一種混合插值算法的露天礦產儲量估算方法及其應用

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

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

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

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

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

1 研究方法

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

圖1 方法流程

1.1 插值算法

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

1.1.1 反距離加權法(IDW)

反距離加權法的表達式[15]為

(1)

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

1.1.2 徑向基函數法

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

(2)

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

1.1.3 Kriging法

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

(3)

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

(4)

(5)

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

(6)

(7)

1.1.4 空間自相關分析

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

圖2 半變異函數結構

1.1.5 交叉驗證與精度評價

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

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

(8)

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

(9)

由于自然界的地質體形狀各異,任何插值方法都無法適應所有的地質體,因此需要對研究區進行具體分析,從而選出最合適的插值方法。本文首先對鉆孔分層數據進行自相關分析,然后選擇相應的插值方法。

1.2 地質體三維建模

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

1.2.1 Revit和Dynamo可視化編程

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

1.2.2 建模流程

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

圖3 建模流程

2 礦山儲量估算系統設計

礦產儲量估算系統采用C/S的結構模式,在Visual Studio 2013 環境下,運用C#語言,在SuperMap平臺上開發,將礦山地質體三維模型與礦產儲量估算系統進行融合。系統包含數據管理、礦山二三維顯示、儲量估算3個模塊。數據管理模塊包括礦區的DOM、DEM及鉆孔數據的管理;礦山二三維顯示模塊包括礦體地質體三維模型、三維地形模型、矢量模型的顯示,以及對模型的交互瀏覽、平移、縮放、旋轉等操作。儲量估算模塊包括對開采量的計算及報表輸出等。系統功能結構如圖4所示。

圖4 系統功能結構

3 試驗與分析

3.1 研究區概況

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

圖5 研究區

3.2 空間數據采集與處理

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

圖6 鉆孔分布

3.3 三維地質體構建

3.3.1 鉆孔插值及精度檢驗

鉆孔數據按照巖性分為白云質灰巖、鮞狀白云巖、白云巖、鮞狀白云巖、灰巖共5層,分別用Z0、Z1、Z2、Z3、Z4表示,并導入數據庫中統一管理。計算各地層數據的變異函數,并繪制變異函數云圖如圖7所示,據此分析各層數據的空間自相關性。

圖7 各地層變異函數云圖

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

表1 Z1、Z4插值結果

表2 Z0、Z2、Z3插值結果

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

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

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

圖8 留一法驗證結果

由Z1和Z4地層的折線圖可以看出,克里金插值結果最貼合紅線,即最接近實測值;反距離加權插值結果誤差最大;規則樣條和張力樣條結果的誤差介于兩者之間。結果表明,當地層符合空間自相關時,克里金插值方法優于其他插值方法。由Z0、Z2、Z3地層的折線圖可以看出,由于3個地層的起伏較大,所有插值方法都與實測值有較大的誤差。Z0和Z3中張力樣條、規則樣條、克里金插值折線大致呈相同趨勢,幾乎貼近于實測值,反距離加權法誤差最為明顯;Z2中張力樣條、規則樣條、克里金插值整體上優于反距離加權插值,但在個別點誤差大于反距離加權插值。

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

表3 交叉驗證結果 m

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

3.3.2 地質體模型生成

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

圖9 地層點模型

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

圖10 地層面模型

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

圖11 地質體模型

圖12 Revit中三維地質模型

3.4 礦山儲量估算系統實現與應用

將三維地質模型導入SuperMap中,與DOM和DEM生成的三維地形模型進行融合,存儲至數據庫中,基于C#進行二次開發,生成礦山儲量估算系統。該系統從2022年初在礦區應用以來,滿足了生產需求,提高了工作效率。全區儲量報表如圖13所示。

圖13 全區儲量報表

4 結 語

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