基于多通道空间光谱全变差的衍射光谱图像复原算法

王 旭 陈 强 孙权森

(南京理工大学计算机科学与工程学院 南京 210094)

Fig. 1 The framework of the proposed MSSTV model图1 MSSTV模型框架

光谱成像技术[1]因其同时具备成像与光谱探测的优点，现已广泛应用于环境监测、军事侦察、医疗诊断、农林业监控、矿物及气体探测等诸多领域.成像光谱仪种类繁多，但是传统的成像光谱仪往往存在构造复杂、光通量小、体积较大以及很难凝视成像等缺点.1995年，美国光量子中心罗姆实验室的Lyons[2]首次提出将衍射光学元件与传统CCD(charge-coupled device)相机结合，为成像系统增加了光谱功能.随后，美国太平洋高技术公司研制了第1台商用成像多光谱检测系统[3]，用于军事目标探测识别.1999年， Hinnrichs[4]研发出了用于油气探测的Sherlock系列中波红外成像光谱仪.2009年，美国陆军研究实验室研制出了工作波段为长波红外的衍射透镜成像光谱仪[5].国内起步较晚，中国科学院长春光机所的于斌[6]、刘英[7]以及浙江大学的梁靖宇等人[8]先后展开了理论研究，并成功研制出原理样机和红外衍射镜头.

由于衍射透镜轴向色散的特性，准焦波长在其焦面所成的准焦像会与相邻波长在该位置处所成的离焦像叠加，同时采集过程中会带入噪声，使得衍射光谱图像严重退化，影响后续应用.因此需要利用图像复原技术对采集的光谱图像进行重构，以提高成像结果的空间分辨率和光谱分辨率.目前，针对衍射光谱图像的三维去卷积问题，国内外学者已提出了最近邻法[2]、逆滤波法[2]、JVC(jan-van cittert)迭代法[9]、Tikhonov正则化方法[9]等图像复原算法.其中大多数方法的抗噪能力较弱，最近邻方法实现简单且能够快速求解，但是复原精度不高；逆滤波法在复原过程中会放大原有噪声或引入额外的噪声，改进后的逆滤波法虽然引入了能够抑制噪声的因子，但需要知道观测图像的信噪比；JVC方法复原效果一般，迭代时间长，求解速度受限且可能引入严重的振铃现象；Tikhonov正则化方法的结果有所提升，能够抑制噪声，但是会导致结果过度光滑，对图像边缘和纹理的处理效果不够理想.此外，这些方法都会产生不同程度的光谱失真现象.

综上，衍射光谱图像重构面临3个问题：1)在去模糊过程中受到噪声的干扰，使得该反卷积问题的求解不稳定，易引入额外的噪声；2)现有的算法只运用了光谱图像的空间信息，忽略了其丰富的光谱信息，复原效果一般，且容易产生光谱失真现象；3)大多数方法耗时长，求解速度受限.基于空间信息的全变差(total variation, TV)正则化[10]可用于此类重构问题，能够较好地保持图像边缘等信息，但是容易平滑掉图像纹理和弱边缘等细节信息从而产生阶梯效应，且没有考虑到光谱先验信息，无法缓解光谱失真问题.

为了解决上述3个问题，本文提出一种基于多通道空间光谱(简称空谱)全变差(multichannel spectral-spatial total variation, MSSTV)的重构模型，充分利用衍射光谱图像的空间信息和光谱信息来抑制噪声和保留图像边缘，使用交替方向乘子法(alternating direction method of multiplier, ADMM)全局一致性框架[11]，通过多次迭代得到清晰图像并加快求解速度，如图1所示.本文还进行了相关的仿真模拟和图像重构实验，采用平均峰值信噪比(mean peak signal-to-noise ratio, MPSNR)、平均结构相似度(mean structural similarity index metric, MSSIM)、平均光谱角距离(mean spectral angle distance, MSAD)等指标评价复原图像的质量[12].

1 相关工作

1.1 衍射光谱成像原理

衍射透镜成像光谱仪由衍射光学元件(diffractive optical element, DOE)、成像镜头和传感器组件、距离控制装置组成，其中衍射光学元件具备特有的轴向色散特点，入射波长与其焦距成反比.同一物体的不同波段图像将沿光轴方向分层成像，可使用传感器沿着光轴方向对固定波段范围内的成像依次扫描，获得光谱数据立方，从而同时完成色散与成像的功能.如图2所示，将探测器移至波长λ0对应的准焦面A处接收其准焦像，此外该传感器还将接收到其他相邻波段在该位置所成的离焦像.移动传感器采集指定波段范围内的衍射光谱图像数据，接着对数据进行三维反卷积处理，从而复原出不同波段各自清晰的光谱图像.

