基于自然梯度的噪声自适应变分贝叶斯滤波算法

2023-10-30 10:13胡玉梅胡振涛

自动化学报 2023年10期

胡玉梅潘泉胡振涛郭振,5

非线性滤波器设计与实现因其普适性及重要性一直是国内外学者研究的热点问题,近年来非线性状态估计理论已成功应用于陆、海、空、天中运动目标的预警与防御,智能交通的精确导航与制导,无人机定位与遥感监测、工业过程监控与故障诊断等众多领域[1-3].考虑到运动目标跟踪系统机动、隐身等复杂多变的人为对抗特征以及非视距、干扰、遮挡等环境因素无法避免,其系统建模、估计与辨识过程中越来越无法回避非线性、非高斯以及参数未知等复杂系统特征的影响[4].在对复杂环境下运动目标系统噪声先验信息进行建模时,建模误差存在于状态演化模型中,并通常假设其满足一定的参数化统计特性[5].然而,在实际工程中,由于先验建模信息的不足导致难以对此类参数进行准确赋值.例如,在现代目标跟踪系统中出现的欺骗、干扰、杂波、未知分布特性的量测噪声和系统噪声等情形,尤其在非合作运动目标强机动场景中,因难以对其运动过程进行精细化建模,常造成目标跟踪航迹间断甚至无法正常起始的现象.

针对系统噪声未知问题,其解决思路通常选取逆伽马分布和逆威沙特(Inverse-Wishart,IW)分布等具有统一参数表达形式的指数族分布作为共轭先验分布.Särkkä 等[6]在变分贝叶斯 (Variational Bayes,VB)推断框架下利用指数族分布构建共轭先验分布近似未知量测噪声的后验概率密度函数(Probability density function,PDF),进而结合贝叶斯滤波机理实现时变噪声方差和目标状态的联合估计,其滤波效果达到与交互式多模型方法相接近的估计精度,并且凭借 VB 便于结合平均场近似解耦理论将高维变量求解转化为多个低维变量计算的特点,使其在解决多未知扰动问题方面更具有优势.随后,文献[7] 采用IW 分布近似多变量噪声后验PDF 的思想,给出变分推断框架下噪声方差估计的一般实现形式.在此基础上,变分贝叶斯方法在未知量测/系统噪声估计方面得到了进一步的发展[8-9],并成功推广应用于多目标跟踪环境[10-11].在文献[12-13]中滤波器的构建均建立在量测噪声和系统噪声统计特性未知的假设条件下,其中,文献[12] 在状态和量测扩维基础上通过采用批处理方式实现对未知噪声矩阵的估计,然而这种扩维方式不可避免造成计算复杂度的急剧增加;文献[13] 给出两类噪声统计特性均未知条件下噪声方差的在线估计方法,并且为避免系统噪声和状态相互独立假设的不合理性,采用独立于当前状态的前一时刻状态预测误差协方差代替系统噪声的统计特性,进而利用平均场理论近似解耦计算边缘后验 PDF 的变分分布,以迭代递推方式求解其估计变量期望的解析解,同时,将状态一步预测误差协方差的先验分布建模为参数化IW 分布,以实现状态后验概率分布的更新.

