克里格法在離子吸附型稀土礦勘查儲量估算中的應(yīng)用

趙 汀，王登紅，王瑞江，鄧茂春，陳為光

(1.中國地質(zhì)科學(xué)院礦產(chǎn)資源研究所，國土資源部成礦作用與資源評價重點實驗室，北京 100037；2.江西省地質(zhì)礦產(chǎn)勘查開發(fā)局贛南地質(zhì)調(diào)查大隊，江西 贛州 341000)

離子吸附型稀土礦是我國獨特的優(yōu)勢礦產(chǎn)資源，近年來取得了很大的找礦進展，但資源消耗也是很快的[1-2]。如何客觀評價我國離子吸附型稀土資源的儲量，乃當務(wù)之急。在我國地質(zhì)礦產(chǎn)勘查工作中，傳統(tǒng)的儲量估算方法有：算術(shù)平均法、地質(zhì)塊段法、多角形法、斷面法(包括垂直剖面法和水平斷面法)及距離平方反比法等[3]。這些傳統(tǒng)方法在我國多年的地質(zhì)找礦工作中已經(jīng)成為國家規(guī)范，絕大部分礦產(chǎn)勘查報告儲量是用傳統(tǒng)方法估算的，但存在的不足之處在于：①把計算儲量的主要參數(shù)(如：品位)視成固定變量，忽略了這類變量的隨機性；②在變量的空間相關(guān)關(guān)系上也僅僅考慮了距離，忽略了變量的地質(zhì)變異規(guī)律性；③無法給出估計精度。克里格插值法的區(qū)域化變量理論及其變差函數(shù)工具，正好有效地彌補了傳統(tǒng)的儲量估算方法的缺點。

1 克里格插值法的基本原理

克里格插值法是法國馬特隆教授以南非礦山地質(zhì)工程師D. G .克里格的名字命名的一種方法[4]。它是以區(qū)域化變量理論為基礎(chǔ)，以變差函數(shù)為主要工具，在保證估計值滿足無偏性條件和最小方差條件的前提下求得估計值。設(shè)區(qū)域化變量f(x)滿足二階平穩(wěn)假設(shè)或本征假設(shè)，則待插點P的估計值為：

式中，fi為n個已知點的函數(shù)值，wi為n個已知點的全系數(shù)。由無偏的條件下，有：

再根據(jù)估計的方差最小的條件:

式中，μ為拉格朗日算子，V(xj-xi)為已知點間的變差函數(shù)值；V(xp-xi)為已知點與待插點間的變差函數(shù)值。求出待插點P的估計值。簡單地說，克里格插值法就是一種特定的滑動加權(quán)平均法[5]。

克里格法是在不斷發(fā)展和完善之中，對于不同情況的條件可以運用相應(yīng)的克里格法。如：當滿足二階平穩(wěn)(或本征)假設(shè)時，若區(qū)域化變量的數(shù)學(xué)期望存在且已知的情況下，用簡單克里格法；若其存在而未知的情況下，用普通克里格法；當原始數(shù)據(jù)服從對數(shù)正態(tài)分布時，用對數(shù)正態(tài)克里格法；在非平穩(wěn)現(xiàn)象中，用泛克里格法；對無分布和非參數(shù)的區(qū)域化變量，用指示克里格法；在計算可采儲量時，用非線性估計量的析取克里格法。此外還有因子克里格法、概率克里格法等等。

2 克里格法在離子吸附型稀土礦儲量估算中的應(yīng)用

南方離子吸附型稀土礦礦體一般呈似層狀，分布于花崗巖風化殼全風化層中。全風化層的全部或部分即是礦體，礦體礦化連續(xù)性好，礦體厚度的穩(wěn)定程度主要受到地形影響。礦體厚度一般以山頂最厚，山脊厚度次之，山坡兩翼及坡腳礦體厚度較薄，礦體品位變化與地形地貌、風化殼厚度及風化殼中礦物的解離程度有著較為密切的關(guān)系。一般地，山腰或山頂偏下部位稀土品位較高，反之則較低；鉆孔顯示品位由上至下逐漸增加，然后在底部再逐漸降低的過程。

在江西贛州離子吸附型稀土礦勘查項目中，儲量計算普遍使用的是水平投影地質(zhì)塊段法，該方法對于礦體厚度、面積、品位的變化控制需要較高的探礦工程網(wǎng)度，需要較高的資金人力投入。本文以江西贛州某稀土礦為例，采用克里格估算方法對礦體頂?shù)装搴推肺贿M行插值計算，最后計算儲量，并與傳統(tǒng)的地質(zhì)塊段法估算的儲量進行對比，以期獲得更好的效果。

2.1 某稀土礦山三維建模

