基于伪距残差的接收机自主完好性监测(receiver autonomous integrity monitoring, RAIM)方法是利用接收机内部观测信息的冗余性进行一致性检验来实现故障的检测与识别^[1-2]，是增强GNSS用户端完好性的主要手段之一，观测信息的冗余度是影响RAIM检测性能的重要因素之一。另外，为满足故障检测的漏检概率和虚警概率需求，卫星与接收机之间的几何构型必须具备一定的质量保证^[3]，通常采用计算保护水平与要求的告警限值进行比较的方法来判断基于残差的RAIM的可用性。随着低轨卫星定轨技术的发展，低轨卫星实时自主定轨的精度可达m级^[4-5]，采用精密定轨方式后精度可达cm级^[6-7]，因此将低轨卫星平台作为导航信号源可以增加导航用户观测信息的冗余度，丰富用户与可见卫星间的几何构型。

目前已建成或提出的低轨导航增强星座有铱星NEXT、Kepler系统及国内的鸿雁星座、微厘空间等^[8-10]。本文通过分析信号波束角对低轨星座覆盖性的影响，基于单故障假设的RAIM保护水平计算模型，将保护水平改善值及保护水平改善比作为评估低轨星座对导航系统完好性增强的性能指标，并对低轨星座在不同高度遮蔽角和等效噪声水平下对北斗导航系统保护水平及RAIM可用性等方面的增强性能进行分析。

1 RAIM算法及其可用性分析 1.1 基于残差的单故障探测模型

线性化后的导航接收机伪距观测方程为：

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{Hx}} + \mathit{\boldsymbol{\varepsilon }} + \mathit{\boldsymbol{b}} $

(1)

式中，y为n×1伪距观测向量；H为n×4观测矩阵；x为4×1待估状态向量；ε为n×1伪距观测噪声向量，表示各伪距观测量的测量误差，ε~N(0, D)，D为观测噪声的协方差矩阵；b为n×1偏差向量，表示各观测量中的测量偏差，在无故障模式下为零向量；n为可见卫星数量。

状态向量的加权最小二乘估计解为：

$ \mathit{\boldsymbol{\hat x}} = {({\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{WH}})^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{Wy}} = \mathit{\boldsymbol{Py}} $

(2)

式中，W为权矩阵，W^-1=diag(σ₁², σ₂², …, σ_n²)。令P=(H^TWH)^-1H^TWy，则有$ \mathit{\boldsymbol{\hat x}}$= Py。

根据状态估计向量计算伪距残差向量：

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{v}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{H\hat x}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{H}}{{({\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{WH}})}^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{Wy}} = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ({\mathit{\boldsymbol{I}}_n} - \mathit{\boldsymbol{H}}{{({\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{WH}})}^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{W}})(\mathit{\boldsymbol{\varepsilon }} + \mathit{\boldsymbol{b}})} \end{array} $

(3)

记S=I_n-H(H^TWH)^-1H^TW，则有 v=S(ε+b)。矩阵S为幂等矩阵，则残差平方和可表示为：

$ \begin{array}{*{20}{c}} {{\rm{SSE}} = {\mathit{\boldsymbol{v}}^{\rm{T}}}\mathit{\boldsymbol{v}} = {{(\mathit{\boldsymbol{\varepsilon }} + \mathit{\boldsymbol{b}})}^{\rm{T}}}{\mathit{\boldsymbol{S}}^2}(\mathit{\boldsymbol{\varepsilon }} + \mathit{\boldsymbol{b}}) = }\\ {{{(\mathit{\boldsymbol{\varepsilon }} + \mathit{\boldsymbol{b}})}^{\rm{T}}}\mathit{\boldsymbol{S}}(\mathit{\boldsymbol{\varepsilon }} + \mathit{\boldsymbol{b}})} \end{array} $

(4)

在无故障模式下，伪距残差向量中各分量为相互独立的正态分布随机误差N(0, σ₀²)，依据统计分布理论，归一化统计量SSE/σ₀²服从自由度为(n-4)的卡方分布，即SSE/σ₀²~χ²(n-4)；在故障模式下，SSE/σ₀²服从自由度为(n-4)的非中心卡方分布，即SSE/σ₀²~χ²(n-4, λ)，λ为非中心化参数。

