CN102193090A - 一种遥感图像混合像元分解方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，具体涉及一种新的基于非负矩阵分解算法的混合像元分解方法。本发明方法根据高光谱图像光谱和丰度的特点，在非负矩阵分解算法的目标函数中引入丰度分离性和平滑性的约束条件，并且在适当的时机移除这些约束条件并继续迭代，从而克服了NMF算法容易陷入局部极小的缺点，可以有效的解决高混合遥感数据的混合像元分解方法。本方法在基于多光谱和高光谱遥感图像的高精度的地物分类以及地面目标的检测和识别方面具有特别重要的应用价值。

Description

一种遥感图像混合像元分解方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种基于非负矩阵分解算法的可解决高混合遥感数据混合像元分解问题的新方法。

背景技术

遥感是本世纪六十年代发展起来的新兴综合技术，与空间、电子光学、计算机、地理学等科学技术紧密相关，是研究地球资源环境的最有力的技术手段之一。近年来，随着成像技术的进步，多波段遥感图像在越来越多的领域得到了广泛应用。由于成像系统空间分辨率的限制和地表的复杂多样，所获得的遥感图像中的一个像元往往包含着多种地物类型，这就形成了混合像元。如何从混合像元广泛存在的多波段遥感图像中准确的提取端元信号，并有效的对混合像元进行分解，已成为了遥感图像定量分析的一个重要研究课题[1]。

目前的混合像元分解算法可以分为两大类。一类是基于几何学的方法，另一类是基于统计学的方法[2]。几何学方法所基于的数学原理是，高光谱图像的所有数据点都位于一个单形体中，单形体的顶点分别对应各端元，因此求解端元光谱相当于找寻这个凸集的顶点。由于几何学方法是从已有的数据集中寻找顶点，因此不适用于没有纯像元的数据集。而基于统计学的方法克服了这个缺点，它充分利用数据的统计特性计算出端元光谱。

非负矩阵分解(Nonnegative Matrix Factorization，NMF)[3]是一种统计学的方法，它是Daniel.D.Lee等在1999年的Nature中首先提出的方法[3]，并用于解决人脸识别和语义分析中的问题。它尝试将一个非负矩阵分解成为两个非负矩阵的乘积。该算法的模型刚好与高光谱图像的线性混合模型一致，因此可用来进行高光谱图像的解混。由于它能够保证非负性并且自动调整迭代步长，近年来在高光谱解混中受到了关注。但是它最大的缺点是其目标函数具有明显的非凸性，因此存在大量的局部极小，如直接应用到高光谱图像，几乎无法得到最优解。

下面介绍与本发明相关的一些概念：

1.线性光谱混合模型

线性混合模型是高光谱遥感图像分析中最常用也最简便的一种数据模型，它假设端元之间没有互相影响。在线性混合模型中，像素的观察值等同于各端元的光谱特征按照它们的丰度进行线性组合。如某一个像素的观察值γ∈i^L×1可被表示为

γ＝Ms+e.                           (1)

其中M＝[m₁，m₂，...，m_P]∈i^L×P被称为端元光谱矩阵，当中的m_j对应着第j个端元的光谱特征，L是波段总数，P是端元个数；s＝(s₁，s₂，...，s_P)^T∈i^P×1为丰度矢量，表示在该像素中各端元所占的比例；e是误差矢量。显然，一个像元必然是由各端元中的若干种组成，且所有端元所占比例之和必须为1。因此丰度必须满足和为一约束(Abundance Sum-to-oneConstraint，ASC)和非负约束(Abundance Nonnegative Constraint，ANC)[4]

$\underset{j}{\mathrm{\Σ}}{s}_j=1,s_j>0.---\left(2\right)$

如果将一幅遥感图像中的每个像素都写成(1)式的形式，并依次排列，可以得到

R＝MS+E                            (3)

这里，R∈i^L×N为高光谱数据矩阵，其中每一行是一个波段，每波段有N个像素。M的定义不变，而S∈i^P×N为丰度矩阵，其中每一行是一个端元的丰度。

2.NMF算法

给定一个非负矩阵V∈i^n×m和一个正整数r＜min(m，n)，NMF的目标是找到两个非负矩阵W∈i^n×r和H∈i^r×m，使其满足

V≈WH.                             (4)

该问题通常由迭代方法解决，最常用的目标函数是欧式距离

$J(W,H)={\left|\right|V-\mathrm{WH}\left|\right|}^2=\underset{\mathrm{ij}}{\mathrm{\Σ}}{(V_{\mathrm{ij}}-{\left(\mathrm{WH}\right)}_{\mathrm{ij}})}^2.---\left(5\right)$

通过迭代使该目标函数达到最小，就可以得到NMF的解.

1999年，Lee和Seung提出了一种进行NMF的乘法迭代方法[3]。该方法通过在传统的梯度下降算法中巧妙地指定变化的步长推导得出。由于在迭代过程中只有加、乘和除操作，只要初始值非负，就能保证结果的非负。首先对式(5)求梯度，有

$\frac{\∂J(W,H)}{\∂W}=-{\mathrm{VH}}^T+{\mathrm{WHH}}^T,---\left(6\right)$

$\frac{\∂J(W,H)}{\∂H}=-W^{T}V+W^T\mathrm{WH}.---\left(7\right)$

那么梯度下降算法的迭代公式为

W←W+η_W.*(VH^T-WHH^T)，           (8)

H←H+η_H.*(W^TV-W^TWH).            (9)

以上的迭代公式并不能保证结果的非负性。如果定义迭代步长[5]

η_W＝W./(WHH^T)，                 (10)

η_H＝H./(W^TWH).                  (11)

将(10)和(11)代入(8)和(9)，得到乘法迭代公式：

W←W.*(VH^T)./(WHH^T)，            (12)

H←H.*(W^TV)./(W^TWH).             (13)

