随着油气勘探开发程度的不断深入, 勘探重点由构造油气藏向岩性油气藏甚至隐蔽性油气藏转变, 而复杂断裂储层作为主要隐蔽性油气藏类型之一, 微小断裂系统的精细识别是复杂储层勘探开发工程中的重要研究内容之一。与构造特征紧密相关的地震属性如构造倾角相干技术、多尺度分频相干技术和多曲率属性综合预测技术等取得了长足的发展, 在火山岩、碎屑岩等多种类型的复杂储层预测中发挥了重要作用^[1-3]。袁晓宇等^[4]采用蚂蚁追踪裂缝检测算法消除噪声对特定频段地震数据的影响, 可以检测不同尺度的断裂, 识别构造裂缝具有较高的精度。各向异性介质理论的发展为叠前裂缝预测提供了理论基础, 在理论模型和实际应用中均取得了一定的研究进展, 为复杂断裂储层的精细预测提供了有效的技术手段, 但是, 这些方法对测井岩石物理微观特征的利用程度不足^[5-7]。近年来, 基于深度学习的断层识别技术得到了快速发展, 通过卷积神经网络、感知归类、残差神经网络等多种方式实现断层面、断层走向等信息的提取, 实现了断层的自动识别与解释^[8-10]。

多属性结合的预测技术能够更加充分挖掘地震资料中丰富的储层发育信息, 建立储层参数与多属性之间的映射关系是利用多属性结合的有效途径。近年来, 人工智能(深度学习和机器学习等)算法在地震勘探中得到了广泛应用, 在岩相识别、属性分析、储层参数预测和油气检测等方面有了深入的发展^[11-13]。支持向量机学习算法(support vector machine, SVM)作为机器学习算法的重要组成部分, 在油气勘探方面, 已应用于岩性特征预测^[14]、烃类识别^[15]、关键储层参数预测^[16-21]。杜炳毅等^[22]从常规测井资料出发建立裂缝敏感曲线, 运用支持向量机算法识别断裂发育特征, 在小样本、非线性数据断裂预测中取得了较好的应用效果。

目前, 叠后复杂断裂储层预测仅仅利用了地震信息, 并未充分有效地利用测井岩石物理等有效信息进行约束, 导致复杂储层的断裂系统识别精度不高。针对上述问题, 我们在充分认识机器学习(支持向量机)理论优势的基础上, 提出了机器学习(支持向量机)的复杂断裂储层预测方法, 并在实际工区进行应用, 取得了理想的效果。

1 基本原理及方法

建立符合储层特征的各向异性岩石物理模型, 通过加入裂缝等复杂缝隙来逼近实际地层, 预测井位处的各向异性参数, 并计算微小断裂系统指示因子。同时, 提取构造类的地震属性(如相干、曲率等), 并进行优化处理, 结合测井敏感曲线优选与微小断裂相关的最优属性集合。选取已知井点处的微小断裂系统属性集和微小断裂系统指示因子作为样本进行机器学习(支持向量机)训练, 建立微小断裂系统属性集与微小断裂系统指示因子曲线之间的非线性映射关系, 最终实现复杂储层微小断裂系统特征的识别, 具体技术流程如图 1所示。

1.1 基于各向异性岩石物理模型的微小断裂系统指示因子曲线计算

本文利用各向异性岩石物理建模^[23]的方法实现对微小断裂系统指示因子曲线的求取, 从储层的本质出发, 充分考虑矿物组分、物性参数(孔隙度、饱和度、泥质含量等)、孔隙类型、各向异性等。

1) 利用Voigt-Reuss-Hill平均计算干岩石混合矿物的骨架模量。

首先, 充分考虑储层的矿物组分, 利用Voigt-Reuss-Hill模型构建干岩石骨架的模量。VOIGT^[24]在各矿物组分沿受力方向平均排列的假设条件下给出了各向同性完全弹性介质各组分平均应力与应变之比, 通过空间体积平均方法求取岩石的等效体积模量; REUSS^[25]假设各矿物组分垂直于受力方向的层状排列方式计算其岩石的等效弹性模量; HILL^[26]认为岩石弹性模量不超过Voigt上边界和Reuss下边界, 结合上下边界求取其算术平均值来近似表达岩石等效弹性模量, 即Voigt-Reuss-Hill平均。

