CN110967706A - 基于改进无迹卡尔曼粒子滤波的raim方法 - Google Patents

基于改进无迹卡尔曼粒子滤波的raim方法 Download PDF

Info

Publication number
CN110967706A
CN110967706A CN201911240626.9A CN201911240626A CN110967706A CN 110967706 A CN110967706 A CN 110967706A CN 201911240626 A CN201911240626 A CN 201911240626A CN 110967706 A CN110967706 A CN 110967706A
Authority
CN
China
Prior art keywords
particle
unscented kalman
filtering
particles
filter
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201911240626.9A
Other languages
English (en)
Inventor
范国超
张依漪
郝创博
薛铸鑫
汪伟峰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Jinghang Computing Communication Research Institute
Original Assignee
Beijing Jinghang Computing Communication Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Jinghang Computing Communication Research Institute filed Critical Beijing Jinghang Computing Communication Research Institute
Priority to CN201911240626.9A priority Critical patent/CN110967706A/zh
Publication of CN110967706A publication Critical patent/CN110967706A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/20—Integrity monitoring, fault detection or fault isolation of space segment

Abstract

本发明属于全球卫星导航检测技术领域，具体涉及一种基于改进无迹卡尔曼粒子滤波的RAIM方法。与现有技术相比较，本发明为实现在非高斯观测噪声环境下有效开展接收机自主完好性监测，通过在粒子滤波方法中引入改进后的无迹卡尔曼滤波，通过改进无迹卡尔曼滤波实现粒子滤波中更合理的状态预测和建议密度计算，能够有效避免粒子退化现象，在保障粒子滤波接收机自主完好性监测方法对非高斯观测噪声环境的适应性的同时，也减少了引入无迹卡尔曼粒子滤波所带来的计算量，保证了算法运算效率。因此，该方案既改善粒子滤波退化问题，同时保障该方法对非高斯观测噪声的适应性，通过对无迹卡尔曼滤波的改进也减少来由其导致的计算量。

Description

