CN105825227B - 一种基于mfocuss和低秩表示的高光谱图像稀疏解混方法 - Google Patents

Abstract

本发明公开了一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法，本发明具体步骤包括：读取原始高光谱数据和已知光谱库数据；构建基于MFOCUSS和低秩表示的稀疏解混模型的目标函数，并将高光谱数据以及已知光谱库作为该目标函数的输入数据和字典，通过求解MFOCUSS和低秩表示模型的目标函数得到整个光谱库的丰度矩阵、剔除光谱库中非真正端元的光谱、将剔除非真正端元后的光谱库反复迭代最终得到真正端元矩阵以及所对应的丰度矩阵；该方法避免了从原始高光谱数据中直接提取端元，并剔除了非真正端元和更新了光谱库，降低了光谱库中端元光谱的自相关性对高光谱解混效果的不利影响，提高丰度估计的精确性。

Description

一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法

技术领域

本发明属于图像处理技术领域，涉及稀疏解混技术领域中的一种基于多点欠定系统聚焦求解(Multiple FOCal Underdetermined System Solver，MFOCUSS)和低秩表示的高光谱图像稀疏解混方法。

背景技术

高光谱遥感图像是由几十乃至数百个连续波段图像组成的三维图像数据，能够将确定地物空间特性的图像与确定地物性质的光谱信息有机结合，可定性定量地对被测对象进行物理分析和识别，因此在国土资源利用、地质勘探、环境监测、海洋资源普查、人工目标识别以及军事探测等领域得到广泛的应用。但受光学器件的限制，高光谱图像中像元通常包含多于一种地物类型，从而形成混合像元。由于混合像元的普遍存在，导致分类和目标探测等高光谱遥感图像应用研究难以达到实用的要求，因此如何有效进行高光谱解混，成为近年来遥感领域的一个研究热点。

在基于线性混合模型的高光谱混合像元分解中，通常采用基于统计的方法和基于几何的方法等实现高光谱解混。近年来，随着压缩感知和稀疏表示理论的蓬勃发展，开始有学者尝试将稀疏性约束加入到高光谱混合像元分解模型中，并用已知的光谱库代替端元集合用于高光谱混合像元分解，形成了第三种方法——稀疏解混方法。稀疏解混方法主要有凸优化算法和贪婪算法等，在凸优化算法方面，Marian-Daniel Iordache,JoséM.Bioucas-Dias和Antonio Plaza在“Sparse Unmixing of Hyperspectral Data”([J].IEEE Transactions on Geoscience Remote Sensing,2011,49(6).)中提出一种基于变量分裂和增广拉格朗日(Sparse Unmixing by variable Splittingand AugmentedLagrangian)的稀疏解混方法，建立了稀疏解混模型，将稀疏性约束加入到高光谱混合像元分解模型中，用已知的光谱库代替端元集合用于高光谱混合像元分解，从而将解混问题转化为求解一个凸优化问题，并且利用交替乘子法求解该稀疏解混模型，但是该方法存在求解速度较慢，而且由于高光谱数据中的实际端元数相对于光谱库中的光谱曲线是非常小的，不能直接从光谱库中提取端元来控制稀疏解混的稀疏性，从而在一定程度上影响了该方法的稀疏解混精度。

发明内容

本发明的目的在于克服上述现有技术的不足，提出一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法，降低光谱库中端元光谱曲线的相关性对高光谱解混效果的不利影响，以提高高光谱图像稀疏解混的精度，克服高光谱图像稀疏解混效率低的问题。

为实现上述目的，本发明的一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法，该方法步骤如下：

步骤1：读取高光谱数据Y∈R^L×N和光谱库A∈R^L×m，其中L表示高光谱数据的波段数，N表示高光谱数据样本的总个数，m表示表示高光谱库中的物质光谱总个数，R表示实数域；

步骤2：对高光谱数据进行解混，构建基于MFOCUSS和低秩表示的稀疏解混模型的目标函数，并将高光谱数据Y以及已知光谱库A作为该目标函数的输入数据和字典，求解得到对应的丰度矩阵X；

步骤3：分别输出丰度矩阵X和剔除非真正端元后的光谱库A，取出剔除非真正端元后的光谱库A中列标号所对应的列，即为真正的端元光谱，取出丰度矩阵中行标号所对应的行，即为真正的端元光谱所对应的丰度图。

