低秩聚类被动压缩鬼成像

低采样率下的高质量鬼成像（GI）对于科学研究和实际应用具有重要意义，为了在低采样率条件下重建高质量图像，提出了一种高质量的被动式压缩鬼成像重构算法（PCGI-LRC）。基于图像的非局域相似块堆叠而成的矩阵具有低秩和稀疏奇异值的假设，从理论和实验上证明了一种对最小二乘问题与非局域相似块低秩近似问题进行联合迭代求解的方法，能够在低采样率（6.25%~50%）条件下实现高质量鬼成像。实验结果表明：与基于稀疏基约束的GI（GI-SBC）和基于全变分约束的GI（GI-TVC）相比，PCGI-LRC在峰值信噪比、结构相似性系数和视觉观测等方面均更优，在抑制重构噪声的同时保持了目标的细节信息，其中PSNR提升效果优于1.1 dB，SSIM提升效果优于0.04。

Abstract

High-quality ghost imaging （GI） at low sampling rate is of great importance for scientific research and practical applications. Therefore， the reconstruction of high-quality images under low sampling rate conditions remains the focus of GI research. In this paper， a high-quality passive compressive ghost image reconstruction algorithm was proposed， called PCGI-LRC. Based on the assumption that the matrices stacked with nonlocal similar blocks of an image have low-rank and sparse singular values， a joint iterative solution of the least squares problem was demonstrated theoretically and experimentally and the low-rank approximation problem of the nonlocal similar blocks can achieve high-quality ghost images under low sampling rate conditions （6.25%-50%）. Moreover， the experimental results show that the proposed algorithm outperforms the GI based on sparse basis constraints （GI-SBC） and GI based on full variational constraints （GI-TVC） algorithms regarding peak signal-to-noise ratio （PSNR）， structural similarity coefficient （SSIM）， and visual observation. Information of the target is preserved while the reconstruction noise is suppressed； the PSNR is improved by more than 1.1 dB and the SSIM improvement is higher than 0.04 dB.

1 引言

鬼成像，又称为关联成像，是一种不同于传统成像方式的新型成像方式，是近年计算光学成像领域的热点之一。鬼成像最早可追溯到著名的HBT强度干涉实验。此后，1995年Pittman、Abouraddy等^［1-2］利用参量下转换产生的动量纠缠光源，通过双光子计数分别实现了“鬼衍射”和“鬼成像”。2002年，Bennink^［3］等用具有互关联的经典光源在相同的双臂光路架构中实现了鬼成像。2006年，Scarcelli^［4］等将激光照射在旋转的毛玻璃上产生赝热光源完成了鬼成像实验。2008年，Shapiro^［5］提出了仅包含物体光路的单光路鬼成像，利用空间光调制器极大地简化了实验系统，并将它命名为计算鬼成像。这类主动式鬼成像通常是利用数字微镜器件（Digital Micromirror Device， DMD）、空间光调制器（Spatial Light Modulator， SLM）等光场调制器件对光源进行调制，通过增加时域上在不同编码光场下的采样数来增加系统的时间自由度数，将场景回波信号与赝热光场进行关联，在时间维上仅通过桶探测器接收目标的回波信号，从而获得目标的图像信息。为了进一步提高鬼成像的实用价值，需要发展被动式鬼成像，与主动式鬼成像不同，自然光和物体自发辐射等被动光源的相干时间在飞秒量级，而测量光场的时空涨落需要探测器响应时间小于光场的相干时间，并且对自然光这类被动光源的光场时空涨落进行实时控制非常困难。2015年，韩申生^［6］提出了一种被动式鬼成像系统，即基于稀疏约束的鬼成像（Ghost Imaging via Sparsity Constraints， GISC）光谱相机。该相机利用空间随机相位调制器实现了光场的相位调制及色散，同时利用随机相位调制器的极近场衍射，将真热光转换成赝热光，解决了真热光相干时间短的测量难题。另外，GISC光谱相机通过将主动式鬼成像中测量成像视场中各像素所对应的光场在时域上的随机涨落转换为空域上的随机涨落，并进行预先测量，从而无需分束器分光就能实时探测光场随时间的无规涨落。

鬼成像相机通过对场景图像进行压缩编码，随后利用重建算法解译编码光场的信息，这一成像方式突破了传统被动式光电系统采样的限制，提高了光学成像系统的信道容量，实现了超分辨率成像^［6］。为了推进鬼成像的实用化进程，必须实现欠采样条件下的高质量被动式鬼成像。重建算法是鬼成像技术中最重要的环节之一，目前常用的鬼成像算法都是基于压缩感知的图像重构算法。压缩感知是一类通过求解欠定逆问题实现信号恢复的方法，鬼成像通常是求解以下两类问题^［7］：（1）求解 $(P_{1})$ 的最小 $l_{1}$ 范数约束优化问题，常使用贪婪迭代优化算法，如OMP算法^［8］、ROMP算法^［9］。贪婪迭代类算法的计算速度较快，但是需要较多的采样次数，对小尺度信号的重建效果较好，但对大尺度信号的重建效果并不理想^［7-8］。为了解决这一问题，BP算法^［10］、FIST算法^［11］等凸优化算法被提出，用于对大尺度稀疏信号的有效重建，这些算法对噪声有较高的稳定性。（2）求解 $(P_{2})$ 的最小化TV或TC问题， $l_{1}$ -Magic^［12-13］，TwIST^{［11，14］}，NESTA^［15］，GAP^［16］以及ADMM^［17］经常被使用。人们还提出了TVAL3算法^［18］，该算法利用图像局部平滑的先验信息，图像重构过程中在保护图像边缘结构信息的同时抑制噪声，但TVAL3忽略了因像素间差异导致的图像梯度差异问题，抗噪性能减弱，并且容易在图像平滑区域产生阶梯效应。为了进一步改进TVAL3算法，2021年Wang等^［19］提出了GISCNL重建算法，该算法在TV约束的基础上考虑了图像的非局域自相似性，并将它作为约束项加入到重建算法中。经过数值仿真论证，GISCNL相比于TVAL3算法有更强的稳定性。为方便论述，本文将求解 $(P_{1})$ 问题的鬼成像称为基于稀疏基约束的鬼成像（Ghost Imaging based on Sparse Basis Constraints， GI-SBC）^［20-22］，将求解 $(P_{2})$ 问题的鬼成像称为基于全变分约束的鬼成像（Ghost Imaging based on Total Variation Constraints， GI-TVC）^［23-25］。

