CN109031356A - 特征斜率加权的最小二乘残差接收机自主完好性监测方法 - Google Patents

Abstract

本发明公开的特征斜率加权的最小二乘残差接收机自主完好性监测方法，属于卫星导航技术领域。本发明实现方法为，构造新的检验统计量，引入与各可见星特征斜率相关的设计权阵，对伪距的最小二乘残差矢量的各分量进行加权，以调节各可见星残差矢量在检验统计量中所占的比重，使检验统计量更准确的反映伪距偏差对定位误差的影响，使检验统计量更准确的反映伪距偏差对定位误差的影响，降低以最小二乘残差LSR方法完成接收机自主完好性监测RAIM时的漏检率和虚警率，为用户提供更为准确的完好性信息。本发明适用于GNSS的接收机自主完好性监测，应用于卫星导航技术领域，提高导航精度。

Description

特征斜率加权的最小二乘残差接收机自主完好性监测方法

技术领域

本发明涉及一种特征斜率加权的最小二乘接收机自主完好性监测方法，属于卫星导航技术领域。

背景技术

接收机自主完好性监测(Receiver Autonomous Integrity Monitoring，RAIM)方法是利用全球导航卫星系统(GNSS)观测信息的一致性实现GNSS完好性监测，对故障星进行检测、识别和隔离的方法。RAIM方法分为两种，快照法和序贯法，前者仅采用当前历元的观测值进行完好性监测，而后者需要联合多个过往历元的观测值。快照法由于计算简单、使用方便而被广泛使用。

快照法包括最小二乘残差(least squares residuals，LSR)法和奇偶矢量法两种，二者具有等效性。最小二乘残差法的原理是利用最小二乘法根据冗余观测量检测和识别卫星故障。选取观测量为伪距，以LSR实现RAIM的具体过程如下：

步骤一、存在观测偏差情况下，建立GNSS伪距观测方程

式[1]中，ρⁱ为第i颗可见星(记作VS_i)(i＝1,…,K,K为可见星总数)和接收机间的伪距；Pⁱ＝(xⁱ,yⁱ,zⁱ)和P^u＝(x^u,y^u,z^u)分别为地心地固坐标(ECEF) 系下VS_i和接收机的三维位置；c为光速；Δt为接收机钟差；ε为观测误差，包括观测噪声和观测偏差。

步骤二、对伪距观测方程进行线性化处理(泰勒展开)

式[2]中，ρⁱ和Rⁱ分别为VS_i和接收机之间的伪距和真距，为接收机的近似位置。

步骤三、联立各可见星与接收机之间的线性化伪距方程得到

Z＝HX+ε [3]

式[3]中，观测矩阵

其中前3个元素是到P^u的三维位置改正数；Z＝[ρ¹-R¹ … ρⁱ-Rⁱ … ρ^K-R^K]^T；ε＝[ε¹ … εⁱ … ε^K]^T，其中εⁱ表示VS_i的观测误差。

步骤四、根据最小二乘定位原理，得到X的最小二乘解为

P^u的最小二乘解为

设结合式[1]和式[5]可知伪距残差矢量ω为

其中是根据式[5]和式[10]由牛顿最小二乘法迭代法计算的。每次迭代被重置为直到为止，s取很小的正数，如10-⁸(的初值设定为0)。在成立的条件下，

步骤五、对ω进行简化，将式[8]代入式[7]得到

步骤六、得到伪距残差ω的平方和

SSE＝ω^Tω [10]

在伪距观测噪声为高斯白噪声(white Gaussian noise，WGN)的假设下，若 GNSS无故障，SSE/σ₀ ²服从自由度为n-4的χ²分布，记作SSE/σ₀ ²～χ²(n-4)；若GNSS有故障，则SSE/σ₀ ²～χ²(n-4,λ)，λ为非中心化参数，σ₀为观测噪声的标准差。

步骤七、构造故障检测的检验统计量Ts

Ts＝SSE/σ₀ ² [11]

以H0代表GNSS无故障，H1代表GNSS有故障，则LSR的故障检测判别式为

式[12]中，Td为根据虚警率(Probability of False Alarm，P_FA)计算的故障检测门限，根据GNSS完好性风险需求，取P_FA＝10^-6。

故障检测的检验统计量Ts由各可见星伪距残差矢量的平方和构造，而故障识别的检验统计量Ts_i由可见星VS_i的伪距残差构造：

式[13]中ωⁱ为式[9]中ω的第i个分量；Q_ii为矩阵Q的第i各对角线元素，Q 为伪距残差矢量的协因数矩阵

Q＝I-H(H^TH)^-1H^T [14]

