Research on penetrating imaging LIDAR based on time-correlated single photon counting (invited)

激光雷达 (Light Detection and Ranging, LIDAR) 是一种重要的主动三维光学成像系统。相比于其他探测手段，激光雷达具有空间分辨率高、抗干扰能力强、集成体积小^[1]等特点，因此在军事战场侦察^[2-3]、大气环境监测^[4-5]、空间交会对接^[6]等领域具有极为重要的应用价值。近年来，随着汽车障碍物探测^[7]与智能驾驶技术^[8]的快速发展，人们对能够快速且精确获取路况信息的传感器技术的需求日益增加。由于激光雷达系统具有高距离/角度分辨率、高信噪比、易集成等特点，可精确获得目标的三维图像^[9-10]，可以完成行驶前方道路识别^[11-12]等任务。因此，其已成为车载视觉领域的热点研究方向。

然而，对于较为复杂的实际环境，如：雾霾或充斥烟尘颗粒的大气，浅海或污染较重的水域等环境，激光雷达仍存在难以解决的一系列问题，包括目标探测距离减小、分辨率和信噪比降低^[13]等。上述复杂环境中激光传播路径存在大量尺寸与激光波长相当的散射颗粒，该尺度下的光散射过程服从Mie散射传播规律。多次Mie散射过程导致激光脉冲在时间和空间维度上产生显著的展宽，脉冲能量也会大幅降低，同时散射介质本身会形成强烈的后向散射。尽管介质后向散射光中携带了大量与散射介质相关的信息，但其会对目标信号的探测形成严重的干扰，除去与大气水云参数^[14]相关的特殊研究领域，其余领域往往通过各种办法将其去除，以实现提高目标信噪比的目的。

目前，提升散射介质中光学传感器穿透能力常用的处理方法包括：被动偏振^[15-16]、直方图均衡^[17-18]、距离选通^[19-20]等，但这些方法适用的环境光学厚度(Optical Thickness, OT) 相对较小。对于OT较大的环境，采用单光子雪崩二极管 (Single Photon Avalanche Diode, SPAD)探测器作为接收装置，统计散射介质中光子飞行时间(Time of Flight, ToF )的时间相关单光子计数(Time-Correlated Single Photon Counting, TCSPC)系统，可显著提高激光雷达系统的穿透深度。为此文中对基于TCSPC的穿透成像激光雷达系统进行了研究，并按照如下步骤展开：

(1) 对比国内外在光学成像领域提升穿透深度所采用的方法，分析从小尺度散射介质（生物组织、毛玻璃）到大尺度散射介质（浑浊的大气与水体）中采用的一般方法。并给出TCSPC技术的原理及优越性。

(2) 分析介质中光传输的规律，并建立散射介质中光场传输的基本模型。优化仿真中的部分流程，得到宽范围适用的散射介质穿透成像仿真模型。借此获取浑浊大气的目标时空信号。

(3) 通过分析仿真模型求解的时域信号，文中提出了一套基于TCSPC穿透激光雷达信号的自适应时域滤波算法。并针对该算法详尽分析了可行性及目标提取效果。

(4) 实验验证以上算法对真实环境的鲁棒性，对比传统相机与该实验系统的成像效果，验证气体环境下基于TCSPC的穿透成像激光雷达系统对散射介质穿透深度的提升效果。

穿透介质成像一直是光学成像领域的一大难题。如何在强烈的散射光背景下提取出微弱的目标信号并重建出高质量的图像是该领域的研究焦点问题。

小尺度探测上，近年来研究焦点主要集中在记忆效应^[21-22]、自相关重建^[21,23]、点扩散函数^[24]、散斑场调控^[25-26]、机器学习^[27]等。但相较于毛玻璃、生物组织等小尺度散射体，浑浊的大气和水域不具备稳定的散射颗粒分布，即介质的相干长度(Coherence Length)为含时的随机函数。因此，应用于小尺度散射介质的重建方法直接移植到大尺度散射场景时效果不够理想。探测尺度增大时，恢复算法需要相应调整。

