基于低秩重建与TV正则的高光谱稀疏解混

徐盈盈，黎建华

（台州学院 电子与信息工程学院，浙江 临海 317000）

0 引言

高光谱影像（hyperspectral imagery）是由卫星等遥感设备上搭载的成像光谱仪获取的三维数据，其中二维图像描述了地表空间特征，而光谱波段揭示了图像每一像元的光谱曲线，这些光谱涵盖了紫外、可见光、近红外和中红外等波段［1］。高光谱影像具有很高的光谱分辨率（通常波段宽度＜10 nm），其丰富的光谱信息可充分反映不同地物的特征，在作物估产、矿物勘探、环境监测、地形分类和军事侦察等领域有着广泛的应用［2］。由于成像传感器通常空间分辨率较低，数据采集过程受大气、光照等干扰，且地物之间存在自然混合等因素，高光谱影像中普遍存在混合像元（mixed pixel），即一个成像单元的光谱是多种物质光谱的混叠信号。混合像元分解，亦称解混（unmixing），其目的是从混合光谱信号中分离纯净端元（endmember）光谱并估计各端元所占丰度（abundance）比例［3］。混合像元分解可以使高光谱影像的应用从像元级提升至亚像元级［4］，对提高遥感应用的精准度、推动定量遥感技术的深入发展具有重要的科学意义。

基于先验光谱库的半监督线性解混方法是近年来的研究热点，利用事先获取的过完备端元光谱库，解混过程等价于从光谱库中挑选端元并求其丰度［5］。由于光谱库所含光谱数目通常远大于待解混的端元个数，该解混方法一般建模为稀疏诱导约束的回归模型，并采用稀疏回归算法进行求解，因此通常称此类方法为稀疏解混。目前稀疏解混的研究重点都在于挖掘待求丰度的先验信息并设计相应的丰度约束，如分别基于对丰度矩阵空间同质性和行稀疏性的考虑，Iordache等人先后提出了经典的SUnSAL-TV模型［6］与 CLSUnSAL 模型［7］。最近几年，Wang 等人［8］提出了 CCSU 模型，通过在协同稀疏解混框架中引入集中约束来更精确地描述丰度图的非局部冗余性。Huang等人［9］提出JSpBLRU模型，同时将影像块结构的联合稀疏与低秩作为丰度估计约束。Li等人［10］提出了一种局部光谱相似性保留约束，以维持局部区域光谱相似性的方式来提高稀疏解混的鲁棒性。Wang等人［11］提出DRSU-TV模型，通过使用双权值来增强光谱和空间域的丰度稀疏性。Zhang等人［12］提出光谱-空间同时加权的稀疏解混模型，更好地保留了丰度图的细节信息。

本文利用稀疏解混模型中重建数据的低秩特性，在解混过程中通过约束重构影像的数据特征来减少噪声等误差对丰度估计的干扰，提出了重建约束与丰度约束协同的稀疏解混模型。与经典解混模型的对比实验表明，稀疏解混模型能有效提高解混精度。

1 稀疏解混建模

现有的稀疏解混方法充分利用了待求丰度的结构特性，但没有从影像数据的角度出发设计模型约束。本文发掘线性解混问题中理想的重建数据应当具有的低秩特性，以此为先验，构造了重建影像的低秩约束项，并考虑丰度矩阵的空间平滑特性，提出结合重建影像低秩约束和丰度TV约束的稀疏解混模型。

1.1 线性解混模型

线性混合假设指像元光谱是若干端元按照各自的丰度比例线性叠加得到的混合信号，该假设符合大部分应用场景且易于建模，其矩阵形式可以表示为：

其中H∈ℝb×n表示高光谱影像数据的矩阵形式，b是光谱的波段数，n是影像包含的像元个数；W∈ℝb×m表示包含了m条纯净物质光谱的端元矩阵或光谱库；X∈ℝm×n是丰度矩阵，N∈ℝb×n为表示模型误差与噪声的扰动项。

