论文推荐 | ​张作宇,廖守亿,孙大为,张合新,王仕成：​稀疏差异先验信息支持的高光谱图像稀疏解混算法

《测绘学报》

构建与学术的桥梁        拉近与权威的距离

复制链接，关注《测绘学报》抖音！

【测绘学报的个人主页】长按复制此条消息，长按复制打开抖音查看TA的更多作品##7NsBSynuc88##[抖音口令]

本文内容来源于《测绘学报》2020年第8期，审图号GS（2020）4062号。

稀疏差异先验信息支持的高光谱图像稀疏解混算法

张作宇¹,廖守亿¹,孙大为²,张合新¹,王仕成¹                            

1.火箭军工程大学, 陕西 西安 710025;
2.火箭军士官学校, 山东 青州 262500

基金项目：国家自然科学基金（61673017；61403398）

摘要：基于光谱库的高光谱稀疏解混技术近年来得到了人们的关注，该技术利用光谱库中光谱样本作为端元，将解混问题转化为稀疏表示问题。然而，由于测量环境的差异，待解混图像的实际端元往往与光谱库中相应光谱信号存在差异。本文提出了一种光谱差异稀疏约束的联合稀疏回归解混算法。首先，假设光谱差异具有稀疏特性，建立了光谱库校正模型，使得在解混过程中可对光谱库进行自适应地调整；然后，将光谱库校正模型与联合稀疏回归解混模型结合，建立了考虑光谱差异的稀疏解混模型；最后，基于交替方向乘子法得到了迭代优化解决方案。分别利用仿真和真实高光谱数据进行了试验验证，结果表明，在光谱库不匹配的情形下，本文方法能够有效提高稀疏解混算法的解混性能。

关键词：高光谱图像    解混    稀疏回归    光谱差异    光谱库校正

引文格式：张作宇, 廖守亿, 孙大为, 等. 稀疏差异先验信息支持的高光谱图像稀疏解混算法. 测绘学报，2020，49(8)：1032-1041. DOI: 10.11947/j.AGCS.2020.20190205.

阅读全文：http://xb.sinomaps.com/article/2020/1001-1595/2020-8-1032.htm
全文概述高光谱成像技术能够在几十或数百个窄波段内对物体和场景进行成像，其波长范围可覆盖可见光和红外的各个波段^[1]。由于具有光谱分辨率高和覆盖波段范围广的优势，高光谱成像技术被广泛应用到各个领域中，如矿物勘探、军事侦察、农作物分析、食品安全等^[2-4]。然而，由于高光谱图像的空间分辨率通常较低，特别是遥感高光谱影像，单个像元内往往包含多种材质，这为实现地面材质的精准分析带来了困难。高光谱解混技术旨在获取图像中每个像元所包含的纯净物质及每种纯净物质在该像素中所占的比例，是对高光谱图像进行分析的重要手段。

按照像元混合方式的不同，光谱混合模型通常分为线性混合模型(linear mixing model, LMM)和非线性混合模型。其中，LMM假定图像中的每个像元都可以被端元线性表示^[5]。而非线性混合模型考虑了光谱混合过程中的非线性因素，如致密混合情况下光线通过微小颗粒产生的散射效应和层级混合情况下光线在物体间发生的多次反射。然而，由于非线性混合模型通常求解复杂且不易扩展，在大多数现实场景中LMM能够满足解混需求，许多高光谱解混工作都是基于LMM来开展^[6-7]。因此，本文同样采用了LMM进行解混。传统的线性解混算法大多采用无监督的方式，主要包括凸几何法^[8-9]、非负矩阵分解^[10-12]、独立成分分析^[13]、贝叶斯方法^[14]等。这些方法往往需要对端元和丰度都进行估算，因此容易受到图像中不存在纯净像元或所估算的端元不具有物理意义的制约。

