CN115792982A - 北斗导航卫星广播星历参数拟合方法及存储介质 - Google Patents

Abstract

本发明的一种北斗导航卫星广播星历参数拟合方法及存储介质，包括首先用户获取北斗的广播星历；其次将根据星历参考时刻的位置和速度求出所需参数，包括星历参考时刻的轨道半长轴相对于参考值的偏差

平近点角

偏心率

近地点幅角

倾角

以及周历元零时刻升交点经度

在此之上设定除参考时间t_oe以外的17个星历参数的求解范围，待求解的轨道参数记为向量形式p，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组；最后使用Levenberg‑Marquardt算法对方程组迭代求解，求出待求星历参数向量p的最优解。本发明的算法复杂度低，能够以更高效率求解出最优解，能够以较高精度拟合北斗卫星广播星历，实现了很好的拟合效果。

Description

北斗导航卫星广播星历参数拟合方法及存储介质

技术领域

本发明涉及卫星导航技术领域，具体涉及一种北斗导航卫星广播星历参数拟合方法。

背景技术

卫星星历可分为精密星历和广播星历两种，其中精密星历提供精度更高的信息，用户获取精密星历主要用于后处理阶段。而广播星历则在用户实时导航定位中有至关重要的作用，用户接收广播星历利用其包含的参数形式的卫星轨道信息，解算出卫星的具体位置。但广播星历的发布时间间隔较长，为了保证动态实时定位性能，用户需要利用精密测轨数据，进行卫星广播星历参数的拟合。现有星历拟合方法有拉格朗日插值法、牛顿插值法、切比雪夫拟合法等，但大多都是针对GPS系统的星历拟合研究，对北斗系统的星历拟合研究相对较少。

北斗卫星导航系统(以下简称北斗)是我国全面自主研制的卫星系统，包括北斗一号、北斗二号和北斗三号系统。与GPS系统构型不同，北斗卫星系统采用混合构型，包括中圆地球轨道(MEO)、地球静止轨道(GEO)、倾斜地球同步轨道(IGSO)三种轨道类型的卫星。北斗系统现已能面向全球提供服务，具备定位导航授时、精密单点定位、短报文通信和国际搜救等多种服务能力。如今全球各国都在开发和完善自己的卫星导航定位系统，开发针对北斗广播星历参数的拟合算法有很大的研究意义。

根据中国卫星导航系统管理办公室发布的北斗卫星导航系统ICD(接口控制文档)，北斗系统播发的B1I和B3I信号采用16参数的广播星历，而B1C、B2a和B2b信号的广播星历包含19个参数，包括1个卫星轨道类型参数和18个准开普勒轨道参数，采用了与GPS相似的广播星历模型。未来，北斗将逐渐从北斗二号过渡到北斗三号系统。而B1C、B2a和B2b信号通过北斗三号中圆地球轨道(MEO)卫星和倾斜地球同步轨道(IGSO)卫星播发，地球静止轨道(GEO)卫星不播发这三个信号。北斗播发的18参数广播星历包括：星历参考时刻t_oe，长半轴的偏差ΔA，长半轴的变化率

平均角速度偏差Δn₀，平均角速度偏差变化率

平近点角M₀，偏心率e、近地点幅角ω，升交点经度Ω₀，轨道面倾角i₀、升交点经度变化率的偏差

轨道面倾角的变化率i₀，轨道倾角的正弦调和改正项振幅C_is和余弦调和改正项的振幅C_ic，轨道半径的正弦调和改正项的振幅C_rs和余弦调和改正项的振幅C_rc，纬度幅角的正弦调和改正项的振幅C_us和余弦调和改正项的振幅C_uc。相对于广播星历16参数，18参数中长半轴、平均角速度和升交点赤经的计算方法不一样，更加体现卫星瞬时变化性。现在大多研究发明围绕16参数的广播星历模型，因此需要针对北斗广播星历参数研究新的拟合算法，也能为北斗系统导航定位研究提供参考。

发明内容

