非線性克里格法在礦體儲量估算中的應(yīng)用

劉修國，尤超，張唯

(中國地質(zhì)大學(xué) 信息工程學(xué)院，湖北 武漢，430074)

在采礦過程中，不可避免地伴隨著一些不確定因素[1]。這些不確定因素將會對礦床儲量品位的估值有嚴(yán)重影響。可靠的礦產(chǎn)資源和礦床儲量的估計對于所有的采礦作業(yè)都是十分重要的，不合理的礦產(chǎn)資源分類對項目的可行性有重大的影響。地質(zhì)統(tǒng)計學(xué)方法通?；陔S機函數(shù)而同時考慮樣本數(shù)據(jù)的空間關(guān)系[2]，所以，它能夠有效地對儲量參數(shù)(品位、厚度和累積量)空間分布進(jìn)行研究。地質(zhì)統(tǒng)計學(xué)方法一般被分為線性克里格，非線性克里格、非參數(shù)克里格和隨機模擬等[3]。由于線性克里格的平滑效應(yīng)，它有時不能準(zhǔn)確地預(yù)測待估點屬性值是否超過某一閾值。相反地，非線性克里格與非參數(shù)克里格方法(多元高斯克里格和指示克里格)則可以利用條件概率分布對這些待估點進(jìn)行較準(zhǔn)確估值[4]。本文作者將選擇某銅礦區(qū)實際數(shù)據(jù)作為實驗數(shù)據(jù)，并利用指示克里格和多元高斯克里格2種非線性克里格方法進(jìn)行研究與應(yīng)用，研究主要包括 2部分內(nèi)容：(1) 通過指示克里格對礦體的邊界進(jìn)行圈定；(2) 利用多元高斯克里格估算礦體的品位值。為了進(jìn)行對比，將應(yīng)用有礦化帶邊界的普通克里格的估值方法和無礦化帶邊界的多元高斯克里格估值方法對同一礦體進(jìn)行品位估值，最后對估值結(jié)果進(jìn)行討論。

1 基本原理

1.1 指示克里格估值方法

指示克里格方法其實質(zhì)就是通過一系列指示閾值，求出相應(yīng)閾值對應(yīng)的條件累積概率分布，并用條件累積概率分布表示位置點Z(u)的不確定性概率模型。指示克里格估值首先需要對數(shù)據(jù)進(jìn)行指示變換，即先選擇1個閾值，一般為礦體的邊界品位Z。其指示化的公式如下：

指示變換后得到一組以0和1形式組織的指示化數(shù)據(jù)。利用這組指示化數(shù)據(jù)，可對待估區(qū)域內(nèi)平均值在某一范圍中的概率進(jìn)行估計，其原理如下。

設(shè)在待估位置u處出現(xiàn)閾值小于或等于Z的概率為F(u,z)，P表示概率，E(u,z)表示在待估位置u處出現(xiàn)閾值小于或等于z的數(shù)學(xué)期望，可得到關(guān)系式如下：

根據(jù)數(shù)學(xué)期望的求取公式可得：

其中：F(u,z)為未知點u處出現(xiàn)品位zu小于或等于閾值Z的累積概率。

1.2 多元高斯克里格估值方法

多元高斯克里格估值方法是對原始樣品數(shù)據(jù)進(jìn)行正態(tài)變換，用變換后的數(shù)據(jù)進(jìn)行線性克里格估值，再對結(jié)果進(jìn)行逆變換。正態(tài)變換后y(ua) =φ(z(ua))被定義為正態(tài)點數(shù)[5]。多元高斯克里格估值方法是建立在多元高斯分布假設(shè)基礎(chǔ)之上的。原始變量在待估位置u處的條件累積分布函數(shù)(ccdf)為

在多元高斯克里格模型下，在位置u處ccdf的平均值和方差分別等同于正態(tài)點數(shù)y(ua)的普通克里格估值的均值(u)和方差(u)。因此，在位置u處的條件累積分布函數(shù)為[6]

