CN111596315A - 一种用于实时监测双频多星座星基增强系统性能的方法 - Google Patents

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

Info

Publication number
CN111596315A
CN111596315A CN202010445180.XA CN202010445180A CN111596315A CN 111596315 A CN111596315 A CN 111596315A CN 202010445180 A CN202010445180 A CN 202010445180A CN 111596315 A CN111596315 A CN 111596315A
Authority
CN
China
Prior art keywords
satellite
monitoring station
sbas
dfmc
time
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.)
Granted
Application number
CN202010445180.XA
Other languages
English (en)
Other versions
CN111596315B (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 CN202010445180.XA priority Critical patent/CN111596315B/zh
Publication of CN111596315A publication Critical patent/CN111596315A/zh
Application granted granted Critical
Publication of CN111596315B publication Critical patent/CN111596315B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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/07—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing data for correcting measured positioning data, e.g. DGPS [differential GPS] or ionosphere corrections
    • G01S19/072—Ionosphere corrections
    • 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
    • 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/10—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing dedicated supplementary positioning signals
    • G01S19/11—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing dedicated supplementary positioning signals wherein the cooperating elements are pseudolites or satellite radio beacon positioning system signal repeaters
    • G01S19/115—Airborne or satellite based pseudolites or repeaters
    • 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/13—Receivers
    • G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements

Landscapes

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

Abstract

本发明提供了一种用于实时监测双频多星座星基增强系统性能的方法，实时监测DFMC SBAS服务性能，确保机场附近区域内DFMC SBAS服务的可用性，利用机场附近监测站的观测数据、GNSS导航电文以及DFMC SBAS增强电文，实时解算定位误差和保护级，通过对比定位误差和保护级来实时监测DFMC SBAS在机场附近的服务性能。本发明具有较强的工程实用性，利用机场附近的监测站作为参考基准，实时解算定位误差和保护级，评估DFMC SBAS在机场的实时服务性能，解决了飞机在进近和着陆过程中因无法实时解算定位误差而导致的完好性风险，保证了DFMC SBAS服务的完好性性能。

Description

一种用于实时监测双频多星座星基增强系统性能的方法
技术领域
本发明涉及卫星导航增强技术领域，是双频多星座(Dual-Frequency Multi-Constellation,DFMC)星基增强系统(Satellite Based Augmentation System,SBAS)中一种实时监测服务性能的方法。
背景技术
目前正在运行的SBAS系统均为单频(Single-Frequency,SF)SBAS。由于电离层异常对服务性能的影响，SF SBAS服务性能尚未达到一类精密进近(CATegory-I,CAT-I)的指标要求。为了消除电离层异常对服务性能的影响，并利用多卫星导航系统的几何布局提高增强星座服务性能，国际民航组织(International Civil Aviation Organization,ICAO)双频多星座星基增强系统标准与建议工作组(DFMC SBAS SARPS working group,DS2)正在研究并制定DFMC SBAS国际标准，以期实现CAT-I指标服务性能。DFMC SBAS播发的增强电文类型如下表所示：
表1 DFMC SBAS增强电文类型
Figure BDA0002505586170000011
中国、美国、欧盟、日本等国家已经根据ICAO发布的DFMC SBAS国际标准草案启动了DFMC SBAS的设计、验证和建设工作，预计2023年前后将向航空用户提供DFMC SBAS服务，并逐渐取代SF SBAS服务成为航空运输的主要导航手段。突发的DFMC SBAS服务中断将可能给航空安全带来严重的后果。为了满足航空用户对基于DFMC SBAS飞行进近的高安全性要求，需要对机场附近区域的DFMC SBAS服务性能进行实时监测，并将实时监测结果发送给机场管制员。如果DFMC SBAS服务出现异常，机场管制员将及时向机场附近的飞机发出告警，停止使用DFMC SBAS服务，改用其他导航手段。
目前，国内外尚未有公开文献对DFMC SBAS服务性能实时监测方法的描述。
发明内容
为了克服现有技术的不足，本发明提供一种用于实时监测双频多星座星基增强系统性能的方法，本发明实时监测DFMC SBAS服务性能，确保机场附近区域内DFMC SBAS服务的可用性，利用机场附近监测站的观测数据、GNSS导航电文以及DFMC SBAS增强电文，实时解算定位误差和保护级，通过对比定位误差和保护级来实时监测DFMC SBAS在机场附近的服务性能。
本发明解决其技术问题所采用的技术方案的步骤为：
步骤一：数据预处理；
机场附近的监测站采集所监测到全球卫星导航系统(Global NavigationSatellite System,GNSS)卫星的观测数据、GNSS导航电文和DFMC SBAS增强电文，监测站i观测到卫星j的双频观测数据如下：
Figure BDA0002505586170000021
Figure BDA0002505586170000022
Figure BDA0002505586170000023
Figure BDA0002505586170000024
其中，
Figure BDA0002505586170000026
和
Figure BDA0002505586170000027
分别为L1和L5频点上的伪距观测量；
Figure BDA0002505586170000028
和
Figure BDA0002505586170000029
分别为L1和L5频点上的载波相位观测量；
Figure BDA00025055861700000210
为监测站i和卫星j间的几何距离；
Figure BDA00025055861700000211
为对流层延迟；bi为监测站接收机时钟与GNSS系统时之间的偏差；Bj为卫星时钟与GNSS系统时之间的偏差；
Figure BDA00025055861700000212
为电离层延迟，对伪距观测量的影响是滞后，对载波相位观测量的影响是超前；
Figure BDA0002505586170000025
f1＝1575.42MHz为载波L1的频率，f5＝1176.45MHz为载波L5的频率；
Figure BDA00025055861700000213
和
Figure BDA00025055861700000214
为伪距观测量上的观测噪声；N1和N5为整周模糊度，λ1＝C/f1和λ5＝C/f5分别为载波L1和L5的波长，C为光速；
Figure BDA0002505586170000037
和
Figure BDA0002505586170000038
为载波相位观测量上的观测噪声；
利用监测站i观测到卫星j的双频伪距观测量和双频载波相位观测量进行数据预处理，i＝1，2，…，M，具体步骤如下：
步骤1.1：周跳探测；
周跳探测利用前5个采样时刻(t-1,t-2,t-3,t-4,t-5)的载波观测量外推出当前时刻的观测量，并将此观测量与当前时刻接收机的载波相位观测量比较，如果超出门限，则认为出现周跳；监测站i对卫星j的周跳监测有下式：
Figure BDA0002505586170000031
Figure BDA0002505586170000032
式(5)和式(6)中，a0、a1、a2为拟合系数，[a0，a1，a2]T＝(FTF)-1FTXL1-L5，
Figure BDA0002505586170000033
F取固定值；
Figure BDA0002505586170000039
为t时刻观测量组合，
Figure BDA0002505586170000034
Figure BDA00025055861700000310
和
Figure BDA00025055861700000311
分别为t时刻L1和L5频点上的载波相位观测量；
Figure BDA00025055861700000312
为由多项式拟合得到的t时刻观测量组合的拟合值；TL1-L5＝0.055为探测门限；
如果L1和L5频点同时出现相同的周跳时，无法检测出周跳，则只利用L5相位观测量进行单频点周跳探测，利用公式(7)、(8)再检测一次：
Figure BDA0002505586170000035
Figure BDA0002505586170000036
其中，b0、b1、b2为拟合系数，[b0，b1，b2]T＝(FTF)-1FTXL5；TL5＝0.35为探测门限；
步骤1.2：双频载波平滑
载波相位观测量通过周跳探测后，则认为没有周跳出现，利用载波相位观测量对伪距观测量进行平滑，首先对载波观测量进行如下变化：
Figure BDA0002505586170000041
Figure BDA0002505586170000042
式(9)、(10)中
Figure BDA0002505586170000043
用
Figure BDA0002505586170000047
平滑伪距观测量中的噪声；
Figure BDA0002505586170000044
其中，Lk表示L1或L5频点，
Figure BDA0002505586170000048
为相应频点的伪距观测量，
Figure BDA0002505586170000049
为相应频点平滑后的伪距观测量，τ＝100为平滑时间；
步骤1.3：消除电离层延迟；
利用L1和L5频点平滑后的伪距观测量消除电离层延迟，消除电离层延迟后的伪距观测量
Figure BDA00025055861700000410
为：
Figure BDA0002505586170000045
步骤二：导航电文解算；
GNSS卫星导航电文中播发的轨道参数为：星历参考时间toe，卫星轨道长半轴αs的平方根，轨道偏心率es，toe时刻的轨道倾角i0，周内时等于0时的轨道升交点赤经Ω0，轨道近地角距ω，toe时刻的平近点角M0，平均运动角速度校正值Δn，轨道倾角变化率i′，轨道升交点赤经变化率
Figure BDA00025055861700000414
升交点角距余弦调和校正振幅Cuc，升交点角距正弦调和校正振幅Cus，轨道半径余弦调和校正振幅Crc，轨道半径正弦调和校正振幅Crs，轨道倾角余弦调和校正振幅Cic，轨道倾角正弦调和校正振幅Cis；利用导航电文播发的轨道参数及卫星位置解算算法，得到卫星星历位置
Figure BDA00025055861700000415
利用卫星星历位置
Figure BDA00025055861700000413
和监测站基准位置[xR，i，yR，i，zR，i]计算星历距离
Figure BDA0002505586170000046
利用导航电文中播发的参考时间toe、参考时刻的卫星时钟偏差αf0、卫星时钟漂移速度αf1和卫星时钟漂移速度的变化率αf2计算t时刻的卫星时钟偏差
Figure BDA0002505586170000051
步骤三：增强电文解算
步骤3.1：差分改正数解算；
星历改正数解算从DFMC SBAS电文32读取星历位置改正信息，参考坐标系为地心地固坐标系，计算公式如下：
Figure BDA0002505586170000052
其中，[δxj δyj δzj]T为卫星j在时刻t的星历改正数；
Figure BDA0002505586170000055
为参考时刻tD卫星j的星历改正数，由电文32播发；
Figure BDA0002505586170000056
为参考时刻tD卫星j的星历改正数变化率，由电文32播发；
钟差改正数解算从DFMC SBAS电文32读取钟差改正信息，计算公式如下：
Figure BDA0002505586170000053
其中，
Figure BDA0002505586170000057
为卫星j在时刻t的钟差改正数，单位：秒；
Figure BDA0002505586170000058
为参考时刻tD卫星j的钟差改正数，由电文32播发；
Figure BDA0002505586170000059
为参考时刻tD卫星j的钟差改正数变化率，由电文32播发；C为光速；
步骤3.2：完好性参数解算；
首先利用电文32播发中播发的协方差矩阵信息比例因子sc和矩阵元素Ex，y，其中，x,y为1,2,3,4，计算卫星j的
Figure BDA00025055861700000510
计算公式如下：
Rj＝SFj·Ej (15)
Figure BDA0002505586170000054
Figure BDA0002505586170000061
其中，
Figure BDA0002505586170000062
Ij为卫星j到监测站的4维方向矢量，前三维是单位方向矢量，第四维是1；
Figure BDA0002505586170000063
Ccovariance由电文37播发；
完好性参数解算利用DFMC SBAS电文32和37播发的信息进行计算，
Figure BDA00025055861700000610
计算公式如下：
Figure BDA0002505586170000064
Figure BDA0002505586170000065
其中，
Figure BDA00025055861700000611
利用电文34、35、36播发的DFREI以及电文37播发的DFREI映射表得出；
Figure BDA00025055861700000612
为电文32第一比特播发时间；RCORR和ICORR由电文37播发，
Figure BDA00025055861700000613
由电文32播发，当t-tCORR≤ICORR时，
Figure BDA0002505586170000066
当t-tCORR＞ICORR时，(RCORR)sv＝RCORR；
Figure BDA00025055861700000614
为向下取整；
步骤四：定位解算；
步骤4.1：对流层延迟估计
使用对流层模型进行修正，对流层延迟估计
Figure BDA00025055861700000615
计算如下：
Figure BDA0002505586170000067
其中，dhyd与dwet分别表示对流层的干分量和湿分量，
Figure BDA0002505586170000068
Figure BDA00025055861700000616
为仰角，
Figure BDA0002505586170000069
b＝acos[cos(φj-φi)×cos(δj-δi)]，φj和δj分别为卫星j所在位置的纬度和经度，φi和δi分别为监测站i所在位置的纬度和经度；
dhyd与dwet由监测站高度信息及五个气象参数的估值计算：
Figure BDA0002505586170000071
Figure BDA0002505586170000072
Figure BDA0002505586170000073
Figure BDA0002505586170000074
其中，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)、水汽变化率λ由监测站的气象传感器提供；
计算P、T、e、β、λ时，将ξ分别替换为P、T、e、β、λ，按表2由下式插值得到：
Figure BDA0002505586170000075
Figure BDA0002505586170000076
表2对流层延迟的气象参数表
Figure BDA0002505586170000081
如果φi≤15或φi＝30或φi＝45或φi＝60或φi≥75，直接利用ξ0(φi)和Δξ(φi)在表1中对应的数值通过式(25)计算，其他情况下，以φi＝40为例，对应的φk＝30，φk+1＝45，利用ξ0(φk+1)、ξ0(φk)、Δξ(φk+1)和Δξ(φk)在表2中对应的数值通过式(26)和(27)计算ξ0(φi)和Δξ(φi)；
步骤4.2：监测站位置解算；
卫星j星历位置
Figure BDA0002505586170000084
经星历改正数[δxj δyj δzj]改正后的位置
Figure BDA0002505586170000085
如下：
Figure BDA0002505586170000082
卫星j时钟偏差
Figure BDA0002505586170000086
经钟差改正数
Figure BDA0002505586170000087
改正后的时钟偏差
Figure BDA0002505586170000088
如下：
Figure BDA0002505586170000083
其中，C为光速；
从平滑后的伪距观测量
Figure BDA00025055861700000910
中消除对流层延迟
Figure BDA00025055861700000911
和卫星时钟偏差
Figure BDA00025055861700000916
得到伪距
Figure BDA00025055861700000913
如下：
Figure BDA0002505586170000091
其中，[xi，yi，zi]为监测站位置，ti为监测接收机时钟偏差；
该伪距方程为非线性方程，要用泰勒级数展开，并取一阶量，将伪距方程转化为线性方程；
Figure BDA0002505586170000092
其中，
Figure BDA0002505586170000093
Figure BDA0002505586170000094
Figure BDA00025055861700000914
为监测站位置估计值，[xi yi zi]为监测站位置与估计值之间的差值，
Figure BDA00025055861700000915
为用户时钟偏差估计值，Δti为监测接收机时钟偏差与估计值之间的差值。
将式(31)进行变换得到：
Figure BDA0002505586170000095
其中，
Figure BDA0002505586170000096
式(32)对应的矩阵形式为：
Z＝HX (33)
其中，
Figure BDA0002505586170000097
X＝[Δxi Δyi Δzi -C·Δti]T，N为监测站观测到的卫星数量；
利用最小二乘法得到：
X＝(HTH)-1HTZ (34)
则监测站位置和时钟偏差为：
Figure BDA0002505586170000098
Figure BDA0002505586170000099
取
Figure BDA0002505586170000105
为[xi yi zi]，通过多次迭代当
Figure BDA0002505586170000106
时，得到监测站位置和时钟偏差；
步骤4.3：监测站定位误差解算
利用定位解算得到的监测站位置[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] (37)
从ECEF坐标系到东北天(East North Up,ENU)坐标系的转换矩阵为：
Figure BDA0002505586170000101
其中，φi和λi分别是监测站所处位置的地理纬度和经度；
得到ENU坐标系下的定位误差如下：
[ΔEi ΔNi ΔUi]＝Pi·[ΔxR，i ΔyR，i ΔzR，i]T (39)
基于式(39)得到水平定位误差(Horizontal Position Error,HPE)和垂直定位误差(Vertical Position Error,VPE)如下：
Figure BDA0002505586170000102
VPE＝ΔUi (41)
步骤五：保护级解算
首先计算监测站与可观测卫星间的观测矩阵G，该矩阵的第j行如下所示：
Figure BDA0002505586170000103
其中，
Figure BDA0002505586170000107
为监测站与卫星j间的仰角；
Figure BDA0002505586170000104
为监测站与卫星j间的方位角；
监测站与第j个可观测的卫星间观测伪距的噪声方差为：
Figure BDA0002505586170000111
其中，
Figure BDA0002505586170000112
Figure BDA0002505586170000116
为卫星j对流层投影系数；
监测站与卫星间观测伪距的协方差矩阵为W，其对角线元素
Figure BDA0002505586170000113
其余元素为0，通过G和W得到：
Figure BDA0002505586170000114
针对航路、终端区、非精密进近(Non-Precision Approach,NPA)飞行阶段，水平保护级(Horizontal Protection Level,HPL)由式(45)计算：
HPL＝KH，NPA·dmajor (45)
其中，
Figure BDA0002505586170000115
KH，NPA＝6.18；
针对一类垂直引导进近(APproach with Vertical guidance I,APV-I)、二类垂直引导进近(APV-II)和CAT-I飞行阶段，HPL和垂直保护级(Vertical Protection Level,VPL)由下式计算：
HPL＝KH，PA·dmajor (46)
VPL＝KVdU (47)
其中，KH，PA＝6.0，KV＝5.33；
步骤六：服务性能评估；
对于针对航路、终端区、NPA飞行阶段，如果HPE≤HPL，表明DFMC SBAS服务正常；如果HPE＞HPL，表明DFMC SBAS服务不能用于导航；
对于APV-I、APV-II和CAT-I飞行阶段，如果HPE≤HPL且VPE≤VPL，表明DFMC SBAS服务正常；如果HPE＞HPL或者VPE＞VPL，则DFMC SBAS服务不能够用于引导飞机进行精密进近，当DFMC SBAS服务不可用时，机场管制员将DFMC SBAS服务不可用的信息告知机场附近准备进近着陆的飞机，飞机需要采用其他导航手段进行进近。
所述步骤4.1中，气象参数气压P(mbar)、温度T(K)、水汽压e(mbar)、温度变化率β(k/m)、水汽变化率λ由当前观测站所在纬度φi和年积日D插值计算，年积日D为当年1月1日起开始计算的天数，插值公式如下：
Figure BDA0002505586170000121
其中，φi为北纬时，Dmin＝28，φi为南纬时，Dmin＝211，ξ0和Δξ分别表示不同纬度的气象参数平均值和季节变化值。
本发明的有益效果在于：
1)提出了实时监测机场附近DFMC SBAS服务性能的方法，给出了明确的处理流程和实施步骤，具有较强的工程实用性，为中国民航SBAS监测与服务系统的建设提供了理论依据和实施思路；
2)利用机场附近的监测站作为参考基准，实时解算定位误差和保护级，评估DFMCSBAS在机场的实时服务性能，解决了飞机在进近和着陆过程中因无法实时解算定位误差而导致的完好性风险，保证了DFMC SBAS服务的完好性性能。
附图说明
图1是本发明DFMC SBAS服务性能实时监测步骤流程图。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
本发明是一种用于实时监测双频多星座星基增强系统性能的方法，具体步骤如图1所示：
步骤一：数据预处理；
机场附近的监测站采集所监测到全球卫星导航系统(Global NavigationSatellite System,GNSS)卫星的观测数据、GNSS导航电文和DFMC SBAS增强电文，监测站i观测到卫星j的双频观测数据如下：
Figure BDA0002505586170000122
Figure BDA0002505586170000123
Figure BDA0002505586170000131
Figure BDA0002505586170000132
其中，
Figure BDA0002505586170000136
和
Figure BDA0002505586170000137
分别为L1和L5频点上的伪距观测量；
Figure BDA0002505586170000138
和
Figure BDA0002505586170000139
分别为L1和L5频点上的载波相位观测量；
Figure BDA00025055861700001310
为监测站i和卫星j间的几何距离；
Figure BDA00025055861700001311
为对流层延迟；bi为监测站接收机时钟与GNSS系统时之间的偏差；Bj为卫星时钟与GNSS系统时之间的偏差；
Figure BDA00025055861700001312
为电离层延迟，对伪距观测量的影响是滞后，对载波相位观测量的影响是超前；
Figure BDA0002505586170000133
f1＝1575.42MHz为载波L1的频率，f5＝1176.45MHz为载波L5的频率；
Figure BDA00025055861700001313
和
Figure BDA00025055861700001314
为伪距观测量上的观测噪声；N1和N5为整周模糊度，由接收机失锁造成；λ1＝c/f1和λ5＝C/f5分别为载波L1和L5的波长，光速C＝299792458m/s；
Figure BDA00025055861700001315
和
Figure BDA00025055861700001316
为载波相位观测量上的观测噪声，该噪声远远小于伪距观测量上的观察噪声。不同时刻的数据会进行标识，未做说明的数据均为t时刻的数据。
利用监测站i(i＝1，2，…，M)观测到卫星j的双频伪距观测量和双频载波相位观测量进行数据预处理，具体步骤如下：
步骤1.1：周跳探测；
周跳探测利用前5个采样时刻(t-1,t-2,t-3,t-4,t-5)的载波观测量外推出当前时刻的观测量，并将此观测量与当前时刻接收机的载波相位观测量比较，如果超出门限，则认为出现周跳；监测站i对卫星j的周跳监测有下式：
Figure BDA0002505586170000134
Figure BDA0002505586170000135
式(5)和式(6)中，α0、α1、α2为拟合系数，[α0，α1，α2]T＝(FTF)-1FTXL1-L5，
Figure BDA0002505586170000141
F取固定值；
Figure BDA0002505586170000142
为t时刻观测量组合，
Figure BDA0002505586170000143
Figure BDA00025055861700001410
和
Figure BDA00025055861700001411
分别为t时刻L1和L5频点上的载波相位观测量；
Figure BDA00025055861700001412
为由多项式拟合得到的t时刻观测量组合的拟合值；TL1-L5＝0.055为探测门限；
如果L1和L5频点同时出现相同的周跳时，无法检测出周跳，则只利用L5相位观测量进行单频点周跳探测，利用公式(7)、(8)再检测一次：
Figure BDA0002505586170000144
Figure BDA0002505586170000145
其中，b0、b1、b2为拟合系数，[b0，b1，b2]T＝(FTF)-1FTXL5；TL5＝0.35为探测门限；
步骤1.2：双频载波平滑
载波相位观测量通过周跳探测后，则认为没有周跳出现，利用载波相位观测量对伪距观测量进行平滑，首先对载波观测量进行如下变化：
Figure BDA0002505586170000146
Figure BDA0002505586170000147
式(9)、(10)中
Figure BDA0002505586170000148
由于
Figure BDA00025055861700001413
前后两个时刻的整周模糊度基本相同，用
Figure BDA00025055861700001414
平滑伪距观测量中的噪声；
Figure BDA0002505586170000149
其中，Lk表示L1或L5频点，
Figure BDA0002505586170000156
为相应频点的伪距观测量，
Figure BDA0002505586170000157
为相应频点平滑后的伪距观测量，τ＝100为平滑时间；
步骤1.3：消除电离层延迟；
利用L1和L5频点平滑后的伪距观测量消除电离层延迟，消除电离层延迟后的伪距观测量
Figure BDA0002505586170000158
为：
Figure BDA0002505586170000151
步骤二：导航电文解算；
GNSS卫星导航电文中播发的轨道参数为：星历参考时间toe，卫星轨道长半轴as的平方根，轨道偏心率es，toe时刻的轨道倾角i0，周内时等于0时的轨道升交点赤经Ω0，轨道近地角距ω，toe时刻的平近点角M0，平均运动角速度校正值Δn，轨道倾角变化率i′，轨道升交点赤经变化率
Figure BDA0002505586170000159
升交点角距余弦调和校正振幅Cuc，升交点角距正弦调和校正振幅Cus，轨道半径余弦调和校正振幅Crc，轨道半径正弦调和校正振幅Crs，轨道倾角余弦调和校正振幅Cic，轨道倾角正弦调和校正振幅Cis；利用导航电文播发的轨道参数及卫星位置解算算法，得到卫星星历位置
Figure BDA0002505586170000152
利用卫星星历位置
Figure BDA0002505586170000153
和监测站基准位置[xR，i，yR，i，zR，i]计算星历距离
Figure BDA0002505586170000154
利用导航电文中播发的参考时间toe、参考时刻的卫星时钟偏差αf0、卫星时钟漂移速度αf1和卫星时钟漂移速度的变化率αf2计算t时刻的卫星时钟偏差
Figure BDA0002505586170000155
步骤三：增强电文解算
步骤3.1：差分改正数解算；
星历改正数解算从DFMC SBAS电文32读取星历位置改正信息，参考坐标系为地心地固坐标系，计算公式如下：
Figure BDA0002505586170000161
其中，[δxj δyj δzj]T为卫星j在时刻t的星历改正数；
Figure BDA0002505586170000166
为参考时刻tD卫星j的星历改正数，由电文32播发；
Figure BDA0002505586170000167
为参考时刻tD卫星j的星历改正数变化率，由电文32播发；
钟差改正数解算从DFMC SBAS电文32读取钟差改正信息，计算公式如下：
Figure BDA0002505586170000162
其中，
Figure BDA0002505586170000168
为卫星j在时刻t的钟差改正数，单位：秒；
Figure BDA0002505586170000169
为参考时刻tD卫星j的钟差改正数，由电文32播发；
Figure BDA00025055861700001610
为参考时刻tD卫星j的钟差改正数变化率，由电文32播发；C为光速；
步骤3.2：完好性参数解算；
首先利用电文32播发中播发的协方差矩阵信息比例因子sc和矩阵元素Ex，y，其中，x,y为1,2,3,4，计算卫星j的
Figure BDA0002505586170000163
计算公式如下：
Rj＝SFj·Ej (15)
Figure BDA0002505586170000164
Figure BDA0002505586170000165
其中，
Figure BDA0002505586170000171
Ij为卫星j到监测站的4维方向矢量，前三维是单位方向矢量，第四维是1；
Figure BDA0002505586170000172
Ccovariance由电文37播发；
完好性参数解算利用DFMC SBAS电文32和37播发的信息进行计算，
Figure BDA0002505586170000177
计算公式如下：
Figure BDA0002505586170000173
Figure BDA0002505586170000174
其中，
Figure BDA0002505586170000178
利用电文34、35、36播发的DFREI以及电文37播发的DFREI映射表得出；
Figure BDA0002505586170000179
为电文32第一比特播发时间；RCORR和lCORR由电文37播发，
Figure BDA00025055861700001710
由电文32播发，当t-tCORR≤ICORR时，
Figure BDA0002505586170000175
当t-tCORR＞ICORR时，(RCORR)sv＝RCORR；
Figure BDA00025055861700001711
为向下取整；
步骤四：定位解算；
步骤4.1：对流层延迟估计
对流层延迟估计需要考虑当地温度、水汽压、高度和气压等的影响，使用对流层模型进行修正，对流层延迟估计
Figure BDA00025055861700001712
计算如下：
Figure BDA0002505586170000176
其中，dhyd与dwet分别表示对流层的干分量和湿分量，
Figure BDA0002505586170000181
Figure BDA00025055861700001810
为仰角，
Figure BDA0002505586170000182
b＝acos[cos(φj-φi)×cos(δj-δi)]，φj和δj分别为卫星j所在位置的纬度和经度，φi和δi分别为监测站i所在位置的纬度和经度；
dhyd与dwet由监测站高度信息及五个气象参数的估值计算：
Figure BDA0002505586170000183
Figure BDA0002505586170000184
Figure BDA0002505586170000185
Figure BDA0002505586170000186
其中，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插值计算，年积日D为当年1月1日起开始计算的天数，插值公式如下：
Figure BDA0002505586170000187
其中，φi为北纬时，Dmin＝28，φi为南纬时，Dmin＝211，ξ0和Δξ分别表示不同纬度的气象参数平均值和季节变化值，计算P、T、e、β、λ时，将ξ分别替换为P、T、e、β、λ，按表2由下式插值得到：
Figure BDA0002505586170000188
Figure BDA0002505586170000189
表2对流层延迟的气象参数表
Figure BDA0002505586170000191
如果φi≤15或φi＝30或φi＝45或φi＝60或φi≥75，直接利用ξ0(φi)和Δξ(φi)在表1中对应的数值通过式(25)计算，其他情况下，以φi＝40为例，对应的φk＝30，φk+1＝45，利用ξ0(φk+1)、ξ0(φk)、Δξ(φk+1)和Δξ(φk)在表2中对应的数值通过式(26)和(27)计算ξ0(φi)和Δξ(φi)；
步骤4.2：监测站位置解算；
卫星j星历位置
Figure BDA0002505586170000194
经星历改正数[δxj δyj δzj]改正后的位置
Figure BDA0002505586170000195
如下：
Figure BDA0002505586170000192
卫星j时钟偏差
Figure BDA0002505586170000196
经钟差改正数
Figure BDA0002505586170000197
改正后的时钟偏差
Figure BDA0002505586170000198
如下：
Figure BDA0002505586170000193
其中，C为光速；
从平滑后的伪距观测量
Figure BDA00025055861700002010
中消除对流层延迟
Figure BDA00025055861700002011
和卫星时钟偏差
Figure BDA00025055861700002016
得到伪距
Figure BDA00025055861700002013
如下：
Figure BDA0002505586170000201
其中，[xi，yi，zi]为监测站位置，ti为监测接收机时钟偏差；
该伪距方程为非线性方程，要用泰勒级数展开，并取一阶量，将伪距方程转化为线性方程；
Figure BDA0002505586170000202
其中，
Figure BDA0002505586170000203
Figure BDA0002505586170000204
Figure BDA00025055861700002014
为监测站位置估计值，[Δxi Δyi Δzi]为监测站位置与估计值之间的差值，
Figure BDA00025055861700002015
为用户时钟偏差估计值，Δti为监测接收机时钟偏差与估计值之间的差值。
将式(31)进行变换得到：
Figure BDA0002505586170000205
其中，
Figure BDA0002505586170000206
式(32)对应的矩阵形式为：
Z＝HX (33)
其中，
Figure BDA0002505586170000207
X＝[Δxi Δyi Δzi -C·Δti]T，N为监测站观测到的卫星数量；
利用最小二乘法得到：
X＝(HTH)-1HTZ (34)
则监测站位置和时钟偏差为：
Figure BDA0002505586170000208
Figure BDA0002505586170000209
取
Figure BDA0002505586170000216
为[xi yi zi]，通过多次迭代当
Figure BDA0002505586170000217
时，得到监测站位置和时钟偏差；
步骤4.3：监测站定位误差解算
利用定位解算得到的监测站位置[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] (37)
从ECEF坐标系到东北天(East North Up,ENU)坐标系的转换矩阵为：
Figure BDA0002505586170000211
其中，φi和λi分别是监测站所处位置的地理纬度和经度；
得到ENU坐标系下的定位误差如下：
[ΔEi ΔNi ΔUi]＝Pi·[ΔxR，i ΔyR，i ΔzR，i]T (39)
基于式(39)得到水平定位误差(Horizontal Position Error,HPE)和垂直定位误差(Vertical Position Error,VPE)如下：
Figure BDA0002505586170000212
VPE＝ΔUi (41)
步骤五：保护级解算
首先计算监测站与可观测卫星间的观测矩阵G，该矩阵的第j行如下所示：
Figure BDA0002505586170000213
其中，
Figure BDA0002505586170000218
为监测站与卫星j间的仰角；
Figure BDA0002505586170000214
为监测站与卫星j间的方位角；
监测站与第j个可观测的卫星间观测伪距的噪声方差为：
Figure BDA0002505586170000215
其中，
Figure BDA0002505586170000221
Figure BDA0002505586170000225
为卫星j对流层投影系数；
监测站与卫星间观测伪距的协方差矩阵为W，其对角线元素
Figure BDA0002505586170000222
其余元素为0，通过G和W得到：
Figure BDA0002505586170000223
针对航路、终端区、非精密进近(Non-Precision Approach,NPA)等飞行阶段，水平保护级(Horizontal Protection Level,HPL)由式(45)计算：
HPL＝KH，NPA·dmajor (45)
其中，
Figure BDA0002505586170000224
KH，NPA＝6.18；
针对一类垂直引导进近(APproach with Vertical guidance I,APV-I)、二类垂直引导进近(APV-II)和CAT-I等飞行阶段，HPL和垂直保护级(Vertical ProtectionLevel,VPL)由下式计算：
HPL＝KH，PA·dmajor (46)
VPL＝KVdU (47)
其中，KH，PA＝6.0，KV＝5.33；
步骤六：服务性能评估；
对于针对航路、终端区、NPA飞行阶段，如果HPE≤HPL，表明DFMC SBAS服务正常；如果HPE＞HPL，表明DFMC SBAS服务不能用于导航；
对于APV-I、APV-II和CAT-I飞行阶段，如果HPE≤HPL且VPE≤VPL，表明DFMC SBAS服务正常；如果HPE＞HPL或者VPE＞VPL，则DFMC SBAS服务不能够用于引导飞机进行精密进近，当DFMC SBAS服务不可用时，机场管制员将DFMC SBAS服务不可用的信息告知机场附近准备进近着陆的飞机，飞机需要采用其他导航手段进行进近。

