基于双稀疏域联合求解的高精度光谱恢复算法

刘世界; 李春来; 徐睿; 唐国良; 吴兵; 徐艳; 王建宇; LIU Shi-Jie; LI Chun-Lai; XU Rui; TANG Guo-Liang; WU Bing; XU Yan; WANG Jian-Yu

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

基于双稀疏域联合求解的高精度光谱恢复算法 PDF

- ORCID：
刘世界 ^1,2
- ORCID：
李春来 ¹
- ORCID：
徐睿 ¹
- ORCID：
唐国良 ^1,2
- ORCID：
吴兵 ^1,2
- ORCID：
徐艳 ^1,2,4
- ORCID：
王建宇 ^1,2,3

1. 中国科学院上海技术物理研究所空间主动光电技术重点实验室，上海　200083； 2. 中国科学院大学，北京　100049； 3. 国科大杭州高等研究院，浙江杭州　310024； 4. 上海科技大学信息科学与技术学院，上海　200020

中图分类号： TP751

最近更新：2021-05-12

DOI：10.11972/j.issn.1001-9014.2021.03.020

摘要

基于压缩感知的光谱成像系统需要合适的算法解码采样数据才能得到最终的光谱成像数据，传统单稀疏域变换算法会带来光谱细节损失等问题。针对该问题，本文提出了利用双稀疏域联合求解的方法（JDSD），将信号分解为低频部分和高频部分，并针对不同频率信号特点分别进行稀疏恢复，进而解码求解以实现高精度恢复信号。在数据验证中，首先利用OMP算法在频域内对光谱信息轮廓进行恢复，利用IRLS算法在空间域内对光谱细节进行补偿，分析了不同稀疏变换对于参数设置的影响，测试了不同算法组合的JDSD对于测试数据的恢复结果。对于500种光谱数据仿真测试表明，双稀疏域联合求解可将光谱恢复保真度大大提升，20%采样率情况下，SAM和GSAM指标由传统方法的0.625和0.515分别提升为0.817和0.659，80%采样率情况下，SAM和GSAM指标由传统方法的0.863和0.808分别提升为0.940和0.897。JDSD算法可以使得光谱吸收峰等细节特征得到高精度保持，对于基于光谱的特征分析、物质识别等应用具有十分重要的意义。

关键词

光谱成像; 光谱特征恢复; 计算成像; 压缩感知

引言

压缩感知是一种新的信息处理理论，指出可以同时对信息压缩和采样^［

1-2］，为数据的存储和传输降低了压力，具有很好应用前景。压缩感知在光谱成像领域内具有广泛的研究，其中一个具有代表性的成像应用是由Rice大学提出的单像素相机SPI（Single Pixel Imager）^{［参考文献 3

百度学术}3］，随后CHSI（Compressive Sensing Hyper Spectral Imager）系统在SPI的基础上增加了分光部件，实现了光谱成像^{［参考文献 4

百度学术}4］，CHSISS（Compressive Hyper Spectral Imaging by projections in both the Spatial and the Spectral domains）又在此基础上提出对空间信息和光谱信息同时进行编码的光谱成像系统^{［参考文献 5

百度学术}5］。Duke大学DISP（Duke Information Spaces Project）研究小组于2007年提出了CASSI（Coded Aperture Snapshot Spectral Imaging）系统，包含单色散结构的SD-CASSI（Single Disperser-CASSI）^{［参考文献 6-7}6-7］和双色散结构的DD-CASSI（Dual Disperser-CASSI）^{［参考文献 8

百度学术}8］。针对CASSI系统的编码问题，Henry Arguello等人提出可选波长探测的编码优化，大大降低了对不必要信息的采样，同时提高了重构精度^{［参考文献 9

百度学术}9］。压缩感知的应用是以信号本身或在变换域内具有稀疏性特点为前提的，包括光谱成像数据在内的绝大多数自然信号，在常用的正交变换域内都具有稀疏性的特点，如傅里叶变换域，小波变换域等，具体表现为在很窄的低频域带宽内，集中了信号的绝大多数信息，而在高频段具有很少甚至可以忽略的信息，这种特性在空间域则体现为数据相关性强、邻域平滑的特点。通常信号中的噪声扰动带来的信息量主要集中在高频段，压缩感知通常是通过感知并保留低频信息，去除高频段信息，因此可以达到去噪效果。但在遥感领域中，针对目标地物的光谱数据，研究者们更关注的是吸收峰等具有指纹特性的细节，这些细节使得光谱曲线变得不平滑，意味着信号的高频信息就变得丰富而不应该被忽略，而传统压缩感知恢复算法的处理方式难以准确地保持该类关键信息，大大降低了光谱应用的可靠性。一般地，压缩感知重构算法仅仅考虑将信号进行一个变换域的变换^{［参考文献 10-13}10-13］，而在单个变换域的条件下，难以将各个频段信息都完整恢复。实际上，一个含有吸收峰的光谱信息，低频段仍占据绝大部分，而吸收峰这些高频信息则在空域更加明显，因此本文提出双稀疏域联合求解（Joint of Double Sparse Domains ，JDSD），将信号首先在频域进行稀疏分解，重构恢复得到光谱曲线轮廓信息，并认为高频信息存在于除去低频采样后的测量残差之中，高频信息在空域更具有稀疏性，因此对该残差在空间域中恢复，得到结果对轮廓信息进行细节补偿，最终得到精确恢复的光谱曲线。实验结果表明，该方法大大提高了光谱恢复精度，为光谱应用打下了更好的基础。

