王业超 陈晓丽 钟晓明 赵海博 张丽莎 苏云

基于全变分正则项的CASSI数据重构算法

王业超1,2陈晓丽1,2钟晓明1,2赵海博1,2张丽莎1,2苏云1,2

（1 北京空间机电研究所，北京 100094） （2 先进光学遥感技术北京市重点试验室，北京 100094）

在编码孔径快照光谱仪（Coded Aperture Snapshot Spectral Imager，CASSI）光谱重构算法设计中，两步迭代收缩阈值（Two-Step Iterative Shrinkage/Threshold，TwIST）算法实现了CASSI的光谱重构，但因正则项选取和引入的噪声随迭代次数而不断放大，导致其光谱重构精度低于80%。文章在TwIST算法基础上，以光谱图像具有空间平滑过渡特性为先验知识，提出两点改进：一是选择全变分正则约束项；二是对每一步迭代的更新项进行全变分去噪处理。为了验证改进后的算法，文章通过计算机仿真CASSI的采集数据，得出仿真数据光谱重构精度为90.93%；并根据CASSI样机采集试验数据，得出试验数据的光谱重构精度为86.56%。改进后的算法可为以后CASSI数据重构提供参考。

两步迭代收缩阈值 压缩感知 编码孔径 全变分 光谱重构 光谱仪

0 引言

编码孔径快照光谱仪（Coded Aperture Snapshot Spectral Imager，CASSI）是在压缩感知理论框架下发展起来的。压缩感知理论指出若信号可压缩或在某个变换域稀疏，则利用与变换基不相关的观测矩阵将高维信号投影至低维空间，可通过求解优化问题从少量投影信息中以高概率重构出原始信号，并证明投影信息包含重构所需的足够信息[1-5]。在压缩感知理论框架下，只要信号满足可压缩或稀疏的条件，即可通过远低于奈奎斯特采样频率采样，重构出原始信号[6-8]。上述两个结论为CASSI算法设计提供理论支撑。基于压缩感知理论，文献[9-13]提出并发展CASSI及其光谱重构算法的设计，但光谱重构精度较低。CASSI系统核心要素和难点是从探测器采样的二维混叠图中重构三维光谱数据立方体。文献[14-20]提出通过硬件与算法联合设计的思路重构CASSI光谱数据。受联合设计启发，研究人员在重构算法中加入全色参照信息提高光谱重构精度，但增加了系统复杂程度[21-22]。目前对重构算法的研究主要包括贪婪迭代算法、凸优化算法和基于贝叶斯框架的重构算法，前两种算法在处理一维信号及二维图像信号中取得了较好的成果，第三种算法由于其时间关联性较强而不适用于光谱图像的重构[23-26]。三维光谱数据具有维度高、矩阵规模大的特点。贪婪迭代算法在求解小规模一维和二维信号优化问题时获得较好结果，但其在每一步都要完成信号或者残差信号在冗余字典每一个原子上的投影计算，计算量成倍增加，重构时间较长，因而贪婪迭代算法不适用于混叠光谱数据重构求解。文献[27]提出两步迭代收缩阈值（Two-Step Iterative Shrinkage/Threshold，TwIST）算法实现CASSI系统光谱重构，但光谱重构精度较低。文献[28]提出结合图像自回归模型的重构算法，因光谱图像模型不同于二维图像，该方法不适用于光谱数据重构。本文在分析TwIST算法基础上，引入两点改进：一是将求解函数中正则项设为光谱图像的全变分正则约束项，使得求解结果最大可能的保留图像本身的高频信息；二是在每次迭代计算前，对更新项进行全变分去噪处理，避免引入的噪声在后续迭代中放大。为了验证基于全变分正则项的TwIST算法对CASSI数据具有较高的光谱重构精度，本文对仿真数据和试验数据进行了重构。

1 CASSI数据构成

编码孔径快照式光谱仪数据可以表示为

编码模板的空间调制使光谱数据损失部分空间信息，棱镜色散后探测器同一像素上采集到不同空间位置光谱的混叠信息，因此重构过程是从不完整的观测中恢复全部光谱数据，需求解严重的病态方程，可通过求解优化问题解决，目标函数表示为

