高回波环境下基于阵列探测的三维计算鬼成像

High-Echo Environment-Based 3D Single-Photon Ghost Imaging via Array Detection

摘要: 传统的单光子三维计算鬼成像(3DSCGI)由于在极低光子通量环境下工作，信噪比(SNR)较低，需要进行大量的重复测量。提高光子通量可以改善这一问题，但会引入堆积效应，导致光子计数直方图与原始脉冲波形相比发生畸变。本研究开发了一种新的单光子三维计算鬼成像技术(3DHCGI)，在高回波光子通量下直接计算64*64单光子相机的光子计数和哈达码矩阵之间的二阶相关函数(SOCF)，显著提高了信噪比，同时避免了堆积效应的影响。此外，利用哈达码矩阵行中平衡的+1和−1分布，我们采用互补探测进一步降低噪声，同时在每个散斑模式仅探测126次的情况下，完成了39.45 m处场景256*256分辨率的三维重建，横向分辨率达到了7.21 mm，显著提升了成像质量。这项工作对于中远距离高分辨率单光子三维成像具有重要意义，提供了一种高效且高质量的成像方法。

Abstract: Typical single-photon 3D computational ghost imaging (3DSCGI) suffers from low signal-to-noise ratios (SNR) due to operating in ultra-low photon flux environment, necessitating numerous repetitive measurements. Enhancing photon flux improves this but introduces pile-up effects, distorting photon counts histogram compared with original pulse waveform. Our study develops a new single-photon 3D computational ghost imaging technology (3DHCGI) computing the second-order correlations function (SOCF) between photon counts of 64*64 single photon camera and Hadamard matrix directly under high echoing photon flux, significantly boosting SNR and avoiding pile-up effect based on the pairwise orthogonality of Hadamard matrix. Additionally, leveraging the balanced +1 and −1 distribution in Hadamard rows, we utilize complementary detection to finish 256*256 3D reconstruction at distance 39.45 m with transverse resolution 7.21 mm and detections of each pattern 126, and further reduce noise and enhance image quality. This work is important for high-resolution single-photon 3D imaging, offering an efficient and high-quality imaging method.

1. 引言

近年来，单光子三维成像技术因其卓越的场景空间结构恢复能力而备受关注[1]-[3]。其中，三维单光子计算鬼成像(3DSCGI)技术通过数字微镜阵列(Digital Micromirror Device, DMD)生成独立散斑模式，并利用单光子雪崩光电二极管(Single Photon Avalanche Diode, SPAD)作为桶探测器，实现了无需扫描的三维成像[4]-[11]。该技术通过计算光子计数与散斑模式之间的二阶关联函数(SOCF)，结合时间关联单光子计数(TCSPC)测量桶信号光子的飞行时间(TOF)，从而重建三维场景($D=c\ast \text{TOF}/2$ ，D表示系统与场景中物体的距离，c为光的传播速度)。此外，压缩感知的引入突破了奈奎斯特采样定理的限制，显著提升了成像质量和效率[12]-[15]。

然而，3DSCGI [16]-[18]仍面临诸多挑战。传统单光子探测器在探测回波信号时，光子计数遵循泊松计数模型[19]-[21]，但由于死区时间限制，仅能记录一个探测周期内首次探测到光子的时间标签。这导致在特定时间间隔内光子事件的发生实际上服从条件概率模型，从而在增加回波光子通量时产生堆积效应[22] [23]。在高回波光子通量下，光子计数时间统计直方图与原始脉冲波形相比会出现明显的前倾现象[24]。这种现象使得每个时间间隔内的光子计数不再与到达探测器的光子数保持固定的线性关系，在鬼成像架构中，这将在成像背景区域引入额外的噪声。为避免堆积效应的影响，3DSCGI工作在极低回波光子通量环境中，从而导致单次探测信噪比很低，对于一种散斑模式往往需要重复测量很多次，一定程度上增加了探测过程的繁琐性，限制其实际应用。

为解决堆积效应的问题，提升探测过程信噪比，研究人员提出了多种方法，以实现从扭曲的时间统计直方图中恢复到达探测器的原始信号[25]，例如基于马尔可夫链的直方图畸变校正算法[26]、光子计数校正方法[27]、基于光子计数分布直方图的恢复方法[28]，以及基于概率图像模型减少范围误差的方法[29]。Liu等人通过信号校正方法在红外单光子压缩成像系统中消除了负像[30] [31]。尽管这些方法在一定程度上恢复了原始信号并提高了检测过程中的信噪比(SNR)，但复杂的恢复过程还是在某种程度上限制了其在实际应用中的广泛推广。

