近几年来, 随着能源需求不断提升和常规油气资源逐步枯竭, 以页岩油气为重点的非常规资源已成为焦点, 通过储层评价有效定位页岩油气的“甜点区”变得尤为重要^[1]。总有机碳含量(total organic carbon, TOC)是用来反映页岩油气富集程度及生烃潜力的重要参数^[2]。目前利用ΔlogR法^[3]进行TOC预测时, 精度偏低且需要人工选定基线, 推广受到一定限制。人工神经网络(artificial neural network, ANN)具有较强的非线性逼近能力, 解决非均质地层非线性映射问题具有较好的应用效果, 正逐渐应用到TOC预测中^[4-8], 但受同步瞬时输入限制, 逼近精度较低。为解决此问题, 何新贵^[9]提出一种连续过程神经元网络(continuous process neural network, CPNN), 已用于故障诊断^[10]、非线性系统建模^[11]等实际问题中。但CPNN存在以下问题:①离散数据拟合时可能得不到合适的解析函数; ②解析函数拟合阶次过高, 易导致龙格震荡现象。综上所述, 受CPNN启发, 为提高TOC的拟合预测精度, 笔者提出一种极限学习离散过程神经网络(ELM discrete process neural network, ELM-DPNN), 输入是与TOC相关性较强的测井曲线离散数据序列, 用向量模拟每一个过程输入, 内部为体现测井曲线随深度变化的时间累积效应, 采用抛物插值数值积分处理离散输入。

1 极限学习机

极限学习机(extreme learning machine, ELM)是由Huang等^[12]提出的一种神经网络快速学习算法, 文献[13]通过理论分析和实验已证明可极大提高学习速度(通常10倍以上), 目前ELM理论体系虽不尽完善, 但处于迅速发展中, 已成为国内外的研究热点^[14-15]。

对于单隐层前馈神经网络和N个任意样本, X=(x₁, x₂, …, x_N)为输入, T =(t₁, t₂, …, t_N)为输出, 则模型输入输出映射关系为

$ \sum\limits_{i = 1}^L {{\beta _i}g\left( {{\mathit{\boldsymbol{a}}_i} \cdot {\mathit{\boldsymbol{x}}_j} + {b_i}} \right)} = {o_j},j = 1, \cdots ,N. $

(1)

式中, g为激励函数; a_i为隐层输入权值; β_i为隐层输出权重。模型学习目标是使误差函数E(w)最小, E(w)为目标输出和期望输出的误差平方和,

$ E\left( w \right) = \sum\limits_{j = 1}^N {{{\left\| {{o_j} - {t_j}} \right\|}^2}} = 0, $

(2)

即存在β_i, a_i和θ_i使得

$ \sum\limits_{i = 1}^L {{\beta _i}g\left( {{\mathit{\boldsymbol{a}}_i} \cdot {\mathit{\boldsymbol{x}}_j} + {b_i}} \right)} = {t_j},j = 1, \cdots ,n. $

(3)

记H为隐层的输出矩阵, β为隐层输出权重向量, 则式(3)可表示为

$ \mathit{\boldsymbol{H\beta }} = \mathit{\boldsymbol{T}}. $

(4)

使用传统的梯度下降法训练神经网络, 不仅计算量大、训练速度慢, 而且对初值选择敏感, 易陷入局部最小。为此Huang提出以下定理:

引理1^[12]  对于N个不同样本, M个隐层节点和一个任意区间无限可导的激活函数g, 则神经网络在a_i∈Rⁿ和b_i∈Rⁿ任意赋值的情况下, 隐层输出矩阵H可逆, 即方程组有精确解, 代价函数E(W)=0。

引理2^[12]  给定任意N个不同样本, 任意小误差e > 0, 及在任意区间无限可导的激活函数g, 总存在一个神经网络, 使a_i∈Rⁿ和b_i∈Rⁿ任意赋值情况下, 误差E(W)≤e。

上述定理表明神经网络只要在权值和阈值赋值的情况下, 可完成对系统的逼近, 从而将神经网络的训练过程转化为式(4)的最小二乘解问题, 其中隐层输出权值为

$ \mathit{\boldsymbol{\beta }} = {\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{H}}} \right)^ - }{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{T}}. $

(5)

2 ELM-DPNN模型