基于改进无迹卡尔曼粒子滤波的RAIM方法
技术领域
本发明属于全球卫星导航检测技术领域，具体涉及一种基于改进无迹卡尔曼粒子滤波的RAIM方法。
背景技术
在非高斯观测噪声环境下，最小二乘残差法、奇偶矢量法、卡尔曼滤波等基于高斯噪声假设的传统接收机自主完好性监测方法已无法准确的实现对卫星故障得到检测和识别。利用粒子滤波本身对非高斯观测噪声有较强适应性的优势，构建粒子滤波接收机自主完好性监测方法，能够在非高斯观测噪声环境下开展接收机自主完好性监测(ReceiverAutonomous Integrity Monitoring，RAIM)，但由于粒子滤波方法本身存在粒子退化等问题，监测结果并不理想。
目前针对粒子滤波在非高斯观测噪声环境下应用于接收机自主完好性监测的研究主要集中在改善粒子滤波中的粒子退化问题，目前提出改善方法虽然达到了一定改善粒子退化现象实现多样性的目的，但都在改进过程中导致其对非高斯观测噪声的适应性降低，无法较好的兼顾粒子退化和非高斯观测噪声适应性这两种问题。
传统的粒子滤波接收机自主完好性监测方法具有明显的粒子退化问题，导致估计结果的不准确。其他改进的粒子滤波接收机自主完好性监测方法能减少粒子退化效果现象，但无法保证改正后的粒子滤波接收机自主完好性监测方法在非高斯观测噪声环境的效果。
发明内容
(一)要解决的技术问题
本发明要解决的技术问题是：如何实现在非高斯观测噪声环境下有效开展接收机自主完好性监测。
(二)技术方案
为解决上述技术问题，本发明提供一种基于改进无迹卡尔曼粒子滤波的RAIM方法，所述方法包括如下步骤：
步骤1：对粒子进行初始化抽样；
步骤2：进行状态预测；
步骤3：对粒子集进行序贯重要性采样；
步骤4：对粒子进行重采样；
步骤5：获得状态估计值；
步骤6：对滤波信息进行一致性检验。
其中，所述步骤1中对粒子进行初始化抽样，首先根据先验概率密度函数p(x0)抽取随机样本构成初始粒子集{x0(i)；i＝1,2,…N}～p(x0)，其中N为粒子数目，每个粒子的权重初始化为1/N，表示为{ω0(i)＝1/N；i＝1,2,…N}。
其中，所述步骤2进行状态预测；
通过改进的无迹卡尔曼滤波对粒子状态进行估计，其过程包括时间更新、生成sigma点，观测更新；该步骤2包括：
步骤21：对粒子
Figure BDA0002306111480000021
进行时间更新；
根据公式(3)、(4)及采样得到的粒子集
Figure BDA0002306111480000022
对粒子进行时间更新；
Figure BDA0002306111480000023
Pk/k-1(i)＝Φk/k-1Pk-1(i)ΦT k/k-1+Qk (4)
其中
Figure BDA0002306111480000024
是
Figure BDA0002306111480000025
的最优滤波估计，Φk/k-1为状态转移矩阵下标；k表示当前观测历元，且k＝1,2,3…，k-1表示上一个相邻观测历元，k/k-1表示两个历元之间的过渡，当k＝1时，
Figure BDA0002306111480000026
为步骤一中得到的粒子集；Pk/k-1为
Figure BDA0002306111480000027
的误差协方差矩阵，Qk为过程噪声矩阵；
步骤22：生成sigma点；
根据公式(5)以及上一步更新后的粒子集
Figure BDA0002306111480000031
分别生成三类sigma点
Figure BDA0002306111480000032
Figure BDA0002306111480000033
其中，n为sigma点的总数，Dk|k-1为
Figure BDA0002306111480000034
的平方根，
Figure BDA0002306111480000035
为
Figure BDA0002306111480000036
的协方差矩阵；
步骤23：更新观测信息；
基于步骤22中得到sigma点信息，通过公式(6)、(7)、(8)进行计算得到新的观测信息
Figure BDA0002306111480000037
误差协方差矩阵
Figure BDA0002306111480000038
误差矩阵
Figure BDA0002306111480000039
Figure BDA00023061114800000310
Figure BDA00023061114800000311
Figure BDA00023061114800000312
其中，
Figure BDA00023061114800000313
和
Figure BDA00023061114800000314
为第j个sigma点的权重，其中
Figure BDA00023061114800000315
可通过求一阶统计特性的权系数确定，
Figure BDA00023061114800000316
可通过求二阶统计特性权系数确定，
Figure BDA00023061114800000317
为
Figure BDA00023061114800000318
的误差协方差矩阵，Rk为过程噪声的方差矩阵，
Figure BDA00023061114800000319
为
Figure BDA00023061114800000320
和
Figure BDA00023061114800000321
的误差矩阵，
Figure BDA00023061114800000322
是对每个sigma点的观测更新；
步骤24：更新滤波信息；
基于新的观测信息，依据公式(9)、(10)、(11)更新各项滤波信息；
Figure BDA0002306111480000041
Figure BDA0002306111480000042
Figure BDA0002306111480000043
其中，Kk为滤波增益矩阵，
Figure BDA0002306111480000044
为
Figure BDA0002306111480000045
的误差协方差矩阵。
其中，所述步骤21中，当k＝1时，使用初始粒子集x0。
其中，所述步骤3对粒子集进行序贯重要性采样；该步骤3包括：
步骤31：通过对粒子集
Figure BDA0002306111480000046
做信息统计，计算粒子集的均值
Figure BDA0002306111480000047
和方差
Figure BDA0002306111480000048
Figure BDA0002306111480000049
Figure BDA00023061114800000410
步骤32：通过粒子集的均值
Figure BDA00023061114800000411
和方差
Figure BDA00023061114800000412
以及建议密度函数采样得到新的粒子集xk(i)：
Figure BDA00023061114800000413
其中，所述步骤4对粒子进行重采样；
根据粒子权重
Figure BDA00023061114800000414
计算有效粒子数
Figure BDA00023061114800000415
然后与门限值Nth比较，如果Neff<Nth，则需要对先验粒子集
Figure BDA00023061114800000416
进行重采样得到新的粒子集
Figure BDA00023061114800000417
其中，所述步骤4中的门限值设置为总粒子数的2/3。
其中，所述步骤5获得状态估计值；
根据重采样的粒子集
Figure BDA0002306111480000051
和权重信息
Figure BDA0002306111480000052
计算当前历元的状态估计值：
Figure BDA0002306111480000053
其中，所述步骤6对滤波信息进行一致性检验；所述步骤6包括：
步骤61：计算一致性检验估计值；
根据粒子集
Figure BDA0002306111480000054
和权重信息
Figure BDA0002306111480000055
计算粒子滤波估计值：
Figure BDA0002306111480000056
Figure BDA0002306111480000057
其中，pq(y)为剔除第q个观测分量的粒子滤波估计值，pM(y)为包含所有观测分量的粒子滤波估计值；
步骤62：计算对数似然比；
根据粒子滤波估计值(pq(y)和pM(y))计算对数似然比：
Figure BDA0002306111480000058
其中，si(y)表示第i颗卫星对应的对数似然比；
步骤63：计算累加对数似然比；
考虑到一致性检验对结果的敏感程度，还会对上述来两个估计值的对数似然比累加后再进行对比分析，累加对数似然比
Figure BDA0002306111480000059
的计算过程为：
Figure BDA00023061114800000510
其中，yi|Yi-1为第i个历元下Y观测向量中的观测信息；
步骤64：一致性检验；
通过对比各观测分量的累加对数似然比估计值，可以实现卫星故障的监测和识别；若系统无故障，则各个分量的累加对数似然比基本一致，且随时间的增加，各曲线整体平稳，没有剧烈波动的情况；若某个卫星发生故障，则在该时刻起，该卫星对应的
Figure BDA0002306111480000061
会发生较大的变化，同时与其他卫星对应的估计值具有很大的区别。
其中，在所述方法执行过程中，该方法与直接基于无迹卡尔曼滤波的粒子滤波方法的主要区别是该方法在sigma点生成过程中只需一次时间更新，而传统的UKF中需要对每个sigma点进行时间更新，应用到粒子滤波接收机自主完好性监测中，共有n+1个粒子滤波器，每个粒子滤波中包含N个粒子，并且每个粒子都需要W个sigma点估计，与单纯运用无迹卡尔曼改善粒子滤波接收机自主完好性监测方法相比较，本方法在时间更新步骤中的运算量将是其运算量的1/[(n+1)×N×W]左右。
(三)有益效果
与现有技术相比较，本发明为实现在非高斯观测噪声环境下有效开展接收机自主完好性监测，通过在粒子滤波方法中引入改进后的无迹卡尔曼滤波，通过改进无迹卡尔曼滤波实现粒子滤波中更合理的状态预测和建议密度计算，能够有效避免粒子退化现象，在保障粒子滤波接收机自主完好性监测方法对非高斯观测噪声环境的适应性的同时，也减少了引入无迹卡尔曼粒子滤波所带来的计算量，保证了算法运算效率。因此，该方案既改善粒子滤波退化问题，同时保障该方法对非高斯观测噪声的适应性，通过对无迹卡尔曼滤波的改进也减少来由其导致的计算量。
附图说明
图1为基于改进无迹卡尔曼粒子滤波接收机自主完好性监测流程示意图。
图2为改进无迹卡尔曼粒子滤波流程示意图。
具体实施方式
为使本发明的目的、内容、和优点更加清楚，下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。
为解决现有技术问题，本发明提供一种基于改进无迹卡尔曼粒子滤波的RAIM方法，所述基于改进无迹卡尔曼粒子滤波的接收机自主完好性监测方法通过引入无迹卡尔曼滤波并加以改进，实现更合理的建议密度计算，能够有效避免粒子退化现象，并且保障粒子滤波接收机自主完好性监测方法对非高斯观测噪声环境的适应性，同时与减少了引入无迹卡尔曼粒子滤波所带来的计算量，保证了运算效率，满足接收机完好性监测对粒子滤波方法的需求。
所述方法包括两个阶段，如图1所示，其中第一部分为改进无迹卡尔曼粒子滤波，用于实现各卫星状态观测量的估计，此部分主要步骤如图2所示；另一部分为一致性检验，通过多颗卫星粒子滤波轨迹结果的比较分析卫星故障，达到完好性监测的目的。
所述方法包括如下步骤：
步骤1：对粒子进行初始化抽样；
步骤2：进行状态预测；
步骤3：对粒子集进行序贯重要性采样；
步骤4：对粒子进行重采样；
步骤5：获得状态估计值；
步骤6：对滤波信息进行一致性检验。
其中，所述步骤1中对粒子进行初始化抽样，首先根据先验概率密度函数p(x0)抽取随机样本构成初始粒子集{x0(i)；i＝1,2,…N}～p(x0)，其中N为粒子数目，每个粒子的权重初始化为1/N，表示为{ω0(i)＝1/N；i＝1,2,…N}。
其中，所述步骤2进行状态预测；
通过改进的无迹卡尔曼滤波对粒子状态进行估计，其过程包括时间更新、生成sigma点，观测更新；该步骤2包括：
步骤21：对粒子
Figure BDA0002306111480000081
进行时间更新；
根据公式(3)、(4)及采样得到的粒子集
Figure BDA0002306111480000082
对粒子进行时间更新；
Figure BDA0002306111480000083
Pk/k-1(i)＝Φk/k-1Pk-1(i)ΦT k/k-1+Qk (4)
其中
Figure BDA0002306111480000084
是
Figure BDA0002306111480000085
的最优滤波估计，Φk/k-1为状态转移矩阵下标；k表示当前观测历元，且k＝1,2,3…，k-1表示上一个相邻观测历元，k/k-1表示两个历元之间的过渡，当k＝1时，
Figure BDA0002306111480000086
为步骤一中得到的粒子集；Pk/k-1为
Figure BDA0002306111480000087
的误差协方差矩阵，Qk为过程噪声矩阵；
步骤22：生成sigma点；
根据公式(5)以及上一步更新后的粒子集
Figure BDA0002306111480000088
分别生成三类sigma点
Figure BDA0002306111480000089
Figure BDA00023061114800000810
其中，n为sigma点的总数，Dk|k-1为
Figure BDA00023061114800000811
的平方根，
Figure BDA00023061114800000812
为
Figure BDA00023061114800000813
的协方差矩阵；
步骤23：更新观测信息；
基于步骤22中得到sigma点信息，通过公式(6)、(7)、(8)进行计算得到新的观测信息
Figure BDA00023061114800000814
误差协方差矩阵
Figure BDA00023061114800000815
误差矩阵
Figure BDA0002306111480000091
Figure BDA0002306111480000092
Figure BDA0002306111480000093
Figure BDA0002306111480000094
其中，
Figure BDA0002306111480000095
和
Figure BDA0002306111480000096
为第j个sigma点的权重，其中
Figure BDA0002306111480000097
可通过求一阶统计特性的权系数确定，
Figure BDA0002306111480000098
可通过求二阶统计特性权系数确定，
Figure BDA0002306111480000099
为
Figure BDA00023061114800000910
的误差协方差矩阵，Rk为过程噪声的方差矩阵，
Figure BDA00023061114800000911
为
Figure BDA00023061114800000912
和
Figure BDA00023061114800000913
的误差矩阵，
Figure BDA00023061114800000914
是对每个sigma点的观测更新；
步骤24：更新滤波信息；
基于新的观测信息，依据公式(9)、(10)、(11)更新各项滤波信息；
Figure BDA00023061114800000915
Figure BDA00023061114800000916
Figure BDA00023061114800000917
其中，Kk为滤波增益矩阵，
Figure BDA00023061114800000918
为
Figure BDA00023061114800000919
的误差协方差矩阵。
其中，所述步骤21中，当k＝1时，使用初始粒子集x0。
其中，所述步骤3对粒子集进行序贯重要性采样；该步骤3包括：
步骤31：通过对粒子集
Figure BDA00023061114800000920
做信息统计，计算粒子集的均值
Figure BDA00023061114800000921
和方差
Figure BDA00023061114800000922
Figure BDA00023061114800000923
Figure BDA0002306111480000101
步骤32：通过粒子集的均值
Figure BDA0002306111480000102
和方差
Figure BDA0002306111480000103
以及建议密度函数采样得到新的粒子集xk(i)：
Figure BDA0002306111480000104
其中，所述步骤4对粒子进行重采样；
根据粒子权重
Figure BDA0002306111480000105
计算有效粒子数
Figure BDA0002306111480000106
然后与门限值Nth比较，如果Neff<Nth，则需要对先验粒子集
Figure BDA0002306111480000107
进行重采样得到新的粒子集
Figure BDA0002306111480000108
其中，所述步骤4中的门限值设置为总粒子数的2/3。
其中，所述步骤5获得状态估计值；
根据重采样的粒子集
Figure BDA0002306111480000109
和权重信息
Figure BDA00023061114800001010
计算当前历元的状态估计值：
Figure BDA00023061114800001011
其中，所述步骤6对滤波信息进行一致性检验；所述步骤6包括：
步骤61：计算一致性检验估计值；
根据粒子集
Figure BDA00023061114800001012
和权重信息
Figure BDA00023061114800001013
计算粒子滤波估计值：
Figure BDA00023061114800001014
Figure BDA00023061114800001015
其中，pq(y)为剔除第q个观测分量的粒子滤波(称为辅助粒子滤波)估计值，pM(y)为包含所有观测分量的粒子滤波(称为主粒子滤波)估计值；
步骤62：计算对数似然比；
根据粒子滤波估计值(pq(y)和pM(y))计算对数似然比：
Figure BDA0002306111480000111
其中，si(y)表示第i颗卫星对应的对数似然比；
步骤63：计算累加对数似然比；
考虑到一致性检验对结果的敏感程度，通常还会对上述来两个估计值的对数似然比累加后再进行对比分析，累加对数似然比
Figure BDA0002306111480000112
的计算过程为：
Figure BDA0002306111480000113
其中，yi|Yi-1为第i个历元下Y观测向量中的观测信息；
步骤64：一致性检验；
通过对比各观测分量的累加对数似然比估计值，可以实现卫星故障的监测和识别；若系统无故障，则各个分量的累加对数似然比基本一致，且随时间的增加，各曲线整体平稳，没有剧烈波动的情况；若某个卫星发生故障，则在该时刻起，该卫星对应的
Figure BDA0002306111480000114
会发生较大的变化，同时与其他卫星对应的估计值具有很大的区别。
此时粒子滤波接收机自主完好性过程如图2。
由于上述过程需要对粒子滤波结果进行一致性检验，则每颗卫星的观测分量都需要一次辅助粒子滤波。假设在k时刻存在M颗可见星，则需要进行M+1次粒子滤波才能完成完好一次完好性监测，这要求应用在完好性监测中的粒子滤波方法需具有较高的计算效率。
在所述方法执行过程中，该方法与直接基于无迹卡尔曼滤波的粒子滤波方法的主要区别是该方法在sigma点生成过程中只需一次时间更新，而传统的UKF中需要对每个sigma点进行时间更新，应用到粒子滤波接收机自主完好性监测中，共有n+1个粒子滤波器，每个粒子滤波中包含N个粒子，并且每个粒子都需要W个sigma点估计，与单纯运用无迹卡尔曼改善粒子滤波接收机自主完好性监测方法相比较，本方法在时间更新步骤中的运算量将是其运算量的1/[(n+1)×N×W]左右。
综上，本发明属于全球卫星导航检测技术领域，具体涉及一种基于改进无迹卡尔曼粒子滤波的RAIM方法。所述方法涉及定位解算、观测数据处理、完好性监测、无迹卡尔曼滤波、粒子滤波等相关技术，主要针对在复杂观测噪声环境下，卫星导航定位服务过程中可能出现的观测数据故障等问题，通过采用改进的无机卡尔曼滤波与粒子滤波进行组合，改进传统粒子滤波粒子退化问题，同时兼顾非高斯观测噪声环境和高效率的应用场景，给出一种有效的全球卫星导航系统接收机自主完好性监测监测方法。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (10)

