BDS是我国自主建设、独立运行的全球卫星导航系统^[1-4]，可为用户提供更多高定位精度和可靠性的导航信息，同时也可为民用航空接收机自主完好性监测(receiver autonomous integrity monitoring，RAIM)提供有力支撑。RAIM对卫星进行故障监测会受到卫星个数和卫星几何分布的影响。当可见卫星数过少或几何分布不理想时，由卫星几何结构引起的定位误差会遮掩卫星故障引起的定位误差，造成完好性监测结果不可信，漏检率增加。因此在进行故障监测之前必须进行可用性判断，以保证不会影响故障监测性能^[5]。RAIM可用性计算中涉及到漏检率和误警率参数，王尔申等^[6]研究发现，RAIM可用性随漏检率和误警率的降低而降低。

目前，基于单星故障的RAIM可用性评估方法主要有水平保护法(horizontal protection level method，HPLM)、最大水平精度因子法(δHDOP_max)和近似径向误差保护法(approximate radial-error protection method，ARPM)，理论上这3种方法具有等价性^[7-8]；考虑到2颗卫星同时发生故障的可能性，倪育德等^[9]从1颗卫星故障的可用性出发，推导出2颗卫星故障的可用性；针对多故障卫星的可用性，陈金平等^[10]提出基于漏检概率和圆概率误差的可用性分析方法。随着GNSS的发展与应用以及飞机更高阶段精密进近的需求，具备垂直导航能力的高级接收机自主完好性监测(advanced RAIM，ARAIM)被引入。ARAIM旨在全球范围内提供垂直导向-200 ft定标性能(LPV-200)级别的航空导航服务，众多学者对此进行了相关算法改进和定位性能方面研究^[11-17]。

现有的双星故障条件下RAIM可用性评估大多沿用已有的矩阵最大特征值方法(matrix maximum eigenvalue method，MMEM)^[18]，该方法可解决2个粗差比值的计算，但该值还可采用高等数学中求极值方法得到。为此，本文提出针对RAIM可用性评估的极大值方法(maxima method，MM)，并基于完好性风险参数和IGMAS中8个GNSS站的实测数据对这2种方法的RAIM可用性进行验证，得到最新的中国境内BDS可用性性能。

1 RAIM方法原理 1.1 RAIM观测方程建立

RAIM主要是基于GNSS中的伪距观测量。若在历元t时接收机r与卫星s进行同步观测，并观测m颗卫星，则观测方程可表示为：

$ \mathit{\boldsymbol{y}}=\mathit{\boldsymbol{H x}}+\mathit{\boldsymbol{\varepsilon}} $

(1)

式中，y为m维伪距观测值与近似卫地距的差向量；H为m行4列设计矩阵，表示各卫星对用户的投影向量；x为测站坐标和接收机钟差构成的四维向量，ε为m维观测伪距的误差向量。

给定观测量权阵为P，则式(1)的最小二乘解为：

$ \begin{gathered} \hat{\boldsymbol{x}}=\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{y}=\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{P}(\boldsymbol{H} \boldsymbol{x}+\boldsymbol{\varepsilon}) \\ =\boldsymbol{x}+\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{\varepsilon} \end{gathered} $

(2)

观测量的改正数向量($ \hat{\boldsymbol{v}}$) 为：

