简介这份资源面向电离层物理、空间天气与无线电传播方向的研究人员和学生提供基于MATLAB的电离层数据处理与可视化脚本。核心文件codeTEC_L.m围绕CODEtec电离层电子密度与总电子含量编码技术展开可用于读取电子密度Ne与TEC观测数据完成清洗、缺失值处理与滤波等预处理并借助IRI等模型进行电离层建模与参数校准。脚本还涉及全球电离层分布的二维、三维绘图、地图投影与颜色映射以及日变化、季节变化和太阳活动关联的时间序列分析同时包含频谱分析与滤波等信号处理思路帮助评估电离层对GPS信号和通信链路的影响。压缩包为zip格式仅含1个m文件体积约1KB结构轻量便于直接放入MATLAB环境运行和二次修改。目前已有487人学习下载适合作为电离层数据分析和建模的入门脚本或算法参考。1. codeTEC_L 与 CODEtec从电离层原始观测到可发布 TEC 产品的落地路径做 GNSS 电离层处理的同行大概率都遇到过这个场景手头攒了一批双频观测文件想算电离层总电子含量TEC但打开 RTKLIB 或自己写的脚本一跑发现卫星弧段短、穿刺点分布稀疏、硬件偏差没扣干净出来的 TEC 序列跳变严重根本没法做后续的电离层建模或闪烁分析。codeTEC_L 和 CODEtec 这两个词本质上指向的就是这条链路——从原始 GNSS 观测到可用的电离层 TEC 产品。CODEtec 是伯尔尼大学 CODE 团队长期维护的电离层 TEC 产品体系codeTEC_L 则更偏向工程实现层面的代码化处理流程。这篇文章不讲教科书里的电离层物理只讲怎么把这条链路跑通数据怎么准备、TEC 怎么提取、硬件偏差怎么处理、平滑和粗差剔除怎么做、最终产品怎么验证。适合已经有一定 GNSS 数据处理基础、想自己搭一套 TEC 提取流程的工程师也适合做电离层建模、空间天气监测、PPP 电离层约束的从业者。2. 电离层 TEC 提取的核心链路从双频观测到 STEC2.1 为什么双频伪距和相位要分开用GNSS 双频观测提取 TEC 的基本原理是利用电离层色散效应不同频率的信号穿过电离层时延迟不同这个差异与 TEC 成正比。伪距观测的差分组合可以直接得到绝对 TEC但精度只有米级噪声大相位观测的差分组合精度高但存在整周模糊度只能得到相对 TEC。工程上的常见做法是用相位平滑伪距既保留伪距的绝对基准又利用相位的高精度。具体来说对于 GPS L1/L2 或北斗 B1/B3 双频数据无几何组合Geometry-Free Combination的伪距和相位分别写为P4 P1 - P2 40.3 * (1/f2^2 - 1/f1^2) * STEC DCB_r DCB_s ε_P L4 L1 - L2 -40.3 * (1/f2^2 - 1/f1^2) * STEC λ1*N1 - λ2*N2 ε_L其中 STEC 是斜路径总电子含量DCB_r 和 DCB_s 分别是接收机和卫星的差分码偏差N 是整周模糊度。伪距组合里 DCB 是必须扣掉的否则 TEC 会有几到几十 TECU 的系统性偏差。2.2 用 Python 实现相位平滑伪距的最小代码下面这段代码演示了从 RINEX 观测文件读取双频数据、计算 P4/L4、做相位平滑伪距的核心步骤。依赖georinex读 RINEXnumpy做数值计算。import georinex as gr import numpy as np # 读取 RINEX 观测文件指定要用的观测类型 obs gr.load(rover0010.24o, useG, tlimNone) # 提取 L1/L2 伪距和相位 P1 obs[C1C].values # GPS L1 C/A 伪距 P2 obs[C2W].values # GPS L2 P(Y) 伪距 L1 obs[L1C].values # GPS L1 相位 L2 obs[L2W].values # GPS L2 相位 f1, f2 1575.42e6, 1227.60e6 # GPS L1/L2 频率 alpha (f1**2 * f2**2) / (f1**2 - f2**2) # 频率因子 # 无几何组合 P4 P1 - P2 L4 L1 - L2 # 相位平滑伪距用相位变化量修正伪距 # 先计算相位组合的弧段均值再逐历元修正 def smooth_p4(P4, L4, window100): 相位平滑伪距window 为平滑窗口历元数 n_epoch P4.shape[0] P4_smooth np.full_like(P4, np.nan) for sat in range(P4.shape[1]): arc_start None for t in range(n_epoch): if not np.isnan(P4[t, sat]) and not np.isnan(L4[t, sat]): if arc_start is None: arc_start t # 弧段起始用伪距作为基准 base P4[t, sat] - L4[t, sat] P4_smooth[t, sat] base L4[t, sat] else: arc_start None return P4_smooth P4_smooth smooth_p4(P4, L4, window100) # 转换为 STEC单位TECU1 TECU 1e16 el/m^2 STEC P4_smooth * f1**2 * f2**2 / (40.3 * (f1**2 - f2**2)) / 1e16这段代码的关键点有三个第一base P4 - L4在弧段起始处计算把伪距的绝对基准和相位的相对变化绑定第二弧段中断失锁、周跳时必须重新初始化 base否则平滑会引入错误第三最终 STEC 的单位换算要统一到 TECU方便和 CODEtec 产品对比。参数方面平滑窗口window一般取 30 到 300 个历元。窗口太短伪距噪声压不住窗口太长电离层本身的变化会被平滑掉尤其在低纬度和赤道异常区TEC 变化率大窗口建议不超过 100。另外截止高度角一般设 15 到 30 度太低会引入多路径和大气折射误差。2.3 硬件偏差 DCB 的处理为什么不能跳过伪距组合 P4 里包含接收机 DCB 和卫星 DCB。卫星 DCB 可以从 IGS 或 CODE 的 DCB 产品里查接收机 DCB 需要自己估计或者用已知基准站做相对标定。如果完全不管 DCBSTEC 会有系统性偏差典型值在 5 到 20 TECU低纬度地区更大。工程上的常见做法是先用 CODE 的月度 DCB 产品扣掉卫星部分接收机 DCB 用最小二乘法在本地估计。具体来说把一天内所有卫星的 STEC 观测按穿刺点纬度做多项式拟合接收机 DCB 作为常数项一并估计。这一步在 codeTEC_L 这类流程里通常是独立模块输入是平滑后的 STEC输出是扣完 DCB 的绝对 TEC。提示如果只是做相对 TEC 变化分析比如电离层闪烁、行波扰动DCB 可以暂时不扣因为它是慢变常数不影响高频变化。但要做绝对 TEC 产品发布DCB 必须处理。3. 从 STEC 到 VTEC投影函数与穿刺点计算3.1 单层模型投影函数的选型STEC 是斜路径上的电子含量要转换成垂直方向的总电子含量VTEC需要投影函数。最常用的是单层模型Single Layer Model, SLM假设电离层集中在某个高度通常 350 到 450 km的薄层上。投影函数为VTEC STEC * cos(z)其中 z 是穿刺点处的天顶角。穿刺点经纬度由接收机位置、卫星方位角和仰角、以及单层高度计算得到。单层高度的选择会影响 VTEC 的精度。常见做法是取 350 km 或 400 km。低纬度地区电离层峰值高度偏高可以取 400 到 450 km中纬度取 350 km 就够。CODEtec 产品通常用 450 km 的单层高度做对比时要注意统一。3.2 穿刺点计算的 Python 实现import numpy as np def pierce_point(lat_r, lon_r, az, el, h_ion450.0, R6371.0): 计算穿刺点经纬度和天顶角 lat_r, lon_r: 接收机纬度和经度度 az, el: 卫星方位角和仰角度 h_ion: 单层高度km R: 地球平均半径km el_rad np.radians(el) az_rad np.radians(az) lat_r_rad np.radians(lat_r) # 地心角 psi np.pi/2 - el_rad - np.arcsin(R/(Rh_ion) * np.cos(el_rad)) # 穿刺点纬度 lat_pp np.arcsin(np.sin(lat_r_rad)*np.cos(psi) np.cos(lat_r_rad)*np.sin(psi)*np.cos(az_rad)) # 穿刺点经度 lon_pp np.radians(lon_r) np.arcsin(np.sin(psi)*np.sin(az_rad)/np.cos(lat_pp)) # 穿刺点天顶角 z_pp np.arcsin(R/(Rh_ion) * np.cos(el_rad)) return np.degrees(lat_pp), np.degrees(lon_pp), np.degrees(z_pp) # 示例接收机在北京39.9N, 116.4E卫星方位角 180 度仰角 60 度 lat_pp, lon_pp, z_pp pierce_point(39.9, 116.4, 180.0, 60.0) VTEC STEC * np.cos(np.radians(z_pp))这段代码里psi是接收机到穿刺点的地心角z_pp是穿刺点天顶角。注意np.arcsin里的参数要保证在 [-1, 1] 范围内否则会出现 NaN。实际处理中仰角低于截止高度角的观测直接丢弃避免投影函数在低仰角处放大误差。参数方面单层高度h_ion建议和你要对比的产品保持一致。如果做 CODEtec 对比用 450 km如果做本地建模350 到 400 km 都可以但整条链路要统一。地球半径R用 6371 km 是平均值精度足够。3.3 多系统融合时的频率差异处理GPS、北斗、GLONASS、Galileo 的频率不同P4/L4 组合的系数也不同。做多系统融合 TEC 时每个系统要单独计算 STEC然后统一投影到 VTEC。GLONASS 是频分多址不同卫星频率不同处理时要逐卫星用对应的频率值不能用一个固定频率。北斗三号 B1C/B2a 和 B1I/B3I 的频率也不一样RINEX 里观测类型要对应清楚。常见做法是先按系统分组每个系统内部按卫星逐颗处理得到各自的 VTEC 序列再按穿刺点位置做空间插值或格网化。如果只是做单站 TEC 时间序列多系统融合可以增加弧段覆盖减少数据空洞。4. 避坑与排查TEC 提取链路上的 5 个血泪教训4.1 周跳没检测平滑后 TEC 出现台阶现象VTEC 序列在某个历元突然跳变几到几十 TECU之后保持在新水平上。原因相位观测发生周跳L4 组合的整周模糊度变了但平滑算法没有重新初始化 base导致伪距基准被错误修正。解决在平滑前做周跳检测常用方法是电离层残差法L1 相位减 L2 相位的变化率或 TurboEdit。检测到周跳就断开弧段重新初始化。代码里可以在smooth_p4函数里加一个判断如果abs(L4[t] - L4[t-1]) threshold就把arc_start置为 None。4.2 DCB 产品时间不匹配TEC 整体偏移现象和 CODEtec 产品对比VTEC 整体偏高或偏低 5 到 15 TECU但变化趋势一致。原因用了错误月份的 DCB 文件或者卫星 DCB 的 PRN 编号和观测文件里的卫星编号不对应。解决下载 DCB 产品时确认时间覆盖观测时段PRN 编号要匹配。GPS 的 PRN 和 SVN 要区分清楚北斗的 C01 到 C60 也要和观测文件里的卫星号对齐。4.3 低仰角数据没剔除VTEC 噪声放大现象VTEC 序列在卫星升起和落下阶段噪声特别大甚至出现负值。原因低仰角时投影函数 cos(z) 很小STEC 的微小误差被放大同时多路径和大气折射误差也大。解决设置截止高度角一般 15 到 30 度。如果要做低仰角电离层研究需要用更精细的投影函数或直接用电离层层析方法不能简单用 SLM。4.4 接收机 DCB 估计时卫星弧段太短现象接收机 DCB 估计值不稳定不同时段差异大。原因参与估计的卫星弧段太短或者卫星几何分布不好导致最小二乘的法方程病态。解决至少用 6 颗以上卫星、弧段长度超过 2 小时的数据做估计如果单站数据不够可以用多站联合估计。另外估计前要把 VTEC 的日变化趋势扣掉只留常数项给 DCB。4.5 单位换算错误TECU 和 el/m^2 混用现象算出来的 TEC 数值和文献或产品差几个数量级。原因TECU 和 el/m^2 的换算搞错了。1 TECU 1e16 el/m^2。P4 的单位是米除以 40.3 和频率因子后得到的是 el/m^2再除以 1e16 才是 TECU。解决在代码里统一用 TECU 作为中间变量最后输出前再确认一次单位。建议在关键步骤打印数值范围正常 VTEC 在 0 到 100 TECU 之间超出这个范围基本就是单位或系数错了。5. 验证与进阶用 CODEtec 产品做交叉验证和精度评估5.1 和 CODEtec 最终产品的对比方法自己提取的 VTEC 要验证最直接的办法是和 CODEtec 的最终产品对比。CODEtec 提供全球电离层 TEC 格网产品时间分辨率有 1 小时和 5 分钟等不同版本。对比时注意几点第一单层高度要一致CODEtec 用 450 km第二时间系统要统一CODEtec 产品通常用 GPS 时第三空间上要把自己的穿刺点插值到 CODEtec 格网上或者反过来把格网值插值到穿刺点位置。对比指标用偏差Bias和均方根误差RMS。正常单站 VTEC 和 CODEtec 格网的偏差在 2 到 5 TECURMS 在 3 到 8 TECU。如果偏差超过 10 TECU大概率是 DCB 或投影函数的问题。低纬度地区偏差会大一些因为电离层梯度大格网插值误差增加。5.2 用双差 TEC 做精度自检如果没有外部产品对比可以用双差 TEC 做自检。同一颗卫星、两个接收机的 VTEC 之差在扣除接收机 DCB 差异后应该接近零短基线或反映真实的电离层梯度长基线。短基线下双差 TEC 的 RMS 可以反映提取精度正常在 0.5 到 2 TECU。如果双差 TEC 噪声很大说明平滑或周跳处理有问题。5.3 进阶技巧用球谐函数做区域 TEC 建模单站 VTEC 只是起点。如果有多站数据可以用球谐函数展开做区域或全球 TEC 建模。电离层 VTEC 在单层模型下可以表示为球谐级数VTEC(φ, λ) Σ Σ P_nm(sinφ) * (C_nm*cos(mλ) S_nm*sin(mλ))其中 P_nm 是勒让德函数C_nm 和 S_nm 是待估系数。用最小二乘求解同时估计接收机 DCB 和卫星 DCB。阶数一般取 2 到 4 阶太高会过拟合。这一步是 CODEtec 这类产品的核心算法自己实现时要注意法方程的条件数必要时加正则化。我自己的习惯是每次跑完 TEC 提取先画三张图——VTEC 时间序列、穿刺点轨迹图、和 CODEtec 的偏差散点图。三张图看完基本就知道这次处理有没有翻车。如果偏差散点图上有明显的系统性倾斜多半是投影函数或 DCB 的问题如果散点图上有离群点回去查周跳和粗差。这套流程跑了几年最深的教训就是不要跳过 DCB不要忽视周跳不要用低仰角数据硬撑。希望帮到你。本文还有配套的精品资源点击获取