进一步的，所述步骤2中的度矩阵X为：

其中，λ和β为拉格朗日参数，其大小根据实际高光谱数据设置；X∈R^m×N为光谱库A的丰度矩阵；MFOCUSS约束项表示取l₂范数的p次方操作；低秩表示项||X||_*为矩阵X的核范数；X≥0是丰度矩阵X的非负性约束；min表示取最小值操作。

更进一步的，所述步骤2具体指：

步骤2.1)初始化丰度矩阵X以及辅助变量V₁，V₂，D₁，D₂，将求解稀疏解混模型的迭代次数k设置为1，求解稀疏解混模型的最大迭代次数设置为600；

步骤2.2)迭代变量X，V₁，V₂，D₁，D₂，交替迭代过程如下：

其中diag(·)为对角矩阵函数，I为单位矩阵，U∑V^T是矩阵的奇异值分解，函数(x)_＋＝max{x,0}，max表示取最大值操作；

步骤2.3)剔除非真正端元，并更新光谱库A：

判断光谱库中每个物质端元A_i(A_i表示光谱库A中第i个端元)所对应的丰度权重系数gamma_i(i＝1,2,...,m)是否小于阈值ω，阈值ω大小根据实际高光谱数据设置，若gamma_i＜ω则认为物质端元A_i为非真正端元，剔除光谱库中A所对应的端元A_i，否则保留光谱库中A所对应的端元A_i，将剔除非真正端元后的光谱库A更新；

端元A_i所对应的丰度权重系数gamma_i计算如下：

步骤2.4)停止迭代条件1：判断丰度矩阵的收敛条件是否满足，若是，则执行步骤3，否则执行步骤2.5；

丰度矩阵的收敛条件为：

取ε＝10^-6；

2.5)停止迭代条件2：将求解稀疏解混模型的迭代次数k加1，判断迭代次数k是否小于最大迭代次数，若是，则执行步骤3，否则，执行步骤2.2。

本发明与现有方法相比具有如下优点：

第一，本发明由于引入了基于MFOCUSS和低秩表示的稀疏解混模型，在高光谱数据稀疏解混过程中能够利用丰度矩阵的全局结构进行高光谱解混，克服了现有技术在高光谱数据稀疏解混过程中，稀疏解混精度低的问题，使本发明具有高光谱数据稀疏解混精度高的优点。

第二，本发明利用现有的光谱库，避免了直接从原始数据中提取端元，并在解混处理的优化迭代过程中，对光谱库进行逐步剔除非真正端元和更新该光谱库，降低了光谱库中端元自相干性对解混效果的不利影响。

附图说明

图1为本发明的流程图；

图2为本发明仿真采用的模拟数据高光谱图像及其所含九种端元的真实丰度图，(a)-(i)分别为九种端元的真实丰度图。

具体实施方式

本发明提供一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法，为使本发明的目的，技术方案及效果更加清楚，明确，以及参照附图并举实例对本发明进一步详细说明。应当理解，此处所描述的具体实施仅用以解释本发明，并不用于限定本发明。

下面结合附图1、2，对本发明的实现步骤做进一步详细描述。

本发明提出一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法，其流程图见图1所示，该解混方法包括如下步骤：

步骤1：用计算机读取数据：用MATLAB R2009b读取高光谱数据Y∈R^L×N和光谱库A∈R^L×m，其中L表示高光谱数据的波段数，N表示高光谱数据样本的总个数，m表示表示高光谱库中的物质光谱总个数，R表示实数域；

步骤2：高光谱数据解混处理。得到高光谱数据Y和光谱库A，对高光谱数据进行解混；假设光谱库A中含有足够丰富的物质光谱，那么高光谱数据Y中含有的端元只占光谱库中物质光谱的极少一部分，因此高光谱数据Y中包含的端元数一般远小于光谱库A中的端元数，高光谱数据Y可以由光谱库A中少数端元通过不同的混合系数线性组合而成，因此高光谱数据中端元在光谱库A中是稀疏的，他们所对应的丰度矩阵也具有稀疏性。因此，构建基于MFOCUSS和低秩表示的稀疏解混模型的目标函数，并将高光谱数据Y以及已知光谱库A作为该目标函数的输入数据和字典，求解得到对应的丰度矩阵X：

