基于一维卷积神经网络的地层智能识别方法

曹茂俊,崔欣锋

(东北石油大学 计算机与信息技术学院,黑龙江 大庆 163318)

0 引 言

如今,随着社会经济的快速发展,石油的消费不断增加。为了满足社会发展需要,需要进一步勘探开发油气资源,而地层智能划分可以为油气资源的持续开发节省人力和物力。

测井曲线包含了随井深变化的储层性质,其直接反映了地下的地质特征。因此,测井曲线的地层划分对于确定地层位置具有重要意义[1]。长期以来,地层划分的人工解释方法以测井解释人员的主观判断为主。人工解释方法需要人工观察测井曲线图形形状的变化以寻找对应地层的半幅点位置,这些方法比较依赖专业领域知识和专家经验[2],而且由于不同专家的经验和熟练程度不同,也会导致地层划分结果的巨大差异,不仅需要大量的人力和物力,而且分层误差也相对较高。随着计算机技术的发展和划分任务的日益繁重,地层划分逐渐由定性识别转向定量识别[3],数理统计、活度函数、小波变换、模糊识别、人工神经网络等方法[4]已经被应用到利用测井曲线进行地层划分的任务中。然而,这些方法大都只考虑个别曲线进行划分,而且大都只考虑点与点之间的关系,而忽略了测井曲线随深度变化上下连续的特点,即局部特征,而深度学习方法为该问题提供了非常具有前景的解决方案。地层划分问题是根据测井曲线的形态识别出不同深度所属的地层,因此可转化为多元序列的分类问题[5]。卷积神经网络(Convolutional Neural Network,CNN)在图像识别领域已经是相对成熟的技术了,由于其强大的图像分类能力也被应用于测井曲线地层的划分[6],但考虑到测井曲线一般为具有一定规律但不规则的形态,对于地层划分来说,输出数据为某一个分层类别,无法实现测井曲线的地层划分。测井曲线序列是一种二维的多元序列,不能直接采用传统的CNN进行地层划分,一维卷积神经网络(1DCNN)在处理序列数据方面效果比较突出[7]。因此,该文选择1DCNN作为地层划分的基础模型,通过对测井曲线进行IMK特征工程处理,即对原始测井曲线分别采用INPEFA技术[8]、中值滤波处理以及K-Means聚类[9]处理,更好地提取了原始曲线的特征,为1DCNN网络提供了更加容易识别的特征,从而增加了模型的地层划分效率。该文结合已有的网络模型及测井曲线的特点,搭建了基于IMK-1DCNN的地层划分网络模型,从而达到对地层划分的智能化方法研究的目的。

1 IMK-1DCNN模型

如图1所示,基于IMK-1DCNN的测井曲线地层识别模型流程大体分为4个部分:

图1 IMK-1DCNN地层识别流程

(1)特征构造。利用INPEFA技术、中值滤波、K-means算法分别构造INPEFA曲线特征、中值滤波特征、聚类特征,简称IMK特征构造,与原有的测井曲线共同构成新的特征集合。

(2)数据集划分。将新的特征集合划分为训练集、验证集和测试集。

(3)模型训练。搭建一维卷积神经网络模型,使用训练集训练IMK-1DCNN模型,使用验证集进行误差计算和权值更新。

(4)模型验证。最终在测试集上验证IMK-1DCNN的分层效率。

1.1 IMK特征构造

图2展示了IMK特征构造的流程,为了处理测井曲线在采集过程中产生的噪声,以及更好地提取原始测井曲线的地层特征,在模型训练之前利用特征程对原始曲线进行了处理,采用IMK特征构造方法构造了新的特征,与原始测井曲线特征共同构成了新的特征矩阵。

图2 IMK特征构造流程

上述特征构造过程中的INPEFA技术是在最大熵频谱分析法(Maximum Entropy Spectral Analysis,MESA)基础上形成的预测误差滤波分析(Prediction Error Filter Analysis,PEFA)通过进一步发展形成合成预测误差滤波分析(Integrated Prediction Error Filter Analysis,INPEFA)技术[10]。其原理是通过最大熵频谱变换,可以由已知点的值在熵最大的原则下推算出下一个点的数值,这个值称为最大熵频谱分析估计值(MESA)。预测误差滤波分析(PEFA)是通过计算每一个深度点的MESA预测值和对应的测井曲线实际值而得到数据差值,如下式所示:

PEFA=RV-MESA

(1)

其中,RV代表实际值,MESA代表最大熵谱分析估计值,得到的较大正值或负值为不连续点,可能为地层不连续点或旋回界面。对PEFA曲线做一个特定积分处理,就得到了合成预测误差滤波分析曲线(INPEFA),其本质是反映泥岩含量的多少且基于测井曲线分析,因此它能够显示通常在原始测井曲线中显示不出来的趋势和模式以及反映沉积的旋回性变化,通常被用来辅助测井解释人员识别隐藏的地层界面[11]。该文运用INPEFA技术利用数据集中的GR曲线构造INPEFA曲线,并将其作为模型的特征之一来构造样本集。

