阅读说明：本技术 一种高精度实时卫星轨道瞬时半长轴的去振算法 (Vibration eliminating algorithm for high-precision real-time satellite orbit instantaneous semi-major axis ) 是由 陈培 程乙洲 贾振俊 耿铖 孙秀聪 于 2021-06-23 设计创作，主要内容包括：本发明提出一种基于离散傅里叶变换和最小二乘拟合的卫星轨道半长轴去振方案,构建适用于卫星轨道相位控制的去振半长轴计算方法,建立了一种结合轨道预报、数值平均和数值拟合的方法,进而设计了一种高精度实时的卫星轨道半长轴的去振算法。本发明采取了数值平均和数值拟合的理念,利用数值平均,给出初次去振半长轴；再利用离散傅里叶变换和数值拟合的方法进一步去振荡。本发明的具体实现方法是利用轨道动力学模型进行轨道预报,从而得到足够长时间段的卫星高精度的绝对位置和速度信息,最后计算得到卫星轨道的瞬时半长轴,利用一定周期内的瞬时半长轴,使用数值平均获取初次去振半长轴；离散傅里叶变换提取残余短周期振荡的固有频率；最后根据残余短周期振荡固有频率,使用数值拟合的方法,计算得到卫星轨道去振平均半长轴。本发明具有导航精度高、实时性好的优点。(The invention provides a satellite orbit semimajor axis vibration eliminating scheme based on discrete Fourier transform and least square fitting, constructs a vibration eliminating semimajor axis calculation method suitable for satellite orbit phase control, establishes a method combining orbit prediction, numerical averaging and numerical fitting, and further designs a high-precision real-time satellite orbit semimajor axis vibration eliminating algorithm. The method adopts the concepts of numerical averaging and numerical fitting, and gives a primary vibration elimination semimajor axis by utilizing the numerical averaging; and further de-oscillating by using a discrete Fourier transform and numerical fitting method. The specific implementation method of the invention is that an orbit dynamics model is utilized to carry out orbit prediction, so as to obtain the high-precision absolute position and speed information of the satellite in a sufficiently long time period, finally the instantaneous semimajor axis of the satellite orbit is obtained by calculation, and the instantaneous semimajor axis in a certain period is utilized to obtain the primary vibration elimination semimajor axis by using numerical value average; extracting the natural frequency of residual short-period oscillation by discrete Fourier transform; and finally, calculating to obtain the satellite orbit vibration-removing average semimajor axis by using a numerical fitting method according to the residual short-period oscillation natural frequency. The invention has the advantages of high navigation precision and good real-time performance.)

一种高精度实时卫星轨道瞬时半长轴的去振算法

技术领域

本发明提供一种高精度实时卫星轨道瞬时半长轴的去振算法，它涉及在一定时间内卫星高精度轨道已知的情况下，通过短时间内的轨道预报加长可利用轨道，轨道瞬时半长轴的分段数值平均初次去振，离散傅里叶变换提取初次去振半长轴残余的短周期振荡固有频率，初次去振半长轴的数值拟合，去掉残余的短周期振荡项，保留长期项和长周期项，最后得到消除短周期振荡后的卫星轨道平均半长轴。属于导航

技术领域

。

背景技术

随着小卫星星座技术的快速发展,小卫星星座受到了世界各航天大国的重视。鸽群卫星(DOVE)星座，由美国Planet Labs公司研制。它使用体积为3U(10cm×10cm×30cm)，重量仅为5.8Kg的CubeSat卫星，组成小卫星编队，可以达到对地3m分辨率的观测，并且可以做到每天遍历拍摄地球一遍，如果使用传统卫星星座造价将非常高昂。星链(StarLink)卫星星座，由著名的SpaceX公司研制，SpaceX公司计划在2024年之前发射超过万颗小卫星，组成卫星星座。根据SpaceX公司，如果一切顺利，星链卫星星座将会提供低成本的互联网连接任务，并且不会受到地理位置的限制，这会极大推动互联网的发展。

但大部分卫星星座都面临着轨道相位保持的问题，一般需要通过轨道相位控制来保持整个星座。假设卫星实际位置为A’，而卫星应该在目标位置为A，它们之间的纬度幅角的夹角为α，卫星需要通过轨道相位控制来保证α角不能太大。

