CN110907953B - 一种卫星故障识别方法、装置及软件接收机 - Google Patents

Abstract

本发明公开了一种卫星故障识别方法、装置及软件接收机，涉及卫星导航与定位技术领域；该方法包括：采集接收机位置数据和钟差，根据接收机位置数据、钟差构建观测矩阵并采用奇偶空间算法对观测矩阵进行分解，生成奇偶空间矩阵；对奇偶空间矩阵进行矩阵运算，得到参与定位的每个卫星的斜率并提取最大斜率值；基于最大斜率值和误警检测值计算定位卫星的相似径向误差保护值，当相似径向误差保护值不小于预设值时，根据奇偶空间矩阵对应的卫星偏差数据计算校验统计量；根据校验统计量、误警检测值和预置的阈值因子检测参与定位的卫星中是否存在故障；本发明具有时间复杂度低、空间复杂度低的特点，能够适应常态飞行、高动态飞行的卫星故障检测和识别。

Description

一种卫星故障识别方法、装置及软件接收机

技术领域

本发明属于卫星导航与定位技术领域，更具体地，涉及一种卫星故障识别方法、装置及软件接收机。

背景技术

在软件接收机的用户位置、速度以及时间信息(Position Velocity Time，PVT)解算 中，需要实时地输出用户经度、纬度、高度和速度等测量信息。为了能够提高接收机的导航精度，必须实时地对卫星进行监测。然而复杂环境下多星座可见卫星的变化容易引 起系统故障，导致出现检测性能不稳定的问题。因此，需要设计一种接收机监测、识别 故障卫星的算法(RAIM，Receiver Autonomous Integrity Monitoring)，当观测信息的误差 超过阈值后，将启动该算法使其能够成功剔除故障卫星，从而保证用户实时接收高精度、 高可靠性的定位信息。

启动RAIM算法至少需要5颗卫星参与定位，进行故障卫星剔除至少需要6颗卫星参与定位。

传统的RAIM算法主要包括最小二乘残差算法、奇偶空间校验的算法，二者在数学上是等价的，故二者具有相同的卫星伪距跳变识别度。但是二者的运算角度具有较大的 差别，最小二乘残差是在欧式距离上进行计算的，计算的空间复杂度高；而奇偶空间校 验是在奇偶空间上进行计算的，计算的时间复杂度高。由于弹载接收机内部DSP资源 有限，因此对算法的时间复杂度、空间复杂度都具有较高要求。

因此，研制一种时间复杂度低、空间复杂度低，并且能够适应常态飞行和高动态飞行，并且具有较高的卫星伪距跳变识度的算法至关重要。

发明内容

针对现有技术的至少一个缺陷或改进需求，本发明提供了一种卫星故障识别方法、 装置及软件接收机，结合最小二乘残差算法和奇偶空间校验算法的各自优点，具有时间 复杂度低、空间复杂度低的特点，并且具有相同的卫星伪距跳变识别度，使其能够适应常态飞行、高动态飞行的卫星故障检测和识别；其目的在于解决现有方法存在的时间复 杂度、空间复杂度高的问题。

为实现上述目的，按照本发明的第一个方面，提供了一种卫星故障识别方法，包括以下步骤：

S1：采集基于最小二乘算法进行卫星定位时收敛的接收机位置数据和接收机钟差， 根据所述接收机位置数据、接收机钟差构建观测矩阵并采用奇偶空间算法对所述观测矩 阵进行分解，生成奇偶空间矩阵；

S2：采用奇偶空间算法对所述标准正交矩阵进行分解，生成奇偶空间矩阵；

S3：对所述奇偶空间矩阵进行矩阵运算，得到参与定位的每个卫星的斜率并提取最 大斜率值；

S4：根据所述最大斜率值和误警检测值计算相似径向误差保护值，当所述相似径向 误差保护值不小于预设值时，根据所述奇偶空间矩阵对应的卫星偏差数据计算校验统计 量；

S5：根据所述校验统计量、误警检测值和预置的阈值因子检测参与定位的卫星中是 否存在故障。

