简介:这份资源面向惯性导航、组合导航方向的学习者与工程人员,提供一套基于扩展卡尔曼滤波(EKF)的INS组合导航MATLAB实现代码,可用于理解姿态、速度与位置估计的完整流程,并作为算法验证与课程设计的参考基础。压缩包共19个文件,约21KB,以cpp与h源码为主,辅以Makefile构建脚本、README说明及gitignore配置,涵盖陀螺仪、加速度计、磁力计、GPS等传感器模块与EKF核心、GPS滤波等实现,结构清晰便于按模块阅读。目前已有468人学习下载,适合希望从代码层面掌握多传感器数据融合与滤波更新的读者,可据此复现仿真、对比估计结果,并在此基础上改进滤波算法或扩展组合导航方案。
1. INS_EKF 组合导航:从一份开源代码到跑通自己的数据
拿到INS_EKF-master.zip这类组合导航代码包的人,多半卡在同一个地方:代码能编译,Demo 能跑,但换成自己的 IMU 和 GNSS 数据就发散。组合导航的核心不是把卡尔曼滤波公式抄对,而是搞清楚状态量怎么定义、噪声参数怎么标、观测在什么条件下可信。INS_EKF 这套代码用扩展卡尔曼滤波把惯性导航的短时高精度和卫星导航的长期稳定性捏在一起,解决的是单一传感器在隧道、城市峡谷、高动态场景下要么漂移要么失锁的问题。它适合做机器人导航、自动驾驶定位、无人机飞控的工程师,也适合刚接触惯性导航、想找一个能改能调的最小系统的学生。下面按“先立住理论、再动手复现、最后避坑”的顺序拆开讲,每一步都落到可执行的命令和参数上。
2. INS_EKF 的状态方程与观测方程:15 维误差状态怎么摆
2.1 为什么用误差状态而不是全状态
组合导航里直接对位置、速度、姿态做 EKF 不是不行,但姿态用四元数或旋转矩阵时,线性化容易出问题,协方差矩阵也难保持正定。工程上更常见的是误差状态卡尔曼滤波,把标称状态和误差状态分开:标称状态用惯性递推走,误差状态用 EKF 估计,最后把误差反馈回去修正标称值。INS_EKF 这类代码通常采用 15 维误差状态,顺序是位置误差、速度误差、姿态误差、加速度计零偏、陀螺仪零偏。这样做的直接好处是姿态误差在小角度下近似为三维旋转向量,线性化误差可控,协方差更新时数值稳定性明显更好。
状态量定义清楚了,接下来要确认坐标系。常见做法是导航系取东北天,机体系取右前上,IMU 测量值在机体系下。如果代码里导航系是北东地,而你的数据是东北天,那姿态解算会直接翻车。我一般会先找到代码里定义导航系的那几行,确认重力向量是[0, 0, -9.81]还是[0, 0, 9.81],这一步花五分钟,能省后面五小时。
2.2 状态转移矩阵的离散化与参数含义
连续时间误差状态方程写成 $\dot{\delta x} = F \delta x + G w$,代码里真正用的是离散形式 $\delta x_{k+1} = \Phi_k \delta x_k + w_k$。INS_EKF 里通常用一阶近似 $\Phi = I + F \Delta t$,在 IMU 频率 100Hz 到 200Hz 下足够。下面这段伪代码展示状态转移矩阵的核心构造逻辑,实际代码里会拆成多个函数。
# 离散化状态转移矩阵构造(示意,对应 15 维误差状态) # 状态顺序: [dp(3), dv(3), dtheta(3), ba(3), bg(3)] def build_F(dt, R, fn, vn): F = np.eye(15) # 位置误差对速度误差的偏导 F[0:3, 3:6] = np.eye(3) * dt # 速度误差对姿态误差的偏导: -R * skew(fn) * dt F[3:6, 6:9] = -R @ skew(fn) * dt # 速度误差对加速度计零偏的偏导: -R * dt F[3:6, 9:12] = -R * dt # 姿态误差对陀螺仪零偏的偏导: -I * dt F[6:9, 12:15] = -np.eye(3) * dt return F这段代码里R是机体系到导航系的旋转矩阵,fn是比力在导航系下的投影,skew是反对称矩阵。参数dt必须和 IMU 实际输出周期一致,如果 IMU 标称 100Hz 但实际有抖动,建议用时间戳差分而不是固定 0.01。F[3:6, 6:9]这一项最容易写错符号,符号错了姿态和速度会耦合发散。F[3:6, 9:12]和F[6:9, 12:15]分别对应零偏的可观测性,如果这两项置零,零偏估计会完全依赖观测更新,收敛极慢。
2.3 观测方程:GNSS 位置速度怎么进滤波器
观测方程 $z = H \delta x + v$ 里,GNSS 通常提供位置和速度。如果只用位置,H 矩阵就是 $[I_3, 0_{3\times12}]$;如果位置速度都用,H 是 $[I_3, 0; 0, I_3, 0]$。量测噪声矩阵 R 一般取对角阵,对角线是 GNSS 位置和速度的方差。这里有个实操细节:很多开源代码把 R 写成固定值,但 GNSS 在城市环境里方差会剧烈变化。我一般会根据卫星数、HDOP 或者 RTK 的 fix 类型动态调 R,固定值在开阔地能用,进城区就等着看位置跳变。
# GNSS 位置+速度观测更新 def update_gnss(P, dx, H, R, z, z_pred): # 新息 y = z - z_pred # 新息协方差 S = H @ P @ H.T + R # 卡尔曼增益 K = P @ H.T @ np.linalg.inv(S) # 状态修正 dx = dx + K @ y # 协方差修正 (Joseph 形式,数值更稳) I_KH = np.eye(15) - K @ H P = I_KH @ P @ I_KH.T + K @ R @ K.T return dx, Pz_pred是用标称状态算出来的预测观测,位置预测就是标称位置,速度预测就是标称速度。S矩阵求逆前建议加一个极小对角项防止奇异。Joseph 形式比简单形式多算一次矩阵乘法,但能保证 P 对称正定,长时间跑下来差别很明显。R 的对角线取值:消费级 GNSS 位置方差可以给 1 到 5 平方米,速度方差给 0.1 到 0.5;RTK 固定解可以降到 0.01 平方米量级。这些值不是拍脑袋,最好用静态数据算一下实际残差。
3. 把 INS_EKF 跑起来:数据接口、时间同步与最小验证
3.1 数据准备:IMU 和 GNSS 的格式对齐
拿到代码后第一件事不是改算法,是让数据能喂进去。INS_EKF 这类工程通常期望 IMU 数据包含时间戳、三轴角速度、三轴加速度,GNSS 数据包含时间戳、经纬高、东北天速度。常见做法是把数据转成 CSV 或直接写一个数据读取类。下面是一个最小数据加载示例,假设 IMU 是 200Hz,GNSS 是 10Hz。
import pandas as pd import numpy as np # IMU: timestamp, gx, gy, gz (rad/s), ax, ay, az (m/s^2) imu = pd.read_csv('imu.csv') imu_ts = imu['timestamp'].values gyro = imu[['gx','gy','gz']].values accel = imu[['ax','ay','az']].values # GNSS: timestamp, lat, lon, alt, ve, vn, vu gnss = pd.read_csv('gnss.csv') gnss_ts = gnss['timestamp'].values gnss_pos = gnss[['lat','lon','alt']].values gnss_vel = gnss[['ve','vn','vu']].values # 检查时间戳单位,常见坑:秒 vs 毫秒 vs 微秒 assert imu_ts[1] - imu_ts[0] < 1.0, "IMU 时间戳单位可能不是秒"时间戳单位是第一个大坑。有的数据集用秒,有的用毫秒,有的用 GPS 周内秒。单位搞错,dt 就错,状态转移矩阵全错。我一般会打印前几个时间戳差值,如果 IMU 标称 200Hz,差值应该在 0.005 秒左右;如果是 5 或者 5000,那就得换算。GNSS 的经纬高要转成局部东北天坐标,常见做法是选第一个 GNSS 点作为原点,用墨卡托或简单球面模型转平面。
3.2 时间同步:IMU 和 GNSS 怎么对齐
IMU 和 GNSS 频率不同,EKF 更新时要么在 GNSS 时刻做观测更新,要么把 GNSS 插值到 IMU 时刻。常见做法是:以 IMU 为主循环,每来一帧 IMU 做一次状态预测,当 IMU 时间戳超过当前 GNSS 时间戳时,用最近一帧 GNSS 做观测更新。这里要注意 GNSS 本身有延迟,有的接收机输出的是天线收到信号的时间,不是解算完成时间。如果延迟没补偿,高速场景下位置会滞后。
gnss_idx = 0 for i in range(len(imu_ts)): dt = imu_ts[i] - imu_ts[i-1] if i > 0 else 0.005 # 状态预测 dx, P = predict(dx, P, gyro[i], accel[i], dt) # 判断是否有 GNSS 更新 while gnss_idx < len(gnss_ts) and gnss_ts[gnss_idx] <= imu_ts[i]: dx, P = update_gnss(P, dx, H, R, gnss_pos[gnss_idx], gnss_vel[gnss_idx]) gnss_idx += 1 # 误差反馈到标称状态 nom_state = feedback(nom_state, dx) dx[:] = 0 # 反馈后误差状态清零feedback函数把位置、速度、姿态、零偏的误差加到标称状态上,姿态误差反馈时要注意用旋转矩阵更新而不是直接加。反馈后误差状态清零是标准做法,不清零的话误差会重复累加。while循环处理的是同一时刻有多帧 GNSS 的情况,实际中很少见,但写了不亏。
3.3 最小验证:静态和直线跑车
代码跑通后不要直接上复杂场景,先做两个最小验证。静态测试:IMU 和 GNSS 都静止,跑十分钟,看位置漂移和姿态收敛。如果位置在几分钟内漂到几十米,说明零偏估计没收敛或者 R 给得太小。直线测试:沿一条已知直线走,看轨迹和 GNSS 轨迹的偏差。这两个测试能暴露大部分参数问题。
| 验证项 | 数据时长 | 关注指标 | 合格参考 |
|---|---|---|---|
| 静态 | 10 分钟 | 位置漂移 | 小于 1 米 |
| 静态 | 10 分钟 | 横滚俯仰收敛 | 收敛到 0.1 度内 |
| 直线 | 2 分钟 | 轨迹偏差 | 小于 2 米 |
| 直线 | 2 分钟 | 速度误差 | 小于 0.2 m/s |
静态测试时 GNSS 位置本身有噪声,所以位置不会完全不动,但应该在 GNSS 噪声范围内。如果位置单调漂移,检查加速度计零偏是否可观测、R 是否过大导致 GNSS 权重太低。直线测试时如果轨迹弯曲,检查航向初始对准和陀螺零偏。
4. INS_EKF 调参避坑:从发散到收敛的 5 个关键点
4.1 现象:滤波器一开始就发散,位置飞到天上
原因通常是初始协方差 P0 给得太小,或者初始姿态误差太大。EKF 是线性化滤波,初始误差超出线性范围就会发散。解决方法是把 P0 的姿态部分放大到 10 度到 30 度的方差,位置部分给 10 米方差,速度给 1 米每秒方差。零偏部分给 0.1 到 1 的方差。先让滤波器“敢修”,收敛后再收紧。
4.2 现象:静态下姿态缓慢旋转
原因是陀螺零偏没有估计,或者零偏的状态转移噪声 Q 太小。陀螺零偏在静态下是可观测的,但如果 Q 给太小,滤波器不相信零偏会变,收敛极慢。我一般把陀螺零偏的过程噪声给到 1e-5 到 1e-4 量级,加速度计零偏给 1e-4 到 1e-3。具体值取决于 IMU 等级,消费级 IMU 零偏稳定性差,Q 要大一些。
4.3 现象:GNSS 更新后位置跳变
原因是 R 给得太小,滤波器过度相信 GNSS。GNSS 在遮挡或多径时位置误差可能到十几米,如果 R 还是 1 平方米,更新后位置会猛跳。解决方法是根据 HDOP 或卫星数动态调 R,HDOP 大于 2 时把位置方差放大到 10 以上。另外可以做新息卡方检验,新息超过阈值时跳过这次更新。
4.4 现象:长时间跑下来协方差矩阵不对称
原因是数值误差累积,简单形式的协方差更新不保证对称正定。解决方法是改用 Joseph 形式,并且每次更新后强制对称P = (P + P.T) / 2。如果还不行,检查 F 矩阵和 H 矩阵的维度是否匹配,维度错了有时不报错但结果完全错。
4.5 现象:换一组数据就完全不能用
原因是参数过拟合到某一组数据。组合导航的参数应该基于传感器手册和实际噪声标定,而不是调参调到某组数据好看。我一般会留一组独立数据做验证,调参只在训练数据上做。如果换数据后误差翻倍,先检查时间戳单位、坐标系定义、单位是否一致,这三个地方最容易出问题。
5. 进阶技巧:用新息序列判断 INS_EKF 有没有真正收敛
5.1 新息白化检验:比看轨迹更早发现问题
轨迹好看不代表滤波器健康。新息序列 $y_k = z_k - H x_k^-$ 在滤波器收敛时应该是零均值白噪声。如果新息有偏或者自相关,说明模型和实际不匹配。实操中我会记录每次 GNSS 更新的新息,跑完后画出来看。位置新息应该在正负两倍标准差内波动,如果持续偏一边,说明零偏或安装角没估计对。
# 新息记录与白化检验 innovations = [] for each gnss update: y = z - z_pred innovations.append(y) # 归一化新息平方 (NIS) nis = y.T @ np.linalg.inv(S) @ y if nis > chi2_threshold: # 例如 3 自由度 99% 置信度约 11.34 skip_update = True # 跑完后检查新息均值 innovations = np.array(innovations) print("位置新息均值:", innovations[:, :3].mean(axis=0)) print("位置新息标准差:", innovations[:, :3].std(axis=0))S是新息协方差,chi2_threshold根据观测维度查卡方分布表。如果 NIS 频繁超阈值,说明 R 给太小或者观测有异常值。这个检验比看轨迹灵敏得多,轨迹在滤波器快发散时可能还看着正常。
5.2 零偏收敛曲线:判断可观测性的直接手段
加速度计和陀螺零偏的估计曲线应该先快速收敛然后稳定。如果零偏一直漂,说明该轴不可观测。静态下水平加速度计零偏和水平陀螺零偏可观测,垂直陀螺零偏不可观测。如果发现某个轴零偏不收敛,不要强行调 Q,先检查运动激励是否足够。车载场景下,转弯和加减速能激励所有轴,直线匀速跑则很多状态不可观测。
5.3 我自己的习惯
每次换传感器,我先跑静态两小时,看零偏收敛和 Allan 方差,把 Q 和 R 定下来,再跑动态。动态数据先跑直线,再跑转弯,最后跑复杂场景。调参记录写在本子上,每次只改一个参数,改完跑同一组数据对比。这套流程看起来慢,但比反复推倒重来快得多。组合导航没有银弹,参数背后是对传感器和场景的理解。希望帮到你。
本文还有配套的精品资源,点击获取