一般而言，卫星是通过控制轨道半长轴来保持轨道相位的。如果一个在理想位置的卫星的轨道半长轴变大，将导致卫星轨道周期增加，这会导致卫星的相位误差越来越大。所以确定卫星轨道半长轴，成为卫星相位控制的重中之重。

卫星的瞬时半长轴，可以分为三部分组成，分别为长期项，长周期项和短周期振荡项。长期项，主要是卫星在地球中心引力场影响下的轨道，是轨道半长轴的主要组成部分，长期项的变化无周期可言，会随着时间逐渐漂移。在地球非球型重力场、地球变形引起的海潮和地固潮、太阳和月球引起的三体引力、大气阻力、太阳光压等摄动力的影响下，轨道半长轴便会产生短周期震荡项和长周期项。卫星轨道半长轴的长周期项是指的是实际轨道半长轴随着时间缓慢地周期性变化，经过一个较长周期后才会比较明显；而卫星轨道半长轴的短周期振荡，是指的轨道半长轴在一个值附近来回快速的振荡，周期相对较短。

目前低轨卫星一般由GNSS定轨技术确定轨道，获得卫星轨道半长轴。GNSS的全称为全球导航卫星系统(Global Navigation Satellite System),包括了美国的GPS、俄罗斯的GLONASS、欧洲的Galileo和中国的北斗等。这是一种能为用户提供空间和时间基准的系统。一种随着GNSS系统的发展，有越来越多的低轨卫星使用了基于GNSS的卫星精密定轨技术。

但是应用GNSS定轨技术，只能获得每个时刻的低轨卫星轨道的瞬时半长轴，这对小卫星相位控制是相当不利的。因为大部分时候，卫星轨道相位控制的目的是为了抑制卫星相位的长期漂移，也就是抑制轨道半长轴的长期漂移，一般不希望对卫星轨道半长轴的短周期振荡做处理，因为短周期振荡项对相位影响较小，且处理起来会浪费很多燃料，这对很多卫星，特别是小卫星而言，是无法接受的。

轨道瞬时半长轴的短周期项，即使对卫星相位的影响较小，但在短期内对瞬时半长轴的影响远大于长期项的漂移部分和长周期项。而卫星的轨道相位控制系统，又恰恰是通过实际半长轴和目标半长轴的差，来确定是否进行轨道控制的。如果不将瞬时半长轴中的短周期振荡项剔除的话，轨道相位控制系统进行控制的次数会远远高于应当控制的次数，这会大大提高燃油损耗，提高卫星的成本。

综上，本发明提出了一种高精度的实时的卫星轨道半长轴的去振算法，它在一定时间段内卫星的高精度轨道已知的情况下，通过短时间的卫星轨道预报获得更长时间段的高精度轨道，之后进行轨道半长轴分段平均(初次去振)，通过离散傅里叶变换确定初次去振后残余的短周期振荡固有频率，最后采取数值拟合的方式，去掉初次去振半长轴中残余的短周期振荡项，保留长期项和长周期项，最后得到消除短周期振荡后的卫星轨道平均半长轴。

发明内容

(一)发明目的

本发明的目的是为了剔除卫星轨道瞬时半长轴的短周期振荡，为轨道相位控制系统提供一个有效的反馈值。本发明提出在一定时间段内卫星的高精度轨道已知的情况下，通过短时间的卫星轨道预报获得更长时间的高精度轨道，之后进行轨道半长轴分段平均(初次去振)，通过离散傅里叶变换确定初次去振后残余的短周期振荡固有频率，最后采取数值拟合的方式，去掉初次去振半长轴中残余的短周期振荡项，保留长期项和长周期项，最后得到消除短周期振荡后的卫星轨道平均半长轴。本发明具有精度高，实时性好的优点，适用于卫星通过轨道半长轴进行轨道相位控制的场景。

(二)技术方案

本发明所述的高精度实时的轨道瞬时半长轴去振算法的实施步骤如下所示：

步骤一：获取一定时间段内的卫星高精度轨道

步骤二：利用轨道动力学预报轨道，获取更长时间段的卫星高精度轨道。

卫星在受到的合外力为F，地球自转角速度为ω，卫星在WGS84坐标系下的位置矢量为r，那么卫星在WGS84坐标系下的动力学方程