设x₁(t_l), x₂(t_l), …, x_n(t_l)是离散过程神经元的离散向量输入, w₁(t_l), w₂(t_l), …, w_n(t_l)为隐层输入权值的离散序列, 则ELM-DPNN结构如图 1所示。

模型映射关系可以表示为

$ y = \sum\limits_{j = 1}^m {{v_j}f\left( {\sum\limits_{i = 1}^n {\int_0^{\rm{T}} {{w_{ij}}\left( t \right){x_i}\left( t \right)dt - \theta } } } \right)} . $

(6)

为完成离散样本输入的时域聚合, 简化连续积分计算, 本文中提出使用抛物插值的数值积分处理离散样本输入。设积分区间[a, b], x_k=a+kh, y_k=f(x_k), $ h = \frac{{b - a}}{n} $, 则f(x)在[a, b]上抛物插值可以表示为

$ \begin{array}{l} I\left( {f\left( x \right)} \right) = \frac{{{h_0}}}{6}\left( {3{f_0} + 3{f_1} - {R_0}} \right) + \\ \frac{1}{6}\sum\limits_{i = 1}^{n - 2} {{h_i}\left( {3{f_i} + 3{f_{i + 1}} - \frac{{{L_i} + {R_i}}}{2}} \right)} + \frac{{{h_{n - 1}}}}{6}\left( {3{f_{n - 1}} + 3{f_n} - {L_{n - 1}}} \right). \end{array} $

(7)

其中

$ {L_i} = \frac{{\delta _i^2}}{{1 + {\delta _i}}}{f_{i - 1}} - {\delta _i}{f_i} + \frac{{{\delta _i}}}{{1 + {\delta _i}}}{f_{i + 1}},i = 1,2, \cdots ,n - 1; $

$ {R_i} = \frac{{\lambda _i^2}}{{1 + {\lambda _i}}}{f_i} - {\lambda _i}{f_{i + 1}} + \frac{{{\lambda _i}}}{{1 + {\lambda _i}}}{f_{i + 2}},i = 0,1, \cdots ,n - 2; $

$ {\lambda _i} = \frac{{{h_i}}}{{{h_{i + 1}}}},h = {x_{i + 1}} - {x_i},{\delta _i} = \frac{{{h_i}}}{{{h_{i - 1}}}},i = 0,1, \cdots ,n - 1. $

根据抛物插值数值积分, 式(6)简化为

$ y = \sum\limits_{j = 1}^m {{v_j}f\left( {\sum\limits_{i = 1}^n {I\left( {{w_{ij}}\left( t \right){x_i}\left( t \right)\Delta t} \right) - {\theta _j}} } \right)} . $

(8)

3 ELM-DPNN学习算法 3.1 ELM-DPNN结构分析

对采用抛物插值数值积分完成离散样本的时域聚合的离散过程神经元进行结构展开, 如图 2所示。

从展开结构中可以看出, 离散过程神经元的输入不再是普通ANN的同步瞬时输入(一维向量), 而是样本采样时间序列x_i(t_l)组成的二维矩阵。同时, 为简化时域聚合运算, 对x_i(t)和w_ij(t)的内积进行数值积分, 体现样本对时间的累积效应。但经过数值积分简化后, 展开的隐层离散过程神经元的结构和普通神经网络基本相同, 是一个规模较为庞大的普通ANN, 若采用传统的梯度下降法, 训练速度慢, 且不易收敛。根据Huang等^[12]提出的定理1和定理2, 在输入权值和隐层阈值随机赋值的情况下, 存在一组输出权值使得网络模型逼近非线性系统, 即可以使用ELM进行离散过程神经网络的训练。

3.2 ELM-DPNN学习算法

(1) 根据试探法, 确定隐层离散过程神经元的数量, 并对阈值和输入权值向量进行随机赋值。

(2) 计算隐层输出矩阵。设(x₁^(k) (t_l), x₂^(k)(t_l), …, x_n^(k) (t_l), d_k)是给定的K个样本, k=1, 2, …, K, l=1, 2, …, L, 为表述方便, 记

$ f\left( {{w_{ij}},x_i^{\left( \mu \right)},{\theta _k}} \right) = f\left( {\sum\limits_{i = 1}^n {\sum\limits_{l = 1}^T {{w_{ij}}\left( t \right)x_i^{\left( \mu \right)}\left( {{t_l}} \right)\Delta {t_l}} + {\theta _k}} } \right). $

