前面的话

地震定位是地震学研究里面比较经典，比较基本，比较重要的问题。这篇文章呢，主要关注的是多台站近震定位的Inglada算法以及和达法，这两个方法算是比较基础简单的了，感觉假设条件挺多的，方法略粗糙简单。

Inglada算法

主要参照万永革教授编著的《地震学导论》一书。其主要思想是利用多个台站的P波、S波的到时信息，同时引入虚波速度的概念，震源到台站的距离就可以表示为虚波速度与P波和S波到时差的乘积。而同时，震源到台站的距离又可以用震源坐标和台站的坐标计算出来，只不过是一个完全平方求和的公式，这样，对于每个台站应用，可以列出一个线性方程组来，解线性方程组即可得到震源的位置，有了震源位置，发震时刻也容易求出来。需要注意的问题是组成的线性方程组，当台站的海拔高度相差不多，使得方程组的系数矩阵具有奇异性，导致求解的震源深度发散，所以可以先求出震源的平面坐标，再求得震源深度，将各个台站求得的震源深度进行平均。我这么说，只有文字，好像也看不懂我说了啥，唉，这个不能插入数学公式，我也木办法，所以，还是参照万永革教授编著的《地震学导论》323页的内容吧。下面是利用Matlab程序实现Inglada算法的程序。

function [FirLocal]=inglada(Vp,Vs,Pg,Sg,Stx,Sty,Stz)
VSpeed=(Vp*Vs)/(Vp-Vs);%虚波的速度
R=VSpeed.*(Sg-Pg);%震源到台站的距离
num=length(Stx);%计算台站的个数

%r的定义
r=zeros(num);
for i=1:num
    r(i)=sqrt(power(Stx(i),2)+power(Sty(i),2)+power(Stz(i),2));
end

%对系数矩阵进行初始化
A=zeros(num-1,3);
%对常数项矩（向量）进行初始化
B=zeros(num-1,1);
%对每个台站计算出来的震源深度矩阵进行初始化
SeisDepth=zeros(num,1);
%对每个台站计算出来的发震时刻矩阵进行初始化
SeisTime=zeros(num,1);

%系数矩阵的定义
for i=1:num-1
    A(i,1)=Stx(1)-Stx(i+1);
end

for i=1:num-1
    A(i,2)=Sty(1)-Sty(i+1);
end

for i=1:num-1
    A(i,3)=Stz(1)-Stz(i+1);
end

%常数项矩阵（向量）的定义
for i=1:num-1
    B(i,1)=1/2*(power(r(1),2)-power(r(i+1),2)+power(R(i+1),2)-power(R(1),2));
end
%求解方程组，pinv用于求解伪逆
FirLocal=pinv(A)*B;

%求解震源深度
for i=1:num
    SeisDepth(i,1)=sqrt(power(R(i),2)-power(FirLocal(1)-Stx(i),2)-power(FirLocal(2)-Sty(i),2))+Stz(i);
end
SeisDepth_M=mean(SeisDepth);
FirLocal(3)=SeisDepth_M;

%当震源深度不符合常理时，强制震源深度为10km，什么鬼，我也不知道为啥是12
if (FirLocal(3)<0 || FirLocal(3)>500)
    FirLocal(3)=12;
end

%计算发震时刻
for i=1:num
    SeisTime(i,1)=Pg(i)-R(i)/Vp;
end

SeisTime_M=mean(SeisTime);
FirLocal(4)=SeisTime_M;

return
    
%输入台站的坐标和P波和S波的到时
Stx=[50 0 50 100 100 100 50 0 0];
Sty=[50 100 100 100 50 0 0 0 50];
Stz=zeros(9);
Pg=[6.4 18.5 11.9 11.9 6.4 11.9 11.9 18.5 15.5];
Sg=[10.7 30.8 19.8 19.8 10.7 19.8 19.8 30.8 25.9];
[FirLocal]=inglada(5,3,Pg,Sg,Stx,Sty,Stz)

和达直线法

主要参照万永革教授编著的《地震学导论》一书。和达直线法的思想是利用在地壳为均匀介质模型下，P波的到时与P波和S波到时差呈直线关系，这个直线的结局就是发震时刻，求这个发震时刻的方法就是和达直线法。我们有P波和S波的到时信息，利用最小二乘法进行直线的拟合得到截距就是发震时刻了。所以程序的核心是最小二乘法的实现。下面是C程序实现和达直线法发震时刻的求解。

/*******************************************************************************/
/***************************LS_Locate.c     2017/10/23**************************/
/***********************First edition written by SeisBird***********************/
/*This code aims to calculate original time of earthquake using Wadachi method */
#include <stdio.h>
#include <math.h>
double a,b; /*the coefficent is b (S-P travel times difference), the intercept(original time of earthquake is a)*/
double Sa,Sb; /*the standard deviation of a(Sa) and b(Sb)*/
double r; /*the correlation coefficent*/
#define length 9 /*the number of the  stations*/
/*******************************************************/
/***********The Least Square Method Function************/
void LeastSquare(double X[], double Y[])
{
    int i;
    double Sy;
    double X_Sum=0,Y_Sum=0,XX_Sum=0,XY_Sum=0,YY_Sum=0;
    double X_Mean,Y_Mean,XX_Mean,XY_Mean,YY_Mean;
    double XX[length],XY[length],YY[length];
    for(i=0;i<length;i++)
    {
    XX[i]=pow(X[i],2);
    YY[i]=pow(Y[i],2);
    XY[i]=X[i]*Y[i];
    }

    for(i=0;i<length;i++)
    {
    X_Sum=X_Sum+X[i];
    Y_Sum=Y_Sum+Y[i];
    XX_Sum=XX_Sum+XX[i];
    XY_Sum=XY_Sum+XY[i];
    YY_Sum=YY_Sum+YY[i];
    }
    X_Mean=X_Sum/length;
    Y_Mean=Y_Sum/length;
    XX_Mean=XX_Sum/length;
    XY_Mean=XY_Sum/length;
    YY_Mean=YY_Sum/length;
    b=(XY_Mean-X_Mean*Y_Mean)/(XX_Mean-pow(X_Mean,2));
    a=Y_Mean-b*X_Mean;
    double SS=0;
    for(i=0;i<length;i++)
    {
    SS=SS+pow(Y[i]-b*X[i]-a,2);
    }
    Sy=sqrt(SS/(length-2));
    Sa=sqrt(pow(X_Mean,2)/(length*(XX_Mean-pow(X_Mean,2))))*Sy;
    Sb=sqrt(1/(length*(XX_Mean-pow(X_Mean,2))))*Sy;
    r=(XY_Mean-X_Mean*Y_Mean)/sqrt((XX_Mean-pow(X_Mean,2))*(YY_Mean-pow(Y_Mean,2)));    
}

/*input P-waves travel time and S-waves travel time*/
/*output a,b and r*/
int main(void)
{
    int i;
    double TP[length]={13.0249,15.1966,17.0544,19.2718,21.0893,23.1923,25.0391,27.4119,29.0323};
    double TS[length]={20.0511,23.2570,26.4523,30.1574,33.2776,36.4374,39.7976,43.2971,46.3749};
    double DT[length];
    for (i=0;i<length;i++)
    {
    DT[i]=TS[i]-TP[i];
    }
    LeastSquare(DT,TP);
    printf("a is %f+-%f, b is %f+-%f\n",a,Sa,b,Sb);
    printf("r is %f\n",r);
    return 0;
}

附程序源码下载

1. Inglada算法地震定位
2. 和达直线法求发震时刻LS_Locate.c