带不确定混合噪声系统的变分贝叶斯期望最大滤波算法

马天力，张 扬，陈超波

（1.西安工业大学电子信息工程学院，西安 710021；2.陕西省自主系统与智能控制国际联合研究中心，西安 710021）

状态估计的目的是从受到噪声污染或者干扰的观测信号中，消除噪声以及干扰的影响，准确获得最优的系统状态[1]。典型的滤波算法是采用卡尔曼滤波器(Kalman filter, KF)，其假设量测噪声已知且服从单一高斯分布。但在实际过程中，量测噪声不仅包含概率特性未知的随机噪声[3]，还可能存在边界未知的有界噪声[4,5]。例如GPS/INS组合导航系统中，GPS解算误差属于非高斯随机噪声，而INS误差是已知其边界的有界噪声。

对于存在未知随机噪声系统的状态估计，主要采用自适应滤波技术[6]，该方法是一种基于贝叶斯理论的递归估计处理策略。需要已知随机变量的概率统计特性，并在递归计算噪声参数的同时进行状态估计。其包括状态增广法[7]、交互式多模型理论[8]以及序贯蒙特卡洛方法[9]。另一类方法是基于期望最大（Expectation Maximization, EM）理论的批处理技术[10]，从包含有隐变量的量测数据中计算最大似然参数。对于包含未知但有界（Unknown-But-Bounded,UBB）噪声系统的状态估计，由于噪声统计特性未知，仅已知噪声所在范围的边界信息，通常采用集员滤波（Set membership filter, SMF）[11]进行求解，与概率方法不同，其以包含系统真实状态的外定界椭球集合为基础，估计得到的是状态变量的可行集而不是后验概率密度，该可行解集由与系统模型、初始状态和UBB噪声假设相一致的所有状态点构成。

上述方法仅针对存在某一种不确定噪声进行计算，对于两种或多种未知噪声同时存在的混合不确定噪声系统，一方面非高斯噪声存在致使集员滤波器选取噪声边界困难，过小的噪声边界造成估计精度下降；另一方面，因集员噪声统计特性为均匀分布，采用相应的贝叶斯滤波算法则无法获得其解析解。解决混合噪声条件下的状态估计问题的核心是如何将两类不确定噪声进行统一描述。现有方法主要包括三类，一种是基于序贯蒙特卡洛采样理论。Zou[12]提出箱粒子滤波算法（Box-Particle filter, BPF），该方法在粒子滤波基础上结合区间分析理论对非高斯噪声干扰下的状态进行估计，但受到箱粒子数目的影响，越大的箱粒子数目，算法运行时间越长。一种是基于随机集合理论[13-15]，将加性集员不确定纳入广义似然函数中。Handbeck[16]利用融合随机集合提出统计与集合理论信息滤波器（Statistical and Set theory Information filter,SSI），传统单一的滤波方法得到的值是SSI估计器的边界情况，即当随机误差消失，其收敛于集合估计；当有界误差为零，其收敛于贝叶斯估计。对于混合噪声，其估计结果是统计状态下具有不确定边界的解集。胡斌[17]将随机观测集看作量化单元，根据量测集识别理论对于量化单元进行计算。另一种策略是非精确概率描述[18,19]，其用概率密度函数的闭型凸集对混合噪声进行构建。Noack[20]将混合不确定噪声表示为概率密度的集合，将系统状态转移与似然密度的信度需求分别扩展到状态转移与似然密度凸集的形式，提出信度状态滤波器（Credal state filter, CSF）。Klumpp[21]对基于随机集合的SSI滤波器和基于集合概率密度的CSF进行了比较，CSF估计结果相比SSI更保守。此外，Ginger[22]结合序贯蒙特卡洛理论提出非精确采样概念，运用混合非精确采样集合对非精确量测进行建模，并在贝叶斯框架下通过集员理论中有界支撑传播非精确信息。但该方法对参数具有较强的敏感度。Henningsson[23]根据控制中观测器的思想，用椭球包含混合噪声中有界集合部分的估计误差，并通过线性矩阵不等式得到最优滤波增益，该算法的性能受到一个权系数的影响，而该系数的取值取决于实际系统总体噪声中随机不确定和有界不确定的某种权重关系，但一般这种权重关系并不明确。江涛[24]提出混合噪声的联合滤波算法（Combined filter, CF），其将统计特性未知但有界的噪声引入到卡尔曼滤波模型中，得到一组包含集合运算的卡尔曼滤波方程。其中UBB噪声应用SMF的思想进行处理，而随机噪声应用卡尔曼滤波的思想进行处理。该算法问题在于统计噪声为参数已知的高斯分布。

