利用單井微地震波形能量反演震源機制

趙 煒 辛 維 毛中華 翟 尚 胡天躍 何 川*

(①北京大學地球與空間科學學院石油與天然氣研究中心，北京 100871； ②中國科學院地質與地球物理研究所，北京 100029； ③中國石化石油工程地球物理有限公司勝利分公司，山東東營 257086)

1 引言

微地震監測作為非常規儲層開發中一種有效的監測手段，可以提供儲層改造過程中巖石破裂產生的微型地震震源位置、發震時刻、震級等大量有價值的信息[1-3]。參考天然地震震源機制相關理論可知，微地震事件的波形中還包含關于震源機制的信息，可以用于求取走向、傾向、滑動角等裂縫破裂面參數，對水力壓裂儲層改造效果的精細評價具有重要的意義[4-8]。

微地震觀測方式一般分為三類，即井下觀測、地面觀測以及井地聯合觀測。地面觀測常采用放射狀布設檢波器，覆蓋范圍廣、觀測方位角多、觀測張角大、信號特征往往較為完備[9]。但實際應用中受到波場傳播路徑以及地面條件的影響，微地震信號十分微弱，有效事件一般局限于少量能量強度相對較大、信噪比較高的事件，大量能量弱、信噪比低的事件被忽略，直接導致震源機制信息相對匱乏，儲層改造效果評價的有效性受到影響。與之相對應，井下觀測的微地震信號信噪比較高，事件數也往往遠多于地面觀測，是目前應用最為普遍的觀測方式。然而出于成本方面的考慮，井下通常采用單井觀測方式，檢波器的數量及空間覆蓋范圍都十分有限，導致觀測張角較小、方位角單一，為震源機制求解帶來極大的不確定性。

Brune[10]提出了基于扁平型裂紋的遠場地震滑動輻射模式理論，是計算地震事件標量地震級的經典理論。該理論適用于剪切型震源(DC)，經眾多學者補充和完善[11-15]，目前被廣泛應用于計算水力壓裂誘發微地震的震級，其局限性在于僅能求出微地震的標量地震矩，無法求解裂隙破裂面的傾向、走向和滑動角。

求解微地震震源機制，可以看作是利用有限樣本(即觀測資料)的某一類或某幾類特征(如振幅、波形、能量、偏振、頻譜等)對模式(即震源機制)進行有效識別[20]。K最近鄰(K-Nearest Neighbor，KNN)分類算法是一種理論上較為成熟的模式識別方法。該方法的思路是：如果一個樣本在特征空間中有k個最相似(特征空間中最鄰近)樣本，并且這些樣本中的大多數屬于某一類別，則該樣本也屬于這個類別。Bailey等[21]提出了一種加權K最近鄰方法，根據近鄰樣本到測試樣本的距離，將較大的權值賦給較近的近鄰，使算法的穩定性大大提高。該方法簡單，無需估計參數，無需訓練，適合對稀有事件分類。

單井監測情況下，現有的微地震震源機制求解方法在弱信噪比條件下往往效果不理想，亟待提出一種更為有效的震源機制求解方法。本文嘗試探討微地震事件P、S波能量與裂縫破裂面參數之間的聯系。以此為出發點，在計算得到震源位置及標量地震矩的前提下，采用微地震事件P、S波能量作為反演特征參數，利用加權K最近鄰算法建立目標函數，反演求解微地震震源機制。為了檢驗方法的效果，分別應用合成微地震數據和實際資料進行測試和分析。

2 方法原理

2.1 微地震震源機制表示方法及激發的地震波

為了描述不同的地震類型，將不同方向的力矩整合為一個張量，被稱之為地震矩張量[22]。地震矩張量可表示為

(1)

式中M0為標量地震矩。采用地震矩張量的9對力矩(圖1)的組合可以描述所有震源類型。

圖1 地震矩張量的9對力矩[23]

