光谱传感器所成的图像中的一个像元通常不止显示一种纯类地物，往往包含多种元素，这样就构成了混合像元^[1]。遥感图像中混合像元的存在不仅影响了基于高光谱图像的地物识别和分类精度，而且已经成为遥感科学向定量化发展的障碍^[2]。现有的数据解混技术大都基于线性光谱混合模型对图像进行分解。高光谱解混分为2个步骤：端元提取和各个像元对应的丰度估计^[3]。已经出现且常用的端元提取方法有纯净像元指数法(pixel purity index, PPI)^[4]、N-FINDR^[5-6]、迭代误差分析法(iterative error analysis, IEA)^[7-8]和顶点成分分析法(vertex component analysis, VCA)^[9]。丰度反演算法常用的是最小二乘算法(least square, LS)。非负矩阵分解^[10](non-negative matrix factorization, NMF)算法是当今较为流行的盲源分离(blind signals separation, BSS)算法之一，能同时进行端元提取和丰度反演，较之前的解混算法更为简便，广泛应用于图像处理。非负矩阵由Paatero等^[10]提出，此后，Lee等^[11]中也分别介绍和推广了NMF的概念。此后，相继提出了多种基于NMF的改进算法，Chen^[12]提出了平滑约束的非负矩阵分解(CNMF)，对非负矩阵模型中的端元矩阵和丰度矩阵分别添加平滑约束，提高了NMF算法的准确性，Miao在^[13]中提出最小体积约束的非负矩阵分解(MVC-NMF)，将高光谱数据中的所有像元包含在由端元形成的单形体中，通过最小化此单形体体积，进行NMF，提高了分解过程的抗噪性能，文献[14]中提出，将NMF扩展到L_1/2稀疏约束，使解混结果更稀疏和精确，孔繁锵等^[15]阐述了非凸稀疏低秩约束的高光谱解混方法，此方法更适用于大信噪比混合图像，Xu等^[16]提出了Pareto优化的稀疏约束的混合像元分解，此方法需要先验信息。以上算法的问题在于：最小体积约束的非负矩阵分解收敛速度非常慢，需要时间过长；稀疏约束的非负矩阵分解忽略了解空间的大小，使得算法在求解过程中造成空间的浪费，拉长了求解时间。本文提出了一种基于自然梯度下降的单形体体积及稀疏约束的非负矩阵分解方法，此算法是将最小体积约束与稀疏约束共同作用于高光谱混合的解混算法中，同时利用自然梯度下降进行求解，提高了计算效率。

1 非负矩阵分解 1.1 线性光谱混合模型

在处理高光谱图像时，通常将图像看成线性混合模型。线性混合模型可以表示为：X= WH + n，其中n表示L维噪声或误差，矩阵X^L×N表示波段为L，像元为N的图像，矩阵W表示端元数量为P的端元矩阵，矩阵H表示P个端元在原图像上所占的的比重(即丰度信息)。丰度信息未知，但须同时满足非负的条件，即W_{i, j}, H_{i, j}≥0，同时丰度矩阵要满足和为一的条件，即$ \sum\limits_{i = 1}^p {{\mathit{\boldsymbol{h}}_i}} = 1, j = 1, 2, \cdots , L$。

1.2 非负矩阵分解

NMF算法由来已久，从1997年Lee等在^[11]中首次提出后就得到了广泛地关注。非负矩阵分解算法能在未知纯净像元的条件下同时得到有意义的端元矩阵及其对应的丰度信息，目前已经在多个领域得到广泛应用。NMF是将一个非负矩阵X分解成2个非负矩阵W和H乘积的形式，W是L×P阶矩阵，H是P×N阶矩阵，P∈min(N，L)。求解非负矩阵最优解的过程目前有2种方式：1)欧氏距离^[17]，建立的目标函数为：$ \parallel \mathit{\boldsymbol{X}} - \mathit{\boldsymbol{WH}}{\parallel ^2} = \sum\limits_{i, j} {{{\left( {{\mathit{\boldsymbol{X}}_{i, j}} - \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{H}}_{i, j}}} \right)}^2}} $；2)K-L散度^[18]，建立目标函数：$ D(\mathit{\boldsymbol{A}}\parallel \mathit{\boldsymbol{WH}}) = \sum\limits_{i, j} {{{\left( {{\mathit{\boldsymbol{X}}_{i, j}} - \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{H}}_{i, j}}} \right)}^2}} $。最方便常用的是利用前者求解最小值来进行非负矩阵分解。W和H的迭代准则为:

$ \mathit{\boldsymbol{W}} \leftarrow \mathit{\boldsymbol{W}}. * \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{H}}^{\rm{T}}}./\mathit{\boldsymbol{WH}}{\mathit{\boldsymbol{H}}^{\rm{T}}} $

