CN113253303B - 一种用于实时监测单频星基增强系统性能的方法 - Google Patents

一种用于实时监测单频星基增强系统性能的方法 Download PDF

Info

Publication number
CN113253303B
CN113253303B CN202110520577.5A CN202110520577A CN113253303B CN 113253303 B CN113253303 B CN 113253303B CN 202110520577 A CN202110520577 A CN 202110520577A CN 113253303 B CN113253303 B CN 113253303B
Authority
CN
China
Prior art keywords
correction
satellite
monitoring station
point
message
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110520577.5A
Other languages
English (en)
Other versions
CN113253303A (zh
Inventor
邵搏
熊帅
原彬
张键
耿永超
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
CETC 20 Research Institute
Original Assignee
CETC 20 Research Institute
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by CETC 20 Research Institute filed Critical CETC 20 Research Institute
Priority to CN202110520577.5A priority Critical patent/CN113253303B/zh
Publication of CN113253303A publication Critical patent/CN113253303A/zh
Application granted granted Critical
Publication of CN113253303B publication Critical patent/CN113253303B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
    • G01S19/08—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing integrity information, e.g. health of satellites or quality of ephemeris data

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Security & Cryptography (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明提供了一种用于实时监测单频星基增强系统性能的方法，确保机场附近区域内SF SBAS服务的可用性，利用机场附近监测站的观测数据、GNSS导航电文以及SF SBAS增强电文，实时解算定位误差和保护级，通过对比定位误差和保护级来实时监测SF SBAS在机场附近的服务性能。本发明提升了SF SBAS增强电文解算效率，优化了SF SBAS增强电文参数使用公式，使得解算流程更为清晰简洁，具有较强的工程实用性，为中国民航SBAS监测与服务系统的建设提供了理论依据和实施思路，解决了飞机在进近和着陆过程中因无法实时解算定位误差而导致的完好性风险，保证了SF SBAS服务的完好性性能。

Description

一种用于实时监测单频星基增强系统性能的方法
技术领域
本发明涉及卫星导航增强技术领域，是单频(Single Frequency,SF)星基增强系统(Satellite Based Augmentation System,SBAS)中一种实时监测SF SBAS服务性能的方法。
背景技术
全球导航卫星系统(Global Navigation Satellite System,GNSS)作为主要导航手段，已经进入了快速发展和应用阶段。为了提升GNSS系统的性能，需要相应增强系统满足不同用户对高完好性和高精度的需求。SBAS由大量分布广泛的监测站对导航卫星进行监测，由地球同步卫星(Geostationary Earth Orbit satellite,GEO)向用户播发改正数信息(星历误差、卫星钟差、电离层延迟)和完好性信息(用户差分距离误差、格网电离层垂直误差)，实现对卫星导航系统定位精度的改进和完好性性能的提高，能够满足从航路、终端区到一类精密进近(CATegory-I,CAT-I)各阶段的导航需求。
北斗全球卫星导航系统(BeiDou Navigation Satellite System,BDS)于2020年7月31日正式建成。作为BDS重要组成部分，北斗星基增强系统(BeiDou Satellite-basedAugmentation System,BDSBAS)按照国际民航组织(International Civil AviationOrganization,ICAO)要求，通过BDS GEO卫星的B1C信号，向中国及周边用户提供SF SBAS服务，播发的增强电文类型如下表所示：
表1 BDSBAS SF SBAS增强电文类型
电文类型 电文内容
0 SF SBAS测试(不可用于生命安全服务)
1 伪随机噪声码(Pseudo Random Noise,PRN)掩码
2-4 快变改正数
6 完好性信息
7 快变改正数降效因子
9 SBAS GEO卫星星历
10 降效参数
17 SBAS GEO卫星历书
18 电离层格网掩码
25 慢变改正数
26 电离层延迟改正数
28 卫星时钟/星历协方差矩阵
63 空白信息
BDSBAS预计2023年前后将向航空用户提供SF SBAS服务，并逐渐取代传统仪表着陆系统成为航空运输的主要导航手段。突发的SF SBAS服务中断将可能给航空安全带来严重的后果。为了满足航空用户对基于SF SBAS飞行进近的高安全性要求，需要对机场附近区域的SF SBAS服务性能进行实时监测，并将实时监测结果发送给机场管制员。如果SF SBAS服务出现异常，机场管制员将及时向机场附近的飞机发出告警，停止使用SF SBAS服务，改用其他导航手段。
目前，国内外尚未有公开文献对SF SBAS服务性能实时监测方法的描述。
发明内容
为了克服现有技术的不足，本发明提供一种用于实时监测单频星基增强系统性能的方法，确保机场附近区域内SF SBAS服务的可用性，利用机场附近监测站的观测数据、GNSS导航电文以及SF SBAS增强电文，实时解算定位误差和保护级，通过对比定位误差和保护级来实时监测SF SBAS在机场附近的服务性能。
本发明解决其技术问题所采用的技术方案的具体步骤如下：
步骤1：数据预处理；
利用监测站观测到卫星的双频伪距观测量和双频载波相位观测量进行数据预处理，预处理的步骤为：
(1)利用前5个采样时刻的载波观测量外推出当前时刻的观测量，并将观测量与当前时刻接收机的载波相位观测量相减并取绝对值，如果该绝对值小于0.055，则认为没有周跳出现，进行后续计算；
(2)利用载波相位观测量对伪距观测量进行平滑，得到平滑后的伪距观测量；
步骤2：导航电文解算；
利用导航电文播发的轨道参数和时钟参数计算卫星星历位置和卫星时钟偏差，利用卫星星历位置和监测站位置计算星历距离；
步骤3：增强电文解算；
步骤3.1：利用SF SBAS增强电文2、3、4中的快变改正信息计算当前时刻的快变改正数；
步骤3.2：利用增强电文25中的慢变改正数信息计算卫星时钟慢变改正数和卫星星历慢变改正数；
步骤3.3：利用增强电文18和26中的电离层格网点信息，推算监测站所在位置电离层改正数；
步骤3.4：基于增强电文2、3、4解算用户差分距离误差，基于增强电文26解算格网电离层垂直误差；
步骤3.5：利用增强电文28中的协方差矩阵信息，解算用户差分距离误差降效降效参数；
步骤3.6：利用增强电文7、电文10、电文25中的相关信息，解算快慢变改正降效参数；
步骤3.7：利用增强电文10和26中的相关信息，解算电离层改正降效参数；
步骤4：定位解算；
(1)基于对流层干湿分量模型，利用监测站所处纬度计算对流层延迟估计；
(2)平滑后的伪距观测量经快变改正数、电离层改正数修正，并消除对流层延迟和卫星时钟偏差后，得到修正后的伪距，对修正后的伪距进行一阶泰勒级数展开，通过迭代法解算监测站位置；
(3)将地心地固坐标系下的定位误差转换到东北天坐标系下，并求解水平和垂直定位误差。
步骤5：保护级解算。利用观测矩阵和观测伪距噪声方差解算水平和垂直保护级。
步骤6：服务性能评估。通过比对定位误差和保护级判断SF SBAS实时服务性能，如果SF SBAS服务出现异常，机场管制员将及时向机场附近的飞机发出告警，停止使用SFSBAS服务。
所述数据预处理的步骤为：
机场附近的监测站采集所监测到全球卫星导航系统卫星的观测数据、GNSS导航电文和SF SBAS增强电文，监测站i观测到卫星j的双频观测数据如下：
其中，和/>分别为L1和L5频点上的伪距观测量；/>和/>分别为L1和L5频点上的载波相位观测量；/>为监测站i和卫星j间的几何距离；/>为对流层延迟；bi为监测站接收机时钟与GNSS系统时之间的偏差；Bj为卫星时钟与GNSS系统时之间的偏差；/>为电离层延迟，对伪距观测量的影响是滞后，对载波相位观测量的影响是超前；/>f1＝1575.42MHz为载波L1的频率，f5＝1176.45MHz为载波L5的频率；/>和/>为伪距观测量上的观测噪声；N1和N5为整周模糊度，由接收机失锁造成；λ1＝C/f1和λ5＝C/f5分别为载波L1和L5的波长，光速C＝299792458m/s；/>和/>为载波相位观测量上的观测噪声，该噪声远远小于伪距观测量上的观察噪声。不同时刻的数据会进行标识，未做说明的数据均为t时刻的数据；
利用监测站i(i＝1，2，…，M)观测到卫星j的双频伪距观测量和双频载波相位观测量进行数据预处理，具体步骤如下：
步骤1.1：周跳探测
周跳探测利用前5个采样时刻(t-1,t-2,t-3,t-4,t-5)的载波观测量外推出当前时刻的观测量，并将此观测量与当前时刻接收机的载波相位观测量比较，如果超出门限则认为出现周跳；监测站i对卫星j的周跳监测有下式：
式(5)和式(6)中，α0、α1、a2为拟合系数，[a0，a1，a2]T＝(FTF)-1FTXL1-L5，F取固定值；/>为t时刻观测量组合，/>和/>分别为t时刻L1和L5频点上的载波相位观测量；/>为由多项式拟合得到的t时刻观测量组合的拟合值；TL1-L5＝0.055为探测门限；
步骤1.2：载波平滑
载波相位观测量通过周跳探测后，则认为没有周跳出现，利用载波相位观测量对伪距观测量进行平滑，首先对载波观测量进行如下变化：
式(9)中
由于前后两个时刻的整周模糊度基本相同，用/>平滑伪距观测量中的噪声；
其中，为t时刻L1频点的伪距观测量，/>为L1频点平滑后的伪距观测量，τ＝100为平滑时间。
所述步骤1.2中，如果L1和L5频点同时出现相同的周跳时，无法检测出周跳，则只利用L5相位观测量进行单频点周跳探测，利用公式(7)、(8)再检测一次：
其中，b0、b1、b2为拟合系数，[b0，b1，b2]T＝(FTF)-1FTXL5；TL5＝0.35为探测门限；
所述步骤二的导航电文解算中，GNSS卫星导航电文中播发的轨道参数为：星历参考时间tce，卫星轨道长半轴as的平方根，轨道偏心率es，toe时刻的轨道倾角i0，周内时等于0时的轨道升交点赤经Ω0，轨道近地角距ω，toe时刻的平近点角M0，平均运动角速度校正值Δn，轨道倾角变化率i’，轨道升交点赤经变化率升交点角距余弦调和校正振幅Cuc，升交点角距正弦调和校正振幅Cus，轨道半径余弦调和校正振幅Crc，轨道半径正弦调和校正振幅Crs，轨道倾角余弦调和校正振幅Cic，轨道倾角正弦调和校正振幅Cis；利用导航电文播发的轨道参数及卫星位置解算通用算法(此处不再赘述)，得到卫星星历位置/>利用卫星星历位置/>和监测站基准位置[xR，i，yR，i，zR，i]计算星历距离
利用导航电文中播发的参考时间toe、参考时刻的卫星时钟偏差af0、卫星时钟漂移速度af1和卫星时钟漂移速度的变化率af2计算t时刻的卫星时钟偏差
所述步骤三的增强电文解算的详细步骤为：
步骤3.1：快变改正数解算
快变改正数信息由增强电文2、3、4播发，计算公式如下：
PRC(t)＝PRC(tof)+RRC(tof)×(t-tof) (11)
如果aii≠0(由电文7播发)，则
其中，PRCcurrent为最新接收到的快变改正数(由电文2、3、4播发)；PRCprevious为PRCcurrent之前接收到的快变改正数(由电文电文2、3、4播发)；Δt＝(tof-tof，previous)；tof为PRCcurrent的参考时刻，为SBAS系统时(SBASNetworkTime，SNT)的整秒开始，SBAS GEO卫星播发电文第一比特的时间；tof，previous为PRCprevious的参考时刻；
如果aii＝0，则RRC＝0；
步骤3.2：慢变改正数解算
慢变改正数信息由增强电文25播发，慢变改正数信息包括卫星时钟慢变改正数和卫星星历慢变改正数信息；
卫星时钟慢变改正数通过下式计算：
δΔtSV(t)＝δaf0+δaf1(t-t0) (13)
其中，δaf0为时钟偏差(电文25)；δaf1为时钟偏差变化率，电文25，速度标识为1；如果速度标识为0该值为0；t0为改正数参考时刻，电文25，速度模式标识为1；
卫星星历慢变改正数通过式(14)计算，采用坐标系WGS-84 ECEF：
其中，[δxk δyk δzk]为时刻t对应的慢变改正数，[δx δy δz]为时刻t0对应的慢变改正数(由电文25播发)，为时刻t0慢变改正数的变化率，由电文25播发；
当速度标识为0时，上式中的速度分量为0；
步骤3.3：电离层改正数解算；
步骤3.3.1：穿刺点经纬度解算
首先计算穿刺点纬度φpp，公式如(15)和(16)所示：
φpp＝sin-1(sinφu cosψpp+cosφu sinψpp cos A) (15)
式中，ψpp为地球中心角；A为方位角；E为高度角；Re为近似地球半径；hI为最大电子密度层的近似高度；
穿刺点经度λpp的计算公式如(17)和(18)所示：
(1)当φu＞70°且tanψppcosA＞tan(π/2-φu)或者φu＜-70°且
tanψpp cos(A+π)＞tan(π/2+φu)时，λpp为：
(2)其他情况下，λpp为：
步骤3.3.2：格网点选择
计算完穿刺点经纬度后，监测站根据接收到的电文18确定使用的电离层格网点，电离层格网点的选择与穿刺点所在的纬度有关，选取策略如下：
(1)穿刺点位于南北纬60°之间
如果穿刺点的周围有4个格网点构成5°×5°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则进入下一步；
如果穿刺点的周围有3个格网点构成5°×5°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则进入下一步；
如果穿刺点的周围有4个格网点构成10°×10°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则进入下一步；
如果穿刺点的周围有3个格网点构成10°×10°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则认为电离层改正数不可用；
(2)穿刺点位于北纬60°～75°之间或者南纬60°～75°之间
如果穿刺点的周围有4个格网点构成5°(纬度)×10°(经度)的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则进入下一步；
如果穿刺点的周围有3个格网点构成5°(纬度)×10°(经度)的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则进入下一步；
如果穿刺点的周围有4个格网点构成10°×10°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则进入下一步；
如果穿刺点的周围有3个格网点构成10°×10°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则认为电离层改正数不可用。
(3)穿刺点位于北纬75°～85°之间或者南纬75°～85°之间
如果两个最近的格网点位于75°，两个最近的格网点位于85°且对应的IGP掩码为1时，通过线性内插得到虚拟格网点构建出10°×10°的四边形区域，则选取这四个网格点，否则认为电离层改正数不可用；
(4)穿刺点位于北纬85°
如果四个格网点位于北纬85°，西经180°、西经90°、0°、东经90°，且对应的IGP掩码为1时，则选取这四个格网点，否则认为电离层改正数不可用；
(5)穿刺点位于南纬85°
如果四个格网点位于南纬85°，西经140°、西经50°、东经40°、东经130°，且对应的IGP掩码为1时，则选取这四个格网点，否则认为电离层改正数不可用；
步骤3.3.3：穿刺点电离层垂直延迟解算
由于增强电文播发的是格网点处的电离层垂直延迟，监测站需要利用周围的格网点来计算穿刺点处的电离层延迟；
(1)周边4个格网点
如果在穿刺点周围选取了4个格网点，穿刺点处的电离层垂直延迟τvpp(φpp，λpp)计算公式如(19)所示：
式中，φpp和λpp分别为穿刺点的纬度和经度；τvi为格网点处的电离层垂直延迟；Wi(xpp，ypp)为格网点的权重函数；
W1(xpp，ypp)＝xppypp (20)
W2(xpp，ypp)＝(1-xpp)ypp (21)
W3(xpp，ypp)＝(1-xpp)(1-ypp) (22)
W4(xpp，ypp)＝xpp(1-ypp) (23)
Δλpp＝λpp-λ1 (24)
Δφpp＝φpp-φ1 (251
如果穿刺点位于北纬85°和南纬85°之间，xpp和ypp计算公式如(26)和(27)所示：
式中，λ1为穿刺点西边格网点的经度；λ2为穿刺点东边格网点的经度；φ1为穿刺点南边格网点的纬度；φ2为穿刺点北边格网点的纬度；
如果位于北纬85°以北或南纬85°以南，xpp和ypp计算公式如(28)和(29)所示：
式中，λ1为穿刺点东边第二个格网点的经度；λ2为穿刺点西边第二个格网点的经度；λ3为穿刺点西边第一个格网点的经度；λ4为穿刺点东边第一个格网点的经度；
监测站利用格网点的计算穿刺点处的/>公式如(30)所示：
(2)周边3个格网点
如果在穿刺点周围选取了3个格网点且位于北纬75°和南纬75°之间，穿刺点处的电离层垂直延迟τvpp(φpp，λpp)计算公式如(31)所示：
式中，φpp和λpp分别为穿刺点的纬度和经度；τvi为格网点处的电离层垂直延迟；Wi(xpp，ypp)为格网点的权重函数；
W1(xpp，ypp)＝ypp (32)
W2(xpp，ypp)＝1-xpp-ypp (33)
W3(xpp，ypp)＝xpp (34)
监测站利用格网点的计算穿刺点处的/>公式如(35)所示：
步骤3.3.4：穿刺点电离层改正数解算
穿刺点处的电离层改正数计算公式如(36)所示：
式中，τspp(φpp，λpp)为穿刺点处电离层倾斜延迟，Fpp为倾斜因子；
穿刺点处的为：
步骤3.4：完好性参数解算
完好性参数包括用户差分距离误差(User Differential Range Error，UDRE)和格网电离层垂直误差(Grid Ionospheric Vertical Error，GIVE)；
UDRE以户差分距离误差索引(UserDifferentialRangeErrorIndicator，UDREI)的形式通过增强电文2、3、4播发，UDREI与UDRE之间的转换关系固定；
GiVE以格网电离层垂直误差索引(GridIonosphericVerticalErrorIndicator，GIVEI)的形式通过增强电文26播发，GIVEI与GIVE之间的转换关系固定；
步骤3.5：UDRE降效参数解算
协方差矩阵信息由增强电文28播发，利用电文中的信息计算上三角矩阵R：
R＝SF·E (39)
其中，SF＝2(scal exponent-5)，
利用矩阵R计算协方差矩阵C：
C＝RT·R (40)
投影到监测站上的UDRE降效参数：
其中，I为卫星到监测站的4维方向矢量，前三维是单位方向矢量，第四维是1；εc＝Ccovariance·SF，Ccovariance由电文10播发，如果电文10中的Ccovariance无效，则该值为0；
步骤3.6：快慢变改正降效参数解算
通过快变改正数和慢变改正数修正后的残差方差通过快慢变降效参数表征，计算公式如下：
其中，RSSUDRE由增强电文10播发；σUDRE按照步骤3.4中的方法计算；δUDRE按照步骤3.5中的方法计算；εfc为快变改正降效参数；εrrc为距离变化改正降效参数；εltc为慢变改正降效参数；εer为应用模式降效参数；
(1)快变改正降效参数解算
快变改正降效参数εfc的计算公式如下：
εfc＝a(t-tu+tlat)2/2 (43)
其中，a为快变改正数降效因子，由增强电文7播发，t为当前时刻，tu为快变改正数的参考时间，如果使用的UDREI是由增强电文2、3、4播发的，则tu为快变改正数的参考时间；如果UDREI是由增强电文6播发并且IODF＝3，则tu也为快变改正数的参考时间；如果IODF≠3，则tu为GEO播发增强电文6第一比特的时间；tlat为系统延迟时间，由增强电文7播发；
(2)距离变化改正降效参数
如果电文7播发的aii＝0，则距离变化改正降效参数εrrc＝0，否则，εrrc将分为两种情况进行计算：
a)当前最新的IODF和前一个IODF都不为3，则：
b)当前最新的IODF和前一个IODF至少有一个为3，则：
计算过程中使用的参数有：a为快变改正数的降效因子，由增强电文7播发；t为当前时刻；Ifc，j为增强电文2、3、4的最短超时间隔，根据电文7播发的a来推算，Brrc由增强电文10播发；IODFcurrent为最新接收的快变改正数信息中的IODF；IODFDrevious为在收到最新信息之前收到的快变改正数信息中的IODF；tof为最新接收的快变改正数信息的参考时间；tof，previous为收到最新信息之前收到的快变改正数信息的参考时间；Δt＝tof-tof，previous；
(3)慢变改正降效参数
与慢变改正相关的降效参数根据电文中是否同时包含偏差和偏差变化量，分为两种情况：
1)增强电文25中速度模式标识为1；
2)增强电文25中速度模式标识为0；
当速度模式标识为1，则：
其中，t为当前时刻；t0为长期改正数的参考时间，由增强电文25播发；Iltc_v1为更新间隔，由增强电文10播发；Cltc_lsb由增强电文10播发；Cltc_v1由增强电文10播发；
当速度模式标识为0，则：
其中，t为当前时刻；tltc为GEO播发长期改正数电文第一比特的时间；Iltc_v0为最小更新间隔，由增强电文10播发；Cltc_v0由增强电文10播发；表示向下取整；
(4)应用模式降效参数
在侧向导航/垂直导航(Lateral NAVigation/Vertical NAVigation，LNAV/VNAV)，导航信标性能(Localizer Performance，LP)，垂直导航信标性能(LocalizerPerformancewithVertical guidance，LPV)等飞行阶段，如果快变改正数或慢变改正数超时，但未超出航路至LNAV的时间限制，需要适用于航路至LNAV的降效参数；
适用于航路至LNAV的应用模式降效参数εer的计算公式如下：
其中，Cer由电文10播发；
步骤3.7：电离层改正降效参数解算
电离层改正降效参数的计算公式如下：
其中，RSSiono由增强电文10播发；σGIVE按照步骤3.4中的方法计算；t为当前时刻；tiono为GEO卫星播发第一比特电离层改正数信息的时间；Ciono_step由增强电文10播发；Ciono_ramp由增强电文10播发；Iiono为电离层改正数信息的最小更新间隔，由增强电文10播发；表示向下取整。
所述步骤四：定位解算；
步骤4.1：对流层延迟估计
对流层模型进行修正，对流层延迟估计计算如下：
其中，dhyd与dwet分别表示对流层的干分量和湿分量， 为仰角，/>b＝acos[cos(φj-φi)×cos(δj-δi)]，φj和δj分别为卫星j所在位置的纬度和经度，φi和δi分别为监测站i所在位置的纬度和经度；
dhyd与dwet由监测站高度信息及五个气象参数的估值计算：
/>
其中，g＝9.80665m/s2，gm＝9.784m/s2，H为监测站海拔，k1＝77.604K/mbar，k2＝382000K2/mbar，Rd＝287.054J/kg/K；
气象参数气压P(mbar)、温度T(K)、水汽压e(mbar)、温度变化率β(k/m)、水汽变化率λ由监测站的气象传感器提供，也可由当前观测站所在纬度φi和年积日D(当年1月1日起开始计算的天数)插值计算，插值公式如下：
其中，Dmin＝28(φi为北纬)，Dmin＝211(φi为南纬)，ξ0和Δξ分别表示不同纬度的气象参数平均值和季节变化值，计算P、T、e、β、λ时，将ξ分别替换为P、T、e、β、λ，按表3由下式插值得到：
如果φi≤15或φi＝30或φi＝45或φi＝60或φi≥75，直接利用ξ0(φi)和Δξ(φi)通过式(32)计算，其他情况下，以φi＝40为例，对应的φk＝30，φk+1＝45，利用ξ0(φk+1)、ξ0(φk)、Δξ(φk+1)和Δξ(φk)通过式(33)和(34)计算ξ0(φi)和Δξ(φi)；
步骤4.2：监测站位置解算；
卫星j星历位置经卫星星历慢变改正数/>改正后的位置/>如下：
卫星j时钟偏差经卫星时钟慢变改正数/>改正后的时钟偏差/>如下：
其中，C为光速；
平滑后的伪距观测量经快变改正数/>电离层改正数/>修正，并消除对流层延迟/>和卫星时钟偏差/>得到伪距/>如下：
/>
其中，[xi，yi，zi]为监测站位置，ti为监测接收机时钟偏差。
该伪距方程为非线性方程，要用泰勒级数展开，并取一阶量，将伪距方程转化为线性方程：
其中，
为监测站位置估计值，[Δxi Δyi Δzi]为监测站位置与估计值之间的差值，/>为用户时钟偏差估计值，Δti为监测接收机时钟偏差与估计值之间的差值；
将式(62)进行变换得到：
其中，
式(63)对应的矩阵形式为：
Z＝HX (51)
其中，
X＝[Δxi Δyi Δzi-C·Δti]T，N为监测站观测到的卫星数量；
利用最小二乘法得到：
X＝(HTH)-1HTZ (52)
则监测站位置和时钟偏差为：
取为[xi yi zi]，通过多次迭代，当/>时，可以得到监测站位置和时钟偏差；
步骤4.3：监测站定位误差解算
利用步骤4.2得到的监测站位置[xi yi zi]，并结合测绘标定得到的真实位置[xR，i，yR，i，zR，i]，得到地心地固(Earth-Centered Earth-Fixed，ECEF)坐标系下的监测站定位误差：
[ΔxR，i ΔyR，i ΔzR，i]＝[xR，i，yR，i，zR，i]-[xi yizi] (55)
从ECEF坐标系到东北天(EastNorthUp，ENU)坐标系的转换矩阵为：
/>
其中，φi和λi分别是监测站所处位置的地理纬度和经度；
得到ENU坐标系下的定位误差如下：
[ΔEi ΔNi ΔUi]＝Pi·[ΔxR，i ΔyR，i ΔzR，i]T (57)
基于式(46)得到如下：
VPE＝ΔUi (72)
得到水平定位误差(Horizontal Position Error，HPE)和垂直定位误差(VerticalPosition Error，VPE)。
所述步骤五的保护级解算的步骤为：
首先计算监测站与可观测卫星间的观测矩阵G，该矩阵的第j行如下所示：
其中，为监测站与卫星j间的仰角；/>为监测站与卫星j间的方位角；
监测站与第j个可观测的卫星间观测伪距的噪声方差为：
其中，按照步骤3.6中的方法计算；/>按照步骤3.3中的方法计算；
按照步骤4.1中的方法计算；
监测站与卫星间观测伪距的协方差矩阵为W，其对角线元素其余元素为0，通过G和W得到：
水平保护级HPL可由下式计算：
HPL＝KH，NPA·dmajor (60)
其中，KH，NPA＝6.18；
针对一类垂直引导进近(APproachwithVertical guidance I，APV-I)、二类垂直引导进近(APV-II)和CAT-I等飞行阶段，HPL和垂直保护级(Vertical Protection Level，VPL)由下式计算：
HPL＝KH，PA·amajor (77)
VPL＝KVdU (78)
其中，KH，PA＝6.0，KV＝5.33。
步骤六：服务性能评估；
对于针对航路、终端区、NPA等飞行阶段，如果HPE≤HPL，表明SF SBAS服务正常；如果HPE＞HPL，表明SF SBAS服务不能用于导航；对于APV-I、APV-II和CAT-I等飞行阶段，如果HPE≤HPL且VPE≤VPL，表明SF SBAS服务正常；如果HPE＞HPL或者VPE＞VPL，则SF SBAS服务不能够用于引导飞机进行精密进近；当SF SBAS服务不可用时，机场管制员需要将SF SBAS服务不可用的信息告知机场附近准备进近着陆的飞机，飞机需要采用其他导航手段进行进近。
本发明的有益效果在于：
1)提升了SF SBAS增强电文解算效率，将一些不需要使用的电文参数进行了删减，优化了SF SBAS增强电文参数使用公式，使得解算流程更为清晰简洁。
2)提出了实时监测机场附近SF SBAS服务性能的方法，给出了明确的处理流程和实施步骤，具有较强的工程实用性，为中国民航SBAS监测与服务系统的建设提供了理论依据和实施思路；
3)利用机场附近的监测站作为参考基准，实时解算定位误差和保护级，评估SFSBAS在机场的实时服务性能，解决了飞机在进近和着陆过程中因无法实时解算定位误差而导致的完好性风险，保证了SF SBAS服务的完好性性能。
附图说明
图1是本发明SF SBAS服务性能实时监测步骤流程图。
图2是本发明电离层穿刺点示意图。
图3是本发明4个格网点与穿刺点分布示意图。
图4是本发明3个格网点与穿刺点分布示意图。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
本发明是一种用于实时监测单频星基增强系统性能的方法，具体步骤如图1所示：
步骤一：数据预处理
机场附近的监测站采集所监测到全球卫星导航系统(Global NavigationSatellite System，GNSS)卫星的观测数据、GNSS导航电文和SF SBAS增强电文，监测站i观测到卫星j的双频观测数据如下：
其中，和/>分别为L1和L5频点上的伪距观测量；/>和/>分别为L1和L5频点上的载波相位观测量；/>为监测站i和卫星j间的几何距离；/>为对流层延迟；bi为监测站接收机时钟与GNSS系统时之间的偏差；Bj为卫星时钟与GNSS系统时之间的偏差；/>为电离层延迟，对伪距观测量的影响是滞后，对载波相位观测量的影响是超前；/>f1＝1575.42MHz为载波L1的频率，f5＝1176.45MHz为载波L5的频率；/>和/>为伪距观测量上的观测噪声；N1和Ns为整周模糊度，由接收机失锁造成；λ1＝C/f1和λ5＝C/f5分别为载波L1和L5的波长，光速C＝299792458m/s；/>和/>为载波相位观测量上的观测噪声，该噪声远远小于伪距观测量上的观察噪声。不同时刻的数据会进行标识，未做说明的数据均为t时刻的数据。
利用监测站i(i＝1，2，…，M)观测到卫星j的双频伪距观测量和双频载波相位观测量进行数据预处理，具体步骤如下：
步骤1.1：周跳探测
周跳探测利用前5个采样时刻(t-1，t-2，t-3，t-4，t-5)的载波观测量外推出当前时刻的观测量，并将此观测量与当前时刻接收机的载波相位观测量比较，如果超出门限则认为出现周跳；监测站i对卫星j的周跳监测有下式：
式(5)和式(6)中，a0、a1、a2为拟合系数，[a0，a1，a2]T＝(FTF)-1FTXL1-L5，F取固定值；/>为t时刻观测量组合，/>和/>分别为t时刻L1和L5频点上的载波相位观测量；/>为由多项式拟合得到的t时刻观测量组合的拟合值；TL1-L5＝0.055为探测门限；
如果L1和L5频点同时出现相同的周跳时，无法检测出周跳，则只利用L5相位观测量进行单频点周跳探测，利用公式(7)、(8)再检测一次：
其中，b0、b1、b2为拟合系数，[b0，b1，b2]T＝(FTF)-1FTXL5；Tl5＝0.35为探测门限；
步骤1.2：载波平滑
载波相位观测量通过周跳探测后，则认为没有周跳出现，利用载波相位观测量对伪距观测量进行平滑，首先对载波观测量进行如下变化：
式(9)中
由于前后两个时刻的整周模糊度基本相同，用/>平滑伪距观测量中的噪声；
其中，为t时刻L1频点的伪距观测量，/>为L1频点平滑后的伪距观测量，τ＝100为平滑时间；
步骤二：导航电文解算
GNSS卫星导航电文中播发的轨道参数为：星历参考时间toe，卫星轨道长半轴as的平方根，轨道偏心率es，toe时刻的轨道倾角i0，周内时等于0时的轨道升交点赤经Ω0，轨道近地角距ω，toe时刻的平近点角M0，平均运动角速度校正值Δn，轨道倾角变化率i′，轨道升交点赤经变化率升交点角距余弦调和校正振幅Cuc，升交点角距正弦调和校正振幅Cus，轨道半径余弦调和校正振幅Crc，轨道半径正弦调和校正振幅Crs，轨道倾角余弦调和校正振幅Cic，轨道倾角正弦调和校正振幅Cis；利用导航电文播发的轨道参数及卫星位置解算通用算法(此处不再赘述)，得到卫星星历位置/>利用卫星星历位置/>和监测站基准位置[xR，i，yR，i，zR，i]计算星历距离/>
利用导航电文中播发的参考时间toe、参考时刻的卫星时钟偏差af0、卫星时钟漂移速度af1和卫星时钟漂移速度的变化率af2计算t时刻的卫星时钟偏差
步骤三：增强电文解算
步骤3.1：快变改正数解算
快变改正数信息由增强电文2、3、4播发，计算公式如下：
PRC(t)＝PRC(tof)+RRC(tof)×(t-tof) (71)
如果aii≠0(由电文7播发)，则
其中，PRCcurrent为最新接收到的快变改正数(由电文2、3、4播发)；PRCprevious为PRCcurrent之前接收到的快变改正数(由电文电文2、3、4播发)；Δt＝(tof-tof，previous)；tof为PRCcurrent的参考时刻，为SBAS系统时(SBAS Network Time，SNT)的整秒开始，SBAS GEO卫星播发电文第一比特的时间；tof，previous为PRCprevious的参考时刻。
如果aii＝0，则RRC＝0。
步骤3.2：慢变改正数解算
慢变改正数信息由增强电文25播发，慢变改正数信息包括卫星时钟慢变改正数和卫星星历慢变改正数信息。
卫星时钟慢变改正数通过下式计算：
δΔtSV(t)＝δaf0+δαf1(t-t0) (73)
其中，δaf0为时钟偏差(电文25)；δaf1为时钟偏差变化率(电文25，速度标识为1；如果速度标识为0该值为0)；t0为改正数参考时刻(电文25，速度模式标识为1)。
卫星星历慢变改正数通过下式计算(坐标系WGS-84ECEF)：
其中，[δxk δyk δzk]为时刻t对应的慢变改正数，[δx δy δz]为时刻t0对应的慢变改正数(由电文25播发)，为时刻t0慢变改正数的变化率(由电文25播发)。
当速度标识为0时，上式中的速度分量为0。
步骤3.3：电离层改正数解算
步骤3.3.1：穿刺点经纬度解算
电离层穿刺点示意图如图2所示，首先计算穿刺点纬度φpp(单位：弧度)，公式如(15)和(16)所示：
φpp＝sin-1(sinφu cosψpp +cosφu sinψpp cos A) (75)
式中，ψpp为地球中心角(单位：弧度)；A为方位角；E为高度角；Re为近似地球半径，计算中取6378.1363km；hI为最大电子密度层的近似高度，计算中取350km。
穿刺点经度λpp(单位：弧度)的计算公式如(17)和(18)所示：
(1)当φu＞70°且tanψppcosA＞tan(π/2-φu)或者φu＜-70°且
tanψppcos(A+π)＞tan(π/2+φu)时，λpp为：
(2)其他情况下，λpp为：
步骤3.3.2：格网点选择
计算完穿刺点经纬度后，监测站需要根据接收到的电文18确定使用的电离层格网点。电离层格网点的选择与穿刺点所在的纬度有关，选取策略如下：
(1)穿刺点位于南北纬60°之间
如果穿刺点的周围有4个格网点构成5°×5°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则采用下面的选择；
如果穿刺点的周围有3个格网点构成5°×5°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则采用下面的选择；
如果穿刺点的周围有4个格网点构成10°×10°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则采用下面的选择；
如果穿刺点的周围有3个格网点构成10°×10°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则认为电离层改正数不可用。
(2)穿刺点位于北纬60°～75°之间或者南纬60°～75°之间
如果穿刺点的周围有4个格网点构成5°(纬度)×10°(经度)的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则采用下面的选择；
如果穿刺点的周围有3个格网点构成5°(纬度)×10°(经度)的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则采用下面的选择；
如果穿刺点的周围有4个格网点构成10°×10°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则采用下面的选择；
如果穿刺点的周围有3个格网点构成10°×10°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则认为电离层改正数不可用。
(3)穿刺点位于北纬75°～85°之间或者南纬75°～85°之间
如果两个最近的格网点位于75°，两个最近的格网点位于85°(使用带9或10的话，经度相距30°，否则相距90°)且对应的IGP掩码为1时，通过线性内插得到虚拟格网点构建出10°×10°的四边形区域，则选取这四个网格点(虚拟网格点)，否则认为电离层改正数不可用。
(4)穿刺点位于北纬85°
如果四个格网点位于北纬85°，西经180°、西经90°、0°、东经90°，且对应的IGP掩码为1时，则选取这四个格网点，否则认为电离层改正数不可用。
(5)穿刺点位于南纬85°
如果四个格网点位于南纬85°，西经140°、西经50°、东经40°、东经130°，且对应的IGP掩码为1时，则选取这四个格网点，否则认为电离层改正数不可用。
步骤3.3.3：穿刺点电离层垂直延迟解算
由于增强电文播发的是格网点处的电离层垂直延迟，监测站需要利用周围的格网点来计算穿刺点处的电离层延迟。
(1)周边4个格网点
如果在穿刺点周围选取了4个格网点，格网点与穿刺点分布如图3所示。穿刺点处的电离层垂直延迟τvpp(φpp，λpp)计算公式如(19)所示：
式中，φpp和λpp分别为穿刺点的纬度和经度；τvi为格网点处的电离层垂直延迟；Wi(xpp，ypp)为格网点的权重函数。
W1(xpp，ypp)＝xppypp (80)
W2(xpp，ypp)＝(1-xpp)ypp (81)
W3(xpp，ypp)＝(1-xpp)(1-ypp) (82)
W4(xpp，ypp)＝xpp(1-ypp) (83)
Δλpp＝λpp-λ1 (84)
Δφpp＝φpp-φ1 (85)
如果穿刺点位于北纬85°和南纬85°之间，xpp和ypp计算公式如(26)和(27)所示：
式中，λ1为穿刺点西边格网点的经度；λ2为穿刺点东边格网点的经度；φ1为穿刺点南边格网点的纬度；φ2为穿刺点北边格网点的纬度。
如果位于北纬85°以北或南纬85°以南，xpp和ypp计算公式如(28)和(29)所示：
式中，λ1为穿刺点东边第二个格网点的经度；λ2为穿刺点西边第二个格网点的经度；λ3为穿刺点西边第一个格网点的经度；λ4为穿刺点东边第一个格网点的经度。
监测站利用格网点的计算穿刺点处的/>公式如(30)所示：/>
(2)周边3个格网点
如果在穿刺点周围选取了3个格网点且位于北纬75°和南纬75°之间，格网点与穿刺点分布如图4所示。穿刺点处的电离层垂直延迟τvpp(φpp，λpp)计算公式如(31)所示：
式中，φpp和λpp分别为穿刺点的纬度和经度；τvi为格网点处的电离层垂直延迟；Wi(xpp，ypp)为格网点的权重函数。
W1(xpp，ypp)＝ypp (92)
W2(xpp，ypp)＝1-xpp-ypp (93)
W3(xpp，ypp)＝xpp (94)
监测站利用格网点的计算穿刺点处的/>公式如(35)所示：
步骤3.3.4：穿刺点电离层改正数解算
穿刺点处的电离层改正数计算公式如(36)所示：
式中，τspp(φpp，λpp)为穿刺点处电离层倾斜延迟，Fpp为倾斜因子。
穿刺点处的为：
步骤3.4：完好性参数解算
完好性参数主要包括用户差分距离误差(User Differential Range Error，UDRE)和格网电离层垂直误差(Grid Ionospheric Vertical Error，GIVE)。
UDRE以户差分距离误差索引(UserDifferentialRangeErrorIndicator，UDREI)的形式通过增强电文2、3、4播发。UDREI与UDRE之间的转换关系是固定的，如表1所示：
表1 UDREI与UDRE转换表
GIVE以格网电离层垂直误差索引(Grid Ionospheric Vertical ErrorIndicator，GIVEI)的形式通过增强电文26播发。GIVEI与GIVE之间的转换关系是固定的，如表2所示：
表2 GIVEI与GIVE转换表
步骤3.5：UDRE降效参数解算
协方差矩阵信息由增强电文28播发，利用电文中的信息计算上三角矩阵R：
R＝SF·E (99)
其中，SF＝2(scaleexponent-5)，
利用矩阵R计算协方差矩阵C：
C＝RT·R (100)
投影到监测站上的UDRE降效参数：
其中，I为卫星到监测站的4维方向矢量(前三维是单位方向矢量，第四维是1)；εc＝Ccovariance·SF，Ccovariance由电文10播发，如果电文10中的Ccovariance无效，则该值为0。
步骤3.6：快慢变改正降效参数解算
通过快变改正数和慢变改正数修正后的残差方差通过快慢变降效参数表征，计算公式如下：
其中，RSSUDRE由增强电文10播发；σUDRE按照步骤3.4中的方法计算；δUDRE按照步骤3.5中的方法计算；εfc为快变改正降效参数；εrrc为距离变化改正降效参数；εltc为慢变改正降效参数；εer为应用模式降效参数。
(1)快变改正降效参数解算
快变改正降效参数εfc的计算公式如下：
εfc＝a(t-tu+tlat)2/2 (103)
其中，a为快变改正数降效因子，由增强电文7播发。t为当前时刻。tu为快变改正数的参考时间，如果使用的UDREI是由增强电文2、3、4播发的，则tu为快变改正数的参考时间；如果UDREI是由增强电文6播发并且IODF＝3，则tu也为快变改正数的参考时间；如果IODF≠3，则tu为GEO播发增强电文6第一比特的时间。tlat为系统延迟时间，由增强电文7播发。
(2)距离变化改正降效参数
如果电文7播发的aii＝0，则距离变化改正降效参数εrrc＝0。否则，εrrc将分为两种情况进行计算。
·当前最新的IODF和前一个IODF都不为3
·当前最新的IODF和前一个IODF至少有一个为3
计算过程中使用的参数有：a为快变改正数的降效因子，由增强电文7播发；t为当前时刻；Ifcj为增强电文2、3、4的最短超时间隔，根据电文7播发的a来推算。Brrc由增强电文10播发；IODFcurrent为最新接收的快变改正数信息中的IODF；IODFprevious为在收到最新信息之前收到的快变改正数信息中的IODF；tof为最新接收的快变改正数信息的参考时间；tof，previous为收到最新信息之前收到的快变改正数信息的参考时间；Δt＝tof-tof，previous。
(3)慢变改正降效参数
与慢变改正相关的降效参数根据电文中是否同时包含偏差和偏差变化量，分为两种情况1)增强电文25中速度模式标识为1；2)增强电文25中速度模式标识为0。
·速度模式标识为1
其中，t为当前时刻；t0为长期改正数的参考时间，由增强电文25播发；Iltc_v1为更新间隔，由增强电文10播发；Cltc_lsb由增强电文10播发；Cltc_v1由增强电文10播发。
·速度模式标识为0
其中，t为当前时刻；tltc为GEO播发长期改正数电文第一比特的时间；Iltc_v0为最小更新间隔，由增强电文10播发；Cltc_v0由增强电文10播发；表示向下取整。
(4)应用模式降效参数
在侧向导航/垂直导航(Lateral NAVigation/Vertical NAVigation，LNAV/VNAV)，导航信标性能(Localizer Performance，LP)，垂直导航信标性能(LocalizerPerformance with Vertical guidance，LPV)等飞行阶段，如果快变改正数或慢变改正数超时，但未超出航路至LNAV的时间限制，需要适用于航路至LNAV的降效参数。
适用于航路至LNAV的应用模式降效参数εer的计算公式如下：
其中，Cer由电文10播发。
步骤3.7：电离层改正降效参数解算
电离层改正降效参数的计算公式如下：
其中，RSSiono由增强电文10播发；σGIVE按照步骤3.4中的方法计算；t为当前时刻；tiono为GEO卫星播发第一比特电离层改正数信息的时间；Ciono_step由增强电文10播发；Ciono_ramp由增强电文10播发；Iiono为电离层改正数信息的最小更新间隔，由增强电文10播发；表示向下取整。
步骤四：定位解算；
步骤4.1：对流层延迟估计
对流层延迟估计需要考虑当地温度、水汽压、高度和气压等的影响，使用对流层模型进行修正，对流层延迟估计计算如下：
其中，dhyd与dwet分别表示对流层的干分量和湿分量，为仰角，/>b＝acos[cos(φj-φi)×cos(δj-δi)]，φj和δj分别为卫星j所在位置的纬度和经度，φi和δi分别为监测站i所在位置的纬度和经度；
dhyd与dwet由监测站高度信息及五个气象参数的估值计算：
其中，g＝9.80665m/s2，gm＝9.784m/s2，H为监测站海拔(单位：米)，k1＝77.604K/mbar，k2＝382000K2/mbar，Rd＝287.054J/kg/K；
气象参数气压P(mbar)、温度T(K)、水汽压e(mbar)、温度变化率β(k/m)、水汽变化率λ由监测站的气象传感器提供，也可由当前观测站所在纬度φi和年积日D(当年1月1日起开始计算的天数)插值计算，插值公式如下：
其中，Dmin＝28(φi为北纬)，Dmin＝211(φi为南纬)，ξ0和Δξ分别表示不同纬度的气象参数平均值和季节变化值，计算P、T、e、β、λ时，将ξ分别替换为P、T、e、β、λ，按表3由下式插值得到：
如果φi≤15或φi＝30或φi＝45或φi＝60或φi≥75，直接利用ξ0(φi)和Δξ(φi)在表3中对应的数值通过式(32)计算，其他情况下，以φi＝40为例，对应的φk＝30，φk+1＝45，利用ξ0(φk+1)、ξ0(φk)、Δξ(φk+1)和Δξ(φk)在表3中对应的数值通过式(33)和(34)计算ξ0(φi)和Δξ(φi)。
表3对流层延迟的气象参数表
步骤4.2：监测站位置解算；
卫星j星历位置经卫星星历慢变改正数/>改正后的位置如下：
卫星j时钟偏差经卫星时钟慢变改正数/>改正后的时钟偏差/>如下：
其中，C为光速。
平滑后的伪距观测量经快变改正数/>电离层改正数/>修正，并消除对流层延迟/>和卫星时钟偏差/>得到伪距/>如下：
其中，[xi，yi，zi]为监测站位置，ti为监测接收机时钟偏差。
该伪距方程为非线性方程，要用泰勒级数展开，并取一阶量，将伪距方程转化为线性方程：
其中，
为监测站位置估计值，[Δxi Δyi Δzi]为监测站位置与估计值之间的差值，/>为用户时钟偏差估计值，Δti为监测接收机时钟偏差与估计值之间的差值。
将式(62)进行变换得到：
其中，
式(63)对应的矩阵形式为：
Z＝HX (124)
其中，
X＝[Δxi Δyi Δzi -C·Δti]T，N为监测站观测到的卫星数量。
利用最小二乘法得到：
X＝(HTH)-1HTZ (125)
则监测站位置和时钟偏差为：
/>
取为[xi yi zi]，通过多次迭代，当/>时，可以得到监测站位置和时钟偏差。
步骤4.3：监测站定位误差解算
利用步骤4.2得到的监测站位置[xi yi zi]，并结合测绘标定得到的真实位置[xR，i，yR，i，zR，i]，得到地心地固(Earth-Centered Earth-Fixed，ECEF)坐标系下的监测站定位误差：
[ΔxR，i ΔyR，i ΔzR，i]＝[xR，i，yR，i，zR，i]-[xi yi zi] (128)
从ECEF坐标系到东北天(East North Up，ENU)坐标系的转换矩阵为：
其中，φi和λi分别是监测站所处位置的地理纬度和经度。
得到ENU坐标系下的定位误差如下：
[ΔEi ΔNiΔUi]＝Pi·[ΔxR，i ΔyR，i ΔzR，i]T (130)
基于式(46)得到水平定位误差(Horizontal Position Error，HPE)和垂直定位误差(Vertical Position Error，VPE)如下：
VPE＝ΔUi (132)
步骤五：保护级解算
首先计算监测站与可观测卫星间的观测矩阵G，该矩阵的第j行如下所示：
其中，为监测站与卫星j间的仰角；/>为监测站与卫星j间的方位角。
监测站与第j个可观测的卫星间观测伪距的噪声方差为：
其中，按照步骤3.6中的方法计算；/>按照步骤3.3中的方法计算；
按照步骤4.1中的方法计算。
监测站与卫星间观测伪距的协方差矩阵为W，其对角线元素其余元素为0，通过G和W得到：
针对航路、终端区、非精密进近(Non-Precision Approach，NPA)等飞行阶段，水平保护级(Horizontal Protection Level，HPL)可由下式计算：
HPL＝KH，NPA·dmajor (136)
其中，
针对一类垂直引导进近(APproachwithVertical guidance I，APV-I)、二类垂直引导进近(APV-II)和CAT-I等飞行阶段，HPL和垂直保护级(Vertical Protection Level，VPL)由下式计算：
HPL＝KH，PA·dmajor (137)
VPL＝KVdU (138)
其中，KH，PA＝6.0，KV＝5.33；
步骤六：服务性能评估；
对于针对航路、终端区、NPA等飞行阶段，如果HPE≤HPL，表明SF SBAS服务正常；如果HPE＞HPL，表明SF SBAS服务不能用于导航。对于APV-I、APV-II和CAT-I等飞行阶段，如果HPE≤HPL且VPE≤VPL，表明SF SBAS服务正常；如果HPE＞HPL或者VPE＞VPL，则SF SBAS服务不能够用于引导飞机进行精密进近。当SF SBAS服务不可用时，机场管制员需要将SF SBAS服务不可用的信息告知机场附近准备进近着陆的飞机，飞机需要采用其他导航手段进行进近。

