低轨卫星J2模型数值外推精度竟然接近STK的HPOP高精度外推模型

最新推荐文章于 2024-04-14 11:22:05 发布

phenixgemini

最新推荐文章于 2024-04-14 11:22:05 发布

阅读量6.6k

点赞数 7

分类专栏：卫星轨道文章标签：算法 python

本文链接：https://blog.csdn.net/phenixgemini/article/details/127161728

版权

卫星同时被 2 个专栏收录

1 篇文章 2 订阅

订阅专栏

轨道

1 篇文章 0 订阅

订阅专栏

1.二体运动方程

卫星在轨运行遵循万有引力定律。在二体模型下，以地心惯性坐标系为参考，卫星所受地球引力为

$\vec{F}=-\frac{GMm}{r^{^{3}}}\vec{r}$

其中

卫星的位置矢量为 $\vec{r}=\left ( x,y,z \right )^{^{T}}$ ；

速度矢量为 $\vec{v}=\left ( {v_{x},v_{y},v_{z}} \right )^{T}$ ；

万有引力常数为G；

地球质量为M；

卫星质量为m，

卫星矢量的模为 $r=\sqrt{x^2+y^2+z^2}$

根据牛顿第二定律，二体问题的微分方程为

$\frac{d\ddot{\vec{r}}}{dt}=\frac{\vec{F}}{m}=-\mu \frac{\vec{r}}{r^{3}}$

其中,

$\mu=GM$ 为引力常量。

由于 $\dot{\vec{r}}=\vec{v}$ , $\dot{\vec{v}}=\ddot{\vec{r}}$ ，则可构成状态方程。将上式写为直角坐标系分量形式如下，则

