CN101866424A - 基于独立分量分析的高光谱遥感图像混合像元分解方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，具体为一种基于独立分量分析的高光谱遥感图像混合像元分解方法。本发明根据高光谱图像的物理特性，在独立分量分析的目标函数中引入丰度和为一约束与丰度非负约束，并且提出一种自适应的丰度建模方法来描述数据的概率分布，从而对各种不同的遥感数据都表现出良好的适用性。本方法可以有效地解决高混合度和各种干扰下的遥感数据混合像元分解问题。在基于多光谱和高光谱遥感图像的高精度的地物分类以及地面目标的检测和识别方面具有重要的应用价值。

Description

基于独立分量分析的高光谱遥感图像混合像元分解方法

技术领域

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

背景技术

遥感是本世纪六十年代发展起来的新兴综合技术，与空间、电子光学、计算机、地理学等科学技术紧密相关，是研究地球资源环境的最有力的技术手段之一。近年来，随着成像技术的进步，多波段遥感图像在越来越多的领域得到了广泛应用。由于成像系统空间分辨率的限制和地表的复杂多样，所获得的遥感图像中的一个像元往往包含着多种地物类型，这就形成了混合像元。如何从混合像元中提取典型地物的光谱，并且求出这些地物所占比例，在实际应用中有着重要意义[1]。

大部分混合像元分解方法都是基于线性混合模型的[2]。基于线性混合模型的像元解混可分为两个步骤，端元提取和丰度解混。端元提取是指从遥感数据中求取端元光谱的过程。丰度解混则是指确定各个端元在混合像元中所占的比例。

独立分量分析(Independent Component Analysis，ICA)[3]把混合像元分解问题视为一种盲分离问题，它们无需光谱信息就可以进行混合像元分解。该方法假设源信号之间相互独立，通过让分离结果彼此独立来恢复出源信号。以观测到的高光谱遥感图像作为混合信号，端元光谱或丰度作为源信号，ICA就可以应用于混合像元的解混，具有运算量小、不易陷入局部极小的优点。但是，直接使用ICA进行光谱解混所面临的一个重要问题是：算法的独立性假设直接与线性混合模型中的和为一约束相矛盾，导致解混结果和真实值之间出现偏差[4]。另一个问题是，ICA假设数据的统计分布规律是一定的，但是真实图像的统计特性各不相同，一般都与事先假设的情况有差别，这就影响了算法的性能。后一问题较难克服，因为ICA是一种盲分离方法，除观测值以外没有其他先验知识，只能假设数据的统计特性为某种确定的分布，然后让结果去逼近这种分布。

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

1线性混合模型

线性混合模型是分析研究高光谱遥感图像的一种简便有效的数据模型。在高光谱数据中，对于每个波段图像中的每一像元，其灰度值可表示为该像元中各端元的光谱特性与各端元在像元中所占比例的线性组合。设高光谱图像中共有L个波段，P个端元，则某一个像素的观察值

可被表示为

x＝As+e.                               (8)

这里，x＝[x₁，x₂，...，x_L]^T；

是端元光谱矩阵，其中

代表第i个端元的光谱向量；

是丰度向量，代表该像素中各端元所占的比例；

表示由于噪声或模型不准等因素导致的误差。从物理实际出发，任意一个像元中各地物所占比例应该是非负的，并且其比例之和等于1。因此，丰度向量s应满足丰度非负约束(Abundance Nonnegative Constraint，ANC)以及和为一约束(Abundance Sum-to-oneConstraint，ASC)

$\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{s}_i=\mathrm{1,0}\≤s_i\≤1,(i=\mathrm{1,2},...,P).---\left(9\right)$

设每个波段的高光谱图像都有N个像素点，则所有这些像素点的观察值可以组成一个矩阵

从而(8)式可表示成矩阵形式

X＝AS+E.                                (10)

这里，X是多通道遥感图像的观测矩阵，其每一列代表单一像元的观测光谱矢量，

是丰度矩阵，其每一行代表一个端元的丰度。

2独立分量分析

在线性模型下，混合像元分解可被视为一种线性解混问题，因而可以通过线性解混算法ICA进行求解。ICA的目标是寻找一种线性变换，使得分离出的源信号彼此相互独立。具体地，我们可以将(8)中的x，A，s分别视为观测信号、混合矩阵、源信号，求出一个解混矩阵W使得

y＝Wx＝Us，                        (11)

这里，y＝[y₁，y₂，...y_P]^T，是源信号s的估计值，ICA通常把y_i视为随机变量；U是每行每列只有一个非0元素的矩阵。

W的求解通过使y的各分量间独立性最大化来进行，衡量独立性的判据有很多，本发明的算法建立在互信息判据的基础上[5][8]。y的各分量之间的互信息的定义为

$I\left(y\right)\stackrel{\mathrm{def}}{=}{D}_{\mathrm{KL}}\left(p\left(y\right)\right|\left|{\mathrm{\Π}}_{i=1}^{P}p\left(y_i\right)\right)=\∫p\left(y\right)\log (p\left(y\right)/{\mathrm{\Π}}_{i=1}^{P}p\left(y_i\right))\mathrm{dy}.---\left(12\right)$

这里，p(y_i)是y_i的概率密度函数，p(y)是y＝[y₁，y₂，...y_p]^T的联合概率密度函数。可见I(y)就是p(y)与

间的K-L散度，也就是两者在概率意义上的距离。互信息总是非负，当且仅当y的各分量互相独立时I(y)＝0。以互信息为目标函数，ICA的目标就可被写作

$\mathrm{Minimize\; I}\left(y\right)=\∫p\left(y\right)\log \left(p\left(y\right)/{\mathrm{\Π}}_{i=1}^{P}p\left(y_i\right)\right)\mathrm{dy},$ 其中y＝Wx.

