尖灭分叉型煤层三维地质建模研究:基于地质统计学和符号矩阵的方法

王艳美, 刘星

(安徽理工大学地球与环境学院, 淮南 232001)

三维地质建模(3D geology modeling),是运用计算机技术,在虚拟三维环境下将空间信息管理、地质解译、空间分析和预测、地学统计、实体内容分析及图形可视化等工具结合起来,用于地质研究的一门技术[1]。为实现煤炭精准开采智慧化,矿井透明地质条件是重要基础[2]。煤矿智能化开采是中国煤炭行业升级转型的必由之路[3],也是煤炭行业未来发展的大趋势,三维地质建模及其可视化技术促进了地质资料信息的管理和共享,在煤炭行业的应用非常广泛且具有极其重要的意义[4]。针对三维地质建模相关研究的各种空间建模技术和方法也不断涌现,根据建模核心过程是否由数学模型控制,可以将其分为显式建模和隐式建模两类[5]。显示建模方法以人机交互为基础,通常基于层面模型[6]来构建,该方法针对区域范围小、构造简单的近平行沉积地层应用效果良好,而针对区域范围大、地层交错复杂、在不同剖面上产状变化大的地质要素进行建模时,需要大量的时间和高要求的专业素养,实际应用程度较低,且所建模型通常与实际地质构造出入较大,模型棱角尖锐,在多z值现象(多逆断层)的地质体中人机都无法对相应地质界线进行连接,无法获得实际的三维地质体,难以进行动态调整且不同的人圈定的地质体差异较大[7]。隐式建模主要核心为空间插值和模拟技术,空间插值的分类方式较多[8-9],常见的空间插值方法有克里金插值、距离反比插值和离散光滑插值等[10]。毕林等[11]提出了一种基于轮廓线和距离场的矿体隐式建模方法,并利用支持向量机和Poisson曲面进行了改进以提高这种方法的自动化程度。邹艳红等[12]提出一套基于隐函数曲面的地质体隐式建模方法,可以快速构建矿体三维模型。Zhang等[13]、唐丙寅等[14]提出一种新的基于不规则三角网与角点网格(triangulate irregular network and corner-point grid,TIN-CPG)的空间数据模型,并基于克里格插值方法实现了三维地层结构与属性模型的一体化展示。郭甲腾等[15]利用全局径向基函数构造的Coons曲面拟合地质特征,并利用布尔操作对模型进行裁剪优化,实现了三维模型的精细构建与系统集成。师素珍等[16]利用地震反演技术进行煤层的分叉合并,证明波阻抗反演对于煤层分叉合并及厚度变化具有较高的精度。空间插值方法针对地层平缓、地质构造简单的地质体建模适用度较高,但对数据量的要求较高,若数据量稀疏容易出现地层穿插现象,且针对复杂地质体的三维地质模型构建过程复杂,难以达到理想状态。因此,不能只单一使用空间插值方法,需要结合其他方法来解决上述问题。

鉴于此,提出一种基于地质统计学和符号矩阵的三维地质建模方法,以适用于大面积区域,不受煤层尖灭或分叉等复杂现象的限制,能够精确定位煤层的位置,并且不限制区域内煤层的层数,使煤层之间的连通性更好。研究成果为煤炭资源精准透明开采设计提供依据。

1 方法

1.1 符号矩阵

符号矩阵是由4个符号学要素组成的显示矩形图式。它是一种不同于数值矩阵的特殊符号表达式[17],其结果可在软件中建模。符号矩阵的计算在许多学科的理论研究与科技计算中都占有重要地位,其一般是二元对立要素。一些简单的符号矩阵计算可手工完成,然而实际的科研工作中往往会产生大量冗繁的符号矩阵或表达式,对此手工计算几乎是不可能的,需要一种不同于数值软件的符号计算工具,对这种符号形式的矩阵进行处理并得到正确结果。主要利用MATLAB软件进行计算,首先对研究区未知位置进行符号距离函数标记,然后再将符号距离函数以{-1, 1}的符号矩阵的形式来标记煤层和非煤层位置,从而有效标识煤层分布情况。

1.2 泛克里金插值方法简介

克里金插值方法是法国G.Matheron教授以南非矿业工程师D.G.Krige(克里格)名字命名的一种对空间分布数据求最优、线性、无偏内插估计量的方法,是地质统计学的重要组成部分,也是地质统计学的核心[18]。由于沉积环境等情况的复杂性,区域化变量在实际应用中很难满足二阶平稳假设,因此选用泛克里金插值方法对待插点煤厚进行预测。在假设非平稳区域化变量Z(x)的变异函数存在,并已知Z(x)的漂移形式为多项式的条件下,讨论各种线性无偏最优估计的问题,称为泛克里金法。区域变量Z(x)在待插值点x0处的真实值Z(x0)是未知的,因此需要用x0影响范围内的测量值的加权平均来估计待插值点的值。