近年来，为了进一步提高鬼成像算法的重建性能，2018年Huang等^［26］提出了一种基于Landweber正则化与引导滤波联合迭代的压缩鬼成像算法，该算法有效降低了欠采样噪声，提高了分辨率。这表明通过分解预重构步骤和去噪步骤来重建原始目标。2019年，Wu等^［27］在鬼成像模型中引入了低秩约束，实验与仿真均表明：对于拥有较强低秩特性的图像，低秩约束下的鬼成像比稀疏约束下的鬼成像有明显的图像质量改善，特别是在欠采样情况下。基于以上两种算法的特点，本文搭建了鬼成像相机实验平台，提出了一种基于低秩聚类的鬼成像算法，并应用于被动式鬼成像系统中，提高了欠采样条件下被动式鬼成像系统的图像重构质量。该算法将图像重构分解为预重构与低秩聚类降噪两部分。在预重构步骤中，利用共轭梯度法求解最小二乘问题以提取更多细节，对图像进行快速重构；在低秩聚类降噪步骤中，将非局域相似的图像块拉成一维向量，并聚合为二维矩阵，利用广义软阈值算法对该矩阵进行低秩矩阵近似，达到抑制重构噪声的目的。在图像重构过程中，算法交替执行预重构步骤与低秩聚类降噪步骤。仿真与实验结果表明，本文提出的算法能够在有限次的联合迭代步骤中提取出散斑场中的场景信息，实现高质量的鬼成像。

2 原理

2.1　成像模型

本文搭建的鬼成像相机原理示意图如图1（a）所示。自然光源照明目标，前置成像镜头将目标的反射光成像至第一成像面上。将第一成像面划分为 $N$ 个小面元，不同位置的小面元出射的宽带热光在自由空间中传播，随后由空间随机相位调制器进行相位调制，并在其后表面附近进行强度叠加形成散斑干涉场。最后，由显微成像镜头将该散斑干涉场成像至像元数为 $M (M ≪ N)$ 的CCD上^［28］。在该鬼成像相机中，选择散射片作为空间随机相位调制器。这一成像过程可以用如下矩阵描述^［28］：

(\begin{matrix} y_{1} \\ y_{2} \\ ⋮ \\ y_{M} \end{matrix}) = (\begin{matrix} φ_{11} & φ_{21} & \dots & φ_{i 1} & \dots & φ_{N 1} \\ φ_{12} & φ_{22} & \dots & φ_{i 2} & \dots & φ_{N 2} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋱ & ⋮ \\ φ_{1 M} & φ_{2 M} & \dots & φ_{i M} & \dots & φ_{N M} \end{matrix}) (\begin{matrix} x_{1} \\ x_{2} \\ ⋮ \\ x_{N} \end{matrix})

图 1. 鬼成像相机原理

Fig. 1. Principle for ghost imaging camera

下载图片查看所有图片

将式（1）简写为矢量形式，得到：

Y = Φ X

其中： $Y = (y_{1}, y_{2}, \dots, y_{M})^{T}$ ，代表由像元数为 $M$ 的探测器采集到的场景的散斑场的图像矢量； $X = (x_{1}, x_{2}, \dots, x_{N})^{T}$ ，代表与原始场景共轭的第一成像面的图像矢量；Φ为鬼成像相机的测量矩阵，第 $i$ 列元素 $Φ_{i} = (Φ_{1}, Φ_{2}, \dots, Φ_{M})^{T}$ ，代表第一成像面上第 $i$ 个面元产生的散斑场的图像矢量。

由式（1）可知，在实际成像之前，需要对空间随机相位调制器进行预先标定以获得测量矩阵Φ。标定方式如图1（b）所示：将第一成像面划分为 $N$ 个面元，第一成像面处的针孔光阑按划分面元的顺序依次移动，由像元数为 $M (M ≪ N)$ 的探测器采集不同面元辐射的散斑场图像，最后将散斑场图像拉成一个列向量，并作为测量阵Φ的某一列。

为了对比不同算法在图1（a）所示的被动式鬼成像系统中的成像效果，需要对系统中的光场调制模块进行数学建模，从而便于鬼成像重建算法的数值仿真。

2.2　光场调制模型

鬼成像相机的光场调制原理如图2所示。第一成像面上一点出射的热光场经过 $z_{1}$ 距离的自由传播到达空间随机相位调制器（本文使用散射片作为空间随机相位调制器），经过光场调制后，再与空间随机相位调制器后表面相距 $z_{2}$ 的平面形成散斑场。

图 2. 光场调制原理

Fig. 2. Schematic diagram of light field modulation

下载图片查看所有图片

由于第一成像面上各点光源近似为空间不相干，因此第一成像面上任意两点 $S_{1} (ξ_{1}, η_{1})$ 以及 $S_{2} (ξ_{2}, η_{2})$ 之间的互强度可表示为：

J (S_{1}, S_{2}) = I (S_{1}) δ (ξ_{1} - ξ_{2}, η_{1} - η_{2})

其中 $δ (\cdot)$ 表示狄拉克函数。第一成像面上点 $S_{1}$ 处的热光场在自由空间中传播 $z_{1}$ 后到达空间随机相位调制器的前表面，成为空间部分相干光。此时，空间随机相位调制器前表面上任意两点 $P_{1} (x_{1}, y_{1})$ ， $P_{2} (x_{2}, y_{2})$ 的互强度为：

