煤矿开采过程中，断层不仅改变了煤(岩)层的埋藏条件，而且使煤(岩)层错断并发生显著位移，这一方面破坏了煤层的连续性和完整性，为煤层开采带来阻力；另一方面，断层处容易发生瓦斯突出、透水、冒顶等地质灾害，严重影响矿区的安全开采。因此，提前探明矿区断层分布是煤矿安全开采的重要内容。目前，地震领域，断层的解释手段主要是人机交互解释，通过解释人员肉眼观察地震波的振幅、相位和时差等特征以确定断层的存在，解释可靠性在一定程度上取决于解释人员对工区有关褶皱、断裂等构造模式的掌握程度，依靠这种方法解释断层具有很大的局限性，并且是一个费时费力的过程^[1]。

为了打破传统断层解释方法的局限性，近年来，不少学者利用属性不连续性来表征断层分布。M. Bahorich等^[2]首先提出了相干体属性，并获得了三维数据中地层的不连续性特征；N. M. Albinhassan等^[3]将霍夫变换应用于时间切面，以增强断层的显示；S. I. Pedersen等^[4]使用蚂蚁跟踪来增强空间不连续性并改善地震数据中的断层特征；F. Admasu等^[5]提出了一种用于三维地震数据的半自动断层跟踪方法，该方法涉及使用log-Gabor滤波器来增强断层振幅的不连续性，并在地震剖面上跟踪断层。这些方法通过完善边缘检测属性来增强断层响应，突出显示断层边缘通常取决于属性的质量。Lu Cai等^[6]开发了一种体绘制技术，融合体渲染多个地震属性，向解释人员显示了三维数据体内部结构的直观视图；孙振宇等^[7]利用支持向量机(Support Vector Machines，简称SVM)进行多属性断层识别。这些研究，通过地震多属性解释很好地避免了对单一属性质量的依赖。在利用多种属性进行地震解释的同时，还需要对其进行数据整合，抓住数据中有效信息，舍弃干扰信息^[8-10]。

主成分分析(Principal Component Analysis，简称PCA)是一种使用广泛的数据整合方法，被广泛应用于信号处理、统计等各个领域。这种方法通过分析提取样本的少量特征来降低空间维数，在降维过程中产生的新特征向量为正交向量，向量之间相互独立，可以帮助我们抓住主要信息，尽可能地消除噪声等干扰^[11-12]。同时，近年来兴起的机器学习方法具有很好的数据泛化能力，可以从已经得到的数据中获取信息，达到分类或者回归的目的。庞大的地震数据量为机器学习在该领域的应用和发展提供了先决条件^[13]。最近邻算法(k-nearest Neighbor Alogorithm，简称kNN)是机器学习领域最常用的算法之一。由于实现简单，理论清晰和分类性能出色，它已在许多领域得到广泛使用^[14]。同时，传统的kNN算法需要大量的存储空间对训练样本进行存储，在多维大数据量的情况下，这意味着庞大的运算量和运算时间^[15]，经过PCA降维，也很大程度上解决了这方面的不足。

笔者采用PCA方法对多种地震属性进行整合，产生新的综合属性，然后将其作为属性输入并利用kNN算法进行断层识别，同时通过交叉验证评估断层识别的准确率，为地震构造解释的研究提供了一种新的思路。

1   基本原理和方法

1.1   主成分分析(PCA)

一般来说，想要获得低维子空间，最常用的方法是对原始的高维空间进行线性变换，给定d维空间中的样本集$ X = \left\{ {{x_1}, {x_2}, \cdots , {x_m}} \right\} \in {R^{d \times m}} $，变换之后得到$ d' $维的空间样本，其中，$ d' \leqslant d $。

