王亞娟，李懷良，2，庹先國，2，3，沈 統(tǒng)

（1.西南科技大學(xué)核廢物與環(huán)境安全國防重點學(xué)科實驗室，四川綿陽621010；2.成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點實驗室，四川成都610059；3.四川輕化工大學(xué)，四川自貢643002）

我國頁巖氣資源豐富，其開發(fā)促進(jìn)了能源結(jié)構(gòu)的調(diào)整和改變[1]。其中，水力壓裂為頁巖氣的勘探開發(fā)提供了重要的技術(shù)手段[2]，而水力壓裂微地震監(jiān)測則為頁巖氣優(yōu)化開發(fā)方案的制定提供了重要的技術(shù)支持。通過微地震監(jiān)測技術(shù)對微地震事件定位、能量大小估計及震源機(jī)理分析等可了解裂縫的形態(tài)分布及發(fā)育過程、裂縫間的相互作用、估算儲層改造體積、進(jìn)行儲層應(yīng)力分析和開展地震災(zāi)害評價[3-5]。

其中，精確的P 震相初至拾取是上述工作實施的重要基礎(chǔ)步驟。雖然此項工作可以手工進(jìn)行，但隨著數(shù)據(jù)量的急劇增加，視覺分析已然成為一項非常耗時和繁瑣的工作[6]。因此，有必要提供一種更高效、快速、客觀、自動的P震相初至拾取方法。

現(xiàn)有的P震相初至拾取方法大部分來源于對天然地震數(shù)據(jù)的處理，例如：能量分析[7-10]、極化分析[11-13]、人工神經(jīng)網(wǎng)絡(luò)[14-15]、自回歸技術(shù)[16-17]、高階統(tǒng)計量[18-22]、小波變換[6]等，這些方法各有特點，在天然地震領(lǐng)域取得了較好的應(yīng)用效果。其中，高階統(tǒng)計量是識別非高斯信號的有效工具，SARAGIOTIS等[19]應(yīng)用基于高階統(tǒng)計函數(shù)的偏度和峰度函數(shù)來識別P震相到達(dá)時間。相對于偏度函數(shù)，峰度法對非高斯信號更加敏感，且對高斯噪聲有一定的抑制作用。當(dāng)信號受到相對較低的噪聲污染時，該方法大多能較好地確定P震相的到達(dá)時。然而，在微地震監(jiān)測領(lǐng)域，由于檢波器自身的技術(shù)局限性以及外部監(jiān)測環(huán)境等因素的影響，實際采集的微地震數(shù)據(jù)不可避免地與高頻噪聲和工頻噪聲混合在一起，信噪比較低且初動不明顯。在這種情況下，峰度法大多效果不佳。

為提高峰度法對微地震信號的拾取精度，本文提出了一種基于特征函數(shù)構(gòu)建的峰度與小波多尺度分解的P震相初至拾取新方法?？紤]到各種噪聲的隨機(jī)性，該方法利用小波多尺度分解提取含強(qiáng)噪聲微地震信號的主成分，將其高頻信號剝離，提高微地震信號的質(zhì)量。在獲得有效信號后，構(gòu)建其特征函數(shù)，并計算該特征函數(shù)序列的峰度值，該峰度值的全局最大斜率即為所拾取的初至點。

1 峰度法

峰度是衡量信號高斯度的四階統(tǒng)計量。服從近似高斯分布的信號，其峰度值接近于0；而非高斯信號（如微地震事件）產(chǎn)生的峰度值較高。因此，峰度可以作為識別非高斯特征信號的有效統(tǒng)計工具[19]。其計算過程如下。

x（n）為微地震信號，其中，n＝1，2，…，為采樣點個數(shù)。采用長度為M的移動時間窗口掃描微地震序列，當(dāng)窗口從M滑動到n（M＜n）時，定義每個滑動窗口的峰度值計算公式為：

其中，

為計算微地震記錄的峰度值序列，需設(shè)置一個移動時間窗口沿時間軸掃描微地震記錄，窗口大小M＝200，移動步長為1，即1～200、2～201、3～202…。如圖1所示，模擬一組信噪比為10 dB 的微地震信號，其信號主頻為120 Hz，采樣間隔為0.1667×10－3s，采樣長度為2000 個采樣點，設(shè)置P 震相初至為1000。在圖1a中，A，B，C分別代表時間窗口移動的不同位置。時窗在沿時間軸移動的過程中，計算每個窗口位置的峰度值，由此得到描述微地震信號非高斯性的峰度曲線，為了與信號的橫軸對齊，峰度曲線的橫軸起點應(yīng)為200，終點為2000。如圖1b所示，該曲線由平穩(wěn)段、陡升段和緩降段組成。在P 波到達(dá)前，信號僅由噪聲構(gòu)成，峰度曲線趨于平穩(wěn)，峰度值為0；隨著時窗的移動，信號成分由純噪聲信號向微地震和噪聲的混合信號變化，此時，峰度值會隨信號成分的變化而發(fā)生變化；P 波到達(dá)一定時間后，信號成分由噪聲與微地震信號的混合信號向純噪聲信號變化，此時體現(xiàn)微地震信號非高斯性的峰度值會逐漸減小，最后趨于0。由上述分析可知，在P波到達(dá)附近峰度曲線會發(fā)生明顯的突變，由此可將峰度曲線的全局最大斜率定義為P震相的初至。圖1a中，黑色虛線表示模擬信號P 震相的初至，它與峰度曲線的最大斜率處重合。由此可以看出，峰度法對具有清晰P 震相起跳的微地震信號有很好的拾取效果。

