研究区为本溪市歪头山尾矿库, 东西长1.87 km, 南北宽0.85 km, 总汇水面积约1.103 km²(见图1)。 位于辽宁省本溪市溪湖区歪头山镇境内, 距离歪头山选矿厂以南2.5 km。 尾矿库原地形横断面呈“ V” 字型, 是一个低山丘陵区的谷地, 在尾矿库边缘每隔5 m左右分布排尾管道, 矿浆经过四周尾矿库边缘流入中心位置, 然后从中心位置向西南方向流入过滤塔。 采样实验选择高分五号卫星过境且天气晴朗的日期进行, 为使采样样品具有代表性, 于2019年6月25日沿排放口及东排放口均匀采集样本77个, 每个采样点利用手持GPS机记录地理坐标, 采集到的样本放入塑封袋密封现场称重带回实验室。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区概况图Fig.1 Overview map of the study area

对原始尾砂样品进行X荧光光谱测试化验成分, 并针对SiO₂、 全铁含量(TFe)及其他主要成分含量进行绘图分析, 结果如图2。 选矿厂排放的尾砂的主要化学元素为SiO₂, Fe, 其次为CaO, MgO和Al₂O₃等, Cu, Pb和Zn等的有色金属含量非常低, 尾砂主要的金属矿物为磁铁矿、 赤铁矿、 其次为褐铁矿、 黄铁矿等, 微量的雌黄铁矿、 毒砂等, 铁是尾砂中含量最多的金属元素。 尾砂的粒度相对较粗, 大部分在74 μ m以上, 小于37 μ m粒径仅占16.84%。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 尾砂样本化学成分含量图Fig.2 Chemical composition of tailings samples

将采集到的样品烘干处理后称重, 对比现场采样称重结果计算含水重量, 将含水重量与样本重量比作为样本实测含水率。 烘干后的样本缓慢注水至饱和, 通过自然蒸发控制含水率, 每隔6 h进行称重计算尾砂含水量, 直至含水率达到稳定, 即风干状态。 室内光谱实验采用美国SVC-HR-1024便携式光谱仪进行测试, 光谱仪波段范围为350~2 500 nm, 共973个波段。 测试在密闭黑暗环境中进行, 将装有样品的黑色容器放在载物台上, 使用功率为1 000 W的卤素灯模拟太阳光作为光源, 卤素灯距样本表面70 cm, 照射方向与水平面的夹角为60° ; 光谱仪镜头垂直距离样本47 cm, 镜头视场角为4° ; 每个含水尾砂样品旋转测试3次光谱后取平均值。 最终得到不同含水率尾砂光谱曲线312条。

高分五号卫星是世界首颗实现对大气和陆地综合观测的全谱段高光谱卫星, 于2019年3月21日正式投入使用。 在可见光及短波红外波段可以获取330个波段的光谱数据。 其中, 可见光近红外波段部分覆盖390~1 090 nm, 光谱分辨率优于5 nm; 短波红外波段光谱覆盖1 000~2 510 nm, 分辨率优于10 nm; 影像空间分辨率为30 m, 幅宽60 km。

本工作使用的数据为与采样时间同期的2019年6月25日的高分五号影像, 预处理在软件ENVI 5.3中进行。 对研究区影像进行辐射定标, 并采用FLAASH模型对图像进行大气校正, 并剔除无反射率信息的水汽吸收波段B193— B200及B246— B262。

高光谱数据波段数据量较大, 并且包含大量冗余信息与背景信号。 竞争性自适应加权重采样法(competitive adaptive reweighted sampling, CARS)以自然界适者生存法则为基础, 能够有效筛选出最优光谱敏感波段。 CARS是一种基于蒙特卡洛采样以及偏最小二乘(PLS)回归系数的波段优选方法, 首先通过蒙特卡洛采样选择部分样本建立PLS模型, 计算各波段反射率的绝对值权重, 通过指数递减函数去除回归系数较小的波段变量; 并对选定的光谱波段子集进行交叉验证, 获得训练集均方根误差最小的模型所对应子集即筛选出的敏感波段集合^[13]。

