基于原始探测光子泊松分布特性的直方图去噪算法研究

摘要

ICESat-2作为新一代多波束激光测高卫星采用了光子计数体制，其探测数据中存在大量噪声，造成光子数据在轨处理和传输时面临巨大挑战。为了在轨高效地对原始探测数据进行去噪处理，以降低星地传输数据量，本文设计了一种基于原始探测光子泊松分布特性的直方图去噪算法。

算法核心

该算法分为光子点云的垂直直方图化和倾斜直方图化。首先，根据光子传输距离将点云数据段划分为二维格网形成垂直直方图，利用直方图箱光子数的均值和标准差计算信噪分离阈值，并对信号光子赋予低、中和高置信度标签来表征信号可靠性；其次，对识别的中、高置信度光子进行线性拟合获取坡度信息，将光子传输距离投影到沿坡面垂直的方向形成倾斜直方图，进行二次信号光子识别且合并置信度标签；最后，对噪声光子进行剔除以实现原始数据的压缩下传。

研究成果

实验结果表明，垂直直方图算法对信号光子召回率R均为1，查准率P和调和F值平均在0.90以上，相较于垂直+倾斜直方图、DBSCAN和OPTICS的运行效率分别达到12倍、3 473倍和1 528倍以上，平均运行时间仅0.048 s，可高效实现数据去噪。垂直直方图算法去噪结果初步满足在轨处理效率要求(<0.25 s)，本文有望为未来星载光子计数激光雷达数据的在轨去噪提供技术参考。

关键词

星载激光雷达 ICESat-2 ATL02 光子计数点云去噪分级直方图空间聚类在轨处理

研究背景

星载激光雷达技术

星载激光雷达作为主动遥感探测技术，可获取高精度的高度信息。目前，国内外的星载激光雷达卫星有冰、云和陆地高程卫星(Ice, Cloud, and Land Elevation Satellite, ICESat)、全球生态系统动力学调查(Global Ecosystem Dynamics Investigation, GEDI)、资源3号卫星(ZY-03)和高分7号卫星(GF-7)等。

激光雷达技术已广泛应用于地形测量、植被冠层参数反演、水深及水下地形测量以及全球碳储量变化监测等方面。

研究挑战

ICESat-2上搭载的先进地形光子计数激光测高仪(Advanced Topographic Laser Altimeter System, ATLAS)是全球首颗采用单光子计数体制的星载激光雷达系统。与传统线性体制激光雷达相比，星载光子计数激光雷达获取的光子点云数据背景噪声率高，对噪声进行高效剔除是光子数据能够科学应用的关键。

随着波束增多，同一探测时刻光子接收数量将出现倍增，且星上处理器在数据处理和传输方面资源有限，如何进一步控制数据量并确保地面信息不丢失将是后续多波束光子计数激光雷达发展的关键问题。

研究问题

噪声问题

光子计数激光雷达数据中存在大量噪声光子，影响数据质量和后续应用

处理效率

星上处理器性能有限，需要高效算法满足在轨实时处理需求

数据传输

大量光子数据对星地传输带宽造成压力，需要有效压缩方法

地形适应性

算法需要适应不同地表类型(城市、植被、海洋等)的光子分布特征

研究方法

研究思路

如图1(a)所示，ATLAS原始探测信息主要包含光子发射时间和传输时间，即ATL02中/delta_time和/ph_tof字段。以数据段中第一个光子为起点，获取每一个光子的相对发射时间和相对传输时间，并结合光子传输速度c(约299 792 458 m/s)和卫星飞行速度v(约7.0 km/s)可以获得未经误差改正(硬件延迟、大气散射、折射和潮汐位移等)的相对传输距离和相对沿轨距离。

针对NASA采集的多波束激光雷达数据(Multiple Altimeter Beam Experimental Lidar, MABEL)，相关学者对其研究分析证实了随机返回的光子独立且遵循泊松分布。在此理论基础上，NASA设计了基于泊松理论的云检测大气直方图算法、距离窗压缩和ATL03置信度分级算法。

图1 基于ATL02光子传输时间的距离转换前后光子点云分布比较

泊松去噪算法

垂直直方图化

ICESat-2上搭载的ATLAS探测器探测频率为50 Hz，即以0.02 s(200次脉冲)间隔对光子回波数据进行记录，除此之外，硬件所能探测的最小脉冲宽度约为20 ns，转换为距离约为3 m。直方图划分时，水平方向以t_min=0.02 s为间隔和垂直方向以z_min=3 m为间隔初始化直方图箱大小，并对光子点云数据直方图化。

P(X=k)=λ^k/k! e^-λ (k=0,1,2,⋯)

式中：X是被ATLAS时间标记的光子事件；P是光子事件X发生k次的概率；!是阶乘；

自然对数函数的底数；参数λ为单位时间内光子事件X平均发生的次数，即在0.02 s内3 m高的直方图箱中的光子数。

图2 光子点云垂直直方图化

背景箱与信号箱识别

根据泊松分布的特性，将光子数低于阈值的点视为噪声点。计算背景噪声阈值α以确定包含背景噪声光子的直方图箱(背景箱)：

α(m)=u(m)+e^aσ(m)

式中：e^a=2.5是经验阈值，其作用是提高找寻背景箱的概率。当H(i, m)<α(m)时，将直方图箱H(i, m)标记为背景箱H_bg(i, m)。

确定信号阈值φ和信噪比Snr，以此找寻直方图中的信号箱：

φ=u_bg+e_mσ_bg

Snr=H(i,m)/u_bg

图3 基于直方图算法的信号箱识别结果

置信度标签分配

计算信号箱的信噪比S_NR，并根据下式为满足S_NR的箱中的每个光子分配低、中、高置信度参数来表征光子作为信号的可信度：