然而，在頁巖氣開發(fā)微地震監(jiān)測過程中，實測的微地震記錄信噪比普遍較低，并常伴隨強(qiáng)尾波、尖峰信號，給P震相的初至拾取造成很大影響。圖2給出了低信噪比模擬微地震信號及其峰度值曲線。圖2a為信噪比為－6 dB 的模擬微地震信號，采樣間隔為0.1667×10－3s，P 震相初至點對應(yīng)第1000個采樣點，如圖中黑色虛線所示；圖2b為微地震信號的峰度值曲線，最大斜率對應(yīng)第1031個采樣點，如圖中黑色虛線所示，拾取結(jié)果滯后了31個采樣點，即延遲了5.1677×10－3s（采樣點個數(shù)×采樣間隔）。由圖2可知，對于低信噪比微地震信號，僅使用峰度法拾取P震相初至?xí)磔^大的拾取誤差，需要改進(jìn)其計算對象及過程，以提高拾取精度。

圖1 高信噪比模擬微地震信號波形（a）及其峰度值曲線（b）（信噪比為10 dB）

圖2 低信噪比模擬微地震信號波形（a）及其峰度值曲線（b）（信噪比為－6 dB）

2 小波多尺度分解及P震相初至拾取

2.1 小波多尺度分解

小波變換作為分析復(fù)雜非平穩(wěn)信號的有效工具，在提高地震資料分辨率和信噪比[20，23]、壓縮地震數(shù)據(jù)[24-25]、表征介質(zhì)奇異結(jié)構(gòu)[26]、地震資料反演[27-28]等方面均得到了廣泛的應(yīng)用。在微地震數(shù)據(jù)處理過程中，小波基的選取以及分解層次的設(shè)定，均會對信號后續(xù)的分析處理產(chǎn)生影響。其中，小波基的選擇必須具有對稱性才能有效防止地震相位失真[29]，而適當(dāng)分解層次的設(shè)定，不僅能減少計算量，還能有效分離低頻信號與高頻信號。經(jīng)過多次數(shù)值模擬實驗，在微地震數(shù)據(jù)處理過程中采用db10小波基最為合適，該小波基具有正交性、對稱性和緊支撐等特點；分解層次設(shè)定為3。

使用db10小波基對微地震信號x（n）進(jìn)行3層小波分解，如圖3所示，可得到尺度空間的近似序列a1，a2，a3和細(xì) 節(jié) 空 間 的 細(xì) 節(jié) 序 列d1，d2，d3。其 中近似序列的重構(gòu)信號對應(yīng)原信號每層分解的低頻部分，細(xì)節(jié)序列的重構(gòu)信號對應(yīng)原信號每層分解的高頻部分。為有效凸顯微地震信號的主成分，降低高頻信號對P 震相初至拾取的影響[29]，本文對近似序列a1，a2，a3進(jìn)行離散小波逆變換，從中選取合適的重構(gòu)近似信號用于P震相初至拾取。

圖3 小波多尺度分解示意

圖4給出了主頻為120 Hz的模擬微地震信號（信噪比為－6 dB）及小波多尺度分解后各層低頻重構(gòu)信號。圖4a為模擬微地震信號；圖4b、圖4c和圖4d分別為近似序列a3，a2和a1通過離散小波逆變換得到的重構(gòu)信號A3，A2，A1。從圖4可以看出，A1，A2，A3分別是濾除部分高頻噪聲后的有效信號，但A1、A2中依舊包含較多高頻成分，而A3在保持相對清晰的P 震相初至信息的同時，剔除了大部分高頻成分。因此，文中將重構(gòu)信號A3用于后續(xù)初至拾取。