Claims (2)

1.一种用于实时监测双频多星座星基增强系统性能的方法，其特征在于包括下述步骤：
步骤一：数据预处理；
机场附近的监测站采集所监测到全球卫星导航系统(Global Navigation SatelliteSystem，GNSS)卫星的观测数据、GNSS导航电文和DFMC SBAS增强电文，监测站i观测到卫星j的双频观测数据如下：
Figure FDA0002505586160000011
Figure FDA0002505586160000012
Figure FDA0002505586160000013
Figure FDA0002505586160000014
其中，
Figure FDA0002505586160000015
和
Figure FDA0002505586160000016
分别为L1和L5频点上的伪距观测量；
Figure FDA0002505586160000017
和
Figure FDA0002505586160000018
分别为L1和L5频点上的载波相位观测量；
Figure FDA0002505586160000019
为监测站i和卫星j间的几何距离；
Figure FDA00025055861600000110
为对流层延迟；bi为监测站接收机时钟与GNSS系统时之间的偏差；Bj为卫星时钟与GNSS系统时之间的偏差；
Figure FDA00025055861600000111
为电离层延迟，对伪距观测量的影响是滞后，对载波相位观测量的影响是超前；
Figure FDA00025055861600000112
f1＝1575.42MHz为载波L1的频率，f5＝1176.45MHz为载波L5的频率；
Figure FDA00025055861600000113
和
Figure FDA00025055861600000114
为伪距观测量上的观测噪声；N1和N5为整周模糊度，λ1＝C/f1和λ5＝C/f5分别为载波L1和L5的波长，C为光速；
Figure FDA00025055861600000115
和
Figure FDA00025055861600000116
为载波相位观测量上的观测噪声；
利用监测站i观测到卫星j的双频伪距观测量和双频载波相位观测量进行数据预处理，i＝1，2，…，M，具体步骤如下：
步骤1.1：周跳探测；
周跳探测利用前5个采样时刻(t-1，t-2，t-3，t-4，t-5)的载波观测量外推出当前时刻的观测量，并将此观测量与当前时刻接收机的载波相位观测量比较，如果超出门限，则认为出现周跳；监测站i对卫星j的周跳监测有下式：
Figure FDA0002505586160000021
Figure FDA0002505586160000022
式(5)和式(6)中，a0、a1、a2为拟合系数，[a0，a1，a2]T＝(FTF)-1FTXL1-L5，
Figure FDA0002505586160000023
F取固定值；
Figure FDA0002505586160000024
为t时刻观测量组合，
Figure FDA0002505586160000025
Figure FDA0002505586160000026
和分别为t时刻L1和L5频点上的载波相位观测量；
Figure FDA0002505586160000028
为由多项式拟合得到的t时刻观测量组合的拟合值；TL1-L5＝0.055为探测门限；
如果L1和L5频点同时出现相同的周跳时，无法检测出周跳，则只利用L5相位观测量进行单频点周跳探测，利用公式(7)、(8)再检测一次：
Figure FDA0002505586160000029
Figure FDA00025055861600000210
其中，b0、b1、b2为拟合系数，[b0，b1，b2]T＝(FTF)-1FTXL5；TL5＝0.35为探测门限；
步骤1.2：双频载波平滑
载波相位观测量通过周跳探测后，则认为没有周跳出现，利用载波相位观测量对伪距观测量进行平滑，首先对载波观测量进行如下变化：
Figure FDA00025055861600000211
Figure FDA00025055861600000212
式(9)、(10)中
Figure FDA00025055861600000213
用
Figure FDA00025055861600000214
平滑伪距观测量中的噪声；
Figure FDA0002505586160000031
其中，Lk表示L1或L5频点，
Figure FDA0002505586160000032
为相应频点的伪距观测量，
Figure FDA0002505586160000033
为相应频点平滑后的伪距观测量，τ＝100为平滑时间；
步骤1.3：消除电离层延迟；
利用L1和L5频点平滑后的伪距观测量消除电离层延迟，消除电离层延迟后的伪距观测量
Figure FDA0002505586160000034
为：
Figure FDA0002505586160000035
步骤二：导航电文解算；
GNSS卫星导航电文中播发的轨道参数为：星历参考时间toe，卫星轨道长半轴as的平方根，轨道偏心率es，toe时刻的轨道倾角i0，周内时等于0时的轨道升交点赤经Ω0，轨道近地角距ω，toe时刻的平近点角M0，平均运动角速度校正值Δn，轨道倾角变化率i′，轨道升交点赤经变化率
Figure FDA00025055861600000310
升交点角距余弦调和校正振幅Cuc，升交点角距正弦调和校正振幅Cus，轨道半径余弦调和校正振幅Crc，轨道半径正弦调和校正振幅Crs，轨道倾角余弦调和校正振幅Cic，轨道倾角正弦调和校正振幅Cis；利用导航电文播发的轨道参数及卫星位置解算算法，得到卫星星历位置
Figure FDA0002505586160000036
利用卫星星历位置
Figure FDA0002505586160000037
和监测站基准位置[xR，i，yR，i，zR，i]计算星历距离
Figure FDA0002505586160000038
利用导航电文中播发的参考时间toe、参考时刻的卫星时钟偏差αf0、卫星时钟漂移速度af1和卫星时钟漂移速度的变化率af2计算t时刻的卫星时钟偏差
Figure FDA0002505586160000039
步骤三：增强电文解算
步骤3.1：差分改正数解算；
星历改正数解算从DFMC SBAS电文32读取星历位置改正信息，参考坐标系为地心地固坐标系，计算公式如下：
Figure FDA0002505586160000041
其中，[δxj δyj δzj]T为卫星j在时刻t的星历改正数；
Figure FDA0002505586160000042
为参考时刻tD卫星j的星历改正数，由电文32播发；
Figure FDA0002505586160000043
为参考时刻tD卫星j的星历改正数变化率，由电文32播发；
钟差改正数解算从DFMC SBAS电文32读取钟差改正信息，计算公式如下：
Figure FDA0002505586160000044
其中，
Figure FDA0002505586160000045
为卫星j在时刻t的钟差改正数，单位：秒；
Figure FDA0002505586160000046
为参考时刻tD卫星j的钟差改正数，由电文32播发；
Figure FDA0002505586160000047
为参考时刻tD卫星j的钟差改正数变化率，由电文32播发；C为光速；
步骤3.2：完好性参数解算；
首先利用电文32播发中播发的协方差矩阵信息比例因子sc和矩阵元素Ex，y，其中，x，y为1，2，3，4，计算卫星j的
Figure FDA0002505586160000048
计算公式如下：
Rj＝SFj·Ej (15)
Figure FDA00025055861600000412
Figure FDA0002505586160000049
其中，
Figure FDA00025055861600000410
Ij为卫星j到监测站的4维方向矢量，前三维是单位方向矢量，第四维是1；
Figure FDA00025055861600000411
Ccovariance由电文37播发；
完好性参数解算利用DFMC SBAS电文32和37播发的信息进行计算，
Figure FDA0002505586160000051
计算公式如下：
Figure FDA0002505586160000052
Figure FDA0002505586160000053
其中，
Figure FDA0002505586160000054
利用电文34、35、36播发的DFREI以及电文37播发的DFREI映射表得出；
Figure FDA0002505586160000055
为电文32第一比特播发时间；RCORR和ICORR由电文37播发，
Figure FDA0002505586160000056
由电文32播发，当t-tCORR≤ICORR时，
Figure FDA0002505586160000057
当t-tCORR＞ICORR时，(RCORR)sv＝RCORR；
Figure FDA0002505586160000058
为向下取整；
步骤四：定位解算；
步骤4.1：对流层延迟估计
使用对流层模型进行修正，对流层延迟估计
Figure FDA0002505586160000059
计算如下：
Figure FDA00025055861600000510
其中，dhyd与dwet分别表示对流层的干分量和湿分量，
Figure FDA00025055861600000511
Figure FDA00025055861600000512
为仰角，
Figure FDA00025055861600000513
b＝acos[cos(φj-φi)×cos(δj-δi)]，φj和δj分别为卫星j所在位置的纬度和经度，φi和δi分别为监测站i所在位置的纬度和经度；
dhyd与dwet由监测站高度信息及五个气象参数的估值计算：
Figure FDA00025055861600000514
Figure FDA00025055861600000515
Figure FDA0002505586160000061
Figure FDA0002505586160000062
其中，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)、水汽变化率λ由监测站的气象传感器提供；
计算P、T、e、β、λ时，将ξ分别替换为P、T、e、β、λ，按表2由下式插值得到：
Figure FDA0002505586160000063
Figure FDA0002505586160000064
表2对流层延迟的气象参数表
Figure FDA0002505586160000065
如果φi≤15或φi＝30或φi＝45或φi＝60或φi≥75，直接利用ξ0(φi)和Δξ(φi)在表1中对应的数值通过式(25)计算，其他情况下，以φi＝40为例，对应的φk＝30，φk+1＝45，利用ξ0(φk+1)、ξ0(φk)、Δξ(φk+1)和Δξ(φk)在表2中对应的数值通过式(26)和(27)计算ξ0(φi)和Δξ(φi)；
步骤4.2：监测站位置解算；
卫星j星历位置
Figure FDA0002505586160000071
经星历改正数[δxj δyj δzj]改正后的位置
Figure FDA0002505586160000072
如下：
Figure FDA0002505586160000073
卫星j时钟偏差
Figure FDA0002505586160000074
经钟差改正数
Figure FDA0002505586160000075
改正后的时钟偏差
Figure FDA0002505586160000076
如下：
Figure FDA0002505586160000077
其中，C为光速；
从平滑后的伪距观测量
Figure FDA0002505586160000078
中消除对流层延迟
Figure FDA0002505586160000079
和卫星时钟偏差
Figure FDA00025055861600000710
得到伪距
Figure FDA00025055861600000711
如下：
Figure FDA00025055861600000712
其中，[xi，yi，zi]为监测站位置，ti为监测接收机时钟偏差；
该伪距方程为非线性方程，要用泰勒级数展开，并取一阶量，将伪距方程转化为线性方程；
Figure FDA00025055861600000713
其中，
Figure FDA00025055861600000714
Figure FDA00025055861600000715
Figure FDA00025055861600000716
为监测站位置估计值，[Δxi Δyi Δzi]为监测站位置与估计值之间的差值，
Figure FDA00025055861600000717
为用户时钟偏差估计值，Δti为监测接收机时钟偏差与估计值之间的差值；
将式(31)进行变换得到：
Figure FDA00025055861600000718
其中，
Figure FDA00025055861600000719
式(32)对应的矩阵形式为：
Z＝HX (33)
其中，
Figure FDA0002505586160000081
X＝[Δxi Δyi Δzi -C·Δti]T，N为监测站观测到的卫星数量；
利用最小二乘法得到：
X＝(HTH)-1HTZ (34)
则监测站位置和时钟偏差为：
Figure FDA0002505586160000082
Figure FDA0002505586160000083
取
Figure FDA0002505586160000084
为[xi yi zi]，通过多次迭代当
Figure FDA0002505586160000085
时，得到监测站位置和时钟偏差；
步骤4.3：监测站定位误差解算
利用定位解算得到的监测站位置[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] (37)
从ECEF坐标系到东北天(East North Up，ENU)坐标系的转换矩阵为：
Figure FDA0002505586160000086
其中，φi和λi分别是监测站所处位置的地理纬度和经度；
得到ENU坐标系下的定位误差如下：
[ΔEi ΔNi ΔUi]＝Pi·[ΔxR，i ΔyR，i ΔzR，i]T (39)
基于式(39)得到水平定位误差(Horizontal Position Error，HPE)和垂直定位误差(Vertical Position Error，VPE)如下：
Figure FDA0002505586160000091
VPE＝ΔUi (41)
步骤五：保护级解算
首先计算监测站与可观测卫星间的观测矩阵G，该矩阵的第j行如下所示：
Figure FDA0002505586160000092
其中，
Figure FDA0002505586160000093
为监测站与卫星j间的仰角；
Figure FDA0002505586160000094
为监测站与卫星j间的方位角；
监测站与第j个可观测的卫星间观测伪距的噪声方差为：
Figure FDA0002505586160000095
其中，
Figure FDA0002505586160000096
Figure FDA0002505586160000097
为卫星j对流层投影系数；
监测站与卫星间观测伪距的协方差矩阵为W，其对角线元素
Figure FDA0002505586160000098
其余元素为0，通过G和W得到：
Figure FDA0002505586160000099
针对航路、终端区、非精密进近(Non-Precision Approach，NPA)飞行阶段，水平保护级(Horizontal Protection Level，HPL)由式(45)计算：
HPL＝KH，NPA·dmajor (45)
其中，
Figure FDA00025055861600000910
KHNPA＝6.18：
针对一类垂直引导进近(APproach with Vertical guidance I，APV-I)、二类垂直引导进近(APV-II)和CAT-I飞行阶段，HPL和垂直保护级(Vertical Protection Level，VPL)由下式计算：
HPL＝KH，pA·dmajor (46)
VPL＝KVdU (47)
其中，KH，PA＝6.0，KV＝5.33；
步骤六：服务性能评估；
对于针对航路、终端区、NPA飞行阶段，如果HPE≤HPL，表明DFMC SBAS服务正常；如果HPE＞HPL，表明DFMC SBAS服务不能用于导航；
对于APV-I、APV-II和CAT-I飞行阶段，如果HPE≤HPL且VPE≤VPL，表明DFMC SBAS服务正常；如果HPE＞HPL或者VPE＞VPL，则DFMC SBAS服务不能够用于引导飞机进行精密进近，当DFMC SBAS服务不可用时，机场管制员将DFMC SBAS服务不可用的信息告知机场附近准备进近着陆的飞机，飞机需要采用其他导航手段进行进近。
2.根据权利要求1所述的一种用于实时监测双频多星座星基增强系统性能的方法，，其特征在于：
所述步骤4.1中，气象参数气压P(mbar)、温度T(K)、水汽压e(mbar)、温度变化率β(k/m)、水汽变化率λ由当前观测站所在纬度φi和年积日D插值计算，年积日D为当年1月1日起开始计算的天数，插值公式如下：
Figure FDA0002505586160000101
其中，φi为北纬时，Dmin＝28，φi为南纬时，Dmin＝211，ξ0和Δξ分别表示不同纬度的气象参数平均值和季节变化值。
CN202010445180.XA 2020-05-23 2020-05-23 一种用于实时监测双频多星座星基增强系统性能的方法 Active CN111596315B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010445180.XA CN111596315B (zh) 2020-05-23 2020-05-23 一种用于实时监测双频多星座星基增强系统性能的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010445180.XA CN111596315B (zh) 2020-05-23 2020-05-23 一种用于实时监测双频多星座星基增强系统性能的方法