考虑到外界环境干扰以及非线性传播等因素的影响,量测数据概率分布往往呈现出重尾和非对称等非高斯特征.例如,在基于电磁信号的距离估计中,障碍物遮挡造成的非视距量测误差往往远大于其他来源服从对称分布的误差,导致量测分布呈现非对称非高斯现象.其次,受传感器精度和灵敏度的限制,当运动目标发生强机动时,雷达极化反射不稳定性将导致的量测随机缺失现象,也将造成量测产生非高斯特征.针对量测噪声非高斯问题的处理,文献[14]采用莱斯分布构建非共轭指数族变换点检测模型,并将其成功应用于雷达目标跟踪系统中.文献[15]则提出一种通过构建状态演化模型和量测模型的条件矩实现非线性非高斯系统的状态估计方法.为了进一步提升滤波器对不同分布形状噪声的鲁棒性,文献[16] 将 Skew-t分布作为非高斯量测似然的近似分布形式,在此基础上设计了一种鲁棒变分贝叶斯估计器.文献[17] 考虑在α散度最小化准则下,利用变分分布实现对后验 PDF 的近似.α散度对异常数据具有较好的抑制作用,但是也打破了对数边缘概率与变分置信下界 (Evidence lower bound,ELBO)、KL 散度(Kullback-Leibler divergence)之和的等式约束关系,其实质是一种伪 VB 推断方法.考虑到学生t分布和 Skew-t分布对因量测异常导致的量测重尾现象和非对称性具有较好鲁棒性[18-19],文献[20] 针对系统噪声和量测噪声均具有重尾特性的线性系统采用学生t分布对状态预测概率密度函数和量测似然函数进行建模,进而提出了一种基于高斯-学生t混合分布的线性滤波器,与传统高斯假设条件下滤波器估计效果相比进一步提升了系统的状态估计的精度和鲁棒性.在系统噪声未知且时变和量测噪声重尾条件下,文献[21] 构建参数化IW 分布作为状态一步预测误差协方差的共轭先验分布,同时,选取参数化学生t分布刻画具有厚尾特点的量测似然函数.然而,上述噪声自适应方法未考虑非线性滤波器自身的优化问题和如何衡量后验分布近似程度的问题.

考虑到在非线性滤波器设计过程中,参数化变分分布对系统状态后验近似的程度是提高估计精度的关键因素之一,基于采样的随机优化和基于拟牛顿、高斯牛顿和梯度上升等确定性优化方法的非线性滤波器相继提出.在随机性优化方面,MCMC(Markov chain Monte Carlo)[22-23]和序贯蒙特卡罗(Sequential Monte Carlo,SMC)利用随机样本来逼近后验概率,并以样本分布近似未知变量后验分布.随机优化方法的优点是在大量样本的条件下能够保证较高精度的估计,但要承担较大的计算负担[24].在确定性优化方面,文献[25] 以最小二乘准则为目标函数,采用牛顿法推导给出一种迭代卡尔曼滤波器.文献[26] 综合梯度法和牛顿法提出一种无噪声条件下的滚动时域估计器,并证明其渐进收敛性和稳定性.文献[27] 针对非线性状态空间的状态估计及其估计误差协方差的迭代更新问题,利用正则化的非线性最小二乘思想提出一种随机增量近端梯度算法.VB 方法通过求解参数化的目标函数,利用参数优化结果逼近后验分布,是一种确定性近似方法.因此,VB 具有确定性优化方法特有的计算量较小的优势,同时由于 VB 便于结合平均场理论近似解耦联合概率分布,进一步减小了计算代价.

非线性状态估计优化的实质是对多维状态后验PDF 的近似逼近,并且其近似程度不能简单地采用欧氏距离进行度量.因此,选取合理度量准则将有利于提高后验 PDF 的近似程度,KL 散度作为分布之间“差异”的度量在后验分布近似性衡量中具有天然优势.文献[17]给出状态估计的变分迭代优化实现形式,获得α散度和 KL 散度下对后验 PDF 更紧密的近似.同时,从信息几何角度出发,概率分布是统计流形上点,在一定条件下两概率分布之间的 KL 散度与作为统计流形度规的 Fisher 信息满足一定的数学关系.基于信息几何理论,Amari[28]利用自然梯度优化方式实现统计流形空间中目标函数的最速梯度下降(或上升).文献[29] 结合自然梯度策略和卡尔曼滤波框架设计一种非线性状态估计方法,在文献[30] 中,该方法进一步推广应用于传感器网络目标跟踪系统.文献[31] 中作者证明了针对非线性状态估计优化的自然梯度方法在克拉美罗下界意义下是渐进最优的.考虑自然梯度优化的优势,以变分置信下界最大化条件下状态估计及其估计误差协方差的自然梯度为切入点,从信息几何角度实现对状态后验概率密度函数的“紧密”逼近,进而提高状态估计精度.