在大尺度探测领域，传统的研究方法包括：主动偏振照明^[15-16]、直方图均衡化^[17-18]等方法，此类方法对于减少雾中目标图像光晕、重现目标细节、提高图像对比度具有一定的可行性。然而，随着环境OT的进一步提高（目标距离增加或散射介质浑浊程度提高），基于目标强度重建目标的方法也逐渐失效。此时目标返回的光子与介质后向散射光子的能量处于同一量级，需升级探测方式才能获取更为理想的目标信号。一些学者在初期探索中采用基于时间选通的激光雷达系统进行探测^[19-20]，实现了约5个衰减长度下的目标物体的图像恢复。文中探讨的基于TCSPC的穿透成像激光雷达系统在一定程度上提高了探测距离。2013年，赫瑞-瓦特大学的McCarthy等人利用InGaAs/InP SPAD探测器，详细分析并搭建了一套1550 nm波段的TCSPC穿透成像激光雷达系统。获得了普通大气环境下，模特目标325 m与910 m的深度图像，及简单多边形目标4.5 km处的深度图像。通过分析不同目标材质与单像素曝光时间对系统空间分辨率的影响，最终获取了低特征目标厘米级分辨率的三维深度图像^[28]。2015年，该团队给出了一种实验室模拟浑浊水体的方式，并利用与参考文献[28]类似的TCSPC激光雷达系统，实验采集了绿、橙、红三种波段下浑浊水体中塑料管目标的TCSPC数据。之后利用像素级的互相关算法清晰还原了1.7 m作用距离、8个衰减长度的光学距离、毫秒级曝光下的目标深度图像^[29]。2018年，麻省理工学院的Satat等人用SPAD相机采集了散射介质下隐藏目标的回波信号。基于时间直方图数据，提出分别拟合散射介质后向散射光与目标弹道光或蛇形光的函数，据此提高了采集信号的信噪比。该团队将其文中系统所得结果与传统相机、光子计数、时间选通的方法进行了比较，认为基于TCSPC系统+算法的采集重建方式可实现更好的探测效果^[30]。2019年，法德圣路易斯研究所的Tobin等人对比讨论了TCSPC的激光雷达系统框架中常用的重建算法如：PwCc、RDI-TV、UA、M-NR3D^[31]，凸显重建算法对目标重建结果的影响。

综上所述，TCSPC激光雷达已成为穿透成像激光雷达领域的研究热点。其具备单光子级别的灵敏度与皮秒量级的时间分辨率，可获得包含目标和散射介质大量信息的波形信号。而针对该波形信号，设计高效的信号提取算法将进一步提升目标重建的质量。鉴于此，文中将围绕时间相关单光子计数的穿透成像激光雷达体制及恢复算法加以展开。

为方便分析散射介质中光传输的规律，文中首先建立散射介质中光传输的基本模型——基于Monte Carlo的Mie散射模型。Mie散射理论是由德国物理学家Gustav Mie于1908年提出的，其主要为解决光传播路径存在尺寸与光波长相当的散射颗粒时，入射光经散射颗粒散射后的散射光场分布^[32]。下文以普适的单球模型为例，对Mie散射理论进行介绍。

如图1所示，设一束光强为${I_0}$的自然光入射至半径为$r$，相对复折射率为$m = {m_r} - {\text{i}}{m_i}$（${m_r}$、${m_i}$分别代表颗粒复折射率的实部与虚部）的各向同性散射颗粒，在距散射体$d$处，散射角为$\theta $（由散射光方向矢量与入射光方向矢量间的夹角确定）的位置，散射光光强$I$（图1中红色部分）为：

Figure 1.  Light field distribution after Mie scattering process

公式（1）中散射振幅函数${S_1}\left( \theta \right)$和${S_2}\left( \theta \right)$具体形式如下：

除此之外，Mie散射理论还给出了消光系数${k_{{\text{ext}}}}$、散射系数${k_{{\text{sca}}}}$、吸收系数${k_{{\text{abs}}}}$等表征散射前后光场能量变化的参数：

