基于自适应噪声方差的卫星定位故障检测法_陈含智.pdf

http:/DOI:10.13700/j.bh.1001-5965.2021.0222基于自适应噪声方差的卫星定位故障检测法陈含智1，孙蕊1，*，邱明1，毛继志2，胡浩亮2，张立东2(1.南京航空航天大学民航学院，南京211106；2.中国航空无线电电子研究所民航空管航空电子技术实验室，上海201109)摘要：不同观测环境的实际观测噪声可能存在较大差别，因此，固定的噪声方差矩阵可能使基于卡尔曼滤波故障检测法的故障检测效果下降。对此，提出自适应噪声方差的卫星定位故障检测算法，利用滑动窗内的历史新息实时估计观测噪声方差矩阵，在此基础上构建故障检测量与识别量对故障进行检测与识别，利用不含故障的新息更新状态向量并解算出定位结果。实验结果表明：所提算法在静态模式下，可 100%检测和识别的最小单阶跃故障为 3m，对持续时间为 100s、斜坡故障增长速率为 0.2m/s 的斜坡故障识别比率为 51.4%，可 100%检测和识别的最小多故障为 4m；在动态模式下，可 100%检测和识别的最小单阶跃故障为 10m，对持续时间为 200s、斜坡故障增长速率为 0.2m/s 的斜坡故障识别比率为 66.25%，可 100%检测和识别的最小多故障为 12m。所提算法性能优于基于最小二乘和卡尔曼滤波的故障检测法。关键词：全球卫星导航系统；全球卫星导航系统质量控制；故障检测；卡尔曼滤波；自适应滤波中图分类号：V249.3文献标志码：A文章编号：1005-5965（2023）02-0406-16全球卫星导航系统(globalnavigationsatellitesys-tem，GNSS)有着广泛的用途，在军事和民用领域发挥着重要作用。民航、自动驾驶等交通领域对 GNSS观测数据的质量有较高的要求1。然而，观测环境容易造成含有较大误差的观测数据，即故障数据，如环境反射的卫星信号造成的多路径效应2。GNSS质量控制旨在保证复杂环境下 GNSS 的观测数据质量。故障检测算法作为 GNSS 质量控制的重要组成部分，通过检测和排除故障观测数据获得高精度可靠的定位结果，保障 GNSS 安全地应用于民航、自动驾驶等交通领域。目前为止，国内外学者已提出多种卫星定位故障检测算法。其中，Lee3提出的伪距比较法、Parkinson 和 Axelrad4提出的最小二乘残差法和 Sturza5提出的奇偶矢量法都是对当前观测值进行一致性检验，被称为“快照法”6。一些学者又提出基于历史信息和当前信息一致性检验的“滤波法”。例如，针对“快照法”难以检测微小误差7，沙海等8提出一种基于滑动窗内多个历元残差矢量累积的抗差卡尔曼滤波方法。Bhatta-charyya 和 Gebre-Egziabher9提出基于卡尔曼滤波的故障检测与保护水平计算方法。针对卡尔曼滤波不适用于非高斯问题10，王尔申等11-12又提出基于粒子滤波的故障检测法。刘江等13提出了适用于非线性观测特征的基于容积卡尔曼滤波的故障检测法。吴云14比较了“快照法”和“滤波法”，通过基于静态和动态数据的实验结果得出，“滤波法”的检测和识别性能优于“快照法”的结论。基于卡尔曼滤波的故障检测法中，观测噪声方差矩阵影响故障检测量和卡尔曼增益矩阵的确定，进一步影响故障检测效果和滤波结果的精度，因此，准确确定观测噪声方差矩阵具有较大的意义。不同观测条件下的观测噪声可能有较大差别，例如收稿日期：2021-05-06；录用日期：2021-08-09；网络出版时间：2021-09-1415：20网络出版地址： J.北京航空航天大学学报，2023，49（2）：406-421.CHEN H Z，SUN R，QIU M，et al.An adaptive noise variance based fault detection algorithm for GNSS positioningJ.Journal of BeijingUniversity of Aeronautics and Astronautics，2023，49（2）：406-421（in Chinese）.2023年2月北京航空航天大学学报February2023第49卷第2期JournalofBeijingUniversityofAeronauticsandAstronauticsVol.49No.2高速无人机上搭载接收机的观测噪声应明显大于静态接收机的观测噪声。因此，按照经验确定观测噪声方差矩阵并将其固定的做法并不合理，有可能造成滤波结果精度下降，甚至故障漏检和误检。为改进基于卡尔曼滤波的故障检测法将观测噪声方差矩阵固定的策略，本文提出基于自适应噪声方差的卫星定位故障检测算法，利用滑动窗内的历史新息序列实时估计观测噪声方差矩阵，再重构故障检测量与识别量，与根据误警率确定的阈值比较来检测和识别故障，完成故障检测后计算卡尔曼增益矩阵，最后利用排除故障后的新息数据进行测量更新并输出滤波的状态向量作为定位结果。1本文算法1.1基于自适应噪声方差的卫星定位故障检测算法框架本文所提算法的框架如图 1 所示：首先初始化滤波参数，预测下一历元的状态向量，再计算新息和观测噪声方差矩阵，接着计算故障检测量并通过与故障检测阈值比较以判断是否存在故障数据，若存在故障数据，则进行故障识别并排除故障数据，完成故障检测与识别后，计算该历元状态向量的最终估计值，输出该历元的定位结果。1.2状态方程和测量方程k设 历元状态向量为Xk=xk,yk,zk,rk,xk,yk,zk,rkT（1）xkykzkkrkkxkykzkkkrkk式中：、和分别为 历元时接收机在地心地固坐标系中的 x 轴、y 轴和 z 轴坐标；为 历元的接收机钟误差与光速的乘积；、和分别为 历元时接收机在地心地固坐标系中的 x 轴、y 轴和 z 轴坐标对时间的导数，即接收机 历元时在 x 轴、y 轴和z 轴方向上的速度分量；为 历元时接收机钟误差与光速的乘积对时间的导数。建立状态方程为Xk+1=k+1|kXk+Wk（2）Xk+1k+1WkWkQk+1|k式中：为历元的状态向量；为过程噪声，假设为高斯白噪声，其协方差矩阵为；为状态转移矩阵，其式为k+1|k=I44tI44044I44（3）I44tkk+1044式中：为 4 阶单位矩阵；为 历元与历元时间间隔；为 4 阶零矩阵。建立观测方程如下：Zk+1=Hk+1Xk+1+Vk+1（4）Zk+1Hk+1Vk+1Vk+1Rk+1式中：为观测向量；为观测矩阵；为伪距观测噪声，假设为高斯白噪声，其协方差矩阵为。观测向量为Zk+1=k+1,1?k+1,1+mk+1,1xk+nk+1,1yk+lk+1,1zkk+1,2?k+1,2+mk+1,2xk+nk+1,2yk+lk+1,2zk.k+1,Nk+1?k+1,Nk+1+mk+1,Nk+1xk+nk+1,Nk+1yk+lk+1,Nk+1zk（5）Nk+1k+1k+1,i(i=1,2,Nk+1)k+1i?k+1,i(i=1,2,Nk+1)(xk+1,yk+1,zk+1)(xk,yk,zk)k+1imk+1,i、nk+1,i、lk+1,i(i=1,2,Nk+1)(xk,yk,zk)k+1i式 中：为历 元 可 见 星 的 数 量；为历元可见卫星 经过误差改正后的伪距；为的近似位置到历元的可见星 的距离；分别为到历元可见星 的方向余弦。观测矩阵如下：Hk+1=Dk+10Nk+14（6）Dk+1=mk+1,1nk+1,1lk+1,11mk+1,2nk+1,2lk+1,21.mk+1,Nk+1nk+1,Nk+1lk+1,Nk+110Nk+14Nk+14式中：；为维零矩阵。计算k+1历元的故障检测量k+1历元的故障检测量大于其阈值?NNNNYYYY计算卫星i的故障识别统计量卫星i故障识别统计量大于其阈值?初始化滤波参数标记卫星i的伪距含故障标记卫星i的伪距不含故障预测际k+1历元的状态向量完成可见星遍历?i=i+1计算新息根据故障检测与识别的结果调整新息计算观测噪声协方差矩阵更新k+1历元的状态向量，解算定位结果k=k+1继续处理下一历元的数据?结束图1本文算法框架Fig.1Frameworkoftheproposedalgorithm第2期陈含智，等：基于自适应噪声方差的卫星定位故障检测法4071.3观测噪声方差矩阵自适应滤波的一个目的是利用观测信息，实时估计和修正噪声统计特性或增益矩阵15。本文通过历史新息序列实时估计观测噪声方差矩阵。k+1根据状态方程，预测历元的状态向量为Xk+1|k=k+1|kXk（7）Xk+1|kk+1Xkk式中：为历元状态向量的一步预测值；为 历元状态向量的滤波估计值。Xk+1|k的误差的协方差矩阵为Pk+1|k=k+1|kPkTk+1|k+Q（8）PkXk式中：为估计值的误差协方差矩阵。计算新息如下：k+1=Zk+1Hk+1Xk+1|k（9）理论上新息的协方差矩阵为Ck+1=Hk+1Pk+1|kHTk+1+Rk+1（10）Ck+1L根据极大似然准则，在长度为 的滑动窗（假设滑动窗内观测的可见星相同）内的最优估计值为16Ck+1=1Lkj=j0jTj（11）j0=kL+1式中：。为使滑动窗内最新的新息对新息协方差的估计起到更大的影响，对式(11)做如下调整：Ck+1=2L(L+1)kj=j0(jk+L)jTj（12）Rk+1于是，的估计值为Rk+1=2L(L+1)kj=j0(jk+L)jTjHk+1Pk+1|kHTk+1（13）Rk+1假设各颗卫星的伪距观测噪声互不相关，同时为了避免估计结果超出合理范围，对估计结果进行如下约束，得到修正后的伪距观测噪声协方差矩阵为Rk+1(p,q)=Rk+1(p,q)2Rk+1(p,q)20 1,p=q120Rk+1(p,q)20 1,p=q220Rk+1(p,q)20 TD于是可以得到故障检测准则：若，则该历元无故障发生；若，则该历元有故障发生。1.5故障识别若检测到有故障，则进一步识别故障。构建故障识别量为k+1,i=k+1(i)Ck+1(i,i)（19）k+1(i)k+1iCk+1(i,i)Ck+1iii=1,2,Nk+1式中：为新息向量的第 个卫星；为新息协方差矩阵的第 行第 列卫星，。k+1,iik+1,i TRik+1,iTRiTRk+1,iTR统计量在卫星 的观测数据无故障时服从标准正态分布。于是可以得到如下故障识别准则：若，则卫星 的伪距不含故障；若,则卫星 的伪距含故障。其中，为统计量的阈值，可根据式（20）确定：+TR12et22dt=PFA2Nk+1（20）e式中：为自然底数。1.6算法步骤本节给出本文所提算法的详细步骤。k=0X0P0、Q、L、PFA、20、1、2步骤1对滤波参数进行初始化：令，确定初始状态向量，根据经验确定。k+1Xk+1|kPk+1|k步骤2根据式(7)预测历元的状态向量，根据式(8)计算。k+1步骤3根据式(9)计算新息。Rk+1k+1步骤4确定观测噪声方差矩阵。若历408北 京 航 空 航 天 大 学 学 报2023年LRk+1Rk+1=20INk+1Nk+1Ck+1元的所有可见星对应滑动窗内的新息序列的长度均等于，根据式(14)计算，否则。根据式(15)计算。步骤5根据式(16)计算故障检测量，并根据故障检测准则判断该历元是否存在故障数据。步骤6若判断该历元不存在故障数据，则直接进入步骤 8；否则根据式(19)计算各颗卫星的故障识别量，并根据故障识别准则识别故障卫星。k+1步骤7对进行如下调整：k+1(i)=Ck+1(i,i)k+1,i TRk+1(i)k+1,i TR（21）k+1k+1(i)k+1ii=1,2,Nk+1式中：为调整后的新息序列；为向量的第 个卫星，。k+1Xk+1步骤8更新历元的状态向量，输出定位结果：增益矩阵为Kk+1=Pk+1|kHTk+1(Ck+1)1（22）k+1于是测量更新后的历元状态向量为Xk+1=Xk+1|k+Kk+1k+1（23）Xk+1计算误差的协方差矩阵Pk+1=(I88Kk+1Hk+1)Pk+1|k（2