CN114440886B - 一种大偏心率轨道高精度轨道计算方法 - Google Patents

Abstract

本申请公开了一种大偏心率轨道高精度计算方法，包含：S1、生成轨道计算区间内的轨道六根数和位置矢量序列；S2、初始化18维的轨道计算参数，前6个参数取为初始六根数，其余参数取为0；S3、判断轨道计算参数的修正次数是否大于5次，若是进入S4，若否则进入S5；S4、结束轨道计算参数的计算；S5、由轨道计算参数生成轨道六根数序列；S6、将轨道六根数序列转换为位置矢量序列，求取位置矢量计算误差；S7、计算位置矢量关于轨道计算参数的雅可比矩阵；S8、最小二乘法求解轨道计算参数修正量；S9、修正轨道计算参数，返回S3。本发明具有复杂度低、覆盖时间长、计算精度高的特点。

Description

一种大偏心率轨道高精度轨道计算方法

技术领域

本申请涉及航天器导航技术领域，具体涉及一种大偏心率轨道高精度轨道计算方法。

背景技术

大偏心率轨道是一类近地点在1000km左右，远地点高度几万公里的轨道。在大偏心率轨道的远地点附近，卫星运动慢，可见时间长，适合对特殊地区的覆盖，尤其是高纬度地区，如苏联的闪电通讯卫星采用的就是大偏心率轨道。为实现高精度对地指向，卫星对地姿态基准需要实时高可靠地获取高精度的卫星位置信息，精度通常要求小于200m。现有的导航方法包括卫星导航、脉冲星导航、紫外导航和上注参数轨道计算等，其中卫星导航是中低轨卫星的主要导航方式，但对于大偏心率轨道，由于其运行高度跨度大，存在卫星导航信号不稳定、不同弧段精度一致性差的问题。上注参数轨道计算的方式基于地面的测定轨结果进行轨道参数拟合，可靠性高，通常作为备份方案保证导航可靠性或者作为主份应用于导航精度要求高且稳定性强的航天器上。现有的参数轨道计算方法主要由以下两种，一种是通过瞬平转换和平根递推的分析法，另一种是基于瞬根谐波分解的拟合法。这两种方法均无法实现大偏心率轨道的高精度轨道计算，分析法的瞬平转换前提是卫星偏心率较小，大偏心率下瞬平转换将引入病态，转换误差大；拟合法的谐波分解要求轨道在不同弧段的特性均一，谐波特性一致，但大偏心率轨道在不同弧段的轨道角速度不同，特性差异大，谐波拟合精度低。

发明内容

为了解决或部分解决相关技术中存在的问题，本申请提供了一种大偏心率轨道高精度轨道计算方法，采用18参数计算轨道六根数，引入地心距谐波修正项、倾角谐波修正项和纬度幅角谐波修正项分别提升轨道计算的径向、法向和切向精度，通过最小二乘法对18参数进行迭代修正，使得轨道计算精度达到卫星导航精度要求，参数维数小、覆盖时间长，可以有效降低地面注数频次。

本申请第一方面提供了一种大偏心率轨道高精度轨道计算方法，其特征在于，包含：

步骤S1、生成轨道计算时间内的轨道六根数序列和轨道位置矢量序列，作为真值；

步骤S2、初始化轨道计算参数，其维数为18，前6个参数初始化为初始时刻轨道六根数，后12个参数初始化为0；

步骤S3、判断轨道计算参数的修正次数是否大于5次，若是，则进入步骤S4，若否，则进入步骤S5；

步骤S4、结束轨道计算参数的迭代求解；

步骤S5、由轨道计算参数和轨道计算时间生成轨道根数序列；

步骤S6、将计算的轨道根数序列转换为位置矢量序列，与真值序列做差，得到位置矢量计算误差；

步骤S7、由生成的轨道根数序列更新位置矢量关于轨道计算参数的雅可比矩阵；

步骤S8、基于位置矢量计算误差和雅可比矩阵，采用最小二乘法求解轨道计算参数的修正量；

