[学习]RTKLib详解：ppp.c与ppp_ar.c

RTKLib详解：ppp.c与ppp_ar.c

本文是 RTKLlib详解 系列文章的一篇，目前该系列文章还在持续总结写作中，以发表的如下，有兴趣的可以翻阅。

[学习] RTKlib详解：功能、工具与源码结构解析
[学习]RTKLib详解：pntpos.c与postpos.c
[学习]RTKLib详解：rtkcmn.c与rtkpos.c
[学习]RTKLib详解：ppp.c与ppp_ar.c

---

Part A： ppp.c

---

一、整体作用与工作流程

该代码实现了GNSS精密单点定位(PPP)算法，基于非差相位和伪距观测值进行卡尔曼滤波状态估计，该代码完整实现了PPP算法的所有核心环节，通过模块化设计将各个误差源独立处理，便于维护和扩展。主要流程包括：

1. 状态时间更新：通过udstate_ppp更新接收机位置、钟差、对流层参数和相位偏差
2. 观测方程构建：res_ppp计算理论距离与观测值的残差
3. 卡尔曼滤波更新：通过filter函数进行状态估计更新
4. 模糊度固定：pppamb尝试进行整周模糊度固定

具体的函数调用关系如下：

---

二、核心函数说明

1. pppos

void pppos(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)

• 功能：PPP主处理函数
• 输入：
  • rtk: RTK控制结构体
  • obs: 观测值数组
  • n: 观测值数量
  • nav: 导航数据
• 处理流程：
  1. 初始化状态向量
  2. 计算卫星位置和钟差
  3. 排除日食卫星观测值
  4. 多次迭代执行残差计算和滤波更新
  5. 更新解算状态和协方差矩阵

2. res_ppp

int res_ppp(int iter, const obsd_t *obs, int n, const double *rs, const double *dts, const double *vare, const int *svh,const nav_t *nav, const double *x, rtk_t *rtk, double *v, double *H, double *R, double *azel)

• 功能：计算观测残差和设计矩阵
• 数学原理：
  • 几何距离计算：$$r=\sqrt{(x_i^s-x{)}^2+(y_i^s-y{)}^2+(z_i^s-z{)}^2}$$
  • 相位观测方程：$$\varphi =\lambda N+r+c(d{t}^r-d{t}^s)+T+I+ϵ$$
• 输出：
  • v: 残差向量
  • H: 设计矩阵
  • R: 观测噪声协方差

3. tide_solid

static void tide_solid(const double *rsun, const double *rmoon, const double *pos, const double *E, double gmst, int opt, double *dr)

• 功能：计算固体地球潮汐位移
• 数学模型：
  $$u=K_2[H_2(\frac{3}{2}{\cos }^2\theta -\frac{1}{2})-3L_2{a}^2]+K_3[H_3(...)]$$
  其中 $K_n=\frac{GM}{r^{n+1}}{R}^n$
• 输入：
  • rsun/rmoon: 日月位置向量
  • pos: 测站坐标
  • E: 坐标转换矩阵

4. prectrop

static double prectrop(gtime_t time, const double *pos, const double *azel,const prcopt_t *opt, const double *x, double *dtdx,double *var)

• 功能：精密对流层延迟模型
• 数学模型：
  $$T=m_h(\theta )ZHD+m_w(\theta )(ZWD+G_n\cot \theta \cos \alpha +G_e\cot \theta \sin \alpha )$$
• 参数：
  • m_h/m_w: 干湿映射函数
  • ZHD: 天顶静力学延迟
  • G_n/G_e: 水平梯度参数

5. corrmeas

static int corrmeas(const obsd_t *obs, const nav_t *nav, const double *pos,const double *azel, const prcopt_t *opt,const double *dantr, const double *dants, double phw,double *meas, double *var, int *brk)

• 功能：电离层和天线校正
• 处理流程：
  1. 无电离层组合计算
  2. 天线相位中心校正
  3. 风化改正
  4. DCB偏差校正

