CN103258330A - 一种高光谱图像端元丰度的估计方法 - Google Patents

Abstract

本发明公开了一种高光谱图像中端元丰度值的估计方法，包括以下步骤：从图像中提取端元；选择混合像元点；线性分解得到对应丰度值；求取端元对应的归一化光谱特征值；在直角坐标系下描点；进行曲线拟合得到二次曲线表达式；通过归一化光谱特征值映射得到剩余点的丰度值；求取估计值与实际值的均方根误差RMSE并判断RMSE是否满足所求精度。由于本发明通过将光谱特征值与端元丰度之间建立一定的关系来快速预测丰度值，从而避免了对所有混合像元点进行线性分解才能得到对应丰度值的弊端，在实际应用过程中，只需要选取极少量分布均匀的像元点进行线性分解，就可以获得所有像元点中对应端元的丰度值，这样就有效地缩短了混合像元分解的时间。

Description

一种高光谱图像端元丰度的估计方法

技术领域

本发明涉及一种遥感图像处理技术，具体涉及一种高光谱图像端元丰度的估计方法。

背景技术

通过遥感技术获得的高光谱图像的波段数可以达到几十甚至上百个，所以数据量非常庞大，这对于其应用领域来说，有利有弊。因其既可以带来丰富的地物信息，又可以带来大量的冗余信息。而且由于遥感仪器和大气传输过程中的混合效应，导致混合像元普遍存在于图像中。因此需要对混合像元进行分解获得其基本的组成成分（即端元）及这些基本组成成分所占的比例（即丰度）。

目前在实际过程中，高光谱图像中端元对应丰度值的获取是通过光谱解混实现的，而光谱解混是通过使计算值与实际值之间的误差达到最小，即最小二乘法来实现的。该算法在实际计算过程中会涉及到大量矩阵的运算，若用在如此大的图像上的处理时间消耗将会很大。用在海洋溢油监测方面，传统的最小二乘法耗时巨大这一特点会降低溢油监测的时效性。本发明建立在线性光谱混合模型之上，利用光谱特征和像元丰度之间的关系，建立了一种快速估计端元对应丰度值的方法，在保证了解混精度的前提下提高了光谱解混的效率。

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

1.线性光谱混合模型

线性混合模型是目前应用最广泛的一种光谱混合模型，它假设各个光子只能“看到”一种物质，其物理意义是光谱仪瞬时场接收到的地面像元的反射率或光谱值是像元中各组分（即端元）反射率的面积比加权和。因此遥感图像中的每个像元可表示为：

$r=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{m}_i{\mathrm{\α}}_i+\mathrm{\ϵ}=\mathrm{M\α}+\mathrm{\ϵ}$

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_i=1$

0≤α_i≤1

其中，m_i是构成像元r的各个端元，M为像元矢量，α_i是各个端元对应的丰度值，α为丰度矢量，ε是误差矢量，N是波段数。显然构成每个像元的各个端元丰度值要满足和为一，且各自大于0小于1的条件。

2.最小二乘法

最小二乘法基于线性混合模型r＝Mα+ε。它的基本思想就是要通过使误差ε＝r-Mα达到最小来求取端元的丰度值，即对min{(r-Mα)^T(r-Mα)}寻求最优解。根据受到约束条件的不同，对应有无约束最小二乘法UCLS，部分约束（和为一约束）最小二乘法SCLS，非负约束最小二乘法NCLS和全约束最小二乘法FCLS。

①在没有任何约束条件下，无约束最小二乘法UCLS为：

${\widehat{\mathrm{\α}}}_{\mathrm{LS}}={\left(M^{T}M\right)}^{-1}{M}^{T}r$

无约束最小二乘法解出的丰度值往往不准确，且有可能为负值，因此不能代表地物的真实丰度。

②在

该条件下，部分约束（和为一约束）最小二乘法SCLS为：

${\widehat{\mathrm{\α}}}_{\mathrm{SCLS}}=(I-\frac{{\left(M^{T}M\right)}^{-1}{\mathrm{II}}^T}{I^T{\left(M^{T}M\right)}^{-1}I}){\widehat{\mathrm{\α}}}_{\mathrm{LS}}+\frac{{\left(M^{T}M\right)}^{-1}I}{I^T{\left(M^{T}M\right)}^{-1}I}$

其中，I为p×p的单位矩阵。

③在0≤α_i≤1的条件下，解混结果称为非负约束最小二乘法NCLS，考虑到该式子为不等式，引入一个p维正数常量c＝(c₁,c₂,...,c_p)来保证非负约束。由此构造拉格朗日函数J：

$J=\frac{1}{2}{(\mathrm{M\α}-r)}^T(\mathrm{M\α}-r)+\mathrm{\λ}(\mathrm{\α}-c)$

其中，α＝c。对该函数求偏导，

${\frac{\∂J}{\∂\mathrm{\α}}|}_{{\widehat{\mathrm{\α}}}_{\mathrm{NCLS}}}=0\⇒M^{T}M{\widehat{\mathrm{\α}}}_{\mathrm{NCLS}}-M^{T}r+\mathrm{\λ}=0$

得到两个迭代方程，即：

${\widehat{\mathrm{\α}}}_{\mathrm{NCLS}}={\left(M^{T}M\right)}^{-1}{M}^{T}r-{\left(M^{T}M\right)}^{-1}\mathrm{\λ}$

$={\widehat{\mathrm{\α}}}_{\mathrm{LS}}-{\left(M^{T}M\right)}^{-1}\mathrm{\λ}$

$\mathrm{\λ}=M^T(r-M{\widehat{\mathrm{\α}}}_{\mathrm{NCLS}})$

通过迭代^[3]得到非负约束最小二乘法NCLS的最优解和拉格朗日乘子向量λ＝(λ₁,λ₂,...,λ_p)^T。

④将两个约束条件同时考虑在内，即可得到全约束最小二乘法FCLS：

