第一次看到“PPP”这三个字母十个人里有八个会先想到那个拨号上网时代的Point-to-Point Protocol。但在卫星导航这个圈子里PPP的另一个含义——Precise Point Positioning精密单点定位——才是真正值得耐心研究的硬骨头。手里一台双频GNSS接收机不架基准站不拉电台和4G差分链路单靠接收机自身把定位精度做到厘米级背后就是一套完整的PPP处理流程和一套不能出错的误差修正公式。这篇内容适合正在做GNSS数据处理、RTK转PPP、或者想搞懂“接收机输出的坐标到底怎么算出来”的朋友。我会把PPP从观测方程到误差修正公式逐项拆开再把从RINEX文件到最终坐标的完整流程走一遍最后夹带一些我实际处理数据时踩过的坑。1. 这个缩写经常被认错PPP到底是什么1.1 和拨号协议无关精密单点定位的前世今生精密单点定位的概念最早是1997年Zumberge等人在JPL提出来的核心思路非常简单直白既然差分定位能消掉卫星钟差、轨道误差、电离层延迟这些共性误差那我能不能不用基准站直接把误差一个一个“算明白”最终靠单台接收机实现精密绝对定位这个想法在当时非常超前因为要做到这一点需要一个前提——必须有高精度的卫星轨道和卫星钟差产品。在IGS国际GNSS服务还没成气候的上世纪九十年代这几乎是奢望。直到IGS的精密星历和精密钟差产品逐渐稳定发布PPP才从论文变成可用的工具。PPP的本质是绝对定位它输出的坐标直接定义在ITRF这样的全球参考框架下。和RTK依赖基准站的相对定位不同PPP没有作用距离的概念只要你头顶上有卫星、手上能拿到精密产品在全球任何地方都能得到一致的定位结果。1.2 在RTK垄断的高精度赛道里PPP凭什么分一杯羹很多人会问RTK精度又高、收敛又快为什么还要用PPP这个问题的答案取决于你身处什么场景。RTK的精度确实好水平能到1-2厘米但它有个天生弱点必须依靠基准站。短基线RTK作用距离一般不超过20公里网络RTK虽然覆盖广也需要运营商布设连续运行参考站网还要解决数据链路的问题。你在海上、沙漠、山区、边境线上作业时很可能根本没有CORS站可用这时候PPP就是唯一的答案。我给一个直观的对比对比项RTKPPP基准站需求必须或依赖CORS网不需要定位类型相对定位绝对定位参考框架随基准站而定ITRF等全球框架典型精度厘米级厘米级到分米级收敛时间秒级到分钟级十到几十分钟作业半径受基线长度限制全球覆盖主要误差处理差分消除逐项建模修正PPP的代价也很明确收敛慢。这是因为单台接收机无法通过差分把模糊度快速固定收敛过程受观测几何、卫星数量和误差模型质量影响很大。但如果你做的是长时段的静态监测比如滑坡监测、海平面监测、形变测量PPP完全可以胜任而且省掉了一堆基准站维护的麻烦。2. 建模思路决定成败观测方程与两种主流实现2.1 无电离层组合用加权差消除一阶电离层PPP中最经典的观测模型是无电离层组合Ionosphere-FreeIF它把双频观测值按频率平方的比例加权组合从而消掉一阶电离层延迟。这个模型从上世纪九十年代沿用至今现在很多业务化PPP软件仍然以它为主打。先说伪距和相位的原始观测方程。对频率$f_i$的伪距$P_i$忽略噪声的情况下可以写成$$ P_i \rho c(dt_r - dt^s) T \frac{I}{f_i^2} B_i $$载波相位观测值$L_i$则多了一个模糊度项并且电离层延迟符号相反$$ L_i \rho c(dt_r - dt^s) T - \frac{I}{f_i^2} \lambda_i N_i b_i $$其中$\rho$是卫星到接收机的几何距离$dt_r$是接收机钟差$dt^s$是卫星钟差$T$是对流层延迟$I/f_i^2$是电离层延迟$\lambda_i$和$N_i$分别是波长和整周模糊度。如果直接用$P_1$和$P_2$做线性组合取$$ P_{IF} \frac{f_1^2 P_1 - f_2^2 P_2}{f_1^2 - f_2^2} $$你会发现电离层延迟项被消掉了类似地载波相位也有$$ L_{IF} \frac{f_1^2 L_1 - f_2^2 L_2}{f_1^2 - f_2^2} $$组合后的PPP观测方程就变成了$$ P_{IF} \rho c(dt_r - dt^s) T \varepsilon_P $$$$ L_{IF} \rho c(dt_r - dt^s) T \lambda_{IF} N_{IF} \varepsilon_L $$这就是经典PPP的核心。无电离层组合最大的优点是模型简单不需要电离层参数直接用双频观测值就能把一阶电离层项“数学消除”。缺点也很明显组合后的观测噪声大约是原始观测值的两到三倍而且模糊度$N_{IF}$不再具有整数特性导致无法直接固定整周模糊度。2.2 非差非组合模型直接估计电离层延迟最近十年越来越多的PPP软件转向非差非组合UncombinedUC模型。这个模型不强行消电离层而是把倾斜电离层延迟当作参数直接估计。它保留了原始伪距和相位观测值模型写起来更贴近物理本质$$ P_i \rho c(dt_r - dt^s) T \gamma_i I_r \varepsilon_{P_i} $$$$ L_i \rho c(dt_r - dt^s) T - \gamma_i I_r \lambda_i N_i \varepsilon_{L_i} $$这里$\gamma_i$是电离层系数$\gamma_i f_1^2 / f_i^2$$I_r$是接收机天顶方向的一阶电离层延迟通常还要用投影函数把它投影到卫星视线方向。我个人的体会是非差非组合模型有四个实打实的好处观测噪声小因为不需要做组合放大三频四频卫导系统扩展起来非常自然直接多写几个方程就行除了定位坐标你还能顺带输出电离层延迟和对流层延迟这对气象和水汽反演很有用模糊度保持原始波长为后续固定整数模糊度提供了更好的基础。代价是状态向量里多了一堆电离层参数卡尔曼滤波的计算量变大而且电离层参数需要合理的随机游走模型否则估计出来会震荡得很厉害。2.3 参数估计里的隐藏变量无论你用IF模型还是UC模型最终都要落到参数估计上。PPP里最关键的部分不是那些显眼的坐标参数而是一堆“隐藏变量”接收机钟差、对流层天顶湿延迟、模糊度参数。以经典IF模型为例单系统PPP的状态向量一般可以写成$$ X \begin{bmatrix} x y z cdt_r ZWD N_{IF}^1 N_{IF}^2 \cdots N_{IF}^n \end{bmatrix}^T $$其中$(x,y,z)$是接收机坐标静态模式里这三个参数过程噪声设为零动态模式则要估计速度和加速度$dt_r$是接收机钟差通常当作白噪声或随机游走过程$ZWD$是对流层天顶湿延迟一般用随机游走建模后面的$N_{IF}$是各颗卫星的无电离层组合模糊度视为常数发生周跳时重置。估计方法主流是卡尔曼滤波。卡尔曼滤波的好处是实时性好适合逐历元处理而且状态预测和观测更新天然支持PPP中“坐标常量钟差白噪声模糊度常量”这种混合模型。如果做后处理也可以用最小二乘批处理但批处理遇到几百个小时的数据时矩阵规模很吓人卡尔曼滤波还是更顺手。3. 误差修正公式逐项拆解从卫星端到地面PPP“精密”二字全靠误差修正撑起来。少了任何一项坐标偏差就可能从厘米级变成分米级甚至米级。我把误差来源分成三类逐个公式拆开讲。3.1 卫星端轨道、钟差、天线相位中心与相对论精密轨道和钟差PPP的卫星位置不能靠广播星历必须用IGS的精密轨道产品。SP3格式的轨道文件通常每5分钟或15分钟一个节点你要用拉格朗日多项式或切比雪夫多项式插值到观测时刻。我自己的经验是用九阶或十阶拉格朗日插值用SP3节点以外的一点时间做收尾比如前后多外推一两个历元防止边缘震荡。精密钟差文件CLK格式给出卫星钟差通常30秒间隔。注意SP3和CLK必须来自同一产品系列而且时间基准要严格对齐否则会出现几厘米到几十厘米的系统差。卫星天线相位中心卫星天线相位中心改正要从igs14.atx天线文件中读取。这个文件里有卫星天线的PCO相位中心偏移和PCV相位中心变化。卫星PCO是星固坐标系下的一个固定矢量但在卫星姿态旋转时它对观测距离的投影会变化。近似的改正公式为$$ \Delta \rho_{ant,sat} \mathbf{e}{rec}^{sat} \cdot \mathbf{R}{body} \cdot \mathbf{r}_{PCO,sat} $$其中$\mathbf{e}{rec}^{sat}$是接收机到卫星的单位矢量$\mathbf{R}{body}$是星固坐标系相对惯性系的旋转矩阵$\mathbf{r}_{PCO,sat}$是卫星天线PCO在星固系下的坐标。PCV则通常按卫星天底角查表。相对论效应卫星钟在轨的高速运动和地球重力场会导致星钟频率偏移。IGS精密钟差产品已经吸收了大部分相对论效应但还遗留了一个周期性相对论项需要自行修正$$ \Delta \rho_{rel} - \frac{2}{c}(\mathbf{r}{sat} \cdot \mathbf{v}{sat}) $$这里的$\mathbf{r}{sat}$和$\mathbf{v}{sat}$是信号发射时刻卫星在地固系中的位置和速度矢量。这项改正虽然量级不大最大几厘米但如果不修正会在坐标时间序列里引入明显的长周期波动尤其在做静态监测时特别刺眼。相位缠绕卫星天线在运行中会不断旋转以保持对日定向这会导致载波相位观测值产生一个与几何无关的额外变化叫相位缠绕Phase Wind-up。动态定位时必须修正静态长时间观测如果忽略它模糊度估计会受污染。相位缠绕的改正量可以写成$$ \Delta \phi \delta \phi 2\pi N $$其中$\delta \phi$由天线体坐标系的等效偶极子方向和接收机天线方向决定$N$由连续跟踪的相位缠绕整数计数决定。大部分软件会自动处理但如果你自己写代码不要漏掉这一项。3.2 传播路径电离层高阶项与对流层湿延迟电离层IF组合已经把一阶电离层项消掉了但二阶和三阶电离层项仍然残留。对绝大多数动态定位场景来说残差只有毫米到厘米级可以不管。如果你做的是最高精度的静态PPP比如地震形变研究就要引入高阶电离层改正公式基于地磁场模型$$ \Delta I_2 \frac{7527c}{f^3} \int N_e B \cos\theta , ds $$二阶项的修正量大约在0-2厘米范围。大多数业务用户忽略它问题不大但你要知道它存在。对流层对流层延迟分干湿两项。干延迟比较稳定可以用Saastamoinen模型直接算$$ ZHD \frac{0.002277 \cdot P}{1 - 0.0026\cos(2\varphi) - 0.00028h} $$其中$P$是测站气压hPa$\varphi$是纬度$h$是海拔高度km。干延迟算出来以后不需要估计当作已知值改正掉就行。湿延迟$ZWD$比较难建模通常作为未知参数用卡尔曼滤波去估计。把天顶总延迟$ZTD$拆成干湿两部分$$ ZTD ZHD ZWD $$无论是干延迟还是湿延迟从天顶方向投影到卫星视线方向都需要映射函数。现在常用的是GMFGlobal Mapping Function$$ mf(\varepsilon) \frac{1 \frac{a}{1\frac{b}{1c}}}{\sin\varepsilon \frac{a}{\sin\varepsilon \frac{b}{\sin\varepsilon c}}} $$其中$\varepsilon$是卫星仰角$a$、$b$、$c$是和测站位置、年积日相关的经验系数。对流层延迟沿视线方向的总改正就是$$ T(\varepsilon) ZHD \cdot mf_{dry}(\varepsilon) ZWD \cdot mf_{wet}(\varepsilon) $$我刚开始做PPP时犯过一个错干湿延迟用了同一个映射函数结果在低仰角卫星上总是出现异常残差。后来换成分别用干湿映射函数情况立刻改善。3.3 地球物理效应固体潮、海潮与极潮很多人在做PPP时容易忽略地球物理效应因为它们不是“导航误差”而是“地球形变”。但PPP是绝对定位坐标框架本身就在随着固体潮、海潮变化所以必须修正。固体潮固体潮是日月引力导致的地壳弹性形变垂直方向最大可达30厘米水平方向也有几厘米是所有地球物理效应里量级最大的绝对不能忽略。IERS2010公约给出了计算公式写成简化形式是$$ \Delta \mathbf{r} \sum_{j1}^{2} \frac{GM_j}{GM} \frac{r^4}{R_j^3} \left[ (3l_2(\hat{R}_j \cdot \hat{r}))\hat{R}_j \left(\frac{3h_2}{2} - \frac{3l_2}{2}\right)(\hat{R}j \cdot \hat{r})^2 \hat{r} \right] \Delta \mathbf{r}{perm} $$其中$GM_j$和$GM$分别是日月和地球的引力常数$r$是地心到测站的距离$R_j$是日月到地心的距离$\hat{R}j$和$\hat{r}$是对应的单位矢量$h_2 \approx 0.6078$和$l_2 \approx 0.0847$是Love数和Shida数。后面还要补一项永久潮汐改正$\Delta\mathbf{r}{perm}$。实际程序里这个过程并不复杂逐历元算一下加到测站坐标上就行。但如果你不修正固体潮静态PPP的垂直分量时间序列会呈现明显的半日波和全日波解算精度会大打折扣。海洋负荷潮汐海洋负荷潮汐是潮汐水体对沿岸地壳的加载形变在离海近的测站上可达几厘米。它需要用海潮模型比如FES2014b、TPXO9输出负荷潮汐系数再结合潮汐谐波参数算改正。如果你的测站离海岸线几十公里以内强烈建议加上海潮改正。否则你会看到坐标时间序列里有一个明显的12小时左右周期分量而这个分量在对流层估计里是消不掉的。极潮极潮是由地球自转轴相对地壳的微小摆动引起的形变最大量级约2.5厘米发生在45度纬度附近。IERS2010给出了简便公式$$ \Delta r_{polar} -33\sin 2\varphi \cdot (m_1 \cos\lambda m_2 \sin\lambda) $$其中$m_1$、$m_2$是极移参数$\lambda$是经度。极潮改正量虽小但长期监测数据处理里加上后坐标残差会明显更平。3.4 与接收机相关的细节接收机天线相位中心同样有PCO和PCV可以通过igs14.atx文件查表改正。如果你用的接收机天线在atx文件里查不到至少要把PCO大致值输进去否则会造成厘米级系统偏差。还有一项容易忽略的是地球自转改正Sagnac效应。信号从卫星传到地面大约要0.07秒这个时间地球转过了大约30米所以计算几何距离时必须在卫星坐标上做地球自转改正或者在方程里显式加入Sagnac项$$ \Delta \rho_{sagnac} \frac{\omega_e}{c}(x^s y_r - y^s x_r) $$其中$\omega_e$是地球自转角速度$(x^s,y^s,z^s)$是卫星坐标$(x_r,y_r,z_r)$是接收机近似坐标。这一项不改正定位结果会直接偏出去几十米属于那种“一步错步步错”的基础操作。4. 从观测值到坐标的完整处理链讲了这么多公式现在把它们串成一条完整的处理链。一套标准的PPP处理流程大致分成四步。4.1 第一步数据采集与RINEX规范化接收机设置上PPP要求至少双频观测值最好能用三频或四频。采样率方面静态测量建议1-30秒动态测量按需求来但要注意数据量和后处理时间之间的平衡。天线架设时天线高的量测要精确到毫米级因为PPP是绝对定位天线高误差会直接进到坐标结果里。周围不能有遮挡物截止仰角建议先放到5度后期处理再根据数据质量调整。采集完原始数据后你需要把接收机厂商的私有格式转成标准RINEX格式。现在各家厂商的转换软件都做得比较友好但要注意RINEX版本。我强烈建议用RINEX 3.04以上版本因为它对多系统、多频率的支持要好得多字段也更规范。4.2 第二步精密产品下载与对齐这一步决定了你能不能解算出高精度坐标。IGS精密产品按时间延迟分几档产品类型轨道精度钟差精度时延IGS最终产品2.5 cm75 ps约14天IGR快速产品2.5 cm75 ps约1天IGU超快速产品3 cm预报5 cm约200 ps实时/预报实时产品RTS/SSR3 cm约100 ps实时广播做高精度事后静态处理首选IGS最终产品。做时效要求高的工程可以用IGR快速产品。实时PPP则要接收SSR改正数通过NTRIP协议从IGS RTS等实时服务商获取。下载完产品之后做“对齐”检查这一步经常被新手忽略检查观测数据的时标是GPST还是UTC和精密产品的时标是否一致检查SP3的参考框架和atx天线文件版本是否匹配比如用igs14的SP3就要配igs14.atx别混用igs08的产品和igs14的天线文件检查CLK文件的卫星钟差是否包含所有你需要解算的卫星系统。4.3 第三步逐历元解算和滤波调参完整流程可以描述成下面的循环读入一个历元的双频伪距和载波相位观测值用历元内所有卫星的观测值做一次标准单点定位SPP给出接收机位置的初始值根据SP3插值出每颗卫星在信号发射时刻的位置根据CLK插值出卫星钟差计算前面章节列出的所有误差修正项相对论、天线相位中心、相位缠绕、固体潮、海潮、极潮、Sagnac效应等形成观测方程代入滤波器的预测状态计算新息执行新息检验如果某颗卫星的残差超过预设阈值给它降权或剔除卡尔曼滤波更新得到当前历元的位置、钟差、ZWD和模糊度进入下一历元重复上述过程。卡尔曼滤波的调参是最考经验的部分。滤波器的过程噪声设置如果太紧参数会被“锁死”跟不上真实变化设置太松解算结果又会被观测噪声带得乱跳。我平时常用的初始参数经验值放在第6章的统一配置清单里。4.4 第四步质量评估与坐标输出处理完成后软件会输出一个坐标时间序列。静态PPP的最终坐标一般取收敛后一段时间比如最后30-60分钟的平均值用这段时间的标准差来评价内符合精度。如果有已知坐标可以计算外符合精度。评价一个PPP结果质量时我会同时看三样东西位置时间序列的收敛曲线是否平滑有没有反复震荡收敛后的标准差水平分量10毫米以内算不错垂直分量20毫米以内算合格模糊度残差和新息序列有没有系统性偏移如果有说明某个误差修正项没过关。5. 收敛时间这个老大难决定因素与压榨技巧5.1 收敛时间为什么跑不短PPP最让用户抓狂的问题就是收敛慢。传统单系统GPS静态PPP收敛到水平10厘米以内顺利的话要30到60分钟不顺利的话两三个小时也正常。收敛慢的根本原因是模糊度与接收机钟差之间的强相关性。单台接收机时接收机钟差和模糊度参数几乎线性相关卡尔曼滤波需要足够多的历元和卫星几何变化才能把它们逐步分离。此外无电离层组合的噪声被放大也拖慢了收敛速度。之前网上搞过不少“快速收敛”的噱头实际上PPP收敛就是个信息量积累的过程。几何分布差、卫星数少、观测噪声大的时候谁来了都没办法。5.2 我实测有效的几招经过大量数据实测我把有效缩短收敛时间的策略按见效程度排个序多系统融合这是效果最明显的一招。GPSBDSGalileoGLONASS四系统融合后卫星数轻松超过25颗几何构型大幅改善收敛时间能比单GPS缩短一半以上。我做一个静态测试时单GPS收敛用了大约42分钟四系统融合大约14分钟就进入了水平10厘米以内。多频观测BDS-3、Galileo都有三频甚至四频信号非差非组合模型可以直接利用多频观测值进一步降低噪声。三频观测值在高仰角下对模糊度解算的帮助非常直观。合理约束ZWD和电离层对流层湿延迟的随机游走谱密度设置要匹配实际气象变化。谱密度设太大ZWD参数会把位置误差吸收掉设太小又跟不上真实的水汽变化。一般静电力可以设在$3\times10^{-8}$到$3\times10^{-7} , \mathrm{m^2/s}$具体值要看数据处理时长。先验坐标约束如果做的是长期监测测站坐标可以先做一个小时左右的快速解把结果作为先验约束反馈进去能明显加速重新收敛。合理设置收敛判定条件不要用“解算到第几分钟就稳定”这种模糊说法。科学的办法是设定一个滑动窗口比如连续10个历元内N、E、U三个方向的坐标变化量都小于某个阈值常见取2-5厘米再判定为收敛。这样判定结果可重复也方便不同批次数据之间做对比。6. 实战中绕不开的坑踩坑记录与参数清单6.1 精密产品时延导致的时间错位我曾经处理一批IGU超快速产品的数据结果N方向整体偏了大约10厘米。排查了好久最后发现是因为IGU产品的钟差时间序列比SP3轨道时间序列晚了整整一个采样间隔插值时没有对齐产生了系统性偏差。这个坑的排查思路很简单在解算前随机选几颗卫星把插值得到的卫星轨道、卫星钟差和广播星历的对应值画在同一个时间轴上看看有没有时间偏移。如果发现错位修正产品的时间标签后再处理。6.2 周跳漏探引发的灾难性模糊度重收敛PPP中周跳一旦漏探等于把一颗卫星的模糊度估计值错误地带到后续所有历元。最坏的情况是收敛完成后一个漏探周跳又把坐标拉偏到分米级需要重新收敛白干半小时。周跳探测我强烈建议用双组合策略。GFGeometry-Free组合对周跳敏感$$ L_{GF} L_1 - L_2 $$它对电离层变化也敏感可能把电离层扰动误判成周跳。MWMelbourne-Wübbena组合则不受几何、钟差和电离层影响$$ L_{MW} \frac{1}{f_1 - f_2}(f_1 L_1 - f_2 L_2) - \frac{1}{f_1 f_2}(f_1 P_1 f_2 P_2) $$MW组合的缺点是伪距噪声大小周跳容易漏。把GF和MW联合起来用再配合TurboEdit算法基本能覆盖绝大部分周跳场景。而且要注意GF组合中的周跳$\Delta N_1$和$\Delta N_2$相等的时候GF组合检测不出来但MW组合几乎是这个特殊场景的唯一救星。6.3 一套可以照抄的PPP配置清单最后给出一套我常用的静态PPP配置参数供你当作初始值参考。不同场景可以在这个基础上微调参数项配置值观测值类型双频或无电离层组合或多频非差非组合截止仰角7度采样间隔30秒后处理静态卫星天线PCO/PCVigs14.atx接收机天线PCO/PCVigs14.atx对流层干延迟Saastamoinen模型 GMF干映射函数对流层湿延迟GMF湿映射函数 随机游走估计固体潮IERS2010海潮FES2014b近海必用极潮IERS2010相位缠绕修正相对论修正周期项卡尔曼滤波过程噪声ZWD$1\times10^{-7} , \mathrm{m^2/s}$接收机钟差过程噪声白噪声或$10^4 , \mathrm{m^2/s}$模糊度初始方差$10^6 , \mathrm{cycle^2}$这套配置不是拍脑袋写的都是实际处理数据过程中一步步调出来的经验值。你在自己数据上跑出来的最优过程噪声可能和我略有差异但至少能保证一个不错的起点。做PPP数据处理三年多下来我最深的体会是这套技术不神秘但也不简单。每一项误差修正单独看都不难难的是把它们全部串起来并在实际数据中排查出那一个让你结果偏几厘米的隐蔽错误。尤其是固体潮、海潮这类地球物理改正看起来和“定位”没什么关系但在静态高精度场景里它们恰恰是决定垂直分量能不能进2厘米的关键。所以如果你正在调PPP解算结果先别急着怀疑算法和软件把你的误差修正项逐项过一遍多数问题都在那里。