馬騰飛

1) 中國北京100081中國地震局地球物理研究所 2) 中國北京100033中國財產(chǎn)再保險有限責(zé)任公司

?

三分量地震記錄的互相關(guān)分析

1) 中國北京100081中國地震局地球物理研究所 2) 中國北京100033中國財產(chǎn)再保險有限責(zé)任公司

地震波場本質(zhì)上是三維矢量波場， 單分量記錄實際上是三維矢量震動在某一方向上的部分投影． 本文基于單分量地震記錄互相關(guān)公式， 提出了一種新的適合三分量地震波形記錄多元互相關(guān)運算的簡易方法， 并以2008年汶川MS8.0地震余震序列波形為例， 對其進行了效果驗證． 結(jié)果表明， 相對單分量互相關(guān)， 該方法可以得到更為合理的全局最優(yōu)結(jié)果， 解決波形識別匹配工作中不同分量間的差異問題． 該方法還可以利用不同分量間的“同源”信息， 有效壓制隨機噪聲， 并從理論上說明其合理性． 其原理及計算過程均較為簡單， 整體運算量較小， 適用于目前“大數(shù)據(jù)”時代的地震數(shù)據(jù)處理．

三分量記錄 多元互相關(guān) 模板匹配方法

引言

波形互相關(guān)技術(shù)是地震學(xué)中一種常用的技術(shù)手段， 目前已經(jīng)在“重復(fù)地震”識別(Schaff， Richards， 2004， 2011; Lietal， 2007， 2011; Maetal， 2014)、 余震事件檢測(Peng， Zhao， 2009; Wuetal， 2014)、 低頻事件觀測與識別(Obara， 2002; Shellyetal， 2007)、 地震精定位(Waldhauser， Ellsworth， 2000; Schaffetal， 2004; Schaff， Waldhauser， 2005)等領(lǐng)域得到了廣泛應(yīng)用． 迄今為止， 地震學(xué)家已經(jīng)發(fā)展了許多先進的數(shù)字記錄處理方法來計算兩個波列間的互相關(guān)系數(shù)， 但這些方法絕大部分是針對單分量的時間序列記錄． 鑒于地震記錄的三維空間屬性， 各個分量僅是質(zhì)點運動在垂直或水平方向的投影， 因此基于各個分量的單分量互相關(guān)系數(shù)計算可能會丟失地震記錄中某些空間相關(guān)信息， 不能反映地震波在傳播過程中波形和震相變化的全貌， 只得到片面性的結(jié)果． 事實上， 只有三分量耦合的空間記錄才能真實地反映實際地震波場所包含的全方位信息．

就地震事件波形而言， 由于其不同震相在不同分量上有較大的運動學(xué)差異， 因此在采用互相關(guān)識別計算時對所用分量均采用統(tǒng)一時間窗口就顯得不大“合適”． 例如： 由于地震信號中噪聲的存在， 當(dāng)初至P波到達時， 垂直分量會產(chǎn)生較明顯的變化， 但水平分量依舊處于信噪比很低的“噪聲模式”； 在S波尾部， 垂向信噪比較低， 水平向卻還有較為明顯的振動． 在這些互相關(guān)運算的時間窗口內(nèi)， 低信噪比“噪聲”的存在無疑會對最終的運算結(jié)果產(chǎn)生較大影響， 使事件的識別與檢測面臨較大困難， 對于震級較小的微震事件來說更是如此． 針對這些問題， 通過采用對不同分量的地震記錄選取不同相關(guān)運算窗口的方法(如垂直分量時間窗口為P波到達后4 s， 水平分量時間窗口為S波到達后4 s)， 便可在一定程度上提高檢測識別的準(zhǔn)確率(Peng， Zhao， 2009; Mengetal， 2012; Wuetal， 2014)． 但是， 這種硬性規(guī)定的不同分量時間窗口難免會“錯殺”一部分不符合這種“標(biāo)準(zhǔn)模式”的地震事件， 從而影響其識別的完整性； 與此同時， 當(dāng)所獲取的各個分量之間的差異性較大時， 如何對所得結(jié)果作出合理的解讀也是一大難題． 對于臺陣記錄我們可以通過各種技術(shù)手段疊加不同臺站的信息來達到壓制噪聲、 提高信噪比的目的(Leonard， Kennett， 1999; Kennett， 2000)， 但對于單臺站地震記錄則無法開展． 因此， 如何充分發(fā)掘不同分量間的“同源”作用， 得到能全面反映地震波三維屬性的相關(guān)信息也是本文將要探討的內(nèi)容．