根據Stein等[24]關于震源模型的假設，本文將所討論的微震震源類型預設為DC型震源，將地震矩中待求解的6個獨立分量轉換為3個獨立的裂縫破裂面參數， DC震源與斷裂面參數的關系如圖2所示，其中D為滑動量，V為斷面的法向。在這一假設的基礎上， DC震源的地震矩張量可以由標量地震矩M0、 斷裂面走向φ、 傾角ψ、 滑動角γ計算[25]。在均勻各向同性介質中DC震源引起的地震波，Aki等[16]給出了相關的公式，經后人整理歸納，將點震源的遠場波場函數與震源機制的關系定義為[26，27]

(2)

(3)

圖2 斷裂面參數與地震矩張量的關系[25]

2.2 S/P波振幅比方法[28，29]

S/P波振幅比方法作為求解天然地震震源機制的主流方法之一，目前被廣泛應用于微地震震源機制求解。以采用S/P波振幅比方法求解震源機制為例，根據圖3所示[30]，將三分量檢波器記錄的遠場微震速度信號積分轉換成位移信號后，沿震源到檢波器路徑分解成徑向與切向分量，拾取P波與橫波(SH與SV波)的最大振幅AP、ASV、ASH，計算其振幅比

(4)

(5)

式中I是檢波器個數。求解一組參數φ、ψ、γ，使理論計算值與實際觀測振幅比的方差和最小，即為觀測記錄的震源機制解。

圖3 三分量檢波器數據旋轉到P波射線坐標系

采用S/P波振幅比方法在保留S/P波振幅比值信息同時，可以有效消除檢波器接收條件對波形振幅絕對值的影響，使反演過程滿足均勻介質假設，具有較強適用性。然而，該方法求取單一時刻振幅比值的同時忽略了震源的強度信息，在觀測方位角度單一、張角很小、有效信號微弱的單井微地震觀測情況下，很容易受到噪聲干擾，導致精度不高、穩定性較差[31]。

2.3 從震源機制角度對比波形能量與S/P波振幅比

微地震震源機制反演算法的關鍵是尋找對震源機制敏感的事件記錄屬性。由于井下觀測檢波器數量以及覆蓋范圍受到監測井軌跡的限制，可獲得的記錄數量有限，采用天然地震中常見的P波初動法求解井下震源機制會出現嚴重的多解性。本文嘗試探討微地震事件波形能量與裂縫破裂面參數之間的關系，并與S/P波振幅比進行對比。

假設地層模型為均勻各向同性介質，P波速度為4500m/s，S波速度為2600m/s，密度為2.5g/cm3，震源位于(0，0，2000m)，觀測系統如圖4所示：觀測井位于(300m，300m)，三分量檢波器共計13級，深度范圍為1680～1800m，間隔為10m，微地震記錄時間采樣間隔為0.25ms。

首先在觀測數據不含任何噪聲的情況下，通過分析圖5所示四種震源機制產生的模擬微地震事件能量與S/P波振幅比(圖6)對目標震源機制參數的收斂情況，對兩者進行對比分析。

圖4 單口觀測井與震源位置示意圖

圖5 不同震源機制的沙灘球

Li等[32]指出，在無其他信息的前提下，將震源機制的(φ、ψ、γ)的搜索范圍分別選為(0°，360°)、(0°，90°)、(0°，180°)，足以包含所有獨立解。根據目標震源機制，ψ角固定在45°，對φ角、γ角所構成的區域進行網格剖分，剖分間隔均選取為5°。依據式(3)在震源位置可以獲得2592(72×36)種不同的剪切震源的觀測記錄，再通過式(4)計算出震源機制模擬觀測記錄的S/P振幅比，并根據模擬記錄求得相應的能量。最后根據下式計算所有震源機制與目標震源模型模擬記錄的能量、S/P波振幅比之差

(6)

式中：E(φn,ψn,γn)和Q(φn,ψn,γn)分別代表搜索網格點處對應事件的能量與S/P波振幅比；E0和Q0為目標震源對應事件的能量與S/P波振幅比；n∈[1，N]，其中N=n1×n2×n3，n1為φ網格剖分個數，n2為ψ剖分個數，n3為γ剖分個數，每個網格點對應一種剪切震源機制。

