星载光子计数激光测距雷达的实时去噪方法

谭崇涛; 于文博; 向雨琰; 李少辉; 余婧; 王倩莹; 李松; TAN Chong-Tao; YU Wen-Bo; XIANG Yu-Yan; LI Shao-Hui; YU Jing; WANG Qian-Ying; LI Song

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

星载光子计数激光测距雷达的实时去噪方法 PDF

- ORCID：
谭崇涛 ¹
✉
- ORCID：
于文博 ¹
- ORCID：
向雨琰 ¹
- ORCID：
李少辉 ²
- ORCID：
余婧 ²
- ORCID：
王倩莹 ²
- ORCID：
李松 ^1,3
✉

1. 武汉大学电子信息学院，湖北武汉 430072； 2. 中国空间技术研究院，北京 100098； 3. 武汉量子技术研究院，湖北武汉 430010

中图分类号： TN958.98

最近更新：2024-05-16

DOI：10.11972/j.issn.1001-9014.2024.02.014

摘要

星载光子计数体制激光测距雷达系统具有高重频、高精度等显著优势，但也面临原始数据量大且噪声数据占比过高的问题。为适应星上数据通道的传输能力，需压缩原始数据量并保障信号光子的查全率，因此必须发展以硬件为主体的实时去噪算法。本文提出一种粗精结合的快速去噪算法，首先基于激光器发射脉宽、系统噪声率、目标特性以及接收光子事件的局部密度信息进行粗去噪，剔除部分噪声光子；再利用直方图统计，对保留的光子事件进行精去噪，确定信号光子区间及最终的信号光子及其时间信息。通过蒙特卡洛仿真和ICESat-2实测数据对算法进行验证，测试结果表明，本算法查全率大于94%、查准率大于93%、调和平均值大于94%，运行效率提高了10%。算法可以实现光子事件的快速实时去噪，为星上硬件实时去噪处理提供了理论基础。

关键词

光子计数; 激光测距; 粗精去噪; 数据密度; 直方图统计

引言

光子计数体制星载激光测距雷达具有高重频、高精度等显著优势，但也面临原始数据量巨大且无用的噪声数据占比过高的问题。为保障对信号光子的查全率，同时压缩原始数据量，以适应星上数据通道的传输能力，必须发展以硬件为主体的实时去噪滤波算法。

光子计数体制的激光测距雷达原始数据去噪主要在光子点云数据的二维时空剖面上进行，通过时空数据密度对噪声和信号进行区分^［

1-7］，算法相对比较复杂，适合在上位机上进行执行。快速实时的信号处理需要在更底层的一维时域上进行去噪，如美国2018年发射的ICESat-2（Ice， Cloud and Elevation Satellite-2）卫星上搭载的光子计数激光雷达系统ATLAS（Advanced Topographic Laser Altimeter System），在光子计数电子卡PCE（Photon Counting Electronics）的FPGA（Field Programmable Gate Array）内部进行直方图统计，依据特定的时间窗口内的光子个数即光子事件的时域密度进行了实时的信号提取和噪声去除^{［参考文献 8-15}8-15］。统计直方图就是在一定的时间窗口范围内，以细分的时间bin宽度作为时间分辨率，落在每个时间bin内的光子事件的数量累计。当统计直方图上出现明显峰值，即可判定信号所在的时间bin，据此对信号和噪声进行区分。统计直方图一般采用固定的时间窗口和时间分辨率。

另外一种时域去噪方法则是依据特定个数光子事件的时间分布范围进行噪声滤除^［

16-18］。该方法本质上与直方图统计一样，都是通过光子时间数据的密度信息对信号光子和噪声光子进行区分。但需要依据系统噪声率和目标反射率设定统计的光子时间数据个数，判定光子时间分布范围是否在一定的阈值范围内。由于信号相对噪声更加密集，信号分布的时间范围比较小，据此可以对信号和噪声进行区分。该方法可以根据目标反射特性灵活调整时间窗口大小，可自适应实现不同信噪比条件下的噪声滤除。

2016年，Rapp J 和Goyal V K对光子计数激光雷达信号和噪声光子的探测概率和虚警概率进行了详细推导^［

16］，得到噪声光子时间数据序列的次序统计量满足β分布的结论。通过设定时间窗口长度，和窗口内满足一定噪声虚警概率的最小光子数的阈值，对噪声进行判定和去除。该算法仍然需要预先确定时间窗口，不利于实时处理。2018年，骆乐、吴长强等人^{［参考文献 17

百度学术}17］利用信号光子和噪声光子在时间轴上不同的分布特性，结合泊松分布的统计规律，通过累积设定的邻域范围内每对相邻光子的时间数据差值，与激光发射脉冲脉宽进行比较，来判定是否为信号光子，并以信号光子时间数据的平均值作为最终到达结果。该方法不需要设定时间窗口大小。2022年，朱文华、汪书潮等人提出一种基于自动选取单像素最佳累积时间的时间自适应扫描方法^{［参考文献 18

百度学术}18］，针对不同反射率目标的实验结果表明，相比单像素固定累积时间的扫描方法，该方法的扫描时间降低了一个数量级，避免了数据累积冗余或不足，一定程度上提高了数据获取效率。

本文提出一种粗精结合的快速去噪算法，首先基于激光器的发射脉宽、系统噪声率、目标反射特性以及接收光子事件的局部密度信息进行粗去噪，剔除部分噪声光子；再利用直方图统计，对保留的光子事件进行精去噪。文章对去噪算法模型进行详细介绍，并通过蒙特卡洛仿真和ICESat-2实测数据对算法进行验证。测试结果表明，本算法可以实现光子事件的快速实时去噪，为星上硬件实时去噪处理提供了理论基础。

1 光子计数体制激光雷达探测概率模型

星载光子计数体制激光雷达采用脉冲飞行时间测距技术^［