以上式子中的″.*″和″./″表示矩阵中的各元素分别相乘和相除。Lee和Seung已证明[5]，在这个迭代公式的作用之下，J(W，H)＝||V-WH||²是单调非增的。根据(3)式，高光谱解混可以看作一个非负矩阵分解问题，因此NMF方法可以被直接应用于高光谱图像，同时还可以满足线性模型所要求的ANC约束。但是NMF的目标函数具有明显的非凸性，因而存在大量局部极小，从而几乎无法得到最优解。

3.K-L散度和SID

K-L散度[6]是一种在信息论中得到广泛应用的信号相似度度量。若两个离散随机信号的概率分布函数分别为P(i)和Q(i)，那么定义Q相对于P的K-L散度为

$D_{\mathrm{KL}}\left(P\right|\left|Q\right)=\underset{i}{\mathrm{\Σ}}P\left(i\right)\log \frac{P\left(i\right)}{Q\left(i\right)}.---\left(14\right)$

式(14)的物理意义是用Q中的元素来表示P中的所有元素所需要的额外信息量的大小。如果这里的log是以2为底的对数，结果单位为bit，如果log是自然对数，结果单位为nat。K-L散度越大，说明完全用Q中的元素来表示P的难度越大，也即两个信号之间的相似度越低。因此，K-L散度本质上是一种非对称的概率意义上的距离。

在K-L散度的基础上，Chang于1999年提出了一种衡量光谱之间相关信息的尺度SID(Spectral Information Divergence)[7]，它被定义为

SID(x，y)＝D(x||y)+D(y||x).               (15)

如果令

$p=\frac{1}{\underset{j}{\mathrm{\Σ}}x\left(j\right)}x,$ $q=\frac{1}{\underset{j}{\mathrm{\Σ}}y\left(j\right)}y,---\left(16\right)$

它们是归一化的x和y，分别表示每个通道所占的比重，或者理解为光谱反射率在通道之间分布的概率函数。根据K-L散度的定义，可以得到

D(x||y)＝D_KL(p||q).                        (17)

可以看出SID实际上是一种对称化的K-L散度，很明显它始终是正数。

与本发明相关的现有技术有如下参考文献：

[1]C-I Chang，Hyperspectral Imaging：Techniques for Spectral Detection and Classification.New York：Plenum，2003.

[2]J.Li and J.M.Bioucas-Dias，“Minimum Volume Simplex Analysis：A Fast Algorithm toUnmix Hyperspectral  Data，”in：IEEE Geoscience and Remote Sensing Symposium(IGARSS’08)，Boston，MA，vol.3，July，2008，pp.250-253.

[3]D.D.Lee and H.S.Seung，“Learning the parts of objects by non-negative matrixfactorization，”Nature，vol.401，Oct.1999，pp.788-791.

[4]N.Keshava，“A survey of spectral unmixing algorithms，”Lincoln Lab.J.，vol.14，no.1，2003，pp.55-73.

[5]D.D.Lee and H.S.Seung，“Algorithms for non-negative matrix factorization，”AdvancedNeural Information Processing Systems，vol.13，2000，pp.556-562.

[6]S.Kullback and R.A.Leibler，“On information and sufficiency，”The Annals ofMathematical Statistics，vol.22，1951，pp.79-86.

[7]C-I Chang，“Spectral information divergence for hyperspectral image analysis，”in：IEEEGeoscience and Remote Sensing Symposium(IGARSS’99)，Hamburg，Germany，vol.1，1999，pp.509-511.

[8]D.C.Heinz and C-I Chang，“ Fully constrained least squares linear spectral mixtureanalysis method for material quantification

[9]J.Nascimento and J.Bioucas-Dias，“Vertex Component Analysis：A Fast Algorithm toUnmix Hyperspectral Data，”IEEE Transactions on Geoscience and Remote Sensing，vol.43，no.4，April，2002，pp.898-910.

[10]S.Jia and Y.Qian，“Constrained Nonnegative Matrix Factorization for HyperspectralUnmixing，”IEEE Transactions on Geoscience and Remote Sensing，vol.47，no.1，Jan.，2009，pp.161-173.

[11]L.Miao and H.Qi，“Endmember extraction from highly mixed data using minimumvolume constrained nonnegative matrix factorization，”IEEE Transactions on Geoscienceand Remote Sensing，vol.45，no.3，Mar.2007，pp.765-777.

[12]D.Landgrebe，Multispectral data analysis：A signal theory perspective.West Lafayette：School of Electrical & Computer Engineering，Purdue University，1998，pp.56-89.

[13]R.N.Clark and G.A.Swayze，“Evolution in Imaging Spectroscopy Analysis and SensorSignal-to-Noise：An Examination of How Far We Have Come，”The 6th Annual JPLAirborne Earth Science Workshop，Mar.，1996.Available：http://speclab.cr.usgs.gov/PAPERS.imspec.evol/aviris.evolution.html

发明内容

本发明的目的是提供一种新的遥感图像处理方法，尤其涉及基于有约束非负矩阵分解的高光谱遥感图像混合像元分解方法。

本发明根据高光谱图像光谱和丰度的特点，在非负矩阵分解算法的目标函数中引入丰度分离性和平滑性的约束条件，在迭代过程中监测目标函数的变化情况，并且在适当的时机移除这些约束条件并继续迭代，从而克服NMF算法容易陷入局部极小的缺点，可以有效的进行高混合遥感数据的混合像元分解。

通过对实际高光谱图像特点的分析，本发明在乘法NMF算法中加入了丰度的分离性和平滑性两个约束条件。所述的约束条件符合高光谱图像的物理实际特性。

由于地物一般都是成片成块分布的，而不会布满整个图像的区域，因此在实际的丰度图上，各种地物分别有其自己的主导区域，它们之间的相关性应该较小。本发明将这样的约束条件称为分离性(Separation)约束。

分离性约束考虑的是不同端元之间相同位置像元的相关关系，但是没有考虑不同像元间的空间关系。平滑性(Smoothness)的加入正是体现了数据的空间关系，这是由于实际地物的分布往往是有规律且在大部分地方都保持连贯性和均匀性，突变只会存在于少数的地物边缘。本发明中平滑性约束的特点在于能够反映更大范围的局部空间关系，并且计算简便。