在过程噪声先验不准确和由于量测随机丢失导致的量测噪声分布非高斯的情况下,针对非线性系统状态估计精度和鲁棒性提升问题,本文提出一种基于自然梯度的噪声自适应变分贝叶斯滤波算法.本文结构如下:第 1 节结合IW 分布和学生t分布分别实现对状态预测误差协方差和量测似然的参数化建模;第 2 节结合平均场近似和坐标上升方法给出了变分贝叶斯框架下未知变量变分分布的迭代更新方式;第 3 节推导给出 ELBO 关于系统状态估计及其误差协方差的自然梯度,进而构建并设计一种基于自然梯度的噪声自适应非线性变分贝叶斯滤波器;第 4 节给出仿真验证与分析;第 5 节总结全文,并展望了后续研究方向.

1 预备知识

假设动态系统隐变量Ξk的量测为zk,根据变分贝叶斯理论可知,在变分贝叶斯框架下以变分分布q(Ξk|ψk)作为桥梁可将难以积分的问题转化为ELBO 的优化问题.

其中,L(ψk)表示具有单调递增特性的变分 ELBO,DKL[q(Ξk|ψk)‖p(Ξk|zk)]表示变分分布q(Ξk|ψk)与后验分布p(Ξk|zk)之间的 KL 散度,其表达式分别为

其中,q(Ξk|ψk)表示以ψk为参数的变分分布.

未知变量的估计实际是对状态后验分布p(Ξk|zk)的近似逼近,当非线性状态后验分布难以获取时,变分贝叶斯方法能够通过构建简单的变分分布实现对p(Ξk|zk)的近似逼近.变分分布选取需考虑先验分布和后验分布具有相同的函数形式,一般情况下具有共轭特性的指数族分布满足这一条件.本文考虑系统噪声未知情况,根据标准卡尔曼滤波实现结构可知,系统噪声的统计特性仅影响状态预测误差协方差,因此可直接对状态预测误差协方差进行先验建模.并且当系统噪声假设服从均值已知但方差未知的多变量高斯分布假设时,可选取IW 分布作为其分布方差矩阵的共轭先验分布;此外,针对量测缺失导致量测噪声出现重尾问题,考虑选取能够有效表征这一现象的学生t分布对量测似然函数进行建模.

1.1 IW 分布共轭先验模型

在贝叶斯滤波框架下,对于系统噪声高斯分布参数未知问题的处理,可转化为状态预测误差协方差的估计问题.这里选取IW 分布作为共轭先验分布,以保证后验分布与先验分布具有相同的函数形式.一个服从IW 分布的n×n维随机对称正定矩阵A的概率密度函数可以表示为

其中,t和T分别表示IW 分布的自由度参数和n×n维逆尺度矩阵,tr(·)表示矩阵的迹,Γn(·)表示n维的伽马函数.进而状态预测误差协方差Pk|k-1的先验分布表示为

对于逆威沙特分布,当tk≥n+1 时,其变量期望即状态预测误差协方差的均值E[Pk|k-1]表示为

令先验分布参数tk表示为

式 (6)进一步转化为

其中,τ表示调谐参数[7].根据式 (5)～(8)计算得到状态预测协方差的先验分布及其分布参数.

1.2 学生t分布量测似然模型

当量测噪声分布由于量测值随机异常或缺失呈现重尾特征时,传统高斯噪声分布模型难以对噪声分布特性进行有效刻画.考虑学生t分布能够更好地体现这一特征,并且对异常值具有较好的鲁棒性.当量测噪声υk满足均值为0、方差为Rk的分布参数时,建立参数化学生t分布模型,即

其中,v表示学生t分布的自由度参数,噪声方差Rk是该分布的尺度矩阵.在此条件下量测似然函数可表示为

通常学生t分布的概率密度函数难以求得封闭解,为了解决这个问题,引入服从伽马分布的辅助随机变量λk,从而进一步将学生t分布转化为服从参数化高斯分布的伽马分布的积分,式 (10)中量测似然函数可转化为

其中,G(λk;αk,βk)表示辅助变量λk服从伽马分布,αk和βk分别为形状参数和逆尺度参数[20].综合上述特点,在传感器量测存在随机异常值情况下量测似然函数表示为如下结构化形式:

并且G(λk;αk,βk)及其均值E[λk]分别表示为

最后,通过对辅助变量λk及其超参数αk和βk更新确定量测似然的表达式.

1.3 平均场近似

根据以上分析建模,式 (2)中的系统隐变量Ξk定义为Ξk=(xk,Pk|k-1,λk).考虑到联合后验分布形式的复杂性以及参数相互耦合的问题,无法直接求得其解析解.在假设各未知变量相互独立的前提下,通过引入平均场近似策略实现联合变分分布的近似解耦,即

其中,q(xk),q(Pk|k-1)和q(λk)分别表示状态xk,状态一步预测协方差Pk|k-1和辅助变量λk的变分分布,通过平均场理论对待估计变量联合分布进行分解,有助于在变分贝叶斯迭代框架结合坐标上升和各类优化方法实现多变量的联合优化.结合上述第1.1 节和第1.2 节中的共轭先验模型和量测似然模型,相应的变分分布分别表示为

2 未知变量的变分分布更新

根据上述构建的共轭先验模型和量测似然模型,变分置信下界中隐变量Ξk包括状态xk,状态一步预测协方差Pk|k-1和辅助变量λk,Ξk可定义为Ξk=(xk,Pk|k-1,λk).采用变分分布近似上述多未知变量的联合后验分布,即

进而,相应变分置信下界转化为

在变分置信下界最大化约束下,结合坐标上升方法,变量xk,Pk|k-1和λk概率分布的对数形式的更新分别表示为

结合建模信息和各变量变分分布的特点,联合后验分布及其对数形式分别表示为

其中,cΞk表示迭代过程中与变量相关的实常数,Dz和n分别表示向量zk和矩阵Pk|k-1的维数,(xk-表示xk|k-1)的缩写形式(下文中类似情况不再一一说明).根据式(26)可知,由于高斯分布、学生t分布和伽马分布均属于指数分布族而具有简单的对数形式,因此便于采用迭代更新的方式计算其解析解.以下给出λk、Pk|k-1以及xk的迭代估计的具体实现原理和步骤.

2.1 λk迭代更新

根据式 (22)和式 (26),在第i+1 次变分迭代更新过程中,λk的变分分布对数形式可表示为

其中,cλk表示xk和Pk|k-1积分后相关的实常数.根据贝叶斯无偏估计理论,,其中和分别表示第i次估计值及其估计偏差.同时,结合局部线性化技术,式 (27)中期望部分可进一步改写为

计算式 (14)中伽马分布的对数表达式,并结合式 (27)和式 (28),不难看出,第i+1 次更新后参数的变分分布中超参数的更新过程可表示为

根据式 (15)中伽马分布的均值表达式,同时结合式 (29)和式 (30),参数λk的第i+1 次变分迭代更新可表示为

由式 (31)可以清晰地看出,辅助随机变量λk的更新与负相关,并且在迭代过程中的变化与上一次迭代的状态估计值及量测函数在此估计值处的局部线性化矩阵有关.根据表达式,第i次迭代的估计值(可理解为第i+1 次迭代估计的先验)越接近真实状态,值越小,根据式 (31)可知则越大;相反,第i次的估计状态距真实状态较远,值越大,则越小.其物理意义是,当发生量测随机缺失时导致滤波更新的新息增大,同时值将随之增大,从而造成式 (11)中实际量测噪声方差增大,这种情况将反馈于滤波过程中状态估计的自适应更新.为了直观地解释这种现象,结合式 (11)和式 (13),采用一维高斯分布和伽马分布说明其相关关系.图1 给出伽马分布参数对量测似然的影响示意图.在图1 中,假设高斯分布的均值和方差分别为25 和2/E[λ],并且E[λ]=α/β,伽马分布参数α和β的值分别在图1(a)和图1(b)中注明.当量测未发生丢失时,β的值较小,高斯分布比较集中(如图1(a)所示);当量测发生随机丢失,β的值增大,高斯分布比较分散(如图1(b)所示).

图1 伽马分布参数对量测似然函数的影响示意图Fig.1 The diagram of the influence of Gamma distribution parameters on likelihood

2.2 Pk|k-1迭代更新

假设第i次迭代状态估计值已知,根据式 (23)和式 (26),在第i+1 次变分贝叶斯迭代更新过程中,Pk|k-1的变分分布的对数形式表示为

