超宽带被动单光子成像 论文解读（Passive Ultra-Wideband Single-Photon Imaging）

Hua Zhouqi, Department of Computer Science & Technology, Tongji University

一、全文概要

文章提出了一种被动超宽带单光子成像方法，通过使用单光子相机和光子检测时间戳，实现了对动态场景在秒到皮秒级时间尺度上的成像，无需大量光线和光源的定时信号，并成功实现了多种前所未见的成像能力。

二、先验知识了解

单光子相机

每个像素只检测一个光子的相机，可以形成非常清晰的图像。主要技术难点在于单光子探测和图像重构。

工作原理为使用一种名为**单光子雪崩二极管（SPAD）**的技术，这是一种能够探测单光子的高灵敏度半导体器件。

Ref：https://zhuanlan.zhihu.com/p/344969590

主动成像

相机和光源使用固定频率的脉冲信号进行同步，

主动成像涉及使用外部光源照亮正在拍摄的场景。 光源发射光子，然后光子从场景中的物体上反射。 主动成像系统检测并测量反射光子以创建图像。 这种方法可以控制照明条件（可以对成像环境进行改变），并可以提供增强的图像质量和深度感知。

被动成像（Passive imaging）

被动成像不依赖外部光源（无需依靠信号反射），它可以捕捉场景中自然存在的环境光或辐射。

被动成像系统工作原理是检测并记录场景中物体发射或反射的光子或辐射。 这种方法通常在光线不足的条件下使用，或者当使用外部光源不可行或不可取时。

超宽带（Ultra-Wideband）

跨越从恒定通量 (0 Hz) 到极端飞行时间尺度 (≥ 10 GHz) 的整个范围的频率带宽，是一种在极端的时间尺度上同时对动态场景进行成像的方法，即几秒到皮秒。

文中的超宽带成像指的是动态视频的成像频率可以有非常大的跨度，从DC（直流频率）到30GHz左右的极高频。

单光子雪崩二极管（SPAD）

固态光电探测器，本文中作为光子接收装置（单光子相机传感器），将模拟光电转为可识别的数字信号。

击穿原理（了解）：强电场下，在晶体中运行的电子和空穴将不断的与晶体原子发生碰撞，通过这样的碰撞可使束缚在共价键中的价电子碰撞出来，产生自由电子-空穴对。新产生的载流子在电场作用下撞出其他价电子，又产生新的自由电子和空穴对。如此连锁反应，使得耗尽层中的载流子的数量雪崩式地增加，流过PN结的电流就急剧增大，这种碰撞电离导致击穿称为雪崩击穿。

Ref：https://zhuanlan.zhihu.com/p/405389221

通量估计技术

用于估计给定场景或图像中光通量或强度的方法。在用于低光水平成像的单光子相机的背景下，通量估计技术在重建光子探测时间戳流中的时变通量方面起着至关重要的作用。

本文中利用随机微积分的见解开发了一种通量探测理论，以实现像素的时变通量的重建，即从单调增加的光子探测时间戳流中估算通量。

非视距成像（NLOS）

能够捕捉由于障碍物或遮挡物而无法直接被摄像机看见的物体或场景的视频片段。它涉及使用间接光或散射光来重建隐藏的场景。

Ref：https://www.zhihu.com/question/387241900/answer/2269213825

异步成像

无需在相机和光源之间进行同步即可捕捉图像或视频的方法。

在传统的成像系统中，相机和光源是同步的，以确保捕获的图像与特定时刻的照明相对应。在异步成像中，相机独立于光源运行，从而在捕捉动态场景时具有更大的灵活性和多功能性。

非齐次泊松过程（Inhomogeneous poisson process）

某事件在任何时刻发生的速率（概率）为 $\lambda$ ，我们假设速率 $\lambda$ 和时间 $t$ 有关，记为 $\lambda (t)$ ，这样的过程我们称为非齐次的泊松过程（Poisson process）。

