一种基于Gibbs采样器的鲁棒滤波方法

本发明属于状态估计技术领域，特别涉及一种鲁棒滤波方法。

背景技术：

kalman滤波器广泛应用于线性状态估计领域。但常规的kalman滤波器仅仅在状态空间模型的过程噪声及观测噪声为高斯分布且噪声统计参数完全已知时才能保证其最优性。而在实际工程应用中，由于使用环境的恶劣且多变，状态空间模型的过程噪声及观测噪声可能在某些环境中出现野值，而在另外一些环境中满足理想的高斯分布，且野值在不同的环境中出现的概率也会发生变化。也就是说模型的过程噪声及观测噪声为非静态厚尾分布。单独的采用固定参数的高斯分布或者厚尾分布来对噪声进行建模都会导致一定的随机建模误差，进而恶化kalman滤波器的性能，甚至可能会导致滤波发散。

高斯-student’st混合(gaussian–student’stmixture，gstm)分布可以较为准确的建模非静态厚尾分布，因而具有较强的工程实用价值。基于gstm分布对状态空间模型噪声进行建模进而设计的鲁棒滤波器都依赖于变分贝叶斯(variationalbayes，vb)近似技术。但目前基于vb近似的状态估计方法通常需要进行联合概率密度的自由因子化近似，得到的只是目标概率密度分布的近似解，其准确性通常难以保证。马尔科夫链蒙特卡洛(markovchainmontecarlo，mcmc)方法是一种随机近似方法，被广泛应用于贝叶斯推论当中。mcmc可以通过产生一条以目标分布为准静态分布的马尔科夫链来实现对复杂分布的采样。当采样数量足够大时，mcmc方法可以保证产生足够精确的估计。将mcmc方法用于鲁棒滤波领域有潜力得到比基于vb近似的方法更好的估计精度，而目前还未出现相关的研究。由于其简单易用，gibbs采样器是mcmc框架中应用最为广泛的方法。

技术实现要素：

本发明的目的是：针对于实际线性状态估计应用中经常出现的线性模型过程噪声及观测噪声为非静态厚尾分布的问题，将模型过程噪声及观测噪声建模为gstm分布，进而提出一种基于gibbs采样器的鲁棒滤波方法，可以在模型噪声为非静态厚尾分布时仍然获得较好的状态估计结果。

本发明的技术方案是：一种基于gibbs采样器的鲁棒滤波方法，首先将线性状态空间模型中的过程噪声及观测噪声建模为gstm分布，且混合分布的混合权重及尺度矩阵未知，将其看作随机变量，先验分布分别建模为beta分布及逆wishart分布；进而，引入满足bernoulli分布的指示变量与满足gamma分布的辅助参数，将gstm分布分解为分层高斯的形式；进而将一步预测概率密度与似然概率密度写成分层高斯形式；在gibbs采样器的框架下，将未知的参数与系统状态同时进行迭代采样。每一个迭代周期的流程为：1)以上一个迭代周期采样出来的过程噪声及观测噪声的尺度矩阵，指示变量以及辅助参数、上一个时间历元的系统状态、以及当前时间历元观测为条件，计算当前时间历元的系统状态的后验分布，并从后验分布中采样当前迭代周期中当前时间历元的系统状态；2)以上一个迭代周期采样出来的过程噪声的尺度矩阵，指示变量以及辅助参数，以及当前迭代周期采样出来的当前时间历元系统状态为条件，计算上一个时间历元的系统状态的后验分布，并从后验分布中采样当前迭代周期中上一个时间历元的系统状态；3)以当前迭代周期采样出来的系统状态、上一个迭代周期采样出来的过程噪声及观测噪声尺度矩阵、过程噪声及观测噪声辅助参数、过程噪声及观测噪声指示变量统计参数、观测变量为条件，计算过程噪声及观测噪声指示变量的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的指示变量；4)以当前迭代周期采样出来的过程噪声及观测噪声的指示变量为条件，计算指示变量统计参数的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的指示变量统计参数；5)以当前迭代周期采样出来的系统状态及指示变量、上一个迭代周期采样出来的过程噪声及观测噪声尺度矩阵、观测变量为条件，计算过程噪声及观测噪声辅助参数的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的辅助参数；6)以当前迭代周期采样出来的系统状态、指示变量、辅助参数、观测变量为条件，计算过程噪声及观测噪声尺度矩阵的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的尺度矩阵。在进行多次迭代之后，选取达到稳态之后的迭代周期状态采样的平均值作为当前时间历元最终的状态估计值。选取达到稳态之后的迭代周期状态采样的方差作为最终的状态方差估计值。选取达到稳态之后的迭代周期尺度矩阵采样的平均值作为当前时间历元最终的尺度矩阵估计值。

