1.   引言

粒子滤波是通过对后验滤波密度的量测更新给出离散状态的贝叶斯近似解。不但适用于非高斯非线性离散系统的状态估计，而且序列重要性重采样是获取样本的主要方法，通过样本与对应似然权值的加权和得到滤波状态^[1]。由于粒子滤波方法结构简单局限性小，因此已被广泛应用到各种实际问题中。

粒子滤波要求系统的模型参数或者量测噪声统计信息保持不变，但在实际应用中可能噪声参数未知且随环境的变化而变化，这时再使用传统粒子滤波算法会使滤波性能下降甚至滤波失效。因此在滤波的同时需要对噪声统计特性或模型结构参数进行更新。针对这个问题，Storvik^[2]给出未知静态参数下的状态和参数同时估计的粒子滤波算法，通过给定模型和未知噪声参数的充分统计量，利用已知的参数递推方法更新参数，得到实时状态估计。文献[3,4]又考虑了统计特性恒定和时变的高斯噪声，用共轭先验设定Gauss分层参数的分布，并进行递推更新。变分贝叶斯(Variational Bayes, VB)方法是一种从结构上近似推断的贝叶斯方法，在假定参数和隐变量独立情况下，能将联合参数估计问题分解为近似的单变量估计问题，依据后验密度和真实密度在KL (Kullback–Leibler)距离最大化准则下表示近似程度^[5]。Piche等人^[6]首次提出用VB结合学生t量测分布模型进行状态估计，从而解决观测过程中出现厚尾非高斯量测，文献[7,8]用学生t分布分别来建模量测模型和过程模型，结合VB方法得到新颖的鲁棒卡尔曼滤波。文献[9]将VB方法与粒子滤波结合用来实现对具有学生t分布的厚尾量测模型进行非线性状态估计，同文献[3,4]一样采用边缘化的方法对状态与参数分别估计，不同之处在于将学生t分布的自由度参数设为待估计参数，同尺度参数与隐变量一道，通过伽马共轭先验分布和VB方法得到参数后验分布表达式。仿真表明，由于自由度可以调节，使学生t分布厚尾形状可变，所以在相同噪声情况下估计性能优于文献[3,4]中所提的方法。粒子滤波与变分贝叶斯方法相结合在文献[10]中也给出一种形式，但它是利用变分贝叶斯在分割的空间内采样粒子，以使粒子滤波适用于高维状态空间的估计，这种形式不同于本文用VB估计参数的研究领域。

从学生$\rm t$构成形式可以发现，决定学生t分布对称轴位置的均值变量与尺度参数、自由度参数具有同等重要地位。而且在实际应用中，如导航^[11]和测距^[4]中，由于传感器偏移使噪声均值非0且随系统运行发生变化，所以对均值变量的估计也十分必要。文献[11,12]通过极大后验和VB方法对线性状态模型的量测噪声均值进行了估计。但对于非线性状态空间模型，使用VB方法对噪声均值的估计还没有涉及。从贝叶斯推断的角度，均值作为变量与尺度参数相关联，因此通过VB推导的均值后验分布表达式不同于以往没有均值参与的，如文献[9]所给形式。在得到量测噪声参数后，对量测噪声时变的噪声同时关联的粒子滤波问题进行探究。尽管文献[13,14]涉及了粒子滤波的噪声关联问题，但是它们不能解决量测噪声时变或有野值的情况。

基于以上讨论，本文用学生t分布对非高斯量测噪声进行建模，同时考虑均值参数的影响，将均值与尺度矩阵联合共轭先验为高斯-伽马分布，利用变分贝叶斯方法结合粒子滤波对时变量测噪声下非线性状态进行估计得到MPF-VBM算法，然后通过获得参数估计研究过程噪声依赖时变量测噪声的粒子滤波算法MPF-VBM-COR。文章的其余部分组织如下，问题公式化在第2节给出，在学生t量测分布下的参数和状态估计算法及噪声关联算法在第3节给出，仿真演示的例子在第4节给出，第5节给出结论。

2.   问题公式化

首先考虑如式(1)、式(2)所示的非线性状态空间模型

