畢麗飛，曾志毅，張建中，黃忠來，芮擁軍，刁 瑞

（1.中國石油化工股份有限公司勝利油田分公司科技處，山東東營257001；2.海洋科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，中國海洋大學(xué)海洋地球科學(xué)學(xué)院，山東青島266100；3.中國石油化工股份有限公司勝利油田分公司物探研究院，山東東營257001）

水力壓裂微震監(jiān)測技術(shù)是一項(xiàng)通過觀測注水壓裂過程中所誘發(fā)的微震事件來監(jiān)測和評估裂縫發(fā)育情況的技術(shù)[1-4]。微震事件識別和到時(shí)拾取是微地震數(shù)據(jù)處理的關(guān)鍵環(huán)節(jié)，直接影響對裂隙成像的可靠性[5]。通常情況下，水力壓裂產(chǎn)生的微震事件能量十分微弱，微震觀測資料的信噪比較低。低信噪比資料的微震事件識別及其到時(shí)拾取是微震監(jiān)測技術(shù)的重要研究課題。

微震監(jiān)測記錄的數(shù)據(jù)量巨大，人工識別并拾取微震事件初至到時(shí)的工作量極大，且容易引入人為誤差。微震事件到時(shí)拾取是根據(jù)有效信號與背景噪聲在能量、偏振、頻譜和統(tǒng)計(jì)等方面的特性差異來實(shí)現(xiàn)的。常見的初至到時(shí)拾取方法基于單道記錄的算法，其中以長短時(shí)窗能量均值比（STA/LTA）算法[6-9]及其改進(jìn)方法為代表[10-12]。該類方法計(jì)算速度快，可以滿足實(shí)時(shí)處理要求，但是對低信噪比微震資料的應(yīng)用效果不佳。近年來，一些學(xué)者綜合利用微震信號與環(huán)境噪聲的多種特征差異，提高了識別和拾取效果。如：宋維琪等[13]利用小波分解與Akaike信息準(zhǔn)則相結(jié)合，有效壓制噪聲引起的局部效應(yīng)，提高微震事件初至拾取的準(zhǔn)確性；吳治濤等[14]聯(lián)合小波變換與偏振分析實(shí)現(xiàn)微震事件的全自動(dòng)識別和拾取P 波初至；譚玉陽等[15]通過綜合地震信號與環(huán)境噪聲在能量、偏振和統(tǒng)計(jì)方面存在的差異，提出一種微震事件初至拾取算法，有效提高了抗噪能力和拾取精度。

GIBBONS等[16]的研究表明，具有相似震源位置及破裂機(jī)理的地震事件通常在記錄上表現(xiàn)出相似的波形特征。一些學(xué)者發(fā)展了基于事件多道記錄相似特征的初至拾取算法[17-19]，如：譚玉陽等[20]通過多道相似系數(shù)實(shí)現(xiàn)微震事件的識別，避免了由于某些道數(shù)據(jù)的信噪比過低而遺漏微震事件的情況；VANDECAR 等[21]利用多道互相關(guān)與最小二乘法相結(jié)合，確定遠(yuǎn)震震相的相對到時(shí)，提升了拾取效率和抗噪性；DE MEERSMAN 等[22]和AKRAM 等[23]分別基于互相關(guān)的迭代算法優(yōu)化了微震事件的初至拾取結(jié)果；魏夢祎等[24]提出一種基于波形互相關(guān)的微震事件自動(dòng)拾取方法，通過疊加獲取高信噪比參考道，提高了初至拾取精度。在拾取微震事件的過程中，以上方法需要對各道微震事件記錄進(jìn)行疊加得到信噪比更高的疊加參考道，但微震事件信號的初至極性與震源機(jī)制相關(guān)，各道微震事件的初至極性會(huì)不一致，即初至極性具有正負(fù)之分[25-28]，致使通過波形直接疊加得到的疊加參考道的信噪比反而降低，影響了對微震事件的識別和到時(shí)拾取精度。

