CN102054273A - 基于单形体三角分解的高光谱遥感图像混合像元分解方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，具体为一种基于单形体三角分解的高光谱遥感图像混合像元分解方法。本发明采用采用线性混合模型，具体步骤包括两部分：端元提取和丰度估计。该方法既是一种单形体类的几何方法，同时又建立在三角分解的代数原理之上。三角分解可采用Cholesky分解和QR分解，能够通过递归操作，在端元提取过程中提高端元的搜索效率。本发明可以有效地提取高光谱遥感数据中的端元，解决相应的混合像元分解问题。在基于高光谱遥感图像的高精度的地物分类以及地面目标的检测和识别方面具有特别重要的应用价值。

Description

基于单形体三角分解的高光谱遥感图像混合像元分解方法

技术领域

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

背景技术

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

大部分混合像元分解方法都是基于线性混合模型的[2]，其像元分解过程可分为两个步骤，端元提取和丰度解混。端元提取是指从遥感数据中求取端元光谱，丰度解混则是指确定各个端元在混合像元中所占的比例。大多数端元提取算法属于几何学方法，这类方法假设高光谱图像的所有数据点都位于一个单形体中，单形体的顶点分别对应各端元，因此求解端元光谱等效于寻找单形体的顶点。这类方法有N-FINDR[3]，Vertex Component Analysis(VCA)[4]，Simplex Growing Algorithm(SGA)[5]等。N-FINDR和SGA都是通过使单形体体积最大化搜索出端元，但是直接计算单形体的体积时会带来较高的运算复杂度，从而导致这两种算法消耗较多的运算时间。

发明内容

本发明的目的在于提供一种运算复杂度低、计算工作量少的高光谱遥感图像混合像元分解方法。

本发明提出的高光谱遥感图像混合像元分解方法，是一种基于矩阵三角分解思想的像元解混方法，该方法将基于单形体的体积运算转换为基于三角分解的乘积计算，从而降低了运算量。我们将该方法命名为基于三角分解的单形体体积分析方法(Simplex VolumeAnalysis based on Triangular Factorization，SVATF)。该方法既是一种单形体类的几何方法，同时又建立在三角分解的代数原理之上，其特点是：SVATF并不是单一的一种方法，而是一类框架，涵盖了不同的具体实现方式。为说明问题，本说明中给出了SVATF的两种具体实现，不论哪种SVATF，都能够通过递归操作，在端元提取过程中提高端元的搜索效率。在提取端元的同时，该类方法还可以估计端元个数。我们还提出了基于QR分解的丰度解混方法，可以同时完成丰度求取和端元光谱校正。实验结果表明，与现有方法相比，SVATF能够以低运算复杂度取得更优的解，而且具有结果唯一的优点，为混合像元分解问题提供了一种有效的解决手段。

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

1线性混合模型

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

可被表示为：

x＝As+n.    (1)

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

是端元光谱矩阵，其中

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

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

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

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

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

忽略误差项n，则从而(1)式可表示成矩阵形式

X＝AS.    (3)

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

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

2单形体体积

n维空间中有P个仿射独立的点e₀，e₁，...，e_P-1，设E＝[e₀，e₁，...，e_P-1]，以这些点为顶点构成的单形体体积为：

$V=\frac{1}{(P-1)!}\left|\det \left(\left[\begin{array}{c}{1}^T\\ E\end{array}\right]\right)\right|---\left(4\right)$

其中1＝[1，1，...，1]^T是全为1的列向量。向量α_i＝e_i-e₀称为该单形体的第i个支撑棱，所有α_i(i＝1，...，P-1)构成支撑棱矩阵：

$\tilde{A}=[{\mathrm{\α}}_1,{\mathrm{\α}}_2,...,{\mathrm{\α}}_{P-1}]---\left(5\right)$

则有

仞等变换不改变行列式的值，因此单形体的体积还可以表示为：

$V=\frac{1}{(P-1)!}\left|\det \left(\tilde{A}\right)\right|=\frac{1}{(P-1)!}{\left|\det \left({\tilde{A}}^T\tilde{A}\right)\right|}^{\frac{1}{2}}---\left(7\right)$

式(4)或式(7)都可以用于单形体体积V的求解。经典的端元提取算法N-FINDR、SGA使用式(4)来衡量V，显然，这里行列式det(*)中必须为方阵才可以进行计算，因此每次搜索都必须先对数据进行降维。这意味着：1.原始信息量的丢失；2.随着行列式维数的逐渐增多，搜索速度将会迅速变慢。而对于式(7)，该体积表达式中

方阵，无需降维处理，因此本发明所提出的算法都将基于式(7)展开。

3矩阵的三角分解

三角分解将一个矩阵分解为正交矩阵与三角矩阵之积(QR分解)，或者分解为一个下三角阵和一个上三角阵之积(LU分解)。矩阵

的QR分解和LU分解分别如式(8)和式(9)所示：

H＝QR    (8)

H＝LU    (9)

这里，Q是正交矩阵，R是上三角矩阵。L是单位下三角，U是上三角阵。完整的QR分解生成m×m大小的正交阵Q和m×n大小的上三角阵R，其中R的后m-n行全为0。大多数情况下，并不需要这些全0的部分，因此在分解时往往只求取Q的前n列和R的前n行。这种简化的QR分解称为经济QR分解，生成m×n大小的列正交阵Q和n×n大小的上三角阵R。若m≥n，则一定存在列正交的矩阵Q ∈R^m×n和上三角阵R∈R^m×n，使得H＝QR。

如果H是对称正定矩阵，将有L＝U^T，(9)式所示的分解进一步表示成：

H＝LL^T    (10)

上式称为Cholesky分解或平方根分解。与普通的LU分解相比，Cholesky分解的计算量大约可以减少一半。当H是对称正定的方阵时，Cholesky分解必然存在。