2) 利用微分等效介质(DEM)模型引入孔隙类型, 计算含孔隙的干岩石骨架的模量。

DEM模型是将包含物逐步添加到背景介质中来模拟实际的双相介质, 等效岩石骨架和包含物的添加次序会影响等效介质的模量大小。

3) 利用Hudson模型和Schoenberg模型考虑各向异性介质的修正模量。

HUDSON^[27]模型假设弹性介质内部的裂缝呈硬币形的椭球缝, 其刚度矩阵表示为:

$ \boldsymbol{c}=\boldsymbol{c}_{b}-\frac{e}{\mu}\left(\begin{array}{cccccc} (\lambda+2 \mu)^{2} U_{11} & \lambda(\lambda+2 \mu) U_{11} & \lambda(\lambda+2 \mu) U_{11} & 0 & 0 & 0 \\ \lambda(\lambda+2 \mu) U_{11} & \lambda^{2} U_{11} & \lambda^{2} U_{11} & 0 & 0 & 0 \\ \lambda(\lambda+2 \mu) U_{11} & \lambda^{2} U_{11} & \lambda^{2} U_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu^{2} U_{33} & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu^{2} U_{33} \end{array}\right)+O\left(e^{2}\right) $

(1)

式中: c_b为各向同性背景岩石的弹性系数矩阵; λ和μ是拉梅常数; U₁₁和U₃₃是由边界条件计算出来的无量纲变量; e是裂缝密度; O(e²)表示裂缝密度的高阶项, 可以忽略。

为了得到裂缝介质的有效参数, SCHOENBERG等^[28]提出了线性滑动理论, 假设介质是无限薄的平面, 满足线性滑动边界, 在各向同性介质中嵌入一组平行的裂缝系统, 其刚度矩阵表示为:

$ \boldsymbol{c}=\boldsymbol{c}_{b}-\left(\begin{array}{cccccc} (\lambda+2 \mu) \varDelta_{\mathrm{N}} & \lambda \varDelta_{\mathrm{N}} & \lambda \varDelta_{\mathrm{N}} & 0 & 0 & 0 \\ \lambda \varDelta_{\mathrm{N}} & \frac{\lambda^{2}}{\lambda+2 \mu} \varDelta_{\mathrm{N}} & \frac{\lambda^{2}}{\lambda+2 \mu} \varDelta_{\mathrm{N}} & 0 & 0 & 0 \\ \lambda \varDelta_{\mathrm{N}} & \frac{\lambda^{2}}{\lambda+2 \mu} \varDelta_{\mathrm{N}} & \frac{\lambda^{2}}{\lambda+2 \mu} \varDelta_{\mathrm{N}} & 0 & 0 & 00 \\ 0 & 0 & 0 & 0 & \mu \varDelta_{\mathrm{T}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \varDelta_{\mathrm{T}} \end{array}\right) $

(2)

式中: Δ_N和Δ_T分别为正向差值和切向差值。

4) 利用Wood方程进行流体混合, 计算混合流体的体积模量。

当混合流体的尺度远远小于地震波长时, 可以用Wood方程计算出含流体混合物的速度。

5) 利用Brown-Korringa模型进行各向异性流体替换。

BROWN等^[29]考虑到各向异性对储层特征的影响, 在Gassmann方程流体替换的基础上提出了各向异性介质中的流体替换理论, 得到了饱和流体的各向异性岩石等效模量:

$ s_{i j k l}^{\mathrm{sat}}=s_{i j k l}^{\mathrm{dry}}-\frac{\left(s_{i j a a}^{\mathrm{dry}}-s_{i j a a}^{0}\right)\left(s_{b b k l}^{\mathrm{dry}}-s_{b b k l}^{0}\right)}{\varphi\left(\beta_{f 1}-\beta_{0}\right)+\left(s_{c c d d}^{\mathrm{dry}}-s_{c c d d}^{0}\right)} $

(3)