步骤S9、对轨道计算参数进行修正，进入步骤S3，开始新的一次参数修正或结束轨道计算参数的求取。

可选地，所述步骤S2中轨道计算参数参数定义如下：

a₀表示半长轴初值，e₀表示偏心率初值，i₀表示倾角初值，Ω₀表示升交点赤经初值，ω₀表示升交点赤经初值，M₀表示平近点角初值，Δn表示平均轨道角速度修正值，表示轨道倾角一次项，/>表示升交点赤经一次项，C_rs表示地心距修正项正弦系数，C_rc表示地心距修正项余弦系数，C_is表示轨道倾角修正项正弦系数，C_ic表示轨道倾角修正项余弦系数，C_us表示纬度幅角修正项正弦系数1，C_uc表示纬度幅角修正项余弦系数1，C_us1表示纬度幅角修正项正弦系数2，C_uc1表示纬度幅角修正项余弦系数2。

可选地，所述步骤S5中由18个轨道计算参数计算轨道六根数的步骤如下：

计算平均轨道角速度：

修正平均轨道角速度：n＝n₀+Δn；

计算当前时刻相对轨道计算起始时刻的时间增量：Δt＝t_k-t₀；

计算当前时刻平近点角初值：M_k＝M₀+nΔt；

初值取E₀＝M_k，采用Newton-Ralphson法，迭代四次得到偏近点角：

计算真近点角初值：

计算当前时刻的纬度幅角初值：u_k0＝f_k0+ω₀；

计算纬度幅角修正项：Δu＝C_uccos 2u_k0+C_ussin 2u_k0+C_uc1cos u_k0+C_us1sin u_k0；

计算地心距修正项：Δr＝C_rccos 2u_k0+C_rssin 2u_k0；

计算轨道倾角修正项：Δi＝C_iccos 2u_k0+C_issin 2u_k0；

修正纬度幅角：u_k＝u_k0+Δu；

计算地心距：

计算轨道倾角：

计算升交点赤经：

修正真近点角：f_k＝u_k-ω₀；

计算半长轴：

更新偏近点角：

计算平近点角：M_k＝E_k-e₀sin E_k；

得到当前时刻轨道根数：σ_k＝[a_k e₀ i_k Ω_k ω₀ M_k]。

可选地，所述步骤S7中位置矢量关于轨道计算参数的雅可比矩阵计算步骤如下：

计算平均轨道角速度：

计算半通径：p＝a₀(1-e₀)；

计算偏近点角，初值取E₀＝M₀，进行四次Newton-Ralphson法迭代：

计算真近点角：

计算纬度幅角：u＝f+ω₀；

计算地心距：

计算位置矢量关于半长轴初值的偏导数：

计算位置矢量关于偏心率初值的偏导数：

计算位置矢量关于轨道倾角初值的偏导数：

计算位置矢量关于升交点赤经初值的偏导数：

计算位置矢量关于近地点幅角初值的偏导数：

计算位置矢量关于平近点角初值的偏导数：

计算位置矢量关于平均轨道角速度修正值的偏导数：

计算位置矢量关于轨道倾角一次项的偏导数：

计算位置矢量关于升交点赤经一次项的偏导数：

计算位置矢量关于地心距修正项正弦系数的偏导数：

计算位置矢量关于地心距修正项余弦系数的偏导数：

计算位置矢量关于轨道倾角修正项正弦系数的偏导数：

计算位置矢量关于轨道倾角修正项余弦系数的偏导数：

计算位置矢量关于纬度幅角修正项正弦系数1的偏导数：

计算位置矢量关于纬度幅角修正项余弦系数1的偏导数：

计算位置矢量关于纬度幅角修正项正弦系数2的偏导数：

计算位置矢量关于纬度幅角修正项余弦系数2的偏导数：

由偏导数组成当前时刻位置矢量关于轨道参数的雅可比矩阵：

合并所有离散时刻的雅可比矩阵：

