談樹成，虎雄崗，金艷珠，蔣順德

（云南大學(xué)地質(zhì)研究所，云南 昆明650091）

固體礦產(chǎn)資源儲量估算是礦產(chǎn)勘查的基本任務(wù)之一，也是礦山生產(chǎn)管理的基礎(chǔ)和礦產(chǎn)地質(zhì)勘查工作成果的總結(jié)，通過對固體礦產(chǎn)資源儲量進(jìn)行估算，可以綜合性地評價礦床的工業(yè)意義，恰當(dāng)?shù)卮_定礦山企業(yè)投資和生產(chǎn)規(guī)模，合理地選擇礦床開采程序和開采方法[1]。常用的固體礦產(chǎn)資源儲量估算方法主要有包括以幾何計算為基礎(chǔ)的傳統(tǒng)方法和以統(tǒng)計學(xué)為基礎(chǔ)的數(shù)學(xué)統(tǒng)計方法兩大類。前者包括算數(shù)平均法、地質(zhì)塊斷法、最近地區(qū)法、等值線法和斷面法等；而后者主要有地質(zhì)統(tǒng)計法、SD儲量估算法和距離平方反比法等[2－3]。地統(tǒng)計學(xué)方法由于考慮了區(qū)域化變量的結(jié)構(gòu)性及隨機(jī)性這二重性特征，在估值精度上與其它方法相比具有明顯的優(yōu)越性，因而獲得了日益廣泛的應(yīng)用，已成為固體礦產(chǎn)資源儲量估算的標(biāo)準(zhǔn)手段之一[4]。國內(nèi)外有關(guān)礦產(chǎn)資源儲量的研究成果較多，在美國，提供地質(zhì)統(tǒng)計學(xué)儲量估算報告已成為礦山企業(yè)獲得銀行貸款的必備條件之一[5]，國內(nèi)學(xué)者陳伯茂以特定礦區(qū)作為研究對象，運(yùn)用地統(tǒng)計學(xué)計算方法進(jìn)行礦產(chǎn)資源儲量分級[6]；孫進(jìn)、張佳榮、侯斌等學(xué)者則詳細(xì)地對比分析了地統(tǒng)計學(xué)與經(jīng)典統(tǒng)計學(xué)在礦產(chǎn)儲量計算中的特點(diǎn)，認(rèn)為應(yīng)用地統(tǒng)計學(xué)方法計算礦產(chǎn)資源儲量可克服傳統(tǒng)計算方法之不足[7]。這些已有的研究成果可為本次研究提供一定的借鑒意義。

然而，由于礦產(chǎn)資源儲量估算涉及到的數(shù)據(jù)量很大，單獨(dú)利用地統(tǒng)計學(xué)方法對固體礦產(chǎn)資源儲量估算的空間數(shù)據(jù)進(jìn)行處理顯得力不從心。近年來，隨著GIS技術(shù)的飛速發(fā)展，使得GIS成為管理、處理分析和顯示輸出空間數(shù)據(jù)之強(qiáng)大工具，其高效的實(shí)時顯示和輸出功能使它廣泛地被運(yùn)用到礦政管理行業(yè)當(dāng)中。因此，運(yùn)用地統(tǒng)計學(xué)與GIS相結(jié)合的方法，探討固體礦產(chǎn)資源儲量的估算，可為國土資源部門對礦產(chǎn)資源儲量的管理提供技術(shù)支撐。

1 研究區(qū)域概況與數(shù)據(jù)源

云南省元江鎳礦礦區(qū)位于唐古拉－昌都－蘭坪－思茅褶皺西，區(qū)域出露地層主要為哀牢山變質(zhì)巖系和三疊系，自酸性至超基性巖均有出露。區(qū)域巖漿活動強(qiáng)烈，礦化帶主要呈北西－南東向展布，整個礦化帶長約2km，礦化較為連續(xù)，工程鎳礦品位1.00%～1.43%，平均品位1.14%。區(qū)內(nèi)礦體呈面型分布，長軸方向長約760～1280m，北部寬760～1280m，南部寬40～200m，礦體規(guī)模為中型，形態(tài)較簡單，厚度一般為1.00～10.80m，最厚為25.60m，平均厚為5.05m。