其中2ω×r是科里奥利加速度，ω×(ω×r)是离心加速度，m为卫星的质量。F为卫星所受到的合外力，包括地球中心引力、地球非球型重力场、地球变形引起的海潮和地固潮、太阳和月球引起的三体引力、大气阻力、太阳光压等摄动力，其中地球中心引力为最主要的力，其他摄动力可看作小量。

地球的角速度矢量在WGS84坐标系下表示为

ω≈[0 0 7.292×10^-5]^Trad/s (2)

通过方程(1)的数值计算，可以将t_i-1时刻的位置速度r_i-1和v_i-1，预报得到卫星t_i时刻的位置速度r_i和v_i。

步骤三：利用瞬时位置速度，求瞬时半长轴。

在t_i时刻的位置速度大小为

r_i＝||r_i|| (3)

v_i＝||v_i|| (4)

在t_i时刻的瞬时半长轴为

其中a_i为t_i时刻的瞬时半长轴，G为万有引力常量，M为中心天体(地球)质量，r_i为t_i时刻的位置大小，v_i为为t_i时刻的速度大小。

步骤四：检查是否有q个N圈的瞬时半长轴数据，如果没有，重复步骤二到三。

整数值N的取值原则：使卫星绕地球旋转N圈的的时间尽可能的接近1天。

整数值q的取值原则：根据所需要的精度调整，一般取3到7。取值越大本算法获得的去振半长轴精度越高，但所需要的原始数据量和计算量越大。

步骤五：瞬时半长轴的分段平均，获取初次去振半长轴。利用各时刻的瞬时半长轴，计算每分段(绕地球N圈)的数值平均后的初次去振半长轴。

注意，各分段(N圈)之间可以有重合的部分。例如可取第一个分段为第1到第N圈，第二个分段取第2到第N+1圈，第三个分段取第3到第N+2圈，依次类推。这样有利于增加数值平均后的初次去振半长轴的数据量。

假设第m个分段(N圈)中瞬时半长轴的采样点有n_m个，那么第m个分段(N圈)的数值平均后的初次去振半场轴为：

其中是第m个分段(N圈)的初次去振半长轴。

共计算k个分段(N圈)的数值平均后的初次去振半长轴 注意这k个分段可以重合但要完全覆盖步骤四中提取的所有瞬时半长轴数据。

步骤六：离散傅里叶变换提取初次去振半长轴残余短周期震荡的固有频率。

根据步骤五中提取的分段初次去振半长轴进行离散傅里叶变换，提取三个主要的残余短周期振荡频率ω₁、ω₂、ω₃。

步骤七：利用数据拟合进一步降低半长轴振荡

拟合函数：

其中τ为自变量，单位为圈。

其中为拟合后的卫星初次去振半场轴，c+k₁τ+k₂τ²为多项式函数，代表半长轴的长周期项和长期项之和，为所需要的平均项；剩余的三角函数部分为短周期振荡项，需要剔除。ω₁、ω₂、ω₃为瞬时半长轴的振荡周期，需要在拟合前用离散傅里叶变换确定它的具体值。公式(7)中的A₁、A₂、A₃、 k₁、k₂、c为拟合所需要确定的参数。

拟合对象：

拟合的方法：采取最小二乘法，使得公式(7)与拟合对象之间的误差最小。

假设拟合函数为

f(τ)＝A₁g₁(τ,ξ)+A₂g₂(τ,ξ)+…+A_lg_l(τ,ξ) (9)

A＝[A₁ A₂ ... A_l]^T和为需要拟合确定的参数，结合公式(8)和公式(9)，确定参数A和w使得

公式(10)达到最小，minF(A,ξ)。

拟合结果：

在拟合多项式中去除三角函数后，剩下的多项式部分就是去除短周期振荡后的半长轴。

最后去振半长轴结果为

(三)优点

本发明提供的高精度实时卫星轨道的去振算法优点在于：

①利用轨道预报增加了可以利用的轨道长度，提高了实时性。

②本发明提出的在一定周期内，利用瞬时半长轴的平均加拟合的去振算法。该算法消除了瞬时半长轴绝大部分的短周期的振荡，保留了半长轴的长期项和长周期项，为卫星的相位控制提供了有利条件。

附图说明

图1是本发明的实施步骤流程图