初始值的选取也是影响结果的重要因素。对于高光谱图像解混，目前通常采用随机初始化的方法。本发明根据随机初始化的弊端，改进了算法的初始化方法。此外，本发明还提出了一种新的算法停止条件，进一步改进了算法的性能。

本发明提出的遥感图像混合像元分解方法，其包括：

1.分离性约束

本发明对K-L散度进行了改进，得到了一种新的函数来衡量相关信息，使其更适合加入NMF中进行迭代，并且能够保证算法的稳定。由于物理含义的改变，该函数已经不能被称为散度。本发明称它为分离性函数。

沿用以上的符号定义，两个矢量x和y的分离性函数被定义为

$\mathrm{Separation}(x,y)=\underset{i}{\mathrm{\Σ}}[p\left(i\right)f\left(\frac{p\left(i\right)}{q\left(i\right)}\right)+q\left(i\right)f\left(\frac{q\left(i\right)}{p\left(i\right)}\right)].---\left(18\right)$

它实质上是用f(x)替换K-L散度中的lnx，f(x)的选取须满足如下几个原则：首先，f(x)在(0，+∞)之间有定义，单调增，f(1)＝0且f(x)+f(1/x)≥0。为了能够体现分离性的物理意义，前三个条件是显然的。而f(x)+f(1/x)≥0可以保证分离性函数始终为正数。其次，f(x)和f′(x)都是有界的，且在x→+∞或x→0时，f′(x)→0。这一点有助于权重的选择，并且可以使算法更稳定，不容易发散。综合以上两点，本发明选择以下的函数：

$f\left(x\right)=1-2^{1-x^2}---\left(19\right)$

它的曲线如图1所示。同样地，这样得到的分离性函数始终为正数，函数值越大说明信号之间的相关性越小。而分离性为0的两个信号则是完全线性相关的。

考虑所解出的一组端元的丰度，在它们之间两两求分离性，并将所有结果的求和，作为整个丰度组的分离性的度量。本发明定义

$J_1\left(S\right)=\frac{1}{P^2}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\underset{j=1}{\overset{P}{\mathrm{\Σ}}}\mathrm{Seperation}(s_i,s_j)---\left(20\right)$

作为分离性约束的目标函数。其中s_i和s_j分别是第i个和第j个端元的丰度矢量。

将该目标函数对S矩阵求导数，可以得到其对S中每个元素的导数为

$\frac{\∂J_1\left(S\right)}{\∂S_{\mathrm{pn}}}=\frac{4\ln 2}{c_p}\left\{\underset{j}{\mathrm{\Σ}}\underset{k}{\mathrm{\Σ}}{Q}_{\mathrm{pk}}\right[\frac{Q_{\mathrm{jk}}^3}{Q_{\mathrm{pk}}^3}{2}^{-\frac{Q_{\mathrm{jk}}^2}{Q_{\mathrm{pk}}^2}}+(\frac{1}{2\ln 2}-\frac{Q_{\mathrm{pk}}^2}{Q_{\mathrm{jk}}^2})2^{-\frac{Q_{\mathrm{pk}}^3}{Q_{\mathrm{jk}}^3}}-\frac{1}{4\ln 2}]---\left(21\right)$

$-\underset{j}{\mathrm{\Σ}}[\frac{Q_{\mathrm{jn}}^3}{Q_{\mathrm{pn}}^3}{2}^{-\frac{Q_{\mathrm{jn}}^2}{Q_{\mathrm{pn}}^2}}+(\frac{1}{2\ln 2}-\frac{Q_{\mathrm{pn}}^2}{Q_{\mathrm{jn}}^2})2^{-\frac{Q_{\mathrm{pn}}^3}{Q_{\mathrm{jn}}^3}}-\frac{1}{4\ln 2}]\}.$

其中

$c_p=\underset{k}{\mathrm{\Σ}}{S}_{\mathrm{pk}},$ $Q_{\mathrm{jk}}=\frac{S_{\mathrm{jk}}}{c_j}.---\left(22\right)$

也就是说，c_p是第p个端元的丰度之和，而Q矩阵即为归一化之后的S矩阵。

2.平滑性约束

分离性约束仅仅是将丰度图中的每个像素点分别独立地考虑，并没有涉及到像素之间的关系。而像素之间的关系事实上表现为空间关系。由于地物总是在缓慢变化的，丰度的变换应该也是平滑的，因此必须引入丰度平滑的约束。为了得到临近像素之间的空间关系，本发明中，暂时改变数据的表示方法。将每个端元的丰度，也即式(3)中S的一行，按照其像素的空间关系整理成对应的一个矩阵。由于丰度平滑约束是对每个端元分别考虑的，本发明只需讨论其中任意一个端元的丰度，设其被整理为S∈i^r×c，其中r和c分别表示图像的行数和列数，并有r×c＝N。

对于矩阵S中的某一个像素S _ij而言，可以用它的值和周围像素的差异来表征该像素附近的平滑程度，考虑某个参考像素S _ij，本发明将除了它本身之外的所有像素划分为八个区域G_k(k＝1，2，...，8)。图2(a)给出了划分的方法，各区域的标号也均被标明。

令S _ijk ^％表示第k个区域的所有像素的加权平均值，也即

${\underset{\‾}{S}}_{\mathrm{ijk}}^{\%}=\underset{m,n,{\underset{\‾}{S}}_{\mathrm{mn}}\∈G_k}{\mathrm{\Σ}}{w}_{\mathrm{mn}}{\underset{\‾}{S}}_{\mathrm{mn}}---\left(23\right)$

其中w_mn为指定位置的权值，可定义参考像素S_ij附近的平滑性函数为

$J_2\left({\underset{\‾}{S}}_{\mathrm{ij}}\right)=\frac{1}{2}\×\frac{1}{8}\underset{k=1}{\overset{8}{\mathrm{\Σ}}}{({\underset{\‾}{S}}_{\mathrm{ijk}}^{\%}-{\underset{\‾}{S}}_{\mathrm{ij}})}^2,---\left(24\right)$

