CN105976404B - 高光谱遥感图像的线性全约束丰度估计方法 - Google Patents

Abstract

本发明公开了一种高光谱遥感图像的线性全约束丰度估计方法，构建原光谱解混问题的对偶问题，通过全约束最小二乘法对丰度进行初始化，结合偏离度ρ因素利用原始对偶内点法优化丰度结果，具体包括以下步骤：S1:建立图像的线性光谱混合模型：S2:对线性光谱混合模型施加约束条件；S3:采用原始对偶内点法和全约束最小二乘丰度估计算法对线性光谱混合模型进行混合像元的丰度估计。该发明有效地克服了全约束最小二乘法精度不高的问题，同时也解决了原始对偶内点法速度慢的问题，可以作为一种解决高光谱遥感图像全约束丰度估计的有效手段。在基于高光谱遥感图像的高精度的解混以及地面目标的检测和识别方面具有重要的应用价值。

Description

高光谱遥感图像的线性全约束丰度估计方法

技术领域

本发明涉及遥感图像处理技术领域，尤其涉及高光谱遥感图像的线性全约束丰度估计方法。

背景技术

高光谱遥感技术能够实现地物空间信息、辐射信息、光谱信息的同步获取，具有波段信息丰富、分辨率高、图谱和一三大特点，因而包含了更多精细的地物信息。高光谱遥感数据最终形成的是一个高光谱数据立方体，这个数据立方体就是一个包含了地物信息的数据宝库，只有对这个数据宝库进行发掘、分析，才能得到有用的信息。

由于受到遥感器分辨率的限制和地物信息复杂性的影响，高光谱遥感数据立方体中包含了大量的混合像元。丰度估计就是对混合像元中端元所占的比例进行估计，所使用的光谱混合模型主要有线性混合模型和非线性混合模型。线性混合模型假设混合像元由几种纯地物按照一定的比例线性混合而成，具有明确的物理意义且容易建模，得到了广泛的应用。

根据实际的物理意义，混合像元的分解需要满足两项约束条件，一是非负约束(Abundance Non-negative Constraint ANC)，即端元丰度不能出现负值；一是和为1约束(Abundance Sum-to-one Constraint,ASC)，即各端元丰度和为1。满足ANC和ASC约束的丰度估计算法即为全约束丰度估计算法。现有技术中采用单一的算法如全约束最小二乘法进行。

发明内容

根据现有技术存在的问题，本发明公开了一种高光谱遥感图像的线性全约束丰度估计方法，构建原光谱解混问题的对偶问题，通过全约束最小二乘法对丰度进行初始化，结合偏离度ρ因素利用原始对偶内点法优化丰度结果，具体包括以下步骤：

S1:建立图像的线性光谱混合模型：

S2:对线性光谱混合模型施加约束条件；

S3:采用原始对偶内点法和全约束最小二乘丰度估计算法对线性光谱混合模型进行混合像元的丰度估计：在该过程中融合偏离度ρ作为丰度估计的一个因素，用偏离度ρ表征混合像元中不能被解析的部分，将端元数量、地形、遥感器、端元纯净性、噪声因素所引起的影响使用解析偏差表示；从而获得高光谱图像中每个混合像元的最佳偏离度ρ_z，得到最佳的丰度估计结果实现原图像的重构。

所述线性光谱混合模型为：

其中：r是混合像元中任意的l维光谱向量，p是端元个数，M＝[m₁,m₂,....m_p]是大小为l*p的端元光谱矩阵，α是p维矢量系数，其各分量元素是混合像元r中端元m_i的丰度,e为l维矢量，表示模型误差项。

将约束条件施加于(1)式上，可以得到ASC约束为：ANC约束为：α_i≥0(i＝1,2,..p)，在最小二乘法算法的基础上加入ASC和ANC约束，则全约束最小二乘法丰度估计算法为：

min_|α∈Ω{(r-Mα)^T(r-Mα)}

S3中将采用原始对偶内点法和全约束最小二乘法相结合时具体为：设公式(2)为原始问题，根据原始对偶理论，引入参数λ，z＝(z₁,z₂,...z_p)^T，则公式(2)的对偶问题如下：

根据KKT条件，公式(2)(3)的最优解等价于求解如下参数方程：