1 算法原理

1.1　基本算法描述

从矩阵计算的角度看，压缩感知的过程可以看作是用一个 $m \times n (m ≪ n)$ 的测量矩阵 $Φ$ 对信号 $x$ 进行线性投影，获取感知测量值 $y$ ： $y = Φ x$ ，其中 $y$ 是m维向量， $x$ 为n维向量，该过程使目标信息从n维降为m维，即感知采样率为 $\frac{m}{n}$ 的少量信息，实现压缩并感知，在测量矩阵 $Φ$ 满足RIP（有限等距性质），并且 $x$ 满足一定稀疏性的前提下，可以准确恢复 $x$ ［1-2］。通常地，需要对信号 $x$ 进行稀疏分解，即 $x = Ψ θ$ ，其中 $Ψ$ 称作稀疏基，可以很多不同的形式，如DCT（Discrete Cosine Transform）正交基，小波变换基，冗余字典等，此时有以下关系：

y = Φ x = Φ Ψ θ = Θ θ

，

（1）

其中 $Θ = Φ Ψ$ 称为压缩感知矩阵， $θ$ 称为 $x$ 的稀疏系数。

压缩感知方程（1）是一个欠定方程，但信号的稀疏性约束可以保证目标解（即稀疏解）具有唯一性。针对该方程的求解具有很多算法，包括贪婪算法，凸优化算法等，但这些算法都是基于对目标信号进行一次稀疏性变换的前提下进行，即当信号稀疏基确定后，认为 $y$ 是 $Φ$ 对信号 $x$ 的投影的同时，也是 $Θ$ 对 $θ$ 的投影，进而对 $θ$ 进行求解。正如前面所说， $θ$ 中的较大的系数集中在信号的低频段，而高频段系数往往会在求稀疏解时因为数值不够大而被丢掉，因此导致类似于光谱数据吸收峰等信息的弱化甚至丢失。

实际上，含有吸收峰细节的光谱曲线看作平滑信号（轮廓）和吸收细节（尖峰）之和，图1给出了ENVI5.1光谱库中巴西棕榈蜡（Carnauba Wax）反射率光谱曲线分解结果。

图1 光谱数据信号可以分为轮廓和细节两部分

Fig.1 Spectrum data can be divided into contour and details

从数学的角度上，信号的压缩感知方程可以分解为：

y = Φ x = Φ (x_{l} + x_{h}) = Φ x_{l} + Φ x_{h} = y_{l} + y_{h}

，

（2）

其中：

y_{l} = Φ x_{l} = Φ Ψ θ_{l} = Θ θ_{l}

，

（3）

y_{h} = Φ x_{h}

（4）

式（2）将原信号 $x$ 分解为 $x_{l}$ 和 $x_{h}$ ，分别对应轮廓（低频）和细节（高频）部分， $y_{l}$ 和 $y_{h}$ 分别为低频采样和高频采样，此时就可以针对信号的高频和低频分别进行稀疏变换并恢复。通常地，轮廓部分（ $x_{l}$ ）空间域比较平滑，在频域 $Ψ$ 中稀疏性极强，因此可以在频域内恢复，即对式（3）求解 $θ_{l}$ ，而尖峰细节（ $x_{h}$ ）在空间域本身就已经表现出良好的稀疏性，不必在进行变换域转换，可以对式（4）直接求解 $x_{h}$ 。

在没有先验条件的前提下，对采样数据 $y$ 进行合适的高低频分解是比较困难的，实际求解中，可以先求解式（1），即对整体信号进行频域稀疏求解，不同的是，此过程中需要合理设置稀疏条件，以对 $y$ 保留一定量的采样残差向量 $y_{r}$ ，其中记：

y_{r} = y - \hat{y}

（5）

式（5）中， $\hat{y} = Φ \hat{x}$ 为进行完全求解后对应的采样值。利用JDSD求解时，令 $y_{l} = \hat{y}$ ，此时可以将 $y_{r}$ 认为是 $y_{h}$ ，进而完成对式（4）的求解。因此从另一个角度讲，JDSD算法也可以认为是对残差向量所包含信息的进一步挖掘，实现对恢复信息的补偿。

1.2　稀疏约束参数的设置

首先这里引入残差值函数 $R (β, Ρ, α)$ ：

R (β, Ρ, α) = \frac{1}{2} {‖β - Ρ α‖}_{2}^{2}

（6）

从式（6）中可以看出， $R (β, Ρ, α)$ 度量的是向量 $β$ 与 $Ρ α$ 的差距， $Ρ α$ 表示矩阵 $Ρ$ 对向量 $α$ 的线性投影值。

