利用背景噪聲研究滇西北波速異常與M≥5.0地震關系

楊建文 葉 泵 高 瓊 陳 佳 王 軍

1 中國地震科學實驗場大理中心，云南省大理市濱海大道，6710002 云南大理滇西北地殼構造活動野外科學觀測研究站，云南省大理市濱海大道，671000

開展地下介質波速變化測量研究對于認識地震孕育機理和探索地震物理預報手段具有重要意義。隨著背景噪聲層析成像技術的發展，經驗格林函數也被應用于地殼內部地震波速度的時空變化特征研究[1]。背景噪聲分布廣且不依賴于特定震源，同時具有可重復性和經濟環保等優勢，非常適合跟蹤地殼內部結構的物性變化[1]。

背景噪聲較天然震源具有更高的時間采樣率，也比人工震源更加經濟。通過背景噪聲波速方法，地震學家觀測到火山區波速的季節性變化[2]，發現火山噴發前數周內波速具有下降特征[3]。部分學者[4-10]對同震和震后地震波速度的變化特征進行了研究，但大部分都基于波形資料來反演強震發生過程及震后的波速變化。

近年來，滇西北地區構造運動強烈，震情形勢嚴峻。本文以滇西北作為實驗地點，基于部分固定臺站記錄的寬頻帶連續波形資料，采用背景噪聲互相關方法研究地殼介質內部地震波速隨時間的變化特征，通過對臺站對走時偏移時間序列進行深加工，提取適合于滇西北地區M≥5.0地震的短臨異常識別指標。該工作的開展對深入研究滇西北地區地下介質的變化過程，探索測震資料及背景噪聲技術在日常震情跟蹤及地震預測中的應用，提高異常識別能力及推進指標體系建設等具有重要意義。

1 數據收集與整理

本文主要研究滇西北地區(25°～28°N，99°～101.5°E)波速變化，選取云南數字地震臺網中5個固定臺站的垂直分量連續波形資料，時間范圍為2012-01～2020-11，該時間段內所有臺站的數據資料連續性都很好，數據基本能夠反映介質波速的變化。研究時段內，研究區共發生6次M≥5.0地震，圖1為5個地震臺站、2個GNSS連續觀測站及6次地震震中分布情況，相關地震參數見表1。

表1 滇西北6次M ≥5.0地震相關參數

圖1 所用地震臺站、GNSS連續觀測站、地震震中分布Fig.1 Distribution of seismic stations, GNSS continuousobservation stations and earthquake epicenter

2 背景噪聲波速測量方法

2.1 單臺站數據預處理

采用滇西北地區5個寬頻帶地震臺記錄的垂直分量連續波形數據計算地震波速度。首先將連續波形數據按天截取[10]，采用1 Hz對波形數據進行抽樣，以降低數據采樣率；然后進行去除儀器響應、去均值、去趨勢、帶通濾波處理[1,8-10]；接著對波形數據進行時域平均，以去除地震事件的影響；最后進行頻譜白化處理[8,10]。以上數據處理的目的是為了消除天然地震和儀器本身異常信號的影響，以獲取高質量的地震背景噪聲。

2.2 背景噪聲互相關及數據疊加

單臺數據預處理后，對由5個臺站構成的10個臺站對每天的背景噪聲進行互相關處理，以獲取各臺站對單日的經驗格林函數[8]。Stehly等[11]研究認為，對任意2個臺站記錄的背景噪聲進行互相關處理，理論上得到的經驗格林函數中有一正、一負2個分支，分別表示臺站對路徑上的因果和非因果信號。但背景噪聲的來源及能量差異，會極大程度地影響經驗格林函數的形態[12]。當臺站兩側的噪聲源分布均勻時，因果信號和非因果信號的到時一致，振幅相同；而當噪聲源分布不均勻時，2個方向信號的到時相同，但振幅不同，在噪聲源能量較強一側產生的信號振幅較大[11-12](圖2)。

圖2 TUS-YOS臺站對經驗格林函數Fig.2 Empirical Green’s function of TUS-YOSstation pair

Bensen等[13]研究認為，對因果和非因果信號進行疊加處理可以等效噪聲源的均勻分布。雖然該處理會損失振幅信息，但可以在一定程度上提高信號的信噪比[12-13]。因此，本文以0為界線，將經驗格林函數的“負支”反號后與“正支”進行對稱疊加，形成“疊加分量”(圖3)。