本发明的提供的高光谱遥感图像混合像元分解方法，具体步骤包括两部分：端元提取和丰度估计，分别描述如下：

1.端元提取

首先介绍端元提取方法

1.1体积公式

对于矩阵设

$Z={\tilde{A}}^T\tilde{A}=\left[\begin{array}{cccc}{\mathrm{\α}}_1^T\·{\mathrm{\α}}_1& {\mathrm{\α}}_1^T\·{\mathrm{\α}}_2& ...& {\mathrm{\α}}_1^T\·{\mathrm{\α}}_{P-1}\\ ...& ...& ...& ...\\ {\mathrm{\α}}_{P-1}^T\·{\mathrm{\α}}_1& {\mathrm{\α}}_{P-1}^T\·{\mathrm{\α}}_2& ...& {\mathrm{\α}}_{P-1}^T\·{\mathrm{\α}}_{P-1}\end{array}\right]=\left[\begin{array}{cccc}{\left|\right|{\mathrm{\α}}_1\left|\right|}^2& {\mathrm{\α}}_1^T\·{\mathrm{\α}}_2& ...& {\mathrm{\α}}_1^T\·{\mathrm{\α}}_{P-1}\\ ...& ...& ...& ...\\ {\mathrm{\α}}_{P-1}^T\·{\mathrm{\α}}_1& {\mathrm{\α}}_{P-1}^T\·{\mathrm{\α}}_2& ...& {\left|\right|{\mathrm{\α}}_{P-1}\left|\right|}^2\end{array}\right]---\left(11\right)$

则式(7)可以表示为：

$V=\frac{1}{(P-1)!}{\left|\det \left(Z\right)\right|}^{\frac{1}{2}}---\left(12\right)$

由于每次端元搜索都要计算单形体体积，如果直接按照(12)式求行列式，运算量会较大。但是如果事先将矩阵进行三角分解，就可以将(12)式的行列式计算转化成乘法计算。因为对于任意一个上三角矩阵L或下三角矩阵U，都有：

det(L)＝l₁₁·l₂₂.....l_PP，det(U)＝u₁₁·u₂₂.....u_PP，    (13)

其中，l_ij和u_ij分别表示矩阵L，U的第i行j列位置处的元素。

理论上，类似于N-FINDR，算法每提取一个端元就要遍历全部观测点进行搜索。对于N-FINDR而言是每次都求新的体积，对于SATF算法则是进行新的三角分解。但是，实际操作中SATF并不需要每次都进行一轮完整的三角分解计算。因为三角分解的公式一般都是递推形式，在不同的端元搜索步骤中，有些量是不变的，可以保留以备下次搜索使用，这是本算法的另一个优势。这一点将在后文的算法实施中进行具体说明。因此本发明可以总结为2点：1.通过最大化单形体体积来搜索端元。2.对矩阵进行三角分解，将(12)式的行列式计算转化成乘法计算。

由于矩阵三角分解的算法多种多样，每一种算法都可以分别衍生出不同的端元提取方法。本发明提出的基于三角分解的单形体分析方法(SATF)，具体选用有代表性的两种三角分解进行推导和阐述。这两种三角分解分别是Cholesky分解和QR分解。

1.2基于Cholesky分解的端元提取算法

显然矩阵Z是对称正定矩阵，所以可以对其使用Cholesky分解：

Z＝LL^T        (14)

因为L矩阵是上三角阵，所以：

|det(L)|＝|l₁₁|·|l₂₂|.....|l_(P-1)(P-1)|    (15)

单形体的体积公式(12)可以进一步表示为：

$V=\frac{1}{(P-1)!}{\left|\det \left({\mathrm{LL}}^T\right)\right|}^{\frac{1}{2}}=\frac{1}{(P-1)!}\left|l_{11}\right|\·\left|l_{22}\right|.....\left|l_{(P-1)(p-1)}\right|---\left(16\right)$

由式(16)可见，可以通过对P进行Cholesky分解，将单形体体积的求解从复杂的行列式计算的转换成简单的乘法计算。这时，最大化单形体体积等效为最大化|l_ii|，(i＝1，2，...，P)，从而使端元搜索效率大为提高。将一个矩阵H∈R^P×P进行Cholesky分解H＝LL^T的步骤如下：

for i＝1，2，....P-1，

for j＝i+1，...，P-1，

$l_{\mathrm{ii}}={(h_{\mathrm{ii}}-\underset{k=1}{\overset{i-1}{\mathrm{\Σ}}}{l}_{\mathrm{ik}}^2)}^{1/2}---\left(17\right)$

$l_{\mathrm{ji}}=(h_{\mathrm{ij}}-\underset{k=1}{\overset{i-1}{\mathrm{\Σ}}}{l}_{\mathrm{ik}}{l}_{\mathrm{jk}})/l_{\mathrm{ii}}---\left(18\right)$

将H＝Z代入式(17)(18)，得

