高角度斷裂約束的方位傅里葉系數裂縫預測方法及在M氣田的應用

范廷恩 杜 昕* 馬淑芳 范洪軍 賀新蔚 樊鵬軍

(①中海油研究總院有限責任公司，北京 100020；②賽吉紀技術服務(北京)有限公司，北京 100020)

0 引言

渤海灣盆地M氣田是全球儲量最大的變質巖潛山凝析氣田，已探明天然氣地質儲量超千億立方米[1-2]。M氣田主要含氣層段發育太古界變質巖潛山裂縫型儲層，勘探階段探井鉆遇儲層品質差異大，存在復雜裂縫。在M氣田開發進程中發現，太古界潛山裂縫型儲層主要受內幕高角度斷裂控制，高角度斷裂區域普遍發育裂縫。因此，裂縫預測是開發階段儲層預測與井位部署面臨的核心問題。2019年，M氣田主體完成寬方位海底電纜(OBC)采集，覆蓋次數為1200，縱橫比達0.72，高品質的寬方位地震數據為疊前定量預測裂縫奠定了基礎。

與疊后裂縫預測方法相比，疊前裂縫預測方法具有定量預測裂縫密度與方向的優勢[3-4]，構建方位反射系數特征方程是后者的理論基礎。M氣田潛山的裂縫多呈近垂直高角度[5]，基于等效介質理論可將裂縫介質抽象為HTI介質[6]。構建HTI介質的方位反射系數近似方程主要有兩類方法：①基于傳統的HTI介質Ruger振幅方位各向異性反射系數近似方程；②基于楊氏模量表征的方位反射系數近似方程[7]。其中方法①應用廣泛[8-9]，并且是多種改進算法的基礎。然而，方法①預測的裂縫方向存在90°的不確定性，并且要求疊前道集的炮檢距與方位角信息豐富，計算效率較低，在一定程度影響了其在三維實際資料中的應用效果[10]。方位傅里葉系數裂縫預測方法[11-12]可緩解方法①的預測裂縫方向不確定性問題，同時解耦振幅隨炮檢距與振幅隨方位角的變化關系，并使用特定炮檢距的分方位疊加數據預測裂縫，顯著提高了計算效率，且預測精度與方法①相當。但是方法①完全基于地震數據，預測過程沒有考慮宏觀地質模式約束，預測結果的穩定性以及與構造背景的吻合程度有待提升。

本文改進了方位傅里葉系數裂縫預測方法。利用f-k濾波、高分辨率Radon變換從疊后地震數據中提取高角度斷裂信息，并映射為方位傅里葉系數裂縫預測中各采樣點的先驗權重，以此建立高角度斷裂約束的方位傅里葉系數裂縫預測反演目標泛函，最后求解得到裂縫密度與方向。新方法的裂縫預測結果與井上解釋的裂縫特征以及各井生產測試情況吻合較好。

1 方法原理

1.1 HTI介質振幅方位各向異性反射系數公式

近炮檢距Ruger縱波方位各向異性反射系數近似公式[9]為

R(φ,θ)=A+[Biso+Banicos2(φ-φsym)]sin2θ

(1)

式中：R(φ,θ)為隨方位角φ與平均入射角θ變化的縱波反射系數；A為振幅截距；Biso為各向同性梯度；Bani為各向異性梯度；φsym為觀測方位角。Hudson理論揭示，各向異性梯度與裂縫密度呈正相關關系，觀測方位則與裂縫走向垂直[13-15]。

方法①通常利用三角函數變換將式(1)改寫為cos2φ與sin2φ的線性組合函數，之后采用最小二乘反演得到裂縫密度與裂縫方位估計。由于式(1)中sin2θ項的影響，導致預測的φsym存在90°的不確定性。

1.2 方位傅里葉反射系數公式

基于式(1)，Shaw等[14]利用傅里葉級數展開給出了隨φsym變化的縱波反射系數公式

RPP(φ,θ)=r0(θ)+r2(θ)cos2(φ-φsym)+

r4(θ)cos4(φ-φsym)

(2)

其中

r0(θ)=A+Bsin2θ+Csin2θtan2θ