在基于先验光谱库的解混问题中，给定的过完备光谱库W一般包含成百上千种物质的光谱，而实测影像H内包含的e种端元通常只有几个到几十个，远远小于光谱库的列维数，即e≪m。因此，与冗余光谱库对应的丰度矩阵X理论上仅有e行非零元素，即X具有稀疏性。此外，鉴于丰度的物理属性，X是满足每列元素之和为1的非负矩阵。基于线性混合模型（1），光谱库先验的解混通常建模为如下的稀疏回归问题：

其中‖X‖0表示用l0范数统计X中的零元素个数，‖·‖F表示用于度量回归误差的Frobenius范数，ϵ表示误差容忍值。式（2）可以用正交匹配追踪等贪婪算法直接求解［13］，或者将l0范数近似替换为更易求解的平滑连续函数或lp范数（0＜p≤1）后再优化求解［14］。

1.2 基于低秩重建和TV正则的解混模型

不同于仅含丰度约束的稀疏解混模型，文中根据重建影像的数据特点构造新的约束项，在SUnSALTV的基础上提出了更鲁棒的扩展模型。Zhao等人［15］已经论证了在无噪声的理想情况下，线性混合的高光谱影像是低秩数据，而实际观测影像由于噪声等其他干扰往往是满秩数据。模型中重建影像WX，本质上正是由光谱库W中的很少量光谱列以丰度X为权重线性混叠而成，因此理想的WX是具有低秩特性的矩阵，且WX的秩通常等于激活端元的个数。据此，将稀疏解混问题构造为如下优化方程：

其中rank(·)表示矩阵的秩，TV(·)表示矩阵的全变分，α,β,γ均为非负参数。

显然式（3）是NP-hard的组合优化问题，为了便于极小化求解，我们对目标约束其进行凸松弛近似，并将模型（3）改写为如下无约束的凸优化问题：

其中‖·‖*表示的核范数是矩阵秩的最小凸包络函数，定义为‖WX‖*=∑iσi(WX)，σi(WX)表示WX矩阵经降序排列的第i个奇异值；‖·‖1是矩阵l0范数的最紧凸松弛l1范数；τ+(X)表示丰度矩阵的非负约束，即X的元素值均为非负时有τ+(X)=0，否则τ+(X)=+∞。PX≡[PhX;PvX]是丰度矩阵的全变分正则的线性算子表示，其中PhX和PvX分别表示丰度图的横向和纵向差分。

模型（4）的核心思想就是通过增加对重建影像的结构约束实现更合理的丰度估计。从理论上分析，目标方程F(X)中对待估丰度矩阵X的全变分约束促使各端元对应的丰度图在相邻像元之间为平滑过渡，符合自然物质分布的局部平滑特性；对重建影像WX的低秩约束，使得该解混模型在数据重建时对噪声有更好的鲁棒性，进一步引导丰度估计朝着有利于无噪声影像重建的方向进行优化。

2 基于ADMM的数值优化算法

优化方程（4）可以采用交替方向乘子法（Alternating Direction Method of Multipliers,ADMM）［16］进行求解。ADMM算法框架的基本思想是通过变量替换的方法，把难以整体优化的原问题分裂为多个易于求解的子问题。我们引入V1～V5共五个辅助变量，将式（4）等价改写为如下带约束的优化问题：

可以将式（5）写成如下的增广拉格朗日形式：

其中参数ρ是大于0的惩罚系数，D/ρ表示与约束GX+BV=0对应的拉格朗日乘子。文献［17］论证了对于形如式（6）的优化问题，在矩阵G满足列满秩且方程g是闭凸函数的条件下，对任意ρ＞0，通过ADMM迭代极小化式（6）可以收敛到存在的解。

为便于迭代公式推导，我们将式（6）所示的增广拉格朗日方程展开如下：

