以伪距、多谱勒频移和载波相位为观测量的全球导航卫星系统/捷联惯导（Global Navigation Satellite System/ strap-down Inertial Navigation System，GNSS/SINS）紧组合系统因其显著的灵活性和稳定性，受到众多大地测量学者的关注。但载波相位观测中，整周模糊度的确定是个难点，为避免其求解，有学者提出了时间差分载波相位（time differenced carrier phase, TDCP），并在实际应用中取得了较好效果^[1-2]。现实导航存在典型的非高斯噪声，如随机模型采用高斯白噪声近似会引起滤波的稳定性和精度下降^[3]。利用高斯混合模型（Gaussian mixture model，GMM）对非高斯噪声建模比传统采用膨大方差的高斯分布近似方法，在图像处理、语音识别和导航完好性监测等方面具有更高的精度^[4-6]。但GMM参数的估计较为复杂，目前GMM参数估计方面的研究主要集中在最大似然估计（maximum likelihood estimate, MLE）、马尔可夫链蒙特卡洛（Markov chain Monte Carlo, MCMC）和期望最大化(expectation maximization, EM)等，其中MLE和EM算法不能自适应估计模型阶数且易出现过拟合现象，MCMC算法计算量大、收敛速度慢，不适合实时性要求较高的应用^[7]。如果GMM参数不准确，会直接影响高斯混合滤波（Gaussian mixture filter, GMF）的估计精度。为此，有学者提出了包含量测噪声的最优贝叶斯滤波，但在实际应用中面临被积函数复杂、变量维数高的问题，且噪声方差阵的引入使其求解更加困难，从而影响其推广使用^[8]。

一种源于机器学习领域的变分贝叶斯(variational Bayesian, VB)以较少运算量进行较高精度的参数估计，可解决贝叶斯估计在实际应用中被积函数计算复杂的问题，适应于高维环境。在参数估计、盲源分离、语音识别、独立成分分析等领域取得了一定的进展^[9-10]。

本文受上述研究启发，考虑到伪距和时间差分载波相位观测量中噪声的非高斯分布特性，提出了一种基于VB的高斯混合自适应滤波，解决了GMM参数的自适应估计问题，进一步优化了GMF算法。通过实验分析验证了该算法在计算精度、实时性和鲁棒性方面的贡献，为其今后在数据处理领域的扩展提供了一定的理论支持。

1 紧组合导航高斯混合滤波模型 1.1 伪距和时间差分载波相位的组合导航

GNSS/SINS紧组合导航的状态向量x为包含姿态、速度、位置、陀螺漂移、加速度计漂移、时钟偏置和时钟漂移的17维列向量，状态方程可建模为：

${\mathit{\boldsymbol{x}}_{t, t - 1}} = f\left( {{\mathit{\boldsymbol{x}}_{t - 1}}} \right) + \mathit{\boldsymbol{gw}}$

(1)

式中，f(·)为非线性传递函数；g为噪声系数阵；w为系统噪声；t为滤波时刻。

量测方程由伪距和TDCP观测量构成，经卫星钟差和大气层各项改正后，整理得：

${\rm{ \mathsf{ δ} }}\rho = c{\rm{ \mathsf{ δ} }}{t_r} + {e_1}{\rm{ \mathsf{ δ} }}x + {e_2}{\rm{ \mathsf{ δ} }}y + {e_3}{\rm{ \mathsf{ δ} }}z + \varepsilon _\rho ^m$

(2)

${\rm{\Delta }}{\varphi ^m} = {\rm{\Delta }}{r^m} + c{\rm{\Delta \mathsf{ δ} }}{t_r} + {\rm{\Delta }}\varepsilon _\varphi ^m$

(3)

式中，δρ为伪距观测量；cδt_r为接收机r的钟差等效距离误差；e为线性化系数；δx、δy和δz为位置估值误差；ε_ρ^m为第m颗卫星到接收机之间由热噪声、多路径和轨道预测残差等引起的噪声残余项；Δφ^m为第m颗卫星的TDCP观测量；cΔδt_r为相邻历元的接收机钟差等效距离差分量；Δr^m为相邻历元间接收机到卫星的几何距离差分量；Δε_φ^m为噪声残余项。

合并式（2）和式（3），得量测方程：

$\mathit{\boldsymbol{z}} = h\left( {{\mathit{\boldsymbol{x}}_{t, t - 1}}} \right) + \mathit{\boldsymbol{v}}$

(4)

式中，z为量测向量；h(·)为量测函数；v为其对应的噪声项。

文献[11-12]通过小波变换和Allan方差对SINS和GNSS噪声进行分析，发现其不符合高斯白噪声特性，是高斯与非高斯分布的混合状态，且差分后增加了噪声之间的相关性，若采用多系统，则其统计特性更为复杂。因此，为改进组合导航滤波的性能，必须对其随机模型采取更为精化的模型进行处理。

1.2 高斯混合滤波

高斯混合滤波采用高斯概率分布逼近的方法，使复杂的非高斯系统可由若干高斯分量以任意精度近似^[4]。对非高斯噪声ε建模，有：

$p\left( \mathit{\boldsymbol{\varepsilon }} \right) = \mathop {\mathop \sum\limits_{i = 1} }\limits^I \mathit{\boldsymbol{\omega }}_{^t}^i \times N\left[ {\mathit{\boldsymbol{\mu }}_t^i, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_t^i} \right]$

(5)

$\mathop {\mathop \sum\limits_{i = 1} }\limits^I \mathit{\boldsymbol{\omega }}_{^t}^i = 1$

(6)

式中，p(ε)表示全局概率密度分布函数; N[]为正态分布；μ_tⁱ和Σ_tⁱ为t时刻第i个分量的均值和协方差；ω_tⁱ为对应的权向量。这样，原滤波器被分解为多个并行滤波器，且每个滤波器都符合高斯白噪声的要求，按最优贝叶斯估计准则，整理得到状态参数的最优后验概率分布：

$p({{\mathit{\boldsymbol{\hat X}}}_t}|{\mathit{\boldsymbol{L}}_t}) = \mathop {\mathop \sum\limits_{i = 1} }\limits^I \mathit{\boldsymbol{\omega }}_{^t}^iN\left[ {{{\mathit{\boldsymbol{\hat X}}}_t} - \mathit{\boldsymbol{\hat X}}_t^i, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_t}}^i} \right]$

(7)

式中，

$\mathit{\boldsymbol{\omega }}_{^t}^i = \frac{{\mathit{\boldsymbol{\omega }}_{t, t - 1}^i\beta _{^t}^i}}{{\mathop {\mathop \sum\limits_{i = 1} }\limits^I \left( {\mathit{\boldsymbol{\omega }}_{t, t - 1}^i\beta _{^t}^i} \right)}} = \frac{{\mathit{\boldsymbol{\omega }}_{t - 1}^i\beta _{^t}^i}}{{\mathop {\mathop \sum\limits_{i = 1} }\limits^I \left( {\mathit{\boldsymbol{\omega }}_{t - 1}^i\beta _{^t}^i} \right)}}$

(8)

$\beta _{^t}^i \approx N[{\mathit{\boldsymbol{L}}_t} - {\mathit{\boldsymbol{H}}_t}{{\mathit{\boldsymbol{\bar X}}}_t}, {\mathit{\boldsymbol{H}}_t}\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\bar X}_t}}^{}\left( {{\mathit{\boldsymbol{H}}_t}{)^{\rm{T}}} + \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{_t}^i} \right]$

(9)

式中，L_t为观测向量；H_t为量测系数矩阵；$\mathit{\boldsymbol{\hat X}}_t^i$为状态估值；X_tⁱ为状态预测值；I为高斯分量总数。由此可见，滤波更新中高斯分量权值可随新息自适应调整，使得算法具有一定的鲁棒性。但GMM初始参数$\{ \mathit{\boldsymbol{\omega }}_t^i, \mathit{\boldsymbol{\mu }}_t^i, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_t^i\} _{i = 1}^I$的确定是个难点，为此本文借鉴VB思路来解决其参数估计问题，平衡计算精度与效率上的矛盾。

2 VB优化高斯混合滤波

VB学习的基本思想是将复杂的多维数值积分问题转化为求解近似分布的优化问题。在贝叶斯估计准则中，观测数据集合为Y，Θ为d个未知模型参数集合Θ={θ₁, θ₂⋯θ_d}，Z为隐变量。当后验分布p(Θ, Z|Y)形式复杂难以得到解析表达式时，可通过最小化Kullback-Leible散度使新分布q(Θ, Z)逼近真实后验分布p(Θ, Z|Y)。由变分泛函数理论，得到通解表达式为：

$q\left( {{\mathit{\boldsymbol{\theta }}_i}} \right) = \frac{{{\rm{exp}}\left\{ {{E_q}\left( {{\mathit{\boldsymbol{\theta }}_{j \ne i}}} \right)\left[ {{\rm{ln}}p\left( {Y, {\rm{\Theta }}, \mathit{\boldsymbol{Z}}} \right)} \right]} \right\}}}{{\mathop \smallint \nolimits^ {\rm{exp}}\left\{ {{E_q}\left( {{\mathit{\boldsymbol{\theta }}_{j \ne i}}} \right)\left[ {{\rm{ln}}p\left( {Y, {\rm{\Theta }}, \mathit{\boldsymbol{Z}}} \right)} \right]} \right\}{\rm{d}}{\mathit{\boldsymbol{\theta }}_i}}}{\rm{}}$

(10)

式中，E_x(·)表示x条件下的期望。共轭指数域函数使得q(θ_i)的分布形式同真实后验分布，其中θ_i为Θ中任意元素。式(10)为隐式解，每个θ_i的计算与θ_j≠i有关，需迭代计算：

${q_n}\left( \mathit{\boldsymbol{Z}} \right) \propto {\rm{exp}}\left\{ {{E_{q\left( \mathit{\Theta } \right)}}\left[ {{\rm{ln}}p\left( {Y, \mathit{\boldsymbol{Z}}, \mathit{\Theta } } \right)} \right]} \right\}$

(11)

${q_n}\left( \mathit{\Theta } \right) \propto {\rm{exp}}\left\{ {{E_{q\left( \mathit{\boldsymbol{Z}} \right)}}\left[ {{\rm{ln}}p\left( {Y, \mathit{\boldsymbol{Z}}, \mathit{\Theta } } \right)} \right]} \right\}$

(12)

n次迭代后，计算收敛^{[7, 9]}。

2.1 高斯混合模型

GMM由K个高斯分布线性组合而成:

$p\left( {Y\left| {\mathit{\boldsymbol{\pi }}, \mathit{\boldsymbol{\mu }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right.} \right) = {\mathit{\boldsymbol{\pi }}_k}N\left( {Y\left| {{\mathit{\boldsymbol{\mu }}_k}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k^{ - 1}} \right.} \right)$

(13)

式中，π=[π₁ π₂ ⋯ π_K]表示权值向量；μ_k和Λ_k分别为第k个高斯分布的均值和协方差。引入隐变量Ζ表征每个D维的观测样本y_u来源于某个特定的高斯分布，则K个高斯分量的隐变量Z可由U×K阶矩阵表示，U为观测样本个数。隐变量Z的条件密度函数为：

$p\left( {\mathit{\boldsymbol{Z}}\left| \mathit{\boldsymbol{\pi }} \right.} \right) = \mathop {\mathop \prod\limits_{u = 1} }\limits^U \mathop {\mathop \prod\limits_{k = 1} }\limits^K \pi _k^{{z_{uk}}}$

(14)

式中，参数z_uk取0或1。

共轭指数分布作为常用的一类参数先验分布，与其后验属同一类分布族，故权向量π的先验分布为狄里克利（Dirichlet）分布，即：

$p\left( \mathit{\boldsymbol{\pi }} \right) = {\rm{Dir}}\left( {\mathit{\boldsymbol{\pi }}\left| {{\mathit{\boldsymbol{\alpha }}^0}} \right.} \right) = C\left( {{\mathit{\boldsymbol{\alpha }}^0}} \right)\mathop {\mathop \prod\limits_{k = 1} }\limits^K \mathit{\boldsymbol{\mu }}_k^{{\alpha _k} - 1}$

(15)

设K个混合分量具有相同先验分布，令α₁=…=α_K=α₀，α⁰是具有K个元素的矢量，α⁰=[α₀ ⋯ α₀]^T。则带有隐变量且每个观测样本独立同分布的GMM概率可表示为：

$\begin{array}{l} p\left( {Y\left| {\mathit{\boldsymbol{Z}}, \mathit{\boldsymbol{\pi }}, \mathit{\boldsymbol{\mu }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right.} \right) = \mathop {\mathop \prod\limits_{u = 1} }\limits^U \mathop {\mathop \prod\limits_{k = 1} }\limits^K \left\{ {\frac{1}{{{{(2{\pi _k})}^{D/2}}}}{{\left| {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}} \right|}^{1/2}}} \right. \times \\ \;\;\;\;\left. {{\rm{exp}}[ - \frac{1}{2}\left( {{\mathit{\boldsymbol{y}}_u} - {\mathit{\boldsymbol{\mu }}_k}{)^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}\left( {{\mathit{\boldsymbol{y}}_u} - {\mathit{\boldsymbol{\mu }}_k}} \right)} \right]} \right\}{^{{z_{uk}}}} \end{array}$

(16)

式中，D为常量。其均值μ和协方差Λ的共轭先验分布为条件高斯分布与威沙特（Wishart）分布的乘积，即：

$\begin{array}{*{20}{c}} {p\left( {\mathit{\boldsymbol{\mu }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right) = p\left( {\mathit{\boldsymbol{\mu }}\left| \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} \right.} \right)p\left( \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} \right) = }\\ {\mathop {\mathop \prod\limits_{k = 1} }\limits^K N\left( {{\mathit{\boldsymbol{\mu }}_k}\left| {{\mathit{\boldsymbol{m}}_0}, {{\left( {{\beta _0}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}} \right)}^{ - 1}}} \right.} \right)W\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}\left| {{\mathit{\boldsymbol{W}}_0}, {v_0}} \right.} \right)} \end{array}$

(17)

式中，W₀为正定阵；m₀为K个元素的矢量；β₀为调节参数；v₀为自由度参数。

由此可见，p(Y|π, μ, Λ)的计算需要对Z、π、Λ和μ这4个参数进行复杂的积分运算。借鉴VB思路，只需引入q(Θ)，在共轭指数域内选取分布函数得到闭合解析式，循环迭代后即可求出最优参数。

2.2 VB优化的模型参数估计

设q(Z, π, Λ, μ)为近似解的联合密度函数，将隐变量与各参数的独立分布因子化，得：

$q\left( {\mathit{\boldsymbol{Z}}, \mathit{\boldsymbol{\pi }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}, \mathit{\boldsymbol{\mu }}} \right) = q\left( \mathit{\boldsymbol{Z}} \right)q\left( \mathit{\boldsymbol{\pi }} \right)\mathop {\mathop \prod\limits_{k = 1} }\limits^K q\left( {{\mathit{\boldsymbol{\mu }}_k}\left| {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}} \right.} \right)q\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}} \right)$

(18)

可由VB思路，通过式（11）计算隐变量Z的近似后验分布：

$\begin{array}{l} q\left( \mathit{\boldsymbol{Z}} \right) \propto {\rm{exp}}\left\{ {\mathop {\mathop \sum\limits_{u = 1} }\limits^U \mathop {\mathop \sum\limits_{k = 1} }\limits^K {z_{uk}}[\frac{1}{2}{E_\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}}\left( {{\rm{ln}}\left| {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}} \right|} \right) - } \right.\frac{D}{2}{\rm{ln}}2{\pi _k} - \\ \left. {\frac{1}{2}{E_{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}, \mathit{\boldsymbol{\mu }}}}[{{({\mathit{\boldsymbol{y}}_u} - {\mathit{\boldsymbol{\mu }}_k})}^T}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}\left( {{\mathit{\boldsymbol{y}}_u} - {\mathit{\boldsymbol{\mu }}_k}} \right)\left] { + {E_{{\pi _k}}}\left( {{\rm{ln}}{\pi _k}} \right)} \right]} \right\} \end{array}$

(19)

方便起见，记式（19）为$q\left( \mathit{\boldsymbol{Z}} \right) \propto \mathop {\mathop \prod\limits_{u = 1} }\limits^U \mathop {\mathop \prod\limits_{k = 1} }\limits^K \gamma _{uk}^{{z_{uk}}}$，其中γ_uk需归一化处理。同理，由式（12）得：