J (P_{1}, P_{2}) = \frac{1}{{\bar{λ}}^{2} R_{1} R_{2}} \iint \iint J (S_{1}, S_{2}) e x p [j \bar{k} (R_{2} - R_{1})] \times K (θ_{1}) K (θ_{2}) d ξ_{1} d η_{1} d ξ_{2} d η_{2}

其中： $\bar{λ}$ 为热光场的中心波长， $\bar{k}$ 为中心波长对应的波矢， $K (θ_{1}), K (θ_{2})$ 为倾斜因子。在傍轴近似条件下有：

K (θ_{1}) K (θ_{2}) \approx 1

R_{1} \approx z_{1} + \frac{(x_{1} - ξ_{1})^{2} + (y_{1} - η_{1})^{2}}{2 z_{1}}

R_{2} \approx z_{1} + \frac{(x_{2} - ξ_{2})^{2} + (y_{2} - η_{2})^{2}}{2 z_{1}}

将式（5）~式（7）代入式（4），并利用狄拉克函数的性质进行积分，化简得：

J (P_{1}, P_{2}) = \frac{e x p {j \bar{k} [(x_{1}^{2} + y_{1}^{2}) - (x_{2}^{2} + y_{2}^{2})] / (2 z_{1})}}{(\bar{λ} z_{1})^{2}} \cdot \iint I (S_{1}) e x p {- j \bar{k} [(x_{1} - x_{2}) ξ + (y_{1} - y_{2}) η] / z_{1}} d ξ d η,

其中： $J (P_{1}, P_{2})$ 表示 $ξ - η$ 平面上任意两点传播到 $x - y$ 平面上得到的互强度分布。假设在第一成像面上的 $(ξ_{0}, η_{0})$ 处放置一个透光直径为 $d$ 的针孔光阑来构造小面元，令：

I (S_{1}) = c i r c (\frac{\sqrt[]{(ξ - ξ_{0})^{2} + (η - η_{0})^{2}}}{d / 2})

其中： $c i r c (\cdot)$ 为圆域函数。将式（9）代入式（8），并计算积分得：

\begin{matrix} J (P_{1}, P_{2}) = \frac{d^{2}}{(2 \bar{λ} z_{1})^{2}} e x p {j \bar{k} [(x_{1}^{2} + y_{1}^{2}) - (x_{2}^{2} + y_{2}^{2})] / (2 z_{1})} \times \\ e x p {- j \bar{k} [(x_{1} - x_{2}) ξ_{0} + (y_{1} - y_{2}) η_{0}] / z_{1}} \times \frac{J_{1} (\frac{\sqrt[]{(x_{1} - x_{2})^{2} + (y_{1} - y_{2})^{2}}}{\bar{λ} z / (π d)})}{\frac{\sqrt[]{(x_{1} - x_{2})^{2} + (y_{1} - y_{2})^{2}}}{2 π \cdot \bar{λ} z / (π d)}}, \end{matrix}

其中 $J_{1} (\cdot)$ 为一阶第一类贝塞尔函数。那么散斑场上任意两点 $Q_{1} (s_{1}, t_{1})$ 和 $Q_{2} (s_{2}, t_{2})$ 的互强度为：

\begin{matrix} J (Q_{1}, Q_{2}) = \frac{1}{{\bar{λ}}^{2} L_{1} L_{2}} \iint \iint J (P_{1}, P_{2}) \times \\ e x p {j \bar{k} (L_{1} - L_{2})} \times e x p {j [ϕ (P_{1}) - ϕ^{*} (P_{2})]} \times \\ K (θ_{3}) K (θ_{4}) d x_{1} d y_{1} d x_{2} d y_{2}, (11) \end{matrix}

其中 $ϕ (P_{1}), ϕ (P_{2})$ 分别表示空间随机相位调制器对点 $P_{1}, P_{2}$ 处的光场调制作用， $K (θ_{3}), K (θ_{4})$ 为倾斜因子。在傍轴近似条件下有：

K (θ_{3}) K (θ_{4}) \approx 1

L_{1} = z_{2} + \frac{(s_{1} - x_{1})^{2} + (t_{1} - y_{1})^{2}}{2 z_{2}}

L_{2} = z_{2} + \frac{(s_{2} - x_{2})^{2} + (t_{2} - y_{2})^{2}}{2 z_{2}}

将式（12）~（14）代入式（11）得：

\begin{array}{l} J (Q_{1}, Q_{2}) = \frac{e x p {j \bar{k} [(s_{1}^{2} + t_{1}^{2}) - (s_{2}^{2} + t_{2}^{2})] / 2 z_{2}}}{(\bar{λ} z_{2})^{2}} \times \iint \iint J (P_{1}, P_{2}) e x p {j [ϕ (P_{1}) - ϕ^{*} (P_{2})]} \times \\ e x p {j \bar{k} [(x_{1}^{2} + y_{1}^{2}) - (x_{2}^{2} + y_{2}^{2}) - 2 (s_{1} x_{1} - s_{2} x_{2}) - 2 (t_{1} y_{1} - t_{2} y_{2})] / (2 z_{2})} d x_{1} d y_{1} d x_{2} d y_{2} . \end{array}

令 $Q_{1} = Q_{2}$ ，并将式（10）代入式（15）得到点 $Q_{1}$ 处的光场强度为：