其中D1～D5分别是对应于式（5）中五个等式约束的拉格朗日乘子。

根据ADMM方法，我们将对式（7）中的每个变量分别优化，即在极小化某一个变量相关的方程时固定其余变量的值，依次交替迭代至达到算法的终止条件。令t表示迭代次数，t=0时所有待优化的变量均初始化为全零矩阵。以t+1轮迭代为例，本文的模型优化被分为以下几个子问题的求解：

1）式（7）对丰度矩阵X求导并令导数等于0，可得X的迭代优化公式为：

其中I表示单位矩阵。

2）重建影像的替换变量V1的优化目标是核范数约束的极小化问题，我们可以套用奇异值阈值（singular value thresholding，SVT）算法［18］，直接给出V1的求解公式：

记矩阵Z的奇异值分解为Z=UΣS*，则根据文献［18］可以定义：

其中矩阵U和S的列分别是矩阵Z经过奇异值分解后的左奇异向量和右奇异向量，S*表示酉矩阵S的共轭转置；矩阵Σ的对角元素σi是矩阵Z的奇异值。

3）通过求导可得用于空间约束的丰度矩阵替换变量V2的更新公式：

其中卷积算子P可以利用离散傅里叶变换进行快速计算。

4）丰度矩阵全变分约束的替换变量V3的目标方程可以用软阈值（soft thresholding）算法［19］进行求解，其更新公式如下：

其中soft(·,η)表示以η为参数的软阈值函数。

5）丰度稀疏约束的替换变量V4的迭代公式同样可用软阈值方法直接得到：

6）丰度非负约束的替换变量V4按如下公式将结果投影至非负象限：

7）拉格朗日乘子D1～D5的更新公式如下：

算法的终止条件是t达到最大迭代次数tmax或者相邻两次迭代的相对误差满在下文的实验中，我们设定tmax=200,ξ=10-4。

3 实验结果与分析

在模拟数据和真实数据上分别测试文中模型与算法的解混效果。通过对比稀疏解混算法与经典解混算法SUnSAL-TV和CLSUnSAL的实验结果，验证了结合重建数据低秩约束与丰度空间约束模型的有效性。实验采用固定其他参数，单独调节其中一个参数的方式交替调参，所有实验方法的模型参数均选自统一参数集合Γ={ }0,0.0005,0.001,0.005,0.01,0.05,0.1,0.5,1,1.5,3,5,10,20，最终取使各算法达到最佳效果的参数组合，以实现公平比较。代码的实现基于MATLAB R2018b平台，硬件配置为英特尔i7-8750H处理器（2.20GHz）和16GB内存。

3.1 模拟数据实验

模拟实验采用的光谱库来自USGS公开的splib06a数字光谱库［20］，该库包含了在波长0.4～2.5微米范围内采集的240种纯净物质的光谱信号，每条光谱记录了224个波段的电磁波信号值。利用开源的模拟丰度图生成代码［21］，可以得到不同混合程度且满足丰度值非负与丰度矩阵列和为1等物理属性的丰度数据集。实验以尺寸为64×64个像素且包含5个端元的模拟丰度为例，随机生成的丰度数据如图1所示，其中每个像素点的数值范围为［0，1］。根据线性混合模型，先从光谱库中随机挑选若干端元光谱，再以随机匹配的丰度图作为权重比例线性叠加，即可合成三维的高光谱模拟混合数据。

图1 模拟混合数据中5个端元对应的真实丰度图

由于模拟数据实验已知真实的丰度图，那么可以通过比较解混算法得到的估计丰度与真实丰度的差距，设计定量指标，给出对解混效果的客观评判。信号与重建比误差（signal-to-reconstruction error，SRE）作为最常用的定量指标，以分贝（dB）为度量单位的SRE定义如下：

其中Xtrue表示真实丰度矩阵，X是解混算法估计的丰度矩阵，E[·]表示数学期望。SRE指标的数值越大，表示解混的丰度估计越精确。

