CN113281716A - 一种光子计数激光雷达数据去噪方法 - Google Patents

Abstract

本发明提供了一种光子计数激光雷达数据去噪方法，属于激光雷达数据处理技术领域。本发明通过对原始光子数据进行四叉空间分割，将原始光子数据构成的空间转化为四叉树结构，然后计算原始光子数据中每个光子的隔离深度，并将隔离深度小于设定的深度阈值的光子标记为噪声光子并去除，从而从原始光子数据中获得信号光子，本发明方法能够从不同信噪比的原始数据中准确提取信号光子，使准确、高效去除光子计数激光雷达数据中的噪声成为可能。

Description

一种光子计数激光雷达数据去噪方法

技术领域

本发明涉及一种光子计数激光雷达数据去噪方法，属于激光雷达数据处理技术领域。

背景技术

遥感是全天时、全天候的对地观测手段，星载激光测高仪已经被证明在获取高精度森林垂直结构和地面高程方面具有独特优势。ICESat-2(Ice Cloud and LandElevation Satellite-2)是NASA于2018年9月发射的新一代激光测高仪，也是全球首颗星载光子计数激光测高仪，其探测性能达到了一个全新的高度，其搭载的ATLAS(theAdvanced Topographic Laser Altimeter System)装载了6束532nm的蓝绿色激光，以10kHz的重复频率对地进行测量。在ICESat-2正式发射之前，NASA进行了包括MABEL和MATLAS在内的多次机载实验，针对具有代表性的地形和测绘条件获取了大量的模拟数据，这些数据的信噪比远低于波形数据，大量噪声广泛分布于数据中，为了生成高质量的数据产品，有必要提出一种准确、有效的噪声去除算法。

目前已经有一些基于光子密度的去噪算法，这些算法均假设信号光子的局部空间分布比噪声光子更加密集。例如，Zhang and Kerekes提出了一种使用水平椭圆形搜索区域的DBSCAN聚类算法，由于光子空间分布与地形相关，该算法能够在平原取得不错的效果，但是在山区的去噪结果不佳；Xie Feng等通过在水平椭圆搜索区域中加入方向参数提出了一种改进DBSCAN算法，增强了地形的鲁棒性，改进DBSCAN算法被最广泛应用于星载光子计数激光雷达数据去噪，然而，搜索区域的尺寸和比例对DBSCAN聚类算法的去噪结果影响比较大，准确计算出这些参数是比较困难的；Bowei Chen等提出了一种使用椭圆搜索区域改进的LOF(local outlier factor)算法，并且也发现了搜索区域的形状对去噪结果的影响很大，该算法虽然证明了水平椭圆搜索区域的效果最优，但是并未给出输入参数的计算方法；Xiaoxiao Zhu提出了一种OPTICS(ordering points to identify the clusteringstructure)的聚类算法并使用水平椭圆搜索区域进行改进，该算法降低了对聚类参数的敏感性，但是搜索区域的半径仍需要人为给定。

综上所述，目前基于光子密度的去噪算法的去噪效果受输入参数的影响较大，不具有输入参数鲁棒性，若输入参数设置的准确就能取得较好的去噪效果，若输入参数设置的不准确去噪效果就不佳，由于在实际应用中输入参数是很难准确计算出来的，导致目前星载光子计数激光雷达数据的去噪效果无法保证。

发明内容

本发明的目的在于提供一种光子计数激光雷达数据去噪方法，能够有效去除光子计数激光雷达数据中的噪声光子。

为了实现上述目的，本发明提供了一种光子计数激光雷达数据去噪方法，该方法包括以下步骤：

(1)获取待处理的光子计数激光雷达数据作为原始光子数据；

(2)利用四叉树方法将原始光子数据所覆盖的空间转化为四叉树：先将原始光子数据所覆盖的空间进行四叉空间分割得到4个相同大小的子空间，根据原始光子数据的坐标将其分配在对应的子空间；然后对每个子空间中包含的光子个数进行判断，若子空间中包含的光子个数＞1，就对该子空间继续进行四叉空间分割，直至分割后的子空间中包含的光子数≤1，若子空间中包含的光子个数≤1，则该子空间停止分割；将原始光子数据所覆盖的空间作为根节点，将每次分割得到的子空间作为子空间节点，将停止分割时子空间中的光子坐标作为叶子结点，从而将原始光子数据所覆盖的空间转化为四叉树；