I(y)对W求负梯度，有

$-\frac{\∂I\left(y\right)}{\∂W}=W^{-T}-f\left(y\right)y^T,---\left(13\right)$

其中，

f(y)＝[-p′(y₁)/p(y₁)，-p′(y₂)/p(y₂)，...，-p′(y_p)/p(y_p)]^T.        (14)

这里p′(y_i)是p(y_i)的导函数，W^-T＝[W^T]^-1表示W的转置的逆矩阵。

根据梯度下降算法，解混矩阵W的迭代公式为

W←W+ηΔW，                    (15)

$\mathrm{\ΔW}=-\frac{\∂I\left(y\right)}{\∂W}.---\left(16\right)$

其中，η是步长因子。当使用自然梯度[6]时，(16)式更新为

ΔW＝(I-f(y)y^T)W.                (17)

这里，

是单位阵。

由(17)式和(14)式可知，求解W需要已知p(y_i)，但实际工作中p(y_i)是很难准确获取的，如果估计取值和真实值相差较大，必然对解混效果造成不良影响。更近一步，根据(11)式可知，ICA并不能保证结果的幅值和源信号的幅值相同，这将严重影响分解结果。本发明针对以上问题，提出了完整的解决方案。

发明内容

本发明的目的在于提出一种鲁棒性好、适应性广、计算效率高的高光谱遥感图像混合像元分解方法。

本发明提出的高光谱遥感图像混合像元分解方法，是一种有约束的独立分量分析方法(Constrained Independent Component Analysis，CICA)，其内容主要包括两个部分：自适应丰度建模(Adaptive Abundance Modeling，AAM)和约束条件。AAM的作用是克服传统ICA只能描述固定分布的缺点，使结果能够更接近于理想值；约束条件用于调整传统ICA算法的目标函数，使算法从物理意义上适用于高光谱混合像元解混问题。

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

1自适应丰度建模

采用Pearson混合模型[7]近似p(y_i)为

$p\left(y_i\right)=(1-{\mathrm{\α}}_i)N({\mathrm{\μ}}_{1\left(i\right)},{\mathrm{\σ}}_i^2)+{\mathrm{\α}}_{i}N({\mathrm{\μ}}_{2\left(i\right)},{\mathrm{\σ}}_i^2),---\left(18\right)$

其中，α_i是权重，0＜α_i＜1，N(μ，σ²)代表均值为μ方差为σ²的正态分布函数，α_i，σ_i，μ_1(i)，和μ_2(i)都是待定参数。这些参数的不同取值决定了函数p(y_i)的形状。

一种经典的ICA算法[8]通过事先指定所有参数的值来确定p(y_i)，也就是说，该方法假设数据的分布规律为某已知的概率分布。这种预定义的p(y_i)和真实地物分布的情况一般是不同的，会导致分解结果有较大偏差。本节提出的AAM方法可以克服这一问题，下面介绍AAM参数学习规则。

(1)参数μ_2(i)。通常，从物理实际出发，遥感图像中的各端元都有各自的分布区域，一般不会充满整幅图像，这种性质称为丰度分布的稀疏性。具体表现为：丰度数据中会出现较多的0值或接近于0的值。该性质可视为y_i在0附近的概率取峰值，因此，为简化问题，我们直接令μ_2(i)＝0，(i＝1，...，P)。图1显示了μ_2(i)＝0时，p(y_i)在区间y_i∈[0，1]的函数曲线(其余参数为α_i＝0.8，μ_1(i)＝1，σ_i＝1/3)。

(2)参数α_i。我们首先给出性质1，详细的证明可见附录A。

性质1.对于概率分布为(18)，且μ_2(i)＝0的信号y_i，其偏度γ₁(y_i)和峭度γ₂(y_i)具有以下性质

a)当0＜α_i＜0.5时，γ₁(y_i)＜0；

b)当0.5＜α_i＜1时，γ₁(y_i)＞0；

c)当

或时，γ₂(y_i)＞0；

d)当

时，γ₂(y_i)＜0。

性质1说明，信号的偏度、峭度与α_i的取值相对应。可以通过α_i的不同取值来建模不同偏度和峭度的信号，反之，也可以根据信号的偏度、峭度确定α_i的值。

实际工作中发现，对于本算法而言，峭度对结果的影响远比偏度重要。为减少计算量，可以只根据峭度γ₂(y_i)测定α_i的值，偏度则设定为γ₁(y_i)＞0恒成立。这时γ₂(y_i)＞0等价于

γ₂(y_i)＜0等价于

上述区间的中值分别是0.8943和0.6643，从而可以给出α_i的调整规则如下：每次迭代时计算当前所得y_i的峭度，如果γ₂(y_i)＞0，令α_i＝0.8943；若γ₂(y_i)＜0，令α_i＝0.6443。

(3)参数μ_1(i)。我们先给出性质2，详细的证明可见附录A。

性质2.对于概率分布为(18)，且满足μ_2(i)＝0的信号y_i，其均值为

mean(y_i)＝(1-α_i)μ_1(i).                (19)

根据性质2，可以给出参数μ_1(i)的调整规则如下：每次迭代中，在确定α_i后计算y_i的均值mean(y_i)，然后令μ_1(i)＝mean(y_i)/(1-α_i)。

(4)参数σ_i。丰度的物理意义决定了必然有y_i∈[0，1]，因此p(0≤y_i≤1)应该尽可能地大。附录B从这点出发做出了σ_i关于取值的分析，取σ_i＝1/3。

至此，AAM中的所有参数都已给出，总结AAM参数学习规则如下：

