CN112099063A - 用于北斗星基增强用户误差最大投影方向快速搜索的方法 - Google Patents

用于北斗星基增强用户误差最大投影方向快速搜索的方法 Download PDF

Info

Publication number
CN112099063A
CN112099063A CN202010861987.1A CN202010861987A CN112099063A CN 112099063 A CN112099063 A CN 112099063A CN 202010861987 A CN202010861987 A CN 202010861987A CN 112099063 A CN112099063 A CN 112099063A
Authority
CN
China
Prior art keywords
error
user
satellite
bdsbas
clock
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
CN202010861987.1A
Other languages
English (en)
Other versions
CN112099063B (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 CN202010861987.1A priority Critical patent/CN112099063B/zh
Publication of CN112099063A publication Critical patent/CN112099063A/zh
Application granted granted Critical
Publication of CN112099063B publication Critical patent/CN112099063B/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/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain
    • 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/20—Integrity monitoring, fault detection or fault isolation of space segment

Abstract

本发明提供了一种用于北斗星基增强用户误差最大投影方向快速搜索的方法，利用当前BDSBAS B1C/B2a增强电文中的轨道/钟差改正数，广播星历，及可由其他机构提供的精密轨道钟差产品，充分考虑误差的几何特性及服务区域限制，达到快速搜索用户误差最大投影方向的目的，并最终基于搜索出的用户误差最大投影方向验证UDRE/DFRE对该误差的包络性能，以评估BDSBAS完好性参数性能。本发明算法执行简单，搜索效率较高，具有较强的工程使用性，将用户误差最大投影方向从服务区域内的面搜索简化成了判断后的边界搜索，大大提高了搜索效率，保证了BDSBAS服务性能监测评估的快速性。

Description

用于北斗星基增强用户误差最大投影方向快速搜索的方法
技术领域
本发明涉及卫星导航增强领域，是北斗星基增强系统(BeiDou Satellite BasedAugmentation System，BDSBAS)性能评估中一种用于用户误差最大投影方向快速搜索的方法。
背景技术
BDSBAS是我国按照国际标准自主建设的星基增强系统(Satellite BasedAugmentation System,SBAS)，通过中国境内分布的监测站，实现对经过我国上空的全球导航卫星系统(Global Navigation Satellite System,GNSS)的完好性监测。
BDSBAS通过地球同步静止卫星(Geosynchronous Earth Orbit,GEO)的B1C频点播发单频慢变改正数(由电文MT24和MT25播发)、快变改正数(由电文MT2～5播发)、电离层格网改正信息(由电文MT18和MT26播发)及用户差分距离误差(UDRE由电文MT2～MT6)、降效参数相关信息(由电文MT7、MT10、MT27和MT28播发)，实现对GPS系统的差分增强，增强服务性能将满足国际民航一类垂直引导进近(APproach with Vertical guidance I，APV-I)指标要求。BDSBAS GEO卫星通过B2a频点提供双频多星座(Dual-Frequency Multi-Constellation,DFMC)星基增强服务，由于在双频定位模式下，用户可通过双频观测量组合自行消除电离层影响，因此B2a不再播发电离层相关的改正数及完好性信息，只播发轨道钟差改正数(由电文MT32播发)、双频测距误差(DFRE由电文MT34～36播发)、降效参数信息(由电文MT37播发)，实现对BDS和GPS的差分增强，DFMC增强服务性能将满足国际民航一类精密进近(CAT-I)性能要求。
BDSBAS B1C单频增强服务中播发的用户差分距离误差(UDRE)和B2a DFMC服务中播发的双频测距误差(DFRE)反映着卫星轨道/钟差经过增强信息改正后的修正误差在服务区域内用户方向上的最大投影值，其应能够以99.9％的概率(3.29σ)对用户误差最大投影进行包络，包络示意图如图1所示。
当前，BDSBAS正处于建设阶段，从用户角度出发对其完好性参数的监测评估成为研究重点。验证完好性参数UDRE/DFRE对服务区域内用户误差最大投影进行99.9％的包络前提便是找到服务区内用户误差投影的最大方向。然而，目前国内外公开文献并没有一种针对BDSBAS用户误差最大投影方向确定的简单、高效方法描述。因此本发明提出了一中实现简单，效率较高的北斗星基增强用户误差最大投影方向快速搜索的方法，以满足BDSBAS建设过程中完好性参数实时监测评估的需求。
发明内容
为了克服现有技术的不足，本发明提供一种用于北斗星基增强用户误差最大投影方向快速搜索的方法，利用当前BDSBAS B1C/B2a增强电文中的轨道/钟差改正数，广播星历，及可由其他机构(International GNSS Service,IGS)提供的精密轨道钟差产品，充分考虑误差的几何特性及服务区域限制，达到快速搜索用户误差最大投影方向的目的，并最终基于搜索出的用户误差最大投影方向验证UDRE/DFRE对该误差的包络性能，以评估BDSBAS完好性参数性能。
本发明解决其技术问题所采用的技术方案的具体步骤是：
步骤1：BDSBAS B1C/B2a增强电文读取；
1)BDSBAS B1C单频快/慢变改正数解算
快变改正数由BDSBAS B1C增强电文MT2-5和MT24播发，计算公式如下：
PRC(t)＝PRCcurrent+RRC(tof)×(t-tof) (1)
若，aii≠0(由电文MT7播发)，则：
Figure BDA0002648445380000021
若，aii＝0，则：
RRC(tof)＝0 (3)
式(1)-(3)中，PRC(t)为当前时刻t的快变改正数；PRCcurrent为最新接收到的快变改正数(电文MT-2-5和MT24播发)；PRCprevious为PRCcurrent之前接收到的快变改正数(电文MT-2-5和MT24播发)；tof为PRCcurrent的参考时刻；RRC(tof)为计算得到的参考时刻tof的距离变化改正数；Δt＝tof-tof，previous；tof，previous为PRCprevious的参考时刻；aii为快变改正数降效因子索引(由电文MT7播发)；由式(1)计算得到的PRC(t)值将直接对当前伪距测量进行相加修正；
慢变改正数由电文MT24和MT25播发，慢变改正数包括卫星时钟慢变改正数和卫星轨道慢变改正数；
卫星时钟慢变改正数计算公式如下：
δΔtSV(t)＝δaf0+δaf1(t-t0)+δafG0 (4)
式(4)中，δΔtSV(t)为解算时刻t的时钟慢变改正数；δαf0为时钟偏差，由电文MT24、25播发；δαf1为时钟偏差变化率，由电文MT24、25播发，速度标识为1；如果速度标识为0该值为0；t0为改正数参考时刻，由电文MT24、25播发，速度模式标识为1；δαfG0为GLONASS卫星改正参数，在电文MT12中播发，针对非GLONASS卫星，该值为0；
卫星星历慢变改正数通过下式计算(坐标系WGS-84ECEF)：
Figure BDA0002648445380000031
当速度标识为0时，式(5)中的速度分量为0，其中
Figure BDA0002648445380000032
为解算时刻t卫星位置三轴改正数；
Figure BDA0002648445380000033
为参考时刻t0卫星位置三轴改正数，由电文MT24、25播发；
Figure BDA0002648445380000034
为参考时刻t0卫星位置三轴改正数变化率，由电文MT24、25播发，速度标示为0时，变化率为0；
2)BDSBAS B2a DMFC轨道钟差改正数解算；
星历改正数解算从BDSBAS B2a增强电文MT32读取星历位置改正信息(参考坐标系WGS-84ECEF)，解算公式如下：
Figure BDA0002648445380000035
其中，
Figure BDA0002648445380000036
为当前时刻t的星历改正数；
Figure BDA0002648445380000037
为参考时刻tD的星历改正数(由电文MT32播发)；
Figure BDA0002648445380000038
为参考时刻tD的星历改正数变化率，由电文MT32播发；
时钟改正数解算公式如下：
Figure BDA0002648445380000039
其中，δΔtSV为当前时刻t的时钟改正数，单位：秒；δB为参考时刻tD的时钟改正数，单位：米(由电文32播发)；
Figure BDA0002648445380000041
为参考时刻tD的时钟改正数变化率，单位：米/秒(由电文32播发)；光速c＝299792458米/秒；
步骤2：BDSBAS B1C/B2a增强轨道/钟差误差求解
BDSBAS B1C/B2a播发的轨道/时钟改正数是对GNSS卫星广播星历解算得到的卫星位置与时钟进行修正，由广播星历计算卫星位置与时钟已有成熟算法，通过广播星历计算得到时刻t的卫星位置与时钟误差分别为
Figure BDA0002648445380000042
Δtsv，bc，则BDSBAS增强后的卫星位置
Figure BDA0002648445380000043
与钟差Δtsv为下式：
Figure BDA0002648445380000044
Δtsv＝Δtsv，bc+δΔtSV (9)
IGS等机构可提供sp3格式与clk格式的GNSS卫星精密轨道与时钟产品文件，其轨道精度可达2cm，时钟精度可达75ps以内，因此以sp3文件与clk文件解算出的GNSS卫星轨道与钟差作为参考真值，计算BDSBAS B1C/B2a增强后的轨道与时钟误差。由sp3文件clk文件参数通过轨道钟差拉格朗日内插方法内插任意时刻GNSS卫星位置与时钟，由sp3产品计算出的卫星位置真值为
Figure BDA0002648445380000045
(坐标参考中心为卫星质心)，由广播星历计算得到的卫星位置
Figure BDA0002648445380000046
坐标参考中心为天线相位中心，在计算轨道误差之前需进行天线相位中心修正，天线相位中心由IGS等机构的ATX文件得到，设由ATX文件得到的天线相位中心偏差(PhaseCenter Offset，PCO)三轴分量为
Figure BDA0002648445380000047
则经过PCO修正后的BDSBAS增强卫星位置如下：
Figure BDA0002648445380000048
式(10)中，T为卫星本体坐标系到ECEF坐标系的3×3姿态转换矩阵，由卫星姿态模型计算得出，可通过GNSS卫星类型、卫星位置速度和UTC时间计算得到，则由式(9)，式(10)和卫星精密位置与钟差可得BDSBAS B1C/B2a增强轨道/钟差误差
Figure BDA0002648445380000051
如下：
Figure BDA0002648445380000052
δΔtsv，sbas＝Δtsv-Δtsv，precise (12)
式(12)中Δtsv，precise为由clk文件解算得到的t时刻卫星精密时钟，δΔtsv，sbas为BDSBAS B1C/B2a增强后的卫星时钟误差；
由式(11)和式(12)可得，在用户位置处，BDSBAS B1C/B2a增强后的卫星轨道/时钟误差的投影Erroruser如下式：
Figure BDA0002648445380000053
式(12)中，euser为1x3用户视线向量，由卫星指向用户所在位置；c为真空中光速；
步骤3：轨道误差矢量方向判断
步骤2中得到GNSS卫星在ECEF坐标系内的位置为
Figure BDA0002648445380000054
轨道误差矢量方向如下式：
Figure BDA0002648445380000055
式(14)中“|| ||”为向量求模运算；eorbit为3x1轨道误差矢量方向；dx，dy，dz分别为矢量三轴分量，由式(14)和GNSS卫星位置，得到轨道误差矢量所在的直线在ECEF坐标系内的表示如下式：
dx·x+dy·y+dz·z+D＝0 (14)
D＝-(dx·xt+dy·yt+dz·zt) (15)
已知地球半径RE＝6378136.6米，则地球球面在ECEF坐标系内的表达式如下：
Figure BDA0002648445380000056
则联立式(14)、式(15)和式(16)，求解误差矢量所在的直线与地球球面是否存在穿刺点及穿刺点的坐标[xp yp zp]T；
穿刺点的求解存在以下两种情况：
情况1：轨道误差矢量所在直线与地球球面存在穿刺点p且穿刺点在BDSBAS服务区域内；
如图3中eorbit1所示，在这种情况下认为矢量eorbit1即是BDSBAS服务区域内用户误差最大投影方向，用户误差最大投影值Erroruser，max如下：
Figure BDA0002648445380000061
算法结束；
情况2：轨道误差矢量所在直线与地球球面不存在穿刺点p，或穿刺点p不在BDSBAS服务区域内；
如图3中eorbit2所示，在该情况下eorbit2虽然为轨道误差最大投影方向，但因其指向位置并不在BDSBAS服务区域内，所以其并不是BDSBAS服务区域内用户误差最大投影方向，在此情况下应当执行步骤4，进行用户误差最大投影方向搜索；
步骤4：BDSBAS服务边界位置遍历
步骤3中的情况2，因为轨道误差所在直线与地球球面穿刺点或穿刺便不再BDSBAS服务区域内，因此需要对BDSBAS服务区域内的位置进行遍历，以寻找服务区域内用户误差投影最大方向，如图4所示。设用户在服务区域内ECEF坐标系下的坐标为[xuser yuser zuser]T，则用户与GNSS卫星的视线向量为euser如下式：
Figure BDA0002648445380000062
则式(13)变换为下式：
Figure BDA0002648445380000063
式(18)中θ为用户实现向量euser与轨道误差矢量eorbit之间的夹角，取值范围为0≤θ≤π；由于在0≤θ≤π的取值范围内，cos(θ)的取值单调递减，因此要使Erroruser取最大值Erroruser，max，则θ取最小值θmin；
如图4所示，设置用户位置在BDSBAS服务区域边界上按固定步长遍历，得到euser1、euser2……euseri个用户视线向量，分别求解处各用户视线向量与轨道误差矢量之间的夹角θ，取其中夹角最小的θmin所对应的视线向量eusermin为用户误差最大投影方向，即可得到用户最误差最大投影误差如下式：
Figure BDA0002648445380000071
通过步骤1-4，即可快速确定北斗星基增强用户误差最大投影方向，计算服务区域范围内的用户最大投影误差，用以评估BDSBAS完好形参数性能。其充分利用了卫星轨道/钟差改正误差的几何特性、BDSBAS服务区域范围限制特性，将用户误差最大投影方向从服务区域内的面搜索简化成了判断后的边界搜索，大大提高了搜索效率，具有较强的工程应用可行性，为BDSBAS服务性能监测提供了新思路、新方法。
本发明的有益效果在于：
1)提出了一种用于北斗星基增强用户误差最大投影方向快速搜索的方法，给出了具体处理流程和实施步骤，算法执行简单，搜索效率较高，具有较强的工程使用性，能够为BDSBAS服务性能监测评估提供理论依据和实施思路。
2)充分利用了卫星轨道/钟差改正误差的几何特性、BDSBAS服务区域范围限制特性，将用户误差最大投影方向从服务区域内的面搜索简化成了判断后的边界搜索，大大提高了搜索效率，保证了BDSBAS服务性能监测评估的快速性。
附图说明
图1为BDSBAS单频UDRE或双频DFRE对用户误差最大投影的包络示意图。
图2为北斗星基增强用户误差最大投影方向快速搜索方法流程图。
图3为BDSBAS B1C/B2a轨道误差矢量与服务区域关系图。
图4为BDSBAS服务边界位置遍历示意图。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
步骤1：BDSBAS B1C/B2a增强电文读取。对实时接收到的BDSBAS B1C/B2a增强电文按照电文类型进行解析，获得当前单频轨道/钟差改正数或双频轨道钟差改正数。
步骤2：BDSBAS B1C/B2a增强轨道/钟差误差求解。以BDSBAS B1C/B2a的轨道/钟差改正数对广播星历轨道钟差进行修正，再与精密轨道/钟差做差求出改正后的B1C/B2a轨道/钟差误差。
步骤3：轨道误差矢量方向判断。对步骤2中的三轴轨道误差所构成的矢量方向进行判断，若该矢量方向延长后与BDSBAS服务区域存在穿刺点，则表明该矢量方向即为服务区域内用户误差最大投影方向，结束；否则，进行步骤4.
步骤4：BDSBAS服务边界位置遍历。以BDSBAS服务区域边界某一点开始，固定步长遍历边界上每个经纬度位置，计算该位置与卫星的视线向量，并计算视线向量与轨道误差矢量方向的夹角，夹角最小时的视线向量即为用户服务区域内误差最大投影方向，结束。
本发明是一种北斗星基增强用户最大误差投影方向快速搜索的方法，具体步骤如图2所示，步骤如下：
步骤1：BDSBAS B1C/B2a增强电文读取
1)BDSBAS B1C单频快/慢变改正数解算
快变改正数由BDSBAS B1C增强电文MT2-5和MT24播发，计算公式如下：
PRC(t)＝PRCcurrent+RRC(tof)×(t-tof) (1)
若，aii≠0(由电文MT7播发)，则：
Figure BDA0002648445380000081
若，aii＝0，则：
RRC(tof)＝0 (3)
式(1)-(3)中，PRC(t)为当前时刻t的快变改正数；PRCcurrent为最新接收到的快变改正数(电文MT-2-5和MT24播发)；PRCprevious为PRCcurrent之前接收到的快变改正数(电文MT-2-5和MT24播发)；tof为PRCcurrent的参考时刻；RRC(tof)为计算得到的参考时刻tof的距离变化改正数；Δt＝tof-tof，previous；tof，previous为PRCprevious的参考时刻；aii为快变改正数降效因子索引(由电文MT7播发)；由式(1)计算得到的PRC(t)值将直接对当前伪距测量进行相加修正。
慢变改正数由电文MT24和MT25播发，慢变改正数包括卫星时钟慢变改正数和卫星轨道慢变改正数；
卫星时钟慢变改正数计算公式如下：
δΔtSV(t)＝δaf0+δaf1(t-t0)+δafG0 (4)
式(4)中，δΔtSV(t)为解算时刻t的时钟慢变改正数；δaf0为时钟偏差(由电文MT24、25播发)；δaf1为时钟偏差变化率(由电文MT24、25播发，速度标识为1；如果速度标识为0该值为0)；t0为改正数参考时刻(由电文MT24、25播发，速度模式标识为1)；δafG0为GLONASS卫星改正参数，在电文MT12中播发，针对非GLONASS卫星，该值为0；
卫星星历慢变改正数通过下式计算(坐标系WGS-84ECEF)：
Figure BDA0002648445380000091
当速度标识为0时，式(5)中的速度分量为0，其中
Figure BDA0002648445380000092
为解算时刻t卫星位置三轴改正数；
Figure BDA0002648445380000093
为参考时刻t0卫星位置三轴改正数(由电文MT24、25播发)；
Figure BDA0002648445380000094
为参考时刻t0卫星位置三轴改正数变化率(由电文MT24、25播发，速度标示为0时，变化率为0)。
2)BDSBAS B2a DMFC轨道钟差改正数解算；
星历改正数解算从BDSBAS B2a增强电文MT32读取星历位置改正信息(参考坐标系WGS-84ECEF)，解算公式如下：
Figure BDA0002648445380000095
其中，
Figure BDA0002648445380000096
为当前时刻t的星历改正数；
Figure BDA0002648445380000097
为参考时刻tD的星历改正数(由电文MT32播发)；
Figure BDA0002648445380000098
为参考时刻tD的星历改正数变化率(由电文MT32播发)。
时钟改正数解算公式如下：
Figure BDA0002648445380000101
其中，δΔtSV为当前时刻t的时钟改正数，单位：秒；δB为参考时刻tD的时钟改正数，单位：米(由电文32播发)；
Figure BDA0002648445380000102
为参考时刻tD的时钟改正数变化率，单位：米/秒(由电文32播发)；光速c＝299792458米/秒；
步骤2：BDSBAS B1C/B2a增强轨道/钟差误差求解
BDSBAS B1C/B2a播发的轨道/时钟改正数是对GNSS卫星广播星历解算得到的卫星位置与时钟进行修正。由广播星历计算卫星位置与时钟已有成熟算法，通过广播星历计算得到时刻t的卫星位置与时钟误差分别为
Figure BDA0002648445380000103
Δtsv，bc，则BDSBAS增强后的卫星位置
Figure BDA0002648445380000104
与钟差Δtsv为下式：
Figure BDA0002648445380000105
Δtsv＝Δtsv，bc+δΔtSV (9)
IGS等机构可提供sp3格式与clk格式的GNSS卫星精密轨道与时钟产品文件，其轨道精度可达2cm，时钟精度可达75ps以内，因此以sp3文件与clk文件解算出的GNSS卫星轨道与钟差作为参考真值，计算BDSBAS B1C/B2a增强后的轨道与时钟误差。由sp3文件clk文件参数通过轨道钟差拉格朗日内插方法内插任意时刻GNSS卫星位置与时钟，由sp3产品计算出的卫星位置真值为
Figure BDA0002648445380000106
(坐标参考中心为卫星质心)，由广播星历计算得到的卫星位置
Figure BDA0002648445380000107
坐标参考中心为天线相位中心，在计算轨道误差之前需进行天线相位中心修正，天线相位中心由IGS等机构的ATX文件得到，设由该文件得到的天线相位中心偏差(PhaseCenter Offset，PCO)三轴分量为
Figure BDA0002648445380000108
则经过PCO修正后的BDSBAS增强卫星位置如下：
Figure BDA0002648445380000111
式(10)中，T为卫星本体坐标系到ECEF坐标系的3×3姿态转换矩阵，由卫星姿态模型计算得出，可通过GNSS卫星类型、卫星位置速度和UTC时间计算得到。则由式(9)，式(10)和卫星精密位置与钟差可得BDSBAS B1C/B2a增强轨道/钟差误差
Figure BDA0002648445380000112
如下：
Figure BDA0002648445380000113
δΔtsv，sbas＝Δtsv-Δtsv，precise (12)
式(12)中Δtsv，precise为由clk文件解算得到的t时刻卫星精密时钟，δΔtsv，sbas为BDSBAS B1C/B2a增强后的卫星时钟误差。
由式(11)和式(12)可得，在用户位置处，BDSBAS B1C/B2a增强后的卫星轨道/时钟误差的投影Erroruser如下式：
Figure BDA0002648445380000114
式(12)中，euser为1x3用户视线向量，由卫星指向用户所在位置；c为真空中光速。
步骤3：轨道误差矢量方向判断
步骤2中得到GNSS卫星在ECEF坐标系内的位置为
Figure BDA0002648445380000116
轨道误差矢量方向如下式：
Figure BDA0002648445380000115
式(14)中“|| ||”为向量求模运算；eorbit为3x1轨道误差矢量方向；dx，dy，dz分别为矢量三轴分量。由式14和GNSS卫星位置，得到轨道误差矢量所在的直线在ECEF坐标系内的表示如下式：
dx·x+dy·y+dz·z+D＝0 (14)
D＝-(dx·xt+dy·yt+dz·zt) (15)
已知地球半径RE＝6378136.6米，则地球球面在ECEF坐标系内的表达式如下：
Figure BDA0002648445380000121
则联立式(14)、式(15)和式(16)可求解误差矢量所在的直线与地球球面是否存在穿刺点，及穿刺点的坐标[xp yp zp]T；
穿刺点的求解存在以下两种情况：
情况1：轨道误差矢量所在直线与地球球面存在穿刺点p且穿刺点在BDSBAS服务区域内；
如图3中eorbit1所示，在这种情况下认为矢量eorbit1即是BDSBAS服务区域内用户误差最大投影方向，用户误差最大投影值Erroruser，max如下：
Figure BDA0002648445380000122
算法结束。
情况2：轨道误差矢量所在直线与地球球面不存在穿刺点p，或穿刺点p不在BDSBAS服务区域内。
如图3中eorbit2所示，在该情况下eorbit2虽然为轨道误差最大投影方向，但因其指向位置并不在BDSBAS服务区域内，所以其并不是BDSBAS服务区域内用户误差最大投影方向，在此情况下应当执行步骤4，进行用户误差最大投影方向搜索。
步骤4：BDSBAS服务边界位置遍历
步骤3中的情况2，因为轨道误差所在直线与地球球面穿刺点或穿刺便不再BDSBAS服务区域内，因此需要对BDSBAS服务区域内的位置进行遍历，以寻找服务区域内用户误差投影最大方向，如图4所示。设用户在服务区域内ECEF坐标系下的坐标为[xuser yuser zuser]T，则用户与GNSS卫星的视线向量为euser如下式：
Figure BDA0002648445380000123
则式(13)变换为下式：
Figure BDA0002648445380000131
式(18)中θ为用户实现向量euser与轨道误差矢量eorbit之间的夹角，取值范围为0≤θ≤π；由于在0≤θ≤π的取值范围内，cos(θ)的取值单调递减，因此要使Erroruser取最大值Erroruser，max，则θ取最小值θmin。
如图4所示，设置用户位置在BDSBAS服务区域边界上按固定步长遍历，得到euser1、euser2……euseri个用户视线向量，分别求解处各用户视线向量与轨道误差矢量之间的夹角θ，取其中夹角最小的θmin所对应的视线向量eusermin为用户误差最大投影方向，即可得到用户最误差最大投影误差如下式：
Figure BDA0002648445380000132
通过步骤1-4，即可快速确定北斗星基增强用户误差最大投影方向，计算服务区域范围内的用户最大投影误差，用以评估BDSBAS完好形参数性能。其充分利用了卫星轨道/钟差改正误差的几何特性、BDSBAS服务区域范围限制特性，将用户误差最大投影方向从服务区域内的面搜索简化成了判断后的边界搜索，大大提高了搜索效率，具有较强的工程应用可行性，为BDSBAS服务性能监测提供了新思路、新方法。

