基于半监督学习的井震联合储层横向孔隙度预测方法

韩宏伟， 刘浩杰， 桑文镜， 魏国华， 韩智颖， 袁三一*

1 中国石化集团公司胜利油田物探研究院， 东营 257000 2 中国石油大学(北京)油气资源与探测国家重点实验室， 北京 102249

0 引言

孔隙度是描述储层岩性、物性和含油气性特征的重要参数之一.准确预测孔隙度有利于开展甜点区域优选、油气储量估算、井网井位设计、油藏数值模拟等工作(Angeleri and Carpi，1982；Doyen，1988；Leite and Vidal，2011).根据不同的地球物理探测和观测技术，孔隙度预测方法主要包括岩石物理实验测量、测井解释和井震联合预测三类.岩石物理实验测量常采用氦气法和高压压汞法等测得岩心孔隙度(韩学辉等，2021).测井解释主要采用经验公式、三孔隙度模型、交会图、岩石物理模型、神经网络等方法构建孔隙度与密度、补偿中子、声波时差、自然伽马、泥质含量、地层电阻率等敏感参数之间的计算模型(Archie，1942；Xu and White，1996；Khaksar and Griffiths，1998；李宏兵等，2019；安鹏等，2019).这两类方法能实现精确测量岩心段或全井段孔隙度，但无法有效外推获得井间孔隙度.

井震联合储层孔隙度预测的理论基础是地震波在地层中的传播速度与岩石孔隙度和泥质含量等有关(Wyllie et al.，1956；Klimentos and McCann，1990).因此，整合地震和测井数据在数据覆盖范围、纵横向分辨率等方面的各自优势，能较好地预测出储层孔隙度的空间分布.井震联合预测方法主要是基于不同的地震属性直接或间接地预测孔隙度.根据利用的地震属性和实现过程的差异，先后发展出叠后单属性预测、叠后多属性融合、叠前多属性融合三类孔隙度预测方法.叠后单属性预测方法一般是基于经验公式和岩石模型等方式进一步线性或非线性地转化叠后波阻抗反演结果(或速度及声波时差)为储层孔隙度(Angeleri and Carpi，1982；de Buyl et al.，1988).或考虑地震垂直走时和储层孔隙度的横向变化的空间相关性，采用协同克里金等地质统计学技术进行储层建模并生成孔隙度模型(Doyen，1988；Moon et al.，2016).叠后多属性融合预测方法采用神经网络、随机模拟等技术直接映射优选的地震属性组合为储层孔隙度.Hampson等(2001)通过多元逐步线性回归分析优选与预测目标参数相关性最好的地震属性组合，并结合两个应用实例说明概率神经网络相比于多元线性回归能更好表达地震属性与孔隙度之间的非线性关系，进而提高预测的可靠性.Leite和Vidal(2011)研究了一种基于神经网络融合自然伽马数据体和递推反演得到的阻抗数据体的三维孔隙度预测方法.

相比于叠后地震数据，叠前地震数据的振幅随偏移距变化(AVO)特征包含着更为丰富的储层岩性、物性与流体信息.因此，理论上叠前数据预测储层物性参数有利于提高预测准确性和可靠性，降低多解性.叠前多属性融合预测方法依托叠前同步反演和多属性融合两大技术，通过叠前反演获得与孔隙度有较好相关性、包含更多反映储层物性特征的叠前弹性参数(如纵横波速度、纵横波速度比、纵横波阻抗、密度、杨氏模量、泊松比等).随后优选并组合叠前AVO属性(如截距、梯度等)、叠后地震属性(如振幅、频率、波形三大类等)和叠前弹性参数中的一类或多类参数，利用神经网络和随机模拟等多属性拟合技术建立基于优选属性的孔隙度预测模型(AlMuhaidib et al.，2012).韩宏伟等(2021)系统地总结了地震储层预测技术表征孔隙度等储层参数的原理与应用案例.