圖6為模擬記錄的能量差、S/P波振幅比差(為了方便對比，將式(6)計算的差值進行了歸一化處理)隨震源機制變化的分布，二者均可以收斂到真實解。采用S/P波振幅比模型解與附近的網格點差值接近，無顯著差別，而采用波形能量，會在模型解附近一定范圍內形成收斂梯度，收斂性優于S/P波振幅比，因此可以認為在單井觀測情況下波形能量較S/P波振幅比對震源機制參數更敏感。

2.4 基于微震事件波形能量的震源機制反演算法

根據本文2.3中微地震事件波形能量對破裂面參數較為敏感這一現象，設計了一種基于微震事件波形能量的震源機制反演算法，具體實現流程如下。

(7)

式中F為正演算子。

(2)求出網格點處模擬波場記錄的波形能量。根據各個檢波器的P波與S波初至信息ti，c(c=1、2，分別表示P波與S波)，求出P波和S波的x、y、z分量在窗口內的能量Un,i,j,c

(8)

式中：ti，c為P波或S波的初至時間；lc為P波或S波的初至時窗長度。

(9)

式中：wi為信噪比權重； KNN為K最近鄰算法算子。

圖6 四個目標震源模型的波形能量差(上)和S/P波振幅比差(下)隨震源機制變化的分布

在實際研究中，不同檢波器上的信噪比可能由于距震源距離或輻射花樣的影響存在顯著差異，因此在確定震源機制時，有必要據此對各個檢波器的作用進行加權。第i個檢波器的權重wi由該道信噪比SNRi確定

(10)

相比常規S/P振幅比方法，本文設計的微震事件波形能量方法在保留微地震事件S/P振幅比信息的同時，通過求取時窗內微地震事件波形能量值，保留了震源的強度信息。在求解震源機制時，輸入參數從3個振幅比增加到6個能量值，增強了對反演目標的約束。此外，單井井中微地震觀測系統的觀測方位角單一、張角較小，所處儲層通常是比較穩定的層狀地層，各檢波器的接收條件差別不大，因此處理實際資料時，檢波器接收條件的差異對波形能量方法影響較小。另外，對時窗內各個分量的P、S波成分求取能量值，無需從波形中提取振幅比，一定程度上減輕了井中微地震觀測資料S波容易受到P波續至干擾、難以提取到合適振幅比的影響，擺脫了對單一時刻振幅比的依賴。

3 模型數據處理

應用無噪聲和加噪模型數據的處理結果，分析波形能量方法與S/P振幅比方法在觀測方向單一、張角較小的直井觀測系統下的適用性。地層模型與觀測系統與2.3中相同，在進行網格剖分時，三個角度的搜索間隔均選取為15°，依據式(7)在震源位置可以得到2184(N=24×7×13)種不同的剪切震源的觀測記錄。

3.1 無噪聲情況下兩種方法反演誤差對比

首先對無噪聲條件下的微地震模擬記錄進行震源機制反演研究。將所有的2184個不同震源機制對應的波場記錄視作待求解震源機制的微震事件，采用本文提出的加權K最近鄰波形能量方法進行震源機制反演。采用波形能量方法φ、γ的反演結果與模型參數存在細微誤差(反演誤差由反演結果與對應模型參數的角度差衡量)，且誤差為0的比例均在97%以上，ψ則完全一致，與采用S/P波振幅比方法反演的誤差分布總體情況類似(圖7)。

采用波形能量方法，在不考慮噪聲情況下，出現誤判的震源機制個數約占全部2184個樣本總數的3%。從分布特征看，誤判主要分為兩類：一類如圖8a和圖8b所示，不同震源機制的ψ角均為0°時，φ角與γ角互為交換，這類震源有著相同的波形，約占全部誤判樣本個數的92%；另一類如圖8c和圖8d所示，不同震源機制的ψ角均為90°，φ角與γ角中一個相同，另一個分別是0°與180°或者180°與0°，引起波形極性反轉，能量值相同，這種情況數量很少，約占模型測試中全部誤判個數的7%。

圖7 無噪聲情況下采用波形能量方法和S/P波振幅比方法的反演誤差分布

圖8 不同震源機制波形記錄的對比