\begin{matrix} I (Q_{1}) = \frac{d^{2}}{(2 {\bar{λ}}^{2} z_{1} z_{2})^{2}} \times \iint \iint e x p {j \bar{k} [(x_{1}^{2} + y_{1}^{2}) - (x_{2}^{2} + y_{2}^{2})] \times [1 / (2 z_{1}) + 1 / (2 z_{2})]} \times \\ e x p {- j \bar{k} [(x_{1} - x_{2}) (ξ_{0} / z_{1} + s / z_{2}) + (y_{1} - y_{2}) (η_{0} / z_{1} + t / z_{2})]} \times \\ \frac{J_{1} (\frac{\sqrt[]{(x_{1} - x_{2})^{2} + (y_{1} - y_{2})^{2}}}{\bar{λ} z / (π d)})}{\frac{\sqrt[]{(x_{1} - x_{2})^{2} + (y_{1} - y_{2})^{2}}}{2 π \cdot \bar{λ} z / (π d)}} d x_{1} d y_{1} d x_{2} d y_{2} . \end{matrix}

设第一成像面上有两个面元分别位于 $(ξ_{0}, η_{0})$ ， $(ξ_{0} + Δ ξ, η_{0} + Δ η)$ 处，对应的散斑场的光场强度分别为 $I (s, t)$ ， $I (s + Δ s, t + Δ t)$ 。由式（16）可知，满足：

Δ ξ / z_{1} + Δ s / z_{2} \approx 0

Δ η / z_{1} + Δ t / z_{2} \approx 0

此时有 $I (s, t) = I (s + Δ s, t + Δ t)$ 。这说明在空间随机相位调制器后表面附近形成的散斑干涉场存在空间相干性。因此，已知第一成像面上某点 $(ξ, η)$ 对应的散斑场的光强分布，可以通过平移得到其近邻点的散斑场的光强分布。

通过测量第一成像面上某一面元对应的散斑场，而后通过平移构造一个矩阵，将这个矩阵作为测量矩阵Φ，即可进行成像数值模拟。根据散斑场的空间相干性，可以大大减少空间随机相位板标定的工作量。在标定过程中等间隔地采集第一成像面处面元产生的散斑场，其近邻面元对应的散斑场由已知的散斑场通过平移得到。

3 算法原理

欠采样条件下的鬼成像实质上是一个不适定问题。为了得到合理的图像估计，需要充分合理利用稀疏/可压缩信号的几何结构。为此本文提出了联合迭代的图像重构算法，代替传统的压缩感知方法来解决最小化问题。在实际应用中，由于采样次数有限，欠采样噪声一直存在，不能忽略。为了提高成像的稳定性，减小欠采样噪声的影响，提高重构图像的信噪比，算法首先利用共轭梯度法来解决不适定问题，得到关联成像的预重构图像，然后基于图像的非局域自相似性，利用低秩聚类降噪算法来抑制图像的预重构噪声，提高成像质量。经过图像预重构和低秩聚类降噪的多次联合迭代，充分利用空间调制光场中所包含的有效物体信息，即可在较低采样率下得到成像效果较佳的重构图像。该算法流程如图3所示。

图 3. 联合迭代的重构算法流程

Fig. 3. Flow chart of reconstruction algorithm with joint iteration

下载图片查看所有图片

3.1　图像预重构

图像预重构是求解一个最小二乘问题。预重构问题描述为：给定测量值 $Y \in R^{M}$ 以及测量矩阵 $Φ_{N \times M}$ ，确定系数 $\hat{X}$ 使得：

\hat{X} = a r g m i n | | Y - Φ \hat{X} | |_{2}

等价于求解如下函数的极值：

F (\hat{X}) = \frac{1}{2} {\hat{X}}^{T} Φ^{T} Φ \hat{X} - {\hat{X}}^{T} Φ^{T} Y

式（20）属于无约束二次规划问题，可以利用共轭梯度法进行快速迭代，使得函数 $F (\hat{X})$ 取得极小值的 $\hat{X}$ ，以达到快速恢复图像细节的目的^［29］。共轭梯度法在最速下降法的基础上结合共轭梯度信息来获得目标函数的最优解，不仅克服了最速下降法收敛速度慢的不足，相比于牛顿法它没有计算Hessen矩阵和求逆过程，也避免了存储空间的占用，是解决大型非线性优化问题非常有效的算法。

3.2　低秩聚类降噪

图像预重构解决的是无约束的最小二乘问题，不可避免地会产生预重构噪声。基于图像的非局域自相似性，通过对含噪的预重构图像进行低秩建模，可以将目标信息从背景噪声中提取出来，并抑制重构噪声。

对于噪声图像中某个参考块 $l$ ，在搜索窗口中选择一组与该含噪图像块 $l$ 相似度较高的图像块，相似度由欧几里得范数定义。不妨设有 $N_{l}$ 个与 $l$ 相似度较高的图像块（包含 $l$ ），将这些相似的块向量堆叠起来，构成一个矩阵：

X_{l} = [x_{1}, x_{2}, \dots, x_{j}, \dots x_{N_{l}}]

由于自然图像的非局域自相似性，由非局域相似块构成的矩阵为低秩矩阵，因此，式（21）中由含噪图像的相似块堆叠而成的矩阵通过低秩建模为：

X_{l} = S_{l} + Γ_{l}

其中： $S_{l}$ 为无噪图像的块矩阵， $Γ_{l}$ 为包含噪声的块矩阵。 $S_{l}$ 可以通过低秩矩阵近似方法获得。低秩矩阵近似的一个解决方法是核范数最小化（Nuclear Norm Minimization， NNM）^［30］，即：

{\hat{S}}_{l} = a r g \underset{S_{l}}{m i n} | | X_{l} - S_{l} | |_{F}^{2} + λ | | S_{l} | |_{*}

其中： $λ$ 为常数， $| | S_{l} | |_{*}$ 为矩阵 $S_{l}$ 的核范数，定义为 $S_{l}$ 奇异值的和，即 $| | S_{l} | |_{*} = \sum_{i} | | σ_{i} | |_{1}$ ， $σ_{i}$ 表示 $S_{l}$ 的第 $i$ 个奇异值。NNM的目标是通过 $S_{l}$ 近似 $X_{l}$ ，同时最小化 $S_{l}$ 的核范数。对于式（23）所描述的问题，Cai^［31］等证明通过对观测矩阵奇异值进行软阈值运算可以很容易地解决基于NNM的低秩矩阵逼近问题，其中保真项由 $F$ 范数描述。因此，式（25）可由奇异值软阈值方法解得：

