GREAT-PVT非差非组合PPP源码阅读

GREAT-PVT 是武汉大学测绘学院 GREAT 团队于 2024 年开源发布的多频多系统 GNSS 精密定位软件。相较于之前的开源 PPP 处理软件，在多系统多频非差非组合观测模型、偏差处理和模糊度固定策略等方面实现得更为完善，可能比 RTKLIB 等老程序更适合作为 GNSS 算法入门后进行多频多系统 UDUC PPP-AR 算法学习和二次开发的平台。

• 一、简单试用
• 1、下载编译
• 2、解算测试
• 二、项目框架
• 1、文件结构
• 2、关键数据结构
• 3、关键函数
• 4、执行流程
• 三、PPP 浮点解及误差改正
• 1、UDUC-PPP 浮点解方程
• 2、周跳探测
• 3、精密误差改正
• 4、多系统多频的偏差处理
• 四、PPP-AR 算法
• 1、UPD 模型
• 2、EWL/WL/NL 级联
• 3、模糊度固定
• 4、质量控制
• 五、总结
• 1、多频算法要点
• 2、容易误导的名称
• 3、未实现算法
• 4、二次开发经验

一、简单试用

官方介绍：为满足高精度卫星导航与定位需求，武汉大学测绘学院 GREAT 团队设计并开发了一款面向全球导航卫星系统（GNSS）用户的开源软件 GREAT-PVT。GREAT-PVT 是 GREAT 软件中的一个重要模块，主要用于精密定位解算。软件采用 C++ 编写，支持用户在 Windows 和 Linux 操作系统上编译安装。

GREAT-PVT 支持多频多系统精密单点定位（PPP）及模糊度固定（PPP-AR），可采用无电离层组合或非差非组合模型。同时，该软件还支持混合频率模式的实时动态定位（RTK）解算。官方论文利用 124 个 MGEX 测站 2023 年第 305～314 天的观测数据开展 PPP-AR 精度评估，并通过动态实验测试 RTK 性能。软件包提供样例数据、说明文档、批处理脚本、绘图脚本及辅助工具，支持多天连续数据处理与分析。项目地址为 GREAT-WHU/GREAT-PVT。

1、下载编译

从 GitHub 下载的项目文件结构如下：

GREAT-PVT/├── src/│   ├── LibGnut/          # 基础库：滤波、模型、数据 I/O、工具│   │   ├── gproc/        # SPP/PPP 基础处理与滤波器│   │   ├── gmodels/      # 误差模型：潮汐、对流层、电离层、相对论、相位缠绕│   │   ├── gset/         # XML 配置解析（gsetproc、gsetgnss 等）│   │   ├── gutils/       # 时间、观测类型、卫星数据、字符串处理│   │   ├── gio/          # 通用 I/O、日志│   │   ├── gcoders/      # 标准格式编解码器（RINEX、SP3、ANTEX 等）│   │   ├── gdata/        # 数据实体：观测、星历、测站/卫星对象│   │   ├── gall/         # 聚合容器：观测、精密产品、偏差、天线、潮汐负荷│   │   ├── gprod/        # 滤波输出产品：坐标、钟差│   │   └── newmat/pugixml/spdlog/  # 第三方库：矩阵运算、XML解析、日志│   ├── LibGREAT/         # 高级库：PPP-AR、UPD、精密偏差、组合模型│   │   ├── gproc/        # PPP 主处理类（gpvtflt、gpppflt）│   │   ├── gmodels/      # 观测方程构建（gcombmodel、gprecisebias）│   │   ├── gset/         # AR 相关 XML 配置解析（gsetamb）│   │   ├── gambfix/      # 模糊度固定：gambiguity、glambda、gbdeci│   │   ├── gcoders/      # GREAT 专用产品编解码器（UPD、IFCB、EOP）│   │   ├── gdata/        # GREAT 专用数据类（UPD、IFCB、极移等）│   │   └── gutils/       # GREAT 专用工具（坐标框架转换）│   ├── app/│   │   └── GREAT_PVT/    # 主程序入口│   └── CMakeLists.txt    # CMake 构建脚本├── doc/                  # 文档与配置模板├── sample_data/          # PPP/RTK 示例数据├── util/                 # 批处理脚本、绘图脚本、极移/UT1 等文件└── bin/                  # 预编译二进制 (Linux/Windows)

doc 文件夹包含中英文使用手册，以及 PPP 和 RTK 的示例 XML 配置文件；src 文件夹是项目源码，其中 GREAT_PVT.cpp 是主程序入口；util 文件夹包含批处理和绘图脚本，使用前需要先解压；sample_data 文件夹提供 PPP 和 RTK 示例数据，使用前同样需要解压。

我的环境是 WSL Ubuntu 20.04，使用 VS Code 调试开发。打开 GREAT-PVT 文件夹后，选择 src/CMakeLists.txt 作为项目的 CMake 入口即可编译；生成的可执行程序位于 ./build_Linux/Bin/GREAT_PVT。不同编译器和依赖版本下的结果可能不同。

2、解算测试

GREAT-PVT 项目包含 PPP 测试数据集 PPPFLT_2023305，其中提供 2023 年第 305 天（11 月 1 日）两个 IGS 测站（GODN、HARB）的观测数据、精密产品、各类 PPP 模式的 XML 配置及示例结果：

PPPFLT_2023305/├── obs/│   └── *.23o                          # RINEX 3.04 观测文件（伪距+载波相位）├── gnss/│   ├── COD0MGXFIN_*.SP3               # CODE 精密星历（SP3-d，5 min 间隔）│   ├── COD0MGXFIN_*.CLK               # CODE 精密钟差（RINEX CLK，30 s 间隔）│   ├── CAS0MGXRAP_*.BSX               # CAS 差分码偏差（DCB，BIAS-SINEX 格式）│   ├── brdc*.23p                      # 混合多系统广播星历（RINEX NAV 3.04）│   └── ifcb_*                         # GPS L5 频间钟偏差（IFCB，GREAT 自定义格式）├── model/│   ├── igs20_*.atx                    # IGS20 天线相位中心改正（ATX/ANTEX）│   ├── jpleph_de405_great             # JPL DE405 行星历表（GREAT 二进制）│   ├── oceanload                      # 海洋负荷潮汐改正（BLQ 格式）│   └── poleut1                        # 极移 / UT1-UTC 参数（ERP 格式）├── upd/│   ├── upd_wl_*_{G,E,C}               # 宽巷（WL）UPD 产品，按 GPS/GAL/BDS 分文件│   ├── upd_nl_*_{G,E,C}               # 窄巷（NL）UPD 产品│   ├── upd_ewl_*_{C,E}                # 超宽巷（EWL）UPD 产品│   ├── upd_ewl24_*_{C,E}              # EWL(2,4) UPD 产品│   └── upd_ewl25_*_{C,E}              # EWL(2,5) UPD 产品├── xml/│   ├── GREAT_PPPFLT_static_*.xml      # 静态模式配置（3 种：DF_Fixed / DF_Float / FF_Fixed）│   └── GREAT_PPPFLT_kinematic_*.xml   # 动态模式配置（3 种：DF_Fixed / DF_Float / FF_Fixed）└── result/    └── *-PPP_kin_DF_Fixed.flt         # 各种 PPP 配置的定位结果（坐标、速度、精度、Ratio 等）

运行程序的命令行格式：GREAT_PVT -x <config.xml>，需要注意以下问题：

1. XML 配置中的路径均以 PPPFLT_2023305 为基准：
• 如果是直接在命令行终端运行，需要切换工作目录到 PPPFLT_2023305；
• 如果用 VSCode 的 launch.json 配置调试运行需要把 cwd 项设置为 PPPFLT_2023305 路径；
2. 结果默认输出到 PPPFLT_2023305/result，可预先重命名示例结果，避免被新结果覆盖。

二、项目框架

1、文件结构

GREAT-PVT 源码包括 app、LibGnut、LibGREAT 三部分，app 里的 GREAT_PVT.cpp 是主函数，调用 LibGnut、LibGREAT 的代码来实现 PPP-AR：

LibGnut 是整个 GREAT-PVT 软件的底层基础设施库，GREAT 团队开源的其它 GNSS 程序（如 GREAT-UPD）也包含这套基础库。它封装了 GNSS 定位所需的基础功能，包括标准数据格式读写（RINEX、SP3、ANTEX 等）、卫星轨道与钟差管理、各类物理模型（对流层、电离层、潮汐、相对论等），以及 Kalman 滤波和平方根滤波通用框架，可被 SPP、PPP、RTK 等定位模式复用。此外，LibGnut 还内嵌 Newmat（稠密矩阵运算）、pugixml（XML 解析）和 spdlog（日志）三个第三方库。

LibGnut 源自 G-Nat 程序包，文件结构和做数据质量分析的Anubis基本上是一样的，XML配置文件也相似，老师几年前的博客做过 G-Nut/Anubis 的介绍：

迷途兄，公众号：不迷途导航程序员G-Nut/Anubis 源码初读小记

LibGREAT 是 GREAT 团队在 LibGnut 之上构建的PPP-AR 专用算法库，实现了从 PPP 浮点解到当历元固定解的精密定位流程。它构建了 IF（消电离层组合）和 UDUC（非差非组合）两种观测方程模型，并实现 LAMBDA 整数搜索、部分模糊度固定和 UPD 改正等 PPP-AR 关键技术。LibGREAT 中的处理类（如 t_gpvtflt）直接继承或调用 LibGnut 的滤波器、数据容器和误差模型等通用能力。