本文提出了一種適于初至極性反轉(zhuǎn)的微震到時(shí)拾取方法。通過計(jì)算任意兩道之間的延時(shí)互相關(guān)，獲得各道微地震的相對到時(shí)；利用相鄰道相乘后再疊加得到參考道，避免了道間波形相似系數(shù)（Semblance）和直接波形疊加得到參考道時(shí)無法有效克服極性反轉(zhuǎn)影響的困難。本文首先介紹了微震事件識別和初至拾取的方法原理和技術(shù)流程，然后分別利用合成數(shù)據(jù)和實(shí)際資料對方法的性能進(jìn)行了測試和分析。

1 方法原理

1.1 相對走時(shí)計(jì)算

取一時(shí)窗內(nèi)的微震記錄，利用互相關(guān)函數(shù)計(jì)算任意兩道微震記錄之間的互相關(guān)系數(shù)序列。時(shí)窗內(nèi)應(yīng)包含同一微震事件在各道記錄的P 波或S波初至，時(shí)窗長度應(yīng)大于微地震波的時(shí)長以及最短和最長偏移距檢波器之間的到時(shí)差，可根據(jù)射孔信號時(shí)長及最短與最長偏移距檢波器的到時(shí)差長度進(jìn)行估計(jì)。互相關(guān)系數(shù)計(jì)算公式為：

式中：xi和xj分別為地震道i和j的微震記錄；W為時(shí)窗長度；k為xi和xj的時(shí)間間隔；ci，j（k）為互相關(guān)系數(shù)?？紤]到初至極性存在反轉(zhuǎn)的情況，對互相關(guān)系數(shù)ci，j（k）取絕對值，其最大值cmaxij對應(yīng)的k為xi和xj的到時(shí)差，記為Δti，j。兩個(gè)地震道i和j包含的同一微震事件的到時(shí)差Δti，j可表示為：

式中：ti和tj分別為地震道i和j的微震初至到時(shí)。對所有地震道記錄進(jìn)行兩兩互相關(guān)處理，可得到兩兩記錄中微震事件的到時(shí)差。以4個(gè)地震道（i＝1，2，3，4）為例，可得到6個(gè)到時(shí)差Δt1，2，Δt1，3，…，Δt3，4。為了避免方程組求解過程的不適定性，增加一個(gè)約束方程∑ti＝0，使相對到時(shí)的均值為0。這樣有下列方程：

若參與計(jì)算的微震記錄有N道，則需利用（1）式進(jìn)行N（N－1）/2次互相關(guān)函數(shù)計(jì)算，得到N（N－1）/2個(gè)兩道間的到時(shí)差數(shù)據(jù)。（3）式寫成矩陣形式為：

式中：A為0，1和－1組成的稀疏系數(shù)矩陣；Δt為到時(shí)差向量；t為到時(shí)向量，也是待求向量。（4）式的最小二乘解為：

利用（5）式求出t后，再用t中的各個(gè)分量對各道記錄進(jìn)行時(shí)差校正。時(shí)差校正后，時(shí)窗內(nèi)各道記錄的微震初至波在時(shí)間上將趨于一致。

1.2 參考道獲取

利用（5）式求出的是相對于時(shí)差校正后各道記錄平均到時(shí)的相對到時(shí)。求微震事件在各道記錄上的絕對到時(shí)，還需要求出各道的平均到時(shí)，即微震事件在疊加道或參考道上的初至到時(shí)。求參考道及其微震事件的絕對到時(shí)的常用做法是：先將時(shí)窗內(nèi)時(shí)差校正后的各道記錄波形直接進(jìn)行疊加，形成高信噪比的參考道，再在參考道上拾取絕對到時(shí)。但是，當(dāng)時(shí)窗內(nèi)微震事件的初至極性存在反轉(zhuǎn)時(shí)，對校正后的各道微震記錄進(jìn)行直接疊加，會(huì)使微震記錄振幅相互抵消，反而降低了參考道記錄的信噪比，從而影響微震事件絕對到時(shí)的拾取精度。有的方法通過對各道的初至極性進(jìn)行識別或通過震源機(jī)制反演，將初至極性校正為相同極性后，再疊加得到參考道[23]，這樣做將會(huì)增加計(jì)算量。為了解決這個(gè)問題，我們利用相鄰道零延遲相乘再疊加的方法獲取參考道（圖1），計(jì)算式為：i＝1

