基于张量分解的光谱图像压缩感知重构 *

赵梓渊 ，唐意东 ，黄树彩

（1. 空军工程大学 信息与导航学院，陕西 西安 710077；2. 中国人民解放军95607 部队，四川 成都 610066；3. 空军工程大学 防空反导学院，陕西 西安 710051）

0 引言

近年来，随着窄带滤光技术和成像光谱技术的快速发展［1-3］，光谱成像已经成为发展天基预警卫星新型探测手段的重要方向［4-6］。光谱图像包含探测场景的空间信息和光谱信息，具有“图谱合一”的特点［7］，提供了丰富详细的物质类别信息，能够大大提高目标检测和分类识别精度，但其庞大的数据量也给数据采集和实时处理带来了巨大挑战。面对这一矛盾，基于压缩感知［8-13］的光谱图像编码，为实现数据的高效采集提供了一条可行途径。现有的压缩采样编码多集中在二维空间域，而光谱维则采用常规的数据压缩方法，忽略了光谱图像的三维空间稀疏性，仍存在资源浪费，其算法存在循环嵌套，时间复杂度亦为较高的平方阶O(n2)。

针对这一问题，本文提出一种基于张量分解的光谱图像压缩感知重构方法，该方法利用光谱图像数据三维空间稀疏性，建立基于三阶张量Tucker 分解的光谱图像重构模型，设计基于正交匹配（orthogonal matching pursuit，OMP）［14］算法的模型求解方法；并将OMP 算法推广到三维空间，设计一种以三阶张量为字典原子的正交匹配追踪算法，在三维空间实现光谱图像数据的压缩采样及解码重构。

1 光谱图像三维空间稀疏性分析

光谱图像数据通常有光谱空间、图像空间和特征空间等3 种表达形式［15］，如图1 所示。

图1 光谱图像数据的表示形式Fig. 1 Representation of spectral image data

图1 中，光谱空间表达为一条连续光谱曲线，包含各像元的光谱信息，反映的是物质辐射特性随波段变化的关系，不同物质具有不同的光谱曲线；图像空间表达能够反映场景中各类物质的空间位置分布和上下文信息；特征空间中，则用多维矢量表示各样本的光谱向量，能够定量反映样本的光谱辐射特性及其变化规律，不同的物质有着不同的统计分布特性。

作为三阶张量，光谱图像包含二维空间数据和一维光谱数据，具有空-谱切片和空-空切片2 种切片分解方式。利用空间域或谱间稀疏性，可以在二维空间域或一维谱域内实现光谱图像压缩采样［16-19］，但实际上光谱图像三维数据块在空谱域三维空间内同样具备稀疏性。以图2a）为例，构造3个不同的离散小波变换基，基于Tucker 分解计算光谱图像的变换域核心张量，其向量形式如图2b）所示。选择其中10%的显著系数重构光谱图像，如图2c）所示，各波段峰值信噪比（peak signal to noise ratio，PSNR）曲线如图2d）所示。可以看出，光谱图像数据在三维空间内具有较强的稀疏性，稀疏分解后，利用部分显著系数能够精确重构原始图像，重构数据各波段的峰值信噪比达到36 dB 左右。可以发现，稀疏分解重构损失了光谱图像的部分高频细节信息，导致重构图像较原始图像更加平滑，不同类型面目标之间的边界有一定程度的模糊。

图2 光谱图像三维空间稀疏性Fig. 2 Sparsity of spectral image in three-dimensional space

2 基于三阶张量分析的压缩感知光谱图像重构

基于光谱图像三维空间稀疏性，直接在三维空间内进行压缩采样和解码重构，可以有效降低计算复杂度，其采样过程如图3 所示。场景辐射通过前光学器件后入射到数字微镜器件（digital mirror device，DMD），空间信息在DMD 上形成一面阵，空间编码图像经过聚光透镜入射到分光计。随后，分光后的空间编码数据入射到随机编码光电传感器阵列，进行光谱维采样编码，最终得到三维采样数据。最后，利用RNG 控制开关1-L的开合进行K次0-1 随机编码，实现光谱维压缩采样编码，此时其光谱维采样率为K/L。

图3 光谱图像三维空间压缩采样Fig. 3 Compressive sensing in three-dimensional space

2.1 TDB-OMP 算法

式中：×n为张量的n-模式积；为变换域系数，即Tucker 分解的核心张量。利用压缩采样矩阵和，分别在光谱图像三维数据的行、列和光谱维进行独立观测：

根据Tucker 分解与Kronecker 积的关系，利用向量化算子vec( · )，将式（2）转化为一维向量形式为

原子更新后，通过求解下列最小二乘问题估计系数向量为

