一种基于最小二乘的非线性相对轨迹预报方法与流程
本发明涉及航天器相对运动轨迹预报技术，具体涉及一种基于最小二乘的非线性相对轨迹预报方法。
背景技术：
航天器相对运动轨迹预报在航天器编队、集群等飞行控制任务中具有重要应用。当两编队航天器飞离测控区域时，需要对其相对运动轨迹进行预报，以判断是否需要进行构型保持、碰撞规避等操作。
用于描述航天器相对运动的方程主要有基于近圆轨道的C-W方程和基于椭圆轨道的T-H方程。但这两个线性方程仅适用于二体假设下近距离的相对轨迹预报，作为改进，不少学者相继研究了考虑二阶非线性项、J2摄动项等的非线性相对运动方程，可以用于更长时间、更高精度的相对轨迹预报。但这些解析的相对运动方程仍然只是对实际飞行力学环境的一定近似，即存在模型误差。
再者，已有研究表明，相对轨迹预报精度不仅与采用的方程有关，还受初始条件(如参考航天器位置、初始相对状态等)影响。实际任务中，一般会有很多历史数据(即预报的初始条件)可以利用。现有相对运动方程都是基于单个数据点进行相对轨迹预报，且没有哪一个方程能对不同的初始点都具有优于其他方程的预报精度。
因此，若直接采用现有相对运动方程进行相对轨迹预报，一方面历史数据没有得到充分利用；另一方面也不知道应该选择哪个初始点和那个相对运动方程，才能获得最高的预报精度。本发明基于最小二乘思想，利用历史数据估计出新的预报初始条件，可以有效解决以上两个问题，最小化相对运动方程中的模型误差与初始条件中的测量误差传播造成的预报误差，以获得最高的相对轨迹预报精度。
技术实现要素：
本发明要解决的技术问题：针对现有技术的上述问题，提供一种考虑了J2摄动项和二阶非线性项，可用于相距较远航天器的长时间、高精度相对轨迹预报，具有设计方法正确合理、对实际工程任务的适用性好的基于最小二乘的非线性相对轨迹预报方法。
为了解决上述技术问题，本发明采用的技术方案为：
一种基于最小二乘的非线性相对轨迹预报方法，步骤包括：
1)将已有的n个历史历元时刻对应的、在主航天器当地轨道坐标系中表示的主航天器和从航天器的相对运动状态历史轨道数据写入观察向量；
2)利用相对运动方程计算各历史历元时刻的预报相对运动状态向量，计算各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量；
3)采用最小二乘方法迭代求解使残差向量的残差和最小的最优初始相对运动状态
4)将所述最优初始相对运动状态带入所述相对运动方程，向前预报到指定的历元时刻tf，获得指定的历元时刻tf的相对运动轨迹预报结果。
优选地，所述步骤1)的详细步骤包括：
1.1)分别指定编队飞行的两航天器为主航天器与从航天器，进行相对轨迹预报的初始时刻为t0，主航天器初始轨道要素为E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]，其中，E为主航天器初始轨道，t0为相对轨迹预报的初始时刻，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近拱点角距，f为真近点角；令从航天器相对主航天器的初始相对运动状态为x(t0)，其中x(t0)在主航天器当地轨道坐标系中表示，该坐标系原点为主航天器质心，x轴沿主航天器地心矢径方向，z轴沿轨道面法向，y轴与x、z轴构成右手坐标系；
1.2)获取n个历史历元时刻从航天器相对主航天器的相对运动状态，记为：[t-n,x(t-n)],[t-(n-1),x(t-(n-1))],…,[t-1,x(t-1)]，其中x(t-n)表示时刻t-n从航天器相对主航天器的初始相对运动状态，x(t-(n-1))表示时刻t-(n-1)从航天器相对主航天器的初始相对运动状态，x(t-1)表示时刻t-1从航天器相对主航天器的初始相对运动状态；
1.3)将n个历史历元时刻从航天器相对主航天器的初始相对运动状态历史轨道数据写入观察向量Z＝[x(t-n),x(t-(n-1)),...,x(t-1)]，其中，x(t-n)表示时刻t-n从航天器相对主航天器的相对运动状态，x(t-(n-1))表示时刻t-(n-1)从航天器相对主航天器的相对运动状态，x(t-1)表示时刻t-1从航天器相对主航天器的相对运动状态。
优选地，所述步骤2)的详细步骤包括：
2.1)利用相对运动方程计算各历史历元时刻的预报相对运动状态向量，得到各历史历元时刻[t-n,t-(n-1),...,t-1]的预报相对运动状态向量其中，表示时刻t-n从航天器相对主航天器的预报相对运动状态，表示时刻t-(n-1)从航天器相对主航天器的预报相对运动状态，表示时刻t-1从航天器相对主航天器的预报相对运动状态；
2.2)根据式(2.2-1)计算各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量；
式(2.2-1)中，e表示各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量，Z表示观察向量，Y表示各历史历元时刻的预报相对运动状态向量，x(t-n)表示时刻t-n从航天器相对主航天器的相对运动状态，x(t-(n-1))表示时刻t-(n-1)从航天器相对主航天器的相对运动状态，x(t-1)表示时刻t-1从航天器相对主航天器的相对运动状态；表示时刻t-n从航天器相对主航天器的预报相对运动状态，表示时刻t-(n-1)从航天器相对主航天器的预报相对运动状态，表示时刻t-1从航天器相对主航天器的预报相对运动状态。
优选地，所述步骤2.1)中利用相对运动方程计算各历史历元时刻的预报相对运动状态向量时，如果仅已知主航天器的轨道半长轴a(t0)或平均轨道角速度其中μ为地心引力常数，a为半长轴，则利用式(2.1-1)所示C-W方程计算各历史历元时刻[t-n,t-(n-1),...,t-1]的预报相对运动状态向量；如果已知主航天器的6个轨道根数E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]，则利用式(2.1-2)所示考虑J2摄动的非线性相对运动方程计算各历史历元时刻[t-n,t-(n-1),...,t-1]的预报相对运动状态向量；
式(2.1-1)中，表示时刻t从航天器相对主航天器的预报相对运动状态，Φ(t,t0)是从t0时刻到t时刻的状态转移矩阵，x0表示t0时刻两航天器的初始相对运动状态，τ＝n(t-t0)，s＝sinτ，c＝cosτ，为主航天器的平均轨道角速度，t0为相对轨迹预报的初始时刻，a为半长轴，μ为地心引力常数；
式(2.1-2)中，表示时刻t从航天器相对主航天器的预报相对运动状态，Φ(t,t0)是从t0时刻到t时刻的状态转移矩阵，为从θ0纬度幅角到θ纬度幅角的状态转移矩阵，x0表示t0时刻两航天器的初始相对运动状态，Ψ(t,t0)为从t0时刻到t时刻的二阶状态转移张量，为从θ0纬度幅角到θ纬度幅角的二阶状态转移张量，T(t)为以纬度幅角θ为自变量的无量纲化坐标到以时间t为自变量的量纲化坐标的转换矩阵，T(t0)为以纬度幅角θ为自变量的无量纲化坐标到以时间t0为自变量的量纲化坐标的转换矩阵，T-1(t0)为转换矩阵T(t0)的逆，为Kronecker张量积。
优选地，所述步骤3)的详细步骤包括：
3.1)定义最小二乘方法的目标函数为式(3.1-1)所示的残差和，使得式(3.1-1)所示残差和最小的条件为式(3.1-2)；
式(3.1-1)中，J表示各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量的残差和，Y表示各历史历元时刻的预报相对运动状态向量，Z表示观察向量；
式(3.1-2)中，J表示各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量的残差和，x0表示t0时刻两航天器的初始相对运动状态，Y表示各历史历元时刻的预报相对运动状态向量，e表示各历史历元时刻的预报结果和历史轨道数据中历史数据的残差向量；
3.2)将相对运动方程、各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量带入式(3.1-2)，取一阶近似可得初始相对运动状态的最优估计值满足式(3.2-1)所示的迭代条件进行迭代求解；
式(3.2-1)中，k为迭代次数，表示第k+1次迭代的结果，表示第k次迭代的结果，的迭代初值x(t0)为时刻t0从航天器相对主航天器的初始相对运动状态，Fk向量的表达式为Fk＝[(Φk)TΦk]-1(Φk)T，Φk表示第k次迭代中的一阶状态转移矩阵，由式(2.1-1)及(2.1-2)可知，Φk在每次迭代中为常量；
3.3)判断是否满足指定的迭代条件，如果不满足指定迭代条件，则跳转执行步骤3.2)；否则当满足指定迭代条件时，将最后一次迭代得到的第k+1次迭代的结果作为最终使得主航天器和从航天器之间相对轨迹预报残差的残差和最小的最优初始相对运动状态
优选地，所述步骤3.3)中指定的迭代条件具体是指第k+1次迭代的结果与第k次迭代的结果之间的变化量小于或等于给定迭代误差、或者迭代次数k达到给定的最大值K。
优选地，所述步骤4)获得指定的历元时刻tf的相对运动轨迹预报结果如式(4-1)或式(4-2)所示；
式(4-1)中，x(tf)表示指定的历元时刻tf的相对运动轨迹预报结果，Φ(tf,t0)是从t0时刻到tf时刻的状态转移矩阵，表示步骤3.3)中获得的最优初始相对运动状态；
式(4-2)中，x(tf)表示指定的历元时刻tf的相对运动轨迹预报结果，Φ(tf,t0)表示是从t0时刻到tf时刻的状态转移矩阵，Ψ(tf,t0)表示从t0时刻到tf时刻的二阶状态转移张量，表示步骤3.3)中获得的最优初始相对运动状态。
本发明基于最小二乘的非线性相对轨迹预报方法的优点如下：
1、本发明基于最小二乘的非线性相对轨迹预报方法利用相对运动方程计算各历史历元时刻的预报相对运动状态向量，计算各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量，考虑了航天器实际飞行环境中的主要摄动因素J2项及二阶非线性项，可用于两相距较远航天器的长时间、高精度相对轨迹预报。
2、本发明利用相对运动方程计算各历史历元时刻的预报相对运动状态向量，计算各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量，采用最小二乘方法迭代求解使残差向量的残差和最小的最优初始相对运动状态将最优初始相对运动状态带入相对运动方程，向前预报到指定的历元时刻tf，获得指定的历元时刻tf的相对运动轨迹预报结果通过解析的相对运动方程进行相对轨迹预报，计算速度快。
附图说明
图1为本发明实施例一方法的基本流程示意图。
图2为本发明实施例方法的相对轨迹预报位置误差对比示意图。
图3为本发明实施例方法的相对轨迹预报速度误差对比示意图。
具体实施方式
实施例一：
如图1所示，本实施例基于最小二乘的非线性相对轨迹预报方法的步骤包括：
1)准备历史轨道数据：
将已有的n个历史历元时刻对应的、在主航天器当地轨道坐标系中表示的主航天器和从航天器的相对运动状态历史轨道数据写入观察向量；
2)计算各历史历元时刻预报结果和历史轨道数据之间的残差：
利用相对运动方程计算各历史历元时刻的预报相对运动状态向量，计算各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量；
3)采用最小二乘方法迭代求解最优初始相对运动状态：
采用最小二乘方法迭代求解使残差向量的残差和最小的最优初始相对运动状态
4)带入最优初始相对运动状态进行相对运动轨迹预报：
将最优初始相对运动状态带入相对运动方程，向前预报到指定的历元时刻tf，获得指定的历元时刻tf的相对运动轨迹预报结果。
本实施例中，步骤1)的详细步骤包括：
1.1)分别指定编队飞行的两航天器为主航天器与从航天器，进行相对轨迹预报的初始时刻为t0，主航天器初始轨道要素为E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]，其中，E为主航天器初始轨道，t0为相对轨迹预报的初始时刻，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近拱点角距，f为真近点角；令从航天器相对主航天器的初始相对运动状态为x(t0)，其中x(t0)在主航天器当地轨道坐标系中表示，该坐标系原点为主航天器质心，x轴沿主航天器地心矢径方向，z轴沿轨道面法向，y轴与x、z轴构成右手坐标系；
1.2)获取n个历史历元时刻从航天器相对主航天器的相对运动状态，记为：[t-n,x(t-n)],[t-(n-1),x(t-(n-1))],…,[t-1,x(t-1)]，其中x(t-n)表示时刻t-n从航天器相对主航天器的初始相对运动状态，x(t-(n-1))表示时刻t-(n-1)从航天器相对主航天器的初始相对运动状态，x(t-1)表示时刻t-1从航天器相对主航天器的初始相对运动状态；
1.3)将n个历史历元时刻从航天器相对主航天器的初始相对运动状态历史轨道数据写入观察向量Z＝[x(t-n),x(t-(n-1)),...,x(t-1)]，其中，x(t-n)表示时刻t-n从航天器相对主航天器的相对运动状态，x(t-(n-1))表示时刻t-(n-1)从航天器相对主航天器的相对运动状态，x(t-1)表示时刻t-1从航天器相对主航天器的相对运动状态。本实施例中，观察向量具体为6n×1的观测列向量，通过前述步骤1.3)最终将已有的n个历史历元时刻对应的、在主航天器当地轨道坐标系中表示的两航天器相对运动状态写入一个6n×1的观测列向量。
本实施例中，步骤2)的详细步骤包括：
2.1)利用相对运动方程计算各历史历元时刻的预报相对运动状态向量，得到各历史历元时刻[t-n,t-(n-1),...,t-1]的预报相对运动状态向量其中，表示时刻t-n从航天器相对主航天器的预报相对运动状态，表示时刻t-(n-1)从航天器相对主航天器的预报相对运动状态，表示时刻t-1从航天器相对主航天器的预报相对运动状态；
2.2)根据式(2.2-1)计算各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量；
式(2.2-1)中，e表示各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量，Z表示观察向量，Y表示各历史历元时刻的预报相对运动状态向量，x(t-n)表示时刻t-n从航天器相对主航天器的相对运动状态，x(t-(n-1))表示时刻t-(n-1)从航天器相对主航天器的相对运动状态，x(t-1)表示时刻t-1从航天器相对主航天器的相对运动状态；表示时刻t-n从航天器相对主航天器的预报相对运动状态，表示时刻t-(n-1)从航天器相对主航天器的预报相对运动状态，表示时刻t-1从航天器相对主航天器的预报相对运动状态。
本实施例中，步骤2.1)中利用相对运动方程计算各历史历元时刻的预报相对运动状态向量时，如果仅已知主航天器的轨道半长轴a(t0)或平均轨道角速度其中μ为地心引力常数，a为半长轴，则利用式(2.1-1)所示C-W方程计算各历史历元时刻[t-n,t-(n-1),...,t-1]的预报相对运动状态向量；如果已知主航天器的6个轨道根数E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]，则利用式(2.1-2)所示考虑J2摄动的非线性相对运动方程计算各历史历元时刻[t-n,t-(n-1),...,t-1]的预报相对运动状态向量；
式(2.1-1)中，表示时刻t从航天器相对主航天器的预报相对运动状态，Φ(t,t0)是从t0时刻到t时刻的状态转移矩阵，x0表示t0时刻两航天器的初始相对运动状态，τ＝n(t-t0)，s＝sinτ，c＝cosτ，为主航天器的平均轨道角速度，t0为相对轨迹预报的初始时刻，a为半长轴，μ为地心引力常数；
式(2.1-2)中，表示时刻t从航天器相对主航天器的预报相对运动状态，Φ(t,t0)是从t0时刻到t时刻的状态转移矩阵，为从θ0纬度幅角到θ纬度幅角的状态转移矩阵，x0表示t0时刻两航天器的初始相对运动状态，Ψ(t,t0)为从t0时刻到t时刻的二阶状态转移张量，为从θ0纬度幅角到θ纬度幅角的二阶状态转移张量，T(t)为以纬度幅角θ为自变量的无量纲化坐标到以时间t为自变量的量纲化坐标的转换矩阵，T(t0)为以纬度幅角θ0为自变量的无量纲化坐标到以时间t0为自变量的量纲化坐标的转换矩阵，反之量纲化坐标到无量纲化坐标的转换矩阵T-1(t0)为转换矩阵T(t0)的逆，为Kronecker张量积。例如，表示用矩阵X的每一个元素与矩阵Y相乘，所得结果为mp×nq的矩阵Z。本实施例中利用相对运动方程将初始相对状态x0向后分别预报到每个历史数据对应的历元时刻，即[t-n,t-(n-1),...,t-1]，获得预报的相对运动状态向量：若仅知道主航天器的轨道半长轴a(t0)(或平均轨道角速度利用C-W方程(2.1-1)预报；若主航天器的6个轨道根数E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]均已知，则利用考虑J2摄动的非线性相对运动方程(2.1-2)预报。
式(2.1-2)可表示为如式(2.1-3)所示分量形式；
式(2.1-3)中，表示t时刻预报相对运动状态的第i个分量，Φi,a(t,t0)表示t0时刻到t时刻的一阶状态转移矩阵的第i行、第a列元素，Ψi,ab(t,t0)表示t0时刻到t时刻的二阶状态转移张量的第i模块、第a行、第b列元素，xa(t0)表示初始t0时刻相对运动状态的第a个分量，xb(t0)表示初始t0时刻相对运动状态的第b个分量，Φij(t,t0)表示t0时刻到t时刻的一阶状态转移矩阵的第i行、第j列元素，Til(t)表示以纬度幅角θ为自变量的无量纲化坐标到以时间t为自变量的量纲化坐标的转换矩阵的第i行、第l列元素，表示从θ0纬度幅角到θ纬度幅角的一阶状态转移矩阵的第l行、第m列元素，Tmj(t0)表示以纬度幅角为自变量的无量纲化坐标到以时间为自变量的量纲化坐标的转换矩阵的第m行、第j列元素，Ψijk(t,t0)表示t0时刻到t时刻的二阶状态转移张量的第i模块、第j行、第k列元素，表示从θ0纬度幅角到θ纬度幅角的二阶状态转移张量的第l模块、第m行、第n列元素，Tnk(t0)表示t0时刻以纬度幅角为自变量的无量纲化坐标到以时间为自变量的量纲化坐标的转换矩阵的第n行、第k列元素。式(2.1-3)采用爱因斯坦求和助记，省略对哑变量的求和符号，即下标a表示状态x(t0)的第a个分量，N＝6为状态变量的维数。哑变量a,b,i,j,k,m,n,l取值均为集合{1,2,…,N}。
状态转移矩阵与二阶状态转移张量可用式(2.1-4)求解。
式(2.1-4)中，表示从θ0纬度幅角到θ纬度幅角的一阶状态转移矩阵的第i行、第j列元素，表示状态转移矩阵的第i行、第l列元素，表示状态转移矩阵的第l行、第j列元素，表示状态转移张量的第i模块、第l行、第m列元素，表示状态转移矩阵的第m行、第k列元素，表示状态转移张量的第i模块、第l行、第m列元素，表示状态转移矩阵的第m行、第j列元素，表示状态转移矩阵的第n行、第k列元素，哑变量a,b,i,j,k,m,n,l取值均为集合{1,2,…,N}。
状态转移矩阵与二阶状态转移张量可用式(2.1-5)求解：
式(2.1-5)中，表示公式推导使用的中间变量，没有物理意义，T(t)表示t时刻以纬度幅角为自变量的无量纲化坐标到以时间为自变量的量纲化坐标的转换矩阵，Π表示相对状态分量次序变换矩阵，Σ(t)表示t时刻轨道要素偏差到相对状态的转换矩阵，表示平均轨道要素偏差到吻切轨道要素偏差的转换矩阵，表示t0时刻平均轨道要素偏差到t时刻平均轨道要素偏差的一阶传递矩阵，表示t0时刻轨道要素偏差的无量纲化矩阵，表示公式推导使用的中间变量的第i模块、第j行、第k列元素，A表示公式推导使用的中间变量，没有物理意义，B表示公式推导使用的中间变量，没有物理意义，Ail表示中间变量A的第i行、第l列元素，Blj表示中间变量B的第l行、第j列元素，Bmk表示中间变量B的第m行、第k列元素，表示t0时刻平均轨道要素偏差到t时刻平均轨道要素偏差的二阶传递张量的第l模块、第j行、第k列元素，Qilm(t)表示t时刻无量纲轨道要素偏差到相对状态的二阶转换张量的第i模块、第l行、第m列元素，表示t时刻轨道要素偏差的无量纲化矩阵。等带上标“bar”的矩阵或张量表示与平均轨道根数对应的值，即将主航天器对应时刻的平均轨道根数带入矩阵D，Γ，H表达式计算获得的值，反之没有上标“bar”的矩阵或张量表示与吻切轨道根数对应的值。以上表达式中有对应初始时刻t0及任意时刻t的量，分别带入t0时刻或t时刻下主航天器的轨道参数，即可计算对应的矩阵及张量。T为以纬度幅角θ为自变量的无量纲化坐标到以时间t为自变量的量纲化坐标的转换矩阵，反之量纲化坐标到无量纲化坐标的转换矩阵T-1为T的逆；Π为相对状态分量次序变换矩阵；Σ为轨道要素偏差到相对状态的转换矩阵；为初始时刻平均轨道要素偏差到任意时刻平均轨道要素偏差的一阶传递矩阵；D为平均轨道要素偏差到吻切轨道要素偏差的转换矩阵；Γ为轨道要素偏差的无量纲化矩阵；P为无量纲轨道要素偏差到相对状态的一阶转换矩阵；Q无量纲轨道要素偏差到相对状态的二阶转换张量；H为初始时刻平均轨道要素偏差到任意时刻平均轨道要素偏差的二阶传递张量。需要说明的是，前述矩阵T，T-1，Π，Γ，P及张量Q、H均为已知矩阵，其详细表达式参见文献[1]:Sengupta P,Vadali S R,Alfriend K T.Second-order state transition for relative motion near perturbed,elliptic orbits[J].Celestial Mechanics and Dynamical Astronomy,2006,97(2)∶101-129。前述矩阵Σ，D同样也均为已知矩阵，其详细表达式参见文献[2]:Gim D W,Alfriend K T.State Transition Matrix of Relative Motion for the Perturbed Noncircular Reference Orbit[J].Journal of Guidance,Control,and Dynamics,2003,26(6)∶956–971。
利用相对运动方程将初始相对状态x0向后分别预报到每个历史数据对应的历元时刻，即[t-n,t-(n-1),...,t-1]，获得预报的相对运动状态向量：若仅知道主航天器的轨道半长轴a(t0)(或平均轨道角速度利用C-W方程(1a)预报；若主航天器的6个轨道根数E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]均已知，则利用考虑J2摄动的非线性相对运动方程(2.2-2)预报。
本实施例中，步骤3)的详细步骤包括：
3.1)定义最小二乘方法的目标函数为式(3.1-1)所示的残差和，使得式(3.1-1)所示残差和最小的条件为式(3.1-2)；
式(3.1-1)中，J表示各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量的残差和，Y表示各历史历元时刻的预报相对运动状态向量，Z表示观察向量；
式(3.1-2)中，J表示各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量的残差和，x0表示t0时刻两航天器的初始相对运动状态，Y表示各历史历元时刻的预报相对运动状态向量，e表示各历史历元时刻的预报结果和历史轨道数据中历史数据的残差向量；
3.2)将相对运动方程、各历史历元时刻的预报相对运动状态向量和观察向量之间的残差向量带入式(3.1-2)，取一阶近似可得初始相对运动状态的最优估计值满足式(3.2-1)所示的迭代条件进行迭代求解；
式(3.2-1)中，k为迭代次数，表示第k+1次迭代的结果，表示第k次迭代的结果，的迭代初值x(t0)为时刻t0从航天器相对主航天器的初始相对运动状态，Fk向量的表达式为Fk＝[(Φk)TΦk]-1(Φk)T，Φk表示第k次迭代中的一阶状态转移矩阵，由式(2.1-1)及(2.1-2)可知，Φk在每次迭代中为常量；
3.3)判断是否满足指定的迭代条件，如果不满足指定迭代条件，则跳转执行步骤3.2)；否则当满足指定迭代条件时，将最后一次迭代得到的第k+1次迭代的结果作为最终使得主航天器和从航天器之间相对轨迹预报残差的残差和最小的最优初始相对运动状态
本实施例中，步骤3.3)中指定的迭代条件具体是指第k+1次迭代的结果与第k次迭代的结果之间的变化量小于或等于给定迭代误差、或者迭代次数k达到给定的最大值K。例如或k≤20；则结束迭代，获得最优初始相对运动状态的估计值：
将基于最小二乘方法计算得到的新初始条件带入原相对运动方程，向前预报到需要的历元时刻，获得具有较高精度的相对运动轨迹预报结果。
本实施例中，步骤4)获得指定的历元时刻tf的相对运动轨迹预报结果如式(4-1)或式(4-2)所示；
式(4-1)中，x(tf)表示指定的历元时刻tf的相对运动轨迹预报结果，Φ(tf,t0)是从t0时刻到tf时刻的状态转移矩阵，表示步骤3.3)中获得的最优初始相对运动状态；
式(4-2)中，xi(tf)表示指定的历元时刻tf的相对运动轨迹预报结果，Φ(tf,t0)表示是从t0时刻到tf时刻的状态转移矩阵，Ψ(tf,t0)表示从t0时刻到tf时刻的二阶状态转移张量，表示步骤3.3)中获得的最优初始相对运动状态。其中，方程(4-1)为基于线性C-W方程的预报，方程(4-2)为考虑J2摄动的非线性相对运动方程预报结果。仿真结果表明，对同样的相对运动方程，采用作为初始条件比直接采用x(t0)作为初始条件的相对轨迹预报精度要高。而对同样的历史数据，基于方程(4-2)的预报结果要比基于方程(4-1)的预报结果精度高。
综上所述，本实施例基于最小二乘的非线性相对轨迹预报方法为一种基于最小二乘的非线性相对轨迹预报精度改进方法，利用历史轨道数据，基于最小二乘原理最小化预报值与历史值的残差和，对相对轨迹预报的初始条件进行估计改进，以提高相对轨迹预报的精度。其中，所用的相对运动方程可根据具体情况选择线性C-W方程或考虑J2摄动的非线性相对运动方程，考虑了J2摄动项和二阶非线性项、可用于相距较远航天器的长时间、高精度相对轨迹预报，具有设计方法正确合理、对实际工程任务的适用性好等优点。
实施例二：
本实施例与实施例一基本相同，其主要不同点为：本实施例中，步骤2.1)中利用相对运动方程计算各历史历元时刻的预报相对运动状态向量时，仅已知主航天器的轨道半长轴a(t0)或平均轨道角速度其中μ为地心引力常数，a为半长轴，且利用式(2.1-1)所示C-W方程计算各历史历元时刻[t-n,t-(n-1),...,t-1]的预报相对运动状态向量；步骤4)获得指定的历元时刻tf的相对运动轨迹预报结果如式(4-1)所示。
本实施例仅根据历史相对运动状态数据、初始相对运动状态及主航天器半长轴(或平均轨道角速度)就可以进行相对轨迹预报。预报过程为解析的线性C-W方程，通过最小二乘估计新的初始条件能最小化预报误差，预报精度得到改进；且基于最小二乘的迭代过程均为解析计算，能够保证迭代过程快速收敛，计算效率高。参见图2所示相对轨迹预报位置误差对比可知，和现有技术采用的CW方程直接预报相比，本实施例(基于CW方程的改进)的相对距离误差x更小；参见图3所示相对轨迹预报速度误差对比可知，和现有技术采用的CW方程直接预报相比，本实施例(基于CW方程的改进)的相对速度误差x更小。
实施例三：
本实施例与实施例一基本相同，其主要不同点为：本实施例中，步骤2.1)中利用相对运动方程计算各历史历元时刻的预报相对运动状态向量时，已知主航天器的6个轨道根数E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]，且利用式(2.1-2)所示考虑J2摄动的非线性相对运动方程计算各历史历元时刻[t-n,t-(n-1),...,t-1]的预报相对运动状态向量；步骤4)获得指定的历元时刻tf的相对运动轨迹预报结果如式(4-2)所示。
本实施例用于初始条件中有历史相对运动状态数据及初始时刻相对运动状态，且主航天器的6个轨道根数E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]均已知的情况。因此在步骤2)及步骤4)中，采用考虑J2摄动的非线性相对运动方程(2.1-2)及(4-2)进行相对轨迹预报。
令编队飞行的两航天器分别为主航天器与从航天器，已知进行相对轨迹预报的初始时刻为t0,主航天器初始轨道要素为E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]，其中μ为地心引力常数，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近拱点角距，f为真近点角。从航天器相对主航天器的初始相对运动状态为x(t0)，其中x(t0)在主航天器当地轨道坐标系中表示，该坐标系原点为主航天器质心，x轴沿主航天器地心矢径方向，z轴沿轨道面法向，y轴与x、z轴构成右手坐标系。
本实施例根据历史相对运动状态数据、初始相对运动状态及主航天器初始轨道根数就可以进行相对轨迹预报。预报过程为解析的考虑J2摄动的非线性相对运动方程，通过最小二乘估计新的初始条件能最小化预报误差，预报精度得到改进；且基于最小二乘的迭代过程均为解析计算，能够保证迭代过程快速收敛，计算效率高。对同样的历史数据及初始条件，本实施例所得精度比实施例二要高。
参见图2所示相对轨迹预报位置误差对比可知，和现有技术采用J2非线性方程直接预报相比，本实施例(采用J2非线性方程的改进)的相对距离误差x更小；参见图3所示相对轨迹预报速度误差对比可知，和现有技术采用J2非线性方程直接预报相比，本实施例(采用J2非线性方程的改进)的相对速度误差x更小。
实施例四：
本实施例用于初始条件中没有历史相对运动状态数据，仅知道初始时刻相对运动状态，及主航天器的6个轨道根数E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]的情况。该情况下，直接采用考虑J2摄动的非线性相对运动方程(2.1-2)进行相对运动轨迹预报。
令编队飞行的两航天器分别为主航天器与从航天器，已知进行相对轨迹预报的初始时刻为t0,主航天器初始轨道要素为E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]，其中μ为地心引力常数，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近拱点角距，f为真近点角。从航天器相对主航天器的初始相对运动状态为x(t0)，其中x(t0)在主航天器当地轨道坐标系中表示，该坐标系原点为主航天器质心，x轴沿主航天器地心矢径方向，z轴沿轨道面法向，y轴与x、z轴构成右手坐标系。
将初始轨道要素为E(t0)＝[a(t0),e(t0),i(t0),Ω(t0),ω(t0),f(t0)]及从航天器相对主航天器的初始相对运动状态为x(t0)直接带入考虑J2摄动的非线性相对运动方程(2.1-2)，即可获得任意时刻的相对运动状态预报结果。
本实施例在没有历史相对运动状态数据的情况下就可以进行相对轨迹预报。预报过程为解析的考虑J2摄动的非线性相对运动方程，计算效率高。对同样的历史数据及初始条件，本实施例所得精度比实施例三要低。
以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
- 还没有人留言评论。精彩留言会获得点赞！