本文针对未知随机不确定噪声和有界不确定噪声共同存在下的状态估计问题，提出基于混合噪声框架的状态估计方法，将随机和有界不确定噪声运用未知参数的混合统计噪声模型进行描述，将混合噪声条件下的状态估计问题描述为未知参数联合统计噪声条件下的参数与状态的估计问题，并设计变分期望最大算法求解其状态及噪声参数。

1 系统模型

假设含有随机噪声和有界噪声的不确定系统的状态方程为

其中xk为k时刻系统状态向量，xk⊂Rn；wk为零均值，方差为Qk的高斯过程噪声，其中Qk为非负正定协方差矩阵。zk为k时刻的系统量测，vk+ek为混合测量噪声，其中vk为零均值非高斯随机量测噪声，ek为UBB噪声，ek∈Ek⊂Rm，Ek为边界未知的有界集合。假设在上述模型中，ek= 0，即量测噪声为典型的随机模型，该模型可以用卡尔曼滤波器进行求解。

随机、集员和混合不确定统计特性模型的概率密度函数表示如图1所示。由于有界和随机并没有严格的划分，两者可以近似转换。对于有界噪声，若利用统计特性进行表述，则高斯白噪声方差σ为有界噪声边界的1/3倍[24]。

图1 随机与集员混合模型Fig.1 Combined stochastic and set-membership model

基于该特性，有界不确定噪声可近似描述为参数未知的高斯噪声，并构建两种不确定噪声的混合高斯模型。

由于高斯混合项中含有隐变量、且每个高斯分布的均值、方差和权重系数均未知。若采用贝叶斯方法估计，则需考虑参数先验分布。对于混合权重系数lα，其先验为Dirichlet分布[10]。

其中λ0是浓度参数。噪声高斯项均值和方差的先验服从Gaussian和Wishart分布。

其中υ0为Wishart分布的自由度，ϖ0为精度矩阵。

令 Θ =[α,μ, Σ,S]，本文的目的是在已知观测集合Z={z1,z2,… ,zN}的情况下，求解最优系统状态X={x0,x1,…xn}和参数集合Θ的最大后验分布。

由状态空间模型可知：

系统的似然函数为：

2 变分EM滤波算法

2.1 变分贝叶斯学习

对于未知参数Θ的估计，若系统量测Z已知，则根据贝叶斯规则：

一般情况下，边缘概率密度函数p(Z)的积分形式复杂，难以求得其解析解，导致难以得到参数后验概率密度函数的解析表达式p(X, Θ|Z)。针对这一问题，根据先验分布满足共轭指数分布族，变分贝叶斯（Variation Bayesian, VB）理论提出构建一个新的分布q(X,Θ)去逼近真实后验分布p(X, Θ|Z)，无需给定其概率密度函数形式，表达形式与p(X, Θ|Z)一致，区别在于参数不同。故边缘似然函数的log形式表示为：