为了模拟不同信噪比的高光谱影像，我们对混合数据添加了不同强度的高斯噪声，得到信噪比分别为SNR=20，30，40dB的三组模拟数据。通过实验参数调优，分别为这三组数据设定模型参数三元组(α,β,γ)为(5,0.05,0.005),(1.5,0.01,0.005),(0.5,0.001,0.001)。三种实验方法在三组模拟数据上得到的解混定量指标结果，如图2所示，可以看到本文方法对不同信噪比数据的解混均取得了高于对比方法的SRE值。对于SNR=20dB这样的大噪声数据，本文方法比SUnSAL-TV算法的结果高出了约3.5dB的SRE值，充分说明了引入低秩重建约束的有效性。在模拟数据的解混实验中，定量指标SRE值的高低已能够客观反映解混结果的好差。如图3所示，以SNR=30dB的模拟数据的解混结果为例，直观展示了文中的丰度估计结果与图1所示真实丰度的相似性。

图2 本文方法与对比算法对不同信噪比模拟数据解混得到的SRE（dB）值

图3 本文算法对模拟数据解混得到的估计丰度图

3.2 真实数据实验

真实数据实验采用公开的常用解混数据集：城市影像（Urban data）及其光谱库［22］。实测城市影像数据的第80个光谱波段的场景，如图4所示，该场景内主要包含四种物质，分别为屋顶、草地、树木和沥青路面。该影像数据的空间分辨率约为2 m×2 m，含有307×307个像元，每个像元记录了162个有效的光谱波段值，实验用的先验光谱库是包含上述四种物质端元光谱的204条光谱集合。

图4 实测高光谱数据

由于真实影像数据中各物质的丰度未知，实验结果可以根据丰度图的视觉效果结合辅助定量指标进行评判。我们用影像的重建误差（reconstruction error,RE）和丰度的稀疏度（sparsity of abundance,SA）作为评价指标，定义：

其中影像重建误差RE能反映用估计丰度和对应端元光谱重建的混合数据与输入影像的误差；稀疏度SA的本意是计算丰度矩阵中0元素的占比，但理论0值在实际优化中通常对应于极小的非零数值解，因此我们用NUM(X＜θ)统计矩阵X中数值小于θ的元素个数，实验中取θ=0.001。好的解混结果应当能同时满足较低的重建误差与较高的丰度稀疏度。由于真实数据实验的两个辅助指标相互影响，那么以RE为首要调参指标，即优先选择使RE值更小的参数组合，在有多组参数使得RE值相同的情况下，再选择使SA值更大的为最优参数。以此为参数调优准则，实验中设定的模型参数(α,β,γ)为（0.01，0.005，0.01）。

三种实验方法在城市影像数据上的解混指标结果见表1所示。总体来看，本文方法取得了最低的重建误差值，丰度稀疏度则介于两个对比方法之间。CLSUnSAL方法的估计丰度虽然有最大的稀疏度，但其影像重建误差远远大于本文方法。如图5所示，从估计丰度的视觉效果上可以看到该城市影像中的四种主要物质的丰度分布，从左到右分别对应屋顶、草地、树木和路面，与图4所示的观测场景基本一致。结合辅助评价指标，可以说明本文方法对提升实测高光谱影像解混是有效性的。

表1 本文方法与对比算法对城市影像解混得到的重建误差（RE）与丰度稀疏度（SA）指标值

4 结语

在基于过完备光谱库的解混问题中，我们利用重建影像具有低秩特性的先验信息，引入重建数据的低秩约束来更好地引导丰度估计，再结合丰度矩阵的空间全变分约束构建了具有良好抗噪性能的稀疏解混模型，在数值优化上采用了基于ADMM的算法实现快速求解。通过与两个经典稀疏解混方法的对比，我们在高光谱模拟数据和实测影像上的解混结果均证实了本文方法能有效提高丰度估计的精度。