式中：μOK(u)為拉格朗日乘數(shù)；普通克里格權(quán)重(u)和誤差(u)可由普通克里格方法計算得到。Y(u)的期望值可表示為

式中的估計值即為地質(zhì)統(tǒng)計學(xué)多元高斯克里格的估計值。

2 研究礦區(qū)概況

研究礦床是以銅為主的銅鋅硫化物礦床。該銅鋅礦共由2個礦體組成，編號分別為Ⅰ和Ⅱ礦體。Ⅰ礦體總體呈南北向展布，總體形態(tài)為透鏡狀，呈隱伏—半隱伏狀產(chǎn)出，地表有鐵帽覆蓋；產(chǎn)于細(xì)碧巖與石英角斑質(zhì)火山碎屑巖之間的接觸面上[7], 由于礦體類型不同，故礦體品位的空間連續(xù)性也大不相同，所以，對于不同礦體應(yīng)分別采用不同模型進(jìn)行儲量估算。由于Ⅰ礦體倒轉(zhuǎn)翼厚度較正常翼大，礦體厚度變化比較穩(wěn)定[8]，故本文以Ⅰ礦體的倒轉(zhuǎn)翼為研究對象。

Ⅰ礦體倒轉(zhuǎn)翼原始取樣數(shù)據(jù)為3 271個，銅品位最小值為0，最大值為14.85%，平均品位為2.37%。圖1所示為原始樣品數(shù)據(jù)的直方圖。從圖1可看出：原始數(shù)據(jù)為單峰分布，并且存在拖尾現(xiàn)象。

圖1 原始數(shù)據(jù)樣品直方圖Fig.1 Histogram of original sample

3 圈定礦體邊界

在進(jìn)行礦體儲量估算之前，首先需要圈定礦體。 在這種情況下，指示概率模型可用于確定礦體的邊界；利用指示概率模型進(jìn)行礦體邊界的確定包括 2個主要步驟：(1) 對數(shù)據(jù)進(jìn)行指示變換，利用變換后的數(shù)據(jù)采用普通克里格方法計算待估點的條件累積概率函數(shù)的估計值，也就是每一個待估點大于指示品位的概率；(2) 根據(jù)已經(jīng)得到的每個待估點大于邊界品位的概率，并同時考慮礦體的連續(xù)性，確定1個用于標(biāo)識礦和非礦的閾值對礦體進(jìn)行圈定。

3.1 指示概率模型

將概率模型應(yīng)用于確定礦體邊界之前需要對數(shù)據(jù)進(jìn)行指示變換，指示變換后得到一組以0和1形式組織的指示化數(shù)據(jù)。利用這組指示化數(shù)據(jù)，可對待估區(qū)域內(nèi)平均值在某一范圍內(nèi)的概率進(jìn)行估計。本研究將邊界品位 0.3%作為指示品位，將小于或等于 0.3%的品位設(shè)定值為1，而對大于0.3%的品位設(shè)定值為0。然后，利用這組指示化數(shù)據(jù)進(jìn)行實驗變差函數(shù)分析。在礦體主要的各向異性方向上，理論變差函數(shù)與實驗變差函數(shù)套合圖形如圖2所示。

3.2 定義閾值

對礦體進(jìn)行塊體劃分，塊段長×寬×高設(shè)置為10 m×10 m×10 m。利用指示化數(shù)據(jù)，根據(jù)得到的實驗變差函數(shù)進(jìn)行普通克里格估值，每個塊體估值后的結(jié)果為0～1之間的值。在指示克里格估值之后，需要確定一個閾值來劃定塊體是否屬于礦化塊。本研究取塊體品位大于0.3%的概率為0.25的水平值作為閾值。將指示值大于0.25的塊體歸入的塊體類；將指示值小于0.25的塊體舍去。此方法圈定的礦體保持了空間分布性與地質(zhì)成礦作用的特征。用指示概率模型得到的礦體邊界三維模型如圖3所示。

4 礦體儲量估算

用指示概率模型確定了礦體邊界之后，再利用多元高斯克里格對礦體儲量進(jìn)行估值。礦產(chǎn)資源儲量估算可按如下步驟進(jìn)行。