本发明提出的一种北斗导航卫星广播星历参数拟合方法，可至少解决上述技术问题之一。

为实现上述目的，本发明采用了以下技术方案：

一种北斗导航卫星广播星历参数拟合方法，其实施步骤如下：

步骤一：准备包含星历参考时刻的一段时间内的卫星位置和速度及对应的时间信息；

步骤二：根据星历参考时刻的位置和速度求出星历参考时刻的轨道半长轴相对于参考值的偏差

平近点角

偏心率

近地点幅角

倾角

以及周历元零时刻升交点经度

步骤三：设定17个星历参数的求解范围，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组；

步骤四：设定带求解参数初始值，使用Levenberg-Marquardt算法对方程组迭代求解。

进一步的，在步骤一中所述的“准备包含星历参考时刻在内的一段时间内的卫星位置和速度及对应的时间信息”，其做法如下：

准备包含星历参考时刻在内的一段时间内的卫星位置和速度及对应的时间信息，星历参考时刻记为t_oe，第i个时刻的卫星位置和速度记为x_i＝(r_i,v_i)，其中，r_i＝(r_x,i,r_y,i,r_z,i)，v_i＝(v_x,i,v_y,i,v_z,i)，为了能够求解出17个轨道参数x_i的数据点数应不小于6；

进一步的，在步骤二中所述的“根据星历参考时刻的位置和速度求出星历参考时刻的轨道半长轴相对于参考值的偏差

平近点角

偏心率

近地点幅角

倾角

以及周历元零时刻升交点经度

”，其做法如下：

S21、将星历参数的求解简化为二体问题，则可以根据星历参考时刻的位置和速度求出星历参考时刻的轨道半长轴相对于参考值的偏差

平近点角

偏心率

近地点幅角

倾角

升交点经度

S22、求周历元零时刻的格林尼治恒星时；

周历元零时刻的儒略日整数部分记为jd，单位为天，小数部分对应的时间记为ut，单位为秒，周历元零时刻格林尼治恒星时GMS计算方法为：

GMS＝24110.54841+8640184.812866·T+0.093104·T²-6.2×10^-6·T³+1.002737909350795·ut

S23、计算周历元零时刻的升交点经度；

式中，mod(GMS,86400)为GMS对86400的余数。

进一步的，在步骤三中所述的“设定17个星历参数的求解范围，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组”，其做法如下：

待求解的轨道参数记为向量形式:

p中元素含义依次为：参考时刻长半轴相对于参考值的偏差、长半轴变化率、参考时刻卫星平均角速度与计算值之差、参考时刻卫星平均角速度与计算值之差的变化率、参考时刻的平近点角、偏心率、近地点幅角、周历元零时刻计算的升交点经度、参考时刻的轨道倾角、升交点赤经变化率、轨道倾角变化率、轨道倾角的正弦调和改正项的振幅、轨道半径的正弦调和改正项的振幅、轨道半径的正弦调和改正项的振幅、轨道半径的余弦调和改正项的振幅、纬度幅角的正弦调和改正项的振幅、纬度幅角的余弦调和改正项的振幅；

时刻t_k，由p到地心地固坐标系下卫星位置r_k的函数关系记为：r_k＝f(p,t_k)，

由x求p的最优化问题：

min||F(p)-b||

s.t. c_i≤p_i≤d_i i＝1,2,…,17

可转换为相应的拉格朗日对偶问题：

s.t. c_i≤p_i≤d_i i＝1,2,…,17

式中，b＝[r₁,r₂,...,r_n]，b为1×3n的矩阵，n为数据点数，F(p)＝[f(p,t₁),f(p,t₂),...,f(p,t_n)]，F(p)为1×3n的矩阵，p_i为p中第i个元素，λ_i为拉格朗日乘子，c_i和d_i为p_i的上下界。

其中，在步骤四中所述的“设定带求解参数初始值，使用Levenberg-Marquardt算法对方程组迭代求解”，其做法如下：

S41、记

取初始点p₀，取λ_i＝10^-6,i＝1,2,...,17，终止控制常数ε，计算

