王 璐,盧小平,李國利,李小雷,康愛峰

(1. 河南理工大學(xué)自然資源部礦山時空信息與生態(tài)修復(fù)重點實驗室,河南 焦作454003;2. 鶴壁恒源礦業(yè)集團有限公司,河南 鶴壁458030)

傳統(tǒng)二維GIS技術(shù)難以描述復(fù)雜環(huán)境地質(zhì)體的三維空間信息,而三維GIS技術(shù)可以更好地解決這一問題[1-2]。近年來,隨著智能化進程的加快和地下空間利用需求的激增,地質(zhì)體三維建模成為三維GIS領(lǐng)域的研究熱點。傳統(tǒng)的地質(zhì)體三維建模中,鉆孔數(shù)據(jù)是最主要的數(shù)據(jù)來源[3]。但是鉆孔數(shù)據(jù)往往成本較高、數(shù)量較少,用少量的鉆孔數(shù)據(jù)表達復(fù)雜的地質(zhì)體內(nèi)部結(jié)構(gòu),往往需要對鉆孔數(shù)據(jù)進行插值加密,生成虛擬鉆孔,以提高建模的精度[4]。

地質(zhì)體三維建模是礦產(chǎn)儲量估算的基礎(chǔ),精確的地質(zhì)體三維模型可以提高儲量估算的精度,在礦區(qū)的生產(chǎn)規(guī)劃和監(jiān)測方面起到指導(dǎo)作用。地質(zhì)體三維建模插值方法大體分為兩類:一類是確定性插值方法,另一類是地統(tǒng)計插值方法。文獻[5]利用局部多項式插值算法進行插值擬合,計算了礦山現(xiàn)開采范圍內(nèi)的土石方量值。文獻[6—8]采用普通克里金和反距離加權(quán)插值法對礦體進行了儲量估算。文獻[9—11]利用克里金插值進行地質(zhì)體三維建模,進而對礦產(chǎn)儲量估算。文獻[12]對比了克里金插值、樣條插值及反距離加權(quán)插值,將礦體表面模型進行疊加,采用礦量計算公式,得到弓長嶺獨木礦區(qū)鐵礦的總儲量。但這些方法大都忽視了地質(zhì)原理對插值方法的要求,對地質(zhì)體采用單一的插值方法難以取得最優(yōu)的插值結(jié)果。

針對以上問題,本文結(jié)合鶴壁市二道莊露天礦的實際情況,提出一種混合插值算法。一方面,對鉆孔分層數(shù)據(jù)進行空間自相關(guān)分析,針對具有空間自相關(guān)性的地層采用克里金、規(guī)則樣條、張力樣條、反距離加權(quán)進行插值,對不具有空間相關(guān)性的地層,則利用規(guī)則樣條、張力樣條、反距離加權(quán)進行插值,通過對比分析選擇出適用于每個地層的插值方法;另一方面,利用Dynamo進行可視化編程,建立基于鉆孔數(shù)據(jù)的地質(zhì)體三維模型,并基于地質(zhì)體三維模型,利用SuperMap進行二次開發(fā),建立礦產(chǎn)儲量估算系統(tǒng)。

1 研究方法

本文研究方法的具體方法流程如圖1所示。

圖1 方法流程

1.1 插值算法

空間插值是將空間連續(xù)的點轉(zhuǎn)化為連續(xù)的數(shù)據(jù)曲面。根據(jù)插值形式的不同,可分為空間外插和空間內(nèi)插[13];根據(jù)插值方法的不同,可分為確定性方法和地質(zhì)統(tǒng)計學(xué)方法[14]。確定性插值方法是基于樣本點之間的相似程度創(chuàng)建一個插值擬合曲面,常見的有反距離加權(quán)法(IDW)、徑向基函數(shù)法(RBF)等。地質(zhì)統(tǒng)計學(xué)方法是根據(jù)樣本點之間的統(tǒng)計規(guī)律與空間自相關(guān)性確定待測點的值,常見的有克里金(Kriging)插值法等。因此,可以根據(jù)已知點之間是否具有空間自相關(guān)性選擇插值方法,從而達到更好的插值效果。

1.1.1 反距離加權(quán)法(IDW)

反距離加權(quán)法的表達式[15]為

(1)

式中,f(x,y)為待測點的屬性值;n為已知點的個數(shù);zi為第i個已知點的屬性值,i=1, 2, …,n;di為(x,y)到(xi,yi)的歐式距離;p為冪指數(shù),通常取值為1～3,但一般p為2時獲得的效果更好。

1.1.2 徑向基函數(shù)法

