Adaptive spatial-temporal correlation depth estimation of photon-counting lidar

近年来，为了以较低的激光发射功率实现对远距离目标的高灵敏探测，基于盖革雪崩光电二极管（Geiger-Mode Avalanche Photodiodes, GM-APD）和时间相关单光子计数技术（Time-Correlated Single-Photon Counting, TCSPC）的光子计数激光雷达成为研究热点，并在地形测绘、自动驾驶和飞行器制导等精密探测领域得到广泛应用^[1-5]。与传统探测体制相比，单光子探测是一种概率探测^[6]，由于探测环境变化而导致的探测概率的下降可以通过延长探测时间进行弥补。采用TCSPC技术记录回波光子的时间信息，将多个探测周期的回波光子数据进行累加，形成光子计数统计直方图，通过时间数据求解目标的距离信息。光子计数激光雷达具有ps量级的时间分辨率和极高的探测灵敏度，可以对单个光子进行响应，尤其在微弱回波信号探测领域具有极大的研究价值，但也正是由于该特点，使得目标探测极易受环境影响，回波数据中通常会包含大量由背景环境引起的噪声光子^[7-9]，例如：太阳背景辐射、激光后向散射等，为信号光子的高质量提取带来不利影响，制约系统整体探测性能的提升。因此，研究如何在强噪声背景下提高光子计数激光雷达的成像质量十分重要。

根据回波光子数据特点，一些利用回波数据时间相关性的滤波算法被提出^[10-12]，这些方法通过深入分析回波光子的时间响应特性，估计得到目标的深度图像，但是针对缺失信号的情况，数据处理效果有待提升。Kirmani等^[13]提出首信号光子法，利用图像的空间相关性，仅使用系统接收到的第一个光子即可得到目标的深度和反射率信息，但成像效果随噪声增大急剧下降。Chen等^[14-15]提出多尺度时间分辨率结合自适应阈值分割的单光子成像方法，实现了噪声和信号的有效区分。Rapp等^[16]提出基于统计模型的解混叠算法，分别构建信号光子和噪声光子的时域模型，并结合空间关系建立超像素进行滤波，提升了图像重构效果。这些方法虽然都在对回波光子时间相关性分析的基础上，将数据的空间关联作为目标重建的重要因素进行考虑，但是噪声不断增强，数据计算复杂度也随之变大。

为了提高目标距离信息的估计精度，提出一种适用于强背景噪声、大数据量的自适应时空关联深度估计算法。所提算法基于回波数据的时间特性和空间关联性，自适应调整时间窗口大小，实现信号光子的粗提取，并结合滑动窗口和自适应均值滤波进行精去噪，得到各像素的飞行时间。通过仿真和室内成像实验对所提算法的滤波效果进行验证，完成待测目标重建。

当使用GM-APD单光子探测器结合TCSPC技术进行微弱回波信号探测时，回波光子数目通常较少，近似认为服从泊松分布，在时刻t像素（i, j）的总回波光子可以表示为^[17]：

式中：$ s\left( \cdot \right) $为激光速率方程；$ {r}_{i,j} $为目标反射率；${{\textit{z}}}_{i,j}$为待测目标的距离；$ {b}_{i,j} $为噪声的期望值。

根据回波光子在时间轴上的分布特性，即：信号光子呈现集中分布，在激光脉宽范围内光子计数值较多，而噪声光子在整个探测时间上通常呈现均匀分布^[11]。利用回波光子在时域上的分布差异性，当待测目标距离值为Z时，回波光子计数模型如图1所示，红色实线表示信号光子累加结果。

Figure 1.  Echo photon counts model

文中所提算法通过3个步骤完成滤波和距离信息估计，分别是自适应调整时间窗口的大小及数据粗提取、滑动窗口求加权时间、自适应均值滤波，如图2（a）所示。

假设像素（i, j）总回波光子个数为N，飞行时间数据$ {T}_{i,j}=\left\{{t}_{i,j}^{1},{t}_{i,j,}^{2},\cdots ,{t}_{i,j}^{N}\right\} $，原始直方图的时间分辨率为$ \mathrm{\tau } $。根据各像素回波光子的时间分布特点，为了在尽可能多保留光子计数数据的同时降低背景噪声带来的干扰，依照从小到大的顺序设置m个小于激光脉冲宽度$ {w}_{p} $的时间分辨率，即W={$ {w}_{1}\mathrm{\tau },{w}_{2}\mathrm{\tau },\cdots ,{w}_{m}\mathrm{\tau } $}。通过构建具有较低时间分辨率的直方图和借邻域数据信息的方式，自适应调整时间窗口大小，初步确定信号光子所在位置，如图2（b）所示。