k＝0,m₀＝10^-3，v＝10，其中，

S42、计算雅克比矩阵J_k，计算

构建增量正规方程

S43、求解增量正规方程得到δ_k；

S44、如果

则令p_k+1＝p_k+δ_k，若||δ_k||<ε_k，停止迭代，输出结果；否则令m_k+1＝m_k/v，转到步骤S42。

S45、如果

则令m_k+1＝v·m_k，重新解正规方程得到δ_k，返回步骤S42。

又一方面，本发明还公开一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如上述方法的步骤。

由上述技术方案可知，本发明的北斗导航卫星广播星历参数拟合方法，包括四个步骤，首先用户获取北斗的广播星历，其中包含星历参考时刻在内的一段时间内的卫星位置和速度信息；其次将根据星历参考时刻的位置和速度求出所需参数，包括星历参考时刻的轨道半长轴相对于参考值的偏差

平近点角

偏心率

近地点幅角

倾角

以及周历元零时刻升交点经度

在此之上设定除参考时间t_oe以外的17个星历参数的求解范围，待求解的轨道参数记为向量形式p，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组；最后使用Levenberg-Marquardt算法对方程组迭代求解，求出待求星历参数向量p的最优解。

通过上述步骤，本发明的一种北斗导航卫星广播星历参数拟合方法通过算法优化，算法复杂度低，能够以更高效率求解出最优解，能够以较高精度拟合北斗卫星广播星历，实现了很好的拟合效果。

依据本发明的设计，本发明实现了一种北斗导航卫星广播星历参数拟合方法，能够很好的拟合北斗卫星广播星历。

依据本发明的设计，本发明实现了一种北斗导航卫星广播星历参数拟合方法，在针对北斗系统广播星历包含19个参数的情况下，满足了实时导航定位的精度要求，对北斗系统的发展和研究进程有重要的参考价值。

附图说明

图1是本发明的星历解算流程；

图2中为使用两小时内240条数据对某颗卫星进行星历拟合后的误差。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。

首先概述本发明实施例的核心原理，B1C、B2a和B2b信号的广播星历包含19个参数，包括1个卫星轨道类型参数和18个准开普勒轨道参数，相对于广播星历16参数，18参数中长半轴、平均角速度和升交点赤经的计算方法不一样，更加体现卫星瞬时变化性。如图1所示，是本次发明算法的流程逻辑框架图，具体实施步骤如下：

第一步：准备包含星历参考时刻在内的一段时间内的卫星位置和速度及对应的时间信息，星历参考时刻记为t_oe，第i个时刻的卫星位置和速度记为x_i＝(r_i,v_i)，其中r_i＝(r_x,i,r_y,i,r_z,i)，v_i＝(v_x,i,v_y,i,v_z,i)，为了能够求解出17个轨道参数x_i的数据个数应不小于6；

第二步：将星历参数的求解简化为二体问题，则可以根据星历参考时刻的位置和速度求出星历参考时刻的轨道半长轴相对于参考值的偏差

平近点角

偏心率

近地点幅角

倾角

升交点经度

计算方法如下：

表1

周历元零时刻的儒略日整数部分记为jd，单位为天，小数部分对应的时间记为ut，单位为秒，周历元零时刻格林尼治恒星时GMS计算方法为：

GMS＝24110.54841+8640184.812866·T+0.093104·T²-6.2×10^-6·T³+1.002737909350795·ut

计算周历元零时刻的升交点经度：

式中，mod(GMS,86400)为GMS对86400的余数。

第三步：设定17个星历参数的求解范围，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组；

其中将待求解的轨道参数记为向量形式

p中元素含义依次为：参考时刻长半轴相对于参考值的偏差、长半轴变化率、参考时刻卫星平均角速度与计算值之差、参考时刻卫星平均角速度与计算值之差的变化率、参考时刻的平近点角、偏心率、近地点幅角、周历元零时刻计算的升交点经度、参考时刻的轨道倾角、升交点赤经变化率、轨道倾角变化率、轨道倾角的正弦调和改正项的振幅、轨道半径的正弦调和改正项的振幅、轨道半径的正弦调和改正项的振幅、轨道半径的余弦调和改正项的振幅、纬度幅角的正弦调和改正项的振幅、纬度幅角的余弦调和改正项的振幅。

