Air and Missile Defense College,Air Force Engineering University,Xi'an 710051,China

Immune adaptive Gaussian mixture particle fi lter for state estimation

Wenlong Huang*,Xiaodan Wang,YiWang,and Guohong Li

Air and Missile Defense College,Air Force Engineering University,Xi'an 710051,China

The particle fi lter(PF)is a fl exible and powerfulsequential Monte Carlo(SMC)technique capable of modeling nonlinear, non-Gaussian,and nonstationary dynamical systems.However, the generic PF suffers from particle degeneracy and sample impoverishment,which greatly affects its performance for nonlinear, non-Gaussian tracking problems.To deal with those issues,an improved PF is proposed.The algorithm consists of a PF that uses an immune adaptive Gaussian mixture model(IAGM)based immune algorithm to re-approximate the posterior density.At the same time,three immune antibody operators are embed in the new fi lter.Instead of using a resample strategy,the newest observation and conditionallikelihood are integrated into those immune antibody operators to update the particles,which can further improve the diversity of particles,and drive particles toward their close localmaximum of the posterior probability.The improved PF algorithm can produce a closed-form expression for the posterior state distribution.Simulation results show the proposed algorithm can maintain the effectiveness and diversity of particles and avoid sample impoverishment,and its performance is superiorto several PFs and Kalman fi lters.

arti fi cial immune,particle fi lter,Gaussian mixture model.

1.Introduction

The particle fi lter(PF)is an attractive estimation procedure for non-linear dynamical systems[1].It is proven to be very successful for solving non-linear and non-Gaussian state estimation problems,including target tracking,speech recognition,fi nancialeconometrics,and mobile robot[2-6].In the generic PF,the algorithms approximate the posterior densities of the hidden states with a setofparticles,and the particles are sampled from the state space and the posterior is updated by propagating the particles and updating their weights based upon the observations'likelihoods.Generally,the performance of the PF is superior to other non-linear fi lters[7].In principle,an infi nite number of particles can exactly represent any given probability density function(PDF).However,the requirementfor the number of particles grows exponentially with dimension,which makes PF sufferdue to“curse of dimensionality”.This problem is exacerbated by the high dimensionality of the state space[8-10].

Moreover,the sample degeneracy phenomenon willunfortunately unavoidablly occur in PF,when there are only a few particles in the vicinity of the true state.Only these particles have signi fi cantweights,and the weights of most of the other particles are near zero.The posterior density cannot be approximated with small sample sets properly, and signi fi cantcomputation is wasted on particles with low weights.In order to reduce the degeneracy effect,several related algorithms have been introduced that incorporate density estimation methods into a particle fi lter[10-14]. In[15],Cham and Rehg employed a piecewise Gaussian modelto representthe measurementfunction and the prediction density.In[16],Vermaak et al.directly modeled the posteriordensity as a Gaussian mixture model(GMM). A kernel-based Bayesian fi lter was introduced in[17].An analytical form for the measurement function p(zk|x1∶k) is estimated from the data using multistage sampling and density interpolation methods.After multiplying the two functions,a smooth approximation for the posterior density is found using a variable-width mean-shift technique. While these methods exploit an analytical form to maintain multiple hypothesesorsimplify the sampling problem, PFs sufferfrom a high computationalcost,especially when large numbers of particles are required to modelthe posteriordistribution.

On the other hand,during the resampling step,particles having large weights are repeatedly selected and theparticles having negligible weights are discarded,which makes multiple samples fall at the same point in the state space[18].As a result,those identicalsamples cannotwell remain the true representative of the posterior distribution.This phenomenon is called the particle impoverishment[19-21].In order to reduce the impoverishmenteffect,several complex resampling strategies are employed or speci fi c prior knowledge aboutthe objects are provided [22-26],in which,the most famous ones are partitioned sampling[25]and unbiased sampling[26].Although those methods can alleviate the impoverishment effect to a certain extent,all of them are not effective in reducing the variance of sample weights which increases overtime.

In order to reduce the impoverishment effect,we propose an immune adaptive GMM(IAGM)based arti fi cial immune algorithm to re-approximate the posteriordensity. The IAGM is obtained over the weighted sample set and each sample is re-weighted via this mixture model.This model is subsequently updated using the new observation to produce a closed-form expression for the posterior state distribution.This approach can provide an analytical expression for the posterior distribution by IAGM that can be used to identify its modes.Instead of employing complex discrete resampling strategies,we can directly sample from the posterior.

