卫星运动对压缩编码孔径光谱成像的影响

杨虹 黄远辉 钟宬 邵晓鹏

(西安电子科技大学技术物理学院，光电成像与图像处理实验室，西安 710071)

1 引言

随着遥感技术的发展，出现了多种光谱成像仪，例如摆扫成像仪、推扫成像仪和可调谐滤波成像仪[1-3]。这些成像仪的设计理念相对简单，但由于大都仅获得数据立方体的一维或二维子集，因而需要对其它维度进行实时扫描以此获取完整的数据立方体; 再者，对于非相干光源，它们的光子收集率都较低，导致信噪比也偏低。一些最新发展起来的光谱成像仪，虽然光子收集率的问题在很大程度上得到解决，但又引入了诸如“失踪三角锥”等一些新问题[4]。

上述这些传统的光谱成像仪均存在一个很重要的特点，那就是获得的观测值的总数大于或等于重建的数据立方体总元素个数。与此不同，Brady等人2006年首次提出了压缩编码孔径光谱成像的概念，这种成像方法基于压缩感知理论，得到的观测值总数远小于数据立方体的元素总数[5]。

压缩编码孔径成像是近年来基于压缩感知理论新兴的一种成像技术，它在可见光成像方面的应用研究相对较多，但在多光谱成像领域的研究相对较少。在保证图像品质的前提下由于压缩编码孔径成像技术能大大减少采样的数据量，使传统遥感光谱成像面临的最普遍的数据量庞大的问题得到解决，所以本文提出将压缩编码孔径成像应用于遥感成像。然而，现有的压缩编码孔径光谱成像系统仅能用于静态成像，而遥感成像是个动态成像的过程，因 此要想将压缩编码孔径应用于遥感光谱成像中，必须分析卫星平台运动对成像过程的影响，进而克服由于卫星平台运动导致的图像品质的下降。

2 压缩编码孔径光谱成像系统原理

2.1 压缩编码孔径原理简介

压缩编码孔径的灵感很大程度上来源于压缩感知这个新兴的领域[6-8]。压缩感知也被称为压缩采样或稀疏采样，是一种寻找欠定线性系统的稀疏解，利用稀疏的或可压缩的信号进行信号重建的方法，是近年来由D. Donoho，E. Candès，J. Romberg和T. Tao等人提出的一种全新的采样理论[9-11]。该理论指出: 只要信号是可压缩的或在某个域上是稀疏的，就可以用一个与变换基不相关的测量矩阵将变换所得高维信号投影到一个低维空间上，使实际采样率远低于 Nyquist采样率，然后通过求解一个优化问题从这些少量的投影中以高概率重构出原信号，这也证明了这样的投影包含了重构信号的足够信息[11]。

编码孔径成像过程分为两步: 第一步是编码过程，光线透过编码板每个小孔后在探测器上都形成投影图像，这些投影图像叠加在一起，在二维平面上形成退化的二维分布信号，即编码图像; 第二步是解码过程，对探测器采集到的图像数据进行滤波和重建，复原出原始目标图像[12]，其基本过程如图1所示:

图1 编码孔径成像过程Fig.1 Imaging procedure of Compressive Coded Aperture

2.2 压缩编码孔径快照光谱成像系统

2008年，Wagadarikar等人改进了他们在2007年提出的单分光元件编码孔孔径光谱成像仪[13]，其结构如图2所示。这是目前世界上发展相对较成熟的一套编码孔径光谱成像系统，由物镜、编码孔径、带通滤波器、F/8中继透镜、双阿米西棱镜和单色CCD探测器构成。物镜将场景成像于编码孔径平面，经编码孔径编码再经过带通滤波器，中继镜头和双阿米西棱镜到达探测器平面。

图2 单分光元件编码孔径快照光谱成像仪光路Fig.2 Optical path of Single Disperser Coded Aperture Snapshot Spectral Imager (SD-CASSI)

3 遥感卫星平台运动建模

随着科技的进步，人们对于光谱图像品质的要求越来越高，消除卫星平台运动引起的光谱混叠对提高光谱图像品质起着十分重要的作用。

3.1 卫星平台运动姿态及其对应像移的计算

卫星平台的运动在空间上通常有3种运动姿态: 俯仰、侧滚和偏航。唐秋艳等人对卫星运动姿态的改变对图像造成的影响进行研究[14]，并且引入了平均掺杂比的概念，架起了卫星平台参数与光谱成像之间的桥梁。在该项研究中，若将探测器阵列上一像素(m,n)处的光谱密度记为Bm(v)，其它像元的混叠光谱统一记作Bm+1(v)，混叠系数为iη，则(m,n)处的光谱密度可记为

