CN108805154A - 一种基于空间聚类的地质断层识别方法 - Google Patents

Abstract

本发明提供一种基于空间聚类的地质断层识别的方法，涉及空间数据挖掘技术领域。一种基于空间聚类进行地质断层识别的方法，首先对勘测数据点进行三角面剖分，计算三角面的坡度，然后通过三角面之间的坡度变化对这些三角面进行空间聚类，坡度变化的大小反映地质连续的破坏程度，根据这些勘测数据点的所在三角面的分类情况，完成对位于不同断层面的勘测数据点划分，从而实现对断层自动识别。本发明提供的基于空间聚类进行地质断层识别的方法突破了现有断层识别方法存在的缺陷，能够对隐藏在大量的离散勘测数据中的断层进行高效的识别和解释。

Description

一种基于空间聚类的地质断层识别方法

技术领域

本发明涉及空间数据挖掘技术领域，尤其涉及一种基于空间聚类的地质断层识别方法。

背景技术

随着我国数字化和信息化进程的加快，产生了大量的地质勘探数据，同时也出现了“数据越丰富，知识越贫乏”的局面，急需对地质勘探数据深入挖掘进行研究。如何在浩瀚的空间数据和人们的知识渴求之间建立一个桥梁的确是一个巨大的挑战，而以聚类、分类等人工智能技术为基础的空间数据挖掘为迎接这个挑战提供了新的支撑技术。

面对不断积累的海量地质勘探数据，为了缩短解释周期，克服断层解释中的主观性，断层的自动或者半自动识别技术逐渐引起了广大学者和科研人员的关注。早期的断层识别方法采用的数据驱动策略，通过局部连续性约束方法来识别断层位置，自2000年以后，许多断层自动或半自动识别方法陆续被提了出来，其中比较著名的方法有：2001年Randen等人提出的“人工蚂蚁”方法来压制不连续属性体中的噪音；2003年Gibson等人提出使用相干体来衡量不连续性，进而使用确定的断层面模型通过最大可信度优先(highestconfidence first,HCF)提取策略提取断层面；2005年Dom和James使用融合了地质先验知识的信号处理技术实现断层自动提取(AFE)；2006年Won-ki Jeong等基于GPU使用交互式操作进行断层识别。Benjamin J等在2008年提出基于水平集(Level Sets)的交互式断层面计算方法。2010年，严哲博士研究了三维地震断层自动识别与智能解释。这些断层识别的方法均集中在对地震波的断层识别和解释及遥感图像的识别，均不适合应用于大量离散的勘测数据。

发明方法

针对现有技术的缺陷，本发明提供一种基于空间聚类的地质断层识别方法，实现矿山勘测及大地勘测离散空间数据点中自动识别地质断层。

一种基于空间聚类的地质断层识别方法，包括以下步骤：

步骤1：输入勘测数据点集和勘测数据点的坡度变化阈值；

将存储勘测数据点的外部文件读入计算机内存，输入勘测数据点集S＝{P₁,P₂,…,P_i,…,P_n}和勘测数据点的坡度变化阈值α，其中，P_i为空间数据点，i∈[1,n]，n为空间数据点总数，P_i由空间三维坐标(x_i,y_i,z_i)构成，三维坐标(x_i,y_i,z_i)所在的三维坐标系的原点根据实际情况人为设定，三维坐标系的X轴和Y轴表示勘测数据点的水平方向信息，Z轴表示勘测数据点的竖直方向信息；通过编程语言编写的数据结构实现对勘测数据点集的存储、组织和管理；

步骤2：按Delaunay准则，从勘测数据点集S中选取三个数据点，生成第一个三角面，计算新生成的三角面的坡度，为三角面添加分类号，并将新生成的三角面加入到三角面集T中，三角面的分类号加入到分类号表C中；

步骤3：按Delaunay准则生成新的三角面，对三角面集T进行扩展，具体方法为：

步骤3.1：以三角面集T中原有的一个三角面为基础三角面，按Delaunay准则，使用勘测数据点集S中其余的数据点生成新的三角面；

步骤3.2：计算新生成的三角面的坡度，同时计算新生成的三角面的坡度变化值θ，并将坡度变化值θ与勘测数据点的坡度变化阈值α进行比较；如果新生成的三角面的坡度变化值θ大于勘测数据点的坡度变化阈值α，则新生产的三角面不属于三角面集T中的分类，为新生成的三角面添加新的分类号，并将新分类号加入到分类号表C中；如果新生成的三角面的坡度变化值θ小于等于勘测数据点的坡度变化阈值α，则新生产的三角面属于三角面集T中的分类，新生成的三角面的分类号与基础三角面的分类号相同；