2、关键数据结构

GREAT-PVT 的数据模型围绕 t_gdata 基类展开，向上聚集为全局容器（gall），向下派生为具体实体（gdata）。

类继承层次：

t_gdata                         ← 所有数据的抽象基类├── t_gobsgnss                  ← 单颗卫星的观测值集合（apply_bias()、MW()、LNL() 等）│   └── t_gsatdata              ← 卫星数据容器（附加轨道、钟差、高度角、PCO/PCV 对象等）├── t_gobj                      ← 天地对象基类（测站 / 卫星的公共属性）│   └── t_grec                  ← 接收机（测站）对象├── t_gallobs                   ← 全局观测容器：管理所有测站、所有历元的原始观测├── t_gallnav                   ← 全局星历容器：管理所有系统的广播/精密星历│   └── t_gallprec              ← 精密产品容器（继承 t_gallnav，附加钟差、ERP 等）├── t_gallbias                  ← 全局偏差容器：管理 BIAS-SINEX/Bernese BIA 产品（DCB/OSB）├── t_gallobj                   ← 全局对象容器：管理所有测站和卫星对象（含天线信息）├── t_gallpcv                   ← 全局天线容器：管理 ANTEX 文件（PCO/PCV 格网）├── t_gallpar                   ← 状态参数向量（独立于 t_gdata，不继承它）├── t_gallproc                  ← 全局处理配置容器└── t_gallprod                  ← 全局产品容器（滤波输出的钟差/坐标等）

关键容器的持有关系：

t_gallobs                      t_gallnav / t_gallprec│  .obs(site, epoch)           │  .eph(sat, epoch)│  返回 vector<t_gsatdata>     │  返回 t_geph / t_gephprec↓                              ↓t_gsatdata ──包含──→ t_gobsgnss ──包含──→ 各频段 t_gobs（码/相位原始值 + LLI 等属性）│                               │├── satcrd() 卫星位置           ├── apply_bias()  施加 DCB/OSB 改正├── clk()    卫星钟差           ├── MW()          Melbourne-Wübbena 组合├── ele()    高度角             ├── LNL()         (2,-1) 窄巷组合├── rho()    几何距离           └── wavelength()  频率对应波长└── pcv()    天线 PCV 对象

滤波器继承层次：

t_gflt                          ← 滤波器抽象基类（update()/add_data()/add_virtual_obs()）├── t_kalman                    ← 标准 Kalman 滤波（Joseph 协方差更新）├── t_SRF                       ← 平方根协方差滤波└── t_SRIF                      ← 平方根信息滤波t_gspp                          ← SPP 基类（虚基类）├── t_gsppflt                   ← SPP 滤波器（_syncSys/_syncIono/_syncIFB/_syncAmb 参数管理）└── t_gppp                      ← PPP 基类（虚继承 t_gspp）    └── t_gpppflt               ← PPP 滤波基础类（多重继承 t_gppp + t_gsppflt + t_gxml）        └── t_gpvtflt           ← PPP-AR/RTK 主处理类（processBatch/ProcessOneEpoch/_amb_resolution）

模型类层次：

t_gbasemodel                    ← 观测方程模型基类├── t_gcombmodel               ← 组合模型（持有 _bias_model）│   ├── t_gcombIF              ← IF 消电离层组合（cmb_equ_IF()）│   ├── t_gcombALL             ← UDUC 非差非组合（cmb_equ()）│   └── t_gcombDD              ← 双差组合（继承 IF + ALL）│t_gbiasmodel                   ← 偏差/设计矩阵模型└── t_gprecisebias             ← 精密偏差计算基类    └── t_gprecisebiasGPP      ← PPP 精密偏差（_prepare_obs_GPP/cmpObs/PCV）

一个历元内的数据流动：

t_gallobs                       t_gallnav                    t_gallbias│ 当前历元观测                   │ 星历+钟差                  │ DCB 产品↓                               ↓                            ↓vector<t_gsatdata> ──预处理──→ _data ──建方程──→ A/P/l ──滤波──→ _param/_Qx                                    │                            │                                    └── AR ──→ _param_fixed ────→ 输出

3、关键函数

4、PPP 及 PPP-AR 执行流程

GREAT-PVT 的 PPP 和 PPP-AR 处理流程以 t_gpvtflt 类为核心，采用逐历元推进的处理模式：程序从 XML 配置中读取起止时间与采样间隔，依次处理每个历元的观测数据，上一历元滤波得到的状态向量作为下一历元的先验信息持续传递。整个流程可划分为四个层次：

1. 首先是观测数据读取和 DCB 改正；
2. 其次是测站预处理与 Kalman 滤波状态预测；
3. 然后是观测方程构建与滤波更新（得到 PPP 浮点解）；
4. 最后是在当历元浮点解基础上进行模糊度固定和约束更新（得到 PPP 固定解）。

下面的调用链从 main() 出发，逐层展开到每个关键环节的具体函数：

main()  └ t_gpvtflt::processBatch()       └ while (time_loop)                      // now 从 _beg_time 步进，到越过 _end_time 结束           ├ _slip_detect()                         [gpvtflt.cpp / gpreproc.cpp]           │   └── t_gpreproc::ProcessBatch()           │       └── MW/EWL/(2,-1)级联周跳探测 + 钟跳检测修复           └ ProcessOneEpoch()               ├ _getData()                         [gpvtflt.cpp]               │   └── 从 t_gallobs 中读取当前历元观测               ├ t_gobsgnss::apply_bias(_gallbias) [gobsgnss.cpp]  // t_gsatdata 继承调用               │   ├── AC == SGG_A: apply_code_phase_bias() // OSB 框架，尚未真正实现               │   └── 其它 AC:      apply_dcb()             // DCB               ├ _timeUpdate()               └ _processEpoch()                   ├ _prepareData()                  [gpvtflt.cpp]                   │   ├── _preprocess()             // 卫星位置、Bancroft 初值、高度角等                   │   └── _gquality_control.processOneEpoch()                   │       ├── excludeBadObs()       // SNR/缺频/码差等检查                   │       └── apply_IGSO_MEO()      // BDS2 SCB                   ├ _predict()                      // 状态同步与时间更新                   │   ├── _syncSys() / _syncIono() / _syncIFB()                   │   ├── _syncAmb()                // 在建方程前增删本历元模糊度参数                   │   ├── _predictCrd()             // 动态时以本历元 SPP 坐标重置先验                   │   ├── _predictClk()             // 钟差白噪声                   │   ├── _predictBias()            // ISB/IFB 随机游走                   │   ├── _predictIono() / _predictTropo()                   │   └── _predictAmb()             // rndwk_amb 可配置，示例配置为 0                   ├ _cmp_equ()                   │   └── for sat in _data: _base_model->cmb_equ()                   │       ├── 组合模型再次调用 apply_bias()，但 _dcb_correct_mark 使其直接返回                   │       ├── _bias_model->cmb_equ() [gprecisebiasGPP.cpp]                   │       │   ├── 固体潮/海潮负荷/极潮/海洋极潮负荷                   │       │   ├── PCO/PCV、对流层、电离层、相位缠绕、相对论                   │       │   └── ISB、IFB 等参数项                   │       ├── IF 模式:   t_gcombIF::cmb_equ_IF()                   │       └── UDUC 模式: t_gcombALL::cmb_equ()                   │           ├── 各频点分别建立 P/L 方程                   │           └── GPS 第 3 频相位按配置施加 IFCB 产品改正                   ├ _filter->update(A, P, l, dx, _Qx) [gflt.cpp]                   │   ├── t_kalman::update()         // Kalman/Joseph 协方差更新                   │   └── 或 t_SRF::update()         // 平方根协方差滤波                   ├ _filter->add_data(...)           // 缓存浮点解和方程，供本历元 AR 使用                   ├ _amb_resolution()                [gpvtflt.cpp]                   │   ├── 设置观测类型、活动模糊度、MW/高度角/SNR                   │   ├── _getSatRef()               // 按弧长×高度角并结合 UPD 可用性选参考星                   │   ├── processBatch("EWL/EWL24/EWL25")                   │   │   └── UPD 改正 → 直接取整筛选 → 独立子集 → EWL 虚拟约束                   │   ├── processBatch("WL")                   │   │   └── UPD 改正 → 直接取整筛选 → 独立子集 → WL 虚拟约束                   │   └── processBatch("NL")                   │       ├── _defineDDAmb() / _applyUpd() / _fixAmbUDUC() 或 _fixAmbIF()                   │       ├── _selectAmb()                   │       ├── _ambSolve() → _lambdaSearch() → LAMBDA4()                   │       │   └── LDLᵀ 分解、Z 变换、整数搜索、Ratio/Bootstrap、部分固定                   │       └── _addFixConstraint() → add_virtual_obs() → update()                   ├ if (_amb_state) _postRes()       // 重组固定解相位残差，不负责判废固定解                   └ for each param: value += dx      // 主状态仍用浮点 dx 更新

单历元内的数据流可以概括为：

上一历元 _param/_Qx    → _predict() 时间更新    → 主滤波 update() 得到浮点 dx/_Qx    → add_data() 复制到 _filter    → _amb_resolution() 在 _filter 上施加整数约束    → _param_fixed 写出当历元 Fixed/Float 结果    → 主 _param 仍执行 value += dx    → 下一历元继续使用浮点 _param/_Qx

