贝叶斯序贯重要性积分滤波器

张宏伟，张小虎，曹 勇

（1. 中山大学航空航天学院，广东广州 510725；2. 中国科学院空间精密测量技术重点实验室，陕西西安 710119；3. 北京东方计量测试研究所，北京 100029）

1 引言

滤波是指从含有噪声的传感测量中估计目标的未知状态或参数，非线性滤波是解决该类不确定信号处理问题的核心技术，已被广泛应用于工程、统计和数学领域［1～8］. 贝叶斯推理为滤波算法提供了一种数学机制，然而在实际应用中，由于非线性动态系统函数的复杂性，难以对方程的积分运算进行闭式求解，通常需利用数值方法近似解算，主要包括数值积分和蒙特卡罗两大类.第一类数值积分方法基于确定性法则. 高斯匹配法［1］近似目标分布简单直接，然而数值精度较为一般. 扩展卡尔曼滤波（Extended Kalman Filter，EKF）算法［2］基于泰勒级数展开，在工程上最易实现且应用最为广泛，然而在复杂传感测量环境中会引入截断误差和基点误差. 无迹卡尔曼滤波（Unscented Kalman Filter，UKF）算法［3］选择固定数量的sigma 点进行无迹变换以近似目标原始分布，在满足协方差矩阵为正定的前提条件下，可达到3阶数值精度. 当状态空间的积分维度适中且满足假定概率模型时，Gauss-Hermite积分算法能够最大程度地减小近似误差［4］. 通常，该类改进卡尔曼滤波（Kalman Filter，KF）算法的应用局限于高斯噪声扰动，因此制约了确定性数值积分方法在强非线性场景中的应用［5］. 第二类蒙特卡罗方法基于随机采样的思想. 该方法能够有效处理非线性非高斯滤波问题，无需对非线性函数进行线性化处理. 马尔可夫链蒙特卡罗（Monte Carlo Markov Chain，MCMC）算法［6］通过模拟与目标分布一致的马尔可夫链生成有效的蒙特卡罗样本，然而在样本较少的情况下难以维持状态为稳态分布，增大数据采集量将导致计算复杂度和时间成本急剧增大. 与此不同，重要性采样（Importance Sampling，IS）算法的优势在于可通过调整样本权重以校正目标重要性函数和目标真实分布之间的偏差，因此，构建能够覆盖有效似然区域的重要性函数至为关键［7］.Doucet等人［8］提出序贯蒙特卡罗采样也称为粒子滤波（Particle Filter，PF）方法，在贝叶斯滤波框架下推导了离散时间动态模型的重要性采样法则. 理论上，该类蒙特卡罗方法的采样误差具有维数不变性，数学上的大数中心极限定理保证其计算收敛性. 然而由维数灾难造成的样本急剧衰减和权值退化仍是亟待解决的技术难题.

从贝叶斯概率论的数学模型分析，复杂非线性高维状态空间中的参数估计属于模糊测度，难以精确解析［9，10］. 其根本原因不在于外在的统计变化量，而是源于复杂非线性动态系统的内在模糊性，主要包括目标模型属性歧义和非线性不稳定传输［11］. 该模糊性制约了许多单个非线性滤波算法的应用，许多研究者提出了结合上述两类数值方法的改进算法. 如通过EKF，UKF 及Gauss-Hermite积分等非线性滤波算法生成重要性函数，即扩展卡尔曼粒子滤波（Extended Kalman Particle Filter，EPF）［12］，无迹卡尔曼粒子滤波（Unscented Kalman Particle Filter，UPF）［13］，高斯-施密特粒子滤波（Gauss-Hermite Particle Filter，GHPF）［14］.Garcia-Fernandez等人［15］证明，当前时刻观测信息较为准确时，改进第一、二阶矩无法直接调制重要性函数和目标分布之间的不匹配问题.

另外，多模型估计器广泛应用于非线性动态系统的状态空间的细化［11］. 其中，交互多模型（Interacting Multiple Model，IMM）算法被成功应用于目标机动跟踪，该算法通过融合有限个单滤波器的滤波输出估计机动目标状态. 事实上，无法通过有限个目标状态模型表征全域的目标机动行为，同时不必要的模型竞争会加大计算复杂度［16］. Hostettler 等人［17］根据Rao-Blackwell理论对状态空间降维，分别采用线性高斯逼近和粒子滤波方法估计系统的状态模型参数和模型状态演化. 王琪等人［18］采用Rao-Blackwell 粒子滤波方法，有效缓解了相控阵导引头指向误差斜率给导弹制导系统带来的寄生回路振荡问题. 该类模型驱动的滤波方法依赖具体问题，对于复杂非线性动态系统的状态估计，难以平衡模型切换的准确性和计算复杂度.