(3)计算原始光子数据中每个光子在所述四叉树中的隔离深度，光子在四叉树中的隔离深度等于从根节点到该光子所在叶子结点所经历的边数；

(4)将隔离深度小于深度阈值的光子标记为噪声光子，通过去除原始光子数据中的噪声光子得到原始光子数据中的信号光子。

本发明是一种不依赖于输入参数的去噪方法，通过对原始光子数据进行四叉空间分割，将原始光子数据所覆盖的空间转化为四叉树结构，然后计算原始光子数据中每个光子的隔离深度，并将隔离深度小于设定的深度阈值的光子标记为噪声光子并去除，从而从原始光子数据中获得信号光子，本发明方法能够从不同信噪比的原始数据中准确提取信号光子，使准确、高效去除光子计数激光雷达数据中的噪声成为可能。

进一步地，为了最大限度的区分信号和噪声，所述深度阈值利用大律法确定。

进一步地，采用大律法确定深度阈值的过程为：将原始光子数据划分为若干设定长度的窗口，将窗口中的光子按照隔离深度降序排列，然后按照下式确定深度阈值：

式中，n是窗口中的光子总数，t是当前假设为信号光子的光子数，ω₀(t)为当前假设的信号比例，ω₁(t)假设的噪声比例，μ₀(t)为信号平均深度，μ₁(t)为噪声平均深度，将类间方差σ²最大时的隔离深度ID(t)作为深度阈值。

附图说明

图1是本发明方法实施例中光子计数激光雷达数据去噪方法流程图；

图2-1是本发明方法实施例中的四叉树原理示意图；

图2-2是本发明方法实施例中的四叉树结构示意图；

图3-1是本发明方法实施例中基于改进DBSCAN方法的山区日间数据噪声去除结果图；

图3-2是本发明方法实施例中基于改进DBSCAN方法的山区夜间数据噪声去除结果图；

图3-3是本发明方法实施例中基于本实施例方法的山区日间数据噪声去除结果图；

图3-4是本发明方法实施例中基于本实施例方法的山区夜间数据噪声去除结果图；

图4-1是本发明方法实施例中基于改进DBSCAN方法、日间、提取地面点所计算出的地表高程和参考DTM高程之间的散点图；

图4-2是本发明方法实施例中基于改进DBSCAN方法、夜间、提取林冠点高程与对应位置DTM高程相减得到的植被高度和CHM高度的散点图；

图4-3是本发明方法实施例中基于本实施例方法、日间、提取地面点所计算出的地表高程和参考DTM高程之间的散点图；

图4-4是本发明方法实施例中基于本实施例方法、夜间、提取林冠点高程与对应位置DTM高程相减得到的植被高度和CHM高度的散点图；

图4-5是本发明方法实施例中基于改进DBSCAN方法、日间、提取地面点所计算出的地表高程和参考DTM高程之间的散点图；

图4-6是本发明方法实施例中基于改进DBSCAN方法、夜间、提取林冠点高程与对应位置DTM高程相减得到的植被高度和CHM高度的散点图；

图4-7是本发明方法实施例中基于本实施例方法、日间、提取地面点所计算出的地表高程和参考DTM高程之间的散点图；

图4-8是本发明方法实施例中基于本实施例方法、夜间、提取林冠点高程与对应位置DTM高程相减得到的植被高度和CHM高度的散点图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。

方法实施例：

本实施例的光子计数激光雷达数据去噪方法(以下简称本实施例方法)，如图1所示，本实施例方法是一种不依赖于输入参数的去噪方法，不同于基于光子密度的去噪算法，本实施例方法对原始光子计数激光雷达数据(以下简称原始光子数据)进行四叉空间分割，将光子空间位置转化为树结构，计算原始光子数据中每个光子的树深度并将深度较小的光子标记为噪声光子并去除，从而从原始光子数据中获得信号光子。

