沉积微相是油气勘探开发工作中至关重要的研究内容之一, 尤其对于明确剩余油富集区产量预测起着重要作用。传统沉积微相划分中, 常以研究前人成果、区域沉积背景为基础, 通过古生物学和沉积学等理论结合关键井岩心相标志分析, 确定研究区内沉积类型; 结合多种测井曲线形态特征, 综合分析岩石厚度、粒度等特征, 划分研究区沉积微相^[1-2]。多种资料交叉分析的过程复杂且繁琐, 很难找出沉积相与各种测井数据之间的定性映射关系。

随着计算机技术的飞速发展, 机器学习在沉积微相上的运用主要方法有: BP神经网络^[3]、支持向量机(SVM)^[4]、模糊聚类神经网络^[5]、KNN算法^[6]和贝叶斯算法^[7], 然而, 这些机器学习算法仍存在自身局限性, 在很多情况下不能满足沉积微相识别要求。BP神经网络虽然实现了输入到输出的映射功能, 具有较强的非线性映射能力, 但是BP神经网络收敛速度慢, 容易陷入局部极小化, 导致网络训练失败, 同时无法观察学习过程, 输出结果难以理解; SVM算法在分类效果上表现很好, 但是, 对大规模训练样本存在分类难的问题, 参数调优过程比较繁琐。模糊聚类分析对噪声和异常值很敏感, 但其结果不稳定; KNN算法计算量过大, 样本不平衡时, 对少数量样本预测准确率低; 贝叶斯算法对多分类问题很有效, 对待测样本预测, 过程简单且速度快, 但是它有分布独立的假设前提, 而现实中预测问题很难完全独立。此外, 这些算法主要依附于特定的地质环境和沉积背景, 构造对自然伽马曲线形态类和物性类参数特征, 而这些定量特征不具有泛化能力。

深度学习具有挖掘数据非线性关系的优势, 目前尚无有效解决沉积微相识别的模型。递归神经网络能够随着时间推移进行顺序处理, 2020年, SONG等^[8]将长短期记忆网络运用于测井曲线形态识别, 仅采用单向序列对钟形、漏斗形、卵形和箱形4类测井相进行识别。经本文实验证明, 用于时间序列预测的时间卷积网络(TCN)同样可用于沉积微相识别, 但是, 长短期记忆网络和时间卷积网络都只能对测井曲线序列在深度上进行单向建模, 当输入的原始信号采样率较低, 利用两者对原始信号进行处理均可以取得较好的效果。但是, 在测井曲线序列中常用的信号大多是高采样率信号, 能够有效反映井沉积特征的信号长度通常大于1000, 甚至更长。而长短期记忆网络和时间卷积网络直接处理含有大量噪声的高采样率原始信号时, 会造成模型参数过大, 使模型难以训练, 容易产生过拟合现象。另外, 普通的长短期记忆网络和时间卷积网络模型只能记忆当前深度之前的曲线变化特征。为了准确反映整口井测井序列的沉积特征, 亟需研究更准确、泛化能力更强的有效识别沉积微相深度学习网络模型。

本文提出了一种DMC-BiLSTM深度学习的沉积微相识别模型。利用DMC特征构造, (D代表地质趋势特征, M代表中值滤波特征, C代表聚类特征), 目的是为了更好提取原始自然伽马曲线特征, 消除噪声对原始信号鲁棒性和泛化能力的影响, 为双向长短期记忆(BiLSTM)网络提供易捕获的特征, 从而提高算法的性能。利用双向长短期记忆网络提取当前测井信号的前后深度序列变化特征, 该变化特征能够更好地反映沉积过程中能量的变化, 预测沉积微相类型。

1 DMC-BiLSTM模型

图 1给出了基于DMC-BiLSTM深度学习的沉积微相智能识别流程, 具体实现过程如下。

1) 数据预处理。主要是构造三大类特征, 包括地质趋势特征, 中值滤波特征和聚类特征, 简称DMC特征构造。生成新的特征集合。

2) 数据划分。将数据划分为训练集和测试集。

