一种适用于卫导系统缺星情况下的完好性监测方法

摘要

本发明公开了一种适用于卫导系统缺星情况下的完好性监测方法，属于卫星导航领域，该方法根据卫导系统完好性的需求，结合惯导精度选取滑动窗窗口大小、确定检测门限；然后利用基于滑动累积的多样本χ2检验方法对卫导系统进行故障检测，判断故障情况和确定故障卫星；本发明有效解决了卫导系统在缺星条件下无法实现故障监测的问题，提高了卫导系统的完好性；应用本发明提供的监测方法，故障偏差在40米之后，故障检测率近似100％；故障幅度在大于100米后，识别率近似100％。

说明书

技术领域

本发明属于卫星导航领域，具体地说，是指一种应用于卫导系统缺星情况下惯导辅助接收机自主完好性监测方法。

背景技术

全球卫星导航系统(Global Navigation Satellite System，GNSS，简称卫导系统)具有在全球范围内提供高精度三维位置和速度信息的能力。目前卫导系统在军用领域和民用领域都有着重要的应用，现已广泛应用于航空、航天、航海、车辆导航、导弹制导、大地测量、野外救生等诸多领域。

完好性是评价导航系统性能的一个重要因素，是系统安全的保障，尤其在航空领域中，完好性直接关系到生命安全。卫导系统在各领域的广泛应用必须要保证其完好性的要求，而完好性监测技术提供了完好性的保障。

接收机自主完好性监测(RAIM)算法提出一定程度地解决故障监测的问题。然而，RAIM算法必须有足够的冗余观测信息进行故障检测和识别，在卫星数少于5颗时RAIM算法不能实现故障监测。

北斗一代系统(BD-1)为我国自主研发的导航系统，目前已经进入老龄化阶段，完好性监测问题显得尤为重要。但由于BD-1卫星数少于5颗，传统的RAIM算法无法实现故障监测。针对BD-1缺星无法实现故障监测的问题，目前已有基于气压高度计/高精度原子钟辅助BD-1系统的完好性监测技术。然而，此技术仅有5维观测信息只能完成故障检测而无法进行故障识别。

此外，美国的GPS(Global Positioning System)、俄罗斯的GLONASS(Global OrbitingNavigation Satellite System)、欧盟的Galileo系统，虽然由于其卫星数多、几何分布好，在多数情况下可以保证5颗以上的可见卫星数，但在卫星信号容易受到干扰或遮挡的复杂环境下，缺星问题仍然存在，完好性无法保障。因此，适用于卫导系统缺星情况下的完好性监测方法具有广阔的应用前景。

发明内容

本发明提出的方法解决了由于卫导系统卫星数少于5颗，传统的RAIM算法无法实现故障监测的问题，该方法基于滑动累积的多样本χ²检验方法通过引入惯导信息对卫导系统卫星多个时刻观测信息进行一致性检验实现故障监测。

本发明提供的方法包括如下步骤：

步骤一：根据卫导系统完好性的需求，结合惯导精度选取滑动窗窗口大小、确定检测门限；

步骤二：利用基于滑动累积的多样本χ²检验方法对卫导系统进行故障检测，判断故障情况；

判断当前滑动窗是否达到设定大小N，若未满足要求，则进行初始化阶段；若已满足要求，则进入故障检测阶段；

所述的初始化阶段，依靠惯导外推给出当前时刻的定位解。

所述的故障检测阶段，在当前时刻k≥N时，滑动窗满足要求达到设定的大小N，此时进行故障监测，具体实施过程如下：

(1)构造统计检测量l_k；

$l_k=\underset{i=k-N+1}{\overset{k}{\Sigma }}{r}_i^T{A}_i^{-1}{r}_i$

其中，A_i为i时刻新息序列r_i的方差，i＝k-N+1，k-N+2，...k；

(2)统计检测量l_k与设置的检测门限T_D作比较，来判断是否有故障；

在故障检测过程中，当滑动窗滑动到k时刻，如果检测k时刻无故障，则利用k-N+1时刻卫导系统的所有观测信息进行滤波处理得到k-N+1时刻的定位解，转步骤四；如果检测k时刻有故障，并且是初次判定有故障，则转(3)；否则转步骤(4)；

(3)利用惯导外推一个时刻，继续进行下一时刻的故障监测；重复执行此过程，直至N个时刻后，如果仍然监测到故障，则进入步骤三实时故障识别算法，以确定故障星。

(4)如果K+N时刻起，故障持续出现，则直接进入步骤三；如果系统判断故障初次出现后又消失，利用惯导外推到一个时刻，继续进行下一时刻的故障判断，重复N次，如N个时刻后仍判断无故障再进入步骤四。