徑向基函數(shù)較適合處理空間各向同性的問題。常見的徑向基函數(shù)有薄板樣條函數(shù)、張力樣條函數(shù)、規(guī)則樣條函數(shù)、高次曲面樣條函數(shù)及反高次曲面樣條函數(shù)[16],其表達式為

(2)

式中,ai為待定系數(shù);di為歐式距離;φ為一種徑向基函數(shù)。

1.1.3 Kriging法

Kriging算法包括普通克里金、簡單克里金、泛克里金、協(xié)同克里金等[17]。半變異函數(shù)是克里金插值的一種模型,用于評估樣本點之間的空間自相關(guān)性[18]。常用的半變異函數(shù)模型有球面模型、指數(shù)模型、高斯模型等。克里金插值算法的數(shù)學(xué)表達式為

(3)

式中,f(x,y)為待測點的屬性值;zi為第i個已知點的屬性值;λi為權(quán)重系數(shù)。能夠滿足點(xo,yo)處的插值f(xo,yo)與真實值zo差值最優(yōu)系數(shù)為

(4)

(5)

同時,滿足無偏估計條件為

(6)

(7)

1.1.4 空間自相關(guān)分析

在對地質(zhì)體插值時,利用變異函數(shù)作為判斷鉆孔數(shù)據(jù)是否空間自相關(guān)的工具[19]。變異函數(shù)主要包括變程、基臺值、塊金值、偏基臺值,如圖2所示。通過繪制半變異函數(shù)云圖表示自相關(guān)分析結(jié)果。半變異云圖反映了采樣點與相鄰采樣點之間的空間關(guān)系,距離越大,半變異值越大,空間相關(guān)性也越好。

圖2 半變異函數(shù)結(jié)構(gòu)

1.1.5 交叉驗證與精度評價

為減少主觀判斷對差值結(jié)果客觀性的影響,選擇留一法進行交叉驗證,并采用平均誤差(ME)和均方根誤差(RMSE)對插值結(jié)果進行精度評價[20]。

平均誤差是預(yù)測值與實際值之間插值的平均值,計算公式為

(8)

均方根誤差是用于描述預(yù)測值與實際值之間的偏離程度,計算公式為

(9)

由于自然界的地質(zhì)體形狀各異,任何插值方法都無法適應(yīng)所有的地質(zhì)體,因此需要對研究區(qū)進行具體分析,從而選出最合適的插值方法。本文首先對鉆孔分層數(shù)據(jù)進行自相關(guān)分析,然后選擇相應(yīng)的插值方法。

1.2 地質(zhì)體三維建模

建筑信息模型(BIM)用于解決復(fù)雜的三維地質(zhì)體的可視化。文獻[21—22]提出了基于BIM信息模型的三維滑坡地質(zhì)災(zāi)害預(yù)測;文獻[23]采用BIM軟件,利用鉆孔數(shù)據(jù)和地質(zhì)剖面數(shù)據(jù)建立了地質(zhì)體三維模型。

1.2.1 Revit和Dynamo可視化編程

近年來,越來越多的人基于BIM技術(shù)在地質(zhì)體三維建模方面進行研究[24-27]。Revit軟件提供了強大的族功能,Dynamo是Revit的一個可視化編程插件,具有強大的節(jié)點庫,包含運算、分析、圖像等功能,節(jié)點可以通過邏輯關(guān)系進行連接實現(xiàn)可視化,也可以采用Python等編程語言達到更多目的[28]。

1.2.2 建模流程

由于研究區(qū)范圍小、地質(zhì)構(gòu)造相對簡單,不涉及裂縫和倒置地層,因此本文基于鉆孔數(shù)據(jù)使用Revit進行地質(zhì)體三維建模。該方法首先通過對分層鉆孔數(shù)據(jù)進行插值處理,選出每層數(shù)據(jù)最佳的插值方法;然后根據(jù)插值結(jié)果生成虛擬鉆孔;最后基于Dynamo生成地質(zhì)模型。具體建模流程如圖3所示。

圖3 建模流程

2 礦山儲量估算系統(tǒng)設(shè)計

礦產(chǎn)儲量估算系統(tǒng)采用C/S的結(jié)構(gòu)模式,在Visual Studio 2013 環(huán)境下,運用C#語言,在SuperMap平臺上開發(fā),將礦山地質(zhì)體三維模型與礦產(chǎn)儲量估算系統(tǒng)進行融合。系統(tǒng)包含數(shù)據(jù)管理、礦山二三維顯示、儲量估算3個模塊。數(shù)據(jù)管理模塊包括礦區(qū)的DOM、DEM及鉆孔數(shù)據(jù)的管理;礦山二三維顯示模塊包括礦體地質(zhì)體三維模型、三維地形模型、矢量模型的顯示,以及對模型的交互瀏覽、平移、縮放、旋轉(zhuǎn)等操作。儲量估算模塊包括對開采量的計算及報表輸出等。系統(tǒng)功能結(jié)構(gòu)如圖4所示。

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

