简介面向GPS电离层研究与导航定位应用人员这份Matlab程序集利用GPS观测数据与精密星历基于多项式或球谐函数模型实现接收机与卫星DCB的同时估计并进一步计算电离层TEC为导航单位提供电离层延迟改正。资源共20个文件、约363KB其中14个.m脚本覆盖主流程与关键算法包括RINEX/SP3读取、穿刺点IPP解算、MDCB/SDCB估计、坐标转换等模块3个.sp3精密星历文件构成完整算例另有.mat参数文件与1个.05i电离层文件便于对照验证与二次开发。已有495人学习适合具有GPS数据处理基础的测绘、导航或空间天气方向学生与工程师。借助这套代码框架用户可以快速开展区域电离层建模实验修改多项式阶数或改用球谐函数模型时也只需调整对应脚本能有效节省从零搭建的时间。1. 用MATLAB把GPS观测、精密星历和DCB估计串成一条TEC产线第一道坎在哪手里有双频GPS观测文件和IGS精密星历想自己算TEC第一道坎就把人卡在原地P2和P1的伪距差不等于电离层延迟里面叠着接收机DCB和卫星DCB两种硬件延迟偏差。不把它们同时估计出来算出的TEC会稳定偏出几个TECU拿去做单频电离层延迟改正结果越改越差。闭环做法是精密星历给位置多项式或球谐函数给VTEC模型接收机DCB、卫星DCB和模型系数放进同一个最小二乘里一起解最后输出TEC并换算成L1/L2延迟改正量。这套MATLAB流程适合做GNSS电离层监测、单频定位改正或DCB标定的研究生和工程人员。2. 从双频伪距到STEC为什么DCB必须和TEC在同一个平差里解出来2.1 双频观测怎么换算成STECDCB在公式里站在哪GPS双频观测提供L1/L2载波相位和P1/P2码伪距。电离层对微波信号是色散介质相延迟超前、群延迟滞后数值上都是40.3×TEC/f^2单位米。把P2减P1几何距离、钟差、对流层这些频率无关项被消掉剩下的是频率相关的电离层项和硬件延迟P2 − P1 40.3·STEC·(1/f2^2 − 1/f1^2) DCB_r DCB^s ε其中f1 1575.42 MHzf2 1227.60 MHz。因为f1 f2括号里是正数所以P2 − P1越大斜向TEC越大。DCB_r是接收机两个频率通道的差分码偏差DCB^s是同一颗卫星的差分码偏差ε则是码噪声和多路径。DCB的单位习惯用纳秒1 ns对应约0.3 m换算到TEC约2~3 TECU。低纬地区白天峰值可能到100 TECU这个偏差占比不大但做全天TEC图或单频改正时它是系统性的不会被平均掉。这就是为什么不能先算TEC再回头修改DCB。STEC里带着一个常数偏置而VTEC模型的第一项也是常数项两者在最小二乘里高度相关。分两步做的结果往往是模型系数把DCB吸收得干干净净TEC图看着平滑DCB却完全错位。所以标题里“同时估计”不是锦上添花是这类处理唯一稳的路径。观测噪声方面码伪距的噪声在分米到米级直接换算STEC会有1~3 TECU的抖动。工程上一般先对P1/P2做载波相位平滑伪距把噪声压到0.2 TECU以内再做差分。Hatch滤波的递推实现很简单单频平滑就可以function P_sm hatchSmooth(P, L, alpha) % P: 原始码伪距序列单位米Nx1 % L: 载波相位观测已乘波长转为米Nx1 % alpha: 平滑系数0.01~0.05越小平滑越强对周跳越敏感 P_sm NaN(size(P)); P_sm(1) P(1); for k 2 : numel(P) P_sm(k) alpha * P(k) (1-alpha) * (P_sm(k-1) (L(k) - L(k-1))); end end这段递推的物理含义是相位历元间差分很干净把它累加在上一历元平滑值上作为当前历元的预测再用当前码观测修正。alpha取0.02时码噪声被压到约1/50但要注意周跳。如果L(k)−L(k−1)跳了整周预测值会突然偏差平滑输出跟着翻车。所以平滑前必须先做周跳探测第5章会专门讲这个坑。2.2 多项式与球谐函数模型参数数量、适用范围和选择逻辑双频观测直接得到的是STEC是倾斜方向电子总含量。要建模、做网格或改正单频用户需要先把STEC投影到垂直方向VTECSTEC VTEC·MF(χ′)χ′是信号在电离层穿刺点处的天顶角。穿刺点地理位置每颗卫星每个历元都在变VTEC分布需要一个函数来表达。“多项式模型”和“球谐函数模型”是两种最常见的表达方式。多项式模型把VTEC写成穿刺点纬度和经度差的高次多项式。以太阳-地磁坐标系里的太阳时角差和地磁纬度差为自变量2阶展开就是VTEC a00 a10·(φ−φ0) a01·(λ−λ0) a20·(φ−φ0)^2 a11·(φ−φ0)·(λ−λ0) a02·(λ−λ0)^26个参数就能描述单站或小区域的TEC空间变化。P阶多项式的参数个数是(P1)(P2)/2P2是6个P3是10个。好处是数值稳定、对观测分布不敏感坏处是大范围区域里真实TEC的空间梯度变化太快多项式会穿出负值或出现边缘翘曲。球谐函数模型的表达是VTEC(φ′, λ′) Σ Σ P̄_nm(sin φ′)·(a_nm·cos(mλ′) b_nm·sin(mλ′))其中P̄_nm是归一化连带勒让德函数φ′和λ′取太阳-地磁坐标系最高阶取N时参数个数是(N1)^2。N4时25个参数N8时81个参数。球谐函数是球面上的正交基函数比多项式更适合做大范围的全局TEC建模代价是和测站分布的均匀性强绑定在高阶项上容易吸收噪声。我一般这样选单站或几百公里范围内的区域网用2~3阶多项式参数少抗差性好跨几个时区或者做全球拼接用4~8阶球谐测站分布不均匀时宁可用4阶球谐加正则化也不要一口气上8阶。模型选择直接影响第4章设计矩阵的尺度和病态程度选型优先级不低。3. 用精密星历做位置和穿刺点准备SP3解析、内插与IPP计算的MATLAB实现3.1 SP3文本解析与拉格朗日内插先把卫星坐标变成逐历元序列算高度角、方位角和穿刺点位置前提是拿到毫米级到厘米级的卫星坐标。广播星历的轨道误差在米级对高度角的影响可以忽略但由高度角算穿透点坐标时存在剖面放大效应所以这里用IGS精密星历。SP3格式是文本文件每个历元一行带*接下来是各颗卫星的PG行function [ephTime, satPos] readSP3(sp3file, prnList) % 读取SP3文件prnList例如[1 3 5 7 9 12 15 18 21 24 27 30] % 返回历元时间序列和对应PRN的地固系坐标(Nx3) fid fopen(sp3file, rt); ephTime []; satPos NaN(5000, 32, 3); % 300秒采样一天约288个历元取5000足够 tIdx 0; while ~feof(fid) tline fgetl(fid); if isempty(tline) continue; end if tline(1) * tIdx tIdx 1; ymdhms sscanf(tline(2:end), %d); ephTime(tIdx) datenum(ymdhms(1), ymdhms(2), ymdhms(3), ... ymdhms(4), ymdhms(5), ymdhms(6)); elseif strncmp(tline, PG, 2) prn str2double(tline(3:4)); if ismember(prn, prnList) x str2double(tline(7:18)); y str2double(tline(20:31)); z str2double(tline(33:44)); satPos(tIdx, prn, :) [x y z]; end end end fclose(fid); endSP3的位置记录是固定列宽第3~4列是PRN号之后每隔14列一个坐标分量。用sscanf按固定列宽提取也可以但str2double配列下标更直观。读取后要把无效历元滤掉SP3文件里可能没有某些PRN的位置记录。得到480秒或300秒采样点的卫星坐标后再根据观测历元内插。线性插值肯定不行卫星在轨道弧段上的运动是光滑曲线线性插值会引入厘米到分米级的误差。常见做法是取插值历元前后各5个历元共11个点做10阶拉格朗日内插。MATLAB里没有内置拉格朗日函数可以这样写function yq lagrange11(t, y, tq) % t: 星历时间序列y: 对应坐标序列tq: 插值时间 % 取最靠近tq的11个点做拉格朗日内插 n numel(t); [~, idx] min(abs(t - tq)); k max(1, min(idx - 5, n - 11)); % 保证窗口内有11个点 tseg t(k : k10); yseg y(k : k10); yq 0; for j 1 : 11 lj 1; for i 1 : 11 if i ~ j lj lj * (tq - tseg(i)) / (tseg(j) - tseg(i)); end end yq yq lj * yseg(j); end end窗口位置由minmax钳制到索引区间[1, n-11]内避免开头结尾取不满11点时基函数发散。插值前要把时间统一成数值型datenum的浮点在相减时精度足够不用担心毫秒级截断。边缘历元的插值振荡仍然存在这点在第5章里要专门处理。3.2 站心系下的高度角与穿刺点IPP纬度经度和映射函数的代码有了卫星位置和接收机坐标先把ECEF差分转成站心ENU。接收机位置从RINEX观测文件头部拿或者用已知测站坐标。ECEF到ENU的旋转矩阵和高度角、方位角计算放在同一个函数里function [el, az] ecef2enuElAz(lat, lon, dr) % dr: 卫星ECEF坐标减接收机ECEF坐标单位米 % lat/lon: 接收机地理纬度/经度单位弧度 slat sin(lat); clat cos(lat); slon sin(lon); clon cos(lon); e -slon .* dr(1) clon .* dr(2); n -slat .* clon .* dr(1) - slat .* slon .* dr(2) clat .* dr(3); u clat .* clon .* dr(1) clat .* slon .* dr(2) slat .* dr(3); el atan2(u, sqrt(e.^2 n.^2)); az atan2(e, n); end高度角低于10°~15°的观测直接扔掉。低高度角信号穿过电离层路径很长映射函数误差和码多路径都放大在DCB估计里很难用。这一步滤完剩下的观测基本上都是干净几何。穿刺点位置要按单层电离层模型来算。设电离层等效高度hIon 350 km地球平均半径Re 6371 km先由高度角el求穿刺点处天顶角χ′χ′ asin( Re / (Re hIon) × cos(el) )然后按球面三角从接收机位置往方位角方向推穿刺点的地心角ψ π/2 − el − χ′穿刺点纬度φ_pp asin(sinφ·cosψ cosφ·sinψ·cosA)经度λ_pp由东向方位角分量决定。代码function [latPp, lonPp, chi] ionoPiercePoint(lat, lon, el, az, hIon) % 单层电离层模型输入输出均为弧度hIon单位米 Re 6371000; psi acos(Re / (Re hIon) * cos(el)); % 地心角 chi pi/2 - el - psi; % 穿刺点天顶角 latPp asin(sin(lat)*cos(psi) cos(lat)*sin(psi)*cos(az)); lonPp lon asin(sin(az)*sin(psi) / cos(latPp)); end映射函数MF 1/cos(χ′)STEC VTEC × MF。这里cos(χ′)在χ′很大时会很小把STEC放大所以低高度角观测在模型估计里权重必须压低。IPP的经度要先把太阳-地磁坐标的转换做好再做球谐展开我这套流程里对多项式模型用地理坐标也能出结果但换成球谐就必须先转坐标系否则高阶项会吸收掉一个系统性的日变化收敛值解释起来很别扭。4. 同一套解算同时给卫星DCB、接收机DCB和TEC模型系数法方程与约束4.1 把观测方程写成A矩阵的一行STEC的映射函数位置别放错现在每一条观测都齐了STEC观测值、穿刺点坐标、高度角、卫星编号、接收机编号。把观测方程写成STEC_i MF(χ′_i) · Σ (模型基函数在IPP处的值 · 对应系数) DCB_r DCB^s ε_i未知数分三块模型系数多项式6~10个或球谐25~81个、接收机DCB每个测站1个、卫星DCB每颗卫星1个。所有未知数放进一个向量x则每个观测对应设计矩阵A的一行。这里有一个很容易放反的细节。STEC是斜向的模型是VTEC因此模型列方向是MF乘以模型基函数不是除以。很多实现把MF放反结果解出来的DCB和TEC整体变形单看残差又看不出来。A矩阵一行的代码示意如下function aRow buildModelRow(MF, base, nModel, iRx, iSv, nRx) % MF: 映射函数值base: 模型基函数在IPP处的取值行向量 % nModel: 模型系数个数iRx/iSv: 接收机/卫星编号nRx: 接收机总数 aRow zeros(1, nModel nRx nSv); aRow(1:nModel) MF .* base; % 注意乘MF不是除 aRow(nModel iRx) 1; % 接收机DCB列 aRow(nModel nRx iSv) 1; % 卫星DCB列 end球谐基函数生成时用MATLAB自带的legendre函数加norm归一化按(n,m)顺序铺开即可function base shBasis(lat, lon, nmax) % 基函数按(0,0),(1,0),(1,1c),(1,1s),(2,0)...顺序排列 % 输出尺寸为(nmax1)^2 × numel(lat) num numel(lat); base zeros((nmax1)^2, num); k 0; for n 0 : nmax pnm legendre(n, sin(lat), norm); for m 0 : n if m 0 k k 1; base(k, :) pnm(m1, :) .* cos(m * lon(:)); else k k 1; base(k, :) pnm(m1, :) .* cos(m * lon(:)); k k 1; base(k, :) pnm(m1, :) .* sin(m * lon(:)); end end end endm0时只有余弦项m0时同一阶次拆成余弦和正弦两项这样球谐基函数总数为(nmax1)^2。所有观测逐行压入A再把对应的STEC观测值按同样的顺序放入b最后调用xHat (A * A) \ (A * b);法方程规模不大模型81个参数、30颗卫星、几十个测站也就是一百多阶直接求逆没问题。这一行解出的是模型系数、卫星DCB和接收机DCB的联合后验估计。4.2 秩亏与基准约束为什么卫星DCB之和为零是常用假设第一次按上面的代码跑xHat几乎必报奇异或给出巨大的DCB值这是秩亏。原因很直观卫星DCB加接收机DCB任意把其中一个整体平移一个常数加进另一个残差完全不变。也就是说“所有卫星DCB整体1所有接收机DCB整体−1”这条方向在法矩阵里没有信息正规矩阵奇异。破解办法是给解空间补一个基准。最常用的是假设一个解算时段内所有参与观测的卫星DCB之和为零Σ DCB^s 0在法方程里增加一行或直接对A做行扩展。扩展方式如下% 在A和b底部追加卫星DCB零和约束 A(end1, nModelnRx1 : end) 1; b(end1, 1) 0; % 给约束行一个比观测小很多的标准差例如1e-6 sigma 1e-6; A(end, :) A(end, :) / sigma; b(end) b(end) / sigma;另一种做法是固定一个已知卫星的DCB例如选一颗健康度好的卫星用IGS发布的产品直接固定这样接收机DCB全部相对该卫星确定。但固定单星会把那颗卫星的DCB产品误差直接引入所有接收机稳定性不如约束所有卫星之和。工程上我一般首选零和约束只有处理单站数据或者卫星数据少时才改用固定单星。参数设置里还有几个和法方程尺度相关的项。观测权重按高度角给常见取sin^2(el)或1/sin²(el)把低高度角观测的误差放大倍数抵消掉。球谐阶数固定时不要一开始就上8阶先4阶跑通再逐级看残差是否有明显下降。正则化方面高阶系数方差不稳时就直接把法矩阵对角加一个1e-6~1e-8的小量虽然破坏严格无偏但能压住震荡。下面是一个我常用的参数表直接复现时可以先照这个设参数取值说明高度角截止10°~15°低高度角映射函数误差大多项式阶数2~3单站或小区域用球谐最高阶4~8测站密时取高阶穿刺点高度350 km也可按450 km平滑系数0.02载波相位平滑伪距卫星DCB约束零和秩亏基准正则化1e-7法矩阵对角加小量这组参数在小范围单站网和大范围IGS网里都跑得通。区别在观测量的数量单站解算时一颗卫星的弧段很短DCB和模型系数的相关性会更高这时约束权重要适当加大。5. 高频避坑指南DCB解算跑飞、TEC跳变的五个真实原因与处理5.1 夜间与低纬度残差不收敛DCB解出来天天飘现象同一台接收机连续两天的DCB估计值相差几个纳秒白天结果看着还行夜间历元残差成片偏大。 原因夜间电离层TEC很低STEC观测值小DCB的相对占比变大加上夜间观测几何变化缓慢法方程对DCB的约束变弱。低纬度地区还有等离子体层和闪烁影响TEC本身不够平滑。 解决夜间弧段单独降低权重或把截止高度角从10°提到15°如果只是要DCB尽量用白天6点到22点的高高度角数据反演夜间数据留给验证用。残差若在固定某颗卫星上成片异常优先怀疑该卫星弧段的相位平滑问题而不是DCB变化。5.2 卫星DCB之和为零的基准在变化连续解出现系统跳变现象连续三天的卫星DCB解出来每颗卫星的DCB当天内部平滑但相邻天之间有整体平移且各颗星平移量相同。 原因正是零和约束自身造成的。零和约束没有把“绝对基准”钉死只是把一天内卫星DCB均值清零。如果某天参与解算的卫星集合不同或某颗卫星数据质量差整体均值会变化所有卫星DCB一起平移。 解决把连续解之间的卫星均值差作为整体偏差先求当天卫星DCB均值再对均值做时间平滑或直接对齐到IGS产品也可以把已知卫星的DCB做弱约束给IGS产品一个0.5 ns的容许量这样基准被软绑定又不会完全继承产品误差。5.3 高阶球谐出现东西向条带或负值现象球谐阶数从4提到8后TEC图出现规律东西向条带低纬区域出现局部负TEC物理上说不通。 原因高阶球谐的自由度增加把测站分布不均匀和观测噪声直接拟合进了模型。低纬地区穿刺点覆盖密集高阶项被噪声激发出现振荡。 解决第一步先检查穿刺点覆盖网格内观测少于一定数量就降阶第二步在法方程对角线加正则化重点对最高阶两阶的参数加大惩罚第三步用多项式模型对比一下如果多项式结果反而更平滑说明观测密度撑不起高阶球谐。5.4 载波相位平滑被周跳污染现象某颗卫星的STEC时间序列出现一个台阶前后两段波形都正常中间跳变残差对应历元出现团状异常。 原因Hatch滤波假设相位历元间差分无整周跳变。周跳发生时L(k)−L(k−1)跳了整周平滑器把伪距改正带偏输出在数分钟内被污染。没有做周跳探测就平滑是所有DCB估计翻车里最隐蔽的一个。 解决平滑前用GF组合或MW组合探测周跳发现跳变就把平滑窗口重置为当前历元。我通常的做法是用L1−L2的多路径组合做粗检测再用电离层残差变化率做精检测周跳处直接重新初始化平滑器。5.5 SP3插值在边缘历元产生假坐标现象测站全天观测数据里每天开头半小时和结尾半小时的DCB解算残差系统性偏大中间时段完全正常。 原因第3章里写的拉格朗日11点插值在数据两端只能取到一侧的点窗口钳位后基函数振荡卫星坐标偏差达到分米量级直接污染高度角和穿刺点。 解决解算时丢弃一天中卫星星历覆盖边缘的时段更稳妥的是用前后两天的SP3文件拼成超长序列再做插值让目标历元始终落在窗口中间。插值误差不直接检查坐标可以在解算后用无电离层组合相位残差做间接验证。6. 残差驱动迭代与验后验证把单天解练到和IGS产品对得上的技巧前面五章跑完能拿到一组可看的DCB和TEC但离稳定可发布还差两步残差清洗和验后质量检查。我现在的习惯做法是两轮加权最小二乘第一步用原始观测跑一次平差保存每个观测的残差第二步按残差绝对值排序剔除超过3倍RMS离差的观测并给剩余观测重新按高度角加权第三步用清洗后的观测重复解算再检查卫星DCB与IGS产品差值的RMS若小于0.5 ns就算稳定第四步用解出的参数回代VTEC按1°×1°网格化输出与同一时段的CODE GIM网格做差RMS在2~4 TECU说明链路基本正常。验证时不只盯DCB还要看TEC图有没有物理上可疑的负值或条带。如果TEC与GIM相差大但DCB差值很小问题多半出在球谐阶数和穿刺点映射上反过来DCB差值大优先怀疑某个接收机的平滑伪距序列还有周跳残留。一个实用的收敛技巧是第一轮用sin²(el)定权第二轮改用残差方差定权测站级重新缩放。伪距噪声水平在测站间差异很大单靠高度角定权会把噪声大的站顶得太高。残差定权后低质量测站的影响被自动削弱DCB的稳定性明显提升。这套流程里最容易忽略的是数据质量标签。RINEX文件里的天线类型、接收机类型和测站概略坐标会直接影响接收机坐标计算坐标差几米穿刺点位置就差几千米低高度角观测全被污染。我早年在这上面吃过亏现在所有测站坐标都会和IGS SINEX比对一遍再开工。希望帮到你。本文还有配套的精品资源点击获取
简介面向GPS电离层研究与导航定位应用人员这份Matlab程序集利用GPS观测数据与精密星历基于多项式或球谐函数模型实现接收机与卫星DCB的同时估计并进一步计算电离层TEC为导航单位提供电离层延迟改正。资源共20个文件、约363KB其中14个.m脚本覆盖主流程与关键算法包括RINEX/SP3读取、穿刺点IPP解算、MDCB/SDCB估计、坐标转换等模块3个.sp3精密星历文件构成完整算例另有.mat参数文件与1个.05i电离层文件便于对照验证与二次开发。已有495人学习适合具有GPS数据处理基础的测绘、导航或空间天气方向学生与工程师。借助这套代码框架用户可以快速开展区域电离层建模实验修改多项式阶数或改用球谐函数模型时也只需调整对应脚本能有效节省从零搭建的时间。1. 用MATLAB把GPS观测、精密星历和DCB估计串成一条TEC产线第一道坎在哪手里有双频GPS观测文件和IGS精密星历想自己算TEC第一道坎就把人卡在原地P2和P1的伪距差不等于电离层延迟里面叠着接收机DCB和卫星DCB两种硬件延迟偏差。不把它们同时估计出来算出的TEC会稳定偏出几个TECU拿去做单频电离层延迟改正结果越改越差。闭环做法是精密星历给位置多项式或球谐函数给VTEC模型接收机DCB、卫星DCB和模型系数放进同一个最小二乘里一起解最后输出TEC并换算成L1/L2延迟改正量。这套MATLAB流程适合做GNSS电离层监测、单频定位改正或DCB标定的研究生和工程人员。2. 从双频伪距到STEC为什么DCB必须和TEC在同一个平差里解出来2.1 双频观测怎么换算成STECDCB在公式里站在哪GPS双频观测提供L1/L2载波相位和P1/P2码伪距。电离层对微波信号是色散介质相延迟超前、群延迟滞后数值上都是40.3×TEC/f^2单位米。把P2减P1几何距离、钟差、对流层这些频率无关项被消掉剩下的是频率相关的电离层项和硬件延迟P2 − P1 40.3·STEC·(1/f2^2 − 1/f1^2) DCB_r DCB^s ε其中f1 1575.42 MHzf2 1227.60 MHz。因为f1 f2括号里是正数所以P2 − P1越大斜向TEC越大。DCB_r是接收机两个频率通道的差分码偏差DCB^s是同一颗卫星的差分码偏差ε则是码噪声和多路径。DCB的单位习惯用纳秒1 ns对应约0.3 m换算到TEC约2~3 TECU。低纬地区白天峰值可能到100 TECU这个偏差占比不大但做全天TEC图或单频改正时它是系统性的不会被平均掉。这就是为什么不能先算TEC再回头修改DCB。STEC里带着一个常数偏置而VTEC模型的第一项也是常数项两者在最小二乘里高度相关。分两步做的结果往往是模型系数把DCB吸收得干干净净TEC图看着平滑DCB却完全错位。所以标题里“同时估计”不是锦上添花是这类处理唯一稳的路径。观测噪声方面码伪距的噪声在分米到米级直接换算STEC会有1~3 TECU的抖动。工程上一般先对P1/P2做载波相位平滑伪距把噪声压到0.2 TECU以内再做差分。Hatch滤波的递推实现很简单单频平滑就可以function P_sm hatchSmooth(P, L, alpha) % P: 原始码伪距序列单位米Nx1 % L: 载波相位观测已乘波长转为米Nx1 % alpha: 平滑系数0.01~0.05越小平滑越强对周跳越敏感 P_sm NaN(size(P)); P_sm(1) P(1); for k 2 : numel(P) P_sm(k) alpha * P(k) (1-alpha) * (P_sm(k-1) (L(k) - L(k-1))); end end这段递推的物理含义是相位历元间差分很干净把它累加在上一历元平滑值上作为当前历元的预测再用当前码观测修正。alpha取0.02时码噪声被压到约1/50但要注意周跳。如果L(k)−L(k−1)跳了整周预测值会突然偏差平滑输出跟着翻车。所以平滑前必须先做周跳探测第5章会专门讲这个坑。2.2 多项式与球谐函数模型参数数量、适用范围和选择逻辑双频观测直接得到的是STEC是倾斜方向电子总含量。要建模、做网格或改正单频用户需要先把STEC投影到垂直方向VTECSTEC VTEC·MF(χ′)χ′是信号在电离层穿刺点处的天顶角。穿刺点地理位置每颗卫星每个历元都在变VTEC分布需要一个函数来表达。“多项式模型”和“球谐函数模型”是两种最常见的表达方式。多项式模型把VTEC写成穿刺点纬度和经度差的高次多项式。以太阳-地磁坐标系里的太阳时角差和地磁纬度差为自变量2阶展开就是VTEC a00 a10·(φ−φ0) a01·(λ−λ0) a20·(φ−φ0)^2 a11·(φ−φ0)·(λ−λ0) a02·(λ−λ0)^26个参数就能描述单站或小区域的TEC空间变化。P阶多项式的参数个数是(P1)(P2)/2P2是6个P3是10个。好处是数值稳定、对观测分布不敏感坏处是大范围区域里真实TEC的空间梯度变化太快多项式会穿出负值或出现边缘翘曲。球谐函数模型的表达是VTEC(φ′, λ′) Σ Σ P̄_nm(sin φ′)·(a_nm·cos(mλ′) b_nm·sin(mλ′))其中P̄_nm是归一化连带勒让德函数φ′和λ′取太阳-地磁坐标系最高阶取N时参数个数是(N1)^2。N4时25个参数N8时81个参数。球谐函数是球面上的正交基函数比多项式更适合做大范围的全局TEC建模代价是和测站分布的均匀性强绑定在高阶项上容易吸收噪声。我一般这样选单站或几百公里范围内的区域网用2~3阶多项式参数少抗差性好跨几个时区或者做全球拼接用4~8阶球谐测站分布不均匀时宁可用4阶球谐加正则化也不要一口气上8阶。模型选择直接影响第4章设计矩阵的尺度和病态程度选型优先级不低。3. 用精密星历做位置和穿刺点准备SP3解析、内插与IPP计算的MATLAB实现3.1 SP3文本解析与拉格朗日内插先把卫星坐标变成逐历元序列算高度角、方位角和穿刺点位置前提是拿到毫米级到厘米级的卫星坐标。广播星历的轨道误差在米级对高度角的影响可以忽略但由高度角算穿透点坐标时存在剖面放大效应所以这里用IGS精密星历。SP3格式是文本文件每个历元一行带*接下来是各颗卫星的PG行function [ephTime, satPos] readSP3(sp3file, prnList) % 读取SP3文件prnList例如[1 3 5 7 9 12 15 18 21 24 27 30] % 返回历元时间序列和对应PRN的地固系坐标(Nx3) fid fopen(sp3file, rt); ephTime []; satPos NaN(5000, 32, 3); % 300秒采样一天约288个历元取5000足够 tIdx 0; while ~feof(fid) tline fgetl(fid); if isempty(tline) continue; end if tline(1) * tIdx tIdx 1; ymdhms sscanf(tline(2:end), %d); ephTime(tIdx) datenum(ymdhms(1), ymdhms(2), ymdhms(3), ... ymdhms(4), ymdhms(5), ymdhms(6)); elseif strncmp(tline, PG, 2) prn str2double(tline(3:4)); if ismember(prn, prnList) x str2double(tline(7:18)); y str2double(tline(20:31)); z str2double(tline(33:44)); satPos(tIdx, prn, :) [x y z]; end end end fclose(fid); endSP3的位置记录是固定列宽第3~4列是PRN号之后每隔14列一个坐标分量。用sscanf按固定列宽提取也可以但str2double配列下标更直观。读取后要把无效历元滤掉SP3文件里可能没有某些PRN的位置记录。得到480秒或300秒采样点的卫星坐标后再根据观测历元内插。线性插值肯定不行卫星在轨道弧段上的运动是光滑曲线线性插值会引入厘米到分米级的误差。常见做法是取插值历元前后各5个历元共11个点做10阶拉格朗日内插。MATLAB里没有内置拉格朗日函数可以这样写function yq lagrange11(t, y, tq) % t: 星历时间序列y: 对应坐标序列tq: 插值时间 % 取最靠近tq的11个点做拉格朗日内插 n numel(t); [~, idx] min(abs(t - tq)); k max(1, min(idx - 5, n - 11)); % 保证窗口内有11个点 tseg t(k : k10); yseg y(k : k10); yq 0; for j 1 : 11 lj 1; for i 1 : 11 if i ~ j lj lj * (tq - tseg(i)) / (tseg(j) - tseg(i)); end end yq yq lj * yseg(j); end end窗口位置由minmax钳制到索引区间[1, n-11]内避免开头结尾取不满11点时基函数发散。插值前要把时间统一成数值型datenum的浮点在相减时精度足够不用担心毫秒级截断。边缘历元的插值振荡仍然存在这点在第5章里要专门处理。3.2 站心系下的高度角与穿刺点IPP纬度经度和映射函数的代码有了卫星位置和接收机坐标先把ECEF差分转成站心ENU。接收机位置从RINEX观测文件头部拿或者用已知测站坐标。ECEF到ENU的旋转矩阵和高度角、方位角计算放在同一个函数里function [el, az] ecef2enuElAz(lat, lon, dr) % dr: 卫星ECEF坐标减接收机ECEF坐标单位米 % lat/lon: 接收机地理纬度/经度单位弧度 slat sin(lat); clat cos(lat); slon sin(lon); clon cos(lon); e -slon .* dr(1) clon .* dr(2); n -slat .* clon .* dr(1) - slat .* slon .* dr(2) clat .* dr(3); u clat .* clon .* dr(1) clat .* slon .* dr(2) slat .* dr(3); el atan2(u, sqrt(e.^2 n.^2)); az atan2(e, n); end高度角低于10°~15°的观测直接扔掉。低高度角信号穿过电离层路径很长映射函数误差和码多路径都放大在DCB估计里很难用。这一步滤完剩下的观测基本上都是干净几何。穿刺点位置要按单层电离层模型来算。设电离层等效高度hIon 350 km地球平均半径Re 6371 km先由高度角el求穿刺点处天顶角χ′χ′ asin( Re / (Re hIon) × cos(el) )然后按球面三角从接收机位置往方位角方向推穿刺点的地心角ψ π/2 − el − χ′穿刺点纬度φ_pp asin(sinφ·cosψ cosφ·sinψ·cosA)经度λ_pp由东向方位角分量决定。代码function [latPp, lonPp, chi] ionoPiercePoint(lat, lon, el, az, hIon) % 单层电离层模型输入输出均为弧度hIon单位米 Re 6371000; psi acos(Re / (Re hIon) * cos(el)); % 地心角 chi pi/2 - el - psi; % 穿刺点天顶角 latPp asin(sin(lat)*cos(psi) cos(lat)*sin(psi)*cos(az)); lonPp lon asin(sin(az)*sin(psi) / cos(latPp)); end映射函数MF 1/cos(χ′)STEC VTEC × MF。这里cos(χ′)在χ′很大时会很小把STEC放大所以低高度角观测在模型估计里权重必须压低。IPP的经度要先把太阳-地磁坐标的转换做好再做球谐展开我这套流程里对多项式模型用地理坐标也能出结果但换成球谐就必须先转坐标系否则高阶项会吸收掉一个系统性的日变化收敛值解释起来很别扭。4. 同一套解算同时给卫星DCB、接收机DCB和TEC模型系数法方程与约束4.1 把观测方程写成A矩阵的一行STEC的映射函数位置别放错现在每一条观测都齐了STEC观测值、穿刺点坐标、高度角、卫星编号、接收机编号。把观测方程写成STEC_i MF(χ′_i) · Σ (模型基函数在IPP处的值 · 对应系数) DCB_r DCB^s ε_i未知数分三块模型系数多项式6~10个或球谐25~81个、接收机DCB每个测站1个、卫星DCB每颗卫星1个。所有未知数放进一个向量x则每个观测对应设计矩阵A的一行。这里有一个很容易放反的细节。STEC是斜向的模型是VTEC因此模型列方向是MF乘以模型基函数不是除以。很多实现把MF放反结果解出来的DCB和TEC整体变形单看残差又看不出来。A矩阵一行的代码示意如下function aRow buildModelRow(MF, base, nModel, iRx, iSv, nRx) % MF: 映射函数值base: 模型基函数在IPP处的取值行向量 % nModel: 模型系数个数iRx/iSv: 接收机/卫星编号nRx: 接收机总数 aRow zeros(1, nModel nRx nSv); aRow(1:nModel) MF .* base; % 注意乘MF不是除 aRow(nModel iRx) 1; % 接收机DCB列 aRow(nModel nRx iSv) 1; % 卫星DCB列 end球谐基函数生成时用MATLAB自带的legendre函数加norm归一化按(n,m)顺序铺开即可function base shBasis(lat, lon, nmax) % 基函数按(0,0),(1,0),(1,1c),(1,1s),(2,0)...顺序排列 % 输出尺寸为(nmax1)^2 × numel(lat) num numel(lat); base zeros((nmax1)^2, num); k 0; for n 0 : nmax pnm legendre(n, sin(lat), norm); for m 0 : n if m 0 k k 1; base(k, :) pnm(m1, :) .* cos(m * lon(:)); else k k 1; base(k, :) pnm(m1, :) .* cos(m * lon(:)); k k 1; base(k, :) pnm(m1, :) .* sin(m * lon(:)); end end end endm0时只有余弦项m0时同一阶次拆成余弦和正弦两项这样球谐基函数总数为(nmax1)^2。所有观测逐行压入A再把对应的STEC观测值按同样的顺序放入b最后调用xHat (A * A) \ (A * b);法方程规模不大模型81个参数、30颗卫星、几十个测站也就是一百多阶直接求逆没问题。这一行解出的是模型系数、卫星DCB和接收机DCB的联合后验估计。4.2 秩亏与基准约束为什么卫星DCB之和为零是常用假设第一次按上面的代码跑xHat几乎必报奇异或给出巨大的DCB值这是秩亏。原因很直观卫星DCB加接收机DCB任意把其中一个整体平移一个常数加进另一个残差完全不变。也就是说“所有卫星DCB整体1所有接收机DCB整体−1”这条方向在法矩阵里没有信息正规矩阵奇异。破解办法是给解空间补一个基准。最常用的是假设一个解算时段内所有参与观测的卫星DCB之和为零Σ DCB^s 0在法方程里增加一行或直接对A做行扩展。扩展方式如下% 在A和b底部追加卫星DCB零和约束 A(end1, nModelnRx1 : end) 1; b(end1, 1) 0; % 给约束行一个比观测小很多的标准差例如1e-6 sigma 1e-6; A(end, :) A(end, :) / sigma; b(end) b(end) / sigma;另一种做法是固定一个已知卫星的DCB例如选一颗健康度好的卫星用IGS发布的产品直接固定这样接收机DCB全部相对该卫星确定。但固定单星会把那颗卫星的DCB产品误差直接引入所有接收机稳定性不如约束所有卫星之和。工程上我一般首选零和约束只有处理单站数据或者卫星数据少时才改用固定单星。参数设置里还有几个和法方程尺度相关的项。观测权重按高度角给常见取sin^2(el)或1/sin²(el)把低高度角观测的误差放大倍数抵消掉。球谐阶数固定时不要一开始就上8阶先4阶跑通再逐级看残差是否有明显下降。正则化方面高阶系数方差不稳时就直接把法矩阵对角加一个1e-6~1e-8的小量虽然破坏严格无偏但能压住震荡。下面是一个我常用的参数表直接复现时可以先照这个设参数取值说明高度角截止10°~15°低高度角映射函数误差大多项式阶数2~3单站或小区域用球谐最高阶4~8测站密时取高阶穿刺点高度350 km也可按450 km平滑系数0.02载波相位平滑伪距卫星DCB约束零和秩亏基准正则化1e-7法矩阵对角加小量这组参数在小范围单站网和大范围IGS网里都跑得通。区别在观测量的数量单站解算时一颗卫星的弧段很短DCB和模型系数的相关性会更高这时约束权重要适当加大。5. 高频避坑指南DCB解算跑飞、TEC跳变的五个真实原因与处理5.1 夜间与低纬度残差不收敛DCB解出来天天飘现象同一台接收机连续两天的DCB估计值相差几个纳秒白天结果看着还行夜间历元残差成片偏大。 原因夜间电离层TEC很低STEC观测值小DCB的相对占比变大加上夜间观测几何变化缓慢法方程对DCB的约束变弱。低纬度地区还有等离子体层和闪烁影响TEC本身不够平滑。 解决夜间弧段单独降低权重或把截止高度角从10°提到15°如果只是要DCB尽量用白天6点到22点的高高度角数据反演夜间数据留给验证用。残差若在固定某颗卫星上成片异常优先怀疑该卫星弧段的相位平滑问题而不是DCB变化。5.2 卫星DCB之和为零的基准在变化连续解出现系统跳变现象连续三天的卫星DCB解出来每颗卫星的DCB当天内部平滑但相邻天之间有整体平移且各颗星平移量相同。 原因正是零和约束自身造成的。零和约束没有把“绝对基准”钉死只是把一天内卫星DCB均值清零。如果某天参与解算的卫星集合不同或某颗卫星数据质量差整体均值会变化所有卫星DCB一起平移。 解决把连续解之间的卫星均值差作为整体偏差先求当天卫星DCB均值再对均值做时间平滑或直接对齐到IGS产品也可以把已知卫星的DCB做弱约束给IGS产品一个0.5 ns的容许量这样基准被软绑定又不会完全继承产品误差。5.3 高阶球谐出现东西向条带或负值现象球谐阶数从4提到8后TEC图出现规律东西向条带低纬区域出现局部负TEC物理上说不通。 原因高阶球谐的自由度增加把测站分布不均匀和观测噪声直接拟合进了模型。低纬地区穿刺点覆盖密集高阶项被噪声激发出现振荡。 解决第一步先检查穿刺点覆盖网格内观测少于一定数量就降阶第二步在法方程对角线加正则化重点对最高阶两阶的参数加大惩罚第三步用多项式模型对比一下如果多项式结果反而更平滑说明观测密度撑不起高阶球谐。5.4 载波相位平滑被周跳污染现象某颗卫星的STEC时间序列出现一个台阶前后两段波形都正常中间跳变残差对应历元出现团状异常。 原因Hatch滤波假设相位历元间差分无整周跳变。周跳发生时L(k)−L(k−1)跳了整周平滑器把伪距改正带偏输出在数分钟内被污染。没有做周跳探测就平滑是所有DCB估计翻车里最隐蔽的一个。 解决平滑前用GF组合或MW组合探测周跳发现跳变就把平滑窗口重置为当前历元。我通常的做法是用L1−L2的多路径组合做粗检测再用电离层残差变化率做精检测周跳处直接重新初始化平滑器。5.5 SP3插值在边缘历元产生假坐标现象测站全天观测数据里每天开头半小时和结尾半小时的DCB解算残差系统性偏大中间时段完全正常。 原因第3章里写的拉格朗日11点插值在数据两端只能取到一侧的点窗口钳位后基函数振荡卫星坐标偏差达到分米量级直接污染高度角和穿刺点。 解决解算时丢弃一天中卫星星历覆盖边缘的时段更稳妥的是用前后两天的SP3文件拼成超长序列再做插值让目标历元始终落在窗口中间。插值误差不直接检查坐标可以在解算后用无电离层组合相位残差做间接验证。6. 残差驱动迭代与验后验证把单天解练到和IGS产品对得上的技巧前面五章跑完能拿到一组可看的DCB和TEC但离稳定可发布还差两步残差清洗和验后质量检查。我现在的习惯做法是两轮加权最小二乘第一步用原始观测跑一次平差保存每个观测的残差第二步按残差绝对值排序剔除超过3倍RMS离差的观测并给剩余观测重新按高度角加权第三步用清洗后的观测重复解算再检查卫星DCB与IGS产品差值的RMS若小于0.5 ns就算稳定第四步用解出的参数回代VTEC按1°×1°网格化输出与同一时段的CODE GIM网格做差RMS在2~4 TECU说明链路基本正常。验证时不只盯DCB还要看TEC图有没有物理上可疑的负值或条带。如果TEC与GIM相差大但DCB差值很小问题多半出在球谐阶数和穿刺点映射上反过来DCB差值大优先怀疑某个接收机的平滑伪距序列还有周跳残留。一个实用的收敛技巧是第一轮用sin²(el)定权第二轮改用残差方差定权测站级重新缩放。伪距噪声水平在测站间差异很大单靠高度角定权会把噪声大的站顶得太高。残差定权后低质量测站的影响被自动削弱DCB的稳定性明显提升。这套流程里最容易忽略的是数据质量标签。RINEX文件里的天线类型、接收机类型和测站概略坐标会直接影响接收机坐标计算坐标差几米穿刺点位置就差几千米低高度角观测全被污染。我早年在这上面吃过亏现在所有测站坐标都会和IGS SINEX比对一遍再开工。希望帮到你。本文还有配套的精品资源点击获取
