接收机自主完好性监测(RAIM)算法源于民航用户对卫星导航系统可靠性的苛刻要求^[1]。目前，卫星导航完好性监测技术主要有3种：卫星自主完好性监测(SAIM)、利用接收机内部冗余信息的RAIM和外部辅助增强系统^[2]。相对于其他2种完好性监测技术，RAIM是利用接收机的超定解对卫星测量值进行一致性校验的技术，其更能保证用户对告警时间的要求。早期的RAIM算法大多基于GPS单星座，可以满足民航航路到非精密进近阶段的完好性要求。目前，全球卫星导航系统(GNSS)发展为多星座多频点，参与定位的卫星成倍增加，接收机具有更多的冗余观测信息，可有效地减少RAIM空洞，也使RAIM技术可应用于一些对完好性性能要求更高的应用，如民航飞行中垂直引导进近(APV)Ⅰ、Ⅱ等阶段。

多星座RAIM算法是指使用2个或2个以上的卫星导航系统的冗余观测信息进行一致性校验，以检测和识别故障。单星座系统中至少观测到5颗卫星才能进行故障检测，考虑到不同星座的时间系统之间的偏差，对于多星座需要更多冗余观测量才能进行故障检测。文献[3]对RAIM的可用性预测进行了分析。文献[4-6]提出了针对多系统多故障情况下的完好性监测算法，其中文献[4]对GPS/Galileo/GLONASS组合导航在APV-Ⅰ、APV-Ⅱ阶段利用RAIM算法进行了分析。Ene等^[7-8]分析了GPS/Galileo组合导航在APV阶段的RAIM可用性。文献[9]对GPS/Galileo组合导航RAIM算法中的截止高度角、用户等效距离误差等问题进行了分析。徐君毅等^[10]对APV阶段的垂直保护水平在全球的分布情况进行了研究。Lee等^[11]提出了最优加权平均解(OWAS)算法，国内的研究人员也对此算法进行了扩展，郭靖等^[12]采用新的分组方法，满足双星座双故障的应用场景，王式太^[13]设计了基于M估计的RAIM粗差探测和剔除算法。其中，OWAS算法能够方便地兼容不同的星座，避免了不同系统间时间系统、坐标系统转换带来的误差影响，具有处理速度快、可用性高的特点，因此，本文采用OWAS算法对BDS/GPS组合导航应用于APV阶段进行了研究。

OWAS算法是将卫星按星座分为2组，监测同时发生在同一星座中的多星故障，但此算法不能有效地对双星座同时发生故障进行检测。本文采用OWAS算法，通过优化加权因子，提高算法的可用性，并对OWAS算法进行了改进。通过实测数据验证了OWAS算法对APV阶段的可用性和双星座双故障的检测。

1.   OWAS RAIM算法

OWAS RAIM算法原理如图 1所示。

图  1  OWAS RAIM算法原理

Figure  1.  OWAS RAIM algorithm principle

下载: 全尺寸图片 幻灯片

加权平均后的定位结果为

式中：加权因子0≤r≤1；x_n(n=1, 2)分别为BDS和GPS的定位结果。

图 1中：d_n(n=1, 2)分别为x_n与x_A的间隔距离。对于加权平均解，要求以0.999的概率满足：|x_A－x|≤|x_A－x_n|+|x_n－x|≤T+a_n。其中：T(n=1, 2)为检验门限，a_n为无故障假设条件下定位误差上限。所以，无故障卫星时，每个导航系统内的保护水平为：PL_n=T+a_n，PL=max(PL₁, PL₂)。OWAS算法在只考虑垂直方向上的定位精度与完好性要求时，可以归纳为如下问题：

式中：VPL为垂直保护级别。

同时，满足以下条件时，RAIM算法可用：

式中：σ_V为垂直方向上的定位误差标准差；σ_V-max为垂直方向上的定位误差标准差的最大值；d_major为最终定位解水平方向上的定位误差标准差；d_major-max为飞行阶段规定的95%水平方向上的定位误差标准差的最大值；VAL为垂直保护门限；HPL为水平保护级别；HAL为水平保护门限^[10]。

σ_V和d_major的标准由以下2个操作需求推导：

1) 95%的导航系统误差(Navigation System Errors，NSE)边界：1.96σ_V=4m，1.96d_major=16m(1.96对应于高斯分布假设的95%分布)。