(1)

$ \mathit{\boldsymbol{H}} \leftarrow \mathit{\boldsymbol{H}}. * {\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{X}}./{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{WH}} $

(2)

每次迭代过程中的H都要满足H_{i, j}≥0，和$\sum\limits_{i = 1}^p {{\mathit{\boldsymbol{h}}_i}} = 1, j = 1, 2, \cdots , L $。

2 本文数据解混算法 2.1 丰度稀疏约束的NMF

单纯的NMF具有非凸性，在求解过程中容易陷入局部最优。不少学者提出了应对策略，其中对NMF添加丰度稀疏约束是有效且常用的一种解决方法。

添加丰度稀疏约束后的NMF能使解混结果更接近端元在真实场景中的分布情况。丰度稀疏约束的NMF的损失函数及约束条件表示为：

$ \left\{ {\begin{array}{*{20}{l}} {f\left( {\mathit{\boldsymbol{W}},\mathit{\boldsymbol{H}}} \right) = \frac{1}{2}\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{WH}}} \right\|_2^2 + \lambda C\left( \mathit{\boldsymbol{H}} \right)}\\ {{\rm{s}}{\rm{.t}}{\rm{. }}\mathit{\boldsymbol{W}} \ge 0,\mathit{\boldsymbol{H}} \ge 0} \end{array}} \right. $

(3)

式中：λ是一个大于0的标量，代表稀疏度量函数C(H)的加权值，在这里稀疏度量函数C(H)可表示为：$ C(\boldsymbol{H})=\|\boldsymbol{H}\|_{1}=\sum\limits_{i} \sum\limits_{j} H_{i j}$。根据文献[4]确定λ的值：

$ \lambda = \frac{1}{{\sqrt L }}\sum\limits_l {\frac{{\sqrt N - {{\left\| {{\mathit{\boldsymbol{y}}_i}} \right\|}_1}/{{\left\| {{\mathit{\boldsymbol{y}}_l}} \right\|}_2}}}{{\sqrt {N - 1} }}} $

(4)

式中：y_l为高光谱图像矩阵X中第l个波段的向量。λ值越高，说明图像中丰度稀疏度越高，像元混合度越低；反之丰度稀疏度越低，像元混合度越高。

2.2 改进的丰度稀疏约束NMF

混合像元的空间分布相对集中在凸面几何体的内部，当只考虑丰度稀疏约束NMF时，就忽略了像元在空间上的实际分布，这时会造成空间的浪费，加入体积最小约束可以有效缩小空间，减短数据解混时间，提高效率。体积在这里指端元到单形体体积中心的欧氏距离，那么体积最小约束项P(W)可以表示为：

$ P\left( \mathit{\boldsymbol{W}} \right) = \sum {\sum {{{\left( {\mathit{\boldsymbol{W}} - \mathit{\boldsymbol{\bar W}}} \right)}^2}} } $

(5)

式中：$\overline{\boldsymbol{W}}=\left(\frac{1}{r} \sum \boldsymbol{w}_{i}\right) \mathbf{1}_{r}^{\mathrm{T}} $是所有端元的平均值矩阵，代表所有端元构成的单形体体积的中心，$ \overline{\boldsymbol{W}} \in \mathbf{R}^{L \times 1}$。

图 1显示了添加体积最小约束前后的NMF的混合像元散点图，星点表示高光谱图像中所有的像元，圆点表示NMF的估计端元，从图 1中可以看出利用单纯NMF方法得出的估计端元形成的单形体所占的空间较大，高光谱图像中的像元并不会充满整个由估计端元形成的单形体，这样就会造成混合数据解混过程中空间的浪费。而添加体积最小约束后的NMF得出的估计端元形成的单形体不仅体积小，还能包含高光谱图像中的所有像元。这是因为添加了体积最小约束后的NMF在求解过程中，端元会逐渐向像元靠近，缩小了解混空间，加快解混速度，节省时间。

基于丰度稀疏和体积最小共同约束的NMF(sparseness and minimum volume constrained NMF，SMVC-NMF)的损失函数可表示为：

$ \left\{ {\begin{array}{*{20}{l}} {f\left( {\mathit{\boldsymbol{W}},\mathit{\boldsymbol{H}}} \right) = \frac{1}{2}\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{WH}}} \right\|_2^2 + \lambda C\left( \mathit{\boldsymbol{H}} \right) + \beta P\left( \mathit{\boldsymbol{W}} \right)}\\ {{\rm{s}}{\rm{.t}}{\rm{. }}\mathit{\boldsymbol{W}} \ge 0,\mathit{\boldsymbol{H}} \ge 0} \end{array}} \right. $

