神经网络作为一种由数据驱动的非线性算法, 因其“端对端”的特点以及强大的学习能力, 广泛应用于图像处理、语音识别、自然语言处理等领域^[1]。根据神经网络所承担任务的数学特点可分为:聚类任务、分类任务、回归任务、重构任务等。

在地球物理领域, 聚类任务主要包括地震相划分^[2], 常使用自组织映射(SOM)神经网络; 分类任务包括储层识别^[3]、初至波识别^[4]、断层识别^[5]等, 常使用概率神经网络(PNN)和具有二维空间适应能力的卷积神经网络(CNN)等; 回归任务包括岩石物理参数预测^[6]、储层参数预测^[7]、测井曲线预测^[8]等, 往往使用前馈神经网络(BP)、径向基神经网络(RBF)、循环神经网络(RNN)等; 重构任务包括叠前、叠后地震数据插值重构^[9]、地震数据去噪^[10]等, 主要使用卷积神经网络(CNN)、深度残差网络(RESnet)以及U型神经网络(U-net)等。

在回归任务中的储层参数预测研究方面, 除神经网络方法外, 树模型类方法^[11]、支持向量机类方法^[12]、多元回归类方法^[13]等方法也被广泛应用。在传统方法基础上, 近年提出了以深度神经决策森林(DNDF)^[14]、自适应神经决策树(ANT)^[15]为代表的多方法融合深度学习模型。虽然这些模型多适用于分类任务, 但只须将模型进行少许变形即可适用于回归任务, 如由深度神经决策森林(DNDF)发展而来的深度回归森林(DRF)^[16]。本质而言, 这些方法或模型均可视为一个复杂的、非线性的函数, 即每个输入都有确定唯一的输出与之对应, 因此这些方法或模型仅适用于复杂唯一性问题。

地震数据反演通常具有非唯一解, 并且受噪声等因素影响较大。由地震数据估计储层参数具有明显的不确定性, 因此直接应用上述方法或模型估计储层参数会存在一定程度的误差。利用概率分布解决不确定性问题效果较好。结合神经网络和概率分布两者优势的混合密度网络是目前解决不确定性问题的较好选择。

混合密度网络(MDN)最早由BISHOP等^[17]和CORNFORD等^[18]提出, 该网络不仅拥有神经网络的学习能力, 还拥有概率分布描述不确定性问题能力, 已经成功应用于具有不确定性的海面风场反演。陈坤堂等^[19]和冯倩等^[20]联合BP神经网络和混合密度网络成功完成海面风场反演, 证明了该网络的有效性。SHAHRAEENI等^[21]将该网络引入岩石物理反演。殷习容等^[22]通过该网络, 利用横纵波速度反演含水饱和度、泥质含量和孔隙度等岩石物理参数, 同时测试该网络抗噪声能力, 证明该网络对岩石物理参数反演的有效性。但混合密度网络在储层参数预测方面的研究文献较少。

由于地震数据与储层参数的不确定性关系, 通过混合密度网络得到的储层参数“观测”概率分布往往呈现“多峰”, 为提升储层参数预测的精确度, 需要对该概率分布进行约束。由于地下介质具有连续分布的、空间横向渐变性的特征, 因而可以假设相邻地震道储层参数之间具有一定的相似性^[23]。以该假设为基础, 可以从已知井点的准确值出发, 建立相邻道的“估计”概率分布来约束“观测”概率分布。相邻道“估计”概率分布的约束, 使得“观测”概率分布的“多峰”性得到一定程度压制。以往研究表明受“过拟合”等影响, 单独使用神经网络得到的结果往往分布杂乱, 空间横向连续性差^[24]。因此反演中需考虑相邻道储层参数, 使得相邻道之间储层参数过渡更为平滑。

