纵波及横波波速均为油气储层非常重要的地球物理参数, 其中纵波波速可通过常规测井方法获得, 而横波波速的获取途径主要包括室内声学测试实验或偶极子声波测井^[1], 这两种方法均具有成本高且耗时长的不足。常规测井曲线具有测试成本较低, 且能够获得连续纵向剖面的特点, 现阶段较为常用的方法是利用常规测井曲线进行纵横波拟合, 得到连续的横波波速剖面, 进而计算储层物性及力学参数。目前, 关于横波波速的研究方法主要以概率统计法及经验公式法为主^[2-3], 这两种方法在碎屑岩储层中均有较好的回归效果。但在泥页岩储层, 尤其是砂泥互层的泥页岩储层中, 受岩性、流体、应力状态等方面的变化, 采用上述两种方法进行拟合, 其相关系数往往较低, 拟合效果不理想。横波波速是求解储层孔隙度、泊松比、杨氏模量等不可或缺的参数, 因此针对泥页岩储层的横波波速的预测工作成为泥页岩储层的重要研究内容之一。人工神经网络是现代数学的一个重要分支, 其中BP(back propagation)神经网络是目前应用最广泛的神经网络模型^[4-6]。BP网络由于具有学习和存贮大量的输入-输出映射关系而无须事前揭示描述这种的映射关系的特点^[7], 在数值预测及识别分析等方面具有广泛的应用。在石油行业, BP神经网络已应用于测井资料处理^[8-9]、岩性识别^[10-15]、孔隙度预测^[16-19]、非常规储层识别^[20-22]、裂缝识别^[23-24]等方面的研究, 但在横波波速预测方面的研究成果目前相对匮乏^[25]。笔者以声学实验结果为基础, 通过训练LM-BP神经网络, 建立考虑多条件的横波波速计算模型。

1 LM-BP神经网络的基本原理

BP神经网络是1986年由Rumelhart和McCelland等^[26-27]提出的一种按误差逆传播算法训练的多层前馈网络, 由信息的正向传播和误差的反向传播两个过程组成。其主要思想为:假设样本向量为P₁, P₂, …, P_q, 均有其对应的输出样本T₁, T₂, …, T_q。神经网络通过收集归纳样本的特征参数, 计算得到实际输出矢量A₁, A₂, …, A_q, 并求取实际输出矢量与输出样本矢量之间的误差, 依据误差修改其权值, 使A_i(i=1, 2, …, q)与输出样本T尽可能地接近。BP神经网络结构模型如图 1所示。

1.1 BP神经网络计算过程

常规BP神经网络的计算过程包括两部分:正向输出和反向误差计算。正向输出过程是依据输入层的样本参数, 通过活化函数分别计算隐含层和输出层单元的数值。反向误差计算则是根据输出层的计算结果, 反向计算输出层和隐含层的误差, 若误差大于设置的最小误差时, 则调整各层权值及阈值, 并将调整后的权值及阈值返回至输入层神经元中, 重新进行正向计算, 直至达到最大训练次数或计算误差小于设置的最小误差为止。BP神经网络反向计算过程的权值调整公式为

$ {\mathit{\boldsymbol{w}}_{i(j + 1)}} = \eta \cdot{\mathit{\boldsymbol{E}}_j}\cdot{\mathit{\boldsymbol{Y}}_k}\cdot{\mathit{\boldsymbol{Y}}_i} + \alpha \cdot\Delta{\mathit{\boldsymbol{w}}_{ij}}. $

(1)

式中, w_i(j+1)为输入层第i个节点到隐含层第j个节点的权值; Δw_ij为输入层到隐含层的误差调整值; η为学习步长; E_j为隐含层各节点的误差; Y_k为输出层第k个节点的输出值; Y_i为输入层第i个节点的输入值; α为动量因子(即动量项参数, α∈[0.5, 0.9])。

