距離冪次反比法(IDW)資源儲量估算流程改進及應用

陳偉康，高幫飛，張書琛，姚 剛，孟 偉

(1.中鐵資源集團有限公司，北京 100039；2.MKM礦業(yè)股份有限公司，科盧韋齊 20043)

資源儲量估算總在特定的估算域(estimated domain)內(nèi)進行，其核心問題是確定哪些已知樣品點數(shù)據(jù)參與估值以及相應的權重。傳統(tǒng)資源儲量估算，一般依據(jù)規(guī)范指定的工業(yè)指標圈定礦體(品位域)，且僅利用礦體內(nèi)的數(shù)據(jù)進行估值，其主要問題在于：一方面，人為分割了地質(zhì)和礦化的空間分布，使得“圈礦”工作更多地關注品位，而忽略影響礦化分布的地質(zhì)規(guī)律；另一方面，僅用工業(yè)品位域內(nèi)的數(shù)據(jù)進行估值，往往會增強高品位數(shù)據(jù)的連續(xù)性，使得整個礦床的平均品位高估，給后期礦山開發(fā)帶來風險[1]。近20年來，隨著三維礦業(yè)軟件在國內(nèi)的普及和應用，距離冪次反比(inversed distance weighting，IDW)和克里格等估值方法不斷引入，尤其是IDW方法操作相對簡單，對空間數(shù)據(jù)構型(configuration)沒有特殊要求，因而廣泛應用于資源儲量估算[2-3]。前人對IDW估值的改進主要集中在算法優(yōu)化上，如參估樣品數(shù)據(jù)的自動優(yōu)化[4]、品位估計公式中參數(shù)確定的遺傳算法[5]、通過計算權重合集的分維數(shù)來確定冪次指數(shù)[6]等。而對綜合利用現(xiàn)有軟件和技術，通過流程改進來實現(xiàn)估值優(yōu)化卻鮮有報道[7]。本文以剛果(金)Kalumbwe銅鈷礦床為例，基于三維礦業(yè)軟件平臺，探討如何利用分形分析和邊界分析方法來確定資源儲量估算的礦化域，同時引入變異函數(shù)分析輔助確定影響IDW估值的距離和冪次參數(shù)，以期通過IDW估值流程優(yōu)化，提高資源儲量估算的準確度，降低后期礦山開發(fā)風險。

1 IDW估值原理及流程

1.1 IDW估值原理及現(xiàn)行估值流程

IDW是一個局部估值工具。未知塊的插值是以未知塊為中心的三維搜索橢球內(nèi)所有樣品點的加權平均值，已知樣品點的權重為距未知礦塊中心點距離冪次的倒數(shù)，計算見式(1)。

i=1,2,…,n

(1)

式中：Z(x)*為未知塊的估計值；z(xi)為搜索橢球內(nèi)的第i個樣品值；w(xi)為對應樣品權重；di為樣品點到未知塊中心點歐氏距離；p為冪次。IDW估值原理及現(xiàn)行估值流程詳見圖1。

圖1 IDW估值原理及現(xiàn)行估值流程Fig.1 Principles and current processes of the IDW estimation method

1.2 存在問題及改進思路

現(xiàn)行IDW估值流程的主要問題包括：①參與資源儲量估算的樣品集為根據(jù)工業(yè)指標圈定礦體內(nèi)的所有數(shù)據(jù)，這一做法人為割裂了區(qū)域化變量在地質(zhì)和統(tǒng)計上的連續(xù)性；②礦體的品位域在工業(yè)品位和特高值之間，估值時容易增強工業(yè)品位數(shù)據(jù)的連續(xù)性，導致估值結果偏高；③影響IDW估值的距離(與橢球半徑有關)和冪次的確定主要依賴于地質(zhì)工程師的認識水平，缺少量化。

本文針對上述問題，提出一種改進的IDW估值工作流程。具體思路和步驟如下所述。①引入“礦化品位域”的概念，域內(nèi)數(shù)據(jù)要具有地質(zhì)和統(tǒng)計意義上的一致性，進一步通過分形分析給出礦化品位域的范圍，進而采用邊界分析對工業(yè)品位域和礦化品位域的邊界條件進行確認；②抓住用哪些數(shù)據(jù)參與估值及其權重分配的中心問題，主要通過變異函數(shù)的井向變異函數(shù)和全向變異函數(shù)分析提供的塊金效應和最優(yōu)變程信息，對IDW估值的距離和冪次參數(shù)進行約束；③所有的應用均基于現(xiàn)有通用三維礦業(yè)軟件平臺，改進的內(nèi)容通用軟件(三維礦業(yè)軟件+Excel)可以實現(xiàn)，不存在技術難度，且技術的可重復性強。

1.3 礦化品位域的確定