以单位权中误差${\hat \sigma } $作为检验统计量进行故障检测：

$ \hat \sigma = \sqrt { {\rm{SSE}} /(n - 4)} $

(5)

根据卡方分布的自由度(n-4)和给定的最大误检概率P_fa，可利用数值积分的方法计算检验门限值T_D：

$ {P_{{\rm{fa}}}} = \int_{T_D^2}^{ + \infty } {{f_{\chi _{n - 4}^2}}} (x){\rm{d}}x $

(6)

通过比较检验统计量与检测门限值的大小来实现对故障的检测，在检测到故障后采用巴尔达数据探测法对故障星进行识别及隔离。

1.2 RAIM可用性判定

RAIM可用性判定包括冗余度和几何构型2个方面，在任一时间段[t₁, t₂]内的RAIM可用性可表示为：

$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{av}}{{\rm{a}}_{{\rm{RAIM}}}} = \\ \frac{{\sum\limits_{t = {t_1},{t_{{\rm{step }}}} = T}^{{t_2}} {{\rm{ bool }}} \left( {{N_{{\rm{sat }}}} \ge {N_{{\rm{req }}}}} \right) \cdot {\rm{ bool }}\left( {{\rm{P}}{{\rm{L}}_t} \le {\rm{AL}}} \right)}}{{\sum\limits_{t = {t_1},{t_{{\rm{step }}}} = T}^{{t_2}} 1 }} \end{array} $

(7)

式中，N_sat和N_req分别为当前时刻的可见卫星数和RAIM需求的卫星数。保护水平(protection level，PL)为当前时刻RAIM的最大不可检测定位误差，通过比较保护水平与给定的告警限值(alert limit，AL)来判断当前时刻的几何构型是否满足完好性指标需求。保护水平包括水平保护水平(HPL)和垂直保护水平(VPL)，常用的表达方式为最大特征斜率与卡方检验最小可检测偏差的乘积：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{HPL}} = {\rm{Hslop}}{{\rm{e}}_{{\rm{max}}}} \times {\rm{pbias}}}\\ {{\rm{VPL}} = {\rm{Vslop}}{{\rm{e}}_{{\rm{max}}}} \times {\rm{pbias}}} \end{array}} \right. $

(8)

式中，pbias为卡方检验的最小可检测偏差，pbias= $ \sqrt \lambda $，λ为卡方检验的非中心化参数，可根据给定的漏检概率P_MD由下式计算得到：

$ 1 - {P_{{\rm{MD}}}} = \int_0^{{T_D}} {{f_{\chi _{n - 4}^2}}} (x,\lambda ){\rm{d}}x $

(9)

Hslope_max和Vslope_max分别为所有可见卫星中最大的水平特征斜率和最大的垂直特征斜率：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{Hslop}}{{\rm{e}}_{{\rm{max}}}} = \mathop {{\rm{max}}}\limits_i \{ {\rm{Hslop}}{{\rm{e}}_i}\} }\\ {{\rm{Vslop}}{{\rm{e}}_{{\rm{max}}}} = \mathop {{\rm{max}}}\limits_i \{ {\rm{Vslop}}{{\rm{e}}_i}\} } \end{array},i = 1,2, \cdots ,n} \right. $

(10)

特征斜率是描述卫星在卫星集合中相对位置关系的参数，卫星的特征斜率越大，其在出现故障时越难以被正确检测出来。特征斜率由卫星在矩阵P和S中对应的分量计算得到，第i颗卫星的特征斜率为：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{Hslope}}{{\rm{ }}_i} = {\sigma _i}\sqrt {(\mathit{\boldsymbol{P}}_{1i}^2 + \mathit{\boldsymbol{P}}_{2i}^2)} /\sqrt {{\mathit{\boldsymbol{S}}_{ii}}} }\\ {{\rm{Vslope}}{{\rm{ }}_i} = {\sigma _i}\sqrt {\mathit{\boldsymbol{P}}_{3i}^2} /\sqrt {{\mathit{\boldsymbol{S}}_{ii}}} } \end{array}} \right. $

(11)

式中，σ_i为第i颗卫星的距离测量误差。

当水平和垂直2个方向同时满足几何条件要求时，认为RAIM在当前几何构型下可用，即

$ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{bool}}({\rm{P}}{{\rm{L}}_t} \le {\rm{AL}}) = }\\ {{\rm{bool(HP}}{{\rm{L}}_t} \le {\rm{HAL}},{\rm{VP}}{{\rm{L}}_t} \le {\rm{VAL}})} \end{array} $

(12)

2 基于低轨星座的完好性监测增强

低轨星座已被国家综合导航定位授时(PNT)体系纳入建设体系^[11]，目前国内外大部分在建的低轨星座也都将导航功能作为星座的重要功能之一。

2.1 低轨导航信号波束角与高度遮蔽角

卫星信号的对地覆盖性取决于卫星的轨道高度和信号波束角，由于低轨卫星的轨道高度较低，其信号波束角会影响用户观测的高度遮蔽角。地球半径R=6 371 km，设卫星的轨道高度为h，则卫星与地球切线和卫星与地心连线之间的夹角为：

$ \alpha = {\rm{arcsin}}\left( {\frac{R}{{R + h}}} \right) $

(13)

对于GNSS导航卫星，MEO卫星的α一般在12°~15°之间，GEO卫星和IGSO卫星的α更小，而GNSS导航信号的半波束角一般为16°~23°。因此，在设置导航卫星的高度遮蔽角时通常仅需考虑信号遮挡等环境因素^[12]。

低轨卫星的轨道高度通常为300~1 500 km，根据式(13)计算得到对应的α约为54°~72°。因此，当低轨卫星信号半波束角小于α时，还需要考虑低轨导航信号波束角对用户观测的高度遮蔽角的限制。

2.2 低轨增强的RAIM可用性

考虑目前低轨卫星可实现精度相对较高的定轨和时间同步，在本文分析中认为加入低轨星座不会引入新的未知状态参数，则低轨增强的线性化观测矩阵为：

$ {\mathit{\boldsymbol{H}}_{{\rm{BDS + LEO}}}} = \left[ {\begin{array}{*{20}{c}} {{h_{1,1}}}&{{h_{1,2}}}&{{h_{1,3}}}&1\\ \vdots & \vdots & \vdots & \vdots \\ {{h_{n,1}}}&{{h_{n,2}}}&{{h_{n,3}}}&1\\ {{h_{(n + 1),1}}}&{{h_{(n + 1),2}}}&{{h_{(n + 1),3}}}&1\\ \vdots & \vdots & \vdots & \vdots \\ {{h_{(n + m),1}}}&{{h_{(n + m),2}}}&{{h_{(n + m),3}}}&1 \end{array}} \right] $

(14)

式中，m为可见的低轨卫星数量。

在加入低轨星座后，保护水平为：

$ \left\{ \begin{array}{l} {\rm{HP}}{{\rm{L}}^{{\rm{BDS}} + {\rm{LEO}}}} = {\rm{HS}}{\mathop{\rm lope}\nolimits} _{\max }^{{\rm{BDS}} + {\rm{LEO}}} \times {\rm{pbia}}{{\rm{s}}^{{\rm{BDS}} + {\rm{LEO}}}}\\ {\rm{VP}}{{\rm{L}}^{{\rm{BDS + LEO}}}} = {\rm{VSlope }}_{\max }^{{\rm{BDS}} + {\rm{LEO}}} \times {\rm{pbia}}{{\rm{s}}^{{\rm{BDS}} + {\rm{LEO}}}} \end{array} \right. $

(15)

由于特征斜率表征的是卫星之间的相对位置关系，因此与DOP值在加入低轨卫星后单调递减的变化趋势存在差异，保护水平在加入低轨卫星后可能变大或变小。定义保护水平改善值δPL为：

$ \left\{ {\begin{array}{*{20}{l}} {\delta {\rm{HPL}} = {\rm{HP}}{{\rm{L}}^{{\rm{BDS}}}} - {\rm{HP}}{{\rm{L}}^{{\rm{BDS + LEO}}}}}\\ {\delta {\rm{VPL}} = {\rm{VP}}{{\rm{L}}^{{\rm{BDS}}}} - {\rm{VP}}{{\rm{L}}^{{\rm{BDS + LEO}}}}} \end{array}} \right. $

(16)

当δPL < 0时，说明加入低轨卫星会使保护水平增大，即可能会增大当前导航系统用户完好性的风险。

定义任一时间段[t₁, t₂]内的保护水平改善比为：

$ {\rm{ rati}}{{\rm{o}}_{{\rm{PL,improvement}}}}{\rm{ }} = \frac{{\sum\limits_{t = {t_1},{t_{{\rm{step }}}} = T}^{{t_2}} {{\rm{ bool }}} \left( {\delta {\rm{P}}{{\rm{L}}_t} > 0} \right)}}{{\sum\limits_{t = {t_1},{t_{{\rm{step }}}} = T}^{{t_2}} 1 }} $

(17)

由式(11)可知，影响δPL的因素主要包括2个方面：1)低轨卫星与北斗导航卫星之间的相对几何构型；2)低轨卫星与北斗导航卫星的等效观测噪声水平β=σ_LEO/σ_BDS，观测噪声包含测距误差、轨道误差和星钟误差等。本文将通过仿真实验来分析高度遮蔽角与等效噪声水平对保护水平和RAIM可用性改善的影响。