其中X＝diag(α₁,α₂,...,α_p),Z＝diag(z₁,z₂,...z_p)，e＝(1₁,1₁,...1_p)^T，u为障碍因子，u→0。公式(4)中第一个式子为原可行性，第二个式子为对偶可行性，第三个式子为互补松弛性；

根据原始对偶理论，对偶间隙为J(α)-D(α)＝α^Tz根据公式(4)中的第三个式子，对偶间隙简化表示为J(α)-D(α)＝pu，p为端元个数，当u→0时，J(α)-D(α)→0，设对偶问题解为W＝(α,λ,υ)，则α、W则分别收敛到原问题(2)和对偶问题(3)的最优解，则W(α)即原始对偶内点法全约束丰度估计的结果；

设第k次迭代的丰度估计的结果为W^k＝(α^k,λ^k,z^k)，丰度估计的步长为ΔW＝(Δα,Δλ,Δz)，ΔX＝diag(Δα₁,Δα₂,...,Δα_p)ΔZ＝diag(Δz₁,Δz₂,...Δz_p),Δz＝(Δz₁,Δz₂,...Δz_p)^TΔα＝(Δα₁,Δα₂,...,Δα_p)^T，采用牛顿法对公式(4)求解，可得：

根据公式(5)可以计算得到ΔW＝(Δα,Δλ,Δz)，记αk＝1-1^Tα^k，dk＝(M^TM)α^k-1λ^k-z^k-M^Tr，zz＝||dk||₂。

在原始对偶内点法算法中加入全约束最小二乘丰度估计算法以及偏离度ρ具体采用如下方式：

Step1运行FCLS算法，以运行结果初始化α⁰。

Step2初始化ρ，需要根据具体情况设置，例ρ＝0.05。

Step3初始化λ⁰，z⁰，由于α⁰的准确性很高，因此λ⁰，z⁰设置为较小的数值，例如1e-8，得到W⁰＝(α⁰,λ⁰,z⁰)。令p为端元个数，同时给定较小的正数0<δ<1；给定终止误差ε>0，置N＝0，最大迭代次数N_s；

Step4当((α^k)^Tz≤ε)&&(αk≤ε)&&(zz≤ε)，则停止丰度的查找，α^k为丰度估计最优值；否则转到Step5。

Step5令μ^k+1＝μ^k(1-δ)，根据公式(8)计算ΔW^k＝(Δα,Δλ,Δz)

Step6W^k+1＝W^k+σΔW^k，N＝N+1。

Step7判断光谱重构误差是否递减，递减则转到Step 8，否则α^k为丰度估计最优值，算法结束。

Step8如果N>＝N_s，算法结束，否则转到Step3。如果将Step8中N_s＝1，即为FCLS+PDIP(N＝1,ρ)算法。经过数据测试发现，该算法分解精度的提高主要集中在迭代的首次。因此，可以通过控制N限制迭代次数，牺牲部分分解精度来较大的降低算法的运行时间。

在混合像元的丰度估计中，将由于端元数量、地形、遥感器、端元纯净性、噪声等一切因素所引起的影响都使用解析偏差统一表示，ρ即表示这种偏离度的大小。对每一个混合像元而言，都有一个最佳的ρ_z，使得该混合像元的分解精度达到最佳。

由于采用了上述技术方案，本发明提供的高光谱遥感图像的线性全约束丰度估计方法，通过构造原解混问题的对偶问题，结合传统的全约束最小二乘丰度估计FCLS方法,将高光谱图像全约束丰度估计问题转化为二次规划的优化问题。进而，该发明引入了丰度的和为一偏离度，对于系统不能解析的噪声等因素进行表示，并将对偶间隙和解混误差同时作为优化目标，提高了算法的解混精度。该发明从全约束线性高光谱混合模型出发，结合最小二乘方法和原始对偶内点法，有效地克服了全约束最小二乘法精度不高的问题，同时也解决了原始对偶内点法速度慢的问题，可以作为一种解决高光谱遥感图像全约束丰度估计的有效手段。在基于高光谱遥感图像的高精度的解混以及地面目标的检测和识别方面具有重要的应用价值。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明中Cuprite数据图的示意图；

图2为本发明中从Cuprite数据中提取的5种矿物质的端元光谱图；

图3为本发明中模拟高光谱数据合成区域图；

图4是模拟高光谱数据彩色展示图；

