基于长短期记忆神经网络补全测井曲线和混合优化XGBoost的岩性识别

潘少伟, 王朝阳, 张 允, 蔡文斌

(1.西安石油大学计算机学院,陕西西安 710065; 2.中国石化石油勘探开发研究院,北京 100083;3.西安石油大学石油工程学院,陕西西安 710065)

准确识别岩性对油田一系列地质工作的开展具有重要意义。目前通过测井曲线识别岩性的方法主要有3种[1-10]:第一种是利用测井曲线和其他岩石物理特征建立交会图等传统的岩性识方法[1];第二种是基于机器学习方法,如朴素贝叶斯(naive Bayes,NB)、决策树(decision tree)、支持向量机(support vector machine,SVM)和随机森林(random forest,RF)等所构建的岩性识别方法[3-9];第三种是基于深度学习(deep learning,DL)方法,如卷积神经网络(convolutional neural networks,CNN)进行岩性识别的方法[10]。利用测井曲线和其他岩石物理特征建立交会图等传统的岩性识别方法可依赖测井工作者的知识、经验和相关地质信息较为准确地进行岩性识别,但此种方法的主观性较强,且耗费时间较长。朴素贝叶斯[3-4]、决策树[5-6]和支持向量机[7-8]的应用在一定程度上提高了岩性识别的效率与准确率,但这些方法本身均存在不足。首先,朴素贝叶斯假定各属性之间是相互独立的,但测井曲线中很多属性是彼此关联的,朴素贝叶斯在解决不同岩性识别问题中带有先天的不足。其次,决策树容易陷入过拟合、忽略数据集中不同属性之间的相关性;支持向量机在处理多分类问题时的训练时间较长,参数选择较困难。随机森林由于对小样本数据的分类效果较差,因此在岩性识别领域的应用受到了限制[9]。此外,通过卷积神经网络进行岩性识别也获得了应用[10],但卷积神经网络模型的训练时间一般较长,且训练集与测试集的构建难度较大。研究区块位于中国南方某油田,地质状况复杂,并且由于某些客观原因,该区块内一些取芯井的测井曲线中存在较多缺失值。笔者首先通过长短期记忆(long short-term memory,LSTM)神经网络建立时间序列预测模型重新生成缺失的测井曲线值,引入XGBoost算法构建不同岩性的识别模型,并利用粒子群优化(particle swarm optimization,PSO)算法和网格搜索(grid search,GS)法对XGBoost的主要超参数进行寻优,最后应用粒子群优化算法和网格搜索法寻优的XGBoost算法完成对测试数据集中不同岩性的识别。

1 原理与方法

1.1 LSTM

循环神经网络(recurrent neural networks,RNN)在计算时间序列中距离较远节点之间的联系时,会涉及雅可比矩阵的多次相乘,导致出现梯度消失或梯度爆炸现象。因此一些学者对RNN进行了改进,包括增加有漏单元、设计门限RNN等。Hochreiter等[11]提出长短期记忆神经网络,即LSTM神经网络。LSTM神经网络具有门(gate)和记忆单元(cell)两种结构。门结构通过抛弃无用信息、增添有用信息的方式保证信息有选择性地通过,这样使得记忆单元状态不断地更新,从而达到控制和保护记忆单元状态的目的。一个LSTM记忆单元内有3种门结构,分别是遗忘门(forget gate)、输入门(input gate)和输出门(output gate)。遗忘门的作用是对记忆单元状态中某些信息进行选择性遗忘,保留有用信息;输入门的作用是控制新信息的输入并将其选择性地记录到记忆单元状态中;输出门的作用是确定记忆单元只输出与当前时间同步的有关信息。这三种门结构在记忆单元中进行矩阵乘法和非线性求和等运算,确保其记忆在迭代过程中不会衰减。

LSTM记忆单元的典型结构如图1所示。由图1可知,在t时刻记忆单元的输入包括前一时刻的隐含层状态变量ht-1、记忆单元的状态变量ct-1和当前时刻的输入信息xt。之后依次通过遗忘门ft、输入门it和输出门ot,记忆单元获得t时刻的隐含层状态变量ht和状态变量ct。最终ht会传入输出层生成LSTM在t时刻的计算结果y,同时ht和ct一起进入到t+1时刻进行新一轮的计算。

