鄭朝治，譚榮建，陳曉軍

(1.昆明理工大學(xué)國(guó)土資源工程學(xué)院，云南昆明 650093;2.云南省文山州國(guó)土資源局土地整理中心，云南文山 663099)

礦體被采出之后，采空區(qū)上覆巖體的應(yīng)力平衡狀態(tài)受到破壞，引起地表變形、移動(dòng)的現(xiàn)象稱為開(kāi)采沉陷［1］。開(kāi)采沉陷是一種作用力集中、危害巨大的地質(zhì)災(zāi)害。它給礦區(qū)的生產(chǎn)、生活造成巨大的破壞，如鐵路架空、道路裂縫、耕地減少、房屋倒塌等。尤其是在人口密集、經(jīng)濟(jì)發(fā)達(dá)的城鎮(zhèn)及其附近，采礦與保護(hù)環(huán)境、建筑物安全的矛盾越來(lái)越突出。據(jù)不完全統(tǒng)計(jì)資料，目前我國(guó)“三下”壓煤約137.9億t［2］，部分礦區(qū)“三下”壓煤量十分巨大，已經(jīng)嚴(yán)重制約了礦區(qū)的發(fā)展。

開(kāi)采沉陷預(yù)測(cè)的基本目的是對(duì)開(kāi)采沉陷及其對(duì)環(huán)境產(chǎn)生的影響做出合理的評(píng)價(jià)，從而避免開(kāi)采對(duì)地表建筑物等造成嚴(yán)重破壞，有助于對(duì)沉陷情況的控制和礦區(qū)項(xiàng)目的決策。

20世紀(jì)90年代以來(lái)，隨著計(jì)算機(jī)技術(shù)的進(jìn)一步發(fā)展，以及地理信息系統(tǒng) (GIS)產(chǎn)業(yè)在全世界普及和發(fā)展，GIS突出的空間分析與三維可視化功能逐漸應(yīng)用于開(kāi)采沉陷問(wèn)題的研究。通過(guò)GIS軟件平臺(tái)的二次開(kāi)發(fā)和定制，可以較容易地實(shí)現(xiàn)預(yù)測(cè)和可視化的程序化，提高預(yù)測(cè)的速度和精度。作者對(duì)以C#為語(yǔ)言工具，以ArcGIS作為開(kāi)發(fā)主要技術(shù)平臺(tái)的開(kāi)采沉陷災(zāi)害預(yù)測(cè)系統(tǒng)做出設(shè)計(jì)與研究。該預(yù)測(cè)系統(tǒng)對(duì)沉陷數(shù)據(jù)進(jìn)行空間分析和可視化，有利于對(duì)地表移動(dòng)和破壞的更直觀和更深刻的理解。

1 開(kāi)采沉陷預(yù)測(cè)的數(shù)學(xué)模型

研究開(kāi)采沉陷作為一門(mén)學(xué)科始于20世紀(jì)20年代，國(guó)內(nèi)外眾多學(xué)者為此做出了許多嘗試和努力。提出的預(yù)測(cè)方法和沉陷模型大致可以分為經(jīng)驗(yàn)方法、剖面函數(shù)法、影響函數(shù)法、理論模型——物理模型。

1.1 概率積分法

概率積分法屬于影響函數(shù)法，因其移動(dòng)和變形預(yù)測(cè)公式中含有概率積分而得名。它認(rèn)為上覆巖土層是由裂隙所切割的大量極小碎塊體所組成，地下開(kāi)采引起的巖土與表層的移動(dòng)規(guī)律和隨機(jī)介質(zhì)模型中碎粒的移動(dòng)規(guī)律宏觀上相似。由于概率積分法算法簡(jiǎn)單可靠，易于計(jì)算機(jī)實(shí)現(xiàn)，目前已成為我國(guó)較為成熟、應(yīng)用最廣的煤礦開(kāi)采沉陷預(yù)測(cè)方法。

