CN112902967A - 一种基于残差卡方-改进序贯概率比的抗欺骗导航方法 - Google Patents

Abstract

本发明公开了一种基于残差卡方‑改进序贯概率比的抗欺骗导航方法。该方法步骤如下：采集捷联惯性导航系统、卫星导航系统、里程计、高度计、磁罗盘的导航信息；以捷联惯性导航系统模型为信息融合模型的基准框架，与其余导航源构成子滤波器，进行信息融合；对所有的子滤波器的输出利用改进的残差χ²检验法进行故障诊断，并针对SINS/GNSS子系统的输出结果利用残差χ²‑改进SPRT的方法进行故障诊断，使得对于该子系统可能发生的欺骗干扰能够快速的识别并隔离；故障诊断与隔离完成后，进行主滤波器的信息融合与分配。本发明提高了多源组合导航系统在干扰环境下的导航精度，特别是增强了系统对于欺骗干扰的容错能力。

Description

一种基于残差卡方-改进序贯概率比的抗欺骗导航方法

技术领域

本发明涉及组合导航技术领域，特别是一种基于残差卡方-改进序贯概率比的抗欺骗导航方法。

背景技术

卫星导航GNSS在军事应用和经济发展中得到日益广泛的应用，重要性日趋明显。目前随着GPS、GLONASS和北斗等卫星导航系统的蓬勃发展，它不仅在民用上产生了巨大的经济价值，同时也在现代高科技战争中为信息化作战系统和精确制导武器提供精确的时空基准，发挥着越来越重要的作用。

然而随着卫星导航重要性的提升，针对卫导的欺骗干扰技术也应运而生，并在战场上发挥了巨大作用。2011年伊朗防空部队在该国东部边境利用欺骗干扰技术俘获的一架美国“RQ-170”无人侦察机。2012年美国国防部组织的对无人直升机GPS接收终端进行欺骗干扰实验，利用低成本的欺骗干扰设备释放欺骗信号，使得本应控制保持固定高度的无人直升机自动驾驶系统不断控制直升机下降。以上案例说明欺骗干扰不只是理论层面的安全隐患，已经成为了卫星导航系统面临的主要现实威胁之一。可见，通过欺骗干扰手段，攻击方可在消耗敌方武装力量的同时，完整保留敌方装备状态，用于提升自身的军事技术水平.这一危害某种意义上远大于以往各种电磁干扰手段。因此当前导航战背景下的卫星导航抗欺骗干扰技术的研究和发展亟需解决。

发明内容

本发明的目的在于提供一种基于残差卡方-改进序贯概率比的抗欺骗导航方法，基于联邦卡尔曼滤波的技术，针对转发式欺骗干扰产生的突变故障和生成式欺骗干扰产生的缓变故障，采用残差χ²-改进SPRT联合故障检测方法进行故障诊断和隔离重构，提高组合导航系统在故障及干扰环境下的容错能力和定位精度。

实现本发明目的的技术解决方案为：一种基于残差卡方-改进序贯概率比的抗欺骗导航方法，包括以下步骤：

步骤1：建立在导航系下的联邦滤波组合导航系统模型，设计联邦滤波器，经过联邦滤波器的最优估计，对导航系统进行校正输出；

步骤2：通过改进的残差χ²检验法对联邦滤波器中的子滤波器进行故障诊断和隔离，在残差χ²计算得到残差向量特征值的基础上对检测阈值进行模糊逻辑和加权处理；

步骤3：针对子滤波器中对导航结果起关键作用的SINS/GNSS子系统，采用改进的SPRT方法进行检测，并将结果和步骤2中的检测结果联合进行故障判决，判断GNSS是否受到欺骗干扰。

本发明与现有技术相比，其显著优点是：(1)在常规的GNSS/SINS组合导航算法的基础上，设计基于联邦卡尔曼滤波算法，以SINS为参考系统，GNSS/SINS、里程计/SINS、高度计/SINS、磁罗盘/SINS四者组合的多源导航系统，提高了组合导航系统的导航定位精度；(2)设计了一种基于改进残差χ²检测的故障检测与故障隔离的方案，当子系统发生故障时能够实时地对子系统进行隔离，使系统具有较好的稳定性和可靠性；(3)设计了一种基于残差χ²-改进SPRT联合故障检测方法，对SINS/GNSS导航子系统进行检测，提高系统对于卫星导航欺骗干扰的及时识别并隔离的能力。

附图说明

图1是本发明的多源组合导航系统结构流程图。

图2是联邦卡尔曼滤波的结构示意图。

图3是残差χ²-改进SPRT联合故障检测模块的结构示意图。

具体实施方式

本发明一种基于残差卡方-改进序贯概率比的抗欺骗导航方法，包括以下步骤：

步骤1：建立在导航系下的联邦滤波组合导航系统模型，设计联邦滤波器，经过联邦滤波器的最优估计，对导航系统进行校正输出；

步骤2：通过改进的残差χ²检验法对联邦滤波器中的子滤波器进行故障诊断和隔离，在残差χ²计算得到残差向量特征值的基础上对检测阈值进行模糊逻辑和加权处理；

步骤3：针对子滤波器中对导航结果起关键作用的SINS/GNSS子系统，采用改进的SPRT方法进行检测，并将结果和步骤2中的检测结果联合进行故障判决，判断GNSS是否受到欺骗干扰。

进一步地，步骤1中所述建立在导航系下的联邦滤波组合导航系统模型，设计联邦滤波器，经过联邦滤波器的最优估计，对导航系统进行校正输出，具体如下：

(1.1)首先在导航坐标系下建立联邦滤波组合导航系统模型；惯性导航系统的误差模型选用通用误差模型，公共误差参考系统选取惯性导航系统，导航坐标系选取东北天地理坐标系，则系统的状态矢量X为