图1 LSTM记忆单元典型结构

为了计算输出值,在t时刻LSTM隐含层的计算过程如下:

ft=σ(Ufxt+Wfht-1+bf),

(1)

it=σ(Uixt+Wiht-1+bi),

(2)

(3)

(4)

ot=σ(Uoxt+Woht-1+bo),

(5)

ht=ottanh(ct).

(6)

其中

式中,ct∈Rm为记忆单元状态变量;xt∈Rp为输入变量;ht∈Rm为输出变量;U(·)∈Rm*p和W(·)∈Rm*p为权重矩阵;b(·)∈Rm为偏置向量,在训练过程中将被不断地优化;ft、it和ot分别代表遗忘门、输入门和输出门的值;σ为sigmod激活函数,取值范围为(0, 1)。tanh(x)的值域是(0, 1)。

1.2 粒子群优化算法

粒子群优化算法[12]是一种随机搜索算法,可以较大的概率收敛并得到全局范围内的最优解。在动态的、多目标的优化问题中,相比传统优化算法,粒子群优化算法具有更快的搜索速度以及避免陷入局部最优的能力,可节省大量的寻优时间,并收敛至最接近的最优解。本文中采用粒子群优化算法对XGBoost算法中的浮点型超参数进行寻优。

1.3 网格搜索法

网格搜索法是一种穷举遍历算法,它首先把所有可能的参数组合在空间中划分成若干个网格,之后遍历网格中的所有交点,对每个参数集合应用交叉验证来计算误差,最后得到误差最小的为全局最优解[13]。网格搜索法可以从诸多参数中获得最优解,但是效率较低,因此本文中仅使用网格搜索法确定XGBoost算法整型超参数的最优值。

1.4 XGBoost算法

梯度提升树(gradient boosting decison tree,GBDT)算法[14]的原理是基于所有先前树,使用梯度下降法生成新树,使得目标函数靠近最小的方向。XGBoost[15]被称为极端梯度提升树,是一种树集成模型,可以用于解决分类和回归问题。XGBoost应用于分类问题时,先不断添加新的回归树,再通过新生成的CART树来拟合前一次模型的残差。训练出的完整模型包括K棵树,每棵树对应的结果之和作为最终的预测值。

式中,xi为第i个数据采样点的特征值集合;R为所有CART树的函数空间;fk(xi)为第i个数据采样点在第k棵树中落在叶子节点处的权重。

正则项的具体表达式为

式中,γ和λ为惩罚系数;T为叶子的个数;w为叶子节点的权重。

XGBoost的目标函数又可改写为

式中,fs(xi)为第s轮训练的新子模型;C为常数项。XGBoost通过不断地寻找fs来优化目标函数,目标函数进行二阶泰勒展开并且移除高阶无穷小项变为

Ωfs+C.

其中

目标函数又可做进一步的化简,具体过程为

(7)

1.5 粒子群和网格搜索寻优的XGBoost算法

采用粒子群优化算法和网格搜索法对XGBoost算法中的7个超参数进行寻优,这7个超参数分别是树的叶子节点进一步分裂所需的最小损失减少量(gamma,记作pg)、构建树时对特征的列采样率(colsample_bytree,记作pcb)、训练样本的采样率(subsample,记作psub)、促使算法更快收敛的参数(scale_pos_weight,记作psca)、每棵树每次节点分裂列采样的比例(colsample_bylevel,记作pcl)、树的最大深度(max_depth,记作pmd)和每棵树权重改变的最大步长(max_delta_step,记作pmds)。把粒子群优化算法和网格搜索法混合寻优的XGBoost算法简记为PSO-GS-XGBoost,其执行的基本流程如下。

