基于长短期记忆神经网络的随钻地层倾角解释方法

孙歧峰，李娜，段友祥，李洪强，唐海全

（1.中国石油大学（华东）计算机科学与技术学院，山东青岛 266580；2.中国石化胜利石油工程有限公司钻井工艺研究院，山东东营 257000；3.中国石化胜利石油工程有限公司测控技术研究院，山东东营 257000）

0 引言

随钻测井与地层评价技术通过实时获取钻头附近的轨迹参数、地质参数和工程参数来实时评价所钻遇地层的情况，进而精确控制钻具运动轨迹以命中最佳地质目标[1]，是高精度地质导向钻井的关键支撑技术。随钻测井技术在分辨地层变化及确定岩性[2-5]、评价泥砂/砂泥岩含量[6]、预测高压地层[7-8]等方面有明显效果，能够提高钻井工程的效益[9]。由于实时数据传输的信息非常有限且解释难度大[10]，将人工智能与随钻测井解释相结合[11-13]以提高实时数据处理的准确性和效率有重要的意义。

针对随钻伽马测井解释过程中面临的问题，结合长短期记忆神经网络（LSTM）及小波变换等[14-16]，本文提出一种实时地层倾角解释方法：首先提取井的层位特征并计算动态阈值，初步判定地层边界位置，然后训练分类器模型，筛选真实的地层边界，最后计算地层相对倾角及方位角。对方法进行详细阐述并开展现场应用。

1 随钻方位伽马数据解释流程

对随钻仪器实时传输上来的钻井液脉冲信号进行解码得到井下各种传感器的测量数据，进行坐标转换、时深转换及校正，然后对数据进行插值得到方位伽马的成像结果。在此基础上，进行随钻方位伽马数据的解释，流程如下。

①地层边界初步划分。首先，采用Lowess（局部加权回归散点平滑法）滤波方法对方位伽马曲线进行过滤，消除与地层变化无关的随机噪声[17]引起的测井曲线中细微尖锐的变化，保留地层变化时的趋势性特征。然后，基于小波变换模极大值计算测井曲线上的极大值点，用Lipschitz指数（Lip）表征测井曲线上极大值点的奇异性以验证各极大值点的真伪，实现测井曲线层位特征的初步提取。

②动态阈值计算。根据伽马图像上不同层序的特征与范围，以熵和类间差为衡量标准计算所有层序的局部阈值。

③地层识别。依据各层序的局部阈值，计算各层序的轮廓点集合。将轮廓点集合作为输入数据，通过LSTM分类模型筛选真假地层边界轮廓。

④倾角计算。用非线性最小二乘法将各轮廓点集合拟合到正弦曲线，得到正弦曲线的振幅、初相位、角速度等关键信息。由倾角计算公式计算得到地层的相对倾角，根据井斜参数得到地层真倾角及方位角。

2 关键技术

2.1 基于小波变换的层位特征提取与动态阈值计算

小波变换作为一种处理信号的重要手段，在层序划分[15]、图像增强[18]等方面有广泛的应用。在信号奇异值计算中，通常使用卷积运算的连续小波变换。

2.1.1 基于小波变换的奇异检测原理

设曲线函数为 f（t），Wf（s,t）为 f（t）在对应母波函数 ψs（t）上的分解。s为尺度因子，在实际应用中，通常取值为 s=2j（j∈Z）[19]，且尺度的选择必须与具体信号结合，通常选定j=4或j=5[20]，Wf（s,t）即为f（t）的卷积型二进制小波变换。

假设 t0是曲线函数 f（t）在尺度 s0下的局部极值点，则有：

如果t0某一邻域内的任意点t都满足如（2）式所示的条件，则称 t0为小波变换的模极大值点，称Wf（s0,t0）为小波变换的模极大值。

基于小波变换模极大值计算曲线函数f（t）的奇异值，通过Lip指数表征f（t）在某点上的奇异性。若曲线函数f（t）在某点处连续可微，导数有界且不连续，则该点的Lip指数为1，该点无奇异性；若f（t）在该点处不连续但有界，则该点的Lip指数为0，该点有奇异性。f（t）在某点处的尖锐程度随Lip指数的减小而增大，通过计算Lip指数可以得知f（t）上某点的奇异程度，辅助校验曲线极大值点的突变特性。

