1、资料内容仅供您学习参考,如有不当或者侵权,请联系改正或者删除。地震定位研究综述作 者 田 玥指导教师 陈晓非摘要 综述了各种地震定位方法的基本原理, 重点介绍了Geiger的经典方法以及在此基础上建立的各种线性方法: 联合定位法, 相对定位法, 和最新的双重残差法; 对每一种方法的应用情况, 特别是国内的工作做了总结; 同时也指出了各种方法的特点, 并进行了相应的比较。另外, 还简要介绍了空间域的定位方法和各种非线性定位方法。关键词: 震定位;线性定位;非线性定位目 录引言11经典定位方法11.1经典方法21.2各种改进方法32经典定位方法42.1震源位置与台站校正的联合反演(JED, JHD
2、)42.2 震源位置与速度结构的联合反演(SSH)52.3 相对定位法( 主事件定位法, ATD) 53空间域内的定位方法台偶时差法64非线性定位方法74.1 牛顿法74.2 全局搜索方法74.3 Bayesian方法85最新定位方法85.1 EHB方法85.2 双重残差法( DDA) 8结论9致谢9参考文献10引言地震定位是地震学中最经典、 最基本的问题之一, 对于研究诸如地震活动构造、 地球内部结构、 震源的几何构造等此类地震学中的基本问题有重要意义。另外, 基于快速准确的地震定位的地震速报, 对于震后的减灾、 救灾工作也是至关重要的。因此, 地震学家一直在不断改进或提出新的定位方法。地震
3、定位问题的提法如下: 根据台站对地震到时的观测资料, 来确定震源的空间坐标和发震时刻, 有时还给出对解的评价。早期的地震定位方法以几何作图法为主1。近三十年来由于计算机技术的飞速发展和广泛应用, 基于科学计算和计算机技术的智能化数值自动定位方法也得到了迅速发展, 并业已成为当前地震定位的主流方法。中国最初的地震定位工作由李善邦先生于1930年在北京鹫峰地震台开创, 1953年开始采用多台站大规模观测数据确定震中, 现在大多使用国际流行的定位方法。本文只介绍当前广泛使用的计算机定位方法, 重点介绍Geiger的经典方法以及在此基础上建立的各种线性方法: 联合定位法, 相对定位法, 和双重残差法,
4、 而且重点总结了国内的有关工作。1经典定位方法1.1经典方法现行的线性定位方法大都源于19 Geiger提出的经典方法2: 设n个台站的观测到时为 求震源及发震时刻, 使得目标函数 (1)最小。其中为到时残差 , (2)为震源到第i个台站的计算走时。 使目标函数取极小值也即 , (3)其中. 为方便, 记 , (4)则由(3)式, 在真解附近任意试探解及其校正矢量满足. (5)也即 . (6)由的定义可得公式(6)的具体表示式. (7) 若偏离真解不大, 则和较小, 可忽略二阶导数项, (7)式被简化为线性最小二乘解: . (8)以矩阵形式表示, 上式为 , (9)其中. 若二阶导数项不可忽略
5、, 则(7)式给出非线性最小二乘解 . (10) 一般各台站的到时数据具有不同的精度, 如果不加以区别, 则具有较低精度的数据将严重干扰结果的精度, 这一问题能够经过引入加权目标函数来解决。设各台站到时残差的方差为, 引入加权目标函数 , (11)按照上述同样的步骤, 经过求(11)式的极小值, 得到如下加权线性最小二乘解, (12)其中为加权方差矩阵: .由方程(9), (10), 或(12)求得后, 以作为新的尝试点, 再求解相应方程。如此重复迭代, 直至足够小( 或满足一定的循环结束条件) , 此时即得估计解。1.2各种改进方法直到20世纪70年代, 随着计算机的迅速兴起, Geiger
6、的思想才被广泛用于地震定位工作。Lee等人连续给出了HYPO71, HYPO7881系列程序3, 至今仍被普遍使用, 中国的赵仲和参与了80、 81版本程序的研制。Backus和Gilbert提出新的反演理论后, Klein提出HYPOINVERSE算法4, Lienert等在此基础上进一步得到HYPOCENTER算法5, Nelson和Vidale也改进了HYPOINVERSE, 提出了三维速度模型下的QUAKE3D方法6。在国内, 经典方法也得到了广泛应用: 赵仲和将HYPO81用于北京台网7, 吴明熙等8和赵卫明等9分别将经典方法用于禄劝地震和灵武地震序列的定位。针对求解基于Geiger
7、方法的线性方程组所遇到的各种问题, 许多学者提出了各种改进方法: 方程(9)的反演可有多种方法。例如当奇异或接近奇异时, 会引起迭代过程的失稳和发散, 此时能够采用奇异值分解( SVD) 求得估计解, 同时还可得到解的分辨率与误差估计; 当矩阵较大时, 能够采用共轭梯度法求解。 为了提高数值计算的稳定性, 一般采用中心化( centering) 、 定标化( scaling) 、 阻尼最小二乘法等方法5。使用最小二乘法( L2准则) 的前提是到时残差遵循Gauss分布, 但这一点常常得不到满足, 此时采用L1准则: , 可降低较大的到时残差的影响10。2经典定位方法经典方法是单事件定位方法。多
8、事件定位法联合定出多个震源以及其它参数( 如台站校正或速度模型) , 旨在解决用简单的速度模型代替复杂的地壳结构所引起的误差, 同时也提高了定位效率。2.1震源位置与台站校正的联合反演(JED, JHD)设有m个事件, n个台站。对每个台站j, 引入”台站校正”, 以弥补由速度模型简化所引起的误差。则对于事件i和台站j( i=1,2,m; j=1,2,n) , 有方程 , (13)其中为观测到时, 为事件i到台站j的计算走时, . 选定初始点和, 将(13)式做一阶Taylor展开, 可得到时残差, (14)设是到时残差的方差, 则可对上式加权: 。这样, 将(14)式用于所有事件和台站, 即
9、可联合反演出m个事件的震源位置及n个台站校正。1967年Douglas最先提出以上理论( JED) 11, 后来Dewey将其扩展成包括震源深度定位的JHD 方法12。为解决由于m, n过大而导致矩阵过大的问题, 1983年Pavlis和Booker提出参数分离的PMLE方法13, 并进一步被Pujol简化14-15。中国王椿镛等16根据昆明台网区域地震初至P波走时资料, 用JHD和参数分离法, 得到各台站P波走时的校正, 而且使定位精度有较大提高。2.2 震源位置与速度结构的联合反演(SSH)1976年, Crosson首次提出该联合反演理论17。由于SSH方法不需要对波速进行校准, 同时还
10、能够获得有关速度结构的很多信息, 是当前被广泛使用的一种定位方法。与JED方法相比, 该方法未引入台站校正, 而是将速度结构作为未知参数与震源同时反演, 由此解决人为构造的速度模型引起的误差。将(13)式改写为 , (15)其中是事件i到台站j的计算走时, 是一维速度模型矢量。给定初值, 将(15)式在该点作一阶Taylor展开可得到: , (16)将(16)式用于所有事件和台站, 即可联合反演出m个震源位置和速度模型。在一维速度结构与震源联合反演的理论基础上, Aki等人将地球内部横向非均匀速度结构网格化, 于1977年提出了三维速度结构与震源联合反演的理论18-19。可是用单一方程组联合反
11、演, 需要巨大的运算量, Pavlis和Booker20, Spencer和Gubbins21用参数分离法进行改进, 使耦合着的速度参数和震源参数分别求解, 大大提高了运算效率。在国内, 赵仲和于1983年建立了一个新的地震波速度模型MDBJ81, 以适应北京地区台网的稀疏分布, 并将该模型用于SSH方法, 提高了北京台网的测定能力22。刘福田引入正交投影算子实现参数分离, 并提出利用矩阵的块结构采取顺序正交三角化的方法, 减轻了运算量23。李强, 刘福田对SSH进一步改进, 应用最新的三维速度结构研究结果, 并考虑方程组的平衡问题以改进震源深度、 发震时刻的测定精度24。另外, 郭贵安等25
12、, 赵燕来等26, 朱元清等27分别将SSH方法用于震源的精确测定工作。2.3 相对定位法( 主事件定位法, ATD) 相对定位法由JED发展而来, 也是一个经典的、 被广泛采用的方法。Spence给出了该理论的详细阐述28。其基本原理是选定一震源位置较为精确的主事件, 计算发生在其周围的一群事件相对于它的位置, 进而计算这群事件的震源位置。设主事件为R, 其震源参数已知; 与R相距很近的待定事件为, 其震源参数为。由JED法列方程: 对事件R: , (17)对事件: , (18)将(18)式在点作一阶Taylor展开, 再与(17)式相减, 得到 . (19)这里引入了到时差( ATD) :
13、 . 由方程式(19)即可反演得到对R的相对位置, 于是可求得其震源参数。相对定位法经过引入到时差, 计算”相对位置”而消除了速度模型引起的误差, 有着独特的优点: 由于与R相距很近, 因此不需要迭代; 对主事件、 待定事件均不需要计算到时残差。该方法所得相对位置与相对到时的误差比经典方法小30%, 但绝对位置与绝对到时依赖于主事件。周仕勇等对该方法作了较大改进29: 定位中避开发震时刻的直接求解, 在确定震源后, 根据地震波的传播速度和距离计算, 而且采用首波到时资料专门确定深度。3空间域内的定位方法台偶时差法上述方法均为时间域内的定位方法, 基于对到时残差的处理, 4个震源参数彼此不完全独
14、立, 定位结果依赖于速度结构和台网分布。为了克服上述缺点, 人们同时提出了空间域内的定位方法, 即用距离残值代替到时残差, 方程只涉及震中位置, 震源深度和发震时刻单独求解, 避免了参数的相互折衷, 定位精度较高。Lomnitz30,Carza等31使用该方法进行远震定位, 震中误差为820km。1957年, Romney提出台偶时差近震定位法32, 利用到时相近、 位置相邻的两个台站( 即台偶) 的到时差和表面平均视速度来建立距离残差方程, 所得方程的条件数低, 易于求解, 而且定位结果对结构的依赖很少; 但对震源深度和发震时刻的确定并没有很好解决。赵珠, 曾融生对此做了改进33: 利用到时
15、曲线的斜率来确定震源深度, 利用到时曲线在时间轴上的截距来确定发震时刻。丁志峰, 曾融生对京津唐地区采用台偶时差法测定了震中, 震源深度的确定使用了不同震相间的到时差34。4非线性定位方法单事件与多事件定位法都是基于Geiger的线性方法, 该方法在很多情况下都会出现问题35。例如省略二阶以上的项不一定合理, 若选择不当, 线性迭代也会使解陷入局部极小点等。 非线性方法是处理这些问题的一个途径。4.1 牛顿法Thurber提出用包含二阶偏导数的非线性牛顿法来处理Geiger的方法所遇到的困难36。在均匀和分层速度模型中, 都存在以下问题: 对于震源深度接近于地表的浅震, 二阶偏导数趋于最大,
16、而一阶偏导数却趋于零, 此时二阶偏导数便极为重要; 另外, 当震源在台网之外时, 二阶偏导数的加入也提高了算法的稳定性。由于深度定位的不确定性来源于线性方法中发震时刻和震源深度的相关性, 而正确变化, 二阶偏导数比一阶偏导数更为敏感, 故它的引入减小了相关性, 提高了算法的稳定性。需要注意的是, 对于多事件定位或三维速度结构, 二阶偏导数的引入大大增加了计算量。Thurber根据牛顿法给出的非线性解为: , (20)此解亦为非线性最小二乘解 (11)。4.2 全局搜索方法非线性最优化理论中的各种全局搜索算法也被广泛的用于地震定位。用这些算法求目标函数的极小值, 能避免解陷入局部极小点, 而且正
17、确形式没有限制, 但计算效率一般都比较低。1988年Prugger和Gendzwill6, 1994年赵珠等37分别将单纯形法用于地震定位。单纯形法算法简单, 不需要求偏导数或逆矩阵, 但它不能给出解的分辨率和误差估计。1979年唐兴国将Powell直接搜索法用于地震定位38。该方法也不需要求偏导数或逆矩阵, 且对迭代初值要求较低, 选取到时最小的台站位置即可。Powell法本身不能给出误差估计, 汪素云等39利用对理论走时作随机扰动的数值实验法给出均方根残差与震源位置误差的关系。另外, 严尊国等40, 汪素云等41分别将此方法应用于三峡地区和青藏高原的重新定位。其它的全局搜索算法还有蒙特卡罗
18、法42, 模拟退火法, 遗传算法 43-44 等, 也都已被用于地震定位工作。4.3 Bayesian方法该定位方法根据Bayesian评估理论形成, 即从统计学的角度看, 模型参数的最佳值使观测数据的概率达到极大。二十世纪80年代, Tarantola和Valette45, Jackson和Matsuura46-47提出了Bayesian定位方法的严格公式及解。5最新定位方法5.1 EHB方法1999年, Engdahl等提出用于全球远震定位的EHB方法48。其中运用了多种震相: P, S, PKiKP, PKPdf, pP, pwP, sP, 而且改进得到了一种适用于后达波震相的全球速度模
19、型, 用以单独确认远震深度震相( pP, pwP, sP) 。与ISC, NEIC的定位结果相比, 用EHB重新定位得到的震中位置的精度明显提高, 而震源深度也得到显著改进。EHB方法可用于日常的快速定位, 也能够用于层析成像及其它地球内部结构的研究。5.2 双重残差法( DDA) 底, Waldhauser和Ellsworth提出了双重残差定位法49, 其基本原理简述如下。对台站k, 引入”事件对”i, j及双重残差 , (21)的计算能够用绝对到时, 也能够用两个事件的到时差。对单事件定位, 我们有 , (22)将(22)式分别用于事件i, j, 并相减, 得到 . (23)将(23)式用
20、于所有台站和事件对, 反演得到震源的绝对位置, 此即双重残差法(DDA).DDA的突出优点在于它能够利用谱域中的互相关分析法读取事件的到时差, 大大提高了到时数据的精确度。与相对定位法不同, 这里的事件正确距离不受限制, 很大程度上提高了该方法的适用性。若使用多种相位的到时差, 定位效果更为显著, 另外, 算法的抗干扰性、 健壮性也较强。从当前看, 这是一种很好的定位方法( 见图5.1) 。图5.1 DDA与其它定位方法对28个相关事件定位结果的比较49。结论从数学上讲, 地震定位问题的实质在于求目标函数的极小值。各种定位方法产生于对目标函数的构造、 处理, 以及求极小值方法的不同。影响地震定
21、位精度的主要因素有: 台网布局, 震相识别, 到时读数, 地壳结构等。在数值计算中, 常遇到下列问题: 走时的计算, 偏导数的计算, 方程的反演求解等。由于台网分布在地表, 给深度定位带来一定的困难。各种定位方法正是针对其中的某几个问题而设, 各有优、 缺点。相对定位所得的震源相对位置精度较高。对于主事件, 能够利用改进后的经典方法进行单事件定位。二者结合将能够得到较好的定位结果。JHD方法中引入的台站校正过于简单, 不足以反映地壳的复杂结构; 而SSH方法中的三维速度模型会带来巨大的运算量。如果我们能够构造一种介于二者之间的校正参数, 比如将台站校正作为有方向的矢量, 进行联合反演, 可能效
22、果更好。在DDA方法中, 当事件对i, j相距较近时, 能够将(23)式化简, 反演得到i, j的相对距离。同时我们能够选取较少的事件, 用联合反演进行绝对定位。将二者结合能够减少运算量, 提高定位效率。致谢 本文是在陈晓非老师的悉心指导下完成的。陈老师优秀的科学修养, 深厚的数理功底, 严谨的治学态度都给我留下了极其深刻的印象, 也成了我努力奋斗的榜样。本文同样凝聚了陈老师的心血, 仅此向陈老师表示深深的谢意。周仕勇博士后以其丰富的理论背景和实践经验, 对本文提出了建设性意见并提供了诸多及时的信息。张海明、 张伟、 邹最红、 曹军等师兄师姐手把手地领我入门, 令我受益匪浅。均在此一并致谢!参
23、考文献1 傅淑芳,刘宝诚. 地震学教程M. 北京:地震出版社,1991,447-480.2 Geiger L. Probability method for the determination of earthquake epicenters from arrival time onlyJ. Bull.St.Louis.Univ, 1912, 8: 60-71.3 Lee W H K, J C Lahr. HYPO71: A computer program for determining hypocenter, magnitude, and first motion pattern of l
24、ocal earthquakesJ. U.S.Geol.Surv. Open-File Rept, 1975, 75-311.4 Klein F W. Hypocenter location program HYPOINVERSE Part I: Users guide to versions 1,2,3 and 4J. U.S.Geol.Surv. Open-File Rept, 1978, 78-694.5 Lienert B R,Berg E, Frazer L N. Hypocenter: An earthquake location method using centered, sc
25、aled, and adaptively damped least squaresJ. Bull.Seism.Soc.Am, 1986,76(3): 771-783 .6 Nelson G D, John E Vidale. Earthquake locations by 3-D finite-difference travel timesJ. Bull.Seism.Soc.Am, 1990, 80(2): 395-410.7 赵仲和.多重模型地震定位程序及其在北京台网的应用J.地震学报,1983,5(2):242254.8 吴明熙,王鸣,孙次昌,等.1985年禄劝地震部分余震的精确定位J.地
26、震学报,1990,12(2):121129.9 赵卫明,金延龙,任庆维. 1988年灵武地震序列的精确定位和发震构造J.地震学报,1992, 14(4):416422.10 Prugger A F, Gendzwill D J. Microearthquake location: A nonlinear approach that makes use of a simplex stepping procedureJ. Bull.Seism.Soc.Am, 1988, 78(2): 799-815.11 Douglas A. Joint epicenter determinationJ. Nat
27、ure, 1976, 215: 45-48.12 Dewey J. Seismicity and tectonics of western VenezuelaJ. Bull.Seism.Soc.Am, 1972, 62(6): 1711-1751.13 Pavlis G, Booker J R. Progressive multiple event location (PMLE)J. Bull.Seism.Soc.Am, 1983,73(6):1753-1777.14 Pujol J. Comments on the joint determination of hypocenter and
28、station correctionsJ. Bull.Seism.Soc.Am, 1988, 78(3):1179-1189.15 Pujol J. Joint event location The JHD technique and applications to data from local seismic networksA. In: Thurber C, N Rabinowitz. Advances in seismic event locationC. Kluwer Academic Publishers, ,163-204.16 王椿镛,王溪莉,颜其中.昆明地震台网多事件定位问题
29、的初步研究J.地震学报,1993,15(2): 136145.17 Crosson R S. Crustal structure modeling of earthquake data,1,Simultameous least squares estimation of hypocenter and velocity parametersJ. J.Geophys.Res,1976,81(17):3036-3046.18 Aki K, Lee W H K. Determination of three-dimensional velocity anomalies under a seismic
30、array using first P arrival times from local earthquakes, part 1 : A homogeneous initial modelJ. J.Geophys.Res,1976,81(23):4381-4399.19 Aki K, et al. Determination of the three-dimensional seismic structure of the lothosphereJ. J.Geophys.Res,1977, 82(2): 277-296.20 Pavlis G, Booker J R. The mixed di
31、screte-continuous inverse problem : Application of the simultaneous determination of earthquake hypocenters and velocity structureJ. J.Geophys. Res, 1980,85(B9):4801-4810.21 Spencer C, Gubbins D. Travel-time inversion for simultaneous earthquake location and velocity structure determination in later
32、ally varying mediaJ. Geophys.J.Roy.Astr.Soc,1980, 63(1):95-116.22 赵仲和.北京地区地震参数与速度结构的联合测定J.地球物理学报,1983,26(2):131139.23 刘福田.震源位置和速度结构的联合反演( I) 理论和方法J.地球物理学报, 1984,27(2): 167175。24 李强,刘福田.一种横向不均匀介质中地震基本参数的测定方法J.中国地震,1991,7(3):54-63.25 郭贵安,冯锐.新丰江水库三维速度结构和震源参数的联合反演J.地球物理学报,1992, 35(3):331341。26 赵燕来,孙若昧,梅
33、世蓉.渤海地区地震参数的修定J.中国地震,1993, 9(2):129137。27 朱元清,范长青,浦小峰.南黄海地震序列时空参数的精细测定和分析J.中国地震, 1995, 11(1):5461。28 Spence W. Relative epicenter determination using P-wave arrival-time differencesJ. Bull. Seism.Soc.Am, 1980,70(1):171-183.29 周仕勇,许忠淮, 韩京等.主地震定位方法分析以及1997年新疆伽师震群高精度定位J.地震学报,1999,21(3):258-265。30 Lomni
34、tze C. A fast epicenter location programJ. Bull.Seism.Soc.Am, 1977, 67(2): 425-431.31 Carza T, Lomnitz C, C Ruiz de velasco. An interactive epicenter location procedure for the RESMAC seismic arrayJ. Bull.Seism.Soc.Am,1979,69(4): 1215-1236.32 Romney C. Seismic waves from the Dixie Valley Fairview Pe
35、ak earthquakesJ. Bull.Seism. Soc.Am,1957, 47(4): 301-319.33 赵珠,曾融生.一种修定震源参数的方法J.地球物理学报,1987,30(4):379-388。34 丁志峰,曾融生.京津唐地区震源深度分布初探J.地震学报,1990,12(3):242-247。35 Lee W H K, Stewart S W. Principles and applications of microearthquake networksM. Academic Press.NewYork,. Adv.Geophys.Suppl.2, 1981.36 Thurb
36、er C H. Nonlinear earthquake location : theory and examplesJ. Bull.Seism.Soc.Am, 1985,75(3):779-790.37 赵珠,丁志峰,易桂喜,等, 西藏地震定位一种使用单纯形优化的非线性方法J.地震学报,1994,16(2):212219.38 唐国兴.用计算机确定地震参数的一个通用方法J.地震学报,1979,1(2):186196.39 汪素云,许忠淮,俞言祥,等.北京西北地区现代微震重新定位J.地震学报,1994,16(1):24-31.40 严尊国,薛军蓉.长江三峡地区弱震重新定位J.中国地震,1987
37、,3(1):5259。41 汪素云,高阿甲,许忠淮,等, 青藏高原东北地区地震重新定位及其活动特征J.地震学报, ,22(3):241-248.42 Billing SD, Sambridge MS, Kennet BLN. Errors in hypocenter location picking, model, and magnitude dependenceJ. Bull.Seism.Soc.Am, 1994, 84(6):1978-1990.43 Xie Z, Spencer TW, Rabinowitz PD, et al. A new regional hypocenter loc
38、ation methodJ. Bull.Seism.Soc.Am, 1996, 86(4):946-958.44 Yong-Ge Wan, Rui-Feng Liu, Hong-Ji Li. The inversion of 3-D crustal structure and hypocenter location in the Beijing-Tianjin-Tangshan-Zhangjiakou area by genetic algorithmJ. Acta Seismological Sinica, 1997, 10(6):769-781.45 Tarantola A, Valett
39、e B. Inverse problem = quest for informationJ. J.Geophys, 1982, 50: 159-170.46 Matsuura M. Bayesian estimation of hypocenter with origin time eliminatedJ. J.Phys. Earth, 1984, 32(6):469-483.48 Jackson D D, Matsuura M. A Bayesian approach to nonlinear inversionJ. J.Geophys. Res, 1985,90(B1):581-591.4
40、9 Engdahl E R, Rob van der Hilst, Raymond Buland. Global teleseismic earthquake relocation with improved travel times and procedures for depth determinationJ. Bull.Seism.Soc.Am, 1998,88(3):722-743.50 Waldhauser F, Ellsworth W L. A double-difference earthquake location algorithm : method and application to the Northern Hayward Fault, CaliforniaJ. Bull.Seism.Soc.Am, , 90 (6):1353-1368.