GPS/BDS双频RTK解算实战：基于NovAtel观测数据的模糊度固定

简介：面向卫星导航算法与程序设计课程实习的源码工程，以诺瓦泰尔接收机为平台，实现GPS与BDS双频RTK解算。压缩包共61个文件，以C++源码为主体，包含24个cpp实现文件、22个h头文件及5个hpp模板文件，并附带Visual Studio解决方案、XML配置等，整体仅93KB，结构紧凑。代码覆盖OEM原始报文解析、双频观测值处理、卫星星历解码、坐标与矩阵运算，以及模糊度解算、粗差探测等核心环节，能够帮助读者串联从观测量获取到定位结果输出的完整RTK流程。所有源码经过严格测试可直接运行，目录按功能清晰划分，适合有卫星导航基础的学生或工程师用于课程设计、算法验证与二次开发。目前已有152人学习，可作为理解双频RTK工程实现的实用参考。

1. 基于NovAtel的GPS/BDS双频RTK解算，这门课设到底在练什么

卫星导航算法课设做到“双频RTK解算”这一步，通常意味着前面伪距单点定位、载波观测量生成、广播星历读取这些环节已经跑通了。标题里把“基于NovAtel”和“GPS/BDS双频RTK”放在一起，本质是要你拿着真实接收机的观测文件，自己写一遍差分定位的完整链路：双频观测值组合、周跳探测、双差方程构建、模糊度浮点解、Lambda固定、固定解输出。这个过程不是调库调参，而是把RTK从原理变成可运行的程序。

选择NovAtel作为数据源有两个现实原因：一是它的接收机在测绘、工程测量领域存量很大，OEM板卡和接收机输出的原始观测数据格式公开，且支持同时输出GPS L1/L2和BDS B1/B2双频数据；二是NovAtel的日志格式相对规整，能直接拿到伪距、载波相位、多普勒和信噪比，省去先从RINEX里猜字段的麻烦。适合的人群很明确：正在做卫星导航课程设计或毕业设计的学生、需要把RTK算法从论文落到代码的工程师。

标题里的zip包大概率是课程实习的工程文件或某次实习的项目归档，但真正值钱的不是包里的成品，而是你能否从原始观测量开始，把双频RTK的每一步误差源、每一个矩阵、每一处参数设置都讲清楚、跑明白。下面这套方案按“理论先行、代码可抄、坑点可见”的思路展开。

2. 从NovAtel原始观测量到可用观测值：协议解析与双频数据组织

2.1 NovAtel日志格式与观测数据字段

NovAtel接收机通过串口或网口输出ASCII或二进制日志，RTK解算中最常用的是RANGEB（二进制格式的原始观测量）和RAWEPHEM（原始星历）。如果接收机配置为输出RINEX格式，也可以直接生成*.obs和*.nav文件，但对课程设计来说，直接从默认日志里抓数据更有“程序设计”的味道。

常见的接收机指令是：

log rangeb ontime 1 log rawephem ontime 60 log bestpos ontime 1

三个命令分别控制观测量、星历和定位结果的输出频率。ontime 1表示每秒输出一次，ontime 60表示每60秒输出一次星历。星历不需要高频输出，因为广播星历的有效期通常是2小时，二次差分和卫星位置计算都基于它。

2.1.1 RANGEB日志的核心字段映射

RANGEB日志以二进制形式输出，每条记录包含接收机时间、卫星编号、信号类型、伪距、载波相位、多普勒、信噪比等字段。以二进制解析时需要注意，NovAtel的二进制消息遵循“头+消息体+CRC”的结构，头部的Message ID字段能区分日志类型。

import struct def parse_rangeb_header(data): # NovAtel二进制消息头为28字节，前4字节为同步字0xAA 0x44 0x12 0x1C sync = data[0:4] if sync != b'\xaa\x44\x12\x1c': raise ValueError("同步字错误") # 头部字段: 消息长度(2字节)、消息ID(2字节)、消息类型(1字节)、端口(1字节) # 空闲(2字节)、接收机时间(8字节)、时间状态(1字节)、消息状态(1字节) msg_len, msg_id = struct.unpack('<HH', data[4:8]) recv_time = struct.unpack('<d', data[14:22])[0] return msg_id, msg_len, recv_time

解析时最重要的三个参数是：字节序为小端（NovAtel使用小端序，代码中<代表little-endian）、消息长度指后续数据长度而不是整条消息长度、接收机时间以秒为单位存储。解析时建议一次性读取整个文件，按同步字切分消息，逐条处理后丢弃CRC校验，不要边读边跳跃定位。