现有的井震联合储层孔隙度预测方法要解决的核心问题是如何优选出与孔隙度相关的地震属性以及提高弹性参数预测精度.人工智能技术可以间接解决上述两个核心问题，为从地震数据中自动提取、优化和组合地震属性来直接预测孔隙度提供了新的可行性途径.近几年来，深度学习在地球物理勘探领域的应用逐步从地震资料处理、解释延伸到弹性参数反演、储层参数预测等方面(Yuan et al.，2018，2022；Sang et al.，2021；Song et al.，2021；王迪等，2021).以卷积神经网络、递归神经网络和生成对抗网络为主的不同类型深度神经网络也开始逐渐应用于储层孔隙度等物性参数预测任务.目前，这类智能建模方法在多测井曲线拟合井点位置的孔隙度方面研究较多，在井震联合直接预测储层横向孔隙度方面研究相对较少.在预测孔隙度曲线方面，考虑到测井数据是反映不同地质时期地层沉积特征的局部相关的时间序列，众多学者采用循环神经网络建立多种敏感的测井参数转换为孔隙度的非线性映射，以实现孔隙度建模与预测(安鹏等，2019；宋辉等，2019；王俊等，2020；Chen et al.，2020).相比于全连接神经网络、多元线性回归分析等方法只考虑来自同一深度点(或时间采样点)的测井曲线之间的相互联系，考虑测井曲线渐变沉积规律和前后关联的循环神经网络具有更高的预测精度和鲁棒性.在井震联合储层横向孔隙度直接预测方面，Das和Mukerji(2020)采用端到端卷积神经网络进行有监督学习，从叠前地震数据预测出孔隙度.Feng等(2020)整合卷积神经网络与褶积模型搭建物理约束的混合网络，模拟地震反演和正演过程开展无监督孔隙度预测.该方法不需要大量的标签数据，适合缓解地震测井数据不均衡引起的油气小样本问题.但是该方法依赖于准确的低频信息和震源子波信息参与混合网络的训练.

为减少低频模型对神经网络预测孔隙度的影响，并缓解基于有监督学习范式下的智能化孔隙度预测结果横向连续性差的问题，本文提出一种基于半监督学习的井震联合储层横向孔隙度预测方法.该方法基于双向门控递归单元(Bi-GRUs)构建由编码网络和解码网络组成的孔隙度预测混合网络(Cho et al.，2014).其中，编码网络模拟地震反演及岩石物理建模过程，转化输入叠后地震数据为孔隙度；解码网络建立数据驱动的正演模型，实现孔隙度非线性映射为生成地震数据.此外，非井旁地震道在每次迭代过程中随机引入混合网络参与网络训练.通过整合井位置的地震数据和(伪)测井数据及非井位置的地震数据开展半监督学习，有望既保证预测孔隙度的精度，同时通过地震波形匹配间接保证孔隙度结果的横向连续性，获得高精度的储层孔隙度模型.

1 方法原理

本小节首先回顾基于有监督学习构建孔隙度预测网络的基本思想.之后介绍本文提出的半监督井震联合孔隙度预测网络的主要原理，包括网络架构、目标函数和实现过程等.最后介绍合成数据集的生成过程，该数据集用于后续测试的网络智能建模与结果定量评价.

1.1 有监督孔隙度预测网络

随着以深度学习为代表的人工智能技术逐渐用于求解高分辨率处理、波阻抗反演、层析反演、全波形反演、弹性参数反演等反问题，储层物性参数智能反演也开始受到重视(Araya-Polo et al.，2018；Kim and Nakata，2018；Alfarraj and AlRegib，2018，2019；Biswas et al.，2019；Gu et al.，2019；Russell，2019；Sun and Alkhalifah，2020；Aleardi and Salusti，2021).目前，通过学习由地震数据和目标参数构成的训练样本对建立有监督储层参数预测网络是一种常用的储层参数智能反演方法(Alfarraj and AlRegib，2019).在有监督学习框架下，本质上孔隙度等物性参数的智能反演与其他储层参数反演类似.其实现的数学基础是理论上神经网络能直接建立地震数据到任意目标参数的非线性映射关系，通过梯度下降迭代优化目标函数而逼近目标解；实现的物理基础是不同类型的地球物理数据从不同角度刻画地下相同目标地质体，它们之间具有内在的物理关联和相似的响应特征.

地层通常遵从渐变的沉积规律，局部地层性质是近似且有关联的.因此，连续记录的地震或测井数据在不同时刻或不同深度获得的响应与其附近响应存在内在关联和局部相似，综合反映了周围地层的性质.基于以上分析，地震记录和测井曲线都可以视为内部关联、局部依赖的序列数据.在语音识别和文本分类等自然语言处理场景中，循环神经网络(RNN)已被广泛应用于解决序列建模问题(Hochreiter and Schmidhuber，1997；Chorowski et al.，2014).在地震勘探领域的弹性参数反演和物性参数预测等方面，不少学者通过比较前馈神经网络和循环神经网络的预测性能发现，后者更适合准确预测储层参数的变化趋势和局部细节(安鹏等，2019；王俊等，2020).相比于前馈神经网络，RNN能记忆上一时刻隐藏层学到的信息，并利用该历史信息影响当前时刻隐藏层的输出.但是，常规RNN由于不同时刻信息的长距离依赖容易出现梯度消失或梯度爆炸的问题, 导致其难以训练获得较精准的网络模型.而双向门控递归单元(Bi-GRUs)通过更新门和重置门灵活决定前一时刻隐层信息的记忆或遗忘，缓解常规RNN模型难以优化更新的问题.此外，其能避免常规RNN信息状态的单向传递，采用两个方向相反的单向门控递归单元(GRU)实现信息的正向和反向传播，使得当前时刻的输出由当前时刻之前和之后的隐层信息共同决定.这一工作机制使得其在预测某一时间点或深度点的孔隙度时能考虑局部相邻的地震或测井模式，获得更能反映地层变化规律的孔隙度结果.因此，本文采用Bi-GRUs作为构建孔隙度预测网络的基本单元.以GRU为例，其内部细胞结构是由重置门、更新门和tanh激活函数三者构成的.重置门和更新门的输出值都与当前时刻的输入和上一时刻的隐层输出有关，且更新门输出值越小，历史状态信息对当前时刻GRU输出结果的影响程度越小.由于Bi-GRUs的工作原理不是本文研究重点，更多的原理介绍参考Cho等(2014)、Alfarraj和AlRegib(2019)和Yuan等(2022).