1.2 塌陷形態(tài)和剖面特點(diǎn)

開(kāi)采剛開(kāi)始時(shí)，隨著采空面積增大，移動(dòng)盆地的面積及地表最大下沉值也增大，盆地呈尖底的“碗狀”，此時(shí)地表的采動(dòng)影響稱為非充分采動(dòng)，如圖1中過(guò)程1和2。隨著開(kāi)采工作面的推進(jìn)，當(dāng)最大下沉值達(dá)到最大值并不再隨著采空面積增大而增大，但盆底仍呈尖底“碗狀”，此時(shí)，地表的采動(dòng)影響稱為充分采動(dòng)，如圖1中過(guò)程3。工作面再繼續(xù)推進(jìn)，下沉最大值不再增大，此時(shí)盆地呈平底的“盤(pán)狀”，此時(shí)的采動(dòng)影響稱為超充分采動(dòng)［3］，如圖1中過(guò)程4。半無(wú)限開(kāi)采，是指在斷面方向采空區(qū)為無(wú)限長(zhǎng)，它也是超充分采動(dòng)的一個(gè)理想狀態(tài)。

沉陷剖面的一側(cè)如圖2所示，其中r是主要影響半徑，銳角β為主要影響角，主要影響角正切tg β=H/r，H為采深。B是計(jì)算時(shí)認(rèn)為的煤壁位置，D為拐點(diǎn)，O為實(shí)際煤壁位置。O到B的距離就是拐點(diǎn)偏距S0的值，當(dāng)B在采空區(qū)那一側(cè)時(shí)S0為正。

圖1 開(kāi)采過(guò)程中塌陷盆地形態(tài)的變化

圖2 沉陷剖面示意圖

1.3 地表移動(dòng)變形計(jì)算

概率積分法中任意點(diǎn) (x，y)的下沉值 W(x，y)可用下式計(jì)算:

式中:m為開(kāi)采厚度;q為下沉系數(shù);a為煤層傾角;W0(x)W0(y)為走向l和傾向l有限開(kāi)采時(shí)主斷面地表下沉值;erf(x)是高斯誤差函數(shù):r、r1、r2分別為走向、下山、上山的主要影響半徑。

地表任意點(diǎn) (x，y)沿方向φ的傾斜值、曲率值、水平移動(dòng)值和水平變形值分別為:

式中:i(x，y，φ)，i0(x)，i0(y)是沿 φ方向、充分采動(dòng)時(shí)點(diǎn) (x，y)對(duì)應(yīng)走向斷面及傾向主斷面的傾斜值;K(x，y，φ)，K0(x)，K0(y)是沿φ方向、充分采動(dòng)時(shí)點(diǎn) (x，y)對(duì)應(yīng)走向主斷面及傾向主斷面的曲率值;U(x，y，φ)，U0(x)，U0(y)是沿 φ方向、充分采動(dòng)時(shí)點(diǎn)(x，y)對(duì)應(yīng)走向斷面及傾向主斷面的水平移動(dòng)值;ε (x，y，φ)，ε0(x)，ε0(y)是沿 φ方向、充分采動(dòng)時(shí)點(diǎn) (x，y)對(duì)應(yīng)走向斷面及傾向主斷面的水平變形值。

1.4 塌陷參數(shù)性質(zhì)和獲取方法

從概率積分法模型得出:走向方向達(dá)到充分采動(dòng)的條件時(shí)，工作面走向長(zhǎng)=l+S1+S2≥2r+S1+S2，即l≥2r;拐點(diǎn)處的下沉值為Wmax/2;拐點(diǎn)處的傾斜和水平移動(dòng)值達(dá)到最大;距拐點(diǎn)左右0.4r處的曲率和水平移動(dòng)值達(dá)到最大。