压缩感知的重构算法的约束条件主要有残差值 $R (y_{r}, Φ, x_{r})$ 和信号的稀疏度 $S_{L} (x)$ ，两者具有一定的等价性，因此可以设置其一。 $R (y_{r}, Φ, x_{r})$ 表示在重构过程中，恢复结果与采样数据的逼近程度，相比较而言，对于 $S_{L} (x)$ 的设置，更能直观体现压缩感知对于稀疏性依赖的特点，因此这里讨论对 $S_{L} (x)$ 的设置。

稀疏度 $S_{L} (x)$ 通常用信号的范数表征：

S_{L} (x) = {‖x‖}_{L}

，

（7）

$S_{0} (x)$ 表示 $x$ 的非零元素的个数，直接表征了信号的稀疏性，常见的算法还有1范数 $S_{1} (x)$ 优化，或者 $p (0 < p < 1)$ 范数 $S_{p} (x)$ 优化。在信号稀疏先验条件下，三者具有等价性。

目标方程要使得 $R (y_{r}, Φ, x_{r})$ 最小化，同时要以 $S (x)$ 大小作为约束，为了充分说明JDSD算法的本质，本文考虑以下两种方式：

m i n R (y, Φ, x) s . t s_{0} (x) \leq σ

（8）

m i n s_{p} (x), (0 < p < 1) s . t R (y, Φ, x) = 0

（9）

式（8）表示在保证信号的稀疏度前提下最小化 $R (y, A, x)$ ， $σ$ 表示信号的稀疏度上限，通常用贪婪算法求解，如MP（Matching Pursuit）和OMP（Orthogonal Matching Pursuit）算法等，其中OMP算法是对MP算法的改进^［

11］，该类算法可以很好地重构出信号的低频信息，因此本文利用OMP算法对式（3）进行求解恢复光谱曲线的轮廓信息。式（9）是子空间追踪（Basis Pursuit，BP）的变形^{［参考文献 12

百度学术}12］，若达到

R (y_{r}, Φ, x_{r}) = 0

的优化效果，可以充分提取采样数据y中所包含信息而得到理想的结果，但实际信号中，噪声会一直存在，此时残差不应该为0，因此一般求解方程为：

m i n s_{p} (x), (0 < p < 1) s . t R (y, Φ, x) \leq ζ

（10）

式（10）中，当信息被准确恢复时，ζ表征噪声大小。本文对式（10）用IRLS（Iterative Reweight Least Squares）算法求解，来恢复光谱曲线的吸收峰细节部分。

一般地，当利用单稀疏域变换求得式（8）的最优解时，σ值要设置比较大以保证足够的信息量可以被恢复，此时尽管 $R (y_{r}, Φ, x_{r})$ 相对较小，但仍然保持了一个远大于ζ的值，意味着还有部分信息并没有得到充分解码，但传统方法会将其认作是噪声而被丢弃，导致光谱细节信息也受到了损失。但本文提出的JDSD算法会将残差向量认作是对高频信息的采样，而对其继续在空域求解，以对光谱恢复进行细节补偿。另一方面，在JDSD算法中，当σ设置较小时，频域恢复的结果与真实信号相差较远，会导致 $R (y_{r}, Φ, x_{r})$ 很大，但此时可以认为高频信息得到足够多的保留，并将其放在高频部分进行恢复，对式（8）的求解可以大大补偿信号的损失。因此，相比较于传统但稀疏域变换算法而言，JDSD对于参数σ的可接受调整范围更大，具有更好的鲁棒性，将在2.2部分进行详细讨论。

综上所述，本文提出的JDSD算法流程如表1所示。

表1 双稀疏域联合求解算法（JDSD）

Table 1 Method Based on Joint of Double Sparse Domains

输入： $y, Φ$

输出： $x$

1：设置σ，求解：

2：设置ζ，求解：

其中：

$y_{h} = y - y_{l}$

$y_{l} = Φ {\hat{x}}_{l} = Φ Ψ {\hat{θ}}_{l} = Θ {\hat{θ}}_{l}$

3、返回：

，其中：

2 实验结果与讨论

2.1　数据来源和评价方法

验证算法的过程中，仍然以ENVI光谱库巴西棕榈蜡（Carnauba Wax）在0.4 μm~2.5 μm波段的反射率光谱曲线作为测试数据（图1），其在1.208 μm，1.388 μm和1.728 μm处有三个明显的吸收特征峰。同时为了说明算法的普适性，本文从该光谱库中随机选取了500种地物的光谱曲线，进行了统计测试，同时为了统一数据处理，对所有数据进行了等间隔重采样，保留512个波段作为原始信号进行仿真。

光谱曲线特征体现在光谱线形中，因此本文涉及所用光谱曲线都进行了归一化。光谱曲线恢复精度具有很多评价方法，其中光谱角度匹配SAM（Spectral Angle Model）可以直观地反映两条线形的匹配程度^［16］：

S A M = \frac{x_{1} x_{2}}{‖x_{1}‖ ‖x_{2}‖}