其中，${{\text{x}}_k} \in {\mathbb{R}^n}$状态向量，量测${{\text{z}}_k} \in {\mathbb{R}^m}$, $f\left( \cdot \right)$和$h\left( \cdot \right)$都是已知非线性函数，过程噪声${{\text{w}}_k} \in {\mathbb{R}^n}$假定服从0均值高斯分布，具有方差${\rm{E}}\left[ {{{\text{w}}_k}{\text{w}}_k^{\rm{T}}} \right] = {{\text{Q}}_k}$，厚尾观测噪声${{\text{v}}_k}$服从对称学生t分布，分布表达式为

其中，$S\left( \cdot \right)$表示学生t分布，$\varGamma \left( \cdot \right)$表示伽马函数，${\text{μ}} ,{\text{Λ}}$和$\nu$分别表示均值，尺度矩阵和自由度参数。应该指出的是学生t分布是一种厚尾的高斯分布，当$\nu \to \infty$时就成为具有相同$\text{μ}$的对应高斯分布。从式(3)可见，即使对于同一个数据点${\text{x}}$，均值的不同也会使密度值不同。所以准确知道均值变量十分重要。为进行利用VB的参数估计，不失一般性，假定${{\text{v}}_k}$的尺度矩阵${\text{Λ} _k}$为对角矩阵，即${\text{Λ} _k} = {\rm diag}({\varLambda _{k,1}},$${\varLambda _{k,{\rm{2}}}},·\!·\!·,{\varLambda _{k,m}})$。对于式(1)、式 (2)所给状态空间模型，依据学生t的量测噪声分布特点，时刻$k$的量测似然可表示为具有隐变量${u_k}$的如式(4)、式(5)的高斯分层形式

其中，${u_k}$是引入的隐变量，$G\left( { \cdot \left| {\alpha ,\beta } \right.} \right)$表示服从参数为$\alpha ,\beta$的伽马分布。由文献[15]所示，均值和尺度矩阵的联合共轭先验服从高斯-伽马分布，并且可分解为各自单独分布的乘积

其中，${\text{η} _{k\left| {k - 1} \right.,j}}$和${\left( {{\beta _{k\left| {k - 1} \right.}}{\varLambda _{k,j}}} \right)^{ - 1}}$分别为${\text{μ} _{k,j}}$的均值和方差，${\beta _{k\left| {k - 1} \right.}}$为精度系数。

自由度${\nu _k}$的共轭先验服从伽马分布

为了更清晰地说明各变量之间的关系，用图1所示的有向图对$k$时刻状态，量测和超参数的贝叶斯分层模型进行表示：

图 1  贝叶斯推断图模型

下载: 全尺寸图片 幻灯片

从图1中可见，相比于在文献[9]，两个额外的超参数${\beta _k}$和${\text{η} _k}$被加入进来，并且${\text{μ} _k},{\text{Λ} _k}$之间具有关联性，相比于文献[6]，自由度变量具有超参数${a_k},{b_k}$。

3.   主要结果

本节主要应用VB对式(1)—式(7)的中参数的变分后验进行推导。VB方法的原理及公式详见文献[7,15]。

首先，状态、量测和参数的联合分布可分解为

其中，${\text{x}}_k^{\left( i \right)}$表示$k$时刻第$i$个粒子，令${\text{Θ} _k} = \left\{ {{\text{μ} _k},{\text{Λ} _k},{\nu _k}} \right\}$表示参数集，利用边缘化思想可分别进行状态与参数的估计，在时刻$k$，状态参数联合后验分布为

而且

其中，$\delta \left( \cdot \right)$表示Dirac-Delta函数。$\omega _k^{\left( i \right)}$是对应第$i$个粒子的重要性权值，式(10)中的似然$p\left( {{{\text{z}}_k}\left| {{\text{x}}_k^{\left( i \right)},\text{Θ} _{k\left| {k - 1} \right.}^{\left( i \right)}} \right.} \right)$服从学生t分布，可见它不但依赖粒子样本值，而且依赖于$k - 1$时刻参数一步预测的样本值。下面给出在量测更新后，$\text{Θ} _k^{\left( i \right)}$依据VB的具体更新过程。

