基于 PPP-AR 的大气提取：三种坐标策略、数学原理与工程实现

约定：本文的 Fixed-coordinate 是“测站坐标作为已知量”，ambiguity-fixed 是“整周模糊度固定”，两者完全不同。长度单位统一为 m，方差单位为 m²。三种策略均以相同的观测、轨钟/偏差产品和 PPP-AR 整数定义为前提。

1. 核心认识：已知精密坐标很有用，但并非必须

GNSS 载波相位、伪距中既包含测站与卫星的几何距离，也包含对流层传播延迟。PPP-AR 利用精密轨道、卫星钟差及相容的码/相位偏差产品处理这些观测；在模型充分、数据质量良好且整数固定正确时，可得到稳定的天顶总延迟（ZTD）或天顶湿延迟（ZWD）。普通浮点 PPP 同样能估计对流层参数，PPP-AR 是改善估计的手段，不是提取 ZTD 的必要前提。[1][2][4] PPP-AR 也不会自动消除对流层误差；整数固定错误反而可能污染结果。

测站高程误差与对流层湿延迟在观测方程中容易相互吸收，这是大气提取时需要特别处理坐标的核心原因。把可靠的精密坐标固定，可以降低这种参数耦合；但若固定坐标本身有误，得到的 ZWD 可能形式方差较小、真实偏差却较大。

三种策略分别是：

• Static PPP-AR：估计整段静态弧段共享的三维坐标，并估计随历元变化的大气和钟差。
• Fixed-coordinate PPP-AR：采用已知三维模型坐标，不再估计坐标改正；其坐标误差不会在解算器的坐标状态中体现。
• Coordinate-constrained PPP-AR：三维坐标仍为未知量，使用精密坐标及可信的协方差施加软约束。

注意：第三种策略的约束方差越小，坐标越难偏离先验；它不能仅凭“给三个约束式设置极小方差”就自动识别短时形变。

2. 先弄清大气提取到底提取什么

2.1 ZTD、ZHD、ZWD 与斜路径延迟

$$\boxed{ZTD(t)=ZHD(t)+ZWD(t)}$$

$$T^s(t)=m_h(e^s)ZHD(t)+m_w(e^s)ZWD(t)+m_g(e^s)[G_N(t)\cos A^s+G_E(t)\sin A^s]$$

符号	物理含义	单位
$t,s$	观测历元和卫星编号	—
$ZTD$	天顶总延迟	m
$ZHD$	天顶静力（通常称干）延迟，常由气压及模型给出	m
$ZWD$	天顶湿延迟，常作为未知参数或先验改正量估计	m
$T^s$	卫星 $s$ 的斜路径中性大气延迟	m
$e^s,A^s$	卫星高度角和方位角；这里方位角从北顺时针	rad（使用三角函数时）
$m_h,m_w,m_g$	静力、湿及梯度映射函数	无量纲
$G_N,G_E$	北向、东向水平梯度参数	m

梯度项不是每套系统都估计；若不估计则省略。梯度映射函数与梯度参数的定义必须成套使用，不能在不同软件或产品间直接混搭。[3]

典型大气产品生成链路是：精密改正及观测预处理 → PPP/PPP-AR 参数估计 → 恢复 ZTD/ZWD → 质量控制 → 区域建模/格网播发。Cui 等（2022）给出了参考站大气提取、区域建模并用于用户 PPP-AR 增强的实例。[7]

2.2 “估计 ZWD”不一定是同一个参数化

状态设计	实际估计量	输出 ZTD 的正确方式
直接估计完整 ZWD	$\widehat{ZWD}$	$ZH{D}_0+\widehat{ZWD}$
估计 ZWD 先验改正	$\widehat{\delta ZWD}$	$ZH{D}_0+ZW{D}_0+\widehat{\delta ZWD}$
直接估计 ZTD 或其改正	$\widehat{ZTD}$ 或 $\widehat{\delta ZTD}$	$\widehat{ZTD}$ 或 $ZT{D}_0+\widehat{\delta ZTD}$

其中下标 $0$ 表示采用的先验值。千万不要把 $\delta ZWD$ 当成完整 ZWD，也不要在已估计完整 ZTD 的情况下再加 ZHD。

ZHD 背景错误、映射函数失配、天线误差或误固定，均可能被 ZWD 状态部分吸收。因此，“PPP 提取的 ZWD”是在所选观测模型下估计的参数，不等于未经模型误差影响的纯净水汽信号。若要从 ZWD 进一步反演可降水量（PWV），还需加权平均气温等气象信息，此处不展开。

2.3 UC/IF 与 PPP-AR 的边界

• UC-PPP 可以联合估计逐星电离层斜延迟，但必须处理码/相位偏差和基准（datum）问题。
• IF-PPP 通过组合消去一阶电离层项；未经相应整数恢复变换的 IF 相位模糊度一般不是整数，不能直接四舍五入固定。[1][4]
• PPP-AR 依赖与所选频点、星座和钟差基准一致的 OSB/UPD 等产品。只有在整周模糊度确实正确固定时，才有条件利用整数约束增强大气估计。错误固定可使 ZTD/ZWD 受到严重污染。

本文重点是中性大气 ZTD/ZWD，不讨论由 UC 模型提取电离层产品的独立可观测性问题。

3. 坐标为什么会影响 ZWD？

3.1 从线性化观测方程看

把精密轨钟、相对论、天线、潮汐及其他已建模项扣除后，一条线性化观测可写成：

$${\ell }_{j,k}=-{\mathbf{u}}_{j,k}^T\delta \mathbf{r}+\delta b_{r,k}+m_{w,j,k}\delta z_k+{\mathbf{h}}_{q,j,k}^T\delta \mathbf{q}+{\epsilon }_{j,k}$$

符号	含义	单位
$j,k$	观测行编号（卫星/频点/观测类型）和历元编号	—
${\ell }_{j,k}$	观测值减去线性化点计算值（OMC）	m
${\mathbf{u}}_{j,k}$	从接收机指向卫星的 ECEF 单位视线向量	无量纲
$\delta \mathbf{r}$	测站 ECEF 坐标改正 $[\delta X,\delta Y,\delta Z{]}^T$	m
$\delta b_{r,k}$	接收机钟差改正的距离形式，即 $c\delta t_{r,k}$	m
$c,\delta t_{r,k}$	真空光速及接收机钟差改正	m/s、s
$m_{w,j,k}$	该卫星方向的湿映射函数	无量纲
$\delta z_k$	选定 ZWD 参数化下的湿延迟改正	m
$\delta \mathbf{q}$	其余待估参数（如模糊度、系统间偏差、梯度、电离层）改正	随参数定义
${\mathbf{h}}_{q,j,k}$	其余参数在该观测行的偏导列向量	对应单位的倒数
${\epsilon }_{j,k}$	噪声和未建模误差	m

此式为理解三种策略的一阶简化模型。坐标变化还会改变卫星高度角、映射函数、ZHD（特别是使用模型气压时）及天线改正；完整 PPP 代码应对实际采用的观测模型求导。若使用本教程的简化雅可比，必须将这些较小的导数视为已忽略的近似，而不是认为它们严格为零。

3.2 为什么重点担心高程？

测站向上或向下移动，会改变各卫星的几何距离；ZWD 增减则通过湿映射函数改变各方向的斜路径延迟。二者对不同高度角观测的影响有可分辨之处，但在可见卫星少、低仰角数据少、观测时间短或权模型不合理时，很容易出现较强相关性。并非所有测站都固定存在相同的“高程到 ZWD”转换系数。[2][8]

固定坐标去掉的是一组待估自由度，而不是把所有对流层建模误差消除。 接收机钟差与 ZWD、梯度和模糊度之间仍可能存在相关性。

4. 三种坐标处理策略

为便于统一表述，定义测站的历元模型坐标为 ${\mathbf{r}}_{0,k}$。它可以由参考历元坐标、速度、已知潮汐/加载改正等构造，并不要求每个历元具有完全相同的 ECEF 数值。

4.1 Static PPP-AR：坐标作为静态未知量

$${\mathbf{r}}_k={\mathbf{r}}_{0,k}+\delta \mathbf{r},\delta \mathbf{r}=\text{同一处理弧段共享的未知常量}$$

$${\mathbf{x}}_k=[\delta X,\delta Y,\delta Z,\delta z_k,b_{r,k},{\mathbf{q}}_k^T{]}^T$$

其中 $\delta z_k$ 是当前历元或当前时间结点的大气状态，$b_{r,k}$ 是接收机钟差距离项，${\mathbf{q}}_k$ 汇集其余参数；这些参数的实际时间相关性由具体解算器定义。

这类模型适合先定坐标、再提取大气或同时进行坐标独立检核。其限制是短弧段坐标—ZWD 相关性可能强。这里的 Static 不是每一历元估计三个全新的 XYZ，更不是完全不做站点位移改正。仅将外部坐标用作非约束的初值仍属于 Static；若把其精度作为先验加入求解，就已属于 Constrained。[1]

4.2 Fixed-coordinate PPP-AR：坐标作为精确已知量

$${\mathbf{r}}_k={\mathbf{r}}_{\mathrm{fix},k},\delta {\mathbf{r}}_k\equiv 0$$

按固定坐标处重新计算的观测方程，不再包含 $\delta X,\delta Y,\delta Z$ 三列。它有利于减少位置—大气耦合，但先验坐标错误可能进入 ZWD、钟差、梯度及残差。Pan 与 Guo（2018）采用精密 IGS 站坐标开展了多系统实时对流层延迟提取。[2]

“完全相信坐标”是估计器的处理假设，并不是声称地面没有潮汐、台站速度或天线更换。 生产解中仍要维持一致的逐历元测站坐标模型。

4.3 Coordinate-constrained PPP-AR：坐标未知且带先验

在当前线性化坐标 ${\mathbf{r}}_{\mathrm{lin}}$ 下，待估改正为 $\delta \mathbf{r}$；外部先验为 ${\mathbf{r}}_{\mathrm{prior}}$。其一般形式是：

$$\underset{{\mathbf{d}}_r}{\underbrace{{\mathbf{r}}_{\mathrm{prior}}-{\mathbf{r}}_{\mathrm{lin}}}}={\mathbf{I}}_3\delta \mathbf{r}+{\mathbf{\eta }}_r,\mathrm{Cov}({\mathbf{\eta }}_r)={\mathbf{C}}_r$$

其中：

