P波初至到時(shí)自動(dòng)拾取技術(shù)研究

王俊 楊馳 詹小艷 朱峰 張金川 錢(qián)婷

摘要：基于粘滯性單自由度振動(dòng)器響應(yīng)下的能量轉(zhuǎn)換理論，提出利用阻尼能量作為目標(biāo)函數(shù)的P波震相到時(shí)拾取方法-SDOF Picker算法。使用該方法對(duì)江蘇及鄰區(qū)2010-2016年實(shí)際記錄的9607組P波初至進(jìn)行到時(shí)自動(dòng)拾取測(cè)試，以地震編目中人工拾取到時(shí)為基準(zhǔn)，與利用AIC算法自動(dòng)抬取的結(jié)果進(jìn)行了系統(tǒng)性對(duì)比分析，結(jié)果顯示：SDOF Picker算法和AIC算法自動(dòng)拾取P波初至的準(zhǔn)確率分別為97.1%、91.8%，中值偏差分別為（0.02±0.61）s、（0.05±0.77）s，方差分別為0.37s2、0.60s2，這表明SDOF Picker算法的在準(zhǔn)確率和拾取精度方面均優(yōu)于AIC算法。

關(guān)鍵詞：粘滯性單自由度振動(dòng)器；阻尼能量；P波初至到時(shí)；AIC算法；SDOF picker算法

中圖分類(lèi)號(hào)：P315.3 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1000-0666（2018）01-0038-08

0 引言

地震定位、震源機(jī)制求解及地震動(dòng)預(yù)測(cè)等許多地震學(xué)應(yīng)用中，P波到時(shí)的精確拾取均十分關(guān)鍵。目前，地震臺(tái)網(wǎng)通常采用半自動(dòng)處理方法，即通過(guò)數(shù)值算法拾取一個(gè)大致位置，再由分析人員進(jìn)行復(fù)核、修正。這是因?yàn)橐话闱闆r下人工拾取的精度更高，尤其是在信噪比較低時(shí)，但人工拾取震相到時(shí)存在主觀性且費(fèi)時(shí)。未來(lái)隨著地震臺(tái)站規(guī)模不斷增大，特別是在地震早期預(yù)警與地震動(dòng)預(yù)測(cè)中，從海量連續(xù)記錄波形數(shù)據(jù)中自動(dòng)、精準(zhǔn)地拾取震相到時(shí)的需求將愈加迫切。

1980年以來(lái)，出現(xiàn)了許多檢測(cè)與拾取不同地震波到時(shí)的方法，大多是通過(guò)尋找地震信號(hào)與背景噪聲（或長(zhǎng)時(shí)平均振動(dòng)）之間能量、極化特征以及頻譜成分的變化特征來(lái)進(jìn)行的（楊配新等，2004；王彩霞等，2013）。通常情況下，地震信號(hào)在某一特定的頻帶范圍內(nèi)會(huì)被強(qiáng)化或沿一定的方向極化（Lomax et al，2012），大多數(shù)檢測(cè)方法通常還會(huì)對(duì)記錄信號(hào)進(jìn)行去除線(xiàn)性趨勢(shì)、濾波等預(yù)處理，以期降低背景噪聲或提高信號(hào)強(qiáng)度。在眾多P波到時(shí)檢測(cè)方法中，STA/LTA方法的應(yīng)用最為普遍（Allen，1978；Earle，Shearer，1994）。該方法是以對(duì)比信號(hào)短時(shí)平均/長(zhǎng)時(shí)平均特征為基礎(chǔ)，通常將隨時(shí)間的絕對(duì)值或均方根作為特征函數(shù)，當(dāng)STA/LTA超過(guò)預(yù)設(shè)或設(shè)定的動(dòng)態(tài)觸發(fā)閾值時(shí)，視為有效檢測(cè)。該方法對(duì)檢測(cè)振幅的變化非常有效，但它的精度依賴(lài)于設(shè)定的檢測(cè)間隔與閾值；另外與地震信號(hào)無(wú)關(guān)的脈沖型噪聲，尤其疊加在微震記錄中時(shí)會(huì)引起誤觸發(fā)檢測(cè)。