由p到r的计算过程如下表所示：

表2

在时刻t_k，由p到地心地固坐标系下卫星位置r_k的函数关系记为：r_k＝f(p,t_k)，由星历时刻的卫星位置和速度向量x＝(r,v)求p的最优化问题：

min||F(p)-b||

s.t. c_i≤p_i≤d_i i＝1,2,…,17

可转换为相应的拉格朗日对偶问题：

s.t. c_i≤p_i≤d_i i＝1,2,…,17

式中，b＝[r₁,r₂,...,r_n]，R为1×3n的矩阵，n为数据点数，F(p)＝[f(p,t₁),f(p,t₂),...,f(p,t_n)]，F(p)为1×3n的矩阵，p_i为p中第i个元素，λ_i为拉格朗日乘子，c_i和d_i为p_i的上下界。上下界的选取如下表所示：

表3

第四步：使用Levenberg-Marquardt算法对方程组迭代求解，其做法如下：

S41、记

取初始点p₀，取λ_i＝10^-6,i＝1,2,...,17，终止控制常数ε，计算

k＝0,m₀＝10^-3，v＝10，其中，

S42、计算雅克比矩阵J_k，计算

构建增量正规方程

其中，雅克比矩阵J_k使用数值法求得，J_k为17×n列的矩阵；J_k(j,i)代表J_k的第j行，第i列；

上式，p_k为p第k次迭代的结果，

为p_k中第i个元素的值在原来的基础上加上Δp，其余元素的值和p_k中保持一致。

S43、求解增量正规方程得到δ_k；

S44、如果

则令p_k+1＝p_k+δ_k，若||δ_k||<ε_k，停止迭代，输出结果；否则令m_k+1＝m_k/v，转到步骤S42。

S45、如果

则令m_k+1＝v·m_k，重新解正规方程得到δ_k，返回步骤S42。

星历拟合误差可用下式评估：

式中，f(p,t_i)为使用求解出的星历计算出的第一步中准备的数据对应时刻的卫星位置，r_i为第1步中准备的卫星位置数据。

若err大于用户预期值，则可以将求解到的轨道参数设为迭代的初始值，重新进行第4步的迭代。

图2中为使用两小时内240条数据对某颗卫星进行星历拟合后的误差，ECEF坐标系下X、Y、Z三个方向的星历拟合误差均小于1.2m，可以满足导航和实时定位的要求。

综上所述，本发明实施例的一种北斗导航卫星广播星历参数拟合方法通过算法优化，算法复杂度低，能够以更高效率求解出最优解，能够以较高精度拟合北斗卫星广播星历，实现了很好的拟合效果。

又一方面，本发明还公开一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如上述任一方法的步骤。

再一方面，本发明还公开一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行如上述任一方法的步骤。

在本申请提供的又一实施例中，还提供了一种包含指令的计算机程序产品，当其在计算机上运行时，使得计算机执行上述实施例中任一方法的步骤。

可理解的是，本发明实施例提供的系统与本发明实施例提供的方法相对应，相关内容的解释、举例和有益效果可以参考上述方法中的相应部分。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一非易失性计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (6)

1.一种北斗导航卫星广播星历参数拟合方法，该方法步骤如下：

步骤一：准备包含星历参考时刻的一段时间内的卫星位置和速度及对应的时间信息；

步骤二：根据星历参考时刻的位置和速度求出星历参考时刻的轨道半长轴相对于参考值的偏差

平近点角

偏心率

近地点幅角

倾角

以及周历元零时刻升交点经度

步骤三：设定17个星历参数的求解范围，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组；

步骤四：设定带求解参数初始值，使用Levenberg-Marquardt算法对方程组迭代求解，完成对卫星广播星历参数的拟合。

