基于信号包络线相似性的特定场地地下爆炸识别方法研究

潘常周 靳 平 徐 雄 王红春 徐恒垒

（中国西安710024西北核技术研究所）

引言

区分地下爆炸与天然地震事件的识别技术是核爆炸地震监测领域的重点研究方向之一.早期的识别技术侧重于远震记录分析，并发展了重要的识别判据，如震级差mb／MS判据（Marshall，Basham，1972）.利用远震记录判据的最大不足是难以用来识别中小震级事件.因为通常情况下中小震级事件仅在区域记录上有较清晰的信号，而远震台站记录的信号则很弱，导致难以利用远震记录可靠地测量信号的特征参数.以mb／MS判据为例，mb＜4.0的绝大多数地震事件在远震记录的面波信号都不发育，从而难以估算面波震级MS.

为了提高对中小震级事件的识别能力，地震事件识别技术研究从20世纪70年代末开始侧重于分析区域记录.例如：Blandford（1982）和Pomeroy等（1982）基于有限的数据归纳了具有较好识别能力的多种区域性识别判据；Taylor等（1989）利用更多地震事件验证了其中部分判据的识别能力，如P／S震相幅值比、震相谱比值等判据.这些区域性识别判据以通常具有较高信噪比的区域震相的特征参数，如不同频带内P波与S波的幅值比、相同震相的不同频带幅值比等为区分事件的特征值.其主要优点在于，即使对低震级事件判据的特征值也容易测量.

由于区域震相的传播路径主要是介质横向不均匀性很强的地壳和上地幔，区域性判据的特征值容易受震相传播路径的影响.因此，针对特定地区的标定和校正判据特征值，是提高判据识别正确率的关键.常用的标定和校正方法是通过拟合特征值与震中距的关系（Patton，Taylor，1995；Walter et al，1995；Hartse et al，1997；Pan et al，2007a），建立一维的校正曲线.而这无法考虑不同方向路径差异的影响，因此一些研究人员（Taylor，Hartse，1998；Fisk，Bottone，2000；Pan et al，2007b）在此基础上利用Kriging（Schultz et al，1998）等技术建立二维的校正曲面，有效地提高了事件识别的正确率.近年来，一些研究（Taylor et al，2002；Fisk et al，2007，2008；Anderson et al，2009；Leidig et al，2010）根据介质品质因子Q值结构模型和地震事件源频谱模型来校正区域震相的频谱，即MDAC（magnitude and distance amplitude correct）方法，进而据此计算判据特征值，取得了更好的识别效果.

尽管如此，判据特征值标定模型的误差仍然是导致上述区域性判据识别错误的重要原因之一.近年来一些研究（Gibbons，Ringdal，2004，2006；Schaff，Waldhauser，2010；王红春等，2012）利用波形互相关方法识别特定场地爆炸事件，避免了识别判据特征参数标定建模和校正的问题，其错误识别率接近于零.虽然这种方法仅仅适用于识别特定场地的爆炸事件，但多数地下核爆炸实际上也是在特定试验场反复多次实施的，因此这种方法也有着重要的应用价值.不过，基于波形互相关的识别方法能够应用的特定场地的空间范围太小，主要原因是波形相关系数随事件间距的增加呈指数衰减（Menke et al，1990；Hough，Field，1996），当事件间距超过1／4波长时波形相关系数一般比较低（Geller，Mueller，1980）.根据Schaff和 Waldhauser（2010）以及王红春等（2012）的研究结果，利用波形互相关法只能识别距离已知爆炸3km以内的事件.然而，目前已知的主要核试验场的面积都比较大，同一试验场很多历史爆炸试验的间距都远大于该距离，达到十几千米甚至数十千米.因此，进一步扩大这种识别方法的适用空间范围，对于提高其在核爆炸监测识别中的实际应用价值具有重要意义.