4.1 正態(tài)點數(shù)變換

圖2 指示變差函數(shù)圖Fig.2 IK variogram map

圖3 礦體邊界三維模型Fig.3 3D domains model

4.2 正態(tài)點數(shù)數(shù)據(jù)的變差函數(shù)分析

正態(tài)點數(shù)數(shù)據(jù)用于連續(xù)性的分析中。首先，從主要的各向異性方性計算正態(tài)點數(shù)數(shù)據(jù)的實驗變差函數(shù)。其中，理論變差函數(shù)與實驗變差函數(shù)的擬合圖形如圖4所示。根據(jù)正態(tài)點數(shù)變換的定義，各方向上的理論變差函數(shù)的基臺值都等于1。

4.3 高斯克里格礦體儲量估算

根據(jù)3個各向異性方向上的理論變差函數(shù)，用高斯克里格估值方法對礦體儲量進(jìn)行估算。同樣，塊段長×寬×高設(shè)置為10 m×10 m×10 m，其中搜索橢球的方位角為110°，傾伏角為20°，傾角為0°。用已知的數(shù)據(jù)進(jìn)行交叉驗證分析，誤差平均值為0.032，標(biāo)準(zhǔn)中誤差為1.011，由此可說明多元高斯克里格估值過程中變差函數(shù)的參數(shù)選擇較合理。

5 結(jié)果與驗證

按照不同的邊際品位，用多元高斯克里格法分別計算礦石儲量，結(jié)果見表1。

為了比較，分別將有礦化帶邊界的情況下用普通克里格方法和無礦化帶邊界的情況下用多元高斯克里格方法對同一礦體進(jìn)行估值計算，得到平均品位—邊際品位曲線對比圖，如圖5所示。綜合分析圖5可以看出：用普通克里格估值方法得到的平均品位明顯低于多元高斯克里格估值方法得到的平均品位(有礦體邊界和無礦體邊界)；而在低邊際品位和高品位位置，普通克里格方法得到的平均品位與多元高斯克里格方法得到的平均品位有更大差異；同時，在邊際品位大于 1.5%時無礦體邊界的多元高斯克里格方法與有礦體邊界的多元高斯克里格方法得到的平均品位較相近。從對比結(jié)果可以看出：由有礦化帶邊界的多元高斯克里格估值方法得到的估值礦體儲量更接近于實際的礦體儲量。

圖4 正態(tài)點數(shù)數(shù)據(jù)各向異性方向上的變差函數(shù)Fig.4 Main anisotropy directions of semi-variogram of normal scores data

表1 多元高斯克里格估值結(jié)果Table 1 Results of multi-Gaussian Kriging

圖5 平均品位—邊際品位曲線對比圖Fig.5 Comparison of result of CutOff-Grade

6 結(jié)論

(1) 指示概率模型是圈定礦體的一種有效方法。

(2) 無礦體邊界的估值會因為無限制的平滑效應(yīng)而增大估值結(jié)果的誤差。

(3) 普通克里格的非平穩(wěn)性和平滑效應(yīng)在某種程度上也會影響儲量估值的準(zhǔn)確性，而有礦化帶邊界的多元高斯克里格方法得到的儲量估值更加接近于實際礦體儲量，從而驗證了多元高斯克里格是一個簡單有效的儲量估值方法。

[1] Bastante F G, Taboada J, Alejano L R. Evaluation of the resources of a slate deposit using indicator Kriging[J].Engineering Geology, 2005, 81(4): 407-418.

[2] Ahmet D A G, Bayram A. Evaluating thickness of bauxite deposit using indicator geostatistics and fuzzy estimation[J].Resource Geology, 2008, 58(2): 188-195.

[3] Glacken I M, Snowden D V. Mineral resource estimation, in mineral resource and ore reserve estimation: The AussIMM guide to good practice[M]. Melbourne: The Australian Institute of Mining and Metallurgy, 2001: 189-198.