3) 模型训练。搭建4层双向长短期记忆网络网络模型, 训练DMC-BiLSTM分类模型。

4) 模型验证。使用测试集验证DMC-BiLSTM分类模型的分类性能。

1.1 DMC特征构造

测井数据在采集的过程中, 由于受到外界干扰或者仪器自身原因, 导致有效信号中夹杂着各种噪声。本文提出的DMC特征构造非常适合测井曲线的信号处理, 在进行沉积微相识别之前, 需要对测井曲线序列进行预处理, 构造新的特征矩阵。

1.1.1 趋势分量

STL算法^[9]为时序分解中的一种常见算法, 将某时刻数据分解为趋势分量、周期分量和余项, 根据对自然伽马曲线的初步分析, 加法模型满足分解要求, 即:

$ x_{t}=T_{t}+S_{t}+R_{t} \quad t=1, \cdots, N $

(1)

式中: x_t为某深度下的曲线数据; T_t为趋势分量; S_t为周期分量; R_t为余项; t为某深度。

1.1.2 中值滤波

中值滤波基本原理为: 通过设定窗口值, 某点的输出信号, 就可以用这个点的窗口范围内所有信号的统计中值来代替。一维中值滤波可定义为:

$ y_{t}=\operatorname{Med}\left(x_{t-N}, \cdots, x_{t}, \cdots, x_{t+N}\right) $

(2)

式中: x_t为深度t时的曲线数据; y_t为深度t时滤波后的数据; t为某深度。

1.1.3 聚类特征

Kmeans是一种无监督学习聚类算法^[10], 具体实现步骤为:

1) 将自然伽马曲线特征数据聚为K类, 选择K个自然伽马数据点作为初始中心点;

2) 按照距离初始中心点最小原则, 将所有曲线数据分到各中心点所在类中;

3) 每类中有数据若干, 计算K个类所有曲线样本数据均值, 作为第二次迭代的K个中心点;

4) 根据该中心点重复第2和第3步, 直到收敛(中心点不再改变或达到指定迭代次数)聚类过程结束。

欧氏距离是便于理解也最常用的一种距离公式。通过欧氏距离可以更好体现微小深度段内自然伽马曲线序列之间的相似程度, 能充分体现自然伽马曲线空间相关性。其定义为:

$ d=\sqrt{\sum\limits_{k=1}^{n}\left(\boldsymbol{x}_{1 k}-\boldsymbol{x}_{2 k}\right)^{2}} $

(3)

式中: x_1k, x_2k为两个n维自然伽马曲线特征向量。

1.2 双向长短期记忆网络原理 1.2.1 基本原理

自然伽马数据可以看作是随深度变化的时间序列, 其形态分类依赖于当前输入和之前状态。假设给定自然伽马序列x={x₁, x₂, …, x_t, …, x_N}, 双向长短期记忆网络^[11]可以同时对每个序列进行正向和反向建模, 因为每个标记编码都包含来自过去和未来的上下文信息, 因此可以更丰富地表示测井序列数据的长期依赖性。双向长短期记忆网络每层由单个长短期记忆网络单元正反向传播构成(图 2)。

长短期记忆网络可以解决循环神经网络(RNN)无法处理长距离的依赖问题。原始循环神经网络的隐藏层只有一个状态, 即h, 它对短期的输入非常敏感。增加一个状态c, 用来保存长期的状态, 称为单元状态。再用3个门, 即遗忘门(f_t)、输入门(i_t)和输出门(o_t), 来控制c。结合构成了长短期记忆网络模型的核心部分“记忆块”(图 3)。总的公式为:

$ \left\{\begin{array}{l} f_{t}=\sigma\left(W_{f} x_{t}+U_{f} h_{t-1}+b_{f}\right) \\ i_{t}=\sigma\left(W_{i} x_{t}+U_{i} h_{t-1}+b_{i}\right) \\ o_{t}=\sigma\left(W_{o} x_{t}+U_{o} h_{t-1}+b_{o}\right) \\ \hat{c}_{t}=\tanh \left(W_{c} x_{t}+U_{c} h_{t-1}+b_{c}\right) \\ c_{t}=f_{t} \otimes c_{t-1}+i_{t} \otimes \hat{c}_{t} \\ h_{t}=o_{t} \otimes \tan h\left(c_{t}\right) \end{array}\right. $

(4)

式中: σ和tanh为两种激活函数; c_t表示t时刻的存储单元; $ \hat{c}_{t}$表示t时刻的候选存储单元; W={W_f; W_i; W_o; W_c; U_f; U_i; U_o; U_C}, 表示权重系数, 由反向传播进行更新; b={b_f; b_i; b_o; b_c}, 表示偏置系数; x_t表示t时刻序列值; h_t表示隐含层输出。

1.2.2 模型配置和训练

本文搭建4层的双向长短期记忆网络单元对自然伽马曲线序列学习, 通过全连接层将学到的“分布式特征表示”映射到样本标记空间, 最后将输出向量输入到包含5个神经单元Softmax层, 执行测井相分类任务。为了避免模型过拟合, 在网络每一层添加系数为0.2的Dropout单元, 通过舍弃一些连接, 增加模型的泛化能力; 采用Adam优化算法, 对梯度进行更新, 加快模型的收敛; 使用交叉熵损失函数作为训练损失, 有:

$ \left\{\begin{array}{l} y=S\left(W_{c} \boldsymbol{x}+b_{c}\right) \\ \theta=-\frac{1}{n} \sum\limits_{\boldsymbol{x}}[y \ln a+(1-y) \ln (1-a)] \end{array}\right. $

(5)

式中: x为输入样本, 本文为20维的自然伽马特征矩阵; y为期望的沉积微相类型输出; S为Softmax函数; W_c为输出节点c的权重系数; b_c为输出节点c的偏置系数; θ为交叉熵损失函数; a为神经元实际输出; n为总类别, 本文为5种, 即分流间湾、前缘席状砂、水下分流河道、河口坝和河道侧缘。

网络层数越多非线性拟合能力越强, 能够识别曲线特征的复杂度越高, 但是, 过深的网络结构容易造成模型过拟合。为了使模型具有最优的准确率和泛化能力, 通过对比不同网络层数在测试集上准确率和损失值的表现, 选择最优的网络层数。如图 4和图 5所示, 网络输入都是20维特征矩阵, 每个模型迭代300轮, 运用早停法, 即模型在20轮中损失值不再下降, 便停止迭代。例如, 图 5中4层双向长短期记忆网络模型在迭代70轮左右后, 损失值在接下来20轮中都趋于平稳, 因此100轮左右后便停止训练, 避免了模型过拟合现象的发生。对比图 4中的曲线, 可看出网络层数从1层增加到4层, 模型在测试集上的准确率不断上升, 在网络层数为4层时, 模型的拟合能力最佳, 准确率最高, 而损失值也收敛得最快, 最终的误差值最小。当网络模型超过4层时, 模型的准确率下降, 并且曲线跌宕幅度变大, 过拟合风险增加。因此, 本文采用了4层的双向长短期记忆网络层数。最终, DMC-BiLSTM沉积微相识别模型如图 6所示, 具体参数如表 1。

2 实验结果与分析 2.1 测井数据来源

本文实验用测井数据来自加州硅谷人工智能公司CrowdAnalytix提供的伽马测井相类型预测(Gamma Log Facies Type Prediction)数据。自然伽马测井特征的变化可以作为粒度的表征, 使其在层序地层学分析中发挥重要作用^[12]。自然伽马测井响应的突变被解释为与不整合和层序边界相关的尖锐的岩性断裂。因此, 测井曲线的基本形态常用来解释沉积旋回和沉积背景^[13]。表 2是三角洲前缘亚相测井曲线类型。

2.2 划分训练集和测试集

共有100口井数据, 每口井的数据共1100条。选取其中80口井的数据作为训练集, 余下20口井的数据作为测试集, 用来验证本文方法对沉积微相的预测能力, 表 3是具体划分的训练集和测试集。

2.3 数据预处理

图 7中给出本文方法的预处理流程, 具体实现过程如下。

1) 对原始自然伽马测井曲线, 采用STL算法, 并使用不同奇数频率窗口进行分解, 去除周期分量和余项, 得到多个地质趋势分量特征; 其次采用与趋势分解相同频率窗口, 对原始曲线进行中值滤波, 得到多个中值滤波特征。

2) 将原始自然伽马曲线和第1步中得到的地质趋势分量和中值滤波特征, 作为Kmeans聚类的输入的特征集合, 进行无监督学习, 聚类时K值的选择要大于沉积微相类别数, 此处选择K=6, 得到具有时空相关性聚类特征。

