CN109785305A - 一种端元可变的高光谱图像光谱混合分析方法及装置 - Google Patents

Abstract

本发明涉及高光谱遥感图像处理领域，公开了一种端元可变的高光谱图像光谱混合分析方法。首先利用微分光谱计算光谱形状稳定性指数，选择波段子集，求各类端元特征光谱；然后，基于类端元特征光谱对像元进行光谱混合分析预确定像元内地物类别；最后，基于枚举法或迭代光谱混合分析，利用像元内所含地物的可变端元法对各混合像元进行精细光谱混合分析。像元内地物类别预确定后结合光谱混合分析精度比较高的枚举法或迭代光谱混合分析法，可以降低光谱混合分析总体时间复杂度，同时提高其精度。因此，本发明提出的端元可变高光谱图像光谱混合分析方法能够提高光谱混合分析的效率。本发明还提供了一种端元可变的高光谱图像光谱混合分析装置。

Description

一种端元可变的高光谱图像光谱混合分析方法及装置

技术领域

本发明涉及图像处理领域，尤其涉及一种端元可变的高光谱图像光谱混合分析方法及装置。

背景技术

随着定量化研究的不断深入，高光谱遥感技术以其光谱分辨率较高的特点受到国内外的广泛关注，在地质制图、植被调查、海洋遥感、农业遥感、大气研究和环境监 测，以及城市应用等方面的研究中取得了引人注目的成果，已经成为当前遥感的一个 重要发展方向。由于空间分辨率的限制，高光谱图像中混合像元普遍存在，光谱混合 分析作为高光谱的一种预处理手段，是实现地物精确分类和识别的重要前提，对发挥 高光谱图像在社会经济和国防建设方面的效益有着重要的现实意义。

国内外众多学者对光谱混合分析做了大量研究，提出了各种有效的方法，其中大多数方法都是基于固定的端元集进行光谱混合分析，即图像中每种地物的纯像元光谱 是唯一的。由于地物的复杂多样性和成像条件的影响，高光谱图像“同物异谱”和“异 物同谱”现象普遍存在，同一种地物的纯像元光谱不一定完全相同，导致对所有像元 用固定的端元集解混精度不高。考虑端元可变性的光谱混合分析是提高高光谱解混精 度的关键技术。

目前光谱可变的光谱混合分析中应用最广泛的是Roberts等人提出的MESMA(Multiple Endmember Spectral Mixture Analysis)方法，该方法需要穷尽所有端元组合， 随着可变端元数增多，计算量将变得十分巨大。为了降低计算量，后人在MESMA 基础上提出了ISMA(Iterative Spectral Mixture Analysis)方法和HMESMA (HierarchicalMultiple Endmember Spectral Mixture Analysis)方法。ISMA算法的基 本原理是以所有地物的所有可变端元构成一个大的初始端元集，对每个像元，通过迭 代剔除丰度最小的端元生成候选端元集，再根据重构误差变化量确定最优端元集解混 相应的像元。HMESMA算法的基本原理是以所有地物的所有可变端元构成一个大 的初始端元集，对每个像元，根据每类端元的最大丰度选择单端元构成其初始单端元 集，进而根据解混误差变化量再分层确定最优端元集。对于包含m种地物的场景，这 两种算法都是假设所有像元内可能包含的地物种类数为1到m，通过迭代确定像元内 实际包含的地物数及最优端元组合。实际场景中，大多数像元内只含一种或少数几种 地物，因此，这两种算法仍然包含一些不必要的丰度求解计算，当地物种类比较多时， 计算量仍然比较大。

发明内容

针对现有技术中的不足，本发明的目的是提供一种端元可变的高光谱图像光谱混合分析方法，通过预确定像元混合特性及混合像元地物组成，再结合解混精度比较高 的MESMA方法或IMSA方法，避免一些不必要的丰度求解运算的同时提高解混效 率，有效克服光谱可变性对光谱混合分析结果的影响。

下面介绍与本发明相关的一些基本变量。

考虑端元可变性，每类地物含多个端元，构成一个端元子集。设m种地物，光谱 波段数为L，n_i为第i类地物类内端元光谱数，总的端元数第i类端元子集 为丰度子集为其中e_ij表示第i类地物内 的第j条光谱，a_ij为其相应丰度。

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

一种端元可变的高光谱图像光谱混合分析方法，包括以下4个步骤：