式中：⊙为Khatri-Rao 积［20］；和分别包含了t= 1，2，…，k次循环中从Φ1，Φ2和Φ3中选择的原子，即vec(Y)为向量化后的压缩采样数据。此时，式（6）的解为

式中：( · )#为矩阵的Moore-Penrose 伪逆。式（7）的计算量特别大，利用Khatri-Rao 积性质可以得到

式中：(n= 1，2，3)为第k次循环中新增入Dn的原子。则第k次循环中Z-1的迭代计算式如下：

最终第k次循环结束时，根据式（8）得到系数向量的计算式为

得到系数向量之后，残差向量为r=y-，将其矩阵化可以得到残差矩阵：

式中：matr( · )为矩阵化算子。

当满足中止判据，则中止迭代，否则转入下一次迭代。经过循环迭代得到变换域系数张量后，利用式（12）可以重构光谱图像数据：

综上所述，TDB-OMP 算法通过n-模式积计算更新原子，采用迭代更新的方式计算Zk，有效降低了计算复杂度，其流程如下，其时间复杂度为线性阶O(n)。

TDB-OMP 算法

输入：感知矩阵{Φ1，Φ2，Φ3}，观测值Y，截断误差ε0，稀疏度K0

初始化：

循环：

步骤 1： 计算残差与感知矩阵的n-模式积，搜索最匹配的原子

步骤 2： 更新标签集和入选原子集

步骤 3：求解最小二乘问题，估计稀疏系数

步骤 4： 更新残差

步骤 5： 更新迭代索引k=k+ 1。

End

2.2 3D-OMP 算法

TDB-OMP 算法利用三阶张量与感知矩阵的Tucker 积搜索最佳匹配原子，并采用迭代递推的方式估计稀疏系数，有效降低了算法的计算复杂度。其不足是系数估计和残差更新仍以向量形式进行，字典中的原子仍是一维向量。为保留光谱图像数据的三维结构信息，进一步降低计算复杂度和存储量需求，本节直接利用感知矩阵各原子的外积作为字典原子，将传统OMP 算法推广到三维空间，即3DOMP，基于光谱图像三维空间稀疏性，在三维空间实现光谱图像压缩采样及解码重构。

式中：wi，j，k为系数三阶张量W的元素，定义k=1，2，…，L)为张量原子，其实质是向量和的外积，即。由此，将观测值Y看成张量原子的线性组合，系数张量W的元素为对应原子的权重系数。观测值Y到字典中各原子的投影为

由式（14），（15）可知，在3D-OMP 算法中，字典包含N2L个原子，每个原子是维数为M×M×K的三阶张量。类似于传统OMP 算法，每次迭代中3D-OMP 首先从字典中搜索最匹配的原子，并将其添加到已选原子集，然后求解所有原子集对应的权重系数。在第t次迭代中，观测值Y为已选t个原子的近似线性组合，此时残差R定义为

权重系数更新的目标是搜索系数向量w=，使得残差R的l2范数最小，其实质为求解下列最优化问题：

式中：R：：k，Y：：k和分别为三阶张量R，Y和Γu正面切片矩阵，则式（17）等效为

式中：tr( · )为矩阵的迹，满足，因此有

可得

矩阵H的元素的计算如下：

令Y=a◦b◦c，则的计算式如下：

得到变换域系数张量后，利用式（1）重构原始光谱图像。综上，3D-OMP 时间复杂度为线性阶O(n)。

其算法如下：

任务：

输入：感知矩阵{Φ1，Φ2，Φ3}，观测值Y，截断误差ε0，稀疏度T0

初始化：

循环：

步骤 1： 搜索最匹配原子

步骤 3： 分别利用式（21）和式（22）计算H和f

步骤 4： 求解最优化问题式（19），估计稀疏系数

步骤 6：t=t+ 1

End

3 仿真校验

为检验TDB-OMP 和3D-OMP 2 种解码重构算法的有效性，将它们应用到光谱图像Our Data［22］，考察它们在不同数据压缩率下的重构性能，并与文献［5］中的2D-CRPG+KL 方法进行比较。式（1）中的稀疏基分别为DWT（discrete wavelet transform）变换矩阵Ψ1=Ψ2和DCT（discrete cosine transform）变换矩阵Ψ3。虽然TDB-OMP 和3D-OMP 可以在行、列和光谱维采用任意的测量矩阵组合来保证数据压缩率，但为确保重构结果的可比性，采用与2D-CRPG 算法相同的测量矩阵组合。给定谱间压缩率为50%，即K= 150，。假设空间采样率为Rspa，则空间域采样矩阵为和，其中。

