矿物混合后光谱值会受到矿物本身颜色、 粒径和比重的综合影响; 混合后光谱特征可分为线性混合和非线性混合。 试验发现, 混合矿物的光谱反射率近似为混合单矿物光谱反射率的线性混合, 除了矿物之间混合非常紧密的, 认为是非线性混合, 其他的均可以视为是线性混合^[13]。 假设某种混合物的高光谱反射率曲线(后文统称光谱曲线)是所含组分光谱曲线按比例的线性混合, 其线性混合模型定义为

mv=As+ε(1)

式(1)中, m_v∈ R^{t× 1}表示一个混合光谱向量, 其中t为光谱波段数; A∈ R^{t× u}表示各个组分光谱向量, u表示组分个数; s∈ R^{u× 1}表示组分丰度矩阵。 ε ∈ R^{t× 1}表示噪声, t维高斯随机噪声或模型误差, 或不确定性因子。

假设物质样本个数为N, 每个样本的混合组分物质类型相同, 但是组分间混合比例不相同。 则N条样本实测光谱曲线就是组分光谱信号按照比例的线性混合, 其混合模型的矩阵形式定义为

Mv=AS+ε(2)

式(2)中, M_v∈ R^{t× N}表示混合光谱向量, 其中t为光谱波段数; A∈ R^{t× u}表示组分光谱基向量, u表示组分个数; S∈ R^{u× N}表示组分丰度矩阵。 ε ∈ R^{t× N}表示噪声, t维高斯随机噪声、 模型误差或不确定性因子。

根据光谱曲线的实际物理意义, 每个组分光谱向量需要满足非负性的约束条件, 一个样本组分混合的丰度矩阵需要同时满足非负性以及和为1的约束条件。

∑Sij=1, ∀i, j, Sij> 0(3)

1.2.1 加权NMF解混目标约束

对于光谱混合模型, 在未知组分源光谱反射曲线和组分比例的条件下进行盲源解混分离, 需要构建一定的解混目标约束条件, 才能实现数学向量分离解混。 光谱线性混合模型的解混过程可以近似认为, 是在目标约束条件下求解源光谱基向量A和丰度S的解析过程。 常用函数有Fast ICA和NMF算法等, 其中非负矩阵分解算法是在ICA算法的基础上进一步优化, 对于非负信号分离检测有较高的分离精度。 基于加权NMF组分光谱解混是从假设未知单一源光谱基向量A和混合比例S两个非负矩阵开始, 经过权函数处理算法, 不断交替迭代趋近最小化目标函数的过程。 设立矿物混合光谱解混的NMF的目标函数J(S, W), 以最小欧氏距离和重加权稀疏建立组合约束条件, 定义为

J(S, W)=minA, S12‖M-v-SW‖2+‖W.* S‖1(4)

式(4)中, W表示控制迭代分解后分量丰度的权重矩阵, ‖ · ‖ ₁表示向量的1范数, .* 表示矩阵对应元素相乘, (· )^T表示向量的转置。 整个解混过程实质就是要找到一个最优方向的W, 使得该方向上的非高斯性最大, 进而求解组分源光谱基向量A。

在进行混合光谱解混过程中需要关注两个步骤, 一个是实测光谱的均值白化预处理, 另一个是解混过程初始权值设定和迭代权值计算。

1.2.2 均值白化预处理

一般在光谱混合模型中, 假设源光谱的各个分量为零均值随机向量。 为了达到分离出的源信号的估计值满足零均值假设, 需要对高光谱数据集M_v进行均值化处理得到矩阵M'_v

M'v=Mv-mean(Mv)(5)

另外, 对输入混合光谱进行白化(Whiting)处理。 白化处理是指通过白化矩阵Q与混合光谱乘积变换实现数据白化。 白化目的是去除数据中的冗余信息, 有效降低解混问题的复杂度, 使解混出组分源光谱的基向量的信号实现非高斯性的最大化, 并且对数据的信息熵不产生影响。 奇异值分解法(singular value decomposition, SVD)是数据白化常用方法, 能比较精确对观测数据进行白化处理。 SVD白化过程是首先计算待解混M'_v矩阵的协方差矩阵T, 然后对T进行奇异值正交化分解, 得到正交化矩阵E和特征值对角矩阵D。 然后计算白化矩阵Q为

