CN110428454B

CN110428454B - 一种高光谱解混方法、装置、电子设备及存储介质

Info

Publication number: CN110428454B
Application number: CN201910743253.0A
Authority: CN
Inventors: 徐翔; 李军; 邹昆; 姚琼
Original assignee: University of Electronic Science and Technology of China Zhongshan Institute
Current assignee: University of Electronic Science and Technology of China Zhongshan Institute
Priority date: 2019-08-13
Filing date: 2019-08-13
Publication date: 2022-02-25
Anticipated expiration: 2039-08-13
Also published as: CN110428454A

Abstract

本申请涉及遥感技术领域，提供一种高光谱解混方法、装置、电子设备及存储介质。其中，高光谱解混方法包括：获取待解混的高光谱图像，该高光谱图像能够表示为端元光谱信号的线性组合，线性组合系数为端元光谱信号的丰度；对第一约束方程进行优化求解，获得端元光谱矩阵以及丰度矩阵；第一约束方程包括利用NMF模型对高光谱图像进行解混时使用的第二约束方程，以及，对丰度矩阵进行稀疏约束的正则化项。该高光谱解混方法实现了对第一约束方程的稳定求解，并且求解过程具有较快的收敛速度，此外，由于第一约束方程中引入的正则化项同时体现了属性多样性和稀疏性，因此方法还可以改善高光谱解混的精度。

Description

一种高光谱解混方法、装置、电子设备及存储介质

技术领域

本发明涉及遥感技术领域，具体而言，涉及一种高光谱解混方法、装置、电子设备及存储介质。

背景技术

高光谱遥感图像通常包含了数以百计的窄波段光谱通道，能识别到地物精细的光谱信息，为定量遥感的应用带来了一种新的有效手段。然而，由于空间分辨率的限制，以及地物分布的复杂多样性，使得高光谱图像中的大量像元仍然以混合像元的形式存在，即单个像元的光谱信号是由像元内各成分物质的光谱信号(称为端元光谱信号)混合而成。混合像元的存在给高光谱图像的应用(例如，地物分类、识别等)造成了不利影响，近年来，国内外研究者提出了一些高光谱解混方法，以期分离出像元内各成分物质的端元光谱信号和对应的信号丰度。

站在统计学角度，高光谱解混问题可以看作是一个盲源分离问题，因此有研究者提出采用非负矩阵分解(Nonnegative Matrix Factorization，NMF)来开展高光谱图像进行解混。然而，现有的基于NMF模型的方法解混结果不甚理想。

发明内容

本申请实施例的目的在于提供一种高光谱解混方法、装置、电子设备及存储介质，以改善上述技术问题。

为实现上述目的，本申请提供如下技术方案：

第一方面，本申请提供一种高光谱解混方法，包括：获取待解混的高光谱图像，所述高光谱图像能够表示为端元光谱信号的线性组合，线性组合系数为所述端元光谱信号的丰度；对第一约束方程进行优化求解，获得所述高光谱图像对应的端元光谱矩阵以及丰度矩阵；其中，所述第一约束方程包括：利用非负矩阵分解NMF模型对所述高光谱图像进行解混时使用的第二约束方程，以及，对所述丰度矩阵进行行稀疏约束的正则化项。

首先，在NMF模型中使用的第二约束方程，由于只引入了非负约束，因此难以得到模型的稳定解，本申请在非负约束的基础上引入了稀疏约束，从而可以对端元光谱矩阵和丰度矩阵进行稳定求解。

其次，发明人长期研究发现，每种地物类别(对应一个端元光谱信号)的丰度图像(对应丰度矩阵的一行)往往具有不同的空间结构属性(如空间分布、几何形态等)，即满足属性多样性。进一步的，运用稀疏表示理论，当给定一个具有表征能力的字典时，每种地物类别的丰度图像可以进行稀疏表示，同时，由于不同的丰度图像之间存在属性的差异，因此在属性多样性的前提下，各丰度图像可以进行不同的稀疏表达(即不共享显著的稀疏系数)，从而可以基于此进行地物类型的判别与分离。在本申请引入第一约束方程的行稀疏约束(即正则化项)中，由于既体现了高光谱图像中空间属性的多样性，又体现了丰度图像在表征上的稀疏性，因此上述方法还可以改善高光谱解混的精度。

在第一方面的一种实现方式中，所述第一约束方程为

其中，所述第二约束方程为

所述正则化项为

Y为所述高光谱图像，M为所述端元光谱矩阵，S为所述丰度矩阵，‖.‖_F表示弗罗贝尼乌斯(Frobenius)范数，‖.‖₁表示L1范数，p为S的行数，s_j为S的第j行，λ_j为s_j对应的阈值，i⁺(.)为非负约束指示函数，定义为

在第一方面的一种实现方式中，对第一约束方程进行优化求解，包括：对由所述第一约束方程拆分而成的第三约束方程以及第四约束方程利用前后向分割(Forward-Backward Splitting，FBS)算法进行交替优化求解；其中，所述第三约束方程为

所述第四约束方程为

在对所述第三约束方程进行求解时，将M固定不变，在对所述第四约束方程进行求解时，将S固定不变。

在第一方面的一种实现方式中，对所述第三约束方程利用FBS算法进行优化求解，包括：对由所述第三约束方程拆分而成的第一约束项

利用梯度下降算法求解，以及，对由所述第三约束方程拆分而成的第二约束项

利用非负软阈值算子求解，得到所述第三约束方程的迭代公式

其中，

为迭代后的S，L_M＝||M^TM||_s′为M的谱范数，λ为S对应的阈值，[.]₊表示正象限投影函数，定义为[x]₊＝max(x,0)。