1.2 数据集划分

经过IMK特征构造之后得到了新的特征集合,数据集的合理划分能够很好地加速神经网络的集成,这里将其按照8∶1∶1的比例划分为训练集、验证集和测试集,其中训练集用来训练模型,验证集用来调整并寻找最佳参数,最后在测试集上进行测试以验证模型的泛化性能及准确率。

1.3 模型训练及验证

在准备好数据集后,需要搭建模型进行训练,该文选用了善于处理序列数据的一维卷积神经网络作为模型的基础架构,处理数据集使之满足模型的输入并进行训练,训练时需要设置网络的批大小、训练轮次、早停策略、Dropout比例等参数以加速网络的训练。在训练完成之后,要在测试集上进行测试来验证模型的识别效果,并通过绘图来对比模型的识别结果与地质专家的标记结果。

2 基于IMK-1DCNN的地层识别方法

2.1 1DCNN网络基本原理

(2)

式中,hi表示输入序列中(i)th个感受野的输出[13]。可以利用1DCNN对测井曲线序列切片做一维卷积运算,提取该深度区间的层位特征信息,从而实现划分地层类型的目的。

2.2 模型总体框架

设计的1DCNN网络架构如图 3所示。网络中各层的含义解释如下。

(1)输入层:将输入的与地层密切相关的测井曲线通过IMK特征构造新的特征矩阵,从而输入到网络中进行特征提取并学习。

(2)卷积层:卷积层使用的是适合处理序列数据的一维卷积(Conv1d),一维卷积常被用以在一维输入信号上进行邻域滤波。为了减少过拟合现象,还在每个卷积层后面引入了BN(Batch Normalization)层[14],它可以加快网络的训练和收敛的速度,控制梯度爆炸且防止梯度消失。

(3)池化层:池化层的本质是一个下采样,因为数据经过卷积之后维度越来越高,而且特征图没有多大改变,在连续多个卷积之后,会产生很大的参数量,不仅会大大地增加网络训练的难度,还容易造成过拟合的现象,所以在卷积层之后放一个池化层对数据进行压缩,降低维度,减少参数量。该文使用了两种池化层,一种是最大池化层(Max Pooling),可以减少卷积层参数误差造成估计均值的偏移,另一种是全局平均池化层(Global Average Pooling),可以减小特征提取时邻域大小受限导致估计值方差增大造成的误差[15]。

(4)Dropout层:为了减少参数量,加快模型的训练速度,缓解过拟合现象,该文引入了Dropout层,通过随机选择其中的一些神经元并将其临时隐藏,然后再进行本次训练和优化。在下一次迭代中,继续随机隐藏一些神经元,如此直至训练结束。由于是随机丢弃,故而每一个mini-batch都在训练不同的网络。设置Dropout层的隐藏比例为0.3(随机将30%的神经元临时隐藏)。

(5)Dense层:即全连接层,在模型中起到“分类器”的作用,它会将学到的“分布式特征表示”映射到样本标记空间。通过将经过网络训练的输出向量输入到包含8个神经元的Dense层,输出最终的地层分类。

2.3 数据归一化

由于不同测井曲线的量纲不同,为了消除量纲的影响,需要对数据进行归一化。如果对未归一化的数据直接进行建模,可能会导致模型对数值大的变量学习过多,而对数值小的变量训练不够充分,往往模型效果不佳。常用的数据归一化的方法有最大最小值归一化、小数定标法、定量特征二值化及均值方差标准化等。

选用的归一化方法是最小值-最大值(Min-Max)方法,顾名思义,该方法是利用数据列中的最大值和最小值进行归一化处理,处理完之后可将数据范围缩放到[0,1]范围内。计算方式为数据列中的数据与该列的最小值作差,再除以其极差。设一条测井曲线为X=(x1,x2,…,xn),则归一化后的数据x'可由式(3)计算得出:

(3)

其中,x'代表数据值x归一化的结果,max(X)和min(X)分别代表测井曲线数据列中的最大值和最小值。

2.4 滑动窗口采样

深度学习需要大量样本进行训练,而单纯地将测井曲线构造成独立的测井图像不能满足深度学习网络需求。在传统的深度学习网络中,往往会采用数据增强的方式来扩充数据集,使得网络能够学习到足够的特征,比如旋转、平移、缩放、随机遮挡、水平翻转和噪声扰动等[16]。但是由于测井曲线具有时序性,其垂向的形态变化往往代表着某种沉积环境,从而导致对测井曲线矩阵利用传统的数据增强方法会破坏其本身所代表的地层和沉积环境的含义。