1 基本方法與理論

1.1 單分量波形互相關(guān)原理

(1)

1.2 三分量多元綜合互相關(guān)系數(shù)

三分量地震記錄有3個獨立分量(垂直分量V， 切向分量T， 徑向分量R)， 不同波列組合后可以形成一個3×n的矩陣. 對于多維矢量矩陣而言， 空間夾角沒有意義， 這種情況下則不能用上述向量相關(guān)的思路來解決多元的相關(guān)問題． 在實際應(yīng)用中， 通常對不同分量兩兩相關(guān)后計算得出3個獨立的相關(guān)系數(shù)， 該相關(guān)系數(shù)矩陣可以用來描述兩矩陣間的相關(guān)關(guān)系． 但在具體工作中我們也會遇到諸如不同分量間最大相關(guān)位置不一致、 各分量間相關(guān)系數(shù)差別較大等問題， 這給我們帶來了較大的挑戰(zhàn)． 事實上， 由于各分量間的的振幅能量、 信噪比(signal noise ratio， 簡寫為SNR)水平均不相同， 上述問題的出現(xiàn)在匹配識別工作中并不罕見． 如果可以找到一種簡單快捷的方法， 能夠綜合考慮各分量的振幅能量水平， 得到全局最優(yōu)結(jié)果， 無疑對此類工作的開展具有重要意義．

一種可行的辦法為矩陣向量化， 即將三分量記錄投影展開到一條直線上， 以便我們能繼續(xù)使用向量相關(guān)的計算公式來處理三分量問題． 需要注意的是， 所采用的變換方式必須使各分量之間滿足等價互易性， 否則所得結(jié)果不唯一．

圖1 三分量地震記錄投影到直線上的示意圖 Fig.1 The diagram of three-component recordings projected onto a line

(2)

其中，

(3)

(4)

(5)

需要注意的是， 式(5)中fi(t)和gi(t)均為歸一化前記錄到的原始數(shù)據(jù)， 其中包含各個分量的絕對振幅信息． 由此可以看出： 這種“拼接”處理的實質(zhì)在于可以將各個分量間不同的振幅及相關(guān)信息置于同一參考系下， 從而得出考慮全局后的整體結(jié)果； 同時將各分量波形置于更大參考系下也可以有效壓制振幅較小、 相位不相關(guān)的噪聲部分， 增大綜合信噪比， 從而提高識別精度． 從形式上看， 式(5)與單分量相關(guān)公式很接近， 且其原理技術(shù)相對簡單， 形式也較為簡潔， 適合大規(guī)模地震數(shù)據(jù)資料的處理計算．

2 應(yīng)用實例

以成都臺(CD2)記錄到的2008年汶川MS8.0地震余震序列中的兩個地震事件波形為例， 詳細分析所得三分量整體相關(guān)公式在實際中的應(yīng)用效果． 該地震事件對的震中距為27.385 km， 為典型的近臺記錄， 目錄參數(shù)引自中國地震臺網(wǎng)中心的《中國地震月報目錄》， 兩個地震事件均屬于微震事件， 震級幾乎相等， 波形數(shù)據(jù)引自中國地震局地球物理研究所“國家數(shù)字測震臺網(wǎng)數(shù)據(jù)備份中心”(Zhengetal， 2010)．

由于成都臺采用甚寬頻地震計， 在較大范圍內(nèi)的頻率響應(yīng)曲線較為平緩， 因此我們采用1—10 Hz四階巴特沃斯(Butterworth)帶通濾波器對去除均值、 線性趨勢后的原始波形進行處理， 這也與前人工作的參數(shù)選取相一致(Lietal， 2007， 2011; Maetal， 2014)． 以P波到時為互相關(guān)計算起點， 所選取的窗口時長為P波與S波走時差的4倍， 滑動時長為±2 s， 這樣便可包含全部的尾波序列， 同時避免后續(xù)噪聲混入影響計算結(jié)果(蔣長勝等， 2008)．

