地震定位研究及应用综述概要.docx

地震定位研究及应用综述概要.docx

《地震定位研究及应用综述概要.docx》由会员分享，可在线阅读，更多相关《地震定位研究及应用综述概要.docx（11页珍藏版）》请在冰豆网上搜索。

地震定位研究及应用综述概要

第25卷第1期

2005年2月地　震　工　程　与　工　程　振　动EARTHQUAKEENGINEERINGANDENGINEERINGVIBRATIONVol.25,No.1Feb.2005收稿日期:

2004-10-12

　　基金项目:

国家科技基础条件平台项目（2003DIB2J099;国家自然科学基金项目（50378086;中国地震局“十五”重点科研项目　　作者简介:

杨文东（1969-,男,高级工程师,主要从事地震监测预报.

文章编号:

100021301（20050120014207

地震定位研究及应用综述

杨文东1,金　星2,3,4,李山有2,马强2

（1.黑龙江省地震局,黑龙江哈尔滨150090;2.中国地震局工程力学研究所,黑龙江哈尔滨150080;

3.南京工业大学岩土工程研究所,江苏南京210009;4.福建省地震局,福建福州350003

摘要:

地震定位是地震学中最经典、最基本的问题之一,提高定位精度也一直是地震学应用研究的重

要课题之一。

本文就目前各种地震定位方法进行了大体分类并概述了基本原理及其应用,介绍了目

前在国内各实时地震观测系统中应用的台网定位方法,同时将部分台网方法应用于黑龙江省绥棱爆

破的定位。

关键词:

地震定位;定位方法;地震台网;计算方法;程序

中图分类号:

P315.916　　　文献标识码:

A

Studyofseisi1,Jin3,Li2,MaQiang2

（1.EarhtquakePr150090,China;2.InstituteofEngineeringMechanics,ChinaEarthquake

China;3.NanjingUniversityofTechnology,Nanjing210009,China;

4.EarthquakeAdministrationofFujianProvince,Fuzhou350003,China

Abstract:

Earthquakelocatingisoneoftheclassicandbasicproblemsinseismology,andtoimprovelocatingpre2cisionisstillanimportantprobleminseismologicalapplicationstudy.Quiddityofearthquakelocatingiscomputingaminimumofobjectivefunction.Therearesomedifferencesbetweenobsevrationtravel2timeandtheorytravel2timeofbothstationsandimaginalhypocenter.Manykindsofmethodsofearthquakelocationaredevelopedby3D2veloc2itymodel.Thispapersortsallearthquakelocatingmethodscurrentlyusedandsummarizesbasictheoryandapplica2tionofthem.Itintroducesnetworklocatingmethodsusedbyeachearthquakereal2timeobservationsysteminChi2na,atthesametimesomeofthesemethodswereusedinblastinglocatinginSuiling,HeilongjiangProvince.Keywords:

earthquakelocation;seismicnetwork;computingmethods;program.

引言

地震定位在地震预报、工程地震以及地球结构研究、地壳应力场分析等诸多方面都是一个基础。

同时,基于快速准确的地震定位的地震速报,对于地震应急、震后减灾和救灾及震后地震趋势预测等至关重要,另外近年发展起来的地震预警技术也需要发展高速准确的地震定位新方法。

地震定位一般指确定地震的震源位置（经度、纬度和深度及发震时刻,并适当给出对定位结果的评价,它的确定受到台网布局、可用定位震相、地震波到时的读数精度、地壳速度模型等诸多因素的影响。

研究地震定位方法和提高地震定位精度,一直是地震科学中的一个重要课题,地震学家在不断改进或提出新的定位方法,期望得到更高的地震定位精　　

度,朱元清等[1]曾对地震定位精度的可能误差进行了较详细的分析。

早期地震定位方法主要是几何作图

法,它的历史可以追溯到地震仪问世的年代,并已形成多种经典的方法。

计算定位起源于1912年,由德国物

理学家Geiger[2]提出,其实质是将非线性方程组线性化,并通过最小二乘原理求解。

随着计算方法和计算机

技术的发展,许多定位方法得到发展,如:

Powell法、联合定位法、双差定位法等,随着多学科的相互交流与渗透,其它学科的最优化方法也得以借鉴,如生物学中的遗传算法和物理学中的模拟退火法等,尤其是近年来基于现代数字地震观测技术和科学计算以及计算机技术的智能化数值自动定位方法得到了迅速发展,并业已成为当前地震定位的主流方法。

随着全球及区域速度结构三维层析成像的研究,在此基础上应用三维结构的地震定位已经被人们所关注,如美国哈佛大学一直将提高定位精度列为主要的研究方向之一,并侧重于

探索使用三维速度模型为主[3]。

本文将对地震定位的一些重要方法进行分析,并说明未来地震定位的发展

方向。

1　传统的几何地震定位方法