优选的，上述卫星故障识别方法，其步骤S3中具体包括：

(1)根据所述奇偶空间矩阵构造第一矩阵和第二矩阵；

B＝I_n-Q_p ^TQ_p

A＝H\B；

其中，Q_p表示奇偶空间矩阵；B表示第一矩阵；I_n表示单位矩阵；A表示第二矩阵； H表示观测矩阵；n表示参与定位的卫星数量；

(2)根据所述第一矩阵和第二矩阵计算参与定位的每个卫星的斜率并取其最大斜率值Slope_max；其中，

优选的，上述卫星故障识别方法，根据预置的正态分布方差、误报警率计算所述误警检测值：

其中，σ_γ表示误警检测值；

表示正态分布方差；P_FA表示误报警率；n表示参与 定位的卫星数量；icdf('ch2',1-P_FA,n-4)表示求取以(n-4)为自由度，累计概率为 (1-P_FA)的χ²分布的临界值。

优选的，上述卫星故障识别方法，其步骤S5中具体包括：

若

则判定参与定位的卫星中存在故障；否则，判定参与定位的卫星中不存在故障；

其中，F_SSE表示校验统计量；σ_γ表示误警检测值；β表示阈值因子；n表示参与定 位的卫星数量。

优选的，上述卫星故障识别方法，还包括排除故障卫星的步骤，具体为：

(1)根据预置的误报警率计算检测门限值；

T_D＝Φ^-1(1-P_FA/2n)

其中，T_D表示检测门限值；P_FA表示误报警率；Φ^-1表示正态分布的标准值；n表 示参与定位的卫星数量；

(2)若r_imax＝max(|P^T·Q_p(:,i)|/(σ₀·norm(Q_p(:,i))))＞T_D，表明第i颗卫星的门限值 大于检测门限值，则判定第i颗卫星存在故障并将其剔除；否则判定参与定位的卫星中 不存在故障卫星。

优选的，上述卫星故障识别方法，所述相似径向误差保护值为最大斜率值与误警检 测值的乘积。

优选的，上述卫星故障识别方法，所述预设值为中国地区飞机非精密进场值，当 n＜＝6,ARP_m(n)＝337；当n＝7,ARP_m(n)＝343；当n＞＝7,ARP_m(n)＝347；

其中，ARP_m(n)表示飞机非精密进场值；n表示参与定位的卫星数量。

优选的，上述卫星故障识别方法，所述校验统计量的计算方法为：

F_SSE＝P^TP，P＝Q_Pε；

其中，F_SSE表示校验统计量；P表示奇偶空间矢量，ε表示伪距噪声矢量。

按照本发明的第二个方面，还提供了一种卫星故障识别装置，包括：

创建单元，用于采集基于最小二乘算法进行卫星定位时收敛的接收机位置数据和接 收机钟差，根据所述接收机位置数据、接收机钟差构建观测矩阵并对所述观测矩阵进行正交三角分解，得到标准正交矩阵；采用奇偶空间算法对所述标准正交矩阵进行分解， 生成奇偶空间矩阵；

第一计算单元，用于对所述奇偶空间矩阵进行矩阵运算，得到参与定位的每个卫星 的斜率并提取最大斜率值；

第一计算单元，用于根据所述最大斜率值和误警检测值计算相似径向误差保护值， 当所述相似径向误差保护值不小于预设值时，根据所述奇偶空间矩阵对应的卫星偏差数 据计算校验统计量；

检测单元，用于根据所述校验统计量、误警检测值和预置的阈值因子检测参与定位 的卫星中是否存在故障。

按照本发明的第三个方面，还提供了一种软件接收机，其包括至少一个处理单元、以及至少一个存储单元，

其中，所述存储单元存储有计算机程序，当所述程序被所述处理单元执行时，使得所述处理单元执行上述任一项卫星故障识别方法的步骤。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益 效果：

(1)本发明提供的卫星故障识别方法、装置及软件接收机，通过结合最小二乘残差算法和奇偶空间校验算法的各自优点，使得新算法不仅与传统最小二乘算法，或者奇 偶空间校验算法相比，具有时间复杂度低、空间复杂度低的特点，并且具有相同的卫星 伪距跳变识别度，使其能够适应常态飞行、高动态飞行的卫星故障检测和识别。

