北斗系统地固系自主定轨算法

常家超，尚 琳，李国通,4，赵璐璐1,，肖 洋

(1. 中国科学院上海微系统与信息技术研究所，上海 200050；2. 中国科学院大学，北京 100049;3. 上海微小卫星工程中心，上海 201203；4. 上海科技大学信息科学与技术学院，上海 201210)

0 引 言

基于星间链路的自主定轨技术是新一代北斗导航系统的关键技术之一，可缩短广播星历的更新周期，弱化对地面系统的依赖，提高导航系统生存能力[1-3]。近年来，许多学者对基于星间链路的自主定轨技术开展了研究，如自主定轨体制研究[4-7]，自主定轨滤波算法研究[2,8]，星载轨道预报算法研究[9-10]，星间链路测量值预处理算法研究[11-12]等。目前，北斗系统采用基于锚固站的分布式自主定轨方案，以减少每颗卫星计算量，提高卫星间独立性[1]，并消除单纯基于星间链路的自主定轨算法中存在的星座整体旋转和时间基准漂移问题[1,7,13]。从2015年3月开始，中国已先后发射了5颗新一代北斗导航卫星并建设了几个锚固站[14]，基于这些星间链路节点的在轨自主定轨试验和基于在轨测量数据的自主定轨研究受到广泛关注。文献[1]中，一颗IGSO卫星基于三个锚固站进行自主定轨，9天在轨自主定轨用户测距误差(User range error, URE)小于6 m，事后处理URE小于2.5 m。文献[15]中，使用4颗卫星和1个锚固站获得的星间链路在轨数据，4颗卫星自主定轨的URE均小于5 m (1σ)。上述结果初步证明北斗系统自主定轨方案和算法的可行性。

但是，上述研究基本都是在惯性系中进行自主定轨，即卫星自主估计并预报其在惯性系中的位置速度。由于导航卫星广播星历是定义在地固系中的，需要将预报的惯性系位置转换到地固系中，才可以通过拟合生成广播星历。另外，使用星地测量数据参与自主定轨，需要将锚固站在地固系中的位置转换到惯性系中。因此，在惯性系中进行自主定轨，将会存在多次惯性系与地固系之间的坐标转换，而每次坐标转换过程中均需要计算岁差章动、格林尼治真恒星时和地极极移等参数，进行9次坐标旋转，程序实现复杂，计算量大，且需要上注地球定向参数(Earth Orientation Parameters, EOP)[16]。

为进一步降低自主定轨算法计算量和星地通信数据量，可以直接在地固系中进行自主定轨。文献[13]提出一种基于广播星历参数的地固系集中式自主定轨方法，但该方法与当前北斗系统分布式自主定轨体制不兼容，需要针对地固系分布式动力学自主定轨算法开展研究。其中，状态转移矩阵是基于卡尔曼滤波的自主定轨算法必须计算的重要数据之一。文献[9,17]分别给出了以第二类无奇点根数和惯性系位置速度为状态量的状态转移矩阵解析计算方法，但它们都在惯性系中完成计算，无法直接应用于以地固系位置速度为状态量的自主定轨算法。针对上述问题，本文首先给出了在地固系中进行自主星历生成的总体框图和地固系自主定轨算法的基本方程，然后详细推导以地固系位置速度为状态量的状态转移矩阵解析计算方法，最后基于仿真数据和真实在轨测量数据对提出的地固系自主定轨算法及其中关键技术进行了验证。

1 地固系自主定轨

1.1 自主星历生成总体框图

自主导航模式与传统运行模式拟实现的功能相同，目标是生成尽可能高精度的广播星历发播给用户。图1给出了地固系自主星历生成的总体框图，主要包括自主定轨和广播星历生成两大部分。根据卡尔曼滤波算法的基本思想，自主定轨部分主要包括时间更新模块和测量更新模块。时间更新模块基于地固系长期预报星历进行一步外推完成状态量时间更新，测量更新模块使用预处理之后的星间链路观测量对预测状态量进行修正获得测量更新状态量。在广播星历生成部分，卫星以测量更新状态量为初值，基于地固系长期预报星历进行短弧段轨道预报，然后对预报轨道进行拟合生成广播星历。

相对于在惯性系中进行星历生成，地固系自主星历生成主要存在两点不同：a) 自主定轨部分，滤波状态量为卫星在地固系中的位置速度，星间/星地互传的各类定轨信息也都定义在地固系中。b) 广播星历生成部分，轨道预报模块预报生成卫星地固系坐标序列，可以直接用于广播星历拟合。因此，在地固系中进行星历生成，不存在惯性系与地固系之间的坐标转换，可进一步节省星上计算资源，并且地固系长期预报星历已隐含使用了预报EOP，与锚固站的星地测量数据可有效消除星座整体旋转误差，不需要上注EOP到卫星进行坐标转换，一定程度上节省了星地数传资源。另一方面，由于滤波状态量和长期预报星历都定义在地固系中，基于卡尔曼滤波的自主定轨算法和基于长期预报星历的轨道预报算法都需要在地固系中完成计算，下面将对地固系自主定轨算法进行理论分析。

1.2 地固系自主定轨理论

自主定轨算法基于卫星轨道动力学和星间链路测量值进行轨道参数估计，离散时间状态方程和测量方程可写为[18]：

Xk=f(Xk-1,k-1)+Wk-1

(1)

Zk=h(Xk,k)+Vk

(2)

式中：Xk为状态向量，自主定轨中取为卫星k时刻位置rk、速度vk，即Xk=(rk,vk)；Zk为观测向量，此处为预处理后的距离测量值；f(·)为非线性状态向量函数；h(·)为非线性观测向量函数；Wk-1为系统动态噪声，Vk为测量噪声，Wk-1和Vk为彼此不相关的零均值白噪声序列。

地面精密定轨系统基于精密动力学模型可对卫星轨道进行长时间预报，并保证预报轨道在一定精度范围内。将状态方程和测量方程围绕长期预报星历作泰勒展开，仅保留线性项得：

(3)

(4)

(5)

(6)

(7)

则式(3)和(4)可重写为：

ΔXk=Φk/k-1ΔXk-1+Wk-1

(8)

ΔZk=HkΔXk+Vk

(9)

由式(8)和(9)知，若取ΔXk为状态量，可按照标准卡尔曼滤波过程对其进行估计，完成自主定轨。

(10)

2 地固系状态转移矩阵解析计算方法

由式(8)～(10)知，基于卡尔曼滤波的自主定轨算法和基于长期预报星历的轨道预报算法中，计算状态转移矩阵Φk/k-1都是其中的关键步骤。本节将在惯性系状态转移矩阵解析计算方法的基础上，通过地固系与第二赤道坐标系之间转换，推导以地固系位置速度为状态量的状态转移矩阵解析计算方法。

2.1 坐标系定义

为推导地固系状态转移矩阵计算方法，定义两个第二赤道坐标系OxI,k-1yI,k-1zI,k-1和OxI,kyI,kzI,k。如图2所示，O为地心，赤道平面为长期预报星历所在地固系赤道平面，A点和B点分别为k-1和k时刻格林尼治子午线与赤道平面交点。OxI,k-1yI,k-1zI,k-1和OxI,kyI,kzI,k定义为：OxI,k-1轴和OxI,k轴分别指向A点和B点，OzI,k-1和OzI,k垂直于赤道平面，且OxI,k-1yI,k-1zI,k-1和OxI,kyI,kzI,k构成右手直角坐标系。

OxI,k-1yI,k-1zI,k-1和OxI,kyI,kzI,k分别与k-1和k时刻地固系的坐标轴指向一致，但不随地球自转而旋转，是一种约束到参考时刻的惯性坐标系。在时间间隔为5min时，岁差章动极移引起的赤道面变化和地球自转速度非均匀性引入的旋转误差小于3×10-9rad，几乎可以忽略，这样OzI,k-1轴和OzI,k轴重合为OzI轴，且OxI,k-1yI,k-1zI,k-1和OxI,kyI,kzI,k之间仅存在饶OzI轴的旋转，若取地球自转角速度为ωe，k与k-1时刻之间的时间间隔为Δt，则旋转角度为：

θ=ωeΔt

(11)

2.2 状态转移矩阵计算方法

忽略系统噪声，标明状态量和状态转移矩阵所在坐标系，状态方程(8)可重写为：

ΔXF,k=ΦF,k/k-1ΔXF,k-1

(12)

式中：ΦF,k/k-1表示地固系状态转移矩阵；ΔXF,k-1和ΔXF,k为地固系中真实轨道与长期预报星历在k-1和k时刻的位置速度误差，根据式(5)得：

(13)

(14)

(15)

式中：I为3×3单位矩阵，0为3×3零矩阵；M为地球自转速度修正矩阵。

(16)

(17)

式中：Rz(θ)为绕OzI轴的旋转矩阵。由式(14)～(17)得：

(18)

2) 将ΔXF,k-1、ΔXF,k转换到OxI,kyI,kzI,k

同式(18)，在OxI,kyI,kzI,k中XF,k-1、XF,k可写为：

XI 2,k-1=R2R1XF,k-1,XI 2,k=R1XF,k

(19)

将式(18)和(19)代入式(13)得：

(20)

式中：ΔXI 2,k-1和ΔXI 2,k分别为ΔXF,k-1和ΔXF,k在OxI,kyI,kzI,k中的表达式。

3) 计算状态转移矩阵ΦF,k/k-1

ΔXI 2,k=ΦI 2,k/k-1ΔXI 2,k-1

(21)

将式(20)代入式(21)，得：

ΔXF,k=R3ΦI 2,k/k-1R2R1ΔXF,k-1

(22)

ΦF,k/k-1=R3ΦI 2,k/k-1R2R1

(23)

式(23)即为以地固系位置速度为状态量的状态转移矩阵解析计算方程。

综上，在惯性系状态转移矩阵解析计算方法的基础上，通过将地固系长期预报星历转换到对应时刻约束的第二赤道坐标系，可以完成地固系状态转移矩阵的计算。该方法相对于惯性系状态转移矩阵计算方法仅增加6次6×6维矩阵乘法，且矩阵中包含多个单位矩阵块和零矩阵块，通过程序优化可进一步降低计算量，仅需增加84次浮点乘法运算。而一次坐标系转换中，仅计算岁差章动、格林尼治真恒星时和地极极移等参数就需要超过500次浮点乘法运算，完整的高精度坐标转换程序计算量更大，且程序实现复杂。因此，本文提出的地固系自主定轨算法以计算状态转移矩阵增加少量计算量为代价，避免了多次计算复杂的坐标转换，进一步降低了星载自主定轨算法的计算量。

3 试验校验

3.1 基于仿真数据的算法校验

3.1.1仿真场景

仿真场景选取Walker 24/3/1星座，种子卫星(PRN01)的轨道根数为：a=27875.8 km，e=10-3，i=55°，Ω=15°，ω=0°，M0=101°。三个锚固站分别位于北京、三亚、喀什。卫星精密轨道由STK软件产生，考虑的摄动力包括：GGM02C 70×70阶地球非球形引力、日月三体引力、地球固体潮汐摄动和太阳光压摄动。卫星钟差使用IGS提供的GPS卫星钟差最终产品，考虑到IGS提供钟差产品的连续可用性，使用GPS 01、02、03、05、07、09、10、11、12、13、15、16、18、19、20、22、23、24、25、27、28、29、30、31号卫星从1914周0秒到1923周0秒(GPST)的钟差数据，依次模拟星座中24颗卫星的钟差参数。实际应用中锚固站时钟始终跟踪北斗时(BDT)，所以锚固站钟差设定为固定值。

基于上述模拟数据，按照北斗星间链路时分空分双向测量体制模拟生成伪距测量值。引入的测量误差包括：对流层延迟、相对论效应改正、相位中心改正、系统误差和测量噪声误差。根据目前星间链路测量水平，星间链路测量噪声设置为0.1 m，考虑到星地复杂的空间环境，星地链路测量噪声设置为0.5 m。星座中每颗卫星与异轨卫星建立8条星间链路，包括4条固定可见链路和4条随机可见链路，同时与可见锚固站建立星地测量链路。

相对于精密轨道动力学模型，长期预报星历动力学模型做如下三种改变：(1) 考虑GGM02C 4×4阶地球非球形引力模型；(2) 轨道初值相对精密轨道初值存在(1 m, 1 m, 1 m, 0.0001 m/s, 0.0001 m/s, 0.0001 m/s )偏差。(3) 由于地面精密定轨系统解算的光压参数存在2%～10%的误差[13]，生成两组长期预报星历，分别在光压参数中引入2%和10%误差。