(3)

(4)

(5)

式(2)實現了振幅隨炮檢距變化關系(AVO)與振幅隨方位角變化關系(AVAZ)的解耦。其中r0(θ)等價表達了第三類AVO變化關系，r2(θ)及r4(θ)為AVAZ項cos2(φ-φsym)與cos4(φ-φsym)的系數。一般情況下，疊前地震道集的最大公共炮檢距存在上限，因此θ≤35°。在此條件下，式(4)中g(g-1)ΔδNtan2θsin2θ項以及高階項r4(θ)可忽略，則式(2)可近似為[15]

RPP(φ,θ)≈r0(θ)+r2(θ)cos2(φ-φsym)

(6)

其中

r2(θ)≈0.5Banisin2θ

(7)

式(6)為擬線性方程組，輸入特定入射角(炮檢距)的分方位疊前道集或部分方位角疊加數據，通過最小二乘反演求解得到Bani(表征裂縫密度)與φsym(表征裂縫方向)。

1.3 基于f-k濾波與高分辨率Radon變換的高角度斷裂信息提取

f-k濾波的理論基礎是二維傅里葉變換，常用于地震資料去噪[16-17]。核心原理為：時—空域的反射波與線性噪聲存在較明顯視速度差異；經過f-k變換反射波與線性噪聲在f-k域的斜率與分布區間顯著不同；在f-k域壓制噪聲，再經過反f-k變換得到去噪后的時—空域信號。

M氣田潛山內幕高角度斷裂在地震剖面上表現為大斜度(60°～70°)、近似呈線性反射形態的幾何學特征(圖1)，本質上反映了潛山內幕高角度斷裂與圍巖的顯著速度差異。因此，可利用對視速度差異敏感性強的f-k濾波從偏移成果數據中提取高角度斷裂信息。

圖1 M氣田典型地震剖面

首先利用二維模型測試f-k濾波提取效果。采用主頻為15Hz的Ricker子波與由模型(圖2a)計算的反射系數(密度為常數)褶積，得到合成地震剖面(圖2b)。利用f-k濾波方法從圖2b中提取的高角度斷裂信息(圖3a)與圖2a基本一致，但分辨率較低，同時存在假象(圖3a的紅框區域)。這是由于紅框區域潛山頂面地層反射傾角較大，并與潛山內幕高角度斷裂的共軛方向相近，通過設計濾波器形態難以有效去除f-k域中的地層響應所致。

圖2 潛山高角度斷裂二維模型

圖3 基于f-k濾波(a)與高分辨率Radon變換(b)提取的高角度斷裂

Radon變換廣泛應用于地震信號處理領域[18-19]，包括多次波壓制、地震反射同相軸識別、數據保幅重建等。三維Radon變換的實質是對輸入數據沿特定空間路徑積分，從而將具有規律排列的信號分解為Radon域內的稀疏散點，以此實現信號識別與分離。通過改變積分路徑，Radon變換能夠追蹤地震數據中不同形態的地層反射信息。三維線性Radon正變換為

(8)

其反變換為

(9)

式中m為三維地震數據體d(x，y，t)在τ-px-py域中的變換結果，x、y、t分別代表主測線、聯絡測線、時間，而px、py、τ分別為主測線方向斜率、聯絡測線方向斜率、時間截距。通過傅里葉變換，可將式(9)寫為矩陣形式

d=Lm

(10)

式中：d為t-x-y域三維地震數據；L為Radon變換算子；m為待求解的τ-px-py域Radon變換系數矩陣。基于最小二乘反演框架，通過正則化約束條件求解m。

高分辨率Radon變換應用稀疏正則化約束項替代傳統最小二乘約束項，使反演結果具有突出強信號、壓制弱信號的稀疏形態特點，較傳統Radon變換的分辨率更高[20]，適合精細提取高角度斷裂信息。高分辨率Radon變換目標泛函為

(11)

