刘瑞华,丁其金

(中国民航大学 电子信息与自动化学院,天津300300)

0 引 言

北斗卫星导航系统(BDS)是我国自主建设、独立运行的卫星导航系统,正在逐步完善,2020年前后,将向全球提供服务[1]。其应用范围也越来越广,在农业、交通及航空等各个行业中发挥着越来越重要的作用。因此,对北斗接收机的需求也越来越大,各导航设备研制单位竞相展开了北斗接收机的研制工作。虽然不同的应用环境,对北斗接收机要求的侧重点不同,但接收机位置解算都是必须考虑的部分。本文对北斗接收机的位置解算过程进行分析研究,并结合MATLAB软件,分别使用北斗空间控制接口文件所提供的算法解算北斗卫星位置和最小二乘法解算北斗接收机位置,并进行精度分析,验证了算法的正确性和可行性。

1 利用广播星历解算卫星位置

北斗卫星位置解算是接收机位置解算中至关重要的一部分,因为卫星位置是否准确直接影响了接收机的定位精度。北斗卫星导航系统采用三种轨道卫星,即地球同步轨道 (GEO)、倾斜地球同步轨道(Inclined IGSO)和中地球轨道(MEO)组成的混合星座[1],因此,北斗卫星位置的解算方法与GPS略有不同。其中IGSO、MEO卫星计算方法相同,GEO卫星因其轨道倾角和偏心率均很小,需采用坐标系旋转的方法增大轨道倾角以提高卫星位置的计算精度[2-3]。北斗官方文件[4]给出了北斗卫星位置解算的星历参数用户算法,这里不再赘述,图1为算法流程框图,表1示出了星历参数及其定义。

表1 星历参数及其定义

本节使用诺瓦泰FlexPark6北斗接收机2018年3月2日2点至23点采集到的实测星历数据,按照文献[4]中的算法进行北斗卫星位置解算,并将解算出的卫星位置分别和精密星历以及诺瓦泰 FlexPark6北斗接收机所提供的卫星位置进行对比,分析所解算出的卫星位置精度。本文使用武汉大学提供的代号为WUM的精密星历,径向误差优于10 cm,比广播星历精度高出1-2个数量级[5];根据诺瓦泰给出的FlexPark6技术指标,该接收机水平定位精度达到1.5 m(RMS),时间精度为20 ns,因此均可以作为真值来评估所解算卫星位置的精度。

按照文献[4]中的算法,进行北斗卫星位置解算,在解算过程中需注意以下几个问题:

1) 在计算时间差tk时,必须考虑周变换的开始或结束。

2) 按照上述算法所解算的卫星位置是卫星在信号发射时刻的位置,而北斗接收机位置解算过程中所使用的却是信号接收时刻的卫星位置。事实上,如果直接将求解的卫星位置和接收机所提供的卫星位置对比,会发现所有类型卫星位置三维坐标互差绝对值的RMS均在350 m左右。原因是信号接收时刻的地心地固坐标系,是由在信号发射时刻的坐标系绕地球自转轴(即坐标系的Z轴)自西向东旋转而成的[6-7]。因此,当以诺瓦泰北斗接收机所提供的卫星位置作为真值进行比较时,要结合伪距观测量,进行地球自转误差的修正。

3) 统一时间系统。利用广播星历所计算的北斗卫星位置采用的是北斗时,而精密星历以及北斗接收机所提供的北斗卫星位置采用的却是GPS时,因此对比时需要统一时间系统。本文将时间统一为GPS时。

4) 精密星历按照15 min的时间间隔给出卫星在空间上的三维坐标,采样间隔太大,不利于对比分析。本文采用滑动拉格朗日插值法对所提供的卫星位置进行插值,以获得较小的采样间隔。文献[8]证明了采用滑动式拉格朗日插值法可以获得mm级的精度,满足研究需求。

在解决了上述问题后,便可以分别与精密星历和北斗接收机所提供的卫星位置进行对比分析。将接收机所提供的卫星位置作为真值,三个坐标互差的绝对值分别设为|dX1|、|dY1|、|dZ1|;将精密星历所提供的卫星位置作为真值,三个坐标互差的绝对值分别设为|dX2|、|dY2|、|dZ2|。数据处理结果如图2至图7及表2所示。

星座类型卫星号|dX1|/|dX2|/m|dY1|/|dY2|/m|dZ1|/|dZ2|/m最大值RMS最大值RMS最大值RMS三维坐标 RMS1/m RMS2/mGEOC012.25/3.731.99/3.002.76/7.132.38/6.170.22/0.920.07/0.453.10/6.88C032.70/5.262.38/4.501.12/3.490.90/2.920.14/1.070.05/0.622.55/5.40IGSOC073.12/1.981.50/0.901.64/1.670.70/0.712.64/2.371.64/1.652.33/2.01C090.86/1.770.58/1.600.13/1.980.10/1.860.92/0.680.67/0.450.89/2.50MEOC114.10/1.261.83/0.694.19/3.031.93/1.612.66/1.511.06/0.662.86/1.88C124.00/1.972.12/0.802.34/1.971.11/1.123.18/2.511.63/1.322.89/1.91

