2015, Vol. 58
2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
We migrate every record of the station backward and take the position where all the imaging results merge relatively as the source location that releases the dominant energy. We propose to use the numerical scheme of varying grid. Firstly we utilize the method of global search to obtain the quick result of the epicenter in the coarse grid condition, and then we determine the final source location just by the local search within the epicenter area with a fine grid. We also employ a 3-D spatial Gaussian filter to obtain the source location and suppress the background noise automatically when the imaging results are scattered or unfocused, getting rid of searching for the major spatial energy points in the conventional least square method.
The results based on the synthetic data show that the method combining the varying-grid 3-D Gaussian beam and spatial Gaussian filter can not only ensure the high efficiency with less steps of ray tracing but also can be more adaptable to the low signal-to-noise ratio data and scattered or unfocused imaging energy. Besides, we also test the proposed method on the data of Zhuolu, Luanxian and Fangshan earthquakes in the Capital Circle region around Beijing, which can further confirm the efficiency and accuracy of the method.
The earthquake location method based on 3-D varying-grid Gaussian beam migration imaging just needs the rough picking of first arrivals, which is more robust to the pickup error than conventional methods. The method implements the unification of imaging accuracy and computational efficiency, and is highly automatic and suitable for noisy data especially in low signal-to-noise ratio cases.
1 引言
地震定位是地震学研究的基本问题之一，对于研究地震活动构造、地球内部结构等地震学问题具有重要意义(蔡明军等，2004；田玥和陈晓非，2002). 其主要应用领域包括天然地震定位、水压致裂检测以及水库的抗震特性研究等方面. 最基本的地震定位方法是盖格法(Geiger，1912)，其后逐渐发展起主事件法(Fitch，1975；Fukao，1972)、双差法(黄耘等，2008；Waldhauser and Ellsworth， 2000；杨智娴等，2003)、交切法(廉超等，2006；周建超和赵爱华，2012)以及网格搜索法(Shearer，1997；许力生等，2013；姚志祥等，2010)等定位方法. 这些传统的地震定位方法都运用了体波震相的到时信息，根据选用体波震相的不同主要分为以下两种：一种是利用P波和S波的到时差计算震源距离，从而确定震源位置(Lay and Wallace， 1995；周建超和赵爱华，2012);另一种是只选择P波，将每个台站的P波观测走时与理论走时残差的绝对值最小作为目标函数，通过一定的反演算法进行震源定位(高星等，2002；田玥和陈晓非，2006). 上述这些常规的地震定位方法对P波和S波的到时拾取具有很强的依赖性，当初至不明显或者被噪声淹没时定位精度很低.
近些年来，利用地震偏移成像技术进行震源定位的方法逐渐发展起来(Baker et al., 2005；Kao and Shan， 2004，2007；McMechan，1982；McMechan et al., 1985；Rentsch et al., 2007)，其优点是不严重 依赖于初至的拾取精度，只需大致获取地震发生时 间附近的波场就可以进行震源定位. McMechan等人(McMechan，1982；McMechan et al., 1985)利用逆时偏移方法对震源进行定位，根据波动方程通过时间的反向延拓获得震源位置，并对1983年美国加州地区长谷火山口的地震事件进行了测试；Kao和Shan(2004，2007)用发展的震源扫描算法(SSA)对2003年美国卡斯卡迪亚地区的地震进行了试算；Baker等人(2005)利用克希霍夫重建法(RKL)反传所有台站不同时间步长的P波振幅在空间网格点进行成像，并成功应用于美国南加州地区的地震定位. 逆时偏移方法在三维情形下的计算量十分庞大，不利于实际应用，克希霍夫法在速度模型存在较强的聚焦效应时精度很差；相比之下，高斯射线束方 法(erven y ′ et al.， 1982； erveny and Peník，1984； erven y ′ 1985，2001； Hale，1992；Hill， 1990，2001；李振春等，2010；Popov，1982；岳玉波等， 2010，2012)是精度和效率的合理折衷，更加适合于非均匀介质中的地震成像问题，且由于高斯束具有一定的射线宽度，与传统克希霍夫法相比可以在同等数据的前提下增加成像空间的照明度. Rentsch等人(2007)利用高斯束偏移对水压致裂问题进行了微震定位研究，并采用极化分析方法有效地节约了计算量.
本文将高斯射线束偏移方法用于天然地震的定位问题. 尽管高斯射线束方法具有较高的计算效率，但在三维大尺度情形下仍然需要进一步的改进才能更好地用于实际地震定位. 为了大幅度地提高三维高斯射线束的执行效率，Rentsch等人(2007)在对水压致裂的微震定位中采用了极化分析方法，通过预估射线方位角的取值范围以节约计算量，其主要缺陷是当信噪比较低时极化分析所确定的方位 角范围可靠性较差. 本文提出了一种基于变网格的高斯射线束地震定位方法，首先利用高斯射线束在粗网格下预估震源的大致位置，然后在置信的方位角范围内再利用高斯射线束在精细网格上进行局部搜索. 这种变网格的高斯射线束方法不但具有很高的计算效率，且其方位角估计比较稳健，主要受粗网格尺度及震中距大小影响，在处理低信噪比资料时较极化分析方法更为可靠和有效. 此外，本文对成像后的结果采用了空间高斯滤波(Oppenheim et al., 1999)，该方法不但可降低噪声对成像结果的影响，而且不需要逐点拾取成像能量的各个交点位置，能够实现震源位置的自动获取. 2 方法
高斯射线束法是波动方程的高频近似解，该方法能够准确刻画波场的运动学和动力学特征，克服了传统射线方法在阴影区和焦散区无法聚焦成像的缺陷，且具有较高的成像精度和计算效率(Hale，1992；Hill，1990)，因此在复杂介质地震偏移成像领域得到了广泛应用. 2.1 三维高斯射线束表达式
erven y' 等(1982)从弹性动力学方程出发，推导出了三维射线中心坐标系下高斯射线束的高频近似解：
，s为沿着中心射线 e t方向的射线路径长度，ω为高斯射线束的频率，τ0为中心射线追踪的初始时刻，n，m为垂直于中心射线平面上 e n和e m方向上的坐标分量. M = PQ －1为一个2×2的复对称矩阵，称为动力学参数矩阵，它表征着高斯射线束的特征，M 的实部特征值表征了高斯射线束的波前曲率，虚部则决定了垂直于中心射线平面上的振幅分布. P，Q 满足：
，矩阵元素vnn，vnm，vmn，vmm分别为速度沿 e n和e m方向的二阶偏导数.图 1为建立在直角坐标系中的射线中心坐标系，e t，e n和e m分别表示沿中心射线和垂直于中心射线的单位矢量. 从震源出发的以中心射线为中心的能量管组成了高斯射线束，射线束的能量分布按偏离中心射线的方向呈高斯型衰减. 射线中心坐标系的三个基矢量 e t，e n，e m采用如下的方式来定义：
其中θ为射线与竖直方向的夹角，为射线出射的方位角.
 | 图 1 射线中心坐标系中的高斯射线束 Fig. 1 Gaussian beam in ray-centered coordinate |