其中

为东北天方向上的姿态角误差；δV_E、δV_N、δV_U为东北天方向上的速度误差；δL、δλ、δH为纬经高方向的位置误差；ε_bx、ε_by、ε_bz为陀螺仪在三个轴向上的随机漂移；ε_rx、ε_ry、ε_rz为陀螺仪在三个轴向上的一阶马尔可夫过程随机噪声；

为加速度计在三个轴向上的常值偏差；

系统的状态方程

为：

其中F(t)为系统的状态转移矩阵，G(t)为误差系数矩阵，W(t)为白噪声随机误差矢量，具体形式如下：

状态转移矩阵F(t)：

式中F_N为系统误差矩阵，F_s为惯性器件的误差转换矩阵，F_M为惯性器件的噪声矩阵，具体如下：

其中I_3×3为3维的单位矩阵，

为姿态转换矩阵；

T_gx、T_gy、T_gz分别为陀螺仪三个轴的相关时间，T_ax、T_ay、T_az分别为加速度计三个轴的相关时间；

误差系数矩阵G(t)为：

白噪声随机误差矢量W(t)为：

W(t)＝[w_gx w_gy w_gz w_rx w_ry w_rz w_ax w_ay w_az]^T (8)

w_gx、w_gy、w_gz为陀螺仪的白噪声误差，w_rx、w_ry、w_rz为陀螺仪扩展误差模型的白噪声误差，w_ax、w_ay、w_az为加速度计的白噪声误差；

系统的量测方程Z为：

Z(t)＝H(t)X(t)+V(t) (9)

其中X(t)为系统的状态矢量，V(t)＝[N_E N_N N_U]^T，N_E、N_N、N_U分别为GPS接收机沿东、北、天方向的位置误差，H(t)根据不同的子滤波器选取；

将系统模型离散化，把状态方程式(2)和量测方程式(9)离散化，得

X_k为k时刻的状态，X_k-1为k-1的状态，W_k-1为离散的白噪声随机误差矢量，Z_k为k时刻的量测值，H_k为离散的量测矩阵，V_k为离散的量测噪声；

式中状态转移矩阵Φ_k|k-1和过程噪声矩阵Γ_k-1定义如下

式中，T为迭代周期，n为大于0的正整数；

(1.2)设计联邦子滤波器

每个子滤波器均采用卡尔曼滤波方式，具体如下：

状态一步预测方程：

式中，

为一步预测的状态估计值，

为k-1时刻的状态估计值；状态估计方程：

式中，

为k时刻的状态估计值，K_k为k时刻的卡尔曼增益；

卡尔曼增益方程：

式中，P_k|k-1为一步预测均方差，R_k为量测噪声的方差矩阵；

一步预测均方差方程：

式中，Q_k-1为k-1时刻的系统噪声的方差矩阵；

均方差估计方程：

式中，P_k|k为状态估计的均方误差；

表示惯性导航系统的位置信息为

表示惯性导航系统的速度信息为

式中λ_I、L_I、H_I为惯性导航系统纬度、经度、高度的测量值，λ_t、L_t、H_t为载体纬度、经度、高度的真值，δλ、δL、δH为惯性导航系统经度、纬度、高度的测量误差值；v_IE、v_IN、v_IU为惯性导航系统在东北天方向上速度的测量值，v_E、v_N、v_U为载体在东北天方向上速度的真值，δv_E、δv_N、δv_U为惯性导航系统在东北天方向上速度的测量误差值；

(1.2.1)设计GNSS/SINS子滤波器

表示卫星导航系统的位置信息为

表示卫星导航系统的速度信息为

式中λ_G、L_G、H_G为卫星导航系统经度、纬度、高度的测量值，λ_t、L_t、H_t为载体经度、纬度、高度的真值，N_E、N_N、N_U为卫星导航系统在东北天方向上位置的测量误差值；v_GE、v_GN、v_GU为卫星导航系统在东北天方向上速度的测量值，v_E、v_N、v_U为载体在东北天方向上速度的真值，M_E、M_N、M_U为惯性导航系统在东北天方向上速度的测量误差值，R_M为子午圈半径，R_N为卯酉圈半径；

根据(1.1)知，GNSS/SINS子滤波器的状态方程为

GNSS/SINS子滤波器的量测方程为

其中，GNSS/SINS子滤波器的量测矩阵H_G和量测噪声矩阵V_G定义如下

V_G(t)＝[N_N N_E N_U M_E M_N M_U]^T (25)

量测噪声当作均值为零的白噪声处理，方差分别为

(1.2.2)设计里程计/SINS子滤波器

表示里程计的速度信息为

式中v_OE、v_ON、v_OU为里程计在东北天方向上速度的测量值，v_E、v_N、v_U为载体在东北天方向上速度的真值，O_E、O_N、O_U为里程计在东北天方向上速度的测量误差值；

由(1.1)得，里程计/SINS子滤波器的状态方程为

里程计/SINS子滤波器的量测方程为

其中，里程计/SINS子滤波器的量测矩阵H_O(t)和量测噪声矩阵V_O(t)定义为

H_O(t)＝[0_3×3 diag[1 1 1]0_3×12] (29)

V_O(t)＝[O_E O_N O_U]^T (30)

量测噪声作为均值为零的白噪声处理；

(1.2.3)设计高度计/SINS子滤波器

表示高度计的高度信息为

h_a＝h_t-A (31)

式中h_a为高度计在高度上的测量值，h_t为载体在天向上位置的真值，A为高度计在高度上的测量误差值；

根据(1.1)知，高度计/SINS子滤波器的状态方程为