下面是几个关键的入口函数，调试程序看代码的时候可以先在这些地方设断点：

• t_gpvtflt::processBatch() 是整个测站/整个弧段的 PPP 批处理入口，它负责初始化时间窗口、按采样率遍历所有历元，并在循环中调用 ProcessOneEpoch()。
• ProcessOneEpoch() 是单历元 PPP 的入口函数，每个历元只执行一次：它先从 t_gallobs 中读取当前历元观测，再调用 t_gobsgnss::apply_bias()（由 t_gsatdata 继承调用），之后调用 _processEpoch()。
• _processEpoch() 是单历元 PPP 的核心处理函数，完成预处理、状态预测、观测方程构建、Kalman 滤波更新等浮点解全过程。
• _amb_resolution() 是单历元 PPP-AR 的入口，它在浮点解之后被调用，依次执行 EWL/WL/NL 级联模糊度固定，并生成当前历元的固定解输出。

三、PPP 浮点解及误差改正

本章详细介绍 UDUC-PPP 浮点解的观测方程、待估参数、设计矩阵、定位模式，以及 PPP 中涉及的周跳探测、精密误差改正和多系统多频偏差处理。这些是理解 PPP 数据预处理与滤波估计的基础。

1、UDUC-PPP 浮点解方程

非差非组合（UDUC）模式直接以每颗卫星、每个频点的原始伪距和相位观测值作为方程输入，不再通过线性组合消除电离层等参数。其优势在于保留了全部原始观测信息，估计得到的大气延迟可用于 PPP-RTK 增强建模或其它科学研究。本小节从数学方程、待估参数和设计矩阵几个方面介绍 UDUC-PPP 观测模型。

1、观测方程

设测站  观测卫星  在频点  上的原始观测值（经 PCO/PCV、潮汐、相对论、相位缠绕等精密改正后），线性化后的观测方程为：

伪距观测：

载波相位观测：

其中：

在精密定位中，观测值加权是平衡不同类型、不同质量数据贡献的手段。伪距观测噪声通常为分米至米级，而载波相位噪声仅为毫米至厘米级，两者的精度差异可达两个数量级，因此必须通过合理的权比设置避免伪距噪声"污染"载波的高精度信息。同时，实际观测环境中存在信号遮挡、多路径、低高度角等不利因素，导致部分观测值质量显著下降，此时需要依据观测质量动态降权或剔除，以确保参数估计的最优性和可靠性。

GREAT-PVT 在 XML 配置中为各 GNSS 系统设定伪距和相位的先验精度，示例配置常用 sigma_C="0.6" sigma_L="0.01"。sigma_C 直接作为米单位的码标准差使用；sigma_L 在 gbiasmodel.cpp 构权时还会乘当前频点波长，因此其数值应按周理解。例如，代码中的 GPS 默认值为伪距 2.0 m、相位 0.02 周，而不是 0.02 m。未配置 obs_weight 时默认采用 SINEL2，示例 XML 则显式选择 SINEL。

2、待估参数

上表顺序对应源码中 _param 的实际构建顺序：t_gsppflt 构造函数先加入 CRD_X/Y/Z、CLK、TRP；随后每历元由 _syncSys()、_syncIono()、_syncIFB()、_syncAmb() 依次动态加入 ISB、SION、IFB、AMB。

• VION 和 ZWD梯度估计没有真正实现；
• 很多的随机模型配置解析都没有在算法中实现，比方说配置解析支持设置 SION 随机游走噪声，但当前代码每个历元都会把 SION 的值和方差重置，后面的随机游走加噪条件因而不会成立。
• 从参数类型和偏导定义看，GLO_IFB、GLO_IFCB、GLO_IFPB 按 GLONASS 卫星设计；但当前 gsppflt.cpp::_syncIFB() 只为 GPS/GAL/BDS/QZS 创建接收机级 IFB，没有创建上述 GLONASS 参数，因此它们尚未进入实际估计流程。

3、设计矩阵

设计矩阵  的每一行对应一个实际建立的伪距或相位方程。设第  颗卫星有  组同时具备码和相位的可用频率，配置的频率上限为 ，则：

• 观测数：；只有所有卫星各频数据完整时才等于 
• 参数数：（省略未启用项）
• 设计矩阵维度：

伪距观测行（系统 、卫星 、频点 ）：

载波相位观测行（系统 、卫星 、频点 ）：

各分块说明：

数量示例：假设GEC系统双频、每系统平均 10 颗可见卫星、双频（L1/L2）、开启 SION。

• 基础参数：3（坐标）+ 1（钟差）+ 1（ZWD）+ 2（ISB）= 7
• 电离层：约 30 颗卫星 = 30
• 模糊度：约 30 颗 × 2 频 = 60
• 单历元总参数约 100 个，状态向量随卫星升降动态变化。

4、定位模式

在 GREAT-PVT 中，坐标估计策略由 <process> 节点下的 crd_constr 和 pos_kin 共同控制，crd_constr 决定坐标是否进入主滤波器估计，pos_kin 影响坐标预测阶段的值更新方式，主要的配置有三种：

钟差在任何模式下都会重置，_predictClk() 不区分 crd_constr 和 pos_kin，每历元都把 CLK 值设为 SPP 估计值。

速度值由多普勒观测计算，受 doppler 配置项控制。PPP 主滤波状态向量中不存在 VEL_X/Y/Z 和 CLK_RAT；它们仅在 doppler=true 时作为多普勒测速子滤波的临时局部参数独立估计，结果只用于文件输出，不参与坐标或钟差的历元间预测。

2、周跳探测

在 GNSS 载波相位观测中，接收机跟踪环路可能因信号遮挡、严重多路径、低信噪比或动态应力等原因导致相位整周计数跳变。精密定位中模糊度被当作时不变常数估计，一旦发生未处理的周跳，滤波器会将跳变后的模糊度当作同一常数继续估计，导致该卫星在所有后续历元中产生累积性偏差。周跳探测的核心思路是借助多频观测值组合消除几何距离、钟差等公共项，将周跳从相位变化中隔离出来：

GREAT-PVT 的周跳探测整体上是两趟扫描的结构，外层入口为 t_gpvtflt::_slip_detect()：

• 第一趟（逐星）：对每颗卫星执行 MW/EWL/NL 级联探测 → _transform() 反解原始频点周跳 → 写 LLI 标记；
• 第二趟（全局）：对所有 GPS 卫星执行钟跳检测 → 若确认钟跳则擦除误标周跳并整体修复相位；

周跳标记传入后续滤波流程后，由 _syncAmb() 依据 LLI 重置对应模糊度，_predictAmb() 仅在此基础上为已有模糊度协方差添加过程噪声。下面逐环节展开。

_slip_detect() 调用 t_gpreproc::ProcessBatch()，内部对每颗可见卫星执行以下步骤：

1. 读取与排序。 获取当前历元与上一历元同颗卫星的观测，用 sort_band() 按频率从高到低排列可用频段（sorted_t1）。

2. 级联探测（_slip()）。 按频段数  循环  级，每级探测相邻两频段间的 WL/EWL 周跳，最后一级额外执行 NL 反推：

• 第 1 级（MW 分支，i=1）：取排序序列尾部最后两个频段构造 MW 组合，历元间差分除以 。若  则记录 WL 周跳到 _v_lcslp[1]。双频系统只有这一级。
• 中间级（EWL 分支，i=2,\dots,d-1）：计算相邻两级宽巷历元间差分之差（dEWL − dWL），扣除前级已探测的周跳贡献后，若  则记录到 _v_lcslp[i]。
• 最后级附带的 NL 分支：当 （即循环最后一级）时，额外用 (2,−1) 组合 LNL() 计算窄巷历元差分，结合 _iono() 估计的 GF 电离层变化率 _dI 和已探测的 WL 周跳，反推窄巷周跳量，记录到 _v_lcslp[i+10]。

3. 反解原始频点周跳（_transform()）。 若某卫星 _v_lcslp 非空，联立各级 WL/EWL/NL 周跳建立线性方程组 ，求解得到各原始频点的周跳值 。只要方程组存在非零解，就对参与求解的全部频点统一写 LLI=1（包含周跳为 0 的频点）。

4. 钟跳检测（_jumps()）。 对所有 GPS 卫星，比较相位平均跳变量  与伪距平均跳变量 。若所有可见 GPS 卫星的  一致且为整数毫秒，则判定发生了接收机钟跳。之所以只用 GPS，是因为 GPS 卫星数量通常足够多， 的一致性可作为可靠的判别依据。

5. 钟跳修复。 若确认钟跳，先调用 _remove_slip() 擦除本历元已记录的全部周跳标记（钟跳导致的相位跳变已被误判为周跳），再调用 _repair() 对所有卫星的所有载波相位观测整体平移  米（ 为钟跳毫秒数），避免大面积模糊度重新初始化。

6. 传递标记。 LLI=1 的观测传入后续滤波流程，由 _syncAmb() 检测到后重置对应模糊度参数并赋予较大初始方差。

值得介绍一下细节：

1. MW 与 EWL 组合的选择逻辑