包络线是信号的重要特征之一，在地震信号自动检测与关联、介质结构反演等研究（Ryzhikov et al，1996；王勤彩等，2009）中均有较好的应用.作者在地下爆炸的分析研究中发现，一些爆炸在间距达到十余千米、波形相似性很弱的情况下，其包络线仍然保持着较高的相似性.因此，本文搜集分析了更多的相关事件，研究基于信号包络线相似性的特定场地的爆炸识别方法，以期更好地应用于地下核爆炸的监测识别.

1 原理方法

一般情况下，特定场地地震事件到相同台站的传播路径大致相同，同一震相的能量衰减及走时都基本一致；当震源机制相同或相近时，在任何空间方位激发的主要震相（P波、S波）的能量相对大小也是基本一致的，因此它们在相同台站上主要震相整体发育特征比较相似.与原始波形相比，波形的能量包络线忽略了局部细微的特征而主要反映更宏观的震相整体发育特征.集中式地下爆炸可以近似认为是膨胀体积源.因此，对于特定场地的地下爆炸，随着爆心距的增大，当它们在相同台站的波形因传播路径局部差异导致相关系数较低的时候，其信号包络线很可能仍然具有较高的相关系数.而天然地震一般是变化复杂的剪切位错源，与地下爆炸有很大差别.其激发主要震相能量的相对大小可能会在个别方位上比较接近，但是在绝大部分方位上均存在较大差异.因此，对于特定场地的地下爆炸与天然地震的波形包络线，有可能在少数台站上具有较强相似性，但在大多数台站上其相似性应该是比较弱的.因此，本文拟利用台网平均信号包络线的相关系数来识别特定场地地下爆炸，具体识别方法如下.

首先，计算事件的信号包络线.信号包络线在许多针对不同应用目的的研究中，其计算方法并不完全相同.针对本文研究目的，需要计算的是反映信号能量的包络线.其计算方法为：对波形序列x（t），分别采用适当带通滤波器滤波后进行希尔伯特变换得到h（t），对h（t）取虚部得到序列y（t），然后按下式计算x（t）的包络线：

并对其进行光滑处理：取长度为1s的滑动汉宁窗，计算窗口内信号的方均根幅值.信号窗口取为从P波开始至Lg波结束（Lg波取2.8—3.7km／s的速度窗口），包含了主要区域震相.记垂直向、东西向、南北向的包络线分别为eZ（t）、eE（t）、eN（t），分别按式（2）和式（3）计算水平向和三分向的包络线为

然后，按相应分向计算已知爆炸与待识别事件在相同台站的信号包络线的相关系数.当两个事件在某台站的第i个频带的信号包络线（分别记为ei1（t）和ei2（t））同时满足一定信噪比要求时，其相关系数为

式中，¯ei1（t）和¯ei2（t）分别为ei1（t）和ei2（t）的平均值.考虑到估算震相到时的误差，我们把已知爆炸信号包络线当作模板，通过在时间域内逐点滑动的方式，在两个事件P波到时对齐前后1s时间范围内滑动计算相关系数，再取该范围内相关系数的最大值.如果选用多个频带分析，则按下式计算不同频带相关系数的平均值：

最后，根据已知爆炸与待识别事件在不同台站的信号包络线的相关系数，按式（5）计算信号包络线相关系数台网平均值，并以此作为筛选识别事件的特征参数，即台网平均信号包络线相关系数大于一定阈值时判别事件为地下爆炸，否则为天然地震.

图1为间距约30km的两次地下爆炸的波形及其包络线以及波形的相关系数和包络线的相关系数，图2为其中一次爆炸与天然地震（间距约18km）的波形及其包络线以及波形的相关系数和包络线的相关系数.可以看出，虽然两次爆炸的间距更大，但它们在波形相似性很弱的情况下，信号包络线相关系数仍然较高，在0.6—0.9之间；而一次爆炸与天然地震的信号包络线相关系数则较低，仅在个别台站接近0.5，多数台站低于0.2.因此，根据与已知爆炸包络线的台网平均相关系数，有望能够有效地区分未知的地下爆炸与天然地震.