本次研究選取礦區(qū)162個工程鉆孔數(shù)據(jù)，均包含礦體厚度（m）和礦體平均品位（%）。根據(jù)鎳礦礦產(chǎn)資源儲量核實(shí)的規(guī)定，將Ni≥1.00%作為最低工業(yè)品位，以厚度≥1m作為最低工業(yè)開采厚度，因此見礦工程中有見工業(yè)礦工程41個，見低品位礦工程114個。根據(jù)工程布置原則和質(zhì)量，圈定了332類和333類資源儲量級別范圍。

2 空間數(shù)據(jù)的處理

2.1 正態(tài)分布檢驗(yàn)

利用直方圖和正態(tài)QQ圖分析功能，可對礦產(chǎn)資源儲量空間樣本數(shù)據(jù)的分布情況進(jìn)行檢驗(yàn)，并可選擇恰當(dāng)?shù)淖儞Q模型對樣本數(shù)據(jù)進(jìn)行計算，使得樣本數(shù)據(jù)在經(jīng)過函數(shù)變換后有更好的正態(tài)分布效果，達(dá)到地統(tǒng)計插值的要求。如圖1所示為沒有經(jīng)過變換、經(jīng)過對數(shù)變換和冪變換后的工業(yè)礦礦體的品位屬性直方圖，通過對比可以發(fā)現(xiàn)，經(jīng)過對數(shù)變換后的樣本數(shù)據(jù)更好的接近正態(tài)分布。

圖1 樣本數(shù)據(jù)直方圖

2.2 趨勢分析

如圖2所示，可以看出礦產(chǎn)資源儲量樣本數(shù)據(jù)在南北方向（Y軸）上表現(xiàn)出一定的倒“U”型，即工業(yè)礦礦體的品位從南向北先變大后變小，而在東西方向（X軸）上呈現(xiàn)出上升的趨勢。因此，在地統(tǒng)計插值時可對其進(jìn)行二階多項(xiàng)式擬合，達(dá)到較好的全局趨勢效果。

3 研究方法概述

3.1 區(qū)域化變量與變異函數(shù)

任何空間上分布的變量均稱為區(qū)域化變量，區(qū)域化變量是以空間點(diǎn)X的三個直角坐標(biāo)（x，y，z）作為自變量的隨機(jī)場Z（X），即Z（X）＝Z（x，y，z）。當(dāng)對其進(jìn)行一次觀測后就得到它的一個現(xiàn)值Z（X），它是一個普通的三元值函數(shù)或是空間點(diǎn)函數(shù)。

變異函數(shù)是一個距離函數(shù)，它描述不同位置變量的相似性，其值越大，相似性越差。通常，隨機(jī)變異函數(shù)值隨距離矢量h的增大而增大，直到h達(dá)到一定值時，變異函數(shù)達(dá)到最大值，然后保持這個常數(shù)不再發(fā)生變化。

3.2 克里格法

通過對研究區(qū)數(shù)據(jù)的分析后，得出樣本數(shù)據(jù)服從一定的正態(tài)分布規(guī)律。因此本次研究擬采用普通克里格、簡單克里格和泛克里格方法對礦體的品位和厚度進(jìn)行空間插值，考慮到礦體品位在空間上存在一定的趨勢，分別利用每種插值方法的二階多項(xiàng)式變換進(jìn)行模擬，從而計算出礦產(chǎn)資源儲量。

圖2 樣本數(shù)據(jù)趨勢分析圖

3.2.1 普通克里格法

當(dāng)區(qū)域化變量滿足二階平穩(wěn)假設(shè)時，通過變異函數(shù)所揭示的空間變異性，估計鄰域內(nèi)變量值來估算某一待估域值，而任一待估塊段V的真值Z（xi）的估算值Z＊v是估計鄰域內(nèi)n個信息值Z（xi）的線性組合，它假設(shè)數(shù)據(jù)變化成正態(tài)分布，認(rèn)為區(qū)域化變量Z的期望值是未知的。插值過程類似于加權(quán)滑動平均，權(quán)重值的確定來自于空間數(shù)據(jù)分析，其計算公式為：

