混合滤波器及其应用

自动化学报

ACTAAUTOMATICASINICA

Vol.39,No.1January,2013

基于高斯和近似的扩展切片高斯混合滤波器及其在

多径估计中的应用

陈杰1

程兰2

甘明刚1

摘要全球卫星导航系统(Globalnavigationsatellitesystem,GNSS)信号的多径估计问题实际上是条件线性状态空间模型下的状态估计问题.根据高斯和理论提出了适用于非高斯噪声环境的扩展切片高斯混合滤波(ExtensionofslicedGaussianmixturefilter,ESGMF)算法.该算法将非高斯噪声的状态概率密度函数(Probabilitydensityfunction,PDF)表示为高斯和的形式,将ESGMF通过一组并行的切片高斯混合滤波器(SlicedGaussianmixturefilter,SGMF)来实现.同时,在ESGMF算法中利用粒子滤波(Particlefilter,PF)中重采样的思想对成指数增加的状态预测PDF的高斯混合个体进行约简,以提高贝叶斯推理的效率.该算法可以获得非高斯噪声下状态PDF的迭代解析表达式.最后,将ESGMF应用于GPS多径参数估计,仿真结果表明,ESGMF算法的估计精度优于基于PF和扩展卡尔曼滤波(ExtendedKalmanfilter,EKF)的算法.关键词

非高斯噪声,高斯和,概率密度函数,切片高斯混合滤波器,多径估计

引用格式陈杰,程兰,甘明刚.基于高斯和近似的扩展切片高斯混合滤波器及其在多径估计中的应用.自动化学报,2013,39(1):1−10DOI

10.3724/SP.J.1004.2013.00001

ExtensionofSGMFUsingGaussianSumApproximationforNonlinear/Non-GaussianModelandItsApplicationin

MultipathEstimation

CHENJie1

CHENGLan2

GANMing-Gang1

AbstractThemultipathestimationofglobalnavigationsatellitesystem(GNSS)signalisactuallythestateestimationofnonlinear/non-Gaussiansystems.TheextensionofslicedGaussianmixturefilter(ESGMF)basedonGaussiansumapproximationisproposedforthestateestimationofnonlinear/non-Gaussianstatespace,andtheprobabilitydensityfunction(PDF)expressionofstatesisderivedrecursivelyforatimevaryingsystem.ResamplingisappliedtothepredictionPDFtoreducethecomplexityofBayesianinference.ThesimulationresultofmultipathestimationwithESGMFshowsthattheESGMFalgorithmperformsbetterinaccuracythanthealgorithmsbasedonparticlefilter(PF)andextendedKalmanfilter(EKF).

KeywordsNon-Gaussiannoise,Gaussiansum,probabilitydensityfunction(PDF),slicedGaussianmixturefilter(SGMF),multipathestimation

CitationJieChen,LanCheng,Ming-GangGan.ExtensionofSGMFusingGaussiansumapproximationfornonlinear/non-Gaussianmodelanditsapplicationinmultipathestimation.ActaAutomaticaSinica,2013,39(1):1−10

状态估计[1]在很多领域具有广泛的应用,如监测、督管和定位[2].对于时变系统而言,迭代的估计方法是非常有用的[3−4],而卡尔曼滤波器是线性高斯噪声系统的最优迭代估计器[5].对于非线性高斯

录用日期2012-03-30

ManuscriptreceivedAugust12,2011;acceptedMarch30,2012国家自然科学基金(61120106010),国家杰出青年科学基金资助(60925011)资助

SupportedbyNationalNaturalScienceFoundationofChina(61120106010)andNationalScienceFundforDistinguishedYoungScholars(60925011)本文责任编委胡小平收稿日期2011-08-12

RecommendedbyAssociateEditorHUXiao-Ping

1.北京理工大学自动化学院复杂系统智能控制与决策教育部重点实验室北京1000812.太原理工大学信息工程学院太原030024

1.KeyLaboratoryofComplexSystemIntelligentControlandDecision(MinistryofEducation),BeijingInstituteofTechnol-ogy,Beijing1000812.CollegeofInformationEngineering,TaiyuanUniversityofTechnology,Taiyuan030024

该文的英文版同时发表在ActaAutomaticaSinica,vol.39,no.1,pp.1−10,2013.

噪声系统的状态估计,研究者提出了扩展卡尔曼滤波器(ExtendedKalmanfilter,EKF).对于非线性非高斯系统,基于EKF的高斯和滤波器(Gaussiansumfilter,GSF)可实现状态估计[6−8].在GSF中,利用测量序列,根据高斯和近似理论以及贝叶斯准则来更新状态的概率密度函数(Probabilitydensityfunction,PDF)[9].然而,当系统的非线性程度较高时,GSF的性能会显著降低[10].

对于非线性非高斯系统,粒子滤波(Particlefilter,PF)是另外一种状态估计方法[11].然而,与GSF一样,PF不能通过确定方式得到状态PDF的迭代解析表达式.而且PF存在样本退化和枯竭的问题.为了克服这些问题,研究者提出了改进PF算法[7,11−12].然而,PF算法的计算复杂度随着状态维数的增加而显著增加[13].对包含有线性子结构的

2自动化非线性系统,边缘PF可有效地降低算法的复杂度.PF算法的共同缺点是不能得到状态PDF的迭代解析表达式[14].

为了在具有线性子结构的非线性系统中得到状态的迭代解析PDF表达式,Klumpp等提出了切片高斯混合滤波(SlicedGaussianmixturefilter,SGMF)算法[5],该算法将非线性子状态的PDF用混合的Dirac函数表示,将线性子状态的PDF用混合的高斯函数表示.SGMF算法在文献[15]中用于协同的无源目标跟踪,但该应用的前提是系统中只存在高斯噪声.

多径干扰是高精度定位系统(如差分全球卫星导航系统(Globalnavigationsatellitesystem,GNSS))的主要误差源[16].多径估计延迟锁定环(Multipathestimationdelaylockloop,MEDLL)利用极大似然估计理论估计多径以抑制其影响.MEDLL实现了从传统延迟锁定环到数字延迟锁定环的转变,是实现数字延迟锁定环的第一次尝试[17].此后,形式各异的多径估计方法被用于多径抑制,如极大似然估计器、EKF估计器[18].然而,这些方法多数只适用于高斯噪声环境,在非高斯噪声环境下其估计性能会显著降低.以PF算法为原型的方法被用于非高斯噪声环境下的多径估计[14,18−19],但时变系统的维数较高时,这类算法的计算复杂度会显著增加.