Step1.置μ_2(i)＝0，σ_i＝1/3。

Step2.每次迭代，按照如下步骤获得α_i和μ_1(i)。

(i).计算当前所得y_i的峭度γ₂(y_i)。

(ii).如果γ₂(y_i)＞0，令α_i＝0.8943；若γ₂(y_i)＜0，令α_i＝0.6443。

(iii).计算y_i的均值mean(y_i)，令μ_1(i)＝mean(y_i)/(1-α_i)。

根据上述调整规则，将μ_2(i)＝0，σ_i＝1/3代入(18)式，可以得出

$-\frac{p^{\′}\left(y_i\right)}{p\left(y_i\right)}=9y_i+\frac{{9\mathrm{\μ}}_{1\left(i\right)}}{[\frac{{\mathrm{\α}}_i}{{\mathrm{\α}}_i-1}\exp (-9{\mathrm{\μ}}_{1\left(i\right)}{y}_i+4.5{\mathrm{\μ}}_{1\left(i\right)}^2)-1]}.---\left(20\right)$

将(20)式代入(17)式，并以符号ΔW_AAM代替ΔW，得出W的更新公式为

ΔW_AAM＝(I+dy^T-9yy^T)W.            (21)

其中，d＝9μ₁·/[1-α·/(α-1)·*exp(-9μ₁·*y+4.5μ₁·*μ₁)]，μ₁＝mean(y)./(1-α)，α＝sign(γ₂)/8+0.7693，mean(y)＝[mean(y₁)，mean(y₂)，...，mean(y_p)]^T，sign是符号函数，γ₂＝[γ₂(y₁)，...，γ₂(y_k)，...，γ₂(y_P)]^T，

上式中的″.*″和″./″表示矩阵元素对应地相乘和相除。

这样就针对单个像素给出了ICA的迭代公式。但是，在实际处理高光谱遥感问题时，都是采用矩阵形式。在矩阵形式下，线性模型表示为(10)式，ICA解混公式表示为

Y＝WX.                           (22)

定义矩阵

的第j列为从而Y的互信息表示为

$I\left(Y\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}I\left(y_j\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{D}_{\mathrm{KL}}\left(p\left(y_j\right)\right|\left|{\mathrm{\Π}}_{i=1}^{P}p\left(y_{\mathrm{ij}}\right)\right).---\left(23\right)$

类似地，通过最小化I(Y)，可给出矩阵形式下W的迭代公式

ΔW_AAM＝(IN+DY^T-9YY^T)W/N.                      (24)

2约束条件

(1)非负约束。线性混合模型的ANC条件要求：丰度向量的各分量都应在[0，1]范围内取值。我们可以通过在算法中引入非负约束项来使结果满足ANC。约束项的设计原则是，当结果满足ANC时约束项等于0，否则大于0。

基于上述原则，我们提出一种非负约束的目标函数如下

$J\_\mathrm{ANC}\left(Y\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}[f\left(y_{\mathrm{ij}}\right)+|f\left(y_{\mathrm{ij}}\right)\left|\right]/2,---\left(25\right)$

其中，f(y_ij)是任意一种可以满足以下条件的函数

$f\left(y_{\mathrm{ij}}\right)\left\{\begin{array}{cc}=0,& y_{\mathrm{ij}}\∈\left[\mathrm{0,1}\right]\\ >0,& y_{\mathrm{ij}}\∉\left[\mathrm{0,1}\right]\end{array}\right..---\left(26\right)$

在本发明中取

式(27)可以保证，当矩阵Y中的负分量出现地越多、负分量的绝对值越大时，J_ANC(Y)的值就越大。这一点有助于加快算法的收敛速度。以符号ΔW_ANC表示J_ANC(Y)对W的负自然梯度，并取正整数b为1，可得

$\mathrm{\ΔW}\_\mathrm{ANC}=-\frac{\∂J\_\mathrm{ANC}\left(Y\right)}{\∂W}{W}^{T}W=-G{X}^T{W}^{T}W.---\left(28\right)$

其中，矩阵G的第i行j列位置处的元素为

$g_{\mathrm{ij}}=\left\{\begin{array}{cc}{y}_{\mathrm{ij}}-0.5,& y_{\mathrm{ij}}\∉\left[\mathrm{0,1}\right]\\ 0,& y_{\mathrm{ij}}\∈\left[\mathrm{0,1}\right]\end{array}\right..---\left(29\right)$

(2)和为一约束。线性混合模型的ASC条件要求：对于任意像素点，所有端元地物在该点的丰度之和必须为1。本发明定义ASC目标函数为

$J\_\mathrm{ASC}\left(Y\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{y}_{\mathrm{ij}}-1)}^{2c}/2c.---\left(30\right)$

这里

当Y满足ASC时，J_ASC(Y)＝0，否则J_ASC(Y)＞0。取c＝1，并令J_ASC(Y)对W求负自然梯度，可得

$\mathrm{\ΔW}\_\mathrm{ASC}=-\frac{\∂J\_\mathrm{ASC}\left(Y\right)}{\∂W}{W}^{T}W=-H{X}^T{W}^{T}W,---\left(31\right)$

其中，矩阵H的第i行j列位置处的元素为

$h_{\mathrm{ij}}=(\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{y}_{\mathrm{ij}}-1).---\left(32\right)$

3.3CICA算法的流程

综合考虑互信息、非负约束以及和为一约束，CICA问题可以被总结为

最小化J(Y)＝I(Y)+η₁J_ANC(Y)+η₂J_ASC(Y)，其中Y＝WX.

参数η₁，η₂用于控制约束条件的权重，这三个参数根据具体计算确定。J(Y)就是CICA算法的总目标函数，根据梯度下降法，可以得到解混矩阵W的迭代公式为

W←W+ηΔW，ΔW＝ΔW_AAM+η₁ΔW_ANC+η₂ΔW_ASC，        (33)

这里，ΔW_AAM，ΔW_ANC，以及ΔW_ASC均已在前面作过说明。

只要给出W的初始值，就可以开始(33)式的算法迭代。选取合理的初始值可简化问题的复杂性。不同于传统ICA的随机初始化，为使算法快速收敛，我们采用如下的初始化方法：使用ATGP[9]从观测图像中选取P个像素点，认为它们是纯像元并以它们的值初始化端元矩阵A，然后以A的伪逆阵作为W的初始值。本算法的具体步骤可以总结如下：

算法基于有约束ICA的遥感图像混合像元分解(CICA)

已知观测矩阵端元个数P。

Step 1.初始化W。

Step 2.开始迭代。

2a)根据式(33)更新W。

2b)根据式(22)计算Y。

