技术领域
[0001] 本
发明属于全球卫星导航检测技术领域，具体涉及一种基于改进无迹卡尔曼粒子滤波的RAIM系统。
背景技术
[0002] 在非高斯观测噪声环境下，最小二乘残差法、奇偶矢量法、卡尔曼滤波等基于高斯噪声假设的传统接收机自主完好性监测方法已无法准确的实现对卫星故障得到检测和识别。利用粒子滤波本身对非高斯观测噪声有较强适应性的优势，构建粒子滤波接收机自主完好性监测方法，能够在非高斯观测噪声环境下开展接收机自主完好性监测(Receiver Autonomous Integrity Monitoring，RAIM)，但由于粒子滤波方法本身存在粒子退化等问题，监测结果并不理想。
[0003] 目前针对粒子滤波在非高斯观测噪声环境下应用于接收机自主完好性监测的研究主要集中在改善粒子滤波中的粒子退化问题，目前提出改善方法虽然达到了一定改善粒子退化现象实现多样性的目的，但都在改进过程中导致其对非高斯观测噪声的适应性降低，无法较好的兼顾粒子退化和非高斯观测噪声适应性这两种问题。
[0004] 传统的粒子滤波接收机自主完好性监测方法具有明显的粒子退化问题，导致估计结果的不准确。其他改进的粒子滤波接收机自主完好性监测方法能减少粒子退化效果现象，但无法保证改正后的粒子滤波接收机自主完好性监测方法在非高斯观测噪声环境的效果。
发明内容
[0005] (一)要解决的技术问题
[0006] 本发明要解决的技术问题是：如何实现在非高斯观测噪声环境下有效开展接收机自主完好性监测。(二)技术方案
[0007] 为解决上述技术问题，本发明提供一种基于改进无迹卡尔曼粒子滤波的RAIM系统，所述系统包括：初始化抽样模
块、状态预测模块、
采样模块、重采样模块、估计值获取模块、一致性检验模块；其中，
[0008] 所述初始化抽样模块用于对粒子进行初始化抽样；
[0009] 所述状态预测模块用于进行状态预测；
[0010] 所述采样模块用于对粒子集进行序贯重要性采样；
[0011] 所述重采样模块用于对粒子进行重采样；
[0012] 所述估计值获取模块用于获得状态估计值；
[0013] 所述一致性检验模块用于对滤波信息进行一致性检验。
[0014] 其中，所述初始化抽样模块对粒子进行初始化抽样，首先根据先验概率
密度函数p(x0)
抽取随机样本构成初始粒子集{x0(i)；i＝1,2,…N}～p(x0)，其中N为粒子数目，每个粒子的权重初始化为1/N，表示为{ω0(i)＝1/N；i＝1,2,…N}。
[0015] 其中，所述状态预测模块进行状态预测；通过改进的无迹卡尔曼滤波对粒子状态进行估计，其过程包括时间更新、生成sigma点，观测更新；该状态预测模块包括：时间更新单元、sigma点生成单元、观测更新单元、滤波信息更新单元；其中，
[0016] (21)所述时间更新单元对粒子 进行时间更新；
[0017] 根据公式(3)、(4)及采样得到的粒子集 对粒子进行时间更新；
[0018]
[0019] Pk/k-1(i)＝Φk/k-1Pk-1(i)ΦTk/k-1+Qk (4)
[0020] 其中 是 的最优滤波估计，Φk/k-1为状态转移矩阵下标；k表示当前观测历元，且k＝1,2,3…，k-1表示上一个相邻观测历元，k/k-1表示两个历元之间的过渡，当k＝1时， 为初始化抽样模块中得到的粒子集；Pk/k-1为 的误差协方差矩阵，Qk为过程噪声矩阵；
[0021] (22)所述sigma点生成单元用于生成sigma点；
[0022] 根据公式(5)以及时间更新单元更新后的粒子集 分别生成三类sigma点[0023]
[0024]
[0025] 其中，n为sigma点的总数，Dk|k-1为 的平方根， 为 的协方差矩阵；
[0026] (23)所述观测更新单元用于更新观测信息；
[0027] 基于sigma点生成单元得到的sigma点信息，通过公式(6)、(7)、(8)进行计算得到新的观测信息 误差协方差矩阵 误差矩阵
[0028]
[0029]
[0030]
[0031] 其中， 和 为第j个sigma点的权重，其中 可通过求一阶统计特性的权系数确定， 可通过求二阶统计特性权系数确定， 为 的误差协方差矩阵，Rk为过程噪声的方差矩阵， 为 和 的误差矩阵， 是对每个sigma点的观测更新；
[0032] (24)所述滤波信息更新单元用于更新滤波信息；
[0033] 基于新的观测信息，依据公式(9)、(10)、(11)更新各项滤波信息；
[0034]
[0035]
[0036]
[0037] 其中，Kk为滤波增益矩阵， 为 的误差协方差矩阵。
[0038] 其中，所述所述时间更新单元工作过程中，当k＝1时，使用初始粒子集x0。
[0039] 其中，所述采样模块对粒子集进行序贯重要性采样；所述采样模块包括：信息统计单元、第一采样单元、
[0040] (31)所述信息统计单元通过对粒子集 做信息统计，计算粒子集的均值 和方差
[0041]
[0042]
[0043] (32)所述第一采样单元通过粒子集的均值 和方差 以及建议密度函数采样得到新的粒子集xk(i)：
[0044]
[0045] 其中，所述重采样模块用于对粒子进行重采样；
[0046] 根据粒子权重 计算有效粒子数 然后与
门限值Nth比较，如果Neff
[0047] 其中，所述重采样模块工作过程中的门限值设置为总粒子数的2/3。
[0048] 其中，所述估计值获取模块获得状态估计值；
[0049] 根据重采样的粒子集 和权重信息 计算当前历元的状态估计值：
[0050]
[0051] 其中，所述一致性检验模块对滤波信息进行一致性检验；所述一致性检验模块包括：一致性检验估计值计算单元、对数似然比计算单元、累加对数似然比计算单元、一致性检验模块；
[0052] (61)所述一致性检验估计值计算单元用于计算一致性检验估计值；
[0053] 根据粒子集 和权重信息 计算粒子滤波估计值：
[0054]
[0055]
[0056] 其中，pq(y)为剔除第q个观测分量的粒子滤波估计值，pM(y)为包含所有观测分量的粒子滤波估计值；
[0057] (62)所述对数似然比计算单元用于计算对数似然比；
[0058] 根据粒子滤波估计值(pq(y)和pM(y))计算对数似然比：
[0059]
[0060] 其中，si(y)表示第i颗卫星对应的对数似然比；
[0061] (63)所述累加对数似然比计算单元用于计算累加对数似然比；
[0062] 考虑到一致性检验对结果的敏感程度，还会对上述来两个估计值的对数似然比累加后再进行对比分析，累加对数似然比Sjk(q)的计算过程为：
[0063]
[0064] 其中，yi|Yi-1为第i个历元下Y观测向量中的观测信息；
[0065] (64)所述一致性检验模块用于进行一致性检验；
[0066] 通过对比各观测分量的累加对数似然比估计值，可以实现卫星故障的监测和识别；若系统无故障，则各个分量的累加对数似然比基本一致，且随时间的增加，各曲线整体平稳，没有剧烈
波动的情况；若某个卫星发生故障，则在该时刻起，该卫星对应的 会发生较大的变化，同时与其他卫星对应的估计值具有很大的区别。
[0067] 其中，在所述系统执行过程中，该系统与直接基于无迹卡尔曼滤波的粒子滤波方法的主要区别是该系统在sigma点生成过程中只需一次时间更新，而传统的UKF中需要对每个sigma点进行时间更新，应用到粒子滤波接收机自主完好性监测中，共有n+1个粒子
滤波器，每个粒子滤波中包含N个粒子，并且每个粒子都需要W个sigma点估计，与单纯运用无迹卡尔曼改善粒子滤波接收机自主完好性监测方法相比较，本系统在时间更新步骤中的运算量将是其运算量的1/[(n+1)×N×W]左右。
[0068] (三)有益效果
[0069] 与
现有技术相比较，本发明为实现在非高斯观测噪声环境下有效开展接收机自主完好性监测，通过在粒子滤波方案中引入改进后的无迹卡尔曼滤波，通过改进无迹卡尔曼滤波实现粒子滤波中更合理的状态预测和建议密度计算，能够有效避免粒子退化现象，在保障粒子滤波接收机自主完好性监测系统对非高斯观测噪声环境的适应性的同时，也减少了引入无迹卡尔曼粒子滤波所带来的计算量，保证了
算法运算效率。因此，该方案既改善粒子滤波退化问题，同时保障该方案对非高斯观测噪声的适应性，通过对无迹卡尔曼滤波的改进也减少来由其导致的计算量。
附图说明
[0070] 图1为基于改进无迹卡尔曼粒子滤波接收机自主完好性监测流程示意图。
[0071] 图2为改进无迹卡尔曼粒子滤波流程示意图。
具体实施方式
[0072] 为使本发明的目的、内容、和优点更加清楚，下面结合附图和
实施例，对本发明的具体实施方式作进一步详细描述。
[0073] 为解决上述技术问题，本发明提供一种基于改进无迹卡尔曼粒子滤波的RAIM系统，所述基于改进无迹卡尔曼粒子滤波的接收机自主完好性监测方法通过引入无迹卡尔曼滤波并加以改进，实现更合理的建议密度计算，能够有效避免粒子退化现象，并且保障粒子滤波接收机自主完好性监测方法对非高斯观测噪声环境的适应性，同时与减少了引入无迹卡尔曼粒子滤波所带来的计算量，保证了运算效率，满足接收机完好性监测对粒子滤波方法的需求。
[0074] 所述方法包括两个阶段，如图1所示，其中第一部分为改进无迹卡尔曼粒子滤波，用于实现各卫星状态观测量的估计，此部分主要步骤如图2所示；另一部分为一致性检验，通过多颗卫星粒子滤波轨迹结果的比较分析卫星故障，达到完好性监测的目的。
[0075] 所述系统包括：初始化抽样模块、状态预测模块、采样模块、重采样模块、估计值获取模块、一致性检验模块；其中，
[0076] 所述初始化抽样模块用于对粒子进行初始化抽样；
[0077] 所述状态预测模块用于进行状态预测；
[0078] 所述采样模块用于对粒子集进行序贯重要性采样；
[0079] 所述重采样模块用于对粒子进行重采样；
[0080] 所述估计值获取模块用于获得状态估计值；
[0081] 所述一致性检验模块用于对滤波信息进行一致性检验。
[0082] 其中，所述初始化抽样模块对粒子进行初始化抽样，首先根据先验概率密度函数p(x0)抽取随机样本构成初始粒子集{x0(i)；i＝1,2,…N}～p(x0)，其中N为粒子数目，每个粒子的权重初始化为1/N，表示为{ω0(i)＝1/N；i＝1,2,…N}。
[0083] 其中，所述状态预测模块进行状态预测；通过改进的无迹卡尔曼滤波对粒子状态进行估计，其过程包括时间更新、生成sigma点，观测更新；该状态预测模块包括：时间更新单元、sigma点生成单元、观测更新单元、滤波信息更新单元；其中，
[0084] (21)所述时间更新单元对粒子 进行时间更新；
[0085] 根据公式(3)、(4)及采样得到的粒子集 对粒子进行时间更新；
[0086]
[0087] Pk/k-1(i)＝Φk/k-1Pk-1(i)ΦTk/k-1+Qk (4)
[0088] 其中 是 的最优滤波估计，Φk/k-1为状态转移矩阵下标；k表示当前观测历元，且k＝1,2,3…，k-1表示上一个相邻观测历元，k/k-1表示两个历元之间的过渡，当k＝1时， 为初始化抽样模块中得到的粒子集；Pk/k-1为 的误差协方差矩阵，Qk为过程噪声矩阵；
[0089] (22)所述sigma点生成单元用于生成sigma点；
[0090] 根据公式(5)以及时间更新单元更新后的粒子集 分别生成三类sigma点[0091]
[0092]
[0093] 其中，n为sigma点的总数，Dk|k-1为 的平方根， 为 的协方差矩阵；
[0094] (23)所述观测更新单元用于更新观测信息；
[0095] 基于sigma点生成单元得到的sigma点信息，通过公式(6)、(7)、(8)进行计算得到新的观测信息 误差协方差矩阵 误差矩阵
[0096]
[0097]
[0098]
[0099] 其中， 和 为第j个sigma点的权重，其中 可通过求一阶统计特性的权系数确定， 可通过求二阶统计特性权系数确定， 为 的误差协方差矩阵，Rk为过程噪声的方差矩阵， 为 和 的误差矩阵， 是对每个sigma点的观测更新；
[0100] (24)所述滤波信息更新单元用于更新滤波信息；
[0101] 基于新的观测信息，依据公式(9)、(10)、(11)更新各项滤波信息；
[0102]
[0103]
[0104]
[0105] 其中，Kk为滤波增益矩阵， 为 的误差协方差矩阵。
[0106] 其中，所述所述时间更新单元工作过程中，当k＝1时，使用初始粒子集x0。
[0107] 其中，所述采样模块对粒子集进行序贯重要性采样；所述采样模块包括：信息统计单元、第一采样单元、
[0108] (31)所述信息统计单元通过对粒子集 做信息统计，计算粒子集的均值 和方差[0109]
[0110]
[0111] (32)所述第一采样单元通过粒子集的均值 和方差 以及建议密度函数采样得到新的粒子集xk(i)：
[0112]
[0113] 其中，所述重采样模块用于对粒子进行重采样；
[0114] 根据粒子权重 计算有效粒子数 然后与门限值Nth比较，如果Neff
[0115] 其中，所述重采样模块工作过程中的门限值设置为总粒子数的2/3。
[0116] 其中，所述估计值获取模块获得状态估计值；
[0117] 根据重采样的粒子集 和权重信息 计算当前历元的状态估计值：
[0118]
[0119] 其中，所述一致性检验模块对滤波信息进行一致性检验；所述一致性检验模块包括：一致性检验估计值计算单元、对数似然比计算单元、累加对数似然比计算单元、一致性检验模块；
[0120] (61)所述一致性检验估计值计算单元用于计算一致性检验估计值；
[0121] 根据粒子集 和权重信息 计算粒子滤波估计值：
[0122]
[0123]
[0124] 其中，pq(y)为剔除第q个观测分量的粒子滤波(称为辅助粒子滤波)估计值，pM(y)为包含所有观测分量的粒子滤波(称为主粒子滤波)估计值；
[0125] (62)所述对数似然比计算单元用于计算对数似然比；
[0126] 根据粒子滤波估计值(pq(y)和pM(y))计算对数似然比：
[0127]
[0128] 其中，si(y)表示第i颗卫星对应的对数似然比；
[0129] (63)所述累加对数似然比计算单元用于计算累加对数似然比；
[0130] 考虑到一致性检验对结果的敏感程度，通常还会对上述来两个估计值的对数似然比累加后再进行对比分析，累加对数似然比 的计算过程为：
[0131]
[0132] 其中，yi|Yi-1为第i个历元下Y观测向量中的观测信息；
[0133] (64)所述一致性检验模块用于进行一致性检验；
[0134] 通过对比各观测分量的累加对数似然比估计值，可以实现卫星故障的监测和识别；若系统无故障，则各个分量的累加对数似然比基本一致，且随时间的增加，各曲线整体平稳，没有剧烈波动的情况；若某个卫星发生故障，则在该时刻起，该卫星对应的 会发生较大的变化，同时与其他卫星对应的估计值具有很大的区别。
[0135] 此时粒子滤波接收机自主完好性过程如图2。
[0136] 由于上述过程需要对粒子滤波结果进行一致性检验，则每颗卫星的观测分量都需要一次辅助粒子滤波。假设在k时刻存在M颗可见星，则需要进行M+1次粒子滤波才能完成完好一次完好性监测，这要求应用在完好性监测中的粒子滤波方法需具有较高的计算效率。
[0137] 其中，在所述系统执行过程中，该系统与直接基于无迹卡尔曼滤波的粒子滤波方法的主要区别是该系统在sigma点生成过程中只需一次时间更新，而传统的UKF中需要对每个sigma点进行时间更新，应用到粒子滤波接收机自主完好性监测中，共有n+1个粒子滤波器，每个粒子滤波中包含N个粒子，并且每个粒子都需要W个sigma点估计，与单纯运用无迹卡尔曼改善粒子滤波接收机自主完好性监测方法相比较，本系统在时间更新步骤中的运算量将是其运算量的1/[(n+1)×N×W]左右。
[0138] 此外，本发明还提供一种基于改进无迹卡尔曼粒子滤波的RAIM方法，所述基于改进无迹卡尔曼粒子滤波的接收机自主完好性监测方法通过引入无迹卡尔曼滤波并加以改进，实现更合理的建议密度计算，能够有效避免粒子退化现象，并且保障粒子滤波接收机自主完好性监测方法对非高斯观测噪声环境的适应性，同时与减少了引入无迹卡尔曼粒子滤波所带来的计算量，保证了运算效率，满足接收机完好性监测对粒子滤波方法的需求。
[0139] 所述方法包括两个阶段，如图1所示，其中第一部分为改进无迹卡尔曼粒子滤波，用于实现各卫星状态观测量的估计，此部分主要步骤如图2所示；另一部分为一致性检验，通过多颗卫星粒子滤波轨迹结果的比较分析卫星故障，达到完好性监测的目的。
[0140] 所述方法包括如下步骤：
[0141] 步骤1：对粒子进行初始化抽样；
[0142] 步骤2：进行状态预测；
[0143] 步骤3：对粒子集进行序贯重要性采样；
[0144] 步骤4：对粒子进行重采样；
[0145] 步骤5：获得状态估计值；
[0146] 步骤6：对滤波信息进行一致性检验。
[0147] 其中，所述步骤1中对粒子进行初始化抽样，首先根据先验概率密度函数p(x0)抽取随机样本构成初始粒子集{x0(i)；i＝1,2,…N}～p(x0)，其中N为粒子数目，每个粒子的权重初始化为1/N，表示为{ω0(i)＝1/N；i＝1,2,…N}。
[0148] 其中，所述步骤2进行状态预测；
[0149] 通过改进的无迹卡尔曼滤波对粒子状态进行估计，其过程包括时间更新、生成sigma点，观测更新；该步骤2包括：
[0150] 步骤21：对粒子 进行时间更新；
[0151] 根据公式(3)、(4)及采样得到的粒子集 对粒子进行时间更新；
[0152]
[0153] Pk/k-1(i)＝Φk/k-1Pk-1(i)ΦTk/k-1+Qk (4)
[0154] 其中 是 的最优滤波估计，Φk/k-1为状态转移矩阵下标；k表示当前观测历元，且k＝1,2,3…，k-1表示上一个相邻观测历元，k/k-1表示两个历元之间的过渡，当k＝1时， 为步骤一中得到的粒子集；Pk/k-1为 的误差协方差矩阵，Qk为过程噪声矩阵；
[0155] 步骤22：生成sigma点；
[0156] 根据公式(5)以及上一步更新后的粒子集 分别生成三类sigma点
[0157]
[0158] 其中，n为sigma点的总数，Dk|k-1为 的平方根， 为 的协方差矩阵；
[0159] 步骤23：更新观测信息；
[0160] 基于步骤22中得到sigma点信息，通过公式(6)、(7)、(8)进行计算得到新的观测信息 误差协方差矩阵 误差矩阵
[0161]
[0162]
[0163]
[0164] 其中， 和 为第j个sigma点的权重，其中 可通过求一阶统计特性的权系数确定， 可通过求二阶统计特性权系数确定， 为 的误差协方差矩阵，Rk为过程噪声的方差矩阵， 为 和 的误差矩阵， 是对每个sigma点的观测更新；
[0165] 步骤24：更新滤波信息；
[0166] 基于新的观测信息，依据公式(9)、(10)、(11)更新各项滤波信息；
[0167]
[0168]
[0169]
[0170] 其中，Kk为滤波增益矩阵， 为 的误差协方差矩阵。
[0171] 其中，所述步骤21中，当k＝1时，使用初始粒子集x0。
[0172] 其中，所述步骤3对粒子集进行序贯重要性采样；该步骤3包括：
[0173] 步骤31：通过对粒子集 做信息统计，计算粒子集的均值 和方差[0174]
[0175]
[0176] 步骤32：通过粒子集的均值 和方差 以及建议密度函数采样得到新的粒子集xk(i)：
[0177]
[0178] 其中，所述步骤4对粒子进行重采样；
[0179] 根据粒子权重 计算有效粒子数 然后与门限值Nth比较，如果Neff
[0180] 其中，所述步骤4中的门限值设置为总粒子数的2/3。
[0181] 其中，所述步骤5获得状态估计值；
[0182] 根据重采样的粒子集 和权重信息 计算当前历元的状态估计值：
[0183]
[0184] 其中，所述步骤6对滤波信息进行一致性检验；所述步骤6包括：
[0185] 步骤61：计算一致性检验估计值；
[0186] 根据粒子集 和权重信息 计算粒子滤波估计值：
[0187]
[0188]
[0189] 其中，pq(y)为剔除第q个观测分量的粒子滤波(称为辅助粒子滤波)估计值，pM(y)为包含所有观测分量的粒子滤波(称为主粒子滤波)估计值；
[0190] 步骤62：计算对数似然比；
[0191] 根据粒子滤波估计值(pq(y)和pM(y))计算对数似然比：
[0192]
[0193] 其中，si(y)表示第i颗卫星对应的对数似然比；
[0194] 步骤63：计算累加对数似然比；
[0195] 考虑到一致性检验对结果的敏感程度，通常还会对上述来两个估计值的对数似然比累加后再进行对比分析，累加对数似然比 的计算过程为：
[0196]
[0197] 其中，yi|Yi-1为第i个历元下Y观测向量中的观测信息；
[0198] 步骤64：一致性检验；
[0199] 通过对比各观测分量的累加对数似然比估计值，可以实现卫星故障的监测和识别；若系统无故障，则各个分量的累加对数似然比基本一致，且随时间的增加，各曲线整体平稳，没有剧烈波动的情况；若某个卫星发生故障，则在该时刻起，该卫星对应的 会发生较大的变化，同时与其他卫星对应的估计值具有很大的区别。
[0200] 此时粒子滤波接收机自主完好性过程如图2。
[0201] 由于上述过程需要对粒子滤波结果进行一致性检验，则每颗卫星的观测分量都需要一次辅助粒子滤波。假设在k时刻存在M颗可见星，则需要进行M+1次粒子滤波才能完成完好一次完好性监测，这要求应用在完好性监测中的粒子滤波方法需具有较高的计算效率。
[0202] 在所述方法执行过程中，该方法与直接基于无迹卡尔曼滤波的粒子滤波方法的主要区别是该方法在sigma点生成过程中只需一次时间更新，而传统的UKF中需要对每个sigma点进行时间更新，应用到粒子滤波接收机自主完好性监测中，共有n+1个粒子滤波器，每个粒子滤波中包含N个粒子，并且每个粒子都需要W个sigma点估计，与单纯运用无迹卡尔曼改善粒子滤波接收机自主完好性监测方法相比较，本方法在时间更新步骤中的运算量将是其运算量的1/[(n+1)×N×W]左右。
[0203] 综上，本发明属于全球卫星导航检测技术领域，具体涉及一种基于改进无迹卡尔曼粒子滤波的RAIM方法。所述方法涉及
定位解算、观测
数据处理、完好性监测、无迹卡尔曼滤波、粒子滤波等相关技术，主要针对在复杂观测噪声环境下，卫星导航定位服务过程中可能出现的观测数据故障等问题，通过采用改进的无机卡尔曼滤波与粒子滤波进行组合，改进传统粒子滤波粒子退化问题，同时兼顾非高斯观测噪声环境和高效率的应用场景，给出一种有效的全球卫星
导航系统接收机自主完好性监测监测方法。
[0204] 以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和
变形，这些改进和变形也应视为本发明的保护范围。
