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

基于自然梯度的噪声自适应变分贝叶斯滤波算法胡玉梅1,2,3潘泉1,2胡振涛4郭振1,2,5摘 要考虑到运动目标跟踪系统机动、隐身等人为对抗特征以及非视距、干扰、遮挡等环境因素,其系统建模、估计与辨识过程中越来越无法回避非线性、非高斯以及参数未知等复杂系统特征的影响.针对过程噪声先验信息不准确以及量测噪声非高斯环境下运动目标的非线性状态估计问题,提出一种基于自然梯度的噪声自适应变分贝叶斯(VariationalBayes,VB)滤波算法.首先,利用指数族分布具有统一表达形式的优势,构建参数化逆威沙特(Inverse-Wishart,IW)分布作为状态一步预测误差协方差的共轭先验分布,同时选取学生 t 分布重构因量测随机缺失导致的具有非高斯特点的似然函数;其次,在变分贝叶斯优化框架下采用平均场理论将状态变量联合后验分布近似分解为独立的变分分布,在此基础上,结合坐标上升方法更新各变量的变分分布参数;进而,结合Fisher信息矩阵推导置信下界最大化关于状态估计及其估计误差协方差的自然梯度,使非线性状态后验分布的近似分布沿梯度下降,以实现对状态后验概率密度函数(Probabilitydensityfunction,PDF)的“紧密”逼近.理论分析和仿真实验表明:相对传统的非线性滤波方法,本文算法对噪声不确定问题具有较好的自适应能力,并且能够获得较高的状态估计精度.关键词非线性滤波,自适应滤波,变分贝叶斯推断,自然梯度,Fisher信息矩阵引用格式胡玉梅,潘泉,胡振涛,郭振.基于自然梯度的噪声自适应变分贝叶斯滤波算法.自动化学报,2023,49(10):20942108DOI10.16383/j.aas.c210964A Novel Noise Adaptive Variational Bayesian Filter Using Natural GradientHUYu-Mei1,2,3PANQuan1,2HUZhen-Tao4GUOZhen1,2,5AbstractConsideringtheincreasingcomplexityandchangeabilityofcharacteristicssuchasmaneuveringandstealthinmovingtargettrackingsystemandtheinfluenceofadversefactorssuchasnon-line-of-sight,interferenceandocclusioninmeasurementenvironment.Stateestimationislikelytobeconfrontedwithcomplexsystemcharac-teristicssuchasnonlinearity,non-Gaussiannoiseandunknownparameters.Aimingatnonlinearadaptivestatees-timationofmovingtargetinasystemwithunknownprocessnoiseandnon-Gaussianmeasurementnoise,anovelnoiseadaptivevariationalBayesian(VB)filterusingnaturalgradientisproposed.Firstly,aparameterizedinverse-Wishart(IW)distributionandastudent stdistributionareconstructedastheconjugatepriordistributionofpre-dictedstateerrorcovarianceandmeasurementlikelihoodrespectively.Then,intheframeworkofvariationalBayesianoptimization,thejointaposterioridistributionofestimationvariablesisapproximatelydecomposedintoindependentvariationaldistributionsbyusingmean-fieldtheory.Onthisbasis,thevariationaldistributionpara-metersofeachvariableareupdatedbycombiningcoordinateascendmethodandthecharacteristicsofexponentialdistributions.Furthermore,undertheconditionofmaximizingevidencelowerbound,thenaturalgradientswithre-specttostateestimationanditserrorcovariancearederivedbycombiningwithFisherinformationmatrix.Sothatthevariationaldistributionofnonlinearstategraduallyapproachestheposterioriprobabilitydensityfunction(PDF)ofstatealongthenaturalgradientdirection.Finally,simulationresultsshowthattheproposedalgorithmhasbetteradaptiveabilitytotheproblemofnoiseuncertaintyandcanobtainhigherestimationaccuracycom-paredtotraditionalalgorithms.Key wordsNonlinearfiltering,adaptivefiltering,variationalBayesianinference,naturalgradient,Fisherinforma-tionmatrixCitationHuYu-Mei,PanQuan,HuZhen-Tao,GuoZhen.AnovelnoiseadaptivevariationalBayesianfilterusingnaturalgradient.ActaAutomaticaSinica,2023,49(10):20942108收稿日期2021-10-13录用日期2022-03-01ManuscriptreceivedOctober13,2021;acceptedMarch1,2022国家自然科学基金(61790552,61976080),西北工业大学博士论文创新基金(CX201915)资助SupportedbyNationalNaturalScienceFoundationofChina(61790552,61976080)andInnovationFoundationforDoctorDis-sertationofNorthwesternPolytechnicalUniversity(CX201915)本文责任编委段海滨RecommendedbyAssociateEditorDUANHai-Bin1.西北工业大学自动化学院西安7100722.信息融合技术教育部重点实验室西安7100723.中国航空工业集团公司西安航空计算技术研究所西安7100654.河南大学人工智能学院郑州4500465.湖北航天技术研究院总体设计所武汉4300401.SchoolofAutomation,NorthwesternPolytechnicalUni-versity,Xi an7100722.KeyLaboratoryofInformationFu-sionTechnology,MinistryofEducation,Xi an7100723.Xi anAeronauticalComputingTechniqueResearchInstitute,AviationIndustryCorporationofChina,Ltd.,Xi an7100654.SchoolofArtificialIntelligence,HenanUniversity,Zhengzhou4500465.Sy-stemDesignInstitute,HubeiAerospaceTechnologyAcademy,Wuhan430040第49卷第10期自动化学报Vol.49,No.102023年10月ACTAAUTOMATICASINICAOctober,2023非线性滤波器设计与实现因其普适性及重要性一直是国内外学者研究的热点问题,近年来非线性状态估计理论已成功应用于陆、海、空、天中运动目标的预警与防御,智能交通的精确导航与制导,无人机定位与遥感监测、工业过程监控与故障诊断等众多领域13.考虑到运动目标跟踪系统机动、隐身等复杂多变的人为对抗特征以及非视距、干扰、遮挡等环境因素无法避免,其系统建模、估计与辨识过程中越来越无法回避非线性、非高斯以及参数未知等复杂系统特征的影响4.在对复杂环境下运动目标系统噪声先验信息进行建模时,建模误差存在于状态演化模型中,并通常假设其满足一定的参数化统计特性5.然而,在实际工程中,由于先验建模信息的不足导致难以对此类参数进行准确赋值.例如,在现代目标跟踪系统中出现的欺骗、干扰、杂波、未知分布特性的量测噪声和系统噪声等情形,尤其在非合作运动目标强机动场景中,因难以对其运动过程进行精细化建模,常造成目标跟踪航迹间断甚至无法正常起始的现象.针对系统噪声未知问题,其解决思路通常选取逆伽马分布和逆威沙特(Inverse-Wishart,IW)分布等具有统一参数表达形式的指数族分布作为共轭先验分布.Srkk等6在变分贝叶斯(VariationalBayes,VB)推断框架下利用指数族分布构建共轭先验分布近似未知量测噪声的后验概率密度函数(Probabilitydensityfunction,PDF),进而结合贝叶斯滤波机理实现时变噪声方差和目标状态的联合估计,其滤波效果达到与交互式多模型方法相接近的估计精度,并且凭借VB便于结合平均场近似解耦理论将高维变量求解转化为多个低维变量计算的特点,使其在解决多未知扰动问题方面更具有优势.随后,文献 7采用 IW 分布近似多变量噪声后验PDF 的思想,给出变分推断框架下噪声方差估计的一般实现形式.在此基础上,变分贝叶斯方法在未知量测/系统噪声估计方面得到了进一步的发展89,并成功推广应用于多目标跟踪环境1011.在文献 1213中滤波器的构建均建立在量测噪声和系统噪声统计特性未知的假设条件下,其中,文献 12在状态和量测扩维基础上通过采用批处理方式实现对未知噪声矩阵的估计,然而这种扩维方式不可避免造成计算复杂度的急剧增加;文献 13给出两类噪声统计特性均未知条件下噪声方差的在线估计方法,并且为避免系统噪声和状态相互独立假设的不合理性,采用独立于当前状态的前一时刻状态预测误差协方差代替系统噪声的统计特性,进而利用平均场理论近似解耦计算边缘后验PDF的变分分布,以迭代递推方式求解其估计变量期望的解析解,同时,将状态一步预测误差协方差的先验分布建模为参数化 IW 分布,以实现状态后验概率分布的更新.tttttt考虑到外界环境干扰以及非线性传播等因素的影响,量测数据概率分布往往呈现出重尾和非对称等非高斯特征.例如,在基于电磁信号的距离估计中,障碍物遮挡造成的非视距量测误差往往远大于其他来源服从对称分布的误差,导致量测分布呈现非对称非高斯现象.其次,受传感器精度和灵敏度的限制,当运动目标发生强机动时,雷达极化反射不稳定性将导致的量测随机缺失现象,也将造成量测产生非高斯特征.针对量测噪声非高斯问题的处理,文献 14 采用莱斯分布构建非共轭指数族变换点检测模型,并将其成功应用于雷达目标跟踪系统中.文献 15 则提出一种通过构建状态演化模型和量测模型的条件矩实现非线性非高斯系统的状态估计方法.为了进一步提升滤波器对不同分布形状噪声的鲁棒性,文献 16将Skew-分布作为非高斯量测似然的近似分布形式,在此基础上设计了一种鲁棒变分贝叶斯估计器.文献 17考虑在 散度最小化准则下,利用变分分布实现对后验PDF的近似.散度对异常数据具有较好的抑制作用,但是也打破了对数边缘概率与变分置信下界(Eviden-celowerbound,ELBO)、KL散度(Kullback-Lei-blerdivergence)之和的等式约束关系,其实质是一种伪VB推断方法.考虑到学生 分布和Skew-分布对因量测异常导致的量测重尾现象和非对称性具有较好鲁棒性1819,文献 20针对系统噪声和量测噪声均具有重尾特性的线性系统采用学生 分布对状态预测概率密度函数和量测似然函数进行建模,进而提出了一种基于高斯学生 混合分布的线性滤波器,与传统高斯假设条件下滤波器估计效果相比进一步提升了系统的状态估计的精度和鲁棒性.在系统噪声未知且时变和量测噪声重尾条件下,文献 21构建参数化 IW 分布作为状态一步预测误差协方差的共轭先验分布,同时,选取参数化学生 分布刻画具有厚尾特点的量测似然函数.然而,上述噪声自适应方法未考虑非线性滤波器自身的优化问题和如何衡量后验分布近似程度的问题.考虑到在非线性滤波器设计过程中,参数化变分分布对系统状态后验近似的程度是提高估计精度的关键因素之一,基于采样的随机优化和基于拟牛顿、高斯牛顿和梯度上升等确定性优化方法的非线性滤波器相继提出.在随机性优化方面,MCMC(MarkovchainMonteCarlo)2223和序贯蒙特卡罗(SequentialMonteCarlo,SMC)利用随机样本来逼近后验概率,并以样本分布近似未知变量后验分布.随机优化方法的优点是在大量样本的条件下能10期胡玉梅等:基于自然梯度的噪声自适应变分贝叶斯滤波算法2095够保证较高精度的估计,但要承担较大的计算负担24.在确定性优化方面,文献 25以最小二乘准则为目标函数,采用牛顿法推导给出一种迭代卡尔曼滤波器.文献 26综合梯度法和牛顿法提出一种无噪声条件下的滚动时域估计器,并证明其渐进收敛性和稳定性.文献 27针对非线性状态空间的状态估计及其估计误差协方差的迭代更新问题,利用正则化的非线性最小二乘思想提出一种随机增量近端梯度算法.VB方法通过求解参数化的目标函数,利用参数优化结果逼近后验分布,是一种确定性近似方法.因此,VB具有确定性优化方法特有的计算量较小的优势,同时由于VB便于结合平均场理论近似解耦联合概率分布,进一步减小了计算代价.非线性状态估计优化的实质是对多维状态后验PDF的近似逼近,并且其近似程度不能简单地采用欧氏距离进行度量.因此,选取合理度量准则将有利于提高后验PDF的近似程度,KL散度作为分布之间“差异”的度量在后验分布近似性衡量中具有天然优势.文献 17 给出状态估计的变分迭代优化实现形式,获得 散度和KL散度下对后验PDF更紧密的近似.同时,从信息几何角度出发,概率分布是统计流形上点,在一定条件下两概率分布之间的KL散度与作为统计流形度规的Fisher信息满足一定的数学关系.基于信息几何理论,Amari28利用自然梯度优化方式实现统计流形空间中目标函数的最速梯度下降(或上升).文献 29结合自然梯度策略和卡尔曼滤波框架设计一种非线性状态估计方法,在文献 30中,该方法进一步推广应用于传感器网络目标跟踪系统.文献 31中作者证明了针对非线性状态估计优化的自然梯度方法在克拉美罗下界意义下是渐进最优的.考虑自然梯度优化的优势,以变分置信下界最大化条件下状态估计及其估计误差协方差的自然梯度为切入点,从信息几何角度实现对状态后验概率密度函数的“紧密”逼近,进而提高状态估计精度.t在过程噪声先验不准确和由于量测随机丢失导致的量测噪声分布非高斯的情况下,针对非线性系统状态估计精度和鲁棒性提升问题,本文提出一种基于自然梯度的噪声自适应变分贝叶斯滤波算法.本文结构如下:第1节结合 IW 分布和学生 分布分别实现对状态预测误差协方差和量测似然的参数化建模;第2节结合平均场近似和坐标上升方法给出了变分贝叶斯框架下未知变量变分分布的迭代更新方式;第3节推导给出ELBO关于系统状态估计及其误差协方差的自然梯度,进而构建并设计一种基于自然梯度的噪声自适应非线性变分贝叶斯滤波器;第4节给出仿真验证与分析;第5节总结全文,并展望了后续研究方向.1预备知识kzkq(k|k)假设动态系统隐变量 的量测为,根据变分贝叶斯理论可知,在变分贝叶斯框架下以变分分布 作为桥梁可将难以积分的问题转化为ELBO的优化问题.logp(zk)=q(k|k)logp(k,zk)q(k|k)dkq(k|k)logp(k|zk)q(k|k)dk=L(k)+DKLq(k|k)p(k|zk)L(k)(1)L(k)DKLq(k|k)p(k|zk)q(k|k)p(k|zk)其中,表示具有单调递增特性的变分ELBO,表示变分分布 与后验分布 之间的KL散度,其表达式分别为L(k)=q(k|k)logp(k,zk)q(k|k)dk(2)DKLq(k|k)p(k|zk)=q(k|k)logp(k|zk)q(k|k)dk(3)q(k|k)k其中,表示以 为参数的变分分布.p(k|zk)p(k|zk)t未知变量的估计实际是对状态后验分布 的近似逼近,当非线性状态后验分布难以获取时,变分贝叶斯方法能够通过构建简单的变分分布实现对 的近似逼近.变分分布选取需考虑先验分布和后验分布具有相同的函数形式,一般情况下具有共轭特性的指数族分布满足这一条件.本文考虑系统噪声未知情况,根据标准卡尔曼滤波实现结构可知,系统噪声的统计特性仅影响状态预测误差协方差,因此可直接对状态预测误差协方差进行先验建模.并且当系统噪声假设服从均值已知但方差未知的多变量高斯分布假设时,可选取IW 分布作为其分布方差矩阵的共轭先验分布;此外,针对量测缺失导致量测噪声出现重尾问题,考虑选取能够有效表征这一现象的学生 分布对量测似然函数进行建模.1.1IW 分布共轭先验模型n nA在贝叶斯滤波框架下,对于系统噪声高斯分布参数未知问题的处理,可转化为状态预测误差协方差的估计问题.这里选取 IW 分布作为共轭先验分布,以保证后验分布与先验分布具有相同的函数形式.一个服从 IW 分布的 维随机对称正定矩阵 的概率密度函数可以表示为2096自动化学报49卷IW(A;t,T)=|T|t2|A|(t+n+1)22nt2n(t2)exp(0.5tr(TA1)(4)tTn ntr()n()nPk|k1其中,和 分别表示 IW 分布的自由度参数和 维逆尺度矩阵,表示矩阵的迹,表示 维的伽马函数.进而状态预测误差协方差 的先验分布表示为p(Pk|k1|zk1)=IW(Pk|k1;tk,Tk)(5)tk n+1EPk|k1对于逆威沙特分布,当 时,其变量期望即状态预测误差协方差的均值 表示为EPk|k1=Tktk n 1(6)tk令先验分布参数 表示为tk=+n+1(7)式(6)进一步转化为Tk=EPk|k1(8)其中,表示调谐参数7.根据式(5)(8)计算得到状态预测协方差的先验分布及其分布参数.t1.2学生 分布量测似然模型tk0Rkt当量测噪声分布由于量测值随机异常或缺失呈现重尾特征时,传统高斯噪声分布模型难以对噪声分布特性进行有效刻画.考虑学生 分布能够更好地体现这一特征,并且对异常值具有较好的鲁棒性.当量测噪声 满足均值为 、方差为 的分布参数时,建立参数化学生 分布模型,即p(k)=St(k;0,Rk,v)(9)vtRk其中,表示学生 分布的自由度参数,噪声方差 是该分布的尺度矩阵.在此条件下量测似然函数可表示为p(zk|xk)=St(zk;hk(xk),Rk,v)(10)tkt通常学生 分布的概率密度函数难以求得封闭解,为了解决这个问题,引入服从伽马分布的辅助随机变量,从而进一步将学生 分布转化为服从参数化高斯分布的伽马分布的积分,式(10)中量测似然函数可转化为p(zk|xk)=(N(zk;hk(xk),Rkk)G(k;k,k)dk(11)G(k;k,k)kkk其中,表示辅助变量 服从伽马分布,和 分别为形状参数和逆尺度参数20.综合上述特点,在传感器量测存在随机异常值情况下量测似然函数表示为如下结构化形式:p(zk|xk,k)=N(zk;hk(xk),Rkk)(12)p(k)=G(k;k,k)(13)G(k;k,k)Ek并且 及其均值 分别表示为G(k;k,k)=kkk(k)k1kekk(14)Ek=kk(15)kkk最后,通过对辅助变量 及其超参数 和 更新确定量测似然的表达式.1.3平均场近似kk=(xk,Pk|k1,k)根据以上分析建模,式(2)中的系统隐变量 定义为.考虑到联合后验分布形式的复杂性以及参数相互耦合的问题,无法直接求得其解析解.在假设各未知变量相互独立的前提下,通过引入平均场近似策略实现联合变分分布的近似解耦,即q(xk,Pk|k1,k)q(xk)q(Pk|k1)q(k)(16)q(xk)q(Pk|k1)q(k)xkPk|k1k其中,和 分别表示状态,状态一步预测协方差 和辅助变量 的变分分布,通过平均场理论对待估计变量联合分布进行分解,有助于在变分贝叶斯迭代框架结合坐标上升和各类优化方法实现多变量的联合优化.结合上述第 1.1 节和第 1.2 节中的共轭先验模型和量测似然模型,相应的变分分布分别表示为q(k)=G(k;k,k)(17)q(Pk|k1)=IW(Pk|k1;tk|k1,Tk|k1)(18)q(xk)=N(xk;xk|k,Pk|k)(19)kkktk|k1Tk|k1Pk|k1xk|kPk|kxk其中,和 为辅助随机变量 的超参数,同理,和 为随机未知变量 的超参数;和 为隐变量 的超参数.在变分贝叶斯优化过程中,可结合平均场理论将联合后验的变分分布近似解耦为多个单变量变分分布,并利用坐标上升方法分别对每个变量进行迭代更新.2未知变量的变分分布更新kxkPk|k1kkk=(xk,Pk|k1,k)根据上述构建的共轭先验模型和量测似然模型,变分置信下界中隐变量 包括状态,状态一步预测协方差 和辅助变量,可定义为.采用变分分布近似上述多未知变量的联合后验分布,即10期胡玉梅等:基于自然梯度的噪声自适应变分贝叶斯滤波算法2097p(xk,Pk|k1,k|zk)q(xk,Pk|k1,k)(20)进而,相应变分置信下界转化为L(k)=(q(xk,Pk|k1,k)logp(xk,Pk|k1,k,zk)q(xk,Pk|k1,k)dxkdPk|k1dk(21)xkPk|k1k在变分置信下界最大化约束下,结合坐标上升方法,变量,和 概率分布的对数形式的更新分别表示为logq(k)=Exk,Pk|k1logp(xk,Pk|k1,k,zk)+cxk,Pk|k1(22)logq(Pk|k1)=Exk,klogp(xk,Pk|k1,k,zk)+cxk,k(23)logq(xk)=EPk|k1,klogp(xk,Pk|k1,k,zk)+cPk|k1,k(24)cPk|k1,kcxk,kcxk,Pk|k1其中,和 分别表示在迭代过程中产生的实常数.结合建模信息和各变量变分分布的特点,联合后验分布及其对数形式分别表示为p(xk,Pk|k1,k|zk)=p(zk|xk,k)p(xk|Pk|k1,zk1)p(Pk|k1)p(k)N(zk;hk(xk),Rkk)N(xk;xk|k1,Pk|k1)IW(Pk|k1;tk|k1,Tk|k1)G(k;k,k)(25)logp(xk,Pk|k1,k|zk)=12k(zk hk(xk)TR1k(zk hk(xk)+(xk xk|k1)TP1k|k1()+(tk|k1+n+2)log|Pk|k1|+tr(Tk|k1P1k|k1)+2kk+(2 2k Dz)logk+ck(26)ckDznzkPk|k1(xkxk|k1)TP1k|k1()(xk xk|k1)TP1k|k1(xkxk|k1)tkPk|k1xk其中,表示迭代过程中与变量相关的实常数,和 分别表示向量 和矩阵 的维数,表示 的缩写形式(下文中类似情况不再一一说明).根据式(26)可知,由于高斯分布、学生 分布和伽马分布均属于指数分布族而具有简单的对数形式,因此便于采用迭代更新的方式计算其解析解.以下给出 、以及 的迭代估计的具体实现原理和步骤.k2.1 迭代更新i+1k根据式(22)和式(26),在第 次变分迭代更新过程中,的变分分布对数形式可表示为logq(i+1)(k)=12Eixk,Pk|k1k(zk hk(xk)TR1k()2kk+(2 2k Dz)logk+ck=12k(2k+tr(R1k)Eixk,Pk|k1(zk hk(xk)T()+(2 2k Dz)logk+ck(27)ckxkPk|k1xk=xik|k+xkxik|k xki其中,表示 和 积分后相关的实常数.根据贝叶斯无偏估计理论,其中 和 分别表示第 次估计值及其估计偏差.同时,结合局部线性化技术,式(27)中期望部分可进一步改写为Eixk,Pk|k1(zk hk(xk)T()(zk Hik(xik|k+xk)T()=(zk Hikxik|k)T()+(Hik)TPik|kHik(28)Hik=hk(xk)xk|xk=xik|kxk=xik|k其中,表示量测函数在 处的雅克比矩阵.i+1i+1ki+1ki+1k计算式(14)中伽马分布的对数表达式,并结合式(27)和式(28),不难看出,第 次更新后参数 的变分分布中超参数 和 的更新过程可表示为i+1k=ik+12Dz(29)i+1k=ik+12tr(R1kik)(30)ik=(zk Hikxik|k)T(zk Hikxik|k)+(Hik)TPik|kHik.其 中,ki+1i+1k根据式(15)中伽马分布的均值表达式,同时结合式(29)和式(30),参数 的第 次变分迭代更新 可表示为i+1k=i+1ki+1k=ik+12Dzik+12tr(R1kik)(31)kikikiki由式(31)可以清晰地看出,辅助随机变量 的更新与 负相关,并且在迭代过程中 的变化与上一次迭代的状态估计值及量测函数在此估计值处的局部线性化矩阵有关.根据 表达式,第 次2098自动化学报49卷i+1iki+1kiiki+1k(zk Hikxik|k)ik252/EE=/迭代的估计值(可理解为第 次迭代估计的先验)越接近真实状态,值越小,根据式(31)可知 则越大;相反,第 次的估计状态距真实状态较远,值越大,则越小.其物理意义是,当发生量测随机缺失时导致滤波更新的新息 增大,同时 值将随之增大,从而造成式(11)中实际量测噪声方差增大,这种情况将反馈于滤波过程中状态估计的自适应更新.为了直观地解释这种现象,结合式(11)和式(13),采用一维高斯分布和伽马分布说明其相关关系.图 1给出伽马分布参数对量测似然的影响示意图.在图 1中,假设高斯分布的均值和方差分别为 和,并且,伽马分布参数 和 的值分别在图 1(a)和图 1(b)中注明.当量测未发生丢失时,的值较小,高斯分布比较集中(如图 1(a)所示);当量测发生随机丢失,的值增大,高斯分布比较分散(如图 1(b)所示).01020304050 x(a)a=10,b=200.020.040.060.080.10概率密度伽马分布高斯分布伽马分布高斯分布01020304050 x(b)a=10,b=600.0050.0100.0150.0200.0250.0300.035概率密度图1伽马分布参数对量测似然函数的影响示意图Fig.1ThediagramoftheinfluenceofGammadistributionparametersonlikelihoodPk|k12.2 迭代更新i假设第 次迭代状态估计值已知,根据式(23)i+1Pk|k1和式(26),在第 次变分贝叶斯迭代更新过程中,的变分分布的对数形式表示为logq(i+1)(Pk|k1)=12Eixk,k(xk xk|k1)TP1k|k1()+(tk|k1+n+2)log|Pk|k1|+tr(Tk|k1P1k|k1)+cPk|k1=12tr(Eixk,k(xk xk|k1)()TP1k|k1)+(tk|k1+n+2)log|Pk|k1|+tr(Tk|k1P1k|k1)+cPk|k1(32)cPk|k1xkkEixk,k(xk xk|k1)(xkxk|k1)T其中,参数 表示 和 积分后的相关性实常数.式(32)中的期望 可转化为Eixk,k(xk xk|k1)()T=Pik|k+(xik|k xk|k1)()T=Pik|k+Kik(zk hk(xk|k1)()T(Kik)T=k(33)KikiitikTikPk|k1其中,表示在第 次变分贝叶斯迭代中卡尔曼滤波的增益.令第 次迭代中 IW 分布参数分别为 和,根据 IW 分布表达式(4),同时结合式(32)和式(33),则状态预测误差协方差 更新后变分分布的超参数分别为ti+1k=tik+1(34)Ti+1k=Tik+k(35)Pk|k1i+1Pi+1k|k1因此,根据式(6)中 均值表达式,同时结合式(34)和式(35),状态预测误差协方差的第 次变分迭代更新 表示为Pi+1k|k1=Ti+1kti+1k n 1=Tik+ktik n(36)zkhk(xk|k1)Pi+1k|k1根据卡尔曼滤波实现原理可知,当过程噪声方差建模不准确时,状态预测误差协方差较大,并且导致量测预测不精确而产生较大的量测新息.同时,由式(33)可知,的迭代更新与量测新息正相关,进而从物理意义上说明状态预测协方差建模的合理性.xk|k2.3 迭代更新i+1根据式(24)和式(26),在第 次迭代更新过程中,系统状态的变分分布对数形式表示为10期胡玉梅等:基于自然梯度的噪声自适应变分贝叶斯滤波算法2099logq(i+1)(xk)=12Ei+1Pk|k1,kk(zk hk(xk)TR1k()+Ei+1Pk|k1,k(xk xk|k1)TP1k|k1()+cxk=12Ei+1Pk|k1,kk(zk hk(xk)TR1k()+(xk xk|k1)TEi+1Pk|k1,kP1k|k1()+cxk(37)cxkPk|k1ki+1Pk|k1k其中,表示与 和 积分后相关性的实常数.在第 次迭代中,依据式(31)和式(36)已经获取状态预测误差协方差 和辅助变量 的估计值Ei+1Pk|k1,kk=i+1k(38)Ei+1Pk|k1,kP1k|k1=(Ei+1Pk|k1,kPk|k1)1=(Pi+1k|k1)1(39)由贝叶斯滤波理论可知,状态后验分布的近似变分分布应具有以下形式q(i+1)(xk)=1Ci+1kp(i+1)(zk|xk)p(i+1)(xk|zk1)=1Ci+1kN(zk;hk(xk),Rki+1k)N(xk;xk|k1,Pi+1k|k1)(40)Ci+1k其中,表示归一化常数Ci+1k=p(i+1)(zk|xk)p(i+1)(xk|zk1)dxk(41)结合高斯分布表达式,式(40)可进一步改写为logq(i+1)(xk)=1Ci+1kvuut|i+1k|(2)Dz+Dx|Rk|Pi+1k|k1|exp12(zkxkhk(xk)xk|k1)TRki+1k00Pi+1k|k11(zkxkhk(xk)xk|k1)(42)xi+1k|kPi+1k|k综合式(37)(42),目标状态后验分布的近似变分分布分别是以 和 为期望和协方差的高斯分布,即q(i+1)(xk)=N(xk;xi+1k|k,Pi+1k|k)(43)i+1Pi+1k|kxi+1k|kPi+1k|k针对非线性系统状态估计问题,现有贝叶斯滤波方法在实现式(43)中状态估计及其误差协方差的迭代更新过程中,往往难以获得与真实后验“紧密”近似的变分分布.虽然在多个变量坐标迭代优化过程中能够弥补部分由于线性化和参数不精确带来的误差,但估计精度提升效果受限,有时难以满足实际工程需求.因此,考虑在获得 和 估计值的基础上,在置信下界最大化条件下结合自然梯度方法和Fisher信息矩阵,推导给出 和 迭代更新解析表达式的具体实现.3基于自然梯度的噪声自适应变分贝叶斯滤波器3.1变分置信下界的自然梯度ikPk|k1ikPik|k1xk第 次变分迭代后,在分别获取参数 和 的估计值 和 的基础上,式(2)可转化为仅包含未知变量 的形式L(k)=q(xk|k)logp(xk,zk)q(xk|k)dxk(44)k=(xk|k,Pk|k)L(k)k定义此时变分参数,通过最大化式(44)中的置信下界即可获得状态估计及其误差协方差.自然梯度通过信息几何方法寻找统计流形空间的目标函数最速下降/上升方向,给出了一种实现式(44)置信下界求解的可行性方案25.下面以置信下界 最大化为目标函数,同时结合Fisher信息矩阵推导关于变分参数 的自然梯度.式(44)的置信下界可以进一步表示为L(k)=Eqlogp(zk|xk)DKLq(xk|k)p(xk)(45)ii+1p(xk)p(xk)q(xk|ik),k假设第 次迭代是第 次迭代的先验,即式(45)中的先验信息 可以近似为 则变分参数 的最优估计值表示为k=argmaxkEqlogp(zk|xk)DKLq(xk|k)q(xk|ik)(46)kk通过计算式(46)等号右边关于 线性化函数,可获取置信下界的自然梯度以及相应的最优.Eqlogp(zk|xk)k=ikDKLq(xk|k)q(xk|ik)k=iki+1定理 1.假设系统状态演化满足高斯分布特性,并且对数似然期望 在 的邻域内一阶可导,同时,在 的邻域内二阶可导,则第 次迭代过2100自动化学报49卷程中变分置信下界的自然梯度为ei+1k=F1ikkEqlogp(zk|xik|k)(47)kk=k ik 0证明.根据变分参数的最优表达式(46)可知,通过对其求关于 的偏导数即可获得置信下界最大化条件下的极值点.由于对数似然期望一阶可导,令,根据泰勒展开上式对数似然期望线性化后可以近似表示为Eqlogp(zk|xk)=Eqlogp(zk|xk)k=ik+kEqlogp(zk|xk)k=ikk+O(k2)(48)kEqlogp(zk|xk)k=ikikO(k2)k2其中,表示对数似然期望在 处的雅克比矩阵,表示 的高阶无穷小.DKLq(xk|k)q(xk|ik)k=ik同时,由于 二阶可导,式(46)中第 2 项在 点处进行泰勒展开后可表示为DKLq(xk|k)q(xk|ik)=12TkFkk+O(k3)(49)Fk其中,为具有半正定性Fisher信息矩阵,即Fk 2kDKLq(xk|k)q(xk|ik)(50)此时,式(46)进一步转化为k=arg maxk0kEqlogp(zk|xk)k12TkFkk(51)ki+1k计算式(51)关于 的偏导,并令其取值为零,可获得在第 次变分迭代中 ELBO关于 的自然梯度ei+1k=F1ikkEqlogp(zk|xik|k)(52)k因此,在置信下界最大化条件下变分参数 的更新表达式为i+1k=ik+ei+1k=ik+F1ikkEqlogp(zk|xk)(53)需要说明的是,VB迭代更新等价于在统计流形空间沿以Fisher信息矩阵作为权重的自然梯度方向移动,从而获得后验分布的最佳近似,并且具有以下性质:FikkEqlogp(zk|xik|k)1)自然梯度推导中,和 q(xk|k)的计算依赖于参数化的变分分布.但是KL散度作为变分分布选择的度量是不变的.k 02)时,式(46)中的第 2 项 KL散度表达式转化为式(49),因此从局部意义上迭代过程中的变分分布与后验分布的近似误差是正定二次型.k3)自然梯度是 Fisher 有效的25,因此上述优化过程是Fisher有效的估计器.特别地,上述 的迭代估计是无偏的.k 0FikkEqlogp(zk|xik|k)q(xk|k)4)在自然梯度的推导