其中的1/8表示在一般情况下，各区域在平滑性中所占的权重相等，均为1/8。

将(24)式推广到所有的像素并求和，得到整个图像区域的平滑性函数

$J_2\left(\underset{\‾}{S}\right)=\frac{1}{2}\×\frac{1}{8}\underset{k=1}{\overset{8}{\mathrm{\Σ}}}{\left|\right|{\underset{\‾}{S}}_k^{\%}-\underset{\‾}{S}\left|\right|}^2,---\left(25\right)$

其中S的定义如本发明上文所述，而S _k ^％是由S _ijk ^％放在对应参考像素的位置所得到的矩阵。

的求解方法由权值w_mn分配的方式决定。分配权值应遵循的原则是，越靠近参考像素的点给予越大的权值。此外，每个区域的所有权值之和应为1。本发明采用图2(b)的方式来进行分配，将最靠近参考像素的位置的权值定为1/2，此后每向外延伸一层，权值减半。对于离参考像素很远的位置，其权值很小，这些值没有实际意义，因此，本发明实际上只采用了参考像素外围最多5层的权值，也即将小于1/32的值指定为0。

根据以上的说明，如果令

$P_1=\left[\begin{array}{ccccc}1& 0& 0& L& 0\\ 1& 0& 0& L& 0\\ \frac{1}{2^2}& \frac{1}{2}& 0& L& 0\\ M& M& O& O& M\\ \frac{1}{2^{N-1}}& \frac{1}{2^{N-2}}& L& \frac{1}{2}& 0\end{array}\right],$ $P_4=\left[\begin{array}{ccccc}0& 0& 0& L& 0\\ \frac{1}{2}& 0& 0& L& 0\\ M& O& O& O& M\\ \frac{1}{2^{N-2}}& L& \frac{1}{2}& 0& 0\\ \frac{1}{2^{N-1}}& L& \frac{1}{2^2}& 1& 1\end{array}\right],---\left(26\right)$

$P_2=\left[\begin{array}{ccccc}0& \frac{1}{2}& \frac{1}{2^2}& L& \frac{1}{2^{N-1}}\\ 0& 0& O& O& M\\ M& M& O& \frac{1}{2}& \frac{1}{2^2}\\ 0& 0& L& 0& 1\\ 0& 0& L& 0& 1\end{array}\right],$ $P_3=\left[\begin{array}{ccccc}1& 1& \frac{1}{2^2}& L& \frac{1}{2^{N-1}}\\ 0& 0& \frac{1}{2}& O& M\\ M& M& O& O& \frac{1}{2^2}\\ 0& 0& L& 0& \frac{1}{2}\\ 0& 0& L& 0& 0\end{array}\right],---\left(27\right)$

则有

${\underset{\‾}{S}}_1^{\%}=P_1\underset{\‾}{S},$ ${\underset{\‾}{S}}_2^{\%}=P_2\underset{\‾}{S},$ ${\underset{\‾}{S}}_3^{\%}=\underset{\‾}{S}{P}_3,$ ${\underset{\‾}{S}}_4^{\%}=\underset{\‾}{S}{P}_4.---\left(28\right)$

如果再令

$P_5=\left[\begin{array}{ccccc}0& 0& 0& L& 0\\ \frac{1}{2}& 0& 0& L& 0\\ \frac{1}{2^2}& \frac{1}{2}& 0& L& 0\\ M& M& O& O& M\\ \frac{1}{2^{N-1}}& \frac{1}{2^{N-2}}& L& \frac{1}{2}& 0\end{array}\right],$ $P_7=\left[\begin{array}{ccccc}0& \frac{1}{2}& \frac{1}{2^2}& L& \frac{1}{2^{N-1}}\\ 0& 0& \frac{1}{2}& L& \frac{1}{2^{N-2}}\\ 0& 0& 0& O& M\\ M& M& O& O& \frac{1}{2}\\ 0& 0& L& 0& 0\end{array}\right],---\left(29\right)$

且P₆＝P₅，P₈＝P₇，T₁与P₅具有相同的形式，T₂与P₇具有相同的形式，则有

${\underset{\‾}{S}}_5^{\%}=P_5\underset{\‾}{S}{T}_2,$ ${\underset{\‾}{S}}_6^{\%}=P_6\underset{\‾}{S}{T}_1,$ ${\underset{\‾}{S}}_7^{\%}=P_7\underset{\‾}{S}{T}_2,$ ${\underset{\‾}{S}}_8^{\%}=P_8\underset{\‾}{S}{T}_1,---\left(30\right)$

以上本发明定义了完整的平滑性函数。对这个约束条件求解其梯度，如果令

P₁₂＝(I-P₁)^T(I-P₁)+(I-P₂)^T(I-P₂)

P₃₄＝(I-P₃)^T(I-P₃)+(I-P₄)^T(I-P₄)

P＝P₅+P₇＝P₆+P₈

T＝T₁+T₂                                   (31)

$A=P_5^T{P}_5+P_7^T{P}_7=P_6^T{P}_6+P_8^T{P}_8$

$B=T_1{T}_1^T+T_2{T}_2^T$

则有

$\frac{\∂J_2\left(\underset{\‾}{S}\right)}{\∂\underset{\‾}{S}}=\frac{1}{8}(P_{12}\underset{\‾}{S}+\underset{\‾}{S}{P}_{34}+4\underset{\‾}{S}+A\underset{\‾}{S}B-2P\underset{\‾}{S}Q).---\left(32\right)$

其中所有符号同本发明文所定义。虽然平滑性目标函数比较复杂，但是在化简后其梯度的计算量很小，(31)式中的所有矩阵都可在算法在初始化时计算好，因此实际迭代过程中每步只需要计算式(32)。和PSNMFSC[10]相比，本发明的平滑性约束可以明显地降低运算复杂度。