自回歸（AR）方法也常角來(lái)檢測(cè)震相。這類(lèi)方法是假設(shè)地震波能被分解為部分穩(wěn)態(tài)的多個(gè)過(guò)程，每個(gè)模型作為一個(gè)AR處理，初至前、后被視為2種具有不同性質(zhì)的狀態(tài)（Sleeman，Van Eek，1999； Rastin et al，2013）。為檢測(cè)到理想的觸發(fā)時(shí)間，需要分離初至前（信號(hào)中僅包含背景噪聲）與初至后（信號(hào)中包含地震波和背景噪聲）特征，使用AR技術(shù)分析不同窗內(nèi)的時(shí)間信號(hào)。例如，一般的背景噪聲通過(guò)相對(duì)低階的AR就可以很好地表征，而地震信號(hào)則通常需要高階的AR才能較好地表征（Leonard，Kennett，1999）。采用AR分析處理時(shí)序信號(hào)時(shí)，赤池信息準(zhǔn)則（AIC）能幫助確定AR處理的階數(shù)，這就意味著模型匹配同樣存在不確定性（Akaike，1974）。AR-AIC是另一種常用的有效方法（Leonard，2000；趙大鵬等，2012），但當(dāng)時(shí)間窗內(nèi)存在多重震相時(shí)，AIC撿拾器則會(huì)拾取振幅最大的震相，另外，它不能辨別時(shí)間窗內(nèi)的信號(hào)是真正的震相到時(shí)，還是由傳感器故障或數(shù)字噪聲引起的干擾信號(hào)（Kalkan，Stephens，2017）。除此之外，神經(jīng)網(wǎng)絡(luò)方法（Dai，MacBeth，1995）、小波變換（高靜懷等，1997；Anant，Dowla，1997；劉希強(qiáng)等，2000）、滑動(dòng)互相關(guān)能量比（Forghani-Arani et al，2013）等方法也被用來(lái)拾取震相到時(shí)，但這些方法都有一個(gè)共性，即需要一個(gè)或多個(gè)經(jīng)驗(yàn)性預(yù)設(shè)，如檢測(cè)間隔、閾值設(shè)定、事件特征、信號(hào)的頻譜成分等。

本文基于粘滯性單自由度（簡(jiǎn)稱(chēng)SDOF）振動(dòng)器響應(yīng)下的能量轉(zhuǎn)換理論，提出利用阻尼能量作為目標(biāo)函數(shù)的P波震相到時(shí)拾取方法——SDOF Pick-er算法。再利用2010-2016年江蘇省地震臺(tái)網(wǎng)記錄到的9607條P波波形進(jìn)行詳細(xì)測(cè)試，并分析了拾取精度。

1 基本原理

在粘滯性單自由度振動(dòng)系統(tǒng)響應(yīng)下，由記錄信號(hào)傳遞給振動(dòng)器的輸入能量會(huì)被轉(zhuǎn)換為彈性應(yīng)變能，隨后由阻尼元件產(chǎn)生阻尼能量進(jìn)行耗散，阻尼力與速度成正比（歐朱光，2003）。圖1給出典型粘滯阻尼單自由度振蕩器的示意圖。

這2個(gè)振動(dòng)器的運(yùn)動(dòng)方程，可以用一個(gè)普通二階微分方程來(lái)表示：

mut+cu+ku=0（1）式中：m為質(zhì)量；c是阻尼系數(shù)；k是剛度系數(shù)；ut是質(zhì)量塊的絕對(duì)（總）位移，且ut=u+ug，ug是地面運(yùn)動(dòng)位移，u是質(zhì)量塊相對(duì)于底座（如地面）的相對(duì)位移。式（1）可以寫(xiě)為：