(6)

将目标函数分别对W、H求偏导，可得：

$ \frac{{\partial f\left( {\mathit{\boldsymbol{W}},\mathit{\boldsymbol{H}}} \right)}}{{\partial \mathit{\boldsymbol{W}}}} = \mathit{\boldsymbol{WH}}{\mathit{\boldsymbol{H}}^{\rm{T}}} - \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{H}}^{\rm{T}}} + \alpha {\mathit{\boldsymbol{W}}^{ - {\rm{T}}}} $

(7)

$ \frac{{\partial f\left( {\mathit{\boldsymbol{W}},\mathit{\boldsymbol{H}}} \right)}}{{\partial \mathit{\boldsymbol{H}}}} = {\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{WH}} - {\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{X}} + \lambda $

(8)

式(7)中出现W^-T，使迭代过程复杂且浪费时间，现利用自然梯度下降法^[19]对原式进行优化，得到利用自然梯度下降的丰度稀疏与体积最小共同约束的NMF(NG-MVSC-NMF)，损失函数表示为：

$ \frac{{\partial f\left( {\mathit{\boldsymbol{W}},\mathit{\boldsymbol{H}}} \right)}}{{\partial \mathit{\boldsymbol{W}}}} = {\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{W}} = \mathit{\boldsymbol{WH}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{W}} - \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^{\rm{T}}}\mathit{\boldsymbol{W}} + \alpha \mathit{\boldsymbol{W}} $

(9)

在黎曼空间中，自然梯度方向被证明是收敛速度最快的^[19]。自然梯度法的使用避免在混合像元分解过程的求逆运算，从而防止式(7)中W变为奇异阵，加快算法的收敛速度。算法收敛，则目标函数数值随迭代次数增加而减小，设置目标函数需要满足的最小值，在实验中函数逐渐收敛，目标函数将在满足条件后停止迭代。

本文NG-MVSC-NMF算法步骤：

输入   非负高光谱混合数据$\mathit{\boldsymbol{X}} \in {{\mathbf{R}}^{L \times N}} $；

输出    2个非负矩阵$ \boldsymbol{W} \in \mathbf{R}^{L \times r}, \boldsymbol{H} \in \mathbf{R}^{r \times N}$，且满足和为1约束 $ \mathbf{1}_{P}^{\mathrm{T}} \boldsymbol{H}=\mathbf{1}_{N}^{\mathrm{T}}$ 。

初始化：设置初始丰度矩阵H，设置迭代次数k，估计端元数量r，设置初始端元矩阵W重复式(9)、式(10)，交替更新W、H，对H进行归一化，直到达到最大迭代次数或满足停止要求后停止更新。

3 实验 3.1 评价指标

利用光谱角距离(spectral angle distance, SAD)^[20]和均方根误差(root mean square error, RMSE)来分析解混结果。光谱角距离：通过比较端元的真实光谱矢量与估计的光谱矢量的夹角大小来判断其相似程度。夹角越小相似度越高，反之越低。利用均方根误差来衡量丰度估计的精度：

$ {\rm{SA}}{{\rm{D}}_i} = {\arccos ^{ - 1}}\left( {\frac{{\mathit{\boldsymbol{W}}_i^{\rm{T}}{{\mathit{\boldsymbol{\hat W}}}_i}}}{{\left\| {\mathit{\boldsymbol{W}}_i^{\rm{T}}} \right\|\left\| {{{\mathit{\boldsymbol{\hat W}}}_i}} \right\|}}} \right) $

(10)

式中：$ {\mathit{\boldsymbol{W}}_i}$和$\hat{\boldsymbol{W}}_{i} $分别表示第i个端元的光谱解混结果和光谱数据库的参考值：

$ {\rm{RMS}}{{\rm{E}}_i} = \sqrt {\sum\limits_{j = 1}^N {\frac{1}{{N - 1}}} \left| {{\mathit{\boldsymbol{H}}_{kj}} - {{\hat H}_{kj}}} \right|_2^2} $

(11)

式中：$ {\mathit{\boldsymbol{H}}_{kj}}$和$ {{\mathit{\boldsymbol{\hat H}}}_{kj}}$分别表示该端元在某一像元的丰度的光谱结果和参考值。

3.2 合成数据实验

本次实验采用人为产生的合成数据，对比L₁-NMF、MVC-NMF和SMC-NMF与本文提出的算法。数据由USGS光谱库中的5种地物作为端元矩阵，分别是钙铁辉石、伊利石、橄榄石、微斜长石、绿脱石；丰度矩阵随机生成，服从Dirichlet分布，且满足“非负”与“归一化”约束；根据高光谱成像原理，添加高斯噪声。将端元矩阵与丰度矩阵相乘，加上高斯噪声得到合成数据。端元光谱信息如图 2。