19-21］进行测距，其探测过程如图1所示。星载激光雷达以一定重频对地发射激光脉冲，经大气传输后入射到目标表面，经由目标反射、二次大气传输后被激光雷达接收系统接收，激光雷达系统同时还会接收到来自大气散射的太阳背景噪声和目标反射的太阳背景噪声。

图1 星载光子计数激光雷达探测过程示意图

Fig. 1 Schematic diagram of the detection process of spaceborne photon counting LiDAR

通常情况下，激光雷达发射系统的激光脉冲的初始分布在时域和空域上可以用高斯分布近似描述。根据Gardner的激光雷达直接探测的相关理论，单程大气传输相当于一次夫琅禾费衍射，高斯脉冲经夫琅禾费衍射后不改变其分布形式^［

22］。对于平面漫反射目标，高斯脉冲经其反射后同样不改变分布形式，因此可以用高斯函数对激光雷达接收回波信号进行近似描述^{［参考文献 23-24}23-24］：

S (t) = \frac{N_{s}}{\sqrt[]{2 π} σ_{s}} e x p [- \frac{{(t - T_{s})}^{2}}{2 σ_{s}^{2}}]

(1)

式中， $N_{s}$ 表示平均信号光子数，可以用激光雷达方程计算， $σ_{s}$ 表示回波信号的均方根脉宽。对于远距离激光测距，当激光入射方向与目标表面法线方向平行时，回波信号的均方根脉宽 $σ_{s}$ 可以用发射的激光脉冲的均方根脉宽 $σ_{t}$ 近似表达，即 $σ_{s} \approx σ_{t}$ 。 $T_{s}$ 为回波信号的时间重心。

光子计数激光雷达的噪声主要包括背景光噪声和单光子探测器本身的暗计数，系统总的噪声光电子数可以表示为^［

25］：

N_{n} = N_{b} + N_{d}

(2)

式中， $N_{n}$ 是总的噪声光电子数， $N_{b}$ 是背景光噪声函数的平均光电子数， $N_{d}$ 是探测器暗计数。

系统探测到的光子个数近似服从泊松分布，在极短时间片 $τ$ 内，单光子探测器探测到 $k$ 个光电子事件的概率为：

P (k, τ) = \frac{N_{τ}^{k} e x p (- N_{τ})}{k!}

(3)

式中， $N_{τ}$ 是极短时间片 $τ$ 的平均光子事件个数。

根据式（3），在时间片 $τ$ 内没有探测到任何光子信号的概率为：

P (k = 0, τ) = e x p (- N_{τ})

(4)

在时间片 $τ$ 内探测到信号的概率为：

P (τ) = 1 - e x p (- N_{τ})

(5)

因此，探测到噪声光子的概率可以表示为：

P_{n} (τ) = 1 - e x p (- N_{τ n})

，

(6)

式中， $N_{τ n}$ 为极短时间片 $τ$ 内的噪声光子数。

回波信号近似为高斯分布，根据高斯分布的“ $3 σ$ ”原则，信号光子在区间 $(T_{s} - 3 σ_{s}, T_{s} + 3 σ_{s})$ 内出现的概率为0.9974，可以作为信号光子所在的区间，而其他的区间近似为只包含噪声。在区间 $(T_{s} - 3 σ_{s}, T_{s} + 3 σ_{s})$ 内信号光子探测概率为：

P_{s} (6 σ_{s}) = 1 - e x p (- N_{s} - 6 σ_{s} f_{n})

(7)

在 $6 σ_{s}$ 脉宽范围内噪声光子的探测概率，也即虚警概率为：

P_{n} (6 σ_{s}) = 1 - e x p (- 6 σ_{s} f_{n})

(8)

式（7）和（8）中的 $N_{s}$ 为 $6 σ_{s}$ 时间范围内的信号和噪声光子事件数， $f_{n}$ 为噪声率。定义 $6 σ_{s}$ 时间范围内的信号光子和噪声光子探测概率的对比度为：

\begin{matrix} C_{6 σ_{s}} = \frac{P_{s} (6 σ_{s})}{P_{n} (6 σ_{s})} = \frac{1 - e x p (- N_{s} - N_{n})}{1 - e x p (- N_{n})} \\ = \frac{1 - e x p (- N_{s} - 6 σ_{s} f_{n})}{1 - e x p (- 6 σ_{s} f_{n})} > 1 \end{matrix}

(9)

在区间 $(T_{s} - 3 σ_{s}, T_{s} + 3 σ_{s})$ 上，以 $6 σ_{s}$ 为时间片大小，目标点的探测概率应大于噪声点。当测量频率较高且目标轮廓表面在小区域内具有一定的连续性和一致性时，经过多次测量，目标点的分布更加集中且具有一定相关性，而噪声点分布相对稀疏，可近似为均匀随机分布，利用目标点和噪声点的这一差异性，可以实现光子计数激光雷达信号光子的提取和噪声去除。

2 去噪算法流程

太阳光产生的背景噪声是星载光子计数激光雷达的主要噪声来源。其不仅对信号探测产生了极大的干扰，无效的噪声还占用了大量的硬件系统资源，影响了对信号的处理能力。采用粗精结合的去噪算法，利用光子事件时间数据局部密度信息，先进行粗去噪，快速筛选出部分噪声光子并剔除，将数据量减少到系统可承受的程度；再进行精去噪，基于对保留的光子事件时间的全局直方图统计进一步去噪，并求取精细化的测距结果，以满足星载光子计数激光雷达实时处理需求。提出的去噪算法流程示意图如图2所示。

图2 去噪算法流程示意图

Fig. 2 Schematic diagram of denoising algorithm flow

具体的算法流程如表1所示。

表1 粗精去噪算法

Table 1 Coarse and fine denoising algorithm

粗精去噪算法

输入：采样到的待去噪的光子事件时间数据的集合 $\{t^{i}\}$