• ${\mathbf{d}}_r$：先验坐标相对线性化坐标的差值（m）；
• ${\mathbf{I}}_3$：三阶单位矩阵；
• ${\mathbf{\eta }}_r$：先验位置误差（m），一般假设均值为零；
• ${\mathbf{C}}_r$：$3\times 3$ 先验坐标协方差（m²），必须正定才能按普通高斯先验做白化。

只有在 ${\mathbf{r}}_{\mathrm{prior}}={\mathbf{r}}_{\mathrm{lin}}$ 时，伪观测右端才是 $0$。 这对迭代最小二乘、更新线性化点或使用不同坐标基准的代码尤为重要。

该模型允许数据提出坐标修正，但其幅度受到先验协方差约束。一个时段共享的静态坐标修正只能描述整段平均偏移，不能直接追踪真实的逐历元形变。若线性化点逐历元包含不同的已知位移改正，应先将先验和估计状态统一为同一参考历元下的常量坐标改正，而不是每个历元把略有不同的 XYZ 当作独立先验重复加入。

5. 三种策略优缺点（工程对比）

对比维度	Static PPP-AR	Fixed-coordinate PPP-AR	Coordinate-constrained PPP-AR
是否估计三维坐标	是，弧段共享	否	是，有先验约束
对外部精密坐标依赖	较低	最高	取决于先验协方差
坐标—ZWD 参数耦合	可能较强	去掉坐标待估项	随先验变强而减弱
ZWD 形式方差	通常较大	正确坐标条件下通常较小	在二者之间（符合第 6 章条件时）
ZWD 真实准确度	受几何与估计稳定性影响	最怕错误先验位置	受先验偏差与权重影响
坐标可诊断性	高	不能直接输出坐标修正	可输出受约束的坐标修正
可否估计时变形变	纯静态模型不能	外部时变模型可以表达，但不能自行估计	纯静态软约束也不能
估计自由度	较多	较少	较多，但额外加入先验信息
实施与运行	定坐标和大气联合估计	简单，关键是坐标准确性	需维护先验方差与相关性
推荐用途	测站坐标复核、无可靠先验	稳定 CORS 站产品生产	先验坐标有可信不确定度、做敏感性分析
核心工程风险	高程与 ZWD 弱可分辨	位置误差被大气等状态吸收	约束失配、重复计权、过强先验

表中的性能是特定线性统计模型下的理论趋势，不是对所有天气、站网或设备的实测排名。 表中“Fixed 形式方差较小”是假定固定坐标精确无误的条件结果。比较不同策略时，应使用相同的观测、改正产品、模糊度固定准则、ZWD 时间模型与质量控制规则。

6. 最小二乘如何严格联系三种策略？

6.1 同一线性化模型

$$\mathbf{y}=\mathbf{G}\delta \mathbf{r}+\mathbf{m}\delta z+\mathbf{B}\mathbf{\beta }+\mathbf{\epsilon },\mathrm{Cov}(\mathbf{\epsilon })=\mathbf{R}$$

其中：

• $\mathbf{y}$：观测减计算向量（m）；
• $\mathbf{G}$：测站三维坐标偏导矩阵（行数为观测数，三列）；
• $\mathbf{m}$：湿延迟参数的偏导列（通常为湿映射函数）；
• $\delta z$：待估 ZWD 改正（m），用于解释耦合的标量例子；
• $\mathbf{B}\mathbf{\beta }$：钟差、模糊度等其他未知参数的线性项；
• $\mathbf{\epsilon },\mathbf{R}$：观测误差及其协方差；
• $\mathbf{W}={\mathbf{R}}^{-1}$：观测权阵。

这里把 $\delta z$ 暂写成一个标量。真实 PPP 处理多个 ZWD 时间结点时，$\mathbf{m}$ 应扩展为矩阵；若 ZWD 具有随机游走、分段平滑或其他先验，须把相应时间约束也加入整体加权系统，才能讨论完整的后验协方差。

6.2 三种最小二乘目标

Static：

$$\underset{\delta \mathbf{r},\delta z,\mathbf{\beta }}{\min }\Vert \mathbf{y}-\mathbf{G}\delta \mathbf{r}-\mathbf{m}\delta z-\mathbf{B}\mathbf{\beta }{\Vert }_{\mathbf{W}}^2$$

Fixed：若线性化点与固定坐标不同，令 ${\mathbf{d}}_f={\mathbf{r}}_{\mathrm{fix}}-{\mathbf{r}}_{\mathrm{lin}}$，则：

$$\underset{\delta z,\mathbf{\beta }}{\min }\Vert \underset{{\mathbf{y}}_f}{\underbrace{\mathbf{y}-\mathbf{G}{\mathbf{d}}_f}}-\mathbf{m}\delta z-\mathbf{B}\mathbf{\beta }{\Vert }_{\mathbf{W}}^2$$

Constrained：

$$\underset{\delta \mathbf{r},\delta z,\mathbf{\beta }}{\min }\Vert \mathbf{y}-\mathbf{G}\delta \mathbf{r}-\mathbf{m}\delta z-\mathbf{B}\mathbf{\beta }{\Vert }_{\mathbf{W}}^2+\Vert \delta \mathbf{r}-{\mathbf{d}}_r{\Vert }_{{\mathbf{C}}_r^{-1}}^2$$

