前面的话
地震定位是地震学研究里面比较经典，比较基本，比较重要的问题。这篇文章呢，主要关注的是多台站近震定位的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;
}
附程序源码下载