2.1.2 层位特征提取与阈值计算

基于小波变换的层位特征提取与动态阈值计算方法流程为：首先对随钻伽马测井曲线进行滤波处理，然后依据小波变换模极大值及 Lip指数对伽马曲线进行层序划分，最后根据动态划分结果计算局部阈值。

2.1.2.1 测井曲线噪声影响及滤波处理

对伽马曲线的平滑要突出不同层岩性的强烈变化，削弱同层内岩性的细微变化，避免产生由伽马测井仪器的放射性测量统计涨落、数据传输时的外界干扰、同层岩性伽马波动等因素[17]导致的曲线小直径波峰现象，使供层位提取使用的测井曲线突变位置更加清晰。

首先采用Lowess滤波方法对方位伽马曲线进行过滤，其基本思路参见文献[21]。以随钻伽马平均曲线为例，各滤波方法对测井曲线平滑处理的结果如图 1所示。图 1表明，与直线型形态学滤波、半圆型形态学滤波相比，Lowess滤波不会丢失重要的曲线突变特征；与Savitzky-Golay滤波相比，Lowess滤波保留的重要突变特征的位置更吻合原曲线形态。

图1 滤波效果对比

2.1.2.2 测井曲线层序划分

基于Lowess滤波后的平均伽马曲线，采用小波变换模极大值计算平均伽马曲线的奇异点，使用Lip指数检测曲线上点的奇异性。平均伽马曲线上的突变点易被检测为奇异点，而突变的位置往往代表可能存在的地层变化，如图 2所示。因此，可以基于小波变换模极大值定位地层可能存在变化的深度，作为计算局部阈值和地层边界轮廓点的依据。

图2 曲线与地层变化关系示意图

结合测井图像特点及地层分布的规律，对于分层点 a0，若在[a0-C，a0+C]内无其他分层点，则 a0的局部阈值计算范围为[a0-C，a0+C]；若在[a0-C，a0]内有其他分层点（设距离最近的分层点为a1），在[a0，a0+C]内无其他分层点，则 a0的局部阈值计算范围为[(a0-a1)ρ，a0+C]；若在[a0-C，a0]内无其他分层点，在[a0，a0+C]内有其他分层点（设距离最近的分层点为a2），则 a0的局部阈值计算范围为[a0-C，(a2-a0)ρ]；若在[a0-C，a0]和[a0，a0+C]内均有其他分层点（设两个区间内距离最近的分层点分别为 a1和 a2），则 a0的局部阈值计算范围为[(a0-a1)ρ，(a2-a0)ρ]。其中，C=δ(D+2H)/λ，为常量，决定分层点a0在查找上下有无其他分层点时的最大范围。δ和ρ分别为与上下无分层点、上下有分层点时的划分范围相关的参数，根据相对倾角在伽马成像图中展开的正弦曲线振幅规律，本文设δ和ρ的值分别为20和0.85。H为伽马图像的成像深度，依据文献[22]研究成果计算后作为常数值处理，本文设为0.035 m。

根据分层点及其局部阈值计算范围，将伽马成像数据划分为不同的层序，然后计算每一层序的局部阈值。

2.1.2.3 局部阈值计算

本文使用熵和类间差作为指标衡量阈值[13]，该指标综合考虑了类内聚和类间离散的特征，局部阈值的计算方法不再赘述。

2.2 基于LSTM的地层识别模型

通过层位特征提取计算的层序受曲线形态、滤波或算法效果等因素的影响势必会导致提取到非真实地层，因此需要再对各层序进行筛选。LSTM模型的记忆门结构使得该模型较适合解决具有序列特征的问题[23]，在多种场景下有广泛的应用[24-25]，而地层的轮廓数据具有序列的特征，因此采用LSTM模型设计地层识别分类器。

2.2.1 地层识别分类器设计

井筒与地层边界相交，反映在随钻伽马成像图上为一条具有正弦特征的曲线，对轮廓点集合的筛选可以抽象为LSTM网络的正弦特征二分类问题。根据计算的轮廓点集合有限点的特征及LSTM设计原则，构建地层识别分类器模型，如图3所示。

图3 基于LSTM的地层识别分类器模型

