CN101221243A - 基于非负矩阵因式分解的遥感图像混合像元分解方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，具体为一种基于非负矩阵因式分解的遥感图像混合像元分解方法。本发明将单形体方法作为非负矩阵分解算法的前端算法，即由单形体方法获得初步端元信号估计结果，然后由其组合成的矩阵作为端元信号矩阵的初始值，代入经修改后的非负矩阵因式分解迭代式进行更新运算，获得最终的分解结果。本发明克服了两种算法各自的缺点，有效的解决了高混合遥感数据的混合像元分解问题。本方法在基于多光谱和高光谱遥感图像的高精度地物分类以及地面目标的检测和识别方面具有重要应用价值。

Description

基于非负矩阵因式分解的遥感图像混合像元分解方法

技术领域

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

背景技术

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

单形体几何学方法[2][3]是目前遥感图像混合像元分解问题中十分重要的一类方法，该类方法物理意义明显，算法较为简单，获得了广泛的应用。然而这类方法均需要一个较强的前提假设，即每一端元在图像中均至少存在一个纯像元，当此前提不满足时，该类方法的精度会大受影响，然而，通常的遥感图像中不存在纯像元的现象又是普遍存在的。N-FINDR算法[4]是单形体几何学方法中一种较为典型的、应用广泛的算法，它也需要上述的前提假设。本文单形体算法的计算主要基于N-FINDR算法，但并不限于N-FINDR算法，其他此类端元提取算法，例如CCA[5]，VCA[6]，SGA[7]等均可同样应用于本文的方法中。

非负矩阵因式分解(Non-negative Matrix Factorization，NMF)是Daniel.D.Lee等在1999年的Nature中首先提出的方法[8]，并用于解决人脸识别和语义分析中的问题。该方法对一个非负的矩阵进行分解，得到两个非负的矩阵[8][9]，其分解模型和混合像元线性分解模型十分相似，但是，若直接将其应用于混合像元分解会存在着局部最小问题.

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

1.线性光谱混合模型

近年的研究中，线性光谱混合模型被广泛的应用于遥感图像中的混合像元分解问题，该模型假设图像中的每个像元都为各个端元像元通过线性混合得到。设V为多通道遥感图像中单一像元的多光谱矢量，W为由各类纯地物信号(端元)的多光谱或高光谱矢量所组成的反射特性矩阵，H为该像元中各类地物所占的百分比(即丰度)，N为模型的误差，则依此模型有如下关系式

V＝WH+N.(1)

若遥感图像有n个通道，其中有m类地物类型，则式中v为n×1的向量，W为n×m的矩阵，H为m×1的向量，N为n×1的向量，对于实际的多通道遥感图像，尤其是高光谱遥感图像，一般有n＞m。

同时，基于混合像元分解问题的实际物理意义，H应满足如下两个约束条件：

1)混合像元中各成分的比例H_i之和应该等于1，即：

$\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{H}_i=1.---\left(2\right)$

2)分解所得各成分的比例H_i应该在[0，1]的范围内，即：

0≤H_i≤1，(i＝1，2，...，m).(3)

2.N-FINDR算法

单形体几何学方法是目前应用广泛的一类遥感图像混合像元分解方法。这类方法将混合像元分解这一代数问题转化为高维空间中的几何问题进行处理。满足上述线性混合模型及约束条件的点集构成了n-1维空间中的凸面单形体，而端元则位于这个凸面单形体的顶点上。例如，在两维空间中，由3个端元线性混合得到的点集构成的凸面单形体即为三角形，如图1(a)所示。单形体几何学方法正是利用遥感图像数据在高维空间中的这一几何特点完成端元的提取工作。本发明中所使用的单形体方法以N-FINDR为例，其他此类算法同样可以用于本发明中。N-FINDR算法是单形体算法中较为典型的一种算法，它通过下述高维空间单形体体积计算公式获得端元光谱，

$E=\left[\begin{array}{cccc}1& 1& ...& 1\\ e_1& e_2& ...& e_m\end{array}\right],---\left(4\right)$