(1)

式(1)中:Z*(x0)为待插值点x0处的煤厚估计值;Z(xi)为已知钻孔的实际测量值;Wi为对应的权重。

求Z*(x0),需要先求权重,其权重不仅取决于待插点与已知点之间的距离h,还与预测区域内的地质环境和变异结构密切相关。变异函数与x无关,只与分隔两点之间的距离h有关,定义为

(2)

式(2)中:h为滞后距;E为方差。

球状模型、指数模型、高斯模型是常用的拟合函数模型[19]。

(1)球状模型的表达式为

(3)

式(3)中:C0为块金效应;C为总基台值;a为变程。

(2)指数模型的表达式为

(4)

(3)高斯模型的表达式为

(5)

泛克里金矩阵表示为

A×W=B

(6)

式(6)中:

其中第n个钻孔到未知待预测第p个钻孔的距离,p为待测区域中未知待预测位置的个数;Xp和Yp为泛克里金估值的位置;a1和a2为待求的偏差参数;Xi、Yi为已知点的坐标,i=1,2,…,n。

求出权重之后,则可利用式(1)求出待插值点处的煤厚。

1.3 反距离加权插值法

反距离加权插值法(inverse distance weighting IDW)是常用的空间插值方法之一,由美国国家气象局提出,是以待插点与已知点之间的距离为权重的插值方法,离待插点越近的已知点被赋予的权重越大,其权重值与距离成反比,其一般表达式为

(7)

(8)

(9)

p越大时较近的已知点将被赋予更高的权重,p越小时权重则更均匀地分配给各已知点。当p=1时称为距离反比法,是一种常用且简便的空间插值方法;当p=2时称为距离平方反比法,选取p=2作为幂指数,来求待插值点的距离函数估值[20]。

2 建模流程

主要利用地质统计学方法和符号矩阵方法对于尖灭、分叉煤层进行三维地质建模,研究的关键是通过已知钻孔数据,利用符号矩阵方法按照一定精度离散煤层和非煤层,通过泛克里金法和反距离加权法对煤层位置进行预测,建立三维地质模型。其流程如图1所示。

图1 模型流程图Fig.1 Model flow chart

2.1 数据准备

通过区域中已知钻孔的实际测量值,求出它的聚煤量和距离函数值,其计算方式如图2所示,聚煤量是指在钻孔处与煤层相交的总厚度,其计算公式为

a、b为钻孔编号;在钻孔a中,该钻孔下有一层煤,T为钻孔a的聚煤量,t1为钻孔a与煤层相交的煤厚,此时T=t1; 在钻孔b中,该钻孔下有四层煤,T为钻孔b的聚煤量,t1为钻孔b与第一层煤层相交的煤厚,t2为钻孔b与第二层煤层相交的煤厚,t3为钻孔b与第三层煤层相交的煤厚,t4为钻孔b与第四层煤层相交的煤厚,此时T=t1+ t2+ t3+t4图2 数据预处理示意图Fig.2 Data preprocessing diagram

(10)

式(10)中:T(a)为在钻孔a处的聚煤量;n(a)为与钻孔a相交的煤层层数;ti(a)为与钻孔a中第i层相交的煤厚。

距离函数与到两域界面的距离有关,两域分别是煤层和非煤层。垂直方向上的栅格点的距离函数为其域到另一个域中最近栅格点的距离,本文定义煤层内距离函数为负值,非煤层内距离函数为正值,则假设两域界面的距离函数为0。每个栅格点是二元对立的,则定义煤层内i(u)=1,非煤层内i(u)=-1。在栅格点u处,确定离它最近的栅格点u′,使得u-u′最小,即为已知点的钻孔在栅格点u处的距离函数F(u),可表示为

(11)

在两域界面处的距离函数为0,非煤层内的距离函数为正,其值在远离界面时增大;煤层内的距离函数为负,其绝对值在远离界面时增大。

2.2 煤厚计算

利用泛克里金法求待插点的煤厚。通过已知钻孔的聚煤量,确定区域化理论变差模型,即可求出变异函数,进而代入式(6),求出权重值,该权重值可用于在每个待插点进行估计,最后利用式(1),求出区域内待插点的煤厚。

2.3 距离函数估值计算

