在油气储量及产量中占据很大比重的碳酸盐岩和致密砂岩储层, 确定裂缝发育的区域、方位和密度具有重大意义^[1-2]。国外早在20世纪70年代就开始了利用地震方法开展裂缝检测的研究与应用, 其中包括利用叠后资料的构造曲率法、张量分析法以及多方位VSP检测技术等。90年代中期, 横波、转换波探测和多分量各向异性分析技术逐渐发展起来, 但是受成本及方法自身的限制, 这些技术都没有实现工业化应用。于是, 很多学者提出了以P波为主、S波为辅的裂缝检测技术, 在垂直定向裂缝的物理模型实验及数值模拟中证实了P波在通过地质裂缝体时表现出很强的方位各向异性特征^[3-4], 并利用该特性在裂缝检测方面取得了比较满意的成果。但是, 目前实际生产中基于该项技术的裂缝检测应用软件较少, 相关的文献资料也不多。个别软件可以利用叠前地震资料提取的信息研究裂缝的特征, 利用与裂缝相关的地质、地球物理、测井和岩石物理等信息实现对裂缝方位和密度的描述。还有个别软件可以为各向异性成像、裂缝检测和油藏描述提供解决方案。商业化的叠前裂缝检测软件还处于日趋成熟和不断完善的阶段, 比如对于不同类型数据、介质的支持方面尚需进一步完善, 在不满足叠前裂缝检测要求的窄方位观测系统下还不能提供更多的解决方案。另外, 在国际油价持续低迷、生产成本与日俱增的经济环境下, 已有的商业软件在采购价格、配套设备、技术支持和软件升级维护费用等方面也是需要考虑的现实问题。因此, 在已经具备成熟方法技术和软件研发能力的条件下, 提高生产效率, 降低生产成本, 自主研发一套功能完备、能够解决裂缝型油气藏勘探实际问题的软件产品很有必要。

1 方法原理及技术流程 1.1 叠前裂缝预测原理

P波在裂缝介质中传播时, 其振幅、频率、速度以及衰减等属性随着观测方位的变化而不同, 而这些变化与产生该变化的裂缝方向和强度相关。一般情况下, P波沿垂直裂缝方向的传播速度比沿平行裂缝方向的传播速度慢^[5-6]。我们提取这些方向和强度特征, 并将其投影到方位平面上, 即可获取当前区域的裂缝走势与强度。对于地震波传播的方向特征, 专家们利用不同的近似解简化了HTI介质P波入射时的AVO方程, 统一的简化算式如下:

$ {R_{\rm{P}}}\left( {\varphi , \theta } \right) = A\left( \theta \right){\rm{ + }}B\left( \theta \right)\cos \left( {2\varphi } \right) $

(1)

式中:R表示P波属性(反射系数、速度、频率等), θ为入射角, φ为相对于当前裂缝方向的方位角, A表示均匀介质下的反射强度, B为固定偏移距下随方位变化的振幅调谐因子。裂缝介质中, 在固定入射角情况下, 地震属性R随方位角φ的变化而变化; 在固定方位角的情况下, 地震属性随偏移距或者入射角的变化而变化。于是, 多个方位的属性变化在平面上的投影呈现类似椭圆的趋势, 而非裂缝型介质(各向同性)则呈现圆形趋势。

根据上述特征, 对给定平面上的一组样本点, 以椭圆方程为模型进行拟合求解, 使得某一椭圆曲线尽量逼近这些样本点(图 1), 并以此椭圆为拟合解, 解析椭圆方程得到裂缝参数。

根据解析几何中椭圆一般方程与椭圆参数的关系, 可以计算椭圆几何参数:

$ A{x^2} + Bxy + C{y^2} + Dx + Ey + F = 0 $

(2a)

$ x = \frac{{BE - 2CD}}{{4AC - {B^2}}}\;\;y = \frac{{BD - 2AE}}{{4AC - {B^2}}} $

(2b)

$ \begin{array}{l} a = 2\sqrt {\frac{{ - 2F}}{{A + C - \sqrt {{B^2} + {{\left( {\frac{{A - C}}{F}} \right)}^2}} }}} \\ b = 2\sqrt {\frac{{ - 2F}}{{A + C + \sqrt {{B^2} + {{\left( {\frac{{A - C}}{F}} \right)}^2}} }}} \end{array} $

