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

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

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

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

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

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

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

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

(1)

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

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

CF1=|X(i)|，

(2)

(3)

(4)

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

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

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

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

(5)

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

(6)

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

(7)

其中，

(8)

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

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

(9)

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

(10)

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

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

圖1 微震P波到時(shí)拾取流程圖

2 模型驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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