中国辽宁省鞍山—本溪地区探明铁矿储量超100亿吨, 工业储量超40亿吨, 居中国首位, 且易于开采, 具有极大的工业价值。 以辽宁省鞍山市鞍千矿露天采场为研究区, 在研究区内不同地点均匀采集鞍山式铁矿石, 对采集后的鞍山式铁矿石进行钻孔、 取芯以及切块处理, 制成185件块状铁矿样品, 如图1所示。 通过化学分析法对实验样本进行化学分析, 获取了实验样本的铁品位含量。 化验结果表明, 实验样本铁品位分布在11.4%~66.83%之间, 平均铁品位为29.85%。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区域位置图 (a): 中国地图; (b): 鞍山市地图; (c): 研究区域地图; (d): 部分样本Fig.1 Location maps of the research area (a): Map of China; (b): Map of Anshan city; (c): Map of study area; (d): Some samples

分别使用美国Resonon公司生产的Pika L和Pika NIR-320高光谱成像仪对样品进行光谱测试, 仪器的具体参数如表1所示。

表1

Table 1

表1(Table 1)

表1 成像光谱仪参数 Table 1 Parameters imaging spectrograph 主要参数 PikaL PikaNIR-320 Spectral range/nm 385.75~1 025.75 886.49~1 718.86 Spectral bands 300 168 FOV/(° ) 17.6 21.7 空间通道数 900 320

表1 成像光谱仪参数 Table 1 Parameters imaging spectrograph

实验过程中, 将样品放置于移动平台上, 以4盏卤素灯为实验光源, 以标准白板为定标背景。 实验结束后使用与高光谱成像仪配套的Spectronon Pro软件将仪器采集的原始DN值转化为反射率数据。 因无法获取单独像元对应矿物的品位数据, 所以圈定每个样品在高光谱图像中所对应的全部像元为感兴趣区(region of interest, ROI), 提取ROI内的平均光谱代表该样品的光谱信息。 剔除传感器波段边缘噪声较大的数据后, VIS-NIR光谱数据保留的波段范围为404~959.89 nm, SWIR光谱数据保留的波段范围为964.48~1 667.87 nm。 高光谱图像采集与平均光谱提取流程如图2所示。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 高光谱图像采集与平均光谱提取流程Fig.2 Process of hyperspectral image acquisition and average spectrum extraction

分别以多元散射校正(multiplicative scatter correction, MSC)与标准正态变量变换(standard normal variate transformation, SNV)对原始光谱数据进行光谱变换, MSC可有效降低散射对光谱数据的影响, 修正光谱数据的基线平移和偏移现象; SNV可以消除固体颗粒间大小不均与表面散射引起的误差, 增强光谱数据中与品位相关的光谱信息。

此外, 高光谱数据包含大量的干扰和无信息变量, 现有研究通常对其进行特征提取后, 再利用机器学习的方法建立成分反演模型^{[16, 19]}。 虽然如局部线性嵌入(locally linear embedding, LLE)、 等距特征映射(isometric feature mapping, ISOMAP)等经典流行降维算法在高光谱数据特征提取中具有较好的效果^[19], 但其并不适用于高光谱图像^[20]。 原因在于流行降维算法具有较高的计算复杂性, 尤其在百万级像素的高光谱影像降维应用中, ISOMAP算法中需要求解全局点间最短路径矩阵并对大型稀疏矩阵进行特征分解, LLE算法虽然是基于局部邻近点的线性重构, 但处理百万级像素高光谱图像时仍需要较大内存来存储光谱间的距离矩阵^[21]。 因此本工作通过蒙特卡洛无信息变量消除(Monte Carlo uninformative variable elimination, MCUVE)、 竞争性自适应重加权采样(competitive adapative reweighted sampling, CARS)、 连续投影算法(successive projections algorithm, SPA)三种特征波段提取算法来减少高光谱图像的多重共线性和冗余, 以提高效率和可靠性。

径向基函数(radial basis function, RBF)神经网络是具有单隐含层的前向神经网络, 隐含层的变换函数是对中心点径向对称衰减的非负非线性函数。 一般采用高斯函数作为隐含层激活函数, 该函数是局部响应函数, 具有学习速度快、 避免出现局部极值等优点, 在解决非线性问题中有着广泛的应用^[22]。

极限学习机(extreme learning machine, ELM)是具有单隐含层的前向神经网络, 它不仅对非线性数据具有较强的拟合能力, 且在训练过程中无需调整隐含层和输入层之间的连接权值, 只需设置隐含层神经元个数即可。 采用稳定性较强的Sigmoid函数作为隐含层激活函数, 并设置隐含层神经元个数为30^[23]。

以决定系数(R²)、 相对分析误差(RPD)作为反演模型准确性的评价指标, 以平均绝对误差(MAE)、 均方根误差(RMSE)作为反演结果精度的评价指标。 R²和RPD越高, MAE与RMSE越小, 说明模型的拟合程度越好。

鞍山式铁矿样品的原始光谱反射率曲线与利用MSC和SNV对原始光谱数据变换后的曲线如图3所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 鞍山式铁矿样品原始及变换后光谱曲线 (a): VIS-NIR; (b): SWIR; (c): VIS-NIR-MSC; (d): SWIR-MSC; (e): VIS-NIR-SNV; (f): SWIR-SNVFig.3 Original and transformed spectral curves of Anshan-type iron ore samples (a): VIS-NIR; (b): SWIR; (c): VIS-NIR-MSC; (d): SWIR-MSC; (e): VIS-NIR-SNV; (f): SWIR-SNV

图3(a)与图3(b)为样本在VIS-NIR与SWIR范围内的光谱曲线, 在VIS-NIR范围内光谱曲线整体呈现缓慢降低的趋势, 在550与600 nm处存在两个较为明显的反射峰, 受Fe³⁺电子跃迁的影响, 样本光谱在700~750 nm区间存在一个微弱的宽缓吸收峰。 在SWIR范围内光谱曲线呈现先上升后下降的趋势, 在1 075与1 175 nm处存在两个较为明显的吸收谷。 由于大气中水汽的影响, 光谱在1 300~1 400 nm区间存在反射谷特征和波段震荡特征。 经MSC与SNV变换后的光谱曲线能有效突出上述光谱特征, 如图3(c)至图3(f)所示。

MCUVE将每个波段视为一个变量, 利用蒙特卡洛抽样法选取不同波段并建立偏最小二乘(partial least squares, PLS)模型, 可以求取每个波段所对应的稳定性, 并以PLS模型的交互验证均方根误差(RMSECV)为评价指标选择特征波段个数。 如图4所示, RMSECV随波段变量数量的增加整体呈现先降低后升高的趋势, 在所选变量数量分别等于18、 15、 17、 5、 18以及16次时, PLS模型达到最低的RMSECV, 即VIS-NIR、 VIS-NIR-MSC、 VIS-NIR-SNV、 SWIR、 SWIR-MSC、 SWIR-SNV六种光谱数据分别取稳定性绝对值最大的18、 15、 17、 5、 18和16个波段为该光谱数据的特征波段。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 RMSECV随选定波段数量的变化(MCUVE) (a): VIS-NIR; (b): VIS-NIR-MSC; (c): VIS-NIR-SNV; (d): SWIR; (e): SWIR-MSC; (f): SWIR-SNVFig.4 RMSECV varies with the number of selected bands(MCUVE) (a): VIS-NIR; (b): VIS-NIR-MSC; (c): VIS-NIR-SNV; (d): SWIR; (e): SWIR-MSC; (f): SWIR-SNV

CARS是一种逐步去除无关变量以保留关键变量的特征波段选择算法, 通过蒙特卡洛抽样方法选取N个波段子集, 然后利用这些波段建立PLSR模型, 以回归系数的绝对值作为每个波段的评价指标, 随着采样次数的增加, 保留波段的个数也随之减少, 在达到最低RRMSECV时仍保留的波段视为经算法选取的特征波段。 图5显示了CARS中RMSECV值随采样次数增加而变化的情况, 由图可知, RMSECV值随着采样次数的增加先呈现缓慢降低的趋势, 这是由于算法消除了大量冗余变量而对模型进行了优化, 而后RMSECV值呈现阶梯式升高的趋势, 是由于一些关键变量被消除而导致模型性能的恶化, 在采样次数分别为27、 26、 26、 26、 22和9次时模型达到最低的RMSECV值, 此时上述六种光谱数据仍保留的特征波段个数为20、 22、 22、 9、 10和52个特征波段。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 RMSECV随采样次数增加而变化的情况(CARS) (a): VIS-NIR; (b): VIS-NIR-MSC; (c): VIS-NIR-SNV; (d): SWIR; (e): SWIR-MSC; (f): SWIR-SNVFig.5 Changes of RMSECV with the increase of sampling runs(CARS) (a): VIS-NIR; (b): VIS-NIR-MSC; (c): VIS-NIR-SNV; (d): SWIR; (e): SWIR-MSC; (f): SWIR-SNV

SPA算法不以最小的RMSECV值为收敛条件, 而是以某点的RMSECV值与最小RMSECV值符合F检验(α =0.25)为截止程序, 在一定程度上避免了过拟合现象。 图6显示了SPA中RMSECV值随所选特征波段数量增加而变化的趋势, 在达到图中红点所示位置后, RMSECV曲线仍呈现下降的趋势, 但随着所选特征博段数量的进一步增加, 改善变得微乎其微, RMSECV曲线变得较为平稳, 且所选点与最低点RMSECV值满足F检验(α =0.25), 因此上述六种光谱数据保留的特征波段数量分别为26、 26、 28、 22、 12和17。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 RMSECV随选定波段数量的变化(SPA) (a): VIS-NIR; (b): VIS-NIR-MSC; (c): VIS-NIR-SNV; (d): SWIR; (e): SWIR-MSC; (f): SWIR-SNVFig.6 RMSECV varies with the number of selected bands(SPA) (a): VIS-NIR; (b): VIS-NIR-MSC; (c): VIS-NIR-SNV; (d): SWIR; (e): SWIR-MSC; (f): SWIR-SNV

VIS-NIR与SWIR光谱数据经不同预处理组合算法提取的特征波段分别如图7与图8所示。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 VIS-NIR光谱数据基于不同预处理方法对选择的特征波段Fig.7 Characteristic bands selected for VIS-NIR spectral data based on different pretreatment methods

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 SWIR光谱数据基于不同预处理方法选择的特征波段Fig.8 Characteristic bands selected for SWIR spectral data based on different pretreatment methods

训练集与测试集按照3:1的比例进行均匀划分, 以139个鞍山式铁矿样本的数据作为模型的输入数据, 使用46个鞍山式铁矿样本的数据检验模型的准确性与精度。 本文以决定系数(R²)、 相对分析误差(RPD)作为反演模型准确性的评价指标, 以平均绝对误差(MAE)、 均方根误差(RMSE)作为反演结果精度的评价指标, 基于RBFNN和ELM建立铁矿品位反演模型准确性与精度的评价指标分别如图9与图10所示。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 基于预处理前后光谱数据的RBFNN预测结果 (a): R2; (b): RPD; (c): RMSE; (d): MAEFig.9 Prediction results of RBFNN based on the spectral data before and after preprocessing (a): R2; (b): RPD; (c): RMSE; (d): MAE

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 预处理前后光谱数据基于ELM的预测结果 (a): R2; (b): RPD; (c): RMSE; (d): MAEFig.10 Prediction results of ELM based on the spectral data before and after preprocessing (a): R2; (b): RPD; (c): RMSE; (d): MAE

由图可知, 原始光谱数据反演精度较低, 因其可能受到实验环境、 矿石表面散射等的影响, 此外, 原始光谱数据含有较多的噪声和冗余信息, 光谱特征不明显。 MSC与SNV两种光谱变换方法可有效突出光谱特征信息, 消除表面散射的影响, 从而提升反演模型的准确性与精度。 MCUVE、 CARS以及SPA通过提取特征波段的方法实现光谱特征提取, 可在有效降低光谱数据冗余的同时较好的保留原始光谱数据特征, 提高光谱质量。 不同的预处理组合算法对反演模型的提升效果不同, 在VIS-NIR范围内MSC-CARS-ELM反演模型效果最优, R²=0.90, RPD=3.02, RMSE=3.27, MAE=2.77。 在SWIR范围内SPA-ELM反演模型效果最优, R²=0.80, RPD=2.23, RMSE=4.42, MAE=3.70。 基于RBFNN与ELM所建模型中皆为VIS-NIR-MSC-CARS组合算法精度最高, 相应的预测值与实测值差异图如图11(a)与图11(b)所示, 其中该组合算法利用RBFNN所建模型R²=0.90, RPD=2.75, RMSE=3.59, MAE=3.04, 模型准确性与ELM所建模型相近, 但反演精度低于ELM所建模型。

图11

Fig.11

	Figure Option ViewDownloadNew Window
	图11 预测值与实测值差异图 (a): VIS-NIR-MSC-MCUVE-RBF; (b): VIS-NIR-MSC-CARS-ELMFig.11 Scatterplots of predicated value vs measured value (a): VIS-NIR-MSC-MCUVE-RBF; (b): VIS-NIR-MSC-CARS-ELM

将感兴趣区中所有像素对应的光谱数据经MSC变换后, 以CARS选择的26个波段(442.73、 467.35、 469.4、 494.15、 496.22、 525.25、 558.63、 585.9、 623.9、 626.02、 672.83、 707.14、 745.98、 809.12、 820.07、 868.52、 892.9、 915.15、 939.71、 941.95、 953.16和959.89 nm)作为输入数据输入至已经建立的ELM模型中, 预测每个像素点的品位数据。 由于样本内存在石英等不包含铁元素的矿物, 且矿石表面在取样切割过程中存在少量破损, 因此将预测值高于100%或低于0%的像素点用黑色表示。 图12为部分鞍山式铁矿样品的品位预测可视化结果, (a)、 (b)和(c)三个样本的品位预测平均值分别为22.34%、 24.79%和42.28%, 品位实测值分别为19.48%、 21.63%和38.59%。

图12

Fig.12

	Figure Option ViewDownloadNew Window
	图12 MSC-CARS-ELM模型得到的鞍山式铁矿品位可视化预测图Fig.12 Visual grade prediction diagram of Anshan type iron ore obtained by MSC-CARS-ELM model

与基于点源高光谱数据所建矿石品位反演模型不同^{[9, 24]}, 本工作以高光谱成像仪采集的室内高光谱图像为数据源, 探究了基于成像高光谱数据实现鞍山式铁矿品位反演的最优预处理组合算法。 通过CARS提取特征波段, 在简化模型输入变量, 降低计算复杂性和提升模型精度的同时, 也有助于开发基于特定波段的多光谱成像系统, 以满足快速在线检测鞍山式铁矿品位的需求。 此外, 成像光谱特征与基于无人机平台的矿区大面积高光谱成像特征更为相近, 本研究结论可以更好的支持矿山大面域成像高光谱数据的品位反演。 然而, 目前研究还处于室内高光谱图像并局限于单一样本的品位反演与可视化阶段, 下一步将基于机载高光谱图像实现矿区大范围的矿岩分类与矿石品位提取工作。