式中: s_ijkl^dry, s_ijkl^sat分别为干岩石和饱和岩石的有效弹性柔度; s_ijkl⁰为矿物的有效弹性柔度; s_ijaa^dry和s_ijaa⁰分别代表刚度矩阵中下标后两项相同干岩石和矿物的有效弹性柔度; s_bbkl^dry和s_bbkl⁰分别代表刚度矩阵中下标前两项相同的干岩石和矿物的有效弹性柔度; s_ccdd^dry和s_ccdd⁰分别代表刚度矩阵中对角线上干岩石和矿物的有效弹性柔度; β_f1和β₀分别为流体和矿物的压缩系数, β_f1=1/K_f1, β₀=1/K₀, K_f1, K₀分别为流体和矿物的体积模量; φ是孔隙度。

按照上述岩石物理建模流程实现对各向异性参数(即Thomsen参数ε^(V), δ^(V), γ)的计算, 进而估算储层各向异性梯度((4)式), 以此作为微小断裂系统指示因子曲线。

$ \varGamma=\frac{1}{2} \varDelta^{(V)}+4 k^{2} \gamma $

(4)

式中: k是横波速度与纵波速度之比。

1.2 地震属性提取与选择 1.2.1 地震属性提取

地震属性分析是断裂、岩性识别与流体检测的重要手段。地震属性是采用多种数学运算从叠后地震资料中求取的具有明确物理意义的特殊数据体, 能够表征地震波几何学、运动学及动力学特征, 能够反映构造、岩性、含油气性和断裂发育特征等地质信息^[30]。对构造特征敏感的地震属性有: 方差、相干、振幅能量梯度、倾角方位角、曲率、构造导向滤波、边缘检测、蚂蚁体、纹理和最大似然体等。

1.2.2 地震属性优选

利用多属性分析提高微小断裂系统识别准确性的关键在于确定最优的属性组合。丁峰等^[31]采用基于灵敏属性分析的相关聚类属性选择方法, 综合分析储层参数与地震属性之间相关性、有效性、符合率等参数, 从而得到最优的地震属性组合。我们利用微小断裂系统指示因子与属性间的相关系数以及属性之间的相关性后进行属性优选, 具体过程如下。

1) 计算微小断裂系统指示因子曲线与各属性间的相关系数向量。

假设有n个属性样本, p个属性变量, 计算微小断裂系统指示因子曲线与每个属性的相关系数并进行降序排列, 构成微小断裂系统指示因子与各属性间的相关系数向量, 公式如下:

$ \boldsymbol{R}_{\mathrm{wm}}=\left[R_{1}, R_{2}, \cdots, R_{p}\right] $

(5)

式中: R_j=$\left[\sum\limits_{i=1}^{n}\left(W_{i}-\bar{W}\right)\left(M_{i}^{j}-\bar{M}^{j}\right)\right]$/$\sqrt{\sum\limits_{i=1}^{n}\left(W_{i}-\bar{W}\right) \sum\limits_{i=1}^{n}\left(M_{i}^{j}-\bar{M}^{j}\right)}$(j=1, 2, …, p)表示第j个属性与微小断裂系统指示因子曲线之间的相关系数, 其中, W_i是微小断裂系统指示因子曲线的第i个样点值, W是微小断裂系统指示因子曲线的样本平均值, M_i^j表示第j个属性的第i个样点值, M^j是第j个属性的样本平均值。

2) 计算地震属性间的相关矩阵。

假设有m个属性样本, l个属性变量, 地震属性间的相关矩阵表示如下:

$ \boldsymbol{R}_{\mathrm{mm}}=\left[\begin{array}{cccc} R_{11} & R_{12} & \cdots & R_{1 l} \\ R_{21} & R_{22} & \cdots & R_{2 l} \\ \vdots & \vdots & & \vdots \\ R_{l 1} & R_{l 2} & \cdots & R_{l l} \end{array}\right] $

(6)

式中: R_ij表示第i个与第j个属性之间的相关系数, 有R_ij=$\left[\sum\limits_{k=1}^{m}\left(M_{k}^{i}-\bar{M}^{i}\right)\left(M_{k}^{j}-\bar{M}^{j}\right)\right]$/$\sqrt{\sum\limits_{k=1}^{m}\left(M_{k}^{i}-\bar{M}^{i}\right) \sum\limits_{k=1}^{m}\left(M_{k}^{j}-\bar{M}^{j}\right)}$, 其中, M_k^j表示第j个属性的第k个样点值, M^j是第j个属性的样本平均值。