采用上述流程， 計算所得的三分量波形互相關(guān)系數(shù)分別為0.81205(E--W)， 0.86067(N--S)和0.76896(U--D)， 如圖2所示． 可以看出， 各分量的波形只在某些不同時段上具有較高的相似性． 由于各分量間的計算結(jié)果差異較大， 對于兩次地震的判定也成為一大難題．

圖3給出了兩次地震在成都臺(CD2)的三分量整體互相關(guān)系數(shù)以及各分量間的互相關(guān)系數(shù)隨時間的變化曲線． 取P波到時前1 s起算， 以窗口中心為計算點， 2 s為互相關(guān)計算滑動窗口長度， 4倍的S-P走時差為窗口長度， 自左向右， 每次移動1個數(shù)據(jù)點(即0.01 s)進行互相關(guān)運算． 為了避免地震波的諧波特性對此處小窗口計算結(jié)果產(chǎn)生較大的影響， 每次進行窗口計算時， 波形只能相對移動±0.01 s (1個數(shù)據(jù)點)． 由圖3可以看出， 南北分量與東西分量在全程運算中均有較高的互相關(guān)值(≥0.8)， 而垂向分量在P波和S波到達時間之外振幅較小部分的互相關(guān)值有較大的波動， 拖累了整體的相關(guān)系數(shù)計算， 因此未達到給定的閾值(≥0.8)． 此外， 不同時段的總體相關(guān)系數(shù)給予不同分量的“權(quán)重”不同， 振幅越大的分量權(quán)重越高， 因此采用三分量相關(guān)可以部分壓制不同分量內(nèi)部低信噪比的“噪聲”部分， 提高識別的效率和準(zhǔn)確率． 在上述例子中， 運用三分量整體相關(guān)公式后可以得出兩個地震事件的互相關(guān)系數(shù)為0.8157， 滿足重復(fù)地震識別互相關(guān)系數(shù)≥0.8的閾值條件，可以視為一對“重復(fù)地震”事件．

圖2 成都(CD2)臺記錄到的2009-07-08(紅色)與2010-09-12(藍色)兩次地震濾波后的歸一化三分量波形對比圖及其互相關(guān)系數(shù)， 圖中綠色豎線表示P波到達時間

圖3 成都(CD2)臺不同分量互相關(guān)系數(shù)隨時間的變化曲線

圖4 不同分量計算得到的互相關(guān)系數(shù)與地震事件三分量綜合后和水平(或垂直)分量的歸一化互相關(guān)系數(shù)

與上述例子類似， 本文詳細統(tǒng)計了成都臺記錄到的一系列不同地震事件之間不同分量的互相關(guān)系數(shù)差異， 得到了采用本文方法后所得到的整體互相關(guān)系數(shù)與各分量互相關(guān)系數(shù)之間的關(guān)系， 如圖4所示． 可以看出， 各分量的相關(guān)系數(shù)總體上與全局相關(guān)系數(shù)呈線性對應(yīng)關(guān)系， 但僅看某一分量有時會出現(xiàn)較大偏離， 因此在這種條件下有必要根據(jù)三分量整體相關(guān)系數(shù)對事件進行合理判斷．

從圖4中也可以看出， 即便是對于某一給定的地震事件而言， 3個分量相關(guān)計算的結(jié)果也會呈現(xiàn)出一定的規(guī)律性， 即水平分量的相關(guān)系數(shù)顯著高于垂直分量， 東西分量的相關(guān)系數(shù)明顯高于南北分量， 因此在只能選擇一個分量作相關(guān)計算時， 垂直分量具有更高的識別可信度． 李宇彤(2012)利用區(qū)域臺網(wǎng)對海城—岫巖地區(qū)“重復(fù)地震”識別的研究也得出類似的結(jié)論， 這也說明了僅用垂直向的波形數(shù)據(jù)進行相關(guān)運算的合理性(Schaff， Richards， 2004， 2011)．

3 討論與結(jié)論

地震波本身為矢量場， 本質(zhì)上為不同特性、 不同類型的振動相互疊加干涉的結(jié)果， 而單分量記錄實際上僅為三維矢量在某一方向上的部分投影， 因此常規(guī)的基于各分量的單分量互相關(guān)計算可能會丟失信號中部分與空間相關(guān)的信息， 不能反映地震波在傳播過程中波形和震相變化的全貌， 從而導(dǎo)致結(jié)果不一致． 本文基于三分量記錄之間的“同源”特性， 提出了一種可以計算三分量記錄總體相關(guān)系數(shù)的簡易方法， 能夠盡可能地利用數(shù)據(jù)的內(nèi)在信息， 壓制隨機噪聲， 并從理論上說明了其合理性． 該方法的原理和計算過程均較為簡單， 整體運算量也較小， 無論從經(jīng)濟上還是技術(shù)上都適用于未來“大數(shù)據(jù)”時代的海量資料處理， 值得在實際工作中推廣應(yīng)用．

由文中的實例可以看出， 采用三分量整體相關(guān)可以解決互相關(guān)系數(shù)在不同分量間的差異以及臨界識別等問題， 所得結(jié)果也不是簡單的三分量單獨相關(guān)運算后的算術(shù)平均值， 而是在各分量間(inter-component)相關(guān)后綜合疊加得到的結(jié)果． 這實際上是一種對信號的壓噪重構(gòu)， 可以達到增強有效信號、 壓制干擾噪聲的目的， 解決了目前“重復(fù)地震”以及類似事件識別工作中遇到的問題． 此外， 采用三分量整體相關(guān)還可以同步三分量地震波形記錄， 避免片面追求各個分量的互相關(guān)系數(shù)單獨最大而造成錯誤時移(這種現(xiàn)象可能是由地震波的諧波特性所導(dǎo)致)， 因此從理論上來講， 也可能會存在整體相關(guān)系數(shù)比3個分量都小的極端情況， 但在實際中由于各種震相混疊、 介質(zhì)不均性等情況的客觀存在， 故難以出現(xiàn)上述情形．

應(yīng)該看到， 這種綜合相關(guān)算法似乎對在有效震相外信噪較低的波形記錄部分的相關(guān)計算效果并不明顯， 因此如何壓制有效震相外的信號噪聲以增加信噪比， 以及充分利用各分量中所包含的地震信息以提高微小地震事件的可探測性也是未來工作的一個可行方向．

吳忠良研究員為本研究進行了分析和指導(dǎo)， 與加州大學(xué)圣克魯茲分校地震學(xué)實驗室的Emily Brodsky教授、 Lian Xue博士、 Stephen Hernandez博士進行了有益討論， 中國地震局地球物理研究所“國家數(shù)字測震臺網(wǎng)數(shù)據(jù)備份中心” (doi:10.7914/SN/CB)為本研究提供了波形數(shù)據(jù)， 審稿專家提出了建設(shè)性的修改意見， 作者在此一并表示誠摯謝意．

蔣長勝， 吳忠良， 李宇彤. 2008. 首都圈地區(qū)“重復(fù)地震”及其在區(qū)域地震臺網(wǎng)定位精度評價中的應(yīng)用[J]. 地球物理學(xué)報， 51(3): 817--827.

Jiang C S， Wu Z L， Li Y T. 2008. Estimating the location accuracy of the Beijing Capital Digital Seismograph Network using repeating events[J].ChineseJournalofGeophysics， 51(3): 817--827 (in Chinese).

李宇彤. 2012. “重復(fù)地震”的若干地震學(xué)問題[D]. 北京: 中國地震局地球物理研究所: 36--40.

Li Y T. 2012.TheSeismologyof‘RepeatingEarthquakes’[D]. Beijing: Institute of Geophysics， China Earthquake Administration: 36--40 (in Chinese).

B?th M. 1974.SpectralAnalysisinGeophysics[M]. Amsterdam: Elsevier Scientific Publishing Company: 87--94.

Fredricksen H， Kessler I. 1977. Lexicographic compositions and deBruijn sequences[J].JCombTheory：SerA， 22(1): 17--30.

Kennett B L N. 2000. Stacking three-component seismograms[J].GeophysJInt， 141(1): 263--269.

Leonard M， Kennett B L N. 1999. Multi-component autoregressive techniques for the analysis of seismograms[J].PhysEarthPlanetInt， 113(1/2/3/4): 247--263.

Li L， Chen Q F， Cheng X， Niu F L. 2007. Spatial clustering and repeating of seismic events observed along the 1976 Tangshan fault， North China[J].GeophysResLett， 34(23): L23309. doi:10.1029/2007GL031594.