为了克服这些限制，本文提出了一种新的三维单光子计算鬼成像技术(3DHCGI)。该技术通过增加回波光子通量显著提高了单次检测SNR，从而减少了重复检测的次数，并实现了高分辨率成像。通过采用各列成对正交的哈达码矩阵作为测量矩阵，在高光子通量下直接计算光子计数与哈达码散斑模式之间的SOCF，有效避免了成像背景区域中堆积效应的影响[32]。此外，我们利用哈达码矩阵每行中平衡的+1和−1分布，采用互补探测进一步降低噪声并提高图像质量[33]。在计算鬼成像中，当散斑场的横向尺寸大小保持稳定时，最终成像的SNR与散斑场中被物体占据的像素数量成反比[34]。因此，我们采用64*64的SPAD相机作为桶探测器，相机阵列上的每个像素作为一个独立子桶探测器探测其所对应立体角范围内的桶信号，从而进一步提高成像质量。

2. 实验原理

Figure 1. Illustration of the 3DHCGI imaging system. The light source 1064 nm Laser emits pulses with pulse width 2 ns and energy of a single pulse being 10 μJ. Predefined patterns of $H^{+}$ and $H^{-}$ are generated on DMD, then alternatively projected into the scenario at the distance D = 39.45 m by the lens $f_1$ = 150 mm. The objects in the scenario are imaged on the camera through lens $f_2$ with adjustable focal length and the camera pixels detect the bucket photons corresponding to their respective solid angle range. And the detecting process is controlled by PC through signal generator (SG)

图1. 3DHCGI成像系统示意图。1064 nm激光器发射脉宽为2 ns、单脉冲能量为10 μJ的脉冲。在DMD上生成与$H^{+}$ 和$H^{-}$ 相应的预定义图案，然后通过焦距为150 mm的透镜$f_1$ 交替投射到距离为39.45 m的场景中。场景中的物体通过具有可调焦距的透镜$f_2$ 成像到相机上，相机像素检测其各自立体角范围对应的桶信号。检测过程由计算机通过信号发生器(SG)控制

图1展示了3DHCGI成像系统的示意图。光源为1064 nm的脉冲激光器，脉宽2 ns，单脉冲能量10 μJ。通过DMD 生成相应于$H^{+}$ 和$H^{-}$ 的预定义散斑模式，并通过焦距为150 mm的透镜$f_1$ 交替投影到距离为39.45 m的场景中。场景中的物体通过焦距可调的透镜$f_2$ 成像到64*64 SPAD相机上，相机像素独立探测其各自立体角范围内的桶信号。整个探测过程由计算机通过信号发生器(SG)控制。

哈达码矩阵由+1和−1组成，各列之间具有成对正交性[31]。在3DHCGI中，我们将哈达码矩阵中所有+1的元素组成矩阵$H^{+}$ ，并将所有−1的元素转换为+1以形成矩阵$H^{-}$ 。然后将$H^{+}$ 和$H^{-}$ 的每一行重组成一个独立的散斑模式，例如，当哈达码矩阵元素组成为4096*4096时，散斑模式的大小为64*64，对应成像的分辨率也为64*64。3DHCGI系统通过控制DMD表面的微镜阵列生成一系列与$H^{+}$ 和$H^{-}$ 对应的散斑模式，并交替投影到场景中。由于不同距离下的点扩散函数在成像距离较远时几乎没有差异，因此可以将不同距离下的散斑场近似视为相同，即散斑场平行穿过场景。随后，SPAD相机根据桶信号光子的TOF区分不同纵向距离，并利用哈达码散斑与不同时间标签范围内的光子计数之间的SOCF重建不同距离处的物体图像，从而完成三维重构。探测过程中每个相机像素作为子桶探测器独立探测其各自立体角范围内的桶信号。

对于4096个独立相机像素中某个i像素，在探测$H^{+}$ 各散斑模式桶信号过程期间，光子事件的发生满足条件概率模型，其概率为[35]：

