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

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

《完整版近震定位及其应用本科毕业设计.docx》由会员分享，可在线阅读，更多相关《完整版近震定位及其应用本科毕业设计.docx（15页珍藏版）》请在冰点文库上搜索。

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

目次

摘要………………………………………………………………………………………Ⅰ

Abstract…………………………………………………………………………………Ⅱ

1前言…………………………………………………………………………………1

2近震的理论基础………………………………………………………………………2

2.1地震学的基本名词和概念……………………………………………………2

2.2地震的分类……………………………………………………………………2

2.3近震的主要震相………………………………………………………………3

3正演计算（模型试算）……………………………………………………………5

3.1基本原理………………………………………………………………………5

3.1.1坐标变换………………………………………………………………5

3.1.2正演计算………………………………………………………………6

3.2正演模型的建立………………………………………………………………7

4反演计算………………………………………………………………………………8

4.1基本原理………………………………………………………………………8

4.1.1当速度V未知的初定方法……………………………………………8

4.1.2当速度V已知的初定方法……………………………………………10

4.2模型反演结果…………………………………………………………………11

4.3误差分析………………………………………………………………………11

5实例剖析………………………………………………………………………………13

6结束语…………………………………………………………………………………17

参考文献…………………………………………………………………………………18

附录近震定位程序………………………………………………………………………19

1前言

四川汶川地震和青海玉树地震发生时的地动山摇、房倒屋塌、生死分离的那种凄凉、悲惨的场面，又一次的让我们感受到自然灾害的威力，更是让灾区的人们心惊肉跳，恐慌难安，甚至谈“地震”色变，这就是地震给人们的最直观的印象。

一次强烈地震的发生，常伴随着地面变形和地层错动，其破坏力是相当大的。

主要表现为：

大型建筑物破坏，普通民房破坏，山崩地裂，人畜伤亡。

如今地震发生越来越频繁，2010年也被“尊称”为国际地震年。

任何事情都具有双面性，地震对我们的生命安全造成了很大的威胁的同时也为专家学者们的研究地球提供了丰富的资料。

目前人类对地球内部结构的认识主要是来自地震学研究，因此地震学也就成为地球科学中的一个重要学科。

由于天然地震具有很大的能量，它所产生的地震波可以穿透很大的深度，传播很远的距离，而且当遇到地球深部的波阻抗差异界面时产生反射波，其中，天然地震的中的近震资料为研究地球深部结构提供了一个重要的信息资源和途径。

利用地震记录进行定位始源于欧洲和日本。

最初使用方位角方法，随后是几何作图法和地球投影法。

《国际地震汇编》（ISS）最先采用最小二乘法进行计算修订震中。

1961年，博尔特和威尔莫合作，改进了计算方案，并首先在ISS使用，随后，国际地震中央局与美国海岸和大地测量局先后使用。

在我国，李善邦先生最初使用方位角和最小二乘进行地震的观测和定位。

1953年，我国开始采用大量观测数据修订震中。

目前我国的地震的定位方法兼作图法和计算法。

其中近震的定位方法主要有石川法、和达法、高桥法、外心方位角法、假定发震时刻定位、等时量板法等十几种方法。

计算机的近震定位方法主要分为速度未知的初定为方法和速度已知的初定为方法。

目前定位精度已达到较高的水平。

如果有合理的台网分布和适用的走时表，对于Δ5km；远震为5—10km。

一般浅源地震（Sn、P*PS*S、P*SS*P、PmPSmS及PmSSmP等，如图2.2所示。

图2.2近震震相及其射线路径示意图

（αc、βc分别为康拉德界面及莫霍面的临界角）

1.PgSg震相

由于地球浅部存在强烈的横向及纵向上的速度变化，尤其是有沉积该层地区的垂向梯度的影响，地震波射线在地壳浅部10公里的范围内发生廻折而到达地表，并在震中附近数公里的范围，这种波束与直达波，被称之为廻折或潜波（DivingWave），对于纵、横波来说，分别命名为Pg和Sg。

2.P*S*震相

P*S*分别为地震波在地壳内部的康拉德界面产生的折射纵、横波，它们的速度分别为6.3~7.1kms及3.6~3.9kms。

应当指出，由于康拉德界面并不是全球普遍存在的界面，在一些地区观测不到。

3.PnSn震相

PnSn分别为地震波在地壳底部的莫霍界面产生的折射纵、横波，它们的速度分别为7.9~8.2kms及4.4~4.8kms。

一般来说，Pn震相在莫霍面的临界反射以远地区总是已初至波的形式出现，与其它折射波一样，波至呈现为线性分布。

4.P*PS*S及P*SS*P震相

P*PS*S分别为地震波在地壳内部的康拉德界面产生的反射纵、横波，P*SS*P分别为地震纵横波入射地壳内部的康拉德界面后发生转换而产生的反射横纵波。