Claims (1)

1.一种用于北斗星基增强用户误差最大投影方向快速搜索的方法，其特征在于包括下述步骤：
步骤1：BDSBAS B1C/B2a增强电文读取；
1)BDSBAS B1C单频快/慢变改正数解算
快变改正数由BDSBAS B1C增强电文MT2-5和MT24播发，计算公式如下：
PRC(t)＝PRCcurrent+RRC(tof)×(t-tof) (1)
若，aii≠0(由电文MT7播发)，则：
Figure FDA0002648445370000011
若，aii＝0，则：
RRC(tof)＝0 (3)
式(1)-(3)中，PRC(t)为当前时刻t的快变改正数；PRCcurrent为最新接收到的快变改正数，由电文MT-2-5和MT24播发；PRCprevious为PRCcurrent之前接收到的快变改正数，电文MT-2-5和MT24播发；tof为PRCcurrent的参考时刻；RRC(tof)为计算得到的参考时刻tof的距离变化改正数；Δt＝tof-tof，previous；tof，previous为PRCprevious的参考时刻；aii为快变改正数降效因子索引，由电文MT7播发；由式(1)计算得到的PRC(t)值将直接对当前伪距测量进行相加修正；
慢变改正数由电文MT24和MT25播发，慢变改正数包括卫星时钟慢变改正数和卫星轨道慢变改正数；
卫星时钟慢变改正数计算公式如下：
δΔtSV(t)＝δaf0+δaf1(t-t0)+δafG0 (4)
式(4)中，δΔtSV(t)为解算时刻t的时钟慢变改正数；δaf0为时钟偏差，由电文MT24、25播发；δaf1为时钟偏差变化率，由电文MT24、25播发，速度标识为1；如果速度标识为0该值为0；t0为改正数参考时刻，由电文MT24、25播发，速度模式标识为1；δafG0为GLONASS卫星改正参数，在电文MT12中播发，针对非GLONASS卫星，该值为0；
卫星星历慢变改正数通过下式在坐标系WGS-84ECEF计算：
Figure FDA0002648445370000021
当速度标识为0时，式(5)中的速度分量为0，其中
Figure FDA0002648445370000022
为解算时刻t卫星位置三轴改正数；
Figure FDA0002648445370000023
为参考时刻t0卫星位置三轴改正数，由电文MT24、25播发；
Figure FDA0002648445370000024
为参考时刻t0卫星位置三轴改正数变化率，由电文MT24、25播发，速度标示为0时，变化率为0；
2)BDSBAS B2a DMFC轨道钟差改正数解算；
星历改正数解算从BDSBAS B2a增强电文MT32读取星历位置改正信息，解算公式如下：
Figure FDA0002648445370000025
其中，
Figure FDA0002648445370000026
为当前时刻t的星历改正数；
Figure FDA0002648445370000027
为参考时刻tD的星历改正数，由电文MT32播发；
Figure FDA0002648445370000028
为参考时刻tD的星历改正数变化率，由电文MT32播发；
时钟改正数解算公式如下：
Figure FDA0002648445370000029
其中，δΔtSV为当前时刻t的时钟改正数，单位：秒；δB为参考时刻tD的时钟改正数，单位：米，由电文32播发；
Figure FDA00026484453700000210
为参考时刻tD的时钟改正数变化率，单位：米/秒，由电文32播发；光速c＝299792458米/秒；
步骤2：BDSBAS B1C/B2a增强轨道/钟差误差求解
BDSBAS B1C/B2a播发的轨道/时钟改正数是对GNSS卫星广播星历解算得到的卫星位置与时钟进行修正，由广播星历计算卫星位置与时钟已有成熟算法，通过广播星历计算得到时刻t的卫星位置与时钟误差分别为
Figure FDA00026484453700000211
Δtsv，bc，则BDSBAS增强后的卫星位置
Figure FDA0002648445370000031
与钟差Δtsv为下式：
Figure FDA0002648445370000032
Δtsv＝Δtsv，bc+δΔtSV (9)
以sp3文件与clk文件解算出的GNSS卫星轨道与钟差作为参考真值，计算BDSBAS B1C/B2a增强后的轨道与时钟误差；由sp3文件clk文件参数通过轨道钟差拉格朗日内插方法内插任意时刻GNSS卫星位置与时钟，由sp3产品计算出的卫星位置真值为
Figure FDA0002648445370000033
坐标参考中心为卫星质心，由广播星历计算得到的卫星位置
Figure FDA0002648445370000034
坐标参考中心为天线相位中心，在计算轨道误差之前需进行天线相位中心修正，天线相位中心由IGS等机构的ATX文件得到，设由ATX文件得到的天线相位中心偏差三轴分量为
Figure FDA0002648445370000035
则经过PCO修正后的BDSBAS增强卫星位置如下：
Figure FDA0002648445370000036
式(10)中，T为卫星本体坐标系到ECEF坐标系的3×3姿态转换矩阵，由卫星姿态模型计算得出，通过GNSS卫星类型、卫星位置速度和UTC时间计算得到，则由式(9)，式(10)和卫星精密位置与钟差可得BDSBAS B1C/B2a增强轨道/钟差误差
Figure FDA0002648445370000037
如下：
Figure FDA0002648445370000038
δΔtsv，sbas＝Δtsv-Δtsv，precise (12)
式(12)中Δtsv，precise为由clk文件解算得到的t时刻卫星精密时钟，δΔtsv，sbas为BDSBASB1C/B2a增强后的卫星时钟误差；
由式(11)和式(12)可得，在用户位置处，BDSBAS B1C/B2a增强后的卫星轨道/时钟误差的投影Erroruser如下式：
Figure FDA0002648445370000041
式(12)中，euser为1x3用户视线向量，由卫星指向用户所在位置；c为真空中光速；
步骤3：轨道误差矢量方向判断
步骤2中得到GNSS卫星在ECEF坐标系内的位置为
Figure FDA0002648445370000042
轨道误差矢量方向如下式：
Figure FDA0002648445370000043
式(14)中“|| ||”为向量求模运算；eorbit为3x1轨道误差矢量方向；dx，dy，dz分别为矢量三轴分量，由式(14)和GNSS卫星位置，得到轨道误差矢量所在的直线在ECEF坐标系内的表示如下式：
dx·x+dy·y+dz·z+D＝0 (14)
D＝-(dx·xt+dy·yt+dz·zt) (15)
已知地球半径RE＝6378136.6米，则地球球面在ECEF坐标系内的表达式如下：
Figure FDA0002648445370000044
则联立式(14)、式(15)和式(16)，求解误差矢量所在的直线与地球球面是否存在穿刺点及穿刺点的坐标[xp yp zp]T；
穿刺点的求解存在以下两种情况：
情况1：轨道误差矢量所在直线与地球球面存在穿刺点p且穿刺点在BDSBAS服务区域内；
矢量eorbit1即是BDSBAS服务区域内用户误差最大投影方向，用户误差最大投影值Erroruser，max如下：
Figure FDA0002648445370000045
算法结束；
情况2：轨道误差矢量所在直线与地球球面不存在穿刺点p，或穿刺点p不在BDSBAS服务区域内；
执行步骤4，进行用户误差最大投影方向搜索；
步骤4：BDSBAS服务边界位置遍历
步骤3中的情况2，对BDSBAS服务区域内的位置进行遍历，以寻找服务区域内用户误差投影最大方向，设用户在服务区域内ECEF坐标系下的坐标为[xuser yuser zuser]T，则用户与GNSS卫星的视线向量为euser如下式：
Figure FDA0002648445370000051
则式(13)变换为下式：
Figure FDA0002648445370000052
式(18)中θ为用户实现向量euser与轨道误差矢量eorbit之间的夹角，取值范围为0≤θ≤π；使Erroruser取最大值Erroruser，max，则θ取最小值θmin；
设置用户位置在BDSBAS服务区域边界上按固定步长遍历，得到euser1、euser2……euseri个用户视线向量，分别求解处各用户视线向量与轨道误差矢量之间的夹角θ，取其中夹角最小的θmin所对应的视线向量eusermin为用户误差最大投影方向，即可得到用户最误差最大投影误差如下式：
Figure FDA0002648445370000053
通过步骤1-4，即可快速确定北斗星基增强用户误差最大投影方向，计算服务区域范围内的用户最大投影误差，用以评估BDSBAS完好形参数性能。
CN202010861987.1A 2020-08-25 2020-08-25 用于北斗星基增强用户误差最大投影方向快速搜索的方法 Active CN112099063B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010861987.1A CN112099063B (zh) 2020-08-25 2020-08-25 用于北斗星基增强用户误差最大投影方向快速搜索的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010861987.1A CN112099063B (zh) 2020-08-25 2020-08-25 用于北斗星基增强用户误差最大投影方向快速搜索的方法