（11）

式（11）中SAM的值表示两个向量 $x_{1}$ 和 $x_{2}$ 的夹角余弦值，取值范围为 $[0,1]$ ，两个向量的形状越是接近，其夹角越小，此时SAM越大，当两条曲线完全相同时，夹角为0，对应的SAM值为1，因此SAM的值越接近于1，两条曲线就越接近。可以通过评价恢复光谱与原始光谱之间的SAM，来评价恢复算法的光谱线形保真度。

另一方面，本文提出算法优势在于可以对光谱吸收峰等细节特征进行补偿，SAM仅仅能体现对于线形匹配度的保持，一般其值会比较高，难以体现对于光谱特征的保持，因此本文用 $G S A M$ （Gradient Spectral Angle Model）梯度光谱角度匹配对光谱特征保真度进行评价，其计算公式为：

G S A M = \frac{< \nabla x_{1}, \nabla x_{2} >}{∥ \nabla x_{1} ∥ ∥ \nabla x_{2} ∥}

（12）

GASM利用光谱曲线的梯度信息来增强光谱特征在比较中的权重，进而客观评价算法对光谱特征保真度的影响。式（12）中， $\nabla x_{1}, \nabla x_{2} \in [- 1,1]$ ，因此 $G S A M \in [- 1,1]$ 。当 $G S A M = - 1$ 时， $\nabla x_{1} 与 \nabla x_{2}$ 的所有相同位置的值幅值相同，符号相反，这种情况在现实光谱比较中基本是不存在的，本文涉及情况下，待比较的两条光谱曲线具有同源性，即在 $\nabla x_{1}, \nabla x_{2}$ 对应相同位置的元素符号相同，此时 $G S A M$ 取值范围为［0，1］，其值越接近1，表明x₁和x₂的相似度越高。

2.2　实验结果与讨论

DCT（Discrete cosine transform）变换基在现实应用中具有良好的稀疏化效果，并能具有代表性地反映数据稀疏化特点，因此首先将其作为基本稀疏基进行算法验证，有关更多的稀疏基对算法的影响将在2.2.2节讨论。OMP（Orthogonal Matching Pursuit）是一种常用的重构算法，可以很好地恢复稀疏信号，其对于信号稀疏性的依赖性具有比较直观地反映，因此为了说明参数的影响，本文首先利用其作为单变换域恢复算法进行重构，进而与本文JDSD算法重构结果进行对比。综上，在基本实验结果分析中，利用DCT作为稀疏基，利用OMP重构轮廓+IRLS补偿细节的组合方式验证所提算法，记作JDSD1。更多地算法组合方式将在2.2.3节讨论。同时指出，本文所有仿真都是在Windows10 （64位）操作系统平台下进行，所使用处理器为Intel（R） Core（TM） i7-8550，内存为8G，仿真软件为MATLAB 2017a （64位版）。

2.2.1　基本实验结果

对于单稀疏域变换而言，稀疏约束决定了恢复结果趋于平滑的程度。用OMP解式（8）的过程中，σ的大小限定了恢复结果非零元素个数，表征信号在频域系数的0范数，对应公式（7）中L=0的情况。对σ的设置直接决定着重构结果的质量。如图2所示，第一列为探测目标（巴西棕榈蜡）经DCT变换的频域组成，σ的值表示取低频系数的个数，第二列给出了σ分别取10、25、250时， OMP算法能恢复的最好结果（100%采样率的情况下），第三列表示恢复结果与原始信号的偏差。从图中可以看出，当σ=10时，只有信号的轮廓得到了恢复，随着σ的增加，恢复结果精度在逐渐提高，但恢复曲线仍然保持一定的平滑特性而导致光谱细节被丢弃或者弱化。本质上，当对σ设定一个值时，恢复结果最多保留前σ大的系数，DCT变换后的频域内，较大的系数基本全部集中在低频部分，因此可以简单地认为，σ设置的大小，直接决定了低频部分保留的多少，当σ为一个较小的数值时，因为保留系数个数不够，而使得重构结果远远偏离原始信号。另一方面，当σ的值足够大而足以涵盖部分高频段时，又因高频系数过小而被丢掉，这就是图中即使σ到达250的时候，光谱细节仍然没有被很好地恢复，实际上，我们在实验中即使将σ设置为最大（512），恢复结果也基本与σ=25时结果相同。因此包含OMP在内的单稀疏域变化算法本身决定了重构结果更偏向于平滑部分的恢复，而细节被严重损失。

图2 分别为10 、25、250时对应信号的高低频分布

Fig. 2 igh and low frequency distribution of corresponding signals when σ is 10, 25, 250

从以上分析可以看出，当仅仅考虑DCT变换域的时候，由于系数分布具有比较单一的特点，而难以同时兼顾高低频信息的恢复。另一方面，当信号的轮廓信息得到恢复后，其残差向量在空间域表现出了良好的稀疏性，只在少数位置具有比较大的值，较小幅值的系数往往可以被认为是噪声而被丢掉。仍使用0范数客观地表征信号的稀疏性，在丢掉小于0.01系数的前提下，图2中三种条件的残差向量对应的0范数分别为52，19和10，而残差向量维度为512，可以明确地看出其良好的稀疏性。因此可以直接对残差向量进行恢复。图3中给出了在设置不同σ的条件下，利用OMP和JDSD1算法恢复结果的对比。