1.一种基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述方法包括如下步骤：
步骤1：对粒子进行初始化抽样；
步骤2：进行状态预测；
步骤3：对粒子集进行序贯重要性采样；
步骤4：对粒子进行重采样；
步骤5：获得状态估计值；
步骤6：对滤波信息进行一致性检验。
2.如权利要求1所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述步骤1中对粒子进行初始化抽样，首先根据先验概率密度函数p(x0)抽取随机样本构成初始粒子集{x0(i)；i＝1,2,…N}～p(x0)，其中N为粒子数目，每个粒子的权重初始化为1/N，表示为{ω0(i)＝1/N；i＝1,2,…N}。
3.如权利要求2所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述步骤2进行状态预测；
通过改进的无迹卡尔曼滤波对粒子状态进行估计，其过程包括时间更新、生成sigma点，观测更新；该步骤2包括：
步骤21：对粒子
Figure FDA0002306111470000011
进行时间更新；
根据公式(3)、(4)及采样得到的粒子集
Figure FDA0002306111470000012
对粒子进行时间更新；
Figure FDA0002306111470000013
Pk/k-1(i)＝Φk/k-1Pk-1(i)ΦT k/k-1+Qk (4)
其中
Figure FDA0002306111470000014
是
Figure FDA0002306111470000015
的最优滤波估计，Φk/k-1为状态转移矩阵下标；k表示当前观测历元，且k＝1,2,3…，k-1表示上一个相邻观测历元，k/k-1表示两个历元之间的过渡，当k＝1时，
Figure FDA0002306111470000021
为步骤一中得到的粒子集；Pk/k-1为
Figure FDA0002306111470000022
的误差协方差矩阵，Qk为过程噪声矩阵；
步骤22：生成sigma点；
根据公式(5)以及上一步更新后的粒子集
Figure FDA0002306111470000023
分别生成三类sigma点
Figure FDA0002306111470000024
Figure FDA0002306111470000025
其中，n为sigma点的总数，Dk|k-1为
Figure FDA0002306111470000026
的平方根，
Figure FDA0002306111470000027
为
Figure FDA0002306111470000028
的协方差矩阵；
步骤23：更新观测信息；
基于步骤22中得到sigma点信息，通过公式(6)、(7)、(8)进行计算得到新的观测信息
Figure FDA0002306111470000029
误差协方差矩阵
Figure FDA00023061114700000210
误差矩阵
Figure FDA00023061114700000211
Figure FDA00023061114700000212
Figure FDA00023061114700000213
Figure FDA00023061114700000214
其中，
Figure FDA00023061114700000215
和
Figure FDA00023061114700000216
为第j个sigma点的权重，其中
Figure FDA00023061114700000217
可通过求一阶统计特性的权系数确定，
Figure FDA00023061114700000218
可通过求二阶统计特性权系数确定，
Figure FDA00023061114700000219
为
Figure FDA00023061114700000220
的误差协方差矩阵，Rk为过程噪声的方差矩阵，
Figure FDA00023061114700000221
为
Figure FDA0002306111470000031
和
Figure FDA0002306111470000032
的误差矩阵，
Figure FDA0002306111470000033
是对每个sigma点的观测更新；
步骤24：更新滤波信息；
基于新的观测信息，依据公式(9)、(10)、(11)更新各项滤波信息；
Figure FDA0002306111470000034
Figure FDA0002306111470000035
Figure FDA0002306111470000036
其中，Kk为滤波增益矩阵，
Figure FDA0002306111470000037
为
Figure FDA0002306111470000038
的误差协方差矩阵。
4.如权利要求3所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述步骤21中，当k＝1时，使用初始粒子集x0。
5.如权利要求3所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述步骤3对粒子集进行序贯重要性采样；该步骤3包括：
步骤31：通过对粒子集
Figure FDA0002306111470000039
做信息统计，计算粒子集的均值
Figure FDA00023061114700000310
和方差
Figure FDA00023061114700000311
Figure FDA00023061114700000312
Figure FDA00023061114700000313
步骤32：通过粒子集的均值
Figure FDA00023061114700000314
和方差
Figure FDA00023061114700000315
以及建议密度函数采样得到新的粒子集xk(i)：
Figure FDA00023061114700000316
6.如权利要求5所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述步骤4对粒子进行重采样；
根据粒子权重
Figure FDA0002306111470000041
计算有效粒子数
Figure FDA0002306111470000042
然后与门限值Nth比较，如果Neff<Nth，则需要对先验粒子集
Figure FDA0002306111470000043
进行重采样得到新的粒子集
Figure FDA0002306111470000044
7.如权利要求6所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述步骤4中的门限值设置为总粒子数的2/3。
8.如权利要求6所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述步骤5获得状态估计值；
根据重采样的粒子集
Figure FDA0002306111470000045
和权重信息
Figure FDA0002306111470000046
计算当前历元的状态估计值：
Figure FDA0002306111470000047
9.如权利要求8所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，所述步骤6对滤波信息进行一致性检验；所述步骤6包括：
步骤61：计算一致性检验估计值；
根据粒子集
Figure FDA0002306111470000048
和权重信息
Figure FDA0002306111470000049
计算粒子滤波估计值：
Figure FDA00023061114700000410
Figure FDA00023061114700000411
其中，pq(y)为剔除第q个观测分量的粒子滤波估计值，pM(y)为包含所有观测分量的粒子滤波估计值；
步骤62：计算对数似然比；
根据粒子滤波估计值(pq(y)和pM(y))计算对数似然比：
Figure FDA0002306111470000051
其中，si(y)表示第i颗卫星对应的对数似然比；
步骤63：计算累加对数似然比；
考虑到一致性检验对结果的敏感程度，还会对上述来两个估计值的对数似然比累加后再进行对比分析，累加对数似然比
Figure FDA0002306111470000052
的计算过程为：
Figure FDA0002306111470000053
其中，yi|Yi-1为第i个历元下Y观测向量中的观测信息；
步骤64：一致性检验；
通过对比各观测分量的累加对数似然比估计值，可以实现卫星故障的监测和识别；若系统无故障，则各个分量的累加对数似然比基本一致，且随时间的增加，各曲线整体平稳，没有剧烈波动的情况；若某个卫星发生故障，则在该时刻起，该卫星对应的
Figure FDA0002306111470000054
会发生较大的变化，同时与其他卫星对应的估计值具有很大的区别。
10.如权利要求9所述的基于改进无迹卡尔曼粒子滤波的RAIM方法，其特征在于，在所述方法执行过程中，该方法与直接基于无迹卡尔曼滤波的粒子滤波方法的主要区别是该方法在sigma点生成过程中只需一次时间更新，而传统的UKF中需要对每个sigma点进行时间更新，应用到粒子滤波接收机自主完好性监测中，共有n+1个粒子滤波器，每个粒子滤波中包含N个粒子，并且每个粒子都需要W个sigma点估计，与单纯运用无迹卡尔曼改善粒子滤波接收机自主完好性监测方法相比较，本方法在时间更新步骤中的运算量将是其运算量的1/[(n+1)×N×W]左右。
CN201911240626.9A 2019-12-06 2019-12-06 基于改进无迹卡尔曼粒子滤波的raim方法 Pending CN110967706A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911240626.9A CN110967706A (zh) 2019-12-06 2019-12-06 基于改进无迹卡尔曼粒子滤波的raim方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911240626.9A CN110967706A (zh) 2019-12-06 2019-12-06 基于改进无迹卡尔曼粒子滤波的raim方法

