克里格法在离子吸附型稀土矿勘查储量估算中的应用

2014-08-01王登红王瑞江邓茂春陈为光

岩矿测试 2014年1期

赵汀，王登红，王瑞江，邓茂春，陈为光

(1.中国地质科学院矿产资源研究所，国土资源部成矿作用与资源评价重点实验室，北京 100037；2.江西省地质矿产勘查开发局赣南地质调查大队，江西赣州 341000)

离子吸附型稀土矿是我国独特的优势矿产资源，近年来取得了很大的找矿进展，但资源消耗也是很快的[1-2]。如何客观评价我国离子吸附型稀土资源的储量，乃当务之急。在我国地质矿产勘查工作中，传统的储量估算方法有：算术平均法、地质块段法、多角形法、断面法(包括垂直剖面法和水平断面法)及距离平方反比法等[3]。这些传统方法在我国多年的地质找矿工作中已经成为国家规范，绝大部分矿产勘查报告储量是用传统方法估算的，但存在的不足之处在于：①把计算储量的主要参数(如：品位)视成固定变量，忽略了这类变量的随机性；②在变量的空间相关关系上也仅仅考虑了距离，忽略了变量的地质变异规律性；③无法给出估计精度。克里格插值法的区域化变量理论及其变差函数工具，正好有效地弥补了传统的储量估算方法的缺点。

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

克里格插值法是法国马特隆教授以南非矿山地质工程师D. G .克里格的名字命名的一种方法[4]。它是以区域化变量理论为基础，以变差函数为主要工具，在保证估计值满足无偏性条件和最小方差条件的前提下求得估计值。设区域化变量f(x)满足二阶平稳假设或本征假设，则待插点P的估计值为：

式中，fi为n个已知点的函数值，wi为n个已知点的全系数。由无偏的条件下，有：

再根据估计的方差最小的条件:

式中，μ为拉格朗日算子，V(xj-xi)为已知点间的变差函数值；V(xp-xi)为已知点与待插点间的变差函数值。求出待插点P的估计值。简单地说，克里格插值法就是一种特定的滑动加权平均法[5]。

克里格法是在不断发展和完善之中，对于不同情况的条件可以运用相应的克里格法。如：当满足二阶平稳(或本征)假设时，若区域化变量的数学期望存在且已知的情况下，用简单克里格法；若其存在而未知的情况下，用普通克里格法；当原始数据服从对数正态分布时，用对数正态克里格法；在非平稳现象中，用泛克里格法；对无分布和非参数的区域化变量，用指示克里格法；在计算可采储量时，用非线性估计量的析取克里格法。此外还有因子克里格法、概率克里格法等等。

2 克里格法在离子吸附型稀土矿储量估算中的应用

南方离子吸附型稀土矿矿体一般呈似层状，分布于花岗岩风化壳全风化层中。全风化层的全部或部分即是矿体，矿体矿化连续性好，矿体厚度的稳定程度主要受到地形影响。矿体厚度一般以山顶最厚，山脊厚度次之，山坡两翼及坡脚矿体厚度较薄，矿体品位变化与地形地貌、风化壳厚度及风化壳中矿物的解离程度有着较为密切的关系。一般地，山腰或山顶偏下部位稀土品位较高，反之则较低；钻孔显示品位由上至下逐渐增加，然后在底部再逐渐降低的过程。

在江西赣州离子吸附型稀土矿勘查项目中，储量计算普遍使用的是水平投影地质块段法，该方法对于矿体厚度、面积、品位的变化控制需要较高的探矿工程网度，需要较高的资金人力投入。本文以江西赣州某稀土矿为例，采用克里格估算方法对矿体顶底板和品位进行插值计算，最后计算储量，并与传统的地质块段法估算的储量进行对比，以期获得更好的效果。

2.1 某稀土矿山三维建模

我国南方离子吸附型稀土矿，似层状展布的离子吸附型稀土矿矿体被雨水切割侵蚀形成沟壑，在探矿工程不足的情况下，块断法估算储量体积的估算误差较大。本文基于我国南方离子吸附型稀土矿赋存的特点[6]设计了一套适合的克里格法计算流程(见图1)，首先搜集地质勘查钻孔数据，按照克里格法所需数据和参数要求建设一个矿床钻孔数据库，数据项要包括平面直角坐标、品位、顶板、底板测量数据、比重等[7-9]。

图1 克里格法储量动态估算流程

常用的矿体三维建模技术包括数字地面模型(DTM)、实体模型、块体模型和直棱柱模型等[10]。对于形态简单的层状矿体多采用DTM、直棱柱等模型对矿体几何形态进行表达。离子吸附型稀土矿矿体受风化壳控制呈似层状分布，故本次估算采用DTM模型作为建模方案。

