CN110516930B - 一种基于Voronoi图的煤层稳定性定量评价方法 - Google Patents

Abstract

本发明公开了一种基于Voronoi图的煤层稳定性定量评价方法，包括以下步骤，S1，确定评价煤层；S2，确定评价区范围；S3，确定评价区煤厚点集；S4，根据离散分布的煤厚点集生成Delaunay三角网；S5，煤厚点点集凸包获取；S6，求取凸包绑定框(C)、评价区边界绑定框(D)及最大绑定框(M)；S7，利用凸包内的点集构建Voronoi多边形；S8，利用凸包上的点集构建Voronoi多边形；S9，用评价区边界裁剪最大绑定框内的Voronoi多边形；S10，计算评价区内每个多边形的面积；S11，确定评价区内的加权平均煤厚；S12，确定基于权值的煤层稳定性评价的其他指标；本发明能更好地评价煤层稳定性，评价过程客观、准确、高效、直观，能够提高煤炭资源评价中信息化、自动化和智能化水平。

Description

一种基于Voronoi图的煤层稳定性定量评价方法

技术领域

本发明涉及煤炭资源地质评价技术领域，特别是一种基于Voronoi图的煤层稳定性定量评价方法。

背景技术

由于煤层厚度、煤层平均厚度及其变化是煤层厚度分级和煤层稳定性评价的重要指标，能否科学、合理、客观地确定煤层厚度及其变化情况，对于客观、准确评价煤层稳定性具有重要意义，目前为止，煤层稳定性定量评价的方法有数理统计法、地质统计法、模糊综合评判法、趋势剩余变异法等。

数理统计法其实质是一种经典的统计法；用它评价煤层稳定性的优点在于煤层厚度平均值、可采性指数反映煤层聚集情况和经济价值；煤层厚度的标准离差反映煤层的绝对误差，而煤层厚度的变异系数则表示煤层单位厚度的波动情况，平均煤层厚度特征数作为划分煤层稳定性和选择勘探网的主要定量数据。但该方法对煤层厚度的空间分布情况无法表示，它把煤层厚度在各点的取值当作随机变量进行处理，忽略了各点取值的相关性及空间分布的特点，因此，在煤层厚度数据点值相同，而煤层厚度的空间分布极不相同的煤层，却得到相同的变异系数。

利用地质统计学，针对反映煤层稳定性的不同变量，构造出其空间分布的理论模型用于定量评价煤层稳定性。显然本方法考虑了煤层厚度区域性变化及其空间分布。但是，地质统计法确定构建半变差函数时的数据选取时，需要考虑煤层厚度的取值范围、煤厚数据分布的角度允许误差限和距离允许误差限，以及为了保证r(h)值的代表性而要求的最低数据配对个数，这些有人为因素参与的活动难免存在主观性；同时建立参数指标时的过程复杂，研究过程太繁琐，不易为大多数人掌握，适合于煤层厚度变化不大情况下使用。

煤层稳定性状况具有一系列的模糊特征，为采用模糊数学方法评价和预测研究提供了可能。目前在煤层稳定性评价过程中，主要采用模糊综合评判法和模式识别，配合以模糊聚类和灰色关联分析等方法。根据煤层稳定性程度具有模糊性的特点，因此，利用模糊综合评判理论，将反映煤层稳定性的指标煤层平均厚度、变异系数和煤层结构复杂程度等有机的结合起来，建立煤层稳定性评价模型，从而，更全面可靠地确定勘查区的煤层稳定性情况。但是此方法存在不足表现在公式繁琐，评价参数不确定，并且现有评价中没有选择代表空间分布特征的指标，不能反映煤层厚度的空间分布。

趋势剩余变异法是结合趋势面法和数理统计法两种方法提出的一种煤层稳定性定量评价方法；其方法的计算步骤为先利用数理统计的方法，算出煤层厚度的平均值、变异系数和可采性指数；再做煤层厚度趋势面分析和剩余值分析，计算剩余值的标准差，求出趋势剩余变异系数，作相应的趋势图和偏差图。这种评价方法考虑了煤厚数据点的空间分布情况，但计算煤厚趋势剩余变异系数基础之一煤厚平均值，仍然是一个算术平均值，没有考虑各煤厚点的空间相关性，同时其他指标仍然采用数理统计法中的相应指标，从而影响评价结果的合理性、准确性；另外利用模糊综合评价方法评价过程及其复杂。