当时间分辨率为$ {w}_{1}\mathrm{\tau } $，若重构直方图中存在单一的光子计数最值，则认为是信号响应，记录最值位置。反之，依次循环所设置的时间分辨率进行直方图重构，直至寻找到单一最值，并对得到的各像素最值位置进行中值滤波，筛选剔除异常最值点。此时，最值位置在原始直方图中对应的时间区间即为所求的时间窗口。若遍历所有预先设置的时间分辨率后，仍然未得到最值位置，则认为无法确定是否包含有效信号，利用各像素数据的空间关系建立邻域数据集，基于该像素3×3邻域的时间信息求解时间窗口，如图3所示。若邻域中无信号像素点数目大于5，则扩大邻域范围至5×5，直至满足条件。将各像素时间窗口内的数据作为后续待处理的回波数据，完成信号光子数据的粗提取，极大地降低了后续计算的数据量。

Figure 2.  (a) Overall flow chart; (b) Adaptive selection of the time window

假设窗口大小为$ {W}_{wind}\mathrm{\tau } $，使用该窗口依次滑过像素（i, j）所提取的光子计数数据$ k\in \left\{\mathrm{1,2},\cdots ,n\right\} $，如图3（b）所示，寻找光子计数最多的窗口位置：

求此时窗口内回波光子的加权时间和$ {T}_{i,j}^{weit} $，即为该像素的飞行时间。

判断在像素（i, j）的n×n邻域内（n一般取3和5），$ {T}_{i,j}^{weit} $与邻域中值$ {T}_{i,j}^{med} $的距离是否小于所设置的阈值$ {T}_{i,j}^{thr} $，若满足，$ {T}_{i,j}^{weit} $即为该像素的飞行时间，反之，认为该像素点不存在有效信号。

Figure 3.  Neighborhood data set (3×3). (a) The time interval of the maximum counts of each pixel in the neighborhood; (b) The time interval of the neighborhood data set and the corresponding photon counts in the original histogram

对经由以上步骤得到的各像素飞行时间进行自适应均值滤波，对不存在有效信号光子的像素点，使用其3×3邻域像素的时间均值进行替代，若邻域数据集中不存在有效的时间信息，则继续扩大邻域范围，直至寻找到有效值。最后，基于估计得到的飞行时间$ \hat{T} $和公式$ \hat{Z}=c\hat{T}/2 $，即可解算出相应的距离信息。

对室内静态目标进行成像探测时，只有物体边缘部分存在深度变化，大部分像素相邻位置的深度几乎不变，待测目标的距离信息整体上不会呈现出起伏变化的趋势。为了验证所提算法在强噪声干扰下，针对深度信息连续变化的目标，是否具有较好的重构效果，使用如图4（a）所示的地形作为仿真对象。考虑探测器死时间的影响，采用蒙特卡洛方法仿真服从泊松分布的光子数据，并通过改变噪声通量的大小，生成包含不同噪声强度的回波数据，文献[14]已经验证该仿真数据的准确性，相关仿真参数见表1。

Figure 4.  (a) The original terrain; (b) Photon statistical histogram (M=10000 and 2 MHz noise)

在${t}_{1} {\text{-}} {t}_{2}$探测时间内，若接收端的平均回波光子数为$ \overline{{M}_{sn}} $，包括信号光子$\overline{{M}_{s}}$和噪声光子$\overline{{M}_{n}}$，则GM-APD探测器产生的平均初始光电子数为：

单次脉冲探测概率^[2]为：

式中：$ {N}_{dead} $为死时间内的光子计数。

结合公式（4）、公式（5）和光子统计直方图估计信号光子数目，当噪声强度约为2 MHz、脉冲积累数目为10 000个时，结果如图4（b）所示，此时，每个脉冲的回波光子数（Photons per Pulse, PPP）约为13.8。

在仿真实验中保持信号通量不变，分别添加强度为0.2、2、3.5、6 MHz的噪声，回波数据的PPP水平约为14.3、13.8、14.6和14.1，定量对比在不同噪声强度下，使用峰值法、Chen算法^[14-15]和所提算法分别进行地形重建的效果，如图5所示，其中所提算法预设的时间分辨率为W={$ 2\mathrm{\tau },\mathrm{ }4\mathrm{\tau },\mathrm{ }8\mathrm{\tau } $}，滑动窗口大小为$ 2\mathrm{\tau } $。使用均方误差（Mean Square Error, MSE）作为评价标准，即：

Figure 5.  Comparison of reconstruction results of different algorithms at different noise intensities. (a) 0.2 MHz noise, 14.3 PPP; (b) 2 MHz noise, 13.8 PPP; (c) 3.5 MHz noise, 14.6 PPP; (d) 6 MHz noise, 14.1 PPP

