石油地球物理勘探 2022 , Vol. 57 Issue (6) : 1436-1444  DOI: 10.13810/j.cnki.issn.1000-7210.2022.06.019 范廷恩, 杜昕, 马淑芳, 范洪军, 贺新蔚, 樊鹏军. 高角度断裂约束的方位傅里叶系数裂缝预测方法及在M气田的应用. 石油地球物理勘探, 2022, 57(6): 1436-1444. DOI: 10.13810/j.cnki.issn.1000-7210.2022.06.019 . FAN Ting'en, DU Xin, MA Shufang, FAN Hongjun, HE Xinwei, FAN Pengjun. Application of high-angle-fault constrained azimuthal Fourier coefficient fracture prediction in M gas field. Oil Geophysical Prospecting, 2022, 57(6): 1436-1444. DOI: 10.13810/j.cnki.issn.1000-7210.2022.06.019 . 杜昕, 北京市朝阳区太阳宫南街6号院中海油研究总院有限责任公司开发地质研究室，100020。Email： duxin7@cnooc.com.cn 。 摘要 ：现有叠前宽方位裂缝预测方法基于数据驱动，地质模式认识无法干预预测过程；渤海湾盆地M气田太古界潜山裂缝性储层主要受内幕带高角度断裂控制。为此，基于M气田宽方位地震资料，提出基于高角度断裂约束的方位傅里叶系数叠前反演裂缝预测方法。首先，利用 f - k 滤波与高分辨率Radon变换从叠后地震资料中提取高角度断裂信息，进一步映射为方位傅里叶系数反演方法中的先验权重，以此建立高角度断裂约束的裂缝预测反演目标泛函，最终求解得到裂缝密度与方位角信息。所提方法在M气田的应用效果表明，裂缝预测结果更具地质意义，与井上解释的裂缝特征以及各井生产测试情况吻合较好，可用于开发井井位设计与优化。 Application of high-angle-fault constrained azimuthal Fourier coefficient fracture prediction in M gas field Abstract : The existing pre-stack wide-azimuth fracture prediction method is data-based, and the geological model cannot affect the prediction process. The Archaean buried-hill fractured reservoirs in the M gas field of Bohai Bay Basin (China) are mainly controlled by inside high-angle faults. To this end, based on the wide-azimuth seismic data of the M gas field, this paper proposes a high-angle-fault constrained azimuthal Fourier coefficient fracture prediction method with pre-stack inversion. Firstly, the f - k filtering and the Radon transformation with a high resolution are used to extract the high-angle fault information from post-stack seismic data, which is then converted as prior weights of the azimuthal Fourier coefficient inversion method, so as to develop the function of high-angle-fault constrained fracture prediction inversion targets and calculate the fracture density and azimuth estimation. The proposed method is applied in the M gas field, and the application effect shows that fracture prediction results are of geological significance and are in good agreement with fracture features interpreted from the electrical imaging well logging and the production and test data of each well. Therefore, the method can be applied to design and optimize the location of production wells.

渤海湾盆地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变换从叠后地震数据中提取高角度断裂信息，并映射为方位傅里叶系数裂缝预测中各采样点的先验权重，以此建立高角度断裂约束的方位傅里叶系数裂缝预测反演目标泛函，最后求解得到裂缝密度与方向。新方法的裂缝预测结果与井上解释的裂缝特征以及各井生产测试情况吻合较好。

R(\phi, \theta)=A+\left[B_{\mathrm{iso}}+B_{\mathrm{ani}} \cos ^2\left(\phi-\phi_{\mathrm{sym}}\right)\right] \sin ^2 \theta

式中： R ( ϕ , θ )为随方位角 ϕ 与平均入射角 θ 变化的纵波反射系数； A 为振幅截距； B _iso 为各向同性梯度； B _ani 为各向异性梯度； ϕ _sym 为观测方位角。Hudson理论揭示，各向异性梯度与裂缝密度呈正相关关系，观测方位则与裂缝走向垂直 ^{[ 13 - 15 ]} 。

方法①通常利用三角函数变换将式(1)改写为cos2 ϕ 与sin2 ϕ 的线性组合函数，之后采用最小二乘反演得到裂缝密度与裂缝方位估计。由于式(1)中sin ² θ 项的影响，导致预测的 ϕ _sym 存在90°的不确定性。

\begin{aligned} R_{\mathrm{PP}}(\phi, \theta)= & r_0(\theta)+r_2(\theta) \cos 2\left(\phi-\phi_{\mathrm{sym}}\right)+ \\ & r_4(\theta) \cos 4\left(\phi-\phi_{\mathrm{sym}}\right) \end{aligned} r_2(\theta)=\frac{1}{2}\left[B_{\mathrm{ani}}+g(g-1) \Delta \delta_{\mathrm{N}} \tan ^2 \theta\right] \sin ^2 \theta r_4(\theta)=\frac{1}{8} g^2\left(\Delta \delta_{\mathrm{T}}-g \Delta \delta_{\mathrm{N}}\right) \sin ^2 \theta \tan ^2 \theta

