一种无输入参数的强噪声背景下ICESat-2点云去噪方法

刘翔，张立华，戴泽源，陈秋，周寅飞

（1 海军大连舰艇学院 军事海洋与测绘系，辽宁 大连 116018）（2 海军大连舰艇学院海洋测绘工程军队重点实验室，辽宁 大连 116018）（3 91001部队，北京 100071）

0 引言

2018年9月15日发射的ICESat-2（The Ice，Cloud and Land Elevation Satellite-2）是美国第二代冰、云和陆地高程卫星，用以执行冰盖高程测量和变化监测、陆地高程测量、全球植被高度测量以及监测云和气溶胶等任务［1-3］。相较于ICESat所搭载的全波形地学激光测高系统（Geoscience Laser Altimeter System，GLAS），ICESat-2搭载的先进地形激光高度计系统（Advanced Topographic Laser Altimeter System，ATLAS）为微脉冲光子计数测高雷达系统，具有高灵敏性、高重频的特性，可以获得光斑更小、密度更高的光子点云数据，增强了对全球三维信息的获取能力［2］。同时这也是全球范围内第一次将该项技术运用于星载平台之上［1］，代表了未来星载激光测高重要发展趋势。但由于ICESat-2光子计数激光雷达发射和探测的激光脉冲都为弱信号，难以区分目标物体表面返回脉冲、大气散射、太阳辐射及仪器自身噪声［4-5］，导致其回波数据中掺杂了大量的背景噪声。在强噪声背景下，较低的信噪比使得信号光子与噪声光子难以区分，如何从背景噪声中提取信号光子已成为ICESat-2数据应用的一项关键而基础的工作。

针对单光子激光测高雷达数据的去噪，国内外学者进行了诸多研究。文献［5-6］通过统计光子与其邻近光子局部距离进行去噪处理，但邻近光子数量为经验设定值；文献［7-11］基于光子滤波核内光子密度值的差异进行去噪，但滤波核范围参数需人为进行设定；文献［12］以DBSCAN聚类方法为基础，将所用圆形滤波核形状改为水平椭圆进行信号光子聚类，但该聚类方法存在对输入参数敏感的问题；文献［13］则以OPTICS方法进行聚类，降低了对输入参数的敏感度，但所采用椭圆滤波核半径仍需通过试错法来得到。以上去噪方法存在参数人为设定较为繁杂或对参数敏感的问题。尽管文献［14-15］将光子分布光栅化为图像，以边缘提取的图像处理方式进行去噪，避免了参数设定过程，但光栅化必然会造成光子点云信息不可逆的损失，影响数据的后续应用。

为了解决单光子激光测高雷达数据去噪对输入参数的依赖性问题，ZHANG Guoping等［16］提出了一种基于四叉树的无参数输入去噪方法，该方法无需人为经验设定的光子在高程和沿轨距离方向上的范围以及周围光子数量，只需利用四叉树将每个光子合理分隔即可通过分隔层值表征光子密度，从而可做到在去噪时无需要参数输入，在ICESat-2的模拟数据MATLAS去噪实验中取得了良好的效果。但对于强噪声背景下的实测ICESat-2数据，四叉树去噪方法处理局部稀疏、相距较近的噪声光子时表征密度不合理，会将较多噪声光子误识别为信号光子，从而影响去噪效果。

为此，本文针对四叉树去噪方法的缺陷，首先以剪枝四叉树的层值表征光子密度，避免强噪声背景下局部稀疏、相距较近的噪声光子被四叉树分隔层次较多而被误表征为密度过大，自适应求取信号光子的密度阈值，完成一级去噪；然后通过箱线图进一步去除少量未能被剪枝四叉树识别局部较密离群噪声光子，完成二级去噪。

1 数据与研究方法

1.1 研究区域概况

本文研究区域及实验所用ICESat-2数据分布情况如图1所示，研究区域1位于美国的北达科他州，该区域为平原地带，地势较为平坦，植被覆盖稀疏；研究区域2位于美国的加利福尼亚州，该区域靠近内华达山脉，地势起伏较大，中等植被覆盖。两个研究区域地势起伏差异较大、地表覆盖类型不同，具有较强的代表性。

图1 研究区域位置及数据分布示意图Fig.1 Location of study area and distribution of data