Claims (8)

1.一种用于实时监测单频星基增强系统性能的方法，其特征在于包括下述步骤：
步骤1：数据预处理；
利用监测站观测到卫星的双频伪距观测量和双频载波相位观测量进行数据预处理，预处理的步骤为：
(1)利用前5个采样时刻的载波观测量外推出当前时刻的观测量，并将观测量与当前时刻接收机的载波相位观测量相减并取绝对值，如果该绝对值小于0.055，则认为没有周跳出现，进行后续计算；
(2)利用载波相位观测量对伪距观测量进行平滑，得到平滑后的伪距观测量；
步骤2：导航电文解算；
利用导航电文播发的轨道参数和时钟参数计算卫星星历位置和卫星时钟偏差，利用卫星星历位置和监测站位置计算星历距离；
步骤3：增强电文解算；
步骤3.1：利用SF SBAS增强电文2、3、4中的快变改正信息计算当前时刻的快变改正数；
步骤3.2：利用增强电文25中的慢变改正数信息计算卫星时钟慢变改正数和卫星星历慢变改正数；
步骤3.3：利用增强电文18和26中的电离层格网点信息，推算监测站所在位置电离层改正数；
步骤3.4：基于增强电文2、3、4解算用户差分距离误差，基于增强电文26解算格网电离层垂直误差；
步骤3.5：利用增强电文28中的协方差矩阵信息，解算用户差分距离误差降效降效参数；
步骤3.6：利用增强电文7、电文10、电文25中的相关信息，解算快慢变改正降效参数；
步骤3.7：利用增强电文10和26中的相关信息，解算电离层改正降效参数；
步骤4：定位解算；
(1)基于对流层干湿分量模型，利用监测站所处纬度计算对流层延迟估计；
(2)平滑后的伪距观测量经快变改正数、电离层改正数修正，并消除对流层延迟和卫星时钟偏差后，得到修正后的伪距，对修正后的伪距进行一阶泰勒级数展开，通过迭代法解算监测站位置；
(3)将地心地固坐标系下的定位误差转换到东北天坐标系下，并求解水平和垂直定位误差；
步骤5：保护级解算；利用观测矩阵和观测伪距噪声方差解算水平和垂直保护级；
步骤6：服务性能评估；通过比对定位误差和保护级判断SF SBAS实时服务性能，如果SFSBAS服务出现异常，机场管制员将及时向机场附近的飞机发出告警，停止使用SF SBAS服务。
2.根据权利要求1所述的用于实时监测单频星基增强系统性能的方法，其特征在于：
所述数据预处理的步骤为：
机场附近的监测站采集所监测到全球卫星导航系统卫星的观测数据、GNSS导航电文和SF SBAS增强电文，监测站i观测到卫星j的双频观测数据如下：
其中，和/>分别为L1和L5频点上的伪距观测量；/>和/>分别为L1和L5频点上的载波相位观测量；/>为监测站i和卫星j间的几何距离；/>为对流层延迟；bi为监测站接收机时钟与GNSS系统时之间的偏差；Bj为卫星时钟与GNSS系统时之间的偏差；/>为电离层延迟，对伪距观测量的影响是滞后，对载波相位观测量的影响是超前；/>f1＝1575.42MHz为载波L1的频率，f5＝1176.45MHz为载波L5的频率；/>和/>为伪距观测量上的观测噪声；N1和N5为整周模糊度，由接收机失锁造成；λ1＝C/f1和λ5＝C/f5分别为载波L1和L5的波长，光速C＝29979245m/s；/>和/>为载波相位观测量上的观测噪声，该噪声远远小于伪距观测量上的观察噪声，不同时刻的数据会进行标识，未做说明的数据均为t时刻的数据；
利用监测站i观测到卫星j的双频伪距观测量和双频载波相位观测量进行数据预处理，i＝1，2，…，M，具体步骤如下：
步骤1.1：周跳探测
周跳探测利用前5个采样时刻(t-1,t-2,t-3,t-4,t-5)的载波观测量外推出当前时刻的观测量，并将此观测量与当前时刻接收机的载波相位观测量比较，如果超出门限则认为出现周跳；监测站i对卫星j的周跳监测有下式：
式(5)和式(6)中，a0、α1、α2为拟合系数，[a0，a1，a2]T＝(FTF)-1FTXL1-L5，F取固定值；/>为t时刻观测量组合，/> 和/>分别为t时刻L1和L5频点上的载波相位观测量；/>为由多项式拟合得到的t时刻观测量组合的拟合值；TL1-L5＝0.055为探测门限；
步骤1.2：载波平滑
载波相位观测量通过周跳探测后，则认为没有周跳出现，利用载波相位观测量对伪距观测量进行平滑，首先对载波观测量进行如下变化：
式(9)中
由于前后两个时刻的整周模糊度基本相同，用/>平滑伪距观测量中的噪声；
其中，为t时刻L1频点的伪距观测量，/>为L1频点平滑后的伪距观测量，τ＝100为平滑时间。
3.根据权利要求2所述的用于实时监测单频星基增强系统性能的方法，其特征在于：
所述步骤1.2中，如果L1和L5频点同时出现相同的周跳时，无法检测出周跳，则只利用L5相位观测量进行单频点周跳探测，利用公式(7)、(8)再检测一次：
其中，b0、b1、b2为拟合系数，[b0，b1，b2]T＝(FTF)-1FTXL5；TL5＝0.35为探测门限。
4.根据权利要求1所述的用于实时监测单频星基增强系统性能的方法，其特征在于：
所述步骤2的导航电文解算中，GNSS卫星导航电文中播发的轨道参数为：星历参考时间toe，卫星轨道长半轴αs的平方根，轨道偏心率es，toe时刻的轨道倾角i0，周内时等于0时的轨道升交点赤经Ω0，轨道近地角距ω，toe时刻的平近点角M0，平均运动角速度校正值Δn，轨道倾角变化率i′，轨道升交点赤经变化率升交点角距余弦调和校正振幅Cuc，升交点角距正弦调和校正振幅Cus，轨道半径余弦调和校正振幅Crc，轨道半径正弦调和校正振幅Crs，轨道倾角余弦调和校正振幅Cic，轨道倾角正弦调和校正振幅Cis；利用导航电文播发的轨道参数及卫星位置解算通用算法(此处不再赘述)，得到卫星星历位置/>利用卫星星历位置/>和监测站基准位置[xR，i，yR，i，zR，i]计算星历距离
利用导航电文中播发的参考时间toe、参考时刻的卫星时钟偏差af0、卫星时钟漂移速度af1和卫星时钟漂移速度的变化率af2计算t时刻的卫星时钟偏差
5.根据权利要求1所述的用于实时监测单频星基增强系统性能的方法，其特征在于：
所述步骤3的增强电文解算的详细步骤为：
步骤3.1：快变改正数解算
快变改正数信息由增强电文2、3、4播发，计算公式如下：
PRC(t)＝PRC(tof)+RRC(tof)×(t-tof) (11)
如果aii≠0(由电文7播发)，则
其中，PRCcurrent为最新接收到的快变改正数(由电文2、3、4播发)；PRCprevious为PRCcurrent之前接收到的快变改正数(由电文电文2、3、4播发)；△t＝(tof-tof,previous)；tof为PRCcurrent的参考时刻，为SBAS系统时(SBAS Network Time,SNT)的整秒开始，SBAS GEO卫星播发电文第一比特的时间；tof,previous为PRCprevious的参考时刻；
如果aii＝0，则RRC＝0；
步骤3.2：慢变改正数解算
慢变改正数信息由增强电文25播发，慢变改正数信息包括卫星时钟慢变改正数和卫星星历慢变改正数信息；
卫星时钟慢变改正数通过下式计算：
δΔtSV(t)＝δaf0+δaf1(t-t0) (13)
其中，δaf0为时钟偏差(电文25)；δαf1为时钟偏差变化率，电文25，速度标识为1；如果速度标识为0该值为0；t0为改正数参考时刻，电文25，速度模式标识为1；
卫星星历慢变改正数通过式(14)计算，采用坐标系WGS-84ECEF：
其中，[δxk δyk δzk]为时刻t对应的慢变改正数，[δx δy δz]为时刻t0对应的慢变改正数(由电文25播发)，为时刻t0慢变改正数的变化率，由电文25播发；
当速度标识为0时，上式中的速度分量为0；
步骤3.3：电离层改正数解算；
步骤3.3.1：穿刺点经纬度解算
首先计算穿刺点纬度φpp，公式如(15)和(16)所示：
φpp＝sin-1(sinφucosψpp+cosφusinψppcosA) (15)
式中，ψpp为地球中心角；A为方位角；E为高度角；Re为近似地球半径；hI为最大电子密度层的近似高度；
穿刺点经度λpp的计算公式如(17)和(18)所示：
(1)当φu>70°且tanψpp cosA>tan(π/2-φu)或者φu<-70°且tanψppcos(A+π)>tan(π/2+φu)时，λpp为：
(2)其他情况下，λpp为：
步骤3.3.2：格网点选择
计算完穿刺点经纬度后，监测站根据接收到的电文18确定使用的电离层格网点，电离层格网点的选择与穿刺点所在的纬度有关，选取策略如下：
(1)穿刺点位于南北纬60°之间
如果穿刺点的周围有4个格网点构成5°×5°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则进入下一步；
如果穿刺点的周围有3个格网点构成5°×5°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则进入下一步；
如果穿刺点的周围有4个格网点构成10°×10°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则进入下一步；
如果穿刺点的周围有3个格网点构成10°×10°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则认为电离层改正数不可用；
(2)穿刺点位于北纬60°～75°之间或者南纬60°～75°之间
如果穿刺点的周围有4个格网点构成5°(纬度)×10°(经度)的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则进入下一步；
如果穿刺点的周围有3个格网点构成5°(纬度)×10°(经度)的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则进入下一步；
如果穿刺点的周围有4个格网点构成10°×10°的四边形区域且对应的IGP掩码为1时，则选取这4个格网点，否则进入下一步；
如果穿刺点的周围有3个格网点构成10°×10°的三角形区域且对应的IGP掩码为1时，则选取这3个格网点，否则认为电离层改正数不可用；
(3)穿刺点位于北纬75°～85°之间或者南纬75°～85°之间
如果两个最近的格网点位于75°，两个最近的格网点位于85°且对应的IGP掩码为1时，通过线性内插得到虚拟格网点构建出10°×10°的四边形区域，则选取这四个网格点，否则认为电离层改正数不可用；
(4)穿刺点位于北纬85°
如果四个格网点位于北纬85°，西经180°、西经90°、0°、东经90°，且对应的IGP掩码为1时，则选取这四个格网点，否则认为电离层改正数不可用；
(5)穿刺点位于南纬85°
如果四个格网点位于南纬85°，西经140°、西经50°、东经40°、东经130°，且对应的IGP掩码为1时，则选取这四个格网点，否则认为电离层改正数不可用；
步骤3.3.3：穿刺点电离层垂直延迟解算
由于增强电文播发的是格网点处的电离层垂直延迟，监测站需要利用周围的格网点来计算穿刺点处的电离层延迟；
(1)周边4个格网点
如果在穿刺点周围选取了4个格网点，穿刺点处的电离层垂直延迟τvpp(φpp，λpp)计算公式如(19)所示：
式中，φpp和λpp分别为穿刺点的纬度和经度；τvi为格网点处的电离层垂直延迟；Wi(xpp,ypp)为格网点的权重函数；
W1(xpp,ypp)＝xppypp (20)
W2(xpp,ypp)＝(1-xpp)ypp (21)
W3(xpp,ypp)＝(1-xpp)(1-ypp) (22)
W4(xpp,ypp)＝xpp(1-ypp) (23)
△λpp＝λpp-λ1 (24)
△φpp＝φpp-φ1 (25)
如果穿刺点位于北纬85°和南纬85°之间，xpp和ypp计算公式如(26)和(27)所示：
式中，λ1为穿刺点西边格网点的经度；λ2为穿刺点东边格网点的经度；φ1为穿刺点南边格网点的纬度；φ2为穿刺点北边格网点的纬度；
如果位于北纬85°以北或南纬85°以南，xpp和ypp计算公式如(28)和(29)所示：
式中，λ1为穿刺点东边第二个格网点的经度；λ2为穿刺点西边第二个格网点的经度；λ3为穿刺点西边第一个格网点的经度；λ4为穿刺点东边第一个格网点的经度；
监测站利用格网点的计算穿刺点处的/>公式如(30)所示：
(2)周边3个格网点
如果在穿刺点周围选取了3个格网点且位于北纬75°和南纬75°之间，穿刺点处的电离层垂直延迟τvpp(φpp,λpp)计算公式如(31)所示：
式中，φpp和λpp分别为穿刺点的纬度和经度；τvi为格网点处的电离层垂直延迟；Wi(xpp,ypp)为格网点的权重函数；
W1(xpp,ypp)＝ypp (32)
W2(xpp,ypp)＝1-xpp-ypp (33)
W3(xpp,ypp)＝xpp (34)
监测站利用格网点的计算穿刺点处的/>公式如(35)所示：
步骤3.3.4：穿刺点电离层改正数解算
穿刺点处的电离层改正数计算公式如(36)所示：
式中，τspp(φpp,λpp)为穿刺点处电离层倾斜延迟，Fpp为倾斜因子；
穿刺点处的为：
步骤3.4：完好性参数解算
完好性参数包括用户差分距离误差(User Differential Range Error,UDRE)和格网电离层垂直误差(Grid Ionospheric Vertical Error,GIVE)；
UDRE以户差分距离误差索引(User Differential Range Error Indicator,UDREI)的形式通过增强电文2、3、4播发，UDREI与UDRE之间的转换关系固定；
GIVE以格网电离层垂直误差索引(Grid Ionospheric Vertical Error Indicator,GIVEI)的形式通过增强电文26播发，GIVEI与GIVE之间的转换关系固定；
步骤3.5：UDRE降效参数解算
协方差矩阵信息由增强电文28播发，利用电文中的信息计算上三角矩阵R：
R＝SF·E (39)
其中，SF＝2(scale exponent-5)，
利用矩阵R计算协方差矩阵C：
C＝RT·R (40)
投影到监测站上的UDRE降效参数：
其中，I为卫星到监测站的4维方向矢量，前三维是单位方向矢量，第四维是1；εc＝Ccovariance·SF，Ccovariance由电文10播发，如果电文10中的Ccovariance无效，则该值为0；
步骤3.6：快慢变改正降效参数解算
通过快变改正数和慢变改正数修正后的残差方差通过快慢变降效参数表征，计算公式如下：
其中，RSSUDRE由增强电文10播发；σUDRE按照步骤3.4中的方法计算；δUDRE按照步骤3.5中的方法计算；εfc为快变改正降效参数；εrrc为距离变化改正降效参数；εltc为慢变改正降效参数；εer为应用模式降效参数；
(1)快变改正降效参数解算
快变改正降效参数εfc的计算公式如下：
εfc＝a(t-tu+tlat)2/2 (43)
其中，a为快变改正数降效因子，由增强电文7播发，t为当前时刻，tu为快变改正数的参考时间，如果使用的UDREI是由增强电文2、3、4播发的，则tu为快变改正数的参考时间；如果UDREI是由增强电文6播发并且IODF＝3，则tu也为快变改正数的参考时间；如果IODF≠3，则tu为GEO播发增强电文6第一比特的时间；tlat为系统延迟时间，由增强电文7播发；
(2)距离变化改正降效参数
如果电文7播发的aii＝0，则距离变化改正降效参数εrrc＝0，否则，εrrc将分为两种情况进行计算：
a)当前最新的IODF和前一个IODF都不为3，则：
b)当前最新的IODF和前一个IODF至少有一个为3，则：
计算过程中使用的参数有：a为快变改正数的降效因子，由增强电文7播发；t为当前时刻；Ifc,j为增强电文2、3、4的最短超时间隔，根据电文7播发的a来推算，Brrc由增强电文10播发；IODFcurrent为最新接收的快变改正数信息中的IODF；IODFprevious为在收到最新信息之前收到的快变改正数信息中的IODF；tof为最新接收的快变改正数信息的参考时间；tof,previous为收到最新信息之前收到的快变改正数信息的参考时间；Δt＝tof-tof,previous；
(3)慢变改正降效参数
与慢变改正相关的降效参数根据电文中是否同时包含偏差和偏差变化量，分为两种情况：
1)增强电文25中速度模式标识为1；
2)增强电文25中速度模式标识为0；
当速度模式标识为1，则：
其中，t为当前时刻；t0为长期改正数的参考时间，由增强电文25播发；Iltc_v1为更新间隔，由增强电文10播发；Cltc_lsb由增强电文10播发；Cltc_v1由增强电文10播发；
当速度模式标识为0，则：
其中，t为当前时刻；tltc为GEO播发长期改正数电文第一比特的时间；Iltc_v0为最小更新间隔，由增强电文10播发；Cltc_v0由增强电文10播发；表示向下取整；
(4)应用模式降效参数
在侧向导航/垂直导航(Lateral NAVigation/Vertical NAVigation,LNAV/VNAV)，导航信标性能(Localizer Performance,LP)，垂直导航信标性能(Localizer Performancewith Vertical guidance,LPV)等飞行阶段，如果快变改正数或慢变改正数超时，但未超出航路至LNAV的时间限制，需要适用于航路至LNAV的降效参数；
适用于航路至LNAV的应用模式降效参数εer的计算公式如下：
其中，Cer由电文10播发；
步骤3.7：电离层改正降效参数解算
电离层改正降效参数的计算公式如下：
其中，RSSiono由增强电文10播发；σGIVE按照步骤3.4中的方法计算；t为当前时刻；tiono为GEO卫星播发第一比特电离层改正数信息的时间；Ciono_step由增强电文10播发；Ciono_ramp由增强电文10播发；Iiono为电离层改正数信息的最小更新间隔，由增强电文10播发；表示向下取整。
6.根据权利要求1所述的用于实时监测单频星基增强系统性能的方法，其特征在于：所述步骤4的定位解算步骤为：
步骤4.1：对流层延迟估计
对流层模型进行修正，对流层延迟估计计算如下：
其中，dhyd与dwet分别表示对流层的干分量和湿分量， 为仰角，/>b＝acos[cos(φj-φi)×cos(δj-δi)]，φj和δj分别为卫星j所在位置的纬度和经度，φi和δi分别为监测站i所在位置的纬度和经度；
dhyd与dwet由监测站高度信息及五个气象参数的估值计算：
其中，g＝9.80665m/s2，gm＝9.784m/s2，H为监测站海拔，k1＝77.604K/mbar，k2＝382000K2/mbar，Rd＝287.054J/kg/K；
气象参数气压P(mbar)、温度T(K)、水汽压e(mbar)、温度变化率β(k/m)、水汽变化率λ由监测站的气象传感器提供，也可由当前观测站所在纬度φi和年积日D(当年1月1日起开始计算的天数)插值计算，插值公式如下：
其中，Dmin＝28(φi为北纬)，Dmin＝211(φi为南纬)，ξ0和Δξ分别表示不同纬度的气象参数平均值和季节变化值，计算P、T、e、β、λ时，将ξ分别替换为P、T、e、β、λ，按表3由下式插值得到：
如果φi≤15或φi＝30或φi＝45或φi＝60或φi≥75，直接利用ξ0(φi)和Δξ(φi)通过式(32)计算，其他情况下，以φi＝40为例，对应的φk＝30，φk+1＝45，利用ξ0(φk+1)、ξ0(φk)、Δξ(φk+1)和Δξ(φk)通过式(33)和(34)计算ξ0(φi)和Δξ(φi)；
步骤4.2：监测站位置解算；
卫星j星历位置经卫星星历慢变改正数/>改正后的位置/>如下：
卫星j时钟偏差经卫星时钟慢变改正数/>改正后的时钟偏差/>如下：
其中，C为光速；
平滑后的伪距观测量经快变改正数/>电离层改正数/>修正，并消除对流层延迟/>和卫星时钟偏差/>得到伪距/>如下：
其中，[xi，yi，zi]为监测站位置，ti为监测接收机时钟偏差；
该伪距方程为非线性方程，要用泰勒级数展开，并取一阶量，将伪距方程转化为线性方程：
/>
其中， 为监测站位置估计值，[Δxi Δyi Δzi]为监测站位置与估计值之间的差值，/>为用户时钟偏差估计值，Δti为监测接收机时钟偏差与估计值之间的差值；
将式(62)进行变换得到：
其中，
式(63)对应的矩阵形式为：
Z＝HX (51)
其中，X＝[Δxi Δyi Δzi-C·Δti]T，N为监测站观测到的卫星数量；
利用最小二乘法得到：
X＝(HTH)-1HTZ (52)
则监测站位置和时钟偏差为：
取为[xi yi zi]，通过多次迭代，当/>时，可以得到监测站位置和时钟偏差；
步骤4.3：监测站定位误差解算
利用步骤4.2得到的监测站位置[xi yi zi]，并结合测绘标定得到的真实位置[xR，i，yR，i，zR，i]，得到地心地固(Earth-Centered Earth-Fixed,ECEF)坐标系下的监测站定位误差：
[ΔxR，i ΔyR，i ΔzR，i]＝[xR，i，yR，i，zR，i]-[xi yi zi] (55)
从ECEF坐标系到东北天(East North Up,ENU)坐标系的转换矩阵为：
其中，φi和λi分别是监测站所处位置的地理纬度和经度；
得到ENU坐标系下的定位误差如下：
[ΔEi ΔNi ΔUi]＝Pi·[ΔxR，i ΔyR，i ΔzR，i]T (57)
基于式(46)得到如下：
VPE＝ΔUi (72)
得到水平定位误差和垂直定位误差。
7.根据权利要求1所述的用于实时监测单频星基增强系统性能的方法，其特征在于：所述步骤5的保护级解算的步骤为：
首先计算监测站与可观测卫星间的观测矩阵G，该矩阵的第j行如下所示：
其中，为监测站与卫星j间的仰角；/>为监测站与卫星j间的方位角；
监测站与第j个可观测的卫星间观测伪距的噪声方差为：
其中，按照步骤3.6中的方法计算；/>按照步骤3.3中的方法计算；按照步骤4.1中的方法计算；
监测站与卫星间观测伪距的协方差矩阵为W，其对角线元素其余元素为0，通过G和W得到：
水平保护级HPL可由下式计算：
HPL＝KH，NPA·dmajor (60)
其中，KH，NPA＝6.18；
针对一类垂直引导进近(APproach with Vertical guidance I,APV-I)、二类垂直引导进近(APV-II)和CAT-I等飞行阶段，HPL和垂直保护级(Vertical Protection Level,VPL)由下式计算：
HPL＝KH，PA·dmajor (77)
VPL＝KVdU (78)
其中，KH，PA＝6.0，KV＝5.33。
8.根据权利要求1所述的用于实时监测单频星基增强系统性能的方法，其特征在于：步骤6：服务性能评估；
对于针对航路、终端区、NPA等飞行阶段，如果HPE≤HPL，表明SF SBAS服务正常；如果HPE＞HPL，表明SF SBAS服务不能用于导航；对于APV-I、APV-II和CAT-I等飞行阶段，如果HPE≤HPL且VPE≤VPL，表明SF SBAS服务正常；如果HPE＞HPL或者VPE＞VPL，则SF SBAS服务不能够用于引导飞机进行精密进近；当SF SBAS服务不可用时，机场管制员需要将SF SBAS服务不可用的信息告知机场附近准备进近着陆的飞机，飞机需要采用其他导航手段进行进近。
CN202110520577.5A 2021-05-13 2021-05-13 一种用于实时监测单频星基增强系统性能的方法 Active CN113253303B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110520577.5A CN113253303B (zh) 2021-05-13 2021-05-13 一种用于实时监测单频星基增强系统性能的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110520577.5A CN113253303B (zh) 2021-05-13 2021-05-13 一种用于实时监测单频星基增强系统性能的方法