步骤1、计算光谱形状稳定性指数，选择光谱形状稳定的波段子集；所述的选 择光谱形状稳定的波段子集包括：

1-1.计算每个端元的微分光谱；

1-2.利用所述微分光谱计算各波段光谱形状稳定性指数；

1-3.求所述光谱形状稳定性指数直方图，确定一个阈值；

1-4.选择光谱形状稳定性指数小于阈值的波段，构成波段子集；

所述的步骤1-1的端元的微分光谱计算公式为

其中为端元e_ij的微分光谱，e_ij表示第i类地物的第j条光谱，i＝1,2,…m，j＝1,2,…n_i，m为地物类别数，n_i为第i类地物类内端元光谱数，L为光谱波段数；

进一步地，步骤1-2的光谱形状稳定性指数计算公式为

其中E_mean,z,i和E_mean,j,i为分别为第z类和第j类端元的第i个波段微分光谱的均值，σ_z,i和 σ_j,i为第z类和第j类端元第i波段微分光谱值的标准差；

所述光谱形状稳定是指同类地物的可变端元光谱形状相似，且不同类地物的端元光谱形状差异明显；

利用微分光谱计算光谱形状稳定性指数，目的是为了确定同种地物可变端元光谱形状相似，且不同种地物的光谱间形状差异大的波段，以保证后续类端元光谱的均值 与同类可变端元光谱的相似性；

步骤2、计算类端元特征光谱，所述类端元特征光谱是指某类地物的所有可变 端元光谱的平均光谱中波段子集对应的部分光谱，其计算过程为：分别计算每类端元 集的平均光谱曲线，从各平均光谱曲线中选波段子集对应的波段，作为各类端元特征 光谱；

步骤3、根据所选的波段子集，利用各类端元特征光谱对像元进行初步光谱混 合分析，快速判定纯像元和混合像元及各像元内所含地物；

步骤4、对每个混合像元，基于像元内所含地物对像元进行精细光谱混合分析。

进一步地，步骤3所述的根据所选的波段子集，利用类端元特征光谱对像元进行初步光谱混合分析，快速判定纯像元和混合像元及各像元内所含地物包括：

3-1.根据选择的波段子集，用类端元特征光谱对每个像元进行非负约束丰度求解；

3-2.对每个像元x，其丰度值大于0的p_x种地物为该像元内所含地物，如果一个 像元内只含一种地物，则该像元为纯像元，置其纯像元标志为1，否则该像元为 混合像元，置其纯像元标志为0。

优选的，所述非负约束丰度求解方法可以采用非负约束的最小二乘法、非负约束的核最小二乘法或渐近波段处理的非负最小二乘法。

进一步地，步骤4所述的对每个混合像元，基于像元内所含地物对像元进行精细光谱混合分析采用枚举法，所述的枚举法包括如下步骤：

(1)从混合像元内所含各地物的可变端元束中分别选一个端元，构成该像元所 有的端元子集；

(2)用各端元子集分别对像元进行全约束光谱混合分析，并计算重构误差；

(3)选择重构误差最小的端元子集，其相应丰度为像元的光谱混合分析结果。

优选的，所述非负约束丰度求解方法可以采用非负约束的最小二乘法、非负约束的核最小二乘法或渐近波段处理的非负最小二乘法。

进一步地，步骤4所述的对每个混合像元，基于像元内所含地物对像元进行精细光谱混合分析采用迭代光谱混合分析法，所述的迭代光谱混合分析法包括如下步骤：

(1)混合像元x内所含p_x类地物的端元构成扩展端元集

(2)用所述扩展端元集对像元x进行全约束丰度求解；

(3)从像元x内所含每类地物的可变端元集分别选取丰度最大的端元，构成像元的初始端元集初始化迭代次数t＝1；

(4)用对像元x全约束丰度求解，得到丰度并计算重构误差r_t；

(5)从中去除丰度最小的端元，得到端元子集更新t＝t+1；

(6)重复步骤(4)和(5)，直到端元数为1；

(7)根据像元x的重构误差变化量，确定最优端元子集，其相应丰度为像元的光 谱混合分析结果。

优选的，所述全约束丰度求解方法可以采用全约束最小二乘法、全约束核最小二乘法或渐近波段处理的全约束最小二乘法。

进一步地，所述的重构误差变化量计算公式为：