为了增加样本量,该文提出一种使用滑动窗口对测井曲线进行采样的方法。设滑动窗口的长度为L,步长为S,数据生成的方式如图4所示,将特征矩阵划分为多个长度为L,宽度为特征数量的切片,每个切片对应的地层标签yi由每个切片中的出现次数最多的标签决定,同时为了避免切片中的地层标签差异过大,步长S不宜设置过大,然后利用样条插值[17]对每个切片进行处理以保证生成的输入矩阵长度大于特征数目。这里使用样条插值是为了在扩充切片长度的同时保证训练曲线的平滑,最大程度地保留原始曲线的地层趋势。

图4 利用滑动窗口生成测井曲线切片示意图

3 实际应用与结果分析

3.1 数据来源

实验使用的数据来自四川盆地川中地区,共包含20口井的测井数据(数据分布见图5)。四川盆地位于四川省龙门山断裂处以东,四面环山,是中新生代以来发展起来的大型沉积盆地。该文研究的主要地层为茅口组和栖霞组,其中茅口组顶部与上二叠统龙潭组呈平行不整合接触,上覆于栖霞组之上。茅口组厚度为119 m～208 m,平均厚237 m,岩性特征为灰色—深灰色亮晶生物碎屑灰岩、泥晶生物碎屑灰岩,泥质灰岩夹硅质结核[18]。栖霞组厚度为23 m～313 m,平均125 m,岩性特征为:下部为深灰色、灰黑色泥质沥青质灰岩,底部常见眼球状微晶灰岩夹炭质沥青质页岩;上部为浅灰色、灰褐色微—亮晶生屑灰岩。

3.2 样本集构造

所用的样本集共20口井,共92 284条数据。将其中18口井作为训练集,剩余2口井作为测试集,验证模型对地层的识别性能。在构造样本集之前要对数据进行预处理。在实际测井中由于测井数据受多因素影响,可能导致数值缺失或出现异常值。因此,对于原始数据,首先要进行缺失值和异常值的处理,该文选择了用众数填充空值及异常值,实验证明该方法简单有效。

结合地质领域知识,以及1.1节的IMK特征构造方法,构造新的特征。最终选取自然伽马(GR)、自然电位(SP)、声波时差(AC)和INPEFA曲线、中值滤波曲线及K-Means聚类曲线共6种测井曲线作为训练特征。其中,GR曲线主要测量地层的放射性,SP曲线主要测量在地层电化学作用下产生的电位,AC曲线反映了地层声波速度,常被用来进行地层划分。另外,在构造中值滤波曲线特征时为了更好地提取测井曲线的边缘信息,使用了20以内的奇数滤波窗口构造了9条新的中值滤波曲线特征。另外以地质专家人工标注地层作为标签,确定了8种地层,分别是石炭系、龙潭组、梁山组、茅三段、茅二段、茅一段、栖二段、栖一段。

图6展示了某单井的样本数据测井解释图,测井数据的采样间隔为0.125 m,包含部分测井曲线和地层标签。为了避免不同量纲的影响,利用2.1节中的最大-最小值方法对数据进行了归一化,将所有测井曲线的数值范围缩放到[0,1]内。

图6 well52井的测井曲线及其测井解释

数据预处理完成后,将归一化后的数据排列到一起形成测井曲线矩阵,然后利用2.4节提到的滑动窗口采样方法将其转化为测井曲线切片,通过对长度L及步长S取3～50内的奇数进行对比实验,发现L为5,S=5时识别效果最佳,故设置长度L为5,即每个曲线切片的实际长度为0.625 m,步长S为5。利用该方法最终得到16 635个长度为5,宽为14的二维灰度测井曲线矩阵切片,经统计,样本地地层平均厚度为160个厚度点,因此对每个切片进行样条插值时,选择160作为目标样本长度,最终将其维度变为160×14以保证生成的输入矩阵长度大于特征数目,其中有14 971张属于训练井,1 664张属于测试井,最终处理完成的部分训练样本如表1所示。其中INPEFA、GR_Med_3、GR_Kmeans分别表示INPEFA曲线、滤波窗口为3的中值滤波曲线以及K-Means聚类曲线,Label代表将原来的8种地层标签映射为0～7的数值。

表1 well52井处理后的训练样本数据

3.3 实验环境和网络参数配置

经过数据集的构造,最终由20口井得到了83 018个实际测井曲线切片及其对应的地层标签,每个切片大小为160×14的二维矩阵,然后将数据集按照9∶1的比例划分为训练集和测试集,得到74 716个训练样本和8 301个测试样本。然后对模型进行训练和验证,并寻找最优参数。

对测井曲线进行地层划分的1DCNN网络模型环境配置为Window10操作系统,使用Python3.6作为编程语言,深度学习框架为Tensorflow,GPU使用GTX1650,显存大小为4 GB。

详细的网络配置和参数设置如表2所示。

表2 网络配置和参数设置

实验设置批次大小(Batch Size)为512,训练轮次(Epoch)为100并运用早停减少模型训练时间,优化器使用AdamOptimizer,学习率设为0.001。

3.4 评价指标

引入多个评价指标,利用这些指标对预测结果做全方位评价,综合反映模型分类效果。包括每一类的混淆矩阵(如表3所示)、整体的准确率(A)、精确率(P)、召回率(R)和F1值。

表3 混淆矩阵

(1)混淆矩阵。

(2)准确率。

准确率反映了对于给定样本,模型预测正确的样本数占总样本的比例:

(4)

(3)精确率。

精确率反映了模型预测为正例的结果中真正例所占的比重:

(5)

(4)召回率。

召回率又称查全率,反映了模型预测中样本的正例有多少被正确预测正确:

(6)

(5)F1值。

F1值为精确率和召回率的调和平均值,由于精确率(P)和召回率(R)是一对相互矛盾的量,所以为了更好地评价分类器的性能,一般使用F1-Score作为评价标准来衡量分类器的综合性能:

(7)

公式(4)～公式(7)中:Tp、FN、TN和FP分别为真正例、假负例、真负例和假正例样本的数量,多分类问题可以转换为二分类,对于正例对应的某一类别,其余的类别均为负例。

3.5 模型训练结果

图7是文中模型训练过程中对应的训练集和验证集的准确率上升曲线和损失值下降曲线。在准确率方面,训练集最高为0.98,验证集最高为0.92;损失值方面,训练集最低为0.17,验证集最低为0.28。表4显示了某测试井的分类结果对应的混淆矩阵,可以直观地看出模型在不同地层上的分类效果,其中最高的精确率为茅三段,达到了1,表示全部识别正确,最低的精确率为梁山组,为0.57。表5显示了各类评价指标。由结果可以看出,文中方法对茅口组地层的识别效果较佳,可能与其灰岩岩性有关,在曲线数值上表现为GR数值较低、SP数值适中,比较易于识别。对于其它地层,各类评价指标略低,可能是由于曲线形态具有相似性,但总体的识别准确率都在0.85以上。总之,模型在验证集上的整体预测能力表现较为优异,对于所有测试井的样本而言,整体平均准确率达到了92.82%。

表4 混淆矩阵

表5 各类评价指标

图7 训练过程中模型准确率-损失值曲线

同时,为了验证文中方法的优越性,将其与双向长短期记忆方法[19]和随机森林方法[20]进行对比实验,实验结果如表6所示。结果表明,文中方法的分类准确率为92.82%,高于双向长短期记忆网络模型的88.23%和随机森林方法的84.43%。

表6 各类地层识别方法结果对比

3.6 盲井的地层识别及划分结果

为了验证模型在未经过模型训练的盲井上的地层识别及划分结果,选取某单井数据进行处理并输入模型验证,该盲井的测试样本构造过程同训练集相同,在经过IMK特征构造、滑动窗口采样(长度L=5,步长S=5)及样条插值后得到1 920个160×14的测试样本作为模型输入。模型的输出为1 920×1的一维向量,其中向量中的每个值代表对应测试样本的预测值。图8给出了模型对某测试井的地层划分可视化结果。前4列为训练特征,其中包含GR自然伽马曲线,INPEFA曲线、AC(声波时差)曲线、SP(自然电位)曲线。后两列的柱状图分别代表人工标定的地层标签以及文中模型对该井的地层划分结果。对比发现,文中方法的地层划分结果,与地质专家人工标定的结果比较相似,具有良好的泛化性能,能够很好地识别地层进而划分出不同的地层。

图 8 某单井的测井解释及模型地层识别结果

4 结束语

提出的基于IMK-1DCNN的地层智能识别方法,构建了3种特征,即INPEFA曲线特征、中值滤波特征和K-Means聚类特征。INPEFA曲线能很好地反映不同地层之间的变化趋势;中值滤波特征在保留有效边缘信息的同时很好地去除了测井曲线中的高频噪声;K-Means聚类特征则能很好地反映测井曲线的时空相关性,区分不同地层之间的边界。将文中方法与双向长短期记忆方法和随机森林方法进行对比,文中方法的分类准确率为92.82%,高于双向长短期记忆网络模型的88.23%和随机森林方法的84.43%。文中方法可以实现测井曲线的端到端学习。实验结果表明,文中方法有助于提取测井曲线序列的隐藏特征,准确地划分不同地层。该方法能有效识别地层,具有较强的泛化能力和鲁棒性,为测井曲线地层划分提供了一种新的途径。