基于多信息融合的层次聚类测井曲线自动分层方法

张景越，肖小玲，，王鹏飞，向家富，张翔

（1.长江大学计算机科学学院，湖北 荆州 434000；2.油气资源与勘探技术教育部重点实验室（长江大学），湖北 武汉 430100）

测井曲线分层是在利用测井资料进行岩性识别、测井相分析、沉积相分析、存储参数等研究工作时，首先要完成的工作，正确合理的分层能够保证后续的研究顺利进行［1］。但每个地区的地质情况不同，没有哪一种方法能够适用于所有地区，需要根据实际情况选择合适的测井曲线和分层方法以提高分层精度。

目前测井曲线分层方法分为两大类： 传统的时域分析方法和新兴的频域分析方法。时域分析方法包括层内差异法、极值方差法、活度分层法，以及人工智能方法中的自组织神经网络、人工神经网络等；频域分析方法包括Hilbert-Huang 变换、小波变换等方法。

时域分析方法基于曲线形态本身进行分层。沈禄银等［2］将多条曲线上的信息整合，形成一条包含多曲线信息的综合曲线，然后使用活度分层法对综合曲线进行分层。极值方差法的指导思想是层内差异小、层间差异大，用求微分、斜率极值点在测井曲线上寻找拐点和半幅点［3-4］。阎辉等［5］提出用自组织神经网络对测井曲线进行自动分层。还有一些使用机器学习法［6-8］，其中人工智能的BP 神经网络方法也被广泛使用［9-12］，但使用BP 神经网络首先需要大量的样本对网络进行训练，并且不同地区需要训练不同的样本。

频域分析方法则通过曲线在频域空间中的特性进行分析。覃瑞东等［13］利用Hilbert-Huang 变换得到测井曲线各个固有模态函数的瞬时频率，将测井曲线的瞬时频率变化剧烈的点作为层界点位置，从而实现测井曲线自动分层。小波变换分层方法是将测井曲线进行小波分解，选择中间尺度的分解波，将曲线瞬时频率变化较大的点作为层界点［14-15］。史清江等［16］提出了一种将小波变换和沃尔什变换联合的测井曲线自动分层方法。Pan 等［17］和Mukherjee 等［18］利用傅里叶变换和小波变换相结合的方式进行分层。频域分析方法目前分层效果较好，其缺点在于将一条曲线分解会有很多的分量，选择不同的分量进行分析会得到不同的分层结果，难以选择合适的分量进行分层，而且往往采用1～2 条测井曲线进行分析，使用的测井信息有限。而本文将多种算法融合、数据集融合［19］，能使用更多的测井信息，在加快分层速度的同时也可提高分层精度。

现有的各种地球物理测井方法，由于其仪器设备的纵向探测范围和对应的纵向分辨率各不相同，因此各种不同的测井曲线在分层能力上是不同的。在选取测井曲线方面，应挑选纵向分辨率高且对分层尤其是薄层敏感的测井曲线，通过Fisher 方法选择特征较为明显的测井曲线；随后采用卡尔曼滤波去噪处理，抑制测井曲线上的毛刺和噪点，曲线滤波后使用多信息融合的方法将相似曲线合并； 然后使用层次聚类方法进行地层划分，将低维空间中难以划分的地层信息转化到高维空间中，属于同一地层的物理性质越相似，在高维空间中的距离就越近，通过距离将其划分为一层。

1 层次聚类分层方法

1.1 基本原理

层次聚类算法是无监督聚类算法中最典型的算法之一，主要任务是把一个数据集分成若干个类或簇，分为凝聚法和分裂法2 种算法［20］。本文主要使用凝聚法中的AGNES 算法进行测井曲线的自动分层方法研究。AGNES 算法是由单个个体开始，把单个个体当成不同的类，然后找出距离最小的2 个类进行合并，不断重复到预期类。

1.2 层次聚类的距离度量对结果的影响

距离是层次聚类方法中一个决定聚类质量的关键因素，层次聚类的距离和规则相似度容易定义而且可以聚类成任意形状。距离度量包含最短距离、 最大距离、平均距离和离差平方和距离等。文中使用的距离为离差平方和距离，致力于最小化集群内的总方差。在每一步中，合并距离最近的2 个集群，以最小化集群的相关损失来建立集群。在合并集群的每一步中，算法考虑每个可能的集群并对其进行组合，选择增加信息损失最小的2 个集群进行合并。这里的信息损失是通过ESS（error sum-of-squares criterion）定义的。每个类的离差平方和用公式可表示为

式中：ESSi为第i 类的离差平方和；xj为第j 个点的值；i为类的序号；N 为类的总个数；n 为族群中包含的点的个数；j 为族群中包含的点的序号。

总离差平方和是所有类的离差平方和之和：