3) 敏感属性相关聚类选择。

应用大量的冗余属性之前, 去除相关性较大的属性, 保留相关性较小的属性。首先对(6)式中的相关系数R_mm设置一个门槛值δ, 并对R_mm的上三角元素R₁₁, R₁₂, …, R_1l, R₂₂, …, R_2l, …, R_ll进行升序排列, 若相关系数R_ij＜δ则保留第i个和第j个属性; 若相关系数R_ij＞δ, 则比较(5)式中R_wm的第i个与第j个相关系数, 保留相关系数较大的属性, 剔除相关系数较小的属性。继续考察剩余的每个R_ij, 符合要求的属性全部组成新的微小断裂系统属性集。值得注意的是, 门槛值的合理设置对微小断裂系统属性集的选择至关重要, 直接影响后期微小断裂系统识别的精度。通常门槛值的设置按照属性优选的要求选择, 若精度要求高, 则门槛值设置应小些; 若精度要求低, 则门槛值设置应大些。

1.3 机器学习(支持向量机)微小断裂系统指示因子估算

支持向量机是VAPNIK^[32]在VC维(vapnik-chervonenkis dimension)和结构风险最小化准则基础上提出的, 是机器学习算法的一种, 对小样本和非线性的函数估计具有很强的适应性, 能够克服传统机器学习中局部最优的不足。其核心思想是通过核函数将低维空间的非线性数据映射到高维空间中, 实现数据的高维空间的数据线性化估计, 实际上是一个凸二次优化问题, 能够保证求取的方程是全局最优的解, 在小样本学习中具有明显的优势。

设已知的样本(微小断裂系统属性集样本和微小断裂系统指示因子样本)映射集为{(x_i, y_i), i=1, 2, …, m}, 其中x_i∈R^N为N维微小断裂系统属性集样本, 为微小断裂系统指示因子曲线的样点值, m为样本个数(已知井点数)。假设非线性映射关系为:

$ f(x)=\boldsymbol{w} \cdot \varphi(x)+b $

(7)

式中: w∈R^N为变换后的特征空间向量; b∈R为偏置量, 是常数项; φ(x)是非线性回归函数。在回归函数中引入松弛变量ξ和ξ^*, 其最优目标函数为:

$ \min \limits_{\omega, b, \xi} \frac{1}{2}\|\boldsymbol{w}\|^{2}+C \sum\limits_{i=1}^{l}\left(\xi_{i}+\xi_{i}^{*}\right) $

(8)

约束条件为:

$ \left\{\begin{array}{l} y_{i}-\boldsymbol{w} \cdot \varphi\left(x_{i}\right)-b \leqslant \varepsilon+\xi_{i} \\ -y_{i}+\boldsymbol{w} \cdot \varphi\left(x_{i}\right)+b \leqslant \varepsilon+\xi_{i}^{*} \quad i=1,2, \cdots, m \\ \xi_{i} \geqslant 0, \xi_{i}^{*} \geqslant 0 \end{array}\right. $

式中: C是惩罚系数; ε是不敏感损失函数。C越大表示超过训练误差ε范围的数据样本惩罚越大, 反映的是误差精度, ε越小表示映射关系与输出误差最小, 最终的估计结果更加精确, 其表示形式为:

$ \begin{array}{ll} |y-f(x)|_{\varepsilon}= & \\ \left\{\begin{array}{cl}0 & |y-f(x)| \leqslant \varepsilon \\ |y-f(x)|-\varepsilon & |y-f(x)|>\varepsilon\end{array}\right. \end{array} $

(9)

(9) 式是高维特征空间的凸二次规划问题, 由满足Mercer条件^[18]的内积核函数K(x, y)=〈φ(x), φ(y)〉替换线性回归函数, 采用拉格朗日乘子法求解满足该约束条件的二次规划问题^[23]:

$ \begin{gathered} \max _{\alpha_{i}, \alpha_{i}^{*}}\left\{L=\sum\limits_{i=1}^{m}\left(\alpha_{i}-\alpha_{i}^{*}\right) y_{i}-\frac{1}{2} \sum\limits_{i=1}^{m}\left(\alpha_{i}-\alpha_{i}^{*}\right) \times\right. \\ \left.\left(\alpha_{i}-\alpha_{i}^{*}\right) K\left(x_{i}, x_{j}\right)+\varepsilon \sum\limits_{i=1}^{m}\left(\alpha_{i}+\alpha_{i}^{*}\right)\right\} \end{gathered} $

(10)

(10) 式应该满足约束条件:

$ \left\{\begin{array}{l} \sum\limits_{i=1}^{l}\left(\alpha_{i}-\alpha_{i}^{*}\right)=0 \\ 0 \leqslant \alpha_{i} \leqslant C \\ 0 \leqslant \alpha_{i}^{*} \leqslant C \end{array}\right. $

其中, α_i≥0, α_i^*≥0(i=1, 2, …, m)是拉格朗日乘子。

求解(10)式, 得到最终的非线性估计函数:

$ f(X)=\sum\limits_{i=1}^{l}\left(\alpha_{i}-\alpha_{i}^{*}\right) K\left(X, X_{j}\right)+b $

(11)

其中, 拉格朗日乘子α_i与α_i^*(i=1, 2, …, m)互不相等。机器学习(支持向量机)预测微小断裂系统特征关键是确定惩罚系数C, 再根据数据特征选定有效的核函数。核函数可将数据空间映射到高维空间, 跳过了非线性映射关系, 直接在输入空间求取回归函数。

2 实际资料应用及效果

将本文方法应用于四川盆地碳酸盐岩实际工区的地震资料进行微小断裂系统的识别应用研究, 其目的层等T₀图如图 2所示, 图 3是该工区内W39井目的层储层的电镜扫描照片和岩心照片。这些资料表明研究区内的微小断裂较为发育。

以地质分析数据以及测井资料为基础, 建立碳酸盐岩储层的各向异性岩石物理模型, 考虑不同的孔隙类型、裂缝以及流体对地震参数的影响, 实现井位处储层参数和各向异性参数的估算, 进而求取断裂敏感因子曲线。具体步骤如下。

1) 利用纵波速度、密度、含水饱和度等通过模拟退火算法对泥质含量和孔隙度曲线进行校正。

2) 从矿物含量曲线出发, 利用Voigt-Reuss-Hill模型计算干燥岩石的等效矿物模量。

3) 利用DEM模型、Hudson模型和Schoernberg模型在干燥岩石中加入孔隙和裂隙得到含孔隙的干燥岩石的等效体积模量。

4) 为了实现储层的骨架和流体混合, 运用Wood方程计算孔隙-裂隙包含流体时的饱和岩石的等效体积模量; 运用Brown-Korringa模型实现该碳酸盐岩储层的各向异性流体替换, 得到Thomsen参数。

5) 利用Thomsen参数计算储层的各向异性梯度, 作为微小断裂系统指示因子曲线。

通过构建碳酸盐岩复杂断裂储层岩石物理模型, 预测W23井、W31井和W39井的Thomsen参数ε^(V), δ^(V)和γ, 然后利用(4)式计算各向异性梯度, 以此作为微小断裂系统指示因子曲线, 估算结果如图 4所示, 每口井显示的曲线从左至右分别是密度(ρ)、自然伽马(γ)曲线和微小断裂系统指示因子(Γ)曲线, 图中黄色阴影部分是目的储层段。由图 4可以发现, 储层段的裂缝敏感因子曲线均有明显的差异响应。

图 5为研究区内过井线的叠后地震剖面, 可以看出, 研究区内存在一套深部大断裂(红粗线所示)一直延伸到浅层, 另外在目的层段附近有一套小断裂(绿粗线所示), 部分微小断裂在剖面中难以识别。

得到微小断裂系统敏感指示因子后, 利用相关聚类方法优选敏感属性集, 其过程如下。

1) 利用三维地震数据体提取属性数据体并作预处理, 得到方差体、构造导向滤波体、相干、倾角、方位角、蚂蚁体、纹理属性、最大正曲率、最大负曲率、边缘检测等多种几何类地震属性。

2) 计算井位处的微小断裂系统敏感指示因子曲线与属性的相关性以及各属性之间的相关系数。

3) 根据步骤2)计算的相关系数通过设定门槛值对输入的属性集合进行相关聚类, 优选得到敏感属性集合。