有监督孔隙度预测网络(图1a，以下简称为有监督网络)一般包括输入层、隐藏层和输出层三个部分.输入层和输出层分别是单道地震记录和网络预测的孔隙度曲线.为与后面的半监督孔隙度预测网络形成对比，中间隐藏层由4个双向门控递归单元和1个全连接层构成.有监督网络主要通过成对的地震记录和孔隙度曲线进行学习和训练，该过程可以视为一种广义的联合叠后地震反演与岩石物理模型的建模流程，即模拟地震数据到弹性参数，弹性参数到物性参数的物理反演过程，实现地震数据直接转化为孔隙度.有监督网络(简记为E)通过反向传播算法不断更新网络参数，减小网络预测结果与孔隙度标签之间的差异，其差异使用如下目标函数LE计算测井监督损失：

(1)

其中，N为训练样本或标签的数量，Si和Pori分别为第i个地震道曲线、第i个孔隙度曲线，θE为有监督网络E的网络参数，E(Si;θE)为网络E预测的孔隙度结果.

通过迭代最小化目标函数，最终得到的有监督孔隙度预测模型能实现从单道地震记录预测出最佳的单道孔隙度曲线.但是，由于该逐道反演过程没有考虑储层的横向变化，特别是在预测的单道孔隙度结果出现较大偏差时，预测孔隙度剖面的横向连续性差.

1.2 半监督孔隙度预测网络

为缓解常见的有监督学习孔隙度预测方法存在着信息利用不足、预测结果连续性差等问题，本文提出一种基于半监督学习的孔隙度预测方法.该半监督孔隙度预测网络(图1b，以下简称半监督网络)是由编码网络和解码网络组成的混合网络.编码网络与上一小节有监督网络的网络结构保持一致(编码网络也简记为E)，也是由4个双向门控循环单元和1个全连接层构成；而解码网络D与编码网络E采用相同的网络模块.E和D的区别在于各自的输入、输出和网络内部权重与偏置不同，E的输入是单道地震记录，输出是预测的孔隙度曲线；D的输入是E预测的孔隙度曲线，输出是生成的地震记录.以编码网络E为例，E中不同隐层的双向门控递归单元通过更新候选隐层状态和隐层状态实现从输入地震数据或上一隐层输出进一步提取高维的多水平特征，不同水平特征最终通过线性回归层加权输出得到期望的孔隙度曲线.预测的孔隙度进一步进入解码网络D正演生成地震数据.

图1 井震联合储层孔隙度预测网络(a) 有监督孔隙度预测网络； (b) 半监督孔隙度预测网络.Fig.1 Reservoir porosity prediction network via seismic and well logs integration(a) Supervised porosity prediction network； (b) Semi-supervised porosity prediction network.

实际工区稀疏分布的测井数据和稠密采集的地震数据导致智能化储层参数预测面临着数据不均衡引起的油气小样本问题.即稀疏的井点位置可以生成地震测井标签样本对；而众多的非井点位置只有无标签的地震数据.为尽可能充分利用地震和测井信息保证井间孔隙度预测效果，本文提出的半监督网络除使用有监督网络中有限的成对井旁地震道和测井标签外，还采用众多的非井旁地震道共同构建孔隙度预测网络的训练集，即联合少量的有标签数据和众多的无标签数据开展半监督学习预测孔隙度.半监督网络在每次迭代优化的过程中，首先利用少量的有标签数据融合窄带地震数据和宽带测井数据这两种不同尺度信息，同时建立地震数据到孔隙度的编码网络和孔隙度到生成地震数据的解码网络，并计算井旁地震道数据匹配损失和测井监督损失.接着进一步在当前迭代周期内随机地从地震剖面上选择选取N个非井位置的地震记录作为无标签数据依次通过编码网络和解码网络而参与网络训练与优化.非井旁地震道由于缺少真实测井标签而无法采用测井监督损失评估其对应的孔隙度预测效果, 因此只采用非井旁地震数据匹配损失间接保证非井位置的孔隙度预测结果准确而且稳定，不易受到地震振幅特征变化而出现较大偏差.井旁地震道数据匹配损失、测井监督损失和非井旁地震数据匹配损失三项构成半监督混合网络的总损失.每次迭代过程计算出总损失后再通过反向传播实现梯度更新并进入到下一次迭代，当训练集和验证集损失都收敛到极小值后，完成半监督网络的训练.由于半监督网络中编码网络与有监督网络结构相同，其目标函数也都是公式(1)；解码网络的目标函数LD为：