在伪距观测噪声为WGN的假设下，如果VS_i无故障，则Ts_i服从标准正态分布，记作Ts_i～N(0,1)；反之，Ts_i～N(δ_i,1)。以H0代表VS_i无故障，H1代表VS_i有故障，则针对VS_i的故障识别判别式为

式[15]中，Td_i为根据P_FA和可见星总数计算的故障识别门限。

根据式[11]和式[13]所构造的伪距残差检验统计量其本质是伪距观测偏差的投影参量，该参量与定位误差为正交关系。图1给出伪距观测偏差在残差域和垂直定位误差域之间的关系。

图1中的横轴为伪距残差域，纵轴为垂直定位误差域，垂直告警限值和检测门限Td将该区域分为正常定位、检出故障、漏检和虚警四部分。垂直告警限值即需告警的垂直定位误差门限。图中斜线向量为伪距偏差向量，其斜率即为伪偏差向量的特征斜率。

当伪距偏差向量出现在左下角时(如图中的A)，定位误差和残差检验统计量均不超过相应的门限，可正常定位；当伪距偏差向量出现在右上角时(如图中的B)，定位误差和残差检验统计量都超过相应门限，此时可正常告警即，即检测并识别故障星；当伪距偏差向量出现在左上角区域时(如图中的C)，残差检验统计量未超出检出门限，但定位误差已超过告警限值，此时的GNSS定位结果是不可用的；当伪距偏差向量出现在右下角区域时(如图中的D)，残差检验统计量已超过检测门限，但定位误差未超过告警限值，此时发生虚警，事实上，此时的GNSS定位结果是可用的。上述四种结果中，正常定位是最理想的，检出故障并排除也可接受，但漏检和虚警都是RAIM监测的异常情况，是设计 RAIM方法时需尽力规避的。

对于每颗卫星，其伪距偏差向量在残差域和垂直定位误差域都存在相应的投影，根据二者关系，伪距偏差向量的特征斜率为

式[16]中，其中，下标3i第3 行第i列的元素，ii表示第i行第i列的元素，H_s是观测矩阵H在站心坐标系中的表示。

根据图1可发现，在伪距偏差向量值一定时，可见星的特征斜率越大，伪距偏差造成的定位误差越大，残差检验统计量反而越小，反之，特征斜率越小，定位误差越小，残差检验统计量反而越大。

发明内容

本发明公开的特征斜率加权的最小二乘残差接收机自主完好性监测方法要解决的技术问题是：通过特征斜率加权的最小二乘残差法实现接收机自主完好性监测，降低以最小二乘残差LSR方法完成接收机自主完好性监测RAIM时的漏检率和虚警率。

本发明的目的是通过下述技术方案实现的。

本发明公开的特征斜率加权的最小二乘残差接收机自主完好性监测方法，实现方法为，构造新的检验统计量，引入与各可见星特征斜率相关的设计权阵，对伪距的最小二乘残差矢量的各分量进行加权，以调节各可见星残差矢量在检验统计量中所占的比重，使检验统计量更准确的反映伪距偏差对定位误差的影响，使检验统计量更准确的反映伪距偏差对定位误差的影响，降低以最小二乘残差LSR方法完成接收机自主完好性监测RAIM时的漏检率和虚警率。为用户提供更为准确的完好性信息。本发明适用于GNSS的接收机自主完好性监测，应用于卫星导航技术领域，提高导航精度。

本发明公开的特征斜率加权的最小二乘残差接收机自主完好性监测方法，又称Slope-LSR方法，包括如下步骤：

步骤一：在某观测历元某地点，共有K颗GNSS可见星。可见星与接收机之间的线性化伪距方程为：

Z＝HX+ε [1]

其中，地心地固ECEF坐标系下的观测矩阵

前3个元素是到P^u的三维位置改正数；Z＝[ρ¹-R¹ … ρⁱ-Rⁱ … ρ^K-R^K]^T；ε＝[ε¹ … εⁱ … ε^K]^T，其中：ρⁱ和Rⁱ分别为第i颗可见星和接收机之间的伪距和真距，为接收机的近似位置， Pⁱ＝(xⁱ,yⁱ,zⁱ)和P^u＝(x^u，y^u，z^u)分别为地心地固ECEF坐标系下第i颗可见星和接收机的三维位置；c为光速；Δt为接收机钟差；ε为观测误差，包括观测噪声和观测偏差。

根据最小二乘定位原理，X的最小二乘解为

P^u的最小二乘解为

设可知伪距最小二乘残差矢量ω为：

步骤二：计算站心坐标系中的观测矩阵H_s

H_s是观测矩阵H在站心坐标系中的表示，即

