詹小艷 王恒知 王俊 繆發(fā)軍 薛瑩瑩 朱升初

摘要：以聲譜圖像模式識別為基本思想，改進了基于地震波波譜包絡分布特征的初至震相檢測方法。改進后的算法，首先將波形記錄特征轉換為譜能量分布，以臺站固有背景噪聲的時頻特征為基準，提取出記錄信號的時頻包絡特征作為事件觸發(fā)檢測的目標函數(shù)。為降低主震尾波對檢測后續(xù)余震的影響，采用非線性標度縮放因子對大于臺站背景噪聲水平的能量進行歸一化處理。通過2組密集余震序列的實測結果表明：相對STA/LTA算法自動檢測結果，新方法正確觸發(fā)率提高約19.6%，漏觸發(fā)率則降低了22.6%，表現(xiàn)出更強的抗干擾能力和觸發(fā)效率，且運行速度快，能夠滿足地震實時處理系統(tǒng)的技術要求。

關鍵詞：余震序列；初至震相；地震波；波譜能量包絡；峰值搜索

中圖分類號：P315.3 文獻標識碼：A 文章編號：1000-0666（2018）02-0258-06

0 引言

完整、可靠的余震序列參數(shù)是中強地震發(fā)生后或密集震群序列活動時，確定主震發(fā)震構造、估計斷層破裂尺度、序列后續(xù)趨勢判定以及開展應急救援等工作關鍵的基礎數(shù)據(jù)。研究表明，余震序列精定位是確定主震發(fā)震構造特征最為有效的途徑之一（黃媛等，2008；朱艾斕等，2008）：利用余震序列的震源機制解還可以進一步反演震源區(qū)構造應力場的分布特征（王勤彩等，2009）。然而，地震學家們通常需要滯后數(shù)天甚至更長的時間才能得到有關強余震序列較為完整、可靠的后續(xù)研究結果。其主要原因是目前地震實時處理系統(tǒng)所采用的STA/LTA算法對余震序列的檢測效率較低，絕大部分余震需要震相分析人員通過查看連續(xù)記錄波形進行查找。隨著臺站密度的日益提高，密集余震序列快速檢測的數(shù)據(jù)量也將進一步增加。因此，研究一種快速、高效的地震事件觸發(fā)檢測算法，具有十分重要的現(xiàn)實意義。

自20世紀70年代以來，地震震相的觸發(fā)檢測一直都是國內(nèi)外學者研究的重點領域，并取得了大量的成果。從最簡單的以振幅大小作為觸發(fā)標準到復雜的圖像模式識別、自適應方法、神經(jīng)網(wǎng)絡法以及小波變換的主成分分析法等（Allen，1978；Joswig，1990；Wang，Teng，1995；Witherset al，1998；Patane et al，1999；Sleeman et al，1999；Leonard et al，1999，2000；劉希強等，2000；Massa et al，2006）。這些震相識別方法是根據(jù)信號與噪聲之間的差異特征而提出，如振幅或能量、頻率、波形的相似性、地震波傳播方向和動力學特征（如偏振、頻譜等），對于單個地震或特定的區(qū)域性地震而言，可取得較好效果。但這些復雜、有效的震相觸發(fā)檢測算法，通常由于運算速度較慢而難以滿足實時處理的技術要求。

