一种基于伪相对运动的航天器实时轨道改进方法与流程

本发明属于航天器轨道确定领域，具体涉及一种基于伪相对运动的航天器实时轨道改进方法。
背景技术：
：航天器轨道确定就是对含有测量误差的航天器运动状态数据，使用统计学原理对航天器轨道进行估值的过程。通过轨道确定可以得到航天器在过去、当前和未来一段时间内任一时刻的运动状态。根据是否采用航天器所受摄动力的力学模型及与力学模型的关系，可以将定轨方法分为运动学定轨、动力学定轨和简化动力学定轨；根据数据处理策略可以分为批处理法和序贯递推法；根据弧段的长度可以分为短弧段法和长弧段法；根据积分方法可以分为单步法和多步法。目前，常用的航天器定轨方法是最小二乘估计方法和各种形式的卡尔曼滤波方法。最小二乘估计方法需要在迭代过程中存储大量的数据以供下次迭代使用；卡尔曼滤波方法会由于数据可观测性差、初轨差、测量数据强非线性等造成滤波定轨结果发散；此外，最小二乘估计方法和卡尔曼滤波方法还存在因可观性造成的法方程矩阵和协方差矩阵病态这一数值问题。由于上述传统轨道确定方法的测量精度是有限的，在观测数据精度没有大幅度改进的情况下，无法使得轨道确定精度得到进一步的提高，因此，需要研究新型的轨道改进方法以提高轨道确定的精度。技术实现要素：本发明的目的在于，为解决现有航天器轨道确定方法的测量精度低的技术问题，提供一种基于伪相对运动的航天器实时轨道改进方法。为实现上述目的，本发明提供的一种基于伪相对运动的航天器实时轨道改进方法，包括：步骤1)将航天器的测量轨道作为参考轨道，以初始时刻测量轨道确定的航天器初始状态作为预报轨道的起点，并依据轨道动力学模型生成该预报轨道；步骤2)将预报轨道与测量轨道之间的差值作为轨道漂移数据，并以该轨道漂移数据所形成的轨迹定义为航天器的伪相对运动；步骤3)分别求解无、有摄动条件下伪相对运动的动力学模型，得到无、有摄动条件下伪相对运动的解析表达式，运用傅里叶变换方法处理轨道漂移数据后，以解算得到的解析表达式反演出航天器的初始状态误差的估计值；步骤4)利用初始状态误差的估计值与参考轨道的初始状态误差值之差，作为航天器轨道改进值。作为上述技术方案的进一步改进，所述步骤3)中的无摄动条件下伪相对运动的解析表达式为：其中，n为轨道平均角运动，即参考轨道角频率，t表示时间，各参数及其关系设置为：其中，(x0,y0,z0)表示三轴初始位置误差，表示三轴初始速度误差；运用傅里叶变换方法处理轨道漂移数据，解算得到常值项xcyc、长期项3xcnt/2和周期项csin(nt+φ)中的各系数xc,b,yc,c,φ，并利用上述解析表达式进行反演，得到初始状态误差的估计值作为上述技术方案的进一步改进，所述步骤3)中的有摄动条件下伪相对运动的解析表达式为：x(t)=A1cos(n1t+θ1)+A2cos(n2t+θ2)+A3y(t)=A4cos(n3t+θ3)+A5cos(n4t+θ4)+A6t+A7z(t)=A8sin(n5t+θ5)]]>其中，各参数及其关系设置为：其中，n为轨道平均角运动，即参考轨道角频率，t表示时间，u为地球引力系数，rref为参考轨道地心距，iref为参考轨道倾角，J2为第二带谐项，Re为地球半径，n1、n2、n3、n4、n5分别为相对运动中各周期分量的角频率；运用傅里叶变换方法处理轨道漂移数据，解算得到常值项A3,A7、长期项A6t和周期项A1cos(n1t+θ1),A2cos(n2t+θ2),A4cos(n3t+θ3),A5cos(n4t+θ4),A8sin(n5t+θ5)中的各系数A1,θ1,A2,θ2,A3,A4,θ3,A5,θ4,A6,A7,A8,θ5，并利用上述解析表达式进行反演，得到初始状态误差的估计值本发明的一种基于伪相对运动的航天器实时轨道改进方法的优点在于：1、高定轨精度，可以在初始定轨技术的基础上将定轨精度提高将近一个量级；2、实时性，可以对时间逆向积分，得到当前时刻航天器的位置矢量和速度矢量；3、普遍适应性，利用CW方程及考虑摄动下的改进方程实现圆/近圆航天器轨道改进，也可利用Lawden方程及其考虑摄动下的变形实现椭圆航天器轨道改进。附图说明图1为本发明的一种基于伪相对运动的航天器实时轨道改进方法流程图。图2为本发明实施例中示出的航天器伪相对运动的轨迹示意图。图3为本发明实施例中航天器轨道径向漂移数据及傅里叶变换处理得到的幅频曲线。图4为本发明实施例中航天器轨道迹向漂移数据及傅里叶变换处理得到的幅频曲线。图5为本发明实施例中航天器轨道法向漂移数据及傅里叶变换处理得到的幅频曲线。具体实施方式下面结合附图和实施例对本发明所述的一种基于伪相对运动的航天器实时轨道改进方法进行详细说明。本发明提供的一种基于伪相对运动的航天器实时轨道改进方法，它能够实现初始轨道确定条件下的航天器轨道改进。本发明通过分别求解无/有摄动条件下的航天器相对运动动力学模型，得到无/有摄动条件下相对运动解析表达式，并与轨道漂移数据相对应，运用傅里叶变换方法解算得到航天器初始位置误差和速度误差，最后便可以反演得到航天器实时位置矢量和速度矢量。具体地，一种基于伪相对运动的航天器实时轨道改进方法，其包括如下步骤：①借鉴航天器编队飞行的思想，将航天器的测量轨道(初定轨轨道)作为参考轨道，基于初定轨状态的预报轨道作为对应的“伪卫星”的轨道，由于测量轨道和预报轨道偏差较小，两者之差可以理解为一种相对运动，即伪相对运动；②将初始时刻初定轨确定的航天器初始状态作为预报轨道的起点，依据轨道动力学模型，生成一条预报轨道，将预报轨道与测量轨道做差，即为轨道漂移数据；③分别求解无、有摄动条件下伪相对运动的动力学模型，在不考虑摄动力情况下，当航天器参考轨道为近地圆/近圆轨道时，参考轨道和预报轨道之间的相对运动满足CW方程，即x(t)=(x0·/n)sinnt-(3x0+2y0·/n)cosnt+4x0+2y0·/ny(t)=(6x0+4y0·/n)sinnt+(2x0·/n)cosnt+y0-2x0·/n-(6x0+3y0·/n)ntz(t)=z0cosnt+(z0·/n)sinnt---(1)]]>其中，(x0,y0,z0)表示三轴初始位置误差，表示三轴初始速度误差(在参考轨道的轨道坐标系上)，n为轨道平均角运动，即参考轨道角频率，t表示时间，这些量均为理论已知量；将CW方程转化为常值项、周期项和长期项组合的形式，即其中各参数及其关系设置如下：④参考轨道和预报轨道之间的轨道漂移数据满足上述CW方程，运用傅里叶变换方法处理该轨道漂移数据，解算得到常值项xcyc、长期项3xcnt/2和周期项csin(nt+φ)系数xc,b,yc,c,φ，并利用式(3)进行反演，得到初始状态误差的估计值⑤反演得到的初始状态误差的估计值与初始状态误差之差，即为基于伪相对运动的航天器实时轨道改进精度值；由于所述的预报轨道是时间正向积分得到的，最后得到的当前时刻前的某一时刻的卫星轨道状态；同样地方式，如果将轨道动力学模型进行时间逆向积分，便可得到当前时刻状态为初始值的预报轨道。如果将预报轨道的起点放在当前时刻之前的某一时刻，最后利用上述步骤①-⑤得到的是那一时刻的轨道状态改进精度值；同理，如果将预报轨道的起点放在当前时刻，最后得到的便是当前时刻的轨道状态改进精度值，从而实现了航天器轨道的实时改进。⑥在考虑J2摄动情况下，航天器的伪相对运动动力学模型表示为：x··-2(nc1)y·-(5c12-2)n2x=-3n2J2Rerref(12-3sin2irefsin2(kt)2-1+3cos2iref8)y··+2(nc1)x·=-3n2J2Rerrefsin2irefsin(kt)cos(kt)---(4)]]>其中，分别为三轴运动量的二阶导数，表示动力学过程。其中各参数及其关系设置如下：n=u/rref3c1=1+sk=nc1+3nJ2Re22rref2cos2irefs=3J2Re28rref2(1+3cos2iref)---(5)]]>其中，n为轨道平均角运动，t表示时间，u为地球引力系数，rref为参考轨道地心距，iref为参考轨道倾角，J2为第二带谐项，Re为地球半径。由于J2摄动下航天器的地心距和轨道倾角变化量很小，并考虑到实际工程应用，假设航天器参考轨道为圆轨道，轨道倾角为28.5°、轨道高度为500km时，J2摄动下的地心距的变化范围为4.48km，相对变化量为0.065％；轨道倾角的变化范围为0.03°，相对变化量为0.11％。，为得到J2摄动下相对运动解析形式，合理假设航天器地心距和轨道倾角为常量，具体值分别取为参考轨道平均半长轴和平均轨道倾角，便得到J2摄动下航天器伪相对运动动力学方程：x··+ay·+bx=32ccos2kty··-ax·=csin2kt---(6)]]>其中各参数及其关系设置如下：a=-2ncb=-(5c2-2)n2c=-3n2J2Resin2iref2rref---(7)]]>求解可得伪相对运动解析表达式为：x=x0·sina2+bta2+b+12k(a2+b)(a2-4k2+b)[(8ak3y0·-a3c-2a3ky0·+2a4kx0-8a2k3x0-abc+4ack2-2abky0·+2a2bkx0)+(a3c+abc+3bck+3a2ck)cos2kt+(-3bck-4ack2-3a2ck-8ak3y0·+2a3ky0·+2b2kx0-8bk3x0+2a2bkx0+2abky0·)cosa2+bt]y=1k2(a2+b)(a2-4k2+b)[(a3k2x0·-4bk4y0-4ak4x0·-4a2k4y0+a4k2y0+b2k2y0+abk2x0·+2a2bk2y0)+(b2ck/2-2bck3-4bk4y0·+b2k2y0·+a2bk2y0·-ab2k2x0-a3bk2x0-a2bck/2+4abk4x0)t+(-b2c/4-a2bc/4+3a3ck/4+bck2+a2ck2+3abck/4)sin2kt+(-a3k2x0·+4ak4x0·-abk2x0·)cosa2+bt]+-12(a2+b)3/2(a2-4k2+b)[a(-2a3y0·-2a2bx0+3ca2-2aby0·+8ak2y0·+4ack-2b2x0+8bk2x0+3bc)sin-a2-bt]---(8)]]>需要说明的是，z轴相对运动分量与xy平面相对运动分量解耦、z轴相对运动实现非超越化需要更多的假设条件、且z轴相对运动定轨精度对整体定轨精度影响很小，因此，z轴相对运动直接取为CW方程中z轴运动形式，即x(t)=x0·sina2+bta2+b+12k(a2+b)(a2-4k2+b)[(8ak3y0·-a3c-2a3ky0·+2a4kx0-8a2k3x0-abc+4ack2-2abky0·+2a2bkx0)+(a3c+abc+3bck+3a2ck)cos2kt+(-3bck-4ack2-3a2ck-8ak3y0·+2a3ky0·+2b2kx0-8bk3x0+2a2bkx0+2abky0·)cosa2+bt]y(t)=1k2(a2+b)(a2-4k2+b)[(a3k2x0·-4bk4y0-4ak4x0·-4a2k4y0+a4k2y0+b2k2y0+abk2x0·+2a2bk2y0)+(b2ck/2-2bck3-4bk4y0·+b2k2y0·+a2bk2y0·-ab2k2x0-a3bk2x0-a2bck/2+4abk4x0)t+(-b2c/4-a2bc/4+3a3ck/4+bck2+a2ck2+3abck/4)sin2kt+(-a3k2x0·+4ak4x0·-abk2x0·)cosa2+bt]+-12(a2+b)3/2(a2-4k2+b)[a(-2a3y0·-2a2bx0+3ca2-2aby0·+8ak2y0·+4ack-2b2x0+8bk2x0+3bc)sin-a2-bt]z(t)=c′sin(nt+φ)---(9)]]>其中：c′cosφ=z0·/nc′sinφ=z0---(10)]]>将公式(9)转化为常值项、周期项、长期项组合的形式：x(t)=A1cos(n1t+θ1)+A2cos(n2t+θ2)+A3y(t)=A4cos(n3t+θ3)+A5cos(n4t+θ4)+A6t+A7z(t)=A8sin(n5t+θ5)---(11)]]>其中：其中，n为轨道平均角运动，即参考轨道角频率，t表示时间，u为地球引力系数，rref为参考轨道地心距，iref为参考轨道倾角，J2为第二带谐项，Re为地球半径，n1、n2、n3、n4、n5分别为相对运动中各周期分量的角频率；⑦在J2摄动下参考轨道与预报轨道之间的轨道漂移数据满足上述运动形式，同样运用傅里叶变换方法处理伪相对运动轨道漂移数据，得到常值项A3,A7、长期项A6t和周期项A1cos(n1t+θ1),A2cos(n2t+θ2),A4cos(n3t+θ3),A5cos(n4t+θ4),A8sin(n5t+θ5)中的各系数A1,θ1,A2,θ2,A3,A4,θ3,A5,θ4,A6,A7,A8,θ5，并利用公式(12)进行反向解算，得到初始状态误差的估计值⑧解算得到的初始状态误差的估计值与初始状态误差之差，即为J2摄动下基于伪相对运动的航天器实时轨道改进精度值；另外，在同时考虑J2摄动和大气阻力摄动情况下，近地航天器轨道因大气阻力摄动引起的位置偏差仅为几米至几十米量级，因此在航天器轨道确定中不必再求解大气阻力摄动下的航天器伪相对运动方程，直接使用J2摄动下的航天器伪相对运动方程(公式9)对航天器轨道进行确定便可，轨道改进精度基本不受影响。实施例一如图1所示，记在初始时刻t0，航天器的真实位置位于A点，初定轨(测量轨道)确定航天器位于B点，由于初定轨误差的存在，实际位置A与测量位置B并不重合，记初定轨误差为(在参考轨道的轨道坐标系上，为假设已知量)，此亦为航天器初始状态误差，将此状态误差在时间、空间上进行延伸，即将预报轨道和测量轨道做差，得到一段时间内的轨道漂移数据，亦即为此段时间内的状态漂移量。由于状态漂移量相对状态量为小量，则预报轨道和测量轨道之间的偏差可以看作一种相对运动，且满足一定的规律。(1)情形一：在不考虑外界摄动力的情况下，处理步骤如下：①航天器轨道为近地圆/近圆轨道，则参考轨道和预报轨道之间的相对运动满足CW方程,即x(t)=(x0·/n)sinnt-(3x0+2y0·/n)cosnt+4x0+2y0·/ny(t)=(6x0+4y0·/n)sinnt+(2x0·/n)cosnt+y0-2x0·/n-(6x0+3y0·/n)ntz(t)=z0cosnt+(z0·/n)sinnt]]>②CW运动方程中包含常值项、周期项、长期项，对CW运动方程进行整理可得：其中各参数及其关系设置如下：③对获得的、满足上述CW方程形式的轨道状态漂移数据进行傅里叶变换处理：F{f(t)}=F(w)=∫-∞∞f(t)e-jwtdt]]>F{af1(t)+bf2(t)}＝aF1(w)+bF2(w)得到相应的幅频响应曲线；④根据幅频响应曲线，提取相应的频率、幅值、相位信息，傅里叶变换处理X轴轨道漂移数据，得到xc,b,n，傅里叶变换处理Y轴漂移数据得到yc，处理Z轴，得到c,φ。利用此信息对初始状态偏差进行辨识，即可得到辨识结果即初始状态误差的估计值，两者之差即为轨道改进精度值。(2)情形二：在考虑外界摄动力的情况下，对于近地航天器，其外界摄动力主要是J2摄动和大气阻力摄动。处理步骤如下：①首先仅考虑J2摄动，此种情况下航天器伪相对运动动力学方程表示为：x··-2(nc1)y·-(5c12-2)n2x=-3n2J2Rerref(12-3sin2irefsin2(kt)2-1+3cos2iref8)y··+2(nc1)x·=-3n2J2Rerrefsin2irefsin(kt)cos(kt)]]>其中各参数及其关系设置如下：n=u/rref3c1=1+sk=nc1+3nJ2Re22rref2cos2irefs=3J2Re28rref2(1+3cos2iref)]]>其中，rref为参考星的实时地心距，iref为参考星的实时轨道倾角。②在J2摄动下，上式中航天器地心距和轨道倾角呈现正弦振荡趋势，导致上式微分方程为一超越方程，无法得到解析解形式，但地心距和轨道倾角变化幅值相比其平均值为小量，因此可以在满足精度的前提下合理假设地心距和轨道倾角为恒定值，其值大小由平均轨道根数求得：rref=a‾iref=i‾]]>由此，在J2摄动下航天器伪相对运动动力学方程由超越方程转化为非超越方程。x(t)=x0·sina2+bta2+b+12k(a2+b)(a2-4k2+b)[(8ak3y0·-a3c-2a3ky0·+2a4kx0-8a2k3x0-abc+4ack2-2abky0·+2a2bkx0)+(a3c+abc+3bck+3a2ck)cos2kt+(-3bck-4ack2-3a2ck-8ak3y0·+2a3ky0·+2b2kx0-8bk3x0+2a2bkx0+2abky0·)cosa2+bt]y(t)=1k2(a2+b)(a2-4k2+b)[(a3k2x0·-4bk4y0-4ak4x0·-4a2k4y0+a4k2y0+b2k2y0+abk2x0·+2a2bk2y0)+(b2ck/2-2bck3-4bk4y0·+b2k2y0·+a2bk2y0·-ab2k2x0-a3bk2x0-a2bck/2+4abk4x0)t+(-b2c/4-a2bc/4+3a3ck/4+bck2+a2ck2+3abck/4)sin2kt+(-a3k2x0·+4ak4x0·-abk2x0·)cosa2+bt]+-12(a2+b)3/2(a2-4k2+b)[a(-2a3y0·-2a2bx0+3ca2-2aby0·+8ak2y0·+4ack-2b2x0+8bk2x0+3bc)sin-a2-bt]z(t)=c′sin(nt+φ)]]>此时，在J2摄动条件下伪相对运动的解析表达式为：x(t)=A1cos(n1t+θ1)+A2cos(n2t+θ2)+A3y(t)=A4cos(n3t+θ3)+A5cos(n4t+θ4)+A6t+A7z(t)=A8sin(n5t+θ5)]]>③对满足上述运动方程形式的轨道状态漂移数据进行傅里叶变换处理：F{f(t)}=F(w)=∫-∞∞f(t)e-jwtdt]]>F{af1(t)+bf2(t)}＝aF1(w)+bF2(w)④对初始状态误差进行辨识，得辨识结果即初始状态误差的估计值，两者之差即为航天器轨道改进精度值。另外，在同时考虑J2摄动和大气阻力摄动情况下，对于常规外形航天器，大气阻力摄动只与航天器面质比、轨道高度有较大关系。在500km的轨道高度下，运行两个周期，大气阻力摄动仅会产生10米左右的位置误差，当轨道高度增加时，大气阻力摄动影响量级呈指数递减，并且大气阻力只对航天器迹向相对运动产生影响，对径向、法向相对运动影响很小，并且引入大气模型会增加计算量，因此可以不考虑大气阻力的影响，直接将J2摄动下的伪相对运动方程作为同时考虑J2摄动和大气阻力摄动时的相对运动方程，然后利用傅里叶变换方法对摄动情况下的相对运动数据处理，最终实现航天器的轨道改进。如表1所示，本实施例所选取的航天器的初始状态矢量、初始状态误差及利用本发明的方法反演得到的初始误差的数据分别为：表1仿真初始数据和仿真结果图3即为相应的航天器轨道径向漂移数据及傅里叶变换处理得到的幅频曲线；由幅频响应曲线中的前2个峰值和起始点，可以提取得到三组频率、幅值和相位信息A1,n1,θ1,A2,n2,θ2,A3。图4即为相应的航天器轨道迹向漂移数据及傅里叶变换处理得到的幅频曲线；由幅频响应曲线中的前2个峰值和起始点，可以提取得到三组频率、幅值和相位信息A4,n3,θ3,A5,n4,θ4,A6,A7。图5即为相应的航天器轨道法向漂移数据及傅里叶变换处理得到的幅频曲线。由幅频响应曲线中的第1个峰值处，可以提取得到一组频率、幅值和相位信息A8,n5,θ5。将图3至图5提取到的频率、幅值和相位信息整合在一起，最后利用公式(12)进行反向计算，便可以得到初始状态误差的估计值(即表1中傅里叶变换反演位置误差和傅里叶变换反演速度误差)，此估计值与初始状态误差(即表1中初始航天器位置误差和初始航天器速度误差)之差即为基于伪相对运动的航天器实时轨道改进精度。最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。当前第1页1 2 3