Fig. 2 The imaging principle diagram of DOE图2 衍射光学元件的成像原理图

不失一般性，本文将衍射透镜成像光谱仪的工作波段按照间距Δλ分为m个波段，并假设在Δλ范围内，对应波段的点扩散函数(point spread func-tion, PSF)固定不变.在波段λi准焦位置处采集到的二维模糊图像可以表示成：

(1)

其中，uj表示波段λj的准焦像，h|i-j|表示波段λj在波段λi准焦位置处的PSF，即离焦量为|i-j|Δλ的PSF，ni表示系统采集图像时产生的噪声，符号“*”表示卷积操作，这里假设噪声为加性噪声，待采集图像gi的大小为n×n.

分别将gi，uj，ni拉伸成一维列向量Gi，Uj，Ni∈Rn2，式(1)可改写成矩阵-向量的形式：

(2)

其中，Hij∈Rn2×n2是点扩散函数卷积核h|i-j|的矩阵表示.基于上述定义，m个波段的成像过程可以集合成一个矩阵运算公式：

G=Hu+N，

(3)

其中

式(3)中，G∈Rmn2表示观测向量，u∈Rmn2表示原始图像，N∈Rmn2表示加性噪声，H∈Rmn2×mn2表示模糊算子.

本文的主要任务是在点扩散函数集合H和观测数据G已知的情况下恢复出干净的衍射光谱图像u.在忽略噪声影响的情况下，只需要对模糊算子H求逆，利用逆滤波算法反向求出原始图像u.然而由于该图像重构过程是一个不适定的反卷积问题，一旦引入噪声，通过逆滤波算法计算得到的解u将是不稳定的,并且对噪声极其敏感.为了保证u的稳定性，需要额外利用一些先验信息，最常见的方法便是在u的重构模型中引入正则项，写成：

(4)

其中等号右边第1项Φfid(u,G)是数据保真项，用来保证理想清晰图像u与观测图像G之间的关系，使得模型的解不会过度失真.对于高斯噪声，数据保真项通常表示为

(5)

其对应u的最大似然估计.式(4)等号右边第2项Φreg(u)为正则项，使用u的先验信息进行建模来约束图像u.正则化参数μ用来平衡数据保真项和正则项之间的权重.

1.2 TV模型

1992年，Rudin等人[13]首次提出全变差正则化模型并用于处理图像去噪问题，后又广泛用于各类图像恢复问题. Tikhonov正则项和TV模型都能够改善不适定问题的病态性，但Tikhonov正则项处理图像边缘时会产生模糊效果，而TV模型却可以有效地保留图像的边缘信息，TV模型为

(6)

式(6)中，Dx，Dy分别表示一阶水平和垂直梯度算子.

对于高光谱图像，传统方法一般是按照式(6)对光谱图像进行逐波段复原操作，波段间相互独立，该逐波段全变差正则项[14]可以表示为

(7)

其中，(Dxu)i,k，(Dyu)i,k分别表示第k波段第i个像素处的水平和垂直差分.

此外，Yang等人[10]提出了MTV(multichannel extension of total variation)模型，该模型对多通道图像逐像素从垂直和水平方向计算全变差，公式为

(8)

MTV模型能够自动调整各波段的去噪强度，噪声强度越强的波段将被平滑得越多.此外，MTV模型还给出了精确最小化TV范数的解决方案，有利于加速求解过程.但是，MTV模型和HTV(hyper-spectral image total variation)模型一样，只考虑了图像的空间信息，在处理光谱图像时，没有施加对光谱域的局部平滑约束，无法缓解光谱失真现象.

1.3 PSF模型

本文采用高斯PSF模型进行仿真，该模型使用二维高斯函数模拟模糊效果，定义为

(9)

其中，高斯分布的方差σ与变量z呈线性变化关系.

2 基于多通道空间光谱全变差的重构模型

2.1 模型的提出

为了充分运用衍射光谱图像的空间信息和光谱信息，本文在MTV模型的基础上加入光谱先验，提出了一个针对多通道的空谱全变差模型，表示为