式中：ESSsum为总离差平方和。

层次聚类会使用不同的距离度量，这个距离度量决定了不同的类之间距离的计算方式。本文使用1 000条具有3 个特征的数据测试不同的距离度量对聚类效果的影响，其结果分为3 个实际类别（见图1，每个类别用不同颜色表示）。不同距离度量对聚类结果影响的测试数据见图2。

图1 实际类别Fig.1 Actual category

图2 不同距离度量的聚类结果Fig.2 Clustering results of different distance metrics

从图2 可以看出： 使用最短距离与平均距离时的聚类结果较差，大部分数据被划分为一类，少部分数据被划分为其余类别；最大距离的划分结果一般；离差平方和距离的聚类效果最好。聚类结果的优劣也可从表1 中各距离度量下的调兰德指数（ARI）看出：ARI的取值范围为［-1，1］，取值越大越好，反映2 种划分的重叠程度；0 附近的值表示重叠的聚类，负值通常表示样本已被分配到错误的集群，其中离差平方和距离的ARI最高，为0.492 5。因此本文选择离差平方和距离作为不同类别合并的距离指标。

表1 各距离度量下的ARITable 1 ARI for each distance metrics

2 测井曲线选择与数据预处理

使用各种分层方法对测井曲线进行分层之前，要先对测井曲线进行特征优选。测井数据中包含很多曲线，并不是每条曲线都适合用作特定地区的分层，而且有些曲线中还有无效数据，做特征优选时要把无效数据用插值法填充［21］。优选后再进行数据的预处理，用滤波方法，去除测井曲线上的毛刺和噪点。

2.1 测井曲线的选择

在实际工作中，若只使用一两条测井曲线的数据会造成信息量较少，从而导致分层结果不准确。使用特征优选的方式，选择测井曲线中包含主要信息的多条曲线，则能更好地描述问题。实际研究中，PCA 方法会改变原始特征量属性，新构建的主特征量物理意义不一定十分明确，因此，使用Fisher 准则进行特征优选［22］。

本文依据Fisher 值选择得分较高的8 条曲线作为特征曲线，分别是铀（U）、补偿中子（CNL）、声波时差（AC）、地层真电阻率（Rt）、冲洗带地层电阻率（RXO）、自然伽马（GR）、无铀伽马（KTH）、密度（DEN）。

2.2 卡尔曼滤波

一般测井信号会受到仪器或是外界环境的影响，在真实值附近无规律波动，因此最终输出的测井信号由地质的真实信息和噪声组成。测井曲线滤波去噪是为了抹除曲线中的小直径波峰、 波谷和一些突变异常值，以减少噪声对曲线特征的影响。

卡尔曼滤波的基本思想是综合利用上一次的物理量状态和测量值对物理量的状态进行预测估计［23］。在使用卡尔曼滤波器时，通过预测偏差来控制滤波效果。预测偏差越小，滤波效果越明显，曲线越光滑；反之，曲线越接近原曲线。偏差过大或过小，滤波效果都不太好，本文测试了3 个预测偏差值对曲线滤波效果的影响，结果见图3。从图中可以看出：曲线太过光滑时，曲线上的薄层信息被抹除；曲线太接近原曲线时，滤波效果不明显。根据结果，本文选择的预测偏差为0.1。

图3 预测偏差对滤波效果的影响Fig.3 Influence of prediction deviation on filtering effect

2.3 测井曲线的相关性分析和加权因子计算

不同的测井曲线从不同的方面反映地质特性，但选出的曲线中有些曲线可能具有高相关性，即2 条曲线之间的形态相近，表现的物理意义相似，含有大量的相似信息，因此需要对测井曲线作相关性分析。各曲线间的相关性如图4 所示。其中地层真电阻率和冲洗带地层电阻率相关性较高，同样，自然伽马与无铀伽马和铀的相关性也比较高，即多条曲线中含有大量相同信息，使用多信息融合的方式，将相似性较高的曲线融合为1 条。曲线融合处理可以消除观测数据中个别参数畸变的局部影响，突出多项观测数据中能量贡献最大的部分［24］，又可以避免多重相关性的影响。

图4 各曲线间的相关系数Fig.4 Correlation coefficient between curves

假设有L 条测井曲线，每条测井曲线有K 个观测点，则测井曲线上的数据可使用归一化矩阵D 表示：

滤波器输出信号OUT：

式中：x 为加权因子；Dl为第l 条测井曲线值。

输出信号能量PS：

噪声能量PN：

信噪比SNR：

以最终输出信号能量和噪声能量之比最大为目标，保留尽可能多的曲线信息。式（4）中的加权因子x为超参数，可使用遗传算法进行求解［25］。根据求解的加权因子，把GR，KTH 和U 融合为1 条曲线CL1，把RXO和Rt融合为1 条曲线CL2。将2 条融合曲线（CL1，CL2）和CNL，AC，DEN 共5 条曲线作为后续层次聚类分层使用的曲线。