在三维各向同性介质中，设 x s=(xs，ys，zs)为震源位置，x r=(xr，yr)为地表地震台站位置，将地表台站数据向下反向延拓至 x =(x，y，z)处的地震波场为u(x，x s;ω)，Hill(2001)根据RayleighⅡ积分公式得到如下表达式：
其中，
，从台站 x r=(xr，yr)反向延拓到空间任一点 x 的三维格林函数可以表示为：
分别为中心射线在x，y，z三个方向上的射线参数. Hill(2001)为协调算法的精度与计算效率，通过引入相位校正因子，此时格林函数可由点 x 附近的射线束中心位置 L =(Lx，Ly，0)出射的高斯束UGB(x，x r，p ;ω)叠加来近似表示：
当 x r和L 距离较远时，公式(6)会存在一定的数值误差. 可以通过对地表观测台站加入一系列重叠的高斯窗来减小误差，高斯窗的中心即为射线束的中心位置，高斯函数的性质如下：
将地表台站的所有记录作为初始条件和边界条件，以一定的射线密度进行反传播，利用三维各向同性介质的高斯射线束叠后偏移成像公式
在处理实际地震数据时存在以下三个主要问题：(1)天然地震的发震时刻是未知的；(2)需要计算的空间尺度较大，直接采用传统的高斯射线束定位的计算效率较低；(3)P波的辐射机制会使各台站数据的振幅与相位随方位角变化，成像能量不一致，且定位处理后能量交点往往不是唯一的. 为了解决上述问题，提出了如图 2所示的变网格高斯射线束定位方法：(1)初步拾取各台站的P波和S波初至，从而可预估一个初始的发震时刻，对该时刻附近某时间窗内的时间点进行循环，对各台站的数据进行归一化，在粗网格下利用高斯射线束通过全局搜索的方法快速定位，确定震源的大致位置和发震时刻；(2)将上述获得的发震时刻和震源位置作为基准，在较为精细的时间窗和细网格情形下进行局部方位角的搜索，再对不同时刻的成像结果进行空间高斯滤波，能量最大的位置被视作最终的震源位置，所对应的时刻即为最终确定的发震时刻.
 | 图 2 利用变网格高斯射线束法定位的流程图 Fig. 2 The flow chart of location using Gaussian beam migration with varying mesh |
