GNSS/SINS组合导航系统由于其优异的性能，已成为导航应用的主要选择，而卡尔曼滤波器作为组合导航系统常用的信息融合方法，其滤波性能严重依赖于准确的系统噪声及测量噪声先验统计特性，利用错误的噪声统计特性将会导致滤波性能下降甚至发散^[1]。然而，在实际应用环境下，由于内部及外部原因，惯性导航传感器存在较大的动态误差，标度因子导致的误差在滤波过程中很难得到补偿，进而导致系统噪声统计特性不准确。当GNSS接收机处于电磁干扰环境下时，测量噪声的统计模型无法准确进行描述^[2]。

自适应滤波方法是解决上述问题的有效途径^[3]。Sage Husa自适应滤波是一种研究较早的方法^[4-6]，然而在实际应用中，该滤波存在无法同时估计系统噪声及测量噪声统计特性、用于组合导航系统时经常存在滤波发散现象等缺陷。交互式多模型滤波器利用多个滤波模型来解决参数变化系统的状态估计问题，并被应用于组合导航测量多模型的切换^[7-9]，但该方法存在运行时间较长、多模型独立参数构成的模型集需涵盖参数范围等缺陷。近年来，基于变分贝叶斯的自适应滤波方法得到长足发展，并被成功应用于组合导航系统中^[10]，以克服Sage Husa自适应滤波在应用中出现的发散现象。但该方法目前只能在线实时估计测量噪声方差的变化，无法克服不准确的系统噪声方差统计特性对滤波精度的影响^[11]。林雪原等^[12]将变分贝叶斯估计方法与抗野值估计方法相结合，以解决测量噪声方差变化、测量值存在硬故障等问题，但对软故障的处理效果不佳。根据卡尔曼滤波方程可以发现，系统噪声方差矩阵会直接影响状态估计的一步预测误差协方差矩阵，且相对于系统噪声方差矩阵，状态估计的一步预测误差协方差矩阵对状态估计的影响更直接，更容易被估计。

为此，针对不准确的系统噪声及测量噪声先验统计，本文基于变分贝叶斯理论，在对系统状态估计的一步预测协方差矩阵及测量噪声的联合后验分布进行近似的基础上，提出GNSS/INS组合导航系统的改进变分贝叶斯自适应滤波方法，以得到复杂环境下组合导航系统的准确滤波结果。

1 组合导航滤波模型

GNSS/INS组合导航系统松组合滤波模型为：

$\left\{\begin{array}{l} \boldsymbol{X}_{k}=\boldsymbol{\varPhi}_{k / k-1} \boldsymbol{X}_{k-1}+\boldsymbol{\varGamma}_{k-1} \boldsymbol{W}_{k-1} \\ \boldsymbol{Z}_{k}=\boldsymbol{H}_{k} \boldsymbol{X}_{k}+\boldsymbol{V}_{k} \end{array}\right.$

(1)

式中，X_k为n×1维状态向量；Φ_k/k－1为k－1时刻到k时刻的状态转移矩阵；Γ_k－1为系统噪声分配矩阵；W_k－1为系统噪声向量，且W_k~N(0, Q_k)；Z_k为m×1维测量向量；H_k为测量矩阵；V_k为测量噪声向量，且V_k~N(0, R_k)。

选取“东北天”导航坐标系，记为n系；“右前上”载体坐标系，记为b系；地心地固系，记为e系；采用15维误差状态向量：

$\boldsymbol{X}=\left[\begin{array}{l} \left(\boldsymbol{\varphi}^n\right)^{\mathrm{T}}\quad\left(\delta \boldsymbol{v}^n\right)^{\mathrm{T}}\quad\left(\delta \boldsymbol{p}^e\right)^{\mathrm{T}}\quad\left(\boldsymbol{\varepsilon}^b\right)^{\mathrm{T}}\quad\left(\boldsymbol{\nabla}^b\right)^{\mathrm{T}} \end{array}\right]^{\mathrm{T}}$

(2)

式中，φⁿ=[φ_E  φ_N  φ_U]^T为n系下的姿态误差角，δvⁿ=[δv_Eⁿ  δv_Nⁿ  δv_Uⁿ]^T为东、北、天3个方向的速度误差，δp^e=[δL  δλ  δh]^T为e系下的纬度、经度和高度3个方向的位置误差，ε^b=[ε_x^b  ε_y^b  ε_z^b]^T为载体系下陀螺仪三轴一阶马尔科夫误差，▽^b=[▽_x^b  ▽_y^b  ▽_z^b]^T为载体系下加速度计三轴一阶马尔科夫误差。

观测向量Z_k取为GNSS与SINS提供的三维位置和三维速度的对应差值。

基于式(1)和式(2)，卡尔曼滤波方程表示为：

$\left\{\begin{array}{l} \hat{\boldsymbol{X}}_{k \mid k-1}=\boldsymbol{\varPhi}_{k, k-1} \hat{\boldsymbol{X}}_{k-1 \mid k-1} \\ \boldsymbol{P}_{k \mid k-1}=\boldsymbol{\varPhi}_{k, k-1} \boldsymbol{P}_{k-1 \mid k-1} \boldsymbol{\varPhi}_{k, k-1}^{\mathrm{T}}+\boldsymbol{\varGamma}_{k-1} \boldsymbol{Q}_{k-1} \boldsymbol{\varGamma}_{k-1}^{\mathrm{T}} \\ \boldsymbol{K}_{k}=\boldsymbol{P}_{k \mid k-1} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\boldsymbol{H}_{k} \boldsymbol{P}_{k \mid k-1} \boldsymbol{H}_{k}^{\mathrm{T}}+\boldsymbol{R}_{k}\right)^{-1} \\ \boldsymbol{P}_{k \mid k}=\left[\boldsymbol{I}-\boldsymbol{K}_{k} \boldsymbol{H}_{k}\right] \boldsymbol{P}_{k \mid k-1} \\ \hat{\boldsymbol{X}}_{k \mid k}=\hat{\boldsymbol{X}}_{k \mid k-1}+\boldsymbol{K}_{k}\left(\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \hat{\boldsymbol{X}}_{k \mid k-1}\right) \end{array}\right.$

(3)

式中，$\hat{\boldsymbol{X}}_{k-1}$为k－1时刻的状态估计值，$\hat{\boldsymbol{X}}_{k \mid k-1}$、P_k|k－1、K_k和P_k分别为k时刻的状态预测值、预测误差协方差矩阵、滤波器增益矩阵和状态协方差矩阵估计值，Q_k－1为系统噪声方差矩阵，R_k为已知的测量噪声方差矩阵。

2 MVBAKF算法 2.1 先验分布选择及时间更新

根据式(3)，Q_k－1的不准确性将直接影响P_k|k－1。所以，当Q_k－1和R_k不准确时，在估计X_k的同时需估计P_k|k－1和R_k，以保证滤波精度。

将P_k|k－1和R_k的先验分布p(P_k|k－1|Z_1:k－1)和p(R_k|Z_1:k－1)选为逆Wishart概率密度函数(probability density function, PDF)^[13]：

$p\left(\boldsymbol{P}_{k \mid k-1} \mid \boldsymbol{Z}_{1: k-1}\right)=\operatorname{IW}\left(\boldsymbol{P}_{k \mid k-1} ; {t_{\hat{k} \mid k-1}}, \hat{\boldsymbol{T}}_{k \mid k-1}\right)$

(4)

$p\left(\boldsymbol{R}_{k} \mid \boldsymbol{Z}_{1: k-1}\right)=\operatorname{IW}\left(\boldsymbol{R}_{k} ; \hat{u}_{k \mid k-1}, \hat{\boldsymbol{U}}_{k \mid k-1}\right)$

(5)

式中，IW(B; λ, Ψ)代表自由度参数为λ和逆尺度矩阵为Ψ的逆Wishart PDF；B和Ψ均为对称正定矩阵，并假设其维数为d×d。

IW(B; λ, Ψ)可以表示为：

$\begin{gathered} \operatorname{IW}(\boldsymbol{B} ; \lambda, \boldsymbol{\varPsi})=\\ \frac{|\boldsymbol{\varPsi}|^{\lambda / 2}|\boldsymbol{B}|^{-(\lambda+d+1) / 2} \exp \left\{-0.5 \operatorname{tr}\left(\boldsymbol{\varPsi} \cdot \boldsymbol{B}^{-1}\right)\right\}}{2^{d \lambda / 2} T_{d}(\lambda / 2)} \end{gathered}$

(6)

式中，|·|和tr(·)分别为行列式运算和求迹运算，T_d(·)是一个d元伽马函数。

如果B~IW(B; λ, Ψ)，当λ>d+1时有：

$E\left[\boldsymbol{B}^{-1}\right]=(\lambda-d-1) \boldsymbol{\varPsi}^{-1}$

(7)

为获取P_k|k－1的先验信息，利用式(7)得到P_k|k－1的均值为：

$\begin{gather*} \frac{\hat{\boldsymbol{T}}_{k \mid k-1}}{\hat{t}_{k \mid k-1}-n-1}=\widetilde{\boldsymbol{P}}_{k \mid k-1}= \\ \boldsymbol{\varPhi}_{k, k-1} \boldsymbol{P}_{k-1 \mid k-1} \boldsymbol{\varPhi}_{k, k-1}^{\mathrm{T}}+\boldsymbol{\varGamma}_{k-1} \widetilde{\boldsymbol{Q}}_{k-1} \boldsymbol{\varGamma}_{k-1}^{\mathrm{T}} \end{gather*}$

(8)

式中，$\widetilde{\boldsymbol{Q}}_{k-1}$为Q_k－1的近似值；并设定：

$\hat{t}_{k \mid k-1}=n+\tau+1$

(9)

式中，τ≥0是一个调整因子。将式(9)代入式(8)，得：

$\hat{\boldsymbol{T}}_{k \mid k-1}=\tau \widetilde{\boldsymbol{P}}_{k \mid k-1}$

(10)

同理，为获取R_k的先验信息，R₀的均值可表示为：

$\widetilde{\boldsymbol{R}}_{0}=\frac{\hat{\boldsymbol{U}}_{0 \mid 0}}{\hat{u}_{0 \mid 0}-m-1}$

(11)

式中，$\hat{\boldsymbol{U}}$和$\hat{u}$可通过一个遗忘因子ρ∈(0, 1)进行预测：

$\hat{u}_{k \mid k-1}=\rho\left(\hat{u}_{k-1 \mid k-1}-m-1\right)+m+1$

(12)

$\hat{\boldsymbol{U}}_{k \mid k-1}=\rho \hat{\boldsymbol{U}}_{k-1 \mid k-1}$

(13)

2.2 基于变分的测量更新

为同时估计X_k、P_k|k－1和R_k，需计算联合后验PDFp(X_k, P_k|k－1, R_k|Z_1:k)，利用变分贝叶斯方法可得^[14-15]：

$p\left(\boldsymbol{X}_{k}, \boldsymbol{P}_{k \mid k-1}, \boldsymbol{R}_{k} \mid \boldsymbol{Z}_{1: k}\right) \approx q\left(\boldsymbol{X}_{k}\right) q\left(\boldsymbol{P}_{k \mid k-1}\right) q\left(\boldsymbol{R}_{k}\right)$

(14)

式中，q(·)代表p(·)的近似后验PDF。

通过最小化近似后验PDF与真实后验PDF之间的KLD距离，得到：

$\log q(\theta) =E_{\varXi^{(-\theta)}}\left[\log p\left(\boldsymbol{\varXi}, \boldsymbol{Z}_{1: k}\right)\right]+c_{\theta}$

(15)

$\boldsymbol{\varXi} \triangleq\left\{\boldsymbol{X}_{k}, \boldsymbol{P}_{k \mid k-1}, \boldsymbol{R}_{k}\right\}$

(16)

式中，联合PDF logp(Ξ, Z_1:k)可以表示为：

$\begin{gather*} p\left(\boldsymbol{\varXi}, \boldsymbol{Z}_{1: k}\right)=\\ p\left(\boldsymbol{Z}_{k} \mid \boldsymbol{X}_{k}, \boldsymbol{R}_{k}\right) p\left(\boldsymbol{X}_{k} \mid \boldsymbol{Z}_{1: k-1}, \boldsymbol{P}_{k \mid k-1}\right) \times \\ p\left(\boldsymbol{P}_{k \mid k-1} \mid \boldsymbol{Z}_{1: k-1}\right) p\left(\boldsymbol{R}_{k} \mid \boldsymbol{Z}_{1: k-1}\right) p\left(\boldsymbol{Z}_{1: k-1}\right)= \\ N\left(\boldsymbol{Z}_{k} ; \boldsymbol{H}_{k} \boldsymbol{X}_{k}, \boldsymbol{R}_{k}\right) N\left(\boldsymbol{X}_{k} ; \hat{\boldsymbol{X}}_{k \mid k-1}, \boldsymbol{P}_{k \mid k-1}\right) \times \\ \operatorname{IW}\left(\boldsymbol{P}_{k \mid k-1} ; \hat{t}_{k \mid k-1}, \hat{\boldsymbol{T}}_{k \mid k-1}\right) \\ \operatorname{IW}\left(\boldsymbol{R}_{k} ; \hat{u}_{k \mid k-1}, \hat{\boldsymbol{U}}_{k \mid k-1}\right) p\left(\boldsymbol{Z}_{1: k-1}\right) \end{gather*}$

(17)

式中，log(·)代表对数函数，E(·)代表期望运算，θ是Ξ中的任一元素，Ξ^(－θ)是除了θ以外所有元素的集合，c_θ代表变量θ的常数，N(·;μ, Σ)代表均值为μ、协方差矩阵为Σ的高斯PDF。

基于式(17)，logp(Ξ, Z_1:k)可表示为：

$\begin{gathered} \log p\left(\boldsymbol{\varXi}, \boldsymbol{Z}_{1: k}\right)=-\\ 0.5\left(\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \boldsymbol{X}_{k}\right){ }^{\mathrm{T}} \boldsymbol{R}_{k}^{-1}\left(\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \boldsymbol{X}_{k}\right)-\\ 0.5\left(m+\hat{\mu}_{k \mid k-1}+2\right) \log \left|\boldsymbol{R}_{k}\right|-0.5\\ \operatorname{tr}\left(\hat{\boldsymbol{U}}_{k \mid k 1} \boldsymbol{R}_{k}^{-1}\right)-0.5\left(n+\hat{t}_{k \mid k-1}+2\right)\\ \log \left|\boldsymbol{P}_{k \mid k-1}\right|-0.5 \operatorname{tr}\left(\hat{\boldsymbol{T}}_{k \mid k-1} \boldsymbol{P}_{k \mid k-1}^{-1}\right)-\\ 0.5\left(\boldsymbol{X}_{k}-\hat{\boldsymbol{X}}_{k \mid k-1}\right)^{\mathrm{T}} \boldsymbol{P}_{k \mid k-1}^{-1}\left(\boldsymbol{X}_{k}-\hat{\boldsymbol{X}}_{k \mid k-1}\right)+c_{\varXi} \end{gathered}$

(18)

令θ=P_k|k－1，且将式(18)代入式(15)，得：

$\begin{gather*} \log q^{(i+1)}\left(\boldsymbol{P}_{k \mid k-1}\right)=- \\ 0.5\left(n+\hat{t}_{k \mid k-1}+2\right) \log \mid \boldsymbol{P}_{k \mid k-1}- \\ 0.5 \operatorname{tr}\left(\left(\hat{\boldsymbol{T}}_{k \mid k-1}+\boldsymbol{A}_{k}^{(i)}\right) \boldsymbol{P}_{k \mid k-1}^{-1}\right)+c_{p} \end{gather*}$

(19)

式中，q⁽ⁱ⁺¹⁾(·)为第i+1次迭代的近似PDF，且A_k⁽ⁱ⁾可表示为：

$\boldsymbol{A}_{k}^{(i)}=\boldsymbol{P}_{k \mid k}^{(i)}+\left(\boldsymbol{X}_{k \mid k}^{(i)}-\boldsymbol{X}_{k \mid k-1}\right)\left(\boldsymbol{X}_{k \mid k}^{(i)}-\boldsymbol{X}_{k \mid k-1}\right)^{\mathrm{T}}$

(20)

基于式(19)，可以更新为逆Wishart PDF：

$q^{(i+1)}\left(\boldsymbol{P}_{k \mid k-1}\right)=\operatorname{IW}\left(\boldsymbol{P}_{k \mid k-1} ; \hat{t}_{k}^{(i+1)}, \hat{\boldsymbol{T}}_{k}^{(i+1)}\right)$

(21)

其中，

$\hat{t}_{k}^{(i+1)}=\hat{t}_{k \mid k-1}+1$

(22)

$\hat{\boldsymbol{T}}_{k}^{(i+1)}=\hat{\boldsymbol{T}}_{k \mid k-1}+\hat{\boldsymbol{A}}_{k}^{(i)}$

(23)

基于式(21)，得到修正的P_k|k－1⁽ⁱ⁺¹⁾为：

$\boldsymbol{P}_{k \mid k-1}^{(i+1)}=\frac{\hat{\boldsymbol{T}}_{k}^{(i+1)}}{\hat{t}_{k}^{(i+1)}-n-1}$

(24)

令θ=R_k，且将式(18)代入式(15)，得：

$\begin{gather*} \log q^{(i+1)}\left(\boldsymbol{R}_{k}\right)=-0.5\left(m+\hat{\mu}_{k \mid k-1}+2\right) \log \left|\boldsymbol{R}_{k}\right|- \\ 0.5 \operatorname{tr}\left(\left(\hat{\boldsymbol{U}}_{k \mid k 1}+\boldsymbol{B}_{k}^{(i)}\right) \boldsymbol{R}_{k}^{-1}\right)+c_{\boldsymbol{R}} \end{gather*}$

(25)

式中，B_k⁽ⁱ⁾可以表示为：

$\begin{gather*} \boldsymbol{B}_{k}^{(i)}=\left(\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \hat{\boldsymbol{X}}_{k \mid k}^{(i)}\right)\left(\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \hat{\boldsymbol{X}}_{k \mid k}^{(i)}\right)^{\mathrm{T}}+ \\ \boldsymbol{H}_{k} \boldsymbol{P}_{k \mid k}^{(i)} \boldsymbol{H}_{k}^{\mathrm{T}} \end{gather*}$

(26)

基于式(25)，q(R_k)更新为逆Wishart PDF：

$q^{(i+1)}\left(\boldsymbol{R}_{k}\right)=\operatorname{IW}\left(\boldsymbol{R}_{k} ; \hat{u}_{k}^{(i+1)}, \hat{\boldsymbol{U}}_{k}^{(i+1)}\right)$

(27)

其中，

$\hat{u}_{k}^{(i+1)}=\hat{u}_{k \mid k-1}+1$

(28)

$\hat{\boldsymbol{U}}_{k}^{(i+1)}=\hat{\boldsymbol{U}}_{k \mid k 1}+\boldsymbol{B}_{k}^{(i)}$

(29)

则修正后的R_k⁽ⁱ⁺¹⁾可以表示为：

$\boldsymbol{R}_{k}^{(i+1)}=\frac{\hat{\boldsymbol{U}}_{k}^{(i+1)}}{\hat{u}_{k}^{(i+1)}-m-1}$

(30)

同理，令θ=X_k，且将式(18)代入式(15)，q(X_k)更新为高斯PDF，得：

$q^{(i+1)}\left(\boldsymbol{X}_{k}\right)=N\left(\boldsymbol{X}_{k} ; \hat{\boldsymbol{X}}_{k \mid k}^{(i+1)}, \boldsymbol{P}_{k \mid k}^{(i+1)}\right)$

(31)

式中，$\hat{\boldsymbol{X}}_{k \mid k}^{(i+1)}$和P_k|k⁽ⁱ⁺¹⁾代表X_k的均值向量和协方差矩阵，可以表示为：

$\boldsymbol{K}_{k}^{(i+1)}=\boldsymbol{P}_{k \mid k-1}^{(i+1)} \boldsymbol{H}_{k}^{\mathrm{T}}\left(\boldsymbol{H}_{k} \boldsymbol{P}_{k \mid k-1}^{(i+1)} \boldsymbol{H}_{k}^{\mathrm{T}}+\boldsymbol{R}_{k}^{(i+1)}\right)^{-1}$

(32)

$\boldsymbol{X}_{k \mid k}^{(i+1)}=\boldsymbol{X}_{k \mid k-1}+\boldsymbol{K}_{k}^{(i+1)}\left(\boldsymbol{Z}_{k}-\boldsymbol{H}_{k} \boldsymbol{X}_{k \mid k-1}\right)$

(33)

$\boldsymbol{P}_{k \mid k}^{(i+1)}=\left(I-\boldsymbol{K}_{k}^{(i+1)} \boldsymbol{H}_{k}\right) \boldsymbol{P}_{k \mid k-1}^{(i+1)}$

(34)

经过N次定点迭代运算后，后验PDF q(X_k)、q(P_k|k－1)和q(R_k)的变分近似可以表示为：

$\begin{gather*} q\left(\boldsymbol{X}_{k}\right) \approx q^{(N)}\left(\boldsymbol{X}_{k}\right)=N\left(\boldsymbol{X}_{k} ; \hat{\boldsymbol{X}}_{k \mid k}^{(N)}, \boldsymbol{P}_{k \mid k}^{(N)}\right)= \\ N\left(\boldsymbol{X}_{k} ; \hat{\boldsymbol{X}}_{k \mid k}, \boldsymbol{P}_{k \mid k}\right) \end{gather*}$

(35)

$\begin{gather*} q\left(\boldsymbol{P}_{k \mid k-1)} \approx q^{(N)}\left(\boldsymbol{P}_{k \mid k-1}\right)=\operatorname{IW}\left(\boldsymbol{P}_{k \mid k-1} ; {t_{\hat{k}}^{( N)}}, \hat{\boldsymbol{T}}_{k}^{(N)}\right)=\right.\\ \operatorname{IW}\left(\boldsymbol{P}_{k \mid k-1} ; \hat{t}_{k \mid k}, \hat{\boldsymbol{T}}_{k \mid k}\right) \end{gather*}$

(36)

$\begin{gather*} q\left(\boldsymbol{R}_{k}\right) \approx q^{(N)}\left(\boldsymbol{R}_{k}\right)=\operatorname{IW}\left(\boldsymbol{R}_{k} ; \hat{u}_{k}^{(N)}, \hat{\boldsymbol{U}}_{k}^{(N)}\right)= \\ \operatorname{IW}\left(\boldsymbol{R}_{k} ; \hat{u}_{k \mid k}, \hat{\boldsymbol{U}}_{k \mid k}\right) \end{gather*}$

(37)

式(35)~(37)表明，经过N次定点迭代运算后，MVBAKF算法的最终滤波值为：$\hat{\boldsymbol{X}}_{k \mid k}=$ $\hat{\boldsymbol{X}}_{k \mid k}^{(N)}、\boldsymbol{P}_{k \mid k}=\boldsymbol{P}_{k \mid k}^{(N)}、\hat{t}_{k \mid k}=\hat{t}_{k}^{(N)}、\hat{\boldsymbol{T}}_{k \mid k}=\hat{\boldsymbol{T}}_{k}^{(N)}、\hat{u}_{k \mid k}=$ $\hat{u}_{k}^{(N)}、\hat{\boldsymbol{U}}_{k \mid k}=\hat{\boldsymbol{U}}_{k}^{(N)}$；并且，测量噪声方差矩阵R_k、状态向量一步预测协方差矩阵P_k|k－1可以表示为：

$\hat{\boldsymbol{R}}_{k} \approx \boldsymbol{R}_{k}^{(N)}=\frac{\hat{\boldsymbol{U}}_{k}^{(N)}}{\hat{u}_{k}^{(N)}-m-1}$

(38)

$\hat{\boldsymbol{P}}_{k \mid k-1} \approx \boldsymbol{P}_{k \mid k-1}^{(N)}=\frac{\hat{\boldsymbol{T}}_{k}^{(N)}}{\hat{t}_{k}^{(N)}-n-1}$

(39)

3 仿真分析 3.1 实验数据及参数设置

为验证本文MVBAKF算法的有效性与优越性，设计了GNSS/SINS组合导航系统仿真实验，并将MVBAKF算法及传统的VBAKF算法进行比较。为了保证仿真实验结果的严谨性，对上述两种滤波器设置了相同的滤波初值、测量噪声、系统噪声及测量数据(即具有相同的初始条件及相同的输入等)。各导航传感器的误差参数设置如表 1所示。

3.2 仿真结果与分析

设置包含加速、减速、爬升、俯冲及转弯等机动过程的航迹作为导航基准，时长为3 600 s，其飞行轨迹如图 1所示。MVBAKF算法中遗忘因子选取为ρ=0.95，VBAKF算法中遗忘因子向量的每个元素均为0.95，以保证两种算法对历史测量信息的保留度相同。仿真时，在GNSS测量噪声方差中引入突变及缓变两个不同阶段，如图 2中黑色实线部分；根据表 1的设置，将滤波器实际的系统噪声方差设置为真实值的9倍，以模拟滤波过程中未知的系统噪声及测量噪声。

对于MVBAKF算法，当循环次数选取为N=5、调整因子选取为τ=6时，图 2给出基于VBAKF算法和MVBAKF算法的测量噪声均方差估计曲线。可以看出，两种算法对测量噪声方差的估计效果基本相同，估计值基本与真实值一致，且随真实值的变化而变化。

图 3给出基于VBAKF算法、MVBAKF算法的组合导航系统的位置误差对比曲线、速度误差对比曲线及姿态误差对比曲线。为更加直观地评估滤波器对状态参数的滤波精度，选取位置和速度的均方根误差(root mean square error，RMSE)为性能指标，定义如下：

$\begin{array}{c} \mathrm{RMSE} \triangleq\\ \sqrt{\frac{1}{M} \sum\limits_{k=1}^{M}\left[\left(\hat{x}_{k}-x_{k}\right)^{2}+\left(\hat{y}_{k}-y_{k}\right)^{2}+\left(\hat{z}_{k}-z_{k}\right)^{2}\right]} \end{array}$

(40)

式中，(x_k, y_k, z_k)代表k时刻真实的位置信息(或速度信息)，($\hat{x}_{k}, \hat{y}_{k}, \hat{z}_{k}$)代表估计的位置信息(或速度信息)。经统计，VBAKF算法和MVBAKF算法得到的位置RMSE分别为10.8 m和9.7 m，速度RMSE分别为0.304 m/s和0.278 m/s。同理，VBAKF算法和MVBAKF算法得到的横滚角RMSE分别为0.043°和0.041°，俯仰角RMSE分别为0.042°和0.035°，航向角RMSE分别为0.26°和0.24°。相对于VBAKF算法，MVBAKF算法可提高位置精度约10%、可提高速度精度约8.6%，同时姿态角精度也得以提升。结合图 2的估计曲线，说明MVBAKF算法可有效降低未知系统噪声对滤波精度的影响。

为评估MVBAKF算法中调整因子τ对滤波精度的影响，在循环次数N=5不变的基础上，图 4给出了τ=2、3、4、5、6、7、8、9、10时的位置均方根误差及速度均方根误差曲线。由图 4可见，当调整因子较小时(如τ=2, 3, 4)，位置均方根误差及速度均方根误差较大；当调整因子较小时(如τ=9, 10)，位置均方根误差及速度均方根误差较小；当调整因子继续增大时(τ>10)，位置均方根误差及速度均方根误差基本稳定。

为评估MVBAKF算法中循环次数N对滤波精度的影响，在调整因子为τ=10不变的基础上，图 5给出了N=1~10时的位置均方根误差及速度均方根误差曲线。

可以看出，当循环次数N=2~4时，位置均方根误差及速度均方根误差达到较小的值，否则位置均方根误差及速度均方根误差较大。

4 结语

在实际环境中，针对具有不准确测量噪声及系统噪声统计特性的线性高斯状态空间GNSS/SINS组合导航系统滤波，提出一种改进的变分贝叶斯自适应滤波算法。在该算法中，将对系统噪声的估计问题变换成对状态向量估计的一步预测误差协方差估计问题，进而将状态估计及其一步预测误差协方差、测量噪声方差通过逆Wishart分布进行推断。

仿真结果表明，与VBAKF算法相比，本文MVBAKF算法具有更好的鲁棒性，且能提高组合导航系统的滤波精度，这是由于MVBAKF算法可以迭代地找到更优的状态估计一步预测误差协方差及测量噪声方差。