利用观测到时确定震源位置和发震时间,已成为现代地震台网的常规任务,但地震定位精度仍然有待提高.准确的地震定位不仅可为地震活动性分析和地球深部结构研究(地震层析成像)提供基础,还可用于推断深部的断裂几何形态和介质物理性质[1-3].确切事件的定位实验表明,地震台网的常规地震定位结果还存在较大的不确定性[4].对于相同的地震序列,定位结果也常常存在着显著的差异[5].地震定位的精度除了受地震台网分布、走时拾取精度、速度模型准确性以及所用的震相种类等因素影响外[3, 6],还与定位方法密切相关[2].为了提高地震定位精度,地震学家发展了多种定位方法,如盖格法[7]、双差法[8]、网格搜索法[9]以及交切法[10-12]等.这些方法各有特点,其中交切法不仅直观而且效率高、稳健性强,在地震台网中有广泛的应用[13-14].交切法利用震源轨迹进行定位,不需要求解方程,即使仅有少量的地震记录,利用交切法也能获得有价值的震中信息.但该方法定位精度较低,特别是震源深度[10].因此,目前该方法主要作为辅助的定位手段.交切法定位精度较低的原因主要是它假定地球模型为均匀或水平均匀介质,而实际上地球内部在径向上和横向上都存在较强的非均匀性[15-20].使用远离实际、过于简化的速度模型必然导致或大或小的定位误差[21].因此,提高传统交切法定位精度的关键是使之基于接近实际的速度模型.
基于复杂速度模型的交切法,其难点在于震源轨迹的计算.当速度结构复杂时,震源轨迹难以表达成解析形式.Zhao & Ding[22]针对二维的复杂速度模型,提出了利用最小走时树射线追踪技术以离散形式准确计算震源轨迹的解决方案.数值实验表明,其改进的交切法可用于复杂速度模型的地震定位.在本项研究中,我们将Zhao & Ding 改进的交切法[22]推广到三维复杂速度模型.我们首先简要介绍三维复杂速度模型的最小走时树射线追踪方法,然后介绍复杂速度模型的地震定位方法.最后通过准确速度模型、扰动速度模型、扰动观测到时及地震在台网外等4种情况地震定位,检验了改进交切法的可行性.
2 最小走时树射线追踪方法最小走时树射线追踪方法基于地震学的惠更斯原理和费玛原理.该方法通过构建最小走时树完成地震波从震源至模型所有节点走时的计算,其基本算法[23]来自网络理论中的Dijkstra方法[24].最小走时树法具有良好的稳健性,适用于复杂速度模型,可计算多种地震波的走时及射线路径.该方法自20世纪80年代中期提出后,不断得到发展和完善[25-36].赵爱华等[37-39]基于动态子波传播区域的思想提出的改进算法不仅计算效率高,而且继承了基本算法稳健、简单的优点.其算法为:
将三维速度模型离散化为大小相等的立方体单元,各离散单元内速度均匀.模型节点位于单元的中心.设N为模型节点的集合,EDP(i)为点i∈ N有效方向点的集合,j∈EDP(i)的可能子波源点集合为PS(i,j),PS(i,j)包括i、i的子波源点s(i)[28]以及s(i)周围的26个节点.子波从i∈ N到点j∈EDP(i)的走时为dt(i,j).初始时,地震波激发点O处的初至波走时t(O)为零.随着地震波在介质中的传播,集合N中的点可分为三部分:走时已达到最小点的集合P,走时尚未达到最小但已更新过点的集合Q和其余点构成的集合R.计算从O点到点i∈ N初至波走时的过程为[39]
(1) 初始化
(2)查询
在Q中查询最小走时t(i)(i∈Q)
(3) 更新
(4) 迭代条件检测
如果Q=Φ,则结束;否则回到(2).其中,A∶=B表示将B赋值于A,Φ 表示空集合.使用不同的地震波速度,可计算相应不同地震波的走时和射线路径.
3 复杂速度模型的地震定位方法复杂速度模型的地震定位方法借鉴了传统交切法的思想,即利用震源轨迹定位.当速度结构复杂时,震源轨迹通常难以解析式表示,为此对其在空间域中以离散方式进行表示.在传统的交切法中[10-14],震源轨迹用观测到时差进行约束,以避开发震时间问题,本文亦采用这种做法[22, 40-41].
假定某地区地下速度分布已知,在地表有M地震台站,Ri(i=1,2,3,…,M).为对区内的地震进行定位,将速度模型离散成大小相等的立方体单元.在模型单元内,地震波速度为常数.利用最小走时树射线追踪方法[39]分别计算在台站Ri(i=1,2,3,…,M)激发的纵波(P波)和横波(S波)到节点x= (x,y,z)的初至波走时,T(x;Ri,P)和T(x;Ri,S)(i=1,2,3,…,M).这样,每2 个走时场可构成1个走时差场:
(1) |
其中,Rj1 和Rj2 分别表示对应第j个走时差场的2个走时场的地震台站;Wj1 和Wj2 分别表示波的类型,即P或S;K为走时差场数.对于区域内的某个地震,震源在x0 =( x0,y0 ,z0),发震时间为t0,台站Ri记录到的P 波、S波的初至到时分别为τRiP和τRiS,其震源轨迹满足如下方程:
(2) |
理论上,震源轨迹将交切于一点,即震源点.但实际上地质模型和初至波走时都会有或多或少的误差,因此震源轨迹将交汇成一个小的区域,该区域的大小反映了定位精度.在本文中,我们将震源定位于轨迹交会最为集中的点x0* = (x0* ,y0* ,z0*),该点使如下的目标函数值最小:
(3) |
其中,I(x)=1.x0* 点可通过扫描整个模型节点空间找到.
震源位置确定后,发震时间t0* 可定为
(4a) |
或
(4b) |
地震定位实践中,通常以理论到时与实际到时之差RAT 评估定位效果,RAT 的表达式为
(5a) |
或
(5b) |
对于交切法,地震定位结果可根据震源轨迹的交汇情况进行定性评估.若震源轨迹在震源附近沿某个方向交汇较多,则意味着在该方向上对震源位置约束较弱.震源轨迹交汇越密集的节点处,总的双重时差(观测到时差与理论到时差之绝对差)RDT越小.因此,具有较小RDT 值的点的空间分布就反映了震源轨迹的交汇情况,据其可定性地估计定位结果的不确定性.为便于叙述,将RDT 值最小的p%的模型节点称为Fp点.
4 模型算例 4.1 速度模型概述
随着地壳结构的地球物理探测,人们可得到较精细的地壳结构模型[17-20].图 1 显示了一个尺度为1km×1km×1km 的复杂速度模型[39].该模型由2层均匀介质与1层非均匀介质组成,其中第一个界面为水平界面:z= - 0.3km; 第二个界面为倾斜界面,其空间位置满足关系x+y+10z+6.0=0;在这两个界面之间嵌有一椭球体E,中心点坐标为(0.5km, 0.5km, -0.5km),其形状满足方程
单位为:km/s.
第二界面以下介质Ⅲ为均匀介质,其地震波速度为vP=3.5km/s, vS=2.1km/s.椭球体内的介质速度为:vP=3.0km/s, vS=1.7km/s.将模型离散成10 m×10 m×10 m 的立方体单元.在地表布设9 个地震台站Ri,其位置为(215+ (m-1)×300 m, 215+ (n-1)×300 m, -5 m),其中i=3(n-1)+m(m,n=1,2,3).假定某个地震发生在点x0(315m, 735m, -655m),即图 1中红色“+"处,发震时间t0=0.0ms.由震源x0 至地震台站Ri的P、S波到时τRiP和τRiS利用改进的最小走时树射线追踪方法[39]计算,子波传播区域半径取1.
对于图 1中的地震,参照Gajewski & Tessmer的做法[40-42],考察如下4种情况对地震定位的影响:(1)准确速度模型、(2)扰动速度模型、(3)扰动观测到时及(4)地震在台网外.激发点在地震台站处的地震波的初至走时使用最小走时树改进算法[39]计算,子波传播区域半径取1.选取RDT 值最小的0.05%的模型节点反映震源位置的约束状况,即Fp点占总的模型节点的百分比p=0.05.
在这个试验中,模型和观测到时都是准确的.我们设计了两个方案A 和B 进行定位.方案A 使用12条以不同台站P 波到时差约束的震源轨迹(I类震源轨迹[40]),对应的走时差场为DT(x;R1,P;Ri,P)、DT(x;R3,P;Ri,P)(i=7、8和9)和DT(x;R5,P;Ri,P)(i=1、3、4、6、7和9).方案B 除了使用方案A 中的震源轨迹外,还使用9条相同台站P、S波到时差约束的震源轨迹(II类震源轨迹[40]),其相应的走时差场为DT(x;Ri,P;Ri,S)(i= 1、2、3…9).采用方案A 和B的定位结果分别显示在图 2(a, b)中.两种方案确定的震源位置(蓝十字)和实际震源位置(红十字)完全一致,计算的发震时间也正确,t0*=0.0ms.在图 2a中,F0.05点的分布形状近似为细长的椭球.这说明I类震源轨迹在震源附近较集中于椭球的长轴方向分布,在该方向上,震源位置约束得不是很好.在图 2b 中,F0.05 点的分布近似为球状,表明综合使用两类震源轨迹能对震源位置有更好的约束.
4.3 扰动速度模型通常,地下速度结构是难以精确获知的.为模拟速度模型有误差的情况,对图 2模型P、S波的速度分别加幅度为其平均值±10% 的随机扰动(扰动利用FORTRAN90 编译器中的随机函数生成),均方速度扰动为5.8%,最大相对速度扰动量为13.3%.这样,图 2模型变为随机介质模型.采用方案A 和B的定位结果分别如图 2(c, d)所示.在图 2c中,确定的震源位置在点(315m, 745m, -635m)处,在x、y和z轴方向的偏差分别为0m、10m 和20m; 发震时间滞后1.1ms, 即t0* =1.1ms.在图 2d中,确定的震源位于点(285m, 765m, -685m)处,定位误差在三个坐标轴方向上均为30m; 确定的发震时间t0* =2.1ms.从F0.05点的分布来看,方案B 的震源约束优于方案A,但后者的定位误差略大于前者.
4.4 扰动观测到时为考察到时拾取误差对震源轨迹分布的影响,对观测到时加幅度为±10 ms的随机扰动(扰动利用FORTRAN90 编译器中的随机函数生成),P、S波的均方到时扰动量分别为7.1ms和6.4ms.使用具有观测误差的到时进行定位,对应方案A和B的定位结果分别如图 2e和2f所示.从图 2e可以看出,方案A 确定的震源位置分别向左、向后和向下偏离实际震源位置80m、和30m 和130m; 相应的发震时间超前-39.2ms, 即t0* =-39.2ms.值得注意的是,定位偏差方向和F0.05 点椭球的长轴方向一致.在图 2f中,确定的震源位置在实际位置的左前下方,三个坐标轴方向上的距离分别为30 m、20 m和20m; 确定的发震时间也更接近实际值,仅超前9.0ms.
在实际中,地震有时在地震台网之外.为模拟这种情况,将台站R1、R2、R4 和R7 去掉.用与方案A、B类似的方案A'、B'进行定位.方案A'使用5 条I类震源轨迹,其对应的走时差场为DT(x;R3,P;R8,P)和DT(x;R5,P;Ri,P)(i=3、6、8和9).方案B'使用10条震源轨迹,除了方案A'中的震源轨迹外,还包括5条II类震源轨迹,其相应的走时差场为DT(x;Ri,P;Ri,S)(i=3、5、6、8和9).采用方案A'、B'的定位结果分别显示在图 2g和2h 中.可以看出,虽然两种方案都能得到正确的震源位置和发震时间,但由于地震在台网外,仅用P 波到时信息很难控制好震源位置,特别是震源深度,如果同时利用P波与S波之间的到时差信息,则震源位置的约束明显改善,这是因为对于同一个台站,P 波与S波的传播路径基本相同,其到时差对于震源位置较为敏感.
5 讨论和结论改进的交切法不再假设震源轨迹为球面或双曲面,而是利用具有良好稳健性的最小走时树射线追踪技术准确计算,因而可用于复杂速度模型的地震定位.在改进的方法中,计算在台站处激发的地震波走时场比较费时,但仅需计算1次.对于速度结构稳定地区、同一地震台网记录的多个地震,改进方法具有较高的效率.总的双重时差值较小的点,其空间分布定性地反映了定位结果的不确定性.
数值模型的定位实验表明:以不同台站P 波到时差约束的震源轨迹由于空间分布方位较窄,对震源的位置约束不是很强,特别是当地震在台网之外时;定位时,利用相同台站P、S 波到时差约束的震源轨迹可显著增强对震源位置的约束.使用多条震源轨迹进行定位,有助于减少由随机因素导致的定位误差;震源位置约束得越好,并不总意味着定位误差会更小.
致谢中国地震局地球物理研究所丁志峰研究员和中国科学院地质与地球物理研究所张中杰研究员对本文初稿进行了认真审校和修改,在此表示衷心感谢.同时,作者也非常感谢评审专家的中肯的评论.
[1] | 田玥, 陈晓非. 地震定位研究综述. 地球物理学进展 , 2002, 17(1): 147–155. Tian Y, Chen X F. Review of seismic location study. Progress in Geophysics (in Chinese) , 2002, 17(1): 147-155. |
[2] | 蔡明军, 山秀明, 徐彦, 等. 从误差观点综述分析地震定位方法. 地震研究 , 2004, 27(4): 314–317. Cai M J, Shan X M, Xu Y, et al. Review of earthquake-locating methods from error. Journal of Seismological Research (in Chinese) , 2004, 27(4): 314-317. |
[3] | 杨文东, 金星, 李山有, 等. 地震定位研究及应用综述. 地震工程与工程振动 , 2005, 25(1): 14–20. Yang W D, Jin X, Li S Y, et al. Study of seismic location methods. Earthquake Engineering and Engineering Vibration (in Chinese) , 2005, 25(1): 14-20. |
[4] | Storchak D A. Results of locating the IASPEI GT (0-5) reference events using the standard ISC procedures. Phys. Earth Planet. Inter. , 2006, 158(1): 4-13. DOI:10.1016/j.pepi.2006.03.003 |
[5] | 黄媛, 吴建平, 张天中, 等. 汶川8.0级大地震及其余震序列重定位研究. 中国科学(D辑) , 2008, 38(10): 1250–1257. Huang Y, Wu J P, Zhang T Z. Relocation of the Ms8.0 Wenchuan earthquake and its aftershock sequence. Science in China (Series D) (in Chinese) , 2008, 38(10): 1250-1257. |
[6] | Billings S D, Sambridge M S, Kennett B L N. Errors in hypocenter location: picking, model, and magnitude dependence. Bull. Seism. Soc. Am. , 1994, 84(6): 1978-1990. |
[7] | Geiger L. Probability method for the determination of earthquake epicenters from the arrival time only. Bull. St. Louis Univ. , 1912, 8: 60-71. |
[8] | Waldhauser F, Ellsworth W L. A double-difference earthquake location algorithm: method and application to the northern Hayward fault, California. Bull. Seismol. Soc. Am. , 2000, 90(6): 1353-1368. DOI:10.1785/0120000006 |
[9] | Shearer P M. Improving local earthquake location using the L1 norm and waveform cross correlation: Application to the Whittier Narrows, California, aftershock sequence. J. Geophys. Res. , 1997, 102(B4): 8269-8283. DOI:10.1029/96JB03228 |
[10] | Udías A. Principles of Seismology. Cambridge: Cambridge University Press, 1999 . |
[11] | Pujol J. Earthquake location tutorial: graphical approach and approximate epicentral location techniques. Seis. Res. Lett. , 2004, 75(1): 63-74. DOI:10.1785/gssrl.75.1.63 |
[12] | 廉超, 李胜乐, 董曼, 等. 球面交切法地震定位. 大地测量与地球动力学 , 2006, 26(2): 99–103. Lian C, Li S L, Dong M, et al. Method of spherical surface shearing in earthquake location. Journal of Geodesy and Geodynamics (in Chinese) , 2006, 26(2): 99-103. |
[13] | 孟玉梅, 赵永, 王斌. 中国地震观测台网地震速报定位偏差的分析. 地震. 2001 : 65 -69. Meng Y M, Zhao Y, Wang B. The analysis on deviation of rapid location by China Seismic Observational Network. Earthquake (in Chinese). 2001 : 65 -69. |
[14] | 宋锐, 顾小红, 王永力, 等. 国家数字地震台网中心实时处理及大地震速报软件系统. 地震 , 2001, 21(4): 47–59. Song R, Gu X H, Wang Y L, et al. The software system for real-time processing and large earthquake rapid determination of NCDSN. Earthquake (in Chinese) , 2001, 21(4): 47-59. |
[15] | 张中杰, 滕吉文, 张霖斌, 等. 地球深部研究的复杂性及其对策. 地球物理学进展 , 1996, 11(3): 10–20. Zhang Z J, Teng J W, Zhang L B, et al. The complexity in the study of earth depth and its countermeasure. Progress in Geophysics (in Chinese) , 1996, 11(3): 10-20. |
[16] | Zhang Z J, Badal J, Li Y K, et al. Crust-upper mantle seismic velocity structure across southeastern China. Tectonophysics , 2005, 395(1-2): 137-157. DOI:10.1016/j.tecto.2004.08.008 |
[17] | Zhang Z J, Klemperer S L. West-east variation in crustal thickness in northern Lhasa block, central Tibet, from deep seismic sounding data. J. Geophys. Res. , 2005, 110(B9): B09403. DOI:10.1029/2004JB003139 |
[18] | Zhang Z J, Yang L Q, Teng J W, et al. An overview of the earth crust under China. Earth Science Reviews , 2011, 104(1-3): 143-166. DOI:10.1016/j.earscirev.2010.10.003 |
[19] | Zhang Z J, Deng Y F, Teng J W, et al. An overview of the crustal structure of the Tibetan plateau after 35 years of deep seismic soundings. Journal of Asian Earth Sciences , 2011, 40(4): 977-989. DOI:10.1016/j.jseaes.2010.03.010 |
[20] | 滕吉文, 王夫运, 赵文智, 等. 鄂尔多斯盆地上地壳速度分布与沉积建造和结晶基底起伏的构造研究. 地球物理学报 , 2008, 51(6): 1753–1766. Teng J W, Wang F Y, Zhao W Z, et al. Velocity distribution of upper crust,undulation of sedimentary formation and crystalline basement beneath the Ordos basin in North China. Chinese J.Geophys. (in Chinese) , 2008, 51(6): 1753-1766. |
[21] | Chen H, Chiu J M, Pujol J, et al. A simple algorithm for local earthquake location using 3D VP and VS models: test examples in the central United States and in central eastern Taiwan. Bull. Seismol. Soc. Am. , 2006, 96(1): 288-305. DOI:10.1785/0120040102 |
[22] | Zhao A H, Ding Z F. An intersection method for locating earthquakes in complex velocity models. Applied Geophysics , 2007, 4(4): 294-300. DOI:10.1007/s11770-007-0036-5 |
[23] | Nakanishi I, Yamaguchi K. A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure. J. Phys. Earth , 1986, 34(2): 195-201. DOI:10.4294/jpe1952.34.195 |
[24] | Dijkstra E W. A note on two problems in connection with graphs. Numer. Math. , 1959, 1(2): 269-271. |
[25] | Moser T J. Shortest path calculation of seismic rays. Geophysics , 1991, 56(1): 59-67. DOI:10.1190/1.1442958 |
[26] | Cao S H, Greenhalgh S. Calculation of the seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm. Geophys. J. Int. , 1993, 114(3): 593-600. DOI:10.1111/gji.1993.114.issue-3 |
[27] | Fischer R, Lees J M. Shortest path ray tracing with sparse graphs. Geophysics , 1993, 58(7): 987-996. DOI:10.1190/1.1443489 |
[28] | Klimes L, Kvasnicka M. 3-D network ray tracing. Geophys. J. Int. , 1994, 116(3): 726-738. DOI:10.1111/gji.1994.116.issue-3 |
[29] | 刘洪, 孟凡林, 李幼铭. 计算最小走时和射线路径的界面网全局方法. 地球物理学报 , 1995, 38(6): 823–832. Liu H, Meng F L, Li Y M. The interface grid method for seeking global minimum travel-time and the correspondent raypath. Chinese J. Geophys. (in Chinese) , 1995, 38(6): 823-832. |
[30] | Cheng N, House L. Minimum traveltime calculation in 3-D graph theory. Geophysics , 1996, 61(6): 1895-1898. DOI:10.1190/1.1444104 |
[31] | 王辉, 常旭. 基于图形结构的三维射线追踪方法. 地球物理学报 , 2000, 43(4): 534–541. Wang H, Chang X. 3-D ray tracing method based on graphic structure. Chinese J. Geophys. (in Chinese) , 2000, 43(4): 534-541. |
[32] | 赵爱华, 张中杰, 王光杰, 等. 非均匀介质中地震波走时与射线路径快速计算技术. 地震学报 , 2000, 13(2): 167–173. Zhao A H, Zhang Z J, Wang G J, et al. A new scheme for fast calculation of seismic traveltimes and ray paths in heterogeneous media. Acta Seismologica Sinica (in Chinese) , 2000, 13(2): 167-173. DOI:10.1007/s11589-000-0006-y |
[33] | 张建中, 陈世军, 徐初伟. 动态网络最短路径射线追踪. 地球物理学报 , 2004, 47(5): 899–904. Zhang J Z, Chen S J, Xu C W. A method of shortest path raytracing with dynamic networks. Chinese J. Geophys. (in Chinese) , 2004, 47(5): 899-904. |
[34] | 张美根, 贾豫葛, 王妙月, 等. 界面二次源波前扩展法全局最小走时射线追踪技术. 地球物理学报 , 2006, 49(4): 1169–1175. Zhang M G, Jia Y G, Wang M Y, et al. A global minimum traveltime raytracing algorithm of wavefront expanding with interface points as secondary sources. Chinese J. Geophys. (in Chinese) , 2006, 49(4): 1169-1175. |
[35] | 张美根, 程冰洁, 李小凡, 等. 一种最短路径射线追踪的快速算法. 地球物理学报 , 2006, 49(5): 1467–1474. Zhang M G, Cheng B J, Li X F, et al. A fast algorithm of shortest path ray tracing. Chinese J. Geophys. (in Chinese) , 2006, 49(5): 1467-1474. |
[36] | Bai C, Greenhalgh S, Zhou B. 3D ray tracing using a modified shortest-path method. Geophysics , 2007, 72(4): T27-T36. DOI:10.1190/1.2732549 |
[37] | 赵爱华, 张中杰, 彭苏萍. 复杂地质模型转换波快速射线追踪方法. 中国矿业大学学报 , 2003, 32(5): 513–516. Zhao A H, Zhang Z J, Peng S P. Fast ray tracing method for converted waves in complex media. Journal of China University of Mining & Technology (in Chinese) , 2003, 32(5): 513-516. |
[38] | Zhao A H, Zhang Z J, Teng J W. Minimum travel time tree algorithm for seismic ray tracing: improvement in efficiency. J. Geophy. Eng. , 2004, 1(4): 245-251. DOI:10.1088/1742-2132/1/4/001 |
[39] | 赵爱华, 张中杰. 三维复杂介质中转换波走时快速计算. 地球物理学报 , 2004, 47(4): 702–707. Zhao A H, Zhang Z J. Fast calculation of converted wave traveltime in 3-D complex media. Chinese J. Geophys. (in Chinese) , 2004, 47(4): 702-707. |
[40] | 赵爱华, 丁志峰, 孙为国, 等. 复杂介质地震定位中震源轨迹的计算. 地球物理学报 , 2008, 51(4): 1188–1195. Zhao A H, Ding Z F, Sun W G, et al. Calculation of focal loci for earthquake location in complex media. Chinese J. Geophys. (in Chinese) , 2008, 51(4): 1188-1195. |
[41] | Zhao A H, Ding Z F. Earthquake location in transversely isotropic media with a tilted symmetry axis. J. Seismology , 2009, 13(2): 301-311. DOI:10.1007/s10950-008-9129-8 |
[42] | Gajewski D, Tessmer E. Reverse modelling for seismic event characterization. Geophys. J. Int. , 2005, 163(1): 276-284. DOI:10.1111/gji.2005.163.issue-1 |