圖3 TUS-YOS臺站對對稱疊加后的經驗格林函數Fig.3 Empirical Green’s function after symmetricsuperposition of TUS-YOS station pair

由于受噪聲源的不均勻分布和噪聲成分的不確定性等因素影響，同時臺站間距離較遠，導致臺站對單日經驗格林函數的信噪比普遍較低，難以分辨地震波形態[10]。為提高信噪比，本文將各臺站對單日前59 d(共60 d)的經驗格林函數進行疊加，作為當天經驗格林函數(圖4(b))。在此基礎上，通過進一步滑動疊加，獲取2012-02-29～2020-11-16當天經驗格林函數(前59 d因疊加數據不足60 d，不作進一步分析)。

為獲得臺站對經驗格林函數隨時間變化的定量信息，需確定各臺站對的參考經驗格林函數[8-10]。Stehly 等[11]研究發現，背景噪聲源的季節性變化會引起由互相關方法重構的面波信號走時變化，因此建議盡量采用較長時間范圍(至少1 a)的資料來確定參考經驗格林函數。本文將各臺站對2012-02-29～2020-11-16共3 184 d的單日經驗格林函數進行疊加，作為各自的參考經驗格林函數(圖4(a))。

陰影部分表示臺站對走時偏移計算中所使用的直達瑞利波信號圖4 TUS-YOS臺站對當天經驗格林函數與參考經驗格林函數互相關計算結果Fig.4 Cross-correlation calculation results of the day empirical Green’s function and reference empirical Green’s function of TUS-YOS station pair

2.3 走時偏移提取

在獲取臺站對當天經驗格林函數和參考經驗格林函數后，對當天經驗格林函數和參考經驗格林函數進行互相關計算[8,10]，獲取相似波形窗口的互相關系數(圖4(c))，互相關系數最大時對應的時間延遲就是這兩段相似波形窗口的走時偏移。對計算結果進行余弦插值以提取更高精度的走時偏移。以互相關系數為約束(互相關系數要求在0.9以上)，通過提取直達瑞利波中部分震相的走時偏移來計算臺站對之間瑞利波的速度。圖5(a)為TUS-YOS臺站對瑞利波走時偏移變化曲線。

圖5 TUS-YOS臺站對瑞利波走時偏移變化曲線Fig.5 The variation curve of Rayleigh wave travel timeoffset of TUS-YOS station pair

從圖5(a)可以看出，TUS-YOS臺站對瑞利波走時偏移變化曲線存在明顯的年變成分。為消除年變成分對走時偏移的影響，采用傅立葉變換以獲取長趨勢變化背景下隱含的波速變化特征。對于去除年變成分后的走時偏移，劃定±1.5倍標準差作為單條曲線的異常指標(以±1.5倍標準差作為異常閾值)。從圖5(c)可以看出，該臺站對的走時偏移在洱源M5.5和M5.0、云龍M5.0、漾濞M5.1等地震前存在較為明顯的超指標異常。

2.4 臺站對走時偏移可靠性分析

由圖1可知，地震臺站YUL和TUS附近分別有云龍和下關2個GNSS連續觀測站，且TUS-YUL臺站對的路徑分布與云龍-下關GNSS基線基本保持平行。圖6為TUS-YUL臺站對瑞利波走時偏移與云龍-下關GNSS基線相對變化的對應關系(GNSS數據獲取時間截至2017-05-20)。從趨勢上來看，TUS-YUL臺站對瑞利波走時偏移與云龍-下關GNSS基線相對變化存在較好的耦合性，即在一段時間內，當TUS-YUL臺站對瑞利波走時偏移縮短(波速加快)時，云龍-下關GNSS基線長度也在縮短(擠壓增強)。由于本文提取的數據為臺站對直達瑞利波走時偏移，反映的是中上地殼物性隨時間的變化特征，而GNSS基線監測的是中上地殼的應力變化，兩者具有較好的耦合性，表明本文提取的瑞利波走時偏移能夠監測到中上地殼介質的物性變化，這也間接表明本文的數據結果具有可靠性。

圖6 TUS-YUL臺站對瑞利波走時偏移與云龍-下關GNSS基線相對變化對應關系Fig.6 Correspondence between the Rayleigh wave travel time offset of TUS-YUL station pair and therelative change of Yunlong-Xiaguan GNSS baseline

3 短臨異常指標提取及效能評價