根据高斯和近似理论,本文提出了适用于非高斯噪声环境的扩展切片高斯混合滤波(ExtensionofslicedGaussianmixturefilter,ESGMF)算法.该算法将噪声PDF近似为高斯函数加权和的形式后,即可将SGMF算法扩展到非高斯噪声环境下.同时,将PF的重采样思想应用于高斯混合个体的约简,目的是提高贝叶斯推理的有效性.ESGMF算法解决了非高斯噪声下的多径估计问题,且其估计精度高于已有的算法.

本文结构组织如下:第1节描述GNSS系统的多径估计问题;第2节介绍了SGMF算法;第3节提出了ESGMF算法,阐述了ESGMF的主要思想,并推导了状态PDF的迭代解析表达式;第4节讨论了ESGMF算法的实现过程;第5节将ESGMF算法用于多径估计;第6节得出结论.

1问题描述

图1是GNSS的信号跟踪结构.与本地载波

混频后,接收信号r(k)与本地产生C/A码的向量c(k−τˆ0−di)=[c(k−τˆ0−d1),···,c(k−τˆ0−

dN)]T

进行相关运算,其中,di(i=1,···,N)是第i路本地码与即时码的相关器间距,N是相关器数目;相关器输出向量yk用于估计多径参数向量Ak和τk;经过后续处理后,如多径重构、多径消除,得到期望信号的时间延迟估计τˆ0,用于控制本地伪码发生器,调整本地产生的即时码c(k−τˆ0),使之与接收信号r(k)的伪码同步.

学报39卷

在多径环境下接收到的GNSS信号可以表示为(M+1)个路径的模型,即一路直接信号与M路多径信号之和,附加噪声信号n(t)[20].不考虑低通滤波器的影响[21−23],图1中即时支路的相关器输出表示为

yk(Ak,0,Ak,i,εk,τk,i)=Ak,0R(ε󰀌Mk)+

Ak,iR(εk+τk,i)+nk(1)

i=1

其中,Ak,0和Ak,i是k时刻直接信号与第i路

多径信号的复合幅度;R(ε)是理想自相关函数;εk=τˆk,0−τk,0,τˆk,0是τk,0的本地估计;τk,0和τk,i分别是k时刻直接信号时间延迟与第i路多径信号的相对时间延迟;nk为k时刻的噪声.通常,nk被认为是高斯噪声,但有时并非如此,如在室内和室外环境下的无线通信中经常遇到的脉冲信号,这种信号的PDF可以近似为有限个高斯个体的混合形式[24−26].

图1信号跟踪结构

Fig.1

Structureofsignaltracking

从式(1)可以看出,相关器输出向量yk=[yk×(Ak,0,Ak,i,ε,τk,i,d1),yk(Ak,0,Ak,i,ε,τk,i,d2),···,yk(Ak,0,Ak,i,ε,τk,i,dN)]T关于幅度向量Ak=[Ak,0,Ak,1,···,Ak,M]T是条件线性的,即给定时间延迟向量τk=[εk,0,τk,1,···,τk,M]T时,yk(·)关于Ak是线性的.所以Ak和τk的估计问题可以转化为包含线性子结构的非线性系统模型的状态估计问题.因此,本文重点考虑混合线性/非线性系统[5]:󰀉

xlk+1=Fkl

(xnk)xlk+Blk(xnk)uk+wlk

xnk+1=Fkn(xn(2)k)+wn

k

yk=Hk(xnln

k)xk+hk(xk)+vk

(3)

其中,xllk∈R和xnk∈Rr

分别表示线性状态子向

量和非线性子状态向量;Fkl(·)和Fkn

(·)分别是线性子空间和非线性子空间的系统矩阵,二者都依赖于

非线性子状态向量xnk;Bl

k(·)表示状态空间模型的非线性部分;uk是输入;wl独立的系统噪声,其PDF分别为k和wnk

分别表示相互p(wlH(·),k)和p(wn

k);k是测量矩阵也依赖于非线性子状态向量xnh(·)k,k是测量模型中的非线性部分;vk是测量噪声,PDF为p(vk).本文的目的是得到k时刻状态xk的PDFp(xk)的准确估计.对于具有状态空间模型的系统而言,估计问题可以通过预测和滤波来完成.

在k+1时刻,预测PDFpP(xk+1)可以通过下

1期陈杰等:基于高斯和近似的扩展切片高斯混合滤波器及其在多径估计中的应用3

式来计算:

󰀍

pP(xk+1)=

pt(xk+1|xk)pe(xk)dxk

(4)

其中,pe(xk)是xk的估计PDF,pt(xk+1|xk)是依赖于状态空间模型的状态转移PDF.

估计PDFpe(xk)可以通过式(5)来计算:

pe(xk)=c·pL(yk|xk)pP(xk)

(5)

其中,pL(yk|xk)是似然函数,似然函数可以认为是对󰀋于状态xk出

现测量值yk的条件概率.c=1/pL(yk|xk)pP(xk)dxk是常数.

除非是在线性高斯噪声情况下,通常不能获得式(4)和式(5)的解析表达式.所以,在非线性非高斯情况下得到PDFp(xk)需要必要的近似处理.

2SGMF

2.1切片高斯混合表示

由于状态空间模型包含线性子结构,为了更加有效地计算式(4)和式(5),可以将全空间的估计问题分解为线性子空间和非线性子空间的估计问题.

切片高斯混合PDFp(xk)将非线性子空间上的状态PDF表示为Dirac混合的形式,将线性子空间上的状态PDF表示为高斯混合的形式[5].

p󰀒xlk

,xn󰀁󰀊k=S󰀃󰀄ωp󰀒x󰀂󰀁k,sδxnk−ξnlk󰀂ξ(6)k,sk,s

s=1

其中,δ(·)是Dirac函数,ωk,s是Dirac函数的非负

权值,S为切片的数目,ξn

k,s∈Rr是每个切片的位置.通常,由式(6)表示的切片高斯混合PDF可以近似估计器中的任意PDF[5].

将非线性子状态xnk的边缘密度表示为Dirac函数加权和的形式,即

p(xn󰀌S󰀃󰀄Sk)

=

ωk,sδxn

n󰀌

k−ξk,s,ωk,s=1

(7)

s=1

s=1

对于给定切片ξn

k,s,xlk的条件概率密度表示为高斯混合的形式:

p

󰀒

xl󰀂󰀂󰀁

󰀌Tk

󰀃kξk,s=

υk,stNxlll

󰀄k−µk,st,Ck,st

(8)

t=1

其中,󰀊Tk

l

l

l

t=1υk,st=1,υk,st,µk,st∈R和Ck,st∈Rl×l分别表示第s个切片的第t个个体的条件权值、均值和方差,Tk是k时刻高斯混合的个体数目.N(µ,C)是均值为µ、方差为C的高斯分布.需要强调的是高斯个体的所有参数都是在给定切片位置

ξnk,s

的条件下计算的[15].2.2SGMF原理

因为估计器中的任意PDF都可以用式(6)表示,所以xk的先验PDF和预测PDF都可以用这种形式来表示.因此,SGMF可以通过迭代的方式获得状态后验PDF的解析表达式,主要处理步骤如图2所示.

图2SGMF处理步骤

Fig.2ProcessingstepsofSGMF

在SGMF的处理过程中,系统噪声ωk和过程噪声νk+1都是高斯噪声,所以SGMF并不适用于非高斯噪声环境.

3ESGMF

3.1ESGMF原理

ESGMF将SGMF扩展到非高斯噪声下,主要处理步骤如图3所示.

图3ESGMF处理步骤

Fig.3ProcessingstepsofESGMF

pnG(wk)和pnG(vk+1)表示p(wk)和p(vk+1)为非高斯概率密度,即ESGMF中系统噪声和测量噪声都为非高斯噪声.为了得到非高斯噪声下状态后验PDF的解析表达式,需要利用高斯和理论.高斯和理论认为任意PDF都可以表示为高斯PDF的加权和,即引理1.

引理1[6].n维向量x的密度函数f(x)可以近

4自动化学报39卷

似为任意精度的加权高斯和,即时刻每个切片的个体数目.

f(x)=

L󰀌i=1

aiN(µi,Bi)

(9)

3.2.2高斯混合的近似

3.2状态概率密度表达式的推导可以看出,预测之后状态的PDF由切片高斯

混合的形式变为高斯混合形式,且高斯个体的数根据引理1,系统噪声和测量噪声可以表示为

pωeeω

目由STk变为P0STk.令ωk式(10)和式(11),其中,λk,p为系统噪声中高斯个+1,g=λk,pωk,sυk,st,

lnng=(p−1)S+(s−1)Tk+t,则式(16)可以体的权值;(mlk,p,Qk,p)和(mk,p,Qk,p)分别对应系

表示为式(17),其中,Nk+1=P0STk.为了能够统高斯个体中线性部分和非线性部分的均值和方差.

υυ

λk,q、mk,q和Rk,q为测量噪声中高斯个体的权值、使状态估计通过迭代的方式实现,需要将高斯混

合形式的预测PDF重新近似为切片高斯混合形均值和方差.P0和Q0分别是系统噪声和测量噪声

pl

式,即将式(17)近似为式(18),其中,υkPDF的高斯个体数目.+1,sg=

ppnpnpn

ωk+1,gN(ξk+1,s−ζk+1,g,Qk+1,g)plpl根据状态转移模型(2)和测量模型(3),可以将󰀊G,µ=µppnpnpnk+1,sgk+1,g,g=1ωk+1,gN(ξk+1,s−ζk+1,g,Qk+1,g)转移密度函数pt(xk+1|xk)和似然函数pL(yk|xk)

pnplpn

Qpl表示为式(12)和式(13).k+1,sg=Qk+1,g,ωk+1,s和ξk+1,s可利用文献

假设初始状态PDF为式(14),其中,S为切片[27]提出的次最优方法来计算.计算出所有切片数目,T0为初始时刻每个切片的高斯个体数目.k的权值和位置之后,对切片权值进行归一化,即󰀊Spnpnpn时刻的状态估计PDF表示为式(15),其中,Tk是kωk+1,s=ωk+1,s/s=1ωk+1,s.

󰀅󰀇󰀈󰀇󰀈󰀆P0P0

lll󰀌󰀌w−mQ0kk,pk,p

,,λω(10)λωp(wk)=k,p=1k,pNnnn

w−m0Qkk,pk,pp=1p=1

p(vk)=

P0

󰀌p=1L

Q0󰀌q=1

T0S󰀌󰀃󰀄󰀃󰀄󰀌󰀒󰀁enelelnenelnlel˜˜0,t,Cpx0,x0=ω0,sυ0,stδx0−ξ0,sNx0−µ0,st

s=1t=1

TkS󰀌󰀄󰀃󰀄󰀃󰀌󰀒󰀁enelelnenelnlel˜˜k,t,Cpxk,xk=ωk,sυk,stδxk−ξk,sNxk−µk,st

s=1t=1

υQ0󰀌q=1

󰀊L

对于整数L,正标量值αi,有i=1ai=1.

根据引理1,系统噪声和测量噪声的非高斯PDF都可以近似为高斯PDF加权和的形式,故而ESGMF可以由一组并行的SGMF来实现,其中预测、近似和滤波过程分别在第3.2.1节到第3.2.3节介绍.

3.2.1预测

用p˜(x)表示x的真实PDF,用p(x)表示x的近似PDF.根据式(4)和多维正态分布函数以及Dirac函数的性质,可以得到k+1时

pn

刻的预测PDF如式(16)所示,其中,ζk+1,pst=󰀒en󰀁pnplnnFkξk,s+mnk,p,Qk+1,pst=Qk,p,µk+1,pst=󰀒en󰀁󰀒en󰀁el

plll

˜k,st+Bkξk,sµξk,suk+mlFk

k,p,Qk+1,pst=󰀒󰀁el󰀒l󰀒en󰀁󰀁T

˜Flξen+Qlk,sCk,stFkξk,sk,p.

k

λυk,qN−

Q0

󰀒󰀁󰀌υ

vk−mk,q,Rk,q,λυk,q=1

q=1

(11)

Qlk,p0