本发明包括以下步骤：

a.将线性状态空间模型中的过程噪声及观测噪声建模为高斯-student’st混合(gaussian–student’stmixture，gstm)分布，且混合分布的混合权重及尺度矩阵未知，将线性状态空间模型的过程噪声及观测噪声看作随机变量，先验分布分别建模为beta分布及逆wishart分布。

b.引入满足bernoulli分布的指示变量与满足gamma分布的辅助参数，将gstm分布分解为分层高斯的形式；进而将一步预测概率密度与似然概率密度写成分层高斯形式。

c.对于每一个时间历元，执行如下步骤：

c1.在gibbs采样器框架下，对线性状态空间模型的过程噪声及观测噪声尺度矩阵、上一个时间历元的系统状态、指示变量的统计参数以及指示变量进行初始采样，并设置辅助参数的初始采样值。

c2.设置总迭代周期数n以及平稳周期数nb，要求满足nb＜n；进行n次迭代，每一次迭代周期的流程为：

c2-1.以上一个迭代周期采样出来的过程噪声及观测噪声的尺度矩阵、指示变量以及辅助参数、上一个时间历元的系统状态、以及当前时间历元观测为条件，计算当前时间历元的系统状态的后验分布，并从后验分布中采样当前迭代周期中当前时间历元的系统状态。

c2-2.以上一个迭代周期采样出来的过程噪声的尺度矩阵，指示变量以及辅助参数，以及当前迭代周期采样出来的当前时间历元系统状态为条件，计算上一个时间历元的系统状态的后验分布，并从后验分布中采样当前迭代周期中上一个时间历元的系统状态。

c2-3.以当前迭代周期采样出来的系统状态、上一个迭代周期采样出来的过程噪声及观测噪声尺度矩阵、过程噪声及观测噪声辅助参数、过程噪声及观测噪声指示变量统计参数、以及观测变量为条件，计算过程噪声及观测噪声指示变量的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的指示变量。

c2-4.以当前迭代周期采样出来的过程噪声及观测噪声的指示变量为条件，计算指示变量统计参数的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的指示变量统计参数。

c2-5.以当前迭代周期采样出来的系统状态及指示变量、上一个迭代周期采样出来的过程噪声及观测噪声尺度矩阵、观测变量为条件，计算过程噪声及观测噪声辅助参数的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的辅助参数。

c2-6.以当前迭代周期采样出来的系统状态、指示变量、辅助参数、观测变量为条件，计算过程噪声及观测噪声尺度矩阵的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的尺度矩阵。

c3.在进行n次迭代之后，选取达到稳态之后的n-nb次迭代周期状态采样的平均值作为当前时间历元最终的状态估计值；选取达到稳态之后的n-nb次迭代周期状态采样的方差作为最终的状态方差估计值；选取达到稳态之后的n-nb次迭代周期尺度矩阵采样的平均值作为当前时间历元最终的尺度矩阵估计值。

具体的，所述步骤a采用以下方法：

对于线性状态空间模型：

xk＝fkxk-1+wk

zk＝hkxk+vk

