测井曲线的自动分层模型.doc

摘 要
在地球物理勘探中，为了了解地下地质情况，以便于对具有不同特点的地层确定研究目标，以及确定将要重点研究的地层，统一不同井号的研究范围，其中测井曲线分层是首先要完成的基础工作。
本文以1号井为标准井，建立数学模型实现了测井曲线的自动分层。在建立模型过程中，对1号井的数据进行了分类：有效值、无效值、过渡值。我们采用零替换的方法处理了题中出现的无效数据，对于其他非正常数据，由于其表现出的无规律性，因此我们采用了中值滤波的处理方法减小了噪声干扰，从而提高了数据质量。
鉴于测井曲线中评价指标过多的情况，首先根据数据的特点进行了初步的筛选，剔除了信息含量少的指标。对于留下的36个指标，又根据信息论的思想，计算出每一项指标的信息量，进一步剔除信息量较低的指标，最终得到22个测井曲线评价指标。该模型对这22个指标进行了主成分分析，得到五个主成分，其累计贡献率达到了80%以上，起到了降维的作用。再根据主成分的方差贡献率确定了每一个主成分的权重，然后将所有主成分加权求和得到一个新的综合指标，从而根据这一综合指标将所有的测井曲线综合为一条测井曲线，利于模型的后续处理。
对于每号井的综合测井曲线，该模型采用matlab软件编程进行了趋势分析，对测井曲线进行了粗分层，确定了分界点的可能位置，然后进行了层界面归并和加权值法命名，达到了测井曲线自动分层的目的。
依据建立的模型对1号井进行了自动分层，根据分层结果论证了该模型的准确性程度，得出该模型有较高的准确性。然后对2至7号井进行了自动分层，通过了人工分层结果进行对比，分析了在测井曲线分层中出现的自动分层模型的准确度问题和人工分层的主观性问题。最后，利用文中建立的模型对8至13号井进行了自动分层，给出了分层结果，并进行了对结果的分析。

关键字：中值滤波；主成分分析；趋势分析