$P_i^{+}={\text{e}}^{-s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}\left(1-{\text{e}}^{-s_1{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}\right)$ (1)

其中，$R\left(X_2\right)$ 表示场景中的物体反射率分布，$\gamma$ 为3DHCGI系统的衰减系数，$s_0$ 和$s_1$ 分别表示属于$0-t_1$ 和$t_1-t_2$ 时间间隔内脉冲激光器发出的光子数。${\eta }_i$ 表示i像素的量子效率，${\Omega }_i$ 表示在$t_1-t_2$ 时间间隔内被i相机像素探测到的物体区域。从公式(1)可以看出，当3DHCGI系统工作在高回波环境中时，$0-t_1$ 时间间隔内未检测到光子事件的概率会变小，而SPAD在一个探测周期中只会产生1个光子计数，因此光子计数的统计直方图会出现前倾，也就是堆积效应的原因，从而导致该时间间隔内光子计数和回波光子数之间不具有确定的线性关系，关联运算后会在背景区域增加额外的噪声。堆积效应导致的信号扭曲问题可以通过光子计数满足的条件概率模型进行解决，如文献[24]中所述，从而消除堆积效应对背景区域的影响。由于3DSCGI工作在极低回波环境中，$0-t_1$ 时间间隔内的光子通量可以近似为0，因此可以将SPAD的每一个时间分辨率内发生的计数过程看成独立的泊松计数过程，从而忽略堆积效应的影响，根据光子事件发生概率的函数形式($1-{\text{e}}^{-x}\approx x$ ，当$x$ 趋近于0时)，可以近似认为在低回波环境下，光子计数正比于每个探测周期内SPAD所接收到的光子数，代替桶信号进行二阶关联运算。但是在3DHCGI中，需要考虑$P_i^{+}$ 表达式中的每一个影响因子。因此，$H^{+}$ 和$P_i^{+}$ 在3DHCGI中的SOCF可以表示为：

$\begin{array}{c}{O}_i^{+}\left(X_1\right)=O_{0i}^{+}\left(X_1\right)-O_{1i}^{+}\left(X_1\right)\\ =\langle H^{+}\left(X_1\right)\rangle {\text{e}}^{-\left(s_0+s_1\right){\eta }_i\gamma \langle H_0\rangle R_{sum}}\left(1-{\text{e}}^{-\left(s_0+s_1\right){\eta }_i\gamma \langle H_0\rangle R\left(X_1\right){\delta }_{1i}}\right)\\ -\langle H^{+}\left(X_1\right)\rangle {\text{e}}^{-s_0{\eta }_i\gamma \langle H_0\rangle R_{sum}}\left(1-{\text{e}}^{-s_0{\eta }_i\gamma \langle H_0\rangle R\left(X_1\right){\delta }_{1i}}\right)\end{array}$ (2)

其中$R_{sum}={\displaystyle {\int }_{{\Omega }_i}R\left(X_1\right)\text{d}{X}_1}$ ，$\langle \cdots \rangle$ 表示系综平均，$H_0$ 表示在散斑场被照亮时单位面积的光子数，${\delta }_{1i}$ 是对应于$H^{+}\left(X_1\right)$ 的散斑场像素与被i相机像素探测到的物体区域的重叠面积，详细证明见附录。此公式也适用

于$H^{-}$ 。从公式 (2)可以看出，成像能量集中在物体区域，因此堆积效应不会在背景区域引入额外噪声。我们进一步利用互补探测消除光子计数中的噪声。由于哈达码矩阵的稀疏率为0.5，因此一对互补图案之间的噪声大致相同。最终的SOCF可以表示为：

$O\left(X_1\right)=2\ast {\left(H^{+}-\langle H^{+}\rangle \right)}^{\text{T}}\left(P_i^{+}-P_i^{-}-\langle P_i^{+}-P_i^{-}\rangle +n^{+}-n^{-}-\langle n^{+}-n^{-}\rangle \right)$ (3)

其中，$n^{+}$ 和$n^{-}$ 分别表示与$H^{+}$ 和$H^{-}$ 对应的互补模式探测过程中引入的噪声，详细证明见附录。这些噪声成分包括环境噪声、暗噪声、统计噪声、串扰噪声以及由于相机在高回波光子通量环境下产生非线性响应引起的噪声。由于在探测过程中，3DHCGI系统和所处环境处于稳定状态，因此每次探测过程中环境噪声的强度可以认为服从均值不变的均匀分布，而暗噪声在每个探测过程中均满足特征量一定的泊松分布，该特征量由探测器本身的属性和工作状态所决定。统计噪声取决于所有探测过程中所发生的随机计数过程，在3DHCGI中，某段时间间隔内光子事件的发生服从条件概率分布，估计该时间间隔内光子计数的统计噪声时可以使用二项分布来进行估计，该二项分布的发生概率取决于每个散斑模式对应的相机入射光子数。相机的串扰概率为固有串扰概率和入射光子数的乘积，由于每个散斑模式的回波光子数均不相同，因此在每个散斑模式的探测过程中，串扰噪声均有其各自不同的统计分布。而高回波环境下相机的非线性响应引起的噪声通常和入射光子通量呈正相关的非线性关系，因此对于不同的散斑模式，其所对应的噪声也是不同的。3DHCGI最大的优势在于通过提升回波光子数，极大程度提升了随机计数过程的SNR，减少了统计噪声所带来的影响，同时互补探测在抑制环境噪声和暗噪声方面非常有效—它们在两种互补图案的探测过程中表现出相同的统计特性，使得环境/暗噪声的方差相减后理论上为零—但在抑制统计噪声、串扰噪声以及相机非线性响应引起的噪声方面效果较差。

本文成像质量通过SNR来评估，表示为：

$\text{SNR}=10\ast {\log }_{10}\frac{{\displaystyle {\int }_{{\Omega }_{object}}O\left(X_1\right)\text{d}{X}_1}/{\displaystyle {\int }_{{\Omega }_{object}}\text{d}{X}_1}}{{\displaystyle {\int }_{{\Omega }_{object}}O\left(X_1\right)\text{d}{X}_1}/{\displaystyle {\int }_{{\Omega }_{object}}\text{d}{X}_1}}$ (4)

其中，${\Omega }_{object}$ 和${\Omega }_{background}$ 分别表示所成图像中物体区域和背景区域的像素数量。

3. 实验配置及结果

3.1. 实验配置

实验配置部分将详细描述3DHCGI的实验装置。如图1所示，使用一个波长为1064 nm的固体激光器作为光源，该激光器单脉冲能量10 μJ，脉宽2 ns，最大脉冲频率10 KHz，本文实验过程中脉冲激光器工作频率为2 KHz。脉冲光束发射后，经反射镜(Mirror, M)反射，通过可变倍数扩束器(Variable Beam Expander, VBE)，形成均匀光斑投射到DMD表面微镜阵列。DMD表面由1024*768个微镜组成，每个微镜尺寸为13.65 μm*13.65 μm，由可编程门控阵列控制，允许每个微镜独立倾斜 ± 12˚，分别对应开(ON)和关(OFF)状态。因此，DMD能够实现对平行光场的0-1调制，从而交替生成相应于$H^{+}$ 和$H^{-}$ 的散斑模式，在经过一个孔径光阑(Adjustable Diaphragm, AD)滤除高频成分后，通过焦距为150 mm的透镜$f_1$ 投射到距离为$D=39.45\text{m}$ 的场景中。实验中，DMD工作频率为2 KHz，每一个散斑模式对应于1个脉冲探测周期。

3DHCGI采用64*64的SPAD相机(型号为GD5551，由重庆四十四所生产)作为桶探测器，其阵列尺寸为3.2 mm*3.2 mm，单个像素尺寸为50 μm*50 μm，感光材料为InGaAs。该相机具有1ns的时间分辨率和最大25 KHz的采样频率。当一定数量的光子到达像素表面时，相机的光子计数遵循泊松计数过程。相机所能设置的最大门宽为4096 μs，如果光子事件发生在时间间隔$t_1-t_2$ 内，读出电路将仅记录一次光子计数，并记录与该光子事件相应的时间标签，前提是在$t_0$ 之前没有光子事件发生并被记录，实际光子计数过程本质上遵循条件概率模型。读出电路集成并封装在相机主体内，具有225 ns的固有链路延迟。在实验过程中，相机在−20℃，2 KHz采集频率，84 V外加电压下工作，平均量子效率为0.2 (在入射光子波长为1064 nm的情况下)，同时设置门宽值133 ns。场景通过相机前的可调焦距透镜$f_2$ 成像于相机阵列，使得每个相机像素独立探测其所对应立体角范围内的桶信号。

相机一旦接收到激光器的外部触发信号，在经过固有链路延迟后即刻开启门宽。同时，读出电路以每纳秒增加1的速率开始计数。在门宽开启时间内到达像素表面并被探测到的光子将被记录其到达时刻相对于脉冲发射时刻所对应的时间标签；如果没有发生光子事件，则记录实验中所设置的门宽值。实验中在1个脉冲探测周期后，相机传输的数据是1*4096的行向量(由64列光子事件时间数据组成)，通过CameraLink协议传输到计算机中进行存储。通过将多次脉冲探测后的数据重组为64*64*L的矩阵(L为探测过程中1个散斑模式所经历的探测脉冲数)，我们可以准确确定在每个散斑模式的探测过程中，每个相机像素所对应的立体角范围内的光子计数，探测过程和数据结构如图2所示。整个3DHCGI系统由计算机通过信号发生器(SG)控制。

(a)

(b)

Figure 2. (a) The detecting process, the pulse is emitted from the 3DHCGI system and arrives at the camera in the period of gate width. After an inherent link delay of 225 ns, the readout circuit starts counting at a rate of 1 ns per increment before the occurring of photon event: $\Delta t_0$ is the inherent link delay, $\Delta t_1$ is the gate width, $\Delta t_{flight}$ is the photon flight time, and T is the detecting cycle duration respectively; (b) data structure, detected row vector 1*4096L is aggregated into a matrix 64*64*L. Different camera pixels are characterized by different color corresponding to different solid angle range of their own and gate width 133 ns will be recorded when no photon event occurs

图2. (a) 探测过程，脉冲由3DHCGI系统发射，并在门宽开启期间到达相机。在固有链路延迟225 ns后，读出电路以每纳秒加1的速率开始计数，直到光子事件发生：$\Delta t_0$ 是固有链路延迟，$\Delta t_1$ 是门宽，$\Delta t_{flight}$ 是光子飞行时间，T是1个脉冲探测周期的持续时间；(b) 数据结构，探测到的1*4096L行向量被重组为一个64*64*L三维矩阵。不同的相机像素以不同颜色表示，对应于它们各自所探测的不同立体角范围，当该像素没有光子事件发生时，将记录设置的门宽值133 ns

3.2. 实验结果

后续部分将比较3DHCGI与3DSCGI的成像结果，突出3DHCGI中所采用的互补探测方法的优势，并展示3DHCGI的三维重建能力。所有实验均在2 KHz的采样频率下进行。在128*128成像实验中，DMD表面每个独立像素由一个4*4的微镜阵列组成，而在256*256成像实验中，每个像素由一个2*2的微镜阵列组成。并且，在整个实验过程中，系统和场景均保持稳定，以确保结果的一致性和可靠性。

3.2.1. 3DHCGI与3DSCGI的对比

图3(a)和图3(b)展示了在超低回波环境(3DSCGI)和高回波环境(3DHCGI)下，使用互补探测方法，并且使用41~45 ns之间的光子计数作为桶信号，每种散斑模式探测的脉冲数为120条件下获得的128*128成像结果。图3(c)和图3(d)显示了这两种环境下的光子计数时间统计直方图。与工作在超低光子通量环境下的3DSCGI相比，3DHCGI的成像质量显著提高。并且由于哈达码矩阵各列向量之间的成对正交性，高回波环境下的堆积效应(图3(d)相对于图3(c)明显的前倾) 并未在成像背景区域引入额外噪声。该实验充分说明了在利用哈达码矩阵的成对正交性避免了堆积效应对背景区域的影响后，3DHCGI相较于传统的3DSCGI有着更强的成像能力，该成像能力的提升，来自于提升回波概率后，随机计数过程产生的光子计数SNR明显变大，统计噪声的影响明显变小。

(a) L = 200 SNR = 16.42 dB (b) L = 200 SNR = 24.48 dB

(c) 3DSCGI光子计数时间统计直方图 (d) 3DHCGI光子计数时间统计直方图

Figure 3. (a), (b) are the 128 *128 complementary detecting imaging results with 41~45 ns camera photon counts and L = 200 pulses per pattern under ultra-low and high echo environments, respectively; (c), (d) are photon counts time histogram of camera under ultra-low and high echo environments, respectively

图3. (a)和(b)分别是超低回波环境和高回波环境下，使用互补探测方法和41~45 ns内的相机光子计数作为桶信号，每种散斑模式经过200个脉冲探测条件下的128*128成像结果；(c)和(d)分别是超低回波环境和高回波环境下相机的光子计数时间统计直方图

3.2.2. 互补检测性能

图4(a)~(j)展示了仅使用$H^{+}$ 散斑模式时的128*128成像结果，将41~45 ns之间的光子计数作为桶信号，每个散斑模式经历的探测脉冲数L从12到120。图4(k)~(t)展示了使用互补探测方法时的结果，也是将41~45 ns之间的光子计数作为桶信号，L从6到60变化。图4(u)比较了两种探测方法的SNR。互补探测方法的SNR始终高于仅使用$H^{+}$ 散斑的方法，且随着L的增加，这种差异变得更大。这是因为，在两种方法选定的探测脉冲数如图4所示的情况下，成像物体的像素值理论上相等，而在仅使用$H^{+}$ 散斑模式的方法中，环境噪声和暗噪声的方差随L成比例增加，在互补探测方法中，这些噪声的方差理论上为零。但是如实验原理部分所论述，由于高回波环境中，互补图案之间的统计噪声、串扰噪声以及相机非线性响应引起的噪声均与每个散斑模式各自对应的回波光子数有关，不具有相同的统计分布，因此即使通过互补探测做为去噪方法，理论上在将互补的散斑模式对应的桶信号作差后，这些噪声所对应方差不为0，成像结果中依然存在一定数量的噪声。

Figure 4. (a)~(j) are the 128*128 imaging results using only $H^{+}$ , with photon counts recorded between 41~45 ns and the number of pulses per pattern L ranging from 12 to 120; (k)~(t) are the 128*128 imaging results using the complementary detection method, with L ranging from 6 to 60; (u) is the performance comparison between the two detection methods

图4. (a)~(j) 是在仅使用$H^{+}$ 的散斑模式情况下的128*128成像结果，将41~45 ns之间的光子计数作为桶信号，每种散斑模式经历的探测脉冲数L从12到120；(k)~(t) 是使用互补探测方法的128*128成像结果，L从6到60；(u) 是两种探测方法的性能对比

3.2.3. 256*256 3D重构结果

图5(a)和图5(b)展示了在高回波环境下，使用互补探测方法获得的256*256成像结果。分别使用35~39 ns和41~45 ns之间的光子计数作为桶信号，每种散斑模式经历的探测脉冲数为126。图5(c)是根据桶信号光子的TOF进行的3D重建。实验环境中，3DHCGI系统的横向分辨率为7.21 mm，可以基于投影透镜$f_1$ 焦距(150 mm)，TOF (38 ns) 以及单个微镜的尺寸(13.65 μm*13.65 μm)计算得出，但是由于堆积效应在光子计数时间统计直方图中并没有消除(图3(c)和图3(d))，因此虽然3DHCGI拥有更高的成像SNR，但是对于场景中物体纵向距离的测量误差更大。这些结果证明了3DHCGI在每种散斑模式仅经历126个脉冲探测周期的情况下，仍具有出色的3D重建能力，相较于传统的3DSCGI，大幅度减少了单个散斑模式所需的重复探测次数，如文献[36]中每种散斑模式需探测100,000次，文献[35]中每种散斑模式需探测50,000次。3DHCGI能够大幅缩减单个散斑模式探测次数的原因在于其工作在高回波环境下，在堆积效应不对成像背景区域产生影响的情况下，极大程度提升了随机计数过程中光子计数的SNR，使得统计噪声的影响很小，并且使用每一个相机像素作为子桶探测器，减少了传统鬼成像架构中桶探测器所对应的散斑模式中物体的像素数，同时，互补探测的应用进一步减少了环境噪声和暗噪声的影响。

(a) L = 126 SNR = 19.06 dB (b) L = 126 SNR = 16.18 dB

(c) 三维重构

Figure 5. (a), (b) are the 256*256 complementary detecting imaging results with photon counts recorded between 35~39 ns and 41~45 ns, respectively, and L = 126 pulses per pattern under high echo environment; (c) is the 3D reconstruction according to the TOF of the bucket photons

图5. (a)和(b)是在高回波环境下，使用互补探测方法获得的256*256成像结果。分别使用35~39 ns和41~45 ns之间的光子计数作为桶信号，每种散斑模式经历的探测脉冲数为126；(c)是根据桶信号光子的TOF进行的3D重建

4. 结论与讨论

本文提出了一种新型的3DHCGI方法，能够在39.45 m距离处实现256*256的三维场景重建，横向分辨率达到7.21 mm，并且每个散斑模式仅需重复探测126次。这一结果突显了3DHCGI在中长距离场景中卓越的三维重建能力。相较于3DSCGI，本文提出的3DHCGI具有以下优势。

1) 工作在高回波光子通量环境中，提升了单次探测的信噪比，减少了光子计数中统计噪声的影响，从而极大程度减少了单个散斑模式重复测量的次数，提升了单光子成像系统的实用性；