(10)

在式(4)的基础上，选择式(5)作为数据保真项，并引入MSSTV正则先验，使得衍射光谱图像重构的结果无论在空间维度还是在光谱维度上都得到改善，具体模型为

(11)

将Pi代入到式(11)中，可得:

(12)

2.2 MSSTV模型求解与优化

由于衍射光谱图像往往是高维的且2.1节所提出的多通道空谱TV模型存在不可微性，快速有效地求解出式(13)的最优值将会比较困难.为了解决这个问题，本文采用交替方向乘子法将复杂的优化问题分割成多个容易求解的子问题[15-16].首先，引入n2个辅助变量wi=Piu，i=1,2,…,n2，式(12)就可以转化成如下的等式约束问题形式：

(13)

然后构建最小化函数(式(13))的增广拉格朗日函数L(u,wi,ρi)，将上述问题转变成无约束优化问题：

(14)

缩放增广拉格朗日函数中的参数ρi，该增广拉格朗日函数可以写成

(15)

接下来，每次仅优化目标函数中的一个变量，固定其他所有变量，然后交替迭代更新每一个待求解变量[17].下面将详细说明每个变量的求解步骤.

步骤1. 固定u和所有的bi，更新wi，只将关于wi的项提出来，因为所有wi之间不会相互影响，可根据ADMM的全局一致性框架，将关于wi的原问题分解成以下n2个互相独立的最小化问题：

(16)

i=1,2,…,n2.

式(16)能够通过二维收缩算法[10]求解：

(17)

i=1,2,…,n2.

步骤2. 固定所有的wi，bi，只保留与u相关的项，优化u，可得到目标函数：

(18)

观察可知，目标函数式(18)只针对变量u的最小化过程其实是个最小二乘问题，其对应的等式可以表示成

(19)

2.1节中提到的Dx，Dy，Dz∈Rmn2×mn2分别表示针对原始图像u的一阶水平、垂直、光谱方向梯度算子的矩阵形式.令P(j)∈Rn2×mn2等于由P1,P2,…,Pn2的第j行按行组成的矩阵，P∈R3mn2×mn2等于由P(1),P(2),…,P(3m)按行组成的矩阵，P(j)(j=1,2,…,3m)分别表示针对每个波段的水平、垂直和光谱方向梯度算子.通过推导可得：

(20)

(21)

步骤3. 固定wi，u，更新bi：

bi=bi+α(Piu-wi),i=1,2,…,n2

,

(22)

其中α为bi的更新步长.

综上所述，我们可以将基于MSSTV模型的衍射光谱图像重构过程总结为算法1中的伪代码：

算法1.使用MSSTV模型的衍射光谱图像重构算法.

输入：点扩散函数矩阵H，观测图像G，拉格朗日乘子bi，模型参数μ>0，β>0，α>0，最大迭代次数Nmax;

输出：重构所得的图像O.

① 初始化：μ=0.125，β=2-6，α=1.618，bi=0，threshold=10-3，Nmax=100;

④ 更新uk+1：

⑥ end while

⑦ 输出O=uk+1.

2.3 时间复杂度分析

根据1.1节定义，所用衍射光谱数据集的波段数为m，每个波段的大小为n×n.步骤1涉及二维收缩算法，其时间复杂度为O(mn2).步骤2对于u的计算主要在于快速傅里叶变换和高斯消去运算，复杂度为O(max(m3n2,mn2log(n2))).步骤3对拉格朗日乘子bi的更新是个线性运算，时间复杂度为O(mn2).综合上述3个步骤的求解复杂度，算法1每一轮迭代的时间复杂度为O(max(m3n2,mn2log(n2))).

3 实验结果及分析

3.1 实验设定

为了验证本文算法的有效性，分别在2组多光谱图像数据上进行实验.第1组是衍射光谱图像重构实验常用的飞机场数据集，包含了410～850 nm共23个波段的多光谱图像数据，每个波段空间大小为256×256，光谱分辨率为20 nm.第2组采用机载多光谱扫描仪数据集Washington DC Mall的部分数据，本文选取了该数据集波段30到波段40的子集进行实验，尺寸为256×256×11.在图像重构前，先根据衍射光谱成像原理模拟衍射透镜成像光谱仪的观测数据，然后加入不同程度的高斯白噪声，仿真所用的点扩散函数模型选取与实测数据模糊程度相当的高斯PSF模型.本文选取改进逆滤波法(IIF)[18]、Tikhonov正则化法[9]、MTV方法[10]作为对比方法，并采用MPSNR，MSSIM，MSAD这3个光谱图像复原质量评价指标作为衡量标准.