高度计/SINS子滤波器的量测方程为

Z_a(t)＝h_I-h_a＝H_a(t)X(t)+V_a(t) (33)

其中高度计/SINS子滤波器的量测矩阵H_a(t)和量测噪声矩阵V_a(t)定义为

H_a(t)＝[0_3×6 diag[0 0 1]0_3×9] (34)

V_a(t)＝A (35)

量测噪声作为均值为零的白噪声处理；

(1.2.4)设计磁罗盘/SINS子滤波器

表示磁罗盘的角度信息为

ψ_c＝ψ_t-δψ (36)

式中ψ_c为磁罗盘测量的磁航向经过磁偏补偿后得到的航向角，ψ_t为载体航向角的真值，δψ为磁罗盘航向角的量测误差；

根据(1.1)知，磁罗盘/SINS子滤波器的状态方程为

磁罗盘/SINS子滤波器的量测方程为

Z_c(t)＝ψ_I-ψ_c＝H_c(t)X(t)+V_c(t) (38)

其中磁罗盘/SINS子滤波器的量测矩阵H_c(t)和量测噪声矩阵V_c(t)定义为

H_c＝[-sinψtanθ -cosψtanθ 1 0_1×12] (39)

V_c(t)＝δψ (40)

式中，ψ为偏航角，θ为俯仰角，量测噪声作为均值为零的白噪声处理；

(1.3)设计联邦滤波主滤波器

联邦滤波器是一种两级滤波结构，公共参考子系统SINS的输出X_k一方面给主滤波器，另一方面给各子滤波器作为量测值；各传感器子系统的输出只给各自的滤波器，各子滤波器的局部估计值X_i及其协方差阵P_i送至主滤波器与主滤波器的估计值一起进行融合，从而得到全局最优估计值；

由主滤波器和子滤波器合成的全局最优估计值

及其相应的协方差阵P_g被放大为

后再反馈给各子滤波器，来重置子滤波器的估计值，即：

主滤波器的协方差阵为

β_m和β_i分别是主滤波器与各子滤波器的信息分配系数，取值根据信息分配原则确定；

采用β_m＝0,β_i＝1/N，无重置结构的联邦卡尔曼滤波，各个局部滤波器独立完成滤波，仅对子滤波器估计值进行全局融合，并将全局最优估计值和协方差矩阵反馈到子滤波器中；

系统信息在主、子滤波器间的分配方法是基于信息分配原则的，且满足：

其中，β_m为主滤波器信息分配系数，β_i为第i个子滤波器对应的信息分配系数；

系统信息在子滤波器中的具体分配方法如下：

其中，

为主滤波器状态估计，P_g、Q_g分别为主滤波器状态协方差阵和系统噪声阵；

为子滤波器状态估计，P_i、Q_i分别为子滤波器状态估计对应的状态估计均方差阵和子滤波器系统噪声；

主滤波器最优信息融合：

式中

为第i个子滤波器的协方差矩阵；

为主滤波器中的协方差矩阵；

为主滤波器和子滤波器融合后的协方差矩阵；

为融合后的状态估计；

为第i个子滤波器的状态估计值；

为主滤波器的状态估计值；N为子滤波器的个数。

进一步地，步骤2所述通过改进的残差χ²检验法对联邦滤波器中的子滤波器进行故障诊断和隔离，在残差χ²计算得到残差向量特征值的基础上对检测阈值进行模糊逻辑和加权处理；具体如下：

(2.1)残差χ²检验

用残差χ²检验对系统的故障进行检测和隔离；由于每个子滤波器均为卡尔曼滤波器，新息向量也即预测残差为：

其中一步预测值

当无故障发生时卡尔曼滤波器的残差r_k是零均值的白噪声，方差为

当系统发生故障时，残差的均值就不再为0，因此对残差的检验能够确定系统是否发生了故障；

对r_k做二元假设：

H₀无故障时：

式中，r_k为子滤波器的预测残差，A_k为残差方差；

H₁有故障时：

式中，μ为常数向量；

定义故障检测函数为：

式中λ_k为故障检测值，是服从自由度为m的χ²分布，即λ_k～χ²(m)，m为量测矩阵Z_k的维数；

故障判定准则为：

其中设置的阈值T_D由误警率P_f确定，给定误警率P_f＝α，则由如下计算公式求出T_D：

式中，λ为均值，n为方差；

(2.2)改进的残差χ²检验法

采用改进的残差χ²检验法，在算出残差向量特征值的基础上对检测阈值进行模糊逻辑和加权处理，能够提高系统对缓变故障的检测能力；

对残差向量进行归一化处理：

式中r_k和A_k与(2.1)中定义相同，

为k时刻归一化后的残差向量；

从当前时刻开始往前选取n个时刻的残差向量组成残差矩阵R：

定义残差的方差矩阵Q为：

Q＝R^TR (52)

根据控制理论知，Q矩阵特征值的最大值有界，即：

λ_max≤S (53)

通过选取最近的N个矩阵最大特征值的平均值进一步提高数据的可靠性，具体的表达式如下：

式中，

为N个矩阵最大特征值的平均值，λ_max(k)为k时刻Q矩阵的最大特征值；

考虑实时性和计算量的问题，建议N与n的值不要取得过大；

此时，故障判断准则为：

(2.3)当某个子滤波器检测出故障时，系统舍弃该子滤波器的信息，主滤波器对剩余的子滤波器进行重构和信息融合，使得剩余的导航子系统能够继续提供导航信息，保证了多源组合导航系统具有一定的容错性能。