$\Vert \mathbf{v}{\Vert }_{\mathbf{W}}^2\equiv {\mathbf{v}}^T\mathbf{W}\mathbf{v}$ 为加权二次范数。Fixed 不是把坐标列删除就结束了：如果 ${\mathbf{d}}_f\ne 0$，还必须更新 OMC。 另一种正确方法是先在固定坐标处完整重新线性化，然后直接删除坐标未知量。

6.3 引入先验时的前提

上面的简单加权相加默认观测误差与外部坐标先验误差相互独立，且坐标先验仅施加一次。如果先验来自同一组观测、同一日 PPP 解，或共享了不可忽略的系统误差，需正确处理交叉相关或改用一致的联合估计，不能当作新的独立信息反复计权。软约束也不意味着“不相信先验”：它是明确规定先验允许偏离的概率尺度。

7. 固定坐标为什么可能降低 ZWD 形式方差？

7.1 必须先处理其余未知量

为避免忽略接收机钟差和模糊度等干扰参数，先从整体线性最小二乘信息矩阵中消去 $\mathbf{B}\mathbf{\beta }$。下列 ${\mathbf{W}}_{\perp }$ 表达式适用于这些干扰参数没有额外独立先验、${\mathbf{B}}^T\mathbf{W}\mathbf{B}$ 可逆的情形；存在钟差、模糊度或 ZWD 时间先验时，应先把对应约束加入完整法方程，再进行分块消元。

$${\mathbf{W}}_{\perp }=\mathbf{W}-\mathbf{W}\mathbf{B}({\mathbf{B}}^T\mathbf{W}\mathbf{B}{)}^{-1}{\mathbf{B}}^T\mathbf{W}$$

$$\mathbf{N}={\left[\begin{array}{cc}\mathbf{G}& \mathbf{m}\end{array}\right]}^T{\mathbf{W}}_{\perp }\left[\begin{array}{cc}\mathbf{G}& \mathbf{m}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{N}}_{rr}& {\mathbf{N}}_{rz}\\ {\mathbf{N}}_{zr}& N_{zz}\end{array}\right]$$

其中 ${\mathbf{W}}_{\perp }$ 是消去干扰参数后的残差信息阵（不必把它误称为原始观测协方差的逆）；${\mathbf{N}}_{rr}$ 是坐标信息块，${\mathbf{N}}_{rz}={\mathbf{N}}_{zr}^T$ 是坐标—ZWD 耦合信息，$N_{zz}$ 是 ZWD 信息标量。若相关矩阵不可逆，则须先处理参数基准和可观测性，不能直接套用逆矩阵。

7.2 Schur 补给出的形式方差

在相同观测数据、随机模型、线性化点、其余参数的先验与消元方式下，并且三种模式采用相同且正确的整数模糊度条件，以下是可辨识标量 ZWD 参数的形式方差：

$$\boxed{P_{z,\mathrm{static}}={\left(N_{zz}-{\mathbf{N}}_{zr}{\mathbf{N}}_{rr}^{-1}{\mathbf{N}}_{rz}\right)}^{-1}}$$

$$\boxed{P_{z,\mathrm{fixed}}=N_{zz}^{-1}}$$

$$\boxed{P_{z,\mathrm{constrained}}={\left[N_{zz}-{\mathbf{N}}_{zr}({\mathbf{N}}_{rr}+{\mathbf{C}}_r^{-1}{)}^{-1}{\mathbf{N}}_{rz}\right]}^{-1}}$$

这里 $P_z$ 表示湿延迟估计的形式方差（m²）；${\mathbf{C}}_r$ 是坐标先验协方差。对可辨识的正常模型，有：

$$P_{z,\mathrm{fixed}}\le P_{z,\mathrm{constrained}}\le P_{z,\mathrm{static}}$$

同时，${\mathbf{C}}_r\to 0$ 的极限对应 Fixed，先验权趋于零对应 Static（各极限需要法方程保持可辨识）。这只是条件形式方差的排序，绝不是三种策略下真实 ZWD RMS 的必然排序。 三次真实 PPP-AR 运行如果固定了不同的整数集合、出现误固定或使用不同线性化状态，就不能直接套用该不等式。

7.3 位置错误如何转移到 ZWD？

若真实坐标相对误固定值的差为 $\delta {\mathbf{r}}_{\mathrm{true}}$，但解算时强制坐标改正为零，则在上述一阶模型下有近似偏差：

$$\mathrm{bias}(\widehat{\delta z})\simeq ({\mathbf{m}}^T{\mathbf{W}}_{\perp }\mathbf{m}{)}^{-1}{\mathbf{m}}^T{\mathbf{W}}_{\perp }\mathbf{G}\delta {\mathbf{r}}_{\mathrm{true}}$$

其中 $\mathrm{bias}(\widehat{\delta z})=E[\widehat{\delta z}]-\delta z_{\mathrm{true}}$，$\delta z_{\mathrm{true}}$ 是真实的湿延迟改正。这里仅推导几何坐标误差项，实际还可能叠加高度相关 ZHD 误差、映射函数误差、潮汐改正遗漏等。此式也说明不存在一个对所有站点、所有卫星几何都固定的“高程误差 → ZWD 偏差”比例系数。