一、问题重述
在地球物理勘探中需要利用测井资料了解地下地质情况，其中测井曲线分层是首先要完成的基础工作。测井曲线分层的目的是为了在今后的研究中，便于对具有不同特点的地层确定研究目标，以及确定将要重点研究的地层，统一不同井号的研究范围。
通常，在一个区域内，通过前期地质研究工作，结合各种测井数据，首先对最早开发的参考井进行详细研究。每一种测井数据，都反映了地质结构的特点和地层的变化，地质人员通过经验，综合各种测井数据反映的地层特点，将井从一定深度开始，对井进行井层划分和命名，如1号井从距井口深368米处开始，依次往下，定名为长31、长32、长33、长41、长42、长61、长62、长63、长71、长72、长73、长81、长82、长91、长92等地层。接着在分析随后开发的2号井时，也根据和1号井分层的特点和规律，依次定名为长31、长32、长33、长41、长42、长61、长62、长63、长71、长72、长73、长81、长82、长91、长92等地层。井的位置不同可能会导致这口井的每一个层位的深度范围也不同，甚至有可能会出现缺失中间某层的现象。如第6号井缺长31、长32层。通常这些工作都是通过人工来进行的，这就是所谓人工分层方法。该方法不仅费时费力，而且分层取值过程中受测井分析人员的经验知识和熟练程度影响较大，主观性较强，也会因为不同的解释人员的个人标准有误差，而造成不同的人员有不同的分层结果。
自动分层的基本思想、实现手段是一个不断发展变化的过程。由人工分层到自动分层，除了计算机工具的引入，各种数据处理技术也被应用于自动分层。随着一个区域开发井的数量增加，我们希望利用已有分层井点数据与变化特点作为控制点，结合每口井丰富的测井曲线数据，如密度 (DEN)、声波 (AC)、中子 (CNL)、自然伽玛 (GR)、自然电位 (SP) 和电阻率 (RT) 等的变化特点，建立合理的数学模型，实现井位分层人工智能处理，也就是实现自动分层。相对于人工分层，自动分层可以避免人为分层的随意性，并可在很大程度上提高工作效率。进行具体的井位分层人工智能处理，这将极大地提高工作效率。另一方面，希望通过自动分层处理，与人工分层的结果进行比较分析，进一步提高分层精度。
下面请完成以下工作：
1. 以1号井为标准井，根据此井的各种测井曲线数据，建立数学模型，对第2号至7号井进行自动分层，并且通过分析，与人工分层结果进行比较分析。考虑是否需要利用你所建立的数学模型，对1号井的分层结果进行说明。
2. 通过前面人工分层与自动分层的比较结果，以及已给的各种测井曲线数据，确定合适的数学模型对第8号井至13号井进行自动分层，并分析你的结论。
二、问题分析
在地球物理勘探中需要利用测井资料了解地下地质情况，其中测井曲线分层是首先要完成的基础工作。传统做法是采用人工分层方法，但是这种方法不仅费时费力，而且受分析人员的经验知识和熟练程度影响较大，主观性较强，存在明显的误差。我们现在希望利用已有分层井点数据与变化特点作为控制点，结合每口井丰富的测井曲线数据，建立合理的数学模型，实现井位分层人工智能处理，也就是实现自动分层。
第一问中，要求以1号井为标准井，根据此井的已知数据，建立模型，对第2至7号井进行自动分层，并且与人工分层结果作比较。第一步，因为在1号井的测井曲线数据中有66项指标，6000多个样本点，所以我们必须剔除掉一些无关紧要的指标，否则由于数据量太大而无法进行求解。我们通过查阅资料，首先人工初步剔除掉30项指标，然后对剩余的36项指标根据每一项中所包含信息量的大小，再次剔除掉部分指标，最终确定出22项指标作为研究对象。第二步，为了消除因仪器设备引起的测量数据扰动，我们采用中值滤波的方法处理数据，得到更有效的数据。第三步，将数据标准化处理后，采用主成分分析法将多个指标综合成少数几个综合指标，即确定出几个主成分。第四步，根据每一个主成分对分层结果的贡献率确定出每个主成分的权重，再将所有主成分综合成一个综合指标，这样每一口井可以得到一条关于这个指标的测井曲线。第五步，我们对测井曲线进行趋势分析，通过计算测井曲线上各点的切线斜率，可以判断层界面的大体位置，从而完成初步分层，然后对初步分层中不必要的小薄层利用层界面归并归并为一层。第六步，我们通过计算每口井每一层的测井值，利用距离判别的思想，分别于1号井作比较，完成对井层的命名。
最后，我们通过作图将第2至7号井自动分层结果与人工分层结果作比较，对比较结果作出分析，并且对1号井的分层结果进行说明。
第二问中，根据第一问中已经建立的模型，结合已有的第8至13号井的测井数据，对第8至13号井进行自动分层，并且对分层的结果作出分析。
三、模型假设
1．测井数据真实可靠，不考虑人为因素引起的误差；
2．不考虑地层变化引起的分层波动；
3．以1号井为标准，它包含所有的地质分层情况。
四、符号说明
符号
含义

第项指标的信息量

第项指标的标准差

第项指标的均值

标准化后第项指标第个观测值

随机向量的均值

随机向量的协方差矩阵

主成分分析中相关系数矩阵

主成分的特征值

主成分的特征向量

五、模型建立的准备
5.1 数据的预处理
我们观察分析每口井的测井数据，发现每项指标测得的数据值可以分为以下几类，第一类，无效值，如-9999；第二类，有效值，如30.885；第三类，过渡值，此类值往往出现在从无效值到有效值的过渡段，如-6196.22。我们考虑到无效值会影响分层计算的结果，而且无效值都是相同的值，为了便于计算，将每口井的测井数据中无效值全部替换为0。而对于过渡值，因为它没有规律可循，无法一一替换，我们将在模型中采用中值滤波的方法消除它带来的噪声干扰。
六、模型的建立与求解
6.1 问题一模型的建立与求解
6.1.1 人工初步剔除指标
在1号井的测井数据中，一共有66项指标，由于数据量太过庞大，计算机无法在短时间内算出结果，所以必须简化指标。我们通过查阅资料，人为初步剔除了30项指标，剔除后的36项指标，如下表所示，
表1 人工初步剔除后的36项指标
DEN
RILD
RILM
RLL8
SP1
R4.0
SP
DEVi
AZIm
GR
AC
RML
RMN
CAL
CNL
R4.0%
cnl%
GR%
AC%
RILD%
RILM%
RLL8%
SP1%
RMN%
SP%
den%
CAL%
RML%
PORT
SW
SH
SXO
RT
CID
CIM
CIL8

