王議迎，丁仁偉,2，3，李建平,4，趙俐紅,2，3，趙 碩，張碩偉

(1. 山東科技大學(xué) 地球科學(xué)與工程學(xué)院，山東 青島 266590; 2. 青島海洋科學(xué)與技術(shù)國家實驗室，山東 青島 266237; 3. 海洋礦產(chǎn)資源評價與探測技術(shù)國家功能實驗室，山東 青島 266237;4. 山東省物化探勘查院，山東 濟(jì)南 221116)

微地震監(jiān)測已經(jīng)成為當(dāng)今油藏地球物理描述和二氧化碳地質(zhì)儲存檢測等方面的關(guān)鍵技術(shù)手段[1]。在非常規(guī)油氣開采過程中，為了提高油氣產(chǎn)量,在儲層附近進(jìn)行高壓注液，從而產(chǎn)生微裂縫、誘發(fā)微地震。實時繪制微震震源，能更好地控制注入過程從而提高儲層開發(fā)的效率[2]。因此，微震的自動檢測對于監(jiān)測誘導(dǎo)水力壓裂過程至關(guān)重要。

隨著微震記錄數(shù)據(jù)規(guī)模的增大，手動拾取微震到時變得更具挑戰(zhàn)性和不可靠，且可能會產(chǎn)生若干假警報或出現(xiàn)到時誤拾現(xiàn)象。因此，各種自動或半自動的方法被提出并廣泛用于微地震震相拾取?，F(xiàn)有的微震震相自動拾取方法主要有長短時窗平均能量比法(short-term average/long-term average, STA/LTA)[3]，該算法簡單、計算效率高、適于實時處理。在過去的幾十年里，許多研究人員對該方法進(jìn)行了改進(jìn)，通過構(gòu)建不同的特征函數(shù)或構(gòu)建能反映信號振幅強度與頻率變化的權(quán)值函數(shù)等方式獲取更準(zhǔn)確的微震到時拾取結(jié)果[4-10]。自回歸(auto regressive, AR)是另一種常用的到時自動拾取方法[11]，該方法基于統(tǒng)計特性的差異將微震記錄分為真實微震信號到達(dá)前和到達(dá)后兩個部分，用以表征微震記錄。在此基礎(chǔ)上，Akaike[12]提出了Akaike信息準(zhǔn)則(Akaike information criterion, AIC)方法，將微震數(shù)據(jù)建模為AR過程，將AIC最小值設(shè)定為微震事件的起始時間。除此之外，高階統(tǒng)計法[13]、神經(jīng)網(wǎng)絡(luò)方法[14]、相關(guān)法[15]以及偏振分析方法[16]等也被廣泛應(yīng)用于微震震相到時拾取。然而，這些方法需要提前得知檢測間隔、人為設(shè)置閾值，并需要根據(jù)信號強度或時間域特征、頻率域特征、綜合特征等不斷調(diào)整拾取參數(shù)，降低了方法的普適性。

對于低信噪比的微震數(shù)據(jù)，微震記錄信號中包含大量隨機(jī)噪聲，這對微震震相的到時拾取是一項極大的挑戰(zhàn)。為了去除隨機(jī)噪聲對拾取震相到時準(zhǔn)確性的影響，常用噪聲壓制的方法來提高原始數(shù)據(jù)的信噪比。Chen等[17]提出基于f-x經(jīng)驗?zāi)B(tài)分解的預(yù)測濾波方法，克服了傳統(tǒng)f-x預(yù)測濾波處理多傾角微震剖面的局限性。Han等[18]提出一種結(jié)合自適應(yīng)閾值的集成經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition, EEMD)微地震降噪方法，通過解決模態(tài)混合問題改進(jìn)了經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)，能夠更好地將包含噪聲的微震數(shù)據(jù)分解為噪聲和有效信號兩部分。Huang等[19]提出數(shù)學(xué)形態(tài)學(xué)濾波(mathematical morphology filter, MMF)方法，在抑制低頻噪聲的同時保留信號的低頻成分。Li等[20]提出一種數(shù)學(xué)形態(tài)學(xué)理論，通過區(qū)分微震信號與背景噪聲來抑制低頻噪聲。這些濾波方法雖取得了較好的降噪效果，但在強隨機(jī)噪聲干擾下未能保留微震數(shù)據(jù)的有效信號形態(tài)。