資源儲量的估算域，既可以是以巖性來劃分的地質(zhì)域，也可以是以品位來圈定的品位域。國內(nèi)通行的做法是圈定品位域用于后續(xù)估值流程[8-11]。為了保持品位域內(nèi)樣品數(shù)據(jù)在地質(zhì)和統(tǒng)計上的一致性，應剔除特異值的影響。特異值既包括特高值，也包括特低值[12]。傳統(tǒng)做法是基于規(guī)范推薦或經(jīng)設計院論證后的邊界品位，作為特低品位的閾值，高于該邊界品位的才會被圈入礦體。而特高值確定的方法比較多，如m+3σ法，變異系數(shù)法，直方圖、CDF、概率圖等[13]，統(tǒng)計分析一般基于礦體內(nèi)數(shù)據(jù)。這些方法存在的問題是：統(tǒng)計分析結果受研究尺度影響；數(shù)據(jù)分組等人為因素干擾；主要針對特高值處理，不能給出特低值。采用分形方法可以有效解決這一問題，具體方法原理和計算流程參見前期的研究成果[14]，本文不再贅述。

2 應用實例概況

Kalumbwe礦床位于剛果(金)銅鈷成礦帶西段，為典型的沉積層狀銅鈷礦床[15]。礦體產(chǎn)出受地層控制，主礦體呈EW走向，傾向南，傾角60°～80°；主礦體長約1 km，厚度30～50 m(圖2(a))，礦床由上部氧化礦體和下部硫化礦體組成。勘查階段主要為鉆孔控制，生產(chǎn)階段采用鉆孔+槽探進行聯(lián)合控制。基本采樣間距2 m，樣品總數(shù)4 768件。剛果(金)政府推薦使用的圈礦工業(yè)指標為銅1.0%，鈷0.1%；礦體最小可采厚度和夾石剔除厚度均為2 m。為了便于比較，利用礦化品位也進行了礦體圈連(圖2(b))。

圖2 剛果(金)Kalumbwe銅鈷礦床鉆孔分布及典型剖面Fig.2 Distribution of boreholes and typical sections of the Kalumbwe Cu-Co Deposit,D R Congo

3 分形分析與礦化品位域的確定

Kalumbwe礦床所有樣品數(shù)據(jù)(Cu=0.01%，為分析下限值)和工業(yè)礦體內(nèi)數(shù)據(jù)(按Cu=1.0%圈礦，未剔除夾石)的直方圖、概率圖以及分形統(tǒng)計結果如圖3所示。由圖3可知，直方圖和概率圖只能給出特高值，無法給出特低值的信息(圖3(a)～圖3(d))。直方圖顯示的特高值與統(tǒng)計數(shù)據(jù)范圍有關，斷點處分別為12%和13%(圖3(a)和圖3(b))。概率圖則存在多個拐點12.05%、14.09%和17.04%，需通過經(jīng)驗判斷哪個為特高值(圖3(c)和圖3(d))。用分形方法計算結果比較穩(wěn)定，特低值為0.3%，特高值為13.9%，分形區(qū)間為0.3%～13.9%(圖3(e)和圖3(f))。這一區(qū)間，即為去掉特異值后的“礦化品位域”。

圖3 樣品數(shù)據(jù)的統(tǒng)計和分形分布特征Fig.3 Statistical and fractal distribution characteristics of the sample data

不同于傳統(tǒng)方法確定礦體的品位域在工業(yè)品位和特高值之間(即所謂的“礦體”或工業(yè)品位域)，分形方法確定的是礦化品位域，域內(nèi)品位變化于礦化品位與特高值之間。如果采用礦化品位域作為最終的估算域，那么部分低于工業(yè)品位的數(shù)據(jù)也將參與資源儲量估算，這取決于工業(yè)品位域的邊界條件。

[7] Leo Suryadinata, “Indonesian Chinese Education: Past and Present”, Indonesia, Vol. 14, No. 34 (1972), p. 70.

4 邊界分析與估算域的確定

邊界分析目的是通過域的邊界條件分析，確定用哪些數(shù)據(jù)進行估值，形成最終的估算域。其方法是考察特定域邊界兩側(cè)一定距離內(nèi)品位數(shù)據(jù)的變化[11]。“硬”邊界(hard boundary)即域內(nèi)外的品位數(shù)據(jù)是一個突變降低的過程，說明域內(nèi)外品位連續(xù)性較差，估值時僅采用域內(nèi)的數(shù)據(jù)。“軟”邊界(soft boundary)即域內(nèi)外的品位數(shù)據(jù)是一個逐步降低的過程，說明域內(nèi)外品位連續(xù)性較好，估值時應采用域外一定范圍內(nèi)數(shù)據(jù)。“半軟”邊界(semi-soft boundary)則位于二者之間，一般也采用邊界外一定距離內(nèi)的數(shù)據(jù)參與估值。

