本文将并行压缩感知采样方法应用于计算层析成像光谱仪的光谱立方体投影测量过程, 在传统的CTIS基础上采用空间光调制器对数据立方体色散投影快速压缩采集, 采用低分辨FPA得到高分辨色散图像, 改进CTIS受FPA性能限制的投影信号采集水平, 从而实现CTIS性能改进, 称为并行压缩感知计算层析光谱成像(blocked compressed sensing computed-tomography imaging spectrometer, BCSCTIS)。 该成像系统原理示意图如图4所示, 由传统CTIS成像系统和色散投影压缩感知成像系统组成, 具体包括以下光学元件组: 成像透镜、 二次成像透镜、 反射式DMD A、 二维光栅、 投影信号并行调制DMD B和面阵探测器CCD, 其中DMD A和DMD B分别用于逐波段逐点标定系统矩阵和并行压缩调制色散投影图像。 来自物体的光经成像透镜在DMD A成像并反射到二次成像系统, 二次成像后经二维光栅分光调制成像在DMD B并反射由CCD收集。 在系统矩阵标定过程中, 连续谱激光器作为光源经滤波器选择波段照射DMD A, 选定区域逐个打开微镜并同步记录CCD上的色散投影, 处理得到成像系统的标定矩阵; 在并行压缩调制过程中, DMD A选定区域微镜设置为全开, 调制矩阵控制DMD B对色散图像并行压缩调制由CCD同步记录压缩调制色散投影图像。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 并行压缩感知计算层析成像光谱仪示意图Fig.4 Schematic setup of block compressed sensing computed-tomography imaging spectrometry system

与传统CTIS相似之处是二次成像系统均由物镜、 光阑和准直器组成, 二维光栅位于二次成像系统与光谱立方体投影面之间。 不同之处是本光学系统使用了两个快速空间光调制器数字微镜阵列(digital micro-mirror device, DMD): DMD A位于成像透镜的像面, 用于CTIS成像系统矩阵的逐点精确标定; DMD B位于色散投影成像面, 用于色散投影图像的高分辨率压缩调制。

如1.1节所述, BCSCTIS成像模型包括传统CTIS系统和色散投影压缩感知成像系统两部分。 CTIS成像模型表示将3D数据立方体(x, y, λ)到2D投影图像的过程, 3D数据立方体和2D投影图像分别为f(x, y, λ)和g。 该投影模型表示为系统矩阵H_m和立方体f线性运算

g=Hmf(x, y, λ)+e(1)

式中, f(x, y, λ)∈R^w×h×b, w×h为立方体空间分辨率, b为立方体波段数; H_m∈R^W×H×R^w×h×b, W×H为二维投影图像的空间分辨率; g∈R^W×H; e表示系统噪声。 系统矩阵H_m表示数据立方体的体素扩散函数(voxel spread function, VSF), 实验中使用逐个波长点光源标定建立。 CTIS光谱成像则是从投影采样数据g中采用重建算法恢复出光谱图像。

色散投影压缩感知成像采用并行压缩感知模型, 它将DMD上N×N有效成像区域分成若干块, 每个块的大小为C×C, 并行成像到阵列探测器FPA的NC×NC像素, FPA上获取的对应第i个观察向量列Y⁽ⁱ⁾表示为式(2)

Y(i)=Φ(i)(Xg(i))+E(i)(2)

式(2)中,Xg(i)为以列优先方式表示空间光调制器中投影图像g的第i分块, Φ⁽ⁱ⁾为第i分块子场景图像的投影算子, E⁽ⁱ⁾为采集到第i分块的随机噪声。

设x∈R^B×B(C≤B≤N)为空间光调制器上形成的目标场景组成块中的其中一块, 并且设ϕ, ϕ^*∈R^B×B为分别显示在空间光调制器上的互补测量二值模式, y, y^*∈RBC×BC为FPA采集到的互补压缩测量, Δy∈RBC×BC是互补压缩测量差分矢量, 则在FPA上获得的每个测量值如式(3)所示

y=ϕ⊗x+e1, y*=ϕ*⊗x+e2(3)

式(3)中, e₁和e₂表示互补压缩测量采集到的两次随机噪声。

由式(4)得到一对互补矩阵的测量结果Δy为

Δy=(y-y*)=(ϕ-ϕ*)⊗x+(e1-e2)(4)