几何定位法是以走时方程为理论依据,根据已知条件,借助平面或球面作图确定震源位置的方法;如:

利用和达法有4个或以上台站的直达波Pg、Sg到时,就可定出震源位置。

假设有2个台站S1、S2,根据均匀介质的直达波方程可得到它们组成的方程组,对其进行整理和变换,最后得到震源轨迹方程

（x-x02+（y-y02+h2=20（1

　　其中（x0,y0,0为球心,R0为半径的球面方程。

这样,4,其解就是震源点。

高桥法、交切法等,都适用均匀介质模型的轨迹法,,,个台不能在一条直线上,4个台分布不能共圆,。

作图法适,位中。

2　线性计算地震定位方法

2.1　绝对计算方法

1912年德国物理学家Geiger提出了一种对理想地球在时间域内的震源计算定位方法,其实质就是将非线性方程组线性化,并通过最小二乘原理求解。

设n个台站的观测到时为t1,t2…,tn,求震源（x0,y0,z0及发震时刻t0,使得目标函数:

<（t0,x0,y0,z0=∑n

i=1r

2i（2

最小。

ri为到时残差:

ri=ti-t0-Ti（x0,y0,z0,Ti为震源到第i个台站的计算走时。

使目标函数取极小值也即:

g（θ=θ<（θ=0

（3由（3式,在真解θ附近任意试探解θ3及其校正矢量δ

θ满足:

g（θ3+Δθg（θ3Tδθ=0（4

由<的定义可得公式（4的具体表达式:

N

i=1i9θji92r1k+2ri9θj9θ3δθj=-ni=Irii9θ3

（5　　若θ3偏离真解θ不大,则ri（θ3和

9θj9较小,忽略二阶导数项式（5简化为:

ni=Ii9θji9=-ni=1ri9θ3（6以矩阵形式表示上式为:

ATAδθ=ATr（7

511期杨文东等:

地震定位研究及应用综述

　　由方程（7求得δθ后,以θ=θ3+δ

θ为新尝试点,再求解相应方程,如此反复迭代,直至<足够小或满足一定的循环结束条件,即得估计解θ。

Geiger法的应用起源于20世纪70年代,用其理论Lee等[4]连续给出了HYPO71~86系列程序现已形

成HYPOELLIPSE定位程序[5],我国赵仲和[6]参与了这项工作;姜燕1986年编写了BLOC86,经过多次修改补充而形成NC91地震定位和震级计算程序。

美国国家地震情报中心（NICE的EDR报告、中国地震局地球物理研究所出版的《中国地震台临时报告》和《中国地震年报》等报告还有各省的地震编目工作仍然采用这种方法。

Geiger法对初始值的依赖性较大,有时会使得定位发散或得不到定位结果。

在实际地震观测中,如果台站较好的包围震中,则可以得到较精确的定位结果;如果台站分布不合理,地震偏在台网的一侧等,则使得方程是一个病态方程,由最小二乘法求解会不稳定,得不到正确的震源位置或解发散。

为了解决这些问题,人们又提出了一些改进方法,如阻尼最小二乘法、迭代收敛法、共轭梯度法以及奇异值分解法等,在一定程度上解决了迭代过程的失稳和发散,提高了数值计算的稳定性。

2.2　相对计算方法

　　相对定位法是一种能有效减小因对地壳结构了解不够精细而引起误差的定位方法,也即2个事件至同一台站的走时差只由2个事件的相对位置以及它们之间的小范围内的波速所决定,与事件到台站的全部路程上的波速无关。

它曾用于研究与水库断层有关的浅源地震活动（Fitch,Muirhead,1974,也曾用于研究中