(2c)

$ \theta = \frac{1}{2}\arctan \frac{B}{{A - C}} $

(2d)

式中:A, B, C, D, E, F为椭圆一般方程参数; x, y为椭圆中心相对坐标; a, b为椭圆长短轴大小; θ为椭圆长轴与正北方向的夹角, 即裂缝的方向; 长短轴之比a/b则反映了裂缝发育的强度(密度)。

1.2 窄方位观测系统优化分析 1.2.1 超道集面元组合分析

方位各向异性椭圆拟合的基础是全方位或者宽方位观测系统采集。对于窄方位观测系统, 在排列横纵比小于0.5的情况下, 不能保证所有方位的数据能够完整采集, 使得椭圆不能获得足够的方位数据信息, 从而影响了拟合参数的准确性。因此, 我们借助于超道集面元统计方法, 将一部分相邻的窄方位面元作为某个宽方位面元的一个方位, 对它们进行组合统计, 使之形成较大面积的宽方位面元, 即超道集面元, 从而获得可以满足叠前裂缝检测需求的覆盖次数和偏移距范围。考虑到工区目的层构造起伏大, 气藏横向变化等因素, 超道集面元的尺寸不宜过大, 需要根据实际情况下单点网格的大小作动态调整, 以便精确反映局部裂缝走向, 提高计算机实时检索的效率。图 2a是单点(未组合)CMP道集, 图 2b是向周边扩展3×3共计9个面元组合后的CMP叠前超道集, 可以看出, 道集组合后, 方位角和偏移距范围都有拓宽, 横纵向各方位的覆盖次数也有明显增加。图 3是组合前后的偏移距-方位角示意图, 可以看出组合前后的覆盖次数均匀性有明显差别。

1.2.2 方位不均匀自动分组分析

在窄方位观测系统条件下, 采用均匀划分方位的方法不能使各个方位的覆盖次数均匀。虽然通过上述超道集组合方式扩展面元范围能在一定程度上达到扩大方位角、偏移距并均衡覆盖次数的目的, 但是从图 3b可以看出, 组合范围有限, 不能完全满足叠前裂缝检测要求。针对该问题, 我们研发了方位不均匀自动分组分析功能, 在超道集面元的基础上, 通过人工调整方位扇区大小和方向, 结合偏移距范围和动态试验, 获得满足要求的划分方案。图 4为软件产品提供的自由方位划分界面。通过该界面的交互分析, 可以实时调整“方位扇区划分方案”表格中每个扇区对应的中心角度和扇区大小, 获取对应的覆盖次数; 也可以调整“参数”选项中的偏移距范围以及最大覆盖次数, 动态获取扇区大小。调整后的方位角、偏移距和对应的覆盖次数可以通过“方位角-偏移距交汇图”和“方位扇区划分图”实时展示出来。在“方位角-偏移距交汇图”中, 不同颜色的柱状图代表不同方位范围在不同偏移距范围的覆盖次数分布。“方位扇区划分图”是偏移距-方位角与覆盖次数在极坐标中的表现形式, 不同方位的扇区大小代表该方位范围对应偏移距区域的覆盖次数分布情况, 扇区的长度代表了该区域覆盖次数在总覆盖次数中所占的比例。从图 4a看出, 按照传统的均匀方位划分方案, 若固定扇区大小为30°, 中心角度为15°, 45°, 75°, 105°, 135°, 165°等6个方位, 该超道集各个方位的覆盖次数严重不均, 受窄方位观测系统设计影响, 75°与105°等方位覆盖次数远低于其他方位, 说明该划分方案明显不利于裂缝检测。

通过人工动态调整该区域(如75°与105°)(超道集)中心角度位置以及扇区大小, 可以使覆盖次数尽量分布均匀(图 4b), 在此基础上, 反复进行不同偏移距-方位角的覆盖次数分析(图 5), 为裂缝检测提供合理的优选参数。从图 5可以看出, 在近偏移距(391~652 m)范围内(红色实线框所示区域), 有效覆盖次数方位角范围在28°~48°与130°~155°区间内, 其他区域无覆盖, 且范围较小; 在远偏移距(3 100~5 600 m)范围内(绿色实线框所示区域), 虽然有效覆盖次数方位角范围比较广, 但在40°~140°区间无覆盖次数分布, 也不满足数据要求; 在中偏移距(800~3 000 m)范围内(黑色实线框所示区域), 有效覆盖次数方位角范围在4°~175°区间内, 且各方位分布比较均匀, 说明该偏移距范围(800~3 000 m)能够满足叠前预测对方位数据的要求。