通常情況下，中、強地震發(fā)生后短時間內(nèi)會在震源區(qū)附近發(fā)生大量的余震，絕大多數(shù)余震發(fā)震時間間隔極短，甚至相互重疊且處于主震的尾波中。例如，據(jù)國家地震臺網(wǎng)中心統(tǒng)一編目結果，2008年汶川M8.0地震后50d內(nèi)，四川地震臺網(wǎng)以及應急流動臺站共記錄到余震總數(shù)超過1.6萬次，余震最密集時1h內(nèi)ML>1.0地震達179次；2014年新疆于田M7.3地震后5d內(nèi)共記錄余震總數(shù)3294次。盡管震相檢測觸發(fā)方面已有不少研究成果，然而針對地震序列快速監(jiān)測的研究成果仍然較少。由于STA/LTA算法（短時平均與長時平均值之比）具有算法簡單、速度快、便于實時處理等特點（Stevenson，1976），在地震事件的實時觸發(fā)檢測中，國、內(nèi)外主要地震機構的處理系統(tǒng)仍普遍采用該算法。STA/LTA算法對于環(huán)境背景噪聲水平低的臺站來說非常有效，但它的抗干擾能力不夠強，若記錄中存在無規(guī)則、隨機大振幅的干擾信號、或脈沖式噪聲以及信噪比較差時，該算法的效率會大大降低。另一方面，在檢測過程中，STA/LTA的觸發(fā)閾值起了決定性的作用，閾值設置得高，小地震則不會觸發(fā)；而當閾值設置得低時，則會出現(xiàn)大量的誤觸發(fā)。盡管在實際計算過程中有2種方案可供選擇，即“鎖定”LTA和連續(xù)更新LTA，但這2種方法仍有優(yōu)缺點：“鎖定”LTA值使記錄系統(tǒng)一直處于觸發(fā)狀態(tài)，而連續(xù)更新LTA值時，會使尾波過早地被截斷（Bor-mann，2006）。此外，采用STA/LTA算法檢測余震序列中的初至時，主震尾波及余震之間的相互疊加就成為干擾STA/LTA算法的主要因素。

本文將針對密集余震序列的記錄特征，以聲譜圖像模式識別方法為基本思想，實驗一種基于地震波波譜包絡分布特征為基礎的初至震相快速檢測方法，并將其應用于2次地震進行驗證。

1 研究方法

Joswing（1990，1993）提出的聲譜圖像模式識別（Sonogram Pattern Recognition）方法，其基本思想是將各種類型的地震信號與臨時噪聲信號定義成不同的模式，在實時檢測過程中，當與定義模式足夠相似時，即被觸發(fā)，該方法是一種正演的決定邏輯觸發(fā)模式。觸發(fā)模式或觸發(fā)閾值的設置不依賴于地噪聲水平，因為它在此階段對地震信號和突發(fā)性噪聲采用相同的處理方式。基本步驟為：（1）將記錄波形進行非線性化的圖像轉換，即通過信號的譜能量構建二維灰度圖像，其中黑度代表潛能量的密度，在處理過程中采用了非線性的標度縮放因子，使它與簡單的譜密度圖像有著很大的差異：（2）基于典型的地震信號和噪聲建立不同的識別模式。如建立不同類型地震信號與交通、工業(yè)噪聲、汽車、火車等大幅度噪聲之間的模式等；（3）實時檢測，即實現(xiàn)檢測閾值和最相似模式的邏輯決定。其中（1）、（2）步中構建波形的非線性圖像和識別模式難度較大且過程十分復雜，這也是該方法難以被廣泛應用的原因所在。

本文是在上述基本思路的基礎上，將其改進為波譜能量時頻包絡線識別模式，基本步驟為：

（1）沿記錄波形連續(xù)滑動地計算一定窗長內(nèi)信號的功率密度譜（PSD），在全頻段內(nèi)以1/8倍頻為單位間隔計算平均功率譜，即計算在短拐角周期Ts和長拐角周期TL=2×Ts內(nèi)的平均功率譜值，賦值給幾何中心周期Tc[Tc=sqrt（Ts×TL）]，記為矩陣Aij：

Aij=psd（fi，tj）（1）式中：fi代表第i個頻率，tj代表j時刻的記錄波形。

為了最大限度地突出地震信號的能量，同時避開噪聲能量最強的頻段（Peterson噪聲模型中微震峰值頻帶1～10s），需要在計算前根據(jù)不同地震信號（P、S波）的卓越頻帶范圍進行相應的濾波處理。

（2）記錄信號的功率譜密度矩陣與該臺穩(wěn)定的噪聲功率譜矩陣nij之比，記為Bij：

Bij=Aij/nij（2）式中：nij是根據(jù)每個臺站背景噪聲的功率譜概率密度函數(shù)計算得到。對于單個臺站的記錄，在某一頻段內(nèi)，將信號中高于該臺背景噪聲功率譜閾值矩陣的部分（即Bij中大于1處），通過非線性的標度縮放因子進行歸一化處理，以降低大振幅信號對能量分布權重的影響，換言之，即提高弱信號（尤其是淹沒在背景噪聲中的信號）在能量分布中的比重。于是可根據(jù)高于臺站背景噪聲水平、歸一化后的能量分布和持續(xù)時間，提取出地震記錄波形在特定頻帶段內(nèi)的譜能量時頻包絡曲線。