H＝[H_1x … H_kx H_1y … H_ky H_1z … H_kz]^T。

可选地，所述步骤S8中，轨道计算参数修正量的计算步骤如下：

对位置矢量关于轨道计算参数的雅可比矩阵进行奇异值分解：

[U,S,V]＝svd(H)；

最小二乘法求解轨道计算参数修正量：

其中，[Δr_x Δr_y Δr_z]^T为轨道计算的位置矢量误差序列。

本申请提供的技术方案可以包括以下有益效果：

本发明具有复杂度低、注数次数少、覆盖时间长等情况的同时，保证轨道计算高精度的特点，为大偏心率轨道高精度导航提供一种有效简便的方法。

应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本申请。

附图说明

为了更清楚地说明本申请专利实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请专利的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本申请实施例中

具体实施方式

下面将参照附图更详细地描述本申请的实施方式。虽然附图中显示了本申请的实施方式，但是应该理解的是，可以以各种形式实现本申请而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本申请更加透彻和完整，并且能够将本申请的范围完整地传达给本领域的技术人员。

在本申请使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本申请。在本申请和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。还应当理解，本文中使用的术语“和/或”是指并包含一个或多个相关联的列出项目的任何或所有可能组合。

应当理解，尽管在本申请可能采用术语“第一”、“第二”、“第三”等来描述各种信息，但这些信息不应限于这些术语。这些术语仅用来将同一类型的信息彼此区分开。例如，在不脱离本申请范围的情况下，第一信息也可以被称为第二信息，类似地，第二信息也可以被称为第一信息。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本申请的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。

下文将结合附图对本申请实施例的技术方案进行详细描述。

请参阅图1，本实施例提供了一种大偏心率轨道高精度轨道计算方法，其特征在于，包含：

步骤S1、生成轨道计算时间内的初始时刻轨道六根数和轨道位置矢量序列，作为真值；

步骤S2、初始化轨道计算参数，其维数为18，前6个参数初始化为初始时刻轨道六根数，后12个参数初始化为0；

步骤S3、判断轨道计算参数的修正次数是否大于5次，若是，则进入步骤S4，结束轨道计算参数的迭代求解，若否，则进入步骤S5；

步骤S5、由轨道计算参数和轨道计算时间计算轨道根数序列；

步骤S6、将计算的轨道根数序列转换为位置矢量序列，与真值做差，得到位置矢量计算误差；

步骤S7、由计算的轨道根数序列更新位置矢量关于轨道计算参数的雅可比矩阵；

步骤S8、基于位置矢量计算误差和雅可比矩阵，采用最小二乘法求解轨道计算参数的修正量；

步骤S9、对轨道计算参数进行修正，进入步骤S3，开始新的一次参数修正或结束轨道计算参数的求取。

可选地，所述步骤S2中轨道计算参数参数定义如下：

可选地，所述步骤S5中由18个轨道计算参数计算轨道六根数的步骤如下：

计算平均轨道角速度：

修正平均轨道角速度：n＝n₀+Δn；

计算当前时刻相对轨道计算起始时刻的时间增量：Δt＝t_k-t₀；

计算当前时刻平近点角初值：M_k＝M₀+nΔt；

初值取E₀＝M_k，采用Newton-Ralphson法，迭代四次得到偏近点角：

计算真近点角初值：

计算当前时刻的纬度幅角初值：u_k0＝f_k0+ω₀；

计算纬度幅角修正项：Δu＝C_uccos 2u_k0+C_ussin 2u_k0+C_uc1cos u_k0+C_us1sin u_k0；

计算地心距修正项：Δr＝C_rccos 2u_k0+C_rssin 2u_k0；

计算轨道倾角修正项：Δi＝C_iccos 2u_k0+C_issin 2u_k0；

修正纬度幅角：u_k＝u_k0+Δu；

计算地心距：

计算轨道倾角：

计算升交点赤经：

修正真近点角：f_k＝u_k-ω₀；

计算半长轴：

更新偏近点角：

计算平近点角：M_k＝E_k-e₀sin E_k；