2) NSE边界满足无故障条件完好性可用性时：K_Vσ_V≤VAL(垂直方向)，K_Hd_major≤HAL(水平方向)，当P_r{HMI}=10^－8/approach(P_r为概率, HMI为危险误导信息, approach为进近)时，K_V=5.73，K_H=6.47^[10]。

σ_V-max和d_major-max的值分别为

1.1   协方差矩阵

OWAS算法为

式中：w₁为GPS定位解在最终解中所占权重；w₂为BDS定位解在最终解中所占权重；，且0≤r≤1。

OWAS算法的位置误差与加权因子的关系如下：

式中：Δx₁为GPS定位解与真实位置之间的误差；Δx₂为BDS定位解与真实位置之间的误差；Δx_A为最终加权平均定位解与真实位置之间的误差。

加权最小二乘定位解的位置误差^[14]表达式如下：

式中：G为伪距误差与接收机位置及钟差的关联矩阵；W为加权矩阵；ΔR为伪距误差矢量；P_n为定位误差的协方差矩阵(n=1, 2分别表示GPS和BDS)。

最优解和各导航系统定位解误差互相关表达式如下：

假设GPS与BDS伪距误差不相关，则E(ΔR₂·ΔR₁^T)=0，E(ΔR₁ΔR₂^T)=0。由此，得

因此，E(Δx_AΔx₁^T)=rP₁，E(Δx_AΔx₂^T)=(1-r)P₂，即

同理，加权平均解x_A的定位误差协方差矩阵P_A如下：

协方差矩阵J_n为

将式(10)、式(11)、式(15)、式(17)代入式(18)中，可得

式中：P_S=P₁+P₂。

同理

1.1.1   加权因子确定保护级别

令

式中：Q^－1为标准单边累加分布函数的逆函数；P_fd为虚警率(参考国际民航不同进近阶段导航性能要求，虚警率与连续性保持一致。虚警率除以4是因为总的允许虚警率在垂直和水平方向上平均分配^[10])；P_md为最大允许漏警率(完好性风险等于漏警率与卫星故障率的乘积)。

令

由式(19)和式(20)可以得出

VPL₁和VPL₂分别为

由图 2可以看出，VPL₁是加权因子的一元线性递减函数；VPL₂是加权因子的一元线性递增函数。在OWAS算法中，VPL=max(VPL₁, VPL₂)，当VPL₁=VPL₂时，得到最小化的VPL，此时可得

图  2  VPL₁、VPL₂与加权因子的函数

Figure  2.  VPL₁ and VPL₂ as a function of weight factor

下载: 全尺寸图片 幻灯片

1.1.2   加权因子确定误差精度σ_V和d_major

OWAS位置误差的协方差矩阵P_A=E(Δx_AΔx_A^T)=r²P₁+(1-r)²P₂。

σ_V和d_major为

由式(3)、式(4)可知，RAIM算法可用要满足：

1.1.3   确定最优加权因子

OWAS最重要的是在降低保护水平的条件下满足精度的要求。

如图 3所示，r_min-VPL表示垂直保护水平最小时的加权因子，但此时不满足σ_V的标准，而r_min-σV满足最小化σ_V，却不满足VPL≤VAL，r_best-int表示在保证精度要求下修正加权因子r_min-VPL来最小化VPL，r_best-acc表示通过修正最小化精度加权因子σ_V保证可用性的要求。

图  3  加权因子的计算

Figure  3.  Determination of weight factor

下载: 全尺寸图片 幻灯片

APV飞行阶段时需更大可能的保证算法的可用性，故采用r_best-int作为最终加权因子，则

式(23)和式(24)代入到式(34)中，可得

由式(35)计算可得2个解，如下：

由于r₁≤r≤r₂，且0≤r≤1，所以

式中：r_{acc, min}=max(r₁, 0)；r_{acc, max}=min(r₂, 1)。

当r_{acc, min}≤r_min-VPL≤r_{acc, max}时，加权因子r如下：

当r_min-VPL≤r_{acc, min}时，加权因子r如下：

当r_{acc, max}≤r_min-VPL时，加权因子r如下：

1.2   故障检测与识别

