基于程函方程的初至P和sPg波走时联合地震定位方法

李孟洋, 刘少林, 潘阳, 杨顶辉, 申文豪,3, 徐锡伟

1 中国科学院大学地球与行星科学学院, 北京 1000492 应急管理部国家自然灾害防治研究院, 北京 1000853 上海佘山地球物理国家野外科学观测研究站, 上海 2016024 清华大学数学科学系, 北京 1000845 中国地质大学(北京)地球科学与资源学院, 北京 100083

0 引言

地震震源位置的精确测定是地震学的基本问题之一,精确的震源位置是研究发震构造和地震危险性的关键参数(崇加军等,2010;Ma and Eaton,2011;韩立波等,2012;李大虎等,2015;Diehl et al.,2021;王未来等,2021;王烁帆等,2022),也是地震层析成像的基础(Zhang and Thurber,2003;Liu et al.,2019).为了得到准确震源位置,前人开展了大量研究(Ma and Eaton,2011).由于初至P波震相易拾取,因此初至P波震相反演被广泛用于地震定位中(曾宪伟等,2019).由于初至P波震相对震源深度敏感程度较低,以至于其难以约束震源深度(Yuan et al.,2020;Rajkumar et al.,2022).尤其在台站稀疏区域,由于数据不足,因此以初至波到时定位反演方法约束震源深度尤为困难(Ma and Eaton,2011;罗艳等,2010).

为了提高震源深度的测定精度,前人提出利用深度震相(例如:sPL、sPg、sPn和sPmP)约束震源深度(Langston,1987,1994;Murphy and Barker,2006;韦生吉等,2009;罗艳等,2010;崇加军等,2010;孙茁等,2014).深度震相是从震源发出的体波先在震中附近界面(例如:地表)反射后再到达观测点的震相(Ma,2010;Yuan et al.,2020).由于深度震相到时对震源深度更敏感,所以被广泛用于测定震源深度(Langston,1987;Bock et al.,1996;Ma and Atkinson,2006;崇加军等,2010;罗艳等,2010;韩立波等,2012).最早Langston(1987)采用sPg和Pg震相到时差研究1968年澳大利亚Meckering近震序列的深度分布.在Langston(1987)的研究中,首先采用层状模型近似研究区域的速度结构,并使用波数积分方法(Wave Number Integration Method;Apsel,1979;Barker,1984)计算震中距为60～95 km范围内台站的合成地震图,然后参考合成地震图拾取sPg和Pg震相到时差.震源深度由sPg和Pg震相到时差除以P波和S波慢度之和得到.Langston(1987)使用的震源深度计算方法,只有在震中距小、地震波近垂直入射时才有效,随着震中距增大,由计算引入的震源深度定位误差会逐渐增大.Ma和Atkinson(2006)提出一种深度震相模拟方法(Regional Depth Phase Modeling,RDPM)测定中小地震(MN≥2.8)的震源深度.在该方法中,不使用近似公式直接计算震源深度,避免了由近似公式计算引入的震源深度定位误差,通过对比理论地震图与实际地震图的sPg与Pg震相到时差获得震源深度,在计算理论地震图时,研究区域内的速度结构使用层状模型近似.

除了sPg震相外,深度震相sPL、sPmP和sPn(图1)也被频繁用于测定震源深度的研究(Ma and Atkinson,2006;崇加军等,2010;罗艳等,2010).不同的深度震相有不同的优势震中范围,定位时应根据震中距合理选择相应的深度震相．当震中距小于50 km时,sPL震相通常是地震图上最清晰的深度震相(Yuan et al.,2020).崇加军等(2010)和罗艳等(2010)通过对比基于层状模型模拟的合成地震图与实际地震图中sPL和P波走时差,确定震中距在50 km以内的目标地震的震源深度.对于震中距约在180～350 km之间的地震,sPmP震相相比于其他深度震相在地震图上发育更好,Bock等(1996)、Saikia等(2001)、Kim等(2006)和Kastrup等(2007)通过对比模拟波形与实际波形中的sPmP震相标定此震中距范围内地震的震源深度.在300 km左右的震中距范围内,还可使用sPn震相测定震源深度,但sPn震相振幅较小,只有中强震可以产生较明显的sPn震相(Ma and Atkinson,2006;崇加军等,2010).孙茁等(2014)和Rajkumar等(2022)采用层状模型中sPn和Pn震相之间的走时差与震源深度的关系计算震源深度.

图1 区域深度震相示意图黑色、红色、绿色、黄色和蓝色曲线分别表示Pg、sPg、sPL、sPmP和sPn震相的射线路径.红色五角星表示震源,黑色三角形表示台站.