双频仅有 L1/L2，无法直接构造既波长足够长、又能同时消去几何与钟差的相位组合，因此需借助 MW：用 L1−L2 宽巷相位结合 P1−P2 伪距组合消去几何和钟差，以探测周跳；而三频可直接用 L2/L5（或 B1I/B3I 等）构造米级超宽巷相位组合EWL，其波长已足够长，对几何与钟差不敏感，无需引入噪声较大的伪距即可直接探测周跳，因而更可靠。

多频情况下，代码并非直接构建一个独立的 EWL 组合，而是计算相邻两级宽巷历元间差分之差（dEWL − dWL），再利用前一级已检测的 WL/EWL 周跳进行改正，从而隔离出当前级的额外周跳。但这种隔离的优势在当前版本的代码中没有利用上，代码中一旦检查到周跳，最后还是对所有频点都重置模糊度。

2. (2,−1) 窄巷组合与电离层改正

并非标准窄巷组合，而是周跳探测中使用的 (2,−1) 相位组合。其定义为：

其中 、 为以米为单位的相位观测值。该组合保留几何距离、钟差和对流层等公共项，等效波长为：

对 GPS L1/L2 约为 15.6 cm。在宽巷周跳  已知后，可配合 (2,−1) 组合求解 L1/L2 上的具体周跳：

其中  即为 (2,−1) 组合的整周模糊度。采用整数系数 (2,−1) 而非标准窄巷 (f_1 N_1+f_2 N_2)/(f_1+f_2)，好处是便于与宽巷建立简单的整数线性关系，从而把 WL/EWL 无法分离的 L1、L2 单独周跳反推出来。

在窄巷 (2,−1) 反推中，_iono() 负责估计并扣除相邻历元间的电离层延迟变化。(2,−1) 组合虽然保留了整周模糊度的整数特性，但它并非电离层无关组合，其历元间差分中仍包含电离层时变项；若不加改正，电离层闪烁或缓慢变化很容易被误判为周跳。_iono() 利用无几何（GF）组合 L1 - L2 的历元间差分进行估计：

k   = (f1 / f2)^2dL  = (L1_t2 - L2_t2) - (L1_t1 - L2_t1)   // GF 组合历元间差分，单位 m_dI = -dL / (k - 1)                          // 电离层变化率，单位 m

得到 _dI 后，NL 反推公式中的 disf 把 _dI 缩放为 (2,−1) 组合中的电离层贡献，再从 dNL 中扣除。剩余差分主要反映真实的 (2,−1) 组合周跳。

3. 组合周跳到原始频点的反解

_transform() 以三频（L1/L2/L5）为例，sorted_t1 = [L1, L2, L5]，d = 3。设 O₁, O₂, O₃ 为待求的原始频率周跳，wlSlp₁ 为 L2-L5 的 WL 周跳，wlSlp₂ 为 L1-L5 的 EWL 周跳，nlSlp 为 L1-L2 的 (2,−1) 组合周跳，则：

求解得到 O₁, O₂, O₃ 后，_transform() 将它们作为各原始频点的周跳记录写入 _gslip/_mslips。一旦 trans=true，代码会给 m_orig 中所有参与求解的频点写入 LLI=1，即便某个频点求得的 O(i) 为 0，后续滤波器依据 LLI 重置相应模糊度。

恢复原始频点模糊度的算法可以继续深入改进：

1. 只重置周跳频点的模糊度：多频组合探测有没有周跳，_transform() 转换周跳到原始频点来判断是哪个频点周跳；
2. 尝试修复周跳：对于有周跳的频点不做模糊度重置，而是给模糊度加上周跳值，并加一个小点的过程噪声；

本节只介绍源码实现；这套周跳探测与钟跳修复算法在电离层闪烁、低成本接收机和城市复杂环境中的实际表现，需通过实测数据进一步评估。

3、UDUC 精密误差改正

精密单点定位之所以能够达到厘米级甚至毫米级的定位精度，关键在于对各类系统误差和物理效应进行完备而精确的改正。GNSS 信号从卫星发射到被地面接收机接收，其间会受到卫星端天线相位中心偏差和变化、接收机端天线相位中心偏差和变化、电离层和对流层延迟、相位缠绕、相对论效应、固体潮、海潮负荷、极潮以及地球定向参数变化等多种因素的影响。

GREAT-PVT 在 gprecisebias.cpp 和 gprecisebiasGPP.cpp 中分不同的步骤实现了上述精密改正：

这里着重介绍一下天线PCO/PCV改正。卫星和接收机天线并非理想点源，信号实际相位中心会随频率、方位角和高度角变化而与几何中心偏离，这种偏差由相位中心偏移（PCO）和相位中心变化（PCV）两部分组成：

• PCO 描述的是天线相位中心相对于天线参考点（ARP）的固定三维偏移，它影响几何距离，因此码和相位都需要改正。
• PCV 描述的是载波相位中心随卫星高度角和方位角的变化量，体现的是”载波相位”等相位面的微小形变，只存在于相位观测中。

IGS 早期采用相对天线校准模型，依赖特定参考天线并通过短基线差分消除公共偏差，无法精确建模绝对相位中心；2011 年全面切换至绝对校准模型后，精度显著提升，但多频点仍共用一套 PCO/PCV 值，或仅对 L1/L2 区分，其他频点通过简单内插或复制获得，导致多频 PPP-AR 中频间偏差（IFB）与天线相位中心混叠，难以分离。

2022 年 11 月 27 日 IGS20 参考框架正式启用后，配套的 igs20.atx 文件实现了逐频点独立校准，为各 GNSS 系统的各信号频点均提供独立的 PCO 和 PCV 值。采用 igs20.atx 后，用户需按频点索引独立读取对应 PCO/PCV 值，替代早期统一或内插的处理方式，确保多频观测值的天线改正相互独立，有效削弱天线频间偏差与 IFB 的混叠。逐频点独立的天线改正能够显著提升 UPD/IFCB 估计的纯净度，改善多频 PPP-AR 的模糊度固定性能。

逐频点天线改正的具体实施包含三个步骤：

1. 首先，依据接收机天线型号和卫星PRN号从igs20.atx文件中匹配对应的天线类型，提取各频点的PCO三维偏移量和PCV改正表；
2. 其次，根据卫星高度角和方位角插值获取当前观测方向的PCV改正值，对卫星端PCO需结合卫星姿态模型计算天线固定坐标系到轨道坐标系的旋转矩阵，将PCO矢量投影至视线方向；
3. 最后，将PCO与PCV改正值叠加至原始观测值，其中PCO通过几何投影转换为距离改正，PCV直接作为附加相位偏差，两者均以米为单位与伪距和载波相位统一量纲后参与后续定位解算。

卫星天线PCO是在卫星本体坐标系（如Z轴指向地心、Y轴沿太阳帆板方向）下定义的，而PPP观测方程需要将其转换到地心地固坐标系（ECEF）或视线方向，这需要卫星姿态，GREAT-PVT 不支持姿态文件读取，只能采用名义姿态模型进行近似处理，即根据太阳-卫星-地球几何关系实时计算理论偏航角，假设卫星完美对日定向。这种方法对GPS Block IIR/IIF、Galileo等标准动态偏航卫星足够近似，但在正午/午夜偏航机动期间会产生显著误差（可达数厘米PCO投影误差），且完全无法处理BDS GEO/IGSO等特殊轨道卫星的姿态异常。如果服务端产品和用户端产品都不用姿态产品，一部分姿态误差可以吸收到服务端估计的偏差参数中；所以说如果服务端和用户端算法一致，影响可以稍微减弱。

GREAT-PVT 中天线 PCO/PCV 数据由 t_gallpcv（管理 ANTEX 文件，位于 gallpcv.cpp）统一读取，并通过 t_gpcv（位于 gpcv.cpp）提供 pcoS()/pcoR()、pcoS_raw()/pcoR_raw() 和 pcvS()/pcvR() 等接口。索引是 GOBSBAND，随后转换到 GFRQ：

// gpcv.cpp / gprecisebiasGPP.cpp 中的核心调用// 公共几何 PCO：IF 模式组合 band1/band2；RAW 模式只取 band1t_gtriple pco(0, 0, 0);if (sat_pcv->pcoS(satdata, pco, lc, band1, band2) > 0)    _crs_sat_pco += t_gtriple(_rot_matrix * pco.crd_cvect());// 接收机端同理t_gtriple pco(0.0, 0.0, 0.0);if (rec_pcv->pcoR(satdata, pco, lc, band1, band2) > 0)    _crs_rec_pco += t_gtriple(_rot_matrix * pco.crd_cvect());

当前代码有两层 PCO 处理。首先，_prepare_obs_GPP() 为公共几何计算 PCO：IF 模式使用双频组合 PCO，RAW_ALL/RAW_MIX 模式的 pcoS()/pcoR() 实际只取传入的 band1（第一频段）。卫星端 PCO 从卫星本体坐标系旋转到 CRS；接收机端 PCO 先由 NEU 转到 TRS，再在 isCRS=true 时转到 CRS：