本文重點分析礦化品位(邊界品位Cu=0.3%，WF=0.3%)和工業(yè)指標(邊界品位Cu=1.0%，WF=1.0%)圈定品位域的邊界條件。WF=0.3%時，域內(nèi)、域外在靠近邊界時分別呈下降和上升趨勢(圖4(a)；WF=1.0%時，也有類似特征(圖4(c))，兩者均為“半軟”邊界條件。不失一般性，工業(yè)品位域為軟邊界或半軟邊界條件時，應采用礦化域數(shù)據(jù)或邊界外一定范圍內(nèi)數(shù)據(jù)進行估值。本文采用WF=0.3%內(nèi)數(shù)據(jù)參與后續(xù)變異函數(shù)分析和估值流程。

圖4 礦化品位域和工業(yè)品位域的邊界分析Fig.4 Boundary analysis for the mineralization and industrial grade domains

5 變異函數(shù)分析與估值參數(shù)的確定

根據(jù)IDW估值原理，進行資源儲量估算前需要確定兩個關鍵參數(shù)，即冪次和橢球搜索半徑。根據(jù)經(jīng)驗，橢球初始搜索半徑按照探明工程間距的1.25倍[16]，而冪次通常主要參考估值域的品位變化系數(shù)[17]。變異函數(shù)是描述區(qū)域化變化空間變異性的工具。研究表明，可以根據(jù)變異函數(shù)分析中塊金效應的大小，來確定IDW冪次取值和橢球搜索半徑相對大小[18]。一般地，塊金效應越小，采用的估值冪次越高，橢球搜索半徑越接近變程。本文通過變異函數(shù)的井向變異函數(shù)和全向變異函數(shù)分析提供的塊金效應和變程信息，對IDW估值的距離和冪次參數(shù)進行約束。

Kalumbwe礦床礦化品位域(WF=0.3%)內(nèi)銅品位數(shù)據(jù)的變異函數(shù)分析結果如圖5所示。由圖5可知，其塊金值為0.35%2，基臺值為3.52%2，計算得到塊金效應為9.07%(圖5(a))；全向變異函數(shù)擬合的變程為a=59 m(圖5(b))。由于變異函數(shù)分析顯示出低塊金效應(<25%)，在后續(xù)估值流程中應采用相對高的冪次(4～5次)，橢球搜索半徑盡量選擇與變程一致(60 m)。

圖5 礦化域內(nèi)銅品位數(shù)據(jù)的井向與全向變異函數(shù)特征Fig.5 The characteristics of downhole directional and omnidirectional variogram of copper grade in the mineralized domain

6 IDW估值結果

基于工業(yè)指標(邊界Cu=1.0%)圈定的線框(WF=1.0%)，用礦化品位域數(shù)據(jù)(邊界品位Cu=0.3%，WF=0.3%)進行估值，考慮取2～5不同冪次的情況，結果如圖6所示。由圖6可知，剖面上礦化整體呈上下部強、中部弱的特征，礦化連續(xù)性方向與礦體產(chǎn)狀基本一致，表明三維橢球參數(shù)選擇適宜。從估值結果看，隨著冪次的增加，礦化的結構性增強，也就是高品位和低品位礦化的連續(xù)性都得到強化。從剖面中部夾石(Cu≤0.3%)的局部驗證情況明顯看出，高冪次(4～5次)時，估值的夾石范圍與原始鉆孔數(shù)據(jù)的夾石范圍基本一致；而低冪次(2～3次)時，估值的夾石范圍明顯較實際鉆孔數(shù)據(jù)范圍偏小。

圖6 不同冪次條件下礦塊模型Fig.6 Block models under different power conditions

7 討 論

7.1 工業(yè)品位域及其邊界條件對估值影響

前已述及，資源儲量估算首先要解決的是用哪些數(shù)據(jù)參與估值。從估值流程來看，是分2個層次或者尺度來實現(xiàn)的，即：①整個礦體范圍內(nèi)哪些數(shù)據(jù)可以用于估值(即估算域)；②對特定待估值塊，具體用哪些局部數(shù)據(jù)進行估值(即搜索橢球范圍)。考慮到礦山設計和開發(fā)，都需要利用一個具有經(jīng)濟性的礦體而不是礦化來進行，本文的討論均以工業(yè)品位域為對象。

由圖7可知，若工業(yè)品位域為硬邊界條件，僅采用工業(yè)品位域內(nèi)數(shù)據(jù)進行估值；但如果是軟邊界，應當適當擴大估算域范圍，而采用礦化品位域內(nèi)數(shù)據(jù)。估算域的范圍將直接影響后續(xù)搜索橢球內(nèi)樣品點的確認。盡管未知塊中心點均位于工業(yè)品位域內(nèi)，但搜索橢球仍可以根據(jù)給定的參與估值數(shù)據(jù)(或者邊界條件)，來確定哪些數(shù)據(jù)點參與其最終的估值(圖8)。

圖7 硬邊界及軟邊界的剖面品位變化Fig.7 Section grade variation of hard and soft boundaries

1-圓從大到小依次代表工業(yè)品位域、礦化品位域和低于礦化品位樣品點；2-參與估值樣品點；3-待估值塊；4-搜索橢球圖8 不同邊界條件對估值影響Fig.8 Impact of different boundary conditions on estimation(資料來源：據(jù)高幫飛[13]修改)

圖9呈現(xiàn)了Kalumbwe礦床WF=1.0%時，在不同邊界條件下，IDW估值結果與實際生產(chǎn)數(shù)據(jù)的對比結果。由圖9可知，軟邊界情況下，估值結果與生產(chǎn)數(shù)據(jù)整體趨勢吻合；但硬邊界情況下，不僅估值結果差距較大，趨勢也與生產(chǎn)數(shù)據(jù)有較大差異。這一結果與上文討論的，本礦床工業(yè)品位域邊界條件為“半軟”邊界的結論是一致的。此外，還可以明顯看出，采用硬邊界條件，估算品位比實際生產(chǎn)數(shù)據(jù)整體偏高。因而，不論在前期礦業(yè)項目投資收購，還是在后期設計開發(fā)階段，都會帶來隱性的經(jīng)營管理風險。

圖9 不同邊界條件估值結果與實際對比分析Fig.9 Comparison and analysis of different boundary conditions

7.2 塊金效應與冪次的確定

為了進一步評價Kalumbwe礦床不同冪次條件下的估值效果，將估值結果與礦山實際生產(chǎn)數(shù)據(jù)進行對比(圖10)。可以看出，低冪次條件下(2～3次)，品位總體被低估；而高冪次(4～5次)，品位總體與生產(chǎn)數(shù)據(jù)接近，分別計算不同冪次條件下的估計方差，見式(2)。

(2)

表1 不同冪次條件下的估值方差比較Table 1 Estimation variance comparison under different power conditions

為進一步評價改進IDW估值流程估值效果，將傳統(tǒng)IDW、改進IDW(冪次為5)、普通克里格三種估值方法估值結果與礦山實際生產(chǎn)數(shù)據(jù)進行對比分析估值效果，其中改進IDW(冪次為5)估值結果最優(yōu)(表2)。

表2 不同估值方法估值方差比較Table 2 Comparison of valuation variances of different estimation methods

7.3 改進的IDW估值流程

綜合上述，本文提出改進IDW估值流程(圖11)。①建立地質(zhì)數(shù)據(jù)庫；②對樣品數(shù)據(jù)進行統(tǒng)計分析；③進行分形分析，剔除特高值及特低值，確定礦化品位域；④分別圈定工業(yè)品位域線框，礦化品位域線框；⑤進行工業(yè)品位域的邊界分析，確定其邊界條件，硬邊界條件下估算域為工業(yè)品位域，軟(半軟)邊界條件下估值域為礦化品位域；⑥進行變異函數(shù)分析，確定估值兩個重要參數(shù)冪次及橢球搜索半徑；⑦進行IDW估值，估值線框選用工業(yè)礦體線框，根據(jù)工業(yè)品位域的邊界條件選擇那些數(shù)據(jù)參與估值，通過變程確定最優(yōu)搜索半徑，通過塊金效應大小確定冪次；⑧輸出資源儲量分類報告。

圖11 改進的IDW估值流程Fig.11 The improved IDW estimation workflow(注：灰色部分為本文修改內(nèi)容)

8 結 論

本文以剛果(金)Kalumbwe銅鈷礦床為例，進行了改進的IDW估值試驗，得出以下結論和認識。

1)分形分析是一種確定資源儲量估算品位域的有效方法。它在不割裂礦化連續(xù)性的基礎上，同時確定特高值和特低值，且不同尺度下均能給出一致的結果。

2)通過工業(yè)品位域的邊界條件分析，明確了參與估值數(shù)據(jù)與域的邊界條件之間的關系，硬邊界條件下全部采用工業(yè)品位域內(nèi)數(shù)據(jù)參與估值，否則采用礦化域數(shù)據(jù)參與估值。

3)引入克里格估值方法中的塊金效應原理作為IDW估值參數(shù)選取依據(jù)，對不同估值方法之間進行融合。改進的IDW估值流程在操作上可行，該流程不僅深化了樣品數(shù)據(jù)的地質(zhì)特征與礦化品位特征之間的聯(lián)系，同時還可以給出接近于實際生產(chǎn)數(shù)據(jù)的估值結果，值得在實際工作中進一步推廣使用。