{\hat{S}}_{l} = U S_{λ} (Σ) V^{T}

其中： $X_{l} = U Σ V^{T}$ 为 $X_{l}$ 的奇异值分解， $S_{λ} (Σ)$ 为带参数 $λ$ 的对角矩阵 $Σ$ 的软阈值函数。对于 $Σ$ 中每个对角线上的元素 $Σ_{i i}$ ，有：

S_{λ} {(Σ)}_{i i} = m a x (Σ_{i i} - λ, 0)

尽管NNM已广泛用于低秩矩阵逼近，但它仍然存在一些问题。为了追求凸性，标准核范数平等地对待每个奇异值，因此，式（24）中的软阈值算子以相同的量 $λ$ 收缩每个奇异值。然而，这忽略了通常对矩阵奇异值的先验知识。例如，矩阵中的列（或行）向量通常位于低维子空间中；较大的奇异值通常与主要的投影方向相关，因此期望较大奇异值能够收缩得少一点，以保留主要的数据分量。显然，NNM及其相应的软阈值算子未能利用这种先验知识。尽管式（23）中的模型是凸的，但它不够灵活，无法处理许多实际问题。为了提高核范数的灵活性，本文采用加权Schatten-p范数进行低秩矩阵近似。矩阵 $S_{l}$ 的加权Schatten-p范数定义为^［32］：

| | S_{l} | |_{w, s} = {(\sum_{i} w_{i} σ_{i} (S_{l})^{p})}^{1 / p}, 0 < p \leq 1

其中： $σ_{i} (S_{l})$ 为 $S_{l}$ 的奇异值， $p$ 为作用于矩阵奇异值 $σ_{i} (S_{l})$ 的幂。对于自然图像， $S_{l}$ 中的较大奇异值比较小奇异值更重要，这是因为它们代表 $S_{l}$ 主要分量的能量。在降噪过程中，奇异值越大，它们应当收缩得越小。因此，分配给 $σ_{i} (S_{l})$ 的权重应与 $σ_{i} (S_{l})$ 成反比，即：

w_{i} = c \sqrt[]{n} / (σ_{i} (S_{l})^{1 / p} + ε)

其中： $c$ 为常数， $n$ 为含噪图像 $X_{l}$ 的相似块数量， $ε$ 为一个避免除零的较小常数。在初始时刻，假设噪声能量在各个特征上分布是均匀的，因此初始化 $σ_{i} (S_{l})$ 为：

σ_{i} (S_{l}) = \sqrt[]{m a x (σ_{i}^{2} (X_{l}) - n σ_{n}^{2}, 0)}

其中： $σ_{i} (X_{l})$ 为 $X_{l}$ 的第 $i$ 个奇异值， $σ_{n}^{2}$ 为 $X_{l}$ 的噪声方差。

将式（26）所示的加权Schatten-p范数应用于式（23），使用噪声方差 $σ_{n}^{2}$ 来归一化 $F$ 范数数据保真项，得到如下目标函数：

{\hat{S}}_{l} = a r g \underset{S_{l}}{m i n} \frac{1}{σ_{n}^{2}} | | X_{l} - S_{l} | |_{F}^{2} + | | S_{l} | |_{w, s}^{p}

其中 $| | S_{l} | |_{w, s}^{p}$ 为 $S_{l}$ 的加权Schatten-p范数的p次方。Xie^［32］证明目标函数（26）的全局最优解可以通过广义软阈值^［33］（Generalized Soft-Thresholding，GST）算法进行求解。因此， $S_{l}$ 的解为：

S_{l} = U \times G S T (Σ, w_{i}, p) \times V

其中 $p = 0.7$ 。对所有含噪图像块执行低秩聚类降噪操作，并将去噪后的图像块聚合得到一幅降噪后的重建图像，随后将降噪后的图像重新返回到3.1进行共轭梯度下降操作。

本文在被动压缩鬼成像问题中引入低秩聚类降噪的方法，得到了一种成像效果较好的被动式鬼成像算法，称该算法为低秩聚类被动压缩鬼成像算法（Passive Compressed Ghost Imaging with Low Rank Clustering， PCGI-LRC）。该算法运行中，最小二乘问题和低秩矩阵近似问题交替求解。当迭代达到一定次数时，停止迭代，得到最终的重构图像。基于低秩聚类联合迭代关联成像方法流程如下：

表 1.

Table 1.

算法1：低秩聚类被动压缩鬼成像算法

输入：测量值 $Y$ ，测量矩阵 $Φ$ ，最大迭代次数 $K$ ，迭代终止阈值 $ε$

1：初始化： $X_{0} = 0$

2：令 $i = 1 : K$

3：采用共轭梯度法求解最小二乘问题，获得预重

构图像 ${\hat{X}}^{(k + 1)}$ ：

${\hat{X}}^{(k + 1)} = a r g \underset{{\hat{X}}^{(k)}}{m i n} 1 / 2 {\hat{X}}^{(k)}^{T} Φ^{T} Φ {\hat{X}}^{(k)} - {\hat{X}}^{(k)}^{T} Φ^{T} y$

4：对图像 ${\hat{X}}^{(k + 1)}$ 进行分块得到局部块 ${\hat{x}}_{l}^{(k + 1)}$

5：寻找局部块 ${\hat{x}}_{l}^{(k + 1)}$ 的相似块矩阵组 ${\hat{X}}_{l}^{(k + 1)}$

6：对每个有噪声的数据矩阵 ${\hat{X}}_{l}^{(k + 1)}$ 进行奇异值分解： $(U_{l}^{(k + 1)}, Σ_{l}^{(k + 1)}, V_{l}^{(k + 1)}) = S V D ({\hat{X}}_{l}^{(k + 1)})$