假设所有样本的波段反射率数据矩阵为X_{m× n}, y_{m× 1}为所有样本目标性质矩阵; W为数据矩阵的组合系数, T为X的得分矩阵, 是X与W的线性组合; c为y与T建立PLS模型的回归系数; e为误差向量, 则有

T=XW(1)

y=Tc+e=XWc+e=Xb+e(2)

式中, b为组合系数与回归系数的结合, 权重ω 定义为

${{\omega }_{i}}=\frac{\left| {{b}_{i}} \right|}{\sum\nolimits_{i=1}^{n}{\left| {{b}_{i}} \right|}},\text{ }\!\!~\!\!\text{ }i=1,\text{ }\!\!~\!\!\text{ }2,\text{ }\!\!~\!\!\text{ }\ldots ,\text{ }\!\!~\!\!\text{ }n$(3)

作为变量优选的评价标准, |b_i|与ω 越大, 该波段变量越重要^[14]。 根据CARS筛选出少量敏感波段, 可以有效降低高光谱数据维数, 达到提高高光谱数据建模效率的目的, 并且便于后续光谱三维变量组合工作的进行。

利用Savizky-Golay对室内可见光-近红外光谱进行平滑去噪处理, 消除光谱信号中的随机噪声, 提高样本信号的信噪比^[15]。 并利用三次样条插值法对室内可见光-近红外光谱进行重采样, 使其波段与高分五号卫星传感器波段相对应。

通过波段组合运算形成包含多波段信息的新的光谱指数参量, 可以消除背景噪声、 以及大气等因素的影响, 指示意义、 提高模型反演地物参数的精度。 三维波段光谱指数法(three-band spectral index, TBI)通过增加维数拓展了传统的二维光谱指数方法, 为挖掘光谱与尾砂属性间更深层的光谱响应关系提供了可能^[16]。 利用尾矿光谱建立五种TBI指数, 如式(4)— 式(8)

TBI1=Ri/(Rj+Rk)(4)

TBI2=(Ri-Rj)/(Rj-Rk)(5)

TBI3=(Ri+Rj)/Rk(6)

TBI4=(Ri-Rj)/[(Ri-Rj)-(Rj-Rk)](7)

TBI5=(Ri-Rj)-(Rj-Rk)(8)

式中, R_i为第i波段反射率值; R_j为第j波段反射率值; R_k为第k波段反射率值。

随机森林(random forests, RF)是近年来广泛应用于高光谱建模反演的经典机器学习方法, 具有参数设置简单、 能够良好挖掘变量间非线性关系的优势。 随机森林由Breiman等于2001年提出, 是一种基于决策树的集成分类器^[17]。 随机森林的决策树构建基于Bootstrap抽样法, 即从原始数据集中有放回地多次随机抽取样本作为子数据集。 在每个决策树的内部节点随机抽取若干个变量作为候补的分枝变量。 最终所有决策树输出结果的均值作为每个样本的属性预测值^[18]。 由于建模数据包含光谱反射率及TBI指数两类数据, 在RF建模过程中能够在每棵决策树节点对输入变量重要性进行进一步评估与选取, 并结合大量非线性决策树, 获得稳健的尾砂含水率回归预测结果。 本研究中随机森林模型采用500的决策树数量, 随机森林回归预测工作通过MATLAB的TreeBagger分类器实现。