$i=1\left\{\begin{array}{c}{l}_{11}=\sqrt{z_{11}}=\left|\right|{\mathrm{\α}}_1\left|\right|\\ l_{21}=z_{21}/l_{11},l_{31}=z_{31}/l_{11},...\end{array}\right.$

$i=2\left\{\begin{array}{c}{l}_{22}=\sqrt{z_{22}-l_{21}^2}=\sqrt{{\left|\right|{\mathrm{\α}}_2\left|\right|}^2-l_{21}^2}\\ l_{32}=(z_{32}-l_{31}{l}_{21})/l_{22},...\end{array}\right.$

$i=3\left\{\begin{array}{c}{l}_{33}=\sqrt{z_{33}-\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{l}_{3k}^2}=\sqrt{{\left|\right|{\mathrm{\α}}_3\left|\right|}^2-\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{l}_{3k}^2}=\sqrt{{\left|\right|{\mathrm{\α}}_3\left|\right|}^2-l_{31}^2-l_{32}^2}\\ l_{43}=l_{\mathrm{ji}}=(z_{43}-\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{l}_{3k}{l}_{4k})/l_{33},...\end{array}\right.---\left(19\right)$

...

由此我们可以得出以下结论：要从遥感图像的N个点里面找出能形成最大单形体体积的端元e₁，e₂，...，e_i，可以不直接计算单形体的体积，而是使用式(19)找出能使l_ii最大的点。

将以上结论应用于端元提取，从而可得出新的端元提取方法，具体流程如下：

算法 基于三角分解的遥感图像端元提取-Cholesky提取(SATF-CH)

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

步骤1初始化。对于所有像元矢量x_n，(n＝1，2，...，N)，

1a)搜索所有像元中矢量中模||x_n||最大的像元，作为端元矢量e₀。即：

$e_0=\underset{x_i}{\arg \max }\left(\right|\left|x_n\right|\left|\right)$

1b)令搜索距离e₀最远的像元作为e₁，即：

$e_1=\underset{x_i}{\arg \max }\left(\right|\left|{\tilde{x}}_n\right|\left|\right),$

并记端元e₁所在的位置为id(1)，即：

$\mathrm{id}\left(1\right)=\underset{n}{\arg \max }\left(\right|\left|{\tilde{x}}_n\right|\left|\right).$

1c)计算a₁＝e₁-e₀， ${\mathrm{\γ}}_1^n=\left|\right|{\tilde{x}}_n\left|\right|,$ ${\mathrm{\η}}_1^n=\left({\tilde{x}}_n^T{a}_1\right)/{\mathrm{\γ}}_1^{\mathrm{id}\left(1\right)}.$

步骤2设已提取的端元为e₀，e₁，...，e_t，这里1≤t＜P。对于所有像元矢量x_n，(n＝1，2，...，N)，

2a)计算 ${\mathrm{\γ}}_{t+1}^n=\sqrt{{\left({\mathrm{\γ}}_t^n\right)}^2-{\left({\mathrm{\η}}_t^n\right)}^2}.$

2b)搜索能够使

最大的e_t+1，即：

$e_{t+1}=\underset{x_n}{\arg \max }\left({\mathrm{\γ}}_{t+1}^n\right),$

并记e_t+1所在的位置为id(t+1)，即：

$\mathrm{id}(t+1)=\underset{n}{\arg \max }\left({\mathrm{\γ}}_{t+1}^n\right).$

2c)计算a_t+1＝e_t+1-e₀。

2d)计算 ${\mathrm{\η}}_{t+1}^n=({\tilde{x}}_n^T{a}_{t+1}-\underset{k=1}{\overset{t}{\mathrm{\Σ}}}{\mathrm{\η}}_k^n{\mathrm{\η}}_k^{\mathrm{id}(t+1)})/{\mathrm{\γ}}_{t+1}^{\mathrm{id}(t+1)}.$

2e)返回2a)继续迭代。

步骤3输出结果A＝[e₀，e₁，...，e_P-1]。

这里，步骤2a)和2d)分别表示在搜索第t+1个端元时，对第n个观测点分别执行Cholesky递推式(17)和式(18)。需要说明的是，Cholesky式(17)可以写成

$l_{\mathrm{ii}}={(h_{\mathrm{ii}}-\underset{k=1}{\overset{i-1}{\mathrm{\Σ}}}{l}_{\mathrm{ik}}^2)}^{1/2}={(h_{\mathrm{ii}}-\underset{k=1}{\overset{i-2}{\mathrm{\Σ}}}{l}_{\mathrm{ik}}^2-l_{i(i-1)}^2)}^{1/2}---\left(20\right)$

该式表示在每次搜索时更新当前观测点的l_ji，然而值得注意的是，(20)式中的

$h_{\mathrm{ii}}-\underset{k=1}{\overset{i-2}{\mathrm{\Σ}}}{l}_{\mathrm{ik}}^2---\left(21\right)$

其实就是当前的观测点在上一次搜索过程的

所以，每次搜索时我们用

将当前求得的l_ii保留，以便于在下一次搜索时直接使用。类似地还有式(18)，求l_ji时所使用的l_ik和l_jk也是以往搜索时求取的l_ji。于是可以在每次搜索结束时用

将算出的l_ji保留，以备以后使用。

1.2基于QR分解的端元提取算法

对于矩阵

若m≥n，则一定存在列正交的矩阵和上三角阵

使得H＝QR。支撑棱矩阵对于高光谱图像一般都有L≥P，故可以对矩阵

进行QR分解：

从而有：

${\tilde{A}}^T=R^T{Q}^T,$ $Z={\tilde{A}}^T\tilde{A}=R^T{Q}^T\mathrm{QR}---\left(22\right)$

因为Q矩阵是列正交的，所以：

Z＝R^TR    (23)

单形体的体积公式(12)可以进一步表示为：

$V=\frac{1}{(P-1)!}{\left|\det \left(R^{T}R\right)\right|}^{\frac{1}{2}}=\frac{1}{(P-1)!}\left|r_{11}\right|\·\left|r_{22}\right|....\left|r_{(P-1)(P-1)}\right|---\left(24\right)$

其中，r_ij表示矩阵R第i行j列位置处的元素。由式(24)可见，可以通过对

进行QR分解，将单形体体积的求解从复杂的行列式计算的转换成简单的乘法计算，即最大化单形体体积等效为最大化|r_ii|，(i＝1，2，...，P-1)，从而使端元搜索效率大为提高。QR分解的具体求法很多，如Gram-Schmidt正交化、household变换，Givens旋转等。其中比较简单易行的是Gram-Schmidt正交化，但是其数字性能不好，这里我们使用一种修正的Gram-Schmidt正交化方法，具体流程如下：

