普通克里格法在草桃背铀矿床资源量估算中的应用

冯巧云 马 涛 邓 宇

(江西省地质局有色地质大队,江西 赣州 341000)

地质统计学是以区域化变量理论作为理论基础,以变异函数作为主要工具,对既具有随机性又具有结构性的变量进行研究的科学,其核心即克里格法,是一种无偏的最小估计方差的估值方法[1-5]。随着应用领域不断扩大,地质统计学理论更加成熟,各种克里格方法、技术不断出现,如参数地质统计学中的普通克里格法、泛克里格法、析取克里格法,非参数地质统计学的指示克里格法、概率克里格法等,其中普通克里格法是其最基本、应用最为广泛的估值方法[6-10]。

本研究借助MinePlan3D矿业软件(原MineSight矿业软件),以隐式建模方法对赣南草桃背铀矿床进行三维地质建模,在充分研究该矿床地质特征及品位变化规律的基础上,制订地质统计学实施方案,以铀边界品位0.01%圈定矿化域,利用普通克里格法对铀品位进行限制性插值,构建矿床资源模型,实现矿床边际品位动态化选择。根据变异函数拟合结果,揭示区域化变量铀品位在垂深方向上比水平方向上矿化连续性更好,为该区后续深部隐伏找矿工作提供有益参考。

1 矿床地质特征

研究区位于武夷山成矿带与南岭成矿带交汇处,大余-会昌Ⅱ级EW向构造与广昌-周田Ⅱ级NE向构造的复合部位,会昌盆地东缘中基性火山岩体带和富城岩体接合处(图1)。区内出露地层主要为震旦-寒武系变余长石石英砂岩、千枚岩、板岩、碳质和硅质板岩,局部可见混合岩化。上白垩系红色砂砾岩、砂岩、粉砂岩不整合在震旦-寒武系之上,中基性火山岩则产于红层下部。富城岩体是区内主要岩体,横跨赣闽二省,东西长约60 km,南北宽18～22 km,总面积为1 095 km2。北部侵入于震旦、寒武系,南部与混合岩呈侵入接触,东部侵入于中侏罗系漳平统,西部被会昌上白垩系红层所覆盖,主要岩性为中粗粒似斑状黑云母花岗岩、细粒似斑状黑云母花岗岩,后期脉岩广泛发育,常见花岗斑岩、石英斑岩、煌斑岩、石英脉等。构造以断裂构造为主,以EW、NE向和环状构造组成构造骨架,断裂性质以压扭性为主。

图1 研究区区域地质特征Fig.1 Regional geological characteristics of the study area

根据矿体特征、构造特征、矿物共生组合、硫同位素测量、沥青铀矿全分析和绝对年龄测定结果,以及铀矿脉之间的穿插关系,可将草桃背铀矿床分为产于火山机构基底的细粒似斑状黑云母花岗岩中岩浆期后热液型(Ⅰ类型)、产于火山震碎带的火山热液型(Ⅱ类型)、产于火山震碎带之紫色带与绿色带界面附近的淋积叠加型(Ⅲ类型),以及产于火山斜坡的火山热液型(Ⅳ类型)4种铀矿化类型。其中,Ⅱ类型铀矿化及Ⅲ类型铀矿化为本矿床主要矿化类型。矿体规模大小不一,矿体形态多为似层状、透镜状,矿体产状多变,受火山机构控制,主要矿体由40m×25m工程间距所控制,矿体产状总体为EW向展布,但与火山口活动关系密切的Ⅱ类型和Ⅲ类型产状受火山构造控制,Ⅰ类型产状为近EW走向,倾向N,Ⅳ类型则主要受NE向构造控制,各矿化类型倾角变化较大,但倾角均较陡,一般为40°～88°,在走向和倾向上矿化都很不均匀,矿体品位变化较大。

2 三维地质建模与样品统计分析

2.1 数据库建立

本研究数据库建立过程中收集了草桃背铀矿床历年来的勘探资料,包括钻探、坑探、槽探等地质资料。收集整理了547个钻孔工程的7 212个样品数据,按照软件对原始信息数据的要求,共建立了孔口坐标表、钻孔测斜表、钻孔化验表与钻孔岩性表4个数据文件。所有数据均经过细致全面的检查校对,可作为矿床建模的基础数据。

2.2 三维地质模型构建

本次三维地质建模构建了岩体模型及矿化域模型。岩体建模采用了隐示建模方法,其原理是基于隐函数对离散的三维点数据进行插值实现表面重建[3]。其中高斯径向基函数(Radial Basis Function,RBF)是常用的隐函数,最早由HARDY引入[4],国际通用三维矿业软件如Leapfrog、MinePlan3D、Micromine、Surpac等也采用该函数作为隐式建模的插值函数[11-15]。MinePlan3D软件的隐式建模功能使用RBF函数生成表面较光滑自然的实体,提供了地表、风化、侵入、断层、不整合、品位壳、岩脉及地层等预设选项以及相互之间的地质时代关系约束功能。草桃背矿床地质数据丰富,通过隐式建模,高效立体展示了火山机构的分布及各类岩体构造的接触关系,弥补了以往只能用二维图件展示地质体导致的认识不足,且建模效率与质量远高于显式建模,为周边类似复杂地质条件的火山口(小富足、上富等)探索提供了有效工具及参考样板。本研究通过隐式建模构建的安山岩实体模型如图2所示。