因此，数理统计法、模糊综合评判法、趋势剩余变异法均涉及到利用简单统计意义上的描述指标：煤层平均厚度、可采性指数、煤厚变异系数，因此在进行具体评价时，均存在评价不够客观、可靠的问题。

发明内容

为了克服上述不足，本发明的目的是要提供一种基于Voronoi图的煤层稳定性定量评价方法，在评价煤层稳定性时，基于地理学第一定律，考虑各个煤厚点及其评价范围的相关性，从而使得评价结果客观、可靠；同时以地学信息化技术为支撑，实现评价过程的信息化，自动化。

为达到上述目的，本发明是按照以下技术方案实施的：

一种基于Voronoi图的煤层稳定性定量评价方法，包括以下步骤：

S1，确定评价煤层；S2，确定评价区范围：根据评价区拐点序列确定，形成多边形区域；S3，确定评价区煤厚点集：利用分布在评价区内及评价区边界外围的钻探、地质调查、山地工程、地球物理勘探等手段确定可靠的煤厚点集；S4，根据离散分布的煤厚点集生成Delaunay三角网及计算各三角形的外心；S5，煤厚点点集凸包获取：利用格雷厄姆扫描算法，从煤厚点集中获取凸包，并将凸包上的点按逆时钟排列；S6，求取凸包绑定框(C)、评价区边界绑定框(D)及最大绑定框(M)，M＝CUD；S7，利用凸包内的点集构建Voronoi多边形；S8，利用凸包上的点集构建Voronoi多边形；S9，用评价区边界裁剪最大绑定框内的Voronoi多边形；S10，计算评价区内每个多边形的面积；S11，确定评价区内的加权平均煤厚；S12，确定基于权值的煤层稳定性评价的其他指标；

具体的，S4中，采用三角网生长算法来构建Delaunay三角网；

具体的，S4中，计算各三角形外心：已知三角形顶点p₁(x₁,y₁)、p₂(x₂,y₂)、p₃(x₃,y₃)，则该三角形外心的坐标为：

具体的，S7中，利用凸包内的点集构建Voronoi多边形的步骤包括：步骤1，创建存储该多边形的数据结构；步骤2，查找与该点共点的所有三角形，将共点三角形的外心加入多边形顶点集；步骤3，将多边形顶点集按逆时钟排序；

具体的，S8中，利用凸包上的点集构建Voronoi多边形的步骤包括：步骤A，创建存储该多边形的数据结构；步骤B，查找与该点共点的所有三角形；若共点三角形的外心位于最大绑定框之上，则将此外心加入多边形顶点集；若共点三角形的外心位于最大绑定框之内：(a)则将此外心加入多边形顶点集；(b)计算凸包每个点与其前、后邻点连线的两条垂直平分线，并计算此二线与最大绑定框的两个交点(由于凸包上的点集按逆时钟方向排列，各邻点连线为逆时钟方向的有向线，所求有向线与最大绑定框的交点是位于该有向线右侧的交点)，然后将此二交点及被分割后相关的弧段上的结点加入多边形顶点集；若共点三角形的外心位于最大绑定框之外：则查找此三角形中顶点位于凸包内的顶点所关联的多边形的两条边，这两条边与最大绑定框相交，并将其中与此共点距离最近交点加入多边形顶点集；同时将被分割后相关的弧段上的结点加入多边形顶点集；步骤C，将多边形顶点集按逆时钟排序；

具体的，S9中，用评价区边界裁剪最大绑定框内Voronoi多边形包括：(1)遍历查找与评价区边界相交的Voronoi多边形；(2)对于每个与评价区边界相交的Voronoi多边形进行求交运算；

具体的，S10中，已知任意多边形具有n个顶点，其中第i个顶点的坐标为(x_i,y_i),i＝1,2…,n，则该多边形的面积为：

上式中的计算结果可能为正，也可能为负，所以要取其绝对值；即最终面积为：A＝Abs(A)；