可以根據(jù)以上性質(zhì)通過(guò)實(shí)測(cè)來(lái)獲取拐點(diǎn)偏距S0和影響半徑r的值，也可以通過(guò)它們來(lái)驗(yàn)證數(shù)學(xué)模型的正確性。另外，一些現(xiàn)有的經(jīng)驗(yàn)數(shù)據(jù)也可以使用，例如:我國(guó)拐點(diǎn)偏距經(jīng)驗(yàn)值為0.05～0.3 H［1］。由于實(shí)際測(cè)量誤差，可能會(huì)造成多種途徑得出的參數(shù)數(shù)據(jù)不完全一致，根據(jù)實(shí)際情況選定合適的參數(shù)值，才能保證系統(tǒng)預(yù)測(cè)結(jié)果的精度。

2 系統(tǒng)設(shè)計(jì)

2.1 ArcEngine組件

ArcEngine組件是ArcObjects組件跨平臺(tái)應(yīng)用的核心集合，它提供了多種開(kāi)發(fā)的接口，可以適應(yīng)NET，Java，VB和C++等多種開(kāi)發(fā)環(huán)境。開(kāi)發(fā)者可以使用這些組件來(lái)開(kāi)發(fā)定制GIS和地圖應(yīng)用。目前，ArcEngine技術(shù)是開(kāi)發(fā)組件式GIS所采用的主流技術(shù)，它逐漸代替了組件式 GIS的其他技術(shù)，如MO、AO等。使用ArcEngine開(kāi)發(fā)出來(lái)的系統(tǒng)不需要安裝桌面版ArcGIS，只需安裝ArcEngine運(yùn)行環(huán)

境即可，這樣既節(jié)省了開(kāi)發(fā)成本，又節(jié)約了系統(tǒng)運(yùn)行所消耗的計(jì)算機(jī)系統(tǒng)資源［4］。

2.2 系統(tǒng)結(jié)構(gòu)和功能

經(jīng)過(guò)需求分析后，對(duì)軟件進(jìn)行初步的整體設(shè)計(jì)。然后逐一實(shí)現(xiàn)數(shù)據(jù)錄入、地表沉陷數(shù)據(jù)的計(jì)算、沉陷數(shù)據(jù)的圖形顯示、數(shù)據(jù)轉(zhuǎn)換接口等各功能部分。調(diào)試運(yùn)行通過(guò)以后，將它們集成組織起來(lái)，形成系統(tǒng)總體的功能模塊 (圖3)。這樣逐一進(jìn)行，最終完成系統(tǒng)的總體開(kāi)發(fā)。

圖3 系統(tǒng)的設(shè)計(jì)方案

參數(shù)輸入模塊。預(yù)測(cè)開(kāi)采沉陷，首先需要確定采區(qū)的基本情況，主要包括:開(kāi)采厚度m，下沉系數(shù)q，煤層傾角α，主要影響角正切tgβ，采深H，拐點(diǎn)偏距S0等。

沉陷數(shù)據(jù)計(jì)算模塊。將基本參數(shù)傳入C#編寫(xiě)的數(shù)學(xué)模型，計(jì)算出所需的高程、距離或者角度數(shù)據(jù)，并通過(guò)ArcEngine組件開(kāi)發(fā)的模塊生成三維空間數(shù)據(jù)。同時(shí)也要支持生成過(guò)程數(shù)據(jù)，以方便對(duì)系統(tǒng)計(jì)算過(guò)程更直觀的觀察與正確性檢驗(yàn)。

基于數(shù)據(jù)庫(kù)應(yīng)用的空間分析模塊。通過(guò)ArcSDE獲得地表地形數(shù)據(jù)，區(qū)縣行政區(qū)劃分?jǐn)?shù)據(jù)，房屋，道路，河流等已有數(shù)據(jù)，并與沉陷空間數(shù)據(jù)進(jìn)行疊加分析，得到沉陷區(qū)地表變化情況，建筑物受損情況，受損設(shè)施分布情況［5］。