从香农信息论的角度分析，充分利用系统的约束信息有助于实现状态空间降维，提高模型预测的可靠性和状态滤波的准确性［10，19，20］. 近年来，约束非线性滤波算法发展迅速，其优势在于可根据应用场景选取辅助条件. 现有公开文献已提出截断观测驱动的约束UKF［19］、融入状态约束的PF［20］、加载深度空间语义特征的PF［21］及基于注意力学习的正则化相关滤波［22］等非线性动态状态估计和目标跟踪算法. 该类算法依赖隐马尔可夫链的动力学模型［23］，在强非线性环境中统计迭代易陷入局部最优，缺乏较强的泛化能力.

本质上，观测量和状态变量之间存在因果映射［11，24］. 然而实际上，由于过程噪声、观测噪声和目标机动运动等不确定因素，在非线性滤波中存在观测似然、重要性函数和真实目标后验分布失配的现象，本文作者在专利中分析了产生该现象的本质原因［25］. 通常，确定的等式或不等式等硬约束条件会破坏常规的高斯分布假设，因此难以对严格的数值边界条件进行闭式解析. 软约束技术能够有效表征非线性动态系统中不确定变量之间的模糊性［10，19，20］. 文献［26］将统计熵的概念推广到不确定信息处理的度量，理论上可通过相关信息熵测度多模态似然对状态估计的影响. 基于此，本文推导了贝叶斯序贯重要性积分滤波（Bayesian Sequential Importance Quadrature Filter，SIQF）算法. 构建系统软空时约束下的有界似然，以消减由非线性动态系统的不可预测性造成的似然分布和目标分布之间的偏差. 通过并行Gauss-Hermite 积分将修正先验和反馈补偿融入状态演化过程，引入最大相关熵调制重要性函数和真实目标分布的一致性，从而增强重要性采样的多样性和预测更新的容错性.

2 问题描述

2.1 系统模型

本文利用空时约束信息构建多模态似然以优化系统状态估计，记为模型集M={M1，…，MK} ，K≥2. 该集合包含K个相互独立的截断先验和后验反馈等元素，将其融入贝叶斯递推过程. 由此，离散非线性动态系统模型可表示为

其中，k∈N 是离散时刻；Xk∈RnX和Zk∈RnZ分别是状态和观测矢量，RnX和RnZ分别是nX维和nZ维空间模型；f和h分别是非线性状态方程和观测方程；υk和ek分别是随机过程噪声和观测噪声，其概率密度函数分别记为pυ(υk)和pe(ek)，假设为相互独立的0 均值白噪声，协方差分别为Συ和Σe；M i k表示第i个模型.

状态变量的初始先验记为p0(Xk)，转移概率为隐马尔可夫序列，似然函数关联系统状态Xk和观测Zk之间的所有信息. 根据贝叶斯全概率公式［1，8］，由观测量推断非线性动态系统中未知状态的后验概率密度为

其中，分母中的归一化常数p(Z1：k)的积分运算难以解算，因此多区域似然的后验分布近似为多高斯分布；和分别是第i个模型下状态分布的均值和方差，分别定义为

在实际应用中，目标状态通常在直角坐标系推演，而观测数据在球面坐标系或者像坐标系传输. 由于类KF 算法受限于高斯假设，因此无法通过提高矩的精度提升整体滤波性能［2，5，11］. 由此可知，在非线性较强的滤波环境中，式（4）和式（5）难以解算. 从算法机理角度出发，重要性采样方法能够有效处理非线性非高斯问题. 然而，对于似然函数、重要性函数偏离目标分布的情况，该算法会发生样本急剧衰减和权值退化，造成滤波发散. 针对该问题，本文构建有界似然并推导能够覆盖多似然区域的多模重要性积分算法.

2.2 提出算法框架

非线性滤波近似的主要困难在于更新步骤产生的误差［5，11，16］，较为常见的例子是若当前测量包含更多有效信息，而过多依赖先验近似目标后验分布会造成较大偏差. 在实际的物理动态系统模型中，传感观测系统提供的测量信息是有限的，因此观测似然应为有界分布. 基于此，通过截断观测噪声概率密度pe(·)构建可行域Λ(Zk)，即

其中，h-1( ·)是式（2）中观测方程h的反函数.

对于存在多域似然的滤波问题，目标状态后验应为多模分布. 因此将其近似为加权混合高斯分布，即

其中，和分别是第i个似然高斯分布的权重、均值和协方差.

综上所述，本文构建了系统空时软约束模型，将有界似然和后验反馈等先验信息融入状态演化过程，以提高整体滤波性能. 考虑动态系统的随机干扰和约束函数的模糊性，为了避免陷入局部最优或鞍点值，采用全局牛顿内点法遍历可行域迭代求解帕累托解，以近似当前有效观测驱动的修正先验. 为了处理非线性非高斯估计问题，通过并行Gauss-Hermite 递推状态以生成混合重要性分布，综合改善重要性采样的多样性和预测协方差的容错性. 图1为算法流程.

