张建峰 沈 军 张昊平

1(哈尔滨工程大学信息与通信工程学院 黑龙江 哈尔滨 150001)2(上海无线电设备研究所 上海 200090)

0 引 言

光谱成像技术通过获取目标的二维空间信息和一维光谱信息来组成图谱结合的数据立方体[1]，它能够帮助我们识别场景中出现的不同物体，具有广阔的应用前景。传统的光谱成像技术，在探测器一次测量中只能获得数据立方体的一个切面，需要反复测量才能获得目标图像，采样过程不仅花费大量时间，而且测量值占大量存储空间，越来越不能满足现实测量的要求。编码孔径光谱成像技术[2-5]是结合压缩感知理论的新型光谱成像方法，很好地打破了传统光谱成像技术的局限性，该系统利用二维探测器捕捉的是一个光谱图像数据立方体的多通道投影。因此，利用编码孔径光谱成像系统在很大程度上减少了采样数据，提高了图像采集效率。另一方面，编码孔径光谱成像系统采集到的图像是高度压缩的[6]，因此从二维测量值中准确高质量重构出三维的光谱图像立方体成为了一个重要的挑战。为了从压缩测量值中更好地重构出光谱图像，人们做了大量的研究。广泛应用的是将该重构问题转化为l1范数最小化问题[3,7-8]，利用相应的重构算法进行求解。Arguello等[7]利用梯度稀疏投影稀疏重建(GPSR)[9]算法，它使用二维小波变换矩阵与一维余弦变换基的克罗内克积构成稀疏基代入重构框架进行求解。Wagadarikar等[3]则利用两步迭代收缩阈值(TwIST)框架进行求解，它主要将l1范数替换成全变分模型作为正则化项[10]，这种方法能快速地恢复出原光谱图像。文献[11]则采用一种近似消息传递框架(AMP)，引入自适应维纳滤波器重构光谱图像。除了像文献[7]中使用三维解析字典外，采用字典学习的方式提高光谱图像重构质量也成为一种有效的方式。文献[12]分别将解析字典和K-SVD离线训练字典代入重构框架，结果显示采用K-SVD离线训练的字典重构的光谱图像质量得到有效提高。文献[13]则在压缩投影和光谱图像间建立神经网络，通过深度学习的方式重构光谱图像。但是构建大量与目标光谱图像相匹配的训练集十分困难，并且学习过程会耗费大量时间。

上述算法虽然在一定程度上实现光谱图像重建，但都没能充分利用光谱图像在特定冗余字典下的稀疏特性。压缩感知理论指出，信号在冗余字典下表示越稀疏，则重构的质量越好[14-15]。可见构造合适的冗余字典对于图像重构至关重要，因此本文提出了一种自适应分裂Bregman迭代的编码孔径光谱图像重构方法。该算法以分块测量方式构建的图像重构框架为基础[16-17]，很好地将分裂Bregman迭代算法[18]与字典学习算法相结合。在迭代的过程中利用前一次图像的估计值构造训练样本，并结合改进的K-SVD算法自适应学习更新字典，此时得到的字典能对光谱图像实现更有效的稀疏表示，多次迭代之后得到待重构图像。相比之下文献[12]则采用其他光谱图像样本学习得到字典，其字典原子不一定与待重构图像相适应。值得注意的是改进的K-SVD算法，其特点是利用K-SVD算法[19]训练好的字典再次对样本进行稀疏表示，提取出误差较大的原子组成新的矩阵并进行一次SVD运算，将变换后前若干奇异值向量加入到原字典当中，进一步使该冗余字典与待重构图像相适应。仿真实验将本文算法与广泛使用的GPSR、TwIST以及OMP重构光谱图像算法进行了比较，通过视觉比较重构光谱图像的细节以及利用峰值信噪比(PSNR)和结构相似性(SSIM)指标作出对比，证明了本文算法的有效性。

1 相关工作

1.1 编码孔径光谱成像系统数学模型

本文主要以编码孔径快照光谱成像系统(CASSI)模型为基础进行研究。CASSI的成像过程主要包括：首先物镜将目标光源聚焦到编码孔径模板上，实现了空间信息采样；再通过色散元器件对不同谱段的投影进行偏移实现编码；然后利用二维的焦平面阵列有探测器测得混叠的图谱信息；最后根据混叠的信息重构出目标图像。图像测量模型如图1所示。

图1 CASSI成像系统抽象模型

假设F0和F1都为M×N×L数据立方体，分别表示原始光谱图像和通过编码孔径模板后的图像。

F1(x,y,λ)=T(x,y)×F0(x,y,λ)

(1)

式中：x、y表示空间坐标，λ表示光谱值;T表示模板，一般用服从伯努利分布的二值矩阵表示，″0″表示不透光，″1″表示透光。原始数据立方体中的所有高光谱片段都使用同样的空间模板来采样。采用光色散器件将光谱图像各个波段的图像沿色散方向彼此错开实现编码，最终在探测器阵列上成像。成像图案表示如下：

(2)

式中：s(λ)代表光谱沿色散方向的移动函数。式(2)中假设沿x轴色散。以参考波长λ为中心，根据棱镜色散系数决定不同波长光谱侧向移动空间距离大小，参考波长以0度入射时假定无偏移产生。假设目标场景的空间维度为M×N，同时光谱维度为L，则总的像素个数为M×N×L。假设比邻的光谱间侧向色散距离为1，在探测器阵列接收过程中光强沿光谱方向叠加，最后得到一张测量图像，其维度为M×(N+L-1)。为了便于计算，将整个过程表示为：

g=Hf+ω

(3)

式中：g是G的向量形式，长度为M×(N+L-1)；f=vec([f0,f1,…,fL-1]) 是F的向量表示，fk为第k个波段的图像数据；w表示测量过程中的噪声；H为观测矩阵，对于CASSI来说，H是由编码孔径形式和色散元器件共同决定的[16]，表示整个采样和色散过程。

(4)

假设f在冗余字典Ψ下的表示是稀疏的，则可以采用压缩感知计算方法根据测量值重构出图谱信息，表示为：

(5)

f=Ψθ

(6)

式中:θ为某冗余字典Ψ下的表示系数;K为稀疏性约束因子。

1.2 分裂Bregman迭代算法

为了方便后续部分的讨论，在这里概要介绍经典的分裂Bregman算法，它是由Goldstein和Osher提出的用于解决l1范数正则化优化问题。文献[18]证明了该算法的收敛性，在迭代的过程中不断逼近最优解。该算法在提高图像处理质量的同时保持高效的计算效率，在图像处理领域应用广泛[20]。

考虑一个约束优化问题如下：

minu∈RN,ν∈RMf(u)+g(ν) s.t.u=Gν

(7)

式中：G∈RN×M，并且f:RN→R,g:RN→R是凸函数。SBI算法对该问题求解的具体步骤如下：

初始化：设置t、b0、u0、ν0的初始值，参数μ>0。

循环执行：

b(t+1)←b(t)-(u(t+1)-Gν(t+1))

t←t+1

直到满足停止迭代条件。

2 分块测量编码孔径光谱成像数学模型

本文主要采用重叠分块测量的方式构造数学模型。分块的好处是可以很好地降低计算复杂度[16-17]，为之后构造训练样本集提供前提条件。分块即在空间维进行分块测量，光谱维保持不变，进行独立的测量和恢复，最后合成整个目标图像。其过程可以表达为：

(8)

(9)

(10)

x=Dα

(11)

式中：D为稀疏基或冗余字典；α为x在D下的稀疏表示系数。测量矩阵H的构造方式不变，只是矩阵的规模变小。

3 自适应分裂Bregman迭代的图像重构算法

在这一部分主要介绍了对原有数学模型进行变换，利用分裂Bregman迭代求解出相应的图像块；以及每次迭代过程中将求解的图像块构造训练样本集，结合改进的K-SVD算法自适应学习更新字典；多次迭代后逐渐逼近求解出待重构图像。

考虑到式(10)中采用的l0范数是非凸函数且是NP难问题，一般采用l1范数来最佳凸逼近l0范数。则上述模型可以改为如下优化问题：

(12)

通过借鉴分裂Bregman算法将式(12)分裂成两个子问题进行最小化求解，这样做的好处是对于子问题的求解难度要远低于对原问题求解，并且经过反复迭代逐渐逼近最优解。

(13)

利用x(t+1)构造训练样本集r(t+1)=x(t+1)-b(t)，利用改进K-SVD算法进行字典更新，即对下式进行求解：

s.t.‖αk‖≤l,∀k

(14)

在更新后的字典下求解稀疏表示系数为：

(15)

b(t+1)=b(t)-(x(t+1)-D(t+1)α(t+1))

(16)

3.1 x求解问题

式(13)可采用最小二乘法得出最优解。首先对x求偏导得到：

HT(y-Hx)+μ(x-D(t)α(t)-b(t))=0

(17)

进而估计出x(t+1)如下：

(18)

在每次迭代过程中字典D和系数α都会进行更新，使得每次估计值不断接近目标图像，最终求解出目标图像。考虑到HTH不具有对角结构，故求逆(HTH+μI)求逆过程计算代价很大，在每次循环的过程中引入泰勒公式二次项变换为：

(19)

此时将式(19)代入式(17)同样采用最小二乘法解得：

(20)

式中：I为单位矩阵。这样求逆的部分变成对角阵，降低了计算的复杂性。

3.2 字典更新

在成像的过程中冗余字典D的构造至关重要，它直接影响到图像稀疏表示，进而影响图像重构质量。本文在分裂Bregman迭代求解的基础上采用一种改进的K-SVD算法进行自适应训练更新冗余字典。主要是利用迭代过程中待重构图像的近似值构造样本训练集结合K-SVD字典学习算法更新字典；之后再次对样本逐列进行稀疏表示，提取误差较大的原子组成新的矩阵，进行一次SVD运算，将前若干列左奇异值向量加入到字典当中。这样做的好处是丰富了字典原子的特征，能够对原始图像进行更好的稀疏表示。具体实现步骤如下：

1) 初始化：设置阈值ε，训练样本集r,空矩阵w,k=0，利用K-SVD算法得到字典D。

2) 重复执行：

(1) 对rk在字典D下进行稀疏表示，得到系数αk；

(3)k=k+1;

(4) 直到遍历整个样本集。

3) 对w进行一次SVD变换，提取前若干列左奇异值向量加入到字典D当中。

3.3 稀疏表示系数α求解问题

对于稀疏系数的求解问题就是对式(14)进行求解，该问题可以等价于求解为下列问题：

s.t.‖α‖0≤L

(21)

对α的求解就变成了常规的稀疏编码问题，该式可以通过正交匹配追踪算法(OMP)进行求解，其中L表示稀疏系数α的稀疏度。

3.4 算法流程

根据以上分析，给出具体的算法流程如下：

1) 输入：分块取样算子R、测量矩阵H、μ、δ、ε、tol以及分块测量值构成的矩阵y。

2) 初始化：冗余字典D0、x(0)、b(0),本文中D0为随机抽取同类型其他光谱图像块组成的字典；

3) 循环执行下列步骤：

(1) 利用式(9)解出x(t+1)、D(t)、α(t)分别表示上一次更新得到的字典以及相应的系数,x(t)为上一次估计的图像值,y表示原始测量值；

(2) 利用改进的K-SVD算法训练学习得到与原始图像相适应的冗余字典D(t+1)；

(3) 利用OMP算法求解式(21),利用更新后的字典再求获得稀疏表示系数α(t+1)；

(4)t←t+1；

(4) 满足循环终止条件,输出x；

(5) 根据式(9)解出目标图像。

4 实 验

实验通过仿真CASSI系统，以及利用本文方法重构光谱图像，比较其成像效果。实验中采用的原始光谱图像来自于哥伦比亚大学光谱图像数据库和ICVL数据集，如图2所示截取图像大小为256×256×8，包括8个波段，波长范围为440～650 nm。通过经典的CASSI系统仿真得到的测量图如图3所示。光谱图像重构即由测量值恢复出原始各个波段的图像。

图2 原始光谱图像

图3 光谱图像压缩测量值

