射线追踪方法定位近震震源空间位置
张潜 尹耿 王玉石 林国良
摘要：介绍了一种基于地震波射线理论、不需要先验速度结构的定位震源空间位置的方法——射线追踪法。该方法利用地震射线会反向汇聚到震源的几何性质，对地面介质分层，并用网格搜索反演方法计算分层速度参数，进而确定震源位置；然后逐步增加分层，重复前面的步骤，对震源位置进行校正，最终得到一个较为精确的定位结果。通过虚拟事件对该方法进行了测试，测试结果表明，该方法理论上能够用于近震空间位置的定位；使用该方法对一个真实震侧进行定位，并和云南地震台网的定位的结果进行了比较，两者较为接近。
关键词：地震定位；无速度结构；射线理论；P波偏振
中图分类号：P3153文献标识码：A文章编号：1000-0666（2017）02-0203-08
0引言
震源空间位置的确定是地震学的经典问题——地震定位的一部分，也是地震学研究的基础。地震定位，无论是对于地球物理学研究，如地球内部的环境构造、地震的发生机制，还是对于震后救援工作的需求，如震害范围和程度的预估、地震趋势的预测以及地震预警，都有着无法替代的重要性。因此，地震定位的方法以及提高地震定位精度和速度的方法，一直都是地震学家的研究目标。
早期的地震定位主要是通过几何作图完成的。Geiger（1912）根据各个观测点到时差最小化的原则，将地震波走时方程组线性化，然后通过最小二乘法来定位地震。这也为以后使用计算机来定位地震提供了基础。20世纪70年代后，伴随着现代计算机技术的发展，地震学家们基于Geiger的理论，给出了一系列地震地位的程序和算法（Lee，Lahr，1975；Klein，1978；Lienert et al，1986；Nelson，Vidale，1990），并根据实际应用中出现的问题，提出了相应的改进方法（Lienert et al，1986；Prugger，Gendzwill，1988）。随后，多事件定位方法（Douglas，1976；Crosson，1976；Spence，1980；Waldhauser，Ellsworth，2000）、基于距离残差的空间域定位方法（Romney，1957；Lomnitze，1977）和非线性的定位方法（Prugger，Gendzwill，1988；Thurber，1985；Tarantola，Valette，1982）也相继出现。震相、振幅，甚至是全波形的信息都被用于震源位置的确定（Engdahl，et al，1998；Warren，Shearer，2005；罗艳等，2013；Zhao，Helmberger，1994；Zhu，Helmberger，1996）。
在地震定位的过程中，影响结果精度的因素大致有：观测系统的结构、速度结构的选取或假设、震相识别以及观测数据（包含观测点的空间位置）的读取和计算过程中产生的误差（朱元清，赵仲和，1997；田玥，陈晓非，2002）。绝大多数时间域和空间域的定位方法，都需要使用某一组先验的速度结构。因此，该种方法中速度结构的选取，将直接影响到结果的误差水平。而通过波形信息来确定震源位置时，震相识别、震源机制解和地球结构模型的可靠性，将决定定位结果的可靠性。
本文将介绍一种基于地震波射线理论，且不需要先验速度结构的定位近震震源位置的方法——射线追踪法，并使用虚拟数据和一个真实震例对该方法进行了检验。[HJ2mm]
1基本原理
假设地球介质是水平均匀且各向同性的分层介质。根据惠更斯原理和地震波射线理论，一束地震波在传播时，射线和分层界面法线的夹角满足方程[KH*2]
[SX（]v1sinθ1[SX）]=……=[SX（]vjsinθj[SX）]=……=[SX（]vnsinθn[SX）]=const[JY]（1）[KH*3/4D]
其中：vj和θj是射线在第j层的波速和夹角（震源所在位置为第n层，地表为第一层上边界）。将地表处的夹角和波速作为已知条件，则[KH*2]
θj=arcsin（[SX（]vjv1[SX）]sinθ1）[JY]（2）[KH*3/4D]
从（2）式可以看出，射線与法线的夹角只与地表入射角θn和当前波速与地表波速的比值vj/v1（记为rj）有关。另一方面，对于一束从震源发出的射线i，其在第j层介质中的传播，符合方程[KH*1]
sij=Δhj/tanθij[JY]（3）[KH*1D]
式中：sij和θij分别代表射线i在第j层介质中的水平位移和与法线的夹角；Δhj表示第j层的厚度。同时，对于射线i，满足方程[KH*1]
θij=arcsin（rjsinθi1）[JY]（4）[KH*1]
结合式（3），可得：[KH*1]
sij=Δhj/tan（arcsin（rjsinθi1））[JY]（5）[KH*1]
因此，震源到射线接收点（即观测点）之间的水平位移为：[KH*2]
[JP2]Si[KG-*4]=[KG-*2]∑[DD（]nj=1[DD）]sij[KG-*4]=[KG-*2]∑[DD（]nj=1[DD）]Δhj/tan（arcsin（rjsinθi1））[JY]（6a）[KH*2D][JP]
或
[JP3]Si[KG-*4]=[KG-*4]∫[DD（]H[KG*4]i0[DD）]si（h）dh[KG-*4]=[KG-*4]∫[DD（]H[KG*4]i0[DD）]1/tan（arcsin（r（h）sinθi1））dh[JP][KH-*4][JX*1][JY]（6b）[KH*1D]
其中：Hi代表震源深度。θi1和Pg波的偏振方向相关，可以通过观测得到；而Si和Hi代表着观测点和震源空间位置的矢量差。因此，我们可以得到关于震源空间位置F（x，y，h）和r（r1，r2，…rn）的方程（组）：[KH*1]
[JB（{][KF（]（xi-y）2+（yi-y）2[KF）]-Si=0[KG*2]（hi-h）-H=0[JB）][JY]（7）[KH*1D]
其中：（xi，yi，hi）为观测点i的空间位置。本文所介绍的方法，就是要找到一组（F，r），使（7）式成立，或者使等式兩边的差值尽量小，从而实现震源的定位。
[KG（0.15mm]由于式（6）和（7）含有三角函数，不能采用线性方程组的方法来求解，而必须寻找一种合适的反演方法来确定（F，r）。网格搜索法是一种不错的方法，但是其计算量会随着参数的增加而呈指数增长。如果待求解的参数越多，则计算所花时间越长。因此，我们采用了一种“分步”的网格搜索法来定位震源，待求参数变少，计算效率得到了提升。图1是对这种方法的一个简单的展示，由于看起来像追溯射线的路径回到震源，我们称该方法为“射线追踪法”，其操作步骤如下：[KG）]
（1）根据各个观测点接收到的Pg波的偏振方向，反向生成一组射线；寻找到各条射线距离最短的一点F1，并将其作为初始定位的结果（F1，r1）（r1=1）。
（2）将介质平均分成2层，每层的厚度为h1/2，r1取步骤（1）的结果。使用网格搜索法，寻找第二层的波速比r2，以及相应的点F2，使F2到根据（r1，r2）做出的反向射线的距离和最小。将（F2，r2）作为最新的结果；
（3）重复步骤（2），直到最后得到的结果（Fn，rn）满足一定的收敛条件。Fn即为最终的定位结果。
通过射线追踪法，我们一步步将介质的分层数从2层增加至n层，且每一次分层都只用搜索新增层的波速比，从而大大减少了计算量。最终，我们会得到震源的位置Fn以及1组和各层相对应的波速比（r1，r2，…，rn）。
2测试数据
为了验证此地震定位方法的准确性，本文虚拟了多组数据对该方法进行测试。测试数据生成[FL）][SD1，1]
[JZ][XC张潜1.TIF][KH*2/3D][HT5K][JZ]
图1射线追踪法步骤的简单演示[JZ]
Fig1Demonstration of ray tracing method[KH*1]
的方法如下：
（1）选择1个一维速度模型作为基准模型，假定震中位置为0，震源深度为10 km。
（2）随机生成一组出射角为（θ1n，…，θmn）的射线，根据步骤1的速度模型和震源位置，计算各条射线到达地表（深度为0）处的角度（θ11，…，θm1）和震中距（S1，…，Sm）。
（3）为（S1，…，Sm）和（θ11，…，θm1）引入误差，组成一组数据，分别对应m个虚拟观测点的坐标和接收到射线的Pg波偏振方向。
对于震中一定范围以外的观测点，Pg波的到达往往会晚于Pn波。此时，Pg波的振幅信息会和其它震相的信息混杂在一起，这就对Pg波偏振方向的确定造成了困难。因此，在生成测试数据时，要注意将θ值限定在一定的范围内。在本研究中，θ的取值范围分为3档，分别为≤ 30°、≤ 45°和≤ 60°，本文将对结果进行比较分析，以观察不同的出射角范围对定位精度的影响。同理，为了解观测误差和观测点数量对于定位精度的影响，本文也选取了相应的参数来进行测试。为了计算上的方便，我们将三维空间简化为二维，即去掉了一个水平维度，变成了如图1所示的情况。这样的简化并不会对测试造成本质性的影响。
本文共虚拟了14万组数据来进行测试，其详细信息见表1。所有的数据按照不同的观测点数量、出射角度范围、误差范围等因素分入14个大组，每个大组包含1万组数据。除观测点数量以外，生成每一组数据所用到的角度和误差都是在限定范围内随机生成的。显然，当误差水平为0，则意味着该组数据没有引入相应的误差。
3测试结果及比较
我们虚构了一个震中位置为0，震源深度为10 km的震源，并根据随机生成的参数生成了观测数据，最终反演得到了14万组定位结果。对这些结果进行比较，以观察不同的条件对于反演结果的影响。图2～4直观的表现了定位结果的差异。
图2展示的是在无误差情况下，观测点的数量和分布对于定位结果的影响。以图4可以看出，在无误差情况下，所有的结果都很接近虚拟的震源位置，尤其是震中的误差，被控制在数米之内；而深度的误差也不超过05 km。观测点的数量和分布都会对震中的反演造成影响，不同的是，对于深度的确定，观测点的分布有着明显的影响，但观测点数量的影响却几乎可以忽略。综合起来看，观测点数量越多，分布的越密集，定位的效果就越好；观测点的分布对于定位结果的影响比观测点数量更加显著。需要注意的是，图2b中，随着观测点的分布越来越稀疏，定位结果中深度的概率峰值很明显的变小。形象的说，观测点的分布越稀疏，得到的震源深度越浅。同样，图2d显示，随着观测点数量的增加，概率密度峰值有轻微的右移，即深度的定位结果轻微的增大。
图3展示的是不同的误差对于定位结果的影响。在引入误差后，反演结果的概率密度曲线与图2相比，明显变得平缓，定位结果的误差水平也显著地增加。但是概率密度峰值对应的震中位置和深度并没有偏离虚拟震源的位置。绝大多数的震中位置的误差在50 m以内，深度的误差水平较高，但是也集中在2 km以内。综合起来看，Pg波偏振方向的误差和观测点坐标的误差对于定位结果有很明显的影响，并且都满足误差越小，定位越准确的规律。相比较而言，观测点坐标误差的影响要更大一些，但并不显著。同时，随着观测点坐标误差的增大，深度概率密度峰值也稍稍左移。并且，坐标误差的增大会导致定位深度的减小。在反演结果中，出现了极少数远远偏离虚拟震源位置的异常结果（主要是深度值异常，图中未显示），在后面的章节中，将对这种情况进行讨论。从反演结果的误差水平上来看，各种误差对于定位结果的影响要远远大于观测点的分布和数量对结果的影响。
图4显示的是在一定误差水平下（Pg波偏振方向的正弦值误差≤ 1%，观测点的坐标误差≤ 1%）的定位结果。同图2相比，定位的精度极大的降低了，同时深度的概率密度峰值随观测点分布和数量的变化趋势都被一定程度的放大了。只是因为误差的影响太大，把这个变化给掩盖住了。同样的，反演结果中也出现了极少的遠远偏离虚拟震源的异常点。
4震例检测
本文选取了2014年8月12日05：[KG-*2]55：[KG-*2]2495发生在云南省大理州宾川县的28级地震对本方法进行检测。在震中20 km范围内，共有15个地震记录仪获取了本次地震的记录。定位的结果见表2。
通过表2可看出，本方法的定位结果与云南地震台网的定位结果很接近，震中位置相差小于2 km，震源深度相差也小于2 km。
5讨论与结论
通过对射线追踪方法的测试，我们认为该方法确实可以确定震源的空间位置。但是，由于需要获得比较“纯净的”Pg波偏振信息，因此，该方法不适用于远震定位。
该方法测试时，发现误差对于定位精确度的影响非常大，远远超过了观测点的分布和数量对于定位精确度的影响。总的来说，降低误差水平，提高观测点的数量和密度都能够提升定位结果的精确度；反之亦然。另外，观测点的密度降低和坐标误差的增大会使深度的定位结果变小；而观测点数量的增加会导致深度的结果变大。
和深度不同，震中的最大概然值并没有随着其它条件的变化而变化，造成两者间区别的原因可能是观测点在不同方向上的分布不同。由于出射角的生成是随机的，因此，观测点在水平方向上的分布应该是比较均匀的，不会产生大部分的观测点都位于震中某一边的情况；而在竖直方向上，所有观测点的值都为0，全部位于震源的上方。在有误差存在的反演结果中，出现了极少数的远远偏离震源的异常值。对于这些异常的产生，我们推测有2种可能的原因：随机误差的不均匀；反演程序自身的缺陷。对于震中距不同的观测点来说，其在反演过程中的权重是不一样的。这一点可以通过图2a、b和图4a、b看出。因此当某个权重较大的观测点获得了较大的误差后，就会对定位过程产生严重的影响。另一方面，如果收敛条件的设定不合适，也可能会导致反演程序给出局部极小值点。
从测试中可以看出，震中定位的精确度是很高的。因此，将射线追踪法（或者其它定位方法）得到的震中作为先验条件，再对震源深度进行单独定位，也许会得到更加可靠的结果。另外，我们注意到，在反演的过程中，我们将地表和震源中间的介质分层，并得到了相应的波速比值r（r1，r2，…，rn）。利用这组数值，我们可以通过线性反演获得“发震时刻”和各层介质的“波速”。但是，由于我们并不清楚这组波速比值是否对应着真实的速度结构，因此，我们也不知道通过该方法得到的“发震时刻”和“波速”是否是正确的。
通过一例真实地震的检验，发现本定位方法的结果和其他研究者的结果较为接近，因此，我们认为射线追踪方法可以用于真实地震的定位。
今后，我们将从5个方面展开工作：改进程序和算法，降低异常结果的产生；提高定位结果，尤其是深度定位结果的精确度；发震时刻的定位；复杂介质条件下，射线追踪法的应用；用更多真实的地震数据来检测该方法，并最终将其应用于实际的地震定位中。
[HTK]在本文的写作和修改过程中，杨润海高工提供了地震资料和指导，颜其中研究员对本文提出了建议，我的硕士导师倪四道教授和崔建文研究员也给予了指导和帮助，在此一并致谢。同时，也感谢审稿专家的修改意见和编辑的辛勤工作。
参考文献：
李聪，雷建设2014滇西南地区地壳速度结构的区域波形反演[J].科学通报，59（34）：3398-3415
罗艳，曾祥芳，倪四道2013震源深度测定方法研究进展[J].地球物理学进展，28（5）：2309-2321
田玥，陈晓非2002地震定位研究综述[J].地球物理学进展，17（1）：147-155
朱元清，赵仲和1997提高地震定位精度新方法的研究[J].地震地磁观测与研究，18（5）：59-67
CROSSON R S1976Crustal structure modeling of earthquake data，1，Simultaneous least squares estimation of hypocenter and velocity parameters[J].J Geophys Res，81（17）：3036-3046
DOUGLAS A1976Joint epicenter determination[J].Nature，215：45-48
ENGDAHL E R，HILST R V D，RAYMOND B1998Global teleseismic earthquake relocation with improved travel times and procedures for depth determination[J].Bull Seism Soc Am，88（3）：722-743
GEIGER L1912Probability method for the determination of earthquake epicenters from arrival time only[J].Bull St Louis Univ，（8）：60-71
KLEIN F W1978Hypocenter location program HYPOINVERSE Part I：Users guide to versions 1，2，3 and 4[R].U S Geol Sur Open-File Rept，78-694
LEE W H K，LAHR J C1975HYPO71：A computer program for determining hypocenter，magnitude，and first motion pattern of local earthquake[R].U S Geol Surv Open-File Rept，75-311
LIENERT B R，BERG E，FRAZER L N1986Hypocenter：An earthquake location method using centered，scaled，and adaptively damped last squares[J].Bull Seism Soc Am，76（3）：771-783
LOMNITZE C1977A fast epicenter location program[J].Bull Seism Soc Am，67（2）：425-431
NELSON G D， VIDALE J E1990Earthquake locations by 3-D finite-difference travel times[J].Bull Seism Soc Am，80（2）：395-410
PRUGGER A F，GENDZWILL D J1988Microearthquake location：A nonlinear approach that makes use of a simplex stepping procedure[J].Bull Seism Soc Am，78（2）：799-815
ROMNEY C1957Seismic waves from the Dixie Valley Fairview Peak earthquakes[J].Bull Seism Soc Am，47（4）：301-319
SPENCE W1980Relative epicenter determination using P-wave arrival-time differences[J].Bull Seism Soc Am，70（1）：171-183
TARANTOLA A，VALETTE B1982Inverse problem = quest for information[J].J Geophys R，50：159-170
THURBER C H1985Nonlinear earthquake location：theory and examples[J].Bull Seism Soc Am，75（3）：779-790
WALDHAUSER F，ELLSWORTH W L2000A double-difference earthquake location algorithm：method and application to the Northern Hayward Fault，California[J].Bull Seism Soc Am，90（6）：1353-1368
WARREN L，SHEARER P2005Using the effects of depth phases on P-wave spectra to determine earthquake depths[J].Bull Seism Soc Am，95（1）：173-184
ZHAO L S，HELMBERGER D V1994Source estimation from broadband regional seismograms[J].Bull Seism Soc Am，84（1）：91-104
ZHU L P，HELMBERGER D V1996Advancement in source estimation techniques using broadband regional seismograms[J].Bull Seism Soc Am，86（5）：1634-1641[ZK）][HJ][FL）][LM]
[STHZ][WT4HZ][JZ]Near Earthquake Locating Based on Ray Tracing Method
[WT5B1][STBZ][JZ]ZHANG Qian1，2，YIN Geng2，WANG Yushi3，LIN Guoliang2
[WT5"B1X][JZ]（1 University of Science and Technology of China，Hefei 230026，Anhui，China）[JZ]
（2 Earthquake Administration of Yunnan Province，Kunming 650224，Yunnan，China）[JZ]
（3 Institute of Geophysics，CEA，Beijing 100081，China）
[WT5HZ][JZ]Abstract[WTB1]
We demonstrate a ray tracing method for earthquake locating without prior velocity structure，which is based on ray theory As the seismic ray focus on the hypocenter，first we divide the layer into sub-layers，and compute the velocities of the sub-layers and the location of the hypocenter by grid-search；then we increase the number of the sub-layers step by step，and repeat the process before correcting the hypocenter location；at last，we get a more accurate result Simulative data have been used to test the method，the results show that，this method can be used for earthquake location in theory We use this method to locate a natural earthquake In contrast with the location result of Yunnan Seismic Network，we find the two results are close
摘要：介绍了一种基于地震波射线理论、不需要先验速度结构的定位震源空间位置的方法——射线追踪法。该方法利用地震射线会反向汇聚到震源的几何性质，对地面介质分层，并用网格搜索反演方法计算分层速度参数，进而确定震源位置；然后逐步增加分层，重复前面的步骤，对震源位置进行校正，最终得到一个较为精确的定位结果。通过虚拟事件对该方法进行了测试，测试结果表明，该方法理论上能够用于近震空间位置的定位；使用该方法对一个真实震侧进行定位，并和云南地震台网的定位的结果进行了比较，两者较为接近。
关键词：地震定位；无速度结构；射线理论；P波偏振
中图分类号：P3153文献标识码：A文章编号：1000-0666（2017）02-0203-08
0引言
震源空间位置的确定是地震学的经典问题——地震定位的一部分，也是地震学研究的基础。地震定位，无论是对于地球物理学研究，如地球内部的环境构造、地震的发生机制，还是对于震后救援工作的需求，如震害范围和程度的预估、地震趋势的预测以及地震预警，都有着无法替代的重要性。因此，地震定位的方法以及提高地震定位精度和速度的方法，一直都是地震学家的研究目标。
早期的地震定位主要是通过几何作图完成的。Geiger（1912）根据各个观测点到时差最小化的原则，将地震波走时方程组线性化，然后通过最小二乘法来定位地震。这也为以后使用计算机来定位地震提供了基础。20世纪70年代后，伴随着现代计算机技术的发展，地震学家们基于Geiger的理论，给出了一系列地震地位的程序和算法（Lee，Lahr，1975；Klein，1978；Lienert et al，1986；Nelson，Vidale，1990），并根据实际应用中出现的问题，提出了相应的改进方法（Lienert et al，1986；Prugger，Gendzwill，1988）。随后，多事件定位方法（Douglas，1976；Crosson，1976；Spence，1980；Waldhauser，Ellsworth，2000）、基于距离残差的空间域定位方法（Romney，1957；Lomnitze，1977）和非线性的定位方法（Prugger，Gendzwill，1988；Thurber，1985；Tarantola，Valette，1982）也相继出现。震相、振幅，甚至是全波形的信息都被用于震源位置的确定（Engdahl，et al，1998；Warren，Shearer，2005；罗艳等，2013；Zhao，Helmberger，1994；Zhu，Helmberger，1996）。
在地震定位的过程中，影响结果精度的因素大致有：观测系统的结构、速度结构的选取或假设、震相识别以及观测数据（包含观测点的空间位置）的读取和计算过程中产生的误差（朱元清，赵仲和，1997；田玥，陈晓非，2002）。绝大多数时间域和空间域的定位方法，都需要使用某一组先验的速度结构。因此，该种方法中速度结构的选取，将直接影响到结果的误差水平。而通过波形信息来确定震源位置时，震相识别、震源机制解和地球结构模型的可靠性，将决定定位结果的可靠性。
本文将介绍一种基于地震波射线理论，且不需要先验速度结构的定位近震震源位置的方法——射线追踪法，并使用虚拟数据和一个真实震例对该方法进行了检验。[HJ2mm]
1基本原理
假设地球介质是水平均匀且各向同性的分层介质。根据惠更斯原理和地震波射线理论，一束地震波在传播时，射线和分层界面法线的夹角满足方程[KH*2]
[SX（]v1sinθ1[SX）]=……=[SX（]vjsinθj[SX）]=……=[SX（]vnsinθn[SX）]=const[JY]（1）[KH*3/4D]
其中：vj和θj是射线在第j层的波速和夹角（震源所在位置为第n层，地表为第一层上边界）。将地表处的夹角和波速作为已知条件，则[KH*2]
θj=arcsin（[SX（]vjv1[SX）]sinθ1）[JY]（2）[KH*3/4D]
从（2）式可以看出，射線与法线的夹角只与地表入射角θn和当前波速与地表波速的比值vj/v1（记为rj）有关。另一方面，对于一束从震源发出的射线i，其在第j层介质中的传播，符合方程[KH*1]
sij=Δhj/tanθij[JY]（3）[KH*1D]
式中：sij和θij分别代表射线i在第j层介质中的水平位移和与法线的夹角；Δhj表示第j层的厚度。同时，对于射线i，满足方程[KH*1]
θij=arcsin（rjsinθi1）[JY]（4）[KH*1]
结合式（3），可得：[KH*1]
sij=Δhj/tan（arcsin（rjsinθi1））[JY]（5）[KH*1]
因此，震源到射线接收点（即观测点）之间的水平位移为：[KH*2]
[JP2]Si[KG-*4]=[KG-*2]∑[DD（]nj=1[DD）]sij[KG-*4]=[KG-*2]∑[DD（]nj=1[DD）]Δhj/tan（arcsin（rjsinθi1））[JY]（6a）[KH*2D][JP]
或
[JP3]Si[KG-*4]=[KG-*4]∫[DD（]H[KG*4]i0[DD）]si（h）dh[KG-*4]=[KG-*4]∫[DD（]H[KG*4]i0[DD）]1/tan（arcsin（r（h）sinθi1））dh[JP][KH-*4][JX*1][JY]（6b）[KH*1D]
其中：Hi代表震源深度。θi1和Pg波的偏振方向相关，可以通过观测得到；而Si和Hi代表着观测点和震源空间位置的矢量差。因此，我们可以得到关于震源空间位置F（x，y，h）和r（r1，r2，…rn）的方程（组）：[KH*1]
[JB（{][KF（]（xi-y）2+（yi-y）2[KF）]-Si=0[KG*2]（hi-h）-H=0[JB）][JY]（7）[KH*1D]
其中：（xi，yi，hi）为观测点i的空间位置。本文所介绍的方法，就是要找到一组（F，r），使（7）式成立，或者使等式兩边的差值尽量小，从而实现震源的定位。
[KG（0.15mm]由于式（6）和（7）含有三角函数，不能采用线性方程组的方法来求解，而必须寻找一种合适的反演方法来确定（F，r）。网格搜索法是一种不错的方法，但是其计算量会随着参数的增加而呈指数增长。如果待求解的参数越多，则计算所花时间越长。因此，我们采用了一种“分步”的网格搜索法来定位震源，待求参数变少，计算效率得到了提升。图1是对这种方法的一个简单的展示，由于看起来像追溯射线的路径回到震源，我们称该方法为“射线追踪法”，其操作步骤如下：[KG）]
（1）根据各个观测点接收到的Pg波的偏振方向，反向生成一组射线；寻找到各条射线距离最短的一点F1，并将其作为初始定位的结果（F1，r1）（r1=1）。
（2）将介质平均分成2层，每层的厚度为h1/2，r1取步骤（1）的结果。使用网格搜索法，寻找第二层的波速比r2，以及相应的点F2，使F2到根据（r1，r2）做出的反向射线的距离和最小。将（F2，r2）作为最新的结果；
（3）重复步骤（2），直到最后得到的结果（Fn，rn）满足一定的收敛条件。Fn即为最终的定位结果。
通过射线追踪法，我们一步步将介质的分层数从2层增加至n层，且每一次分层都只用搜索新增层的波速比，从而大大减少了计算量。最终，我们会得到震源的位置Fn以及1组和各层相对应的波速比（r1，r2，…，rn）。
2测试数据
为了验证此地震定位方法的准确性，本文虚拟了多组数据对该方法进行测试。测试数据生成[FL）][SD1，1]
[JZ][XC张潜1.TIF][KH*2/3D][HT5K][JZ]
图1射线追踪法步骤的简单演示[JZ]
Fig1Demonstration of ray tracing method[KH*1]
的方法如下：
（1）选择1个一维速度模型作为基准模型，假定震中位置为0，震源深度为10 km。
（2）随机生成一组出射角为（θ1n，…，θmn）的射线，根据步骤1的速度模型和震源位置，计算各条射线到达地表（深度为0）处的角度（θ11，…，θm1）和震中距（S1，…，Sm）。
（3）为（S1，…，Sm）和（θ11，…，θm1）引入误差，组成一组数据，分别对应m个虚拟观测点的坐标和接收到射线的Pg波偏振方向。
对于震中一定范围以外的观测点，Pg波的到达往往会晚于Pn波。此时，Pg波的振幅信息会和其它震相的信息混杂在一起，这就对Pg波偏振方向的确定造成了困难。因此，在生成测试数据时，要注意将θ值限定在一定的范围内。在本研究中，θ的取值范围分为3档，分别为≤ 30°、≤ 45°和≤ 60°，本文将对结果进行比较分析，以观察不同的出射角范围对定位精度的影响。同理，为了解观测误差和观测点数量对于定位精度的影响，本文也选取了相应的参数来进行测试。为了计算上的方便，我们将三维空间简化为二维，即去掉了一个水平维度，变成了如图1所示的情况。这样的简化并不会对测试造成本质性的影响。
本文共虚拟了14万组数据来进行测试，其详细信息见表1。所有的数据按照不同的观测点数量、出射角度范围、误差范围等因素分入14个大组，每个大组包含1万组数据。除观测点数量以外，生成每一组数据所用到的角度和误差都是在限定范围内随机生成的。显然，当误差水平为0，则意味着该组数据没有引入相应的误差。
3测试结果及比较
我们虚构了一个震中位置为0，震源深度为10 km的震源，并根据随机生成的参数生成了观测数据，最终反演得到了14万组定位结果。对这些结果进行比较，以观察不同的条件对于反演结果的影响。图2～4直观的表现了定位结果的差异。
图2展示的是在无误差情况下，观测点的数量和分布对于定位结果的影响。以图4可以看出，在无误差情况下，所有的结果都很接近虚拟的震源位置，尤其是震中的误差，被控制在数米之内；而深度的误差也不超过05 km。观测点的数量和分布都会对震中的反演造成影响，不同的是，对于深度的确定，观测点的分布有着明显的影响，但观测点数量的影响却几乎可以忽略。综合起来看，观测点数量越多，分布的越密集，定位的效果就越好；观测点的分布对于定位结果的影响比观测点数量更加显著。需要注意的是，图2b中，随着观测点的分布越来越稀疏，定位结果中深度的概率峰值很明显的变小。形象的说，观测点的分布越稀疏，得到的震源深度越浅。同样，图2d显示，随着观测点数量的增加，概率密度峰值有轻微的右移，即深度的定位结果轻微的增大。
图3展示的是不同的误差对于定位结果的影响。在引入误差后，反演结果的概率密度曲线与图2相比，明显变得平缓，定位结果的误差水平也显著地增加。但是概率密度峰值对应的震中位置和深度并没有偏离虚拟震源的位置。绝大多数的震中位置的误差在50 m以内，深度的误差水平较高，但是也集中在2 km以内。综合起来看，Pg波偏振方向的误差和观测点坐标的误差对于定位结果有很明显的影响，并且都满足误差越小，定位越准确的规律。相比较而言，观测点坐标误差的影响要更大一些，但并不显著。同时，随着观测点坐标误差的增大，深度概率密度峰值也稍稍左移。并且，坐标误差的增大会导致定位深度的减小。在反演结果中，出现了极少数远远偏离虚拟震源位置的异常结果（主要是深度值异常，图中未显示），在后面的章节中，将对这种情况进行讨论。从反演结果的误差水平上来看，各种误差对于定位结果的影响要远远大于观测点的分布和数量对结果的影响。
图4显示的是在一定误差水平下（Pg波偏振方向的正弦值误差≤ 1%，观测点的坐标误差≤ 1%）的定位结果。同图2相比，定位的精度极大的降低了，同时深度的概率密度峰值随观测点分布和数量的变化趋势都被一定程度的放大了。只是因为误差的影响太大，把这个变化给掩盖住了。同样的，反演结果中也出现了极少的遠远偏离虚拟震源的异常点。
4震例检测
本文选取了2014年8月12日05：[KG-*2]55：[KG-*2]2495发生在云南省大理州宾川县的28级地震对本方法进行检测。在震中20 km范围内，共有15个地震记录仪获取了本次地震的记录。定位的结果见表2。
通过表2可看出，本方法的定位结果与云南地震台网的定位结果很接近，震中位置相差小于2 km，震源深度相差也小于2 km。
5讨论与结论
通过对射线追踪方法的测试，我们认为该方法确实可以确定震源的空间位置。但是，由于需要获得比较“纯净的”Pg波偏振信息，因此，该方法不适用于远震定位。
该方法测试时，发现误差对于定位精确度的影响非常大，远远超过了观测点的分布和数量对于定位精确度的影响。总的来说，降低误差水平，提高观测点的数量和密度都能够提升定位结果的精确度；反之亦然。另外，观测点的密度降低和坐标误差的增大会使深度的定位结果变小；而观测点数量的增加会导致深度的结果变大。
和深度不同，震中的最大概然值并没有随着其它条件的变化而变化，造成两者间区别的原因可能是观测点在不同方向上的分布不同。由于出射角的生成是随机的，因此，观测点在水平方向上的分布应该是比较均匀的，不会产生大部分的观测点都位于震中某一边的情况；而在竖直方向上，所有观测点的值都为0，全部位于震源的上方。在有误差存在的反演结果中，出现了极少数的远远偏离震源的异常值。对于这些异常的产生，我们推测有2种可能的原因：随机误差的不均匀；反演程序自身的缺陷。对于震中距不同的观测点来说，其在反演过程中的权重是不一样的。这一点可以通过图2a、b和图4a、b看出。因此当某个权重较大的观测点获得了较大的误差后，就会对定位过程产生严重的影响。另一方面，如果收敛条件的设定不合适，也可能会导致反演程序给出局部极小值点。
从测试中可以看出，震中定位的精确度是很高的。因此，将射线追踪法（或者其它定位方法）得到的震中作为先验条件，再对震源深度进行单独定位，也许会得到更加可靠的结果。另外，我们注意到，在反演的过程中，我们将地表和震源中间的介质分层，并得到了相应的波速比值r（r1，r2，…，rn）。利用这组数值，我们可以通过线性反演获得“发震时刻”和各层介质的“波速”。但是，由于我们并不清楚这组波速比值是否对应着真实的速度结构，因此，我们也不知道通过该方法得到的“发震时刻”和“波速”是否是正确的。
通过一例真实地震的检验，发现本定位方法的结果和其他研究者的结果较为接近，因此，我们认为射线追踪方法可以用于真实地震的定位。
今后，我们将从5个方面展开工作：改进程序和算法，降低异常结果的产生；提高定位结果，尤其是深度定位结果的精确度；发震时刻的定位；复杂介质条件下，射线追踪法的应用；用更多真实的地震数据来检测该方法，并最终将其应用于实际的地震定位中。
[HTK]在本文的写作和修改过程中，杨润海高工提供了地震资料和指导，颜其中研究员对本文提出了建议，我的硕士导师倪四道教授和崔建文研究员也给予了指导和帮助，在此一并致谢。同时，也感谢审稿专家的修改意见和编辑的辛勤工作。
参考文献：
李聪，雷建设2014滇西南地区地壳速度结构的区域波形反演[J].科学通报，59（34）：3398-3415
罗艳，曾祥芳，倪四道2013震源深度测定方法研究进展[J].地球物理学进展，28（5）：2309-2321
田玥，陈晓非2002地震定位研究综述[J].地球物理学进展，17（1）：147-155
朱元清，赵仲和1997提高地震定位精度新方法的研究[J].地震地磁观测与研究，18（5）：59-67
CROSSON R S1976Crustal structure modeling of earthquake data，1，Simultaneous least squares estimation of hypocenter and velocity parameters[J].J Geophys Res，81（17）：3036-3046
DOUGLAS A1976Joint epicenter determination[J].Nature，215：45-48
ENGDAHL E R，HILST R V D，RAYMOND B1998Global teleseismic earthquake relocation with improved travel times and procedures for depth determination[J].Bull Seism Soc Am，88（3）：722-743
GEIGER L1912Probability method for the determination of earthquake epicenters from arrival time only[J].Bull St Louis Univ，（8）：60-71
KLEIN F W1978Hypocenter location program HYPOINVERSE Part I：Users guide to versions 1，2，3 and 4[R].U S Geol Sur Open-File Rept，78-694
LEE W H K，LAHR J C1975HYPO71：A computer program for determining hypocenter，magnitude，and first motion pattern of local earthquake[R].U S Geol Surv Open-File Rept，75-311
LIENERT B R，BERG E，FRAZER L N1986Hypocenter：An earthquake location method using centered，scaled，and adaptively damped last squares[J].Bull Seism Soc Am，76（3）：771-783
LOMNITZE C1977A fast epicenter location program[J].Bull Seism Soc Am，67（2）：425-431
NELSON G D， VIDALE J E1990Earthquake locations by 3-D finite-difference travel times[J].Bull Seism Soc Am，80（2）：395-410
PRUGGER A F，GENDZWILL D J1988Microearthquake location：A nonlinear approach that makes use of a simplex stepping procedure[J].Bull Seism Soc Am，78（2）：799-815
ROMNEY C1957Seismic waves from the Dixie Valley Fairview Peak earthquakes[J].Bull Seism Soc Am，47（4）：301-319
SPENCE W1980Relative epicenter determination using P-wave arrival-time differences[J].Bull Seism Soc Am，70（1）：171-183
TARANTOLA A，VALETTE B1982Inverse problem = quest for information[J].J Geophys R，50：159-170
THURBER C H1985Nonlinear earthquake location：theory and examples[J].Bull Seism Soc Am，75（3）：779-790
WALDHAUSER F，ELLSWORTH W L2000A double-difference earthquake location algorithm：method and application to the Northern Hayward Fault，California[J].Bull Seism Soc Am，90（6）：1353-1368
WARREN L，SHEARER P2005Using the effects of depth phases on P-wave spectra to determine earthquake depths[J].Bull Seism Soc Am，95（1）：173-184
ZHAO L S，HELMBERGER D V1994Source estimation from broadband regional seismograms[J].Bull Seism Soc Am，84（1）：91-104
ZHU L P，HELMBERGER D V1996Advancement in source estimation techniques using broadband regional seismograms[J].Bull Seism Soc Am，86（5）：1634-1641[ZK）][HJ][FL）][LM]
[STHZ][WT4HZ][JZ]Near Earthquake Locating Based on Ray Tracing Method
[WT5B1][STBZ][JZ]ZHANG Qian1，2，YIN Geng2，WANG Yushi3，LIN Guoliang2
[WT5"B1X][JZ]（1 University of Science and Technology of China，Hefei 230026，Anhui，China）[JZ]
（2 Earthquake Administration of Yunnan Province，Kunming 650224，Yunnan，China）[JZ]
（3 Institute of Geophysics，CEA，Beijing 100081，China）
[WT5HZ][JZ]Abstract[WTB1]
We demonstrate a ray tracing method for earthquake locating without prior velocity structure，which is based on ray theory As the seismic ray focus on the hypocenter，first we divide the layer into sub-layers，and compute the velocities of the sub-layers and the location of the hypocenter by grid-search；then we increase the number of the sub-layers step by step，and repeat the process before correcting the hypocenter location；at last，we get a more accurate result Simulative data have been used to test the method，the results show that，this method can be used for earthquake location in theory We use this method to locate a natural earthquake In contrast with the location result of Yunnan Seismic Network，we find the two results are close