（3）通過設定譜能量的觸發(fā)閾值對信號進行觸發(fā)檢測，這樣觸發(fā)檢測就完全取決于臺站隨頻率的噪聲水平和各個頻點上高于背景噪聲水平的持續(xù)時間，在持續(xù)時間內(nèi)的信號，即認為是有效的地震信號。

2 實測結果

2.1 江蘇高郵M4.9地震余震序列

2012年7月20日20時11分，在江蘇高郵與寶應交界處發(fā)生M4.9地震，這是近30年來江蘇省陸地上發(fā)生的最大地震。距離主震震中最近的是寶應地震臺，震中距約為36.8km。據(jù)江蘇省地震臺網(wǎng)記錄，震后4h內(nèi)的寶應臺共記錄余震44次，其中震后1h內(nèi)的余震為29次（表1）。為驗證本文提出的方法，選取了震后1h內(nèi)的連續(xù)記錄波形進行回溯性實測。在提取波譜能量包絡線之前，先確定被檢測臺站的各頻點上背景噪聲水平nij。本文中采用功率譜概率密度函數(shù)法計算nij（McNamara，Boaz，2005；王俊等，2009）。取寶應臺背景噪聲加速度功率譜概率密度函數(shù)80%的分布概率代表其背景噪聲水平，如圖1中黃色線條所示，在1～20Hz頻帶范圍內(nèi)的平均值約為-145.7dB（噪聲水平的有效值約為5.16×10-8m/s）。

在震中距小于100km的地方震記錄中，P波的周期約為0.05～0.2s，同時為避開背景噪聲微震峰值的影響，在計算前對連續(xù)記錄進行1Hz以上的Butterworth高通濾波處理。首先，將垂直向記錄波形中高于各頻點背景噪聲水平的能量進行歸一化處理，如圖2b中紅色區(qū)域所示，白色曲線即為根據(jù)時頻能量分布提取出的特征包絡曲線。對比實際記錄波形（圖2a）可以看出，在經(jīng)過歸一化后的能量包絡曲線中，由于大振幅信號在能量分布中的權重被削弱了，除了可以清晰地展現(xiàn)出記錄數(shù)據(jù)中震級較大（ML>3.0）的地震信號外，重疊在主震以及較大余震尾波中的微小余震信號也變得十分清晰，如在主震后600s時間窗內(nèi)有多組清晰的峰值曲線。如圖2c所示，當把臺站背景噪聲水平nij人為提高5dB時，可以看到記錄信號中的大部分余震淹沒在背景噪聲中，導致無法分辨。這表明記錄信號的時頻能量特征分布主要取決于記錄臺站在頻域內(nèi)的背景噪聲水平。事實上對于一個固定的地震臺站而言，一般情況下背景噪聲水平的源是相對穩(wěn)定的（王俊等，2013）。因此，在準確求得臺站的背景噪聲水平后，采用以背景噪聲水平基準的特征包絡曲線來反映記錄信號的特征是客觀、可靠的。

隨后，利用峰值搜索法對時頻包絡線進行觸發(fā)檢測，具體方法是：在設定的觸發(fā)能量閾值Fc（圖2b中黑色線條）線以上，根據(jù)峰值的持續(xù)時間（Tc）、坡度（Pc，即峰值的高度除以持續(xù)時間的1/2）閾值，進行峰值搜索①。圖2b中黑色箭頭表示在Fc=10Hz、Tc=20s、Pc=2時的觸發(fā)情況，除主震外，還檢測到29次觸發(fā)事件。將觸發(fā)結果與正式編目結果中寶應臺的Pg波震相到時進行對比，以檢驗觸發(fā)準確率，將觸發(fā)時間誤差在一定范圍內(nèi)的事件視為同一事件，本文中誤差取值為2s。由于波譜能量時頻包絡線法檢測到的位置是地震信號中能量分布最強的位置，對于地方震而言它對應的位置一般是S波的振幅最大位置處（圖2b，c）。以人工地震編目中Pg震相到時為基準，寶應臺實際記錄的S-P到時差約為5s，因此將到時誤差在4～6s范圍內(nèi)的事件視為同一事件。將準確觸發(fā)事件數(shù)與觸發(fā)總數(shù)之比，定義為準確率；未觸發(fā)的事件數(shù)與實際事件數(shù)之比，定義為漏觸發(fā)率。統(tǒng)計結果顯示：采用波譜能量包絡方法準確觸發(fā)的事件共有26次，準確率約為83.4%，漏觸發(fā)率約10.3%。相比而言，采用STA/LTA算法在觸發(fā)閾值為3時的常規(guī)狀態(tài)下，共檢測出23次，但只有16次觸發(fā)是有效的，準確觸發(fā)率約為69.5，漏觸發(fā)率達44.8%。詳細的對應關系見表1。