(1)利用粒子群优化算法对XGBoost算法中的pg、pcb、psub、psca和pcl5个浮点型超参数进行寻优,确定pmd和pmds2个整型超参数的取值范围,具体步骤为:①初始化粒子的种群规模,局部学习因子c1、全局学习因子c2、惯性权重ω、最大迭代次数和适应度函数func()等;②应用初始化的参数进行迭代训练,计算每个粒子的适应度;③ 若迭代到一定次数,整个粒子群收敛,则此时pg、pcb、psub、psca和pcl5个浮点型超参数的值即为它们的最优值;如果整个粒子群不收敛,且没有达到最大迭代次数,则返回②继续进行训练。

(2)通过网格搜索法,结合粒子群优化算法寻优获得的pg、pcb、psub、psca和pcl的最优值,进一步确定pmd和pmds的最优值。

(3)由XGBoost算法中5个浮点型超参数的最优值和(2)中所获得的pmd和pmds2个整型超参数的最优值,构建PSO-GS-XGBoost岩性识别模型。

1.6 评价指标

为评价LSTM神经网络对缺省测井参数的预测效果,选取均方根误差(记作prmse)、平均绝对误差(记作pmae)和平均绝对百分比误差(记作pmape)作为评价指标。计算公式为

2 基于PSO-GS-XGBoost的岩性识别

2.1 研究区概况

研究区位于中国南方某油田,面积约为2.5 km2,自1995年正式投入开发,主要产油层段为古近系阜宁组一段(简称阜一段)和阜宁组二段(简称阜二段)。其中阜一段为三角洲前缘亚相沉积,分为前缘席状砂、河口坝、水下分流河道和水下分流间湾微相;以灰色粉砂岩、棕色含油粉砂岩和灰色泥岩不等厚互层为特征,三角洲相砂体是主要储集层之一;自然伽马负异常比较明显,电阻率曲线值普遍较高,自然电位曲线表现为钟形、漏斗形、箱形、直线形及其复合形态。阜二段为滨浅湖亚相沉积,发育砂坝、生物滩、鲕粒滩、灰质滩和滨浅湖泥微相,鲕粒滩和生物滩的含油性较好,生物滩挖掘潜力较大;为一套湖相沉积地层,发育灰色、深灰色粉砂岩、生物灰岩、鲕粒灰岩和泥岩;自然伽马出现中幅低值异常,底部电阻率曲线呈尖刀状,称“山字形”段,自然电位曲线多呈直线形、钟形、箱形及其复合形态。研究区内现有H1～H8共8口取芯井。

2.2 基于LSTM神经网络的测井曲线补全

在测井工作中测井仪器发生故障或其他客观原因导致测井曲线中的一种或多种参数在某一深度段没有采集到数据或采集到异常数据。前人在解决此类问题时,或采用某种物理模型直接反演得到缺失或异常的测井参数[16-17],或采用交会图和多元回归方法生成缺失或异常的测井参数[18-19]。这两种方法在早期的测井曲线补全工作中发挥了积极作用,但均存在不足:物理模型反演得到测井曲线的方法简化了真实的地层状况;交会图和多元回归方法则忽略了测井曲线中不同参数之间复杂的映射关系。

测井曲线表征的地层按照时间先后顺序沉积而成,并且测井曲线之间的采样间隔较小。可见,测井曲线的采样点之间存在一定的时间相关性,可将其看作时间序列数据。LSTM神经网络在当前时间序列数据预测中得到了广泛应用[20-23]。本文中将LSTM神经网络应用于取芯井缺失值处或异常值处的测井参数补全中。以H1井测井曲线中的中感应电阻率(RILM)为例,研究基于LSTM神经网络的测井曲线补全过程。

H1井测井曲线中共有3 032条数据记录(表1),前2 801条数据记录非常完整,但其后231条数据记录的中感应电阻率中出现了空值。综合地质研究表明,取芯井不同测井曲线参数之间存在较强的非线性映射关系。将H1井测井曲线中前2 801条数据记录作为LSTM神经网络的数据集,并按照7∶3的比例将其划分为训练数据集和测试数据集。把H1井的自然电位、自然伽马、深感应电阻率、浅侧向电阻率和自然电位作为LSTM神经网络的输入,中感应电阻率作为输出,时间步长设置为30,利用训练数据集构建补全中感应电阻率的LSTM神经网络预测模型,通过测试数据集对该预测模型进行测试,测试合格后将其应用于H1井中感应电阻率缺失值处新数据的生成中。

