简介针对RINEX 2.11格式的双频GPS观测数据这份MATLAB工具包实现了总电子含量TEC的完整计算流程适合从事GNSS电离层研究、大地测量及相关课程设计的工程师与学生。程序内置ProcessTECCalculation.m主脚本支持输出垂直TEC、倾斜TEC、带接收机/卫星DCB的STEC以及ROTI指数并包含周跳修正、DCB文件解析、卫星位置计算等配套函数。资源共84个文件以46个DCB差分码偏差文件、21个m格式MATLAB源码、可执行工具及mexw32/64编译组件为主压缩包大小约28.89MB结构清晰便于直接调用。目前已有873人浏览学习适合需要快速上手TEC计算与电离层监测分析的MATLAB开发者。1. 在 MATLAB 里从 RINEX 2.11 算 TEC为什么值得自己写一遍GNSS 数据处理里总电子含量 TEC 是最容易拿到、又最容易被算错的一个物理量。手里有一台双频接收机导出的原始观测是 RINEX 2.11 格式想在 MATLAB 里按历元算出每颗卫星方向的 TEC第一反应往往是去找现成工具箱。但 MATLAB 自带的 RINEX 读取函数对 2.x 老格式的支持并不完整遇到删帧、非标空白、观测值跨行时会静默出错。自己按 RINEX 2.11 标准写解析再把双频伪距和载波相位组合成无几何量两百行左右就能在 MATLAB 里跑通完整流程。这里说的 TEC是电离层总电子含量 Total Electron Content检索时注意别和热电制冷器的 TEC 混在一起。这条路径适合做电离层研究、GNSS 数据后处理以及刚接触精密测量的工程师整条链路从文件到 TECU 曲线都是可控的。2. RINEX 2.11 观测文件的结构与 MATLAB 解析实现2.1 头文件段里决定 TEC 计算成败的观测类型字段RINEX 2.11 观测文件由头文件段和数据段构成头文件段以END OF HEADER结束。与 TEC 计算强相关的头文件记录是# / TYPES OF OBSERV它按顺序列出数据段中每颗卫星的观测量类型。这个顺序一旦读错后续列映射全部错位而且程序不会报错只会输出一套看起来合理、实际完全不可用的 TEC。双频接收机最常见的输出类型是 C1、P1、P2、L1、L2。计算 TEC 使用的是无几何组合要求两个频率的伪距和相位同时存在。C1 与 P1 的区别在码类型消费级接收机很多只输出 C1 和 P2此时用 C1-P2 组合同样能算 TEC只是常数偏置与 P1-P2 组合不同。下面是这些观测类型在 TEC 流程中的角色。观测类型含义单位TEC 计算中的角色C1L1 频率 C/A 码伪距米与 P2 组成无几何伪距组合P1L1 频率 P 码伪距米C1 缺失时替代 C1P2L2 频率 P 码伪距米无几何组合的低频端L1L1 载波相位周相位平滑与周跳检测L2L2 载波相位周相位平滑与周跳检测S1 / S2L1 / L2 信噪比dB-Hz数据剔除与权重分配不值得在解析阶段花太多精力处理 S1/S2TEC 精度主要由伪距和相位决定。信噪比可以在后续质量控制里做阈值筛选但大多数单站 TEC 任务不会用到它。2.2 用固定列宽逐行解析的 read_rinex211_obs 函数RINEX 2.11 的历元行和数据行是固定列宽格式不能用空格切割因为观测值字段之间允许出现可变数量的空格。常见做法是按字符位置切片这也是下面这段代码选择的方式。function obs read_rinex211_obs(filepath) % read_rinex211_obs: RINEX 2.11 双频观测文件读取 % 输出 obs 为结构体数组, 字段: epoch, prn, C1, P2, L1, L2 % 缺失观测量置 NaN, 不满足双频条件的历元自动剔除 fid fopen(filepath, r); if fid -1 error(无法打开文件: %s, filepath); end % ---- 头文件段解析 ---- obs_types {}; while true line fgetl(fid); if ~ischar(line) error(文件缺少 END OF HEADER); end if numel(line) 61 label strtrim(line(61:end)); if strcmp(label, END OF HEADER) break; end if strcmp(label, # / TYPES OF OBSERV) obs_types regexp(strtrim(line(11:60)), \S, match); end end end % 确定所需观测类型在数据列中的位置 need {C1, P2, L1, L2}; col zeros(1, 4); for k 1:4 idx find(strcmp(obs_types, need{k}), 1); if ~isempty(idx) col(k) idx; elseif strcmp(need{k}, C1) % C1 缺失时退回 P1 idx find(strcmp(obs_types, P1), 1); if ~isempty(idx), col(k) idx; end end end if any(col 0) error(观测文件缺少双频必需类型: C1/P1, P2, L1, L2); end % ---- 数据段历元循环 ---- obs struct(epoch, {}, prn, {}, C1, {}, P2, {}, L1, {}, L2, {}); n 0; n_obs_types numel(obs_types); while true line fgetl(fid); if ~ischar(line), break; end if numel(line) 32, continue; end % 历元行固定宽度: 2位年, 2位月, 2位日, 2位时, 2位分, 11位秒, 标志, 卫星数 yy str2double(line(2:3)); mo str2double(line(5:6)); dd str2double(line(8:9)); hh str2double(line(11:12)); mi str2double(line(14:15)); ss str2double(line(16:26)); flag str2double(line(29)); nsat str2double(line(30:32)); if ~isnan(flag) flag ~ 0 flag ~ 1 continue; % 事件历元不参与 TEC 计算 end % 卫星列表: 从第 33 列起每 3 字符一颗卫星 sats cell(1, nsat); for k 1:nsat st 32 3 * (k - 1) 1; s strtrim(line(st:min(st 2, numel(line)))); if numel(s) 2 % GPS 卫星未带系统标识时补 G s [G, s]; end sats{k} s; end % 逐颗卫星读取观测值 for k 1:nsat sat_line fgetl(fid); if ~ischar(sat_line), break; end vals nan(1, n_obs_types); % 一行最多 5 个观测值, 超过部分继续读后续行 for row_idx 0:ceil(n_obs_types / 5) - 1 if row_idx 0 sat_line fgetl(fid); if ~ischar(sat_line), break; end end for j 1:min(5, n_obs_types - row_idx * 5) start_col 3 14 * (row_idx * 5 j - 1) 1; val str2double(sat_line(start_col:min(start_col 13, numel(sat_line)))); if val 0, val nan; end vals(row_idx * 5 j) val; end end c1 vals(col(1)); p2 vals(col(2)); l1 vals(col(3)); l2 vals(col(4)); if all(~isnan([c1, p2, l1, l2])) n n 1; obs(n).epoch datetime(yy 2000, mo, dd, hh, mi, ss); obs(n).prn sats{k}; obs(n).C1 c1; obs(n).P2 p2; obs(n).L1 l1; obs(n).L2 l2; end end end fclose(fid); end这段代码的逻辑重点是固定列宽切片和numel(line) 32的防御性判断。RINEX 2.11 中有些老文件的历元行会以空格结尾直接按空格拆分会把时间字段和卫星数拆散所以按位置切片是唯一可靠的方式。line(61:end)解析头文件标签RINEX 格式规定标签从第 61 列开始但实际文件可能有截断因此先判numel(line) 61。观测类型超过 5 个时同一颗卫星的观测值会跨行内层循环按 5 个一组继续读后续行。这里最容易出错的点是跨行后的续行没有卫星号但前 3 列仍然保留读取起始列仍要从第 4 列开始。2.3 解析后转成时间序列的高效做法解析出的结构体数组通常要按卫星分组转成时间序列方便后续按弧段处理。prn_list unique({obs.prn}); gps_prn prn_list(startsWith(prn_list, G)); first_prn gps_prn{1}; sel strcmp({obs.prn}, first_prn); t [obs(sel).epoch]; C1 [obs(sel).C1]; P2 [obs(sel).P2]; L1 [obs(sel).L1]; L2 [obs(sel).L2];卫星命名统一成 G 加两位 PRN 后与导航星历文件、精密星历文件的匹配会省去很多麻烦。first_prn只是演示实际处理时应该把所有可见卫星都循环一遍。这里要特别提醒不要把 RINEX 数据行交给 AI 辅助工具自动“优化”成按空格拆分比如用 codex 改写后很容易引入strsplit(line)这类实现遇到观测值缺失的空白字段时整行解析就断了。3. TEC 解算公式双频伪距无几何组合的推导与 MATLAB 实现3.1 为什么双频组合能抵消全部几何项电离层对微波信号是色散介质。伪距测距时信号穿过电离层产生与频率平方成反比的附加延迟载波相位则出现等量的相位提前。对同一颗卫星同一时刻伪距观测量里除了几何距离、钟差、对流层延迟剩下的频率相关项只有电离层延迟和硬件延迟。把 L1 和 L2 两个伪距做差几何距离、卫星钟差、接收机钟差、对流层延迟全部消掉剩下的就是电离层项的差。这是双频 TEC 计算能成立的根本原因。单频接收机必须依赖外部电离层模型或者格网产品本质上是在猜 TEC双频接收机则是直接测量 TEC这也是标题里强调双频接收器的意义。伪距无几何组合的基本关系是P2 - P1 40.3 × TEC × (1/f2² - 1/f1²)其中 TEC 是斜路径上的电子总量单位是电子/平方米P2 与 P1 的单位是米f1 与 f2 是信号频率。注意这个公式没有考虑硬件延迟所以真实数据里 TEC 结果会带一个常数偏置第 5 章会讲如何处理。3.2 简化的比例系数 9.524 是怎么来的GPS 的 L1 频率为 1575.42 MHzL2 频率为 1227.60 MHz代入后整理比例系数TEC (P2 - P1) × f1²f2² / (40.3 × (f1² - f2²))把频率值代入除以 1e16 把单位换成 TECU得到 TECU 9.524 × (P2 - P1)。这个系数在文献里经常直接给出没有推导过程。自己算一遍的好处是能确认符号。如果某篇参考资料用的是 P1-P2系数就要取负号。实际计算时还有一个容易忽略的单位问题。RINEX 2.11 文件里 L1、L2 的单位是周必须先乘波长转换成米才能与伪距组合放在同一个尺度上。3.3 直接计算的 MATLAB 代码与符号约定% 双频 TEC 计算: 伪距无几何组合 % C1, P2 单位: 米; TEC_raw 单位: TECU TEC_raw 9.524 * (P2 - C1); % 载波相位无几何组合, 单位: 米 lambda1 299792458 / 1575.42e6; lambda2 299792458 / 1227.60e6; L_geom L1 * lambda1 - L2 * lambda2; % 相位 TEC 差分变化, 用于周跳检测和后继平滑 dTEC_phase 9.524 * diff(L_geom);这里的符号约定是 P2-P1 为正时 TEC 为正。L_geom表示 L1 相位距离减去 L2 相位距离其模糊度项是常数差分后只剩 TEC 变化和噪声。diff(L_geom)如果出现明显跳变就是周跳信号后面做平滑时需要在这里断开弧段。直接算出来的TEC_raw噪声很大。伪距噪声通常在分米到米级对应 TEC 误差约为几个到几十个 TECU。把TEC_raw和cumsum(dTEC_phase)画在同一张图中两条曲线形状应该一致只是相位版本带常数偏置。如果形状对不上优先检查头文件里# / TYPES OF OBSERV的列顺序是否解析正确。4. 提高 TEC 精度的关键处理相位平滑、周跳检测与 VTEC 投影4.1 伪距噪声与相位模糊度的互补关系TEC_raw噪声大直接用于电离层建模不够。载波相位观测噪声只有毫米到厘米级其无几何组合能精确刻画 TEC 的短时变化趋势但含有未知模糊度常数。把两者结合就是经典的分段平滑思路伪距 TEC 决定常量水平相位 TEC 决定变化细节。另外双频接收机算出来的 TEC 是从接收机到卫星的斜路径电子总量。做单站电离层监测时通常还要投影到垂直方向得到 VTEC。投影需要卫星仰角而 RINEX 观测文件里没有仰角需要结合接收机概略坐标和卫星位置计算。接收机概略坐标可以直接读 RINEX 头文件里的APPROX POSITION XYZ卫星位置则要由导航文件星历计算这部分本身是另一个主题这里假设仰角序列已经通过星历解算得到记为elev_deg。4.2 基于 Hatch 滤波的平滑实现对每一颗卫星的连续弧段使用如下递推公式PS(k) w × P_raw(k) (1 - w) × [PS(k-1) L_geom(k) - L_geom(k-1)]其中 P_raw(k) P2(k) - C1(k)L_geom(k) L1(k)×λ1 - L2(k)×λ2w 是平滑权重取值在 0.01 到 0.05 之间。w 越大平滑越弱w 越小平滑越强但对周跳越敏感。function [TEC_smooth, tec_seg] smooth_tec(C1, P2, L1, L2, w, cycle_threshold) % smooth_tec: 载波相位平滑伪距的无几何组合 % C1,P2,L1,L2: 单颗卫星的时间序列 % w: 平滑权重; cycle_threshold: 周跳阈值(米) lambda1 299792458 / 1575.42e6; lambda2 299792458 / 1227.60e6; L_geom L1 * lambda1 - L2 * lambda2; P_raw P2 - C1; n numel(P_raw); TEC_smooth nan(1, n); tec_seg zeros(1, n); % 弧段编号 seg_id 0; ps_prev nan; for k 1:n if isnan(P_raw(k)) || isnan(L_geom(k)) ps_prev nan; % 数据中断, 重置滤波 continue; end if k 1 abs(L_geom(k) - L_geom(k-1)) cycle_threshold ps_prev nan; % 周跳, 重置弧段 end if isnan(ps_prev) ps P_raw(k); % 新弧段用伪距初始化 seg_id seg_id 1; else ps w * P_raw(k) (1 - w) * (ps_prev L_geom(k) - L_geom(k-1)); end ps_prev ps; TEC_smooth(k) ps * 9.524; tec_seg(k) seg_id; end endps_prev置为nan的两种情况分别是数据缺失和周跳它们在物理上都代表信号链路被打断模糊度常数发生随机跳变必须重置弧段。tec_seg用来标记弧段编号后续按弧段统计 TEC 均值时有用。周跳阈值的选取要结合环境。L_geom 这个组合的噪声在厘米级周跳在该组合上至少产生几厘米到几十厘米的跳变。静态测量场景取 0.05 米比较合适车载等动态环境卫星信号遮挡频繁可放宽到 0.10 米否则频繁误判周跳会导致弧段过短、平滑失效。4.3 薄层电离层模型下的 VTEC 投影垂直投影采用薄层电离层模型假设自由电子集中在一个距地面高度 H 的薄球壳上。设接收机处卫星仰角为 E穿刺点天顶角 χ 满足sin(χ) Re / (Re H) × sin(90° - E)其中 Re 取 6371 kmH 取 350 km 或 450 kmVTEC STEC × cos(χ)。function VTEC stec_to_vtec(STEC, elev_deg) % stec_to_vtec: 斜路径 TEC 转垂直 TEC % STEC: 斜路径 TEC (TECU); elev_deg: 卫星仰角 (度) Re 6371.0; H 350.0; % 电离层薄层高度, 常用 350 或 450 km chi deg2rad(90 - elev_deg); sin_chi_p Re / (Re H) * sin(chi); cos_chi_p sqrt(1 - sin_chi_p^2); VTEC STEC * cos_chi_p; end这个映射函数在低仰角时放大倍数很大因此必须先做截止角筛选。单站 TEC 计算通常取 10° 到 15° 截止角低于截止角的数据投影误差和伪距噪声都太大算出来的 VTEC 没有使用价值。H 取 350 还是 450 km对高仰角卫星几乎没有差别主要影响低仰角弧段的投影幅度建议固定一个值并写进处理日志。参数推荐值影响平滑权重 w0.01~0.05越小越平滑但收敛越慢周跳阈值0.05~0.10 米过小误判噪声为周跳薄层高度 H350 或 450 km影响低仰角投影幅度截止角10°~15°低于此值数据噪声大5. 验证 TEC 结果是否可信的三个实用技巧5.1 相位与伪距 TEC 的斜率一致性检查TEC 计算代码跑通后先不要急着看绝对数值先检查形状。对同一颗卫星把伪距 TEC 与相位差分累积 TEC 画在同一张图里两条曲线的斜率在连续弧段内应该一致。如果斜率不一致最可能的原因是观测列错位比如把 L1 和 L2 的顺序调换了或者 C1 和 P2 的符号被写反。figure; plot(t, TEC_raw, ., MarkerSize, 4); hold on; plot(t, TEC_smooth, -, LineWidth, 1.2); legend(伪距 TEC, 平滑 TEC); ylabel(TECU);这一步还能检查周跳标记是否合理。正常情况下平滑后的 TEC 是一条连续曲线如果出现大量锯齿状分段说明cycle_threshold设得太小把观测噪声误判成了周跳。5.2 与 IGS GIM 网格 TEC 对比的边界条件与 IGS 的全球电离层格网产品对比是最常用的外部验证方式。IGS GIM 以 IONEX 格式发布空间分辨率约 2.5°×5°时间分辨率 1 小时。对比时的常见误区是直接拿单站 30 秒采样 TEC 与 GIM 逐历元比较这没有意义因为 GIM 本身抹平了短时电离层变化。正确做法是先对本地 VTEC 做 10 到 15 分钟的滑动平均再插值到 GIM 对应时刻和穿刺点位置。中纬度白天两者相差 3 到 10 TECU 属于正常范围太阳活动高年会更大。需要明确一点GIM 也是模型值不是真值它更适合用来检查量级和趋势不能作为逐历元真值。5.3 输出带弧段号的 CSV 并检查 DCB 痕迹多颗卫星的 TEC 序列放在一起时不同弧段之间会因为卫星和接收机硬件延迟出现常数阶梯差。这个差值就是 DCB 的表现。单站双频处理无法单独分离接收机与卫星 DCB但可以把它当作常数偏置处理。如果发现某条弧段整体比其他卫星高几十 TECU不一定是计算错误先确认是不是 DCB 导致的。落盘时把弧段号保存下来后续分析会省很多事。result table(t, repmat(first_prn, numel(t), 1), TEC_raw, TEC_smooth, tec_seg, ... VariableNames, {epoch, prn, TEC_raw, TEC_smooth, seg_id}); writetable(result, tec_result.csv);文件里保留seg_id字段后按弧段统计 TEC 均值、筛选连续观测时长、做 DCB 估计都无需再对周跳做二次判断。这也是把处理流程固化成脚本后最低成本保存完整语义的输出方式。本文还有配套的精品资源点击获取
MATLAB实现RINEX 2.11双频观测TEC解算与验证
简介针对RINEX 2.11格式的双频GPS观测数据这份MATLAB工具包实现了总电子含量TEC的完整计算流程适合从事GNSS电离层研究、大地测量及相关课程设计的工程师与学生。程序内置ProcessTECCalculation.m主脚本支持输出垂直TEC、倾斜TEC、带接收机/卫星DCB的STEC以及ROTI指数并包含周跳修正、DCB文件解析、卫星位置计算等配套函数。资源共84个文件以46个DCB差分码偏差文件、21个m格式MATLAB源码、可执行工具及mexw32/64编译组件为主压缩包大小约28.89MB结构清晰便于直接调用。目前已有873人浏览学习适合需要快速上手TEC计算与电离层监测分析的MATLAB开发者。1. 在 MATLAB 里从 RINEX 2.11 算 TEC为什么值得自己写一遍GNSS 数据处理里总电子含量 TEC 是最容易拿到、又最容易被算错的一个物理量。手里有一台双频接收机导出的原始观测是 RINEX 2.11 格式想在 MATLAB 里按历元算出每颗卫星方向的 TEC第一反应往往是去找现成工具箱。但 MATLAB 自带的 RINEX 读取函数对 2.x 老格式的支持并不完整遇到删帧、非标空白、观测值跨行时会静默出错。自己按 RINEX 2.11 标准写解析再把双频伪距和载波相位组合成无几何量两百行左右就能在 MATLAB 里跑通完整流程。这里说的 TEC是电离层总电子含量 Total Electron Content检索时注意别和热电制冷器的 TEC 混在一起。这条路径适合做电离层研究、GNSS 数据后处理以及刚接触精密测量的工程师整条链路从文件到 TECU 曲线都是可控的。2. RINEX 2.11 观测文件的结构与 MATLAB 解析实现2.1 头文件段里决定 TEC 计算成败的观测类型字段RINEX 2.11 观测文件由头文件段和数据段构成头文件段以END OF HEADER结束。与 TEC 计算强相关的头文件记录是# / TYPES OF OBSERV它按顺序列出数据段中每颗卫星的观测量类型。这个顺序一旦读错后续列映射全部错位而且程序不会报错只会输出一套看起来合理、实际完全不可用的 TEC。双频接收机最常见的输出类型是 C1、P1、P2、L1、L2。计算 TEC 使用的是无几何组合要求两个频率的伪距和相位同时存在。C1 与 P1 的区别在码类型消费级接收机很多只输出 C1 和 P2此时用 C1-P2 组合同样能算 TEC只是常数偏置与 P1-P2 组合不同。下面是这些观测类型在 TEC 流程中的角色。观测类型含义单位TEC 计算中的角色C1L1 频率 C/A 码伪距米与 P2 组成无几何伪距组合P1L1 频率 P 码伪距米C1 缺失时替代 C1P2L2 频率 P 码伪距米无几何组合的低频端L1L1 载波相位周相位平滑与周跳检测L2L2 载波相位周相位平滑与周跳检测S1 / S2L1 / L2 信噪比dB-Hz数据剔除与权重分配不值得在解析阶段花太多精力处理 S1/S2TEC 精度主要由伪距和相位决定。信噪比可以在后续质量控制里做阈值筛选但大多数单站 TEC 任务不会用到它。2.2 用固定列宽逐行解析的 read_rinex211_obs 函数RINEX 2.11 的历元行和数据行是固定列宽格式不能用空格切割因为观测值字段之间允许出现可变数量的空格。常见做法是按字符位置切片这也是下面这段代码选择的方式。function obs read_rinex211_obs(filepath) % read_rinex211_obs: RINEX 2.11 双频观测文件读取 % 输出 obs 为结构体数组, 字段: epoch, prn, C1, P2, L1, L2 % 缺失观测量置 NaN, 不满足双频条件的历元自动剔除 fid fopen(filepath, r); if fid -1 error(无法打开文件: %s, filepath); end % ---- 头文件段解析 ---- obs_types {}; while true line fgetl(fid); if ~ischar(line) error(文件缺少 END OF HEADER); end if numel(line) 61 label strtrim(line(61:end)); if strcmp(label, END OF HEADER) break; end if strcmp(label, # / TYPES OF OBSERV) obs_types regexp(strtrim(line(11:60)), \S, match); end end end % 确定所需观测类型在数据列中的位置 need {C1, P2, L1, L2}; col zeros(1, 4); for k 1:4 idx find(strcmp(obs_types, need{k}), 1); if ~isempty(idx) col(k) idx; elseif strcmp(need{k}, C1) % C1 缺失时退回 P1 idx find(strcmp(obs_types, P1), 1); if ~isempty(idx), col(k) idx; end end end if any(col 0) error(观测文件缺少双频必需类型: C1/P1, P2, L1, L2); end % ---- 数据段历元循环 ---- obs struct(epoch, {}, prn, {}, C1, {}, P2, {}, L1, {}, L2, {}); n 0; n_obs_types numel(obs_types); while true line fgetl(fid); if ~ischar(line), break; end if numel(line) 32, continue; end % 历元行固定宽度: 2位年, 2位月, 2位日, 2位时, 2位分, 11位秒, 标志, 卫星数 yy str2double(line(2:3)); mo str2double(line(5:6)); dd str2double(line(8:9)); hh str2double(line(11:12)); mi str2double(line(14:15)); ss str2double(line(16:26)); flag str2double(line(29)); nsat str2double(line(30:32)); if ~isnan(flag) flag ~ 0 flag ~ 1 continue; % 事件历元不参与 TEC 计算 end % 卫星列表: 从第 33 列起每 3 字符一颗卫星 sats cell(1, nsat); for k 1:nsat st 32 3 * (k - 1) 1; s strtrim(line(st:min(st 2, numel(line)))); if numel(s) 2 % GPS 卫星未带系统标识时补 G s [G, s]; end sats{k} s; end % 逐颗卫星读取观测值 for k 1:nsat sat_line fgetl(fid); if ~ischar(sat_line), break; end vals nan(1, n_obs_types); % 一行最多 5 个观测值, 超过部分继续读后续行 for row_idx 0:ceil(n_obs_types / 5) - 1 if row_idx 0 sat_line fgetl(fid); if ~ischar(sat_line), break; end end for j 1:min(5, n_obs_types - row_idx * 5) start_col 3 14 * (row_idx * 5 j - 1) 1; val str2double(sat_line(start_col:min(start_col 13, numel(sat_line)))); if val 0, val nan; end vals(row_idx * 5 j) val; end end c1 vals(col(1)); p2 vals(col(2)); l1 vals(col(3)); l2 vals(col(4)); if all(~isnan([c1, p2, l1, l2])) n n 1; obs(n).epoch datetime(yy 2000, mo, dd, hh, mi, ss); obs(n).prn sats{k}; obs(n).C1 c1; obs(n).P2 p2; obs(n).L1 l1; obs(n).L2 l2; end end end fclose(fid); end这段代码的逻辑重点是固定列宽切片和numel(line) 32的防御性判断。RINEX 2.11 中有些老文件的历元行会以空格结尾直接按空格拆分会把时间字段和卫星数拆散所以按位置切片是唯一可靠的方式。line(61:end)解析头文件标签RINEX 格式规定标签从第 61 列开始但实际文件可能有截断因此先判numel(line) 61。观测类型超过 5 个时同一颗卫星的观测值会跨行内层循环按 5 个一组继续读后续行。这里最容易出错的点是跨行后的续行没有卫星号但前 3 列仍然保留读取起始列仍要从第 4 列开始。2.3 解析后转成时间序列的高效做法解析出的结构体数组通常要按卫星分组转成时间序列方便后续按弧段处理。prn_list unique({obs.prn}); gps_prn prn_list(startsWith(prn_list, G)); first_prn gps_prn{1}; sel strcmp({obs.prn}, first_prn); t [obs(sel).epoch]; C1 [obs(sel).C1]; P2 [obs(sel).P2]; L1 [obs(sel).L1]; L2 [obs(sel).L2];卫星命名统一成 G 加两位 PRN 后与导航星历文件、精密星历文件的匹配会省去很多麻烦。first_prn只是演示实际处理时应该把所有可见卫星都循环一遍。这里要特别提醒不要把 RINEX 数据行交给 AI 辅助工具自动“优化”成按空格拆分比如用 codex 改写后很容易引入strsplit(line)这类实现遇到观测值缺失的空白字段时整行解析就断了。3. TEC 解算公式双频伪距无几何组合的推导与 MATLAB 实现3.1 为什么双频组合能抵消全部几何项电离层对微波信号是色散介质。伪距测距时信号穿过电离层产生与频率平方成反比的附加延迟载波相位则出现等量的相位提前。对同一颗卫星同一时刻伪距观测量里除了几何距离、钟差、对流层延迟剩下的频率相关项只有电离层延迟和硬件延迟。把 L1 和 L2 两个伪距做差几何距离、卫星钟差、接收机钟差、对流层延迟全部消掉剩下的就是电离层项的差。这是双频 TEC 计算能成立的根本原因。单频接收机必须依赖外部电离层模型或者格网产品本质上是在猜 TEC双频接收机则是直接测量 TEC这也是标题里强调双频接收器的意义。伪距无几何组合的基本关系是P2 - P1 40.3 × TEC × (1/f2² - 1/f1²)其中 TEC 是斜路径上的电子总量单位是电子/平方米P2 与 P1 的单位是米f1 与 f2 是信号频率。注意这个公式没有考虑硬件延迟所以真实数据里 TEC 结果会带一个常数偏置第 5 章会讲如何处理。3.2 简化的比例系数 9.524 是怎么来的GPS 的 L1 频率为 1575.42 MHzL2 频率为 1227.60 MHz代入后整理比例系数TEC (P2 - P1) × f1²f2² / (40.3 × (f1² - f2²))把频率值代入除以 1e16 把单位换成 TECU得到 TECU 9.524 × (P2 - P1)。这个系数在文献里经常直接给出没有推导过程。自己算一遍的好处是能确认符号。如果某篇参考资料用的是 P1-P2系数就要取负号。实际计算时还有一个容易忽略的单位问题。RINEX 2.11 文件里 L1、L2 的单位是周必须先乘波长转换成米才能与伪距组合放在同一个尺度上。3.3 直接计算的 MATLAB 代码与符号约定% 双频 TEC 计算: 伪距无几何组合 % C1, P2 单位: 米; TEC_raw 单位: TECU TEC_raw 9.524 * (P2 - C1); % 载波相位无几何组合, 单位: 米 lambda1 299792458 / 1575.42e6; lambda2 299792458 / 1227.60e6; L_geom L1 * lambda1 - L2 * lambda2; % 相位 TEC 差分变化, 用于周跳检测和后继平滑 dTEC_phase 9.524 * diff(L_geom);这里的符号约定是 P2-P1 为正时 TEC 为正。L_geom表示 L1 相位距离减去 L2 相位距离其模糊度项是常数差分后只剩 TEC 变化和噪声。diff(L_geom)如果出现明显跳变就是周跳信号后面做平滑时需要在这里断开弧段。直接算出来的TEC_raw噪声很大。伪距噪声通常在分米到米级对应 TEC 误差约为几个到几十个 TECU。把TEC_raw和cumsum(dTEC_phase)画在同一张图中两条曲线形状应该一致只是相位版本带常数偏置。如果形状对不上优先检查头文件里# / TYPES OF OBSERV的列顺序是否解析正确。4. 提高 TEC 精度的关键处理相位平滑、周跳检测与 VTEC 投影4.1 伪距噪声与相位模糊度的互补关系TEC_raw噪声大直接用于电离层建模不够。载波相位观测噪声只有毫米到厘米级其无几何组合能精确刻画 TEC 的短时变化趋势但含有未知模糊度常数。把两者结合就是经典的分段平滑思路伪距 TEC 决定常量水平相位 TEC 决定变化细节。另外双频接收机算出来的 TEC 是从接收机到卫星的斜路径电子总量。做单站电离层监测时通常还要投影到垂直方向得到 VTEC。投影需要卫星仰角而 RINEX 观测文件里没有仰角需要结合接收机概略坐标和卫星位置计算。接收机概略坐标可以直接读 RINEX 头文件里的APPROX POSITION XYZ卫星位置则要由导航文件星历计算这部分本身是另一个主题这里假设仰角序列已经通过星历解算得到记为elev_deg。4.2 基于 Hatch 滤波的平滑实现对每一颗卫星的连续弧段使用如下递推公式PS(k) w × P_raw(k) (1 - w) × [PS(k-1) L_geom(k) - L_geom(k-1)]其中 P_raw(k) P2(k) - C1(k)L_geom(k) L1(k)×λ1 - L2(k)×λ2w 是平滑权重取值在 0.01 到 0.05 之间。w 越大平滑越弱w 越小平滑越强但对周跳越敏感。function [TEC_smooth, tec_seg] smooth_tec(C1, P2, L1, L2, w, cycle_threshold) % smooth_tec: 载波相位平滑伪距的无几何组合 % C1,P2,L1,L2: 单颗卫星的时间序列 % w: 平滑权重; cycle_threshold: 周跳阈值(米) lambda1 299792458 / 1575.42e6; lambda2 299792458 / 1227.60e6; L_geom L1 * lambda1 - L2 * lambda2; P_raw P2 - C1; n numel(P_raw); TEC_smooth nan(1, n); tec_seg zeros(1, n); % 弧段编号 seg_id 0; ps_prev nan; for k 1:n if isnan(P_raw(k)) || isnan(L_geom(k)) ps_prev nan; % 数据中断, 重置滤波 continue; end if k 1 abs(L_geom(k) - L_geom(k-1)) cycle_threshold ps_prev nan; % 周跳, 重置弧段 end if isnan(ps_prev) ps P_raw(k); % 新弧段用伪距初始化 seg_id seg_id 1; else ps w * P_raw(k) (1 - w) * (ps_prev L_geom(k) - L_geom(k-1)); end ps_prev ps; TEC_smooth(k) ps * 9.524; tec_seg(k) seg_id; end endps_prev置为nan的两种情况分别是数据缺失和周跳它们在物理上都代表信号链路被打断模糊度常数发生随机跳变必须重置弧段。tec_seg用来标记弧段编号后续按弧段统计 TEC 均值时有用。周跳阈值的选取要结合环境。L_geom 这个组合的噪声在厘米级周跳在该组合上至少产生几厘米到几十厘米的跳变。静态测量场景取 0.05 米比较合适车载等动态环境卫星信号遮挡频繁可放宽到 0.10 米否则频繁误判周跳会导致弧段过短、平滑失效。4.3 薄层电离层模型下的 VTEC 投影垂直投影采用薄层电离层模型假设自由电子集中在一个距地面高度 H 的薄球壳上。设接收机处卫星仰角为 E穿刺点天顶角 χ 满足sin(χ) Re / (Re H) × sin(90° - E)其中 Re 取 6371 kmH 取 350 km 或 450 kmVTEC STEC × cos(χ)。function VTEC stec_to_vtec(STEC, elev_deg) % stec_to_vtec: 斜路径 TEC 转垂直 TEC % STEC: 斜路径 TEC (TECU); elev_deg: 卫星仰角 (度) Re 6371.0; H 350.0; % 电离层薄层高度, 常用 350 或 450 km chi deg2rad(90 - elev_deg); sin_chi_p Re / (Re H) * sin(chi); cos_chi_p sqrt(1 - sin_chi_p^2); VTEC STEC * cos_chi_p; end这个映射函数在低仰角时放大倍数很大因此必须先做截止角筛选。单站 TEC 计算通常取 10° 到 15° 截止角低于截止角的数据投影误差和伪距噪声都太大算出来的 VTEC 没有使用价值。H 取 350 还是 450 km对高仰角卫星几乎没有差别主要影响低仰角弧段的投影幅度建议固定一个值并写进处理日志。参数推荐值影响平滑权重 w0.01~0.05越小越平滑但收敛越慢周跳阈值0.05~0.10 米过小误判噪声为周跳薄层高度 H350 或 450 km影响低仰角投影幅度截止角10°~15°低于此值数据噪声大5. 验证 TEC 结果是否可信的三个实用技巧5.1 相位与伪距 TEC 的斜率一致性检查TEC 计算代码跑通后先不要急着看绝对数值先检查形状。对同一颗卫星把伪距 TEC 与相位差分累积 TEC 画在同一张图里两条曲线的斜率在连续弧段内应该一致。如果斜率不一致最可能的原因是观测列错位比如把 L1 和 L2 的顺序调换了或者 C1 和 P2 的符号被写反。figure; plot(t, TEC_raw, ., MarkerSize, 4); hold on; plot(t, TEC_smooth, -, LineWidth, 1.2); legend(伪距 TEC, 平滑 TEC); ylabel(TECU);这一步还能检查周跳标记是否合理。正常情况下平滑后的 TEC 是一条连续曲线如果出现大量锯齿状分段说明cycle_threshold设得太小把观测噪声误判成了周跳。5.2 与 IGS GIM 网格 TEC 对比的边界条件与 IGS 的全球电离层格网产品对比是最常用的外部验证方式。IGS GIM 以 IONEX 格式发布空间分辨率约 2.5°×5°时间分辨率 1 小时。对比时的常见误区是直接拿单站 30 秒采样 TEC 与 GIM 逐历元比较这没有意义因为 GIM 本身抹平了短时电离层变化。正确做法是先对本地 VTEC 做 10 到 15 分钟的滑动平均再插值到 GIM 对应时刻和穿刺点位置。中纬度白天两者相差 3 到 10 TECU 属于正常范围太阳活动高年会更大。需要明确一点GIM 也是模型值不是真值它更适合用来检查量级和趋势不能作为逐历元真值。5.3 输出带弧段号的 CSV 并检查 DCB 痕迹多颗卫星的 TEC 序列放在一起时不同弧段之间会因为卫星和接收机硬件延迟出现常数阶梯差。这个差值就是 DCB 的表现。单站双频处理无法单独分离接收机与卫星 DCB但可以把它当作常数偏置处理。如果发现某条弧段整体比其他卫星高几十 TECU不一定是计算错误先确认是不是 DCB 导致的。落盘时把弧段号保存下来后续分析会省很多事。result table(t, repmat(first_prn, numel(t), 1), TEC_raw, TEC_smooth, tec_seg, ... VariableNames, {epoch, prn, TEC_raw, TEC_smooth, seg_id}); writetable(result, tec_result.csv);文件里保留seg_id字段后按弧段统计 TEC 均值、筛选连续观测时长、做 DCB 估计都无需再对周跳做二次判断。这也是把处理流程固化成脚本后最低成本保存完整语义的输出方式。本文还有配套的精品资源点击获取
