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

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

0 引言

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

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

1 研究區(qū)數(shù)據(jù)基礎(chǔ)

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

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

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

2.1 地質(zhì)統(tǒng)計學(xué)基礎(chǔ)

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

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

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

2.2 空間數(shù)據(jù)分析

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

2.2.1 礦井?dāng)?shù)據(jù)分布直方圖和正態(tài)Q Q圖檢驗

采用直方圖法和正態(tài)QQ圖法對110個樣本點進行了分析，比較樣本點數(shù)據(jù)對數(shù)變化前后的正態(tài)分布。由圖2、圖3可得，樣本點數(shù)據(jù)經(jīng)過對數(shù)變換后更符合正態(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ū)域的表面在各處出現(xiàn)漸變時，將該表面與采樣點擬合，得到漸進趨勢的平滑表面[11]。使用ArcGIS軟件地質(zhì)統(tǒng)計學(xué)模塊中的趨勢分析工具對數(shù)據(jù)進行處理，結(jié)果如圖4所示，數(shù)據(jù)在X和Y方向上都呈二階分布。更加直觀地了解整個數(shù)據(jù)在區(qū)域范圍內(nèi)有較強的趨勢變化。

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

2.2.3 空間自相關(guān)

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

圖5 半變異函數(shù)云圖Fig.5 Semi-variogram cloud

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

3.1 有效性評價

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

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

3.2 分布趨勢剔除

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

表1 趨勢面階數(shù)誤差結(jié)果Tab.1 Order error results of the trend surface

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

3.3.1 普通克里金法插值分析

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

表2 普通克里金誤差統(tǒng)計Tab.2 Statistics of the Kriging error

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

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

3.3.2 簡單克里金法插值分析

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

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

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

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

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

3.3.3 兩種方法插值結(jié)果比較

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

表4 高斯模型步長數(shù)為6時誤差統(tǒng)計Tab.4 Error statistics of step number 6 in the Gaussian model

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

4 三維可視化圖

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

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

圖8 三維地質(zhì)模型Fig.8 Three-dimensional geological model

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

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

表5 品位值區(qū)間與顏色對應(yīng)表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ū)域。根據(jù)表5得到，該礦體某金屬元素品位值范圍在30%～60%之間，其中品位值在40%～50%的區(qū)域范圍最大，為采礦提供了參考。

5 結(jié)論

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

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

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

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