2017, Vol. 60
2. 中国科学院大学, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
近场震源定位在实际应用中需求广泛，如爆破监测、车辆追踪等.传统的震源定位方法都是以地震波信号为基础，通过求解目标函数的极小值.不同的定位方法区别在对于目标函数的构造、处理以及优化方法上(田玥和陈晓非, 2002).面对复杂的波速环境，国内外许多学者提出了不同的定位算法，例如：联合反演法(Crosson, 1976)，主事件定位法(Spence, 1980)，EHB法(Engdahl et al., 1998)，单纯形法(赵珠等, 1994)，SAMS方法(高星等, 2002)，双差定位法(Waldhauser and Ellsworth, 2000)，反投影法(Walker et al., 2005)，DE算法(李志伟等，2006)，Bayesian法(MATSU, apos, Ura, 2009)，逆时偏移法(Nakata and Beroza, 2016)、基于自适应量子遗传算法(范建柯等, 2016)等等.这些方法在反演过程中需要一个相对可靠的速度模型.但是在近场地表震源定位中，获得高质量的近地表二维速度模型的成本较高，且反演计算耗时相对较大.因此为有效解决近场震源定位中的精度与计算速度的问题，需要一种新的方法.
对于非稳态信号的定位，一个最基本的几何约束思路是通过计算地震仪对之间信号的到时差，震源则位于以这一对地震仪为离心点的双曲线的一支上，所有的双曲线相交的点或区域则为定位结果.到时差可以通过波形互相关(Carter, 1981)或者自适应滤波器(Doh-Hyoung and Youngjin, 2005; Ho et al., 1993; Klee et al., 2006)计算得到.该方法要求不同地震仪接收的信号必须具有很高的互相关性.另一种方法是通过拾取初至波的绝对到时，然后相减获取到时差.常见的自动、高质量拾取到时的方法包括AIC法(Maeda, 1985)，PAI-S/K法(Saragiotis et al., 2002)，峰度拾取算法(刘劲松等, 2013)，短时信号小波分析法(Thanasopoulos and Avaritsiotis, 2011)以及连续小波分解法(Bogiatzis and Ishii, 2015).
这种几何约束的框架均暗含了一个假设：介质的波速在各个方向上都是均匀，并且波在接收点与震源点之间的传播没有阻碍.实际工作中这种条件是几乎无法满足.因此，该假设条件会给定位测量带来误差.
为解决由介质的非均匀性及不连续性带来的震源定位难题，Baxter等(2007)提出了到时差分布图(Delta T Mapping，后称DTM)方法，并率先在具有复杂几何结构的金属材料板的声波定位中取得了较高的定位精度.该方法利用人工震源获得研究区域的到时差“结构”，这些“结构”图可进一步用来进行新震源的高精度定位.Eaton等(2012)证明在具有各向异性的材料上，DTM也能得到非常好的结果.在数据处理手段方面，Al-Jumaili等(2016)在前人研究的基础上，提出全自动的处理技术，一方面减少了由主观分析在挑选高质量到时数据时引入的误差，另一方面大大提高了数据处理效率.Hensman等(2010)不同于他人利用线性插值获得位置到时差的映射，而是利用高斯过程回归(Gaussian Process Regression，后简称GPR)得到震源定位，同样取得了较高的定位精度.
DTM方法被证明在小尺度各向异性、复杂几何结构的材料板声波定位中具有较高精度，且无需波速度信息，本文首次将其引入近场震源定位中.本文首先介绍了DTM方法实现的具体步骤与两种不同定位方法.然后介绍了本次实际定位实验的具体情况，包括点位布设方案、地震仪器、仪器布设方案与炸点实验过程.随后分析了实验测得的近场炸点地震数据，说明了走时拾取的方法和有效数据的挑选.最后，通过上述方法得到了定位结果，分析了定位误差，比较了不同方法在应用上的优缺点，并提出了实验的改进方案，最后结合震源扫描法，DTM方向具有对多目标定位问题的潜力.结果表明，DTM方法是一种可靠、高精度、具有近实时性的近场震源定位方法.
2 定位方法DTM方法定位的根本原理是，首先利用人工炸点实验对区域的“到时差结构”进行刻画，然后对于新的爆破事件，利用它的到时差数据代入之上的结构中，寻求最匹配的一点，即为定位结果.所谓到时差模型，是指一种对区域内到时差与位置之间的关系的描述.
2.1 DTM方法定位流程DTM方法用于近场震源包括到时差模型训练与事件定位两个主要环节，其中，到时差模型的获取包括震源与地震仪点位设计、人工震源实施、到时差模型计算三个环节(Baxter et al., 2007).
(1) 震源与地震仪点位设计：首先选定测区，并进行网格划分，震源点通常设计位于网格结点处，将用来得到相应格点处的到时差结构，称为模型点位.地震仪点位设计在模型节点的四周.
(2) 人工震源实施：在网格点位置进行人工震源爆破.震源能量要求在距源点最远处的仪器能接受到有效的爆破信号，以获得更多有效的地震记录.
(3) 到时差模型计算：获取震源的波形数据后，对于每一个震源事件，分别计算该震源事件的每一个地震仪对之间的到时差.假定有4个地震仪(S01—S04)，共获得6个到时差数据，分布是S01—S02、S01—S03、S01—S04、S02—S03、S02—S04、S03—S04.
(4) 事件定位：对于新的震源事件，计算每一个地震仪对的到时差，代入已建立的到时差模型(见下一节)，进行新的事件的定位.对于多个事件，若模型训练时的台阵位置不改变，该模型可重复使用.
值得注意的是，地震仪点位之所以分布在模型节点的四周，是为了获得更大范围的地震仪对方位角差，来更好地约束定位结果.仪器布设的偏移距一般应该保证两个条件：一是偏移距不能太大，以致仪器不能接受到有效的爆破信号；二是偏移距不能太小，绝大多数地震仪对之间的到时差应可测，而不是湮没于采样间隔之间而无法区分.
对于到时差的计算，一是通过拾取绝对到时，即Δt=ta-tb，ta、tb分别为同一爆破信号在不同记录中的到时；二是通过波形互相关，来提取由同一震源产生的相关的爆破信号进而计算获得，即Δt=
在第三步计算DTM后，获得了每个模型点上、对应每一地震仪对的到时差.这些模型数据包含了该测区的“到时差”结构，是该测区炸点定位的依据.为了达到高精度的定位要求，需进一步利用已有的模型数据构建更好的定位模型，该过程中有两种方法.
(1) 网格搜索法
对于每一个地震仪对，每个模型点位上已经分别获得了一个到时差的值，即每个模型点位置(X, Y)都对应一个到时差ΔT.为了获得除模型点位外、区域内其他位置的到时差的分布，可以通过线性散点插值获得更高密度的分布图，后续定位坐标值有效位数不小于插值间隔.对新震源事件定位，采用最小差法(Al-Jumaili et al., 2016)，计算下式：
|
(1) |
其中Ts是新震源事件相对某一地震仪对的到时差，Tmap是同一地震仪对的到时差分布图.绝对值表示在该地震仪对下到时差的残差，残差最小的区域即为新震源事件的所在位置；sum表示所有地震仪对的累积到时差，n代表地震仪对的个数.对全部有效模型进行线性(平权)叠加，累积到时差最小的位置作为新震源事件的定位结果，该位置实际为全部有效到时差模型残差最小区域的交集.
(2) 统计定位法
与上述方法获得位置到时差，即(X, Y)→ΔT不同是，这里可以利用高斯回归过程(Gaussian Process Regression，简称GPR)(Rasmussen, 2006)获得到时差与位置坐标的映射关系，通过计算新的到时差与原有的到时差的相似程度，结合训练样本的坐标来预测相对应的新位置坐标.
GPR是为了通过训练样本数据，估计在新样本空间上的后验概率分布函数(Rasmussen, 2006).考虑一组训练样本D={(x1, y1), (x2, y2), …, (xn, yn)}，它由含噪声的某一过程f实现：
|
(2) |
其中xi是
它对样本的要求是，不同样本点上的函数值是互相关的，函数值之间的协方差，例如f(xp)和f(xq)之间是由xp和xq所决定.这种决定关系称为核函数，记为k(xp, xq).最广泛使用的核函数为高斯函数：
|
(3) |
其中，σ2f是输入样本的方差.l为核函数尺度，表示两个样本点间互相关随之距离的衰减变化.这两个参数控制了GPR的光滑程度.
实际实验中，观测值并不是真正的函数值，而是带有一定的噪声(假设为高斯噪声).因此相应的协方差函数应为
|
(4) |
其中σ2n为观测得到的高斯噪声方差，δpq=1当且仅当p=q，其余情况取0.对于完整的输入样本X，对应的目标值 y的协方差为
|
(5) |
其中K是n×n的协方差矩阵，K[p, q] =k(xp, xq).协方差与观测值y之间相互独立.
在样本空间中，y ~
|
(6) |
其中μx*=kT*(K +σn2 I)-1 y，σ2x*=k(x*, x*)-kT*(K +σ2n I)-1k*.
训练GPR涉及到3个超参数的选择：σ2f，l以及σ2n.超参数的选择对模型表达有着极大的影响，合适的超参数才能提供高的定位精度(Ferris, Haehnel, Fox, 2006).超参数的优化，一是可以利用广义交叉验证(Hensman et al., 2010; Nelder Mead, 1965)，二是可以寻求目标值的最大似然(Rasmussen, 2006).
3 实验设置实验场所位于北京市顺义区某农场中，西面紧靠公路、工厂、学校.震源定位测试的区域主要包括荆棘地、道路、树林、蔬菜种植地，树林部分区域布设了地下管网.地震波的主要传播路径在地表的土壤介质中.
本实验采用QS-05A便携式微功耗宽频带地震仪(图 2a)，频带宽度为5 s~250 Hz，系统噪音低于1 μV，内置GPS实时定位授时，采样率500 Hz，数据记录格式为MiniSeed.本次实验一共布置了22台地震仪(图 1a)，布设采取填埋方式(图 2b)，仪器布放位置由手机定位导航，与实际有一定偏差.由于用于模型训练的数据与测试事件的数据是在同一台阵下记录得到，因此地震仪的位置不影响定位分析.另外，本次地震仪均采用连续实时GPS授时系统，保证了单个仪器的时间准确性.
|
图 1 实验点位布设图 (a)台站与模型点位分布图，原点处为M1；(b)测试区域爆破点分布图. Fig. 1 Layout of the experiment setting (a) Layout for stations and model grids, the origin is at M1; (b)Layout for explosions. |
|
图 2 (a) QS-05A便携式微功耗宽频带地震仪；(b)野外布设图 Fig. 2 (a) Seismic sensor of type QS-05A; (b) The scene in the experiment |
爆破点定位采用差分GPS(DGPS)定位，定位精度达厘米级.本次共用24炮(以M标记)作为模型训练使用，10炮(以T标记)作为测试，以检验定位精度.模型网格的横向(东西向)间距为28 m，纵向(南北向)间距为30 m，模型横向长度为140 m，纵向长度为90 m(图 1b).
爆破采用普通爆竹，第一种为单响，威力大(编号B01)；第二种为间隔约2.5 s的两响，威力小(编号B02).对于B02型爆竹，布埋时使其头向下，这样保证在第一炮响后，它仍旧能大约留在原地而不是在空中发生第二次爆破.在24个模型网格中，除7、10、15、18、19号点用的是B02型爆竹，其余均采用B01型.在10个测试事件中，除7号点用B01型外，其余全部采用B02型.
4 数据记录及分析 4.1 地震记录本实验仅采用垂直分量的地震记录.图 3a与图 3b分别对应B01型和B02型爆竹的原始地震记录.可以看出，B01型爆竹的数据信噪比大于B02型.对于偏移距较大的地震仪(S2、S3、S4等)，B02型爆竹的信号几乎湮灭在背景噪声中.为得到信号的有效频段，对原始记录进行了小波分析(图 5).分析结果表明爆破信号优势频带为64~190 Hz，环境噪声的优势频带为2~33 Hz，即有效信号与环境噪声的频带可明显区分开.在后续计算到时差的过程中，仅采用了对应中心频率为64~190 Hz的小波系数分量.
|
图 3 M1记录的B01型爆竹原始地震记录(垂直分量) 横坐标为UTC+08:00时间，纵坐标为台站编号，图中为振幅归一化波形. Fig. 3 The raw record (vertical component) at M1 for B01 explosions The horizontal ordinate denotes time (UTC+08:00) while vertical one the index of seismometers. Waves′ amplitudes are normalized. |
|
图 4 M7记录的B02型爆竹原始地震记录(垂直分量) 横坐标为UTC+08:00时间，纵坐标为台站编号，图中为振幅归一化波形. Fig. 4 The raw record (vertical component) at M7 for B02 explosions The horizontal ordinate denotes time (UTC+08:00) while vertical one the index of seismometers. Inside is the waves with normalized amplitude. |
|
图 5 S1地震仪记录的B01型(a)和B02型(b)爆竹信号连续小波变换结果 Fig. 5 Continuous wavelet transformations for B01 (a) and B02 (b) explosions |
无论是提取精确的走时，或是进行波形互相关，首先要获得正确的时间与波形.本次实验中，仪器内部采用实时GPS进行授时，会实时将采样时间误差记录在仪器内的数据中(MiniSeed头段中的TimeCorrection).除了单个采集单元中的时钟误差外，为了让不同仪器的数据同步在相同的时间轴上，还存在同步误差.这两者之和即为将需要校正的钟差值.一般钟差校正的方法如下：
首先对时间序列进行偏移以校正大于采样间隔的钟差值，然后对时间序列进行等间隔插零，使得原信号的采样点位于插值后的采样点中.接着对信号进行FFT，去掉高于原Nyquist频率的成份，再通过IFFT得到插值后的信号.最后对信号进行等间隔降采样，得到校正后的波形.
在GPS实时授时情况下，仪器本身钟差极小，几乎都为0(仪器钟差最小可分辨值为0.1 ms).须要校正的钟差主要是与统一时间轴间的偏差.图 6为钟差较正的一个例子，从图中可看出其对信号的偏差的影响基本也很小.后续的计算均使用校正后的信号.
|
图 6 钟差较正的例子 信号为M12震源在S7处的记录，蓝色为原始信号，橙色为校正后信号. Fig. 6 An example of time correction The signal is recorded for M12 at S7. The blue denotes original signal while the orange the one after correction. |
AIC方法可以快速、自动拾取地震波到时.在地震波或声波走时拾取中最常见的AIC公式为
|
(7) |
其中var(x)代表时间序列x的方差，x[1:t]代表时间序列x从时间1到t的切片，T为所选取的信号的总时长.
该方法一般需要引入一种信号的特征函数，使得对应的AIC的变化更为敏感，提高拾取精度(Allen, 1982).拾取到时若仅采用宽时窗下的特征信号，将所得AIC的全局最小值处作为到时可能存在偏差(图 7a)，例如S7地震台站对M12震源的提取偏差约为15 ms.AIC方法提取到时非常依赖波形特征.理想情况下，AIC计算的范围应开始于非有效信息部分，终结于有效信息结束处.Sedlak等(2008)提出了双步AIC法，在第一步计算信号特征函数下的AIC(图 7b)，得到一个初步的到时估计点后，第二步集中在这个估计点的一个较小的邻域再重新计算特征函数的AIC(图 7c)，所得的全局最小值处作为改进的到时估计.本文采用的特征函数为CF(i)=x(i)2+R*(x(i)-x(i-1))2.通过数据测试，本文取R为10.第一次AIC计算的范围选为波形最大振幅值前后分别0.32 s和0.18 s.第二次AIC计算的范围选为第一次拾取到的到时前后分别0.12 s和0.08 s.
|
图 7 拾取绝对到时图 (a)归一化、高通滤波后、震源为M12、接受点为S7的原始信号，红色实线代表第一次拾取结果，虚线代表第二次，其结果标注在图上TOA，灰色部分是第二步AIC所考虑的特征函数时间范围；(b)第一次AIC计算结果，红色虚线代表其值最小处；(c)第二次AIC计算结果，红色实线代表其值最小处. Fig. 7 Illustration for picking arrivals (a) is the normalized signal recorded for M12 at S7 and only retains coefficients at high scale. Red dashed line denotes the first picking while the red solid line the second one. The grey rectangular indicates the time window for calculating the second AIC function.(b) is the AIC firstly calculated while the red dashed line denotes its minimum position; (c) is the AIC secondly calculated while the red solid line denotes its minimum position. |
高质量模型的建立需要高质量的数据(Al-Jumaili et al., 2016; Baxter et al., 2007; Eaton et al., 2012; Hensman et al., 2010).低质量的数据，如信噪比较低，或者有其他的干扰信号甚至假信号的地震记录，在自动挑选走时过程中会提供虚假的到时信息，在波形互相关中将引入额外的互相关分量，给到时差模型带来误差.如果模型不能尽可能准确地刻画测区的到时差结构，在此基础上的定位精度也会受影响.
信号的信噪比可以作为我们自动挑选有效数据的一个准则.在信号噪比计算中，信号时窗从到时开始至其后0.08 s，噪声时窗从到时至其前0.2 s.信号噪比定义为这两个时窗中信号的均方根之比.其中信噪比高于6(该值由试验选定)的地震记录自动标注为有效，其后计算得到有效的到时差数据，而由其余低信噪比数据得到的到时差则标注为无效值，用于后续过滤.
对有效数据进行挑选后，会出现部分到时差数据的缺失，处理方法如下：在网格搜索法中，对于某一地震对，如果初始模型网格中某一个结点处没有有效的到时差数据，则该点的值由周围有效的结果插值得到.如果一个模型中有效的结果过少而无法进行散点插值时，定位时便跳过对该模型的引入；在使用GPR时，如果对于某一对地震对没有有效的到时差数据时，一种方法是参考机器学习中如何处理数据丢失情况下的预测方法(Smola, Vishwanathan, Hofmann, 2005)，二是训练若干组GPR进行联合预测(Hensman et al., 2010).本文为简便起见，直接剔除了低信噪比数据所在维度下所有的输入样本，用剩余维度的样本来进行训练.
5 定位结果及分析 5.1 网格搜索法定位结果本实验到时差模型线性散点插值间距为0.5 m. 图 8为S13与S19间的到时差模型.由于两种方法所得到的到时差数据不同，因此插值后的到时差分布图也有一定差异，进而影响后续的定位效果.图 9分别为对应该上述两种DTM模型对T7的定位结果.定位结果处的累积到时差分别为0.8636 s与0.2342 s，定位结果与实际震源位置分别相差1.44 m与1.63 m.最终定位结果见图 10，对应误差汇总见表 1.
|
图 8 S13与S19之间到时差模型图. 到时差计算方法(a)AIC法；(b)波形互相关法 Fig. 8 Delta T model between S13 and S19 calculated from (a) AIC method, (b) cross-correlation |
|
图 9 T7的定位结果 使用的DTM模型由(a)AIC法；(b)波形互相关法. Fig. 9 Locating result for T7 DTM used here are calculated by (a) AIC method, (b) cross-correlation method. |
|
图 10 网格搜索法定位结果汇总 CC代表到时差由采用波形互相关得到，AIC代表由AIC法拾取绝对到时进而计算得到. Fig. 10 Locating results by grid search method CC denotes cross-correlation while AIC the difference of absolute arrivals picked by two-step AIC method. |
|
表 1 网格搜索法的定位误差汇总 Table 1 Localization errors by grid search method |
本方法中，利用上述由两种到时差方法得到的22个台站下、24个模型点位的数据训练GPR.在每种模型中，分别只单独考虑位置的横、纵坐标，将多变量预测分解成两个单变量预测.在得到两个GPR训练模型后，代入测试事件的到时差，分别得到相应的横、纵坐标数值.该方法采用了MATLAB中的统计与机器学习工具箱.图 11是利用GPR得到的定位结果汇总图，表 2是对应的误差汇总.
|
图 11 统计法定位结果汇总 CC代表到时差由采用波形互相关得到，AIC代表由AIC法拾取绝对到时进而计算得到. Fig. 11 Locating results by statistical method CC denotes cross-correlation while AIC the difference of absolute arrivals picked by two-step AIC method. |
|
|
表 2 利用统计定位法得到的定位误差汇总 Table 2 Localization errors by statistical locating method |
DGPS的定位误差以及经纬度转墨卡托坐标的偏差是系统性误差，在本实验的尺度下，对定位结果的影响可忽略不计.而真正影响定位的误差主要来自两个方面.
第一个来源于到时差的测量不够精确.对于波形互相关法，外部相干信号的干扰会使其产生误差.对于双步AIC法拾取走时，时窗的设定不是对每一台都是最优的.这使得在不同的高频噪声水平下，不能最优地区分所有波形的信息与非信息部分，因此在某些地震记录中所拾取的到时会有些许偏差，进而将误差传递、累积至到时差.
第二个来源于模型建立不够准确.在网格搜索法中，采用了线性散点插值法来推测非模型点位处的到时差结构，这与实际情况可能是不相符的，例如地表介质有断裂，波速是非均匀的，或者由地形起伏造成波的传播时长变化.因此在匹配测试事件的到时差结构时产生了误差.在统计定位法中，从样本数目上看，在精确刻画到时差与位置之间映射关系时存在不足，因而存在误差.此外该方法是求基于输入样本的位置的概率分布，而最大概率处是最可能的，与实际比并不一定是最准确的.
5.3.2 方法对比对比两种到时差计算的方法，在本实验中，利用波形互相关计算到时差更为准确.从计算角度而言，波形互相关更为简便快捷，而利用双步AIC法拾取走时需要相当的计算量与参数设置.而从实时处理角度而言，互相关法需要一定时段内的所有数据，而走时拾取法中，可以在仪器中通过设定长短时窗平均的阀值，仅传输爆破信号一段数据，减缓数据传输压力.
对比两种模型的构建，网格搜索法较为直接、方便，并且可以直观地分析测区的到时差结构.但当测区扩大且要保证相当的精度时，该方法需要储存大量的模型数据，并且模型残差的叠加耗时会增长.而对统计定位法来说，尽管训练回归模型需要相当的时间，但定位阶段计算耗时少，适合大尺度范围要求下的实时定位处理.但是该方法需要考虑到缺失数据下的预测，须进一步改进计算方法.而在网格搜索法中，遇到缺失的数据可以略去与之有关的模型的叠加，处理比较简便.总体而言，无论使用何种计算策略，DTM方法具有近实时性的优点.
5.3.3 实验设计优化考虑T7使用的是B01型爆竹，但其定位精度并不明显高于其他测试事件.尽管对于T7，有效的数据更多，但是近处的有效震源记录已经能够很好地约束定位区域，说明实验中DTM是有冗余的，即不须要这么多的接受点.
图 12a与图 12b分别是网格搜索法与统计定位法中，台站数对T7的定位误差的影响，其中到时差由波形互相关得到.台数的选取方式为由近及远(具体顺序为：S19—S17—S21—S15—S22—S6— S10—S14—S13—S9—S11—S16—S20—S5—S7— S18—S1—S2—S3—S4—S8—S12).可以看出前者方法中，台站数目仅需7台，定位结果即可收敛.在后者方法中，台站数目也仅需8台即可达到相应级别的精度.因此说明了定位精度的改善不在于盲目地增加台站数目.
|
图 12 (a)网格搜索法；(b)统计定位法中，对T7定位的误差变化图.其中到时差计算使用波形互相关法 Fig. 12 The change of position error for T7 in (a) grid search method, (b) statistical method. The time difference here is calculated by cross-correlation |
事实上，提高定位精度，更在于增加区域内到时差标定的密度，一方面在网格搜索法中减少线性插值的误差，另一方面在统计定位法中增加学习样本数目.考虑到人工震源实验的成本较高，结合地震波传播路径的可逆性，可通过互换模型点位与接收点，来大幅提高区域内到时差的精度.具体地，在前期到时差标定过程中，将地震仪布置在模型点位处，在测区外围实施少数几次爆破来做标定.对于某一爆破而言，相当于在此放置了一个地震仪，在测区中的各个地震仪位置，分别进行了一次爆破标定.在定位阶段，将相应于爆破标定数的地震仪分别布置在之前的爆破位置，用于接收测试事件的地震记录.值得注意的是，此实验设计中所记录的是针对同一台站(标定期间的炮点)、不同模型点位的地震信号.为了获得针对同一模型点位、不同台站(标定期间的炮点)的到时差，由于此时地震记录并不是在同一时间下接收到的，因此无法使用波形互相关来获得到时差，而需要记录下各标定爆破的精确到时，然后相减获得到时差.因此在前期对到时差标定过程中，须要获得精确的标定爆破时刻，这对建立精确的到时差模型至关重要.
如上述分析，实际用于接受的地震仪数目并不多，因此用于标定的爆破数目也降低了许多，大大节省了人力、财力成本，同时获得了测区内高精度的到时差结构.此外由于台站数减少，相应的定位计算过程中的数据量也大大减少.结合流动布阵，可以很方便地拓展测区.本文实验中，所采用的是普通的爆竹，精确记录爆破时刻是非常困难的.但在大型的工程实验中，这种实验设计，相比传统的DTM方法，拥有上述的优点，非常值得应用.
5.3.4 结合震源扫描法震源扫描法(source-scanning algorithm，后简称SSA)(Kao, Shan, 2004)是利用已知的地震波形，在空间、时间网格进行搜索震源的最优分布的一种到时定位方法.其中，在计算要用于叠加的时间窗的时候，涉及到理论到时的计算，这以往需要使用速度模型.而在本实验中的网格搜索定位法中，已经使用线性散点插值获得了区域内的到时差模型，可以代替速度模型.具体来说，以某一台站的某一时窗的波形为参考，其余任意一台站、在任意一点处，根据到时差模型数据，平移相应的时间量便可直接选取到将要叠加的波形.之后进行叠加、平均，得到亮度值.
本实验中将时间的搜索点设置为S1台站中可分辨的爆破震相最大振幅处，空间搜索步长为0.5 m，计算时波形的邻域选为最大振幅处前后分别为0.1 s和0.2 s，权重函数选为中心点位于最大振幅处的Hanning窗.图 13是T7结合SSA的定位结果，图 14则为定位结果汇总，表 3为定位误差汇总.可以看出结合到时差模型的SSA定位效果也相当不错.
|
图 13 结合SSA的T7定位结果 到时差由(a)波形互相关；(b)AIC方法得到. Fig. 13 Locating result for T7 combining SSA DTM used here are calculated by (a) cross-correlation method, (b)AIC method. |
|
图 14 结合SSA定位结果汇总 CC代表到时差由采用波形互相关得到，AIC代表由AIC法拾取绝对到时进而计算得到. Fig. 14 Locating results by combining SSA CC denotes cross-correlation while AIC the difference of absolute arrivals picked by two-step AIC method. |
|
|
表 3 结合SSA得到的定位误差汇总 Table 3 Localization errors by combining SSA |
SSA方法从计算速度来看远不如之前的两种定位方法，但是这种方法最大的优点是有助于解决多目标定位问题(Liao et al., 2012).之前两种方法在定位单目标的前提，是要得到对应该震源的到时，或者是仅含有该震源信号的波形.而当有多个、连续的震源存在时，各个地震仪上记录到的爆破信号的顺序是未知的，甚至可能混叠在一起，进而无法正确区分，也就不能得到定位结果.而结合SSA则可以一定程度地规避这个问题.
Liao等(2012)中图 1展示了相关的合成数据实验结果，并验证了SSA在多目标定位中的潜力.由于本文实验中各炮点是依次相隔较长时间引爆，对所有台站而言，记录均无混叠、乱序，因而无法给出实际定位效果与更深入的讨论.
6 结论本文运用DTM进行了近场震源的定位实验.测试区域为140 m×90 m.用于模型训练的爆破一共24个，按横向间距28 m，纵向间距30 m均匀分布在测试区域内.用于验证测试的爆破一共10个，随机选择在测试区域.到时差的计算分别运用波形互相关和双步AIC法.定位方法分别采用网格搜索法与统计定位法.不同计算策略下的定位误差约0.5至5.1 m.利用DTM方法并不需要获得测区的速度结构信息，因此适合具有不同地表、地质结构的测区.另外不需要了解地震仪的具体位置，减少了外部测量误差因素的影响.并且在台阵不发生变化的情况下，一个测区的DTM可反复使用.此外还可以添加新的模型数据以加强对测区到时差结构的描述精度.DTM方法并不需要太多的台站.定位精度的改善在于提高测区到时差标定的密度.可以利用地震波传播路径的可逆性，通过交换模型点位与接受点，大幅减少标定操作，提高测区的到时差结构的精度，并减少数据处理量.进一步结合SSA具有解决多目标定位问题的潜力.
这些特点使得DTM方法非常适合于不同野外环境下对震源事件进行高精度、近实时监测.该方法的缺限在于，当台阵发生变化时，原有的DTM便无法再使用；或者当气候、水文等环境条件发生变化，地表介质波速发生相应改变，以致到时差结构变化，此时需要重新学习到时差结构以适应当前情况下的定位作业.
致谢感谢中国科学院地质与地球物理研究所徐锡强、赵剑阳、李少卿提供地震仪以完成此实验，华通信安宇天鹰、唐军在耗材准备、野外作业上的鼎力相助.此外，非常感谢两位审稿专家对文章的润色、分析错误的指出、以及思路的拓展.
| Al-Jumaili S K, Pearson M R, Holford K M, et al. 2016. Acoustic emission source location in complex structures using full automatic delta T mapping technique. Mechanical Systems and Signal Processing, 72-73, 513-524. doi:http://dx.doi.org/10.1016/j.ymssp.2015.11.026. | |
| Allen R. 1982. Automatic phase pickers: their present use and future prospects. Bulletin of the Seismological Society of America, 72(6B): S225-S242. | |
| Asgari S, Stafsudd J Z, Hudson R E, et al. 2015. Moving source localization using seismic signal processing. Journal of Sound and Vibration, 335: 384-396. DOI:10.1016/j.jsv.2014.09.027 | |
| Baxter M G, Pullin R Holford, K M, Evans S L. 2007. Delta T source location for acoustic emission. Mechanical Systems and Signal Processing, 21(3): 1512-1520. DOI:10.1016/j.ymssp.2006.05.003 | |
| Bogiatzis P, Ishii M. 2015. Continuous wavelet decomposition algorithms for automatic detection of compressional-and shear-wave arrival times. Bulletin of the Seismological Society of America, 105(3): 1628-1641. DOI:10.1785/0120140267 | |
| Carter G. 1981. Time delay estimation for passive sonar signal processing. IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(3): 463-470. DOI:10.1109/TASSP.1981.1163560 | |
| Crosson R S. 1976. Crustal structure modeling of earthquake data: 1. Simultaneous least squares estimation of hypocenter and velocity parameters. Journal of Geophysical Research, 81(17): 3036-3046. | |
| Doh-Hyoung K, Youngjin P. 2005. Explicitly-adaptive time delay estimation for wide-band signals. IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, 88(11): 3214-3217. | |
| Eaton M J, Pullin R, Holford K M. 2012. Acoustic emission source location in composite materials using Delta T Mapping. Composites Part A: Applied Science and Manufacturing, 43(6): 856-863. DOI:10.1016/j.compositesa.2012.01.023 | |
| Engdahl E R, Der Hilst R D V, Buland R P. 1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination. Bulletin of the Seismological Society of America, 88(3): 722-743. | |
| Ferris B, Haehnel D, Fox D. 2006. Gaussian processes for signal strength-based location estimation.//Proc. of Robotics Science and Systems. | |
| Hensman J, Mills R, Pierce S, et al. 2010. Locating acoustic emission sources in complex structures using Gaussian processes. Mechanical Systems and Signal Processing, 24(1): 211-223. DOI:10.1016/j.ymssp.2009.05.018 | |
| Ho K, Chan Y T, Ching P. 1993. Adaptive time-delay estimation in nonstationary signal and/or noise power environments. IEEE Transactions on Signal Processing, 41(7): 2289-2299. DOI:10.1109/78.224240 | |
| Kao H, Shan S J. 2004. The source-scanning algorithm: Mapping the distribution of seismic sources in time and space. Geophysical Journal International, 157(2): 589-594. DOI:10.1111/gji.2004.157.issue-2 | |
| Klee U, Gehrig T, McDonough J. 2006. Kalman filters for time delay of arrival-based source localization. EURASIP Journal on Advances in Signal Processing, 2006(1): 1-15. | |
| Liao Y C, Kao H, Rosenberger A, et al. 2012. Delineating complex spatiotemporal distribution of earthquake aftershocks: An improved source-scanning algorithm. Geophysical Journal International, 189(3): 1753-1770. DOI:10.1111/gji.2012.189.issue-3 | |
| Maeda N. 1985. A method for reading and checking phase times in auto-processing system of seismic wave data. Zisin, 38(3): 365-379. DOI:10.4294/zisin1948.38.3_365 | |
| Matsu, Apos, Ura M. 2009. Bayesian estimation of hypocenter with origin time eliminated. Journal of Physics of the Earth, 32(6): 469-483. | |
| Nakata N, Beroza G C. 2016. Reverse time migration for microseismic sources using the geometric mean as an imaging condition. Geophysics, 81(2): KS51-KS60. DOI:10.1190/geo2015-0278.1 | |
| Nelder J A, Mead R. 1965. A simplex method for function minimization. The Computer Journal, 7(4): 308-313. DOI:10.1093/comjnl/7.4.308 | |
| Rasmussen C E. 2006. Gaussian processes for machine learning. Massachusetts: The MIT Press. | |
| Saragiotis C D, Hadjileontiadis L J, Panas S M. 2002. PAI-S/K: A robust automatic seismic P phase arrival identification scheme. IEEE Transactions on Geoscience and Remote Sensing, 40(6): 1395-1404. DOI:10.1109/TGRS.2002.800438 | |
| Sedlak P, Hirose Y, Enoki M, et al. 2008. Arrival time detection in thin multilayer plates on the basis of Akaike information criterion. Journal of Acoustic Emission, 26: 182-189. | |
| Smola A J, Vishwanathan S V N, Hofmann T. 2005. Kernel methods for missing variables.//The Tenth International Workshop on Artificial Intelligence & Statistics. | |
| Spence W. 1980. Relative epicenter determination using P-wave arrival-time differences. Bulletin of the Seismological Society of America, 70(1): 171-183. | |
| Thanasopoulos I, Avaritsiotis J. 2011. Wavelet analysis of short range seismic signals for accurate time of arrival estimation in dispersive environments. IET Science, Measurement & Technology, 5(4): 125-133. | |
| Waldhauser F, Ellsworth W L. 2000. A double-difference earthquake location algorithm: Method and application to the northern Hayward fault, California. Bulletin of the Seismological Society of America, 90(6): 1353-1368. DOI:10.1785/0120000006 | |
| Walker K T, Ishii M, Shearer P M. 2005. Rupture details of the 28 March 2005 Sumatra MW8.6 earthquake imaged with teleseismic P waves. Geophysical Research Letters, 32(24): L24303. DOI:10.1029/2005GL024395 | |
| 范建柯, 丁志峰, 徐小明, 等. 2016. 基于自适应量子遗传算法的地震重定位及其在山东地区的应用. 地球物理学报, 59(11): 4075–4088. DOI:10.6038/cjg20161112 | |
| 高星, 王卫民, 姚振兴. 2002. 用于地震定位的SAMS方法. 地球物理学报, 45(4): 525–532. | |
| 李志伟, 胥颐, 郝天珧, 等. 2006. 利用DE算法反演地壳速度模型和地震定位. 地球物理学进展, 21(2): 370–378. | |
| 刘劲松, 王赟, 姚振兴. 2013. 微地震信号到时自动拾取方法. 地球物理学报, 56(5): 1660–1666. DOI:10.6038/cjg20130523 | |
| 田玥, 陈晓非. 2002. 地震定位研究综述. 地球物理学进展, 17(1): 147–155. | |
| 赵珠, 丁志峰, 易桂喜, 王建格. 1994. 西藏地震定位──一种使用单纯形优化的非线性方法. 地震学报(2): 212–219. | |