$ \begin{gathered} \hat{\boldsymbol{v}}=\boldsymbol{y}-\boldsymbol{H} \hat{\boldsymbol{x}}=\boldsymbol{H} \boldsymbol{x}+\boldsymbol{\varepsilon}-\boldsymbol{H}\left(\boldsymbol{x}+\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{P H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{\varepsilon}\right)\\ =\left(\boldsymbol{I}-\boldsymbol{H}\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{P}\right) \boldsymbol{\varepsilon} \end{gathered} $

(3)

令

$ \left\{\begin{array}{l} \boldsymbol{Q}=\boldsymbol{P}^{-1}-\boldsymbol{H}\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \\ \boldsymbol{A}=\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \\ \boldsymbol{B}=\boldsymbol{H}\left(\boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{H}\right)^{-1} \boldsymbol{H}^{\mathrm{T}} \boldsymbol{P} \end{array}\right. $

(4)

则定位误差(Δx)及$ \hat{\boldsymbol{v}}$的加权平方和($ \hat{\boldsymbol{v}}^{\mathrm{T}} {\boldsymbol{P}}\hat{\boldsymbol{v}}$) 为：

$ \left\{\begin{array}{l} \Delta \boldsymbol{x}=\hat{\boldsymbol{x}}-\boldsymbol{x}=\boldsymbol{A} \boldsymbol{\varepsilon} \\ \hat{\boldsymbol{v}}=(\boldsymbol{I}-\boldsymbol{B}) \boldsymbol{\varepsilon} \\ \hat{\boldsymbol{v}}^{\mathrm{T}} {\boldsymbol{P}}\hat{\boldsymbol{v}}=\boldsymbol{\varepsilon}^{\mathrm{T}} \boldsymbol{P}(\boldsymbol{I}-\boldsymbol{B}) \boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}^{\mathrm{T}} \boldsymbol{P Q P \boldsymbol { \varepsilon }} \end{array}\right. $

(5)

1.2 检验统计量构建

假定ε服从0均值整体正态随机分布。设分量方差为σ₀²，依据统计理论，$ {\mathit{\boldsymbol{\hat v}}^{\rm{T}}}\mathit{\boldsymbol{P\hat v}}$/σ₀²服从自由度为(m－4)的卡方分布χ²(m－4)；当观测量中存在粗差时，$ {\mathit{\boldsymbol{\hat v}}^{\rm{T}}}\mathit{\boldsymbol{P\hat v}}$/σ₀²服从自由度为(m－4)的非中心卡方分布χ²(m－4, λ)，其中λ为$ {\mathit{\boldsymbol{\hat v}}^{\rm{T}}}\mathit{\boldsymbol{P\hat v}}$/σ₀²分布的非中心参数。根据$ {\mathit{\boldsymbol{\hat v}}^{\rm{T}}}\mathit{\boldsymbol{P\hat v}}$/σ₀²的分布特征可进行故障卫星探测。

无故障发生时，概率密度用f_χm－4²(x)表示。当系统处于正常状态时，无伪距故障发生，如果系统发出监测报警，则为虚警(false alarm，FA)，虚警率(P_FA)可表示为P_FA=P($ \frac{{{{\mathit{\boldsymbol{\hat v}}}^{\rm{T}}}\mathit{\boldsymbol{P\hat v}}}}{{\sigma _0^2}}$＞α²|无故障)，α²为χ²分布的分位点，其值由P_FA和自由度(m－4)确定，应有以下等式成立：

$ P_{\mathrm{FA}}=P\left(\frac{\hat{\boldsymbol{v}}^{\mathrm{T}} {\boldsymbol{P}} \hat{\boldsymbol{v}}}{\sigma_0^2}>\alpha^2\right)=\int_{a^2}^{\infty} f_{\chi_{m-4}^2}(x) \mathrm{d} x $

(6)

存在故障时，概率密度用f_{χm－4, λ²}(x)表示。当发生伪距故障时，如果系统未发出监测报警，则为漏检(missed detection，MD)，漏检率(P_MD)可表示为P_MD=P($ \frac{{{{\mathit{\boldsymbol{\hat v}}}^{\rm{T}}}\mathit{\boldsymbol{P\hat v}}}}{{\sigma _0^2}}$＜α²|故障)。存在故障时，λ由α²、P_MD和自由度(m－4)确定，应有以下等式成立：

$ P_{\mathrm{MD}}=P\left(\frac{\hat{\boldsymbol{v}}^{\mathrm{T}} {\boldsymbol{P}} \hat{\boldsymbol{v}}}{\sigma_0^2}<\alpha^2\right)=\int_0^{a^2} f_{\chi_{m-4, \lambda}^2}(x) \mathrm{d} x $

(7)

一般取检验统计量(T)为：

$ T=\sqrt{\hat{\boldsymbol{v}}^{\mathrm{T}} {\boldsymbol{P}} \hat{\boldsymbol{v}} /(m-4)} $

(8)

在导航过程中，将实时解算的T值与监测门限(T_D)进行比较，T_D计算可参考文献[19]。若T＜T_D，则系统无故障，反之，系统存在故障。

2 RAIM可用性分析方法

在进行RAIM之前，首先根据性能指标对当前可见星的几何分布进行判断，分析是否适合进行完好性监测，即判断RAIM算法是否可用。

2.1 监测统计量与定位误差的关系

RAIM可用性判断时需要构建T与定位误差之间的数学关系。由于RAIM不具备垂直导航能力，因此本文只讨论水平保护级别(horizontal protection level，HPL)。

图 1中横轴表示T，纵轴表示水平径向误差(horizontal radial error，HRE)，将这2个变量取比值后得到图中斜线斜率(SLOPE)，即

$ \mathrm{SLOPE}=\frac{\mathrm{HRE}}{T} $

(9)

从图 1可以看出，T与HRE呈线性关系，沿此倾斜线，当伪距观测量中存在偏差时，T与HRE线性增加。偏差较小时，HRE小于HPL且T小于T_D，系统处于正常状态(图 1左下角)，P_MD较低；偏差较大时，HRE大于HPL且T大于T_D，系统能正确监测偏差(图 1右上角)，P_MD也较低；当偏差处于两者之间时，P_MD可能达到最大(图 1左上角)。对应相同T值时，SLOPE越大，卫星发生故障时MD概率越高。对应卫星的SLOPE_s为：

$ \mathrm{SLOPE}_{\mathrm{s}}=\frac{\mathrm{HPL}_{\mathrm{s}}}{T} $

(10)

在m个观测量中，HPL可表示为：

$ \mathrm{HPL}=\max \left(\mathrm{SLOPE}_{\mathrm{s}}\right) \cdot T=\mathrm{SLOPE}_{\max } \cdot T $

(11)

RAIM可用性分析就是利用SLOPE_s在各历元的最大值SLOPE_max进行判定，依据SLOPE_max求出HPL，并与水平告警限值(horizontal alarm limit，HAL)进行比较，如果HPL＜HAL，说明RAIM算法可用，否则不可用。

2.2 双星故障时MM求解HPL

在忽略观测伪距随机误差的情况下，当观测误差ε在第i个和第j个观测值中存在粗差b_i和b_j时，则有：

$ \boldsymbol{\varepsilon}=\left[\begin{array}{lllllll} 0 & \cdots & b_i & \cdots & b_j & \cdots & 0 \end{array}\right]^{\mathrm{T}} $

(12)

根据式(5)可以得到定位解的误差向量为：

$ \begin{array}{c} \Delta \hat{\boldsymbol{x}}_{i j}=\left[\begin{array}{ccccccc} A_{11} & \cdots & A_{1 i} & \cdots & A_{1 j} & \cdots & A_{1 m} \\ A_{21} & \cdots & A_{2 i} & \cdots & A_{2 j} & \cdots & A_{2 m} \\ A_{31} & \cdots & A_{3 i} & \cdots & A_{3 j} & \cdots & A_{3 m} \\ A_{41} & \cdots & A_{4 i} & \cdots & A_{4 j} & \cdots & A_{4 m} \end{array}\right]\left[\begin{array}{c} 0 \\ \vdots \\ b_i \\ \vdots \\ b_j \\ \vdots \\ 0 \end{array}\right]\\ =\left[\begin{array}{l} A_{1 i} b_i+A_{1 j} b_j \\ A_{2 i} b_i+A_{2 j} b_j \\ A_{3 i} b_i+A_{3 j} b_j \\ A_{4 i} b_i+A_{4 j} b_j \end{array}\right] \end{array} $

(13)

则有：

$ \mathrm{HRE}_{i j}^2=\left(A_{1 i} b_i+A_{1 j} b_j\right)^2+\left(A_{2 i} b_i+A_{2 j} b_j\right)^2 $

(14)

将式(12)代入式(5)可得：

$ \begin{array}{c} \left(\hat{\boldsymbol{v}}^{\mathrm{T}} {\boldsymbol{P}}\hat{\boldsymbol{v}}\right)_{i j}=\boldsymbol{\varepsilon}^{\mathrm{T}} \boldsymbol{P}(\boldsymbol{I}-\boldsymbol{B}) \boldsymbol{\varepsilon}=\left[\begin{array}{ll} b_i P_i & b_j P_j \end{array}\right] \\ {\left[\begin{array}{cc} 1-B_{i i} & -B_{i j} \\ -B_{i j} & 1-B_{j j} \end{array}\right]\left[\begin{array}{l} b_i \\ b_j \end{array}\right]=b_i^2 P_i\left(1-B_{i i}\right)+} \\ b_j^2 P_j\left(1-B_{j j}\right)-b_i b_j P_i B_{i j}-b_i b_j P_j B_{i j} \end{array} $

(15)

令$ k=\frac{b_i}{b_j} $, 则：

$ \begin{gathered} \operatorname{SLOPE}_{i j}^2=\frac{\operatorname{HRE}_{i j}^2 \cdot(m-4)}{\left(\hat{\boldsymbol{v}}^{\mathrm{T}} \boldsymbol{P} \hat{\boldsymbol{v}}\right)_{i j}}= \\ \frac{\left[\left(k A_{1 i}+A_{1 j}\right)^2+\left(k A_{2 i}+A_{2 j}\right)^2\right] \cdot(m-4)}{k^2 P_i\left(1-B_{i i}\right)+P_j\left(1-B_{j j}\right)-k P_i B_{i j}-k P_j B_{i j}} \end{gathered} $

(16)

以往对k值主要采用矩阵最大特征值方法(MMEM)^{[9, 18, 20]}进行求解，但其也可采用高等数学中求极值的方法获取，具体如下：

为求式(16)的极大值，两边取对数，并对式中k求偏导，可得一元二次方程式：

$ a k^2+b k+c=0 $

(17)

设k₁和k₂分别为式(17)两根，将k₁、k₂分别代入式(16)，则当前历元的(SLOPE_ij)_max可表示为：

$ \begin{gathered} \left(\mathrm{SLOPE}_{i j}\right)_{\max }= \\ \max \left(\left(\mathrm{SLOPE}_{i j}\right)_{k_1}, \left(\mathrm{SLOPE}_{i j}\right)_{k_2}\right) \end{gathered} $

(18)

即可得到双星故障时的HPL为^[19]：

$ \mathrm{HPL}=\left(\mathrm{SLOPE}_{i j}\right)_{\max } \cdot \sigma_0 \cdot \sqrt{\lambda} / \sqrt{(m-4)} $

(19)

2.3 双星故障时MMEM求解HPL

目前求解HPL主要采用MMEM方法，该方法主要是构建关于粗差向量b的二次型矩阵，根据矩阵的最大值为矩阵最大特征值这一定律求解HPL，具体计算过程可参考文献[18]。

3 数据处理 3.1 实验数据和处理策略

实验数据采用国际GNSS监测评估系统(international GNSS monitoring and assessment system，IGMAS)提供的中国境内BJF1、CHU1、GUA1、KUN1、LHA1、SHA1、WUH1、XIA1共8个连续跟踪站2020-09-06的观测数据和星历文件。通过自编程序统计这些测站在截止高度角5°以上的BDS卫星可见性、DOP值和基于2种方法的双星故障HPL值，比较并分析中国境内BDS可用性和2种双星故障RAIM可用性的差异。具体数据情况及处理策略为：1)BJF1和CHU1站观测时段为00:00:00~09:59:30；GUA1、KUN1、LHA1、SHA1站观测时段为00:00:00~15:59:30；WUH1和XIA1站观测时段为00:00:00~16:59:30；2)考虑到每个测站不完全包含双频或三频观测值，因此采用BDS B1I频率的单频观测值；3)全球星历文件中每隔1 h给出用于计算每个历元的BDS卫星轨道参数、卫星钟差参数和卫星状态信息，其1 d的全球星历文件可从MGEX网站下载；4)电离层延迟采用Klobuchar模型进行改正，对流层延迟采用Saastamoinen模型进行改正；5)式(2)中P采用已知的高度角定权法^[21]，定权时考虑3种异构卫星之间的差异，将GEO和IGSO+MEO卫星的权比关系设置为1 ∶5^[22]。

