高角度断裂约束的方位傅里叶系数裂缝预测方法及在M气田的应用

2022-12-09范廷恩马淑芳范洪军贺新蔚樊鹏军

石油地球物理勘探 2022年6期

范廷恩杜昕* 马淑芳范洪军贺新蔚樊鹏军

(①中海油研究总院有限责任公司，北京 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为考虑高角度断裂约束的风化带裂缝方向预测结果。可见，风化带裂缝发育方向整体以近东西向为主，局部存在北东、南西向，这与电成像测井揭示的裂缝发育方向较吻合。