首先，对于均值变量${\text{μ} _k}$和尺度矩阵${\text{Λ} _k}$，通过基本贝叶斯公式$p(x|y) = {{p(x,y)} / {p(y)}}$, $\left( {{\text{μ} _k},{\text{Λ} _k}} \right)$联合变分后验的对数函数表示为

其中，${{\rm{E}}_x}\left[ \cdot \right]$表示括号中的概率分布对变量$x$求期望，${C_{\text{μ}} ,{\text{Λ}}}$是独立于$\text{μ}$和$\text{Λ}$的常数项，${\left( \cdot \right)_j}$表示括号中向量的第$j$个分量。${\text{μ} _k}$边缘后验对数函数为

由于式(12)是${\text{μ} _k}$的2次式，所以变分后验$q\left( {{\text{μ} _k}\left| {{\text{Λ} _k}} \right.} \right)$服从正态分布${N}\left( {{{\text{μ}} _k}\left| {{{\text{η}} _{k\left| k \right.}},{{\left( {{\beta _{k\left| k \right.}}{{\text{Λ}} _{k\left| k \right.}}} \right)}^{ - 1}}} \right.} \right)$，其中参数的递推公式为

${\text{Λ} _k}$的变分后验分布对数式由式(11)减式(12)得到

将式(13)，式(14)代入式(15)，通过表达式化简和共轭先验的概念，得到${\varLambda _{k,j}}$变分后验服从伽马分布$q\left( {{\varLambda _{k,j}}} \right) = {{G}}\left( {{\varLambda _{k,j}}\left| {{c_{k\left| k \right.,j}},{d_{k\left| k \right.,j}}} \right.} \right)$，其中参数

其中，${\left[ {\text{A}} \right]_{jj}}$表示矩阵${\text{A}}$主对角线第$jj$个元素。

隐变量${u_k}$变分后验$q\left( {{u_k}} \right)$的对数函数表示为

对式(18)指数化运算，$q\left( {{u_k}} \right)$服从伽马分布$q\left( {{u_k}} \right) = G({u_k}\left| {\nu _k^1,\nu _k^2} \right.)$，其中

这里，在${{\rm{E}}_{{\text{μ} _k} {\text{Λ} _k}}}\left[ {{\left( {{{\text{z}}_k} - h\left( {{\text{x}}_k^{\left( i \right)}} \right) - {\text{μ} _k}} \right)}^{\rm{T}}}{\text{Λ} _k}\left( {{\text{z}}_k} - \right.\right.$$\left. h\left( {{\text{x}}_k^{\left( i \right)}} \right)- {\text{μ} _k} ) \right]$推导中，使用了${\text{μ} _k}$的如式(20)的2阶矩

自由度${\nu _k}$的近似后验为

使用Stirling的近似$\ln \varGamma \left(\dfrac{{{\nu _k}}}{2}\right) \approx \dfrac{{{\nu _k} - 1}}{2}\ln \dfrac{{{\nu _k}}}{2}$$- \dfrac{{{\nu _k}}}{2}$, $q\left( {{\nu _k}} \right)$为伽马分布$q\left( {{\nu _k}} \right) = G\left( {{\nu _k}\left| {{a_{k\left| k \right.}},{b_{k\left| k \right.}}} \right.} \right)$，其中

在以上每一个参数推导过程中，都需要对其余参数进行期望操作，根据分布期望公式，参数期望如式(24)给出

其中，$\varPsi ( \cdot )$为di-gamma函数。

尽管以上参数更新相互耦合，但这一问题解决可通过对任一参数初始化实现，然后依次迭代更新。同时可见，如果均值变量为0，涉及均值的超参数${\text{η} _k}$和${\beta _k}$的递推式(13)、式(14)以及在式(17)和式(19)含有的${\text{η} _k}$和${\beta _k}$都消失了，算法退化为文献[9]所提算法。因此所提算法可视为文献[9]算法的丰富和扩展。

噪声可能为时变的情形，因此引入遗忘因子$\rho \in \left( {0,1} \right)$标识噪声波动程度，噪声参数时间更新模型为

表1给出所提算法1的流程。