训练集由两部分组成：第 1部分为人工测井解释并标注的样本数据；第 2部分为计算机根据地层边界特征生成的样本数据。规定正样本是具有正弦特征的边界点集合，负样本是不具有正弦特征的点集合。训练集总共包含300组正样本和300组负样本，每组的数据点数量为 N（N∈[50，100]），N的范围与伽马成像预处理时的方位伽马插值数量相关。从后续的分类结果可见，样本数量设置为600组即可达到预期精度。设计正样本训练数据时，应结合正弦曲线表达式及地层边缘的特征，在此基础上产生的训练数据集示例如图4所示。模型最终输出结果为1，0标签，分别表示数据点集合正弦特征的有、无。

图4 建模产生的训练数据集示例

LSTM 模型隐藏层神经元的选择对训练结果有较大的影响，通常依据经验公式确定神经元个数[26]：

采用交叉熵函数作为损失函数计算分类器结果与样本标签的误差，通过误差反向传播进行训练，在交叉熵平稳变化后，用 Adam梯度优化算法继续调整，与其他梯度优化器相比，Adam优化器在LSTM模型的实际训练中效果更优[27]。

2.2.2 模型验证与分析

结合（3）式设定LSTM模型中不同的隐藏层神经元节点数，对比不同神经元节点数的交叉熵值，选取最小误差的节点数作为地层分类器的参数。本文不同神经元节点数的交叉熵值如图5所示。

图5 不同神经元节点数的交叉熵

由图5可见，节点数在12附近时训练效果最优，在4，6，8处效果较差。因此，本文选用12作为隐藏层的神经元节点数，在该节点数下的地层识别分类器模型训练结果如图6所示。该模型训练在epoch=80附近收敛，epoch=120附近时平均损失值下降至1.2×10-4左右，模型测试集的平均损失值约为2.3×10-4，说明该模型能够正确地对有正弦特征的轮廓点数据和无正弦特征的轮廓点数据进行分类。

图6 模型收敛曲线

在实际地层识别中，需要先计算各层序的局部阈值，得到各层序内的轮廓点集合：假设某层序的阈值为th1，遍历该层序中的所有数据点，若f(xi,yi)=th1，则标记(xi,yi)为轮廓点；若 f(xi,yi)-th1与 f(xi+1,yi)-th1结果符号相反，则其中结果较小的数据点也将标记为轮廓点。依此类推，计算所有层序的轮廓点集合并保存至边界集合，边界集合即由各轮廓点集合组成。运用地层识别分类器模型对边界集合中的各轮廓点集合进行正弦性检查，去除无正弦特征的轮廓点集合，使边界集合能反映真实的地层边界。

通过上述方法对某深度上的18组轮廓点数据进行筛选，结果如表1所示。结果表明，训练的基于LSTM的地层识别分类器具备筛选真、假地层的能力。

表1 轮廓点集合识别结果

2.3 倾角计算

通过拟合边界集合中经过正弦性检查后的所有轮廓点集合，可获得地层边界轮廓对应的正弦信息。拟合方法为：结合非线性最小二乘法，通过重复迭代寻找函数的局部最小值，从而使目标函数取得最优解。根据拟合后得到的正弦信息（见（4）式）及文献[17]中的方法计算地层相对倾角和方位角，公式分别如（5）式、（6）式所示，具体推导过程不再赘述。

3 方法验证与应用

为验证本文方法的准确性和可行性，从方位伽马数据解释和随钻实时数据处理两方面进行应用及分析。

3.1 方位伽马数据解释

计算解释的数据选用中国某油田X井井场8道方位伽马数据（2 100～2 380 m井段），该井方位伽马数据前期处理完善、伽马成像结果清晰。通过对真实的井数据进行相关地层解释，验证本文方法的准确性，包括算法的地层识别和倾角解释能力。

本文方法的分层结果、拟合结果、倾角识别结果及与全局阈值法[13]的对比如图 7所示。可以看出，本文方法有较好的分层及轮廓拟合效果，可以准确地描述地层的变化、刻画地层边界轮廓曲线；与全局阈值法相比，本文方法能解释出更多的地层倾角。局部井段的地层倾角本文方法解释结果与人工解释结果对比如表 2所示，可见本文方法解释结果与人工解释结果相当。全部井段长度为280 m，测点22 400个，人工解释层数为50，本文方法解释层数为46，本文方法在较长井段、较多测点数的情况下，层数识别率为92%，倾角解释准确度达到96.6%，证明本文方法准确性较好。