具体的，S11中，由于每个Voronoi多边形内有且仅有一个煤厚点，设评价区共划分为n个多边形，第i个多边形的煤厚值为M_i，面积为A_i，则评价区的加权平均煤厚为：

具体的，S12中，所述基于权值的煤层稳定性评价的其他指标包括加权煤层可采性指数、加权煤厚变异系数，所述加权煤层可采性指数K_wm为：

其中，n′为煤层厚度大于或等于可采厚度的煤厚点数，n表示参与评价的煤厚点数；A_i为第i个煤层厚度大于或等于可采厚度的煤厚点所在多边形面积，A_j为第j个煤厚点所在多边形面积；

所述加权煤厚变异系数为；

式中，γ_w——加权煤厚变异系数，M_i——每个煤厚点煤层厚度，

——煤矿(或分区)的加权平均煤层厚，S_w——加权煤厚均方差值。

与现有技术相比，本发明基于计算机及信息科学技术，针对煤炭资源评价中煤层稳定性评价现行方法存在问题，提出了基于Voronoi图的煤层稳定性定量评价方法，能更好地评价煤层稳定性，评价过程客观、准确、高效、直观，能够提高煤炭资源评价中信息化、自动化和智能化水平；主要优点如下：

(1)评价区煤层稳定性评价的重要指标——平均煤层厚度的计算中，采用了加权平均煤层厚度

提高了该指标的客观性、可靠性。现行方法中，平均煤层厚度的计算采用简单的算术平均值，没有考虑煤厚点的空间分布与每个煤厚点相关空间范围；本方法根据最大相关性原则，将空间离散分布的煤厚点相关空间区域按照Voronoi多边形剖分，以严格的数学理论为基础，保证各煤厚点在其相关多边形中的唯一性，以及煤厚点与其相关多边形内任意点的距离最近，从而使煤厚点能够代表相关多边形区域内任意点煤厚。因此，以评价区内各个多边形的面积作为权值计算评价区加权平均煤厚，对于表征评价区煤层的总体厚度来说显得更加客观、真实。

(2)以多边形面积为权值的加权煤层可采性指数K_wm，可以更客观、真实地反映评价区煤层的可采情况。现行方法中，可采性指数采用评价区内达到或超过最低可采厚度的煤厚点数的百分数来表示。由于煤厚点在空间上分布位置、分布密度不同，使得不同煤厚点能表征的空间范围存在差异，简单的百分数表示法则忽略了这种差异。采用以多边形面积为权值的加权煤层可采性指数，能更客观、真实反映评价区的煤层的可采性。

(3)基于加权平均煤厚的加权煤厚变异系数，能够更合理地表达评价区内煤层厚度的变化情况。以加权平均煤厚为基础的加权煤厚变异系数。

(4)评价过程的高度自动化、信息化、可视化。相比以往评价方法，本方法只需要将评价区边界点集、煤厚点集、最低可采厚度数据按一定格式存储在数据库或Excel表中，评价过程会自动按算法实现煤厚点集凸包提取、Delaunay三角剖分、以评价区边界为约束条件的Voronoi图生成、多边形面积计算、评价指标的生成及评价结果的展示等全过程。

(5)评价结果的唯一性、不变性。对于特定的数据集，本方法中采用的基于Delaunay三角剖分的Voronoi图生成方法(间接法)，从算法上保证了Delaunay三角网及其对偶Voronoi图生成具有唯一性。即重复本方法，生成的与各煤厚点相关的多边形都是唯一的、具有几何不变性，与生成的时间因素无关。

(6)评价区Voronoi图维护的高效。对于需要增加或删除特定煤厚点的情况，本方法仅需要从算法上对增、删部位的Voronoi图进行局部调整重构，不会影响其他部位，所以维护容易、高效。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的凸包上点的Voronoi多边形顶点图；

图2为本发明的凸包内点关联的多边形与最大绑定框相交图；

图3为本发明的交点凹凸性判断图；

图4为本发明的交点凹凸性判断图；

图5为本发明交点大于1时评价区边界多边形图；

图6为本发明交点为1时除交点之外两多边形没有公共区域图；

图7为本发明交点为1时一个多边形被另一多边形包围图；

图8为本发明实施例山东某井田范围及其Voronoi剖分图；