3.1.2仿真结果

基于上述模拟星间/星地测量数据和长期预报星历，每颗卫星独立地解算本星轨道，解算周期为5 min。以URE评估自主定轨算法输出轨道精度。不包含钟差误差的URE定义为[5]：

e(i,k)=

(24)

式中：e(i,k)为卫星i在k时刻的URE；ΔR(i,k)、ΔT(i,k)、ΔN(i,k)为卫星i在k时刻的径向、切向、法向轨道误差。每颗卫星自主定轨n次，则卫星i的URE均方根定义为：

(25)

仿真60天，自主定轨结果如图4和图5所示。星座中各卫星自主定轨结果相似，图4仅给出了PRN01卫星基于含有10%光压误差的长期预报星历自主定轨60天的URE；图5给出了星座中各卫星基于两组不同长期预报星历进行自主定轨的URE均方根值。

由图4可看出，基于3个锚固站自主定轨60天，URE小于2.5 m，且误差没有增大趋势，其波动是由各时刻星地可见性、星座拓扑结构和长期预报星历精度等因素不同引起的，存在约7天的周期变化趋势，与星间链路网络回归周期一致。对比图5中两组基于不同长期预报星历的自主定轨结果，可以看出长期预报星历对自主定轨精度影响较大。基于包含10%光压误差的长期预报星历，卫星自主定轨URE均方根小于0.9 m；基于包含2%光压误差的长期预报星历，卫星自主定轨URE均方根小于0.4 m，两组定轨结果都满足设计要求。将其与文献[9]中惯性系自主定轨结果和文献[13]中基于广播星历的地固系集中式自主定轨结果进行对比。由表1知，本文提出的地固系分布式自主定轨算法与惯性系自主定轨算法和地固系集中式自主定轨算法定轨精度相当，证明了该算法及其中地固系状态转移矩阵计算方法、轨道预报算法等关键技术的可行性。

表1 不同自主定轨算法定轨精度对比Table 1 The orbit determination accuracy comparison with different autonomous orbit determination algorithms

3.2 基于在轨数据的算法校验

3.2.1测试场景

基于新一代北斗导航卫星获得的星间链路实测数据，对本文提出的地固系自主定轨算法进行进一步的验证。考虑获取数据的完整性，使用2016年7月12日到2016年7月18日4颗新一代北斗卫星(Sate1、Sate2、Sate3、Sate4)和北京锚固站(StatBJ)共5个星间链路节点之间的双向测量数据开展自主定轨试验。表2给出了4颗卫星初始时刻的轨道根数，它们在惯性空间的位置关系如图6所示。

北斗星间链路时分测量体制使得双向测距时刻不同，与不同节点之间的测量时刻也不同，且测距值中含有相位中心改正、设备收发时延、对流层延迟、相对论效应改正等误差。因此，自主定轨软件需要

表2 卫星轨道根数Table 2 The orbit elements of satellites

对原始测量数据进行历元归算并消除其中的各类误差[11-12,14]。其中，相位中心改正、设备收发时延和对流层延迟计算方法对于地固系和惯性系自主定轨是相同的，可使用相应的误差模型计算[16]。按照式(14)修正卫星地固系速度后，计算卫星位置矢量与速度矢量内积并乘以相应系数便可得到相对论效应改正。基于预报轨道和钟差计算双向伪距在观测历元相对于目标历元的卫星距离和钟差变化，并将其补偿到伪距值中，便可实现高精度厉元归算[14]。但是，由于观测厉元信号收发时刻不同，预报的收发卫星地固系位置实际上定义在不同的坐标系中。因此，在计算观测厉元收发卫星之间的距离时，需要将不同时刻地固系中的卫星位置统一到相同坐标系中，进而完成伪距厉元归算。

实际运行过程中，北斗主控站使用数值积分方法基于精密动力学模型进行长时间轨道预报获得长期预报星历，并上注卫星作为自主定轨输入参数。此处，使用从主控站获取的长期预报星历进行试验。

