馬騰飛

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

?

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

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

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

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

引言

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

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

1 基本方法與理論

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

(1)

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

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

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

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

(2)

其中，

(3)

(4)

(5)

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

2 應(yīng)用實(shí)例

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

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

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

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

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

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

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

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

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

3 討論與結(jié)論

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

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

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

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

蔣長(zhǎng)勝， 吳忠良， 李宇彤. 2008. 首都圈地區(qū)“重復(fù)地震”及其在區(qū)域地震臺(tái)網(wǎng)定位精度評(píng)價(jià)中的應(yīng)用[J]. 地球物理學(xué)報(bào)， 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é)問(wèn)題[D]. 北京: 中國(guó)地震局地球物理研究所: 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

國(guó)家留學(xué)基金委和中國(guó)地震局“地震科技青年骨干人才培養(yǎng)項(xià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é)報(bào), 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.