对三维区域进行栅格化,为三维网格中的每个栅格点进行距离函数估值。指定一个垂直搜索窗口来限制使用哪些栅格点对当前栅格点进行估计,并且不限制水平搜索窗口,直接使用煤厚计算求出的权重值,将分配给已知钻孔的泛克里金权重值应用于该钻孔中位于垂直搜索窗口内的每个栅格点。可以直接使用该权重值,也可以利用反距离加权法对泛克里金权重值进行加权,使其在地质环境上更接近被估计栅格点的数据,无论如何权重值之和都为初始分配给已知钻孔的泛克里金权重值。其距离估值的计算公式为

E*(u)=

(12)

式(12)中:E*(u)为待求栅格点的距离估值,F(ui,j)为第i个钻孔垂直搜索窗口中第j个栅格点的距离函数;n为已知钻孔个数;m为垂直搜索窗口中栅格点的个数;d(ui,j)为待求栅格点到第i个钻孔垂直搜索窗口中第j个栅格点的距离。

2.4 根据距离函数估值匹配煤厚

通过计算三维网格中每个栅格点的距离函数可以确定煤层的位置。已知三维网格各列中每个栅格点的深度和距离函数,通过这两个矢量即可确定煤层顶底板。为此,需要考虑距离函数阈值,距离函数估值位于距离函数阈值的相反侧,煤层与非煤层的边界位于阈值的位置。考虑距离函数阈值为0,其示例如图3(a)所示,煤顶在4～4.5 m之间,其距离函数估值为0.35和-0.14;煤底在8～8.5 m,其距离函数估值为-0.27和0.23,则煤厚估计值是4 m。为了使煤厚精确度更高,需对煤顶和煤底进行精细化,如图3(b)所示,初始垂直分辨率为0.5 m,煤厚估计值是4 m,令垂直分辨率为0.05 m,煤顶在4.35～4.4 m,其距离函数估值为0.01和-0.04;煤底在8.25～8.3 m,其距离函数估值为-0.02和0.03,则煤厚估计值是3.9 m,与泛克里金法预测的煤厚更加相近,但距离函数阈值为0 时的煤厚估值与泛克里金预测的煤厚并不总是相近,为此考虑距离函数阈值为0.5 m,其示例如图4所示,煤顶在3.8～3.85 m,其距离函数估值为0.55和0.49;煤底在8.75～8.8 m,其距离函数估值为0.49和0.54,则煤厚估计值是4.95 m。由此可知,距离函数阈值与煤厚估计值成正比,增加距离函数阈值煤厚估计值也增加,减小距离函数阈值煤厚估计值也相应减小。其在多煤层中也同样适用。将距离函数阈值应用于距离函数估值这一列中可以识别不同的煤层和不同的煤层厚度,其示例如图5所示,低距离函数阈值识别

图3 距离函数阈值为0Fig.3 Distance function threshold is 0

图4 距离函数阈值为0.5Fig.4 Distance function threshold is 0.5

图5 距离函数阈值图Fig.5 Distance function threshold diagram

出两个煤层,高距离函数阈值识别出一个煤层,聚煤量也在增加。

2.5 构建三维地质模型

可以绘制三维地质模型的软件有很多,如RMS、GOCAD、EVS、ArcGIS软件等,选择GOCAD作为建模软件[21],首先在该软件上对三维地质模型进行栅格化,生成S-Grid网格模型[22],其为M×N×K个子单元的规则拓扑结构模型,其每个单元体有8个栅格点,其栅格点的坐标可以单独储存,将其以Object的形式导出,对符号矩阵进行计算,其计算公式为

(13)

式(13)中:Classi,j,k为第(i,j,k)个子单元的符号矩阵值,i=1,2,…,M;j=1,2,…,N;k=1,2,…,K,其构建出K个M×N的符号矩阵,将符号矩阵中的值对properties进行重新赋值,重新导入即可看到煤层分布情况。

3 应用实例

3.1 研究区地质背景

潘二煤矿隶属于安徽省淮南市潘集区,矿区范围为116°45′00″E～116°56′15″E,32°45′00″N～32° 50′00″N,其主要含煤地层为石炭系、二叠系,石炭系煤层发育差,无工业价值,二叠系煤层主要为山西组、下、上石河子组[23]。由于研究区范围广,数据量庞大,选取潘集煤矿西北部进行三维地质建模,主要研究6、7、8号煤层群在地下的分布状况,其主要煤层为7-1煤层,厚0～4.95 m,可采指数为0.9,其煤层由于环境变化、后期侵蚀等出现分叉、尖灭情况,煤层结构较为复杂,含暗煤、亮煤,煤层顶底板均以泥岩为主。研究区共有21个钻孔勘探点,其二维分布如图6所示。

图6 研究区钻孔分布图Fig.6 Borehole distribution map of the study area