Moreover,an immune clone and mutation operator is incorporated into the new fi lter,induced by the newestobservation information and the likelihood function.The immune clone and mutation iteration can drive particles toward the local modes of the posterior probability effectively.At the same time,this procedure can perform implicitly iterative sampling through immune antibody diversity operators,thus we can obtain more diversity of particles sampled to alleviate the particle impoverishmentphenomenon greatly.Then the posterior density of the state is re-approximated by IAGMwith the immune diversity samples to maintain the effectiveness and diversity of samples.

2.IAGMdensity estimate

2.1 Particle impoverishment problem

During the resampling period of PF,the particle impoverishment problem always arises,because when we try to approximate a continuous distribution with a discrete one, there is a strictly positive probability that two samples drawn from the distribution will be identical.On the contrary,with a continuous distribution,no two samples drawn from the distribution will be identical.Thus we approximate the fi ltering distribution with a continuous GMM density estimate function in an analytic form.The advantage of the analytic representation is that it provides a globalview of the landscape of the likelihood function and thus enables ef fi cient sample placement.The GMM obtained by an adaptive global optimal algorithm can minimize a certain distance measure between the true fi ltering distribution and the GMM density estimate function.

2.2 Immune adaptive Gaussian mixture model

According to[27],we can directly modelthe posteriordensity as a GMM.However these models are optimized by the expectation maximization(EM)algorithm and adapta fi xed number of components instead of fl exible number of componentsin each subnet.As we know,the EMalgorithm only produces solutions that are locally optimal and may not achieve the globally optimal solution,thus it is sensitive to initialization.Moreover,the numberofcomponents of the mixture modelmustbe known in advance.This paper introduces a novelglobalsearch mechanism based on the arti fi cialimmune clonalselection algorithm into GMM techniques.The novelalgorithm can improve the effectiveness in estimating the parameters and determining simultaneously the optimalnumber of components automatically in each mixture modelsubnetwhich allows us to approximate non-Gaussian multimodaldensities.

The immune system(IS)can be considered to be a remarkably ef fi cient and powerful information processing system which operates in a highly paralleland distributed manner[28-30].Clonal selection is an adaptive search technique designed to fi nd outnear optimalsolutions to the large-scale optimization problem with multiple localmaxima.In this paper,we integrate the process of EMinto the clonalselection algorithm,thus obtain a method which is able to simultaneously performestimation ofthe modelparameters and determine the optimalnumberofcomponents in each subnet GMMautomatically.

2.2.1 Basis conception

In IAGM,antigen represents a problem,and antibodies represent candidates of the problem.Each subnet GMM is coded as an antibody to represent a possible solution, and each antibody gene is composed of two parts.The fi rst gene segmentis binary,which is used to encode the number of components.Each of these bits is related to a particular component.The second gene segment consists of fl oating point values representing parameters of the models.The parametersforeach componentinclude the mixing proportionπk,the mean vectorμkand the covariance matrixΣk.

The minimum description length(MDL)criterion is used as the antibody-antigen af fi nity for model selection. The best individual is the one that has the lowest MDL value.The MDL criterion is a consistent estimator ofmodelorder thatis expressed as

The fi tness of the antibody

where MDLmaxand MDLminare the highest and the lowest MDL values respectively.The fi tness of an antibody is evaluated by invoking an R-iteration EM algorithm.Starting from the mixture modelindicated by an antibody,the algorithm runs the E-step and the M-step for R iterations.In our implementations R is setto 3.

2.2.2 IAGMdensity estimate algorithm

The IAGMprocess can be visualized as follows:Step 1Randomly generate n antibodies:

Step 2Perform R EM steps on each antibodyCompute the antibody-antigen af fi nities MDL ofallthe antibodies.Antibody clone:

Step 3Mutation

The fi rstmutation operator inverses the binary value of each antibody gene in the fi rstpartwith the mutation prob-which modi fi es the numberof cluster K.tion probability which varies with the antibody fi tness,and prandis a uniform random variable over[0,1].