本文采用分块采样的方式，使用6×6的窗口进行重叠取块，移动步长为4，分块可以有效降低计算复杂度。仿真得到测量值矩阵，然后代入式(12)求解框架中，利用提出的自适应分裂Bregman算法进行求解，重构出目标光谱图像。该框架中的测量矩阵H是由采样模板T以式(4)构成，而仿真中采样模板T主要为服从伯努利分布的随机矩阵，仿照大量文献中分布参数设置为0.5。实验中利用本文算法进行光谱图像重构，设置初始参数μ=τ=0.05，ε=tol=103。最后与文献[7]采用GPSR算法、文献[3]采用的TwIST算法以及压缩感知中常用的OMP算法重构的光谱图像进行比较。比较的方法从主观和客观上考虑，主观评价是图像的成像效果，客观评价首先以常用的峰值信噪比来进行衡量。考虑到光谱图像的特殊性，在计算时分别求出每个波段对应的PSNR值，计算如下式：

(22)

(23)

式中:λ对应第λ个波段的图像，I和K分别表示重构后的图像以及原始图像，MAXIλ表示第λ个波段上图像的最大值。

表1展示了重构的各个波段图像的峰值信噪比。可以看出，利用本文方法重构的光谱图像峰值信噪比明显高于其他算法，而GPSR算法和OMP算法总体上好于TwIST算法。其中：本文算法与GPSR算法和OMP算法求解框架都是以稀疏系数的l1范数作为正则化条件，引入稀疏基或冗余字典；而TwIST是以图像空间信息平滑性作为正则化条件进行约束。可以看出通过利用光谱图像在某一冗余字典下的稀疏性作为约束更加符合解决编码孔径光谱图像重构的问题。本文在分裂Bregman算法迭代求解过程中利用改进的K-SVD自适应训练更新字典，使得冗余字典能够对待重构图像进行更好的稀疏表示，充分利用光谱图像在特定冗余字典下具有稀疏性这一特征，进而提高图像的质量。

表1 重构各个波段图像的峰值信噪比 dB

图4展示了光谱图像压缩测量值在加入不同程度的噪声下，光谱图像的重构效果。为了便于比较，这里对每个波段图像的峰值信噪比进行统计平均值然后比较。可以看出，利用本文方法重构的光谱图像各个波段的峰值信噪比平均值总体在其他三种方法之上，但是在噪声影响较大的情况下图像重构质量提高较小，而当噪声影响较小时，图像重构质量提高较大。

图4 不同程度噪声影响下重构图像质量

图5-图7是利用三组光谱图像进行测试的结果，分别选取其中一个波段的图像并截取一部分进行显示，通过观察图5刻度尺部分、图6建筑物边缘以及图7的路灯部分清晰可见，利用OMP和GPSR重构出的光谱图像在边缘部分不平滑，出现很多相互交错的噪声，总体效果要好于TwIST算法恢复出的图像；相比之下TwIST恢复出的图像虽然很平滑但是图像细节部分丢失严重，整体较为模糊；而利用本文方法恢复出的图像细节更加清晰，而且图像边缘部分也恢复得很好。

(a) OMP算法 (b) GPSR算法

(c) TwIST算法 (d) 本文方法图5 测试图像1

(a) OMP算法 (b) GPSR算法

(c) TwIST算法 (d) 本文方法图6 测试图像2

(a) OMP算法 (b) GPSR算法

(a) TwIST算法 (b) 本文方法图7 测试图像3

利用结构相似性(SSIM)指标进一步衡量图像质量，计算公式为：

(24)

式中：μ表示图像的均值；σ表示图像方差或协方差；c1和c2表示常数。这里的图像均指某一特定波段上的图像。同样为了方便比较，采用平均结构相似性(MSSIM)，即对每个波段图像的SSIM求统计平均。表2展示了上述三幅测试图像重构以后的MSSIM值，MSSIM值越大则表示重构图像与原图像越接近，可以看出本文方法与其他三种算法相比具有一定的优势。

表2 重构光谱图像的MSSIM值

5 结 语

为了解决编码孔径光谱成像技术中，利用传统方法造成光谱图像重构质量不高等问题，本文提出了一种自适应分裂Bregman迭代的图像重构算法。该算法很好地将分裂Bregman迭代算法与字典学习相结合，在每次迭代的过程中自适应训练更新字典。充分利用光谱图像在特定冗余字典下具有稀疏性这一特征，提高图像重构质量。实验通过与传统光谱图像重构过程中广泛使用的GPSR、OMP以及TwIST算法进行比较，结果表明本文算法性能更好。同时，与上述三种方法相比本文算法也存在重构时间长的问题，下一步将进一步改善字典学习的过程，缩短重构时间。