在第一方面的一种实现方式中，对所述第四约束方程利用FBS算法进行优化求解，包括：对由所述第四约束方程拆分而成的第三约束项

利用梯度下降算法求解，以及，由所述第四约束方程拆分而成的第四约束项g(M)＝i⁺(M)利用正象限投影函数进行投影，得到所述第四约束方程的迭代公式

其中，

为迭代后的M，L_S＝||SS^T||_s′为S的谱范数，[.]₊表示正象限投影函数，定义为[x]₊＝max(x,0)。

以上三种实现方式简要描述了优化求解第一约束方程的一些可能的方式。由于第一约束方程是非凸问题，因此在一种实现方式中可以将其拆分为第三约束方程以及第四约束方程两个凸问题，然后进行交替优化求，以获得第一约束方程的稳定解。

进一步的，在一种实现方式中，对于第三约束方程可以分割为第一约束项和第二约束项进行求解。其中，第一约束项是凸可微的，因此可以直接用梯度下降算法求解，而第二约束项是凸不可微的，因此可以用非负软阈值算子进行近似计算。在求解时，还可以采用FBS算法加快收敛速度，FBS算法融合了对第一约束项和第二约束项的求解过程，得到最终的更新公式

类似的，在一种实现方式中，对于第四约束方程可以分割为第三约束项和第四约束项进行求解。其中，第三约束项是凸可微的，因此可以直接用梯度下降算法求解，而第四约束项是凸不可微的，因此可以用正象限投影函数进行投影。在求解时，还可以采用FBS算法加快收敛速度，FBS算法融合了对第三约束项和第四约束项的求解过程，得到最终的更新公式

在第一方面的一种实现方式中，在优化求解的每一轮迭代过程中减小λ的取值。

在第一方面的一种实现方式中，在优化求解的每一轮迭代过程中减小λ的取值，包括：在优化求解的每一轮迭代后对残差res＝Y-MS的标准差进行估计，并根据获得的估计值减小λ的取值。

以上两种实现方式在优化求解的迭代过程中不断缩小阈值λ，其中，包括利用对残差的标准差的估值来动态地调整λ的方式。这些实现方式的意义在于：刚开始迭代时，通过设置较大的阈值针对高光谱图像中的显著特征进行考虑，随着迭代轮次的增加，逐渐设置较小的阈值以考虑高光谱图像中更细小的特征，从而达到快速调节模型参数，提高方法的收敛速度以及求解稳定性的目的。

第二方面，本申请实施例提供一种高光谱解混装置，包括：图像获取模块，用于获取待解混的高光谱图像，所述高光谱图像能够表示为端元光谱信号的线性组合，线性组合系数为所述端元光谱信号的丰度；解混模块，用于对第一约束方程进行优化求解，获得所述高光谱图像对应的端元光谱矩阵以及丰度矩阵；其中，所述第一约束方程包括：利用非负矩阵分解NMF模型对所述高光谱图像进行解混时使用的第二约束方程，以及，对所述丰度矩阵进行行稀疏约束的正则化项。

第三方面，本申请实施例提供一种电子设备，包括存储器以及处理器，所述存储器中存储有计算机程序指令，所述计算机程序指令被所述处理器读取并运行时，执行第一方面或第一方面的任意一种可能的实现方式提供的方法的步骤。

第四方面，本申请实施例提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序指令，所述计算机程序指令被处理器读取并运行时，执行第一方面或第一方面的任意一种可能的实现方式提供的方法的步骤。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对本申请实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为线性光谱混合模型的示意图；

图2为本申请实施例提供的一种高光谱解混方法的流程图；

图3为模拟数据下估计端元与地物参考端元之间的平均SAD；

图4为模拟数据下估计丰度与地物参考丰度之间的RMSE；

图5为模拟数据下不同高光谱解混方法在低混以及高混下的平均SAD；

图6为HYDICE城市影像数据下估计端元与地物参考端元之间的SAD；

图7为HYDICE城市影像数据下估计端元与参考端元光谱曲线的对比结果；

图8为Cuprite矿区影像数据下估计端元与地物参考端元之间的SAD；

图9为本申请实施例提供的一种高光谱解混装置的功能模块图；

图10为本申请实施例提供的一种电子设备的结构示意图。

具体实施方式

下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行描述。应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

本申请实施例在现有技术中基于NMF模型的高光谱解混方法的基础上进行了改进，提供一种基于广义属性成分分析模型的高光谱解混方法，该方法在NMF模型中使用的第二约束方程的基础上，通过引入丰度矩阵的行稀疏约束形成第一约束方程，并进一步在一些实现方式中提出了可用于求解第一约束方程的迭代算法。该高光谱解混方法求解过程稳定高效，并且解混结果精度较高。

为便于理解发明内容，在介绍本申请实施例提供的高光谱解混方法前，先对高光谱图像的线性光谱混合模型和NMF模型作简单介绍。

线性混合模型

线性光谱混合模型基于以下假设：在宏观尺度下，高光谱图像的混合像元内不同成分物质(即端元)的空间分布是离散均质分布，因此每个瞬时光线只与一种物质进行交互，而不存在不同成分物质之间的多重散射现象。在该假设条件下，光谱混合只是发生在光谱仪设备内部，即由于空间分辨率的限制，使得尽管从每种成分物质反射得到的信号是各自独立的，但是经过光谱仪处理后的像元光谱信号是由各成分物质的光谱反射信号线性混合而来。