2) 利用哈达码矩阵列向量之间的成对正交性，有效避免了堆积效应对成像背景区域的影响；

3) 采用互补探测的方式，在相同数量探测周期的情况下，削弱了环境噪声和暗噪声的影响，进一步提升了成像SNR。

4) 不再受限于超低回波的实验条件，3DHCGI的应用场景更加灵活。

尽管3DHCGI实现了较高的成像质量，但堆积效应仍然是一个挑战，时间统计直方图依然会出现前倾，如图3(c)和图3(d)所示，从而降低时间分辨率。除了本文讨论的噪声外，相机在运行过程中阵列的不稳定性也可能在成像过程中引入额外噪声。然而，本文所展示的实验结果充分证明了3DHCGI优异的成像能力。我们的工作对于中长距离单光子成像应用具有重要意义，特别是在高回波光子通量可以实现的场景中。未来的工作将致力于进一步优化系统，以解决剩余的挑战，如堆积效应和相机不稳定性，从而提高3DHCGI在实际成像场景中的整体性能和准确性。

致 谢

谢谢导师上海理工大学王春芳教授的指导，表示衷心的感谢。

附 录

实验原理部分证明

以下是公式(2)的证明：

本文中$t_j$ 到$t_{j+1}$ 时间间隔内，脉冲激光器发出的光子数为：