如果固定坐标仍有外部不确定度，解算器给出的“固定坐标条件下 ZWD 方差”通常没有自动包含这部分坐标不确定度；需要做误差传播或敏感性检验。若外部位置误差视为随机量、协方差为 ${\mathbf{C}}_r$，且与观测噪声独立，则一阶传播项为：

$${\mathrm{Var}}_r(\widehat{\delta z})\simeq {\mathbf{K}}_r{\mathbf{C}}_r{\mathbf{K}}_r^T,{\mathbf{K}}_r=({\mathbf{m}}^T{\mathbf{W}}_{\perp }\mathbf{m}{)}^{-1}{\mathbf{m}}^T{\mathbf{W}}_{\perp }\mathbf{G}$$

其中 ${\mathbf{K}}_r$ 是位置误差到 ZWD 估值的线性灵敏度矩阵（$1\times 3$）；它只表示坐标不确定度的一个传播分量，不能在存在相关性或参数重估条件变化时机械与其他方差相加。评估实际准确度还应考虑偏差平方项（MSE = 方差 + 偏差平方）。

8. 如何设置坐标软约束的协方差？

8.1 先验标准差不是调小就更好

在三个坐标分量互不相关的教学假设下：

$${\mathbf{C}}_r=\mathrm{diag}({\sigma }_X^2,{\sigma }_Y^2,{\sigma }_Z^2),{\mathbf{W}}_r={\mathbf{C}}_r^{-1}$$

假设的每轴先验标准差	约束强度	含义
1 mm	很强	仅当先验坐标、参考框架和未建模位移确实达到相应水平时才合理
5 mm	较强	允许小幅静态偏移，但不保证能检测逐历元形变
20 mm	相对较弱	先验容忍度较大，坐标—ZWD 耦合也可能增强
$\sigma \to 0$	无限强的极限	数学上逼近 Fixed
$\sigma \to \infty$	无外部先验的极限	数学上逼近 Static

这些数值仅用于理解权重的量级，不是 GNSS 参考站的通用推荐配置。真实 ${\mathbf{C}}_r$ 应根据坐标解的协方差、速度传播、参考框架差异、台站环境位移和观测相关性建立。先验的标准差是对真实不确定性的量化，不是为了让结果更平滑而调节的经验旋钮。

8.2 在 ENU 下设置更符合物理认知

若水平和垂向先验精度不同，可先在 ENU 中设置：

$${\mathbf{C}}_{ENU}=\mathrm{diag}({\sigma }_E^2,{\sigma }_N^2,{\sigma }_U^2),{\mathbf{C}}_{XYZ}={\mathbf{T}}^T{\mathbf{C}}_{ENU}\mathbf{T}$$

其中 $\mathbf{T}$ 是从 ECEF 到 ENU 的正交旋转矩阵，${\sigma }_E,{\sigma }_N,{\sigma }_U$ 是东、北、天向先验标准差（m）。ENU 协方差不能未经旋转就当作 XYZ 协方差使用。 若有完整先验协方差，应保留非对角相关项。

8.3 是否必须添加三个坐标约束方程？

不必须。 三条 XYZ 伪观测表示对三维位置施加联合先验；如果已有可靠的水平方向基准、但只希望控制或检验高程，可以只对某一方向添加约束。例如只约束局部 Up 方向时：

$$d_U={\mathbf{t}}_U^T\delta \mathbf{r}+{\eta }_U,\mathrm{Var}({\eta }_U)={\sigma }_U^2$$

其中 ${\mathbf{t}}_U^T$ 是 ECEF→ENU 旋转矩阵的第三行（从 XYZ 改正得到 U 改正），$d_U$ 为先验与线性化坐标之间的 Up 方向差（m），${\eta }_U$ 是 Up 先验误差，${\sigma }_U$ 是该方向先验标准差（m）。仅约束 U 不会自动消除 E/N 与其余参数的耦合；如果问题要求完整三维先验，应使用三维全协方差。

8.4 强约束是否等于限制了毫米级形变？

不等于。先验协方差 ${\mathbf{C}}_r$ 描述的是先验坐标的不确定度，不直接规定坐标随时间变化的物理规律。若 $\delta \mathbf{r}$ 在整段解算中共享，它只能估计一个静态偏移。即使 ${\sigma }_U$ 取 5 mm，也不会自动识别某一小时内突然出现、下一小时消失的高程形变。

9. 时变坐标：先验约束与过程噪声是两种不同的信息

如需估计短时真实位移，必须把坐标定义为时变状态，例如：

$$\delta {\mathbf{r}}_k=\delta {\mathbf{r}}_{k-1}+{\mathbf{w}}_{r,k},\mathrm{Cov}({\mathbf{w}}_{r,k})={\mathbf{Q}}_{r,k}$$

符号	含义
$k$	当前历元编号
$\delta {\mathbf{r}}_k$	当前历元相对参考模型坐标的三维修正（m）
${\mathbf{w}}_{r,k}$	历元间未知坐标变化（m）
${\mathbf{Q}}_{r,k}$	该时间间隔内允许坐标变化的过程噪声协方差（m²）

${\mathbf{C}}_r$ 决定“多相信先验位置”，${\mathbf{Q}}_{r,k}$ 决定“允许位置怎样随时间变化”。 这两个参数的职责不同。对于潮汐、台站长期速度、已知震后阶跃等，可预先改正的运动应优先采用明确的物理/分段模型，而非任由随机游走吸收。[6]