进一步地，步骤3中所述针对子滤波器中对导航结果起关键作用的SINS/GNSS子系统，采用改进的SPRT方法进行检测，并将结果和步骤2中的检测结果联合进行故障判决，判断GNSS是否受到欺骗干扰，具体如下：

(3.1)残差χ²检测

利用(2.1)中阐述的残差χ²检测法对SINS/GNSS子滤波器的估计结果进行故障诊断；

(3.2)改进SPRT检测

假设某未知正态随机变量x的k次序贯独立样本为{x_i|i＝1,2,…,k}，根据概率统计理论，近似有x～

且

式中x_i为第i次样本值，

为样本均值，既是有故障时的均值，也为无故障时的均值；

为样本方差；

假定状态变量x的量测实际值为x^*，则定义如下两个事件

其中，x₀为故障状态下的真实值；

量测序列x₁,x₂,…,x_k属于两个样本类H₀-正常类和H₁-故障类之一，概率密度函数为

式中，p(x_i|H₀)和p(x_i|H₁)分别为两个样本类的概率密度；

则得似然比L(k)：

对式(59)取自然对数，得到对数似然比λ(k)：

将式(60)简化，得：

由式(61)知，λ(k)不可能为负值，当系统处于正常工作状态时，随着k的增加，样本均值

将逐渐逼近x₀，因此λ(k)→0；当系统发生故障时，随着k的增加，样本均值

依然会逐渐逼近故障情况下的真值，但真值不为x₀；

综上，得到似然比λ(k)的递推计算公式：

式中，

为k-1个样本的均值，

为此时的样本方差；

若变量x的先验方差为σ²，则用σ²替代上式中的

若不确定变量x的先验方差，则用实时计算的样本方差替代；

Wald序贯概率算法中给出了两个检测阈值，分别为T(H₀)和T(H₁)；对于上述改进的SPRT算法，则仅需设定一个检测阈值T(H₁)，该值的计算方式与Wald算法一致，即由误警率P_f和漏检率P_m确定：

则判定准则为

(3.3)故障判决

上述(3.1)和(3.2)两种故障检测方法是并行计算的；对SINS/GNSS子滤波器的残差值分别使用残差χ²和改进SPRT两种检测法进行故障检测，得到两个检测结果，并且将当前时刻的改进SPRT检测值与上一时刻的作差得到增量；把增量和两个检测结果一同输入到系统故障综合决策规则中，进行最终故障判决，从而得出最终的检测结果；系统故障综合决策规则见表1，表中用“是”与“否”表示是否发生故障；

表1系统故障综合决策规则

由此得到对SINS/GNSS子滤波器估计结果进行故障检测的结果。

下面结合附图及具体实施例对本发明作进一步详细说明。

实施例

结合图1，本实施例为一种基于残差χ²-改进SPRT联合故障检测的抗欺骗导航方法，包括以下步骤：

步骤1：采集来自SINS、GNSS、里程计、高度计、磁罗盘等导航信息源经过如图1中解算并输出的导航信息。

步骤2：如图1和图2中，以SINS为公共信息参考系统，与其他传感器系统形成各个子滤波器，分别进行状态估计。

(2.1)建立子系统的状态方程

以惯性导航系统的误差模型为通用误差模型，公共误差参考系统选取惯性导航系统，导航坐标系选取东北天地理坐标系，则系统的状态矢量X为

其中

为东北天方向上的姿态角误差；δV_E、δV_N、δV_U为东北天方向上的速度误差；δL、δλ、δH为纬经高方向的位置误差；ε_bx、ε_by、ε_bz为陀螺仪在三个轴向上的随机漂移；ε_rx、ε_ry、ε_rz为陀螺仪在三个轴向上的一阶马尔可夫过程随机噪声；

为加速度计在三个轴向上的常值偏差。

系统的状态方程

为：

其中F(t)为系统的状态转移矩阵，G(t)为误差系数矩阵，W(t)为白噪声随机误差矢量，具体形式如下：

状态转移矩阵F(t)：

式中F_N为系统误差矩阵，F_s为惯性器件的误差转换矩阵，F_M为惯性器件的噪声矩阵，具体如下：

其中

为姿态转换矩阵；

误差系数矩阵G(t)为：

白噪声随机误差矢量W(t)为：

W(t)＝[w_gx w_gy w_gz w_rx w_ry w_rz w_ax w_ay w_az]^T (1.8)

系统的量测方程Z为：

Z(t)＝H(t)X(t)+V(t) (1.9)

其中X(t)为系统的状态矢量，V(t)＝[N_N N_E N_U]^T，H(t)根据不同的子滤波器选取。

将系统模型离散化，把状态方程式(1.2)和量测方程式(1.9)离散化，得

式中

式中，T为迭代周期。

(2.2)建立各子系统的量测方程

(2.2.1)GNSS/SINS子系统的量测方程

根据(2.1)知，GNSS/SINS子系统的状态方程为

GNSS/SINS子系统的量测方程为

其中

V_G(t)＝[N_N N_E N_U M_E M_N M_U]^T (1.16)

式中λ_G、L_G、H_G为卫星导航系统经度、纬度、高度的测量值，λ_t、L_t、H_t为载体经度、纬度、高度的真值，N_E、N_N、N_U为卫星导航系统在东北天方向上位置的测量误差值；v_GE、v_GN、v_GU为卫星导航系统在东北天方向上速度的测量值，v_E、v_N、v_U为载体在东北天方向上速度的真值，M_E、M_N、M_U为惯性导航系统在东北天方向上速度的测量误差值。

量测噪声作为均值为零的白噪声处理，其方差分别为

(2.2.2)里程计/SINS子系统的量测方程

根据(2.1)知，里程计/SINS子系统的状态方程为

里程计/SINS子系统的量测方程为

其中

