风化壳淋积型稀土矿（又称离子型稀土矿或离子吸附型稀土矿）富含中重稀土元素,是我国重要的战略资源^[1-3]。现阶段,主要以硫酸铵、硫酸镁及其相关的混合电解质作为浸矿剂,采用原地浸矿工艺开采离子型稀土矿^[4-5]。浸矿剂用量是开采风化壳淋积型稀土矿的一个重要参数,若浸矿剂用量过多,增加开采成本、加重环境污染风险,若浸矿剂用量过少,浸矿不充分,造成资源浪费^[6-8]。离子型稀土矿原地浸出开采技术规范（报批稿）建议硫酸铵溶液的消耗体积与矿石体积比值为0.33∶1,秦磊等^[9]建议采用稀土储量来确定浸矿剂用量。因此,矿石量或稀土储量是稀土开采的关键参数。采用赣南钻进行生产勘探是估算稀土储量和矿石量的主要方式^[10],一般地,勘探网度越密,矿体稀土储量和矿石量估算越准确,但是过小勘探间距,使得勘探孔数量过多,勘探工作量大。同时,勘探孔是潜在的优先流通道,过多的勘探孔增加注液时浸矿剂溶液串孔的风险,增加注液盲区,导致浸矿不均匀,造成资源浪费^[11]。因此,合理布置勘探网度,在保证勘探精度的基础上减少勘探孔数量,对于风化壳淋积型稀土矿的开采具有重要意义。

针对勘探网度优化方法,国内外学者开展了大量研究。梅晓仁等^[12]基于普通克里金法,提出了露天煤矿勘探网度的优化方法。HOSSEIN等^[13]通过估计品位、矿石厚度和估计误差,以确定勘探钻孔的最佳位置。ANAM等^[14]使用考虑地层不确定性的预测优化模型来确定最佳钻孔位置,并通过假设检验验证了预测模型的拟合优度。文献[15-17]还将智能算法引入,实现勘探网度的优化。上述对勘探网度优化重在优化勘探孔数量,确定数量后,在开采区域均匀布置勘探孔。风化壳淋积型稀土矿稀土品位具有明显的空间变异性^[18-19],均匀布置勘探孔对于稀土品位变化较大的区域存在勘探孔数量不足的问题,导致该区域的储量等参数估算不准确,对于稀土品位变化较小的区域,勘探孔数量布置过多,增加了不必要的勘探工作量。AHSAN等^[20]引入了一种自适应取样方法,以期获得基础地质的精确局部模型,同时最大限度地减少钻孔的数量。但该方法假设数据是无测试误差的,可是实际应用中无法避免测试误差。

为了结合风化壳淋积型稀土矿稀土品位空间变异程度,自适应布置勘探孔,本文假设矿体的单位面积稀土质量分布与地形高程的变化剧烈程度存在正相关关系,以矿体地形高程的特征点作为勘探孔的布置位置,再结合勘探数据,对局部区域进行加密勘探孔,建立风化壳淋积型稀土矿勘探孔优化的特征点法,通过实际矿体的勘探数据讨论方法的合理性。

1   特征点提取方法

矿体的单位面积稀土质量分布与地形高程的变化剧烈程度存在正相关关系,即可以用地形的特征点近似作为单位面积稀土质量分布的特征点。设置勘探线间距和角度（正东依逆时针方向至勘探线转动的角度为α）,在开采边界内布置平行勘探线（图1）,采用张鸿刚等^[21]提出的曲线特征点提取方法,提取每条勘探线上地形高程的特征点,以地形高程特征点作为稀土资源预勘探孔的孔口位置,根据勘探间距的限值对预勘探孔进行修正。计算每条勘探线上勘探孔孔口的间距,取所有勘探间距的均值或中位数作为勘探线的间距,重新确定稀土资源的勘探孔孔口坐标。勘探后,确定单位面积稀土质量分布图,对单位面积稀土质量分布变化大而勘探孔较少的区域进行加密,对加密位置进行再次勘探,综合两次勘探结果,采用储量估算方法确定开采矿体的稀土储量。具体的步骤如下（筛选勘探点的流程如图2所示）：

图  1  开采矿体地形及边界

Figure  1.  Topography and boundary of the ore body to be mined

下载: 全尺寸图片 幻灯片

图  2  根据地形筛选勘探孔孔口坐标的流程示意

Figure  2.  Flowchart for screening the coordinates of exploration hole openings based on topography

下载: 全尺寸图片 幻灯片