建立新的向量N和s，其中 $N=\left[\begin{array}{c}\mathrm{\δM}\\ l^T\end{array}\right],$ $s=\left[\begin{array}{c}\mathrm{\δr}\\ 1\end{array}\right],$

δ控制和为一的约束条件。将N和s分别代替全约束最小二乘法NCLS算法中的M和r，通过迭代得到的丰度矩阵即为全约束最小二乘法FCLS。

3.光谱角距离

光谱角(SAM，spectral angle mapper)，也被称为波谱角，具体是以所要求取的像元光谱值与样本参考光谱值之间的夹角来定义的。在一个N维的样本空间中，某两个影像像元之间的光谱角的数学表达公式可表示如下：

$\mathrm{SAM}={\cos }^{-1}\frac{A\·B}{\left|A\right|*\left|B\right|}={\cos }^{-1}\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i*B_i}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i^2}*\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{B}_i^2}}$

式中，A＝(A₁,A₂,…,A_N)和B＝(B₁,B₂,…,B_N)分别代表样本空间中的两个影像像元的光谱值，可以是两个混合像元，也可以是一个端元、一个混合像元，N为维数，也是波段数，α即为光谱角，其取值范围由α的余弦值的范围决定，变化范围集中在

代表了两个像元或光谱曲线之间的相似程度：同种像元之间的光谱角很小，甚至可以接近于0，故对应的余弦值接近于1；相反，不同种的像元则会使α值变的很大，因此对应的余弦值可以接近于0。

4.光谱信息散度

光谱信息散度(SID,Spectral Information Divergence)是用来表征两个像元光谱特征之间概率分布的不同的。设像元矢量X＝(x₁,x₂,…,x_L)^Τ对应的概率矢量为P＝(p₁,p₂,…,p_L)^Τ，Y＝(y₁,y₂,…,y_L)^Τ是另一个不同于X的像元矢量，其对应的概率矢量为Q＝(q₁,q₂…,q_L)^Τ，则

$P_j=x_j/\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{x}_l,q_j=y_j/\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{y}_l$

根据信息论的知识，可获得像元X和Y对应于波段B_l的自信息量为：

I_l(x)＝-logp_l    I_l(y)＝-logq_l

则Y关于X的相对熵和X关于Y的相对熵为：

$D\left(X\right|\left|Y\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l{D}_l\left(X\right|\left|Y\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l(I_l\left(Y\right)-I_l\left(X\right))=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l\log \left(\frac{p_l}{q_l}\right)$

$D\left(Y\right|\left|X\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l{D}_l\left(Y\right|\left|X\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l(I_l\left(X\right)-I_l\left(Y\right))=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l\log \left(\frac{q_l}{p_l}\right)$

因此光谱信息散度SID的定义公式为：

SID(X,Y)＝D(X||Y)+D(Y||X)

5.光谱相关系数

两像元之间的光谱相关系数(SCC,Spectral Correlation Coefficient)的定义公式为：

$\mathrm{SCC}(X,Y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(x_i-{\mathrm{\μ}}_x)*(y_i-{\mathrm{\μ}}_y)}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-{\mathrm{\μ}}_x)}^2}*\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(y_i-{\mathrm{\μ}}_y)}^2}}$

其中，x_i和y_i分别为像元X和Y在各个波段上对应的光谱值，μ_x和μ_y分别为其平均值，N是像元空间的维度，或者说是图像的波段数。

该式的取值范围介于[-1,1]之间。若两个光谱矢量反向成比例，则其相关系数为-1，此时称两者为完全负相关，反之，若两个光谱矢量相等或者正向成比例，则其相关系数为1，此时称两者为完全正相关。

6.混合像元分解的精度评价

6.1已知真实数据的精度评价

对于模拟合成的图像，由于事先知道各个混合像元中的真实丰度值，因此某个端元的丰度结果可以用估计出的结果和真实值之间的均方根误差RMSE(root mean square error)来衡量，其公式为：

$\mathrm{RMSE}=\sqrt{\frac{1}{n}\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{(\mathrm{\α}\left(p\right)-\widehat{\mathrm{\α}}\left(p\right))}^2}$

其中，α(p)和分别是对应于图像中第p个像元的实际丰度和估计丰度，n是像元点的个数。

6.2未知真实数据的精度评价

对于真实的遥感图像，由于事先不知道各个混合像元中的真实丰度值，因此其分解结果可以用真实光谱值和估计光谱值之间的重建误差RE(reconstruction error)来衡量，其公式为：

$\mathrm{RE}=\sqrt{\frac{1}{\mathrm{nL}}\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|\hat{Y}\left(p\right)-Y\left(p\right)\left|\right|}^2}$

其中，Y(p)是实际的光谱值向量，

是估计的光谱值向量，L是波段数，n是像元点的个数。

估计光谱值向量可以由估计出的各个端元丰度值与端元对应的各个波段上的光谱值相乘后各自相加得到，具体按照公式（19）求得：