由于平台的稳定系数一般是10-3(°)/s ～10-1(°)/s，即平台运动引起的混叠通常都不会超过一个像元间隔，所以可以将Bm+1(v)直接视为相邻像元的光谱密度，于是式(1)可改写为

由光谱混叠模型中混叠光谱各成分的比例等于各地物面积之比可知[14-15]，平均混叠权重系数由像移路径包围的总面积决定。

俯仰会在探测器的x方向产生像移，像移为Δx(t)=ftan[Δα(t)]; 侧滚在探测器的y方向产生像移，由图可得Δy(t)=ftan[Δβ(t)]; 偏航在水平和垂直方向均产生像移，但由于垂直方向像移较小，可近似为0，其像移主要集中在水平方向，且位移大小与偏航角度相关Δx(t)=y·Δθ(t)。其中,f为系统焦距; Δα(t)，Δβ(t)为t时刻的振动角度，由于卫星平台的稳定度较高，俯仰和侧滚的像移可分别近似表示成Δx(t)=f·Δα(t)，Δy(t)=f·Δβ(t)。

3.2 平均混叠权重系数矩阵的生成

在实际成像过程中，由于运动造成的像移往往是上述3种像移的叠加且卫星的稳定性导致像移不会超出一个像素范围，因此混叠光谱可简化为原始光谱与其周围8个相邻像素光谱的叠加[15]，如图3所示:

图3 8邻域像素光谱混叠示意Fig.3 Spectrum aliasing from the 8 neighbor pixels of the center pixel

观察可得，探测器上的像移有4种可能的情况，可描述为相邻8个光谱中的3个参与混叠，即光谱2,5,9参与混叠，光谱4,5,8参与混叠，光谱3,4,7参与混叠，光谱2,3,6参与混叠。

假设卫星平台的运动是一个随机过程，即在一个运动周期内4种光谱混叠情况等概率分布，且每种混叠情况的平均像移取最大像移量与最小相移量的平均值，设混叠光谱像元2～9所占平均混叠权重系数分别对应，由图3可得:

其中:

式(2)可改写为:

其中，B(m,n)(v)表示不存在光谱混叠时的核心像素处的光谱信息;B(m,n)′(v)表示存在光谱混叠后的核心像素处的光谱信息:η(p,q)，B(p,q)(v)分别表示核心像素8个邻域像素对应的混叠权重系数和光谱信息，将式(8)代入式(7)，即可获得每个像素存在混叠后的光谱。

4 仿真算法与结果分析

4.1 结合运动模糊模型的TwIST重建算法

在重建过程中，采用了Bioucas-Dias and Figueiredo提出的两步迭代收缩阈值算法(Two-step Iterative Shrinkage/Thresholding , TwIST)[16]。此算法结合了迭代重权值收缩算法(Iterative Re-weighted Shrinkage ,IRS)和迭代收缩阈值算法(Iterative Shrinkage/Thresholding , IST)的优点，既能解决方程严重欠定的问题又具有较好的去噪效果。

TwIST算法解的形式为:

其中:

Γλ表示算子操作;α，β表示权重系数;K表示编码矩阵。

TwIST算法的流程大致可分为以下4部分:

1) 获取参数β，α，和τ;

2) 赋初始值x0;

3) 迭代计算:

4) 当估计值与原图的差值小于某一阈值时，结束迭代，得到的x即为所求。

在本文中，此算法将求数据立方体的解等效为求如下这个非线性无约束最小化问题的最优解的过程:

式中g表示理想图像;f表示真实图像;H表示模糊矩阵;Φ(f)表示正则项;τ表示正则系数。

其中，采用全变分的方法进行正则化，即所谓的TV范数:

对于连续和离散的情况，TV范数分别有如下形式:

4.2 基于CASSI系统的仿真试验设计

整个试验采用MATLAB进行仿真和重建，其流程大致如下: 首先用编码孔径模板卷积选取多光谱原图，然后再用运动模糊算子卷积编码后的图像，以此来模拟实际遥感编码孔径的成像过程; 在获得编码和加入模糊的图像后，采用TwIST算法对图像进行重建，并对重建结果进行了分析和总结。

在仿真过程中，采用了由Brady等人用单分光元件编码孔径快照光谱成像仪所拍摄的光谱数为 33的多光谱图像作为源图像，它对应的成像场景的彩色图像如图4[13]。