Q=ED-12ET(7)

混合光谱白化处理后得到 M-v

M-v=QM'v=Q[Mv-mean(Mv)](8)

1.2.3 权重矩阵确定

加权NMF光谱解混模型是在合适的权值向量下, 通过加权正则化迭代计算获得更稀疏的解。 为了提高算法解混的分离精度, 消除数据不确定性对总体解混结果的影响, 模型从初始权确定方法和迭代权函数两个方面进行了优化。 一般解混模型初始权会选取单位矩阵进行迭代。 文中初始权值W₀为了提高模型收敛速度, 初始权为混合光谱与组分光谱基向量光谱角余弦值, 光谱角输出结果值越小表示两个光谱相似度越高, 其W₀对应初始权值越接近于1; 反之, 表示两个光谱相似度越小, 其W₀对应初始权值越接近于0。 如果未知组分源光谱的情况下, 可以设定每个组分的初始权值相等。

W0=cos(RAD)+δ=MTvA^‖MTv‖‖A^‖+δ(9)

式(9)中, M_v指某个混合物实测光谱向量, A^指标准波谱数据库的参考高光谱向量。 式RAD度量是两个向量之间的光谱角距离, ‖ · ‖ 表示向量的2范数, (· )^T表示向量的转置, δ 是一个正小常数, 防止权重出现0, 导致解算异常。

混合光谱解算过程的迭代权函数, 是根据上一次解算丰度矩阵来计算下一次迭代的权重,

Wi, j{k+1}=1/(|Si, j{k}|+eps)(10)

式(10)中, Si, j{k}表示第k次迭代的丰度矩阵, i、 j分别表示第i个混合光谱中在第j个组分中所占的比例, 表示元素相除。 eps表一个正的常数, 防止分母为0。

NMF盲源解混方法的最大特点在于它可以充分考虑源信号的统计独立性、 稀疏性、 时空无关性等特性来估计不同信号源光谱, 是一种稳健高效的算法。 但是解算结果存在两种不确定性, 一个是分离出光谱基向量排列(permutation)顺序的不确定, 即无法了解所恢复的光谱基向量, 对应A中的哪个组分源光谱; 另一个是解算出的光谱基向量具有信号尺度(scaling)不确定性, 即无法恢复出光谱基向量丰度或能量矩阵S, 与已知源光谱反射率之间存在数值差异和符号不确定性。 应用解混模型时, 通常预知了混合物中主要的组分构成, 实测了已知组分源光谱或者在标准波谱库组分光谱特性曲线。 然后与NMF算法解算出的组分光谱基向量进行识别排序计算。 根据X_SAD大小来识别分解光谱基向量与源光谱的对应关系, X_SAD表示真实光谱基向量与解混恢复光谱基向量之间的相似度, 其定义为

XSAD=ATA^‖AT‖‖A^‖(11)

式(11)中, A∈ R^{t× u}表示解混未知排列组分的光谱基向量, u表示组分个数, t为光谱波段数; A^∈ R^{t× u}表示已知排列顺序标准光谱数据库组分的光谱基向量; ‖ · ‖ 表示向量的2范数, (· )^T表示向量的转置。 X_SAD∈ (0, 100%), 当X_SAD≥ 99.5%时, 认为两个基向量为同一类光谱, 重新光谱排序。

使用本文提出的加权NMF高光谱组分盲源解混算法的流程(见图1)如下:

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 基于NMF的混合矿物高光谱曲线的盲源解混流程Fig.1 Blind source demodulation process of NMF-based hyperspectral curve of mixed minerals

Step1: 得用式(8)对混合物高光谱数据集M_v均值白化处理, 计算待解混矩阵 M-v;

Step2: 得用式(9)进行M_v与参考光谱的相关性分析, 计算解混的初始权阵W₀;