东地震的深度（Jackson,Fitch,1979,国内马秀芳（1982最早使用相对定位法,周仕勇等[7该方法做了改进。

为避免主地震误差传递到待定地震中,Gotetal（,Deichmann（1992等对于时空丛集的一群地震,,把这一方

法称为改进的主事件定位法。

Waldhauser等[8]（DD法,应

用于加州北海沃德断层上,还有Zhang,[9]。

JHD1971、HDC（EngdahlandBergman,

2000;Jordanand1983,它们都是处理多事件定位的方法,并都和illiamRodi等,2002。

Clifford等（2003在DD定位法发展了DD层析成像定位法,在实际应用中得到的结果优于DD法定位。

2.2.1　　　双重残差地震定位法由Waldhauser和Ellsworth于2000年提出,由每2个相邻地震的观测走时差减去

理论计算值的走时差得到的残差（“双差”

确定其相对位置:

ik9mmi-jk9m

mj=drijk（8式中,drijk是“双差”

drijk=（tik-tjkobj-（tjk-tjk

cal（9即第i个地震至第k个地震台站的地震波的走时tjk与第j个地震至第k个地震台站的地震波的走时tk之差

的观测值（tik-tjkobj与理论计算值（tik-tikcal的残差,Δmi=（Δxi,Δyi,Δzi,Δ

τiT,i=1,2,……,N构成的维数为4N的矢量,W是一用以对每个方程加权的对角线矩阵。

将由所有地震（i,j=1,2,…,N,所有台站（k=1,2,…得到的形式如式（8的方程用矩阵形式表示,便得到方程:

WGm=Wd

（10式中,G是一个M×4N的偏微商矩阵,M是双差观测的数目,N是地震数,d是双差资料矢量,m是由待定的

震源参数改变量（Δxi,Δyi,Δzi,Δ

τiT,i=1,2,…,N构成的维数为4N的矢量,W是一用以对每个方程加权的对角线矩阵。

反演中,引进了一个表示所有的地震经重新定位后其平均“位移”为0（也即其“矩心”不动的

约束条件:

∑N

i=IΔmi=0,以阻尼最小二乘法求解方程（10,此时问题归结为:

Wm=W（11

　　λ为阻尼因子,I为单位矩阵。

由正则方程可以得到方程（11的解为:

^m=（GTw-1G-1GTW-1d或者对

于数目不是很大的地震丛集,可以利用奇异值分解法得到方程（11的解:

^m=VpΛ-1pUTpd:

UP和VP分别为矩

阵G的两个正交奇异矢量矩阵G的2个正交奇异矢量矩阵,Λp是由G的奇异值构成的对角线矩阵。

实际计算采用共轭梯度法求解方程（11,得到阻尼最小二乘解,同时将奇异值分解法应用于部分资料以获得有

61　　　　　　地　震　工　程　与　工　程　振　动　　　　　　　　　　　　　　25卷

关模型参数的误差、分辨率等信息。

主事件法和双差法都是地震的相对定位方法,但它们的成立条件和应用范围是有区别的。

主事件法定位结果依赖于主事件的位置;双差法不要求主事件,结果不依赖于初始的震源位置。

主事件法只适用于较小空间范围的地震精确定位;双差法可用于大空间范围的地震精确定位。

杨智娴等[11]

利用双差定位法做了大量工作,并得到了很好的重定位结果;白玲等（2003.11利用多重相对定位法对1996年12月16日发生的顺义地震的余震震群进行重新定位;杨智娴等（2002.7还利用主事件定位法对1998年张北-尚义地震序列重新定位。

3　非线性计算法

　　非线性定位方法能够很好地解决在求目标函数极小值时避免其陷入局部极小点的问题。

3.1　Powell法

　　Powell方法[12]

是一种直接搜索目标函数极小值的有效方法,是一种改进的共扼梯度法。

如第Ti个地震

台记录到某次地震的P波到时的观测值为Ti,由给定的（初值或迭代值震源参数X（<、λ、h、t求出的理论

到时为τi（x,由此定义的目标函数为二次型函数

<（x=∑ni=I[T

i-τi（x]2（12

　　限值。

Powell方法不要计算导数,对初值的要求低,,它在台网的自动快速测报中有一定的优势。

3.2　遗传算法

　　遗传算法[13,14],。

具体应用:

参数（x,y,z,t可由相应的网格坐标（ix,iy,iz,it:

X=Xmin+iXΔX,这里Xmin为X参数搜索的下限,ΔX为X参数的网格单位,其余参数类推。

因此,一组整数（ix,iy,iz,it就可以表示为一震源位置。

首先在搜索范围内随机产生一组个体

（种群（xi,yi,zi,ti,i=1,2,…,Q,Q是种群中的个体总数,逐个计算每个模型对应各台观测资料的拟合差:

F（j=sta

i=IObi（i-Cal（i2　　j=1,2,…,Q（13

F（j是种群中第j个个体的拟合差,sta是观测值总个数,Obi（i是观测到时,Cal（i是计算到时,运用拟合差确定各模型的生存概率,个体生存概率将随拟合差的增大而减小,优胜劣汰。

根据个体的生存概率可采用轮盘赌方法,再选择Q个个体组成新的种群。

交配从种群中随机地选择两两一对的双亲,分别以交配概率随机地交换部分基因进行交配,由此产生2个新的后代代替原来的双亲。

以一定的变异概率使个体变异为新的个体,这样可以增加种群的多样性,从而能避免局部极值。

不断重复上述,最终得到满意的震源位置。

3.3　GMEL方法

　　WilliamRodi和E.RobertEngdah等[15,16]人介绍了一种新的网格多事件定位方法———GMEL法,并对其于JHD、HDC、PMEL及DD方法进行了理论比较和在1999年Izmit/Duzce地震序列的两个毗邻震群的实际应用上进行了比较研究,GMEL法是网格搜索求解,适合于多个台站和震相来自多个事件的地震到时、方位和观测的任意数据组。

一般只考虑包含到时的数据组,就是m个地震事件和n个地震台站到时,并分别用tj和xj表示第j个事件的发震时刻和震源位置。

假设每一对事件和台站至多只有一个观测到时,并且使用dij表示第i个台站和第j个事件,并可写成表达式:

dij=tj+Ti（xj+ci+eij（14

式中,Ti是假想地球模型下第i个台站的已知走时函数,ci是对真实地球和假想地球模型之间未校正的误差

函数,末知数xj、tj和ci。

ci有时使用台站校正代替路径校正;eij是一个观测误差;多事件定位问题就是解方程（14,求得带有路径校正的m个事件的定位参数（xj,tj。

非线性定位方法对初始位置的依赖性弱,不需要求偏导数或逆矩阵,避免解陷入局部极小点,但是运算量比较大,效率相对较低。

此方法还有网格搜索法、单纯形法、模拟退火法、蒙特卡罗法等。

这些定位方法很

711期杨文东等:

地震定位研究及应用综述

多作者都进行了实际应用,取得了一定效果;但是大部分方法还都在研究分析中,在实际地震台网观测定位中应用的主要有Powell法和遗传算法。

4　应用三维地壳结构的地震定位研究

以往人们研究地震定位问题都是基于一维的地壳模型,并通过采用不同的参数反演算法和对数据进行相应的加权处理等来提高地震定位的精度。

一维地壳模型由最初的J-B表发展到iasp91走时表和sp6模

型,Kennett等在此基础上提出了ak135模型,Engah[17]应用其对全球1900年~1999年的地震进行了重新定

位,明显改善了定位结果。

随着地震层析成像反演近年来的研究发展,三维地壳模型的建立（Nataf等的

3SMAC模型,Mooney等[18]的CRUST5.1模型,等等,为三维地震定位提供了基础。

Smith等[19]对已知的26

个核爆事件和83个有良好定位结果的地震事件,通过全球三维P波速SP12/WM13模型和一维的J-B、PREM、iasp91模型的定位对比分析,表明三维定位结果明显好于一维模型结果,定位偏差减少约40%。

An2tolik等也进行了三维结构模型地震定位的应用研究。

三维地震定位在地震成核、精细地质结构刻画、核查和震源位置的精度提高以及地震破裂的精细过程、EIDC建立的称为“地表真实”事件的核爆破目录和地震目录等方面都展示了相当大的应用研究潜力。

尽管三维地震定位研究在全球尺度目前还局限于应用探索阶段,区域尺度体现在三维结构的联合反演和特定震源区的应用研究中,应用于日常地震定位已为期不远。

5　台网观测系统实用地震定位法

　　的,,尤其是地震事件的智能化数值自动定位,已,为快速的地震参数报告、震后的地震应急和救灾、地震趋势的快速判断。

在台网定位中应用最广的是基于Geiger法发展起来的计算定位法,如HYPO系列和LOC系列。

目前国内有多种数字台网处理系统,在这些系统中都包含多种计算定位法,用于多种震相不同地震的定位。

5.1　国家数字地震台网中心的实时地震处理及大震速报系统

[20]国家数字台网中心大地震速报定位使用目标函数最小值搜索法和单台方位角法进行地震定位。

尤以搜索法为主,即在全空间范围内寻找一个点（x,y,z,根据最小二乘法,如果该点发生地震,地震波到达各台站的计算到时与实际观测到时差的平方和最小（与其他点相比,该点即为发生地震的位置。

目标函数如下:

F（x,y,z,t=

∑n

i=I（Ti-T′i（x,y,z,t2（15Ti为i地震台的观测到时,T′i为i地震台的计算到时。

目标函数使用J-B表计算理论到时,用P、S等多种

震相和7个深度的走时表。

搜索法中将全球分成网格,扫描网格点,寻找目标函数值最小的点,算出震中后再在地图上标出其位置,同时图上画出相应的一组双曲线或圆,通过观察震中位置与几何点的符合程度,确认计算定位结果的正确性,这就将计算法的精确和几何法的直观相结合,提高了初始定位的精度和可信度。

该定位系统可根据使用的震相或震中距选用不同的目标函数,这样可提高地震定位精度,另外需要时可增加使用其他震相的目标函数。

表1列出了使用不同震相的计算定位和对应的几何定位。

表1　地震定位方法和所用震相对照

使用震相

PS-PS-P,P方位角S-P,方位角定位

方

法搜索法∑ni=I∑nj=i+1（ΔTpij-Δtpiu2∑ni=I∑nj=i+1（ΔTspij-Δtspi2∑ni=1∑nj=i+I（ΔTpij-Δtpiu+ΔTspi-ΔTspi2几何法双曲线相交圆交切双曲线、圆相交直线相交圆和直线相交81　　　　　　地　震　工　程　与　工　程　振　动