CN108388863A - 一种高光谱遥感图像混合像元分解方法 - Google Patents

Abstract

本发明公开了一种高光谱遥感图像混合像元分解方法。本发明将近似稀疏约束的多尺度非负矩阵分解算法应用于高光谱遥感图像解混，针对现有高光谱遥感图像稀疏解混技术中L0模型求解难的问题，利用一种新的近似稀疏模型替代L0模型进行求解，并在此基础上，考虑到高光谱遥感影像复杂的多尺度空间几何结构，增加了全变差空间约束，使得混合像元的分解精度和性能明显提高。

Description

一种高光谱遥感图像混合像元分解方法

技术领域

本发明涉及遥感图像处理技术，尤其涉及一种高光谱遥感图像混合像元分解方法。

背景技术

由于高光谱数据能同时记录同一个场景中数百个光谱段信息，因此被广泛应用于诸多领域。但是由于高光谱传感器的低空间分辨率和地物的复杂多样性，使得空间相邻物质的光谱不可避免地融合在一起，导致高光谱数据产生混合现象，这种混合像元现象大大影响了地物识别和区分的精度。高光谱解混技术可以对这些混合像元进行分解，使之分解为典型的地物光谱(端元)和对应光谱所占比例(丰度)。

线性高光谱解混技术是一种光谱解混的标准技术，它可以把高光谱中的混合像元分解为光谱端元和端元丰度，其模型为端元的线性组合与端元丰度的加权的形式，并且受到现实情况的约束，需要满足端元和丰度数据的非负性与丰度和为一的条件。由于其模型简单、效率较高、物理含义明确(可参见文献[张绍泉.基于高光谱遥感影像稀疏解混的水域变化检测[D].南昌工程学院,2015.])，并且在一般的情形下能够得到比较满意的效果，是目前混合像元分解中的主流。

由于线性高光谱模型与非负矩阵分解(NMF)[D.D.Lee andH.S.Seung.Algorithmsfor non-negative matrix factorization[J].Advances in Neural InformationProcessing Systems,2001,13:556-562.]模型十分相似，且非负矩阵分解算法已经非常成熟，因此许多学者把非负矩阵分解与高光谱解混结合起来，形成了NMF高光谱解混方法。

稀疏线性高光谱解混技术是利用已知的光谱库代替端元集合，在线性高光谱解混基础上添加了端元的稀疏特性，用尽可能少的非零系数对原始信号进行表示的过程，这能让我们更容易的获取信号中所包含的信息并且简化信号处理问题的求解。

稀疏正则化NMF解混技术把NMF优化分解方法与稀疏正则化计算合二为一，使得高光谱数据在处理性能与精度方面得到了很大的优化。以至于后期，很多学者在此基础上进行了许多改进，其中最典型的是L0稀疏非负矩阵分解方法(可参考文献[高钛.基于近似稀疏约束非负矩阵分解的高光谱图像混合像元分解[D].江西科技师范大学,2014])。

然而，L0范数稀疏性好，但是其高光谱解混存在一个NP问题，求解困难，不适用与实际问题(可参考文献[徐晨光,邓承志.稀疏正则化的非负矩阵分解高光谱解混算法比较与分析[J].南昌工程学院学报,2017,36(6):78-87])，并且现有各类解混模型主要集中于利用高光谱的光谱特性，其像素邻域之间的空间信息却没有得到充分有效的利用。大部分遥感高光谱图片都是针对真实地物进行采样获取的图片，其相邻像素的丰度间具有一定的相似性。如果能对L0稀疏进行改进，并且把空间信息应用到高光谱解混中来，将对其解混结果产生很大的改善作用。

发明内容

本发明所要解决的技术问题在于克服现有技术不足，提供一种高光谱遥感图像混合像元分解方法，针对现有高光谱遥感图像稀疏解混技术中L0模型求解难的问题，利用一种新的近似稀疏模型替代L0模型进行求解；并在此基础上，考虑到高光谱遥感影像复杂的多尺度空间几何结构，增加了多尺度全变差空间约束，使之解混精度和性能明显提高。