步骤1,给定勘探线间距（D）和角度（α）,在边界内布置勘探线。

步骤2,基于地形数据,采用泛克里金法^[22] （漂移函数为二次函数,变异函数为球状函数）确定各勘探线上对应的地形高程。；

步骤3,采用曲线特征点提取方法^[21],提取各勘探线的特征点作为勘探点,第i条勘探线上,各勘探点依次记为p_i1, p_i2, p_i3, …,勘探孔的孔口坐标分别记为（x_i1, y_i1, z_i1）, （x_i2, y_i2, z_i2）, （x_i3, y_i3, z_i3）, …。

步骤4,一条勘探线上,相邻勘探点间距（d）需满足最大间距（d_max）和最小间距（d_min）的要求^[23],本文d_min和d_max分别取15.00 m和50.00 m。对于第i条勘探线,令p_temp = p_i1,计算p_temp和p_i2的间距d_i,2,若d_i,2<d_min,则删除p_i2；若d_i,2>d_min,保留p_i2,令p_temp = p_i2；若d_i,2>d_max,在p_i1和p_i2中间增加一个勘探孔,勘探孔孔口的水平和垂直坐标分别为（x_i1+ x_i2）/2和（y_i1+ y_i2）/2,通过插值算法确定孔口的高程,令p_temp = p_i2；继续计算p_temp和p_i3的间距d_i,3,判断d_i,3与d_min的大小；如此重复,筛选出满足最小间距的勘探点。

步骤5,计算每条勘探线上相邻勘探孔孔口的间距,取所有间距的平均值d_mean与中位数d_median较小值（式（1））作为勘探线间距,重复步骤1—步骤4,确定稀土资源勘探孔孔口坐标。

式（1）中：min( ∙ )为取最小值运算。

步骤6,按照上述设计结果进行勘探,根据勘探孔的品位信息,采用泛克里金法确定单位面积稀土质量分布图,计算上述勘探线位置的单位面积稀土质量,按照步骤3和步骤4确定单位面积稀土质量分布的特征点,计算同一条勘探线上,每个特征点与该线上所有勘探点的间距,若均大于d_min,则保留特征点,否则删除该特征点,将所有保留的特征点作为加密勘探孔,进行再次勘探。

步骤7,综合2次勘探结果,采用地质块段法^[24]确定矿体的稀土储量。

2   稀土配分测试试验

试验场地位于福建省某风化壳淋积型稀土矿区,采用赣南钻^[25]对目标矿体进行勘探,勘探网度为（15~20） m×（15~20） m,约每1 m取一个样品,测试样品的稀土品位,测量取样段起点和终点与孔口的距离。当遇到半风化层或者含水率过高无法取矿样时,结束该孔的勘探。4个矿体分别记为KT1、KT2、KT3和KT4,勘探孔数量分别为175、188、122、223个,采用GPS-RTK测量勘探孔孔口的坐标。矿体的地形、勘探边界和勘探孔位置如图3所示。采用ICP-MS测试4个矿体的离子相稀土配分。

图  3  矿体地形、勘探边界以及勘探孔的位置：（a） KT1；（b） KT2；（c） KT3；（d） KT4

Figure  3.  Topography, exploration boundaries and location of exploration holes： （a） KT1； （b） KT2； （c） KT3； （d） KT4

下载: 全尺寸图片 幻灯片

3   结果与讨论

3.1   稀土配分

4个矿体的稀土配分如表1所列,轻稀土（La-Eu）占比分别为50.1%、49.2%、49.6%和49.1%,重稀土（Ga-Lu+Y）占比分别为49.9%、50.8%、50.4%和50.9%,钇含量均大于25%,铕含量均在0.5%~1.0%范围内,结合《风化壳淋积型稀土矿化工冶金》^[26]和《风化壳离子吸附型稀土矿产地质勘查规范》（2019）^[25]对风化壳淋积型稀土矿的分类原则,4个矿体均可视为富钇富铕重稀土矿。对于重稀土矿,边界品位和最低工业品位分别取0.020%和0.035%,最小可采厚度为1.00~2.00 m,夹石剔除厚度为2.00~4.00 m^[25],本文最小可采厚度和夹石剔除厚度分别取1.00 m和2.00 m。

表  1  矿体的稀土配分

Table  1.  Rare earth partitioning of ore body

矿体	La2O3	CeO2	Pr6O11	Nd2O3	Sm2O3	Eu2O3	Gd2O3	Tb4O7	Dy2O3	Ho2O3	Er2O3	Tm2O3	Yb2O3	Lu2O3	Y2O3
KT1	17.50	4.15	5.29	18.51	3.89	0.72	3.92	0.97	3.90	0.93	2.60	0.22	2.03	0.05	35.32
KT2	20.09	2.38	4.29	17.61	4.03	0.77	4.96	0.67	5.29	0.91	2.62	0.33	2.31	0.33	33.42
KT3	18.61	5.88	4.56	15.87	4.08	0.59	4.34	0.76	5.00	0.93	2.78	0.38	2.33	0.34	33.53
KT4	17.74	3.09	5.44	17.72	4.28	0.80	4.78	0.81	4.67	0.95	3.08	0.40	2.18	0.30	33.76

下载: 导出CSV 

| 显示表格

3.2   稀土品位

4个矿体的取样点数分别为2644、3746、764、3660个,每个取样点品位数据如图4所示。根据勘探信息,可以得到4个矿体的勘探深度范围分别为4.00~33.00 m、6.00~35.10 m、1.20~18.20 m和4.55~38.30 m,平均勘探深度分别为17.17 、21.79 、7.25 m和18.73 m,见矿深度范围分别为1.00~13.00 m、1.00~13.00 m、1.00~12.00 m和1.00~16.00 m,平均见矿深度分别为4.61 m、4.69 m、3.39 m和4.45 m,4个矿体的见矿深度均较浅。采用地质块段法^[24]计算矿体的稀土储量,考虑“穿鞋戴帽”,“穿鞋”和“戴帽”的最大厚度均取4.00 m,得到稀土储量（称为勘探储量）分别为181.43、278.90、62.69、204.22 t。

图  4  矿体的品位信息：(a) KT1；(b) KT2；(c) KT3；(d) KT4

Figure  4.  Ore body grade information: (a) KT1; (b) KT2; (c) KT3; (d) KT4

下载: 全尺寸图片 幻灯片

3.3   勘探孔网优化

当勘探线与正东的方向夹角（α）为0°时,采用特征点法优化KT1矿体的勘探孔网,得到设计勘探孔的结果如图5所示。由于设计勘探孔的孔口坐标与实际勘探孔的孔口坐标存在差异,为了计算设计勘探孔对应的资源量,取与设计勘探孔距离最近的实际勘探孔作为设计勘探孔,采用地质块段法计算设计勘探孔对应的稀土储量（称为设计储量）。勘探孔数量和设计储量随α的变化情况如图6所示。由图6(a)可知,随着α的增加,4个矿体的设计勘探孔数量呈现先减小后增加的趋势,相较于实际勘探孔,4个矿体设计勘探孔数的压缩率（压缩率为勘探孔设计值与实际值的比值）变化范围分别为[38.86%, 66.86%]、[47.87%, 69.68%]、[49.18%, 73.77%]和[47.53%, 82.96%],在α = 90°附近,4个矿体的设计勘探孔数量和压缩率均取得最小值。

图  5  KT1矿体勘探线方向角为0°时勘探孔网优化结果

Figure  5.  Optimization results of KT1 orebody when the α = 0°

下载: 全尺寸图片 幻灯片

图  6  勘探孔数量和设计储量随勘探线方向夹角α的变化情况：(a) 勘探孔数量；(b) 稀土储量

Figure  6.  Number of exploration holes and designed reserves change with direction angle α of exploration line: (a) exploration holes; (b) rare earth reserves

下载: 全尺寸图片 幻灯片

随着α的增加,设计储量的计算结果相对比较随机,在一定范围内波动。定义设计储量的计算误差为设计储量与勘探储量的相对误差,如式(2)所示^[27-28]。4个矿体设计储量的计算误差范围分别为[0.39%, 32.68%]、[0.08%, 17.91%]、[0.16%, 27.20%]和[0.14%, 21.84%],若选择合理的α,设计储量的计算误差可以小于0.20%,若α选择不合理,设计储量的计算误差超过15.00%,由此可见,合理选择α对于提高特征点法设计勘探孔网的精度十分重要。

式（2）中：η为设计储量的误差；M_design和M_explorate分别为设计储量和勘探储量。

3.4   α角的选定

为了合理确定α,以不同α对应的勘探孔数量为横坐标、以设计储量与勘探储量的相对误差为纵坐标绘制数据图,结果如图7所示。由图7可知,随着勘探孔数量的增加,相对误差的变化区间越窄,即设计储量的结果越稳定。对于本文4个矿体,当α分别取170°、5°、178°和2°时,勘探孔数量达到最大值,分别为117、131、90、185个,设计储量的计算误差分别为6.86%、4.75%、7.84%和7.49%。因此,从结果的稳定性和精度两方面考虑,以勘探孔数量最大值作为目标函数（式（3））^[29],确定相应的α。确定勘探线的方向角后,采用提出的特征点法设计矿体勘探孔的空间位置。

图  7  设计储量与勘探储量的相对误差（η）随勘探孔数量的变化情况：（a） KT1；（b） KT2；（c） KT3；（d） KT4

Figure  7.  Relative error of designed reserves and explored reserves （η） varies with the number of exploration holes: （a） KT1; （b） KT2; （c） KT3; （d） KT4

下载: 全尺寸图片 幻灯片

式（3）中：OF为目标函数；N_α为勘探线与正东的方向夹角为α时,采用特征点法设计得到的勘探孔数量；max（ ∙ ）为取最大值运算。

3.5   优化结果分析

当α分别取170°、5°、178°和2°时,4个矿体勘探孔优化结果的压缩率分别为74.52%、78.92%、73.77%和82.96%,即相比实际勘探孔网,采用提出的优化方法可以减少约20.00%的勘探工作量,稀土储量的计算结果的偏差在8.00%以内。

通过式（4）可以计算出每个实际勘探孔位置的单位面积稀土质量^[30],分别以x和y方向作为x和y坐标,以单位面积稀土质量为z坐标,得到单位面积稀土质量的分布情况（图8）,同样的,以x和y方向作为x和y坐标,以勘探孔孔口高程为z坐标,得到勘探孔孔口坐标的分布情况。采用平均曲率^[31-32]量化单位面积稀土质量和勘探孔孔口高程随水平坐标变化的剧烈程度,以单位面积稀土质量的平均曲率为横坐标、以勘探孔孔口高程的平均曲率为纵坐标绘制数据图,结果如图9所示。采用相关系数（式（5））量化单位面积稀土质量与勘探孔孔口高程的相关性^[33-34],结果分别为0.390、0.403、0.882和0.556；同时采用t检验^[35-36]进行显著性检验,计算结果如表2所列,4个矿体p值均在0.001以下,满足小于0.05的显著性检验要求^[37-38],表明单位面积稀土质量与勘探孔孔口高程的平均曲率呈正相关关系,即基本假设“矿体的单位面积稀土质量分布与地形的高程的变化剧烈程度存在正相关关系”是合理的。

图  8  实际勘探结果的单位面积稀土质量分布云图：（a） KT1；（b） KT2；（c） KT3；（d） KT4

Figure  8.  Cloud map of rare earth mass distribution per unit area of actual exploration results（a） KT1; （b） KT2; （c） KT3; （d） KT4

下载: 全尺寸图片 幻灯片

图  9  单位面积稀土质量与勘探孔孔口高程平均曲率的关系：（a） KT1；（b） KT2；（c） KT3；（d） KT4

Figure  9.  Relationship between rare earth mass per unit area and elevation of exploration holes: （a） KT1; （b） KT2; （c） KT3; （d） KT4

下载: 全尺寸图片 幻灯片

表  2  矿体单位面积稀土质量与勘探孔孔口高程显著性的t检验

Table  2.  A t-test for the significance of rare earth mass per unit area of ore body versus exploration hole aperture elevation

矿体		方差方程的Levene检验		均值方程的t检验
		F	显著性	t	df	显著性（双侧）	均值差值	标准误差值	差值的95%置信区间
									下限	上限
KT1	假设方差相等	5.032	0.026	-4.535	172.000	<0.001	-0.216	0.048	-0.310	-0.122
	假设方差不相等	—	—	-4.535	161.981	<0.001	-0.216	0.048	-0.310	-0.122
KT2	假设方差相等	9.849	0.002	-4.301	186.000	<0.001	-0.212	0.049	-0.309	-0.115
	假设方差不相等	—	—	-4.363	174.768	<0.001	-0.212	0.049	-0.308	-0.116
KT3	假设方差相等	25.360	<0.001	-9.952	120.000	<0.001	-0.757	0.076	-0.908	-0.607
	假设方差不相等	—	—	-10.327	100.073	<0.001	-0.757	0.073	-0.903	-0.612
KT4	假设方差相等	23.272	<0.001	-8.944	215.000	<0.001	-0.202	0.025	-0.246	-0.157
	假设方差不相等	—	—	-9.075	159.521	<0.001	-0.202	0.022	-0.245	-0.158