$\mathrm{Vol}\left(E\right)=\frac{1}{(m-1)!}\mathrm{abs}\left(\right|E\left|\right),---\left(5\right)$

其中e_i为表征第i个端元的列向量，Vol是由这m个端元所构成的单形体的体积，|·|为行列式运算符。算法将图像中的所有点均代入E进行运算，当得到的Vol最大时即找到了代表m个端元所在的点，然后再通过最小二乘方法进行解混。

可以看出，N-FINDR这类基于单形体几何学的方法必须假设在遥感图像中，每一种端元地物都必须存在至少一个纯像元，否则端元的提取便会发生偏差，从而使解混的精度也大受影响。但是，由于算法本身的计算过程，即使某些端元成分在图像中不存在纯像元，算法最终计算结束时依然会找到最接近于其纯像元的混合像元，如图1(b)所示。事实上，这种现象在实际的遥感图像中是普遍存在的。

3.NMF算法

NMF算法使用非负约束条件，通过迭代运算将一个非负矩阵(向量)V分解为两个非负矩阵W和H的乘积

V＝W×H(6)

其中W为n×r阶矩阵，H为r×m阶矩阵，r须事先给定，通常，r应小于m和n。

该算法是在W和H是非负矩阵(向量)的约束条件下使目标函数：

${|V-W\×H|}^2=\underset{i,j}{\mathrm{\Σ}}{(V_{i,j}-{(W\×H)}_{i,j})}^2---\left(7\right)$

取极小，将上式分别对W和H求导数，并取迭代步长为：

${\mathrm{\η}}_{i,j}=\frac{H_{i,j}}{{\left({\mathrm{WW}}^{T}H\right)}_{i,j}},(i=\mathrm{1,2},...,n;j=\mathrm{1,2},...,m),---\left(8\right)$

由梯度下降法就可以得到如下的NMF算法的迭代公式

$H_{i,j}\←H_{i,j}\frac{{\left(W^{T}V\right)}_{i,j}}{{\left(W^T\mathrm{WH}\right)}_{i,j}},(i=\mathrm{1,2},...,n;j=\mathrm{1,2},...,m),---\left(9\right)$

$W_{i,j}\←W_{i,j}\frac{{\left({\mathrm{VH}}^T\right)}_{i,j}}{{\left({\mathrm{WHH}}^T\right)}_{i,j}},(i=\mathrm{1,2},...,n;j=\mathrm{1,2},...,m).---\left(10\right)$

上式的更新规则已被证明是收敛的。在算法中，只用到了乘、加运算，保证了结果的非负性。另外，算法无需选择学习速率，增强了算法的适用性。关于NMF算法更详尽的分析可见参考文献[8][9]。

可以看出NMF算法的分解模型和混合像元线性分解模型是几乎一样的，同时可以自动满足非负的约束条件，因此我们考虑将其应用于遥感图像混合像元分解问题中。然而，NMF算法最初提出是应用于人脸识别等领域，该算法的迭代式只能保证收敛于目标函数的局部最小点，因此若直接将其应用于定量的遥感图像混合像元分解问题进行盲分解，则容易因陷入局部极小点而无法获得精确的分解结果。

发明内容

本发明的目的在于提出一种基于单形方法和非负矩阵分解算法的遥感图像混合像元分解方法，以便从高混合遥感数据中提取准确的端元光谱信号并获得满足全约束的丰都估计结果。

本发明提出的遥感图像混合像元分解方法，基本内容如下：

为了使原有的NMF算法能同时满足混合像元分解模型的和为一约束条件，对NMF算法做出如下修改，保持W矩阵的更新公式不变，而在每步对H的迭代运算中，令

$\stackrel{\‾}{V}=\left[\begin{array}{c}V\\ {\mathrm{\δ}1}_N^T\end{array}\right]$ $\stackrel{\‾}{W}=\left[\begin{array}{c}W\\ {\mathrm{\δ}1}_m^T\end{array}\right],---\left(11\right)$

然后使用下式更新H：

$H_{i,j}\←H_{i,j}\frac{{\left({\stackrel{\‾}{W}}^T\stackrel{\‾}{V}\right)}_{i,j}}{{\left({\stackrel{\‾}{W}}^T\stackrel{\‾}{W}H\right)}_{i,j}}.---\left(12\right)$