前人常在层状模型中利用深度震相测定震源深度,这种做法能简化深度震相走时计算,但是忽略了模型横向非均匀的影响(Yuan et al.,2020).当层状模型与真实地下结构之间存在较大差异时,可能会导致测定震源位置与实际位置之间出现较大偏差(Bock et al.,1996;Ma and Eaton,2011;Yuan et al.,2020;Diehl et al.,2021).随着地球物理观测数据的累积以及成像技术的发展,地球物理成像得到的三维非均匀模型越来越接近于实际地下结构.根据Diehl等(2021)的研究,区域地下介质模型不准能导致4 km以上的地震定位误差,而使用三维非均匀模型能减小模型误差对地震精定位的影响,因此有必要在三维非均匀结构模型之中开展地震定位研究.

为了实现三维非均匀介质模型中的地震精定位,本文通过快速行进法(Fast Marching Method,FMM)求解程函方程,得到区域非均匀模型中的初至P和sPg震相走时.联合P和sPg震相走时实现非均匀模型中的地震重定位.以2022年四川泸定MS6.8地震为例,检验本文地震定位方法对于提高地震定位精度的有效性.

1 P和sPg波走时联合地震定位方法

本文利用初至P波震相与深度震相sPg(图1)到时联合测定震源精确位置,为表述方便,后文将此地震定位方法称为联合定位方法.联合定位方法可归纳为两步.步骤1:利用实际的初至P波到时,对给定的初始震源位置进行重定位(Geiger,1912);步骤2:以步骤1的定位结果为初始位置,利用深度震相到时进一步提高地震定位精度.

关于初至P波地震重定位,其观测方程可表示为:

(1)

其中i表示地震编号,j表示台站编号,上标obs表示观测走时,cal表示计算走时,Δφi、Δλi和Δhi分别为地震i的初始震源位置在纬度、经度、震源深度方向的扰动,ΔTi为发震时刻扰动.将公式(1)中的走时差表示为向量形式b,震源位置及发震时刻的扰动由向量X表示,并将公式(1)中震源参数扰动系数置于灵敏度矩阵A中对应位置(Liu et al.,2019).然后,使用共轭梯度法迭代求解如下带阻尼的最小二乘问题(Puspito et al.,1993;Tong et al.,2014),实现基于初至P波到时的地震重定位:

(2)

公式(2)中Cd为先验数据协方差矩阵,Cm为先验模型协方差矩阵,η为阻尼因子.

在实现初至P波到时重定位后,本文利用sPg波到时进一步提高地震定位精度.利用sPg震相到时提高地震定位精度可分为两步.第一,在初至P波到时定位结果的周围一定范围内设置尝试网格点,网格点范围和网格点间距根据初至P波定位结果确定.网格点范围为水平方向±2.0 km、震源深度方向±10 km,对于地震定位,该范围可满足要求(罗艳等,2010).第二,计算这些网格点到台站的理论sPg震相走时,对比理论与实际sPg波走时并找出两者差别最小的网格点,此点即为最终的震源位置.其中理论sPg波走时的计算可分为三步.第一,将震中与台站之间的地表划分为等间距网格点.第二,分别计算震源到网格点的S波走时和台站到网格点的P波走时.第三,将所有的S波走时与P波走时相加,找出相加后的最小走时,此即为sPg波走时.

求解公式(2)和计算深度震相sPg走时都需要计算地震射线,本研究使用快速行进法(FMM)求解程函方程得到走时场和射线路径.FMM结合迎风差分格式和快速排序技术,可计算非均匀性较强介质(如:速度异常体尺度大于2倍波长尺度、相对扰动大于15%;Koulakov et al.,2016)中的地震波走时场和射线路径(Sethian and Povovici,1999;Rawlinson et al.,2006).非均匀介质中的程函方程可表示为:

(3)

=S(i,j,k),

(4)

(5)

(6)

(7)

(8)

(9)

(10)

其中δr、δθ和δφ分别表示r、θ和φ方向的离散步长,ri表示网格点(i,j,k)到球心的距离.

2 数据

2022年9月5日12时52分,在鲜水河断裂(磨西段)附近的四川省甘孜州泸定县(29.59°N,102.08°E)发生了MS6.8地震,震源深度16.0 km(由中国地震台网中心测定给出).在该地震的震中距约250 km范围内分布着24个固定台站,固定台站方位覆盖良好,平均台间距约为50 km(地震科学数据中心,2007;郑秀芬等,2009;图2).这些台站为本文的地震定位研究提供了宝贵数据.利用这些数据,本文开展相关合成测试和实际数据定位研究,验证本文提出的联合定位方法对地震定位精度提升的有效性.