3.2 卫星可见数和DOP值

利用8个IGMAS站的BDS观测数据和全球星历文件，计算出每个测站每个历元对应卫星的高度角，并统计可见卫星数、平均卫星数、GDOP值及其平均值。结果表明：1)中国境内8个IGMAS站的BDS可见卫星数最小值为17颗，完全能够满足RAIM可用性评估；2)GDOP最小值、平均值和最大值的均值分别为1.117、1.437、2.186，说明BDS卫星在中国境内空间分布情况很好。

3.3 2种方法计算双星故障时HPL比较

目前BDS在中国区域的定位精度为5 m左右^[4]。保守起见，将σ₀设置为8 m，P_MD设置为10^－3，P_FA设置为10^－5[19]。利用2种方法分别计算BDS双星故障时HPL值，图 2~5分别为8个IGMAS站的HPL坐标时间序列，表 1(单位m) 为8个IGMAS站的HPL最小值、最大值和平均值。通过分析可知：

1) 从整体上看，MMEM计算的HPL小于MM计算的HPL，表明MMEM得到的RAIM可用性高。利用MM求解的8个IGMAS站的HPL最小值、平均值、最大值的均值分别为124.921 m、233.141m、3 020.005 m；MMEM求解的HPL最小值、平均值和最大值的均值分别为219.093 m、515.110 m、1 302.518 m。