式中1_N和1_m分别为一个N维和一个m维元素全为1的列矢量。δ为一个拉格朗日乘子，它的大小控制丰度和为1约束的严格程度，在实际应用中，使用者可根据数据维度的大小及所要求满足约束的严格程度，自行选取合适的δ值。经过上述处理，使得原始的NMF算法在应用于遥感图象混合像元分解问题中时，能自动在每步迭代运算中约束丰度矩阵列归一，从而满足混合像元分解丰都和为一约束条件。

其次，为了获得高混合遥感数据的准确分解结果，本发明将N-FINDR算法作为NMF算法的初始值选取方法，即利用最大化单形式获得初步端元信号估计结果e₁，e₂，...，e_m，然后将其组合成的矩阵做为W初始值，代入经修改后的NMF迭代式进行更新运算，利用NMF算法可以修正端元信号矩阵的能力获得最终的分解结果。

对于一幅由线性光谱混合模型描述的n波段遥感图像，在无噪声环境下，其所有像元在n维空间中正好构成了一个m-1维的单形体(我们称这个单形体所在的m-1维子空间为数据空间)，而端元则位于这个单形体的顶点上，如图1(a)所示。单形体几何学的方法必须假设在遥感图像中，每一种端元地物都必须存在至少一个纯像元，否则端元的提取便会发生偏差，从而使解混的精度大受影响。但是，由于算法本身的计算过程，即使某些端元成分在图像中不存在纯像元，算法最终计算结束时依然会找到最接近于其纯像元的混合像元，如图1(b)所示。事实上，这种现象在实际的遥感图像中是普遍存在的。本发明将单形方法作为前端方法使用，可以从高混合遥感数据中提取出最接近光谱准确值的端元光谱矩阵。

根据上述内容，本发明方法的具体步骤可归纳如下：

1.建立线性光谱混合模型

设V为多通道遥感图像中单一像元的多光谱矢量，W为由各类纯地物信号(端元)的多光谱或高光谱矢量所组成的反射特性矩阵，H为该像元中各类地物所占的百分比(即丰度)，N为模型的误差，则依此模型有如下关系式：

V＝WH+N.(1)

若遥感图像有n个通道，其中有m类地物类型，则式中V为n×1的向量，W为n×m的矩阵，H为m×1的向量，N为n×1的向量，对于实际的多通道遥感图像，尤其是高光谱遥感图像，一般有n＞m。

同时，基于混合像元分解问题的实际物理意义，H应满足如下两个约束条件：

(1)混合像元中各成分的比例H_i之和应该等于1，即：

$\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{H}_i=1---\left(2\right)$

(2)分解所得各成分的比例H_i应该在[0，1]的范围内，即：

0≤H_i≤1，(i＝1，2，...，m)；.(3)

2.用单形体方法对端元信号进行预估计

这里以N-FINDR为例，其他此类算法同样可以用于本发明中。N-FINDR算法是单形体算法中较为典型的一种算法，它通过下述高维空间单形体体积计算公式获得端元光谱，

$E=\left[\begin{array}{cccc}1& 1& ...& 1\\ e_1& e_2& ...& e_m\end{array}\right],---\left(4\right)$

$\mathrm{Vol}\left(E\right)=\frac{1}{(m-1)!}\mathrm{abs}\left(\right|E\left|\right),---\left(5\right)$

其中e_i为表征第i个端元的列向量，Vol是由这m个端元所构成的单形体的体积，|·|为行列式运算符。将图像中的所有点均代入E进行运算，当得到的Vol最大时即找到了代表m个端元所在的点e₁，e₂...e_m，此时令

W＝[e₁，e₂...e_m]，

作为步骤(3)迭代计算的初始值，

3使用修改的非负矩阵分解算法获的最终分解结果迭代运算模型(1)中的：

V＝W×H

NMF算法在W和H是非负矩阵(向量)的约束条件下使目标函数：

