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

矢志航天的阿洪

已于 2024-04-10 08:02:07 修改

阅读量1.4k

点赞数 5

文章标签：算法

于 2024-04-09 22:52:58 首次发布

本文链接：https://blog.csdn.net/weixin_57997461/article/details/137403155

版权

本文探讨了在卫星定轨中考虑J2摄动时的状态转移方程，通过线性化和数值积分方法（如四阶龙格库塔法）进行求解。与STK预报器的对比显示，状态转移方程方法精度较低，需改进。

摘要由CSDN通过智能技术生成

一、背景介绍

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

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

$\ddot{\boldsymbol{r}}=\boldsymbol{g}=\textup{grad}\boldsymbol{U}=\bigtriangledown\boldsymbol{ U}\\ U=\frac{\mu _e}{r}[1-\frac{J_2}{2}(\frac{a_e}{r})^2(3\sin^2(\varphi)-1)]$

式中， $\boldsymbol{U}$ 为地球的势函数， $\mu_e$ 为地球引力常数， $r$ 为该位置离地球质心的距离， $a_e$ 为地球赤道半径， $\varphi$ 为该位置对应的纬度

根据状态转移方程，设卫星的状态空间为 $\boldsymbol{X}=[\boldsymbol{r},\boldsymbol{v},\boldsymbol{q}]^T$ ，其中 $\boldsymbol{q}$ 为系统观测误差向量，则系统运动方程满足：

$\boldsymbol{\dot{X}}=[\boldsymbol{\dot{r}},\boldsymbol{\dot{v}},\boldsymbol{\dot{q}}]^T=f(\boldsymbol{X},t)=[v,\bigtriangledown U,0]^T$

得到的状态转移方程的微分形式 $f(\boldsymbol{X},t)$ 为

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

$\dot{{\boldsymbol{\Phi}}}(t,t_0)=\boldsymbol{A}(t)\boldsymbol{\Phi}(t,t_0),\boldsymbol{\Phi}(t_0,t_0)=\boldsymbol{I},\boldsymbol{A}(t)=\frac{\partial{f(\boldsymbol{X},t)}}{\partial{\boldsymbol{X}}}$

$(\dot{x},\dot{y},\dot{z},\dot{v_{x0}},\dot{v_{y0}},\dot{v_{z0}})^{T}=\begin{bmatrix} 0&0 &0 & 1 & 0 &0 \\ 0& 0 &0 & 0 & 1&0 \\ 0& 0 & 0 &0 & 0& 1\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_x}{\partial y} & \frac{\partial a_x}{\partial z} & 0 & 0 &0 \\ \frac{\partial a_y}{\partial x}& \frac{\partial a_y}{\partial y}& \frac{\partial a_y}{\partial z} & 0& 0& 0\\ \frac{\partial a_z}{\partial x} & \frac{\partial a_z}{\partial y} & \frac{\partial a_z}{\partial z} &0 &0 & 0 \end{bmatrix}(x_0,y_0,z_0,v_{x0},v_{y0},v_{z0})^{T}$

需要求出 $t_0$ 时刻后的状态位置 $(x_i,y_i,z_i,v_{xi},v_{yi},v_{zi})$ ，以及从初始状态到该末状态的状态转移方程，对于以上非线性动力学系统，通常采用数值的方法来求解。

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

$\boldsymbol{\dot{{x}}}(t)=\boldsymbol{a}(\boldsymbol{x},t),\boldsymbol{x}(t_0)=\boldsymbol{x}_0$

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

$\boldsymbol{\dot{{x}}}(t_{k+1})=\boldsymbol{x}(t_k)+\frac{h}{6}[\Delta \boldsymbol{x}_k^{(1)}+2\Delta \boldsymbol{x}_k^{(2)}+2\Delta \boldsymbol{x}_k^{(3)}+\Delta \boldsymbol{x}_k^{(4)}]$

伴随向量为

$\Delta \boldsymbol{x}_k^{(1)}=\boldsymbol{a}(\boldsymbol{x}(t_k),t_k)\\ \Delta \boldsymbol{x}_k^{(2)}=\boldsymbol{a}(\boldsymbol{x}(t_k)+\frac{1}{2}\Delta \boldsymbol{x}_k^{(1)},t_k+\frac{1}{2}h)\\ \Delta \boldsymbol{x}_k^{(3)}=\boldsymbol{a}(\boldsymbol{x}(t_k)+\frac{1}{2}\Delta \boldsymbol{x}_k^{(2)},t_k+\frac{1}{2}h)\\ \Delta \boldsymbol{x}_k^{(4)}=\boldsymbol{a}(\boldsymbol{x}(t_k)+\Delta \boldsymbol{x}_k^{(3)},t_k+h)$

为了更好的理解这些公式，可以将上式应用于齐次定常系统

${\boldsymbol{\dot{x}}(t)}=\boldsymbol{A}\cdot \boldsymbol{x}(t)\\ x(t_{k+1})\approx\boldsymbol{ \phi}(h)\boldsymbol{x}(t_k)\\ \boldsymbol{ \phi}(h)=\boldsymbol{E}+\boldsymbol{A}h+\frac{1}{2}\boldsymbol{AA}h^2+\frac{1}{6}\boldsymbol{AAA}h^3+\frac{1}{24}\boldsymbol{AAAA}h^4$

推得单步误差次数为 $h^5$ 。一步法中最重要的是增量 $h$ ，它对积分误差有直接的影响。特别是，我们可以通过一更小的增量计算来提高结果的准确性

二、具体实现

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摄动的 $f(\boldsymbol{X},t)$ ，写成如下函数

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;