0 引言

曲流河砂体是陆相沉积中最主要的油气储集体之一。伴随河道形态变化和反复沉积、侵蚀，平原区曲流河河谷内残留点坝复合沉积体^[1-3]。识别曲流河点坝对于研究曲流河演化特征、指导油气藏开发等具有重要意义^[4-5]。然而，曲流河往往经历多期改造而形成复合微相单元，如何利用地震资料准确预测砂体厚度的变化并识别复合体内部的单一点坝，是曲流河储层预测的关键。围绕上述问题形成了一系列的地震储层预测技术(如地层切片^[6-10]、地震属性分析^[11-13]及地震反演^[14]等技术)，其中地震属性分析是最简单、有效且发展迅速的方法。新的单一属性类型的研发和多属性融合是地震属性分析技术发展的两大方向。近年来，以分频属性融合技术^[15]、基于机器学习的多属性分析^[16-18]等为代表的地震属性融合分析技术广泛应用于河流相储层的研究。

极限学习机(ELM)是一种基于前馈神经网络构建的机器学习方法，它优化了隐层节点的调整，采用随机神经元和随机连接，在学习速度上显著优于支持向量机(SVM)和深度学习方法^[19-20]，而且泛化能力强，预测误差小^[21]。陈芊澍等^[22]和Wang等^[23]利用极限学习机算法进行地震多属性分析，在煤层及裂缝带预测方面取得了良好的效果，但是在储层研究领域的其他方面的应用尚未见文献报道。

本文选取渤海湾盆地济阳坳陷孤东油田馆陶组高弯度曲流河，将极限学习机算法应用于地震多属性融合中以预测砂体分布，按照“井震结合，动(生产动态资料)静(地震、地质等资料)结合”的思路识别曲流带内复合点坝沉积体中的单一点坝，揭示高弯度曲流河内部结构。

1 研究区概况

孤东油田紧邻黄河入海口，是渤海湾盆地济阳坳陷沾化凹陷内的一个大型披覆背斜构造整装油藏，主要的油气层为新近系馆陶组。研究区发育典型的高弯度曲流河沉积^{[9, 24-25]}，但近年来的新井钻探和注水开发成果揭示曲流河点坝沉积微相内部具有强非均质性，其中单一点坝认识不清成为当前制约油气开发的关键因素。研究区面积为2.45km²，共完钻井60余口(图 1)。三维地震资料覆盖全区，道间距为25m。目的层平均深度约为1260m，地震纵波平均速度为2150m/s。目的层段地震资料主频为35Hz，垂向分辨率(按照1/4λ计算)约为15.4m。

2 方法原理

每种地震属性都包含了复杂的地质信息，但是在理论上地震属性与地质参数之间缺乏确定性对应关系。从多种地震属性中获得需要的地质信息是地质和地球物理人员共同的目标。在这种大信息量、非线性的分析中，神经网络具有独特的优势。本文采用基于前馈神经网络构建的极限学习机算法进行地震多属性分析，实现多属性融合，预测砂体厚度及分布。

2.1 极限学习机算法

极限学习机是Huang等^[26]提出的一种前馈神经网络，由输入层、隐含层和输出层组成(图 2)。

在利用极限学习机进行多属性分析时，将来自不同井点的数据形成样本集合。例如对于N个不同样本形成集合{(X_j, t_j)}，其中，j=1, 2, …, N，X_j=x_j1, x_j2, …, x_jn^T∈ Rⁿ，t_j= t_j1, t_j2, …，t_jm^T∈ R^m。当隐含层节点数量为L个时，激励函数为g(x)的前馈神经输出O_j可表示为

$ \sum\limits_{i = 1}^L {{\mathit{\boldsymbol{\beta }}_i}g({\mathit{\boldsymbol{W}}_i} \bullet {\mathit{\boldsymbol{X}}_j} + {b_i}) = {\mathit{\boldsymbol{O}}_j}} $

(1)

式中：W_i= w_i1, w_i2, …w_in^T表示输入层到第i个隐含层节点的输入权值；b_i为第i个隐含层节点的偏差；β_i= β_i1, β_i2, …β_im^T是连接第i个隐含层节点与输出层的权值。

网络的学习目标是使输出的误差最小，即$ \sum\limits_{i = 1}^L {\left\| {{\mathit{\boldsymbol{O}}_j} - {\mathit{\boldsymbol{t}}_j}} \right\| = 0} $，则式(1)可表示为

$ \sum\limits_{i = 1}^L {{\mathit{\boldsymbol{\beta }}_i}g({\mathit{\boldsymbol{W}}_i} \bullet {\mathit{\boldsymbol{X}}_j} + {b_i}) = {\mathit{\boldsymbol{t}}_j}} $

