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.