Publications (1)

Publication Number Publication Date
CN110967706A true CN110967706A (zh) 2020-04-07

Family

ID=70033303

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911240626.9A Pending CN110967706A (zh) 2019-12-06 2019-12-06 基于改进无迹卡尔曼粒子滤波的raim方法

Country Status (1)

Country Link
CN (1) CN110967706A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112558115A (zh) * 2020-11-30 2021-03-26 中航机载系统共性技术有限公司 一种基于自适应bfo-pso改进粒子滤波的卫星raim监测方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105004351A (zh) * 2015-05-14 2015-10-28 东南大学 基于自适应upf的sins大方位失准角初始对准方法
CN110455287A (zh) * 2019-07-24 2019-11-15 南京理工大学 自适应无迹卡尔曼粒子滤波方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105004351A (zh) * 2015-05-14 2015-10-28 东南大学 基于自适应upf的sins大方位失准角初始对准方法
CN110455287A (zh) * 2019-07-24 2019-11-15 南京理工大学 自适应无迹卡尔曼粒子滤波方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
冉星浩 等: "基于无迹卡尔曼滤波和权值优化的改进粒子滤波算法", 《探测与控制学报》 *
彭雅奇 等: "基于鲁棒扩展卡尔曼粒子滤波的RAIM算法", 《系统工程与电子技术》 *
朱立新 等: "一种基于UKF改进算法的组合导航系统故障的检测方法", 《计算机应用研究》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112558115A (zh) * 2020-11-30 2021-03-26 中航机载系统共性技术有限公司 一种基于自适应bfo-pso改进粒子滤波的卫星raim监测方法
CN112558115B (zh) * 2020-11-30 2023-10-10 中航机载系统共性技术有限公司 一种基于自适应bfo-pso改进粒子滤波的卫星raim监测方法