(1)令 ${\mathrm{\α}}_i^{\left(1\right)}={\mathrm{\α}}_i(i=\mathrm{1,2},...,P-1)$

for i＝1，2，...，P-2，

${\mathrm{\β}}_i={\mathrm{\α}}_i^{\left(i\right)},$ r_ii＝||β_i||，    (25)

q_i＝β_i/||β_i||，    (26)

$r_{\mathrm{ik}}=q_i^H{\mathrm{\α}}_k^{\left(i\right)},(k=i+1,...,P-1),---\left(27\right)$

${\mathrm{\α}}_k^{(i+1)}={\mathrm{\α}}_k^i-r_{\mathrm{ik}}{q}_i,(k=i+1,...,P-1),---\left(28\right)$

(2)计算 $r_{(p-1)(p-1)}=\left|\right|{\mathrm{\α}}_{(p-1)}^{(p-1)}\left|\right|,$ $q_{p-1}={\mathrm{\α}}_{p-1}^{(p-1)}/r_{(p-1)(p-1)}.$

让单形体体积最大化等效为使三角阵对角元r_ii(i＝1，2，...，P-1)(也即||β_i||)最大化。由QR分解的定义知这里β₁，β₂，...，β_P-1是两两正交的，称为正交基。将H＝Z代入式(25)～(28)，可得三角阵正交基β_i的表达式如下

i＝1    β₁＝α₁

$\begin{array}{cc}i=2& {\mathrm{\β}}_2={\mathrm{\α}}_2^{\left(2\right)}={\mathrm{\α}}_2-\frac{{\mathrm{\α}}_2\·{\mathrm{\β}}_1}{\left|\right|{\mathrm{\β}}_1\left|\right|}{\mathrm{\β}}_1\end{array}$

$\begin{array}{cc}i=3& \left\{\begin{array}{c}{\mathrm{\α}}_3^{\left(2\right)}={\mathrm{\α}}_3-\frac{{\mathrm{\α}}_3^T{\mathrm{\β}}_1}{\left|\right|{\mathrm{\β}}_1\left|\right|\·\left|\right|{\mathrm{\β}}_1\left|\right|}{\mathrm{\β}}_1\\ {\mathrm{\β}}_3={\mathrm{\α}}_3^{\left(3\right)}={\mathrm{\α}}_3^{\left(2\right)}-\frac{{\mathrm{\α}}_3^{\left(2\right)}\·{\mathrm{\β}}_2}{\left|\right|{\mathrm{\β}}_2\left|\right|\·\left|\right|{\mathrm{\β}}_2\left|\right|}{\mathrm{\β}}_2\end{array}\right.\end{array}$

…

即有：

${\mathrm{\α}}_{i+1}^{\left(i\right)}={\mathrm{\α}}_{i+1}^{(i-1)}-\frac{{\mathrm{\α}}_{i+1}^{(i-1)}\·{\mathrm{\β}}_{i-1}}{\left|\right|{\mathrm{\β}}_{i-1}\left|\right|\·\left|\right|{\mathrm{\β}}_{i-1}\left|\right|}{\mathrm{\β}}_{i-1}$ (29)

${\mathrm{\β}}_{i+1}={\mathrm{\α}}_{i+1}^{(i+1)}={\mathrm{\α}}_{i+1}^{\left(i\right)}-\frac{{\mathrm{\α}}_{i+1}^{\left(i\right)}\·{\mathrm{\β}}_1}{\left|\right|{\mathrm{\β}}_i\left|\right|\·\left|\right|{\mathrm{\β}}_i\left|\right|}{\mathrm{\β}}_i$

类似于Cholesky端元提取，我们可得新的端元提取方法，具体流程如下。

算法基于三角分解的遥感图像端元提取-QR提取(SATF-QR)

已知观测矩阵

端元个数P。

步骤1初始化。对于所有像元矢量x_n，(n＝1，2，...，N)，

1a)搜索所有像元中矢量中模||x_n||最大的像元，作为端元矢量e₀。即

$e_0=\underset{x_i}{\arg \max }\left(\right|\left|x_n\right|\left|\right).$

1b)令搜索距离e₀最远的像元作为e₁，即

$e_1=\underset{x_i}{\arg \max }\left(\right|\left|{\tilde{x}}_n\right|\left|\right),$

并记端元e₁所在的位置为id(1)，即

$\mathrm{id}\left(1\right)=\underset{n}{\arg \max }\left(\right|\left|{\tilde{x}}_n\right|\left|\right).$

1c)计算a₁＝e₁-e₀，

取β₁＝α₁。

步骤2设已提取的端元为e₀，e₁，...，e_t，(1≤t＜P)。对于所有像元矢量x_n，(n＝1，2，...，N)，

2a)计算 ${\mathrm{\θ}}_n^{t+1}={\mathrm{\θ}}_n^t-\frac{{\mathrm{\θ}}_n^t\·{\mathrm{\β}}_t}{{\mathrm{\β}}_t\·{\mathrm{\β}}_t}{\mathrm{\β}}_t,$

2b)搜索能够使最大的e_t+1，即：

$e_{t+1}=\underset{x_n}{\arg \max }\left(\right|\left|{\mathrm{\θ}}_n^{t+1}\right|\left|\right),$

记e_t+1所在的位置为id(t+1)，即：

$\mathrm{id}(t+1)=\underset{n}{\arg \max }\left(\right|\left|{\mathrm{\θ}}_n^{t+1}\right|\left|\right).$

同时记录下此时的即：