3 試驗與分析

3.1 研究區(qū)概況

研究區(qū)為鶴壁市二道莊礦區(qū),面積約為5.2 km2, 海拔高度為150～440 m,如圖5所示。

圖5 研究區(qū)

3.2 空間數(shù)據(jù)采集與處理

利用無人機采集空間數(shù)據(jù),經(jīng)內(nèi)業(yè)處理得到礦區(qū)的DEM與DOM。共收集40個工程地質(zhì)鉆孔,鉆孔深度為40.89～121.56 m。原始鉆孔柱狀圖中包含鉆孔編號、孔口坐標、開孔日期、終孔日期、終孔深度、分層情況、巖性描述等信息。鉆孔數(shù)據(jù)坐標采用CGCS2000。得到鉆孔數(shù)據(jù)空間分布如圖6所示。

圖6 鉆孔分布

3.3 三維地質(zhì)體構(gòu)建

3.3.1 鉆孔插值及精度檢驗

鉆孔數(shù)據(jù)按照巖性分為白云質(zhì)灰?guī)r、鮞狀白云巖、白云巖、鮞狀白云巖、灰?guī)r共5層,分別用Z0、Z1、Z2、Z3、Z4表示,并導(dǎo)入數(shù)據(jù)庫中統(tǒng)一管理。計算各地層數(shù)據(jù)的變異函數(shù),并繪制變異函數(shù)云圖如圖7所示,據(jù)此分析各層數(shù)據(jù)的空間自相關(guān)性。

圖7 各地層變異函數(shù)云圖

由圖7可知,Z0、Z2、Z3地層的云圖點較散亂,變異函數(shù)值隨已知點之間距離的增大呈先增大、后減小的趨勢,不符合空間自相關(guān)的云圖特點,即不具備空間自相關(guān)性;Z1、Z4的變異函數(shù)值隨距離的增大有明顯的增大趨勢,符合空間自相關(guān)的特性。因此,對Z1、Z4采用克里金插值、徑向基函數(shù)插值及反距離加權(quán)插值;對Z0、Z2、Z3則采用徑向基函數(shù)插值和反距離加權(quán)插值。其中,克里金插值選取普通克里金插值,變異函數(shù)選擇球面函數(shù),徑向基函數(shù)插值選取規(guī)則樣條函數(shù)和張力樣條函數(shù)兩種插值方法。插值結(jié)果對比見表1、表2。

表1 Z1、Z4插值結(jié)果

表2 Z0、Z2、Z3插值結(jié)果

由表1可知,Z1、Z4地層較平坦,高低點起伏不大,整體呈上(北)高、下(南)低的情況。反距離加權(quán)插值得到的結(jié)果為,在高點和低點周圍出現(xiàn)明顯的凹包和凸包,“牛眼”現(xiàn)象較為明顯,基本符合反距離加權(quán)插值的特點?？死锝?、規(guī)則樣條、張力樣條3種插值得到的結(jié)果大致相同,但克里金插結(jié)果較平滑,在高低點處的過渡較平緩,優(yōu)于規(guī)則樣條和張力樣條插值;規(guī)則樣條和張力樣條插值結(jié)果,在Z1和Z4地層上部出現(xiàn)輕微的“牛眼”現(xiàn)象。綜上分析可知,在符合空間自相關(guān)時,克里金插值優(yōu)于其他3種插值方法。

由表2可知,Z0、Z2、Z3地層的地形起伏較大,其反距離加權(quán)插值都出現(xiàn)較明顯的“牛眼”現(xiàn)象,特別是較高和較低區(qū)域的“牛眼”更加突顯;規(guī)則樣條和張力樣條插值法在Z0和Z3地層的插值結(jié)果都較平滑,在Z2層兩種方法都出現(xiàn)輕微的“牛眼”,這是由于該層局部起伏較大。

對插值結(jié)果進行對比后,對鉆孔數(shù)據(jù)采用留一法驗證各插值方法的精度,結(jié)果如圖8所示,其中數(shù)據(jù)點到紅色實線的距離表示各插值方法的誤差值。

圖8 留一法驗證結(jié)果