(2)

可简化为

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

(3)

其中

$ \begin{array}{c} \mathit{\boldsymbol{H}}({\mathit{\boldsymbol{W}}_1}, \cdots {\mathit{\boldsymbol{W}}_L}, {b_1}, \cdots , {b_L}, {\mathit{\boldsymbol{X}}_1}, \cdots {\mathit{\boldsymbol{X}}_L}) = \\ {\left( {\begin{array}{*{20}{c}} {g({\mathit{\boldsymbol{W}}_1} \bullet {\mathit{\boldsymbol{X}}_1} + {b_1})}& \cdots &{g({\mathit{\boldsymbol{W}}_L} \bullet {\mathit{\boldsymbol{X}}_1} + {b_L})}\\ \vdots &{}& \vdots \\ {g({\mathit{\boldsymbol{W}}_1} \bullet {\mathit{\boldsymbol{X}}_N} + {b_1})}& \cdots &{g({\mathit{\boldsymbol{W}}_L} \bullet {\mathit{\boldsymbol{X}}_N} + {b_L})} \end{array}} \right)_{N \times L}} \end{array} $

$ \mathit{\boldsymbol{\beta }} = [\mathit{\boldsymbol{\beta }}_1^{\rm{T}} \cdots \mathit{\boldsymbol{\beta }}_L^{\rm{T}}]_{L \times m}^{\rm{T}} $

$ \mathit{\boldsymbol{T}} = [\mathit{\boldsymbol{T}}_1^{\rm{T}} \cdots \mathit{\boldsymbol{T}}_L^{\rm{T}}]_{N \times m}^{\rm{T}} $

在极限学习机中，输入权重W_i和隐含层的偏差b_i(由算法随机给定)，使隐含层的输出矩阵H被唯一确定。因此，训练问题转化为求解输出权值矩阵的最小二乘解的问题，输出权值矩阵为

$ \mathit{\boldsymbol{\hat \beta }} = {\mathit{\boldsymbol{H}}^{† }}\mathit{\boldsymbol{T}} $

(4)

式中${\mathit{\boldsymbol{H}}^{† }} $为H的Moore-Penrose广义逆。

在计算过程中，只需调整隐含层神经元结点数，不需要迭代。因此，相比传统基于梯度下降的迭代方法，极限学习机算法能够解决收敛速度慢、容易陷入局部极小值的问题，具有极端“快速”的特点^[27]，可明显缩短学习时间；同时，所得解为唯一最优解，网络的泛化性能得到了保证^[28]。

2.2 基于极限学习机的多属性融合

将地震属性分析技术与机器学习方法相融合，可预测砂体厚度。其中有两个关键问题：一是如何选取地震属性；二是如何融合地震属性预测砂体。为此，本文设计的解决路径包含四个部分(图 3)，即地震属性提取与优选(流程A)、特征变量与目标变量的统计(流程B)、地震属性与地质参数关系的建立(流程C)及砂体厚度的地质解释(流程D)。其中流程A和B针对第一个关键问题，解决地震属性的选取问题；流程C针对第二个关键问题，解决属性的融合与地质参数的预测问题。

需要特别说明的是，在流程A中，需要对数据进行预处理和属性优选，具体分为以下四步。

(1) 单一地震属性提取。为了减少属性类型选取的盲目性，需要根据研究区砂体发育特点(厚度、组合模式、沉积相类型等)和地震属性的物理意义进行选择。需要强调的是，充分利用前人的经验是一个有效的途径。

(2) 地震属性标准化处理。不同地震属性物理意义和值域存在差异，开展属性融合分析之前，需要对训练集和测试集中的地震属性进行标准化处理，以避免数值较大的地震属性对融合结果产生过度影响。

(3) 单一属性与井点砂岩厚度之间的相关性分析及筛选。计算每种地震属性与井点砂岩厚度之间的相关性，筛除相关性差的地震属性类型。

(4) 属性聚类分析进行再次筛选。一些在物理意义上相近的地震属性与砂岩厚度之间相关性都较高时，容易形成信息重复。为了避免这一情况，实现数据有效降维，需要对上一步筛选后的地震属性进行聚类分析，在聚类结果中选择与砂岩厚度相关性最高的属性。

经过上述属性处理与优选后，使用极限学习机算法，构建多属性融合预测砂岩厚度的网络模型，最后将全区的地震属性分析数据作为特征变量输入最优模型，预测砂体厚度及分布。

3 应用实例 3.1 模型正演

