基于VTK的地下矿开采沉陷预计研究与实现

朱忠华，王李管，毕林，钟德云

(中南大学资源与安全工程学院，中南大学数字矿山研究中心，湖南长沙 410083)

0 引言

地下矿开采破坏了采区周围岩体的完整性和应力平衡，使上覆岩层断裂、移动．随着开采的进行，这种移动逐渐向上延伸，不可避免地使地表产生移动、变形、裂缝、沉降、塌陷，造成不同程度的破坏，对生态环境构成严重影响，甚至威胁人民群众的生产、生活．波兰在该国和具托姆车站下开采，采厚20 m，车站普遍下沉达3 m［1］;德国柏留克城曾因开采导致地面突然塌陷，毁坏了31座房屋．我国是开采沉陷的重灾区，地下矿开采导致矿区地表下沉，造成严重的影响，如安徽淮北市，因采煤沉陷，全市耕地面积逐年减少，30多年来已造成地表沉陷超1．13万hm2［2］．此外，我国“三下”压煤量巨大，达133亿多吨［1］，其中兖州矿区“三下”压煤占其可采储量51．1%［3］，其开采环境复杂，这部分资源的合理利用也有赖于开采沉陷的合理预计．因此，方便准确地进行开采沉陷预计具有重要的意义．

科研人员和矿山工程技术人员对开采沉陷预计工作投入了越来越多的研究．胡振琪［2］对采煤沉陷的土地资源管理与复垦进行了研究，但未给出有效的沉陷预计方法;谢和平等［4］将有限差分软件Flac3d应用到开采沉陷预计．沉陷预计的随机介质理论，首先由波兰学者李特威尼申(J．Litwiniszyn)于20世纪50年代引入岩层移动研究，后由我国学者刘宝琛、廖国华等发展为概率积分法，经过多年的研究，该方法目前已成为我国较成熟的、应用最为广泛的预计方法之一［5］．国内学者对此做了一些研究［6－8］，然而离矿山企业的实际需求还有一定差距．

VTK(visualization toolkit)是美国Kitware公司的大型可视化工具包．VTK提供了开源面向对象的可视化类库，其具有较好的底层数据结构，方便建立可视化管线，在众多领域尤其是医疗可视化得到了较好应用［9－12］，但在矿业和岩土工程中的应用尚处于起步阶段［13－14］．本文以VTK可视化工具包为平台，采用C++编程语言扩展VTK，将开采沉陷预计算法写入VTK过滤器，实现地下矿开采沉陷预计;通过扩展VTK过滤器来拓展其专业应用．

1 VTK工具包简介

VTK内核由C++构建，它将可视化过程中常用的算法封装成多个类，可方便地对数据集进行各种变换和操作，实现三维数据场的可视化．VTK主要由C++类库和解译包装层两个基本子系统组成，可在C++、Tcl/Tk、Java、Python语言环境下使用，其内部采用流水管线机制，如图1所示．

VTK可视化管线由数据对象(data objects)和过程对象(process objects)组成．数据对象代表数据信息;过程对象在输入数据上进行操作，产生输出数据．过程对象按对象开始，继续和终止可视化数据流可以分为源对象(source objects)、过滤器对象(filter objects)和映射器对象(mapper objects)．过滤器对象用来对数据进行各种转换和处理，要有一个或者多个输入数据，产生一个或者多个输出数据．映射器对象通常用来将数据转换为基本图元，也可以将数据存入文档，或者用来建立与其它软件系统或设备的接口．

用户根据实际情况，可以建立自己需要的可视化管线，也可以重写或者扩展过滤器对象(filter objects)．

图1 VTK可视化流水管线Fig．1 The visualization pipeline of VTK

2 开采沉陷预计理论基础及VTK实现流程

2．1 开采沉陷预计理论基础

地下矿开采引起周围岩体和地表移动，开采沉陷预计就是定量地研究受开采影响的岩层和地表移动在时间和空间上的分布规律．开采沉陷预计的常用方法有典型曲线法、剖面函数法和概率积分法［5］．本文采用基于随机介质理论的概率积分法沉陷预计模型，下沉量为:

式中:W0为最大下沉量，W0=mq cosα，m为煤层厚，q为下沉系数，cosα为煤层倾角余弦;D1s和D3为开采工作面范围尺寸;r为主要影响半径，r=H/tanβ，H为地面待计算任意点A(x，y)与煤层上微元B(s，t)的高程差，tanβ为主要影响角正切．