Step3: 以W₀为输入权阵, 利用式(4)对 M-v进行加权NMF迭代分解;

Step4: 利用式(10)计算迭代权重, 并以此作为Step3的输入权阵;

Step5: 重复执行Step3— Step4, 直到满足终止条件;

Step6: 输出光谱基向量, 利用式(11)判断组分基向量并重析排序;

Step7: 以重新排序的A为输出端W估算各个组分丰度。

为了分析基于加权NMF的光谱解混模型的实际解混效果, 选取了3组物质开展光谱解混分离的检测试验。 1组化合物为化学纯的氧化铜和氧化亚铜(图中氧化铜为A、 氧化亚铜为B), 2组化合物为化学纯的氢氧化铜和碱式碳酸铜(图中氢氧化铜为C、 碱式碳酸铜为D), 3组化合物为化学纯的蓝铜矿和孔雀石(图中蓝铜矿为E、 孔雀石为F)。 试验样本制作按照粒径相同、 总摩尔质量不变的原则, 设计混合物两种物质比例, 见表1。 样本制作时, 采用千分位精密电子天平称取一定重量两种物质粉末, 置入量杯用电动搅拌棒充分混合后, 倒入黑色纸盒, 震荡至样本表面平整。

表1

Table 1

表1(Table 1)

表1 试验1组和2组样本的物质摩尔质量配比 Table 1 Molar ratios of substances in samples 1 and 2

表1 试验1组和2组样本的物质摩尔质量配比 Table 1 Molar ratios of substances in samples 1 and 2

采用ASD Field Spec4高光谱仪测量样本的高光谱反射曲线。 光谱仪测量波段范围为350~2 500 nm, 采样间隔为1.4 nm/350~1 000 nm, 2 nm/1 000~2 500 nm。 因为矿物混合时会受到密度不同、 震荡不完全、 混合不充分等因素的影响, 即使同一个样本的光谱值也会存在微小差异。 为了消除采样光谱的差异, 一个样本需要测量多个点光谱值, 将多点平均值作为最终的光谱值。 光谱测量在正常光照的室内进行, 每个样本采用“ 蛇形” 选取60测量点, 每个测量点利用光谱矿物探头垂直于样本并且保持2 mm左右的距离进行测量, 样本光谱值是60点采样曲线的平均值。 实测光谱数据经过断点修复, 均值平滑等预处理后, 成为NMF高光谱解混的原始数据集。 不同比例混合物样本的实测光谱曲线都需要进行去均值和白化预处理, 见图2所示。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 (a) CuO与Cu2O混合光谱曲线; (b) Cu(OH)2与Cu2(OH)2CO3混合光谱曲线; (c)蓝铜矿与孔雀石混合光谱曲线Fig.2 (a) spectra of CuO and Cu2O mixtures; (b) spectra of Cu(OH)2 and Cu2(OH)2CO3 mixtures; (c) spectra of malachite and Azurite mixtures

图2(a)为CuO与Cu₂O混合物的原始光谱曲线。 图中红色曲线为Cu₂O光谱曲线, 在550 nm波段处有明显Cu⁺吸收谷; 灰色曲线为CuO光谱曲线, 在850 nm时Cu²⁺吸收位置; 两种化合物光谱曲线都在2 200 nm处羟基吸收谷。 图2(b)Cu(OH)₂与Cu₂(OH)₂CO₃混合物的原始光谱曲线。 图中灰色的曲线是Cu(OH)₂光谱曲线, 在480 nm有铜离子吸收波峰, 1 450 nm有较深水分子吸收谷; 红色曲线为Cu₂(OH)₂CO₃光谱曲线, 在550 nm有铜离子吸收波峰; 两种化合物光谱曲线在1 940~2 000 nm有多处水和羟基的吸收谷, 在2 200~2 400 nm有碳酸根振动吸收。 图2(c)蓝铜矿与孔雀石按照表2的摩尔质量配比混后物光谱曲线。 其主要吸收位置与Cu(OH)₂和Cu₂(OH)₂CO₃混合物光谱曲线相近, 但对应吸收位置峰值有差别。