2.1.2 双频观测值的识别规则

GPS和BDS的观测值在RANGEB中通过“信号类型”字段区分。GPS L1对应1C或L1C，L2对应2L或L2W；BDS B1对应2I（老接口文档中为B1I），B3对应6I。如果接收机支持BDS B1/B2双频，则B2对应7I。这里有一个容易踩的坑：NovAtel在不同固件版本中对BDS信号类型的命名不完全一致，较老固件用2I和7I表示B1I和B2I，新固件可能增加1P等新信号类型。

def is_gps_l1(sig_type): return sig_type in ('1C', '1L', '1S') def is_gps_l2(sig_type): return sig_type in ('2L', '2W', '2S') def is_bds_b1(sig_type): # 老固件为2I，新固件可能出现1P等 return sig_type in ('2I', '1P') def is_bds_b2(sig_type): return sig_type in ('7I', '5P')

判断信号类型的依据是NovAtel文档中的Signal Type字段，课程设计代码里一般只要保留GPS L1/L2和BDS B1/B2四类信号即可。每个历元下同一颗卫星如果有多个信号频率，要分别存为独立观测值，不能用单一结构体存储。

2.2 伪距与载波相位的数据清洗

原始观测值不能直接用，要做三步清洗。第一，删除信噪比过低的观测值，一般阈值取30 dB-Hz，城市环境下可降到25；第二，删除伪距或载波相位为0或异常浮点值的记录；第三，检查载波相位的周跳标记，NovAtel的RANGEB中周跳标记是一个bitmask，但更可靠的做法是用GF组合（Geometry-Free组合）自行探测。

2.2.1 周跳探测的GF组合实现

GF组合是双频载波相位之差，表达式为：

GF = L1 - L2

其中L1和L2是以米为单位的载波相位观测值。由于双频路径延迟差异通常小于0.05米，GF组合在无周跳时是一条平滑曲线，相邻历元的突变量超过阈值即判定为周跳。阈值的经验取值是0.04米，在电离层活跃期可放宽到0.08米。

def detect_cycle_slip(prev_gf, curr_gf, threshold=0.04): """ 用GF组合检测周跳 prev_gf: 上一历元GF组合值 curr_gf: 当前历元GF组合值 返回True表示发生周跳 """ if prev_gf is None: return False return abs(curr_gf - prev_gf) > threshold

若周跳被检测到，该卫星当前历元的载波相位需在后续双差处理中做重置处理，即不再参与连续弧段构建，可以重新初始化模糊度。GF组合不能区分L1和L2哪个发生了周跳，但对于RTK解算来说，只需知道该弧段中断即可。

2.2.2 按PRN和频率组织观测集合

完成清洗后，把观测值组织成epoch -> prn -> frequency -> observation的四层字典结构。数据结构设计直接决定后续双差构建的代码复杂度。

observations = { "2024-03-15 08:00:00.000": { "G01": {"L1": {"pseudo_range": 20456890.123, "carrier_phase": 107253189.235, "snr": 42.5}, "L2": {"pseudo_range": 20456891.345, "carrier_phase": 83643218.456, "snr": 38.2}}, "C01": {"B1": {"pseudo_range": 20456892.234, "carrier_phase": 107253190.567, "snr": 40.1}} } }

这里G01是GPS卫星PRN1，C01是BDS卫星PRN1。按这个结构组织后，后续构建双差方程时可以在单个历元内按PRN遍历，也可以跨历元按连续弧段遍历。存储时保留全部原始值，不要预先扣除任何误差项，误差处理放在观测方程构建阶段。

3. 双频RTK解算的数学框架：双差方程与模糊度浮点解

3.1 为什么要用双差模型

RTK的核心是消除或削弱误差源。单差是站间差分，可以消除卫星钟差；双差是站间再星间差分，可以消除接收机钟差。对于短基线（小于20公里），双差后电离层延迟和对流层延迟基本被消除，剩余主要未知量是坐标增量和整周模糊度。

双差载波相位观测方程写成矩阵形式为：

y = A * x + B * N + ε

其中y是双差载波相位观测值向量，x是基线向量增量（三维），N是双差模糊度向量（每个频率每颗卫星一个），A是几何矩阵，B是模糊度系数矩阵（一般为单位阵）。浮点解的目标是估计x和N的实数解及其协方差。

3.2 选参考星与构建双差观测值