(2)本发明提供的卫星故障识别方法、装置及软件接收机，在相同的误警率和漏警率的前提下，通过在卫星故障检测和识别时设置阈值因子，增大了卫星的最小伪距跳 变识别度。

附图说明

图1是本发明提供的卫星故障识别方法的流程图；

图2是本发明实施例提供的卫星故障识别方法的算法流程图；

图3是本发明实施例提供的软件接收机通过接收卫星信号模拟器发射的卫星信号， 通过模拟器上位机，控制GPS卫星均伪距不跳变的定位效果图；

图4是本发明实施例提供的软件接收机通过接收卫星信号模拟器发射的卫星信号， 通过模拟器上位机，控制BD卫星均伪距不跳变的定位效果图；

图5是本发明实施例提供的软件接收机通过接收卫星信号模拟器发射的卫星信号， 通过模拟器上位机，控制GPS+BD卫星均伪距不跳变的组合定位效果图；

图6是本发明实施例提供的软件接收机通过接收卫星信号模拟器发射的卫星信号， 通过模拟器上位机，控制GPS的27#卫星伪距跳变45m，使用卫星故障识别算法剔除故障卫星的效果图。

图7是本发明实施例提供的软件接收机通过接收卫星信号模拟器发射的卫星信号， 通过模拟器上位机，控制BD的103#卫星伪距跳变45m，使用卫星故障识别算法剔除故障卫星的效果图。

图8是本发明实施例提供的软件接收机通过接收卫星信号模拟器发射的卫星信号， 通过模拟器上位机，控制GPS+BD卫星组合定位时，GPS伪距跳变45m，使用卫星故 障识别算法剔除故障卫星的效果图。

图9是本发明实施例提供的软件接收机通过接收卫星信号模拟器发射的卫星信号， 通过模拟器上位机，控制GPS+BD卫星组合定位时，BD的103#卫星伪距跳变45m，使 用卫星故障识别算法剔除故障卫星的效果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征 只要彼此之间未构成冲突就可以相互组合。

本发明所提供的一种卫星故障识别方法，联合最小二乘残差算法和奇偶空间校验算 法，减少算法的时间复杂度和空间复杂度；另外通过设置阈值因子，使识别算法具有较高的故障卫星伪距跳变识别度。

图1是本发明提供的卫星故障识别方法的流程图，参见图1所示，该方法包括以下步骤：

S1：采集基于最小二乘算法进行卫星定位时收敛的接收机位置数据和接收机钟差， 根据所述接收机位置数据、接收机钟差构建观测矩阵H并对所述观测矩阵H进行正交三角分解，得到标准正交矩阵Q；

假设某时间内有n颗卫星参与定位，若n＜5，则不启动RAIM算法；若n＞5，则 启动RAIM算法，同时构建采用最小二乘算法进行定位时收敛的观测矩阵H；对观测矩 阵H进行正交三角分解(QR分解)，得到标准正交矩阵Q和上三角矩阵R。

S2：采用奇偶空间算法对所述标准正交矩阵Q进行分解，生成奇偶空间矩阵Q_p；

本实施例中，提取标准正交矩阵Q的转置矩阵Q^T的5～n行数据，生成奇偶空间矩阵Q_p。

S3：根据所述奇偶空间矩阵Q_p计算参与定位的每个卫星的斜率并提取最大斜率值Slope_max；具体包括：

(1)根据奇偶空间矩阵Q_p构造第一矩阵B，使得B＝I_n-Q_p ^TQ_p；然后构造第二矩 阵A，使得A＝H\B；

(2)根据第一矩阵B和第二矩阵A计算参与定位的每个卫星的斜率并取其最大值Slope_max；

在软件接收机中，最难检测的参与卫星定位的故障卫星就是其能够产生最大斜率的 卫星。因此，本发明采用相似径向误差保护值ARP作为定位卫星在空间的几何分布作为衡量标准；首先计算参与定位的每颗卫星的斜率，并取其最大值Slope_max；

