基于VTK的地下礦開采沉陷預(yù)計(jì)研究與實(shí)現(xiàn)

朱忠華，王李管，畢林，鐘德云

(中南大學(xué)資源與安全工程學(xué)院，中南大學(xué)數(shù)字礦山研究中心，湖南長(zhǎng)沙 410083)

0 引言

地下礦開采破壞了采區(qū)周圍巖體的完整性和應(yīng)力平衡，使上覆巖層斷裂、移動(dòng)．隨著開采的進(jìn)行，這種移動(dòng)逐漸向上延伸，不可避免地使地表產(chǎn)生移動(dòng)、變形、裂縫、沉降、塌陷，造成不同程度的破壞，對(duì)生態(tài)環(huán)境構(gòu)成嚴(yán)重影響，甚至威脅人民群眾的生產(chǎn)、生活．波蘭在該國(guó)和具托姆車站下開采，采厚20 m，車站普遍下沉達(dá)3 m［1］;德國(guó)柏留克城曾因開采導(dǎo)致地面突然塌陷，毀壞了31座房屋．我國(guó)是開采沉陷的重災(zāi)區(qū)，地下礦開采導(dǎo)致礦區(qū)地表下沉，造成嚴(yán)重的影響，如安徽淮北市，因采煤沉陷，全市耕地面積逐年減少，30多年來(lái)已造成地表沉陷超1．13萬(wàn)hm2［2］．此外，我國(guó)“三下”壓煤量巨大，達(dá)133億多噸［1］，其中兗州礦區(qū)“三下”壓煤占其可采儲(chǔ)量51．1%［3］，其開采環(huán)境復(fù)雜，這部分資源的合理利用也有賴于開采沉陷的合理預(yù)計(jì)．因此，方便準(zhǔn)確地進(jìn)行開采沉陷預(yù)計(jì)具有重要的意義．

科研人員和礦山工程技術(shù)人員對(duì)開采沉陷預(yù)計(jì)工作投入了越來(lái)越多的研究．胡振琪［2］對(duì)采煤沉陷的土地資源管理與復(fù)墾進(jìn)行了研究，但未給出有效的沉陷預(yù)計(jì)方法;謝和平等［4］將有限差分軟件Flac3d應(yīng)用到開采沉陷預(yù)計(jì)．沉陷預(yù)計(jì)的隨機(jī)介質(zhì)理論，首先由波蘭學(xué)者李特威尼申(J．Litwiniszyn)于20世紀(jì)50年代引入巖層移動(dòng)研究，后由我國(guó)學(xué)者劉寶琛、廖國(guó)華等發(fā)展為概率積分法，經(jīng)過(guò)多年的研究，該方法目前已成為我國(guó)較成熟的、應(yīng)用最為廣泛的預(yù)計(jì)方法之一［5］．國(guó)內(nèi)學(xué)者對(duì)此做了一些研究［6－8］，然而離礦山企業(yè)的實(shí)際需求還有一定差距．

VTK(visualization toolkit)是美國(guó)Kitware公司的大型可視化工具包．VTK提供了開源面向?qū)ο蟮目梢暬悗?kù)，其具有較好的底層數(shù)據(jù)結(jié)構(gòu)，方便建立可視化管線，在眾多領(lǐng)域尤其是醫(yī)療可視化得到了較好應(yīng)用［9－12］，但在礦業(yè)和巖土工程中的應(yīng)用尚處于起步階段［13－14］．本文以VTK可視化工具包為平臺(tái)，采用C++編程語(yǔ)言擴(kuò)展VTK，將開采沉陷預(yù)計(jì)算法寫入VTK過(guò)濾器，實(shí)現(xiàn)地下礦開采沉陷預(yù)計(jì);通過(guò)擴(kuò)展VTK過(guò)濾器來(lái)拓展其專業(yè)應(yīng)用．

1 VTK工具包簡(jiǎn)介