图3 同采样率在不同稀疏度约束下的恢复结果 (a) 采样率：40%，σ: 10，(b) 采样率：40%，σ: 25，(c) 采样率：40%，σ: 250，(d) 采样率：80%，σ: 10，(e) 采样率：80%，σ: 25，(f) 采样率：80%，σ: 250

Fig.3 Recovery results of different sampling rates under different sparsity constraints (a) Sampling rate: 40%, σ: 10, (b) Sampling rate: 40%, σ: 25, (c) Sampling rate: 40%, σ: 250,(d) Sampling rate: 80%, σ: 10, (e) Sampling rate: 80%, σ: 25, (f) Sampling rate: 80%, σ: 250

从图3中可以看出，JDSD1算法对OMP算法恢复结果很好地进行了细节补偿，在相同采样率下，JDSD1恢复的结果精度明显高于OMP恢复结果，同时可以看出，在相同采样率下，三种σ取值下JDSD1的恢复结果相差不大，说明该算法对于稀疏度的限制更小，具有更好的鲁棒性。图4中给出了ENVI光谱库中随机选取了500种地物的光谱曲线，进行的实验统计结果。图4（a）和图4 （b）分别给出了当σ=25时，两种算法在不同采样率下的恢复结果的SAM和GSAM对比，可以明显地看出，JDSD1比OMP具有更好的算法性能。特别地，在低采样率的条件下仍然保持了较好的恢复结果，如20%采样率条件下SAM可以达到0.85，GSAM由OMP的0.47提升到0.56，80%采样率条件下，光谱恢复结果SAM可以达到0.96以上，GSAM达到了0.90以上，表明光谱线形保真度和光谱特征保真度都得到了较大的提升。图4（c）和图4（d）给出了40%采样率条件下，不同σ的取值对于两种算法恢复结果的对比，从曲线的斜率变化程度上看，不同σ取值对于JDSD1的影响远远小于对OMP的影响，进一步证明了JDSD1对于参数的容忍度更强，具有良好的鲁棒性。

图4 对500种样品测试的结果 (a)（b）:不同采样率下恢复精度对比（σ：25），（c）（d）: σ 对恢复精度的影响（采样率：40%）

Fig.4 Test results on 500 samples (a) (b): Comparison of recovery accuracy at different sampling rates (σ: 15) ,(c) (d):σ’s impact on recovery accuracy (Sampling rate: 40%)

2.2.2　稀疏基的影响

尽管JDSD1比OMP对于参数σ的容忍度更强，但是σ的设置仍然需要根据信号的稀疏特点设定。信号在不同的稀疏基下具有不同的稀疏特点，这其中包括稀疏系数的幅值大小和位置分布。如图5（a）所示，为巴西棕榈蜡光谱曲线对应DCT变换和小波变换（DWT：db8）的系数分布，其中小波是3个不同层数分解对应的结果。可以看出，不同的稀疏分解结果中，“大系数”的集中程度具有很大的差别。一方面，不同的稀疏基与对应测量矩阵的相关性匹配不同，另一方面，稀疏系数的分布影响着JDSD算法中对参数σ的设置，这都将直接影响重构结果。图5（b）和5（c）分别给出了σ设置为10和100情况下，信号重构过程中的分布情况，显然，当σ较小时，DWT频域内信息仅仅得到了少量恢复，导致残差较大，已经不具备空域稀疏的特点，无法再按照上述过程进行细节补偿。

图5 巴西棕榈蜡光谱曲线对应不同稀疏变换的信息分布 (a):不同稀疏变换系数分布，（b）：σ =10时信号重构过程中的分布情况，(c):σ =100时信号重构过程中的分布情况

Fig.5 Information distribution of carnauba wax spectrum corresponding to different sparse transforms (a) Distribution of different sparse transform coefficients， (b): Distribution during signal reconstruction at σ=10， (c) Distribution during signal reconstruction at σ=100

(b)	OMP 全采样恢复
	残差向量
(c)	OMP 全采样恢复
(c)	残差向量
		DCT	DWT(db8-5)	DWT(db8-4)	DWT(db8-3)

