企业官网建设流程全解析

简介：这份资源面向电离层物理、空间天气与无线电传播方向的研究人员和学生，提供基于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 提取的核心链路：从双频观测到 STEC

2.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读 RINEX，numpy做数值计算。

import georinex as gr import numpy as np # 读取 RINEX 观测文件，指定要用的观测类型 obs = gr.load('rover0010.24o', use='G', tlim=None) # 提取 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, window=100): """相位平滑伪距，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, window=100) # 转换为 STEC（单位：TECU，1 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 需要自己估计或者用已知基准站做相对标定。如果完全不管 DCB，STEC 会有系统性偏差，典型值在 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_ion=450.0, R=6371.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/(R+h_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/(R+h_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 TECU，RMS 在 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，不要忽视周跳，不要用低仰角数据硬撑。希望帮到你。

本文还有配套的精品资源，点击获取