󰀈󰀆0

Qnk,p

(12)

󰀅󰀇Nλω

k,p

pt(xk+1|xk)=

xlk+1

llln

Fk(xnk)xk−Bk(xk)uk

nnn

xnk+1−Fk(xk)−mk,p

−mlk,p

󰀈󰀇,

p(yk|xk)=

󰀄󰀃

nlnυ

λk,qNyk−Hk(xk)xk−hk(xk)−mk,q,Rk,q

(13)

(14)

(15)

󰀍

p˜P(xk+1)=

pt(xk+1|xk)·pe(xk)dxk=

󰀈󰀇󰀈󰀆󰀅󰀇TkP0󰀌S󰀌plpll󰀌Qk+1,pst0xk+1−µk+1,pstenel

λωωυN,pnpnk,sk,stnk,p

0Qx−ζk+1k+1,pstk+1,pstp=1s=1t=1

(16)

1期陈杰等:基于高斯和近似的扩展切片高斯混合滤波器及其在多径估计中的应用

elυk+1,st

5

重新近似过程可以看作是非线性子状态的更

新过程,在这个过程中线性子状态保持不变.所以重新近似过程中,近似算法的性能直接决定非线性子状态的估计精度.同时,预测后每个切片具

pl

有相同的高斯个体数,这是因为µplk+1,sg=µk+1,g,

pl

Qplk+1,sg=Qk+1,g,这表明预测后的高斯个体数目决定了近似后每个切片的高斯个体数目.

=

󰀊Nk+1󰀊Q0

g=1

plυ

γsgqυk+1,sgλk+1,q

q=1

plυγsgqυk+1,sgλk+1,q

(t=(g−1)Q0

3.2.3滤波

k+1时刻的状态估计PDFpe(xk+1)可以通过

滤波步来计算.根据式(5)可以得到式(19).对式(19)中每个切片的高斯个体权值归一化可得:

Nk+1

+q)值得注意的是,归一化是在每个切片下分别进行的,而不是对所有的个体权值进行归一化.

从式(19)的推导过程可以看出,估计过程实际上是在给定的非线性子状态下对线性子状态的更新过程.每个切片有Nk+1Q0个相同的高斯个体,对应Nk+1Q0个不同的权值.所以式(19)可以重新表示为式(20),其中Tk+1=Nk+1Q0.

至此,根据图2完成了k+1时刻状态后验PDF的迭代解析表达式的推导,而k+1时刻状态的估计值由式(21)计算.−µplk+1,g

pn−ζk+1,g

󰀈󰀇,Qplk+1,g

0

󰀈󰀆0Qpnk+1,g

(17)

p˜P(xk+1)=

󰀌

g=1

󰀅󰀇p

ωkN+1,g

xlk+1

xnk+1

pP(xk+1)=

SNk+1

󰀌󰀌s=1g=1

󰀒n󰀁󰀒l󰀁pnpnplplpl

ωkυδx−ξNx−µ,Qk+1k+1,sk+1+1,sk+1,sgk+1,sgk+1,sg

(18)

pe(xk+1)=cpL(xk+1)pP(xk+1)=

c

Q0

󰀌q=1SNk+1󰀌󰀌s=1g=1

pnpnplplplnl

ωk+1,sυk+1,sgδ(xk+1−ξk+1,s)N(xk+1−µk+1,g,Qk+1,sg)=SNk+1

󰀌󰀌s=1g=1

nlnυ

λυk+1,qN(yk+1−Hk+1(xk+1)xk+1−hk+1(xk+1)−mk+1,q,Rk+1,q)×

c

Q0

󰀌q=1

λυk+1,qpnpnpln

ωkυδ(x−ξk+1k+1,s)×+1,sk+1,sg

plplllnυ

N(yk+1−Hk+1(xnk+1)xk+1−hk+1(xk+1)−mk+1,q,Rk+1,q)N(xk+1−µk+1,sg,Qk+1,sg)=󰀐󰀏󰀎󰀑ConditionallylineardimensionS󰀌s=1

Nk+1Q0

pnn

ωk+1,sδ(xk+1

c

−ξpnk+1,s)

󰀌󰀌

g=1q=1

plυl˜el˜elγstqυkk+1,sgq,Ck+1,sgq)=+1,sgλk+1,qN(xk+1−µ

c

S󰀌s=1

Nk+1Q0

enennωk+1,sδ(xk+1−ξk+1,s)

󰀌󰀌

g=1q=1

plυl˜el˜elγsgqυkk+1,sgq,Ck+1,sgq)+1,sgλk+1,qN(xk+1−µ

(19)

其中

pnen

ωk+1,s=ωk+1,s,

ξpnξenk+1,s=󰀅󰀅k+1,s,

󰀆󰀆󰀒pn󰀁pl

󰀒󰀒pn󰀁󰀁T󰀒pn󰀁plyk+1−Hk+1ξk+1,sµk+1,sg

󰀒pn󰀁HQγsgq=N,R+Hξk+1ξk+1,sk+1,qk+1k+1,sgk+1,s

−hk+1ξk+1,s−mυ

k+1,q󰀃󰀒󰀒pn󰀁󰀁T󰀄󰀁󰀁󰀒󰀒󰀒pn󰀁plTpnpl

Rk+1,q+Hk+1ξk+1,sQk+1,sgHk+1ξk+1,sKsgq=Qk+1,sgHk+1ξk+1,s

1

c=Q0SNk+1󰀊󰀊󰀊plenυ

ωk+1,sγsgqυk+1,sgλk+1,q

s=1g=1q=1

6自动化学报39卷

从推导过程可以看出,非线性子状态和线性子状态可以分别进行估计;非线性子状态的估计精确度由将高斯混合近似为切片高斯混合的近似算法决定;线性子状态的估计精度由非线性子状态的估计精度和卡尔曼滤波的精度共同决定,这是因为卡尔曼滤波是在给定非线性子状态下进行的.

4ESGMF的实现

4.1基于重采样的高斯混合约简

从式(17)可以看出,在k+1时刻,预测PDF的高斯个体数目为Nk+1=P0STk.预测PDF被重新近似为切片高斯混合的形式之后,有S个切片,每个切片有Nk+1个高斯个体.同时,由于非高斯测量噪声也被建模为具有Q0个个体的高斯和的形式,这样经过滤波之后,后验PDF具有S个切片,每个切片有Tk+1=Nk+1Q0个个体.所以,在k+2时刻预测后高斯混合个体数目增加到Nk+2=P0STk+1.根据递推的方法可知,k时刻预测后的高斯个体数

(k−1)

目为Nk=(P0S)kQ0T0.然而,在实际应用中对成指数增加的高斯个体进行处理是很困难的,这也增加了贝叶斯推理的复杂度.所以,应该限制高斯个体的数目.

从式(17)可以看出,个体的重要性是通过其权值来体现的.对于给定的预测PDF,只需要少数的个体就可以近似整个预测PDF.通常,认为权值大的个体在近似过程中产生较大的作用.这可以通过式(18)来解释,在式(18)中权值越大的的个体的

pn

切片位置ξk+1,s与高斯个体中非线性部分的均值ζpnk+1,g之间的似然度越大,其高斯个体对应的权值也越大.所以,在实际应用中,可以通过G个个体的加权和来近似预测PDF,这样就剔除了权值较小的个体.k+1时刻,经过约简后的预测PDF表示为式(22),其中pPr(xk+1)是约简后的预测PDF.显然,对于给定的误差代价函数,G依赖于预测PDF

S󰀌s=1S󰀌s=1G󰀌g=1

G󰀌g=1

p

的形状和权值ωk+1,g.假设在k+1时刻,预测PDF可以通过G个高斯个体的加权和来近似.在实际应用中,在每一时刻都保持个体数目恒定是很有利的.本文利用粒子滤波中重采样的思想保持个体数目稳定.重采样可以剔除权值较小的个体,保留权值大的个体.算法1给出了一种简单的重采样方法.

算法1.高斯混合约简子算法

p

对Nk+1个体进行降序排列.1)根据权值ωk