$\begin{array}{*{20}{c}} {q\left( {\mathit{\boldsymbol{\pi }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}, \mathit{\boldsymbol{\mu }}} \right) \propto {\rm{exp}}[\mathop {\mathop \sum\limits_{u = 1} }\limits^U \mathop {\mathop \sum\limits_{k = 1} }\limits^K {E_{{z_{uk}}}}\left[ {{z_{uk}}} \right] \cdot }\\ {\left( {{\rm{ln}}N\left( {{\mathit{\boldsymbol{y}}_u}\left| {{\mathit{\boldsymbol{\mu }}_k}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k^{ - 1}} \right.} \right) + {\rm{ln}}{\pi _k}} \right) + {\rm{lnDir}}\left( {\mathit{\boldsymbol{\pi }}\left| {{\mathit{\boldsymbol{\alpha }}^0}} \right.} \right) + }\\ {{\rm{ln}}\left( {N\left( {{\mathit{\boldsymbol{\mu }}_k}\left| {{\mathit{\boldsymbol{m}}_0}, {{({\beta _0}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k})}^{ - 1}}} \right.} \right)W\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}\left| {{\mathit{\boldsymbol{W}}_0}, {v_0}} \right.} \right)} \right)]} \end{array}$

(20)

因z_uk是离散型变量，式（11）为Z的近似分布，则每个z_uk的期望为：

${E_{{z_{uk}}}}\left[ {{z_{uk}}} \right] = 1 \times {\gamma _{uk}} = {\gamma _{uk}}$

(21)

式中, 分量权值、分量协方差和分量均值3个参数相互独立。对于权值参数π_k，将其当作另一个Dirichlet分布：

$q\left( \mathit{\boldsymbol{\pi }} \right) = C\mathop {\mathop \prod\limits_{k = 1} }\limits^K \pi _k^{\left( {{\alpha _0} + \mathop {\mathop \sum\limits_{u = 1} }\limits^U {\gamma _{uk}} - 1} \right)} = {\rm{Dir}}\left( {\mathit{\boldsymbol{\pi }}\left| {{\mathit{\boldsymbol{\alpha }}^0}} \right.} \right)$

(22)

式中，C为常量。为方便起见，令：

则

${\alpha _k} = {\alpha _0} + {N_k}$

(23)

参照式（18），分别处理q(μ_k|Λ_k)与q(Λ_k)推导分量协方差和分量均值的近似分布：

$q\left( {{\mathit{\boldsymbol{\mu }}_k}\left| {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}} \right.} \right) \propto N({\mathit{\boldsymbol{\mu }}_k}\left| {{\mathit{\boldsymbol{m}}_k}} \right., \left( {{\beta _k}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}{)^{ - 1}}} \right)$

(24)

式中，β_k=β₀+N_k；${\mathit{\boldsymbol{m}}_k} = \frac{{{\beta _0}{\mathit{\boldsymbol{m}}_0} + {N_k}{{\mathit{\boldsymbol{\bar y}}}_k}}}{{{\beta _0} + {N_k}}}$。

整理得到新的Wishart分布：

$q\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}} \right) \propto W\left( {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}\left| {{v_k}, {\mathit{\boldsymbol{W}}_k}} \right.} \right)$

(25)

式中，v_k=v₀+N_k；$\mathit{\boldsymbol{W}}_k^{ - 1} = \mathit{\boldsymbol{W}}_0^{ - 1} + {N_k}{\mathit{\boldsymbol{S}}_k} + \frac{{{\beta _0}{N_k}}}{{{\beta _0} + {N_k}}}\left( {{{\mathit{\boldsymbol{\bar y}}}_k} - {\mathit{\boldsymbol{m}}_0}} \right){({{\mathit{\boldsymbol{\bar y}}}_k} - {\mathit{\boldsymbol{m}}_0})^{\rm{T}}}$。最终可由高斯分布、Dirichlet分布和Wishart分布得式（19）中E_{Λ, μ}[(y_u-μ_k)^TΛ_k(y_u-μ_k)]、E_Λ(ln|Λ_k|)、E_πk(lnπ_k)各项为：

$\begin{array}{*{20}{c}} {{E_{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}, \mathit{\boldsymbol{\mu }}}}[\left( {{\mathit{\boldsymbol{y}}_u} - {\mathit{\boldsymbol{\mu }}_k}{)^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}\left( {{\mathit{\boldsymbol{y}}_u} - {\mathit{\boldsymbol{\mu }}_k}} \right)} \right]{\rm{}} = }\\ {{v_k}{{({\mathit{\boldsymbol{y}}_u} - {\mathit{\boldsymbol{m}}_k})}^{\rm{T}}}{\mathit{\boldsymbol{W}}_k}\left( {{\mathit{\boldsymbol{y}}_u} - {\mathit{\boldsymbol{m}}_k}} \right) + D\beta _k^{ - 1}} \end{array}$

