[学习]RTKLib详解：pntpos.c与postpos.c

RTKLib详解：pntpos.c与postpos.c

Part A: pntpos.c

一、概述

pntpos.c 是 RTKLIB 中的一个核心模块，用于实现单点定位（Single Point Positioning, SPP）。其主要功能是根据 GNSS 观测数据（伪距和多普勒观测值）和导航星历数据，估计接收机的位置、速度和时钟偏差（PVT）。该模块实现了完整的单点定位功能，涵盖了从卫星轨道计算到接收机位置/速度估计的全过程。其核心是基于伪距观测和多普勒观测的最小二乘估计，并结合电离层、对流层延迟改正及RAIM机制提高解的精度和可靠性，支持多种卫星系统（GPS、GLONASS、Galileo、BeiDou、QZSS 和 SBAS），并实现了包括电离层改正、对流层改正、周跳检测、RAIM（Receiver Autonomous Integrity Monitoring）等功能。

---

二、整体工作流程

pntpos.c 的主函数是 pntpos()，它负责调用其他子函数完成整个定位过程。以下是整体的工作流程：

1. 初始化参数：
   • 设置处理选项（如是否使用广播电离层模型等）。
2. 计算卫星位置、速度、钟差：
   • 调用 satposs() 函数获取所有可见卫星的轨道信息。
3. 估计接收机位置：
   • 调用 estpos() 函数，通过最小二乘法迭代求解接收机坐标和钟差。
4. 残差计算与设计矩阵构建：
   • 在 rescode() 函数中计算伪距残差，并构建观测方程的设计矩阵。
5. 误差建模：
   • 包括电离层延迟 (ionocorr())、对流层延迟 (tropcorr())、伪距测量噪声 (varerr()) 等。
6. 质量控制与验证：
   • 使用 valsol() 进行解的质量检验（GDOP、残差卡方检验）。
7. RAIM（接收机自主完好性监测）：
   • 如果启用 RAIM-FDE（Failure Detection and Exclusion），则调用 raim_fde() 排除异常卫星。
8. 估计接收机速度：
   • 利用多普勒观测值，在 estvel() 中进行速度估计。

---

三、主要函数说明

1. pntpos()

• 功能：主函数，协调整个定位流程。
• 输入参数：
  • obs: 观测数据数组
  • n: 观测数据数量
  • nav: 导航数据结构
  • opt: 处理选项结构体
• 输出参数：
  • sol: 解算结果（位置、速度、时间）
  • azel: 卫星方位角和仰角
  • ssat: 卫星状态信息
  • msg: 错误信息
• 调用关系：
  • satposs()
  • estpos()
  • raim_fde()
  • estvel()

---

2. satposs()

• 功能：计算卫星在地心坐标系中的位置、速度、钟差。
• 数学原理：
  • 根据星历数据（广播星历或精密星历）和当前时间，利用开普勒轨道模型或数值积分方法计算卫星轨道。
• 调用关系：
  • 内部调用 ephpos() 或 peph2pos() 等函数。

---

3. estpos()

• 功能：估计接收机位置和钟差，采用迭代最小二乘法。
• 数学原理：
  • 建立伪距观测方程：
    $$P_i={\rho }_i+c(\delta t_r-\delta t^s)+I_i+T_i+{\epsilon }_i$$
    其中：
    • $ \rho_i $: 几何距离
    • $ c $: 光速
    • $ \delta t_r, \delta t^s $: 接收机和卫星钟差
    • $ I_i, T_i $: 电离层和对流层延迟
    • $ \varepsilon_i $: 测量误差
• 调用关系：
  • rescode()
  • lsq() (最小二乘求解)
  • valsol() (解质量评估)

---

4. rescode()

• 功能：计算伪距残差及其对应的观测矩阵。
• 关键步骤：
  • 计算几何距离、方位角、仰角
  • 电离层改正 (ionocorr())
  • 对流层改正 (tropcorr())
  • 构建观测方程设计矩阵
• 输出：
  • v: 残差向量
  • H: 设计矩阵
  • var: 每个观测的方差

---

5. prange()

• 功能：计算经过码偏差修正后的伪距值。
• 关键步骤：
  • 双频数据：电离层自由组合（IFLC）
  • 单频数据：使用广播电离层模型或TGD改正
• 调用关系：
  • testsnr() (信噪比掩码检查)
  • gettgd() (获取TGD参数)

---

6. ionocorr()

• 功能：计算电离层延迟（L1频率上的值）。
• 支持模型：
  • 广播模型 (IONOOPT_BRDC)
  • SBAS模型 (IONOOPT_SBAS)
  • IONEX TEC模型 (IONOOPT_TEC)
• 输出：
  • ion: L1上的电离层延迟
  • var: 方差

---

7. tropcorr()

• 功能：计算对流层延迟。
• 支持模型：
  • Saastamoinen模型 (TROPOPT_SAAS)
  • SBAS模型 (TROPOPT_SBAS)
• 输出：
  • trp: 对流层延迟
  • var: 方差

---

8. valsol()

• 功能：解的质量评估。
• 评估内容：
  • GDOP（几何精度因子）是否过大
  • 残差平方和是否符合卡方分布
• 输出：
  • msg: 错误信息

---

9. raim_fde()

• 功能：RAIM（接收机自主完好性监测）的FDE算法，用于排除故障卫星。
• 逻辑：
  • 逐颗剔除一颗卫星，重新解算位置
  • 比较不同剔除情况下的RMS残差
  • 选择RMS最小的情况作为最终解
• 调用关系：
  • estpos() (多次调用以剔除不同卫星)

---

10. estvel()

• 功能：利用多普勒观测值估计接收机速度。
• 数学原理：
  • 多普勒观测方程：
    $$f_D=-\frac{1}{\lambda }({\vec{v}}_s-{\vec{v}}_r)\cdot \vec{e}$$
    其中：
    • $ f_D $: 多普勒频移
    • $ \lambda $: 波长
    • $ \vec{v}_s, \vec{v}_r $: 卫星和接收机速度
    • $ \vec{e} $: 卫星到接收机的视线方向单位向量
• 调用关系：
  • resdop() (计算多普勒残差)

---

四、函数调用关系图（Mermaid）

---

五、数学推导补充

1 伪距观测方程

伪距观测值 $ P_i $ 可表示为：

$$P_i={\rho }_i+c(\delta t_r-\delta t^s)+I_i+T_i+{\epsilon }_i$$

其中：

• $ \rho_i $: 卫星与接收机之间的几何距离
• $ c $: 光速
• $ \delta t_r $: 接收机钟差
• $ \delta t^s $: 卫星钟差
• $ I_i $: 电离层延迟
• $ T_i $: 对流层延迟
• $ \varepsilon_i $: 测量噪声

2 最小二乘解

定义残差向量 $ \mathbf{v} $ 和设计矩阵 $ \mathbf{H} $，最小二乘解为：

$$\hat{\mathbf{x}}=({\mathbf{H}}^T\mathbf{WH}{)}^{-1}{\mathbf{H}}^T\mathbf{Wv}$$

其中 $ \mathbf{W} $ 是权矩阵（由各观测的方差决定）。

---

Part B: postpos.c

一、整体作用

postpos.c 是 RTKLIB 的后处理定位核心模块，实现多模式 GNSS 数据后处理定位功能，支持以下特性：

• 支持单点定位(SPP)、差分定位(DGPS)、实时动态定位(RTK)、精密单点定位(PPP)等多种定位模式
• 支持前向/反向/组合处理模式
• 支持多频多系统观测数据处理
• 集成SBAS/PPP-RTK/SSR等增强系统改正数处理
• 提供多种坐标格式输出

注意
该模块通过模块化设计实现了复杂GNSS数据处理功能，但全局变量的使用可能影响多线程处理能力。

二、工作流程

三、关键函数分析

1. postpos

功能: 主处理函数，协调整个后处理流程
输入参数:

• ts/te: 处理起止时间
• ti/tu: 时间间隔和单位
• popt/sopt/fopt: 各类处理选项
• infile/n: 输入文件列表及数量
• rov/base: 流动站/基准站列表
  输出参数: 处理状态码

2. execses_b

功能: 基准站循环处理入口
输入参数:

• ts/te: 当前时段起止时间
• popt/sopt/fopt: 处理选项
• flag: 控制标志
• infile/index/n: 文件信息
• rov/base: 站点列表
  输出参数: 处理状态码