图1 两次地下爆炸的波形及其包络线相似性比较（a）2010-01-28爆炸波形及其包络线；（b）2009-04-09爆炸波形及其包络线；（c）波形相关系数；（d）包络线相关系数Fig.1 Comparison of waveforms of two explosions，as well as correlation between their envelopes（a）Waveforms and envelopes of the explosion on 28January 2010；（b）Waveforms and envelopes of the explosion on 9April 2009；（c）Correlations between pairs of waveforms；（d）Correlations between pairs of envelopes

2 应用效果分析

新疆维吾尔自治区托克逊县及其附近地区有多个采矿爆破场，同时该地区的地震活动性也较强，有较多位置相近的地下爆炸与天然地震.徐雄（2011）曾对该地区2009—2010年ML≥1.5事件进行聚类分析，并确定了部分人工爆炸.我们利用该地区地震事件来评估本文提出方法的识别效果.本文选用新疆地震局2009—2010年地震目录中以2010-04-15爆炸震中为中心、45km为半径、ML≥1.5的事件，共527个（图3）.根据徐雄（2011）的分析结果，这些事件中有418个被确认为地下爆炸（最大震级ML2.9），另外有109个为未经核实分析的未知事件（最大震级为ML3.5）.本文选取2010-04-15爆炸作为已知爆炸，分别按垂直向、水平向和三分向计算它与其余526个事件的信号包络线相关系数，对这些事件进行识别，并统计误识别率.

为了尽量准确地评估有效识别事件的间距范围，我们利用双差定位法（Waldhauser，Ellswort，2000）对上述418次地下爆炸进行了相对定位，结果如图4所示.待识别分析的爆炸事件与2010-04-15爆炸的最大间距约为30km.识别分析所用的台站共17个（图3）.这些台站对事件形成了较好的包围效果，震中距（相对于2010-04-15爆炸）约为70—350km.事件和台站分布总体上满足本文所提出的识别方法的基本要求，即监测台站震中距应远大于事件间距，而且对事件形成较好的包围效果.不过对于具体事件而言，能够记录到信号的台站未必满足上述基本要求.为此，本文要求进行筛选识别的事件需满足以下条件：①包络线P波信噪比大于2；②震中距为事件间距的5倍以上；③符合前两条要求的台站在3个以上并且能够对事件形成较好包围效果 （这里采用反映台站缺失最大方位范围的最大间隙角作为限定，要求最大间隙角小于240°）.

图3 识别分析的地震事件及台站分布Fig.3 Distribution of seismic events for identification and stations Triangles denote stations，asterisk denotes the explosion on 15April 2010，red dots denote other explosions，and blue dots denote unknown events

图4 地下爆炸事件的相对位置Fig.4 Relative location of known explosions The asterisk denotes the explosion on 15April 2010，and blue dots denote other explosions to be identified

根据上述要求，符合筛选识别条件的事件共有317个，约占事件总数的60%；未达到识别条件的主要是震级较小的事件，约87%为ML＜2.0的事件.识别结果如图5所示，可以看出：垂直向的识别效果最好，绝大部分爆炸事件与2010-04-15爆炸的包络线相关系数较高，特别是事件间距在15km以内的爆炸，多数的相关系数大于0.5；同时大多数未知事件与2010-04-15爆炸的相关系数均较低.图5中相关系数较高的两个异常未知事件，很有可能确实是地下爆炸.原因有两方面：一是其波形主要震相发育特征与2010-04-15爆炸比较相似，具有P波较强S波较弱的爆炸波形特征（图6）；二是这两个事件所属的重复事件群（总共12个事件，最大震级为ML1.8，其它事件信号未达到识别条件）中所有事件的发震时间都集中在新建地区经常进行工业爆破的时段内（北京时间16—21时），说明不大可能是天然地震.不过，要确认这两个事件为地下爆炸，还需要更多的证据支持.

图5 台网平均信号包络线相关系数与事件间距的关系Fig.5 Relationship between envelope correlation coefficient of net average and inter-distance of corresponding event-pair.The dots denote known explosions，and the squares denote unknown events

图6 异常的未知事件与已知爆炸波形（4—8Hz）的对比Fig.6 Comparison of waveforms of known explosion with those of abnormal unknown events（4—8Hz band-pass filtered）

本文从保守评估判据识别能力的角度考虑，把所有未知事件都当作天然地震.根据垂直向的结果统计不同相关系数判别该阈值条件下的误警率（天然地震中被识别为爆炸的比例）和漏检率（地下爆炸中未能正确识别的比例）.假设以台网平均的包络线相关系数c作为判别阈值，即相关系数大于c的事件判别为地下爆炸，小于c的事件判别为天然地震，则可分别根据式（6）和（7）统计误警率Rf和漏检率Rm.式中，Ne和Nx分别为天然地震和地下爆炸的总数，Nxe为被识别为地下爆炸的天然地震数，Nex为被识别为天然地震的地下爆炸数.结果如图7所示，可以看出：当相关系数阈值约为0.38时，误警率约为6.3%；间距在15km以内爆炸漏检率仅为0.5%，间距在15—30 km之间的爆炸漏检率约为5.0%.

图7 误警率和漏检率与识别阈值的关系Fig.7 Relationship between wrong alert rate／missing rate and correlation threshold

3 结论与讨论

本文分析结果表明，利用信号包络线相似性能够较好地识别特定场地的地下爆炸.以1次已知历史地下爆炸为模板，根据它与其它事件信号包络线的台网平均相关系数对事件进行识别.对于距离模板事件15km以内的地下爆炸，仅有1次被错误识别，漏检率约为0.5%；对于距离模板事件15—30km的地下爆炸，有5次被错误识别，漏检率约为5.0%.所有未经核实性质的事件中，有2个被识别为爆炸.根据波形特征和发震时间等特点分析，这两个事件很有可能确实属于地下爆炸.即使这两个事件是属于被错误识别的天然地震，其误警率也仅为6.3%.而P／S幅值比等判据的误警率和漏检率分别约为10%和5%（Pan et al，2007a）.波形互相关识别方法的误警率虽然低至接近于零，但其只能识别事件间距小于3km的爆炸（Schaff，Waldhauser，2010；王红春等，2012）.因此，与波形互相关方法相比，本文识别方法大幅提高了有效识别特别场地爆炸事件的空间范围，对实际应用具有重要意义.虽然本文方法错误识别率略为偏高，但仍然低于P／S幅值比等常用识别判据.

同时，在上述研究中也发现了一些值得讨论的问题.第一个问题是关于利用信号包络线相似性识别特定场地爆炸事件的有效空间范围.本文分析结果表明，30km以内识别效果是比较好的，但对更大间距事件的识别效果尚不清楚.因为本文所分析的爆炸震级都很小，震中距较大的台站一般没有记录信号，只能利用震中距较小的台站，不适合分析间距更大的事件.而且，本文只针对一个场地进行了分析.而对于不同场地，有效识别空间范围可能受地质条件（如地形起伏等）的影响.另外，本文分析的是只有1个已知爆炸的情况，在实际监测中，敏感场地可能有多个分散的已知爆炸.在这种情况下，综合利用这些已知爆炸，可以对特定场地内较大区域范围的事件进行筛选识别，而且可能会取得更好的识别效果.