6.1.2 根据信息量的大小再次剔除指标
在信息论中，如果某项指标的变化幅度越大，我们就认为这项指标所携带的信息量越大，对结果的影响越大，也就是说这项指标越重要。而在1号井的测井数据中，测井开始和测井结束时每项指标的数据都是无效的，这些数据对分层不能提供有效地信息，所以我们取从长31的底深（294）到长92的底深（857.8）这一段测井数据作为每项指标的有效观测值，然后我们定义某项指标所有有效观测值的标准差与均值的比值大小为该指标的信息量，即
.…………………………………… (1)
其中，表示第项指标的信息量，表示第项指标的标准差，表示第项指标的均值。
我们在Excel中利用Average函数和STDEVP函数分别求出1号井人工初步剔除后每项指标的均值和标准差，然后代入(1)式，求得每项指标的信息量，作出如下图：

图1 36项指标的信息量
图中横坐标表示每一项指标，纵坐标表示每一项指标对应的信息量，从中可以看出有些指标的信息量非常小，完全可以剔除掉，经过再次剔除后剩余22项指标，它们的信息量如下图：

图2 22项指标的信息量
从图2中可以看出，这22项指标的信息量基本都大于0.5，我们认为这些指标携带着影响分层的主要信息，是比较重要的指标。
最终选取的22项指标如下表所示：
表2 根据信息量大小再次剔除后的指标
RILD
RILM
RLL8
R4.0
DEVi
AZIm
RML
RMN
CNL
R4.0%
cnl%
RT
RILD%
RILM%
RLL8%
CID
CIM
CIL8
SH
RMN%
RML%
PORT

6.1.3 中值滤波
为了消除因仪器设备产生的测量数据扰动，我们采用中值滤波对数据进行处理。
假设有数据样本序列，当滤波窗口为，其中值滤波的步骤为：
（1）取以第个数据为中心的个数据行进行排序(顺序或逆序均可)；
（2）取排序后的中间值，即第个数据作为第点的滤波值；
（3）自上而下迭代计算离散曲线上的各点。
中值滤波流程如图3所示

图3 中值滤波流程图
输入;

调用快速排序，得

输出
开始
结束

否
是

滤波窗口中的的的视目的不同而定。时，滤波器失效；较小（如）时可保持曲线的幅度值；较大时可只保持曲线的变化趋势。在这里，我们没有找到很好的办法确定的大小，只能经过对不同值的尝试，最终确定当时，比较合理。
6.1.4 标准化处理数据
为了消除各项指标之间不同量纲的影响，对每项指标的观测值进行标准化处理，即：
.…………………………………… (2)
式中，，，表示第项指标的第个观测值，表示第项指标所有观测值的均值，且，表示第项指标所有观测值的标准差，且，表示标准化处理后第项指标的第个观测值。
根据式（2），我们对所有指标的观测值进行标准化处理。
6.1.5 主成分分析
主成分分析法是利用降维的思想，在力保数据信息损失最少的原则下，把多个指标转化为少数几个综合指标的一种对多变量数据进行最佳综合简化的多元统计方法。
设它的数据有个指标，,分别用表示，这个指标构成的维的随机向量为。设随机向量的均值为，协方差矩阵为。对进行线性变换，形成新的综合变量,满足下式：

其中