VTK內(nèi)核由C++構(gòu)建，它將可視化過(guò)程中常用的算法封裝成多個(gè)類，可方便地對(duì)數(shù)據(jù)集進(jìn)行各種變換和操作，實(shí)現(xiàn)三維數(shù)據(jù)場(chǎng)的可視化．VTK主要由C++類庫(kù)和解譯包裝層兩個(gè)基本子系統(tǒng)組成，可在C++、Tcl/Tk、Java、Python語(yǔ)言環(huán)境下使用，其內(nèi)部采用流水管線機(jī)制，如圖1所示．

VTK可視化管線由數(shù)據(jù)對(duì)象(data objects)和過(guò)程對(duì)象(process objects)組成．?dāng)?shù)據(jù)對(duì)象代表數(shù)據(jù)信息;過(guò)程對(duì)象在輸入數(shù)據(jù)上進(jìn)行操作，產(chǎn)生輸出數(shù)據(jù)．過(guò)程對(duì)象按對(duì)象開始，繼續(xù)和終止可視化數(shù)據(jù)流可以分為源對(duì)象(source objects)、過(guò)濾器對(duì)象(filter objects)和映射器對(duì)象(mapper objects)．過(guò)濾器對(duì)象用來(lái)對(duì)數(shù)據(jù)進(jìn)行各種轉(zhuǎn)換和處理，要有一個(gè)或者多個(gè)輸入數(shù)據(jù)，產(chǎn)生一個(gè)或者多個(gè)輸出數(shù)據(jù)．映射器對(duì)象通常用來(lái)將數(shù)據(jù)轉(zhuǎn)換為基本圖元，也可以將數(shù)據(jù)存入文檔，或者用來(lái)建立與其它軟件系統(tǒng)或設(shè)備的接口．

用戶根據(jù)實(shí)際情況，可以建立自己需要的可視化管線，也可以重寫或者擴(kuò)展過(guò)濾器對(duì)象(filter objects)．

圖1 VTK可視化流水管線Fig．1 The visualization pipeline of VTK

2 開采沉陷預(yù)計(jì)理論基礎(chǔ)及VTK實(shí)現(xiàn)流程

2．1 開采沉陷預(yù)計(jì)理論基礎(chǔ)

地下礦開采引起周圍巖體和地表移動(dòng)，開采沉陷預(yù)計(jì)就是定量地研究受開采影響的巖層和地表移動(dòng)在時(shí)間和空間上的分布規(guī)律．開采沉陷預(yù)計(jì)的常用方法有典型曲線法、剖面函數(shù)法和概率積分法［5］．本文采用基于隨機(jī)介質(zhì)理論的概率積分法沉陷預(yù)計(jì)模型，下沉量為:

式中:W0為最大下沉量，W0=mq cosα，m為煤層厚，q為下沉系數(shù)，cosα為煤層傾角余弦;D1s和D3為開采工作面范圍尺寸;r為主要影響半徑，r=H/tanβ，H為地面待計(jì)算任意點(diǎn)A(x，y)與煤層上微元B(s，t)的高程差，tanβ為主要影響角正切．

此外，開采沉陷預(yù)計(jì)使用的開采影響傳播角θ與煤層傾角α相關(guān)，θ=90°－κα，κ取0．6～0．8［5］;拐點(diǎn)偏距s0表示由于煤壁上方煤巖層的懸臂作用引起的拐點(diǎn)的偏移距離，在式(1)中未能體現(xiàn)，需在計(jì)算中考慮．

由文獻(xiàn)［5］可知，點(diǎn)A(x，y)在沿φ方向的傾斜i(x，y，φ)和曲率k(x，y，φ)分別為:

而水平移動(dòng)U(x，y，φ)和水平變形ε(x，y，φ)分別為:

式(2)～式(5)中:相關(guān)參數(shù)右上角的上標(biāo)“°”表示有限開采的移動(dòng)變形值;參數(shù)b為水平移動(dòng)系數(shù)．

2．2 基于VTK的開采沉陷預(yù)計(jì)流程