对统计结果进行分析可以看出:1) 将北斗接收机所提供的北斗卫星位置作为真值,经过对地球自转误差的修正,三维坐标RMS均在5 m以内。即在位置解算过程中,如果不考虑地球自转的影响,将会产生较大的误差;2) 将精密星历所提供的北斗卫星位置作为真值,三维坐标RMS均在10 m以内;以C07、C09号卫星为代表的IGSO卫星和以C11、C12号卫星为代表的MEO卫星三维RMS小于以C01和C03号卫星为代表的GEO卫星,即IGSO和MEO卫星精度好于GEO卫星。以上结果证明了所解算的卫星位置的正确性,同时验证了星历参数用户算法的正确性。

2 北斗接收机位置解算

本文使用基于线性迭代的最小二乘法解算接收机位置,因为该方法的位置解算仅使用当前数据,因此复杂度低,容易实现,是北斗接收机中最常用的方法[9]。其主要思想是根据当前参数矩阵迭代当前的解,直到误差小于用户预定的门限值,其解算过程如下:

1) 伪距观测量修正

从诺瓦泰接收机中所获得的伪距ρi(i表示不同的卫星)包含了卫星钟差、电离层延时及对流层延时等误差[10],需对其进行修正。即:

ρi=ri+cδtu-cδts+I+T+ερi,

(1)

(2)

2) 构建伪距观测量方程

假设可以得到n颗卫星的伪距观测量,则伪距观测量方程为

(3)

3) 线性化伪距观测方程

假设用户位置和接收机钟差的初始值为(x0,y0,z0,δtu0),基于这个初始值将式(3)进行一阶泰勒级数展开得到:

(4)

定义如下矢量:

hi

(5)

dx0[(xu-x0),(yu-y0),(zu-z0)]T,

(6)

式(5)中前三列由卫星位置与用户估计位置之间的方向余弦组成,第四列为常值1;dx0为线性化后初始点和真实位置之间的修正量。将式(5)、(6)带入到式(4)并略去高阶项整理得:

δρi=hidx0+ερi

(7)

(8)

将式(8)整理为矩阵的形式:

δρ=Hdx0+ερ.

(9)

4) 求解修正量

dx0=(HTH)-1HTδρ.

(10)

5) 更新接收机位置

x1=x0+dx0

(11)

6) 判断牛顿迭代的收敛性

以更新后的接收机位置和钟差为初始位置,重复步骤2)-5),对于第k次更新来说,

(12)

判断‖dxk-1‖是否小于预定的门限值,预定的门限值取决于用户的精度需求,本文将门限值设为0.001。如果已经小于门限值,则可以将当前这一次迭代计算后的更新值(即xk)作为接收机的位置和钟差。即:

xk=xk-1+dxk-1.

(13)

为了验证上述算法的正确性,本文选取天津某已知基准点2018年3月2日2点至23点诺瓦泰北斗接收机接收到的伪距信息结合第2节所解算的北斗卫星位置进行接收机位置解算,并以诺瓦泰接收机所提供的位置信息作为真值进行精度分析。上文提到诺瓦泰FlexPark6接收机水平定位精度达到1.5 m(RMS),时间精度为20 ns,因此可以作为真值来评估所解算接收机位置的精度。数据处理结果如图8至图10及表3所示。

X轴误差/mY轴误差/mZ轴误差/m 最大值1.813.472.15 平均值0.060.170.09 RMS0.100.270.20 三维RMS0.35

由图8和图9可以看出所计算的接收机位置分布较为密集,三维方向上波动不大。图10和表3为对比结果,可以看出: 1) 三个方向上误差最大均不超过5 m,误差平均值均小于0.5 m,三维RMS也小于0.5 m; 2) 通过图10可以看出,三个方向上误差都集中在1 m以下,个别时间段出现较大误差,经分析发现,出现较大误差的原因是卫星个数减少,DOP值变大。结果表明,定位精度均优于北斗官方性能规范中要求的精度[11],证明了所解算接收机位置的正确性,从而证明了最小二乘法用于北斗接收机位置解算的可行性。

3 结束语

本文利用实测数据验证了北斗接收机位置解算中所涉及到的算法的正确性,并指出卫星位置及接收机位置解算过程中需要注意的问题。数据计算结果显示,利用文献[4]中星历参数用户算法解算出的北斗卫星位置以及利用最小二乘法解算出的接收机位置精度均满足北斗官方性能规范要求。接收机位置解算作为接收机基本功能之一,所得出的位置精度直接体现了接收机的性能好坏。因此,对位置解算过程中涉及到的相关算法进行验证,为北斗接收机的设计和研发提供了更可靠的理论依据,具有重要的现实意义。