输出：去噪后的光子事件时间数据集合 $\{t_{s f}\}$ 和探测目标的距离值 $Z$

粗去噪

1：确定单次去噪对应的发射脉冲数量 $N$ 和 $N$ 发脉冲接收到的光子事件数量 $N_{c}$

2： $\{t^{(i)}\} = S o r t (t^{i}), i = \{1,2, \dots, N_{c}\}$ ， $\{t^{(i)}\}$ 为 $\{t^{i}\}$ 排序后的数据

3：在 $\{t^{(i)}\}$ 中每次连续选取 $n$ 个数据 $\{t_{n}^{(i)}\}$ 进行计算

$t_{m a x} = m a x \{t_{n}^{(i)}\}$

$t_{m i n} = m i n \{t_{n}^{(i)}\}$

$t_{∆} = t_{m a x} - t_{m i n}$

$\{t_{s c}\}$ =｛ $t_{n}^{(i)} |d^{i} = \frac{t_{∆}}{n - 1} < T_{p}, T_{p} = 6 σ_{s}|$ ， $\{t_{s c}\}$ 为可能的信号光子事件时间数据集合， $d^{i}$ 为光子事件时间数据密度， $T_{p}$ 为发射激光脉冲脉宽

4：采用重叠选取的方式选择下一次的n个数据，直至遍历完 $\{t^{(i)}\}$ 中的所有数据，连续两次选取数据重叠数量为n-1。为保证查全率，判定为噪声光子事件的数据可在其他样本中改判为信号光子事件，反之，判定为信号光子事件的数据不可改判为噪声光子事件。

精去噪

5： $t_{m a x}^{'} = h i s t o g r a m \{t_{s c}\}$ ， $t_{m a x}^{'}$ 为直方图统计峰值对应的时间值

6： $\{t_{s f}\} = \{t_{s c} | (t_{m a x}^{'} - T_{p}) \leq t_{s c} \leq (t_{m a x}^{'} + T_{p})\}$ ， $\{t_{s f}\}$ 为最终去噪后保留下来的信号光子事件时间数据

7： $t_{m e a n} = m e a n \{t_{s f}\}$ ， $t_{m e a n}$ 为筛选出的信号光子时间数据的平均值。

8： $Z = \frac{c}{2} t_{m e a n}$ ， $c$ 为光速， $Z$ 为距离值。

3 阈值参数确定

上述实时去噪算法核心参数为累计光子事件个数N和每次取出的样本个数n。N和n的选取需要综合考虑发射激光脉宽、系统噪声率和探测概率等因素，同时确保在满足系统探测指标的同时，系统硬件资源开销最小，以达到最优化设计。

3.1　 $N$ 的确定

N的取值需满足探测概率和虚警概率要求。为保证探测的有效性，在区间 $(T_{s} - 3 σ_{s}, T_{s} + 3 σ_{s})$ 内的探测概率和虚警概率均需处于合理的取值范围。本算法选取探测概率不小于90%、虚警概率不大于10%以确保对信号的有效识别，即

\begin{array}{l} P_{s} (6 σ_{s}) = 1 - e x p (- N_{s} \times N - N_{n} \times N) \geq 0.9, \\ N \geq 1 且 N \in Z \\ P_{n} (6 σ_{s}) = 1 - e x p (- N_{n} \times N) \leq 0.1, \\ N \geq 1 且 N \in Z \end{array}

，

(10)

可以得到：

\begin{array}{l} N \times N_{n} = N \times 6 σ_{s} f_{n} \leq 0.11, N \geq 1 且 N \in Z \\ N \times (N_{s} + N_{n}) \geq 2.3, N \geq 1 且 N \in Z \end{array}

(11)

星载光子计数体制激光雷达的噪声包括太阳背景光噪声和探测器件的暗计数噪声。由于暗计数噪声比太阳背景光噪声低2~3个数量级，一般可以忽略。太阳背景光噪声可由下式计算得到^［

24］：

\begin{array}{l} f_{A} = \int_{0}^{\infty} N_{λ}^{0} Δ λ \frac{β θ_{r}^{2} L^{2} π T^{s e c θ_{s}} A_{r} η_{r} η_{q} T_{h o r} (l)}{4 h v \cdot 4 π L^{2}} d l \\ = N_{λ}^{0} Δ λ \frac{β θ_{r}^{2} T^{s e c θ_{s}} A_{r} η_{r} η_{q}}{16 h v} \int_{0}^{\infty} T_{h o r} (l) d l \end{array}

(12)

图3为以ATLAS系统参数作为输入，在不同太阳天顶角条件下，随地表反射率变化的系统噪声率曲线。

图3 ATLAS噪声率与太阳天顶角和地表反射率理论关系图^［

24］

Fig. 3 Theoretical relationship between ATLAS noise rate and solar zenith angle and surface reflectance^［

24］

根据公式（12）预估得到系统噪声率，通过激光雷达方程计算得到平均信号光子数，进而可以确定去噪算法中的参数N。

3.2　 $n$ 的确定

定义随机变量 $X_{1}, X_{2}, \dots, X_{n}$ 次序统计量 $X_{(1)}, X_{(2)}, \dots, X_{(n)}$ 的极差为 $R_{n} = X_{(n)} - X_{(1)}$ 。

提出的粗去噪方法其实就是在排好序的 $N$ 个数据里面，每次求 $n$ 个次序统计量的极差，判定极差是否在设定的范围内，即排序后的 $n$ 个光子事件时间数据的最大值和最小值的差值，是否满足一定阈值要求。

为了确定 $n$ 值，不妨假设 $n$ 个光子全部为噪声光子，噪声光子事件时间数据满足均匀分布 $U (a, b)$ ，其中 $a$ 为所有噪声光子事件时间数据最小值， $b$ 为所有噪声光子事件时间数据最大值，将数据归一化到 $(0,1)$ 的范围内后，均匀分布 $U (0,1)$ 的极差为 $β$ 分布 $β (n - 1,2)$ ^［