7：利用广义软阈值算法进行低秩矩阵近似

8：图像更新，对所有块进行加权平均，得到噪声

抑制后的图像 $q^{(k + 1)}$

9：若 $(q^{(k + 1)} - {\hat{X}}^{(k)}) / {\hat{X}}^{(k)} < ε$ ，则停止循环，设置 $X = q^{(k + 1)}$ ；否则设置 ${\hat{X}}^{(k)} = q^{(k + 1)}$ ，联合迭代继续

10：若迭代次数 $i > K$ ，则停止循环

输出：重构图像 $X$

查看所有表

4 仿真结果与分析

为了测试PCGI-LRC的成像效果，首先进行数值仿真，将它与GI-TVC、离散余弦域下的GI-SBC两类算法进行对比。为了客观地评价该方法的性能，用峰值信噪比（Peak Signal-to-Noise Ratio， PSNR）和结构相似性系数（Structural Similarity Index， SSIM）来定量地衡量图像重建质量。这里选取OMP算法^［8］和FISTA算法^［11］作为GI-SBC的重建算法，选取TVAL3^［18］和GISCNL^［19］作为GI-TVC的重建算法。

本文选择灰度图像“house”（图像尺寸为128×128）作为仿真中需要重建的目标场景。基于式（18）和式（19）的结论，将第一成像面划分为128×128块面元，随后等间隔地选取第一成像面上各面元对应的散斑场，其余面元对应的散斑场通过平移获得，最后将所有面元的散斑场拉成一维向量并聚合成二维矩阵作为被动式鬼成像系统的测量矩阵。

在非局域相似块低秩近似的过程中，图像块的尺寸会影响成像性能。若图像块尺寸过大，成像结果会更平滑；如果图像块尺寸过小，则无法达到抑制重构噪声的目的。在鬼成像过程中，需要在保持图像细节和抑制重构噪声之间权衡。在4，5，6之间选取图像块的尺寸，并根据仿真或实验得到的最佳成像结果调整图像块的尺寸。从图像的左上角选择第一个参考块，并按照先列后行的顺序从一个参考块切换到下一个。数值仿真是在Intel（R） Core （TM） CPU i7-11800H处理器（2.3 GHz），16 G内存电脑上的MATLAB R2021a版本上运行的。

不同采样率下，GI-SBC，GI-TVC和PCGI-LRC（图4（b）中蓝色曲线）对“house”图像的重构结果如图4（a）所示，对应的数值仿真结果计算的PSNR曲线如图4（b）所示（彩图见期刊电子版）。对比结果可知：利用不同算法在不同采样率条件下进行鬼成像，本文算法的PSNR值均优于GI-SBC和GI-TVC，且目测结果与原始图像更接近，重建伪影更少。特别当采样率为0.5时，GI-SBC，GI-TVC和PCGI-LRC的PSNR值分别为22.744 0，22.913 0和28.307 4 dB。从这些重构图像和评价曲线上可以看出，PCGI-LRC在主观和客观上都优于其他方法。

图 4. 灰度图像仿真结果

Fig. 4. Simulation results of gray images

下载图片查看所有图片

如上所述，PCGI-LRC是一个最小二乘问题和非局域相似块低秩近似问题的联合迭代求解过程。为了验证该方法每一步骤的有效性，在采样率为0.25的情况下，“house”图像在不同迭代次数 $i$ 下的均方误差（Mean Squared Error，MSE）的曲线如图5所示。由图5可知，随着迭代次数的增加，MSE曲线呈单调递减趋势，最终趋于平稳，表明本文方法具有较好的稳定性。还可以从中观察到，约10次迭代即可得到最佳的重构图像。当MSE基本稳定或迭代次数 $i$ 达到最大迭代次数 $K$ 时，停止迭代，得到最终的重构图像。

图 5. 采样率为0.25的条件下本文算法在不同迭代次数下的MSE变化曲线

Fig. 5. MSE variation curve of PCGI-LRC under different iterations at sampling rate of 0.25

下载图片查看所有图片

5 测量实验与结果

为验证算法的实际性能，本文搭建了被动式压缩鬼成像实验平台，在采样率分别为0.062 5，0.125，0.25以及0.5的情况下进行鬼成像实验，实验光路如图6所示。采用LED灯（KM-2FL400400， 400 mm×400 mm）作为照明光源，场景“GI”辐射出来的光经过准直镜头（Collimating Lens， LBMT-2，f=400 mm）准直后经由分束镜（Beam Splitter， BS1455-A， 50：50）进行分光：其中一束光由监视相机（Monitor， TRI050S-MC，3 μm×3 μm）进行采集，监视相机将采集到的原始场景图像作为实验中的对比图像，另一束光由被动压缩鬼成像系统（Passive Compressed Ghost Imaging）采集（图中虚线框所标注）。在被动式鬼成像系统中，前置成像镜头（LBMT-1，f=200 mm）将光场聚焦在第一成像面处，空间随机相位板（GCL-201123， 1500）将光场信息随机编码至干涉场中，最后由显微物镜（KP-10A035X，β=22）将干涉场成像至CCD（TRI050S-MC， 3 μm×3 μm）中。其中，散射片与第一成像面相距12.70 mm，散斑场与散射片后表面相距1.5 mm。

图 6. 被动式压缩鬼成像实验光路

Fig. 6. Experimental optical path of passive compressed ghost imaging

下载图片查看所有图片