图2 2022年泸定MS6.8地震震中及周围台站分布图红色五角星和圆表示地震,深红色实线为研究区域内主要断裂(Deng et al.,2003;An et al.,2023).倒三角表示台站,红色倒三角形为后续研究中用于计算sPg震相走时的台站.其中BHB:巴颜喀拉块体;SCB:四川盆地;WSB:川西块体;ANHF:安宁河断裂;LDF:理塘—德巫断裂;LMF:龙门山断裂;XSHF:鲜水河断裂;KD:康定地震;WC:汶川地震;LS:芦山地震;LD:泸定地震.图中的历史地震及泸定地震的震源机制解信息来源于中国地震台网中心(http:∥data.earthquake.cn).左下角插图显示研究区域的位置.

3 合成数据测试

Liu等(2021)联合大量体波和面波走时数据,反演了青藏高原东南缘地壳和最上层地幔速度结构,并给出三维区域非均匀模型SWChinaCVM-2.0.在地壳和地幔中该模型水平方向的分辨率分别约为20 km和40 km.后文中,我们采用三维模型SWChinaCVM-2.0开展合成测试研究,以验证本文基于程函方程的联合定位方法对于提升震源位置测定精度的有效性.合成测试研究主要讨论两种情况下的地震定位,第一,地震重定位时不存在模型误差的影响;第二,重定位时存在模型误差的影响.

首先展示重定位时不存在模型误差影响的情况.此情况下的合成测试分为两部分,初至P波到时重定位和sPg震相到时提高地震定位精度.初至P波到时重定位可归纳为如下三步:第一,在SWChinaCVM-2.0模型中采用FMM计算泸定地震震源位置(29.59°N,102.08°E,16 km)到台站(图2)的合成初至P波走时,并将标准差为 0.1 s 的高斯噪声加入至合成走时数据中,用于模拟含噪声的实际地震走时数据.第二,在泸定地震震源位置的经度、纬度和深度方向上分别加入标准差为7 km的高斯噪声获得一系列用于重定位的初始震源位置(表1).第三,使用合成走时数据和初至P波到时定位方法对表1所示的震源进行重定位,重定位时使用的模型为SWChinaCVM-2.0,重定位结果如表2所示.由表2可知,当重定位所使用的速度模型为真实模型时,初至P波到时定位方法较好地恢复了震源位置.重定位的结果与实际震源位置之间水平方向的误差平均值为0.1 km,深度方向的误差平均值为0.36 km.

表1 联合地震重定位方法所使用的初始震源位置

表2 初至P波到时重定位结果

使用sPg震相进一步约束震源位置.使用sPg震相到时约束震源位置可归纳为如下两步:第一,在SWChinaCVM-2.0中,使用FMM计算泸定地震震源到台站(图2)的合成sPg震相走时,将标准差为0.1 s的高斯噪声加入至合成走时数据中模拟实际sPg震相走时所含有的噪声.第二,将第一步中计算的合成走时当作实际sPg震相走时,然后在SWChinaCVM-2.0模型利用第1节中所述的sPg震相到时重定位方法对初至P波到时定位结果(表2)进行再次定位,重定位结果如表3所示.值得指出的是,将地表离散为等间距网格点,用于计算sPg震相走时,考虑计算时间,将网格间距设置为1 km.事实上,进一步加密网格,减小网格间距,虽然能提升定位精度,但实际测试表明进一步减小网格间距(例如,将网格间距降低至0.4 km)对定位精度提升有限.对比表3和表2可知,加入sPg震相到时约束震源位置的联合定位方法提高了震源深度的测定精度,深度方向定位误差的平均值从0.36 km减小到0.02 km.在水平方向上,联合定位方法相对于初至P波定位方法定位精度并未提高,这可能是因为在不存在模型误差影响情况下,初至P波定位方法已经获得了较准确的震中位置.

表3 联合方法的重定位结果