1.3 叠前裂缝检测流程

目前方位叠前地震数据大致可分为CMP、方位CRP、方位角度道集、方位叠加道集等4种类型, 其中CMP道集中炮检坐标为实际数据, 方位信息可实时计算, 其他几种类型为固定方位信息。不同类型的数据需采用不同的处理流程, 如图 6所示。

首先通过常规地震保幅处理, 包括道编辑、带通滤波、剩余静校正、真振幅恢复、地表振幅一致性补偿及动静校等, 为裂缝检测提供优质的数据支持; 然后针对叠前道集同相轴不平的现象, 采用沿层自动追踪方法提取地震波形的形态特征, 在相似的波形上提取显著点作为层位点, 提高方位各向异性椭圆的指示精度。图 7a为沿时间线拾取的结果, 将方位振幅沿时间线拾取并投影到极坐标系后, 椭圆拟合解析的裂缝方向解析为148.84°; 图 7b为沿层自动追踪拾取的结果, 自动追踪有效拾取了目的层波峰数据, 裂缝方向解析为128.23°, 且解析的裂缝密度也不相同。另外, 针对不同类型的数据, 在处理流程中集成了不同的道集处理方法以满足叠前裂缝检测的要求。

1) 道集归一化, 用于对叠前方位道集和方位角度道集进行归一化处理。理论上各个方位的道集在任一时间点的零偏移距振幅值相同, 但实际上不同方位的道集在零偏移距时同一时间点处的振幅值并不收敛到相同值^[7-8](图 8a)。此时需要将超道集中所有方位道集提取的垂直入射地震道叠加, 计算叠加道的包络X_st, 并用叠加包络X_st的每个采样除以超道集面元中各个方位角度道集的垂直入射地震道包络X_env, 得到每个方位道集的归一化标定地震道X_nor, 即X_nor=X_st/X_env, 最后将每个方位的归一化标定地震道乘以该方位道集中的各道, 即可实现归一化。归一化后的各方位道集反映了垂直入射时来自不同方位的反射都趋于同一个反射点(图 8b)。

2) P, G拟合剖面。即对方位角度道集数据进行各个方位的P, G拟合。在入射角θ＜25°的条件下, 来自两个弹性介质之间的平面反射纵波振幅与sin(2θ)呈近似线性关系。

$ R\left( \theta \right){\rm{ = }}P{\rm{ + }}G \times \sin \left( {2\theta } \right) $

(3)

式中:P为线性方程的截距; G为方程的斜率或梯度。从某种意义上来说, 这种拟合得到的P波剖面比普通叠加剖面更接近于真正的零炮检距剖面, 或者说更接近于垂直入射的振幅估计^[9-11]。G值所组成的剖面称作梯度剖面, 它反映了各个方位反射振幅随入射角的变化率以及变化趋势。通常情况下, 垂直于裂缝方向的G值比平行于裂缝方向的G值变化更剧烈, 对裂缝的发育较为敏感。

3) 方位时差计算。P波从不同方位穿过裂缝介质时, 除了振幅变化外, 同一层位反射的子波走时亦有不同。在某个方位道集的同相轴上, 零偏移距到达时间与远道偏移距的到达时间之差被称为时差。若不同方位的时差为周期变化, 说明存在速度各向异性^[12-13], 可用其表征裂缝的发育程度和方位。由于各个方位子波的走时在垂直入射情况下都相等, 因此必须先将各方位零偏移距处的同相轴校正到同一个时间, 再计算方位的时差值。

4) 方位包络处理。反演过程中若是以P波振幅随方位角的变化来判断裂缝发育的强度和方位, 则所用的数据为各方位的入射角部分叠加振幅, 这种部分叠加道仍然存在正负值和过零点。袁刚等^[14]指出, 使用各方位叠加道的包络(瞬时振幅, 图 9蓝线所示)具有更为合适的能量意义, 且能够克服波形振幅的零点带来的纵向不连续现象。

2 软件架构设计 2.1 面向集成环境的框架结构