还要避免把同一坐标先验作为全新的独立观测每历元重复注入。对静态共享坐标，重复使用同一先验会虚假累积信息。若采用动态状态并逐历元使用位置观测，应建立其时间相关结构。

10. C++11 + Eigen：三种策略的最小二乘实现

下面代码只演示同一个线性化观测模型中的坐标策略差异，不包含轨钟改正、ZHD 计算、OSB/UPD 处理、整数固定及状态时间更新。因此它不是可以直接产出 ZTD 的完整 PPP-AR 程序。

与常见简化例子相比，这里特别处理了两个工程细节：

1. Fixed：即使固定坐标不等于线性化坐标，也要把对应观测改正移到右端；
2. Constrained：允许先验坐标不等于线性化坐标，并使用完整的 $3\times 3$ 先验协方差。

#include <Eigen/Dense>
#include <stdexcept>

enum class CoordMode { Static, Fixed, Constrained };

// H 的前三列依次为 dX,dY,dZ；其余为 dZWD、clock 等待估参数。
// y: 在同一个线性化点建立的 OMC (m)；obsStd: 各行观测标准差 (m)。
// fixedDelta: 固定坐标 - 线性化坐标 (m)，仅 Fixed 使用。
// priorDelta: 先验坐标 - 线性化坐标 (m)，仅 Constrained 使用。
// priorCov:   先验坐标的 XYZ 协方差 (m^2)，仅 Constrained 使用。
// 返回完整参数改正向量；Fixed 的前三个元素为指定坐标差。
Eigen::VectorXd solveCoordMode(
    const Eigen::MatrixXd& H, const Eigen::VectorXd& y,
    const Eigen::VectorXd& obsStd, CoordMode mode,
    const Eigen::Vector3d& fixedDelta,
    const Eigen::Vector3d& priorDelta,
    const Eigen::Matrix3d& priorCov)
{
    const Eigen::Index n = H.rows();
    const Eigen::Index p = H.cols();
    if (p < 4 || n < 1 || y.size() != n || obsStd.size() != n ||
        !H.allFinite() || !y.allFinite() || !obsStd.allFinite() ||
        (obsStd.array() <= 0.0).any())
        throw std::invalid_argument("invalid observation input");

    const bool fixed = (mode == CoordMode::Fixed);
    const bool constrained = (mode == CoordMode::Constrained);
    const Eigen::Index offset = fixed ? 3 : 0;
    const Eigen::Index q = p - offset;

    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n + (constrained ? 3 : 0), q);
    Eigen::VectorXd b = Eigen::VectorXd::Zero(A.rows());

    if (fixed && !fixedDelta.allFinite())
        throw std::invalid_argument("invalid fixed coordinates");

    Eigen::VectorXd rhs = y;
    if (fixed) rhs.noalias() -= H.leftCols(3) * fixedDelta;
    for (Eigen::Index i = 0; i < n; ++i) {
        A.row(i) = H.block(i, offset, 1, q) / obsStd(i);
        b(i) = rhs(i) / obsStd(i);
    }

    if (constrained) {
        if (!priorDelta.allFinite() || !priorCov.allFinite())
            throw std::invalid_argument("invalid coordinate prior");
        Eigen::LLT<Eigen::Matrix3d> llt(priorCov);
        if (llt.info() != Eigen::Success)
            throw std::invalid_argument("prior covariance must be positive definite");
        const Eigen::Matrix3d Linv =
            llt.matrixL().solve(Eigen::Matrix3d::Identity());
        A.bottomRows(3).leftCols(3) = Linv;
        b.tail(3) = Linv * priorDelta;
    }

    Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(A);
    if (qr.rank() < q)
        throw std::runtime_error("rank-deficient design matrix");

    const Eigen::VectorXd result = qr.solve(b);
    Eigen::VectorXd dx = Eigen::VectorXd::Zero(p);
    if (fixed) {
        dx.head(3) = fixedDelta;
        dx.tail(q) = result;
    } else {
        dx = result;
    }
    return dx;
}

代码对应的 LaTeX 公式（与实现紧邻，方便逐行核对）：

观测白化矩阵、观测和设计矩阵为：

$$\mathbf{D}=\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_n),{\mathbf{A}}_0={\mathbf{D}}^{-1}\mathbf{H},{\mathbf{b}}_0={\mathbf{D}}^{-1}\mathbf{y}$$

其中 $n$ 为观测数、${\sigma }_i$ 为第 $i$ 行观测标准差、$\mathbf{D}$ 单位 m；$\mathbf{H}$ 为设计矩阵，$\mathbf{y}$ 是 OMC 向量。

Fixed 模式删除坐标三列，并改正右端：

$${\mathbf{A}}_f={\mathbf{D}}^{-1}{\mathbf{H}}_{[:,4:p]},{\mathbf{b}}_f={\mathbf{D}}^{-1}(\mathbf{y}-\mathbf{G}{\mathbf{d}}_f)$$

这里 $\mathbf{G}={\mathbf{H}}_{[:,1:3]}$ 为坐标三列，${\mathbf{H}}_{[:,4:p]}$ 表示余下参数列，${\mathbf{d}}_f$ 为固定坐标相对线性化点的差值。列下标在数学表述中从 1 开始，C++ 的数组下标从 0 开始。