通过上述步骤, 最终优选了相干、最大负曲率、最大正曲率和倾角属性作为微小断裂系统属性集, 目的层的属性集合平面图如图 6所示。

以已知的W23井、W39井和W31井井位处的微小断裂指示因子采样值和敏感属性集作为训练数据, 在目的层以上20ms和目的层以下10ms的有效层段内, 结合井震标定结果, 按照2ms的采样率选取样本, 每口井16个采样点, 共选择48个采样点处的微小断裂指示因子曲线和敏感属性集数据对。在计算过程中, 根据微小断裂指示因子与属性集合样本之间的数据特征, 选择高斯基作为支持向量机的核函数。高斯核函数的具体表示形式为:

$ \begin{gathered} K(x, y)=\exp \left(-\frac{\|x-y\|^{2}}{2 \sigma^{2}}\right)= \\ \exp \left(-g\|x-y\|^{2}\right) \end{gathered} $

(12)

式中: x是训练样本; y是核函数中心; g=1/(2σ²), 且g＞0, 称为核函数参数。

运用交叉验证的方式优选支持向量机的惩罚系数C和核函数参数g, 交叉验证的主要步骤可参考文献[33]和文献[34]。通过交叉验证方法得到误差最小时的惩罚系数C=1、核函数参数g=50, 在此基础上开展支持向量预测, 微小断裂系统指示因子预测结果如图 7所示。根据支持向量机计算过程设置测试数据评判支持向量机预测结果的准确性, 已知结果与预测结果的平均吻合率为89.6%, 说明支持向量机预测的可靠性。另外, 预测结果表明, 微小断裂敏感因子预测平面图均能正确反映微小断裂发育带, 这是相干属性和倾角属性中不易揭示的; 最大正曲率和最大负曲率虽然能够反映这些断裂, 但反映的微小断裂特征比较杂乱, 不利于识别真实的断裂信息。图 5和图 7中的蓝色箭头分别是对应的断裂特征在剖面和预测平面中的反映, 可以看出, 两者吻合程度较高。W39井的测井综合解释结果(图 8)表明, 在深度4410~4428m处为优质白云岩储层, 测试产气达到24.9×10⁴m³/d, 同时W39井的电镜扫描照片、岩心照片均显示该层段断裂比较发育, 预测结果表明W39井附近微小断裂发育, 与已知的井点信息吻合程度很高。通过上述分析, 可以认为本文机器学习(支持向量机)的预测结果能够实现复杂储层断裂特征的精细识别。同时, 结合地震剖面中的断层解释, 可有效区分深部(断层如果从深部延展到目的层就认为是深部断层)和浅部断层(断层只是目的层内分布), 进而指导微小断裂特征的准确解释(图 9), 为微小断裂系统识别提供一定的技术支持。

3 结论与认识

结合多元信息(地震、测井)的优势, 在各向异性岩石物理建模的基础上开展机器学习(支持向量机)复杂储层的微小断裂系统识别, 预测结果与已有的地质认识和测井解释结果相吻合, 形成了新的微小断裂系统预测技术体系。

1) 微小断裂指示因子曲线的计算综合考虑矿物组分、储层特征、孔隙类型以及各向异性等, 构建符合复杂储层地质特征的各向异性岩石物理模型, 获取能够准确反映井位处各向异性特征的微小断裂指示因子曲线。许多构造类地震属性能够在一定程度上反映储层的断裂特征, 但是, 需要从众多属性中优选具有代表性的属性进行机器学习, 保证储层预测结果的有效性与准确性。

2) 与常规的属性预测相比, 机器学习(支持向量机)的微小断裂系统识别方法有效地解决了已知井点少、非线性特征的储层预测难题, 可以提供高横向分辨率的微小系统断裂识别结果。实际资料的应用结果表明, 机器学习(支持向量机)预测结果与已知地质认识、测井解释匹配程度较高, 能够为复杂储层的微小断裂系统精确识别与综合预测提供可靠的地球物理依据。

3) 本文方法在微小断裂系统识别方面取得了较好的应用效果, 近几年智能物探技术得到了快速发展, 为解决复杂储层的断裂识别提供了有效途径, 这是我们下一步的重点研究内容。