本文采用以下方法粗略地确定地震的发震时刻. 从地震记录中初步地拾取直达P波和S波震相到时，根据和达法有：
图 3为变网格高斯射线束定位方法的示意图. 图中黑色三角表示地震台站的位置，黑色五角星的位置为粗网格下通过全局搜索的方法快速计算的震中位置. 由数值测试的结果可知，在速度模型准确 的前提下，利用高斯射线束定位的数值误差不会超过网格大小的一半，选取 2 倍的粗网格大小作为误差区间即可保证震中位于绿色区域内. 然而在处理实际数据时，需要考虑不精确的速度模型带来的误差，一般取10倍的粗网格大小作为可信的震中区域. 基于确定的震中区域，再利用细网格对该区域进行局部搜索，仅需搜索大小的方位角即可，同理，我们也可以较为准确地约束倾角的取值范围. 同传统极化分析约束方位角与倾角的方法相比，该方法在处理信噪比较低数据时更为有效与可靠. 变网格的高斯束定位方法可以大幅度提高计算效率，提高的效率约为粗细网格比的三次方，本文在处理实际数据时采用的粗细网格比为5，故较直接利用三维高斯束定位的效率可提高约100倍.
 | 图 3 变网格高斯射线束定位法的示意图，粗线为粗网格，细线为细网格 Fig. 3 Sketch map of Gaussian beam migration with varying mesh |
本文为了解决成像空间能量不能汇聚的问题，引入空间高斯滤波函数(Oppenheim，1999)：
其中，σx，σy，σz决定了高斯函数在x，y，z三个方向的宽度. 设I(x，y，z;ti)为某个时刻P波记录在三维空间的成像结果，对该成像结果采用空间高斯滤波，设J(x，y，z;ti)为经过高斯滤波之后的成像结果，则有：
对时间窗口内的所有时刻获得的高斯滤波三维成像结果J(x，y，z;ti)(i=1，2，3，…，n)进行遍历，当某个时间点处的能量取得最大值时，则将该点定为最终的震源位置，相应的时间ti即为地震的发震时间.
空间高斯滤波不但可以使成像空间的能量汇聚于一点，而且可以去除噪声能量的干扰，为程序自动提取震中位置创造了条件. 图 4a为二维空间中一系列台站成像的能量交点，左上角五个较为集中的点表示可信度较高的能量交点，右下角的孤立点表示可信度较差的能量交点. 图 4b为对图 4a的二维空间进行高斯滤波得到的能量分布结果；图 4c中“十”字表示通过高斯滤波获得的有效能量中心，三角形表示使用最小二乘法获得的有效能量中心. 从图中可以看出，传统的最小二乘法易受到噪声点的干扰，而采用高斯滤波则可有效地消除噪声能量干扰，使空间能量更为集中. 此外，高斯滤波直接对二维空间进行操作，避免了传统最小二乘法需要对空间主要能量点的搜索与选取，因而自动化程度较高. 通过选取合适的σx，σy和σz大小可以有效地控制高斯滤波窗口的宽度，进而在压制噪声和提高空间分辨率之间达到平衡.
 | 图 4 高斯函数滤波示意图 (a)二维空间离散的能量点;(b)经过高斯函数滤波的能量分布;(c)高斯函数滤波和最小二乘法结果对比. Fig. 4 Sketch of Gaussian-function filter (a)Discrete energy points in 2-D space；(b)Energy distribution after processed by Gaussian-function filter； (c)Results comparison between the Gaussian-function filter and the least square method. |
为了验证高斯射线束偏移成像用于震源定位的准确性，我们采用理论地震图进行数值试验. 选取符合天然地震学基本假设的层状凹陷模型，如图 5a所示为速度模型的X-Z剖面，以该剖面的垂向中心为对称轴旋转360°构建三维速度模型. 模型三维尺寸(长×宽×高)为20 km×20 km×15 km. 震源为 双力偶源，震源时间函数采用主频为10 Hz的 Ricker子波，震源位置为x=10.0 km，y=10.0 km，z=10.0 km. 高斯射线束选用的参数为：(1)粗网格时，ω0=2π×10 rad·s-1，ωh=2π×50 rad·s-1，wl=400 m；(2)细网格时，ω0=2π×40 rad·s-1，ωh=2π×200 rad·s-1，wl=200 m. 如图 5b，5c所示为地震震中和地表台站的分布情况以及相应的Z分量波形记录，其中五角星代表震中位置，三角形为台站的分布位置.
 | 图 5(a)速度模型的X-Z剖面；(b)地表台站及震中分布;(c)理论波形记录 Fig. 5(a)X-Z profile of the velocity model;(b)Distribution of the stations and epicenter;(c)Synthetic waveform records |