步骤三：进行故障识别，确定故障卫星；

选取故障识别算法确定故障卫星，将故障卫星进行隔离，剔除识别出的故障星后同样进行滤波处理得到k-N+1时刻的定位解，完成了k时刻的故障检测。

步骤四：滑动窗向前滑动，进行下一时刻的故障监测。

本发明提供的方法有以下优点：

(1)本发明有效解决了卫导系统在缺星条件下无法实现故障监测的问题，提高了卫导系统的完好性。

(2)应用本发明提供的监测方法，故障偏差在40米之后，故障检测率近似100％；故障幅度在大于100米后，识别率近似100％。

附图说明

图1是本发明完好性监测方法流程图；

图2是基于滑动累积的多样本χ²检验方法流程图；

图3是实施例中故障检测率曲线图；

图4是实施例中故障识别率曲线图。

具体实施方式

本发明是一种应用于卫导系统的完好性监测的方法，应用本发明提供的方法可以有效解决卫导系统在缺星条件下的故障监测问题。本发明采用滑动累积的方法，通过开辟一定大小窗口并对窗口内多个时刻卫星观测信息按照一定的检验准则进行一致性检验，实现对卫导系统的完好性监测，监测方法流程如图1所示，具体通过如下步骤实现：

步骤一：根据卫导系统完好性的需求，结合惯导精度选取滑动窗窗口大小、确定检测门限；

选取滑动窗窗口大小为N，则滑动窗内可存储N个时刻的卫星观测信息。这N个时刻的卫星观测信息包含当前时刻(k时刻)观测信息和先前N-1个时刻的观测信息，其中当前时刻观测信息存入滑动窗最前端，而前一时刻(k-1时刻)至最早时刻(k-N+1时刻)共N-1个观测信息依次向后存入滑动窗。

滑动窗的窗口大小N的选取需要考虑几方面的因素：

(1)N的大小必须在卫导系统告警时限内；

(2)随着N的增大，故障检测性能会因所检测信号信噪比的增加而提升；N选择过小会导致检测性能较差；

(3)N的选取与惯导精度有重要关系，N的增大使惯导外推时间加长，这不但会导致惯导系统定位精度下降使其超出定位告警极限造成完好性风险，而且会抑制信噪比的增加，尤其是对于较低等级的惯导系统，N的增大会导致惯导系统精度的明显下降。综合上述的三个因素，N的选取需要对卫导系统完好性所规定的告警极限、告警时限和惯导系统精度几方面因素综合考虑。通常N的选取范围可选为5～10之间。

在选定N的大小后，检测门限T_D可通过限定虚警率P_f＝α来确定，根据Neyman-Pearson准则，有：

P_f＝P[l＞T_D|H₀]＝α                                                 (1)

式(1)表示在无故障时统计检测量l大于检测门限T_D的发生概率，即虚警率。式(1)中统计检测量l为服从χ²分布的随机变量，H₀表示无故障假设，则虚警率P_f具体可写为：

$P_f=1-\underset{0}{\overset{T_D}{\int }}\frac{1}{2^{n/2}\Gamma (n/2)}{l}^{\frac{n}{2}-1}{e}^{-\frac{1}{2}}\mathrm{dl}=\alpha ---\left(2\right)$

式(2)中：为伽马函数，n为统计检测量l所服从χ²分布的自由度，α为卫导系统最大虚警率，根据航空领域中对卫导系统所规定的工作性能最低标准选取。通过对式(2)数值解算可以获得检测门限T_D。

步骤二：利用基于滑动累积的多样本χ²检验方法对卫导系统进行故障检测，判断故障情况；

判断当前滑动窗是否达到设定大小N。若未满足要求，即此时滑动窗大小是小于N，则进行初始化阶段；若已满足要求，即此时滑动窗大小等于N，则进入故障检测阶段。具体如下：

情况一：初始化阶段；

卡尔曼滤波器达到稳定状态需要经历一段时间，在达到稳定状态前由于卫导和惯导组合系统的精度较差而无法保证故障检测的有效性。以卫导和惯导组合系统达到稳定状态时刻为0时刻，在当前时刻k＜N时，由于滑动窗还未开辟到N大小，不能积累足够的观测信息，无法实现故障监测；所以此时只能依靠惯导外推给出当前时刻的定位解，即利用惯导将上一时刻定位解(k-1时刻)外推一个时刻给出当前时刻的定位解。

情况二：正常工作状态；

在当前时刻k≥N时，滑动窗满足要求达到设定的大小N，此时卫导和惯导组合系统进入正常工作状态可以正常的对卫导系统进行故障监测，具体实施过程如图2所示：

(1)构造统计检测量l_k。

卫导系统与惯导组合系统的离散化卡尔曼滤波器方程如下：

$\left(\begin{array}{c}{X}_k=F_{k-1}{X}_{k-1}+W_{k-1}\\ Z_k=H_k{X}_k+V_k\end{array}\right)$

式中，X_k为k时刻的状态向量，即由惯导三维姿态、速度、位置误差项，陀螺仪、加速度计误差项，卫星钟差、频差项组成。F_k-1为k-1时刻到k时刻的卫导与惯导组合系统一步转移矩阵；W_k-1为k-1时刻的卫导与惯导组合系统噪声协方差矩阵；Z_k为表示k时刻的观测量，即由卫导系统卫星的观测伪距和相应的惯导估计伪距之差组成；H_k表示k时刻的观测矩阵；V_k为k时刻的量测噪声协方差矩阵。卫导与惯导组合系统噪声和量测噪声是互不相关的零均值的白噪声序列。

卡尔曼滤波k时刻的新息序列r_k定义为：

$r_k=Z_k-H_k{\hat{X}}_k^{-}$

式中，表示由k-1时刻的滤波值得到的状态一步预测值。表示由观测矩阵和状态一步预测值得到的量测值的预测值。

本发明利用多个时刻新息序列构造统计检测量，具体方法如下：

在k时刻，组合系统从k-N+1时刻至k时刻共N个观测量(Z_k-N+1，Z_k-N+2，…Z_k-1，Z_k)，根据k-N时刻的定位解，利用惯导外推得到的N个观测量的预测值(i＝k-N+1，k-N+2，...k)，对i时刻进行新息计算，得到先前N-1个时刻的新息序列r_i，即：

$r_i=Z_i-H_i{\hat{X}}_i^{-},$

其中，i＝k-N+1，k-N+2，...k。

则构造统计检测量l_k为：

$l_k=\underset{i=k-N+1}{\overset{k}{\Sigma }}{r}_i^T{A}_i^{-1}{r}_i$

式中，A_i为i时刻新息序列r_i的方差，i＝k-N+1，k-N+2，...k。

(2)统计检测量l_k与设置的检测门限T_D作比较，来判断是否有故障。

故障检测可以归结为一个二元检测问题：即通过由样本信息构造出的统计检测量l_k与设置的检测门限T_D作比较，来判定H₁(有故障)和H₀(无故障)两种假设哪一个成立。

设定检测门限为T_D，并制定故障判定准则为：

在故障检测过程中，当滑动窗滑动到k时刻，如果检测k时刻无故障，则利用k-N+1时刻卫导系统的所有观测信息进行滤波处理得到k-N+1时刻的定位解，转步骤四；如果检测k时刻有故障，并且是初次判定有故障(即从无故障到故障出现这一过程)，则转(3)；否则转步骤(4)。

(3)故障有可能只是在k时刻或先前几个时刻后才出现，即k-N+1时刻以后和k时刻之前一段时间的故障情况无法确定。在这种情况下组合系统并不立即进行故障识别和滤波处理，而是利用惯导外推一个时刻，继续进行下一时刻的故障监测。重复执行此过程，即利用惯导外推N个时刻，直至N个时刻后，如果仍然监测到故障，则进入步骤三实时故障识别算法，以确定故障星。

(4)如果K+N时刻起，故障持续出现，则直接进入步骤三；如果系统判断故障初次出现后又消失，即从有故障到判断无故障这一过程，也不立即进行滤波处理，而是利用惯导外推到一个时刻，继续进行下一时刻的故障判断，重复N次，如N个时刻后仍判断无故障再进入步骤四。

步骤三：选用合适的故障识别算法，以确定故障卫星；

如果经步骤三检测到k时刻有故障则进行故障隔离，即经过步骤三的故障识别算法，剔除识别出的故障星后同样进行滤波处理得到k-N+1时刻的定位解。此时，即完成了k时刻的故障检测。

判断出故障后可以用(m-1)法或比较新息分量的方法来识别故障星。不论是(m-1)法还是比较新息分量的方法，都需要对N个时刻的新息分量综合考虑，而不是只对故障时刻的新息分量进行处理。本实施例中采用比较新息分量的方法进行故障识别，即将滑动窗口中的N个时刻的新息序列相加，比较此时得到新的变量中的p(p即为卫导系统卫星数)个分量，即比较中的p个分量的大小，如图2，最大分量对应的卫星视为故障星。

步骤四：滑动窗向前滑动，进行下一时刻的故障监测。