式中： B 、 C 为由Thomsen各向异性参数组成的各向异性系数；Δ δ _T 、Δ δ _N 分别为切向弱度与法向弱度； g = β ² / α ² ， α 与 β 分别为反射界面附近纵波速度 α 、横波速度 β 的平均值。

式(2)实现了振幅随炮检距变化关系(AVO)与振幅随方位角变化关系(AVAZ)的解耦。其中 r ₀ ( θ )等价表达了第三类AVO变化关系， r ₂ ( θ )及 r ₄ ( θ )为AVAZ项cos2( ϕ - ϕ _sym )与cos4( ϕ - ϕ _sym )的系数。一般情况下，叠前地震道集的最大公共炮检距存在上限，因此 θ ≤35°。在此条件下，式(4)中 g ( g -1)Δ δ _N tan ² θ sin ² θ 项以及高阶项 r ₄ ( θ )可忽略，则式(2)可近似为 ^{[ 15 ]}

R_{\mathrm{PP}}(\phi, \theta) \approx r_0(\theta)+r_2(\theta) \cos 2\left(\phi-\phi_{\mathrm{sym}}\right)

f - k 滤波的理论基础是二维傅里叶变换，常用于地震资料去噪 ^{[ 16 - 17 ]} 。核心原理为：时—空域的反射波与线性噪声存在较明显视速度差异；经过 f - k 变换反射波与线性噪声在 f - k 域的斜率与分布区间显著不同；在 f - k 域压制噪声，再经过反 f - k 变换得到去噪后的时—空域信号。

M气田潜山内幕高角度断裂在地震剖面上表现为大斜度(60°~70°)、近似呈线性反射形态的几何学特征( 图 1 )，本质上反映了潜山内幕高角度断裂与围岩的显著速度差异。因此，可利用对视速度差异敏感性强的 f - k 滤波从偏移成果数据中提取高角度断裂信息。

m\left(p_x, p_y, \tau\right)=\sum\limits_x \sum\limits_y d\left(x, y, t=\tau+p_x x+p_y y\right)

其反变换为

d(x, y, t)=\sum\limits_{p_x} \sum\limits_{p_y} m\left(p_x, p_y, \tau=t-p_x x-p_y y\right)

式中 m 为三维地震数据体 d ( x ， y ， t )在 τ - p _x - p _y 域中的变换结果， x 、 y 、 t 分别代表主测线、联络测线、时间，而 p _x 、 p _y 、 τ 分别为主测线方向斜率、联络测线方向斜率、时间截距。通过傅里叶变换，可将式(9)写为矩阵形式

式中： d 为 t - x - y 域三维地震数据； L 为Radon变换算子； m 为待求解的 τ - p _x - p _y 域Radon变换系数矩阵。基于最小二乘反演框架，通过正则化约束条件求解 m 。

高分辨率Radon变换应用稀疏正则化约束项替代传统最小二乘约束项，使反演结果具有突出强信号、压制弱信号的稀疏形态特点，较传统Radon变换的分辨率更高 ^{[ 20 ]} ，适合精细提取高角度断裂信息。高分辨率Radon变换目标泛函为

J(\boldsymbol{m})=\|\boldsymbol{d}-\boldsymbol{L} \boldsymbol{m}\|_2^2+\lambda\|\boldsymbol{m}\|_1

式中：‖ d － Lm ‖ ₂ ² 为最小二乘数据误差项；‖ m ‖ ₁ 为正则化约束项； λ 为约束项权重超参数。对式(11)的最小化问题一般采用迭代重加权最小二乘算法(IRSL)或交替方向乘子法(ADMM)等最优化算法求解 ^{[ 21 ]} 。

针对 f - k 滤波方法提取高角度断裂信息存在假象的不足，应用三维高分辨率Radon变换方法从 f - k 滤波结果( 图 3a )中进一步提取高角度断裂信息，通过在 τ - p _x - p _y 域中压制非高角度反射信号，然后将只保留高角度反射信号的 τ - p _x - p _y 域阈值结果利用反Radon变换转换至 t - x - y 域( 图 3b )，以此优化高角度断裂提取结果。可见，与 图 3a 相比， 图 3b 的分辨率更高，同时能够较好地分离高角度断裂反射与地层反射。

根据M气田裂缝发育程度主要受潜山内幕高角度断裂控制的认识，将 f - k 滤波与高分辨率三维Radon变换提取的高角度断裂信息映射为三维地震数据各采样点参与AVAZ裂缝预测的先验权重，使高角度断裂发育区域的数据样点对裂缝预测结果的影响更大。另一方面，加入先验权重本身也是一种平滑机制，可提升裂缝预测结果的横向稳定性。