得到当前时刻轨道根数：σ_k＝[a_k e₀ i_k Ω_k ω₀ M_k]。

可选地，所述步骤S7中位置矢量关于轨道计算参数的雅可比矩阵计算步骤如下：

计算平均轨道角速度：

计算半通径：p＝a₀(1-e₀)；

计算偏近点角，初值取E₀＝M₀，进行四次Newton-Ralphson法迭代：

计算真近点角：

计算纬度幅角：u＝f+ω₀；

计算地心距：

计算位置矢量关于半长轴初值的偏导数：

计算位置矢量关于偏心率初值的偏导数：

计算位置矢量关于轨道倾角初值的偏导数：

计算位置矢量关于升交点赤经初值的偏导数：

计算位置矢量关于近地点幅角初值的偏导数：

计算位置矢量关于平近点角初值的偏导数：

计算位置矢量关于平均轨道角速度修正值的偏导数：

计算位置矢量关于轨道倾角一次项的偏导数：

计算位置矢量关于升交点赤经一次项的偏导数：

计算位置矢量关于地心距修正项正弦系数的偏导数：

计算位置矢量关于地心距修正项余弦系数的偏导数：

计算位置矢量关于轨道倾角修正项正弦系数的偏导数：

计算位置矢量关于轨道倾角修正项余弦系数的偏导数：

计算位置矢量关于纬度幅角修正项正弦系数1的偏导数：

计算位置矢量关于纬度幅角修正项余弦系数1的偏导数：

计算位置矢量关于纬度幅角修正项正弦系数2的偏导数：

计算位置矢量关于纬度幅角修正项余弦系数2的偏导数：

由偏导数组成当前时刻位置矢量关于轨道参数的雅可比矩阵：

合并所有离散时刻的雅可比矩阵：

H＝[H_1x … H_kx H_1y … H_ky H_1z … H_kz]^T。

可选地，所述步骤S8中，轨道计算参数修正量的计算步骤如下：

对位置矢量关于轨道计算参数的雅可比矩阵进行奇异值分解：

[U,S,V]＝svd(H)；

最小二乘法求解轨道计算参数修正量：

其中，[Δr_x Δr_y Δr_z]^T为轨道计算的位置矢量误差序列。

以上所述，仅为本申请的实施例而已，并非用于限定本申请的保护范围。凡在本申请的精神和范围之内做出的任何修改、等同替换和改进等，均包含在本申请的保护范围之内。

Claims (3)

1.一种大偏心率轨道高精度轨道计算方法，其特征在于，包含：

步骤S1、生成轨道计算时间内的轨道六根数序列和位置矢量序列，作为真值；

步骤S2、初始化轨道计算参数，其维数为18，前6个参数初始化为初始时刻轨道六根数，后12个参数初始化为0，

所述轨道计算参数参数定义如下：

a₀表示半长轴初值，e₀表示偏心率初值，i₀表示倾角初值，Ω₀表示升交点赤经初值，ω₀表示升交点赤经初值，M₀表示平近点角初值，Δn表示平均轨道角速度修正值，表示轨道倾角一次项，/>表示升交点赤经一次项，C_rs表示地心距修正项正弦系数，C_rc表示地心距修正项余弦系数，C_is表示轨道倾角修正项正弦系数，C_ic表示轨道倾角修正项余弦系数，C_us表示纬度幅角修正项正弦系数1，C_uc表示纬度幅角修正项余弦系数1，C_us1表示纬度幅角修正项正弦系数2，C_uc1表示纬度幅角修正项余弦系数2；

步骤S3、判断轨道计算参数的修正次数是否大于5次，若是，则进入步骤S4，若否，则进入步骤S5；

步骤S4、结束轨道计算参数求解；

步骤S5、由轨道计算参数和轨道计算时间序列生成相应的轨道根数序列，

由18个轨道计算参数和轨道计算时间生成轨道六根数的步骤如下：

计算平均轨道角速度：

修正平均轨道角速度：n＝n₀+Δn；