(9)

则ELM-DPNN的隐层输出矩阵H如下:

$ \mathit{\boldsymbol{H}} = {\left[ {\begin{array}{*{20}{c}} {f\left( {{w_{i1}},x_i^{\left( 1 \right)},{\theta _1}} \right)}& \cdots &{f\left( {{w_{im}},x_i^{\left( 1 \right)},{\theta _m}} \right)}\\ \vdots&\vdots&\vdots \\ {f\left( {{w_{i1}},x_i^{\left( K \right)},{\theta _1}} \right)}& \cdots &{f\left( {{w_{im}},x_i^{\left( K \right)},{\theta _m}} \right)} \end{array}} \right]_{K \times m}}. $

(10)

H中第k行是第k个样本x_i^(k) (t)在每个隐层离散过程神经元的输出向量; 第j列是对于输入样本中的每个样本在第j个隐层离散过程神经元的输出向量。根据式(4), ELM-DPNN输入输出映射关系可写为

$ \left[ {\begin{array}{*{20}{c}} {f\left( {{w_{ij}},x_i^L,{\theta _k}} \right)}& \cdots &{f\left( {{w_{ij}},x_i^L,{\theta _k}} \right)}\\ \vdots&\vdots&\vdots \\ {f\left( {{w_{ij}},x_i^L,{\theta _k}} \right)}& \cdots &{f\left( {{w_{ij}},x_i^L,{\theta _k}} \right)} \end{array}} \right] \times \left[ {\begin{array}{*{20}{c}} {{v_1}}\\ {{v_2}}\\ \vdots \\ {{v_m}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{d_1}}\\ {{d_2}}\\ \vdots \\ {{d_k}} \end{array}} \right]. $

(11)

为方便描述，记O=[d₁, d₂, …, d_k]_k×1^T和V=[v₁, v₂, …, v_m]_m×1^T, O和V分别为期望输出向量和隐层输出权值向量。

(3) 计算隐层输出权值。根据ELM-DPNN的隐层输出权值向量如下:

$ \mathit{\boldsymbol{V}} = {\left( {{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{H}}} \right)^ - }{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{O}}. $

(12)

记A=H^TH, 在求解(H^TH)^-时, 若A为非奇异方阵, 则可以根据|A|直接求解A^-。但A如果是奇异方阵, 无法根据|A|求解A^-。本文中提出使用SVD求解A的Moore-Penrose广义逆矩阵A^†, 适用于A为非奇异方阵和奇异方阵的两种情况。按照SVD分解法, 将A分解如下:

$ \mathit{\boldsymbol{A}} = \mathit{\boldsymbol{P}}{\bf{diag}}\left( {{\sigma _1},{\sigma _2}, \cdots ,{\sigma _r},0, \cdots ,0} \right){\mathit{\boldsymbol{Q}}^{\rm{T}}}. $

(13)

则Moore-Penrose广义逆矩阵A^†为:

$ {\mathit{\boldsymbol{A}}^\mathit{\dagger }} = \mathit{\boldsymbol{Q}}{\bf{diag}}\left( {{\sigma _1},{\sigma _2}, \cdots ,{\sigma _r},0, \cdots ,0} \right){\mathit{\boldsymbol{P}}^{\rm{T}}}. $

(14)

式中, P为m×m的正交矩阵; Q为n×n的正交矩阵, 矩阵A的秩为r。记矩阵A的特征值为λ₁≥λ₂≥…≥λ_r>λ_r+1=…=λ_n=0, 则奇异值$ {\sigma _i} = \sqrt {{\lambda _i}} $(i=1, 2, …, r)。

记Σ=diag(σ₁, σ₂, …, σ_r), σ₁≥σ₂≥…≥σ_r > 0, 设Q₁和Q₂分别为Q的前r列和后n-r列, 即Q =(Q₁, Q₂), 根据下式求解Q₁和Q₂:

$ \mathit{\boldsymbol{Q}}_1^{\rm{T}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{Q}}_1}\mathit{\boldsymbol{ = }}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^2},\;\;{\left( {\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{Q}}_2}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{Q}}_2}} \right) = \mathit{\boldsymbol{O}}. $

