高光谱遥感技术是近年来非常流行的技术之一，通过遥感平台形成的图像具有高光谱分辨率、低空间分辨率的特点，这样就造成了图像中的一个像元中不仅包含一种纯粹地物，往往含有多种地物，称为混合像元^[1][2]。非负矩阵分解(non-negative matrix factorization, NMF)^[3]算法是基于线性光谱混合模型(linear spectral mixture model，LSMM)^[4]进行展开的，是当今较为流行的盲源分离(blind signals separation, BSS)算法之一，广泛应用于图像处理。NMF最早由Paatero和Tapper^[3]发现，此后，Lee和Seung^[4]分别介绍和推广了它的概念。NMF是基于数据非负的条件下进行的，这样能更好的贴合实际物理意义，使研究更具现实意义。此后，多种基于NMF的改进算法相继提出，2005年Chen Z^[5]提出了平滑约束的非负矩阵分解，在端元矩阵和丰度矩阵加入平滑约束(CNMF)，提高了NMF算法的精度，2007年Miao Li等^[6]将最小体积约束加入到非负矩阵分解(MVC-NMF)，提高了分解过程的抗噪性能，2009年Qian Y等^[7]提出L_1/2稀疏约束的NMF，使解混结果更稀疏和精确。2014年魏一苇等^[8]提出改进的MVC-NMF。有些真实地物信息通过科学测量可以获得，但以上方法都没有利用真实地物信息，浪费数据资源。Tang等^[9]将端元矩阵分成两部分，一部分代表已知端元，另一部分代表未知端元，但缺点是求解时间长且结果不稀疏。本文利用已经获得的真实地物信息作为先验信息对非负矩阵分解进行优化，改善利用先验信息求解时间长的缺陷。

1 非负矩阵分解算法 1.1 LSMM

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

1.2 非负矩阵分解算法及改进

非负矩阵分解算法能在未知纯净像元的条件下同时得到有意义的端元矩阵及其对应的丰度信息，目前已经在多个领域得到广泛应用。非负矩阵分解是将一个非负矩阵X分解成两个非负矩阵W和H乘积的形式，W是L×P阶矩阵，H是P×N阶矩阵，P∈min(N，L)。求解非负矩阵最优解的过程目前有两种方式：1)欧氏距离^[9]，建立的目标函数为${\left\| {\mathit{\boldsymbol{X}}{\bf{ - }}\mathit{\boldsymbol{WH}}} \right\|^2} = \sum\limits_{i, j} {{{\left({{\mathit{\boldsymbol{X}}_{i, j}} - \mathit{\boldsymbol{W}}{\mathit{\boldsymbol{H}}_{i, j}}} \right)}^2}} $；2)K-L散度^[10]，建立目标函数为$D\left({\mathit{\boldsymbol{A|WH}}} \right) = \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}} \cdot \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{H}}^{\rm{T}}} \cdot /\mathit{\boldsymbol{WH}}{\mathit{\boldsymbol{H}}^{\rm{T}}} $

(1)

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

(2)

每次迭代过程中的H都要满足H_{i, j}≥0和$\sum\limits_{i = 1}^P {{h_i}} = 1$, j=1, 2, …, L。

随着高光谱数据解混研究的不断深入，真实高光谱图像中的某些特征的先验知识是可以获得的。在矿区等多个场景中，铁、云母石等材料通过现场调查极易识别，因此它们的光谱特征可以从光谱库中获得。由于现今的大多数分解算法在建立模型时并没有考虑这些真实端元的已知条件，造成数据浪费及解混结果精度偏低。本文提出基于已知端元约束的非负矩阵分解的高光谱解混方法则可充分利用已知端元信息。

本文算法根据文献[10]建立基本模型为：

$f(\mathit{\boldsymbol{W}}, \mathit{\boldsymbol{H}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{X}}{\bf{ - }}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{H}}_1}{\bf{ - }}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{H}}_2}} \right\|_2^2 $

(3)