1.2 数据来源

1.2.1 ICESat-2/ATLAS数据

2019年5 月美国国家冰雪数据中心NSIDC（National Snow and Ice Data Center）开始向全球用户公开发布ICESat-2/ATLAS数据［4］。本文使用其中的Level-2级数据产品ATL03作为本次去噪的实验数据，ATL03为全球定位光子数据，它记录了光子精确的经纬度和高程信息，是生成诸如海冰高程（ATL07）［17］、陆地植被高程（ATL08）［18］等其他产品的数据基础［19］。

1.2.2 机载激光雷达数据

为验证本文方法去噪效果，实验选择美国国家生态观测站网络NEON（National Ecological Observatory Network）发布的机载激光雷达高程数据产品作为验证数据，包括数字表面模型（Digital Surface Model，DSM）和数字地面模型（Digital Terrain Model，DTM），二者分辨率皆为1 m。同时为保证验证的准确性，机载激光雷达数据观测时间与ICESat-2数据观测时间均为同年同月，以减少因时相变化所造成的不同数据产品之间的高程误差。

1.3 研究方法

1.3.1 一级去噪

在ICESat-2光子点云中，信号光子与噪声光子空间分布差异往往较大，具体表现为信号光子较噪声光子分布更加密集［20-21］，这种密度的差异是当前ICESat-2点云去噪的基本依据。因此，选取合适的密度阈值是正确分类信号光子和噪声光子的关键。本文根据信号光子和噪声光子密度不同的特点，挖掘Otsu法在自适应求取灰度特性有明显差别的两类目标分割阈值上的优势，通过剪枝四叉树的层值来定量表示光子密度，以密度替代灰度，求取信号光子和噪声光子密度最大类间方差作为分类阈值，对两类光子进行分类。

1.3.1.1 四叉树层值表征光子密度

四叉树层值表征光子密度是利用四叉树方法对光子分隔，即在高程和沿轨距离二维空间D上按照光子的位置将空间进行四象限递归划分，依照光子最终所处四叉树象限叶节点的层值来定量光子的的密度。分隔过程为首先根据ICESat-2光子数据的范围给定初始根节点空间D0，记录D0中光子个数为sum(D0)，当sum(D0)超过设定的空间光子最大容纳值maxlimit，则将当前空间划分为右上（Ⅰ）、左上（Ⅱ）、左下（Ⅲ）、右下（Ⅳ）四个相同的子象限空间，计算各子空间内光子数，对超过maxlimit的子空间再次进行划分，递归执行以上划分，直至各子空间内光子数在最大容纳值之内，此时光子处于各叶节点所对应的划分空间中。空间划分准则可归纳为

式中，D′k-1为第k-1层中可进行划分的象限集合；Dk-1为第k-1层象限集合；X为Dk-1中某个象限空间；sum(X)为象限空间X内光子数；为保证去噪精度，通常需对每个光子进行分隔，因此maxlimit一般设为1。

图2为某区域利用四叉树分隔阴影部分光子的示意图，图3为相应的四叉树树状结构图，从中可以看出由于噪声光子分布较信号光子稀疏，因此在分隔时只需要较少的划分次数即可被分隔开，处于的叶节点象限空间层值k较小，而信号光子因分布密集所处叶节点象限空间层值k较高，层值k与密度存在着明显的正相关关系，因此采用层值来对信号和噪声光子密度进行表征，将层值作为光子的密度标签。

图3 四叉树分隔光子树状结构图Fig.3 Tree structure of photons separated by quadtree

1.3.1.2 剪枝四叉树层值表征光子密度

由于噪声光子的分布呈现随机性，在利用四叉树分隔ICESat-2光子时，会出现局部稀疏情形下某些噪声光子相距较近而需要较多的划分次数方可分隔的现象，如图4红框标识的象限所示。此时利用四叉树方法表征光子密度时，这类噪声光子由于所在叶节点象限层值较高，接近信号光子的层值，如图5所示，导致在去噪过程中被误分为信号光子。特别是在强噪声背景下，信噪比较低，光子点云中会大量散布此类噪声光子，进而导致四叉树方法难以达到较好的去噪效果。

图5 四叉树误分隔光子树状结构图Fig.5 Tree structure of photons incorrectly separated by quadtree

进一步分析光子密度的表征过程，光子层值需要反映的是光子与周围光子而不是与最邻近光子之间的疏密程度，在图4红框标识的象限中，内部光子虽然距离较近，但是象限进行一次四等划分后并不能将其分隔，从其与周围光子分布的稀疏程度来看其密度仍然较小，所以此时象限空间所在层对应的层值就代表了其内部光子的密度。基于此，为减少噪声光子的误识别，本文引入数据搜索过程中的“剪枝”思想改进原有四叉树生成过程，并采用剪枝四叉树层值来表征光子密度。数据搜索算法中的“剪枝”是通过一定的策略对树结构进行优化，删除部分无效节点，提升算法效率。而本文所述的“剪枝”过程，则是先对将要进行划分的某象限进行分析判断，若划分一次后没有光子被分隔开，则此划分不进行，同时此象限的划分相应终止，以此来剪掉由于不合理的划分而生成的枝状结构。划分准则为

图4 四叉树误分隔光子示意图Fig.4 Photons incorrectly separated by quadtree

式中，nsum[split(X)]≥1表示为对象限空间X四等划分后，子象限空间内光子数不小于1的空间个数，即存在光子的子象限空间个数。

对四叉树划分改进后，划分准则不再依据maxlimit，而是顾及了象限空间光子的分布情况。按照划分准则，对象限空间X的划分存在以下两种情况：

1）sum(X)≤1，即象限空间X不存在或只存在一个光子，对此象限不进行划分，此象限划分终止；

2）sum(X)≥2，即象限空间X存在多个光子，先对象限空间X进行分析判断，若进行一次四等划分后光子在同一象限，则不进行此划分，同时此象限划分终止；若不在同一象限，则进行划分，对划分后的子象限空间按照其内光子数迭代执行以上两种划分情况，直至划分终止。

图6为剪枝四叉树划分后阴影部分光子的分隔情况，相对于图4来说，红框象限按照划分准则不再进行划分。从所对应的树状结构图图7可以看出，通过上述划分方式，红色阴影部分的枝状结构被剪掉，噪声光子所处层级降低，层值变小，避免了其被误识别为信号光子。

图6 剪枝四叉树分隔光子示意图Fig.6 Photons separated by pruned quadtree

图7 剪枝四叉树分隔光子树状结构图Fig.7 Tree structure of photons separated by pruned quadtree

1.3.1.3 Otsu法求取光子密度阈值

由于受地表反射率、太阳高度角、大气散射等因素的影响，沿轨方向上不同区域内光子的噪声密度会有所差异［11］。为适应这种沿轨方向上信号光子和噪声光子密度对比的不同，在沿轨方向划分100 m等距离窗口，分窗口利用Otsu法自适应求取各窗口光子的密度阈值，依据密度阈值进行噪声光子和信号光子的分类。

设一个窗口内所有光子层值取值范围为[0，k]，假定密度阈值为t，将窗口内的光子分为层值为[0，t-1]和[t，k]两类，类间方差为

式中，

式中，pi为窗口内层值为i的光子所占比，ω1、ω2分别为层值小于t和大于等于t的光子所占比，μ1、μ2分别为层值小于t和大于等于t的光子层值均值，μ为窗口内所有光子的层值均值。

当部分噪声光子被错分为信号光子或者部分信号光子被错分为噪声光子时，就会导致类间方差变小，因此类间方差最大时意味着两类光子错分的概率最小［22-23］，所以窗口内光子的最佳密度阈值t*为

即t从(0，k)中取值，t*为使类间方差σ2(t)最大时所对应的t值。此时窗口内层值小于t*的光子分类为噪声光子，其他为信号光子，依据分类结果将噪声光子去除完成本文的一级去噪。

1.3.2 二级去噪

经过剪枝四叉树去噪以后，绝大部分噪点得到了消除。虽然噪声光子在整体上比信号光子分布稀疏，但依然存在孤立噪声光子簇，此时利用剪枝四叉树层值表征密度时，此类局部噪声光子由于分布密集，因此其层值偏大，被误识别为信号光子。对于这种少量的离群噪声光子，采用箱线图的方法进行去除。由于箱线图是通过光子高程四分位数来设置上下阈值进行离群点的判定，因此为顾及高程变化对箱线图去噪的影响，同样采用分窗口的形式进行去噪，窗口宽度在地形复杂区域和平坦区域分别取100 m和50 m。其中箱线图的上下限为［24-25］

式中，Q1、Q3分别为窗口内所有光子高程的下四分位数和上四分位数。如图8为箱线图去噪的示意图，蓝点为高程值大于Upperlimit或小于Lowerlimit的光子，其被识别为离群噪点。

图8 箱线图离群点检测Fig.8 Outlier detection by box-plot

1.3.3 去噪效果验证

为有效地验证去噪效果，采取定性和定量相结合的方式与机载Lidar生成的DTM和DSM进行对比。定性验证采用将四叉树方法、本文所提方法去噪后光子拟合的地表和冠层顶部曲线，与相应位置处DTM和DSM剖面的高程曲线进行比对的方式。在植被区地表光子和冠顶光子分别位于信号光子点云底部和顶部，由此采用移动窗口的方式，将窗口内高程值最小和最大的信号光子分别作为地表和冠层的种子光子，然后利用三次样条插值函数拟合得到地表和冠层曲线［9］，对比的真实地表曲线为DTM剖面高程曲线，真实冠顶曲线为DSM剖面高程曲线。定量验证为计算种子光子高程和验证高程数据之间均方根误差RMSE和决定系数R2，其中RMSE和R2表示为

式中，N为计算的光子点数，hi为ICESat-2光子高程值，ĥi为光子所处位置对应的DTM或DSM高程值。

2 试验结果与分析

本文在北达科他州研究区域1和加利福尼亚州研究区域2选择两条ICESat-2条带开展去噪实验。图9和图10分别为两区域光子在沿轨方向和高程的二维分布图，中间密集部分为信号光子，大量的噪声光子较为稀疏地分布在信号光子带的两侧，光子数据的信噪比皆比较低，处于强噪声背景下。

图9 研究区域1原始点云Fig.9 Initial point cloud in study area 1

图10 研究区域2原始点云Fig.10 Initial point cloud in study area 2

2.1 研究区域1去噪结果

在研究区域1开展去噪实验，图11和图12分别为利用四叉树方法和剪枝四叉树方法的去噪结果，其中四叉树方法和剪枝四叉树方法采用了相同的窗口划分策略和光子密度阈值求取方法，去噪结果中蓝色光子为去除的噪声光子，红色光子为去噪后得到的信号光子。从中，可以明显看出，在强噪声背景下，剪枝四叉树方法相比四叉树方法误识别的噪声光子数量得到减少，去噪效果更好。将上述两种方法经过二次去噪得到如图13和14所示结果图，可以看到四叉树方法去噪结果在经过二次去噪后仍存在一定量的噪声光子，而本文所提剪枝四叉树方法和二级去噪则基本去除了所有噪声光子。这是由于二次去噪所采用的箱线图法的去噪性能会随着离群噪点的增多而降低，而剪枝四叉树方法减少了误识别噪声光子的数量，所以使剪枝四叉树结合二级去噪在去噪效果上要优于四叉树方法结合二级去噪。

图11 四叉树方法在研究区域1去噪结果Fig.11 Denoising result of quadtree method in study area1

图12 剪枝四叉树方法在研究区域1去噪结果Fig.12 Denoising result of pruned quadtree method in study area1

图13 四叉树方法和二级去噪在研究区域1去噪结果Fig.13 Denoising result of quadtree method and second-level denoising in study area 1

图14 剪枝四叉树方法和二级去噪在研究区域1去噪结果Fig.14 Denoising result of pruned quadtree method and second-level denoising in study area 1

为进一步验证本文所提方法相较于四叉树方法具有更优的去噪性能以及去噪结果的准确程度，本文以机载激光雷达生成的高程数据产品为标准进行对比试验。由于此区域植被覆盖稀疏，因此在研究区域1只使用DTM作为对比数据。分别用四叉树方法提取的信号光子和本文所提方法提取的信号光子得到地表种子光子，将拟合的地表曲线与DTM剖面高程曲线进行对比，结果如图15所示。可以看到，利用四叉树方法提取的信号光子所拟合的地表曲线与DTM剖面的真实地表曲线之间存在着较大偏差，无法正确反映地表走势，而本文所提方法提取的信号光子拟合的地表曲线则与DTM剖面的地表曲线基本一致。进一步地，定量计算两种方法所得地表种子光子高程与DTM高程之间的均方根误差与决定系数，其结果如表1所示。可以看出，本文方法的RMSE为0.91 m，优于四叉树方法所对应的10.51 m，R2为0.997，非常接近理想值1，计算光子的高程与分辨率为1 m的机载激光雷达高程数据产品高程基本一致，较四叉树方法对应的0.713，有显著提高。

图15 研究区域1地表曲线结果Fig.15 Ground curves in study area 1

表1 研究区域1地表精度评价Table 1 Ground accuracy evaluation in study area 1

2.2 研究区域2去噪

图16和图17分别为利用四叉树方法和剪枝四叉树方法在研究区域2的去噪结果，可以看出，四叉树方法误识别的噪声光子较多，而剪枝四叉树方法明显减少了误识别噪声光子的数量。图18为四叉树方法和二级去噪的结果，仍然残留了一定量的噪声光子，图19所示的剪枝四叉树方法和二级去噪则基本去除了所有的噪声光子。

图16 四叉树方法在研究区域2去噪结果Fig.16 Denoising result of quadtree method in study area 2

图17 剪枝四叉树方法在研究区域2去噪结果Fig.17 Denoising result of pruned quadtree method in study area 2

图18 四叉树方法和二级去噪在研究区域2去噪结果Fig.18 Denoising result of quadtree method and second-level denoising in study area 2

图19 剪枝四叉树方法和二级去噪在研究区域2去噪结果Fig.19 Denoising result of pruned quadtree method and second-level denoising in study area 2

同样为进一步验证本文所提方法相较于四叉树方法具有更好的去噪性能以及去噪结果的准确程度，对于研究区域2，从四叉树方法提取的信号光子和本文所提方法提取的信号光子得到冠顶和地表种子点，将拟合的冠顶和地表曲线分别与DSM和DTM剖面高程曲线进行对比，如图20和21所示。从图中可以看出，利用四叉树方法提取的信号光子拟合的冠顶和地表曲线与DSM和DTM剖面的冠顶和地表曲线都存在较大的偏差，而利用本文所提方法提取的信号光子拟合的曲线则与DSM和DTM剖面曲线基本一致。从表2和表3可以看出，研究区域2本文方法所得信号光子中提取的冠层和地表种子光子高程与DSM和DTM高程之间的RMSE分别为3.56 m和2.47 m，明显优于四叉树方法对应的90.17 m和90.92 m，R2为0.998和0.999，接近理想值1，计算光子的高程与分辨率为1 m的机载激光雷达高程数据产品高程基本一致，与四叉树对应的0.4和0.156相比，有显著提高。

图20 研究区域2冠顶曲线结果Fig.20 Canopy top curves in study area 2

图21 研究区域2地表曲线结果Fig.21 Ground curves in study area 2

表2 研究区域2冠顶精度评价Table 2 Canopy top accuracy evaluation in study area 2

表3 研究区域2地表精度评价Table 3 Ground accuracy evaluation in study area 2

3 结论

本文提出了一种无输入参数强噪声背景下ICESat-2点云去噪方法，以剪枝四叉树层值表征光子密度，自适应求取信号光子密度阈值，进行一级去噪；然后通过箱线图进一步去除少量离群噪声光子，进行二级去噪。本文所提方法避免了强噪声背景下局部稀疏、相距较近的噪声光子被四叉树分隔层次较多而被误表征为密度过大，自适应求取信号光子的密度阈值，并对相聚密集噪声光子进行了有效去除，与四叉树方法相比，具有更优的去噪性能。本文所提方法去噪后的信号光子所拟合的冠顶和地表曲线与分辨率为1 m的机载激光雷达高程数据产品高程剖面曲线基本一致。

当然，限于所掌握的公开数据有限等原因，本文只选取了较为典型的区域进行了比对实验，更多区域的实验验证还有待以后进一步进行。同时，本文去噪方法中采用了划分窗口的方式，可能会存在不同窗口去噪结果相邻边界不一致的现象，后续研究中将结合窗口拓展等方式来降低此类边缘效应。