表2

Table 2

表2(Table 2)

表2 试验样本蓝铜矿与孔雀石摩尔质量配比 Table 2 Molar ratio of Azurite to malachite in samples

表2 试验样本蓝铜矿与孔雀石摩尔质量配比 Table 2 Molar ratio of Azurite to malachite in samples

从图中可以看出, 两种化合物按照不同比例混合后的光谱曲线, 其反射率值会在这两种化合物光谱曲线之间进行波动, 整体呈现明显的按比例递增的变化规律。

光谱解混试验的结果采用三种指标进行评价: 解混性能指数PI, 分解后光谱相似度的光谱角距离指标D_SAM和解混光谱曲线精度指标X_RMSE, 指标具体数学公式如下。

(1)解混性能指数PI是表示混合光谱分离性能指标, 当解混完全时PI=0。 实际上当PI接近10^-2时说明该算法分离性能好。

PI=1n(n-1)∑i=1n{∑i=1n|gik|maxj|gij|-1+∑k=1n|gki|maxj|gji|-1}(12)

式(12)中, g_ij为解混后全局矩阵C=AS的元素; maxj|g_ij|表示C的第i行中绝对值最大值; maxj|g_ji|表示第j列元素绝对值最大值。

(2)光谱角距离D_SAM表示为解混光谱基向量与已知标准光谱之间相似程度。 当输出结果值越小表示两个光谱越匹配, 相似度越高; 反之, 表示两个光谱距离越大, 相似度越小。

DSAM=arccosAA^(ATA^)1/2(ATA^)1/2(13)

式(13)中, A某种物质已知标准光谱基向量, A^(λ )为解混后某种物质光谱基向量, (· )^T表示向量的转置。

(3)均方根误差X_RMSE是用来评价实际光谱与解混光谱之间的误差水平。 其X_RMSE数值越接近于0, 表明光谱解混丰度的误差越小, 精度越高; 反之, 表明解混的误差越大, 精度越低。

XRMSE=1N∑i=1N[A(i)-A^(i)]2(14)

式(14)中, N为光谱采样的总波段数, i为采样波段, A(i)波段区间的某种物质已知标准波谱反射率, A^(i)为解混之后波段区间的某种物质波谱反射率。

图3是加权NMF的光谱解混的试验结果, 是3组试验的源光谱曲线和解混后光谱曲线。 从图中可以看出, 解混后光谱曲线总体上能够保持源光谱的吸收位置特征, 但是吸收深度、 吸收宽度与源光谱存在一定的差异。 从图3(a)解混的两条光谱曲线基本保持了源光谱典型的吸收位置和吸收峰, 光谱整体变化趋势大体一致。 CuO解混后光谱与源光谱符合度很高, 两条曲线基本重合。 Cu₂O解混后光谱与源光谱总体变化相似, 两者反射率数值相差较大; 图3(b)中解混后Cu(OH)₂、 Cu₂(OH)₂CO₃光谱与源光谱的典型吸收峰变化相同, 但是两者的吸收位置存在微小偏移; 图3(c)中解混后蓝铜矿、 孔雀石光谱与源光谱典型吸收峰的变化相同, 两者吸收位置偏移较大, 如孔雀石源光谱550 nm吸收峰, 在解混后光谱吸收峰发生左偏至510 nm。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 (a) CuO与Cu2O混合解混光谱曲线; (b) Cu(OH)2与Cu2(OH)2CO3混合解混光谱曲线; (c) 蓝铜矿与孔雀石混合解混光谱曲线Fig.3 (a) spectra of de-mixing result of CuO and Cu2O mixtures; (b) spectra of de-mixing result of Cu(OH)2 and Cu2(OH)2CO3 mixtures; (c) spectra of de-mixing result malachite and Azurite mixtures

综合分析3组光谱解混的结果可以得出, 加权NMF解混模型的盲源解混效果十分明显。 在未知混合光谱的先验条件基础上, 仍然可以准确分离出基源光谱。 解混后光谱与源光谱的曲线整体变化趋势相同, 但是具体的吸收位置和吸收峰存在一定偏移, 对应吸收位置的反射率也存在明显的差异。