H_s＝HR [6]

其中，R为可见星与接收机之间的真距。

步骤三：计算各可见星的伪距偏差特征斜率SlopeV_i，其中：i＝1,…,K。

对于每颗GNSS可见星，其伪距偏差向量在残差域和垂直定位误差域都存在相应的投影，根据二者关系，伪距偏差向量的特征斜率为

其中，其中，下标3i第3行第 i列的元素，ii表示第i行第i列的元素。

在伪距偏差向量值一定时，可见星的特征斜率越大，伪距偏差造成的定位误差越大，残差检验统计量反而越小，反之，特征斜率越小，定位误差越小，残差检验统计量反而越大。

步骤四：根据各可见星的特征斜率设计权阵W。

其中，w_i是与SlopeV_i有关的设计权值，

式中，函数med(·)代表取中位数，T_Slope是伪距偏差向量在残差域和垂直定位误差域投影分别等于故障检测门限和垂直告警限值时的特征斜率，即

T_Slope＝HAL/Td [10]

式中，HAL和Td分别代表垂直告警限值和故障检测门限。

步骤五：构造故障检测的检验统计量Ts_w。

步骤六：利用构造故障检测的检验统计量Ts_w进行故障检测。

Ts_w＜Td，不存在故障星；

Ts_w≥Td，存在故障星。

若Ts＜Td，则转入步骤九；Ts≥Td，则转入步骤七。

步骤七：构造故障识别的检验统计量

步骤八：利用构造故障识别的检验统计量识别故障星星号。

VSi无故障；

VSi有故障。

其中，下标i即为对应的故障星星号。

步骤九：根据步骤六故障检测和步骤八的故障识别结果，输出接收机自主完好性监测RAIM监测结果。

根据步骤六进行故障检测，如果未检测到故障，则监测结果为GNSS不存在故障星；如果监测到故障，根据步骤八识别故障星星号，则监测结果为GNSS 存在故障星，并输出故障星的星号，即实现输出接收机自主完好性监测RAIM 监测结果。

有益效果：

1、本发明公开的特征斜率加权的最小二乘残差接收机自主完好性监测方法，通过特征斜率加权的最小二乘残差法实现接收机自主完好性监测，即通过步骤五构造故障检测的检验统计量Ts_w，该统计量中引入与各可见星特征斜率相关的设计权阵，以调节各可见星残差矢量在检验统计量中所占的比重，使检验统计量更准确的反映伪距偏差对定位误差的影响，降低以最小二乘残差LSR方法完成接收机自主完好性监测RAIM时的漏检率和虚警率。

2、本发明公开的特征斜率加权的最小二乘残差接收机自主完好性监测方法，利用构造故障识别的检验统计量识别故障星星号，能够在发现可见星故障的情况下准确的输出对应的故障星星号，便于使用过程中的快速修复及解决，进一步提高了RAIM结果的可靠性。

附图说明

图1为伪距观测偏差投影图；

图2为本发明特征斜率加权的最小二乘残差RAIM方法流程图；

图3为对G17可见星注入45m观测偏差的仿真结果，其中，图3a为检验统计量结果，图3b为垂直定位误差结果；

图4为对G09可见星注入45m观测偏差的仿真结果,其中，图4a为检验统计量结果，图4b为垂直定位误差结果。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。

下面以具体数值进一步说明本发明的有益效果。

本实施例在数学仿真环境中进行，设定仿真开始时间为UTC时间2015-5-2 12:00:00，仿真步长为1s，仿真持续时间为400s，即共有400个观测历元。接收机初始位置为纬度39°，经度116°，高度0m，在ECEF坐标系中的位置是 [-2175778.59,4461007.21,3992315.37]m。伪距观测量是根据GPS的RINEX星历和接收机的真实位置计算所得的真实距离加上标准差为5m的WGN，符合实际情形中排除电离层、对流层和多径误差后伪距观测量的真实情况。GPS卫星是否可见的高度判别角为15°。在第50～150历元，分别在可见星G09和G17的伪距观测量中注入45m(9σ₀)观测偏差，以Slope-LSR方法实现RAIM监测，具体实施步骤如下，分别以不加伪距偏差和加入伪距观测偏差元给出具体说明：

实施例1：

不加伪距观测偏差

步骤1：计算伪距最小二乘残差矢量ω和ECEF坐标系下的观测矩阵H。在第1个观测历元共有10颗可见星，各可见星的位置和伪距参见表1。

表1.第1个历元GPS可见星的三维位置和伪距

结合表1所示的卫星位置、伪距以及接收机的位置，计算伪距最小二乘残差矢量ω和观测矩阵H。