3.2 结果分析

根据衍射成像模型(式(1))，利用飞机场数据合成模糊图并在每个通道中添加高斯白噪声，每个通道的噪声方差相同，且均值都为0.选取方差为3和方差为23两种情况，比较本文方法和其他方法在不同噪声程度下的复原效果.为了直观地展示本文方法和其他算法在衍射光谱图像重构中的表现，图3和图4分别给出了噪声方差为23时波段510 nm和噪声方差为3时波段810 nm的对比结果.图3和图4的第2行展示了第1行中框出内容的放大图像.尽管改进逆滤波法和Tikhonov正则化法去除了大部分模糊，但是局部仍有残留且噪声较大，见图3(c)和图3(d).比较图3(e)和图3(f)，可以发现本文的MSSTV模型使用光谱一致性有助于保留更多的细节信息，且对噪声的抑制作用更加明显.

Fig. 3 The reconstructed results of simulated Airport dataset(band 510 nm)图3 飞机场仿真数据集复原结果(波段510 nm)

Fig. 4 The reconstructed results of simulated Airport dataset(band 810 nm)图4 飞机场仿真数据集复原结果(波段810 nm)

在Washington DC Mall数据集的每个通道中加入均值为0、方差分别为11和1的高斯白噪声.图5展示噪声方差为11时波段32的复原结果.在光谱分辨率较小且光谱响应普遍较低时，4种方法的重构效果都比较好，由于本文的MSSTV方法同时利用空间相关信息和光谱局部光滑性约束，从图5(f)中可以看出该方法的结果在抑制噪声方面更具有优势，且能更好地保留图像边缘和去除模糊.

Fig. 5 The reconstructed results of simulated Washington DC Mall dataset(band 32)图5 Washington DC Mall仿真数据集复原结果(波段32)

Table 1 The Comparison of Evaluation Indexes of Reconstructed Results with Different Algorithms
表1 不同方法复原结果的评价指标比较

DatasetVarianceIndexDegradedIIFTikhonovMTVMSSTVMPSNR16.0427.2327.2428.5430.41Airport23MSSIM0.64620.78900.78920.91010.9434MSAD0.98410.99220.99230.99250.9936MPSNR15.9629.2929.5730.1131.46Airport3MSSIM0.65010.95030.95090.97490.9758MSAD0.98410.99230.99260.99270.9934MPSNR14.3235.5135.5235.8136.75Washington DC Mall11MSSIM0.58060.95850.95850.96600.9745MSAD0.99110.99380.99390.99530.9964MPSNR14.3941.8042.7042.8043.01Washington DC Mall1MSSIM0.58700.99430.99480.99490.9951MSAD0.99130.99640.99700.99680.9972

定量分析如表1所示，表1中列出了2种数据集在不同噪声程度下的退化图像以及4种方法复原结果的定量评估，所用的3种评价指标均是数值越高，重构效果越好.显然，在2组数据集上不同噪声程度下MSSTV方法的MPSNR，MSSIM，MSAD与其他方法相比均保持最高.当噪声方差为23时，本文方法的MSSIM指标比改进逆滤波和Tikhonov正则化提升了0.15左右，和MTV相比也略有提高，说明本文模型采用的TV正则项在提升复原效果的同时有效地保留了图像结构信息.此外，本文方法的MSAD指标始终高于复原效果较好的MTV方法，这归功于在光谱方向施加了较强的光谱平滑约束，相邻波段提供额外的互补信息从而避免某一波段的单独退化，但是当光谱一致性得不到保证时，可能会引入不必要的噪声.

Fig. 6 Special pixels in the Airport image图6 机场图像中指定像素

对于光谱维度恢复好坏的评价除了使用MSAD指标外，我们还绘制了指定像素的光谱曲线.如图6所示，选择飞机场数据集噪声方差为23时A,B处的像素进行分析，具体比较见图7和图8，横轴表示波段，纵轴表示光谱反射率并归一化至[0,1].3种对比方法在中间波段区域均发生了不同程度的光谱