面向对象(OOP)的集成方法是近十几年来比较流行且实用的软件工程设计思路, 该方法以工作对象(object)为基本元素, 对数据与逻辑行为进行封装, 大量运用类实例对象模拟现实行为, 其代码具有一定的可重用性和可继承性, 在规模较小的软件体系中能发挥较强的控制作用^[15]。但是对于功能层次多、后续扩展需求强、模块(module)关系复杂的中大型软件工程, 纯面向对象的设计架构往往会缺乏对整体系统的详细划分, 出现对象设计结构不合理、各部分关系失调甚至对象继承无限放大的代码冗余现象, 使得后续维护升级困难。我们在软件架构的设计中使用面向集成环境的分层组件(component)模型, 在传统面向对象设计的基础上建立一种组件设计模式, 其兼具组件与模块的特色(图 10)。与面向对象的编程思想相比, 一个类也可以视为一个组件, 但更多情况下, 组件提供了更为高层的系统视角, 可以形象地描述如下:组件如同一片树林, 对象是树, 模块是某些相同种类树木的组合。基于构件的软件工程(CBSE), 也被称为基于组件的开发(CBD), 是一种强调技术重用的方法, 即将松散耦合的独立组件定义、实现和组合到系统中^[16]。一个独立的组件可以是一个软件包、web服务、系统资源或者封装了特定功能函数的模块。这样独立出来的组件可以单独维护和升级而不会影响到其他的组件。而模块化的目的在于将一个技术组件按功能做拆分, 形成相互独立的功能, 以便于每个模块只包含与其功能相关的内容, 模块之间通过接口调用。将一个大的系统模块化之后, 每个模块都可以被高度复用。组件与模块密不可分, 组件化的目的是为了解耦, 将系统拆分成多个组件, 分离组件边界和责任, 便于独立升级和维护。

在控制服务层面, 将已有的成熟控件(图形平台)作为新架构的核心集成环境, 与地震数据平台管理结合作为主服务组件, 形成一套基础的图形化地震数据服务中心。在应用层设计中, 将应用主控界面(用户界面)、应用框架(交互应用框架、批量模块应用框架)和应用服务(平面绘图框架、剖面绘图框架、算法库)设计为单独的集成组件, 利用核心集成环境提供的公共接口集成在平台上。在后台服务支持层面, 提供了跨平台的数据访问和文件I/O机制, 通过系统配置实现对多种数据库的支持; 封装了通用的数据库访问对象模型, 支持常用的数据格式文件解析和读写, 为应用软件提供统一的数据访问接口(图 11)。

2.2 模块分类与开发

在用户应用层面, 将算法、工具(交汇图、并行框架等)、交互视图封装成可以单独使用的功能模块, 具体如下。

1) 观测系统显示。全工区或选定区域炮检关系、覆盖次数分布、偏移距方位角统计。

2) 实时部分叠加。选定区域范围叠前道集分方位、分角度、分偏移距实时部分叠加。

3) 叠前叠后联动。叠前道集与叠后剖面的联合解释与联动分析, 实时方位划分等。

4) AVAZ正演模拟。分析地震振幅属性的AVAZ特征, 用以指示方位各向异性椭圆。

5) 叠前道集分析。常规的叠前道集处理。

6) 叠前裂缝反演。裂缝密度与方位参数反演。

7) 裂缝综合显示。全工区或选定区域裂缝信息统计与显示。

通过制定层次关系与调用次序, 各模块之间相互衔接, 内部降低耦合程度, 在继承共享平台服务的同时, 实现整个框架的集成组建。软件运行时, 各功能模块间除数据关联外, 一般不产生依赖性。每个模块都可与控制服务组件结合, 编译生成定制功能的执行程序独立运行, 这样不仅可以有效实现交互与批处理的分离, 避免大数据量运算过程中实时交互操作导致的不稳定性, 还可以在软件研发生产线上实现模块与组件的分工维护与测试, 减小软件回归测试的范围, 提高集成测试与错误修正的效率。

3 应用测试 3.1 工区概况

