CN115825920B - 一种顾及冰川形态的ICESat-2光子去噪方法 - Google Patents

一种顾及冰川形态的ICESat-2光子去噪方法 Download PDF

Info

Publication number
CN115825920B
CN115825920B CN202310092563.7A CN202310092563A CN115825920B CN 115825920 B CN115825920 B CN 115825920B CN 202310092563 A CN202310092563 A CN 202310092563A CN 115825920 B CN115825920 B CN 115825920B
Authority
CN
China
Prior art keywords
photon
photons
density
neighborhood
ransac
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202310092563.7A
Other languages
English (en)
Other versions
CN115825920A (zh
Inventor
黄荣刚
常瑞杰
江利明
张文典
汪汉胜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Precision Measurement Science and Technology Innovation of CAS
Original Assignee
Institute of Precision Measurement Science and Technology Innovation of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Precision Measurement Science and Technology Innovation of CAS filed Critical Institute of Precision Measurement Science and Technology Innovation of CAS
Priority to CN202310092563.7A priority Critical patent/CN115825920B/zh
Publication of CN115825920A publication Critical patent/CN115825920A/zh
Application granted granted Critical
Publication of CN115825920B publication Critical patent/CN115825920B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Photometry And Measurement Of Optical Pulse Characteristics (AREA)

Abstract

本发明公开了一种顾及冰川形态的ICESat‑2光子去噪方法，包括以下步骤：多尺度RANSAC确定最佳光子椭圆邻域；加权密度的估算，包含连续性密度、相似性密度以及混合密度的估算；沿轨分区段自适应阈值的确定。本发明提出的光子去噪新方法，更好的解决了冰川形态尺度不一、方向变化等因素导致附着在地表附近的噪声光子密度与信号光子密度相似而难以分离的问题，能够有效提高信号光子的提取质量，进而服务于冰川表面高程探测。

Description

一种顾及冰川形态的ICESat-2光子去噪方法
技术领域
本发明属于星载激光雷达光子点云处理领域，具体涉及一种顾及冰川形态的ICESat-2光子去噪方法。
背景技术
激光测高卫星ICESat-2（The Ice, Cloud and Land Elevation Satellite-2）搭载了先进地形激光测高系统（Advanced Topography Laser Altimeter System，ATLAS），使用了全新的微脉冲多波束光子计数激光雷达，可以捕获三维高精度和高密度的光子点云。这项任务的一个主要科学目标是监测冰冻圈的表面形态及其变化。ATL03为全球地理定位光子数据产品，在晴朗的大气条件下，沿轨分辨率可达0.7米，为测量细小的冰川表面形态提供了前所未有的机会，如小裂隙、冰川湖泊等。然而，光子计数激光雷达拥有较高的灵敏度，容易受到太阳背景噪声、大气散射和仪器暗计数的影响，导致ATL03数据产品包含大量的噪声光子，且噪声光子在空间中随机广泛分布，很难从回波信号光子中剔除。因此，有必要进一步开发去噪方法，以便从ATL03数据产品中检索出具有高分辨率的冰川表面光子。
现有去噪方法基本假设目标表面光子密度比背景噪声密度大，以此提取目标表面光子。这些方法在强信噪比和简单形态区域质量较高，而在复杂形态（如表面裂隙）和低信噪比处理方面的稳健性和质量均还有待提高。此外，现有的去噪算法不能有效滤除信号光子附近的噪声光子。根据原始光子数据分布特点，我们发现：除了光子密度外，在局部空间分布上还存在其他差异，如冰川表面空间连续、密度尺度差异及方向变化等，对其深入挖掘将有助于提高光子去噪质量和适应能力。
发明内容
本发明的目的在于针对现有技术中存在的复杂形态区域光子去噪困难问题，提供一种顾及冰川形态的ICESat-2光子去噪方法。利用本发明可以最大限度的剔除附着在地表附近的噪声光子，且可以较为完整的保留复杂形态区域的信号光子，有助于在冰川变化过程的研究中提供高分辨率的表面特征，同时也为未受扰动的冰面提供了表面高度，为评估冰盖范围内的表面高程变化和质量变化提供了基础。
为了解决上述技术问题，本发明采用如下的技术方案：
一种顾及冰川形态的ICESat-2光子去噪方法，包括以下步骤：
步骤1、计算每个光子的K近邻距离；
步骤2、选定当前中心光子，使用基于圆的邻域进行多尺度RANSAC线性拟合；
步骤3、确定当前中心光子的最佳的圆邻域尺度以及椭圆邻域的长轴大小和方向；
步骤4、计算当前光子中心点的连续性密度；
步骤5、计算当前光子中心点的相似性密度；
步骤6、根据连续性密度和相似性密度计算混合加权密度；
步骤7、将沿轨光子分成不相交的沿轨区段；
步骤8、对每个沿轨区段的光子高程进行直方图统计，提取候选信号光子的海拔区间；
步骤9、根据区段内的归一化混合加权密度平均值确定自适应密度阈值；
步骤10、根据自适应密度阈值将信号光子与噪声光子进行分离。
如上所述步骤2包括以下步骤：
步骤2.1、将多尺度的集合设为
Figure SMS_1
，其中，
Figure SMS_2
为各个尺度，且每个尺度都对应着一个圆的半径，以每一个光子作为当前中心光子，搜索各个尺度下圆邻域内的光子；
步骤2.2、利用高斯混合模型以光子的K近邻距离为分类依据对当前中心光子每个尺度的邻域光子进行初步分类，以获得每个尺度下可能的信号光子集合
Figure SMS_3
；
步骤2.3、使用RANSAC算法对当前中心光子每个尺度下可能的信号光子进行RANSAC线性拟合，得到当前中心光子在每个尺度下的方向、RANSAC拟合线以及内点集合
Figure SMS_4
。
如上所述步骤3包括以下步骤：
最佳圆邻域尺度
Figure SMS_5
为步骤2中每个尺度下RANSAC拟合得到的内点数量与可能的信号光子数量之比的最大值对应的圆邻域尺度；
若当前中心光子在步骤2中得到的内点集合中，则被判定为RANSAC线性拟合有效，椭圆长轴大小为
Figure SMS_6
，椭圆长轴方向为当前中心光子在最佳圆邻域尺度下的方向，椭圆短轴为设定的固定值；
若当前中心光子不在步骤2中得到的内点集合，则被判定为RANSAC线性拟合无效，椭圆长轴方向被定义为水平向，椭圆长轴大小为设定的固定值，椭圆短轴为设定的固定值。
如上所述步骤4包括以下步骤：
当RANSAC线性拟合有效时，根据椭圆邻域内光子到RANSAC拟合线的距离分配连续性权重；
当RANSAC线性拟合无效时，光子的连续性权重分配为零，
统计椭圆邻域中所有光子点的连续性权重之和作为当前光子中心点的连续性密度。
如上所述步骤5包括以下步骤：
步骤5.1、为当前中心光子的椭圆邻域内的光子分配相似性权重：
为椭圆邻域内的光子搜索
Figure SMS_7
个最近的光子；计算椭圆邻域内的光子与其
Figure SMS_8
个最近的光子之间的距离的平均值作为光子间距
Figure SMS_9
；评估椭圆邻域内的光子和当前中心光子之间的光子间距之差，为椭圆邻域内的光子分配一个相似性权重；
步骤5.2、统计椭圆邻域中所有光子点的相似性权重之和作为当前光子中心点的相似性密度。
如上所述步骤6中的混合加权密度基于以下公式获得：
Figure SMS_11
式中：
Figure SMS_13
和
Figure SMS_17
分别是连续性密度系数和相似性密度系数，
Figure SMS_12
为连续性密度，
Figure SMS_15
为相似性密度，
Figure SMS_16
为当前中心光子的椭圆邻域的椭圆长轴，
Figure SMS_18
为当前中心光子的椭圆邻域的椭圆短轴，
Figure SMS_10
是当前中心光子
Figure SMS_14
的归一化混合加权密度。
如上所述步骤8中候选信号光子的海拔区间为
Figure SMS_19
，其中，
Figure SMS_20
为光子高程直方图最大频率区间的最小高程值，
Figure SMS_21
为光子高程直方图最大频率区间的最大高程值，
Figure SMS_22
为扩展值。
如上所述步骤9中自适应密度阈值基于以下公式：
Figure SMS_23
式中：
Figure SMS_24
为自适应密度阈值，
Figure SMS_25
为区段内的归一化混合加权密度平均值，
Figure SMS_26
为密度校正值，
Figure SMS_27
为区段内的RANSAC线性拟合的有效拟合比例，
Figure SMS_28
为RANSAC线性拟合的有效拟合比例阈值。
如上所述步骤4中椭圆邻域内光子的连续性权重基于以下公式获得：
Figure SMS_30
式中：
Figure SMS_32
为椭圆邻域内光子的连续性权重，
Figure SMS_34
为第一高斯核函数标准差，
Figure SMS_31
为光子
Figure SMS_33
到RANSAC拟合线的距离，
Figure SMS_35
是RANSAC线性拟合的有效拟合指标，若能有效拟合，则
Figure SMS_36
，若为无效拟合，则
Figure SMS_29
。
如上所述步骤5中椭圆邻域内光子的相似性权重基于以下公式获得：
Figure SMS_37
其中，
Figure SMS_38
为椭圆邻域内光子的相似性权重，
Figure SMS_39
为第二高斯核函数标准差，
Figure SMS_40
为椭圆邻域内的光子与当前中心光子的空间分布差异。
与现有技术相比，本发明具有以下有益效果：
（1）高精度性：本发明通过多尺度RANSAC确定光子点的最佳椭圆邻域，并采用一维高斯核进行连续性密度、相似性密度以及混合密度的估算，更好的解决了冰川形态尺度不一、方向变化等因素导致附着在地表附近的噪声光子密度与信号光子密度相似而难以分离的问题，能够有效提高信号光子的提取质量，进而服务于冰川表面高程探测；
（2）实用性：本发明可以实现不同冰川地貌形态的高精度光子去噪，且可以获取大多数场景的地表坡度变化，为表面裂隙三维特征的提取提供基础。
附图说明
图1为本发明的流程示意图。
图2（a）为雅各布港冰川表面的复杂形态区域一的光子去噪结果图，其中
Figure SMS_42
和
Figure SMS_45
分别为弱波束原始数据图和强波束原始数据图，
Figure SMS_47
和
Figure SMS_43
分别为弱波束混合加权密度图和强波束混合加权密度图，
Figure SMS_44
和
Figure SMS_46
分别为弱波束去噪结果图和强波束去噪结果图，
Figure SMS_48
和
Figure SMS_41
分别为弱波束去噪结果精度评估图和强波束去噪结果精度评估图；
图2（b）为雅各布港冰川表面的复杂形态区域二的光子去噪结果图，其中
Figure SMS_50
和
Figure SMS_53
分别为弱波束原始数据图和强波束原始数据图，
Figure SMS_54
和
Figure SMS_51
分别为弱波束混合加权密度图和强波束混合加权密度图，
Figure SMS_52
和
Figure SMS_55
分别为弱波束去噪结果图和强波束去噪结果图，
Figure SMS_56
和
Figure SMS_49
分别为弱波束去噪结果精度评估图和强波束去噪结果精度评估图；
图2（c）为雅各布港冰川表面的斜坡区域的光子去噪结果图，其中
Figure SMS_59
和
Figure SMS_61
分别为弱波束原始数据图和强波束原始数据图，
Figure SMS_63
和
Figure SMS_58
分别为弱波束混合加权密度图和强波束混合加权密度图，
Figure SMS_60
和
Figure SMS_62
分别为弱波束去噪结果图和强波束去噪结果图，
Figure SMS_64
和
Figure SMS_57
分别为弱波束去噪结果精度评估图和强波束去噪结果精度评估图；
图2（d）为雅各布港冰川表面的平缓区域的光子去噪结果图，其中
Figure SMS_66
和
Figure SMS_68
分别为弱波束原始数据图和强波束原始数据图，
Figure SMS_71
和
Figure SMS_67
分别为弱波束混合加权密度图和强波束混合加权密度图，
Figure SMS_69
和
Figure SMS_70
分别为弱波束去噪结果图和强波束去噪结果图，
Figure SMS_72
和
Figure SMS_65
分别为弱波束去噪结果精度评估图和强波束去噪结果精度评估图。
图3（a）为复杂形态区域一的去噪结果分别与RANSAC算法、改进的DBSCAN算法和KNN算法的去噪结果的对比图，其中
Figure SMS_75
和
Figure SMS_76
分别为弱波束RANSAC去噪结果图和强波束RANSAC去噪结果图，
Figure SMS_78
和
Figure SMS_74
分别为弱波束改进的DBSCAN去噪结果图和强波束改进的DBSCAN去噪结果图，
Figure SMS_77
和
Figure SMS_79
分别为弱波束KNN去噪结果图和强波束KNN去噪结果图，
Figure SMS_80
和
Figure SMS_73
分别为弱波束本发明的去噪结果图和强波束本发明的去噪结果图；
图3（b）为复杂形态区域二的去噪结果分别与RANSAC算法、改进的DBSCAN算法和KNN算法的去噪结果的对比图，其中
Figure SMS_82
和
Figure SMS_85
分别为弱波束RANSAC去噪结果图和强波束RANSAC去噪结果图，
Figure SMS_86
和
Figure SMS_83
分别为弱波束改进的DBSCAN去噪结果图和强波束改进的DBSCAN去噪结果图，
Figure SMS_84
和
Figure SMS_87
分别为弱波束KNN去噪结果图和强波束KNN去噪结果图，
Figure SMS_88
和
Figure SMS_81
分别为弱波束本发明的去噪结果图和强波束本发明的去噪结果图；
图3（c）为斜坡区域的去噪结果分别与RANSAC算法、改进的DBSCAN算法和KNN算法的去噪结果的对比图，其中
Figure SMS_90
和
Figure SMS_93
分别为弱波束RANSAC去噪结果图和强波束RANSAC去噪结果图，
Figure SMS_95
和
Figure SMS_91
分别为弱波束改进的DBSCAN去噪结果图和强波束改进的DBSCAN去噪结果图，
Figure SMS_92
和
Figure SMS_94
分别为弱波束KNN去噪结果图和强波束KNN去噪结果图，
Figure SMS_96
和
Figure SMS_89
分别为弱波束本发明的去噪结果图和强波束本发明的去噪结果图；
图3（d）为平缓区域的去噪结果分别与RANSAC算法、改进的DBSCAN算法和KNN算法的去噪结果的对比图，其中
Figure SMS_99
和
Figure SMS_100
分别为弱波束RANSAC去噪结果图和强波束RANSAC去噪结果图，
Figure SMS_102
和
Figure SMS_98
分别为弱波束改进的DBSCAN去噪结果图和强波束改进的DBSCAN去噪结果图，
Figure SMS_101
和
Figure SMS_103
分别为弱波束KNN去噪结果图和强波束KNN去噪结果图，
Figure SMS_104
和
Figure SMS_97
分别为弱波束本发明的去噪结果图和强波束本发明的去噪结果图。
图4为本发明的去噪结果与ATL03置信度光子、标准陆冰沿轨高度产品ATL06的对比图。其中，（a）和（c）分别为复杂形态区域一弱波束ATL03置信度光子图和强波束ATL03置信度光子图，（b）和（d）分别为复杂形态区域一弱波束本发明的去噪结果与ATL06的对比图和强波束本发明的去噪结果与ATL06的对比图，（e）和（g）分别为复杂形态区域二弱波束ATL03置信度光子图和强波束ATL03置信度光子图，（f）和（h）分别为复杂形态区域二弱波束本发明的去噪结果与ATL06的对比图和强波束本发明的去噪结果与ATL06的对比图。
具体实施方式
为了便于本领域普通技术人员理解和实施本发明，下面结合实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。
实施例
一种顾及冰川形态的ICESat-2光子去噪方法，具体包括以下步骤：
步骤一：计算光子的K近邻距离。针对ICESat-2 ATL03光子数据，采用KdTree构建光子空间索引，并计算每个光子的K近邻距离（KNN）。
Figure SMS_105
公式（1）
式中：
Figure SMS_108
表示光子数据中的第
Figure SMS_112
个光子到周围最近的
Figure SMS_113
个光子的总距离，
Figure SMS_109
和
Figure SMS_111
分别表示第
Figure SMS_116
个光子和第
Figure SMS_117
个光子的沿轨距离，
Figure SMS_106
和
Figure SMS_110
分别表示第
Figure SMS_114
个光子和第
Figure SMS_115
个光子的光子高程。
Figure SMS_107
是计算
Figure SMS_118
的重要参数，即选取近邻光子的个数，实质上是一个平滑因子，其作用是平滑同类别光子内部的空间异质性。若
Figure SMS_119
过小，则平滑效果不明显，信号光子与噪声光子的距离和会比较相近，难以进行区分；若
Figure SMS_120
过大，则将噪声光子与信号光子之间的差异过度平滑，使得光子去噪不完全。
步骤二：多尺度RANSAC拟合。求取了每个光子的K近邻距离后，再进行多尺度RANSAC线性拟合，由于使用了基于圆的邻域，每一个光子的步骤如下。
首先，将多尺度的集合设为
Figure SMS_121
，其中，
Figure SMS_122
为各个尺度，且每个尺度都对应着一个圆的半径。以光子数据中的每一个光子作为当前中心光子，搜索其各个尺度下圆邻域内的光子。
然后，为了减小背景噪声对拟合结果的影响，利用高斯混合模型（GMM）以光子的K近邻距离为分类依据对当前中心光子每个尺度的邻域光子进行初步分类，以获得每个尺度下可能的信号光子集合
Figure SMS_123
。
最后，使用RANSAC算法对当前中心光子每个尺度下可能的信号光子进行RANSAC线性拟合，得到当前中心光子在每个尺度下的方向、RANSAC拟合线以及内点集合
Figure SMS_124
。
步骤三：确定当前中心光子的最佳的圆邻域尺度以及椭圆邻域的长轴大小和方向。最佳圆邻域尺度
Figure SMS_125
为步骤二中每个尺度下RANSAC拟合得到的内点数量与可能的信号光子数量之比的最大值对应的圆邻域尺度，比值越大，则该尺度下RANSAC拟合得到的内点为信号光子的概率越大，拟合结果越好。考虑到信号光子在水平方向密集分布，为抑制垂直方向上噪声光子对后续信号光子密度估算的影响，将当前中心光子的圆邻域变换为椭圆邻域。
若当前中心光子在步骤二中得到的内点集合中，则被判定为RANSAC线性拟合有效，椭圆长轴大小为
Figure SMS_126
，椭圆长轴方向为当前中心光子在最佳圆邻域尺度下的方向。
若当前中心光子不在步骤二中得到的内点集合，则被判定为RANSAC线性拟合无效，在这种情况下，椭圆长轴方向被定义为水平向，椭圆长轴大小被设定为一个设定的固定值（如20m）。此外，无论中心光子是否在步骤二中得到的内点集合，椭圆短轴都设为固定值。
Figure SMS_127
公式（2）
式中：
Figure SMS_130
为
Figure SMS_132
尺度上拟合得到的内点
Figure SMS_133
的数量，
Figure SMS_129
为
Figure SMS_131
尺度上高斯混合模型得到的可能的信号光子
Figure SMS_134
的数量，
Figure SMS_136
为尺度序号,
Figure SMS_128
为具有最高拟合比例的尺度，也即
Figure SMS_135
最大对应的尺度，argmax为求函数参数的函数。
步骤四：当前光子中心点的连续性密度估算。当前中心光子的椭圆邻域是根据步骤三自适应确定的，当RANSAC线性拟合有效时，可以认为椭圆邻域中局部表面可能的信号光子是平滑的、连续的和线性的，利用一维高斯核函数，根据椭圆邻域内光子到RANSAC拟合线的距离分配连续性权重，距离越小，权重越大。此外，当RANSAC线性拟合无效时，局部表面不光滑或不存在，光子的连续性权重就直接分配为零，如公式（3）。
Figure SMS_137
公式（3）
式中：
Figure SMS_138
表示椭圆邻域内光子被分配到的连续性权重，
Figure SMS_139
为第一高斯核函数标准差，
Figure SMS_140
为光子
Figure SMS_141
到RANSAC拟合线的距离，
Figure SMS_142
是RANSAC线性拟合的有效拟合指标，若能有效拟合，则
Figure SMS_143
，若为无效拟合，则
Figure SMS_144
。
在给椭圆邻域中所有的光子分配了连续性权重后，基于连续性的加权密度可以根据公式（4）进行计算，即对于每个当前中心光子，统计椭圆邻域中所有光子点的连续性权重之和。
Figure SMS_145
公式（4）
式中：光子
Figure SMS_146
来自于当前中心光子
Figure SMS_147
椭圆邻域内的光子集合
Figure SMS_148
，
Figure SMS_149
为椭圆邻域内光子点
Figure SMS_150
的连续性权重，
Figure SMS_151
为当前光子中心点
Figure SMS_152
的连续性密度。
步骤五：当前光子中心点的相似性密度估算。在当前中心光子的椭圆邻域内，光子的属性可能是不一致的，是信号光子和噪声光子的混合体。这样一来，信号光子和附近噪声光子的密度估计就会受到严重影响，它们的密度差异就会变小。因此，在密度估计中给不同类别的光子分配不同的权重，将改善信号和附近噪声光子之间的密度差。一般来说，具有相同类别属性的光子的空间分布在一个局部区域是近似一致的。也就是说，一个信号光子的光子间距与相邻的信号光子的光子间距大致相同，相邻的噪声也应遵守同样的规律。因此，把光子间距作为相似性指标，用当前中心光子的椭圆邻域内的光子与当前中心光子的光子间距之差来评估权重。对于一个当前中心光子
Figure SMS_153
，为其椭圆邻域内的光子分配相似性权重的步骤如下。1）为椭圆邻域内的光子
Figure SMS_154
搜索
Figure SMS_155
个最近的光子；2）计算椭圆邻域内的光子
Figure SMS_156
与其
Figure SMS_157
个最近的光子之间的距离的平均值作为光子间距
Figure SMS_158
；3）评估椭圆邻域内的光子和当前中心光子
Figure SMS_159
之间的光子间距之差，并通过一维高斯核函数为椭圆邻域内的光子分配一个相似性权重，如公式（5）所示。
Figure SMS_160
公式（5）
式中：
Figure SMS_161
表示椭圆邻域内光子被分配到的相似性权重，
Figure SMS_162
为第二高斯核函数标准差，其中，第二高斯核函数标准差与第一高斯核函数标准差仅在数值上有差异，
Figure SMS_163
为椭圆邻域内的光子
Figure SMS_164
与当前中心光子
Figure SMS_165
的空间分布差异,即光子间距之差。
在给椭圆邻域中所有的光子点分配了相似性权重后，可以根据公式（6）计算出当前中心光子
Figure SMS_166
基于相似性的加权密度。
Figure SMS_167
公式（6）
式中：光子
Figure SMS_168
来自于当前中心光子
Figure SMS_169
椭圆邻域内的光子集合
Figure SMS_170
，
Figure SMS_171
为椭圆邻域内光子点
Figure SMS_172
的相似性权重，
Figure SMS_173
为当前光子中心点
Figure SMS_174
的相似性密度。
步骤六：混合加权密度估算。连续性密度反映了局部表面的形状特征，而相似度密度则表达了光子点间的空间分布差异。因此，将两类加权密度结合起来作为混合加权密度。混合加权密度有利于放大信号光子和噪声光子之间的差异。此外，由于椭圆邻域的长轴大小对于不同区域的光子可能不同，混合加权密度应当进行归一化处理。
Figure SMS_175
公式（7）
式中：
Figure SMS_178
和
Figure SMS_179
分别是连续性密度系数和相似性密度系数，
Figure SMS_181
为当前中心光子
Figure SMS_177
的椭圆邻域的椭圆长轴，
Figure SMS_180
为当前中心光子
Figure SMS_182
的椭圆邻域的椭圆短轴，
Figure SMS_183
是当前中心光子
Figure SMS_176
的归一化混合加权密度。
步骤七：沿轨光子的分段。为了克服信号光子和噪声光子之间的区域密度差异对阈值选择的影响，将沿轨光子分成不相交的沿轨区段。这样一来，信号光子密度、噪声光子密度和两种密度的差值将在一个沿轨区段内保持较好的一致性。此外，局部表面形状将相对简单，信号光子的海拔也将更加集中。
步骤八：候选信号光子的提取。沿轨光子进行分段后，对每个沿轨区段的光子高程进行直方图统计，以提取候选信号光子的海拔区间。对于一个沿轨区段，首先生成一个基于光子高程的直方图，其中，横坐标为沿轨距离，纵坐标为光子高程。然后，选择光子高程区间
Figure SMS_184
的光子作为候选信号光子，其中，
Figure SMS_185
为光子高程直方图最大频率区间的最小高程值，
Figure SMS_186
为光子高程直方图最大频率区间的最大高程值，最大频率区间即为光子数量最多的区间。为了确保所有的信号光子都包含在海拔区间内，设定一个扩展值
Figure SMS_187
对海拔区间进行扩展，最终候选信号光子的海拔区间为
Figure SMS_188
。
步骤九：自适应密度阈值的确定。一般来说，每一个沿轨区段中的信号光子和噪声光子都存在密度差，密度阈值
Figure SMS_190
可以设置为沿轨区段内所有光子密度的平均值。然而，对于不同的沿轨区段，信号光子和噪声光子之间的密度差是不同的。为了提高
Figure SMS_194
的鲁棒性，在密度阈值中加入一个校正值
Figure SMS_197
，以适应不同数据质量的沿轨区段，如公式（8）所示。在本发明中，当大部分光子和其椭圆邻域内的光子能够满足步骤三中的RANSAC有效拟合时，该数据就具有良好的质量。因此，在沿轨区段内计算有效拟合比率
Figure SMS_191
，用于描述沿轨区段中有多少光子与其椭圆邻域内的光子可以被有效拟合。当
Figure SMS_192
大于阈值
Figure SMS_195
时，说明该沿轨区段的数据质量较好，为了更好地去除信号光子周围的噪声，密度阈值可以设置为一个较高的值，因此加上一个校正值
Figure SMS_198
。相反，当
Figure SMS_189
小于阈值
Figure SMS_193
时，数据质量不被认为是好的，而且信噪比可能很低。因此，为了保证信号光子的充分提取，应该设置一个较低的密度阈值，将密度阈值减去一个校正值
Figure SMS_196
。
Figure SMS_199
公式（8）
式中：
Figure SMS_200
为自适应密度阈值，
Figure SMS_201
为区段内的归一化混合加权密度平均值，
Figure SMS_202
为密度校正值，
Figure SMS_203
为区段内的RANSAC线性拟合的有效拟合比例，
Figure SMS_204
为有效拟合比例判定数据质量情况的RANSAC线性拟合的有效拟合比例阈值。
步骤十：信号光子与噪声光子的分离。在得到每一个沿轨区段的自适应密度阈值后，将每一个沿轨区段内光子密度大于自适应密度阈值的光子判定为信号光子；将光子密度小于自适应密度阈值的光子判定为噪声光子。
为说明本发明的有效性，选择了雅各布港冰川表面复杂形态区域、斜坡区域、平缓区域的强弱信号数据进行测试。表1为参数设置说明表。
表1参数设置说明表
Figure SMS_205
本发明提供了精度更高的去噪结果，且能够处理不同质量的数据，在多种地形结构的F1分数均可以达到0.9以上。不同区域去噪结果总体评价如表2所示。
表2本发明去噪结果精度评估表
Figure SMS_206
相比较于现有方法，本方法对附着在地表附近的噪声光子拥有更好的剔除效果，这将提供更为精确的高分辨率地表光子，从而获取更为准确的地面高度信息。表3-5分别为基于RANSAC算法、改进的DBSCAN算法、KNN算法的去噪结果的精度评估表。
表3基于RANSAC算法的去噪结果的精度评估表
Figure SMS_207
表4基于改进的DBSCAN算法的去噪结果的精度评估表
Figure SMS_208
表5基于KNN算法的去噪结果的精度评估表
Figure SMS_209
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (4)

1.一种顾及冰川形态的ICESat-2光子去噪方法，其特征在于，包括以下步骤：
步骤1、计算每个光子的K近邻距离；
步骤2、选定当前中心光子，使用基于圆的邻域进行多尺度RANSAC线性拟合；
步骤3、确定当前中心光子的最佳的圆邻域尺度以及椭圆邻域的长轴大小和方向；
步骤4、计算当前光子中心点的连续性密度；
步骤5、计算当前光子中心点的相似性密度；
步骤6、根据连续性密度和相似性密度计算混合加权密度；
步骤7、将沿轨光子分成不相交的沿轨区段；
步骤8、对每个沿轨区段的光子高程进行直方图统计，提取候选信号光子的海拔区间；
步骤9、根据区段内的归一化混合加权密度平均值确定自适应密度阈值；
步骤10、根据自适应密度阈值将信号光子与噪声光子进行分离，
所述步骤2包括以下步骤：
步骤2.1、将多尺度的集合设为
Figure QLYQS_1
，其中，
Figure QLYQS_2
为各个尺度，且每个尺度都对应着一个圆的半径，以每一个光子作为当前中心光子，搜索各个尺度下圆邻域内的光子；
步骤2.2、利用高斯混合模型以光子的K近邻距离为分类依据对当前中心光子每个尺度的邻域光子进行初步分类，以获得每个尺度下可能的信号光子集合
Figure QLYQS_3
；
步骤2.3、使用RANSAC算法对当前中心光子每个尺度下可能的信号光子进行RANSAC线性拟合，得到当前中心光子在每个尺度下的方向、RANSAC拟合线以及内点集合
Figure QLYQS_4
，
所述步骤3包括以下步骤：
最佳圆邻域尺度
Figure QLYQS_5
为步骤2中每个尺度下RANSAC拟合得到的内点数量与可能的信号光子数量之比的最大值对应的圆邻域尺度；
若当前中心光子在步骤2中得到的内点集合中，则被判定为RANSAC线性拟合有效，椭圆长轴大小为
Figure QLYQS_6
，椭圆长轴方向为当前中心光子在最佳圆邻域尺度下的方向，椭圆短轴为设定的固定值；
若当前中心光子不在步骤2中得到的内点集合，则被判定为RANSAC线性拟合无效，椭圆长轴方向被定义为水平向，椭圆长轴大小为设定的固定值，椭圆短轴为设定的固定值，
所述步骤4包括以下步骤：
当RANSAC线性拟合有效时，根据椭圆邻域内光子到RANSAC拟合线的距离分配连续性权重；
当RANSAC线性拟合无效时，光子的连续性权重分配为零，
统计椭圆邻域中所有光子点的连续性权重之和作为当前光子中心点的连续性密度，
所述步骤4中椭圆邻域内光子的连续性权重基于以下公式获得：
Figure QLYQS_8
式中：
Figure QLYQS_10
为椭圆邻域内光子的连续性权重，
Figure QLYQS_12
为第一高斯核函数标准差，
Figure QLYQS_9
为光子
Figure QLYQS_11
到RANSAC拟合线的距离，
Figure QLYQS_13
是RANSAC线性拟合的有效拟合指标，若能有效拟合，则
Figure QLYQS_14
，若为无效拟合，则
Figure QLYQS_7
，
所述步骤5包括以下步骤：
步骤5.1、为当前中心光子的椭圆邻域内的光子分配相似性权重：
为椭圆邻域内的光子搜索
Figure QLYQS_15
个最近的光子；计算椭圆邻域内的光子与其
Figure QLYQS_16
个最近的光子之间的距离的平均值作为光子间距
Figure QLYQS_17
；评估椭圆邻域内的光子和当前中心光子之间的光子间距之差，为椭圆邻域内的光子分配一个相似性权重；
步骤5.2、统计椭圆邻域中所有光子点的相似性权重之和作为当前光子中心点的相似性密度，
所述步骤5中椭圆邻域内光子的相似性权重基于以下公式获得：
Figure QLYQS_18
其中，
Figure QLYQS_19
为椭圆邻域内光子的相似性权重，
Figure QLYQS_20
为第二高斯核函数标准差，
Figure QLYQS_21
为椭圆邻域内的光子与当前中心光子的空间分布差异。
2.根据权利要求1所述一种顾及冰川形态的ICESat-2光子去噪方法，其特征在于，所述步骤6中的混合加权密度基于以下公式获得：
Figure QLYQS_24
式中：
Figure QLYQS_27
和
Figure QLYQS_29
分别是连续性密度系数和相似性密度系数，
Figure QLYQS_23
为连续性密度，
Figure QLYQS_26
为相似性密度，
Figure QLYQS_28
为当前中心光子的椭圆邻域的椭圆长轴，
Figure QLYQS_30
为当前中心光子的椭圆邻域的椭圆短轴，
Figure QLYQS_22
是当前中心光子
Figure QLYQS_25
的归一化混合加权密度。
3.根据权利要求2所述一种顾及冰川形态的ICESat-2光子去噪方法，其特征在于，所述步骤8中候选信号光子的海拔区间为
Figure QLYQS_31
，其中，
Figure QLYQS_32
为光子高程直方图最大频率区间的最小高程值，
Figure QLYQS_33
为光子高程直方图最大频率区间的最大高程值，
Figure QLYQS_34
为扩展值。
4.根据权利要求3所述一种顾及冰川形态的ICESat-2光子去噪方法，其特征在于，所述步骤9中自适应密度阈值基于以下公式：
Figure QLYQS_35
式中：
Figure QLYQS_36
为自适应密度阈值，
Figure QLYQS_37
为区段内的归一化混合加权密度平均值，
Figure QLYQS_38
为密度校正值，
Figure QLYQS_39
为区段内的RANSAC线性拟合的有效拟合比例，
Figure QLYQS_40
为RANSAC线性拟合的有效拟合比例阈值。
CN202310092563.7A 2023-02-10 2023-02-10 一种顾及冰川形态的ICESat-2光子去噪方法 Active CN115825920B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310092563.7A CN115825920B (zh) 2023-02-10 2023-02-10 一种顾及冰川形态的ICESat-2光子去噪方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310092563.7A CN115825920B (zh) 2023-02-10 2023-02-10 一种顾及冰川形态的ICESat-2光子去噪方法

Publications (2)

Publication Number Publication Date
CN115825920A CN115825920A (zh) 2023-03-21
CN115825920B true CN115825920B (zh) 2023-05-05

Family

ID=85520954

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310092563.7A Active CN115825920B (zh) 2023-02-10 2023-02-10 一种顾及冰川形态的ICESat-2光子去噪方法

Country Status (1)

Country Link
CN (1) CN115825920B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116243273B (zh) * 2023-05-09 2023-09-15 中国地质大学（武汉） 一种针对植被冠层提取的光子计数激光雷达数据滤波方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP5520786B2 (ja) * 2010-11-11 2014-06-11 株式会社パスコ レーザ密度分布推定装置、レーザ密度分布推定方法、及びプログラム
CN109799494B (zh) * 2017-11-17 2020-09-29 中国林业科学研究院资源信息研究所 一种星载光子计数激光雷达数据快速去噪滤波方法
CN111665517B (zh) * 2020-05-29 2022-11-18 同济大学 一种基于密度统计的单光子激光测高数据去噪方法及装置
CN112986964B (zh) * 2021-02-26 2023-03-31 北京空间机电研究所 基于噪声邻域密度的光子计数激光点云自适应去噪方法
CN114355367B (zh) * 2022-01-10 2024-10-01 中国人民解放军61540部队 一种基于星载单光子激光雷达数据测量浅海水深的方法
CN114779215A (zh) * 2022-04-21 2022-07-22 昆明理工大学 一种植被覆盖区星载光子计数激光雷达数据去噪方法
CN115222949B (zh) * 2022-09-21 2022-12-06 自然资源部第一海洋研究所 一种基于激光卫星数据的浅海区域光子去噪方法

Also Published As

Publication number Publication date
CN115825920A (zh) 2023-03-21

Similar Documents

Publication Publication Date Title
CN108957453B (zh) 一种基于多目标跟踪的高精度动目标成像及识别方法
CN111830506B (zh) 一种基于K-means聚类算法的海面风速方法
US7787657B2 (en) SAR ATR treeline extended operating condition
CN108681525B (zh) 一种基于车载激光扫描数据的路面点云强度增强方法
US7116265B2 (en) Recognition algorithm for the unknown target rejection based on shape statistics obtained from orthogonal distance function
CN111665517B (zh) 一种基于密度统计的单光子激光测高数据去噪方法及装置
CN108171193B (zh) 基于超像素局部信息度量的极化sar舰船目标检测方法
CN115825920B (zh) 一种顾及冰川形态的ICESat-2光子去噪方法
CN112669333B (zh) 一种单木信息提取方法
CN106952242A (zh) 一种基于体素的渐进不规则三角网点云滤波方法
CN109100697B (zh) 一种基于地面监视雷达系统的目标凝聚方法
CN113256990B (zh) 基于聚类算法的雷达采集道路车辆信息的方法及系统
CN113466827B (zh) 一种基于改进局部稀疏算法的去噪方法
CN113281716B (zh) 一种光子计数激光雷达数据去噪方法
CN113034569A (zh) 一种基于点云数据的船舶超限预警方法及系统
CN105138992A (zh) 一种基于区域主动轮廓模型的海岸线检测方法
CN108983194B (zh) 一种基于地面监视雷达系统的目标提取及凝聚方法
CN116626685B (zh) 基于机器学习的河道底泥实时监测方法及系统
CN117437147A (zh) 一种浅水域单光子点云的去噪方法
CN111080647B (zh) 基于自适应滑动窗口滤波和fcm的sar图像分割方法
CN116165635A (zh) 多级滤波算法的日间条件下不同波束光子云数据去噪方法
CN115436966A (zh) 一种激光雷达参考水深控制点批量提取方法
CN114779235A (zh) 基于车载毫米波雷达的道路边界检测方法和系统
Zhao et al. Automatic extraction of floating ice at Antarctic continental margin from remotely sensed imagery using object-based segmentation
Ma et al. Ocean SAR image segmentation and edge gradient feature extraction

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant