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

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

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

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

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

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

1 改進的STA/LTA方法

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

(1)

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

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

CF1=|X(i)|，

(2)

(3)

(4)

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

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

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

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

(5)

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

(6)

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

(7)

其中，

(8)

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

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

(9)

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

(10)

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

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

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

2 模型驗證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結論

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

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

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

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