常规BP神经网络的主要计算过程是通过不断调整权值, 实现计算样本与实际样本的误差逐渐减小至设定的最小误差。然而在实际计算中, 须事先给出动量因子α及学习步长η, 目前这两个参数的赋值依据主要是通过实际经验给出, 存在较大的随机性, 计算结果往往存在较大误差。

1.2 LM算法的反向误差计算

利用LM(Levenberg-Marquardt)优化算法, 能够对BP神经网络的反向误差计算过程进行优化。LM算法是最优化算法中的一种, 它是梯度下降法与高斯-牛顿法的结合, 既有高斯-牛顿法的局部收敛性, 又具有梯度法的全局特性^[28]。设BP神经网络正向计算得到的计算值Y_k, 代入下式计算输出误差:

$ \mathit{\boldsymbol{E}}\left( \mathit{\boldsymbol{w}} \right) = \frac{1}{2}\mathop \sum \limits_{k = 1}^p {\left\| {{\mathit{\boldsymbol{Y}}_k} - \mathit{\boldsymbol{Y}}{' _k}} \right\|^2} = \frac{1}{2}\mathop \sum \limits_{k = 1}^p \mathit{\boldsymbol{e}}_k^2\left( \mathit{\boldsymbol{w}} \right). $

(2)

式中, E(w)为误差指标函数; Y_k为正向计算的神经网络输出向量; Y′_k为实际的输出样本向量; p为样本个数; w为神经网络的权值组成的向量; e_i(w)为误差。

若w_n表示第n次迭代的权值和阈值所组成的向量, w_n+1为调整后的权值和阈值所组成的向量, 则有

$ {\mathit{\boldsymbol{w}}_{n + 1}} = {\mathit{\boldsymbol{w}}_n} + \Delta\mathit{\boldsymbol{w}}. $

(3)

其中Δw为权值增量, 其表达式为

$ \Delta\mathit{\boldsymbol{w}} = {[{\mathit{\boldsymbol{J}}^{\rm{T}}}\left( \mathit{\boldsymbol{w}} \right)\mathit{\boldsymbol{J}}\left( \mathit{\boldsymbol{w}} \right) + \mu \mathit{\boldsymbol{I}}]^{ - 1}}{\mathit{\boldsymbol{J}}^{\rm{T}}}\left( \mathit{\boldsymbol{w}} \right)\mathit{\boldsymbol{e}}\left( \mathit{\boldsymbol{w}} \right){\rm{ }}. $

(4)

其中

$ \mathit{\boldsymbol{J}}\left( \mathit{\boldsymbol{w}} \right)=\left[ \begin{matrix} \frac{\partial {{\mathit{\boldsymbol{e}}}_{1}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{1}}}&\frac{\partial {{\mathit{\boldsymbol{e}}}_{1}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{2}}}&\cdots &\frac{\partial {{\mathit{\boldsymbol{e}}}_{1}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{n}}} \\ \frac{\partial {{\mathit{\boldsymbol{e}}}_{2}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{1}}}&\frac{\partial {{\mathit{\boldsymbol{e}}}_{2}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{2}}}&\cdots &\frac{\partial {{\mathit{\boldsymbol{e}}}_{2}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{n}}} \\ \vdots &\vdots &\ddots &\vdots \\ \frac{\partial {{\mathit{\boldsymbol{e}}}_{N}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{1}}}&\frac{\partial {{\mathit{\boldsymbol{e}}}_{N}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{2}}}&\cdots &\frac{\partial {{\mathit{\boldsymbol{e}}}_{N}}\left( \mathit{\boldsymbol{w}} \right)}{\partial {{\mathit{\boldsymbol{w}}}_{n}}} \\ \end{matrix} \right]. $

式中, I为单位矩阵; μ为用户定义的学习率; J(w)为Jacobian矩阵, 即权值w的雅克比矩阵。

利用LM算法训练BP神经网络, 首先须给出训练误差的允许值ε, 训练学习率μ₀和调整常数β(0 < β < 1), 并对神经网络权值及阈值向量进行初始化处理, 令首次迭代次数n=0, 学习率μ=μ₀。再根据公式(2)计算神经网络输出向量Y_i及其误差指标函数E(w), 根据公式(3)和公式(4)分别计算输出权值w的雅克比矩阵J(w)和权值增量Δw。若E(w) < ε, 则计算结束。否则, 根据公式(3)计算新的网络权值和阈值向量。若E(w_n+1) < E(w_n), 则令计算次数n=n+1, μ=μβ (若E(w_n+1)> E(w_n), 则μ=μ/β), 重复上述计算步骤, 直到达到最大训练次数或输出误差小于ε为止(图 2)。

与常规BP神经网络相比, LM-BP神经网络具有训练时长短、收敛速度快、数值稳定性较好的特点。由公式(2)~(4)可知, LM-BP神经网络的反向误差计算过程不考虑网络初始权值w₀和动量项系数α, 从而很好地解决了常规BP神经网络中存在的初始权值及动量项系数难调整的问题。同时, 由公式(4)可以看出, 当μ=0时, 公式(4)则为高斯-牛顿法; 若学习率较高, 则反向算法近似梯度下降法, 每完成一次迭代, 则μ减小一些, 这样在接近设定的最小误差的过程中, 实现由梯度下降法逐渐转变为高斯-牛顿法。另外由于[J^T(w)J(w)+μI]为正定矩阵, 因此公式(4)的解总是存在的, 保证了计算结果的收敛性, 计算效率及精度得到提高。

2 计算实例

本文中以新场气田须五段泥页岩地层为研究对象。新场气田地区位于中国四川省川西坳陷中部, 其中须五段地层埋深约为2 000~2 800 m, 地层岩性以泥岩为主, 夹杂炭质页岩、粉砂质泥岩、粉砂岩和细砂岩的等厚或不等厚互层组合^[29-30], 孔隙度平均为4.26%, 渗透率为(0.000 5~0.004 4)×10^-3 μm², 属于致密页岩气藏^[31-33]。须五段早期的研究中曾利用砂泥岩纵横波速的统计关系分别建立了砂泥岩横波时差的预测模型, 但由于考虑的条件较少, 导致横波时差拟合的相关系数较低, 线性拟合方法并不适用(图 3)。本文中从岩石机械波理论出发, 综合考虑影响岩石声波波速的因素, 分别测试了在不同岩性、饱和状态、纵波波速、围压及轴压条件下的岩石声学测试^[34-35]。根据测试结果, 建立LM-BP神经网络模型, 将上述实验条件及纵波波速作为输入样本, 以横波波速为目标输出样本, 通过训练LM-BP神经网络, 提高横波波速预测的准确率。

2.1 岩石机械波理论

若岩石为各项同性介质, 不考虑体力, 机械波在固体介质中传播的拉梅运动方程^[36]为

$ \rho \frac{{{\partial ^2}u}}{{\partial {t^2}}} = \left( {\lambda + \mu } \right)\nabla {\theta _{\rm{t}}} + \mu {\nabla ^2}\mathit{\boldsymbol{u}}{\rm{ }}. $

(5)

式中, ρ为岩石密度, kg/m³; u为位移矢量场, m; θ_t为体积应变, m³/m³; μ为剪切模量, GPa; λ为拉梅常数, GPa。其中, 拉梅常数的表达式为

$ \lambda = \frac{E}{{\left( {1 + \upsilon } \right)\left( {1 - 2\upsilon } \right)}}\upsilon {\rm{ }}. $

(6)

式中, E为动弹性模量, GPa; υ为动泊松比。

由波动方程可得到的纵横波速计算方程为

$ \left\{ \begin{array}{l} {\mathit{v}_{\rm{p}}} = \sqrt {\frac{{\lambda + 2\mu }}{\rho }}, \\ {\mathit{v}_{\rm{s}}} = \sqrt {\mu /\rho } . \end{array} \right. $

(7)

式中, v_p和v_s分别为纵波波速和横波波速, m/s。

由公式(7)可以得到各向同性介质中的纵横波波速的数学关系:

$ v_\text{p}^2 = \frac{\lambda }{\rho } + 2v_{\rm{s}}^2. $

(8)

将拉梅常数的计算公式代入式(8)中, 经化简可推出:

$ v_{\rm{s}}^2 = \frac{{E\upsilon }}{{2\rho \left( {1 + \upsilon } \right)\left( {1 - 2\upsilon } \right)}} + \frac{{v_{\rm{p}}^2}}{2}. $

(9)

根据动弹性模量E、动泊松比υ的定义公式, 公式(9)可变形为

$ v_{\rm{s}}^2 = \frac{{\sigma {\varepsilon _{\rm{a}}}}}{{2\rho ({\varepsilon _{\rm{l}}} + {\varepsilon _{\rm{a}}})({\varepsilon _{\rm{l}}} - 2{\varepsilon _{\rm{a}}})}} + \frac{{v_\text{p}^2}}{2}. $

(10)

式中, σ为施加在岩石介质的载荷, MPa; ε_l和ε_a分别为岩石介质的轴向应变和横向应变, 无量纲。

令$b = \frac{{\sigma {\varepsilon _a}}}{{2\rho ({\varepsilon _l} + {\varepsilon _a})({\varepsilon _l} - 2{\varepsilon _a})}}$, 则式(10)简化为

$ v_{\rm{s}}^2 = v_\rm{p}^2/2 + {\rm{b}}. $

(11)

从公式(10)可以看出, 岩石类型或饱和状态不同, 密度ρ不同, 岩石应力状态改变, 其应变量也会随之变化, 这些对b值均有影响, 导致纵横波波速关系发生改变。若仅利用纵横波平方关系线性统计拟合, 可能会得到近似线性关系, 但往往相关系数较低, 很难保证拟合精度。鉴于此, 针对影响纵横波关系的多项参数, 设计了在不同岩性、饱和状态、围压及轴压条件下的岩石声学实验, 得到了各条件下纵横波波速特征。

2.2 陆相泥页岩声学响应特征

岩石的声学实验采用成都理工大学油气藏地质及开发国家重点实验室多功能岩石物理测试系统, 根据理论推导结果, 采用等侧向压缩三轴实验分别测试了不同岩性在不同条件下的岩石纵波及横波波速, 测试样本点共计363个。根据测试结果, 样品岩性及测试条件改变, 纵横波的传播速度均有不同程度的变化。

2.2.1 岩性

岩性变化对纵横波的传播速度的影响较为明显, 在相同围压及相同饱和状态的条件下, 砂岩波速整体大于泥岩波速(图 4)。

2.2.2 饱和状态

饱水后样品的纵波值明显升高, 砂岩升高幅度为7.93%, 泥岩升高幅度为2.11%(图 5), 而饱水后的样品横波波速值与饱气条件下相比, 大部分样品略降, 砂岩整体降低幅度为3.31%, 泥岩整体降低幅度为6.73%。

2.2.3 围压

围压的改变对波速的影响整体表现为正相关关系, 即随着围压的增加, 纵横波波速值均增加(图 6)。但增加到较高围压条件下, 单位围压下纵波波速的增幅已经不太明显, 并开始出现下降的趋势; 而横波波速的变化趋势明显, 随着围压的增加, 波速阶段增量及所占的比例逐渐下降。

2.2.4 轴压

轴压对纵横波波速的影响(图 7)主要体现在随着轴压的增加, 起初会有一个初始压密段, 表现为起初纵、横波波速, 特别是纵波波速会出现明显增加。当轴压增大到较大值时, 逐渐产生微裂纹, 导致纵横波波速稳定或降低。由于轴压的增加是导致环向应变的主导因素, 因此对于横波来说, 变化不大或出现降低, 而纵波则会出现一定程度升高。

2.3 LM-BP神经网络计算过程

依据声学实验过程中的各项条件参数, 选取岩性、纵波波速、轴压、围压、饱和状态为输入层样本, 横波波速为输出层样本, 样本数总计363个。岩性以砂泥岩区分, 设砂岩值为1, 泥岩为0;饱和状态设饱水状态为1, 饱气为0。为消除样本数据在量纲及数量级上的差异, 首先对样本数据进行标准化处理, 标准化公式为

$ x{' _{ij}} = \frac{{{x_{ij}} - {{\overline x }_j}}}{{{S_j}}}. $

(12)

式中, x′_ij为样本标准化参数; x_ij为样本参数(i=1, 2, …p, j=1, 2, …n); x_j为样本第j列平均值; S_j为样本第j列标准差。

计算模型采用3层BP神经网络模型, 活化函数为S函数(Sigmoid函数), 训练学习率为0.1, 隐含层神经元阈值设置为(-0.5, 0.5)区间内的随机值, 误差反馈算法为LM算法, 训练次数为10 000次, 学习总误差为0.000 1, 输入层神经元为5个, 隐含层节点数10个, 输出层节点为1个, 建立横波波速预测的LM-BP神经网络结构。首先正向计算横波波速, 再利用公式(2)~(4)得到横波正向计算结果与样本参数的误差指标函数E(w)及权值增量Δw, 并依据计算误差, 调整隐含层和输出层的权值w_i及阈值T_i, 再重复正向计算结果, 直到误差指标函数E(w)小于设定的学习总误差0.000 1为止。表 1为神经网络样本标准化后的计算参数。

通过训练BP神经网络, 得到横波波速的计算值, 统计计算结果表明, 平均相对误差为2.22%, 相对误差的最大值为10.79%, 相对误差在5%以内的样本数为336个, 占样本总数的92.56%(图 8)。与纵横波的统计回归方法相比, LM-BP神经网络由于在计算过程中已对岩性进行了区分, 因此在利用测井资料计算横波时差时, 无须再对计算结果按岩性进行分类即可得到连续的横波时差计算结果, 降低了计算横波时差剖面及岩石动力学参数的工作量。同时, 该方法综合了岩性、饱和状态、轴压及围压等对纵横波波速的影响, 横波波速计算结果的准确率得到提高(图 9)。

本文中利用LM-BP神经网络计算横波波速, 共需要5个实验参数:有效围压、轴压、纵波波速、饱和状态和岩性参数。在实际计算中, 岩性参数可根据不同岩性分别以数字代替; 纵波波速则来自于常规测井曲线声波时差测井数据计算得到; 饱和状态可以通过测井解释及生产测试资料获取; 轴压及有效围压可以根据井点地应力测试资料计算得到, 其中轴压为第一主应力σ₁, 而有效围压可以通过其定义公式求取。有效围压可通过下式估算得到:

$ {p_{\rm{e}}} = p - {p_{\rm{f}}} $

式中, p_e为有效围压, MPa; p为小主应力σ₃测试结果, MPa; p_f为地层流体压力, MPa。

3 结论

(1) 利用LM算法能够有效回避网络初始权值、动量项参数等取值方面的困难, 并通过将高斯-牛顿算法和梯度下降法的结合, 达到提高BP神经网络的运算速度和精度的效果。

(2) 泥页岩地层纵横波波速受多个条件影响, 若忽略这些条件, 拟合结果精度较低。利用机械波理论对纵横波波速进行分析, 并通过多条件岩石声学实验确定影响纵横波波速的条件主要包括岩性、饱和状态及应力状态。

(3) 将影响波速的多个条件及纵波波速作为输入样本, 横波波速作为预期输出样本。利用LM-BP神经网络预测横波时差, 相比经验公式法及统计回归法, 能够更全面地考虑影响纵横波波速的条件, 横波时差的预测精度有所提高。