步骤3.3：将新生成的三角面加入到三角面集T中，对三角面集T进行扩展；

步骤3.4：判断是否所有的勘测数据点已全部加入到三角面集T中，如果所有的勘测数据点已全部加入到了三角面集T中，则生成三角面集T＝{t₁,t₂,…,t_k,…,t_m}及分类号表C＝{1,2,…,N}，其中，t_k为第k个三角面，m为最终生成的三角面集中三角面的总个数，N为总的分类号，然后执行步骤4，否则返回步骤3.1，以新生成的三角面为基础三角面，继续生成新的三角面，扩展三角面集T；

步骤4：根据勘测数据点所在三角面的分类号，为勘测数据点集S中的点添加分类号，对勘测数据点集S中的数据点进行分类，具体方法为：

勘测数据点集S＝{P₁,P₂,…,P_i,…,P_n}中的每个数据点至少作为三角面集T＝{t₁,t₂,…,t_k,…,t_m}中一个三角面的顶点，其分类由所在三角面的分类情况决定；如果勘测数据点所在的所有三角面的分类号均相同，则将该勘测点的分类号设置为所在三角面的分类号，分类号相同的勘测数据点为一类；如果勘测数据点所在的多个三角面有不同的分类号，则该勘测数据点处在断层分界面上，将该勘测数据点设置为分界点；

步骤5：将聚类后的三角面集T和添加了分类号的勘测数据点集S输出，作为勘探钻孔断层识别结果，并将三角面集T和添加了分类号的勘测数据点集S存储在数据库中。

进一步地，步骤1所述使用编程语言编写的数据结构包括勘测点数据结构、边数据结构和三角面数据结构三类数据结构，每类数据结构包括的具体信息分别如下：

勘测点数据结构，该数据结构中包括勘测点的位置坐标及该点的分类号；

边数据结构，该数据结构包括组成边的顶点的编号、对应的邻接三角面的编号和边所在三角面的编号，其中，邻接三角面的编号描述三角面间的拓扑关系；

三角面数据结构，该数据结构包括三角面的三个顶点、三条边、三角面的梯度、三角面的方向向量和三角面的分类号。

进一步地，步骤2所述新生成的三角面的坡度为该三角面的方向向量与三维坐标系Z轴的正方向的夹角，其计算公式如下：

式中，β_k为新生成的三角面t_k的坡度，V_k为新生成的三角面t_k的方向向量，p_i、p_i-1和p_i-2为新生成的三角面t_k的三个顶点，z为Z轴正方向的单位向量。

进一步地，步骤3.2所述新生成的三角面的坡度变化值θ为新生成的三角面t_k的坡度与基础三角面t_k-1的坡度的差的绝对值，计算公式如下所示：

θ(k,k-1)＝|β_k-β_k-1|

其中，θ(k,k-1)为新生成的三角面t_k与基础三角面t_k-1的坡度变化值，β_k为新生成的三角面t_k的坡度，β_k-1为基础三角面t_k-1的坡度。

由上述技术方案可知，本发明的有益效果在于：本发明提供的一种基于空间聚类的地质断层识别方法，对勘测数据点进行三角面剖分，计算三角面的坡度，通过三角面之间的坡度变化对这些三角面进行空间聚类，坡度变化大小反映地质连续的破坏程度，根据这些勘测数据点的所在三角面的分类情况，完成对位于不同断层面的勘测数据点划分，从而实现对断层自动识别，突破了现有断层识别方法存在的缺陷，能够对隐藏在大量的离散勘测数据中的断层进行高效的识别和解释。

附图说明

图1为本发明实施例提供的一种基于空间聚类的地质断层识别方法的流程图；

图2为本发明实施例提供的勘探钻孔的分布俯视示意图；

图3为本发明实施例提供的勘探钻孔的空间位置分布示意图；

图4为本发明实施例提供的勘测数据点构成的多个三角面的坡度示意图；

图5为本发明实施例提供的勘探钻孔聚类结果的西南等轴测图；

图6为本发明实施例提供的勘探钻孔聚类结果的平面示意图；

图7为本发明实施例提供的离散勘探钻孔数据断层识别结果的立体示意图。

图中：1、第一类勘测数据点；2、第二类勘测数据点；3、第三类勘测数据点；4、第四类勘测数据点。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