Δr_t＝|r_t+1-r_t|/r_t

本发明还提供了一种端元可变高光谱图像光谱混合分析装置，所述装置包括波段选择模块、类端元特征光谱计算模块、像元内地物类别预确定模块和混合像元精细光 谱混合分析模块。

所述波段选择模块，用于选择异类端元形状差异明显而同类端元形状相似性高的稳定波段子集；

所述类端元特征光谱计算模块，用于计算每类端元的平均光谱曲线，并从中选择波段子集对应的波段，得到各类端元特征光谱；

所述像元内地物类别预确定模块，用各类平均端元光谱对像元进行初步光谱混合分析，判定混合像元内地物组成，同时筛选纯像元；

所述混合像元精细光谱混合分析模块，用于精确求取混合像元内所含地物的丰度；

进一步地，所述高光谱图像光谱混合分析装置还包括应用模块，所述应用模块使用光谱混合分析结果进行高光谱图像亚像元定位，实现混合像元内地物的具体空间分 布。

进一步地，所述的波段选择模块具体步骤包括：

(1)计算每个端元光谱的微分光谱

其中为端元e_ij的微分光谱，e_ij表示第i类地物的第j条光谱，i＝1,2,…m， j＝1,2,…n_i，m为地物类别数，n_i为第i类地物类内端元光谱数，L为光谱波段数；

(2)利用所述微分光谱计算光谱形状稳定性指数

(3)根据光谱形状稳定性指数直方图，确定一个阈值；

(4)选择光谱形状稳定性指数小于阈值的波段，构成波段子集；

进一步地，所述的混合像元精细光谱混合分析模块选用迭代光谱混合分析具体步骤 包括：

(1)混合像元x内所含p_x类地物的端元构成扩展端元集

(2)用所述扩展端元集对像元x进行全约束丰度求解；

(3)从像元x内所含每类地物的可变端元集分别选取丰度最大的端元，构成像 元的初始端元集初始化迭代次数t＝1；

(4)用对像元x全约束丰度求解，得到并计算重构误差r_t；

(5)从中去除丰度最小的端元，得到端元子集更新t＝t+1；

(6)重复步骤(4)和(5)，直到中端元数为1。

(7)

(8)计算重构误差变化量Δr_t＝|r_t+1-r_t|/r_t，t＝1,2,…p_x-1，根据重构误 差变化量，确定最优端元子集，其相应丰度为像元的光谱混合分析结果。

本发明的有益效果在于：

针对现有技术中对所有像元用同样的方法进行光谱混合分析导致运算量大的问题，本发明提出像元混合特性及像元内地物类别预确定，可以降低光谱混合分析总体 时间复杂度同时提高其精度；

本发明提出利用微分光谱计算光谱形状稳定性指数选择波段后的平均端元光谱进行光谱混合分析可以提高像元内地物类别确定的准确性；

本发明提出结合枚举法或迭代光谱混合分析法进行精细光谱混合分析可以解决像元内地物类别预确定可能引起的误差问题，进一步提高光谱混合分析精度。

附图说明

图1为本发明高光谱图像光谱混合分析方法实施例的基本步骤流程图。

图2为本发明高光谱图像光谱混合分析装置的结构示意图。

图3为真实地物类别标签图。

图4为6种地物的真实丰度分布图。

图5为USGS光谱库中6类地物的可变端元光谱曲线。

图6为仿真的高光谱图像

图7为根据步骤一得到的所选波段光谱曲线。

图8为经过光谱混合分析得到的丰度分布图。

图9为三种算法求得的丰度用于亚像元定位的结果。

图10是EE_FCLS、ISMA、MGP_ESMA和MGP_ISMA四种方法的丰度均方根 误差及亚像元定位精度。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面结合具体的实施例详细说明本发明。下面描述了具体实施例以简化本发明。但是需要认识到，本发明不局限于 所说明的实施例，并且在不脱离基本原理的前提下，本发明的各种修改是可能的，这 些等价形式同样落于本申请所附权利要求书所限定的范围。

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提 下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，为所发明的高光谱图像光谱混合分析方法实施例的基本步骤流程图。如图2所示，为所发明的高光谱图像光谱混合分析装置的结构示意图。