$\hat{Y}\left(p\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{M}_j{\mathrm{\α}}_j$

其中，M_j为第j个端元的光谱值向量，α_j为该端元对应的丰度值，N为端元的个数。

从最小二乘法的计算公式可以看出，不论是无约束还是有约束的结果，均需要经过矩阵的多次运算，并且为了使得结果准确，一般采用全约束的最小二乘法FCLS，而该方法中结果的获取需要经过矩阵的多次迭代，这无疑会使得端元丰度值的获取经历过长的时间，无法适应在对实时性要求强的场合使用。

与本发明密切相关的参考文献如下：

[1]童庆禧，张兵，郑兰芬.高光谱遥感——原理、技术与应用[M].北京：高等教育出版社.2006-6.

[2]Rasmum Bro and Sijmen De Jong.A Fast Nonnegativity constrained LeastSquares Algorithm[J].Journal of Chemometrics,1997,11(5):393-401.

[3]Chein-I-Chang,D.C.Heinz.Constrained Subpixel Target Detection forRemotely Sensed Imagery[J].IEEE Transactions on Geoscience and Remote Sensing,2000,38(3):1144-1159.

[4]D.C.Heinz,Chein-I-Chang.Fully Constrained Least Squares Linear SpectralMixture Analysis Method for Material Quantification in Hyperspectral Imagery[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(3):529-544.

[5]Chein-I-Chang.An Information-Theoretic Approach to Spectral Variability,Similarity,and Discrimination for Hyperspectral Image Analysis[J].IEEETRANSACTIONS ON INFORMATION THEORY,2000.VOL.46,NO.5:1927-1932.

[6]A.Halimi,Y.Altmann,N.Dobigeon,and J.-Y.Tourneret,“Nonlinearunmixing of hyperspectral images using a generalized bilinear model,”IEEE Trans.Geosci.Remote Sensing,vol.49,NO.11:4153-4162,Nov.2011.

发明内容

为了解决现有技术存在的上述问题，本发明要设计一种能够在保证分解精度的前提下，提高混合像元分解效率的高光谱图像中端元丰度值的估计方法。

为了实现上述目的，本发明的技术方案如下：

一种高光谱图像中端元丰度值的估计方法，包括以下步骤：

A、从图像中提取端元

对于模拟图像，根据相关性系数公式（1），随机模拟几组相关性小的数据组成模拟端元e1，e2，…，en，其中两像元之间的光谱相关系SCC（X,Y）的定义公式为：

$\mathrm{SCC}(X,Y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(x_i-{\mathrm{\μ}}_x)*(y_i-{\mathrm{\μ}}_y)}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-{\mathrm{\μ}}_x)}^2}*\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(y_i-{\mathrm{\μ}}_y)}^2}}---\left(1\right)$

其中，x_i和y_i分别为像元X和Y在各个波段上对应的光谱值，μ_x和μ_y分别为其平均值，N是像元空间的维度，或者说是图像的波段数；

对于实际溢油图像，将油膜最厚的点即纯油区域选为石油端元，将纯海水区域选为海水端元；

B、选择k个混合像元点

对于模拟图像，按照随机的比例即随机的丰度值混合成k个混合像元，从而获得模拟混合图像；对于实际图像，选择图像中k个分布均匀的混合像元点；

C、线性分解得到对应丰度值

对于模拟图像，由于在混合成混合像元时各端元的比例是随机的，因此丰度值是事先知道的；对于实际高光谱图像，需要对选择出的混合像元点进行线性分解；

D、求取端元对应的归一化光谱特征值

所述的归一化光谱特征值包括归一化光谱角距离NSAM和归一化光谱信息散度NSID，归一化光谱角距离NSAM的计算公式为：

${\mathrm{NSAM}}_i=\frac{{\mathrm{SAM}}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{SAM}}_j}---\left(2\right)$

其中，NSAM_i在这里代表某个特定端元与第i个混合像元之间的归一化光谱角距离，SAM_i代表某个特定端元与第i个混合像元之间的光谱角距离，n代表所有混合像元点的个数，j的取值范围为1到n；在一个N维的样本空间中，某两个影像像元之间的光谱角的数学表达公式可表示如下：

$\mathrm{SAM}={\cos }^{-1}\frac{A\·B}{\left|A\right|*\left|B\right|}={\cos }^{-1}\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i*B_i}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i^2}*\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{B}_i^2}}---\left(2\right)$

式中，A＝(A₁,A₂,…,A_N)和B＝(B₁,B₂,…,B_N)分别代表样本空间中的两个影像像元的光谱值，可以是两个混合像元，也可以是一个端元、一个混合像元，N为维数，也是波段数，α即为光谱角；

归一化光谱信息散度NSID的计算公式为：

${\mathrm{NSID}}_i=\frac{{\mathrm{SID}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{SID}}_j}---\left(4\right)$

其中，NSID_i在这里代表某个特定端元与其中一个混合像元之间的归一化光谱信息散度，SID_i代表某个特定端元与其中一个混合像元之间的光谱信息散度，n代表所有混合像元点的个数；光谱信息散度SID的定义公式为：

SID(X,Y)＝D(X||Y)+D(Y||X)                               （5）

$P_j=x_j/\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{x}_l,q_j=y_j/\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{y}_l---\left(6\right)$

I_l(x)＝-logp_l    I_l(y)＝-logq_l                            （7）

