The general location problem

Setting up the problem

The model parameters

The earthquake location can be defined by its hypocenter (latitude, longitude, and depth). The earthquake's origin time is also unknown. We therefore need to solve for four different model parameters.
Define the vector x as the estimate of the source parameters, where T₀ is the event origin time (s), h is the depth (km), Q is latitude (°N), and F is longitude (°E):

x = [T₀, h, Q, F]^T

Now define a small perturbation to these source parameters:

dx = [dT₀, dh, dQ, dF]^T

The data

T_i^obs = arrival time observed at the i-th station.

How to relate the data to the model parameters

T_i^calc = predicted arrival time = T₀ + travel time calculated from source to the i-th station.
For this lab, the theoretical travel times are obtained by interpolations of travel time tables that represent the 1D Earth model ak135.

T_i^obs = T_i^calc + dT_i/dT⁰ dT⁰ + dT_i/dh dh + dT_i/dQ dQ + dT_i/dF dF + higher order terms

We can thus relate the vector with arrival time residuals to the vector with model perturbations:

dT_i = T_i^obs - T_i^calc = A_ijdx_j + higher order terms

where A_ij = dT_i/dx_j,   i = 1, ...N; j = 1,...4
A_i1 = 1
A_i2 = -cos(i_i)/v
A_i3 = -dT/dDcos(z_i) = -p_icos(z_i)
A_i4 = --dT/dDsin(z_i)cos(Q_i) = -p_isin(z_i)cos(Q_i)

where i_i is the take-off angle at the source for the i-th station (angle between vertical and ray), and z_i is the azimuth (angle measured clockwise from North to the source-to-receiver ray direction).

Solving the problem

If the initial event location estimate is good enough, we can discard the higher order terms and linearize the problem. We solve this set of linear equations by performing a 'least-squares inversion'. We'll be finding the vector dx that minimizes the difference between observed and calculated arrival time residuals. Because the location problem is nonlinear, we have to do this iteratively. For every new location we find, we have to recalculate the partial derivatives (the take-off angle and azimuth will change slightly when the location is different).

dT = Adx

A^TdT = A^T A dx

A^TA is a symmetric matrix and it thus has an inverse. Its inverse can be evaluated using scilab command 'inv()'. Premultiply both sides of the previous equation by (A^TA)^-1 to obtain dx:

dx = (A^TA)^-1 A^T dT

This is the least squares solution to the linear system of equations.
Now, add this perturbation to the previous estimate of the model parameters to obtain a new solution:

x^new = x + dx

You can iterate this procedure until you find an acceptable solution. One can use several criteria for this. We use the criterium that the variance of the arrival time residuals has to be minimized.
e = S_i (T_i^obs - T_i^calc)² / (ndata - nmodel)

ndata = number of data points (observations)
nmodel = number of model parameters (3 or 4, depending on whether you keep the depth fixed)

In fact, we minimize a quantity called Root Mean Square residual, the square root of the variance. It gives a direct indication of how well we match the observed arrival times.