+1,g

2)保留前G个个体,将G个权值和个体分别pr

记为ωk+1,g和

󰀈󰀆󰀈󰀇󰀅󰀇

prlprll

0Qk+1,gxk+1−µk+1,g

,Nprnprn

0Qkxnk+1−ζk+1,g+1,g

g=1,···,G.

3)定义权值阈值ωthreshold

pr

当(min(ωk+1,g<ωthreshold))forg=1到G

prpr

根据G个个体的权值ωk+1,1,···,ωk+1,G,

prnprlrl

从1,···,G中选取j.令{µpk+1,g,ζk+1,g,Qk+1,g,

prnprlprlprnprrn

Qpk+1,g}={µk+1,j,ζk+1,j,Qk+1,j,Qk+1,j},ωk+1,g=1/G.

endforend.

算法1中,权值阈值ωthreshold依赖于期望的精度,即期望精度越高ωthreshold越小.经过大量的仿真结果发现当ωthreshold<0.05时,可以得到期望的精度.经过高斯混合约简后,可以得到式(23).

重采样可以平衡个体数目G和这些个体权值之间的关系,使计算量不会浪费在产生不重要的个体上.而且,通过复制权值大的个体,可以降低贝叶斯推理的复杂度.经过高斯混合约简,Nk个高斯个体由G个个体来近似,这就保证了在每一时刻预测PDF的高斯个体数目保持不变.

Tk+1

p(xk+1)=

e

ennωk+1,sδ(xk+1−ξenk+1,s)

󰀌

t=1

ell˜el˜elυk+1,stN(xk+1−µk+1,st,Ck+1,st)

(20)

xnk+1

=

enen

ωk+1,sξk+1,s,xlk+1xlk+1

xnk+1

=

S󰀌s=1

Tk+1

en

ωk+1,s

󰀌

t=1

el

˜elυk+1,stµk+1,st

(21)

󰀈󰀆

󰀅󰀇

pr

ωk+1,gN

pPr(xk+1)=

prl

−µk+1,g

prn−ζk+1,g

󰀈󰀇,

prlQk+1,g

0

prn

Qk+1,g

00

rn

Qpk+1,g

(22)

󰀅󰀇

pr

ωk+1,gN

pPr(xk+1)=

xlk+1

xnk+1

rl

−µpk+1,g

prn−ζk+1,g

󰀈󰀇,

rl

Qpk+1,g

󰀈󰀆,

G󰀌g=1

0

pr

ωk+1,g=1

(23)

1期陈杰等:基于高斯和近似的扩展切片高斯混合滤波器及其在多径估计中的应用7

4.2ESGMF算法

ESGMF算法的实现过程如算法2所示.在ESGMF算法中,初始化设置已知的参数.迭代通过递归的方式来估计状态,每次迭代都进行预测、高斯混合约简、切片高斯混合近似、更新以及状态估计.高斯混合约简之后,对权值进行归一化,以保证约简后的函数仍然具有PDF的性质.基于同样的原因,对预测PDF用切片高斯混合的形式近似之后,也对每个切片的高斯个体权值分别进行归一化.

算法2.ESGMF算法//初始化

S,Ten0,ωen0,s,υel0,st,ξ0,s,µ˜el0,st,C˜0el,st,P0,λω0,p,ml0,p,

mn0,p,Ql0,p,Qn0,p,Q0,λυ

υ0,q,m0,q,R0,q.

//迭代

fork=1toN//1)预测

对于p=1,···,P0,s=1,···,S1,···,Tp,tpn

=

0,根据式(18)计算ωk,pst,ζpnk,pst,Qk,pst,

µ˜plk,pst,Q˜plk,pst.

//2)高斯混合约简定义高斯个体数目G.

根据高斯混合约简算法,只保留G个个体,记

其权值、均值和方差为ωprprnrn

µprlprlk,g,ζk,g,Qpk,g,k,g,Qk,g,g=1,···,G,并对G个个体的权值进行归一化.

//3)切片高斯混合近似

根据文献[24]中Dirac混合近似的次最优算法计算ωpplplpl

k,s,ξpnk,s,υk,sg,µk,sg,Ck,sg,s=1,···,S,g=1,···,G.并对每个切片的高斯个体权值分别进行归一化.

//4)更新

对于s=1,···,S,g=1,···,G,q=1,···,Q0,根据式(19)计算状态后验PDF的参

数ωenenk,s,ξk,s,υelk,sg,µ

˜el˜el//5)k,sg,Ck,sg.状态估计