在南方某三维工区, 针对目标层开展了叠前方位各向异性裂缝预测与分析。该工区总面积约970 km², 整体构造为由南往北逐渐变缓的背斜, 越往北部构造越深越平缓。该工区资料采集观测系统排列纵横比约1.0：0.8, 是典型的宽方位观测系统, 为叠前裂缝检测提供了良好的数据支持。工区内已有2口探井A井和B井, 其中A井较深, 钻遇震旦系灯影组, B井钻遇寒武系筇竹寺组。该区目标层为寒武系筇竹寺组的页岩层, 依据钻井显示, 从上至下有1~4共4套页岩层。B井钻井结果表明, 本区主要目标层为1号和4号页岩层, 目标层段测试获气5×10⁴ m³/d以上, 表明该区具有较好的勘探潜力。

3.2 数据准备

该区前期已经对叠前地震数据做了相应的处理, 包括原始CMP道集动校正、去噪(面波/斜干扰/大值等)、地表一致性振幅补偿、超道集面元组合、分方位叠加、振幅补偿(反Q滤波)、去随机噪声等, 形成了从15°~145°(30°间隔)的6个方位叠加数据体(图 12)。数据准备主要包括:①对各方位叠加数据进行去异常值预处理, 消除野值带来的异常与压道现象; ②对目标层位不同区域做单点方位各向异性分析, 初步掌握全测网裂缝的发育与分布情况, 为后续的批量处理提供准确的反演参数; ③进行分方位振幅和波阻抗数据椭圆拟合, 批量计算裂缝属性。与此同时, 为验证反演结果的准确性, 采用叠后裂缝检测方法开展地层不连续性检测, 了解大中断裂与构造的展布特征, 最后结合叠后裂缝属性和叠前裂缝检测结果, 进行裂缝综合预测与对比。限于篇幅, 本文仅针对4号页岩层进行预测与分析。

3.3 单点联动分析

单点联动分析是指裂缝检测批量运算前利用图形用户交互联动功能, 实时结合叠前、叠后道集协助用户了解目标层位某CDP范围(局部)内细微裂缝的方位各向异性特征。单点联动分析有3种方法:①沿时间线分析, 将各个方位数据沿相同时间线方向提取并拟合椭圆; ②沿时间线自动追踪分析, 在一定时窗范围内对波峰或波谷进行自动搜索, 提取对应的波峰或波谷数据再拟合椭圆, 该方法在道集同相轴不平时能较好地改善分析质量; ③瞬时振幅包络分析, 通过提取各个方位振幅的包络数据, 拟合包络椭圆, 分析裂缝属性。

图 13为三种方法椭圆拟合的界面, 可以看出三种方法拟合结果基本相似, 但沿时间线自动追踪的拟合结果(图 13b)更接近于真实值, 各个方位的数据点基本都分布在椭圆线上。对应地, 沿时间线分析(图 13a)以及瞬时振幅包络分析(图 13c)两种方法拟合的椭圆均没有完全贴合每个数据点。当然, 也不是所有的数据都会出现上述情况, 不同的数据可以动态地选择不同的分析方法, 得到较好的拟合结果。从研究区南、北部随机抽取的4号页岩层方位拟合结果来看, 可以初步判断工区南部(图 14)裂缝发育程度比工区北部(图 15)强。

3.4 叠前裂缝检测 3.4.1 方位振幅AVAZ正演模拟

根据叠前地震资料属性数据的方位各向异性特征, 一般情况下沿裂缝走向方向振幅衰减慢, 而垂直裂缝走向方向振幅衰减快^[17]。但方位振幅不同于速度、频率、波阻抗及衰减等属性数据, 由于地层界面上下岩性及岩性组合的变化, 往往会出现与上述规律相反的情况, 其结果是方位各向异性椭圆在裂缝方位指示上出现了90°的偏差, 影响了解释成果的准确性。针对这种现象, 需要实时分析地下裂缝型储层的地震响应规律^[18], 利用反、透射系数与弹性常数关系的AVO简化方程^[19], 模拟地震波在已知裂缝模型中传播时各个方位的地震响应, 以确定椭圆长短轴的正确指示方向。

