J2摄动条件下的轨道转移方程

J2摄动条件下的轨道转移方程

引用

CSDN

1.

https://blog.csdn.net/weixin_57997461/article/details/137403155

一、背景介绍

在进行精密定轨的过程中，通常涉及到求解卫星自身的状态转移方程，根据考虑不同的摄动力会对应着不同的状态转移方程，在本文中仅考虑J2项摄动影响，讨论卫星的状态转移方程。

根据牛顿第二定律以及J2项摄动情况下的势函数，可以得到下式

式中，
为地球的势函数，
为地球引力常数，
为该位置离地球质心的距离，
为地球赤道半径，
为该位置对应的纬度

根据状态转移方程，设卫星的状态空间为
，其中
为系统观测误差向量，则系统运动方程满足：

得到的状态转移方程的微分形式
为

由于动力学模型的非线性，系统状态转移方程
不能够解析给出，但给出其满足的微分方程。

需要求出
时刻后的状态位置
，以及从初始状态到该末状态的状态转移方程，对于以上非线性动力学系统，通常采用数值的方法来求解。

对于数值积分，首先将输入变量的时间函数插入到状态方程中，给出
向量微分方程

有几种数值积分方法可用，这已被广泛描述，如Butcher。问题的求解取决于初始状态
，因此被称之为初值问题，用于解决初值问题的各种方法可以细分为一步法、多步法和外推法。著名的是四阶龙格库塔法，取步长为
伴随向量为

为了更好的理解这些公式，可以将上式应用于齐次定常系统
推得单步误差次数为
。一步法中最重要的是增量
，它对积分误差有直接的影响。特别是，我们可以通过一更小的增量计算来提高结果的准确性

二、具体实现

2.1 线性化过程

根据文章求出的状态微分方程，明显是一个非线性的函数，因此，需要对其进行线性化，这个过程通过MATLAB的syms系统进行求解，可以直接求出线性化后的状态方程，公式如下所示

syms x y z J2 ae miu
r=sqrt(x^2+y^2+z^2)
dvx=-miu*x/r^3-1.5*J2*ae^2*miu*x/r^5*(1-5*(z/r)^2);
dvy=-miu*y/r^3-1.5*J2*ae^2*miu*y/r^5*(1-5*(z/r)^2);
dvz=-miu*z/r^3-1.5*J2*ae^2*miu*z/r^5*(3-5*(z/r)^2);
dvxdx=simplify(diff(dvx,x));
dvxdy=simplify(diff(dvx,y));
dvxdz=simplify(diff(dvx,z));
dvydx=simplify(diff(dvy,x));
dvydy=simplify(diff(dvy,y));
dvydz=simplify(diff(dvy,z));
dvzdx=simplify(diff(dvz,x));
dvzdy=simplify(diff(dvz,y));
dvzdz=simplify(diff(dvz,z));  