$D\left(X\right|\left|Y\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l{D}_l\left(X\right|\left|Y\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l(I_l\left(Y\right)-I_l\left(X\right))=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l\log \left(\frac{p_l}{q_l}\right)---\left(8\right)$

$D\left(Y\right|\left|X\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l{D}_l\left(Y\right|\left|X\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l(I_l\left(X\right)-I_l\left(Y\right))=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l\log \left(\frac{q_l}{p_l}\right)---\left(9\right)$

其中，像元矢量X＝(x₁,x₂，…,x_L)^Τ对应的概率矢量为P＝(p₁,p₂，…,p_L)^Τ，Y＝(y₁,y₂,…,y_L)^Τ是另一个不同于X的像元矢量，其对应的概率矢量为Q＝(q₁,q₂,…,q_L)^Τ，I_l(x)和I^l(y)是像元X和Y对应于波段B_l的自信息量，D(Y||X)和D(X||Y)是Y关于X的相对熵和X关于Y的相对熵；

根据上述公式求取当前端元与合成的混合像元的归一化光谱角距离NSAM或归一化光谱信息散度NSID；

E、将归一化光谱特征值与对应的丰度值在直角坐标系下描点

对于模拟图像，由于在合成混合像元时各个端元所占的比例是事先知道的，即端元丰度值是已知的；而对于实际高光谱图像，在步骤C中已经线性分解得到对应的丰度值，因此根据步骤D计算出的该端元与所选混合像元之间的归一化光谱角距离NSAM及归一化光谱信息散度NSID在直角坐标系下描点：一是将归一化光谱角距离NSAM作为横坐标、丰度值作为纵坐标在直角坐标系上进行描点，得到NSAM-丰度值分布图；二是将归一化光谱信息散度NSID作为横坐标，丰度值作为纵坐标在直角坐标系上进行描点，得到NSID-丰度值分布图；

F、进行曲线拟合得到二次曲线表达式y＝ax²+bx+c

对NSAM-丰度值分布图和NSID-丰度值分布图分别采用最小二乘拟合法进行拟合，通过拟合，得到线性二次曲线，线性二次曲线的表达式分别为：

y₁＝a₁x₁ ²+b₁x₁+c₁                                          （10）

y₂＝a₂x₂ ²+b₂x₂+c₂                                         （11）

其中，x₁代表NSAM，y₁代表丰度，x₂代表NSID，y₂代表丰度；

G、通过归一化光谱特征值映射得到剩余点的丰度值

由拟合出的基于光谱角SAM和光谱信息散度SID的二次曲线表达式y₁＝a₁x₁ ²+b₁x₁+c₁和y₂＝a₂x₂ ²+b₂x₂+c₂，通过计算光谱特征值NSAM或NSID快速映射获取各个端元在混合图像中所有混合像元点中的估计丰度值；

H、求取估计值与实际值的均方根误差RMSE

将步骤G快速映射出的估计丰度值与实际的丰度值计算均方根误差RMSE，并测量丰度估计所需的时间，以估计其高效性；

对于模拟合成的图像，由于事先知道各个混合像元中的真实丰度值，因此某个端元的丰度结果用估计出的结果和真实值之间的均方根误差RMSE来衡量，其公式为：

$\mathrm{RMSE}=\sqrt{\frac{1}{n}\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{(\mathrm{\α}\left(p\right)-\widehat{\mathrm{\α}}\left(p\right))}^2}---\left(12\right)$

其中，α(p)和

分别是对应于图像中第p个像元的实际丰度和估计丰度，n是像元点的个数；

对于真实的遥感图像，由于事先不知道各个混合像元中的真实丰度值，因此其分解结果用真实光谱值和估计光谱值之间的重建误差RE来衡量，其公式为：

$\mathrm{RE}=\sqrt{\frac{1}{\mathrm{nL}}\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|\hat{Y}\left(p\right)-Y\left(p\right)\left|\right|}^2}---\left(13\right)$

其中，Y(p)是实际的光谱值向量，

是估计的光谱值向量，L是波段数，n是像元点的个数；

估计光谱值向量

由估计出的各个端元丰度值与端元对应的各个波段上的光谱值相乘后各自相加得到，具体按照公式（14）求得：

$\hat{Y}\left(p\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{M}_j{\mathrm{\α}}_j---\left(14\right)$

其中，M_j为第j个端元的光谱值向量，α_j为该端元对应的丰度值，N为端元的个数；

I、判断RMSE或RE是否满足所求精度

计算出均方根误差RMSE或RE之后，看所求精度是否满足要求，若不满足，则可能是拟合点选取的个数太少，或者是混合像元点选取的不好，那么就跳回步骤B，重新选择混合像元点，或增加混合像元点的个数，继续往下执行，直到所求均方根误差RMSE或RE满足所求精度为止。

与现有技术相比，本发明具有以下有益效果：

1、本发明命名为基于光谱特征的丰度快速估计方法AEMSC（AbundanceEstimation Method based on Spectral Characteristics）。由于光谱特征不是只包含这两种，因此本发明并不是单一的一种类型，而是一类框架，涵盖了不同的实现方式。为了说明本发明的有效性，本说明中给出了AEMSC的两种实现方式SAM和SID，在具体实施方式中将会看到，不论哪一种构成方式，其获取结果所需时间均很短，因此均能有效地提高丰度估计的效率。实验结果表明，与直接对遥感图像进行最小二乘混合像元分解相比，本发明能够在保证分解精度的前提下，大大提高混合像元分解的效率，为快速获取端元丰度提供了一种新方法。

2、本发明建立在线性光谱混合模型基础上，将光谱特征与丰度值之间建立起了关系，因此可以将丰度值的求取简单转化为对光谱特征值的求取。

3、由于本发明通过将光谱特征值与端元丰度之间建立一定的关系来快速预测丰度值，从而避免了对所有混合像元点进行线性分解才能得到对应丰度值的弊端，在实际应用过程中，只需要选取极少量分布均匀的像元点进行线性分解，就可以获得所有像元点中对应端元的丰度值，这样就有效地缩短了混合像元分解的时间。

4、由于本发明选取的光谱特征求取简单，因此大大降低了获取结果所需的时间。

附图说明

本发明共有附图27张，其中：

图1是本发明的处理流程图。

图2是模拟混合图像的第一波段示意图。

图3是模拟混合图像的第二波段示意图。

图4是模拟混合图像的第三波段示意图。

图5是端元1基于NSAM的丰度估计曲线。

图6是端元1基于NSID的丰度估计曲线。

图7是端元2基于NSAM的丰度估计曲线。

图8是端元2基于NSID的丰度估计曲线。

图9是端元3基于NSAM的丰度估计曲线。

图10是端元3基于NSID的丰度估计曲线。

图11是油膜模拟混合图像第一波段示意图。

图12是油膜模拟混合图像第二波段示意图。

图13是油膜模拟混合图像第三波段示意图。

图14是4个点的NSAM丰度估计曲线。

图15是6个点的NSAM丰度估计曲线。

图16是9个点的NSAM丰度估计曲线。

图17是13个点的NSAM丰度估计曲线。

图18是不同拟合点的RMSE变化趋势示意图。

图19是蓬莱溢油部分图像。

图20是实验用遥感图像。

图21是实验用遥感图像。

图22是对应于图20的NSAM丰度估计曲线。

图23是对应于图21的NSAM丰度估计曲线。

图24是对应于图21的最小二乘法石油丰度分布灰度图。

图25是对应于图21的AEMSC法石油丰度分布灰度图。

图26是对应于图22的最小二乘法石油丰度分布灰度图。

图27是对应于图22的AEMSC法石油丰度分布灰度图。

具体实施方式

下面结合附图对本发明进行进一步说明。下面分别根据模拟数据、实验数据和真实高光谱图像数据对本发明的具体步骤进行描述。

一、按照模拟数据建立AEMSC

如图1所示，一种高光谱图像中端元丰度值的估计方法，具体步骤如下：

A、根据相关性系数公式（1），随机模拟几组相关性小的数据组成模拟端元，其中，

B、再按照随机的比例即随机的丰度值混合成混合像元，从而获得模拟混合图像，如图2-4所示，图像的大小为12*12个像素，图像的四周分别为四个端元。

C、对于模拟图像，模拟丰度值事先知道，因此不需要对其进行线性分解。

D、根据公式（2）和公式（4）可以获得各个端元与混合像元的归一化光谱角距离NSAM和归一化光谱信息散度NSID。

E、将获得的归一化光谱特征值与丰度值在直角坐标系下描点。

F、采用线性二次最小二乘法拟合，获得的分别基于NSAM和NSID的AEMSC曲线分别如图5-10所示。其中：

图5是端元1基于NSAM的丰度估计曲线；图中，横坐标代表NSAM，纵坐标代表丰度，拟合出的曲线表达式为y＝4381.0x²-140.7198x+1.1448；

图6是端元1基于NSID的丰度估计曲线；图中，横坐标代表NSID，纵坐标代表丰度，拟合出的曲线表达式为y＝1284.5x²-64.0676x+0.6811；

图7是端元2基于NSAM的丰度估计曲线；图中，横坐标代表NSAM，纵坐标代表丰度，拟合出的曲线表达式为y＝4734.6x²-139.2751x+1.1094；

图8是端元2基于NSID的丰度估计曲线；图中，横坐标代表NSID，纵坐标代表丰度，拟合出的曲线表达式为y＝3075.2x²-96.0542x+0.8228；

图9是端元3基于NSAM的丰度估计曲线；图中，横坐标代表NSAM，纵坐标代表丰度，拟合出的曲线表达式为y＝2714.1x²-95.4297x+0.8344；

图10是端元3基于NSID的丰度估计曲线；图中，横坐标代表NSID，纵坐标代表丰度，拟合出的曲线表达式为y＝635.3557x²-40.1158x+0.5058；

G、分别通过基于NSAM和NSID的AEMSC得到原始模拟图中各个端元对应的丰度值。

H、按照公式（12）可以分别获得对应各个端元的丰度估计均方根误差如表1所示，由于4个端元对应像元点中丰度值的和为一，因此端元4的丰度值可以由1减去三个端元的丰度值而得到。其中，表中所测时间由曲线拟合时间加上丰度估计时间之和测得。

表1基于不同光谱特性的AEMSC方法的均方根误差和时间

由此表可以发现，对于基于NSAM和NSID的丰度估计方式，无论哪一种方式，其所需时间均很短，因此均能有效地提高丰度估计的效率。进一步对比发现：对于任意一个端元，基于NSAM的丰度估计方法的误差均要比基于NSID的方法小一些，因此其估计端元丰度的性能要更好。虽然基于NSID的丰度估计方法也能获得比较满意的结果，但为了获得丰度估计更加精确的结果，本发明以后均采用基于NSAM的丰度估计方法。

I、根据获得的各个端元的丰度估计均方根误差，可以看出该方法的估计结果均在所需要的精度范围内。

二、按照实验数据验证AEMSC

为进一步地验证本发明的鲁棒性，本发明设计了室内油膜光谱曲线采集实验，使用的仪器是Avantes公司生产的AvaSpec-2024光谱仪，该仪器的性能为：采集波长范围是300nm～1100nm，测量精度为0.6nm，在采集波长范围内可产生1335个采样，即1335个波段。实验时间为2012年3月下旬上午10点，在实验过程中分别在自来水和海水存在的前提下，倒入不同量的汽油和柴油来获取所需要的数据。在使用清洁自来水的实验过程中，由于该实验水槽为一直径为11cm的圆桶，在往水槽中滴入汽油的过程中，当滴入汽油的量达到200ml时就可以将该水槽平面盖满，因此可以认为此时获得的结果为不含有水的纯汽油的光谱曲线。具体验证步骤如下：

A、根据实验数据组成模拟端元并合成混合像元

在实验验证过程中，根据实验获取的纯油和纯水的光谱值组成相应的端元，这里，从汽油的量达到200ml时测得的光谱曲线和所测得的纯水光谱曲线中分别选取R、G、B三个波段700.57nm、546.63nm、435.88nm分别对应的反射率，得到含有三个波段的汽油和纯水的端元。

B、将模拟端元合成混合像元

按照随机的比例，即随机的丰度值混合成混合像元，最后也得到模拟混合图像，将测得的这两个端元按照随机的比例混合成20个混合像元，从而组成一幅模拟混合图像，如图11-13所示。

C、对于合成的图像，丰度值事先知道，因此不需要对其进行线性分解。

D、求取模拟端元与混合像元的归一化光谱特征值

由于基于光谱角SAM的丰度估计结果精度更高一些，因此为了提高精度，这里采用基于光谱角SAM的丰度估计方式。

根据公式（2）求取当前端元与合成的混合像元的归一化光谱角距离NSAM。

E、将归一化光谱特征值和其相应的端元丰度值在坐标轴上描点

由于在利用实验获取的值合成混合像元时各个端元所占的比例是事先知道的，即端元丰度值是已知的，步骤D已计算出该端元与所选混合像元之间的归一化光谱角距离NSAM，选取合适的点，将其归一化光谱角距离NSAM作为横坐标，丰度值作为纵坐标在直角坐标系上进行描点。

F、利用二次曲线拟合得到基于光谱角SAM的丰度估计二次曲线表达式

对归一化光谱角距离NSAM和丰度值在坐标轴上描出的点，采用最小二乘拟合法进行拟合，通过拟合，得出的曲线也为线性二次曲线，因此利用拟合出的二次曲线表达式，就可以通过计算光谱特征值来快速映射得到丰度值。

G、通过像元间光谱特征值的简单求取，快速预测端元对应的丰度值

由拟合出的基于光谱角SAM的二次曲线表达式，快速映射来获取各个端元在实验模拟混合图像中剩余混合像元点中的模拟丰度值。

H、计算估计出的丰度值与实际模拟丰度值的RMSE

利用公式（12），将步骤G快速映射出的丰度值与步骤C中实际的模拟丰度值计算均方根误差RMSE，评估基于光谱角SAM的丰度估计结果。

I、再次选取不同模拟混合像元点重复E、F、G步骤，计算其快速映射出的丰度值与实际的模拟丰度值的均方根误差RMSE，得出不同数量的混合像元点的选取对应的均方根误差变化曲线。在实际使用中可以根据所能接受的误差范围合理选择一定数量的混合像元点，即拟合点。

分别选取不同个数的像元点应用AEMSC方法，得到的丰度估计方式及RMSE变化趋势分别如图14-18所示。其中：

图14是4个点的NSAM丰度估计曲线；其中，横坐标代表NSAM，纵坐标代表丰度，拟合出的曲线表达式为y＝43.747x²-14.766x+0.97888。

图15是6个点的NSAM丰度估计曲线，其中，横坐标代表NSAM，纵坐标代表丰度，拟合出的曲线表达式为y＝41.58x²-14.604x+0.97649。

图16是9个点的NSAM丰度估计曲线，其中，横坐标代表NSAM，纵坐标代表丰度，拟合出的曲线表达式为y＝41.629x²-14.602x+0.97621。

图17是13个点的NSAM丰度估计曲线，其中，横坐标代表NSAM，纵坐标代表丰度，拟合出的曲线表达式为y＝41.698x²-14.601x+0.97609。

图18是不同拟合点的RMSE变化趋势，其中，横坐标代表拟合点个数，纵坐标代表RMSE；

图14-17中AEMSC的基本走势是一定的，只是在拟合曲线的系数上有微量的差别，因此该方法是固定的。图18表明随着拟合点个数的增加，误差值逐渐减小，但是降到一定值之后就会变得平缓，因此在实际应用过程中可以根据实际需要选取合适的拟合点个数。

三、根据真实高光谱图像数据验证AEMSC

本实验采用的数据来自于蓬莱-3C平台溢油遥感图像，如图19所示，该图像的总波段数目为258，每个波段的图像大小为3904*521，其中部分图像如图20-21所示。对其进行对数残差预处理后选取实验所用的图像，图20：大小为189*150，选取三个波段700.9600nm,500.3900nm,439.9400nm；图21：大小为258*173，选取三个波段703.3100nm,500.3900nm,442.1800nm。具体验证步骤如下：

A、从高光谱图像中选择端元

为了快速得到估计结果，对于溢油图像，可以粗略在油膜最厚的区域选择石油端元，而在纯海水区域选择海水端元；

B、选择高光谱图像中分布均匀的几个典型像元

对于实际高光谱图像，由于事先不知道混合像元点中对应端元的丰度，需要获取少量混合像元点中对应端元的丰度值，为了使拟合曲线对整幅图像都具有普遍性，则混合像元点的选取需要选择那些在图像中分布均匀的若干个典型像元，以提高丰度值的估计精度。为了提高精度，图20中选取了16个分布均匀的典型混合像元进行了线性分解；同理，图21中选取了12个分布均匀的典型混合像元进行了线性分解。

C、利用线性混合像元分解的方法计算每个像元中各端元的丰度值

选择好分布均匀的典型混合像元后，采用线性混合像元分解的方法对选择好的混合像元点进行分解，具体采用全约束最小二乘法FCLS。

D、求取所选端元与混合像元的归一化光谱特征值

求取当前端元与所选的若干个分布均匀的混合像元的归一化光谱角距离，得到基于SAM的归一化光谱角距离NSAM。

E、将归一化光谱角距离NSAM和相应的端元丰度值在坐标轴上描点；

通过线性分解获得了端元在所选若干个混合像元中的丰度值，并且也获得了该端元与所选混合像元之间的归一化光谱角距离NSAM，将归一化光谱角距离NSAM和对应的丰度值在直角坐标系上进行描点。

F、利用高次曲线拟合得到基于SAM的丰度估计二次曲线表达式

对坐标轴上描出的点进行曲线拟合，采用最小二乘拟合，通过拟合，得出的曲线也为线性二次曲线，因此利用拟合出的二次曲线表达式，就可以通过计算光谱特征值来快速映射得到丰度值；对于图20，拟合出的NSAM丰度估计曲线如图22所示，其中，横坐标表示归一化光谱角距离NSAM，纵坐标代表的是石油端元对应的丰度值，拟合出的二次曲线方程为y＝2574700x²-8879.9x+1.0054。对于图21，拟合出的NSAM丰度估计曲线如图23所示，其中，横坐标表示归一化光谱角距离NSAM，纵坐标代表的是石油端元对应的丰度值，拟合出的二次曲线方程为y＝13557000x²-17185x+0.98799。

G、通过像元间光谱值的简单求取，快速预测端元对应的丰度值

由拟合出的基于光谱角SAM的二次曲线表达式，快速映射来获取各个端元在图像中所有混合像元点中的丰度估计值。对于图20，AEMSC方法与最小二乘法得到的石油丰度分布图分别如图24-25所示；对于图21，AEMSC方法与最小二乘法得到的石油丰度分布图分别如图26和27所示。其中黑色的区域为原油覆盖的区域，黑色越深代表原油油膜厚度越深，即原油端元在该像元中所占的丰度值越大；白色区域为海水所在的区域，颜色越白代表海水的纯度越高，其它的灰色代表的区域为原油和海水以不同比例混合而成的区域。

H、求取实际光谱值和估计光谱值的重构误差RE

对于真实的遥感图像，由于事先不知道各个混合像元中的真实丰度值，因此可以利用估计出的丰度值与端元光谱值反演出估计光谱值，这样就可以求取实际光谱值和估计光谱值的重构误差RE。

I、如果该误差不在允许范围内，可以重新选择混合像元点，从步骤B开始重复，一直到计算出的误差在允许范围内。

利用公式（13）计算其重建误差RE，并测量所消耗的时间分别如表2所示。对于图20，可以看到本发明的AEMSC方法获得的石油丰度估计值的重建误差RE和全约束最小二乘法FCLS的解混结果的重建误差RE是在一个数量级上的，都保持在0.6%左右，因此可以说丰度估计模型AEMSC保持了原有算法的精度，但是在结果获取的时间上却减少了，减少的时间比例约为8.82%，有效地提高了丰度结果获取的效率，为高光谱图像中端元丰度值的快速获取提供了一种新的方法。对于图21，本发明获得的石油丰度估计值的重建误差RE和全约束最小二乘法FCLS的解混结果的重建误差RE是在一个数量级上的，都保持在0.7%左右，因此可以说本发明的丰度估计模型保持了原有算法的精度，但是在结果获取的时间上却减少了，减少的时间比例约为33.55%，有效地提高了丰度结果获取的效率，再次证明了本发明为高光谱图像中端元丰度值的快速获取提供了一种新方法的结论的正确性。

表2  AEMSC方法与最小二乘法的RE与时间比较

以上是分别根据模拟数据、实验数据和真实高光谱图像数据对本发明的具体步骤进行的描述。首先，本发明通过计算归一化光谱角距离NSAM和归一化光谱信息散度NSID来完成方法的搭建。利用合成的模拟数据，先计算NSAM和NSID值，这样就可以将其和对应的丰度值在坐标轴上进行描点，然后根据高次曲线拟合方式对两者的关系进行拟合，这样就获得了基于光谱特征的丰度估计模型AEMSC。最后所测得的时间以及所求的精度均在所需范围内，因此证明了本发明的高效性，并在此基础上分析了两种丰度估计方法的优劣。

其次，为了测试和验证所得方法的正确性，本发明进行了油膜光谱曲线采集实验。通过实验获取的数据合成了混合像元，然后选取不同数量的拟合样本进行丰度估计方法的搭建，每一次搭建都得到一个误差值，这样就获得了不同拟合样本所对应的误差变化曲线。分析每次所求得的误差值，再次证明了本发明能够保证所求结果的精度。

最后，本发明应用在蓬莱19-3C平台实际高光谱溢油图像上。完成对其的预处理之后，只需要将其中的一小部分混合像元点进行分解就可以完成丰度估计方法的搭建。利用本发明可以获得图像中所有混合像元点中对应石油的丰度值，然后将其获得的结果与直接用最小二乘法的计算结果进行对比。

模拟和实验以及实际图像均表明，本发明提出的方法在分解混合像元的效率上优于对图像直接采用最小二乘法，特别地，在效率上可以提高33.55%，因此可以满足实时结果的需求。

Claims (1)

1.一种高光谱图像中端元丰度值的估计方法，其特征在于：包括以下步骤：

A、从图像中提取端元

对于模拟图像，根据相关性系数公式（1），随机模拟几组相关性小的数据组成模拟端元e1，e2，…，en，其中两像元之间的光谱相关系SCC（X,Y）的定义公式为：

$\mathrm{SCC}(X,Y)=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(x_i-{\mathrm{\μ}}_x)*(y_i-{\mathrm{\μ}}_y)}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(x_i-{\mathrm{\μ}}_x)}^2}*\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{(y_i-{\mathrm{\μ}}_y)}^2}}---\left(1\right)$