式中：N是相鄰地震道數(shù)的配對數(shù)；W是時(shí)窗長度；w是時(shí)窗內(nèi)采樣點(diǎn)序號；x′i（w）表示時(shí)差校正后的地震道i在采樣點(diǎn)w處的幅值；X（w）為獲得的參考道在采樣點(diǎn)w處的幅值。圖1a中，時(shí)差校正后微震初至被“拉平”，其中左邊記錄（黑色地震道）與右邊記錄（紅色地震道）的初至極性相反。對時(shí)差校正后的微震記錄進(jìn)行相鄰道的零延遲相乘，初至極性相同（同為正或負(fù)極性）的相鄰道記錄的零延遲相乘后的幅值均為正值；而初至極性相反（一正一負(fù)）的情況僅出現(xiàn)在不同初至極性區(qū)域的分界處，如圖1b中第8道（紅色）所示。幅值為負(fù)的道的數(shù)量很少，對整個(gè)相乘道的疊加值影響很小，而隨機(jī)噪聲的極性無一定規(guī)律，疊加后振幅減弱，從而得到信噪比更高的參考道，如圖1c所示。

圖1 相鄰道零延遲相乘再疊加得到參考道示意a存在初至極性相反的地震道；b 相鄰道零延遲相乘道；c參考道

1.3 絕對到時(shí)拾取

由相鄰道零延遲相乘再疊加得到的參考道，雖然沒有完整的微震信號波形，但是在微震信號初至?xí)r刻附近有較大的振幅或能量變化，因此，本文利用STA/LTA 算法判斷參考道是否含有有效微震信號以及確定其初至到時(shí)。STA/LTA 算法的計(jì)算式為：

式中：X（w）為參考道記錄；WL和WS分別為長、短時(shí)窗的長度；m為采樣時(shí)刻；STA 和LTA 分別為短時(shí)窗和長時(shí)窗內(nèi)的能量均值。在初至到時(shí)附近，STA 比LTA 變化更快，得到時(shí)窗內(nèi)的STA/LTA會(huì)出現(xiàn)明顯的極大值，一般近似地將STA/LTA 曲線的極大值點(diǎn)作為初至到時(shí)點(diǎn)。若時(shí)窗內(nèi)的STA/LTA 的極大值超過所設(shè)閾值[9]，則認(rèn)為該時(shí)窗內(nèi)含微震事件，并且同時(shí)將STA/LTA 曲線的極大值點(diǎn)所對應(yīng)的時(shí)刻取為該微震事件在參考道的初至到時(shí)，記為T0；否則認(rèn)為該時(shí)窗內(nèi)不含微震事件。本文以各個(gè)時(shí)窗的長短窗能量比值的均值為依據(jù)，將該均值的某個(gè)倍數(shù)（以下稱為長短窗能量比值均值倍數(shù)）作為微震事件的識別閾值，這樣，就可在每個(gè)時(shí)窗內(nèi)自適應(yīng)地設(shè)置閾值，從而避免固定閾值所帶來的識別錯(cuò)誤和遺漏。

獲得參考道的初至到時(shí)T0，再結(jié)合（5）式中各道的相對到時(shí)ti，則各道微震事件的初至絕對到時(shí)Ti可由下式得到：

1.4 算法流程

本文的微震事件初至拾取算法流程如圖2所示，關(guān)鍵步驟概括如下。

1）確定滑動(dòng)時(shí)窗大小及滑動(dòng)步長，設(shè)置用于判別微震時(shí)間的長短窗能量比值均值的倍數(shù)b。

2）提取時(shí)窗內(nèi)的微震記錄，計(jì)算該時(shí)窗內(nèi)的道間互相關(guān)系數(shù)，確定微震事件的到時(shí)差Δti，j，并利用（5）式計(jì)算各道微震事件的相對到時(shí)ti。

3）利用得到的相對到時(shí)，對各道微震信號進(jìn)行時(shí)差校正，對時(shí)差校正后的微震記錄進(jìn)行相鄰道零延遲相乘再疊加處理，得到參考道。

4）利用STA/LTA 算法得到參考道的長短時(shí)窗能量比曲線，并計(jì)算能量比均值與設(shè)置的倍數(shù)b的乘積，即微震事件識別閾值R。判斷該時(shí)窗能量比極大值F是否大于微震事件識別閾值R，若是，進(jìn)行步驟5）；否則返回步驟2），進(jìn)行下個(gè)時(shí)窗的處理。