本实施例以某煤矿的钻孔数据作为勘测数据的来源，以钻孔的位置坐标(x,y)、煤层的底板标高z作为输入数据，通过一种基于空间聚类的地质断层识别方法实现对钻孔数据点中的地质断层进行识别。

一种基于空间聚类的地质断层识别方法，如图1所示，包括以下步骤：

步骤1：输入勘测数据点集和勘测数据点的坡度变化阈值；

将存储勘测数据点的外部文件读入计算机内存，输入勘测数据点集S＝{P₁,P₂,…,P_i,…,P_n}和勘测数据点的坡度变化阈值α，其中，P_i为空间数据点，i∈[1,n]，n为空间数据点总数，P_i由空间三维坐标(x_i,y_i,z_i)构成；本实施例使用C++编程语言编写的数据结构实现对勘测数据点集的存储、组织和管理；

使用C++编程语言编写的数据结构包括勘测点数据结构、边数据结构和三角面数据结构三类数据结构，每类数据结构包括的具体信息分别如下：

勘测点数据结构，该数据结构中包括勘测点的位置坐标及该点的分类号；

边数据结构，该数据结构包括组成边的顶点的编号、对应的邻接三角面的编号和边所在三角面的编号，其中，邻接三角面的编号描述三角面间的拓扑关系；

三角面数据结构，该数据结构包括三角面的三个顶点、三条边、三角面的梯度、三角面的方向向量和三角面的分类号。

本实施例以勘探钻孔的位置、煤层的底板标高为空间数据点，建立勘探钻孔成果数据库，对空间数据点，进行存储和管理。数据库中存储的244个勘探钻孔的坐标及17#煤的底板高程的部分数据如表1所示，根据现场实际情况设定坡度变化角阈值α为65°。勘探钻孔的分布如图2和图3所示，图2将实施例用到的244个样本勘探钻孔绘制到图形中，图中显示的是钻孔编号、17#煤层的底板标高和钻孔的平面位置，具体数值如图2中局部放大图所示；图3将绘制的勘探钻孔在XYZ三维坐标系中显示，三个坐标轴分别表示X、Y和Z坐标，从图中能直观地看出这些样本钻孔在空间的位置分布。

表1输入的部分数据表

钻孔号	(x，y)	17#底板标高(m)
			26-6	(38528.5，-67898.31)	878.08
25-95-1	(38885.67，-67831.09)	861.87
			…	…	…
25-观6	(39493.84，-68476.62)	878.60

步骤2：按Delaunay准则，从勘测数据点集S中选取三个数据点，生成第一个三角面，计算新生成的三角面的坡度，为新生成的三角面添加分类号，并将新生成的三角面加入到三角面集T中，表示为T＝{t₁}，三角面的分类号加入到分类号表C中，表示为C＝{1}；

新生成的三角面的坡度，如图4所示，为该三角面的方向向量与三维坐标系Z轴的正方向的夹角，其计算公式如下：

式中，β_k为新生成的三角面的坡度，V_k为新生成的三角面t_k的方向向量，p_i、p_i-1和p_i-2为新生成的三角面t_k的三个顶点，z为Z轴正方向的单位向量。

根据三角面的坡度计算公式及已给的空间数据点值，得到本实施例中新生成的第一个三角面的坡度为β₁＝11.24°。

步骤3：按Delaunay准则生成新的三角面，对三角面集T进行扩展，具体方法为：

步骤3.1：以三角面集T中原有的一个三角面为基础三角面，按Delaunay准则，使用勘测数据点集S中其余的数据点生成新的三角面；

步骤3.2：计算新生成的三角面的坡度，同时计算新生成的三角面的坡度变化值θ，并将坡度变化值θ与勘测数据点的坡度变化阈值α进行比较；如果新生成的三角面的坡度变化值θ大于勘测数据点的坡度变化阈值α，则新生产的三角面不属于三角面集T中的分类，为新生成的三角面添加新的分类号，并将新分类号加入到分类号表C中；如果新生成的三角面的坡度变化值θ小于等于勘测数据点的坡度变化阈值α，则新生产的三角面属于三角面集T中的分类，新生成的三角面分类号与基础三角面的分类号相同；

新生成的三角面的坡度变化值θ，如图4所示，为新生成的三角面t_k的坡度与基础三角面t_k-1的坡度的差的绝对值，计算公式如下所示：

θ(k,k-1)＝|β_k-β_k-1|

其中，θ(k,k-1)为新生成的三角面t_k的坡度变化值，β_k为新生成的三角面t_k的坡度，β_k-1为作为新生成三角面生成基础的三角面t_k-1的坡度；

步骤3.3：将新生成的三角面加入到三角面集T，对三角面集T进行扩展；

步骤3.4：判断是否所有的勘测数据点加入到三角面集T中，如果所有的勘测数据点已全部加入到了三角面集T中，则生成三角面集T＝{t₁,t₂,…,t_k,…,t_m}及分类号表C＝{1,2,…,N}，其中，t_k为第k个三角面，m为最终生成的三角面集中三角面的总个数，N为总的分类号，然后执行步骤4，否则返回步骤3.1，以新生成的三角面为基础三角面，继续生成新的三角面，扩展三角面集T；

本实施例以第一个三角面t₁为例对三角面集T进行扩展。以第一个三角面t₁为基础三角面，按Delaunay准则，使用勘测数据点集S中其余数据点对三角面集T进行扩展，得到第二个三角面t₂，使用三角面t₂的顶点数据进行计算，得到第二个三角面的坡度β₂＝13.57°，则新生成的三角面t₂的坡度变化可由θ(2,1)表示，根据坡度变化值计算公式得，θ(2,1)＝|β₂-β₁|＝|13.57°-11.24°|＝2.33°，将新生成的三角面的坡度变化值θ(2,1)与勘测数据的坡度变化阈值65°进行比较，θ(2,1)＝2.33°＜65°，坡度变化值θ(2,1)小于勘测数据的坡度变化阈值，所以新生成的三角面t₂的分类号与基础三角面t₁的分类号相同。

步骤4：根据勘测数据点所在三角面的分类号，为勘测数据点集S中的点添加分类号，对勘测数据点集S中的数据点进行分类，具体方法为：

勘测数据点集S＝{P₁,P₂,…,P_i,…,P_n}中的每个数据点至少作为三角面集T＝{t₁,t₂,…,t_k,…,t_m}中一个三角面的顶点，其分类由所在三角面的分类情况决定；如果勘测数据点所在的所有三角面的分类号均相同，则将该勘测点的分类号设置为所在三角面的分类号，分类号相同的勘测数据点为一类；如果勘测数据点所在的多个三角面有不同的分类号，则该勘测数据点处在断层分界面上，将该勘测数据点设置为分界点；

本实施例中勘测数据点集S的聚类结果如表2所示，得到第一类数据点、第二类数据点、第三类数据点和第四类数据点共四类，每类数据点的边界为断层分界面，作为分界点的勘测数据点处在该分界面上；将这些勘测数据点的聚类结果用如图5和图6表示，图中不同的分类点用不同形状进行表示，第一类勘测数据点用圆形表示，第二类勘测数据点用十字形表示，第三类勘测数据点用三角形表示，第四类勘测数据点用正方形表示。

表2聚类结果部分数据表

钻孔号	(x，y)	17#底板标高(m)	分类号
				26-6	(38528.5，-67898.31)	878.08	1
25-95-1	(38885.67，-67831.09)	861.87	1
				26-5	(38376.35，-68097.94)	898.18	4
33-6	(36690.60，-66871.80)	742.76	3
				…	…	…	…
25-观6	(38000.03，-68799.89)	879.84	2

步骤5：将聚类后的三角面集T和添加了分类号的勘测数据点集S输出，作为勘探钻孔断层识别结果，并将三角面集和添加了分类号的勘测数据点集S存储在数据库中。

本实施例除了进行数据库存储输出以外，还进行了三维可视化输出，如图7所示，作为对离散勘探钻孔数据断层识别结果的立体显示，从图中可以看出虚线表示的两个断层面将勘测数据点分为四类，分别为第一类勘测数据点1、第二类勘测数据点2、第三类勘测数据点3和第四类勘测数据点4。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (4)

1.一种基于空间聚类的地质断层识别方法，其特征在于：包括以下步骤：

步骤1：输入勘测数据点集和勘测数据点的坡度变化阈值；

将存储勘测数据点的外部文件读入计算机内存，输入勘测数据点集S＝{P₁,P₂,…,P_i,…,P_n}和勘测数据点的坡度变化阈值α，其中，P_i为空间数据点，i∈[1,n]，n为空间数据点总数，P_i由空间三维坐标(x_i,y_i,z_i)构成，三维坐标(x_i,y_i,z_i)所在的三维坐标系的原点根据实际情况人为设定，三维坐标系的X轴和Y轴表示勘测数据点的水平方向信息，Z轴表示勘测数据点的竖直方向信息；通过编程语言编写的数据结构实现对勘测数据点集的存储、组织和管理；