其中，自由能量函数 F (q(X,Θ)是logp(Z)的下界函数[6]。K L(q(X, Θ) ||p(X, Θ |Z))为q(X,Θ)和p(X,Θ|Z)之间的Kullback-Leibler散度，其用来衡量两个概率分布之间的差异。因KL散度非负，由式(11)可知，只需使KL散度最小，求最大化自由能量函数max F (q(X,Θ) )，即可获得最佳近似分布。为了使模型参数可求，VB采用均值域逼近理论将多变量联合概率分布近似逼近为各个变量边缘分布的乘积[25]。q(X,Θ)的表达式如下：

对于每一个参数的近似分布q( Θi)对KL散度

求偏导，即可求得q(Θi)的通解表达式

2.2 变分EM估计策略

由于在联合后验概率分布计算过程中，所求参数相互耦合。故本文将EM理论中参数迭代计算后验参数思想引入VB理论，设计参数耦合更新策略，通过反复迭代求解最大期望。在变分贝叶斯期望（Variational Bayesian expectation, VBE）步骤中，利用前向-后向递推算法计算隐变量的后验，并在变分贝叶斯最大化（Variational Bayesian maximization, VBM）步骤中对模型超参数进行更新。

A 变分贝叶斯期望步骤

根据文献[25]，利用前向递推算法求解隐变量的充分统计量，通过输入隐马尔可夫模型参数矩阵以及观测序列计算该序列发生的概率。假设ak(xk) =p(xk|z1:k)为隐变量的后验概率密度，结合公式(14)可知状态xk的边缘后验为：

其中ζk=p(zk|zk-1)为归一化因子。a(xk)中关于xk-1边缘概率密度函数满足高斯分布形式，其均值为，方差为经过简化：

将上式带入式(15)中，并对式(15)进行边缘化处理，使其只含有隐状态xk项，便能够计算高斯分布的均值θk和方差Ωk，如下所示

与前项递推滤波过程不同，后向滤波根据观测量zk+1:N估计现有状态xk的后验，采用并行递推的形式实现。假设bk(xk) =p(zk+1:N|xk)，xk满足终止条件bN(xN) = 1。则：

与ak(xk)类似，b(xk)同样满足高斯分布的形式，则：

近似分布q(xk)的后验概率密度计算采用前向-后向联合滤波，如下所示：

则系统状态后验的充分统计量为：

B 变分贝叶斯最大化步骤

一旦获得隐变量的最大后验分布，在VBM步骤中需要计算系统模型的超参数。由于参数变量空间维度高导致其计算困难，所提算法利用近似边缘分布去分别逼近其统计特性。根据共轭指数分布的性质[26]，高斯混合模型参数的近似边缘分布定义为：

对于高斯混合模型系数αl的推理，假设其变分后验分布与先验分布形式相同，均满足Dirichlet分布，根据式(13)(14)，对联合概率密度函数的近似分布取对数后为：

由于混合系数的后验概率密度形式与先验相同，经过变分推断：

同理，对于混合噪声高斯项均值μl和方差Σl的联合概率密度函数近似分布取对数：

其中Tr(·)表示矩阵的迹，根据式(6)，q(μl, Σl)的后验分布为高斯威沙特分布的形式。首先考虑均值μl：

上式展开后，合并关于μl的同类项，得到新的高斯分布q(μl|Σl) =N(μl|ml,(βlΣl)-1)，其中：

与混合系数的计算类似，上式为递推形式。由式(6)可知，lμ和Σl相互耦合，无法利用式(30)直接进行求解。因此，采用两者联合概率密度函数减去lμ的概率密度函数。

对上式中迹矩阵进行合并简化，其中高斯混合项方差的后验近似分布为，则：

至此，式(24)-(36)为VBEM算法的状态及参数更新方程。式(24)(25)中状态均值与方差的计算与超参数有关，式(33)(35)中超参数的更新与状态均值、方差同样相关，相互耦合。在VBEM中，需要迭代求解状态与参数，直至达到收敛。该方法与EM算法类似，其均是寻求能量函数关于隐变量分布最大值的方法。不同之处在于，EM算法在期望最大化步骤中，并未赋予参数任何先验信息，而是作为未知随机变量进行计算。而VBEM算法考虑参数近似后验分布q(Θ)，对参数的超参数进行更新计算。若将所求参数分布假设为Dirichlet函数，则VBEM对概率模型参数进行点估计求解，所提出的VBEM算法简化为EM算法。

3 仿真验证

为了验证所提算法的有效性，分别采用卡尔曼滤波器（KF），联合滤波（CF），箱粒子滤波算法（BPF）以及VBEM算法对GPS/INS松耦合系统开展仿真实验。取东北天地理坐标系为导航坐标系，以INS误差作为状态变量，GPS和INS的位置差作为观测变量，系统状态方程为：

其中x=[φE,φN,φU,δVE,δVN,δVU,δL,δλ,δh,εE,εN,εU,∇E,∇N,∇U]T，φE、φN和φU为姿态误差角；δVE、δVN和δVU为东北天向的速度误差；δL、δλ和δh为经纬度及高度误差，εE、εN与εU为陀螺的一阶马尔科夫漂移误差；∇E、∇N与 ∇U为加速度计零点漂移误差。F为系统状态矩阵。

系统量测方程为：

其中z=[VEG-VEI,VNG-VNI,VUG-VUI,LG-LI,λG-λI,hG-hI]T；LG、LI、λG、λI、hG、hI分别为GPS和INS观测的位置信息；VEG、VEI、VNG、VNI、VUG、VUI分别为GPS和INS获得的速度分量，Hk=[03×3,I6×6,03×3]为观测矩阵。

模拟车辆在二维空间运动，其行驶轨迹如图2所示。车辆初始位置为北纬34.246 °，东经108.909 °，初始速度为10 m/s；初始姿态为方位角0 °，俯仰角为0 °，横滚角为0 °。仿真时间为600 s。vk为服从Student-T分布的随机噪声序列，均值为0，自由度为1，协方差矩阵δL=δλ= 10-5(°)。ek=[ek,L,ek,λ]T为 量 测 UBB噪 声 ，ek,L,ek,λ⊂[-1 ,1]× 1 0-4(°)。

图2 车辆真实轨迹Fig.2 Trajectory of the vehicle

图3为四种算法100次蒙特卡洛仿真实验下的经纬度的误差结果。从图中可以看出，对于同时存在未知参数统计噪声和未知边界的有界噪声的组合导航系统，变分EM算法的经纬度均方根误差明显低于其他三种算法。主要原因是VBEM将混合噪声建模为未知参数的高斯混合模型，进行联合估计。而KF算法仅针对量测噪声参数已知的高斯分布模型。因此，结果存在较大偏差，估计性能不足。另外，CF算法与BPF算法均是采用区间分析策略，该策略需要已知噪声的上、下界。由于噪声模型有偏，不准确的噪声边界设定造成估计结果出现误差。表1为四种算法的经纬度均方根误差(RMSE)统计表。

图3 不同算法的经纬度误差Fig.3 The latitude and longitude errors from different method

表1 四种算法的经纬度误差统计表Tab.1 Statistical table of latitude and longitude error

从表中可以看出，VBEM算法相比于KF，CF和BPF在纬度的平均均方根误差分别下降了83.32%、80.81%和79.84%。在经度方面，VBEM算法相比于KF，CF和BPF平均均方根误差分别下降了61.28%、60.97%、56.78%。

本文对算法的运算时间进行了分析，表2为四种算法的相对运算仿真时间对比，令卡尔曼滤波算法运算时间为单位1。由表2可知，本文所提算法虽然相比于KF，CF具有较长的运算时间，但明显低于基于蒙特卡洛理论的BPF算法。

表2 四种算法的相对运行时间Tab.2 The running time of different method

为了验证VBEM算法中相关参数对系统的影响，分别从初始量测噪声均值以及精度矩阵对系统影响角度进行了分析。图4为不同初始噪声均值mL，mλ下的经纬度误差。从图中可以看出，在噪声初始均值不相同的情况下，采用变分EM算法均可使系统收敛。初始噪声均值越贴近真实值，其经纬度误差越小。图5为不同初始精度矩阵Λ下的滤波结果。可以看出，初始噪声方差越接近真实值，经纬度的均方根误差越小，主要是因为VBEM参数递归中需要引入初始噪声参数进行计算，较大误差的先验信息会对滤波系统造成干扰。

图4 不同初始噪声均值下的经纬度误差Fig.4 The latitude and longitude errors from different initial noise means

图5 不同初始精度矩阵下的经纬度误差Fig.5 The latitude and longitude errors from different initial precision

为了说明所提算法中递归循环次数对系统的影响，进行了相关验证。采用不同迭代次数分别计算整个滤波过程的均方根误差。从图6中可以看出，所提算法的收敛效果受到迭代次数的影响。迭代次数越低，相对应滤波算法的误差越大。迭代次数为20时，在递归计算过程中自由能量函数并没有达到最大，造成估计结果有偏。

图6 不同迭代次数下的经纬度误差Fig.6 The latitude and longitude errors from different iteration times

3 结 论

针对混合不确定量测噪声条件下的状态估计问题，提出基于变分贝叶斯期望最大理论的鲁棒滤波算法。由于随机噪声统计特性与有界噪声边界参数均未知，若仅采用基于统计信息或基于区间分析的策略，易造成误差偏移。因此，本文将随机不确定噪声与未知有界噪声运用未知参数的高斯混合模型进行统一表示，通过变分期望算法对模型隐变量进行推理，同时运用变分最大化算法更新模型参数。最终获得状态估计结果。仿真试验结果表明，在含有混合不确定噪声的情况下，所提算法具有较高的估计精度。鉴于变分EM具有较强的自适应性，未来考虑非线性系统模型以及具有复杂时变特性的随机及有界噪声模型。