(2)

其中，S′i为第i个非井位置的地震记录曲线，θD为解码网络D的网络参数，D(E(Si;θE);θD)和D(E(S′i;θE);θD)分别为Si、S′i通过混合网络后生成的地震数据.解码网络目标函数中的第一、二项分别用于计算井旁地震道和非井旁地震道的地震数据匹配损失.此外，不同迭代次数下的非井旁地震道是从原始地震数据中随机选择的，以保证网络较好地适应不同波形特征预测出精确的孔隙度.混合网络的总目标函数Ltotal为：

Ltotal=LE+λLD，

(3)

其中，λ为正则化参数，用于控制地震数据匹配程度对孔隙度预测效果的影响.当λ=0，公式(3)与公式(1)等价，即有监督学习范式是半监督学习范式的一种特例.本文提出的半监督孔隙度预测方法摈弃了弹性参数反演、属性优选等复杂过程，编码网络和解码网络在地震数据匹配损失和孔隙度估计损失的联合监督下同时迭代优化网络参数，实现从地震数据直接定量预测地下储层孔隙度参数的空间分布，而且保证预测孔隙度能通过解码网络“正演”生成逼近真实的地震数据.

1.3 合成数据集生成

为证明本文提出的方法预测孔隙度的有效性，选取常用的Marmousi模型(Martin et al.，2006)进行测试.图2a为通过速度模型导出的孔隙度模型，其模拟地层的岩性以泥岩为主，还含有部分砂岩、泥灰岩和盐岩.模型中部从浅到深依次经过三条大断裂、一个泥灰岩背斜、一个不整合面、一个楔形真空盐岩层和一个泥岩背斜地层.图2b为阻抗导出的反射系数与震源子波正演生成的无噪合成地震数据.震源子波的主频为35 Hz，采样间隔为1 ms.图2c为在图2b基础上添加随机噪声得到的含噪合成地震数据，噪信比为0.3(噪声与地震信号的能量之比).

图2 用于孔隙度预测的合成数据集(a) 孔隙度模型； (b) 无噪合成地震数据； (c) 含噪合成地震数据.Fig.2 The synthetic dataset for porosity prediction(a) Porosity model； (b) Noise-free synthetic seismic data； (c) Noisy synthetic seismic data.

选取孔隙度模型和无噪合成地震剖面中的部分数据构建训练集和验证集，无噪和含噪合成地震剖面作为测试数据，孔隙度模型作为参考标准评价有监督和半监督孔隙度预测方法的优劣.评价指标采用相对误差(R)和相关系数(CC)两种.相对误差定义为：

(4)

其中，X和X′分别为真实孔隙度模型和预测孔隙度结果，m和n分别为孔隙度模型或预测孔隙度结果的时间采样点数和道数.相关系数定义为：

(5)

其中，X为真实孔隙度模型(或观测地震数据)，X′为预测孔隙度结果(或生成地震数据)，Cov(X,X′)、Var[X]和Var[X′]分别为X与X′之间的协方差、X的方差和X′的方差.

2 例子

2.1 数值模型测试

本小节采用Marmousi模型进行方法测试与评价.训练有监督网络时，训练集由随机抽取的8条地震记录及其相应的孔隙度曲线(图2b和图2a中的黑线)构成，验证集由等间隔抽取的6条地震记录及其对应的孔隙度曲线(图2b和图2a中的红线)组成，测试集由无噪地震数据(图2b)、含噪地震数据(图2c)和真实孔隙度模型(图2a)形成.三个数据集中每条地震记录或孔隙度曲线的大小都为750×1.训练半监督网络时，验证集和测试集与有监督网络使用的保持一致，训练集为成对的地震记录和孔隙度曲线(图2b和图2a中的黑线，即有监督网络的训练集)及非井旁地震道.在半监督网络的不同迭代训练周期内，其训练集中的非井旁地震道是利用不同种子数随机生成的，而有标签数据对始终维持不变，即用于训练半监督网络的无标签数据远多于有标签数据.此外，不同数据集中的地震数据和孔隙度都采用Z-score归一化进行预处理：