图5是FCLS、PDIP和优化FCLS+PDIP(N＝1,ρ)算法的SRMSE和ρ值关系示意图；

图6是真实高光谱图像Hydice数据第120波段示意图。

具体实施方式

为使本发明的技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：

如图1所示的一种高光谱遥感图像的线性全约束丰度估计方法，构建原光谱解混问题的对偶问题，通过全约束最小二乘法对丰度进行初始化，结合偏离度ρ因素利用原始对偶内点法优化丰度结果，具体包括以下步骤：

S1:建立图像的线性光谱混合模型：

S2:对线性光谱混合模型施加约束条件；

S3:采用原始对偶内点法和全约束最小二乘丰度估计算法对线性光谱混合模型进行混合像元的丰度估计：在该过程中融合偏离度ρ作为丰度估计的一个因素，用偏离度ρ表征混合像元中不能被解析的部分，将端元数量、地形、遥感器、端元纯净性、噪声因素所引起的影响使用解析偏差表示；从而获得高光谱图像中每个混合像元的最佳偏离度ρ_z，得到最佳的丰度估计结果实现原图像的重构。

实施例：分别用仿真数据和实际遥感图像数据为例，说明本发明的具体实施方式。我们用重构误差和运行时间来衡量全约束解混算法的优劣，运行环境为Windows7系统下的Matlab7.0平台。

重构误差采用光谱均方根误差(Spectral Root Mean Square Error,SRMSE)作为丰度评价标准，如公式(6)所示。SRMSE值越小表明丰度估计的结果越准确，在混合像元分解中，可以代表混合像元分解的精度，是一个绝对评价标准。具体计算公式如下：

其中，p是波段个数，N是像素的个数，r_ij是第i个混合像元在第j波段的光谱亮度值，M为端元矩阵，α_i是第i个混合像元的丰度估计值。

模拟高光谱数据Cuprite数据是一个公开的高光谱数据，因此以Cuprite图像数据为基础，提取5种矿物质A、B、C、D、E，分别是Kaolinite,Calcite,Alunite，Buddingtonite，Muscovite,如图1所示，它们的辐射亮度值曲线如图2所示。

以5种矿物质为端元合成波段189个波段的32*32大小的模拟图像。将图像划分为A，B，C，D,e,h,g,h,o共9个区域，如图3所示。将Calcite作为基础物质，在9个区域内都添加。以在A区域内每个像元添加80％的Alunite和20％的Calcite；在B区域内添加80％的Buddingtonite和20％的Calcite；在C区域内添加80％的Muscovite和20％的Calcite；在D区域内添加80％的Kaolinite和20％的Calcite；在e区域内添加40％的Alunite，40％的Buddingtonite以及20％的Calcite；在f区域内添加40％的Alunite，40％的Muscovite以及20％的Calcite；在g区域内添加40％的Muscovite，40％的Kaolinite以及20％的Calcite；在h区域内添加40％的Kaolinite，40％的Buddingtonite以及20％的Calcite；在o区域内添加20％的Alunite，20％的Buddingtonite,20％的Calcite，20％的Kaolinite以及20％的Muscovite；同时为每个像元中加入30:1的信噪比，完成模拟混合像元数据。从189个波段中提取3个波段合成彩色图4所示。

仿真实验中的测试数据共有4组，本别依据FCLS、PDIP、FCLS+PDIP(ρ)、FCLS+PDIP(ρ,N＝1)产生。PDIP、FCLS+PDIP(ρ)、FCLS+PDIP(ρ,N＝1)中的参数ε设置为1e-8；0<δ<1，取δ＝0.1；更新步长σ＝0.9。

实验1 测试重构误差与偏离度ρ之间的关系。在这个实验中，我们改变偏离度ρ的值，分别测试各算法进行全约束解混时的重构误差。表1、2分别给出

了FCLS、PDIP的数据结果。表1FCLS算法在不同ρ值下的丰度估计结果

表1 FCLS算法在不同ρ值下的丰度估计结果

表2 PDIP算法在不同ρ值下丰度估计结果

从表1、2可以看出，适当的增加ρ值可以进一步的降低SRMSE值，也就是提高了分解精度，说明ρ预设值在向最佳ρ_z靠近。如果继续增大ρ值，SRMSE值反而会升高，也就是降低了分解精度，说明ρ的预设值在偏离最佳ρ_z。因此，要想获得最佳的丰度估计结果，就要找到最佳ρ_z，而且是找到每个混合像元的最佳ρ_z，才能实现最佳的丰度估计。