// gprecisebiasGPP::cmb_equ() 内部// 1. 获取卫星和接收机的 PCV 对象shared_ptr<t_gpcv> sat_pcv = sat_obj->pcv(crt_epo);shared_ptr<t_gpcv> rec_pcv = rec_obj->pcv(crt_epo);// 2. PCO 改正（最终统一到 CRS）if (sat_pcv) {t_gtriple pco(0, 0, 0);    sat_pcv->pcoS(_crt_obs, pco, _observ, band1, band2);    _crs_sat_pco += t_gtriple(_rot_matrix * pco.crd_cvect());}if (rec_pcv) {t_gtriple pco(0, 0, 0);    rec_pcv->pcoR(_crt_obs, pco, _observ, band1, band2);    _crs_rec_pco += t_gtriple(_rot_matrix * pco.crd_cvect());}// 3. cmpObs() 中按当前观测频段再次计算 PCO 投影；相位观测还计算 PCVmodel += PCV(_corrt_sat_pcv, _corrt_rec_pcv, epoch, ...);

PCV() 这个函数名略有迷惑性：它对码和相位都返回当前频段的收发端 PCO 投影，只有相位观测才额外插值卫星端和接收机端 PCV。

4、多系统多频的偏差处理

在多系统多频 PPP-AR 中，需要同时处理多个 GNSS 系统、多个频率的观测值。不同系统采用不同的时间基准，不同频率信号在卫星和接收机硬件中的传播延迟也存在差异，这些差异表现为各种形式的偏差。本节逐一介绍 GREAT-PVT 中的偏差处理方式，包括 DCB（产品改正）、ISB/IFB（参数估计）、IFCB（模型/产品改正）、BDS2 SCB（模型改正）等。

GREAT-PVT 将这些偏差分为两大类：

• 外部产品改正类：通过 BIAS-SINEX、UPD、IFCB 等文件提供，直接改正到原始观测值上；
• 滤波器参数估计类：在滤波器中将偏差作为待估参数，随观测数据动态更新。

t_gallbias 只管理从 BIAS-SINEX/Bernese BIA 等文件读入的产品偏差（如 DCB/OSB）。ISB、IFB 不是存放在 t_gallbias 中的产品，而是 t_gallpar 状态向量里的滤波参数。观测值偏差在 t_gobsgnss::apply_bias() 中按配置选择 OSB 框架或 DCB 路径（t_gsatdata 继承该接口）；状态偏差则在时间更新和观测方程中处理。对于 BDS2 卫星端伪距偏差和 GPS L5 的 IFCB，代码另有专门模型函数。

1. 配置读取   └── gsetproc.cpp / gsetinp.cpp       ├── <bias> 文件路径       ├── <upd> 文件路径       ├── <ifcb_model> COR       └── <bds_code_bias_corr> true/false2. 偏差产品加载   └── GREAT_PVT.cpp       └── t_biasinex / t_biabernese           └── t_gallbias::_mapBias3. 观测预处理   └── gpvtflt.cpp::_prepareData()       └── gqualitycontrol::processOneEpoch()           ├── excludeBadObs()           └── t_gbds_codebias_cor::apply_IGSO_MEO()  [BDS2 SCB]4. DCB/OSB 改正   ├── gpvtflt.cpp::ProcessOneEpoch()   │   └── t_gobsgnss::apply_bias(_gallbias)  // 由 t_gsatdata 继承调用   │       ├── SGG_A 模式：apply_code_phase_bias() [OSB 框架]   │       └── 其它模式：apply_dcb()               [DCB]   └── gcombmodel.cpp::cmb_equ() 再次调用 apply_bias()       └── _dcb_correct_mark 已置位，直接返回，不重复改正5. 参数动态创建   └── gsppflt.cpp       ├── _syncSys() → GLO/GAL/BDS/QZS ISB       └── _syncIFB() → GPS/GAL/BDS/QZS IFB   参数预测   └── gpvtflt.cpp::_predictBias()       └── 为 ISB 和 IFB 初始化方差或增加相应系统的过程噪声6. 观测方程构建   └── gprecisebiasGPP.cpp::cmpObs()       ├── isbDelay()  → ISB + GLO_IFB       ├── ifbDelay()  → GPS/GAL/BDS/QZS IFB       ├── ionoDelay() → SION       └── windUp() / tropoDelay() / PCV()7. 模式 IFCB   └── gcombmodel.cpp       └── ifcbDelay() [GPS L5 IFCB]8. UPD 产品加载   └── GREAT_PVT.cpp       └── t_upd           └── t_gupd::_upd9. 模糊度固定中的 UPD 改正   └── gambiguity.cpp       ├── _getSatRef()      [参考星 UPD 质控]       ├── _applyUpd()       [RAW_ALL/RAW_MIX 模式]       ├── _applyWLUpd()     [WL/EWL 单独应用]       └── _getSingleUpd()   [单星 UPD 查询]

下文逐项介绍各类偏差的处理方式。

1、DCB

DCB（Differential Code Bias）表示两个信号之间的码偏差差值。在 GREAT-PVT 中，apply_dcb() 将 DCB 改正分为两步：频内 DCB 改正（同一频点内不同调制码信号间的偏差）和 频间 DCB 改正（不同频点间的偏差）。两者都基于 t_gallbias::get(epo, sat, obs1, obs2) 查询到的差分偏差值。

ProcessOneEpoch() 和 IF/UDUC 组合模型都会调用 apply_bias()，第一次完成改正后 apply_dcb() 将 _dcb_correct_mark 置为 true，后续调用检查到该标志后直接返回，因此粗差迭代中不会重复累计 DCB。

频内 DCB 改正是针对 GPS 和 GLONASS 的特定码信号，将其改正到与精密钟差产品一致的参考信号。例如，精密钟差产品通常基于 C1W 和 C2W，而实际观测可能是 C1C、C2L、C2X、C5Q 等，因此需要加上这些信号与参考信号之间的 DCB。GPS 的频内 DCB 改正逻辑如下：