在上面各式中，表示第个主成分，表示第个主成分中，第个指标的载荷。
　　主成分分析的基本步骤如下：
第一步：求协方差矩阵。对于1号井，有22项指标需要处理，则先由原始数据得到协方差矩阵。
第二步：指标的标准化处理。利用5.1.4中标准化处理后的结果，可得到新的数据矩阵。
第三步：计算相关系数矩阵。其中，为原来变量与的相关系数，其计算公式为：

因为是实对称矩阵（即），所以只需计算其上三角元素或下三角元素即可。
第四步：计算特征值与特征向量。解特征方程求特征值。因为为正定矩阵，所以其特征值都为正数，将其按大小顺序排序，即。特征值是各主成分的方差，它的大小反映了各个主成分在描述被评价对象上所起的作用；然后根据方程可确定特征向量的矩阵。
第五步：计算主成分贡献率及累计贡献率，确定主成分个数。主成分的贡献率为，累计贡献率为。一般取累计贡献率达80％至95％的特征值所对应的个主成分。
我们通过Matlab编程按照上述步骤求解，最终得到的每个主成分中各个指标的系数矩阵如下表所示：
表3 每个主成分中各个指标的系数
主成分
指标
1
2
3
4
5
CIL8
-0.3133
0.0688
-0.015
0.0681
-0.0251
CIM
-0.3035
0.0884
-0.022
0.1307
0.0451
CID
-0.2208
0.1906
-0.0754
0.1907
0.1135
RT
-0.1613
0.2042
-0.071
0.3788
0.0157
SH
0.0106
-0.0507
-0.67
-0.181
0.0165
PORT
0.0358
-0.0074
-0.6848
-0.076
0.1351
RMN%
-0.2505
-0.0996
-0.0182
0.2945
0.0617
RML%
-0.2485
-0.1373
0.0274
0.253
0.0546
RLL8%
0.0721
0.4417
-0.089
0.1052
0.1163
RILM%
-0.0731
0.1732
0.014
0.1389
-0.0008
RILD%
0.0733
0.4413
-0.0925
0.1206
0.091
cnl%
-0.2355
0.1872
0.057
-0.3493
-0.2522
R4.0%
-0.2237
0.214
0.0631
-0.3794
-0.233
CNL
-0.1587
0.2643
-0.0012
-0.3353
-0.2717
RMN
-0.2198
0.0124
0.0846
-0.2713
0.5864
RML
-0.2151
-0.015
0.1267
-0.2969
0.574
AZIm
-0.075
-0.367
-0.0282
-0.0276
-0.0166
DEVi
0.1457
0.368
0.0402
0.0877
0.0444
R4.0
-0.3133
0.0688
-0.015
0.0681
-0.0251
RLL8
0.2811
0.1003
0.0762
-0.0005
0.1826
RILM
0.2785
0.1223
0.0913
-0.0306
0.1659
RILD
0.2938
0.0856
0.0759
-0.1003
0.0918
在表3中，每一行表示同一个指标在不同主成分中的系数，每一列表示同一个主成分中不同指标的系数，因为主成分是原始评价指标的的线性组合，在这个线性组合中各变量的系数有大有小、有正有负。一般而言，线性组合中系数绝对值大的指标表明其对该主成分的属性做出了较大的贡献，若几个指标系数相当时，则应认为这一主成分是这几个指标的综合。
6.1.6 综合曲线
确定出综合评分函数，即各个主成分的加权线性组合。

其中，为第个主成分的权值，表示第个主成分，表示综合评分值。
计算解得权向量。
6.1.7 利用趋势分析划分井层
趋势分析是一种既直观又有效的分层方法，特别是对于层中心对称的曲线有很好的分层效果，通过计算测井曲线上各点的切线斜率，可以判断层界面的大体位置。
6.1.7.1 斜率的计算
设离散的测井曲线上有对样本

表示测井样本的实测值，表示样本的深度，则曲线上点的切线斜率可近似为：

为便于计算判断，令，对于等间距的采样点，为定值，可设，所以。
可见，当时，曲线处于拐点处；当较大时，曲线的斜率较小。给定门限值，当时，可认为实测值有显著变化，它一般处于层界面上。
6.1.7.2 趋势分析的流程

图4 趋势分析的流程图
输入;；