同时，利用MATLAS数据和G-LiHT数据测试并评价本实施例方法的性能，此外，还将本实施例方法与改进DBSCAN方法进行了对比，实验结果表明，本实施例方法能够从不同信噪比的原始数据中准确提取信号光子，而且表现得比改进DBSCAN方法更好。总得来说，本实施例方法使准确、高效去除光子计数激光雷达数据中的噪声成为可能。

下面对本实施例中用到的性能测试数据进行简单介绍：

(1)MATLAS数据

在ICESat-2发射前，NASA针对代表性地形，在不同时间获取了一些机载模拟数据用于早期科学分析。MABEL(Multiple Altimeter Beam Experimental LiDAR)是NASA针对ICESat-2任务获取的模拟数据之一，而MATLAS则根据ATLAS的仪器参数对MABEL数据中的信号光子和噪声光子数量进行了调整。MATLAS数据与ICESat-2数据类似，但是在测量条件上进行了精心设计，其数据信噪比更有代表性。

因此，本实施例决定使用MATLAS数据而非ICESat-2数据来评估本实施例方法在不同信噪比下的表现，评估时，8条西海岸的MATLAS数据被使用，包括4条日间数据和4条夜间数据。

(2)G-LiHT数据

G-LiHT系统是NASA为了解陆地生态系统的组成、结构和作用而设计的机载测量系统，可以同时获取地表三维点云、超光谱影像和热感影像。G-LiHT的点云数据使用Riegl的VQ-480激光雷达获取，点云密度大于10个/㎡，可以用于生产高精度的地形和植被高程产品。由于G-LiHT的部分数据轨迹与MATLAS数据存在重合，由点云生成的(DTM)和(CHM)被用作参考数据。其中，DTM为数字地形模型(Digital Terrain Model)，CHM为林冠高度模型(Canopy Height Model)。

下面对本实施例方法进行详细介绍：

(1)将原始光子数据所覆盖的空间转化为四叉树

与基于光子密度的去噪算法假设信号光子的空间分布比噪声光子更加密集不同，本实施例方法是一种基于孤立的去噪算法，这种算法假设噪声光子疏散的空间分布使其比信号光子更容易被孤立，基于这个假设，本实施例方法首先将原始光子数据所覆盖的空间进行四叉空间分割得到4个相同大小的子空间，并按照他们的位置命名为LL(lower left)、LR(lower right)、UL(upper left)、UR(upper right)，然后，从原始光子数据中随机选取光子并根据坐标将其分到对应的子空间中；然后，对每个子空间中包含的光子个数进行判断，若子空间中包含的光子个数＞1，就对该子空间继续进行四叉空间分割，直至分割后的子空间中包含的光子数≤1，若子空间中包含的光子个数≤1，则该子空间停止分割。将原始光子数据所覆盖的空间作为根节点，将每次分割得到的子空间作为子空间节点，将停止分割时子空间中的光子坐标作为叶子结点，从而将原始光子数据所覆盖的空间转化为四叉树，此时，每个光子都被隔离在叶子结点中。

定义隔离深度是指光子从根节点遍历四叉树，直到遍历结束于该光子所在叶子结点所经历的边数。

由于四叉树是光子空间位置的树结构表示，每个非叶子结点代表了一组光子的坐标范围，而叶子结点表示了一个光子的空间坐标，所以每个光子的隔离深度与该光子被孤立所需要的树分裂次数是一致的。

图2-1和图2-2分别展示了一组光子数据集中不同光子被孤立的情况及其在四叉树中的位置。图中，Noise(Noise Photon)代表噪声光子，Signal(Signal Photon)代表信号光子，ROOT为根节点，Node为节点。对于光子数据，由于噪声光子距离其他光子较远，四叉树只需要更少次分裂就可以将其孤立；而信号光子的情况则相反。在树结构中，信号光子需要更多次树分裂，位于树的底部；噪声光子需要更少次树分裂，位于树的上部。尽管只是一个例子，但是仍然不难得出结论：信号光子的隔离深度大于噪声光子。

(2)利用大律法(即类间方差最大法)确定深度阈值进行噪声去除

在构造的四叉树中，噪声光子容易被分隔，隔离深度较小。相比之下，信号光子的隔离深度往往较大。因此，隔离深度可以作为发现噪声光子的一种依据。

