马天力， 张 扬， 刘 盼， 高 嵩
(西安工业大学电子信息工程学院，西安，710021)
目标跟踪是在已知受噪声干扰的量测信息条件下，利用相关滤波算法对目标状态(位置，速度等)进行估计的过程，其在弹道导弹防御、空中目标监控和水下目标探测等领域应用广泛[1-2]。然而系统的量测数据与目标动态参数间的关系大多为非线性，一般采用点估计和概率密度估计[3]的方法对此类系统状态进行计算。点估计方法是通过非线性逼近策略输出状态均值和协方差矩阵获得系统状态。根据逼近策略可分为函数逼近、统计量逼近和随机模型逼近；
其代表算法分别为扩展卡尔曼滤波[4]、无迹卡尔曼滤波[5]和容积卡尔曼滤波(cubature kalman filter, CKF)[6]。概率密度估计则是通过构建状态的概率密度表达，计算状态后验概率密度函数，并根据后验信息获得系统状态[7]，该方法需要大量状态对后验概率密度函数进行描述，计算量较大，因此，实际系统大多采用点估计器进行状态估计。
一般来说，跟踪系统中量测噪声的统计特性往往是未知且时变的。针对该类问题，文献[8]提出改进Sage-Husa自适应Kalman滤波算法，运用协方差匹配判据来判断滤波发散趋势，并引入自适应衰减因子修正预测误差协方差，从而调整滤波增益阵，抑制滤波发散。文献[9]提出基于最大后验概率的自适应无迹卡尔曼滤波算法，利用输出的量测信息在线更新噪声均值和协方差。Chang等[10]通过Huber代价函数重新表述量测噪声，无需对非线性方程进行线性化处理，提高了滤波算法的估计精度。彭美康等[11]提出基于鲁棒M估计的自适应容积卡尔曼滤波算法，通过修正观测向量的误差协方差，将损失函数最小化，在观测噪声污染率较高的情况下能够自适应抑制野值对系统的干扰。Sarkka等[12]提出了变分贝叶斯自适应卡尔曼滤波算法，将变分贝叶斯方法与卡尔曼滤波结合，对模型状态和时变量测噪声参数进行联合递归估计，一定程度上降低时变量测噪声方差对系统的影响。虽然上述方法可以有效解决时变量测噪声条件下非线性系统状态估计问题，但均需要满足量测噪声服从高斯分布的假设，具有一定局限性[13]。
解决时变非高斯噪声条件下的非线性滤波问题主要核心是对非高斯噪声进行描述。Fan等[14]运用混合高斯模型对非高斯量测噪声进行表示，利用集成卡尔曼滤波方法估计系统状态；
董祥祥等[15]运用Pearson Type VII分布对重尾噪声进行建模，结合鲁棒容积卡尔曼滤波器，实现状态向量、模型参数的联合估计；
Zhu等[16]利用Student-T分布描述非高斯噪声模型，采用变分学习算法对系统模型和噪声参数进行联合估计。在此基础上，Huang等[17]将Student-T 分布加权积分分解为球面积分和径向积分，并利用球面径向变换规则进行积分数值化求解，提出基于Student-T分布的随机容积滤波器。然而在实际目标跟踪过程中，受传感器设备老化、故障或电磁干扰等因素影响，导致量测噪声出现不对称重尾性，利用混合高斯模型，Student-T分布模型难以对量测噪声进行准确表述，不准确的模型描述将严重降低非线性目标跟踪系统的滤波性能。Wang等[18]用Skew-T模型对非对称量测噪声进行建模，提出基于Skew-T分布扩展卡尔曼滤波算法(Skew-T extended kalman filter, SkewEKF)，假设噪声参数已知，一旦先验噪声信息与实际值存在较大偏差，难以保证准确的滤波效果。
因此，本文针对不确定重尾量测噪声条件下非线性系统状态估计问题，提出基于变分推理的鲁棒容积卡尔曼滤波(variational inference robust cubature kalman filter, VIRCKF)算法，该方法运用Skew-T分布对不对称重尾噪声进行表示，然后基于三阶球面径向规则计算状态和量测的容积点，并通过采样所得的容积点逼近系统状态和协方差，同时结合变分推理学习估计得到近似后验分布，并对系统状态进行更新，从而获得准确的目标状态。
考虑不对称重尾量测噪声下的目标跟踪问题, 建立系统模型：
xk=Axk-1+wk-1
(1)
zk=h(xk)+ek
(2)
式中：xk为k时刻系统状态向量；
A为状态转移矩阵；
wk-1为过程噪声，假设其服从零均值、协方差为Qk-1的高斯分布；
h(·)为量测方程；
zk是k时刻的传感器接收到的量测值；
ek为不对称重尾量测噪声，且Cov(wk,ek)=0。
(3)
式中：ST(·)表示Skew-T分布；
[ek]i表示k时刻噪声的第i个元素；
μ、σ和δ分别是位置参数，扩展参数和模型偏斜系数，v表示自由度，σ、δ、v>0。
以一维Skew-T分布为例，当i=1时，式(3)可表示为[19]：
(4)
且
(5)
(6)
式中：t(·;μ,σ2,v)为Student-T分布的概率密度函数；
T(·;0,1,v)是Student-T分布的累积分布函数；
Γ(·)表示伽马函数。
Skew-T分布、Student-T分布和高斯分布的概率密度函数如图1所示。当δ趋于0时，Skew-T分布接近Student-T分布；
当δ趋于0且v趋于∞时，Skew-T分布接近正态分布。由图1可以看出，Skew-T分布偏斜更明显，更适合对不对称重尾噪声进行模拟。
图1 3种分布的概率密度函数
2.1 变分推断
假设系统状态向量先验服从高斯分布，则式(1)表示为p(xk|xk|k-1)=N(xk;Axk|k-1,Qk-1)，系统一步预测的概率密度函数p(xk|z1:k-1)为：
p(xk|z1:k-1)=N(xk;xk|k-1,Pk|k-1)
(7)
结合式(2)和(3)得到量测似然函数：
p(zk|xk)∶ST(zk;h(xk),R,Δ,v)
(8)
由于式(8)中似然函数为Skew-T分布，其不满足贝叶斯推断中参数属于共轭指数分布簇的假设条件[20]，造成参数难以进行求解。因此，将式(8)的似然函数表示为：
(9)
(10)
(11)
式中：NH(·)为半正态分布，G(·)为伽马分布。
根据贝叶斯规则[21]，参数集合的后验概率密度函数计算为：
(12)
式(12)的难点在于边缘似然函数p(z1:k)的计算，故引入变分贝叶斯理论，考虑寻找一个简单近似分布q(xk,uk,Λk)去逼近真实后验分布p(xk,uk,Λk|z1:k)，取p(z1:k)的对数形式：
KL(q(xk,uk,Λk)‖p(xk,uk,Λk|z1:k))
(13)
KL(q(xk,uk,Λk)‖p(xk,uk,Λk|z1:k)) 为KL散度，其用来衡量近似分布q(xk,uk,Λk)与真实后验分布p(xk,uk,Λk|z1:k)之间的差异[22]，由于KL散度非负，只需令变分下界F(q(xk,uk,Λk))达到最大，即可获得后验概率最优近似分布。
系统联合概率密度函数为：
(14)
根据均值域逼近理论[23]，后验概率密度函数的近似因子分解形式为：
p(xk,uk,Λk|z1:k)≈qx(xk)qu(uk)qΛ(Λk)
(15)
2.2 变分贝叶斯容积卡尔曼滤波
非线性滤波问题可以归结为一个积分计算问题，被积函数为非线性函数和密度乘积的形式。容积卡尔曼滤波算法是将积分转换为球面径向形式，采用球面径向规则进行数值积分。由于包含Skew-T分布的后验概率密度函数形式复杂，考虑将变分推断引入CKF算法。在状态预测阶段，因状态向量服从高斯分布，可用三阶球面径向规则计算标准的高斯加权积分，即通过一组容积点数值化计算系统的状态均值和协方差。
(16)
ϑi,k-1|k-1=Sk-1|k-1ξi+xk-1|k-1
(17)
依据式(1)和式(7)，状态向量的先验均值xk|k-1和协方差矩阵Pk|k-1可表示为：
(18)
(19)
对容积点ϑi,k-1|k-1进行更新，
ϑi,k|k-1=Aϑi,k-1|k-1
(20)
则状态向量先验均值和协方差矩阵可表示为：
(21)
(22)
在量测更新阶段，量测服从Skew-T分布，通过变分推断计算其近似形式，用更新的容积点计算系统量测均值和协方差。其量测容积点计算如下：
(23)
ϑi,k|k-1=Sk|k-1ξi+xk|k-1
(24)
结合式(7)和(8)，量测向量的先验概率密度p(zk|z1:k-1)为：
(25)
量测向量先验均值zk|k-1和协方差矩阵Pzz,k|k-1表示为：
(26)
(27)
状态向量和量测的互协方差矩阵Pxz,k|k-1为:
(28)
其中：
p(xk,zk|z1:k-1)=p(zk|xk)p(xk|z1:k-1)=
ST(zk;h(xk),R,Δ,v)N(xk;xk|k-1,Pk|k-1)
(29)
考虑Skew-T分布用其近似分布式(9)代替，式(27)～(28)可以表示为式(32)～(33) ：
zi,k|k-1=h(ϑi,k|k-1)
(30)
(31)
(32)
(33)
系统状态的后验概率密度函数为：
p(xk,uk,Λk|z1:k)=
p(xk|z1:k-1)p(Λk)·p(uk|Λk)p(zk|xk,uk,Λk)
(34)
由于后验概率密度函数形式复杂，根据式(15)，通过求解logqx(xk)、logqu(uk)和logqΛ(Λk)可得到系统的状态和协方差，具体计算过程如下：
logqx(xk)=
(35)
由于变分贝叶斯近似适用于所有共轭指数域中的分布函数，那么后验可由其先验唯一确定[24]，所以其后验概率密度函数的形式为：
qx(xk)=N(xk;xk|k,Pk|k)
(36)
式中：更新后的卡尔曼滤波增益Kk、系统状态xk|k和协方差矩阵Pk|k分别为：
(37)
xk|k=xk|k-1+Kk(zk-k|k-1-Δūk)
(38)
(39)
同理，对于隐变量uk，有:
logqu(uk)=
(40)
其后验概率密度函数的形式为：
qu(uk)=NT(uk;uk|k,Uk|k)
(41)
式中：NT(·)表示截断正态分布。qu(uk)的后验概率密度参数的更新方程为：
uk|k=Ku(zk-h(xk|k))
(42)
(43)
(44)
对于似然函数参数Λk，有：
logqΛ(Λk)
=
(45)
其后验概率密度函数的形式为：
(46)
其中：
φk=R-1((zk-h(xk|k))(zk-h(xk|k))T+Pzz,k|k)+
(47)
式中：式(35)、(40)、(45)中的cx、cu和cΛ是与变量xk、uk和Λk相关的常数。
为了验证所提算法的有效性，考虑时变Skew-T分布噪声模型，分别采用CKF、SkewEKF以及VIRCKF算法对目标跟踪系统进行仿真实验。本文主要通过3种算法的均方根误差(RMSE)来说明算法的滤波性能。令目标运动模型为匀速模型，目标运动方程为：
xk=Axk-1+wk-1
(48)
(49)
式中：采样时间T=1 s，过程噪声协方差矩阵Q=diag([5×10-2,5×10-3,5×10-2,5×10-3])。
量测方程为：
(50)
式中：zk=[rkθk]T，rk和θk分别表示k时刻目标径向距离与方位角。量测噪声参数设置如下：
假设量测噪声协方差初始值R0=I2，令目标的初始位置为x0|0=[500 m,25 m/s,500 m,-20 m/s]T，状态协方差矩阵初始值P0=diag([2,10-2,2,10-2])。仿真时间为200 s。对3种算法分别进行100次Monte Carlo实验。
图2、图3分别为3种算法在X、Y方向的位置和速度均方根误差。从图2可以看出，当t>100 s时，由于量测噪声的突然增大，CKF算法的均方根误差出现波动；
当t>150 s时，量测噪声再次增大，CKF算法的RMSE也随之增大。SkewEKF算法的RMSE明显小于CKF，因为SkewEKF算法的本质是利用泰勒级数展开将非线性系统进行线性化，存在一定误差。而本文所提算法通过构建近似分布，迭代计算，直至逼近系统真实后验概率密度函数，具有自适应性。从图3可以看出，在150 s噪声增大时，3种算法速度RMSE均增大，但VIRCKF算法RMSE更小。
图2 不同方向的位置均方根误差
图3 不同方向的速度均方根误差
通过计算VIRCKF算法的变分下界函数，对算法的收敛性进行分析。图4为实验过程中系统状态的变分下界函数。可以看出，该算法在X方向位置变分下界稳定在[-3.40,-3.00]，在Y方向最终稳定在[-4.41,-3.65]，表明算法稳定收敛。
图4 变分下界函数
图5为量测噪声的近似分布拟合真实分布的过程，可以看出，随着参数的迭代优化，近似分布会不断逼近真实分布。
图5 量测噪声的真实分布与近似分布的拟合过程
图6 变分下界函数随值的变化
图7 变分下界函数随值的变化
分别对3种算法的计算复杂度进行分析，CKF算法的时间复杂度为O(6Nn)，SkewEKF的时间复杂度为O(Nm)，VIRCKF算法的时间复杂度为O(8Nm+6Nnm+10N)。其中，N为总运行时间，2n表示容积点数，m为变分迭代次数。虽然VIRCKF的时间复杂度更高，但综合3种算法的滤波效果和鲁棒性分析，本文所提的VIRCKF算法滤波精度更高，性能更加优越。
本文针对不确定重尾量测噪声条件下非线性系统目标跟踪问题，提出基于变分推理的鲁棒容积卡尔曼滤波算法。将Skew-T分布的概率密度函数用近似形式代替，然后利用球面径向规则对函数进行数值积分求解，同时结合变分推理学习估计得到系统的近似联合后验分布。实验结果表明，在量测噪声为重尾分布且参数不确定时，本文算法具有较高的估计精度和较好的滤波效果。未来可以探索更有效的估计算法，解决不确定过程噪声与量测噪声下的机动目标跟踪问题。
猜你喜欢概率密度函数变分后验求解变分不等式和不动点问题的公共元的修正次梯度外梯度算法数学物理学报(2022年5期)2022-10-09幂分布的有效估计*广西民族大学学报(自然科学版)(2022年1期)2022-05-18一类传输问题的自适应FEM-BEM方法山西大学学报（自然科学版）(2021年4期)2021-08-31一种基于折扣因子D的贝叶斯方法在MRCT中的应用研究*中国卫生统计(2020年3期)2020-06-28基于贝叶斯理论的云模型参数估计研究统计与决策(2019年6期)2019-04-22自反巴拿赫空间中方向扰动的广义混合变分不等式的可解性Acta Mathematica Scientia(English Series)(2018年6期)2018-03-01已知f（x）如何求F（x）当代旅游(2018年8期)2018-02-19基于变构模型的概率密度函数的教学探索数学学习与研究(2018年2期)2018-02-09一种基于最大后验框架的聚类分析多基线干涉SAR高度重建算法雷达学报(2017年6期)2017-03-26基于变分水平集方法的数字图像分割研究中国市场(2016年45期)2016-05-17本文来源：http://www.triumph-cn.com/fanwendaquan/gongwenfanwen/2023/0822/99386.html