图6给出了四种稀疏基下对500种样本测试的恢复结果，分别对应低采样率（40%）和高采样率（80%）两种情况。结合小波变换理论，在合理的范围内，一般较大层数的小波分解，可以使得系数较大程度的集中，这将增加信号的稀疏度，更有利于信号重构，因此这里仅针对最大层数分解的小波变换进行对比研究。本文选择db4，sysm4，coif4小波进行分析，这三种小波具有不同的特点，是现实信号处理中常用的稀疏基。从图中更可以看出，除了DCT变换外，其它稀疏变换下，σ值的设定对JDSD1的重构结果具有很大的影响，尤其是在低采样率条件下，如果σ值设置过小，会大大降低重构结果的光谱线形保真度（SAM）和光谱特征保真度（GSAM）。另一方面，可以明显看出，当σ大于某个阈值时，重构精度会急剧上升，并在大于该阈值的范围内，对σ变化的敏感度会下降，因此通常可以设置σ为一个较大的值来提高光谱恢复精度，对于本文测试数据，当设置 $σ \geq 6$ 4时，对于四种稀疏基，40%采样率条件下，SAM和GSAM分别都可以达到0.89和0.71以上，80%采样率条件下SAM和GSAM分别可以达到0.95和0.90以上。这是因为当σ足够大时，可以涵盖全部低频系数而使得光谱轮廓信息可以得到充分恢复，进而可以利用JDSD1算法进行更好的进行细节补偿。

图6 σ对不同稀疏分解恢复结果的影响 (a)（b）采样率为40%条件下，（c）（d）采样率为80%条件下

Fig.6 Effect of σ on the recovery results of different sparse decompositions. (a) (b) at a sampling rate of 40%， (c) (d) at a sampling rate of 80%

另一方面，尽管可以通过为σ设定为一个较大的值来提高光谱恢复精度，但是当σ值较大时，会大大降低JDSD算法的速度，如图7所示，为测试样本在不同σ条件下恢复一条512维光谱曲线所需要的平均时间，可以看出随着σ的增大，JDSD1算法速度大大降低。这是因为 OMP作为JDSD1的组成部分之一，它是一种贪婪算法，在σ较大时，需要恢复的系数个数增多，随着系数的增多，该类算法速度将大大降低。因此在现实应用中，需要根据实际需求来设定参数σ以保证恢复速度和精度的平衡。

图7 σ对不同稀疏分解恢复速度的影响(a)（b）采样率为40%条件下，（c）（d）采样率为80%条件下

Fig.7 Effect of σ on recovery speed of different sparse decompositions. (a) (b) at a sampling rate of 40%，(c) (d) at a sampling rate of 80%

2.2.3　不同重构组合的JDSD算法

除了OMP和IRLS算法外，传统单稀疏域重构算法还有很多，如两步软阈值迭代阈值（TwIST， Two step iterative shrinkage/thresholding）^［

12］，梯度投影（GPSR，Gradient Projection for Sparse Reconstruction：）^{［参考文献 13

百度学术}13］，这些方法两两组合也可以形成JDSD算法。在不同算法中，对于稀疏约束具有不同的表现形式，对应的参数设置也具有一定差别，但相互之间都具有等价性，本质上都是保留适量测量残差，利用其进行细节补偿。表2给出了几种单稀疏变换域重构算法结果，以及它们分别与IRLS组合的JDSD恢复结果，表中SAM和GSAM同样是对500种样本恢复的平均值。如表所示，20%采样率条件下，四种单稀疏域变换重构算法平均恢复结果SAM和GSAM分别为0.625和0.515，而三种JDSD平均重构结果对应的值分别为0.817和0.659，40%和80%采样率情况下，JDSD算法同样大大提高了光谱恢复精度，充分证明了JDSD算法的优越性。另一方面，从整体上看，JDSD1的恢复结果要好于JDSD2和JDSD3，这是因为OMP善于恢复平滑信号，IRLS善于弥补残差细节，两者达到了很好地优势互补，而TwIST和GPSR是针对信号整体进行重构，不能对高低频信号进行很好的分解，使得残差所含信息具有较大的噪声，直接降低了IRLS的细节补偿效果。

表2 不同组合的JDSD算法恢复结果对比

Table 2 Comparison of JDSD algorithm recovery results of different combinations

算法名称	采样率	DCT		DWT(db4)		DWT(sysm4)		DWT(coif4)		Average
算法名称	采样率	SAM	GSAM	SAM	GSAM	SAM	GSAM	SAM	GSAM	SAM	GSAM
OMP	20%	0.654	0.479	0.663	0.482	0.677	0.579	0.664	0.445	0.665	0.496
	40%	0.723	0.574	0.733	0.597	0.734	0.677	0.733	0.554	0.738	0.600
	80%	0.918	0.792	0.926	0.823	0.921	0.893	0.883	0.807	0.912	0.829
IRLS	20%	0.662	0.481	0.670	0.583	0.680	0.493	0.631	0.443	0.614	0.500
	40%	0.716	0.592	0.726	0.697	0.734	0.607	0.700	0.561	0.719	0.614
	80%	0.898	0.772	0.902	0.879	0.920	0.882	0.858	0.747	0.895	0.820
TwIST	20%	0.596	0.556	0.586	0.577	0.622	0.572	0.583	0.544	0.597	0.562
	40%	0.668	0.647	0.679	0.628	0.718	0.672	0.651	0.633	0.679	0.645
	80%	0.843	0.830	0.881	0.866	0.903	0.850	0.832	0.827	0.865	0.843