本发明具体采用以下技术方案解决上述技术问题：

一种高光谱遥感图像混合像元分解方法，首先估计出高光谱遥感图像中端元的个数；然后基于所估计出的端元个数，利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解，得到高光谱遥感图像的端元矩阵和丰度矩阵；所述近似稀疏约束的多尺度非负矩阵分解方法中，近似稀疏约束模型为：

其中，λ_A0为正则化参数，S_l(pn)表示丰度矩阵在第l尺度的分解尺度因子S_l中第p行第n列的值，σ为近似参数，P为估计出的高光谱遥感图像中端元的个数，N表示高光谱遥感图像中的像素个数。

进一步地，在利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解的过程中，至少在其中一个尺度的第l尺度非负矩阵分解的目标函数中引入了全变差空间约束项λ_TVTV(S_l)，其中，TV(S_l)＝∑_{i,j}∈ε||S_l(i)-S_l(j)||₁，λ_TV为非零的全变差正则化参数， S_l表示丰度矩阵在第l尺度的分解尺度因子，S_l(i)、S_l(j)分别表示在S_l中任意一个元素的同一邻域中位置i和位置j的丰度值，ε表示所述邻域的取值范围。

优选地，所述l层非负矩阵分解的目标函数具体如下：

其中，λ、λ_A0为正则化参数，X_l、A_l分别为混合像元、端元矩阵在第l层非负矩阵分解的分解值，S_l(pn)表示丰度矩阵在第l尺度的分解尺度因子S_l中第p行第n列的值，σ为近似参数，P为估计出的高光谱遥感图像中端元的个数，N表示高光谱遥感图像中的像素个数，上标T为矩阵转置符号。

优选地，在利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解的过程中，其余各尺度非负矩阵分解的目标函数与第l尺度非负矩阵分解的目标函数相同。

优选地，在利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解的过程中，每一尺度非负矩阵分解均使用交替迭代算法对所述目标函数进行求解。

进一步地，所述交替迭代算法的停止条件为：迭代次数达到预设的最大迭代次数，或者，前后两次迭代的误差不超过预设误差容忍度的情形连续出现次数达到预设次数。

优选地，在进行第一尺度非负矩阵分解时，使用VCA-FCLS算法进行端元矩阵和丰度矩阵的初始化；在其余各尺度非负矩阵分解时，使用随机方法进行端元矩阵和丰度矩阵的初始化。

优选地，正则化参数λ的计算公式如下：

式中，t表示循环迭代中迭代的次数，λ₀和τ是用来调节λ的影响因数；

正则化参数λ_A0、λ_TV的计算公式如下：

式中，L为所述近似稀疏约束的多尺度非负矩阵分解方法中非负矩阵分解的尺度数，x_l表示高光谱遥感图像中第l维的光谱数据。

优选地，在进行S_l更新时，通过对X_l和A_l进行下列扩展操作来实现丰度的和唯一约束：

其中，参数δ是调节矩阵与和为一约束的强弱项。

优选地，使用最小误差的高光谱信号端元数确定Hysime方法估计出高光谱遥感图像中端元的个数。

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

本发明将近似稀疏约束的多尺度非负矩阵分解算法应用于高光谱遥感图像解混，针对现有高光谱遥感图像稀疏解混技术中L0模型求解难的问题，利用一种新的近似稀疏模型替代L0模型进行求解，并在此基础上，考虑到高光谱遥感影像复杂的多尺度空间几何结构，增加了全变差空间约束，使得混合像元的分解精度和性能明显提高。

附图说明

图1a～图1c为Jasper Ridge高光谱数据。

具体实施方式

针对现有高光谱遥感图像稀疏解混技术中L0模型求解难的问题。本发明的思路是利用一种新的近似稀疏模型替代L0模型进行求解，并在此基础上，考虑到高光谱遥感影像复杂的多尺度空间几何结构，增加了全变差空间约束，使得混合像元的分解精度和性能明显提高。