分析这种差异现象存在的原因有: (1) 源光谱是在实验室理想条件下测量获得的, 而解混光谱的基向量是数学解析的结果, 两种光谱获取条件不同, 存在一定的差异是难免的。 由于NMF解混模型中, 输入光谱数据是经过去均值和白化处理后的, 输入数据的本身存在一定的失真形变, 解混结果的反射率自然会产生相应的差异; (2) NMF模型解混原理是寻找一个基向量组对混合光谱数据进行线性逼近, 也就是解混基向组只能够代表数据间潜在相对的结构关系, 不是绝对差异的体现; (3) NMF算法解混要求数据具有非负性, 且一般为超高斯或亚高斯信号。 本试验数据是350~2 500 nm全波段反射光谱, 未对数据进行分块高斯性判定, 在一定程度上会导致解混后光谱精确度, 产生一定偏移和差值。

选用分离精度和均方差评价试验解混的效果, 见表3。 利用式(12)实际计算的性能指数PI来衡量解混效果, PI计算值是分离矩阵W与实际产混合矩阵A的乘积C=WA(C常被称为全局矩阵)与单位矩阵I之间的差值。 分离精度PI数值越接近0, 分离效果越准确。 均方差是解混光谱与源光谱对应波段反射率差值的均方差, 数值大小说明两者反射率差异。 表中3组样本分离精度PI都为0.15左右, 说明NMF模型光谱解混分离较好, 能够正确解混出与源光谱相似的基向量。 其中氧化铜和氧化亚铜的解混精度PI值最小, 达到了0.121 3, 分离度最高。 但是3组样本分离光谱的均方差数值较大, 说明解混光谱与源光谱反射率存在明显差异。 其中CuO解混光谱与源光谱相似度好, 均方差为0.026 8, 结果与图2中两者基本重合图形相吻合; Cu₂O解混均方差为0.369 4, 图2中解混光谱与源光谱反射率确实差异较大。 但是表中光谱角数值差异较多, 说明解混光谱和源光谱存在相似, 但还有一定偏移。

表3

Table 3

表3(Table 3)

表3 光谱解混分离效果的性能指标 Table 3 Performance indicators of spectral de-mixing and separation effect 混合物 分离精度PI XRMSE 光谱角距离 氧化铜 氧化亚铜 0.121 3 0.026 8 0.369 4 0.025 0 0.118 3 氢氧化铜 碱式碳酸铜 0.147 7 0.209 5 0.128 9 0.141 5 0.093 6 蓝铜矿 孔雀石 0.137 5 0.107 1 0.188 2 0.051 2 0.076 8

表3 光谱解混分离效果的性能指标 Table 3 Performance indicators of spectral de-mixing and separation effect

由于实际生产现场的检测物质不一定是化学纯化合物的混合, 还会夹杂少量未知物质。 因此生产现场检测目标的光谱曲线, 是多种已知化合物的光谱曲线与夹杂少量未知物质光谱的线性混合。 假设少量的未知物质光谱值为光谱数据的干扰噪声, 试验时将实测混合光谱加入一定比例的噪声, 模拟复杂环境下的混合光谱数据。 添加噪声的方法是向已知光谱数据直接添加某个信噪比(SNR)的噪声。 SNR是光谱信号强度除以噪声的强度, 信噪比越大表明噪声越小, 见式(15)

SNR=10lgE(s)E(n)=20×lg1ε(15)

式(15)中, E(s)和E(n)分别是信号和噪声的能量, ε 为噪声的百分比。