function GAMMA=trans(x,y,z,J2,ae,miu)
GAMMA=zeros(3);
GAMMA(1,1)=-(2*miu*(x^2 + y^2 + z^2)^3 - 6*miu*x^2*(x^2 + y^2 + z^2)^2 + 3*J2*ae^2*miu*(x^2 + y^2 + z^2)^2 + 105*J2*ae^2*miu*x^2*z^2 - 15*J2*ae^2*miu*x^2*(x^2 + y^2 + z^2) - 15*J2*ae^2*miu*z^2*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
GAMMA(1,2)=(6*miu*x*y*(x^2 + y^2 + z^2)^2 - 105*J2*ae^2*miu*x*y*z^2 + 15*J2*ae^2*miu*x*y*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
GAMMA(1,3)=(6*miu*x*z*(x^2 + y^2 + z^2)^2 - 105*J2*ae^2*miu*x*z^3 + 45*J2*ae^2*miu*x*z*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
GAMMA(2,1)=(6*miu*x*y*(x^2 + y^2 + z^2)^2 - 105*J2*ae^2*miu*x*y*z^2 + 15*J2*ae^2*miu*x*y*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
GAMMA(2,2)=-(2*miu*(x^2 + y^2 + z^2)^3 - 6*miu*y^2*(x^2 + y^2 + z^2)^2 + 3*J2*ae^2*miu*(x^2 + y^2 + z^2)^2 + 105*J2*ae^2*miu*y^2*z^2 - 15*J2*ae^2*miu*y^2*(x^2 + y^2 + z^2) - 15*J2*ae^2*miu*z^2*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
GAMMA(2,3)=(6*miu*y*z*(x^2 + y^2 + z^2)^2 - 105*J2*ae^2*miu*y*z^3 + 45*J2*ae^2*miu*y*z*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
GAMMA(3,1)=(6*miu*x*z*(x^2 + y^2 + z^2)^2 - 105*J2*ae^2*miu*x*z^3 + 45*J2*ae^2*miu*x*z*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
GAMMA(3,2)=(6*miu*y*z*(x^2 + y^2 + z^2)^2 - 105*J2*ae^2*miu*y*z^3 + 45*J2*ae^2*miu*y*z*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
GAMMA(3,3)=-(2*miu*(x^2 + y^2 + z^2)^3 - 6*miu*z^2*(x^2 + y^2 + z^2)^2 + 9*J2*ae^2*miu*(x^2 + y^2 + z^2)^2 + 105*J2*ae^2*miu*z^4 - 90*J2*ae^2*miu*z^2*(x^2 + y^2 + z^2))/(2*(x^2 + y^2 + z^2)^(9/2));
end  

根据线性化结果和龙格库塔法的关系，分别使用状态转移的方法和数值积分的方法，对J2轨道进行预报。第一个代码为J2摄动的
，写成如下函数

function x1=RK_J2(x0,t)
x1=zeros(6,1)
J2=1.08263e-3;
ae=6378.137;
miu=3.986e5;
x=x0(1);
y=x0(2);
z=x0(3);
vx=x0(4);
vy=x0(5);
vz=x0(6);
r=sqrt(x^2+y^2+z^2);
x1(1)=vx;x1(2)=vy;x1(3)=vz;
x1(4)=-miu*x/r^3-1.5*J2*ae^2*miu*x/r^5*(1-5*(z/r)^2);
x1(5)=-miu*y/r^3-1.5*J2*ae^2*miu*y/r^5*(1-5*(z/r)^2);
x1(6)=-miu*z/r^3-1.5*J2*ae^2*miu*z/r^5*(3-5*(z/r)^2);
end  

使用四阶龙格库塔法数值积分，步长为10s

y=[23763.0112742573,-14408.2177617449,-2541.17345408073,0.911182669796365, 1.04987004575796,2.33432272229561];
h=10;
dx1=RK_J2(m,0)
dx2=RK_J2(m+dx1,0.5*h);
dx3=RK_J2(m+dx2,0.5*h);
dx4=RK_J2(m+dx3,h);
mmm=m+h/6*(dx1+2*dx2+2*dx3+dx4)  

使用状态转移方程求解

GM=3.986e5;
Re=6378.137;
A=zeros(6);
A(1,4)=1;
A(2,5)=1;
A(3,6)=1;
gamma=trans(m(1),m(2),m(3),J_2,Re,GM);
for i=1:3
    for j=1:3
        A(3+i,j)=gamma(i,j);
    end
end
% 取步长为10s
E=eye(6);h=10
Phi=E+A.*h+1/2*A*A.*h^2+1/6*A*A*A.*h^3+1/24*A*A*A*A.*h^4;
m1=Phi*m;  

三、结果对比

使用STK的J2pertubation预报器进行对比，其中1为STK预报器，2为状态转移方程方法，3为四阶龙格库塔法数值积分，下表为对初始位置10s后位置和速度的预报

比较上式可以知道，状态转移方法相比于四阶龙格库塔数值积分法，精度存在明显的不足，需要进一步提高精度。