齐次泊松过程的特点是平稳增量。但是当事件的频度不恒定，产生非齐次泊松过程，不具有平稳增量。

Ref：https://zhuanlan.zhihu.com/p/28785318

奈奎斯特定理（Nyquist theorem）

又称奈奎斯特-香农采样定理，信号处理和通信理论中的一个基本概念，是主动光子直方成像（Active histogram-based imaging）作为通量采集情况下最大可重建频率的主导原理。

主要内容：为了准确地重建来自其样本的连续信号，采样速率必须至少是信号中存在的最高频率分量的两倍。

对应采样频率称为奈奎斯特极限（Nyquist limit）。

傅立叶谱（Fourier spectrum）

傅里叶频谱是通过对通量函数执行傅里叶变换获得的，该函数将信号从时域转换为频域。生成的频谱提供有关通量函数中存在的不同频率分量的振幅和相位的信息。

Ref：

[1] https://zhuanlan.zhihu.com/p/19759362

[2] https://zhuanlan.zhihu.com/p/34989414

函数内积（inner product）

在傅立叶级数知识中，函数的内积用于描述两个函数之间的关系（相关性）。现规定两函数 $f(x)$ 与 $g(x)$ 与区间 $[a,b]$，且两函数在该区间上可积且平方可积。积分 ${\int }_a^{b}f(x)g(x)dx$ 记作函数的内积，记作 $<f(x),g(x)>$。

鞅（martingale）

鞅是一种连续时间的随机过程，满足：已知过去某一时刻 $s$ 以及之前所有时刻的观测值，若某一时刻 $t$ 的观测值的条件期望等于过去某一时刻 $s$ 的观测值，则称这一随机过程是鞅。

概率里，其表示的则是一类既无向上趋势，又无向下趋势的随机过程(Stochastic processes)，文中用作可加零均值噪声。

傅立叶基函数（Fourier basis functions）和 傅立叶变换

傅立叶分析中使用的一组函数，用于将周期信号表示为正弦函数之和。这些函数以下方程定义：
$$p_f(t)=e^{-j2\pi ft}$$
其中 $f$ 是频率，$t$ 是时间。

在复平面中，可以理解为以原点为圆心，$1/f$ 为周期，半径为1的顺时针圆形运动。

傅里叶基函数构成了表示频域信号的完整正交基础。它们使我们能够将信号分解为其组成频率，并确定每个频率分量的振幅和相位。

傅立叶变换中通常使用信号函数 $g(t)$ 与基函数相乘后对时间积分，将整个函数转为关于频率 $f$ 的函数，从而根据图像可以得到实部的（多个）峰值，进而可以进行不同频率信号的分离：
$$\hat{g}(f)={\int }_{-\infty }^{+\infty }g(t)e^{-2\pi ift}dt$$

Ref：https://www.bilibili.com/video/BV1pW411J7s8/?spm_id_from=333.999.0.0&vd_source=b4f15dea13e1086ebea8e3a493d35c8c

傅立叶混淆（Aliasing）

简单来说就是一种由于采样不足而导致信号的高频分量被错误地表示为低频分量的现象。

在成像环境中，当尝试从使用高于一定极限的频率获得的测量结果中重建信号时，可能会出现傅里叶混叠混乱。该极限由奈奎斯特定理和所使用的采样率确定。如果正在探测的频率超过此限制，则生成的测量值将被混叠，并且高频信息将在重建过程中失真或丢失。

定时分辨率（Timing resolution）

计时系统或设备可以精确测量或解析的最小时间间隔。它表示可以测量或记录时间的精度。

文中按照相机单光子雪崩二极管（SPAD）的时序分辨率，整个采集间隔实际上被划分为时间箱（bin）。

CFAR（Constant false alarm rate）

恒定虚警概率下的检测器，作用为在含有噪声的情况下确定信号存在还是不存在。