具体地，本发明高光谱遥感图像混合像元分解方法，首先估计出高光谱遥感图像中端元的个数；然后基于所估计出的端元个数，利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解，得到高光谱遥感图像的端元矩阵和丰度矩阵；所述近似稀疏约束的多尺度非负矩阵分解方法中，近似稀疏约束模型为：

其中，λ_A0为正则化参数，S_l(pn)表示丰度矩阵在第l尺度的分解尺度因子S_l中第p行第 n列的值，σ为近似参数，P为估计出的高光谱遥感图像中端元的个数，N表示高光谱遥感图像中的像素个数。

进一步地，在利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解的过程中，至少在其中一个尺度的第l尺度非负矩阵分解的目标函数中引入了全变差空间约束项λ_TVTV(S_l)，其中，TV(S_l)＝∑_{i,j}∈ε||S_l(i)-S_l(j)||₁，λ_TV为非零的全变差正则化参数，S_l表示丰度矩阵在第l尺度的分解尺度因子，S_l(i)、S_l(j)分别表示在S_l中任意一个元素的同一邻域中位置i和位置j的丰度值，ε表示所述邻域的取值范围。

上述技术方案中，可以将现有稀疏解混技术中的目标函数与所述的新的稀疏约束模型和多尺度全变差空间约束项相结合以得到引入了多尺度全变差空间约束的近似稀疏目标函数；本发明优选采用以下的目标函数：

其中

TV(S_l)＝∑_{i，j}∈ε||S_l(i)-S_l(j)||₁ (2)

在公式(1)中λ、λ_A0和λ_TV为正则化参数，X_l、A_l、S_l、为混合像元、端元矩阵与丰度矩阵的第l次分解值，S_l(pn)表示矩阵(即丰度矩阵在第l尺度的分解尺度因子)S_l的第p 行第n列的值，σ为近似参数，σ的值越小的话近似函数的性质越接近L0范数。TV(S_l) 是S_l全变差空间约束项，表示在相同一端元的邻域像素的平滑过渡，由公式(2)表示， (S_l(i)-S_l(j))表示S_l中任意一个元素的同一邻域中位置i和位置j的丰度差值，ε表示所述邻域的取值范围。

在实际应用中，可以通过设置λ_TV值得大小来控制其空间平滑影响权重，而当参数λ_TV为零时，公式(1)变成普通的稀疏约束方法。为了计算方便，我们把公式(2)进行如下优化：

其中H_hS_l＝[h₁，h₂，...，h_n]，h_i＝x_i-x_ih，x_i和x_ih分别表示矩阵中数值与其横向邻域的差值。 H_vS_l＝[v₁，v₂，...，v_n]，v_i＝x_i-x_iv，x_i和x_iv分别表示矩阵中数值与其纵向邻域的差值。因此我们可以把公式(1)写成

其中

把公式(3)的约束项和公式融合到一起，这样方便后续的计算：

其中i_-R₊(S_l)项代表S_l≥0，i_S(S_l)项代表i_-R₊(A_l)代表A_l≥0。

由于限定项中S_l≥0，所以|S_l(pn)|/σ²项必大于0，这样，我们可以把公式(5)简化成：

对于多尺度的非负矩阵分解算法，在每一尺度的计算，都相当于进行一次普通的非负矩阵分解算法。因此在多尺度非负矩阵分解计算中，可以采用相同的目标函数和相同的迭代规则，也可以采用完全不同的目标函数和相应的迭代规则来进行求解。为了计算方便和程序易于编写，本发明优选在各尺度非负矩阵分解算法中采用相同的目标函数和相同的迭代规则来进行求解。

交替迭代算法(ADMM)[M.Afonso,J.Bioucas-Dias,and M.Figueiredo,“Anaugmented Lagrangian approach to the constrained optimization formulation ofimaging inverse problems,”IEEE Trans.Image Process.,vol.20,no.3,pp.681–695,Mar.2011.]整合了许多经典优化思路，结合现代统计学习所遇到的问题，是一个比较好实施的分布式计算框架。因此本发明优选采用ADMM作为各尺度非负矩阵分解算法中的迭代算法。

把公式(6)利用ADMM算法进行优化，并且把约束项全部列入到公式中去，得到如下公式：

其中：

