(完整版)近震定位及其应用本科毕业设计.doc

1、目次摘要 ,Abstract,1前言 ,12近震的理论基础 ,22.1地震学的基本名词和概念 ,22.2地震的分类 ,22.3近震的主要震相 ,33正演计算（模型试算） ,53.1基本原理 ,53.1.1坐标变换 ,53.1.2正演计算 ,63.2正演模型的建立 ,74反演计算 ,84.1基本原理 ,84.1.1当速度 V 未知的初定方法 ,84.1.2当速度 V 已知的初定方法 ,104.2模型反演结果 ,114.3误差分析 ,115实例剖析 ,136结束语 ,17参考文献 ,18附录 近震定位程序 ,191 前 言四川汶川地震和青海玉树地震发生时的地动山摇、房倒屋塌、生死分离的那种凄凉、悲

2、惨的场面，又一次的让我们感受到自然灾害的威力，更是让灾区的人们心惊肉跳，恐慌难安，甚至谈“地震”色变，这就是地震给人们的最直观的印象。一次强烈地震的发生，常伴随着地面变形和地层错动，其破坏力是相当大的。主要表现为：大型建筑物破坏，普通民房破坏，山崩地裂，人畜伤亡。如今地震发生越来越频繁， 2010 年也被“尊称”为国际地震年。任何事情都具有双面性，地震对我们的生命安全造成了很大的威胁的同时也为专家学者们的研究地球提供了丰富的资料。目前人类对地球内部结构的认识主要是来自地震学研究，因此地震学也就成为地球科学中的一个重要学科。由于天然地震具有很大的能量，它所产生的地震波可以穿透很大的深度，传播很远

3、的距离，而且当遇到地球深部的波阻抗差异界面时产生反射波，其中，天然地震的中的近震资料为研究地球深部结构提供了一个重要的信息资源和途径。利用地震记录进行定位始源于欧洲和日本。最初使用方位角方法，随后是几何作图法和地球投影法。国际地震汇编（ISS）最先采用最小二乘法进行计算修订震中。1961年，博尔特和威尔莫合作，改进了计算方案，并首先在ISS使用，随后，国际地震中央局与美国海岸和大地测量局先后使用。在我国，李善邦先生最初使用方位角和最小二乘进行地震的观测和定位。1953年，我国开始采用大量观测数据修订震中。目前我国的地震的定位方法兼作图法和计算法。其中近震的定位方法主要有石川法、和达法、高桥法、

4、外心方位角法、假定发震时刻定位、等时量板法等十几种方法。计算机的近震定位方法主要分为速度未知的初定为方法和速度已知的初定为方法。目前定位精度已达到较高的水平。如果有合理的台网分布和适用的走时表，对于5km；远震为 5 10km。一般浅源地震（ Sn、P*PS*S、P*SS*P、 PmPSmS及 PmSSmP等，如图 2.2 所示。图 2.2近震震相及其射线路径示意图（c、c 分别为康拉德界面及莫霍面的临界角）1. PgSg 震相由于地球浅部存在强烈的横向及纵向上的速度变化，尤其是有沉积该层地区的垂向梯度的影响，地震波射线在地壳浅部10 公里的范围内发生廻折而到达地表，并在震中附近数公里的范围，

5、这种波束与直达波，被称之为廻折或潜波（Diving Wave ），对于纵、横波来说，分别命名为Pg 和 Sg。2. P*S* 震相P*S* 分别为地震波在地壳内部的康拉德界面产生的折射纵、横波，它们的速度分别为 6.37.1kms 及 3.63.9kms 。应当指出， 由于康拉德界面并不是全球普遍存在的界面，在一些地区观测不到。3. PnSn 震相PnSn 分别为地震波在地壳底部的莫霍界面产生的折射纵、横波，它们的速度分别为 7.98.2kms 及 4.44.8kms 。一般来说， Pn 震相在莫霍面的临界反射以远地区总是已初至波的形式出现，与其它折射波一样，波至呈现为线性分布。4. P*PS

6、*S 及 P*SS*P 震相P*PS*S 分别为地震波在地壳内部的康拉德界面产生的反射纵、横波，P*SS*P分别为地震纵横波入射地壳内部的康拉德界面后发生转换而产生的反射横纵波。3 正演计算（模型试算）测定震源位置及发震时刻最常用的方法是将非线性问题化为线性问题，然后利用最小二乘法原理迭代，全主元消元法等求解。通过反复修定，以得到震源参数的最佳值。此方法对近震经纬度的确定比较准确的。近震的定位要确定四个参数，即震源的空间坐标（， >4，震中距小于170km，已知直达纵波 P 到时，而且台站方位和震中距分布合理的情况下。设接收到某次地震直达波到时的台站有n 个（n4）。其中任一台站的坐标为