H_O(t)＝[0_3×3 diag[1 1 1]0_3×12] (1.19)

V_O(t)＝[O_E O_N O_U]^T (1.20)

式中v_OE、v_ON、v_OU为里程计在东北天方向上速度的测量值，v_E、v_N、v_U为载体在东北天方向上速度的真值，O_E、O_N、O_U为里程计在东北天方向上速度的测量误差值，量测噪声作为均值为零的白噪声处理。

(2.2.3)高度计/SINS子系统的量测方程

根据(2.1)知，高度计/SINS子系统的状态方程为

高度计/SINS子系统的量测方程为

Z_a(t)＝h_I-h_a＝H_a(t)X(t)+V_a(t)(1.22)

其中

H_a(t)＝[0_3×6 diag[0 0 1]0_3×9] (1.23)

V_a(t)＝A (1.24)

式中h_a为高度计在高度上的测量值，h_t为载体在天向上位置的真值，A为高度计在高度上的测量误差值，量测噪声作为均值为零的白噪声处理。

(2.2.4)磁罗盘/SINS子系统的量测方程

根据(2.1)知，磁罗盘/SINS子系统的状态方程为

磁罗盘/SINS子系统的量测方程为

Z_c(t)＝ψ_I-ψ_c＝H_c(t)X(t)+V_c(t) (1.26)

其中

H_c＝[-sinψtanθ -cosψtanθ 1 0_1×12] (1.27)

V_c(t)＝δψ (1.28)

式中ψ_c为磁罗盘测量的磁航向经过磁偏补偿后得到的航向角，ψ_t为载体航向角的真值，δψ为磁罗盘航向角的量测误差，量测噪声作为均值为零的白噪声处理。

(2.3)设计联邦子滤波器：

每个子滤波器均采用卡尔曼滤波方式，具体如下：

状态一步预测方程：

状态估计方程：

卡尔曼增益方程：

一步预测均方差方程：

均方差估计方程：

步骤3：故障检测模块

(3.1)对除SINS/GNSS子系统外的其他子系统进行残差χ²检验(故障检测1)

用残差χ²检验对系统的故障进行检测和隔离；由于每个子滤波器均为卡尔曼滤波器，其新息向量(也称为预测残差)为：

其中一步预测值

当无故障发生时，卡尔曼滤波器的残差r_k是零均值的白噪声，其方差为

当系统发生故障时，残差的均值就不再为0，因此对残差进行检验可以确定系统是否发生了故障。

对r_k做二元假设，H₀无故障时：

H₁有故障时：

定义故障检测函数为：

式中λ_k是服从自由度为m的χ²分布，即λ_k～χ²(m)，m为量测矩阵Z_k的维数。

故障判定准则为：

其中设置的阈值T_D由误警率P_f确定。给定误警率P_f＝α，则由如下计算公式求出T_D：

(3.2)对SINS/GNSS子系统进行残差χ²-改进SPRT联合故障检测(故障检测2)

(3.2.1)残差χ²检测

利用(3.1)中阐述的残差χ²检测法对SINS/GNSS子滤波器的估计结果进行故障诊断。

(3.2.2)改进SPRT检测

假设某未知正态随机变量x的k次序贯独立样本为{x_i|i＝1,2,…,k}，根据概率统计理论，近似有x～

且

式中

为样本均值，既是有故障时的均值，也为无故障时的均值，

为样本方差。

假定状态变量x的量测实际值为x^*，则可以定义如下两个事件

其中，x₀为故障状态下的真实值。

量测序列x₁,x₂,…,x_k属于两个样本类H₀(正常类)和H₁(故障类)之一，其概率密度函数为:

则可得似然比：

对式(1.43)取自然对数，得到对数似然比：

将式(1.44)简化，可得：

由式(1.45)可知，λ(k)不可能为负值，当系统处于正常工作状态时，随着k的增加，样本均值

将逐渐逼近x₀，因此λ(k)→0；当系统发生故障时，随着k的增加，样本均值

依然会逐渐逼近故障情况下的真值，但真值不为x₀。

综上，得到似然比λ(k)的递推计算公式：

若变量x的先验方差为σ²，则用其替代上式中的

若不确定变量x的先验方差，则可用实时计算的样本方差替代。

对于上述改进的SPRT算法，需设定一个检测阈值T(H₁)，该值的计算方式与Wald算法一致，即由误警率P_f和漏检率P_m确定：

则判定准则为:

(3.2.3)故障判决

上述(3.2.1)和(3.2.2)两种故障检测方法是并行计算的，如图3所示。

对SINS/GNSS子滤波器的残差值分别使用残差χ²和改进SPRT两种检测法进行故障检测，得到两个检测结果，并且将当前时刻的改进SPRT检测值与上一时刻的作差得到增量。把增量和两个检测结果一同输入到系统故障综合决策规则中，进行最终故障判决，从而得出最终的检测结果。系统故障综合决策规则见下表，表中用“是”与“否”表示是否发生故障。

表1系统故障综合决策规则

由此可以得到对SINS/GNSS子滤波器的估计值进行故障检测的结果。

步骤4：当某个子滤波器检测出故障时，系统舍弃该子滤波器的信息。主滤波器对剩余的子滤波器进行重构和信息融合。

如图2所示，公共信息参考系统SINS的输出X_k一方面给主滤波器，另一方面给各子滤波器作为量测值；各传感器子系统的输出只给各自的滤波器，各子滤波器的局部估计值X_i及其协方差阵P_i送至主滤波器与主滤波器的估计值一起进行融合，从而得到全局最优估计值。

由主滤波器和子滤波器合成的全局最优估计值

及其相应的协方差阵P_g被放大为