mu+cu+ku=-mug（2）式（2）的右邊是與地動(dòng)加速度（ug）有關(guān)的激發(fā)函數(shù)。對(duì)式（1）、（2）中u進(jìn)行積分將得到2種不同定義的輸入能量。對(duì)式（1）中的u進(jìn)行積分，給出了線(xiàn)彈性粘滯阻尼SDOF振蕩器運(yùn)動(dòng)（圖1a）的絕對(duì)（總）能量公式，其遭受的運(yùn)動(dòng)加速度（ug）可以表示為：式中：t代表時(shí)間。式（3）還能寫(xiě)成更一般的表達(dá)式，用能量來(lái)表示：

Ek+Eζ+ES=EI（4）式中：Ek是絕對(duì)動(dòng)能；Eζ是阻尼能量；ES是彈性應(yīng)變能。這些項(xiàng)相加就等于傳遞給SDOF振蕩器的能量，即絕對(duì)輸入能量EI。

對(duì)式（2）進(jìn)行積分，可得出基準(zhǔn)固定情況下（圖1b）線(xiàn)彈性粘滯阻尼SDOF振蕩器的相對(duì)能量公式：式（5）可寫(xiě)成：

Ek'+Eζ+ES=EI'（6）式中：Ek'是相對(duì)動(dòng)能；EI'是相對(duì)輸入能量。式（4）中的EI表示外力（mut）在質(zhì)量塊上的作用已經(jīng)完成，相當(dāng)于總基底-剪切力對(duì)地面運(yùn)動(dòng)位移上的功已經(jīng)完成。另一方面，式（6）中的EI'代表等效側(cè)力在基準(zhǔn)固定情況下振蕩器的功已完成，從而去除了剛體轉(zhuǎn)換效應(yīng)。2個(gè)能量公式之間EI'和EI的差異導(dǎo)致了2個(gè)動(dòng)能Ek'和Ek之間的差異，然而，阻尼能量和應(yīng)變能量項(xiàng)在這兩種定義中均被保留。在粘滯性阻尼SDOF振蕩器中，單位重量的阻尼能量耗散的功率，通過(guò)對(duì)Eζ進(jìn)行時(shí)間的微分，可表示為：式中：ωD是角頻率；ζ是阻尼比，用2π/TD來(lái)表示；TD是阻尼振動(dòng)的固有周期，與無(wú)阻尼振動(dòng)的固有周期（Tn）相關(guān)，可表示為：

2 P波初至到時(shí)拾取算法

在地震波激勵(lì)情況下，能量的變化可采用底座固定狀態(tài)下線(xiàn)彈性粘滯阻尼SDOF振動(dòng)器（圖1b）的相對(duì)能量公式來(lái)表述。在時(shí)域里，由式（2）可得到質(zhì)量塊的相對(duì)速度u。在u已知的情況下，利用式（5）便可計(jì)算出不同的能量值。由式（5）和式（7）可知，阻尼能量是累積的函數(shù)，且正比于質(zhì)量塊的相對(duì)速度，這使得阻尼能量是一個(gè)隨時(shí)間變化的光滑包絡(luò)線(xiàn)函數(shù)，在記錄開(kāi)始至P波到達(dá)之前接近于零（或平坦），P波到達(dá)后迅速增大（Erol，2016），本文提出的方法便是利用阻尼能量的平滑特性來(lái)進(jìn)行P波初至到時(shí)的檢測(cè)，簡(jiǎn)稱(chēng)為SDOF Picker算法。由于背景噪聲與地震波信號(hào)之間具有不同的頻率分布特征，因此阻尼能量函數(shù)在信號(hào)開(kāi)始時(shí)即會(huì)發(fā)生顯著變化，于是可將其作為跟蹤并檢測(cè)P波初至到時(shí)的有效度量。固有周期小的振動(dòng)器，由于不會(huì)產(chǎn)生延遲和震相以外的響應(yīng)，因此對(duì)于這種變化具有更高的敏感性。