图4 给出了空间采样率为40%～90%时，运用2D-CRPG+KL，TDB-OMP 和3D-OMP 算法重构光谱图像的曲线。从图4 可以看出，在相同的数据压缩率下，TDB-OMP 的重构性能在多数情况下优于2DCRPG+KL，而3D-OMP 算法的重构性能总是优于其他2 种算法，特别是在信噪比较低的波段（波段150～300），其优势更加明显。这主要是因为，2DCRPG+KL 对各波段图像进行独立压缩采样和解码重构，忽略了波段间的相关信息。TDB-OMP 虽然在三维空间实现光谱图像重构，使得重构性能有所提升，但其系数和残差更新仍以向量形式进行，忽略了光谱图像三维数据的结构信息。尽管递推求解方式在一定程度上减小了计算量，但TDB-OMP 需要的存储量仍然较大。而3D-OMP 不仅在三维空间重构光谱图像，而且直接以三阶张量作为字典原子，有效保留了光谱图像的结构信息，进一步提升了重构性能。同时，由于2D-CRPG+KL 方法在进行谱域KL 变换编码时，无论谱域压缩率大还是小都需要计算空间采样数据的协方差矩阵，并进行特征值分解，使得其计算量和需要的存储空间都较大。为保证重构质量，3D-OMP 能够灵活地设置空间采样率和谱域压缩率对光谱图像三维数据进行观测，而对于2D-CRPG+KL 来说，则需要分析空间采样率和谱域压缩率的关系，以及它们对重构性能的影响来优化整体的数据采样率，如何调整空间采样率和谱间压缩率仍然有待解决。

图4 不同采样率下重构光谱图像各波段PSNRFig. 4 PSNRs of reconstructed spectral image with different sampling rate

图5 给出了2D-CRPG+KL，TDB-OMP，3D-OMP 3 种算法重构光谱图像的均峰值信噪比（average PSNR，APSNR）与空间采样率的关系。可以看出，当空间采样率较小时，TDB-OMP 和3D-OMP 的重构性能明显优于2D-CRPG+KL，而随着空间采样率的增长，3 种算法的APSNR 都逐渐增大。其中TDBOMP 和3D-OMP 都呈近似线性增长，但前者的增长速度小于后者，两者之间的差距不断扩大；而2DCRPG+KL 的增长速度较快，当空间采样率较大时，其重构性能超过TDB-OMP，直至赶上3D-OMP。

图5 APSNR 随采样率变化Fig. 5 Change of APSNR over sampling rate

图6，7 分别给出了空间采样率较小（40%）和较大（90%）时，运用3 种算法重构光谱图像的假色图。从图6 可以看出，当空间采样率较小时，2D-CRPG+KL 重构图像存在明显的模糊，而TDB-OMP 重构图像也与原始图像有比较明显的区别，3D-OMP 能够更好地保留原始图像的结构信息，其APSNR 略高于TDB-OMP。从图7 可以看出，当空间采样率为90%时，不同算法重构光谱图像的APSNR 均超过39 dB，且2D-CRPG+KL 的APSNR 值最大。这说明，当采样资源充足，空间采样率接近全采样时，3种重构算法均能较好地保留图像的结构信息和细节特征，三维空间压缩采样的优势并不明显。

图6 空间采样率为40%时不同方法重构光谱图像Fig. 6 Reconstructed spectral images with 0.4 sampling rate

图7 空间采样率为90%时不同方法重构光谱图像Fig. 7 Reconstructed spectral images with 0.9 sampling rate

图8 ，9 给出了空间采样率为40%和90%时，不同算法得到的重构光谱曲线，选取Our Data 中比较典型的3 类样本，其空间位置分别为（80，45），（170，70）和（151，151）。从图9 可以看出，当空间采样率达到90%时，3 种算法均能比较完整地保留典型样本的光谱信息。而从图8 可以看出，当采样率较低时，相比于2D-CRPG+KL 方法，TDB-OMP 和3D-OMP 算法能够更好地保留各类典型样本的光谱信息，以更低的资源消耗提升光谱探测信息的应用精度。

图8 空间采样率为40%时3 种方法的重构光谱Fig. 8 Reconstructed spectrum with 0.4 sampling rate

图9 空间采样率为90%时3 种方法的重构光谱Fig. 9 Reconstructed spectrum with 0.9 sampling rate

4 结束语

本文针对现有光谱图像压缩采样方法仍然存在的高资源消耗问题，利用光谱图像的三维空间稀疏性和数据块结构信息，在三维空间建立基于压缩感知理论的光谱图像压缩采样及重构模型，提出两种基于三阶张量分解的压缩感知光谱图像重构算法，有效降低了重构算法计算复杂度，改善了光谱图像的重构效果，增强了空间采样率和谱间压缩率设置的灵活性。