其中：下标k表示第k个时间历元的变量；为系统状态向量，为系统观测向量，为状态转移矩阵，为系统观测矩阵，为过程噪声，为观测噪声；

将系统过程噪声及观测噪声建模为gstm分布，即：

p(wk|τk)＝τkn(wk；0,qk)+(1-τk)st(wk；0,qk,ωk)

其中：n(x；a,a)表示随机向量x满足以a为均值，以a为方差的高斯分布；st(x；μ,σ,ω)表示以μ为均值向量、σ为尺度矩阵、ω为自由度的满足student’t分布的随机向量x；τk与分别表示第k个时间历元的混合权重，为未知随机变量，其先验分布满足beta分布：

p(τk)＝beta(τk；e,1-e)

其中：beta(x；a,b)表示以a,b为参数的满足beta分布的随机变量x；e和f为调制参数；将过程噪声与观测噪声gstm分布的参数qk与rk记作尺度矩阵，将其均看作未知的，其先验分布建模为共轭先验：逆wishart分布，如下：

其中：iw(x；a,b)表示随机矩阵x满足以a为自由度，以b为逆尺度矩阵的逆wishart分布；下标i|j代表以第j个时间历元之前及第j个时间历元的系统观测为条件的第i个时间历元的变量估计值；设定每一个时间历元过程噪声及观测噪声尺度矩阵名义值分别为与对于每一个时间历元，令：

得出：

其中：ρu及ρt为调制参数；考虑到实际应用中线性空间模型的过程噪声及观测噪声统计参数通常为慢时变的，定义当前时间历元过程噪声尺度矩阵及观测噪声尺度矩阵的名义值为上一个时间历元的对应后验估计值，即：

具体的，所述步骤b采用以下方法：

引入满足bernoulli分布的随机变量sk与tk，对应的统计参数分别为混合权重τk与即：

其中：bern(x；p)表示以p为参数，满足bernoulli分布的离散随机变量x，且满足x∈{0,1}；过程噪声与观测噪声概率密度函数写作：

进一步，基于student’t分布的分层高斯性质，引入满足gamma分布的辅助参数将其分解为高斯分布与gamma分布的分层组合，即：

其中：gamma(x；a,b)表示以a为形状参数，b为速率参数的满足gamma分布的随机变量x；进而过程噪声与观测噪声概率密度函数写作：

其中：

p(γk)＝gamma(γk；ωk/2,ωk/2)

p(λk)＝gamma(λk；νk/2,νk/2)

根据过程噪声及观测噪声的概率密度函数以及线性空间模型，得到系统一步预测概率密度及似然概率密度分别为：

写成分层概率模型形式为：

具体的，所述步骤c1采用以下方法：

根据所述步骤a得出的过程噪声及观测噪声尺度矩阵先验分布参数直接从中采样得出过程噪声尺度矩阵及观测噪声尺度矩阵初始采样，即：

其中：a^(j)表示随机变量a在第j次迭代当中的采样值；根据上一时间历元所输出的后验状态估计及后验方差矩阵，直接从中采样得出上一个时间历元的系统状态初始采样，即：

根据调制参数e和f，直接从beta(τk；e,1-e)和中采样指示变量的统计参数，即：

进而，根据所采样出来的与从及中采样得出指示变量的初始采样，即：

设置辅助参数的初始采样值为：

具体的，所述步骤c2-1采用以下方法：

在前一迭代周期过程噪声尺度矩阵采样值观测噪声尺度矩阵采样值过程噪声指示变量采样值观测噪声指示变量采样值过程噪声辅助参数采样值观测噪声辅助参数采样值以及上一个时间历元系统状态采样值已知的情况下，采用如下方法计算当前时间历元状态后验分布参数：

1)计算等效过程噪声方差矩阵与观测噪声方差矩阵：

2)计算更新增益：

3)计算当前时间历元后验状态估计值：

4)计算后验方差：