图9为本发明实施例山东某井田评价区边界裁剪后的内部Voronoi图。

具体实施方式

下面结合附图以及具体实施例对本发明作进一步描述，在此发明的示意性实施例以及说明用来解释本发明，但并不作为对本发明的限定。

一种基于Voronoi图的煤层稳定性定量评价方法，包括以下步骤：

S1，确定评价煤层。

S2，确定评价区范围：根据评价区拐点序列确定，形成多边形区域。

S3，确定评价区煤厚点集：利用分布在评价区内及评价区边界外围的钻探、地质调查、山地工程、地球物理勘探等手段确定可靠的煤厚点集。

S4，根据离散分布的煤厚点集生成Delaunay三角网及计算各三角形的外心；由于一般评价区煤厚点数目有限，采用三角网生长算法构建。三角网生长算法逻辑简单，易于编程实现，主要操作是在点集中搜寻给定基线符合要求的邻域点。本算法在最坏情况下的时间复杂度为O(n²)。存储各个三角形及其构成此三角形的边、顶点；计算各三角形外心：已知三角形顶点p₁(x₁,y₁)、p₂(x₂,y₂)、p₃(x₃,y₃)，则该三角形外心的坐标为：

S5，煤厚点点集凸包获取：利用格雷厄姆扫描(Graham Scan)算法，从煤厚点集中获取凸包，并将凸包上的点按逆时钟排列。

S6，求取凸包绑定框(C)、评价区边界绑定框(D)及最大绑定框(M)，M＝CUD。

S7，利用凸包内的点集构建Voronoi多边形；包括：步骤1，创建存储该多边形的数据结构；步骤2，查找与该点共点的所有三角形，将共点三角形的外心加入多边形顶点集；步骤3，将多边形顶点集按逆时钟排序。

S8，利用凸包上的点集构建Voronoi多边形；包括：步骤A，创建存储该多边形的数据结构；步骤B，查找与该点共点的所有三角形；若共点三角形的外心位于最大绑定框之上，则将此外心加入多边形顶点集；如图1所示，与M、N共点的三角形MNQ的外心K位于最大绑定框上，则将K点分别加入与M、N关联的两个多边形顶点集中；若共点三角形的外心位于最大绑定框之内：(a)则将此外心加入多边形顶点集；(b)计算凸包每个点与其前、后邻点连线的两条垂直平分线，并计算此二线与最大绑定框的两个交点(由于凸包上的点集按逆时钟方向排列，各邻点连线为逆时钟方向的有向线，所求有向线与最大绑定框的交点是位于该有向线右侧的交点)，然后将此二交点及被分割后相关的弧段上的结点加入多边形顶点集；如图1所示，对于凸包上的点G，计算其前、后点H、F连线的垂直平分线VT、WS与最大绑定框的交点T、S，T和S点分别位于有向线GH、FG的右侧，同时T、S点在最大绑定框分割出弧段SUT，其中该弧段结点U也位于有向线GH、FG的右侧，将点S、T及U加入到G点关联的多边形顶点集中；若共点三角形的外心位于最大绑定框之外：则查找此三角形中顶点位于凸包内的顶点所关联的多边形的两条边，这两条边与最大绑定框相交，并将其中与此共点距离最近交点加入多边形顶点集；同时将被分割后相关的弧段上的结点加入多边形顶点集；如图2所示，对于凸包上的点E，其共点三角形EFP的外心C位于最大绑定框外，该三角形的顶点P位于凸包内，P所关联的多边形有两条边分别与最大绑定框相交于A、B两点，将与E点最近的边与最大绑定框的交点A加入与E关联的多边形顶点集中；同时将相关的弧段上的结点D加入多边形顶点集；步骤C，将多边形顶点集按逆时钟排序。

S9，用评价区边界裁剪最大绑定框内的Voronoi多边形；实质是评价区边界任意多边形(凹或凸)与Voronoi凸多边形的求交问题。解决方法：基于交点的凹凸性实现任意多边形求交运算，其优点是算法简单，适用于任意多边形。前提条件：Voronoi凸多边形上顶点、评价区边界任意多边形顶点均按逆时钟排列，这样多边形上的边即为有向边(矢量直线)。矢量直线与矢量直线交点的凹凸性：如图3，矢量