表1 H1井测井曲线参数

生成中感应电阻率缺失值的LSTM神经网络预测模型的训练过程如图2所示。由图2可看出,当训练轮次达到20后,LSTM神经网络预测模型的损失值基本趋于稳定,训练轮次达到100时,训练完成,此时的损失值介于0和0.01之间。利用测试数据集对LSTM神经网络预测模型进行测试,测试结果使用prmse、pmae、pmape3个评价指标表征。再从测试数据集中截取300个数据点,将中感应电阻率的真实值与LSTM神经网络预测模型的预测值进行对比,所得结果如图3所示。

图2 LSTM神经网络损失函数随训练轮次变化

生成中感应电阻率缺失值的LSTM神经网络预测模型在测试数据集上产生的prmse为0.024 5,pmae为0.020 6,pmape为4.004 0,均极小。由图3可看出,中感应电阻率的真实值曲线与LSTM神经网络预测模型的预测值曲线除在极值处存在一定差异外,在其他数据点处均拟合得很好。这说明本文中构建的中感应电阻率缺失值的LSTM神经网络预测模型是准确的、合理的。将该模型应用到H1井测井曲线后231条数据记录中即可预测获得中感应电阻率缺失的231个数据值。

图3 中感应电阻率真实值与预测值对比

2.3 数据集构建

鉴于研究区面积较小,内部井网较密,且不同井内发育的地层非常相似,因此从8口取芯井中随机抽取5口取芯井H1、H3、H5、H6和H8来构建数据集。根据阜一段和阜二段的综合地质特征和测井数据特点,采用岩心分析孔隙度数据与测井声波时差数据计算机自动相关对比的方法,将岩心深度归位到测井深度上,使岩心数据与测井数据相匹配。之后把阜一段和阜二段的全部岩性划分为灰岩、泥质砂岩、砂岩和泥岩4类,分别用1、2、3和4表示,结合岩心数据完成对H1、H3、H5、H6和H8井的测井曲线岩性标签的确定,由H1、H3、H5和H6井共获得11 495条数据记录(表2)。表2中第8列为标签值,不同标签值代表了不同的岩性。

表2 由H1、H3、H5、H6和H8井获得的测井曲线数据集

2.4 数据集预处理

在应用多维度的测井曲线建立不同岩性识别模型时,由于不同测井参数的量纲不一致,且同一种测井参数数值间的差距也较大,需要对测井曲线数据集进行预处理,以提高预测准确率,减少训练时间。使用标准差归一化公式来完成预处理,公式为

2.5 基于粒子群和网格搜索XGBoost算法的超参数寻优

把预处理之后的测井曲线数据集按7∶3的比例划分成训练数据集和测试数据集。基于训练数据集,利用粒子群优化算法对XGBoost模型中pg、pcb、psub、psca和pcl5个浮点型超参数和pmd和pmds2个整型超参数进行寻优。经过多次调试后,设置粒子的种群规模为50,局部学习因子c1为0.2,全局学习因子c2为0.5,惯性权重ω为0.5,最大迭代次数为50,适应度函数func()为待优化XGBoost算法在训练数据集上进行岩性识别的准确率。最终得到浮点型超参数的寻优结果如表3所示。

表3 XGBoost浮点型超参数的寻优结果

同时获得整型超参数pmd的最优值范围为[3,10],pmds的最优值范围为[0,5]。保持5个浮点型超参数的最优值不变,利用网格搜索法对整型超参数pmd与pmds进行寻优,最终寻优得到pmd为6,pmds为0。

2.6 PSO-GS-XGBoost岩性识别模型构建与测试

确定pg、pcb、psub、psca、pcl、pmd、pmds的取值分别为5.413、0.467、0.743、3.711、0.707、6和0时,即成功建立了PSO-GS-XGBoost岩性识别模型。