(6)

在训练数据、优化器、目标函数和反向传播算法等的共同作用下，有监督网络通过预先设定最大迭代次数完成初步学习与训练，通过绘制损失曲线观察训练损失和验证损失是否同时收敛到极小值，若没有则调整最大迭代次数直至二者都满足收敛条件.最终确定最大迭代次数设置为300次时可以获得最优的有监督孔隙度预测模型，对应的损失曲线见图3.之后推广到归一化后的无噪地震数据进行测试，网络的输出结果采用公式(6)中的均值和方差反归一化处理化后得到最终预测的孔隙度(图4a).孔隙度预测结果与真实孔隙度模型(图2a)的相关系数和相对误差分别为0.9616和2.59%，表明预测结果整体精度较高.网络单道反演预测得到的孔隙度在浅部横向连续性较好，这是因为孔隙度模型在浅部(小于0.25 s的部分)地层起伏较小，横向变化较为缓慢，且浅部倾斜地层、褶皱和断层内的孔隙度整体变化也较小.而孔隙度模型在大于0.25 s的中深部地层起伏较大，发育有陡倾角的断层和背斜等，且孔隙度在中深度部的局部差异较大.导致有监督网络在背斜等复杂构造位置推广能力相对较差，因而预测结果的横向连续性会变差，局部出现明显的预测偏差.其原因可能是对应位置的低频成分难以准确估计.可以考虑学习相对孔隙度、加入孔隙度初始模型参与网络训练、多道地震数据预测单道孔隙度等方式缓解该问题.此外，若采用不同随机种子为有监督网络设置不同的网络参数初始化状态，以此建立的多个有监督网络预测出多个孔隙度模型，取不同结果的均值作为统计性孔隙度预测结果也可以缓解该问题.但是，该方式的计算成本较高.采用传统模型驱动类方法可以获得速度或阻抗参数，再根据孔隙度与速度或阻抗的岩石物理关系可以进一步获得预测的孔隙度结果.参考She等(2018)的工作，图4b为采用基于全变分约束阻抗反演方法导出的孔隙度结果.从定量预测角度来看，有监督方法预测出更符合真实模型的孔隙度结果，其横向连续性较传统方法预测效果要差一些.其原因是传统方法使用了初始模型提供测井先验信息.但是，该传统方法由弹性参数到孔隙度的二次转换及其方法本身对地下地质结构的块状假设导致预测的孔隙度整体出现较大偏差，表现为预测孔隙度值整体偏小.

图3 有监督孔隙度预测网络的损失曲线Fig.3 Loss curves of the supervised porosity prediction network

图4 智能方法与传统方法孔隙度预测结果对比(a) 基于有监督孔隙度预测网络得到的孔隙度结果(迭代300次)； (b) 基于叠后波阻抗反演方法导出的孔隙度结果.Fig.4 Comparison of porosity prediction results between the intelligent method and the traditional method(a) Porosity prediction result via supervised porosity prediction network (300 iterations)； (b) Porosity prediction result derived from the post-stack acoustic impedance inversion method.

在与有监督网络的最大迭代次数保持一致的情况下，使用提前准备的数据集训练半监督网络，其使用(伪)测井数据和(伪)井位置及非(伪)井位置的地震数据匹配共同监督和质控网络的训练过程.目标函数(式(3))中地震数据匹配项的正则化参数λ会影响半监督网络的迭代优化方向和孔隙度预测效果.下面对比分析不同正则化参数情况下训练得到的半监督网络推广到无噪地震数据(图2b)预测的孔隙度结果与有监督网络预测孔隙度结果之间的差异，并讨论λ大小对孔隙度预测及地震数据匹配效果的影响.当λ设置为0时，半监督网络预测的孔隙度结果(图5a)与有监督网络预测的孔隙度结果(图4a)完全一致，表明此时两种网络建立的孔隙度预测模型是基本相同的.由于此时的地震数据匹配项完全没有起到任何作用，所以预测孔隙度结果(图5a)的精度及横向连续性并没有得到改善，并且其进一步通过半监督网络中的解码部分生成的地震数据(图5f)与观测地震数据(图2b)无法匹配，二者的残差(图5k)极大.当λ设置较小时(λ=0.1)，半监督网络预测得到的孔隙度结果(图5b)相较于图3的横向连续性得到一定改善，但是地震数据匹配损失占总损失的比例较小，导致半监督网络生成的地震数据(图5g)与观测地震数据(图2b)之间的残差仍然较大(图5l).当λ设置较大时(λ=1.5)，地震数据匹配过多地影响孔隙度预测效果，使得孔隙度结果(图5c)中出现了与观测地震数据(图2b)形态相似的预测假象，并且在浅部(0～0.1 s)表现尤为明显.此时，生成地震数据(图5h)相比于图5g进一步接近观测地震数据.当λ设置合理时(λ=0.37)，在地震波形控制和测井监督的共同作用下，不仅可以获得比有监督方法预测精度更高且横向连续性更好的孔隙度预测结果(图5d)，而且预测孔隙度通过数据驱动的“正演模型”(即解码网络)可以进一步实现准确的地震数据匹配(图5i).以上4组对照测试与有监督网络设置的迭代次数一致，都为300次.在确定出合适的正则化参数后，继续增加迭代次数(如增大到1000次)可以进一步提高孔隙度预测结果(图5e)的质量，同时进一步减小生成地震数据(图5j)与观测地震数据之间的差异(图5o).对比图5a、b及图4a可以看到，引入地震数据匹配约束会改善孔隙度预测结果的横向连续性，但需要设置合理的正则化参数才能保证预测准确性和连续性的同时提升.对比图5f—j及图5k—o不难发现，地震数据匹配的好坏可以在一定程度上反映预测孔隙度剖面的精度，并且通过生成地震数据也可以间接判断预测孔隙度剖面的横向连续性.这是因为当孔隙度预测结果出现较大偏差时也会导致生成的地震记录不准确，即预测孔隙度剖面与生成地震剖面的横向连续性具有较强的一致性(如图5a、f，图5b、g等).