為進(jìn)一步確認(rèn)上述小波基選擇、分解層次設(shè)定以及信號選取的有效性，圖5給出了合成微地震信號與重構(gòu)信號A3的時頻分析結(jié)果。從圖5a可以看出，相對于有效微地震信號，噪聲信號頻率較高。圖5b為所選重構(gòu)信號A3及其時頻分析結(jié)果。由圖5b可見，小波多尺度分解可以實現(xiàn)高頻噪聲和有效信號的分離，信噪比可以得到極大提高，而且重構(gòu)信號A3與原始信號的能量信息在時頻譜上能保持一致，證明了A3可以代替原始數(shù)據(jù)進(jìn)行初至拾取。

2.2 特征函數(shù)構(gòu)建

為進(jìn)一步提高本文方法的初至拾取精度，重新構(gòu)建了特征函數(shù)序列。常見的特征函數(shù)通常包括以下4種[30-31]：

特征函數(shù)（4）和特征函數(shù)（5）在時域上表現(xiàn)較為突出，體現(xiàn)信號的幅值變化；特征函數(shù)（6）和特征函數(shù)（7）體現(xiàn)信號的頻率變化。LI等[31]指出僅用（6）式作為信號的特征函數(shù)時，若背景噪聲的頻率高于地震信號的頻率，則該特征函數(shù)無效。因此，特征函數(shù)的好壞直接影響初至拾取的結(jié)果。為充分反映信號頻率與幅值的變化，并兼顧微地震信號的低信噪比和多震相特點，本文構(gòu)建了如下特征函數(shù)：

根據(jù)峰度法原理求特征函數(shù)序列的峰度值：

^CFi和δi的表達(dá)式為：

圖4 低信噪比模擬微地震信號波形及小波分解后的各層低頻重構(gòu)信號（黑色實線表示模擬信號的初至點）

圖5 合成微地震信號（a）與重構(gòu)信號A 3（b）的時頻分析結(jié)果

2.3 P震相初至拾取

本文方法的步驟可概括如下。

1）對于低信噪比微地震信號序列x（n），利用小波多尺度分解將其高頻噪聲與有效信號分離，選取重構(gòu)信號A3用于初至拾取。其中所選小波基為db10，分解層次為3。

2）通過（8）式構(gòu)建重構(gòu)信號A3的特征函數(shù)序列。

3）將峰度值計算應(yīng)用于所構(gòu)建的特征函數(shù)序列，選取峰度值的全局最大斜率作為P 震相初至。其中，時窗的設(shè)置對拾取精度和效率有至關(guān)重要的影響[20-21]，時窗過長不僅會增加計算量，而且會造成信息冗余，影響拾取精度。經(jīng)大量數(shù)值模擬實驗，本文方法的時窗長度設(shè)置為200。

圖6對比了峰度法、本文方法、小波分解與高階統(tǒng)計量聯(lián)合方法[20]對低信噪比（信噪比為－6 dB）模擬微地震數(shù)據(jù)P震相的拾取結(jié)果。圖6a為模擬微地震信號波形及上述3種方法的初至拾取結(jié)果，圖中黑色實線表示預(yù)設(shè)的P 震相初至，紅色實線表示本文方法的拾取結(jié)果，綠色實線表示小波分解與高階統(tǒng)計量聯(lián)合方法的拾取結(jié)果，藍(lán)色實線表示峰度法的拾取結(jié)果；圖6b為圖6a的局部放大圖。從圖6可以清晰看出，本文方法的拾取效果更好，其余兩種方法的拾取誤差均大于本文方法。

圖6 不同方法的模擬微地震數(shù)據(jù)P震相初至拾取結(jié)果（信噪比為－6 dB；P震相初至在第1 000個樣點處）

3 應(yīng)用效果分析

3.1 模擬微地震數(shù)據(jù)實驗

圖7 模擬微地震信號（黑色實線表示初至點在第1 000個樣點處）

為了驗證本文方法拾取P震相初至的可行性和穩(wěn)定性，分別采用本文方法、峰度法、小波分解與高階統(tǒng)計量聯(lián)合方法對模擬微地震信號進(jìn)行拾取精度與效率對比。圖7所示為模擬微地震信號，該信號是由衰減的正弦波產(chǎn)生，主頻為120 Hz，采樣間隔為0.1667×10－3s，數(shù)據(jù)長度為2000個采樣點，P 震相初至點預(yù)設(shè)在第1000個樣點處。

原始信號與加入噪聲的能量比值為信噪比[32]：

式中：T為信號長度；S，N分別表示原始信號與加入的噪聲信號。采用公式（12）構(gòu)造信噪比為－5，－6，－7，－8，－9，－10 dB 的含噪微地震信號，對應(yīng)每個信噪比值構(gòu)造1000組模擬微地震信號，分別采用上述3種方法對其進(jìn)行初至拾取，進(jìn)而計算相同信噪比模擬微地震信號拾取誤差的平均值，結(jié)果如表1 所示。對比3種方法的拾取結(jié)果可見，本文方法的拾取精度更高。實驗表明，本文方法的拾取誤差能夠控制在0.0302×10－3～1.3002×10－3s。