我國南方離子吸附型稀土礦，似層狀展布的離子吸附型稀土礦礦體被雨水切割侵蝕形成溝壑，在探礦工程不足的情況下，塊斷法估算儲量體積的估算誤差較大。本文基于我國南方離子吸附型稀土礦賦存的特點[6]設(shè)計了一套適合的克里格法計算流程(見圖1)，首先搜集地質(zhì)勘查鉆孔數(shù)據(jù)，按照克里格法所需數(shù)據(jù)和參數(shù)要求建設(shè)一個礦床鉆孔數(shù)據(jù)庫，數(shù)據(jù)項要包括平面直角坐標、品位、頂板、底板測量數(shù)據(jù)、比重等[7-9]。

圖1 克里格法儲量動態(tài)估算流程

常用的礦體三維建模技術(shù)包括數(shù)字地面模型(DTM)、實體模型、塊體模型和直棱柱模型等[10]。對于形態(tài)簡單的層狀礦體多采用DTM、直棱柱等模型對礦體幾何形態(tài)進行表達。離子吸附型稀土礦礦體受風化殼控制呈似層狀分布，故本次估算采用DTM模型作為建模方案。

變差函數(shù)設(shè)定是克里格法估算最重要的部分，因為變差函數(shù)反映的是礦化現(xiàn)象的空間相關(guān)性，稀土礦礦體在水平和垂向方向上的礦化現(xiàn)象存在一定的規(guī)律性，制定對應(yīng)的變差函數(shù)能對品位模擬起到關(guān)鍵作用。在常規(guī)儲量計算工作之后，由于克里格儲量估算法有可視、快速、靈活的特點，這樣通過設(shè)定不同的市場價格對價格儲量進行敏感性分析非常方便。

為了對克里格法估算結(jié)果進行驗證，選擇一個工作程度較高的稀土礦山，同時采用克里格法和塊斷法兩種方法進行儲量估算比較，以驗證新方法的正確性。

此稀土礦地處南嶺東西向復(fù)雜褶皺帶東段之南雄—周田一級構(gòu)造帶和新華夏系之廣昌—定南北北東向構(gòu)造帶復(fù)合部位，大地構(gòu)造位置屬華夏板塊之“武夷褶皺帶”的西坡南端。礦區(qū)內(nèi)大面積出露中生代酸性花崗巖——寨背巖體，寨背巖體為復(fù)式巖體，出露面積約254 km2，屬燕山早期第一階段的花崗巖。寨背巖體屬酸性富堿質(zhì)鋁過飽和巖石，具有高硅、富鉀特征，巖石副礦物主要有獨居石、磷釔礦、螢石、磁鐵礦和磷灰石等；稀土總量為343.82 μg/g，屬輕稀土富集型[11]。不同期次的花崗巖遭受長期的風化淋積作用，在其風化殼中形成規(guī)模不等的離子吸附型稀土礦床。

從圖2礦山DTM模型看出勘查區(qū)工作區(qū)屬低山丘陵地帶，高低起伏，坡度較緩，高程為400～500 m，坡度一般在40°左右，在 450～500 m之間共布置18個鉆孔。

圖2 某礦山克里格數(shù)字高程DTM模型

2.2 克里格法與塊斷法估算儲量的結(jié)果比較

2.2.1 稀土氧化物儲量估算方法

依據(jù)《稀土礦產(chǎn)地質(zhì)勘查規(guī)范》(DZ/T 0204—2002)，本區(qū)離子吸附型輕稀土礦床一般工業(yè)指標為：邊界品位(TR2O3) 0.05%，最低工業(yè)品位(TR2O3)0.08%，最小可采厚度1 m，夾石剔除厚度2 m。離子吸附型稀土礦稀土氧化物儲量的計算公式為：

Q=S×H×D

P=Q×C

式中，Q為礦石量(t)；S為礦塊面積(m2，計算機讀取)；H為礦體厚度(m)；D為礦石體重(t/m3)；C為TR2O3品位(%)；P為TR2O3(t)。

2.2.2 稀土氧化物儲量計算結(jié)果

結(jié)合各礦山2007年及2008年度儲量動態(tài)檢測報告提供的開采數(shù)據(jù)[11]，本區(qū)礦石體重為1.539 t/m3。通過三維礦體建模和克里格法計算，估計出整個礦體的品位分布(見圖3)。通過普通克里格法計算儲量，計算結(jié)果是：礦體體積322萬m3，比塊段法增加了34萬m3，變化率為11.8%；儲量估算結(jié)果為4332 t，比塊段法增加586 t，變化率為15%。與實際勘探值相比較，克里格法的計算結(jié)果基本合理，充分體現(xiàn)出克里格法快速、準確、方便的特點。

圖3 稀土離子吸附相品位等值線和礦體頂?shù)装鍒D

3 克里格法在儲量估算中的動態(tài)評價

3.1 稀土礦原地浸礦工藝和價格變化需要對邊界品位的動態(tài)評估