具体实施方式

本发明所述的高精度实时的轨道瞬时半长轴去振算法的实施步骤如下所示：

步骤一：获取一定时间段内的卫星高精度轨道

步骤二：利用轨道动力学预报轨道，获取更长时间段的卫星高精度轨道。

卫星在受到的合外力为F，地球自转角速度为ω，卫星在WGS84坐标系下的位置矢量为r，那么卫星在WGS84坐标系下的动力学方程

其中2ω×r是科里奥利加速度，ω×(ω×r)是离心加速度，m为卫星的质量。F为卫星所受到的合外力，包括地球中心引力、地球非球型重力场、地球变形引起的海潮和地固潮、太阳和月球引起的三体引力、大气阻力、太阳光压等摄动力，其中地球中心引力为最主要的力，其他摄动力可看作小量。

地球的角速度矢量在WGS84坐标系下表示为

ω≈[0 0 7.292×10^-5]^Trad/s (13)

通过方程(1)的数值计算，可以将t_i-1时刻的位置速度r_i-1和v_i-1，预报得到卫星t_i时刻的位置速度r_i和v_i。

步骤三：利用瞬时位置速度，求瞬时半长轴。

在t_i时刻的位置速度大小为

r_i＝||r_i|| (14)

v_i＝||v_i|| (15)

在t_i时刻的瞬时半长轴为

其中a_i为t_i时刻的瞬时半长轴，G为万有引力常量，M为中心天体(地球)质量，r_i为t_i时刻的位置大小，v_i为为t_i时刻的速度大小。

步骤四：检查是否有q个N圈的瞬时半长轴数据，如果没有，重复步骤二到三。

整数值N的取值原则：使卫星绕地球旋转N圈的的时间尽可能的接近1天。

整数值q的取值原则：根据所需要的精度调整，一般取3到7。取值越大本算法获得的去振半长轴精度越高，但所需要的原始数据量和计算量越大。

步骤五：瞬时半长轴的分段平均，获取初次去振半长轴。利用各时刻的瞬时半长轴，计算每分段(绕地球N圈)的数值平均后的初次去振半长轴。

注意，各分段(N圈)之间可以有重合的部分。例如可取第一个分段为第1到第N圈，第二个分段取第2到第N+1圈，第三个分段取第3到第N+2圈，依次类推。这样有利于增加数值平均后的初次去振半长轴的数据量。

假设第m个分段(N圈)中瞬时半长轴的采样点有n_m个，那么第m个分段(N圈)的数值平均后的初次去振半场轴为：

其中是第m个分段(N圈)的初次去振半长轴。

共计算k个分段(N圈)的数值平均后的初次去振半长轴 注意这k个分段可以重合但要完全覆盖步骤四中提取的所有瞬时半长轴数据。

步骤六：离散傅里叶变换提取初次去振半长轴残余短周期震荡的固有频率。

根据步骤五中提取的分段初次去振半长轴进行离散傅里叶变换，提取三个主要的残余短周期振荡频率ω₁、ω₂、ω₃。

步骤七：利用数据拟合进一步降低半长轴振荡

拟合函数：

其中τ为自变量，单位为圈。

其中为拟合后的卫星初次去振半场轴，c+k₁τ+k₂τ²为多项式函数，代表半长轴的长周期项和长期项之和，为所需要的平均项；剩余的三角函数部分为短周期振荡项，需要剔除。ω₁、ω₂、ω₃为瞬时半长轴的振荡周期，需要在拟合前用离散傅里叶变换确定它的具体值。公式(7)中的A₁、A₂、A₃、 k₁、k₂、c为拟合所需要确定的参数。

拟合对象：

拟合的方法：采取最小二乘法，使得公式(7)与拟合对象之间的误差最小。

假设拟合函数为

f(τ)＝A₁g₁(τ,ξ)+A₂g₂(τ,ξ)+…+A_lg_l(τ,ξ) (20)

A＝[A₁ A₂ ... A_l]^T和为需要拟合确定的参数，结合公式(8)和公式(9)，确定参数A和w使得

公式(10)达到最小，minF(A,ξ)。

拟合结果：

在拟合多项式中去除三角函数后，剩下的多项式部分就是去除短周期振荡后的半长轴。

最后去振半长轴结果为