图1为线性光谱混合模型的示意图。参照图1，光谱反射在场景中的三种成分物质m1、m2、m3上同时发生，最终得到的像元光谱信号y来自于对每种成分物质反射信号的加权平均，权值设定为每种成分物质在该像元中的丰度，即α1、α2、α3。

假设

代表高光谱图像的第n个混合像元的光谱信号(向量形式)，其中，n＝1,2,…,t，t为高光谱图像中混合像元的总数，b为波段数。基于线性光谱混合模型，x(n)可以表示为：

x(n)＝Ms(n)+w(n) (1)

其中，

为端元光谱矩阵，M的每一列对应一个端元光谱信号，p为端元总数，

为丰度向量，s(n)中的每个元素表示x(n)中相应端元所占的比例，

为模型误差向量。

丰度向量s(n)通常需要满足两个约束条件：

a.非负约束(Abundance Non-negativity Constraint,ANC)，即s_i(n)≥0，i＝1,2,…,p；

b.和为1约束(Abundance Sum-to-one Constraint,ASC)，即

对高光谱图像中所有混合像元的光谱信号都采用(1)进行表示，则(1)中的线性光谱混合模型可以写成矩阵形式：

Y＝MS+W (2)

其中，

为高光谱图像，每列表示一个混合像元的光谱信号，

为丰度矩阵，每列表示一个混合像元的丰度向量，

为模型误差矩阵。

当然，也可以从另一个角度来理解(2)，M的每一列对应一个端元光谱信号，丰度矩阵S＝[s₁ ^T,s₂ ^T,…,s_p ^T]^T的每一行表示一个端元光谱信号对应的丰度图像，从而在线性光谱混合模型中，高光谱图像可以视为端元光谱信号的线性组合，线性组合系数为端元光谱信号的丰度。

NMF模型

NMF模型具有与(2)相同的公式表达，但在表示的含义上有所区别。NMF模型原本用于盲源分离领域，在NMF模型中一个混合信号可以表示成多个源信号的线性组合。将NMF模型应用到高光谱解混领域时，高光谱图像Y视为上述混合信号，丰度矩阵S的每一行(即一个丰度图像)视为一个上述源信号，端元光谱矩阵M用于对源信号进行线性混合，因此也称为混合矩阵。

为了应用NMF模型对(2)进行求解，可以假定W为高斯随机噪声，基于最大似然估计，可以得到NMF模型进行高光谱解混时使用的约束方程：

其中，‖.‖_F表示弗罗贝尼乌斯范数，‖.‖₁表示L1范数，1_p和1_t表示两个元素全为1的列向量，M≥0且S≥0对应线性光谱混合模型中的非负约束，

对应线性光谱混合模型中的和为1约束。

对(3)进行优化求解，计算出的M和S即为高光谱解混结果，但目前研究者提出的求解算法速度较慢，并且收敛性也难以保证。

图2为本申请实施例提供的一种高光谱解混方法的流程图。参照图2，该方法包括：

步骤S100：获取待解混的高光谱图像。

步骤S110：对第一约束方程进行优化求解，获得高光谱图像对应的端元光谱矩阵以及丰度矩阵。

高光谱图像可以通过，但不限于光谱仪等设备进行采集。本申请实施例在NMF模型的基础上进一步提出广义属性成分分析模型用于高光谱解混的分析求解，在高光谱解混中应用该模型的前提和NMF模型类似，都要满足(2)中对高光谱图像的假设。

广义属性成分分析模型

广义属性成分分析模型基于稀疏性和属性多样性先验，它可以作为对属性成分分析的一种泛化。在属性成分分析中，假定原始图像是多个属性成分的线性组合(例如，分块平滑与纹理成分)，其中每个属性成分能够在特定字典下具有稀疏的表达，而且该字典无法对其他属性成分提供稀疏的表达，从而利用字典来完成不同属性成分的分离。具体到高光谱解混问题，属性成分可以是指高光谱图像中隐含的一些高层空间语义信息，例如纹理，结构、形态、分布等特性。

广义属性成分分析模型进一步扩展了属性成分分析框架，并应用于盲源分离目的。该模型假定原始图像Y是由多个源图像通过混合矩阵M线性混合形成(这里由于是在介绍盲源分离，因此Y、M可以作更广义理解，不限于高光谱图像和端元光谱矩阵)，每个源图像又可以进一步建模为多个属性成分的线性组合形式，其中每个属性成分可以用对应的字典进行稀疏表示，如下所示：

其中，s_j表示一个源图像，K是组成该源图像的属性成分数量，s_jk表示源图像s_j的第k个属性成分，α_jk是属性成分s_jk在字典D_k下的稀疏系数。定义联合字典D：＝(D₁,D₂,…,D_k)，则D扮演了不同属性成分的判别器，决定了源图像分离的质量。

能够实现从原始图像中分离源图像的关键在于源图像的稀疏表达和原始图像中属性成分的多样性。所谓稀疏表达，就是说选择合适的字典后，每个源图像在字典下都可以只用少数几个显著的稀疏系数来表示(其他稀疏系数都很小，例如接近于零)。所谓属性成分的多样性(以下简称属性多样性)，就是说存在大量的属性，从而两个源图像属性完全相同的概率极小，即不同的源图像之间存在属性差异，进而它们的稀疏表达系数也不相关，可以认为属性多样性在一定程度上也是形成稀疏性的原因。综合以上两点，可以基于源图像稀疏表达系数的不同，对它们进行分离。相应的约束方程如下：

其中，α为属性成分的稀疏表达系数矩阵，正则项

为属性成分的稀疏约束，λ为阈值。对(5)进行优化求解，计算出的M和α即为盲源分离结果