式中：无因次直径($\alpha = {{{\text{π }}r} \mathord{\left/ {\vphantom {{{\text{π }}r} \lambda }} \right. } \lambda }$)由粒径$r$和入射光波长$\lambda $确定；Mie系数中包含的${a_n}$和${b_n}$是由Ricatti-Bessel函数定义的，与散射颗粒的复折射率$m$有关。${\pi _n}$与${\tau _n}$则用来表征颗粒的散射强度，是关于$\cos \theta $的勒让德函数和一阶缔合勒让德函数。

Ricatti-Bessel函数一般按无穷级数的形式给出，这给数值计算Mie散射系数带来了一定的困难。不仅如此，公式(2)、(3)中${a_n}$和${b_n}$包含的复杂递推关系会导致Mie散射系数在计算时产生剧烈震荡^[33]。考虑到Ricatti-Bessel函数会随项数增加快速收敛，因此不断有学者提出Ricatti-Bessel函数适当截断项经验公式^[34-35]以及快速计算Mie散射系数的优化方法^[36]。文中采用向前递推与向后递推相结合的方式，权衡计算耗时与精度，建立了一套宽范围高精度Mie散射系数的计算方法。部分计算结果如下所示：

表1中，计算Mie散射系数的平均用时约为$41.812\;{\text{ms}}$。文中算法相比Dave、Bohren算法的精度更高、适用范围更广，相比Lentz算法的耗时更低。验证了文中算法具备宽范围高精度的Mie散射系数计算能力，为后续散射介质传输模型的建立奠定了基础。

对于大气环境，散射颗粒的浓度大致为0.5~50 mg/L^[39]。在此类介质中传输时，绝大多数的光会被多次散射（即散射光），仅有极少的光子未被散射或只被散射了几次（即弹道光与蛇形光）。由于散射传输中存在大量随机且无序的散射事件，仅用单次散射光场分布模型已无法精确描述。基于此，文中选择Monte Carlo方法来模拟散射介质中的光传输过程。如图2所示仿真过程，首先引入光子团的概念，光子团是指：一定数量状态完全相同的光子集合，且该集合内光子传播的行为当且仅当与散射颗粒作用时才会产生差异^[40]。基于图1所示的单次散射模型，计算单次反照率（光子团由散射引起的能量损失率）${\omega _0}$，以及散射自由程$l$、散射角$\theta $、方位角$\;\beta $的概率分布^[36,41-42]，便可利用Monte Carlo方法来模拟散射介质中的光传输过程。

Figure 2.  Simulation procedure of Mie scattering model for light field transmission based on Monte Carlo^[43]

由于穿透激光雷达系统旨在获取散射介质中目标的时空信息，文中对Mie散射模型做出了改进，在散射介质中置入待测目标，并精确统计光子团的飞行时间。过程中需计算光子团与目标的作用坐标，即光子飞行路径所在直线的参数方程与目标法向表面微元的交点，接下来按照光子团入射坐标及方向余弦抽样，光子脉冲抽样，光子自由程抽样，光子团散射角、方位角抽样，最后按照目标碰撞假定、光子团能量、光子团飞行路程累加的流程完成Monte Carlo模拟。

由于该仿真模型中所用抽样的精度依赖于Mie散射系数的计算能力，因而仿真模型同样具有宽范围散射介质模拟、高仿真精度的特点。借此仿真模型，可对不同种类，不同强度的散射介质；不同深度，不同材质的目标；不同范围，不同空间采样率的激光雷达进行探测模拟，并获取目标的强度图、深度图、时间直方图，对后续研究具有显著的指导作用。

明确光子在散射介质中的传输方式后，利用上节建立的仿真模型，对基于TCSPC的穿透成像激光雷达系统进行了探测模拟。仿真中发射与接收信号通过模拟发射具有特定波长、能量、数量的光子团；捕获接收孔径内光子团并记录光子团的能量与飞行时间来实现。探测场景在仿真中由散射颗粒的复折射率、尺寸、单位体积浓度、以及物体的目标函数共同确定。利用该模型模拟的雷达接收信号与目标强度图像如图3所示（仿真规模：光子团数$1 \times {10^7}$）。

