河流相储层是碎屑岩储层的典型代表, 约占我国近43%的碎屑岩地质储量, 此类储层通常由多期多级次的复合砂体构成, 控制着地下储层空间分布与渗流机理^[1]。地质历史时期, 由于河流或水体频繁的迁移、改道和截弯取直等作用, 因而储层多级次复合砂体单元叠置, 关系复杂, 同时砂体单元之间普遍发育泥质构型界面, 界面的低渗或阻隔特性阻碍了流体运移, 形成了储层的非均质特性^[2]。随着油田进入开发中后期, 产油量下降, 含水率升高, 剩余油的开发是进一步挖潜增产的关键。储层构型界面的阻隔特性容易造成空间上不同的剩余油富集区, 富集区内部具有相对高孔隙高渗透率的储集特性, 但是富集区之间流体被阻隔, 这种非连通或半连通状态导致后期注采效果不佳, 难以有效开采剩余油^[1]。因此, 开展储层单元的连通性分析, 尤其是储层构型界面的空间展布规律研究, 对于中后期的老油田二次开发、认识储层渗流机理、预测剩余油分布以及提高采收率都具有重要的指导意义。

剩余油分布预测可从储层内部结构剖析出发^[3-6], 如纵向上识别旋回层序, 平面上划分储层单元, 分析可能的韵律水淹特征(曲流河点坝与复合点坝、辫状河道与心滩和高含水期水淹规律等)。通常基于测井和地震数据来确定剩余油分布, 预测方法包括地震正演与岩石物理结合的方法^[7]、时移地震方法^[8]、基于开发数据的人工智能方法^[9]等。这些分析或预测方法的适用性与数据质量、储层非均质性和开发阶段等因素关系密切。

陆上油田开发井网通常相对密集, 结合探地雷达、地震与油藏动态等资料, 学者们发展了“井震结合、以井为主”的储层构型研究思路, 但这种思路不能满足海上油田“稀疏井网大井距”开发要求, 海上油气田井距通常为300~500 m, 这种尺度在横向上大于多数砂体构型单元的尺度, 故导致难以描述井间砂体单元和构型界面。

为满足海上油气田的精细表征砂体单元与构型界面的需求, 除测井信息之外, 还需要充分挖掘丰富的高品质三维地震数据信息, 研究如何利用地震信息构建砂体构型界面模型, 为进一步研究砂体连通性、储层渗流机理和剩余油空间分布奠定模型基础。本文以海上油气田复合砂体构型界面为研究对象, 从高品质三维地震数据出发, 提取与复合砂体构型界面相关的地震属性, 使用图形处理相关算法, 将构型界面抽象到界面建模, 为分析油水在储层复合砂体系统中的渗流机理和剩余油的空间分布奠定了技术基础。

1 复合砂体构型界面建模的理论基础 1.1 复合砂体构型与构型界面理论

复合砂体构型是三角洲和河流相砂体储层中结构单元的大小、形状、规模和单元间叠置关系的总称。构型界面中的“构型”(Architecture)概念源于第一届国际河流沉积学会议上ALLEN的报告, 该报告明确描述了河流相中河道与溢岸沉积的几何形态与内部组合关系^[10]。之后, MIALL^[11]通过分析美国科罗拉多高原侏罗-白垩纪河流相, 提出了河流相沉积单元分类分析方法。1994年, WILLIS等^[12]利用河流相砂体构型理论详细地解释了北巴基斯坦的中新世时期溢岸沉积单元。随后构型相关理论在国内成为研究热点, 并有学者进行了大量的河道复合砂体沉积构型的研究。曲流河沉积的构型研究是近年来的研究热点之一^{[1, 13]}, 胡光义等^[1]提出了一套完整的河流相复合砂体构型概念体系和表征方法, 对当今河流相构型单元和构型界面的研究具有重要的指导意义, 通过定义复合砂体的概念, 明确指出复合砂体的单元组合性和级次性, 并提出每一级次复合砂体均是由次一级砂体单元和构型界面组合而成。

1.2 砂体构型界面