从图5的仿真结果可以看出：当噪声强度在0.2 ~ 6 MHz之间，所提算法均能重建出地形信息，随着噪声强度增加，均方误差由0.04 增加至0.51 ，较其他方法进行地形重建的均方误差最小，尤其在强噪声数据中展现出较好的效果。当噪声在0.2~3.5 MHz之间时，峰值法重构的均方误差从0.11 增大至1.30 ，且在同一噪声强度下，其重构效果较其他方法略差，当噪声高于3.5 MHz时，失去重建能力。Chen算法在噪声强度低于2 MHz时，与文中算法的处理效果基本相近，但对于强噪声数据，处理效果急剧下降，在噪声强度高于3 MHz时无法准确提取出地形高度信息。在噪声强度相同时，所提算法较峰值法和Chen算法进行地形重建的均方误差降低了至少20%。仿真实验表明，所提算法针对深度信息起伏变化的目标具有良好的重建效果。

实验系统如图6所示，激光的中心波长为1064 nm，脉宽为3.5 ns，重频为6.5 kHz。发射激光脉冲通过分束器后，少部分光作为同步信号，大部分光到达光学发射系统后，通过控制快反镜的旋转角度改变激光传播方向，在目标区域形成激光点阵，完成目标扫描。在接收端，回波激光和背景杂散光激发单光子探测器的响应，光子探测效率为2.8%。通过TCSPC模块记录回波光子的飞行时间，最小时间分辨率τ为64 ps。

对室内静态目标进行成像探测，待测目标由不同距离处的2个纸盒组成，如图7（a）所示。在暗室环境中，设置各像素的采集时间为200 ms，图像大小为60 pixel×60 pixel像素，使用峰值法和中值滤波法对回波数据进行处理，得到目标深度重建的真值图像，如图7（b）所示。

保持激光发射功率不变，在接收端引入波长为1064 nm的连续激光光源，通过调整连续激光的光功率来模拟不同强度的噪声环境。设置各像素的脉冲积累个数为20，采集包含不同噪声强度的回波光子飞行时间数据，噪声强度估计值分别为0.23 MHz、1.21 MHz、2.17 MHz、3.02 MHz和5.08 MHz。使用不同算法对以上数据进行处理，结果如图8所示。其中，所提算法设置的重构直方图时间分辨率W={$ 20\mathrm{\tau }, \mathrm{ }30\mathrm{\tau },\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }36\mathrm{\tau },40\mathrm{\tau },\mathrm{ }50\mathrm{\tau } $}，滑动窗口大小为$ 2\mathrm{\tau } $。

Figure 6.  (a) Block diagram of photon counting lidar system; (b) Physical picture of device

Figure 7.  (a) Visible band image; (b) Reconstructed truth depth image of target

Figure 8.  Reconstruction results of target. (a) 0.23 MHz noise, 21.6 PPP; (b) 1.21 MHz noise, 21.2 PPP; (c) 2.17 MHz noise, 21.8 PPP; (d) 3.02 MHz noise, 22.5 PPP; (e) 5.08 MHz noise, 23.9 PPP

从图8可看出，所提算法在0.23 MHz、1.21 MHz、 2.17 MHz、3.02 MHz和5.08 MHz噪声强度下，重建的均方误差分别为0.009、0.011、0.013、0.017 和0.017。从图8可以得出，在噪声强度较低时，采用峰值法可以恢复出目标的部分深度信息，但是针对强噪声数据，峰值法失去效果。在0.23 ~ 5.08 MHz噪声强度范围内，Chen算法和所提算法均能够完成目标重建，前者的均方误差由0.010 增加至5.74 ，后者的均方误差由0.009 增加至0.017 ，使用Chen算法重建质量下降较为明显。随着噪声强度逐渐增加，峰值法几乎无法提取出目标信息，Chen算法重建得到的目标深度图中，噪声点数明显增多。在同种噪声强度下，使用所提算法较其他两种方法进行目标重建的效果更好，重建的均方误差至少降低了约91%，进一步验证了所提算法对强噪声数据的鲁棒性。

文中提出一种光子计数激光雷达的自适应时空关联深度估计算法，通过重构具有较低时间分辨率的直方图，并结合邻域像素的时间信息，自适应寻找信号光子所在的时间窗口，对窗口内的数据进行滑窗后，结合自适应均值滤波，完成距离信息估计。既解决了传统峰值法存在多个最值的问题，又极大降低了待处理的数据量。针对起伏变化的地形和室内静态目标成像探测，所提算法较峰值法、Chen算法具有更好的目标重建效果，均方误差至少降低了约20%。未来的工作将在此基础之上，研究稀疏回波光子数据中的目标重建。