从表1、2中还可以看出，PDIP算法可以事先精确的控制每个混合像元的ρ值，而FCLS算法则做不到。FCLS算法得到的ρ值范围在预设值上下波动，从表1中可以看出，ρ均值小于预设值。从分解结果上看，FCLS算法的SRMSE值略小于PDIP算法，也就是分解精度较高，说明FCLS算法比PDIP算法在噪声处理上较强，对ρ值有一定的自适应性，有自动寻找最佳ρ_z的趋势。

PDIP算法可以事先精确的控制ρ值，这是该算法的特点，但是在丰度估计中并不是优点，因为每个混合像元都有自己最佳ρ_z，因此事先精确控制ρ_z并不是一个好的选择。但是，如果以FCLS算法分解结果的ρ值作为标准，事先设置在PDIP算法上，经过数据测试可以看到，PDIP算法得到的丰度值、SRMSE值与FCLS算法一致。这也说明了最佳的丰度估计在于对噪声的处理，在于对最佳ρ_z的寻找。

由于事前并不知道最佳ρ_z，因此ρ的设置上具有盲目性，如果设置的不好，不仅起不到提高分解精度的目的，反而会使分解经过更差。例如根据数据的模拟过程，信噪比为30:1，因此如果设置ρ＝1/30，不管是FCLS算法还是PDIP算法，SRMSE值反而高于ρ＝0时的分解结果。如果用NSCLS算法做对比，NSCLS算法是实现了ASC约束，但是在实现ANC约束时，采用标准化的方式进行，打破了目标函数最小化的约束，因此是一个分解效果比较差的全约束丰度估计算法。NSCLS算法分解结果的SRMSE值为126.09609195，FCLS、PDIP算法的SRMSE在ρ＝1/30时已经比较接近该值了，说明FCLS、PDIP算法在ρ＝1/30时已经退化成了一个比较差的算法了。因此，在FCLS、PDIP算法上如果添加ρ值不当，会严重影响算法的分解精度，也说明了在FCLS、PDIP算法上直接添加ρ值，可行性不高。

然后，通过结合FCLS与PDIP，加入解混误差标准，使得该算法对ρ值有了较强的自适应性，提高了分解精度，而且大幅度降低了在ρ值设置不合理的情况下分解精度的误差，避免了ρ值设置不合理的危险性。表3是优化FCLS+PDIP(N＝1,ρ)算法在不同ρ值下的丰度估计结果。图5是优化FCLS+PDIP(N＝1,ρ)算法与FCLS、PDIP算法分解结果的比较。

表3 优化FCLS+PDIP(N＝1,ρ)算法在不同ρ值下丰度估计结果

从图5以及表1、2、3可以看到，优化FCLS+PDIP(N＝1,ρ)算法在分解精度上有了进一步的提高，而且在ρ值设置较大时，并不会降低分解的精度，相对于原始FCLS、PDIP算法，分解精度仍有提高。因此，本发明算法对ρ值有了一定的自适应性，强化了ρ值的自动寻找最佳ρ_z的能力。

实验2 测试各算法的运行时间，每种算法的每组参数设置都在测试数据上运行20次后取平均值。PDIP算法花费时间较长，FCLS和PDIP算法想结合之后，虽然相对于PDIP能够降低算法的部分运行时间，但是仍然较长。由于该算法分解精度的提高主要集中在首次，因此可以通过牺牲部分分解精度，控制迭代次数(N＝1)来有效的降低算法的运行时间。以ρ＝0.003为例，迭代次数、SRMSE值、时间如表4所示。

表4 FCLS+PDIP(ρ)算法在ρ＝0.003下的丰度估计

从表4可以看出，优化FCLS+PDIP(ρ)在分解精度上主要集中在首次FCLS+PDIP(N＝1,ρ)，牺牲之后的迭代次数可以有效的降低分解时间，使得分解时间靠近FCLS算法。针对不同的ρ值，3种算法的运行时间结果如表5所示。

表5 各算法在不同ρ值下的运行时间