在广义属性成分分析中，当不同的源图像之间具有一定的不相关性时，基于属性多样性的稀疏表达能增强源图像之间的分离，而且对源图像的表示越稀疏，源图像之间将越容易分离。广义属性成分分析对噪声具有鲁棒性，因为稀疏的源图像通常只有少数显著的稀疏系数，而噪声通常是非稀疏的，因此更容易在稀疏域中分离出来。

在将广义属性成分分析模型应用到高光谱解混中时，高光谱图像Y视为上述原始图像，丰度矩阵S的一行s_j(即一个丰度图像)视为一个上述源图像，端元光谱矩阵M视为上述混合矩阵。发明人长期研究发现，每种地物类别(对应一个端元光谱信号)的丰度图像往往具有不同的空间结构属性(如空间分布、几何形态等)，即满足属性多样性。进一步的，运用稀疏表示理论，当给定一个具有表征能力的字典时，每种地物类别的丰度图像可以进行稀疏表示，同时，由于不同的丰度图像之间存在属性的差异，因此在属性多样性的前提下，各丰度图像可以进行不同的稀疏表达(即不共享显著的稀疏系数)，从而可以基于属性多样性的稀疏表达分离不同地物对应的丰度图像，即实现高光谱解混。在上述分析的基础上，基于广义属性成分分析模型进行高光谱解混的第一约束方程可以定义为：

(6)可以视为在利用NMF模型进行高光谱解混时使用的第二约束方程：

的基础上引入了表征对丰度矩阵S进行行稀疏约束的正则化项

为s_j对应的阈值。

其中，i⁺(.)为非负约束指示函数，定义为

x可以是数值、向量、矩阵等，因此i⁺(M)和i⁺(S)可以等效为对M和S的非负约束，从而(7)和(3)实质上是相同的。注意，(7)中虽然没有体现和为1约束，但这个约束也可以在求解出S后再进行归一化来体现。

对(6)进行优化求解可以得到M和S，即实现了高光谱解混。对比NMF模型中使用的(7)，(7)中由于只引入了非负约束，因此难以得到稳定的解，本申请中的方法在非负约束的基础上引入了稀疏约束，从而可以对端元光谱矩阵和丰度矩阵进行稳定求解。

另外，(6)也可以视为(5)在应用到高光谱解混这一具体问题时产生的一种特殊形式，之前已经分析过，每种地物对应的丰度图像往往具有各自不同的空间结构属性，即同时满足稀疏性和属性多样性，因此将(5)的正则化项中对稀疏系数的约束转化为(6)中对丰度图像的约束是合理的。并且(5)中的稀疏系数α_jk是在变换域，‖Y-MS‖是在直接域，求解起来比较困难，(6)中的||s_j||也是在直接域，求解起来相对容易。

进一步的，发明人研究发现，在不考虑物理意义时，对行稀疏约束做累加和于对整个丰度矩阵做稀疏约束并没有明显的区别，因此在一些实现方式中，为简化求解过程，可以将(6)中的所有λ_j设置为相同的阈值λ，得到：

对(8)进行优化求解可以得到M和S，同样可以实现高光谱解混。

在使用广义属性成分分析模型进行高光谱解混时综合考虑了属性多样性和稀疏性，并将对丰度矩阵的行稀疏约束通过正则化项的形式引入到模型优化求解所使用的第一约束方程中，由于正则化项中既体现了高光谱图像中空间属性的多样性，又体现了丰度图像在表征上的稀疏性，因此该方法可以改善高光谱解混的精度。

下面介绍第一约束方程(7)的优化求解方法，注意以下仅仅提供了一些可能的求解方法而不是唯一的求解方法。

由于第一约束方程是非凸问题，因此在一些实现方式中可以将其拆分为第三约束方程以及第四约束方程两个凸问题，然后进行交替优化求，以获得第一约束方程的稳定解。其中，第三约束方程为：

第四约束方程为：

优化求解过程需要进行多轮次的迭代，直至满足一定的条件结束。交替优化求解即是说在每轮迭代中M和S是交替进行更新的，即在对第三约束方程中的S进行求解时，将M固定不变，在对第四约束方程中的M进行求解时，将S固定不变。

进一步的，第三约束方程中的优化目标是一个不可微的凸函数，可以基于近似分割原理，将其分割为第一约束项和第二约束项两个凸子函数，然后交替优化求解。其中，第一约束项为

第二约束项为

第一约束项是凸可微的，因此可以直接用梯度下降算法求解，而第二约束项是凸不可微的，因此可以用非负软阈值算子进行近似计算。非负软阈值的定义为：

[ST_λ(S)]₊→sign(S)[|S|-λ]₊ (11)

表示对g(S)指向的近似估计。

在求解时，还可以采用FBS算法以加快收敛速度，FBS算法融合了对第一约束项和第二约束项的求解过程，并采用了(8)中对行稀疏约束的简化表示，得到求解第三约束方程的迭代更新公式：

其中，

为迭代后的S，L_M＝||M^TM||_s′为M的谱范数，M^T(MS-Y)为f(S)的梯度计算结果，λ为S对应的阈值，[.]₊表示正象限投影函数，定义为[x]₊＝max(x,0)。

与求解第三约束方程类似的，第四约束方程也分割为第三约束项和第四约束项两个凸子函数。其中，第三约束项为

第四约束项为g(M)＝i⁺(M)，第三约束项是凸可微的，因此可以直接用梯度下降算法求解，而第四约束项是凸不可微的，因此可以利用正象限投影函数进行投影。

在求解时，还可以采用FBS算法以加快收敛速度，FBS算法融合了对第三约束项和第四约束项的求解过程，得到求解第四约束方程的迭代更新公式：

