基于鉆孔數據進行礦體克里金插值分析

劉亞靜 ，趙秀梅 ，2，甘德清 ，2

（1.華北理工大學 礦業(yè)工程學院，河北 唐山 063210；2.河北省礦業(yè)開發(fā)與安全工程實驗室，河北 唐山 063210）

0 引言

我國鎢資源儲量和生產量都居于世界第一，隨著對鎢需求量的增加，大規(guī)模的礦山開采給生態(tài)環(huán)境帶來巨大的影響[1-2]。為了保障礦區(qū)的可持續(xù)發(fā)展，準確模擬礦區(qū)地下巖層的分布形態(tài)[3-4]，李海港等[5]學者對鎢礦進行三維建?？梢暬Ｓ捎陔x散鉆探數據建立的曲面不能準確反映巖層的分布，需要確定精確的空間插值方法，利用最優(yōu)插值方法，對礦區(qū)金屬元素進行儲量計算。

研究基于某礦區(qū)鉆孔數據，使用GIS空間地質統計學和相關軟件，分析數據特征[6-8]，通過交叉驗證方法，得到最優(yōu)插值方法，為礦體儲量計算提供了依據。以期為該礦山地理信息領域的進一步研究提供參考。

1 研究區(qū)數據基礎

研究礦井樣本點經整理共110個樣本點數據，將結果整理為EXCEL表格數據。為了將EXCEL數據轉換為屬性特征表，使用ArcGIS的添加x，y數據工具進行轉換，利用ArcGIS定義投影工具，建立相應的數據系統[9-10]，得到鉆孔的空間位置分布，如圖1所示。

圖1 礦區(qū)點分布Fig.1 Distribution of the mining area

2 礦區(qū)鉆孔數據空間特征分析

2.1 地質統計學基礎

克里金插值法又稱空間插值，是一種基于變量理論和結構分析，對有限區(qū)域內的區(qū)域變量進行無偏最優(yōu)估計的方法。采用已知樣本點之間的距離和空間方位來反映空間相關性，對未知點進行線性無偏最優(yōu)預計[11]。

克里金法插值公式如式（1）所示。

其中：λi為第i個樣本點的未知權重；N是已知點的數目；Z（xi）為位置i的樣本值；Z（x0）為未知樣本值[13]。

2.2 空間數據分析

ArcGIS軟件的地質統計學模塊提供了許多功能，其中直方圖法、正態(tài)QQ圖法，趨勢分析和半方差/協方差函數云可分析數據[12]。當獲得樣本點時，首先需要對數據進行分析，分析數據是否服從正態(tài)分布，是否存在趨勢，檢驗數據是否合乎實際情況，剔除明顯差異點。

2.2.1 礦井數據分布直方圖和正態(tài)Q Q圖檢驗

采用直方圖法和正態(tài)QQ圖法對110個樣本點進行了分析，比較樣本點數據對數變化前后的正態(tài)分布。由圖2、圖3可得，樣本點數據經過對數變換后更符合正態(tài)分布。

圖2 直方圖Log變化后分布Fig.2 Histogram Log distribution after the change

圖3 正態(tài)QQ圖Log變化后分布Fig.3 Log distribution after the change in the normal QQ chart

2.2.2 趨勢分析

感興趣區(qū)域的表面在各處出現漸變時，將該表面與采樣點擬合，得到漸進趨勢的平滑表面[11]。使用ArcGIS軟件地質統計學模塊中的趨勢分析工具對數據進行處理，結果如圖4所示，數據在X和Y方向上都呈二階分布。更加直觀地了解整個數據在區(qū)域范圍內有較強的趨勢變化。

圖4 礦區(qū)地質云數據趨勢分布圖Fig.4 Distribution of geological cloud data in the mining area

2.2.3 空間自相關

利用半變異函數的協方差云圖對所得數據進行分析。其中，水平坐標反映的是兩點之間的距離，而垂直坐標反映的是-1～+1之間的半變異值，區(qū)域化變量距離越近，相似性越大。結果如圖5顯示，半變異值主要集中在-0.45～+1之間，空間自相關較強。因此，可以得出該數據可用于克里金插值的結論。

圖5 半變異函數云圖Fig.5 Semi-variogram cloud

3 比較不同空間插值選取最優(yōu)插值方式

3.1 有效性評價

交叉驗證是指從樣本數據集中刪除一點，使用其他樣本值估計此點屬性，并計算與實際值之間誤差，重復以上步驟直到遍歷所有采樣點。在多種插值方法，多個參數之間選取最優(yōu)方法，采用交叉驗證來比較優(yōu)劣。

交叉驗證的參數包含多個誤差，在本文采用了標準均值（SEM）、標準化均方根誤差（RMSSE）、平均標準誤差（ASE）、均方根誤差（RMSE）和平均誤差（ME）進行評估[14-15]。平均標準誤差越接近于零，準確度就越高[14]。同時，RMSSE誤差和SEM誤差越小，精度越高。

3.2 分布趨勢剔除

采用最常用的普通克里金插值法和簡單克里金插值法，比較了常數、一階和二階趨勢面與數據符合情況[15-17]。其中變量參數統一設置為高斯變異函數模型，步長數為12，模型類型為穩(wěn)定，計算偏基臺值為true，搜索鄰域類型為標準工具，表1為誤差結果。觀察表1發(fā)現，普通克里金插值和簡單克里金插值中，一階趨勢面ME、SEM、RMSSE誤差較小，效果較好。因此，本研究基于一階趨勢面進行以下數據分析。