3.2 模型匹配

已知其研究区内钻孔的聚煤量及距离函数,将该区域内的已知钻孔数据集称之为区域化变量,对其进行模拟,求得变异函数如图7所示,其理论模型为高斯模型, 块金效应C0=0.001,总基台值C=0.954,变程a=2 251.666。利用泛克里金法求出待插点的煤厚,并结合反距离加权法求待插点处的各栅格点的距离函数估值,选取距离函数阈值为0,垂直窗口的大小为5时,其光滑度更高,与泛克里金法求出的煤厚值最为匹配。

图7 变异函数图Fig.7 Variation function diagram

3.3 三维地质建模构建

根据研究区边界范围建立SGrid模型,利用符号矩阵方法,对该研究区内钻孔数据进行预处理,按照一定精度离散煤层和非煤层,进行区域空间建模,通过泛克里金法确定二维网格煤厚,反距离加权法确定三维网格中的距离函数估值,局部窗口对煤层进行精准定位,确定与煤层厚度相匹配的距离函数截距,对各六面体单元属性进行预测及结果可视化,对其进行剖切显示,得到研究区三维地质模型如图8所示,由于该区域不含复杂地质结构,整个区域煤层按一定趋势平缓变化,其煤层分布比较集中,且煤层连通性也比较好,众多薄煤层也可以精准定位,并在适当位置分叉、尖灭,所建三维模型符合实际规律。如图9(a)所示,钻孔14-2与钻孔12-4附近有尖灭现象,在图9(b)中钻孔L4-1和钻孔13-1中煤层附近有分叉现象,其过渡自然,符合地质规律。

图8 研究区三维地质模型栅栏图(z方向拉伸5倍)Fig.8 Three-dimensional geological model fence map of the study area(stretch 5 times in z direction)

图9 三维地质模型中煤层尖灭、分叉剖面与实际钻孔媒段对比Fig.9 Comparison between the coal seam pinch-out and bifurcation sections in the three-dimensional geological model and the actual borehole media section is drawn

创建了煤层剖面示意图,水平250 m,垂直20 m,包括两个相距150 m的钻孔,每个钻孔与4 m厚的煤层相交。图10说明改变垂直窗口个数的效果,在每种情况下,两个钻孔相交接缝之间的垂直偏移和煤层位置都是相同的。初始煤层呈断裂状态,随着垂直窗口大小的增加,接缝的连接变得平滑。结果表明:垂直窗口较小时,煤层会出现不连续现象,实际应用时,选择垂直窗口为7时满足建模要求,图11展示了改变煤层位置的效果,在每种情况下,垂直窗口大小保持不变,煤层连通性随着两钻孔相交煤层的位置的距离增加而变化。初始煤层呈平滑状态,接缝似乎一直在拉伸,在两钻孔相交煤层的位置为10 m时,煤层开始断裂,直到最后以约6 m的间距“断裂”。两钻孔相交煤层距离在10 m之内,煤层不会出现断裂现象,在超过10 m时出现断裂,与地质规律一致。

图10 垂直窗口变化后对煤层连通性影响Fig.10 Influence of vertical window change on coal seam connectivity

图11 煤层位置变化后对煤层连通性影响Fig.11 Influence of coal seam connectivity after the change of coal seam position

3.4 模型检验

在研究区内随机抽取3个参与建模钻孔,将原始钻孔数据与模型预测数据对比,其结果如表1所示。由表1可知,二者具有良好的一致性,精度约为0.3 m,满足开采要求[24]。

表1 原始钻孔数据与模型预测数据精度对比Table 1 Comparison of accuracy between original borehole data and model prediction data

4 结论

(1)以钻孔数据为基础,提出一种基于地质统计学和符号矩阵的三维地质建模方法,可以快速建立分叉、尖灭连接项,精度可靠,误差小于0.3 m,满足开采要求。

(2)本文方法使用符号矩阵方法来构建煤层几何模型。利用符号矩阵方法计算出的距离函数被局部截断匹配正确的煤层厚度,接缝自然向外扩张并分叉,局部窗口对对距离函数模型进行精细化,以实现高精度的单个煤层的标高,且跟传统方法相比,本文方法减少人工交互,自动化程度提高,避免层面建模需要人工经验,先解译再连接的难题。

(3)实践表明,通过在潘集煤矿中使用该方法,成功地预测出7号煤在地层下的分布情况,符合地质规律,该方法中的算法运行速度和精确定位煤层和非煤层边界的能力使其成为一个强大的建模工具,它可以模拟地下的煤层分布,使地质体三维可视化,为煤炭资源精准透明开采设计提供依据。