${|V-W\×H|}^2=\underset{i,j}{\mathrm{\Σ}}{\left(V_{i,j}{(W\×H)}_{i,j}\right)}^2---\left(7\right)$

取极小，将上式分别对W和H求导数，并取迭代步长为：

${\mathrm{\η}}_{i,j}=\frac{H_{i,j}}{{\left({\mathrm{WW}}^{T}H\right)}_{i,j}},(i=\mathrm{1,2},...,n;j=\mathrm{1,2},...,m),---\left(8\right)$

保持W矩阵的更新公式不变，而在每步对H的迭代运算中，令：

$\stackrel{\‾}{V}=\left[\begin{array}{c}V\\ {\mathrm{\δ}1}_N^T\end{array}\right]$ $\stackrel{\‾}{W}=\left[\begin{array}{c}W\\ {\mathrm{\δ}1}_m^T\end{array}\right],---\left(11\right)$

然后使用下式更新H：

$H_{i,j}\←H_{i,j}\frac{{\left({\stackrel{\‾}{W}}^T\stackrel{\‾}{V}\right)}_{i,j}}{{\left({\stackrel{\‾}{W}}^T\stackrel{\‾}{W}H\right)}_{i,j}}.---\left(12\right)$

式中1_N和1_m分别为一个N维和一个m维元素全为1的列矢量。δ为一个拉格朗日乘子，10≤δ≤200，它的大小控制丰度和为1约束的严格程度，在本文的实验中，δ取为15，在实际应用中，使用者可根据数据维度的大小及所要求满足约束的严格程度，自行选取合适的δ值。经过上述处理，使得原始的NMF算法在应用于遥感图象混合像元分解问题中时，能自动在每步迭代运算中约束丰度矩阵列归一，从而满足混合像元分解丰都和为一约束条件。

将步骤2中获得的W作为这里W的初始值，并对W和H使用经修改后的NMF算法进行迭代运算，获得最终的分解结果。

这样做的意义在于，虽然对于高混合数据N-FINDR无法获得精确的端元信号，然而由于其给出的结果为在数据集中对各端元信号最接近的估计，因此将其作为NMF算法的初始值等价于将NMF算法的梯度下降初始点选定于遥感图像数据信息中最靠近代价函数全局最优点的位置，从该点出发，进行经修改后满足全约束的迭代运算，NMF将有效的避开局部最小，其修正后的结果可以更加准确的逼近真实的端元信号值。

本发明的优点

本发明为一种基于单形体方法和非负矩阵因式分解算法的遥感图像混合像元分解方法。其优点在于，单形体方法的预处理可以获得遥感数据中最接近端元信号准确值的预估值；采用经修改过的非负矩阵因式分解算法更新式可以同时满足混合像元分解的两个约束条件；将单形体方法获得的预估值作为非负矩阵因式分解算法的初始值，可以有效避开遥感图像混合像元分解的局部最小，从而获得更加准确的分解结果。本发明在基于多光谱和高光谱遥感图像的高精度的地物分类以及地面目标的检测和识别方面具有重要意义。

附图说明：

图1 N-FINDR算法性能示意图(a)未缺失各端元的纯像元的情况，(b)缺失端元a的纯像元的情况。

图2波段1和波段2端元提取结果对比。

图3 LANDSAT上海地区图像(a)所截取的上海出海口区域，(b)遮蔽地表所用掩膜。

图4作为参考的地物分布图(a)叶绿素，(b)泥沙。

图5叶绿素和泥沙端元提取结果(a)叶绿素，(b)泥沙。

图6 N-FINDR采用最小二乘解混结果示意图(a)叶绿素，(b)泥沙。

图7本发明解混结果示意图(a)叶绿素，(b)泥沙。

图8实验所用AVIRIS数据。

图9 MNF变换特征值排列图。

图10端元光谱矢量结果对比(a)Kaolinte，(b)Alunite。

图11 NMF得到的丰度解混结果(a)Kaolinte，(b)Alunite。

具体实施方式

1.首先，为了使原有的NMF算法能同时满足混合像元分解模型的和为一约束条件，我们对NMF算法做出如下修改，保持W矩阵的更新公式不变，而在每步对H的迭代运算中，令