计算
输出各层界面
开始
结束
与异号

标志作为新界面的始点

否
否
否
是

6.1.8 层界面归并
在5.1.3中，我们为了消除因仪器设备造成的测量数据扰动，采用中值滤波对数据进行处理，但中值滤波也无法保证滤掉所有由仪器测量误差所产生的尖峰，这将在分层计算中分出一些不必要的小薄层，这时需要利用归并薄层界面来消除。
薄层归并一般采用扫描迭代的方法。将曲线看成一深度序列，相邻两点的深度差为：

对于给定的层厚最小值，当时，将此分层并入相邻的顶层（或相邻底层）。
6.1.8 井层命名
在确定出层界面之后，接下来就是要计算层测井值。即又一层内的多个采样点的测井值确定出能反映该层物理特征的一个测井值。然后根据此测井值，与标准井1号井的测井值作比较，从而对井层命名。
我们采用取加权值的方法，具体的做法是，找出层内的个极值，以极值间的距离（厚度）作为权，进行加权平均，即：

取作为本层测井值。
在命名中，首先计算出1号井每一层的测井值，再计算出其余任一口井每一层的测井值，然后与1号井的测井值相比较，运用距离判别的思想，判断出待命名井各层与1号井各层的对应关系，即完成命名。
综合上述所有的模型，并结合已有的测井曲线，对第2号井至7号井进行自动分层。
最终的分层结果如下表所示：
表4 第2至7号井分层结果

2号井
3号井
4号井
5号井
6号井
7号井
长31底深
318.57
298.52
302.33
307.35
300.48
305.98
长32底深
357.15
308.91
321.4
314.62
315.41
349.57
长33底深
383.5
327.1
344.8
322.77
334.93
416.75
长41底深
429.57
378.6
383.2
350.3
388.36
464.33
长42底深
514.77
422.92
408.78
368.2
406.43
487.78
长61底深
544.94
450.7
464.95
384.1
453.08
519.98
长62底深
574.4
464.45
494.07
397.35
522.3
549
长63底深
637.2
529.94
517.4
407.1
554.43
583.88
长71底深
667.96
549.15
549.9
447.1
577.24
651.1
长72底深
690.85
573.03
662.6
451.9
604.76
674.4
长73底深
722.7
637.53
688.3
469.6
629.18
711.4
长81底深
779.35
680.3
706.3
849.8
656.61
754.72
长82底深
800.6
722.93
735.17
895.2
717.97
793
长91底深
847.4
771.96
803.82
900.2
764.93
850.8
长92底深
874
801.64
854.65
955.7
806.88
893.55
6.1.9 对1号井的分层结果进行说明
我们为了检验模型的正确性，采用直观的柱形图对一号井人工分层的结果与自动分层的结果作比较，如下图5所示：

图5
在图中，横坐标表示井层，纵坐标表示每一井层的底深。从图中可以看出，人工分层和自动分层的结果差异不大，这说明我们建立的模型可靠，准确性高。
6.1.10 对比人工分层与自动分层结果
我们将第2至7号井的人工分层结果与自动分层结果用直观的柱形图进行对比，如图6—图11所示：

图6 图7

图8 图9

图10 图11
上面每一幅图中，横坐标均表示井层，纵坐标均表示每一井层的底深。从上面几幅图中可以看出，人工分层结果与自动分层结果大体相同，但部分井仍有明显的差异，这一差异应该源自于两方面，
第一，人工分层受测井分析人员的经验知识和熟练程度影响较大，主观性较强，也会因为不同的解释人员的个人标准有误差，而造成不同的人员有不同的分层结果，所以人工分层的结果不是完全可靠。
第二，自动分层模型本身存在缺陷导致模型的准确性有限，因为在建立模型时，所采用的每一种基本模型、算法，都存在缺陷，所以最终导致模型的准确性有限。
正是上面两方面原因导致人工分层结果与自动分层结果存在差异。
6.2 问题二模型的建立与求解
6.2.1 对第8至13号井分层
在问题二中，根据问题一中确定的模型，以及已给的各种测井曲线，对第8号井至13号井进行自动分层。
最终分层的结果如下表所示：