双差需要先确定一颗参考星。参考星的选择标准是高度角最高且没有周跳的卫星。GPS和BDS分别选择各自的参考星，不能混用，因为GPS和BDS的星间差分不能跨系统合并，否则会产生系统间偏差项。

构建双差的代码逻辑：

def form_double_difference(base_obs, rover_obs, ref_sat, other_sats, freq): """ base_obs: 基准站观测值字典 rover_obs: 流动站观测值字典 ref_sat: 参考星PRN other_sats: 其他卫星PRN列表 freq: 频率标识，如'L1' 返回双差观测值列表 """ dd_list = [] base_ref = base_obs[ref_sat][freq]["carrier_phase"] rover_ref = rover_obs[ref_sat][freq]["carrier_phase"] sd_base_ref = base_ref + base_obs[ref_sat][freq]["pseudo_range"] * 0 sd_rover_ref = rover_ref + rover_obs[ref_sat][freq]["pseudo_range"] * 0 # 单差：流动站减基准站 sd_ref = rover_ref - base_ref for sat in other_sats: if sat not in base_obs or sat not in rover_obs: continue sd_sat = rover_obs[sat][freq]["carrier_phase"] - base_obs[sat][freq]["carrier_phase"] dd = sd_sat - sd_ref dd_list.append((sat, dd)) return dd_list

这段代码中，伪距乘以0是保留格式的占位写法，实际双差伪距需要另外计算。载波相位双差以米为单位参与计算，而模糊度是以周为单位，两者之间的关系是相位（米）等于波长乘以相位（周），所以后面的设计矩阵中要乘以波长。

3.3 最小二乘估计浮点解

不做卡尔曼滤波时，用最小二乘就能得到浮点解。观测方程线性化后，构建法方程：

N = A^T * P * A U = A^T * P * y x_hat = inv(N) * U

这里的P是权阵，通常取为高度角相关的对角阵。高度角越低，噪声越大，权重越小。加权策略采用简单的高度角正弦模型，权值为sin(elevation)^2。

import numpy as np def ls_solve(A, y, weights): """ 加权最小二乘求解浮点解 A: 设计矩阵 (m x n) y: 观测值向量 (m x 1) weights: 权阵对角元素 (m x 1) 返回 (参数估计, 协方差矩阵) """ P = np.diag(weights) N = A.T @ P @ A U = A.T @ P @ y try: cov = np.linalg.inv(N) x_hat = cov @ U return x_hat, cov except np.linalg.LinAlgError: print("法方程奇异，检查卫星数和几何构型") return None, None

设计矩阵A的构建按标准RTK流程进行。A的前三列是流动站到卫星的单位向量差，后几列对应各颗卫星的模糊度参数。坐标参数的初值可以用伪距单点定位结果或接收机输出的BESTPOS结果，通常精度在米级，迭代两三次即可收敛。

3.4 Lambda整数模糊度固定

浮点解得到模糊度的实数估计和协方差矩阵后，用Lambda方法搜索整数解。Lambda的核心思想是先对模糊度做Z变换降相关，再在变换域内搜索，最后还原到原域。

def lambda_search(amb_float, amb_cov, num_candidates=2): """ 简化版Lambda搜索 实际使用建议调用成熟的quad-rtk或RTKLIB的lambda实现 """ # Z变换降相关（此处示意，实际需实现整数Gauss变换） Z = np.eye(len(amb_float)) L = np.linalg.cholesky(amb_cov).T transformed_mean = np.linalg.inv(Z) @ amb_float # 在降相关域内做整数最小二乘搜索 # 搜索空间由卡方阈值确定 chi2 = 5.0 candidates = [] # ... 实际搜索逻辑较复杂，需要维护候选集并不断收缩搜索半径 return candidates

这里不贴完整Lambda实现，因为全部展开超出篇幅，但必须说清楚两点：第一，Lambda搜索要设置ratio阈值来验证固定解可靠性，通常ratio大于3才接受固定结果；第二，BDS参与解算时模糊度维数会明显增加，搜索耗时可能翻倍，降相关做得不好时尤其明显。

4. GPS/BDS双频融合解算的实现步骤与参数调优

4.1 系统间时间基准的统一

GPS时和BDS时存在14秒的整秒差（BDS时领先GPS时14秒），同时两个系统还有微小的闰秒偏差。NovAtel接收机输出的接收机时间通常统一为GPS时，但卫星钟差改正需要按各自系统的时间基准计算。处理方式是在读取星历时分别标记系统类型，计算卫星位置时用对应系统的时间参数，不强行把BDS时间转换到GPS时。