在P波初至到時(shí)的拾取中為防止共振，將SD-OF振動(dòng)器的固有周期設(shè)置為0.01s，因此具有相對(duì)較高的諧振頻率，該頻率高于地震波信號(hào)中大部分頻率；阻尼比取值為0.6，這類(lèi)似于短周期地震計(jì)。在0.6的阻尼水平或更高的狀態(tài)下，頻率響應(yīng)接近于巴特沃斯（Butterworth）濾波器所能保留的最大的平坦幅度和相位角響應(yīng)，這意味著輸入信號(hào)的頻譜能完整地被傳遞到振蕩器的響應(yīng)中。在這樣的固有周期和阻尼比狀態(tài)下，振蕩器能迅速回到平衡位置，且沒(méi)有自由振蕩也不影響輸入能量特性。圖2是2017年2月24日江蘇金湖ML1.5地震寶應(yīng)臺(tái)垂直分向的記錄波形及各能量計(jì)算結(jié)果，震中距約為11.7km，從圖中看到，相對(duì)輸入能量是相對(duì)動(dòng)能（圖2c）、彈性應(yīng)變能量（圖2d）、阻尼能量（圖2e）的總和，大部分作為應(yīng)變能被保留和釋放，隨后由阻尼器進(jìn)行能量耗散，振動(dòng)器的動(dòng)能可忽略不計(jì)，因?yàn)橘|(zhì)量塊相對(duì)于基準(zhǔn)位置幾乎沒(méi)有運(yùn)動(dòng)。最后，通過(guò)檢測(cè)阻尼能量功率（圖2f）來(lái)進(jìn)行震相到時(shí)拾取。

為驗(yàn)證SDOF Picker算法的拾取精度，本文將其與AIC算法、人工拾取的結(jié)果進(jìn)行對(duì)比分析。圖3是2017年江蘇金湖ML1.5地震寶應(yīng)臺(tái)三分向記錄的檢測(cè)結(jié)果，記錄波形采樣率為100Hz，記錄中BHE分向的噪聲水平顯著高于其它2個(gè)分向，這3種方法拾取的P波到時(shí)位置用不同顏色的豎直線(xiàn)進(jìn)行標(biāo)注（圖3a）。在這次地震中SDOF Picker算法在BHZ，BHE、BHN三分向上的撿拾結(jié)果與人工拾取結(jié)果之間偏差分別僅為0.01s，0.02s，0.01s，與AIC算法的拾取結(jié)果之間偏差則分別為0.02s、2.83s、0.02s；從圖3中可看出，盡管AIC算法在BHZ、BHN分向上的檢測(cè)結(jié)果與人工分析結(jié)果也十分接近，但在具有高噪聲干擾的BHE分向上，AIC算法則未正確拾取出P波到時(shí)，拾取出的位置為S波震相到時(shí)。為更進(jìn)一步檢驗(yàn)SDOF Picker算法在低信噪比下的拾取性能，定義信噪比為人工初至拾取震相到時(shí)前、后各1s的振幅平方之比（Hildyard et al，2008），圖5給出了P波段信噪比分別為1.16、3.24、5.81時(shí)江蘇泗洪臺(tái)（SH）、徐州臺(tái)（XZ）及連云港臺(tái)（LYG）垂直向記錄的檢測(cè)實(shí)例，其中SH臺(tái)、XZ臺(tái)的儀器頻帶范圍為60s～50Hz，LYG臺(tái)的儀器頻帶范圍為120s～50Hz，記錄采樣率均為100Hz，到時(shí)偏差對(duì)比結(jié)果見(jiàn)表1。

從表1中的3組實(shí)測(cè)結(jié)果看，隨著信噪比的提高，2種算法的結(jié)果與人工分析結(jié)果之間的偏差均逐漸減小，但SDOF Picker算法的檢測(cè)結(jié)果與人工分析結(jié)果之間偏差更小，意味著即使是在信噪比極低的情況下，阻尼能量的功率也能精準(zhǔn)地反映出背景噪聲與地震波信號(hào)之間的差異特征，而從基于原始記錄波形的AIC曲線(xiàn)結(jié)果可以看出（圖5a-2、b-2），在低信噪比下的寬頻帶記錄中容易受到記錄中長(zhǎng)周期固有信號(hào)的“誤導(dǎo)”，致使震相到時(shí)拾取偏差較大。

3 系統(tǒng)性測(cè)試