26］。当

n

取不同的值时，

β

分布曲线如图4所示。

图4 β分布曲线：（a） $β$ 分布的概率密度函数曲线；（b） $β$ 分布的累积分布函数曲线

Fig. 4 β distribution curve：（a） the probability density function curve of β distribution；（b） the cumulative distribution function curve of β distribution

图4（a）为 $β$ 分布概率密度分布函数（Probability Density Function， PDF）曲线，横坐标为极差值，纵坐标为概率密度。从曲线可以看出，随着 $n$ 不断增大，PDF整体分布不断往右移动，峰值点右移。在同样的极差值条件下（图中黄色竖直线所示），PDF越小，即噪声越强，噪声分布范围小于设定阈值的概率（即噪声被误判为信号的概率）也越小，但是同时需要处理的数据量也不断增加。

图4（b）为 $β$ 分布累积概率分布函数（Cumulative Distribution Function， CDF）曲线，横坐标为极差值，纵坐标为累积概率值，对应着虚警概率。以10%的虚警概率作为阈值（图（5）中黄色水平线所示），随着 $n$ 不断增大，极差值也越来越大，即在同样的虚警概率要求下，噪声光子数越多，允许的噪声光子事件的分布范围也越大，星上系统需要存储和处理的数据也会越大。为了实时存储和处理的方便， $n$ 取3比较合适。此时均匀分布 $U (0,1)$ 的极差阈值为0.1958，即归一化的3个噪声光子事件时间数据的最大值和最小值的差值不能超过0.1958，对应均匀分布 $U (a, b)$ 的极差阈值为 $0.195 8 \times (b - a)$ ，其中 $b - a$ 为所有噪声光子事件时间数据的最大分布范围，即实际的3个噪声光子事件时间数据的最大值和最小值的差值不能超过噪声时间数据分布范围的0.1958倍。一般而言，实际系统噪声的时间分布范围很宽， $n$ 取3可以满足虚警概率要求。

假设 $n$ 个光子全部为信号光子，由于其满足高斯分布，信号光子事件时间数据的极差理论上为无穷大，但是如果只考虑有效的 $(- 3 σ_{s}, 3 σ_{s})$ 范围，则极差最大值为 $6 σ_{s}$ ，探测概率为0.9974，满足探测概率大于0.9的要求。即当 $n$ 个信号光子事件时间数据的最大值和最小值的差值为激光发射脉冲的脉宽时，都可以满足探测概率的要求。综合噪声和信号光子事件的时间分布特性， $n$ 取3可以满足去噪要求。

4 蒙特卡罗仿真

蒙特卡罗方法是一种以抽样和随机数的产生为基础的随机性方法，其基本原理是通过数字模拟实验，得到所要求解的出现某种事件的概率^［

27］。本文采用了蒙特卡罗仿真验证上述算法可行性及阈值选择的有效性。具体仿真流程如图5所示。

图5 蒙特卡罗仿真过程流程图

Fig. 5 Process flow diagram of Monte Carlo Simulation

（1）根据有效信号回波模型，对回波信号时域分布以极小时间片τ进行离散化，并生成探测概率分布曲线和探测概率密度曲线，根据探测概率密度曲线利用随机数生成器生成信号光子事件数据序列 $t_{1}$ ；

（2）根据噪声模型，设置噪声率参数与选通门宽度，即约束噪声光子数据量与数据范围，利用均匀分布随机数生成器，构造均匀分布的噪声光子事件数据序列 $t_{2}$ ；

（3）将高斯分布和均匀分布的光子事件时间数据序列 $t_{1}$ 和 $t_{2}$ 进行随机混合，构造实际光子事件时间数据序列 $t_{3}$ ；

（4）利用去噪算法对光子事件时间数据序列 $t_{3}$ 中的噪声进行剔除，得到去噪后的光子事件时间数据序列 $t_{4}$ ；

（5）求取去噪后的光子事件时间数据序列 $t_{4}$ 的均值 $t_{m e a n}$ ，换算成目标的距离值Z $= \frac{1}{2} c t_{m e a n}$ ，其中 $c$ 为光速，通常取 $3 \times 10^{8} m / s$ 。

在蒙特卡罗仿真过程中，参考ATLAS系统的参数设计，激光重频设为 $10 k H z$ ，相邻两次激光脉冲发射时间间隔为0.1 ms，高斯分布的信号均方根脉宽 $σ_{s}$ 设为 $0.67 n s$ ，选取平均信号光子数 $N_{s}$ 为3个。ATLAS每次对200次激光脉冲发射即0.02 s内的回波探测信号进行直方图统计，本文构造10次激光发射脉冲，信号光子总数为30个，时间分布范围 $6 σ_{s}$ 为4 ns，以5 000 ns为真值基准，信号均值 $T_{s}$ 设为5 000 ns。针对不同的地表类型，实际的仿真过程和结果如下。

4.1　陆地

取陆地反射率为0.2~0.3，噪声率一般小于3 MHz。将噪声率设为3 MHz，即1 s有 $3 \times 10^{6}$ 个噪声光子，1 ns有 $3 \times 10^{- 3}$ 个噪声光子。10次激光发射重复探测时间内，噪声光子数为 $3 \times 10^{2}$ 个，均匀分布在 $10^{4}$ ns的时间范围内，蒙特卡罗仿真结果如图6所示。

图6 陆地仿真结果：（a）原始信号和噪声光子事件时间数据分布图；（b）去噪后光子事件时间数据分布曲线图，方框中是图表的图例说明，横线表示的是去噪后光子事件时间数据的均值，绿色圆圈表示的是信号，红色圆圈表示的是噪声