其中，λ和β为拉格朗日参数，其大小根据实际高光谱数据设置；X∈R^m×N为光谱库A的丰度矩阵；MFOCUSS约束项表示取l₂范数的p次方操作；低秩表示项||X||_*为矩阵X的核范数；X≥0是丰度矩阵X的非负性约束；min表示取最小值操作。

将上式转化为无约束优化问题：

其中，表示取F范数的平方操作，l_R+(·)为丰度矩阵的非负约束函数，定义为||X||_*为矩阵X的核范数，其具体定义如下：

其中σ_i(1≤i≤min{m,N})是矩阵X的奇异值。

针对上述目标函数，采用拉格朗日交替乘子法求解上述函数极值，从而得到上述优化模型的最优解——最优丰度矩阵X；

引入辅助变量V₁和V₂，并且V₁＝X,V₂＝X，上式目标函数可以转化为：

转化为无约束形式为：

其中μ为正则化参数，其大小根据实际高光谱数据设置。引入误差反馈变量D₁和D₂，上述优化问题转化为迭代求解过程，如下：

具体求解步骤如下：

2.1)初始化丰度矩阵X以及辅助变量V₁，V₂，D₁，D₂，将求解稀疏解混模型的迭代次数k设置为1，求解稀疏解混模型的最大迭代次数设置为600；

2.2)迭代变量X，V₁，V₂，D₁，D₂，交替迭代过程如下：

其中diag(·)为对角矩阵函数，I为单位矩阵，U∑V^T是矩阵的奇异值分解，函数(x)₊＝max{x,0}，max表示取最大值操作。

2.3)剔除非真正端元，并更新光谱库A：

判断光谱库中每个物质端元A_i(A_i表示光谱库A中第i个端元)所对应的丰度权重系数gamma_i(i＝1,2,...,m)是否小于阈值ω，阈值ω大小根据实际高光谱数据设置，若gamma_i＜ω则认为物质端元A_i为非真正端元，剔除光谱库中A所对应的端元A_i，否则保留光谱库中A所对应的端元A_i，将剔除非真正端元后的光谱库A更新。端元A_i所对应的丰度权重系数gamma_i计算如下：

2.4)停止迭代条件1：判断丰度矩阵的收敛条件是否满足，若是，则执行步骤3，否则执行步骤2.5。丰度矩阵的收敛条件为：

这里取ε＝10^-6。

2.5)停止迭代条件2：将求解稀疏解混模型的迭代次数k加1，判断迭代次数k是否小于最大迭代次数600，若是，则执行步骤3，否则，执行步骤2.2。

步骤3：输出解混结果：分别输出丰度矩阵X和剔除非真正端元后的光谱库A。取出剔除非真正端元后的光谱库A中列标号所对应的列，即为真正的端元光谱，取出丰度矩阵中行标号所对应的行，即为真正的端元光谱所对应的丰度图。

下面结合仿真实验对本发明的效果做进一步的说明。

1.仿真条件

本发明所用的实验数据为文献“Total variation spatial regularization forsparse hyperspectral unmixing”([J].IEEE Transactions on Geoscience RemoteSensing,2012,50(11).)中构造的模拟高光谱图像，该模拟高光谱图像包含100×100个像素点，224个波段，波长范围为0.4到2.5微米。采用信号重构误差公式评估本发明的性能，高光谱的丰度矩阵的重构误差越大

光谱库A是从美国地质勘探局USGS的光谱库splib06中选择498条光谱曲线构成。实验硬件环境:CPU(Intel(R)Core(TM)i5-3470，3.2GHz)，内存(2GB)，硬盘(7200转)；软件环境Matlab R2009a。

2.仿真内容：

仿真内容主要包括本发明与现有几种方法在模拟数据上的，现有的几种方法包括：稀疏法SUnSAL，稀疏结合全变差空间正则SUnSAL-TV这两种方法；本发明的基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法缩写为MFOCUSS-LRR。

3.仿真实验结果分析：

表1展示的是模拟数据在SUnSAL，SUnSAL-TV和本发明MFOCUSS-LRR在信噪比分别为20dB、30dB、40dB时的丰度矩阵信号重构误差(SRE)。定义用dB表示为：SRE(dB)≡10log₁₀(SRE)，其中X为真实的丰度矩阵，为通过解混算法求得的丰度矩阵，E[·]为求期望操作。