选取Z分量中的P波数据进行试算，定位结果如图 6所示. 图 6a—6c为利用常规的高斯射线束方法进行定位获得的过震源中心处三个方向的成像 切片，确定的震源位置为(10.0 km，10.0 km，10.0 km)，与理论测试的震源位置高度一致. 图 6d—6f为单独利用变网格的高斯射线束进行定位的结果，有效地减少了射线束的搜索数量，并使成像能量集中在有效范围内；图 6g—6i 为在变网格法基础上加入高斯滤波后获得的最终计算结果，与常规的高斯射线束定位方法相比，不但成像结果更为集中，有效地避免了成像结果中聚焦点不唯一的情况，而且大幅度提高了计算效率，在大尺度的天然地震定位中的优势会更加明显.
 | 图 6 理论合成数据的成像切片图
(a)—(c)为常规高斯射线束的成像结果;(d)—(f)为单独利用变网格法高斯射线束的成像结果;(g)—(i)为在变网格法基础上采用高斯滤波的成像结果. Fig. 6 Imaging slices with synthetic data (a)—(c)Imaging results calculated by conventional Gaussian beams;(d)—(f)Imaging results calculated by varying-grid Gaussian beams method alone;(g)—(i)Imaging results calculated by varying-grid method joint Gaussian filter. |
 | 图 7 加入白噪声的理论合成数据 Fig. 7 Synthetic data with white noise |
如图 7所示为加入噪声的合成数据，在每个台站的数据中均加入了幅值为该台站P波震相最大幅值30％的随机噪声，其他条件保持不变. 由于数据的噪声水平较高，P波的初至基本完全被淹没在噪声之中，此时，利用极化分析的手段获的方位角信息极易受到时窗选择的影响，会直接影响到震源位置定位的精度和可靠性. 由图 8可见，采用变网格的高斯束定位方法并对成像结果进行高斯滤波，依然可以获得准确的震源位置，并能够有效地去除噪声能量的干扰.
 | 图 8 含噪声数据的成像切片 (a)—(c)为常规高斯射线束的成像结果;(d)—(f)为改进的高斯射线束的成像结果. Fig. 8 Imaging slices with noise-mixed data (a)—(c)Imaging results calculated by conventional Gaussian beams;(d)—(f)Imaging results calculated by modified Gaussian beams. |