式（2）右边第一项为系统保真项，衡量优化结果与系统观测之间的误差；第二项为正则项，一般根据求解目标的固有性质对求解过程进行约束。

正则项对光谱重构精度有较大影响。稀疏重构的梯度投影算法以场景信息在空间上分段平滑作为先验知识，实现优化算法的重构[29]。基于全变分的正则约束项可以表示为

式中 （，）为离散后的空间坐标；为离散后的谱段；F,n,k，F+1,n,k，F,n+1,k分别为空间坐标为（,），（1,），（,+1）谱段的光谱数据。

本文采用全变分正则项作为求解优化问题的约束，在抑制噪声的同时防止损失过多细节信息来保证光谱重构精度。

2 基于全变分正则项的TwIST算法

式（2）中的优化问题可由TwIST算法求解。本文通过改进TwIST算法，提高光谱重构精度。

2.1 TwIST算法

TwIST算法由迭代收缩阈值算法（Iterative Shrinkage/Thresholding，IST）和迭代加权收缩算法（Iterative Re-Weighted Shrinkage，IRS）结合而来，核心思想是利用前两个估计值来更新当前值，即“两步”迭代，主要迭代更新过程如下：

2.2 改进后的TwIST算法

全变分正则项基于光谱图像空间分段平滑过渡的先验知识[30]，对病态方程的求解施加约束。本文对TwIST算法提出两点改进：一是将目标函数的约束项指定为光谱图像的全变分正则约束项，使最终重构近似解保留图像更多高频细节，同时实现较快的收敛；二是在每次迭代计算前，对更新项进行全变分去噪处理，以免引入的噪声随迭代计算而放大。

式中（•）为光谱数据的目标函数。

图1 改进的TwIST算法流程图

3 仿真数据的重构及分析

3.1 建立仿真模型

CASSI系统数学模型为

式中 (,)表示连续的空间坐标；(,)为获取的图像；(,,)表示输入的原始光谱信息；(,,)表示编码模板空间调制和色散棱镜光谱调制；(,)表示系统引入的各种噪声。

离散形式可以表示为

通过Matlab实现CASSI成像的仿真建模。以随机的0/1矩阵建立编码模板的空间调制，采用理想线性色散模型。以高光谱数据集作为输入，获得的仿真图像分辨率667像素×520像素，共30个谱段，谱段范围426~670nm。图2为仿真模拟采集结果。

图2 仿真模拟采集结果

图2（a）仿真目标场景包含黑白棋盘格、标准色卡以及其他目标物；图2（b）仿真采集中光谱沿竖直方向色散，获取含光谱信息的编码混叠图。

图3为等效观测矩阵仿真结果，共30个谱段，谱段范围426~670nm。等效观测矩阵在CASSI系统中起到空间调制作用，在算法中作为等效观测矩阵参与光谱数据立方体重构。

图3 等效观测矩阵仿真结果

3.2 仿真光谱数据重构结果

图4是以文献[27]中的TwIST算法重构的仿真光谱数据。在426~496nm和587~670nm的光谱切片信息因噪声等影响而损失严重，光谱数据整体的重构品质较差[13,27]。经计算平均光谱重构精度为73.21%。综合不同仿真数据的重构结果，TwIST算法的光谱重构精度低于80%。

图4 TwIST算法重构光谱数据

图5为改进后的基于全变分正则约束的TwIST算法重构出的30个光谱切片。

图5 基于全变分正则项的TwIST算法仿真数据重构光谱

对重构后的光谱切片进行分析可以看出，图2（a）右侧的色卡在图5的光谱切片中清晰可见，色卡上不同块的亮暗表示当前谱段下不同颜色灰度值的强度；图2（a）左侧的黑白色棋盘格在图5不同光谱切片中边缘清晰可见，重构的光谱数据较好地保留了空间细节信息，验证了基于全变分正则项的TwIST算法对CASSI数据重构的可行性。