复合砂体构型界面是多级次复合砂体单元中同级次单元间低渗透阻隔界面, 是储层中的“渗流屏障”。复合砂体可分为多个级次, 同级次砂体单元间为低渗透构型界面, 以表 1为例, 说明同级次砂体单元与构型界面的几何位置关系^[1-2]。如表 1所示, 复合砂体可分为4个代表性的级次: 复合河道带级、单河道带级、复合点坝级和单点坝砂级。而4个级次可进一步分为7个子类。以第1级次为例, 复合河道带级次的构型单元为单河道或河道砂, 单河道间的构型界面为泥质含量高的河道间(或废弃河道)沉积。从表 1可以看出, 构型界面在空间上具有面和线的特性, 这是本文对构型界面几何抽象化的基础。

为了描述构型单元与构型界面的相对关系以及流体在构型单元中的流动机理, 我们建立了界面地质模型。传统的构型界面建模方法是基于“以井为主”的研究思路, 先利用井网的测井资料和岩心等资料, 结合沉积模式、级次约束, 解剖构型, 之后使用确定性方法建模^[14-17]。本文着重提出了另一套思路, 即不受限于密井网资料的要求, 充分发挥海上油田高品质三维地震资料的优势, 突破传统上简单的沉积模式的限制, 利用地震数据预测地下真实的构型界面分布, 建立了基于地震数据的界面建模方法。

1.3 地震属性的“膨胀”、“腐蚀”与“骨架化”

砂体构型界面的抽象过程是基于数学形态学中的“膨胀”和“腐蚀”的图像处理方法。“膨胀”和“腐蚀”过程可以用图 1来描述, A经过B的“扩边”处理形成“膨胀”的结果(图 1a), 相反, A经过B的“削边”处理形成“腐蚀”的结果(图 1b)。理论上砂体构型单元内部地震波速度偏高, 而含泥界面地震波速度偏低, 速度差异导致砂体与界面间形成地震波阻抗差, 地震强反射特征明显, 平面上形成条带或块状的高亮属性。将地震强振幅属性转换成图像, 通过“膨胀”过程连接某些反射能量, 再通过“腐蚀”过程, 形成线状的几何特征(图 2)。通过进一步“骨架化”实现对块状属性的完全“线性化”, 提取线段信息, 将线段面垂直插入到油藏模型中, 实现含构型界面建模。

“膨胀”和“腐蚀”作为数学形态学图像处理的基本处理方法, 对图像特征提取和识别有着重要影响^[18]。通常, 结构元素(对应图 1中的B)用元素为1的矩阵表示, 同时标明结构元素的原点。“膨胀”过程中只要结构元素与图像(对应图 1中A)边缘有交集, 保留平移结构元素的原点即可实现该点的“膨胀”, 遍历全部图像区域, 实现图像目标区边缘的“膨胀”。同样, “腐蚀”过程则是当结构元素与图像目标区完全覆盖时, 才保留结构元素的原点位置, 这样非原点的某些位置被剔除, 即可实现目标区的“收缩”。其中, 如果二值图像表示为A, 进行运算的结构元素表示为B, 可定义如下操作过程。

A被B“膨胀”、“腐蚀”分别表示为:

$ A \oplus B = \left\{ {z|\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over B} } \right) \cap A \ne \phi } \right\} $

(1)

$ A\mathit{\Theta} B = \left\{ {z|{{\left( B \right)}_z} \cap {A^C} \ne \phi } \right\} $

(2)

式中: 集合A, B称为膨胀结构元素, 存在于所有形态学的操作中。${\left( {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over B} } \right)}$={c|c=-b, b∈B}, 表示对集合B做关于原点的映射; (B)_z={c|c=z+b, b∈B}, 表示集合B沿z方向的平移; A^C={a|a∉A}, 表示集合A的补集是不包含集合A的所有元素组成的集合。

在图像处理应用中, 更多地以各种组合来应用“膨胀”和“腐蚀”。其中最常用的3种方式为开操作, 闭操作, 击中或击不中变换。

A被B开操作表示为A$ \circ $B, 定义为A被B“腐蚀”, 然后再用B“膨胀”的结果:

$ A \circ B = \left( {A\mathit{\Theta }B} \right) \oplus B $

(3)

A被B闭操作则表示为, 是指先“膨胀”再“腐蚀”:

$ A \cdot B = \left( {A \oplus B} \right)\mathit{\Theta B} $

(4)

对于击中或击不中变换在匹配一副图像中像素的特定结构是非常有用的。A被B击中或击不中变换表示为$A \otimes B$:

$ A \otimes B = \left( {A\mathit{\Theta }{B_1}} \right) \cap \left( {{A^C}\mathit{\Theta }{B_2}} \right) $

(5)

数学形态学的“骨架化”过程, 主要是针对二值图而言, 将输入的具有一定宽度的图像轮廓用逐次去掉边缘的方法最终变为宽度仅为一个像素的骨架。设想在同一时刻, 将目标的边缘线都点燃, 火的前沿以匀速向内部蔓延, 当前沿相交时火焰熄灭, 火焰熄灭点的结合就是骨架。对于地震属性, 可使用骨架化的方法取得最终界限在平面内的线段特征。

1.4 霍夫变换

在砂体构型界面腐蚀形成“骨架”之后, 需要对图像的“骨架”数字化形成可存取的线段, 即记录下线段的两个端点, 利用这两个端点信息即可定义界面, 该过程可以通过人工拾取实现, 但如果“骨架”线段太多, 可使用霍夫变换的方法实现线段的自动拾取。

霍夫变换是一种寻找并连接图像中直线线段的方法。假设在直角坐标系内有一条直线定义为:

$ x\cos \theta + y\sin \theta = \rho $

(6)

我们通常认为θ和ρ为已知量, 水平直线的θ=0, ρ为正的x截距, 类似地, 当水平直线的θ=90°时, ρ为正的y截距。转换为法线方程后, 在参数平面内以正弦曲线显示上述两条直线。每一条正弦曲线表示通过特定点(x_i, y_i)的一族直线, 交叉点对应于通过不同点的线, 图 3展示了霍夫变换的变化过程。

霍夫变换将参数空间再细分为累加单元, 由此可以提升霍夫变换的精度, 其中的细分数目决定了这些点共线的精度。

需要注意的是, 霍夫变换对线段的自动拾取存在拾取精度不高, 且在一套直线上的两个线段不能被分别记录的不足。为解决霍夫变换对输出参数敏感的问题, 可参考文献[19]提出的方法, 手动修改霍夫变换生成结果, 使之更为合理。

2 复合砂体构型界面建模

为了研究含构型界面的曲流河油藏中流体的动态特征, 需要从地震数据中提取构型界面的几何信息。本文对构型界面的研究从最简单的垂直界面入手, 使得实现的构型界面在油藏模型中是垂直面, 从油藏的平面俯视图看构型界面是由一系列相连或断开的线段组成, 后续的研究将考虑从三维地震数据中更直接地提取符合实际界面产状信息, 如界面倾斜或弯曲等。

从地震数据上看, 在储层内, 砂体单元与构型界面由于岩性差异导致地震波速度差异明显, 形成高值反射系数, 进而形成地震强反射信号。除此之外, 导致信号增强的因素很多, 包括断层、裂缝、构型界面等。对于由非构型界面引起的强信号, 可以通过断层和裂缝解释将非界面的因素剔除。下面根据图 4介绍复合砂体构型界面的建模流程及其关键步骤。

1) 在三维地震成像数据体中确定目的层油藏时间窗口, 提取窗口内地震数据, 计算窗口内地震数据的平面属性, 此时由窗口内三维地震数据沿z方向融合转换成二维平面地震属性, 这里地震属性是指信号能量属性或一致性属性(Coherence)等。

2) 将二维平面地震属性转换成彩色图像, 便于实现图像处理, 转换后, 只保留和构型界面有关的地震属性信息, 如果图像中还有井位/断层/刻度等信息单元, 需提前剔除。

3) 将彩色图像转换成黑白二值图像。

