企业官网建设流程全解析

简介：针对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 流程中的角色。

不值得在解析阶段花太多精力处理 S1/S2，TEC 精度主要由伪距和相位决定。信噪比可以在后续质量控制里做阈值筛选，但大多数单站 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 MHz，L2 频率为 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)×λ2，w 是平滑权重，取值在 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 end

ps_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 km，H 取 350 km 或 450 km，VTEC = 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，对高仰角卫星几乎没有差别，主要影响低仰角弧段的投影幅度，建议固定一个值并写进处理日志。

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 估计都无需再对周跳做二次判断。这也是把处理流程固化成脚本后，最低成本保存完整语义的输出方式。

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