贛南稀土礦山多采用原地浸礦采礦工藝，在實際生產(chǎn)過程中如果收液溝槽位置過高，而不是位于礦體最低部位，殘留礦量較大可能造成資源浪費問題；或者收液溝槽布設(shè)不合理，使部分母液流失。而稀土的價格近幾年上漲較快，勘查時確定的邊界品位隨著價格提高的發(fā)生變化，原來一些不經(jīng)濟的含稀土的風化殼也可以合理利用，這樣通過動態(tài)快速的礦體圈定方法可以方便地輔助礦山調(diào)整采礦設(shè)計。

3.2 價格-邊界品位敏感性分析輔助克里格法快速動態(tài)評價資源

我國稀土供應(yīng)在世界市場上處于絕對領(lǐng)導(dǎo)地位，但是對于國際稀土價格長期喪失定價權(quán)，國內(nèi)稀土礦山對國際市場價格反映非常敏感，稀土價格像過山車一樣，變化起伏[12-14]，而邊界品位還沿用以往的數(shù)據(jù)，這將浪費大量的礦產(chǎn)資源。根據(jù)礦石的開采成本、可選性和市場價格，經(jīng)濟合理地確定礦石在不同價格下的動態(tài)邊界品位，動態(tài)調(diào)整礦山的資源量，提高資源的利用率，延長礦山服務(wù)年限，可實現(xiàn)礦山的可持續(xù)發(fā)展[15-17]，同時對科學(xué)開發(fā)我國稀土資源具有重要的戰(zhàn)略意義。

由上可知，本文需要建立邊界品位和價格之間的量化關(guān)系才能進行分析，以下利用盈虧平衡的思想建立了一個價格-邊界品位估算模型。

投入：礦山的成本以其他的投資金額。即：

β· (Cα+IαF)

產(chǎn)出：礦山的收益，銷售額除去稅金部分。即：

a·ρ·ε0·(Vc-Tc)

由于需要盈虧平衡，投入=產(chǎn)出，即：

β·(Ca+IaF)=a·ρ·ε0·(Vc-Tc)

所以，得到邊界品位(α，%)的計算模型為：

式中：β—精礦品位(%)；Cα—折算成本(元/噸)；Iα—年產(chǎn)一噸礦石的單位投資指標(元/噸)；F—資金還原系數(shù)；ρ—采礦回采率(%)；ε0—選礦回收率(%)；Vc—精礦產(chǎn)品價格(元/噸)；Tc—每噸精礦的稅金(元/噸)。

將計算模型視為變量α和Vc的函數(shù)關(guān)系式，其他量視為定量參數(shù)(可根據(jù)各自礦山的情況獲得相應(yīng)確定值)，由此可通過此模型進行價格敏感性分析。

通過克里格法三維儲量估算軟件，按照市場價格靈活設(shè)定邊界品位，圈定礦體邊界，設(shè)置采礦工程位置，使得資源利用效率最大化。例如，當前市場上精礦價格為10萬元/噸時，運用價格-邊界品位敏感性分析模塊,求出企業(yè)凈利潤為零時的邊界品位為0.08%。圖4中等值線為“0.08%”的等值線圈定區(qū)域是此價格下經(jīng)濟可采的礦體，隨著精礦的市場價格上漲到12萬元/噸，邊界品位降低至0.07%時，礦體圈定面積增加了15%，探明資源儲量相應(yīng)增加，而價格隨市場發(fā)生下降時，邊界品位上升到0.09%，礦體圈定面積則相應(yīng)減少了10%。

圖4 價格-邊界品位敏感性分析

贛南離子吸附型稀土礦山開采方法普遍采用原地浸礦工藝，收集稀土母液的導(dǎo)流孔組成的收液工程控制面要與礦體底部構(gòu)成的空間曲面處于同一位置或略低位置，這樣才能達到最優(yōu)的采礦設(shè)計方案。在以克里格法為基礎(chǔ)的三維估算系統(tǒng)中動態(tài)設(shè)置邊界品位，采用礦山生產(chǎn)盈虧平衡法作價格敏感性分析，可以根據(jù)當時的稀土價格靈活地圈定礦體邊界，并能迅速地計算出所圈定礦塊的品位分布、體重、礦石量、金屬量等，幫助礦山選擇合理的采礦工程布置，保證了礦產(chǎn)資源的合理利用。

4 結(jié)語

根據(jù)贛南某離子吸附型稀土礦的地質(zhì)特征，研究了離子吸附型稀土礦礦體展布、品位、厚度、穩(wěn)定性等特點；根據(jù)礦山勘查資料建立了礦體三維模型，使得礦體的形態(tài)更加直觀，方便后期的勘查工程的部署。在三維模型的基礎(chǔ)上，應(yīng)用克里格法對礦山地質(zhì)經(jīng)濟參數(shù)(品位、厚度、分布)進行了無偏線形估計，采用地質(zhì)統(tǒng)計學(xué)軟件求取了礦山資源儲量。與傳統(tǒng)的地質(zhì)塊段方法相比較，克里格法具有計算速度快、精度高的特點，更加符合實際。

