LXX
读完需要
速读仅需 1 分钟
前言：
ICESat-2 数据经过处理后，得到大致的地面点数据，对地面点数据进行过滤拟合得出最终的地面点，一般有三角网加密或是滤波等方法。因此，本文分享一种多次迭代和滤波最终提取出地面点的方法。
希望各位同学点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
1
在遥感和地理信息系统（GIS）分析中，特别是当处理地形数据时，常常需要识别并提取地面点（ground points）以准确描述地形表面。这些地面点的提取对于生成数字高程模型（DEM）、分析水流路径和进行生态环境监测等任务至关重要。然而，遥感数据中经常包含各种噪声和非地面点（如植被和建筑物反射的光子），这些干扰需要被有效地滤除。
此处为一种迭代滤波方法，结合多项式拟合、高斯平滑和百分位数滤波来提取地面点。这种方法可以有效地处理和分析沿轨道（Along-Track）和高程（Height）数据。
主要包括以下方法：
（1）多项式拟合：多项式拟合通过最小二乘法估计多项式系数，使得拟合曲线能够尽可能准确地描述数据的整体趋势。多项式的阶数（如线性拟合、二次拟合等）可以根据数据的复杂程度进行选择。
（2）高斯平滑：高斯平滑是一种信号处理技术，通过将高斯核函数与数据进行卷积，达到平滑信号的效果。高斯核函数的标准差（σ）决定了平滑的程度，较大的σ会导致更平滑的曲线。
（3）残差计算：残差是衡量数据点偏离拟合曲线程度的指标，通过计算残差可以识别出拟合效果较差的数据点（可能是噪声或异常值）。
（4）百分位数滤波：百分位数滤波通过计算数据分布的百分位数，确定一个合理的范围，保留主要数据点，滤除极端值。这种方法对处理非正态分布的数据尤其有效。
（5）迭代处理：迭代处理是一种逐步逼近的过程，每次迭代都会在前一次的基础上进一步优化结果。通过多次迭代，可以逐步去除噪声和异常值，最终提取出准确的地面点。
2
（1）数据读取：读取包含沿轨道和高程数据的 CSV 文件，转换为 Numpy 数组，以便进行后续处理。
import warningsimport pandas as pdimport numpy as npimport matplotlib.pyplot as pltfrom scipy.ndimage import gaussian_filterfrom xy.class.polyfit_functions import polyfit_numpydata = pd.read_csv('ground_points.csv')xy = data[['Along-Track (m)', 'Height (m MSL)']].to_numpy()xy_loop = np.array(xy, copy=True)
（2）多项式拟合：使用多项式拟合方法对数据进行拟合，以捕捉数据的整体趋势。多项式拟合是通过最小二乘法来估计多项式系数，使得多项式曲线尽可能接近数据点。
（3）高斯平滑：对拟合结果进行高斯平滑，以减少噪声和局部波动。高斯平滑通过卷积操作，将高斯核函数应用于数据，从而实现平滑效果。
for i in range(20):print(f'Loop {i+1}')print(f'Number of photons: {len(xy_loop)}')coeffs, pfit_elev = polyfit_numpy(xy_loop, 150, 1)pfit_elev_smooth = gaussian_filter(pfit_elev, sigma=5)pfit_elev_smooth[np.abs(xy_loop[:,1]-pfit_elev_smooth) > 150] = np.nanprint(f'Number of elements in smoothed polyfit: {np.sum(~np.isnan(pfit_elev_smooth))}')
（4）残差计算：计算数据点与平滑后的拟合曲线之间的残差，以评估拟合的质量。残差是数据点的实际值与拟合值之间的差异。
with warnings.catch_warnings():warnings.simplefilter("ignore", category=RuntimeWarning)residuals = xy_loop[:,1] - pfit_elev_smoothprint(f'Max residual at index {np.nanargmax(residuals)} with {np.nanmax(residuals):.2f} m')print(f'Min residual at index {np.nanargmin(residuals)} with {np.nanmin(residuals):.2f} m')print(f'StDev of residuals: {np.nanstd(residuals):.2f}')
（5）百分位数滤波：按 30 米区间对数据进行分组，并计算每组残差的百分位数（2%-98%）。通过滤除残差超出此范围的数据点，保留主要数据。这一步有助于进一步减少噪声和异常值。
（6）迭代处理：多次迭代上述步骤，每次迭代后更新数据，逐步去除非地面点和噪声，保留地面点。每次迭代都会过滤掉一些离群点，使得剩余的数据更加接近真实地面点。
split_at = xy_loop[:, 0].searchsorted(np.arange(0, int(max(xy_loop[:,0])), 30))xy_bins = np.split(xy_loop, split_at)residuals_bins = np.split(residuals, split_at)new_xy_bins = []for j in range(len(xy_bins)):filtered_ind, = np.where((residuals_bins[j] > np.nanpercentile(residuals_bins[j], 2)) &(residuals_bins[j] < np.nanpercentile(residuals_bins[j], 98)))new_xy_bins.append(xy_bins[j][filtered_ind])xy_loop = np.concatenate(new_xy_bins)np.save(f'xy_loop_{i}.npy', xy_loop)print(f'No of ground photon candidates: {xy_loop.shape[0]}\n')
（7）结果保存和可视化：将提取的地面点保存为新的 CSV 文件，并通过可视化展示原始数据和地面点，以便于检查和验证结果。
ground_points = pd.DataFrame({'Along-Track (m)': xy_loop[:, 0], 'Height (m MSL)': xy_loop[:, 1]})ground_points.to_csv('ground_points.csv', index=False)plt.scatter(xy[:, 0], xy[:, 1], color='blue', label='All Photons')plt.scatter(xy_loop[:, 0], xy_loop[:, 1], color='red', label='Ground Photons')plt.xlabel('Along-Track (m)')plt.ylabel('Height (m MSL)')plt.legend()plt.show()
3
在地形起伏较大区域效果不是太好，仍需后续改进......
