尧图网页设计 
简介这份资源面向卫星遥感、海洋学与地球物理领域的科研人员及技术开发者围绕国产风云三号E星FY-3E搭载的GNOS-II仪器复现基于星载GNSS-R的海面高度反演模型。内容融合传统物理模型与机器学习方法采用随机森林和卷积神经网络对比北斗与GPS反射信号的反演精度并给出数据预处理、物理模型实现、模型训练评估的完整代码与解释适合希望快速上手国产卫星测高研究的读者。资源包共1个PDF文件约889KB集中呈现论文复现思路、代码实现与误差分析、空间分布可视化等关键环节。目前已有124人学习可作为从理论到实践的技术参考与实现框架。1. 从一张 FY-3E 的 DDM 图说起这套 GNSS-R 海面高度反演复现包到底能干什么海面高度这东西传统上要么靠验潮站要么靠卫星雷达高度计。验潮站是点测量出了近岸就抓瞎雷达高度计精度高但重访周期长、成本高而且近岸波形容易被陆地污染。GNSS-R 走的是另一条路——把导航卫星当成免费的 L 波段信号源用低轨卫星上的接收机接海面反射信号从反射信号里反推海面高度。FY-3E 作为国产极轨气象卫星搭载了 GNSS-R 载荷能同时接收 GPS、BDS、Galileo 的反射信号这就给海面高度反演提供了大量国产数据源。这份复现包的核心是把物理模型和机器学习两条路线揉在一起做海面高度反演。物理模型负责给出可解释的几何关系机器学习负责吃掉那些物理模型搞不定的残差和噪声。包里包含完整的代码、数据预处理脚本、模型训练流程和结果评估适合做遥感、GNSS-R、海洋测绘方向的研究生和工程师尤其是想拿国产卫星数据发论文、但又卡在“数据怎么处理、模型怎么搭、结果怎么验证”这三步上的人。下面我按“数据怎么进、模型怎么搭、坑在哪、怎么验证”的顺序把这份资源拆开讲。2. FY-3E GNSS-R 数据解析与 DDM 预处理从原始 L1 级数据到可训练样本2.1 为什么 FY-3E 的 DDM 不能直接拿来训练GNSS-R 的核心观测量是 DDMDelay-Doppler Map也就是反射信号在时延和多普勒两个维度上的功率分布。FY-3E 的 L1 级数据里DDM 是按镜面反射点附近的一小块区域存储的每个 DDM 对应一个镜面点附带几何参数卫星高度、入射角、反射角、接收机状态和辅助的海洋气象数据。问题在于原始 DDM 里混着相干分量和非相干分量海面粗糙度、风浪、电离层延迟都会影响功率分布直接拿原始功率值去回归海面高度模型学到的多半是噪声。常见做法是先做 DDM 的归一化和特征提取。归一化是为了消除接收机增益和距离衰减的影响特征提取则是把二维 DDM 压缩成几个物理意义明确的标量比如峰值功率、峰值时延、峰值多普勒、时延扩展、多普勒扩展、以及峰值周围的功率梯度。这些特征里峰值时延和镜面点几何关系直接挂钩是反演海面高度的主特征时延扩展和多普勒扩展反映海面粗糙度是修正项。我一般会先把 DDM 转成 dB 尺度再做峰值检测然后按镜面点位置对齐。这一步如果对齐做错后面所有特征都是歪的。import numpy as np from scipy.ndimage import maximum_filter def ddm_to_features(ddm_power, delay_bins, doppler_bins): ddm_power: 2D array, 原始 DDM 功率 (线性尺度) delay_bins: 1D array, 时延轴 (单位: 码片) doppler_bins: 1D array, 多普勒轴 (单位: Hz) 返回: 特征字典 # 转 dB 尺度避免动态范围过大 ddm_db 10 * np.log10(ddm_power 1e-12) # 找峰值位置 peak_idx np.unravel_index(np.argmax(ddm_db), ddm_db.shape) peak_delay delay_bins[peak_idx[0]] peak_doppler doppler_bins[peak_idx[1]] peak_power ddm_db[peak_idx] # 时延扩展: 峰值周围 -3dB 范围内的时延宽度 threshold peak_power - 3 mask ddm_db threshold delay_extent delay_bins[mask.any(axis1)].max() - delay_bins[mask.any(axis1)].min() # 多普勒扩展: 同理 doppler_extent doppler_bins[mask.any(axis0)].max() - doppler_bins[mask.any(axis0)].min() # 峰值梯度: 峰值点沿时延方向的一阶差分 delay_grad np.gradient(ddm_db, axis0)[peak_idx] return { peak_delay: peak_delay, peak_doppler: peak_doppler, peak_power: peak_power, delay_extent: delay_extent, doppler_extent: doppler_extent, delay_grad: delay_grad }这段代码的关键参数是threshold peak_power - 3也就是 -3dB 门限。门限设太高扩展量算出来偏小海面粗糙度信息丢失设太低噪声进来扩展量虚高。我试过 -2dB 到 -5dB-3dB 在 FY-3E 数据上比较稳。delay_bins和doppler_bins的单位要注意FY-3E 的 L1 产品里时延轴单位是码片多普勒轴是 Hz换算成物理距离和多普勒速度时要乘对应的系数。2.2 镜面点几何解算与高度初值计算DDM 特征只是输入真正要反演的是海面高度。物理模型给出的初值来自镜面点几何发射机导航卫星、接收机FY-3E、镜面反射点三者构成一个椭圆镜面点到接收机的路径时延和直达信号时延之差可以换算成镜面点相对于参考椭球的高度。这个高度是“几何高度”不是最终的海面高度但它是机器学习模型最重要的输入特征之一。具体步骤是先从 L1 数据里读出导航卫星和 FY-3E 的精密轨道位置再迭代求解镜面点位置然后算路径时延差。FY-3E 的精密轨道星历一般用 TLE 或 SP3 格式包里带了一个简化版的 SGP4 传播器精度够用但如果你要做高精度反演建议换成官方精密星历。from sgp4.api import Satrec import numpy as np def compute_specular_point(tle_line1, tle_line2, t, receiver_pos): tle_line1/2: TLE 两行根数 t: 时间 (datetime) receiver_pos: FY-3E 位置 (ECEF, km) 返回: 镜面点 ECEF 坐标和几何高度初值 sat Satrec.twoline2rv(tle_line1, tle_line2) jd t.timestamp() / 86400.0 2440587.5 e, r, v sat.sgp4(jd, 0.0) if e ! 0: raise RuntimeError(fSGP4 error: {e}) transmitter_pos np.array(r) # km, ECEF # 迭代求镜面点: 初始猜测为收发中点在地球表面的投影 midpoint (transmitter_pos receiver_pos) / 2.0 specular midpoint / np.linalg.norm(midpoint) * 6371.0 # 地球平均半径 for _ in range(10): # 更新镜面点: 满足入射角等于反射角 vec_t transmitter_pos - specular vec_r receiver_pos - specular normal specular / np.linalg.norm(specular) cos_i np.dot(vec_t, normal) / np.linalg.norm(vec_t) cos_r np.dot(vec_r, normal) / np.linalg.norm(vec_r) # 沿法线方向微调 specular specular normal * (cos_i - cos_r) * 10.0 specular specular / np.linalg.norm(specular) * 6371.0 # 几何高度初值: 路径时延差换算 range_t np.linalg.norm(transmitter_pos - specular) range_r np.linalg.norm(receiver_pos - specular) range_direct np.linalg.norm(transmitter_pos - receiver_pos) delay_diff (range_t range_r - range_direct) / 299792.458 # 秒 h_geo delay_diff * 299792.458 / 2.0 # 简化几何高度 return specular, h_geo这里迭代次数设 10 次是因为镜面点方程收敛快一般 5 到 7 次就稳定了。specular的初始猜测用收发中点投影在低仰角时偏差大但迭代能拉回来。h_geo是简化几何高度实际物理模型里还要考虑地球曲率和大气折射包里有一个atmospheric_correction函数做这件事参数用的是标准大气模型如果你有实测气象数据可以替换。提示FY-3E 的 TLE 更新频率不高做长时序反演时轨道误差会累积。我一般每 6 小时重新拉一次 TLE或者直接用 SP3 精密星历插值。3. 物理模型与机器学习融合双分支网络结构设计与训练策略3.1 物理约束怎么嵌进网络不是加个损失函数就完事物理模型和机器学习的融合常见做法有三种一是物理模型算初值机器学习做残差修正二是物理约束作为正则项加在损失函数里三是把物理方程嵌进网络结构做成可微分模块。这份复现包用的是第一种加第二种的混合物理模型给出几何高度初值和路径时延差作为额外输入特征同时损失函数里加了一项物理一致性约束让网络预测的高度和几何高度之间的偏差不要太大。为什么不用纯物理模型因为物理模型在近岸、高海况、低仰角时误差大这些恰恰是机器学习能补的地方。为什么不用纯机器学习因为纯数据驱动需要大量标注数据而海面高度的真值验潮站、雷达高度计在时空上稀疏纯 ML 容易过拟合到特定海域。网络结构是双分支一个分支吃 DDM 特征峰值时延、扩展、梯度等另一个分支吃几何特征入射角、镜面点纬度经度、几何高度初值、路径时延差。两个分支各接两层全连接然后拼接再接三层全连接回归出海面高度。激活函数用 ReLU输出层线性。import torch import torch.nn as nn class FusionNet(nn.Module): def __init__(self, ddm_feat_dim6, geo_feat_dim5, hidden64): super().__init__() # DDM 特征分支 self.ddm_branch nn.Sequential( nn.Linear(ddm_feat_dim, hidden), nn.ReLU(), nn.Linear(hidden, hidden), nn.ReLU() ) # 几何特征分支 self.geo_branch nn.Sequential( nn.Linear(geo_feat_dim, hidden), nn.ReLU(), nn.Linear(hidden, hidden), nn.ReLU() ) # 融合回归头 self.regressor nn.Sequential( nn.Linear(hidden * 2, hidden), nn.ReLU(), nn.Linear(hidden, hidden // 2), nn.ReLU(), nn.Linear(hidden // 2, 1) ) def forward(self, ddm_feat, geo_feat): ddm_out self.ddm_branch(ddm_feat) geo_out self.geo_branch(geo_feat) fused torch.cat([ddm_out, geo_out], dim-1) return self.regressor(fused).squeeze(-1) def physics_informed_loss(pred, target, h_geo, lambda_phys0.1): pred: 网络预测高度 target: 真值高度 h_geo: 物理模型几何高度初值 lambda_phys: 物理约束权重 mse nn.functional.mse_loss(pred, target) # 物理一致性: 预测值不应偏离几何高度太远 phys nn.functional.mse_loss(pred, h_geo) return mse lambda_phys * physlambda_phys这个权重是调参重点。设太大网络被物理初值拽住学不到残差设太小物理约束形同虚设。我在 FY-3E 数据上试过 0.01 到 1.00.1 左右比较平衡。hidden64是经验值数据量大可以加到 128但 FY-3E 单颗星的 DDM 样本量有限太大容易过拟合。3.2 训练集划分与数据增强别让同一片海域同时出现在训练和测试里遥感数据做机器学习最大的坑是时空自相关。同一片海域相邻时刻的样本高度相关如果随机划分训练集和测试集测试集里的样本可能在训练集里有“邻居”评估结果虚高。正确做法是按时间或空间划分比如用前 70% 时间的样本训练后 30% 测试或者按海域分太平洋训练大西洋测试。包里默认用的是时间划分split_by_time函数按时间戳排序后切分。数据增强方面DDM 特征可以加高斯噪声模拟接收机噪声几何特征可以做小幅度旋转模拟轨道误差。但要注意海面高度真值不能动增强只加在输入特征上。def split_by_time(features, labels, timestamps, train_ratio0.7): 按时间排序后切分避免时空泄漏 idx np.argsort(timestamps) n_train int(len(idx) * train_ratio) train_idx idx[:n_train] test_idx idx[n_train:] return (features[train_idx], labels[train_idx], features[test_idx], labels[test_idx]) def augment_ddm(ddm_feat, noise_std0.02): DDM 特征加高斯噪声 noise np.random.normal(0, noise_std, ddm_feat.shape) return ddm_feat noisenoise_std0.02是归一化后的尺度如果特征没归一化这个值要按特征量级调。时间划分的train_ratio我一般用 0.7但如果数据时间跨度短比如只有一个月建议用 0.8否则测试集样本太少评估不稳定。注意FY-3E 的 GNSS-R 数据在极区和高纬度海域覆盖密低纬度稀疏。按时间划分时如果训练期和测试期跨越的季节不同海况差异会导致模型泛化变差。我一般会检查训练集和测试集的海况分布有效波高、风速确保两者不要差太远。4. 避坑与排查FY-3E GNSS-R 反演里最容易翻车的五个地方4.1 现象训练损失降得很快但验证损失不降反升原因最常见的是时空泄漏。训练集和测试集里有同一片海域相邻时刻的样本模型记住了“这片海域的高度大概是多少”而不是学到反演关系。另一个可能是 DDM 特征归一化用了全局统计量而全局统计量里包含了测试集信息。解决按时间或空间划分数据集归一化统计量只用训练集算。包里normalize_features函数有一个fit_on_train_only参数默认 True但如果你手动改过数据流程要检查这个参数有没有被覆盖。4.2 现象反演高度在近岸区域偏差特别大远海还行原因近岸的 DDM 里混入了陆地反射相干分量被破坏物理模型的镜面点假设也不成立。机器学习模型如果训练集里近岸样本少就会把近岸样本当成噪声处理。解决在特征里加入“距离海岸线距离”和“陆地掩膜比例”两个辅助特征让模型知道当前样本是不是近岸。包里有一个coast_distance函数用的是 GSHHG 海岸线数据分辨率 1km够用。另外近岸样本在训练时可以适当加权sample_weight参数设成远海样本的 2 到 3 倍。4.3 现象物理约束损失项一直很大降不下去原因lambda_phys设太大或者几何高度初值本身误差大。FY-3E 的 TLE 轨道精度有限几何高度初值在低仰角时可能偏几十米物理约束硬拽着网络往错误方向走。解决先检查几何高度初值的误差分布如果初值本身偏差大把lambda_phys降到 0.01 以下或者改用“软约束”——只约束预测值和几何高度的差值不要超过某个阈值而不是强制接近。包里physics_informed_loss有一个huber_delta参数设成 10 米超过 10 米的偏差用线性损失避免梯度爆炸。4.4 现象模型在某个特定海域比如南海表现特别差原因训练集里这个海域的样本太少或者这个海域的海况比如内波、强流和训练集其他海域差异大。GNSS-R 对海面粗糙度敏感内波会改变粗糙度分布DDM 特征跟着变。解决做海域分层抽样确保每个海域在训练集里都有足够样本。如果某个海域样本实在少用迁移学习先在样本多的海域预训练再在目标海域微调。包里finetune脚本支持这个流程freeze_layers参数控制冻结哪几层。4.5 现象推理时单样本耗时太长没法做业务化原因双分支网络虽然不大但如果逐样本推理Python 循环开销大。另外DDM 特征提取里的maximum_filter和np.gradient在 CPU 上逐样本跑也慢。解决把特征提取和网络推理都批量化。maximum_filter可以换成torch.nn.functional.max_pool2d在 GPU 上批处理。网络推理用torch.no_grad()加batch_size256单样本耗时能从几十毫秒降到几毫秒。包里batch_inference函数做了这件事但要注意显存batch_size太大反而会 OOM。5. 验证与进阶用交叉验证和残差分析判断模型是不是真的学到了东西5.1 别只看 RMSE分海况、分仰角的误差拆解RMSE 是整体指标但海面高度反演的误差在海况和仰角上分布不均匀。我一般会把测试集按有效波高分成三档低、中、高按仰角分成三档低、中、高分别算 RMSE 和 bias。如果模型在低仰角或高海况下误差明显大说明特征里缺少对应信息或者训练集里这类样本太少。包里error_breakdown函数输出一个表格按海况和仰角交叉分组。下面是一个典型结果的结构海况档位仰角档位RMSE (m)Bias (m)样本数低低0.42-0.051200低中0.310.023500低高0.280.012800中低0.67-0.12900中中0.45-0.032600中高0.380.002100高低1.12-0.25400高中0.78-0.081100高高0.61-0.02800从这张表能看出低仰角加高海况是误差最大的组合RMSE 超过 1 米。如果业务上要求全海况 RMSE 小于 0.5 米这个模型还不够需要补高海况样本或者加海况修正项。5.2 残差分析模型没学到的部分长什么样残差预测值减真值如果只是随机噪声说明模型学到了主要关系如果残差和某个特征有相关性说明模型漏掉了这个特征的信息。我一般会画残差 vs 几何高度初值、残差 vs 有效波高、残差 vs 镜面点纬度三张图。包里residual_analysis函数会输出残差和这些特征的相关系数。如果残差和有效波高的相关系数超过 0.3说明模型对海况的修正不够可以考虑在特征里加有效波高的二次项或者把网络加宽。如果残差和纬度相关可能是轨道误差或电离层延迟有纬度依赖需要加对应的修正特征。def residual_analysis(pred, target, aux_features, feature_names): pred, target: 预测值和真值 aux_features: 辅助特征矩阵 (n_samples, n_features) feature_names: 特征名列表 residual pred - target print(fResidual mean: {residual.mean():.4f}, std: {residual.std():.4f}) for i, name in enumerate(feature_names): corr np.corrcoef(residual, aux_features[:, i])[0, 1] print(fCorr(residual, {name}) {corr:.4f}) return residual这个函数跑完如果某个特征的相关系数绝对值大于 0.2就值得深挖。我遇到过残差和镜面点经度相关系数 0.35 的情况后来发现是某个经度区间里 BDS 反射信号的信噪比偏低DDM 特征质量差加了一个信噪比阈值过滤后相关系数降到 0.1。5.3 交叉验证时间序列交叉验证比随机 K 折靠谱随机 K 折在时空数据上会高估性能我一般用时间序列交叉验证把时间轴分成 K 段每次用前 K-1 段训练最后一段测试然后滑动。包里time_series_cv函数实现了这个逻辑n_splits5是默认值。如果数据时间跨度小于 3 个月n_splits设 3 就行再多每折样本太少。交叉验证的另一个用处是调参。lambda_phys、hidden、learning_rate这三个参数对结果影响最大我一般先用时间序列 CV 粗扫一遍再在最优附近细调。粗扫范围lambda_phys取 [0.01, 0.05, 0.1, 0.5]hidden取 [32, 64, 128]learning_rate取 [1e-4, 5e-4, 1e-3]。细调时固定两个调一个。从那以后我每次做 GNSS-R 反演都强制先跑一遍时间序列 CV 和残差分析再去看整体 RMSE。整体指标好看但残差有结构的模型上了业务迟早翻车。希望帮到你。本文还有配套的精品资源点击获取 FY-3E GNSS-R海面高度反演：物理与机器学习融合复现包
简介这份资源面向卫星遥感、海洋学与地球物理领域的科研人员及技术开发者围绕国产风云三号E星FY-3E搭载的GNOS-II仪器复现基于星载GNSS-R的海面高度反演模型。内容融合传统物理模型与机器学习方法采用随机森林和卷积神经网络对比北斗与GPS反射信号的反演精度并给出数据预处理、物理模型实现、模型训练评估的完整代码与解释适合希望快速上手国产卫星测高研究的读者。资源包共1个PDF文件约889KB集中呈现论文复现思路、代码实现与误差分析、空间分布可视化等关键环节。目前已有124人学习可作为从理论到实践的技术参考与实现框架。1. 从一张 FY-3E 的 DDM 图说起这套 GNSS-R 海面高度反演复现包到底能干什么海面高度这东西传统上要么靠验潮站要么靠卫星雷达高度计。验潮站是点测量出了近岸就抓瞎雷达高度计精度高但重访周期长、成本高而且近岸波形容易被陆地污染。GNSS-R 走的是另一条路——把导航卫星当成免费的 L 波段信号源用低轨卫星上的接收机接海面反射信号从反射信号里反推海面高度。FY-3E 作为国产极轨气象卫星搭载了 GNSS-R 载荷能同时接收 GPS、BDS、Galileo 的反射信号这就给海面高度反演提供了大量国产数据源。这份复现包的核心是把物理模型和机器学习两条路线揉在一起做海面高度反演。物理模型负责给出可解释的几何关系机器学习负责吃掉那些物理模型搞不定的残差和噪声。包里包含完整的代码、数据预处理脚本、模型训练流程和结果评估适合做遥感、GNSS-R、海洋测绘方向的研究生和工程师尤其是想拿国产卫星数据发论文、但又卡在“数据怎么处理、模型怎么搭、结果怎么验证”这三步上的人。下面我按“数据怎么进、模型怎么搭、坑在哪、怎么验证”的顺序把这份资源拆开讲。2. FY-3E GNSS-R 数据解析与 DDM 预处理从原始 L1 级数据到可训练样本2.1 为什么 FY-3E 的 DDM 不能直接拿来训练GNSS-R 的核心观测量是 DDMDelay-Doppler Map也就是反射信号在时延和多普勒两个维度上的功率分布。FY-3E 的 L1 级数据里DDM 是按镜面反射点附近的一小块区域存储的每个 DDM 对应一个镜面点附带几何参数卫星高度、入射角、反射角、接收机状态和辅助的海洋气象数据。问题在于原始 DDM 里混着相干分量和非相干分量海面粗糙度、风浪、电离层延迟都会影响功率分布直接拿原始功率值去回归海面高度模型学到的多半是噪声。常见做法是先做 DDM 的归一化和特征提取。归一化是为了消除接收机增益和距离衰减的影响特征提取则是把二维 DDM 压缩成几个物理意义明确的标量比如峰值功率、峰值时延、峰值多普勒、时延扩展、多普勒扩展、以及峰值周围的功率梯度。这些特征里峰值时延和镜面点几何关系直接挂钩是反演海面高度的主特征时延扩展和多普勒扩展反映海面粗糙度是修正项。我一般会先把 DDM 转成 dB 尺度再做峰值检测然后按镜面点位置对齐。这一步如果对齐做错后面所有特征都是歪的。import numpy as np from scipy.ndimage import maximum_filter def ddm_to_features(ddm_power, delay_bins, doppler_bins): ddm_power: 2D array, 原始 DDM 功率 (线性尺度) delay_bins: 1D array, 时延轴 (单位: 码片) doppler_bins: 1D array, 多普勒轴 (单位: Hz) 返回: 特征字典 # 转 dB 尺度避免动态范围过大 ddm_db 10 * np.log10(ddm_power 1e-12) # 找峰值位置 peak_idx np.unravel_index(np.argmax(ddm_db), ddm_db.shape) peak_delay delay_bins[peak_idx[0]] peak_doppler doppler_bins[peak_idx[1]] peak_power ddm_db[peak_idx] # 时延扩展: 峰值周围 -3dB 范围内的时延宽度 threshold peak_power - 3 mask ddm_db threshold delay_extent delay_bins[mask.any(axis1)].max() - delay_bins[mask.any(axis1)].min() # 多普勒扩展: 同理 doppler_extent doppler_bins[mask.any(axis0)].max() - doppler_bins[mask.any(axis0)].min() # 峰值梯度: 峰值点沿时延方向的一阶差分 delay_grad np.gradient(ddm_db, axis0)[peak_idx] return { peak_delay: peak_delay, peak_doppler: peak_doppler, peak_power: peak_power, delay_extent: delay_extent, doppler_extent: doppler_extent, delay_grad: delay_grad }这段代码的关键参数是threshold peak_power - 3也就是 -3dB 门限。门限设太高扩展量算出来偏小海面粗糙度信息丢失设太低噪声进来扩展量虚高。我试过 -2dB 到 -5dB-3dB 在 FY-3E 数据上比较稳。delay_bins和doppler_bins的单位要注意FY-3E 的 L1 产品里时延轴单位是码片多普勒轴是 Hz换算成物理距离和多普勒速度时要乘对应的系数。2.2 镜面点几何解算与高度初值计算DDM 特征只是输入真正要反演的是海面高度。物理模型给出的初值来自镜面点几何发射机导航卫星、接收机FY-3E、镜面反射点三者构成一个椭圆镜面点到接收机的路径时延和直达信号时延之差可以换算成镜面点相对于参考椭球的高度。这个高度是“几何高度”不是最终的海面高度但它是机器学习模型最重要的输入特征之一。具体步骤是先从 L1 数据里读出导航卫星和 FY-3E 的精密轨道位置再迭代求解镜面点位置然后算路径时延差。FY-3E 的精密轨道星历一般用 TLE 或 SP3 格式包里带了一个简化版的 SGP4 传播器精度够用但如果你要做高精度反演建议换成官方精密星历。from sgp4.api import Satrec import numpy as np def compute_specular_point(tle_line1, tle_line2, t, receiver_pos): tle_line1/2: TLE 两行根数 t: 时间 (datetime) receiver_pos: FY-3E 位置 (ECEF, km) 返回: 镜面点 ECEF 坐标和几何高度初值 sat Satrec.twoline2rv(tle_line1, tle_line2) jd t.timestamp() / 86400.0 2440587.5 e, r, v sat.sgp4(jd, 0.0) if e ! 0: raise RuntimeError(fSGP4 error: {e}) transmitter_pos np.array(r) # km, ECEF # 迭代求镜面点: 初始猜测为收发中点在地球表面的投影 midpoint (transmitter_pos receiver_pos) / 2.0 specular midpoint / np.linalg.norm(midpoint) * 6371.0 # 地球平均半径 for _ in range(10): # 更新镜面点: 满足入射角等于反射角 vec_t transmitter_pos - specular vec_r receiver_pos - specular normal specular / np.linalg.norm(specular) cos_i np.dot(vec_t, normal) / np.linalg.norm(vec_t) cos_r np.dot(vec_r, normal) / np.linalg.norm(vec_r) # 沿法线方向微调 specular specular normal * (cos_i - cos_r) * 10.0 specular specular / np.linalg.norm(specular) * 6371.0 # 几何高度初值: 路径时延差换算 range_t np.linalg.norm(transmitter_pos - specular) range_r np.linalg.norm(receiver_pos - specular) range_direct np.linalg.norm(transmitter_pos - receiver_pos) delay_diff (range_t range_r - range_direct) / 299792.458 # 秒 h_geo delay_diff * 299792.458 / 2.0 # 简化几何高度 return specular, h_geo这里迭代次数设 10 次是因为镜面点方程收敛快一般 5 到 7 次就稳定了。specular的初始猜测用收发中点投影在低仰角时偏差大但迭代能拉回来。h_geo是简化几何高度实际物理模型里还要考虑地球曲率和大气折射包里有一个atmospheric_correction函数做这件事参数用的是标准大气模型如果你有实测气象数据可以替换。提示FY-3E 的 TLE 更新频率不高做长时序反演时轨道误差会累积。我一般每 6 小时重新拉一次 TLE或者直接用 SP3 精密星历插值。3. 物理模型与机器学习融合双分支网络结构设计与训练策略3.1 物理约束怎么嵌进网络不是加个损失函数就完事物理模型和机器学习的融合常见做法有三种一是物理模型算初值机器学习做残差修正二是物理约束作为正则项加在损失函数里三是把物理方程嵌进网络结构做成可微分模块。这份复现包用的是第一种加第二种的混合物理模型给出几何高度初值和路径时延差作为额外输入特征同时损失函数里加了一项物理一致性约束让网络预测的高度和几何高度之间的偏差不要太大。为什么不用纯物理模型因为物理模型在近岸、高海况、低仰角时误差大这些恰恰是机器学习能补的地方。为什么不用纯机器学习因为纯数据驱动需要大量标注数据而海面高度的真值验潮站、雷达高度计在时空上稀疏纯 ML 容易过拟合到特定海域。网络结构是双分支一个分支吃 DDM 特征峰值时延、扩展、梯度等另一个分支吃几何特征入射角、镜面点纬度经度、几何高度初值、路径时延差。两个分支各接两层全连接然后拼接再接三层全连接回归出海面高度。激活函数用 ReLU输出层线性。import torch import torch.nn as nn class FusionNet(nn.Module): def __init__(self, ddm_feat_dim6, geo_feat_dim5, hidden64): super().__init__() # DDM 特征分支 self.ddm_branch nn.Sequential( nn.Linear(ddm_feat_dim, hidden), nn.ReLU(), nn.Linear(hidden, hidden), nn.ReLU() ) # 几何特征分支 self.geo_branch nn.Sequential( nn.Linear(geo_feat_dim, hidden), nn.ReLU(), nn.Linear(hidden, hidden), nn.ReLU() ) # 融合回归头 self.regressor nn.Sequential( nn.Linear(hidden * 2, hidden), nn.ReLU(), nn.Linear(hidden, hidden // 2), nn.ReLU(), nn.Linear(hidden // 2, 1) ) def forward(self, ddm_feat, geo_feat): ddm_out self.ddm_branch(ddm_feat) geo_out self.geo_branch(geo_feat) fused torch.cat([ddm_out, geo_out], dim-1) return self.regressor(fused).squeeze(-1) def physics_informed_loss(pred, target, h_geo, lambda_phys0.1): pred: 网络预测高度 target: 真值高度 h_geo: 物理模型几何高度初值 lambda_phys: 物理约束权重 mse nn.functional.mse_loss(pred, target) # 物理一致性: 预测值不应偏离几何高度太远 phys nn.functional.mse_loss(pred, h_geo) return mse lambda_phys * physlambda_phys这个权重是调参重点。设太大网络被物理初值拽住学不到残差设太小物理约束形同虚设。我在 FY-3E 数据上试过 0.01 到 1.00.1 左右比较平衡。hidden64是经验值数据量大可以加到 128但 FY-3E 单颗星的 DDM 样本量有限太大容易过拟合。3.2 训练集划分与数据增强别让同一片海域同时出现在训练和测试里遥感数据做机器学习最大的坑是时空自相关。同一片海域相邻时刻的样本高度相关如果随机划分训练集和测试集测试集里的样本可能在训练集里有“邻居”评估结果虚高。正确做法是按时间或空间划分比如用前 70% 时间的样本训练后 30% 测试或者按海域分太平洋训练大西洋测试。包里默认用的是时间划分split_by_time函数按时间戳排序后切分。数据增强方面DDM 特征可以加高斯噪声模拟接收机噪声几何特征可以做小幅度旋转模拟轨道误差。但要注意海面高度真值不能动增强只加在输入特征上。def split_by_time(features, labels, timestamps, train_ratio0.7): 按时间排序后切分避免时空泄漏 idx np.argsort(timestamps) n_train int(len(idx) * train_ratio) train_idx idx[:n_train] test_idx idx[n_train:] return (features[train_idx], labels[train_idx], features[test_idx], labels[test_idx]) def augment_ddm(ddm_feat, noise_std0.02): DDM 特征加高斯噪声 noise np.random.normal(0, noise_std, ddm_feat.shape) return ddm_feat noisenoise_std0.02是归一化后的尺度如果特征没归一化这个值要按特征量级调。时间划分的train_ratio我一般用 0.7但如果数据时间跨度短比如只有一个月建议用 0.8否则测试集样本太少评估不稳定。注意FY-3E 的 GNSS-R 数据在极区和高纬度海域覆盖密低纬度稀疏。按时间划分时如果训练期和测试期跨越的季节不同海况差异会导致模型泛化变差。我一般会检查训练集和测试集的海况分布有效波高、风速确保两者不要差太远。4. 避坑与排查FY-3E GNSS-R 反演里最容易翻车的五个地方4.1 现象训练损失降得很快但验证损失不降反升原因最常见的是时空泄漏。训练集和测试集里有同一片海域相邻时刻的样本模型记住了“这片海域的高度大概是多少”而不是学到反演关系。另一个可能是 DDM 特征归一化用了全局统计量而全局统计量里包含了测试集信息。解决按时间或空间划分数据集归一化统计量只用训练集算。包里normalize_features函数有一个fit_on_train_only参数默认 True但如果你手动改过数据流程要检查这个参数有没有被覆盖。4.2 现象反演高度在近岸区域偏差特别大远海还行原因近岸的 DDM 里混入了陆地反射相干分量被破坏物理模型的镜面点假设也不成立。机器学习模型如果训练集里近岸样本少就会把近岸样本当成噪声处理。解决在特征里加入“距离海岸线距离”和“陆地掩膜比例”两个辅助特征让模型知道当前样本是不是近岸。包里有一个coast_distance函数用的是 GSHHG 海岸线数据分辨率 1km够用。另外近岸样本在训练时可以适当加权sample_weight参数设成远海样本的 2 到 3 倍。4.3 现象物理约束损失项一直很大降不下去原因lambda_phys设太大或者几何高度初值本身误差大。FY-3E 的 TLE 轨道精度有限几何高度初值在低仰角时可能偏几十米物理约束硬拽着网络往错误方向走。解决先检查几何高度初值的误差分布如果初值本身偏差大把lambda_phys降到 0.01 以下或者改用“软约束”——只约束预测值和几何高度的差值不要超过某个阈值而不是强制接近。包里physics_informed_loss有一个huber_delta参数设成 10 米超过 10 米的偏差用线性损失避免梯度爆炸。4.4 现象模型在某个特定海域比如南海表现特别差原因训练集里这个海域的样本太少或者这个海域的海况比如内波、强流和训练集其他海域差异大。GNSS-R 对海面粗糙度敏感内波会改变粗糙度分布DDM 特征跟着变。解决做海域分层抽样确保每个海域在训练集里都有足够样本。如果某个海域样本实在少用迁移学习先在样本多的海域预训练再在目标海域微调。包里finetune脚本支持这个流程freeze_layers参数控制冻结哪几层。4.5 现象推理时单样本耗时太长没法做业务化原因双分支网络虽然不大但如果逐样本推理Python 循环开销大。另外DDM 特征提取里的maximum_filter和np.gradient在 CPU 上逐样本跑也慢。解决把特征提取和网络推理都批量化。maximum_filter可以换成torch.nn.functional.max_pool2d在 GPU 上批处理。网络推理用torch.no_grad()加batch_size256单样本耗时能从几十毫秒降到几毫秒。包里batch_inference函数做了这件事但要注意显存batch_size太大反而会 OOM。5. 验证与进阶用交叉验证和残差分析判断模型是不是真的学到了东西5.1 别只看 RMSE分海况、分仰角的误差拆解RMSE 是整体指标但海面高度反演的误差在海况和仰角上分布不均匀。我一般会把测试集按有效波高分成三档低、中、高按仰角分成三档低、中、高分别算 RMSE 和 bias。如果模型在低仰角或高海况下误差明显大说明特征里缺少对应信息或者训练集里这类样本太少。包里error_breakdown函数输出一个表格按海况和仰角交叉分组。下面是一个典型结果的结构海况档位仰角档位RMSE (m)Bias (m)样本数低低0.42-0.051200低中0.310.023500低高0.280.012800中低0.67-0.12900中中0.45-0.032600中高0.380.002100高低1.12-0.25400高中0.78-0.081100高高0.61-0.02800从这张表能看出低仰角加高海况是误差最大的组合RMSE 超过 1 米。如果业务上要求全海况 RMSE 小于 0.5 米这个模型还不够需要补高海况样本或者加海况修正项。5.2 残差分析模型没学到的部分长什么样残差预测值减真值如果只是随机噪声说明模型学到了主要关系如果残差和某个特征有相关性说明模型漏掉了这个特征的信息。我一般会画残差 vs 几何高度初值、残差 vs 有效波高、残差 vs 镜面点纬度三张图。包里residual_analysis函数会输出残差和这些特征的相关系数。如果残差和有效波高的相关系数超过 0.3说明模型对海况的修正不够可以考虑在特征里加有效波高的二次项或者把网络加宽。如果残差和纬度相关可能是轨道误差或电离层延迟有纬度依赖需要加对应的修正特征。def residual_analysis(pred, target, aux_features, feature_names): pred, target: 预测值和真值 aux_features: 辅助特征矩阵 (n_samples, n_features) feature_names: 特征名列表 residual pred - target print(fResidual mean: {residual.mean():.4f}, std: {residual.std():.4f}) for i, name in enumerate(feature_names): corr np.corrcoef(residual, aux_features[:, i])[0, 1] print(fCorr(residual, {name}) {corr:.4f}) return residual这个函数跑完如果某个特征的相关系数绝对值大于 0.2就值得深挖。我遇到过残差和镜面点经度相关系数 0.35 的情况后来发现是某个经度区间里 BDS 反射信号的信噪比偏低DDM 特征质量差加了一个信噪比阈值过滤后相关系数降到 0.1。5.3 交叉验证时间序列交叉验证比随机 K 折靠谱随机 K 折在时空数据上会高估性能我一般用时间序列交叉验证把时间轴分成 K 段每次用前 K-1 段训练最后一段测试然后滑动。包里time_series_cv函数实现了这个逻辑n_splits5是默认值。如果数据时间跨度小于 3 个月n_splits设 3 就行再多每折样本太少。交叉验证的另一个用处是调参。lambda_phys、hidden、learning_rate这三个参数对结果影响最大我一般先用时间序列 CV 粗扫一遍再在最优附近细调。粗扫范围lambda_phys取 [0.01, 0.05, 0.1, 0.5]hidden取 [32, 64, 128]learning_rate取 [1e-4, 5e-4, 1e-3]。细调时固定两个调一个。从那以后我每次做 GNSS-R 反演都强制先跑一遍时间序列 CV 和残差分析再去看整体 RMSE。整体指标好看但残差有结构的模型上了业务迟早翻车。希望帮到你。本文还有配套的精品资源点击获取