2c)如果算法没有收敛，返回2a)继续迭代。

Step 3.输出结果。

如果要获取端元矩阵A的值，可以在求出丰度矩阵Y后，通过解非负最小二乘问题来实现。也可以使用NMF[10]来实现，即以Y作为NMF的初始值，进行迭代求解。

本发明的优点

本发明为一种新的基于ICA算法的混合像元分解方法。本发明方法是根据高光谱图像的特点，提出了一种自适应丰度建模(AAM)方法，并向ICA算法中引入约束条件，从而服了ICA算法独立性假设问题和概率分布统计不变问题，是一种可以有效的解决高混合遥感数据的混合像元分解方法。其中的AAM的学习规则根据图像的物理特性设计，可以按照遥感数据的特点进行自动调整，使算法对于多种数据都有良好的适用性。另外，本发明引入了ANC与ASC约束条件，它们具有明确的物理意义且计算简便。新方法在基于多光谱和高光谱遥感图像的高精度的地物分类以及地面目标的检测和识别方面具有特别重要的应用价值。

仿真实验结果表明，所提出的算法克服了传统ICA方法的问题，能得出更优的解，对纯像元缺失和受到干扰的数据，都表现出了较好的鲁棒性。同时，本方法的一个显著特点是，在端元数估计错误的情况下仍能得到正确结果，具有重要的实际意义。对于实际高光谱遥感数据实验，该方法也得到了理想的结果，进一步证实了方法的有效性和对于各种不同数据的适用性。

附图说明

图1为p(y_i)函数曲线。

图3为算法的抗噪声性能。其中，(a)RMSE，(b)SAD。

图2为丰度信号图。

图4为像素数对算法的影响。(a)RMSE，(b)SAD。

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

图6为时的理想结果。

图7为端元数估计误差对算法的影响。

图8为自适应概率模型对算法的影响。(a)RMSE，(b)SAD。

图9为Indiana数据的伪彩色图。

图10为Indiana数据的分解结果。(a)干草堆，(b)大豆，(c)草地，(d)玉米，(e)树林，(f)人工建筑。

图11为Cuprite数据的伪彩色图。

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

图13Urban数据的伪彩色图。

图14Urban数据的丰度分解结果。(a)水泥路面，(b)沥青路面，(c)树木，(d)草地，(e)屋顶1，(f)屋顶2。

具体实施方式

下面，分别用仿真数据和实际遥感图像数据为例说明本发明的具体的实施方式：

1仿真数据

我们将本发明提出的CICA方法与以下3种典型算法进行对比分析：HOS-ICPA[11]，VCA[12]，以及MVCNMF[13]。其中VCA算法只能得到光谱矩阵，因此在解出光谱后用FCLS[14]求出丰度矩阵，并将这种方法记为VCA-FCLS。用仿真数据测试以上所有算法的性能，以SAD和RMSE两个指标衡量所有算法解出的结果和真实参考值之间的差异。

光谱角距离(Spectral Angle Distance，SAD)用于衡量算法解出的光谱与已知的参考光谱间的差异程度，第i个端元的真实光谱矢量a_i＝[a_i1，...，a_ik，...，a_iL]^T与其对应的解混结果

之间的SAD定义为

${\mathrm{SAD}}_i=\arccos \frac{a_i^T{\tilde{a}}_i}{\left|\right|a_i\left|\right|\·\left|\right|{\tilde{a}}_i\left|\right|},(i=1,...P).---\left(34\right)$

均方根误差(Root Mean Square Error，RMSE)是衡量丰度解混的指标。设第i个端元的真实丰度矢量为s_i＝[s_i1，...，s_ij，...，s_iN]，它与其对应的解混结果y_i＝[y_i1，...，y_ij，...，y_iN]之间的RMSE定义为

${\mathrm{RMSE}}_i=\sqrt{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(y_{\mathrm{ij}}-s_{\mathrm{ij}})}^2},(i=1,...P).---\left(35\right)$