Publications (2)

Publication Number Publication Date
CN113253303A CN113253303A (zh) 2021-08-13
CN113253303B true CN113253303B (zh) 2023-11-10

Family

ID=77181537

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110520577.5A Active CN113253303B (zh) 2021-05-13 2021-05-13 一种用于实时监测单频星基增强系统性能的方法

Country Status (1)

Country Link
CN (1) CN113253303B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114152962A (zh) * 2021-11-14 2022-03-08 中国电子科技集团公司第二十研究所 一种星基增强系统服务范围确定的方法
CN117335899B (zh) * 2023-10-09 2024-04-19 中国人民解放军32021部队 一种北斗星基增强服务降效程度评估方法
CN117579391B (zh) * 2024-01-16 2024-04-05 中国人民解放军战略支援部队航天工程大学 一种基于sbas认证服务的钟差电文降频优化方法和系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2642316A1 (de) * 2012-03-21 2013-09-25 Astrium GmbH Bestimmung von Grid Ionospheric Vertical Delay Paramentern für ein Satelliten-gestützes Erweiterungssystem
CN110988929A (zh) * 2019-12-21 2020-04-10 中国电子科技集团公司第二十研究所 电离层影响下的gbas系统性能评估方法及装置
CN111596315A (zh) * 2020-05-23 2020-08-28 中国电子科技集团公司第二十研究所 一种用于实时监测双频多星座星基增强系统性能的方法
CN112034489A (zh) * 2020-07-20 2020-12-04 中国科学院空天信息创新研究院 一种基于多源数据融合的全球电离层格网生成方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9945954B2 (en) * 2014-11-20 2018-04-17 Honeywell International Inc. Using space-based augmentation system (SBAS) grid ionosphere vertical error (GIVE) information to mitigate ionosphere errors for ground based augmentation systems (GBAS)

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2642316A1 (de) * 2012-03-21 2013-09-25 Astrium GmbH Bestimmung von Grid Ionospheric Vertical Delay Paramentern für ein Satelliten-gestützes Erweiterungssystem
CN110988929A (zh) * 2019-12-21 2020-04-10 中国电子科技集团公司第二十研究所 电离层影响下的gbas系统性能评估方法及装置
CN111596315A (zh) * 2020-05-23 2020-08-28 中国电子科技集团公司第二十研究所 一种用于实时监测双频多星座星基增强系统性能的方法
CN112034489A (zh) * 2020-07-20 2020-12-04 中国科学院空天信息创新研究院 一种基于多源数据融合的全球电离层格网生成方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
北斗卫星导航系统完好性参数研究;崔瑞云;倪育德;王琳琳;刘逸;;现代导航(01);第17-22页 *