按照：

对当前时间历元的系统状态进行采样；

所述步骤c2-2采用以下方法：

在前一迭代周期过程噪声尺度矩阵采样值指示变量采样值辅助参数采样值以及当前迭代周期采样出来的当前时间历元系统状态采样值已知的情况下，采用如下方法计算上一时间历元状态后验分布参数：

1)计算等效过程噪声方差矩阵：

2)计算更新增益：

3)计算上一时间历元后验状态估计值：

4)计算后验方差：

按照：

对上一时间历元的系统状态进行采样；

所述步骤c2-3采用以下方法：

在当前迭代周期系统状态采样值及上一个迭代周期过程噪声尺度矩阵采样值观测噪声尺度矩阵采样值过程噪声辅助参数采样值观测噪声辅助参数采样值过程噪声指示变量统计参数采样值以及观测噪声指示变量统计参数采样值已知的情况下，采用如下方法计算过程噪声及观测噪声指示变量的后验分布参数：

根据贝叶斯规则：

其对数形式为：

其中：为bernoulli分布，其对数形式为：

根据运动学模型及gstm分布的分层高斯形式将似然概率密度对数形式写作：

其中：

进而得到：

即，sk的后验分布满足bernoulli分布且：

其中：

类似的，对于观测噪声指示变量，其后验分布对数形式满足：

其中:为bernoulli分布，其对数形式为：

根据观测模型及gstm分布的分层高斯形式将似然概率密度对数形式写作：

其中：

进而得到：

即,tk的后验分布满足bernoulli分布且：

其中：

从：

对当前迭代周期过程噪声及观测噪声指示变量进行采样；

所述步骤c2-4采用以下方法：

在当前迭代周期过程噪声指示变量采样值以及观测噪声的指示变量采样值已知的情况下，采用如下方法计算指示变量统计参数的后验分布参数：

根据贝叶斯规则：

其对数形式为：

其中：p(τk)为beta分布，其对数形式为：

logp(τk)＝(e-1)logτk-elog(1-τk)+logc1

为bernoulli分布，其对数形式为：

则：

即，过程噪声指示变量统计参数后验分布满足beta分布：

同样，对于观测噪声指示变量统计参数，有：

且：

得到：

即，过程噪声指示变量统计参数后验分布满足beta分布：

从

对当前迭代周期过程噪声及观测噪声指示变量统计参数进行采样；

所述步骤c2-5采用以下方法：

在当前迭代周期系统状态采样值及指示变量采样值及上一个迭代周期过程噪声尺度矩阵采样值以及观测噪声尺度矩阵采样值已知的情况下，采用如下方法计算过程噪声及观测噪声辅助参数的后验分布参数：

根据贝叶斯规则：

其对数形式为：

其中：p(γk)＝gamma(γk；ωk/2,ωk/2)，其对数形式为：

根据运动学模型及gstm分布的分层高斯形式得到的对数形式为：

其中：

则：

其中：

即，过程噪声辅助参数后验分布满足gamma分布：

同样，对于观测噪声辅助参数，有：

且：

其中：

得到：

其中：

即，观测噪声辅助参数后验分布满足gamma分布：

从

对过程噪声及观测噪声辅助参数进行采样；

所述步骤c2-6采用以下方法：

在当前迭代周期系统状态采样值及指示变量采样值及辅助参数采样值及已知的情况下，采用如下方法计算过程噪声及观测噪声尺度矩阵的后验分布参数：

根据贝叶斯规则：

其对数形式为：

过程噪声尺度矩阵建模为逆wishart分布，其对数形式为：

根据运动学模型及gstm分布的分层高斯形式得到的对数形式为：

计算得到：

其中：

也就是说，过程噪声尺度矩阵的后验概率密度为以为自由度，以为逆尺度矩阵的逆wishart分布；

类似的，对于观测噪声尺度矩阵，有：

其中：

计算得到：

其中：