Constrained 模式采用 ${\mathbf{C}}_r=\mathbf{L}{\mathbf{L}}^T$（Cholesky 分解）以及 ${\mathbf{d}}_r={\mathbf{r}}_{\mathrm{prior}}-{\mathbf{r}}_{\mathrm{lin}}$：

$$\underset{{\mathbf{A}}_c}{\underbrace{\left[\begin{array}{c}{\mathbf{D}}^{-1}\mathbf{H}\\ {\mathbf{L}}^{-1}{\mathbf{S}}_r\end{array}\right]}}\delta \mathbf{x}\simeq \underset{{\mathbf{b}}_c}{\underbrace{\left[\begin{array}{c}{\mathbf{D}}^{-1}\mathbf{y}\\ {\mathbf{L}}^{-1}{\mathbf{d}}_r\end{array}\right]}},{\mathbf{S}}_r=[{\mathbf{I}}_{3}0]$$

这里 $\mathbf{L}$ 是下三角矩阵，${\mathbf{S}}_r$ 从完整参数改正中抽取前三个坐标分量；${\mathbf{A}}_c,{\mathbf{b}}_c$ 是增加坐标先验后的白化方程。Static 仅使用 ${\mathbf{A}}_0,{\mathbf{b}}_0$。

本代码对观测采用彼此独立的标准差，但对位置先验允许完整协方差。真实 PPP 的观测噪声相关、跨历元先验、模糊度固定、非线性迭代及偏差产品一致性须在外层处理。QR 秩检查也不能代替条件数、残差、周跳及误固定风险检验。特别注意：固定坐标时若使用的 $\mathbf{y}$ 已在固定坐标处重新计算，传入的 fixedDelta 应为零；不能重复施加几何改正。若固定坐标与线性化点差异较大，应在固定坐标处完整重算观测模型及偏导数，而不是依赖一次线性近似。代码还假设 priorCov 是对称正定的可信协方差，实用工程应先验证它的对称性及数值尺度。

11. 完整处理流程（Mermaid）

上图表达模块关系，不限定具体软件的严格执行顺序。实时 PPP-AR 通常在每历元循环执行估计、周跳处理、整数固定及状态质量控制。

12. 决定大气产品质量的其他关键因素

因素	为什么重要	工程检查点
坐标参考框架、参考历元与速度	与精密轨道不一致会形成几何误差	ITRF/IGS 框架、历元与站点速度是否统一
站点真实位移	潮汐/加载位移会进入观测	固体潮、极潮、海潮负荷按同一规范处理 [6]
天线参考点	标石、ARP、相位中心不是同一个位置	天线高、PCO、PCV、radome、ANTEX 与站点日志 [5]
ZHD 及压力	ZHD 的误差可被 ZWD 吸收	站点气压或背景模型、海拔与垂直气压外推
映射函数、低仰角与梯度	几何分离和误差传播都依赖它们	高度角截止、权函数、梯度状态及模型一致性
轨钟与 OSB/UPD	PPP-AR 整数特性依赖产品基准一致	星座、频点、信号、钟差与偏差产品配套 [4]
相位质量与模糊度固定	周跳/误固定可能造成 ZWD 台阶或尖峰	固定前后残差、固定率、误固定检验
大气状态时间模型	过强平滑可能丢失快速天气变化	ZWD 更新间隔、过程噪声和输出延迟
多系统钟差与 ISB	未估计偏差可能与大气互相吸收	接收机钟差、系统间偏差和信号间偏差配置
先验相关性	独立性假设错误会低估不确定度	先验坐标来源、先验使用次数、交叉协方差

IGS 天线工作组明确要求地面天线、卫星天线的绝对相位中心模型及参考框架保持一致；因此精密坐标不仅是“三个 XYZ 数字”。[5]

13. 如何验证三种策略，避免只比较形式标准差？

推荐采用同一测站、相同轨钟/OSB、相同截止高度角、相同 ZWD 时间模型进行对照。

评价对象	建议指标	避免的误判
ZTD/ZWD 实测准确度	相对可信外部产品的 Bias、STD、RMS	形式标准差小就认为真实误差小
时间连续性	历元差分、跳变率、断流后重收敛	平滑导致的迟滞被当成稳定
位置与约束	$\delta E,\delta N,\delta U$、先验归一化残差	先验错误被 ZWD 吸收而不自知
PPP-AR 整数固定	固定率、误固定风险、固定前后变化	把错误固定带来的小协方差当成功
区域大气产品	留站验证、空间交叉验证、边界/高程敏感性	训练站自评过于乐观
用户增强结果	独立用户 PPP-RTK 收敛时间和位置精度	站点 ZTD 的微小差异未必改善用户结果

一个有用的敏感性检查：在统一其他条件的前提下，比较 Static、Fixed 以及不同有依据先验标准差的 Constrained 结果；若 Fixed 与软约束的 ZTD 出现稳定差异，优先核查坐标框架、参考点、ZHD、潮汐、天线改正，再分析真实大气变化。

对比时应给出处理弧段和收敛/稳定段的定义，并区分“条件形式精度”“外部真值误差”和“用户增强效果”。必要时先用长弧段确定坐标，再在独立观测时段提取大气，以减少同源坐标先验与原始观测的相关性；若重复使用同一数据，须明确其并非独立验证。