5）拾取STA/LTA 能量比曲線的極大值對應(yīng)時(shí)刻T0，再利用（8）式求取時(shí)窗內(nèi)各道微震事件初至的絕對到時(shí)Ti。

圖2 微震事件初至到時(shí)拾取算法流程

2 合成數(shù)據(jù)測試

首先利用合成數(shù)據(jù)來檢驗(yàn)本文方法的應(yīng)用效果。圖3 為合成微震記錄，時(shí)間長度為3s，采樣率為1 ms。圖3a為無噪聲微震記錄，含有3個(gè)有效微震事件，從第1個(gè)到第3個(gè)微震事件的振幅不斷減小，其中第1個(gè)微震事件初至極性在第8～16道上發(fā)生反轉(zhuǎn)。圖3b為添加均值為0，標(biāo)準(zhǔn)差為0.4的高斯白噪的微震記錄，其中第3個(gè)微震事件信號已經(jīng)完全淹沒在噪聲中。分別采用常規(guī)方法和本文方法進(jìn)行測試，選用的滑動(dòng)時(shí)窗長度和滑動(dòng)步長分別為120個(gè)采樣點(diǎn)和15個(gè)采樣點(diǎn)，取各個(gè)時(shí)窗參考道的長短窗能量比均值的3.5倍，作為微震事件的識別閾值。

圖4a和圖4b分別為采用各道記錄波形線性疊加參考道和本文的相鄰道零延遲相乘疊加參考道得到的STA/LTA 能量比曲線（藍(lán)色曲線）及自適應(yīng)閾值曲線（紅色曲線）?？梢钥闯?，圖4a所示的微地震記錄直接疊加的STA/LTA 能量比曲線只有1個(gè)明顯的峰值（紅色箭頭所指），指示了第2個(gè)微震事件，而對能量最強(qiáng)但存在初至極性反轉(zhuǎn)的第1個(gè)微震事件和信噪比較低的第3個(gè)微震事件卻沒反映出來；圖4b所示的曲線有3 個(gè)明顯的峰值（紅色箭頭所指），3個(gè)微震事件都被識別出來?？梢?，本文方法不但適于初至極性存在反轉(zhuǎn)的情況，對低信噪比資料也有很強(qiáng)的適應(yīng)性。

圖3 合成微震記錄

圖4 采用不同方法識別參考道微震事件的STA/LTA 能量比曲線（藍(lán)色）及自適應(yīng)閾值曲線（紅色）

下面用本文方法對微震事件初至到時(shí)進(jìn)行拾取，關(guān)鍵步驟的處理結(jié)果如圖5 所示。圖5a 為根據(jù)圖3b中3個(gè)微震信號所在的時(shí)窗分別提取出來的經(jīng)時(shí)差校正后的含微震時(shí)間記錄，與理論微震事件的到時(shí)一致，說明微震事件相對到時(shí)計(jì)算的準(zhǔn)確性。圖5b是相鄰道零延遲相乘結(jié)果，在一定程度上壓制了不相干的背景噪聲。圖5c是對圖5b各道進(jìn)行疊加的結(jié)果，即參考道，具有明顯的振幅極值，消除了初至極性反轉(zhuǎn)對參考道的不利影響。圖5d 為利用STA/LTA 算法得到的參考道的STA/LTA 能量比曲線及其極大值對應(yīng)的參考道的初至到時(shí)T0（紅色虛線所示）。

圖5 用本文方法對微震事件初至到時(shí)拾取過程的相關(guān)結(jié)果（左、中、右分別對應(yīng)圖3中的第1、第2和第3個(gè)微震事件）