3.2.1 算法在不同噪声情况下的解混效果

在抗噪性能实验中，使用的是分辨率为60×60，端元数目为5个，添加高斯噪声分别为25、35、45、55、65 dB的合成数据进行对比实验，验证在不同噪声情况下，NG-MVSC-NMF与其他对比算法的抗噪性能。实验均以50次迭代的方式结束。解混所得端元矩阵的光谱角距离SAD均值和丰度矩阵的均方根误差RMSE的均值计算结果以及算法所用时间如表 1所示。

在这个实验中，分别进行NG-MVSC-NMF与其他对比算法对添加不同噪声的合成高光谱图像的解混，此实验中，端元矩阵和丰度矩阵均已知，所以实验得出的SAD和RMSE结果精度较高。表 1显示了合成图像在添加不同噪声情况下利用本文算法解混与其他对比算法的SAD、RMSE以及所需时间的比较。

由表 1可以看出，各个算法在添加不同噪声情况下各自得到的解混结果大致相同，说明噪声的变化对各个算法分解过程的影响并没有很明显；在对比算法中，SMC-NMF得到的SAD结果最好，但RMSE却没有其他2种算法好，而本文算法得出的SAD和RMSE结果均比其他算法好，这说明NG-MVSC-NMF的抗噪性能是最好的。在相同分辨率、添加相同噪声情况下，本文提出的算法解混所需时间是最短的。

3.2.2 算法在不同像元情况下的解混效果

在本次实验中，使用的是添加噪声为60 dB，端元数目为5个，不同像元个数分别为20×20、40×40、60×60、80×80、100×100的合成数据进行对比实验，验证在不同像元情况下，NG-MVSC-NMF与其他对比算法的解混效果。实验均以50次迭代的方式结束。解混所得端元矩阵的光谱角距离SAD均值和丰度矩阵的均方根误差RMSE的均值计算结果以及算法所用时间如表 2所示。

在这个实验中，分别进行NG-MVSC-NMF与其他对比算法对不同像元合成高光谱图像的解混，此实验中，端元矩阵和丰度矩阵均已知，所以实验得出的SAD和RMSE结果精度较高。表 2显示了合成图像在不同像元情况下利用本文算法解混与其他对比算法的SAD、RMSE以所需时间的比较。

由表 2可以看出，SAD数值随着像元的增加而减小，这说明本文算法更适用于像元多的高光谱图像分解，与其他算法对比，相同像元情况下的SAD值最小，即解混得到的估计端元与真实端元信息最为靠近；在40×40像元的情况下，本文算法得出的解混效果显示出明显的误差，所以不做比较；而在其他数量的像元情况下，NG-MVSC-NMF算法得到对比其他算法有明显优势的解混效果；本文算法在相同像元数量下对混合像元分解所需要的时间更短，像元数越多，优势就越明显。

3.3 真实数据实验

为了验证本文算法的性能，还进行了关于真实高光谱数据的解混实验。本次实验采用1997年AVIRIS在美国内华达州采集的Cuprite区域的高光谱图像作为真实数据进行实验，如图 3所示。去除不利于实验的波段(1~2, 104~113, 148~167, 221~224)，剩余188个波段作为实验波段，图像共有250×191个像元。本组数据在USGS库中选取真实端元。本次实验选取了8个端元作为真实端元信息，但由于缺乏现场真实丰度数据信息，因此通过显示算法解混得到的地物分量图对丰度进行评估。

根据文献[14]，稀疏正则化参数λ取0.38，经过3种方法得到的解混结果如表 3及图 3所示，可以看出NG-MVSC-NMF解混所用时间比其他3种方法所用时间更少，且给出了更好的结果。

本次实验总共选取了8类端元进行实验：高岭石、明矾石、水铵长石、玉髓、黄钾铁矾、蒙脱石、白云母、榍石。表 3显示了各方法解混后8种估计端元与真实端元之间的SAD及各算法解混所需时间指标。从SAD平均值可以看出，NG-MVSC-NMF算法得到的估计端元与MVC-NMF结果相差不多，但比其他对比算法得到的估计端元更为接近真实端元信息。图 3显示解混后的各端元对应的丰度图，各种物质都清晰可见。

4 结论

1) 相对于以往常用NMF模型，本文采用基于稀疏及体积约束的NMF模型，介绍了一种基于自然梯度的求解算法。实验可看出该方法在解混精度上有所提高，得出的地物分量图十分清晰，所需时间显著减少。

2) 由于实际情况下在高光谱混合图像中很有可能存在异常点，或者噪声的影响使得顶点分析与真实端元相差甚远，因此如何保证算法的鲁棒性值得进一步研究。