因此,根据式 (6)中Pk|k-1均值表达式,同时结合式 (34)和式 (35),状态预测误差协方差的第i+1 次变分迭代更新表示为

根据卡尔曼滤波实现原理可知,当过程噪声方差建模不准确时,状态预测误差协方差较大,并且导致量测预测不精确而产生较大的量测新息zkhk(xk|k-1).同时,由式 (33)可知,的迭代更新与量测新息正相关,进而从物理意义上说明状态预测协方差建模的合理性.

2.3 xk|k迭代更新

根据式 (24)和式 (26),在第i+1 次迭代更新过程中,系统状态的变分分布对数形式表示为

其中,cxk表示与Pk|k-1和λk积分后相关性的实常数.在第i+1 次迭代中,依据式 (31)和式 (36)已经获取状态预测误差协方差Pk|k-1和辅助变量λk的估计值

由贝叶斯滤波理论可知,状态后验分布的近似变分分布应具有以下形式

结合高斯分布表达式,式 (40)可进一步改写为

综合式 (37)～(42),目标状态后验分布的近似变分分布分别是以为期望和协方差的高斯分布,即

针对非线性系统状态估计问题,现有贝叶斯滤波方法在实现式 (43)中状态估计及其误差协方差的迭代更新过程中,往往难以获得与真实后验“紧密”近似的变分分布.虽然在多个变量坐标迭代优化过程中能够弥补部分由于线性化和参数不精确带来的误差,但估计精度提升效果受限,有时难以满足实际工程需求.因此,考虑在获得λi+1和估计值的基础上,在置信下界最大化条件下结合自然梯度方法和 Fisher 信息矩阵,推导给出迭代更新解析表达式的具体实现.

3 基于自然梯度的噪声自适应变分贝叶斯滤波器

3.1 变分置信下界的自然梯度

第i次变分迭代后,在分别获取参数λk和Pk|k-1的估计值的基础上,式 (2)可转化为仅包含未知变量xk的形式

定义此时变分参数ψk=(xk|k,Pk|k),通过最大化式(44)中的置信下界即可获得状态估计及其误差协方差.自然梯度通过信息几何方法寻找统计流形空间的目标函数最速下降/上升方向,给出了一种实现式 (44)置信下界求解的可行性方案[25].下面以置信下界L(ψk)最大化为目标函数,同时结合Fisher 信息矩阵推导关于变分参数ψk的自然梯度.式 (44)的置信下界可以进一步表示为

假设第i次迭代是第i+1 次迭代的先验,即式(45中的先验信息p(xk)可以近似为,则变分参数ψk的最优估计值表示为

通过计算式 (46)等号右边关于ψk线性化函数,可获取置信下界的自然梯度以及相应的最优.

定理 1.假设系统状态演化满足高斯分布特性,并且对数似然期望Eq[logp(zk|xk)]在的邻域内一阶可导,同时,在的邻域内二阶可导,则第i+1 次迭代过程中变分置信下界的自然梯度为

证明.根据变分参数的最优表达式 (46)可知,通过对其求关于ψk的偏导数即可获得置信下界最大化条件下的极值点.由于对数似然期望一阶可导,令,根据泰勒展开上式对数似然期望线性化后可以近似表示为

其中,Fψk为具有半正定性 Fisher 信息矩阵,即

此时,式 (46)进一步转化为

计算式(51)关于ψk的偏导,并令其取值为零,可获得在第i+1 次变分迭代中ELBO 关于ψk的自然梯度

因此,在置信下界最大化条件下变分参数ψk的更新表达式为

需要说明的是,VB 迭代更新等价于在统计流形空间沿以 Fisher 信息矩阵作为权重的自然梯度方向移动,从而获得后验分布的最佳近似,并且具有以下性质:

2)Δψk →0 时,式 (46)中的第2 项KL 散度表达式转化为式 (49),因此从局部意义上迭代过程中的变分分布与后验分布的近似误差是正定二次型.

3)自然梯度是Fisher 有效的[25],因此上述优化过程是 Fisher 有效的估计器.特别地,上述ψk的迭代估计是无偏的.