此外，开采沉陷预计使用的开采影响传播角θ与煤层倾角α相关，θ=90°－κα，κ取0．6～0．8［5］;拐点偏距s0表示由于煤壁上方煤岩层的悬臂作用引起的拐点的偏移距离，在式(1)中未能体现，需在计算中考虑．

由文献［5］可知，点A(x，y)在沿φ方向的倾斜i(x，y，φ)和曲率k(x，y，φ)分别为:

而水平移动U(x，y，φ)和水平变形ε(x，y，φ)分别为:

式(2)～式(5)中:相关参数右上角的上标“°”表示有限开采的移动变形值;参数b为水平移动系数．

2．2 基于VTK的开采沉陷预计流程

开采沉陷对地表的影响通过地表移动和变形量来表征和衡量．通过开采沉陷预计需要得到的也是地表的移动和变形量，主要是:水平移动、垂直移动(下沉)、倾斜变形、曲率变形和水平变形，因此开采沉陷预计需要实现这些移动变形量的计算，以及对预计结果的提取、处理和对沉陷地表动态更新．

开采沉陷预计将开采工作面和地表作为输入，通过预计模型计算得到沉陷结果，主要流程如图2所示．

图2 地表开采沉陷预计流程Fig．2 Flow chart of surface subsidence prediction

地表建模在VTK现有的类基础上扩展完成，得到的地表DTM(digital terrain model)作为开采沉陷预计模型的输入．开采工作面参数用两个结构体struc_wpborder和struc_wpPara传入，其中struc_wpborder存放工作面边界极值、走向和中心点位置，struc_wpPara存放煤层开采厚度、下沉系数、倾角、主要影响角正切、水平移动系数及平均深度．沉陷预计根据开采沉陷预计模型，在VTK中开发对应的功能，得到沉陷预计的结果;后处理则是将沉陷预计的结果进行提取和处理，由VTK现有的类完成．

3 基于VTK开采沉陷预计实现

3．1 基于VTK的地表建模

地表开采沉陷预计的输入为地表模型DTM，DTM模型由矿山测量得到的点和线数据经过三角化算法得到．Delaunay三角化是实际应用中使用最多的三角剖分算法［15－17］．VTK现有的三角化由vtkDelaunay2D类提供，该类可以对输入的点线数据进行Delaunay三角化(delaunay triangulation)及限制Delaunay三角化(consrained delaunay triangulation)．但是Delaunay三角化存在边界消失的不足，而限制Delaunay三角化虽解决了边界消失的问题，但仍存在单元形状不佳的情况［17］，为克服这些不足，更好地反应地表沉陷，需要扩展VTK三角化算法．本文在vtkDelaunay2D类的基础上，将相容Delaunay三角化(conforming delaunay triangulation)引入到VTK中，命名为vtkConformingDelaunay2D类．在使用vtkConformingDelaunay2D建模得到的二维模型基础上，给每个顶点赋以高程，即可完成三维地表建模．

3．1．1 相容Delaunay三角化

相容的Delaunay三角化在消失的线段上插入新的节点来恢复消失的线段，新增点称为史坦纳点(Steiner point)［17］，相容的Delaunay三角化图例见图3．

图3 相容的Delaunay三角化图例Fig．3 Conforming Delaunay triangulation

3．1．2 vtk ConformingDelaunay2D 类进行地表建模

vtkConformingDelaunay2D类依据相容Delaunay三角化规则，由vtkDelaunay2D类扩展而来，限于篇幅，在此不对其详述实现过程．

矿山原始的点线数据，设置到数据集vtkPolyData对象，作为相容Delaunay三角化的输入，经过vtk-ConformingDelaunay2D对象三角化后可以得到高质量的地表模型．使用vtkConformingDelaunay2D进行地表建模的部分代码如下:

vtkPolyData*profile=vtkPolyData::New();

profile－ ＞SetPoints(points);//设置点

profile－＞SetPolys(poly);//设置线的拓扑

vtkConformingDelaunay2D*del=vtkConformingDelaunay2D::New();

del－ ＞SetInput(profile);

del－ ＞Excute();

vtkPolyData*profileDTM=NULL;

profileDTM=del－ ＞GetOutput();

3．1．3 可视化管线