图6为重构结果分析。从30个谱段中选取红绿蓝波段合成彩色图，见图6（a），与图2（a）仿真目标场景相比，图像的色彩和空间信息均具有较好的一致性。图6（b）为场景中重构光谱曲线与原始光谱曲线的对比，重构光谱曲线与参照数据较好的吻合，验证了基于全变分正则项的TwIST算法具有较好的光谱重构精度。

图6 重构结果分析

选取图2（a）中色卡作为分析目标，色卡有24个不同颜色，编号如图7。

图7 色卡的分块编号

计算色卡各分块的光谱重构精度

式中Q为重构数据谱段在空间位置的灰度值与总灰度值比值；F为重构数据谱段在空间位置的灰度值；D表示谱段在空间位置的灰度值与总灰度值比值；C表示谱段在空间位置的灰度值。

根据式（13）、（14）、（15）计算图5各光谱切片色卡分块的光谱重构精度，见表1。

表1 仿真色卡分块重构精度

Tab.1 The reconstruction accuracy of simulation color blocks

根据表1可计算得出平均光谱重构精度为90.93%，优于TwIST算法的光谱重构精度，与输入数据具有较好的一致性。

4 试验数据的重构及分析

为了进一步验证改进后的TwIST算法对CASSI数据的光谱重构精度，需要搭建CASSI试验样机，并采集数据。图8为试验采集结果，图8（a）试验目标场景，包含黑白棋盘格、标准色卡以及其他目标物，图8（b）CASSI系统试验采集图，试验采集中光谱沿竖直方向色散，展示了编码后光谱的混叠结果。

图8 试验采集结果

图9为等效观测矩阵采集结果。图9（a）为467nm谱段等效观测矩阵的采集结果，图9（b）为等效观测矩阵局部细节。

图9 等效观测矩阵采集结果

图10为基于全变分正则项的TwIST算法重构出的光谱数据切片，共有34个谱段，谱段范围415~670nm。

图10 基于全变分正则项的TwIST算法试验数据重构光谱

对重构出的光谱切片进行分析可以看出，图8（a）右侧的色卡在图10中不同光谱切片中清晰可见，色卡分块的亮暗表示当前波长不同颜色灰度值的强度；图8（a）左侧的黑白色棋盘格在图10中不同光谱切片中边缘清晰可见，重构的光谱数据较好的保留了空间细节信息，这些验证了基于全变分正则项的TwIST算法对CASSI实测数据重构的有效性。

按式（12）、（13）、（14）计算图10各光谱切片中色卡分块光谱重构精度，见表2。

表2 试验色卡分块光谱重构精度

Tab.2 The reconstruction accuracy of experimental color blocks

根据表2可计算得到平均光谱重构精度为86.56%，优于TwIST算法的光谱重构精度，与试验实测数据具有较好的一致性。

5 结束语

本文基于全变分正则项，对TwIST算法提出两点改进，并通过对仿真数据和试验数据的重构，验证了改进的基于全变分正则项的TwIST算法具有较好的光谱重构精度。改进的基于全变分正则项的TwIST算法对仿真数据光谱重构精度为90.93%，对试验数据的光谱重构精度为86.56%，均优于TwIST算法的光谱重构精度。该优化设计的光谱重构算法可为CASSI光谱重构精度提升提供参考。

[1] DONOHO D L. Compressed Sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306.

[2] CANDES E J, WAKIN M B. An Introduction to Compressive Sampling[J]. IEEE Signal Processing Magazine, 2008, 25(2): 21-30.

[3] KUTYNIOK G. Compressed Sensing: Theory and Applications[J]. Corr, 2012, 52(4): 1289-1306.

[4] ELAD M. Optimized Projections for Compressed Sensing[J]. IEEE Transactions on Signal Processing, 2007, 55(12): 5695-5702.

[5] BLUMENSATH T, DAVIES M E. Iterative Hard Thresholding for Compressed Sensing[J]. Applied & Computational Harmonic Analysis, 2009, 27(3): 265-274.

[6] FIGUERIREDO M A T, NOWAK R D, WRIGH S J. Gradient Projection for Sparse Reconstruction: Application to Compressed Sensing and Other Inverse Problems[J]. IEEE Journal of Selected Topics in Signal Processing, 2008, 1(4): 586-597.