式中：為估算點(diǎn)的值；λi為權(quán)重系數(shù)，是各已知樣品在Z（xi）估算影響大小的系數(shù)，而估算Z＊V的可靠性主要取決于怎樣計算或選擇權(quán)重系數(shù)λi。xi表示估算點(diǎn)的位置。

普通克里格插值是滿足假設(shè)的區(qū)域化變量的線性估計，它假設(shè)數(shù)據(jù)的變化呈正態(tài)分布，假設(shè)區(qū)域化變量Z的期望值是未知的，如果為已知的常數(shù)，則為簡單克里格方法。當(dāng)樣本數(shù)據(jù)在空間上存在明顯的趨勢，普通克里格方法將不再適用，此時需要考慮泛克里格插值方法。

3.2.2 簡單克里格法

簡單克里格是區(qū)域化變量的線性估計，它假設(shè)數(shù)據(jù)變化成正態(tài)分布，認(rèn)為區(qū)域化變量Z的期望值為已知的某一常數(shù)。與普通克里格法不同的是，公式Z（s）＝μ＋ε（s）中的u是一個已知常量，而并非一個未知常量，如圖3所示。因?yàn)榧僭O(shè)確切已知常量μ，那么也確切已知數(shù)據(jù)位置上的ε（s），而對于普通克里金法，如果估算了μ，也會估算ε（s）。因此，在ε（s）已知的情況下，可以比ε（s）未知時更好地估算自相關(guān)。

圖3 Z（s）＝μ＋ε（s）模型示意圖

3.2.3 泛里格法

泛克里格假設(shè)數(shù)據(jù)中存在主導(dǎo)趨勢，且該趨勢可以用一個確定的函數(shù)或多項(xiàng)式來擬合。在進(jìn)行泛克里格分析時，首先分析數(shù)據(jù)中存在的變化趨勢，獲得擬合模型；其次，對殘差數(shù)據(jù)進(jìn)行克里格分析；最后，將趨勢面分析和殘差分析的克里格結(jié)果加和，得到最終結(jié)果。由此可見，克里格方法明顯優(yōu)于趨勢面分析，泛克里格的結(jié)果也要優(yōu)于普通克里格的結(jié)果。

3.3 留一交叉驗(yàn)證法

留 一 交 叉 驗(yàn) 證 法 （Leave－One－Out Cross Validation，簡記為LOO－CV），是交叉驗(yàn)證的重要方法之一。該方法假設(shè)有N個采樣點(diǎn)，對于第i個訓(xùn)練樣本，將其取出，對剩下N－1個樣本進(jìn)行訓(xùn)練，得到?jīng)Q策函數(shù)，并用其測試第i個訓(xùn)練樣本，該過程重復(fù)N次，并得到N個模型，用這N個模型最終的驗(yàn)證集的分類準(zhǔn)確率的平均數(shù)作為此下LOOCV插值模型的精度指標(biāo)，用此方法求出的誤差對于實(shí)際中的測試誤差來說幾乎是無偏的。

4 插值結(jié)果與檢驗(yàn)分析

4.1 插值結(jié)果

為了使工程樣本數(shù)據(jù)空間插值達(dá)到最好的效果，本次研究選取多種插值方法并對其插值結(jié)構(gòu)進(jìn)行比較，以期選取出高效、最佳的插值模型。主要有：普通克里格無變換插值、普通克里格二階多項(xiàng)式插值、簡單克里格無變換插值、簡單克里格二階多項(xiàng)式插值、泛克里格無變換插值和泛克里格二階多項(xiàng)式插值6種模型。

4.2 交叉驗(yàn)證

為了充分比較和評價以上6種插值方法的結(jié)果。利用AML語言進(jìn)行編程，實(shí)現(xiàn)留一交叉驗(yàn)證的計算機(jī)輔助自動化處理?；舅悸肥牵好看闻R時從某時刻降雨量觀測值中移出一個觀測值，然后基于余下的數(shù)據(jù)采用各種插值算子來重新估計它的值，通過計算和比較采樣點(diǎn)預(yù)測值與真實(shí)值之間的誤差大小來評價插值模型的精度。其評價原則為：接近0的平均誤差、較小的均方根預(yù)測誤差、與均方根預(yù)測誤差相似的平均標(biāo)準(zhǔn)誤差、接近0的標(biāo)準(zhǔn)平均值預(yù)測誤差和較大的預(yù)測值與真實(shí)值之間的相關(guān)系數(shù)。表1所示為6種插值方法的結(jié)果誤差比較。