$s_j=I_0{\displaystyle {\int }_{t_j}^{t_{j+1}}S\left(t\right)\text{d}t}$ (5)

$S\left(t\right)$ 和$I_0$ 分别为脉冲的时间波形以及脉冲激光器单脉冲所包含总的光子数。

在$H^{+}$ 散斑模式的探测过程中，相机阵列上i像素的回波概率如公式(1)所示，因此$H^{+}$ 散斑模式和$P_i^{+}$ 的SOCF可以写为：

$\begin{array}{c}{O}_i^{+}\left(X_1\right)=\langle H^{+}\left(X_1\right)P_i^{+}\rangle -\langle H^{+}\left(X_1\right)\rangle \langle P_i^{+}\rangle \\ =\langle H^{+}\left(X_1\right){\text{e}}^{-s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}\rangle -\langle H^{+}\left(X_1\right)\rangle \langle {\text{e}}^{-s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}\rangle \\ -\langle H^{+}\left(X_1\right){\text{e}}^{-(s_0+s_1){\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}\rangle +\langle H^{+}\left(X_1\right)\rangle \langle {\text{e}}^{-s\left(s_0+s_1\right){\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}\rangle \\ =O_{0i}^{+}\left(X_1\right)-O_{1i}^{+}\left(X_1\right)\end{array}$ (6)

关于参数的描述如正文中所述。

由于3DHCGI工作在高回波环境中，因此需要考虑${\text{e}}^{-s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}$ 和${\text{e}}^{-\left(s_0+s_1\right){\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}$ 的高阶Tailor展开项。

$\begin{array}{l}{\text{e}}^{-s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}}\\ =1-s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}+\frac{{\left(s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}\right)}^2}{2!}\\ -\frac{{\left(s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}\right)}^3}{3!}+\cdots \end{array}$ (7)

因此$H^{+}$ 散斑模式和一阶Tailor展开项的SOCF为：

$\begin{array}{l}\langle H^{+}\left(X_1\right)s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}\rangle -\langle H^{+}\left(X_1\right)\rangle \langle s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}\rangle \\ =s_0{\eta }_i\gamma \langle H^{+}\left(X_1\right)\rangle \langle H_0\rangle R\left(X_1\right){\delta }_{1i}\end{array}$ (8)

$H_0$ 表示散斑场中某像素亮起时，单位面积入射光子数，其他参数如正文中所述。

$H^{+}$ 散斑模式和二阶Tailor展开项的SOCF为：

$\begin{array}{l}\langle H^{+}\left(X_1\right){\left(s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}\right)}^2\rangle -\langle H^{+}\left(X_1\right)\rangle \langle {\left(s_0{\eta }_i\gamma {\displaystyle {\int }_{{\Omega }_i}{H}^{+}\left(X_1\right)R\left(X_1\right)\text{d}{X}_1}\right)}^2\rangle \\ ={\left(s_0{\eta }_i\gamma \right)}^2{\displaystyle {\int }_{{\Omega }_i}{\displaystyle {\int }_{{\Omega }_i}\left(\langle H^{+}\left(X_1\right)H^{+}\left(X_2\right)H^{+}\left(X_3\right)\rangle -\langle H^{+}\left(X_1\right)\rangle \langle H^{+}\left(X_2\right)H^{+}\left(X_3\right)\rangle \right)R\left(X_2\right)R\left(X_3\right)\text{d}{X}_2\text{d}{X}_3}}\\ ={\left(s_0{\eta }_i\gamma \right)}^2\langle H^{+}\left(X_1\right)\rangle {\langle H_0\rangle }^2{\displaystyle {\int }_{{\Omega }_i}{\displaystyle {\int }_{{\Omega }_i}(\langle \left(H\left(1\right)+1\right)\left(H\left(2\right)+1\right)\left(H\left(3\right)+1\right)\rangle }}\\ -\langle H\left(1\right)+1\rangle \langle \left(H\left(2\right)+1\right)\left(H\left(3\right)+1\right)\rangle )R\left(X_2\right)R\left(X_3\right)\text{d}{X}_2\text{d}{X}_3\\ ={\left(s_0{\eta }_i\gamma \right)}^2\langle H^{+}\left(X_1\right)\rangle {\langle H_0\rangle }^2{\displaystyle {\int }_{{\Omega }_i}{\displaystyle {\int }_{{\Omega }_i}\langle H\left(1\right)H\left(2\right)H\left(3\right)\rangle +\langle H\left(1\right)H\left(2\right)\rangle +\langle H\left(1\right)H\left(3\right)\rangle \text{d}{X}_2\text{d}{X}_3}}\\ ={\left(s_0{\eta }_i\gamma \right)}^2\langle H^{+}\left(X_1\right)\rangle {\langle H_0\rangle }^2\left[{\left(R\left(X_1\right){\delta }_{1i}+R_{sum}\right)}^2-R_{sum}^2\right]\end{array}$ (9)