GPSR	20%	0.624	0.504	0.636	0.544	0.624	0.482	0.615	0.482	0.625	0.503
	40%	0.697	0.547	0.717	0.583	0.699	0.547	0.686	0.533	0.700	0.553
	80%	0.770	0.692	0.798	0.752	0.791	0.632	0.764	0.684	0.781	0.690
JDSD1 (OMP+IRLS)	20%	0.851	0.561	0.887	0.865	0.891	0.661	0.889	0.875	0.880	0.741
	40%	0.895	0.712	0.932	0.919	0.942	0.880	0.941	0.900	0.928	0.853
	80%	0.962	0.871	0.942	0.932	0.952	0.892	0.963	0.921	0.955	0.904
JDSD2 (OMP+TwIST)	20%	0.856	0.663	0.863	0.667	0.756	0.663	0.756	0.593	0.808	0.647
	40%	0.932	0.752	0.941	0.764	0.832	0.732	0.872	0.792	0.894	0.760
	80%	0.972	0.891	0.980	0.923	0.882	0.879	0.942	0.894	0.944	0.897
JDSD3 (OMP+ GPSR)	20%	0.831	0.561	0.827	0.661	0.731	0.541	0.661	0.591	0.762	0.589
	40%	0.865	0.700	0.853	0.779	0.865	0.765	0.765	0.724	0.837	0.742
	80%	0.942	0.847	0.936	0.907	0.902	0.837	0.902	0.877	0.921	0.891

2.2.4　计算光谱成像实验

如图8（a）为利用编码板作为编码器件搭建的光谱成像系统，其中编码板是在镀铬层的玻璃基板上光刻编码制成，色散部件是SpecIm公司用在V10E成像光谱仪上的PGP组件。当编码板只保留单条狭缝时，系统就变成了传统的推扫式高光谱成像系统（PHI，Push-broom Hyperspectral Imager），当编码按照高斯随机排列时，系统就变成了单色散编码孔径光谱成像系统（SD-CASSI）。图8（b）为80%采样率条件下SD-CASSI对外场景成像结果合成的彩色图，所使用恢复算法为本文提出的JDSD。图8（c）为在相同条件下PHI 的成像结果，两者采样时间相差10s左右。取两图相同位置均匀区域（位置A）计算信噪比，图8（b）为39.19db，而图8（c）为39.17db，因此可以得出结论，欠采样条件下的SD-CASSI系统利用JDSD1算法已经取得了和传统系统全采样条件下相当的成像效果。位置B为一个具有明显光谱特征（红色）的目标，可以作为光谱恢复评价的测试点。JDSD1比JDSD2和JDSD3具有更好的恢复效果，这里仅仅用其结果与OMP恢复结果对比。图8（d）~（g）分别为20%，40%， 80%和100%采样率条件下对测试点恢复的光谱曲线结果，并给出了与PHI采样结果的对比。可以看出，JDSD1比OMP可以更好的恢复光谱细节，比如在40%采样条件下，利用OMP算法时，760nm的氧气吸收峰没有得到很好的恢复，大大降低了对于光谱探测的可信度，而JDSD1通过对细节的补偿，可以很好的保持光谱特征。随着采样率的上升，JDSD1的恢复结果明显优于OMP算法，在760nm的氧气吸收带具有更深的深度，以及目标在红色光谱段（620~780nm）恢复出了更强的光谱特征信息。从光谱恢复结果上同样可以看出，80%采样时，SD-CASSI利用JDSD1算法已经取得了和PHI相当的结果，100%采样时，SD-CASSI系统的光谱深度优于PHI系统，进一步说明了JDSD1对于光谱细节补偿的作用。

图8 实验室验证装置和验证结果（a）CASSI系统，（b）80%采样恢复真彩色图，(c)PHI 成像结果（d）~（g）20%，40%，80%，100%采样率条件下两种算法对位置B的恢复结果

Fig.8 Laboratory verification equipment and verification results (a) CASSI system, (b) recovered true color image by 80% sampling, (c) PHI imaging results, (d)~（g） Recovered results by two algorithms under 20%，40%，80% and 100%sampling

3 结论

提出了基于双稀疏域联合求解（JDSD）的方法，首先利用OMP算法对欠采样信号在频域内求解，然后针对因光谱吸收峰的存在带来的残差向量具有空间稀疏性的特点，利用IRLS进一步解码，达到细节补偿的目的，使得重构结果精度大大提高。JDSD算法参数σ的设置是决定重构结果的直接因素，影响其设置的主要因素是稀疏基的选择，不同的稀疏分解决定了信号的系数分布特点，进而决定了应该设置σ的大小。一般地，σ需要选择较大的值，以保证有足够多的轮廓信息被恢复。但是随着σ的增大，JDSD算法的速度急剧下降，这是因为贪婪算法作为JDSD的组成部分，直接影响了恢复速度。因此在实际应用中，需要根据实际需求来设定σ以平衡算法速度和恢复精度。不同的算法可以组合成更多的JDSD形式，本文选择了三种JDSD算法形式，并对选定的大量实验样本进行了测试，利用SAM和GSAM分别评价了恢复算法的光谱线形保真度和光谱特征保真度。测试结果表明，在20%采样率条件下，传统单稀疏域恢复算法的平均恢复结果的SAM和GSAM分别为0.625和0.515，而三种JDSD算法平均重构结果对应的值分别为0.817和0.659，随着采样率的提升，JDSD的恢复结果都远远优于传统算法。最后，利用搭建的编码孔径成像系统进行了对外实验， JDSD算法可以得到更好的光谱特征恢复，为实际光谱应用中提供了更好的基础。另一方面，针对JDSD算法的性能定量化分析和工程应用还需要进一步研究，主要有以下几个方面。