表1 幾種克里格插值后的參數(shù)對比

4.2.1 平均誤差分析

從表1中可得到：普通克里格無變換插值、普通里格二階多項(xiàng)式插值和泛克里格無變換插值誤差較高，誤差絕對值均大于0.006。而簡單克里格無變換插值和簡單克里格二階多項(xiàng)式插值所得到的誤差均較小，其中以簡單克里格二階多項(xiàng)式插值平均值誤差最小，僅為－0.0004，因此，就平均誤差而言，簡單克里格二階多項(xiàng)式插值效果無疑是最好的。

4.2.2 均方根誤差分析

從表1中可得到：泛克里格二階多項(xiàng)式插值的均方根誤差最大，達(dá)到0.8919。其它方法的均方根誤差絕對值均小于0.5，其中以簡單克里格二階多項(xiàng)式插值的均方根誤差最小，為0.4214，說明在Z（s）＝μ＋ε（s）中常量u已知的情況下進(jìn)行礦體品位的二階模擬可以得到很好的反映出其真實(shí)情況，得到很好的效果，因此，從均方根誤差這一指標(biāo)來看，簡單克里格二階多項(xiàng)式插值為效果最好，精度最高的一種方法。

4.2.3 其它評價指標(biāo)分析

除了利用平均誤差和均方根誤差進(jìn)行檢驗(yàn)外，本次研究還利用了標(biāo)準(zhǔn)平均值誤差、標(biāo)準(zhǔn)均方根誤差和平均標(biāo)準(zhǔn)誤差等其余3項(xiàng)指標(biāo)進(jìn)行對比分析?？梢园l(fā)現(xiàn)，除了利用普通克里格二階多項(xiàng)式插值得到最小的標(biāo)準(zhǔn)平均值誤差外，剩下的兩項(xiàng)指標(biāo)均以簡單克里格二階多項(xiàng)式插值誤差最小。從而得到簡單克里格二階多項(xiàng)式插值為本次研究中估算工程品位精度最高、效果最好的空間插值方法。

運(yùn)用同樣的方法對工業(yè)礦礦體厚度進(jìn)行地統(tǒng)計插值，不同的是，工業(yè)礦礦體厚度在空間分布上沒有表現(xiàn)出一定的空間趨勢，因此，在地統(tǒng)計插值前不需要對其進(jìn)行變換處理。本次研究選取普通克里格、簡單克里格和泛克里格3種插值模型對工業(yè)礦礦體厚度進(jìn)行空間插值，通過比較后選擇出最優(yōu)模型。

經(jīng)過簡單克里格二階多項(xiàng)式變換插值后的工業(yè)礦礦體品位預(yù)測圖和普通克里格無變換插值后工業(yè)礦礦體厚度預(yù)測圖如圖4所示。

5 礦產(chǎn)資源的估算及其結(jié)果分析

礦產(chǎn)資源儲量（礦石量和金屬量）估算公式為：

式中：Q為塊段礦石量，t；S為塊段面積，m2；H為塊段平均厚度，m；D為塊段平均體重，t／m3；P為塊段金屬量，t；C為塊段平均品位，%。

根據(jù)上式中的計算方法，將礦區(qū)工業(yè)礦礦體品位及厚度數(shù)據(jù)柵格化后，再進(jìn)行空間數(shù)據(jù)的疊加，可分別得到332類和333類工業(yè)礦鎳礦礦石量和鎳礦金屬量，再利用礦區(qū)鎳礦資源儲量核實(shí)報告中的332類和333類工業(yè)礦鎳礦礦石量和鎳礦金屬量對本次研究進(jìn)行檢驗(yàn)，從而確定研究方法和技術(shù)的可靠性，其中：

變幅（%）＝（本次研究所估算的量－資源儲量核實(shí)結(jié)果）／資源儲量核實(shí)結(jié)果×%