其中，

为迭代后的M，L_S＝||SS^T||_s′为S的谱范数。

对于(12)中使用的阈值λ，在迭代求解的过程中可以始终保持不变。或者，也可以在刚开始迭代时，通过设置较大的阈值针对高光谱图像中的显著特征进行考虑，并在每一轮迭代后，通过减小该阈值以考虑高光谱图像中更细小的特征，从而达到快速调节模型参数，提高方法的收敛速度以及求解稳定性的目的。

例如，在一些实现方式中，可以在优化求解的每一轮迭代后，计算解混形成的残差，定义为：res＝Y-MS，其中，M和S取本轮迭代后的值。然后对该残差的标准差res_std进行估计，并根据获得的估计值减小λ的取值，这种方式下由于对残差的标准差的估计值是不断变换的，所以λ是动态下降的。

下面通过一些基于Matlab语法的伪代码简单总结高光谱解混方法的流程：

输入：

高光谱图像Y；

端元数量p；％可以通过估计获得

阈值计算参数σ

外层最大迭代数MAX_ITER

内层FBS迭代数FBS_ITER

输出：

端元光谱矩阵M和丰度矩阵S

方法主流程：

第一步：初始化M和S，获得M和S合理的初始值

％M的初始值置为Y的p个主成分向量，可以通过主成分分析获得

M＝[PC_1,…,PC_p]；

％初始化过程仅执行两次循环，当然也可以是其他次数

For k＝1:2

％非约束更新

S＝M\Y；

M＝Y/S；

％FBS约束更新

S＝FBS_S(M,S,Y,0,FBS_ITER)；

M＝FBS_M(M,S,Y,0,FBS_ITER)；

End

第二步：计算初始阈值λ

M＝Normalize_M(M,S)；％首先对M执行正则化，并将尺度赋值给S

λ＝‖M^T(MS-Y)‖_∞％计算初始阈值λ，

第三步：进入迭代主循环

For i＝1:MAX_ITER

％基于FBS优化算法对(9)进行优化求解，得到新的S。

S＝FBS_S(M,S,Y,λ)；

％基于FBS优化算法对(10)进行优化求解，，得到新的M。

S＝Normalize_S(M,S)；％对S执行正则化，并将尺度赋值给MM＝FBS_M(M,S,Y,0)；％更新M时，阈值λ设置为0

％对λ做阈值下降

M＝Normalize_M(M,S)；％对M执行正则化，并将尺度赋值给S

λ＝updateLambda(M,S,Y,i)；％阈值下降更新

end

％FBS_S算法的主体，FBS_M类似，就不写出代码了

％前后向分割算法包括两个步骤，首先是一个前向步骤，即gradient算子的工作，

％该步骤在上一次迭代的S上执行梯度下降，即沿着前次S的负梯度方向做下降移动。

％然后是一个proximal近似估计步骤，或者称为后向梯度下降步骤。

％这个近似步骤也是可以看成是在做梯度下降，但是不是在上一次迭代结果上做梯度下降，

％而是在估计得到的最终结果上做梯度下降，所以叫后向步骤。

高光谱解混的结果可以用于矿物填图、亚像元级目标探测、森林资源清查、水体污染与监测、土壤土质调查、农业监测等方面。以亚像元目标探测为例，高光谱遥感影像具有光谱分辨率高、图谱合一的特点，可以提供区分不同物质的诊断性光谱信息，因此在地物目标探测领域具有独特的优势。由于地物分布情况复杂、传感器空间分辨率的限制、目标数目少尺寸小等原因,待探测目标通常与其他地物共同组成混合像元，这时目标以信号比较弱小的亚像元形式存在，问题转化为在像元内部确定目标信号存在性的亚像元小目标探测问题，此时可以通过解混实现对这些物质的探测，即判断混合像元中是否存在某种物质。显然，由于本申请实施例提供的高光谱解混方法提高了解混精度、解混效率以及解混的稳定性，因此有利于高光谱解混技术在上述方面的推广应用，例如，可以更快更精确地实现小目标探测。

为验证本申请实施例提供的高光谱解混方法相对于一些主流解混方法的优越性，发明人进行了大量实验，在后文中予以简要说明。为了简化描述，将本申请实施例提供的高光谱解混算法称为GACA-HU(Generalized Attribute Component AnalysisHyperspectral Unmixing)方法，其具体实现可以参考上述伪代码中的实现方式。

模拟数据实验

模拟数据实验是指通过生成模拟数据对各种高光谱解混方法进行定量评估。比较方法包括顶点成分分析(Vertex Component Analysis，VCA)+完全约束最小二乘法(FullyConstrained Least Squares，FCLS)、最小体积约束NMF(Minimum Volume ConstrainedNMF，MVC-NMF)以及R-CoNMF。VCA+FCLS是经典的二步解混方法，其中，VCA算法用来识别端元光谱，然后利用全约束最小二乘算法来进行丰度估计。MVC-NMF和R-CoNMF均为基于NMF模型的解混方法。其中，MVC-NMF通过施加最小单形体体积约束到端元集合上，而R-CoNMF则施加了一个协同稀疏约束到丰度矩阵上。

实验中，使用了两个测度来进行定量评估：第一个测度是光谱角距离SAD，定义为估计端元

与地物参考端元m的光谱角度，公式如(14)所示，光谱信号越相似，则SAD值越小。

第二个测度是均方根误差RMSE，定义为估计丰度

与地物参考丰度s的均方根，公式如下所示：

在以下的内容中，首先，简单介绍模拟数据的生成方式；然后，给出相应的实验设置；最后，对实验结果进行探讨。

模拟数据生成

模拟数据的生成步骤如下：首先，从USGS数字光谱库中随机挑选一批地物光谱信号。设置模拟图像大小为256×256和221个波段，然后，将图像均匀分割成256个大小为16×16的子块，并给每个子块中的所有像元随机设置一个端元信号。之后，使用邻域窗口大小为33×33的空间低通滤波器对每个图像像元进行滤波处理，从而生成具有局部平滑丰度的混合像元。这里，需要强调的是，滤波窗口越大，像元的混合程度越高。进一步的，为了从图像中完全移除纯净像元，在实验中还对每个像元的丰度进行阈值控制，即如果某个像元的丰度大于参数θ(θ<1)，则将该像元的丰度设置为1/p，其中，p表示生成模拟数据时用到的端元数目。最后，不同信噪比的零均值高斯白噪声被加入模拟数据中来得到最终的模拟实验数据。

实验设置

针对模拟数据实验，发明人从三个角度进行了定量评估：1)噪声鲁棒性；2)像元混合程度；3)端元数目。在考虑噪声鲁棒性时，在实验中生成了不同信噪比的模拟数据，包括10dB、30dB、50dB、70dB以及无噪声几种情况。考虑到高光谱解混方法对像元混合程度的敏感性，发明人选取了两种不同的θ取值(0.6、0.8)，来对每个像元的丰度进行阈值约束。考虑到算法对端元数目的敏感性，发明人分别使用了不同数量的端元来生成模拟图像，包括4、6、8、10。

需要强调的是，由于在GACA-HU方法的迭代优化过程中采用了阈值下降策略，因此方法只有唯一的调节参数用来计算每一次迭代的阈值。即每次迭代时，阈值来自于参数σ乘上当前迭代所估计的噪声标准偏差。发明人研究发现，当信噪比较低时，σ取值往往小于1，从而起到在分解结果中抑制噪声的效果；而当信噪比较高时，σ取值通常在范围[1，10]之间(信噪比越高，σ取值通常越大)。

此外，GACA-HU方法的其他参数在实验中均设置为固定值。具体包括：外层迭代的最大迭代数设置为500，两个内层迭代的最大迭代数均设置为80。M的初始值设置为Y的前p个主成分，S初始值设置为M^TY。然后，通过执行两次交替最小二乘更新和约束FBS算法来初始化M和S。其中，初始化时FBS算法中的正则化参数λ设置为0。然后进入主循环后，初始λ值设置为‖M^T(MS-Y)‖_∞，从而使得S的系数在第一次迭代时保持非增，然后逐渐下降到σ×res_std。

其他解混方法的相关参数设置如下：

对于VCA+FCLS方法，VCA只需要提供端元数目p以及原始高光谱图像Y即可执行，没有其他要设置的参数；FCLS则利用VCA得到的端元集以及原始影像，来反演出丰度。

对于MVC-NMF方法，有一个初始化步骤，采用VCA来得到初始端元，然后采用FCLS来得的初始丰度矩阵。对于MVC-NMF的迭代主循环，则设置80次迭代，10^-6的循环终止条件和0.015的模拟退火温度，以得到一个相对鲁棒的结果。

对于R-CoNMF方法，其参数包括两个正则化参数α和β，两个近似迭代优化子问题中的正则化参数λ_t和μ_t，以及一个代表最小单形体的极限端元集。由于R-CoNMF假定预先不知道端元数目，因此首先使用最小误差高光谱信号辨识(HySime)算法进行端元数目p的估计，然后采用主成分分析保留原始高光谱图像Y的前p个主成分，最后使VCA进行初始端元的估计以及使用FCLS进行初始丰度的估计。对于R-CoNMF的迭代主循环，具体参数设置为：β＝10^-1,α＝10^-8，λ_t＝1，μ_t＝1。此外，针对体积约束采用了VCA得到的端元集作为最小单形体的边界。

最后，整个实验过程是在相同的软硬件平台上运行，其软硬件配置如下：MatlabR2015a桌面版，Intel Core i7 CPU 3.6GHz，32GB内存。

实验结果分析

1)噪声鲁棒性

该实验的目的是评估GACA-HU方法在不同信噪比下的性能。图3和图4分别示出了在采用GACA-HU方法时，基于不同信噪比、不同端元数目以及不同混合程度条件(θ＝0.8为较低混合程度，后文简称低混，θ＝0.6为较高混合程度，后文简称高混)下的SAD和RMSE实验结果。图3的实验结果表明，尽管当信噪比逐渐增大时，SAD逐渐减小，但是实验结果的差别并不明显，特别是当端元数目和混合程度均较低时。随着端元数目的增加，在小于30dB信噪比时的结果与大于30dB的结果差别不大。因此，该实验结果充分说明了本申请提供的GACA-HU方法具有较好的噪声鲁棒性。图4的实验结果具有和图3的实验结果类似的结论，即表明了GACA-HU方法对噪声具有鲁棒性，此处不再重复分析。

为了进一步展现GACA-HU方法的优越性，图5(包括图5(A)至图5(H))显示了各种高光谱解混方法在不同信噪比、端元数目以及不同像元混合程度下的SAD实验结果。实验结果表明，所有方法在信噪比增加时，都得到了逐渐下降的SAD结果，当信噪比较高时，R-CoNMF方法在低混合程度下获得了最好的解混结果。然而，当信噪比逐渐降低时，可以发现噪声对MVC-NMF方法和R-CoNMF方法的影响更大。与之相反，GACA-HU方法表现出更好的噪声鲁棒性。这是因为MVC-NMF方法和R-CoNMF方法均采用了单形体体积约束，该约束对噪声更敏感。

