简介：全球导航卫星系统（GNSS）技术中，载波相位测量是实现高精度定位的核心。RTK（实时动态差分）通过站间和星间双差消除误差，结合卡尔曼滤波估计状态，再经LAMBDA算法固定整周模糊度，最终获得厘米级坐标。这一技术广泛应用于测绘、自动驾驶和精密农业。在工程实践中，处理NovAtel等接收机的原始观测日志，需理解RANGEB数据格式、双频组合及GPS/BDS联合解算中的系统间偏差。本文基于课程实习案例，系统拆解从二进制日志解析到模糊度固定的完整链路，并分析常见调参与验证方法，帮助开发者提升固定率与定位精度。
1. 卫星导航算法课程实习：基于NovAtel双频RTK解算，从数据到厘米级固定解
拿到这个压缩包的人，大概率是正在做卫星导航算法与程序设计课程实习的同学。包里最常见的组合是：一组NovAtel接收机采集的原始观测日志、一份基线长度已知的静态数据、还有一个要求你从零写RTK解算程序的任务说明。NovAtel双频RTK解算这件事，听起来像是个黑匣子——打开软件就有厘米级坐标，但自己动手写一遍才发现，从二进制日志里抠出伪距和载波相位只是开始，后面还有周跳探测、双差方程、模糊度固定一连串关卡。这篇笔记就按课程实习里最常走的路径，把这条链路完整拆开，让新手能照着搭出可运行的解算流程，也让已经跑通基本流程的同学知道问题出在哪、参数怎么调、结果怎么验证。
2. 从NovAtel原始日志到双频双系统观测值：解析这一步决定后面所有事
RTK解算的输入不是坐标，不是单点定位结果，而是原始观测值——伪距、载波相位、多普勒和信噪比。NovAtel接收机输出的日志格式里，这几个量分布在不同的LOG消息中，解析这一步如果搞错单位或者搞错频点对应关系，后面整个滤波器都会跟着翻车。常见做法是让接收机同时输出原始观测量、星历和定位结果三类日志，然后写一个解析器把它们转成统一的内部结构。
2.1 先认清楚NovAtel的日志类型：RANGEB和BESTPOS各管什么
NovAtel接收机的日志体系里，和RTK解算直接相关的主要有三类。第一类是原始观测值日志，常见的有RANGEB和OBS，RANGEB按卫星输出伪距、载波相位、多普勒和信噪比，OBS则是按历元输出的通用观测值；第二类是星历日志，常用的是EPHEM或RAWEPHEM，里面有广播星历参数，用于计算卫星位置和卫星钟差；第三类是定位结果日志，比如BESTPOS，它只能用来做参考对比，不能作为解算输入。
RANGEB日志里每个卫星对应一条记录，字段包含卫星号、伪距、载波相位（单位是周期）、多普勒频移和信号载噪比等。伪距单位是米，载波相位单位是周期，这两个量后面要组合使用。课程实习里最常遇到的坑是：把NovAtel输出的载波相位直接当成米来用，或者搞混了相位和相位率。解析前先确认单位，这是最便宜的一步。
2.2 解析RANGEB：伪距、载波相位、多普勒的提取与单位换算
RANGEB的二进制格式看起来复杂，但实际解析时只需要按字节偏移读出几个关键字段。接收机手册里会给出每个字段的长度和单位，我一般会先把一条RANGEB消息完整打印出来对照一遍，确认卫星号、伪距、相位这几个字段的位置，再写正式解析器。
下面是一段按常见RANGEB字段定义写出的解析示意，实际字段偏移以接收机手册为准：
typedef struct { uint32_t sat_id; /* 卫星号，高8位是系统标识 */ double pseudo_range; /* 伪距，单位米 */ double carrier_phase; /* 载波相位，单位周期 */ double doppler; /* 多普勒，单位Hz */ double cn0; /* 载噪比，单位dB-Hz */ } novatel_rangeb_t; int parse_rangeb(const uint8_t *buf, int len, novatel_rangeb_t *out) { int offset = 0; memcpy(&out->sat_id, buf + offset, 4); offset += 4; memcpy(&out->pseudo_range, buf + offset, 8); offset += 8; memcpy(&out->carrier_phase, buf + offset, 8); offset += 8; memcpy(&out->doppler, buf + offset, 8); offset += 8; memcpy(&out->cn0, buf + offset, 4); offset += 4; /* 根据手册继续读冗余字段 */ return 0; }这段逻辑的关键是字段顺序和长度要按接收机的输出定义来。伪距直接就能用于伪距观测方程，载波相位乘以波长就能得到以米为单位的相位观测值：
phase_m = carrier_phase * wavelengthwavelength = 光速 / 频率。GPS L1约0.1903米，BDS B1I约0.1920米，混用会直接导致一米级偏差，这是常见错误里最隐蔽的一个。多普勒在RTK里不参与定位解算，但可以用来辅助周跳探测和速度估计，解析后建议先保留下来，后面会用得上。
2.3 把观测值组织成解算结构：一个可用的结构体与解析片段
解析出来的数据要按历元组织。RTK解算是以历元为单位的，每个历元包含基准站和流动站各自的观测快照。我一般会用两个结构体分别表示单站观测和双站观测：
typedef struct { double t_rcv; /* 接收机时刻 */ int sat_count; /* 当前历元卫星数 */ int sys[32]; /* 系统标识：1=GPS, 2=BDS */ int prn[32]; /* 卫星PRN号 */ double pr[32]; /* 伪距 */ double cp[32]; /* 载波相位，单位米 */ double cn0[32]; } obs_t; typedef struct { obs_t base; /* 基准站观测 */ obs_t rover; /* 流动站观测 */ double dt; /* 两站接收机时间差 */ } epoch_obs_t;组织好这个结构之后，后面的双差构建、周跳探测、滤波器更新都只需要遍历这个数组，不用再回原始日志里翻字段。这一步值得多花时间测试：先用两个接收机静止时的数据跑一遍，检查每个历元卫星数是否稳定、伪距是否在合理范围、相位是否连续，再进RTK主流程。解析层至少承载了后面一半的调试工作量。
3. RTK解算核心算法链：双差方程、双频组合、滤波与模糊度固定
观测值准备好之后，进入RTK算法主体。RTK的本质是利用两个接收机的观测值做差分，消除卫星钟差、接收机钟差和大部分大气延迟，然后通过载波相位的高精度测量值加上整周模糊度解算，最终得到厘米级相对位置。整个算法链按顺序可以做四件事：构建双差观测方程、选择双频组合、设计滤波器状态、固定模糊度并输出固定解。
3.1 双差观测方程：时钟误差是怎么消掉的
单差是站间差分，把同一颗卫星在基准站和流动站的观测值相减，接收机钟差被消掉。双差是在单差的基础上再在卫星间做一次差分，把卫星钟差也消掉。剩下的量是双差几何距离、双差电离层延迟、双差对流层残差和双差模糊度。
以GPS L1和BDS B1I为例，双差伪距和双差相位观测方程可以写成：
∇ΔP = ∇Δρ + ∇ΔI + ∇ΔT + ε_P ∇ΔL = ∇Δρ - ∇ΔI + ∇ΔT + λ·∇ΔN + ε_L其中∇Δρ是双差几何距离，∇ΔI是双差电离层延迟，∇ΔT是双差对流层残差，∇ΔN是双差整周模糊度。相位观测值此时是米单位，所以方程里直接乘了波长。由于双差模糊度天然是整数，这为后面的模糊度固定提供了基础。
构建双差方程时需要一个参考卫星。每颗卫星的观测值都要和同一颗参考卫星做差分，参考卫星通常选高度角最高的那颗，因为高度角高意味着大气延迟残差小、周跳概率低。参考卫星一旦发生周跳，所有双差观测都要重组，这个细节在实习代码里经常被忽略。
3.2 状态向量与卡尔曼滤波：位置、钟差、模糊度一起估
RTK解算的滤波器状态向量一般包含流动站位置、基准站和流动站之间的相对天顶对流层延迟、以及所有卫星的双差模糊度。实际代码中，单差模糊度和双差模糊度都有人用，熟练的做法是用双差模糊度直接当状态，因为滤波结束后做整数约束更方便。
卡尔曼滤波的预测和更新可以这样组织：
# 状态向量: [x, y, z, dTrop, N1, N2, ..., Nn] # 观测向量: 双差伪距 + 双差相位 def rtk_filter(epoch_obs, x_pri, P_pri, x_sat, x_base): H = build_dd_matrix(epoch_obs, x_sat, x_base) # 双差几何矩阵 R = build_obs_noise(epoch_obs) # 伪距/相位噪声 K = P_pri @ H.T @ inv(H @ P_pri @ H.T + R) x_post = x_pri + K @ (z_obs - H @ x_pri) P_post = (I - K @ H) @ P_pri return x_post, P_post这里的关键是H矩阵的列对应状态向量，前四列是位置和天顶对流层的偏导数，后面每一列对应一个双差模糊度，取值是±1或0，取决于该卫星在双差组合中与参考卫星的符号关系。噪声矩阵R的对角元素根据伪距和相位的先验方差设置，通常伪距方差设为0.3^2米，相位方差设为0.003^2米。
滤波器跑起来之后，模糊度的浮点解会逐步收敛。但浮点解精度不够，必须把双差模糊度约束为整数才能得到固定解。这一步常见做法是把浮点模糊度及其协方差矩阵交给LAMBDA算法，通过整数最小二乘搜索得到最优整数组合。
3.3 模糊度固定：LAMBDA搜索与固定解输出
LAMBDA是模糊度固定的标准方法。它先把模糊度向量做整数变换，让搜索空间更接近球形，再用树搜索找到使得残差平方和最小的整数组合。实习代码里不重新造轮子，而是调用现成的LAMBDA函数接口。
int lambda(int n, int m, const double *a, const double *Q, double *F, double *s, double *s2); /* n: 模糊度个数 m: 固定解的候选个数 a: 浮点模糊度向量 Q: 浮点模糊度协方差矩阵 F: 输出的固定模糊度向量 s/s2: 残差平方和 */固定之后要做一个比值检验，通常用次优解和最优解的残差平方和之比来判断固定是否可靠。实践中这个比值阈值设在2到3之间，低于阈值时宁可用浮点解，也不要强行固定。固定解的位置更新则是把模糊度代入法方程回代得到。
双频数据的优势在模糊度固定阶段非常明显。利用双频观测值可以构建宽巷组合，宽巷波长约86厘米，远大于单频波长，更容易固定宽巷模糊度，再回代求解窄巷模糊度。这一步能显著提高固定率，也是课程实习想让同学体会到双频价值的关键点。
4. GPS与BDS联合解算：双系统组合的模型细节与代码落点
单GPS系统的RTK在城市环境里经常因为卫星数不足而无法固定。加入BDS之后，可用卫星数增加，几何构型改善，固定率和收敛速度都会有明显提升。但双系统不是简单地把观测值堆在一起，两个系统的频点不同、信号结构不同，还有一个被忽略但影响巨大的系统间偏差参数需要处理。
4.1 GPS L1/L2与BDS B1I/B3I的频点差异怎么处理
GPS L1的频率是1575.42 MHz，L2是1227.60 MHz；BDS B1I是1561.098 MHz，B3I是1268.52 MHz。RTK解算用到的是双频消电离层组合或窄巷组合，频率不同，波长和电离层延迟系数都不同。代码里的波长表必须按系统分别设置，不能共用一个常数。
我一般会在代码里维护一张频点表：
/* 载波频率表，单位Hz */ static const double FREQ_GPS_L1 = 1575.42e6; static const double FREQ_BDS_B1I = 1561.098e6; static const double FREQ_BDS_B3I = 1268.52e6; double get_wavelength(int sys, int freq_band) { if (sys == SYS_GPS && freq_band == 1) return C / FREQ_GPS_L1; if (sys == SYS_BDS && freq_band == 1) return C / FREQ_BDS_B1I; return 0.0; }双频组合使用时，GPS的组合系数用L1/L2频率计算，BDS用B1I/B3I频率计算，两者不能混用。课程实习里如果把GPS的L2波长套到BDS B3I上，电离层残差和模糊度全都会出问题，且这类错误非常隐蔽，因为解算不会崩溃，只是精度和固定率上不去。
4.2 系统间偏差ISB：联合解算不能忽略的额外未知数
GPS和BDS联合解算时，双差消除的是单系统内的钟差。跨系统做差分时，两个系统的信号延迟基准不同，会引入一个随时间缓慢变化的系统间偏差。这个偏差如果不处理，相当于把两个系统的观测值强行拼在一起，会导致模糊度固定失败、残差异常。
常见做法是在状态向量里增加两个ISB参数，一个对应GPS到BDS的伪距偏差，一个对应相位偏差。很多实习代码用常数模型，静态场景下这是够用的；动态场景或接收机冷启动时，ISB会有跳变，建议用随机游走模型驱动。
# 状态向量扩展为: # [x, y, z, dTrop, ISB_p, ISB_c, N1, N2, ...] # ISB随时间变化: 方差按过程噪声缓慢放大把ISB估计出来之后，GPS和BDS的观测值才能真正融合成一个解算系统。这一步是双系统联合解算和单系统解算在算法上的本质区别，也是课程实习评分里拉开差距的地方。
4.3 联合解算的观测方程组织：双系统怎么放进同一组法方程
双系统观测方程的组织方式是把GPS和BDS的双差观测分别构建，然后拼进同一个H矩阵和同一个观测向量。每一颗卫星的模糊度状态是独立的，但位置和天顶对流层状态是共享的，ISB是BDS独有的。
下面是一段构建双系统H矩阵的示意：
int build_H_dd(epoch_obs_t *ep, double *pos_rover, double *pos_base, double *x_sat_gps, double *x_sat_bds, double *H, double *z, int *n_obs) { int n = 0; /* GPS双差方程 */ for (int i = 0; i < ep->gps_count; i++) { /* 计算双差几何矩阵行 */ H[n * state_dim] = dd_geometric(...); H[n * state_dim + 3] = trop_coef; /* 天顶对流层 */ H[n * state_dim + 4 + i] = 1.0; /* GPS模糊度 */ z[n] = dd_phase - dd_geom; n++; } /* BDS双差方程，注意ISB列 */ for (int j = 0; j < ep->bds_count; j++) { H[n * state_dim] = dd_geometric(...); H[n * state_dim + 4] = 1.0; /* ISB伪距/相位偏差 */ H[n * state_dim + 5 + gps_amb_count + j] = 1.0; /* BDS模糊度 */ z[n] = dd_phase - dd_geom; n++; } return n; }这里关键点是ISB参数的列位置要与模糊度的列位置错开，否则法方程会秩亏。整理好H矩阵后，后面的卡尔曼滤波更新和单系统完全一致。联合解算的好处在于模糊度固定的候选卫星更多、PDOP值更小，实践中BDS占比越高，窄巷固定成功率越快。
5. 实习中最容易踩的坑：数据质量、周跳漏检与低固定率的排查
RTK解算程序跑通一遍不难，难的是结果稳定。实习阶段最常见的翻车场景就集中在这几个现象上，每一个都有明确的排查路径。
5.1 现象一：固定解突然跳变回浮点解
接收机静止、环境良好，固定解跑了一阵子，突然跳到浮点解，过几秒又恢复。这个现象通常不是滤波发散，而是某个参考卫星发生了周跳。
原因是参考卫星的载波相位出现滑周，双差模糊度整体被打断，滤波器只能重新收敛。解决方法是加入周跳探测，最常用的就是GF组合：
GF = L1 - L2GF组合消除了几何距离和大部分大气延迟，只剩电离层残差和模糊度差。当GF在连续历元之间变化超过阈值，比如4厘米，就判定该卫星发生周跳，将其模糊度重置。
5.2 现象二：加了北斗反而定位精度变差
单GPS能固定，加入BDS之后固定率降低或者结果抖动变大。这通常是ISB没有估计造成的。
原因是GPS和BDS的观测值在系统间存在信号延迟差异，如果忽略这个偏差，相当于把两个系统的所有误差强行混在一起。解决方法是按前文所述引入ISB参数，或者退一步，在实习初期先用单系统跑通，再扩展双系统，不要一上来就同时调试两个系统的代码。
5.3 现象三：模糊度固定率只有百分之六七十
固定率上不去，先检查观测值本身的精度。低高度角卫星的多路径误差大，模糊度搜索空间被污染，固定率自然低。常见做法是把截止高度角从5度提到10度或15度，代价是卫星数略减，但固定率通常会有明显提升。
另外要检查伪距噪声的方差设置是否合理。伪距噪声设置过大会导致浮点模糊度收敛慢，过小会导致滤波过快锁死到错误值。我一般会把BDS GEO卫星的伪距噪声设得比GPS略大，因为GEO卫星位于赤道上空，在中高纬度地区的高度角偏低，多路径更严重。
5.4 现象四：残差图有明显的斜条纹
双差残差随时间呈线性或周期性漂移，这类问题基本可以锁定到两个原因。一是对流层残差没有被充分建模，基线超过10公里时尤其明显；二是模糊度被固定到了错误整数上，比值检验通过但残差内部互相吸收。
解决方法是先看基线和采样间隔，基线长则必须估计天顶对流层并在状态向量里建随机游走模型；如果残差条纹出现在短基线上，优先检查双差模糊度的取值符号是否和H矩阵一致，符号反了会导致残差按几何距离的倍数周期性变化。
6. 评估定位结果的验证技巧：固定率统计、残差绘图与参数复盘
解算程序写完，不能只看NMEA输出里的坐标就完事。课程实习的验收重点一般有三个：固定率、定位精度和收敛时间，这三个指标都需要一个独立的评估脚本。
固定率指固定解历元数占总解算历元数的比例，一般要求静态场景能达到95%以上。定位精度用固定解坐标和已知坐标的差值统计RMS来评估。收敛时间则看从开始解算到第一次固定解出现的时间间隔。写一个简单的Python脚本，逐历元读入解算结果，统计这些指标：
import numpy as np # 假设输出文件每行为: 时间,sys,pos_x,pos_y,pos_z,amb_state,ratio data = np.loadtxt('rtk_result.txt', delimiter=',') # 固定率 fixed_mask = data[:, 5] > 100 # 固定解标记 fixed_rate = fixed_mask.sum() / len(data) # 与真值对比 true_pos = np.array([x, y, z]) # 已知坐标 rms = np.sqrt(((data[fixed_mask, 2:5] - true_pos) ** 2).sum(axis=1).mean()) # 首次固定时间 first_fix_idx = np.argmax(fixed_mask) first_fix_time = data[first_fix_idx, 0] - data[0, 0] print(f'固定率: {fixed_rate:.2%}, 位置RMS: {rms:.3f}m, 首次固定: {first_fix_time:.1f}s')参数复盘也很关键。同一个数据，把噪声方差、过程噪声、截止高度角分别放大和缩小，固定率和RMS的变化趋势能直接反映你的滤波器参数是否合理。我第一次完整跑通双系统RTK时，固定率只有85%，反复检查代码没发现问题。后来把伪距方差从0.3改到0.5，固定率升到97%。看似细节的参数，其实决定了整个实习项目的最终评价。
最后给个我自己的习惯：每改一个参数，只动一个变量，跑完记录固定率和RMS，不要同时改多个参数。这样复盘才有依据，答辩时被问到“这个参数为什么这样设”也有底气。希望帮到你。
本文还有配套的精品资源，点击获取