4)在自然梯度的推导过程中,在无限小的邻域内 (即Δψk →0 )给出 Fisher 信息矩阵与 KL 散度之间的近似关系,并且的计算依赖于参数化的变分分布q(xk|ψk),因此,随着变分迭代的进行,Fihser 信息不断变化.

3.2 基于自然梯度的滤波更新

系统状态的变分分布包含两个超参数,即当前时刻状态估计xk|k及其误差协方差Pk|k,为实现其迭代更新需分别计算相应的自然梯度.定理 2 及其证明给出了变分置信下界分别关于xk|k和Pk|k的自然梯度及其推导过程.

定理 2.假设第i次迭代更新的状态估计及其估计误差协方差分别为,在满足定理 1 的条件下,对数似然期望Eq[logp(zk|xk)]分别在xk=的邻域内一阶可导;同时,KL散度分别在和的邻域内二阶可导;并且在已经分别获取第i次迭代状态预测误差协方差以及辅助变量估计值的前提下,变分置信下界关于xk和Pk|k的自然梯度分别表示为

证明.首先推导在定理 2 中假设条件下 KL 散度关于xk|k和Pk|k的二阶导数.若高斯分布q1～N(ξ1;µ1,C1)和q2～N(ξ2;µ2,C2)具有相同维数d,则两分布之间的KL 散度表示为

因此,迭代过程中变分分布之间的 KL 散度表示为

结合式 (50),关于xk|k的 Fisher 信息矩阵表示为

关于Pk|k的 Fisher 信息矩阵的推导如下:

其中,哈达玛积◦表示取矩阵的对角线元素构成单位阵.忽略误差协方差矩阵非对角线元素的影响,式(58)可以近似简化为

由矩阵求导知识可知,矩阵X逆函数X→X-1的导数表示为,其中,m和p分别表示矩阵X的不同行,n和l分别表示其不同列.因此,在迭代优化过程中关于估计误差协方差的 Fisher 信息矩阵表示为

其中,符号⊗表示矩阵克罗内克积运算.

结合服从高斯分布的似然函数表达和局部线性化方法,式(61)可转化为

同时,根据 Price 定理可知

由于Eq[zk-hk(xk)]=0,局部线性化处理后,式(63)可表示为

综合上述过程,在第i+1 次迭代中变分置信下界关于xk|k和Pk|k的自然梯度分别表示为

结合式 (53)中变分参数迭代更新以及式 (66)和式(67)中自然梯度的具体表达式,第i+1 次迭代时状态估计值及其误差协方差表示为

综合上述式(68)和式 (69)的构建过程,本文给出了一种基于自然梯度的噪声自适应变分贝叶斯滤波算法,为便于理解算法整体实现流程,下面以伪代码的形式给出算法具体实现步骤,详见算法1.

算法1.基于自然梯度的噪声自适应变分贝叶斯滤波算法

由算法1 中状态估计以及估计误差协方差的更新过程可知,算法1 具有以下特点和优势.

3)自然梯度的方向是目标函数在统计流形空间中 ELBO 的最速上升方向,即 ELBO 沿基于Fisher 信息矩阵的自然梯度向使其增大的方向移动,以从信息几何角度获得状态后验分布的最优近似.

为实现迭代优化过程的自适应,采用以下相对估计误差er作为判断迭代是否终止的条件.

其中,ϵ表示一个较小的正实数,其取值大小可根据实际应用场景的精度需求合理设定.

4 实验验证

4.1 仿真算例