3 仿真结果与分析

仿真实验所用数据为北斗三号卫星的轨道数据，包括29颗已发射卫星2020-03-26 00:00:00~2020-03-27 00:00:00的实际轨道数据和1颗未发射卫星的设计轨道数据，其中高度遮蔽角设为5°，低轨星座的轨道参数参照微厘空间的星座构型，包括120颗倾斜轨道卫星(Walker 120/12/0，轨道倾角55°，轨道高度980 km)和30颗极地轨道卫星(Walker 30/3/0，轨道倾角85°，轨道高度1 250 km)，仿真时长为24 h，采样间隔为60″，给定的虚警概率为3.3×10^-7，漏检概率为1×10^-3，先验观测噪声均方根为6 m。

将全球区域划分为6°×6°网格点，单北斗系统的保护水平平均值和RAIM可用性见图 1。由图可知，IGSO卫星和GEO卫星会增大北斗系统对东半球区域的覆盖性，东半球的保护水平和RAIM可用性明显优于西半球，高纬度和低纬度地区的保护水平和RAIM可用性优于中纬度地区。

分别将低轨卫星的高度遮蔽角设为5°、10°、15°和20°，加入低轨星座后HPL和VPL改善值的平均值和RAIM可用性见图 2和3。

由图 2和3可以看出，加入低轨星座后保护水平相对较差的中纬度地区得到明显改善，低轨卫星的高度遮蔽角越小，对地覆盖性越好，对保护水平和RAIM可用性的改善效果越好。表 1为不同高度遮蔽角下HPL和VPL改善值。

将低轨卫星的高度遮蔽角设为5°，等效噪声水平分别设为2/3和4/3，对应的HPL和VPL改善值和RAIM可用性见图 4。由图可知，等效噪声水平越低，低轨星座的增强效果越好。

表 2为不同等效误差水平下HPL和VPL改善值的数值特征，相比于单北斗系统，HPL改善值的平均值分别降低37.21%、35.54%和30.08%，VPL改善值的平均值分别降低33.14%、31.66%和26.92%。

图 5为不同等效误差水平下低轨星座的保护水平改善比，由图可知，加入低轨星座后东半球中低纬度地区保护水平改善比相对较低，且等效误差水平越高，保护水平改善比越低。

4 结语

保护水平和RAIM可用性取决于卫星与接收机之间的几何构型及导航信号的等效误差水平，在导航星座的基础上加入适当的低轨星座可以增加导航用户可见卫星数，改善接收机与卫星间的几何构型。根据仿真结果可以得出以下结论：

1) 与导航卫星不同，低轨卫星的高度遮蔽角同时受环境遮挡和信号波束角限制，在分析低轨星座对地覆盖性时，应根据天线的波束角设置低轨卫星的高度遮蔽角；

2) 低轨卫星的高度遮蔽角越大，对地覆盖性越差，对保护水平和RAIM可用性的增强效果越差，在低轨导航增强星座的设计中应予以考虑；

3) 低轨导航信号的等效误差水平越高，对保护水平和RAIM可用性产生负面影响的比例越大，尤其是对保护水平原本相对较好的地区。因此，保证低轨卫星具备较高的轨道和星钟精度对提升低轨增强性能至关重要。