Figure 3.  Detective simulation results of TCSPC lidar system. (a) ${\text{OT}} = 2.016\;3$; (b) ${\text{OT}} = 4.032\;5$

图3所示仿真结果预示：若散射介质光学厚度进一步增加（${\text{OT}} \geqslant 4$时），携带目标信息的弹道光在回波信号中的占比将远小于1，此时常规光学成像手段将难以提取混合在后向散射信号中的目标信息。不仅如此，单位距离内散射事件的增多可严重展宽时空点扩散函数，即雷达接收信号的方向性将随介质${\text{OT}}$的增加大幅下降。散射介质这种退化目标强度与方向性的特性，严重影响了激光雷达系统的成像效果。为此，引入与探测手段相匹配的算法将成为进一步提高探测系统穿透深度的关键。

对比图3(a)、(b)得出：即使目标在强度分布的特征难以辨别，但在时间信号相应位置处仍存在突出的目标波形。此时引入目标距离信息，再借助时间选通等方式便能较好的提取目标特征。然而，这种提取手段的前提是散射介质的后向散射与目标峰形不产生混叠，或目标峰远高于临近通道的光子计数。在实际环境中，散射介质${\text{OT}}$的提升往往伴随着目标深度或散射颗粒浓度的增加。若目标深度不变，所处环境的散射颗粒浓度增加，目标峰则可能与介质的后向散射信号产生混叠，如图4所示。

Figure 4.  Simulation results of aliasing between the target peak and the backscattered signal of the medium

对于散射程度较小的大尺度探测环境，常规算法通常能够有效地提取目标信息。但该类算法的恢复效果较依赖于信号中弹道光的信噪比，且通常不考虑散射介质的传输特性，难以应用于散射更强的环境。散射介质中散射颗粒浓度提高会增强介质的后向散射；抑制目标峰的幅值；目标峰最终淹没于后向散射信号难以被采集。对于存在混叠的散射场景，常规算法便难以适用。此时需从散射介质的传输特性入手，以有效提取目标信息。对此，文中提出一种基于信号时域分布的自适应目标提取算法：

首先，设采集信号中介质后向散射满足的分布为${S_{\text{b}}}$，目标峰满足的分布为${S_{\text{o}}}$，则采集信号$S$满足：

式中：${c_{\text{b}}}$、${c_{\text{o}}}$代表叠加权重；$ {s_{\text{n}}}\left( t \right) $表示采集信号中的随机噪声。$ {c_{\text{b}}}{S_{\text{b}}} $由Mie理论给出，散射介质中光子团的散射自由程应满足：

式中：$\bar l$表示当前散射环境的平均自由程。由于碰撞次数相同的光子团在接收信号中应服从相同的分布，且Monte Carlo过程的光子团互相独立。因而碰撞次数相同的光子团的叠加应服从伽马分布，即总计发生$i$次碰撞的接收信号满足时间分布：

式中：$\bar t$为散射环境中两次散射事件的平均时间间隔，满足$\bar t = {{\bar l} \mathord{\left/ {\vphantom {{\bar l} {\text{c}}}} \right. } {{c}}}$（${{c}}$为光速）。

借助Monte Carlo仿真模型，可以统计出接收信号中不同碰撞次数的光子团满足的频数序列（$N = {n_1},{n_2}, \cdots ,{n_i}, \cdots $，式中$i$表示光子团发生碰撞的次数，${n_i}$表示总计发生$i$次碰撞的光子团个数，且$i \to \infty $）。实际过程中的单次反照率令光子团能量将随碰撞次数指数递减，若选择最大碰撞次数$a$作为次数统计的上限，可以最大似然的估计出发生$i$次碰撞的光子团在接收信号中出现的频率为：$ {p_i} = {{{n_i}} \mathord{\left/ {\vphantom {{{n_i}} {\sum\nolimits_{i = 1}^a {{n_i}} }}} \right. } {\sum\nolimits_{i = 1}^a {{n_i}} }} $，那么采集信号中后向散射的分布${S_{\text{b}}}$应为^[30]：

