一种面向民航的北斗卫星轨迹定位拟合方法

罗银辉，郑迦馨

（中国民用航空飞行学院计算机学院，广汉 618307）

0 引言

广播式自动相关监视（automatic dependent surveillance-broadcast，ADS-B）系统是国外的设备，系统的卫星信号源来自美国的GPS，一旦出现GPS信号精度降低甚至不可用的情况，民航飞行、训练保障等方面就会受到很大影响，探索研究将部分航空器的ADS-B机载设备的GPS信号源改成北斗信号源，并在实际运行中检验北斗系统的运行精度、可用性、完好性和连续性。事实证明，基于北斗的ADS-B其导航经度可用性、完好性和连续性并不比GPS差。卫星定位［1］是ADS-B技术的两大核心之一，因此，通过卫星轨迹确定的北斗信号源十分重要。当前，常用的获取轨道位置的方法是通过对精密星历进行多种插值法，包括拉格朗日、牛顿、三角多项式插值方法［2-4］，但是精密星历难以获取，且主要用于事后处理；对于广播星历的插值，沙海等人［5］使用Powell法对卫星轨迹进行拟合，但仅针对2小时之内的数据，时间范围过短；F Cao等人［6］通过构造一个等于脉冲力的恒定力来建立持续增加的加速度来对GEO卫星进行轨道拟合，但是结果表明操纵前的观测结果确实会影响轨道确定和预测，且精度为米级，不够精确。本文针对上述方法的缺点，提出基于最小二乘法［7］原理对三种卫星分别进行长时间的轨迹拟合，并就拟合值再与真实数值进行比较，验证方法的有效性。

1 北斗卫星广播星历及卫星位置计算

1.1 北斗卫星广播星历

北斗卫星系统采用类GPS开普勒根数的广播星历，具有计算简单、外推能力强和高精度的特点，具体参数包含历元时的开普勒根数I、ω、M0、Ω），长期项改正参数（Δn、i0、Ω̇e）以及（Cuc、Crs、Cis、Crc、Cic、Cuc）短周期改正项振幅3种不同类型的参数项［8］。

在GNSS领域中普遍采用与接收机独立交换 格 式（receiver independent exchange format，RINEX）［9］进行数据交换。RINEX格式例子解析如图1所示，第一行“3.03”表示版本号，“N”表示广播星历，第二行表示创建此文件的程序和机构“IONOSPHERIC CORR”表示的是电离层的校正参数，“END OF HEADER”表示表头结束。之后的每八行记录对应卫星在特定历元下的广播星历参数。各个参数所在位置以及所代表的含义如表1所示。

图1 BDS导航信息文件

表1 BDS卫星星历各个参数所在位置以及所代表的的含义

续表1

1.2 卫星在历元t时的位置计算

由于ADS-B主要的位置数据信息来源于北斗卫星导航系统（beidou navigation satellite system，BDS）［10］，因 此 BDS 提 供 的 定 位 信 息 精 确 度 与ADS-B的运行能力息息相关。此时卫星的位置计算就显得尤为重要。在通常情况下，卫星的位置数据是根据接收机接收到的卫星导航电文中的广播星历参数按照直接法由固定的公式来计算，具体步骤如下：

计算卫星运行时的平均角速度：

其中，μ为地球的引力常数。

计算纬度幅角、径向、轨道倾角摄动改正项：

计算IGSO卫星、MEO卫星在ECEF坐标系下的升交点经度：

需要注意的是tk应为两个时间历元之间的实际的全部时间差，必须对星期交接时的起始和结束进行考虑，也就是说如果大于（或小于）302400 s，则要相应的减去（或加上）604800 s。卫星钟差［11］可由下式计算（用表示卫星的编号）：

式中，时间变量的单位是s，计算的钟差单位为10-6s。

2 北斗卫星广播星历拟合

为解决直接法过程繁琐、效率低的问题，本文以时间为自变量，采用最小二乘法对广播星历进行拟合，这样当接收机接收到新的参数时可以直接带入最小二乘法进行计算。