The second mutation operator modi fi es the parameters of the mixture model and it only works on the antibody genes.The number of the genes willoffer mutation for the antibody?ai(t),which is shown as

We adoptthe Gaussian mutation to generate a new antibody as follows

where g(0,σ)is a Gaussian random variable,ωis a variable,andω=0.1.

Step 4Af fi nity maturation:Perform R EM steps on antibody population Z(t).Compute the antibodies fi tness of Z(t).Step 5will die according to the death probability pd.We setξ=1e?4.

Step 7t=t+1:If the stop criterion is achieved,then stop;otherwise go to Step 2.

In IAGM-PF,IAGM density estimation can form continuous estimates of the posterior as follows:

where p(xk|θj)is the conditionalprobability density function corresponding to the j th component.

3.Immune antibody operator

In the implementation of general particle fi ltering,the resampling step is crucialbecause withoutit,the variance of the particle weights quickly increases,i.e.very few normalized weights are substantial.However,the resampling strategy results in the particle impoverishment problem, i.e.the loss of diversity for the particles.In order to deal with particle impoverishment,three immune antibody diversity operators are embedded in the new method.The iterative procedure of immune clone can redistribute particles to the localmodesofthe observation effectively.Atthe same time,this procedure can derive the particles towards high likelihood regions to relieve the particle impoverishmentphenomenon greatly.

3.1 Antibody generating operator

We obtain the antibody population Ab(k)at the k th step from the continuousestimate ofthe posteriorapproximated by IAGM.

where SC(k)(·)is a sampling operator,C(k)is decided bynormalized particle weights as follows:

where nc＞n is a given value related to the clone scale, Int(·)is a supremum function,and N is the number of particles.

3.2 Antibody mutating operator

In order to incorporate the mostcurrentmeasurements,we mutate every generated antibody by sigma points scaled unscented transformation.

where i=1,2,...,C(k),and Kkcan be obtained by sigma points scaled unscented transform[13].

Calculate sigma points:

Propagate the particle with time update:

Incorporate new observation,i.e.measurementupdate: where nxis the dimension of the state x,λ=α2(nx+κ),+(1?α2+β),

The scaled unscented transformation parameters are set toα=1,β=0 andκ=2.

3.3 Antibody selection operator

Compute the conditional likelihood1,2,...,C(k))for every antibody,then select Abi?(k)

4.Improved particle filter

Consider the evolution of the state sequence{xk,k∈N} given by

where f is a possibly nonlinearfunction of the state xk?1, {vk?1,k∈N}is an i.i.d.process noise sequence,and N is the step numbers.

The measurementofthe state sequence is

where h is a possibly nonlinear function,{nk?1,k∈N} is an i.i.d.measurement noise sequence.In particular,we seek fi ltered estimates of xkbased on the set of all available measurements y1∶kup to time k.

In this paper,IAGMdensity estimation is used to form continuous estimates of the posterior.This model is subsequently updated using the immune antibody operators which operate the newestobservation to produce a closedform expression for the posterior state distribution.The corresponding algorithm is described below.

Step 1Initializationfrom the initialization state with

the density p(x0);

(ii)Form the initial GMMwith the EMalgorithm.

Step 2Generating the proposaldistribution

The posterior state density is approximated by the IAGM as(6).Then the IAGM is used as the proposaldistribution in the following measurementupdate step.

Step 3Importance samplingate the unnormalized particle weights:

Step 4Immune antibody operators,obtain the new particle and its corresponding weights with(12).

Step 5Estimate the state with

Step 6Fitthe particles using the immune adaptive EM algorithm and attain the particle representation based on the Gaussian mixture.

Step 7Set k=k+1 and go back to Step 2.

5.Experimentalresults

In order to evaluate the performance of the proposed immune adaptive Gaussian mixture particle fi lter algorithm (IAGM PF),a number of simulations are performed.The simulation results are compared with other fi lter algorithms in terms of state estimation accuracy,including extended Kalman fi lter(EKF),unscented Kalman fi lter (UKF),generic particle fi lter(Generic PF),Monte Carlo particle fi lter(MCMC PF)[13],kernel-based Bayesian fi lter(KBBF)[12],partitioned sampling particle fi ler(PS PF) [13]and deterministic resampling unbiased sampling particle fi ler(DR-PF)[14].Allof the experiments are carried outon a CPUIntelCore(TM)3.4 GHz PCwith 4 GB memory.