ω＝[0.54,1.63,-6.13,10.87,0.18,3.46,0.13,-2.82,-8.55,0.69]^T

步骤2：计算站心坐标系中的观测矩阵H_s＝HR。

将纬度39°，经度116°代入公式计算可得：

步骤3：计算各可见星的伪距偏差特征斜率SlopeV_i(i＝1,…,K)。

根据观测矩阵H_s可计算各可见星的特征斜率，在第1个观测历元，各可见星的特征斜率值参见表2。

步骤4：根据各可见星的特征斜率设计权阵W。

在第1个观测历元，特征斜率设计权阵W的各对角线元素参见表2。

表2.第1个历元GPS可见星的特征斜率值

序号	PRN	slopeVi	wi
				1	G01	0.86	1.70
2	G03	0.52	1.03
				3	G07	0.41	0.82
4	G08	0.30	0.60
				5	G09	0.10	0.20
6	G11	0.30	0.60
				7	G17	1.29	2.55
8	G20	0.22	0.43
				9	G28	0.62	1.23
10	G31	0.79	1.57

步骤5：构造故障检测的检验统计量Ts_w。

第1个观测历元计算的检验统计量Ts_w＝8.30；检验统计量取χ²分布的分位数(自由度为6)，即Td＝α，P_α＝10-⁶，其值为Td＝38.26

步骤6：检测故障。

结果Ts_w＜Td，转入步骤9。

步骤9：输出RAIM监测结果：GPS无故障。

实施例2：

加入伪距观测偏差

以在G17中注入45m观测偏差为例，给出具体执行步骤。

步骤1：计算伪距最小二乘残差矢量ω和ECEF坐标系下的观测矩阵H。 (H与第1观测历元中的H大致相等，此处不再给出具体数值)

ω＝[6.06,5.24,0.00,-4.36,-17.94,8.54,18.59,11.46,-6.26,-21.33]^T

步骤2：计算站心坐标系中的观测矩阵H_s＝HR。(H_s与第1观测历元中的 H_s大致相等，此处不再给出具体数值)

步骤3：计算各可见星的伪距偏差特征斜率SlopeV_i(i＝1,…,K)。

根据观测矩阵H_s可计算各可见星的特征斜率，在第50个观测历元，各可见星的特征斜率值参见表3。

步骤4：根据各可见星的特征斜率设计权阵W。

在第55个观测历元，特征斜率设计权阵W的各对角线元素参见表3。

表3.第50个历元GPS可见星的特征斜率值和设计权值

序号	PRN	slopeVi	wi
				1	G01	0.85	1.71
2	G03	0.53	1.06
				3	G07	0.42	0.84
4	G08	0.30	0.60
				5	G09	0.10	0.20
6	G11	0.30	0.60
				7	G17	1.27	2.56
8	G20	0.23	0.46
				9	G28	0.63	1.27
10	G31	0.80	1.60

步骤5：构造故障检测的检验统计量Ts_w。

第55个观测历元计算的检验统计量Ts_w＝43.98；检验统计量Td＝38.26。

步骤6：检测故障。

Ts_w≥Td，转入步骤7。

步骤7：构造故障识别的检验统计量

计算各可见星的故障识别检验统计量，在第55个观测历元，各可见星的参见表4。

表5.第55个历元各GPS可见星的故障识别检验统计量

故障识别门限Td_i由可见星总数和正态分布规律共同确定。10颗可见卫星共得到10个故障识别检验统计量，给定系统虚警概率为P_α＝10^-6，则每个故障识别检验统计量的虚警概率为Pα/K。Td_i为标准正态分布分位数，即Td＝α， P_α＝10-⁶/2K。在第150个观测历元，共有10颗可见星，Td_i＝5.33。

步骤8：识别故障。

分别比较表5中各可见星的和Td_i可知，G17的检验统计量超过门限，其余可见星的未超过门限，故G17为故障星，其余为无故障星。

步骤9：输出RAIM监测结果。

监测结果：GPS存在故障星，故障星为G17。

本实施例中，以现有的最小二乘残差(LSR)RAIM方法为对比，证明 Slope-LSR方法可降低漏检率和虚警率。仿真结果如图3和图4所示。

图3和图4是在第50-150个观测历元分别对G17和G09两颗可见星注入 45m伪距偏差所得的Slope-LSR仿真结果，包括每个观测历元的检验统计量和垂直定位误差，其中检验统计量以LSR做对比。

如图3所示，由于G17的特征斜率大(见表3)，对其注入45m观测偏差后垂向定位误差明显增加。在第50-150历元，Slope-LSR计算所得检验统计量超出故障检测门限的个数(42个)明显多于LSR(19个)，故Slope-LSR与LSR 相比漏检率明显下降。