由于需要獲取短臨異常信息，地震事件對應的時間間隔應不超過90 d。對10個臺站對的瑞利波走時偏移時間序列進行傅立葉滑動去年變處理，并設定±1.5倍標準差作為異常閾值，結合研究區6次地震對各臺站對走時偏移時間序列的映震能力進行效能評價。現階段應用最為廣泛的預報效能評價方法為R值評分法[14]，其主要依據預測指標實際預報地震的有震報準率c和預報占時率(預報時間占有率)b[15]，可反映預測方法或指標與地震的相關程度。R值計算公式[15]為：

R=c-b=有震報準率-預報占時率=

(1)

式中，R=1表示全報準；R=0表示預報未起作用。R值越大，預報效果越好。根據R值評分規則，在97.5%置信度下，當預報效能R>R0(預報效能檢驗的最低值)時[16]，表明所評估的預測方法可通過統計檢驗，高于隨機預測的預報效能。

表2為10個臺站對瑞利波走時偏移時間序列的預報效能檢驗結果。

表2 各臺站對預報效能檢驗結果

由表2可知，在10個臺站對中，有7個臺站對通過了檢驗，3個臺站對未通過檢驗，但通過檢驗的7個臺站對的對應率普遍偏低(各臺站對均存在較高的虛報)。考慮到單個臺站對的虛報率偏高，且可能存在不確定性等因素，需要將通過檢驗的7個臺站對進行綜合分析，以獲取穩定可靠的地震短臨異常識別指標(綜合指標)。本文采用自適應加權綜合預測方法獲取綜合概率指標，計算公式為[17]：

(2)

式中，n為預報因子總數；Kti為到時間t為止的第i項因子動態值(當第i項因子出現異常后，其所對應的地震危險時段K=1，其余時段K=0)；Pti為到時間t為止的第i項因子歷史映震率；Nti為到時間t為止的第i項因子歷史映震次數；Pt為時間t的綜合概率值[17-18]。

從理論和實際預測中均可發現，在各項因子與地震的相關性及預測貢獻中，最重要的表征值為映震率Pti。另外，各因子的映震次數是因子與地震相關性置信度評價最主要的因素，通過置信度檢驗的相關因子，映震次數Nti越大，可信度就越高[18]。

圖7為采用自適應加權綜合預測方法對通過檢驗的7個臺站對進行綜合分析后獲取的地震短臨異常識別指標(綜合指標)。

圖7 滇西北地區地震短臨異常識別指標Fig.7 The identification index of short-term andimminent earthquake anomaly in northwestern Yunnan

在獲取滇西北地區的綜合指標后，設定2倍均值線作為異常閾值(此處為0.60)，當綜合概率值Pt超過2倍均值線時，表示出現異常，認為滇西北地區未來90 d內存在發生M≥5.0地震的危險。

在提取滇西北地區的綜合指標和劃定2倍均值作為異常閾值后，需要結合對應的地震事件對綜合指標進行效能評價(效能評價采用R值評分法)。預報效能檢驗結果表明，利用該綜合指標對滇西北地區2012年以來發生的6次M≥5.0地震進行90 d短臨預報時，異常指標共出現8次(2016年異常指標超過90 d，故劃分為2次)，其中準確預報地震5次，漏報1次，虛報4次，預報效能評分R為0.692，R0為0.475，R>R0，效能檢驗通過。該綜合指標的地震對應率為62.50%，概括率為83.33%，表3為綜合指標異常出現時間與地震發生時間的對應情況。

表3 綜合指標異常出現與地震發生對應情況

由表3可知，大部分地震的發生與異常指標的出現時間間隔均不超過45 d，表明該綜合指標在時間上具有較好的短臨預報能力。結合圖1中臺站對路徑及地震震中分布可知，該綜合指標對臺站對覆蓋范圍內的洱源M5.5、洱源M5.0、云龍M5.0和漾濞M5.1地震都作出了準確的短臨預報，不存在漏報現象，從而檢驗了該指標的預測能力。此外，該指標對臺站對范圍之外的昌寧M5.0地震也可作出一定預報，表明該指標對臺站對覆蓋范圍及鄰區的地震均有較好的預報能力。

4 討 論