极限学习机(extreme learning machine, ELM)是一种基于单隐含层前馈神经网络(single hidden layer feedforward neural network, SLFN)的机器学习方法。 与传统的SLFN不同, ELM通过随机生成或人为给定输入层权重, 且不需要更新, 具有训练参数少、 泛化能力强的特点, 且能够克服训练速度慢以及局部最优解的问题^{[19, 20]}。 对基于影像的大范围尾砂含水率预测, 能够得到良好应用。 但由于随机给定权值与偏差, ELM通常需要较多隐含层节点才能获得理想精度。 粒子群优化算法(particle swarm optimization, PSO)是一种进化计算技术, 源于对鸟群捕食的行为研究^[21]。 粒子群算法通过无质量的粒子模拟篷群中的鸟, 每个粒子具有两个属性: 位置与速度, 速度代表粒子的移动快慢, 位置决定粒子移动方向。 每个粒子在搜索空间中单独搜索最优解, 并且对比整个粒子群的个体极值获得当前全局最优解。 其余粒子们追随当前的最优粒子在解空间中搜索(即迭代过程), 粒子通过自己找到的个体极值与当前全局最优解调整自己的速度与位置, 直至全局最优解满足结束条件^[22]。 相比于ELM算法, PSO-ELM算法能够依靠较少的隐含层节点获得更高的精度。 本工作所建立PSO-ELM模型参数如图3, v_max为最大粒子速度: v_min为最小粒子速度; c₁为个体学习因子加速度常数; c₂为社会学习因子加速度常数; M_num为最大迭代次数。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 建模方法及参数设置Fig.3 Modeling method and parameter setting

深度学习是一种数据驱动的方法, 能够通过多层神经网络来提取具有多个抽象级别的复杂特征。 其中, 卷积神经网络(convolutional neural network, CNN)作为深度学习的代表算法之一, 近年来在光谱分析中的应用逐渐增多^[23]。 卷积神经网络是一种包含卷积层的深度神经网络; 卷积层的作用是自动提取局部特征, 激活函数层通过为神经网络添加非线性映射, 从而使深度神经网络具备分层的非线性映射学习能力。 通过多层卷积层, CNN算法能够不断放大输入变量中的有效特征、 抑制无用特征对于建模工作的干扰, 并且在建模数据数量多时能够取得更好的建模效果, 适用于类似本研究建模输入变量复杂、 光谱数据较多的情况。 由于卷积层输出的映射特征量往往非常大, 通常需要添加池化层来降低输出向量维数。 全连接层则负责建立提取出的特征与输出目标间的映射^[24]。 具体网络结构与参数如图3, f_num为卷积核数量, k_size为卷积核大小, u_num为神经元个数, 各层神经网络所使用激活函数均为ReLU。

利用以上三种反演方法对尾砂光谱及TBI指数数据集进行建模。 实验获得了不同含水率尾砂光谱曲线312条, 其中208条作为训练集样本, 104条作为验证集样本。 选取决定系数R²、 均方根误差(root mean squared error, RMSE)、 相对分析误差(residual prediction deviation, RPD)作为模型精度评价指标。 R²越高、 RMSE越低, 则证明模型反演精度越高。 RPD的定义为预测样本标准差与均方根误差的比值, 当RPD≥ 2时表明模型具有极好的预测能力, 当1.4≤ RPD< 2时表明模型可对样品作粗略估测, 而RPD< 1.4则表明模型无法对样品进行预测^[25]。

不同含水率尾砂的光谱如图4, 反射率随着尾砂样品含水率的减小, 对应光谱曲线呈现规律性的变化, 光谱反射率随含水率的减少而上升; 尾砂光谱在670、 912、 1 451、 1 921和2 350 nm处存在吸收谷, 其中912、 1 451和1 921 nm处吸收谷较为明显, 912 nm处吸收谷为Fe³⁺的宽吸收带, 其吸收深度不随含水率的减少而变化, 1455nm和1930nm处为O— H吸收谷, 吸收深度大小随着含水率的减少而减小^{[26, 27]}。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同含水率尾砂光谱Fig.4 Spectra of tailings with different moisture contents

光谱反射率随含水率的增加整体呈下降趋势, 且越向近红外波段发展, 光谱反射率随含水率的变化越为明显, 不同含水率光谱的分离性更好。 当含水率趋于饱和时(接近22.19%), 尾砂光谱反射率随含水率的降低速度逐渐变慢, 直至达到最低值。 当含水率小于15.07%时, 2 100~2 500 nm波段范围内光谱曲线较为平坦, 反射率值高于1 100~1 850 nm波段区域反射率; 当含水率大于15.07%时, 2 100~2 500 nm波段光谱曲线出现明显的向下弯曲形态, 反射率相较1 100~1 850 nm波段更低。

