CN105044738A - 一种接收机自主完好性监视的预测方法及预测系统 - Google Patents

Abstract

本发明涉及一种接收机自主完好性监视(RAIM)的预测方法及预测系统，包括以下步骤：采集三维集地形数据、飞行计划数据、GPS卫星历书数据；根据上述计算出飞行计划数据中某个时刻的GPS卫星位置坐标、飞机位置坐标、航段地形遮蔽角和卫星相对飞机几何分布；根据地形遮蔽角、卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星；分析卫星伪距离测量误差，并构建权系数矩阵；根据加权水平保护值(HPL)算法，计算出水平保护值；将水平保护值和水平告警门限值比较，获取接收机自主完好性监视结果。在筛选可见卫星的步骤中，增加地形遮蔽角的筛选步骤，提高了接收机自主完好性监视(RAIM)预测的准确性和可信度。

Description

一种接收机自主完好性监视的预测方法及预测系统

技术领域

本发明涉及民用航空飞行运行签派放行技术支持领域，特别是一种接收机自主完好性监视的预测方法及系统。

背景技术

接收机自主完好性监视(RAIM)是一种仅利用GNSS伪距测量信息进行完好性监视的技术。RAIM算法内嵌在接收机内部，接收机在接收GNSS测量信号进行定位解算的同时执行RAIM算法。RAIM算法的完好性监视基本原理为检测各个卫星测量值之间的一致性，因此需要有冗余观测信息。通常只有当接收机能获取五个以上卫星测量值时，RAIM算法才能有效执行，实现故障检测(FD)，即判断接收到的卫星测量值是否存在异常。当获取六个以上卫星测量值时，RAIM算法能够进行故障识别(FDE)，判断哪个测量值存在异常。RAIM不需要借助外界其他设备和信息，检测速度快，方便易行，对于高动态运行的民航飞行运行而言，能够满足其实时监视的需求。另外，RAIM技术应用成本低，易实现，在全球范围内均可使用，成为目前使用最为广泛的一种完好性监视手段。

GNSS卫星信号属于视距信号，对于非透明障碍物穿透能力很弱。飞机在运行中，特别是在地形复杂的机场终端区域运行时，飞行高度低，部分GNSS信号会受到附近山体或者障碍物的遮挡，导致信号丢失。传统RAIM预测系统大都未考虑地形遮蔽的影响，通常以机载接收机规定的仰角门限预测观测到的可见卫星情况。我国地貌起伏不平，截止2014年底拥有着全世界42个高高原机场(海拔超过2438米)中的15个，并且今后还会有更多的西部山区高高原机场建设开通。这些特殊机场地形复杂，飞机通常沿着河谷实施离场和进近，地形遮蔽情况严重，有些机场的地形遮蔽角甚至达到20度以上，严重影响机载GNSS接收机有效地接收GNSS卫星信号。如果RAIM预测算法不对地形遮蔽加以考虑，仅仅以接收机门限仰角作为可见卫星的仰角门限，其预测结果不能准确直观的反映实际飞行中的RAIM可用性情况，存在安全隐患。

发明内容

本发明提供一种接收机自主完好性监视的预测方法及系统，以解决上述预测结果不能准确直观的反映实际飞行中的RAIM可用性情况，存在安全隐患的技术问题。

为了解决上述技术问题，本发明提供一种接收机自主完好性监视的预测方法，包括以下步骤：

步骤S1，采集三维集地形数据、飞行计划数据、GPS卫星历书数据；

步骤S2，根据所述三维地形数据、所述飞行计划数据和所述GPS卫星历书数据计算出所述飞行计划数据中某个时刻的GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角；

根据所述GPS卫星位置坐标、所述飞机位置坐标和所述航段地形遮蔽角计算卫星相对飞机几何分布；

步骤S3，根据所述地形遮蔽角、所述卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星；

步骤S4，根据卫星相对飞机几何分布分析卫星伪距离测量误差，并根据所述伪距离测量误差构建权系数矩阵；

步骤S5，将所述观测矩阵和权系数矩阵代入加权HPL算法，计算出水平保护值；

步骤S6，将所述水平保护值和水平告警门限值比较，获取所述接收机自主完好性监视结果。

本发明的有益效果是：通过在筛选可见卫星的步骤中，增加地形遮蔽角的筛选步骤，提高了接收机自主完好性监视预测的准确性，特别是在地形复杂的机场终端区域运行时，飞行高度低，部分GNSS信号会受到附近山体或者障碍物的遮挡，这样增加地形遮蔽角的筛选步骤，去除山体或者障碍物的遮挡的影响，使得接收机自主完好性监视预测的准确性和可信度更高。

进一步，所述步骤S2中计算飞机位置坐标步骤包括，构建基于性能的导航航线模型，所述基于性能的导航航线模型是基于WGS-84基准椭球体的大椭圆航线；

根据所述大椭圆航线模型以及距离计算算法计算某个时刻的飞机的位置坐标(x、y、z)；

其中，所述大椭圆航线的模型为：