表1 不同方法在不同信噪比时的SRE(dB)比较

从表1可以看出，本发明在20～40dB信噪比下都要比SUnSAL和SUnSAL-TV方法能够获得更高的SRE值。可以看出随着SNR的提高，各种算法的信号重构误差SRE逐步提高。由于模拟数据的相邻像素之间具有较强的空间相关性，SUnSAL-TV方法因为加入了全变差正则项约束考虑到了空间相关性，SUnSAL-TV方法比SUnSAL方法在性能上有了较大程度的提高，而本发明在字典更新中逐步剔除非真正端元，消除了非真正端元对真正端元的影响，获得了比SUnSAL-TV方法更高的SRE值。

表2 不同方法在不同信噪比时的运行时间(s)比较

从上表中可以看出，本发明在不同信噪比下的运行时间都要比现有技术中的SUnSAL和SUnSAL-TV技术的运行时间短，说明了本发明采用交替乘子法求解基于MFOCUSS和低秩表示的稀疏解混模型，使得本发明具有高光谱数据稀疏解混耗时短、效率高的优点。

Claims (3)

1.一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法，其特征在于：该方法步骤如下：

步骤1：读取高光谱数据Y∈R^L×N和光谱库A∈R^L×m，其中L表示高光谱数据的波段数，N表示高光谱数据样本的总个数，m表示高光谱库中的物质光谱总个数，R表示实数域；

步骤2：对高光谱数据进行解混，构建基于MFOCUSS和低秩表示的稀疏解混模型的目标函数，并将高光谱数据Y以及已知光谱库A作为该目标函数的输入数据和字典，求解得到对应的丰度矩阵X；剔除光谱库中非真正端元的光谱、将剔除非真正端元后的光谱库重新代入MFOCUSS和低秩表示模型的目标函数计算所对应的丰度矩阵；

步骤3：分别输出丰度矩阵X和剔除非真正端元后的光谱库A，取出剔除非真正端元后的光谱库A中列标号所对应的列，即为真正的端元光谱，取出丰度矩阵中行标号所对应的行，即为真正的端元光谱所对应的丰度图。

2.根据权利要求1所述的一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法，其特征在于：

所述步骤2中的丰度矩阵X为：

其中，λ和β为拉格朗日参数，其大小根据实际高光谱数据设置；X∈R^m×N为光谱库A的丰度矩阵；MFOCUSS约束项0＜p≤1，表示取l₂范数的p次方操作；低秩表示项||X||_*为矩阵X的核范数；X≥0是丰度矩阵X的非负性约束；min表示取最小值操作。

3.根据权利要求2所述的一种基于MFOCUSS和低秩表示的高光谱图像稀疏解混方法，其特征在于：

所述步骤2具体指：

步骤2.1)初始化丰度矩阵X以及辅助变量V₁，V₂，D₁，D₂，将求解稀疏解混模型的迭代次数k设置为1，求解稀疏解混模型的最大迭代次数设置为600；步骤2.2)迭代变量X，V₁，V₂，D₁，D₂，交替迭代过程如下：

其中diag(·)为对角矩阵函数，I为单位矩阵，U∑V^T是矩阵的奇异值分解，函数(x)₊＝max{x,0}，max表示取最大值操作；μ为正则化参数；

步骤2.3)剔除非真正端元，并更新光谱库A：

判断光谱库中每个物质端元A_i，A_i表示光谱库A中第i个端元，所对应的丰度权重系数gamma_i，i＝1,2,...,m，是否小于阈值ω，阈值ω大小根据实际高光谱数据设置，若gamma_i＜ω则认为物质端元A_i为非真正端元，剔除光谱库中A所对应的端元A_i，否则保留光谱库中A所对应的端元A_i，将剔除非真正端元后的光谱库A更新；

端元A_i所对应的丰度权重系数gamma_i计算如下：

步骤2.4)停止迭代条件1：判断丰度矩阵的收敛条件是否满足，若是，则执行步骤3，否则执行步骤2.5)；

丰度矩阵的收敛条件为：

取ε＝10^-6；

步骤2.5)停止迭代条件2：将求解稀疏解混模型的迭代次数k加1，判断迭代次数k是否小于最大迭代次数，若是，则执行步骤3，否则，执行步骤2.2)。