Publications (2)

Publication Number Publication Date
CN112099063A true CN112099063A (zh) 2020-12-18
CN112099063B CN112099063B (zh) 2022-12-27

Family

ID=73754608

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010861987.1A Active CN112099063B (zh) 2020-08-25 2020-08-25 用于北斗星基增强用户误差最大投影方向快速搜索的方法

Country Status (1)

Country Link
CN (1) CN112099063B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113391334A (zh) * 2021-07-14 2021-09-14 武汉大学 基于北斗全球短报文的空间信号测距误差改正数编码方法
CN114152962A (zh) * 2021-11-14 2022-03-08 中国电子科技集团公司第二十研究所 一种星基增强系统服务范围确定的方法
CN114966758A (zh) * 2022-05-24 2022-08-30 中国科学院国家授时中心 基于地球椭球的低轨卫星最差投影瞬时轨道ure的解析方法
CN117335899A (zh) * 2023-10-09 2024-01-02 中国人民解放军32021部队 一种北斗星基增强服务降效程度评估方法
CN117579391A (zh) * 2024-01-16 2024-02-20 中国人民解放军战略支援部队航天工程大学 一种基于sbas认证服务的钟差电文降频优化方法和系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100123622A1 (en) * 2008-11-17 2010-05-20 Neil Harper System and method for determining the location of a mobile device
CN106468774A (zh) * 2016-09-09 2017-03-01 北京航空航天大学 一种应用于星基增强系统的星历星钟改正参数及空间信号完好性参数方法
CN107229061A (zh) * 2017-07-18 2017-10-03 武汉大学 一种基于低轨卫星的星地差分实时精密定位方法
CN110133689A (zh) * 2019-05-24 2019-08-16 中国科学院国家授时中心 自适应用户自主完好性监测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100123622A1 (en) * 2008-11-17 2010-05-20 Neil Harper System and method for determining the location of a mobile device
CN106468774A (zh) * 2016-09-09 2017-03-01 北京航空航天大学 一种应用于星基增强系统的星历星钟改正参数及空间信号完好性参数方法
CN107229061A (zh) * 2017-07-18 2017-10-03 武汉大学 一种基于低轨卫星的星地差分实时精密定位方法
CN110133689A (zh) * 2019-05-24 2019-08-16 中国科学院国家授时中心 自适应用户自主完好性监测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
倪育德 等: "基于WAAS监测数据的精度和完好性参数解算", 《空间科学学报》 *
倪育德 等: "基于WAAS监测数据的精度和完好性参数解算", 《空间科学学报》, 31 December 2019 (2019-12-31), pages 502 - 511 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113391334A (zh) * 2021-07-14 2021-09-14 武汉大学 基于北斗全球短报文的空间信号测距误差改正数编码方法
CN113391334B (zh) * 2021-07-14 2022-10-14 武汉大学 基于北斗全球短报文的空间信号测距误差改正数编码方法
CN114152962A (zh) * 2021-11-14 2022-03-08 中国电子科技集团公司第二十研究所 一种星基增强系统服务范围确定的方法
CN114966758A (zh) * 2022-05-24 2022-08-30 中国科学院国家授时中心 基于地球椭球的低轨卫星最差投影瞬时轨道ure的解析方法
CN117335899A (zh) * 2023-10-09 2024-01-02 中国人民解放军32021部队 一种北斗星基增强服务降效程度评估方法
CN117335899B (zh) * 2023-10-09 2024-04-19 中国人民解放军32021部队 一种北斗星基增强服务降效程度评估方法
CN117579391A (zh) * 2024-01-16 2024-02-20 中国人民解放军战略支援部队航天工程大学 一种基于sbas认证服务的钟差电文降频优化方法和系统
CN117579391B (zh) * 2024-01-16 2024-04-05 中国人民解放军战略支援部队航天工程大学 一种基于sbas认证服务的钟差电文降频优化方法和系统

