前言：
全球海平面变化是气候变化最敏感的指示器之一。传统验潮站和卫星高度计提供了多年的连续记录，但近年来星载激光雷达（LiDAR）——尤其是ICESat-2卫星上的ATLAS（先进地形激光测高系统）——凭借其光子计数技术，能够以厘米级精度获取海面及地物的三维信息。然而，原始光子数据中混杂着大量大气散射、太阳背景噪声以及陆地回波，如何从中准确分离出属于海面的光子群，是提取海平面高度的关键步骤。
因此，本文分享一个简洁高效的海平面提取算法，并完整展示从数据读取、正态拟合到可视化输出的全过程。
希望各位同学点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
1
数据
ICESat-2于2018年9月发射，搭载的ATLAS系统采用多波束微脉冲光子计数技术，以10 kHz的重复频率向地面发射6束激光（分为强/弱波束，间隔3.3 km）。其第三代数据产品ATL03提供了所有光子事件的经纬度、高度、时间及质量标记，是后续ATL12（海洋高度）、ATL06（陆地冰高度）等高级产品的基础。
本案例使用的原始数据来自于ATL03产品的一次过境（轨道号 1777559014556），观测时间为2025年9月11日13时29分（UTC），大致覆盖了某片开阔海洋区域。原数据为HDF5格式，但为了方便教学，已提前提取了关键字段并导出为CSV文本文件，文件名为 photon_2025-09-11_t1329_1777559014556.csv。
2
研究方法
在开阔海面上，噪声光子在高度方向上的分布几乎是均匀的（来自大气散射和暗电流），而海面回波由于镜面反射外加少量水体散射，会形成明显的峰值。因此，我们将所有光子的高度直方图视为一个 高斯分布 + 均匀背景噪声的混合模型。如果只关心海面平均高度，最简单的策略就是找到分布的主峰位置。然而，实际数据可能存在偶尔的稀疏极高/极低点（比如海浪破碎、前向散射的云层光子），它们会拉高方差、偏移均值。为此，我们引入递归截断方法：先拟合一个初始正态分布，然后截掉距均值超过某个倍数标准差（n_sigmas）的离群点，再对剩余点重新拟合，迭代一次或多次，直到分布稳定。代码中实现了两个函数：
get_normal_distribution(dataX, n_sigmas, n)
输入：待拟合的一维数组dataX，截断宽度系数n_sigmas（默认0.5），递归次数n。
过程：
（1）使用scipy.stats.norm.fit估计当前数据的均值mu和标准差sigma。
（2）如果n > 0，则将数据截断到 [mu - n_sigmas*sigma, mu + n_sigmas*sigma] 区间内，然后递归调用自身，并将n减1。
（3）如果n == 0，返回最终的mu, sigma。
经验表明，海面光子的自然分布宽度（sigma_sea）大约对应0.3~0.5 m（视有效波高而定），而噪声光子的标准差则要大得多（可达几十米）。用0.5倍初始标准差截掉大部分噪声后，剩余点几乎全是海面光子，再拟合一次即可得到准确的海面参数。迭代更多次意义不大，因为第二次拟合后分布区间已非常稳定。
get_sea_level(df, index, n_sigmas)
直接调用get_normal_distribution提取光子高度列，返回(mu, (mu - n_sigmas*sigma, mu + n_sigmas*sigma))。
其中sea_range是一个二元组，用于后续提取海面光子。
2.2 海面光子提取
得到海面平均高度mu和标准差sigma后，定义一个区间 [mu - 0.5*sigma, mu + 0.5*sigma]。这一区间的宽度约为1倍标准差（因为0.5+0.5=1），而对于真实海面高斯分布，约68%的光子应落在此区间内。实际计算中，我们设置n_sigmas=0.5，区间总宽度为sigma。代码中的函数：def get_water_surface_points(sea_level, sea_range, df, index): sea_surface_points = df[(df[index] >= sea_range[0]) & (df[index] <= sea_range[1])] return sea_surface_points
该函数返回一个DataFrame，即被判定为“海面”的光子子集，后续可用于单独分析。3
输出与结果
原始数据：
直方图呈现出典型的“尖峰+平底”结构：在mu附近（约25.80米）处出现一个极高峰值，两侧快速下降。红色虚线标出的提取区间将大部分噪声光子排除在外。
结果可视化：
可以看出海面提取效果较好。
4
完整代码
看到最后，如果有用请大家点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
最后，祝大家五一快乐！！！
# -*- coding: utf-8 -*-"""@Time: 2026/4/30 22:33@Author: LXX@File: 海平面提取.py@IDE: PyCharm@Motto: ABC(Always Be Coding)"""import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport matplotlib as mplfrom scipy import stats# ================== Nature 风格绘图配置 ==================mpl.rcParams['font.family'] = 'sans-serif'mpl.rcParams['font.sans-serif'] = ['Arial', 'Helvetica', 'DejaVu Sans']mpl.rcParams['axes.linewidth'] = 1.2mpl.rcParams['axes.labelsize'] = 12mpl.rcParams['xtick.labelsize'] = 10mpl.rcParams['ytick.labelsize'] = 10mpl.rcParams['xtick.major.width'] = 1.2mpl.rcParams['ytick.major.width'] = 1.2mpl.rcParams['xtick.direction'] = 'in'mpl.rcParams['ytick.direction'] = 'in'mpl.rcParams['legend.fontsize'] = 10mpl.rcParams['legend.frameon'] = Falsempl.rcParams['figure.dpi'] = 300# ================== 用户核心函数 ==================def get_normal_distribution(dataX: np.ndarray, n_sigmas: float = 0.5, n: int = 0): """递归截断离群点的高斯拟合""" mu, sigma = stats.norm.fit(dataX) if n > 0: n -= 1 mask = (dataX > mu - n_sigmas * sigma) & (dataX < mu + n_sigmas * sigma) dataX = dataX[mask] return get_normal_distribution(dataX, n_sigmas, n) return mu, sigmadef get_sea_level(df: pd.DataFrame, index: str = "photon height", n_sigmas: float = 0.5): """提取海平面高度及提取区间""" assert index in df.columns, f"索引 {index} 不在数据中" heights = df[index].values mu, sigma = get_normal_distribution(heights, n_sigmas, n=1) return mu, (mu - n_sigmas * sigma, mu + n_sigmas * sigma)def get_water_surface_points(sea_level: float, sea_range: tuple, df: pd.DataFrame, index: str = "photon height") -> pd.DataFrame: """根据区间筛选海面光子""" if sea_level is None or sea_range is None: sea_level, sea_range = get_sea_level(df, index) sea_surface_points = df[(df[index] >= sea_range[0]) & (df[index] <= sea_range[1])] return sea_surface_points# ================== 数据加载 ==================file_path = "photon_2025-09-11_t1329_1777559014556.csv"df = pd.read_csv(file_path)# ================== 海平面提取 ==================mu, sea_range = get_sea_level(df, index="photon height", n_sigmas=0.5)df_water = get_water_surface_points(mu, sea_range, df, index="photon height")print(f"Sea Level (mu): {mu:.4f} m")print(f"Sea Range: {sea_range[0]:.4f} m to {sea_range[1]:.4f} m")print(f"Total points: {len(df)}")print(f"Water surface points: {len(df_water)}")# ================== 可视化部分 ==================# --- 图1: 直方图及高斯拟合 ---fig1, ax1 = plt.subplots(figsize=(6, 4))heights = df["photon height"].valuesp1, p99 = np.percentile(heights, 1), np.percentile(heights, 99)mask = (heights >= p1) & (heights <= p99)ax1.hist(heights[mask], bins=100, density=True, alpha=0.6, color='#7CAE00', label='Photon Height')xmin, xmax = ax1.get_xlim()x = np.linspace(xmin, xmax, 100)sigma_reconstructed = (mu - sea_range[0]) / 0.5p = stats.norm.pdf(x, mu, sigma_reconstructed)ax1.plot(x, p, 'k', linewidth=2, label=f'Fit: $\mu={mu:.2f}$')ax1.axvline(mu, color='#F8766D', linestyle='--', linewidth=2, label='Sea Level')ax1.axvspan(sea_range[0], sea_range[1], color='#F8766D', alpha=0.2, label='Extraction Range')ax1.set_xlabel('Photon Height (m)')ax1.set_ylabel('Density')ax1.legend(loc='upper right')plt.tight_layout()fig1.savefig('fig1_histogram.png')plt.show()plt.close(fig1)# --- 图2: 原始+海面光子叠加散点图 ---fig2, ax2 = plt.subplots(figsize=(8, 4))plot_df = df.sample(50000) if len(df) > 50000 else dfax2.scatter(plot_df['latitude'], plot_df['photon height'], s=1, c='#BEBEBE', alpha=1, label='Raw Photons')ax2.scatter(df_water['latitude'], df_water['photon height'], s=1, c='#00BFC4', alpha=1, label='Water Surface')ax2.axhline(mu, color='#F8766D', linestyle='-', linewidth=1, label='Estimated Sea Level')ax2.axhline(sea_range[0], color='k', linestyle='--', linewidth=1)ax2.axhline(sea_range[1], color='k', linestyle='--', linewidth=1)ax2.set_xlabel('Latitude (°)')ax2.set_ylabel('Photon Height (m)')y_center = muy_window = 30ax2.set_ylim(y_center - y_window, y_center + y_window)ax2.legend(loc='upper right', markerscale=5)plt.tight_layout()fig2.savefig('fig2_overlay.png')plt.show()plt.close(fig2)# --- 图3: 原始光子点云---fig3, ax3 = plt.subplots(figsize=(8, 4))ax3.scatter(plot_df['latitude'], plot_df['photon height'], s=2, c='#C0C0C0', alpha=1, label='Raw Photons', edgecolors='none')ax3.set_xlabel('Latitude (°)')ax3.set_ylabel('Photon Height (m)')ax3.set_ylim(y_center - y_window, y_center + y_window)ax3.legend(loc='upper right', markerscale=4)plt.tight_layout()plt.savefig('fig3_raw_photons.png', dpi=300, bbox_inches='tight')plt.show()print("Images saved.")
10个月宝宝每天需要喝多少奶粉？