在进行实验测试时，以所有端元的平均SAD和平均RMSE作为评价标准。本实验使用的仿真数据含有五个端元，其丰度数据遵循以下原则：分布是连续的、渐变的，丰度值逐渐从最大值过渡到最小值，满足和ANC和ASC。对于光谱数据，我们从美国地质勘测局(USGS)公布的矿物光谱库(Available at http://speclab.cr.usgs.gov/spectral.lib04/lib04-AVIRIS.html)中选取，共有224个波段。这样产生的仿真数据，其特性能够尽可能地接近真实遥感数据，有利于评价算法性能。

具体的丰度产生方法是：为5种端元分别产生5幅大小为r×r像素的丰度图，丰度图上各点的灰度值代表该点的丰度，其取值范围是[0，1]。每幅图上各选一个点作为“中心点”，该点的值设为数据的纯度。中心点以外的点，其值逐渐递减。最后把所有端元的丰度信号之和归一化以满足ASC。这里，纯度的取值范围是[0.5，1]。图2显示了令r＝100，纯度为1时，得到的5个地物的丰度图。其中(a)-(e)各代表一种端元的丰度分布，都是100×100的正方形区域。这里丰度图既有超高斯分布又有亚高斯分布：(a)是亚高斯，(b)-(e)是超高斯，这样可以近似多种类型的分布情况，使仿真数据在评价算法时更具有代表性。

通过四个实验评价本发明所提出方法的有效性，并与其他算法进行比较。四个实验分别测试噪声强度、像元个数、纯像元缺失、端元数估计误差对算法性能的影响。实验中CICA算法的参数统一设置为：η＝0.1/P，η₁＝5/N，η₃＝20/N。

实验1抗噪声干扰的性能。我们按照(8)式把高斯白噪声加到仿真数据中，然后测试不同算法的解混结果。一般说来，在真实遥感图像中，通常的情形是：有些地物存在纯像元，有些地物没有纯像元，为了让仿真数据近似这种情况，在实验中指定3个端元的纯度为0.7(无纯像元)，其余端元的纯度为1(有纯像元)。图3(a)和(b)分别给出了在不同的信噪比下，各算法的RMSE和SAD的变化情况。可见，随着信噪比的降低，所有算法的整体效果在都在逐渐变差。CICA算法的效果始终是最好的，其次是VCA-FCLS和MVCNMF，HOS-ICPA则基本没有解出正确结果。

实验2算法性能与像素个数的关系。这个实验中，改变仿真图像的大小，以此来比较像素个数对各算法的影响。根据图像大小的不同，从50×50到100×100像素，总共分为5档。图4(a)和(b)分别给出了在不同的像元个数下，各算法的RMSE和SAD的变化情况。可见，随着像元个数的增加，所有方法的性能都有小幅改善。VCA是基于单形体的方法，像素数的增加有利于对单形体顶点的寻找；CICA和HOS-ICPA是基于统计的方法，样本数的增加意味着能够获取更多的统计信息；MVCNMF是基于NMF是方法，随着像素个数的增加可以获得更多信息，因而效果就随之提高。但是NMF算法具有局部极小，所以图中也出现了随着像素个数提高而效果下降的情况，这在一定程度上反应了该方法的不稳定性。

为衡量各算法的计算复杂度，以运算耗时(以“秒”为计时单位)作为衡量指标，在表1中给出了各算法的耗时情况。实验中所使用的计算机软硬件环境如下：CPU为Xeon(TM)3.60GHZ，内存3.00GB，所有程序在MATLAB中开发实现。可见，VCA，CICA和HOS-ICPA的计算耗时相差无几，而MVCNMF的耗时相当多，大约是其他算法的1000倍。这是因为MVCNMF基于NMF，计算量较大，而且该方法还涉及单形体体积的计算，进一步增加了算法复杂度。VCA-FCLS方法中，VCA算法迭代次数和端元数相同，FCLS算法解最小二乘问题，二者的运算量都很小。HOS-ICPA和CICA都是ICA类方法，经过降维等预处理后简化了问题复杂度，耗时都较少。本实验说明，算法能够以低运算耗时取得理想的结果。

表1同像素个数下的运算时间(时间单位：秒)

实验3纯像元缺失情况下的适应性。一般而言，实际遥感图像中并非所有地物都存在纯像元。因此，纯像元缺失时能否得到良好的分解结果是衡量算法性能优劣的一个重要指标。本实验中，让纯度从1改变到0.5，比较各算法对纯像元缺失情况的适应性。图5(a)和图5(b)分别给出了在不同的纯像元缺失程度下，各算法的RMSE和SAD的变化情况。可见，随着像元纯度的降低，CICA和VCA-FCLS的性能都在逐渐变差，但是CICA在效果上始终优于VCA-FCLS。而MVCNMF和HOS-ICPA则基本上没有解出正确的结果，RMSE和SAD两项指标都始终处于较高的水平。

实验4抗端元数估计误差的鲁棒性。现有的混合像元解混算法必须精确已知图像中存在的端元的个数P，这就需要根据先验知识或者使用专门的算法来进行端元数估计。但是，由于数据集本身的复杂性，端元个数的估计一般是不准确的。事实上，对于同样的实际遥感数据，不同的端元数估计算法求出的端元个数一般是不同的，同一种估计算法在不同的参数下给出的端元个数也是不同的[15]。因此，能否在不准确的端元数估计下依然给出正确的结果具有重要的实际意义，本实验对所提议算法抗端元数估计误差的性能进行分析。

仿真数据的端元数P＝5，假设端元数被错误地估计为一个大于5的数P′。例如P′＝8，设纯度为1，这时所提议算法的解混结果如图6所示。进一步地，令P′从5～10变化，以此模拟不同程度的端元数估计误差。我们比较各种算法的解混结果和理想结果间的RMSE和SAD，如图7所示。这里没有放上VCA-FCLS的结果，因为P＞5时，它的RMSE太大和其他方法不具有可比性。可见，随着端元数估计误差的增大，CICA几乎没有受到影响，总是能得出正确的结果，而其它几种方法的都没能得出正确解。本实验说明，在端元数估计偏大的情况下，CICA不但能求出实际存在的端元丰度，同时还能以全0的形式给出被错误地多估计的端元的丰度，从而指出正确的端元个数。这一优点是其它算法所不具备的。

实验5AAM对结果的影响。自适应的概率模型是CICA算法的一个重要组成部分，因此，有必要考察该模型的引入是否提高了算法的性能。我们将CICA的模型(18)的参数固定，其余操作不变，从而得到一种新算法，记为CICA2。比较CICA2和CICA的解混结果，以此验证自适应模型的有效性。类似于文献[8]，本实验CICA2的模型(18)中的四个参数被固定为α_i＝0.5，μ_1(i)＝1，μ_2(i)＝-1，σ_i＝1。图8(a)和图8(b)分别给出了在不同的纯像元缺失程度下，两种算法相应的RMSE和SAD的变化情况。可见，在纯度大于0.9时CICA的分解效果和CICA2的效果差别不大；当纯度低于0.9时，CICA明显优于CICA2。这是因为，数据纯度的变化也意味着概率分布的变化，而自适应的概率模型能够比固定参数更好地描述不同形式的分布，所以可以在一定程度上保证混合像元分解的精度。

2实际遥感数据

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

A.Indiana数据实验

本实验采用了美国印第安纳州西拉斐特(West Lafayette)地区的高光谱遥感数据(Available at http://cobweb.ecn.purdue.edu/～biehl/MultiSpec/)。该数据由机载可见光及红外成像光谱仪(Airborne Visible/Infrared Imaging Spectrometer，AVIRIS)拍摄于1992年6月，成像区域是Purdue大学的一块遥感试验田，图像大小为145×145，波长范围从0.4～2.5μm，光谱分辨率为10nm，空间分辨率为17m，共有224个波段。取第70，86，136波段分别作为R，G，B分量合成伪彩色图如图9所示。由Purdue大学的实地勘测报告[16]可知，该地区的地表主要由各种农作物和人工建筑所覆盖，典型地物包括玉米、大豆、干草、树林、草地、公路、石塔以及一些房屋。

数据的第1～4，78～82，103～115，148～166以及211～220波段由于为水吸收波段或信噪比太低而被事先移除，本发明使用剩下的169个波段进行实验。为定量衡量算法的性能，我们根据地物真实报告所提供的分布情况对端元进行手动提取，在图像中为各地物分别选择多个较纯像元，通过求出每类像元观测值的平均值，以此作为各端元的参考光谱。本数据中总共提取了六个端元光谱，分别对应人工建筑、玉米、大豆、草地、树林和干草堆。

采用本发明的方法对该数据的分解结果如图10所示，其中各参数的选取均与仿真数据相同。解出的六个端元分别是人工建筑、玉米、大豆、草地、树林和干草堆，通过比较可以发现，丰度解混的结果与真实分布[16]非常吻合。为比较CICA，VCA，HOS-ICPA和MVCNMF四种方法的解混结果，我们将所有算法解出的光谱与手动选取的参考光谱进行比较，求出其光谱角，结果如表2所示。可见，CICA表现出了最好的性能，其余三种方法则误差过大。

B.Cuprite数据实验

本实验采用的数据由机载可见光及红外成像光谱仪(AVIRIS)拍摄于美国内华达州Cuprite地区(Available at http://aviris.jpl.nasa.gov/html/aviris.freedata.html)。数据成像于1997年6月19日，波长范围0.37～2.48μm，光谱分辨率10nm，共有224个波段。实验所用图像大小为250×191，图11显示了该数据的伪彩色图。该数据已被广泛应用于高光谱图像解混算法的评价中，文献[17]给出了其地物真实分布的报告。Cuprite区域的地表主要为裸露的矿物，除了个别地物的分布比较突出外，其他物质的混合度都比较高。

在算法进行解混前，去数据的除低信噪比波段、水吸收波段(包括波段1～2，104～113，148～167以及221～224)，使用剩下的188个波段进行实验。使用CICA所得到的丰度解混结果如图12所示。与实地勘测地物分布相比较，可以确定所求出的这些矿物，具体包括：明矾石Alunite，多水高岭石Halloysite，白云母Muscovite，玉髓Chalcedony，蒙脱石Montmorillonite，高岭石Kaolinite，皂石Nortronite，黄甲铁石Jarosite以及沙漠地表Desert Vanish。为定量衡量算法性能，以美国地质勘测局(USGS)公布的矿物光谱库作为参考光谱，求出几种不同算法的解混结果与参考值之间的光谱角，结果如表3所示。通过比较可以看出，本发明的方法总体上给出了最好的结果。

表2Indiana数据的光谱分解结果与参考值之间的光谱角(SAD)

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

C.Urban数据实验

本实验采用的数据是超光谱图像收集实验仪器(Hyperspectral Digital ImageryCollection Experiment，HYDICE)拍摄的Urban数据，成像于1997年6月19曰，该数据是在4430m的飞行高度上拍摄的，空间分辨率2m，光谱分辨率10nm，共有210个波段(Available at http://www.tec.army.mil/Hypercube)。图像大小为307×307。图13给出了以第12,99和171波段分别作为R、G、B分量合成的伪彩色图。Urban数据的成像区域位于美国得克萨斯州胡德堡(Fort Hood，TX)附近的科珀勒斯科夫镇(Copperas Cove)，典型地物有路面，房屋，草坪和树木等。其中，根据质地的不同，路面可以分成沥青路面和水泥路面，二者的光谱反射特性是不同的；房顶根据反射率的高低则可以分成较亮屋顶和较暗屋顶两种。因此可以认为本数据中存在6个端元，分别是：沥青路面，水泥路面，屋顶1，屋顶2，草地，树木。

在解混前，第1～4，76，87，101～111，136～153以及198～210波段由于信噪比太低或为水吸收波段而被移除，剩下的162个波段用于进行算法的验证实验。使用CICA算法得出的丰度解混结果如图14所示，可见，分解结果符合地物真实分布情况。为进一步衡量算法性能，和Indiana数据实验类似，我们参考实际地物分布情况，通过对图像中的端元进行手动提取得出相应的参考光谱。表给出了所有算法的解混结果和参考光谱之间的光谱角，显然本发明的方法表现出了最好的综合性能。

表4Urban数据的光谱分解结果与参考值之间的光谱角(SAD)

参考文献

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

[2]T.M.Lillesand and R.W.Kiefer，Remote Sensing and Image Interpretation，4thed.New York：John Wiley & Sons，Inc.，2000.

[3]A.Hyvarinen，J.Karhunen，and E.Oja，Independent Component Analysis，New York：Wiley，2001.

[4]J.Nascimento and J.Bioucas-Dias，“Does Independent Component Analysis Playa Role in Unmixing Hyperspectral Data？”IEEE Transactions on Geoscience andRemote Sensing，vol.43，no.4，pp.175-187，April，2005.

[5]H.H.Yang，and Amari，S.Adaptive Online Learning Algorithms for BlindSeparation：Maximum Entropy and Minimum Mutual Information，Neural Computation，vol.9，no.7，pp.1457-1482，1997.

[6]Shun-Ichi Amari，“Natural gradient works efficiently in learning，”NeuralComputation，vol.10 no.2，pp.251-276，Feb.15，1998.

[7]K.Pearson，Contributions to the theory of mathematical evolution.Phil.Trans.R.Soc.Lond.，185，pp.71-110，1894.

[8]M.Girolami，“An alternative perspective on adaptive independent componentanalysis algorithms，”Neural Computation，vol.10，pp.2103-2114，1998.

[9]H.Ren and C.-I Chang，“Automatic spectral target recognition in hyperspectralimagery，”IEEE Trans.Aerosp.Electron.Syst.，vol.39，no.4，pp.1232-1249，Oct.2003.

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

[11]J.Wang and C-I Chang，“Applications of independent component analysisin endmember extraction and abundance quantification for hyperspectralimagery，”IEEE Transactions on Geoscience and Remote Sensing，vol.44，no.9，pp.2601-2616，Sep.2006.

[12]J.Nascimento and J.Bioucas-Dias，“Vertex Component Analysis：A FastAlgorithm to Unmix Hyperspectral Data，”IEEE Transactions on Geoscience andRemote Sensing，vol.43，no.4，pp.898-910，April，2002.

[13]L.Miao and H.Qi，“Endmember extraction from highly mixed data usingminimum volume constrained nonnegative matrix factorization，”IEEETransactions on Geoscience and Remote Sensing，vol.45，no.3，pp.765-777，Mar.2007.

[14]D.C.Heinz and C-I Chang，“Fully constrained least squares linear spectralmixture analysis method for material quantification in hyperspectralimagery，”IEEE Transactions on Geoscience and Remote Sensing，vol.39，no.3，pp.529-545，Mar.2001.

[15]C.-I Chang and Q.Du，“Estimation of number of spectrally distinct signalsources in hyperspectral imagery，”IEEE Trans.Geosci.Remote Sens.，vol.42，no.3，pp.608-619，Mar.2004.

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

[17]R.N.Clark and G.A.Swayze，“Evolution in Imaging Spectroscopy Analysisand Sensor Signal-to-Noise：An Examination of How Far We Have Come，”The 6thAnnual JPL Airborne Earth Science Workshop，Mar.，1996.[Online].Available：http://speclab.cr.usgs.gov/PAPERS.imspec.evol/aviris.evolution.html

附录

附录A：性质1，2的证明

下面我们利用高阶累积量证明性质1和性质2。独立成分分析经常直接或间接地使用到高阶累积量[3]，设y是一个连续标量随机变量，其概率密度函数为p(y)，则y的矩生成函数(第一特征函数)为

其中，i²＝-1，每个概率分布由它的特征函数唯一决定，反之亦然.

y的累积量生成函数(第二特征函数)定义为

y的n阶累积量定义为

${\mathrm{\κ}}_n\left(y\right)={(-i)}^n\frac{d^n\mathrm{\φ}\left(w\right)}{{\mathrm{dw}}^n}{|}_{w=0}.---\left(38\right)$

本发明使用的概率模型是

根据(36)式和(37)式，可求出y_i的累积量生成函数为

$\mathrm{\φ}{\left(w\right)}_i=\ln \left\{(1-{\mathrm{\α}}_i)\exp \right[i{\mathrm{\μ}}_{1\left(i\right)}w-\frac{{\mathrm{\σ}}_i^2{w}^2}{2}]+{\mathrm{\α}}_i\exp [i{\mathrm{\μ}}_{2\left(i\right)}w-\frac{{\mathrm{\σ}}_i^2{w}^2}{2}\left]\right\}.---\left(39\right)$

代入(38)式，可求出y_i的前4个累积量分别是

κ₁(y_i)＝(1-α_i)μ_1(i)+α_iμ_2(i)                          (40)

κ₂(y_i)＝α_i(1-α_i)(μ_1(i)-μ_2(i))²+σ²，                    (41)

κ₃(y_i)＝α_i(1-α_i)(2α_i-1)(μ_1(i)-μ_2(i))³，                (42)

κ₄(y_i)＝α_i(1-α_i)(μ_1(i)-μ_2(i))⁴(6α_i ²-6α_i+1)，          (43)

将μ₂＝0代入(40)-(43)式，可得

κ₁(y_i)＝(1-α_i)μ_1(i)，                                     (44)

${\mathrm{\κ}}_2\left(y_i\right)={\mathrm{\α}}_i(1-{\mathrm{\α}}_i){\mathrm{\μ}}_{1\left(i\right)}^2+{\mathrm{\σ}}^2,---\left(45\right)$

${\mathrm{\κ}}_3\left(y_i\right)={\mathrm{\α}}_i(1-{\mathrm{\α}}_i)(2{\mathrm{\α}}_i-1){\mathrm{\μ}}_{1\left(i\right)}^3,---\left(46\right)$

${\mathrm{\κ}}_4\left(y_i\right)={\mathrm{\α}}_i(1-{\mathrm{\α}}_i){\mathrm{\μ}}_{1\left(i\right)}^4(6{{\mathrm{\α}}_i}^2-6{\mathrm{\α}}_i+1).---\left(47\right)$

所以我们有

a)当0＜α_i＜0.5时，κ₃(y_i)＜0；

b)当0.5＜α_i＜1时，κ₃(y_i)＞0；

c)当

或

时，κ₄(y_i)＞0；

d)当