输入：图像的N个像元x_k，k＝1,2,…,N，第i类端元子集其中m为地物的种类数，n_i为第i类地物类内端 元光谱数，总的端元数用e_ij表示第i类地物内的第j条光谱，j＝1,2,…,n_i， N，n_i，m都为正整数；

步骤一：基于光谱形状稳定性指数选择波段子集。

(1)计算每个端元光谱的微分光谱

其中为端元e_ij的微分光谱，e_ij表示第i类地物的第j条光谱，i＝1,2,…m， j＝1,2,…n_i，m为地物类别数，n_i为第i类地物类内端元光谱数，L为光谱波段数；

(2)计算各个波段的光谱形状稳定性指数

其中E_mean,z,i和E_mean,j,i为分别为第z类和第j类端元的第i个波段微分光谱的均值，σ_z,i和 σ_j,i为第z类和第j类端元第i波段微分光谱值的标准差。

(3)求光谱形状稳定性指数直方图hist(SSI)，根据直方图设定SSI的阈值T。

(4)求波段子集B＝∪{i}，其中i∈{1,2,…L-1}，并且SSI_i<T。

步骤二：计算类端元特征光谱

求每类端元的平均光谱根据波段子集B从中 取部分数据构成端元特征光谱

步骤三：用类端元特征光谱对像元进行光谱混合分析，确定像元内所含地物

对每个像元x_k，k＝1,2,…N，N为像元个数；

各类端元特征光谱按列向量扩展构成特征光谱矩阵对像元x_k的子集x_k(B)非负约束求丰度a_k＝[a_k1,a_k2,…,a_km]；

求x_k内含地物类别索引矩阵I_k：如果a_kj＞0，则I_k＝[I_k,j]；

如果I_k中只有一个元素，则Flag_k＝1，a_k＝0,否则Flag_k＝0；

非负约束丰度求解选用非负约束的最小二乘法，可用matlab函数lsqnonneg 实现。

步骤四：对每个混合像元，基于像元内混合地物对像元进行精细光谱混合分析。

(1)对于Flag_k＝0的混合像元x_k，计算索引矩阵I_k元素个数p_k，从所有地物的 可变端元集中选p_k种地物可变端元集并合并为扩展端元集

(2)用所述扩展端元集对像元x_k进行全约束丰度求解，

(3)从像元x_k内所含每类地物的可变端元集分别选取丰度最大的端元，构成像 元的初始端元集即:

(4)初始化迭代次数t＝1；

(5)用对像元x_k全约束丰度求解，得到并计算重构误差

(6)令从中去除第q个端元，得到端元子集 更新t＝t+1；

(7)重复步骤(5)和(6)，直到中端元数为1；

(8)

(9)计算重构误差变化量Δr_t＝|r_t+1-r_t|/r_t，t＝1,2,…p_k-1，根据重构误差 变化量，确定最优端元子集，其相应丰度为像元的光谱混合分析结果。

全约束丰度求解选用全约束最小二乘方法。已知含k个端元的矩阵M和像元 r，记1^T＝(1,1,...1)，全约束最小二乘求丰度的步骤 如下：

(1)初始化P⁽⁰⁾＝{1,2,...,k}，并令迭代次数t＝0；

(2)计算并令S_cls＝S_ls；

(3)假如S_cls的全部分量值全部是非负的，则停止迭代；否则，t＝t+1；

(4)将集合P^(t-1)中与S_cls负分量对应的元素移到集合R^(t-1)中，得到新的索 引集合P^(t)与R^(t)，并建立了另一个新集合N^(t)＝P(^t)；

(5)从S_cls中取出下标和N^(t)中的元素相等的值，记作中删除 所有行号或列号为P^(t)中元素的行与列，得到新的矩阵

(6)计算如果λ^(t)的分量值全部为负值，则直接跳转到步骤 (10)；

(7)求记在λ^(t)中的序号为q，将R^(t)中第q个元 素移到集合P^(t)；

(8)删除中列号为P^(t)中元素的所有列，得到一个新的矩阵并 计算

(9)将集合P^(t)中与负分量对应的元素移到集合R(^t)中，跳转到步骤(5)；

(10)删除中列号为P^(t)中元素的所有列，得到一个新的矩阵

(11)计算的值，跳转到步骤(3)。

下面以仿真高光谱图像为例说明具体的实施方式。

