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

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

Info

Publication number
CN112099063B
CN112099063B CN202010861987.1A CN202010861987A CN112099063B CN 112099063 B CN112099063 B CN 112099063B CN 202010861987 A CN202010861987 A CN 202010861987A CN 112099063 B CN112099063 B CN 112099063B
Authority
CN
China
Prior art keywords
error
user
satellite
clock
bdsbas
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
CN202010861987.1A
Other languages
English (en)
Other versions
CN112099063A (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 GDA0003807829210000021
若，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 GDA0003807829210000031
当速度标识为0时，式(5)中的速度分量为0，其中
Figure GDA0003807829210000032
为解算时刻t卫星位置三轴改正数；
Figure GDA0003807829210000033
为参考时刻t0卫星位置三轴改正数，由电文MT24、25播发；
Figure GDA0003807829210000034
为参考时刻t0卫星位置三轴改正数变化率，由电文MT24、25播发，速度标示为0时，变化率为0；
2)BDSBAS B2a DMFC轨道钟差改正数解算；
星历改正数解算从BDSBAS B2a增强电文MT32读取星历位置改正信息(参考坐标系WGS-84ECEF)，解算公式如下：
Figure GDA0003807829210000035
其中，
Figure GDA0003807829210000036
为当前时刻t的星历改正数；
Figure GDA0003807829210000037
为参考时刻tD的星历改正数(由电文MT32播发)；
Figure GDA0003807829210000038
为参考时刻tD的星历改正数变化率，由电文MT32播发；
时钟改正数解算公式如下：
Figure GDA0003807829210000039
其中，δΔtSV为当前时刻t的时钟改正数，单位：秒；δB为参考时刻tD的时钟改正数，单位：米(由电文32播发)；
Figure GDA0003807829210000041
为参考时刻tD的时钟改正数变化率，单位：米/秒(由电文32播发)；光速c＝299792458米/秒；
步骤2：BDSBAS B1C/B2a增强轨道/钟差误差求解
BDSBAS B1C/B2a播发的轨道/时钟改正数是对GNSS卫星广播星历解算得到的卫星位置与时钟进行修正，由广播星历计算卫星位置与时钟已有成熟算法，通过广播星历计算得到时刻t的卫星位置与时钟误差分别为
Figure GDA0003807829210000042
Δtsv，bc，则BDSBAS增强后的卫星位置
Figure GDA0003807829210000043
与钟差Δtsv为下式：
Figure GDA0003807829210000044
Δtsv＝Δtsv，bc+δΔtSV (9)
IGS等机构可提供sp3格式与clk格式的GNSS卫星精密轨道与时钟产品文件，其轨道精度可达2cm，时钟精度可达75ps以内，因此以sp3文件与clk文件解算出的GNSS 卫星轨道与钟差作为参考真值，计算BDSBAS B1C/B2a增强后的轨道与时钟误差。由sp3文件clk文件参数通过轨道钟差拉格朗日内插方法内插任意时刻GNSS卫星位置与时钟，由sp3产品计算出的卫星位置真值为
Figure GDA0003807829210000045
(坐标参考中心为卫星质心)，由广播星历计算得到的卫星位置
Figure GDA0003807829210000046
坐标参考中心为天线相位中心，在计算轨道误差之前需进行天线相位中心修正，天线相位中心由IGS等机构的ATX文件得到，设由ATX 文件得到的天线相位中心偏差(PhaseCenter Offset，PCO)三轴分量为
Figure GDA0003807829210000047
则经过 PCO修正后的BDSBAS增强卫星位置如下：
Figure GDA0003807829210000048
式(10)中，T为卫星本体坐标系到ECEF坐标系的3×3姿态转换矩阵，由卫星姿态模型计算得出，可通过GNSS卫星类型、卫星位置速度和UTC时间计算得到，则由式(9)，式(10)和卫星精密位置与钟差可得BDSBAS B1C/B2a增强轨道/钟差误差
Figure GDA0003807829210000051
如下：
Figure GDA0003807829210000052
δΔtsv，sbas＝Δtsv-Δtsv，precise (12)
式(12)中Δtsv，precise为由clk文件解算得到的t时刻卫星精密时钟，δΔtsv，sbas为BDSBAS B1C/B2a增强后的卫星时钟误差；
由式(11)和式(12)可得，在用户位置处，BDSBAS B1C/B2a增强后的卫星轨道/时钟误差的投影Erroruser如下式：
Figure GDA0003807829210000053
式(13)中，euser为1x3用户视线向量，由卫星指向用户所在位置；c为真空中光速；
步骤3：轨道误差矢量方向判断
步骤2中得到GNSS卫星在ECEF坐标系内的位置为
Figure GDA0003807829210000054
轨道误差矢量方向如下式：
Figure GDA0003807829210000055
式(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 GDA0003807829210000056
则联立式(14)、式(15)和式(16)，求解误差矢量所在的直线与地球球面是否存在穿刺点及穿刺点的坐标[xp yp zp]T；
穿刺点的求解存在以下两种情况：
情况1：轨道误差矢量所在直线与地球球面存在穿刺点p且穿刺点在BDSBAS服务区域内；
如图3中eorbit1所示，在这种情况下认为矢量eorbit1即是BDSBAS服务区域内用户误差最大投影方向，用户误差最大投影值Erroruser，max如下：
Figure GDA0003807829210000061
算法结束；
情况2：轨道误差矢量所在直线与地球球面不存在穿刺点p，或穿刺点p不在BDSBAS服务区域内；
如图3中eorbit2所示，在该情况下eorbit2虽然为轨道误差最大投影方向，但因其指向位置并不在BDSBAS服务区域内，所以其并不是BDSBAS服务区域内用户误差最大投影方向，在此情况下应当执行步骤4，进行用户误差最大投影方向搜索；
步骤4：BDSBAS服务边界位置遍历
步骤3中的情况2，因为轨道误差所在直线与地球球面穿刺点或穿刺便不再BDSBAS服务区域内，因此需要对BDSBAS服务区域内的位置进行遍历，以寻找服务区域内用户误差投影最大方向，如图4所示。设用户在服务区域内ECEF坐标系下的坐标为[xuser yuserzuser]T，则用户与GNSS卫星的视线向量为euser如下式：
Figure GDA0003807829210000062
则式(13)变换为下式：
Figure GDA0003807829210000063
式(18)中θ为用户实现向量euser与轨道误差矢量eorbit之间的夹角，取值范围为 0≤θ≤π；由于在0≤θ≤π的取值范围内，cos(θ)的取值单调递减，因此要使Erroruser取最大值Erroruser，max，则θ取最小值θmin；
如图4所示，设置用户位置在BDSBAS服务区域边界上按固定步长遍历，得到 euser1、euser2......euseri个用户视线向量，分别求解处各用户视线向量与轨道误差矢量之间的夹角θ，取其中夹角最小的θmin所对应的视线向量eusermin为用户误差最大投影方向，即可得到用户最误差最大投影误差如下式：
Figure GDA0003807829210000071
通过步骤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 GDA0003807829210000081
若，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-84 ECEF)：
Figure GDA0003807829210000091
当速度标识为0时，式(5)中的速度分量为0，其中
Figure GDA0003807829210000092
为解算时刻t卫星位置三轴改正数；
Figure GDA0003807829210000093
为参考时刻t0卫星位置三轴改正数(由电文MT24、25播发)；
Figure GDA0003807829210000094
为参考时刻t0卫星位置三轴改正数变化率(由电文MT24、25播发，速度标示为0时，变化率为0)。
2)BDSBAS B2a DMFC轨道钟差改正数解算；
星历改正数解算从BDSBAS B2a增强电文MT32读取星历位置改正信息(参考坐标系WGS-84ECEF)，解算公式如下：
Figure GDA0003807829210000095
其中，
Figure GDA0003807829210000096
为当前时刻t的星历改正数；
Figure GDA0003807829210000097
为参考时刻tD的星历改正数(由电文MT32播发)；
Figure GDA0003807829210000098
为参考时刻tD的星历改正数变化率(由电文MT32播发)。
时钟改正数解算公式如下：
Figure GDA0003807829210000101
其中，δΔtSV为当前时刻t的时钟改正数，单位：秒；δB为参考时刻tD的时钟改正数，单位：米(由电文32播发)；
Figure GDA0003807829210000102
为参考时刻tD的时钟改正数变化率，单位：米/秒(由电文32播发)；光速c＝299792458米/秒；
步骤2：BDSBAS B1C/B2a增强轨道/钟差误差求解
BDSBAS B1C/B2a播发的轨道/时钟改正数是对GNSS卫星广播星历解算得到的卫星位置与时钟进行修正。由广播星历计算卫星位置与时钟已有成熟算法，通过广播星历计算得到时刻t的卫星位置与时钟误差分别为
Figure GDA0003807829210000103
Δtsv，bc，则BDSBAS增强后的卫星位置
Figure GDA0003807829210000104
与钟差Δtsv为下式：
Figure GDA0003807829210000105
Δtsv＝Δtsv，bc+δΔtSV (9)
IGS等机构可提供sp3格式与clk格式的GNSS卫星精密轨道与时钟产品文件，其轨道精度可达2cm，时钟精度可达75ps以内，因此以sp3文件与clk文件解算出的GNSS 卫星轨道与钟差作为参考真值，计算BDSBAS B1C/B2a增强后的轨道与时钟误差。由 sp3文件clk文件参数通过轨道钟差拉格朗日内插方法内插任意时刻GNSS卫星位置与时钟，由sp3产品计算出的卫星位置真值为
Figure GDA0003807829210000106
(坐标参考中心为卫星质心)，由广播星历计算得到的卫星位置
Figure GDA0003807829210000107
坐标参考中心为天线相位中心，在计算轨道误差之前需进行天线相位中心修正，天线相位中心由IGS等机构的ATX文件得到，设由该文件得到的天线相位中心偏差(PhaseCenter Offset，PCO)三轴分量为
Figure GDA0003807829210000108
则经过 PCO修正后的BDSBAS增强卫星位置如下：
Figure GDA0003807829210000111
式(10)中，T为卫星本体坐标系到ECEF坐标系的3×3姿态转换矩阵，由卫星姿态模型计算得出，可通过GNSS卫星类型、卫星位置速度和UTC时间计算得到。则由式(9)，式(10)和卫星精密位置与钟差可得BDSBAS B1C/B2a增强轨道/钟差误差
Figure GDA0003807829210000112
如下：
Figure GDA0003807829210000113
δΔtsv，sbas＝Δtsv-Δtsv，precise (12)
式(12)中Δtsv，precise为由clk文件解算得到的t时刻卫星精密时钟，δΔtsv，sbas为BDSBAS B1C/B2a增强后的卫星时钟误差。
由式(11)和式(12)可得，在用户位置处，BDSBAS B1C/B2a增强后的卫星轨道/时钟误差的投影Erroruser如下式：
Figure GDA0003807829210000114
式(13)中，euser为1x3用户视线向量，由卫星指向用户所在位置；c为真空中光速。
步骤3：轨道误差矢量方向判断
步骤2中得到GNSS卫星在ECEF坐标系内的位置为
Figure GDA0003807829210000115
轨道误差矢量方向如下式：
Figure GDA0003807829210000116
式(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 GDA0003807829210000121
则联立式(14)、式(15)和式(16)可求解误差矢量所在的直线与地球球面是否存在穿刺点，及穿刺点的坐标[xp yp zp]T；
穿刺点的求解存在以下两种情况：
情况1：轨道误差矢量所在直线与地球球面存在穿刺点p且穿刺点在BDSBAS服务区域内；
如图3中eorbit1所示，在这种情况下认为矢量eorbit1即是BDSBAS服务区域内用户误差最大投影方向，用户误差最大投影值Erroruser，max如下：
Figure GDA0003807829210000122
算法结束。
情况2：轨道误差矢量所在直线与地球球面不存在穿刺点p，或穿刺点p不在BDSBAS服务区域内。
如图3中eorbit2所示，在该情况下eorbit2虽然为轨道误差最大投影方向，但因其指向位置并不在BDSBAS服务区域内，所以其并不是BDSBAS服务区域内用户误差最大投影方向，在此情况下应当执行步骤4，进行用户误差最大投影方向搜索。
步骤4：BDSBAS服务边界位置遍历
步骤3中的情况2，因为轨道误差所在直线与地球球面穿刺点或穿刺便不再BDSBAS服务区域内，因此需要对BDSBAS服务区域内的位置进行遍历，以寻找服务区域内用户误差投影最大方向，如图4所示。设用户在服务区域内ECEF坐标系下的坐标为[xuser yuserzuser]T，则用户与GNSS卫星的视线向量为euser如下式：
Figure GDA0003807829210000123
则式(13)变换为下式：
Figure GDA0003807829210000131
式(18)中θ为用户实现向量euser与轨道误差矢量eorbit之间的夹角，取值范围为 0≤θ≤π；由于在0≤θ≤π的取值范围内，cos(θ)的取值单调递减，因此要使Erroruser取最大值Erroruser，max，则θ取最小值θmin。
如图4所示，设置用户位置在BDSBAS服务区域边界上按固定步长遍历，得到 euser1、euser2......euseri个用户视线向量，分别求解处各用户视线向量与轨道误差矢量之间的夹角θ，取其中夹角最小的θmin所对应的视线向量eusermin为用户误差最大投影方向，即可得到用户最误差最大投影误差如下式：
Figure GDA0003807829210000132
通过步骤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 FDA0003807829200000011
若，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 FDA0003807829200000021
当速度标识为0时，式(5)中的速度分量为0，其中
Figure FDA0003807829200000022
为解算时刻t卫星位置三轴改正数；
Figure FDA0003807829200000023
为参考时刻t0卫星位置三轴改正数，由电文MT24、25播发；
Figure FDA0003807829200000024
为参考时刻t0卫星位置三轴改正数变化率，由电文MT24、25播发，速度标示为0时，变化率为0；
2)BDSBAS B2a DMFC轨道钟差改正数解算；
星历改正数解算从BDSBAS B2a增强电文MT32读取星历位置改正信息，解算公式如下：
Figure FDA0003807829200000025
其中，
Figure FDA0003807829200000026
为当前时刻t的星历改正数；
Figure FDA0003807829200000027
为参考时刻tD的星历改正数，由电文MT32播发；
Figure FDA0003807829200000028
为参考时刻tD的星历改正数变化率，由电文MT32播发；
时钟改正数解算公式如下：
Figure FDA0003807829200000029
其中，δΔtSV为当前时刻t的时钟改正数，单位：秒；δB为参考时刻tD的时钟改正数，单位：米，由电文32播发；
Figure FDA00038078292000000210
为参考时刻tD的时钟改正数变化率，单位：米/秒，由电文32播发；光速c＝299792458米/秒；
步骤2：BDSBAS B1C/B2a增强轨道/钟差误差求解
BDSBAS B1C/B2a播发的轨道/时钟改正数是对GNSS卫星广播星历解算得到的卫星位置与时钟进行修正，由广播星历计算卫星位置与时钟已有成熟算法，通过广播星历计算得到时刻t的卫星位置与时钟误差分别为
Figure FDA00038078292000000211
Δtsv，bc，则BDSBAS增强后的卫星位置
Figure FDA0003807829200000031
与钟差Δtsv为下式：
Figure FDA0003807829200000032
Δtsv＝Δtsv，bc+δΔtsv (9)
以sp3文件与clk文件解算出的GNSS卫星轨道与钟差作为参考真值，计算BDSBAS B1C/B2a增强后的轨道与时钟误差；由sp3文件clk文件参数通过轨道钟差拉格朗日内插方法内插任意时刻GNSS卫星位置与时钟，由sp3产品计算出的卫星位置真值为
Figure FDA0003807829200000033
坐标参考中心为卫星质心，由广播星历计算得到的卫星位置
Figure FDA0003807829200000034
坐标参考中心为天线相位中心，在计算轨道误差之前需进行天线相位中心修正，天线相位中心由IGS机构的ATX文件得到，设由ATX文件得到的天线相位中心偏差三轴分量为
Figure FDA0003807829200000035
则经过PCO修正后的BDSBAS增强卫星位置如下：
Figure FDA0003807829200000036
式(10)中，T为卫星本体坐标系到ECEF坐标系的3×3姿态转换矩阵，由卫星姿态模型计算得出，通过GNSS卫星类型、卫星位置速度和UTC时间计算得到，则由式(9)，式(10)和卫星精密位置与钟差可得BDSBAS B1C/B2a增强轨道/钟差误差
Figure FDA0003807829200000037
如下：
Figure FDA0003807829200000038
δΔtsv，sbas＝Δtsv-Δtsv，precise (12)
式(12)中Δtsv，precise为由clk文件解算得到的t时刻卫星精密时钟，δΔtsv，sbas为BDSBASB1C/B2a增强后的卫星时钟误差；
由式(11)和式(12)可得，在用户位置处，BDSBAS B1C/B2a增强后的卫星轨道/时钟误差的投影Erroruser如下式：
Figure FDA0003807829200000041
式(13)中，euser为1x3用户视线向量，由卫星指向用户所在位置；c为真空中光速；
步骤3：轨道误差矢量方向判断
步骤2中得到GNSS卫星在ECEF坐标系内的位置为
Figure FDA0003807829200000042
轨道误差矢量方向如下式：
Figure FDA0003807829200000043
式(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 FDA0003807829200000044
则联立式(14)、式(15)和式(16)，求解误差矢量所在的直线与地球球面是否存在穿刺点及穿刺点的坐标[xp yp zp]T；
穿刺点的求解存在以下两种情况：
情况1：轨道误差矢量所在直线与地球球面存在穿刺点p且穿刺点在BDSBAS服务区域内；
矢量eorbit1即是BDSBAS服务区域内用户误差最大投影方向，用户误差最大投影值Erroruser，max如下：
Figure FDA0003807829200000045
算法结束；
情况2：轨道误差矢量所在直线与地球球面不存在穿刺点p，或穿刺点p不在BDSBAS服务区域内；
执行步骤4，进行用户误差最大投影方向搜索；
步骤4：BDSBAS服务边界位置遍历
步骤3中的情况2，对BDSBAS服务区域内的位置进行遍历，以寻找服务区域内用户误差投影最大方向，设用户在服务区域内ECEF坐标系下的坐标为[xuser yuser zuser]T，则用户与GNSS卫星的视线向量为euser如下式：
Figure FDA0003807829200000051
则式(13)变换为下式：
Figure FDA0003807829200000052
式(18)中θ为用户实现向量euser与轨道误差矢量eorbit之间的夹角，取值范围为0≤θ≤π；使Erroruser取最大值Erroruser，max，则θ取最小值θmin；
设置用户位置在BDSBAS服务区域边界上按固定步长遍历，得到euser1、euser2......euseri个用户视线向量，分别求解处各用户视线向量与轨道误差矢量之间的夹角θ，取其中夹角最小的θmin所对应的视线向量eusermin为用户误差最大投影方向，即可得到用户最误差最大投影误差如下式：
Figure FDA0003807829200000053
通过步骤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 CN112099063A (zh) 2020-12-18
CN112099063B true 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)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113391334B (zh) * 2021-07-14 2022-10-14 武汉大学 基于北斗全球短报文的空间信号测距误差改正数编码方法
CN114966758B (zh) * 2022-05-24 2023-05-26 中国科学院国家授时中心 基于地球椭球的低轨卫星最差投影瞬时轨道ure的解析方法
CN117579391B (zh) * 2024-01-16 2024-04-05 中国人民解放军战略支援部队航天工程大学 一种基于sbas认证服务的钟差电文降频优化方法和系统

Citations (3)

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

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8125377B2 (en) * 2008-11-17 2012-02-28 Andrew Llc System and method for determining the location of a mobile device

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
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 (1)

* Cited by examiner, † Cited by third party
Title
基于WAAS监测数据的精度和完好性参数解算;倪育德 等;《空间科学学报》;20191231;第502-511页 *

Also Published As

Publication number Publication date
CN112099063A (zh) 2020-12-18

Similar Documents

Publication Publication Date Title
CN112099063B (zh) 用于北斗星基增强用户误差最大投影方向快速搜索的方法
CN112711048B (zh) 基于北斗三号rdss短报文的ssr传输方法及高精度定位系统
CN103823222B (zh) 通过扩展sps轨道信息进行定位的方法和装置
CN105182384A (zh) 一种双模实时伪距差分定位系统和伪距改正数据生成方法
CN111856534B (zh) 智能终端的双模gnss载波精密单点定位方法及系统
US9933523B2 (en) Systems and methods to enhance reliability of measured position data
Santerre et al. Single point positioning using GPS, GLONASS and BeiDou satellites
CN102426017A (zh) 一种基于星敏感器确定载体相对于地理坐标系姿态的方法
CN112285749B (zh) 全球导航卫星系统原始观测数据处理方法、装置及存储介质
US6844847B2 (en) Method and device for instantaneous determination of orientation, based on satellite positioning signals
CN101419274B (zh) 电离层延迟误差的获取方法及系统
CN110988934A (zh) 多模式接收机星基增强技术装置及处理方法
CN112596077B (zh) 一种针对终端载体为低轨卫星的卫星导航信号仿真方法
CN115951378B (zh) 一种基于北斗星基增强信息的自适应信息融合定位方法
CN106814376B (zh) 一种快速精确厘米级单点定位方法
KR102057547B1 (ko) Lte 기반 이동 통신 기지국을 이용한 이동국의 위치 보정 방법
US11693120B2 (en) System and method for providing GNSS corrections
CN107505645B (zh) 一种导航定位系统及方法
JP2010112759A (ja) 移動体位置測位装置
JP3611526B2 (ja) 衛星測位システム、その地上局及び地上端末
Yang et al. Research on real-time precise point positioning method based on short message communication
CN113176596A (zh) 气压高程约束定位方法
JP2017129555A (ja) 衛星航法システムにおける測位誤差の補正方法及びその装置
Januszewski How the ionosphere affects positioning solution using terrestrial and satellite navigation systems?
CN116256782B (zh) 一种基于双天线gnss单差算法的多路径误差消除方法

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