本文在混合密度网络的基础上, 以相邻道横向渐变为约束条件, 提出了基于空间近似概率约束的混合密度网络储层参数预测方法。以砂体厚度预测为例, 分别利用本文方法对理论模型和实际资料进行了试验, 同时与BP神经网络和混合密度网络的预测结果进行了对比, 并结合预测结果分析了3种方法在地震储层参数预测方面的效果。

1 算法原理 1.1 混合密度网络 1.1.1 混合模型

将多个不同频率的正弦或余弦信号进行组合可以得到一个复杂信号。同理, 将多个简单概率分布进行组合也可以得到一个复杂概率分布, 即混合概率分布模型。假设简单概率分布为高斯概率分布(N), 在地震数据反演中地震属性(s)已知, 则储层参数(r)的概率分布可表述如下:

$ P(r|s) = \sum\limits_{g = 1}^G {{w_g}} (s)N(r|{\mu _g}(s),\sigma _g^2(s)) $

(1)

式中:P(r|s)为复杂概率分布; w_g(s), μ_g(s)和σ_g²(s)分别为第g个高斯概率分布的权值、均值和方差。

1.1.2 混合密度网络结构

混合密度网络数学模型来自于混合概率分布模型。与BP神经网络相似, 混合密度网络一般由4部分组成:输入层, 隐藏层, 预处理层和输出层(图 1), 但其输出层输出的值并不是储层参数, 而是一系列与混合概率分布模型相关的参数^[25]。

对于上述混合密度网络结构, 其前向过程可分为:输入层-隐藏层、隐藏层-预处理层以及预处理层-输出层3个部分。

输入层-隐藏层与BP神经网络一致, 采取全连接方式计算:

$ {h_j} = f(\sum\limits_{i = 1}^I {w_{ji}^h} {x_i} + b_j^h) $

(2)

式中:h_j和x_i分别为隐藏层第j个节点值和输入层第i个节点值; f为激活函数; I为输入层节点的数目; w_ji^h为全连接权值; b_j^h为偏置。

隐藏层-预处理层也采取全连接方式计算, 但不使用激活函数:

$ z_k^c = \sum\limits_{j = 1}^J {w_{kj}^c} {h_j} + b_k^c $

(3)

式中:z_k^c为预处理层c类第k个节点值; c包括权值w、方差σ和均值μ三类; J为隐藏层节点的数目。

根据高斯混合复杂概率分布的限制条件, 对预处理层-输出层不同类别节点采取不同处理方式:各高斯分布权值为正数, 权值之和为1以及方差为正数。则由混合密度网络得到的“观测”概率分布(P^obs)可表述如下:

$ \left\{ {\begin{array}{*{20}{l}} {{P^{{\rm{obs}}}}(y|x) = \sum\limits_{k = 1}^K {{w_k}} N(y|{\mu _k},\sigma _k^2)}\\ {{w_k} = \frac{{{\rm{exp}}(z_k^w)}}{{\sum\limits_{k = 1}^K {{\rm{exp}}} (z_k^w)}}}\\ {\sigma _k^2 = {\rm{exp}}(z_k^\sigma )}\\ {{\mu _k} = z_k^\mu } \end{array}} \right. $

(4)

根据最大似然估计, 混合密度网络的误差函数L一般为:

$ L = - \sum\limits_{d = 1}^D {{\rm{ln}}} {P^{{\rm{ obs }}}}({Y_d}|{X_d}) $

(5)

式中:X_d和Y_d分别为训练集中第d个训练样本的输入和输出; D为参与训练样本数目。

根据误差函数, 使用BP算法进行训练。混合密度网络可视为BP神经网络的一种概率上的推广^[25], 在单独使用混合密度网络时, 往往将预测区间内概率最大处作为预测输出。

1.2 相邻道估计模型

以往神经网络仅仅单独应用于每一地震道, 未考虑地层介质在横向上具有的渐变性, 由此可假设相邻地震道储层参数之间具有一定相似性^[23], 在此基础上建立相邻道“估计”数学模型。对于二维地震数据, 每一地震道均具有唯一的一维坐标n, 其中n为该地震道的inline号。若第n道储层参数r_n未知, 而第n-1道储层参数r_n-1已知, 则r_n的“估计”概率分布(P^est)可表述如下:

$ {P^{{\rm{est}}}}({r_n}|{r_{n - 1}},D_{n,n - 1}^2) = N({r_n}|{r_{n - 1}},D_{n,n - 1}^2) $

(6)

式中:D_{n, n-1}²为方差, 其值大小反映了横向变化的剧烈程度, 可通过相邻道之间的地震属性差异进行定量描述。该差异有多种选择, 其中马氏距离计算过程中考虑整个样本集, 更能反映两道之间的差异, 所以通过马氏距离表示D_{n, n-1}²:

$ D_{n,n - 1}^2 = ({\mathit{\boldsymbol{s}}_n} - {\mathit{\boldsymbol{s}}_{n - 1}}){\mathit{\boldsymbol{C}}^{ - 1}}{({\mathit{\boldsymbol{s}}_n} - {\mathit{\boldsymbol{s}}_{n - 1}})^{\rm{T}}} + {\rm{const}} $

(7)

式中:s_n和s_n-1为第n道和第n-1道的地震属性; C为整个工区地震属性组成矩阵的协方差矩阵; const为常数, 表示整个工区的基础差异。

当第n-1道储层参数r_n-1不确定, 仅知其概率分布为bel(r_n-1)时, 则r_n的“估计”概率分布(P^est)可表述如下:

$ \begin{array}{*{20}{c}} {{P^{{\rm{ est }}}}({r_n}|{r_{n - 1}},D_{n,n - 1}^2) = \int {{\rm{bel}}({r_{n - 1}})} \cdot }\\ {N({r_n}|{r_{n - 1}},D_{n,n - 1}^2){\rm{d}}{r_{n - 1}}} \end{array} $

(8)

实际地震数据多为三维数据, 对此每一地震道均具有唯一的二维坐标(n, m), 其中n和m分别代表其inline号和crossline号。仿照卡尔曼滤波相邻估计方法, 将其相邻估计方法应用于二维空间(图 2), 其中蓝色表示第一轮相邻估计, 红色表示第二轮相邻估计。卡尔曼滤波相邻估计多用于一维时间域, 是单向、单点且不交叉的(图 2a), 而二维空间域相邻估计是多向、集体且存在交叉的(图 2b)。

1.3 储层参数预测

对于储层参数未知的(n, m)道, 使用混合密度网络可以得到该道储层参数r_{n, m}的“观测”概率分布(P^obs), 同时根据相邻道“估计”关系((6)式和(8)式)得到(n, m)道储层参数r_{n, m}的“估计”概率分布(P^est)。由于二维空间域相邻估计具有交叉的特点(图 2b), 为简便表达, 设(n, m)道的相邻道隶属于集合w, 将两种概率分布进行融合, 得:

$ \begin{array}{l} {P^{{\rm{ merge }}}}({r_{n,m}}|{s_{n,m}},{s_w},{r_w}) = \left\{ {\left[ {\prod\limits_{({n^\prime },{m^\prime }) \in w} {{P^{{\rm{est}}}}} (\left. {{r_{n,m}}} \right|} \right.} \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {\left. {{r_{{n^\prime },{m^\prime }}},D_{({n^\prime },{m^\prime }),(n,m)}^2)} \right]{P^{{\rm{obs}}}}({r_{n,m}}|{s_{n,m}})} \right\}\eta \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{{k^\prime } = 1}^{{K^\prime }} {\tau _{{k^\prime }}^{{\rm{ merge }}}} N({r_{n,m}}|\mu _{{k^\prime }}^{{\rm{ merge }}},{(\sigma _{{k^\prime }}^{{\rm{ merge }}})^2}) \end{array} $

(9)

式中; η为归一化参数; K′为融合后概率分布的高斯分布数目。

显然, P^merge拥有的高斯分布数目远大于P^est和P^obs拥有的高斯分布数目。为运算简便, 采取近似方法, 即在P^merge中选取具有最大的K″个权值的高斯分布, 重新归一化后组合, 得到(n, m)道储层参数r_{n, m}的概率分布bel(r_{n, m}):

$ {\rm{bel}} ({r_{n,m}}) = \sum\limits_{k = 1}^{{K^{\prime \prime }}} {w_{{k^{\prime \prime }}}^{{\rm{bel}}}} N({r_{n,m}}|\mu _{{k^{\prime \prime }}}^{{\rm{bel}}},{(\sigma _{{k^{\prime \prime }}}^{{\rm{bel}}})^2}) $

(10)

以井位处作为起始位置, 根据二维空间外推方案(图 2b), 以及公式(6)、公式(8)、公式(9)及公式(10)求取相邻道储层参数概率分布, 再以这些相邻道为起始位置重复以上过程, 直到完成全三维工区地震道的计算, 即完成了全工区储层参数概率预测。单井储层参数概率预测流程如图 3所示。

实际工区有多口井, 往往会产生多个不同初始点的全工区储层参数概率预测结果。概率理论下的相邻道估计过程中会不断积累误差, 为了在减少误差的同时融合不同初始点的反演结果, 根据预测道与各初始位置的距离进行加权平均, 得到加权后的预测道储层参数概率分布bel^avg(r_{n, m}):

$ {\rm{be}}{{\rm{l}}^{{\rm{arg}}}}({r_{n,m}}) = \frac{{\sum\limits_{i = 1}^{{N_{{\rm{log}}}}} {{{\rm{e}}^{ - d_{w_{{\rm{log}}(n,m)}^i}^2}}} {\rm{be}}{{\rm{l}}^{w_{{\rm{log}}}^i}}({r_{n,m}})}}{{\sum\limits_{i = 1}^{{N_{{\rm{log}}}}} {{{\rm{e}}^{ - d_{w_{{\rm{log}}(n,m)}^i}^2}}} }} $

(11)

式中:N_log为井数目; d_{wlogⁱ(n, m)}为第i口井与(n, m)道之间的距离; belw_logⁱ(r_{n, m})为以第i口井作为初始位置外推得到的(n, m)道储层参数概率分布。

选取(n, m)道加权后概率分布bel^avg(r_{n, m})的期望作为(n, m)道最终的储层参数r_{n, m}^predict:

$ r_{n,m}^{{\rm{ predict }}} = E({\rm{be}}{{\rm{l}}^{{\rm{avg}}}}({r_{n,m}})) $

(12)

2 理论模型试验

设计的理论模型如图 4a所示, 包括砂泥岩薄互层以及上、下两套围岩, 其中目的层厚度50m。砂体分为类菱形、类平行四边形以及凹凸形3类, 地震波速度2350m/s, 密度2.1g/cm³; 泥岩地震波速度2420m/s, 密度2.2g/cm³; 围岩地震波速度2500m/s, 密度2.0g/cm³。图 4b为采用主频为50Hz的零相位雷克子波进行正演模拟的地震数据。

对于该地震数据(图 4b), 首先根据顶、底界面提取了振幅类、瞬时类、频谱类等八大类60种地震属性^[26], 同时采用粗糙集(RS)^[27]、三参数敏感属性分析(SDC)^[28]等方法, 联合挑选出7种敏感地震属性。然后随机选取其中30%, 20%, 10%以及5%道数据作为训练数据, 最终得到如图 5至图 8所示的不同训练样本数下砂体厚度(d)的预测结果及绝对误差曲线。对比本文方法(红色曲线)、BP神经网络(蓝色曲线)以及混合密度网络(绿色曲线)这3种方法的预测结果可以看出, 随着训练样本的减少, 3种方法预测曲线越来越偏离实际厚度曲线。但是在相同条件下, 本文方法与实际厚度曲线更为接近, 横向过渡更为连续。

考虑到单次随机选取训练数据不具有代表性, 进行50次训练数据的Monte Carlo随机挑选, 进行3种方法对比预测试验。选择预测砂体厚度曲线与实际砂体厚度曲线的平均误差和相关系数作为预测效果的评价指标。采用3种方法分别计算50次预测试验的平均误差和平均相关系数, 结果如图 9所示。

由图 9a可见, 随着训练样本的减少, BP神经网络预测的平均误差由0.6m上升到3.1m;混合密度网络平均误差由1.3m上升到3.9m;本文方法平均误差由0.2m上升到2.1m。由图 9b可见, 随着训练样本的减少, BP神经网络平均相关系数由0.98下降至0.83;混合密度网络平均相关系数由0.89下降至0.67;本文方法平均相关系数由0.99下降至0.92。综合图 5至图 9可见, 在相同训练样本下, 单纯混合密度网络预测效果最差, BP神经网络较好, 本文方法最好。

3 实际数据应用

基于上述实验结果, 将本文方法应用于X工区某实际地震数据。从工区提取的多种属性中, 同样通过粗糙集(RS)、三参数敏感属性分析(SDC)等方法联合优选出8种敏感属性(图 10至图 13)。将本工区的56口井分为41口训练井(图 10至图 13中含Tr前缀井)和15口测试井(图 10至图 13中含Te前缀井)。

在训练井位置处抽取地震道敏感属性与训练井组合成训练集, 训练BP神经网络和混合密度网络, 将训练完成的网络应用于整个工区, 得到工区的砂体厚度分布预测结果(图 14a和图 14b)。再采用本文方法对整个工区进行预测, 得到工区的砂体厚度分布(图 14c)。从图 14可以看出, BP神经网络和混合密度网络预测结果砂体厚度分布较为杂乱, 横向不连续, 而本文方法预测的砂体厚度横向更为连续, 分布集中。

选取该工区测试井位处砂体预测厚度与已知的测试井位处的砂体厚度值进行对比, 其结果如图 15所示。统计各个测试井数据, BP神经网络预测的平均误差为1.08m, 误差大于1m的井有7口, 最大误差为3.6m, 位于井Te2处; 混合密度网络平均误差为1.59m, 误差大于1m的井有7口, 最大误差为5.8m, 位于井Te13处; 本文方法平均误差为0.7m, 误差大于1m的井有3口, 最大误差为2.5m, 位于Te3处。由此可见本文方法预测的误差最低。

综上所述, 综合判断本文方法预测结果能够反映工区砂体厚度的分布情况。

4 结论

传统神经网络可视为一个复杂的、非线性的函数, 即每个输入都有确定的输出与之对应。但是利用地震数据估计储层参数具有明显的不确定性, 因而传统神经网络不能很好适应。对地震预测的不确定性问题, 最好使用概率分布来间接解决。结合了神经网络和概率分布优势的混合密度网络, 虽然从理论上能很好地解决多解性, 但实际上的表现却不如BP神经网络, 这是由于混合密度网络需要训练的参数过多, 以及概率分布本身所具有的不确定性造成的。

本文提出的基于空间近似概率约束的混合密度网络砂体厚度预测方法, 结合了混合密度网络和横向渐变假设, 借鉴地质统计学从井出发的观点, 实现了概率分布下的储层参数最优估计。通过理论数据和实际资料的应用测试, 对比了BP神经网络、混合密度网络以及本文方法的预测效果。结果表明, 本文方法具有更高的精确性和较大的应用潜力。

本文方法中, 地震道储层参数概率分布需要由各已知井外推得到, 而各已知井外推是相互独立的运算, 所以需要较大次数的迭代运算才能实现所有已知井的外推, 如何实现多井同步外推是下一步研究的重点。