（1）现实应用中还存在更多的重构算法和稀疏表达方法，对信号进行更好的特征挖掘以寻求更好的JDSD组合形式，以及对该类算法速度的提升的是未来研究的重点方向之一。

（2） JDSD方法利用采样残差向量来表示高频信号，但是噪声与吸收峰同时属于高频信息，在实际应用中需要进行有效的区分。如图3（a）所示，在较低采样率条件下，JDSD算法恢复的结果线形更接近真值，但同时引入了一些噪声信息，因此需要更深入的研究来区分和处理有效信息和噪声信息。

（3）从应用角度上，本文是基于光谱数据相关性进行的算法，实际上，光谱成像领域中，还存在二维空间图像数据，在未来应用研究中可以进一步挖掘空间和光谱之间的相关性，进而提高算法的适用性和鲁棒性。

References

Donoho D L. Compressed sensing［J］. IEEE Transactions on Information Theory， 2006， 52（4）： 1289-1306. [百度学术]

Candes E J， Wakin M B. An Introduction to Compressive Sampling［J］. IEEE Signal Processing Magazine， 2008， 25（2）： 21-30. [百度学术]

Duarte M F， Davenport M A， Takhar D， et al. Single-Pixel Imaging Via Compressive Sampling［J］. IEEE Signal Processing Magazine， 2008， 25（2）： 83–91. [百度学术]

Sun T，Kelly K. Compressive sensing hyper spectral imager［C］// Computational Optical Sensing and Imaging， October 13-15，2009，San Jose， California ，United States Washington： Optical Society of America， 2009：Ctua5. [百度学术]

Yitzhak A， Chaim V， Yair R. et.al， Compressive hyper spectral imaging by random separable projections in both the spatial and the spectral domains［J］.Applied Optics， 2013， 52（10）：D46-D54. [百度学术]

Wagadarikar A， John R， Willett R， et al. Single Disperser Design for Coded Aperture Snapshot Spectral Imaging［J］. Applied Optics. 2008， 47（10）： B44–B51. [百度学术]

Kittle D， Choi K， Wagadarikar A，et.al. Multiframe image estimation for coded aperture snapshot spectral imagers［J］， Applied Optics. 2010， 49（36）：6824-6833. [百度学术]

Gehm M E， John R， Brady D J， et.al， Single-Shot Compressive Spectral Imaging with a Dual-Disperser Architecture［J］， Optics Express. 2007， 15（21）： 14013-14027. [百度学术]

Henry Arguello， Gonzalo Arce ，Code Aperture Design For Compressive Spectral Imaging［C］//18th European Signal Processing Cofference（EUSIPCO-2010）Aalborg， Denmark， August 23-27，2010：1434-1438. [百度学术]

Tropp J， Gilbert A. Signal Recovery from Random Measurements via Orthogonal Matching Pursuit［J］. IEEE Transactions on Information Theory， 2007， 53（12）： 4655-4666. [百度学术]

Chartrand R， Yin W T ， Iteratively Reweighted Algorithms For Compressive Sensing［C］. 2008 IEEE International Conference on Acoustics， Speech and Signal Processing.Las Vegas， NV， USA，2008：1-4. [百度学术]

Bioucas J D and Figueiredo M， A new TwIST： Two step iterative shrinkage/thresholding algorithms for image restoration［J］ IEEE Trans. Image Process， 2007，16（12）： 2992–3004. [百度学术]

Figueiredo A.T ， Robert D N， Wright S J， Gradient Projection for Sparse Reconstruction： Application to Compressed Sensing and Other Inverse Problems， IEEE Journal Of Selected Topics In Signal Processing， 2007，1（4）：586-597.. [百度学术]

Lu Xu-Liang， JIA Qi，RONG Xian-hui，et.al， Application of degree of grey incidence on similarity of spectral curves［J］， Journal of PLA University of Science and Technology（Natural Science Edition）（吕绪良，贾其，荣先辉等，灰色关联度在光谱曲线相似度分析中的应用，解放军理工大学学报（自然科学版））.2011， 12（5）：496-500. [百度学术]

您是第12687457 位访问者

主管单位：中国科学院

主办单位：中国科学院上海技术物理研究所，中国光学学会

地址：上海市玉田路500号电话：021-25051553

51La

首页

学报简介

编委会

征稿简则

版权声明

开放获取

出版道德声明

相关下载

联系我们

常见问题

English Version

基于双稀疏域联合求解的高精度光谱恢复算法 PDF

摘要

关键词

引言