常见的轨道外推模型

最近研究卫星轨道需要用到一些推导外推模型，我一般会使用STK的HPOP模型进行高精度轨道外推，但因为一些原因不能使用STK，在此分享几个比较常见的轨道外推模型：

1. 轨道外推模型

1.1 拉格朗日系数法

该函数我是从一个类中直接取出来的，如果要使用需要提前将卫星的位置速度转换至轨道六根数输入。具体求法可以参考国防科大张洪波老师出版的《航天器轨道力学理论与方法》书中P85-P93。

def calculation_spacecraft_status(self, t):
        """
        :param t: 航天器轨道状态信息外推时间
        :return: 航天器轨道状态信息
        """
        def delta_E_equation_t(delta_E):
            # 根据当前时刻的平近点角差确定偏近点角差()
            eq = delta_E + sigma_t * (1 - np.cos(delta_E)) / np.sqrt(self.a_t) - \
                 (1 - self.r_t / self.a_t) * np.sin(delta_E) - np.sqrt(self.u / self.a_t**3) * t
            return eq

        def delta_E_equation_c(delta_E):
            # 根据当前时刻的平近点角差确定偏近点角差()
            eq = delta_E + sigma_c * (1 - np.cos(delta_E)) / np.sqrt(self.a_c) - \
                 (1 - self.r_c / self.a_c) * np.sin(delta_E) - np.sqrt(self.u / self.a_c**3) * t
            return eq

        def Lagrange_factor_c(sigma, delta_E):
            # 拉格朗日系数,外推轨道位置速度
            r = self.a_c + (self.r_c - self.a_c) * np.cos(delta_E) + sigma * np.sqrt(self.a_c) * np.sin(delta_E)
            F = 1 - self.a_c * (1 - np.cos(delta_E)) / self.r_c
            G = self.a_c * sigma * (1 - np.cos(delta_E)) / np.sqrt(self.u) + \
                self.r_c * np.sqrt(self.a_c / self.u) * np.sin(delta_E)
            F_c = -np.sqrt(self.u * self.a_c) * np.sin(delta_E) / (r * self.r_c)
            G_c = 1 - self.a_c * (1 - np.cos(delta_E)) / r

            return F, G, F_c, G_c

        # 偏近点角差
        sigma_t = np.dot(self.R0_t, self.V0_t) / np.sqrt(self.u)
        sigma_c = np.dot(self.R0_c, self.V0_c) / np.sqrt(self.u)
        # 当前时刻处的偏近点角差
        delta_E_t = fsolve(delta_E_equation_t, 0)
        delta_E_c = fsolve(delta_E_equation_c, 0)
        # print(self.R0_t, self.V0_t, self.R0_c, self.V0_c)
        # print(delta_E_t, sigma_c, self.a_c, self.r_c / self.a_c, np.sqrt(self.u / self.a_c ** 3) * t)
        #
        # 航天器位置速度外推(惯性系下)
        F, G, F_t, G_t = self.Lagrange_factor(sigma_t, delta_E_t)
        R_i_t = F * self.R0_t + G * self.V0_t
        V_i_t = F_t * self.R0_t + G_t * self.V0_t

        F, G, F_c, G_c = Lagrange_factor_c(sigma_c, delta_E_c)
        R_i_c = F * self.R0_c + G * self.V0_c
        V_i_c = F_c * self.R0_c + G_t * self.V0_c

        return np.array([R_i_c, V_i_c]).ravel(), np.array([R_i_t, V_i_t]).ravel()

1.2 四阶龙格库塔算法

import numpy as np

# 常数定义
G = 6.67430e-11  # 万有引力常数 (m^3 kg^−1 s^−2)
M_earth = 5.972e24  # 地球质量 (kg)
time_step = 1.0  # 时间步长，单位为秒

# 计算引力加速度的函数
def acceleration(r):
    r_magnitude = np.linalg.norm(r)  # 计算当前位置的距离（米）
    return -G * M_earth / r_magnitude ** 3 * r  # 返回加速度（m/s^2）


# 四阶龙格-库塔方法进行轨道推演
def runge_kutta_orbit(r0, v0, time_step, total_time):
    r = r0
    v = v0

    # 记录轨道数据
    positions = [r]
    velocities = [v]

    # 计算总时间步数
    num_steps = int(total_time / time_step)

    for _ in range(num_steps):
        # 计算四阶龙格-库塔法的k值
        k1_v = acceleration(r)
        k1_r = v

        k2_v = acceleration(r + 0.5 * time_step * k1_r)
        k2_r = v + 0.5 * time_step * k1_v

        k3_v = acceleration(r + 0.5 * time_step * k2_r)
        k3_r = v + 0.5 * time_step * k2_v

        k4_v = acceleration(r + time_step * k3_r)
        k4_r = v + time_step * k3_v

        # 计算新的位置和速度
        v = v + (time_step / 6) * (k1_v + 2 * k2_v + 2 * k3_v + k4_v)
        r = r + (time_step / 6) * (k1_r + 2 * k2_r + 2 * k3_r + k4_r)

        # 存储轨道数据
        positions.append(r)
        velocities.append(v)

    return np.array(positions), np.array(velocities)