3正演计算（模型试算）

测定震源位置及发震时刻最常用的方法是将非线性问题化为线性问题，然后利用最小二乘法原理迭代，全主元消元法等求解。

通过反复修定，以得到震源参数的最佳值。

此方法对近震经纬度的确定比较准确的。

近震的定位要确定四个参数，即震源的空间坐标（，，>4，震中距小于170km，已知直达纵波P到时，而且台站

方位和震中距分布合理的情况下。

设接收到某次地震直达波到时的台站有n个（n≥4）。

其中任一台站的坐标为（xi,yi）相应的直达波到时为Ti。

待求的震中位置为x,y,震源深度z,发震时刻为T（图4.1）；按照均匀地壳模型，可列出下列方程：

图4.1震源计算机定位原理图

（4.1）

i=1,2,…，n.

式中为地震波在地壳中的平均传播速度，因区域不同而异。

式（4.1）为非线性方程组，为便于直接求解，首先应将其线性化，变为线性方程组。

为此将（4.1）式展开得

（4.2）

i=1,2,…,n.

对于第一个台，（i=1）可以写出

（4.3）

将（4.3）式中i分别代以2，3，…，n与（4.3）式相减得

（4.4）

共有n-1个线性方程式。

（4.4）式以消去了深度参数z，只有x、y、T、V四个未知数。

由于n-1一般都大于4，故（4.4）式所表示的线性方程个数多于未知数个数，即构成所谓的超定方程组

………………

（4.5）

为了表达简洁起见，将此超定方程组用矩阵形式表达

（4.6）

其中，对于n个台站的数据,

其中U=V2，W=V2T。

当地震台站数目n小于5时，上述方程组为不定方程，有多解；大于5时，方程组为超定方程，通过变换

（A’为A的转置）

使得超定方程组变为正定方程组，从而计算出震源参数，这就是最小二乘法原理。

只要解出此方程组（4.6），相应震源参数也可求得。

即：

震中位置:

（x，y）

介质的波速:

V＝

发震时刻:

T=WU.

又由：

任取一台站（i）数据代入此式便得震源深度：

4.1.2当波速V已知的初定为方法

因为所研究区域位于青藏高原，其地壳的厚度较大，所以近震的深度基本在上地

壳的范围之内，其速度大约为v=5.8kms.当速度为已知值后，在上述方法中的线性方程组就会减少一个未知数，所以（3.5）变化之后，

（4.7）

其中

W=V2T。

只要解出此方程组（4.7），相应震源参数也可求得。

即：

震中位置:

（x，y）

发震时刻:

T=WV2.

同理，再由

求出震源的深度。

这两种方法的主要流程为：

图4.2震源计算机定位流程图

4.2模型反演结果

为了更好的检验程序，在正演t的结果上均加了6s，下面为各个反演方法所得到的结果：

1.当介质的速度V未知时，两个模型反演出的地震参数分别为：

AnIntroductiontoSeismologyEarthquakes,andEarthStructure.BlackwellPublishing,2003.

[2]李善邦，中国地震.地震出版社,1981.

[3]单娜琳,程志平,刘云祯.工程地震勘探.冶金工业出版社，2006.

[4]姜枚,王有学,钱辉.造山的高原—青藏高原及其邻区的宽频地震探测与地壳上地幔结构.地质出版社，2009.

[5]熊章强,方根显.浅层地震勘探.地震出版社，2002.

[6]大港油田科技丛书编委会.地震勘探资料处理和解释技术.石油工业出版社，1999.

[7]时振梁,张少泉,赵荣国等.地震工作手册.地震出版社，1990.

[8]傅淑芳,刘宝城,李文艺.地震学教程（上、下册）.地震出版社，1980.

[9]彭国伦.Fortran95程序设计.中国电力出版社，2002.

附录程序

programtest

implicitnone

integer:

:

i,j,k,nstn,istn

real*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,,MaxLoncharacter:

:

filename*100

write（*,*）"inputfilename"

read（*,'（a）'）filename

!

readthepicksforallstations

open（1,file=filename）

read（1,*）nstn

x0=0;y0=0

doistn=1,nstn

read（1,*）x（istn）,y（istn）,t（istn）

x（istn）=x（istn）*atan（1.00）45

y（istn）=y（istn）*atan（1.00）45

x0=x0+x（istn）

y0=y0+y（istn）

enddo

la0=x0nstn

ph0=y0nstn

close

（1）

goto10

!

calculatingsytheticdata

la0=99*atan（1.00）45

ph0=27*atan（1.00）45

open（2,file='mod13.dat'）

write（2,'（I2）'）nstn

v=5.0

=1,nstn

lai=x（istn）

phi=y（istn）