Fig. 7 The reflectance of special pixel A(228,180) in Airport dataset图7 机场数据集指定像素A(228,180)的反射率

失真，光谱曲线呈锯齿状，而MSSTV方法的光谱曲线在这段区域则比较光滑.这是因为前3种方法只使用了衍射光谱图像的空间信息，而MSSTV利用了2个局部空间平滑约束和1个局部谱间平滑约束，可以在去除模糊和噪声的同时尽可能缓解锯齿状光谱失真.此外，由于该重构问题自身的病态性和边缘波段缺乏足够相邻波段的辅助信息，前几个波段的复原效果都表现不佳.

为了对本文算法和其他算法的时间复杂度做定性比较，本文测试了各个算法在相同环境下的运行时间.所有算法均在配备Intel i5 CPU和8 GB内存的笔记本上运行，软件环境均为Windows10系统下的Matlab2017a.表2显示了飞机场数据集在高斯白噪声方差为23的情况下各算法运行10次并取平均值得到的运行时间.从表2中可以看出MSSTV的运行速度远超于IIF方法，但稍微落后于MTV.此外，由于Tikhonov能够闭式求解，而MSSTV需要多次迭代，因此MSSTV此类迭代方法与Tikhonov方法的运行时间会相差较多.总体而言，MSSTV能够在保证求解速度的情况下有效提高复原效果.

Table 2 The Comparison of Running Time of Different Methods表2 不同方法的运行时间比较 s

3.3 收敛性分析

为了更深入地分析本文方法的收敛性，我们给出飞机场数据在噪声方差为23的情况下每次迭代结果的MPSNR和MSSIM值来表示算法的收敛速度，如图9所示，MSSTV方法在不到10次迭代的情况下就逐渐达到收敛状态.

Fig. 9 MPSNR and MSSIM values versus the iteration number of MSSTV algorithm in simulated Airport dataset图9 在飞机场仿真数据集上MSSTV算法每次迭代对应的MPSNR和MSSIM值

Fig. 10 NSDE versus the iteration number of MSSTV algorithm in simulated Airport dataset图10 在飞机场仿真数据集上MSSTV算法每次迭代对应的NSDE值

图10展示了在飞机场数据集上选取噪声方差为23时MSSTV算法每次迭代对应的NSDE值变化，在7次迭代之后NSDE值趋于收敛.

3.4 参数敏感性分析

最后，本文还讨论了在大噪声情况下正则化参数μ和惩罚参数β的选择问题，图11展示了飞机场数据集在噪声方差为23的情况下参数μ和β与评价指标MPSNR和MSSIM之间的关系，μ选择从2-7～26，相邻参数间隔2倍，同理β选择从2-10～2-2.从图11可以看出，参数μ对本文算法性能的影响很大，当μ的取值在0.05～1之间会得到相对好的结果.另外，该算法对参数β有较好的鲁棒性，当β的取值在2-10～2-4之间都能产生很好的复原效果.

Fig. 11 Sensitivity analysis of regularization parameters μ, β图11 正则化参数μ和β的敏感度分析

4 结 论

本文利用局部空间和谱间平滑性约束，提出了一种基于多通道空谱全变差的衍射光谱图像重构模型MSSTV.该模型采用空间全变差抑制噪声并尽可能保留图像边缘信息，同时利用谱间全变差施加光谱一致性约束，起到保持光谱局部光滑和进一步抑制噪声的作用，并将两者集成到统一的模型中.在此基础上，通过ADMM全局一致性框架将问题划分为多个更容易求解的子问题，有效地解决最优解求值问题，并与快速傅里叶变换算法结合提升框架计算速度.本文选取了2组数据集和2种噪声情况进行仿真实验，与多种最先进的衍射图像复原方法比较，实验结果表明，本文方法可以有效地去除模糊和噪声，同时尽可能地保留边缘等细节信息.此外，MSSTV还可以有效地抑制谱间噪声并减缓锯齿状光谱失真.

在未来的工作里，MSSTV模型还有很大的改进空间，可以考虑参考更多的光谱信息对求解模型进行约束，比如根据衍射光谱图像的相关性使用全局低秩约束.此外，还可以考虑根据空间结构、谱间距离以及噪声强度自动计算空谱权重，调节光谱一致性约束的强度，进一步提高MSSTV的鲁棒性和稳定性.