3) 原始自然伽马曲线、地质趋势分量、中值滤波特征和聚类特征组成新的特征组合, 对新特征集合进行归一化处理。

实验均在Intel(R)Core(TM)i5-8300CPU@2.3GHz、8GB RAM设备上进行。图 8以9号井为例, 纵坐标表示曲线数值大小, 横坐标表示深度(需要特别说明的是, 原始数据中无深度数据, 每口井有1100条数据, 文中以每口井的第1个点到最后1个点的长度作为理论深度), 其中分解、滤波频率窗口值均为19, 对自然伽马曲线进行趋势分量特征、中值滤波特征、Kmeans聚类特征构造。能够看出新的特征集合不仅能够去除曲线中的高频噪声, 保留地质趋势特征及有效的边缘形态特征, 还能充分体现曲线序列的时空相关性。在原始数据中为了使模型拟合能力更强, 取数值范围为3~20所有奇数滤波窗口(奇数滤波相比偶数滤波有中心点, 对测井曲线整体形态以及曲线边缘特征更加敏感, 可以更有效地提取边缘信息)进行自然伽马曲线特征构造。最终新特征与原始自然伽马曲线数据, 共同生成20维特征矩阵。

2.4 实验结果

表 3反映出数据中分流间湾类型明显多于其余4个微相类型, 存在数据不均衡现象, 因此, 采用单一评价指标并不科学。本文引入多个评价指标, 包括每一类的混淆矩阵(如表 4所示)、精确率、召回率、F1值和整体的准确率^[14], 利用这些指标对预测结果做全方位评价, 综合反映模型分类效果。

1) 混淆矩阵。

2) 准确率。

模型预测正确的结果占总结果的比重为:

$ A=\frac{T_{P}+T_{N}}{T_{P}+T_{N}+F_{P}+F_{N}} $

(6)

3) 精确率。

模型预测为正例的结果中的真正实例所占比重为:

$ P=\frac{T_{P}}{T_{P}+F_{P}} $

(7)

4) 召回率。

模型预测中样本的正例有多少被预测正确。

$ R=\frac{T_{P}}{T_{P}+F_{N}} $

(8)

5) F1值。

精确率和召回率的调和平均值:

$ F_{1}=2 * \frac{P * R}{P+R} $

(9)

公式(6)~公式(9)中: A为准确率; P为精确率; R为召回率; F₁为F1值; T_P为真正例; F_P为假正例; T_N为真负例; F_N为假负例。

从混淆矩阵及各类评价指标(表 5)中可以看出, 本文方法在测试集上分流间湾、前缘席状砂、水下分流河道、河口坝、河道侧缘分类精确率分别达到了96%, 81%, 93%, 87%和88%。其中, 对呈低幅微齿状的分流间湾划分效果最佳, 可能与分流间湾主要以泥岩沉积为主, 表现为高自然伽马值有关, 易识别; 其次, 是对箱形水下分流河道的识别, 精确率、召回率以及F1值都达到90%以上; 在其它曲线形态识别上, 各类指标值略低, 可能因为三者在自然伽马曲线形态上区分不大, 只是在沉积时水流进积和退积的区别, 但在预测的精确率也都达到80%以上。总而言之, 模型在测试集上的整体预测能力表现优异, 平均精确率、平均召回率和平均F1值都在85%以上, 准确率也达到了91.69%。