7、（ xi ,y i ）相应的直达波到时为Ti 。待求的震中位置为x,y, 震源深度 z, 发震时刻为 T（图 4.1 ）；按照均匀地壳模型，可列出下列方程：图 4.1震源计算机定位原理图(x xi ) 2( y yi ) 2z2V 2 (T Ti ) 2（4.1 ）i=1,2, ， n.式中为地震波在地壳中的平均传播速度，因区域不同而异。式（4.1 ）为非线性方程组，为便于直接求解，首先应将其线性化，变为线性方程组。为此将（ 4.1 ）式展开得x 22xi x x 2y 22yi y yi2z2V 2 (Ti22TTiT 2 )(4.2)i=1,2,n.对于第一个台，（ i=1 ）可以写出x2

8、2x1 x x2y 22 y1 y y12z2V2(T22TTT 2 )(4.3)11将(4.3) 式中 i 分别代以 2，3，, ， n 与 (4.3) 式相减得(xix1 )x ( yiy1 ) y (ti2t12 )V 2 / 2 (tit1 )TV 2( xi2x12yi2y12 ) / 2（4.4 ）共有 n-1 个线性方程式。 (4.4) 式以消去了深度参数z，只有 x、 y、T、V 四个未知数。由于 n-1 一般都大于 4，故 (4.4) 式所表示的线性方程个数多于未知数个数，即构成所谓的超定方程组( x2x1 )x ( y2y1 ) y (t22t12 )V 2( x2x2y

9、2y 2 ) / 22121( x3x1 ) x ( y3y1 ) y (t32t12 )V 2( x2x2y 2y 2 ) / 23131,( xnx1 ) x ( y ny1 ) y (t n2t12 )V 2( x2x 2y 2y 2 ) / 2n1n1为了表达简洁起见，将此超定方程组用矩阵形式表达(4.6)其中，对于 n 个台站的数据 ,/ 2(t2t1 )TV 2/ 2(t 3t1 )TV 2/ 2 (t nt1 )TV 2（4.5 ）x2x1y2y1(t22t12 ) / 2t1t 2x3x1y3y1(t32t12 ) / 2t1t 3Axnx1yny1(tn2t12 ) / 2t

10、1t n其中2U=V，2W=VT。当地震台站数目n小于 5时，上述方程组为不定方程，有多解；大于5时，方程组为超定方程，通过变换（A为 A的转置）使得超定方程组变为正定方程组，从而计算出震源参数，这就是最小二乘法原理。只要解出此方程组 (4.6) ，相应震源参数也可求得。即：震中位置 :（ x， y）介质的波速 :V发震时刻 :T=WU.又由：( x xi )2( y yi )2z2V 2 (T Ti ) 2任取一台站 (i) 数据代入此式便得震源深度：zV 2 (TTi ) 2( xxi )2( yyi )24.1.2当波速 V 已知的初定为方法因为所研究区域位于青藏高原，其地壳的厚度较大，

11、所以近震的深度基本在上地壳的范围之内， 其速度大约为 v=5.8kms. 当速度为已知值后， 在上述方法中的线性方程组就会减少一个未知数，所以（3.5 ）变化之后，( xnx1 ) x ( yny1 ) y (tnt1 )TV 2( xn2x12yn2y12 ) / 2 (tn2（4.7 ）t12 )V 2 / 2其中x2x1y2y1t1t 2Ax3x1y3y1t1t3xnx1yny1t1tnx22x12y22y12t22t121x32x12y32y12t32t12B2xn2x12yn2y12tn2t122W=VT。只要解出此方程组 (4.7) ，相应震源参数也可求得。即：震中位置 :（ x，

12、 y）发震时刻 :T=W V2.同理，再由 zV 2 (TTi )2(xxi ) 2( yyi ) 2 求出震源的深度。这两种方法的主要流程为：开始台站经纬度转成直角坐标值计算方程组系数矩阵A和常数矩阵B求解方程组AX=B依 方 程 组 的 解 进 一 步 求 得 震 源 参 数 (x,y,z,T)作出震源与台站相对位置图形直角坐标值转回经纬度表示输出结束图 4.2震源计算机定位流程图4.2模型反演结果为了更好的检验程序， 在正演 t 的结果上均加了6s，下面为各个反演方法所得到的结果：1. 当介质的速度 V 未知时，两个模型反演出的地震参数分别为：An IntroductiontoSeism

13、ologyEarthquakes,andEarthStructure.BlackwellPublishing, 2003.2 李善邦，中国地震 . 地震出版社 ,1981.3 单娜琳 , 程志平 , 刘云祯 . 工程地震勘探 . 冶金工业出版社， 2006.4 姜枚 , 王有学 , 钱辉 . 造山的高原青藏高原及其邻区的宽频地震探测与地壳上地幔结构 . 地质出版社， 2009.5 熊章强 , 方根显 . 浅层地震勘探 . 地震出版社， 2002.6 大港油田科技丛书编委会 . 地震勘探资料处理和解释技术 . 石油工业出版社， 1999.7 时振梁 , 张少泉 , 赵荣国等 . 地震工作手册 .