利用得到的各道相對到時(shí)ti和參考道初至到時(shí)T0并用（8）式計(jì)算出微震事件在各道的初至到時(shí)Ti。圖6為采用本文方法拾取微震事件初至到時(shí)的誤差（拾取到時(shí)與理論到時(shí)之差）。從圖6可以看出，整體拾取誤差較小，說明采用本文方法可以得到較為準(zhǔn)確的微震事件初至的絕對到時(shí)。同時(shí)我們采用基于波形線性疊加參考道的識別方法拾取第2個(gè)微震事件初至到時(shí)，其相關(guān)結(jié)果和拾取誤差如圖7所示，與采用本文方法得到的第2個(gè)微震事件初至到時(shí)拾取結(jié)果僅相差1 ms（1個(gè)采樣點(diǎn)長度），說明本文方法的拾取精度與常用拾取方法的拾取精度相當(dāng)。但是，采用常用方法只能識別出不存在極性反轉(zhuǎn)的第2個(gè)微震事件，未能識別出含有極性反轉(zhuǎn)的第1個(gè)微震事件和信噪比很低的第3個(gè)微震事件，而采用本文方法可以識別出這3個(gè)微震事件?？梢?，采用本文方法不僅可以較為準(zhǔn)確地識別和拾取低信噪比微震事件，且能較好地克服初至極性反轉(zhuǎn)所帶來的不利影響，具有更強(qiáng)的適用性。

圖6 采用本文方法拾取微震事件初至到時(shí)的誤差（拾取到時(shí)減去理論到時(shí)）

圖7 采用基于波形線性疊加參考道的識別方法拾取第2個(gè)微震事件到時(shí)的相關(guān)結(jié)果及拾取誤差

3 實(shí)際資料處理

圖8 實(shí)際野外檢波器布設(shè)位置

采用本文方法對某區(qū)實(shí)際地面微震監(jiān)測資料進(jìn)行處理和分析。野外布設(shè)42個(gè)檢波器，位置如圖8所示。圖9為一段經(jīng)過帶通濾波后的地面微震監(jiān)測記錄，時(shí)間長度為4 s，時(shí)間采樣間隔為1 ms。采用本文方法和常用方法處理該段實(shí)際資料時(shí)選取的參數(shù)為：滑動(dòng)時(shí)窗為200個(gè)采樣點(diǎn)，滑動(dòng)步長為15個(gè)采樣點(diǎn)，判別微震事件的閾值設(shè)置為該時(shí)窗內(nèi)長短窗能量比均值的3.5倍。圖10a和圖10b分別為采用微地震記錄波形線性疊加方法和本文方法得到的參考道STA/LTA 能量比曲線（藍(lán)色）及自適應(yīng)閾值曲線（紅色），可以看出，微地震記錄波形線性疊加參考道的STA/LTA 能量比曲線沒有出現(xiàn)明顯的極大值，而本文方法的STA/LTA 能量比曲線出現(xiàn)了2個(gè)超過識別閾值的明顯的峰值（圖10b中紅色箭頭所示），檢測出該段微震記錄中的2個(gè)微震事件。

圖9 野外實(shí)際微震監(jiān)測記錄

利用本文方法對該段微震記錄的2個(gè)微震事件的初至到時(shí)進(jìn)行了拾取。圖11為實(shí)際微震信號的拾取過程的相關(guān)結(jié)果。圖11a為時(shí)差校正后的微震信號記錄，可以看出，微震事件波形基本上被校正拉平，其中第1個(gè)微震事件中的第5～12道和第35～40道的微震初至極性發(fā)生了反轉(zhuǎn)，波形線性疊加得到的參考道振幅會(huì)相互抵消，這是利用常規(guī)波形線性疊加參考道方法無法識別出微震事件的主要原因。第2個(gè)微震事件記錄的信噪比較低，個(gè)別道初至同相軸未被完全拉平，也導(dǎo)致波形線性疊加得到的參考道振幅較小，漏識了該微震事件。本文方法利用相鄰道零延遲相乘再疊加得到參考道，突出了被拉平的各道初至同相軸的作用，受部分未拉平道的影響相對較小。圖11b為利用相鄰道零延遲相乘的結(jié)果，有效壓制了背景噪聲，提高了信噪比。圖11c為由相鄰道零延遲相乘再疊加形成的參考道。圖11d為利用STA/LTA 算法對參考道進(jìn)行處理得到的STA/LTA 能量比曲線，以及由該曲線的極大值點(diǎn)得到參考道的初至到時(shí)T0（圖中紅色虛線所示）。最后把各道的相對到時(shí)ti和參考道的絕對到時(shí)T0相加，得到了微震事件在各道的初至絕對到時(shí)Ti，如圖12所示（紅點(diǎn)為各道微震事件的初至到時(shí)時(shí)刻）。