为了验证网格间距对高斯射线束定位精度的影响，进行如下数值测试. 选取如图 5a所示的 速度模型，将震源位置放在x=8.1 km，y=9.2 km，z=9.5 km处，网格间距分别选取为0.05 km，0.1 km，0.2 km，0.3 km，0.4 km.
如表 1所示，当震源的深度及震中位置坐标恰好位于网格点上，即dx=0.05 km和0.1 km时，利用高斯射线束偏移成像获得的震源位置与理论位置完全吻合，且利用不同网格间距的定位结果相同，这表明在一定程度上本方法对于粗网格情形下的定位依然具有较高的计算精度. 当震源位置不刚好位于网格点上，即dx=0.2 km，0.3 km，0.4 km时，利用高斯射线束偏移成像定位的结果会出现不超过半个网格间距的数值定位误差，这在大尺度的天然地震 定位中是可以接受的. 需要特殊说明的是，本数值测试的目的是在假设速度模型准确的前提下，验证方法本身的可靠性，但在处理实际数据时，由于速度模型是未知的，此时本文方法与其他地震定位方法面临同样的问题，即不能解决由速度模型的不准确所带来的误差问题.
|
表 1 不同网格间距的计算结果对比 Table 1 Comparison of the calculation results with different grid spacing |
为了验证数据的噪声水平对极化分析的影响，并针对计算效率将极化分析法与变网格法作对比，本文进行了如下数值测试. 利用如图 5b所示的观测系统，采用均匀半无限空间速度模型，其中 P波的速度为6.0 km·s-1，S波速度为3.5 km·s-1，介质密度为2.7×103 kg·m-3，震源深度为10.0 km，利用广义反射透射系数矩阵法(Yao and Harkider， 1983)计算出三分量理论地震图. 采用极化分析的方法分别对理论数据及含30％白噪声的数据进行极化分析，同时利用变网格法搜索含噪声数据的方位角与倾角信息，并给出相应的计算结果，如表 2所示.
|
|
表 2 方位角与倾角范围结果对比 Table 2 Comparison of the determination of the azimuth and incidence angles |
由表 2可知，当存在较强的噪声情况下，极化分析方法不能有效地约束方位角和倾角的信息，采用粗网格进行全局搜索可将角度信息约束在一个较小的范围内，表明利用变网格的高斯束定位方法在计算效率和成像精度上实现了较好的统一. 4 应用实例
根据中国地震台网中心(CENC)公布的2005年至2010年间首都圈地震数据，震源深度基本都在25 km以内，属于浅震，一般位于上地壳，因此本文 使用30 km深的ISPA91速度模型(表 3)对河北涿鹿、滦县以及北京房山等地区的三次地震事件进行 试算，并与地震局系统采用的JOPENS自动定位软件(利用“单纯型法”)的定位结果 作对比. 高斯射线束选用的参数为：(1)粗网格情况下，ω0=2π×5 rad·s-1，ωh=2π×10 rad·s-1，wl=2000 m；(2)细网格情况下，ω0=2π×20 rad·s-1，ωh=2π×80 rad·s-1，wl=400 m.
|
|
表 3 ISPA91速度模型和PREM速度模型 Table 3 ISPA91 velocity model and PREM velocity model |
如表 4所示，选择2013年1月份发生在涿鹿、滦县以及房山的三次地震事件，从150 km左右的震中距范围内挑选台站数据，通过高斯变换将经纬度投影为平面坐标，对计算区域进行网格化，利用变网格的三维高斯射线束定位. 本文以河北涿鹿地震(2013-01-04)为例，给出其台站分布(图 9)、地震波 形记录(图 10)以及利用本方法进行定位的结果(图 11).
|
|
表 4 实际地震数据信息及相关参数 Table 4 Information of the real seismic data and related parameters |
 | 图 9 涿鹿地震的台站分布 Fig. 9 Distribution of stations in selected Zhuolu earthquake |
 | 图 10 涿鹿地震实际地震资料 图中虚线表示利用和达法确定的发震时刻范围. Fig. 10 Real seismic data in Zhuolu earthquake Dashed lines in the picture indicate the range of the origin time of earthquake with Heda method. |
 | 图 11 涿鹿地震成像结果切片 (a)Z=19.2 km处，X-Y切片上的成像结果；(b)X=90.0 km处，Y-Z切片上的成像结果；(c)Y=89.6 km处，X-Z切片上的成像结果；(d)在震源位置的X-Y切片上，当X=90.0 km时，Y方向上的能量分布；(e)在震源位置的X-Y切片上，当Y=89.6 km时，X 方向上的能量分布；(f)在震源位置的X-Z切片上，当X=90.0 km时，Z方向上的能量分布. Fig. 11 Slice of imaging space in Zhuolu earthquake (a)at Z=19.2 km，the imaging result on the X-Y slice;(b)at X=90.0 km，the imaging result on the Y-Z slice;(c)at Y=89.6 km，the imaging result on the X-Z slice;(d)On the X-Y slice of the source location，the distribution of the energy in Y direction at X=90.0 km;(e)On the X-Y slice of the source location，the distribution of the energy in X direction at Y=89.6 km;(f)On the X-Z slice of the source location，the distribution of the energy in Z direction at X=90.0 km. |
通过计算，得到涿鹿地震的发震时刻为:2013年1月4日，07 ∶ 12 ∶ 13.6. 由图 11a—11c可见，震中位置为(90.0 km，89.6 km)，相应的经纬度为115.437°E，40.247°N，震源深度为19.2 km.
为了获取震源位置的最优估计，我们对成像空间采用高斯滤波，而高斯滤波函数的应用使成像空间各个能量点的位置误差呈正态分布. 根据统计模型的基本假设，在传统地震定位中，若理论走时和实际走时之间的误差满足正态分布，则可用置信区间的范围来描述定位误差(Flinn，1965)，我们沿用该基本假设，以最大能量点的95%作为误差容许范围来评价定位结果. 图 11d—11f分别为震源定位结果在X、Y、Z三个方向上的误差分布图. 横虚线表示95%的置信区间，两条竖虚线之间的距离表示置信区间内的误差范围. 可见，Y方向上的最大能量点位置为Y=89.6 km，误差为±2.6 km；X方向上 的最大能量点位置为X=90.0 km，误差为±2.2 km； Z方向上的最大能量点位置为Z=19.2 km，误差为±3.2 km.
为了验证本文方法对涿鹿、滦县以及房山三次地震事件定位的适用性与准确性，表 5给出了JOPENS利用“单纯型法”和本文方法的定位结果对比. 结果表明，两种定位方法的结果相近：发震时刻相差1.0 s左右，震中相差4.5 km以内，深度相差 2 km以内，总体上来讲，两种方法的计算结果具有 较好的一致性，虽然二者的震源参数在空间位置上存在一定的差异，但位于合理范围内. 本文的方法选取P波和S波最大振幅作为P波和S波的到时，而JOPENS自动定位系统基于传统的P波与S波初至到时. P波初至与P波最大振幅到时之间的不一致性会使二者的计算结果存在一定的差异. 本文采用的方法和JOPENS定位系统都需要对地震事件有一定的识别度. 当噪声水平较高时，拾取P波和S波初至有较大的难度，而本方法采用P波最大振幅来确定震源位置，不易受背景噪声对初至拾取 的影响，因而更加适合于低信噪比的微地震定位问题.
|
|
表 5 两种方法计算的三次地震的震源参数比较 TableTable 5 Comparison of seismic source parameters in three earthquakes with two methods |
速度结构是影响定位准确性的一个重要因素，无论是盖戈法还是双差法，其定位结果都强烈地依赖所使用的速度模型(许力生等，2013). 由于真实的速度模型往往是未知的，为了验证速度模型存在误差时，高斯射线束方法的定位能力，仍以河北涿鹿地震(2013-01-04)为例，给出采用PREM速度模型(参见表 3)的定位结果，如图 12所示.
 | 图 12 基于PREM模型的涿鹿地震成像结果切片 (a)Z=22.0 km处，X-Y切片上的成像结果；(b)X=90.4 km 处，Y-Z切片上的成像结果；(c)Y=89.6 km处，X-Z切片上 的成像结果. Fig. 12 Slice of imaging space based on PREM velocity model in Zhuolu earthquake (a)At Z=22.0 km，the imaging result on the X-Y slice;(b)At X=90.4 km，the imaging result on the Y-Z slice;(c)At Y=89.6 km，the imaging result on the X-Z slice. |