Also Published As

Publication number Publication date
CN113253303A (zh) 2021-08-13

Similar Documents

Publication Publication Date Title
US10690775B2 (en) Crowdsourcing atmospheric correction data
CN113253303B (zh) 一种用于实时监测单频星基增强系统性能的方法
CN111596315B (zh) 一种用于实时监测双频多星座星基增强系统性能的方法
US10078140B2 (en) Navigation satellite system positioning involving the generation of advanced correction information
AU2009333478B2 (en) Methods and systems to increase accuracy in the navigation of single frequency receivers
Yang et al. Chinese navigation satellite systems
US11378699B2 (en) System and method for determining GNSS positioning corrections
CN110007326B (zh) 一种用于星基增强系统的双频测距误差参数生成方法
US11740364B1 (en) System and method to reduce PPP filter convergence time using LEO frequency band signals
Jokinen Enhanced ambiguity resolution and integrity monitoring methods for precise point positioning
Abdelazeem et al. MGR-DCB: a precise model for multi-constellation GNSS receiver differential code bias
CN101419274B (zh) 电离层延迟误差的获取方法及系统
CN105044733B (zh) 一种高精度的导航卫星tgd参数标定方法
Boon et al. Precise aircraft positioning by fast ambiguity resolution using improved troposphere modeling
CN113267793B (zh) 一种基于外部增强信息的gbas对流层参数生成方法
Xing et al. Analysis of RDSS positioning accuracy based on RNSS wide area differential technique
Kee Wide Area Differential GPS (WADGPS)
Thompson et al. Utilization and Validation of DSS-17 on the CAPSTONE Lunar Mission
US11714198B2 (en) Single-epoch pseudo-range positioning under varying ionosphere delays
Aleshkin et al. Calculation of navigation corrections for a single-frequency GNSS receiver based on satellite radio occultation data
Baba et al. Estimation Of global positioning system measurement errors For GAGAN applications
Zheng Generation of network-based differential corrections for regional GNSS services
Hu et al. Research on the performance of multi-GNSS medium length baseline RTK with LEO-augmented
Kim et al. Preliminary Test Results of RTK-aided Conical Domain Model for SBAS Ionospheric Correction
CN105572712B (zh) 一种北斗系统载波相位整周模糊度实时确定方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant