简介这份资源面向从事GNSS电离层分析与GPS定位精度研究的学习者和工程人员提供一套基于RINEX文件的电离层总电子含量TEC计算工具。它能够解析RINEX 2.11观测数据提取GPS双频伪距与相位观测值通过差分处理消除对流层延迟与钟差影响并结合DCB改正、周跳探测与相位平滑等步骤最终输出TEC随时间变化的曲线与分布图可用于空间天气监测、电离层特性研究及定位误差校正等场景。压缩包共16个文件全部为m脚本文件整体约16KB涵盖卫星位置计算、DCB读取、周跳改正、IPP映射、TEC计算与绘图等完整流程模块结构紧凑、便于二次开发与移植。目前已有325人学习下载。读者可借此快速搭建从RINEX读取到TEC成图的分析链路理解双频差分与电离层建模的关键实现并在此基础上开展自己的电离层实验与算法验证。1. 从一份 RINEX 观测文件算出 TEC这个标题到底在解决什么问题电离层总电子含量TEC是 GNSS 信号传播路径上自由电子密度的积分量单位是 TECU1 TECU 等于 10¹⁶ 个电子每平方米。做定位的人关心它是因为电离层延迟是单频 GNSS 最大的误差源之一做空间天气的人关心它是因为 TEC 的时空变化直接反映电离层活动。而Calculate Tec_GPSRINEX_TEC_这个标题说的就是一件很具体的事拿一份 GPS 的 RINEX 观测文件从双频伪距和载波相位里把 TEC 算出来。听起来不复杂但真正动手你会发现几个绕不开的问题RINEX 版本不同字段布局不一样GPS 的 TEC 公式里有个卫星和接收机之间的差分码偏差DCB必须扣掉载波相位算出来的 TEC 精度高但有整周模糊度伪距算出来的没有模糊度但噪声大。怎么把两者平滑到一起怎么处理 DCB怎么把斜向 TEC 映射到垂直 TEC这些才是决定你算出来的东西能不能用的关键。这篇文章面向的是手里已经有 RINEX 文件、想自己算 TEC 的工程师和学生。我会从文件解析讲到 TEC 计算公式再到 DCB 改正和平滑算法最后给出可复现的代码和参数设置。不依赖任何特定软件包用 Python 从零搭一遍你能看到每一步在做什么、参数为什么这么设。2. RINEX 观测文件解析与双频数据提取2.1 RINEX 版本差异与观测码选择RINEX 观测文件目前常见的有 2.11 和 3.x 两个大版本。2.11 的观测记录按卫星逐行排列每行固定 80 字符观测类型在文件头的# / TYPES OF OBSERV里定义。3.x 改成了按历元、按卫星的 XML 风格观测码用C1C、L1C、C2W、L2W这种三字符标识含义更明确但解析逻辑完全不同。算 TEC 需要的是双频伪距和双频载波相位。GPS 上常用的是 L1 和 L2伪距对应C1C/C1W和C2W/C2C载波相位对应L1C/L1W和L2W/L2C。如果你拿到的文件里 L2 只有C2W没有L2W那载波相位那一路就没法用只能退回伪距单路算 TEC精度会差一个量级。我一般拿到文件先看文件头确认版本、观测类型、采样率、测站坐标。这几项决定了后面所有参数怎么设。采样率 1Hz 和 30s 的平滑窗口长度完全不同测站坐标用来算卫星高度角和方位角进而做斜向到垂直的映射。# 读取 RINEX 文件头提取版本、观测类型、测站坐标 def parse_rinex_header(filepath): header { version: None, obs_types: [], approx_xyz: [0.0, 0.0, 0.0], interval: 30.0 } with open(filepath, r) as f: for line in f: label line[60:].strip() if RINEX VERSION / TYPE in label: header[version] float(line[0:9]) elif # / TYPES OF OBSERV in label: # 2.11 格式每行最多 9 个观测码可能跨多行 n int(line[0:6]) codes line[6:60].split() header[obs_types].extend(codes) elif APPROX POSITION XYZ in label: header[approx_xyz] [ float(line[0:14]), float(line[14:28]), float(line[28:42]) ] elif INTERVAL in label: header[interval] float(line[0:10]) elif END OF HEADER in label: break return header这段代码只处理了 2.11 的头解析。3.x 的头结构类似但观测码定义在SYS / # / OBS TYPES里按系统分组。实际用的时候我会先判断版本号再走不同分支。approx_xyz是测站在 WGS84 下的近似坐标精度到米级就够用后面算高度角不需要太准。2.2 按历元提取双频伪距与相位解析完头之后正文部分要按历元逐条读。2.11 的格式里每个历元第一行是时间戳和卫星列表后面每颗卫星的观测值按观测类型顺序排列每个观测值占 14 字符加一个 LLI 标志和一个信号强度。这里最容易翻车的地方是卫星列表超过 12 颗时会续行续行没有时间戳解析时要把续行的卫星拼回去。# 按历元提取指定卫星的 L1/L2 伪距和载波相位 def extract_dual_freq(filepath, header, target_prnG01): epochs [] with open(filepath, r) as f: # 跳过文件头 for line in f: if END OF HEADER in line: break current_epoch None sat_list [] while True: line f.readline() if not line: break # 历元头行年份在 1-3 列 if line[0:2].strip().isdigit() or line[0:1] : try: year int(line[1:3]) except ValueError: continue # 解析卫星列表处理续行 n_sat int(line[29:32]) sats [] idx 32 while len(sats) n_sat: sats.extend([line[idxi*3:idxi*33] for i in range(3) if line[idxi*3:idxi*33].strip()]) if len(sats) n_sat: line f.readline() idx 32 sat_list sats current_epoch {time: line[0:26], data: {}} epochs.append(current_epoch) # 读取每颗卫星的观测值 for sat in sat_list: obs_line f.readline() obs {} for i, code in enumerate(header[obs_types]): val_str obs_line[i*16:i*1614].strip() if val_str: obs[code] float(val_str) current_epoch[data][sat] obs return epochs这段代码的核心逻辑是先读历元头拿到卫星列表再按卫星数逐行读观测值。i*16这个偏移量是 2.11 的固定格式每个观测值占 16 字符14 数值 1 LLI 1 信号强度。如果你用的是 3.x 文件这个偏移量要改成按观测码动态计算。提取到C1C、C2W、L1C、L2W之后就可以进入 TEC 计算了。这里有个细节伪距和相位的单位。伪距是米相位是周算 TEC 之前相位要乘以波长转成米。L1 波长约 0.1903 米L2 约 0.2442 米这两个常数在公式里反复出现建议定义成模块级常量。3. TEC 计算公式与 DCB 改正的落地实现3.1 双频伪距与相位算 TEC 的公式推导GPS 双频观测里L1 和 L2 的电离层延迟与频率平方成反比。利用这个关系可以把两个频率上的伪距观测组合起来消掉大部分非电离层项剩下的就是电离层延迟。伪距组合算 TEC 的公式是TEC_p (f1² · f2²) / (40.3 · (f1² - f2²)) · (P2 - P1)其中 f1 1575.42 MHzf2 1227.60 MHz40.3 是常数P1 和 P2 是 L1 和 L2 的伪距观测值单位米。算出来的 TEC 单位是 TECU。载波相位组合类似但相位观测里包含整周模糊度所以算出来的 TEC 有一个未知的常数偏差TEC_φ (f1² · f2²) / (40.3 · (f1² - f2²)) · (L1·λ1 - L2·λ2) const伪距 TEC 没有模糊度但噪声大典型噪声在 1-3 TECU相位 TEC 噪声小两个量级但有个未知常数。工程上的标准做法是用相位 TEC 的變化量去平滑伪距 TEC 的绝对值这就是载波相位平滑伪距的思路。import numpy as np F1 1575.42e6 # L1 频率 Hz F2 1227.60e6 # L2 频率 Hz C 299792458.0 # 光速 m/s LAM1 C / F1 # L1 波长 LAM2 C / F2 # L2 波长 K 40.3e16 # 电离层常数配合 TECU 使用 def tec_from_pseudorange(P1, P2): 双频伪距算斜向 TEC单位 TECU coef (F1**2 * F2**2) / (K * (F1**2 - F2**2)) return coef * (P2 - P1) def tec_from_phase(L1, L2): 双频载波相位算斜向 TEC 的相对值单位 TECU coef (F1**2 * F2**2) / (K * (F1**2 - F2**2)) # L1/L2 是周数先转米 return coef * (L1 * LAM1 - L2 * LAM2)K 40.3e16这个写法是为了让结果直接落在 TECU 量级。如果你把 40.3 直接写进去算出来的单位是电子每平方米数值会大到没法看。这个常数在文献里有两种写法注意区分。3.2 DCB 改正卫星和接收机偏差怎么扣上面算出来的 TEC 里包含卫星和接收机的差分码偏差DCB。卫星 DCB 是卫星端 P1 和 P2 硬件延迟的差接收机 DCB 是接收机端的对应差值。这两个偏差如果不扣TEC 会有几个到十几个 TECU 的系统性偏移做绝对 TEC 分析时完全不能用。DCB 改正的公式是TEC_corrected TEC_raw - (DCB_sat DCB_recv) / (K · (f1² - f2²) / (f1² · f2²))实际工程里卫星 DCB 可以从公开的 DCB 产品文件里查接收机 DCB 一般用最小二乘估计或者用已知的 TEC 模型反算。如果你只是做相对变化分析不扣 DCB 也能看出趋势但要做绝对 TEC 或者和测高仪对比DCB 必须扣。def apply_dcb(tec_raw, dcb_sat, dcb_recv): 扣 DCB 改正 dcb_sat / dcb_recv 单位米P1-P2 的硬件延迟差 coef (F1**2 * F2**2) / (K * (F1**2 - F2**2)) return tec_raw - coef * (dcb_sat dcb_recv)这里dcb_sat和dcb_recv的单位是米不是 TECU。很多 DCB 产品给的是纳秒要乘以光速转成米。我见过有人直接把纳秒代进去结果 TEC 偏了几百个 TECU排查了半天才发现是单位问题。3.3 载波相位平滑伪距的实现与窗口选择相位平滑伪距的基本思路是用相位 TEC 的历元间变化量去修正伪距 TEC 的绝对值。具体做法是在第一个历元用伪距 TEC 作为基准之后每个历元用相位 TEC 的变化量累加同时用递推公式把伪距的噪声压下去。经典的 Hatch 平滑公式是TEC_smooth(k) (1/N) · TEC_p(k) ((N-1)/N) · (TEC_smooth(k-1) ΔTEC_φ(k))其中 N 是平滑窗口长度ΔTEC_φ(k) 是当前历元相对上一历元的相位 TEC 变化量。N 越大平滑效果越好但相位失锁时会引入偏差。工程上 N 一般取 30 到 100对应 30s 采样率下 15 到 50 分钟。def smooth_tec(tec_p, tec_phi, window50): 载波相位平滑伪距 TEC tec_p: 伪距 TEC 序列 tec_phi: 相位 TEC 序列相对值 window: 平滑窗口长度 n len(tec_p) tec_smooth np.zeros(n) tec_smooth[0] tec_p[0] for k in range(1, n): d_phi tec_phi[k] - tec_phi[k-1] w min(k, window) tec_smooth[k] (1.0/w) * tec_p[k] ((w-1.0)/w) * (tec_smooth[k-1] d_phi) return tec_smooth窗口长度的选择要看采样率和电离层活动程度。采样率 1Hz 时窗口可以取 100 到 30030s 采样率取 30 到 60 比较合适。电离层活跃期相位失锁频繁窗口要短一些否则一次失锁会污染后面很长时间的平滑结果。判断失锁可以用相位观测的 LLI 标志或者检测相邻历元相位 TEC 变化量是否超过阈值。4. 斜向 TEC 到垂直 TEC 的映射与精度验证4.1 投影函数的选择与参数设置从观测文件直接算出来的是斜向 TECSTEC也就是信号传播路径上的积分。要做电离层建模或者和别的测站对比需要转成垂直 TECVTEC。转换靠投影函数VTEC STEC / F(z)F(z) 是投影函数z 是卫星在测站处的高度角。最常用的是单层模型假设下的余弦投影F(z) cos(arcsin(R/(RH) · sin(z)))其中 R 是地球半径约 6371 kmH 是电离层单层高度一般取 350 km 到 450 km。H 的选择会影响 VTEC 的绝对值350 km 和 450 km 算出来的 VTEC 能差 5% 到 10%。做长期分析时 H 要固定不能这次用 350 下次用 450。import math def mapping_function(elev_rad, h_ion350.0, r_earth6371.0): 单层模型投影函数elev_rad 是高度角弧度 z math.pi/2 - elev_rad sin_z math.sin(z) arg r_earth / (r_earth h_ion) * sin_z # 防止数值越界 arg max(-1.0, min(1.0, arg)) return math.cos(math.asin(arg)) def stec_to_vtec(stec, elev_rad, h_ion350.0): return stec * mapping_function(elev_rad, h_ion)高度角从测站坐标和卫星坐标算。卫星坐标可以从广播星历或者精密星历里取。如果你只做单站分析用广播星历算高度角精度就够误差在 0.1 度以内对投影函数的影响可以忽略。4.2 用已知 TEC 模型做交叉验证算出来的 VTEC 对不对最直接的验证方法是和公开的电离层模型对比。常见的有 IRI 模型、CODE 的 GIM 产品。对比时要注意时间系统和坐标系统要统一GIM 产品一般是 UTC 时间、地磁纬度你的 VTEC 要转过去才能比。我一般会挑几个电离层平静日做对比看 VTEC 的日变化曲线形状对不对。平静日 VTEC 应该在地方时正午前后达到峰值夜间降到最低峰值在 20 到 80 TECU 之间取决于纬度和季节。如果你的曲线峰值出现在半夜或者峰值超过 150 TECU那大概率是 DCB 没扣对或者投影函数用反了。# 简单的日变化检查找峰值对应的地方时 def check_diurnal_pattern(vtec_series, times_utc, lon_deg): 检查 VTEC 峰值是否在地方时正午附近 local_times [(t lon_deg/15.0) % 24 for t in times_utc] peak_idx np.argmax(vtec_series) peak_local local_times[peak_idx] peak_val vtec_series[peak_idx] print(f峰值 VTEC: {peak_val:.1f} TECU, 地方时: {peak_local:.1f}h) if not (10 peak_local 14): print(警告峰值不在正午附近检查 DCB 和投影函数) return peak_local, peak_val这个检查看起来简单但能抓住大部分低级错误。我见过有人把 L1 和 L2 的观测码搞反了算出来的 TEC 是负的日变化完全反过来了。也见过 DCB 符号搞反的VTEC 整体偏了 20 多个 TECU。4.3 多站对比与精度评估指标单站验证只能看趋势要做定量评估需要多站对比。选几个纬度不同的测站同一天算 VTEC看同一时刻不同纬度的 VTEC 是否符合电离层赤道异常的特征——赤道附近峰值高中纬度低高纬度更低。如果赤道站的 VTEC 反而比中纬度低那说明某个站的 DCB 或者坐标有问题。定量指标可以用均方根误差RMS和偏差Bias。和 GIM 对比时RMS 在 2 到 5 TECU 之间算正常Bias 在 ±2 TECU 以内可以接受。如果 RMS 超过 10 TECU先检查 DCB再检查投影函数最后检查高度角截止——低于 15 度的高度角数据噪声大一般要剔除。def evaluate_vtec(vtec_calc, vtec_ref): 计算 RMS 和 Bias diff np.array(vtec_calc) - np.array(vtec_ref) rms np.sqrt(np.mean(diff**2)) bias np.mean(diff) print(fRMS: {rms:.2f} TECU, Bias: {bias:.2f} TECU) return rms, bias高度角截止我一般设 20 度。低于 20 度时投影函数对高度角误差很敏感而且多路径效应严重伪距噪声能到几米算出来的 TEC 波动很大。设 20 度会损失一些低高度角的数据但换来的精度提升值得。5. 避坑与排查算 TEC 时最容易翻车的五个地方5.1 现象TEC 出现负值或日变化反转原因L1 和 L2 的观测码搞反了或者伪距组合的符号写反了。GPS 的 L1 频率高于 L2电离层延迟在 L1 上小于 L2所以 P2 - P1 应该是正值。如果算出来是负的先检查观测码映射。解决打印前几个历元的 P1、P2 原始值确认 P2 大于 P1。如果 P2 小于 P1检查 RINEX 头里的观测码顺序确认你取的C1C和C2W没有对调。另外注意有些接收机用C2C而不是C2W取错了会拿到空值。5.2 现象VTEC 整体偏高或偏低几十个 TECU原因DCB 没扣或者 DCB 单位搞错了。卫星 DCB 产品一般给纳秒要乘以光速转米。接收机 DCB 如果用了错误的估计值会引入系统性偏差。解决先不扣 DCB 算一遍再扣 DCB 算一遍看差值是否在合理范围一般 5 到 20 TECU。如果差值超过 50 TECU检查 DCB 的单位和符号。接收机 DCB 可以用平静日的夜间 VTEC 应该接近零这个先验来估计——夜间电离层电子密度低VTEC 一般在 5 TECU 以下。5.3 现象相位平滑后 TEC 出现台阶状跳变原因载波相位失锁导致整周跳变相位 TEC 的变化量里混入了周跳平滑结果被污染。解决在平滑之前做周跳检测。简单的方法是看相邻历元相位 TEC 变化量是否超过阈值比如 0.5 TECU/s超过就标记为周跳重置平滑窗口。更稳的方法是用 LLI 标志但有些文件里 LLI 不可靠还是要靠数值检测。def detect_cycle_slip(tec_phi, threshold0.5): 检测周跳返回布尔序列 slips np.zeros(len(tec_phi), dtypebool) for k in range(1, len(tec_phi)): if abs(tec_phi[k] - tec_phi[k-1]) threshold: slips[k] True return slips5.4 现象低高度角数据波动剧烈拉高整体 RMS原因低高度角时多路径效应和大气折射误差大伪距噪声显著增加投影函数对高度角误差的放大作用也更强。解决设高度角截止一般 15 到 20 度。截止之后数据量会减少但精度提升明显。如果必须用低高度角数据可以加大平滑窗口或者用相位 TEC 为主、伪距 TEC 只做绝对基准。5.5 现象不同测站算出的 VTEC 在同一时刻差异很大原因测站坐标错误、时间系统不统一、或者某个站的接收机 DCB 估计有误。解决先检查测站坐标用文件头里的近似坐标算高度角和精密星历对比。再检查时间标签RINEX 2.11 的时间是 GPS 时还是 UTC 要看文件头里的TIME SYSTEM标识。最后做多站对比时把每个站的 VTEC 减去自身夜间均值看相对变化是否一致这样能排除 DCB 的影响。6. 把 TEC 算准之后还能做什么从单站到区域建模的进阶路径单站 TEC 算通之后下一步自然是多站组网。把同一区域多个测站的 VTEC 放在一起可以做电离层梯度分析、行波检测甚至反演电子密度剖面。但这一步的前提是每个站的 DCB 都扣干净了否则站间的系统性偏差会淹没真实的电离层信号。我自己的习惯是每次拿到新数据先跑一遍单站日变化检查确认峰值在正午、夜间接近零、没有负值。这三条过了再进入多站流程。多站的时候我会先做站间交叉验证选两颗共同可见的卫星比较两个站在同一时刻、同一卫星的 VTEC差异应该在 2 TECU 以内。如果超过 5 TECU先查接收机 DCB再查坐标。区域建模常用的是多项式拟合或者球谐函数。多项式拟合简单适合小区域几百公里把 VTEC 表示为纬度和经度的二次多项式用最小二乘解系数。球谐函数适合大区域甚至全球但需要的数据量大单站数据做不了。我一般从多项式开始阶数取 2 到 3阶数太高会过拟合阶数太低拟合残差大。def fit_vtec_poly(lat, lon, vtec, order2): 二次多项式拟合 VTEC lat/lon 单位度 from numpy.linalg import lstsq # 构造设计矩阵 cols [] for i in range(order1): for j in range(order1-i): cols.append((lat**i) * (lon**j)) A np.column_stack(cols) coef, _, _, _ lstsq(A, vtec, rcondNone) return coef def predict_vtec_poly(coef, lat, lon, order2): cols [] for i in range(order1): for j in range(order1-i): cols.append((lat**i) * (lon**j)) return np.dot(np.array(cols), coef)拟合完要看残差分布。如果残差在空间上有明显结构比如某个方向系统性偏大说明阶数不够或者该区域有电离层梯度没被模型抓住。这时候要么升阶要么改用局部插值方法比如克里金。还有一个进阶方向是结合精密星历做更精确的高度角计算。广播星历的卫星位置误差在米级对高度角的影响在 0.01 度量级对投影函数的影响可以忽略。但如果你要做实时 TEC 或者高精度应用精密星历的延迟和精度优势就体现出来了。最后说一个我踩过的坑TEC 的单位。TECU 是 10¹⁶ el/m²但有些文献用 10¹⁶ el/m² 有些用 10¹⁷公式里的常数会差一个量级。我一般会在代码里把常数定义成模块级变量算完之后打印几个典型值看看量级对不对——平静日 VTEC 在 10 到 50 TECU电离层暴时能到 100 以上。如果算出来是 0.1 或者 1000先查常数。这套流程从 RINEX 解析到 VTEC 验证代码量不大但每一步都有容易翻车的地方。我的习惯是每加一个改正项就打印一次中间结果和上一步对比看变化量是否合理。这样出问题的时候能快速定位到是哪一步引入的。希望帮到你。本文还有配套的精品资源点击获取