4) 根据数学形态学的处理方法, 将黑白二值图像进行“膨胀、腐蚀和骨架化”处理, 形成骨架图像。

5) 使用霍夫变换和人工拾取相结合的方法, 对骨架线段进行线段数字化。

这里需要注意的是, 在得到腐蚀后的骨架图像之后, 需要通过人工解释的方法, 剔除可能由油藏边界、断层或裂缝导致的骨架线段。一旦数值化构型界面之后, 我们可以将垂直界面插入到油藏模型中, 设置油藏和构型界面的孔隙度、渗透率、饱和度等物性参数, 为分析油水在储层复合砂体系统中的渗流机理和剩余油的空间分布奠定基础。

3 实际应用

以渤海Q油田为例, 该油田正处于开发的中后期, 油水关系复杂, 井间储层预测砂体分布不确定性高, 砂体连通性关系不明, 需要将井震信息有效地集成到油藏模型中, 深入开展砂体单元的叠置关系和构型界面空间展布的研究, 分析油水渗流机理和剩余油的分布规律。渤海Q油田面积约10 km², 油藏深度为1.1 km, 属于典型的河流相沉积, 各级砂体单元(复合河道、单河道、复合点坝和单点坝)与构型界面发育广泛。

以渤海Q油田某油层R1单元为例, 具体说明构型界面抽象和建模过程。R1油层单元平均厚度26 m, 为一套开发单元, 砂体遍布整个工区(面积约10 km²), 连通性好, 但是总体生产中动静态矛盾突出, 注采收效情况复杂, 推断和普遍发育的构型泥质界面(泥质隔夹层)可能存在一定关系。图 5是R1单元对应砂岩储层内地震振幅能量的平面显示, 高亮属性表示砂层内部岩性差异导致的强振幅反射, 高亮位置和界面位置相关, 开发阶段注采位置与效果分析确定了这种相关性的存在, 尤其是某些位置相邻注采井井距小, 但是注采效果不明显, 而测井资料显示泥质界限的存在阻隔了流体流动, 水驱受限。从图 5可以看出几乎所有井都位于低亮地震属性区域, 进一步说明界面与地震高亮反射能量的相关性。

本文基于地震属性构建构型界面模型, 以简化的垂直界面刻画实际构型界面特征(后续研究将考虑更复杂的界面特征, 本文不做详细讨论), 以研究其几何属性和物性对剩余油开采的影响。界面建模的过程如下:

1) 首先, 提取地震振幅属性, R1单元沉积层序简单, 综合解释成果发现没有断层和裂缝发育, 我们认为高亮部分表示由砂体构型界面所引起的信号增强;

2) 剔除多余信息, 对图 5的地震振幅能量剔除外在因素, 得到图 6, 剔除的对象包括井位井名、标尺、背景图等;

3) 将图 6的彩色图像转换成图 7的黑白二值图像;

4) 将图 7的黑白二值图像“膨胀”后, 得到图 8, 局部断开的地震属性连通性得到了增强;

5) 对图 8进行“腐蚀”和“骨架化”处理, 得到图 9;

6) 对图 9进行自动(霍夫变换)和手动修改结合的线段提取, 得到如图 10所示的线段化的砂体构型界面平面分布结果。

至此, 我们已经完成从地震数据中提取砂体构型界面的几何信息的任务。接下来, 我们将这些信息插入R1油组单元的油藏模型(图 11), 实现构型界面建模, 为进一步的油水流动模拟和剩余油分布的研究提供了基础。将数值化的线段经过的油藏网格用于表示数值化后的构型界面, 构型界面油藏网格具有低孔隙度、低渗透率的物性, 在设置完所有网格物性后, 此模型可用于油藏流体数值模拟, 以及研究构型界面的存在对流体渗流和剩余油分布的影响。

