蔡音飛，李曉靜，鄧偉男

(1.太原理工大學 礦業(yè)工程學院，山西 太原 030024；2.洛林大學 南錫高等礦業(yè)學院 ，法國 南錫 54000；

3.山西農(nóng)業(yè)大學 資源與環(huán)境學院，山西 太谷 030801；4.天地科技股份有限公司 開采設計事業(yè)部，北京 100013)

?

基于Mathematica的顧及地形的開采沉陷三維可視化

蔡音飛1，2，李曉靜3，鄧偉男4

(1.太原理工大學 礦業(yè)工程學院，山西 太原 030024；2.洛林大學 南錫高等礦業(yè)學院 ，法國 南錫 54000；

3.山西農(nóng)業(yè)大學 資源與環(huán)境學院，山西 太谷 030801；4.天地科技股份有限公司 開采設計事業(yè)部，北京 100013)

[摘要]利用Mathematica軟件，考慮原始地表數(shù)據(jù)和由顧及地形的開采沉陷預計方法計算得到的沉陷數(shù)據(jù)，實現(xiàn)了原始地表、沉陷和受沉陷影響的地表的三維可視化。顧及地表的開采沉陷預計需使用2個偏斜影響函數(shù)分別描述單元開采的下沉和水平移動，進而按經(jīng)典方法進行積分。應用Mathematica的插值函數(shù)和繪圖函數(shù)，繪制所需的各種三維圖形，可直觀反映沉陷后的地表形態(tài)。

[關鍵詞]Mathematica；地形影響；影響函數(shù)法；插值；三維可視化

采礦引起的地表沉陷對礦區(qū)的土地、建筑、道路等的功能和安全造成了不同程度的影響[1-2]。很多礦區(qū)的地表是非水平的，結合地表形態(tài)因素，以三維圖示的方法表現(xiàn)沉陷和沉陷后的地形，可以更直觀、準確地對沉陷區(qū)土地、地物進行損毀分析。單獨沉陷數(shù)據(jù)的三維顯示在很多研究中都有報告，但只能表現(xiàn)礦區(qū)地表為水平面時的開采沉陷狀況，而無法顧及礦區(qū)原有的地表特征，因此尚不能夠真實反映開采后的地表情況[3]。有些研究工作已實現(xiàn)了沉陷影響下的礦區(qū)地形三維可視化，或采用DEM法編程實現(xiàn)[4-5]，或使用GIS軟件[6-11]、南方CASS[12]等配合沉陷數(shù)據(jù)實現(xiàn)。這些方法雖然在原始地形上疊加了沉陷數(shù)據(jù)以實現(xiàn)受開采影響的地形三維顯示，但是在計算沉陷時均未能考慮原始的非水平地表對沉陷的影響，并且沉陷預計與其后的三維顯示成圖均不是在同一軟件中實現(xiàn)的，操作和數(shù)據(jù)轉換有諸多不便。

本文使用考慮地形變化的開采沉陷預計程序，計算非水平地表條件下的沉陷，然后與原始地形數(shù)據(jù)相疊加，獲得受開采影響后的地表數(shù)據(jù)，并實現(xiàn)三維顯示，為礦山開采前的設計、規(guī)劃，開采后的評估、復墾提供可視化支持和依據(jù)。

本文沉陷預計和數(shù)據(jù)疊加均在Mathematica中完成。Mathematica是由沃爾夫勒姆研究公司開發(fā)的一款被廣泛使用的科學計算軟件，擁有強大的數(shù)值計算和符號運算能力，是科學計算的完全集成環(huán)境、交互式的問題解決系統(tǒng)[13]，其提供的高等數(shù)學函數(shù)可以直接用于本研究的開采沉陷預計程序的開發(fā)，極大地簡化了編程過程；其提供的插值函數(shù)和繪圖函數(shù)也為三維顯示輸出提供了便利。

1顧及地表形態(tài)的開采沉陷預計方法

理論上，水平地表條件下，沉陷與地形疊加后的形態(tài)與僅考慮沉陷不疊加地表的形態(tài)是完全一致的。本文研究對象是非水平地表條件下的開采沉陷三維可視化，首先需將地表變化的影響納入開采沉陷預計計算中。

本文使用的開采沉陷預計程序是基于影響函數(shù)法[14-15]，在Mathematica環(huán)境下開發(fā)的。筆者在文獻[16]中詳細說明了考慮地形因素的開采沉陷預計程序的基本原理，即采用偏斜正態(tài)分布的概率密度函數(shù)，如式(1)、式(2)所示。作為影響函數(shù)，用以描述單元開采在地表傾向方向上的下沉和水平移動，按照經(jīng)典影響函數(shù)法，采用積分的方法求取任意地表點的移動變形量。