結(jié)果對比如表2所示。通過表2中的結(jié)果對比可以得出，本次研究結(jié)果與資源儲量核實(shí)結(jié)果基本一致，變幅不到5%。其中，通過本次研究得到的332類鎳礦礦石量與332類鎳礦金屬量分別為578210t和6920t，而332類資源儲量核實(shí)結(jié)果中礦石量為572100t，金屬量為6835t，變幅分別為1.07%和1.24%；通過本次研究得到的333類鎳礦礦石量與333類鎳礦金屬量分別為128678t和1469t，而333類資源儲量核實(shí)結(jié)果中礦石量為123200t，金屬量為1407t，變幅分別為4.45%和4.41%。這主要是由于332類和333類鎳礦資源工程布置質(zhì)量的差異。根據(jù)國土資源部發(fā)布的《銅、鉛、鋅、銀、鎳、鉬礦地質(zhì)勘查規(guī)范》，此鎳礦資源儲量勘查中332類資源量的網(wǎng)格間距確定為100m×100m，而333類資源量的網(wǎng)格間距確定為200m×200m，從而造成了333類資源量在本次研究與資源儲量核實(shí)結(jié)果中存在較大的差異。

圖4 工業(yè)礦品位和厚度克里格插值結(jié)果

表2 研究結(jié)果與鎳礦礦產(chǎn)資源儲量核實(shí)對比分析表

6 結(jié)論

通過對工程鉆孔數(shù)據(jù)的空間分析和變換處理，利用地統(tǒng)計學(xué)方法對礦體的品位與厚度進(jìn)行檢驗(yàn)分析，尋求礦體的空間分布規(guī)律，結(jié)合GIS空間數(shù)據(jù)的管理與分析功能，估算出研究區(qū)域工業(yè)礦332類資源儲量礦石量和金屬量以及333類資源儲量礦石量和金屬量。結(jié)果顯示，利用地統(tǒng)計學(xué)與GIS技術(shù)相結(jié)合的方法對固體礦產(chǎn)資源儲量進(jìn)行估算，其結(jié)果與固體礦產(chǎn)資源儲量核實(shí)結(jié)果基本一致，332類資源儲量變幅不超過1.5%，而333類資源儲量變幅在4%～5%之間，充分論證了研究結(jié)果的準(zhǔn)確性。研究表明，方法更為恰當(dāng)、技術(shù)更為嚴(yán)密、結(jié)果更為可靠等特點(diǎn)，充分顯示了基于地統(tǒng)計學(xué)與GIS對固體礦產(chǎn)資源儲量估算進(jìn)行研究的優(yōu)越性。因此，研究方法可為國土資源部門對礦產(chǎn)資源儲量的管理提供技術(shù)支撐，成果可為其他礦區(qū)礦產(chǎn)資源儲量的估算提供重要的參考價值。

[1]陳國旭.支持傳統(tǒng)方法的礦產(chǎn)資源儲量估算軟件研究進(jìn)展[J].國土資源信息化，2009（2）：8－12.

[2]袁勇，王志龍.儲量計算克立格法與SD法的對比分析[J].內(nèi)蒙古石油化工，2009（4）：42－43.

[3]王正帥，顧和和.數(shù)字化礦山資源儲量計算方法[J].煤田地質(zhì)與勘探，2010，38（3）：6－9.

[4]萬昌林，金銅標(biāo).地質(zhì)統(tǒng)計學(xué)在礦產(chǎn)資源信息研究中的應(yīng)用[J].有色冶金設(shè)計與研究，2002（4）：18－21.

[5]馬洪濱，熊俊楠.基于地質(zhì)統(tǒng)計學(xué)的儲量估算系統(tǒng)[J].煤炭學(xué)報，2007，32（3）：267－271.

[6]陳伯茂.地質(zhì)統(tǒng)計學(xué)在礦產(chǎn)儲量分級中的應(yīng)用[J].地質(zhì)與勘探，1986，22（2）：24－30.

[7]孫進(jìn)，張佳榮，侯斌.礦產(chǎn)儲量計算經(jīng)典統(tǒng)計與地質(zhì)統(tǒng)計學(xué)方法的對比分析[J].采礦技術(shù)，2005，5（2）：80－81.