Similar Documents

Publication Publication Date Title
Cai et al. Cycle slip detection and repair for undifferenced GPS observations under high ionospheric activity
CN110007317B (zh) 一种选星优化的高级接收机自主完好性监测方法
US10197678B1 (en) H-ARAIM system of optimizing a horizontal protection level
CN110646820B (zh) Rtk定位数据的质检方法、装置、设备和存储介质
Zhou et al. Critical issues on Kalman filter with colored and correlated system noises
CN104375157B (zh) 短基线下惯导辅助的北斗单频整周模糊度求解方法
CN110006427B (zh) 一种低动态高振动环境下的bds/ins紧组合导航方法
Langley An optimized least-squares technique for improving ambiguity resolution and computational efficiency
JP2011522269A (ja) 異常な擬似距離測定値から無線ナビゲーション受信機ユーザを保護するための方法
CN107505636A (zh) 海基jpals的定位域mrcc方法及装置
CN108761498B (zh) 一种针对高级接收机自主完好性监测的位置估计优化方法
CN113670337A (zh) 一种用于gnss/ins组合导航卫星缓变故障检测方法
Zhu et al. Extended Kalman filter (EKF) innovation-based integrity monitoring scheme with C/N 0 weighting
CN110967706A (zh) 基于改进无迹卡尔曼粒子滤波的raim方法
CN101630000B (zh) 一种评估干扰信号对gps性能影响的系统
CN105511481B (zh) 一种星载定轨优化方法
Marnissi et al. An auxiliary variable method for Langevin based MCMC algorithms
Green et al. Position-domain integrity analysis for generalized integer aperture bootstrapping
CN110988921A (zh) 基于改进无迹卡尔曼粒子滤波的raim系统
CN106506008B (zh) 一种基于结构化测量矩阵的块稀疏信号恢复方法
CN116495204A (zh) 一种基于异动概率的轨道异动检测方法及系统
CN104035112B (zh) 利用虚拟高程模型辅助城市环境下卫星导航定位的方法
Hwang et al. NIORAIM integrity monitoring performance in simultaneous two-fault satellite scenarios
Khanafseh et al. Integrity risk of cycle resolution in the presence of bounded faults
Renaux et al. The Bayesian Abel bound on the mean square error

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication

Application publication date: 20200407

RJ01 Rejection of invention patent application after publication