<p>I think that it is possible with fours stations yes, not with three.</p>

<p><strong>With the assumption of (roughly) homogeneous velocities we get the following forward problem</strong>. Assuming that <span class="math-container">$v_p$</span> and <span class="math-container">$v_s$</span> are constants between the Earthquake and all your receivers yields the following system:
<span class="math-container">\begin{align}
    d &amp;= \sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2}=\text{3 knowns, 3 unknowns},\quad \text{for} \left\{ \begin{array}{rl} \Delta x &amp;= x_\text{earthquake}-x_\text{receiver} \\ \Delta y &amp; = y_\text{earthquake} - y_\text{receiver} \\ \Delta z &amp; = z_\text{earthquake} - z_\text{receiver}  \end{array} \begin{split} \end{split}\right.,\\
\Delta t &amp; = t_s - t_p = \text{measured},\\
p &amp;= \left(\frac{1}{v_s} - \frac{1}{v_p} \right) = \text{1 unknown}.
\end{align}</span>
Restating the equation of your question then gives the forward problem:
<span class="math-container">$$ \Delta t = p d.$$</span>
For example, for <span class="math-container">$v_p=2000$</span> m/s and <span class="math-container">$v_s=1200$</span> m/s, we find <span class="math-container">$p=1/3000$</span> and can compute e.g. with <span class="math-container">$d = 600$</span> that <span class="math-container">$\Delta t = 600/3000 =0.2$</span> s.</p>

<p><strong>Now, we must set=up the inverse problem</strong>. You generally need at least as many equations as unknowns. In this case we have 4 unknowns (the velocity and the <span class="math-container">$(x,y,z)$</span> earthquake location). That means that we need 4 (unique) equations. I think that the simplest way of obtaining 4 unique equations is to combine 4 measurements:
<span class="math-container">$$ p\begin{pmatrix} d_1 \\ d_2 \\ d_3 \\ d_4 \end{pmatrix} = \begin{pmatrix} \Delta t_1 \\ \Delta t_2 \\ \Delta t_3 \\ \Delta t_4 \end{pmatrix}, $$</span>
Denoting with <span class="math-container">$x_{1,2,3,4}$</span> and <span class="math-container">$y_{1,2,3,4}$</span> the latitude and longitude of the receivers, and assuming measurements at the surface, we obtain:
<span class="math-container">$$ p\begin{pmatrix} \sqrt{(x_\text{earthquake}-x_1)^2 + (y_\text{earthquake}-y_1)^2 + z_\text{earthquake}^2} \\ \sqrt{(x_\text{earthquake}-x_2)^2 + (y_\text{earthquake}-y_2)^2 + z_\text{earthquake}^2} \\ \sqrt{(x_\text{earthquake}-x_3)^2 + (y_\text{earthquake}-y_3)^2 + z_\text{earthquake}^2 } \\ \sqrt{(x_\text{earthquake}-x_4)^2 + (y_\text{earthquake}-y_4)^2 + z_\text{earthquake}^2} \end{pmatrix} = \begin{pmatrix} \Delta t_1 \\ \Delta t_2 \\ \Delta t_3 \\ \Delta t_4 \end{pmatrix}. $$</span></p>

<p><strong>Finally, we must solve the inverse problem</strong>. Unfortunately, the system is non-linear (the unknowns are 'hidden' in the root of the square and can't be separated). You must therefore use some non-linear solver. For example, in MATLAB you can do this</p>

<pre><code>% True earthquake location
x_earthquake=0;
y_earthquake=0;
z_earthquake=1000;

% True station locations
x_1 = 1000; y_1 = 1000;
x_2 =  500; y_2 = -300;
x_3 = -400; y_3 = -100;
x_4 =  -10; y_4 = 800;
X = [x_1;x_2;x_3;x_4];
Y = [y_1;y_2;y_3;y_4];

% True velocity structure
p0 = 1/3000;

% Forward equation (F(1)=x_eartquake, F(2)=y_earthquake, F(3)=z_earthquake, F(4)=p. 
t =@(F) F(4) * sqrt( (F(1)-X).^2 + (F(2)-Y).^2 + F(3).^2);

% True recordings
measured_times = t([x_earthquake,y_earthquake,z_earthquake,p0]);

% Misfit function (=0 at optimum)
G =@(F) t(F) - measured_times;

% Invert
%options = optimoptions('fsolve','FiniteDifferenceType','central'); % using these will give better results in MATLAB
%F_inv = fsolve( G, [1000,1000,1000,0.1],options)
F_inv = fsolve( G, [1000,1000,1000,0.1])
x_earthquake_inv = F_inv(1)
y_earthquake_inv = F_inv(2)
z_earthquake_inv = F_inv(3)
p_inv            = F_inv(4)
</code></pre>

<p>If you copy this into <a href="https://octave-online.net/" rel="nofollow noreferrer">https://octave-online.net/</a>, for example, it will find an earthquake at (0,0,1000) and p=1/3000, with some small errors due to the non-linear solver.</p>

CC BY-SA 4.0
2019-03-20T08:39:01.660Z
2019-03-20T12:03:16.053Z
2019-03-20T12:03:16.053Z
Erik
<p>I'm trying to write a program that would determine the magnitude and epicenter of an earthquake. From <a href="http://www.dartmouth.edu/~ears5/handouts/EQ_distance.html" rel="nofollow noreferrer">http://www.dartmouth.edu/~ears5/handouts/EQ_distance.html</a>, the distance of the epicenter can be calculated with </p>

<p><span class="math-container">$d = (ts-tp)/(1/vs)-(1/vp)$</span></p>

<p>However, I'm not really sure as to what the values Vs and Vp are. From the travel time graph Vs and Vp varies with the distance of the sensor from the earthquake. Given that I have 3 sensors with given latitude and longitude, from which I am able to get the lag time, is there anyway I can estimate the Vs and Vp from these variables?</p>

<p>Otherwise, is there any other way I could get the distance of the sensors from the epicenter?</p>

<earthquakes>
How do I get the velocity of the P-wave or S-wave without using the travel time graph?
nelaaam
# How do I get the velocity of the P-wave or S-wave without using the travel time graph?

(Asked by: username_0 on 2019-03-20)

I'm trying to write a program that would determine the magnitude and epicenter of an earthquake. From <http://www.dartmouth.edu/~ears5/handouts/EQ_distance.html>, the distance of the epicenter can be calculated with

$d = (ts-tp)/(1/vs)-(1/vp)$

However, I'm not really sure as to what the values Vs and Vp are. From the travel time graph Vs and Vp varies with the distance of the sensor from the earthquake. Given that I have 3 sensors with given latitude and longitude, from which I am able to get the lag time, is there anyway I can estimate the Vs and Vp from these variables?

Otherwise, is there any other way I could get the distance of the sensors from the epicenter?

**Tags:** earthquakes


---

## 1 answers

### Answer by username_1 on 2019-03-20. Score: 3

I think that it is possible with fours stations yes, not with three.

**With the assumption of (roughly) homogeneous velocities we get the following forward problem**. Assuming that $v\_p$ and $v\_s$ are constants between the Earthquake and all your receivers yields the following system:
\begin{align}
d &= \sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2}=\text{3 knowns, 3 unknowns},\quad \text{for} \left\{ \begin{array}{rl} \Delta x &= x\_\text{earthquake}-x\_\text{receiver} \\ \Delta y & = y\_\text{earthquake} - y\_\text{receiver} \\ \Delta z & = z\_\text{earthquake} - z\_\text{receiver} \end{array} \begin{split} \end{split}\right.,\\
\Delta t & = t\_s - t\_p = \text{measured},\\
p &= \left(\frac{1}{v\_s} - \frac{1}{v\_p} \right) = \text{1 unknown}.
\end{align}
Restating the equation of your question then gives the forward problem:
$$ \Delta t = p d.$$
For example, for $v\_p=2000$ m/s and $v\_s=1200$ m/s, we find $p=1/3000$ and can compute e.g. with $d = 600$ that $\Delta t = 600/3000 =0.2$ s.

**Now, we must set=up the inverse problem**. You generally need at least as many equations as unknowns. In this case we have 4 unknowns (the velocity and the $(x,y,z)$ earthquake location). That means that we need 4 (unique) equations. I think that the simplest way of obtaining 4 unique equations is to combine 4 measurements:
$$ p\begin{pmatrix} d\_1 \\ d\_2 \\ d\_3 \\ d\_4 \end{pmatrix} = \begin{pmatrix} \Delta t\_1 \\ \Delta t\_2 \\ \Delta t\_3 \\ \Delta t\_4 \end{pmatrix}, $$
Denoting with $x\_{1,2,3,4}$ and $y\_{1,2,3,4}$ the latitude and longitude of the receivers, and assuming measurements at the surface, we obtain:
$$ p\begin{pmatrix} \sqrt{(x\_\text{earthquake}-x\_1)^2 + (y\_\text{earthquake}-y\_1)^2 + z\_\text{earthquake}^2} \\ \sqrt{(x\_\text{earthquake}-x\_2)^2 + (y\_\text{earthquake}-y\_2)^2 + z\_\text{earthquake}^2} \\ \sqrt{(x\_\text{earthquake}-x\_3)^2 + (y\_\text{earthquake}-y\_3)^2 + z\_\text{earthquake}^2 } \\ \sqrt{(x\_\text{earthquake}-x\_4)^2 + (y\_\text{earthquake}-y\_4)^2 + z\_\text{earthquake}^2} \end{pmatrix} = \begin{pmatrix} \Delta t\_1 \\ \Delta t\_2 \\ \Delta t\_3 \\ \Delta t\_4 \end{pmatrix}. $$

**Finally, we must solve the inverse problem**. Unfortunately, the system is non-linear (the unknowns are 'hidden' in the root of the square and can't be separated). You must therefore use some non-linear solver. For example, in MATLAB you can do this

```
% True earthquake location
x_earthquake=0;
y_earthquake=0;
z_earthquake=1000;

% True station locations
x_1 = 1000; y_1 = 1000;
x_2 =  500; y_2 = -300;
x_3 = -400; y_3 = -100;
x_4 =  -10; y_4 = 800;
X = [x_1;x_2;x_3;x_4];
Y = [y_1;y_2;y_3;y_4];

% True velocity structure
p0 = 1/3000;

% Forward equation (F(1)=x_eartquake, F(2)=y_earthquake, F(3)=z_earthquake, F(4)=p. 
t =@(F) F(4) * sqrt( (F(1)-X).^2 + (F(2)-Y).^2 + F(3).^2);

% True recordings
measured_times = t([x_earthquake,y_earthquake,z_earthquake,p0]);

% Misfit function (=0 at optimum)
G =@(F) t(F) - measured_times;

% Invert
%options = optimoptions('fsolve','FiniteDifferenceType','central'); % using these will give better results in MATLAB
%F_inv = fsolve( G, [1000,1000,1000,0.1],options)
F_inv = fsolve( G, [1000,1000,1000,0.1])
x_earthquake_inv = F_inv(1)
y_earthquake_inv = F_inv(2)
z_earthquake_inv = F_inv(3)
p_inv            = F_inv(4)

```

If you copy this into <https://octave-online.net/>, for example, it will find an earthquake at (0,0,1000) and p=1/3000, with some small errors due to the non-linear solver.


---