注：“—”代表无此数据。

下载: 导出CSV 

| 显示表格

式（4）、式（5）中：ω_i为第i个勘探孔位置的单位面积稀土质量,kg/m²；ε_ij和L_ij分别为第i个勘探孔中第j个样品的品位和取样长度,L_ij的单位为m；ρ_b为矿体的密度,本文取ρ_b = 1520 kg/m³；κ_explorate和κ_altitude分别为单位面积稀土质量和勘探孔孔口高程的平均曲率,${\overline{\kappa }}_{\mathrm{a}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{u}\mathrm{d}\mathrm{e}}=\sum {\kappa }_{\mathrm{a}\mathrm{l}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{u}\mathrm{d}\mathrm{e}}/\left(I{J}_i\right)$,${\overline{\kappa }}_{\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{l}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}}=\sum {\kappa }_{\mathrm{e}\mathrm{x}\mathrm{p}\mathrm{l}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}}/\left(I{J}_i\right)$。

为了进一步讨论特征点法局部勘探精度,绘制实际勘探和优化设计的单位面积稀土质量的绝对误差分布云图^[39-40],结果如图10所示。单位面积稀土质量的绝对误差主要分布在-0.7~0.7 kg/m²范围内,局部区域的绝对误差的绝对值超过1.3 kg/m²。

图  10  实际勘探和优化设计的单位面积稀土质量的绝对误差分布云图：（a） KT1；（b） KT2；（c） KT3；（d） KT4

Figure  10.  Cloud plot of absolute error distribution of rare earth mass per unit area for actual exploration and optimized design: （a） KT1; （b） KT2; （c） KT3; （d） KT4

下载: 全尺寸图片 幻灯片

由图11可知,对于A-1区域:5个设计勘探孔的矿层的平均品位在0.036%~0.044%范围内变化,矿层厚度的范围为2.0~12.9 m,单位面积稀土质量在1.42~8.01 kg/m²范围内波动,2个未纳入设计的勘探孔的矿层平均品位、矿层厚度和单位面积稀土质量分别大于0.059%、16.8 m和13.4 kg/m²,矿层平均品位和矿层厚度均明显大于设计勘探孔的数据,说明A-1区域的未设计勘探孔范围存在矿层平均品位和矿层厚度的突变,导致单位面积稀土质量的估算误差偏大。A-2和A-3区域单位面积稀土质量估算误差偏大的原因与A-1的相同。

图  11  A-1、B-2和D-2区域勘探的基本参数：（a）平均品位；（b）矿层厚度；（c）单位面积稀土质量

Figure  11.  Basic parameters for exploration in A-1, B-2 and D-2 areas: （a） average grade; （b） ore layer thickness; （c）mass of rare earths per unit area

下载: 全尺寸图片 幻灯片

对于B-2区域,设计勘探孔和未纳入设计的实际勘探孔的矿层平均品位相近,但是未纳入设计的实际勘探孔的矿层厚度显著大于设计勘探孔,B-1、C-2和D-2区域的情况与之类似。对于这些区域,单位面积稀土质量估算误差偏大的原因是未设计勘探孔范围存在矿层厚度突变的情况。

对于D-1区域,d-198的矿层平均品位为0.088%,约是该区域设计勘探孔最大矿层平均品位的2倍,而矿层厚度与该区域设计勘探孔最大矿层厚度相当,说明该区域单位面积稀土质量估算误差偏大的原因是未设计勘探孔范围存在矿层平均品位突变的情况。C-1区域的情况与之相同。

由此可见,单位面积稀土质量估算误差偏大的主要原因是未设计勘探孔范围存在矿层平均品位和矿层厚度突变。

4   结 论

1）假设矿体的单位面积稀土储量分布与地形高程的变化剧烈程度存在正相关关系,设置勘探线,提取每条勘探线上的地形特征点作为设计勘探孔,结合实际勘探孔,对单位面积稀土质量分布剧烈的区域进行加密勘探,以此建立了风化壳淋积型稀土矿勘探网度优化的特征点法。

2）随着勘探线与正东方向夹角的变化,勘探孔数量呈现先增加后减小的趋势,在α = 90°附近,勘探孔数量取得最小值。采用勘探孔数量取最大值时对应的α进行勘探时,能保证勘探结果的精度和稳定性。

3）单位面积稀土质量分布与相应高程的平均曲率的相关系数大于0.390,表明两者的平均曲率呈正相关关系,即论证了基本假设的合理性。