xn

󰀊Sk=ωenenk,sξk,s,xl=󰀊SGQs=1

kωen󰀊0s=1

k,sυelt=1

k,stµ˜elk,st.//6)为下次迭代做准备

T0=GQ0,ωen

0,s=ωenenenk,s,ξ0,s=ξk,s,υel0,st=υelk,st,µ˜el0,st=µ˜elk,st,C˜0el,st=C˜k,stel.end

5多径估计

根据式(1),令xT

T

k=[(Ak),(τk)T],xk∈

R2(M+1)×1.假设多径参数在观测时间内保持不变,则系统模型可以表示为一阶马尔科夫模型

xk+1=Fkxk+wk(24)

wk表示非高斯系统噪声.令xlk=Ak,xn

k=τk,则

模型(24)可以表示为如下形式:

󰀇xl󰀈󰀇󰀈󰀇󰀈

Fklxlkwlk+1k

xn=F(25)

k+1knxn

+kwn

k

其中,Fkl∈RM+1和Fn

∈RM+1分别为线性状态和非线性状态的转移矩阵k

,wl

和wn是相互独立的.为了估计多径参数,需要k

N个相关器输出值k,相关器输出向量可以表示为

y󰀒xn󰀁l

k=Rk

xk+vk(26)其中,yk=[yk(Ak,0,Ak,i,ε,τk,i,d1),yk(Ak,0,Ak,i,ε,

τk,i,d2),···,yk(Ak,0,Ak,i,ε,τk,i,dN)]T.vk为非高斯测量噪声.

从式(25)和式(26)可以看出,如果H(xnR(xnk)=

k),Blk(xnk)uk=0,hk(xn

k)=0,非高斯噪声下的多径估计问题可以转化为非高斯噪声下包含线性子结构的非线性状态空间模型的状态估计问题.

不失一般性,仿真中采用GPS的C/A码信号,假设只有直接信号和一路多径信号,即M=1.考虑最坏的情形,即直接信号与多径信号同相[26],且此后的仿真只考虑这种情形.设定归一化的幅度A0=1,A1=0.6.设定时间延迟τ0=10Tc,τ1=0.5Tc,Ts=Tc/10,其中,Tc是C/A码的码元持续时间,Ts是采样间隔.τ0的本地估计是τˆ0=10.2Tc,所以ε=τˆ0−τ0=0.2Tc.本文的目的是估计τ0、A1和τ1,这里τ0=τˆ0−ε.

在系统模型(24)中,考虑Fk=I4×4,即认为多径参数在观测周期内保持不变,wk和vk是可以由高斯混合模型表示的非高斯噪声.系󰀅󰀇统噪声的PDF󰀈󰀇设定为p(ω0)=󰀊l󰀈󰀆

Pλ0−ml0ωω0,pQl0p=10,pNωn,0,p

,0−mn0,p0Qn

0,p

其中,P0=2,λω0,1=0.6,λω0,2=0.4,ml

0,1=

[0,0]T,mlT0,2=[0.2,0.2],mnT0,1=[0,0],mn

0,2=

[0.2,0.2]T

,Ql0,p=Qn0,p=diag{0.01,0.01}.

󰀊测量噪声的PDF可以表示为p(νk)=Q0υ

k−mυk,q,Rk,q),其中,Q0=2,λυ

q=1λk,qN(νk,1=

0.7,λυk,1=0.3,mυk,1=[0,0,0,0,0]T,mυ

k,2=[0.1,0.1,0.1,0.1,0.1]T,Rk,1=Rk,2=diag{0.01,0.01,0.01,0.01,0.01}.

为了说明ESGMF算法的可行性,在初始状态的估计PDF中考虑了两个高斯个体:第一个高斯个体的均值是真值,第二个高斯个体的均值在一定程度上偏离真值.两个个体的方差是根据经验设定的.仿真中初始状态的估计pe󰀒󰀃PDFxl0,xn󰀁󰀊G0󰀃设为

󰀄0=ω0jNxl0−µ˜elj=10,j,C˜0el,j

×󰀃Nxn

en0−µ

˜0,j,C˜󰀄󰀄

0en,j这里,G0=2,ω0,1=0.6,ω0,1=0.4,

8自动化µ˜el0,1=[10.6]T,µ˜en0,1=[0.20.5]T

,µ˜el0,2=

[0.90.5]T,µ˜en0,2=[0.150.45]T,C˜0el,1

=C˜0en,2=diag{0.01,0.01}.初始时刻状态PDF的切片高斯混合形式可以由第3.2.3节的方法来计算.

将预测PDF近似为切片高斯混合的形式时要计算次边缘积分,这需要考虑状态的初始积分域[24].仿真中将状态的初始积分域设为

󰀇󰀇󰀈󰀇󰀈󰀈󰀇󰀇󰀈󰀇󰀈󰀈

xl∈00,1−0.50.5

1,xn∈0,

1

这是基于如下考虑:归一化的直接信号幅度和

多径信号幅度小于“1”.直接信号的时间延迟估计偏差小于0.5Tc,多径信号相对于直接信号的时间延迟小于Tc.

将文献[7]中提出的III型高斯和滤波器(Gaus-siansumfilterIII,GSF3)和III型高斯和粒子滤波器(GaussiansumparticlefilterIII,GSPF3)以及文献[8]中提出的GSF作为对比算法.状态的均方根误差(Rootmeansquareerror,RMSE)作为状态估计的性能指标.直接信号时间延迟和多径相关输出的估计结果如图4和图5所示,其中多径相关输出等于A1R(ε−τ1),ESGMF算法中预测PDF的高斯个体数目限制为G=10(滤波后GSF、GSF3和GSPF3的高斯个体数目也限制为G).从图4可以看出,ESGMF算法具有最好的估计性能.而且随着切片数目S的增加,ESGMF算法的估计性能显著提高,而增加GSPF3和GSF3算法的粒子数目Np二者性能改善并不明显,如图5所示.

从图4和图5可以看出,与GSPF3、GSF3以及GSF相比,ESGMF具有较好的估计性能.这是因为GSPF3是基于PF的,而PF是根据重要性密度函数产生的随机样本进行状态估计的,通常我们并不能获得最优的重要性密度函数;GSF3是基于EKF估计和PF预测的,GSF是只基于EKF的,EKF需要对状态方程进行线性化,线性化过程带来了误差;而ESGMF是基于确定性样本选择策略的,其后验PDF可以通过式(20)获得迭代的解析表达式.

不难看出RMSE越大定位误差越大,因此定义如下性能指标:

In=M1×293m+M2×7.5m

GPSC/A码的一个码元持续时间对应的空间距离为Tc×c=977.5ns×(3×108m/s)≈293m,GPS载波一个周期对应的空间距离为c/fc=(3×108m/s)/(10.23MHz)≈30m.而最大载波跟踪误差为四分之一周期[28],所以载波跟踪误差引起的最大定位偏差为30m/4=7.5m.其中,M1为直接信号时间延迟的RMSE,M1为多径相关输出的RMSE.在性能指标In下不同算法的引起的定位偏差如表1所示,可以看出在该指标下得出的结论与

学报39卷

从图4和图5中得出的结论是吻合的.

图4S=20,Np=50时的估计结果

Fig.4EstimationresultswhenS=20,Np=50

图5S=50,Np=100时的估计结果

Fig.5EstimationresultswhenS=50,Np=100

1期陈杰等:基于高斯和近似的扩展切片高斯混合滤波器及其在多径估计中的应用

表1不同算法的定位偏差(m)

Positioningerrorofdifferentalgorithms(m)

GSF345.07543.5783

GSPF324.1922.4192

ESGMF21.33507.6250

9

Table1

算法

GSF43.572541.4329

S=20,Np=50S=50,Np=100

虽然ESGMF性能比其他三种算法好,但仍存在估计偏差,这主要是由以下原因造成的:

1)用于确定切片位置的算法为了降低计算复杂度,本文采用了可以快速计算切片位置的次最优算法,这种算法的复杂度与状态的维数成线性关系.同时切片的数目也影响状态估计的性能,这是因为非线性子状态的PDF是由一系列的Dirac函数来近似的.