高斯拉普拉斯(Laplacian of Gaussian, LoG)算子源于Marr和Hildreth[21]提出的視覺理論，自提出后就受到了廣泛關(guān)注。在圖像處理領(lǐng)域，LoG算子濾波器被廣泛應(yīng)用于邊緣檢測等眾多圖像處理領(lǐng)域。Saad等[22]提出改進(jìn)的高斯拉普拉斯(modified Laplacian of Gaussian, MLoG)算子濾波方法，利用高斯平滑濾波器和拉普拉斯銳化濾波器，通過平滑背景噪聲來檢測微震信號的突變。為準(zhǔn)確拾取微震P波到時，本研究利用MLoG算子濾波方法提高微震數(shù)據(jù)的信噪比，然后通過新提出的特征函數(shù)，在滑動時窗內(nèi)分別計算STA和LTA的值，繪制STA/LTA曲線，曲線的峰值記為微震P波到時。

1 改進(jìn)的STA/LTA方法

STA/LTA最早由Stevenson[3]用于拾取微震震相到時，通過STA與LTA之比反映微震信號振幅和能量的變化。STA主要反映微震信號能量的平均值，LTA主要反映背景噪聲能量的平均值。當(dāng)微震信號到達(dá)時，STA比LTA變化快，相應(yīng)的STA/LTA值會出現(xiàn)明顯的增加，當(dāng)某一時刻其比值大于預(yù)先設(shè)定的觸發(fā)閾值時，即可判定該時刻為P波的到時，從而達(dá)到自動檢測和拾取P波到時的目的。STA/LTA方法計算公式為[6]：

(1)

式中：STAi和LTAi(i=1，2，…，N)分別為微震信號在i時刻的短時窗平均能量值和長時窗平均能量值；NSTA和NLTA分別為對應(yīng)的短時窗長度和長時窗長度；CFj(j=1，2，…，M)為信號在j時刻的特征函數(shù)值，表征微震數(shù)據(jù)的振幅、頻率與能量的變化情況。

STA/LTA方法基于特征函數(shù)的突變進(jìn)行識別。不同特征函數(shù)對信號變化的靈敏度不同，具有不同的識別能力。選取適當(dāng)?shù)奶卣骱瘮?shù)對準(zhǔn)確拾取微震到時極為重要，常用的特征函數(shù)有[7-8]：

CF1=|X(i)|，

(2)

(3)

(4)

式中X(i)代表原始微震事件在第i個采樣點的輸入信號。特征函數(shù)CF1和CF2僅對振幅變化響應(yīng)；特征函數(shù)CF3對振幅變化和頻率變化都有響應(yīng)。劉晗等[9]指出特征函數(shù)CF2和CF3的檢測效果較好，但極易受到噪聲的干擾，對微震P波到時的自動拾取帶來一定誤差，從而影響拾取結(jié)果的準(zhǔn)確性。

(3)紙漿洗滌過程t時刻的操作參數(shù)為q(t)=[q1(t), q2(t), q3(t)]。其中，q1(t)為上漿濃度，q2(t)為上漿流量，q3(t)為清水加入量。

為了增強算法對振幅和頻率變化的敏感度，Allen[4]、劉曉明等[8]對傳統(tǒng)STA/LTA方法進(jìn)行了改進(jìn)并用來進(jìn)行微震事件的檢測和震相到時拾取。

Allen[4]改進(jìn)的特征函數(shù)為：

(5)

其中：X(i)為原始微震事件在第i個采樣點的輸入信號；k為權(quán)重因子，是根據(jù)輸入信號X(i)的采樣頻率和噪聲特性，對信號振幅和頻率的權(quán)重進(jìn)行分配：

(6)

劉曉明等[8]在公式(5)的基礎(chǔ)上進(jìn)行改進(jìn)，計算整個微震事件序列X(i)的權(quán)值。改進(jìn)后的公式為：

(7)

其中，

(8)

文獻(xiàn)[4，8]根據(jù)輸入信號X(i)的采樣頻率和噪聲特性，引入權(quán)重因子k對信號振幅和頻率的權(quán)重進(jìn)行分配，構(gòu)建了一種新的特征函數(shù)，克服了傳統(tǒng)方法不能適應(yīng)復(fù)雜多變的微震信號到時拾取的缺點。然而，這些方法在低信噪比的微震數(shù)據(jù)應(yīng)用中，由于輸入信號X(i)受到噪聲的影響較大，往往得不到準(zhǔn)確的到時結(jié)果。

本方法根據(jù)信號包含的能量、振幅、頻率變化等特性，結(jié)合特征函數(shù)CF2，增強真實信號與噪聲的信噪比，以微震時間序列的振幅平方X(i)2，即特征函數(shù)CF2作為輸入信號。鑒于Allen[4]、劉曉明等[8]的權(quán)值構(gòu)建思想，構(gòu)建了一種新的特征函數(shù)：

(9)