表 1  基于VB带有噪声均值估计的边缘粒子滤波(MPF-VBM)算法

从分布${ {\text{N} } }\left( { { {\text{x} }_0}\left| { { {\text{m} }_{0\left| 0 \right.} },{ {\text{P} }_{0\left| 0 \right.} } } \right.} \right)$采样粒子${\text{x}}_0^{\left( i \right)}$$i = 1,2, ·\!·\!· ,N$，并且设置权值$\omega _0^{\left( i \right)} = {1 / N}$；初始化超参数$a_0^{\left( i \right)},\;\;b_0^{\left( i \right)},\;\;c_0^{\left( i \right)},\;\;d_0^{\left( i \right)},\nu _0^{\left( i \right)},\text{η} _0^{\left( i \right)}$和$\beta _0^{\left( i \right)}$；计算初 始参数${\text{μ} _0}$, ${\text{Λ} _0}$和${\nu _0}$期望。对时刻$k = 1,2, ·\!·\!· ,K$

对每一个粒子$i = 1,2, ·\!·\!· ,N$

(1) 使用式(25)做噪声参数的时间更新； 　(2) 从状态传递方程$p\left({\text{x} }_k^{\left( i \right)}\left| {\text{x} }_{k - 1}^{\left( i \right)} \right.\right)$做粒子的一步预测；

(3) 通过新量测${{\text{z}}_k}$用式(15)更新重要性权值；

(4) 必要的话，粒子重采样；

(5) 使用式(13)，式(14)，式(16)，式(17)，式(18)，式(19)，式(22)，式(23)，并利用重采样粒子做噪声参数后验更新， $\varTheta _{k\left| k \right.}^{\left( i \right)} = T\left( {\varTheta _{k\left| {k - 1} \right.}^{\left( i \right)},x_k^{\left( i \right)},{z_k}} \right)$；其中$T\left( \cdot \right)$代表参数充分统计量；

得出当前的噪声参数期望值及对应状态值；

在$k = = K$前，进行下一次循环。

下载: 导出CSV 

| 显示表格

依据表1，接着考虑过程噪声依赖时变量测噪声的粒子滤波算法。将学生t量测似然式(4)改写成高斯分层形式

其中，正态分布的方差表达式为${{{{\text{R}}_k}} / {{u_k}}}$, ${{\text{R}}_k}$表示式(4)中${\left( {{\text{Λ} _k}} \right)^{ - 1}}$。假定过程噪声依赖量测噪声的关联矩阵为${{\text{S}}_k}$，即

同时假定${{\text{w}}_k}$, ${{\text{v}}_k}$各自单独序列仍然是相互独立的。下面利用解关联框架进行噪声解关联操作，根据式(2), ${{\text{z}}_k} - h\left( {{{\text{x}}_k}} \right) - {{\text{v}}_k} = 0$，那么状态传递方程式(1)修改为

这里，${f^*}\left( {{{\text{x}}_k}} \right) = f\left( {{{\text{x}}_k}} \right) - {{\text{J}}_k}h\left( {{{\text{x}}_k}} \right)$, ${\text{w}}_k^* = {{\text{w}}_k}-$${{\text{J}}_k}{{\text{v}}_k}$, ${{\text{J}}_k}$是要确定的辅助参数。${\text{w}}_k^*$可被看成新的过程噪声，由于是高斯分布的线性组合，所以仍然服从高斯分布，且它与量测噪声${{\text{v}}_k}$的协方差为

令协方差${\rm{E}}\left[ {{\text{w}}_k^*{\text{v}}_k^{\rm{T}}} \right] = 0$，得到辅助矩阵${{\text{J}}_k}$为

式(29)，式(2)连同式(31)一起给出噪声同时关联的等价状态空间模型。通过式(31)，等价传递模型的高斯过程噪声${\text{w}}_k^*$的均值与方差为