表1 趨勢面階數誤差結果Tab.1 Order error results of the trend surface

3.3 克里金不同插值方法插值比較

3.3.1 普通克里金法插值分析

普通克里金法是一種線性克里金法，區(qū)域變量滿足二階平穩(wěn)且期望為未知常數，較符合實際情況，應用范圍較廣。本研究的初始步長為6，間隔為3，最終步長為21?；谝浑A趨勢面，對樣本數據進行對數變換，選擇高斯、球面和指數模型。結果如表2，球面函數模型和指數函數模型的均方根誤差和平均標準誤差較大；高斯函數模型的平均標準誤差最接近均方根誤差，均方根誤差最低。相比之下，高斯模型的均方根誤差和平均標準誤差都很小，因此選擇高斯模型作為最優(yōu)參數。

表2 普通克里金誤差統計Tab.2 Statistics of the Kriging error

對高斯模型的不同步長數進行對比來確定最優(yōu)步長數，要求各個誤差值越小越好，結果如圖6曲線所示。從圖6對比結果可知，誤差值都較低時，步長數為6，故確定步長數6為最優(yōu)參數。

圖6 普通克里金法插值誤差曲線Fig.6 Interpolation error curve of ordinary Kriging method

3.3.2 簡單克里金法插值分析

簡單克里金法，區(qū)域化變量滿足二階平穩(wěn)，期望值通常很難準確估計，估計值精度低。逐項比較誤差的精度評定是一個復雜的過程，且含有不確定因素，因此定義4個評價指數S，V，W，Z，計算公式為：

其中：n 代表步長數，S，V，W，Z 四個值越小代表插值精度越高。根據以上公式，簡單克里金誤差評價指數見表3。

表3 簡單克里金誤差評價指數Tab.3 Error evaluation indices of simple Kriging method

綜合比較，高斯球面函數模型，V和W評價指數相對較小，分別為3.979和3.187，精度更高，固選取高斯函數為最優(yōu)模型。分析圖7，當步長數為6，簡單克里金法插值的誤差相對較小，故確定步長數6為最優(yōu)參數。

圖7 簡單克里金法插值誤差曲線Fig.7 Interpolation error curve of simple Kriging method

3.3.3 兩種方法插值結果比較

分析插值結果得到，當函數模型選取高斯，步長數為6，采用普通克里金和簡單克里金兩種方法進行插值結果分析，結果如表4。

表4 高斯模型步長數為6時誤差統計Tab.4 Error statistics of step number 6 in the Gaussian model

根據誤差統計表分析，普通克里金法的ME、RMSE、SEM和ASE誤差值都為最低，普通克里金方法插值效果最好。因此，文章選取高斯球面函數，步長數為6，普通克里金方法的空間插值方法最好。

4 三維可視化圖

基于此數據，針對某礦體通過Surpac軟件建立三維地質模型[18-20]，有以下幾個步驟：

（1）根據勘探線剖面圖繪制各條勘探線的刨面；（2）將線狀剖面圖依次連接，形成三維礦體模型；（3）根據模型的坐標范圍，建立三維模型的塊體模型，建成塊體模型如圖8所示；（4）通過三角網對塊體模型進行約束，并保持礦體的內部局限性，形成礦體約束模型；（5）對約束后的塊體模型進行賦值；（6）賦值后的模型進行克里金插值，將不同品位區(qū)間賦予不同顏色。

圖8 三維地質模型Fig.8 Three-dimensional geological model

在克里金插值中，考慮了信息樣本的形狀和大小、信息樣本與待估塊的空間分布位置以及品位的空間結構等主要幾何特征。為了實現無偏、線性和精確的估計，對每個品位值進行賦值，然后用加權平均法估計區(qū)塊品位。

將部分礦體模型進行克里金插值，不同品位值區(qū)間賦予不同顏色如表5，結果如圖9所示。

表5 品位值區(qū)間與顏色對應表Tab.5 Grade value range and color coreespondence table

圖9 礦體某金屬元素品位分布Fig.9 Grade distribution of a metal element in ore body

從圖9可以看到，礦體某金屬元素品位值顏色分布情況，首先淡藍色區(qū)域最大，接著是深紫色區(qū)域，然后是淡紫色區(qū)域。根據表5得到，該礦體某金屬元素品位值范圍在30%～60%之間，其中品位值在40%～50%的區(qū)域范圍最大，為采礦提供了參考。

5 結論

研究基于實際鉆孔數據，采用相關軟件對數據進行了分析，得到最優(yōu)插值方法，探索得到結論。

（1）對110個鉆孔數據進行空間插值分析，從空間自相關的半變異函數云圖中得到，數據集中在-0.45～+1的范圍內，局部變量的間隔越近，相似度越大，說明數據可以用于克里金插值。

（2）從插值結果來看，當采用普通克里金，高斯函數模型，步長設置為6時，更能準確地反映地層的空間分布。

（3）通過建模，對模型進行克里金插值，將插值后的模型進行不同屬性著色，更加直觀地看出礦體元素品位值的分布。