在卫星平台参数方面，采用了美国宇航局戈达德航天中心提供的LAND-SAT-4卫星在轨测量的振动数据。卫星平台的振动主要来源于太阳能电池阵列驱动所产生的频率为1Hz、振幅为100μrad的振动，以及卫星反作用轮基波的二次谐波所产生的频率为100Hz、振幅为4μrad和频率为200Hz、振幅为0.6μrad的振动，卫星振动可认为是这3个谐波振动的合成[15]:

图4 成像场景的可见光图像Fig.4 Visible image of the imaging scene

其中φ1,φ2,φ3是随机生成的初始相位角，取值范围为[0, 2π]。在实际仿真中将这3个初始相位角取为[φ1,φ2,φ3]=[0, π/6, π/3]，其周期远大于光谱仪的曝光时间。

在试验中，假设在一个振动周期内，平台抖动引起的光谱混叠等概率分布于4种混叠方式中，而每一块像元的混叠面积取4种混叠情况下，该像元在一个振动周期内可能出现的最大混叠面积与最小混叠面积的平均值的求和，然后与原始像元面积求比值，即得到该像元对原始像元的混叠权重系数，核心处的原始像元混叠权重系数可通过用1减去其它8个邻近像素的混叠权重系数的和得到，最终生成混叠权重系数矩阵，在焦距与像元边长约为10∶1时，混叠权重系数矩阵如图5和图6所示。

图5 俯仰与侧滚综合作用的混叠权重系数矩阵Fig.5 Aliasing coefficient matrix of pitch and roll combined

图6 偏航产生的混叠权重系数矩阵Fig.6 Aliasing coefficient matrix of yaw

比较两个权重系数矩阵，不难发现第二个矩阵元素远小于第一个矩阵，可忽略。因此在实际操作中，仅采用由俯仰和侧滚产生的混叠权重系数矩阵作为图像的模糊算子，用此矩阵对探测器上获得数据立方体进行模糊操作。

在此需强调一点，当平台运动引起的像移超出一个像素范围时，光谱混叠的范围就不再局限于原始像素的8邻域，而会向外扩展到更大的范围，此时需要重新获得相应的更大尺寸的混叠系数矩阵。

4.3 结论分析

实验结果表明，卫星平台的运动会导致严重的光谱混叠，若不针对其进行去模糊操作，复原的图像不能去除光谱混叠，与原图差距很大; 而在重建中加入模糊算子进行去模糊操作后，得到的图像光谱混叠情况去除良好，分辨率得到很大提升。

图7 场景在33个波段上的原始光谱图像Fig.7 33 original spectral images of the scene

图8 加入运动模型后重建出的图像Fig.8 Recovered images applied the motion model

图7为图4场景下的33通道多光谱图像，图8为加入运动模糊后的多光谱图像。对比两幅图，很明显地看出与原图相比，加入运动模型后的图像在483nm, 505～537nm, 637nm和650nm波段处图像强度提升很多，即发生了严重的光谱混叠。

图9为在原图7的基础上加入运动模型，并在图像重建时引入相应模糊算子进行去模糊操作后得到的图像。为了更好地体现去模糊重建的效果，截取了第21波段图像的局部图进行对比说明。

从图10中不难看出，加入运动模糊后的图像在边界处(香蕉与苹果的分界处和苹果与背景的边界处)产生了严重的伪轮廓，并且丢失了很多细节信息(例如苹果的柄在(b)中是看不见的)。经过去模糊重建后，成功消除了边界区域的伪轮廓，并且恢复了如苹果柄这样的细节信息。

图9 进行去模糊重建后的图像Fig.9 Recovered image with deblured

图10 第21谱段局部对比图Fig.10 Local contrast of the 21st spectrum

图11为含去模糊操作的重建结果和不含去模糊操作的重建结果的PSNR(峰值信噪比)对比图。红线表示去模糊过的图像的PSNR值，蓝线表示未去模糊的图像的PSNR值。可见，去模糊的图像PSNR值普遍远高于未去模糊的图像。

图12为原图、含去模糊操作的重建结果和不含去模糊操作的重建图像的平均梯度对比图。可见，含去模糊重建的图像与原图更一致，更接近真实值，未去模糊的图像则相对有很大偏离。

图11 重建结果PSNR对比Fig.11 PSNR of recovered images

图12 重建结果平均梯度对比Fig.12 Average gradient of recovered images

通过以上对比，可以看出在多光谱遥感中，由于平台抖动引起的光谱混叠是不可忽视的，而重建过程中是否包含相应的去模糊操作对复原图像的品质有至关重要的影响，因此对平台抖动进行建模，再根据抖动模型得到相应模糊算子进行去模糊重建对提高多光谱图像品质是十分有意义的。

5 结束语