图1 算法流程框架

3 贝叶斯序贯重要性积分滤波算法的实现

3.1 近似有界似然先验

根据式（6）定义的可行域Λ(Zk)可构建当前观测驱动的有界似然. 同时，在贝叶斯全概率公式中引入辅助变量集a，目标状态的后验概率密度函数修正为

其中，plik(·)和p1(·)分别是似然函数和修正先验的概率密度 函数；ς为 常数；是指示函数，若取为1，否则为0.

在近似修正先验的数值规划过程中，需分解目标状态的标量和矢量，选择标量分量构成凸的目标函数.实际应用中，可通过最小二乘和多视图对极几何约束等测量手段将目标初始位置X0定位在可行域Λ(Zk)内. 从几何上看，可行域的中心处观测误差应为最小.因此，构建目标函数为

在可行域内，通过backtracking 回溯算法求解真实目标状态分布的均值，即

其中，s和dk分别是迭代步长和搜索方向，该数值规划优化方法的详细步骤请参考本文作者之前的工作［19］. 根据该统计中心计算可行域的第二阶矩，即方差.

综上所述，在第k时刻，可行域的中心或修正先验概率密度函数记为其中，下标“1”用以区分原始先验概率. 相应地，原始先验概率分布记为这里请注意，可通过扩维方法获取目标状态矢量集［3］. 为了覆盖多似然区域，通过加权p1( ·)和p0( ·)修正目标状态先验为

其中，α是权重因子，用来衡量不同区域似然先验对目标状态估计的影响. 其求解过程将在下一节详细描述.

3.2 并行Gauss-Hermite 生成重要性分布

根据Gauss-Hermite 积分法则［5，6］，可通过有限个一维高斯分布迭代近似多维高斯分布的积分，即

其中，I表示单位矩阵是积分点及其权值；l是Hermite 多项式的特征根个数；是方差的Cholesky解，即

3.2.1 预测

在第k-1 时刻，从修正先验分布选取2l-1个nX维积分点为

状态预测的均值和协方差为

测量预测均值为

3.2.2 更新

计算测量的预测协方差和测量的互协方差为

计算卡尔曼滤波增益为

滤波输出均值为

滤波输出协方差为

至此，采用最大相关信息熵度量本文2个先验分布之间的相似性. 构建最大似然目标函数，以求解权重因子α，即

其中，gi( ·)是第i个似然分布的核函数，定义为

其中，exp( ·)是以自然常数e为底的指数函数.

3.3 重要性采样更新

对应式（13）的多模先验分布，根据多模重要性分布采样重要性样本，即

计算样本权重为

归一化重要性样本的权值为

根据蒙特卡罗法则，通过狄拉克函数δ( ·)加权近似目标状态的后验概率密度函数为

在已知目标状态后验分布的情况下，加权重要性样本估计真实目标状态分布的均值和方差为

4 仿真实验与分析

本节研究以下2 个例子. 第一个例子是一维单变量非平稳增长模型（Univariate Nonstationary Growth Model，UNGM），将SIQF 算法与文献［13］和文献［19］中的扩展卡尔曼滤波（EKF）、无迹卡尔曼滤波（UKF）、传统粒子滤波（Generical Particle Filter，GPF）、扩展卡尔曼粒子滤波和无迹卡尔曼滤波算法比较；第二个例子采用我国某通用民航公司实测民航航迹作为目标数据，观测传感器为雷达. 由于测量环境干扰等不确定因素，部分观测量之间时间间隔达数十秒，故称之为稀疏观测. 在该稀疏观测场景中，EKF 和UKF 跟踪误差过大甚至发散［19］. 因此，将SIQF 算法与文献［11］、文献［16］和文献［18］中的交互多模型扩展卡尔曼滤波（Interacting Multiple Model Extended Kalman Filter，IMMEKF）、交互多模型无迹卡尔曼滤波（Interacting Multiple Model Unscented Kalman Filter，IMMUKF）以及基于Rao-Blackwell 的多模型粒子滤波（Multiple Model Rao-Blackwell Particle Filter，MMRBPF）算法比较.

4.1 一维单变量非平稳增长模型

一维单变量非平稳增长模型的状态方程和观测方程模型分别为

其中，系数为已知常数，a=0.5，b=2.5，c=8，φ1=0.2，φ2=0.5；状态初值X0在[0，1]范围内平均取值；υk是服从参数为[2，3]的伽马分布的过程噪声；ek是均值为0且方差为0.5的观测噪声.