式(4)中, ⊗表示元素乘积。

最终, BCSCTIS成像模型表示为式(1)和式(2), 完成对数据立方体二维色散投影的并行压缩测量。 数据立方体重建过程为: 从测量值Y⁽ⁱ⁾通过压缩感知重建算法恢复得到立方体投影g(g={Xg(i)|i=1, …,NC×NC}), 再由g重建数据立方体f(x, y, λ), 见式(5):

f(x, y, λ)=Hm-1{(Φ(i))-1(Y(i))i=1, …, NC×NC}(5)

通过仿真实验验证1.2节BCSCTIS成像模型, 处理流程包括采样数据立方体色散投影、 并行压缩测量色散投影、 重建高分辨色散投影和光谱数据立方体, 对比了直接使用低分辨探测器测量结果和基于并行压缩测量的高分辨率重建两种方法的性能。

仿真实验用数据立方体分辨率为128×128×12(128×128空间分辨率和12个波长), 按相邻波长间2像素色散距离生成436×436高分辨率2D色散投影和190 096×196 608规模的色散投影变换矩阵。 光谱色散投影仿真使用109×109低分辨探测器和4×4分块并行压缩感知采集。 以两种方法分别进行光谱重建。 其中, 低分辨采集光谱成像使用低分辨投影变换矩阵重建, 矩阵大小为118 811×196 608; 并行压缩感知光谱成像使用Hadamard测量矩阵和BCSCTIS系统模型测量和重建。

数据立方体仿真物体选自本·古里安大学跨学科计算视觉实验室(interdisciplinary computational vision laboratory, iCVL)高光谱数据集^[15](BGU iCVL Hyperspectral Image Dataset, plt_0411-1155), 包括430~650 nm间隔20 nm的12个波段, 原始光谱信息如图5所示。 图6(a)为仿真得到的投影图像低分辨率探测器直接测量结果, 分辨率为109×109像素, (b)为在并行压缩感知测量系统中利用同样分辨率探测器得到的测量结果, (c)为通过BCS重建得到的436×436像素高分辨色散投影图像, 采样次数为11次(采样率11/16=68.75%)。 仿真实验光谱重建质量采用数据立方体原图与重建结果的均方误差MSE评价, 计算公式如式(7)

MSE=1m∑i=1m(yi-ŷi)2(7)

式(7)中, y_i是原始图像像素灰度值,ŷi是重建图像像素灰度值, m是图像中的像素数。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 原始光谱信息 (a): 各波长光谱图像; (b): 多波段彩色图像Fig.5 Original spectral information (a): Multi-band image; (b): Multispectral composite colored image

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 (a) 109×109色散投影图像直接测量结果; (b) 109×109压缩测量结果; (c) 436×436重建色散投影图像Fig.6 (a) 109×109 pixels dispersion projection image directly; (b) Compressed sampling image in 109×109 pixels; (c) Reconstructed high-resolution dispersion projection image in 436×436 pixels