2.2 新疆和靜M6.6地震余震序列

為了進一步驗證本文方法，選取1組更為復雜的序列進行實測。2012年6月30日5時7分，在新疆和靜（43.4°N，84.8°E）發(fā)生M6.6地震，根據(jù)新疆地震臺網(wǎng)記錄結果，震后1h內(nèi)共記錄余震37次，距離震中最近的石場臺（SCH臺）震中距約為100km。石場臺背景噪聲的加速度功率譜概率密度分布如圖3所示，1～20Hz頻帶范圍內(nèi)的平均值約為-143.3dB（噪聲水平的有效值為6.82×10-8m/s）。圖4是震后1h內(nèi)石場臺的記錄波形和時頻能量分布結果。

同樣，對石場臺的時頻包絡特征曲線進行檢測，共檢測到39次觸發(fā)結果。石場臺記錄中的S-P到時差約為16s，將自動檢測到時與人工編目中Pg波震相到時差在15～17s范圍內(nèi)的事件視為同一事件，結果顯示：正確觸發(fā)的事件共有31次，正確率約83.8%，漏觸發(fā)率約為16.2%。而采用STA/LTA自動算法共檢測出觸發(fā)事件45次，但其中正確的觸發(fā)事件數(shù)僅有27次，正確率約為60.0%，漏觸發(fā)率約為27.0%。通過以上結果可以看出，本研究提出的檢測方法，在這2個指標上也明顯優(yōu)于STA/LTA算法的檢測效率。為了對這2種算法的運算速度進行比較，將2組記錄波形從0s處開始，以窗長為180s的間隔，漸進式增加記錄數(shù)據(jù)長度進行計算，2種算法實際計算耗時對比如圖5所示。對于處理3600s采樣率為100Hz的數(shù)據(jù)，在數(shù)據(jù)記錄長度小于2000s的階段，時頻包絡特征算法的計算耗時與STA/LTA算法基本一致。隨著數(shù)據(jù)長度的增加，時頻包絡算法的計算耗時的增加幅度逐漸略大于STA/LTA算法，這主要是由于隨著記錄數(shù)據(jù)長度的增加，計算時頻能量分布所需耗時逐漸增加，但這種影響在實時處理過程中可以忽略不計，因為通常情況下實時檢測的數(shù)據(jù)長度為s級或數(shù)十s，在這種情況下時頻包絡特征算法的計算耗時可以達到s級。

3 結論與討論

本文針對密集余震序列的記錄特征，研究了一種基于地震波波譜包絡分布特征的快速觸發(fā)檢測方法。在繼承聲譜圖像模式識別的正演邏輯觸發(fā)模式的基礎上，將其改進為以臺站固有的地噪聲水平為觸發(fā)的判據(jù)。改進后的算法，將構建二維灰度圖像轉換為譜能量分布圖，并且舍棄了原來需要構建不同信號之間的波形非線性圖像和復雜識別模式的過程，在初至震相的觸發(fā)效率上僅取決于臺站記錄頻帶內(nèi)背景噪聲固有的時頻特征，這樣不僅克服了由單一頻率、無規(guī)則的瞬時大振幅而引起的誤觸發(fā)，還很好地避開了STA/LTA算法中設置觸發(fā)閾值的問題。通過頻譜分析精確分離出記錄信號中與背景噪聲有差異的信號后，繼續(xù)沿用原方法中采用非線性標度縮放因子對大于臺站背景噪聲水平的能量進行歸一化處理的方式，顯著提高了記錄信號中弱信號的信噪比，使得特征包絡曲線的提取更加精確，因此在序列地震的觸發(fā)檢測中大大降低了主震尾波對后續(xù)地震檢測的影響。