2.根据权利要求1所述的一种北斗导航卫星广播星历参数拟合方法，其特征在于：所述步骤一具体过程如下：

准备包含星历参考时刻的一段时间内的卫星位置和速度及对应的时间信息，星历参考时刻记为t_oe，第i个时刻的卫星位置和速度记为x_i＝(r_i,v_i)，其中，r_i＝(r_x,i,r_y,i,r_z,i)，v_i＝(v_x,i,v_y,i,v_z,i)，为了能够求解出17个轨道参数x_i的数据点数应不小于6。

3.根据权利要求1所述的一种北斗导航卫星广播星历参数拟合方法，其特征在于：所述步骤二具体过程如下：

S21、将星历参数的求解简化为二体问题，则根据星历参考时刻的位置和速度求出星历参考时刻的轨道半长轴相对于参考值的偏差

平近点角

偏心率

近地点幅角

倾角

升交点经度

S22、求周历元零时刻的格林尼治恒星时；

周历元零时刻的儒略日整数部分记为jd，单位为天，小数部分对应的时间记为ut，单位转为秒，周历元零时刻格林尼治恒星时GMS计算方法为：

GMS＝24110.54841+8640184.812866·T+0.093104·T²-6.2×10^-6·T³+1.002737909350795·ut

S23、计算周历元零时刻的升交点经度；

式中，mod(GMS,86400)为GMS对86400的余数。

4.根据权利要求1所述的一种北斗导航卫星广播星历参数拟合方法，其特征在于：所述步骤三具体过程如下：

设定17个星历参数的求解范围，将求解星历参数的最优化问题转为等价的拉格朗日对偶问题，构建方程组；

待求解的轨道参数记为向量形式:

p中元素含义依次为：参考时刻长半轴相对于参考值的偏差、长半轴变化率、参考时刻卫星平均角速度与计算值之差、参考时刻卫星平均角速度与计算值之差的变化率、参考时刻的平近点角、偏心率、近地点幅角、周历元零时刻计算的升交点经度、参考时刻的轨道倾角、升交点赤经变化率、轨道倾角变化率、轨道倾角的正弦调和改正项的振幅、轨道半径的正弦调和改正项的振幅、轨道半径的正弦调和改正项的振幅、轨道半径的余弦调和改正项的振幅、纬度幅角的正弦调和改正项的振幅、纬度幅角的余弦调和改正项的振幅；

时刻t_k，由p到地心地固坐标系下卫星位置r_k的函数关系记为：r_k＝f(p,t_k)，

由x求p的最优化问题：

min||F(p)-b||

s.t.c_i≤p_i≤d_i i＝1,2,…,17

可转换为相应的拉格朗日对偶问题：

s.t.c_i≤p_i≤d_i i＝1,2,…,17

式中，b＝[r₁,r₂,...,r_n]，b为1×3n的矩阵，n为数据点数，F(p)＝[f(p,t₁),f(p,t₂),...,f(p,t_n)]，F(p)为1×3n的矩阵，p_i为p中第i个元素，λ_i为拉格朗日乘子，c_i和d_i为p_i的上下界。

5.根据权利要求1所述的一种北斗导航卫星广播星历参数拟合方法，其特征在于：所述步骤四具体过程如下：

S41、记

取初始点p₀，取λ_i＝10^-6,i＝1,2,...,17，终止控制常数ε，计算

k＝0,m₀＝10^-3，v＝10；，其中，

S42、计算雅克比矩阵J_k，计算

构建增量正规方程

S43、求解增量正规方程得到δ_k；

S44、如果

则令p_k+1＝p_k+δ_k，若||δ_k||<ε_k，停止迭代，输出结果；否则令m_k+1＝m_k/v，转到步骤S42；

S45、如果

则令m_k+1＝v·m_k，重新解正规方程得到δ_k，返回步骤S42。

6.一种计算机可读存储介质，存储有计算机程序，所述计算机程序被处理器执行时，使得所述处理器执行如权利要求1至5中任一项所述方法的步骤。

Images