壓裂進(jìn)行時，施工噪音較大，壓裂產(chǎn)生的微地震事件通常被強的背景噪聲淹沒，遠(yuǎn)高于系統(tǒng)分辨微地震事件的水平，且難以從記錄波形中被人工識別出來[10]。因此，針對低信噪比的微震數(shù)據(jù)，在應(yīng)用STA/LTA算法前，使用Saad等[22]提出的MLoG濾波方法去除背景噪聲，用來提高微震數(shù)據(jù)的信噪比。MLoG算子濾波方法的表達(dá)式為：

(10)

式中：σ為高斯濾波器的標(biāo)準(zhǔn)偏差，n為高斯指數(shù)，C為濾波階數(shù)。

通過MLoG算子濾波器對背景噪聲進(jìn)行平滑處理，增強了微震數(shù)據(jù)細(xì)節(jié)；在保留原始數(shù)據(jù)特征的基礎(chǔ)上進(jìn)行去噪，使微震有效信號與背景噪聲之間的差值最大化，增強信噪比。利用本方法進(jìn)行微震P波到時拾取的基本流程如圖1所示。

圖1 微震P波到時拾取流程圖

2 模型驗證

本節(jié)首先通過模擬數(shù)據(jù)來驗證改進(jìn)的STA/LTA算法的自動拾取效果。模擬數(shù)據(jù)的采樣頻率設(shè)置為1 000 Hz，信噪比設(shè)置為20，選取的采樣點數(shù)為3 500。

圖2(a)為模擬微震數(shù)據(jù)記錄。對該數(shù)據(jù)分別采用特征函數(shù)CF2和CF4得到的STA/LTA進(jìn)行到時拾取，對應(yīng)的長時窗NLTA設(shè)置為50 ms，短時窗NSTA設(shè)置為20 ms。圖2(b)為圖2(a)模擬微震數(shù)據(jù)的特征函數(shù)CF2和CF4響應(yīng)結(jié)果；圖2(c)為對圖2(a)分別利用特征函數(shù)CF2和CF4進(jìn)行計算的STA/LTA到時拾取結(jié)果?？梢钥闯?，基于兩種特征函數(shù)的STA/LTA拾取的P波到時與模擬的數(shù)據(jù)記錄完全對應(yīng)，證明了改進(jìn)的特征函數(shù)CF4理論正確并且可行。

為驗證算法在弱信噪比情況下微震事件的識別能力，對圖2(a)中的模擬微震數(shù)據(jù)記錄添加了90 dB高斯白噪聲，見圖3(a)。圖3(b)為特征函數(shù)CF2和CF4對加噪后微震數(shù)據(jù)的響應(yīng)結(jié)果。圖3(c)為對應(yīng)的STA/LTA到時拾取結(jié)果。可以看出，在90 dB噪聲干擾下，特征函數(shù)CF4對應(yīng)的STA/LTA曲線在微震有效信號到來時刻的幅值變化更加明顯，對應(yīng)的P波到時與利用CF2得到的STA/LTA拾取結(jié)果一致，表明基于特征函數(shù)CF4的STA/LTA的微震到時拾取方法具有一定的抗噪性和準(zhǔn)確性。

為進(jìn)一步驗證算法的抗噪聲干擾能力，對圖2(a)中的合成數(shù)據(jù)記錄添加100 dB高斯白噪聲，見圖4(a)。圖4(b)為直接采用特征函數(shù)CF2和CF4對添加100 dB高斯白噪聲合成數(shù)據(jù)的響應(yīng)結(jié)果。對該合成數(shù)據(jù)采用特征函數(shù)CF2和CF4的STA/LTA進(jìn)行到時拾取時，無法準(zhǔn)確拾取P波到時，見圖4(c)。造成這一現(xiàn)象的原因是噪聲干擾強度太大，真實的微震信號被強噪聲淹沒。因此，為了有效解決這一問題，利用MLoG低通濾波方法對圖4(a)中的微震數(shù)據(jù)進(jìn)行濾波處理，相應(yīng)的參數(shù)σ為10,n為32，濾波結(jié)果見圖5(a)。與圖2(a)相比，MLoG濾波后的微震數(shù)據(jù)雖然降低了微震能量幅值，但有效保留了原始信號形態(tài)特征。圖5(b)為特征函數(shù)CF2和CF4對圖4(a)濾波后微震數(shù)據(jù)的響應(yīng)結(jié)果。圖5(c)為濾波后微震數(shù)據(jù)CF2和CF4的STA/LTA的到時拾取結(jié)果。對比發(fā)現(xiàn)，特征函數(shù)CF4的STA/LTA結(jié)果中只有一個比較明顯的峰值，與真實的微震到時一致。而特征函數(shù)CF2的STA/LTA結(jié)果中存在多個比較明顯的峰值，且只有一個記錄與真實的微震到時一致，其他均為錯誤微震事件。因此，針對強噪聲干擾下的微震數(shù)據(jù)，利用MLoG濾波方法可以保留原始微震數(shù)據(jù)形態(tài)特征，聯(lián)合改進(jìn)的特征函數(shù)CF4的STA/LTA方法能有效提高微震數(shù)據(jù)的抗噪性和到時拾取的準(zhǔn)確性。