$\stackrel{\‾}{V}=\left[\begin{array}{c}V\\ {\mathrm{\δ}1}_N^T\end{array}\right]$ $\stackrel{\‾}{W}=\left[\begin{array}{c}W\\ {\mathrm{\δ}1}_m^T\end{array}\right],$

然后使用下式更新H：

$H_{i,j}\←H_{i,j}\frac{{\left({\stackrel{\‾}{W}}^T\stackrel{\‾}{V}\right)}_{i,j}}{{\left({\stackrel{\‾}{W}}^T\stackrel{\‾}{W}H\right)}_{i,j}}.$

式中1_N和1_m分别为一个N维和一个m维元素全为1的列矢量。δ为一个拉格朗日乘子，它的大小控制丰度和为1约束的严格程度，在本文的实验中，δ取为15。经过上述处理，使得原始的NMF算法在应用于遥感图象混合像元分解问题中时，能自动在每步迭代运算中约束丰度矩阵的列归一，从而满足混合像元分解丰度和为一约束条件。

其次，为了获得高混合遥感数据的准确分解结果，新方法将N-FINDR算法作为NMF算法的初始值选取方法，即利用最大化单形体体积获得初步端元信号估计结果e₁，e₂，...，e_m，然后将其组合成的矩阵作为W的初始值，代入经修改后的NMF迭代式进行更新运算，利用NMF算法可以修正端元信号矩阵的能力获得最终的分解结果。这样做的意义在于，虽然对于高混合数据N-FINDR无法获得精确的端元信号，然而由于其给出的结果为在数据集中对各端元信号最接近的估计，因此将其作为NMF算法的初始值等价于将NMF算法的梯度下降初始点选定于遥感图像数据信息中最靠近代价函数全局最优点的位置，从该点出发，进行经修改后满足全约束的迭代运算，NMF将有效的避开局部最小，其修正后的结果可以更加准确的逼近真实的端元信号值。

下面，我们分别以模拟和实际遥感图像数据为例说明具体的实施方式：

1.模拟遥感图像数据

实验所用模拟数据为5波段3端元模拟图像，其端元光谱特性矩阵W取为：

$\left[\begin{array}{ccc}35& 230& 135\\ 40& 90& 210\\ 200& 30& 85\\ 120& 45& 55\\ 60& 100& 40\end{array}\right]$

混合矩阵H产生办法如下，首先，随机产生一3维向量，做归一处理，然后检查，若其中有2元素大于0.30，则认为是合格的丰度向量，将其作为H的一列，否则舍弃，按此方法产生3行10000列混合矩阵H。N为添加的信噪比SNR为30db的高斯白噪声，则V＝W×H+N为所需要的5波段3端元100×100模拟图像，且其中每一像元所含端元丰度至少有2个大于0.30。

所提议的方法(组合N-FINDR和NMF的方法，这里简称为NF-NMF)为，首先用N-FINDR方法对模拟图像进行预处理，找出其中的3个端元，然后将结果作为NMF算法的初值进行迭代获得最终分解结果.

A.端元分解结果对比

这里，我们分别引入相关系数、端元信号角度均方根误差和端元信号光谱信息偏差3种指标[7]作为评测端元提取精度标准。

1)相关系数(Correlation coefficients)  基于相关系数的本发明方法与原N-FINDR算法提取端元准确度对比结果如表1所示。

2)端元信号角度均方根误差(Spectral Angel Mapper(SAM))第i个端元信号m_i和其估计值之间的夹角θ_i定义为

${\mathrm{\&theta;}}_i\&equiv;\left(\arccos \frac{<m_i,\widehat{m_i}>}{\left|\right|m_i\left|\right|\left|\right|\widehat{m_i}\left|\right|}\right),$

则端元信号角度均方根误差ε_θ可以定义为

${\mathrm{\&epsiv;}}_{\mathrm{\&theta;}}={\left(\frac{1}{p}E\right[{\left|\right|\mathrm{\&theta;}\left|\right|}^2\left]\right)}^{1/2}.$