采用印第安树高光谱图像的真实标签和从USGS光谱库中的光谱数据用于 图像仿真。截取(1:60,70:144)大小的标签图像，如图3所示，按3×3比例缩小 构造混合比例，得多图4所示的真实丰度分布图；从USGS光谱库中选择6类地 物光谱，每类地物取3条光谱数据构成端元集，如图5所示；从每类端元中随机 取一个端元构成端元矩阵，并按照混合比例仿真一个像元，重复此过程，完成所 有像元数据仿真；最后，加30db的噪声，得到如图6所示的高光谱图像。

根据实施例步骤一，得到图7所示的所选波段光谱曲线图。图8为最后得到 的结果丰度分布图。为进一步验证本发明专利的应用效果，将本专利方法和其他 方法得到的结果丰度用于亚像元定位。便于描述，将直接用所有端元集全约束求 丰度进行光谱混合分析的方法记为EE_FCLS，将本专利提出的混合地物类别预 确定结合枚举光谱混合分析的算法记为MGP_ESMA(Emulative Spectral Mixture Analysis after Mixed GroundPredetermined)，把本专利提出的混合地物类 别预确定结合迭代光谱混合分析的算法记为MGP_ISMA(Iterated Spectral Mixture Analysis after Mixed GroundPredetermined)，图9为EE_FCLS、 MGP_ESMA和MGP_ISMA三种方法求得的丰度用于亚像元定位的结果，表1 为EE_FCLS、ISMA、MGP_ESMA和MGP_ISMA四种方法的丰度均方根误差 及亚像元定位精度。

图7-9和表1(图10)的结果表明：(1)基于光谱形状稳定性指数选择的波 段可以满足保持同类端元形状相似和不同类端元光谱差异的要求；(2)混合地物 类别预确定可以提高光谱混合分析精度；(3)MGP_ESMA精度最高但耗时最长， MGP_ISMA可以同时提高精度和降低时间复杂度，因此性能最好。

Claims (9)

1.一种端元可变的高光谱图像光谱混合分析方法，其特征在于包括以下4个步骤：

步骤1、计算光谱形状稳定性指数，选择光谱形状稳定的波段子集；所述的选择光谱形状稳定的波段子集包括：

1-1.计算每个端元的微分光谱；

1-2.利用所述微分光谱计算各波段光谱形状稳定性指数；

1-3.求所述光谱形状稳定性指数直方图，确定一个阈值；

1-4.选择光谱形状稳定性指数小于阈值的波段，构成波段子集；

所述的步骤1-1的端元的微分光谱计算公式为

其中为端元e_ij的微分光谱，e_ij表示第i类地物的第j条光谱，i＝1,2,…m，j＝1,2,…n_i，m为地物类别数，n_i为第i类地物类内端元光谱数，L为光谱波段数；

所述的步骤1-2的光谱形状稳定性指数计算公式为

其中E_mean,z,i和E_mean,j,i为分别为第z类和第j类端元的第i个波段微分光谱的均值，σ_z,i和σ_j,i为第z类和第j类端元第i波段微分光谱值的标准差；

所述光谱形状稳定是指同类地物的可变端元光谱形状相似，且不同类地物的端元光谱形状差异明显；

步骤2、计算类端元特征光谱，所述类端元特征光谱是指某类地物的所有可变端元光谱的平均光谱中波段子集对应的部分光谱，其计算过程为：分别计算每类端元集的平均光谱曲线，从各平均光谱曲线中选波段子集对应的波段，作为各类端元特征光谱；

步骤3、根据所选的波段子集，利用各类端元特征光谱对像元进行初步光谱混合分析，快速判定纯像元和混合像元及各像元内所含地物；

步骤4、对每个混合像元，基于像元内所含地物对像元进行精细光谱混合分析。

2.如权利要求1所述的端元可变高光谱图像光谱混合分析方法，其特征在于，步骤3所述的根据所选的波段子集，利用各类端元特征光谱对像元进行初步光谱混合分析，快速判定纯像元和混合像元及各像元内混合地物包括：

3-1.根据选择的波段子集，用各类端元特征光谱对每个像元进行非负约束丰度求解；

3-2.对每个像元x，其丰度值大于0的p_x种地物为该像元内所含地物，如果一个像元内只含一种地物，则该像元为纯像元，置其纯像元标志为1，否则该像元为混合像元，置其纯像元标志为0。