时，κ₄(y_i)＜0。

信号y_i的均值定义为

E(y_i)＝mean(y_i)＝κ₁，                  (48)

偏度定义为

${\mathrm{\γ}}_1\left(y_i\right)=\frac{{\mathrm{\κ}}_3\left(y_i\right)}{{\mathrm{\κ}}_2^{3/2}\left(y_i\right)},---\left(49\right)$

峭度定义为

${\mathrm{\γ}}_2\left(y_i\right)=\frac{{\mathrm{\κ}}_4\left(y_i\right)}{{\mathrm{\κ}}_2^2\left(y_i\right)},---\left(50\right)$

当γ₁(y_i)＞0时，信号y_i为右偏，否则为左偏。

当γ₂(y_i)＞0时，信号y_i为超高斯型，否则为亚高斯型。

由(44)式和(48)式可得

mean(y_i)＝(1-α_i)μ_1(i)，                  (51)

因而性质2得证。

由于0＜α_i＜1，故

因此偏度γ₁(y_i)，峭度γ₂(y_i)的符号分别与κ₃(y_i)，κ₄(y_i)的符号相同。我们有

a)当0＜α_i＜0.5时，γ₁(y_i)＜0；

b)当0.5＜α_i＜1时，γ₁(y_i)＞0；

c)当