${\mathrm{\β}}_{t-1}=\underset{{\mathrm{\θ}}_n^{t+1}}{\arg \max }\left(\right|\left|{\mathrm{\θ}}_n^{t+1}\right|\left|\right)$

2c)返回2a)继续迭代。

步骤3输出结果：端元光谱A＝[e₀，e₁，...，e_P-1]和端元的位置索引[id(1)，id(2)，...，id(P)]。

该算法中，端元e₀，e₁，...，e_P-1是逐个依次找出的。搜索新的端元要遍历每个观测点，通过QR递推式(29)计算出对当前观测点所对应的三角阵正交基(对角元)，如果当前观测点的正交基是所有观测点中最大的，则该观测点就是新的端元。这里，我们通过步骤2a)实现式(29)的QR正交基计算，步骤2b)找出最大对角元，然后记录下新的端元。对于同一个观测点而言，由于正交基计算公式(29)是递推式，每一次端元搜索时计算该式所需要的量

β_i都已在上一次搜索时求出：是该观测点在上次端元搜索时所求出的正交基，β_i是上次端元搜索时找到的最大正交基。所以，我们在步骤2中分别以

和β_t+1将当前求得的

和β_i+1保留，以便于在下一次搜索时直接使用。

2丰度估计

在提取端元后，可以使用最小二乘方法来求解丰度信号。但是对于某些受噪声干扰、纯象元缺失等因素影响的数据，将导致算法端元提取有偏差，进而影响到丰度求取的精度。针对这一情况，本发明提出一种基于QR分解的丰度估计，在其中引入了线性混合模型的ANC约束和ASC约束。该方法并不是单纯的丰度估计方法，而是可以同时完成进行端元光谱信息校正和丰度估计，对于端元提取出现偏差的数据也能适用。在说明该方法前，我们首先介绍相关的代数原理。

2.1基于QR分解的线性系统求解

将光谱矩阵A做(经济)QR分解A＝QR，同时忽略误差e，则线性混合模型(1)可以写成如下形式：

x＝QRs    (30)

由于为正交矩阵，所以式(30)变形可以得到下式：

Q^Tx＝Rs

由式(31)可以看出，任意一点所对应的丰度信号s都是该式所示的线性方程组的解，而该方程组的求解非常简单，因为未知数s的系数矩阵R是上三角矩阵。在线性混合模型的矩阵形式(3)下，丰度信号的求解可以直接写成：

S＝inv(R)Q^TX    (32)

由X＝AS可得X^T＝S^TA^T，因而类似地，已知丰度矩阵S和观测矩阵X也可以借助QR分解来求解A。对S^T进行QR分解S^T＝Q_SR_S，然后得：

$A={\left(\mathrm{inv}\left(R_S\right)Q_S^T{X}^T\right)}^T---\left(33\right)$

2.2ASC约束条件

线性混合模型的ASC条件要求：对于任意像素点，所有端元地物在该点的丰度之和必须为1。为满足该条件，我们使用以下方法进行扩展：

$A_1\←\left[\begin{array}{c}A\\ 1^T\end{array}\right],$ $X_1\←\left[\begin{array}{c}X\\ 1^T\end{array}\right].---\left(34\right)$

其中1＝[1，1，...，1]^T是元素全为1的列矢量。

2.3算法实现

基于QR分解的丰度估计(QR-based Abundance Quantification Algorithm，QR-AQA)

步骤1、按式(34)进行扩展。

步骤2、执行QR分解A₁＝Q_AR_A，求取丰度矩阵S：

$S=\mathrm{inv}\left(R_A\right)Q_A^T{X}_1---\left(35\right)$

步骤3、找出丰度矩阵S中不满足ASC条件的元素

计算其在总点数中的比例，如果

$\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{\hat{s}}_{\mathrm{ij}}/N/P<\mathrm{\&zeta;}---\left(36\right)$

就视为算法收敛，跳转至0，其中ζ是一预定义的小量。

步骤4、对S进行限幅并进行并归一化

4a)S中所有大于1的元素置1，小于0的置0；

4b)对于s_j＝[s_1j，s_2j，...，s_Pj]^T，进行归一化：

${\tilde{s}}_j=\frac{s_j}{{\mathrm{\&Sigma;}}_{i=1}^P{s}_{\mathrm{ij}}},j=1,...,N---\left(37\right)$

将

按列排列，构造出矩阵

步骤5、执行QR分解

并校正光谱矩阵A₁：

$A_1=\mathrm{\&lambda;}{\left(\mathrm{inv}\left(R_S\right)Q_S^T{X}_1^T\right)}^T+(1-\mathrm{\&lambda;})A_1---\left(38\right)$

这里λ是遗忘因子，满足0＜λ＜1。

步骤6、将矩阵A₁和X₁的最后一行重新置零：

A₁(end，：)＝1^T，X₁(end，：)＝1^T    (39)

步骤7、转至步骤1继续迭代。

步骤8、输出丰度矩阵S和光谱矩阵A，A由A₁去掉其最后一行而得到。

显然我们直接根据式(31)就可以进行丰度估计，但是其估计精度会受端元矩阵A的精度影响。针对这一问题，本发明在每次迭代中都校正端元矩阵，并且把丰度的ASC和ANC都引入了求解过程。光谱矩阵的反馈校正在步骤5的式(38)执行，式(38)包括2部分，一是由在上一步中求出的丰度矩阵来反推光谱，根据式(33)引申得到；另一部分则是旧光谱矩阵A，因为无法保证由式(33)得到的光谱就是最优解，所以通过使用遗忘因子逐渐加入新信息削弱旧光谱，以防止校正过程出现剧烈波动。另外，针对丰度的约束条件，其中ANC约束要求丰度信号大于0小于1，本发明通过4a)限定信号必须在[0，1]范围内取值。对于ASC约束，则涉及了两种实现方法，一种是在步骤1中按照式(34)进行扩展，第二种则是在步骤4按式(37)中将丰度信号进行归一化。与式(34)不同之处在于，式(37)是一种启发式方法(heuristic approach)，其目的是通过让信号立即满足ASC以加快收敛速度。同时，式(37)的这种让丰度信号立即满足ASC的行为又不会影响算法的稳定性，因为我们的输出丰度信号是S而不是被强行归一化的