优选的，所述非负约束丰度求解方法可以采用非负约束的最小二乘法、非负约束的核最小二乘法或渐近波段处理的非负最小二乘法。

3.如权利要求1所述的端元可变高光谱图像光谱混合分析方法，其特征在于，步骤4所述的对每个混合像元，基于像元内所含地物对像元进行精细光谱混合分析采用枚举法，所述的枚举法包括如下步骤：

(1)从混合像元内所含各地物的可变端元束中分别选一个端元，构成该像元所有的端元子集；

(2)用各端元子集分别对像元进行全约束光谱混合分析，并计算重构误差；

(3)选择重构误差最小的端元子集，其相应丰度为像元的光谱混合分析结果。

优选的，所述非负约束丰度求解方法可以采用非负约束的最小二乘法、非负约束的核最小二乘法或渐近波段处理的非负最小二乘法。

4.如权利要求1所述的端元可变高光谱图像光谱混合分析方法，其特征在于，步骤4所述的对每个混合像元，基于像元内所含地物对像元进行精细光谱混合分析采用迭代光谱混合分析法，所述的迭代光谱混合分析法包括如下步骤：

(1)混合像元x内所含p_x类地物的端元构成扩展端元集

(2)用所述扩展端元集对像元x进行全约束丰度求解；

(3)从像元x内所含每类地物的可变端元集分别选取丰度最大的端元，构成像元的初始端元集初始化迭代次数t＝1；

(4)用对像元x全约束丰度求解，得到丰度并计算重构误差r_t；

(5)从中去除丰度最小的端元，得到端元子集更新t＝t+1；

(6)重复步骤(4)和(5)，直到端元数为1；

(7)根据像元x的重构误差变化量，确定最优端元子集，其相应丰度为像元的光谱混合分析结果。

优选的，所述全约束丰度求解方法可以采用全约束最小二乘法、全约束核最小二乘法或渐近波段处理的全约束最小二乘法。

5.如权利要求7所述的端元可变高光谱图像光谱混合分析方法，其特征在于，所述的重构误差变化量计算公式为：

Δr_t＝|r_t+1-r_t|/r_t

6.一种端元可变高光谱图像光谱混合分析装置，其特征在于，包括

波段选择模块，用于选择异类端元形状差异明显而同类端元形状相似性高的稳定波段子集；

类端元特征光谱计算模块，用于计算每类端元的平均光谱曲线，并从中选择波段子集对应的波段，得到各类端元特征光谱；

像元内地物类别预确定模块，用各类端元特征光谱对像元进行初步光谱混合分析，判定混合像元内地物组成，同时筛选纯像元；

混合像元精细光谱混合分析模块，用于精确求取混合像元内所含地物的丰度；

亚像元定位模块，用于根据光谱混合分析结果，实现混合像元内地物的具体空间分布。

7.根据权利要求6所述的端元可变高光谱图像光谱混合分析装置，其特征在于，所述装置还包括应用模块，所述应用模块使用光谱混合分析结果来进行高光谱图像亚像元定位。

8.根据权利要求6所述的端元可变高光谱图像光谱混合分析装置，其特征在于，

所述的波段选择模块具体步骤包括：

(1)计算每个端元光谱的微分光谱

(2)利用所述微分光谱计算光谱形状稳定性指数

(3)根据光谱形状稳定性指数直方图，确定一个阈值；

(4)选择光谱形状稳定性指数小于阈值的波段，构成波段子集。

9.根据权利要求6所述的端元可变高光谱图像光谱混合分析装置，其特征在于，所述的混合像元精细光谱混合分析模块选用迭代光谱混合分析的具体步骤包括：

(1)混合像元x内所含p_x类地物的端元构成扩展端元集

(2)用所述扩展端元集对像元x进行全约束丰度求解；

(3)从像元内所含每类地物的可变端元集分别选取丰度最大的端元，构成像元的初始端元集初始化迭代次数t＝1；

(4)用对像元x全约束丰度求解，得到并计算重构误差r_t；

(5)从中去除丰度最小的端元，得到端元子集更新t＝t+1；

(6)重复步骤(4)和(5)，直到中端元数为1。

(7)

(8)计算重构误差变化量Δr_t＝|r_t+1-r_t|/r_t，t＝1,2,…p_x-1，根据重构误差变化量，确定最优端元子集，其相应丰度为像元的光谱混合分析结果。