也就是说，观测噪声方差矩阵的后验概率密度为以为自由度，以为逆尺度矩阵的逆wishart分布；按照：

采样与r⁽ⁱ⁺¹⁾，得到当前迭代周期的过程噪声尺度矩阵及观测噪声尺度矩阵采样值。

具体的，所述步骤c3采用以下方法：

经历n次迭代采样后，得到当前时间历元状态变量采样集合为选取达到稳态之后的n-nb次迭代周期采样的平均值作为最终的状态估计值，即：

状态方差估计值的计算方法为：

选取当前时间历元过程噪声尺度矩阵及观测噪声尺度矩阵估计值为：

有益效果：本发明通过将线性空间模型当中的过程噪声及观测噪声建模为gstm分布来处理实际状态估计问题中经常出现的非静态厚尾噪声。通过引入满足bernoulli分布的指示变量及满足gamma分布的辅助参数将gstm分布分解为分层gauss形式。采用gibbs采样器来对线性状态空间模型当中未知的系统状态和参数进行随机采样，得到基于gibbs采样器的鲁棒滤波方法。所提出的方法可以在模型噪声为非静态厚尾分布时仍然取得较好的状态估计结果，对传感器的观测野值及噪声特性的变化具有一定的鲁棒性。

附图说明

图1为本发明的流程图；

图2为本发明即gib-gstm与kf、gib-kf、vb-gstm、vb-rst五种方法的位置估计均方根误差比较；

图3为本发明即gib-gstm与kf、gib-kf、vb-gstm、vb-rst五种方法的速度估计均方根误差比较。

具体实施方式

实施例1，参见附图1，一种基于gibbs采样器的鲁棒滤波方法，包括以下步骤：

a.将线性状态空间模型中的过程噪声及观测噪声建模为gstm分布，且混合分布的混合权重及尺度矩阵未知，将线性状态空间模型的过程噪声及观测噪声看作随机变量，先验分布分别建模为beta分布及逆wishart分布。

对于线性状态空间模型：

xk＝fkxk-1+wk

zk＝hkxk+vk

其中：下标k表示第k个时间历元的变量；为系统状态向量，为系统观测向量，为状态转移矩阵，为系统观测矩阵，为过程噪声，为观测噪声；

将系统过程噪声及观测噪声建模为gstm分布，即：

p(wk|τk)＝τkn(wk；0,qk)+(1-τk)st(wk；0,qk,ωk)

其中：n(x；a,a)表示随机向量x满足以a为均值，以a为方差的高斯分布；st(x；μ,σ,ω)表示以μ为均值向量、σ为尺度矩阵、ω为自由度的满足student’t分布的随机向量x；τk与分别表示第k个时间历元的混合权重，为未知随机变量，其先验分布满足beta分布：

p(τk)＝beta(τk；e,1-e)

其中：beta(x；a,b)表示以a,b为参数的满足beta分布的随机变量x；e和f为调制参数；将过程噪声与观测噪声gstm分布的参数qk与rk记作尺度矩阵，将其均看作未知的，其先验分布建模为共轭先验：逆wishart分布，如下：

其中：iw(x；a,b)表示随机矩阵x满足以a为自由度，以b为逆尺度矩阵的逆wishart分布；下标i|j代表以第j个时间历元之前及第j个时间历元的系统观测为条件的第i个时间历元的变量估计值；设定每一个时间历元过程噪声及观测噪声尺度矩阵名义值分别为与对于每一个时间历元，令：

得出：

其中：ρu及ρt为调制参数；考虑到实际应用中线性空间模型的过程噪声及观测噪声统计参数通常为慢时变的，定义当前时间历元过程噪声尺度矩阵及观测噪声尺度矩阵的名义值为上一个时间历元的对应后验估计值，即：

b.引入满足bernoulli分布的指示变量与满足gamma分布的辅助参数，将gstm分布分解为分层高斯的形式；进而将一步预测概率密度与似然概率密度写成分层高斯形式。