Publications (2)

Publication Number Publication Date
CN111596315A true CN111596315A (zh) 2020-08-28
CN111596315B CN111596315B (zh) 2022-07-22

Family

ID=72189279

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010445180.XA Active CN111596315B (zh) 2020-05-23 2020-05-23 一种用于实时监测双频多星座星基增强系统性能的方法

Country Status (1)

Country Link
CN (1) CN111596315B (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112462394A (zh) * 2020-11-10 2021-03-09 中国科学院国家授时中心 一种基于gnss全视比对的远程时间频率配送方法
CN113253303A (zh) * 2021-05-13 2021-08-13 中国电子科技集团公司第二十研究所 一种用于实时监测单频星基增强系统性能的方法
CN113568020A (zh) * 2021-09-27 2021-10-29 长沙学院 一种顾及硬件频间差的卫星导航定位误差修正方法和装置
CN114047526A (zh) * 2022-01-12 2022-02-15 天津七一二通信广播股份有限公司 基于双频双星座gbas的电离层异常监测方法及装置
CN114244545A (zh) * 2020-09-08 2022-03-25 中国科学院空天信息创新研究院 Sbas电文认证方法、装置及接收机
CN116718195A (zh) * 2023-08-03 2023-09-08 中国科学院空天信息创新研究院 基于双频定位的飞行导航方法、装置、设备和存储介质

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110184595A1 (en) * 2010-01-27 2011-07-28 Airbus Operations (S.A.S.) Method And Device For Aiding The Piloting Of An Aircraft During A Final Approach Phase
US20130079958A1 (en) * 2011-09-22 2013-03-28 Ecole Nationale De L'aviation Civile (E.N.A.C.) Method and system for determining the position of an aircraft during its approach to a landing runway
US20140035778A1 (en) * 2012-08-03 2014-02-06 Thales Method of monitoring the integrity of radio-navigation stations in a satellite based augmentation system
CN105068088A (zh) * 2015-06-29 2015-11-18 北京航空航天大学 双频卫星导航星基增强系统可用性预测方法
CN106468774A (zh) * 2016-09-09 2017-03-01 北京航空航天大学 一种应用于星基增强系统的星历星钟改正参数及空间信号完好性参数方法
CN109001767A (zh) * 2018-08-11 2018-12-14 西北工业大学 一种利用低轨卫星增强多基准一致性检测的方法
CN109542084A (zh) * 2018-11-19 2019-03-29 北京航空航天大学 一种星基增强系统完好性故障仿真方法
CN110007326A (zh) * 2019-04-15 2019-07-12 中国电子科技集团公司第二十研究所 一种用于星基增强系统的双频测距误差参数生成方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110184595A1 (en) * 2010-01-27 2011-07-28 Airbus Operations (S.A.S.) Method And Device For Aiding The Piloting Of An Aircraft During A Final Approach Phase
US20130079958A1 (en) * 2011-09-22 2013-03-28 Ecole Nationale De L'aviation Civile (E.N.A.C.) Method and system for determining the position of an aircraft during its approach to a landing runway
US20140035778A1 (en) * 2012-08-03 2014-02-06 Thales Method of monitoring the integrity of radio-navigation stations in a satellite based augmentation system
CN105068088A (zh) * 2015-06-29 2015-11-18 北京航空航天大学 双频卫星导航星基增强系统可用性预测方法
CN106468774A (zh) * 2016-09-09 2017-03-01 北京航空航天大学 一种应用于星基增强系统的星历星钟改正参数及空间信号完好性参数方法
CN109001767A (zh) * 2018-08-11 2018-12-14 西北工业大学 一种利用低轨卫星增强多基准一致性检测的方法
CN109542084A (zh) * 2018-11-19 2019-03-29 北京航空航天大学 一种星基增强系统完好性故障仿真方法
CN110007326A (zh) * 2019-04-15 2019-07-12 中国电子科技集团公司第二十研究所 一种用于星基增强系统的双频测距误差参数生成方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
DIDIER FLAMENT等: "RAIM in dual frequency/multi constellation APV/LPV operations in aeronautics", 《2010 5TH ESA WORKSHOP ON SATELLITE NAVIGATION TECHNOLOGIES AND EUROPEAN WORKSHOP ON GNSS SIGNALS AND SIGNAL PROCESSING (NAVITEC)》 *
JUNJIE BAO等: "Ionospheric Anomaly Detection to Support the BD SBAS", 《IEEE ACCESS 》 *
伍维甲等: "基于多接收机的本地监测系统监测算法研究", 《科学技术与工程》 *
杜娟: "星基增强系统互操作及其关键技术研究", 《中国优秀博硕士学位论文全文数据库(博士)基础科学辑(月刊)》 *
牛飞: "GNSS完好性增强理论与方法研究", 《中国优秀博硕士学位论文全文数据库(博士)基础科学辑(月刊)》 *
耿永超等: "地基区域完好性监测系统发展探讨", 《2008年导航定位技术进步与创新专家论坛论文集》 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114244545A (zh) * 2020-09-08 2022-03-25 中国科学院空天信息创新研究院 Sbas电文认证方法、装置及接收机
CN114244545B (zh) * 2020-09-08 2023-11-14 中国科学院空天信息创新研究院 Sbas电文认证方法、装置及接收机
CN112462394A (zh) * 2020-11-10 2021-03-09 中国科学院国家授时中心 一种基于gnss全视比对的远程时间频率配送方法
CN112462394B (zh) * 2020-11-10 2023-08-25 中国科学院国家授时中心 一种基于gnss全视比对的远程时间频率配送方法
CN113253303A (zh) * 2021-05-13 2021-08-13 中国电子科技集团公司第二十研究所 一种用于实时监测单频星基增强系统性能的方法
CN113253303B (zh) * 2021-05-13 2023-11-10 中国电子科技集团公司第二十研究所 一种用于实时监测单频星基增强系统性能的方法
CN113568020A (zh) * 2021-09-27 2021-10-29 长沙学院 一种顾及硬件频间差的卫星导航定位误差修正方法和装置
CN114047526A (zh) * 2022-01-12 2022-02-15 天津七一二通信广播股份有限公司 基于双频双星座gbas的电离层异常监测方法及装置
CN116718195A (zh) * 2023-08-03 2023-09-08 中国科学院空天信息创新研究院 基于双频定位的飞行导航方法、装置、设备和存储介质
CN116718195B (zh) * 2023-08-03 2023-11-14 中国科学院空天信息创新研究院 基于双频定位的飞行导航方法、装置、设备和存储介质

Also Published As

Publication number Publication date
CN111596315B (zh) 2022-07-22

Similar Documents

Publication Publication Date Title
CN111596315B (zh) 一种用于实时监测双频多星座星基增强系统性能的方法
CN113253303B (zh) 一种用于实时监测单频星基增强系统性能的方法
CN110007326B (zh) 一种用于星基增强系统的双频测距误差参数生成方法
CA2250196C (en) Gps/irs global position determination method and apparatus with integrity loss provisions
US11378699B2 (en) System and method for determining GNSS positioning corrections
Yang et al. Chinese navigation satellite systems
CN105044747B (zh) 一种基于多星共视和滤波的时间同步装置及其方法
CN101395443B (zh) 混合定位方法和设备
CN104316943B (zh) 一种伪距离和多普勒组合差分定位系统及方法
CN106324629A (zh) 一种bds_gps_glonass融合精密单点定位方法
CN104656108A (zh) 一种顾及高程差异的稀疏参考站网络天顶对流层延迟建模方法
CN111913203B (zh) 一种动态基线定位域监测方法
McGraw Tropospheric error modeling for high integrity airborne GNSS navigation
CN105510946A (zh) 一种bds卫星载波相位整周模糊度快速解算方法
CN116577810A (zh) 一种卫星导航高精度服务完好性监测方法及装置
CN113267793B (zh) 一种基于外部增强信息的gbas对流层参数生成方法
CN103389502B (zh) 基于多个地面基站高精度确定载体加速度的方法
Basile et al. Multi-frequency precise point positioning using GPS and Galileo data with smoothed ionospheric corrections
Khojasteh et al. Introduction to global navigation satellite systems and its errors
Kablak Procedure for determining tropospheric delays in the ZAKPOS/UA-EUPOS network of active reference stations
Jgouta et al. Usage of a correction model to enhance the evaluation of the zenith tropospheric delay
Aleshkin et al. Calculation of navigation corrections for a single-frequency GNSS receiver based on satellite radio occultation data
El-Mowafy Using Multiple Reference Station GPS Networks for Aircraft Precision Approach and Airport Surface Navigation
Parveen et al. Position Error Calculations for IRNSS System Using Pseudo Range Method
CN114430808A (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