式中：W₁、H₁分别代表已知端元矩阵和丰度矩阵；W₂、H₂分别代表未知端元矩阵和丰度矩阵。添加约束后模型为

$\begin{array}{*{20}{c}} {f(\mathit{\boldsymbol{W}}, \mathit{\boldsymbol{H}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{X}}{\bf{ - }}{\mathit{\boldsymbol{W}}_1}{\mathit{\boldsymbol{H}}_1}{\bf{ - }}{\mathit{\boldsymbol{W}}_2}{\mathit{\boldsymbol{H}}_2}} \right\|_2^2{\bf{ + }}}\\ {\frac{\lambda }{2}C\left( {{\mathit{\boldsymbol{W}}_1}} \right) + \frac{\beta }{2}S\left( {{\mathit{\boldsymbol{H}}_{\rm{1}}}} \right) + \frac{\gamma }{2}S\left( {{\mathit{\boldsymbol{H}}_2}} \right)} \end{array} $

(4)

式中：C(W₁)表示单形体体积最小约束；S(H₁)表示丰度稀疏约束；λ、β、γ分别代表体积最小约束和丰度稀疏约束的权值，均大于0。体积约束权值由读者自行设置。根据文献[12]中提出的稀疏度求解方法，可将丰度稀疏约束权值用公式$\beta、\gamma = \frac{1}{{\sqrt L }}\frac{{\sqrt N - {{\left\| {{x_l}} \right\|}_1}/{{\left\| {{x_l}} \right\|}_2}}}{{\sqrt {N - 1} }}$来定义，其中x_l代表高光谱数据矩阵X的第l个波段向量，由于是在同一图像内分别计算丰度信息，因此上述两个稀疏约束权值β、γ相等。

由于端元不能随意确定，当现场的测量端元与估计的端元相差太大时也会产生极大的误差，因此本文将真实端元与初始化端元之间范数相对小的数据作为已知端元，保证解混结果的精度。

本文算法步骤如下：

输入  高光谱图像数据X∈R^L×N

利用顶点成分分析(vertex component analysis，VCA)和全约束最小二乘(full constraint least square，FCLS)对端元和丰度进行初始化，计算真实端元与初始端元的范数，选出已知端元；根据式(4)处理数据；

输出  W, H

2 实验及结果分析 2.1 评价指标

在高光谱解混中，常见的评价指标有光谱角距离(spectral angle distance, SAD)和均方根误差(root mean square error, RMSE)^[13]。SAD表示估计端元与初始端元之间的光谱角距离，RMSE表示端元对应的估计丰度信息与真实丰度之间的均方根误差。光谱角距离越小，均方根误差值越小，说明解混效果越好。在合成数据中可以同时利用以上两种评价指标来评判解混效果，在真实高光谱数据中，由于缺乏现场真实地物的分布图，不能利用RMSE进行衡量。

SAD表示为：