$\left\{\begin{matrix}\dot{x}=v_x \\\dot{y}=v_y \\\dot{z}=v_z \\\dot{v}_x=-\mu\frac{x}{r^3} \\\dot{v}_y=-\mu\frac{y}{r^3} \\\dot{v}_z=-\mu\frac{z}{r^3} \end{matrix}\right.$

上式为六个未知数，六个方程，若给定初始条件为 $\vec{r}_0,\vec{v}_0$ ，恰好可以满足方程可解。

2.J2模型下的摄动方程

低轨卫星所处的空间环境很复杂，要承受地球非球形引力，大气阻力，日月引力及太阳光压等多种摄动因素的影响。考虑J2摄动下的加速度，则上述状态方程变为

$\left\{\begin{matrix}\dot{x}=v_x \\\dot{y}=v_y \\\dot{z}=v_z \\\dot{v}_x=-\mu\frac{x}{r^3}+\Delta g_x \\\dot{v}_y=-\mu\frac{y}{r^3}+\Delta g_y \\\dot{v}_z=-\mu\frac{z}{r^3}+\Delta g_z \end{matrix}\right.$

其中

$\left\{\begin{matrix}\Delta g_x=-\frac{3}{2}J_2{R_e}^{2}\mu\frac{x}{r^5} \\ \Delta g_y=-\frac{3}{2}J_2{R_e}^{2}\mu\frac{y}{r^5} \\ \Delta g_z=-\frac{3}{2}J_2{R_e}^{2}\mu[3-5(\frac{z}{r})^2] \end{matrix}\right.$

3.龙格库塔算法

当 $a\leq t\leq b$ ，对初值问题

$\left\{\begin{matrix}\dot{r}=f(t,r) \\ r(a)=r_0 \end{matrix}\right.$

常用四阶龙格库塔算法，直接给出算法如下

$\left\{\begin{matrix}r_{m+1}=r_m+\frac{h}{6}\left ( K_1+2K_2+2K_3+K_4 \right ) \\ K_1=f(t_m,r_m) \\ K_2=f(t_m+\frac{h}{2},r_m+\frac{h}{2}K_1) \\ K_3=f(t_m+h,r_m+\frac{h}{2}K_2) \\ K_4=f(t_m+h,r_m+hK_3) \end{matrix}\right.$

4.STK轨道

本文为了验证数值外推与STK生成轨道的精度，选取轨道参数为

时间：2022-10-01 00:00:00（UTCG）

半长轴a:6878.137km

偏心率e:0.001

倾角i:60.0

升交点赤经RAAN:12

近地点幅角omega:0

平近点角M:300

在J2模型和HPOP模型下的参数设置如下。除六根数设置外，其他均为默认参数。至于HPOP和J2模型的区别，简而言之，可认为HPOP模型是高精度外推模型就好了。

将此参数输入STK，在J2模型下，得到卫星位置和速度的一组初值

$\vec{r}_0=\begin{Bmatrix}3971.676026 \\ -2202.172866 \\ -5161.178823 \end{Bmatrix}, \vec{v}_0=\begin{Bmatrix}6.059801 \\3.231769 \\3.293050 \end{Bmatrix}$

在HPOP模型下，得到卫星位置和速度的一组初值仍然如下值

$\vec{r}_0=\begin{Bmatrix}3971.676026 \\ -2202.172866 \\ -5161.178823 \end{Bmatrix}, \vec{v}_0=\begin{Bmatrix}6.059801 \\3.231769 \\3.293050 \end{Bmatrix}$

分别以初值作为方程的初值进行数值计算。利用Python进行编程计算。代码见第五章。

可见以同样的初值 $\vec{r}_0,\vec{v}_0$ 进行外推，利用J2模型进行数值外推，一天内的误差与STK的J2模型外推误差竟然远大于STK的HPOP模型外推误差。这与常规理解相左。理论上不论什么工具，同样的J2外推误差应该近似，但是在数值外推的情况下，只考虑J2项摄动进行数值外推，竟然与STK的高精度模型HPOP外推得到的结果类似。至于为什么出现这样的情况，尚未查明。

J2数值外推与STK 的J2、HPOP模型外推位置和速度误差如下所示。

STK轨道模型	Rerror	Verror	备注
J2	<514km	<0.558km/s	外推1天的误差
HPOP	<18.1km	<0.0202km/s	外推1天的误差

研究卫星轨道时通常采用STK的HPOP模型进行高精度轨道外推。但在一些小型仿真案例中，如果不能使用STK，根据上面的仿真结果，可用J2模型数值外推，根据需要，在短时间内仍然可以得到与STK可比拟的精度。

5.python代码

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams["font.sans-serif"]="SimHei"
plt.rcParams["axes.unicode_minus"]=False
mu=398600
Re=6378.137# 地球半径
J2=0.00108263# J2项

##状态方程
def StateEq(t,RV):
    x=RV[0]
    y=RV[1]
    z=RV[2]
    vx=RV[3]
    vy=RV[4]
    vz=RV[5]
    r=np.sqrt(x**2+y**2+z**2)
    gx=-mu * x / r ** 3
    gy=-mu * y / r ** 3
    gz=-mu * z / r ** 3
    dgx = -3 / 2 * J2 * Re ** 2 * mu * x / r ** 5 * (1 - 5 * (z / r) ** 2)
    dgy = -3 / 2 * J2 * Re ** 2 * mu * y / r ** 5 * (1 - 5 * (z / r) ** 2)
    dgz = -3 / 2 * J2 * Re ** 2 * mu * z / r ** 5 * (3 - 5 * (z / r) ** 2)
    f=np.array([vx,vy,vz,gx+dgx,gy+dgy,gz+dgz])
    return f
## 龙格库塔算法
def RungeKutta(t0,r0,h):
    K1=StateEq(t0,r0)
    K2=StateEq(t0 + 2 / h, r0 + h / 2 * K1)
    K3=StateEq(t0 + h, r0 + h / 2 * K2)
    K4 = StateEq(t0 + h, r0 + h * K3)
    r1=r0 + h / 6 * (K1 + 2 * K2 + 2 * K3 + K4)
    return r1

h=1# 步长
N=86400 # 外推时间，一天86400秒
# 初始位置和速度
RV0=np.array([3971.676026,-2202.172866, -5161.178823,6.059801,3.231769,3.293050])
data=np.array(np.zeros((np.int32(N/h+1),6)))
data[0]=RV0
for i in range(np.int32(N/h)):
    data[i+1]=RungeKutta(i*h,RV0,h)
    RV0=data[i+1]
##
df=pd.read_csv(r"J2_J2000_Position_Velocity.csv")
error=df.iloc[:,1:]-data
plt.subplot(2,2,1)
plt.plot(error.iloc[:,0:3],label=['Xerror','Yerror','Zerror'])
plt.xlabel('时间/s')
plt.ylabel('位置误差/km')
plt.legend()
plt.grid()
plt.title('J2数值外推与STK的J2模型外推的位置误差')
plt.subplot(2,2,2)
plt.plot(error.iloc[:,3:6],label=['VXerror','VYerror','VZerror'])
plt.xlabel('时间/s')
plt.ylabel('速度误差/(km/s)')
plt.legend()
plt.grid()
plt.title('J2数值外推与STK的J2模型外推的速度误差')

df=pd.read_csv(r"HPOP_J2000_Position_Velocity.csv")
error=df.iloc[:,1:]-data
plt.subplot(2,2,3)
plt.plot(error.iloc[:,0:3],label=['Xerror','Yerror','Zerror'])
plt.xlabel('时间/s')
plt.ylabel('位置误差/km')
plt.legend()
plt.grid()
plt.title('J2数值外推与STK的HPOP模型外推的位置误差')
plt.subplot(2,2,4)
plt.plot(error.iloc[:,3:6],label=['VXerror','VYerror','VZerror'])
plt.xlabel('时间/s')
plt.ylabel('速度误差/(km/s)')
plt.legend()
plt.grid()
plt.title('J2数值外推与STK的HPOP模型外推的速度误差')