结合之前的分析，优化FCLS+PDIP(N＝1,ρ)算法在分解精度上有了较大的提高，而且为了兼顾分解时间，控制N＝1，使得优化FCLS+PDIP(N＝1,ρ)算法成为了一个比较高效的算法。在模拟高光谱数据集上，令ρ＝0.01，比较本发明算法和原有算法的分解精度，以及本发明算法的分解优势，如表6所示。

表6 优化FCLS+PDIP(N＝1,ρ＝0.01)算法相对于其他算法在分解精度上的提高比例。

2、实际遥感高光谱数据

真实高光谱数据采用Hydice数据进行测试，共169波段，光谱空间分辨率为1.5米。图像大小为64*64，图6显示的是第120波段的图像。Hydice数据主要包含interferer,grass,tree,road及其他地物。本发明分别从中提取3、6、9、12、18、24、36个端元进行全约束丰度估计测试。

实验1 测试各算法的解混精度。在真实高光谱数据上，由于地物环境的复杂性，噪声通常较大，认为5％的噪声是一个合理的预设值，因此设置ρ＝0.05。验证不同端元下的丰度估计结果，其重构误差和相他算法的精度提高比例如表7、8所示。

表7 全约束丰度估计算法的SRMSE值(ρ＝0.05)

表8 优化FCLS+PDIP(N＝1，ρ＝0.05)算法相对于其他算法在分解精度上的提高比例

在表7中，从纵向数据比较来看，每种算法都随着端元个数的增多，SRMSE值逐渐减小，也就是说分解精度都在提高。从横向数据比较来看，FCLS算法SRMSE值相对于PDIP算法只有微弱的降低，而本发明算法较大幅度的降低了SRMSE值。结合表8可以看到，本发明算法相对于FCLS、PDIP算法在分解精度上，提高比例基本都在10％以上，提高的幅度较大。因此，可以看到本发明算法在端元确定的情况下，不用通过增加端元个数，也能够较大幅度的提高分解精度。

FCLS+PDIP(N＝1,ρ)算法虽然能够一定程度的自动寻找最佳ρ_z，但仍然不能找到最佳ρ_z，可以通过设置不同的ρ值，找到一个比较精确的丰度估计结果。表9是在不同ρ的预设值下的丰度估计结果。

表9 优化FCLS+PDIP(N＝1,ρ)算法丰度估计结果SRMSE值

可以看到在ρ的预设值在0.1时，丰度估计的精度比ρ＝0.05时更好。在端元个数为3、6时，ρ＝0.15则SRMSE值进一步减小，因此对于Hydice数据，在不同的端元个数下，都有其对应的最佳ρ_z，要针对每一种情况进行特定的分析，才能得到较好的丰度估计结果。如果能做到对每个混合像元都做到个性化分析，则能够得到更为精确的丰度估计结果。

实验2 测试各算法的运行时间，结果如表10所示。

表10 全约束丰度估计算法的运行时间

从表10中可以看出，优化FCLS+PDIP(N＝1,ρ)算法随着端元数目的增多，相对于N不限时节省的时间也越来越多，端元数目从9个开始，FCLS+PDIP(N＝1,ρ)算法的执行时间逐步靠近FCLS算法。因此，本发明FCLS+PDIP(N＝1,ρ)算法是一个兼顾了分解精度与分解时间的比较高效的一个算法。

综上所述，可以得出两点结论：一是本发明提出了混合像元分解精度和ρ值之间深刻关系，只有找到高光谱图像中每个混合像元的最佳ρ_z，才能得到最佳的丰度估计结果；二是优化FCLS+PDIP(N＝1,ρ)算法虽然不能完全实现最佳ρ_z的自动查找功能，但是避免了ρ值设置不合理的危险性，而且在一定程度上实现ρ值的自动查找功能，较大程度的提高了混合像元的分解精度，而且兼顾了算法的运行时间，算法运行也较快，因此本发明仍然是一个比较高效的算法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (5)

1.一种高光谱遥感图像的线性全约束丰度估计方法，其特征在于：构建原光谱解混问题的对偶问题，通过全约束最小二乘法对丰度进行初始化，结合偏离度ρ因素利用原始对偶内点法优化丰度结果，具体包括以下步骤：

S1:建立图像的线性光谱混合模型：

S2:对线性光谱混合模型施加约束条件；