使用x构造具有n+1列和以的长度为行数的Vandermonde矩阵［12］并生成线性方程组。

推导过程：

也就是说X*A=Y，那么A=(X′*X)-1*X′Y，得系数矩阵A。由于Vandermonde矩阵中的列是向量x的幂，因此条件数X对于高阶拟合来说通常较大，生成一个奇异系数矩阵。若处于这些情况下，中心化和缩放可改善系统的数值属性使得拟合更加可靠。

3 仿真实验及结果分析

3.1 实验数据及数据处理

本文利用中国卫星导航系统管理办公室测试评估中心2021年的广播星历（B1I/B3I）数据，选取了从2021年1月1日到2021年7月10日的数据（共172813条数据）。

由于数据包含误差和错误信息，需对数据进行预处理［13］。首先，判断缺失数据和错误数据，并进行清洗。然后，判断格式错误数据，并进行数据校正，形成规格化数据。最后，将数据格式转换为特征和索引格式。流程如图1所示。

图1 数据预处理流程

3.2 轨迹分析

取t为当前历元的GPS时间，对模型进行仿真实验，得到每个时刻的卫星轨迹坐标的曲线图，形成卫星轨迹图，可以看到不同类型卫星的轨迹大不相同。

（1）以2021年C01号卫星为例的GEO卫星在地固坐标系下的轨迹如图2（a）所示。可以看到，GEO卫星的轨迹重叠度很高，基本符合随地球自转而运行的特点。

（2）以2021年C07号卫星为例的IGSO卫星在地固坐标系下的轨迹如图2（b）所示。IGSO卫星运动轨迹垂直向下在地球表面的投影呈现“8”字形，实验中发现卫星轨迹呈螺旋上升又螺旋下降的重复运动，符合星下点［14］轨迹。

（3）以2021年C12号卫星为例的MEO卫星在地固坐标系下的轨迹如图2（c）所示。MEO卫星的运行轨迹垂直向下在地球表面的投影呈现波浪线，实验中发现该卫星运行中仅保证了x、y方向的坐标运行成圆形，而z方向上呈波浪线形，符合星下点轨迹。

图2 2021年卫星运行轨迹

3.3 拟合分析

取12小时内时间间隔为60 min的时间序列t为横坐标，对模型进行仿真实验，得到其对应的实际运行曲线和在时由最小二乘法拟合后的曲线如图3所示，分别用星号、菱形和上三角形表示拟合后的轨迹，用黑色线条表示三条原始曲线，可以看出曲线基本能被拟合，且精度较高。

根据卫星运行轨迹图可以看出对于GEO卫星三个方向上的位置变动较大，而对于MEO和IGSO卫星，相比较于x、y方向的坐标变化，z方向上的变化可以忽略不计。因此在图3（b）～（c）中z方向上的曲线相比于x、y方向几乎呈直线。

图3 卫星轨迹拟合曲线图（续）

图3 卫星轨迹拟合曲线

3.4 误差分析

图4为BDS各类卫星在12小时内经最小二乘法拟合后得到的每个时刻的曲线拟合位置与实际位置之间的误差e，在n+1时，可以看到，GEO卫星、IGSO卫星和MEO卫星的拟合精度分别为10-5m、10-6m和10-4m，以星号、加号和直线表示对于x、y、z坐标的误差值，可以看到误差呈现出中间较稳定两端误差更大的情况，而相比较于x和z坐标的变化，y的误差变化在两端更大，但这对于精度为米级的广播星历来说都可以忽略不计。

图4 2021年卫星坐标误差

图4 2021年卫星坐标误差（续）

4 结语

本文通过北斗卫星广播星历计算出卫星轨迹并可视化，通过使用最小二乘法对12小时内的卫星轨迹进行拟合，拟合精度高，对于以米为精度的广播卫星可以忽略不计。本文提出的拟合方法可以实现轨道的实时运算，并且计算简单，提高了拟合效率，为获取卫星位置信息提供了可靠可行的方法和理论依据，同时也为ADS-B系统的运行提供了获取北斗卫星信号源的方法。