2210 仪器仪表学报 第30卷 其中，y1、y2、y3、4、5、6为当前时刻的测量值. 粒子集中的粒子数。同时，在保证一定的估计精度下， 粒子滤波器MCMC-P℉JJ=B,C,D,E,F,G的输入为测量子 适当减少粒子的个数以及时间窗的长度。 集。{，pOJ=BCD,E,F.G为粒子滤波器给出的状态 估计和似然概率密度。详细的算法为： 41 算法验证与实验结果分析 1)初始化，0： 根据先验概率密度(X。)产生N个初始化粒子 1)实验条件 (K”:i=12,N-pX)，并设置初始化权值为 利用实验室的NewStar100GPS/GNSS接收机研发平 ,=1/N。 台提取实验数据，观测的数据包括用于接收机位置解算 2)当k=1,2,…,执行下面的步骤： 的卫星位置坐标和卫星到用户之间的测量伪距。静态采 ①状态预测：根据前面产生的粒子以及状态方程产 集剥量数据3508，采集数据的时间为第498周中的第 生新的粒子{X-(:i=l2N~pX1X)。 293249~293599s(2009年3月11日17时27分30秒)，在 ②计算似然概率密度：在获得测量值y:后，根据测量 此期间，可见的卫星为7颗，跟踪解算卫星为6颗。卫星 值构成测量全集和测量子集。首先，计算状态观测似然 D号分别为5、9、12、18、22、30.对应的测量值分别为 概率密度： Y=(yl,y2,y3,y4,y5,y6)。同时，并行利用u-blox公司的 w(i)=p(y IXu-()).(i=1..N:j=1..m) RCB-4H模块监测得到此时这些卫星无故障，处于正常状 然后，分别求出各个MCMC粒子滤波器的似然联合 态。为了验证某颗卫星发生故障时，算法是否能够对故 概率密度。假设观测时刻k测量值之间彼此统计独立，因 障卫星进行正确检测，人为地加入偏差，这里，对用于解 此，各个MCMC粒子滤波器的似然联合概率密度可以表 算的第5号卫星在=180的观测时刻注入偏差，即： 示为测量子集中各个测量值的状态测量似然概率密度的 0,t=0-180(s) b= (4) 乘积，即为： 5m,t=180~350(s) p(Y IXu()),J=A.B.C.D.E.F.G 在实验中，粒子数选为900个，计算检验统计量的窗 3)计算对数似然比： 口长度选为40。实验所用数据是从实际的接收机中采集 R5)=hpW1Xa✉) 得到的，其中观测噪声服从高斯核拉普拉斯分布 pY1X-()' J=B,C.D.E,F.G (Gaussian core-Laplacian)例。 4)累加对数似然比： 2)实验结果与分析 S-R)()J-B.C.D.E.F.G 图2给出了利用卡尔曼滤波和MCMC PF粒子滤波的 检验统计量，图3给出了正常条件下的MCMC粒子滤波完 式中：W是累加时间窗长。 好性监测实验结果，图4给出了人为加入故障后的MCMC 5)计算检验统计量： 粒子滤波完好性监测实验结果。 Y=max Sum_LLR(J),J=B.C,D.E,F,G 5004 6)故障检测与隔离： 450叶 一是鼓 如果Y>T,T为检测门限，则给出告警，同 400 时，通过故障观测子集给出故障卫星并进行隔 350 300 离：否则，系统工作正常。 7)重采样： 200 依据重要性权重()=w'() /乞0对粒子 华150 100 进行重采样，采样后的粒子集为(XO:i=l,N)。 50 8)M-H抽样： 0 30100150200250300350 利用MCMC方法对式(7)中得到的粒子进行M-H采 时间s 样，得到新的粒子集(XO,1/N:i=1,,W)。 图2卡尔曼滤波和MCMC粒子滤波检验统计量 为了减少粒子滤波算法的运算量，算法中采用状态变 Fig.2 Test statistics for Kalman filter and MCMC PF 量的转移概率密度函数作为重要性概率密度函数，并减少 万方数据2210 仪器仪表学报 第30卷 其中，y1、y2、)，3、y4、y5、))6为当前时刻的测量值。 粒子滤波器MCMC．PF J(J=B，C,D，E，EG)的输入为测量子 集。∥，pJ()，)；‘，=且c，D,E,F，G1为粒子滤波器给出的状态 估计和似然概率密度。详细的算法为： 1)初始化，k=O： 根据先验概率密度p(X。)产生Ⅳ个初始化粒子 {Xo“’：f=1．2’…，N}-烈Xo)，并设置初始化权值为 ％‘‘’=I／N。 2)当k=l，2，…，执行下面的步骤： ①状态预测：根据前面产生的粒子以及状态方程产 生新的粒子(鼍¨∞：f=l'2'…，Ⅳ卜p(邑I邑_)。 ②计算似然概率密度：在获得测量值岁。后，根据测量 值构成测量全集和测量子集。首先，计算状态观测似然 概率密度： 魄砸)=p魄u’I五¨∞)，“=卜·，Ⅳ；．『=l，o o功 然后，分别求出各个MCMC粒子滤波器的似然联合 概率密度。假设观测时刻忌测量值之间彼此统计独立，因 此，各个MCMC粒子滤波器的似然联合概率密度可以表 示为测量子集中各个测量值的状态测量似然概率密度的 乘积，即为： p(L¨’I Xklk-I(f))，J=A，E c’D，EFG 3)计算对数似然比： 峭‘(肛h面p(Yk硒CJ)l Xk,k_,(i))，．，=E c，D，E只G 4)累加对数似然比： Sum—LLR(J)=∑叫‘(．，)，JffiB，C，D，E，F，G 式中：w是累加时间窗长。 5)计算检验统计量： 圪=maxSum—L酞(．，)，J=B，C，D，E。F，G 6)故障检测与隔离： 如果托>T，哟检测门限，则给出告警，同 时，通过故障观测子集给出故障卫星并进行隔 离；否则，系统工作正常。 7)重采样： 依据重要性权重矾’(f)=峨．f(f)／∑'．，。砸)对粒子 进行重采样，采样后的粒子集为{夏。“’：汪l'…，N}。 8)M—H抽样： 利用MCMC方法对式(7)中得到的粒子进行M．H采 样，得到新的粒子集{X。“’，llN；i=1，…，N}。 为了减少粒子滤波算法的运算量，算法中采用状态变 量的转移概率密度函数作为重要性概率密度函数，并减少 粒子集中的粒子数。同时，在保证一定的估计精度下， 适当减少粒子的个数以及时间窗的长度。 4算法验证与实验结果分析 1)实验条件 利用实验室的NewStarl00 GPS／GNSS接收机研发平 台提取实验数据，观测的数据包括用于接收机位置解算 的卫星位置坐标和卫星到用户之间的测量伪距。静态采 集测量数据350 s，采集数据的时间为第498周中的第 293 249～293 599 s(2009年3月11日17时27分30秒)，在 此期间，可见的卫星为7颗，跟踪解算卫星为6颗。卫星 ID号分别为5、9、12、18、22、30，对应的测量值分别为 y=(yl，y2，y3，y4，y5，y6)。同时，并行利用U-blox公司的 RCB．4H模块监测得到此时这些卫星无故障，处于正常状 态。为了验证某颗卫星发生故障时，算法是否能够对故 障卫星进行正确检测，人为地加入偏差，这里，对用于解 算的第5号卫星在t=180的观测时刻注入偏差，即： ， 10，t=0—180 (s) … 驴lsm，：180。3505m t 180 350(s) 。 竹’ 。= 一 (s) ～ 在实验中，粒子数选为900个，计算检验统计量的窗 口长度选为40。实验所用数据是从实际的接收机中采集 得到的，其中观测噪声服从高斯核拉普拉斯分布 (Gaussian corc．Laplacian)[91。 2)实验结果与分析 图2给出了利用卡尔曼滤波和MCMC PF粒子滤波的 检验统计量，图3给出了正常条件下的MCMC粒子滤波完 好性监测实验结果，图4给出了人为加入故障后的MCMC 粒子滤波完好性监测实验结果。 窦 蠕 蠡 时间，暑 图2卡尔曼滤波和MCMC粒子滤波检验统计量 Fig．2 Test statistics for Kalman filter and MCMC PF 万方数据