式中：$ Z \in {R^{d \times m}} $是样本在新空间中的投影；$ {\boldsymbol{W}} \in {R^{d \times d'}} $是变换矩阵。

基于线性变换来降维的方法都符合式(1)的基本形式，不同的线性降维方法之间的差别在于对低维子空间的性质有不同的要求，也就是对$ {\boldsymbol{W}} $施加了不同的约束。

主成分分析作为一种常用的线性降维方法，其要求低维子空间对样本具有最大可分性。这种方法最初是由K. Pearson于1901年对非随机变量引入的。1933年，H. Hotelling进一步完善了PCA的数学基础，将其推广到随机向量^[16-17]。其核心思想是通过坐标旋转，寻找新的正交基，从而将数据投影到使数据方差最大的n维坐标轴上，得到数据在新坐标系中的表示以消除原数据空间的多重共线性，从而达到数据降维的目的^[18]。这种方法的主要流程如图 1所示。

图  1  PCA算法流程

Figure  1.  PCA algorithm process

下载: 全尺寸图片 幻灯片

降维后低维空间的维数$ d' $通常是由用户事先指定，根据实际问题，研究人员常常根据累计方差贡献率大于80%这一指标来确定主成分个数。在这一过程中，特征值最小的$ d - d' $个特征向量被舍弃了，这是降维导致的结果。通过舍弃这部分信息，一方面能够使样本的有效信息密度增大；另一方面，当数据受到噪声影响时，最小的特征值所对应的特征向量往往与噪声有关，将其舍弃能在一定程度上起到去噪的效果^[19]。

1.2   最近邻算法(kNN)

最近邻算法(k-Nearest Neighbor)是T. Cover和P. Hart于1967年提出一种常用的监督学习算法^[20]。kNN的工作机制是：给定测试样本，基于某种距离度量找出训练集中与其最靠近的k个训练样本，然后基于这k个“邻居”的信息来进行预测。通常，在分类任务中可使用“投票法”，即选择这k个样本中出现最多的类别标记作为预测结果。

图 2是最近邻算法的原理。从图中可以看出，k是一个重要参数，当k取值不同时，分类结果会显著不同；另一方面，若采用不同距离的计算方式，则找出的“近邻”可能有显著差别，从而也会导致分类结果的不同^[21]。

图  2  最近邻算法的原理

Figure  2.  Principle of the k-nearest neighbor algorithm

下载: 全尺寸图片 幻灯片

常用的距离计算方法有欧式距离、曼哈顿距离和闵可夫斯基距离等^[22]。本文采用欧氏距离计算方法：设特征空间$ X $是$ d $维实数向量空间$ {{\boldsymbol{R}}_d} $，$ {x}_{i}, {x}_{j}\in X，{x}_{i}={({x}_{i}^{(1)}, {x}_{i}^{(2)}, \cdots , {x}_{i}^{(d)})}^{\rm{T}}, {x}_{j}=({x}_{j}^{(1)}, {x}_{j}^{(2)}, \cdots , x_j^{(d)}{)^{\rm{T}}}, {x_i}、{x_j}$的距离$ {L_0} $可以用下式计算^[23]。

从式(2)中可以看出，特征空间$ X $的维度$ d $越大，则该公式计算的繁琐程度越高，而测试样本在利用kNN进行预测时，需要与每个训练样本都进行距离计算，然后通过排序找到前k个最近邻，这意味着极大运算量和运算时间。利用主成分分析的方法可以将$ d $维向量降维为$ d' $，使得每次距离计算的计算量大大减少，弥补了kNN这方面的不足。

2   数据集构建

2.1   研究区概况

峰峰矿区羊东煤矿位于中国河北邯郸。本矿区内含煤地层为石炭系上统本溪组、石炭–二叠系太原组和二叠系下统山西组，井田含煤地层平均总厚约210 m，含可采煤层6层，即2、4、6、7、8、9号煤，本次研究的目的煤层为2号煤层。2号煤层顶板为粉砂岩，底板为含泥质及炭质粉砂岩，属稳定可采煤层。

研究靶区的勘探面积为3 km²，工作面内已有3条巷道、15口钻井，其在矿区内的分布情况如图 3所示。根据巷道和钻井提供的坐标位置与构造特征(是否存在断层)对应关系信息，在2号煤层确定139个点的坐标及该位置所对应的构造特征(表 1)。

图  3  羊东煤矿2号煤层已有钻井、巷道分布情况

Figure  3.  Distribution of existing drilling wells and roadways of No.2 coal seam in Yangdong Coal Mine

下载: 全尺寸图片 幻灯片

表  1  139个点位置坐标及其对应的构造信息

Table  1.  Position coordinates and their corresponding geological information of the 139 points

序号	X坐标	Y坐标	是否存在断层	信息来源		序号	X坐标	Y坐标	是否存在断层	信息来源
1	525992	4039808	是	巷道1		128	524834	4039530	是	井1402
2	525992	4039808	否	巷道1		129	525672	4039327	否	井1404
$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $		130	525337	4039388	否	井148
44	525976	4039789	是	巷道2		131	524960	4039741	否	井141
45	525944	4039817	否	巷道2		132	525224	4039824	是	井144
$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $		133	525772	4039573	否	井145
82	525958	4039766	是	巷道3		134	525386	4039710	否	井149
83	525927	4039795	是	巷道3		135	525475	4040141	否	井1301
$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $		136	525665	4040023	否	井136
125	524887	4039379	是	井155		137	525054	4039979	是	井1303
126	525257	4039147	是	井156		138	525225	4040225	是	井130
127	525592	4039078	否	井157		139	525285	4040421	否	井1204

下载: 导出CSV 

| 显示表格

2.2   特征提取

特征提取的目的是获取研究对象的尽可能多的信息。地震属性是由叠前或叠后地震数据，经过数学变换而得出的有关地震波的几何学、运动学、动力学或统计学特征，不同的地震属性代表不同的物理特性。因此，地震属性是非常合适的特征提取对象^[24]。以方差体为例，方差体的计算是通过求取加权移动的方差值，得到三维数据体中每个时间样点的方差，其计算如式(3)所示。方差值越大，说明相似性越差，常被用于表征地层的不连续变化及检测地下断层。对三维地震数据体提取方差属性，包含了层位不连续信息。

式中：$ {\delta ^2} $为当前采样点的方差值；$ t $为当前采样点的时刻值；$ {\omega _{j - t}} $为三角形权重函数；$ {m_{ij}} $表示第$ i $道$ j $时刻采样点的振幅值；$ \overline {{m_j}} $为所有$ i $道数据在$ j $时刻的平均振幅值；$ {\omega _{ij}} $为第$ i $道$ j $时刻采样点对应的三角形权重函数；$ L $为时窗的长度；$ I $为参与计算的地震数据的道数。

对数据体分别提取了10种地震属性，包括：方差、衰减系数、走向曲率、反射强度、瞬时相位、最大振幅、瞬时频率、倾角偏差、倾角连续性、混沌体，这10种属性均能用于表征断层^[7]。然后，将已知构造信息的139个点与提取得到的地震属性根据坐标位置对应起来，构建羊东煤矿2号煤层已知数据集，见表 2。其中，非断层数据97组，用标签‘0’表示；断层数据42组，用标签‘1’表示。

表  2  羊东煤矿2号煤层已知数据集

Table  2.  Known data on geological information of the No.2 coal seam in Yangdong Coal Mine

方差	衰减系数	走向曲率	反射强度	瞬时相位	最大振幅	瞬时频率	倾角偏差	倾角连续性	混沌体	label
0.03	0.05	0.02	535.46	110.37	220 983.1	43.39	0.02	0.08	0.05	1
0.01	0.02	0.01	1 211.85	106.04	72 728.63	46.09	0	0.01	0.01	0
$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $
0.02	0.05	0	410.21	116.25	108 717.8	37.72	0	0.08	0.01	0

下载: 导出CSV 

| 显示表格

2.3   特征选择

对于从地震数据体提取得到的属性数据集，需要利用PCA对其中的原始地震数据进行整合，整合得到的综合属性将作为特征选择的结果，用于kNN进行断层识别。本文利用SPSS软件对数据进行主成分分析，具体处理过程如下。

输入的原始地震属性数据集由47 442个样本点组成，每个样本点有10种地震属性值，则输入样本集$ X = {({x_{ij}})_{47\;442 \times 10}} $，其中，$ x_{ij} $表示第$ i $个样本点的第$ j $个属性值。由于从线性变换的本质来说，PCA是在线性空间做一个旋转，然后取低维子空间上的投影点来代替原样本点，以达到降维的目的，所以首先要保证原本空间里的点是以原点为中心分布的，同时考虑到不同的地震属性之间具有不同的量纲，数量级差异较大，因此，PCA的第一步是对每组地震属性分别进行标准化处理，其计算原理如式(4)所示。标准化处理后，各组地震属性$ z_j $的均值为0，方差为1，构成了以原点为中心分布且不同维度的属性具有同等重要性的地震属性簇。

式中：$ z_{ij} $为标准化处理后输出的地震属性值；$ {x_{ij}} $为输入的原始地震属性值；$ \overline {{x_j}} $为$ {x_j} $的均值；$ {\sigma _{{x_j}}} $为$ {x_j} $的标准差。

对原始数据标准化后，得到了输出矩阵$ \boldsymbol{Z} = {({z_{ij}})_{47\;442 \times 10}}$。为了得到方差变化最大的方向，需要计算矩阵$ \boldsymbol{Z} $的各个列向量之间的协方差，根据式(5)可以求得其协方差矩阵$ \boldsymbol{C} $。

其中， ${c_{mn}} = \frac{1}{{47\;442}}\sum\nolimits_{k = 1}^{47\;442} {{z_{mk}}{z_{kn}}(m、n = 1, 2, \cdots , 10)} $。

利用特征方程$ \left| {C - \lambda E} \right| = 0 $，可以求出协方差矩阵$ \boldsymbol{C}$的特征值$ {\lambda _j}(j = 1, 2, \cdots , 10) $。同时，定义第$ j $个特征值对应方差贡献率为$ a_j $，则：

各主成分特征值λ_j和及其对应的方差贡献率a_j的计算结果见表 3。计算得到的协方差矩阵的特征值就等于对应主成分的方差，其大小反映了第j个主成分所包含原始数据全部信息的权重，也反映了各主成分贡献的大小。根据计算结果，确定主成分的数目为6，累计方差贡献率达80.503%，包含了大部分有效信息的同时也达到了降维的目的。

表  3  主成分特征值及其方差贡献率

Table  3.  Principal component eigenvalue and its variance contribution rate

主成分编号	特征值	方差贡献率/%	累计方差贡献率/%		主成分编号	特征值	方差贡献率/%	累计方差贡献率/%
1	2.559	25.585	25.585		6	0.916	9.161	80.503
2	1.486	14.862	40.447		7	0.736	7.356	87.859
3	1.141	11.407	51.854		8	0.552	5.517	93.376
4	0.999	9.993	61.847		9	0.493	4.927	98.303
5	0.950	9.495	71.343		10	0.170	1.697	100.000

下载: 导出CSV 

| 显示表格

确定完主成分后，利用主成分的特征值可以求出对应的特征向量$ {\mathit{\boldsymbol{\gamma }}_j} = ({\gamma _{j1}}, {\gamma _{j2}}, \cdots , {\gamma _{jp}})(p = 1, 2, \cdots , 6) $。然后以主成分特征向量的分量值为权数，根据式(7)，对标准化后的属性数据进行加权就得到了第$ p $个主成分的综合属性$ Y_p $，结果见表 4。

表  4  综合属性数据集

Table  4.  Integrated attribute data set

样本点	Y1	Y2	Y3	Y4	Y5	Y6
1	0.623 25	–0.003 02	–0.798 82	1.624 09	1.256 43	–1.289 87
2	0.604 85	0.264 97	–0.836 58	0.084 20	–0.408 08	–1.456 64
3	0.475 82	0.193 16	–0.979 45	0.866 78	–0.336 68	–1.396 93
$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $	$ \vdots $
47 442	–0.169 09	1.711 36	–0.007 00	–3.548 22	0.643 57	–0.643 87

下载: 导出CSV 

| 显示表格

3   模型构建与断层自动识别

3.1   交叉验证的思想

在对kNN算法分类准确率进行评估时，通常的方法是将已知数据集随机分为两部分：训练集和测试集。这种方法的样本选取虽然具有随机性，但不同的分类方法训练的结果可能有差异。为了消除这种影响，本文在样本选取时采用了十折交叉验证的方法：将用于分类的139组数据，分成10组，其中每组数据数量如图 4所示，每次实验，依次分配训练集9份验证集1份，10次结果的均值作为算法精度的估计^[23]。该方法的优点是：所有的样本都进行了训练，每一份样本又各自作为测试集进行验证，提高了预测的精确度；同时，也降低了分组不同的影响。

图  4  每组数据数量

Figure  4.  The number of data set in each group

下载: 全尺寸图片 幻灯片

3.2   模型构建与评估

k值的选取是构建kNN分类模型的重要参数，在过拟合和欠拟合之间保持平衡是选择k值的关键，k值的合理选取可以尽可能减少噪声对输出类别的影响^[25]。根据经验规则，k一般取奇数且小于训练样本数的平方根，考虑到样本集数据量的大小，本文k的取值1、3、5、7、9、11，然后依据验证集准确率大小随k取值的变化规律，确定模型的最佳k值。

为了更好地观察PCA提取特征前后对kNN分类效果的影响，分别将PCA处理前后的2组数据集作为输入，并利用十折交叉验证的思想和图 4的分组方式对2组数据进行分组，然后利用kNN对2组数据进行分类。验证集分类准确率随k取值的变化如图 5所示。

图  5  各组验证集分类准确率随k取值的变化

Figure  5.  Variation of classification accuracy of each verification set with the value of k

下载: 全尺寸图片 幻灯片

以分组1为例，从图 5中可以看出：在该种分组下，未经过PCA提取主成分的数据集在k取值3、5、9、11时，验证集分类准确率最高，为84.61%；同样的分组下，经过PCA提取主成分的数据集在k取值1、9、11时，验证集分类准确率最高，为92.31%。同样，其他各组均选择验证集准确率最高时的k值作为当前分组下kNN模型中k的最佳取值，并对其对应取值下的验证集准确率进行统计。统计结果显示，未经过PCA选择特征的kNN分类平均准确率为87.75%，经过PCA选择特征的kNN分类平均准确率为89.23%。

为了观察出k的取值对分类准确率的影响。对不同分组下同一k值的验证集准确率进行统计并求平均值，得到的变化规律曲线如图 6中PCA-kNN(a)和kNN(a)所示。从图 6中可以看出，当k的取值为9时，对应的验证集平均准确率最高，即构建模型最优，是该组数据的最佳k值。

图  6  验证集平均准确率随k取值的变化规律

Figure  6.  Variation of the average accuracy of the test set with the value of k

下载: 全尺寸图片 幻灯片

另外，还增加了一组对比实验。在该次实验中，对已知数据集中选取训练集和测试集的比例做出调整，每次选择3组作为训练集，其余7组作为测试集，被选取为训练集的小组依次为：1、2、3组，2、3、4组，⋯，10、1、2组，共进行了10次选取测试。对该对比实验的结果进行分析可知，没有经过PCA特征选择时，各分组下最佳k值对应的分类平均准确率为73.79%；经过PCA特征选择时，分类平均准确率为71.63%。同时，将测试得到的同一k值的验证集准确率求平均，结果如图 6中曲线PCA-kNN(b)和kNN(b)所示。对比前面的准确率计算结果和曲线a、b可以看出，训练样本数量的减少，使得构建的模型泛化能力降低，模型预测的准确率明显下降，这一方面反映了测井和巷道提供已知信息的重要性，训练样本越多，构建的模型预测效果越好。

3.3   断层分布预测

利用交叉验证法确定了PCA-kNN的最优参数取值后，对羊东煤矿2号煤层的断层分布进行了预测，预测结果如图 7所示。

图  7  kNN断层预测分布与人工解释结果对比

Figure  7.  Comparison of kNN fault prediction distribution and artificial interpretation results

下载: 全尺寸图片 幻灯片

从图 7中可以看出：PCA-kNN预测的断层分布与钻井、巷道揭露的断层分布吻合程度较高，断层走向基本一致，这表明利用PCA-kNN的方法可以实现煤矿断层的分布预测；通过对比可以发现，使用kNN模型获得的断层分布显然比人工断层解释结果更为宽松，这是由于PCA-kNN利用综合属性预测断层区域时，它仅对综合属性值的异常区域产生断层响应，因此，预测结果是相对松散分布的断层点；图中的A处和F处均有断层揭露，人工解释的结果却并未解释出这2个区域的断层，这是因为该处的断层落差较小，地震解释剖面上的断层特征不明显，人工解释常常难以解释出这种落差较小的断层，但是PCA-kNN将这2个区域均预测为断层区，这表明PCA-kNN模型相较于人通过肉眼来解释响应更为灵敏，在小断层的预测方面具有一定的优势；同理，可以推测B、C两处，PCA-kNN的预测结果与人工解释结果相反，极有可能是该区域存在小断层分布，但是人工并未解释出这些断层；另外，图中的D、E处被人为地解释为断层分布区，但是PCA-kNN的预测结果显示这些区域的断层扩展长度很小，这表明这些区域断层解释的可靠性相对较低。总体来看，PCA-kNN可以实现工区的断层分布预测，并且相较于人工断层解释而言，这一方法具有快速、直观，能够更好地识别小断层的特点。

4   结论

a. 本文在开展研究过程中，随机形成了2种类型的数据集：其中，数据集类型1的训练数据集与测试数据集的比率为9∶1，用于表示训练数据比测试数据多的情况；数据集类型2的训练数据集与测试数据集的比率为3∶7，用于表示训练数据比测试数据少的情况。基于以上2种数据集，分析了kNN和PCA-kNN的自动断层识别的准确率，并在峰峰矿区羊东煤矿进行了基于PCA-kNN算法的断层自动识别。

b. 基于kNN和PCA-kNN，数据集类型1的断层识别精度高于数据集类型2的断层识别精度。这表明断层识别的精度与数据集数量密切相关。当训练数据集的数量大于测试数据集的数量时，断层识别的准确性更高。

c. 基于数据集类型1和数据集类型2，分别开展了kNN和PCA-kNN方法的断层识别。结果表明，基于PCA-kNN模型的断层识别准确率要高于单纯基于kNN模型的断层识别准确率。特别是在数据集类型2中，当训练数据集的数量少于测试数据集的数量时，基于PCA-kNN模型的断层识别精度超过了基于kNN模型的断层识别精度。这表明由PCA形成的单个综合属性比单个地震属性具有更高的有效信息密度，并且具有进一步探查可以表征断层的地震属性的能力。

d. 在实际应用中，使用所构建的模型来预测实际的断层分布。这种情况类似于本研究中的数据集类型2，即训练数据小于测试数据。为了获得比测试数据更多的训练数据，有必要在勘探区域中整合多个已知的地质数据，以进行实际的大数据分析，从而有助于形成与数据集类型1类似的情况。如果发生数据集类型2，可以考虑采用降维方法(例如PCA)来提高断层自动识别的准确性。