图5 基于半监督孔隙度预测网络测试得到的孔隙度结果(a)—(e) 孔隙度预测结果； (f)—(j) 孔隙度预测结果通过解码网络D生成的地震数据； (k)—(o) 生成地震数据(f)—(j)与观测地震数据(图2b)之间的残差. (a)—(e)对应的正则化参数分别为0、0.1、1.5、0.37、0.37.此外，(a)—(d)对应的迭代次数都设置为300次，(e)对应的迭代次数设置为1000次.Fig.5 Porosity results via the semi-supervised porosity prediction network(a)—(e) Porosity prediction results； (f)—(j) Generated seismic data by inputting the predicted porosity of (a)—(e) into the decoder D，respectively； (k)—(o) Data residual between the generated seismic data of (f)—(j) and the observed seismic data of Fig.2b，respectively. The regularization parameter corresponding to (a)—(e) is set to 0，0.1，1.5，0.37 and 0.37， respectively. In addition, the iterations for (a)—(d) are set to 300，and for (e) are set to 1000.

在直观对比图4和图5中的预测孔隙度剖面差异基础上，进一步对比有监督与半监督两种网络在盲井位置预测的单道孔隙度曲线的差异.图6每个子图中的黑线、蓝线和红线分别表示从真实孔隙度模型(图2a)、有监督网络预测的孔隙度剖面(图4a)和半监督网络预测的孔隙度剖面(图5d)中提取出的来自同一CDP位置的孔隙度曲线.图6a—c对应的CDP位置分别是16、499和688.在经过楔形体及下伏地层(图6a的椭圆内)时，孔隙度曲线出现了急剧增大到固定值后又急剧减小的现象.两种网络都没有很好地适应这种纵向跨度比较大的孔隙度突变，但是半监督网络相比于有监督网络预测效果略好.在经过断层及背斜等(图6b、c椭圆内)区域时，半监督网络比有监督网络能更好地适应局部孔隙度的复杂变化，说明前者更适用于实际非均质性强的复杂储层孔隙度预测.从图6可以看出，半监督网络相比于有监督网络能提高不同位置的单道孔隙度预测精度，最终在孔隙度预测剖面上表现为减小数据驱动的逐道反演引起的横向不连续性.

图6 有监督与半监督网络盲井测试结果对比基于第16道(a)、第499道(b)、第668道(c)地震记录预测的孔隙度曲线.Fig.6 Comparisons of blind well test via supervised network and semi-supervised networkPredicted porosity curves based on the 16th (a)，499th (b)， and 668th (c) seismic record. Black，blue and red lines represent true supervised network inverted，and semi-supervised network inverted porosity curves，respectively.

实际地震数据在野外采集时往往由于人为或环境因素的干扰而受到噪声污染，为进一步说明半监督网络相比于有监督网络具有更好的抗噪性.采用与图4和图5d分别对应的无噪数据训练得到的有监督网络和半监督网络直接推广到含噪合成地震数据(图2c，即模拟的实际数据).此时，半监督网络预测孔隙度剖面(图7b)比有监督网络预测孔隙度剖面(图7a)更为干净和连续，说明前者比后者更适用于低信噪比地震资料.表1统计了图4—图7中所有测试预测孔隙度剖面与真实孔隙度剖面之间的相对误差和相关系数，及半监督网络生成地震数据与对应的输入地震数据之间的相关系数.