Fig. 6 Land simulation results：（a） distribution map of the original signal and noise photon event time data；（b） distribution curve of the photon event time data after denoising， the diagram is illustrated in the box， the horizontal line represents the mean of the photon event time data after denoising， the green circle represents the signal， and the red circle represents the noise

图6（a）横坐标为采样序号，按照到达时刻顺序采样；纵坐标为光子事件时间数据值，绿色为信号光子，红色为噪声光子。可以看到，信号光子事件时间数据聚集在高斯信号均值 $T_{s}$ 附近，噪声事件时间数据均匀分布。图6（b）为去噪后光子事件时间数据分布曲线图，横坐标为去噪后的光子事件序号，纵坐标为去噪后的光子事件飞行时间数据，绿色为信号光子，红色为噪声光子。上述说明在以下图示中保持一致。由于去噪过程中按照光子事件时间数据大小进行了排序，所以去噪后光子事件时间数据分布曲线为一条递增的曲线，从图6（b）中可以看出，经过去噪后，信号光子被提取了出来，且信号光子分布在 $T_{s}$ 附近，误差100 ps量级，实现了点云去噪。

4.2　海洋

海洋噪声率一般小于5 MHz，设为5 MHz，即1 s有 $5 \times 10^{6}$ 个噪声光子，1 ns有 $5 \times 10^{- 3}$ 个噪声光子。10次激光发射重复探测时间内，噪声光子数为 $5 \times 10^{2}$ 个，均匀分布在 $10^{4}$ ns的时间范围内，蒙特卡罗仿真结果如图7所示。

图7 海洋仿真结果：（a）原始信号和噪声光子事件时间数据分布图；（b）去噪后光子事件时间数据分布曲线图，方框中是图表的图例说明，横线表示的是去噪后光子事件时间数据的均值，绿色圆圈表示的是信号，红色圆圈表示的是噪声

Fig. 7 Ocean simulation results：（a） distribution map of the original signal and noise photon event time data；（b） distribution curve of the photon event time data after denoising， the diagram is illustrated in the box， the horizontal line represents the mean of the photon event time data after denoising， the green circle represents the signal， and the red circle represents the noise

从图7（a）中可以看出，相对图6（a），信号光子分布不变，但是噪声光子更加密集。从图7（b）中可以看出，信号光子被提取了出来，且信号光子分布在 $T_{s}$ 附近，误差为100 ps量级，达到了去噪的效果。

4.3　陆地冰

陆地冰噪声率按8 MHz计算，则1 s有 $8 \times 10^{6}$ 个噪声光子，1 ns有 $8 \times 10^{- 3}$ 个噪声光子。10次激光发射重复探测时间内，噪声光子数为 $8 \times 10^{2}$ 个，均匀分布在 $10^{4} n s$ 的时间范围内，蒙特卡罗仿真结果如图8所示。

图8 陆地冰仿真结果：（a）原始信号和噪声数据分布图；（b）去噪后数据分布曲线图，方框中是图表的图例说明，横线表示的是去噪后光子事件时间数据的均值，绿色圆圈表示的是信号，红色圆圈表示的是噪声

Fig. 8 Land ice simulation results：（a） distribution map of the original signal and noise photon event time data；（b） distribution curve of the photon event time data after denoising， the diagram is illustrated in the box， the horizontal line represents the mean of the photon event time data after denoising， the green circle represents the signal， and the red circle represents the noise

相对图7（a），图8（a）中信号分布不变，噪声分布更加密集，数据量进一步增多。图8（b）显示去噪后有少量噪声光子没有被剔除，但是有效信号都被提取了出来，达到去噪效果。

4.4　海洋冰

海冰反射率为0.8~0.9，噪声率设为10 MHz，即1 s有 $10^{7}$ 个噪声光子，1 ns有 $10^{- 2}$ 个噪声光子。10次激光发射重复探测时间内，噪声光子数为 $10^{3}$ 个，均匀分布在 $10^{4}$ ns的时间范围内，蒙特卡罗仿真结果如图9所示。

图9 海洋冰仿真结果：（a）原始信号和噪声数据分布图；（b）去噪后数据分布曲线图，方框中是图表的图例说明，横线表示的是去噪后光子事件时间数据的均值，绿色圆圈表示的是信号，红色圆圈表示的是噪声

Fig. 9 Ocean ice simulation results：（a） distribution map of the original signal and noise photon event time data；（b） distribution curve of the photon event time data after denoising， the diagram is illustrated in the box， the horizontal line represents the mean of the photon event time data after denoising， the green circle represents the signal， and the red circle represents the noise

图9结果显示，算法能有效实现去噪功能。

为对去噪算法的效果进行定量评价，引入如下评价指标：

查准率：

P = \frac{T P}{T P + F P}

(13)

查全率：

R = \frac{T P}{T P + F N}

(14)

调和平均值：

F = \frac{2 P R}{P + R}

，

(15)

其中 $T P$ ：信号光子被判定为信号的个数； $F P$ ：噪声光子被判定为信号的个数； $T N$ ：噪声光子被判定为噪声的个数； $F N$ ：信号光子被判定为噪声的个数。

在选定的系统参数和阈值条件下，对四种目标进行多次重复仿真实验，实验结果如表2所示。

表2 粗精去噪算法的实验结果

Table 2 Experimental results of coarse and fine denoising algorithm

序号	目标	TP	FP	TN	查全率R	查准率P	调和平均值F	数据压缩比
1	陆地	30	0.2	299.8	1	0.993 377 5	0.996 677 7	10.927 152 32
2	海洋	30	0.9	499.1	1	0.970 873 8	0.985 221 7	17.152 103 56
3	陆地冰	30	0.4	799.6	1	0.986 842 1	0.993 377 5	27.302 631 58
4	海洋冰	30	0.6	999.4	1	0.980 392 2	0.990 099	33.660 130 72

表3 常见的直方图统计去噪算法实验结果