计算当前时刻相对轨道计算起始时刻的时间增量：Δt＝t_k-t₀；

计算当前时刻平近点角初值：M_k＝M₀+nΔt；

初值取E₀＝M_k，采用Newton-Ralphson法，迭代四次得到偏近点角：

计算真近点角初值：

计算当前时刻的纬度幅角初值：u_k0＝f_k0+ω₀；

计算纬度幅角修正项：Δu＝C_uc cos2u_k0+C_us sin2u_k0+C_uc1 cosu_k0+C_us1 sinu_k0；

计算地心距修正项：Δr＝C_rc cos2u_k0+C_rs sin2u_k0；

计算轨道倾角修正项：Δi＝C_ic cos2u_k0+C_is sin2u_k0；

修正纬度幅角：u_k＝u_k0+Δu；

计算地心距：

计算轨道倾角：i_k＝i₀+Δi+iΔt；

计算升交点赤经：

修正真近点角：f_k＝u_k-ω₀；

计算半长轴：

更新偏近点角：

计算平近点角：M_k＝E_k-e₀sin E_k；

得到当前时刻轨道根数：σ_k＝[a_k e₀ i_k Ω_k ω₀ M_k]；

步骤S6、将生成的轨道根数序列转换为位置矢量序列，与真值序列做差，得到位置矢量计算误差；

步骤S7、基于生成的轨道根数序列，更新位置矢量关于轨道计算参数的雅可比矩阵；

步骤S8、基于位置矢量计算误差和雅可比矩阵，采用最小二乘法求解轨道计算参数的修正量；

步骤S9、对轨道计算参数进行修正，进入步骤S3，开始新的一次参数修正或结束轨道计算参数的求取。

2.如权利要求1所述的大偏心率轨道高精度轨道计算方法，其特征在于，所述步骤S7中位置矢量关于轨道计算参数的雅可比矩阵计算步骤如下：

计算平均轨道角速度：

计算半通径：p＝a₀(1-e₀)；

计算偏近点角，初值取E₀＝M，进行四次Newton-Ralphson法迭代：

计算真近点角：

计算纬度幅角：u＝f+ω₀；

计算地心距：

计算位置矢量关于半长轴初值的偏导数：

计算位置矢量关于偏心率初值的偏导数：

计算位置矢量关于轨道倾角初值的偏导数：

计算位置矢量关于升交点赤经初值的偏导数：

计算位置矢量关于近地点幅角初值的偏导数：

计算位置矢量关于平近点角初值的偏导数：

计算位置矢量关于平均轨道角速度修正值的偏导数：计算位置矢量关于轨道倾角一次项的偏导数：/>

计算位置矢量关于升交点赤经一次项的偏导数：

计算位置矢量关于地心距修正项正弦系数的偏导数：

计算位置矢量关于地心距修正项余弦系数的偏导数：

计算位置矢量关于轨道倾角修正项正弦系数的偏导数：

计算位置矢量关于轨道倾角修正项余弦系数的偏导数：

计算位置矢量关于纬度幅角修正项正弦系数1的偏导数：

计算位置矢量关于纬度幅角修正项余弦系数1的偏导数：

计算位置矢量关于纬度幅角修正项正弦系数2的偏导数：

计算位置矢量关于纬度幅角修正项余弦系数2的偏导数：

由偏导数组成当前时刻位置矢量关于轨道参数的雅可比矩阵：

合并所有离散时刻的雅可比矩阵：

H＝[H_1x … H_kx H_1y … H_ky H_1z … H_kz]^T。

3.如权利要求1所述的大偏心率轨道高精度轨道计算方法，其特征在于，所述步骤S8中，轨道计算参数修正量的计算步骤如下：对位置矢量关于轨道计算参数的雅可比矩阵进行奇异值分解：

[U,S,V]＝svd(H)；

最小二乘法求解轨道计算参数修正量：

其中，[Δr_x Δr_y Δr_z]^T为轨道计算的位置矢量误差序列。