构建合适的目标重要性函数是重要性采样算法的关键步骤，它与实际目标分布的匹配程度直接影响整体滤波性能和样本数目及计算复杂度.

图2～5 分别是采样10，200，300，500 个样本的结果. 从滤波误差趋势的定性比较分析可知：相比EKF和UKF 算法，GPF 显示出处理非高斯噪声的优势.EPF 和UPF 算法分别由EKF 和UKF 算法的非线性估计生成重要性函数.Gauss-Hermite 方法估计精度高于EKF，且状态方差有较重的尾部，因此SIQF滤波性能优于UPF.

图2 10个样本

图3 200个样本

图4 300个样本

图5是样本增多到500 个，单独比较基于蒙特卡罗采样的数值方法的滤波性能. 随着时间推进，在状态演化过程中，EPF，UPF 和GPF 滤波算法的估计误差在第30 s 后明显增大，出现这个现象的主要原因是重要性函数与实际目标的真实分布失配引起样本多样性衰减和权值退化. 相比UPF 滤波算法，SIQF 算法在无需牺牲计算复杂度的情况下将平均误差减小了63%. 这主要是因为SIQF 融合截断先验和后验反馈构建覆盖多似然的目标重要性函数，从根本上增强了重要性样本的多样性和预测准确性.

图5 500个样本

4.2 稀疏观测机动目标跟踪

图6 是我国某通用民航公司广播式自动相关监视（Automatic Dependent Surveillance Broadcast，ADS-B）系统监测的某民航飞机航迹信息. 其中，黑色曲线是由全球定位系统（Golobal Positional System，GPS）接收的航迹，红色曲线是监测雷达系统的观测数据. 目标状态变量由位置、速度和转弯率组成，表示为Xk=，其中ωk为第k时刻目标的转弯率.目标状态转移矩阵和过程噪声分别为

图6 目标GPS定位和观测航迹

其中，T为雷达的采样时间间隔，取为1 s；过程噪声标准差为σv=0.1 km/s2，σω=0.1 rad/s2. 第k时刻的雷达测距rk和测角θk可由下式计算：

其中，(Xs，Ys)是雷达的位置坐标，观测噪声标准差分别为σr=0.1 km和σθ=3 mrad.

图7、图8 和图9 分别给出了IMMEKF，IMMUKF，MMRBPF 和本文提出的SIQF 滤波算法在二维空间、X轴和Y轴运动方向上跟踪机动目标的位置均方根误差. 从该4 种滤波算法的均方根误差趋势的定性比较可知：相比传统交互多模型滤波算法和MMRBPF 滤波算法，SIQF 算法的跟踪效果在滤波精度方面表现出较大优势. 这主要是因为SIQF 算法采用最大相关信息熵准则衡量截断先验和后验反馈约束信息，在多模型参数下构建的混合高斯目标重要性函数覆盖了多似然信息，有效减小了目标重要性函数和实际目标状态的真实分布之间的偏差，从而提高了系统状态在动态演化过程中对随机噪声的抗干扰能力.

图7 二维空间的均方根误差

图8 X轴方向的均方根误差

图9 Y轴方向的均方根误差

另外，为了分析SIQF算法的计算复杂度. 表1统计了IMMEKF，IMMUKF，MMRBPF 和SIQF 这4 种滤波器执行100 轮蒙特卡罗实验的误差均值和所需的平均执行时间. 从这2 个参数的定量比较可知：相比传统的交互多模型滤波算法，SIQF 滤波算法的蒙特卡罗实验执行时间有所增加，主要原因在于近似有界似然和构建重要性函数过程中的迭代约束优化. 相比MMRBPF 滤波算法，SIQF 滤波算法取得与之相当的滤波精度，而计算时间降低了一个数量级. 这主要是因为引入最大相关信息熵调制多模似然和目标真实分布之间的匹配程度，而不是依靠增多模型个数以表征目标机动行为，从而避免了多模型之间不必要的竞争. 同时，引入特征辅助变量加强了该算法学习有效似然先验知识的自适应性，在理论上实现了状态空间降维.

表1 100轮蒙特卡罗轮实验的滤波误差均值平均执行时间

5 结束语

本文分析了非线性滤波中观测似然、重要性函数与目标真实分布失配的根本原因，为了消减根源于模型歧义和预测有偏的2种模糊性，根据软约束理论修正有界似然，引入了最大相关信息熵准则构建覆盖多似然模态的贝叶斯序贯重要性积分滤波算法. 仿真实验证明，该算法充分融合Gauss-Hermite 的高精度和重要性采样算法的灵活性，提高了贝叶斯滤波器抗击多峰似然和非高斯噪声干扰的能力. 结合本文作者的课题研究，后续工作将基于空间目标星图和光度曲线信息，研究分布式贝叶斯重要性积分算法跟踪和反演空间目标姿态.