LXX
读完需要
速读仅需 1 分钟
前言：
ICESat-2 数据经过几次去噪，绝大多数的噪声点已经被剔除，但由于仪器和大气环境的影响，光子点云的噪声可能也会出现聚集的情况，而通过上述两级去噪会把这类情况也误判为信号光子点，这对于后续的分析影响是较大的。本文分享一种箱线图去噪方法，是对某论文的方法复现。
希望各位同学点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
1
分位数/箱线图作为数据挖掘中常用的方法之一，可以很好地检测到数据中的异常值和离群点。其对于二次去噪后对光子点云中可能残留的噪声光子点的剔除是适用的。因此通过箱线图的方法对光子点云数据进行再次去噪，其中分为地形相对平坦状况下的分位数去噪和地形起伏状况下的分位数去噪。具体步骤如下：
（1）当光子点云所处地形相对平坦时，可以直接使用分位数方法。首先按以下公式 3-8 计算光子高程方向上的上四分位数（Q=75）和下四分位数（Q=25），同时计算两者的间距 QR：
然后，通过计算箱线图的上下最值确定范围，以确定去噪的阈值。根据阈值，将不在阈值范围内的噪声点剔除，计算方法按公式示：
式中：Point（j）为第 j 个光子点的判断，为第 j 个光子点的高程信息。noise 和 signal 分别为噪声点和信号点。基于此便可以去除平坦地形状况下的可能残留的噪声光子点。
（2）当光子点云所处的地形有一定的坡度时，运用上述的分位数方法便会出现无法监测到噪声点的情况。这是因为在平坦地区时存留在信号光子点上下的噪声点比较明显，较容易去除；而在地形有起伏的区域，在沿轨方向上噪声点的分布是非常离散的，使得箱线图计算的上下四分位数的间距很大，此时运用分位数无法取得很好的效果。下图显示了运用箱线图去噪的方法，光子点处于地形起伏区域，噪声光子点的分布较为离散，此时的箱线图无法很好地检测出噪声光子点。
但是，从图可以看出，光子点云在沿坡度方向更加聚集，类似于水平方向的光子点云聚集形态。因此可以将有坡度的光子点云进行旋转，使其处于水平方向，然后再用箱线图方法进行去噪。光子点云的旋转按照如下公式所示：
式中，h2 为光子点云经过旋转 x 度所得到的高程值。Std 用来计算标准差。al_t 为光子沿轨方向的距离。上式是通过计算旋转后光子点云高程的标准差是否达到最小，若达到，则经过该旋转角度转换后的光子点云沿水平方向分布。下图为旋转点云进行箱线图去噪的方法，左图为旋转后呈水平方向的光子点云，右图为检测到的噪声点箱线图。因此通过旋转点云再进行去噪的这种方法是可行的。
2
1.导入库
import osimport pandas as pdimport numpy as npimport matplotlib.pyplot as pltimport tkinter as tkfrom tkinter import filedialogimport seaborn as snsimport time
2.加载数据
# 读取CSV数据文件root = tk.Tk()root.withdraw()file_path = filedialog.askopenfilename(filetypes=[('CSV Files', '*.csv')])# 读取CSV文件data = pd.read_csv(file_path, header=0)# 获取原文件路径、文件名和后缀file_dir, file_name = os.path.split(file_path)file_prefix, file_suffix = os.path.splitext(file_name)
3.将光子数据进行旋转至水平方向
def rotate_to_horizontal(segment):angles = []x = segment['Along-Track (m)'].valuesy = segment['Height (m HAE)'].valuesother_cols = {col: segment[col].values for col in segment.columns if col not in ['Along-Track (m)', 'Height (m HAE)']}for phi in range(0, 180, 2):x_prime = x * np.sin(np.deg2rad(phi)) + y * np.cos(np.deg2rad(phi))std = np.std(x_prime)angles.append((phi, std))rotation_angle = min(angles, key=lambda x: x[1])[0]x_prime = x * np.sin(np.deg2rad(rotation_angle)) + y * np.cos(np.deg2rad(rotation_angle))y_prime = -x * np.cos(np.deg2rad(rotation_angle)) + y * np.sin(np.deg2rad(rotation_angle))rotated_segment = pd.DataFrame({'Along-Track (m)': x_prime, 'Height (m HAE)': y_prime, **other_cols})return rotated_segment, rotation_angle
4.计算光子点云高程Height方向的上下四分位数,以及上下四分位间距IQR
def calculate_iqr(segment):height = segment['Height (m HAE)'].valuesq1 = np.percentile(height, 25)q3 = np.percentile(height, 75)iqr = q3 - q1return q1, q3, iqr
5.计算箱线图上下限去噪阈值，将阈值之外的点作为噪声点将其剔除
def denoise_segment(segment, q1, q3, iqr):......return segment.dropna()
6.将去噪结果旋转至原方向
def rotate_to_original(segment, rotation_angle):x_prime = segment['Along-Track (m)'].valuesy_prime = segment['Height (m HAE)'].valuesother_cols = {col: segment[col].values for col in segment.columns if col not in ['Along-Track (m)', 'Height (m HAE)']}x = x_prime * np.sin(np.deg2rad(rotation_angle)) - y_prime * np.cos(np.deg2rad(rotation_angle))y = x_prime * np.cos(np.deg2rad(rotation_angle)) + y_prime * np.sin(np.deg2rad(rotation_angle))original_segment = pd.DataFrame({'Along-Track (m)': x, 'Height (m HAE)': y, **other_cols})return original_segment
.......
3
（创作不易，完整代码后台发送“箱线图去噪”）
