Gaussian Mixture Filter Based on Variational Bayesian Learning Optimization and Its Application to Integrated Navigation
-
摘要: 针对时间差分载波相位/捷联惯导紧组合系统在非高斯噪声环境工作时，采用高斯混合滤波遇到的混合模型参数估计问题，提出了一种变分贝叶斯学习优化的高斯混合自适应滤波算法。该算法借鉴变分学习理论，准确高效地实现了高斯混合模型参数的自适应估计，进一步精化了滤波算法中的随机模型，能够显著提高估计精度，降低计算负担，改善滤波性能。实验结果表明，相比传统滤波算法，该算法的估计精度得到了进一步改善，运算耗时仅与拓展卡尔曼滤波相当。Abstract: We present a novel Gaussian mixture filter (GMF) improved by variational Bayesian learning. The method is mainly used in time-differenced carrier phase/strap-down Inertial Navigation System integrated navigation to deal with parameter estimation of Gaussian mixture model in non-Gaussian noise environment. By means of variational Bayesian learning theory, this algorithm accurately and efficiently estimates parameters of GMF, and further refines the stochastic model. And it improves accuracy of GMF, reduces the computational complexity, and improves the effectiveness. The experimental results illustrate that the proposed filter substantially outperforms the existing algorithm in terms of estimation accuracy, and is computationally much more efficient. It provides theoretical support for integrated navigation data fusion strategy based on GMF.
-
以伪距、多谱勒频移和载波相位为观测量的全球导航卫星系统/捷联惯导（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δtr为接收机r的钟差等效距离误差；e为线性化系数；δx、δy和δz为位置估值误差；ερm为第m颗卫星到接收机之间由热噪声、多路径和轨道预测残差等引起的噪声残余项；Δφm为第m颗卫星的TDCP观测量；cΔδtr为相邻历元的接收机钟差等效距离差分量；Δrm为相邻历元间接收机到卫星的几何距离差分量；Δεφ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[]为正态分布；μti和Σti为t时刻第i个分量的均值和协方差；ωti为对应的权向量。这样，原滤波器被分解为多个并行滤波器，且每个滤波器都符合高斯白噪声的要求，按最优贝叶斯估计准则，整理得到状态参数的最优后验概率分布：
$$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) 式中，Lt为观测向量；Ht为量测系数矩阵；$\mathit{\boldsymbol{\hat X}}_t^i$为状态估值；Xti为状态预测值；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个未知模型参数集合Θ={θ1, θ2⋯θ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) 式中，Ex(·)表示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) 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) 式中，π=[π1 π2 ⋯ πK]表示权值向量；μk和Λk分别为第k个高斯分布的均值和协方差。引入隐变量Ζ表征每个D维的观测样本yu来源于某个特定的高斯分布，则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) 式中，参数zuk取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个混合分量具有相同先验分布，令α1=…=αK=α0，α0是具有K个元素的矢量，α0=[α0 ⋯ α0]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) 式中，W0为正定阵；m0为K个元素的矢量；β0为调节参数；v0为自由度参数。
由此可见，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) 因zuk是离散型变量，式（11）为Z的近似分布，则每个zuk的期望为：
$${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为常量。为方便起见，令：
$$\mathop {\mathop \sum\limits_{u = 1} }\limits^U {\gamma _{uk}} = {N_k}$$ $${{\mathit{\boldsymbol{\bar y}}}_k} = \frac{1}{{{N_k}}}\mathop {\mathop \sum\limits_{u = 1} }\limits^U {\gamma _{uk}}{\mathit{\boldsymbol{y}}_u}$$ $${\mathit{\boldsymbol{S}}_k} = \frac{1}{{{N_k}}}\mathop {\mathop \sum\limits_{u = 1} }\limits^U {\gamma _{uk}}\left( {{\mathit{\boldsymbol{y}}_u} - {{\mathit{\boldsymbol{\bar y}}}_k}} \right){({\mathit{\boldsymbol{y}}_u} - {{\mathit{\boldsymbol{\bar y}}}_k})^{\rm{T}}}$$ 则
$${\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=β0+Nk；${\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) 式中，vk=v0+Nk；$\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Λ, μ[(yu-μk)TΛk(yu-μ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 变分优化的高斯混合EKF滤波算法框架Figure 1. Flowchart of VB Based on Gaussian Mixture Extended Kalman Filter1）根据状态方程式（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，设置初始参数α0=1，β0=1，v0=1，m0=2，W0中各元素设置为2I，I=2，K=4, 收敛阈值设为ΔF < 1×10-6[9]，并将权值小于1%的分量删减。
通过p(v)=εN(μ1, R1)+(1-ε)N(μ2, R2)建模非高斯混合噪声，其中v为混合噪声项，ε为高斯分量的比例系数，μ和R为对应的均值和协方差。
3.1.1 仿真实验1
据文献[13]可观测性分析，加速运动或高机动状态的载体姿态角可观性较好。模拟载体在初始位置北纬34.4°、东经113.5°匀速运动300 s，并在50~250 s发生加速和转弯的过程。将非高斯噪声置于GNSS观测量中，取ε=0.1，μ1=μ2=0，R1=3R2=R，R=diag（0.32, 0.32…0.32），并每隔20历元附加2.5 m的粗差。结果如图 2所示，用均方根误差（root mean square error，RMSE）做定量比较，如表 1所示。
表 1 仿真1姿态误差结果（RMSE）Table 1. Attitude Error of Simulation 1(RMSE)方案 俯仰角/（°） 横滚角/（°） 航偏角/（°） 每历元耗时/s 1 0.582 9 0.673 8 0.872 7 0.084 2 0.314 5 0.378 4 0.734 7 0.092 3 0.227 6 0.339 9 0.583 1 0.206 4 0.164 8 0.099 8 0.472 7 0.530 5 0.197 4 0.107 3 0.474 6 0.108 3.1.2 仿真实验2
模拟载体以60 m/s匀速运动120 s，并发生两段时长30 s的卫星失锁状况，非高斯噪声设置同仿真实验1，结果如图 3和表 2所示。
表 2 卫星失锁状态最大定位误差结果Table 2. Maximum Positioning Error Under Satellite Losslock方案 失锁20~50 s 失锁70~100 s 每历元耗时/s 水平误差/m 竖直误差/m 水平误差/m 竖直误差/m 1 60.49 2.68 60.78 2.75 0.098 2 58.72 2.01 58.33 2.03 0.114 3 56.86 1.81 56.17 1.72 0.240 4 51.40 1.43 52.93 1.47 0.608 5 51.41 1.48 52.83 1.46 0.124 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所示。
表 3 惯性器件精度Table 3. Inertia Device Accuracy惯性器件 刻度系数 零偏 随机游走 MEMS 陀螺仪 0.2% 2 (°)/$\sqrt {\rm{h}} $ 0.1 (°)/s IMU 加速度计 0.3% 36 μg/$\sqrt {\rm{Hz}} $ 0.1 mg RLG 陀螺仪 5×10-6 0.002 (°)/$\sqrt {\rm{h}} $ 0.003 (°)/s IMU 加速度计 100×10-6 8 μg/$\sqrt {\rm{Hz}} $ 25 μg 表 4 动态数据姿态误差结果（RMSE）Table 4. Attitude Error of Dynamic Data(RMSE)方案 俯仰角/（°） 横滚角/（°） 航偏角/（°） 每历元耗时/s 1 0.425 9 0.182 7 3.686 1 0.086 5 0.389 9 0.102 1 2.589 6 0.108 对比实验结果，方案5的姿态角精度整体得以改善，但基于EKF的滤波框架中，高阶截断误差会引起部分跳跃点。方案5的俯仰角误差在少部分历元时精度比方案1差，是由于受载体机动影响，GMM参数学习对初值敏感，进而影响结果。从表 4中发现，本文算法作为一种高效的自适应滤波，可以明显改善姿态角精度，尤其是航偏角，提高了约1.8倍，每历元耗时为0.1 s，能保持实时性。由此证明，高机动状态下，观测条件较差时，使用基于VB优化的GMF，可较好地精化随机模型，适合TDCP/SINS的动态导航数据处理。
4 结语
为了进一步提高时间差分载波相位观测量的GNSS/SINS数据处理性能，平衡滤波在算法精度和计算效率上的矛盾，本文通过对非高斯噪声影响下的随机模型进行精化处理，利用VB理论对GMM参数进行估计，构建了一种基于VB学习的混合高斯自适应滤波算法。实验结果表明，相比传统GMF算法，本文算法有效实现了混合模型参数的自适应估计，在非高斯噪声条件下有着更好的鲁棒性，同时借鉴VB学习优化了算法过程，计算负担小，具有较好的工程应用价值。
-
![]()
图 1 变分优化的高斯混合EKF滤波算法框架
Figure 1. Flowchart of VB Based on Gaussian Mixture Extended Kalman Filter
表 1 仿真1姿态误差结果（RMSE）
Table 1 Attitude Error of Simulation 1(RMSE)
方案 俯仰角/（°） 横滚角/（°） 航偏角/（°） 每历元耗时/s 1 0.582 9 0.673 8 0.872 7 0.084 2 0.314 5 0.378 4 0.734 7 0.092 3 0.227 6 0.339 9 0.583 1 0.206 4 0.164 8 0.099 8 0.472 7 0.530 5 0.197 4 0.107 3 0.474 6 0.108 
下载: 导出CSV
表 2 卫星失锁状态最大定位误差结果
Table 2 Maximum Positioning Error Under Satellite Losslock
方案 失锁20~50 s 失锁70~100 s 每历元耗时/s 水平误差/m 竖直误差/m 水平误差/m 竖直误差/m 1 60.49 2.68 60.78 2.75 0.098 2 58.72 2.01 58.33 2.03 0.114 3 56.86 1.81 56.17 1.72 0.240 4 51.40 1.43 52.93 1.47 0.608 5 51.41 1.48 52.83 1.46 0.124
下载: 导出CSV
表 3 惯性器件精度
Table 3 Inertia Device Accuracy
惯性器件 刻度系数 零偏 随机游走 MEMS 陀螺仪 0.2% 2 (°)/$\sqrt {\rm{h}} $ 0.1 (°)/s IMU 加速度计 0.3% 36 μg/$\sqrt {\rm{Hz}} $ 0.1 mg RLG 陀螺仪 5×10-6 0.002 (°)/$\sqrt {\rm{h}} $ 0.003 (°)/s IMU 加速度计 100×10-6 8 μg/$\sqrt {\rm{Hz}} $ 25 μg
下载: 导出CSV
表 4 动态数据姿态误差结果（RMSE）
Table 4 Attitude Error of Dynamic Data(RMSE)
方案 俯仰角/（°） 横滚角/（°） 航偏角/（°） 每历元耗时/s 1 0.425 9 0.182 7 3.686 1 0.086 5 0.389 9 0.102 1 2.589 6 0.108
下载: 导出CSV
-
[1] Han S, Wang J. Integrated GPS/INS Navigation System with Dual-Rate Kalman Filter [J]. GPS Solutions, 2012, 16(3): 389-404 doi: 10.1007/s10291-011-0240-x
[2] Zhao Y, Becker M, Becker D, et al. Improving the Performance of Tight-Coupled GPS/INS Navigation by Using Time-Differenced GPS-Carrier-Phase Measurements and Low-Cost MEMS IMU [J]. Gyroscopy and Navigation, 2015, 6(2):133-142 doi: 10.1134/S2075108715020108
[3] Chang G. Robust Kalman Filtering Based on Mahalanobis Distance as Outlier Judging Criterion [J]. Journal of Geodesy, 2014, 88(4): 391-401 doi: 10.1007/s00190-013-0690-8
[4] Alspach D L, Sorenson H W. Nonlinear Bayesian Estimation Using Gaussian Sum Approximations [J]. IEEE Transactions on Automatic Control, 1972, 17(4): 439-448 doi: 10.1109/TAC.1972.1100034
[5] 曹轶之.非高斯/非线性滤波算法研究及其在GPS动态定位中的应用[D].郑州: 信息工程大学, 2012 http://cdmd.cnki.com.cn/Article/CDMD-90005-1013161239.htm Cao Yizhi. Research of Non-Gaussian/Nonlinear Filtering Algorithms and Its Applications in GPS Kinematic Positioning [D]. Zhengzhou: Information Engineering University, 2012 http://cdmd.cnki.com.cn/Article/CDMD-90005-1013161239.htm
[6] 唐炉亮, 杨雪, 靳晨, 等.基于约束高斯混合模型的车道信息获取[J].武汉大学学报·信息科学版, 2017, 42(3): 341-347 http://ch.whu.edu.cn/CN/abstract/abstract5684.shtml Tang Luliang, Yang Xue, Jin Chen, et al. Traffic Lane Number Extraction Based on the Constrained Gaussian Mixture Model [J]. Geomatics and Information Science of Wuhan University, 2017, 42(3): 341-347 http://ch.whu.edu.cn/CN/abstract/abstract5684.shtml
[7] Ishikawa Y, Takeuchi I, Nakano R. Multi-directional Search from the Primitive Initial Point for Gaussian Mixture Estimation Using Variational Bayes Method [J]. Neural Networks, 2010, 23(3): 356-364 doi: 10.1016/j.neunet.2009.08.003
[8] Lim K L, Wang H. MAP Approximation to the Variational Bayes Gaussian Mixture Model and Application [J]. Soft Computing, 2017 (3): 1-13 http://cn.bing.com/academic/profile?id=1e46675d931c561e10ea5fa8b913bda4&encoded=0&v=paper_preview&mkt=zh-cn
[9] Vrettas M D, Cornford D, Opper M. Estimating Parameters in Stochastic Systems: A Variational Bayesian Approach [J]. Physica D Nonlinear Phenomena, 2011, 240(23): 1877-1900 doi: 10.1016/j.physd.2011.08.013
[10] Miao Z Y, Lv Y L, Xu D J, et al. Analysis of a Variational Bayesian Adaptive Cubature Kalman Filter Tracking Loop for High Dynamic Conditions [J]. GPS Solutions, 2017, 21(1): 111-122 doi: 10.1007/s10291-015-0510-0
[11] Stepanov O A, Chelpanov I B, Motorin A V. Accuracy of Sensor Bias Estimation and Its Relationship with Allan Variance [J]. Gyroscopy and Navigation, 2017, 8(1): 51-57 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=725670abb056ae1e19b5d91071445a17
[12] Niu X J, Chen Q J, Zang Q. Using Allan Variance to Analyze the Error Characteristics of GNSS Positioning [J]. GPS Solutions, 2014, 18(2): 231-242 doi: 10.1007/s10291-013-0324-x
[13] Niu X, Li Y, Zhang Q, et al. Observability Analysis of Non-holonomic Constraints for Land-Vehicle Navigation Systems [J]. Journal of Global Positioning Systems, 2012, 11(1): 80-88 doi: 10.5081/jgps
-
期刊类型引用(6)
1. 林雪原，刘丽丽，董云云，陈祥光，杨海利. 改进的GNSS/SINS组合导航系统自适应滤波算法. 武汉大学学报(信息科学版). 2023(01): 127-134 . 
百度学术
2. 王勋，崔先强，高天杭. 动力学模型自适应滤波算法研究. 武汉大学学报(信息科学版). 2023(05): 741-748 . 百度学术
3. 孙玉梅，潘新龙，王美春. 基于VB-AUKF的GNSS/SINS组合导航信息融合算法. 弹箭与制导学报. 2023(03): 1-8 . 百度学术
4. 王思琦，刘江，蔡伯根，赵阳. 基于VBGMM-DCNN的列车卫星定位欺骗干扰检测方法. 导航定位与授时. 2023(04): 58-68 . 百度学术
5. 朱璐瑛，孙炜玮，刘成铭，孙兆玮. 多传感器组合导航系统的联邦UKF算法研究. 电子测量与仪器学报. 2022(07): 91-98 . 百度学术
6. 张大为，沈勇. 基于PLC和机器视觉的O型密封圈尺寸检测系统设计. 电子设计工程. 2021(14): 163-167 . 百度学术
其他类型引用(3)
计量
- 文章访问数: 1830
- HTML全文浏览量: 174
- PDF下载量: 134
- 被引次数: 9