6. udbias_ppp

static void udbias_ppp(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)

• 功能：相位偏差状态更新
• 处理步骤：
  1. LLI检测周跳
  2. 几何自由组合检测周跳
  3. 相位偏差初始化
  4. 相位-码一致性校正

---

三、数学原理补充

卡尔曼滤波更新公式：

1. 预测：

   • 状态转移：${\hat{x}}_{k|k-1}^{-}={\Phi }_{k-1}{x}_{k-1}$
   • 协方差预测：$P_{k|k-1}^{-}={\Phi }_{k-1}{P}_{k-1}{\Phi }_{k-1}^T+Q_{k-1}$

2. 更新：

   • 卡尔曼增益：$K_k=P_{k|k-1}^{-}{H}_k^T(H_k{P}_{k|k-1}^{-}{H}_k^T+R_k{)}^{-1}$
   • 状态更新：${\hat{x}}_k={\hat{x}}_{k|k-1}^{-}+K_k(z_k-H_k{\hat{x}}_{k|k-1}^{-})$
   • 协方差更新：$P_k=(I-K_k{H}_k)P_{k|k-1}^{-}$

Saastamoinen对流层模型：
$$ZHD=0.002277\frac{P}{1-0.00266\cos (2\varphi )-0.00028H}$$
其中P为地面气压，φ为纬度，H为高程

---

四、代码特点

1. 支持多系统(GPS/GLO/GAL/SBS)
2. 提供多种电离层处理模式(广播星历/STEC/无模型)
3. 包含完整的误差修正模型(潮汐/对流层/天线偏差)
4. 实现了PPP-AR模糊度固定功能
5. 支持动态/静态模式切换

---

Part B: ppp_ar.c

---

一、整体作用与工作流程分析

这段代码实现了精密单点定位（PPP）中的整数模糊度解析（AR）功能，主要目标是通过宽巷（Wide-Lane, WL）和窄巷（Narrow-Lane, NL）组合观测值，固定载波相位模糊度以提高定位精度。整体工作流程如下：

1. LC（线性组合）观测值预处理：计算三频观测值的线性组合（如宽巷、窄巷组合），并进行噪声统计。
2. 宽巷模糊度固定：通过双频组合的宽巷模糊度解析，利用几何无关模型和置信度验证。
3. 窄巷模糊度固定：基于宽巷结果，采用四舍五入法（ROUND）或整数最小二乘法（ILS）进一步固定窄巷模糊度。
4. 状态更新：将固定后的模糊度作为约束条件，更新系统状态和协方差矩阵。

详细函数调用关系如下所示。

---

二、函数功能与参数说明

1. lam_LC(i, j, k)

• 功能：计算线性组合的波长（m）。
• 输入：
  • i, j, k：频率系数（如 1,-1,0 表示 L1-L2 组合）。
• 输出：波长 ${\lambda }_{\text{LC}}=\frac{c}{i{f}_1+j{f}_2+k{f}_5}$。

2. L_LC(i, j, k, L)

• 功能：计算载波相位线性组合（m）。
• 输入：
  • L[0], L[1], L[2]：L1/L2/L5 载波相位观测值（周）。
• 输出：组合值 $L_{\text{LC}}=\frac{i{f}_1{L}_1+j{f}_2{L}_2+k{f}_5{L}_5}{i{f}_1+j{f}_2+k{f}_5}$。

3. P_LC(i, j, k, P)

• 功能：计算伪距线性组合（m）。
• 输入：
  • P[0], P[1], P[2]：L1/L2/L5 伪距观测值（m）。
• 输出：组合值 $P_{\text{LC}}=\frac{i{f}_1{P}_1+j{f}_2{P}_2+k{f}_5{P}_5}{i{f}_1+j{f}_2+k{f}_5}$。

4. var_LC(i, j, k, sig)

• 功能：计算线性组合的噪声方差（m²）。
• 输入：
  • sig：原始观测值噪声（m）。
