蔡音飛，李曉靜，鄧偉男

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

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

?

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

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

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

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

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

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

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

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

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

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

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

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

(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)為水平移動(dòng)在地表傾向方向上的影響函數(shù)；sm，α，μ，σ為地形相關(guān)的修正系數(shù)，對于ifv(x)，ifh(x)其取值是不同的，其中，sm用于修正單元開采下沉和水平移動(dòng)的數(shù)值大小、α用于修正他們的偏斜度、μ用于修正他們在地表傾向方向上相對于開采單元的位置、σ用于修正他們的影響范圍。

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

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

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

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

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

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

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

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

2.1原始地表

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

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

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

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

“InterpolationOrder → 1”表示不進(jìn)行平滑處理，該值越大，插值后的函數(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”是圖形顯示相關(guān)附加設(shè)置的集合(如：顏色、字體、坐標(biāo)軸、比例、觀察點(diǎn)等)。

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

圖2　原始地表三維示意

2.2開采沉陷計(jì)算

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

圖3　下沉三維示意

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

2.3受沉陷影響的地表

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

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

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

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

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

3結(jié)論

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

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

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

[參考文獻(xiàn)]

[1]鄧偉男，張華興，徐乃忠.開采損害影響范圍的多方法比較鑒定模式及其應(yīng)用[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ū)地表動(dòng)態(tài)沉陷模擬與復(fù)墾耕地率分析[J].煤炭科學(xué)技術(shù)，2013，41(8)：126-128.

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

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

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

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

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

[9]王京衛(wèi)，欒紅，汝續(xù)偉.顧及礦區(qū)地貌特征的開采沉陷三維可視化實(shí)現(xiàn)方法[J].山東科技大學(xué)學(xué)報(bào)(自然科學(xué)版)，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]盧志剛，劉興權(quán)，唐義宏.基于GIS和FLAC3D礦山地表沉陷可視化[J].現(xiàn)代礦業(yè)，2012 (11)：44-46.

[12]易四海，戴華陽，廉旭剛，等.顧及地表沉陷的礦區(qū)地貌可視化實(shí)現(xiàn)方法[J].湖南科技大學(xué)學(xué)報(bào)(自然科學(xué)版)，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]何國清，楊倫，凌賡娣，等.礦山開采沉陷學(xué)[M].徐州：中國礦業(yè)大學(xué)出版社，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.

[責(zé)任編輯：徐乃忠]

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

[文獻(xiàn)標(biāo)識碼]A

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

[作者簡介]蔡音飛(1983-)，男，浙江嘉興人，博士，主要研究方向?yàn)椤叭隆辈擅杭夹g(shù)及其理論的研究與應(yīng)用。

[收稿日期]2015-07-19

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

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