OWAS RAIM算法是以GPS、BDS到x_A的距离作为故障检测的检验统计量^[15-16]，按照国际民航组织(ICAO)要求的性能标准得到检验门限。当双星座同时出现故障时，2个检验统计量均会在故障产生的时刻超过检验门限。当单星座产生故障时，故障星座的定位解使得OWAS与无故障星座的距离增加，容易产生误警，此时采用加权最小二乘残差的故障检测对OWAS定位域的故障检测进行辅助，完成单星座故障的检测。

1.2.1   故障检测

GNSS伪距观测线性化方程为

式中：ε为观测噪声矢量。

对式(42)运用加权最小二乘法求解，得到的解为

式中：B=(G^TWG)^－1G^TW。

伪距残差矢量v为

矩阵C的协因数矩阵为

伪距残差平方和为

在卫星导航系统正常工作的情况下，各卫星伪距残差都较小，伪距残差平方和也较小；当某个测量伪距存在较大偏差时，伪距残差平方和会变大。假设卫星导航系统正常工作时，观测噪声ε中的各个分量是均值为0、方差为σ²的正态分布随机误差，且它们相互独立。依据统计分布理论，若ε的均值为0，SSE_WLS服从自由度为(n－4)的χ²分布；若ε的均值不为0，则SSE_WLS服从自由度为(n－4)的非中心化χ²分布。对SSE_WLS作二元假设：

1) 无故障假设H₀:E(ε)=0，则SSE_WLS~χ²(n－4)。

2) 有故障假设H₁：E(ε)≠0，则SSE_WLS~χ²(λ, n－4)。

由式(46)可以解得t_WLS值。

1.2.2   故障识别

根据残差和观测误差的关系式，定义故障识别的检验统计量为

式中：i=1, 2, …, n，并且n为可见卫星数目；v_i为利用加权最小二乘法求得的伪距残差；Q_ii为矩阵Q的第i行第i列；d_i~N(0, 1)。

对统计量d_i做二元假设：

1) H₀(无故障):E(ε_i)=0，则d_i~N(0, 1)。

2) H₁(有故障):E(ε_i)≠0，则d_i~N(δ_i, 1)。

其中，δ_i为统计量偏移参数。给定总体的误警率P_FA，可计算出对应的检验门限T_1WLS。计算如下：

由式(48)可解得门限值T_1WLS。

设为检验统计量，对应的检验门限值为T_1WLS。比较检验统计量d_i与检验门限T_1WLS进行卫星故障识别。

2.   算法验证与结果分析

采用IGS(International GPS Service)2016年7月31日0时的BDS、GPS系统RINEX数据对本文算法进行仿真验证。

在GPS可见星中的第16颗卫星、BDS可见星中的第7颗卫星中，分别在500~900s时刻人为加入30、50、70m伪距偏差。

图 4(a)、(b)分别为未加伪距偏差时GPS和BDS的垂向定位域检验统计量。图 5(a)、(b)分别为在500 s开始加入30、50、70 m伪距偏差时GPS、BDS的检验统计量。可知，当加入伪距偏差时，检验统计量值会发生跳变; 未加入伪距偏差时，检验统计量值波动较小。

图  4  BDS/GPS未加入伪距偏差的检验统计量

Figure  4.  BDS/GPS test statistics with no pseudo-range bias

下载: 全尺寸图片 幻灯片

图  5  BDS/GPS加入不同伪距偏差的检验统计量

Figure  5.  BDS/GPS test statistics with different pseudo-range bias

下载: 全尺寸图片 幻灯片

双故障卫星检测结果如图 6所示。

图  6  BDS/GPS故障卫星检测

Figure  6.  BDS fault satellite detection

下载: 全尺寸图片 幻灯片

由图 6可知，在GPS和BDS系统的第16颗卫星和第7颗卫星中加入50m伪距偏差，检验统计量大于检验门限，可以对故障进行检测。

双星座双故障算法识别率结果如图 7所示。可知，当伪距偏差值较小时，不能识别出故障卫星；当故障卫星伪距偏差在30m以上时，BDS/GPS能识别全部故障。

图  7  故障识别率

Figure  7.  Fault identification rate

下载: 全尺寸图片 幻灯片

3.   结论

本文研究了一种BDS/GPS组合导航RAIM算法，对算法进行了详细的理论分析，并利用IGS网站数据对算法进行验证。

1) 本文算法对BDS/GPS双卫星故障检测与识别是可用的。

2) 当伪距偏差在30m以上时可有效地检测与识别。

研究结果对多星座组合导航接收机自主完好性监测算法研究有一定的参考价值。