開采沉陷對(duì)地表的影響通過(guò)地表移動(dòng)和變形量來(lái)表征和衡量．通過(guò)開采沉陷預(yù)計(jì)需要得到的也是地表的移動(dòng)和變形量，主要是:水平移動(dòng)、垂直移動(dòng)(下沉)、傾斜變形、曲率變形和水平變形，因此開采沉陷預(yù)計(jì)需要實(shí)現(xiàn)這些移動(dòng)變形量的計(jì)算，以及對(duì)預(yù)計(jì)結(jié)果的提取、處理和對(duì)沉陷地表動(dòng)態(tài)更新．

開采沉陷預(yù)計(jì)將開采工作面和地表作為輸入，通過(guò)預(yù)計(jì)模型計(jì)算得到沉陷結(jié)果，主要流程如圖2所示．

圖2 地表開采沉陷預(yù)計(jì)流程Fig．2 Flow chart of surface subsidence prediction

地表建模在VTK現(xiàn)有的類基礎(chǔ)上擴(kuò)展完成，得到的地表DTM(digital terrain model)作為開采沉陷預(yù)計(jì)模型的輸入．開采工作面參數(shù)用兩個(gè)結(jié)構(gòu)體struc_wpborder和struc_wpPara傳入，其中struc_wpborder存放工作面邊界極值、走向和中心點(diǎn)位置，struc_wpPara存放煤層開采厚度、下沉系數(shù)、傾角、主要影響角正切、水平移動(dòng)系數(shù)及平均深度．沉陷預(yù)計(jì)根據(jù)開采沉陷預(yù)計(jì)模型，在VTK中開發(fā)對(duì)應(yīng)的功能，得到沉陷預(yù)計(jì)的結(jié)果;后處理則是將沉陷預(yù)計(jì)的結(jié)果進(jìn)行提取和處理，由VTK現(xiàn)有的類完成．

3 基于VTK開采沉陷預(yù)計(jì)實(shí)現(xiàn)

3．1 基于VTK的地表建模

地表開采沉陷預(yù)計(jì)的輸入為地表模型DTM，DTM模型由礦山測(cè)量得到的點(diǎn)和線數(shù)據(jù)經(jīng)過(guò)三角化算法得到．Delaunay三角化是實(shí)際應(yīng)用中使用最多的三角剖分算法［15－17］．VTK現(xiàn)有的三角化由vtkDelaunay2D類提供，該類可以對(duì)輸入的點(diǎn)線數(shù)據(jù)進(jìn)行Delaunay三角化(delaunay triangulation)及限制Delaunay三角化(consrained delaunay triangulation)．但是Delaunay三角化存在邊界消失的不足，而限制Delaunay三角化雖解決了邊界消失的問(wèn)題，但仍存在單元形狀不佳的情況［17］，為克服這些不足，更好地反應(yīng)地表沉陷，需要擴(kuò)展VTK三角化算法．本文在vtkDelaunay2D類的基礎(chǔ)上，將相容Delaunay三角化(conforming delaunay triangulation)引入到VTK中，命名為vtkConformingDelaunay2D類．在使用vtkConformingDelaunay2D建模得到的二維模型基礎(chǔ)上，給每個(gè)頂點(diǎn)賦以高程，即可完成三維地表建模．

3．1．1 相容Delaunay三角化

相容的Delaunay三角化在消失的線段上插入新的節(jié)點(diǎn)來(lái)恢復(fù)消失的線段，新增點(diǎn)稱為史坦納點(diǎn)(Steiner point)［17］，相容的Delaunay三角化圖例見圖3．

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

3．1．2 vtk ConformingDelaunay2D 類進(jìn)行地表建模

vtkConformingDelaunay2D類依據(jù)相容Delaunay三角化規(guī)則，由vtkDelaunay2D類擴(kuò)展而來(lái)，限于篇幅，在此不對(duì)其詳述實(shí)現(xiàn)過(guò)程．

礦山原始的點(diǎn)線數(shù)據(jù)，設(shè)置到數(shù)據(jù)集vtkPolyData對(duì)象，作為相容Delaunay三角化的輸入，經(jīng)過(guò)vtk-ConformingDelaunay2D對(duì)象三角化后可以得到高質(zhì)量的地表模型．使用vtkConformingDelaunay2D進(jìn)行地表建模的部分代碼如下:

vtkPolyData*profile=vtkPolyData::New();

profile－ ＞SetPoints(points);//設(shè)置點(diǎn)

profile－＞SetPolys(poly);//設(shè)置線的拓?fù)?/p>

vtkConformingDelaunay2D*del=vtkConformingDelaunay2D::New();

del－ ＞SetInput(profile);

del－ ＞Excute();

vtkPolyData*profileDTM=NULL;

profileDTM=del－ ＞GetOutput();

3．1．3 可視化管線

開采沉陷預(yù)計(jì)VTK可視化數(shù)據(jù)集vtkDataSet選用多邊形數(shù)據(jù)vtkPolyData，它由點(diǎn)、線、多段線、三角形、四邊形、多邊形、三角條帶等基本單元構(gòu)成，大多數(shù)過(guò)濾器都可以處理這種類型的數(shù)據(jù)．地表建模可視化實(shí)現(xiàn)的部分代碼如下:

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開采沉陷預(yù)計(jì)過(guò)濾器設(shè)計(jì)

開采沉陷預(yù)計(jì)屬于礦業(yè)及巖土工程的專業(yè)方向，是本文的核心算法，VTK本身不帶有該類算法，因此需要重新開發(fā)．本文將開采沉陷預(yù)計(jì)算法寫入VTK過(guò)濾器(Filter)，命名為vtkSubsidenceFilter，以實(shí)現(xiàn)地表的沉陷預(yù)計(jì)．

vtkSubsidenceFilter輸入地表采用3．1．2節(jié)建模得到的地表profileDTM，沉陷預(yù)計(jì)功能用void CacuSubsidence()方法實(shí)現(xiàn)．由本文2．1節(jié)可知，沉陷預(yù)計(jì)模型實(shí)質(zhì)上為二重積分，參照文獻(xiàn)［18－19］，本文沉陷預(yù)計(jì)采用變步長(zhǎng)Simpson二重積分計(jì)算，即將該積分對(duì)2．1節(jié)的開采沉陷預(yù)計(jì)模型的求解用C++編程語(yǔ)言寫入CacuSubsidence()函數(shù)，作為vtkSubsidenceFilter類的成員函數(shù)，則可以計(jì)算地表模型的沉陷值．用vtkSubsidenceFilter類進(jìn)行沉陷預(yù)計(jì)的部分代碼如下:

vtkSubsidenceFilter*filter=vtkSubsidenceFilter::New();

filter－ ＞ SetInput(profileDTM);//設(shè)置點(diǎn)

filter－＞SetParameter1(wpborder);//傳入工作面參數(shù)，wpborder為struc_wpborder結(jié)構(gòu)體變量

filter－＞SetParameter2(wpPara);//傳入工作面參數(shù)，wpPara為struc_wpPara類型變量

filter－＞CacuSubsidence();

vtkPolyData*profileDTM=NULL;

profileDTM=del－ ＞GetOutput();

vtkSubsidenceFilter類實(shí)現(xiàn)了沉陷值的計(jì)算，同時(shí)將沉陷值返回到地表，得到沉陷后的地表，從而實(shí)現(xiàn)了地表的動(dòng)態(tài)更新，故該類的輸出是沉陷后的三維地表模型．

3．3 屬性數(shù)據(jù)提取與后處理

3．3．1 屬性數(shù)據(jù)提取與vtkDataArray類

屬性數(shù)據(jù)用于描述數(shù)據(jù)集的點(diǎn)或者單元屬性信息，在VTK中可以用vtkDataArray類存儲(chǔ)屬性數(shù)據(jù)，屬性數(shù)據(jù)可以是標(biāo)量、矢量、張量及紋理坐標(biāo)等形式，每個(gè)屬性數(shù)據(jù)分別與vtkDataSet中的點(diǎn)和單元相關(guān)聯(lián)．該類定義了很多應(yīng)用程序接口．