在实际成像前需要对空间随机相位调制器进行预先标定，标定过程如图1（b）所示。在第一成像面处放置针孔光阑作为标定过程中的点光源，针孔光阑固定在电位移台上，可以通过计算机调整针孔光阑的位置。空间随机相位调制器对透过针孔光阑的光进行相位调制，并在CCD的探测面上形成具有强度起伏的光强分布，即散斑场。标定过程中，将第一成像面划分为 $100 \times 100$ 个面元，第一成像面处的针孔光阑按照面元的划分顺序依次进行移动。为提高相位板的标定效率，针孔光阑在第一成像面上每移动3个面元的距离，计算机则控制CCD采集一次散斑场，并作为鬼成像系统中参考臂的赝热光场。其余面元对应的散斑场则由已知的散斑场通过双线性插值法计算得到。标定完成后，移除针孔光阑与电位移台便可直接对场景进行成像。实际成像中，在图6所示的准直镜头的焦平面处放置靶标，用于模拟无穷远处的场景。鬼成像相机收集靶标的散斑场信息并送入上位机（Intel Core i7-11800H RAM 16GB），利用鬼成像算法进行图像重构。原始场景的图像尺寸为 $100 \times 100$ ，如图7所示。首先，按不同采样率对CCD采集到的散斑场图像进行空间随机采样。例如，对于0.062 5的空间采样率，需要随机选取散斑场图像中 $0.062 5 \times 100 \times 100 = 625$ 个像素点作为式（4）中的列向量 $Y$ 。随后基于鬼成像原理，利用重构算法对散斑场进行反演，最终得到场景的高分辨率图像。

图 7. 目标图像及其散斑场

Fig. 7. Target image and its scattering field

下载图片查看所有图片

最后，将PCGI-LRC与GI-SBC，GI-TVC两类成像算法进行对比，成像目标为刻有字母“SITP”的靶标。不同鬼成像算法对场景的成像效果以及客观评价指标分别如图8和表1所示。

图 8. 不同算法的成像效果对比

Fig. 8. Comparison of imaging effects of different algorithms

下载图片查看所有图片

10.37188/OPE.20243201.0012.T002

表1

不同鬼成像算法的重建效果对比

Tab.1Comparison of reconstruction effects of different ghost imaging algorithms

Algorithm	Rate=0.062 5			Rate=0.125			Rate=0.25			Rate=0.5
Algorithm	SSIM	PSNR/dB	MSE	SSIM	PSNR/dB	MSE	SSIM	PSNR/dB	MSE	SSIM	PSNR/dB	MSE
GI-SBC（OMP）	0.118 4	16.396 0	0.149 1	0.138 9	16.892 5	0.133 0	0.141 3	17.012 1	0.128 7	0.147 2	17.863 2	0.106 3
GI-SBC（FISTA）	0.079 0	14.572 2	0.226 9	0.092 6	15.535 3	0.181 8	0.105 6	16.448 2	0.147 3	0.123 3	17.752 4	0.109 1
GI-TV（TVAL3）	0.091 7	17.247 6	0.122 6	0.092 4	17.388 6	0.118 6	0.104 6	17.673 9	0.110 9	0.106 0	17.284 9	0.121 5
GI-TV（GISCNL）	0.104 1	17.662 1	0.111 4	0.113 4	17.731 0	0.109 6	0.115 0	17.950 1	0.104 3	0.115 8	18.026 8	0.102 4
PCGI-LRC	0.147 2	18.751 5	0.086 8	0.152 6	18.784 9	0.086 0	0.156 3	19.178 4	0.078 6	0.153 5	19.164 6	0.078 8

结合图8的成像效果与表1所示的数据可以看出，GI-SBC类算法能够较好地恢复目标的细节信息，但重构噪声较大；GI-TVC类算法虽然相比于GI-SBC类算法拥有更低的重构噪声，但对于图像细节的重构能力不足；本文提出的算法能够较好地恢复出目标的细节，同时有效地抑制重构噪声。显然，本文提出的鬼成像算法在PSNR，SSIM以及MSE 3项指标上均优于另外两类成像算法。

6 结论

为了在低采样率条件下减小图像重构噪声的影响，提高图像的重构精度，本文提出了一种被动式压缩鬼成像方案，该方法对最小二乘问题和非邻域相似块低秩近似问题进行联合迭代求解。在联合迭代过程中，首先采用共轭梯度法求解最小二乘问题来获得预重构图像，随后利用图像非局域自相似性的特性以滑动窗口的方式处理预重构图像，其中滑动窗口尺寸固定。将这些相似块的阵列堆叠成低秩矩阵，通过广义软阈值算法收缩衰减预重构噪声。仿真和实验结果均表明，被动式鬼成像体制下，本文算法在PSNR、SSIM、视觉观测和数值计算等方面均优于GI-SBC以及GI-TVC算法，未来有望应用于目标检测、遥感等领域。

参考文献

[1] STREKALOV D V, SERGIENKO A V, KLYSHKO D N, et al. Observation of two-photon “ghost” interference and diffraction[J]. Physical Review Letters, 1995, 74(18): 3600-3603.

[2] PITTMAN T B, SHIH Y H, STREKALOV D V, et al. Optical imaging by means of two-photon quantum entanglement[J]. Physical Review A, 1995, 52(5): R3429-R3432.

[3] BENNINK R S, BENTLEY S J, BOYD R W. “two-photon” coincidence imaging with a classical source[J]. Physical Review Letters, 2002, 89(11): 113601.

[4] SCARCELLI G, BERARDI V, SHIH Y. Quantum interference approach to two-photon correlation phenomena of chaotic light[J]. Journal of Modern Optics, 2006, 53(16/17): 2279-2292.

[5] SHAPIRO J H. Computational ghost imaging[J]. Physical Review A, 2008, 78(6).

[6] LIU Z T, TAN S Y, WU J R, et al. Spectral camera based on ghost imaging via sparsity constraints[J]. Scientific Reports, 2016, 6: 25718.