(1)

(2)

(3)

(4)

式(1)～式(4)中：G(x)為正態(tài)分布(高斯分布)的概率密度函數(shù)；erfc(x)為高斯誤差函數(shù)，用于偏斜G(x)；G(x)和erfc(x)的積即為偏斜正太分布的概率密度函數(shù)；ifv(x)為下沉在地表傾向方向上的影響函數(shù)；ifh(x)為水平移動在地表傾向方向上的影響函數(shù)；sm，α，μ，σ為地形相關的修正系數(shù)，對于ifv(x)，ifh(x)其取值是不同的，其中，sm用于修正單元開采下沉和水平移動的數(shù)值大小、α用于修正他們的偏斜度、μ用于修正他們在地表傾向方向上相對于開采單元的位置、σ用于修正他們的影響范圍。

相對于經(jīng)典影響函數(shù)法，本方法最主要的特點為：

(1)使用了2個影響函數(shù)分別描述單元開采的下沉和水平移動。研究發(fā)現(xiàn)，在非水平地表條件下，水平移動不再和傾斜(下沉的導數(shù))相似，即水平移動無法再由下沉計算，因此需要定義2個影響函數(shù)。其他的地表移動變形量(傾斜、水平變形、曲率)仍可用經(jīng)典方法計算。

(2)地表高程數(shù)據(jù)需要作為已知數(shù)據(jù)輸入計算。上述2個影響函數(shù)中的地形影響修正系數(shù)都可以表示為地表相對礦體高程和地表傾角的函數(shù)。實際上，地表傾角也可以由地表高程求導得到，所以對于本方法而言，只有地表高程數(shù)據(jù)(更準確地說是地表相對礦體高程數(shù)據(jù))是較經(jīng)典方法新增的輸入數(shù)據(jù)。

通過使用按照上文原理開發(fā)的開采沉陷預計程序，可以在沉陷計算中考慮非水平地表的影響，相對于原來的預計方法，能夠得到相對更精確的預計結果。

2受沉陷影響的地表三維可視化

Mathematica提供了一個插值函數(shù)(Interpolation)，適用于任意維數(shù)數(shù)據(jù)的插值(本文涉及的地形數(shù)據(jù)、沉陷數(shù)據(jù)均為三維)，插值的方法可選擇樣條曲線(Spline)插值法或埃爾米特(Hermite)插值法，得到的結果是一個函數(shù)，可以和其他任何數(shù)學函數(shù)一樣使用。應用該插值函數(shù)，配合Mathematica提供的三維繪圖函數(shù)(Plot3D)可以方便地繪制地表和沉陷的三維圖形，包括受沉陷影響后的地表形態(tài)。