为了测试重定位使用的地下结构模型存在误差时,联合定位方法对震源精确位置的约束能力,本文基于AK135(Kennett et al.,1995)构建合成地壳速度梯度模型(崇加军等,2010),即速度随深度线性增加,为了论述方便后文中将此速度梯度模型记为M0.模型M0从地上5 km到地下35 km范围内,P波速度从5.8 km·s-1均匀增大到6.5 km·s-1,S波速度从3.46 km·s-1增大到3.85 km·s-1.合成测试时,首先在SWChinaCVM-2.0模型中计算初至P波走时和sPg波走时,将标准差为0.1 s的高斯噪声加入计算的走时数据中.然后将计算的走时数据当作实际观测数据,在M0模型中利用初至P波定位方法和本文联合定位方法开展地震重定位.初至P波重定位结果展示在表4中,联合方法定位结果展示在表5中.根据表4可知,当重定位使用的地下结构模型存在误差时(模型相对误差平均值为9.77%,标准差为10.22%;图3),使用初至P波到时定位方法仍能较好的约束震中位置,但得到的震源深度与真实情况偏离较大.重定位后得到的震源位置与实际位置在震中方向误差平均值约为2.3 km,深度方向的误差平均值约为19 km.联合定位方法改善了此情况下初至P波到时定位方法对震源深度约束不足问题,将定位误差平均值从19 km减小至3.19 km.此外,联合定位方法对震中位置的定位精度也有微弱的提高,震中位置的误差平均值从2.38 km减小到2.10 km.

表4 存在模型误差影响情况下初至P波到时重定位结果

表5 存在模型误差影响情况下使用联合方法的重定位结果

图3 模型相对误差图该图由M0模型与SWChinaCVM-2.0模型之差除以SWChinaCVM-2.0模型获得.红色数字为模型误差的平均值和误差的标准差.

为了进一步定量研究模型误差对本文联合地震定位方法的定位结果的影响,对SWChinaCVM-2.0模型分别施加5%、10%、15%和20%的随机误差,即对模型离散网格点上的速度值附加高斯噪声,获得误差模型M1、M2、M3和M4,并基于模型M1至M4开展重定位研究.重定位时,首先在SWChinaCVM-2.0模型中计算P与sPg波走时,并将此走时作为实际观测走时,然后分别在M1、M2、M3和M4中利用联合定位方法开展地震重定位,重定位结果相对于实际震源位置的误差展示在图4中.根据图4可知,地震定位的误差随着模型误差的增大而增大.震源深度方向的定位误差与模型误差成近似线性关系,模型误差每增加5%,定位误差的平均值约增加1.0 km.在震中位置方向定位误差与模型误差之间的线性关系不明显.

图4 地震定位误差图红色和黑色直方图分别表示震中方向和震源深度方向的定位误差平均值.

4 实际数据定位

本节选择2022年9月5日四川泸定MS6.8地震为研究对象,检验联合地震定位方法在实际地震定位中的有效性.首先根据图2所示台站记录的地震波形数据手动拾取初至P波(图5;表6)与sPg波(图6)到时,然后在区域非均匀模型SWChinaCVM-2.0之上开展联合方法地震重定位,初始震源位置为(102.08°E,29.59°N,16 km),重定位结果如表7所示.根据表7可知,初至P波到时重定位结果与联合方法重定位结果震中位置相差约为3.0 km,震源深度相差约为2.5 km. 此外,为了检验联合定位方法相比于前人深度震相定位方法的准确性,本文还使用Ma和Atkinson(2006)提出的深度震相模拟方法对泸定地震开展地震重定位.开展深度震相模拟方法重定位时,使用的层状模型如表8所示(Wang et al.,2007;王椿镛等,2010;孙茁等,2014),地震发震断层走向为163°、倾角为77°、滑动角为-5°(张喆等,2023).采用频率-波数(FK)方法计算震源深度6 km至11 km的理论地震图(Zhu and Rivera,2002),如图7所示.根据图7可知,随着震源深度增加,sPg与P波的到时差逐渐增大,实际地震图的sPg与P波到时差位于7 km和8 km理论地震图的到时差之间,因此,由深度震相模拟方法得出的泸定地震震源深度为7.5 km.

表6 根据图2所示台站记录的波形数据拾取的初至P波到时

表7 2022年泸定MS6.8地震重定位结果

表8 计算理论地震图使用的层状模型

图5 初至P波到时拾取图图中的黑色竖线标记初至P波到时,台站名显示在图的右上角,波形记录为垂直分量.左下角的数字表示纵坐标的倍数.

图6 sPg震相到时拾取图图中红色竖线标记初至P和sPg震相到时,台站名标记在图的右上角,震中距标记在图的右下角.在每个子图中,波形从上到下分别为径向、横向和垂向分量.

图7 台站SC.GZA处的合成与实际波形图图中红色曲线为不同震源深度的合成地震波形,黑色曲线为实际地震波形,所有波形记录根据初至P波到时对齐,波形记录均为垂直分量.实际地震波形记录与7 km和8 km的合成波形记录拟合较好.