引入满足bernoulli分布的随机变量sk与tk，对应的统计参数分别为混合权重τk与即：

其中：bern(x；p)表示以p为参数，满足bernoulli分布的离散随机变量x，且满足x∈{0,1}；过程噪声与观测噪声概率密度函数写作：

进一步，基于student’t分布的分层高斯性质，引入满足gamma分布的辅助参数将其分解为高斯分布与gamma分布的分层组合，即：

其中：gamma(x；a,b)表示以a为形状参数，b为速率参数的满足gamma分布的随机变量x；进而过程噪声与观测噪声概率密度函数写作：

其中：

p(γk)＝gamma(γk；ωk/2,ωk/2)

p(λk)＝gamma(λk；νk/2,νk/2)

根据过程噪声及观测噪声的概率密度函数以及线性空间模型，得到系统一步预测概率密度及似然概率密度分别为：

写成分层概率模型形式为：

c.对于每一个时间历元，执行如下步骤：

c1.在gibbs采样器框架下，对线性状态空间模型的过程噪声及观测噪声尺度矩阵、上一个时间历元的系统状态、指示变量的统计参数以及指示变量进行初始采样，并设置辅助参数的初始采样值。

根据所述步骤a得出的过程噪声及观测噪声尺度矩阵先验分布参数直接从中采样得出过程噪声尺度矩阵及观测噪声尺度矩阵初始采样，即：

其中：a^(j)表示随机变量a在第j次迭代当中的采样值；根据上一时间历元所输出的后验状态估计及后验方差矩阵，直接从中采样得出上一个时间历元的系统状态初始采样，即：

根据调制参数e和f，直接从beta(τk；e,1-e)和中采样指示变量的统计参数，即：

进而，根据所采样出来的与从及中采样得出指示变量的初始采样，即：

设置辅助参数的初始采样值为：

c2.设置总迭代周期数n以及平稳周期数nb，要求满足nb＜n；进行n次迭代，每一次迭代周期的流程为：

c2-1.以上一个迭代周期采样出来的过程噪声及观测噪声的尺度矩阵、指示变量以及辅助参数、上一个时间历元的系统状态、以及当前时间历元观测为条件，计算当前时间历元的系统状态的后验分布，并从后验分布中采样当前迭代周期中当前时间历元的系统状态。

在前一迭代周期过程噪声尺度矩阵采样值观测噪声尺度矩阵采样值过程噪声指示变量采样值观测噪声指示变量采样值过程噪声辅助参数采样值观测噪声辅助参数采样值以及上一个时间历元系统状态采样值已知的情况下，采用如下方法计算当前时间历元状态后验分布参数：

1)计算等效过程噪声方差矩阵与观测噪声方差矩阵：

2)计算更新增益：

3)计算当前时间历元后验状态估计值：

4)计算后验方差：

按照：

对当前时间历元的系统状态进行采样。

c2-2.以上一个迭代周期采样出来的过程噪声的尺度矩阵，指示变量以及辅助参数，以及当前迭代周期采样出来的当前时间历元系统状态为条件，计算上一个时间历元的系统状态的后验分布，并从后验分布中采样当前迭代周期中上一个时间历元的系统状态。

在前一迭代周期过程噪声尺度矩阵采样值指示变量采样值辅助参数采样值以及当前迭代周期采样出来的当前时间历元系统状态采样值已知的情况下，采用如下方法计算上一时间历元状态后验分布参数：

1)计算等效过程噪声方差矩阵：

2)计算更新增益：

3)计算上一时间历元后验状态估计值：

4)计算后验方差：

按照：

对上一时间历元的系统状态进行采样。

c2-3.以当前迭代周期采样出来的系统状态、上一个迭代周期采样出来的过程噪声及观测噪声尺度矩阵、过程噪声及观测噪声辅助参数、过程噪声及观测噪声指示变量统计参数、以及观测变量为条件，计算过程噪声及观测噪声指示变量的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的指示变量。