Table 3 Experimental results of common histogram statistical denoising algorithms

序号	目标	TP	FP	TN	查全率R	查准率P	调和平均值F	数据压缩比
1	陆地	30	0.2	299.8	1	0.993 377 5	0.996 677 7	10.927 152 32
2	海洋	30	0.9	499.1	1	0.970 873 8	0.985 221 7	17.152 103 56
3	陆地冰	30	0.5	799.5	1	0.983 606 6	0.991 735 6	27.213 114 75
4	海洋冰	30	0.6	999.4	1	0.980 392 2	0.990 099	33.660 130 72

从实验结果来看，提出的粗精去噪算法与常见的直方图统计去噪算法的结果基本一致，这主要是因为提出的算法精去噪部分也是用的直方图统计方法。从各项指标来看，提出的去噪算法能够保证查全率，即能够找出所有的信号光子；查准率在0.97~0.99之间，即还有一部分噪声光子被误当作信号光子提取出来，但是整体而言，调和平均值F在0.98~0.99之间，属于正常水平。同时可以看到，去噪算法将数据量压缩到了10倍量级，大大减轻了后续数据处理和存储的压力。

常见的直方图统计去噪算法需要预先设定一段统计的时间窗口范围，对于不同反射率的目标采集探测时间是一样的，所以统计的目标信息会出现不足和冗余的情况。提出的粗精去噪算法根据目标反射特性灵活变动时间窗口范围，依据系统噪声率和目标反射率确定系统统计需要的最少光子时间数据个数进行处理，将有限的资源按需分配到每种目标，在减少处理的数据量和处理时间的同时，保证每种目标的探测精度。ATLAS每次直方图统计需要200个激光脉冲，实际上对于不同的目标，需要的数据量并不一样，反射率高的目标只需要几十个甚至几个激光发射脉冲。本文提出的粗精去噪算法，通过粗去噪对目标的差异性进行适配，使得统计的数据量最小，同时又保留了直方图统计所具有的精度。

5 实测数据测试

本文选取ICESat-2在海洋、海冰区域的数据对提出的算法进行验证，选取了海洋区域的一轨数据和海冰区域的一轨数据，同时使用同一轨数据的强波束数据和弱波束数据进行验证，对应不同平均信号光子数情况。ICESat-2所搭载的激光雷达ATLAS激光器重频为10 kHz，发射激光均方根脉宽为0.67 ns。根据实测数据的噪声率和平均信号光子数，选择合适的N值，验证本文所提出的实时去噪算法。去噪效果如图10所示。

图10 ICESat-2数据去噪结果，（a）、（b）、（c）、（d）为强波束数据；（e）、（f）、（g）、（h）为弱波束数据；（a）、（b）、（e）、（f）为海冰数据；（c）、（d）、（g）、（h）为海洋数据；（a）、（c）、（e）、（g）噪声率为1 MHz；（b）、（d）、（f）、（h）噪声率为5 MHz，蓝色为噪声光子，绿色为信号光子，红色为去噪后的信号光子

Fig. 10 ICESat-2 data denoising results，（a），（b），（c） and （d） are strong beam data；（e），（f），（g） and （h） are weak beam data；（a），（b），（e） and （f） are sea ice data；（c），（d），（g） and （h） are ocean data； the noise rate of （a），（c），（e） and （g） is 1 MHz； the noise rate of （b），（d），（f） and （h） is 5 MHz， blue represents noise photons， green represents signal photons， and red represents denoised signal photons

表4 ICESat-2数据去噪结果

Table 4 ICESat-2 data denoising results

目标	噪声率（MHz）	平均信号光子数（个）	N	查全率（%）	查准率（%）	调和平均值（%）
海洋	5	1.24	5	97.65	93.56	95.56
海洋	5	8.81	5	99.16	99.36	99.26
海洋	1	1.24	50	99.98	99.11	99.54
海洋	1	8.81	50	98.86	99.88	99.36
海冰	1	1.1	50	96.42	97.73	97.07
海冰	1	4.42	50	94.38	99.55	96.9
海冰	5	1.1	5	94.84	94.44	94.64
海冰	5	4.42	5	95.82	97.92	96.86

从ICESat-2实测数据的实验结果来看，本文提出的实时去噪算法的实际查全率在0.94~0.99之间，查准率在0.93~0.99之间，整体调和平均值在0.94~0.99之间。即该实时去噪算法可以实现星载光子计数体制激光雷达的去噪功能。

基于栅格的点云数据去噪方法是从人卫激光测距数据处理领域发展而来的，是一种基于概率探测的数据处理方法^［

25］。由于其计算方法较为简单，栅格去噪法常用于光子点云的快速去噪。选取了海冰和海洋区域，从查全率、查准率、调和平均值、计算时间等四个方面将本文提出的去噪算法与栅格法比对，结果如图11所示。

图11 ICESat-2数据栅格法去噪结果，（a）、（b）、（c）、（d）为强波束数据；（e）、（f）、（g）、（h）为弱波束数据；（a）、（b）、（e）、（f）为海冰数据；（c）、（d）、（g）、（h）为海洋数据；（a）、（c）、（e）、（g）噪声率为1 MHz；（b）、（d）、（f）、（h）噪声率为5 MHz，蓝色为噪声光子，绿色为信号光子，红色为去噪后的信号光子

Fig. 11 The denoising results of ICESat-2 data grid method，（a），（b），（c） and （d） are strong beam data；（e），（f），（g） and （h） are weak beam data；（a），（b），（e） and （f） are sea ice data；（c），（d），（g） and （h） are ocean data； the noise rate of （a），（c），（e） and （g） is 1 MHz； the noise rate of （b），（d），（f） and （h） is 5 MHz， blue represents noise photons， green represents signal photons， and red represents denoised signal photons

表5 ICESat-2数据栅格法去噪结果

Table 5 Denoising results of ICESat-2 data grid method

目标	噪声率（MHz）	平均信号光子数（个）	查全率（%）	查准率（%）	调和平均值（%）
海洋	5	1.24	98.72	79.56	88.11
海洋	5	8.81	92.90	90.05	91.45
海洋	1	1.24	98.72	94.31	96.47
海洋	1	8.81	92.90	92.19	92.55
海冰	1	1.1	99.62	96.00	97.78
海冰	1	4.42	99.21	98.26	98.73
海冰	5	1.1	99.42	84.97	91.63
海冰	5	4.42	98.97	94.60	96.74

表6 去噪计算时间对比

Table 6 Comparison of denoising calculation time

目标	噪声率（MHz）	平均信号光子数（个）	本文算法计算时间（s）	栅格法计算时间（s）
海洋	5	1.24	0.642 2	0.716 0
海洋	5	8.81	0.754 1	0.871 8
海洋	1	1.24	0.320 0	0.346 7
海洋	1	8.81	0.438 3	0.473 3
海冰	1	1.1	0.325 8	0.354 7
海冰	1	4.42	0.342 0	0.378 0
海冰	5	1.1	0.634 1	0.738 6
海冰	5	4.42	0.643 2	0.798 9

结果对比可看到，相比于快速去噪的栅格点云去噪法，本文提出的实时去噪方法的查准率明显优于栅格法，查全率与栅格法相当。随着信噪比的降低，本文算法和栅格去噪法的调和平均值均有所下降，但本文算法的效果明显优于栅格法。从去噪时间对比，本文提出的算法运行时间均小于栅格法运行时间，计算时间降低了10%~20%，且随着信噪比的降低，运行效率对比差异更明显。

6 结论

为了将星载单光子激光雷达海量原始数据快速提炼为高价值小数据量的信息，大幅降低下传数据量，并将激光探测有效信息及时发送至用户，降低用户信息获取延时，并辅助减轻地面处理压力，本文提出一种粗精结合的星载光子计数激光测距实时去噪算法，利用光子事件数据的局部密度信息，初步快速识别出部分噪声光子并剔除，再对保留的光子事件进行全局直方图统计，求取精细化的测距结果。相比于传统的直方图统计和其他的时域去噪算法，本方法能够充分联合并利用局部密度和全局分布信息对噪声和信号进行快速实时去噪，满足卫星实时处理需求。

References

ZHENG Xiang-Yang， DING Yu-Xing， WANG Hai-Wei， et al. Processing of point cloud data for photon counting laser range finder system under movement condition［J］. INFRARED郑向阳，丁宇星，王海伟，等.运动条件下光子计数激光测距系统点云数据的处理［J］.红外， 2016， 37（6）：24. [百度学术]

XIE Feng，YANG Gui，SHU Rong，et al. An adaptive directional filter for photon counting LiDAR point cloud data［J］.Journal of Infrared and Millimeter Waves谢锋，杨贵，舒嵘，等.方向自适应的光子计数激光雷达滤波方法［J］.红外与毫米波学报，2017，36（1）：107-113. [百度学术]

LI Song， WANG Yue， HUANG Ke， et al. An improved clustering based denoising method for single photon point cloud data： CN201811424017.4 ［P］. CN109344812A.李松，王玥，黄科，等.一种改进的基于聚类的单光子点云数据去噪方法： CN201811424017.4［P］. CN109344812A. [百度学术]

LI Kai，ZHANG Yong-Sheng，TONG Xiao-Chong，et al. Research on De-noising and Filtering Algorithm of Single Photon LiDAR Data.［J］.Navigation and Control李凯，张永生，童晓冲，等.单光子激光雷达数据去噪与滤波算法［J］. 导航与控制， 2020， 19（1）： 67-76. [百度学术]

JIAO Huihui， XIE Junfeng， LIU Ren， et al. Discussion on Denoising Method of Photon Counting LiDAR for Satellite Ground Observation［J］. Spacecraft Recovery＆Remote Sensing焦慧慧，谢俊峰，刘仁，等. 星载对地观测光子计数激光雷达去噪方法浅析［J］. 航天返回与遥感，2021，42（5）： 140-150. [百度学术]

HE Guang-Hui，WANG Hong，FANG Qiang，et al. Spaceborne photon counting lidar point cloud denoising method with the adaptive mountain slope［J］. Journal of Infrared and Millimeter Waves何光辉，王虹，方强，等.山地坡度自适应星载光子计数激光雷达点云去噪方法［J］.红外与毫米波学报， 2023， 42（2）： 250-259. [百度学术]

Wang Zhen-Hua， Chen Shi-Xian， Kong Wei， et al. Comparison and Analysis of Denoising for Photon-Counting LiDAR Data［J］. Laser & Optoelectronics Progress王振华，陈诗贤，孔伟，等.光子计数激光雷达中光子点云滤波方法的比较与分析［J］.激光与光电子学进展， 2023， 60（6）： 334-341. [百度学术]

MCGARRY J L. ATLAS Flight Science Receiver Algorithms［J/OL］. 2019［2022-10-25］. https：//ntrs. nasa. gov/citations/20190031952. [百度学术]

ICE， CLOUD， and Land Elevation Satellite （ICESat-2） Project Algorithm Theoretical Basis Document （ATBD） For ATLAS Level 1A Processing （John DiMarzio，David Hancock，2017） https：//icesat-2.gsfc.nasa.gov/sites/default/files/files/ATL01_13Nov2017.pdf [百度学术]

Martino A. J.， M. R. Bock， C. Gosmeyer， et al. Ice， Cloud， and Land Elevation Satellite （ICESat-2） Project Algorithm Theoretical Basis Document （ATBD） for ATL02 （Level-1B） Data Product Processing， Version 6（2022）. ICESat-2 Project， DOI： 10.5067/7LI3JNHLHB6X. https：//icesat-2.gsfc.nasa.gov/sites/default/files/page_files/ICESat2_ATL02_ATBD_r006.pdf [百度学术]