vtkDataArray類是所有數(shù)據(jù)數(shù)組對(duì)象的抽象超類，用來(lái)存儲(chǔ)連續(xù)的同一類型的數(shù)據(jù)，如char、int、float等．在定義數(shù)據(jù)數(shù)組的過(guò)程中必須動(dòng)態(tài)為其分配內(nèi)存空間．vtkDataArray中的元素由組元(tuple)組成，每一個(gè)組元是一個(gè)結(jié)構(gòu)體，該結(jié)構(gòu)體由N個(gè)數(shù)據(jù)類型相同的分量組成．?dāng)?shù)據(jù)數(shù)組可以存儲(chǔ)的數(shù)據(jù)包括幾何數(shù)據(jù)和屬性數(shù)據(jù)(矢量、向量、張量)等．

開采沉陷預(yù)計(jì)的屬性數(shù)據(jù)有沉降值、水平移動(dòng)、水平變形、曲率變形和傾斜變形等，均為矢量，數(shù)據(jù)類型均為float，其中沉降值的方向垂直向下，可視為標(biāo)量．水平移動(dòng)、水平變形、曲率變形及傾斜變形為矢量，可以將其分解為X和Y方向的分量分別計(jì)算和提取，方法和沉降值的計(jì)算與提取類似．

3．3．2 開采沉陷預(yù)計(jì)后處理

開采沉陷預(yù)計(jì)后處理可用VTK現(xiàn)有的vtkMarchingContourFiltter、vtkCutter、vtkThresholdd等過(guò)濾器．vtkMarchingContourFiltter的輸入可以是繼承于vtkDataSet的任何數(shù)據(jù)類型，用于產(chǎn)生等值線或等值面;vtkCutter過(guò)濾器用于對(duì)數(shù)據(jù)集進(jìn)行切片;而vtkThresholdd過(guò)濾器的功能是提取那些滿足閾值范圍的cells．

結(jié)果后處理的可視化管線建立與3．1．3節(jié)類似．通過(guò)后處理，可從計(jì)算結(jié)果數(shù)據(jù)中提取等值線、生成結(jié)果云圖以及完成沉陷地表的屬性配色，使結(jié)果表達(dá)更豐富和形象．

4 應(yīng)用實(shí)例

某礦區(qū)地質(zhì)采礦條件:采煤方法為綜采放頂煤，頂板管理為全部垮落法;開采煤層為6煤，該煤層傾角為α=12°，開采深度為202～230 m，開采厚度8 m．該礦區(qū)概率積分法參數(shù)為:q=0．76，tanβ=2．0，b=0．36，θ=82．8°，s0=13m．

初始地表和對(duì)應(yīng)開采工作面如圖5所示，開采工作面從上到下分別為602、604、606、608、610、612工作面及6202、6204和6206工作面，開采沉陷預(yù)計(jì)結(jié)果如圖6，開采沉陷影響面積見表1．

表1 開采沉陷影響范圍Tab．1 Scope of mining subsidence

圖5 輸入地表和開采工作面Fig．5 The input DTM and working face model

圖6 開采沉陷預(yù)計(jì)結(jié)果Fig．6 The calculation result of mining subsidence

從圖6可知，沉陷預(yù)計(jì)最大值為5．22 m，與最大實(shí)測(cè)值5．07 m較接近;沉陷最大值發(fā)生在606工作面，預(yù)計(jì)傾斜最大值發(fā)生在沉陷等值線最密處，與實(shí)際情況相符;由表1可知，沉陷區(qū)面積與實(shí)際較吻合，計(jì)算具有較好的可靠性．

5 結(jié)語(yǔ)

基于VTK設(shè)計(jì)了vtkConformingDelaunay2D類，實(shí)現(xiàn)相容Delaunay三角化，得到高質(zhì)量的地表模型;根據(jù)沉陷預(yù)計(jì)模型，設(shè)計(jì)開采沉陷預(yù)計(jì)過(guò)濾器vtkSubsidenceFilter，實(shí)現(xiàn)沉陷計(jì)算;基于VTK實(shí)現(xiàn)沉陷預(yù)計(jì)結(jié)果的實(shí)時(shí)獲取和沉陷地表的動(dòng)態(tài)更新．該方法突破了以往VTK主要用于三維可視化的局限，通過(guò)擴(kuò)展VTK過(guò)濾器，拓寬了其專業(yè)應(yīng)用．研究成果可用于地下煤礦與金屬礦開采沉陷預(yù)計(jì)、開采沉陷的土地資源管理與復(fù)墾、以及自然崩落法放礦的礦巖接觸面模擬．應(yīng)用實(shí)例表明，這種應(yīng)用取得了較好效果．

［1］侯長(zhǎng)祥，馮濤，熊仁欽，等．礦床“三下一上”開采［M］．北京:煤炭工業(yè)出版社，2001．

［2］胡振琪．采煤沉陷地的土地資源管理與復(fù)墾［M］．北京:煤炭工業(yè)出版社，1996．

［3］趙經(jīng)徹，高延法，張懷新．兗州礦區(qū)開采沉陷控制的研究［J］．煤炭學(xué)報(bào)，1997，22(3):248－252．

［4］謝和平，周宏偉，王金安，等．FLAC3D在煤礦開采沉陷預(yù)測(cè)中的應(yīng)用及對(duì)比分析［J］．巖石力學(xué)與工程學(xué)報(bào)，1999，18(4):397－401．

［5］何國(guó)清，楊倫，凌賡娣，等．礦山開采沉陷學(xué)［M］．徐州:中國(guó)礦業(yè)大學(xué)出版社，1991．

［6］陳雨，張振文，張彥敏．基于概率積分法的地面沉陷災(zāi)害預(yù)測(cè)［J］．遼寧工程技術(shù)大學(xué)學(xué)報(bào)，2007，26(11):143－145．

［7］顧葉，宋振柏，張勝偉．基于概率積分法的開采沉陷預(yù)計(jì)研究［J］．山東理工大學(xué)學(xué)報(bào):自然科學(xué)版，2011，25(1):33－36．

［8］蔡美峰，李春雷，謝謨文，等．北洺河鐵礦開采沉陷預(yù)計(jì)及地表變形監(jiān)測(cè)與分析［J］．北京科技大學(xué)學(xué)報(bào)，2008，30(2):109－114．

［9］陳曉軍，閻世舉，吳軼群，等．基于VTK的口腔種植手術(shù)導(dǎo)航系統(tǒng)的圖像配準(zhǔn)［J］．計(jì)算機(jī)工程，2006，32(23):4－6．

［10］王延華，洪飛，吳恩華．基于VTK的醫(yī)學(xué)圖像處理子系統(tǒng)設(shè)計(jì)和實(shí)現(xiàn)［J］．計(jì)算機(jī)工程與應(yīng)用，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］劉偉寧，陳耀武，張磊．海底聲吶實(shí)時(shí)可視化顯示系統(tǒng)設(shè)計(jì)［J］．計(jì)算機(jī)工程，2010，36(23):249－251．

［13］畢林，王李管，陳建宏，等．基于VTK的礦體三維可視化研究與實(shí)現(xiàn)［J］．計(jì)算機(jī)工程與應(yīng)用，2008，44(10):78－81．

［14］歐耿鑫，陳喜，佘超，等．基于VTK的地下水?dāng)?shù)值模擬三維可視化開發(fā)應(yīng)用［J］．水文，2009，29(1):17－20．

［15］譚云蘭，李光耀，夏潔武，等．Delaunay三角網(wǎng)高效構(gòu)建及地形仿真應(yīng)用［J］．計(jì)算機(jī)工程，2012，38(22):287－290．

［16］陳學(xué)工，黃晶晶．Delaunay三角剖分中的約束邊嵌入算法［J］．計(jì)算機(jī)工程，2007，33(16):56－58．

［17］李智欣．Triangle與TetGen網(wǎng)格產(chǎn)生器于有限元素程式之使用［M］．臺(tái)南:國(guó)立成功大學(xué)，2004．

［18］李慶楊，王能超，易大義．?dāng)?shù)值分析［M］．武漢:華中科技大學(xué)出版社，2001．

［19］蔣和理．二重積分優(yōu)化梯形與辛卜生算法［J］．合肥工業(yè)大學(xué)學(xué)報(bào)，1987，9(6):34－36．