图7 有监督与半监督网络含噪测试结果对比(迭代300次)(a) 有监督网络基于含噪地震数据(图2c)得到的孔隙度预测结果； (b) 半监督网络基于含噪地震数据(图2c)得到的孔隙度预测结果.Fig.7 Comparisons of noise-tolerance test results via supervised network and semi-supervised network (300 iterations)(a) Predicted porosity result by inputting noisy seismic data of Fig.2c into the supervised network； (b) Predicted porosity result by inputting noisy seismic data of Fig.2c into the semi-supervised network.

表1 不同方法预测孔隙度结果对比Table 1 Comparisons of porosity prediction results using different methods

Marmousi模型数据测试结果表明：(1)有监督或半监督孔隙度预测方法相比于传统的联合地震反演与岩石物理模型的反演方法能避免弹性参数到物性参数的二次转换引起误差累积和放大，实现直接从地震数据预测出孔隙度，且预测结果更为准确.(2)相比于有监督网络，半监督网络通过非井位置的地震数据匹配间接保证井间孔隙度的预测精度.井间地震数据匹配程度一定程度上可以间接反映孔隙度预测质量.当非井位置的孔隙度预测结果也较好而不出现较大偏差时，其通过解码网络才能合成与观测地震数据接近的生成地震数据.半监督网络通过学习反映局部地层变化的大量地震波形信息而增强了刻画薄层孔隙度细节及薄厚层孔隙度突变的能力，提高了孔隙度预测的精度和横向连续性.(3)在使用的有标签数据一致的情况下，半监督网络与有监督网络的预测效果存在一定关联.即当公式(3)中的正则化参数越小时，二者预测效果越接近.半监督网络需要引入合适的地震数据匹配约束数据驱动反演过程，以提高孔隙度预测精度.

2.2 实际数据应用

最后，采用来自中国东部某地区的叠后地震数据对本文提出的半监督孔隙度预测方法进行检验.由于工区直井缺少实际解释的孔隙度曲线，本文首先采用工区内已解释的斜井孔隙度曲线与不同的孔隙度敏感参数进行交会分析.如图8所示，在统计工区111口斜井在储层段(1300～1400 m)的孔隙度与速度关系发现，使用线性回归即可较好地从速度曲线拟合出孔隙度曲线，拟合结果与已知孔隙度值的相关系数较高，为0.8385.在明确孔隙度曲线生成方法的基础上，下面针对某连井剖面(图9a)开展储层孔隙度预测.该连井剖面有141道，每道地震记录有111个时间采样点，道间距为12.5 m，时间采样间隔为2 ms.连井剖面的时间范围为1.1～1.32 s，所在储层埋深较浅，自上而下依次为曲流河沉积和辫状河沉积，主要的沉积微相包括心滩坝、辫状河道、河道边缘和河间滩等.储层物性较好，平均孔隙度为30%左右，各小层的平均渗透率在(600～1100)×10-3μm2左右.由于河流频繁改道及多条大断层的发育使得储层非均质性强，岩性纵横向变化大，呈现出多个砂泥岩薄互层纵向交替叠置.在地震剖面上表现为整体较破碎，地震同相轴不连续且变化快.

图8 储层段孔隙度与速度的交会图Fig.8 Cross-plot of porosity versus velocity in the reservoir intervals

如图9a所示，该连井地震剖面从左到右依次经过6口井(命名为A1—A6).6口井经过井震标定合成的地震数据与对应的井旁地震道匹配程度高，平均相关系数为0.79.利用图8的岩石物理关系拟合出的6条孔隙度曲线，选取前5条孔隙度曲线及其对应的地震数据训练有监督网络和半监督网络.训练的迭代次数都设置为3000次，正则化参数设置为1，其他设置与合成数据例子相同.训练好的有监督和半监督孔隙度预测模型推广到整个地震剖面(图9a)，预测孔隙度剖面分别为图9b、c.自然电位(SP)曲线在该工区能很好地刻画河道砂体，通过插入井位置的SP曲线(图9b、c中的虚线)可以看出，整体上有监督网络和半监督网络预测的孔隙度具有一定的相似性，但是后者在井旁预测的孔隙度与SP曲线变化趋势更为吻合，更好地表征出不同河道砂体孔隙度的差异.有监督或半监督网络预测的孔隙度模型在剖面上部(1.1～1.25 s)以相对低孔地层包裹高孔地层为主，符合曲流河出现的“泥包砂”现象；在剖面下部(1.25～1.32 s)相对高孔地层包裹低孔地层为主，符合辫状河呈现的“砂包泥”现象.但是，半监督方法比有监督方法预测出的孔隙度结果纵横向分辨率更高，连续性更好.此外，半监督孔隙度预测结果(图9c)对应的生成地震数据(图9d)与实际地震数据(图9a)接近，二者残差小(图9e)，表明预测的孔隙度剖面是满足地震数据匹配的.