广播星历计算GPS和BDS卫星位置时，两者都采用开普勒轨道参数，但BDS的GEO卫星需要额外处理轨道倾角的小偏差。对课程设计来说，直接按ICD文档定义计算即可，注意BDS的轨道参数单位与GPS一致，但时间变量需要用BDS时elapsed seconds。

4.2 观测方程融合：同历元联合平差

GPS和BDS双频融合时，观测方程纵向拼接。假设一个历元有8颗GPS卫星和6颗BDS卫星，每颗卫星双频共2个观测值，总观测数约为28个（减去参考星后）。未知参数包括3个坐标增量，加GPS模糊度和BDS模糊度，每个频率独立模糊度，总计可能超过20个。

联合平差的好处是显著改善几何构型，尤其是在GPS卫星数不足5颗的场景下，BDS的加入能维持解算。实现时注意设计矩阵的拼接：

def build_design_matrix(gps_dd, bds_dd, base_pos, rover_pos_approx, sat_positions): """ gps_dd: GPS双差观测值列表 [(sat, phase_dd)] bds_dd: BDS双差观测值列表 [(sat, phase_dd)] """ rows = len(gps_dd) + len(bds_dd) num_amb = len(gps_dd) + len(bds_dd) A = np.zeros((rows, 3 + num_amb)) row_idx = 0 # GPS观测值的几何矩阵部分 for sat, _ in gps_dd: los = sat_positions[sat] - rover_pos_approx distance = np.linalg.norm(los) unit_vec = los / distance A[row_idx, 0:3] = -unit_vec # 流动站到卫星的单位向量取负 # 模糊度系数列 A[row_idx, 3 + row_idx] = 1.0 row_idx += 1 # BDS观测值类似，模糊度列偏移GPS模糊度个数 offset = len(gps_dd) for sat, _ in bds_dd: los = sat_positions[sat] - rover_pos_approx distance = np.linalg.norm(los) unit_vec = los / distance A[row_idx, 0:3] = -unit_vec A[row_idx, 3 + offset + (row_idx - offset)] = 1.0 row_idx += 1 return A

每个频率独立构建方程组，但坐标参数共用，也就是说L1频率和L2频率的观测方程里，前三列坐标参数是相同的，拼接时坐标列要重叠。另一种做法是直接把双频观测值按波长缩放后全部放入一个方程中解算。

4.3 关键参数设置与调优

RTK解算的核心参数表如下：

高度角截止的权衡是：太低会引入多径和大气残余误差，太高会减少可用卫星数导致模糊度不可固定。课程设计中建议手动调整并记录固定率变化，这是一个很好的报告分析点。

4.4 失败处理：卫星数不足与解算中断

BDS和GPS融合也不是万能的。常见失败场景是双系统参考星选择冲突、某一频率在基准站或流动站缺失、以及观测文件时间不同步。处理策略是：

def check_common_satellites(base_obs, rover_obs, freq): """ 检查同一历元基准站和流动站共视卫星 返回共视卫星列表 """ base_sats = set(base_obs.keys()) rover_sats = set(rover_obs.keys()) common = base_sats & rover_sats valid = [] for sat in common: if freq in base_obs[sat] and freq in rover_obs[sat]: valid.append(sat) return valid

如果有效卫星数少于4颗，本历元无法解算，直接跳过并记录。如果少于2颗，连双差都构建不起来。对连续解算的要求是：前一个历元的模糊度固定值可以作为当前历元的虚拟观测值，提高后续历元的固定成功率，这就是后续章节要讲的时间传递约束。

5. 基于实测数据验证解算结果：固定率分析与精度评估

5.1 数据准备与静基座验证方案

课程实习的常见数据来源有两种：一是NovAtel接收机采集的真实静态数据，通常基准站和流动站架设在已知点，基线长度从几米到几公里；二是公开数据集或实验室已有的存档数据。无论哪种，验证RTK解算的最直接方法是静态基线，因为基线真值已知，可以精确计算每个历元的定位误差。

标准操作是让接收机静态采集至少30分钟数据，流动站坐标设为已知点（或用PPK软件解算值作为参考），然后对比算法输出坐标序列的均值与真值。统计指标包括：固定率（固定解历元数占总历元数的比例）、RMS误差（三维方向）、模糊度固定后坐标稳定性。

5.2 可视化输出坐标时间序列与固定状态

解算完成后，把坐标输出为CSV文件，包含时间、X/Y/Z或经纬高、解算状态（固定/浮点/单点）、PDOP值。可视化时重点关注三件事：固定解是否连续、固定解坐标是否围绕真值波动、浮点解与固定解的跳变量级。