針對f-k濾波方法提取高角度斷裂信息存在假象的不足，應用三維高分辨率Radon變換方法從f-k濾波結果(圖3a)中進一步提取高角度斷裂信息，通過在τ-px-py域中壓制非高角度反射信號，然后將只保留高角度反射信號的τ-px-py域閾值結果利用反Radon變換轉換至t-x-y域(圖3b)，以此優化高角度斷裂提取結果。可見，與圖3a相比，圖3b的分辨率更高，同時能夠較好地分離高角度斷裂反射與地層反射。

根據M氣田裂縫發育程度主要受潛山內幕高角度斷裂控制的認識，將f-k濾波與高分辨率三維Radon變換提取的高角度斷裂信息映射為三維地震數據各采樣點參與AVAZ裂縫預測的先驗權重，使高角度斷裂發育區域的數據樣點對裂縫預測結果的影響更大。另一方面，加入先驗權重本身也是一種平滑機制，可提升裂縫預測結果的橫向穩定性。

權重計算公式為

(12)

式中：wx,y,t為各采樣點參與AVAZ裂縫預測的先驗權重；SGauss為高斯平滑算子；rx、ry、rt分別為主測線、聯絡測線、時間方向的權重計算半徑；ax,y,t為對dx,y,t以cthreshold作為門檻值的閾值結果，而dx,y,t為由f-k濾波與高分辨率Radon變換綜合提取的高角度斷裂信息。進一步對wx,y,t歸一化處理

(13)

1.4 基于高角度斷裂約束的方位傅里葉系數裂縫預測

將式(6)與式(13)結合，建立高角度斷裂約束的方位傅里葉系數裂縫預測公式

wR=w(r0+Fφr2)

(14)

式中：w為先驗權重；R為方位各向異性反射系數；r0為方位各向同性項；r2為方位各向異性項，代表裂縫發育密度信息；Fφ為方位角函數項cos2(φ-φsym)組成的角度矩陣算子。基于式(14)利用最小二乘反演建立高角度斷裂約束的方位傅里葉系數裂縫預測目標泛函

(15)

式中μ為最小二乘約束項的超參數。式(15)的解析解為

(16)

式中I為單位矩陣。

可通過IRSL、共軛梯度等算法求解式(16)。基于高角度斷裂約束的方位傅里葉系數裂縫預測技術流程如圖4所示。

圖4 基于高角度斷裂約束的方位傅里葉系數裂縫預測技術流程

2 實例驗證

M氣田位于渤中凹陷渤中19構造脊上，該構造脊同其西側的渤中13構造脊及東側的渤中21-22構造脊共同構成面積近440km2的古潛山群[1]，整體發育在太古界基底之上并受郯廬走滑斷裂帶后期切割改造，形成了大型復雜構造脊(圖5)。

圖5 M氣田潛山頂面裂縫預測先驗權重平面圖

M氣田裂縫型儲層受潛山內幕的高角度斷裂控制，沿高角度斷裂走向呈“漏斗狀”分布。準確預測高角度斷裂的空間展布是裂縫型儲層描述突破的關鍵。圖6為MM′地震剖面。可見：A、B井附近高角度斷裂發育程度一般，由A、B井的測井響應解釋的內幕帶儲層含氣性較差；C、D井附近的高角度斷裂較發育，與C、D井的測井響應解釋的儲層含氣性較好相吻合。

圖6 MM′地震剖面

圖7為基于f-k濾波提取的潛山內幕高角度斷裂信息。可見，f-k濾波結果分辨率較低，同時提取的高角度反射信號連續性較差。圖8為應用三維高分辨率Radon變換對圖7進一步提取的高角度斷裂信息。可見，提取效果得到一定改善，分辨率較高，高角度反射更連續。

圖7 基于f-k濾波提取的潛山內幕高角度斷裂信息

圖8 應用三維高分辨率Radon變換對

基于式(12)與式(13)計算各采樣點參與裂縫預測的先驗權重，得到高角度斷裂先驗約束權重剖面(圖9)。可見，高角度斷裂帶發育區域對應更大的約束權重值，對式(15)、式(16)的干預作用更強。此外，潛山頂面裂縫預測先驗權重平面圖(圖5)表明：D、C、A、B四口井中，C、D井先驗權重值較大，A、B井先驗權重值較小，這與圖6反映的高角度斷裂發育規律相吻合，驗證了先驗權重估計的合理性。