本文中使用CFAR的工作原理是：设置所需的误报概率 $\alpha$，并将每个频率的估计傅里叶探测能量 $|p_f(T){|}^2$ 与从非中心卡方分布（non-central ${\chi }^2$ distribution ）得出的阈值进行比较。

• 如果频率的能量测量超过阈值（$CD{F}_{{\chi }^2}^{-1}(1-\alpha )N(t)/2t^2$），则该频率被视为具有统计学意义并被检测为有效频率。
• 如果能量测量值低于阈值，则该频率被视为有噪音并从分析中删除。

关于阈值的设定，CFAR 探测器会考虑探测到的光子总数（$N(t)$）和通量函数中特定频率的主导性，以调整探测到其他频率的概率。

三、介绍

在高速光源、快速相机等硬件加持下，高速成像已经可以实现，但是有以下的问题还没有被解决：

• 当前高速成像的基本法则是：场景变化越快，就需要越多的光线来精确成像才不会有过多的噪声和运动模糊，所以在低通量环境无法实现。
• 需要操作相机和信号源之间的同步频率，使用相同的重复频率成像。
• 捕捉超快事件的同时无法同时捕捉较慢事件，因为同步周期是无法准确捕获事件的边界。时间在同步周期结束，在较长时间跨度内发生的任何事情都会被模糊化。

简单来说，在低光条件下，现有的单光子相机的光通量估计技术无法适用于秒到皮秒级时间尺度的成像。

因此本文开发了一种基于随机微积分（stochastic calculus）和傅里叶级数分解（Fourier series decomposition）的通量探测理论，使用基于统计学和时间序列分析的数学框架，演示了一种将光效和超宽带宽结合在通量函数估计中的新方法。

四、被动超宽带成像基本理论

本文的成像仪器使用单光子雪崩二极管SPAD，它可以以极高的时间精度标记单个光子的到达时间。但是具有以下非理想性质：

1. 量子效率（quantum efficiency）：SPAD一像素检测到的光子的频率/实际光子频率。根据波长不同，量子效率可能远低于1。
2. 死区时间（dead time）：即SPAD的最小的光子检测间隔。如果光子密集到来，落在死区时间则不能被检测到。
3. 时间戳量化（timestamp quantization）：SPAD 具有有限的时间戳分辨率（通常为几皮秒范围内），这会导致时间戳量化。
4. 抖动（jitter）：SPAD 由于计时电子的不稳定性，导致探测光子具有时间延迟变化。

而本文的研究具有以下特征：

• 无源无同步成像：无法改变成像场景，没有电子定时信号来同步。
• 通量函数具有时变性：将一个像素的入射光表示为未知的时变函数 $\varphi (t)$，单位为光子数每秒。是一个连续函数，频率以 $f_{min}$ 和 $f_{max}$ 为界。
• 超宽频通量（Ultra-wideband flux）：在不事先了解成像对象频谱的情况下，重建跨越宽频率范围的通量函数。
• 非齐次泊松光子到达模型：单光子成像中光子到达之间时间跨度不可忽略，故 $\varphi (t)$ 是控制光子到达的非齐次泊松过程的速率函数：
  • $\overline{\varphi (t)}$：平均值表示观测间隔 $[0,t]$ 内的平均通量（单位：光子/秒）
  • $Tavg=1/\overline{\varphi (t)}$：连续光子到达的平均时间跨度（单位：秒）
• 低通量光子检测模型：由于SPAD具有量子效率和死区时间的问题，所以本文设定的环境是低通量光照模型，即 $Tavg>>T_{Dead}$
• 绝对时间戳检测流（stream of absolute detection timestamps）：由于不使用外部同步定时信号，设定光子探测时间戳会根据SPAD的内部时钟单调增加，从而产生绝对时间戳流。（而不是以外部同步频率为时间度量，即“相对”时间戳）

传统技术

1 - 被动光子间成像（Passive inter-photon imaging）