根据河道的迁移特征，曲流河道间的横向组合可以划分为两种类型、三种模式，即点坝起始端与废弃河道叠加、废弃河道与废弃河道叠加两大类。其中前一类又可进一步细分为废弃河道泥岩保留和不保留两种模式。

对上述三种模型开展正演模拟(图 4)。根据井点信息，设置点坝主体位置砂岩厚度为12m，废弃河道处砂岩厚度为7m，砂岩速度为2160m/s，泥岩速度为2500m/s，子波为35Hz雷克子波。对正演结果进行-90°相位转换。

模式一(图 4a)为两期河道同向迁移叠加，后期河道起始端紧邻前一期并部分侵蚀前一期废弃河道，但前期废弃河道泥岩保留。在点坝叠置部位砂体厚度快速减小，正演模拟结果可见振幅快速减弱。

模式二(图 4b)同样为两期河道同向迁移叠加，但后期河道完全侵蚀前一期的废弃河道，甚至侵蚀部分河道砂岩，前一期废弃河道细粒沉积部分没有残留。在这种情况下，点坝叠置位置砂体厚度减小速度慢、程度弱，振幅较两侧相对减弱。

模式三(图 4c)为距离较远的两期河道相向迁移，后期河道侵蚀了前期河道而在两条河道之间保留了后期的废弃河道。点坝叠置部分砂体厚度由两侧向中间快速减小，振幅快速减弱至低值。

三种叠置模式的正演响应特征的相同点是在点坝叠置部位振幅减弱，不同点是砂体减薄方式导致的振幅减弱程度的差异和波形错动的位置。从振幅变化可以看出，振幅与砂岩厚度变化有良好的对应关系，这就意味着振幅能够较好地反映砂体的厚度变化。

3.2 多属性融合过程

对目的层段提取振幅类、频谱类、相位类等共计16种常见的地震属性，计算各属性与砂体厚度的相关系数，初步优选出6种与砂体厚度相关系数大于0.5的地震属性(表 1)。对6种初选地震属性聚类分析，结果显示可为三大类(图 5)。平均能量(AE)、能量自相关(EOAC)、振幅平方差(VA)、均方根振幅(RMS)四者距离在2以下，为第一类；弧长(ARC)和最小振幅值(MV)分别作为第二类和第三类，它们与第一类距离大于5且彼此距离也大于5。选取每类中与砂体厚度相关系数最高的地震属性，即最终选择RMS、ARC和MV属性作为参与下一步机器学习的特征变量来源。

3.3 预测效果 3.3.1 砂体厚度预测

选取拟合优度(R²)和均方根误差(RMSE)作为评价预测结果准确性的参数。拟合优度越接近于1则模型预测效果越好；均方根误差越接近于0则说明预测准确性越高。

训练集多属性融合预测与实际砂体厚度的R²为0.8237，RMSE为0.5813m(图 6a)；单一属性(RMS)预测与实际砂体厚度的R²为0.6562，RMSE为0.9546m(图 6b)。测试集预测与实际砂体厚度R²为0.8275，RMSE为0.5921m(图 6c)；单一属性(RMS)预测与实际砂体厚度的R²为0.6787，RMSE为0.9906m(图 6d)。上述结果表明，基于极限学习机的多属性融合预测结果比单一属性预测结果准确度更高。

3.3.2 异常区预测

研究区大部分钻井揭示，均方根振幅单一地震属性与砂岩厚度之间具有较高的相关性，但局部存在与实钻不吻合的异常区，这降低了单一属性预测砂岩厚度的可靠性，因而需要采用多属性预测砂岩厚度。

研究区36-255井区的6口井存在地震属性与砂岩厚度关系异常(图 7红色虚线框内)。钻井揭示6口井钻遇点坝主体，砂体厚度均为10m左右。单一地震属性表现为低值，根据单一地震属性预测为砂岩不发育区(图 7b)，这与实钻结果不符。采用极限学习机多属性融合预测的6口井的砂岩厚度与实际砂岩厚度的R²为0.6639，RMSE为0.2961m(图 8a)。而单一的均方根振幅属性预测的砂岩厚度与实际砂岩厚度R²为0.4104，RMSE为1.0987m(图 8b)。

如上所述，极限学习机驱动的多属性预测方法明显优于单属性预测方法。

3.3.3 不同机器学习方法对比

选取研究区资料齐全的30口井建立对照组，比较ELM、BP神经网络和SVM算法进行多属性融合预测的效果。以砂岩厚度预测误差10%作为标准，预测误差小于10%为吻合，大于10%为不吻合，统计三种不同方法的砂岩厚度预测单井吻合率。结果(表 2)显示，ELM方法的单井吻合率为93.3%，SVM算法为86.7%，而BP神经网络预测单井吻合率为83.3%，这表明极限学习机驱动的多属性融合方法(本文方法)效果更好。