表2 采用本文方法及人工方法对真实方位伽马数据进行解释后得到的地层倾角解释结果对比

图7 本文算法与全局阈值法计算结果对比

3.2 随钻实时数据处理

Y井是中国东部油田某采油厂完善采油井网中的一口井，是为挖潜小层高部位剩余油、进一步提高储量动用程度、改善开发效果而钻探的一口开发水平井。地质导向服务采用了中国自主研制的近钻头地质导向系统+随钻方位伽马+一体化随钻测量解释系统组成的新型地质导向系统。地质导向技术服务井段从A靶前开始，从井深2 184 m钻进到井深3 075 m完钻，累计进尺891 m，井下工作82 h。系统不仅实现常规MWD（随钻测量）井斜、方位、工具面、温度等轨迹姿态测量参数上传，同时上传的近钻头井斜、两道伽马（上、下）测量曲线也为实时地质导向提供了及时准确的决策依据。

使用该井数据验证本文方法在实时钻井生产中插值成像、地层边界解释效果。在该场景下本文方法的处理结果如图 8所示。可以看出，本文方法基于数据插值的方法得到了实时伽马数据的清晰成像；基于地层划分及局部阈值计算的方法得到了符合实际情形的地层划分；基于LSTM分类器模型得到了准确的地层边界，特别是图8中井深2 440～2 450 m位置，本文方法通过分类器模型筛查出了多条无正弦特征的轮廓点集合，并结合快速正弦拟合的方法，将轮廓点拟合到了正确的位置。由表 3可见，本文方法在真实井实时数据中仍能够保持可接受的误差，与人工解释的结果大致相同，该段耗时约为0.27 s，比人工解释效率高。通过对图8与表3所示结果的分析，证明了本文方法在实际生产中的效果。

表3 采用本文方法及人工方法对随钻实时数据进行解释后得到的地层倾角解释结果对比

图8 井场随钻实时数据的解释结果

以上计算结果与分析表明，本文方法在自动地层识别方面有较高的准确率且计算速度较快，具备处理实时随钻测井数据的能力。实际应用中，利用近钻头井斜测量确保轨迹准确入靶并在靶体中延伸，借助近钻头上下方位伽马及电阻率测量数据进行处理分析，保证了目的层钻遇率。

4 结论

相较于全局阈值法和固定阈值法，基于小波变换的动态阈值法使地层界面解释的准确率显著提升，且在地层岩性变化稳定时减少了无用计算，并适用于地质构造变化较快的情况。

运用LSTM设计了基于随钻伽马数据特征的地层边界分类器模型，建立了地层边界特征数据集，在训练过程中调整神经元节点数及学习率等参数避免导致欠/过拟合的问题。分析对比发现，该分类器模型能准确识别出具有正弦特征的轮廓，从而筛选出真实的地层边界。

以人工智能为基础的随钻测井解释方法能够为钻井工程技术人员的决策提供辅助，为实时地质导向应用提供技术支持。

符号注释：

a0，a1，a2——分层点，即层位提取步骤中的小波变换模极大值点位置；A——振幅；D——钻孔直径，m；f(·)——曲线函数；H——伽马图像的成像深度，m；i——曲线上点的序号；j——与尺度因子相关的系数；K——常数，K∈[0,10]；m，n——输入层和输出层的节点数；M——神经元个数；N——每组样本中数据点数量；s——尺度因子；th1——某层序的局部阈值；t——曲线 f(t)上的某点；t0——曲线函数 f(t)在尺度 s0下的局部极值点；Wf(s,t)——f(t)在对应母波函数ψs(t)上的分解；x——直角坐标系x轴坐标，m；xi，yi——曲线上某点的坐标，m；y0——偏距，m；Z——整数集；α，β——地层相对倾角和方位角，（°）；δ——与上下无分层点时的划分范围相关的参数；λ——测井仪器测点的间距，m；ρ——与上下有分层点时的划分范围相关的参数；φ——初相；ψs(t)——母波函数；ω——角速度。