在当前迭代周期系统状态采样值及上一个迭代周期过程噪声尺度矩阵采样值观测噪声尺度矩阵采样值过程噪声辅助参数采样值观测噪声辅助参数采样值过程噪声指示变量统计参数采样值以及观测噪声指示变量统计参数采样值已知的情况下，采用如下方法计算过程噪声及观测噪声指示变量的后验分布参数：

根据贝叶斯规则：

其对数形式为：

其中：为bernoulli分布，其对数形式为：

根据运动学模型及gstm分布的分层高斯形式将似然概率密度对数形式写作：

其中：

进而得到：

即，sk的后验分布满足bernoulli分布且：

其中：

类似的，对于观测噪声指示变量，其后验分布对数形式满足：

其中:为bernoulli分布，其对数形式为：

根据观测模型及gstm分布的分层高斯形式将似然概率密度对数形式写作：

其中：

进而得到：

即,tk的后验分布满足bernoulli分布且：

其中：

从：

对当前迭代周期过程噪声及观测噪声指示变量进行采样。

c2-4.以当前迭代周期采样出来的过程噪声及观测噪声的指示变量为条件，计算指示变量统计参数的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的指示变量统计参数。

在当前迭代周期过程噪声指示变量采样值以及观测噪声的指示变量采样值已知的情况下，采用如下方法计算指示变量统计参数的后验分布参数：

根据贝叶斯规则：

其对数形式为：

其中：p(τk)为beta分布，其对数形式为：

logp(τk)＝(e-1)logτk-elog(1-τk)+logc1

为bernoulli分布，其对数形式为：

则：

即，过程噪声指示变量统计参数后验分布满足beta分布：

同样，对于观测噪声指示变量统计参数，有：

且：

得到：

即，过程噪声指示变量统计参数后验分布满足beta分布：

从

对当前迭代周期过程噪声及观测噪声指示变量统计参数进行采样。

c2-5.以当前迭代周期采样出来的系统状态及指示变量、上一个迭代周期采样出来的过程噪声及观测噪声尺度矩阵、观测变量为条件，计算过程噪声及观测噪声辅助参数的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的辅助参数。

在当前迭代周期系统状态采样值及指示变量采样值及上一个迭代周期过程噪声尺度矩阵采样值以及观测噪声尺度矩阵采样值已知的情况下，采用如下方法计算过程噪声及观测噪声辅助参数的后验分布参数：

根据贝叶斯规则：

其对数形式为：

其中：p(γk)＝gamma(γk；ωk/2,ωk/2)，其对数形式为：

根据运动学模型及gstm分布的分层高斯形式得到的对数形式为：

其中：

则：

其中：

即，过程噪声辅助参数后验分布满足gamma分布：

同样，对于观测噪声辅助参数，有：

且：

其中：

得到：

其中：

即，观测噪声辅助参数后验分布满足gamma分布：

从

对过程噪声及观测噪声辅助参数进行采样。

c2-6.以当前迭代周期采样出来的系统状态、指示变量、辅助参数、观测变量为条件，计算过程噪声及观测噪声尺度矩阵的后验分布，并从后验分布当中采样当前迭代周期过程噪声及观测噪声的尺度矩阵。

在当前迭代周期系统状态采样值及指示变量采样值及辅助参数采样值及已知的情况下，采用如下方法计算过程噪声及观测噪声尺度矩阵的后验分布参数：

根据贝叶斯规则：

其对数形式为：

过程噪声尺度矩阵建模为逆wishart分布，其对数形式为：

根据运动学模型及gstm分布的分层高斯形式得到的对数形式为：

计算得到：

其中：

也就是说，过程噪声尺度矩阵的后验概率密度为以为自由度，以为逆尺度矩阵的逆wishart分布；

类似的，对于观测噪声尺度矩阵，有：