变差函数设定是克里格法估算最重要的部分，因为变差函数反映的是矿化现象的空间相关性，稀土矿矿体在水平和垂向方向上的矿化现象存在一定的规律性，制定对应的变差函数能对品位模拟起到关键作用。在常规储量计算工作之后，由于克里格储量估算法有可视、快速、灵活的特点，这样通过设定不同的市场价格对价格储量进行敏感性分析非常方便。

为了对克里格法估算结果进行验证，选择一个工作程度较高的稀土矿山，同时采用克里格法和块断法两种方法进行储量估算比较，以验证新方法的正确性。

此稀土矿地处南岭东西向复杂褶皱带东段之南雄—周田一级构造带和新华夏系之广昌—定南北北东向构造带复合部位，大地构造位置属华夏板块之“武夷褶皱带”的西坡南端。矿区内大面积出露中生代酸性花岗岩——寨背岩体，寨背岩体为复式岩体，出露面积约254 km2，属燕山早期第一阶段的花岗岩。寨背岩体属酸性富碱质铝过饱和岩石，具有高硅、富钾特征，岩石副矿物主要有独居石、磷钇矿、萤石、磁铁矿和磷灰石等；稀土总量为343.82 μg/g，属轻稀土富集型[11]。不同期次的花岗岩遭受长期的风化淋积作用，在其风化壳中形成规模不等的离子吸附型稀土矿床。

从图2矿山DTM模型看出勘查区工作区属低山丘陵地带，高低起伏，坡度较缓，高程为400～500 m，坡度一般在40°左右，在 450～500 m之间共布置18个钻孔。

图2 某矿山克里格数字高程DTM模型

2.2 克里格法与块断法估算储量的结果比较

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