后再反馈给各子滤波器，来重置子滤波器的估计值，即：

主滤波器的协方差阵为

β_i根据信息分配原则确定。

采用β_m＝0,β_i＝1/N，无重置结构的联邦卡尔曼滤波，各个局部滤波器独立完成滤波，仅对子滤波器估计值进行全局融合，并将全局最优估计值和协方差矩阵反馈到子滤波器中。

系统信息在主、子滤波器间的分配方法是基于信息分配原则的，且满足：

其中，β_m为主滤波器信息分配系数，β_i为第i个子滤波器对应的信息分配系数。

系统信息在子滤波器中的具体分配方法如下：

其中，

为主滤波器状态估计，P_g、Q_g分别为主滤波器状态协方差阵和系统噪声阵；

为子滤波器状态估计，P_i、Q_i分别为子滤波器状态估计对应的状态估计均方程阵和子滤波器系统噪声。

主滤波器最优信息融合：

式中

为主滤波器中的协方差矩阵；

为主滤波器和子滤波器融合后的协方差矩阵；

为融合后的状态估计。

综上所述，本发明在常规的GNSS/SINS组合导航算法的基础上，设计基于联邦卡尔曼滤波算法，以SINS为参考系统，GNSS/SINS、里程计/SINS、高度计/SINS、磁罗盘/SINS四者组合的多源导航系统，提高了组合导航系统的导航定位精度；设计了一种基于改进残差χ²检测的故障检测与故障隔离的方案，当子系统发生故障时能够实时地对子系统进行隔离，使系统具有较好的稳定性和可靠性；设计了一种基于残差χ²-改进SPRT联合故障检测方法，对SINS/GNSS导航子系统进行检测，提高系统对于卫星导航欺骗干扰的及时识别并隔离的能力。

Claims (4)

1.一种基于残差卡方-改进序贯概率比的抗欺骗导航方法，其特征在于，包括以下步骤：

步骤1：建立在导航系下的联邦滤波组合导航系统模型，设计联邦滤波器，经过联邦滤波器的最优估计，对导航系统进行校正输出；

步骤2：通过改进的残差χ²检验法对联邦滤波器中的子滤波器进行故障诊断和隔离，在残差χ²计算得到残差向量特征值的基础上对检测阈值进行模糊逻辑和加权处理；

步骤3：针对子滤波器中对导航结果起关键作用的SINS/GNSS子系统，采用改进的SPRT方法进行检测，并将结果和步骤2中的检测结果联合进行故障判决，判断GNSS是否受到欺骗干扰。

2.根据权利要求1所述的基于残差卡方-改进序贯概率比的抗欺骗导航方法，其特征在于，步骤1中所述建立在导航系下的联邦滤波组合导航系统模型，设计联邦滤波器，经过联邦滤波器的最优估计，对导航系统进行校正输出，具体如下：

(1.1)首先在导航坐标系下建立联邦滤波组合导航系统模型；惯性导航系统的误差模型选用通用误差模型，公共误差参考系统选取惯性导航系统，导航坐标系选取东北天地理坐标系，则系统的状态矢量X为

其中

为东北天方向上的姿态角误差；δV_E、δV_N、δV_U为东北天方向上的速度误差；δL、δλ、δH为纬经高方向的位置误差；ε_bx、ε_by、ε_bz为陀螺仪在三个轴向上的随机漂移；ε_rx、ε_ry、ε_rz为陀螺仪在三个轴向上的一阶马尔可夫过程随机噪声；

为加速度计在三个轴向上的常值偏差；

系统的状态方程

为：

其中F(t)为系统的状态转移矩阵，G(t)为误差系数矩阵，W(t)为白噪声随机误差矢量，具体形式如下：

状态转移矩阵F(t)：

式中F_N为系统误差矩阵，F_s为惯性器件的误差转换矩阵，F_M为惯性器件的噪声矩阵，具体如下：

其中I_3×3为3维的单位矩阵，

为姿态转换矩阵；

T_gx、T_gy、T_gz分别为陀螺仪三个轴的相关时间，T_ax、T_ay、T_az分别为加速度计三个轴的相关时间；

误差系数矩阵G(t)为：

白噪声随机误差矢量W(t)为：

W(t)＝[w_gx w_gy w_gz w_rx w_ry w_rz w_ax w_ay w_az]^T (8)

w_gx、w_gy、w_gz为陀螺仪的白噪声误差，w_rx、w_ry、w_rz为陀螺仪扩展误差模型的白噪声误差，w_ax、w_ay、w_az为加速度计的白噪声误差；

系统的量测方程Z为：

Z(t)＝H(t)X(t)+V(t) (9)

其中X(t)为系统的状态矢量，V(t)＝[N_E N_N N_U]^T，N_E、N_N、N_U分别为GPS接收机沿东、北、天方向的位置误差，H(t)根据不同的子滤波器选取；

将系统模型离散化，把状态方程式(2)和量测方程式(9)离散化，得

X_k为k时刻的状态，X_k-1为k-1的状态，W_k-1为离散的白噪声随机误差矢量，Z_k为k时刻的量测值，H_k为离散的量测矩阵，V_k为离散的量测噪声；

式中状态转移矩阵Φ_k|k-1和过程噪声矩阵Γ_k-1定义如下

式中，T为迭代周期，n为大于0的正整数；

(1.2)设计联邦子滤波器

每个子滤波器均采用卡尔曼滤波方式，具体如下：

状态一步预测方程：

式中，

为一步预测的状态估计值，

为k-1时刻的状态估计值；

状态估计方程：

式中，

为k时刻的状态估计值，K_k为k时刻的卡尔曼增益；

卡尔曼增益方程：

式中，P_k|k-1为一步预测均方差，R_k为量测噪声的方差矩阵；

一步预测均方差方程：

式中，Q_k-1为k-1时刻的系统噪声的方差矩阵；

均方差估计方程：

式中，P_k|k为状态估计的均方误差；

表示惯性导航系统的位置信息为

表示惯性导航系统的速度信息为

式中λ_I、L_I、H_I为惯性导航系统纬度、经度、高度的测量值，λ_t、L_t、H_t为载体纬度、经度、高度的真值，δλ、δL、δH为惯性导航系统经度、纬度、高度的测量误差值；v_IE、v_IN、v_IU为惯性导航系统在东北天方向上速度的测量值，v_E、v_N、v_U为载体在东北天方向上速度的真值，δv_E、δv_N、δv_U为惯性导航系统在东北天方向上速度的测量误差值；

(1.2.1)设计GNSS/SINS子滤波器

表示卫星导航系统的位置信息为

表示卫星导航系统的速度信息为

式中λ_G、L_G、H_G为卫星导航系统经度、纬度、高度的测量值，λ_t、L_t、H_t为载体经度、纬度、高度的真值，N_E、N_N、N_U为卫星导航系统在东北天方向上位置的测量误差值；v_GE、v_GN、v_GU为卫星导航系统在东北天方向上速度的测量值，v_E、v_N、v_U为载体在东北天方向上速度的真值，M_E、M_N、M_U为惯性导航系统在东北天方向上速度的测量误差值，R_M为子午圈半径，R_N为卯酉圈半径；

根据(1.1)知，GNSS/SINS子滤波器的状态方程为

GNSS/SINS子滤波器的量测方程为

其中，GNSS/SINS子滤波器的量测矩阵H_G和量测噪声矩阵V_G定义如下

V_G(t)＝[N_N N_E N_U M_E M_N M_U]^T (25)

量测噪声当作均值为零的白噪声处理，方差分别为

(1.2.2)设计里程计/SINS子滤波器

表示里程计的速度信息为

式中v_OE、v_ON、v_OU为里程计在东北天方向上速度的测量值，v_E、v_N、v_U为载体在东北天方向上速度的真值，O_E、O_N、O_U为里程计在东北天方向上速度的测量误差值；

由(1.1)得，里程计/SINS子滤波器的状态方程为

里程计/SINS子滤波器的量测方程为

其中，里程计/SINS子滤波器的量测矩阵H_O(t)和量测噪声矩阵V_O(t)定义为

H_O(t)＝[0_3×3 diag[1 1 1] 0_3×12] (29)

V_O(t)＝[O_E O_N O_U]^T (30)

量测噪声作为均值为零的白噪声处理；

(1.2.3)设计高度计/SINS子滤波器

表示高度计的高度信息为

h_a＝h_t-A (31)

式中h_a为高度计在高度上的测量值，h_t为载体在天向上位置的真值，A为高度计在高度上的测量误差值；

根据(1.1)知，高度计/SINS子滤波器的状态方程为

高度计/SINS子滤波器的量测方程为

Z_a(t)＝h_I-h_a＝H_a(t)X(t)+V_a(t) (33)

其中高度计/SINS子滤波器的量测矩阵H_a(t)和量测噪声矩阵V_a(t)定义为

H_a(t)＝[0_3×6 diag[0 0 1] 0_3×9] (34)

V_a(t)＝A (35)

量测噪声作为均值为零的白噪声处理；

(1.2.4)设计磁罗盘/SINS子滤波器

表示磁罗盘的角度信息为

ψ_c＝ψ_t-δψ (36)

式中ψ_c为磁罗盘测量的磁航向经过磁偏补偿后得到的航向角，ψ_t为载体航向角的真值，δψ为磁罗盘航向角的量测误差；

根据(1.1)知，磁罗盘/SINS子滤波器的状态方程为

磁罗盘/SINS子滤波器的量测方程为

Z_c(t)＝ψ_I-ψ_c＝H_c(t)X(t)+V_c(t) (38)

其中磁罗盘/SINS子滤波器的量测矩阵H_c(t)和量测噪声矩阵V_c(t)定义为

H_c＝[-sinψtanθ -cosψtanθ 1 0_1×12] (39)

V_c(t)＝δψ (40)

式中，ψ为偏航角，θ为俯仰角，量测噪声作为均值为零的白噪声处理；

(1.3)设计联邦滤波主滤波器

联邦滤波器是一种两级滤波结构，公共参考子系统SINS的输出X_k一方面给主滤波器，另一方面给各子滤波器作为量测值；各传感器子系统的输出只给各自的滤波器，各子滤波器的局部估计值X_i及其协方差阵P_i送至主滤波器与主滤波器的估计值一起进行融合，从而得到全局最优估计值；

由主滤波器和子滤波器合成的全局最优估计值

及其相应的协方差阵P_g被放大为

后再反馈给各子滤波器，来重置子滤波器的估计值，即：

主滤波器的协方差阵为

β_m和β_i分别是主滤波器与各子滤波器的信息分配系数，取值根据信息分配原则确定；

采用β_m＝0,β_i＝1/N，无重置结构的联邦卡尔曼滤波，各个局部滤波器独立完成滤波，仅对子滤波器估计值进行全局融合，并将全局最优估计值和协方差矩阵反馈到子滤波器中；

系统信息在主、子滤波器间的分配方法是基于信息分配原则的，且满足：

其中，β_m为主滤波器信息分配系数，β_i为第i个子滤波器对应的信息分配系数；

系统信息在子滤波器中的具体分配方法如下：

其中，

为主滤波器状态估计，P_g、Q_g分别为主滤波器状态协方差阵和系统噪声阵；