利用测试数据集对该PSO-GS-XGBoost岩性识别模型进行测试,所得结果如图4所示。图4中共有6幅子图,依次对应PSO-GS-XGBoost、随机森林、未经优化的XGBoost算法、决策树、支持向量机和朴素贝叶斯在同一测试数据集上的测试结果。

由图4可看出:在6种岩性识别方法中,本文中提出的PSO-GS-XGBoost算法对灰岩、泥质砂岩、砂岩的识别准确率均最高,分别为75.00%、99.40%、99.90%,朴素贝叶斯对泥岩的识别准确率最高,为99.20%;而支持向量机对灰岩和泥岩的识别准确率最低,分别为0和86.20%,朴素贝叶斯对泥质砂岩和砂岩的识别准确率最低,分别为89.70%和80.10%;随机森林、未经优化的XGBoost算法和决策树对灰岩、泥质砂岩、砂岩和泥岩4种不同岩性的识别准确率介于支持向量机和朴素贝叶斯与PSO-GS-XGBoost算法之间。在4种不同岩性中,6种不同方法对灰岩的识别准确率最高为75.00%,最低为0;对泥质砂岩,砂岩和泥岩识别的最高、最低准确率分别为99.40%、89.70%,99.90%、80.10%和99.20%、86.20%。

图4 6种岩性识别方法对测试数据集中不同岩性识别结果的统计直方图

综合上述分析可以发现:在对灰岩、泥质砂岩、砂岩和泥岩4种不同岩性的识别中,PSO-GS-XGBoost算法的表现最好,而朴素贝叶斯和支持向量机的表现最差,这说明基于树的机器学习方法在经过超参数寻优之后,在不同岩性识别中要优于未进行超参数寻优的同类机器学习方法和朴素贝叶斯、支持向量机、决策树与随机森林这些经典机器学习方法;6种岩性识别方法对灰岩的识别准确率均不高,这是由于灰岩的训练数据集较少所致,但相对而言,PSO-GS-XGBoost算法对灰岩的识别准确率最高,这充分体现出其在岩性识别方面的先进性与优越性。

2.7 PSO-GS-XGBoost岩性识别模型的验证

从H2、H4和H7井中随机抽出H7井对PSO-GS-XGBoost岩性识别模型进行验证。同样把H7井测井曲线中的自然电位、自然伽马、深感应电阻率、中感应电阻率、浅侧向电阻率和自然电位作为PSO-GS-XGBoost岩性识别模型的输入,计算输出相应深度下岩性的标签值(1、2、3和4)。把这些标签值投影到H7井的测井曲线图中,所得的结果如图5所示(图中包含H7井的真实岩性剖面和PSO-GS-XGBoost识别的岩性剖面)。

由图5可看出,本文中PSO-GS-XGBoost算法对不同岩性的识别结果与真实的岩性类别基本吻合。但由于不同的分类标准,本文中PSO-GS-XGBoost算法仅能够对砂岩、泥岩和灰岩实现比较准确地识别,而对同属于砂岩类别的灰质砂岩和粉砂岩则不能进行识别;并且由于训练数据的缺失,在对灰岩识别时出现了误判,将1 653.0 m处的一小部分灰岩识别为砂岩。

3 结束语

为克服已有岩性识别方法中存在的不足,提出一种利用粒子群优化算法和网格搜索法对XGBoost算法进行混合寻优的岩性识别方法。与随机森林、未经优化的XGBoost算法、决策树、支持向量机和朴素贝叶斯相比,在不同岩性识别中PSO-GS-XGBoost算法的表现最为优异。由于不同的分类标准,提出的PSO-GS-XGBoost岩性识别方法在对实际岩性识别时会出现一定的误差,且当训练数据不足时,该方法对已有岩性识别的准确率也会降低。此外,在进行测井数据采集时,受制于油田的实际情况未能采集到更多的测井曲线数据,如中子密度测井曲线等。在后续的研究工作中,将采集更多的测井曲线数据,进一步细化基于PSO-GS-XGBoost的岩性识别模型,提高识别准确率。