S3:采用原始对偶内点法和全约束最小二乘丰度估计算法对线性光谱混合模型进行混合像元的丰度估计：在该过程中融合偏离度ρ作为丰度估计的一个因素，用偏离度ρ表征混合像元中不能被解析的部分，其中偏离度是对丰度的“和为1”约束的松弛，ρ将端元数量、地形、遥感器、端元纯净性、噪声因素所引起的影响使用解析偏差表示；从而获得高光谱图像中每个混合像元的最佳偏离度ρ_z，得到最佳的丰度估计结果实现原图像的重构。

2.根据权利要求1所述的高光谱遥感图像的线性全约束丰度估计方法，其特征还在于：所述线性光谱混合模型为：

其中：r是混合像元中任意的l维光谱向量，p是端元个数，M＝[m₁,m₂,....m_p]是大小为l*p的端元光谱矩阵，α是p维矢量系数，其各分量元素是混合像元r中端元m_i的丰度,e为l维矢量，表示模型误差项。

3.根据权利要求2所述的高光谱遥感图像的线性全约束丰度估计方法，其特征还在于：将约束条件施加于(1)式上，可以得到ASC约束为：ANC约束为：α_i≥0(i＝1,2,..p)，在最小二乘法算法的基础上加入ASC和ANC约束，则全约束最小二乘法丰度估计算法为：

min_|α∈Ω{(r-Mα)^T(r-Mα)}

4.根据权利要求3所述的高光谱遥感图像的线性全约束丰度估计方法，其特征还在于：S3中将采用原始对偶内点法和全约束最小二乘法相结合时具体为：设公式(2)为原始问题，根据原始对偶理论，引入参数λ，z＝(z₁,z₂,...z_p)^T，则公式(2)的对偶问题如下：

根据KKT条件，公式(2)(3)的最优解等价于求解如下参数方程：

其中X＝diag(α₁,α₂,...,α_p),Z＝diag(z₁,z₂,...z_p)，e＝(1₁,1₁,...1_p)^T，u为障碍因子，u→0；公式(4)中第一个式子为原可行性，第二个式子为对偶可行性，第三个式子为互补松弛性；

J(α)为原始问题； 根据原始对偶理论，对偶间隙为J(α)-D(α)＝α^Tz根据公式(4)中的第三个式子，对偶间隙简化表示为J(α)-D(α)＝pu，p为端元个数，当u→0时，J(α)-D(α)→0，设对偶问题解为W＝(α,λ,z)，则α、W则分别收敛到原问题(2)和对偶问题(3)的最优解，则W(α)即原始对偶内点法全约束丰度估计的结果；

设第k次迭代的丰度估计的结果为W^k＝(α^k,λ^k,z^k)，丰度估计的步长为ΔW＝(Δα,Δλ,Δz)，ΔX＝diag(Δα₁,Δα₂,...,Δα_p)ΔZ＝diag(Δz₁,Δz₂,...Δz_p),Δz＝(Δz₁,Δz₂,...Δz_p)^TΔα＝(Δα₁,Δα₂,...,Δα_p)^T，采用牛顿法对公式(4)求解，可得：

根据公式(5)可以计算得到ΔW＝(Δα,Δλ,Δz)，记αk＝1-1^Tα^k，dk＝(M^TM)α^k-1λ^k-z^k-M^Tr，zz＝||dk||₂。

5.根据权利要求4所述的高光谱遥感图像的线性全约束丰度估计方法，其特征还在于：在原始对偶内点法算法中加入全约束最小二乘丰度估计算法以及偏离度ρ具体采用如下方式：

Step1 运行FCLS算法，以运行结果初始化α⁰；

Step2 初始化ρ，根据实际情况设置ρ的数值；

Step3 初始化λ⁰，z⁰，得到W⁰＝(α⁰,λ⁰,z⁰)，令p为端元个数，设0<δ<1；给定终止误差ε>0，置N＝0，最大迭代次数N_s；

Step4 当((α^k)^Tz≤ε)&&(αk≤ε)&&(zz≤ε)，则停止丰度的查找，α^k为丰度估计最优值；否则转到Step5；

Step5 令μ^k+1＝μ^k(1-δ)，根据公式(5)计算ΔW^k＝(Δα,Δλ,Δz)；

Step6 W^k+1＝W^k+σΔW^k，N＝N+1；

Step7 判断光谱重构误差是否递减，递减则转到Step 8，否则α^k为丰度估计最优值，算法结束；

Step8 如果N>＝N_s，算法结束，否则转到Step3。