2) 2种方法得到的HPL时间序列中，大部分历元的HPL值一致性很好，HPL平均值的差值约为281.969 m，主要原因可能为这2种方法的计算公式中所用参数个数和精度存在差异。

3) 在中国区域，双星故障条件下HPL值的空间分布存在差异，表明RAIM可用性与地理位置相关。

3.4 计算时间比较

表 2(单位s)为2种方法的计算时间，从表中可以看出，MM和MMEM的平均计算时间分别为11.733 s和21.387 s，MM比MMEM计算时间短。

3.5 2种方法双星故障时可用性统计

飞机在飞行阶段包含非精密进近(non-precision approach，NPA)、终端、本土航路和远洋航路4个飞行阶段^[23]，4个阶段的HAL值分别为555.6 m、1 852 m、5 556 m和7 408 m。基于2种方法得到双星故障条件下HPL值，根据每个测站计算出观测历元中各阶段的可用性百分比，再计算8个测站的可用性比值均值，得到每个阶段的可用性结果(表 3)。从表中可以看出：1)MM计算的可用性在航路和远洋阶段达到100%，而MMEM计算的RAIM可用性除在NPA阶段外，其他阶段的可用性均达到100%；2)结合所计算的HPL值，虽然MMEM在NPA阶段的可用性不及MM，但整体上MMEM可用性优于MM。

4 结语

本文从双星存在故障的前提条件出发，在总结已有的最大特征值计算RAIM可用性方法的基础上，提出RAIM可用性评估的极大值方法，并以2020-09-06中国区域8个IGMAS连续跟踪站的BDS观测数据为例，比较和分析中国境内HPL变化及RAIM可用性情况，得到以下结论：

1) 中国境内可见卫星数和GDOP完全能够满足可用性计算需求，中国境内BDS目前至少可观测到17颗可见卫星，GDOP值最大为2.186，可为RAIM可用性计算和故障监测提供丰富的数据源。

2) 虽然MM计算的RAIM可用性在NPA阶段高于MMEM，且耗时较少，但总体来说，MMEM计算HPL的时间序列小于MM，即对应的RAIM可用性更高。

3) 2种方法计算的HPL存在差异，主要原因可能为各方法所用参数个数及本身精度存在差异，且RAIM可用性计算中涉及到完好性风险参数，不同的参数初值会对RAIM可用性评估产生不同影响，因此参数取值应引起重视。