为了检验本文联合方法定位结果的可信度,本文挑选了32个方位覆盖较好且独立的台站开展数值测试,台站分布如图8所示.根据图8所示的32个台站记录的地震波形数据,我们拾取了实际初至P波到时(表9).然后使用FMM分别计算联合定位方法给出的震源位置、初始震源位置、初至P波定位给出的震源位置和深度震相模拟方法给出的震源位置(102.08°E,29.59°N,7.5 km)到台站的合成走时,并计算实际走时与合成走时之间的走时差,走时差分布如图9所示.由图9可知,初始震源位置的走时差平均值和方差分别为0.22 s和0.82 s,经过联合方法重定位后走时差的平均值和方差分别减小为-0.07 s和0.67 s.从均值和方差的对比可知,联合方法定位结果比初始震源位置更加接近实际震源位置.此外,深度震相模拟方法定位结果的走时差均值和方差分别为-0.06 s和0.83 s,误差大于联合方法,联合定位方法给出的震源位置比深度震相模拟方法给出的震源位置更接近实际震源位置.

表9 根据图8所示台站记录的波形数据拾取的初至P波到时

图9 走时差分布直方图(a)、(b)、(c)和(d)由实际走时减去合成走时获得,但合成走时分别根据联合方法给出的震源位置、重定位初始震源位置、初至P波定位给出的震源位置和区域深度震相模拟给出的震源位置计算得到.右上角数字为走时差的平均值和方差.

5 结论与讨论

为了提高地震定位的精度,本文构造了基于程函方程的P和sPg波走时联合地震定位方法.在联合地震定位方法中,首先利用初至P波震相到时数据对给定初始震源位置进行重定位,然后再根据区域深度震相sPg到时进一步提高地震定位的精度.由于sPg震相是震源深度敏感震相,本文联合定位方法能较准确的约束震源深度.在联合定位过程中,本文采用快速行进法求解程函方程,得到非均匀模型中的P和sPg波走时,由此可以实现任意非均匀模型中的高精度地震定位.当非均匀模型接近实际地下结构时,联合定位方法能消除模型误差对定位结果的影响.通过合成测试分析,证实了基于程函方程的联合地震定位方法对提高地震定位的精度的有效性.利用本文的联合地震定位方法,对2022年四川泸定MS6.8地震开展实际数据的重定位研究,重定位结果显示泸定地震的震源位置为(102.094°E,29.601°N,13.87 km).根据走时差数据的分析,验证了重定位后的震源位置比初始震源位置(102.08°E,29.59°N,16.0 km)更接近实际震源位置.

虽然本文联合地震定位方法能有效提升地震定位的精度,但目前仍存在两方面的不足.第一,联合定位方法计算量较大,需要消耗较多计算时间.相比于初至P波走时定位,联合定位方法需要更多次的正演模拟,原因在于计算震源周围尝试网格点到台站的sPg波震相走时过程中,对于每个尝试网格点都需要一次正演模拟.对于本文区域地震定位而言,一次正演模拟的计算时间约为4 min,尝试网格点的个数约为1000,本研究中使用具有40个CPU计算核心的台式计算机,联合地震定位耗时约2 h.该计算时间远大于初至P波走时定位时间(约10 min).因此联合定位计算时间上并不占优势.但值得指出的是,对于每个尝试网格点走时场的数值求解具有天然并行性,可以通过增加CPU计算核心数目或者通过GPU(Graphics Processing Unit, 图像处理器)加速将有效减少计算时间,例如,在采用1000个计算核心计算时,联合定位计算时间将降至10 min以内.第二,联合定位方法难以将初至P和sPg波到时同时反演,而是通过两步反演的方式实现联合反演.在一阶线性近似条件下,初至P波走时扰动和震源参数扰动具有公式(1)所示的线性关系.对于sPg波,震源扰动会引起地表反射点变化,导致难以给出sPg波走时扰动与震源参数扰动之间的线性关系.因此,难以建立包含P和sPg波走时统一的反演观测方程用于震源参数同时反演.如何建立sPg波走时扰动与震源参数扰动之间的线性关系值得进一步研究.

本文所有数据和联合定位源代码可从https:∥doi.org/10.6084/m9.figshare.22336492.v1自由下载,定位代码可无偿用于学术用途使用.

致谢感谢审稿专家提出的建设性意见,这些意见对本文质量提升至关重要.中国地震局地球物理研究所“国家数字测震台网数据备份中心”和中国地震局地球物理研究所地震科学数据中心为本研究提供地震波形数据．本文所有图件由GMT(Generic Mapping Tools;Wessel et al.,2013)绘制.