或

时，γ₂(y_i)＞0；

d)当

时，γ₂(y_i)＜0。

从而性质1得证.

附录B：σ_i取值的分析

丰度信号的物理意义决定了y_i∈[0，1]，因此我们希望p(0≤y_i≤1)应该足够大。

由(18)式知，p(y_i)由

和

两部分组成，其中，参数α_i和μ_1(i)是已知量，它们可以根据前面说明的调整规则确定(调整规则不涉及σ_i)。我们简化问题，假设已经求出α_i＝0.5，μ_1(i)＝1，这时正态分布

和

关于y_i＝0.5对称，让p(0≤y_i≤1)足够大就可以近似成对于正态变量

p(0≤z≤1)应该足够大。

正态变量的一个著名性质是

p(μ_2(i)-3σ_i＜z＜μ_2(i)+3σ_i)≈99.74％，                    (52)

即变量z的取值落在μ_2(i)-3σ_i＜y_i＜μ_2(i)+3σ_i内几乎是一定的，这就是所谓的“3σ法则”。又因为μ_2(i)＝0，所以p(-3σ_i＜z＜3σ_i)≈99.74％。正态分布z的概率密度函数关于z＝0偶对称，进一步有

p(0≤z≤3σ_i)≈50％，                    (53)

p(3σ_i≤z≤1)≈0.                        (54)

