基于信號包絡線相似性的特定場地地下爆炸識別方法研究

潘常周 靳 平 徐 雄 王紅春 徐恒壘

（中國西安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.