由于在实际算法中，是对S而非S进行迭代，因此在求出每个端元的梯度之后，必须再将它们分别重新整理成矢量，再合并作为梯度矩阵

3.ASSNMF迭代公式

在线性混合模型中，任意像素在所有端元上的丰度之和必须为1。这一点必须也作为约束条件被加入到算法中。这里本发明直接采用了文献[8]中的方法。

在算法初始化的时候进行下面的扩展：

$M_A\←\left[\begin{array}{c}M\\ \mathrm{\δ}{1}^T\end{array}\right],$ $R_A\←\left[\begin{array}{c}R\\ {\mathrm{\δ}1}^T\end{array}\right].---\left(33\right)$

其中1^T是元素全为1的矢量，而δ控制着ASC影响的大小。用M_A和R_A替代R和M进行算法的迭代，并在每次迭代完M之后，从M_A和R_A中取出原本属于M和R的部分，并重新进行一次式(33)的扩展即可。这样在算法收敛的同时也能满足ASC。

将所有的约束条件都加到NMF算法中去，即可得到最后的迭代公式。所述的ASSNMF的问题可被归纳为：

最小化并满足以下条件：M≥0，S≥0。

其中λ₁和λ₂分别是分离性和平滑性约束条件的权重，而其他符号已经作过说明。分别用R_A和M_A取代R和M，可以得到梯度下降算法的迭代公式：

M_A←M_A+η_M.*(R_AS^T-M_ASS^T)，                      (34)

$S\←S+{\mathrm{\η}}_S.*({M_A}^T{R}_A-{M_A}^T{M}_{A}S+{\mathrm{\λ}}_1\frac{\∂J_1\left(S\right)}{\∂S}-{\mathrm{\λ}}_2\frac{\∂J_2\left(S\right)}{\∂S}).---\left(35\right)$

与(10)和(11)类似地定义迭代步长

η_M＝M_A./(M_ASS^T)，                              (36)

η_S＝S./(M_A ^TM_AS).                               (37)

将(36)和(37)代入(34)和(35)之后得到乘法迭代公式：

M_A←M_A.*(R_AS^T)./(M_ASS^T)，               (38)

$S\←S.*({M_A}^T{R}_A+{\mathrm{\λ}}_1\frac{\∂J_1\left(S\right)}{\∂S}-{\mathrm{\λ}}_2\frac{\∂J_2\left(S\right)}{\∂S})./\left({M_A}^T{M}_{A}S\right).---\left(39\right)$

4.约束条件的移除

本发明中，所述的几个约束条件处于互相牵制的状态，其中分离性和平滑性作为防止局部极小的手段，在算法检测到分离性目标函数收敛之后，移除分离性和平滑性约束，从而改善欧式距离以及ASC约束的收敛情况，同时也可降低算法复杂度，加快运算速度。

本发明中，用每次迭代与前一次迭代之间J₁(S)的相对变化

$\mathrm{chg}=\frac{J_1{\left(S\right)}^t-J_1{\left(S\right)}^{t-1}}{J_1{\left(S\right)}^{t-1}}---\left(40\right)$

作为判断的标准，当连续N次迭代后(本发明中，取N＝20)该值均小于某个阈值η时，可以认为分离性目标函数已经收敛，此时将λ₁和λ₂置为0，移除分离性和平滑性约束，继续进行只有ASC约束的NMF，直至整体算法的收敛。

5.算法的初始化

由于初始值能够影响算法最终能否收敛到最优解以及所需要的迭代时间，因此进行有效的初始化是整个算法中至关重要的一步。

端元矩阵M的初始化一般是采用从图像中随机选取指定个数的像素点，认为它们是纯像元并将它们的值作为端元的方法，但是如果有几个随机点落到同一个很小的区域内，就会出现两个或更多端元的初始值完全一样的极端情况，实验表明，这种情况很容易让算法陷入局部极小，因此本发明对随机方法作了改进，当每产生一个新的随机像元，就通过式(15)所定义的SID比较它与已经产生的端元之间的差异，当这个差异大于指定阈值时，才将该像元作为初始的端元，否则就重新寻找另一个随机像元。

在得到M后，采用下式进行S的初始化：

S＝(M^TM)^-1M^TR.                      (41)

根据NMF的特点，初始值必须非负才能保证结果非负，否则在迭代过程中可能发散。因此在产生S之后，还需要将它归一化到[0，1]之间，才能开始进行迭代。

本发明经实验发现，对结果影响较大的参数只有λ₁和λ₂两个，分别是丰度分离性和平滑性的约束。实验结果表明，在一定的范围内，算法的结果只和这两个参数之比λ₁/λ₂有关，因而事实上只需要调节其中一个。该比值在一定范围内都能得到正确的解。

本发明具有如下优点：

本发明提供了一种新的基于非负矩阵分解(Nonnegative Matrix Factorization，NMF)算法的混合像元分解方法。所述方法根据高光谱图像光谱和丰度的特点，在非负矩阵分解算法的目标函数中引入丰度分离性和平滑性的约束条件，并且在适当的时机移除这些约束条件并继续迭代，从而克服了NMF算法容易陷入局部极小的缺点，可以有效的解决高混合遥感数据的混合像元分解方法。其中的分离性约束是在K-L散度的基础上改进而来，可以保证算法的稳定。而平滑性约束可以对整个图像范围进行全局考虑，并且计算简便。本方法在基于多光谱和高光谱遥感图像的高精度的地物分类以及地面目标的检测和识别方面具有特别重要的应用价值。

附图说明：

图1分离性函数中的f(x)。

图2丰度平滑性示意.(a)像素周围区域的划分，(b)权值的分配方式。

图3仿真数据示意.(a)丰度，(b)光谱，(c)信噪比为20dB时第185波段的数据。

图4算法的抗噪声性能.(a)RMSE，(b)SAD。

图5纯像元缺失程度对算法的影响.(a)RMSE，(b)SAD.。

图6像元个数对算法的影响.(a)RMSE，(b)SAD。