由Z1和Z4地層的折線圖可以看出,克里金插值結(jié)果最貼合紅線,即最接近實測值;反距離加權(quán)插值結(jié)果誤差最大;規(guī)則樣條和張力樣條結(jié)果的誤差介于兩者之間。結(jié)果表明,當?shù)貙臃峡臻g自相關(guān)時,克里金插值方法優(yōu)于其他插值方法。由Z0、Z2、Z3地層的折線圖可以看出,由于3個地層的起伏較大,所有插值方法都與實測值有較大的誤差。Z0和Z3中張力樣條、規(guī)則樣條、克里金插值折線大致呈相同趨勢,幾乎貼近于實測值,反距離加權(quán)法誤差最為明顯;Z2中張力樣條、規(guī)則樣條、克里金插值整體上優(yōu)于反距離加權(quán)插值,但在個別點誤差大于反距離加權(quán)插值。

對不同插值方法的插值精度用平均誤差(ME)和均方根誤差(RMSE)進行評價,結(jié)果見表3。張力樣條插值法在Z0地層ME和RMSE最小,精度最好;Z1地層克里金插值法的ME和RMSE最小;Z2選擇規(guī)則樣條插值法的ME和RMSE最小;Z3和Z4交叉驗證結(jié)果較特殊,Z3規(guī)則樣條插值法的RMSE最小,張力樣條插值法的ME最小,Z4克里金插值法的RMSE最小,規(guī)則樣條插值法的ME最小。

表3 交叉驗證結(jié)果 m

為了選擇適用于各地層最優(yōu)的插值方法,根據(jù)插值結(jié)果、交叉驗證與精度評價可知:Z0選擇張力樣條插值法,Z1和Z4選擇克里金插值法,Z2和Z3選擇規(guī)則樣條插值法。然后根據(jù)規(guī)則格網(wǎng)生成點模型,以便后邊生成地質(zhì)體模型。

3.3.2 地質(zhì)體模型生成

將插值加密得到的虛擬鉆孔點數(shù)據(jù)整理成Excel數(shù)據(jù)并導(dǎo)入Dynamo中,通過Topography.ByPoints節(jié)點生成虛擬鉆孔點模型。生成的點模型如圖9所示。

圖9 地層點模型

將生成的地層點模型通過Topodraphy.PolySurface節(jié)點生成地層模型,如圖10所示。通過Surface.PerimeterCurves節(jié)點生成地形曲面的邊界,并將邊界向上向下投影到一定的平面。

圖10 地層面模型

將生成的地層面模型及兩個投影面進行Loft放樣,生成體模型,通過Geometry.SplitByTools節(jié)點與地形曲面進行裁剪,得到地質(zhì)體模型,如圖11所示。通過Springs.FamilyInstance.ByGeometry節(jié)點將生成的地質(zhì)體導(dǎo)入Revit軟件,可在Revit軟件中更改各地層的屬性,如圖12所示。

圖11 地質(zhì)體模型

圖12 Revit中三維地質(zhì)模型

3.4 礦山儲量估算系統(tǒng)實現(xiàn)與應(yīng)用

將三維地質(zhì)模型導(dǎo)入SuperMap中,與DOM和DEM生成的三維地形模型進行融合,存儲至數(shù)據(jù)庫中,基于C#進行二次開發(fā),生成礦山儲量估算系統(tǒng)。該系統(tǒng)從2022年初在礦區(qū)應(yīng)用以來,滿足了生產(chǎn)需求,提高了工作效率。全區(qū)儲量報表如圖13所示。

圖13 全區(qū)儲量報表

4 結(jié) 語

為滿足露天礦山生產(chǎn)管理的實際需求,本文提出了一種基于混合插值算法的露天礦地質(zhì)體三維建模方法,并據(jù)此進行礦山儲量估算,為礦山企業(yè)實現(xiàn)生產(chǎn)科學(xué)管理提供了數(shù)據(jù)支撐。得出結(jié)論如下:①將鉆孔數(shù)據(jù)分層,結(jié)合空間自相關(guān)分析,選出每個地層最適合的插值方法,不局限于一個區(qū)域的插值只選取單一的插值方法,提高了鉆孔插值的精度,進而提高了地質(zhì)體三維建模的精度;②通過Revit與Dynamo,采用“點—面—體”的建模思路構(gòu)建三維地質(zhì)體模型,基于Revit平臺構(gòu)建的模型可以與其他BIM及GIS軟件進行交互。雖然本文方法及研究成果在露天礦山得到應(yīng)用,但對于地質(zhì)結(jié)構(gòu)復(fù)雜情況的露天礦產(chǎn)儲量估算仍需進一步研究。