$\left\{\begin{array}{c}px+qy+\mathrm{\ξ}z=0\\ \frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{b^2}=1\end{array}\right.$

其中，坐标原点为大椭圆航线模型中心，XY平面为赤道平面，其中X轴指向0°经度，Y轴指向东经90°，Z轴指向地理北极极点；A和B为所述大椭圆航线模型中的两点，A点的空间直角坐标记为(x_A，y_A，z_A)，B点的空间直角坐标记为(x_B，y_B，z_B)；a为基准椭球的长半轴，b为基准椭球的短半轴；p、q和ζ满足如下公式：

$p=\left|\begin{array}{cc}{y}_A& z_A\\ y_B& z_B\end{array}\right|,q=\left|\begin{array}{cc}{x}_B& z_B\\ x_A& z_A\end{array}\right|,\mathrm{\ξ}=\left|\begin{array}{cc}{x}_A& y_A\\ x_B& y_B\end{array}\right|.$

采用上述进一步方案的有益效果是：通过建立基于WGS-84基准椭球体的大椭圆航线，基于性能的导航(PBN)摒弃了传统的台到台飞行模式，采用点到点飞行，即从一个航路点飞往另一个航路点，航路点严格基于WGS-84坐标系定义。而是基于WGS-84基准椭球体的大椭圆航线，这样飞机的位置预测更为准确。

进一步，所述步骤S2中地形遮蔽角的计算包括：在所述三维集地形数据设置高程提取点的经纬度坐标，根据所述经纬度坐标提取相应的高程信息；根据所述高程信息和所述飞机位置坐标计算出某个时刻的地形遮蔽角。

采用上述进一步方案的有益效果是：通过提取三维集地形数据中高程点的经纬度坐标，计算出飞机位置坐标计算出某个时刻的地形遮蔽角，提高了地形遮蔽角的准确性。

进一步，所述步骤S2中计算卫星相对飞机几何分布具体包括：以观测步长对某个时刻上所述飞机位置为中心，建立站心坐标系，计算某个时刻所述卫星相对于所述飞机的卫星观测仰角、卫星观测方位角和卫星观测距离。

采用上述进一步方案的有益效果是：通过建立站心坐标系，能够准确计算出某个时刻所述卫星相对于所述飞机的卫星观测仰角、卫星观测方位角和卫星观测距离。

进一步，所述步骤S3中筛选可见卫星包括：根据所述卫星观测仰角、所述地形遮蔽角和接收机门限仰角，筛选除所述卫星观测仰角小于所述接收机门限仰角的卫星和所述卫星观测仰角小于所述地形遮蔽角的卫星。

采用上述进一步方案的有益效果是：将卫星观测仰角小于地形遮蔽角的卫星筛除，确保预测到的卫星在任意状态下都是可见的，提高了系统预测的准确性。

进一步，根据所述卫星相对飞机几何分布构建观测矩阵G，所述观测矩阵G为： $G=\left[\begin{array}{cccc}cos\left({\mathrm{\θ}}^{\left(1\right)}\right)cos\left({\mathrm{\α}}^{\left(1\right)}\right)& cos\left({\mathrm{\θ}}^{\left(1\right)}\right)\sin \left({\mathrm{\α}}^{\left(1\right)}\right)& sin\left({\mathrm{\θ}}^{\left(1\right)}\right)& 1\\ cos\left({\mathrm{\θ}}^{\left(2\right)}\right)\cos \left({\mathrm{\α}}^{\left(2\right)}\right)& cos\left({\mathrm{\θ}}^{\left(2\right)}\right)\sin \left({\mathrm{\α}}^{\left(2\right)}\right)& \sin \left({\mathrm{\θ}}^{\left(2\right)}\right)& 1\\ ...& ...& ...& ...\\ cos\left({\mathrm{\θ}}^{\left(N\right)}\right)\cos \left({\mathrm{\α}}^{\left(N\right)}\right)& cos\left({\mathrm{\θ}}^{\left(N\right)}\right)\sin \left({\mathrm{\α}}^{\left(N\right)}\right)& sin\left({\mathrm{\θ}}^{\left(N\right)}\right)& 1\end{array}\right]$

其中，α^(N)和θ^(N)分别表示N个可见卫星相对于飞机的卫星观测方位角和卫星观测仰角。

采用上述进一步方案的有益效果是：通过建立观测矩阵G，便于准确计算系统的水平保护值。

进一步，所述步骤S4中卫星伪距离测量误差σi满足如下公式：

${\mathrm{\σ}}_i^2={\mathrm{\σ}}_{i,URA}^2+{\mathrm{\σ}}_{i,iono}^2+{\mathrm{\σ}}_{i,tropo}^2+{\mathrm{\σ}}_{i,mp}^2+{\mathrm{\σ}}_{i,tcvr}^2$

其中，σ_i，URA为卫星时钟和星历误差，σ_i，iono为电离层延时误差，σ_i，tropo为对流层延时误差，σ_i，mp为接收机多路径误差，σ_i，tcvr为接收机热噪声误差，i为1至N；

所述权系数矩阵为：

$W=\left[\begin{array}{cccc}1/{\mathrm{\σ}}_1& & & \\ & 1/{\mathrm{\σ}}_2& & \\ & & ...& \\ & & & 1/{\mathrm{\σ}}_N\end{array}\right].$

采用上述进一步方案的有益效果是：通过建立权系数矩阵，便于准确计算系统的水平保护值。

进一步，所述步骤S5中所述水平保护值的计算公式为：

$HPL=\underset{i}{\max }\sqrt{(A_{w,1i}^2+A_{w,2i}^2)/S_{w,ii}{C}_{ii}}\×\sqrt{\mathrm{\λ}}$

其中，A_W，1i和A_W，2i分别代表矩阵A_W第1行第i列和第2行第i列数值；所述矩阵A_W为：

A_w＝(G^TCG)^-1G^TC

S_W，1i是矩阵S_W的第i行第i列数值，矩阵S_W为：

S_w＝I_N-G(G^TCG)^-1G^TC

C_ii为矩阵C的第i行第i列元素，λ是满足漏检概率要求的非中心分布χ²(N-4，λ)中的非中心参量，N为可见卫星的数量，I_N为N阶单位矩阵，矩阵C表示权系数矩阵W的协方差，矩阵C为：

C＝W^TW。

采用上述进一步方案的有益效果是：通过水平保护值(HPL)的计算公式，准确计算出水平保护值(HPL)，便于将计算出的水平保护值和水平告警门限值比较，准确获取接收机自主完好性监视结果。

本发明还提供一种接收机自主完好性监视的预测系统，其特征在于，包括数据采集模块、第一计算模块、可见卫星筛选模块、权系数矩阵构建模块、水平保护值计算模块和检测结果获取模块；

所述数据采集模块，用于采集三维集地形数据、飞行计划数据、GPS卫星历书数据；

所述第一计算模块，用于根据所述三维地形数据、所述飞行计划数据和所述GPS卫星历书数据计算出所述飞行计划数据中某个时刻的GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角；根据所述GPS卫星位置坐标、所述飞机位置坐标和所述航段地形遮蔽角计算卫星相对飞机几何分布；

所述可见卫星筛选模块，用于根据所述地形遮蔽角、所述卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星；

所述权系数矩阵构建模块，用于根据卫星相对飞机几何分布分析卫星伪距离测量误差，并根据所述伪距离测量误差构建权系数矩阵；

所述水平保护值计算模块，用于将所述观测矩阵和权系数矩阵代入加权HPL算法，计算出水平保护值；

所述检测结果获取模块，用于将所述水平保护值和水平告警门限值比较，获取所述接收机自主完好性监视结果。

本发明接收机自主完好性监视的预测系统的有益效果是：有在可见卫星筛选模块中，根据地形遮蔽角、所述卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星；提高了接收机自主完好性监视预测的准确性，特别是在地形复杂的机场终端区域运行时，飞行高度低，部分GNSS信号会受到附近山体或者障碍物的遮挡，这样增加地形遮蔽角的筛选步骤，去除山体或者障碍物的遮挡的影响，使得接收机自主完好性监视预测的准确性和可行度更高。

进一步，所述预测系统还包括航班计划调整模块，用于根据所述检测结果获取模块获取的接收机自主完好性监视结果，调整航班计划。

采用上述进一步方案的有益效果是：通过航班计划调整模块，充分利用接收机自主完好性监视预测的结果。

附图说明

图1是本发明的预测方法的实施方式控制流程图，

图2是本发明的预测系统的结构框图，

图3是基于WGS-84基准椭球体的大椭圆航线模型图，

图4是内插点坐标计算流程图，

图5是地形遮蔽角计算流程图，

图6是高程提取点选取图，

图7是高程信息提取示意图，

图8是程高信息保存示意图，

图9是各方位地形遮蔽角分布示意图，

图10是飞机位置为原点的站心坐标系示意图，

图11是站心坐标系图，

图12是可见卫星判断流程图，

图13是电离层延时示意图，

图14是GPS各测量误差分布图，

图15是不同卫星数目下的非中心参量分布图。

附图中，各标号所代表的部件列表如下：

1、卫星伪距总测量误差，2、卫星自身误差，3、电离层误差，4、对流层误差，10、数据采集模块，20、第一计算模块，30、可见卫星筛选模块，40、权系数矩阵构建模块，50、水平保护值计算模块，60、检测结果获取模块，70、航班计划调整模块

具体实施方式

下面结合附图和实施方式对本发明作进一步的说明。

本发明的预测方法的实施方式控制流程图参见图1，接收机自主完好性监视的预测方法，包括以下步骤。

步骤S1，采集三维集地形数据、飞行计划数据、GPS卫星历书数据。

步骤S2，根据三维地形数据、飞行计划数据和GPS卫星历书数据计算出飞行计划数据中某个时刻的GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角。

根据GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角计算卫星相对飞机几何分布。

步骤S3，根据地形遮蔽角、卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星。

步骤S4，根据卫星相对飞机几何分布分析卫星伪距离测量误差，并根据伪距离测量误差构建权系数矩阵。

步骤S5，将观测矩阵和权系数矩阵代入加权HPL算法，计算出水平保护值。

步骤S6，将水平保护值和水平告警门限值比较，获取接收机自主完好性监视结果。

通过在筛选可见卫星的步骤中，增加地形遮蔽角的筛选步骤，提高了接收机自主完好性监视预测的准确性，特别是在地形复杂的机场终端区域运行时，飞行高度低，部分GNSS信号会受到附近山体或者障碍物的遮挡，这样增加地形遮蔽角的筛选步骤，去除山体或者障碍物的遮挡的影响，使得接收机自主完好性监视预测的准确性和可信度更高。

本发明的预测系统的结构框图参见图2，得接收机自主完好性监视的预测系统包括数据采集模块10、第一计算模块20、可见卫星筛选模块30、权系数矩阵构建模块40、水平保护值计算模块50、检测结果获取模块60和航班计划调整模块70。

数据采集模块10，用于采集三维集地形数据、飞行计划数据、GPS卫星历书数据。

第一计算模块20，用于根据三维地形数据、飞行计划数据和GPS卫星历书数据计算出飞行计划数据中某个时刻的GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角；根据GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角计算卫星相对飞机几何分布。

可见卫星筛选模块30，用于根据地形遮蔽角、卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星。

权系数矩阵构建模块40，用于根据卫星相对飞机几何分布分析卫星伪距离测量误差，并根据伪距离测量误差构建权系数矩阵。

水平保护值计算模块50，用于将观测矩阵和权系数矩阵代入加权HPL算法，计算出水平保护值。

检测结果获取模块60，用于将水平保护值和水平告警门限值比较，获取接收机自主完好性监视结果。

预测系统还包括航班计划调整模块70，用于根据检测结果获取模块获取的接收机自主完好性监视结果，调整航班计划。

本发明接收机自主完好性监视的预测系统的有益效果是：有在可见卫星筛选模块中，根据地形遮蔽角、卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星；提高了接收机自主完好性监视预测的准确性，特别是在地形复杂的机场终端区域运行时，飞行高度低，部分GNSS信号会受到附近山体或者障碍物的遮挡，这样增加地形遮蔽角的筛选步骤，去除山体或者障碍物的遮挡的影响，使得接收机自主完好性监视预测的准确性和可行度更高。

本发明接收机自主完好性监视的预测方法具体操作如下。

根据所述三维地形数据、飞行计划数据和GPS卫星历书数据计算出所述飞行计划数据中某个时刻的GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角。

其中，计算飞行计划数据中某个时刻的GPS卫星位置坐标的实施过程如下：

美国http：//www.navcen.uscg.gov/网站上每周六的24时更新一次历书参数，即发布新一期的历书。历书参数的有效期为一个星期，其中包含的星期数，用以帮助判断历书的有效性。有效期内的历书参数的计算精度较高，能够满足飞行前的RAIM预测的需要。

下表为下载的历书文件中某颗卫星的一套历书参数。

下面给出了某一颗卫星在未来t时刻的位置坐标计算步骤。

步骤一、计算规划时间t_k

卫星历书中包含的轨道参数以历书参考时间t_oe为基准，为了计算卫星在t时刻的位置，首先必须求出t时刻与参考时间t_oe之间的差异，即

t_k＝t-t_oe

时刻t以秒为单位，从历书发布的那个星期周六晚上12点开始计时。这里必须注意，参考时间t_oe是GPS时间，为了保持时钟基准一致，时刻t也必须换算成GPS时间，UTC与GPST之间的换算关系可参照下表内容。

步骤二、计算卫星的平均角速度

$n_0=\sqrt{\frac{GM}{{a_s}^3}}$

其中GM为地球引力常数，参照下表中的内容。

基本参数	参数值
		椭球体长半轴a	6378.137km
椭球体短半轴b	6356.7523142km
		椭球体极偏率f	1/298.257223561
地球自转角速度Ωe	7.292115146×10-5rad/s
		地球引力常数u＝GM	3.986005×1014m3/s2

步骤三、计算t时刻的平近点角M_k

将参考时间的平近点角M₀代入以下的线性公式：

M_k＝M₀+n₀t_k

M_k应在0到2π之间，可以根据计算结果加减2π，以符合要求。

步骤四、计算t时刻的偏近点角E_k

偏近点角E_k和平近点角M_k关系式如下：

M_k＝E_k-e_ssin(E_k)

其中历书参数e_s是轨道偏心率，这里可以采用迭代的方法，将E_k从以上的开普勒方程中求解出来。

步骤五、计算t时刻的真近点角v_k

真近点角v_k和偏近点角E_k存在以下关系

$v_k=arctan\left(\frac{\sqrt{1-{e_s}^2}sin\left(E_k\right)}{cos\left(E_k\right)-e_s}\right)$

v_k值在-π和π之间。

步骤六、计算t时刻的升交点角距Φ_k

将卫星历书给出的近地点角距w代入下式：

Φ_k＝v_k+w步骤七、计算在卫星运行轨道平面的卫星矢径长度r_k

卫星矢径长度r_k和偏近点角E_k存在以下关系：

r_k＝a_s(1-e_scos(E_k))

步骤八、计算t时刻卫星在轨道平面的位置(x_k’，y_k’)

$\left\{\begin{array}{c}{x}_k^{\′}=r_{k}cos\left({\mathrm{\Φ}}_k\right)\\ y_k^{\′}=r_{k}sin\left({\mathrm{\Φ}}_k\right)\end{array}\right.$

步骤九、计算t时刻的升交点赤经Ω_k

升交点赤经的线性模型如下：

Ω_k＝Ω₀+(Ω-Ω_e)t_k-Ω_et_oe

其中，升交点赤经变化率Ω和升交点赤经Ω₀均由历书提供，Ω_e为地球自转角速度，参考下表中的内容。

基本参数	参数值
		椭球体长半轴a	6378.137km
椭球体短半轴b	6356.7523142km
		椭球体极偏率f	1/298.257223561
地球自转角速度Ωe	7.292115146×10-5rad/
		地球引力常数u＝GM	3.986005×1014m3/s2

步骤十、计算t时刻卫星在WGS-84地心地固直角坐标系中的坐标(x_k，y_k，z_k)，计算公式如下：

$\left\{\begin{array}{c}{x}_k=x_k^{\′}\cos \left({\mathrm{\Ω}}_k\right)y_k^{\′}\cos \left(i_0\right){\mathrm{sin\Ω}}_k\\ y_k=x_k^{\′}\sin \left({\mathrm{\Ω}}_k\right)-y_k^{\′}\cos \left(i_0\right){\mathrm{cos\Ω}}_k\\ z_k=y_k^{\′}\sin \left(i_0\right)\end{array}\right.$

其中

i₀＝δ_i+0.3π

以上便是整个GPS卫星位置计算算法，最终得到各颗卫星在某时刻的位置分布。

其中，飞行计划数据中某个时刻的飞机位置坐标的实施过程如下：

根据飞行计划制定的飞行路径以及飞机在各个飞行阶段的地速和预计到达时刻等信息，推测各个阶段飞机运行的四位维航迹，即三维位置坐标和对应的一维时间矢量，从而实现飞机位置预测。

PBN运行摒弃了传统的台到台飞行模式，采用点到点飞行，即从一个航路点飞往另一个航路点，航路点严格基于WGS-84坐标系定义。因此PBN航线模型不是传统的等角航线或者大圆航线，而是基于WGS-84基准椭球体的大椭圆航线。

大椭圆航线为航路起点和航路终点以及地心构成的空间平面与地球基准椭球面相交而形成的大椭圆上的一段劣弧。

建立如图3所示的直角坐标系，坐标原点为WGS-84椭球模型中心，XY平面为赤道平面，其中X轴指向0°经度，Y轴指向东经90°，Z轴指向地理北极极点。A点和B点分别为PBN航路上的两个航路点，它们之间在基准椭球面上的最短连线所在的曲线模型C，即为大椭圆航线。E点为大椭圆航线和赤道平面的一个交点。

WGS-84基准椭球面表示为：

$\frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{b^2}=1$

其中a为基准椭球的长半轴，b为基准椭球的短半轴，WGS-84坐标系及旋转椭球体基本参数定义下表中内容所列。

基本参数	参数值
		椭球体长半轴a	6378.137km
椭球体短半轴b	6356.7523142km
		椭球体极偏率f	1/298.257223561
地球自转角速度Ωe	7.292115146×10-5rad/s
		地球引力常数u＝GM	3.986005×1014m3/s2

A点的空间直角坐标记为(x_A，y_A，z_A)，B点的空间直角坐标记为(x_B，y_B，z_B)。基于WGS-84坐标系的直角坐标(x，y，z)和大地坐标(φ，λ)的转换公式如下

其中，N为卯酉圆曲率半径，e为椭球偏心率，它与WGS-84基准椭球体的长半轴a和短半轴b存在以下关系：

$e^2=\frac{a^2-b^2}{a^2}$

$N=\frac{a}{\sqrt{1-e^2{\sin }^2\left(\mathrm{\φ}\right)}}$

A点和B点以及地心O构成一个空间平面，该平面与WGS-84基准椭球面的交线就是飞机待飞的大椭圆航线。该平面过坐标系原点(基准椭球体中心)，平面方程可表示如下：

px+qy+ξz＝0

代入A点和B点空间直角坐标，得：

$p=\left|\begin{array}{cc}{y}_A& z_A\\ y_B& z_B\end{array}\right|,q=\left|\begin{array}{cc}{x}_B& z_B\\ x_A& z_A\end{array}\right|,\mathrm{\ξ}=\left|\begin{array}{cc}{x}_A& y_A\\ x_B& y_B\end{array}\right|$

因此，由航路点A、B所构建的PBN航路模型表示如下：

$\left\{\begin{array}{c}px+qy+\mathrm{\ξ}z=0\\ \frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{b^2}=1\end{array}\right.$

由PBN航路模型表达的大椭圆航线方程易知，大椭圆C的长半轴即为赤道半径a，短半轴b’为

$b^{\′}=\sqrt{b^2+\frac{{\mathrm{\ξ}}^2(a^2-b^2)}{p^2+q^2+{\mathrm{\ξ}}^2}}$

大椭圆C和赤道平面的一个交点E的直角坐标为

使用向量夹角公式，可求得矢径和的夹角α₁，矢径和的夹角为α₂，通过椭圆圆弧积分即可得到大椭圆上的劣弧AB的长度S。

$S={\∫}_{{\mathrm{\α}}_1}^{{\mathrm{\α}}_2}\sqrt{a^2+(b^{\′2}-a^2){\sin }^2\mathrm{\θ}}d\mathrm{\θ}$

如图3所示，位置点P处于大椭圆航线C上，距离起航路点A的距离为S₀。假设P点大地坐标为(φ，λ)，由于P点位于大椭圆C上，那么点P、A、B、O共面，即点P满足公式

px+qy+ξz＝0

所表达的平面方程。

根据坐标转换公式

将大地坐标转换成直角坐标，代入平面方程，得到该大椭圆航线上的位置点经纬度关系如下：

联立距离公式：

$S={\∫}_{{\mathrm{\α}}_1}^{{\mathrm{\α}}_2}\sqrt{a^2+(b^{\′2}-a^2){\sin }^2\mathrm{\θ}}d\mathrm{\θ}$

两组非线性方程解算两个未知数，利用matlab软件进行迭代运算，求得内插点坐标，具体算法流程如图4所示。利用matlab软件计算内插点坐标的核心思想是不停地进行迭代运算，直至符合要求达到收敛。其中，门限值G0根据需要的精度要求设定，在机载导航数据库中，位置点的经纬度坐标通常精确到0.01秒，对应的弧度为4.85×10^-8rad，计算得到的距离最小分辨率约为0.3米，因此设定的判决门限G₀不能低于0.3米。另外，在设置内插点P的初始经度时，为了减少迭代运算量，可以采用差值法赋予初始经度值。

本发明的中的地形数据来自Googleearth，通过地形遍历，提取相应的高程信息，GoogleEarth地形数据库数据来源于高分辨率的卫星影像和航拍，其中航拍影像的分辨率最多能够达到1米，而地球探测卫星系统则能够提供最高为0.61米的高精度影像。GoogleEarth通过访问Key-hole的航拍和卫星影像扩展数据库，用户浏览的地球三维地貌有效分辨率至少为100m，在我国大陆能达到30m，针对特定区域(如大城市、著名风景区等)会提供分辨率为0.6m-1m的高精度影像。

所述三维集地形数据设置高程提取点的经纬度坐标，根据所述经纬度坐标提取相应的高程信息；根据所述高程信息和所述飞机位置坐标计算出某个时刻的地形遮蔽角。计算飞机位置坐标计算出某个时刻的地形遮蔽角流程如图5所示。

第一步，设置高程提取点的经纬度坐标；如图6所示，以某个时刻飞机的位置P为中心，在每条方位线上等间隔取高程提取点，步长为d，共取n个点，方位线夹角为α，则需要提取高程数据的方位线共(360/α)条。根据P点位置坐标值，计算各个提取点的经纬度坐标。

第二步，将计算得到的所有提取点经纬度坐标导入excel文档，方便高程提取软件调用，程序如下：

xlswrite(′A.xls′，A，1，′A2′)；

A为存储提取点经纬度坐标值得矩阵，导入到A.xls文档sheet1表，从A列第2行开始导入，如下下表格内容所示。

经度	纬度
		100.1461	29.18865
100.15	29.19089
		100.1539	29.19314
100.1578	29.19538
		100.1616	29.19763
100.1655	29.19988
		100.1694	29.20212
100.1733	29.20437
		100.1772	29.20661
100.1811	29.20886
		100.185	29.2111
100.1889	29.21335

第三步，GoogleEarth高程提取软件调用excel表中存储的提取点经纬度信息，并提取相应的高程信息，如图7所示。

第四步，所有提取点的高程信息提取完毕之后，将高程数据保存到txt文档，如图8中内容所示。

第五步，将txt文档中的高程数据导入到MATLAB软件，计算地形遮蔽角，程序如下：

fid＝fopen(′高程.txt′，′r′)；％打开存储高程数据的txt文档“高程.txt”

H＝fscanf(fid，′％f′)；％将提取到的高程数据存储到矩阵H中

fclose(fid)；％关闭txt文档

某一位置点三维坐标为(E100.1422N29.18644500)，提取点步长D＝500m，方位线夹角为10°，每条方位线取20个点，仿真结果如图13所示。

图9中仿真结果显示，在该位置点70°方位线上的地形遮蔽角最大，为23.9345°，可以认为地形对于该位置点的最大遮蔽角θ_g＝23.9345°。

根据所述GPS卫星位置坐标、所述飞机位置坐标和所述航段地形遮蔽角计算卫星相对飞机几何分布的实施过程如下：

为了更直观的反应GPS卫星星座相对于飞机的分布情况，需要建立一个以飞机所在位置点为坐标原点的站心坐标系。三个坐标轴互相垂直，满足右手定则要求，分别指向正东东方向、正北方向和天顶方向，因此站心坐标系亦称东北天(ENU)坐标系，如图10所示。在站心坐标系中，用卫星相对于飞机的仰角、方位角和距离三要素来表达卫星的位置信息。

某一时刻飞机所在位置点P的大地坐标为(λ，φ，h)，首先需要将其转换成地心地固坐标(x，y，z)，转换方法可参照式(4.2)。对应时刻某颗卫星的位置点S坐标为(x^(s)，y^(s)，z^(s))，则从从用户到该卫星的观测向量为

$\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\\ \mathrm{\Δ}z\end{array}\right]=\left[\begin{array}{c}{x}^{\left(s\right)}\\ y^{\left(s\right)}\\ z^{\left(s\right)}\end{array}\right]-\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$

该卫星在点P处的单位观测量I^(m)为

$I^{\left(m\right)}=\frac{1}{\sqrt{{\mathrm{\Δx}}^2+{\mathrm{\Δy}}^2+{\mathrm{\Δz}}^2}}\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\\ \mathrm{\Δ}z\end{array}\right]$

观测向量[ΔxΔyΔz]^T可等效地表达在以飞机位置点P为中心的站心坐标系中的向量[ΔeΔnΔu]^T，变换关系为

$\left[\begin{array}{c}\mathrm{\Δ}e\\ \mathrm{\Δ}n\\ \mathrm{\Δ}u\end{array}\right]=S\·\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\\ \mathrm{\Δ}z\end{array}\right]$

其中Δe为东向相对距离，Δn为北向相对距离，Δu为天向相对距离。坐标变换矩阵S为

式中，λ和φ分别指的是飞机位置点P的大地坐标经度值和纬度值。坐标变换矩阵S是一个单位正交矩阵，无论是从地心地固直角坐标系转换到站心坐标系，还是反过来，卫星观测矢量的长度保持不变。建立的站心坐标系如图11所示。

得到飞机位置处P点的卫星观测向量[ΔeΔnΔu]^T，接下来可以计算该卫星相对于飞机位置点P的方位与仰角。

卫星仰角θ是观测矢量高出由东向和北向两轴所构成的水平面的角度，即

$\mathrm{\θ}=arcsin\left(\frac{\mathrm{\Δ}u}{\sqrt{{\left(\mathrm{\Δ}e\right)}^2+{\left(\mathrm{\Δ}n\right)}^2+{\left(\mathrm{\Δ}u\right)}^2}}\right)$

卫星仰角θ值在-90°到90°之间，由于地球是个非透明体，因此对于可见GPS卫星来讲，仰角θ的值在大多数情况下大于0°，它的最大值为90°，指卫星位于正头顶方向。

卫星的方位角α定义为北向顺时针转到观测矢量在水平面内的投影方向上的角度，即

$\mathrm{\α}=\arctan \left(\frac{\mathrm{\Δ}e}{\mathrm{\Δ}n}\right)---\left(4.30\right)$

根据所述地形遮蔽角、所述卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星的实施过程如下：

对于接收机而言，卫星是否可见跟该卫星的观测仰角有直接关系，仰角超过观测仰角门限值，即为可见。考虑到过低仰角的GPS卫星伪距测量误差很大，为了保证GPS卫星信号质量，减小定位误差，接收机内部软件中设定了可用仰角门限，低于这个可用仰角的卫星视为被遮蔽，其测量信息将不会参与定位算法。FAATSOC129接收机内置的可用仰角为7.5°，值得注意的是，FAATSOC129接收机将SA政策引起的伪距测量误差进行了编码并写入了定位算法中。FAATSOC146接收机规定的仰角门限为5°，目前大多数接收机采用5°或者更低的遮蔽角。

在民航中，传统的RAIM预测系统对于可见卫星的预测，通常以接收机遮蔽角作为判断标准，比如我国空管部门的RAIM预测系统以5°仰角为可用仰角。在PBN航线上，大部分区域的地形遮蔽角一般不会超过5°，然而在一些地形复杂的山区机场，飞机在实施终端区进近时，由于飞行高度低，地形遮蔽有可能超过接收机仰角门限。事实上，我国地形复杂的特殊机场众多，某些机场所在区域的地形遮蔽角甚至达到20°。因此，飞机在这些区域运行时，仅以接收机内置仰角门限作为可见卫星的判决标准，RAIM预测结果可信度不高，容易造成漏检，存在安全隐患。

可见卫星判断流程如图12所示，其中θ_i是第i颗卫星相对于接收机的观测仰角，θ_R为接收机内置仰角门限，θ_g为飞机所在位置的地形遮蔽角。首先统一根据接收机仰角门限剔除那些卫星仰角低于此门限值的卫星，然后根据飞机位置的地形遮蔽角θ_g对剩下的卫星进行第二次筛选，均满足要求的即为可见卫星。

卫星信号从产生到被接收整个过程所产生的伪距误差主要来自三个部分：与卫星有关误差、与信号传播有关的误差和与接收机有关的误差。其中与卫星有关误差来自两个方面：卫星时钟误差和卫星星历误差，共约产生6米的测量误差；与信号传播有关的误差通常由电离层延时误差、对流层延时误差和多路径误差构成^[30]。卫星伪距测量误差σ_i表达为

${\mathrm{\σ}}_i^2={\mathrm{\σ}}_{i,URA}^2+{\mathrm{\σ}}_{i,iono}^2+{\mathrm{\σ}}_{i,tropo}^2+{\mathrm{\σ}}_{i,mp}^2+{\mathrm{\σ}}_{i,tc\mathrm{\ν}r}^2$

对中的各个误差分量的相应定义下表内容所列。

其中与卫星相关误差σ_i，URA和接收机热噪声误差σ_i，tcvr是定值，对流层延时误差σ_i，tropo和接收机多路径误差σ_i，mp取决于卫星观测仰角θ。电离层延时误差σ_i，iono模型中T_iono是电离层校正值，仅来自卫星星历数据，由于本发明使用历书参数，所以，σ_i，iono定义为

σ_i，iono＝F_pp·τ_vert

其中

${\mathrm{\&tau;}}_{vert}=\left\{\begin{array}{cc}9m& 0\&le;\left|{\mathrm{\&eta;}}_m\right|\&le;20\\ 4.5m& 20<\left|{\mathrm{\&eta;}}_m\right|\&le;55\\ 6m& 55<\left|{\mathrm{\&eta;}}_m\right|\end{array}\right.$

η_m为电离层刺穿点在地面投影点的磁纬度，F_PP为倾斜因子。民用导航用户使用L1频率，单频GPS卫星信号传播过程中，穿过电离层的示意图如图13所示。

电离层刺穿点指的是信号与电离层平均高度面的交点P，它的卫星相对于刺穿点P的天顶角为ζ’，卫星相对于接收机位置点R的天顶角为ζ。天顶角指的是观测矢量和天顶方向的夹角，它与仰角的关系为

$\mathrm{\&zeta;}=\frac{\mathrm{\&pi;}}{2}-\mathrm{\&theta;}8$

倾斜因子F_PP表达式为

F_pp＝1/cosζ′8

根据三角形的正弦定理，可知

$\frac{{\mathrm{sin\&zeta;}}^{\&prime;}}{R_e}=\frac{sin\mathrm{\&zeta;}}{R_e+h_i}8$

其中R_e为地球平均半径，h_i为电离层平均高度，联立表达式

F_pp＝1/cosζ′

可解得倾斜因子F_PP为

$F_{pp}={\&lsqb;1-{\left(\frac{R_{e}sin\mathrm{\&zeta;}}{R_e+h_i}\right)}^2\&rsqb;}^{-1/2}---\left(4.40\right)$

地球平均半径R_e约为6368km，电离层平均高度h_i约为350km，代入上式即可得到仅与卫星观测仰角相关的倾斜因子F_PP表达式。

为了减小计算量，这里给出了一个计算F_PP的近似公式：

$F_{pp}=1+16\&times;{(0.53-\frac{\mathrm{\&theta;}}{\mathrm{\&pi;}})}^3$

式中，θ为卫星观测仰角，在这里以弧度为单位。

另外P点在地面的投影点A磁纬度值η_m，可根据卫星方位角以及A点和R点在地球表面的圆弧角w进行计算。

假如某一时刻，接收机位置点R坐标为N30E100，卫星方位角为80°，不同卫星仰角下的伪距测量误差仿真情况如图14所示。

从图19中可以看出，卫星伪距总测量误差1与卫星仰角息息相关，随着卫星仰角减小，其伪距总测量误差1在逐渐变大，并且当卫星仰角越靠近0°时候，所对应的伪距测量误差急剧增大。另外，电离层误差2、卫星自身误差3和对流层误差4是测量误差的主要来源。卫星自身误差3是一个固定值，因此影响测量误差变化的主要因子是电离层误差2，图中电离层误差曲线分布趋势几乎与总测量误差曲线保持一致。GPS系统采用伪距进行定位，伪距总测量误差1越大引起的最后定位误差也会越大。因此，接收机在设计过程中，为了保证定位精度，会在接收机软件中内置一个卫星可用仰角门限，低于该门限值的卫星信号将不被接收。

GPS接收机在定位解算中，会使用加权最小二乘法解算接收机位置信息。加权最小二乘法的目的就是让更准确的卫星伪距测量值在定位解算中有着更高的参与度，卫星测量值准确度取决于卫星测量误差。加权最小二乘法通常以卫星伪距测量误差均方差的倒数作为该卫星测量值的权重，即

$W=\left[\begin{array}{cccc}1/{\mathrm{\&sigma;}}_1& & & \\ & 1/{\mathrm{\&sigma;}}_2& & \\ & & ...& \\ & & & 1/{\mathrm{\&sigma;}}_N\end{array}\right].$

根据所述卫星相对飞机几何分布构建观测矩阵G的实施过程如下：

假设某一时刻GPS接收机能接收到N颗可见卫星信号，即产生N个观测量，那么GPS测量方程为

y＝Gx+v

其中y是N×1维观测矢量；x是N×1维待测用户状态矢量，包括用户接收机自身真实的三维位置坐标和接收机相对于卫星的时钟偏差；v是N×1维测量误差矢量，在没有卫星故障时，v向量的各个元素均服从(0，σ²)的高斯分布；G是N×4维观测矩阵，它由3列观测卫星到接收机的方向余弦和全部是1的第四列组成^[18]，表示如下：

$G=\left[\begin{array}{cccc}cos\left({\mathrm{\&theta;}}^{\left(1\right)}\right)cos\left({\mathrm{\&alpha;}}^{\left(1\right)}\right)& cos\left({\mathrm{\&theta;}}^{\left(1\right)}\right)\sin \left({\mathrm{\&alpha;}}^{\left(1\right)}\right)& sin\left({\mathrm{\&theta;}}^{\left(1\right)}\right)& 1\\ cos\left({\mathrm{\&theta;}}^{\left(2\right)}\right)\cos \left({\mathrm{\&alpha;}}^{\left(2\right)}\right)& cos\left({\mathrm{\&theta;}}^{\left(2\right)}\right)\sin \left({\mathrm{\&alpha;}}^{\left(2\right)}\right)& \sin \left({\mathrm{\&theta;}}^{\left(2\right)}\right)& 1\\ ...& ...& ...& ...\\ cos\left({\mathrm{\&theta;}}^{\left(N\right)}\right)\cos \left({\mathrm{\&alpha;}}^{\left(N\right)}\right)& cos\left({\mathrm{\&theta;}}^{\left(N\right)}\right)\sin \left({\mathrm{\&alpha;}}^{\left(N\right)}\right)& sin\left({\mathrm{\&theta;}}^{\left(N\right)}\right)& 1\end{array}\right]$

其中，α^(N)和θ^(N)分别表示可N个见卫星相对于用户接收机的方位角和仰角。

通过最小二乘法可以得到用户状态矢量的最小二乘解为

$\stackrel{\&OverBar;}{x}={\left(G^{T}G\right)}^{-1}{G}^{T}y$

将计算得到的最小二乘解代入方程式

y＝Gx+v

的右侧，并将得到的计算伪距矢量与实际的伪距观测相比较，他们的差值即为伪距残差D，数学上推导如下：

$D-y-\stackrel{\&OverBar;}{y}=y-G\stackrel{\&OverBar;}{x}=y-G{\left(G^{T}G\right)}^{-1}{G}^{T}y$

将方程式

y＝Gx+v

中观测矢量y的表达式代入上，得到

D＝(Gx+v)-G(G^TG)^-1G^T(Gx+v)

＝Gx+v-[Gx+G(G^TG)^-1G^Tv]

＝[I_N-G(G^TG)^-1G^T]v＝Sv

其中，矩阵S定义为

S＝I_N-G(G^TG)^-1G^T。

将所述观测矩阵和权系数矩阵代入加权HPL算法，计算出水平保护值的具体实施过程为：

对于GPS测量方程，通过加权最小二乘法，即分别在观测矩阵G和观测矢量y之前同时乘上权重矩阵W，得到加权最小二乘解为

$\begin{array}{c}\stackrel{\&OverBar;}{x}={\left({\left(WG\right)}^{T}WG\right)}^{-1}{\left(WG\right)}^{T}Wy\\ ={\left(G^{T}CG\right)}^{-1}{G}^{T}Cy\end{array}$

其中，C＝W^TW，W为权重矩阵，表示如下

$W=\left[\begin{array}{cccc}{w}_1& & & \\ & w_2& & \\ & & ...& \\ & & & w_N\end{array}\right]$

因此，伪距残差D_W为

$\begin{array}{c}{D}_W=y-\stackrel{\&OverBar;}{y}=y-G{\left(G^{T}CG\right)}^{-1}{G}^{T}Cy\\ =\&lsqb;I_N-G{\left(G^{T}CG\right)}^{-1}{G}^{T}C\&rsqb;(Gx+v)\\ =\&lsqb;I_N-G{\left(G^{T}CG\right)}^{-1}{G}^{T}C\&rsqb;v\\ =S_{w}v\end{array}$

式中，投影矩阵S_W为

S_w＝I_N-G(G^TCG)^-1G^TC

令A_w＝(G^TCG)^-1G^TC，所以

S_w＝I_N-GA_w

那么，伪距残差平方和F_W为

F_w＝(WD)^TWD＝D^TCD)

＝v^TCS_wv

根据最小二乘法原理，选择统计检测量为

$r_w=\sqrt{F_w/(N-4)}$

同理HPL方法，如果第i颗卫星存在故障，引起伪距测量误差为b_i，那么因此而导致的水平定位误差是

$E_{w,i}=b_i\sqrt{A_{w,1i}^2+A_{w,2i}^2}$

其中，A_W，1i和A_W，2i分别代表矩阵A_W第1行第i列和第2行第i列数值。忽略观测噪声，相应的统计检测量为

$r_{w,i}=b_i\sqrt{S_{w,ii}{C}_{ii}/(N-4)}$

式中，S_W，1i是矩阵S_W的第i行第i列数值，C_ii为矩阵C的第i行第i列元素。

那么该卫星的特征斜率K_i为

$K=E_{w,i}/r_{w,i}=\sqrt{(A_{w,1i}^2+A_{w,2i}^2)(N-4)/S_{w,ii}{C}_{ii}}$

此时，满足漏警率要求的检测门限值r_D’为

$r_D^{\&prime;}=\sqrt{\mathrm{\&lambda;}/(N-4)}$

λ是满足漏检概率要求的非中心分布χ²(N-4，λ)中的非中心参量，同HPL方法一样，非中心参量λ值仅与当前可见卫星数目以及系统要求的虚警率和漏警率性能指标有关。在实际应用中，SA取消后的虚警率和漏警率要求与SA开启时一样。因此，在当时可见卫星数目及几何布局一定时，加权HPL方法中的非中心参量λ值分布与非加权HPL结果一致，如图15中表格内容所示。

根据HPL的定义，最终得出加权HPL的表达式为

$HPL=K_{\max }{r}_D^{\&prime;}=\underset{i}{\max }\sqrt{(A_{w,1i}^2+A_{w,2i}^2)/S_{w,ii}{C}_{ii}}\&times;\sqrt{\mathrm{\&lambda;}}$

从上式中可以看出，加权HPL和非加权HPL结果主要区别在于卫星最大水平斜率不同。在SA政策取消后，卫星自身误差和对流层误差(即信号大气传播延时误差)成为伪距测量误差的主要来源。其中，信号在穿过大气过程中引起的延时误差，与信号传播的角度密切相关，因此相对于接收机，不同观测仰角的卫星信号传播误差不尽相同。加权HPL算法充分考虑到这一特性，不再将卫星伪距测量误差设置为一个固定值，而是根据各个误差校正模型确定相应的误差分量，并由此为各个卫星测量伪距设置了相应的权重系数，与实际GPS定位算法(加权最小二乘法)原理上保持一致，使得RAIM预测结果可信度更高。另外，加权HPL算法中的测量误差整体上要比非加权HPL算法低，所以，理论上其RAIM预测结果，即RAIM可用性更好。这一点在实际应用以及实验仿真中已经得到验证。所以，本发明将加权HPL算法作为RAIM预测的核心算法，

本发明以××机场的RNPAR进近程序作为仿真对象。由于该机场属于高高原机场，地形条件极其复杂，进近程序周围高大山体较多，对于GPS信号遮挡较为严重。具体仿真数据如下表内容所列。

基于WGS-84基准椭球体构建航路大椭圆航线模型，计算得到该进近航迹的四维坐标信息，如下表内容所列。

从上表内容中可以发现，在该RNP进近阶段，随着飞机越来越接近跑道入口，地形遮蔽角越来越大，最大达到8.315°。分别令可见卫星的仰角门限为5°和10°，进行RAIM预测，仿真结果如下表内容所列。

根据各个航段要求的水平告警门限值，如下表内容所示，RNPAR最后进近阶段的HAL值为40m。两次不同可用仰角下的RAIM预测结果显示，5°仰角下，RAIM预测结果均小于40m，即RAIM可用，而在10°仰角下，预测点所在位置的可见卫星数量有所减少，RAIM预测结果显示，部分预测点的HPL值超过40m限值，判定RAIM不可用。

由此可见，RAIM预测算法设置的遮蔽角不同，得到的RAIM预测结果也会不同。通过高程提取，计算得到该进近航段地形遮蔽角最大超过8°，如果不考虑地形遮蔽影响，设置的可见卫星仰角门限为通常意义上的5°，会导致RAIM预测结果与实际不相符，因为飞机在实际进近过程中，部分时间接收不到那些卫星仰角低于8°的卫星测量信号，存在安全隐患。因此，将地形遮蔽角因子纳入到RAIM预测理论是由必要的，尤其是在执行RNP最后进近时，飞行高度低，可操纵余度小，航空器高度依赖GNSS提供水平航迹引导，要求GNSS的定位导航精度绝对可靠，所以RAIM预测过程必须要和实际情况相一致，才能保证RAIM预测结果可信，从而保障飞行安全。

以上对本发明自接收机自主完好性监视的预测方法及预测系统进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述。以上实施过程的说明只是用于帮助理解本发明的核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种接收机自主完好性监视的预测方法，其特征在于，包括以下步骤：

步骤S1，采集三维集地形数据、飞行计划数据、GPS卫星历书数据；

步骤S2，根据所述三维地形数据、所述飞行计划数据和所述GPS卫星历书数据计算出所述飞行计划数据中某个时刻的GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角；

根据所述GPS卫星位置坐标、所述飞机位置坐标和所述航段地形遮蔽角计算卫星相对飞机几何分布；

步骤S3，根据所述地形遮蔽角、所述卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星；

步骤S4，根据卫星相对飞机几何分布分析卫星伪距离测量误差，并根据所述伪距离测量误差构建权系数矩阵；

步骤S5，将所述观测矩阵和权系数矩阵代入加权HPL算法，计算出水平保护值；

步骤S6，将所述水平保护值和水平告警门限值比较，获取所述接收机自主完好性监视结果。

2.根据权利要求1所述的接收机自主完好性监视的预测方法，其特征在于，所述步骤S2中计算飞机位置坐标步骤包括，构建基于性能的导航航线模型，所述基于性能的导航航线模型是基于WGS-84基准椭球体的大椭圆航线；

根据所述大椭圆航线模型以及距离计算算法计算某个时刻的飞机的位置坐标(x、y、z)；

其中，所述大椭圆航线的模型为：

$\left\{\begin{array}{c}px+qy+\mathrm{\&xi;}z=0\\ \frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{b^2}=1\end{array}\right.$

其中，坐标原点为大椭圆航线模型中心，XY平面为赤道平面，其中X轴指向0°经度，Y轴指向东经90°，Z轴指向地理北极极点；A和B为所述大椭圆航线模型中的两点，A点的空间直角坐标记为(x_A，y_A，z_A)，B点的空间直角坐标记为(x_B，y_B，z_B)；a为基准椭球的长半轴，b为基准椭球的短半轴；p、q和ζ满足如下公式：

$p=\left|\begin{array}{cc}{y}_A& z_A\\ y_B& z_B\end{array}\right|,q=\left|\begin{array}{cc}{x}_B& z_B\\ x_A& z_A\end{array}\right|,\mathrm{\&xi;}=\left|\begin{array}{cc}{x}_A& y_A\\ x_B& y_B\end{array}\right|.$

3.据权利要求2所述的接收机自主完好性监视的预测方法，其特征在于，所述步骤S2中地形遮蔽角的计算包括：在所述三维集地形分辨率数据设置高程提取点的经纬度坐标，根据所述经纬度坐标提取相应的高程信息；根据所述高程信息和所述飞机位置坐标计算出某个时刻的地形遮蔽角。

4.根据权利要求2或3所述的接收机自主完好性监视的预测方法，其特征在于，所述步骤S2中计算卫星相对飞机几何分布具体包括：以观测步长对应某个时刻上所述飞机位置为中心，建立站心坐标系，计算某个时刻所述卫星相对于所述飞机的卫星观测仰角、卫星观测方位角和卫星观测距离。

5.根据权利要求4所述的接收机自主完好性监视的预测方法，其特征在于，所述步骤S3中筛选可见卫星包括：根据所述卫星观测仰角、所述地形遮蔽角和接收机门限仰角，筛选除所述卫星观测仰角小于所述接收机门限仰角的卫星和所述卫星观测仰角小于所述地形遮蔽角的卫星。

6.根据权利要求5所述的接收机自主完好性监视的预测方法，其特征在于，根据所述卫星相对飞机几何分布构建观测矩阵G，所述观测矩阵G为：

$G=\left[\begin{array}{cccc}\cos \left({\mathrm{\&theta;}}^{\left(1\right)}\right)\cos \left({\mathrm{\&alpha;}}^{\left(1\right)}\right)& \cos \left({\mathrm{\&theta;}}^{\left(1\right)}\right)\sin \left({\mathrm{\&alpha;}}^{\left(1\right)}\right)& \sin \left({\mathrm{\&theta;}}^{\left(1\right)}\right)& 1\\ \cos \left({\mathrm{\&theta;}}^{\left(2\right)}\right)\cos \left({\mathrm{\&alpha;}}^{\left(2\right)}\right)& \cos \left({\mathrm{\&theta;}}^{\left(2\right)}\right)\sin \left({\mathrm{\&alpha;}}^{\left(2\right)}\right)& \sin \left({\mathrm{\&theta;}}^{\left(2\right)}\right)& 1\\ ...& ...& ...& ...\\ \cos \left({\mathrm{\&theta;}}^{\left(N\right)}\right)\cos \left({\mathrm{\&alpha;}}^{\left(N\right)}\right)& \cos \left({\mathrm{\&theta;}}^{\left(N\right)}\right)\sin \left({\mathrm{\&alpha;}}^{\left(N\right)}\right)& \sin \left({\mathrm{\&theta;}}^{\left(N\right)}\right)& 1\end{array}\right]$

其中，α^(N)和θ^(N)分别表示N个可见卫星相对于飞机的卫星观测方位角和卫星观测仰角。

7.根据权利要求6所述的接收机自主完好性监视的预测方法，其特征在于，所述步骤S4中卫星伪距离测量误差σi满足如下公式：

${\mathrm{\&sigma;}}_i^2={\mathrm{\&sigma;}}_{i,URA}^2+{\mathrm{\&sigma;}}_{i,iono}^2+{\mathrm{\&sigma;}}_{i,tropo}^2+{\mathrm{\&sigma;}}_{i,mp}^2+{\mathrm{\&sigma;}}_{i,tcvr}^2$

其中，σ_i，URA为卫星时钟和星历误差，σ_i，iono为电离层延时误差，σ_i，tropo为对流层延时误差，σ_i，mp为接收机多路径误差，σ_i，tcvr为接收机热噪声误差，i为1到N；

所述权系数矩阵为：

$W=\left[\begin{array}{cccc}1/{\mathrm{\&sigma;}}_1& & & \\ & 1/{\mathrm{\&sigma;}}_2& & \\ & & ...& \\ & & & 1/{\mathrm{\&sigma;}}_N\end{array}\right].$

8.根据权利要求7所述的接收机自主完好性监视的预测方法，其特征在于，所述步骤S5中所述水平保护值的计算公式为：

$HPL=\underset{i}{\max }\sqrt{(A_{w,1i}^2+A_{w,2i}^2)/S_{w,ii}{C}_{ii}}\&times;\sqrt{\mathrm{\&lambda;}}$

其中，A_w，1i和A_W，2i分别代表矩阵A_W第1行第i列和第2行第i列数值；所述矩阵A_W为：

A_w＝(G^TCG)^-1G^TC

S_W，1i是矩阵S_W的第i行第i列数值，矩阵S_w为：

S_w＝I_N-G(G^TCG)^-1G^TC

C_ii为矩阵C的第i行第i列元素，λ是满足漏检概率要求的非中心分布χ²(N-4，λ)中的非中心参量，N为可见卫星的数量，I_N为N阶单位矩阵，矩阵C表示权系数矩阵W的协方差，矩阵C为：

C＝W^TW。

9.一种采用权利要求1至8任一所述预测方法的接收机自主完好性监视的预测系统，其特征在于，包括数据采集模块、第一计算模块、可见卫星筛选模块、权系数矩阵构建模块、水平保护值计算模块和检测结果获取模块；

所述数据采集模块，用于采集三维集地形数据、飞行计划数据、GPS卫星历书数据；

所述第一计算模块，用于根据所述三维地形数据、所述飞行计划数据和所述GPS卫星历书数据计算出所述飞行计划数据中某个时刻的GPS卫星位置坐标、飞机位置坐标和航段地形遮蔽角；根据所述GPS卫星位置坐标、所述飞机位置坐标和所述航段地形遮蔽角计算卫星相对飞机几何分布；

所述可见卫星筛选模块，用于根据所述地形遮蔽角、所述卫星相对飞机几何分布和接收机门限仰角筛选出可见卫星；

所述权系数矩阵构建模块，用于根据卫星相对飞机几何分布分析卫星伪距离测量误差，并根据所述伪距离测量误差构建权系数矩阵；

所述水平保护值计算模块，用于将所述观测矩阵和权系数矩阵代入加权HPL算法，计算出水平保护值；

所述检测结果获取模块，用于将所述水平保护值和水平告警门限值比较，获取所述接收机自主完好性监视结果。

10.根据权利要求9所述的接收机自主完好性监视的预测系统，其特征在于，所述预测系统还包括航班计划调整模块，用于根据所述检测结果获取模块获取的接收机自主完好性监视结果，调整航班计划。