近些年，人们花费了大量人力物力对不同材质的样本进行了光谱测量，形成了适用于各种用途的光谱库。例如，美国USGS光谱库^[15]包含1300多种光谱样本，覆盖了矿物、岩石、液体、人造物、植被等多种类型。这些丰富的光谱信息可有效地用于解混。假设场景中的所有端元都被光谱库覆盖，依据LMM，图像中每一个像元都可被光谱库中的一小部分光谱线性表示。将光谱库作为字典，利用光谱库进行解混的过程可看作一种稀疏表示问题。因此，许多稀疏表示的理论和工具可应用到解混中^[16-17]。文献[18]将高光谱解混看作单个像元的稀疏回归问题，提出了基于变量分离和增广拉格朗日的稀疏解混算法(sparse unmixing algorithm via variable splitting and augmented Lagrangian，SUnSAL)。为联合空间信息，文献[19]将全变差(total variation，TV)约束引入SUnSAL中，提出了SUnSAL-TV算法。由于稀疏解混中端元通常只占光谱库很小的一部分，且通常光谱库具有较强的自相关性，文献[20]利用丰度的联合稀疏特性提出了联合SUnSAL算法(collaborative SUnSAL, CLSUnSAL)。

当前，大多稀疏解混算法假设场景中的实际端元与光谱库中相关光谱保持一致。然而，由于成像尺度、光学系统误差、环境因素(如光照，温度)等方面的不同，图像中端元很可能与光谱库中相关光谱存在差异，这将影响算法的解混效果。最近，一些学者致力于研究光谱库与实际场景不匹配的问题，取得了一些有益的研究成果。文献[21]将光谱库校正变量引入稀疏表示项中并假设光谱差异服从高斯分布，提出了字典调整的非凸增强稀疏回归(dictionary-adjusted nonconvex sparsity-encouraging regression，DANSER)算法。文献[22]采用冗余光谱代表高光谱图像中的光谱差异成分，使冗余光谱不参与解混过程以减轻光谱差异带来的影响，提出了稀疏冗余解混(sparse redundant unmixing，SpaRedU)算法。DANSER引入了光谱库校正变量，有效地刻画了光谱变异过程，然而该算法假设光谱差异服从高斯分布，这可能与实际的光谱变异特性不符。SpaRedU剔除了冗余光谱项，使存在光谱差异情况下的丰度估计得到了一定改善，然而该方法并不能得到变异后的端元。

光谱变异可由与光谱库信号相关的光谱变异和与光谱库信号不相关的光谱变异两部分组成。其中，与光谱库信号相关的光谱变异成分通常由光照等环境因素产生，主要表现为信号整体的强度变化。这部分光谱变异对解混结果产生的影响可采用在解混过程中忽略丰度和为1的条件来缓解。而与光谱库信号不相关的光谱变异形成原因较为复杂。在地质过程中，矿物分子中的特定原子会在特定条件下被化学性质相近的其他原子替代，这种现象被称为原子替换。原子替换是产生不相关光谱变异的重要因素之一，在不同矿物之间发生较为普遍。由于原子替换只发生在特定的原子上，由此产生的矿物分子光谱特性的改变通常与替换后原子的光谱吸收特性相关。由于每个原子都具有特定的吸收波峰，因此变异前后的光谱主要表现为吸收波峰的变化，由此产生的光谱差异具有较强的波段属性。

本文重点考虑了产生光谱变异的原子替换过程，假设光谱库与实际端元只在特定波段范围内存在光谱差异，提出了一种光谱差异稀疏约束的联合稀疏回归(spectral difference sparse constrained collaborative sparse regression，SDSCSR)解混算法。引入了光谱库校正变量来描述光谱差异，利用了光谱差异的稀疏特性，并基于交替方向乘子法(alternating direction method of multipliers，ADMM)^[23]得到了迭代优化解决方案。由于原子替换现象主要存在于地表不同矿物之间，因此本文方法更适用于遥感中的矿物勘探应用。

1  联合稀疏解混模型联合稀疏解混模型在光谱库自相关性较强的情况下具有较好的解混性能，在许多研究中得到了应用^{[2, 6]}。本节对LMM和联合稀疏回归进行简要介绍。

1.1  LMM

LMM假设高光谱图像中任意像元的光谱响应可由端元线性组合而成。假定表示高光谱图像中具有L个波段的像素矢量，A=为具有M个端元的端元矩阵，LMM可表示为

(1)

式中，为y对应的丰度矢量；为误差项。在LMM中，丰度通常还要满足两个条件，即丰度非负约束(abundance nonnegativity constraint，ANC)和丰度和为1约束(abundance sum-to-one constraint，ASC)，可用式(2)表示