其中，

S4：根据最大斜率值Slope_max和误警检测值σ_γ计算相似径向误差保护值ARP，当相似径向误差保护值ARP不小于预设值时，根据奇偶空间矩阵Q_p计算校验统计量F_SSE；

首先，根据预置的正态分布方差

误报警率P_FA计算误警检测值σ_γ：

其中，icdf('ch2',1-P_FA,n-4)表示求取以(n-4)为自由度，累计概率为(1-P_FA)的χ²分布的临界值。

相似径向误差保护值ARP为最大斜率值Slope_max与误警检测值σ_γ的乘积，即ARP＝Slope_max·σ_γ；

然后进行RAIM算法可用性判断，若ARP＝Slope_max·σ_r≥ARP_m(n)，则RAIM算法 可用；否则不可用。其中ARP_m(n)是中国地区飞机非精密进场值，当n＜＝6,ARP_m(n)＝337。 当n＝7,ARP_m(n)＝343。当n＞＝7,ARP_m(n)＝347；

当判定RAIM算法可用时，根据奇偶空间矩阵Q_p计算校验统计量F_SSE，具体的：

F_SSE＝P^TP，P＝Q_Pε；

其中，P表示奇偶空间矢量，表示卫星的偏差数据；ε表示伪距噪声矢量。

S5：根据校验统计量F_SSE、误警检测值σ_γ和预置的阈值因子检测参与定位的卫星中是否存在故障具体包括：

若

则判定参与定位的卫星中存在故障；否则，判定参与定位的卫星中不存在故障；其中，β是根据具体的接收机性能设置的阈值因子，本实施例 中的取值范围为0＜β＜20。

S6：排除故障卫星；

根据预置的误报警率P_FA计算检测门限值T_D；若第i颗卫星的门限值大于检测门限值T_D，则判定第i颗卫星存在故障并将其剔除；否则判定参与定位的卫星中不存在故障 卫星；具体的，若r_imax＝max(|P^T·Q_p(:,i)|/(σ₀·norm(Q_p(:,i))))＞T_D，则表明第i颗卫星 存在故障并将其剔除；

其中，T_D＝Φ^-1(1-P_FA/2n)，Φ^-1表示正态分布的标准值；

考虑到软件接收机快速定位的需求，通常将可见卫星数n存储在数组中，接收机DSP程序在烧写到FLASH之前，就可计算出误警检测值数组σ_r'[n-4]、检测门限数组 T_D'[n-4]中，有效避免重复计算，提高定位速度。

下面结合附图和具体的实例对本发明提供的卫星故障识别方法作更进一步的详细 说明。

实施例一

图2较为详细地描述了本发明提出的卫星故障识别方法的具体的算法流程，参见图 2所示，具体流程如下：

(1)软件接收机通过天线接收卫星信号后，首先启动接收机内FLASH中的存储程序，实现对卫星信号的位同步、帧同步后，解析有效卫星星历信息，进行卫星导航电文 译码，最后获取卫星编号、卫星伪距等信息。通过上述信息计算出卫星的俯仰角、方位 角，并删除俯仰角小于10度的卫星。

(2)待选取有效卫星后，采用最小二乘算法进行粗略定位后，选用收敛后的接收机位置x、y、z以及接收机钟差d_t构建观测矩阵H；软件接收机上电一段时间，得到某一 时刻的观测矩阵H如下所示：

(3)从观测矩阵H可知，当前参与定位的卫星数n＝7，启动卫星故障识别算法，同时对观测矩阵H进行QR分解，得到Q、R矩阵，获取奇偶空间矩阵Q_p；

(4)构造矩阵B，使得B＝I_n-Q_p ^TQ_p；然后构造矩阵A，使得A＝H\B；

(5)计算参与定位的每颗卫星的斜率，并取其最大值Slope_max，其中

(6)给定正态分布的方差

误警报率P_FA＝1e-7，计算误警检测值σ_r和检测门限值T_D；

T_D＝Φ^-1(1-P_FA/2n)