Neumann T. A.， A. Brenner， D. Hancock， et al. Ice， Cloud， and Land Elevation Satellite （ICESat-2） Project Algorithm Theoretical Basis Document （ATBD） for Global Geolocated Photons ATL03 （2022）， Version 6. ICESat-2 Project， DOI： 10.5067/GA5KCLJT7LOT. https：//icesat-2.gsfc.nasa.gov/sites/default/files/page_files/ICESat2_ATL03_ATBD_r006.pdf [百度学术]

Kwok R.， Markus T.， Morison J.， et al. Profiling Sea Ice with a Multiple Altimeter Beam Experimental Lidar （MABEL）［J］. Journal of Atmospheric and Oceanic Technology， 2014， 31（5）： 1151-1168. [百度学术]

Jasinski M. F.， Stoll J. D.， Cook W. B.， et al. Inland and Near-Shore Water Profiles Derived from the High-Altitude Multiple Altimeter Beam Experimental Lidar （MABEL）［J］. Journal of Coastal Research， 2016， 76： 44-55. [百度学术]

Aabdalati W.， Zwally H. J.， Bindschadler R.， et al. The ICESat-2 Laser Altimetry Mission［J］. Proc. IEEE， 2010， 98（5）： 735-751. [百度学术]

Dabney P.， Harding D.， Abshire J.， et al. The Slope Imaging Multi-Polarization Photon-Counting Lidar： Development and Performance Results［C］. 2010 IEEE International Geoscience and Remote Sensing Symposium. Honolulu， HI， USA： IEEE， 2010： 653-656. [百度学术]

RAPP J， GOYAL V K. A Few Photons Among Many： Unmixing Signal and Noise for Photon-Efficient Active Imaging［J］. IEEE Transactions on Computational Imaging， 2017， 3（3）： 445-459. [百度学术]

LUO Le， WU Chang-Qiang， LIN Jie， et al. Time-domain denoising based on photon-counting LiDAR［J］. Optics and Precision Engineering骆乐，吴长强，林杰，等.基于光子计数激光雷达的时域去噪［J］.光学精密工程， 2018， 26（5）： 1175-1180. [百度学术]

Zhu Wen-Hua， Wang Shu-Chao， Wang Kai-Di， et al. Adaptive acquisition time scanning method for photon counting imaging system［J］. Acta Physica Sinica朱文华，汪书潮，王凯迪，等.光子计数成像系统的自适应累积时间扫描方法［J］.物理学报， 2022， 71（15）： 314-322. [百度学术]

ZHU Lei， HUANG Geng-Hua， OUYANG Jun-Hua， et al. Study on time interval measurement system in photon counting imaging LIDAR［J］.Journal of Infrared and Millimeter Waves朱磊，黄庚华，欧阳俊华，等.光子计数成像激光雷达时间间隔测量系统研究［J］.红外与毫米波学报， 2008， 27（6）： 461-464. [百度学术]

HOU Li-Bing， GUO Ying， HUANG Geng-Hua， et al. A Time-to-digital Converter Used in Photon-counting LIDAR［J］. Journal of Infrared and Millimeter Waves侯利冰，郭颖，黄庚华，等.光子计数激光雷达时间-数字转换系统［J］.红外与毫米波学报， 2012， 31（3）： 243-247. [百度学术]

HOU Li-Bing， HUANG Geng-Hua， KUANG Yao-Wu， et al.Research of Photon Counting Laser Ranging Technology［J］. Science Technology and Engineering侯利冰，黄庚华，况耀武，等.光子计数激光测距技术研究［J］.科学技术与工程， 2013， 13（18）： 5186-5190. [百度学术]

Gardner C. S.. Target Signatures for Laser Altimeters： An Analysis［J］. Applied Optics， 1982， 21（3）： 448-453. [百度学术]

HUANG Ke， LI Song， MA Yue， et al. Detection Probability Model of Single-Photon Laser Altimetry and Its Range Accuracy［J］. Chinese Journal of Lasers黄科，李松，马跃，等.单光子模式激光测高探测概率模型与精度分析［J］.中国激光， 2016， 43（11）： 229-234. [百度学术]

LI Song， ZHANG Zhi-Yu， MA Yue， et al. Detection Theory， Link Simulation， and Data Processing of Spaceborne Photon Counting Lidar ［M］. Science Press， 2022.李松，张智宇，马跃，等.星载光子计数激光雷达探测理论、链路仿真与数据处理［M］. 科学出版社， 2022. [百度学术]

Degnan J J. Photon-counting multikilohertz microlaser altimeters for airborne and spaceborne topographic measurements［J］. 2002， 34（3）：503-549. [百度学术]

XIONG Jia-Bing， WANG Zhi-Xiang. Properties of Uniform Distribution Order Statistics ［J］. [百度学术]

Studies in College Mathematics熊加兵，王志祥.均匀分布次序统计量的性质［J］.高等数学研究， 2010， 13（1）： 17-19. 10.3969/j.issn.1008-1399.2010.01.006 [百度学术]

Coyac A， Hespel L， Riviere N， et al. Performance assessment of simulated 3D laser images using Geiger-mode avalanche photo-diode： tests on simple synthetic scenarios： SPIE Security + Defence， 2015［C］. [百度学术]

您是第12692036 位访问者

主管单位：中国科学院

主办单位：中国科学院上海技术物理研究所，中国光学学会

地址：上海市玉田路500号电话：021-25051553

51La

首页

学报简介

编委会

征稿简则

版权声明

开放获取

出版道德声明

相关下载

联系我们

常见问题

English Version

星载光子计数激光测距雷达的实时去噪方法 PDF

摘要

关键词

引言

1 光子计数体制激光雷达探测概率模型

2 去噪算法流程