3.3.4 单一点坝识别

应用本文方法预测的砂岩厚度及分布(图 9)可以确定在泥岩沉积背景下发育的条带状河道，厚度大值区对应的点坝复合体被砂体不发育的“月牙形”(泥岩发育带)分割。结合钻井标定可以识别出五个单一点坝，即点坝A~点坝E。

测线LL′(图 10a)上可见侧积方向相同的两个单期叠置点坝(点坝D和点坝B)，前一期废弃河道完全被后期河道侵蚀，并且没有残留，属于模式二(图 4)。37-254井两侧反射同相轴连续性较好，但是该井东侧的反射振幅更强(图 10a)。钻井显示两期砂体叠置(图 10b)，与地震反射的振幅变化一致。

测线MM′(图 11a)上可见两个侧积方向相对的单期叠置点坝，即点坝C与点坝A(属于图 4中的模式三)。38-2294井位于点坝叠置区(图 9)，处于弱振幅区，反射同相轴略向下凹，显示了时间厚度减小的特点(图 11a)。这是由于该井所在位置为废弃河道叠置区、砂岩厚度减小造成的。从连井剖面(图 11b)可见，38-2294井砂岩厚度减小且砂体顶面低于两侧其他井，说明了这是河道废弃造成的。

测线NN′(图 12a)上可见三个侧积方向大致相同的单期点坝(点坝D、点坝B和点坝A)，36-275井和38-295井处于单期点坝的叠置区。测井资料揭示发育废弃河道残留的特征(图 12b)，该点坝组合属于模式一(图 4)。前述两井点处的振幅变化以及同相轴揭示的高程关系表明，自东向西三个单一点坝发育时间由老逐渐变新，这与钻井结果一致。

综上所述，研究区发育高弯度曲流河的点坝复合体，北部和西南部发育泛滥平原，点坝复合体内部被残留的废弃河道将曲流带分为五个单一点坝(图 13a)。其中，点坝C与点坝A侧积方向相对；点坝A、点坝B和点坝D侧积方向相近，点坝B末期的废弃河道泥岩部分保留。点坝C、点坝D和点坝E的废弃河道在平面上连续，发育稳定，结合现代河流沉积特征(图 13b)可判断为同期河道沉积的三个连续点坝。

新钻38-C295井是38-295井的侧钻井，在F靶点钻遇目的层，在G靶点出目的层。本文砂岩厚度预测及单一点坝解释结果表明，该井目的层段处于废弃河道内(图 14a)。随钻自然伽马曲线(图 14b)显示38-C295井目的层顶部发育厚度3.00m的泥质沉积，砂岩厚度为8.10m，而平面上与之紧邻的37-274井砂岩厚度为10.24m，从而揭示了38-C295井在目的层顶部钻遇废弃河道，验证了单一点坝的地震解释结果的合理性。

开发动态资料可以反映井间砂体的连通性，可以检验储层内部结构的预测结果。图 14a所示的三口井位于同一个开发井组，其中37-274井为注水井，38-275井和38-295井为开发井，二者与注水井37-274井的距离均约为250m。根据本文方法单一点坝识别结果，37-274井与38-275井位于点坝B内，而38-295井则位于点坝A(图 9，图 14a)。注采资料显示，37-274井注水量的变化与38-275井的产液量有良好的对应性。例如，1999年4月至2000年2月，随着37-274井注水量的降低，对应的采油井38-275井的液量发生了降低；而2001年2月至2002年5月，37-274井注水量经历了先调减后调增的措施变化，38-275井的产液量出现了一致性的变化，这两口注采对子井的变化趋势和时间具有一致性。而同期注水井37-274井与邻近的采油井38-295井的产液量变化并无对应性(图 15)。从生产动态信息来看，37-274井与38-275井的井间砂体连通性要远远好于37-274井与38-295井，这从砂体井间连通性方面验证了本文方法单一点坝识别结果的准确性。

4 结论

(1) 基于极限学习机的多属性融合技术比单一属性对砂体厚度的预测更准确，在实例区的单井吻合率达到93.3%，高于支持向量机和BP神经网络算法。

(2) 孤东油田馆陶组曲流带发育三种单一点坝组合模式，即相向迁移点坝叠置模式、同向迁移废弃河道残留模式和同向迁移废弃河道无残留模式，三种模式的地震反射差异主要为横向连续性的变化和振幅幅度的差异。