(2)

ANC和ASC使得LMM更具物理意义。然而，由于光照和环境的影响，真实高光谱图像中往往存在较大的光照强度变化，ASC约束饱受争议，因此，本文忽略ASC约束。假设包含N个像素的高光谱图像表示为，丰度矩阵表示为，误差矩阵表示为，则LMM可表示为如下矩阵形式

(3)

1.2  联合稀疏回归

在传统盲解混中，通常A和X都依据Y来获取。而当利用光谱库进行解混时，将光谱库作为字典，光谱库中每一种材质光谱作为一个原子，解混过程可以看成用字典中原子对高光谱图像进行线性表示的过程，参与表示的原子即为端元。由于参与表示的原子只占光谱库很小的一部分，上述问题是一个稀疏回归问题。定义包含K个光谱样本的光谱库为，字典辅助的光谱混合模型可表示为

(4)

式中，为行稀疏矩阵，C的非零行构成了丰度矩阵X。联合稀疏回归的核心思想是，充分利用C的行稀疏属性，提高端元识别的准确性和丰度估计的精度。联合稀疏回归解混方法的目标函数可表示为

(5)

式中，λ>0为常数，代表C的每个行向量的2范数之和，即。该约束可有效提升C的行稀疏度，使C中大部分行向量的值为零。由此可看出，联合稀疏回归的目标是寻找一个非负行稀疏的矩阵C，使得Y≈DC。

2  SDSCSR解混算法传统稀疏解混算法通常直接采用光谱库中的光谱信号作为端元。然而，现实中很难找到一个光谱库正好匹配当前图像，实际端元很可能与光谱库中对应光谱信号存在差异。在这一节中，为减轻光谱库不匹配对解混性能造成的影响，提出了SDSCSR解混算法。

2.1  稀疏差异先验信息支持的联合稀疏解混模型

在地质过程中，环境中的实际端元通常会由于原子替换而发生轻微变化。在特定条件下，分子中的一些离子会被其他离子替代。这种替换的发生通常是少量且缓慢的，因此并不影响材质的基本属性，变化前后的材质依然可以认为属于相同类型。然而，原子替换会造成分子震动频率和电子能级跃迁的变化，从而影响材质对光线吸收波峰的变化，使材质光谱发生变化。由于分子震动和电子能级跃迁均对特定频率的电磁波敏感，由此引起的光谱变化具有较强的波段属性。图 1展示了两种Alunite变体的光谱曲线及吸收波峰位置。对于Alunite，K⁺离子可以被替换为Na⁺，Alunite GDS84 Na03中的K⁺离子含量要远高于Alunite GDS82 Na82中的含量。从图 1可看出，Alunite GDS84 Na03在1.5、1.8、2.2和2.3 μm附近的光谱反射率明显低于Alunite GDS82 Na82。这主要由于原子替换造成了K⁺离子含量的不同。研究表明，原子替换是造成光谱变异的重要原因之一^[22]。由于原子替换主要在材质光谱的特定波段上产生影响，因此可认为由其产生的光谱差异是稀疏的。

图 1                         Alunite GDS84 Na03和Alunite GDS82 Na82光谱反射率图及吸收波峰 Fig. 1     Reflectance of Alunite GDS84 Na03 and Alunite GDS82 Na82 as well as their absorption peaks      

图选项

基于以上分析，本文建立了光谱库校正模型，假设图像中的任一实际端元a_n可由式(6)表示

(6)

式中，d_kn为端元a_n对应的光谱库光谱信号；r_kn为光谱库校正变量。

在模型中，光谱库并不是无限制地变化，而是在一定的范围内变化。实际端元依然需要被光谱库覆盖。为充分利用光谱差异的稀疏特性，对r_kn进行如下约束

(7)

式中，τ>0为常数，代表光谱变化的容忍度。

由于光谱库通常包含几十或上百种光谱样本，光谱样本间很可能具有较强的相关性。联合稀疏解混在光谱库自相关性较强的情况下相比于普通稀疏解混方法性能更好。本文以联合稀疏解混为基础，依据建立的光谱库校正模型，并对式(7)进行松弛化处理，提出了新的解混模型

(8)

式中，α、β > 0为约束项系数；D^′为校正后光谱库；R为光谱库校正矩阵。相比于原始联合稀疏解混模型，新模型主要添加了光谱库校正项以减轻由光谱库不匹配造成的不良影响。其主要思想是，允许光谱库中光谱信号在一定范围内进行局部调整，并充分利用光谱差异的稀疏特性，使得调整后光谱库更加适应当前场景。

2.2  基于ADMM的迭代优化解决方案

为求解式(8)的最小化问题，利用ADMM设计迭代优化算法。由于C、D^′和R都是未知的，该问题是非凸的。采用交替优化的求解策略，即一次只针对一个变量进行优化，保持其他变量不变。为方便起见，式(8)可被重写为如下形式

(9)

式中，ι_R+(C)为指示函数，当C中所有元素均为非负时，ι_R+(C)的值为0，否则为+∞。

首先，对C进行优化，对于第t+1次迭代，针对C的子目标函数为

(10)

可以看出，式(10)与联合稀疏回归目标函数相同，可采用类似的解决方案。式(10)等价于如下形式

(11)

为更加简洁，式(11)可被重写为

(12)

式(12)的增广拉格朗日函数为

(13)

式中，μ > 0；W=(W₁, W₂)；W/μ代表拉格朗日乘子。

分别针对C、V和W进行优化可得到式(10)的解决方案。对于C的更新，需要对式(14)进行求解

(14)

其解决方案为

(15)

对于V₁的更新，需要求解式(16)

(16)

可采用向量软阈值^[24]的方法来求解式(16)，独立地对整行变量进行更新

(17)

式中，ζ=C^t+1-W₁^t；vect-soft(b, θ)代表向量软阈值函数，采用式(18)计算

(18)

对于V₂的更新，需要求解式(19)

(19)

可通过将V₂的解投影到非负坐标轴上来考虑ι_R+项，得到V₂的更新公式

(20)

接下来对D^′进行优化，考虑式(21)的子目标函数

(21)

其解决方案为

(22)

最后，对R进行更新，考虑式(23)子目标函数

(23)

可采用软阈值法^[25]对式(23)进行求解

(24)

式中，ψ=D^′t+1－D；soft(b, θ)代表软阈值函数，采用下式计算

(25)

式中，.*和./分别代表矩阵对应元素的乘和除操作。

算法的整体流程可总结如下。

输入：高光谱数据Y，光谱库D，约束项系数λ、α、β。

输出：校正后光谱库D^′，丰度矩阵C。

(1) 初始化迭代次数t=0，初始化D^′0，C⁰，R⁰，V⁰，和W⁰。

(2) 利用式(15)更新C^t+1。

(3) 利用式(17)更新V₁^t+1。

(4) 利用式(20)更新V₂^t+1。

(5) 利用式(22)更新D^′t+1。

(6) 利用式(24)更新R^t+1。

(7) 更新拉格朗日乘子。

(8) W₁^t+1=W₁^t-C^t+1+V₁^t+1。

(9) W₂^t+1=W₂^t-C^t+1+V₂^t+1。

(10) 更新迭代次数：t=t+1。

(11) 返回步骤(2)直到满足停止条件。

在以上算法步骤中，采用最小二乘法对C进行初始化，将D^′初始化为D，将R初始化为零矩阵。为使算法尽可能充分收敛又不至于运行时间过长，迭代停止条件被设置为两个条件，一是达到最大迭代次数，二是满足，其中σ > 0为误差容忍度，以上任意一个条件达到，则停止迭代。在算法中，参数μ对收敛速度有较大影响，本文采用在迭代中不断地更新μ值，使得初次残差和二次残差的比例保持在一定范围内。考虑到本文算法的计算复杂度，在算法中最耗费时间的步骤是C和D^′的更新，其中计算复杂度最高的运算分别为D^′TY和YC^T，它们的计算复杂度均为。算法中其他步骤均具有更低的计算复杂度。因此，本文算法的总体计算复杂度是。

3  试验结果与讨论为验证本文算法的有效性，分别采用模拟数据和真实数据进行了试验。在这些试验中，将本文算法与其他3种解混算法进行对比，分别是SUnSAL、CLSUnSAL和DANSER。采用信号与重构误差比(signal-to-reconstruction error，SRE)来评估算法的性能，具体定义为

(26)

式中，X为真实丰度矩阵；C为行稀疏的真实丰度矩阵；为解混算法输出的估计丰度矩阵。SRE值越大，丰度估计结果就越准确。

3.1  模拟数据试验

在模拟数据试验中，从USGS光谱库中随机选取40种物质的光谱信号作为光谱库，这些光谱信号包含224个波段并均匀地分布在0.4~2.5 μm。选取两个包含相同类别材质光谱的光谱库用于试验，分别为和，它们在相同材质光谱之间存在一定的光谱差异。采用球形高斯场随机生成大小为100×100像素的仿真丰度图像，如图 2所示。随机从D₁和D₂中选取M个不同种类物质的光谱作为端元，其中从D₂中选取的端元个数为J，依据端元和丰度生成高光谱数据。为模拟实际高光谱图像中的噪声，将具有一定信噪比(signal-to-noise ratio，SNR)的高斯白噪声添加到仿真数据中。在解混过程中，采用D₁作为光谱库进行解混，将来自D₂中的端元视为变异端元，从而模拟光谱库不匹配现象。在接下来的试验中，除非特殊说明，设置端元总数M为8，不匹配端元个数与端元总数之比J/M为0.5，SNR为30 dB。

图 2                         6个端元对应的仿真丰度图像 Fig. 2     The simulated abundance maps of six endmembers      

图选项

首先，讨论参数设置对算法性能的影响。本文算法主要包括3个参数：λ、α和β。对于λ，由于λ控制联合稀疏约束的强弱，本文参考CLSUnSAL对λ进行设置。α控制光谱库校正项的强弱，类似于μ、α对算法的收敛速度和稳定性有影响，本文将α设置为固定值1000。图 3展示了不同β值对算法性能的影响。可以看出，算法在β∈[3, 100]区间取得了较好的解混结果，过小的β值会使R的变化不受约束从而造成算法性能的大幅下降，而过大的β值可能造成无法对光谱库进行有效调整。

图 3 不同β值对应的试验结果 Fig. 3     The experimental results for different β values      

图选项

为充分验证算法的解混性能，生成了不同M、J和SNR值情况下的模拟数据。在接下来的仿真中，无论仿真数据如何变化，解混算法的参数保持不变，对比算法采用原文献中推荐的参数设置。为SDSCSR、SUnSAL、CLSUnSAL、DANSER的稀疏约束项系数λ分别设置为0.1、0.1、0.1、0.5。此外，为DANSER设置参数为p=0.75和μ=10⁵，为SDSCSR设置参数为(α, β)=(1000, 10)，将最大迭代次数设置为1000次。为提高试验结果的可靠性，每组试验重复进行20次，取其平均值作为最终结果。

图 4展示了算法在不同端元个数下的SRE结果。在这个试验中，M的取值为4~12均匀取值，并保持J/M的值为0.5。可以看出，在所有的M取值中，SDSCSR均取得了最好的解混结果，只有当M为12时，SDSCSR的优势较小。本试验验证了SDSCSR对端元个数的抗差性。

图 4 不同方法在不同端元个数下对应的试验结果 Fig. 4     The experimental results for different methods under different endmember numbers      

图选项

图 5展示了算法在不同不匹配端元个数下的解混性能。在这个试验中，保持端元总数为8个，通过变化J的取值控制光谱库不匹配的程度。可以看出，在存在不匹配端元的情形下，SDSCSR的解混结果均优于其他对比算法。在不匹配端元个数为0的情形下，SDSCSR取得了与CLSUnSAL相当的解混结果。特别地，随着J取值的增大，SDSCSR相对于CLSUnSAL的优势更加显著。

图 5 不同方法在不同不匹配端元个数下对应的试验结果 Fig. 5     The experimental results for different methods under different mismatched endmember numbers      

图选项

图 6展示了不同噪声水平下算法的解混性能。在这个试验中，仿真图像的SNR被设置为在15~40 dB均匀变化。可以看出，当SNR高于25 dB时，SDSCSR比其他算法有更好的解混结果。在低SNR的情形，SDSCSR和DANSER的性能下降明显，这主要由于较高的噪声水平会对它们的光谱库校正模型产生较大影响，使得噪声成分被误认为是光谱变异。

图 6 不同方法在不同噪声水平下对应的试验结果 Fig. 6     The experimental results for different methods under different noise levels      

图选项

表 1统计了各个算法的运行时间。SDSCSR比CLSUnSAL花费更多的时间收敛，这主要是由于SDSCSR的计算复杂度更高。

表 1 不同算法的运行时间Tab. 1 Time cost of different algorithmss

algorithm	SUnSAL	CLSUnSAL	DANSER	SDSCSR
time	9.557	9.958	70.251	26.029

表选项

3.2  真实数据试验

采用Cuprite数据对SDSCSR及其对比算法进行测试。该数据由AVIRIS高光谱成像仪于加拿大Nevada矿区上空拍摄，具有224个波段均匀地分布在0.4~2.5 μm，采集时间为1997年。Cuprite数据被广泛地用于高光谱研究，场景中一些材质对应的丰度图像已得到大家的普遍认可。为减少噪声的影响，去除数据中的低SNR波段(1—2，104—113，148—167，221—224)，利用剩余的188个波段进行试验。由于原始Cuprite数据较大，试验截取了包含250×191像素的子图像，其伪彩色合成图如图 7(a)所示。试验中采用的光谱库包含来自于USGS光谱库的240个光谱。为便于试验，光谱库中对应的低SNR波段也被去除。研究指出，Cuprite数据中的实际光谱与USGS光谱库中对应光谱存在一定的校正误差^[21]。图 7(b)展示了USGS于1995年制作的该地区矿物分布图。需要注意的是，公开的AVIRIS Cuprite数据采集于1997年，而矿物分布图制作于1995年，因此他们之间可能存在一定的误差。然而，由于真实的丰度图像无法获得，该分布图可以作为参考从而对算法估计的丰度图像进行定性分析。

图 7                         Cuprite数据 Fig. 7     The Cuprite dataset      

图选项

试验中，对比算法依然采用原文献中推荐的参数设置，为SUnSAL和CLSUnSAL分别设置稀疏约束参数λ为0.001和0.01，为DANSER设置参数(p, λ, μ)=(0.5, 3, 1000)，为SDSCSR设置参数(λ, α, β)=(0.01, 1000, 10)。图 8展示了不同算法解混得到的4种矿物对应的丰度图像以及由Tricorder 3.3软件^[26]得到的矿物分类图。这4种矿物被认为在Cuprite数据中占主要成分。从图 8可以看出，SDSCSR得到的丰度估计结果与分类图具有较高的一致性。然而，SDSCSR估计的丰度图像与对比算法存在一些不同。对于Buddingtonite和Chalcedony，SDSCSR估计的丰度显著区域在视觉上更加接近于分类图，而对于Montmorillonite，SDSCSR估计的丰度显著区域与DANSER保持一致。因此，本文算法能够使光谱库较好地适应实际场景，在存在光谱差异的情况下具有较好的解混性能。

图 8                         Cuprite数据矿物分类图及不同解混解混算法得到的矿物丰度图 Fig. 8     Classification maps of the Cuprite dataset and abundance maps estimated by different unmixing algorithms      

图选项

4  结论

本文提出了一种稀疏差异先验信息支持的联合稀疏解混模型，该模型假设光谱差异具有稀疏特性，更加符合实际存在原子替换的情况。仿真和真实数据试验表明，本文算法在端元与光谱库不匹配情况下较对比算法能够提供更有效的解混结果。考虑到试验中本文算法在低SNR下的解混性能下降，未来将针对光谱差异特性研究更加细致的光谱库校正模型，提高算法对噪声的抗差性。此外，如何自动地选择算法参数也是将来需要研究的重要问题。

作者简介

第一作者简介：张作宇(1992-), 男, 博士生, 研究方向为高光谱图像混合像元分解与融合。E-mail:[e-mail]

权威 | 专业 | 学术 | 前沿

微信、抖音小视频投稿邮箱 | [e-mail]

微信公众号中搜索「测绘学报」，关注我们，长按上图二维码，关注学术前沿动态。

欢迎加入《测绘学报》作者QQ群： 751717395

进群请备注：姓名+单位+稿件编号