由于常规分层方法只能使用1 条曲线，单一测井曲线包含信息较少，容易产生串层现象，因此本文采用多信息融合的方法［26］，将8 条测井曲线根据相关性进行分类融合，然后将2 条融合曲线（CL1，CL2）同CNL，AC，DEN 一起使用遗传算法求得新的权值。将CL1，CL2，CNL，AC 和DEN 融合成为一条多信息的综合测井曲线，对综合曲线进行归一化、滤波处理，最终将处理后的综合曲线作为活度分层等方法的分层曲线。

3 分层结果对比

3.1 不同测井曲线

对同一井段使用不同测井曲线的分层结果如图5 所示（红线表示划分的层界面）。由图可以看出：单独使用自然伽马曲线进行分层得到的层界点较少，在5 050～5 150 m 的许多小层发生了窜层的情况；使用铀曲线又会划分出过多的薄层。由于单条曲线包含的地层信息较少，使用常规分层方法无法对地层进行有效划分。运用多信息融合方法可将不同曲线上的岩石物理响应特征都保存在一条综合曲线上，在使用综合曲线分层时，能够充分利用多条测井曲线上的地层信息，从而取得不错的划分结果。

3.2 不同分层方法

图6 是本文方法与其他分层方法的分层结果对比。从图中可以看出，仅使用综合曲线的活度分层法和极值方差法能准确划分出一些地层，但当综合曲线某个区间的变化为“急—缓—急”时，活度分层曲线就会出现“极大值—极小值—极大值”的情况，在这个区间内会划分出过多的薄层。深度在5 240～5 270 m 的区间内，5 250 m 处的活度较小，两侧的活度较大，因此在曲线上可划出2 个层界点。

图6 其他分层方法与本文方法分层结果对比Fig.6 Comparison of stratification results of other stratification methods and those of the stratification method in this paper

极值方差法的分层依据和活度分层法较为相似，因此极值方差法也会出现类似的情况（见图6b），而这样的层界点并不符合实际的地层划分。

本文使用多信息融合的层次聚类测井曲线分层方法，将相关性较高的测井曲线融合，既避免了曲线间多重相关性的影响，又保存了地层信息，还减少了计算量。由于属于同一地层各个点的物理性质相似，且高维空间中的距离相近，因此使用聚类方法将同一地层的数据划分为一类，以达到划分地层的目的。由于层次聚类法是在高维空间中使用多条测井曲线进行划分，有些层界面在1 条曲线中无法体现，如5 150～5 200 m 的区间内，使用本文方法可以划分出有效地层 （见图6c），常规分层方法则无法做到这一点。本文方法在曲线走势复杂的地方划分出的地层较多，平缓处划分出的地层较少，符合实际地层划分规律。相对于常规分层方法，本文方法能划分出更多的有效地层，且准确率更高。

各分层方法的分层结果如表2 所示。自动分层（活度分层法、极值方差法、本文方法）的结果与地质分层（人工分层）的结果并不是一一对应的。这是由于自动分层是基于岩石物理响应差异进行的，而地质分层主要基于岩性的差异。不同的岩性可能具有相同的岩石物理响应，而相同的岩性也可能具有不同的岩石物理响应。

由表2 可知，人工分层结果共24 个层界点。活度分层法对应17 个层界点，准确率为70%；极值方差法有58%的层界点与人工分层相近；本文提出的多信息融合的层次聚类法的准确率则达到了83%，表明本文方法的分层结果与人工分层结果较为接近。

3.3 方法验证

某地区以砂岩和页岩为主要的油储层，岩性以泥页岩为主，具有较多的薄互层。图7 是在部分井段应用本文方法的分层结果。多口井的应用表明，本文方法可以有效地提取和保护薄层信息，并且能够在主要的曲线突变处设置分层，地层划分准确率最高可以达到88.89%。

4 结束语

常规测井曲线分层方法，通常采用单一测井曲线进行地层划分，使用的地层信息有限，容易发生窜层现象。本文在数据预处理时首先使用Fisher 方法选取纵向分辨率较高的测井曲线，然后提出多信息融合方法，用滤波器将多条测井曲线进行融合，使多条曲线上的不同岩石物理响应能汇聚在一条综合曲线上，有效地提取和保护了薄层信息。相比于常规分层方法，本文方法能够综合多条测井曲线进行分层工作，使用更多的地层信息，划分出单条曲线上岩石物理响应不明显的地层，有利于克服窜层问题。从应用结果来看，本文方法分层效果较好，地层划分准确率可以达到88.89%，可用于辅助测井曲线的分层工作。