為更系統(tǒng)地測(cè)試SDOF Picker算法的可靠性，本文對(duì)2010-2016年江蘇及鄰區(qū)1444次M0.1～4.9地震（圖6），選取其中震中距小于200km地震的9607個(gè)P波初至到時(shí)進(jìn)行自動(dòng)拾取，以江蘇省測(cè)震臺(tái)網(wǎng)人工正式編目結(jié)果中的P波到時(shí)為基準(zhǔn)，將自動(dòng)拾取到時(shí)與人工到時(shí)之間偏差大于2s的結(jié)果視為誤拾取，2種自動(dòng)算法拾取結(jié)果的準(zhǔn)確率與精度對(duì)比見(jiàn)表2。

從表2的9607組自動(dòng)拾取結(jié)果中可以看到，相比于AIC算法，SDOF Picker算法拾取的準(zhǔn)確率提高了5.3%，中值偏差減小了0.03s，方差減小了0.23s2，中值偏差和方差的減小意昧著與人工分析結(jié)果之間系統(tǒng)性偏差更小，自動(dòng)拾取精度優(yōu)于AIC算法；從絕對(duì)偏差統(tǒng)計(jì)分布結(jié)果中可以清晰地看到，SDOF Picker算法拾取結(jié)果的正態(tài)分布特征也更加顯著（圖7）。根據(jù)到時(shí)偏差與信噪比之間的統(tǒng)計(jì)結(jié)果（圖8）：在AIC算法自動(dòng)檢測(cè)結(jié)果中，與人工分析結(jié)果相比到時(shí)偏差小于0.5s的比例約為總樣本數(shù)的56.5%，而在SDOF Picker算法的自動(dòng)拾取結(jié)果中，這一比例則高達(dá)76.7，相對(duì)于AIC算法拾取效率提高了20.2%。當(dāng)信噪比范圍設(shè)定在較低的1～5時(shí)，震相樣本數(shù)為4 015組，此時(shí)AIC算法自動(dòng)拾取結(jié)果中偏差小于0.5s比例約占樣本數(shù)的38.9%，而在SDOF Picker算法的結(jié)果中這一比例則高達(dá)65.8%，相對(duì)于AIC算法拾取效率提高的幅度更大，約為26.9%。這表明SDOF Picker算法拾取性能更優(yōu)，尤其是在低信噪比下能正確拾取出的弱信號(hào)更多，拾取效率更高。

4 結(jié)論

本文基于粘滯性單自由度振動(dòng)器響應(yīng)下的能量轉(zhuǎn)換理論，提出利用阻尼能量作為目標(biāo)函數(shù)的P波初至拾取算法SDOF Picker算法。該算法將記錄信號(hào)轉(zhuǎn)換到粘滯性單自由度振動(dòng)器響應(yīng)下，由于阻尼能量函數(shù)是與時(shí)間相關(guān)的光滑一包絡(luò)線(xiàn)函數(shù)，在P震相到達(dá)時(shí)會(huì)迅速增大可以清晰地拾取出來(lái)，因此相比于傳統(tǒng)的檢測(cè)方法，如STA/LTA等，SDOFPicker算法不需要指定任何檢測(cè)間隔或設(shè)定觸發(fā)閾值，也不需要對(duì)記錄臺(tái)站的背景噪聲或被檢測(cè)事件等進(jìn)行先驗(yàn)研究，因此便于廣泛地進(jìn)行實(shí)際應(yīng)用轉(zhuǎn)化。利用2010-2016年期間江蘇地震臺(tái)網(wǎng)實(shí)際記錄的9607組P波初至測(cè)試了其性能，測(cè)試結(jié)果顯示：SDOF Picker算法拾取的準(zhǔn)確率和精度整體均優(yōu)于AIG算法，特別是在低信噪比記錄情況下，其拾取效率更高。