3. execses_r

功能: 流动站循环处理入口
输入参数:

• ts/te: 当前时段起止时间
• popt/sopt/fopt: 处理选项
• flag: 控制标志
• infile/index/n: 文件信息
• outfile: 输出文件
  输出参数: 处理状态码

4. execses

功能: 核心处理会话执行
关键处理步骤:

1. 初始化观测/导航数据
2. 设置天线参数
3. 读取潮汐参数
4. 定位解算
5. 结果输出

5. procpos

功能: 单次定位解算核心
数学原理:

• 使用RTK滤波算法进行状态估计
• 支持前向/反向处理模式
• 动态模式下实时更新SBAS/LEX/RTCM改正数

6. combres

功能: 组合解算结果
数学原理:

• 卡尔曼平滑器组合正反向解
• 质量检验公式：
  $$dr[i{]}^2\le 16\times ({\text{var}}_f[i]+{\text{var}}_b[i])$$
  其中 $dr[i]$ 为位置差，${\text{var}}_f[i]$ 和 ${\text{var}}_b[i]$ 为前向/反向解的方差。

7. inputobs

功能: 输入观测数据及改正数
关键处理:

• 双站观测数据同步
• SBAS/LEX/RTCM改正数更新
• 支持前向/反向数据处理模式

8. readobsnav

功能: 读取观测和导航数据
处理流程:

1. 解析RINEX文件
2. 排序观测数据
3. 去除重复星历
4. 设置时间窗口

9. antpos

功能: 天线位置计算
坐标转换公式:
$$\Delta \mathbf{r}=\mathbf{R}(\varphi ,\lambda )\cdot \delta$$
其中：

• $\mathbf{R}(\varphi ,\lambda )$ 为坐标旋转矩阵
• $\varphi ,\lambda$ 为地理经纬度
• $\delta$ 为局部ENU偏移

10. setpcv

功能: 设置天线相位中心参数
处理流程:

1. 加载卫星天线参数
2. 计算接收机天线相位中心修正
3. 更新处理选项中的天线参数

四、函数调用关系

五、数学原理补充

1. 坐标系转换

• ECEF到LLH转换迭代算法：
  $$\begin{array}{rl}a& =\text{semi-major axis}\\ e^2& =1-{\left(\frac{b}{a}\right)}^2\\ V& =\frac{a}{\sqrt{1-e^2{\sin }^2\varphi }}\end{array}$$
• ENU坐标计算：
  $$\mathbf{R}=\left[\begin{array}{ccc}-\sin \varphi \cos \lambda & -\sin \lambda & \cos \varphi \cos \lambda \\ -\sin \varphi \sin \lambda & \cos \lambda & \cos \varphi \sin \lambda \\ \cos \varphi & 0& \sin \varphi \end{array}\right]$$

2. 改正数组合

组合解算使用加权最小二乘：
$$\begin{array}{rl}\hat{\mathbf{x}}& =({\mathbf{P}}_f^{-1}+{\mathbf{P}}_b^{-1}{)}^{-1}({\mathbf{P}}_f^{-1}{\mathbf{x}}_f+{\mathbf{P}}_b^{-1}{\mathbf{x}}_b)\\ \hat{\mathbf{P}}& =({\mathbf{P}}_f^{-1}+{\mathbf{P}}_b^{-1}{)}^{-1}\end{array}$$

3. 时间同步

采用GPST时间系统，时间转换关系：
$$\begin{array}{rl}\text{GPST}& =\text{UTC}+\text{leap\_seconds}\\ \text{UTC}& =\text{GPST}-\text{leap\_seconds}\end{array}$$

4. 野值剔除

通过QF值质量检验：
$$\text{QF}=\frac{{\text{残差}}^2}{\text{方差阈值}}$$

六、内存管理

• 观测数据使用动态数组管理
• 解算结果采用双缓冲技术
• 注意全局变量(obss, navs等)的生命周期控制

七、特殊处理机制

1. 野值剔除: 通过QF值质量检验
2. 多路径抑制: 使用历元间差分检测
3. 周跳修复: 在rtkpos中实现LAMBDA方法
4. 电离层延迟: 双频组合消除一阶项

---

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

---