CARS方法对变量的特征优选基于达尔文进化论的“ 适者生存” 理论。 将通过蒙特卡洛采样筛选出的多个数据集模型中交叉验证标准误差最小的模型对应波长变量作为优选后的敏感波段^[28]。 利用CARS方法对尾砂的含水率与光谱进行敏感波段筛选, 共筛选出含水率敏感波段18个, 波段分布位置如图 5(a)。 对筛选出的尾砂水分敏感波段与尾砂含水率进行相关性分析, 结果如图5(b)。 尾砂的含水率敏感波段主要分布于三个区域, 可见光波段的400~500 nm密集分布5个敏感波段, 相关系数达-0.53~-0.55; 700~1 290 nm范围内共7个敏感波段, 分布于912 nm处的Fe³⁺宽吸收带两端, 相关系数达-0.57~-0.62; 2 000 nm后的近红外波段共有6个敏感波段, 分布于1 930 nm的O— H吸收带右肩周围, 该部分波段与尾砂含水率相关性表现最佳, 达-0.71~-0.75, 与尾砂含水率呈现强相关关系。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 尾砂水分敏感波段CARS筛选结果 (a): 敏感波段位置; (b): 敏感波段相关性分析Fig.5 The tailings moisture sensitive bands screened by CARS (a): Sensitive bands; (b): Correlation analysis of sensitive bands

利用筛选后的18个敏感波段的反射率逐个组合进行TBI指数计算, 探索不同波段组合方法与含水率间的相关性。 并引入皮尔森相关系数衡量TBI指数与尾砂含水率间的相关程度。 最终得到的各波段组合相关系数如图6所示。 图6展示了不同角度的TBI指数相关性分布, 以及最佳相关性系数所在位置切片, 每个元素的x, y, z坐标表示用于计算该指数的三个波段, 将TBI指数相关性矩阵进行旋转展示, 矩阵在水平及垂直方向的旋转角度标注于矩阵左上角。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 基于CARS筛选后敏感波段的TBI指数相关性分布Fig.6 Three-dimensional spectral index correlation distribution based on sensitive bands screened by CARS

每种TBI指数相关性最高值及波段组合方式列于表1。 可见5种TBI指数的相关性最大值分布于0.74~0.84间, 均与尾砂含水率间达到强相关或极强相关^[29]。 如图6, 由三个波段差值组成的TBI5指数表现最佳, 图上呈红色的指数相关系数达到0.60以上, 与含水率呈强相关关系, 达到红色的TBI5波段组合指数共有1 527个。 并且TBI5指数与尾砂含水率的相关系数最高达到0.84, 组成该TBI指数的三个波段分别位于反射率最低的可见光波段437.39 nm处、 铁离子吸收带附近的1 097.47 nm波段以及处于O— H吸收带的1 990.67 nm波段。 对比敏感波段相关性分析结果, TBI指数显著地提高了光谱变量与尾砂含水率间的相关性, 通过三个波段相减计算得到的TBI指数, 能够放大尾砂光谱中水含量相关波段的反射率变化, 更深层地挖掘光谱与尾砂含水率的响应关系, 从而与含水率间达到较高的相关性系数。

表1

Table 1

表1(Table 1)

表1 最佳三维指数波段组合 Table 1 Band combination of the best three-dimensionam indices 指数 名称 波段组合 最佳相关 系数 TBI 1 1097.47 nm/(1 990.67 nm-1 097.47 nm) 0.76 TBI2 (1 291.45 nm-2 041.18 nm)/(2 041.18 nm-1 990.67 nm) 0.81 TBI3 (2 285.71 nm+1 015.88 nm)/ 1 990.67 nm 0.74 TBI4 (1 990.67 nm-2 176.16 nm)/ [(1 990.67 nm-2 176.16 nm)- (2 176.16 nm-963.55 nm)] 0.82 TBI5 (1 097.47 nm-1 990.67 nm)- (1 990.67 nm-437.39 nm) 0.84

表1 最佳三维指数波段组合 Table 1 Band combination of the best three-dimensionam indices

根据前文的相关性分析结果, 计算所有样本的5种最佳TBI指数, 并结合CARS筛选后的18个含水率敏感波段建立新的联合数据集(共23维)。 分别利用RF, PSO-ELM与CNN三种机器学习方法, 对比18个波段原始反射率数据集与TBI指数联合数据集的含水率反演精度, 结果如表2。 原始反射率数据集与PSO-ELM方法联合建立的含水率模型验证精度相对较高, 达到R²=0.84, RMSE=2.75%, RF方法和CNN方法分别达到了R²=0.78, RMSE=3.18%和R²=0.82, RMSE=3.00%。 在三种建模方法下, 在原始的敏感波段数据集中加入了TBI指数后均提高了尾砂光谱的含水率预测精度。 RF模型验证精度由R²=0.78提升至R²=0.85; PSO-ELM模型验证精度由R²=0.84提升至R²=0.88; CNN方法的精度提升最为显著, 结合了TBI指数数据集与CNN方法的尾砂含水率反演模型达到了最高精度R²=0.92, RMSE=1.99%, 能够实现尾砂含水率的高精度预测。

表2

Table 2

表2(Table 2)

表2 不同数据集与建模方法精度 Table 2 Accuracies of different datasets and modeling methods 建模方法 数据集 建模精度 验证精度 R2 RMSE/% RPD R2 RMSE/% RPD RF 敏感波段(18个) 0.97 1.37 4.98 0.78 3.18 2.15 敏感波段(18个)+TBI指数(5个) 0.97 1.19 5.72 0.85 2.67 2.56 PSO-ELM 敏感波段(18个) 0.86 2.50 2.53 0.84 2.75 2.50 敏感波段(18个)+TBI指数(5个) 0.88 2.34 2.73 0.88 2.39 2.89 CNN 敏感波段(18个) 0.90 2.20 3.09 0.82 3.00 2.27 敏感波段(18个)+TBI指数(5个) 0.94 1.72 3.94 0.92 1.99 3.43

表2 不同数据集与建模方法精度 Table 2 Accuracies of different datasets and modeling methods

经过对比实验, 基于TBI指数联合数据集的CNN尾砂含水率反演模型取得了最佳精度。 因此, 提取采样同期高分五号遥感影像上的同波段反射率数据并计算得到TBI指数数据, 利用建立好的CNN模型逐像素进行含水率反演, 并制作尾矿库含水率分布图, 结果如图7。 尾矿库表面含水率呈现明显的变化分布, 其中含水率在5%以下的干砂主要集中在尾矿中间区域, 远离尾矿边缘排放口及中间汇水区域, 约占总面积的34.62%; 含水率处于5%~12%的半干尾砂分布最广, 占总面积的45.85%, 在尾矿库中均匀分布; 含水率处于12%~20%的湿润及潮湿尾砂占总面积的18.31%, 离散分布于尾矿库边缘; 而含水率大于20%的潮湿尾砂出现较少, 仅占总面积的1.22%, 出现在尾矿边缘区域及中间汇水区域。 尾砂含水率的空间预测结果与尾矿周围排放口位置呈现一致规律。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 歪头山尾矿库含水率分布反演图Fig.7 Inversion map of moisture content distribution in Waitoushan tailings-dam

从采样点中随机挑选9个对比采样实测含水率及影像CNN模型反演结果进行精度验证, 采样点实测含水率与影像预测含水率间关系如图8所示。 影像的预测含水率与实测值间呈现较为一致的变化关系, 对比实测值与预测值的验证精度R²达到0.79, RMSE为1.02%, RPD为2.20。 证明本研究所建立数据集及模型能够实现较高精度地对尾矿库表层含水率进行反演。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 歪头山尾矿库含水率反演结果实地验证精度Fig.8 Field accuracy verification of moisture content inversion for Waitoushan tailings-dam