图7移除约束条件前后的收敛曲线对比.(a)欧式距离，(b)ASC，(c)RMSE，(d)SAD。

图8 Indiana数据的伪彩色图。

图9 Indiana数据的分解结果.(a)人工建筑，(b)小麦，(c)玉米，(d)大豆，(e)植被，(f)干草堆。

图10 Urban数据的伪彩色图。

图11 Urban数据的丰度分解结果.(a)屋顶1，(b)树木，(c)水泥路面，(d)屋顶2/阴影，(e)草地，(f)沥青路面/停车场。

图12 Cuprite数据的伪彩色图。

图13 Cuprite数据的丰度分解结果，(a)明矾石Alunite，(b)榍石Sphene，(c)高岭石Kaolinite#1，(d)玉髓Chalcedony，(e)铵长石Buddingtonite，(f)皂石Nortronite，(g)高岭石Kaolinite#2，(h)白云母Muscovite，(i)黄甲铁石Jarosite，(j)高岭石Kaolinite#3，(k)蒙脱石Montmorillonite，(l)沙漠地表Desert Vanish。

具体实施方式

下面，分别用仿真数据和实际遥感图像数据为例进一步描述本发明。

实施例1

1.仿真数据

人工产生的仿真数据被用于测试算法的性能。本发明提出的ASSNMF与以下几种高光谱图像的解混算法进行比较：VCA[9]、带有光谱和丰度的平滑性以及丰度稀疏性约束的PSNMFSC[10]以及带有最小体积约束的MVCNMF[11]。其中VCA只能得到光谱矩阵，其他方法都可以直接从数据中解出光谱和丰度。对于VCA，本发明在其解出的光谱的基础上使用FCLS算法[8]，得到对应的丰度，该方法记为VCA-FCLS。

本发明使用光谱角距离(Spectral Angel Distance，SAD)和均方根误差(Root MeanSquare Error，RMSE)两个指标衡量解混结果的优劣。该两个指标分别被用来衡量光谱与丰度解混结果与参考值的近似程度。对于第k个端元而言，其SAD被定义为[9]

${\mathrm{SAD}}_k={\cos }^{-1}\frac{\⟨{\hat{m}}_k,m_k\⟩}{\left|\right|{\hat{m}}_k\left|\right|\left|\right|m_k\left|\right|}---\left(42\right)$

其中m_k和

分别表示该端元的光谱解混结果和参考光谱。而丰度与参考值之间的RMSE被定义为[10]

${\mathrm{RMSE}}_k=\sqrt{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(S_{\mathrm{kj}}-{\hat{S}}_{\mathrm{kj}})}^2}---\left(43\right)$

其中S_kj和分别表示该端元在某一指定像元的丰度解混结果和参考值。在进行仿真数据的测试时，本发明将所有端元的平均SAD以及平均RMSE作为比较标准。

在仿真数据的产生上，本发明遵循的原则是使它的丰度特性尽可能接近真实数据，具有成片分布、连续变化、有一定轮廓的特点。而光谱，本发明采用了美国地质调查局(USGS)光谱库中的矿物光谱，这样产生的仿真数据和真实的遥感图像有一定的相似之处，在评价算法的性能时更具代表性。

本发明中，丰度产生方法是，首先在一个100×100的正方形区域(第三个实验除外)的四个角上各放置一种端元，然后在在正方形区域的中间再放置第5种端元。它们的丰度值都是从中心往外随着距离的增加而递减直至为0。最后对所有丰度之和进行归一化，利用和为一约束，将不同端元混合到一起。根据所指定的各端元中心的丰度值的不同，可以产生纯像元缺失程度不同的图像。图3(a)中给出了纯像元没有缺失的情况。从USGS矿物光谱库中选取5种224波段的光谱，如图3(b)所示，对应的地物根据图例从1到5分别是明矾石(Alunite_GDS83)、铵长石(Buddingtonite_NHB231)、石榴石(Pyrope)、皂石(Nontronite_GDS41)和沙漠地表(Desert_Vanish_GDS78A)，分别作为各端元的反射特性，与丰度矩阵相乘并加上不同程度的高斯白噪声，就得到了实验用的仿真数据。其大小为224×10000(第三个实验除外)。作为本发明的一个实例，图3(c)中给出了由图3(a)中的丰度形成的仿真数据加上噪声后，第185波段的图像，其信噪比为20dB。

本发明通过四个实验，说明所提出方法的有效性，并与其他类似方法进行比较。在第一个实验中，改变往图像中所加的噪声等级，评判算法的抗噪声性能。第二个实验通过改变丰度数据的纯像元缺失程度，研究算法对于无纯像元数据的适应性。第三个实验研究算法性能和像元个数的关系。第四个实验是验证去掉分离性和平滑性约束项之后继续迭代的方法是否能够对分解结果起到改善作用。在上述实验中，本发明算法ASSNMF的参数均设置为λ₁＝0.28×N，λ₂＝5/P，δ＝0.02×L。PSNMFSC算法的参数按照文献[10]的标准，设置为α＝0.5，β＝0.1，γ_m＝0.01，γ_s＝0.5。此外，对于所有的算法，均认为端元个数是已知的。

实验1测试算法的抗噪声性能：改变往仿真数据中所加的噪声，以此来比较各算法的抗噪声性能。根据信噪比的不同，从∞(无噪声)、35dB∶10dB，每隔5dB为一档。

实验显示，本发明的算法和VCA-FCLS解出了正确的丰度轮廓与光谱，而MVCNMF和PSNMFSC则没有解出正确的结果，光谱和丰度都与参考值相差甚远。图4(a)和图4(b)分别给出了几种算法的RMSE指标和SAD指标随着信噪比变化的情况。结果显示，随着信噪比的降低，所有算法的效果在整体上都逐渐变差。在所有方法中，本发明的方法综合表现最好，当SNR＝∞，也即没有噪声的时候，VCA-FCLS给出了最好的结果，两项评价指标均为零。事实上由于VCA正确地寻找到纯像元并将其作为端元，该结果与标准值完全相同。然而随着噪声的增加，VCA寻找的端元逐渐开始偏移准确值，SAD和RMSE逐渐增加，此时本发明相对于VCA-FCLS的优势开始凸显。MVCNMF和PSNMFSC的两项指标都始终处于较高的水平。此外值得注意的是，由于所有这些算法都不能去除噪声，当信噪比降到15dB以下时，由于噪声方差过大，所有方法的RMSE都不理想，但是此时本发明的算法解出的SAD是最准确的，并且求解出丰度轮廓上与标准丰度(图3(a)所示)最为接近。