[7] 孙鸣捷，闫崧明，王思源. 鬼成像和单像素成像技术中的重建算法［J］. 激光与光电子学进展， 2022， 59（2）： 0200001. doi: 10.3788/LOP202259.0200001SUNM J， YANS M， WANGS Y. Reconstruction algorithm in ghost imaging and single pixel imaging technology［J］. Laser & Optoelectronics Progress， 2022， 59（2）： 0200001.（in Chinese）. doi: 10.3788/LOP202259.0200001

[8] TROPP J A, GILBERT A C. Signal recovery from random measurements via orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655-4666.

[9] NEEDELLD， VERSHYNINR. Greedy signal recovery and uncertainty principles［C］.Proc SPIE 6814， Computational Imaging VI， 2008， 6814： 139-150. doi: 10.1117/12.776996

[10] CHEN S S, DONOHO D L, SAUNDERS M A. Atomic decomposition by basis pursuit[J]. SIAM Review, 2001, 43(1): 129-159.

[11] BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM Journal on Imaging Sciences, 2009, 2(1): 183-202.

[12] CANDES E J, TAO T. Near-optimal signal recovery from random projections： universal encoding strategies？[J]. IEEE Transactions on Information Theory, 2006, 52(12): 5406-5425.

[13] CANDES E J, 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.

[14] BIOUCAS-DIASJ M， FIGUEIREDOM A T. Two-step algorithms for linear inverse problems with non-quadratic regularization［C］.2007 IEEE International Conference on Image Processing. September 16-October 19， 2007. San Antonio， TX， USA. IEEE， 2007， 1： I-105-I-108. doi: 10.1109/icip.2007.4378902

[15] BECKER S, BOBIN J, CANDÈS E J. NESTA： a fast and accurate first-order method for sparse recovery[J]. SIAM Journal on Imaging Sciences, 2011, 4(1): 1-39.

[16] YUANX. Generalized alternating projection based total variation minimization for compressive sensing［C］.2016 IEEE International Conference on Image Processing （ICIP）. September25-28， 2016. Phoenix， AZ， USA. IEEE， 2016： 2539-2543.

[17] WAHLBERG B, BOYD S, ANNERGREN M, et al. An ADMM algorithm for a class of total variation regularized estimation problems[J]. IFAC Proceedings Volumes, 2012, 45(16): 83-88.

[18] LI C B, YIN W T, JIANG H, et al. An efficient augmented Lagrangian method with applications to total variation minimization[J]. Computational Optimization and Applications, 2013, 56(3): 507-530.

[19] WANG P W, WANG C L, YU C P, et al. Color ghost imaging via sparsity constraint and non-local self-similarity[J]. Chinese Optics Letters, 2021, 19(2).

[20] KATZ O, BROMBERG Y, SILBERBERG Y. Compressive ghost imaging[J]. Applied Physics Letters, 2009, 95(13): 131110.

[21] BROMBERG Y, KATZ O, SILBERBERG Y. Ghost imaging with a single detector[J]. Physical Review A, 2009, 79(5).

[22] LI W W, TONG Z S, XIAO K, et al. Single-frame wide-field nanoscopy based on ghost imaging via sparsity constraints[J]. Optica, 2019, 6(12): 1515-1523.

[23] LIU S, LIU Z, HU C, et al. Spectral ghost imaging camera with super-Rayleigh modulator[J]. Optics Communications, 2020, 472: 126017.

[24] CHU C Y, LIU S Y, LIU Z T, et al. Spectral polarization camera based on ghost imaging via sparsity constraints[J]. Applied Optics, 2021, 60(16): 4632-4638.

[25] WANG T, CHEN M Y, WU H, et al. Underwater compressive computational ghost imaging with wavelet enhancement[J]. Applied Optics, 2021, 60(23): 6950-6957.

[26] HUANG H, ZHOU C, TIAN T, et al. High-quality compressive ghost imaging[J]. Optics Communications, 2018, 412: 60-65.

[27] WU G H, LI T H, LI J H, et al. Ghost imaging under low-rank constraint[J]. Optics Letters, 2019, 44(17): 4311-4314.

[28] 吴建荣，沈夏，喻虹，等. 基于相位调制的单次曝光压缩感知成像［J］. 光学学报， 2014， 34（10）： 1011005. doi: 10.3788/aos201434.1011005WUJ R， SHENX， YUH， et al. Snapshot compressive imaging by phase modulation［J］. Acta Optica Sinica， 2014， 34（10）： 1011005.（in Chinese）. doi: 10.3788/aos201434.1011005

[29] 穆哲. 电容层析成像图像重建算法研究［D］. 兰州：西北师范大学， 2021.MUZH. Research on Image Reconstruction Algorithm of Electrical Capacitance Tomography［D］. Lanzhou： Northwest Normal University， 2021. （in Chinese）

[30] CANDÈS E, RECHT B. Exact matrix completion via convex optimization[J]. Communications of the ACM, 2012, 55(6): 111-119.

[31] CAI J F, CANDÈS E J, SHEN Z W. A singular value thresholding algorithm for matrix completion[J]. SIAM Journal on Optimization, 2010, 20(4): 1956-1982.

[32] XIE Y, GU S H, LIU Y, et al. Weighted schatten p-norm minimization for image denoising and background subtraction[J]. IEEE Transactions on Image Processing, 2016, 25(10): 4842-4857.

[33] ZUOW M， MENGD Y， ZHANGL， et al. A generalized iterated shrinkage algorithm for non-convex sparse coding［C］.2013 IEEE International Conference on Computer Vision. December1-8， 2013. Sydney， Australia. IEEE， 2013： 217-224.

雷腾, 张义民, 马一哲, 丁学专, 吴滢跃, 王世勇. 低秩聚类被动压缩鬼成像[J]. 光学精密工程, 2024, 32(1): 12. Teng LEI, Yiming ZHANG, Yizhe MA, Xuezhuan DING, Yingyue WU, Shiyong WANG. Passive compressive ghost imaging with low rank clustering[J]. Optics and Precision Engineering, 2024, 32(1): 12.