其中，x_i和y_i分别为像元X和Y在各个波段上对应的光谱值，μ_x和μ_y分别为其平均值，N是像元空间的维度，或者说是图像的波段数；

对于实际溢油图像，将油膜最厚的点即纯油区域选为石油端元，将纯海水区域选为海水端元；

B、选择k个混合像元点

对于模拟图像，按照随机的比例即随机的丰度值混合成k个混合像元，从而获得模拟混合图像；对于实际图像，选择图像中k个分布均匀的混合像元点；

C、线性分解得到对应丰度值

对于模拟图像，由于在混合成混合像元时各端元的比例是随机的，因此丰度值是事先知道的；对于实际高光谱图像，需要对选择出的混合像元点进行线性分解；

D、求取端元对应的归一化光谱特征值

所述的归一化光谱特征值包括归一化光谱角距离NSAM和归一化光谱信息散度NSID，归一化光谱角距离NSAM的计算公式为：

${\mathrm{NSAM}}_i=\frac{{\mathrm{SAM}}_i}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{SAM}}_j}---\left(2\right)$

其中，NSAM_i在这里代表某个特定端元与第i个混合像元之间的归一化光谱角距离，SAM_i代表某个特定端元与第i个混合像元之间的光谱角距离，n代表所有混合像元点的个数，j的取值范围为1到n；在一个N维的样本空间中，某两个影像像元之间的光谱角的数学表达公式可表示如下：

$\mathrm{SAM}={\cos }^{-1}\frac{A\·B}{\left|A\right|*\left|B\right|}={\cos }^{-1}\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i*B_i}{\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{A}_i^2}*\sqrt{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{B}_i^2}}---\left(3\right)$

式中，A＝(A₁,A₂，…,A_N)和B＝(B₁,B₂,…,B_N)分别代表样本空间中的两个影像像元的光谱值，可以是两个混合像元，也可以是一个端元、一个混合像元，N为维数，也是波段数，α即为光谱角；

归一化光谱信息散度NSID的计算公式为：

${\mathrm{NSID}}_i=\frac{{\mathrm{SID}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{SID}}_j}---\left(4\right)$

其中，NSID_i在这里代表某个特定端元与其中一个混合像元之间的归一化光谱信息散度，SID_i代表某个特定端元与其中一个混合像元之间的光谱信息散度，n代表所有混合像元点的个数；光谱信息散度SID的定义公式为：

SID(X,Y)＝D(X||Y)+D(Y||X)                                     （5）

$P_j=x_j/\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{x}_l,q_j=y_j/\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{y}_l---\left(6\right)$

I_l(x)＝-logp_l    I_l(y)＝-logq_l                              （7）

$D\left(X\right|\left|Y\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l{D}_l\left(X\right|\left|Y\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l(I_l\left(Y\right)-I_l\left(X\right))=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{p}_l\log \left(\frac{p_l}{q_l}\right)---\left(8\right)$

$D\left(Y\right|\left|X\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l{D}_l\left(Y\right|\left|X\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l(I_l\left(X\right)-I_l\left(Y\right))=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{q}_l\log \left(\frac{q_l}{p_l}\right)---\left(9\right)$

其中，像元矢量X＝(x₁,x₂,…,x_L)^Τ对应的概率矢量为P＝(p₁,p₂,…,p_L)^Τ，Y＝(y₁,y₂,…,y_L)^Τ是另一个不同于X的像元矢量，其对应的概率矢量为Q＝(q₁,q₂,…,q_L)^Τ，I_l(x)和I_l(y)是像元X和Y对应于波段Bl的自信息量，D(Y||X)和D(X||Y)是Y关于X的相对熵和X关于Y的相对熵；

根据上述公式求取当前端元与合成的混合像元的归一化光谱角距离NSAM或归一化光谱信息散度NSID；

E、将归一化光谱特征值与对应的丰度值在直角坐标系下描点

对于模拟图像，由于在合成混合像元时各个端元所占的比例是事先知道的，即端元丰度值是已知的；而对于实际高光谱图像，在步骤C中已经线性分解得到对应的丰度值，因此根据步骤D计算出的该端元与所选混合像元之间的归一化光谱角距离NSAM及归一化光谱信息散度NSID在直角坐标系下描点：一是将归一化光谱角距离NSAM作为横坐标、丰度值作为纵坐标在直角坐标系上进行描点，得到NSAM-丰度值分布图；二是将归一化光谱信息散度NSID作为横坐标，丰度值作为纵坐标在直角坐标系上进行描点，得到NSID-丰度值分布图；

F、进行曲线拟合得到二次曲线表达式y＝ax²+bx+c

对NSAM-丰度值分布图和NSID-丰度值分布图分别采用最小二乘拟合法进行拟合，通过拟合，得到线性二次曲线，线性二次曲线的表达式分别为：

y₁＝a₁x₁ ²+b₁x₁+c₁                                       （10）

y₂＝a₂x₂ ²+b₂x₂+c₂                                       （11）

其中，x₁代表NSAM，y₁代表丰度，x₂代表NSID，y₂代表丰度；

G、通过归一化光谱特征值映射得到剩余点的丰度值

由拟合出的基于光谱角SAM和光谱信息散度SID的二次曲线表达式y₁＝a₁x₁ ²+b₁x₁+c₁和y₂＝a₂x₂ ²+b₂x₂+c₂，通过计算光谱特征值NSAM或NSID快速映射获取各个端元在混合图像中所有混合像元点中的估计丰度值；

H、求取估计值与实际值的均方根误差RMSE

将步骤G快速映射出的估计丰度值与实际的丰度值计算均方根误差RMSE，并测量丰度估计所需的时间，以估计其高效性；

对于模拟合成的图像，由于事先知道各个混合像元中的真实丰度值，因此某个端元的丰度结果用估计出的结果和真实值之间的均方根误差RMSE来衡量，其公式为：

$\mathrm{RMSE}=\sqrt{\frac{1}{n}\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{(\mathrm{\α}\left(p\right)-\widehat{\mathrm{\α}}\left(p\right))}^2}---\left(12\right)$

其中，α(p)和

分别是对应于图像中第p个像元的实际丰度和估计丰度，n是像元点的个数；

对于真实的遥感图像，由于事先不知道各个混合像元中的真实丰度值，因此其分解结果用真实光谱值和估计光谱值之间的重建误差RE来衡量，其公式为：

$\mathrm{RE}=\sqrt{\frac{1}{\mathrm{nL}}\underset{p=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|\hat{Y}\left(p\right)-Y\left(p\right)\left|\right|}^2}---\left(13\right)$

其中，Y(p)是实际的光谱值向量，

是估计的光谱值向量，L是波段数，n是像元点的个数；

估计光谱值向量

由估计出的各个端元丰度值与端元对应的各个波段上的光谱值相乘后各自相加得到，具体按照公式（14）求得：

$\hat{Y}\left(p\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{M}_j{\mathrm{\α}}_j---\left(14\right)$

其中，M_j为第j个端元的光谱值向量，α_j为该端元对应的丰度值，N为端元的个数；

I、判断RMSE或RE是否满足所求精度

计算出均方根误差RMSE或RE之后，看所求精度是否满足要求，若不满足，则可能是拟合点选取的个数太少，或者是混合像元点选取的不好，那么就跳回步骤B，重新选择混合像元点，或增加混合像元点的个数，继续往下执行，直到所求均方根误差RMSE或RE满足所求精度为止。

Images