2組余震序列的實測結果表明：本文所提出的方法，相對STA/LTA算法自動檢測結果，正確觸發(fā)率提高約19.6%，漏觸發(fā)率則降低了22.6%。這表明基于地震波波譜能量分布的檢測方法具有更強的抗干擾能力。由于本文方法是通過提取包絡特征曲線來進行初至震相的觸發(fā)檢測，因此具有與STA/LTA算法（即短時平均與長時平均值之比的曲線）相似的特點，即運行速度快，能夠滿足地震實時處理系統(tǒng)的技術要求。

初至震相的識別是建立在識別信號和噪聲差異的基礎之上，如振幅、頻率、偏振、功率譜等。由于記錄信號常常會受到很多干擾因素的影響，如震級大小、震中距、記錄臺站相對于震源的方位甚至記錄系統(tǒng)故障等。因此，沒有任何單一的方法可以保證初至震相的識別完全可靠；但從本文的實測結果來看，基于地震波波譜包絡分布特征的檢測方法，在密集余震序列的初至震相檢測中是具有明顯優(yōu)勢的。

參考文獻：

黃媛，吳建平，張?zhí)熘?2008.汶川8.0級大地震及其余震序列重定位研究[J].中國科學：地球科學，38（10）：1242-1240.

劉希強，周蕙蘭，沈萍，等.2000.用于三分向記錄震相識別的小波變換方法[J].地震學報，22（2）：125-131

王俊，徐戈，孫業(yè)君.2009.江蘇省區(qū)域地表背景噪聲特性的分析[J].地震研究，32（2）：155-161.

王俊，鄭定昌，詹小艷，等.2013.基于背景噪聲互相關格林函數(shù)的單臺時間誤差枯計[J].地震學報，35（6）：888-901

王勤彩，陳章立，鄭斯華.2009.汶川大地震余震序列震源機制的空間分段特征[J].科學通報，54（16）：2348-2354

朱艾斕，徐錫偉，刁桂荃，等.2008.汶川MS8.0地震部分余震重新定位及地震構造初步分析[J].地震地質(zhì)，30（3）：759-767.

Allen R.1978.Automatic earthquake recognition and timing from singletrace[J].BSSA，68（5）：1521-1532.

Bormadn P.2012.NewManual of Seismological Observatory Practice 2（NMSOP-2）Potsdam：Deutsches GeoForschungsZentrum GFZ，pp.1-65，doi：http：//doi.org/10.2312/GFZ.NMSOP-2_ch8

Joswig M.1990.Pattern recognition for earthquake detection[J].BSSA，80（1）：170-186.

Joswig M.1993.Single-trace detection and array-wide confidence asso-ciation of local earthquakes and explosions[J].Computers and Geo-sciences，19（2）：207-221.

Leonard M，Kennett B L N.1999.Multi-component autoregressive tech-niques for the analysis of seismograms[J].Phys Earth Planet Interi，113（1-4）：247-264.

Leonard M.2000.Comparison of manual and automatic onset time picking[J].BSSA，90（6）：1384-1390.

Massa M，F(xiàn)erretti G，Spallarossa D.2006.Improving automatic locationprocedure by waveform similarity analysis：An application in theSouth Western Alps（Italy）[J].Phys Earth Planet Interi，154（1）：18-29.

McNamam D E，Boaz R I.2005.Seismic Noise Analysis System UsingPower Spectral Density Probability Density Functions：A Stand-A-lone Software Package[R].Virginia U S Geological Survey.Patane D，F(xiàn)errari F，F(xiàn)ernicci F.1999.First application of ASDP software：A case study at Mt.Etna volcano and in the Acri region（Southern It-aly）[J].Phys Earth Planet lnteri，113（1-4）：75-88.

Sleeman R，Van E，Torild.1999.Robust automatic P-phase picking：Anon-line implementation in the analysis of broadband seismogram re-cordings[J].Phys Earth Planet Interi，113（1-4）：265-275.

Stevenson P R.1976.Microearthquakes Flathead Lake，Montana：A studyusing automatic earthquake processing[J].BSSA，66（1）：61-79.

Wang J，Teng T L.1995.Artificial neural network-based seismic detec-tor[J].BSSA，85（1）：308-319.

Withers M，Aster R，Young C，Beiriger J，et at，1998.A comparison of se-lect trigger algorithMS for automated global seismic phase and eventdetection[J].BSSA，88（1）：95-106.