本節(jié)結合法國洛林地區(qū)某礦實例[16]進行說明。洛林曾是法國主要的礦產(chǎn)區(qū)，當?shù)卮蠖嗖捎梅恐ㄩ_采，很多礦山在停采多年后面臨礦柱逐漸失穩(wěn)，地表沉陷頻發(fā)的問題。本文研究的2個開采區(qū)塊近年來有失穩(wěn)的跡象，其拐點坐標見表1(文中圖表所用坐標系均為RGF93 / Lambert-93，為便于表述，所有平面坐標均統(tǒng)一減去一定值，經(jīng)過現(xiàn)場實測，這2個區(qū)塊造成的最大地表下沉值約為0.5m，影響角約為30°。

表1　失穩(wěn)的開采區(qū)塊拐點坐標

2.1原始地表

任何途徑獲得的原始地表數(shù)據(jù)均可以用三維數(shù)組表示：surfdata = {{x1，y1，z1}，{x2，y2，z2}，…,{xn，yn，zn}}，其中x,y為地表點的平面坐標，z為對應點的高程。開采區(qū)塊上方的若干離散地表點高程數(shù)據(jù)已知，這些點的位置繪于圖1(圖中，每個點相連的豎線底端位置是該點的平面位置，豎線長度為相對最低點的高程)，其高程約介于+330m到+390m之間，最大高差約60m。

圖1　有高程數(shù)據(jù)的地表點位置示意

以x，y為自變量，使用樣條曲線插值法，地表點高程z的插值函數(shù)(此處記為SurfZ[x，y])可以用以下命令在Mathematica中求?。?/p>

SurfZ[x，y]= Interpolation [surfdata，Method → “Spline”，InterpolationOrder → 1]

“InterpolationOrder → 1”表示不進行平滑處理，該值越大，插值后的函數(shù)圖形越光滑。

在Mathematica中使用Plot3D函數(shù)(僅需提供三維函數(shù)作為輸入數(shù)據(jù)即可輸出圖形)，可將原始地表函數(shù)SurfZ(x,y)輸出為三維圖形：

Plot3D [SurfZ[x，y]，{x，0，600}，{y，0，600}，UserSettings]

“{x，0，600}”和“{y，0，600}”用于指定輸出范圍；“UserSettings”是圖形顯示相關附加設置的集合(如：顏色、字體、坐標軸、比例、觀察點等)。

輸出的原始地表三維圖見圖2。需要說明的是(對圖3，圖4亦同，不再贅述)：本文將地表的顏色設置為高程的函數(shù)并用灰度區(qū)分，實際上，Mathematica提供了很多顏色函數(shù)，支持彩色輸出；觀察點位置可以設置，或者直接在程序中拖動圖形進行三維旋轉；傳統(tǒng)的平面等值線、色暈圖與三維輸出方法類似，因非本文內(nèi)容，在此不展開說明；圖中地表的網(wǎng)格線是為了更清楚地顯示，實際操作中并不需要求網(wǎng)格點處的函數(shù)值；為了便于將開采區(qū)塊與地表等一起顯示，圖中開采區(qū)塊的高程是經(jīng)過調(diào)整的，各點實際高程以表1為準。

圖2　原始地表三維示意

2.2開采沉陷計算

應用上文考慮地表形態(tài)的開采沉陷預計方法(已在Mathematica中編程實現(xiàn))，計算由開采區(qū)塊1，2礦柱失穩(wěn)造成的沉陷，獲得下沉、水平移動數(shù)據(jù)。與處理地表數(shù)據(jù)類似，可以在Mathematica下，用插值函數(shù)算得下沉、x方向水平移動、y方向水平移動的函數(shù)，記為SubsV[x，y]、SubsHx[x，y]、SubsHy[x，y]。以下沉為例，其三維圖如圖3所示。

圖3　下沉三維示意

在Mathematica中，插值函數(shù)與一般初等函數(shù)一樣，可以通過賦予自變量值得到函數(shù)值；也可以進行四則運算、求導、積分等。因此用下沉、水平移動的函數(shù)求偏導可以求得相應方向上的傾斜、水平變形、曲率。

2.3受沉陷影響的地表

在求取受下沉影響的地表時，無需如傳統(tǒng)用列表的方法一樣考慮數(shù)據(jù)點位置對應問題，直接將原始地表函數(shù)SurfZ[x，y]與下沉函數(shù)SubsV[x，y]相加(下沉記為負值)，即可得到考慮下沉的地表函數(shù)。水平移動與原始地表疊加的方法類似，但稍微復雜，需將SubsHx[x，y]，SubsHy[x，y]與對應的平面坐標相加。由于沉陷量(下沉0.5m左右)與地表高差(60m左右)相比較小，按照實際沉陷量與地表疊加后與原始地表對比不明顯。本文將沉陷量放大20倍后，疊加沉陷(下沉和水平移動)的地表三維示意見圖4。

圖4　受沉陷影響的地表三維示意

2.4數(shù)據(jù)、圖形輸出

Mathematica可以將數(shù)據(jù)以各種自定義的列表形式輸出到Text或者Excel等常用軟件，與其他程序配合使用無障礙；圖形也可以以圖片格式(常用均可)或者矢量格式(包括eps，svg，dxf等)輸出，便于進一步處理。

本文為說明Mathematica在圖形顯示方面的便利和優(yōu)越性，所有圖片均僅在Mathematica中通過調(diào)整顯示參數(shù)(UserSettings)直接生成，未在其他軟件中處理。

3結論

(1)使用考慮地形影響的開采沉陷疊加原始地表，實現(xiàn)了受沉陷影響的地表三維可視化。

(2)沉陷計算和受沉陷影響的地表三維可視化均在同一個軟件中實現(xiàn)，操作更方便，無需考慮數(shù)據(jù)轉換問題。

(3)使用Mathematica提供的高等數(shù)學函數(shù)、插值函數(shù)、繪圖函數(shù)，編程更簡練，效率更高。

[參考文獻]

[1]鄧偉男，張華興，徐乃忠.開采損害影響范圍的多方法比較鑒定模式及其應用[J].煤礦開采，2013，18(2)：105-107.

[2]Li Xiao-jing，Hu Zhen-qi，Li Shuang-cheng，et al.Anomalies ofmountainousmining paddy in western China[J].Soil and Tillage Research，2015(145)：10-19.

[3]肖武，胡振琪，李太啟.采區(qū)地表動態(tài)沉陷模擬與復墾耕地率分析[J].煤炭科學技術，2013，41(8)：126-128.

[4]柴華彬，鄒友峰，劉景艷.DTM在開采沉陷可視化預計中的應用[J].遼寧工程技術大學學報(自然科學版)，2004，23(2)：171-174.

[5]陳秋計，胡振琪，劉昌華，等.DEM在礦區(qū)土地復墾中的應用研究[J].金屬礦山，2006 (2)：67-68.

[6]戴華陽.基于傾角變化的開采沉陷模型及其GIS可視化應用研究[J].巖石力學與工程學報，2002，21(1)：148-148.

[7]劉立民，劉漢龍，連傳杰.基于GIS的礦山塌陷損害評價系統(tǒng)及可視化方法[J].防災減災工程學報，2003，23(1)：69-73.

[8]孟峰.礦區(qū)地表沉陷預計三維可視化研究[J].陜西煤炭，2006 (2)：10-12.

[9]王京衛(wèi)，欒紅，汝續(xù)偉.顧及礦區(qū)地貌特征的開采沉陷三維可視化實現(xiàn)方法[J].山東科技大學學報(自然科學版)，2007，26(1)：8-11.

[10]Dai Hua-yang，Ren Li-yan，Wang meng，et al.Water distribution extracted frommining subsidence area using Kriging interpolation algorithm[J].Transactions of Nonferrousmetals Society of China，2011(21)：723-726.

[11]盧志剛，劉興權，唐義宏.基于GIS和FLAC3D礦山地表沉陷可視化[J].現(xiàn)代礦業(yè)，2012 (11)：44-46.

[12]易四海，戴華陽，廉旭剛，等.顧及地表沉陷的礦區(qū)地貌可視化實現(xiàn)方法[J].湖南科技大學學報(自然科學版)，2008，23(4)：81-84.

[13]Wolfram Stephen.Themathematica book[M].Wolframmedia，Incorporated，1996.

[14]National Coal Board.Subsidence engineers’handbook[M]. London：National Coal Board，1975.

[15]何國清，楊倫，凌賡娣，等.礦山開采沉陷學[M].徐州：中國礦業(yè)大學出版社，1991.

[16]Cai Yin-fei，Verdel Thierry，Deck Olivier.On the topography influence on subsidence due to horizontal undergroundmining using the influence functionmethod[J].Computers and Geotechnics，2014，61(6)：328-340.

[責任編輯：徐乃忠]

3-D Visualization of Mining Subsidence Considering Terrain Based on Mathematica

CAI Yin-fei1，2,LI Xiao-jing3,DENG Wei-nan4

(1.College of Mining Engineering,Taiyuan University of Technology,Taiyuan 030024,China;2.Nancy Higher Mining Institute,Lorraine University,Nancy 54000,France;3.College of Resources and Environment,Shanxi Agricultral University,Taigu 030801;4.Coal Mining & Designing Department,Tiandi Science & Technology Co.,Ltd.,Beijing 100013,China)

Abstract:Applying mathematica software,3-d visualization of original surface,subsidence,and surface influenced by subsidence was realized with original surface data and subsidence data calculated by mining subsidence prediction method. The mining subsidence prediction considering terrain need to use deflection influence functions respectively describing sink and horizontal displacement of elements,and farther finish integral with classic method. Using interpolation function and drawing function,various 3-d graphs were drawn which could directly reflect surface configuration after subsidence.

Key words:Mathematica;terrain influence;influence function method;interpolation;3-d visualization

[中圖分類號]TD327

[文獻標識碼]A

[文章編號]1006-6225(2016)01-0080-04

[作者簡介]蔡音飛(1983-)，男，浙江嘉興人，博士，主要研究方向為“三下”采煤技術及其理論的研究與應用。

[收稿日期]2015-07-19

[DOI]10.13532/j.cnki.cn11-3677/td.2016.01.022

[引用格式]蔡音飛，李曉靜，鄧偉男.基于Mathematica的顧及地形的開采沉陷三維可視化[J].煤礦開采，2016，21(1)：80-83.