圖10 采用不同方法得到的參考道STA/LTA 能量比曲線（藍(lán)色）及自適應(yīng)閾值曲線（紅色）

圖11 實(shí)際微震事件初至到時(shí)拾取過程的相關(guān)結(jié)果（左列對應(yīng)915～1 115 ms時(shí)窗記錄；右列對應(yīng)2 565～2 765 ms時(shí)窗記錄）

為進(jìn)一步說明本文方法對實(shí)際資料中微震事件的識別率及初至拾取精度，以震源掃描定位算法[28]得到微震事件數(shù)目作為標(biāo)準(zhǔn)，分別采用本文方法和常用的基于波形線性疊加的微震事件識別方法對同一段壓裂的微地震監(jiān)測數(shù)據(jù)進(jìn)行識別和拾取，結(jié)果如表1所示。從表1可以看出，采用本文方法的誤拾率和漏拾率都小于采用波形線性疊加參考道方法的結(jié)果，尤其是漏拾率遠(yuǎn)小于后者，說明采用本文方法能識別出更多的有效事件，具有更高的微震事件識別能力。

圖12 采用本文方法對某段監(jiān)測記錄的微震事件到時(shí)的拾取結(jié)果（紅點(diǎn)為初至到時(shí)時(shí)刻）

由于采用掃描定位方法確定的初至到時(shí)受到速度模型精度的影響，因此本文在采用掃描定位法確定到時(shí)的基礎(chǔ)上進(jìn)行人工拾取，再以人工拾取結(jié)果為準(zhǔn)確值，計(jì)算本文方法和波形線性疊加參考道方法拾取微震初至到時(shí)的誤差。波形線性疊加參考道方法共識別出26個(gè)微震事件，這里分別計(jì)算了采用本文方法和波形線性疊加參考道法拾取這26個(gè)微震事件的到時(shí)誤差，統(tǒng)計(jì)的拾取誤差的頻率如圖13所示。從圖13可以看出，采用波形線性疊加參考道方法的拾取誤差主要分布在2～21 ms，而采用本文方法的拾取誤差主要集中在2～12 ms?？傮w上，采用本文方法的拾取誤差要低于采用波形線性疊加參考道方法的誤差。

表1 采用不同方法識別微震事件的結(jié)果

圖13 采用不同方法拾取初至到時(shí)的誤差率

4 結(jié)論

常用的多道相似系數(shù)法由于在識別有效微震事件時(shí)受初至極性反轉(zhuǎn)的影響，因而會(huì)導(dǎo)致微震事件的遺漏。本文根據(jù)微震事件的波形相似特征，并用相鄰道零延遲相乘再疊加的方法取代常用的波形線性疊加方法，獲取更高信噪比的參考道，提出了一種從多道微震監(jiān)測記錄中識別和拾取微震到時(shí)的方法。該方法不但能夠壓制不相干噪聲，從低信噪比監(jiān)測資料中獲取有效微震事件信息，而且能夠克服初至極性反轉(zhuǎn)的不利影響。合成數(shù)據(jù)和實(shí)際資料應(yīng)用結(jié)果表明，本文方法具有較強(qiáng)的微震事件識別能力和較高的微震到時(shí)的拾取精度，對低信噪比資料和存在極性反轉(zhuǎn)資料有較強(qiáng)的適應(yīng)性。

在判別參考道各時(shí)窗是否存在微震事件時(shí)，以各個(gè)時(shí)窗的長短窗能量比均值為依據(jù)，將該均值的某個(gè)倍數(shù)作為判別的閾值，這樣自適應(yīng)地設(shè)置閾值，避免了固定閾值帶來的識別誤差。

采用本文方法能較好地區(qū)別具有波形相似性的微震記錄與不規(guī)則噪聲干擾。但在實(shí)際資料中，除背景噪聲外，還可能包括一些相干噪聲，如反射波和折射波及續(xù)至波等，這些噪聲會(huì)對互相關(guān)時(shí)差估計(jì)和參考道產(chǎn)生影響，導(dǎo)致對微震事件識別困難及對到時(shí)拾取產(chǎn)生誤差，需要通過分析壓裂區(qū)域微震事件的初至曲線特征來進(jìn)一步研判和去除錯(cuò)誤的拾取結(jié)果。