成果的輸出模塊。系統(tǒng)的最終成果是多樣性的，包括興趣點(diǎn)單點(diǎn)沉陷查詢、任一點(diǎn)的移動(dòng)形變的變化圖、受災(zāi)區(qū)域的面數(shù)據(jù)等。還可以加入簡(jiǎn)單的經(jīng)濟(jì)損失模型以計(jì)算出開(kāi)采塌陷損失。

3 系統(tǒng)開(kāi)發(fā)的關(guān)鍵問(wèn)題

3.1 數(shù)值的使用方法

計(jì)算模塊根據(jù)開(kāi)采沉陷模型計(jì)算得到的是一系列的數(shù)據(jù)點(diǎn)陣，其中包括橫縱坐標(biāo)與點(diǎn)的下沉、曲率、水平移動(dòng)和變形等值。這些數(shù)據(jù)點(diǎn)都是離散的，不能構(gòu)成所需要的等值線結(jié)果。但是由于塌陷模型的特殊性，得到的沉陷面一定是關(guān)于地理坐標(biāo)連續(xù)且平滑的。所以，可以通過(guò)插值的方法用這些離散點(diǎn)構(gòu)建出很近似于真實(shí)情況的曲面。插值方法的實(shí)質(zhì)是在每2個(gè)已知點(diǎn)間插入多個(gè)點(diǎn)，然后將這些點(diǎn)用直線段依次連接起來(lái)［6］。當(dāng)加入了大量?jī)?nèi)插點(diǎn)時(shí)，用直線段連接后就可近似看作是一條光滑曲線。插值方法有許多種，其中3次樣條插值法已經(jīng)廣泛應(yīng)用于繪制等高線等解決實(shí)際問(wèn)題的多種工程軟件中，有較高的可靠性和較成熟的算法。所以本系統(tǒng)采用3次樣條作為模擬曲面的算法。圖4為沉陷數(shù)據(jù)模塊并插值后得到的沉陷盆地效果圖。

圖4 矩形開(kāi)采區(qū)的沉陷效果

3.2 地物的查詢分析

模擬預(yù)測(cè)礦區(qū)的受災(zāi)影響，需要首先建立該區(qū)域關(guān)于受災(zāi)參數(shù)的等值線，例如豎直位移，水平位移的等值線，從而得到已達(dá)到災(zāi)害程度的面狀要素。然后，通過(guò)進(jìn)一步空間查詢才能得到所影響的占有一定面積的地面單位，例如房屋，煙囪，道路等?；诳臻g位置的查詢，是根據(jù)要素與要素之間的空間關(guān)系進(jìn)行的。它用到 ArcEngine中的SpatialFilter組件類(lèi)，其中SpatialRel用于獲取或者設(shè)置查詢時(shí)的空間關(guān)系。設(shè)置空間關(guān)系為eariSpatialRelIntersects(相交)和 eariSpatialRel(包含)時(shí)，用受災(zāi)面對(duì)其他地面設(shè)施面進(jìn)行查詢，就可以得到受災(zāi)的地物。除此之外還有另一種受災(zāi)情況需要考慮，那就是管網(wǎng)類(lèi)的設(shè)施受損，比如公路、溝渠。此時(shí)必須得出受損的部分，僅僅查詢計(jì)算并不能滿足需求，必須對(duì)其進(jìn)行空間分析，具體方法為:用上面的方法查詢得到管線類(lèi)面或線，然后利用 ItopologicalOperator接口的Intersection方法得到與受災(zāi)區(qū)域的交集部分。

3.3 數(shù)據(jù)的疊加分析

由于空間數(shù)據(jù)三維問(wèn)題的復(fù)雜性，現(xiàn)有的軟件大多數(shù)停留在生成沉陷等值線這一步。為了使系統(tǒng)得到的模型更加直觀地表達(dá)，還需要建立塌陷前后3D模型的對(duì)比。使用 ArcEngine中 Scene加載DEM數(shù)據(jù)，使用 ApplyTin3Dproperties應(yīng)用3D屬性到遙感柵格圖像［7］，這樣就可以把影像覆蓋在建立好的高程模型上。

系統(tǒng)先將輸入的采空區(qū)數(shù)據(jù)進(jìn)行計(jì)算，得到離散的點(diǎn)陣列數(shù)據(jù)，進(jìn)行等高線繪制等操作。但是把豎直下陷數(shù)據(jù)和地表數(shù)據(jù)疊加時(shí)，必須重新自動(dòng)獲取離散點(diǎn)信息。這樣在3次樣條模擬的曲面上進(jìn)行離散點(diǎn)陣列的取樣，則會(huì)造成精度的損失。所以本系統(tǒng)使用最初計(jì)算模塊采用的點(diǎn)陣列，經(jīng)過(guò)坐標(biāo)轉(zhuǎn)換后直接與地表模型的對(duì)應(yīng)點(diǎn)高程相加，然后更新地表DEM數(shù)據(jù)。這樣解決數(shù)據(jù)疊加的問(wèn)題就更加準(zhǔn)確合理。

以上方法同樣適用于多個(gè)工作面引起沉陷的疊加，多個(gè)工作面在計(jì)算模塊中使用相同的離散點(diǎn)陣采樣，先把離散點(diǎn)陣上變形值相加，再與地面空間數(shù)據(jù)疊加以獲得相對(duì)較好的精度。

4 小結(jié)與展望

沉陷災(zāi)害預(yù)測(cè)是礦山災(zāi)害防治的重要內(nèi)容，其依據(jù)的基礎(chǔ)信息大多數(shù)與空間位置密切相關(guān)，使用GIS技術(shù)進(jìn)行評(píng)估，不但提高了評(píng)估的精度，還將沉陷地表形變、沉陷災(zāi)害評(píng)估結(jié)果可視化，可以更好地輔助開(kāi)采方案的設(shè)計(jì)、采動(dòng)損毀預(yù)防措施的制訂和實(shí)施。

本研究方案以ArcEngine組件為工具使用GIS技術(shù)結(jié)合礦區(qū)開(kāi)采沉陷的專(zhuān)業(yè)預(yù)計(jì)模型，采用先進(jìn)的數(shù)據(jù)處理手段、可視化的表達(dá)與多元化的分析方法，一體化地對(duì)礦區(qū)沉陷災(zāi)害進(jìn)行預(yù)測(cè)預(yù)報(bào)［8］。但軟件的設(shè)計(jì)還存在著不足，如沉陷模型較單一，只適用于比較理想化的工作面等。

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

［2］ 譚志祥，鄧喀中.建筑物下采煤理論與實(shí)踐 ［M］.徐州:中國(guó)礦業(yè)大學(xué)出版社，2006:485-488.

［3］ 徐永圻.煤礦開(kāi)采學(xué)［M］.徐州:中國(guó)礦業(yè)大學(xué)出版社，1999:417-418.

［4］ 邱洪鋼，張青蓮，陸紹強(qiáng).ArcEngine從入門(mén)到精通 ［M］.北京:人民郵電出版社，2010.

［5］ 張麗軍，劉亮，潘家偉.基于GIS的地震災(zāi)害損失快速評(píng)估系統(tǒng)開(kāi)發(fā)［J］.科技創(chuàng)新導(dǎo)報(bào)，2009(25):128-129.

［6］ 李慶揚(yáng)，王能超，易大義.數(shù)值分析［M］.武漢:華中科技大學(xué)出版社，2006:22-23.

［7］ 劉興權(quán)，盧賽飛.基于ArcEngine的地物三維可視化實(shí)現(xiàn)［J］.地理空間信息，2008，6(6):5-8.

［8］ 王來(lái)強(qiáng)，姚頑強(qiáng).基于GIS對(duì)煤礦開(kāi)采沉陷系統(tǒng)的研究［J］.技術(shù)與創(chuàng)新管理，2009，30(1):96-98.