在相同計算平臺上，采用信噪比為－6 dB 的1000組模擬微地震信號對上述3種方法的計算效率進(jìn)行比較，統(tǒng)計結(jié)果如表2所示。由表2可見，相比于傳統(tǒng)峰度法，本文方法的計算效率有所降低，這是由于小波多尺度分解與特征函數(shù)的構(gòu)建增加了計算量，時間消耗有所提高。但效率的降低隨之而來的是拾取精度的提高。本文方法中小波基選取、分解層次的設(shè)定及拾取準(zhǔn)則均優(yōu)于小波變換與高階統(tǒng)計量聯(lián)合方法，因此本文方法的計算效率略高。由表1 和表2可知，本文方法相比于小波分解與高階統(tǒng)計量聯(lián)合方法在拾取精度和計算效率方面都有所提高。

表1_3種方法對不同信噪比模擬微地震記錄的P震相初至拾取誤差統(tǒng)計

3.2 應(yīng)用實例

為檢測本文方法的實用性，分別采用峰度法、本文方法對實測微地震記錄的P 震相初至進(jìn)行拾取，并與人工拾取結(jié)果進(jìn)行了對比。

圖8給出了采用不同方法對6組微地震信號的P震相初至拾取結(jié)果。其中，各圖上圖為實測微地震記錄的振幅-時間曲線，其采樣間隔為0.1667×10－3s；下圖為實測微地震信號的重構(gòu)信號A3的局部放大圖，圖中黑色實線為人工初至的拾取結(jié)果，紅色實線為本文方法對原信號初至的拾取結(jié)果，藍(lán)色實線表示傳統(tǒng)峰度法對原信號初至的拾取結(jié)果。從圖中可以看出，相比于傳統(tǒng)方法，本文方法拾取結(jié)果與人工拾取結(jié)果更為接近。其中，圖8a為P 震相初至相對清晰的微地震信號初至拾取結(jié)果，由圖8a可見，與傳統(tǒng)峰度法相比，本文方法的拾取精度提高了1.002×10－3s；圖8b、圖8c為具有強(qiáng)尾波的微地震信號初至拾取結(jié)果，由圖8b和圖8c可見，與傳統(tǒng)方法相比，本文方法可有效克服尾波信號的干擾，圖8b中本文方法的拾取精度比傳統(tǒng)峰度法提高了12.6692×10－3s，圖8c中本文方法的拾取精度比傳統(tǒng)峰度法提高了9.0018×10－3s；圖8d為伴隨有尖峰的微地震信號初至拾取結(jié)果，圖8d中本文方法的拾取精度比傳統(tǒng)峰度法提高了133.36×10－3s；圖8e、圖8f為受環(huán)境噪聲影響較強(qiáng)的微地震信號初至拾取結(jié)果，圖8e和圖8f中本文方法的拾取精度比傳統(tǒng)峰度法分別提高了7.5015×10－3s和5.6678×10－3s。綜上所述，本文方法的拾取結(jié)果與人工拾取結(jié)果更相近，實現(xiàn)了更準(zhǔn)確的P震相初至拾取。

圖8 不同方法對6組微地震信號的P震相到達(dá)時拾取結(jié)果

4 結(jié)論

針對頁巖氣微地震監(jiān)測資料信噪比低、初至難以拾取的問題，本文提出了一種基于特征函數(shù)構(gòu)建的峰度與小波多尺度分解的P震相自動拾取方法。相比于傳統(tǒng)峰度法，該方法既保持了小波分解的多尺度分析特性以及對隨機(jī)噪聲的抑制特性，又克服了傳統(tǒng)峰度法對低信噪比、強(qiáng)尾波和尖峰的微地震信號拾取不準(zhǔn)確的缺點。將該算法用于不同信噪比模擬微地震信號的測試，并與傳統(tǒng)峰度法、小波分解與高階統(tǒng)計量聯(lián)合方法進(jìn)行比較，結(jié)果表明，本文方法能夠?qū)Φ托旁氡任⒌卣饠?shù)據(jù)的P 震相進(jìn)行準(zhǔn)確拾取，且誤差為0.0302×10－3～1.3002×10－3s，具有一定的魯棒性和準(zhǔn)確性。將本文方法應(yīng)用于實際微地震記錄的P震相初至拾取，進(jìn)一步驗證了本文方法的有效性和實用性。本文只是針對P 波初至拾取進(jìn)行研究，然而，在微地震數(shù)據(jù)分析過程中僅對P 波進(jìn)行拾取是遠(yuǎn)遠(yuǎn)不夠的，適當(dāng)?shù)腟波模型和可靠的S波到達(dá)時為震源定位、地震層析成像均提供了重要的約束。因此我們下一步將研究如何實現(xiàn)P波、S波的同時拾取。