受外界环境影响，同一条光子数据的信噪比会发生变化，低信噪比区域中光子的隔离深度显然比高信噪比中光子的隔离深度大。因此，为了最大限度地区分信号光子和噪声光子，将数据按照沿轨道距离(沿轨道距离指卫星沿着轨道运行的距离)划分为若干100m(即设定长度，可根据实际需要设置)的窗口，将窗口中光子按照隔离深度降序排列。使用大津法计算隔离深度阈值，将光子自适应地划分为信号光子和噪声光子。

大津法计算隔离深度阈值的公式如下所示。

式中，n是窗口中的光子总数，t是当前假设为信号光子的光子数，ω₀(t)为当前假设的信号比例，ω₁(t)假设的噪声比例，μ₀(t)为信号平均深度，μ₁(t)为噪声平均深度。

将类间方差σ²最大时的隔离深度ID(t)作为深度阈值。对于任一光子，若其隔离深度大于ID(t)，将其标记为信号光子；否则，标记为噪声光子。

(3)性能评价

由于MATLAS数据中没有给出光子的信号标签，所以需要从侧面评价算法的准确性。首先使用移动窗口法提取地面点和林冠点并基于此计算地表高程和植被高度。然后，计算地表高程和植被高度和参考数据之间的差异。Coefficient of determination(R²)和root-mean-square error(RMSE)被计算，被用作反映提取结果和实测数据之间误差的大小。

除此之外，为了更加客观得评价本实施例方法，改进DBSCAN方法被作为对比算法引入实验。

如图3-1至图3-4所示，分别为不同时间获取的山区光子计数激光雷达数据的噪声去除结果，其中，图3-1为基于改进DBSCAN方法的山区日间数据的噪声去除结果，图3-2为基于改进DBSCAN方法的山区夜间数据的噪声去除结果，图3-3为基于本实施例方法的山区日间数据的噪声去除结果，图3-4为基于本实施例方法的山区夜间数据的噪声去除结果。图中，纵坐标Elevations是高程，横坐标Along Track Distance是沿轨道距离，Noise是噪声，Canopy photons是植被光子，Ground photons是地面光子，TOC photons是林冠光子，TOC是林冠，Reference DTM是参考DTM，Reference TOC是参考林冠曲线。

结合图3-1至图3-4可以看出：本实施例方法和改进DBSCAN方法都能够提取信号光子。两种方法在处理夜间数据时取得了相似的结果。然而，在改进DBSCAN方法处理日间数据的结果中有少量信号周围的噪声光子被错分类，而本实施例方法去除了绝大多数位于地面和林冠周围的噪声光子。总得来说，由于信噪比提高，夜间数据的处理效果要优于日间数据的处理效果，本实施例方法略优于改进DBSCAN方法，去除了大多数噪声光子。

为进一步验证本实施例方法的有效性，通过提取地面点和林冠点计算地表高程和植被高度来进一步评价本实施例方法的表现，计算结果与参考数据之间的误差如图4-1至图4-8所示。

图4-1、图4-3、图4-5、图4-7是提取地面点所计算出的地表高程和参考DTM高程之间的散点图，图中，纵坐标Derived Ground Elevations是计算得到的地表高程，横坐标Referenced DTM Elevations是参考DTM高程。通过分析发现：通过算法提取的地面点计算出的地表高程和参考DTM高程具有强相关性，而且两种算法处理的结果精度均优于1m。这说明，两种算法都能够准确提取信号光子。此外，对于同一种算法，夜间结果的精度要优于日间结果的精度。这说明数据的信噪比对算法的表现有影响，当信噪比下降时，算法性能也下降。对于相同信噪比的数据来说，本实施例方法的精度均优于改进DBSCAN精度。结合图3-1的分析，这可能是因为改进DBSCAN难以确定准确的输入参数，去除那些分布在地面光子周围，与地面点具有相同邻域密度的噪声光子。尽管已经将算法优化为使用方向自适应的椭圆形邻域，但是该算法对输入参数的敏感性和统一的密度阈值仍然导致部分噪声光子被错误标记为信号。与此相对比的是，本实施例方法并没有输入参数，只依据隔离深度来标记噪声，而且在噪声去除前使用移动窗口将数据划分为小段并使用大津法自适应得提取阈值。这种做法有效避免了与改进DBSCAN一样的问题，使本实施例方法所恢复的地面高程精度更高。