其中p为端元个数，θ＝[θ₁，θ₂，θ₃，…θ_p]，E[]表示求期望。该指标越接近于0，说明提取的端元与实际端元矢量越为接近。基于角度均方根误差的本发明方法与原N-FINDR算法提取端元准确度对比如表2所示。

3)端元信号光谱信息偏差(Spectral Information Divergence(SID))  端元信号光谱信息定义为：

${\mathrm{SID}}_{m_i,\widehat{m_i}}\&equiv;D<m_i|\widehat{m_i}>+D<\widehat{m_i}|m_i>,$

式中

表示第i个端元信号矢量m_i对其估计值的相对熵 $\mathrm{SID}={\left(\frac{1}{p}E\right[{\left|\right|\mathrm{\&phi;}\left|\right|}^2\left]\right)}^{1/2}$

$D<m_i|{\hat{m}}_i>\&equiv;\underset{j=1}{\overset{L}{\mathrm{\&Sigma;}}}{p}_j\log \left(\frac{p_j}{q_j}\right),$

其中 $\mathrm{\&phi;}=[{\mathrm{SID}}_{m_1,\widehat{m_1}},{\mathrm{SID}}_{m_2,\widehat{m_2}},...{\mathrm{SID}}_{m_p,\widehat{m_p}}],$ $p_j=m_{\mathrm{ij}}/\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{m}_{\mathrm{ik}},$ $q_j=\widehat{m_{\mathrm{ij}}}/\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}\widehat{m_{\mathrm{ik}}},$ m_ij表示m_i的第j个分量。该指标越接近于0，说明提取的端元与实际端元矢量越为接近。基于光谱信息偏差的本发明方法与原N-FINDR算法提取端元准确度对比如表3所示。

表1基于相关系数的本发明方法与原    表2基于角度均方根误差的本发明方法与

N-FINDR算法提取端元准确度对比      原N-FINDR算法提取端元准确度对比

表3基于光谱信息偏差的本发明方法与  表4本发明方法与原N-FINDR算法

原N-FINDR算法提取端元准确度对比    采用最小二乘丰度解混准确度对比

B.丰度解混结果对比

本发明方法与原N-FINDR算法采用最小二乘丰度解混准确度对比结果如表4所示。从上述图示及表格结果可以看出，当图像中并不存在纯像元时，本发明提出的方法在端元提取和丰度解混方面的精度都要优于原N-FINDR算法结合最小二乘算法的方法。

2.实际遥感图像数据

2.1上海地区LANDSAT图像数据

选用2000年7月14日Landsat7 ETM+拍摄的上海崇明岛地区的多光谱遥感图像，其中以第1-5、7通道的6幅图像(截取350×350)做为实验数据(如图3(a)所示)，对算法可行性进行验证。提取如图3(b)所示的掩膜，将地物遮蔽后仅对水体进行PCA分析得到特征值(^*1.0e+004)由大至小排列为：4.7215，0.6540，0.1112，0.0160，0.0073，0.0039。由PCA结果可知道该地区水体遥感图像端元数目取3较为合适，这与海洋二类水体光谱反射特性主要由水体、泥沙、叶绿素三种物质决定的实际相符合。

首先根据相关文献资料[10]对该较大区域进行端元提取和丰度解混，将结果作为进行比较的标准值。然后，进一步截取出海口附近混合较为严重，纯像元明显缺失的较小区域(图3(a)中的方框区域)对本发明提出的方法在该情况下性能优越性进行验证。所截取区域叶绿素和泥沙2种端元作为参考的准确分布如图4所示。

下面，仅从对该区域出发，分别用原N-FINDR算法和本发明提出方法完成端元提取和丰度解混工作，其结果对比如图3所示。

A.端元光谱矢量

叶绿素和泥沙端元提取结果如图5所示。对比叶绿素(波长550nm以下存在吸收峰，550nm以上反射率逐渐增大)和泥沙(分别在600nm和800nm波长附近存在双反射峰)光谱特性及ETM波段数据可知，所提取的端元很好的匹配了对应物质的光谱反射特性。同时由结果可以看出，本发明提出的方法较仅用N-FINDR算法在纯像元缺失时端元提取的准确性有了明显提高。