为了研究构型界面对剩余油分布的影响, 对渤海Q油田某油层(对应图 11中间区域)的机理模型进行模拟试验。试验展示了两个机理模型(无界面存在), 并模拟出基于反九点法井网方案开采若干年后剩余油的空间分布特征。图 12中油层平均厚度为26 m, 其中蓝色网格代表构型界面, 在三维空间中构型界面为三维垂直网格面, 界面从油层界面的顶到底完全穿透, 渗透率为0, 即流体无法穿过界面, 但可以从界面两端绕过。在本部分的数值模拟研究中, 利用无构型界面(初始模型)和含构型界面(界面模型)两个模型进行对比研究, 通过对比分析两个模型在相同的注采条件下22年后的剩余油空间分布特征, 定量分析构型界面的发育如何影响剩余油的空间分布。

图 12a和图 12b中蓝色单元渗透率几乎为0, 表示界面单元, 相对图 12a来说, 图 12b的构型界面完全阻隔流体, 流体只能从无界面网格通过, 其它位置由于流体通道变小, 流体流动不畅, 界面网格的渗透率几乎为0, 但界面网格位置的孔隙度和初始网格相同。

图 12c和图 12d是两个模型油藏在注采22年之后的含油饱和度分布情况。对比发现, 开发22年后, 剩余油饱和度的分布差异显著。无构型界面的图 12c含油饱和度分布平缓, 从中间注水井到其它生产井过度平滑; 而含界面的图 12d含油饱和度受界面控制, 界面网格单元位置渗透率几乎为0, 并且含油饱和度因渗透率几乎为0不发生变化(同开采前初始模型的含油饱和度一致), 图 12d的红色网格表示高含油饱和度的网格, 其中多数是界面网格, 界面位置剩余油完全没有动用, 还有一些是受界面隔断的网格, 剩余油被阻隔, 动用程度较低。另外, 由于水驱油的通道多经过界面之间的可渗透区域, 故在相同的注采条件下, 图 12d多数非界面位置剩余油动用更多。

图 12e和12f对应于初始模型和含构型界面模型的剩余储量分布。由于界面的发育, 导致含构型界面模型的红色区域位置(图 12f)剩余储量依然较高, 动用程度低, 难以形成有利注采条件。

表 2是含构型界面模型数模开发生产动态数据表, 可以看出: 经过22年的注采, 相对无构型界面模型, 含构型界面模型下, 石油少采了1.8%, 少采了68 328.0 m³, 平均单井少采8 841.0 m³, 每口生产井每年少采1 105.0 m³; 假设2000年的含水率为0, 2008年, 含构型界面模型的含水率相对增加了6.5%, 2022年, 含水率相对增加了2.8%。

总的来说, 本试验表明, 构型界面的发育导致了开采效果明显变差, 受界面的遮挡, 剩余油难以实现有效开采, 且分布更加复杂, 含水率上升。

4 结论

相对于传统海上复合砂体构型界面建模多采用的沉积模式与多井分析方式, 本文提出了一种新的基于纯地震数据构型界面建模的方法。该方法从地震数据的平面属性出发, 引入“数学形态学”方法处理平面属性, 形成表征构型界面的直线线段数据, 然后将线段作为边界构建垂直插入油层的垂直界面, 最后实现含构型界面模型建模, 以供油藏模拟和剩余油分析与预测使用。这种方法的最大优点是, 利用三维地震数据的丰富空间信息实现界面的快速建模, 简化了构型界面的建模过程。该方法开辟了一种地震约束河流相构型界面油藏建模的新思路, 作为一种界面建模新的尝试, 对曲流河储层预测具有重要的启示和指导意义, 为后续研究含构型界面的油藏建模、流体流动机理和剩余油分布规律奠定了基础。

本文以渤海Q油田为例, 阐述了曲流河构型界面建模的详细过程, 并且给出了一个构型机理模型的数值模拟试验, 定量分析了非渗透性构型界面对剩余油分布的影响, 为本方法的应用提供了参考。

作为基于地震数据的首次构型界面建模尝试, 本文只讨论了三维构型界面垂直的情况, 尚未讨论实际油层中三维界面可能倾斜(或曲面)的情况, 后续的研究将探索直接从地震数据中提取更具代表性的倾斜界面(或曲面)的过程和方法, 实现更具实际参考意义的构型界面建模和剩余油空间分布机理研究。