第二个问题是关于利用信号包络线相似性是否能够约束天然地震的震源机制.上述分析结果表明：由于震源机制相同，多次地下爆炸之间信号包络线具有较好的相似性；而由于震源机制不同，则地下爆炸与天然地震之间信号包络线相似性较差.那么，震源机制相同或相近的天然地震之间信号包络线是否也具有较好相似性，而震源机制不同的天然地震之间信号包络线相似性也较差呢？如果答案是肯定的，说明信号包络线相似性对震源机制有较强的约束力，则可以利用信号包络线反演震源机制.这对于中小震级地震来说非常重要，因为中小震级地震信号的优势频率相对比较高，利用波形反演震源机制就需要高精度介质结构模型来支持计算高频理论地震图.而目前甚至在将来很长一段时间内，介质结构模型的精度都很难满足这样的需求.如果关注的是反映主要震相整体发育特征的包络线，则计算理论地震图对介质结构模型精度要求就没那么高.当然，与地下爆炸相比，天然地震情况要复杂得多.除了震源机制外，地震信号包络线相似性还可能受其它因素的影响，如震源深度、断层破裂方向等.因此，能否利用信号包络线反演震源机制，还需要开展更深入的研究.

王勤彩，陈章立，Asano Y，郑斯华，Hasegawa A.2009.利用尾波包络线反演方法研究伽师强震群区地壳的非均匀结构［J］.地球物理学报，52（1）：90-98.

Wang Q C，Chen Z L，Asano Y，Zheng S H，Hasegawa A.2009.Imaging crustal heterogeneity in Jiashi strong earthquake swarm region by coda envelope inversion analysis［J］.Chinese Journal of Geophysics，52（1）：90-98（in Chinese）.

王红春，靳平，何燕.2012.基于三分向台站波形的重复地下爆炸相关检测［J］.地球物理学报，55（3）：937-943.doi：10.6038／j.issn.0001-5733.2012.03.023

Wang H C，Jin P，He Y.2012.Cross-correlation detection of repeating underground explosions using three-component stations［J］.Chinese Journal of Geophysics，55（3）：937-943.doi：10.6038／j.issn.0001-5733.2012.03.023（in Chinese）.

徐雄.2011.地震数据库中重复地震事件搜索方法研究及应用［D］.西安：西北核技术研究所：74.

Xu X.2011.Study of Algorithms for Repeat Events Exploration in Seismic Database and Its Application［D］.Xi’an：Northwest Institute of Nuclear Technology：74（in Chinese）.

Anderson D N，Walter W R，Fagan D K，Mercier T M，Taylor S R.2009.Regional multistation discriminants：Magnitude，distance，and amplitude corrections，and sources of error［J］.Bull Seismol Soc Am，99（2A）：794-808.

Blandford R R.1982.Seismic event discrimination［J］.Bull Seismol Soc Am，72（6B）：S69-S87.

Fisk M D，Bottone S.2000.Regional seismic event characterization using a Bayesian Kriging approach［C］∥22nd Annual Seismic Research Symposium.Planning for Verification of and Compliance with the Comprehensive Nuclear-Test-Ban Treaty （CTBT）.New Orleans：Department of Defense and U.S.Department of Energy：01-03.

Fisk M D，Taylor S R，Patton H J，Walter W R.2007.Next-generation MDAC discrimination procedure using multidimensional spectral analyses［C］∥Proceedings of the 29th Monitoring Research Review：Ground-Based Nuclear Explosion Monitoring Technologies.Denver，Colorado：National Nuclear Security Administration，LA-UR-07-5613，1：551-560.

Fisk M D，Taylor S R，Patton H J，Walter W R.2008.Applications of a next-generation MDAC discrimination procedure using two-dimensional grids of regional P／S spectral ratios［C］∥Proceedings of the 30th Monitoring Research Review：Ground-Based Nuclear Explosion Monitoring Technologies.Portsmouth，Virginia：National Nuclear Security Administration，LA-UR-08-05261，1：583-592.

Geller R J，Mueller C S.1980.Four similar earthquakes in central California［J］.Geophys Res Lett，7（10）：821-824.

Gibbons S J，Ringdal F.2004.A waveform correlation procedure for detecting decoupled chemical explosions［C］∥NORSAR Scientific Report：Semiannual Technical Summary （2）.Kjeller：Norwegian Seismic Array：41-50.

Gibbons S J，Ringdal F.2006.The detection of low magnitude seismic events using array-based waveform correlation［J］.Geophys J Int，165（1）：149-166.doi：10.1111／j.1365-246X.2006.02865.x.

Hartse H E，Taylor S R，Phillips W S，Randall G E.1997.A preliminary study of regional seismic discrimination in central Asia with emphasis on western China［J］.Bull Seismol Soc Am，87（3）：551-568.

Hough S E，Field E H.1996.On the coherence of ground motion in the San Fernad Valley［J］.Bull Seismol Soc Am，86（6）：1724-1732.

Leidig M，Stump B，Walter W R，Stroujkova A，Yang X，Bonner J.2010.Examination of P／S spectral ratios for small explosions at local distances and interpretation of moment tensors estimated from near-source data［C］∥Proceedings of the 32th Monitoring Research Review：Ground-Based Nuclear Explosion Monitoring Technologies.Orlando，Florido：National Nuclear Security Administration，LA-UR-10-05578，1：416-426.

Marshall P D，Basham P W.1972.Discrimination between earthquakes and underground explosions employing an improved MSscale［J］.Geophys J R astr Soc，28（5）：431-458.

Menke W，Lerner-Lam A L，Dubendorff B，Pacheco J.1990.Polarization and coherence of 5to 30Hz seismic wave fields at a hard-rock site and their relevance to velocity heterogeneities in the crust［J］.Bull Seismol Soc Am，80（2）：430-449.

Pan C Z，Jin P，Wang H C.2007a.Applicability of P／S amplitude ratios for the discrimination of low magnitude seismic events［J］.Acta Seismologica Sinica，20（5）：553-561.

Pan C Z，Jin P，Xiao W G.2007b.Calibration of P／S amplitude ratios for seismic events in Xinjiang and its adjacent areas based on a Bayesian Kriging method［J］.Acta Seismologica Sinica，20（6）：664-674.

Patton H J，Taylor S R.1995.Analysis of Lg spectral ratios from NTS explosions：Implications for the source mechanisms of spall and the generation of Lg waves［J］.Bull Seismol Soc Am，85（1）：220-236.

Pomeroy P W，Best W J，McEvilly T V.1982.Test ban treaty verification with regional data：A review［J］.Bull Seismol Soc Am，72（6B）：S89-S129.

Ryzhikov G A，Biryulina M S，Husebye E S.1996.A novel approach to automatic monitoring of regional seismic events［J］.IRIS News Letter，15（1）：12-14.

Schaff D P，Waldhauser F.2010.One magnitude unit reduction in detection threshold by cross correlation applied to Parkfield（California）and China seismicity［J］.Bull Seismol Soc Am，100（6）：3224-3238.

Schultz C A，Myers S C，Hipp J，Young C J.1998.Nonstationary Bayesian Kriging：A predictive technique to generate spatial correction for seismic detection，location，and identification［J］.Bull Seismol Soc Am，88（5）：1275-1288.

Taylor S R，Denny M D，Vergino E S，Glaser R E.1989.Regional discrimination between NTS explosions and western U.S.earthquakes［J］.Bull Seismol Soc Am，79（4）：1142-1176

Taylor S R，Hartse H E.1998.A procedure for estimation of source and propagation amplitude corrections for regional seismic discriminants［J］.J Geophys Res，103（B2）：2781-2789.

Taylor S R，Velasco A A，Hartse H E，Phillips W S，Walter W R，Rodgers A J.2002.Amplitude corrections for regional seismic discriminants［J］.Pure Appl Geophys，159（4）：623-650.

Waldhauser F，Ellswort W L.2000.A double-difference earthquake location algorithm：Method and application to the Northern Hayward fault，California［J］.Bull Seismol Soc Am，90（1）：1353-1368.

Walter W R，Mayeda K M，Patton H J.1995.Phase and spectral ratio discrimination between NTS earthquakes and explosion，PartⅠ：Empirical observations［J］.Bull Seismol Soc Am，85（4）：1050-1067.