依据《稀土矿产地质勘查规范》(DZ/T 0204—2002)，本区离子吸附型轻稀土矿床一般工业指标为：边界品位(TR2O3) 0.05%，最低工业品位(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 稀土氧化物储量计算结果

结合各矿山2007年及2008年度储量动态检测报告提供的开采数据[11]，本区矿石体重为1.539 t/m3。通过三维矿体建模和克里格法计算，估计出整个矿体的品位分布(见图3)。通过普通克里格法计算储量，计算结果是：矿体体积322万m3，比块段法增加了34万m3，变化率为11.8%；储量估算结果为4332 t，比块段法增加586 t，变化率为15%。与实际勘探值相比较，克里格法的计算结果基本合理，充分体现出克里格法快速、准确、方便的特点。

图3 稀土离子吸附相品位等值线和矿体顶底板图

3 克里格法在储量估算中的动态评价

3.1 稀土矿原地浸矿工艺和价格变化需要对边界品位的动态评估

赣南稀土矿山多采用原地浸矿采矿工艺，在实际生产过程中如果收液沟槽位置过高，而不是位于矿体最低部位，残留矿量较大可能造成资源浪费问题；或者收液沟槽布设不合理，使部分母液流失。而稀土的价格近几年上涨较快，勘查时确定的边界品位随着价格提高的发生变化，原来一些不经济的含稀土的风化壳也可以合理利用，这样通过动态快速的矿体圈定方法可以方便地辅助矿山调整采矿设计。

3.2 价格-边界品位敏感性分析辅助克里格法快速动态评价资源

我国稀土供应在世界市场上处于绝对领导地位，但是对于国际稀土价格长期丧失定价权，国内稀土矿山对国际市场价格反映非常敏感，稀土价格像过山车一样，变化起伏[12-14]，而边界品位还沿用以往的数据，这将浪费大量的矿产资源。根据矿石的开采成本、可选性和市场价格，经济合理地确定矿石在不同价格下的动态边界品位，动态调整矿山的资源量，提高资源的利用率，延长矿山服务年限，可实现矿山的可持续发展[15-17]，同时对科学开发我国稀土资源具有重要的战略意义。

由上可知，本文需要建立边界品位和价格之间的量化关系才能进行分析，以下利用盈亏平衡的思想建立了一个价格-边界品位估算模型。

投入：矿山的成本以其他的投资金额。即：

β· (Cα+IαF)

产出：矿山的收益，销售额除去税金部分。即：

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

由于需要盈亏平衡，投入=产出，即：

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

所以，得到边界品位(α，%)的计算模型为：

式中：β—精矿品位(%)；Cα—折算成本(元/吨)；Iα—年产一吨矿石的单位投资指标(元/吨)；F—资金还原系数；ρ—采矿回采率(%)；ε0—选矿回收率(%)；Vc—精矿产品价格(元/吨)；Tc—每吨精矿的税金(元/吨)。

将计算模型视为变量α和Vc的函数关系式，其他量视为定量参数(可根据各自矿山的情况获得相应确定值)，由此可通过此模型进行价格敏感性分析。

通过克里格法三维储量估算软件，按照市场价格灵活设定边界品位，圈定矿体边界，设置采矿工程位置，使得资源利用效率最大化。例如，当前市场上精矿价格为10万元/吨时，运用价格-边界品位敏感性分析模块,求出企业净利润为零时的边界品位为0.08%。图4中等值线为“0.08%”的等值线圈定区域是此价格下经济可采的矿体，随着精矿的市场价格上涨到12万元/吨，边界品位降低至0.07%时，矿体圈定面积增加了15%，探明资源储量相应增加，而价格随市场发生下降时，边界品位上升到0.09%，矿体圈定面积则相应减少了10%。

图4 价格-边界品位敏感性分析

赣南离子吸附型稀土矿山开采方法普遍采用原地浸矿工艺，收集稀土母液的导流孔组成的收液工程控制面要与矿体底部构成的空间曲面处于同一位置或略低位置，这样才能达到最优的采矿设计方案。在以克里格法为基础的三维估算系统中动态设置边界品位，采用矿山生产盈亏平衡法作价格敏感性分析，可以根据当时的稀土价格灵活地圈定矿体边界，并能迅速地计算出所圈定矿块的品位分布、体重、矿石量、金属量等，帮助矿山选择合理的采矿工程布置，保证了矿产资源的合理利用。

4 结语

根据赣南某离子吸附型稀土矿的地质特征，研究了离子吸附型稀土矿矿体展布、品位、厚度、稳定性等特点；根据矿山勘查资料建立了矿体三维模型，使得矿体的形态更加直观，方便后期的勘查工程的部署。在三维模型的基础上，应用克里格法对矿山地质经济参数(品位、厚度、分布)进行了无偏线形估计，采用地质统计学软件求取了矿山资源储量。与传统的地质块段方法相比较，克里格法具有计算速度快、精度高的特点，更加符合实际。

在克里格地质统计学方法基础上结合技术经济评价模型，并编写软件具体实现，这是本研究的一个特色，对于采用原地浸矿工艺的矿山，其工程部署也要适当依据市场价格调整才能更加充分地利用资源，帮助矿山企业及时掌握矿山资源储量动态价值，通过研究矿山生产盈亏平衡条件,设置价格敏感性分析模型条件，输入最新稀土价格灵活地计算出最低边界品位，并在以克里格法为基础的三维估算系统中动态设置边界品位，从而迅速地动态圈定出经济可供开采的矿块及其品位分布、矿石量、金属量等。可见，本文建立的以克里格法为基础的三维估算系统能够充分满足矿山动态管理的需要。

致谢：江西省地质矿产勘查开发局赣南地质调查大队提供的野外数据对此研究帮助极大，在此深表感谢！

5 参考文献

[1] 王登红，王瑞江，李建康，赵芝，于扬，代晶晶，陈郑辉，李德先，屈文俊，邓茂春，付小方，孙艳，郑国栋.中国三稀矿产资源战略调查研究进展综述[J].中国地质，2013,40(2)：361-370.

[2] 袁忠信，李建康，王登红，郑国栋，娄德波，陈郑辉，赵芝，于扬.中国稀土矿床成矿规律[M].北京：地质出版社，2013：1-116.

[3] DZ/T 0204—2002，稀土矿产地质勘查规范[M].2003.

[4] 王仁铎，胡光道.线性地质统计学[M].北京：地质出版社，1988：4-10.

[5] 陈欢欢，李星，丁文秀.Surfer 8.0等值线绘制中的十二种插值方法[J].工程地球物理学报，2007，4(1)：52-57.

[6] 邓茂春，王登红，曾载淋，张永忠，赵芝，邹新勇，陈斌锋.风化壳离子吸附型稀土圈矿方法评价[J].岩矿测试，2013，32(5)：803 -809.

[7] 陈郑辉，肖克炎，吴健生，朱裕生.阿舍勒铜锌矿床三维模型的研究[J].地质论评，2000,46(增刊): 319-324.

[8] 张新宇，肖克炎，刘光胜，揣媛媛.阿舍勒铜矿可视化储量计算的指示克里格法应用研究[ J].吉林大学学报(地球科学版)，2006，36(1)：305-308.

[9] 杨东来.地质体三维建模方法与技术指南[M].北京：地质出版社，2007：65-69.

[10] 刘海英，刘修国，李超岭.基于地质统计学法的三维储量估算系统研究与应用[J].吉林大学学报(地球科学版)，2009，39(3)：541-546.