与

的交点为O，设由A出发，

沿逆时钟方向旋转到

的旋转角α，若α≤180°，则O为凸点，如图3；若180°<α≤360°，则O为凹点，如图4；具体实现时，旋转角α的取值范围确定可通过

来实现。用评价区边界裁剪最大绑定框内Voronoi多边形包括：(1)遍历查找与评价区边界相交的Voronoi多边形；(2)对于每个与评价区边界相交的Voronoi多边形进行求交运算；得到的多边形可能是凸多边形，也可能是凹多边形；若交点个数大于1(如图5)，设P为评价区边界多边形，Q为与边界相交的Voronoi多边形。从P中的某一交点出发，根据P中以该交点为终点的矢量与多边形Q中以该交点为起点的矢量判定该交点的凹凸性。若为凸点，则转向沿Q多边形遍历；若为凹点，则继续沿P多边形遍历。如此，每遇到一个凸点，转向另一个多边形，否则继续遍历当前多边形，直至回到出发点。若尚有未曾遍历的交点，则重复上述过程，直至遍历所有交点。如此遍历路径上的点集形成的多边形即为两多边形的交集。

若交点个数为1(如图6)，这时有两种情况：①除交点之外两多边形没有公共区域；②一个多边形被另一多边形包围。这两种情况可通过交点计数检验法来判别：通过一个多边形P、Q上的任一异于交点的顶点向任意方向作射线(限制条件：射线不得通过另一多边形的端点或切点。若通过时，则必须调整射线方向，避开端点或切点)，若两条射线与另一多边形的交点个数均为偶数(含0)，则P、Q的关系属于情况①，此时P∩Q即为交点O，如图6；若一奇一偶，则P、Q的关系属于情况②，此时P∩Q即为奇数交点射线的发出多边形Q，如图7。

S10，计算评价区内每个多边形的面积；已知任意多边形具有n个顶点，其中第i个顶点的坐标为(x_i,y_i),i＝1,2…,n，则该多边形的面积为：

上式中的计算结果可能为正，也可能为负，所以要取其绝对值；即最终面积为：A＝Abs(A)。

S11，确定评价区内的加权平均煤厚；由于每个Voronoi多边形内有且仅有一个煤厚点，设评价区共划分为n个多边形，第i个多边形的煤厚值为M_i，面积为A_i，则评价区的加权平均煤厚为：

S12，确定基于权值的煤层稳定性评价的其他指标；所述基于权值的煤层稳定性评价的其他指标包括加权煤层可采性指数、加权煤厚变异系数，所述加权煤层可采性指数K_wm为：

其中，n′为煤层厚度大于或等于可采厚度的煤厚点数，n表示参与评价的煤厚点数；A_i为第i个煤层厚度大于或等于可采厚度的煤厚点所在多边形面积，A_j为第j个煤厚点所在多边形面积；

所述加权煤厚变异系数为；

式中，γ_w——加权煤厚变异系数，M_i——每个煤厚点煤层厚度，

——煤矿(或分区)的加权平均煤层厚，S_w——加权煤厚均方差值。

山东某井田开采的主要煤层为3_上煤，评价区范围如图8所示，参与煤厚点共51个，煤厚0.38～3.55m，属中厚煤层。煤厚数据如下表示。

通过本算法先获取凸包，凸包内进行Delaunay三角剖分，并构建Voronoi多边形，求取评价区边界绑定框与凸包绑定框并集，得到最大绑定框。位于凸包上的煤厚点的Voronoi多边形需要通过增加由凸包上的相应边的垂直平分线与最大绑定框的交点来构建。最终得到最大绑定框内的Voronoi图如图8中区域构成。经评价区边界裁剪后的内部Voronoi图如图9所示。根据本方法计算得到：平均煤层厚度：1.72m；加权平均煤层厚度：1.60m；可采性指数：0.84；加权可采性指数：0.82；变异系数：27％；加权变异系数：30％；加权方法煤层稳定性(主要指标)：较稳定；加权方法煤层稳定性(辅助指标)：较稳定。