[4] Emery X. Uncertainty modeling and spatial prediction by multi-Gaussian Kriging: Accounting for an unknown mean value[J]. Computers & Geosciences, 2008, 34(11): 1431-1442.

[5] Goovaerts P. Comparative performance of indicator algorithms for modeling conditional probability distribution functions[J].Mathematical Geology, 1994, 26(3): 389-411.

[6] Emery X. Ordinary multigaussian Kriging for mapping conditional probabilities of soil properties[J]. Geoderma, 2006,132(1/2): 75-88.

[7] 張新宇, 肖克炎, 劉光勝, 等. 阿舍勒銅礦可視化儲量計算的指示克里格法應(yīng)用研究[J]. 吉林大學(xué)學(xué)報: 地球科學(xué)版, 2006,36(2): 305-308.

ZHANG Xinyu, XIAO Keyan, LIU Guangsheng, et al. The applied research of indicator Kriging method in visualized reserve computing of Ashele copper mine[J]. Journal of Jilin University: Earth Science Edition, 2006, 36(2): 305-308.

[8] 李章林, 張夏林, 翁正平. 指示克里格在礦體儲量計算方面的研究與應(yīng)用[J]. 礦業(yè)快報, 2008, 20(1): 11-15.

LI Zhanglin, ZHANG Xialin, WENG Zhengping. Research and application of indicator Kriging method in ore body reserves calculation[J]. Express Information of Mineing Industry, 2008,20(1): 11-15.

[9] 黃競先, 侯景儒. 泛克立格法和指示克立格法在地球化學(xué)探礦中的應(yīng)用[J]. 地球科學(xué): 中國地質(zhì)大學(xué)學(xué)報, 1994, 19(3):321-328.

HUANG Jinxian, HOU Jinru. Application of universal Kriging and indicator Kriging for geochemical prospecting[J]. Earth Sciende: Journal of China University of Geosciendes, 1994,19(3): 321-328.

[10] 劉修國, 李超嶺, 倪平澤, 等. 可視化固體礦產(chǎn)儲量估算方法研究[J]. 中國地質(zhì), 2007, 34: 407-412.

LIU Xiuguo, LI Chaolin, NI Pingze, et al. Study on visualization of solid mineral reserve estimation[J]. Geology in China, 2007,34: 407-412.

[11] Abzalov M Z, Humpherys M. Resource estimation of structurally complex and discontinuous mineralization using non-linear geostatitisics: Case study of a nesothermal gold deposit in northern Canada[J]. Explor Mining Geol, 2002,11(1/2/3/4): 19-29.

[12] 孫玉建. 地質(zhì)統(tǒng)計學(xué)在固體礦產(chǎn)資源評價中的若干問題研究[D]. 北京: 中國地質(zhì)大學(xué)地球科學(xué)與資源學(xué)院, 2008: 51-68.

SUN Yujian. A study on several issues on application of geostatistics in solid mineral resources estimation[D]. Beijing:China University of Geosciences. School of Earth Sciences and Resources, 2008: 51-68.

[13] 王全明. 我國銅礦勘查程度及資源潛力預(yù)測[D]. 北京: 中國地質(zhì)大學(xué)地球科學(xué)與資源學(xué)院, 2005: 99-129.

WANG Quanming. Degrees of copper ore exploration and resources potential of copper in lower degree areas in China[D].Beijing: China University of Geosciences. School of Earth Science and Resources, 2005: 99-129.

[14] 侯景儒, 黃競先. 地質(zhì)統(tǒng)計學(xué)在固體礦產(chǎn)資源/儲量分類中的應(yīng)用[J]. 地質(zhì)與勘探, 2001, 37(6): 61-66.

HOU Jingru, HUANG jingxian. Application of geosatistics in classification for resource/reserves of solid fuels and mineral commodities[J]. Geology and Prospecting, 2001, 37(6): 61-66.

[15] Bastante F G, Ordonez C, Taboada J, et al. Comparison of indicator Kriging, conditional indicator simulation and multiple-point statistics used to model slate deposits[J].Engineering Geology, 2008, 98(1/2): 50-59.