其中：

计算得到：

其中：

也就是说，观测噪声方差矩阵的后验概率密度为以为自由度，以为逆尺度矩阵的逆wishart分布；按照：

采样与r⁽ⁱ⁺¹⁾，得到当前迭代周期的过程噪声尺度矩阵及观测噪声尺度矩阵采样值。

c3.在进行n次迭代之后，选取达到稳态之后的n-nb次迭代周期状态采样的平均值作为当前时间历元最终的状态估计值；选取达到稳态之后的n-nb次迭代周期状态采样的方差作为最终的状态方差估计值；选取达到稳态之后的n-nb次迭代周期尺度矩阵采样的平均值作为当前时间历元最终的尺度矩阵估计值。

经历n次迭代采样后，得到当前时间历元状态变量采样集合为选取达到稳态之后的n-nb次迭代周期采样的平均值作为最终的状态估计值，即：

状态方差估计值的计算方法为：

选取当前时间历元过程噪声尺度矩阵及观测噪声尺度矩阵估计值为：

实施例2，实施例1所实现的伪代码为：

实施例3，利用实施例1所述的鲁棒滤波方法通过仿真数据进行验证。

仿真案例为二维目标跟踪，定义状态向量xk＝[xkykvx,kvy,k]^t，其中，xk，yk，vx,k以及vy,k分别表示x与y方向的位置以及x与y方向的速度。观测变量为x与y方向的位置。系统运动学模型及观测模型分别为：

其中：δt为离散时间间隔，在本仿真中选择为1秒。总仿真时间t选择为200秒。设置真实的过程噪声以及观测噪声统计特性为非静态厚尾的：

其中：

名义过程噪声方差矩阵及观测噪声方差矩阵分别设置为其对应的真实值，即仿真调制参数设置为ρt＝5,ρu＝5,e＝0.85,f＝0.85,ω0＝0.5,ν0＝0.5，总迭代周期n＝5000，平稳周期为nb＝1000。

作为比较，本实施例同时展示了以名义噪声参数解算的kalman滤波器(kf)、基于噪声参数未知的gauss分布，采用gibbs采样器进行求解的自适应kalman滤波器(gib-kf)、基于噪声gstm分布，采用vb近似求解的鲁棒滤波器(vb-gstm)以及基于噪声student’st分布，采用vb近似求解的鲁棒滤波器(vb-rst)的估计结果。本发明的方法简记作“gib-gstm”。

500独立的次montecarlo仿真的结果被用来验证本发明所提出的方法。位置和速度的均方根误差(rmse)及平均均方根误差(armse)被用来评价不同方法的性能。几种评价指标的计算方式如下：

其中：及为第i次montecarlo仿真中航行器在k时刻真实的坐标；及为第i次montecarlo仿真中航行器在k时刻估计出的位置坐标；及为第i次montecarlo仿真中航行器在k时刻真实的速度；及为第i次montecarlo仿真中航行器在k时刻估计出的速度；m＝500表示总的仿真次数。

附图2为五种方法的位置估计均方误差比较，附图3为五种方法的速度均方误差比较。五种方法的平均均方根误差如表1所示。

表1

根据图2、图3以及表1，可以看出相比于kf、gib-kf、vb-gstm与vb-rst，实施例1所提出的gib-gstm可以在噪声特性为非静态厚尾的情况下获得更好的状态估计精度，包括位置估计以及速度估计。在系统噪声满足gauss分布时，实施例1所提出的方法估计性能略逊于kf、gib-kf、vb-gstm与vb-rst，而在系统噪声出现野值时，实施例1所提出的方法估计性能明显优于kf、gib-kf、vb-gstm与vb-rst。总体来说，实施例1所提出的gib-gstm可以适应非静态噪声特性环境，且可以在每一种噪声特性下都取得较为理想的估计性能，有较好的实际应用潜力。

虽然，上文中已经用一般性说明及具体实施例对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。