5.1 Univariate nonstationary growth model

The univariate modelis a non-stationary,highly non-linear modelthatis also bimodalin nature.Itis given by the following setof equations:

where uk?1and vkare zero mean Gaussian random variables with variancesσvandσu,respectively,i.e.uk?1∈In this experiment,the following values for the parameters are considered:α=0.5, β=25,γ=8,σv=1,the totalnumber of time steps is setas 50.

We consider 50 particles for all of the four algorithms. Fig.1 shows the true trajectory of the target along with the fi ltering estimates of one exemplar run,and Fig.2 shows observations corresponding estimate error.In order to quantify the performance,the simulation is run 100 times with re-initialization to generate statisticalaverages. We use the traditionalmeasure ofthe performance:the root mean squared error(RMSE)and variance to evaluate those fi lters.The performance of all the different fi lters is summarized in Table 1,including the following aspects:the means,variances of the RMSE of the state estimates and the average executeion time.From these results,we can conclude thatthe average accuracy of ouralgorithm is better than EKF,Generic PF and MCMC PF.

Table 1 Experimental results of state estimation of UNGMwith 50 particles

Fig.1 Single run result oftrue trajectory and filtering solutions for 50 particles

Fig.2 Distribution of corresponding estimate error for Fig.1

For the sake of completeness,we set different parameters to simulate the four algorithms,considering the number of particles equal to 100 and 500,withσu=10 and σu=1.The comparison results can be seen in Table 2.

Table 2 Experimentalresults ofstate estimation of UNGMwith different parameters

From these results we can see our algorithm is more robust compared with other nonlinear fi lters,where the observation is severely contaminated by the noise.As shown in Table 2,our algorithm expends more time than the other fi lters,since the immune adaptive step is added for the optimization mixture model.

5.2 Nonlinear and non-Gaussian model

For this expreiment,a time-series is generated by the following process model:

where vkis a Gamma Ga(k,θ)random variable modeling the process noise,andω=4e?2 andφ1=0.5 are scalar parameters.A non-stationary observation model is as follows:

whereφ2=0.2 andφ3=0.5.The observation noise nkis drawn from a Gaussian distribution N(0,0.001).Obviously,the process model is non-Gaussian,and the observation modelis nonlinear when k?30.

Different fi lters are used to estimate the underlying state sequence xkgiven the noisy observations yk.All of the particle fi lters use 100 particles and residual resampling. The time step is set as 60.Fig.3 shows the tracking performance of one exemplar run.The superior performance of the IAGMPF is clearly evident.From Fig.4,itis obvious thatthe UKF cannottrack the state wellwhen k?30, since the observation model is nonlinear,meanwhile the particle class fi lters can work better,especially IAGMPF.

Fig.3 Single run result of true trajectory and filtering solutions for 100 particles

Fig.4 Distribution of corresponding estimate error

The experiment is run 100 times with random reinitialization for each run to generate statistical averages with the parameters setabove.Table 3 summarizes the performance ofthe different fi lters.The table shows the means and variances of RMSE of the state estimates.It is obvious that our algorithm is more robustto the non-Gaussian noise distributions,compared with other nonlinear fi lters, where the observation is severely contaminated by the non-Gaussian noise.

Table 3 Experimental results of state estimation of nonlinear and non-Gaussian model

5.3 Two dimension single object tracking with cluttered measurements

Forthe fi rstexample we shallconsidera scenario,in which we are tracking a single target in two dimensions having cluttered measurements.The dynamic model of the objectis described by the discretized Wienerprocess velocity model,where the state of the targetcan be written as

where(xk,yk)is the object's position and(˙xk,˙yk)is the velocity in two dimensionalcartesian coordinates.The discretized dynamics can be expressed with a linear,timeinvariantplantequation

where qk?1is discrete Gaussian white process noise having moments E[qk?1]=0,

The size of the time step is set toΔt=0.1,and the power spectraldensity of process noise to q=0.1.

The likelihood of clutter measurements is de fi ned to be uniformin space[?5,5]×[?4,4].The measurementmodel for the actualtargetis linear with additive Gaussian noise. Thus,we can express the jointmeasurementlikelihood as