2)状态的初始积分范围为了降低数值积分的计算量,本文选择了状态PDF的主要部分对应的状态范围,忽略了边缘区域,这样根据次边缘积分计算的切片位置可能偏离根据真实的边缘积分获得的切片位置.

3)高斯混合约简算法本文只保留了权值较大的G个个体,并通过重采样对G个个体中权值较小的个体进行剔除,并用权值大的个体进行替换以保持个体数目不变.这虽然可以有效地提高贝叶斯推理的效率,但也一定程度上造成了误差.而且正如前文所说,G的大小也影响状态估计的精度.

6结论

本文根据多径下GNSS信号的特点,将非高斯噪声环境下的多径估计问题转化为具有线性子结构的非线性系统的状态估计问题.并在非高斯噪声环境下,针对混合线性/非线性系统提出了ESGMF算法,扩展了SGMF算法.利用高斯和理论将ES-GMF通过一组并行的SGMF来实现.将估计器中的状态PDF表示为切片高斯混合的加权和形式,可以在非高斯噪声下得到状态后验PDF的闭环迭代解析表达式.而且ESGMF算法将重采样思想用于高斯混合约简,提高了贝叶斯推理的复杂度.最后,将ESGMF算法用于多径估计,仿真结果表明所提出算法的多径估计精度优于对比算法.

References

1AlfiA.PSOwithadaptivemutationandinertiaweightanditsapplicationinparameterestimationofdynamicsystem.ActaAutomaticaSinica,2011,37(5):541−5492HanebeckUD,BriechleK,RauhA.ProgressiveBayes:anewframeworkfornonlinearstateestimation.In:Proceed-ingsoftheInternationalSocietyforOpticalEngineering.Orlando,FL,USA:SPIE,2003,5099:256−2673AnDe-Xi,LiangYan,ZhouDong-Hua.Arobustadaptivefilterbasedononlinefilterparameteridentification.ActaAutomaticaSinica,2004,30(4):560−566

(安德玺,梁彦,周东华.一种基于滤波参数在线辨识的鲁棒自适应滤波器.自动化学报,2004,30(4):560−566)

4SongQi,HanJian-Da.AnadaptiveUKFalgorithmforthestateandparameterestimationsofamobilerobot.ActaAu-tomaticaSinica,2008,34(1):72−79

(宋崎,韩建达.一种用于移动机器人状态和参数估计的自适应UKF算法.自动化学报,2008,34(1):72−79)5KlumppV,SawoF,HanebeckUD,FrankenD.TheslicedGaussianmixturefilterforefficientnonlinearestimation.In:Proceedingsofthe11thInternationalConferenceonInfor-mationFusion.Cologne,Germany:IEEE,2008.1−86SorensonHW,AlspachDL.RecursiveBayesianestimationusingGaussiansums.Automatica,1971,7(4):465−4797KotechaJH,DjuricPM.Gaussiansumparticlefilter-ing.IEEETransactionsonSignalProcessing,2003,51(10):2602−26128LinQing,YinJian-Jun,ZhangJian-Qiu,HuBo.Gaus-siansumfilteringmethodsfornonlinearnon-Gaussianmod-els.SystemsEngineeringandElectronics,2010,32(12):2494−2499

(林青,尹建军,张建秋,胡波.非线性非高斯模型的高斯和滤波算法.系统工程与电子技术,2010,32(12):2494−2499)9CaputiMJ,MooseRL.AmodifiedGaussiansumapproachtoestimationofnon-Gaussiansignals.IEEETransactionsonAerospaceandElectronicSystems,1993,29(2):446−45110AlspachDL,SorensonHW.NonlinearBayesianestimation

usingGaussiansumapproximations.IEEETransactionsonAutomaticControl,1972,17(4):439−44811ZhuZhi-Yu.ParticleFilterAlgorithmandItsApplication.

Beijing:SciencePress,2010.27−32

(朱志宇.粒子滤波算法及其应用.北京:科学出版社,2010.27−32)12SahaN,RoyD.AGirsanovparticlefilterinnonlinear

engineeringdynamics.PhysicsLettersA,2009,373(6):627−63513SchrempfOC,HanebeckUD.Recursivepredictionof

stochasticnonlinearsystemsbasedonoptimaldiracmixtureapproximations.In:Proceedingsofthe2007AmericanCon-trolConference.NewYork,USA:IEEE,2007.1768−177414SchonT,GustafssonF,NordlundPJ.Marginalizedparticle