Also Published As

Publication number Publication date
CN112099063B (zh) 2022-12-27

Similar Documents

Publication Publication Date Title
CN112099063B (zh) 用于北斗星基增强用户误差最大投影方向快速搜索的方法
CN112711048B (zh) 基于北斗三号rdss短报文的ssr传输方法及高精度定位系统
CN103823222B (zh) 通过扩展sps轨道信息进行定位的方法和装置
US8493266B2 (en) Global positioning system (GPS) receiver and method of determining location of GPS receiver
CN102565814B (zh) 卫星导航系统的信号精度和定位服务可用性的评估方法
CN111856534B (zh) 智能终端的双模gnss载波精密单点定位方法及系统
CN1645162A (zh) 确定卫星定位系统中的时间的方法和装置
CN102426017A (zh) 一种基于星敏感器确定载体相对于地理坐标系姿态的方法
CN112285749B (zh) 全球导航卫星系统原始观测数据处理方法、装置及存储介质
Li et al. Improving availability and accuracy of GPS/BDS positioning using QZSS for single receiver
Wang et al. Integrity monitoring of PPP-RTK positioning; Part II: LEO augmentation
Gioia et al. The Galileo High Accuracy Service (HAS): decoding and processing live corrections for code-based positioning
RietDorf et al. Precise positioning in real-time using navigation satellites and telecommunication
CN107505645B (zh) 一种导航定位系统及方法
Anghileri et al. Reduced navigation data for a fast first fix
CN113176596B (zh) 气压高程约束定位方法
Yang et al. Research on real-time precise point positioning method based on short message communication
van Diggelen Assisted GNSS
CN117782080B (zh) 基于PPP-B2b/INS的实时空基导航系统及方法
CN113835108B (zh) 一种定位数据的处理方法及装置、计算设备
CN117055087B (zh) 一种雾霾影响区域的gnss实时定位解算方法
CN116256782B (zh) 一种基于双天线gnss单差算法的多路径误差消除方法
Solé et al. Enhancement study of GPS snapshot techniques
CN117782080A (zh) 基于PPP-B2b/INS的实时空基导航系统及方法
Park et al. PPP-RTK service for port navigation using state space representation in DGNSS medium frequency wave

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