B.丰度解混结果

原N-FINDR算法采用最小二乘解混以及本发明提出算法得到的2种端元分布结果分别如图6和图7所示。可以看出，由于此时N-FINDR算法提取的端元实际已为混合较为严重的混合像元，端元光谱矩阵误差较大，使得最小二乘法已无法准确分解出2种端元成分的分布，尤其对于叶绿素成分。由于N-FINDR找到的对应端元像素误差较大且与图像中上部像素比较接近，造成解混时最小二乘法将这部分像素丰度倾向于完全归类于叶绿素，结果造成估计的叶绿素丰度大幅偏高(灰度分布图中表现为亮度极高)。而本发明提出算法由于改善了提取端元的精度，使得获得的丰度分布图准确性有了很大提高。

2.2 AVIRIS数据

实验使用ENVI软件自带的Cuprite地区的AVIRIS数据图8，大小为400×300，波长范围是1.99-2.48μm，共有172-221波段间的50个波段数据。对该高光谱遥感数据首先做MNF变换(如图9所示)，由MNF变换结果知该地区遥感图像端元个数应取为10-13个。本实验取端元数为11，并进行端元提取和丰度解混。其中网上提供了5种矿物：Alunite，Buddingtonite，Calcite，Kaolinite，Muscovite实地探测的空间分布图。本文以Kaolinite、Alunite两种分布较为广泛的典型矿物结果为例进行说明。实验选取Kaolinte、Alunite两种该地区典型矿物作为特征端元，首先用N-FINDR找出图像中2种矿物各自的纯像元，然后全部用混合像元予与屏蔽。对本发明所提议方法在纯像元缺失时的优越性进行验证。

用N-FINDR算法初步得到的端元矢量和经NMF迭代运算后得到的端元矢量分别与纯像元光谱矢量比较结果如下所示。Kaolinte的端元提取结果如图10(a)所示，本发明方法与原N-FINDR算法提取Kaolinite端元准确度对比结果如表5所示。Alunite的端元提取结果如图10(b)所示，本发明方法与原N-FINDR算法提取Alunite端元准确度对比结果如表6所示。另外，丰度解混结果如图11所示，对比http://speclab.cr.usgs.gov/cuprite.html网上提供的实地探测结果，也有很好的吻合。

表5本发明方法与原N-FINDR算法提取Kaolinite端元准确度对比

方法	N-FINDR	NF-NMF
			相关系数	0.9489	0.9932
角度均方根误差	0.3264	0.1459

表6本发明方法与原N-FINDR算法提取Alunite端元准确度对比

方法	N-FINDR	NF-NMF
			相关系数	0.9025	0.9935
角度均方根误差	0.2700	0.0702

从数据指标和图示中可以看出，当纯像元并不存在时，经过NMF迭代运算后得到的端元光谱矢量准确性有了明显的提高。

参考文献

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

[2]J.Boardman，“Automating spectral unmixing of AVIRIS data using convex geometryconcepts，”in Summaries 4th Annu.JPL Airborne Geoscience Workshop，vol.1，1993，JPLPub.93-26，pp.11-14.

[3]A.Ifarraguerri and C.-I Chang，“Multispectral and hyperspectral image analysis with convexcones，”，IEEE Trans.Geosci.Remote Sens.，vol.37，no.2，pp.756-770，Mar.1999.

[4]M.E.Winter，“N-FINDR：An algorithm for fast autonomous spectral end-memberdetermination in hyperspectral data，”in Proc.SPIE Conf.Imaging Spectrometry V，1999，pp.266-275.

[5]M.L.Mavrovouniotis，A.M.Harper，and A.Ifarraguerri，“Classification of pyrolysis massspectra of biological agents using convex cones，”J.Chemometrics，vol.8，pp.305-333，1994.

[6]JoséM.P.Nascimento，and JoséM.Bioucas Dias，“Vertex Component Analysis：A FastAlgorithm to Unmix Hyperspectral Data，”IEEE Transaction on Geoscience and RemoteSensing，Vol.43，No.4，April 2005