基于PREM速度模型，得到涿鹿地震的发震时 刻为: 2013年1月4日，07 ∶ 12 ∶ 13.6. 由图 12可见，震中位置为(90.4 km，89.6 km)，相应的经纬度为115.440°E，40.245°N，震源深度为22.0 km. 与基于ISPA91速度模型获得的震源时空参数相比，二 者的发震时刻相同，震中位置相差0.4 km，震源深度相差2.2 km. 由此可见，当速度模型存在误差时，高斯射线束也能较好进行定位，其中震中受到的影响相对较小，震源深度对速度模型的误差较为敏感，但均位于误差允许的范围内. 尽管如此，定位结果的准确性与速度模型的精确程度仍然息息相关，高斯射线束的定位结果仍需要一个尽可能准确的速度模型作为保证. 5 讨论和结论
本文旨在初步探讨利用高斯射线束进行区域地震定位的可行性，而一个准确可靠的定位结果同时受观测资料、方程系统、速度模型以及求解方法的制约(许力生等，2013).通过数值算例和实际资料的处理结果可知，在可信度较高的速度模型基础上，本文提出的方法具有较好的成像精度. 本文在初步确定发震时刻时，同时使用P波和S波到时信息，当S波的到时信息难以获取时，可以单纯使用P波，借助常规方法初步地确定震源时空参数，在此基础上再利用细网格的高斯射线束定位，结果将更为准确可靠，在实际问题中的应用将更为便利. 本文在偏移成像定位过程中只采用了三分量中Z分量的P波数据，下一步可以综合利用三分量的弹性波数据进行P-S波联合定位，有望获得更多的震源信息与更高的定位精度.
在计算效率方面，变网格的三维高斯射线束在 对近震进行地震定位时具有较高的计算效率，但在处理远台站数据时，高斯射线束方法具有先天的不足，即高斯射线束在传播过程中的有效宽度会随着 距离的增大而急剧增加(erven y ′ et al.， 1982)，将会同时影响到计算效率和成像精度，因此如何使高斯射线束快速高精度地远震定位将是未来的主要研究方向之一. 由于高斯射线束在处理复杂介质时比传统射线方法具有更好的适应性，比逆时偏移方法具有更高的计算效率，在下一步的工作中，将主要针对具体小区域的水压致裂或矿山微破裂等定位问题，根据测井数据或偏移成像中的速度分析手段，将会获得高精度的速度模型，此时高斯射线束在三维复杂速度模型中的优势将会得到更好的体现. 与此同时，在处理天然地震数据时，采用相对简单的高斯滤波函数即可实现较好的震源位置估计，在处理三维复杂速度模型时，高斯滤波函数的选取需要进一步考虑速度和台站分布等因素的影响.
本文提出了一种利用变网格的三维高斯射线束对天然地震震源进行定位的方法，同传统的基于极化分析的方法(Rentsch et al., 2007)相比，在保证具有较高计算效率的同时提高了对低信噪比数据的适应性. 通过对预估发震时刻附近的时窗扫描可以获得相对准确的发震时刻. 提出了空间高斯滤波法对各台站成像交点散乱分布情形下的震源位置估计，既能极大地提高方法的自动化程度，自动剔除偏离过大的交点，又能有效压制背景噪声对常规初至拾取带来的影响. 该方法只需对波形初至到时进行粗略拾取，受到初至拾取误差的干扰明显小于常规方法，因而压制噪声的能力更强. 理论模型的测试和涿鹿、滦县、房山等地区实际地震资料的处理证实了方法的有效性.
致谢 感谢中国地震局地球物理研究所提供的首都圈地震台网数据，感谢广东省智源工程抗震科技公司开发的JOPENS软件为本文的方法验证提供了参考，感谢王卫民、赵连锋和郝金来与作者的有益讨论和建议.| [1] | Baker T, Granat R, Clayton R W. 2005. Real-time earthquake location using Kirchhoff reconstruction. Bull. Seism. Soc. Am., 95(2): 699-707. |
| [2] | 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. |
| [3] | erveny' V, Popov M M, Peník L. 1982. Computation of wave fields in inhomogeneous media: Gaussian beam approach. Geophys. J. R. Astr. Soc., 70(1): 109-128. |
| [4] | erveny' V, Peník L. 1984. Gaussian beams in elastic 2-D laterally varying layered structures. Geophys. J. R. Astr. Soc., 78(1): 65-91. |
| [5] | erveny' V. 1985. Gaussian beam synthetic seismograms. J. Geophys., 58: 44-72. |
| [6] | erveny' V. 2001. Seismic Ray Theory. Cambridge: Cambridge University Press. |
| [7] | Fitch T J. 1975. Compressional velocity in source regions of deep earthquakes: an application of the master earthquake technique. Earth Planet. Sci. Lett., 26(2): 156-166. |
| [8] | Flinn E A. 1965. Confidence regions and error determinations for seismic event location. Rev. Geophys., 3(1): 157-185. |
| [9] | Fukao Y. 1972. Source process of a large deep-focus earthquake and its tectonic implications—the western Brazil earthquake of 1963. Physics of the Earth and Planetary Interior, 5: 61-76. |
| [10] | Gao X, Wang W M, Yao Z X. 2002. Hypocentral determination using simulated annealing method. Chinese J. Geophys. (in Chinese), 45(4): 525-532. |
| [11] | Geiger L. 1912. Probability method for the determination of earthquake epicenters from arrival time only. Bull. St. Louis. Univ., 8: 60-71. |
| [12] | Hale D. 1992. Migration by the Kirchhoff, slant stack and Gaussian beam methods. Colorado School of Mines Center for Wave Phenomena Report, 121. |
| [13] | Hill N R. 1990. Gaussian beam migration. Geophysics, 55(11): 1416-1428. |
| [14] | Hill N R. 2001. Prestack Gaussian-beam depth migration. Geophysics, 66(4): 1240-1250. |
| [15] | Huang Y, Li Q H, Zhang Y S, et al. 2008. Relocation of earthquakes in Jiangsu and neighboring areas, China and analysis of structural features. Chinese J. Geophys. (in Chinese), 51(1): 175-185. |
| [16] | Kao H, Shan S J. 2004. The source-scanning algorithm: Mapping the distribution of seismic sources in time and space. Geophys. J. Int., 157(2): 589-594. |
| [17] | Kao H, Shan S J. 2007. Rapid identification of earthquake rupture plane using Source-Scanning Algorithm. Geophys. J. Int., 168(3): 1011-1020. |
| [18] | Lay T, Wallace T C. 1995. Modern Global Seismology. New York: Academic Press Inc.. |
| [19] | Li Z C, Yue Y B, Guo C B, et al. 2010. Gaussian beam common angle preserved-amplitude migration. Oil Geophysical Prospecting (in Chinese), 45(3): 360-365. |
| [20] | 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. |
| [21] | McMechan G A. 1982. Determination of source parameters by wavefield extrapolation. Geophys. J. R. Astr. Soc., 71(3): 613-628. |
| [22] | McMechan G A, Luetgert J H, Mooney W D. 1985. Imaging of earthquake sources in Long Valley Caldera, California, 1983. Bull. Seism. Soc. Am., 75(4): 1005-1020. |
| [23] | Oppenheim A V, Schafer R W, Buck J R. 1999. Discrete-Time Signal Processing. 2nd ed. Upper Saddle River, NJ: Prentice-Hall. |
| [24] | Popov M M. 1982. A new method of computation of wave fields using Gaussian beams. Wave Motion, 4(1): 85-97. |
| [25] | Rentsch S, Buske S, Luth S, et al. 2007. Fast location of seismicity: A migration-type approach with application to hydraulic-fracturing data. Geophysics, 72(1): S33-S40. |
| [26] | Shearer P M. 1997. Improving local earthquake locations using the L1 norm and waveform cross correlation: Application to the Whittier Narrows, California, aftershock sequence. J. Geophys. Res.: Solid Earth, 102(B4): 8269-8283. |
| [27] | Tian Y, Chen X F. 2002. Review of seismic location study. Progress in Geophysics (in Chinese), 17(1): 147-155, doi: 10.3969/j.issn.1004-2903.2002.01.022. |
| [28] | Tian Y, Chen X F. 2006. Simultaneous inversion of hypocenters and velocity using the quasi-Newton method and trust region method. Chinese J. Geophys. (in Chinese), 49(3): 845-854. |
| [29] | Waldhauser F, Ellsworth W L. 2000. A double-difference earthquake location algorithm: method and application to the northern Hayward fault, California. Bull. Seism. Soc. Am., 90(6): 1353-1368. |
| [30] | Xu L S, Du H L, Yan C, et al. 2013. A method for determination of earthquake hypocentroid: time-reversal imaging technique Ⅰ—Priniple and numerical tests. Chinese J. Geophys. (in Chinese), 56(4): 1190-1206, doi: 10.6038/cjg20130414. |
| [31] | Yang Z X, Chen Y T, Zheng Y J, et al. 2003. Accurate relocation of earthquakes in central-western China using double difference earthquake location algorithm. Science in China Series D: Earth Sciences, 46(S2): 181-188. |
| [32] | Yao Z X, Harkider D G. 1983. A generalized reflection-transmission coefficient matrix and discrete wavenumber method for synthetic seismograms. Bull. Seism. Soc. Am., 73(6A): 1685-1699. |
| [33] | Yao Z X, Wang C Y, Lou H. 2010. Ascertaining the structure parameters of Kunlun fault zone using the grid searching method based on trapped wave correlation. Chinese J. Geophys. (in Chinese), 53(5): 1167-1172, doi: 10.3969/j.issn.0001-5733.2010.05.018. |
| [34] | Yue Y B, Li Z C, Zhang P, et al. 2010. Prestack Gaussian beam depth migration under complex surface conditions. Applied Geophysics, 7(2): 143-148. |
| [35] | Yue Y B, Li Z C, Qian Z P, et al. 2012. Amplitude-preserved Gaussian beam migration under complex topographic conditions. Chinese J. Geophys. (in Chinese), 55(4): 1376-1383, doi: 10.6038/j.issn.0001-5733.2012.04.033. |
| [36] | 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. |
| [37] | 蔡明军, 山秀明, 徐彦等. 2004. 从误差观点综述分析地震定位方法. 地震研究, 27(4): 314-317. |
| [38] | 高星, 王卫民, 姚振兴. 2002. 用于地震定位的SAMS方法. 地球物理学报, 45(4): 525-532. |
| [39] | 黄耘, 李清河, 张元生等. 2008. 江苏及临区地震重新定位和构造特征分析. 地球物理学报, 51(1): 175-185. |
| [40] | 李振春, 岳玉波, 郭朝斌等. 2010. 高斯波束共角度保幅深度偏移. 石油地球物理勘探, 45(3): 360-365. |
| [41] | 廉超, 李胜乐, 董曼等. 2006. 球面交切法地震定位. 大地测量与地球动力学, 26(2): 99-103. |
| [42] | 田玥, 陈晓非. 2002. 地震定位研究综述. 地球物理学进展, 17(1): 147-155, doi: 10.3969/j.issn.1004-2903.2002.01.022. |
| [43] | 田玥, 陈晓非. 2006. 利用拟牛顿法和信赖域法联合反演震中分布与一维速度结构. 地球物理学报, 49(3): 845-854. |
| [44] | 许力生, 杜海林, 严川等. 2013. 一种确定震源中心的方法: 逆时成像技术(一)——原理与数值实验. 地球物理学报, 56(4): 1190-1206, doi: 10.6038/cjg20130414. |
| [45] | 杨智娴, 陈运泰, 郑月军等. 2003. 双差地震定位法在我国中西部地区地震精确定位中的应用. 中国科学D辑: 地球科学, 33(增刊): 129-134. |
| [46] | 姚志祥, 王椿镛, 楼海. 2010. 利用基于围陷波波形相关的网格搜索法确定昆仑山断裂带结构参数. 地球物理学报, 53(5): 1167-1172, doi: 10.3969/j.issn.0001-5733.2010.05.018. |
| [47] | 岳玉波, 李振春, 钱忠平等. 2012. 复杂地表条件下保幅高斯束偏移. 地球物理学报, 55(4): 1376-1383, doi: 10.6038/j.issn.0001-5733.2012.04.033. |
| [48] | 周建超, 赵爱华. 2012. 三维复杂速度模型的交切法地震定位. 地球物理学报, 55(10): 3347-3354, doi: 10.6038/j.issn.0001-5733.2012.10.017. |