不同震源機制波形相同或彼此極性反轉，引起波形能量值相同，導致波形能量方法求解相應震源機制時出現誤判，這兩種誤判情況分別是由于DC震源假設以及算法本身引起的。從前者看，DC震源假設導致震源機制求解本身具有多解性，S/P波振幅比方法同樣面臨這個問題，在實際情況中可以采用先驗信息約束。從算法本身看，基于提取能量值的震源機制反演方法導致極性信息缺失，無法從極性角度對求解進行約束。但總的來說，誤判情況占比很低(3%)，對本文所提出方法影響較小。

3.2 不同噪聲情況下兩種方法反演誤差對比

利用噪聲數據分別對波形能量方法與S/P波振幅比方法進行抗噪性分析。噪聲模型數據由式(7)模擬的地震波與不同幅度的隨機噪聲疊加而成。以事件S波初至為中心，向前、向后延拓一個短時窗，計算其信噪比。

(11)

式中：An表示隨機噪聲均方根振幅；As+n表示加噪后模擬地震記錄S波信號均方根振幅。

分別采用波形能量與S/P波振幅比的方法，求出不同信噪比條件下2184個震源機制的傾向、走向、滑動角，并計算反演誤差分布(反演誤差由反演結果與對應模型參數的角度差衡量)。當信噪比為40dB時，兩種方法的反演誤差分布類似，都保持著很高的精度，說明在高信噪比情況下，兩種方法的反演精度趨近(圖9a)。隨著信噪比下降，二者誤差開始上升，S/P波振幅比方法總體誤差上升尤為顯著(圖9b～圖9d)。比較而言，波形能量方法的反演誤差在各個信噪比條件下均低于S/P波振幅比方法，當信噪比下降到10dB時，波形能量方法誤差為0°區間的樣本比例仍大于50%，此時S/P波振幅比方法已經降到20%左右，且分布相對分散。可見，在同等噪聲情況下，波形能量方法反演精度優于S/P波振幅比方法。

圖9 不同信噪比下兩種方法反演的誤差分布

4 實際資料處理結果

實際資料來自于油田一口水平井第7壓裂段的微地震監測數據。觀測系統為一組15級井下檢波器，級間距為10m，時間采樣間隔為0.5ms。該觀測系統布設在距第7壓裂段射孔點平均距離624.34m處的一口直井中，2h15min的微地震記錄中共計獲得66個有效的微地震事件，逐道計算各道記錄的S波信噪比(圖10)，信噪比集中分布在0到15dB之間。有效微地震事件經過前期初至拾取、標量地震矩反演等處理，其震源定位結果如圖11所示。

圖10 實際資料S波信噪比分布

采用K最近鄰波形能量震源機制反演結果如圖12、圖13所示。分析圖12可見，裂縫走向集中分布在135°左右，與壓裂產生的主裂縫方向一致；裂縫傾角主要分布在15°到30°之間，為低角度縫；裂縫滑動角集中在90°附近。從沙灘球圖(圖13)可見，微震事件震源位置集中分布在距離射孔位置較近的東側，并且震源裂縫破裂面走向垂直于井壁，與主裂縫延伸方向基本一致，東側距離射孔點較遠的位置存在少量滑動破裂裂縫，可能是儲層壓裂改造后被再次激活的天然裂縫。射孔位置的西南側存在少量裂縫，裂縫破裂面走向大致順著井壁。

圖11 震源定位結果俯視圖

圖12 實際監測數據反演得到的震源走向(a)、傾角(b)和滑動角(c)分布

圖13 震源沙灘球分布俯視圖及其局部放大顯示

5 結束語

研究發現，微地震事件波形能量對裂縫破裂面參數較為敏感。在此基礎上，提出了一種利用波形能量改進型微地震震源機制反演方法。在單井觀測系統條件下，無論模型數據還是實際數據均取得了較好的效果，證明了該方法的有效性。下一步將研究增加極性約束，提高反演準確率。

感謝中石化地球物理公司勝利分公司于靜、馮剛為本文研究提供了實際監測資料，感謝北京大學地球與空間科學學院喻志超博士為論文提供了許多寶貴的建議。