图7对比展示了直接测量重建光谱和压缩感知测量重建光谱, 结果显示压缩感知重建光谱MSE远优于直接测量, 其全波段细节更清晰, 说明了BCSCTIS模型的有效性, 验证了低分辨探测器获取高分辨二维色散投影和高质量光谱重建的可行性。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 (a) 低分辨采集重建光谱(低分辨率图像重建光谱); (b) BCSCTIS重建结果(压缩感知重建图像重建光谱); (c) 全波段合成的彩色物体图像(上为低分辨重建合成, 下为压缩感知重建合成)Fig.7 (a) Reconstructed spectral image from low-resolution sampling; (b) Reconstructed spectral image by BCSCTIS; (c) Colored multiple bands spectral image (the image by conventional low-resolution dispersion above, the image by BCSCTIS bellow

光学实验使用连续谱激光器和卤素灯作为光源分别对数字物体和胶片^[16]透射彩色物体进行成像光谱实验, 其中, 数字物体选用中国科学院国家空间科学中心缩写“NSSC”, 物体光谱信息包括“N”: 490 nm、 左侧“S”: 570 nm、 右侧“S”: 530 nm、 “C”: 610 nm, DMD依次显示字母并选择相应波段照射成像合成为多光谱物体; 胶片透射彩色物体印有压缩感知缩写“CS”胶片(“C”部分使用绿色滤光片遮挡, 规格为柯达Filter Wratten Red 25; “S”部分使用红色滤光片遮挡, 规格为柯达Filter Wratten Green 58), 彩色滤光片光谱透过率曲线如图8所示。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 柯达滤光片透过率曲线[16] 红色滤光片: 25; 绿色滤光片: 58Fig.8 Transmittance curve of Kodak color filters No.25 is red filter, No.58 is green filter

系统标定使用连续谱激光器和DMD对430, 450, 470, 490, 510, 530, 550, 570, 590, 610, 630和650 nm共12个波段进行逐点精确标定, 图4所示包含了标定光路, 标定时激光器使用声光可调谐滤波器(acousto-optic tunable filter, AOTF)多通道可调滤波器逐个选波段照射DMD, 逐个翻转DMD中心区域128×128微镜, CCD同步采集每个点的色散投影。 其中, 激光器型号为YSL Supercontinuum Source SC-PRO-7; 滤波器为SuperK Select UV/nIR1; DMD为TI DLP7100, 分辨率为1 024×768, 使用中心有效区域128×128像素, 每个微镜大小为13.68 μm×13.68 μm; CCD为IMPERX B1620, 1 600×1 200像素, 使用中心有效区域950×976像素, 像素尺寸7.4 μm。

光学实验利用了空间光调制器DMD的高空间分辨率、 快速翻转和区域可编程性, 结合连续谱激光器和AOTF多通道可调滤波器进行逐点、 逐波长系统矩阵的精确快速标定, 即完成DMD中心区域128×128像素12个波段数据立方体到950×976像素色散投影图像的标定, 将记录到的逐点色散投影图集合组织为系统矩阵, 维度为(950×976)×(128×128×12), 并以稀疏形式存储。

实验中数据立方体标定共有128×128×12个点, 受探测器帧频限制(实验中IMPERX B1620帧频通常为35 fps), 完成数据立方体中全部体素耗时很长, 严重影响实验和应用效率。 为提高标定效率, 改进系统整体性能, 本文提出并设计实现了目标区域逐点并行标定方法, 采用该方法提高了并行度、 降低了数据加载次数和加载量。 图9为同时标定P1, P2, P3和P4四个点的9个方向投影色散投影并行标定示意图。 标定算法流程如表1所示。 通过该算法, 系统标定时间可以降低到原来的四分之一。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 并行标定示意图Fig.9 Schematic diagram of system matrix calibration in parallel

表1

Table 1

表1(Table 1)

表1 并行标定算法 Table 1 The system matrix calibration parallelly algorithm 已知:标注到光学系统波段i的±1级色散距离为2ω i通过逐波段点光源色散图像测量得到。初始并行标定区域位置的起始坐标为(xi, yi), 区域大小为d×d和每个区域初始标定矩阵为Mi, 进行多个区域同时迭代提取标定数据 1. 标定信号抽取的Maski=δ(ω , ε ), ε 为标定位置对应色散点的松弛量, λi∈{430, 450, 470, 490, 510, 530, 550, 570, 590, 610, 630, 650} 2. for k← 1~d 3. for k← 1~d 4. Wndi← Ψ (Maski, k)确定迭代掩膜滑动窗口位置 5. τ← imread(k)读取对应位置的并行采样帧 6. Mi← Wndi *τ*表示卷积运算 7. end for 8. end for 9. Save Mi as HHi存储并行标定区域的标定矩阵 10. 按照区域列优先从{HHi}以区域列像素数量为步长并行抽取合并为全区域标定矩阵HH 11. norm(HH)←x-min(x)max(x)-min(x), x∈HH采用最大-最小归一化方法对进行归一化处理

表1 并行标定算法 Table 1 The system matrix calibration parallelly algorithm

图10—图12列出多光谱数字物体“NSSC”和透射彩色物体“CS”的低分辨测量和并行压缩测量与成像光谱重建结果对比, 其中低分辨测量采用122×122, 为体现并行压缩感知优势采用更低分辨率61×61进行成像测量, 结果表明并行压缩感知计算层析成像光谱展现出更准确的光谱重建。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 (a) NSSC/CS色散投影图像(122×122); (b) NSSC/CS压缩测量结果(61×61); (c) NSSC/CS重建色散投影图像(976×976)Fig.10 (a) NSSC/CS Low-resolution dispersion projection image in 122×122 pixels; (b) NSSC/CS Compressed sampling in 61×61 pixels; (c) NSSC/CS Reconstructed high-resolution dispersion projection image in 976×976 pixels

图11

Fig.11

Figure Option ViewDownloadNew Window

图11 数字物体光谱重建结果, 波长单位为nm (a): 低分辨重建; (b): 并行压缩感知重建; (c): 全波段合成物体(上为低分辨重建合成, 下为压缩感知重建合成)Fig.11 Reconstructed spectral object of digital chart, nm as unit (a): Reconstructed multispectral image using conventional low-resolution dispersion; (b): Reconstructed multispectral image using BCSCTIS; (c): Multi-spectral composite colored image(the image by conventional low-resolution dispersion above, the image by BCSCTIS bellow)

图12

Fig.12

Figure Option ViewDownloadNew Window

图12 数字物体光谱重建结果, 波长单位为nm (a): 低分辨重建; (b): 并行压缩感知重建; (c): 全波段合成物体(上为低分辨重建合成, 下为压缩感知重建合成)Fig.12 Reconstructed spectral transmitting colored object, nm as unit (a): Reconstructed multispectral image using conventional low-resolution dispersion; (b): Reconstructed multispectral image using BCSCTIS; (c): Multi-spectral composite colored image(the image by conventional low-resolution dispersion above, the image by BCSCTIS bellow

图10中(a)-NSSC/CS分别为探测器采集到的“NSSC”和“CS”物体的色散投影图像, 分辨率为122×122像素; (b)-NSSC/CS为使用DMD按照16×16分块并行压缩感知测量到的“NSSC”和“CS”物体图像, 相应的探测器成像分辨率为61×61, 采样率为0.687 5; (c)-NSSC/CS列为对压缩采样低分辨率图(b)-NSSC/CS重建得到的高分辨色散投影图, 分辨率为976×976。

图11展示了数字物体“NSSC”光谱重建结果对比。 (a) 为122×122低分辨重建结果, 每个光谱图像都出现大量信号干扰: 在450, 470, 510和550 nm光谱图像中出现了不应出现的信号; 490 nm光谱物体“N”、 530 nm光谱物体“S”和570 nm光谱物体“S”较为准确重建, 但细节出现不同程度模糊, 仍不如BCSCTIS重建结果; 610 nm光谱图像则出现了比较严重的干扰。 (b) 为利用61×61像素的低分辨测量结果进行并行压缩感知光谱重建结果, 数字物体的每个字母在对应的波段均得到准确重建: 光谱物体“N”出现在490 nm、 左侧“S”出现在570 nm、 右侧“S”出现在530 nm、 “C”出现在610 nm, 和多光谱物体一致; (c) 为全波段成像结果, BCSCTIS重建彩色多光谱物体细节也更为清晰、 生动、 可辨。

图12展示了透射彩色物体“CS”光谱重建结果对比。 (a) 为122×122低分辨重建结果, 重建效果较差: 470和490 nm光谱中错误出现了“S”信号; 510, 530, 550和570 nm对应光谱呈现比较模糊的物体“C”, 同时还出现了信号“S”串扰, 610 nm处出现强度高于630和650 nm光谱的“S”信号, 这也不符合红色滤光片的透过率曲线特征。 (b) 为61×61的低分辨并行压缩感知光谱重建结果: 绿色字母物体“C”信号明显出现在波长510, 530, 550和570 nm, 对应于图8中绿色滤光片58号光谱透过曲线范围, 红色字母物体“S”信号明显出现在波长630和650 nm, 对应于图10中红色滤光片25号光谱透过曲线范围, 其他波段没有出现干扰信号。 (c) 为全波段成像结果, BCSCTIS重建彩色物体细节也更为准确。

本节通过计算层析成像光学实验结合光谱色散投影并行压缩感知仿真实验验证了并行压缩感知计算层析成像光谱模型的正确性, 实验中BCSCTIS重建出更高分辨率的色散投影图像和更为准确的光谱立方体。 实验验证了依据压缩感知信号采样理论采用调制测量方式, 可以改善传统计算层析成像光谱系统中探测器分辨率对重建质量的制约问题。