[7] SHI G M, LIU D H, GAO D H, et al. Advances in Theory and Application of Compressed Sensing[J]. Acta Electronica Sinica, 2009, 37(5): 1070-1081.

[8] 石光明, 刘丹华, 高大化, 等. 压缩感知理论及其研究进展[J]. 电子学报, 2009, 37(5): 1070-1081. SHI Guangming, LIU Dianhua, GAO Dahua. Advances in Theory and Application of Compressed Sensing[J]. ACTA Electronica Sinica, 2009, 37(5): 1070-1081. (in Chinese)

[9] WAGADARIKAR A, JOHN R, WILLETT R, et al. Single Disperser Design for Coded Aperture Snapshot Spectral Imaging[J]. Applied Optics, 2008, 47(10): 44-51.

[10] WAGADARIKAR A A, PITSIANIS N P, SUN X, et al. Single-shot Compressive Spectral Imaging with A Dual-disperser Architecture[J]. Optics Express, 2007, 15(21): 14013-14027.

[11] ARCE G R, BRADY D J, CARIN L, et al. Compressive Coded Aperture Spectral Imaging: An Introduction[J]. IEEE Signal Processing Magazine, 2014, 31(1): 105-115.

[12] WAGADARIKAR A A, PITSIANIS N P, SUN X, et al. Video Rate Spectral Imaging Using A Coded Aperture Snapshot Spectral Imager[J]. Optics Express, 2009, 17(8): 6368-6388.

[13] WAGADARIKAR A A, PITSIANIS N P, SUN X, et al. Spectral Image Estimation for Coded Aperture Snapshot Spectral Imagers[J]. Proceedings of SPIE - The International Society for Optical Engineering, 2008, 7076(36): 6824-6833.

[14] LLULL P, LIAO X, YUAN X, et al. Coded Aperture Compressive Temporal Imaging[J]. Optics Express, 2013, 21(9): 10526-10545.

[15] ARGUELLO H, ARCE G R. Colored Coded Aperture Design by Concentration of Measure in Compressive Spectral Imaging[J]. IEEE Transactions on Image Processing, 2014, 23(4): 1896-1908.

[16] KITTLE D. Design and Fabrication of an Ultraviolet-visible Coded Aperture Snapshot Spectral Imager[J]. Optical Engineering, 2012, 51(7): 071403.

[17] WANG Lizhi, LI Huan, GAO Dahua, et al. A Flexible Design for Coded Aperture Snapshot Spectral Imager[C]// SPIE Vol.9273, Optoelectronic Imaging and Multimedia Technology III, International Society for Optics and Photonics, 2014.

[18] 石光明, 刘丹华, 高大化. 随机编码感知的高分辨率遥感光谱计算成像[J]. 航天返回与遥感, 2011, 32(5): 60-66. SHI Guangming, LIU Danhua, GAO Dahua. High-resolution Computational Spectral Imaging of Remote Sensing Based on Coded Sensing[J]. Spacecraft Recovery & Remote Sensing, 2011, 32(5): 60-66. ( in Chinese)

[19] 李国. 基于双通道的压缩光谱成像及其重构算法GPU实现[D]. 西安: 西安电子科技大学, 2014. LI Guo. The Double-channel Compressed Spectral Imaging and Reconstruction Algorithm Implemented on GPU[D]. Xi'an: Xidian University. 2014. (in Chinese)

[20] RUEDA H, ARGUELLO H, ARCE G R. Compressive Spectral Imaging Based on Colored Coded Apertures[C]// 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), May 4-9, 2014, Florence, Italy. IEEE, 2014.

[21] WANG Lizhi，GAO Dahua，SHI Guangming，et al. Double-channel Compressive Spectral Imaging via Complementary Code Patterns[C]// 2013 IEEE International Conference on Signal Processing，Communication and Computing (ICSPCC), Aug 5-8, 2013, Kunming, China. IEEE, 2013: 1-5.