为了进一步评估所提出的DMC-BiLSTM模型, 还使用其它两种适用于序列预测的深度学习模型(长短期记忆网络和时间卷积网络)与本文方法在每个类别的识别精度、整体准确率和运行时长等方面进行对比。需要特别说明的是, 在本文实验中, 3种神经网络除了基础算法本质上的不同外, 其余参数都尽量保持了一致, 神经网络隐含层的激活函数均设置为ReLU函数, 输出层的激活函数设置为Softmax函数。在模型的训练过程中, 都使用了小批量梯度下降和Adam优化器, 对梯度进行更新, 且运用了早停法对3个模型均训练300轮。DMC-BiLSTM模型采用了4层的双向长短期记忆网络单元, 长短期记忆网络模型采用了4层网络单元, 时间卷积网络模型采用了4层的残差单元。模型的输入均为20维自然伽马曲线特征, 长度为1100的测井曲线序列。

最终结果如表 6所示。在运行的实效性方面, 此处以秒(s)作为运行时长的单位。本文方法的双向长短期记忆网络单元是由正向和反向传播的长短期记忆网络单元构成, 模型的运行时长大概是普通长短期记忆网络单元的两倍左右; 而时间卷积网络可以对测井序列进行并行运算, 降低了运行时长。虽然本文方法在运行时长上相比其它两种深度学习算法不具有优势, 但明显缩短了人工划分的工作时长。本文方法无论在每一类识别的精确率上, 还是对沉积微相识别的整体准确率上效果都是最佳的, 最终准确率达到91.69%。与其它两种模型相比, 高出长短期记忆网络单元模型7%左右, 高出时间卷积网络模型12%左右。证明了本文方法在沉积微相分类准确度方面具有明显的优势。

图 9中给出了3种模型对97号井沉积微相识别的可视化结果。前4列是4种曲线特征, 自然伽马列表示原始自然伽马曲线, 地质趋势和中值滤波列都是以滤波窗口19为例对应的两种曲线特征, 聚类特征列就是前3种特征集合, 经过Kmeans无监督学习聚类后的曲线特征; 后4个柱状图, 分别代表着人工标定的沉积微相, 以及DMC-BiLSTM、长短期记忆网络和时间卷积网络3种模型对97号井沉积微相识别的结果。对比发现, 本文方法的沉积微相识别结果, 与地质专家人工标定的结果吻合率更高, 能够很好地识别出不同曲线形态类别。不足之处在于, 当水动力发生改变, 同时产生的能量变化微弱, 那么曲线的顶底接触只发生微小变化, 但沉积微相过渡成另外一种类型, 模型在识别这种情况下的不同沉积微相之间的顶底接触往往会出现误差(图 9中红色箭头所示)。当水动力能量变化较大时, 本文方法能够很好地识别不同沉积微相之间的顶底界限。而其它两种方法识别结果是紊乱的, 不仅不能够很好地体现出整口井的沉积规律, 而且在识别不同类别之间的界限上表现较差。

3 结论

本文提出的基于DMC-BiLSTM的沉积微相智能识别方法, 构建了3种特征, 即地质趋势特征、中值滤波特征和聚类特征。地质趋势特征能够很好地反映一定深度内曲线的变化趋势; 中值滤波特征既能去除测井曲线中的高频噪声, 又能保留有效的边缘信息; 聚类特征则能更好地反映测井曲线的时空相关性, 区分不同曲线形态之间的边界。将本文方法与仿真得到的时间卷积网络模型和长短期记忆网络模型进行比较, 本文方法的分类准确率为91.69%, 高于时间卷积网络模型的79.21%和长短期记忆网络模型的83.99%。本文方法可以实现测井曲线的端到端学习, 对分流间湾、前缘席状砂、分流河道、河口坝和河道侧缘等各类沉积微相的识别精确率分别为96%, 81%, 93%, 87%和88%。实验结果表明, 本文方法有助于提取测井曲线序列的隐藏特征, 准确地划分不同沉积微相的分类界限。该方法能有效识别沉积微相, 具有较强的泛化能力和鲁棒性, 为沉积微相的识别提供了一种新途径。