一种激光雷达回波检测海底地形信号异常点去除方法

1.本发明属于信号处理领域，尤其涉及一种激光雷达回波检测海底地形信号异常点去除方法。


背景技术：

2.美国国家航空航天局(nasa)于2018年9月15日发射了ice-sat-2，其携带了第一颗星载光子计数激光雷达atlas。由于探测器灵敏，在atlas光学特性的影响下，原始光子数据包含大量太阳背景噪声光子，在采用一定的方法检测出海底反射的光子信号后，依然残留大量噪声信号靠近海表面以及海底信号周围。
3.离群点检测的目的是过滤每个数据集中的异常数据，因此离群点检测的主要目标是发现与正常情况明显不同或不一致的离群点或噪声。由于激光雷达检测到的海底光子信号仍然包含少量噪声光子，为了解决这些问题，需要针对当前场景下的异常点处理方法，将多余噪声点去除，以便进行后续信号处理。


技术实现要素：

4.为了去除激光雷达回波检测海底地形信号中的噪声点，本发明提供一种激光雷达回波检测海底地形信号异常点去除方法。
5.本发明的目的通过如下的技术方案来实现：
6.一种激光雷达回波检测海底地形信号异常点去除方法，其特征在于，该方法包括如下步骤：
7.s1：读取激光雷达检测海底地形结果信号数据，横轴为沿轨距离，纵轴为高程，当前数据段中当前数据段最小沿轨距离和最大沿轨距离分别是x
min
和x
max
，当前数据段最小和最大高程分别是y
min
和y
max
；
8.s2：对激光雷达检测海底地形结果信号数据进行小波滤波，小波滤波能对信号进行初步的噪声去除，根据小波分解各层信号进行噪声水平估计，采用硬阈值进行优化；
9.s3：针对距离特征近似的信号簇进行分类处理，沿沿轨距离方向采用聚类方法将信号聚类；
10.s4：对每一类数据，设置沿轨距离滑动窗，计算每个滑动窗中沿轨距离数据向量a的绝对偏差中位值mad1，将沿轨数据向量上大于mad对应的数据点去除；其中mad1的计算公式为：
11.mad1＝c
×
median(|a-median(a)|)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
12.其中erfcinv()是逆补余误差函数，median()是取中位数函数；
13.s5：将s4得到的所有类别的结果数据汇总后，沿沿轨距离方向升序重新排列；
14.s6：对s5获得的数据设置高程滑动窗，计算每个高程滑动窗中高程数据向量b的mad2，将沿高程方向上大于mad2的对应的数据点去除
15.mad2＝c
×
median(|b-median(b)|)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
16.进一步地，常见聚类分类方法有k-means,但所述s3中采用的k-medoids算法中，将从当前簇中选择这样一个点作为中心点，该点是从当前簇到所有其他点的距离的最小和，这使得每一个簇大小不会有很大变化，能够更好分类，使用该方法将聚类数据为3类。
17.进一步地，由于沿轨方向的异常点信号多于高程方向上的异常点信号，滑动窗越短，绝对偏差中位值越小，异常点范围越大，因此所述沿轨距离滑动窗的窗长小于所述高程滑动窗的窗长。
18.进一步地，所述沿轨距离滑动窗的窗长为n1，高程滑动窗的窗长为n2:
19.n2＝n1
×
0.1
×
(x
max-x
min
)/(y
max-y
min
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
20.本发明的有益效果是：
21.由于激光雷达检测到的海底光子信号仍然包含少量噪声光子，主要集中在海面和海底信号附近，本发明的方法针对当前场景下，采用横坐标与纵坐标结合的异常点处理方法，将多余噪声点去除，提高了后续信号处理精度。
附图说明
22.图1是本发明方法实施例的流程图；
23.图2是采用本发明方法后的结果对比图，其中左图为原始信号结果，右图为采用本发明所述方法后信号结果。
具体实施方式
24.下面根据附图和优选实施例详细描述本发明，本发明的目的和效果将变得更加明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
25.如图1所示，本发明的激光雷达回波检测海底地形信号异常点去除方法，包括如下步骤：
26.s1：读取激光雷达检测海底地形结果信号数据，横轴为沿轨距离，纵轴为高程，当前数据段最小沿轨距离和最大沿轨距离分别是x
min
和x
max
，当前数据段最小和最大高程分别是y
min
和y
max
；
27.s2：对激光雷达检测海底地形结果信号数据进行小波滤波，小波滤波能对信号进行初步的噪声去除，根据小波分解各层信号进行噪声水平估计，采用硬阈值进行优化；
28.s3：针对距离特征近似的信号簇进行分类处理，沿沿轨距离方向采用聚类方法将信号聚类；
29.s4：对每一类数据，设置沿轨距离滑动窗，计算每个滑动窗中沿轨距离数据向量a的绝对偏差中位值mad1，将沿轨数据向量上大于mad对应的数据点去除；其中mad1的计算公式为：
30.mad1＝c
×
median(|a-median(a)|)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
31.其中erfcinv()是逆补余误差函数，median()是取中位数函数；
32.s5：将s4得到的所有类别的结果数据汇总后，沿沿轨距离方向升序重新排列；
33.s6：对s5获得的数据设置高程滑动窗，计算每个高程滑动窗中高程数据向量b的mad2，将沿高程方向上大于mad2的对应的数据点去除
34.mad2＝c
×
median(|b-median(b)|)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
35.常见聚类分类方法有k-means等,但该实施例中，s3采用的k-medoids算法，将从当前簇中选择这样一个点作为中心点，该点是从当前簇到所有其他点的距离的最小和，这使得每一个簇大小不会有很大变化，能够更好分类，使用该方法将聚类数据为3类。
36.由于沿轨方向的异常点信号多于高程方向上的异常点信号，滑动窗越短，绝对偏差中位值越小，异常点范围越大，因此，该实施例中，优选沿轨距离滑动窗的窗长小于所述高程滑动窗的窗长。
37.所述沿轨距离滑动窗的窗长为n1，高程滑动窗的窗长为n2:
38.n2＝n1
×
0.1
×
(x
max-x
min
)/(y
max-y
nin
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
39.图2中左图是基于icesat-2 atl03_20190221014318_08340207_004_01.h5原始文件提取出来的st.thomas岛屿附近的水深，右图为采用本发明的方法进行异常点去除后的地形图，从图中可以看到，原始文件提取出来的地形回波信号中包含了大量出了海底地形外的噪声信号，如果保留会干扰有效信号判断，采用本发明的方法处理后，噪声光子几乎都被去除，只剩下有效地形信号，对后续回波信号处理具有极大益处。
40.本领域普通技术人员可以理解，以上所述仅为发明的优选实例而已，并不用于限制发明，尽管参照前述实例对发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在发明的精神和原则之内，所做的修改、等同替换等均应包含在发明的保护范围之内。