只是为调整端元光谱A而需要的一个中间量而已。在实际应用中，各种混合像元分解算法都无法保证求取的结果总是能够满足ANC或ASC，因此有些算法就采用了取绝对值或者类似于这里的步骤4的方法来对所求的结果做出强制改变，这样做的缺点是降低算法的稳定性，可能会对结果带来不可预测的影响，而本发明所提出的上述步骤则避免了这一点。

本发明的优点

本发明为一种新的基于三角分解算法的混合像元分解方法：Simplex Volume Analysisbased on Triangular Factorization，SVATF。其特点是通过对单形体体积公式进行三角分解，将求体积的行列式运算等效成简单的数乘运算。由于三角分解的方法多种多样，不同的方法将引申出不同的具体实现方式，所以SVATF本质上是涵盖了多种具体方法的算法框架。不论是哪种实现方式，三角分解几乎都以递推方式实现，因此本发明的端元搜索也是一种递归的过程，可以在每次搜索时可以利用上一次已获取的端元光谱信息，因而具有比现有方法更高的效率。此外，针对大多数基于单形体体积的方法需要纯象元假设的缺点，本发明还提出了一种新的基于QR分解的丰度估计方法。我们通过递归迭代，在求出丰度的同时对不精确的光谱信号进行矫正。SVATF是一种稳态算法，其解混方式决定了对于同一数据每次运行都能获得唯一的结果。仿真实验结果表明，所提出的算法对纯像元缺失和受到各种干扰的数据，都表现出良好的鲁棒性，以更低的运算耗时得出了更优的解。进一步地，本发明对于实际高光谱遥感数据实验也得到了理想的结果，证实了算法的有效性和对于各种不同数据的适用性。

附图说明：

图1 丰度信号图。

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

图3 纯象元缺失对算法的影响。其中，(a)RMSE，(b)SAD。

图4 Cuprite数据的伪彩色图。

图5 Cuprite数的丰度分解结果。其中，(a)白云母，(b)沙漠地表，(c)明矾石，(d)高岭石，(e)蒙脱石，(f)黄甲铁石，(g)铵长石，(h)高岭石#2，(i)皂石，(j)玉髓，(k)高岭石#3，(1)钙钛硅酸盐。

图6所提议算法提取端元光谱信号(虚线)与USGS库标准光谱信号(实线)对比。其中，(a)白云母，(b)沙漠地表，(c)明矾石，(d)高岭石，(e)蒙脱石，(f)黄甲铁石，(g)铵长石，(h)高岭石#2，(i)皂石，(j)玉髓，(k)高岭石#3，(1)钙钛硅酸盐。

具体实施方式

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

1仿真数据

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

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

之间的SAD定义为：

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

均方根误差(Root Mean Square Error，RMSE)是衡量丰度解混的指标。设第i个端元的真实丰度矢量为[s_i1，s_i2，...，s_iN]，它与其对应的解混结果

之间的RMSE定义为：

${\mathrm{RMSE}}_i=\sqrt{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{s}}_{\mathrm{ij}}-s_{\mathrm{ij}})}^2},(i=1,...P)---\left(41\right)$

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

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

通过两个实验评价本文所提出方法的有效性，并与其他算法进行比较。两个个实验分别测试噪声强度、纯像元缺失这两种情况对算法性能的影响。SVATF进行丰度估计时的参数统一设置为：遗忘因子λ＝0.1，终止条件ζ＝0.1。

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

实验2纯像元缺失情况下的适应性。纯像元缺失时能否得到良好的分解结果是衡量算法性能优劣的一个重要指标。本实验中，我们让纯度从1改变到0.5，比较各算法对纯像元缺失情况的适应性。图3(a)和(b)分别给出了在不同的纯像元缺失程度下，各算法的RMSE和SAD的变化情况。可见，随着像元纯度的降低，SVATF和VCA-FCLS的性能都在逐渐变差，但是SVATF在效果上始终优于VCA-FCLS。而MVCNMF和HOS-ICPA的结果则误差较大，RMSE和SAD两项指标都始终处于较高的水平。

2实际遥感数据

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

在算法进行解混前，去数据的除低信噪比波段、水吸收波段(包括波段1～2，104～113，148～167以及221～224)，使用剩下的188个波段进行实验。使用SVATF所得到的丰度解混结果如图5所示。与实地勘测地物分布0相比较，可以确定所求出的这些矿物，具体包括：白云母Muscovite#1，沙漠地表Desert_Varnish，明矾石Alunite，高岭石Kaolinite#1，蒙脱石Montmorillonite，黄甲铁石Jarosite，铵长石Buddingtonite，高岭石Kaolinite#2，皂石Nontronite，玉髓Chalcadony，高岭石Kaolinite#3，钙钛硅酸盐Sphene。为了进一步对目视判别结果进行验证，我们将所提取出的端元信号与美国标准数字光谱库(USGS)公布的对应矿物光谱进行了比对，结果如图6所示。其中，实线为标准地物光谱，虚线为所提议算法提取结果，可以看出，本方法提取的结果很好地吻合了标准地物光谱。

参考文献

