企业官网建设流程全解析

简介：一套基于MATLAB的GPS双频TEC计算工具包，面向GNSS数据处理、电离层监测等领域的研究人员和学生，用于从RINEX 2.11格式的观测与导航文件计算总电子含量。程序内置完整解算流程，可输出垂直总电子含量（VTEC）、倾斜总电子含量（STEC）、带接收机或卫星差分码偏差（DCB）校正的STEC，以及由码伪距和载波相位分别导出的STEC和ROTI指数，并保留仰角参数，适合分析单日电离层变化。代码依赖Cygwin环境运行，包内附Cygwin安装程序和MEX加速模块，便于有MATLAB与GPS基础的用户直接运行或深度修改。压缩包共84个文件、约28.89MB，含21个MATLAB脚本、46个DCB偏差文件、6个RAR分卷、4个MEX编译模块、3个exe及3个说明文档，并配有1个20n导航电文实测数据。MEX模块用于加速RINEX文件读取，DCB文件用于卫星和接收机偏差校正，脚本覆盖从数据读取、周跳修正到绘图输出的完整流程。目前已有873人学习下载，适合需要快速搭建TEC计算流程、验证电离层监测算法的研究人员参考。

1. 双频 TEC 与 RINEX 2.11 的适配关系

GPS 卫星同时播发 L1 1575.42 MHz 与 L2 1227.60 MHz 两个频率，电离层对这两路信号的群延迟差，恰好正比于信号路径上的总电子含量 TEC。双频接收机把伪距和载波相位按 RINEX 2.11 格式落盘后，计算 TEC 就变成一次文件解析加一次线性组合。RINEX 2.11 的观测文件按 80 字符定长块组织，卫星号、观测值、失锁标识全部挤在同一行里，解析的坑远大于公式本身。这篇文章按「头部类型表 → 逐历元切行 → 几何无关组合 → 相位平滑 → 与 IGS 网格比对」的顺序走完整条链路，适合刚拿到双频接收机原始数据、想在 MATLAB 里快速画出 TEC 时间序列的定位算法工程师和空间物理方向的研究生。

2. RINEX 2.11 观测文件解析：头部观测类型表、历元行与 MATLAB 读取

2.1 # / TYPES OF OBSERV 记录里的双频观测码

RINEX 2.11 的观测文件（O 文件）以RINEX VERSION / TYPE行声明版本，版本号直接写2.11。真正决定怎么切数据的是头部里的# / TYPES OF OBSERV记录，它按顺序列出每个历元后面跟哪些观测量。2.11 版本沿用两位字符的观测码，以下几个是双频 TEC 计算里最常见的组合：

要注意 2.11 里C2的语义在不同厂家的接收机里不完全一致，有些指 L2 C/A，有些指 L2C。如果 O 文件头里同时出现P2和C2，优先用P2参与 TEC 计算，因为 P 码的测距精度和 DCB 标定都更稳定。部分新接收机输出的是 3 位字符的C2S、C2L或C5Q，那已经不是 RINEX 2.11 的约定，直接按 2.12 或 3.x 处理，不要硬塞进这套解析逻辑里。

2.2 按 80 字符块切分卫星号、观测值与 LLI

RINEX 2 系列的数据区每行固定 80 个字符。历元行以>开头，后面依次是年、月、日、时、分、秒，接着是接收机钟差标志和卫星数。观测行前 3 个字符是卫星号，第一位是系统标识（G为 GPS，R为 GLONASS，E为 Galileo，S为 SBAS），后两位是 PRN 号。从第 4 个字符开始，每个观测值占 16 个字符，其中前 14 位是浮点数值，接着 1 位是 LLI 失锁标识，最后 1 位是信噪比等级。

LLI 为 1 或 2 时，代表载波相位发生了失锁或半周模糊度跳变，这一历元的相位观测不能直接进平滑滤波。下面这段 MATLAB 代码先读头部，把观测类型列表存下来，再按类型数量切分观测行：