若发射系统的输入功率为${P_{{\text{in}}}}$，则${c_{\text{b}}}{S_{\text{b}}}$为：

式中：${\omega _0}$为散射颗粒的单次散射反照率。

在已知$ {c_{\text{b}}}{S_{\text{b}}} $的基础上，若通过滤波对信号$S$的噪声加以滤除，就能获得目标信号$ {c_{\text{o}}}{S_{\text{o}}} $，例如采取小波变换等方式处理信号中的${s_{\text{n}}}$部分。为考察该算法是否有效，文中首先利用仿真数据进行验证，下图为一组时间直方信号的去噪结果，包含后向散射信号${S_{\text{b}}}$分布统计拟合结果与差分信号的目标峰提取结果：

据公式(7)，对去噪信号与${c_{\text{b}}}{S_{\text{b}}}$进行差分便可得到${c_{\text{o}}}{S_{\text{o}}}$。但由于小波变换等去噪算法不能理想的剔除系统噪声，差分后的结果一般为$ {c_{\text{o}}}{S_{\text{o}}} + \delta \left( {{s_{\text{n}}}} \right) $，这里$ \delta \left( {{s_{\text{n}}}} \right) $表示经由小波变换后残余的噪声。

考虑目标信噪比趋于$1$的情况：目标峰处于后向散射信号拖尾、且强度较低的部分，差分过程对目标峰型引入的变形较小。因$\delta \left( {{s_{\text{n}}}} \right)$与${S_{\text{o}}}$存在分布差异，可采用相关的方法提取${c_{\text{o}}}{S_{\text{o}}}$，并得出目标峰位及标准差，完成目标的提取（图5 (b)）。

Figure 5.  Denoising results of time histogram signal. (a) Distribution fitting result of backscattering signal${S_{\text{b}}}$; (b) Target peak extraction result for differential signal

综上，该算法能够自适应地提取出目标特征，算法平均耗时为$1.023{\text{s}}$，对目标强度图像的信背比（目标强度与背景平均强度的比）与边缘锐度带来较大提升。为进一步验证该算法的适用范围，文中对两种介质的仿真信号进行了特征提取，结果如下：

图6、图7所示为散射介质中激光雷达系统的仿真探测信号的强度图与时间直方图。两组图的仿真参数如下：图6中，目标距离接收孔径$0.2\;{\text{m}}$、散射颗粒不对称因子$0.924$、从图(a)~(c)可知目标处${\text{OT}}$分别为：$4.993\;6$、$6.048\;5$、$6.912\;3$，算法平均耗时为$1.085\;{\text{s}}$；图7中，目标距离接收孔径1 m、散射颗粒不对称因子$0.524$、从图(a)~(c)可知目标处${\text{OT}}$分别为：$1.920\;3$、$3.072\;6$、$3.936\;4$，算法平均耗时为$1.117\;{\text{s}}$。对比两图结果，尽管在${\text{OT}}$较大时该算法恢复目标的强度较低，但仍较好地保留了目标特征。相比于恢复前的强度分布，目标信背比与目标特征的留存度都有了较大的提升。通过仿真数据的验证，证明了该算法的去散射能力与信号峰的信噪比直接相关。该算法的实用性依靠实验进一步验证。

Figure 6.  Feature extraction results for target in medium 1

Figure 7.  Feature extraction results for target in medium 2

为验证真实场景下该目标提取算法的效果，文中搭建了一套结合摆镜扫描的非共焦TCSPC穿透激光雷达系统，该系统的原理图如图8所示。

Figure 8.  Non-confocal TCSPC penetrating lidar system based on galvanometer scanning