为子滤波器状态估计，P_i、Q_i分别为子滤波器状态估计对应的状态估计均方差阵和子滤波器系统噪声；

主滤波器最优信息融合：

式中

为第i个子滤波器的协方差矩阵；

为主滤波器中的协方差矩阵；

为主滤波器和子滤波器融合后的协方差矩阵；

为融合后的状态估计；

为第i个子滤波器的状态估计值；

为主滤波器的状态估计值；N为子滤波器的个数。

3.根据权利要求2所述的基于残差卡方-改进序贯概率比的抗欺骗导航方法，其特征在于，步骤2所述通过改进的残差χ²检验法对联邦滤波器中的子滤波器进行故障诊断和隔离，在残差χ²计算得到残差向量特征值的基础上对检测阈值进行模糊逻辑和加权处理；具体如下：

(2.1)残差χ²检验

用残差χ²检验对系统的故障进行检测和隔离；由于每个子滤波器均为卡尔曼滤波器，新息向量也即预测残差为：

其中一步预测值

当无故障发生时卡尔曼滤波器的残差r_k是零均值的白噪声，方差为

当系统发生故障时，残差的均值就不再为0，因此对残差的检验能够确定系统是否发生了故障；

对r_k做二元假设：

H₀无故障时：

式中，r_k为子滤波器的预测残差，A_k为残差方差；

H₁有故障时：

式中，μ为常数向量；

定义故障检测函数为：

式中λ_k为故障检测值，是服从自由度为m的χ²分布，即λ_k～χ²(m)，m为量测矩阵Z_k的维数；

故障判定准则为：

其中设置的阈值T_D由误警率P_f确定，给定误警率P_f＝α，则由如下计算公式求出T_D：

式中，λ为均值，n为方差；

(2.2)改进的残差χ²检验法

采用改进的残差χ²检验法，在算出残差向量特征值的基础上对检测阈值进行模糊逻辑和加权处理，能够提高系统对缓变故障的检测能力；

对残差向量进行归一化处理：

式中r_k和A_k与(2.1)中定义相同，

为k时刻归一化后的残差向量；

从当前时刻开始往前选取n个时刻的残差向量组成残差矩阵R：

定义残差的方差矩阵Q为：

Q＝R^TR (52)

根据控制理论知，Q矩阵特征值的最大值有界，即：

λ_max≤S (53)

通过选取最近的N个矩阵最大特征值的平均值进一步提高数据的可靠性，具体的表达式如下：

式中，

为N个矩阵最大特征值的平均值，λ_max(k)为k时刻Q矩阵的最大特征值；

考虑实时性和计算量的问题，建议N与n的值不要取得过大；

此时，故障判断准则为：

(2.3)当某个子滤波器检测出故障时，系统舍弃该子滤波器的信息，主滤波器对剩余的子滤波器进行重构和信息融合，使得剩余的导航子系统能够继续提供导航信息，保证了多源组合导航系统具有一定的容错性能。

4.根据权利要求3所述的基于残差卡方-改进序贯概率比的抗欺骗导航方法，其特征在于，步骤3中所述针对子滤波器中对导航结果起关键作用的SINS/GNSS子系统，采用改进的SPRT方法进行检测，并将结果和步骤2中的检测结果联合进行故障判决，判断GNSS是否受到欺骗干扰，具体如下：

(3.1)残差χ²检测

利用(2.1)中阐述的残差χ²检测法对SINS/GNSS子滤波器的估计结果进行故障诊断；

(3.2)改进SPRT检测

假设某未知正态随机变量x的k次序贯独立样本为{x_i|i＝1,2,…,k}，根据概率统计理论，近似有

且

式中x_i为第i次样本值，

为样本均值，既是有故障时的均值，也为无故障时的均值；

为样本方差；

假定状态变量x的量测实际值为x^*，则定义如下两个事件

其中，x₀为故障状态下的真实值；

量测序列x₁,x₂,…,x_k属于两个样本类H₀-正常类和H₁-故障类之一，概率密度函数为

式中，p(x_i|H₀)和p(x_i|H₁)分别为两个样本类的概率密度；

则得似然比L(k)：

对式(59)取自然对数，得到对数似然比λ(k)：

将式(60)简化，得：

由式(61)知，λ(k)不可能为负值，当系统处于正常工作状态时，随着k的增加，样本均值

将逐渐逼近x₀，因此λ(k)→0；当系统发生故障时，随着k的增加，样本均值

依然会逐渐逼近故障情况下的真值，但真值不为x₀；

综上，得到似然比λ(k)的递推计算公式：

式中，

为k-1个样本的均值，

为此时的样本方差；

若变量x的先验方差为σ²，则用σ²替代上式中的

若不确定变量x的先验方差，则用实时计算的样本方差替代；

Wald序贯概率算法中给出了两个检测阈值，分别为T(H₀)和T(H₁)；对于上述改进的SPRT算法，则仅需设定一个检测阈值T(H₁)，该值的计算方式与Wald算法一致，即由误警率P_f和漏检率P_m确定：

则判定准则为

(3.3)故障判决

上述(3.1)和(3.2)两种故障检测方法是并行计算的；对SINS/GNSS子滤波器的残差值分别使用残差χ²和改进SPRT两种检测法进行故障检测，得到两个检测结果，并且将当前时刻的改进SPRT检测值与上一时刻的作差得到增量；把增量和两个检测结果一同输入到系统故障综合决策规则中，进行最终故障判决，从而得出最终的检测结果；系统故障综合决策规则见表1，表中用“是”与“否”表示是否发生故障；

表1 系统故障综合决策规则

由此得到对SINS/GNSS子滤波器估计结果进行故障检测的结果。

Images