其中，n＝7；

(7)卫星故障识别算法可用性判断：ARP＝Slope_max·σ_r≥ARP_m(7)＝343，则卫星故障识别算法可用；

(8)根据步骤(3)中的奇偶空间矩阵Q_p，计算校验统计量F_SSE＝P^TP；其中，P代表 奇偶空间矢量，该奇偶空间矢量表示卫星的偏差数据；

(9)故障检测：由于

则表明参与定位的卫星中存在故障；其中，β为给定的阈值因子；

(10)故障排除：r_imax＝max(|P^T·Q_p(:,i)|/(σ₀·norm(Q_p(:,i))))＞T_D，则说明第1颗卫 星存在故障。

(11)在剔除参与定位的故障卫星后，采用最小二乘算法重新利用剩余卫星进行定位， 输出精确的用户位置和用户钟差；

(12)采用VS2010编写MFC程序，将接收机的位置x、y、z，速度信息v_x、v_y、v_z， UTC时间，参与定位卫星的伪距、伪距率、速度、跟踪状态、多普勒计数等有用信息输 出显示给用户。

图6～图9中的实线框表示在满足高动态场景下，卫星伪距发生跳变时，使用本发明 提供的卫星故障识别方法成功将其剔除，并且给出0X42F8的标志；状态掩码0X42F8 即表示识别出此颗卫星是故障卫星。

实施例二

本实施例提供了一种卫星故障识别装置，该装置包括：

创建单元，用于采集基于最小二乘算法进行卫星定位时收敛的接收机位置数据和接 收机钟差，根据上述接收机位置数据、接收机钟差构建观测矩阵H并对该观测矩阵H进行正交三角分解，得到标准正交矩阵Q；采用奇偶空间算法对上述标准正交矩阵Q进行 分解，生成奇偶空间矩阵Q_p；

第一计算单元，用于对上述奇偶空间矩阵Q_p进行矩阵运算，得到参与定位的每个卫星的斜率并提取最大斜率值Slope_max；

第一计算单元，用于根据上述最大斜率值Slope_max和误警检测值σ_γ计算相似径向误 差保护值ARP，当该相似径向误差保护值ARP不小于预设值时，根据奇偶空间矩阵Q_p对应的卫星偏差数据计算校验统计量F_SSE；

检测单元，用于根据上述校验统计量F_SSE、误警检测值σ_γ和预置的阈值因子检测参与定位的卫星中是否存在故障。

关于卫星故障识别装置及其各功能单元的具体限定可以参见上文中对于卫星故障 识别方法的限定，在此不再赘述。上述卫星故障识别装置的各个功能单元可全部或部分通过软件、硬件及其组合来实现。上述各功能单元可以以硬件形式内嵌于或独立于计算 机装置的处理器中，也可以以软件形式存储于计算机装置的存储器中，以便于处理器调 用执行以上各个单元对应的操作。

实施例三

本实施例提供了一种软件接收机，其包括至少一个处理器、以及至少一个存储器，其中，存储器中存储有计算机程序，当计算机程序被处理器执行时，使得处理器执行上 述卫星故障识别方法的步骤。处理器和存储器的类型不作具体限制，例如：处理器可以 是微处理器、数字信息处理器、片上可编程逻辑系统等；存储器可以是易失性存储器、 非易失性存储器或者它们的组合等。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限 制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种卫星故障识别方法，其特征在于，包括以下步骤：

S1：采集基于最小二乘算法进行卫星定位时收敛的接收机位置数据和接收机钟差，根据所述接收机位置数据、接收机钟差构建观测矩阵并采用奇偶空间算法对所述观测矩阵进行分解，生成奇偶空间矩阵；

S2：对所述奇偶空间矩阵进行矩阵运算，得到参与定位的每个卫星的斜率并提取最大斜率值；

S3：基于所述最大斜率值和误警检测值计算定位卫星的相似径向误差保护值，当所述相似径向误差保护值不小于预设值时，根据所述奇偶空间矩阵对应的卫星偏差数据计算校验统计量；