矿化域模型是矿床成因及其在三维空间上变化特征的具体反映,通过矿化域模型不仅能在三维视角下展示矿化现象的几何空间形态,同时也为后续品位模型插值提供了限制约束域。本研究采用剖面解译法进行建模,在各个剖面对矿化边界进行界定,按照几何法的圈连和外推原则,以铀边界品位0.01%在剖面上圈定低品位矿化域,根据矿化域剖面在水平面的投影分布圈定各水平分段。

2.3 样品统计特征分析

2.3.1 原始样品分析统计

利用地质统计学方法进行矿床资源储量评估时,有必要对样品数据进行统计特征分析。只有在样品数据的统计特征确定以后,方可有效选择地质统计学插值方法[2]。根据建立的数据库,对原始样品铀品位进行了统计分析,结果如图3所示。通过图3可以看出原始样铀品位服从偏正态分布,故本研究选用普通克里格法进行插值计算。

图3 含特高品位的原始样的铀品位直方图Fig.3 Uranium grade histogram of the original sample containing extremely high grade

2.3.2 特高品位处理

对于应用普通克里格法插值的原始数据需要对特高品位进行处理,以保证参数估值精度。目前,常用的特高品位判断方法有置信区间判断法、图形判别法、公式判别法等[2],本研究采用图形判别法确定特高品位。从原始样铀品位累计概率分布曲线图(图4)可以看出,在99.9%分位数上,其直线分布出现明显的转折拐点,根据统计学方法对特高品位的判断与替换方法,该99.9%分位数上对应的品位值1.66%即为特高品位临界值,故以1.66%品位值作为上限值代替原始样铀品位中的特高品位。

图4 原始样铀品位累计概率分布曲线Fig.4 Cumulative probability distribution curve of original uranium grade

2.3.3 样品组合与统计分析

地质统计学要求有效数据必须建立在固定长度的支撑上,因此,有必要对样长不等的钻孔、坑探品位数据按一定的长度进行组合计算并进行分析[1]。组合样长度一般根据矿山开采方法、夹石剔除厚度、矿床类型及品位在垂深方向的变异程度等确定;通常要求一个采矿分段(台阶)高度内需有1个或若干个组合样;组合样的统计数字特征需与原始样相近,不取过大的组合样长度是为了降低样品组合导致的品位平均化程度,并确保组合样数量满足地质统计学计算变异函数需有较多的样本个体为基础的要求。本研究建模选取1.25 m样品组合长度进行统计分析。组合样与原始样的统计特征值比较见表1。由表1可知:无论是否处理了特高品位,二者的均值、标准差及变化系数等都相近。组合样品位分布直方图(图5)分析表明:组合样铀品位服从正态分布,且与原始样品位直方图分布近似,因而取1.25 m作为组合样长度较适宜。经过特高品位处理后,组合样品位的变化系数明显降低,说明数据的稳定性得到加强,为后续变异函数分析奠定了基础。

表1 1.25 m组合样铀品位统计及其与原始样的统计比较Table 1 Uranium grade statistics of 1.25 m composite sample and its statistical comparison with original sample

图5 处理特高品位后的组合样铀品位直方图Fig.5 Uranium grade histogram of combined sample after processing ultra-high grade

3 变异函数计算及结构分析

3.1 变异函数模型拟合与结构分析

矿区勘查控制工程基本网度间距为40 m(走向)×25 m(倾向),钻孔垂深方向铀品位变异大,因此在计算变异函数时,基本滞后距取目前最密工程间距的1/2即12.5 m,搜索角为15°,进行0°、15°、30°、…、165°共12个水平方向以及垂直方向(倾角分别设为0°、15°、30°、…、75°)的变异函数计算。

根据铀的试验变异函数等值线图(图6)分析,结合矿体产状,选择0°与90°这一组试验变异函数进行拟合。以90°方向为例,得出变异函数及其拟合结果如图7所示,各方向上的参数取值见表2。

表2 组合样铀品位球状模型变异函数参数Table 2 Parameters of variogram of spherical model of uranium grade for composite samples

图6 铀的试验变异函数等值线(水平方向)Fig.6 Isoline of test variogram of uranium(horizontal direction)

图7 组合样铀品位试验变异函数曲线(水平0°方位、倾角30°)Fig.7 Variable function curve of combined sample uranium grade test (horizontal 0° azimuth,dip 30°)

铀的试验变异函数结果拟合表明:铀品位变异函数具有几何异向性特点,表现为在水平0°/倾角30°、水平90°/倾角30°和垂深方向上具有相同的块金常数值和基台值,但变程不同。据上述铀变异函数几何异向性参数,矿业软件可自动进行套合。

3.2 搜索椭球体参数确定

搜索椭球体用于克里格法插值对邻域样品搜索时,搜索半轴(半径)一般不宜超过变异函数的变程[2]。根据变异函数拟合确定的变程等参数,估值的普通克里格法估计邻域的半径长度分别为长轴方向40 m,短轴方向25 m,垂直方向1.8 m;搜索椭球体全域扇区数8个,每个待估块段参与克里格法估计的最少样品数有2件,最多15件。

3.3 交叉验证

试验得到的变异函数模型有必要进一步进行理论检验。交叉验证结果见表3。由表3可知:实际值与估计值的均值偏差近于0;且残差分布图服从正态分布,符合交叉验证的2个理论判据,表明用该组变异函数参数对组合样铀品位进行估计是无偏的,满足区域化变量的内蕴假设[2],即通过地质统计学方法进行品位插值建模具有合理性。

表3 交叉验证结果Table 3 Cross-validation results

4 块段模型及克里格估值

待估块段尺寸确定需要综合考虑矿床类型、矿体规模、夹石剔除厚度、探矿工程间距等因素。该矿主矿体厚度一般为2～6 m,倾向N,倾角30°左右,地下采矿作业不允许过度贫化,块段尺寸在Z方向不宜超过3 m。由于本矿床勘查以40 m(走向)×25 m(倾向)工程间距求控制资源量,故本研究取该投影网度最小间距的1/2,即10 m×10 m作为待估块段水平X、Y方向的尺寸。因此,模型块段尺寸设置为10 m(东向)×10 m(北向)×2.5 m(垂直方向)。

设置好搜索椭球体参数,并在矿化域边界限制的约束前提下对铀品位进行普通克里格法估值,通过估值得到每个块段的克里格品位及块段的克里格方差值,模型块段估值的三维效果如图8所示。图8显示铀矿化非常明显地环绕火山口分布,表明铀矿化严格受火山机构控制,且铀矿体主要分布在火山口西部。

图8 模型块段估值的三维效果Fig.8 Three-dimensional effect of model segment estimation

5 资源量估算及验证对比

品位插值后,可快速统计出不同边际品位下的资源量。草桃背矿床模型以铀克里格边际品位0.047%估算矿床资源量,与原用地质块段法(边界品位0.03%,最低工业品位0.05%)估算结果进行了对比,发现两者估算的矿石量相对误差为0.34%,金属氧化物量相对误差为3%,平均品位相对误差为3.2%,均在允许误差范围内。草桃背矿床模型资源量估算结果的铀克里格边际品位与矿石量的关系如图9所示。由图9可知:在矿石量a4与a5之间的靠近a4处,品位-吨位曲线出现明显拐点,其对应的边际品位低于0.05%,约0.047%。

图9 草桃背矿床铀矿块边界品位与矿石量关系Fig.9 Relationship between uranium block boundary grade and ore quantity in Caotaobei Deposit

对本次地质统计学法资源量估算结果进行全局验证,首先分别计算统计出间隔2.5 m的每个水平(中段)的组合样品平均值及该水平块段模型克里格法估计平均值,分别绘制连线图;然后在垂直方向上对比相同水平两条连线的吻合程度,结果如图10所示。

图10 草桃背矿床铀品位全局验证曲线Fig.10 Global verification curves of uranium grade in Caotaobei Deposit

由图10可知:组合样铀品位曲线与克里格法铀品位曲线吻合度好,反映出克里格法估计结果可靠性好。

6 结 论

(1)本研究使用普通克里格法估算草桃背铀矿床资源量,并与原用地质块段法的估算结果进行了对比,两种方法估算的矿石量相对误差为0.34%,金属氧化物量相对误差为3%,平均品位相对误差为3.2%,根据品位-吨位关系图,提出了用铀克里格品位0.047%作为边际品位的单指标方案,为矿山今后适应市场经济需要调整边际品位奠定了基础。

(2)采用隐式建模方法构建了该矿床岩体三维地质模型,弥补了以往仅有二维认识的不足,直观地表征了火山机构内各岩石的分布形态,对红盆底部的构造变异部位和紫色蚀变带分布有了更清晰的了解,为进一步探索红盆及开展火山岩型铀矿深部找矿工作提供了地质依据。

(3)根据试验变异函数拟合得到的铀品位球状模型理论变异函数,其水平方向主轴变程为63 m,次轴变程为56 m,钻孔方向变程为75 m,即垂深方向变程比水平方向大,说明垂深方向上的矿化比水平方向上矿化连续性更好。