圖2 合成數(shù)據(jù)記錄自動到時拾取結(jié)果

圖3 添加90 dB高斯白噪聲后的到時自動拾取結(jié)果

圖4 添加100 dB高斯白噪聲后的到時自動拾取結(jié)果

圖5 圖4(a)濾波后的到時自動拾取結(jié)果

3 實測數(shù)據(jù)應(yīng)用

某礦區(qū)因歷經(jīng)多期構(gòu)造運動，地質(zhì)條件復(fù)雜，斷裂發(fā)育。區(qū)內(nèi)發(fā)育有固鎮(zhèn)—長豐、光武—固鎮(zhèn)等正斷層，這些斷層具有良好的隔水作用，為煤層氣聚集儲存提供了有利條件。所選用的壓裂井是布置在該礦區(qū)祁南塊段上的一口煤層氣試采取芯井。為監(jiān)測該壓裂井附近水力壓裂裂縫發(fā)育、分布規(guī)律等情況，在該區(qū)域建立微地震監(jiān)測系統(tǒng)，記錄因連續(xù)壓裂引起的儲層物理特性改變而產(chǎn)生的微地震活動，為礦區(qū)地面水力壓裂工程成果的考察檢驗提供指導(dǎo)性建議。

微地震監(jiān)測系統(tǒng)的采樣頻率為1 000 Hz，采集到的三分量微地震監(jiān)測數(shù)據(jù)見圖6和圖7中黑線。利用人工拾取的某一個P波到時為16 259 ms。在實測數(shù)據(jù)應(yīng)用中，因輸入的是三分量微震數(shù)據(jù)，所以到時拾取結(jié)果取三分量微震數(shù)據(jù)的均值。針對礦區(qū)微震數(shù)據(jù)，首先利用MLoG濾波方法對其進(jìn)行濾波處理，對應(yīng)的σ為2.5,n為10，然后分別選取特征函數(shù)CF2和CF4計算STA/LTA，STA/LTA曲線的峰值記為微震P波到時。結(jié)果見圖6和圖7。

圖6 以CF2作為特征函數(shù)的STA/LTA監(jiān)測結(jié)果

圖7 以CF4作為特征函數(shù)的STA/LTA監(jiān)測結(jié)果

對比兩種結(jié)果可知，當(dāng)微震有效信號首次到來時，采用CF2和CF4作為特征函數(shù)計算的STA/LTA曲線幅值相應(yīng)地存在較大的變化，實現(xiàn)了對微地震有效事件的準(zhǔn)確識別。兩種特征函數(shù)對應(yīng)的STA/LTA在X、Y、Z分量上拾取的到時相同，分別為16 262、16 261、16 256 ms，計算3個分量到時的均值為16 259 ms，與人工拾取結(jié)果一致，證明了該方法在實際應(yīng)用中的精確性和有效性。

4 結(jié)論

本研究提出一種聯(lián)合改進(jìn)STA/LTA與MLoG算子的微震P波震相到時自動拾取方法，將其應(yīng)用于模型數(shù)據(jù)與實際微震數(shù)據(jù)進(jìn)行驗算，得出如下結(jié)論：

1) 在強背景噪聲影響下，本研究改進(jìn)的特征函數(shù)CF4采用高階冪級數(shù)放大微震事件的有效信息，提高了P波到時拾取結(jié)果的準(zhǔn)確性。

2) 針對弱信噪比的微震數(shù)據(jù)，在強噪聲干擾下引入MLoG濾波算子對背景噪聲進(jìn)行平滑處理，增強信噪比，且在去除隨機(jī)噪聲干擾的同時最大程度保留有效信號形態(tài)，聯(lián)合改進(jìn)的特征函數(shù)CF4的STA/LTA方法能有效提高微震數(shù)據(jù)的抗噪性。

3) 通過模型數(shù)據(jù)和某地區(qū)微震實測數(shù)據(jù)試算，證明了本方法相對于傳統(tǒng)P波到時自動拾取方法有較強的抗噪能力，具有更高的準(zhǔn)確性。