开采沉陷预计VTK可视化数据集vtkDataSet选用多边形数据vtkPolyData，它由点、线、多段线、三角形、四边形、多边形、三角条带等基本单元构成，大多数过滤器都可以处理这种类型的数据．地表建模可视化实现的部分代码如下:

vtkPolyDataMapper*mapper=vtkPolyDataMapper::New();

mapper－ ＞ SetInput(profileDTM);

vtkActor*acror=vtkActor::New();

acror－＞SetMapper(mapper);

Render－ ＞AddActor(actor);

基于VTK地表建模得到的地表模型见图4．

图4 初始地表模型Fig．4 The architecture of initial DTM

3．2 VTK开采沉陷预计过滤器设计

开采沉陷预计属于矿业及岩土工程的专业方向，是本文的核心算法，VTK本身不带有该类算法，因此需要重新开发．本文将开采沉陷预计算法写入VTK过滤器(Filter)，命名为vtkSubsidenceFilter，以实现地表的沉陷预计．

vtkSubsidenceFilter输入地表采用3．1．2节建模得到的地表profileDTM，沉陷预计功能用void CacuSubsidence()方法实现．由本文2．1节可知，沉陷预计模型实质上为二重积分，参照文献［18－19］，本文沉陷预计采用变步长Simpson二重积分计算，即将该积分对2．1节的开采沉陷预计模型的求解用C++编程语言写入CacuSubsidence()函数，作为vtkSubsidenceFilter类的成员函数，则可以计算地表模型的沉陷值．用vtkSubsidenceFilter类进行沉陷预计的部分代码如下:

vtkSubsidenceFilter*filter=vtkSubsidenceFilter::New();

filter－ ＞ SetInput(profileDTM);//设置点

filter－＞SetParameter1(wpborder);//传入工作面参数，wpborder为struc_wpborder结构体变量

filter－＞SetParameter2(wpPara);//传入工作面参数，wpPara为struc_wpPara类型变量

filter－＞CacuSubsidence();

vtkPolyData*profileDTM=NULL;

profileDTM=del－ ＞GetOutput();

vtkSubsidenceFilter类实现了沉陷值的计算，同时将沉陷值返回到地表，得到沉陷后的地表，从而实现了地表的动态更新，故该类的输出是沉陷后的三维地表模型．

3．3 属性数据提取与后处理

3．3．1 属性数据提取与vtkDataArray类

属性数据用于描述数据集的点或者单元属性信息，在VTK中可以用vtkDataArray类存储属性数据，属性数据可以是标量、矢量、张量及纹理坐标等形式，每个属性数据分别与vtkDataSet中的点和单元相关联．该类定义了很多应用程序接口．

vtkDataArray类是所有数据数组对象的抽象超类，用来存储连续的同一类型的数据，如char、int、float等．在定义数据数组的过程中必须动态为其分配内存空间．vtkDataArray中的元素由组元(tuple)组成，每一个组元是一个结构体，该结构体由N个数据类型相同的分量组成．数据数组可以存储的数据包括几何数据和属性数据(矢量、向量、张量)等．

开采沉陷预计的属性数据有沉降值、水平移动、水平变形、曲率变形和倾斜变形等，均为矢量，数据类型均为float，其中沉降值的方向垂直向下，可视为标量．水平移动、水平变形、曲率变形及倾斜变形为矢量，可以将其分解为X和Y方向的分量分别计算和提取，方法和沉降值的计算与提取类似．

3．3．2 开采沉陷预计后处理

开采沉陷预计后处理可用VTK现有的vtkMarchingContourFiltter、vtkCutter、vtkThresholdd等过滤器．vtkMarchingContourFiltter的输入可以是继承于vtkDataSet的任何数据类型，用于产生等值线或等值面;vtkCutter过滤器用于对数据集进行切片;而vtkThresholdd过滤器的功能是提取那些满足阈值范围的cells．

结果后处理的可视化管线建立与3．1．3节类似．通过后处理，可从计算结果数据中提取等值线、生成结果云图以及完成沉陷地表的属性配色，使结果表达更丰富和形象．

4 应用实例

某矿区地质采矿条件:采煤方法为综采放顶煤，顶板管理为全部垮落法;开采煤层为6煤，该煤层倾角为α=12°，开采深度为202～230 m，开采厚度8 m．该矿区概率积分法参数为:q=0．76，tanβ=2．0，b=0．36，θ=82．8°，s0=13m．