[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，4^th ed.NewYork：John Wiley&Sons，Inc.，2000.

[3]M.E.Winter，“N-find：an algorithm for fast autonomous spectral endmemberdetermination in hyperspectral data，”in：Proc.of the SPIE conference on Imaging SpectrometryV，vol.3753，pp.266-275，1999.

[4]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，pp.898-910，April，2002.

[5]C-I Chang，C-C Wu，W.Liu，and Y-C Ouyang，“A new growing method forsimplex-based endmember extraction algorithm ，”IEEE Transactions on Geoscience and RemoteSensing，vol.44，no.10，pp.2804-2819，2006.

[6]J.Wang and C-I Chang，“Applications of independent component analysis inendmember extraction and abundance quantification for hyperspectral imagery，”IEEETransactions on Geoscience and Remote Sensing，vol.44，no.9，pp.2601-2616，Sep.2006.

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

[8]D.C.Heinz and C-I Chang，“Fully constrained least squares linear spectral mixtureanalysis method for material quantification in hyperspectral imagery，”IEEE Transactions onGeoscience and Remote Sensing，vol.39，no.3，pp.529-545，Mar.2001.

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 JPL AirborneEarth Science Workshop，Mar.，1996.[Online].Available：http://speclab.cr.usgs.gov/PAPERS.imspec.evol/aviris.evolution.html。

Claims (3)

1.一种基于单形体三角分解的高光谱遥感图像混合像元分解方法，采用线性混合模型，具体步骤包括两部分：端元提取和丰度估计，分别描述如下：

端元提取

n维空间中有P个仿射独立的点e₀，e₁，...，e_P-1，设E＝[e₀，e₁，...，e_P-1]，向量α_i＝e_i-e₀称为该单形体的第i个支撑棱，所有支撑棱构成矩阵

$\tilde{A}=[{\mathrm{\&alpha;}}_1,{\mathrm{\&alpha;}}_2,...,{\mathrm{\&alpha;}}_{P-1}]---\left(1\right)$

对于支撑棱矩阵设：

$Z={\tilde{A}}^T\tilde{A}=\left[\begin{array}{cccc}{\mathrm{\&alpha;}}_1^T\&CenterDot;{\mathrm{\&alpha;}}_1& {\mathrm{\&alpha;}}_1^T\&CenterDot;{\mathrm{\&alpha;}}_2& ...& {\mathrm{\&alpha;}}_1^T\&CenterDot;{\mathrm{\&alpha;}}_{P-1}\\ ...& ...& ...& ...\\ {\mathrm{\&alpha;}}_{P-1}^T\&CenterDot;{\mathrm{\&alpha;}}_1& {\mathrm{\&alpha;}}_{P-1}^T\&CenterDot;{\mathrm{\&alpha;}}_2& ...& {\mathrm{\&alpha;}}_{P-1}^T\&CenterDot;{\mathrm{\&alpha;}}_{P-1}\end{array}\right]=\left[\begin{array}{cccc}{\left|\right|{\mathrm{\&alpha;}}_1\left|\right|}^2& {\mathrm{\&alpha;}}_1^T\&CenterDot;{\mathrm{\&alpha;}}_2& ...& {\mathrm{\&alpha;}}_1^T\&CenterDot;{\mathrm{\&alpha;}}_{P-1}\\ ...& ...& ...& ...\\ {\mathrm{\&alpha;}}_{P-1}^T\&CenterDot;{\mathrm{\&alpha;}}_1& {\mathrm{\&alpha;}}_{P-1}^T\&CenterDot;{\mathrm{\&alpha;}}_2& ...& {\left|\right|{\mathrm{\&alpha;}}_{P-1}\left|\right|}^2\end{array}\right]---\left(2\right)$

则单形体的体积表示为：

$V=\frac{1}{(P-1)!}{\left|\det \left(Z\right)\right|}^{\frac{1}{2}}---\left(3\right)$

通过最大化单形体体积搜索端元，即对所有观测点计算其构成的单形体体积，搜索能形成最大体积的观测点作为端元；通过对体积矩阵Z或支撑棱矩阵

进行三角分解，将求体积时所用的(3)式的行列式计算转化成三角阵对角元的乘法计算；

丰度估计

在提取端元后，采用基于QR分解的丰度估计，其步骤为：

步骤1、按线性混合模型的ASC条件，对于任意像素点，所有端元地物在该点的丰度之和必须为1，为满足该条件，使用以下方法进行扩展：

$A_1\&LeftArrow;\left[\begin{array}{c}A\\ 1^T\end{array}\right],$ $X_1\&LeftArrow;\left[\begin{array}{c}X\\ 1^T\end{array}\right].---\left(4\right)$

其中1＝[1，1，...，1]^T是元素全为1的列矢量；

步骤2、执行QR分解A₁＝Q_AR_A，求取丰度矩阵S：

$S=\mathrm{inv}\left(R_A\right)Q_A^T{X}_1---\left(5\right)$

步骤3、找出丰度矩阵S中不满足ASC条件的元素计算其在总点数中的比例，如果：

$\underset{j}{\mathrm{\&Sigma;}}\underset{i}{\mathrm{\&Sigma;}}{\hat{s}}_{\mathrm{ij}}/N/P<\mathrm{\&zeta;}---\left(6\right)$

视为算法收敛，跳转至0，其中ζ是一预定义的小量；

步骤4、对S进行限幅并进行并归一化：

4a)S中所有大于1的元素置1，小于0的置0；

4b)对于s_j＝[s_1j，s_2j，...，s_Pj]^T，进行归一化：

${\tilde{s}}_j=\frac{s_j}{{\mathrm{\&Sigma;}}_{i=1}^P{s}_{\mathrm{ij}}},j=1,...,N---\left(7\right)$

将

按列排列，构造出矩阵

步骤5、执行QR分解

并校正光谱矩阵A₁：

$A_1=\mathrm{\&lambda;}{\left(\mathrm{inv}\left(R_S\right)Q_S^T{X}_1^T\right)}^T+(1-\mathrm{\&lambda;})A_1---\left(8\right)$