通過背景噪聲互相關計算可獲取臺站對的經驗格林函數。研究表明，噪聲源往往呈不均勻分布，且具有一定的季節性變化，這與理論假設存在差異[12]。該差異使得通過互相關計算獲取的經驗格林函數并不完全收斂于臺站間的理論格林函數，并存在季節性變化，從而導致獲取的走時偏移(波速變化)受季節變化影響存在明顯的年變特性。Wang等[19]基于日本臺站數據對地震波的季節性變化特征進行研究發現，地震波速度存在夏季低、冬季高的特點。根據波速和走時偏移的計算關系[3]dv/v=-dt/t可知，本文計算得到的瑞利波走時偏移dt存在夏季高、冬季低(對應波速dv為夏季低、冬季高)的變化特征，與Wang等[19]的結果一致。為消除季節性變化對走時偏移的影響，本文對走時偏移時間序列進行傅立葉滑動去年變處理，該方法理論基礎嚴密，是一種線性、無偏、最優的周期估計方法，利用該方法去除年變成分的可行性已得到認可[20]。傅立葉滑動去年變方法能最大限度地去除年變特性對走時偏移的影響，因此認為本文提取的短臨異常識別指標(綜合指標)具有可靠性。

本文基于6個震例，采用R值評分法對臺站對走時偏移時間序列的映震能力進行檢驗發現，有7個臺站對通過檢驗，3個未通過檢驗。究其原因，可能與強震發生區域的波速變化存在一定的地域性，且不同臺站對之間的面波傳播路徑存在差異有關。具體原因還有待進一步研究。

2021-05-21在本文臺站對覆蓋范圍內發生漾濞6.4級地震，為驗證本文提取的短臨異常識別指標(綜合指標)是否對該地震的發震時刻具有指示意義，對綜合指標進行數據更新，結果如圖8所示。

圖8 更新后的滇西北地震短臨異常識別指標Fig.8 The updated identification index of short-term and imminent earthquake anomaly in northwest Yunnan

由圖8可知，該短臨異常識別指標(綜合指標)在2021-04-20(距漾濞6.4級地震發震時刻32 d)出現明顯的超指標異常，因此漾濞6.4級地震的發生是對本文提取的短臨異常識別指標(綜合指標)的很好檢驗。

在上述分析的基礎上，對2021-04-20異常指標出現時的異常臺站對進行提取，以探討異常臺站對的分布與漾濞6.4級地震震中位置的對應關系，結果如圖9所示。從圖中可以看出，2021-04-20綜合指標出現異常時，震中附近DAY-TUS、DAY-YUL和YUL-ZOD臺站對的走時偏移出現超指標異常，漾濞6.4級地震發生在3個異常臺站對的附近區域。由此可見，基于綜合指標提取的異常臺站對分布對未來地震的發震地點具有一定的指示意義。

圖9 漾濞6.4級地震震中與異常臺站對分布對應情況Fig.9 Correspondence of the epicenter distribution of Yangbi MS6.4 earthquake and abnormal station pairs

5 結 語

本文基于地震臺站記錄的寬頻帶連續波形資料，采用背景噪聲互相關及傅立葉變換等方法，獲取臺站對瑞利波走時偏移。通過對臺站對走時偏移時間序列進行深加工處理，進而提取出適合于滇西北地區的地震短臨異常識別指標(綜合指標)。結果表明，該指標物理意義明確，具有對背景噪聲波速測量結果高度綜合的特點，可提高對異常的識別能力，對滇西北地區M≥5.0地震的發震時間具有較好的指示意義。但有幾點需要注意：

1)在對每個臺站對瑞利波走時偏移時間序列進行效能評價時，本文選取±1.5倍標準差作為異常閾值，主要是為提高信噪比，閾值設置過低會帶入部分噪聲。

2)在綜合異常指標中，選取2倍均值線作為異常閾值，大部分指標在超過閾值指標線后較短時間內發生地震，因此在后續預測實踐中，及時處理數據尤為重要。另外，該閾值可根據實際需要進行動態調整。

3)本文研究是將背景噪聲波速測量結果應用于地震短臨預報的一次很好的嘗試。研究發現，單個臺站對的瑞利波走時偏移時間序列對整個云南地區M≥5.0地震的發震時刻均有較好的指示意義，因此本文方法可拓展至整個云南地區。基于云南地區固定臺站的連續波形資料，提取適合于云南地區M≥5.0地震的短臨異常識別指標，對整個云南地區的震情跟蹤及預測預報工作都具有重要意義。

致謝:本文研究所用連續波形資料來源于云南數字地震臺網，GNSS基線數據來源于云南省地震局信息中心，背景噪聲互相關計算程序由中國地震局地球物理研究所地球物理先導技術研究室王偉濤研究員提供，綜合指標基于云南省保山市地震局李宗興高級工程師開發的前兆處理軟件計算得到，在此表示衷心感謝。