初始地表和对应开采工作面如图5所示，开采工作面从上到下分别为602、604、606、608、610、612工作面及6202、6204和6206工作面，开采沉陷预计结果如图6，开采沉陷影响面积见表1．

表1 开采沉陷影响范围Tab．1 Scope of mining subsidence

图5 输入地表和开采工作面Fig．5 The input DTM and working face model

图6 开采沉陷预计结果Fig．6 The calculation result of mining subsidence

从图6可知，沉陷预计最大值为5．22 m，与最大实测值5．07 m较接近;沉陷最大值发生在606工作面，预计倾斜最大值发生在沉陷等值线最密处，与实际情况相符;由表1可知，沉陷区面积与实际较吻合，计算具有较好的可靠性．

5 结语

基于VTK设计了vtkConformingDelaunay2D类，实现相容Delaunay三角化，得到高质量的地表模型;根据沉陷预计模型，设计开采沉陷预计过滤器vtkSubsidenceFilter，实现沉陷计算;基于VTK实现沉陷预计结果的实时获取和沉陷地表的动态更新．该方法突破了以往VTK主要用于三维可视化的局限，通过扩展VTK过滤器，拓宽了其专业应用．研究成果可用于地下煤矿与金属矿开采沉陷预计、开采沉陷的土地资源管理与复垦、以及自然崩落法放矿的矿岩接触面模拟．应用实例表明，这种应用取得了较好效果．

［1］侯长祥，冯涛，熊仁钦，等．矿床“三下一上”开采［M］．北京:煤炭工业出版社，2001．

［2］胡振琪．采煤沉陷地的土地资源管理与复垦［M］．北京:煤炭工业出版社，1996．

［3］赵经彻，高延法，张怀新．兖州矿区开采沉陷控制的研究［J］．煤炭学报，1997，22(3):248－252．

［4］谢和平，周宏伟，王金安，等．FLAC3D在煤矿开采沉陷预测中的应用及对比分析［J］．岩石力学与工程学报，1999，18(4):397－401．

［5］何国清，杨伦，凌赓娣，等．矿山开采沉陷学［M］．徐州:中国矿业大学出版社，1991．

［6］陈雨，张振文，张彦敏．基于概率积分法的地面沉陷灾害预测［J］．辽宁工程技术大学学报，2007，26(11):143－145．

［7］顾叶，宋振柏，张胜伟．基于概率积分法的开采沉陷预计研究［J］．山东理工大学学报:自然科学版，2011，25(1):33－36．

［8］蔡美峰，李春雷，谢谟文，等．北洺河铁矿开采沉陷预计及地表变形监测与分析［J］．北京科技大学学报，2008，30(2):109－114．

［9］陈晓军，阎世举，吴轶群，等．基于VTK的口腔种植手术导航系统的图像配准［J］．计算机工程，2006，32(23):4－6．

［10］王延华，洪飞，吴恩华．基于VTK的医学图像处理子系统设计和实现［J］．计算机工程与应用，2003，39(8):205－207．

［11］Schroeder W，Martin K，Lorensen B．The visualization toolkit:an object－ oriented to 3D graphics［M］．3rd ．［s．l．］:Kitware Inc，2002．

［12］刘伟宁，陈耀武，张磊．海底声呐实时可视化显示系统设计［J］．计算机工程，2010，36(23):249－251．

［13］毕林，王李管，陈建宏，等．基于VTK的矿体三维可视化研究与实现［J］．计算机工程与应用，2008，44(10):78－81．

［14］欧耿鑫，陈喜，佘超，等．基于VTK的地下水数值模拟三维可视化开发应用［J］．水文，2009，29(1):17－20．

［15］谭云兰，李光耀，夏洁武，等．Delaunay三角网高效构建及地形仿真应用［J］．计算机工程，2012，38(22):287－290．

［16］陈学工，黄晶晶．Delaunay三角剖分中的约束边嵌入算法［J］．计算机工程，2007，33(16):56－58．

［17］李智欣．Triangle与TetGen网格产生器于有限元素程式之使用［M］．台南:国立成功大学，2004．

［18］李庆杨，王能超，易大义．数值分析［M］．武汉:华中科技大学出版社，2001．

［19］蒋和理．二重积分优化梯形与辛卜生算法［J］．合肥工业大学学报，1987，9(6):34－36．