3.2.2测试结果

基于实测星间链路数据和主控站生成的长期预报星历，4颗卫星自主定轨原型软件并行估计本星轨道，将输出的定轨结果与精密轨道进行比较评估自主定轨精度。此处，以主控站融合星间链路双向测量数据和监测站L波段伪距相位数据联合定轨结果为精密轨道，其径向误差小于0.1 m，URE小于0.2 m[14,19]。图7为4颗卫星自主定轨7天URE，图8给出了4颗卫星自主定轨URE均方根。

由图7可看到，基于真实星间/星地测量数据和主控站生成的长期预报星历，本文提出的地固系自主定轨算法定轨7天URE小于0.9 m，且误差没有增大的趋势，受星间链路拓扑结构和实测数据量等因素变化影响，定轨误差存在波动，但波动不具有周期性。由图8知，各卫星自主定轨7天的URE均方根在0.26 m～0.40 m之间，与仿真场景中基于包含2%光压误差长期预报星历的自主定轨结果相当，说明主控站生成的长期预报星历精度比较高。上述基于真实数据的定轨结果，进一步证明了所提地固系分布式自主定轨算法及其中关键技术的可行性，同时也验证了星载自主定轨原型软件的有效性。

4 结 论

本文针对地固系分布式自主定轨算法及其中的关键技术进行研究，给出了在地固系中进行自主星历生成的总体框图和地固系自主定轨算法的基本方程，推导了以地固系位置速度为状态量的状态转移矩阵解析计算方法。基于仿真数据和真实在轨测量数据的测试结果表明，所提自主定轨算法与传统算法定轨精度相当，证明该算法及其中的关键技术是可行的，也说明星载自主定轨原型软件是有效的。所提自主定轨算法以计算状态转移矩阵增加少量计算量为代价，避免了多次计算复杂的坐标转换，且地固系长期预报星历已隐含使用了预报EOP，无需另外上注EOP进行坐标转换，可进一步降低星载自主定轨算法计算量和星地通信数据量，具有很强的工程应用价值。

参 考 文 献

[1] Chang J C, Shang L, Li G. The research on system error of inter-satellite-link (ISL) measurements for autonomous navigation of Beidou system [J]. Adv. Space Res., 2017, 60(1): 65-81.

[2] Shang L, Liu G H, Zhang R, et al. An information fusion algorithm for integrated autonomous orbit determination of navigation satellites [J]. Acta Astronautica, 2013, 85(4): 33-40.

[3] 王冬霞, 辛洁, 薛峰, 等. GNSS星间链路自主导航技术研究进展及展望 [J]. 宇航学报, 2016, 37(11): 1279-1288. [Wang Dong-xia, Xin Jie, Xue Feng, et al. Development and prospect of GNSS autonomous navigation based on inter-satellite link [J]. Journal of Astronautics, 2016, 37(11): 1279-1288.]

[4] 宋小勇, 毛悦, 贾小林, 等. 基于星间测距的分布式自主星历更新算法[J]. 武汉大学学报·信息科学版, 2010, 35(10): 1161-1164. [Song Xiao-yong, Mao Yue, Jia Xiao-lin, et al. The distributed processing algorithm for autonomously updating the ephemeris of navigation satellite by inter-satellite links [J]. Geomatics and Information Science of Wuhan University, 2010, 35(10): 1161-1164.]

[5] 林益明, 秦子增, 初海彬, 等. 基于星间链路的分布式导航自主定轨算法研究 [J]. 宇航学报, 2010, 31(9): 2088-2094. [Lin Yi-ming, Qin Zi-zeng, Chu Hai-bin, et al. A satellite cross link-based GNSS distributed autonomous orbit determination algorithm [J]. Journal of Astronautics, 2010, 31(9): 2088-2094.]

[6] 蔡志武, 韩春好, 陈金平,等. 导航卫星长期自主定轨的星座旋转误差分析与控制 [J]. 宇航学报, 2008, 29(2): 522-528. [Cai Zhi-wu, Han Chun-hao, Chen Jin-ping, et al. Constellation rotation error analysis and control in long-term autonomous orbit determination for navigation satellites [J]. Journal of Astronautics, 2008, 29(2): 522-528.]