filtersformixedlinear/nonlinearstate-spacemodels.IEEETransactionsonSignalProcessing,2005,53(7):2279−228915HorstJ,SawoF,KlumppV,HanebackUD,FrankenD.

ExtensionoftheslicedGaussianmixturefilterwithapplica-tiontocooperativepassivetargettracking.In:Proceedingsofthe12thInternationalConferenceonInformationFusion.Seattle,WA,USA:IEEE,2009.587−594

10自动化16ChenA,ChaboryA,EscherAC,MacabiauC.Compar-isonsofmultipathmodelingstrategiesfortheestimationofGPSpositioningerror.In:Proceedingsofthe3rdEuropeanConferenceonAntennasandPropagation.Berlin,Germany:IEEE,2009.2300−230417VanNeeDJR.Themultipathestimationdelaylockloop.

In:Proceedingsofthe2ndIEEEInternationalSymposiumonSpreadSpectrumTechniquesandApplications.Yoko-hama,Japan:IEEE,1992.39−4218LentmaierM,KrachB,RobertsonP.Bayesiantimedelay

estimationofGNSSsignalsinmultipathdynamicenviron-ments.InternationalJournalofNavigationandObservation,2008.ArticleID:37265119WangXian,YuanHong.GPSmultipathestimationbased

particlefiltering.ControlandDecision,2010,25(8):1139−1148

(王宪,袁洪.基于粒子滤波的GPS多径估计.控制与决策,2010,25(8):1139−1148)20SahmoudiM,AminMG.RobusttrackingofweakGPSsig-nalsinmultipathandjammingenvironments.SignalPro-cessing,2009,89(7):1320−133321ChengLan,ChenJie,GanMing-Gang.Multipatherror

analysisofcarriertrackingloopinGPSreceiver.In:Pro-ceedingsofthe29thChineseControlConference.Beijing,China:IEEE,2010.4137−4141

(程兰,陈杰,甘明刚.GPS接收机载波跟踪多径误差分析.第29届中国控制会议论文集.北京,中国:IEEE,2010.4137−4141)22ChengLan,ChenJie,GanMing-Gang.Modelingandanal-ysisofmultipatheliminatingtechniqueofPNcoderang-ingreceiver.TransactionsofBeijingInstituteofTechnology,2011,31(3):318−323

(程兰,陈杰,甘明刚.伪码测距接收机中MET技术误差建模与分析.北京理工大学学报,2011,31(3):318−323)23ChenJ,ChengL,GanMG.ModelingofGPScodeand

carriertrackingerrorinmultipath.ChineseJournalofElec-tronics,2012,21(1):78−8424SenguptaD,KayS.Efficientestimationofparametersfor

non-Gaussianautogressiveprocesses.IEEETransactionsonAcoustics,SpeechandSignalProcessing,1989,37(6):785−79425SteinDWJ.DetectionofrandomsignalsinGaussianmix-turenoise.IEEETransactionsonInformationTheory,1995,41(6):1788−180126SoubielleJ,FijalkowI,DuvautP,BibautA.GPSposition-inginamultipathenvironment.IEEETransactionsonSig-nalProcessing,2002,50(1):141−15027KlumppV,HanebeckUD.Diracmixturetreesforfast

suboptimalmulti-dimensionaldensityapproximation.In:Proceedingsofthe2008IEEEInternationalConferenceonMultisensorFusionandIntegrationforIntelligentSystems.Seoul,Korea:IEEE,2008.593−600

学报39卷

28SunLi,WangYin-Feng,ZhangQi-Shan.Multipatheffects

onGPScarrierphaseobservable.ActaAeronauticaetAs-tronauticaSinica,1998,19(S1):76−78

(孙礼,王银锋,张其善.多径效应对GPS载波相位观测量的影响.航空学报,1998,19(S1):76−78)

陈杰北京理工大学自动化学院教授.分别于1986年、1993年、2000年获得北京理工大学控制理论与控制工程专业学士、硕士、博士学位.1989年至1990年为美国加利福利亚州立大学的的访问学者.1996年至1997年为英国伯明翰大学E&E学院的研究员.主要研究方向为多指标优化,多目标决策与控制,受限

非线性控制.E-mail:chenjie@bit.edu.cn

(CHENJieProfessorattheSchoolofAutomation,Bei-jingInstituteofTechnology.Hereceivedhisbachelor,mas-terandPh.D.degreesincontroltheoryandcontrolengi-neeringfromBeijingInstituteofTechnologyin1986,1993,and2000,respectively.From1989to1990,hewasavisit-ingscholarinCaliforniaStateUniversity,USA.From1996to1997,hewasaresearchfellowinSchoolofE&E,theUniversityofBirmingham,UK.Hisresearchinterestcov-ersmulti-indexoptimization,multi-objectivedecisionandcontrol,intelligentcontrol,andconstrainednonlinearcon-trol.)

程兰太原理工大学信息工程学院讲师.分别于2005年和2008年获得太原理工大学士和硕士学位,2012年获得北京理工大学博士学位.主要研究方向为GPS和北斗接收机的伪码和载波跟踪以及多径缓和.本文通信作者.E-mail:taolan1983@126.com

(CHENGLanLectureratTaiyuan

UniversityofTechnology.ShereceivedherbachelorandmasterdegreesfromTaiyuanUniversityofTechnologyin2005and2008,respectively.ShereceivedherPh.D.degreeinpatternrecognitionandintelligentsystemsfromBeijingInstituteofTechnologyin2012.Herresearchinterestcov-ersthePNcodetrackingandcarriertrackingofGPSandBeidoureceiver,andmultipathmitigation.Correspondingauthorofthispaper.)

甘明刚北京理工大学自动化学院讲师.分别于2002年和2007年获得北京理工大学控制理论与控制工程专业学士和博士学位.主要研究方向为智能信息处理和图像处理.E-mail:agan@bit.edu.cn(GANMing-GangLectureratBei-jingInstituteofTechnology.Here-ceivedhisbachelorandmasterdegrees

incontroltheoryandcontrolengineeringfromBeijingIn-stituteofTechnologyin2002and2007,respectively.Hisresearchinterestcoversintelligentinformationprocessingandimageprocessing.)