import matplotlib.pyplot as plt def plot_position_results(csv_file): """ 绘制坐标时间序列 固定解用蓝点，浮点解用红点 """ import pandas as pd df = pd.read_csv(csv_file) fix_mask = df['status'] == 'FIX' plt.figure(figsize=(12, 6)) plt.plot(df['gps_week_second'][fix_mask], df['east_error'][fix_mask], 'b.', markersize=2, label='固定解') plt.plot(df['gps_week_second'][~fix_mask], df['east_error'][~fix_mask], 'r.', markersize=2, label='浮点解') plt.axhline(y=0, color='k', linestyle='--', linewidth=0.8) plt.xlabel('时间 (秒)') plt.ylabel('东向误差 (米)') plt.legend() plt.grid(True, alpha=0.3) plt.show()

误差计算：把解算坐标与真值坐标都转到站心坐标系（ENU），才能直观看到水平精度和高程精度的差异。静态场景下固定解的水平RMS通常优于2厘米，高程RMS约3~4厘米，高程略差是因为卫星几何中天顶方向约束较弱。

5.3 固定率偏低时优先排查的三个方向

固定率如果低于70%，先检查观测值质量。用多普勒值和伪距变化量做对比，伪距的变化率应该接近多普勒乘以波长，偏差超过1米/秒说明伪距存在粗差。再检查基准站和流动站是否用同一颗参考星，如果某一系统参考星在其中一个站周跳频繁，就切换到另一颗。最后检查卫星位置计算是否正确，特别是BDS GEO卫星，轨道计算稍有小误差就会导致几厘米到几分米的系统性偏差，直接影响模糊度固定。

6. 进阶技巧：利用先验约束提升双频RTK的固定率

在常规解算跑通之后，可以加入历元间的模糊度连续约束。原理是：如果卫星在连续多个历元没有周跳，其整周模糊度应保持不变。将上一历元固定的模糊度作为当前历元的虚拟观测方程加入，能显著减少待估参数，增强法方程强度。

def add_ambiguity_constraint(A, y, weights, fixed_amb_map): """ 向观测方程追加模糊度约束 fixed_amb_map: 已固定模糊度字典 {param_index: (value, variance)} """ rows = len(fixed_amb_map) if rows == 0: return A, y, weights extra_cols = A.shape[1] A_new = np.vstack([A, np.zeros((rows, extra_cols))]) y_new = np.hstack([y, np.zeros(rows)]) w_new = np.hstack([weights, np.zeros(rows)]) for i, (idx, (val, var)) in enumerate(fixed_amb_map.items()): A_new[A.shape[0] + i, idx] = 1.0 y_new[A.shape[0] + i] = val w_new[A.shape[0] + i] = 1.0 / var return A_new, y_new, w_new

加入模糊度约束后有一个需要避免的坑：如果上一历元固定的模糊度本身就是错误的，约束会把后续历元都带偏。因此，只有当某个模糊度连续固定超过10个历元且ratio值持续超过3.0时，才将其加入约束集合。一旦该卫星被检测到周跳，立即从约束集合中移除。

另一个实用技巧是针对双频观测值做电离层残差的交叉验证。在短基线场景中，L1和L2的模糊度应该满足窄巷和宽巷的线性关系。利用L1模糊度 - L2模糊度 = 宽巷模糊度的约束，可以剔除部分模糊度搜索空间中的错误组合。这个方法在基线长度超过10公里时尤其有效，因为电离层残差会同时影响L1和L2，但组合约束可以抵消大部分影响。

对于BDS双频，B1/B2频率间隔比GPS L1/L2更大，电离层误差影响也更明显，所以在BDS双频中使用宽巷/窄巷约束的收益比GPS更高。实际操作中，先在宽巷域做一次固定，把宽巷模糊度取整固定后，再回代到原始观测方程求解窄巷模糊度，这样两步固定策略能把BDS模糊度的固定率提高10%到20%。

最后补充一个差异化对比的思路：分别用GPS单系统、BDS单系统和GPS+BDS联合解算同一组数据，对比三种模式的固定率和RMS。通常GPS+BDS联合解算在卫星数超过20颗时会明显优于单系统，但在卫星数较少时，联合解算因为几何构型差反而可能无法固定。用这个实验来写课设报告，既展示了算法实现能力，也体现了对系统间差异的深入理解。

本文还有配套的精品资源，点击获取