LXX
读完需要
速读仅需 1 分钟
前言:
好久没更新 ICESat-2 的内容了,最近忙里偷闲分享一个 ICESat-2 的 ATL03 数据去噪的方法,是对文献《基于统计直方图两步法的星载单光子数据去噪》中部分精去噪方法的复现。
希望各位同学点个关注,点个小赞,这将是更新的动力,不胜感激❥(^_-)
1
由于 ICESat-2/ATLAS 原始光子点云条带跨度长,研究区范围较大,且在光子点云沿轨距离和高程的二维剖面中,存在着大量的噪声点和信号点,如果直接对光子点云进行去噪,势必会增加去噪多余的工作和运算时间,效率不高且受到大量噪声点的影响使得去噪结果不佳。因此在进行精去前,基于信号光子点和噪声光子点在高程方向上密度存在差异的特性,采用基于高程统计直方图的粗去噪方法,旨在确定原始光子点云中信号光子点的高程范围,为后续的去噪工作降低数据量。
方法主要是对光子点云条带进行格网划分,每个格网窗口进行光子点云数据的高程统计,并建立每个窗口内的光子高程统计直方图,然后利用高斯滤波来找到每个格网中高程统计直方图的峰值位置,由于信号光子点云在空间中更加聚集且在统计中有着更大的频数,因此将高程统计直方图的峰值确定为信号光子的中心位置。
具体方法见前期内容
2
粗去噪后,光子数据中仅含少量噪声,剩余有效信号光子之间距离小于有效信号光子与噪声光子或噪声光子与噪声光子之间的距离。为显著区分噪声光子,本文基于式(5)计算每个光子的 K 邻域距离平方的均值作为衡量标准。
将统计数据集 D 内所有光子的 K 邻域距离平方的均值,按其范围划分为 dn 个间隔,统计每个间隔内的光子数量,绘制距离平方频数直方图,该直方图近似呈单侧正态分布,故同样满足正态分布特性。以 1*σ 置信度作为信号与噪声的分界点,对应有效信号光子置信区间如式(6)所示:
3
本文给出大家完整代码,如果对您有用的话还请点个小赞,多谢呀!
# -*- coding: utf-8 -*-
"""
@Time: 2024/3/11 11:13
@Author: LXX
@File: 距离平方统计直方图去噪.py
@IDE:PyCharm
@Motto:ABC(Always Be Coding)
"""
import tkinter as tk
from tkinter import filedialog
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import KDTree
# 打开文件选择弹窗并选择文件
def open_file_dialog():
root = tk.Tk()
root.withdraw() # 隐藏主窗口
file_path = filedialog.askopenfilename(title="Select ICESat-2 CSV File", filetypes=[("CSV files", "*.csv")])
return file_path
# 计算K邻域内的平均距离平方
def calculate_dist_ave(along_track, height, K=5):
dist_ave = np.zeros(len(along_track))
for i in range(len(along_track)):
neighbors_idx = np.arange(max(0, i - K // 2), min(len(along_track), i + K // 2 + 1))
neighbors = np.sqrt(
(along_track[i] - along_track[neighbors_idx]) ** 2 + (height[i] - height[neighbors_idx]) ** 2)
dist_ave[i] = np.mean(neighbors)
return dist_ave
# 去噪处理函数
def remove_noise(df, avg_distances, threshold=1):
# 根据平均距离和阈值判断噪声
signal_mask = avg_distances < threshold
df_valid = df[signal_mask] # 有效信号
df_noise = df[~signal_mask] # 噪声点
valid_points = df_valid.index # 有效点的索引
return df_valid, df_noise, avg_distances, valid_points
# 去噪处理
def denoise_data(df, along_track, height, threshold=2, K=5):
avg_distances = calculate_dist_ave(along_track, height, K)
df_valid, df_noise, avg_distances, valid_points = remove_noise(df, avg_distances, threshold)
return df_valid, df_noise, avg_distances
# 主函数
def main():
# 打开文件
file_path = open_file_dialog()
# 读取CSV文件
data = pd.read_csv(file_path)
along_track = data['Along-Track (m)'].values
height = data['Height (m HAE)'].values
# 去噪处理
df_valid, df_noise, avg_distances = denoise_data(data, along_track, height)
# 画图
plt.figure(figsize=(12, 8))
# 原始数据图
plt.subplot(3, 1, 1)
plt.scatter(along_track, height, s=0.01, c='blue', label='Original Data')
plt.xlabel('Along-Track (m)')
plt.ylabel('Height (m HAE)')
plt.title('Original ICESat-2 Data')
plt.legend()
# 平均距离图 - 使用颜色条来表示 Average Distance
plt.subplot(3, 1, 2)
# 使用scatter, c=avg_distances 表示用颜色条映射平均距离的值,cmap是颜色映射
sc = plt.scatter(along_track, avg_distances, c=avg_distances, cmap='magma', s=0.1, label='Average Distance')
# 添加颜色条,表示平均距离的数值
cbar = plt.colorbar(sc)
cbar.set_label('Average Distance')
# 坐标轴和标题
plt.xlabel('Along-Track (m)')
plt.ylabel('Average Distance')
plt.title('Average Distance of K Neighbors')
# 添加图例
plt.legend()
'''
# 平均距离图
plt.subplot(3, 1, 2)
plt.plot(along_track, avg_distances, c='orange', label='Average Distance (K Neighbors)')
plt.xlabel('Along-Track (m)')
plt.ylabel('Average Distance')
plt.title('Average Distance of K Neighbors')
plt.legend()
'''
# 去噪后数据图
plt.subplot(3, 1, 3)
plt.scatter(df_valid['Along-Track (m)'], df_valid['Height (m HAE)'], s=1, c='red', label='Denoised Data')
plt.xlabel('Along-Track (m)')
plt.ylabel('Height (m HAE)')
plt.title('Denoised ICESat-2 Data')
plt.legend()
# 显示图形
plt.tight_layout()
plt.show()
if __name__ == '__main__':
main()
4
可以看到去噪结果尚可,但存在一定的离群光子点云尚未去除,可采用离群点检测等方法进行去除;同时该方法也需要手动设置参数,自适应阈值供大家探索开发。
下面是几个去噪代码,创作不易大家有需要再看。
改进DBSCAN进行ICESat-2光子点云去噪(完整代码)
【参考】
[1]焦慧慧,谢俊峰,刘仁,等.基于统计直方图两步法的星载单光子数据去噪[J].太赫兹科学与电子信息学报, 2023, 21(3):384-391.