实验2测试算法对无纯像元情况的适应性：改变仿真数据的纯像元缺失程度，以此来比较各算法对无纯像元情况的适应性。本发明用所有端元中出现的最高的丰度值定量表示纯像元缺失的程度，从1(有纯像元)∶0.5，每隔0.1为一档，而信噪比被固定在30dB。实验显示，只有本发明的算法和VCA-FCLS解出了正确的丰度轮廓与光谱，而其他方法都未能解出正确的结果。图5(a)和图5(b)分别给出了在不同的纯像元缺失程度下，各算法的RMSE和SAD的变化情况。实验显示，随着像元纯度的降低，本发明的方法和VCA-FCLS的效果都在逐渐变差，而其他几种方法由于没有解出正确的结果，故其评价参数变化不大。本发明的算法在RMSE和SAD上优于VCA-FCLS。值得注意的是，当像元纯度降低到0.6时，所有的端元已经混合得非常严重，因而VCA-FCLS和本发明的方法在SAD上相继出现了一个大幅的变化。

实验3研究算法性能和像元个数的关系：改变仿真数据的大小，研究当像元数改变时算法性能的变化。实验的范围是从20×20到100×100，纯度被固定在0.9，信噪比被固定在30dB。

图6(a)和图6(b)分别给出了在像元数不同时，各算法的RMSE和SAD的变化情况。可以看出随着像元数的增加，所有算法的性能都有小幅度的改善，对于VCA而言，这是因为更多的像元能形成更为完整的单形体，使顶点的寻找更为精确。而对于其他算法而言，这是因为能从观察矩阵中得到更丰富的信息。另外，与前两个实验一样，只有本发明的算法和VCA-FCLS解出了正确的丰度轮廓与光谱，本发明算法给出了最好的结果。

实验4研究移除约束条件对结果的影响：为了验证移除分离性和平滑性约束的方法的有效性，在每步迭代完后，检查当前的欧式距离和ASC的收敛情况，并根据参考的标准丰度和光谱，计算出每一步的结果与标准值之间的RMSE和SAD两个参数，绘制成收敛曲线。此外在分离性目标函数收敛的时候，并不立即移除约束条件，而是设定让算法继续迭代一段时间后再移除，以更清楚地说明移除约束条件是否会对结果起到改善作用。本实验中使用有纯像元且无噪声的数据。实验结果如图7所示，所设定的移除约束条件的时机是第5000次迭代完成后。结果显示，在第5000次迭代之前，所有参数都已经收敛，但是移除约束条件后，它们开始继续下降，直至收敛到更低的值。

实施例2 实际数据实验

分别使用三个实际的高光谱遥感图像数据集对所提出算法的性能进行测试。

第一个数据集是由机载可见光及红外成像光谱仪(Airborne Visible/InfraredImaging Spectrometer，AVIRIS)拍摄的Indiana数据。它成像于1992年6月，成像区域为美国印第安那州的派恩遥感测试点，该数据有220个波段，波长范围从0.4：2.5μm，光谱分辨率为10nm，空间分辨率为17m。实验所用的图像大小为145×145。该数据已被广泛地用于遥感图像的混合像元分解算法的研究比较。Purdue大学已经给出一份关于该区域的实地调查报告[12]。该地区是位于印第安纳州西部城市西拉斐特(West Lafayette，IN)西北方向约10km处的一片农田，覆盖该区域的主要是各种农作物(大约占三分之二，包括玉米、小麦、大豆、干草堆)和天然植被(大约占三分之一，由树林、草地等组成)。除此之外，还有一些人工用地，包括区域顶部的一条双向高速公路(U.S.52和U.S.231)和一条铁路、中间的一条2级公路(杰克逊高速公路)、区域上方的一个无线电发射塔以及一些零星的房屋。取第70、86、136波段分别作为R、G、B分量合成伪彩色图，如图8所示。

在进行分解之前，该数据的第1∶4，78∶82，103∶115，148∶166以及211∶220波段由于信噪比太低或为水吸收波段而被移除，剩下169个波段被用于进一步处理。为了定量评价算法的性能，本发明在分解之前参考地物真实分布情况[12]对图像中的端元进行了人工提取，以此来得到各端元的参考光谱。在本图像中共有72个较纯的像元被选取出来，其中对应玉米、小麦、植被、人工建筑、干草和大豆的像元数分别为5，10，16，13，15和13。本实验求取每一类像元观察值的均值，作为该端元的参考光谱。

采用本发明的方法对Indiana数据的丰度解混结果如图9所示，其中各参数的选取均与仿真数据相同。通过比较结果显示，丰度解混结果与实地调查结果[12]非常吻合。另外，本发明将解出的光谱与手动选取的参考光谱比较其光谱角，以定量衡量算法的性能。与上述仿真数据一样，本发明同时比较了ASSNMF，VCA，PSNMFSC和MVCNMF四种方法的解混结果，结果显示，ASSNMF得到了最好的结果，其他的算法则误差较大。

表1是Indiana数据的光谱分解结果与参考值之间的光谱角。

表1

第二个数据集是由超光谱图像收集实验仪器(Hyperspectral Digital Imagery CollectionExperiment，HYDICE)拍摄于1995年10月的Urban数据，该数据共有210个波段，光谱分辨率为10nm，空间分辨率为2m，是在4430m的飞行高度上拍摄的。成像区域位于美国得克萨斯州胡德堡(Fort Hood，TX)附近的科珀勒斯科夫镇(Copperas Cove)，画面范围内的主要地物有位于顶部的一条高速公路(U.S.190)，公路旁的一家大型购物中心以及门前的停车场，一些小公路，整齐排列的小房屋，草坪和树木等。另外，由于拍摄时太阳角度较低，树木和房屋都在地面留下了阴影。

取第12、99和171分别作为R、G、B分量合成伪彩色图，如图10所示。从中可以看出，屋顶主要有两种不同的颜色，其中一种较亮，而另一种与地上的阴影很接近。公路有沥青路面和水泥路面两种，沥青路面反射率较低，水泥路面反射率较高，另外停车场也属于沥青路面。而草坪和树木的光谱特性也有所差别。因此可以认为，该数据中有6个端元，分别是沥青路面/停车场，屋顶1，屋顶2/阴影，草地，树木，水泥路面。

在进行分解之前，该数据的第1∶4，76，87，101∶111，136∶153以及198∶210波段由于信噪比太低或为水吸收波段而被移除，剩下162个波段被用于进一步处理。为了定量评价算法的性能，本发明在分解之前参考以上所描述的地物分布情况，对图像中的端元进行了人工提取，以此来得到各端元的参考光谱。在本图像中，每个端元均有15个较纯的像元被选取出来，共90个。本实验求取每一类像元观察值的均值，作为该端元的参考光谱。

使用本发明所提出的方法对Urban数据的丰度分解结果如图11所示。可以看出，分解结果很符合实际地物的分布情况。表2中比较了几种不同方法所得到的光谱解混结果与参考光谱之间的差异，用光谱角表示。很明显，除了屋顶1之外，本发明的方法均得到了令人满意的结果。另外在所有方法中，本发明的方法总体上表现出了最好的综合性能。

表2 Urban数据的光谱分解结果与参考值之间的光谱角

第三个数据集是由AVIRIS拍摄于1997年6月19日的美国内华达州Cuprite地区数据，该数据的波长范围是0.37～2.48μm，光谱分辨率为10nm，共有224个波段。该数据已被广泛应用于高光谱图像解混算法的评价中，并且Swayze和Clark等人已经给出了关于该地区的地物真实分布的报告[13]。该地区位于内华达州南部城市拉斯维加斯(LasVegas，NV)西北方向约200km的沙漠中，区域内有一条高速公路(U.S.95)，地表基本无植物覆盖，主要为裸露的矿物，并且给出了各种矿物的大致分布情况。另外，该地区各种矿物之间的混合现象较为普遍，很适合用来检验算法对高混合度数据的适应能力。

本发明从图像中截取了大小为250×191的一块进行实验，其伪彩色图如图12所示。在进行分解之前，第1～2，104～113，148～167以及221～224波段由于信噪比过低或者为水吸收波段而被舍弃掉，留下188个波段进行进一步的处理。根据已有的对该数据的研究错误！未找到引用源。，本发明设置端元数目为12，使用本发明提出的方法所得到的丰度解混结果如图13所示。与实地勘测地物分布图[13]相比较，可以确定出这些端元各自对应的矿物种类。为了进一步衡量算法性能，本实验将美国地质调查局(USGS)库中的对应矿物光谱作为参考光谱，并求取解混结果与它们之间的光谱角。几种方法的光谱角的定量对比如表3所示。可以看出，在几种方法中，本发明所提出的方法给出了最好的结果。

表3 Cuprite数据的光谱分解结果与参考值之间的光谱角

Claims (8)

1.一种遥感图像混合像元分解方法，其特征在于，在非负矩阵分解算法的目标函数中引入丰度分离性和平滑性的约束条件，在迭代过程中监测目标函数的变化情况，在适当的时机移除所述的约束条件并继续迭代，进行高混合遥感数据的混合像元分解。

2.按权利要求1所述的遥感图像混合像元分解方法，其特征在于，其包括步骤：

1)分离性约束

对K-L散度进行改进，得分离性函数，

其中，用f(x)替换K-L散度中的lnx，f(x)的选取满足如下原则：f(x)在(0，+∞)之间有定义，单调增，f(1)＝0且f(x)+f(1/x)≥0；

2)平滑性约束

将每个端元的丰度，按照其像素的空间关系整理成对应的一个矩阵；

3)ASSNMF迭代公式

将所有的约束条件都加到NMF算法中去，得到最后的迭代公式；

4)约束条件的移除

在算法检测到分离性目标函数收敛之后，移除分离性和平滑性约束；

5)算法的初始化

将产生的新的随机像元，与已经产生的端元比较其间的差异，当差异大于指定阈值时，将该新像元作为初始的端元。

3.按权利要求1所述的遥感图像混合像元分解方法，其特征在于，其中所述的分离性函数为：

4.按权利要求3所述的遥感图像混合像元分解方法，其特征在于，其中所述的分离性函数始终为正数。

5.按权利要求3所述的遥感图像混合像元分解方法，其特征在于，其中步骤1)中，将所有结果的求和，作为整个丰度组的分离性的度量，定义为：

$J_1\left(S\right)=\frac{1}{P^2}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\underset{j=1}{\overset{P}{\mathrm{\Σ}}}\mathrm{Seperation}(s_i,s_j)$

作为分离性约束的目标函数，其中s_i和s_j分别是第i个和第j个端元的丰度矢量。

6.按权利要求1所述的遥感图像混合像元分解方法，其特征在于，步骤2)中，将任意一个端元的丰度，设其被整理为S∈i^r×c，其中r和c分别表示图像的行数和列数，并有r×c＝N。

7.按权利要求1所述的遥感图像混合像元分解方法，其特征在于，步骤2)中，用矩阵中的某一个像素的值和周围像素的差异表征该像素附近的平滑程度。

8.按权利要求1所述的遥感图像混合像元分解方法，其特征在于，步骤2)中，将某个参考像素除其本身之外的所有像素划分为八个区域G_k(k＝1，2，...，8)，各区域在平滑性中所占的权重相等。

Images