在克里格地質(zhì)統(tǒng)計學(xué)方法基礎(chǔ)上結(jié)合技術(shù)經(jīng)濟評價模型，并編寫軟件具體實現(xiàn)，這是本研究的一個特色，對于采用原地浸礦工藝的礦山，其工程部署也要適當依據(jù)市場價格調(diào)整才能更加充分地利用資源，幫助礦山企業(yè)及時掌握礦山資源儲量動態(tài)價值，通過研究礦山生產(chǎn)盈虧平衡條件,設(shè)置價格敏感性分析模型條件，輸入最新稀土價格靈活地計算出最低邊界品位，并在以克里格法為基礎(chǔ)的三維估算系統(tǒng)中動態(tài)設(shè)置邊界品位，從而迅速地動態(tài)圈定出經(jīng)濟可供開采的礦塊及其品位分布、礦石量、金屬量等。可見，本文建立的以克里格法為基礎(chǔ)的三維估算系統(tǒng)能夠充分滿足礦山動態(tài)管理的需要。

致謝：江西省地質(zhì)礦產(chǎn)勘查開發(fā)局贛南地質(zhì)調(diào)查大隊提供的野外數(shù)據(jù)對此研究幫助極大，在此深表感謝！

5 參考文獻

[1] 王登紅，王瑞江，李建康，趙芝，于揚，代晶晶，陳鄭輝，李德先，屈文俊，鄧茂春，付小方，孫艷，鄭國棟.中國三稀礦產(chǎn)資源戰(zhàn)略調(diào)查研究進展綜述[J].中國地質(zhì)，2013,40(2)：361-370.

[2] 袁忠信，李建康，王登紅，鄭國棟，婁德波，陳鄭輝，趙芝，于揚.中國稀土礦床成礦規(guī)律[M].北京：地質(zhì)出版社，2013：1-116.

[3] DZ/T 0204—2002，稀土礦產(chǎn)地質(zhì)勘查規(guī)范[M].2003.

[4] 王仁鐸，胡光道.線性地質(zhì)統(tǒng)計學(xué)[M].北京：地質(zhì)出版社，1988：4-10.

[5] 陳歡歡，李星，丁文秀.Surfer 8.0等值線繪制中的十二種插值方法[J].工程地球物理學(xué)報，2007，4(1)：52-57.

[6] 鄧茂春，王登紅，曾載淋，張永忠，趙芝，鄒新勇，陳斌鋒.風化殼離子吸附型稀土圈礦方法評價[J].巖礦測試，2013，32(5)：803 -809.

[7] 陳鄭輝，肖克炎，吳健生，朱裕生.阿舍勒銅鋅礦床三維模型的研究[J].地質(zhì)論評，2000,46(增刊): 319-324.

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

[9] 楊東來.地質(zhì)體三維建模方法與技術(shù)指南[M].北京：地質(zhì)出版社，2007：65-69.

[10] 劉海英，劉修國，李超嶺.基于地質(zhì)統(tǒng)計學(xué)法的三維儲量估算系統(tǒng)研究與應(yīng)用[J].吉林大學(xué)學(xué)報(地球科學(xué)版)，2009，39(3)：541-546.

[11] 鄧茂春，艾生根，劉明發(fā).江西省定南縣嶺北稀土礦核查礦區(qū)資源儲量核查報告[R].贛州：贛南地質(zhì)大隊，2010.

[12] 王登紅，趙芝，于揚，趙汀，李建康，代晶晶，劉新星，何晗晗.離子吸附型稀土資源研究進展存在問題及今后研究方向[J].巖礦測試，2013，32(5)：796-802.

[13] 章文靖.稀土資源管理政策的國際比較及對中國的啟示[D].上海：東華大學(xué)，2012：22-63.

[14] 李皓，管宏平.中國建立稀土戰(zhàn)略儲備制度的國際戰(zhàn)略意義[J].河北青年管理干部學(xué)院學(xué)報，2009(5)：89-92.

[15] 汪朝，趙勇，賈明濤，肖英才，李寧.盈虧平衡的思想建立邊界品位估算模型[J].稀土，2011，35(4)：22-25.

[16] 李章林，張夏林，翁正平.指示克里格法在礦體圈定中的應(yīng)用研究[J].礦業(yè)快報，2008(465)：11-15.

[17] 姚筱煌，周雪珊.礦產(chǎn)資源經(jīng)濟評價及其程序[J].鈾礦冶,1996，10(2)：1-7.