根据(53)式，要让p(0≤z≤1)尽可能大应有3σ_i≤1，即σ_i≤1/3。根据(54)式，z的取值几乎不会落在[0，1]的子区间[3σ_i，1]上。一般而言我们希望z在[0，1]的各个子区间都能有分布，因此σ_i不可以取的过小，应有σ_i≥1/3。

综合以上分析，取σ_i＝1/3，(i＝1，...，P)。

需要指出的是，虽然这里给出的σ_i＝1/3是近似值，但是σ_i精确与否不会对算法带来决定性的影响。因为推导a_i和μ_1(i)的调整规则无需考虑参数σ_i取值多少，对我们的算法而言σ_i只是一种辅助性的参数而非决定性参数。

Claims (1)

1.一种基于独立分量分析的高光谱遥感图像混合像元分解方法，其特征在于具体步骤如下：

已知观测矩阵

端元个数P，

步骤1、初始化解混矩阵W；

步骤2、开始迭代：

a)根据解混矩阵W的迭代公式更新W：

W←W+ηΔW，ΔW＝ΔW_AAM+η₁ΔW_ANC+η₂ΔW_ASC，    (1)

参数η₁，η₂为用于控制约束条件的权重，η为迭代步长，

b)根据ICA解混公式计算Y：

Y＝WX.               (2)

c)如果算法没有收敛，返回a)继续迭代；

步骤3、输出结果；

其中，

ΔW_AAM＝(I+dy^T-9yy^T)W.                    (3)

式中，d＝9μ₁·/[1-α·/(α-1)·*exp(-9μ₁·*y+4.5μ₁·*μ₁)]，μ₁＝mean(y)./(1-α)，α＝sign(γ₂)/8+0.7693，mean(y)＝[mean(y₁)，mean(y₂)，...，mean(y_p)]^T，sign是符号函数，γ₂＝[γ₂(y₁)，...，γ₂(y_k)，...，γ₂(y_P)]^T，

符号″.*″和″./″表示矩阵元素对应地相乘和相除；

式中，矩阵G的第i行j列位置处的元素为

式中，矩阵H的第i行j列位置处的元素为

。

Images