根据ADMM算法，可以得到如下

其中：

在公式(9)中μ为正常数，表示归一化拉格朗日乘子。d_(i+1)是d_(i)与对应优化项误差的差值。根据交替迭代算法(AMDD)算法，如下表:

表1交替迭代算法(AMDD)算法伪代码

我们利用公式(9)并结合表1中公式进行求解可得：

当把基于近似稀疏约束的多尺度非负矩阵分解算法直接应用到高光谱图像的混合像元分解之前，需要对算法进行初始化和相关参数设置。

对高光谱图像进行混合像元的分解，必须先估计出高光谱图像中端元的个数。在基于非负矩阵分解的高光谱图像混合像元分解中运用的最多的两种算法是虚拟维度方法(VD)，最小误差的高光谱信号端元数确定方法(Hysime)。本发明优选采用最小误差的高光谱信号端元数确定(Hysime)来估计高光谱图像中端元的个数。

在基于非负矩阵分解的高光谱图像混合像元分解算法中用的比较多的两种初始化方法是随机方法(随机产生两个大小为R^L×P的端元矩阵和R^P×N的丰度矩阵，其中两个矩阵中的值都在0和1之间)和VCA-FCLS方法(用VCA算法来计算出端元矩阵A，然后用FCLS算法来求解出丰度矩阵S作为初始值)。对于这两种方法而言，优选采用 VCA-FCLS做初始化，从而可以加快算法的计算速度，减少计算时间，并且效果更好。

由于本发明方法需要做多尺度的非负矩阵分解，因此需要多次进行初始化。本发明优选的方法是：第一尺度的非负矩阵分解算法中用VCA-FCLS来对矩阵A₁和矩阵S₁进行初始化；而对于其他尺度的计算我们选择随机方法进行初始化。

根据以上公式，我们可以得出本发明方法(简称AL₀TV-MLNMF)的算法伪代码，如下表

表2AL₀-MLNMF算法

为了便于公众理解本发明技术方案并验证其技术效果，下面以真实的高光谱数据Jasper Ridge为例来对本发明方法进行详细说明。

高光谱数据Jasper Ridge像素的大小为512×614，每个像素包含了224个波段，波长范围是380nm到2500nm，光谱的分辨率为9.46nm。考虑到高光谱数据过于复杂，利用了原始图像中从(105,269)像素位置开始的100×100的像素大小进行分解，并且考虑到密集的水蒸气和大气效应吸收的波段，过滤掉(1-3,108-112,154-166,220-224)波段，这样最后留下198个波段做下一步处理。在高光谱数据中包含了四个端元，分别是道路、泥沙、水、树木。图1a～图1c为Jasper Ridge的地面真实数据，其中图1a为真实图片，图1b为各端元丰度分布图，图1c为各端元光谱信息图。

本发明的解混步骤具体如下：

步骤1、估计端元的数目P

在做高光谱图像的混合像元分解工作前，第一步就是要预估出被测地物的端元数目 P。现有基于非负矩阵分解的高光谱混合像元分解常用的几种估算方法为虚拟维度(VD)法和高光谱最小误差端元确定(Hysime)法。本实例中使用高光谱最小误差端元确定法(Hysime)来对端元数目P进行预估。

步骤2、初始化

预估完端元数目，就要对端元矩阵与丰度矩阵进行初始化工作。目前通常基于非负矩阵分解的高光谱混合像元解换常用的几种初始化方法为VCA-FCLS初始化方法和随机初始化方法。VCA-FCLS初始化方法是首先采用VCA方法对端元矩阵进行预估，然后利用预估出来的端元矩阵和FCLS方法求出丰度矩阵，作为初始值。随机初始化方法是随机产生一个R^L×P大小的端元矩阵和一个R^P×N大小的丰度矩阵，两个矩阵中值为0 到1中随机值。这两种方法，随机初始化方法虽然简单但是效果没有VCA-FCLS好，因此本实例中采用VCA-FCLS方法作为初始化方法，不但可以提高计算速率，并且计算出来的数据效果好。