Li L， Chen Q F， Niu F L， Su J. 2011. Deep slip rates along the Longmen Shan fault zone estimated from repeating microearthquakes[J].JGeophysRes， 116(B9): B09310. doi:10.1029/2011JB008406.

Ma X J， Wu Z L， Jiang C S. 2014. ‘Repeating earthquakes’associated with the WFSD-1 drilling site[J].Tectonophy-sics， 619/620: 44--50.

Meng X F， Yu X， Peng Z G， Hong B. 2012. Detecting earthquakes around Salton Sea following the 2010MW7.2 El Mayor-Cucapah earthquake using GPU parallel computing[J].ProcCompSci， 9: 937--946.

Obara K. 2002. Nonvolcanic deep tremor associated with subduction in Southwest Japan[J].Science， 296(5573): 1679--1681. doi:10.1126/science.1070378.

Peng Z G， Zhao P. 2009. Migration of early aftershocks following the 2004 Parkfield earthquake[J].NatGeosci， 2(12): 877--881.

Schaff D P， Bokelmann G H R， Ellsworth W L， Zanzerkia E， Waldhauser F， Beroza G C. 2004. Optimizing correlation techniques for improved earthquake location[J].BullSeismolSocAm， 94(2): 705--721.

Schaff D P， Richards P G. 2004. Repeating seismic events in China[J].Science， 303(5661): 1176--1178.

Schaff D P， Waldhauser F. 2005. Waveform cross-correlation-based differential travel-time measurements at the Northern California Seismic Network[J].BullSeismolSocAm， 95(6): 2446--2461.

Schaff D P， Richards P G. 2011. On finding and using repeating seismic events in and near China[J].JGeophysRes， 116(B3): B03309. doi:10.1029/2010JB007895.

Shelly D R， Beroza G C， Ide S. 2007. Non-volcanic tremor and low-frequency earthquake swarms[J].Nature， 446(7133): 305--307.

Waldhauser F， Ellsworth W L. 2000. A double-difference earthquake location algorithm: Method and application to the northern Hayward fault， California[J].BullSeismolSocAm， 90(6): 1353--1368.

Wu C Q， Meng X F， Peng Z G， Ben-Zion Y. 2014. Lack of spatiotemporal localization of foreshocks before the 1999MW7.1 Düzce， Turkey， earthquake[J].BullSeismolSocAm， 104(1): 560--566.

Zheng X F， Yao Z X， Liang J H， Zheng J. 2010. The role played and opportunities provided by IGP DMC of China National Seismic Network in Wenchuan earthquake disaster relief and researches[J].BullSeismolSocAm， 100(5B): 2866--2872. doi:10.1785/0120090257.

Cross-correlation analysis of three-component seismic recordings

1)InstituteofGeophysics，ChinaEarthquakeAdministration，Beijing100081，China2)ChinaProperty&CasualtyReinsuranceCompanyLtd.,Beijing100033,China

Seismic wave is a three dimensional vector wavefield， the single component recordings are actually the projection of particle motion along certain directions. Based on the single component seismogram cross-correlation formula， this paper presents a novel simple solution which is suitable for the calculation of three-component seismogram cross-correlation， and the effectiveness of this new approach is verifiedviaa practical case from the aftershock sequence of 2008 WenchuanMS8.0 earthquake. Compared with single component seismogram cross-correlation， the new approach can obtain a global optimized result more reasonably and erase the discrepancy between different components in the work of template waveform matching. Also this new formula can take advantage of the congenerous between different components， suppress the ambient seismic noise effectively， and its rationality was demonstrated in theory. This new approach requires rather small computations as its simplicity in principles and procedures， which is suitable for the seismic data processing in the current era of “big data”.

three-component seismogram; multiple cross-correlation; template matching method

國家留學(xué)基金委和中國地震局“地震科技青年骨干人才培養(yǎng)項目”(20105007)聯(lián)合資助.

2015-02-17收到初稿， 2015-05-07決定采用修改稿.

e-mail: matengfei@cpcr.com.cn

10.11939/jass.2016.01.009

P315.63

A

馬騰飛. 2016. 三分量地震記錄的互相關(guān)分析. 地震學(xué)報, 38(1): 96--102. doi:10.11939/jass.2016.01.009.

Ma T F. 2016. Cross-correlation analysis of three-component seismic recordings.ActaSeismologicaSinica, 38(1): 96--102. doi:10.11939/jass.2016.01.009.