2018, Vol. 61
地震定位是地震学最基本的问题之一.准确的地震定位对地震活动性分析、核爆监测和地球深部结构研究等都有重要意义(曾融生, 1991；Bondár and North, 1999；陈运泰等，2009).震源参数(震源位置和发震时间)通常根据观测到时确定，确定的震源参数使理论计算的和观测的到时相差最小(Thurber, 1986; Moser et al., 1992).因而，计算理论到时的速度模型对定位结果有直接的影响(Billings et al., 1994).如果速度模型远离实际，即使对地震进行相对定位也可能会产生较大的误差(Michelini and Lomax, 2004).大量研究结果表明地球内部、特别是在地壳和上地幔广泛存在较强的纵向和横向非均匀性(Goff and Holliger, 2003；滕吉文等，2004；Zhang et al., 2011)，但实际中的地震定位仍较普遍基于一维速度模型(Bondár et al., 2015).这在很大程度上与可用的三维速度模型较为稀缺有关，不过这种状况正不断得到改善.以中国为例，21世纪以来，使用分布密集的固定和/或流动地震台站观测记录，获得了华北(齐诚等，2006；Ding et al., 2009；于湘伟等，2010)、山东(苏道磊等, 2016)及川滇(见王椿镛等，2015)等很多地区的地壳高分辨率三维速度图像.三维精细速度模型的建立，使高精度的地震定位成为可能(陈棋福等，2001；Lomax et al., 2001; Chen et al., 2006；DeShon et al., 2007; de Kool and Kennett, 2014；Myers et al., 2015).
地震定位的精度除了受速度模型影响外，还与所拾取到时(或波形)的准确性、地震台网分布、所用震相种类及定位方法等因素有关(蔡明军等, 2004；杨文东等，2005; Bai et al., 2006).为更好地确定震源参数，地震学家提出了多种定位方法(田玥和陈晓非, 2002；蔡明军等, 2004；Richards et al., 2006)：盖革法(Geiger, 1912)、双差法(Waldhauser and Ellsworth, 2000)、干涉走时法(王璐琛等，2016)、网格搜索法(Shearer, 1997)、单纯形法(Prugger and Gendzwill, 1988)、模拟退火法(Billings，1994)、遗传算法(Kennett and Sambridge, 1992)、邻域算法(Sambridge and Kennett, 2001)、交切法(Pujol, 2004; 廉超等, 2006；Zhao and Ding, 2007，2009)及逆时成像法(许力生等, 2013a, 2013b)等.对于大多数方法而言，三维模型定位的主要困难在于走时或波形的正演模拟，但对于交切法则不只如此.交切法使用震源轨迹确定震中或震源位置，具有稳健直观的优点，在地震台网中有广泛的应用(孟玉梅等, 2001; 宋锐等, 2001).传统交切法通常基于均匀介质模型，震源轨迹以到时差约束，形状为球面或双曲面.当速度模型复杂时，震源轨迹则难以表达成解析形式.因此，将传统交切法扩展至三维速度模型还需要解决震源轨迹计算问题.对于该问题，目前主要有两种解决方案.一是以绝对残差较小的点表示震源轨迹(Zhou，1994；Font et al., 2004；Theunissen et al., 2012).这种方法简单易行，但得到的震源轨迹较为粗略，即使速度模型和到时都准确，震源轨迹通常也不会交于一点，而是交会成一个区域.第二种方案基于射线追踪技术：震源轨迹以地震波在绝对残差场中从最小或局部最小点至震源轨迹代表点的射线路径进行表示(赵爱华等, 2008, 2015).这种方法可精细计算震源轨迹，但还仅限于二维模型，同时兼顾计算轨迹的精细性与完整性较为困难.为精细且完整地计算震源轨迹，提出了一种类似削去层层外皮获得果核的算法从绝对残差较小的点中选取震源轨迹所在模型单元的节点(为便于叙述，将其称为震源轨迹节点)作为轨迹代表点.由于选取的轨迹代表点仅是部分轨迹节点，因此这种算法仅适用于二维模型；对于三维模型，计算的震源轨迹存在较多“空洞”，完整性较差.此外，目前计算震源轨迹的算法，还不太适应大批量数据的自动处理.例如，对于由多段组成的或稳定性较弱的震源轨迹，为防止产生虚假轨迹或保证轨迹的完整性，尚需人工精心设置轨迹代表点的选取范围和射线路径的计算区域.现代地震台网记录的地震事件及观测数据量都是巨大的，大量事件的高精度定位，要求震源轨迹的计算不仅准确而且效率高.
本文尝试在已有工作的基础上(赵爱华等, 2008, 2015)，发展适用于三维模型的震源轨迹计算方法.所期求的方法，除了计算的震源轨迹既精细且较完整外，还具有良好的稳健性、较强的自适应性和较高的计算效率，从而适用于大批量数据的自动处理.
1 震源轨迹方程设研究区地下的三维速度分布已知，在地表有M个地震台站Ri(i=1, 2, 3, …, M).将模型离散成大小相等的立方体单元，在单元内速度为常数，单元的中心点为模型节点.利用射线追踪方法，例如最小走时树方法(Klimeš and Kvasni ka, 1994；王辉和常旭，2000；Bai et al., 2007)、走时插值方法(Zhang et al., 2011, 2015)或程函方程有限差分方法(Vidale, 1990; Sethian and Popovic, 1999; de Kool et al., 2006)等，分别计算在台站Ri处激发的纵波(P波)和横波(S波)到点x =(x, y, z)的初至波走时T(x; Ri, P)和T(x; Ri, S).对于区内震源在x0=(x0, y0, z0)、发震时间为t0的地震，台站Ri记录到的P波、S波的初至到时分别记为τRiP和τRiS.根据震源-台站的对称性，可构建震源点所满足的空间方程.方程所表示的曲线(二维)或曲面(三维)，通常被称为震源的轨迹(参见傅淑芳和刘宝诚，1991)；相应的方程，本文称之为震源轨迹方程.以到时差和到时为约束条件的震源轨迹方程，构建如下(赵爱华等，2015).
1.1 到时差约束的震源轨迹每两个观测到时可组成1个到时对.设组成到时对j的到时分别为τRj1Wj1和τRj2Wj2，其中Rj1和Rj2为记录到时的台站，Wj1和Wj2为震相类型，即P或S.到时τRj1Wj1与τRj2Wj2之差约束的震源轨迹满足如下方程：
|
(1) |
在方程(1)中，左边为理论的到时差，右边为观测的到时差.到时差残差(RDT)场FRDT为：
|
(2) |
其中，I(x)为单位场，即其值处处等于1.
1.2 到时约束的震源轨迹根据观测到时可构建事件的发震时间场：
|
(3) |
其中，Gj和Wj分别为第j个观测到时τGjWj的台站和震相类型(P或S波)，K为观测到时个数.这样，到时τGjWj约束的震源轨迹满足如下方程：
|
(4) |
到时残差(RAT)场FRAT为：
|
(5) |
其中，I(x)为单位场.
2 轨迹计算对于复杂速度模型，震源轨迹方程(1)或(4)难以给出解析解.但从式(2)或(5)可知，震源轨迹在绝对残差场|FRDT|或|FRAT|中具有接近于0的低值.绝对残差场中数值较小的区域依震源轨迹的几何分布自适应地划分成1个或多个连通区域.若在每个连通区域内最小点处激发地震波，地震波将沿震源轨迹快速传播，而在震源轨迹两侧传播较慢.因此，若选取部分模型节点作为震源轨迹的代表点，则代表点至与其同一连通区域内地震波初始点的射线路径可近似地表示震源轨迹.为使计算的震源轨迹精细，震源轨迹代表点应是震源轨迹所在模型单元的节点(简称震源轨迹节点).在三维模型中，震源轨迹为曲面.在曲面上连接两点的路径不止一条.根据费玛原理可知，地震波沿走时最小的路径传播，会尽可能绕开绝对残差较大的轨迹节点.若这些轨迹节点没有被选为轨迹代表点，最后计算的震源轨迹就有可能因没有与其相连的射线路径而存在“空洞”，如图 1所示.因此，为精细且完整地计算震源轨迹，选取的震源轨迹代表点应不仅是而且包括全部的轨迹节点.基于这种认识，提出了适用于三维模型的震源轨迹计算方法.
|
图 1 计算的震源轨迹(蓝色线)存在空洞(红色椭圆)情形的示意图 棕色点为选取的震源轨迹参考点; 绿色十字为实际震源轨迹所经过模型单元的节点. Fig. 1 Diagram of a calculated hypocentral locus (blue lines) with holes (red ellipses) Brown dots are the selected focal locus reference points, and the green crosses indicate the nodes of the model elements passed by the real hypocentral locus. |
震源轨迹代表点通常在绝对残差场中选取，选取的依据主要是(1)轨迹节点具有很小的值和(2)在轨迹的法线方向上轨迹节点的数值局部最小.轨迹节点至实际轨迹的距离有远有近，其绝对残差有大有小.当震源轨迹的稳定性较弱时，绝对残差很小的点离震源轨迹也可能较远.单单根据绝对残差大小，仅能大致确定震源轨迹代表点的选取范围.在绝对残差场中，确定震源轨迹的法线方向较为困难，这是因为轨迹法线方向并不总是最大梯度方向，而存在局部最小点的最大梯度方向也不全是轨迹法线方向.选取局部最小点或最大梯度方向上局部最小点作为轨迹代表点，都难以使选取的代表点仅是或包含全部震源轨迹节点.
和绝对残差场相比，残差场包含更丰富的信息.图 2显示了均匀介质模型中虚拟事件的震源轨迹及其残差场.震源轨迹来自文献(赵爱华等, 2015, 图 1)，模型以3 km×3 km的网格剖分.从图 2可以看出，无论是仅有1段的震源轨迹(a)还是有多段的震源轨迹(b)，都有如下特征：(1)震源轨迹(图中紫色实线)位于残差正(红色)负(蓝色)节点之间；(2)在震源轨迹法线方向上残差变化最快，即绝对梯度最大；(3)在震源轨迹法线方向上，离震源轨迹越近残差的绝对值越小.根据以上三个特征，本文提出了一种在残差场中选取震源轨迹代表点的方法.
|
图 2 二维均匀介质中事件的震源轨迹及其残差场 (a)和(b)中的震源轨迹对应同一P波到时，但其发震时刻场分别以P到时单独、P和S波到时联合构建.金色十字和棕色点分别为“削皮”算法和本文方法选取的震源轨迹代表点.紫色和绿色线分别为震源轨迹的理论解和计算结果. Fig. 2 Hypocentral loci in their residual fields of an event in a two-dimensional homogeneous medium Hypocentral loci respond to the same P-wave arrival time but their origin time fields are constructed respectively with (a) only P-wave arrival times and (b) both P- and S-wave data. The golden crosses and the brown dots indicate the hypocentral locus reference points selected with the peeling method and the method in this study, respectively. The green and the purple lines respectively show the calculated and theoretical hypocentral loci. |
选取轨迹代表点的步骤如下：
(1) 组成点对
将每一个模型节点和残差正负极性与其不同的邻点组成点对.对于二维模型，考虑节点上下、左右4个方向的邻点；对于三维模型，考虑节点上下、左右、前后6个方向的邻点.
(2) 确定法线点对
对于同一个模型节点，选取绝对梯度最大的点对为法线点对.
(3) 选取轨迹代表点
选取法线点对中残差绝对值较小的点为震源轨迹代表点.
在上述算法中，组成点对的邻点残差正负极性不同，即其残差相乘之积小于或等于零，表明点对中的两点分别位于震源轨迹的两侧.对于某个模型节点，若其位于震源轨迹另外一侧的邻点不止一个，则该点与其邻点可组成多个两点残差正负极性不同的点对.绝对梯度最大的点对(即法线点对)所在方向最接近轨迹法线方向，两者相差不超过45°.在轨迹法线方向上，离震源轨迹越近残差绝对值越小.由此推知，法线点对中残差绝对值较小的点更靠近震源轨迹.这样，选取的轨迹代表点就是在包含震源轨迹的1个模型单元范围内最接近震源轨迹的节点，即震源轨迹节点.图 2中的金色十字和棕色点分别为“削皮”算法(赵爱华等，2015)和本文方法选取的震源轨迹参考点.可以看出，本文方法选取的轨迹参考点不仅是而且包含全部震源轨迹节点，而“削皮”算法仅选取部分轨迹节点为轨迹参考点.值得指出的是，当震源轨迹在两个邻点间穿越2次或更多偶数次时，由于邻点的残差正负极性相同，本文方法会漏选其轨迹节点，但这种情况很少.
2.2 射线路径计算选取的震源轨迹代表点不仅是而且包括了全部或几乎全部的震源轨迹节点，可直接用于地震定位.若将这些离散的模型节点用射线路径连接起来则可更直观地展示震源轨迹的几何形态，这对评估定位效果、分析速度模型等因素对定位结果的影响是有利的.连接轨迹代表点的射线路径利用改进的最小走时树方法(Zhao et al, 2004; 赵爱华和张中杰，2004)计算.和其他射线追踪方法相比，最小走时树方法具有很好的稳健性，适用于复杂速度模型，并且可同时计算地震波从震源至其他所有模型节点的走时和射线路径(Nakanishi and Yamaguchi, 1986; Moser, 1991；Cao and Greenhalgh, 1993).为防止产生虚假震源轨迹或使轨迹的完整性受到损害，和二维模型的做法一样(赵爱华等，2015)，射线追踪区域也被限制在绝对残差场中数值小于RFL的区域，区域中的每个连通区域完整包含且仅包含1段震源轨迹；以每个连通区域内绝对残差最小的点为地震波初始点.注意到震源轨迹代表点包含了全部或几乎全部轨迹节点，这里将RFL取为比轨迹代表点绝对残差最大值RFLRPmax略大(例如大0.01)的值.对于二维和三维模型，计算震源轨迹射线路径的算法具有相同的框架，具体如下：
(1) 查找地震波初始点
在震源轨迹代表点集合Z中查询绝对残差最小的点s, 将s从Z中移出.
(2) 射线路径计算
以s为地震波初始点，在绝对残差场中计算地震波走时及射线路径直至最小走时树算法终止；若节点j∈Z被计算过，则将其从Z中移出.
(3) 迭代条件检测
如果Z成为空集则结束；否则回到(1).
轨迹代表点至其各自地震波初始点的射线路径有很多部分是重复的.为去掉冗余射线，在统计表示震源轨迹的射线路径时，仅考虑射线路径经过次数为1的震源轨迹代表点.对于这些代表点，其射线路径重复部分也仅统计1次.图 2中震源轨迹的计算结果如图中绿线所示.可以看出，计算的震源轨迹和理论结果(图中紫色实线)符合得很好.
3 算例 3.1 速度模型
地球内部的速度结构是复杂的，既存在较强的速度间断面，也存在速度在空间上的连续变化(例如, 白志明和王椿镛，2004).为此，选用了一个较复杂的三维速度模型来检验本文方法的实用性.速度模型来自文献(赵爱华和张中杰，2004)，如图 3所示.水平界面R1深度为0.3 km, 倾斜界面R2满足方程x+y+10z+6.0=0.在这两个界面之间嵌有一椭球体E，中心点坐标为(0.5 km, 0.5 km, -0.5 km)，其形状满足方程
|
|
图 3 三维复杂速度结构模型 三角形表示地震台站；红色符号“+”为实际震源位置. Fig. 3 A 3-D complex velocity model Triangles refer to the seismic stations and the red signs of "+" to the true hypocenter. |
界面R2以下介质Ⅲ为均匀介质，其地震波速度为vP=3.5 km·s-1, vS=2.1 km·s-1.椭球体内的介质速度为：vP=3.0 km·s-1, vS=1.7 km·s-1.
图 3模型用10 m×10 m×10 m的网格剖分.和周建超和赵爱华(2012)的做法一样，将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(315 m, 735 m, -655 m)，即图 2中红色“+”处，发震时间t0=0.0 ms.由震源x0至地震台站Ri的P、S波到时τRiP和τRiS利用三维的最小走时树射线追踪方法(赵爱华和张中杰，2004)计算，子波传播区域半径取1.
3.2 震源轨迹对于图 3中的事件，9个台站的P波到时两两组合，总计可产生C92=36个到时差.这些到时差约束的震源轨迹有很多在几何形态上非常接近.根据几何形态，本文选取了其中有代表性的3条显示在图 4a—4c中.轨迹旁边的数字对“i-j”表示该轨迹以到时差τRiP-τRjP约束.可以看到，计算的震源轨迹不仅很精细，为恰好经过震源的曲面而非包含震源的板带，并且相当完整，不存在“空洞”；受介质非均匀性的影响，轨迹的形状不再是双曲面；轨迹在震源附近以垂直向分布为主，对震中有较好的约束.
|
图 4 图 3中事件以到时差约束的震源轨迹 左图以不同台站P波到时差约束；右图以相同台站的P波和S波到时差约束. Fig. 4 Hypocentral loci constrained with arrival time differences for the event in Fig. 3 The left hypocentral loci are constrained with arrival time differences between P waves at different stations; the right those with ones between P and S waves at the same stations. |
同一台站P波和S波到时差约束的震源轨迹总计有9条.根据几何形态，从其中选取了有代表性的3条显示在图 4d—4f中.轨迹旁边的数字“i”表示该轨迹以到时差τRiS-τRiP约束.震源轨迹不像在均匀介质中那样为球面形.图 4d中的轨迹受介质非均匀性的影响最为显著，该轨迹存在明显的内凹形变.距离震中较近的台站，其P-S波到时差约束的震源轨迹(例如图 4f)在震源附近有较大的水平分量，可较好地约束震源深度.
到时约束的震源轨迹涉及发震时刻场构建.构建发震时刻场的到时个数越多，震源轨迹理论求解越困难；使用多个观测到时也有利于减小随机干扰的影响.为此，这里考虑全部观测台站的情况.全部台站P波到时约束的震源轨迹有9条，P波和S波到时联合约束的震源轨迹有18条.根据几何形态，从这两种轨迹中各选出有代表性的3条显示在图 5中.轨迹旁边的数字“i”表示该轨迹对应到时τRiP.可以看到，对于同一到时，发震时刻场不同，震源轨迹的几何形态也不相同，甚至相差较大，例如图 5a与图 5d.和到时差约束的震源轨迹相比，到时约束的震源轨迹在几何形态上更具有多样性.例如，图 5c中的震源轨迹总体外凸，局部内凹，既不是球面也不是双曲面；图 5f中的震源轨迹由分开的上下两部分组成，上部为桶状，下部为向外、向上伸展的不规则曲面.到时约束的震源轨迹几何形态复杂，意味着其在空间上可能交会于多个点，而每个交会点均是总绝对到时残差函数的极小值点.因此，当以盖革法(以总的到时残差为目标函数)定位时，所给初始震源参数应离真值较近.
|
图 5 图 3中事件以到时约束的震源轨迹 左图和右图震源轨迹都对应P波到时，但其发震时刻场分别以P波到时单独、P和S波到时联合构建. Fig. 5 Hypocentral loci constrained with arrival times for the event in Fig. 3 The left and right hypocentral loci respond to P waves but their origin time fields are constructed respectively with only P-wave arrival times and both P- and S-wave data. |
在三维模型中震源轨迹通常为连续的曲面.图 5f中的轨迹在其底部存在孤立的突起(红色箭头所指)，显得异常.为此，以图 5c和图 5f两个几何形态较为复杂的震源轨迹作为代表，对计算结果进行了检验.
图 3模型速度结构较为复杂，震源轨迹解析求解较为困难.由于震源轨迹计算的正确与否关键取决于轨迹代表点选取的是否恰当，因此计算结果可通过考查所选震源轨迹代表点在残差场中的分布进行检验.为此，对图 5c和图 5f中轨迹的残差场分别于震源点(315 m, 735 m, -655 m)和近模型中心点(505 m, 505 m, -505 m)处做垂直于坐标轴的3个切片.残差场切片如图 6所示，棕色点为选取的震源轨迹代表点.可以看出，所选轨迹代表点恰好沿震源轨迹(白色部分)分布，而且没有间断.由此推断，基于这些代表点计算的震源轨迹应和实际相符.
|
图 6 图 5c(左图)和5f(右图)中震源轨迹所在到时残差场的切片 棕色点为选取的震源轨迹代表点；红色“+”为震源在切片上的投影. Fig. 6 Slices of arrival time residual fields containing the hypocentral loci in Figs. 5c (left) and 5f (right) The brown dots are the selected hypocentral locus reference points. The signs of red "+" indicate the projective points of the true hypocenter onto the slices. |
在图 6d和6e中，可以看到在震源轨迹左下方有一小段残差为正、但数值很小的区域(黑色箭头所指).结合图 5f中突起的孤立震源轨迹来看，该低值部分可能也包含震源轨迹.若如此，该段轨迹被漏掉的原因估计是轨迹在两个相邻节点之间穿越了两次或更多偶数次致使轨迹两侧节点的残差正负极性相同.和整个震源轨迹相比，该部分所占比例很小，并且离震源较远；更重要的是，这种情况很少.由于实际定位通常使用多条震源轨迹，因此该部分即使真是震源轨迹而没被计算上也不会对定位结果有较大影响.
需要说明的是，实际中的速度模型和拾取的到时，受多种因素影响，总会存在或大或小的误差，从而计算的震源轨迹通常不会恰好经过实际震源.这样，同一事件的震源轨迹将不会交会于一点，而是交会成一个区域.将交会区内的哪个点(例如轨迹经过最多的点、交会区的中心点或加权中心点)确定为震源更好，这还是一个需要研究的问题.
4 讨论 4.1 算法对弱稳定性震源轨迹的适用性震源轨迹对震源位置的约束能力，不仅与轨迹在震源附近的方位分布有关，而且还与轨迹的稳定性有关(赵爱华等，2015).震源轨迹越稳定，受扰动影响越小，对震源的约束越强；反之亦然.将绝对残差场看作地形图，弱稳定性震源轨迹所在的“峡谷”较为宽缓.通常来讲，两个相距较近的远台P波到时差约束的震源轨迹和远台P波到时约束的震源轨迹，稳定性较弱；但包括近台和远台的以到时约束的震源轨迹也有稳定性弱的.定位实践中，通常使用所有可能的震源轨迹确定震源位置(Font et al., 2004; Theunissen et al., 2012)，而这就可能涉及弱稳定性震源轨迹的计算.
图 7显示了一条具有弱稳定性的震源轨迹及其绝对残差场.震源轨迹来自文献(赵爱华等, 2015, 图 8a)：速度模型以云南地区为背景，对遮放—宾川速度剖面(白志明和王椿镛，2004)略加平滑而成，以1 km×1 km的网格剖分；虚拟震源位于点(130.5 km, -27.5 km)处，如红色“+”所示.为清晰显示震源轨迹，对绝对残差做了取对数处理.易见，震源轨迹由分开的左右两段组成，右半段稳定性较弱.用本文方法对该条轨迹计算.可以看到，选取的轨迹代表点(棕色点)恰好沿震源轨迹曲线分布，即使对于稳定性较弱的右半段轨迹亦是如此.在轨迹形状急剧变化的部位(红色圆圈所示)，轨迹代表点的分布似乎有所间断，其原因估计是(1)轨迹穿越同一模型单元两次，而轨迹节点仅有1个；或者(2)轨迹近乎平行地穿越两个相邻模型单元，单元节点正好组成法线点对，尽管它们均是轨迹节点，但根据规则仅能选取其中之一作为轨迹代表点.由于轨迹形状变化剧烈的情况较少，轨迹节点漏选的范围较小，而震源参数的确定通常使用多条震源轨迹，因此轨迹代表点的这种间断对地震定位的结果影响不大.计算的震源轨迹总体上和实际符合得很好，仅在红色圆圈所示部分略有不符：射线路径不是沿轨迹曲线将轨迹代表点相连，而是以直线将其连接.这是由于震源轨迹在该部分稳定性较弱、形状变化剧烈，地震波沿着轨迹之间区域中的直线比沿着震源轨迹曲线传播得更快.该算例表明，即使对稳定性较弱的震源轨迹，本文方法也能较好地计算.
|
图 7 弱稳定性震源轨迹及其绝对残差场 棕色点为选取的震源轨迹代表点；蓝色实线为计算的震源轨迹. Fig. 7 A weakly stable hypocentral locus and its absolute residual field Brown dots indicate the selected FLRPs; blue solid lines are the calculated hypocentral locus. |
|
图 8 图 2a中震源轨迹代表点(实线)及其初选点(虚线)的百分比 Fig. 8 Percentages of primary (dashed line) and final (solid line) reference points for the hypocentral locus in Fig. 2a |
震源轨迹计算包括三个步骤：(1)地震波走时场计算、(2)震源轨迹代表点选取和(3)震源轨迹射线路径计算.其中，第一步最为费时，尤其是三维速度模型.这一步仅涉及速度模型和地震台站、与地震事件无关.地球内部结构的变化通常是很缓慢的，短时间内可看作是稳定的.这样，地震波走时场仅需计算1次.
震源轨迹代表点从整个模型空间选取，即使震源轨迹由分开的多段组成或稳定性较弱，选取的轨迹代表点不仅是而且包含了全部或几乎全部震源轨迹节点.这部分仅涉及简单的运算，计算速度很快，并且有较大的提升空间.注意到震源轨迹节点具有较小的绝对残差值，因此可从残差绝对值较小的点中选取震源轨迹代表点.为便于叙述，将残差绝对值小于正实数RFLRP的点称为初选轨迹代表点.图 8显示了均匀模型算例中图 2a震源轨迹其初选(虚线)、实际(实线)代表点占所有模型节点百分比随参数RFLRP的变化，RFLRP的步长取0.005 s.可以看出，自RFLRP大于0.15 s(确切地说是0.120 s, 轨迹节点的最大绝对残差为0.118 s)之后，震源轨迹代表点的选取结果都恒定不变.这表明，为保证计算轨迹的完整性，参数RFLRP只要大于震源轨迹节点的最大绝对残差值即可.以RFLRP取较大的0.3 s为例，初选震源轨迹代表点不到模型单元节点总数的25%，从这些初选点中选取轨迹代表点可大大节省计算量.
震源轨迹射线路径的计算限制在绝对残差小于正实数RFL的区域. RFL取为比轨迹代表点的最大绝对残差值RFLRPmax略大(例如大0.01)的值.这样，既可保证计算轨迹的完整性，也可防止虚假震源轨迹的产生，同时也使得计算区域较小、大大提高了计算效率. RFL依据选取的轨迹代表点设置，这样轨迹射线路径的计算就不需要人工干预，从而适于计算机的自动处理.
根据以上分析可知，对于速度结构稳定的地区，当处理同一台网记录的多个事件时，震源轨迹计算具有较高的效率；对于震源轨迹的组成段数，算法没有限制，即使对弱稳定性的震源轨迹也能很好地计算——这意味着，对于多个事件，在构建其震源轨迹时，可用相同方式对观测到时进行组合；计算参数在很宽的范围内取值，计算的震源轨迹都精细且较完整.综合这些特点来看，本文方法对大批量数据的自动处理有较好的适应性.
5 结论本文提出的震源轨迹计算方法基于最小走时树射线追踪技术，对三维复杂速度模型具有良好的稳健性.表征震源轨迹的代表点可在整个模型空间选取，对轨迹的组成段数没有限制，选取的轨迹代表点不仅是而且包含全部或几乎全部震源轨迹节点，即使震源轨迹稳定性较弱.射线追踪仅限于绝对残差很小的区域，区域根据选取的轨迹代表点自适应地设置，便于计算机的自动处理.最为费时的地震波走时场计算仅需1次，对于批量事件，震源轨迹计算效率较高.
数值模型计算表明：计算的震源轨迹不仅精细而且较完整，可用于高精度的地震定位(包括震源位置的确定、定位效果的评估及速度模型等因素的影响分析等).介质的非均匀性对震源轨迹的几何形态有显著影响.和到时差约束的震源轨迹相比，到时约束的震源轨迹在几何形态上更具有多样性.
致谢 中国地震局地球物理研究所吴庆举研究员和丁志峰研究员对本项工作给予了支持，作者同时也感谢评审专家的中肯评述.
Bai C, Greenhalgh S, Zhou B. 2007. 3D ray tracing using a modified shortest-path method. Geophysics, 72(4): T27-T36. DOI:10.1190/1.2732549 |
Bai L, Wu Z L, Zhang T Z, et al. 2006. The effect of distribution of stations upon location error:Statistical tests based on the double-difference earthquake location algorithm and the bootstrap method. Earth Planets Space, 58(2): e9-e12. DOI:10.1186/BF03353364 |
Bai Z M, Wang C Y. 2004. Tomography research of the Zhefang-Binchuan and Menglian-Malong wide-angle seismic profiles in Yunnan province. Chinese J. Geophys (in Chinese), 47(2): 257-267. |
Billings S D. 1994. Simulated annealing for earthquake location. Geophys. J. Int., 118(3): 680-692. DOI:10.1111/gji.1994.118.issue-3 |
Billings S D, Sambridge M S, Kennett B L N. 1994. Errors in hypocenter location:picking, model, and magnitude dependence. Bull. Seismol. Soc. Am., 84(6): 1978-1990. |
Bondár I, Myers S C, Engdahl E R. 2015. Earthquake location.//Beer M, Kougioumtzoglou I A, Patelli E, et al eds. Encyclopedia of Earthquake Engineering. Heidelberg: Springer, 661-676.
|
Bondár I, North R G. 1999. Development of calibration techniques for the Comprehensive Nuclear-Test-Ban Treaty (CTBT) international monitoring system. Phys. Earth Planet. Inter., 113(1-4): 11-24. DOI:10.1016/S0031-9201(99)00033-3 |
Cai M J, Shan X M, Xu Y, et al. 2004. Review of earthquake-locating methods from error. Journal of Seismological Research (in Chinese), 27(4): 314-317. |
Cao S H, Greenhalgh S. 1993. Calculation of the seismic first-break time field and its ray path distribution using a minimum traveltime tree algorithm. Geophys. J. Int., 114(3): 593-600. DOI:10.1111/gji.1993.114.issue-3 |
Chen H, Chiu J M, Pujol J, et al. 2006. 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., 96(1): 288-305. DOI:10.1785/0120040102 |
Chen Q F, Zhang Y Q, Zhou J. 2001. Review of earthquake location with three dimensional earth model for digital seismic observation. Earthquake (in Chinese), 21(2): 29-40. |
Chen Y T. 2009. Earthquake prediction:Retrospect and prospect. Sci. China Ser. D-Earth Sci (in Chinese), 39(12): 1633-1658. |
de Kool M, Kennett B L N. 2014. Practical earthquake location on a continental scale in Australia using the AuSREM 3D velocity model. Bull. Seismol. Soc. Am., 104(6): 2755-2767. DOI:10.1785/0120140105 |
de Kool M, Rawlinson N, Sambridge M. 2006. A practical grid-based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media. Geophys. J. Int., 167(1): 253-270. DOI:10.1111/gji.2006.167.issue-1 |
DeShon H R, Thurber C H, Rowe C. 2007. High-precision earthquake location and three-dimensional P wave velocity determination at Redoubt Volcano, Alaska. J. Geophys. Res., 112(B7): B07312. DOI:10.1029/2006JB004751 |
Ding Z F, Zhou X F, Wu Y, et al. 2009. Tomographic imaging of P wave velocity structure beneath the region around Beijing. Earthquake Science, 22(4): 403-408. DOI:10.1007/s11589-009-0403-9 |
Font Y, Kao H, Lallemand S, et al. 2004. Hypocentre determination offshore of eastern Taiwan using the Maximum intersection method. Geophys. J. Int., 158(2): 655-675. DOI:10.1111/gji.2004.158.issue-2 |
Fu S F, Liu B C. 1991. A Tutorial of Seismology (in Chinese). Beijing: Seismological Press.
|
Geiger L. 1912. Probability method for the determination of earthquake epicenters from the arrival time only. Bull. St. Louis Univ., 8: 60-71. |
Goff J A, Holliger K. 2003. Heterogeneity in the Crust and Upper Mantle: Nature, Scaling, and Seismic Properties. New York: Kluwer Academic/Plenum Publishers.
|
Klimeš L, Kvasnička M. 1994. 3-D network ray tracing. Geophys. J. Int., 116(3): 726-738. DOI:10.1111/gji.1994.116.issue-3 |
Kennett B L N, Sambridge M S. 1992. Earthquake location-genetic algorithms for teleseisms. Phys. Earth Planet. Inter., 75(1-3): 103-110. DOI:10.1016/0031-9201(92)90121-B |
Lian C, Li S L, Dong M, et al. 2006. Method of spherical surface shearing in earthquake location. Journal of Geodesy and Geodynamics (in Chinese), 26(2): 99-103. |
Lomax A, Zollo A, Capuano P, et al. 2001. Precise, absolute earthquake location under Somma-Vesuvius volcano using a new three-dimensional velocity model. Geophys. J. Int., 146(2): 313-331. DOI:10.1046/j.0956-540x.2001.01444.x |
Meng Y M, Zhao Y, Wang B. 2001. The analysis on deviation of rapid location by China Seismic Observational Network. Earthquake (in Chinese), 21(3): 65-69. |
Michelini A, Lomax A. 2004. The effect of velocity structure errors on double-difference earthquake location. Geophys. Res. Lett., 31(9): L09602. DOI:10.1029/2004GL019682 |
Moser T J. 1991. Shortest path calculation of seismic rays. Geophysics, 56(1): 59-67. DOI:10.1190/1.1442958 |
Moser T J, van Eck T, Nolet G. 1992. Hypocenter determination in strongly heterogeneous Earth models using the shortest path method. J. Geophys. Res., 97(B5): 6563-6572. DOI:10.1029/91JB03176 |
Myers S C, Simmons N A, Johannesson G, et al. 2015. Improved regional and teleseismic P-wave travel-time prediction and event location using a global 3D velocity model. Bull. Seismol. Soc. Am., 105(3): 1642-1660. DOI:10.1785/0120140272 |
Nakanishi I, Yamaguchi K. 1986. A numerical experiment on nonlinear image reconstruction from first-arrival times for two-dimensional island arc structure. J. Phys. Earth, 34(2): 195-201. DOI:10.4294/jpe1952.34.195 |
Prugger A F, Gendzwill D J. 1988. Microearthquake location:A nonlinear approach that makes use of a simplex stepping procedure. Bull. Seismol. Soc. Am., 78(2): 799-815. |
Pujol J. 2004. Earthquake location tutorial:graphical approach and approximate epicentral location techniques. Seismol. Res. Lett., 75(1): 63-74. DOI:10.1785/gssrl.75.1.63 |
Qi C, Zhao D P, Chen Y, et al. 2006. 3-D P and S wave velocity structures and their relationship to strong earthquakes in the Chinese capital region. Chinese J. Geophys (in Chinese), 49(3): 805-815. |
Richards P G, Waldhauser F, Schaff D, et al. 2006. The applicability of modern methods of earthquake location. Pure and Applied Geophysics, 163(2-3): 351-372. DOI:10.1007/s00024-005-0019-5 |
Sambridge M S, Kennett B L N. 2001. Seismic event location:nonlinear inversion using a neighbourhood algorithm. Pure and Applied Geophysics, 158(1-2): 241-257. |
Sethian J A, Popovic A M. 1999. 3-D traveltime computation using the fast marching method. Geophysics, 64(2): 516-523. DOI:10.1190/1.1444558 |
Shearer P M. 1997. Improving local earthquake location using the L1 norm and waveform cross correlation:Application to the Whittier Narrows, California, aftershock sequence. J. Geophys. Res., 102(B4): 8269-8283. DOI:10.1029/96JB03228 |
Song R, Gu X H, Wang Y L, et al. 2001. The software system for real-time processing and large earthquake rapid determination of NCDSN. Earthquake (in Chinese), 21(4): 47-59. |
Su D L, Fan J K, Wu S G, et al. 2016. 3D P wave velocity structures of crust and their relationship with earthquakes in the Shandong area. Chinese J. Geophys (in Chinese), 59(4): 1335-1349. DOI:10.6038/cjg20160415 |
Teng J W, Zhang Z J, Bai W M, et al. 2004. Physics of Lithosphere (in Chinese). Beijing: Science Press.
|
Theunissen T, Font Y, Lallemand S, et al. 2012. Improvements of the maximum intersection method for 3D absolute earthquake locations. Bull. Seismol. Soc. Am., 102(4): 1764-1785. DOI:10.1785/0120100311 |
Tian Y, Chen X F. 2002. Review of seismic location study. Progress in Geophys (in Chinese), 17(1): 147-155. |
Thurber C H. 1986. Analysis methods for kinematic data from local earthquakes. Rev. Geophys., 24(4): 793-805. DOI:10.1029/RG024i004p00793 |
Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions. Geophysics, 55(5): 521-526. DOI:10.1190/1.1442863 |
Waldhauser F, Ellsworth W L. 2000. A double-difference earthquake location algorithm:method and application to the northern Hayward fault, California. Bull. Seismol. Soc. Am., 90(6): 1353-1368. DOI:10.1785/0120000006 |
Wang C Y, Yang W C, Wu J P, et al. 2015. Study on the lithospheric structure and earthquakes in North-South Tectonic Belt. Chinese J. Geophys (in Chinese), 58(11): 3867-3901. DOI:10.6038/cjg20151101 |
Wang H, Chang X. 2000. 3-D ray tracing method based on graphic structure. Chinese J. Geophys (in Chinese), 43(4): 534-541. |
Wang L, Chang X, Wang Y B. 2016. Locating micro-seismic events based on interferometric traveltime inversion. Chinese J. Geophys (in Chinese), 59(8): 3037-3045. DOI:10.6038/cjg20160826 |
Xu L S, Du H L, Yan C, et al. 2013a. A method for determination of earthquake hypocentroid:time-reversal imaging technique I-Principle and numerical tests. Chinese J. Geophys (in Chinese), 56(4): 1190-1206. DOI:10.6038/cjg20130414 |
Xu L S, Yan C, Zhang X, et al. 2013b. A method for determination of earthquake hypocentroid:time-reversal imaging technique Ⅱ-An examination based on people-made earthquakes. Chinese J. Geophys (in Chinese), 56(12): 4009-4027. DOI:10.6038/cjg20131207 |
Yang W D, Jin X, Li S Y, et al. 2005. Study of seismic location methods. Earthquake Engineering and Engineering Vibration (in Chinese), 25(1): 14-20. |
Yu X W, Chen Y T, Zhang H. 2010. Three-dimensional crustal P-wave velocity structure and seismicity analysis in Beijing-Tianjin-Tangshan Region. Chinese J. Geophys (in Chinese), 53(8): 1817-1828. DOI:10.3969/j.issn.0001-5733.2010.08.007 |
Zeng R S. 1991. The problems of lithospheric structure and geodynamics. Advances in Earth Science (in Chinese), 6(2): 1-10. |
Zhang J Z, Huang Y Q, Song L P, et al. 2011a. Fast and accurate 3-D ray tracing using bilinear traveltime interpolation and the wave front group marching. Geophys. J. Int., 184(3): 1327-1340. DOI:10.1111/gji.2011.184.issue-3 |
Zhang J Z, Shi J J, Song L P, et al. 2015. Linear traveltime perturbation interpolation:a novel method to compute 3-D traveltimes. Geophys. J. Int., 203(1): 548-552. DOI:10.1093/gji/ggv316 |
Zhang Z J, Yang L Q, Teng J W, et al. 2011b. An overview of the earth crust under China. Earth-Science Reviews, 104(1-3): 143-166. DOI:10.1016/j.earscirev.2010.10.003 |
Zhao A H, Ding Z F. 2007. An intersection method for locating earthquakes in complex velocity models. Applied Geophysics, 4(4): 294-300. DOI:10.1007/s11770-007-0036-5 |
Zhao A H, Ding Z F. 2009. Earthquake location in transversely isotropic media with a tilted symmetry axis. J. Seismol., 13(2): 301-311. DOI:10.1007/s10950-008-9129-8 |
Zhao A H, Ding Z F, Bai Z M. 2015. Improvement of the ray-tracing based method calculating hypocentral loci for earthquake location. Chinese J. Geophys (in Chinese), 58(9): 3272-3285. DOI:10.6038/cjg20150922 |
Zhao A H, Ding Z F, Sun W G, et al. 2008. Calculation of focal loci for earthquake location in complex media. Chinese J. Geophys (in Chinese), 51(4): 1188-1195. |
Zhao A H, Zhang Z J. 2004. Fast calculation of converted wave traveltime in 3-D complex media. Chinese J. Geophys (in Chinese), 47(4): 702-707. |
Zhao A H, Zhang Z J, Teng J W. 2004. Minimum travel time tree algorithm for seismic ray tracing:improvement in efficiency. J. Geophys. Eng., 1(4): 245-251. DOI:10.1088/1742-2132/1/4/001 |
Zhou H W. 1994. Rapid three-dimensional hypocentral determination using a master station method. J. Geophys. Res., 99(B8): 15439-15455. DOI:10.1029/94JB00934 |
Zhou J C, Zhao A H. 2012. An intersection method for locating earthquakes in 3-D complex velocity models. Chinese J. Geophys (in Chinese), 55(10): 3347-3354. DOI:10.6038/j.issn.0001-5733.2012.10.017 |
白志明, 王椿镛. 2004. 云南遮放-宾川和孟连-马龙宽角地震剖面的层析成像研究. 地球物理学报, 47(2): 257-267. DOI:10.3321/j.issn:0001-5733.2004.02.012 |
蔡明军, 山秀明, 徐彦, 等. 2004. 从误差观点综述分析地震定位方法. 地震研究, 27(4): 314-317. DOI:10.3969/j.issn.1000-0666.2004.04.005 |
陈棋福, 张跃勤, 周静. 2001. 数字观测时代的全球三维结构与地震定位研究. 地震, 21(2): 29-40. |
陈运泰. 2009. 地震预测:回顾与展望. 中国科学D辑:地球科学, 39(12): 1633-1658. |
傅淑芳, 刘宝诚. 1991. 地震学教程. 北京: 地震出版社.
|
廉超, 李胜乐, 董曼, 等. 2006. 球面交切法地震定位. 大地测量与地球动力学, 26(2): 99-103. |
孟玉梅, 赵永, 王斌. 2001. 中国地震观测台网地震速报定位偏差的分析. 地震, 21(3): 65-69. |
齐诚, 赵大鹏, 陈顒, 等. 2006. 首都圈地区地壳P波和S波三维速度结构及其与大地震的关系. 地球物理学报, 49(3): 805-815. DOI:10.3321/j.issn:0001-5733.2006.03.024 |
宋锐, 顾小红, 王永力, 等. 2001. 国家数字地震台网中心实时处理及大地震速报软件系统. 地震, 21(4): 47-59. |
苏道磊, 范建柯, 吴时国, 等. 2016. 山东地区地壳P波三维速度结构及其与地震活动的关系. 地球物理学报, 59(4): 1335-1349. DOI:10.6038/cjg20160415 |
滕吉文, 张中杰, 白武明, 等. 2004. 岩石圈物理学. 北京: 科学出版社.
|
田玥, 陈晓非. 2002. 地震定位研究综述. 地球物理学进展, 17(1): 147-155. DOI:10.3969/j.issn.1004-2903.2002.01.022 |
王椿镛, 杨文采, 吴建平, 等. 2015. 南北构造带岩石圈结构与地震的研究. 地球物理学报, 58(11): 3867-3901. DOI:10.6038/cjg20151101 |
王辉, 常旭. 2000. 基于图形结构的三维射线追踪方法. 地球物理学报, 43(4): 534-541. DOI:10.3321/j.issn:0001-5733.2000.04.014 |
王璐琛, 常旭, 王一博. 2016. 干涉走时微地震震源定位方法. 地球物理学报, 59(8): 3037-3045. DOI:10.6038/cjg20160826 |
许力生, 杜海林, 严川, 等. 2013a. 一种确定震源中心的方法:逆时成像技术(一)-原理与数值实验. 地球物理学报, 56(4): 1190-1206. DOI:10.6038/cjg20130414 |
许力生, 严川, 张旭, 等. 2013b. 一种确定震源中心的方法:逆时成像技术(二)-基于人工地震的检验. 地球物理学报, 56(12): 4009-4027. DOI:10.6038/cjg20131207 |
杨文东, 金星, 李山有, 等. 2005. 地震定位研究及应用综述. 地震工程与工程振动, 25(1): 14-20. DOI:10.3969/j.issn.1000-1301.2005.01.003 |
于湘伟, 陈运泰, 张怀. 2010. 京津唐地区地壳三维P波速度结构与地震活动性分析. 地球物理学报, 53(8): 1817-1828. DOI:10.3969/j.issn.0001-5733.2010.08.007 |
曾融生. 1991. 大陆岩石圈构造与地球动力学. 地球科学进展, 6(2): 1-10. |
赵爱华, 丁志峰, 白志明. 2015. 基于射线追踪技术计算地震定位中震源轨迹的改进方法. 地球物理学报, 58(9): 3272-3285. DOI:10.6038/cjg20150922 |
赵爱华, 丁志峰, 孙为国, 等. 2008. 复杂介质地震定位中震源轨迹的计算. 地球物理学报, 51(4): 1188-1195. DOI:10.3321/j.issn:0001-5733.2008.04.029 |
赵爱华, 张中杰. 2004. 三维复杂介质中转换波走时快速计算. 地球物理学报, 47(4): 702-707. DOI:10.3321/j.issn:0001-5733.2004.04.023 |
周建超, 赵爱华. 2012. 三维复杂速度模型的交切法地震定位. 地球物理学报, 55(10): 3347-3354. DOI:10.6038/j.issn.0001-5733.2012.10.017 |