因此，由式(1)，式(2)组成具有噪声关联的系统转化为由式(29)，式(2)组成不含噪声关联的等价系统滤波问题。算法流程与表1算法相同，就是在第(2)步中增加对${\text{w}}_k^*$噪声参数的计算式(32)和式(33)，状态转移模型由${p_{{{\text{w}}_k}}}\left( {{\text{x}}_k^{\left( i \right)}\left| {{\text{x}}_{k - 1}^{\left( i \right)}} \right.} \right)$替换为${p_{{\text{w}}_k^*}}\left( {\text{x}}_k^{\left( i \right)}\left| \right.\right.$$\left. {{\text{x}}_{k - 1}^{\left( i \right)},{{\text{z}}_{k - 1}}} \right)$，这样得到时变量测噪声下噪声同时关联的粒子滤波算法，记作MPF-VB-COR-1(Marginalized Particle Filter with VB for Correlation-1)。

4.   仿真验证

这里考虑在文献[4,9]广泛使用的如下非线性模型来验证上面所提两种算法的性能

算法1　初始状态${{\text{x}}_0} \sim {N}\left( {0,5} \right)$，过程噪声分布$p\left( {{w_k}} \right)=N\left( {{w_k}\left| {0,5} \right.} \right)$，量测噪声${v_k}$服从学生t分布，初始超参数为：${a_0} = {b_0} = 0.12$, ${c_0} = 2$, ${d_0} = 5$, ${\eta _{{0}}} = 1$和${\beta _{{0}}} = 2$，遗忘因子选取根据文献[9]的遗忘因子与超参数收敛时间的比较说明，设为$\rho = 1 - \exp \left( { - 4} \right)$。仿真时间长度设为1000。根据文献[9]的粒子数与状态均方根误差平均值(ARMSE)的关系，考虑程序运行的时间因素，粒子数设为100。Monte Carlo仿真运行30次，所提的MPF-VBM算法与另外两种算法MPF-VBS (Student)^[9]和MPF-CP (Conjugate Prior)^[3]在相同仿真状况下进行均方根误差(RMSE)^[8]及(ARMSE)^[9]比较。仿真设置3种不同情况噪声，噪声及对应的状态ARMSE及时间如表2所示。

表 2  对应3种噪声的3种算法的均方根误差平均值(ARMSE)

量测噪声	MPF-VBM	MPF-VBS	MPF-CP
N(0,1)+20%U(–20, 20)野值	5.5673	5.5468	8.5571
N(6, 1)无野值	4.6170	8.1141	8.8132
N(6,5)+20%U(20, 60)野值	6.3353	8.6288	8.7588
运行时间(s)	0.01067	0.007187	0.005587

下载: 导出CSV 

| 显示表格

${\rm U}\left( \cdot \right)$表示均匀分布，野值以随机形式进入噪声分布。第1种标准高斯噪声加入20%服从均匀分布野值，噪声均值为0；第2种高斯噪声无野值加入，均值为6；第3种高斯噪声均值为6有20% U(20, 60)的野值加入。MPF-VBM与MPF-CP量测噪声均值估计比较图如图2所示。

图 2  算法噪声均值估计比较图

下载: 全尺寸图片 幻灯片

从图2可见，本文MPF-VBM算法对3种噪声情况的均值都具有稳健而准确的估计能力，相比较而言，MPF-CP的均值估计易受到野值干扰，跟踪真值时的震荡性较强。但从图2(b)估计图所示，在没有野值干扰时，MPF-CP均值估计具有收敛速度快优点。

图3—图5给出了3种算法对应3种噪声的RMSE的比较图，为了清晰起见，这里只列出了300～330 s时间范围内的RMSE的比较曲线。从图3—图5可见，所提MPF-VBM算法对3种噪声都始终处于估计性能最优越的算法行列。3种算法对应3种噪声ARMSE如表2所示，结果验证了RMSE图的显示。在运行时间上，MPF-VBM近似为MPF-VBS时间的1.5倍，为MPF-CP时间的2倍。量级为ms级。

图 3  3种算法对应第1种噪声RMSE比较图

下载: 全尺寸图片 幻灯片

图 5  3种算法对应第3种噪声RMSE比较图

下载: 全尺寸图片 幻灯片

图 4  3种算法对应第2种噪声RMSE比较图

下载: 全尺寸图片 幻灯片