这里λ是遗忘因子，满足0＜λ＜1；

步骤6、将矩阵A₁和X₁的最后一行重新置零：

A₁(end，：)＝1^T，X₁(end，：)＝1^T    (9)

步骤7、转至步骤1继续迭代；

步骤8、输出丰度矩阵S和光谱矩阵A，A由A₁去掉其最后一行而得到。

2.根据权利要求1所述的基于单形体三角分解的高光谱遥感图像混合像元分解方法，其特征在于所述通过最大化单形体体积搜索端元，其中通过对体积矩阵Z进行三角分解，进行端元提取-Cholesky提取，具体步骤为：

已知观测矩阵

端元个数P；

步骤1、初始化，对于所有像元矢量x_n，(n＝1，2，...，N)，

1a)搜索所有像元中矢量中模||x_n||最大的像元，作为端元矢量e₀，即：

$e_0=\underset{x_i}{\arg \max }\left(\right|\left|x_n\right|\left|\right).$

1b)令

搜索距离e₀最远的像元作为e₁，即：

$e_1=\underset{x_i}{\arg \max }\left({\tilde{x}}_n\right),$

记端元e₁所在的位置为id(1)，即：

$\mathrm{id}\left(1\right)=\underset{n}{\arg \max }\left(\right|\left|{\tilde{x}}_n\right|\left|\right).$

1c)计算 $a_1=e_1-e_0,{\mathrm{\&gamma;}}_1^n=\left|\right|{\tilde{x}}_n\left|\right|,$ ${\mathrm{\&eta;}}_1^n=\left({\tilde{x}}_n^T{a}_1\right)/{\mathrm{\&gamma;}}_1^{\mathrm{id}\left(1\right)};$

步骤2、设已提取的端元为e₀，e₁，...，e_t，这里1≤t＜P，对于所有像元矢量x_n，(n＝1，2，..，N)，

2a)计算 ${\mathrm{\&gamma;}}_{t+1}^n=\sqrt{{\left({\mathrm{\&gamma;}}_t^n\right)}^2-{\left({\mathrm{\&eta;}}_t^n\right)}^2};$

2b)搜索能够使

最大的e_t+1，即：

$e_{t+1}=\underset{x_n}{\arg \max }\left({\mathrm{\&gamma;}}_{t+1}^n\right),$

记e_t+1所在的位置为id(t+1)，即：

$\mathrm{id}(t+1)=\underset{n}{\arg \max }\left({\mathrm{\&gamma;}}_{t+1}^n\right).$

2c)计算a_t+1＝e_t+1-e₀；

2d)计算 ${\mathrm{\&eta;}}_{t+1}^n=({\tilde{x}}_n^T{a}_{t+1}-\underset{k=1}{\overset{t}{\mathrm{\&Sigma;}}}{\mathrm{\&eta;}}_k^n{\mathrm{\&eta;}}_k^{\mathrm{id}(t+1)})/{\mathrm{\&gamma;}}_{t+1}^{\mathrm{id}(t+1)};$

2e)返回2a)继续迭代；

步骤3、输出结果A＝[e₀，e₁，...，e_P-1]。

3.根据权利要求1所述的基于单形体三角分解的高光谱遥感图像混合像元分解方法，其特征在于所述通过最大化单形体体积搜索端元，其中通过对支撑棱矩阵

进行三角分解，进行端元提取-QR提取，具体步骤为：

已知观测矩阵

端元个数P，

少骤1、仞始化，对于所有像元矢量x_n，(n＝1，2，...，N)，

1a)搜索所有像元中矢量中模||x_n||最大的像元，作为端元矢量e₀，即：

$e_0=\underset{x_i}{\arg \max }\left(\right|\left|x_n\right|\left|\right).$

1b)令

搜索距离e₀最远的像元作为e₁，即：

$e_1=\underset{x_i}{\arg \max }\left(\right|\left|{\tilde{x}}_n\right|\left|\right),$

记端元e₁所在的位置为id(1)，即：

$\mathrm{id}\left(1\right)=\underset{n}{\arg \max }\left(\right|\left|{\tilde{x}}_n\right|\left|\right).$

1c)计算a₁＝e₁-e₀，

取β₁＝α₁；

步骤2、设已提取的端元为e₀，e₁，...，e_t，(1≤t＜P)，对于所有像元矢x_n，(n＝1，2，...，N)，

2a)计算 ${\mathrm{\&theta;}}_n^{t+1}={\mathrm{\&theta;}}_n^t-\frac{{\mathrm{\&theta;}}_n^t\&CenterDot;{\mathrm{\&beta;}}_t}{{\mathrm{\&beta;}}_t\&CenterDot;{\mathrm{\&beta;}}_t}{\mathrm{\&beta;}}_t,$

2b)搜索能够使

最大的e_t+1，即：

$e_{t+1}=\underset{x_n}{\arg \max }\left(\right|\left|{\mathrm{\&theta;}}_n^{t+1}\right|\left|\right),$

记e_t+1所在的位置为id(t+1)，即：

$\mathrm{id}(t+1)=\underset{n}{\arg \max }\left(\right|\left|{\mathrm{\&theta;}}_n^{t+1}\right|\left|\right).$

同时记录下此时的即：

${\mathrm{\&beta;}}_{t+1}=\underset{{\mathrm{\&theta;}}_n^{t+1}}{\arg \max }\left(\right|\left|{\mathrm{\&theta;}}_n^{t+1}\right|\left|\right)$

2c)返回2a)继续迭代；

步骤3、输出结果：端元光谱A＝[e₀，e₁，...，e_P-1]和端元的位置索引[id(1)，id(2)，...，id(P)]。

Images