[7]Chein-I Chang，Chao-Cheng Wu，Wei-min Liu，and Yen-Chieh Ouyang，“A New GrowingMethod for Simplex-Based Endmember Extraction Algorithm”，IEEE Transaction onGeoscience and Remote Sensing，Vol.44，No.10，October 2006

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

[9]D.D.Lee and H.S.Seung，““Algorithms for non-negative matrix factorization，”inProceedings of Neural Information Processing Systems，2001，vol.13，pp.556-562.

[10]恽才兴主编，海岸带及近海卫星遥感综合应用技术(第4章第2节，悬浮泥沙遥感监测)，海洋出版社，北京，2005年1月。

Claims (1)

1.一种基于非负矩阵因式分解的遥感图像混合像元分解方法，其特征在于具体步骤如下：

(1)、建立线性光谱混合模型

设V为多通道遥感图像中单一像元的多光谱矢量，W为由各类纯地物信号的多光谱或高光谱矢量所组成的反射特性矩阵，H为该像元中各类地物所占的百分比，N为模型的误差，则依此模型有如下关系式：

V＝WH+N，(1)

若遥感图像有n个通道，其中有m类地物类型，则式中V为n×1的向量，W为n×m的矩阵，H为m×1的向量，N为n×1的向量，n＞m；

基于混合像元分解问题的实际物理意义，H应满足如下两个约束条件：

①混合像元中各成分的比例H_i之和应该等于1，即：

$\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{H}_i=1---\left(2\right)$

②分解所得各成分的比例H_i应该在[0，1]的范围内，即：

0≤H_i≤1，(i＝1，2，...，m)；(3)

(2)、用单形体方法对端元信号进行预估计

采用N-FINDR算法，通过下述高维空间单形体体积计算公式获得端元光谱：

$E=\left[\begin{array}{cccc}1& 1& ...& 1\\ e_1& e_2& ...& e_m\end{array}\right],---\left(4\right)$

$\mathrm{Vol}\left(E\right)=\frac{1}{(m-1)!}\mathrm{abs}\left(\right|E\left|\right),---\left(5\right)$

其中e_i为表征第i个端元的列向量，Vol是由这m个端元所构成的单形体的体积，|·|为行列式运算符；将图像中的所有点均代入E进行运算，当得到的Vol最大时即找到了代表m个端元所在的点e₁，e₂...e_m，此时令

W＝[e₁，e₂...e_m]，

作为步骤(3)迭代计算的初始值；

(3)、使用修改的非负矩阵分解算法获的最终分解结果

迭代运算模型(1)中的：

V＝W×H

NMF算法在W和H是非负矩阵的约束条件下使目标函数：

${|V-W\&times;H|}^2=\underset{i,j}{\mathrm{\&Sigma;}}{(V_{i,j}-{(W\&times;H)}_{i,j})}^2---\left(7\right)$

取极小，将上式分别对W和H求导数，并取迭代步长为：

${\mathrm{\&eta;}}_{i,j}=\frac{H_{i,j}}{{\left({\mathrm{WW}}^{T}H\right)}_{i,j}},i=\mathrm{1,2},...,n;j=\mathrm{1,2},...,m,---\left(8\right)$

保持W矩阵的更新公式不变，而在每步对H的迭代运算中，令：

$\stackrel{\&OverBar;}{V}=\left[\begin{array}{c}V\\ {\mathrm{\&delta;}1}_N^T\end{array}\right]$ $\stackrel{\&OverBar;}{W}=\left[\begin{array}{c}W\\ {\mathrm{\&delta;}1}_m^T\end{array}\right],---\left(11\right)$

然后使用下式更新H：

$H_{i,j}\&LeftArrow;H_{i,j}\frac{{\left({\stackrel{\&OverBar;}{W}}^T\stackrel{\&OverBar;}{V}\right)}_{i,j}}{{\left({\stackrel{\&OverBar;}{W}}^T\stackrel{\&OverBar;}{W}H\right)}_{i,j}}.---\left(12\right)$

式中1_N和1_m分别为一个N维和一个m维元素全为1的列矢量，δ为一个拉格朗日乘子，10≤δ≤200。

Images