${\rm{SA}}{{\rm{D}}_i} = {\cos ^{ - 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) $

(5)

式中W_i和Ŵ_i分别表示第i个端元的光谱解混结果和光谱数据库的参考值。

RMSE表示为：

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

(6)

式中H_kj和Ĥ_kj分别表示该端元在某一像元的丰度的光谱结果和参考值。

2.2 合成高数据的实验

从美国地质勘探局(USGS)测量的光谱数据库选取不同端元合成模拟数据，此光谱数据库含224个波段，501个端元，波长范围为0.38~2.5 μm，光谱分辨率达10 nm。本文选取其中5种端元作为端元矩阵；随机生成相应的丰度矩阵，同时满足“非负”与“归一化”双重约束，且服从Dirichlet分布。合成数据中含有真实端元信息和真实端元分布情况，因此在对解混结果进行评价时采用SAD和RMSE。

2.2.1 不同像元的合成数据的实验

为了对比不同像元个数的高光谱数据对解混结果的影响，在本次实验中，利用不同像元个数分别为20×20、40×40、60×60、80×80、100×100的合成数据进行对比实验。当实验达到最高迭代次数停止运算。本文给出了解混所得端元矩阵的光谱角距离SAD均值和丰度矩阵的均方根误差RMSE的均值计算结果。

在这个实验中，分别进行本文算法与其他对比算法对不同像元合成高光谱数据的解混，此实验中，端元矩阵和丰度矩阵均已知，因此实验得出的SAD和RMSE结果精度较高。表 1和表 2分别给出了合成数据在不同像元情况下利用本文算法解混与其他对比算法的SAD和RMSE比较。

由表 1可以看出，本文算法的SAD数值始终比其他对比算法要低得多，即解混得到的估计端元与真实端元信息最为靠近，在80×80像元的情况下，本文算法得出的解混效果最好；由表 2可以看出本文算法得出的丰度信息更为接近真实丰度信息，解混效果最好；由表 3可以看出本文提出算法相比其他算法对数据解混所需时间是最少的，但是比单纯NMF算法时间长，这是因为NMF算法模型简单，且容易陷入局部最优。

2.2.2 不同噪声的合成数据的实验

为了判断本文算法的抗噪性能，实验采用添加不同噪声的合成数据进行解混，并与其他算法的解混效果进行对比。合成数据分别加入25、35、45、55、65 dB的噪声，当实验达到要求或到达最大迭代次数停止迭代，本文给出了解混所得端元矩阵的光谱角距离SAD均值和丰度矩阵的均方根误差RMSE的均值计算结果。

在这个实验中，分别进行本文算法与其他对比算法对添加不同噪声的合成高光谱数据的解混，端元矩阵和丰度矩阵均已知，因此实验得出的SAD和RMSE结果精度较高。表 4和表 5给出了合成数据在添加不同噪声情况下利用本文算法与其他对比算法解混结果的SAD和RMSE比较。由表 4可以看出，本文算法的SAD数值始终比其他对比算法要低得多，即解混得到的估计端元与真实端元信息最为靠近，在添加噪声为65 dB时，本文算法得出的解混效果最好；由表 5可以看出随着噪声增大，本文算法的解混效果整体呈上升趋势，得到的估计丰度信息更为接近真实丰度信息，解混效果最好；由表 6可以看出本文提出算法比其他对比约束NMF算法对数据解混所需时间更少，但是比单纯NMF算法时间长，这是因为NMF算法模型简单，容易陷入局部最优。

2.3 真实高光谱数据的实验

真实数据选自1997年AVIRIS在美国内达华州采集的Cuprite地区的高光谱图像数据。滤除对实验不利或无用的波段后，高光谱数据为X∈R^L×N, L=180, N=250×191，即共有188个波段，47 750个像元。根据实际测量的地物信息，选取其中8种地物作为真实端元进行解混。由于缺乏真实场景下的端元分布信息，利用SAD和丰度图对解混结果进行评价。真实高光谱解混后丰度如图 1。

图 2是本文算法和其他对比算法对高光谱解混后的SAD比较。由图 1、2可以看出本文算法得到的解混效果更好。

通过真实数据解混实验得到各个算法所需时间分别是42.082 3 s(NMF)、2 412.759 820 s(MVC-NMF)、913.002 5 s(L_1/2-NMF)、37.708 0 s(本文算法)。相比可知，本文算法对真实数据解混所需时间更短。

3 结论

本文在NMF方法来解决高光谱解混问题的基础上对其进行优化，添加先验信息约束，使解混结果更接近真实情况。根据实验数据分析，NMF及MVC-NMF算法收敛效果差，L_1/2-NMF算法对空间信息利用不足，本文算法克服了以上的缺点，使计算过程中有较好的收敛效果，同时很大程度地利用了先验信息。仿真实验结果表明，相对于MVC-NMF算法和L_1/2-NMF算法，本文提出的算法具有更好的抗噪性能，即使在多个像元情况下，也有较好的解混效果。如何保证算法的稳定性值得进一步研究。