在多尺度非负矩阵分解的算法中，每一尺度算法都是在做一次简单的非负矩阵分解。而在简单的非负矩阵分解中都需要进行初始化工作，因此在多尺度非负矩阵分解算法中，每个尺度既可以使用相同方法进行初始化，也可以使用不同的方法进行初始化。本实例中为了即提高精度又使得计算速度加快，使用混合的初始化方法。即在第一尺度使用VCA-FCLS方法，之后其他尺度的初始化选择随机初始化方法。

步骤3、相关参数设置

对本发明AL₀TV-MLNMF算法的模型公式(7)中三个正则化参数λ、λ_A0、λ_TV进行设置。

通常参数λ设置为如下：

公式中t表示循环迭代中程序迭代的次数，λ₀和τ用来调节λ影响因数。

公式中参数λ_A0和λ_TV，本发明为增强其算法稳定性采用自适应参数。依靠丰度矩阵稀疏度来决定参数λ_A0和λ_TV的值。而由于丰度矩阵又不可预知，采用高光谱数据的图像数据来替代。对参数λ_A0和λ_TV，我们采用了L₁范数与L₂范数相结合的方式来度量其稀疏度。具体公式如下：

公式中x_l表示了被测地物高光谱数据第l维的光谱数据。

在本实例中设置的多尺度非负矩阵分解总尺度数L为50，每层的最大迭代次数T为150。

步骤4、非负数约束、和为一约束条件设置

在高光谱混合像元分解中存在着两个重要的约束条件：非负约束(ANC)和丰度和为一约束(ASC)。由于非负矩阵分解的特点就是分解出两个非负的矩阵，其非负性可以保证。因此只需要保证丰度的和唯一约束。本实例采用了文献[M.P.Bertsekas, ConstrainedOptimization and Lagrange Multiplier Methods.New York:Academic,1982.]中提到的方法进行简单有效的约束，该方法广泛的应用于高光谱混合像元分解中。其数学表达公式如(23)，在对公式(12)第l尺度进行迭代矩阵S_l更新时，对X_l和A_l进行下列扩展操作：

参数δ是调节矩阵与和为一约束的强弱项，其值越大时，其结果效果越好，但其越难收敛；当其值越小时，效果会变差，但其越容易收敛。因此为了使分解效果与收敛度达到平衡，在本实例中设置δ值为20。

步骤5、停止条件设置

通常迭代方法采用的停止条件有两种：第一种是通过设定的迭代次数来控制的，程序如果运行超过了设定值就停止；第二种是通过设定前后两次迭代误差值来控制，在程序运行时不断比较前后两次的误差，如果小于给定的停止条件值足够次数的话就停止迭代，其前后两次误差计算表达公式如下：

||f(A_l,S_l)_new-f(A_l,S_l)_old||＜ε (24)

公式中f(A_l,S_l)_new和f(A_l,S_l)_old分别表示在第t次的迭代目标函数值与第t-1次的迭代目标函数值，ε表示误差的容忍度。

本实例中采用了这两种相结合的方式，设置最大迭代次数为500，以达到最大迭代次数，或者连续10次前后两次误差不超过误差容忍度10^-4为迭代结束条件。

步骤6、迭代

对于多尺度的非负矩阵分解算法，在每一尺度的计算，都相当于进行一次普通的非负矩阵分解算法。因此在多尺度非负矩阵分解计算中，可以采用相同的目标函数和相同的迭代规则也可以采用完全不同的目标函数和相应的迭代规则来进行求解。为了计算方便和程序易于编写，在本实施例中采用相同的目标函数和相同的迭代规则来进行求解。采用交替迭代算法(ADMM)算法进行求解，求解步骤：

1：令X₁＝X

2：当l＝1到L时进行3-5步循环

3：对参数进行初始化：设置的初始值。

4：重复

(a)利用公式(11)对矩阵进行更新

(b)利用公式(12)对矩阵进行更新，并且利用步骤4方法设置和为一非负约束

(c)利用公式(13)-(20)对矩阵进行更新；

(d)利用公式(10)对矩阵D进行更新；

直到达到步骤5的停止条件

5：令X_l+1＝S_l