权重计算公式为

\left\{\begin{array}{l} w_{x, y, t}=S_{\mathrm{Gauss}}\left(\sum\limits_{x-r_x}^{x+r_x} \sum\limits_{y-r_y}^{y+r_y} \sum\limits_{t-r_t}^{t+r_t} a_{x, y, t}\right) \\ a_{x, y, t}= \begin{cases}d_{x, y, t} & d_{x, y, t} \geqslant c_{\text {threshold }} \\ 0 & d_{x, y, t}<c_{\text {threshold }}\end{cases} \end{array}\right.

式中： w _{x , y , t} 为各采样点参与AVAZ裂缝预测的先验权重； S _Gauss 为高斯平滑算子； r _x 、 r _y 、 r _t 分别为主测线、联络测线、时间方向的权重计算半径； a _{x , y , t} 为对 d _{x , y , t} 以 c _threshold 作为门槛值的阈值结果，而 d _{x , y , t} 为由 f - k 滤波与高分辨率Radon变换综合提取的高角度断裂信息。进一步对 w _{x , y , t} 归一化处理

\bar{w}_{x, y, t}=\frac{w_{x, y, t}-\min \left(w_{x, y, t}\right)}{\max \left(w_{x, y, t}\right)-\min \left(w_{x, y, t}\right)} w \boldsymbol{R}=\boldsymbol{w}\left(\boldsymbol{r}_0+\boldsymbol{F}_\phi \boldsymbol{r}_2\right)

式中： w 为先验权重； R 为方位各向异性反射系数； r ₀ 为方位各向同性项； r ₂ 为方位各向异性项，代表裂缝发育密度信息； F _ϕ 为方位角函数项cos2( ϕ - ϕ _sym )组成的角度矩阵算子。基于式(14)利用最小二乘反演建立高角度断裂约束的方位傅里叶系数裂缝预测目标泛函

J\left(\boldsymbol{r}_2\right)=\left\|\boldsymbol{w}\left(\boldsymbol{R}-\boldsymbol{r}_0-\boldsymbol{F}_\phi \boldsymbol{r}_2\right)\right\|_2^2+\mu\left\|\boldsymbol{r}_2\right\|_2^2

式中 μ 为最小二乘约束项的超参数。式(15)的解析解为

\boldsymbol{r}_2=\left(\boldsymbol{F}_\phi^{\mathrm{T}} \boldsymbol{w}^{\mathrm{T}} \boldsymbol{w} \boldsymbol{F}{ }_\phi+\mu \boldsymbol{I}\right)^{-1} \boldsymbol{F}_\phi^{\mathrm{T}} \boldsymbol{w}^{\mathrm{T}} \boldsymbol{w}\left(\boldsymbol{R}-\boldsymbol{r}_0\right)

式中 I 为单位矩阵。

可通过IRSL、共轭梯度等算法求解式(16)。基于高角度断裂约束的方位傅里叶系数裂缝预测技术流程如 图 4 所示。

将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井附近的裂缝密度差异并不大。

徐长贵, 于海波, 王军, 等. 渤海海域渤中19-6大型凝析气田形成条件与成藏特征[J]. 石油勘探与开发, 2019, 46(1): 25-38.
XU Changgui, YU Haibo, WANG Jun, et al. Formation conditions and accumulation characteristics of Bozhong 19-6 large condensate gas field in oil offshore Bohai Bay Basin[J]. Petroleum Exploration and Development, 2019, 46(1): 25-38. 薛永安, 王德英. 渤海湾油型湖盆大型天然气藏形成条件与勘探方向[J]. 石油勘探与开发, 2020, 47(2): 260-271.
XUE Yong'an, WANG Deying. Formation conditions and exploration direction of large natural gas reservoirs in the oil-prone Bohai Bay Basin, East China[J]. Petroleum Exploration and Development, 2020, 47(2): 260-271. 姜晓宇, 张研, 甘利灯, 等. 花岗岩潜山裂缝地震预测技术[J]. 石油地球物理探, 2020, 55(3): 694-704.
JIANG Xiaoyu, ZHANG Yan, GAN Lideng, et al. Seismic techniques for predicting fractures in granite buried hills[J]. Oil Geophysical Prospecting, 2020, 55(3): 694-704. 陈怀震, 印兴耀, 高成国, 等. 基于各向异性岩石物理的缝隙流体因子AVAZ反演[J]. 地球物理学报, 2014, 57(3): 968-978.
CHEN Huaizhen, YIN Xingyao, GAO Chengguo, et al. AVAZ inversion for fluid factor based on fracture anisotropic rock physics theory[J]. Chinese Journal of Geophysics, 2014, 57(3): 968-978. 张志军, 肖广锐, 李尧. 渤中19-6油田变质岩潜山内幕裂缝地震响应特征及预测技术[J]. 石油地球物理勘探, 2021, 56(4): 845-852.
ZHANG Zhijun, XIAO Guangrui, LI Yao. Seismic response characteristics and prediction of fractured inside metamorphic buried hill of Bozhong 19-6 oilfield[J]. Oil Geophysical Prospecting, 2021, 56(4): 845-852. 王建花, 张金淼, 吴国忱. 宽方位杨氏模量反演和裂缝预测方法及应用[J]. 石油地球物理勘探, 2021, 56(3): 593-602.
WANG Jianhua, ZHANG jinmiao, WU Guochen. Wide-azimuth Young's modulus inversion and fracture prediction: An example of H structure in Bozhong sag[J]. Oil Geophysical Prospecting, 2021, 56(3): 593-602. 宗兆云, 印兴耀, 张峰, 等. 杨氏模量和泊松比反射系数近似方程及叠前地震反演[J]. 地球物理学报, 2012, 55(11): 3786-3794.
ZONG Zhaoyun, YIN Xingyao, ZHANG Feng, et al. Reflection coefficient equation and pre-stack seismic inversion with Young's modulus and Poisson ration[J]. Chinese Journal of Geophysics, 2012, 55(11): 3786-3794. 王康宁, 孙赞东, 侯昕晔. 基于傅里叶级数展开的纵波方位各向异性裂缝预测[J]. 石油物探, 2015, 54(6): 755-761.
WANG Kangning, SUN Zandong, HOU Xinye. Fracture prediction by P-wave azimuthal anisotropy based on Fourier series decomposition[J]. Geophysical Prospecting for Petroleum, 2015, 54(6): 755-761. 宋维琪, 徐月森, 张云银, 等. 基于傅里叶级数分析的各向异性参数估计及裂缝预测[J]. 石油物探, 2020, 59(1): 108-113.
SONG Weiqi, XU Yuesen, ZHANG Yunyin, et al. Anisotropy parameter estimation and fracture prediction based on Fourier series analysis[J]. Geophysical Prospecting for Petroleum, 2020, 59(1): 108-113. PAN X P, ZHANG G, YIN X Y. Amplitude variation with offset and azimuth inversion for fluid indicator and fracture weakness in an oil-bearing fractured re-servoir[J]. Geophysics, 2019, 84(3): N41-N53. ZHANG G Z, LI L, PAN X P, et al. Azimuthal Fourier coefficient inversion for horizontal and vertical fracture characterization in weakly orthorhombic media[J]. Geophysics, 2020, 85(6): C199-C214. RUGER A. P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry[J]. Geophysics, 1997, 62(3): 713-722. HUDSAON J A. Wave speeds and attenuation of elastic waves in material containing cracks[J]. Geophysical Journal of the Royal Astronomical Society, 1981, 64(1): 133-150. SHAW R K, SEN M K. Use of AVOA data to estimate fluid indicator in a vertically fractured medium[J]. Geophysics, 2006, 71(3): C15-C24. DOWNTON J, RUSSELL H. Azimuthal Fourier coefficients: A simple method to estimate fracture parameters[C]. SEG Technical Program Expanded Abstracts, 2011, 30: 269-273. 李庆忠. 地震信号内插与噪音剔除(二)[J]. 地球物理学报, 1988, 31(3): 329-341.
LI Qingzhong. Seismic signal interpolation and noise deletion(Ⅱ)[J]. Chinese Journal of Geophysics, 1988, 31(3): 329-341. 谭军, 宋鹏, 李金山, 等. 基于同相轴追踪的三维地震资料多次波压制方法[J]. 石油地球物理勘探, 2017, 52(5): 895-905.
TAN Jun, SONG Peng, LI Jinshan, et al. 3D multiple suppression based on event tracing[J]. Oil Geophysical Prospecting, 2017, 52(5): 895-905. SABBIONE J I, SACCHI M D. Restricted model domain time Radon transforms[J]. Geophysics, 2016, 81(6): A17-A21. GHOLAMI A, ZAND T. Three-parameter Radon transform based on shifted hyperbolas[J]. Geophy-sics, 2018, 83(1): V39-V48. 杜昕, 范廷恩, 董建华, 等. 基于多层感知机网络的薄储层预测[J]. 石油地球物理勘探, 2020, 55(6): 1178-1187.
DU Xin, FAN Ting'en, DONG Jianhua, et al. Characterization of thin sand reservoirs based on a multi-layer perceptron deep neural network[J]. Oil Geophysical Prospecting, 2020, 55(6): 1178-1187. DU X, LI G F, ZHANG M, et al. Multichannel band-controlled deconvolution based on a data-driven structural regularization[J]. Geophysics, 2018, 83(5): R401-R411.