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

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

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

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

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

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

1 峰度法

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

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

其中，

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

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

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

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

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

2.1 小波多尺度分解

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

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

圖3 小波多尺度分解示意

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

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

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

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

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

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

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

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

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

2.3 P震相初至拾取

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

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

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

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

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

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

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

3.1 模擬微地震數(shù)據(jù)實(shí)驗(yàn)

圖7 模擬微地震信號(hào)（黑色實(shí)線表示初至點(diǎn)在第1 000個(gè)樣點(diǎn)處）

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

原始信號(hào)與加入噪聲的能量比值為信噪比[32]：

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

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

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

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

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

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

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

4 結(jié)論

針對(duì)頁(yè)巖氣微地震監(jiān)測(cè)資料信噪比低、初至難以拾取的問(wèn)題，本文提出了一種基于特征函數(shù)構(gòu)建的峰度與小波多尺度分解的P震相自動(dòng)拾取方法。相比于傳統(tǒng)峰度法，該方法既保持了小波分解的多尺度分析特性以及對(duì)隨機(jī)噪聲的抑制特性，又克服了傳統(tǒng)峰度法對(duì)低信噪比、強(qiáng)尾波和尖峰的微地震信號(hào)拾取不準(zhǔn)確的缺點(diǎn)。將該算法用于不同信噪比模擬微地震信號(hào)的測(cè)試，并與傳統(tǒng)峰度法、小波分解與高階統(tǒng)計(jì)量聯(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)用于實(shí)際微地震記錄的P震相初至拾取，進(jìn)一步驗(yàn)證了本文方法的有效性和實(shí)用性。本文只是針對(duì)P 波初至拾取進(jìn)行研究，然而，在微地震數(shù)據(jù)分析過(guò)程中僅對(duì)P 波進(jìn)行拾取是遠(yuǎn)遠(yuǎn)不夠的，適當(dāng)?shù)腟波模型和可靠的S波到達(dá)時(shí)為震源定位、地震層析成像均提供了重要的約束。因此我們下一步將研究如何實(shí)現(xiàn)P波、S波的同時(shí)拾取。