# 初始状态
r0 = np.array([7378137, 0, 0])  # 卫星初始位置 (m)
v0 = np.array([0, 7350.1, 0])  # 卫星初始速度 (m/s)

# 推演卫星轨道
total_time = 100
positions, velocities = runge_kutta_orbit(r0, v0, time_step, total_time)

# 打印结果
print("卫星轨道位置（m）：", positions[-1])
print("卫星速度（m/s）：", velocities[-1])

1.3 欧拉法

import numpy as np

# 常数定义
G = 6.67430e-11  # 万有引力常数 (m^3 kg^−1 s^−2)
M_earth = 5.972e24  # 地球质量 (kg)
m_satellite = 1.0  # 卫星质量 (kg)
time_step = 1.0  # 时间步长，单位为秒

def satellite_orbit(r0, v0, time_step, total_time):
    # 初始位置和速度
    r = r0
    v = v0

    # 记录轨道数据
    positions = [r]
    velocities = [v]

    # 计算总时间步数
    num_steps = int(total_time / time_step)

    for _ in range(num_steps):
        # 计算引力方向
        r_magnitude = np.linalg.norm(r)
        acceleration = -G * M_earth / r_magnitude ** 3 * r  # 引力加速度 (m/s^2)

        # 更新速度和位置
        v = v + acceleration * time_step
        r = r + v * time_step

        # 存储轨道数据
        positions.append(r)
        velocities.append(v)

    return np.array(positions), np.array(velocities)

# 初始状态
r0 = np.array([7378137, 0, 0])  # 卫星初始位置 (m)
v0 = np.array([0, 7350.1, 0])  # 卫星初始速度 (m/s)

# 推演卫星轨道
total_time = 100 # 推演时间
positions, velocities = satellite_orbit(r0, v0, time_step, total_time)

print("卫星轨道位置（m）：", positions[-1])
print("卫星速度（m/s）：", velocities[-1])

1.4 CW方程

import numpy as np

# 常数定义
G = 6.67430e-11  # 万有引力常数 (m^3 kg^−1 s^−2)
M_earth = 5.972e24  # 地球质量 (kg)

# 计算引力加速度的函数
def acceleration(r):
    r_magnitude = np.linalg.norm(r)  # 计算当前位置的距离（米）
    return -G * M_earth / r_magnitude ** 3 * r  # 返回加速度 (m/s^2)

# Cowell方法轨道推演
def cowell_orbit(r0, v0, total_time, time_step=1.0):
    r = r0
    v = v0

    # 存储轨道数据
    positions = [r]
    velocities = [v]

    # 计算总时间步数
    num_steps = int(total_time / time_step)

    for _ in range(num_steps):
        # 计算加速度
        a = acceleration(r)

        # 更新速度和位置（使用欧拉法）
        v = v + a * time_step  # 更新速度
        r = r + v * time_step  # 更新位置

        # 存储轨道数据
        positions.append(r)
        velocities.append(v)

    return np.array(positions), np.array(velocities)

# 初始状态
r0 = np.array([7378137, 0, 0])  # 卫星初始位置 (m)
v0 = np.array([0, 7350.1, 0])  # 卫星初始速度 (m/s)

# 推演卫星轨道
total_time = 100  # 推演1小时
time_step = 1.0  # 时间步长，单位为秒
positions, velocities = cowell_orbit(r0, v0, total_time, time_step)

# 打印结果
print("卫星轨道位置（m）：", positions[-1])
print("卫星速度（m/s）：", velocities[-1])

1.5 比较

使用上面四个模型和STK在相同输入数据以及相同外推时间下的结果对比：

拉格朗日：[ 7.34155608e+06  7.33794867e+05  0.00000000e+00 -7.31013347e+02 7.31365808e+03  0.00000000e+00]
龙格库塔：[7341557.06851583  733794.90050119     0 -730.99346743 7313.65907585    0]
CW：[7341189.75673622  733794.93144015       0  -731.04732227 7314.01962396    0]
欧拉法：[7341189.75673622  733794.93144015       0    -731.04732227 7314.01962396    0 ]

STK：[7341493.332, 733973.096, -0.544, -732.087, 7313.590, -0.011]

可以看出来，误差是有的，但是误差在短时间的外推下还是比较小。