滑动窗向前滑动一个窗口，此时滑动窗剔除末端(k-N+1时刻)的观测信息并将最新时刻(k+1时刻)的观测信息存储到最前端。此时即可实施最新时刻的故障监测。

本发明通过滑动累积的多样本χ²检验方法实现了对卫导系统的故障监测，虽然存在一定的检验延迟，但仍能保证其在所要求性能之内。可见，本发明能有效地解决卫导系统缺星条件下的完好性监测问题。

实施例1

由于BD-1只有3颗卫星，对于完好性监测缺星为其固有问题，且目前已经进入老龄化阶段，完好性监测问题显得尤为重要。因此，以BD-1作为研究实例进一步说明本发明提供的方法，有以下具体实施步骤：

步骤一、根据卫导系统完好性性能要求，结合惯导精度选取滑动窗口大小、确定检测门限；

本实例中惯导精度较高并经过气压高度计阻尼，其中陀螺仪零偏0.02度/小时，加速度计零偏100mg。

完好性性能参数选取，根据国际航空无线电技术委员会(Radio Technical Commission forAeronautics，RTCA)所制定的对于GPS在民用航空飞行时完好性保障方面的最低工作标准(Minimum Operational Performance Standards，MOPS)，如下表：

表1RTCA GPS最低工作状态标准(MOPS)

  飞行阶段  水平定位告警极限  (m) 最大虚警率(最大允许  告警率)(hr) 告警时限(s)  越洋  7412.8  0.002  30  航路  3706.4  0.002  30  终端  1853.2  0.002  10  非精密进近  555.96  0.002  10

根据实际定位观测间隔(本实例为1s)，最大虚警率为α＝0.002/3600≈5.56*10^-7。针对表1中非精密进近飞行阶段的情况，其告警时限为10s，在10s的范围内滑动窗窗口N选择过小会导致检测性能较差；而对于精度较好的惯导，N可以选取较大的值(但必须小于10，过大的N值会系统导致精度下降从而超出定位告警极限)，本实例选取滑动窗大小N＝6。此时，统计检测量l在无故障时服从自由度为n＝N·p＝18(p＝3)的中心χ²分布，即l～χ²(0，18)。根据，P_f＝P[l＞T_D|H₀]＝α，解得门限T_D＝63.47。

步骤二：利用基于滑动累积的多样本χ²检验方法对卫导系统进行故障检测，判断故障情况；

系统正常工作状态时，BD-1卫星的故障检测。假定BD-1的中星(位于BD-1三颗同步卫星中间的称为中星)在k时刻出现故障，在其正常伪距基础上加入偏差其余卫星无故障。

情况一：加入较小偏差40米；

在故障发生时刻k，统计检测量l_k＝29.19；k+1时刻，l_k+1＝53.31；k+2时刻，l_k+2＝59.44，此三个时刻l＜T_D，都未检测到故障；而在k+3时刻，l_k+3＝86.84，l_k+3＞T_D＝63.47，即检测到故障。这是由于在k+3时刻统计检测量中包含了三个时刻故障卫星信息，增加了检测信号的样本数从而实现了故障检测。在此后时刻，均有l＞T_D即都能正常的检测到故障。

虽然上述监测方法中有3秒的故障检测延迟，然而由于在MOPS所要求的告警时限(10s)之内，所以故障检测算法仍然有效。

情况二：加入较大偏差100米；

在故障发生时刻k，统计检测量l_k＝138.49。由于l_k＞T_D＝63.47，即检测到故障。同样在此后时刻，均有l＞T_D，故障检测正常。

需要说明的是：对于情况一、情况二，在故障监测的过程中，若故障首次出现，则应利用惯导外推一个时刻，继续进行下一时刻的故障监测。重复执行此过程，即利用惯导外推N个时刻，直至N个时刻后，如果仍然监测到故障，即仍有l＞T_D则进行故障识别算法。

步骤三、比较新息分量的方法来识别故障星。

将k-N+1时刻至k时刻共N个时刻的新息序列相加，得到一个新的变量，其三个分量的绝对值分别为43.44、139.22、39.53。其中最大的分量139.22对应中星，即识别中星为故障星。由于中星设置为故障星，即比较新息分量的方法实现了故障的识别。

实施例2

选择BD-1中星为故障卫星，设置故障偏差由0m以10m为步长增长到150m。用户位置设置为北京(东经116°，北纬40°)附近，进行一段5小时的仿真，得到故障检测率曲线如图3所示，故障识别率曲线如图4所示。

由图3可以看出，故障偏差在40米之后，故障检测率近似100％。由图4可以看出，故障幅度在大于100米后，识别率近似100％。即在故障幅度大于100米后，本算法能很好的完成故障监测。