(26)

$\begin{array}{*{20}{l}} {{E_\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}}\left( {{\rm{ln}}\left| {{\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}_k}} \right|} \right) = \mathit{\Psi }\left( {\frac{{{v_k} + 1 - i}}{2}} \right) + }\\ {{\rm{}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;D{\rm{ln}}2 + {\rm{ln}}\left( {{\mathit{\boldsymbol{W}}_k}} \right)} \end{array}$

(27)

${E_{{{\rm{ \mathsf{ π} }}_k}}}\left( {{\rm{ln}}{\pi _k}} \right) = \mathit{\Psi }\left( {{\alpha _k}} \right) - \mathit{\Psi }\left( {\mathop {\mathop \sum\limits_{k = 1} }\limits^K {\alpha _k}} \right)$

(28)

式中，Ψ(·)为双伽马函数。

2.3 判断收敛下限的方法

将式（18）代入F(q(Θ, Z))得：

$\begin{array}{*{20}{c}} {F\left( {q\left( {\mathit{\boldsymbol{Z}}, \mathit{\boldsymbol{\pi }}, \mathit{\boldsymbol{\mu }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right)} \right) = E[{\rm{ln}}p\left( {Y\left| {\mathit{\boldsymbol{Z}}, \mathit{\boldsymbol{\mu }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right.} \right) + }\\ {{\rm{ln}}p\left( {\mathit{\boldsymbol{\mu }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right) + {\rm{ln}}p\left( {\mathit{\boldsymbol{Z}}\left| \mathit{\boldsymbol{\pi }} \right.} \right) + {\rm{ln}}p\left( \mathit{\boldsymbol{\pi }} \right)] - }\\ {{\rm{E}}\left[ {{\rm{ln}}q\left( \mathit{\boldsymbol{\pi }} \right) + {\rm{ln}}q\left( \mathit{\boldsymbol{Z}} \right) + {\rm{ln}}q\left( {\mathit{\boldsymbol{\mu }}, \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right)} \right]} \end{array}$

(29)

式中，F(·)为判断收敛性下界函数，通过式（26）至式（28）的性质，计算式（29）中的各项期望。其余期望代入每个近似分布的负熵值计算，熵计算公式为：

${H_{p\left( \mathit{\Theta } \right)}} = - \mathop \smallint \nolimits^ p\left( \mathit{\Theta } \right){\rm{ln}}p\left( \mathit{\Theta } \right){\rm{d}}\mathit{\Theta }$

(30)

每次迭代后计算下界变化值ΔF，当该值低于预先设定的阈值时，则算法收敛，即求得逼近原分布的近似分布。

2.4 VB优化的GMF算法实现

具体滤波算法实现如图 1所示，其中EKF表示扩展卡尔曼滤波(extended Kalman filter)。算法步骤为：

1）根据状态方程式（1）和量测方程式（4）分别对其噪声项进行GMM建模。

2）贝叶斯学习阶段，由式（19）更新隐变量，由式（20）更新参数集，通过式（29）判断收敛，离线计算GMM参数。

3）滤波更新阶段，利用式（5）至式（9）完成自适应滤波估计，得到后验密度。

4）借鉴文献[6]中拟贝叶斯近似法对混合分量缩减，控制更新后高斯分量数。

3 计算与分析 3.1 模拟算例

TDCP/SINS紧组合系统仿真参数设置为：陀螺仪刻度系数0.3%，零偏2°/$\sqrt {\rm{h}} $，随机游走0.1°/s；加速度计刻度系数0.3%，零偏30 μg/$\sqrt {\rm{Hz}} $，随机游走0.1 mg；初始姿态误差协方差20″，初始速度协方差10 cm/s, 初始位置误差协方差3 m；接收机钟差和钟漂协方差分别为10 m和1 m/s。计算平台为Inter Core i5-7200U, RAM 4 GB, Matlab 2014ra。采用5种方案进行50次仿真测试。

1) 方案1：EKF算法。

2) 方案2：文献[10]中的抗差自适应EKF算法。

3) 方案3：基于MCMC的GMEKF（Gaussian mixture extended Kalman filter）算法，样本容量为1 000。

4) 方案4：基于MCMC的GMEKF，样本容量为5 000。

5) 方案5：基于VB的GMEKF，设置初始参数α₀=1，β₀=1，v₀=1，m₀=2，W₀中各元素设置为2I，I=2，K=4, 收敛阈值设为ΔF < 1×10^-6[9]，并将权值小于1%的分量删减。

通过p(v)=εN(μ₁, R₁)+(1-ε)N(μ₂, R₂)建模非高斯混合噪声，其中v为混合噪声项，ε为高斯分量的比例系数，μ和R为对应的均值和协方差。

3.1.1 仿真实验1

据文献[13]可观测性分析，加速运动或高机动状态的载体姿态角可观性较好。模拟载体在初始位置北纬34.4°、东经113.5°匀速运动300 s，并在50~250 s发生加速和转弯的过程。将非高斯噪声置于GNSS观测量中，取ε=0.1，μ₁=μ₂=0，R₁=3R₂=R，R=diag（0.3², 0.3²…0.3²），并每隔20历元附加2.5 m的粗差。结果如图 2所示，用均方根误差（root mean square error，RMSE）做定量比较，如表 1所示。

3.1.2 仿真实验2

模拟载体以60 m/s匀速运动120 s，并发生两段时长30 s的卫星失锁状况，非高斯噪声设置同仿真实验1，结果如图 3和表 2所示。

3.1.3 仿真实验3

为验证算法的一般性，设置高斯混合分布的不同比例系数ε，仿真实验1和仿真实验2的最大定姿误差和最大定位误差分别如图 4(a)、4(b)所示。

3.1.4 实验结果分析

1）仿真实验1表明，除方案1外，其余方案均能抑制受粗差和非高斯噪声污染时随机模型对定姿的影响。方案2采用自适应因子和观测信息等价协方差矩阵，虽具有一定抗差性，但仍基于高斯白噪声假设，限制了估计精度。方案3、4和5通过GMM对非高斯噪声建模，精化了观测量随机模型，从而进一步提高了滤波量测更新精度。

2）仿真实验2表明，没有GNSS辅助，SINS单独工作导致定位误差累积快。而方案3、4和5采用GMM对动力学模型中的非高斯噪声建模，改善了滤波时间更新精度。

3）前两组仿真实验中，方案3的精度逊色于方案4，因MCMC估计受样本容量限制，当容量少时，不能较好地拟合GMM参数，导致GMF精度降低，而容量大又会影响实时性。方案5的精度接近于方案4，耗时明显减少，可见本文算法通过VB能有效平衡滤波精度和效率上的矛盾。

4）仿真实验3中，当ε取0时，随机模型呈高斯白噪声，方案1误差较小。当ε取值发生变化时，方案1结果恶化，而方案4、5无论在何种取值下都能较好地估计GMM参数，使GMF取得较好效果。

3.2 动态数据算例

机动载体采用微惯导（micro electro mechanical system interial measurement unit，MEMS IMU）与激光惯导（ring laser gyro IMU，RLG IMU），其标称精度见表 3，将后者输出结果作为参考值。沿途有建筑物、树林和水面，数据采集环境较为复杂，可导致明显的多路径效应并使系统失去高斯白噪声特性。为验证算法性能，截取一段1 000 s的数据，用§3.1中的方案1和方案5进行分析，结果如图 5和表 4所示。

对比实验结果，方案5的姿态角精度整体得以改善，但基于EKF的滤波框架中，高阶截断误差会引起部分跳跃点。方案5的俯仰角误差在少部分历元时精度比方案1差，是由于受载体机动影响，GMM参数学习对初值敏感，进而影响结果。从表 4中发现，本文算法作为一种高效的自适应滤波，可以明显改善姿态角精度，尤其是航偏角，提高了约1.8倍，每历元耗时为0.1 s，能保持实时性。由此证明，高机动状态下，观测条件较差时，使用基于VB优化的GMF，可较好地精化随机模型，适合TDCP/SINS的动态导航数据处理。

4 结语

为了进一步提高时间差分载波相位观测量的GNSS/SINS数据处理性能，平衡滤波在算法精度和计算效率上的矛盾，本文通过对非高斯噪声影响下的随机模型进行精化处理，利用VB理论对GMM参数进行估计，构建了一种基于VB学习的混合高斯自适应滤波算法。实验结果表明，相比传统GMF算法，本文算法有效实现了混合模型参数的自适应估计，在非高斯噪声条件下有着更好的鲁棒性，同时借鉴VB学习优化了算法过程，计算负担小，具有较好的工程应用价值。