图4-2、图4-4、图4-6、图4-8是提取林冠点高程与对应位置DTM高程相减得到的植被高度和CHM高度的散点图，图中，横坐标Derived Canopy Heights是参考林冠高程模型高度，纵坐标Referenced CHM Heights是计算得到的林冠高度。通过分析发现：和预期一致的是，植被高度的计算结果相对于地表高程的计算结果误差更大，相关性更弱。这是因为植被高度随着自身生长和外界影响不断变化，在不同季节、风速和时间下获取的数据的差异会比较大。分析处理结果发现，改进DBSCAN处理日间数据的精度最低，而本实施例方法处理夜间数据的精度最高。对于植被光子提取结果，尽管两种算法都受到了信噪比的影响，但是本实施例方法在不同信噪比中的表现均优于改进DBSCAN。这种现象可能由两方面原因共同导致。第一，植被光子的空间特征与其周围的噪声光子十分相似。改进DBSCAN使用椭圆邻域计算光子密度时，不准确的椭圆邻域会导致这些噪声光子与植被光子的密度是一样或相似的。然而，确定准确的椭圆邻域是困难的。与此对应的是，尽管本实施例方法同样难以准确区分植被光子及其周围的噪声，但是其计算结果并不受邻域形状的影响，算法本身并不会引入误差。第二，植被光子的空间分布与具体的植被种类、大小相关，其密度是不固定的。改进DBSCAN使用全局密度阈值，并不能适应这种变化；本实施例方法分别对每个100m窗口使用大津法计算分类阈值，能够在一定程度上改善这个问题。

综上所述：本实施例提出了一种新的光子计数激光雷达数据去噪算法，在这个基于孤立的算法中，一颗四叉树被基于光子坐标构建，以获得隔离深度。另外，为了更好地去除噪声，数据被划分为若干小段并被使用大津法计算噪声阈值。基于模拟ICESat-2数据的实验表明本实施例方法在主观和客观评价方面都优于改进DBSCAN。而且，本实施例方法无需输入参数，具有较好的鲁棒性。

Claims (3)

1.一种光子计数激光雷达数据去噪方法，其特征在于，该方法包括以下步骤：

(1)获取待处理的光子计数激光雷达数据作为原始光子数据；

(2)利用四叉树方法将原始光子数据所覆盖的空间转化为四叉树：先将原始光子数据所覆盖的空间进行四叉空间分割得到4个相同大小的子空间，根据原始光子数据的坐标将其分配在对应的子空间；然后对每个子空间中包含的光子个数进行判断，若子空间中包含的光子个数＞1，就对该子空间继续进行四叉空间分割，直至分割后的子空间中包含的光子数≤1，若子空间中包含的光子个数≤1，则该子空间停止分割；将原始光子数据所覆盖的空间作为根节点，将每次分割得到的子空间作为子空间节点，将停止分割时子空间中的光子坐标作为叶子结点，从而将原始光子数据所覆盖的空间转化为四叉树；

(3)计算原始光子数据中每个光子在所述四叉树中的隔离深度，光子在四叉树中的隔离深度等于从根节点到该光子所在叶子结点所经历的边数；

(4)将隔离深度小于深度阈值的光子标记为噪声光子，通过去除原始光子数据中的噪声光子得到原始光子数据中的信号光子。

2.根据权利要求1所述的光子计数激光雷达数据去噪方法，其特征在于，所述深度阈值利用大律法确定。

3.根据权利要求2所述的光子计数激光雷达数据去噪方法，其特征在于，采用大律法确定深度阈值的过程为：将原始光子数据划分为若干设定长度的窗口，将窗口中的光子按照隔离深度降序排列，然后按照下式确定深度阈值：

式中，n是窗口中的光子总数，t是当前假设为信号光子的光子数，ω₀(t)为当前假设的信号比例，ω₁(t)假设的噪声比例，μ₀(t)为信号平均深度，μ₁(t)为噪声平均深度，将类间方差σ²最大时的隔离深度ID(t)作为深度阈值。

Images