6：结束循环

步骤7、输出解混结果

迭代之后，根据迭代数据可以得到解混结果。端元矩阵为A＝A₁A₂…A_L，丰度矩阵为S＝S_L。其中A_x为多尺度分解的结果。S_L为最后一个尺度分解的结果。将采用上述方法对JasperRidge进行高光谱数据混合像元分解后得到解混结果与几种现有算法进行比较(L_1/2-NMF，AL₀-NMF，AL₀-MLNMF)，可以发现，相比现有技术，本发明的 AL₀TV-MLNMF算法具有较明显的优势。

Claims (10)

1.一种高光谱遥感图像混合像元分解方法，首先估计出高光谱遥感图像中端元的个数；然后基于所估计出的端元个数，利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解，得到高光谱遥感图像的端元矩阵和丰度矩阵；其特征在于，所述近似稀疏约束的多尺度非负矩阵分解方法中，近似稀疏约束模型为：

其中，λ_A0为正则化参数，S_l(pn)表示丰度矩阵在第l尺度的分解尺度因子S_l中第p行第n列的值，σ为近似参数，P为估计出的高光谱遥感图像中端元的个数，N表示高光谱遥感图像中的像素个数。

2.如权利要求1所述方法，其特征在于，在利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解的过程中，至少在其中一个尺度的第l尺度非负矩阵分解的目标函数中引入了全变差空间约束项λ_TVTV(S_l)，其中，TV(S_l)＝∑_{i,j}∈ε||S_l(i)-S_l(j)||₁，λ_TV为非零的全变差正则化参数，S_l表示丰度矩阵在第l尺度的分解尺度因子，S_l(i)、S_l(j)分别表示在S_l中任意一个元素的同一邻域中位置i和位置j的丰度值，ε表示所述邻域的取值范围。

3.如权利要求2所述方法，其特征在于，所述第l尺度非负矩阵分解的目标函数具体如下：

其中，λ、λ_A0为正则化参数，X_l、A_l分别为混合像元、端元矩阵在第l层非负矩阵分解的分解值，S_l(pn)表示丰度矩阵在第l尺度的分解尺度因子S_l中第p行第n列的值，σ为近似参数，P为估计出的高光谱遥感图像中端元的个数，N表示高光谱遥感图像中的像素个数，上标T为矩阵转置符号。

4.如权利要求3所述方法，其特征在于，在利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解的过程中，其余各尺度非负矩阵分解的目标函数与第l尺度非负矩阵分解的目标函数相同。

5.如权利要求4所述方法，其特征在于，在利用近似稀疏约束的多尺度非负矩阵分解方法进行混合像元分解的过程中，每一尺度非负矩阵分解均使用交替迭代算法对所述目标函数进行求解。

6.如权利要求5所述方法，其特征在于，所述交替迭代算法的停止条件为：迭代次数达到预设的最大迭代次数，或者，前后两次迭代的误差不超过预设误差容忍度的情形连续出现次数达到预设次数。

7.如权利要求4所述方法，其特征在于，在进行第一尺度非负矩阵分解时，使用VCA-FCLS算法进行端元矩阵和丰度矩阵的初始化；在其余各尺度非负矩阵分解时，使用随机方法进行端元矩阵和丰度矩阵的初始化。

8.如权利要求4所述方法，其特征在于，正则化参数λ的计算公式如下：

式中，t表示循环迭代中迭代的次数，λ₀和τ是用来调节λ的影响因数；

正则化参数λ_A0、λ_TV的计算公式如下：

式中，L为所述近似稀疏约束的多尺度非负矩阵分解方法中非负矩阵分解的尺度数，x_l表示高光谱遥感图像中第l维的光谱数据。

9.如权利要求4所述方法，其特征在于，在进行S_l更新时，通过对X_l和A_l进行下列扩展操作来实现丰度的和唯一约束：

其中，参数δ是调节矩阵与和为一约束的强弱项。

10.如权利要求1～9任一项所述方法，其特征在于，使用最小误差的高光谱信号端元数确定Hysime方法估计出高光谱遥感图像中端元的个数。