(15)

记P=(p₁, p₂, …, p_m), 其中列向量p_i为A^TA的正交单位特征向量, 且满足$ {\mathit{\boldsymbol{p}}_i} = \frac{{\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_i}}}{{{\sigma _i}}} $。最后利用Q和P求解A^†, 按式(12)计算隐层输出权值向量。

4 TOC预测实例

本文中以A地区B1井和B2井为例, 其中B1井在深度2 930~3 098 m的泥页岩TOC值为2%~15%, 属于TOC中高值情况, B2井在深度3 055~3 140 m的TOC值小于3%, 属于TOC低值情况, 根据这两种情况分别验证所提方法有效性。常规测井资料包括自然电位SP、自然伽马GR、深侧向电阻率RLLD、密度DEN、补偿中子CNL、声波时差AC等。

实验对比包括两部分:①模型学习能力对比。与文献[5]的RBF-ANN、文献[16]的正交基展开连续PNN(orthogonal basis-based PNN, OBPNN)、BP-DPNN进行对比, 检验本文的ELM-DPNN的学习能力是否提高; ②TOC的拟合效果和预测能力对比。将ELM-DPNN与ΔlogR法和文献[6]的RBF-ANN进行对比, 检验TOC的拟合预测精度是否提高。

4.1 数据样本选取

利用神经网络进行TOC预测, 测井曲线的选择很重要。以B1井为例, 图 3为各条测井曲线和TOC的相关性分析, 可以看出DEN、CNL和AC与TOC的相关性要明显优于SP、GR和RLLD, 因此为更好地进行TOC和测井曲线的非线性映射, 选择DEN、CNL、AC作为模型的输入。同理, 经过相关性分析, B2井也选取DEN、CNL、AC三条测井曲线。

在样本选择方面, B1井共选取180组代表性样本, 由DEN、CNL、AC测井曲线和TOC实测值构成, 其中135组作为训练样本, 其余45组为测试样本, 对神经网络模型和ΔlogR法预测的TOC进行效果验证。同理, B2共选取46组样本, 38组作为训练样本, 8组为测试样本用来验证。

由于每条测井曲线量纲不同, 数据量级相差较大, 为了让神经网络模型能够有效反映每条测井曲线与TOC的非线性映射关系, 在输入模型前必须对每条测井曲线的离散数据进行归一化处理, 如下:

$ {X^ * } = \frac{{X - {X_{\min }}}}{{{X_{\max }} - {X_{\min }}}}. $

式中, X为原始的测井曲线数值; X_min和X_max为该条测井曲线离散数据样本的最小值和最大值; X^*为归一化后的测井曲线数值, 即X^*∈[0, 1]。

4.2 TOC预测结果对比 4.2.1 模型训练情况对比

实验环境说明:所有训练实验都在Intel i7-3520、2.90 GHz主频、8G内存、win 7-64位、Matlab2012上完成。最大迭代步数maxGen=5 000, 误差精度ε=0.000 5, 模型的拓扑结构为3-11-1。其中OBPNN的正交基展开函数选用勒让德基函数, 正交基函数项数L=21, OBPNN和BP-DPNN的学习速度α=β=γ=0.5, 惯性系数η=0.1。

针对B1和B2井的训练样本, 各独立进行10次训练实验, 对结果数据进行统计, 对比情况如表 1所示。结果表明, 在训练速度方面, ELM-DPNN的训练速度较BP-DPNN和OBPNN有了很大的提升, 约5~8倍, 有效解决了过程神经网络学习训练速度较慢的问题。

4.2.2 B1井TOC预测情况

B1井埋深2 930~3 098 m的泥页岩TOC值属于中高值情况。将训练好的RBF-ANN、ELM-DPNN进行45组TOC测试样本的拟合预测, 与ΔlogR法进行对比, 结果如表 2所示, 其中ELM-DPNN的训练样本和测试样本拟合预测的TOC相关性交互图如图 4所示。