算法2　所提噪声关联的MPF-VB-COR与算法MPF-VBM及文献[13]中的算法PF-COR在相同噪声关联情况下进行性能比较。Monte Carlo仿真次数为30次。过程噪声与量测噪声超参数设置与算法1相同。设过程噪声${{\text{w}}_k}$量测噪声${{\text{v}}_k}$有关系式${{\text{w}}_k} = {{\text{T}}_k}{{\text{v}}_k}$，其中${{\text{T}}_k}$为常矩阵，此处设为$\left[ {{{\text{T}}_k}} \right] = 0.8$，得到关联矩阵${{\text{S}}_k} \!=\! {\rm{E}}\left[ {{{\text{w}}_k}{\text{v}}_k^{\rm{T}}} \right] \!=\! {{\text{T}}_k}{\rm{E}}\left[ {{{\text{v}}_k}{\text{v}}_k^{\rm{T}}} \right]$$= {{{{\text{T}}_k}{\rm E}\left[ {{{\text{R}}_k}} \right]} / {{\rm E}\left[ {{u_k}} \right]}}$。为了更清晰表示算法性能，图6—图8给出单次Monte Carlo状态的ARMSE的曲线比较图，表3给出设置的量测噪声及30次Monte Carlo的ARMSE。其中，第2种噪声方差$R$如下变化，当$t \in \left[ {1\;{\rm s}:300\;{\rm s}} \right]$, $R = 8$；当$t \in \left[ {301\;{\rm s}:700\;{\rm s}} \right]$, $R = 20$；当$t \in \left[ {701\;{\rm s}:}\right.$$\left.{1000\;{\rm s}} \right]$, $R = 2$。

图 6  3种算法对第1种关联噪声的ARMSE比较图

下载: 全尺寸图片 幻灯片

图 8  3种算法对第3种关联噪声的状态ARMSE比较图

下载: 全尺寸图片 幻灯片

表 3  对应3种噪声的3种算法的均方根误差平均值(ARMSE)

量测噪声	MPF-VBM	MPF-VB-COR	PF-COR
${{N} }\left( {0,\;{5 / 2} } \right)$	3.8662	3.4025	3.1384
${{N} }\left( {0,\;R} \right)$ $R$时变	5.5178	4.5536	5.0392
${{N} }\left( {0,\;{5 / 2} } \right)$+20% Unif (–5, 5)野值	8.5892	8.7574	11.0393

下载: 导出CSV 

| 显示表格

图 7  3种算法对第2种关联噪声的状态ARMSE比较图

下载: 全尺寸图片 幻灯片

从图6—图8及表3可见，当第1种量测噪声是参数确定的高斯噪声时，PF-COR的估计性能最好，本文的MPF-VB-COR-1与PF-COR估计性能相近，两者都好于没有考虑关联性的MPF-VBM，从ARMSE曲线看，含有变分估计方法的估计一致性更好；第2种量测噪声是方差时变的高斯噪声，MPF-VB-COR的估计性能好于PF-COR，两者的性能都优于MPF-VBM，仍然说明考虑了噪声关联情况粒子滤波算法此时更好。在第3种具有野值的高斯噪声，MPF-VB-COR与MPF-VBM的性能基本相同，好于PF-COR，说明此时具有抑制野值干扰的滤波算法具有更好性能，而所提MPF-VB-COR也是具有这一能力的。

5.   结论

本文讨论了具有学生t量测噪声的鲁棒粒子滤波算法。量测噪声被设计成服从学生t分布，利用变分贝叶斯方法对t分布的噪声包括均值在内的全部参数进行了估计。先给出了参数的共轭先验，然后推导出参数变分后验表达式。通过边缘化思想给出估计参数值对应的状态值。将此算法与已有学生t量测噪声的粒子滤波方法比较，对于给出的噪声情况，所提粒子滤波展现出其有效性和优越的性能。接着基于所提鲁棒粒子滤波框架，探究了量测噪声时变的噪声关联的粒子滤波，仿真表明当量测噪声时变或有野值存在时，所提噪声关联的粒子滤波算法优于已有噪声关联粒子滤波算法。所提算法运行时间一般为${10^{ - 2}}$ s的量级。