进一步局部放大图9a—c，比较有监督方法与半监督方法描述辫状河储层段孔隙度的差异.图10a—c分别为储层段的过井地震剖面，有监督和半监督方法预测的储层段孔隙度.其中，图10a—c中2条红线表示储层段上下界面位置的两个目标层位，图10b—c中的黑框部分为井位置的真实孔隙度，井间为预测孔隙度.相比于有监督方法，半监督方法预测的孔隙度横向变化在盲井位置(A6)附近更为合理.从图10b—c中的黑框可以看到，A5井位置的孔隙度值比其余5口井的要大.而有监督方法建立的地震到孔隙度的非线性映射关系很大程度依赖于测井标签，导致A5井参与有监督网络训练后，A4井与A5井之间孔隙度预测结果(图10b的白框)出现成片的高孔特征，这与辫状河横向砂体变化快的地质认识不符.半监督方法比有监督方法提高了地震数据在孔隙度预测中的参与程度，缓解了特殊测井标签(A5)对网络建模的影响.通过地震波形匹配约束降低了数据驱动逆过程的多解性，间接保证不同空间位置孔隙度的差异性.因此，在A4井与A5井之间的半监督孔隙度结果(图10c的白框)更符合地质认识，表明半监督方法在储层段刻画不同尺度河道砂体的孔隙度的能力更强.最后，对比一口训练井(A1)和一口测试盲井(A6)的单道孔隙度预测结果可以看到，半监督孔隙度预测方法(图11中的红线)和有监督孔隙度预测方法(图11中的蓝线)都能较好地实现训练井位置的孔隙度曲线拟合，但是前者的盲井拟合效果更符合实际孔隙度曲线(图11b中的黑线)的整体趋势，刻画出更为准确的小层孔隙度变化细节.在盲井位置，半监督方法和有监督方法预测的孔隙度结果与真实孔隙度的相关系数分别为为0.8151和0.6949.

图11 基于有监督与半监督网络的训练井与测试井预测孔隙度结果对比(a) 训练井(A1)位置的孔隙度预测结果； (b) 测试井(A6)位置的孔隙度预测结果. (a)、(b)中的黑线、蓝线和红线分别代表真实、有监督网络预测、半监督网络预测的孔隙度曲线.Fig.11 Comparisons of predicted porosity curves between supervised network and semi-supervised network at the training well and test well locations(a) Predicted porosity curves at the training well (A1)； (b) Predicted porosity curves at the test well (A6). Black，blue and red lines in (a) and (b) represent true，supervised network inverted，and semi-supervised network inverted porosity curves，respectively.

实际数据测试结果表明：半监督方法比有监督方法更适用于表征实际非均质性储层的孔隙度.它通过构建地震测井双约束目标函数降低了测井标签对孔隙度建模的决定性作用，地震数据波形匹配控制缓解了特殊测井数据对建立的孔隙度预测网络的负面影响，并间接保证横向孔隙度预测精度，实现更充分地利用地震测井多尺度信息估计更为准确的孔隙度模型.

3 结论

本文通过模拟传统的地震正反演与岩石物理建模过程，提出了一种基于双向门控递归单元神经网络的半监督井震联合储层横向孔隙度预测方法.通过充分利用少量的(伪)测井数据及对应的井旁地震道和大量的非井旁地震道建立地震数据到孔隙度、孔隙度到地震数据的闭环映射，实现地震数据直接预测孔隙度.合成数据和实际数据预测结果表明，相比于仅使用成对的地震数据和测井标签建立的有监督孔隙度预测方法，半监督孔隙度预测方法具有更好的泛化性和鲁棒性.其在保证井位置预测精度的同时通过非井旁地震数据匹配进一步提高井间孔隙度预测精度，进而改善预测孔隙度剖面的横向连续性和准确性.通过学习不同位置不同信噪比的实际地震数据，半监督方法更好地适应地震波形变化，并降低对测井标签的过度依赖，缓解特殊测井数据引起的孔隙度建模误差.因此，半监督孔隙度预测方法更适用于实际非均质性储层，实现合理地刻画不同辫状河河道砂体形态及其内部的孔隙度差异.

地震数据匹配项的正则化参数直接影响半监督孔隙度预测方法的应用效果，如何自动优选出合理的正则化参数有待进一步研究.此外，半监督孔隙度预测方法本身受到地震测井数据质量的固有限制，提高地震数据的信噪比和分辨率，制作数量充足、质量可靠、更贴近实际的(相控)孔隙度曲线也是下一步需要研究的工作.