[7] Yi H, Xu B, Gao Y T, et al. Long-term semi-autonomous orbit determination supported by a few ground stations for navigation constellation [J]. Sci. China Phys. Mech. Astron., 2011, 54(7): 1342-1353.

[8] 尚琳, 刘国华, 刘善伍, 等. 利用分步Kalman滤波器的自主定轨信息融合算法 [J]. 宇航学报, 2013, 34(3): 333-339. [Shang Lin, Liu Guo-hua, Liu Shan-wu, et al. An information fusion algorithm of autonomous orbit determination based-on step Kalman filter [J]. Journal of Astronautics, 2013, 34(3): 333-339.]

[9] Wang H H, Chen Z G, Zheng J J, et al. A new algorithm for onboard autonomous orbit determination of navigation satellites [J]. Journal of Navigation, 2011, 64(S1): S162-S179.

[10] 王海红, 陈忠贵, 初海彬, 等. 导航卫星星载自主轨道预报技术[J]. 宇航学报, 2012, 33(8): 1019-1026. [Wang Hai-hong, Chen Zhong-gui, Chu Hai-bin, et al. On-board autonomous orbit prediction algorithm for navigation satellites [J]. Journal of Astronautics, 2012, 33(8): 1019-1026.]

[11] 毛悦, 宋小勇, 贾小林, 等. 星间链路观测数据归化方法研究[J]. 武汉大学学报·信息科学版, 2013, 38(10): 1201-1206. [Mao Yue, Song Xiao-yong, Jia Xiao-lin, et al. Naturalization method research on inter-satellite link observation data [J]. Geomatics and Information Science of Wuhan University, 2013, 38(10): 1201-1206.]

[12] 宋小勇, 毛悦, 冯来平, 等. BD卫星星间链路定轨结果及分析[J]. 测绘学报, 2017, 46(5): 547-553. [Song Xiao-yong, Mao Yue, Feng Lai-ping, et al. The preliminary result and analysis for BD orbit determination with inter-satellite link data [J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(5): 547-553.]

[13] 毛悦, 胡小工, 宋小勇, 等. 基于广播星历参数的卫星自主导航算法[J]. 中国科学: 物理学 力学 天文学, 2015 45(7): 079512. [Mao Yue, Hu Xiao-gong, Song Xiao-yong, et al. Satellite autonomous navigation algorithm analysis based on broadcast ephemeris parameters [J]. Sci Sin-Phys. Mech. Astron., 2015, 45(7): 079512.]

[14] 唐成盼, 胡小工, 周善石, 等. 利用星间双向测距数据进行北斗卫星集中式自主定轨的初步结果分析 [J]. 中国科学: 物理学 力学 天文学, 2017, 47(2): 029501. [Tang Chen-pan, Hu Xiao-gong, Zhou Shan-shi, et al. Centralized autonomous orbit determination of Beidou navigation satellites with inter-satellite link measurements: preliminary results [J]. Sci Sin-Phys. Mech. Astron., 2017, 47(2): 029501.]

[15] Wang H H, Chen Q L, Jia W S, et al. Research on autonomous orbit determination test based on BDS inter-satellite-link on-orbit data [C]. The 8th China Satellite Navigation Conference, Shanghai, May 24-25, 2017.

[16] 文援兰, 廖瑛, 梁加红, 等. 卫星导航系统分析与仿真技术[M]. 北京: 中国宇航出版社, 2009: 69-74.

[17] Liu L, Zhang Q, Liao X H. Problem of algorithm in precision orbit determination [J]. Science in China (Series A), 1999, 42(5): 552-560.

[18] 秦永元, 张洪钺, 汪叔华. 卡尔曼滤波与组合导航原理[M]. 第3版. 西安: 西北工业大学出版社, 2015: 199-212.

[19] 陈金平, 胡小工, 唐成盼, 等. 北斗新一代试验卫星星钟及轨道精度初步分析 [J]. 中国科学: 物理学 力学 天文学, 2016, 46(11): 119502. [Chen Jin-ping, Hu Xiao-gong, Tang Chen-pan, et al. Orbit determination and time synchronization for new-generation Beidou satellites: Preliminary results [J]. Sci Sin-Phys. Mech. Astron., 2016, 46(11): 119502.]