集团总工程师
正高级工程师
郝治朝
(同煤集团地质勘测处,山西,大同,037003)
摘要:本文引入空间插值方法将回采区域格网化,并以此格网为基础,结合生产实际揭露的煤层底板标高,借助数学分析方法,预测回采工作面前方的地质构造。
关键字: 空间插值 格网 构造 预测
由于矿井小断层构造发育规律认识不清,探查手段落后,构造预测准确性较差,给煤矿生产和建设带来了很多困难。同时,因对小断层构造认识不清而导致底板突水,甚至淹没矿井的现象也时有发生。因此,必须对矿井小断层构造井进行探查和预测,并尽量提高预测的准确性。利用数学地质中统计分析的方法,通过对已被揭露的小断层构造的分析研究,经数学计算,建立一种数学模型,用它来拟合矿井小断层构造的出现和展布规律,可以达到对矿山小断层构造定量的认识。
随着GIS技术的不断发展以及研究中对空间高质量数据的要求,空间数据插值应用越来越广。通常大家通过曲面拟合来描述煤层底板,这在实际应用中受到各种地质条件的限制,本文采用了更为灵活的格网插值方式,对于断层影响下的不连续曲面达到了更符合客观的拟合。以此为基础,可以借助数据分析方法更为科学地预测回采工作面前方的地质构造。
所谓三角形网格化也就是对平面点集的计算其三角剖分,具体实现过程是将平面点集各点用互不相交的直线段连接起来,构成一个三角形网格。在这种情况下,三角网中各三角形顶点仍为原始数据点。在追踪等值线时,以其中的每个三角形为单元,进行等值点的计算及等值线的追踪。
一般而言,三角形网格构成方式有多种,不同的三角化方式有不同的标准。目前被普遍接受和广泛采用的三角剖分有最小权三角剖分、贪婪三角剖分和Delaunay三角剖分。最小权三角剖分是指该离散点集可以够成的所有三角剖分中的边总长最小的三角剖分;而贪婪三角剖分是指一次加入三角剖分的一条边,直至产生要求的边数构成的三角剖分;Delaunay三角剖分需要满足这样特性即在三角剖分时获得的任意三角形的外接圆内不包含原始离散点集中的其它点。
绘制等值线时,为了获得较好的等值点内插精度,总希望最临近的数据点构成三角形,即三角网中三角形的边长之和最小。最小权三角剖分满足这一准则,但其求解算法问题至今尚未解决。而贪婪三角剖分是常以局部情况为基础作最优选择而不追求最优解,它只希望得到较为满意解的方法,而不考虑各种可能的整体情况。因此从理论上来讲Delaunay三角剖分被认为是可实现的最优剖分算法。
所谓的规则格网即依据一定的插值原理将原始离散数据点内插为按一定的行距和列距分布在相互垂直的两组平行线的交点上的点数据集。其插入原理是:假设有n个测量数据点(Z1,Z2…Zn),其中Zi=(xi,yi,zi),把这些数据点投影到XOY平面上,得到一个相应的点集,记为{Zi(xi,yi) | i=1,2…n},进而以此点集为依据在XOY平面上按照选定的插值方法建立一定数量行列的规则矩形网格,每一个网格节点即为一个插入点。
对于矩形格网建立采用的插值方法已比较成熟,包括反距离加权、克里金、趋势面法、改进谢德别法、移动平均法等。虽然矩形网格法绘制等值线图具有算法简单的优点,但它要求原始数据点分布较为均匀或均匀。当计算格网中左上部分节点特征值时,由于参考数据较为邻近结果也较为精确,而计算格网右下方向部分节点特征值时,由于原始数据点稀疏,参考点与当前计算点相关度较低导致结果可能与实际相差较大。
通过以上对三角网和规则格网的介绍可知,就对原始数据点的要求而言,三角网格法不仅适合于规则分布的原始数据点,也适合于不规则、甚至畸形分布的原始数据点。而矩形网格法要求采集的原始离散数据点分布较为均匀,被研究的空间变量连续分布,不如三角网格法灵活。就格网化效率而言,在实际网格化过程中,由于一般采集的原始离散数据点分布不均匀,因此使用三角网格法可以省去将原始离散数据点变换成矩形网格节点数据的处理过程,格网化效率较高。就插值准确度而言,由于三角网格的网格节点即是原始数据点,在格网点的插值过程中也是直接依据原始数据计算,因此在原始数据点分布比较均匀的情况下可以生成格网点分布规律且数值精度较高的格网。更重要的是相对于规则矩形格网,三角网构建较为灵活,可以方便地将断层线等线性约束条件嵌入到三角网中,对绘制支持断层等地质构造有重要意义。
根据逐点插入法及优化算法的思想,生成Delaunay三角网的基本步骤如下:
(1)检索数据集查找到数据点的最小与最大横纵坐标分别记录为xmin、xmax、ymin、ymax,从而确定包纳所有的边界点和散点的超级三角形顶点坐标,并将该三角形作为Delaunay三角网的第一条记录。
超级三角形顶点的记录顺序是顺时针的,而且采用逐点插入法生成的Delaunay三角网中所有的三角形将均以顺时针方向记录顶点。这在断层识别中有重要意义。
(2)从离散数据中取一点,遍历已有三角网,判断该离散点是否位于当前三角形外接圆内。如果是,将当前三角形各边放入优化队列,并从三角网内删除该三角形;当遍历完毕,如果散点不在任何一个三角形中,则表明点落在了计算区域之外,应删除此点而不添加到三角网中。反之,按如下方法进行优化:
(a)从优化队列中取出一条边,开始优化此边;
(b)如果此边属于边界边,则此边不用优化;
(c)如果此边所在的正在被优化的三角形的外接形的除此边外的另一个顶点。在此情况下应将当前的优化边从优化队列中删除,同时将该三角形的另外的两边加入到优化队列中删除后送入优化后的队列,便于建立新的三角形;
(d)重复(a)、(b)、(c)步,直到优化队列为空。
(3)重复执行2,直到所有的离散点处理完毕,得到要绘制的Delaunay三角网。
在步骤2执行过程中,最重要的问题是判断点与三角形外接圆的关系。解决方法如图3.1所示,当判断点Q与三角形ABC外接圆关系时,首先计算三角形的外接圆圆心O点坐标,进而比较OQ距离与半径R的大小。
图3.1点与三边形外接圆关系示意图
嵌入断层即将断层线作为约束条件加入已生成的三角网,实现所有断层线线段均作为三角网中某三角形一边,以便于计算等值点和追踪等值线。计算方法如下:
(1)遍历Delaunay三角网,检索并记录作为三角形边的所有断层线线段,对比之下得到断层线段中非三角网某边的断层线线段并为记录;
(2)计算线段中点,并将各中点逐次插入原三角网中,并将该中点按其位置记入断层线;
(3)执行1、2直到断层线中各线段均为Delaunay三角网的一边;
本文断层线节点的赋值思想是将断层线节点特征值与其所处的三角形直接关联,即在包含断层节点的三角形中直接记录断层线节点值。此时,以三角网为基础追踪等值线时可以采用通用方法完成包含断层线的等值线追踪。如在图3.2中红色的线表示断层线,红色的点代表断层线节点,而黑色的点表示其它数据点。当为断层线节点B点赋值时,需在△ABE、△ABJ、△BCE、△BCK、△BKJ中分别记录B点特征值。
解决断层线节点赋值的关键问题在于搜索合理的参考点以计算特征值,由于断层两侧数据的无关性,当为某个断层线节点赋值时需要分别以该断层线两侧数据点集作为参考计算特征值,并将所得的特征值分别赋予其所在三角形。求解断层线与等值线交点关键是为了显示等值线分布趋势而精度要求并非十分精确,因此本文中断层线节点特征值计算时选择参考点为与其直接相关的离散数据点或矩形格网节点,即在三角网中与其共存于同一三角形中的其它数据点。如图3.2中,当计算断层线节点A的特征值时, D、E、G、I、J五点为参考点。
图3.2 含断层的约束三角网示意图
获得参考数据点集后,需要按其与断层位置关系进行分类。为了便于为三角形中断层线节点赋值,本文采用的方法是对当前断层线节点为其顶点的三角形进行分类,进而求解每组三角形中断层线节点的特征值。如图3.2 中,当计算B点特征值时,首先检索出与包含B点的三角形分别为△ABE、△ABJ、△BCE、△BCK、△BKJ。 进而按其与断层线的位置关系将△ABE和△BCE分为一组,△ABJ、△BCK、△BKJ划为另一组。在第一组中特征值计算参考的数据点为J、K点, 而第二组中参考的数据点为E点。获取参考数据点后,即可采用一定插值方法(如克里金、反距离加权)计算特征值并赋予对应的三角形顶点。
本文分类算法的关键思想是当两三角形的公用顶点为非断层线节点时,此三角形即划为一类。该方法的正确性用反证法即可证明,若两三角形的公用顶点为非断层线节点时,且这两个位于断层线两侧,则其中一个三角形的该非断层线节点与其第三个顶点的连线必与断层线相交,这与三角网的基本特性矛盾。应用此方法需要注意一种特殊情况如图3.3所示,同样红色的线表示断层线,红色的点代表断层线节点,而黑色的点表示其它数据点。当计算断层线节点B点时,需△ABE、△ABJ、△BCE、△BCK、△BKJ分组,其中△ABE为独立为一组且该三角形顶点均为断层线节点,此时参考的数据点为与该三角形公用当前断层线节点对边的三角形中的非断层线节点,在图3.3中所示为D点。
图3.3 断层线节点特征值计算特殊情况图示
总结断层线节点计算步骤如下:
(1)检索Delaunay三角网,将所有包含断层线节点的三角形记录成集;
(2)在断层线节点中取出一点;
(3)检索1中的三角形集,记录并存储与当前断层线节点所处的三角形;
(4)将与相前断层线节点相关的三角形按位于断层线方位进行分类并检索获取每一类的参考点集;
(5)分别计算两侧的断层节点特征值,并赋予相应的三角形;
(6)重复执行2-5,直到所有的断层线节点处理完毕。
首先通过整个或局部格网确立一个平面,进而分析实际揭露各点煤层底板标高到该平面的距离在某一方向(纵向或横向)上的数据分布的相关关系,如果呈现出增加到峰值然后降低或减少到谷底然后增加,且近似对称分布,则是褶曲出现的标志。
平面法线式参照方程:
Xcosα+ycosβ+zocsγ-p0=0
平面的确定要求与煤层产状面拟合较好,否则预测准确率较低。因此可以选用局部格网作为拟合平面确立依据,以增加拟合精度。与此同时,需确定一个判断指标以评价其拟合准确度。
煤强度(P)相对而言降低和节理发育强度(J)相对增大,是采掘前方出现断层的两个重要征兆。
据此预测步骤如下:
1)对已知断层的观测数据建立经验公式
Sq=K1H+b1
Sj=K2H+b2
2)连续监测
当P≥A或J≤B时,视为正常征兆;
当P<A或J>B时,则为异常征兆,此时记录Pm1和Jn1点在井下的位置或三维坐标,特别要注意Pm2、Pm3和Jn1、Jn2的观测。
3)动态预测
f0=(P≥A)∩(J≤B)
f1= Pm2∩Jn1∪Pm1∩Jn2
f2= Pm1∩Jn3∪Pm2∩Jn2
f3= Pm2∩Jn3∪Pm3∩Jn2
f4= Pm3∩Jn3
f真值表
模糊函数 |
意义 |
赋值 |
f0 |
不太可能 |
10% |
f1 |
有可能 |
30% |
f2 |
很可能 |
50% |
f3 |
极有可能 |
70% |
f4 |
几乎肯定 |
90% |
半隐伏及完全隐伏断层发育规模的预测依据如下公式:
α——岩层真倾角;
β——断层真倾角;
θ——断层倾向与剖面方位夹角;
ω——岩层走向与断层走向的夹角;
φ——岩层倾向与剖面方位夹角;
ρ——断层走向线与断煤交线在断层面内的夹角;
γ——断层面综合擦痕伏角。
同煤集团同发东周窑矿地质构造极其复杂,在首采面顺槽掘进的过程中多次遇到断层、陷落柱等地质构造,对正常的生产建设造成很大影响,甚至一度因掘全岩而建设进度缓慢。本文研究以此矿资料为基础。
格网化依据的数据来源于两方面,一是勘探过程中的钻孔资料,二是来源于巷道掘进过程中实测的煤层底板标高。首先以此基础数据生成经格网化后生成格网,进而在此基础上,借助上述数据分析方法预测掘进、回采前方的地质构造,收到了较好的效果。
表5.1 工作面两巷揭露标高
点编号 |
横坐标 |
纵坐标 |
标高 |
1 |
52629200 |
4431700 |
878.2 |
2 |
52629300 |
4431700 |
879.3 |
3 |
52629400 |
4431700 |
880.2 |
4 |
52629500 |
4431700 |
877.8 |
5 |
52629600 |
4431700 |
878.2 |
6 |
52629700 |
4431700 |
881.2 |
7 |
52629200 |
4431500 |
870.2 |
8 |
52629300 |
4431500 |
872.2 |
9 |
52629400 |
4431500 |
878.2 |
10 |
52629500 |
4431500 |
873.2 |
11 |
52629600 |
4431500 |
878.2 |
12 |
52629700 |
4431500 |
883.2 |
如上表所示为工作面两巷揭露的部分煤层底板标高,以此数据为基础生成格网,进而绘制煤层底板等高线如图5.1所示。在巷道掘进过程中揭露的煤层底板标高与理论值差异如图5.2所示,进而可以预测该区域断层的存在。
图5.1 以格网为基础绘制的底板等高线图
图5.2 掘进揭露与理论拟合的差异
本文首先提出了采用格网化的方式取代曲面拟合来描述煤层底板面,进而研究了断层等地质构造存在条件下的格网化赋值方式。并以此格网为基础,借助数据分析方法,预测回采工作面前方褶曲、断层等构造的存在及其规模,为地质构造定量化描述和预测提供了依据。为完成上述成果,本文在一些方面做了探索。如以断层线结点赋值的方式,考虑了格网化过程中已知断层对数值分布的影响,为进一步分析提供了更为可靠的基础格网。本文在褶曲的预测中,基于局部格网选择拟合平面,避免了以整个区域拟合平面造成的误差较大的问题。
1.崔洪庆.煤矿生产地质保障技术[M].北京:煤炭工业出版社,2005
2.乔美萍,周荣福等.基于GIS的空间插值法在矿井地质构造预测中的应用研究[J].矿业工程,2010(2)
3.李志华,窦林名,张卫东.矿井地质构造预测与应用[J].煤炭工程.2005(2)
4王国华.煤矿中小型地质构造预测方法[J].水力采煤与管道运输.2006(3)
Research on geological structures in working face of coal mine
Abstract:This paper introduced spatial interpolation method into regional grid transformation, and then based on fruiting grid,combined with the floor level of production practice and using mathematical analysis method, to predict the geological structure of front mining face.
Key words:Spatial interpolation;Grid;geological structure; prediction