我们利用A井纵横波速度和密度数据与提取的子波数据, 结合分析层段的各向异性参数与角度信息, 模拟了HTI介质AVAZ地震道集(图 16)。图 16同时显示了正演模拟结果中各个方位角道集数据在某时间线上的振幅随入射角变化曲线以及20°入射角处各个方位振幅的正演模拟椭圆。可以看出, 沿裂缝方向(0)振幅随入射角增大而增大; 垂直于裂缝方向(90°)振幅随入射角增大而减小。模拟结果表明:对于该目标层, 地震振幅椭圆拟合的长轴代表了裂缝的走向, 裂缝方向解释为“长轴”。在裂缝属性计算时, 需选择椭圆长轴作为裂缝方向。若条件允许, 可以多做几口井的正演模拟实验来相互验证。

3.4.2 方位振幅裂缝属性反演

图 17为分方位地震叠加数据裂缝反演参数界面, 可以设置裂缝属性计算范围(Inline、Xline以及时窗)和计算参数(振幅类型、拟合方式、椭圆长轴指示), 利用沿时间线分析、沿时间线自动追踪及瞬时振幅包络三种方法进行裂缝密度和方向批量运算。

3.4.3 叠前裂缝综合分析

图 18a是用沿时间线分析方法计算得到的4号页岩层裂缝发育密度沿层切片, 可以看出, 裂缝发育密度分布具有一定的规律性, 在断裂发育部位, 裂缝分布更为集中, 南部裂缝较北部更为发育。图 18b是用沿时间线自动追踪方法计算得到的沿层切片, 可以看出, 裂缝发育密度分布与图 18a基本一致, 但规律性较好, 北部构造平缓区裂缝密度值较小, 裂缝不发育, 而南部由于受断裂的影响, 裂缝较发育。图 18c是用瞬时振幅包络方法计算得到的沿层切片, 可以看出, 裂缝发育密度分布与图 18b相似, 差异不明显。综合分析可知, 由于沿时间线自动追踪分析方法采用了沿波峰或波谷自动追踪的算法, 克服了各方位地震同相轴在时间上的微小差异, 使裂缝预测结果更加精确, 更符合地下实际地质情况。

为了更加直观地了解研究区断裂-裂缝发育情况, 本文结合叠后地震不连续性裂缝属性参数, 进行了叠前、叠后断裂-裂缝综合分析。根据地震勘探原理, 连续性较好的同相轴一般对应连续性较好的地层; 相反, 连续性较差的同相轴对应连续性较差的地质体, 如各种断层构造。利用叠后裂缝检测技术可以描述地震数据的空间连续性, 常用的方法有相干分析、曲率分析、倾角检测、应力分析等。这些方法通常以叠后地震数据为基础, 研究的是整个网状断裂系统中较大尺度骨架与轮廓的断裂(在平面上呈线性分布, 数量较少, 结构清晰), 也用于断层、河道、尖灭线的刻画。而对于较小尺度、数量更庞大的小型或微断裂, 利用叠后裂缝检测技术往往不能进行详细刻画, 但可以进行叠合对比, 验证叠前检测技术反演的全区大断裂带发育概貌与走势。图 19分别是研究区最大正曲率、叠前裂缝密度与发育方向以及三者的叠合显示结果。通过对比可以看出:叠前反演的裂缝发育密度和方向与最大正曲率所刻画的断裂-裂缝发育特征非常吻合, 其准确性与可靠性得到了进一步验证, 且在微裂缝的刻画上更为细致。

裂缝密度与方向的叠合显示不仅较好地展现了全区裂缝发育的相对强度, 而且能直观地显示裂缝发育的方向。对不同区域范围裂缝发育方向进行统计(图 20), 可以更细致地了解该区域内裂缝发育的优势方向, 为该区页岩气开发方案实施及井位部署提供准确的基础资料。

4 结束语

本文研究了叠前方位各向异性裂缝检测的方法和流程, 介绍了满足实际需求、具备自主知识产权的软件产品。该软件采用面向集成环境的组件式模型架构, 有效弥补了传统面向对象的软件设计思想在研发规模、开发效率、升级维护方面的不足, 使软件模块的可重用性得到进一步提高。实际应用及对比结果表明, 该软件可以为后续页岩气开发方案的部署提供准确的参考资料, 为裂缝型油气藏勘探开发提供有效的技术支持, 指导勘探井位的部署。通过进一步的扩展开发, 该软件可以为叠前弹性参数反演提供必要的数据支持, 满足油气储层叠前流体识别的需求。

致谢: 感谢中国石油化工股份有限公司西南油气分公司勘探院领导、专家在软件推广与应用测试中提供的大力支持和帮助！