公式(9)中，$R_{sum}={\displaystyle {\int }_{{\Omega }_i}R\left(X_1\right)\text{d}{X}_1}$ 。

同理，$H^{+}$ 散斑模式和三阶Tailor展开项的SOCF为：

${\left(s_0{\eta }_i\gamma \right)}^3\langle H^{+}\left(X_1\right)\rangle {\langle H_0\rangle }^3\left[{\left(R\left(X_1\right){\delta }_{1i}+R_{sum}\right)}^3-R_{sum}^3\right]$ (10)

因此，我们可以得到$O_{0i}^{+}\left(X_1\right)$ 的表达式为：

$O_{0i}^{+}\left(X_1\right)=\langle H^{+}\left(X_1\right)\rangle {\text{e}}^{-s_0{\eta }_i\gamma \langle H_0\rangle R_{sum}}\left({\text{e}}^{-s_0{\eta }_i\gamma \langle H_0\rangle R\left(X_1\right){\delta }_{1i}}-1\right)$ (11)

$O_{1i}^{+}\left(X_1\right)$ 的计算过程和$O_{0i}^{+}\left(X_1\right)$ 一致，因此，$O_i^{+}\left(X_1\right)$ 的表达式如公式 (2)所示。

以下是公式 (3)的证明：

在3DHCGI系统中，$H^{+}$ 和$H^{-}$ 的散斑模式被交替投影到场景中，最终的二阶关联运算过程可以表述为：

$\begin{array}{l}{O}_i^{+}\left(X_1\right)={\left(H^{+}-\langle H^{+}\rangle \right)}^{\text{T}}\left(P_i^{+}-\langle P_i^{+}\rangle \right)\\ O_i^{-}\left(X_1\right)={\left[\left(1-H^{+}\right)-\langle 1-H^{+}\rangle \right]}^{\text{T}}\left(P_i^{-}-\langle P_i^{-}\rangle \right)\end{array}$ (12)

我们可以使组合桶信号取$P_i^{+}-P_i^{-}$ 的形式完成互补探测，进一步减少光子计数中的噪声影响，从而获得$O_i\left(X_1\right)$ 的SOCF表达式为：

$O_i\left(X_1\right)=2\ast {\left(H^{+}-\langle H^{+}\rangle \right)}^{\text{T}}\left(p_i^{+}-p_i^{-}-\langle p_i^{+}-p_i^{-}\rangle +n^{+}-n^{-}-\langle n^{+}-n^{-}\rangle \right)$ (13)

关于噪声成分及其影响的分析如正文中所述。