简介这份资源面向卫星遥感、海洋学与地球物理领域的科研人员及技术开发者围绕国产风云三号E星FY-3E搭载的GNOS-II仪器复现基于星载GNSS-R的海面高度反演模型。内容融合传统物理模型与机器学习方法采用随机森林和卷积神经网络对比评估北斗与GPS反射信号的反演性能并给出数据预处理、特征提取、模型训练与误差分析等完整实现思路。资源包为1个PDF文件约889KB集中呈现论文复现所需的代码与解释便于读者对照理解物理模型公式推导与机器学习建模流程。目前已有124人学习适合希望掌握国产卫星测高技术、开展GNSS-R海面测高研究或搭建反演实验框架的读者参考可从中获取从理论到代码落地的系统指导。1. FY-3E 星载 GNSS-R 海面高度反演从国产卫星数据到可复现的建模路径FY-3E 是国产极轨气象卫星里第一颗搭载 GNSS-R 载荷的业务星它接收导航卫星直射信号和海面反射信号通过两者的路径差反演海面高度、风速、海冰等参数。海面高度反演这件事本质是把一段时延-多普勒图DDM里的镜面反射点位置、波形前沿、相干/非相干特性映射成一个厘米到分米级的测高结果。传统做法靠几何模型加波形重跟踪物理可解释但误差项多纯机器学习做法拟合能力强却容易在训练集外的海况下翻车。这篇要讲的就是把物理模型和机器学习融合起来在 FY-3E 数据上跑通一套海面高度反演模型从数据读取、特征构造、模型搭建到误差评估每一步都给出可抄的代码和参数解释。适合做遥感反演、GNSS-R 信号处理、以及想把机器学习落地到国产卫星数据的从业者。2. FY-3E GNSS-R 数据怎么读从 L1 级 DDM 到反演可用的观测量2.1 先搞清楚 FY-3E GNSS-R 的数据层级和文件结构FY-3E 的 GNSS-R 产品一般分 L1 和 L2 两级。L1 是经过定标和几何定位的时延-多普勒图附带镜面反射点的经纬度、入射角、方位角、接收机高度等辅助信息L2 是已经反演好的海面风速、海面高度等地球物理参数。做反演模型设计通常从 L1 入手因为你要自己控制特征和标签的构造过程。文件格式常见的是 HDF5内部按轨道分段存储每个段包含若干积分周期通常 1ms 或 1s 量级的 DDM 矩阵。读数据第一步不是急着建模而是确认三件事DDM 的时延轴和多普勒轴怎么定义、镜面反射点坐标怎么算出来的、每个 DDM 对应的接收机高度和入射角是否已经做过几何改正。这三件事决定了你后面构造的特征有没有物理意义。我一般会先写一个探查脚本把单轨数据里所有 DDM 的维度、时延分辨率、多普勒分辨率、辅助字段名打印出来确认没有字段缺失再往下走。import h5py import numpy as np # 打开单轨 FY-3E GNSS-R L1 文件 f h5py.File(FY3E_GNOS_L1_20230101.h5, r) # 打印顶层组结构确认数据组织方式 def print_structure(name, obj): if isinstance(obj, h5py.Dataset): print(fDataset: {name}, shape{obj.shape}, dtype{obj.dtype}) elif isinstance(obj, h5py.Group): print(fGroup: {name}) f.visititems(print_structure) # 读取一个积分周期的 DDM 和辅助几何参数 ddm f[/Science/DDM][0] # 形状通常为 (delay_bins, doppler_bins) sp_lat f[/Science/SP_Lat][0] # 镜面反射点纬度 sp_lon f[/Science/SP_Lon][0] # 镜面反射点经度 inc_angle f[/Science/IncidenceAngle][0] # 入射角 rx_height f[/Science/RxHeight][0] # 接收机高度 print(ddm.shape, sp_lat, sp_lon, inc_angle, rx_height)这段代码先遍历文件结构再取第一个积分周期的 DDM 和几何参数。参数说明delay_bins是时延维通常对应 0.25 个码片或更细doppler_bins是多普勒维分辨率取决于相干积分时间。IncidenceAngle和RxHeight是后面几何改正和特征构造的关键输入如果这两个字段有缺失整条反演链路都要重新考虑。2.2 从 DDM 里提取海面高度相关的观测量DDM 本身是一个二维功率图直接扔给机器学习模型不是不行但维度高、冗余大而且物理意义被淹没。常见做法是先做波形重跟踪提取镜面反射点的时延偏移再结合几何关系换算成高度。具体来说海面高度反演的核心观测量是镜面反射点相对于参考椭球面的高度它由接收机高度、入射角、以及反射信号相对于直射信号的额外路径延迟共同决定。我一般会从 DDM 里提取三类特征第一类是波形前沿特征包括前沿 10%、50%、90% 功率对应的时延位置以及前沿斜率第二类是峰值特征包括峰值功率、峰值时延、峰值多普勒第三类是形状特征包括 DDM 的等效宽度、对称性、以及相干分量占比。这三类特征加起来大概 15 到 20 维既保留了物理可解释性又不会让模型过拟合。def extract_waveform_features(ddm, delay_axis): 从单个 DDM 提取波形前沿和峰值特征 # 沿多普勒维取最大值得到一维时延波形 waveform np.max(ddm, axis1) waveform_norm waveform / (np.max(waveform) 1e-12) # 峰值位置和峰值功率 peak_idx np.argmax(waveform) peak_power waveform[peak_idx] peak_delay delay_axis[peak_idx] # 前沿 10%、50%、90% 功率对应的时延 def crossing_level(level): idx np.where(waveform_norm level)[0] return delay_axis[idx[0]] if len(idx) 0 else np.nan delay_10 crossing_level(0.1) delay_50 crossing_level(0.5) delay_90 crossing_level(0.9) # 前沿斜率用 10% 到 90% 的时延差近似 leading_edge_slope (0.9 - 0.1) / (delay_90 - delay_10 1e-12) # 等效宽度用波形二阶矩 power_sum np.sum(waveform_norm) mean_delay np.sum(delay_axis * waveform_norm) / (power_sum 1e-12) equiv_width np.sqrt( np.sum((delay_axis - mean_delay) ** 2 * waveform_norm) / (power_sum 1e-12) ) return { peak_power: peak_power, peak_delay: peak_delay, delay_10: delay_10, delay_50: delay_50, delay_90: delay_90, leading_edge_slope: leading_edge_slope, equiv_width: equiv_width }这段函数把二维 DDM 压成一维波形再提取七个标量特征。参数说明delay_axis是时延轴的实际物理值单位通常是码片或米必须和 DDM 的列索引对应上否则算出来的时延偏移没有物理意义。crossing_level里用idx[0]取第一个超过阈值的点这是前沿检测的常用做法但如果波形有噪声毛刺建议先做一次 3 点滑动平均再检测。2.3 几何改正把时延偏移换算成海面高度有了峰值时延和前沿时延下一步是几何改正。GNSS-R 测高的基本几何关系是反射路径比直射路径多走的距离等于反射信号相对直射信号的额外时延乘以光速。这个额外路径延迟再投影到天底方向就得到镜面反射点相对于参考面的高度。实际处理中还要考虑地球曲率、大气延迟、以及接收机姿态误差。我一般会先算一个几何高度初值再用经验模型做大气和潮汐改正。几何高度初值的公式是海面高度 接收机高度 - (额外路径延迟 / 2) * cos(入射角)。这里额外路径延迟取峰值时延相对于直射信号时延的偏移量。注意这个公式假设镜面反射点正好在接收机和导航卫星的连线与地球表面的交点附近实际处理中需要用迭代法修正镜面反射点位置。def geometric_height_correction(peak_delay, direct_delay, inc_angle, rx_height): 几何改正从时延偏移换算海面高度初值 c 299792458.0 # 光速单位 m/s # 额外路径延迟单位米 extra_path (peak_delay - direct_delay) * c # 投影到天底方向 height_geom rx_height - (extra_path / 2.0) * np.cos(np.deg2rad(inc_angle)) return height_geom # 假设 direct_delay 已经从辅助数据里读到 height_geom geometric_height_correction( peak_delaypeak_delay, direct_delaydirect_delay, inc_angleinc_angle, rx_heightrx_height ) print(f几何高度初值: {height_geom:.3f} m)这段代码给出几何高度初值。参数说明direct_delay是直射信号的时延参考通常从辅助数据里读取如果没有可以用 DDM 的时延轴零点近似。inc_angle单位是度rx_height单位是米。算出来的height_geom只是初值后面还要用机器学习模型做残差修正把大气延迟、海况偏差、以及几何近似误差一起吸收掉。3. 物理模型和机器学习怎么融合残差建模与特征工程3.1 为什么纯物理模型不够纯机器学习也不稳纯物理模型的问题在于它把海面当成一个理想镜面但实际海面有波浪、有泡沫、有风生粗糙度这些因素会让反射信号的前沿展宽、峰值偏移物理模型很难把所有误差项都写清楚。纯机器学习的问题在于它不知道物理规律如果训练集里缺少高海况样本模型在高风速下会给出离谱的高度值。我见过最典型的翻车场景是模型在低风速下误差只有几厘米一到高风速就跳到几十厘米因为训练集里高风速样本太少模型把噪声当成了信号。融合的思路是物理模型负责给出一个可解释的初值机器学习负责学习残差。残差 真实高度 - 物理模型高度。这样模型只需要学习物理模型没覆盖的那部分误差学习难度大大降低而且外推能力比纯机器学习好。常见做法有两种一种是直接把物理模型输出作为特征之一和 DDM 特征一起送进模型另一种是先算残差再单独训练一个残差模型。我一般用第二种因为残差的分布更集中模型更容易收敛。3.2 特征工程把物理量、几何量和 DDM 特征拼成一张表特征表的设计决定了模型的上限。我一般会把特征分成四组第一组是 DDM 波形特征就是上一章提取的那七个第二组是几何特征包括入射角、方位角、接收机高度、镜面反射点经纬度第三组是物理模型输出包括几何高度初值、大气改正量、潮汐改正量第四组是辅助特征包括导航卫星高度角、信号载噪比、以及海况指示量比如风速初值。这四组特征加起来大概 25 到 30 维。注意经纬度不要直接送进模型因为经纬度是周期性变量直接送进去会让模型误以为经度 179 和 -179 差很远。常见做法是把经纬度转成 sin/cos 编码或者直接去掉经纬度改用海区标识。我一般会保留 sin/cos 编码因为不同海区的海况差异确实会影响反演精度。import pandas as pd import numpy as np def build_feature_table(ddm_features, geom_features, phys_features, aux_features): 把四组特征拼成一张宽表 df pd.DataFrame() # 第一组DDM 波形特征 for k, v in ddm_features.items(): df[k] v # 第二组几何特征经纬度做 sin/cos 编码 df[inc_angle] geom_features[inc_angle] df[azimuth] geom_features[azimuth] df[rx_height] geom_features[rx_height] df[sp_lat_sin] np.sin(np.deg2rad(geom_features[sp_lat])) df[sp_lat_cos] np.cos(np.deg2rad(geom_features[sp_lat])) df[sp_lon_sin] np.sin(np.deg2rad(geom_features[sp_lon])) df[sp_lon_cos] np.cos(np.deg2rad(geom_features[sp_lon])) # 第三组物理模型输出 df[height_geom] phys_features[height_geom] df[atmo_corr] phys_features[atmo_corr] df[tide_corr] phys_features[tide_corr] # 第四组辅助特征 df[sat_elev] aux_features[sat_elev] df[cn0] aux_features[cn0] df[wind_speed_init] aux_features[wind_speed_init] return df这段代码把四组特征拼成一张表。参数说明sp_lat和sp_lon是镜面反射点经纬度做 sin/cos 编码后每个变成两维避免周期性断裂。height_geom是上一章算的几何高度初值atmo_corr和tide_corr是大气和潮汐改正量如果没有现成数据可以先用经验模型算一个近似值。wind_speed_init是海面风速初值可以从 L2 产品里读也可以从 DDM 特征里估。3.3 残差建模用梯度提升树还是神经网络残差建模的模型选择我一般先试梯度提升树比如 XGBoost 或 LightGBM再试一个浅层神经网络做对比。梯度提升树的优势是训练快、对特征尺度不敏感、可解释性好能看特征重要性神经网络的优势是能捕捉特征之间的非线性交互但需要更多调参和正则化。在 FY-3E 这种样本量不算特别大的场景下梯度提升树通常是更稳的选择。训练残差模型时标签是真实高度减去物理模型高度。真实高度从哪来常见做法是用验潮站数据、或者用其他测高卫星比如 Jason 系列的交叉点做参考。如果没有外部参考也可以用 L2 产品里的海面高度作为标签但要注意 L2 产品本身也有误差这时候你训练出来的模型是在拟合 L2 的误差而不是真实误差。我一般会优先用验潮站数据因为它的精度最高但空间覆盖有限需要做时空匹配。import lightgbm as lgb from sklearn.model_selection import train_test_split from sklearn.metrics import mean_absolute_error, root_mean_squared_error # 假设 df 是特征表label 是残差标签 X df.drop(columns[residual]) y df[residual] X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, random_state42 ) # LightGBM 参数树数量、学习率、叶子数、最小叶子样本数 params { objective: regression, metric: mae, n_estimators: 800, learning_rate: 0.03, num_leaves: 63, min_child_samples: 20, subsample: 0.8, colsample_bytree: 0.8, reg_alpha: 0.1, reg_lambda: 0.1, random_state: 42 } model lgb.LGBMRegressor(**params) model.fit( X_train, y_train, eval_set[(X_test, y_test)], eval_metricmae, callbacks[lgb.early_stopping(50), lgb.log_evaluation(100)] ) y_pred model.predict(X_test) print(fMAE: {mean_absolute_error(y_test, y_pred):.4f} m) print(fRMSE: {root_mean_squared_error(y_test, y_pred):.4f} m)这段代码训练一个 LightGBM 残差模型。参数说明n_estimators是树的数量配合early_stopping用实际训练时可能不到 800 就停了learning_rate是学习率0.03 是比较稳的值太大容易震荡太小训练慢num_leaves是叶子数63 是中等复杂度如果过拟合就降到 31min_child_samples是最小叶子样本数20 是防止过拟合的常用值。subsample和colsample_bytree是行采样和列采样比例0.8 是经验值。reg_alpha和reg_lambda是 L1 和 L2 正则化系数0.1 是轻量正则。3.4 模型评估不能只看 MAE还要看海况分层误差评估残差模型时只看整体 MAE 和 RMSE 是不够的因为海面高度反演的误差和海况强相关。我一般会按风速分层看低风速5 m/s、中风速5-10 m/s、高风速10 m/s三档的误差分别是多少。如果高风速档的误差明显大于低风速档说明模型在高海况下外推能力不足需要补充高风速样本或者调整特征。另一个要看的指标是误差的空间分布。把测试集的误差按镜面反射点经纬度画出来看有没有明显的区域聚集。如果某个海区的误差特别大可能是该海区的海况特殊或者训练集里该海区样本太少。这时候可以考虑按海区分组做交叉验证或者给不同海区加样本权重。# 按风速分层评估 df_test X_test.copy() df_test[y_true] y_test df_test[y_pred] y_pred df_test[abs_error] np.abs(df_test[y_true] - df_test[y_pred]) bins [0, 5, 10, 100] labels [low_wind, mid_wind, high_wind] df_test[wind_bin] pd.cut(df_test[wind_speed_init], binsbins, labelslabels) for label in labels: subset df_test[df_test[wind_bin] label] if len(subset) 0: print(f{label}: MAE{subset[abs_error].mean():.4f} m, f样本数{len(subset)})这段代码按风速分层算 MAE。参数说明wind_speed_init是特征表里的风速初值分档阈值 5 m/s 和 10 m/s 是 GNSS-R 海面高度反演里常用的分界点。如果高风速档样本数太少算出来的 MAE 统计意义不强需要结合置信区间一起看。4. 避坑与排查FY-3E GNSS-R 反演里最容易翻车的五个地方4.1 现象模型在测试集上 MAE 很小但实际应用时误差突然变大原因训练集和测试集来自同一轨或相邻轨数据分布太接近模型没有见过不同海况或不同几何构型的样本。这是典型的过拟合到轨道特征而不是学到物理规律。解决按轨道做交叉验证而不是随机划分。具体做法是把数据按轨道号分组每次留一整轨做测试其余轨做训练。这样能真实反映模型在新轨道上的表现。如果按轨道交叉验证的误差明显大于随机划分说明模型泛化能力不足需要增加训练轨数量或者简化模型。4.2 现象DDM 波形特征里出现大量 NaN模型训练时报错原因某些积分周期的 DDM 功率太低波形前沿检测找不到超过阈值的点crossing_level返回 NaN。这种情况在低入射角或高海况下比较常见。解决在特征提取阶段加一个有效性判断如果峰值功率低于某个阈值比如噪声底的 3 倍直接标记该样本无效不送进训练。对于部分特征缺失的样本可以用中位数填充但要加一个缺失指示特征让模型知道这个值是被填充的。4.3 现象几何高度初值和真实高度偏差很大残差模型学不动原因几何改正公式里的direct_delay取错了或者入射角没有做地球曲率改正。GNSS-R 的几何关系在低入射角下对误差很敏感入射角差 1 度高度可能差几十厘米。解决先单独验证几何改正链路拿几个已知高度的验潮站数据做对比看几何高度初值的偏差是否在合理范围内。如果偏差大检查direct_delay的符号和单位检查入射角是否已经做过地球曲率改正。常见做法是用迭代法重新算镜面反射点位置直到收敛。4.4 现象残差模型的特征重要性里经纬度编码占了主导原因经纬度 sin/cos 编码虽然解决了周期性问题但如果训练集里某些海区样本特别多模型会倾向于用经纬度来记忆海区而不是学物理规律。这会导致模型在新海区表现很差。解决要么去掉经纬度特征改用海区标识要么给不同海区的样本加权重让每个海区的总权重相近。我一般会先试去掉经纬度看误差有没有明显变化。如果去掉后误差变化不大说明经纬度特征本来就是冗余的如果误差明显变大说明海区差异确实重要这时候应该保留但加权重。4.5 现象模型训练时损失下降很快但验证集损失早早停止下降原因模型复杂度太高或者特征里有泄漏。特征泄漏的常见来源是把真实高度相关的量不小心放进了特征表比如 L2 产品里的海面高度。另一个来源是时间泄漏比如用未来时刻的数据预测过去时刻的高度。解决检查特征表里有没有和标签直接相关的量。所有从 L2 产品里读的字段都要仔细审查确认它不是标签的另一种表达。时间泄漏的检查方法是按时间顺序划分训练集和测试集而不是随机划分。如果按时间划分的误差明显大于随机划分说明有时间泄漏。5. 把模型跑稳的进阶技巧从单轨验证到多轨业务化5.1 用交叉点做独立验证不依赖训练集里的标签残差模型的标签来自验潮站或 L2 产品但验潮站空间覆盖有限L2 产品本身有误差。更独立的验证方法是用交叉点同一颗卫星在不同轨道上经过同一海区时海面高度应该一致。把交叉点上的两个反演结果做差差值就是反演误差的一个估计。这个方法不需要外部参考数据适合做业务化前的自检。具体做法是先找出所有轨道交叉点然后对每个交叉点取两条轨道上最近时刻的反演结果算差值。如果差值的标准差在几厘米量级说明模型一致性不错如果超过 10 厘米说明模型在不同几何构型下表现不稳定需要检查特征里有没有和轨道相关的偏差。def cross_over_validation(df, time_col, lat_col, lon_col, height_col, max_time_diff300): 交叉点验证找时空接近的样本对算高度差 df df.sort_values(time_col).reset_index(dropTrue) diffs [] for i in range(len(df)): for j in range(i 1, len(df)): dt abs(df.loc[j, time_col] - df.loc[i, time_col]) if dt max_time_diff: break dlat abs(df.loc[j, lat_col] - df.loc[i, lat_col]) dlon abs(df.loc[j, lon_col] - df.loc[i, lon_col]) if dlat 0.1 and dlon 0.1: diffs.append(df.loc[j, height_col] - df.loc[i, height_col]) diffs np.array(diffs) print(f交叉点数量: {len(diffs)}) print(f高度差均值: {diffs.mean():.4f} m) print(f高度差标准差: {diffs.std():.4f} m) return diffs这段代码做交叉点验证。参数说明max_time_diff是时间窗口单位秒300 秒是常用值dlat和dlon是空间窗口0.1 度大约是 11 公里。注意这个双重循环在样本量大时会很慢实际使用时可以用 KD 树做空间索引加速。5.2 模型更新策略什么时候该重新训练GNSS-R 反演模型不是训练一次就能一直用的。卫星载荷状态会变海况有季节性变化训练集里的样本分布会逐渐偏移。我一般会监控两个指标一是交叉点验证的高度差标准差二是残差模型在最近一个月数据上的 MAE。如果这两个指标连续两个月上升超过 20%就该考虑重新训练了。重新训练时不要直接把新数据加进去就完事。要先检查新数据的特征分布和旧数据有没有明显差异如果有说明数据生成过程变了需要重新做特征工程。另外新模型上线前一定要用交叉点验证做一次独立评估确认不比旧模型差再切换。5.3 一个具体技巧用分位数回归给出误差范围海面高度反演不光要给一个值还要给一个误差范围。常见做法是用分位数回归训练三个模型分别预测 10%、50%、90% 分位数这样就能给出一个 80% 置信区间。LightGBM 支持分位数回归只需要把objective改成quantilealpha设成对应的分位数。def train_quantile_models(X_train, y_train, X_test, quantiles(0.1, 0.5, 0.9)): 训练分位数回归模型给出预测区间 models {} preds {} for q in quantiles: params { objective: quantile, alpha: q, metric: quantile, n_estimators: 600, learning_rate: 0.05, num_leaves: 31, min_child_samples: 30, random_state: 42 } model lgb.LGBMRegressor(**params) model.fit(X_train, y_train) models[q] model preds[q] model.predict(X_test) return models, preds models, preds train_quantile_models(X_train, y_train, X_test) interval_width preds[0.9] - preds[0.1] print(f80% 置信区间平均宽度: {interval_width.mean():.4f} m)这段代码训练三个分位数模型。参数说明alpha是目标分位数0.1 对应下界0.9 对应上界num_leaves降到 31 是因为分位数回归比均值回归更容易过拟合需要更简单的模型。interval_width是 80% 置信区间的宽度如果这个宽度在高风速下明显变大说明模型在高海况下的不确定性更高这是符合物理直觉的。我自己做 FY-3E 反演时最大的教训是不要迷信单一指标。MAE 小不代表模型能用交叉点验证稳、分层误差均衡、置信区间合理这三条都满足才敢往业务化推。希望帮到你。本文还有配套的精品资源点击获取