calllp2xy（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（istn）+6

enddo

write（2,'（4f20.6）'）la0*45.00atan（1.00）,ph0*45.00atan（1.00）,v,=1,nstn

lai=x（istn）

phi=y（istn）

calllp2xy（lai,phi,la0,ph0,xx,yy）

x（istn）=xx

y（istn）=yy

enddo

!

calculatingcoefficentsforthelinearequations

doistn=2,nstn

coef（istn,1）=x（istn）-x

（1）

coef（istn,2）=y（istn）-y

（1）

coef（istn,3）=（t（istn）**2-t

（1）**2）2.00

coef（istn,4）=-t（istn）+t

（1）

delt（istn）=（x（istn）**2-x

（1）**2+y（istn）**2-y

（1）**2）2.00

enddo

!

leastsquar-rootmethodtogetthesolution

!

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（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）

doi=1,4

doj=1,4

solut（i,j）=0

dok=1,nstn

solut（i,j）=solut（i,j）+coef（k,i）*coef（k,j）

enddo

enddo

enddo

doi=1,4

solut（i,5）=0

dok=1,nstn

solut（i,5）=solut（i,5）+coef（k,i）*delt（k）

enddo

enddo

callgaussian（solut,4,5）

xx=solut（1,5）

yy=solut（2,5）

v0=sqrt（solut（3,5））

t0=solut（4,5）solut（3,5）

=1,nstn

callxy2lp（xx,yy,la0,ph0,lai,phi）

print*,'------laiphiv0t0（1.00）,phi*45.00atan（1.00）,v0,t0,（A,N1,M1）

implicitnone

integer:

:

n1,m1,k1,k2,ik,i,j,l,i1

REAL*8:

:

A（n1,m1）

REAL*8:

:

BMAX,T,EPS

EPS=0.

DOK1=1,N1

BMAX=0.

DOI=K1,N1

IF（BMAX-ABS（A（I,K1））.LT.0）THEN

BMAX=ABS（A（I,K1））

L=I

ENDIF

ENDDO

IF（BMAX.LT.EPS）STOP4444

IF（L.NE.K1）THEN

DOJ=K1,M1

T=A（L,J）

A（L,J）=A（K1,J）

A（K1,J）=T

ENDDO

ENDIF

T=1.A（K1,K1）

K2=K1+1

DOJ=K2,M1

A（K1,J）=A（K1,J）*T

DOI=K2,N1

A（I,J）=A（I,J）-A（I,K1）*A（K1,J）

ENDDO

ENDDO

ENDDO

DOIK=2,N1

I=M1-IK

I1=I+1

DOJ=I1,N1

A（I,M1）=A（I,M1）-A（I,J）*A（J,M1）

ENDDO

ENDDO

RETURN

ENDSUBROUTINEGaussian

subroutinelp2xy（lai,phi,la0,ph0,xx,yy）

implicitnone

real*8:

:

ph0,la0,phi,lai,ee,a,r1,r2,r3,dl,dp,sinxcos,xx,yy

a=6378.160

sinxcos=sin（ph0）*cos（ph0）

R1=sqrt（1-ee*sin（ph0）**2.00）

r2=aR1

r3=a*（1-ee）R1**3

dl=lai-la0

dp=phi-ph0

xx=r2*dl*cos（ph0）

yy=r3*dp+0.50*r2*dl**2.00*sinxcos

return

endsubroutinelp2xy

subroutinexy2lp（xx,yy,la0,ph0,lai,phi）

implicitnone

real*8:

:

ph0,la0,phi,lai,ee,a,r1,r2,r3,sinxcos,xx,yy

a=6378.160

sinxcos=sin（ph0）*cos（ph0）

R1=sqrt（1-ee*sin（ph0）**2.00）

r2=aR1

r3=a*（1-ee）R1**3

lai=xx（r2*cos（ph0））+la0

phi=yyr3-0.50*xx**2*tan（ph0）（r2*r3）+ph0

return

endsubroutinexy2lp

当速度给定时，在求解方程参数和未知数的量稍有变化。

即

doistn=2,nstn

coef（istn,1）=x（istn）-x

（1）

coef（istn,2）=y（istn）-y

（1）

coef（istn,3）=（-t（istn）+t

（1））*v**2

delt（istn）=（x（istn）**2-x

（1）**2+y（istn）**2-y

（1）**2）2.00-（t（istn）**2-t

（1）**2）*v**22.00

enddo

doi=1,3

doj=1,3

solut（i,j）=0

dok=1,nstn

solut（i,j）=solut（i,j）+coef（k,i）*coef（k,j）

enddo

enddo

enddo

doi=1,3

solut（i,4）=0

dok=1,nstn

solut（i,4）=solut（i,4）+coef（k,i）*delt（k）

enddo

enddo