where the measurementmodelmatrices are

The system is simulated 240 time steps.Fig.5 shows the simulated(real and clutter)measurements,where the targetis given slightly randomized accelerations such that itachieves a curved trajectory.

Fig.5 True trajectory and distribution of observations and its noises

Fig.6 Distribution of particles of IAGMPF

Fig.7 Filtering results of IAGM PF algorithm using only 30 particles

The distribution of particles of IAGMPF is displayed in Fig.6,where there are 30 particles used for state estimation.In Fig.7 we plot the fi ltering results of different PF algorithms using 30 particles.Moreover,the experimentis run 50 times with random re-initialization for each run to generate statistical averages.Table 4 summarizes the performance ofthe different fi lters.Itcan be seen thatthe estimates ofouralgorithm are signi fi cantly more accurate than other fi lters,even under the strong cluttered measurements and nonlinear objecttracking systems.Meanwhile the executing time of IAGM PF is acceptable,which is a little more than thatof the others.

Table 4 Experimental results of state estimation of two dimension single object tracking

6.Conclusions

In this paper,a hybrid PF is proposed to avoid the impoverishment and degeneracy issue of the generic PF.The improved algorithm consists ofa particle fi lter thatuses an immune adaptive Gaussian mixture algorithm based immune clone algorithm to re-approximate the posteriordensity,and get a closed-form expression for the posterior state distribution,which can effectively reduce the impoverishment effect.Moreover three immune antibody operators are embedded in the new fi lter to substitute for the resample strategy,in which the newest observation and conditionallikelihood are integrated.We validate the algorithm by computersimulations in severalstate estimate applications.Our algorithm outperforms the generic PF and several KFs,with a little more computational complexity than generic PF.In the future,it will be interesting to investigate how to reduce the time expense meanwhile notto prejudice its performance.

[1]W.M.Zhang,G.Du,S.Zhang,et al.Study of nonlinear fi lter methods:particle fi lter.Journal of Systems Engineering and Electronics,2006,17(1):1-5.

[2]M.G.S.Bruno.Sequential Monte Carlo methods for nonlinear discrete-time fi ltering.Synthesis Lectures on Signal Processing,2013,6(1):1-99.

[3]R.H.Zhan,L.Wang,J.W.Wan,etal.Passive targettracking using marginalized particle fi lter.JournalofSystems Engineering and Electronics,2007,18(3):503-508.

[4]P.Chavali,A.Nehorai.Hierarchicalparticle fi ltering formultimodaldata fusion with application to multiple-targettracking. Signal Processing,2014,97(1):207-220.

[5]H.F.Lopes,R.S.Tsay.Particle fi lters and Bayesian inference in fi nancial econometrics.Journal of Forecasting,2011,30: 168-209.

[6]M.Zajac.Online faultdetection of a mobile robotwith a parallelized particle fi lter.Neurocomputing,2014,126(1):151-165.

[7]M.S.Arulampalam,S.Maskell,N.Gordom,etal.A tutorial on particle fi lters for online nonlinear/non-Gaussian Bayesian tracking.IEEE Trans.on Signal Processing,2002,50(2):8-15.

[8]T.Li,S.Sun,T.P.Sattar.Adapting sample size in particle fi lters through KLD-resampling.Electronics Letters,2013, 46(12):740-742.

[9]S.H.Lee,M.West.Convergence of the Markov chain distributed particle fi lter(MCDPF).IEEE Trans.on Signal Processing,2013,61(4):801-812.

[10]A.Banerjee,P.Burlina.Ef fi cient particle fi ltering via sparse kernel density estimation.IEEE Trans.on Image Processing, 2010,19(9):2480-2490.

[11]W.Yi,M.R.Morelande,L.Kong,et al.A computationally ef fi cient particle fi lter for multitarget tracking using an independence approximation.IEEE Trans.on Signal Processing, 2013,61(4):843-856.

[12]N.Whiteley,S.Singh,S.Godsill.Auxiliary particle implementation of probability hypothesis density fi lter.IEEE Trans. on Aerospace and Electronic Systems,2010,46(3):1437-1454.

[13]J.Liu,W.Wang,F.Ma.A regularized auxiliary particle fi ltering approach for system state estimation and battery life prediction.SmartMaterials and Structures,2011,20(7):1-9.

[14]P.M.Stano,Z.Lendek,R.Babuska.Saturated particle fi lter: almost sure convergence and improved resampling.Automatica,2013,49(1):147-159.

[15]T.Cham,J.M.Rehg.Amultiple hypothesis approach to fi gure tracking.Proc.of the IEEE Conference on Computer Vision Pattern Recognition,1999:239-245.

[16]J.Vermaak,A.Doucet,P.Perez.Maintaining multi-modality through mixture tracking.Proc.of the IEEE International Conference on Computer Vision,2003:1110-1116.

[17]B.Han,Y.Zhu,D.Comaniciu,et al.Kernel-based Bayesian fi ltering for object tracking.Proc.of the IEEE Conference Computer Vision and Pattern Recognition,2005:227-234.

[18]Y.R.Wang,X.L.Tang,Q.Cui.Dynamic appearance model for particle fi lter based visual tracking.Pattern Recognition, 2012,45(12):4510-4523.

[19]W.Hassan,N.Bangalore,P.Birch,etal.An adaptive sample count particle fi lter.Computer Vision and Image Understanding,2012,116(12):1208-1222.

[20]Y.Q.Xia,Z.H.Deng,L.Li,etal.A new continuous-discrete particle fi lterfor continuous-discrete nonlinearsystems.Information Sciences,2013,242:64-75.

[21]J.Read,K.Achutegui,J.Miguez.A distributed particle fi lter for nonlinear tracking in wireless sensor networks.Signal Processing,2014,98(1):121-134.

[22]X.Fu,Y.Jia.An improvementon resampling algorithm ofparticle fi lters.IEEE Trans.on Signal Processing,2010,58(10): 5414-5420.

[23]J.Y.Zuo,Y.Liang,Y.Z.Zhang,etal.Particle fi lterwith multimode sampling strategy.Signal Processing,2013,93(11): 3192-3201.

[24]T.Li,S.Sun,T.P.Sattar.Adapting sample size in particle fi lters through KLD-resampling.Electronics Letters,2013, 46(12):740-742.

[25]S.Sarkka.Bayesian filtering and smoothing.Cambridge: Cambridge University Press,2013.

[26]T.Li,T.P.Sattar,S.Sun.Deterministic resampling:unbiased sampling to avoid sample impoverishment in particle fi lters. Signal Processing,2012,92(7):1637-1645.

[27]F.Lindsten.An ef fi cient stochastic approximation EM algorithm using conditional particle fi lters.Proc.of the IEEE International Conference on Acoustics,Speech and Signal Processing,2013:6274-6278.

[28]W.Zhang,G.G.Yen,Z.He.Constrained optimization via arti fi cial immune system.IEEE Trans.on Cybernetics,2014, 44(2):185-198.

[29]W.S.Dong,G.M.Shi,L.Zhang.Immune memory clonalselection algorithms fordesigning stack fi lters.Neurocomputing, 2007,70(4):777-784.

[30]H.Han,Y.S.Ding,K.R.Hao,et al.An evolutionary particle fi lter with the immune genetic algorithm for intelligent video targettracking.Computers and Mathematics with Applications,2011,62(7):2685-2695.

Biographies

Wenlong Huangwas born in 1973.He received his Ph.D.degree from Xidian University,China,in 2009.He is now a postdoctor and lecturer of Air Force Engineering University.His research interests include intelligent information processing and arti fi cial immune system.

E-mail:wljy rghymt@sohu.com

Xiaodan Wangwas born in 1966.She received her Ph.D.degree from Northwestern Polytechnical University,China,in 2000.She is now a professor of Air Force Engineering University.Her research interests include machine learning,pattern recognition and data mining.

E-mail:afeu w@163.com

10.1109/JSEE.2015.00095

Manuscriptreceived April14,2014.

*Corresponding author.

This work was supported by the National Natural Science Foundation of China(61273275;61402517;61573375),the Open Research Fund of State Key Laboratory of Astronautic Dynamics(2012ADL-DW0202),the Natural Science Foundation of ShaanxiProvince of China(2013JQ8035), and the China Postdoctoral Science Foundation(2013M542331).