ph_conf(m)= ⎧
2 S_NR(m)<20
3 20≤S_NR(m)<50
4 50≤S_NR(m) ⎫

式中：ph_conf是光子置信度标签，数值越大，光子作为信号的可信度越高。

为了防止单个直方图中信号光子识别不全而导致信号光子丢失，根据信号光子相对离散程度R确定信号探测概率阈值r_limit，以此作为查寻率r的约束：

R(m)=(M_ph(m)-u(m))/∑H(i,m)

r(m)=∑(ph_conf(m)>0)/∑H(i,m)

当R≥10^-3时，结合ATLAS系统对目标探测概率确定r_limit取值为0.2，反之取值为0.1。

倾斜直方图化

为了保证在陡峭地形条件下信号光子不遗漏，需相邻直方图间沿轨方向时间重叠率在50%以上。拟合坡度，对原始光子点云倾斜处理，将光子传输距离投影到沿坡面垂直的方向形成倾斜直方图，进行二次信号光子识别且合并置信度标签。

y=c₀+c₁t

φ=arctan(c₁)

d'(m)=d(m)×cosφ(m)-t×sinφ(m)

图4 倾斜直方图化过程中原始光子传输距离倾斜示意图

实验与分析

实验数据

ICESat-2/ATLAS飞行轨道高度约500 km，ATLAS系统发射的532 nm激光脉冲对水体具有穿透性，其最大穿透深度达40 m。本文主要利用ATL02数据进行实验，选取了不同地理特征的实验区，数据可在NASA网站获取。

图5 试验区地理位置分布与ICESat-2星下点轨迹图

表1 各实验区数据详细信息统计

地表类型	数据编号	地理位置
植被区	Data1	大兴安岭
高山区	Data2	祁连山脉
沙漠区	Data3	塔克拉马干沙漠
海洋区	Data4	印度洋
极地区	Data5	格陵兰岛
陆海区	Data6	南海
湖泊区	Data7	哈德森湾
城市区	Data8	芝加哥

性能评估指标

召回率(R)

TP/(TP+FN)

被正确识别的信号光子占所有真实信号光子的比例

查准率(P)

TP/(TP+FP)

被识别为信号的光子中真正是信号光子的比例

调和平均(F)

2PR/(P+R)

召回率和查准率的调和平均值，综合评估算法性能

实验结果

使用本文直方图去噪算法分别对8个试验区数据进行去噪和置信度分级，同时将ATL03中信号光子映射回ATL02数据，对比分析信号光子分布差异。

图10和图11分别为城市区与陆海区去噪结果。城市中建筑物的密集分布和多层次结构会导致光子在不同高度和位置上的分布不均匀。532 nm的激光脉冲具有水体穿透能力，在近海岸水下地形测绘方面得到了广泛应用。

图6 城市区结果对比

图7 陆海区结果对比

算法性能比较

表3 基于不同算法的实验区ATL02光子数据去噪结果评估指标统计比较

实验区	垂直直方图			垂直+倾斜直方图			DBSCAN			OPTICS
实验区	R	P	F	R	P	F	R	P	F	R	P	F
植被区	1.000	0.915	0.956	1.000	0.916	0.956	0.998	1.000	0.999	1.000	0.994	0.997
高山区	1.000	0.832	0.908	1.000	0.850	0.919	1.000	0.786	0.880	1.000	0.753	0.859

运行效率分析

为了评估本文直方图算法在星上进行数据去噪压缩的可行性，分析了不同方法的去噪效果以及算法复杂度。卫星在空间中的运行速度约为7 km/s，通过ATL02数据提供的飞行时间参数(/tof字段)可粗略计算沿轨飞行距离，截取7 km（约0.064°）的数据段利用4种算法进行去噪实验。

垂直直方图相较于垂直+倾斜直方图、DBSCAN和OPTICS运行效率分别达到12倍、3 473倍和1 528倍以上，平均运行时间仅0.048 s，初步认为满足星上处理效率要求(<0.25 s)。

图8 基于不同算法的试验区ATL02光子数据去噪运行效率统计比较

结论与展望

研究结论

通过对8种地表类型数据实验结果表明，直方图去噪算法具有较高的地形适应性和较好的去噪压缩效果。垂直直方图和垂直+倾斜直方图信号召回率均达到了1，F值平均高于0.90。

本文直方图算法设计考虑了硬件参数与探测器性能，降低了由于参数经验阈值取值不合理对去噪效果的影响。

基于密度的空间聚类算法精度虽明显较高，但其较低的运行效率无法满足星上处理器性能要求。垂直直方图平均运行时间仅0.048 s，初步满足星上处理效率要求(<0.25 s)。

未来展望

考虑到星上处理器与日常计算机性能、程序语言以及运行环境差异等因素影响，本文评估结果可能存在偏差，有必要进一步展开研究。

由于星载光子计数激光雷达可接收500 km高度范围内返回的任何光子，考虑到ATL02数据为星上压缩处理后的数据产品，因此，后续需要通过模拟最原始的探测数据进行距离窗压缩和本文的直方图去噪实验。

在本文目前研究中，主要关注实验区内的地理空间特征差异，但尚未充分考虑四季更替所带来的地理特征变化对直方图去噪算法的影响。

研究价值

本文设计的基于原始探测光子泊松分布特性的直方图去噪算法，通过利用光子数的泊松分布特性，优先识别并剔除噪声光子，可以有效提高信号光子的提取效率，并降低误差，这种方法能够更好地适应不同数据分布情况下的去噪需求。研究结果验证了直方图算法的适用性和有效性，
有望为光子计数星载激光雷达光子点云在轨去噪提供理论参考。