2)像元混合程度

该实验的目的是评估GACA-HU方法在不同像元混合程度下的性能。通常，混合程度越低，SAD值将越小。观察图5(A)、图5(C)、图5(E)以及图5(G)可知，尽管GACA-HU方法在低混情况下没有获得最好的解混结果(R-CoNMF方法在大多数情况下结果最好)，然而，观察图5(B)、图5(D)、图5(F)以及图5(H)可知，当处在高混情况下时，GACA-HU方法通常能得到最好的解混结果。此外，从图3中的定量结果可以观察到GACA-HU方法对混合程度不敏感，特别是在低信噪比情况下的高混数据上的结果，均好于其他算法。总的来说，该实验结果表明GACA-HU方法对混合程度具有更好的泛化性能。考虑到一些真实高光谱图像场景中，往往存在大量高度混合的像元，因此GACA-HU方法的这种泛化能力将带来非常显著的应用价值。

进一步的，可以通过实验获得各种高光谱解混方法在不同混合程度下所得到的丰度图像，例如，该丰度估计结果可以来自于四个端元的合成图像。结果表明，无论是在低混还是高混下，GACA-HU方法得到的丰度图像能更好地保持空间局部同质性，以及地物分布形态特性，而其他方法在空间分布和形态特性上则存在一定程度的模糊。

3)端元数目

该实验的目的是评估GACA-HU方法对不同端元数目的性能。通常，端元数目越大，解混过程越困难。然而，图5的实验结果表明，GACA-HU方法所得到的SAD结果受不同端元数目的影响较小。例如，在8端元和10端元情况下，GACA-HU算法仍然能够得到较好的解混结果。类似的，R-CoNMF方法也表现出对端元数目的稳定性(例如，8端元下的解混结果R-CoNMF方法更好，而10端元下GACA-HU方法则具有更好的解混结果)。此外，在高度混合的情况下，GACA-HU在大端元数目条件下仍然具有最优的SAD值。因此，实验结果说明了GACA-HU方法对端元数目的稳定性要优于其他现有方法。

真实高光谱图像实验

两幅真实的高光谱影像数据被用来进行各种高光谱解混方法的性能评估，一幅图像是采用HYDICE传感器采集的城市场景影像，另一幅图像是采用AVIRIS传感器采集的Cuprite矿区影像。这里所称的高光谱解混方法仍然包括GACA-HU、VCA+FCLS、MVC-NMF以及R-CoNMF几种。

HYDICE城市影像数据

该幅城市影像具有307×307的空间尺寸，以及210个光谱波段，采集自1995年美国德克萨斯州科波拉斯湾。影像的空间分辨率是4平方米，光谱分辨率10纳米，覆盖波段范围从0.4到2.5微米。去除噪声和水汽波段后，保留了162个波段。该幅影像包括4个典型地物端元类：道路、屋顶、树木和草地。其中，树木和草地具有非常相似的光谱信号特征，这为它们的分离带来了较大困难。

图6给出了该幅影像上的SAD实验结果，除了包含平均SAD外，还给出了每类地物的SAD值。从图6可以观察到，GACA-HU方法能获得非常好的SAD结果，特别是对树木和草地的区分能力，而其他方法对这两类地物仍然存在一定程度的偏差。这是因为GACA-HU方法充分利用了地物空间分布的属性多样性特征，能更好地分离光谱相似、但空间分布和几何结构不相似的地物类型。

此外，还可以通过实验获得各种高光谱解混方法得到的端元光谱信号和相应的丰度图。从这些丰度图中可以看出，GACA-HU方法所获得的地物丰度图呈现出非常清晰的空间分布特征。例如，道路表现出独特的现状结构，草地表现出局部块状结构。而屋顶和树木，则以更小的块分散在整个场景中。总的来说，所得到的丰度图像较好地展现了不同地物类别的空间分布和结构特性。

进一步的，图7(包括图7(A)至图7(D))给出了采用GACA-HU方法时估计端元与参考端元光谱曲线的对比结果。通过该结果可以看出，两种曲线匹配良好(即曲线形态接近)，尤其是树木和草地这两类地物。而对于道路和屋顶，一些光谱的扰动仍然存在。

总体来说，实验结果表明了GACA-HU方法在HYDICE城市影像数据的解混中的优势。

Cuprite矿区影像数据

该幅影像采集自1997年的Cuprite矿区，影像具有250×190的空间尺寸，以及224个光谱波段，波段覆盖范围从0.4到2.5微米，达到10纳米的光谱分辨率和20米空间分辨率。去除水吸收波段和低信噪比波段后，保留了188个通道用于实验。Cuprite数据主要包含了各种矿物分布，这些矿物的光谱信号都包含在USGS数字光谱库中。

考虑到该幅影像中矿物分布的复杂性，端元数目非常难以准确地获取。因此，实验采用了12种代表性的矿物端元的情况，这种情况在一些公开算法中已被采用。

图8给出了该幅影像上的SAD实验结果。从图8的结果中可以得出以下结论，对于水铵长石矿物，GACA-HU方法和R-CoNMF方法得到了比MVC-NMF方法和VCA+FCLS方法更小的SAD值，而对于蓝宝石、黄钾铁矾、榍石等矿物类型，GACA-HU方法获得了最小的SAD值。平均SAD的结果也表明了GACA-HU优于其他对比方法。但是也需要承认，由于该场景的复杂性，GACA-HU方法的优势并不明显，这是因为该矿区的各种矿物的混合程度更高，而且更多的是物质之间的精细混合。导致每种矿物的丰度图像呈现出碎片化，空间一致性较弱。这种情况下，空间形态多样性在解混过程中的优势难以充分发挥。