14、 地震出版社， 1990.8 傅淑芳 , 刘宝城 , 李文艺 . 地震学教程（上、下册） . 地震出版社， 1980.9 彭国伦 .Fortran 95 程序设计 . 中国电力出版社， 2002.附录程序program testimplicit noneinteger : i, j, k, nstn, istnreal*8 : x0, y0, v0, t0, xx, yy, coef(1000,5)real*8 : delt(1000), solut(4,5), x(1000), y(1000), t(1000)real*8 : ph0, la0, phi, lai ,v, MaxLoncha

15、racter:filename*100write(*,*) "input filename"read(*,'(a)') filename!read the picks for all stationsopen(1,file=filename)read(1,*) nstnx0=0;y0=0do istn=1,nstnread(1,*) x(istn),y(istn), t(istn)x(istn)=x(istn)*atan(1.00)45y(istn)=y(istn)*atan(1.00)45x0=x0+x(istn)y0=y0+y(istn)end dola

16、0=x0nstnph0=y0nstnclose(1)goto 10!calculating sythetic datala0=99*atan(1.00)45ph0=27*atan(1.00)45open(2,file='mod13.dat')write(2,'(I2)') nstnv=5.0=1,nstnlai=x(istn)phi=y(istn)call lp2xy(lai,phi,la0,ph0,xx,yy)t(istn)=sqrt(xx*2.00+yy*2.00+)*45.00atan(1.00),y(istn)*45.00atan(1.00), t(is

17、tn)+6end dowrite(2,'(4f20.6)')la0*45.00atan(1.00),ph0*45.00atan(1.00),v,=1,nstnlai=x(istn)phi=y(istn)call lp2xy(lai,phi,la0,ph0,xx,yy)x(istn)=xxy(istn)=yyend do!calculating coefficents for the linear equationsdo istn=2,nstncoef(istn,1)=x(istn)-x(1)coef(istn,2)=y(istn)-y(1)coef(istn,3)=(t(ist

18、n)*2-t(1)*2)2.00coef(istn,4)=-t(istn)+t(1)delt(istn) =(x(istn)*2-x(1)*2+y(istn)*2-y(1)*2)2.00end do!least squar-root method to get the solution! dt(1)de(1,1)de(1,2) de(1,3) de(1,4) de(1,5)! dt(2)de(2,1)de(2,2) de(2,3) de(2,4) de(2,5)! .! .! .! dt(n)de(n,1)de(n,2) de(n,3) de(n,4) de(n,5)! de(1,1) de(

19、2,1) . de(n,1)! de(1,2) de(2,2) . de(n,2)! de(1,3) de(2,3) . de(n,3)! de(1,4) de(2,4) . de(n,4)! de(1,5) de(2,5) . de(n,5) do i=1,4do j=1,4 solut(i,j)=0 do k=1,nstnsolut(i,j)=solut(i,j)+coef(k,i)*coef(k,j)end doend doend dodo i=1,4solut(i,5)=0do k=1,nstnsolut(i,5)=solut(i,5)+coef(k,i)*delt(k)end doe

20、nd docall gaussian(solut,4,5)xx=solut(1,5)yy=solut(2,5)v0=sqrt(solut(3,5)t0=solut(4,5)solut(3,5)=1,nstncall xy2lp(xx,yy,la0,ph0,lai,phi)print*,'-laiphiv0t0(1.00),phi*45.00atan(1.00),v0,t0, (A,N1,M1)implicit noneinteger :n1,m1,k1,k2,ik,i,j ,l, i1REAL*8 : A(n1,m1)REAL*8 : BMAX,T,EPSEPS=0.DO K1=1,N

21、1BMAX=0.DO I=K1,N1IF(BMAX-ABS(A(I,K1).LT.0) THENBMAX=ABS(A(I,K1)L=IEND IFEND DOIF(BMAX.LT.EPS) STOP 4444IF(L.NE.K1) THENDO J=K1,M1T=A(L,J)A(L,J)=A(K1,J)A(K1,J)=TEND DOEND IFT=1.A(K1,K1)K2=K1+1DO J=K2,M1A(K1,J)=A(K1,J)*TDO I=K2,N1A(I,J)=A(I,J)-A(I,K1)*A(K1,J)END DOEND DOEND DODO IK=2,N1I=M1-IKI1=I+1D

22、O J=I1,N1A(I,M1)=A(I,M1)-A(I,J)*A(J,M1)END DOEND DORETURNEND SUBROUTINE Gaussiansubroutine lp2xy(lai,phi,la0,ph0,xx,yy)implicit nonereal*8 : ph0, la0, phi, lai, ee, a, r1,r2, r3, dl, dp, sinxcos, xx,yya=6378.160sinxcos=sin(ph0)*cos(ph0)R1=sqrt(1-ee*sin(ph0)*2.00)r2=aR1r3=a*(1-ee)R1*3dl=lai-la0dp=phi-ph0xx=r2*dl