步骤2：按Delaunay准则，从勘测数据点集S中选取三个数据点，生成第一个三角面，计算新生成的三角面的坡度，为三角面添加分类号，并将新生成的三角面加入到三角面集T中，三角面的分类号加入到分类号表C中；

步骤3：按Delaunay准则生成新的三角面，对三角面集T进行扩展，具体方法为：

步骤3.1：以三角面集T中原有的一个三角面为基础三角面，按Delaunay准则，使用勘测数据点集S中其余的数据点生成新的三角面；

步骤3.2：计算新生成的三角面的坡度，同时计算新生成的三角面的坡度变化值θ，并将坡度变化值θ与勘测数据点的坡度变化阈值α进行比较；如果新生成的三角面的坡度变化值θ大于勘测数据点的坡度变化阈值α，则新生产的三角面不属于三角面集T中的分类，为新生成的三角面添加新的分类号，并将新分类号加入到分类号表C中；如果新生成的三角面的坡度变化值θ小于等于勘测数据点的坡度变化阈值α，则新生产的三角面属于三角面集T中的分类，新生成的三角面的分类号与基础三角面的分类号相同；

步骤3.3：将新生成的三角面加入到三角面集T中，对三角面集T进行扩展；

步骤3.4：判断是否所有的勘测数据点已全部加入到三角面集T中，如果所有的勘测数据点已全部加入到了三角面集T中，则生成三角面集T＝{t₁,t₂,…,t_k,…,t_m}及分类号表C＝{1,2,…,N}，其中，t_k为第k个三角面，m为最终生成的三角面集中三角面的总个数，N为总的分类号，然后执行步骤4，否则返回步骤3.1，以新生成的三角面为基础三角面，继续生成新的三角面，扩展三角面集T；

步骤4：根据勘测数据点所在三角面的分类号，为勘测数据点集S中的点添加分类号，对勘测数据点集S中的数据点进行分类，具体方法为：

勘测数据点集S＝{P₁,P₂,…,P_i,…,P_n}中的每个数据点至少作为三角面集T＝{t₁,t₂,…,t_k,…,t_m}中一个三角面的顶点，其分类由所在三角面的分类情况决定；如果勘测数据点所在的所有三角面的分类号均相同，则将该勘测点的分类号设置为所在三角面的分类号，分类号相同的勘测数据点为一类；如果勘测数据点所在的多个三角面有不同的分类号，则该勘测数据点处在断层分界面上，将该勘测数据点设置为分界点；

步骤5：将聚类后的三角面集T和添加了分类号的勘测数据点集S输出，作为勘探钻孔断层识别结果，并将三角面集T和添加了分类号的勘测数据点集S存储在数据库中。

2.根据权利要求1所述的一种基于空间聚类的地质断层识别方法，其特征在于：步骤1所述使用编程语言编写的数据结构包括勘测点数据结构、边数据结构和三角面数据结构三类数据结构，每类数据结构包括的具体信息分别如下：

勘测点数据结构，该数据结构中包括勘测点的位置坐标及该点的分类号；

边数据结构，该数据结构包括组成边的顶点的编号、对应的邻接三角面的编号和边所在三角面的编号，其中，邻接三角面的编号描述三角面间的拓扑关系；

三角面数据结构，该数据结构包括三角面的三个顶点、三条边、三角面的梯度、三角面的方向向量和三角面的分类号。

3.根据权利要求1所述的一种基于空间聚类的地质断层识别方法，其特征在于：步骤2所述新生成的三角面的坡度为该三角面的方向向量与三维坐标系Z轴的正方向的夹角，其计算公式如下：

式中，β_k为新生成的三角面t_k的坡度，V_k为新生成的三角面t_k的方向向量，p_i、p_i-1和p_i-2为新生成的三角面t_k的三个顶点，z为Z轴正方向的单位向量。

4.根据权利要求3所述的一种基于空间聚类的地质断层识别方法，其特征在于：步骤3.2所述新生成的三角面的坡度变化值θ为新生成的三角面t_k的坡度与基础三角面t_k-1的坡度的差的绝对值，计算公式如下所示：

θ(k,k-1)＝|β_k-β_k-1|

其中，θ(k,k-1)为新生成的三角面t_k与基础三角面t_k-1的坡度变化值，β_k为新生成的三角面t_k的坡度，β_k-1为基础三角面t_k-1的坡度。