[22] WANG Lizhi, XIONG Zhiwei, GAO Dahua. et al. Dual-camera Design for Coded Aperture Snapshot Spectral Imaging[J]. Applied Optics, 2015, 54(4): 848-858.

[23] TROPP J A, GILBERT A C. Signal Recovery from Random Measurements Via Orthogonal Matching Pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655-4666.

[24] BOYD, VANDENBERGHE, FAYBUSOYICH. Convex Optimization[J]. IEEE Transactions on Automatic Control, 2006, 51(11): 1859-1859.

[25] CHAE B G, LEE S. Sparse-View CT Image Recovery Using Two-step Iterative Shrinkage-thresholding Algorithm[J]. ETRI Journal, 2015, 37(6): 1251-1258.

[26] LI Y, HU B L, WANG S, et al. Coded Aperture Spectral Data Restoration Based on Two-step Iterative Shrinkage/Threshold[J]. Spectroscopy and Spectral Analysis, 2014, 34 (3): 847-850.

[27] BIOUCAS-DIA J M, FIGUEIREDO M A T. A New TwIST: Two-step Iterative Shrinkage/Thresholding Algorithms for Image Restoration[J]. IEEE Transactions on Image Processing, 2007, 16(12): 2992-3004.

[28] WU Xiaolin, ZHANG Xiangjun, WANG Jia. Model-guided Adaptive Recovery of Compressive Sensing[C] //2009 Data Compression Conference (DCC 2009), March 16-18, 2009, Snowbird, UT, USA. IEEE, 2009: 123-132.

[29] XU H K. Averaged Mappings and the Gradient-projection Algorithm[J]. Journal of Optimization Theory and Applications, 2011, 150(2): 360-378.

[30] JIA R Q, ZHAO H. A Fast Algorithm for The Total Variation Model of Image Denoising[J]. Advances in Computational Mathematics, 2010, 33(2): 231-241.

A Reconstruction Algorithm of CASSI Data Based on Total Variation Regular Terms

WANG Yechao1,2CHEN Xiaoli1,2ZHONG Xiaoming1,2ZHAO Haibo1,2ZHANG Lisha1,2SU Yun1,2

（1 Beijing Institute of Space Mechanics & Electricity, Beijing 100094, China） （2 Key Laboratory for Advanced Optical Remote Sensing Technology of Beijing, Beijing 100094, China）

In designing spectral reconstruction algorithm for Coded Aperture Snapshot Spectral Imager (CASSI), the two-step iterative shrinkage/threshold (TwIST) algorithm can realize spectral reconstruction of CASSI, but the accuracy of spectral reconstruction is lower than 80% due to selection of regular terms and amplification of the introduced noise with the increase of iteration number. In this paper, based on the study of the traditional TwIST algorithm, two improvements are proposed, selecting the total variation regular constraint terms and denoising the updated terms in each iteration. In order to verify the improved algorithm, the data of CASSI are simulated by computer, with the spectral reconstruction accuracy 90.93%. The CASSI prototype is built to collect the experimental data, and the spectral reconstruction accuracy of the experimental data is 86.56%. The research results provide a reference for the following CASSI data reconstruction.

two-step iterative shrinkage/threshold (TwIST); compressed sensing; coded aperture; total variation; spectral reconstruction; spectral imager

TP751

A

1009-8518(2020)01-0091-11

10.3969/j.issn.1009-8518.2020.01.011

2019-10-16

北京市科技计划课题（Z181100003018003）

王业超, 陈晓丽, 钟晓明, 等. 基于全变分正则项的CASSI数据重构算法[J]. 航天返回与遥感, 2020, 41(1): 91-101.

WANG Yechao, CHEN Xiaoli, ZHONG Xiaoming, et al. A Reconstruction Algorithm of CASSI Data Based on Total Variation Regular Term[J]. Spacecraft Recovery & Remote Sensing, 2020, 41(1): 91-101. (in Chinese)

王业超，男，1993年生，2016年获北京航空航天大学校飞行器设计专业学士学位，现在中国空间技术研究院飞行器设计专业攻读硕士学位。研究方向为计算光谱成像技术。E-mail：cast_wangyc_508@163.com。

（编辑：王丽霞）