• 输出：方差 ${\sigma }_{\text{LC}}^2=\frac{i^2{f}_1^2+j^2{f}_2^2+k^2{f}_5^2}{(i{f}_1+j{f}_2+k{f}_5{)}^2}{\sigma }^2$。

5. conffunc(N, B, sig)

• 功能：计算模糊度整数假设的置信度。
• 输入：
  • N：候选整数模糊度。
  • B：浮点解模糊度。
  • sig：模糊度标准差。
• 输出：置信度$P={\sum }_{i=1}^7\left[\text{erfc}\left(\frac{i-x}{\sqrt{2}\sigma }\right)-\text{erfc}\left(\frac{i+x}{\sqrt{2}\sigma }\right)\right]$，其中 $x=|B-N|$。

6. average_LC(rtk, obs, n, nav, azel)

• 功能：对LC观测值进行滑动窗口平均。
• 输入：
  • rtk：RTK状态结构体。
  • obs：观测数据数组。
  • azel：卫星方位角仰角数组。
• 输出：更新 rtk->ambc 中的LC均值和方差。

7. fix_amb_WL(rtk, nav, sat1, sat2, NW)

• 功能：固定宽巷模糊度。

• 输入：

  • sat1, sat2：卫星编号。

• 输出：

  • NW：固定后的整数宽巷模糊度。

• 公式：
  $$B_W=\frac{\Delta L{C}_0^{(1)}-\Delta L{C}_0^{(2)}}{{\lambda }_{WL}}+\delta b_{WL}^{(1)}-\delta b_{WL}^{(2)}\text{(无REV\_WL\_FCB)}$$

  其中 $\delta b_{WL}$ 为宽巷小数周偏差（FCB）。

8. fix_amb_ROUND(rtk, sat1, sat2, NW, n)

• 功能：通过四舍五入法固定窄巷模糊度。
• 输入：
  • NW：宽巷模糊度数组。
• 输出：更新状态向量 rtk->x 和协方差矩阵 rtk->P。

9. fix_amb_ILS(rtk, sat1, sat2, NW, n)

• 功能：通过整数最小二乘（ILS）固定窄巷模糊度。
• 数学原理：
  定义窄巷模糊度 $B_{NL}=\frac{c_1{\lambda }_1{N}_1+c_2{\lambda }_2{N}_2}{{\lambda }_{NL}}$，构造变换矩阵 $ D $，求解整数最小二乘问题：
  $$\underset{N\in {\mathbb{Z}}^n}{\min }\Vert D^{-1}(B_{\text{float}}-N){\Vert }_Q^2$$
  使用 LAMBDA 算法降维后进行比值检验（ratio-test）。

10. pppamb(rtk, obs, n, nav, azel)

• 主函数：协调整个模糊度解析流程。
• 流程：
  1. 计算LC观测值均值。
  2. 枚举卫星对，固定宽巷模糊度。
  3. 根据 modear 选择ROUND或ILS方法固定窄巷模糊度。

---

三、数学推导补充

1.宽巷与窄巷组合

• 宽巷组合（如 L1-L2）：
  $${\lambda }_{WL}=\frac{c}{f_1-f_2}\approx 0.86\text{m}$$
• 窄巷组合（如 L1+L2）：
  $${\lambda }_{NL}=\frac{c}{f_1+f_2}\approx 0.11\text{m}$$

2.整数最小二乘（ILS）

定义模糊度浮点解 $\hat{N}\in {\mathbb{R}}^n$ 及其协方差 $Q$，优化问题：
$$\underset{N\in {\mathbb{Z}}^n}{\min }\Vert \hat{N}-N{\Vert }_Q^2$$
通过 Cholesky 分解降维后搜索最近整数解，并验证比值 $\frac{\Vert \hat{N}-N_2{\Vert }_Q^2}{\Vert \hat{N}-N_1{\Vert }_Q^2}$ 是否超过阈值。

---

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）

---