实验系统主要包括：超窄脉冲激光器、PIN光电探测器、单像素单光子探测器、扫描振镜、TCSPC模块、以及光路系统。其主要性能参数如下：激光脉冲的波长为$532\;{\text{nm}}$、脉宽为$5\;{\text{ps}}$、单光子计数器时间分辨率为$ {\text{64 ps}} $，对不同种类的散射环境均能实现快速采集。系统的工作状态如下：激光器发出的脉冲光被分束器分为两路，一路作为发射信号对散射场景进行照明。另一路经PIN光电探测器转换为TCSPC模块的起始信号，SPAD探测到回波光子输出电脉冲作为TCSPC模块的停止信号，记录得到的时间差即为光子的飞行时间。相比于传统的激光雷达系统，TCSPC激光雷达具备单光子级别的灵敏度与皮秒量级的时间分辨率，可获得包含目标和散射介质大量细节信息的波形信号，更适合研究通过信号处理来抑制背景、提取信号的方法。文中设计的这套TCSPC激光雷达系统具有发射信号脉宽短、时间分辨率高的特点，有利于在散射较强场景下实现瞬态光场的波形精确提取。同时采用非共焦扫描接收方式，可以有效避免系统内部的散射干扰，有助于在大发射功率下实现对目标回波的快速扫描获取，提高系统空间分辨能力。

实验放置“S”字型目标（目标条宽$5\;{\text{cm}}$，所占区域面积约为30 mm×30 mm）于距探测器空间距离1.31 m处的雾仓内。为方便对比，图9给出了消费级相机对不同浓度散射介质的成像结果。

Figure 9.  Imaging results of different concentrations of scattering media with a consumer-grade camera

实验时雾仓内填充高浓度的气溶胶液滴来模拟浓雾环境。该环境下时域信号滤波与信号拟合的目标提取算法对目标特征的提取结果如图10、图11所示。

Figure 10.  Statistical fitting results of signal distribution in a scattering scene

结合图10 、图11的实验结果（算法平均耗时为$1.118\;{\text{s}}$），证明该算法对于不同强度的散射环境均有一定程度的提取作用。在图11(c)中，该算法的特征提取效果受到了较大影响，此时该算法主要起自适应目标定位的作用。考虑此时采集信号的信背比（目标信号与介质散射信号的峰值比）约为 $0.025$、信噪比约为 $1.209$，文中提出的算法对目标横向特征的提取具有较为显著的效果。

最后给出该系统的景深特性与目标重构算法对不同距离下多目标鲁棒性的测试实验，实验结果如下：

图12 （算法平均耗时为$1.643\;{\text{s}}$）中，$1\;{\text{m}}$ 处的${\text{OT}} \sim 2.655\;7$，两目标的纵向间隔约为 $0.55\;{\text{m}}$。实验结果证明：对雾中纵深较长、且目标信号与后向散射信号混叠严重的复杂目标，文中提出的算法可较为完整地提取出目标特征，针对不同距离下的多目标信号仍有较强的鲁棒性。

Figure 11.  Feature extraction results of target "S". (a) ${\text{OT}} \text{-} 1.541\;6$; (b) ${\text{OT}} \text{-} 2.792\;4$; (c) ${\text{OT}} \text{-} 4.140\;1$

Figure 12.  Results of multi-objective experiments. (a) Scene graph; (b) Received signals and fitting results; (c) The intensity image after algorithm processing (bottom right is the raw intensity image); (d) The depth image after algorithm processing (unit: m)

为克服传统激光雷达系统难以在散射环境下成像的问题。文中建立了散射介质下的瞬态光场传输模型，并以此模型的统计规律提出了一套基于自适应的时域信号滤波与信号拟合的目标提取算法。通过仿真数据与实验数据的验证，证明了该算法能够从信背比为$0.025$、信噪比$1.209$ 的回波信号中提取出目标特征，同时证明了仿真模型的可靠性。文中搭建的非共焦TCSPC穿透成像激光雷达系统结合该滤波算法，可实现${\text{OT}} > 4$的散射环境中目标横向特征的快速提取，以及${\text{OT}} > 2.6$环境下多目标物体的纵向特征的快速提取。该算法对于提高激光雷达在散射介质中的穿透深度具有较强作用，同时为雷达接收信号的时域分析提供了新的思路。