圖9 高角度斷裂先驗權重剖面

通過AVAZ特性分析以及對比不同入射角道集疊加地震剖面，優選裂縫預測的最佳部分入射角疊加數據。A井風化帶頂面AVAZ分析結果表明，原始疊前道集最大入射角約為35°，AVAZ特性隨著入射角增大而更明顯(圖10)。進一步采用不同入射角范圍數據疊加形成部分入射角疊加地震數據(圖11)，其中20°～35°疊加地震數據的高角度斷裂特征明顯(圖11c)。根據上述分析，確定以20°～35°疊加數據作為特定入射角(平均入射角為28°)的疊加地震數據，開展方位傅里葉系數裂縫預測。

圖10 A井風化帶頂面AVAZ特性分析(左)及30°入射角的分方位道集(右)

圖11 部分入射角道集疊加地震剖面

將20°～35°入射角道集資料以30°為間隔分為6組分方位角(0°～30°、30°～60°、60°～120°、120°～150°、150°～180°)數據并疊加，作為疊前裂縫預測的基礎輸入數據，基于式(15)實現高角度斷裂約束的方位傅里葉系數裂縫預測，得到裂縫密度與裂縫方向預測數據體。圖12為風化帶裂縫密度預測結果。可見：①與未考慮高角度斷裂約束的方法(圖12b)相比，考慮高角度斷裂約束的方法預測的裂縫密度(圖12a)與電成像測井解釋結果整體較吻合，均表現為A、C井區裂縫最發育，G、H井區次之。②圖12b的裂縫密度在研究區東側偏高，如F井(圖12黃框區域)的測井解釋裂縫密度(2.4)明顯低于C井(4.8)，但F井與C井附近的裂縫密度差異并不大。

圖12 風化帶裂縫密度預測結果

圖13為預測裂縫密度連井剖面。可見：①考慮高角度斷裂約束的方法預測的裂縫密度與測井流體解釋結果一致(圖13a)，井上解釋氣藏較發育部位的裂縫密度也較高；同時由A井生產測試解釋的干層對應較低的裂縫密度值(圖13a藍色方框)，結果較合理。②未考慮高角度斷裂約束的方法在G井與D井附近預測的裂縫密度相對偏低(圖13b)。圖14為考慮高角度斷裂約束的風化帶裂縫方向預測結果。可見，風化帶裂縫發育方向整體以近東西向為主，局部存在北東、南西向，這與電成像測井揭示的裂縫發育方向較吻合。

圖13 不同方法預測裂縫密度連井剖面

圖14 考慮高角度斷裂約束的風化帶裂縫方向預測結果

進一步利用已投產開發井驗證所提方法的裂縫預測結果(圖15)。目前各開發井鉆遇的儲層普遍存在裂縫，而開發井實鉆軌跡與預測裂縫密度疊合表明，各井第一靶點基本著陸于裂縫密度較高區域，說明預測結果較可靠。

圖15 開發井實鉆軌跡與預測裂縫密度疊合圖

將新方法獲得的裂縫預測成果應用于M氣田開發井井位設計與優化。圖16為設計井軌跡與預測裂縫密度、方向以及潛山頂面等T0圖的三維疊合顯示。可見，設計開發井靶點坐標位于裂縫密度較高區域，同時井軌跡與裂縫方向呈一定角度。后期將基于裂縫預測成果進一步建立M氣田潛山立體縫網，持續評價與優化設計井位及軌跡。

圖16 設計井軌跡與預測裂縫密度、方向以及潛山頂面等T0圖的三維疊合顯示

3 結束語

本文提出了基于高角度斷裂約束的方位傅里葉系數裂縫預測方法，考慮了裂縫沿著高角度斷裂呈帶狀分布的宏觀地質認識，增強了預測結果的地質意義。運用該方法預測M氣田太古界潛山裂縫型儲層，探井及開發井數據驗證了所提方法的可靠性，為具有類似高角度構造特征的油氣田的裂縫型儲層預測提供了新思路。