如图3所示，由于G09的特征斜率大(见表3)，对其注入45m观测偏差后垂向定位误差没有明显变化。在第50-150历元，Slope-LSR计算所得检验统计量超出故障检测门限的个数(33个)明显少于LSR(99个)，故Slope-LSR与 LSR相比虚警率明显下降。

综上所述，与传统LSR相比，Slope-LSR可有效降低漏检率和虚警率。

本发明的Slope-LSR法对传统LSR法进行了改进，降低了利用最小二乘伪距残差实现RAIM时的漏检率和虚警率。将本方法运用到实际导航过程中，主要是航空领域，可降低危险误导信息和虚警信息的出现概率，为用户提供更为准确的完好性信息。本发明所提供的方法适用于GNSS的接收机自主完好性监测。

需要特别说明的是，由于航空领域对垂向导航服务性能要求较高，本发明的设计权阵中选用的特征斜率为垂向特征斜率。对于其它导航应用场景，如果对水平方向导航服务性能要求较高，可将设计权阵中的特征斜率选为水平特征斜率。

Claims (2)

1.特征斜率加权的最小二乘残差接收机自主完好性监测方法，又称Slope-LSR方法，其特征在于：包括如下步骤，

步骤一：在某观测历元某地点，共有K颗GNSS可见星；可见星与接收机之间的线性化伪距方程为：

Z＝HX+ε [1]

其中，地心地固ECEF坐标系下的观测矩阵

前3个元素是到P^u的三维位置改正数；Z＝[ρ¹-R¹ … ρⁱ-Rⁱ … ρ^K-R^K]^T；ε＝[ε¹ … εⁱ … ε^K]^T，其中：ρⁱ和Rⁱ分别为第i颗可见星和接收机之间的伪距和真距，为接收机的近似位置，Pⁱ＝(xⁱ,yⁱ,zⁱ)和P^u＝(x^u,y^u,z^u)分别为地心地固ECEF坐标系下第i颗可见星和接收机的三维位置；c为光速；Δt为接收机钟差；ε为观测误差，包括观测噪声和观测偏差；

根据最小二乘定位原理，X的最小二乘解为

P^u的最小二乘解为

设可知伪距最小二乘残差矢量ω为：

步骤二：计算站心坐标系中的观测矩阵H_s

H_s是观测矩阵H在站心坐标系中的表示，即

H_s＝HR [6]

其中，R为可见星与接收机之间的真距；

步骤三：计算各可见星的伪距偏差特征斜率SlopeV_i，其中：i＝1,…,K；

对于每颗GNSS可见星，其伪距偏差向量在残差域和垂直定位误差域都存在相应的投影，根据二者关系，伪距偏差向量的特征斜率为

其中，其中，下标3i第3行第i列的元素，ii表示第i行第i列的元素；

在伪距偏差向量值一定时，可见星的特征斜率越大，伪距偏差造成的定位误差越大，残差检验统计量反而越小，反之，特征斜率越小，定位误差越小，残差检验统计量反而越大；

步骤四：根据各可见星的特征斜率设计权阵W；

其中，w_i是与SlopeV_i有关的设计权值，

式中，函数代表取中位数，T_Slope是伪距偏差向量在残差域和垂直定位误差域投影分别等于故障检测门限和垂直告警限值时的特征斜率，即

T_Slope＝HAL/Td [10]

式中，HAL和Td分别代表垂直告警限值和故障检测门限；

步骤五：构造故障检测的检验统计量Ts_w；

步骤六：利用构造故障检测的检验统计量Ts_w进行故障检测；

Ts_w＜Td，不存在故障星；

Ts_w≥Td，存在故障星；

若Ts＜Td，则转入步骤九；Ts≥Td，则转入步骤七；

步骤七：构造故障识别的检验统计量

步骤八：利用构造故障识别的检验统计量识别故障星星号；

VSi无故障；

VSi有故障；

其中，下标i即为对应的故障星星号；

步骤九：根据步骤六故障检测和步骤八的故障识别结果，输出接收机自主完好性监测RAIM监测结果；

根据步骤六进行故障检测，如果未检测到故障，则监测结果为GNSS不存在故障星；如果监测到故障，根据步骤八识别故障星星号，则监测结果为GNSS存在故障星，并输出故障星的星号，即实现输出接收机自主完好性监测RAIM监测结果。

2.如权利要求1所述的特征斜率加权的最小二乘残差接收机自主完好性监测方法，其特征在于：应用于卫星导航技术领域，提高导航精度。