本方法得到的加权平均煤层厚度为1.60m，比算术平均法得到的结果小0.12m；加权可采性指数为0.82，比原方法的结果小0.02；加权变异系数为30％，比原方法结果大3％。

结合加权方法煤层稳定性评价的主要指标、辅助指标的评价结果，本评价区煤层的稳定性为较稳定煤层，符合实际情况。评价过程快速、高效、自动化水平高。

本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。

Claims (7)

1.一种基于Voronoi图的煤层稳定性定量评价方法，其特征在于，包括以下步骤：

S1，确定评价煤层；

S2，确定评价区范围：根据评价区拐点序列确定，形成多边形区域；

S3，确定评价区煤厚点集：利用分布在评价区内及评价区边界外围的钻探、地质调查、山地工程、地球物理勘探手段确定可靠的煤厚点集；

S4，根据离散分布的煤厚点集生成Delaunay三角网及计算各三角形的外心；

S5，煤厚点点集凸包获取：利用格雷厄姆扫描算法，从煤厚点集中获取凸包，并将凸包上的点按逆时钟排列；

S6，求取凸包绑定框C、评价区边界绑定框D及最大绑定框M，M＝C∪D；

S7，利用凸包内的点集构建Voronoi多边形；

S8，利用凸包上的点集构建Voronoi多边形；

S9，用评价区边界裁剪最大绑定框内的Voronoi多边形；

S10，计算评价区内每个多边形的面积；

S11，确定评价区内的加权平均煤厚，由于每个Voronoi多边形内有且仅有一个煤厚点，设评价区共划分为n个多边形，第i个多边形的煤厚值为M_i，面积为A_i，则评价区的加权平均煤厚为：

S12，确定基于权值的煤层稳定性评价的其他指标，所述基于权值的煤层稳定性评价的其他指标包括加权煤层可采性指数、加权煤厚变异系数，所述加权煤层可采性指数K_wm为：

其中，n′为煤层厚度大于或等于可采厚度的煤厚点数，n表示参与评价的煤厚点数；A_i为第i个煤层厚度大于或等于可采厚度的煤厚点所在多边形面积，A_j为第j个煤厚点所在多边形面积；

所述加权煤厚变异系数为；

式中，γ_w——加权煤厚变异系数，M_i——每个煤厚点煤层厚度，

——煤矿的加权平均煤层厚，S_w——加权煤厚均方差值。

2.根据权利要求1所述的基于Voronoi图的煤层稳定性定量评价方法，其特征在于，S4中，采用三角网生长算法来构建Delaunay三角网。

3.根据权利要求2所述的基于Voronoi图的煤层稳定性定量评价方法，其特征在于，S4中，计算各三角形外心：已知三角形顶点p₁(x₁,y₁)、p₂(x₂,y₂)、p₃(x₃,y₃)，则该三角形外心的坐标为：

4.根据权利要求3所述的基于Voronoi图的煤层稳定性定量评价方法，其特征在于，S7中，利用凸包内的点集构建Voronoi多边形的步骤包括：步骤1，创建存储该多边形的数据结构；步骤2，查找与该点共点的所有三角形，将共点三角形的外心加入多边形顶点集；步骤3，将多边形顶点集按逆时钟排序。

5.根据权利要求4所述的基于Voronoi图的煤层稳定性定量评价方法，其特征在于，S8中，利用凸包上的点集构建Voronoi多边形的步骤包括：步骤A，创建存储该多边形的数据结构；步骤B，查找与该点共点的所有三角形；步骤C，将多边形顶点集按逆时钟排序。

6.根据权利要求5所述的基于Voronoi图的煤层稳定性定量评价方法，其特征在于，S9中，用评价区边界裁剪最大绑定框内Voronoi多边形包括：(1)遍历查找与评价区边界相交的Voronoi多边形；(2)对于每个与评价区边界相交的Voronoi多边形进行求交运算。

7.根据权利要求6所述的基于Voronoi图的煤层稳定性定量评价方法，其特征在于，S10中，已知任意多边形具有n个顶点，其中第i个顶点的坐标为(x_i,y_i),i＝1,2…,n，则该多边形的面积为：

上式中的计算结果可能为正，也可能为负，所以要取其绝对值；即最终面积为：A＝Abs(A)。

Images