fid = fopen('obs.rnx', 'r'); % 跳过头部，记录观测类型列表 obsTypes = {}; while true line = fgetl(fid); if contains(line, 'END OF HEADER') break; end if contains(line, '# / TYPES OF OBSERV') parts = regexp(line, '\S+', 'match'); nTypes = str2double(parts{1}); obsTypes = [obsTypes, parts(2:end)]; % 2.11 允许观测类型跨多行，继续读直到凑满 nTypes while numel(obsTypes) < nTypes line = fgetl(fid); parts = regexp(line, '\S+', 'match'); obsTypes = [obsTypes, parts(:)']; end end end

这段代码用regexp按空白切分，不依赖固定列位置，能同时兼容不同厂家对头部空格填充的差异。nTypes决定后续每个历元要读几个观测值，obsTypes的顺序就是观测行里每 16 字符块对应的物理量顺序。读取头部时常见的坑是# / TYPES OF OBSERV一行放不下全部类型，规范允许续行，所以必须用while循环累计，不能只读一行就完事。

2.3 按历元循环、按卫星索引的 MATLAB 数据组织

头部读完以后进入数据区，逐行读取历元行和观测行。观测行里取数值部分用固定列切片，因为str2double对 14 位定宽的纯数字字符串同样适用：

% 历元行解析 epochLine = fgetl(fid); if epochLine(1) ~= '>' error('期望历元行，实际读到其他内容'); end % 2.11 历元行：秒段占 F11.7，之后是钟差标志和卫星数 nSats = str2double(epochLine(31:33)); epochTime = datetime(epochLine(3:28), 'InputFormat', ... 'yy MM dd HH mm ss.SSSSSSS'); % 逐颗卫星读观测行 for s = 1:nSats obsLine = fgetl(fid); prn = obsLine(1:3); obsValues = nan(1, numel(obsTypes)); valid = true; % LLI 标记是否正常 for k = 1:numel(obsTypes) col = 4 + (k - 1) * 16; token = obsLine(col:col + 13); if ~all(isspace(token)) obsValues(k) = str2double(token); end % col+14 是 LLI，非 0 说明相位失锁 if obsLine(col + 14) ~= '0' valid = false; end end % 存入容器：data.(prn).time、data.(prn).C1 等 end

这里的核心是把obsTypes的顺序映射到数值数组的下标。等所有历元读完后，直接用data.G12.C1就能取出某颗卫星的双频伪距序列。LLI的判断放在读行的循环里，一旦发现失锁就把valid置为false，后续相位平滑时这一段的递推状态要重置。

3. 几何无关组合计算 TEC：公式、常量与 MATLAB 实现

3.1 P4 与 L4 的物理意义

电离层对伪距的群延迟与载波相位的相延迟方向相反、数值近似相等：伪距被拉长，相位被提前。取 L1 和 L2 的伪距差：

P4 = P2 - P1 = 40.3 × TEC × (1 / f2² - 1 / f1²)

这里 TEC 的单位是电子数每平方米，频率单位是 Hz，P4 的单位是米。把 L1、L2 的频率代进去，得到工程上最常用的换算关系：

TEC（TECU）≈ 9.52 × (P2 - P1)（米）

1 TECU = 1e16 电子每平方米。这个系数 9.52 是由两个频率的平方倒数差算出来的，换用不同频率对（比如 L1/L5）时系数会变，不能套用。

载波相位的组合 L4 = L2 - L1（均换算成米）方向和 P4 相反，因为相位超前导致观测距离偏短。L4 的噪声只有毫米级，但包含整周模糊度常数项，只能提供 TEC 的相对变化量。把伪距的绝对量级和相位的平滑变化结合起来，就是相位平滑伪距的基本思路。

3.2 伪距 TEC 与相位 TEC 的 MATLAB 代码

在 MATLAB 中按卫星逐历元计算 TEC，先提取已经解析好的 C1 和 P2 序列：

f1 = 1575.42e6; % L1 频率 f2 = 1227.60e6; % L2 频率 K = 1 / (40.3 * (1/f2^2 - 1/f1^2)); % 约 9.52e16，单位 TECU/m % 假设已经按历元对齐了 C1、P2、L1、L2 P4 = P2 - C1; % 群延迟组合，单位米 L1m = L1_l * 0.190293672798365; % L1 相位周数转米 L2m = L2_l * 0.244210213424568; % L2 相位周数转米 L4 = L2m - L1m; % 相位组合，单位米 TEC_pseudo = K * P4; % 伪距 TEC，含 DCB 和噪声 TEC_phase = -K * L4; % 相位 TEC，含模糊度常数

两点说明。第一，K值在代码里保留高精度，不要用 9.52 这种四舍五入后的数去做逐历元计算，否则 100 TECU 的量级上会积累出约 0.5 TECU 的系统差。第二，L1_l是从 RINEX 文件里读出的原始相位周数，必须乘上对应频率的波长才能和伪距做组合。GPS L1 波长为 0.1903 米，L2 波长为 0.2442 米，代码里这两个常量要写全。

3.3 DCB 与多路径：误差预算表

伪距 TEC 里最麻烦的系统误差是差分码偏差 DCB，来自卫星和接收机的硬件群延迟差异。双频伪距组合无法消除 DCB，只能靠外部产品改正。误差预算大致如下：

1 纳秒的群延迟差对应约 2.86 TECU，所以接收机 DCB 是伪距 TEC 最大的污染源。IGS 发布的 DCB 产品按卫星和接收机分别给出 P1-P2 偏差，使用时从对应日期的 DCB 文件里取值，把TEC_pseudo整体平移即可。多路径的影响具有低频特性，在静态观测环境下用相位平滑伪距能把大部分多路径误差压下去，但 DCB 不受平滑影响，必须单独处理。

4. MATLAB 相位平滑伪距与绝对 TEC 弧段拟合

4.1 Hatch 滤波的权重选择

相位平滑伪距的经典实现是 Hatch 递推。设当前历元索引为 k，平滑后的 P4 记为 Ps：

Ps(k) = w × P4(k) + (1 - w) × [Ps(k-1) - (L4(k) - L4(k-1))]

注意 L4 项前面的负号。因为 P4 随 TEC 增大而增大，L4 随 TEC 增大而减小，递推时要把 L4 的变化量取反才能和 P4 对齐。权重 w 有两种取法：一是 w = 1/k，从弧段起点开始渐近收敛，适合长时间连续跟踪；二是固定窗口，w = 1/N，N 取 60 到 120 个历元，对 30 秒采样间隔就是 30 到 60 分钟。静态测站上多路径的周期通常在几分钟到几十分钟，N 太小压制不彻底，N 太大会把电离层真实变化也平滑掉。我一般先看数据弧段的连续长度，短弧段用 w = 1/k，长弧段用固定 N = 100。

4.2 递推平滑代码与周跳重置

下面这段代码按单颗卫星的连续弧段做递推，遇到数据中断或 LLI 异常就重置滤波器：

function TEC_sm = smoothP4(P4, L4, LLI, N) n = numel(P4); TEC_sm = nan(size(P4)); w = 1 / N; % 固定窗口权重 isInit = true; % 当前弧段是否刚刚开始 for k = 1:n if isnan(P4(k)) || isnan(L4(k)) || LLI(k) ~= 0 isInit = true; % 数据缺失或失锁，重置 continue; end if isInit TEC_sm(k) = P4(k); prevL4 = L4(k); prevPs = P4(k); isInit = false; else pred = prevPs - (L4(k) - prevL4); TEC_sm(k) = w * P4(k) + (1 - w) * pred; prevPs = TEC_sm(k); prevL4 = L4(k); end end end

权重w的作用是决定伪距观测对新估计的信任程度。w太大，平滑结果里伪距噪声残留多；w太小，伪距里的多路径和 DCB 偏移会被慢慢当成真实 TEC 吸收，产生低频漂移。LLI(k) ~= 0的判断放在循环开头，保证周跳后的第一历元只用当前伪距初始化，不让跳动后的模糊度污染后续值。

4.3 用夜间参考值消除常数偏置

相位平滑只能消除伪距噪声和大部分多路径，平滑后的结果仍然带着 DCB 和相位模糊度的常数偏移。绝对 TEC 的恢复通常借助夜间电离层极小值：午夜到凌晨 4 点之间，垂直总电子含量降到个位数 TECU，可以把这段时间平滑 TEC 的中位数当作常数偏置，从全天数据里减掉。

% 假设 t 是 UTC 时间数组，TEC_sm 是平滑后的 TEC nightIdx = (hour(t) >= 0) & (hour(t) < 4); bias = median(TEC_sm(nightIdx), 'omitnan'); TEC_abs = TEC_sm - bias;

夜间低电离层假设在中纬度地区基本成立，但在低纬赤道异常区和太阳活动高年，夜间 TEC 仍可能残留 10 TECU 左右的基线，这时用夜间中位数会低估全天 TEC。更稳妥的做法是用当天所有卫星的 TEC 序列做最小二乘：把常数偏置和 DCB 合并成一个待估参数，同时引入卫星 DCB 产品固定卫星端偏差，只估接收机偏差。这套做法在 MATLAB 里用\运算符直接解超定方程即可，比单纯减夜间中位数多一层对 DCB 的显式建模。

5. 与 IGS TEC 网格比对：计算结果的三种验证做法

5.1 广播星历解算卫星位置，转穿刺点坐标

验证前先明确一个边界：上面计算得到的 TEC 是斜路径上的 STEC，而 IGS 发布的全球 TEC 网格是垂直方向 VTEC。要对比，必须先把斜 TEC 投影到垂直方向。投影需要卫星与接收机的几何关系，也就是卫星位置和仰角。常见做法是利用同一时段的广播星历文件，在 MATLAB 里按 ICD-GPS-200 的星历算法解算卫星 ECEF 坐标，然后以测站坐标和卫星坐标计算仰角，再用投影函数：

VTEC = STEC × cos(asin(Re × sin(z) / (Re + h_iono)))

其中 Re 取 6371 公里，电离层薄层高度 h_iono 取 350 公里到 450 公里之间。IGS 网格产品采用的等效薄层高度是 450 公里，做严格比对应保持同一取值。投影函数的误差在中低仰角下显著，建议比对时只取仰角大于 30 度的卫星。

5.2 IONEX 网格插值比对

IGS 的 IONEX 文件按经纬度 5 度 × 2.5 度给出 VTEC 网格，时间分辨率常见为 2 小时或 1 小时，MATLAB 读取后可以用interpn插值到接收机穿刺点位置：

% latVec, lonVec, tVec 为 IONEX 网格坐标 % vtecGrid 为 3 维数组，顺序对应 lat/lon/time vtecIgs = interpn(latVec, lonVec, tVec, vtecGrid, ... latIPP, lonIPP, tIPP, 'linear');

插值完成后按卫星、按历元画散点图，横轴是本文计算的垂直 TEC，纵轴是 IGS 网格值，理想情况下点云围绕斜率 1 的对角线分布。误差来源要分别看：IGS 网格本身在中纬度的精度约 1 到 3 TECU，低纬和高纬地区会更大；穿刺点投影函数的薄层高度假设会带来系统偏差；接收机 DCB 改正不彻底则表现为整体平移。用表格记录不同卫星的统计结果，比只画一个总图更容易定位问题。

5.3 逐卫星残差诊断的脚本化

收尾的实用技巧是把验证做成可复现脚本，输出三列关键指标：每颗卫星的 VTEC 均值、与 IGS 网格的差值中位数、差值标准差。差值中位数非零说明 DCB 或投影函数有系统偏差，差值标准差偏大则指向周跳处理不干净或多路径严重。把这套指标按天滚动输出，能直接看出接收机 DCB 的日稳定性，也可以反过来验证平滑窗口 N 选得是否合理。

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