为验证本文提出算法的可行性和有效性,采用二维笛卡尔坐标系下三雷达量测的目标跟踪系统进I2表示2 维单位阵.初始估计值和其行500 次 Monte Carlo 仿真实验,目标做转弯率未知的匀速转弯运动,并假设过程噪声先验信息不准确且缓慢变化[13],同时三个雷达传感器均存在量测随机缺失现象.本文所提算法记为 VBAKF-NG,其三个传感器分布式估计误差协方差融合方法记为DVBAKF-NG.仿真实验中选取对比算法分别为:结合本文构建的逆威沙特分布的共轭先验分布模型和服从学生t分布的量测似然模型,分别利用扩展卡尔曼滤波器 (Extended Kalman filter,EKF)和无迹卡尔曼滤波 (Unscented Kalman filter,UKF)滤波器实现状态的预测与更新,其中变分扩展卡尔曼自适应滤波算法记为 EKF-VB,其对应的三个传感器分布式估计误差协方差融合方法记为DEKF-VB,变分无迹卡尔曼自适应滤波方法记为 UKF-VB,采用传统自适应结构的无迹卡尔曼滤波方法记为 UKF-FN,基于虚拟量测采样的集合卡尔曼滤波器记为 EnKF,其对应的三个传感器分布式估计误差协方差融合方法分别记为 DUKF-VB,DUKF-FN 和 DEnKF;同时,将采用真实量测噪声和过程噪声的 EKF 和 UKF 实现方式分别记为DEKF-TN 和 DUKF-TN.

4.1.1 场景设置

在k时刻目标状态表示为,其中,[xk yk]T和分别表示k时刻目标在水平方向和竖直方向的位置与速度,=-0.052 rad/s表示转弯率.相应地,线性化后的状态转移矩阵Fk表示为

其中,T=1 表示采样周期.过程噪声ωk随时间缓慢变化且先验知识不准确,假设系统噪声协方差的变化满足方程

其中,q1=0.3 m2/s2和q2=1×10-4m2/s2分别为噪声方差的参数,ts表示仿真的总时长,误差协方差分别为x0|0=[2 100 m,0 m/s,1 800 m,30 m/s,-π/50 rad/s]和P0|0=diag{150,10,150,10,0.0001}.直角坐标系内三个雷达的位置分别为(0 m,0 m)、(500 m,0 m)和(0 m,500 m),其量测方程表示为

仿真实验中,迭代次数Niter=5,传感器扫描周期T=1 s,其余仿真实验参数取值如表1 所示.

表1 仿真参数Table 1 Simulation parameters

4.1.2 结果与分析

图2 和图3 分别给出基于传感器 1 量测的状态预测误差协方差对角线上位置分量和速度分量的估计值,从图2 和图3 可以看出,基于本文建模方法的3 种滤波算法均能对预测误差协方差进行有效估计,避免了因过程噪声先验建模不准对状态估计精度造成的不利影响.同时,VBAKF-NG 在更新阶段采用自然梯度优化获得更高精度的估计和相应的误差协方差,并作为下一时刻的先验信息,因而其预测误差协方差最小,并且更有利于在时序上形成一个“当前时刻高精度估计⇄下一时刻高精度预测”的良性循环.

图2 传感器1 位置状态预测误差协方差的估计值Fig.2 The expectation of the position state prediction error covariance from Sensor 1

图3 传感器1 速度状态预测误差协方差的估计值Fig.3 The expectation of the velocity state prediction error covariance from Sensor 1

图4～7 分别给出传感器 2 和传感器 3 的状态预测误差协方差对角线上位置量和速度量的估计值,可以看出,基于传感器 2 和传感器 3 量测的3种算法对状态预测误差协方差估计的曲线与传感器1 相似,即VBAKF-NG 预测误差协方差位置分量和速度分量的估计值曲线始终处于最优的位置且更加平滑,体现出相对于 EKF-VB 和 UKF-VB 能够获取更高的估计精度和鲁棒性.

图4 传感器2 位置状态预测误差协方差的估计值Fig.4 The expectation of the position state prediction error covariance from Sensor 2

图5 传感器2 速度状态预测误差协方差的估计值Fig.5 The expectation of the velocity state prediction error covariance from Sensor 2

图6 传感器3 位置状态预测误差协方差的估计值Fig.6 The expectation of the position state prediction error covariance from Sensor 3

图7 传感器3 速度状态预测误差协方差的估计值Fig.7 The expectation of the velocity state prediction error covariance from Sensor 3