本文通过对卫星平台的运动建模，获得了对应的模糊算子，进而用此算子对多光谱图像进行模糊处理，来模拟编码孔径遥感成像的过程; 然后再将此模糊算子应用于图像重建过程中，获得了很好的重建效果。除此之外，还进行了对比试验——重建时不加入模糊算子，证明了模糊算子对重建图像像质的影响。通过整个仿真和重建过程，可知卫星遥感平台的抖动对光谱成像过程具有不容忽视的影响，获得准确的运动模型将对后期的图像重建工作十分重要。

压缩编码孔径多光谱成像是一项基于压缩感知(CS)理论发展起来的新技术，其在遥感领域的应用仍有很长的路要走。本文仅对将编码孔径运用到遥感多光谱成像所遇到的第一道难题——严重的光谱混叠提出了初步的解决方案，对于文中所提及的混叠权重系数矩阵与焦距和像元尺寸的关系、重建的优化算法及系统模型优化等问题，仍需进一步研究。

References)

[1] Green R O, Eastwood M L, Sarture C M, et al. Imaging Spectroscopy and the Airborne VisibleInfrared Imaging Spectrometer(AVIRIS) [J]. Remote Sens. Environ, 1998, 65(3): 227–248.

[2] Herrala E, Okkonen J T, Hyvarinen T S, et al. Imaging Spectrometer for Process Industry Applications[J].Proc. SPIE, 1994,2248(33): 33–40.

[3] Morris H, Hoyt C, Treado P.Imaging Spectrometers for Fluorescence and Raman Microscopy: Acousto-optic and Liquid Crystal Tunable Filters[J]. Appl. Spec- trosc, 1994, 48(7): 857–866.

[4] Descour M R, Volin C E, Dereniak E L, et al.Demonstration of a High-speed Nonscanning Imaging Spectrometer[J]. Opt. Lett.,1997, 22(16): 1271–1273.

[5] Brady D J, Gehm M E.Compressive Imaging Spectrometers Using Coded Apertures[J]. Proc. SPIE, 2006, 6246: 62460A.

[6] Candès E, Romberg J, Tao T.Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[J]. IEEE Trans. Inf. Theory, 2006, 52(2): 489–509.

[7] Candès E, Tao T. Near Optimal Signal Recovery From Random Projections: Universal Encoding Strategies[J]. IEEE Trans. Inf.Theory, 2006, 52(12): 5406–5425.

[8] Donoho D. Compressed Sensing[J]. IEEE Trans. Inf. Theory, 2006, 52(4): 1289–1306.

[9] Donoho D L. Compressed Sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289–1306.

[10] Candès E, Romberg J, Tao T. Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information[J]. IEEE Transactions on Information Theory, 2006, 52(2): 489–509.

[11] Candès E. Compressive Sampling[C].In: Proceedings of International Congress of Mathematicians. Zürich, Switzerland: European Mathematical Society Publishing House, 2006: 1433–1452.

[12] 许晓峰.光学压缩编码成像及其复原算法研究[D]. 西安电子科技大学硕士研究生学位论文, 2011.XU Xiaofeng. Optical Compressive Coded Aperture Imaging and Research on the Recovery Algorithm [D] Master thesis of Xidian Univ., 2011. (in Chinese)

[13] Wagadarikar A A, Pitsianis N P, SUN Xiaobai,et al. Spectral Image Estimation for Coded Aperture Snapshot Spectral Imagers[J]. Proc. SPIE, 2008, 7076: 707602.

[14] 唐秋艳, 唐义, 曹玮亮, 等. 卫星平台复杂振动引起的光谱成像退化仿真研究[J]. 物理学报,2012, 61(7): 070202-1-070202-9.TANG Qiuyan, TANG Yi, CAO Weiliang, et al. Simulation of Imaging Spectrometers Degraded by Satellite Vibrations[J].Acta Physica Sinica, 2012, 61(7): 070202-1-070202-9.(in Chinese)

[15] 王晓燕, 唐义, 唐秋艳, 等. 色散型推扫成像光谱仪在轨动态成像建模与分析[J]. 北京理工大学学报, 2012, 32(3):291–295.WANG Xiaoyan, TANG Yi, TANG Qiuyan, et al. Dynamic Imaging Modeling and Analysis of Dispersive Spectrometer with Pushbroom Imaging on Orbit [J] Journal of Beijing Institute of Technology, 2012, 32(3): 291–295. (in Chinese)

[16] Bioucas-Dias J M, Figueiredo M A. A New TwIST: Two-step Iterative Shrinkage /Thresholding Algorithms for Image Restoration[J].IEEE Trans Image Process, 2007, 16(12): 2992–3004.