if (gsys == GSYS::GPS){switch (gobs_type.gobs())    {case C1C:  dcb_value = allbias->get(gepo, gsat, C1W, C1C); break;case C1P:case C1Y:case C1W:  dcb_value = 0.0; break;  // 本身就是参考信号case C2C:  dcb_value = allbias->get(gepo, gsat, C2W, C2C); break;case C2L:  dcb_value = allbias->get(gepo, gsat, C2W, C2L); break;case C2X:  dcb_value = allbias->get(gepo, gsat, C2W, C2X); break;case C2P:case C2Y:case C2W:  dcb_value = 0.0; break;  // 本身就是参考信号case C5Q:  dcb_value = allbias->get(gepo, gsat, C5X, C5Q); break;case C5X:  dcb_value = 0.0; break;  // 参考信号default:   dcb_value = 0.0; break;    }}

GLONASS 的频内 DCB 改正逻辑类似：

elseif (gsys == GSYS::GLO){switch (gobs_type.gobs())    {case C1C:  dcb_value = allbias->get(gepo, gsat, C1P, C1C); break;case C1P:  dcb_value = 0.0; break;case C2C:  dcb_value = allbias->get(gepo, gsat, C2P, C2C); break;case C2P:  dcb_value = 0.0; break;default:   dcb_value = 0.0; break;    }}

对于 Galileo、BDS、QZSS，代码中没有单独的频内 DCB 分支，其频内差异（如 E1 上的 C1C/C1X、B1I 上的 C2I/C2X）没做处理。

我的理解：之所以处理 GPS 和 GLONASS 的频内偏差是因为它们最主流的伪距观测数据（C1C）和精密钟差基准所用通道不一致，必须需要改。BDS、Galileo 如果用非精密钟差基准的通道，也是需要改频内DCB消除伪距观测值和精密星历之间的差异。

频间 DCB 改正是将不同频点的码观测统一到与精密钟差产品绑定的两个参考频点组合上。apply_dcb() 的核心流程为：

1. 确定参考频点 band1（第一频率）和 band2（第二频率，BDS 特殊处理为第三频率 B3I）；
2. 计算无电离层组合系数 alfa12、beta12；
3. 获取频间 DCB inter_dcb_12 = interFreqDcb(*allbias, band1, band2)；
4. 逐码观测应用改正：
• band1 信号：-beta12 * inter_dcb_12
• band2 信号：+alfa12 * inter_dcb_12
• 其他任意频点 band_x：alfa12 * DCB(band1,band_x) + beta12 * DCB(band2,band_x)

if (gobs_type.band() == band1)    dcb_value += -beta12 * inter_dcb_12;elseif (gobs_type.band() == band2)    dcb_value += +alfa12 * inter_dcb_12;else{    GOBSBAND band_x = gobs_type.band();    inter_dcb_1x = this->interFreqDcb(*allbias, band1, band_x);    inter_dcb_2x = this->interFreqDcb(*allbias, band2, band_x);    dcb_value += alfa12 * inter_dcb_1x + beta12 * inter_dcb_2x;}obs_value += dcb_value;this->resetobs(obs_type, obs_value);

当处理 Galileo 或 BDS 的四频、五频数据时，band_x 可以是第 4 或第 5 频点。apply_dcb() 通过 interFreqDcb() 动态查询 band1-band_x 和 band2-band_x 的 DCB 值，然后以 alfa12/beta12 加权组合，从而将任意额外频点的码偏差改正到与精密钟差产品一致的参考频点组合上。若外部 DCB 产品未提供对应频间 DCB，则该部分改正量为 0，残差由后续 IFB 参数估计吸收。

卫星端频间 DCB 也可以不改正，由电离层参数吸收，并进一步联合 VTEC 建模估计 卫星 DCB。GREAT-PVT 当前会施加频间 DCB，因此直接用其输出电离层参数反估同一 卫星 DCB 时，结果会接近零。

示例数据的DCB改正是用CAS0产品，CAS0没有做逐频点的天线改正，和用户端算法不完全匹配，如果服务端和用户端都用做了天线改正的CAS1产品应该效果会更好（CAS团队有论文说用CAS1产品定位精度优于CAS1，他们的钟差现在是输入CAS1-DCB算来的）。

2、UPD

UPD（Uncalibrated Phase Delay）是卫星和接收机端相位硬件延迟中不能被钟差、轨道产品完全吸收的小数部分。在 PPP-AR 中，浮点解估计的模糊度包含了这些未校准延迟，直接取整固定会导致系统性偏差。通过引入外部 UPD 产品，可以将星间单差模糊度的小数部分预先扣除，使模糊度恢复整数特性。

UPD 文件由 src/LibGREAT/gcoders/upd.cpp 解码，数据存储在 t_gupd 的 _upd[UPDTYPE][epoch][prn] 结构中。XML 配置为 <ambiguity><upd_mode>upd</upd_mode></ambiguity>，在 GREAT_PVT.cpp 中通过 IFMT::UPD_INP 加载并通过 Add_UPD() 注入各测站滤波器。具体用法在第四章中详细介绍。

3、OSB

OSB（Observable-specific Signal Bias）直接给出每个信号相对于系统时的绝对偏差。当输入武大 SGG_A 产品时，apply_code_phase_bias() 按观测类型逐个获取 WHU 相位偏差产品并计算无电离层组合偏差。需要注意的是，当前代码中组合偏差仅通过 std::cout 输出，实际未写入观测值（resetobs 仍使用原始值），该函数可视为 OSB 改正的预留框架。若想拓展实现基于 OSB 的 PPP-AR，可以在此完善。

4、IFCB

GPS Block IIF 及后续卫星的 L5 频段卫星钟差与 L1/L2 不一致，存在约 0.1～0.5 m 的系统性偏差，称为 IFCB。由于精密钟差产品通常基于 L1/L2 无电离层组合估计，使用 L5 观测时需要额外改正其相对 L1/L2 钟差的系统偏移。BDS/GAL 中这种偏差量级小，所以通常只对GPS做改正。GREAT-PVT 通过产品处理 IFCB，XML 配置项 <ifcb_model>COR</ifcb_model> 用于启用改正。

<ifcb_model>EST</ifcb_model> 模式尚未实现。

IFCB 产品通过 src/LibGREAT/gcoders/ifcb.cpp 解码，文件格式与 UPD 类似，包含 EPOCH-TIME 和每颗卫星的 value/sigma/npoint。数据存储在 t_gifcb 中。gprecisebias.cpp::ifcbDelay() 负责查询产品并返回改正量：

doublet_gprecisebias::ifcbDelay(t_gsatdata &satdata, t_gifcb *ifcb, OBSCOMBIN obscombin){double L = 0.0, ifcb_GPS = 0.0;if (ifcb == nullptr && _gifcb == nullptr) returnfalse;if (ifcb == nullptr) ifcb = _gifcb;    t_gtime epoch = satdata.epoch();string sat = satdata.sat();    GOBSBAND b1 = t_gsys::band_priority(satdata.gsys(), FREQ_1);    GOBSBAND b3 = t_gsys::band_priority(satdata.gsys(), FREQ_3);double l1 = satdata.wavelength(b1);double l3 = satdata.wavelength(b3);double c2 = l1 / l3;double c1 = 1.0 / (1.0 - c2 * c2);    c2 = c2 * c2 * c1;  // L1/L3 无电离层组合中 L3 系数的幅值（真实 L3 系数为 -c2）if (satdata.gsys() == GSYS::GPS)    {        one_epoch_ifcb epoch_ifcb = ifcb->get_epo_ifcb(epoch);if (epoch_ifcb.find(sat) == epoch_ifcb.end()) returnfalse;if (epoch_ifcb[sat]->sigma > 0.2 || epoch_ifcb[sat]->npoint <= 2)returnfalse;        ifcb_GPS = epoch_ifcb[sat]->value;if (obscombin == OBSCOMBIN::IONO_FREE)            L += ifcb_GPS;              // 直接加到 IF 组合观测残差elseif (obscombin == OBSCOMBIN::RAW_ALL)            L += ifcb_GPS / (-c2);      // 转换到 RAW L3 相位层面    }return L;}

RAW_ALL 模式下，ifcbDelay() 将 L1/L3 IF 组合层面的 IFCB 偏差，按该组合中 L3 的系数 (-c_2) 进行缩放，等效归化到原始第三频相位 (L_3) 上作为改正量：

// IONO_FREE 模式：构建 L1/L3 IF 组合相位时直接加入if (_ifcb_model == IFCB_MODEL::COR && _frequency == 3 &&    gobs1.is_phase() && _freq_index[gsys][b2] == FREQ_3 && obsdata.gsys() == GPS){    l_IF += dynamic_cast<t_gprecisebias *>(_bias_model.get())->ifcbDelay(obsdata, nullptr, OBSCOMBIN::IONO_FREE);}// RAW_ALL 模式：对第 3 频点（L5）的原始相位改正if (_ifcb_model == IFCB_MODEL::COR && _frequency >= 3 && freq == FREQ_3){    tempL.l[0] += dynamic_cast<t_gprecisebias *>(_bias_model.get())->ifcbDelay(obsdata, nullptr, OBSCOMBIN::RAW_ALL);}

5、BDS2 SCB

研究发现，BDS2 卫星存在一种可能由多路径引起的卫星端的伪距观测值系统偏差，称为北斗卫星伪距偏差，也称北斗伪距多路径延迟。该伪距偏差只存在于伪距观测值，和卫星高度角密切相关，可基于观测值的 MP 组合得到：

MP 组合是无电离层无几何组合观测值，对其平滑后可分离出伪距多路径噪声。Wanninger 基于卫星高度角，采用分段线性拟合的方法建立了节点间隔 10 度的北斗二号 IGSO 和 MEO 卫星伪距偏差的经验模型，如下表所示：

for (auto obs_type : obs_vec){if (!t_gobs(obs_type).is_code()) continue;    GOBSBAND b = t_gobs(obs_type).band();double obs_P = data->obs_C(obs_type);    obs_P += Band_cor[b];    data->resetobs(obs_type, obs_P);    data->setrangestate("bds_code_bias", true);}

6、ISB

ISB 是不同 GNSS 系统时间基准与公共接收机钟参数之间的差异。_syncSys() 对观测到的 GLO/GAL/BDS/QZS 分别管理 GLO_ISB、GAL_ISB、BDS_ISB、QZS_ISB；若数据中只有某一个非 GPS 系统，则不会为该系统创建 ISB。多系统情况下，只要某非 GPS 系统出现且不是唯一系统，代码就会创建相应 ISB，并不显式要求 GPS 必须存在。在常见的 GPS+多系统解算中，可以把 GPS 钟理解为基准、其它系统 ISB 相对它估计；若没有 GPS 却同时处理多个非 GPS 系统，当前创建逻辑仍会给各系统都加 ISB，需要另行检查参数秩亏问题。

ISB 参数在 gsppflt.cpp::_syncSys() 中创建和动态管理，在 gpvtflt.cpp::_predictBias() 中按随机游走模型预测。首次创建或重新初始化时状态置零并赋初始方差；之后在首轮迭代（_cntrep == 1）按系统增加 _gloStoModel->getQ()、_galStoModel->getQ() 等过程噪声。

7、IFB

IFB 是同一接收机不同频点间伪距硬件延迟的差异。当出现额外频率码观测时，gsppflt.cpp::_syncIFB() 按 _freq_index 的逻辑频率序号动态创建参数：GPS/QZSS 只为 FREQ_3 创建 IFB_GPS/IFB_QZS，Galileo 和 BDS 则可为 FREQ_3/4/5 创建三组 IFB。具体对应哪个物理频段由 XML 频段映射决定，不应固定写成某一组信号名。在设计矩阵中，相应伪距的 IFB 偏导为 1.0，相位偏导为 0。

IFB 并非创建后保持常数。gpvtflt.cpp::_predictBias() 首次为它设置很大的初始方差 3000²；后续在 _cntrep == 1 时，分别用 GPS/GAL/BDS/QZS 对应随机模型的 getQ() 增加协方差。因此当前实现也是随机游走式预测。

8、GLONASS

GLONASS 采用 FDMA（频分多址）体制，各卫星的载波频率互不相同，因此接收机端频间偏差与卫星 PRN 相关，无法通过星间差分完全消除。理论上可额外估计卫星级 GLO_IFB/IFCB/IFPB 参数，或由同型号接收机的 UPD/DCB 产品吸收。GREAT-PVT 在 gpar.h 中定义了 GLO_IFB、GLO_IFCB、GLO_IFPB 三种参数类型，gpar.cpp 和 gbiasmodel.cpp 也给出了相应的偏导计算；但当前 _syncIFB() 中没有为 GLONASS 创建这些参数的代码，因此尚未形成完整可用的估计流程。

四、PPP-AR 算法

PPP-AR（PPP with Ambiguity Resolution）在浮点解的基础上，通过引入外部 UPD 产品恢复模糊度的整数特性，并利用 LAMBDA 整数搜索和级联固定策略获得固定解。本章依次介绍 UPD 改正模型、EWL/WL/NL 级联固定策略、模糊度固定的具体流程（构建双差、预筛选、整数搜索、固定检验、虚拟观测约束）以及质量控制体系。

1、UPD 模型

载波相位观测值中不仅包含整周模糊度 ，还混杂了卫星端和接收机端的相位硬件延迟（分别记为  和 ）。这些硬件延迟来源于信号在卫星/接收机内部电路和天线中的传播时延，量级可达纳秒级，且无法被绝对校准。精密钟差产品通常基于无电离层组合伪距估计，吸收的是伪距硬件延迟的 IF 组合值，因此相位硬件延迟中与伪距不同的部分会被浮点解系数，导致其偏离整数。这部分未被钟差产品吸收的相位小数延迟，统称为 UPD（Uncalibrated Phase Delay，非校准相位延迟）。

设测站  观测卫星  的浮点模糊度（已扣除 PCO/PCV、潮汐等已知改正）可分解为：

其中  为真实整数模糊度， 为接收机端 UPD（与卫星无关）， 为卫星端 UPD（与测站无关）。该方程同时包含  和  两类未知量，直接求解秩亏。PPP-AR 的处理策略是服务端通过网解估计并分离卫星 UPD，向用户播发以消除卫星端相位小数偏差；用户端采用星间单差进一步消除接收机 UPD，从而实现模糊度整数解恢复：

具体推导如下。对于同系统、同频率的两颗卫星  和 ，分别写出非差浮点模糊度：

第一步：星间单差，消除接收机端 UPD。两式相减，公共的  被抵消：

第二步：卫星端 UPD 产品改正。从 UPD 产品文件中读取  和 ，将差分改正量  加到浮点单差模糊度上：

改正后的  即为星间单差整数模糊度，可以直接取整固定（EWL/WL）或进入 LAMBDA 搜索（NL）。

UPD 产品在代码中的时间索引方式如下：

2、EWL/WL/NL 级联

级联模糊度固定是 PPP-AR 中提高固定成功率、缩短收敛时间和增强可靠性的经典策略。其基本思想是按照波长从长到短的顺序逐层固定：首先固定波长最长、最易固定的超宽巷 EWL 模糊度，然后利用已固定的 EWL 约束来固定宽巷 WL 模糊度，最后利用已固定的 WL 约束来固定窄巷 NL 模糊度。EWL 波长可达数米，WL 波长通常为 0.75～0.86 m，二者一般比 NL 更容易固定；NL 波长约为分米量级，对噪声和残余模型误差更敏感，但是固定后精度更高。

对于 BDS-3、Galileo 等四频/五频系统，GREAT-PVT 在三频 EWL（L2–L3）基础上还引入了 EWL24（L2–L4）和 EWL25（L2–L5）两级超宽巷，三者均以 L2 为公共锚点。_amb_resolution() 按 EWL → EWL24 → EWL25 → WL → NL 的调用顺序执行，每固定一层就通过虚拟观测约束更新 _filter，使下一层获得更紧的浮点解初值与协方差。代码里的 "EWL" 不带后缀时特指 L2–L3 组合，"EWL24"/"EWL25" 才分别对应 L2–L4、L2–L5；它们的实际波长与具体系统的频段配置有关，但总体上远长于 WL/NL，因而更容易固定。

下表给出 GREAT-PVT 样例数据  PPPFLT_2023305  UPD 产品配套使用的 GEC 系统信号映射：

几点说明：

• 示例数据中 BDS 没有用 B1C 而用到了 B2ab 频点；
• 公共锚点是 FREQ_2：所有 EWL/EWL24/EWL25 都以 FREQ_2 为公共锚点，这样逐级固定后才能把整数约束传递到 WL/NL。
• 名称不代表波长顺序：EWL / EWL24 / EWL25 是按频率序号命名，不是按波长从长到短排序。例如 Galileo 中 EWL24 (~19.5 m) 最长，而 EWL25 (~2.93 m) 反而比 EWL (~9.77 m) 短。
• 与 UPD 产品绑定：UPD 产品（upd_ewl_*、upd_ewl24_*、upd_ewl25_*、upd_wl_*、upd_nl_*）是依据上述特定信号组合生成的，XML 中的 <band>/<freq> 配置必须与 UPD 产品一致，不能随意更改（这是用UPD产品做PPP-AR的一大缺点，用OSB就自由多了）。

GREAT-PVT 在 _amb_resolution() 中依次调用 processBatch("EWL")、processBatch("WL") 和 processBatch("NL")，并在 processBatch() 内部完成星间双差模糊度构建、UPD 改正、预筛选、独立子集选择以及虚拟观测约束等步骤（在下一节详细介绍）。每一层 processBatch() 都在同一个 _filter 构造虚拟观测：EWL/EWL24/EWL25 分别约束相应频率对，WL 约束宽巷组合，随后调用 gflt->update() 更新完整状态和协方差；下一层再从更新后的 _filter 读取参数与协方差。

if ((_observ == RAW_ALL || _observ == RAW_MIX) && _frequency >= 3) {    _ambfix->processBatch(_epoch, _filter, "EWL");}if ((_observ == RAW_ALL || _observ == RAW_MIX) && _frequency >= 4) {    _ambfix->processBatch(_epoch, _filter, "EWL24");}if ((_observ == RAW_ALL || _observ == RAW_MIX) && _frequency >= 5) {    _ambfix->processBatch(_epoch, _filter, "EWL25");}if (_observ == RAW_ALL || (_observ == RAW_MIX && _frequency >= 2)) {    _ambfix->processBatch(_epoch, _filter, "WL");}int nlfix_valid = _ambfix->processBatch(_epoch, _filter, "NL");

3、模糊度固定

本节深入 t_gambiguity 类的内部实现，从 processBatch() 调度入口开始，逐项展开星间双差构建、UPD 改正、预筛选、LAMBDA 搜索、固定检验以及虚拟观测约束回带的完整流程。

1. t_gambiguity::processBatch

t_gambiguity::processBatch() 是 GREAT-PVT  PPP-AR 级联固定的核心调度入口，由 _amb_resolution() 按 逐层调用；它接收当前历元、滤波器工作副本和固定模式，依次完成星间双差模糊度构造、UPD 差分改正应用、近整数预筛选、独立子集选择，并在 NL 层调用 LAMBDA 整数搜索，若固定成功则把整数解以强权重虚拟观测形式回注滤波器，从而把当前层的整数约束传递给下一层，实现多频模糊度的逐级固定。流程如下：

1. 构建星间双差模糊度：将滤波器中的模糊度参数两两配对，组成同一测站、同系统、同频点的两颗卫星之间的双差模糊度（消除接收机端 UPD）。代码通过遍历所有 par_type::AMB_Lx 参数，要求同系统、同测站、同频点、共视时间满足 _min_common_time，构造 t_dd_ambiguity 对象。
2. 应用 UPD 改正：
• EWL/WL 模式调用 _applyWLUpd()，扣除卫星端 WL/EWL UPD；
• NL 模式调用 _applyUpd()，扣除 NL UPD，同时恢复已固定的 WL/EWL 结果。
3. 预筛选：
• UDUC 的 EWL/WL：_fixAmbWL(mode) 直接判断 |rwl − round(rwl)| < maxdev（EWL 使用相应 rewl）。函数中虽有 _fix_epo_num[sat1/2] > 20 且 _lock_epo_num[sat1/2] > 200 时放宽到 maxdev × 1.1 的分支，但 EWL/WL 不运行 LAMBDA，其 amb_cmn Ratio 保持 0；在示例的正 Ratio 阈值下，这两层不会进入递增 _fix_epo_num 的成功分支，因此该放宽条件通常不可达。
• NL（UDUC）：_fixAmbUDUC() 继承已保存的 EWL/WL 标志，再对目标 NL 量执行同类 maxdev 近整数判断。该路径不调用 t_gbdeci，也不使用 maxsig/alpha 综合决策；t_gbdeci 出现在 IF 等其它分支。
4. 选择独立子集：_selectAmb() 按弧长降序排列，剔除线性相关的 DD（通过 _checkAmbDepend 正交基检验），确保进入 LAMBDA 的模糊度相互独立。
5. LAMBDA 搜索（仅 NL）：调用 _ambSolve() 进行整数搜索；EWL/WL 跳过 LAMBDA。
6. 虚拟观测约束：固定成功后，NL 调用 _addFixConstraint()（p0 = 1E9），WL/EWL 调用 _addFixConstraintWL()（p0 = 1E6，EWL 为 1E4），将整数解作为强约束注入滤波器。

2. 模糊度预筛选

在进入 LAMBDA 搜索之前，会先对模糊度做一次预筛选，目的是剔除明显不可固定的候选、降低后续搜索空间。代码中定义了三种模糊度预筛选阈值：

• Taper = ，反映浮点值靠近整数的程度，随偏差增大而减小，随 sigma 增大而减小。
• C 为基于互补误差函数的置信度积分，反映浮点值偏离整数的概率。

在 UDUC 模式下只用到了 maxdev，EWL/WL 由 _fixAmbWL(mode) 按与最近整数的偏差直接判定，NL 由 _fixAmbUDUC() 预筛，然后进入 _selectAmb() 与 _ambSolve() 的整数搜索。

3. LAMBDA

LAMBDA 是 Teunissen 提出的整数最小二乘估计方法，通过整数去相关显著压缩实际搜索空间、提高搜索效率。其核心流程包括：

1. 对协方差矩阵进行 LDLᵀ 分解（FMFAC6()）；
2. 通过整数高斯变换和条件交换执行 Z 变换，以减少搜索次数（SRC1i()）；
3. 在变换后的空间中由 FI71() 进行整数候选搜索；pBootStrapping() 计算成功率，不是搜索函数；
4. 最后将整数解逆变换回原始空间；

4. 部分模糊度固定

观测值质量有好有差，尤其是当卫星几何分布不佳或大气条件复杂时，并非所有模糊度都能可靠固定。部分模糊度固定策略通过剔除低置信度模糊度，仅固定高置信度子集，在保证可靠性的同时提高固定率。GREAT-PVT 的 _ambSolve() 实现了三级剔除策略：

每次剔除后重新运行 LAMBDA 搜索，直到满足 Ratio 阈值；如果剩余数量小于等于 _part_fix_num 时仍不满足 Ratio 阈值，则判定为模糊度固定失败。

5. 模糊度固定检验

模糊度整数搜索会产生多组候选解，需要可靠的检验指标来判断最优候选是否可信。GREAT-PVT 采用 Ratio Test 作为主要检验手段，辅以 Bootstrapping 成功率 作为二次验证。两者核心实现均在 src/LibGREAT/gambfix/glambda.cpp 的 t_glambda 类及 src/LibGREAT/gambfix/gambiguity.cpp 中。

Ratio Test（比率检验）

Ratio Test 是 LAMBDA 方法中最经典的模糊度确认检验。其基本原理是：LAMBDA 整数搜索会找出按残差二次型从小到大排序的多组候选整数解，如果最优解显著优于次优解，说明最优解具有足够的置信度；反之，若最优解与次优解的残差接近，则可能存在多组竞争解，最优解不可信。

设最优整数候选的残差二次型为 、次优为 ，Ratio 值定义为：

当  大于预设阈值（如 2.0 或 3.0）时，认为最优解显著占优，固定成功；否则拒绝固定该组模糊度。

在 GREAT-PVT 的实现中，_lambdaSearch() 调用 LAMBDA4() 时指定 maxcan = 2，只搜索最优和次优共两个候选：

// gambiguity.cpp::_lambdaSearch()constint maxcan = 2;lambda->LAMBDA4(maxcan, namb, Q, fbias, &ncan, &ipos, cands, disall, boot);// 计算 Ratioif (double_eq(disall[1], 0.0))    ratio = disall[0];       // 仅 1 个候选时，Ratio = Ω₁ (一般 > 阈值则通过)else    ratio = disall[1] / disall[0];  // Ω₂ / Ω₁

其中 disall[0] 和 disall[1] 由 FI71() 整数搜索与 COLLECTs()/STOREs() 候选收集后排序得出。阈值 _ratio 从 XML 配置读取（<lambda_ratio>），示例配为 2.0。若 ratio < _ratio，则清空候选并判定固定失败。

在部分模糊度固定循环中，每轮剔除低置信度模糊度后重新执行 LAMBDA4 并计算 Ratio，直到 Ratio 满足阈值或剩余模糊度数不大于 _part_fix_num。

Bootstrapping 成功率

Bootstrapping 成功率（Bootstrapping Success Rate）用于评估模糊度整数固定的可靠性，作为 Ratio Test 的附加验证。未配置 <boot> 时返回 -0.001，等效于不启用该门限；配置正阈值后程序初始化 _bootfile。不过 boot 文件只在 NL 搜索最终同时通过 Ratio 与 Bootstrap 时写入，并非每个历元都写。

在 LAMBDA 的 Z 变换（降相关） 之后，模糊度协方差矩阵被近似对角化。设变换后第  维的方差为 （即  分解中的对角元 ），则该维度的 Bootstrapping 成功概率 为：

其中  为标准正态分布 CDF。

整体 Bootstrapping 成功率是各维度成功概率的 连乘积：

逻辑含义：逐维按"最近整数"取整时，所有维度均不翻边的联合概率。

在标准 LAMBDA 流程中，完成 LDLᵀ 分解 (FMFAC6) 和 Z 变换 (SRC1i) 后，立即计算 Bootstrapping：

// glambda.cpp:534-539*boot = 1.0;for (i = 0; i < iN; i++){    dt = sqrt(pdD[i]);           // Z变换后第i维的标准差    *boot *= pBootStrapping(dt); // 查表并累乘}

• pdD[i]：Z 变换后协方差矩阵的对角元（方差）
• dt：该维度的标准差 

程序采用 预计算查找表 + 线性插值，避免运行时重复计算正态积分：

// glambda.cpp:651-721doublet_glambda::pBootStrapping(constdouble &sig){constdouble BootSt_2Phi_sigRecip_1[461] = {1.00000000000, 1.00000000000, 1.00000000000, 1.00000000000,1.00000000000, 0.99999999999, 0.99999999999, 0.99999999999,0.99999999999, 0.99999999998, /* ... 共461个 ... */0.90441929545    };// ...if (sig >= 0.3)        dRes = 0.9;elseif (sig <= 0.07)        dRes = 1.0;else    {        dt = (sig - 0.07) * 2000.0;  // 步长 0.0005        i  = int(dt);        j  = i + 1;        dRes = (j - dt) * BootSt_2Phi_sigRecip_1[i]             + (dt - i) * BootSt_2Phi_sigRecip_1[j];    }return dRes;}

6、固定解回带

模糊度固定成功后，程序将整数模糊度作为强约束虚拟观测加入滤波器，重新求解得到当前历元的固定坐标、钟差和其它参数。

代码中通过 _addFixConstraint() 构建固定模糊度虚拟观测，NL 使用 p0 = 1E9；_addFixConstraintWL() 中 WL 使用 1E6、EWL 使用 1E4。按 IF、RAW_ALL/RAW_MIX 及 GLONASS/非 GLONASS 情况计算 Ba、Bb、integer 和 flt：

// src/LibGREAT/gambfix/gambiguity.cppbool t_gambiguity::_addFixConstraint(t_gflt *gflt){double dl = 0, flt = 0, integer = 0, Ba = 1, Bb = -1;double p0 = 1E9;  // NL 约束使用极大权重for (auto itdd = _DD.begin(); itdd != _DD.end(); itdd++)    {// 跳过尚未固定的 DDbool wl_fix = ((_obstype == gnut::OBSCOMBIN::RAW_MIX && !itdd->isSngleFreq) ||                       _obstype != gnut::OBSCOMBIN::RAW_MIX);if ((!itdd->isWlFixed && wl_fix) || !itdd->isNlFixed)continue;// ... 按 IF / RAW_ALL / RAW_MIX 及 GLONASS 区分计算 integer 和 flt        dl = integer - flt;vector<pair<int, double>> B;        B.push_back(make_pair(index_sat1, Ba));        B.push_back(make_pair(index_sat2, Bb));        Matrix B_mat;        SymmetricMatrix P_mat;        ColumnVector l_mat;        t_gfltEquationMatrix virtual_equ;        t_gobscombtype type;        virtual_equ.add_equ(B, p0, dl, _site, ..., false);        virtual_equ.chageNewMat(B_mat, P_mat, l_mat, gflt->npar_number());        gflt->resetQ();                    // 每次加虚拟观测前重置        gflt->add_virtual_obs(B_mat, P_mat, l_mat);    }    gflt->update();returntrue;}

这些约束更新的是 _filter 中缓存的当历元解。_processEpoch() 在 AR 前已经把浮点参数、浮点改正数和浮点协方差复制进 _filter；AR 在这个对象上施加约束，固定成功后用 _ambfix->getFinalParams() 形成 _param_fixed 并写出本历元结果。随后主状态更新仍使用 AR 前的局部浮点改正数 dx：

_filter->add_data(_param, dx, _Qx, _sig_unit, Qsav); // 保存当历元浮点参数/改正数/协方差_filter->add_data(A, P, l);                          // 加入当历元观测方程_filter->add_data(vtpv, nobs_total, npar_number);    // 加入 vtpv / 自由度_amb_resolution();                                   // 约束并输出当历元固定解if (_amb_state) _postRes(A, P, l, dx);              // 固定成功时重组固定解残差for (unsignedint i = 0; i < _param.parNumber(); i++)    _param[i].value(_param[i].value() + dx(_param[i].index)); // 主状态仍加浮点 dx

因此，固定后的 _param_fixed 和 _filter->Qx() 没有覆盖持久化的 _param、_Qx，不能说固定状态在下一历元直接延续或加快下一历元主滤波收敛。

processBatch() 每次开始都会 _DD.clear() 并重新构造当前历元的卫星间差分模糊度。_DD_previous[mode] = _DD 虽然在固定成功时被赋值，但全工程没有读取 _DD_previous 的代码；_DD_save = _DD 保存的是当前历元集合，随后由 _postRes() 查找 L2 对应项，并不是下一历元初始化来源。EWL/WL 整数映射、固定标志和连续固定计数等 AR 辅助量会跨历元保留，但这不等于固定参数与协方差回写。

4、质量控制

GREAT-PVT 的质量控制贯穿整个 PPP-AR 流程，从观测数据预处理、浮点解后验残差剔除，到 UDUC 模糊度近整数预筛和 NL 整数搜索检验，形成多层递进的质量保障体系。