表5 第8至13号井分层结果

8号井
9号井
10号井
11号井
12号井
13号井
长31底深
303.13
299.45
302.67
316.15
303.25
305.7
长32底深
316.8
315.1
341.86
350.59
326.37
322.13
长33底深
383.8
329.3
365.18
389.73
369.68
338.25
长41底深
427.45
352.3
395.43
437.68
396.58
373.7
长42底深
464.1
376.6
419.18
460.93
446.1
390.3
长61底深
485.8
447.2
456.74
480.3
484.05
411.9
长62底深
505.1
471.3
491.84
509.36
512.5
441.4
长63底深
516.8
498.7
522.72
573.55
529.5
472.38
长71底深
550.65
551.7
588.11
612.3
548.7
501.25
长72底深
585.8
574.15
635.09
648.55
606.3
532.35
长73底深
674.3
602.6
655.55
670.45
641.35
612.15
长81底深
696.3
629.4
680.99
692.61
661.7
731.8
长82底深
716.3
654.9
703.05
766.93
691.6
760.05
长91底深
730
710.2
763.51
794.28
734
790.2
长92底深
824
810.95
803.9
851.57
800.2
843.1
6.2.2 对第8至13号井分层结果的分析
从表5的分层结果可以看出，由于第8至13号井基本处在同一区域，所以它们每一层的分层结果相近，但个别井有明显的不同。从表中还可以看出，每一口井不同层的层深在一定范围内，这一特征与标准井是符合的，所以是合理的。

六、模型的评价
6.1模型的优点
1. 合理的简化指标数量。由于每一口井的数据量都非常庞大，无法对全部进行分析计算，所以只能剔除掉一些无关紧要的指标，以减少数据量。我们首先人工初步剔除掉30项指标，再根据各个指标携带的信息量大小，剔除掉信息量小的指标，最终筛选出22项指标作为研究对象。
2. 综合所有主成分得到一个综合指标。通过主成分分析法得到若干个主成分后，我们将每个主成分对结果的贡献率作为各自的权重，再将所有主成分加权组合成一个新的综合指标，而每一口井都能得到唯一一条关于这项指标的测井曲线，然后只需根据这条综合曲线进行分层。
6.2模型的缺点
1. 在测井数据的处理中，只是对无效值替换为0，对于过渡值没有找到很好的处理方法，所以无法消除过渡值对分层结果的影响。
2. 在对数据做中值滤波处理时，需要确定滤波窗口大小，我们只能定性分析出窗口大小对结果的影响，但是不能定量的确定出窗口大小，只能通过若干次试验选取，缺乏依据。
3．在对综合曲线做趋势分析时，无法确定出门限值，只能根据1-7号井的已有分层对门限值进行估计。

参考文献
[1] 纪荣艺，樊洪海等，测井曲线自动分层模型设计与实现，石油钻探技术，2007年第35卷第2期
[2] 谢忠怀，李保利，测井地层自动化分的边缘检测最优分割方法，新疆石油地质，2006年27卷第2期。
[3] 陈文鹏，杨绍贵，武帅杰，测井曲线自动分层
[4] 肖波，韩学辉，周开金，支乐菲，测井曲线自动分层方法回顾与展望，地球物理学进展，2010年25卷第5期。

附录
中值滤波及归一化程序代码：
clear;
clc;
x=xlsread('11号井.xls','sheet1','T2:T5522');
c=medfilt2(x,[4,1]);
a=mean(c);
sig=std(c);
for i=1:length(c)
b(i)=(x(i)-a)/sig;
end
[a,b]=xlswrite('11号井.xls',’b’,'sheet2','T2:T5522');

粗分层层界点选择程序代码：
clear ;
clc;
global N
N=input('请输入门限值N=');
x=xlsread('13口井的综合测井曲线.xls','sheet1','I4:I7722');
j=1;
std(x)
for i=2:length(x)
if (abs(x(i))>N)&(x(i)*x(i-1)<0 )
b(j)=i*0.125+294.05;
j=j+1;
end
end
b'
plot(b);