顾名思义就是只检测到个别光子，而假设光子间的通量是恒定的。上图倒数第二行中轴上标记的就是光子探测的时间，最后一行图像呈现了“光子间”的通量假设是相同的。

由于存在“时间跨度内通量没有变化”的假设，意味着该时间跨度内的光子探测速率是通量最大频率的上限。如果在考虑的时间跨度内通量变化显著，则光子探测速率可能无法准确反映通量的最大频率。因此，采用这种方式的SPAD仅限于慢速成像（ $f_{max}=x\ast kHz$ ）。

该种成像方式没有同步信号，因此使用的是绝对时间戳。

2 - 主动直方图成像（Active histogram-based imaging）

使用了已知的同步信号 $f_{sync}$ 进行成像，采用相对时间戳检测。同步频率通常在低MHz范围内，在信噪比与光子丢失或延迟的风险之间取得平衡。

SPAD的检测分辨率构成一个个时间箱（bin），而 $f_{sync}$ 构成一个个检测周期，一个周期一般含有固定个数个时间箱。经过一次次检测周期，光子在时间箱中叠加。

因此可以适用于高速成像（和被动光子间成像相反），但是一个个检测周期时间没有形成关系，所以不适用与低频成像。

但是要求入射通量必须是一个周期函数满足 $T=1/f_{sync}$，如果来自以非倍数 $f_{sync}$ 发射的光子会导致直方图伪影和额外的光子噪声，如上图${\varphi }_2$中的5Hz和${\varphi }_3$中的7Hz的成像结果。

提出新技术：被动超宽带成像机制（The passive ultra-wideband imaging regime）

上面的被动单光子间成像和主动直方图成像都有成像频率的限制，可以说主动直方图法不能检测低频的原因是“同步信号”的存在，导致大时间跨度被分割为若干个不相干的检测周期。

所以本文考虑其极限情况——同步频率降为0，即检测周期就是整个采样的时间跨度。这会导致由于没有周期之间的叠加，每一个时间箱中的光子只可能是0个或者1个，时间戳（由于没有同步信号）变为绝对，每个光子都落在其特有的时间箱中：

至关重要的事，根据奈奎斯特定理所定义的极限（$f_{max}\le 1/(2Q)$），这种方法对极限内所有频率（$DC\to 1/(2Q)$）和任何光源都有重构潜能。

至此，由于没有直方图的协助，任务的主要难点在于：根据全0/1的绝对时间戳构建通量函数的表达式。

五、通量函数的构建

光子计数方程

由于需要在“一个像素上检测到的绝对时间戳流”和“产生他们的通量函数”之间产生联系，所以首先要做的就是构建一个光子计数的理论框架 $N(t)$，表示直到时间 $t$ 之前的光子探测次数：

• 第一个项是通量 $\varphi (u)$ 的积分，直到时间 $t$，它表示到该点之前探测到的光总量。
• 第二个项是鞅 $M(t)$，是一个马丁格尔噪声项，它表示光子探测中的随机波动。

这个累加结果并不尽如人意，高度不连续会导致 $N(t)$ 的傅立叶域中产生很多密集且虚假的噪声。但是这可以由后续的CFAR进行过滤。

通量探测方程

$$p(\tau )=⟨p,\varphi ⟩+M_p(t)$$

方程中 $p$ 和 $\varphi$ 应用函数内积运算，衡量两者相关程度。

$p(\tau )$ 是探测的衡量度，对探测函数在绝对时间戳上求和：
$$p(\tau )=\sum _{t\in \tau }p(t)$$
这个衡量度近似服从正态分布，均值为 $⟨p,\varphi ⟩$，方差为 $⟨p^2,\varphi ⟩$

傅立叶基函数

$$p_f(t)=e^{-2\pi ift}$$

傅立叶探测值近似服从复正态分布（complex normal distribution），均值和协方差矩阵如下：
$$\mu =\left[\begin{array}{cccc}& ⟨\cos (2\pi ft),\varphi (t)⟩& ⟨-\sin (2\pi ft),\varphi (t)⟩& \end{array}\right]$$

$$\sum =\left[\begin{array}{cccc}& ⟨{\cos }^2(2\pi ft),\varphi (t)⟩& 0& \\ & 0& ⟨{\sin }^2(2\pi ft),\varphi (t)⟩& \end{array}\right]$$

可以使用傅里叶探测测量来估计特定通量频率的精度。借助傅立叶探测的归一化能量（也就是傅立叶结果）：
$$p_f^{ϵ}(\tau )=Re[p_f(\tau )/\sqrt{\sum 1,1}{]}^2+Im[p_f(\tau )/\sqrt{\sum 2,2}{]}^2$$
遵循非中心卡方分布。

CFAR频率预警

使用恒定误报率（CFAR）探测器根据所需的误报概率识别和消除噪声频率：
$$|p_f(\tau ){|}^2\ge CD{F}_{{\chi }^2}^{-1}(1-\alpha )N(t)/(2t^2)$$
探测到频率的概率与探测到的光子总数成正比，如果特定频率在通量函数中占主导地位，则探测到频率的概率会降低。

算法流程

关于以上数学计算的用途，可以总结为：

光子计数函数提供了对时变通量的离散和噪声测量，通量探测方程用于重建光子探测时间戳中的通量，并使用CFAR从重建的通量中去除噪声。

而计算通量探测方程均值的作用是：

通过在通量探测方程中加入平均值，即使在低光子数或不可忽视的死区时间条件下，重建算法也可以准确估计一段时间内的平均通量。在统计分析和噪声建模中也起着作用，允许使用CFAR等技术识别和消除噪声频率。

算法实现：

扫描超宽频率带，通过傅立叶重建求出突出频率，再使用CFAR进行阈值判断，分离噪声，最后得出分离的有效频率。

六、成像和实验结果

实验1：多频率分离和重建

左下角为实验环境，环境内有多个不同频率的光源，SPAD自由运动成像后经过傅立叶域重建，可以得出四种有效频率（中间栏和右下角），通过对投影仪频率的图像成像，可以得到右上角结果，使用少量光子重建的成像效果和真实效果不相上下。

实验2：高速风扇捕捉

左边三张图片为新技术实现，1kfps可以捕捉风扇叶片（低频），250Gfps可以捕捉到光的运动（中间两张对比，高频），且两者都无需同步。而传统的主动直方图成像方法无法恢复未同步的风扇旋转。

实验3：多光源场景中分离成像和非视距视频采集

第一行左边展示多光源实验装置，右边展示在具有多个光源的场景中非视距视频成像的设置。第二行图像中 B、C 和 D 点的三个峰值对应于进入瓶子 (B)、传播到瓶盖 © 并向后反射的光脉冲 (D)。

实验结果：

• 展示明了一种前所未有的能力，即能够在没有同步的情况下使用多个光源对场景进行成像。
• 展示了被动 NLOS 视频采集，它允许对由于遮挡或反射而无法直接看到的物体进行成像。

实验4：被动非视距视频采集

使用光栅扫描激光投影仪以被动非视线（NLOS）方式采集场景图像。重建的图像是指使用研究人员开发的通量探测理论和傅里叶域通量重建算法重建的场景图像，可以展示出即使只有2000～3000个光子被捕捉（非视距+低通量），依然可以实现很好的重构效果。

实验5：SPAD阵列技术对比

SPAD阵列有两个常用方法，左边展示是重建每像素通量，右边是分段恒定处理（假设每个像素探测到的光量在一定时间内保持恒定）。

实验结果：

• 即使对于简单的场景，分段恒定通量假设也不成立。
• 每个像素检测到的光量可能会随着时间的推移而变化，在重建图像时需要考虑这种变化。