进一步的，还可以通过实验获得采用GACA-HU方法时识别的端元光谱信号与相应的USGS库中最匹配的光谱信号的对比结果。从该结果中可以看到，尽管由于一些矿物类型的复杂性和变异性，但是GACA-HU方法仍然能够较好地识别出这些代表性矿物的光谱信号，从而再次证明了GACA-HU方法的有效性。

图9本申请实施例提供的高光谱解混装置200的功能模块图。参照图9高光谱解混装置200包括：图像获取模块210，用于获取待解混的高光谱图像，所述高光谱图像能够表示为端元光谱信号的线性组合，线性组合系数为所述端元光谱信号的丰度；解混模块，用于对第一约束方程进行优化求解，获得所述高光谱图像对应的端元光谱矩阵以及丰度矩阵；其中，所述第一约束方程包括：利用NMF模型对所述高光谱图像进行解混时使用的第二约束方程，以及，对所述丰度矩阵进行行稀疏约束的正则化项。

在高光谱解混装置200的一些实现方式中，所述第一约束方程为

其中，所述第二约束方程为

所述正则化项为

Y为所述高光谱图像，M为所述端元光谱矩阵，S为所述丰度矩阵，‖.‖_F表示弗罗贝尼乌斯范数，‖.‖₁表示L1范数，p为S的行数，s_j为S的第j行，λ_j为s_j对应的阈值，i⁺(.)为非负约束指示函数，定义为

在高光谱解混装置200的一些实现方式中，对第一约束方程进行优化求解，包括：对由所述第一约束方程拆分而成的第三约束方程以及第四约束方程利用前后向分割FBS算法进行交替优化求解；其中，所述第三约束方程为

所述第四约束方程为

在高光谱解混装置200的一些实现方式中，对所述第三约束方程利用FBS算法进行优化求解，包括：对由所述第三约束方程拆分而成的第一约束项

其中，

在高光谱解混装置200的一些实现方式中，对所述第四约束方程利用FBS算法进行优化求解，包括：对由所述第四约束方程拆分而成的第三约束项

其中，

在高光谱解混装置200的一些实现方式中，在优化求解的每一轮迭代过程中减小λ的取值。

在高光谱解混装置200的一些实现方式中，在优化求解的每一轮迭代过程中减小λ的取值，包括：在优化求解的每一轮迭代后对残差res＝Y-MS的标准差进行估计，并根据获得的估计值减小λ的取值。

本申请实施例提供的高光谱解混装置200，其实现原理及产生的技术效果在前述方法实施例中已经介绍，为简要描述，装置实施例部分未提及之处，可参考方法施例中相应内容。

图10本申请实施例提供的电子设备300的结构示意图。参照图10电子设备300包括：处理器310、存储器320以及通信接口330，这些组件通过通信总线340和/或其他形式的连接机构(未示出)互连并相互通讯。

其中，存储器320包括一个或多个(图中仅示出一个)，其可以是，但不限于，随机存取存储器(RandomAccess Memory，RAM)，只读存储器(Read Only Memory，ROM)，可编程只读存储器(Programmable Read-Only Memory，PROM)，可擦除只读存储器(ErasableProgrammable Read-Only Memory，EPROM)，电可擦除只读存储器(Electric ErasableProgrammable Read-Only Memory，EEPROM)等。处理器310以及其他可能的组件可对存储器320进行访问，读和/或写其中的数据。

处理器310包括一个或多个(图中仅示出一个)，其可以是一种集成电路芯片，具有信号的处理能力。上述的处理器310可以是通用处理器，包括中央处理器(CentralProcessing Unit，CPU)、微控制单元(Micro Controller Unit，MCU)、网络处理器(NetworkProcessor，NP)或者其他常规处理器；还可以是专用处理器，包括数字信号处理器(DigitalSignalProcessor，DSP)、专用集成电路(Application Specific IntegratedCircuits，ASIC)、现场可编程门阵列(Field Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。

通信接口330包括一个或多个(图中仅示出一个)，可以用于和其他设备进行直接或间接地通信，以便进行数据的交互。通信接口330可以是以太网接口；可以是移动通信网络接口，例如3G、4G、5G网络的接口；还是可以是具有数据收发功能的其他类型的接口。

在存储器320中可以存储一个或多个计算机程序指令，处理器310可以读取并运行这些计算机程序指令，以实现本申请实施例提供的高光谱解混方法的步骤以及其他期望的功能。

可以理解，图10示的结构仅为示意，电子设备300还可以包括比图10所示更多或者更少的组件，或者具有与图10示不同的配置。图10所示的各组件可以采用硬件、软件或其组合实现。于本申请实施例中，电子设备300可以是，但不限于专用设备、台式机、笔记本电脑、平板电脑、智能手机、智能穿戴设备、车载设备等实体设备，或者虚拟机等虚拟设备。电子设备300可以是一台设备，例如单台服务器；也可以是多台设备的组合，例如服务器集群。

本申请实施例还提供一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序指令，所述计算机程序指令被计算机的处理器读取并运行时，执行本申请实施例提供的高光谱解混方法的步骤。例如，计算机可读存储介质可以实现为图10电子设备300中的存储器320。

在本申请所提供的实施例中，应该理解到，所揭露装置和方法，可以通过其他的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其他的形式。

另外，作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

再者，在本申请各个实施例中的各功能模块可以集成在一起形成一个独立的部分，也可以是各个模块单独存在，也可以两个或两个以上模块集成形成一个独立的部分。

以上所述仅为本申请的实施例而已，并不用于限制本申请的保护范围，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。