用于北斗星基增强用户误差最大投影方向快速搜索的方法与流程
本发明涉及卫星导航增强领域，是北斗星基增强系统(beidousatellitebasedaugmentationsystem，bdsbas)性能评估中一种用于用户误差最大投影方向快速搜索的方法。
背景技术：
bdsbas是我国按照国际标准自主建设的星基增强系统(satellitebasedaugmentationsystem,sbas)，通过中国境内分布的监测站，实现对经过我国上空的全球导航卫星系统(globalnavigationsatellitesystem,gnss)的完好性监测。
bdsbas通过地球同步静止卫星(geosynchronousearthorbit,geo)的b1c频点播发单频慢变改正数(由电文mt24和mt25播发)、快变改正数(由电文mt2～5播发)、电离层格网改正信息(由电文mt18和mt26播发)及用户差分距离误差(udre由电文mt2～mt6)、降效参数相关信息(由电文mt7、mt10、mt27和mt28播发)，实现对gps系统的差分增强，增强服务性能将满足国际民航一类垂直引导进近(approachwithverticalguidancei，apv-i)指标要求。bdsbasgeo卫星通过b2a频点提供双频多星座(dual-frequencymulti-constellation,dfmc)星基增强服务，由于在双频定位模式下，用户可通过双频观测量组合自行消除电离层影响，因此b2a不再播发电离层相关的改正数及完好性信息，只播发轨道钟差改正数(由电文mt32播发)、双频测距误差(dfre由电文mt34～36播发)、降效参数信息(由电文mt37播发)，实现对bds和gps的差分增强，dfmc增强服务性能将满足国际民航一类精密进近(cat-i)性能要求。
bdsbasb1c单频增强服务中播发的用户差分距离误差(udre)和b2adfmc服务中播发的双频测距误差(dfre)反映着卫星轨道/钟差经过增强信息改正后的修正误差在服务区域内用户方向上的最大投影值，其应能够以99.9％的概率(3.29σ)对用户误差最大投影进行包络，包络示意图如图1所示。
当前，bdsbas正处于建设阶段，从用户角度出发对其完好性参数的监测评估成为研究重点。验证完好性参数udre/dfre对服务区域内用户误差最大投影进行99.9％的包络前提便是找到服务区内用户误差投影的最大方向。然而，目前国内外公开文献并没有一种针对bdsbas用户误差最大投影方向确定的简单、高效方法描述。因此本发明提出了一中实现简单，效率较高的北斗星基增强用户误差最大投影方向快速搜索的方法，以满足bdsbas建设过程中完好性参数实时监测评估的需求。
技术实现要素：
为了克服现有技术的不足，本发明提供一种用于北斗星基增强用户误差最大投影方向快速搜索的方法，利用当前bdsbasb1c/b2a增强电文中的轨道/钟差改正数，广播星历，及可由其他机构(internationalgnssservice,igs)提供的精密轨道钟差产品，充分考虑误差的几何特性及服务区域限制，达到快速搜索用户误差最大投影方向的目的，并最终基于搜索出的用户误差最大投影方向验证udre/dfre对该误差的包络性能，以评估bdsbas完好性参数性能。
本发明解决其技术问题所采用的技术方案的具体步骤是：
步骤1：bdsbasb1c/b2a增强电文读取；
1)bdsbasb1c单频快/慢变改正数解算
快变改正数由bdsbasb1c增强电文mt2-5和mt24播发，计算公式如下：
prc(t)＝prccurrent+rrc(tof)×(t-tof)(1)
若，aii≠0(由电文mt7播发)，则：
若，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)：
当速度标识为0时，式(5)中的速度分量为0，其中
2)bdsbasb2admfc轨道钟差改正数解算；
星历改正数解算从bdsbasb2a增强电文mt32读取星历位置改正信息(参考坐标系wgs-84ecef)，解算公式如下：
其中，
时钟改正数解算公式如下：
其中，δδtsv为当前时刻t的时钟改正数，单位：秒；δb为参考时刻td的时钟改正数，单位：米(由电文32播发)；
步骤2：bdsbasb1c/b2a增强轨道/钟差误差求解
bdsbasb1c/b2a播发的轨道/时钟改正数是对gnss卫星广播星历解算得到的卫星位置与时钟进行修正，由广播星历计算卫星位置与时钟已有成熟算法，通过广播星历计算得到时刻t的卫星位置与时钟误差分别为
δtsv＝δtsv，bc+δδtsv(9)
igs等机构可提供sp3格式与clk格式的gnss卫星精密轨道与时钟产品文件，其轨道精度可达2cm，时钟精度可达75ps以内，因此以sp3文件与clk文件解算出的gnss卫星轨道与钟差作为参考真值，计算bdsbasb1c/b2a增强后的轨道与时钟误差。由sp3文件clk文件参数通过轨道钟差拉格朗日内插方法内插任意时刻gnss卫星位置与时钟，由sp3产品计算出的卫星位置真值为
式(10)中，t为卫星本体坐标系到ecef坐标系的3×3姿态转换矩阵，由卫星姿态模型计算得出，可通过gnss卫星类型、卫星位置速度和utc时间计算得到，则由式(9)，式(10)和卫星精密位置与钟差可得bdsbasb1c/b2a增强轨道/钟差误差
δδtsv，sbas＝δtsv-δtsv，precise(12)
式(12)中δtsv，precise为由clk文件解算得到的t时刻卫星精密时钟，δδtsv，sbas为bdsbasb1c/b2a增强后的卫星时钟误差；
由式(11)和式(12)可得，在用户位置处，bdsbasb1c/b2a增强后的卫星轨道/时钟误差的投影erroruser如下式：
式(12)中，euser为1x3用户视线向量，由卫星指向用户所在位置；c为真空中光速；
步骤3：轨道误差矢量方向判断
步骤2中得到gnss卫星在ecef坐标系内的位置为
式(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坐标系内的表达式如下：
则联立式(14)、式(15)和式(16)，求解误差矢量所在的直线与地球球面是否存在穿刺点及穿刺点的坐标[xpypzp]t；
穿刺点的求解存在以下两种情况：
情况1：轨道误差矢量所在直线与地球球面存在穿刺点p且穿刺点在bdsbas服务区域内；
如图3中eorbit1所示，在这种情况下认为矢量eorbit1即是bdsbas服务区域内用户误差最大投影方向，用户误差最大投影值erroruser，max如下：
算法结束；
情况2：轨道误差矢量所在直线与地球球面不存在穿刺点p，或穿刺点p不在bdsbas服务区域内；
如图3中eorbit2所示，在该情况下eorbit2虽然为轨道误差最大投影方向，但因其指向位置并不在bdsbas服务区域内，所以其并不是bdsbas服务区域内用户误差最大投影方向，在此情况下应当执行步骤4，进行用户误差最大投影方向搜索；
步骤4：bdsbas服务边界位置遍历
步骤3中的情况2，因为轨道误差所在直线与地球球面穿刺点或穿刺便不再bdsbas服务区域内，因此需要对bdsbas服务区域内的位置进行遍历，以寻找服务区域内用户误差投影最大方向，如图4所示。设用户在服务区域内ecef坐标系下的坐标为[xuseryuserzuser]t，则用户与gnss卫星的视线向量为euser如下式：
则式(13)变换为下式：
式(18)中θ为用户实现向量euser与轨道误差矢量eorbit之间的夹角，取值范围为0≤θ≤π；由于在0≤θ≤π的取值范围内，cos(θ)的取值单调递减，因此要使erroruser取最大值erroruser，max，则θ取最小值θmin；
如图4所示，设置用户位置在bdsbas服务区域边界上按固定步长遍历，得到euser1、euser2……euseri个用户视线向量，分别求解处各用户视线向量与轨道误差矢量之间的夹角θ，取其中夹角最小的θmin所对应的视线向量eusermin为用户误差最大投影方向，即可得到用户最误差最大投影误差如下式：
通过步骤1-4，即可快速确定北斗星基增强用户误差最大投影方向，计算服务区域范围内的用户最大投影误差，用以评估bdsbas完好形参数性能。其充分利用了卫星轨道/钟差改正误差的几何特性、bdsbas服务区域范围限制特性，将用户误差最大投影方向从服务区域内的面搜索简化成了判断后的边界搜索，大大提高了搜索效率，具有较强的工程应用可行性，为bdsbas服务性能监测提供了新思路、新方法。
本发明的有益效果在于：
1)提出了一种用于北斗星基增强用户误差最大投影方向快速搜索的方法，给出了具体处理流程和实施步骤，算法执行简单，搜索效率较高，具有较强的工程使用性，能够为bdsbas服务性能监测评估提供理论依据和实施思路。
2)充分利用了卫星轨道/钟差改正误差的几何特性、bdsbas服务区域范围限制特性，将用户误差最大投影方向从服务区域内的面搜索简化成了判断后的边界搜索，大大提高了搜索效率，保证了bdsbas服务性能监测评估的快速性。
附图说明
图1为bdsbas单频udre或双频dfre对用户误差最大投影的包络示意图。
图2为北斗星基增强用户误差最大投影方向快速搜索方法流程图。
图3为bdsbasb1c/b2a轨道误差矢量与服务区域关系图。
图4为bdsbas服务边界位置遍历示意图。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
步骤1：bdsbasb1c/b2a增强电文读取。对实时接收到的bdsbasb1c/b2a增强电文按照电文类型进行解析，获得当前单频轨道/钟差改正数或双频轨道钟差改正数。
步骤2：bdsbasb1c/b2a增强轨道/钟差误差求解。以bdsbasb1c/b2a的轨道/钟差改正数对广播星历轨道钟差进行修正，再与精密轨道/钟差做差求出改正后的b1c/b2a轨道/钟差误差。
步骤3：轨道误差矢量方向判断。对步骤2中的三轴轨道误差所构成的矢量方向进行判断，若该矢量方向延长后与bdsbas服务区域存在穿刺点，则表明该矢量方向即为服务区域内用户误差最大投影方向，结束；否则，进行步骤4.
步骤4：bdsbas服务边界位置遍历。以bdsbas服务区域边界某一点开始，固定步长遍历边界上每个经纬度位置，计算该位置与卫星的视线向量，并计算视线向量与轨道误差矢量方向的夹角，夹角最小时的视线向量即为用户服务区域内误差最大投影方向，结束。
本发明是一种北斗星基增强用户最大误差投影方向快速搜索的方法，具体步骤如图2所示，步骤如下：
步骤1：bdsbasb1c/b2a增强电文读取
1)bdsbasb1c单频快/慢变改正数解算
快变改正数由bdsbasb1c增强电文mt2-5和mt24播发，计算公式如下：
prc(t)＝prccurrent+rrc(tof)×(t-tof)(1)
若，aii≠0(由电文mt7播发)，则：
若，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)：
当速度标识为0时，式(5)中的速度分量为0，其中
2)bdsbasb2admfc轨道钟差改正数解算；
星历改正数解算从bdsbasb2a增强电文mt32读取星历位置改正信息(参考坐标系wgs-84ecef)，解算公式如下：
其中，
时钟改正数解算公式如下：
其中，δδtsv为当前时刻t的时钟改正数，单位：秒；δb为参考时刻td的时钟改正数，单位：米(由电文32播发)；
步骤2：bdsbasb1c/b2a增强轨道/钟差误差求解
bdsbasb1c/b2a播发的轨道/时钟改正数是对gnss卫星广播星历解算得到的卫星位置与时钟进行修正。由广播星历计算卫星位置与时钟已有成熟算法，通过广播星历计算得到时刻t的卫星位置与时钟误差分别为
δtsv＝δtsv，bc+δδtsv(9)
igs等机构可提供sp3格式与clk格式的gnss卫星精密轨道与时钟产品文件，其轨道精度可达2cm，时钟精度可达75ps以内，因此以sp3文件与clk文件解算出的gnss卫星轨道与钟差作为参考真值，计算bdsbasb1c/b2a增强后的轨道与时钟误差。由sp3文件clk文件参数通过轨道钟差拉格朗日内插方法内插任意时刻gnss卫星位置与时钟，由sp3产品计算出的卫星位置真值为
式(10)中，t为卫星本体坐标系到ecef坐标系的3×3姿态转换矩阵，由卫星姿态模型计算得出，可通过gnss卫星类型、卫星位置速度和utc时间计算得到。则由式(9)，式(10)和卫星精密位置与钟差可得bdsbasb1c/b2a增强轨道/钟差误差
δδtsv，sbas＝δtsv-δtsv，precise(12)
式(12)中δtsv，precise为由clk文件解算得到的t时刻卫星精密时钟，δδtsv，sbas为bdsbasb1c/b2a增强后的卫星时钟误差。
由式(11)和式(12)可得，在用户位置处，bdsbasb1c/b2a增强后的卫星轨道/时钟误差的投影erroruser如下式：
式(12)中，euser为1x3用户视线向量，由卫星指向用户所在位置；c为真空中光速。
步骤3：轨道误差矢量方向判断
步骤2中得到gnss卫星在ecef坐标系内的位置为
式(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坐标系内的表达式如下：
则联立式(14)、式(15)和式(16)可求解误差矢量所在的直线与地球球面是否存在穿刺点，及穿刺点的坐标[xpypzp]t；
穿刺点的求解存在以下两种情况：
情况1：轨道误差矢量所在直线与地球球面存在穿刺点p且穿刺点在bdsbas服务区域内；
如图3中eorbit1所示，在这种情况下认为矢量eorbit1即是bdsbas服务区域内用户误差最大投影方向，用户误差最大投影值erroruser，max如下：
算法结束。
情况2：轨道误差矢量所在直线与地球球面不存在穿刺点p，或穿刺点p不在bdsbas服务区域内。
如图3中eorbit2所示，在该情况下eorbit2虽然为轨道误差最大投影方向，但因其指向位置并不在bdsbas服务区域内，所以其并不是bdsbas服务区域内用户误差最大投影方向，在此情况下应当执行步骤4，进行用户误差最大投影方向搜索。
步骤4：bdsbas服务边界位置遍历
步骤3中的情况2，因为轨道误差所在直线与地球球面穿刺点或穿刺便不再bdsbas服务区域内，因此需要对bdsbas服务区域内的位置进行遍历，以寻找服务区域内用户误差投影最大方向，如图4所示。设用户在服务区域内ecef坐标系下的坐标为[xuseryuserzuser]t，则用户与gnss卫星的视线向量为euser如下式：
则式(13)变换为下式：
式(18)中θ为用户实现向量euser与轨道误差矢量eorbit之间的夹角，取值范围为0≤θ≤π；由于在0≤θ≤π的取值范围内，cos(θ)的取值单调递减，因此要使erroruser取最大值erroruser，max，则θ取最小值θmin。
如图4所示，设置用户位置在bdsbas服务区域边界上按固定步长遍历，得到euser1、euser2……euseri个用户视线向量，分别求解处各用户视线向量与轨道误差矢量之间的夹角θ，取其中夹角最小的θmin所对应的视线向量eusermin为用户误差最大投影方向，即可得到用户最误差最大投影误差如下式：
通过步骤1-4，即可快速确定北斗星基增强用户误差最大投影方向，计算服务区域范围内的用户最大投影误差，用以评估bdsbas完好形参数性能。其充分利用了卫星轨道/钟差改正误差的几何特性、bdsbas服务区域范围限制特性，将用户误差最大投影方向从服务区域内的面搜索简化成了判断后的边界搜索，大大提高了搜索效率，具有较强的工程应用可行性，为bdsbas服务性能监测提供了新思路、新方法。
- 还没有人留言评论。精彩留言会获得点赞！