前面的話
地震定位是地震學(xué)研究里面比較經(jīng)典陵霉，比較基本说搅，比較重要的問題。這篇文章呢律歼，主要關(guān)注的是多臺站近震定位的Inglada算法以及和達法民镜，這兩個方法算是比較基礎(chǔ)簡單的了，感覺假設(shè)條件挺多的险毁，方法略粗糙簡單制圈。
Inglada算法
主要參照萬永革教授編著的《地震學(xué)導(dǎo)論》一書。其主要思想是利用多個臺站的P波畔况、S波的到時信息鲸鹦，同時引入虛波速度的概念狮惜，震源到臺站的距離就可以表示為虛波速度與P波和S波到時差的乘積塑悼。而同時，震源到臺站的距離又可以用震源坐標和臺站的坐標計算出來辩稽，只不過是一個完全平方求和的公式吵瞻，這樣葛菇，對于每個臺站應(yīng)用，可以列出一個線性方程組來橡羞，解線性方程組即可得到震源的位置眯停，有了震源位置，發(fā)震時刻也容易求出來卿泽。需要注意的問題是組成的線性方程組莺债，當臺站的海拔高度相差不多，使得方程組的系數(shù)矩陣具有奇異性签夭，導(dǎo)致求解的震源深度發(fā)散齐邦，所以可以先求出震源的平面坐標，再求得震源深度覆致，將各個臺站求得的震源深度進行平均侄旬。我這么說肺蔚，只有文字煌妈，好像也看不懂我說了啥，唉宣羊，這個不能插入數(shù)學(xué)公式璧诵，我也木辦法，所以仇冯，還是參照萬永革教授編著的《地震學(xué)導(dǎo)論》323頁的內(nèi)容吧之宿。下面是利用Matlab程序?qū)崿F(xiàn)Inglada算法的程序。
function [FirLocal]=inglada(Vp,Vs,Pg,Sg,Stx,Sty,Stz)
VSpeed=(Vp*Vs)/(Vp-Vs);%虛波的速度
R=VSpeed.*(Sg-Pg);%震源到臺站的距離
num=length(Stx);%計算臺站的個數(shù)
%r的定義
r=zeros(num);
for i=1:num
r(i)=sqrt(power(Stx(i),2)+power(Sty(i),2)+power(Stz(i),2));
end
%對系數(shù)矩陣進行初始化
A=zeros(num-1,3);
%對常數(shù)項矩（向量）進行初始化
B=zeros(num-1,1);
%對每個臺站計算出來的震源深度矩陣進行初始化
SeisDepth=zeros(num,1);
%對每個臺站計算出來的發(fā)震時刻矩陣進行初始化
SeisTime=zeros(num,1);
%系數(shù)矩陣的定義
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
%常數(shù)項矩陣（向量）的定義
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
%計算發(fā)震時刻
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)
和達直線法
主要參照萬永革教授編著的《地震學(xué)導(dǎo)論》一書等缀。和達直線法的思想是利用在地殼為均勻介質(zhì)模型下枷莉，P波的到時與P波和S波到時差呈直線關(guān)系，這個直線的結(jié)局就是發(fā)震時刻尺迂，求這個發(fā)震時刻的方法就是和達直線法笤妙。我們有P波和S波的到時信息，利用最小二乘法進行直線的擬合得到截距就是發(fā)震時刻了噪裕。所以程序的核心是最小二乘法的實現(xiàn)蹲盘。下面是C程序?qū)崿F(xiàn)和達直線法發(fā)震時刻的求解。
/*******************************************************************************/
/***************************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;
}
附程序源碼下載