自動(dòng)識(shí)別和檢測(cè)震相對(duì)于地震早期預(yù)警和快速處理主震之后的大量余震非常關(guān)鍵，因此本文所描述的檢測(cè)算法對(duì)于高精度自動(dòng)拾取P波初至到時(shí)具有實(shí)際意義。目前，SDOF Picker算法已經(jīng)在江蘇省地震臺(tái)網(wǎng)進(jìn)行實(shí)際應(yīng)用，主要用于系統(tǒng)性自動(dòng)識(shí)別前震時(shí)間窗和基于自動(dòng)拾取的P波初至到時(shí)進(jìn)行地震定位。未來(lái)還將在此理論基礎(chǔ)上進(jìn)一步研究S波震相到時(shí)的自動(dòng)拾取技術(shù)。

參考文獻(xiàn)：

高靜懷，汪文秉，朱光明，等.1997.小波變換與信號(hào)瞬時(shí)特征分析[J].地球物理學(xué)報(bào)，40（6）；821-832.

劉希強(qiáng)，周蕙蘭，沈萍，等.2000門(mén)用于三分向記錄震相識(shí)別的小波變換方法[J].地震學(xué)報(bào)，22（2）：125-131.

歐朱光.2003.工程振動(dòng)[M].武漢：武漢大學(xué)出版社，17-55.

王彩霞，白超英，王馨.2013.地震震相初至自動(dòng)檢測(cè)技術(shù)綜述[J].地球物理學(xué)進(jìn)展，28（5）：2363-2375.

楊配新，鄧存華，劉希強(qiáng)，等.2004.數(shù)字化地震記錄震相自動(dòng)識(shí)別的方法研究[J].地震研究，27（4）308-313.

趙大鵬，劉希強(qiáng)，李紅，等.2012.峰度和AIC方法在區(qū)域地震事件和直達(dá)P波初動(dòng)自動(dòng)識(shí)別方面的應(yīng)用[J].地震研究，35（2）；220-225.

Akaike H.1974.A new look at the statistical model identification[J]Trans Automat Contr，19（6）：716-723.

Allen R V.1978.Automatic earthquake recognition and timing from singletraces[J].Bull Seismol Soc Am，68（5）：1521-1532.

Anant K S，Dowla F U.1997.Wavelet transform methods for phase identi-fication in three-component seismograms[J].Bull Seismol SocAm，87（6）：1598-1612.

Dai H，MacBeth C.1995.Automatic picking of seismic arrivals in localearthquake data using an artificial neural network[J].Geophys JInt，120（3）：758-774.

Earle P，Shearer P.1994.Characterization of global seismograms using anautomatic picking algorithm[J].Bull Seismol Soc Am，84（2）：366-376.

Erol K.2016.An Automatic P-Phase Arrival-Time Picker[J].BullSeismol Soc Am，106（3）：1-16.

Forghani-Arani F，Behura J，Haines S S，et al.2013.Anautomated cross-correlation based event detection technique and its application to asurface passive data set[J].Geophys Prospect，71（4）：778-787.

Hildyard M W，Nippress S E，Rietbrock A.2008.Event detection andphase picking using a time-domain estimate of predominate periodTpd[J].Bull Seismol Soc Am，98（6）：3025-3032.

Kalkan E，Stephens C.2017.Systematic Comparisons Between PRISMVersion 1.0.0，BAP，and CSMIP Ground-Motion Processing[R].U S Geol Surv Open-File Rept.

Leonard M，Kennett B L N.1999.Multi-component autoregressive tech-niques for the analysis of seismograms[J].Phys Earth Planet In，113，247-263.

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

Lomax A，Satriano C，Vassallo M.2012.Automatic picker developmentsand optimization；Filter Picker-A robust，broadband picker for real-time seismic monitoring and earthquake early warning[J].SeismolRes Lett，83（3）：531-540.

Rastin S J，Unsworth C P，BenitesR，et al.2013.Using real and syntheticwaveforms of the Matata swarm to assess the performance of NewZealand GeoNet phase pickers[J].Bull Seismol Soc Am，103（4）：2173-2187.

Sleeman R，Van Eck T.1999.Robust automatic P-phase picking；An on-line implementation in the analysis of broadband seismogram re-cordings[J].Phys Earth Planet In，113（1-4）：265-275.