14. 工程选型建议

• 坐标长期稳定、精密坐标与框架/位移模型可靠的 CORS 站：优先评估 Fixed-coordinate PPP-AR，并同时跑 Constrained 做独立敏感性测试。
• 已有精密坐标，但存在可量化的不确定度：采用 Constrained PPP-AR；给出 ${\mathbf{C}}_r$ 的来源，而不是为了得到“漂亮的 ZWD STD”任意减小方差。
• 无可靠精密先验，或正在重新确定测站坐标：采用 Static PPP-AR；适当增加观测弧长，并用独立坐标或后处理产品交叉检查。
• 确实存在逐历元形变或设备变动：先修正站点日志和天线/坐标模型；必要时引入分段坐标或时变状态模型，不应让一个静态软约束同时承担定坐标和形变识别的任务。

最终判断标准不是哪套方法让 ZWD 形式方差最小，而是谁能在一致的坐标与误差模型下输出实际偏差小、稳定、连续、可独立验证的大气产品。

15. 六个容易混淆的概念

1. 固定坐标 = 固定模糊度？ 不是。坐标固定控制的是位置未知量；PPP-AR 的模糊度固定控制的是整数约束，它们可任意组合。

2. Static PPP-AR 会忽略地球潮汐吗？ 不会。坐标状态可表示相对逐历元测站模型坐标的共享常量改正。

3. 三条 XYZ 软约束总是 0=δX,δY,δZ 吗？ 只有先验等于当前线性化点时才是零右端；一般应使用先验与当前线性化坐标之差。

4. 三条软约束方差越小，就越能允许毫米级形变？ 恰恰相反。方差越小越接近固定坐标；描述随时间变化的位移需要时间模型。

5. Fixed 的 ZWD 标准差最小，所以 ZTD 一定最准？ 不成立。固定错误坐标可以把位置误差转入 ZWD，而该误差未必体现在形式协方差中。

6. PPP-AR 大气提取必须估计完整 ZWD 吗？ 不一定。可以估计完整 ZWD、对 ZWD 先验的改正量或直接估计 ZTD；关键是状态定义与最后产品恢复公式一致。

参考文献

[1] ESA Navipedia. PPP Fundamentals. https://gssc.esa.int/navipedia/index.php/PPP_Fundamentals
用于 PPP 静态参数、接收机钟差、对流层和相位模糊度的基本参数化；不作为三种坐标策略优劣排序的直接实测证据。

[2] Pan, L., & Guo, F. (2018). Real-time tropospheric delay retrieval with GPS, GLONASS, Galileo and BDS data. Scientific Reports, 8, Article 17067. https://doi.org/10.1038/s41598-018-35155-3
用于多系统 PPP 实时 ZTD 提取、采用 IGS 精密站坐标固定位置的实例。

[3] Pacione, R., & Douša, J. (2020). SINEX_TRO: Solution (Software/Technique) INdependent EXchange Format for TROpospheric and meteorological parameters, Version 2.00 (December 2020). IGS format specification. https://files.igs.org/pub/data/format/sinex_tro_v2.00.pdf
用于 ZTD、ZHD、ZWD、梯度、单位和产品元数据的定义；此文是格式规范，不是本文坐标软约束算法的来源。

[4] International GNSS Service. Precise Point Positioning with Ambiguity Resolution (PPP-AR) Working Group. https://www.igs.org/wg/ppp-ar/
用于 PPP-AR 星钟/码相位偏差产品一致性与互操作性要求。

[5] International GNSS Service. Antenna Committee. https://www.igs.org/wg/antenna/
用于 ANTEX、天线相位中心改正、ARP/PCO/PCV 与参考框架一致性。

[6] Petit, G., & Luzum, B. (eds.). (2010). IERS Conventions (2010). IERS Technical Note No. 36, Chapter 7: Displacement of reference points. https://www.iers.org/iers/en/publications/technicalnotes/tn36
用于测站参考点位移模型、固体潮、极潮及相关加载改正的规范背景。

[7] Cui, B., Wang, J., Li, P., Ge, M., & Schuh, H. (2022). Modeling wide-area tropospheric delay corrections for fast PPP ambiguity resolution. GPS Solutions, 26, Article 56. https://doi.org/10.1007/s10291-022-01243-1
用于参考站 PPP-AR 大气提取、区域建模与用户增强的处理链路。

[8] Brack, A., & He, S. (2026). Inequality-constrained integer least-squares and its application to troposphere-bounded GNSS PPP. GPS Solutions, 30, Article 23 (online: 12 November 2025). https://doi.org/10.1007/s10291-025-01991-w
用于 ZWD 与 Up 参数相关性及合理物理区间约束的研究背景。该文约束的是ZWD 的不等式边界，并未提出本教程的 XYZ 高斯软约束三策略。

引用范围：上述论文和规范分别提供 PPP/PPP-AR 基础、固定精密坐标提取 ZTD 的研究实例、ZWD 区域增强、SINEX_TRO 产品规范、参考点位移与天线基准；三策略统一推导、Schur 补形式方差排序以及 C++11 示例为本文基于加权最小二乘构造的教学表达，不应误认为某篇论文提出了这套同名分类。