为了测试加权NMF光谱解混模型抗噪性能或者算法适应性, 以表2摩尔质量配比的蓝铜矿与孔雀石混合物光谱曲线为基础, 分别加入5%~15%的高斯噪声制作成新的试验数据。 其中5%高斯白噪声, 信噪比SNR为26.2 dB; 10%高斯白噪声, 信噪比SNR为19.8 dB; 15%高斯白噪声, 信噪比SNR为16.5 dB, 图4为添加噪声后的蓝铜矿与孔雀石混合光谱曲线。 从图中可以看出, 不同比例蓝铜矿与孔雀石混合后光谱加入噪声后, 光谱数据仍然呈现一定的规律性。 添加大噪声的光谱曲线的反射率值浮动明显, 添加小噪声反射率值只有轻微波动。 添加噪声后光谱数据与源数据的光谱曲线的整体变化趋势相同, 加入噪声没有明显影响光谱数据的吸收位置, 但对吸收峰值影响较大, 曲线间界线变得不清晰。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 (a) 蓝铜矿与孔雀石加5%噪声混合光谱曲线; (b) 蓝铜矿与孔雀石加10%噪声混合光谱曲线; (c) 蓝铜矿与孔雀石加15%噪声混合光谱曲线Fig.4 (a) spectra of malachite and Azurite mixtures +5% noise; (b) spectra of malachite and Azurite mixtures +10% noise; (c) spectra of malachite and Azurite mixtures +15% noise

图5是采用加权NMF方法光谱解混后的蓝铜矿与孔雀石光谱曲线。 从图中可以看出, NMF解混模型可以分解出组分源光谱的基向量, 解混后光谱曲线与源光谱曲线变化特征总体一致, 还存在一定扰动噪声。 解混后光谱曲线保持了铜离子吸收波峰、 1 940~2 000 nm有多处水和羟基的吸收谷, 在2 200~2 400 nm处有碳酸根吸收特征, 但吸收深度、 吸收宽度存在差异。 加入10%的噪声后蓝铜矿和孔雀石的解混光谱与源光谱特征接近, 依旧可以识别。 当噪声达到15%时, 铜离子和碳酸根吸收特征仍然明显, 而水和羟基的吸收谷特征已被噪声淹没, 说明加权NMF的光谱解混算法适应性有一定抗噪上限。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 (a) 蓝铜矿与孔雀石加5%噪声混合解混光谱曲线; (b) 蓝铜矿与孔雀石加10%噪声混合解混光谱曲线; (c) 蓝铜矿与孔雀石加15%噪声混合解混光谱曲线Fig.5 (a) spectra of de-mixing of malachite and Azurite mixtures +5% noise; (b) spectra of de-mixing of malachite and Azurite mixtures +10% noise; (c) spectra of de-mixing of malachite and Azurite mixtures +15% noise

表4为加噪后蓝铜矿与孔雀石混合光谱解混的指标。 从表中可以看出, 该算法整体分离效果较好, 解混抗噪能力较强。 随着混合光谱噪声的增大, 解混分离精度随之缓慢降低, 具有较强的噪声适应性。 当无噪声时, 分离精度PI为0.137 5, 当噪声达到15%时, 解混分离精度PI为0.181 4, 分离指标总体下降不到0.03。 另外, 随着混合光谱噪声的增大, 解混后光谱与源光谱的均方根误差数值波动很小。 加入噪声解混出孔雀石光谱均方根误差波动未超过0.02, 蓝铜矿光谱均方根误差波动未超过0.05, 总体差异不突出, 说明NMF解混算法具有较好鲁棒性。

表4

Table 4

表4(Table 4)

表4 光谱解混分离效果中性能指标 Table 4 Performance indicators of spectral de-mixing and separation effect 蓝铜矿和孔雀石 分离精度 PI 光谱角距离 XRMSE 孔雀石 蓝铜矿 孔雀石 蓝铜矿 混合原始光谱 0.137 5 0.076 8 0.051 2 0.188 2 0.107 1 混合光谱+5%噪 0.142 6 0.160 3 0.094 6 0.191 2 0.144 2 混合光谱+10%噪 0.162 3 0.222 4 0.096 2 0.181 3 0.151 3 混合光谱+15%噪 0.181 4 0.175 1 0.125 0 0.173 1 0.133 1

表4 光谱解混分离效果中性能指标 Table 4 Performance indicators of spectral de-mixing and separation effect