S4：根据所述校验统计量、误警检测值和预置的阈值因子检测参与定位的卫星中是否存在故障：

若

则判定参与定位的卫星中存在故障；否则，判定参与定位的卫星中不存在故障；

其中，F_SSE表示校验统计量；σ_γ表示误警检测值；β表示根据接收机性能设置的阈值因子，0<β<20；n表示参与定位的卫星数量；

F_SSE＝P^TP，P＝Q_Pε；

其中，F_SSE表示校验统计量；P表示奇偶空间矢量，Q_p表示奇偶空间矩阵，ε表示伪距噪声矢量。

2.如权利要求1所述的卫星故障识别方法，其特征在于，步骤S2中具体包括：

(1)根据所述奇偶空间矩阵构造第一矩阵和第二矩阵；

B＝I_n-Q_p ^TQ_p

A＝H\B；

其中，B表示第一矩阵；I_n表示单位矩阵；A表示第二矩阵；H表示观测矩阵；n表示参与定位的卫星数量；

(2)根据所述第一矩阵和第二矩阵计算参与定位的每个卫星的斜率并取其最大斜率值Slope_max；其中，

3.如权利要求1所述的卫星故障识别方法，其特征在于，根据预置的正态分布方差、误报警率计算所述误警检测值：

其中，σ_γ表示误警检测值；

表示正态分布方差；P_FA表示误报警率；n表示参与定位的卫星数量；icdf('ch2',1-P_FA,n-4)表示求取以(n-4)为自由度，累计概率为(1-P_FA)的χ²分布的临界值。

4.如权利要求1所述的卫星故障识别方法，其特征在于，还包括排除故障卫星的步骤，具体为：

根据预置的误报警率计算检测门限值，若第i颗卫星的门限值大于检测门限值，则判定第i颗卫星存在故障并将其剔除；

T_D＝Φ^-1(1-P_FA/2n)

其中，T_D表示检测门限值；P_FA表示误报警率；Φ^-1表示正态分布的标准值；n表示参与定位的卫星数量。

5.如权利要求1所述的卫星故障识别方法，其特征在于，所述相似径向误差保护值为最大斜率值与误警检测值的乘积。

6.如权利要求1所述的卫星故障识别方法，其特征在于，所述预设值为中国地区飞机非精密进场值，当n<＝6,ARP_m(n)＝337；当n＝7,ARP_m(n)＝343；当n>＝7,ARP_m(n)＝347；

其中，ARP_m(n)表示飞机非精密进场值；n表示参与定位的卫星数量。

7.一种卫星故障识别装置，其特征在于，包括：

创建单元，用于采集基于最小二乘算法进行卫星定位时收敛的接收机位置数据和接收机钟差，根据所述接收机位置数据、接收机钟差构建观测矩阵并采用奇偶空间算法对所述观测矩阵进行分解，生成奇偶空间矩阵；

第一计算单元，用于对所述奇偶空间矩阵进行矩阵运算，得到参与定位的每个卫星的斜率并提取最大斜率值；

第一计算单元，用于根据所述最大斜率值和误警检测值计算定位卫星的相似径向误差保护值，当所述相似径向误差保护值不小于预设值时，根据所述奇偶空间矩阵对应的卫星偏差数据计算校验统计量；

检测单元，用于根据所述校验统计量、误警检测值和预置的阈值因子检测参与定位的卫星中是否存在故障：

若

则判定参与定位的卫星中存在故障；否则，判定参与定位的卫星中不存在故障；

其中，F_SSE表示校验统计量；σ_γ表示误警检测值；β表示根据接收机性能设置的阈值因子，0<β<20；n表示参与定位的卫星数量；

F_SSE＝P^TP，P＝Q_Pε；

其中，F_SSE表示校验统计量；P表示奇偶空间矢量，Q_p表示奇偶空间矩阵，ε表示伪距噪声矢量。

8.一种软件接收机，其特征在于，包括至少一个处理单元、以及至少一个存储单元，

其中，所述存储单元存储有计算机程序，当所述程序被所述处理单元执行时，使得所述处理单元执行权利要求1～6任一项所述方法的步骤。

Images