图8给出传感器 1 在径向距和方位角上的量测噪声方差.图8 中黑色跳变表示该传感器此时出现量测丢失情况,从图8 中可以看出,UKF-VB 和 VBAKF-NG 能够准确地识别出量测噪声方差的跳变,以便及时调整分布参数更新.同时,从图8 中可以清晰地看出,VBAKF-NG 对噪声方差跳变的识别精度更高.图9 和图10 分别给出传感器 2 和传感器 3 在径向距和方位角上的量测噪声方差.结果同样表明与 UKF-VB 和 EKF-VB 相比,VBAKFNG 对噪声方差跳变具有更高的识别精度.实际上,图2～7 中 VBAKF-NG 之所以在滤波精度和鲁棒性上优于 EKF-VB 和UKF-VB,其原因正是由于VBAKF-NG 可以有效辨识量测随机缺失现象产生,从而使得算法具有较好的容错特性.

图8 传感器1 量测噪声方差的估计值Fig.8 The expectation of the measurement noise variance from Sensor 1

图9 传感器2 量测噪声方差的估计值Fig.9 The expectation of the measurement noise variance from Sensor 2

图10 传感器3 量测噪声方差的估计值Fig.10 The expectation of the measurement noise variance from Sensor 3

图11和图12 分别给出上述6 种方法状态估计的径向距均方根误差 (Root mean squared error,RMSE)和径向速度均方根误差的对比,从图11 和图12 中可以看出,采用IW 分布和学生t分布分别建模预测误差协方差和量测噪声的 DEKF-VB,DUKF-VB 和 DVBAKF-NG 三种方法均优于传统自适应滤波算法 DUKF-FN 和集合卡尔曼滤波算法 DEnKF.同时,由于 DVBAKF-NG 在更新过程中采用自然梯度的方法更新目标状态,与 DEKFVB,DUKF-VB 相比,取得了相对更优的估计精度.需要说明的是,DUKF-TN 在滤波过程中采用真实的过程噪声,因此具有最高的精度.

图11 径向距估计RMSE 的对比Fig.11 The RMSE comparison of radical range estimation

图12 径向速度估计RMSE 的对比Fig.12 The RMSE comparison of radical velocity estimation

4.2 毫米波雷达数据验证

实验场景为采用 ARS408-21 毫米波雷达对地面沿直线运动的行人目标进行跟踪.在量测过程中考虑量测噪声非高斯特点,行人分别在 17 时刻和27 时刻被遮挡,导致出现两次量测丢失.同时考虑系统噪声方差的未知时变情况,行人的速度大小不断变化.在不考虑杂波和数据关联的情况下,假设一个时刻返回一个量测,量测数据和滤波结果如图13 所示.

图13 毫米波雷达地面行人跟踪实验Fig.13 Tracking a pedestrian by using a millimeter wave radar

需要说明的是,由于真实系统噪声未知,在此实验中采用真实系统噪声方差的 UKF-TN 算法未参与对比.从图13 中可以看出,由于系统噪声未知,在滤波初始阶段,EnKF 和 UKF-FN 算法的估计结果与量测偏离较大,采用IW 分布建模的 EKFVB、UKF-VB 和 VBAKF-NG 与量测偏离较小.当发生量测缺失时,对比算法均存在不同程度的发散现象,VBAKF-NG 能够较好地抑制滤波发散.因此,毫米波雷达地面行人跟踪实验验证了本文所提 VBAKF-NG 算法的优越性.

5 结束语

针对过程噪声时变且先验信息建模不准确和量测噪声非高斯问题,首先,引入参数化的IW 分布作为状态预测误差协方差矩阵的共轭分布建立先验信息模型,同时构建参数化学生t分布模型以近似由于量测随机丢失造成的非高斯量测似然函数;继而,结合平均场理论近似解耦状态和多参数的联合概率密度函数以获取相应的边缘变分分布,并采用坐标上升的迭代方法更新状态和各参数的边缘变分分布;在此基础上,结合 Fisher 信息矩阵推导给出置信下界关于状态估计及其误差协方差的自然梯度,进而给出一种基于自然梯度的噪声自适应变分贝叶斯滤波算法,实现过程噪声时变和由于量测随机缺失导致的量测噪声非高斯条件下非线性系统的高精度估计.本文算法处理的对象局限在单目标跟踪系统,如何将其拓展到机动多目标跟踪系统,在变分贝叶斯框架下,结合自然梯度策略有效地解决模型自适应选取以及多目标跟踪中数据关联问题是下一步将开展的工作.