通过均方根误差RMSE、相关系数R等预测指标可发现, 对于B1井TOC拟合预测精度从高到低依次为ELM-DPNN、RBF-ANN和ΔlogR法。对于ELM-DPNN模型, 135组训练样本的实测TOC和预测TOC的相关系数可达到R=0.901 0(判定系数R²=0.811 9), 45组测试样本的实测TOC和预测TOC的相关系数R=0.862 3(判定系数R²=0.743 6), 对于B1井2 930~3 098 m井段的TOC预测, ELM-DPNN的预测精度比ΔlogR有明显提高, 相关系数提高约13%, 比文献[5]的RBF高出约6%, 更好地体现了TOC的整体变化趋势。图 5为三种方法在B1井的TOC预测情况对比。

4.2.3 B2井TOC预测情况

B2井埋深3 055~3 140 m的泥页岩TOC值属于TOC低值的情况。本次实验的ΔlogR、RBF-ANN、ELM-DPNN三种方法预测TOC的相关系数分别为0.856 3、0.885 6、0.888 9, 均较好地完成TOC的预测, ELM-DPNN比ΔlogR略有提高, 约为3%, 其他TOC拟合预测指标对比见表 2, 三种方法在B2井的TOC预测情况对比如图 6所示。

4.3 实验结果分析 4.3.1 方法适用性分析

首先, ΔlogR法的前提是假定地层岩性基本不变, 声波和孔隙度的差异变化仅由有机质含量变化引起, 使用声波和电阻率的幅度差计算TOC, 因此该方法本质上是假定TOC与电阻率、声波之间存在线性映射关系。但在A地区, 储层段岩性发生改变, 电阻和声波曲线间的差异不能单纯反映地层TOC含量的变化。当TOC较低时, 对声波、电阻率的变化影响较小, 使用ΔlogR法计算TOC存在一定的偏差。当TOC较高时, 通过对B1井进行岩心TOC相关性分析, 电阻率与TOC间的相关性较差, ΔlogR法的测井响应线性假定与实际不符, 导致TOC预测精度偏低。

其次, TOC与测井曲线之间应属于一种非线性映射关系, TOC的变化在测井曲线上均有一定的响应。由于神经网络具有很强的映射能力, 选择与TOC相关性较高的多条测井曲线作为神经网络的特征参数, 考虑了多种测井方法对TOC的综合响应, 避免了由于复杂岩性变化导致的TOC预测模型与测井响应不匹配的问题, 利用神经网络的特性可完成复杂非线性系统的映射逼近, 因此TOC的拟合预测精度要优于ΔlogR法。

最后, 本文中所提的ELM-DPNN模型通过测井曲线的离散数据向量模拟过程式输入, 以矩阵作为模型整体输入, 隐层利用抛物插值数值积分的时间聚合算子, 减小了测井曲线过程信号的时间累积误差, 因在非线性系统的逼近能力上优于普通的神经网络, TOC的拟合预测精度较高。

综合所述理论分析和实验结果数据可知, 当电阻率与TOC相关性较差时, 对于低值、中高值的TOC拟合预测问题, ELM-DPNN具有较好的适用性。

4.3.2 学习速度分析

ELM-DPNN在学习速度上明显优于BP-DPNN和OBPNN, 是因为ELM-DPNN对输入权值、阈值进行随机赋值, 通过使用Moore-Penrose广义逆计算出隐层的输出权值, 无需计算梯度信息, 学习速度较BP-DPNN提升约5倍。文献[11]提出的OBPNN在训练过程中须对离散输入进行拟合和正交基展开, 计算量大, 训练时间长, 训练时间约为ELM-DPNN的8倍。同时由于OBPNN在拟合过程中可能存在解析函数不存在、项数不易确定等问题, 导致模型的训练收敛能力劣于ELM-DPNN, 仅达到约60%。这些问题在ELM-DPNN中有明显改善。

5 结束语

为提高TOC拟合预测精度, 减小测井曲线随深度变化的时间累积误差, 提出ELM-DPNN预测模型, 并根据极限学习理论完成模型参数求解, 最后以B1和B2井为例进行TOC预测。从实验结果可以看出, 电阻率与TOC的线性关系不明显时, 模型预测效果优于ΔlogR法, 从而为TOC预测提供了一种新的智能方法。但ELM-DPNN也存在一定不足, 为完成模型训练, 所需实测的TOC样本较多, 否则易导致模型训练欠拟合情况发生, 影响预测精度, 但在理论上训练后的模型可以用于该地区其他井的TOC预测。