Python 卫星通信模拟：低轨星座的轨道力学计算
摘要
本文详细介绍了使用 Python 进行低轨卫星星座轨道力学计算的完整模拟方法。我们将从基础理论出发，逐步构建一个完整的轨道模拟系统，涵盖卫星轨道参数计算、星座构型设计、轨道摄动模型以及可见性分析等关键方面。
1. 引言
低地球轨道（LEO）卫星星座已成为全球通信、地球观测和导航系统的关键基础设施。星座由数十至数千颗卫星组成，如 Starlink、OneWeb 等系统。精确的轨道力学计算对于星座设计、卫星管理和通信链路维护至关重要。
Python 凭借其强大的科学计算库（如 NumPy、SciPy）和可视化工具（如 Matplotlib），成为轨道力学模拟的理想选择。本文将构建一个完整的轨道模拟框架，涵盖从基础理论到实际应用的各个方面。
2. 轨道力学基础理论
2.1 牛顿万有引力定律
轨道力学的核心是牛顿万有引力定律：
F = G * m1 * m2 / r^2
对于地球与卫星系统：
a = -mu * r / r^3
其中 mu = GM 是地球引力常数，约为 3.986×10^14 m^3/s^2。
2.2 运动方程
卫星在地球引力场中的运动方程为：
d^2r/dt^2 = -mu * r / r^3 + apert
其中 apert 表示各种摄动加速度。
3. 轨道参数与坐标系统
3.1 经典轨道要素
卫星轨道通常用六个轨道要素描述：
- 半长轴 (a)：轨道尺寸
- 偏心率 (e)：轨道形状
- 轨道倾角 (i)：轨道平面倾斜度
- 升交点赤经 (Ω)：轨道平面方向
- 近地点幅角 (ω)：轨道椭圆方向
- 真近点角 (ν)：卫星在轨道上的位置
3.2 坐标系统转换
轨道计算涉及多个坐标系统：
- 地心惯性坐标系 (ECI)
- 地固坐标系 (ECEF)
- 轨道平面坐标系
4. 二体问题与开普勒轨道计算
4.1 开普勒方程
对于椭圆轨道，开普勒方程为：
M = E - e * sin(E)
其中 M 是平近点角，E 是偏近点角。
4.2 位置与速度计算
根据轨道要素计算卫星位置和速度：
import numpy as np
from math import sin, cos, sqrt, radians, degrees
class KeplerOrbit:
"""开普勒轨道计算类"""
def __init__(self, a, e, i, raan, arg_peri, nu, mu=3.986004418e14):
"""
初始化轨道要素
参数:
a: 半长轴 (m)
e: 偏心率
i: 轨道倾角 (度)
raan: 升交点赤经 (度)
arg_peri: 近地点幅角 (度)
nu: 真近点角 (度)
mu: 地球引力常数 (m^3/s^2)
"""
self.a = a
self.e = e
self.i = radians(i)
self.raan = radians(raan)
self.arg_peri = radians(arg_peri)
self.nu = radians(nu)
self.mu = mu
def calculate_position_velocity(self):
"""计算卫星在地心惯性坐标系中的位置和速度"""
# 1. 计算轨道平面内的位置和速度
# 半正焦弦 p
p = self.a * (1 - self.e**2)
# 轨道平面内的极坐标位置
r = p / (1 + self.e * cos(self.nu))
# 位置在轨道平面内的分量
r_perifocal = np.array([r * cos(self.nu), r * sin(self.nu), 0.0])
# 速度在轨道平面内的分量
v_perifocal = np.array([
-sin(self.nu) * sqrt(self.mu / p),
(self.e + cos(self.nu)) * sqrt(self.mu / p),
0.0
])
# 2. 转换矩阵：从轨道平面到地心惯性坐标系
# 绕 Z 轴旋转 (-raan)
R3_raan = np.array([
[cos(-self.raan), -sin(-self.raan), 0],
[sin(-self.raan), cos(-self.raan), 0],
[0, 0, 1]
])
# 绕 X 轴旋转 (-i)
R1_i = np.array([
[1, 0, 0],
[0, cos(-self.i), -sin(-self.i)],
[0, sin(-self.i), cos(-self.i)]
])
# 绕 Z 轴旋转 (-arg_peri)
R3_arg_peri = np.array([
[cos(-self.arg_peri), -sin(-self.arg_peri), 0],
[sin(-self.arg_peri), cos(-self.arg_peri), 0],
[0, 0, 1]
])
# 完整的转换矩阵 Q
Q = R3_raan @ R1_i @ R3_arg_peri
# 3. 转换到地心惯性坐标系
r_eci = Q @ r_perifocal
v_eci = Q @ v_perifocal
return r_eci, v_eci
def calculate_orbital_period(self):
"""计算轨道周期"""
return 2 * np.pi * sqrt(self.a**3 / self.mu)
def propagate_orbit(self, delta_t):
"""传播轨道到新时间点
参数:
delta_t: 时间增量 (秒)
"""
# 计算平近点角变化
n = sqrt(self.mu / self.a**3)
# 平均运动
M0 = self.eccentric_to_mean(self.true_to_eccentric(self.nu))
M = M0 + n * delta_t
# 解开普勒方程求偏近点角
E = self.solve_kepler(M, self.e)
# 计算新的真近点角
nu_new = self.eccentric_to_true(E)
# 更新轨道要素
return KeplerOrbit(
self.a, self.e, degrees(self.i), degrees(self.raan),
degrees(self.arg_peri), degrees(nu_new), self.mu
)
@staticmethod
def true_to_eccentric(nu, e):
"""真近点角转偏近点角"""
return 2 * np.arctan(np.sqrt((1 - e) / (1 + e)) * np.tan(nu / 2))
@staticmethod
def eccentric_to_true(E, e):
"""偏近点角转真近点角"""
return 2 * np.arctan(np.sqrt((1 + e) / (1 - e)) * np.tan(E / 2))
@staticmethod
def eccentric_to_mean(E, e):
"""偏近点角转平近点角"""
return E - e * np.sin(E)
@staticmethod
def solve_kepler(M, e, tol=1e-12, max_iter=100):
"""解开普勒方程 E - e*sin(E) = M
参数:
M: 平近点角
e: 偏心率
tol: 容差
max_iter: 最大迭代次数
返回:
E: 偏近点角
"""
# 初始估计
if e < 0.8:
E = M
else:
E = np.pi
# 牛顿 - 拉弗森迭代
for i in range(max_iter):
f = E - e * np.sin(E) - M
f_prime = 1 - e * np.cos(E)
delta = f / f_prime
E -= delta
if abs(delta) < tol:
break
return E
5. 轨道摄动模型
实际轨道受多种摄动影响，主要考虑：
5.1 J2 摄动（地球扁率）
地球非球形导致的主要摄动项：
U_J2 = (mu * J2 * Re^2) / (2 * r^3) * (3 * sin^2(phi) - 1)
其中 J2 ≈ 1.08263×10^-3，Re 为地球赤道半径。
class PerturbationModel:
"""轨道摄动模型"""
def __init__(self, J2=1.08263e-3, R_e=6378137.0):
self.J2 = J2
self.R_e = R_e
def J2_acceleration(self, r_eci):
"""计算 J2 摄动加速度
参数:
r_eci: 卫星在地心惯性坐标系中的位置向量 (m)
返回:
a_J2: J2 摄动加速度向量 (m/s^2)
"""
x, y, z = r_eci
r = np.linalg.norm(r_eci)
# 地球引力常数
mu = 3.986004418e14
# J2 摄动加速度分量
factor = 1.5 * self.J2 * mu * self.R_e**2 / r**5
a_x = factor * x * (5 * z**2 / r**2 - 1)
a_y = factor * y * (5 * z**2 / r**2 - 1)
a_z = factor * z * (5 * z**2 / r**2 - 3)
return np.array([a_x, a_y, a_z])
def atmospheric_drag(self, r_eci, v_eci, satellite_params):
"""大气阻力摄动
参数:
r_eci: 位置向量 (m)
v_eci: 速度向量 (m/s)
satellite_params: 卫星参数字典
返回:
a_drag: 大气阻力加速度 (m/s^2)
"""
# 提取卫星参数
Cd = satellite_params.get('drag_coefficient', 2.2)
A = satellite_params.get('cross_sectional_area', 1.0) # m^2
m = satellite_params.get('mass', 100.0) # kg
# 计算高度（简化模型）
r = np.linalg.norm(r_eci)
h = r - self.R_e
# 大气密度模型（指数模型）
rho0 = 1.225 # 海平面密度 kg/m^3
H = 8500.0 # 标高 m
rho = rho0 * np.exp(-h / H)
# 相对速度（忽略大气旋转）
v_rel = np.linalg.norm(v_eci)
# 阻力加速度
a_drag_magnitude = -0.5 * Cd * (A/m) * rho * v_rel**2
# 阻力方向与速度方向相反
if v_rel > 0:
a_drag = a_drag_magnitude * (-v_eci / v_rel)
else:
a_drag = np.zeros(3)
return a_drag
def solar_radiation_pressure(self, r_eci, satellite_params):
"""太阳辐射压力
参数:
r_eci: 位置向量 (m)
satellite_params: 卫星参数字典
返回:
a_srp: 太阳辐射压力加速度 (m/s^2)
"""
# 太阳辐射常数
P_sun = 4.56e-6 # N/m^2
# 卫星参数
Cr = satellite_params.get('reflectivity_coefficient', 1.3)
A = satellite_params.get('cross_sectional_area', 1.0) # m^2
m = satellite_params.get('mass', 100.0) # kg
# 简化：假设太阳在 X 方向（春分点方向）
sun_direction = np.array([1, 0, 0])
# 卫星指向太阳的方向
sat_to_sun = sun_direction - r_eci / np.linalg.norm(r_eci)
sat_to_sun = sat_to_sun / np.linalg.norm(sat_to_sun)
# 辐射压力加速度
a_srp_magnitude = -Cr * (A/m) * P_sun
return a_srp_magnitude * sat_to_sun
5.2 数值积分器
对于包含摄动的精确轨道计算，需要使用数值积分：
class NumericalIntegrator:
"""轨道数值积分器"""
def __init__(self, perturbation_model=None):
self.perturbation_model = perturbation_model or PerturbationModel()
self.mu = 3.986004418e14
def ode_equations(self, t, y, satellite_params):
"""轨道运动微分方程
参数:
t: 时间 (s)
y: 状态向量 [x, y, z, vx, vy, vz]
satellite_params: 卫星参数
返回:
dydt: 状态向量导数
"""
# 位置和速度
r = y[:3]
v = y[3:]
# 二体问题加速度
r_norm = np.linalg.norm(r)
a_two_body = -self.mu * r / r_norm**3
# 总加速度（初始化为二体加速度）
a_total = a_two_body
# 添加摄动加速度
if self.perturbation_model:
# J2 摄动
a_J2 = self.perturbation_model.J2_acceleration(r)
a_total += a_J2
# 大气阻力（仅适用于低轨）
r_norm_km = r_norm / 1000
if r_norm_km < 2000: # 高度小于 2000km
a_drag = self.perturbation_model.atmospheric_drag(r, v, satellite_params)
a_total += a_drag
# 太阳辐射压力
a_srp = self.perturbation_model.solar_radiation_pressure(r, satellite_params)
a_total += a_srp
# 状态导数
dydt = np.zeros(6)
dydt[:3] = v # 位置导数 = 速度
dydt[3:] = a_total # 速度导数 = 加速度
return dydt
def integrate_orbit(self, initial_state, t_span, dt, satellite_params, method='rk4'):
"""积分轨道
参数:
initial_state: 初始状态 [x, y, z, vx, vy, vz] (m, m/s)
t_span: 时间范围 [t0, tf] (s)
dt: 时间步长 (s)
satellite_params: 卫星参数
method: 积分方法 ('rk4' 或 'adams')
返回:
t_array: 时间数组
state_history: 状态历史
"""
t0, tf = t_span
n_steps = int((tf - t0) / dt) + 1
# 初始化数组
t_array = np.linspace(t0, tf, n_steps)
state_history = np.zeros((n_steps, 6))
state_history[0] = initial_state
if method == 'rk4':
# 龙格 - 库塔 4 阶积分
for i in range(n_steps - 1):
t = t_array[i]
y = state_history[i]
k1 = self.ode_equations(t, y, satellite_params)
k2 = self.ode_equations(t + dt/2, y + dt/2 * k1, satellite_params)
k3 = self.ode_equations(t + dt/2, y + dt/2 * k2, satellite_params)
k4 = self.ode_equations(t + dt, y + dt * k3, satellite_params)
state_history[i+1] = y + dt/6 * (k1 + 2*k2 + 2*k3 + k4)
elif method == 'adams':
# 亚当斯 - 巴什福斯方法（4 阶）
# 使用 RK4 初始化前 4 步
for i in range(3):
t = t_array[i]
y = state_history[i]
k1 = self.ode_equations(t, y, satellite_params)
k2 = self.ode_equations(t + dt/2, y + dt/2 * k1, satellite_params)
k3 = self.ode_equations(t + dt/2, y + dt/2 * k2, satellite_params)
k4 = self.ode_equations(t + dt, y + dt * k3, satellite_params)
state_history[i+1] = y + dt/6 * (k1 + 2*k2 + 2*k3 + k4)
# 亚当斯 - 巴什福斯主循环
for i in range(3, n_steps - 1):
t = t_array[i]
y = state_history[i]
# 计算过去 4 个点的导数
f_0 = self.ode_equations(t_array[i], state_history[i], satellite_params)
f_1 = self.ode_equations(t_array[i-1], state_history[i-1], satellite_params)
f_2 = self.ode_equations(t_array[i-2], state_history[i-2], satellite_params)
f_3 = self.ode_equations(t_array[i-3], state_history[i-3], satellite_params)
# 4 阶亚当斯 - 巴什福斯公式
state_history[i+1] = y + dt/24 * (55*f_0 - 59*f_1 + 37*f_2 - 9*f_3)
return t_array, state_history
6. 低轨星座构型设计
6.1 Walker 星座
Walker 星座是最常用的星座构型，由三个参数描述：
- 总卫星数 (T)
- 轨道面数 (P)
- 相位因子 (F)
class WalkerConstellation:
"""Walker 星座设计类"""
def __init__(self, T, P, F, altitude, inclination, raan0=0):
"""
初始化 Walker 星座
参数:
T: 卫星总数
P: 轨道面数
F: 相位因子 (0 到 P-1 之间的整数)
altitude: 轨道高度 (km)
inclination: 轨道倾角 (度)
raan0: 第一个轨道面的升交点赤经 (度)
"""
self.T = T
self.P = P
self.F = F
self.altitude = altitude # km
self.inclination = inclination # 度
self.raan0 = raan0 # 度
# 计算每个轨道面的卫星数
self.S = T // P # 每个轨道面的卫星数
# 计算轨道半长轴
R_e = 6371.0 # 地球半径 (km)
self.a = (R_e + altitude) * 1000 # 转换为米
# 生成星座中的所有卫星
self.satellites = self.generate_constellation()
def generate_constellation(self):
"""生成 Walker 星座中的所有卫星"""
satellites = []
# 每个轨道面
for p in range(self.P):
# 轨道面升交点赤经
raan = self.raan0 + p * (360.0 / self.P)
# 每个轨道面内的卫星
for s in range(self.S):
# 轨道面内卫星的相位
phase = (360.0 / self.T) * (s * self.P + p * self.F)
# 创建卫星轨道要素
# 假设圆轨道 (e=0)
# 真近点角 = 相位角
nu = phase
# 创建轨道对象
orbit = KeplerOrbit(
a=self.a, e=0.0, i=self.inclination,
raan=raan, arg_peri=0.0, # 圆轨道，近地点幅角为 0
nu=nu
)
satellites.append({
'plane': p,
'position_in_plane': s,
'phase': phase,
'orbit': orbit
})
return satellites
def get_satellite_position_velocity(self, satellite_idx, time=0):
"""获取指定卫星在给定时间的位置和速度
参数:
satellite_idx: 卫星索引
time: 时间 (秒)
返回:
position: 位置向量 (m)
velocity: 速度向量 (m/s)
"""
satellite = self.satellites[satellite_idx]
# 如果时间不为 0，传播轨道
if time != 0:
orbit = satellite['orbit'].propagate_orbit(time)
else:
orbit = satellite['orbit']
# 计算位置和速度
return orbit.calculate_position_velocity()
def calculate_coverage_statistics(self, ground_points, min_elevation=10):
"""计算星座对地面点的覆盖统计
参数:
ground_points: 地面点列表 [(lat, lon), ...] (度)
min_elevation: 最小仰角 (度)
返回:
coverage_stats: 覆盖统计字典
"""
coverage_stats = {
'total_points': len(ground_points),
'min_elevation': min_elevation,
'coverage_percentage': [],
'average_gap': [],
'max_gap': []
}
for lat, lon in ground_points:
# 模拟 24 小时覆盖
time_step = 60 # 1 分钟
total_time = 24 * 3600 # 24 小时
n_steps = total_time // time_step
coverage_mask = np.zeros(n_steps, dtype=bool)
for step in range(n_steps):
t = step * time_step
# 检查是否有卫星可见
visible = False
for sat in self.satellites:
# 传播轨道到时间 t
orbit = sat['orbit'].propagate_orbit(t)
r_eci, _ = orbit.calculate_position_velocity()
# 计算仰角
elevation = self.calculate_elevation(r_eci, lat, lon, t)
if elevation >= min_elevation:
visible = True
break
coverage_mask[step] = visible
# 计算统计量
coverage_percentage = np.sum(coverage_mask) / n_steps * 100
# 计算覆盖间隔
gaps = self.find_coverage_gaps(coverage_mask, time_step)
if gaps:
avg_gap = np.mean(gaps)
max_gap = np.max(gaps)
else:
avg_gap = 0
max_gap = 0
coverage_stats['coverage_percentage'].append(coverage_percentage)
coverage_stats['average_gap'].append(avg_gap)
coverage_stats['max_gap'].append(max_gap)
return coverage_stats
@staticmethod
def calculate_elevation(r_eci, lat, lon, t):
"""计算卫星相对于地面点的仰角
参数:
r_eci: 卫星位置 (ECI 坐标系，m)
lat: 地面点纬度 (度)
lon: 地面点经度 (度)
t: 时间 (秒，用于地球旋转)
返回:
elevation: 仰角 (度)
"""
# 地球自转角速度 (rad/s)
omega_e = 7.2921150e-5
# 将地面点转换到 ECI 坐标系
lat_rad = np.radians(lat)
lon_rad = np.radians(lon)
# 地球半径
R_e = 6378137.0 # 赤道半径 (m)
# 考虑地球旋转的地面点经度
lon_rad_rotated = lon_rad + omega_e * t
# 地面点 ECI 坐标
r_ground_eci = np.array([
R_e * np.cos(lat_rad) * np.cos(lon_rad_rotated),
R_e * np.cos(lat_rad) * np.sin(lon_rad_rotated),
R_e * np.sin(lat_rad)
])
# 卫星相对于地面点的向量
r_rel = r_eci - r_ground_eci
# 地面点处的单位法向量
n_ground = r_ground_eci / np.linalg.norm(r_ground_eci)
# 计算仰角
r_rel_norm = np.linalg.norm(r_rel)
sin_elevation = np.dot(r_rel, n_ground) / r_rel_norm
elevation = np.degrees(np.arcsin(sin_elevation))
return elevation
@staticmethod
def find_coverage_gaps(coverage_mask, time_step):
"""找出覆盖间隔
参数:
coverage_mask: 布尔数组，True 表示有覆盖
time_step: 时间步长 (秒)
返回:
gaps: 覆盖间隔列表 (秒)
"""
gaps = []
current_gap = 0
for covered in coverage_mask:
if not covered:
current_gap += time_step
else:
if current_gap > 0:
gaps.append(current_gap)
current_gap = 0
# 最后一个间隔
if current_gap > 0:
gaps.append(current_gap)
return gaps
7. Python 实现：轨道计算库
7.1 完整轨道计算框架
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
from datetime import datetime, timedelta
class SatelliteOrbitSimulator:
"""卫星轨道模拟器"""
def __init__(self, constellation_config=None):
"""
初始化轨道模拟器
参数:
constellation_config: 星座配置字典
"""
self.constellation_config = constellation_config
self.constellation = None
self.integrator = NumericalIntegrator()
# 地球参数
self.R_e = 6378137.0 # 赤道半径 (m)
self.mu = 3.986004418e14 # 地球引力常数
# 时间系统
self.jd2000 = 2451545.0 # 2000 年 1 月 1 日 12:00 UT 的儒略日
def create_walker_constellation(self, T, P, F, altitude, inclination):
"""创建 Walker 星座"""
self.constellation = WalkerConstellation(T, P, F, altitude, inclination)
return self.constellation
def simulate_constellation(self, duration_hours=24, time_step=60):
"""模拟星座运行
参数:
duration_hours: 模拟时长 (小时)
time_step: 时间步长 (秒)
返回:
simulation_data: 模拟数据
"""
if self.constellation is None:
raise ValueError("未创建星座，请先调用 create_walker_constellation")
total_time = duration_hours * 3600 # 转换为秒
n_steps = int(total_time / time_step) + 1
# 初始化数据存储
simulation_data = {
'time': np.linspace(0, total_time, n_steps),
'satellite_positions': [],
'satellite_velocities': [],
'ground_track': []
}
# 为每颗卫星分配存储
n_satellites = len(self.constellation.satellites)
for i in range(n_satellites):
simulation_data['satellite_positions'].append(np.zeros((n_steps, 3)))
simulation_data['satellite_velocities'].append(np.zeros((n_steps, 3)))
simulation_data['ground_track'].append(np.zeros((n_steps, 2))) # 经纬度
# 时间循环
for step_idx, t in enumerate(simulation_data['time']):
for sat_idx in range(n_satellites):
# 获取卫星位置和速度
r_eci, v_eci = self.constellation.get_satellite_position_velocity(sat_idx, t)
# 存储数据
simulation_data['satellite_positions'][sat_idx][step_idx] = r_eci
simulation_data['satellite_velocities'][sat_idx][step_idx] = v_eci
# 计算并存储星下点
lat, lon = self.eci_to_geodetic(r_eci, t)
simulation_data['ground_track'][sat_idx][step_idx] = [lat, lon]
return simulation_data
def eci_to_geodetic(self, r_eci, t):
"""将 ECI 坐标转换为地理坐标
参数:
r_eci: ECI 坐标 (m)
t: 时间 (秒)
返回:
lat: 纬度 (度)
lon: 经度 (度)
"""
x, y, z = r_eci
# 地球自转角度
omega_e = 7.2921150e-5 # rad/s
theta = omega_e * t
# 旋转到地固坐标系
x_rot = x * np.cos(theta) + y * np.sin(theta)
y_rot = -x * np.sin(theta) + y * np.cos(theta)
z_rot = z
# 计算经度
lon = np.degrees(np.arctan2(y_rot, x_rot))
# 计算地心纬度
r = np.sqrt(x_rot**2 + y_rot**2 + z_rot**2)
lat_geocentric = np.degrees(np.arcsin(z_rot / r))
# 简化：假设地球为球体
lat = lat_geocentric
return lat, lon
def calculate_coverage(self, ground_stations, duration_hours=24, min_elevation=10):
"""计算对地面站的覆盖
参数:
ground_stations: 地面站列表 [(lat, lon, name), ...]
duration_hours: 模拟时长 (小时)
min_elevation: 最小仰角 (度)
返回:
coverage_data: 覆盖数据
"""
if self.constellation is None:
raise ValueError("未创建星座")
total_time = duration_hours * 3600
time_step = 60 # 1 分钟
n_steps = int(total_time / time_step)
coverage_data = {
'ground_stations': ground_stations,
'time': np.linspace(0, total_time, n_steps),
'visibility': np.zeros((len(ground_stations), n_steps)),
'elevation': np.zeros((len(ground_stations), len(self.constellation.satellites), n_steps))
}
# 计算每个时间步的可见性
for t_idx, t in enumerate(coverage_data['time']):
for gs_idx, (lat, lon, _) in enumerate(ground_stations):
best_elevation = -90
for sat_idx, sat in enumerate(self.constellation.satellites):
# 获取卫星位置
r_eci, _ = self.constellation.get_satellite_position_velocity(sat_idx, t)
# 计算仰角
elevation = self.constellation.calculate_elevation(r_eci, lat, lon, t)
coverage_data['elevation'][gs_idx, sat_idx, t_idx] = elevation
# 更新最佳仰角
if elevation > best_elevation:
best_elevation = elevation
# 检查可见性
if best_elevation >= min_elevation:
coverage_data['visibility'][gs_idx, t_idx] = 1
return coverage_data
7.2 高级轨道分析工具
class OrbitAnalyzer:
"""轨道分析工具"""
@staticmethod
def calculate_orbit_elements(r_eci, v_eci, mu=3.986004418e14):
"""从位置和速度计算轨道要素
参数:
r_eci: 位置向量 (m)
v_eci: 速度向量 (m/s)
mu: 引力常数
返回:
elements: 轨道要素字典
"""
# 位置和速度的模
r = np.linalg.norm(r_eci)
v = np.linalg.norm(v_eci)
# 角动量向量
h = np.cross(r_eci, v_eci)
h_norm = np.linalg.norm(h)
# 节点向量
n = np.cross([0, 0, 1], h)
n_norm = np.linalg.norm(n)
# 偏心率向量
e_vec = ((v**2 - mu/r) * r_eci - np.dot(r_eci, v_eci) * v_eci) / mu
e = np.linalg.norm(e_vec)
# 轨道能量
energy = v**2/2 - mu/r
# 半长轴
if e != 1:
a = -mu/(2*energy)
else:
a = float('inf')
# 轨道倾角
i = np.arccos(h[2]/h_norm)
# 升交点赤经
if n_norm != 0:
raan = np.arccos(n[0]/n_norm)
if n[1] < 0:
raan = 2*np.pi - raan
else:
raan = 0
# 近地点幅角
if n_norm != 0 and e > 0:
arg_peri = np.arccos(np.dot(n, e_vec)/(n_norm*e))
if e_vec[2] < 0:
arg_peri = 2*np.pi - arg_peri
else:
arg_peri = 0
# 真近点角
if e > 0:
nu = np.arccos(np.dot(e_vec, r_eci)/(e*r))
if np.dot(r_eci, v_eci) < 0:
nu = 2*np.pi - nu
else:
nu = 0
# 平近点角
if e < 1:
E = 2 * np.arctan(np.sqrt((1-e)/(1+e)) * np.tan(nu/2))
M = E - e*np.sin(E)
else:
M = 0
return {
'a': a, # 半长轴 (m)
'e': e, # 偏心率
'i': np.degrees(i), # 倾角 (度)
'raan': np.degrees(raan), # 升交点赤经 (度)
'arg_peri': np.degrees(arg_peri), # 近地点幅角 (度)
'nu': np.degrees(nu), # 真近点角 (度)
'M': M, # 平近点角 (rad)
'period': 2*np.pi*np.sqrt(a**3/mu) if a>0 else float('inf')
}
@staticmethod
def analyze_constellation_coverage(coverage_data):
"""分析星座覆盖性能
参数:
coverage_data: 覆盖数据
返回:
analysis: 分析结果
"""
n_stations = len(coverage_data['ground_stations'])
n_time_steps = len(coverage_data['time'])
analysis = {
'station_coverage': {},
'global_metrics': {}
}
# 每个地面站的覆盖统计
for gs_idx, (lat, lon, name) in enumerate(coverage_data['ground_stations']):
visibility = coverage_data['visibility'][gs_idx]
coverage_percentage = np.sum(visibility) / n_time_steps * 100
# 计算覆盖间隔
gaps = []
current_gap = 0
for v in visibility:
if v == 0:
current_gap += 1
else:
if current_gap > 0:
gaps.append(current_gap)
current_gap = 0
if current_gap > 0:
gaps.append(current_gap)
avg_gap = np.mean(gaps) if gaps else 0
max_gap = np.max(gaps) if gaps else 0
analysis['station_coverage'][name] = {
'latitude': lat,
'longitude': lon,
'coverage_percentage': coverage_percentage,
'average_gap_minutes': avg_gap,
'max_gap_minutes': max_gap,
'gaps': gaps
}
# 全局指标
all_visibility = coverage_data['visibility']
global_coverage = np.mean(all_visibility) * 100
# 计算同时覆盖的地面站数量
simultaneous_coverage = np.sum(all_visibility, axis=0)
analysis['global_metrics'] = {
'global_coverage_percentage': global_coverage,
'min_simultaneous_coverage': np.min(simultaneous_coverage),
'max_simultaneous_coverage': np.max(simultaneous_coverage),
'avg_simultaneous_coverage': np.mean(simultaneous_coverage)
}
return analysis
8. 轨道可视化与仿真
8.1 3D 可视化
class OrbitVisualizer:
"""轨道可视化工具"""
def __init__(self):
self.fig = None
self.ax = None
def plot_3d_orbit(self, simulation_data, title="卫星轨道", show_earth=True):
"""绘制 3D 轨道图"""
self.fig = plt.figure(figsize=(12, 10))
self.ax = self.fig.add_subplot(111, projection='3d')
# 绘制地球
if show_earth:
self._plot_earth()
# 绘制每颗卫星的轨道
colors = plt.cm.rainbow(np.linspace(0, 1, len(simulation_data['satellite_positions'])))
for idx, positions in enumerate(simulation_data['satellite_positions']):
# 提取坐标
x = positions[:, 0] / 1000 # 转换为 km
y = positions[:, 1] / 1000
z = positions[:, 2] / 1000
# 绘制轨道
self.ax.plot(x, y, z, color=colors[idx], linewidth=1.0, alpha=0.7, label=f'Sat {idx+1}')
# 设置图形属性
self.ax.set_xlabel('X (km)')
self.ax.set_ylabel('Y (km)')
self.ax.set_zlabel('Z (km)')
self.ax.set_title(title)
# 设置坐标轴比例
max_range = max([np.ptp(x) for x in [self.ax.get_xlim(), self.ax.get_ylim(), self.ax.get_zlim()]] ) / 2.0
mid_x = (self.ax.get_xlim()[0] + self.ax.get_xlim()[1]) * 0.5
mid_y = (self.ax.get_ylim()[0] + self.ax.get_ylim()[1]) * 0.5
mid_z = (self.ax.get_zlim()[0] + self.ax.get_zlim()[1]) * 0.5
self.ax.set_xlim(mid_x - max_range, mid_x + max_range)
self.ax.set_ylim(mid_y - max_range, mid_y + max_range)
self.ax.set_zlim(mid_z - max_range, mid_z + max_range)
plt.legend(loc='upper right', bbox_to_anchor=(1.15, 1))
plt.tight_layout()
return self.fig, self.ax
def _plot_earth(self):
"""绘制地球"""
# 地球半径 (km)
R_e = 6378.137
# 创建球体
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = R_e * np.outer(np.cos(u), np.sin(v))
y = R_e * np.outer(np.sin(u), np.sin(v))
z = R_e * np.outer(np.ones(np.size(u)), np.cos(v))
# 绘制地球
self.ax.plot_surface(x, y, z, color='lightblue', alpha=0.3, edgecolor='none')
# 绘制赤道
theta = np.linspace(0, 2*np.pi, 100)
x_eq = R_e * np.cos(theta)
y_eq = R_e * np.sin(theta)
self.ax.plot(x_eq, y_eq, 0, 'b--', alpha=0.5, linewidth=0.5)
def plot_ground_tracks(self, simulation_data, title="星下点轨迹"):
"""绘制星下点轨迹"""
fig, axes = plt.subplots(1, 2, figsize=(16, 8))
# 设置地图背景
for ax in axes:
ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.set_xlabel('经度 (度)')
ax.set_ylabel('纬度 (度)')
ax.grid(True, alpha=0.3)
# 添加经纬度网格
ax.set_xticks(np.arange(-180, 181, 30))
ax.set_yticks(np.arange(-90, 91, 30))
# 绘制每颗卫星的星下点
colors = plt.cm.rainbow(np.linspace(0, 1, len(simulation_data['ground_track'])))
for idx, ground_track in enumerate(simulation_data['ground_track']):
lons = ground_track[:, 1]
lats = ground_track[:, 0]
# 处理经度不连续
lons_wrapped = np.where(lons > 180, lons - 360, lons)
# 绘制轨迹
axes[0].plot(lons_wrapped, lats, color=colors[idx], linewidth=0.5, alpha=0.7)
# 绘制散点图
axes[1].scatter(lons_wrapped[::10], lats[::10], # 每隔 10 个点取一个
color=colors[idx], s=1, alpha=0.5)
axes[0].set_title('星下点轨迹 - 连线图')
axes[1].set_title('星下点轨迹 - 散点图')
plt.suptitle(title)
plt.tight_layout()
return fig, axes
def plot_coverage_analysis(self, coverage_data, analysis_results, station_name=None):
"""绘制覆盖分析图"""
if station_name:
# 绘制单个地面站的覆盖
for name, data in analysis_results['station_coverage'].items():
if name == station_name:
station_data = data
break
else:
raise ValueError(f"未找到地面站：{station_name}")
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
# 1. 可见性时间序列
ax = axes[0, 0]
for gs_idx, (lat, lon, name_gs) in enumerate(coverage_data['ground_stations']):
if name_gs == station_name:
visibility = coverage_data['visibility'][gs_idx]
time_hours = coverage_data['time'] / 3600
ax.plot(time_hours, visibility, 'b-', linewidth=1)
ax.fill_between(time_hours, 0, visibility, alpha=0.3)
ax.set_xlabel('时间 (小时)')
ax.set_ylabel('可见性')
ax.set_title(f'{station_name} 可见性 ({station_data["coverage_percentage"]:.1f}%)')
ax.set_ylim(-0.1, 1.1)
ax.grid(True, alpha=0.3)
break
# 2. 仰角时间序列
ax = axes[0, 1]
for gs_idx, (lat, lon, name_gs) in enumerate(coverage_data['ground_stations']):
if name_gs == station_name:
n_satellites = coverage_data['elevation'].shape[1]
time_hours = coverage_data['time'] / 3600
for sat_idx in range(n_satellites):
elevation = coverage_data['elevation'][gs_idx, sat_idx, :]
ax.plot(time_hours, elevation, linewidth=0.5, alpha=0.5)
ax.set_xlabel('时间 (小时)')
ax.set_ylabel('仰角 (度)')
ax.set_title('各卫星仰角')
ax.grid(True, alpha=0.3)
break
# 3. 覆盖间隔直方图
ax = axes[1, 0]
if station_data['gaps']:
gaps_minutes = np.array(station_data['gaps'])
ax.hist(gaps_minutes, bins=20, edgecolor='black', alpha=0.7)
ax.set_xlabel('覆盖间隔 (分钟)')
ax.set_ylabel('频次')
ax.set_title(f'覆盖间隔分布 (平均：{station_data["average_gap_minutes"]:.1f} 分钟)')
ax.grid(True, alpha=0.3)
# 4. 全天仰角分布
ax = axes[1, 1]
for gs_idx, (lat, lon, name_gs) in enumerate(coverage_data['ground_stations']):
if name_gs == station_name:
all_elevations = coverage_data['elevation'][gs_idx, :, :].flatten()
all_elevations = all_elevations[all_elevations > -90]
ax.hist(all_elevations, bins=30, edgecolor='black', alpha=0.7)
ax.set_xlabel('仰角 (度)')
ax.set_ylabel('频次')
ax.set_title('全天仰角分布')
ax.grid(True, alpha=0.3)
break
plt.suptitle(f'地面站覆盖分析：{station_name}', fontsize=16)
plt.tight_layout()
else:
# 绘制所有地面站的覆盖统计
fig, axes = plt.subplots(2, 2, figsize=(15, 10))
# 1. 覆盖百分比条形图
ax = axes[0, 0]
station_names = []
coverage_percentages = []
for name, data in analysis_results['station_coverage'].items():
station_names.append(name)
coverage_percentages.append(data['coverage_percentage'])
y_pos = np.arange(len(station_names))
bars = ax.barh(y_pos, coverage_percentages, alpha=0.7)
ax.set_yticks(y_pos)
ax.set_yticklabels(station_names)
ax.set_xlabel('覆盖百分比 (%)')
ax.set_title('各地面站覆盖百分比')
ax.grid(True, alpha=0.3, axis='x')
# 2. 同时覆盖地面站数量
ax = axes[0, 1]
time_hours = coverage_data['time'] / 3600
all_visibility = coverage_data['visibility']
simultaneous_coverage = np.sum(all_visibility, axis=0)
ax.plot(time_hours, simultaneous_coverage, 'b-', linewidth=1)
ax.fill_between(time_hours, 0, simultaneous_coverage, alpha=0.3)
ax.set_xlabel('时间 (小时)')
ax.set_ylabel('同时覆盖地面站数量')
ax.set_title(f'同时覆盖地面站数量 (平均：{analysis_results["global_metrics"]["avg_simultaneous_coverage"]:.1f})')
ax.grid(True, alpha=0.3)
# 3. 全球覆盖统计
ax = axes[1, 0]
metrics = analysis_results['global_metrics']
metric_names = ['全球覆盖', '最小同时覆盖', '最大同时覆盖', '平均同时覆盖']
metric_values = [
metrics['global_coverage_percentage'],
metrics['min_simultaneous_coverage'],
metrics['max_simultaneous_coverage'],
metrics['avg_simultaneous_coverage']
]
x_pos = np.arange(len(metric_names))
bars = ax.bar(x_pos, metric_values, alpha=0.7)
ax.set_xticks(x_pos)
ax.set_xticklabels(metric_names, rotation=45, ha='right')
ax.set_ylabel('数值')
ax.set_title('全球覆盖统计')
ax.grid(True, alpha=0.3, axis='y')
# 添加数值标签
for bar, value in zip(bars, metric_values):
height = bar.get_height()
ax.text(bar.get_x() + bar.get_width()/2., height, f'{value:.1f}', ha='center', va='bottom')
# 4. 经纬度覆盖热图
ax = axes[1, 1]
# 创建经纬度网格
lats = []
lons = []
coverage_values = []
for name, data in analysis_results['station_coverage'].items():
lats.append(data['latitude'])
lons.append(data['longitude'])
coverage_values.append(data['coverage_percentage'])
# 简单的散点图表示
scatter = ax.scatter(lons, lats, c=coverage_values, cmap='RdYlGn', s=100, edgecolor='black')
ax.set_xlabel('经度 (度)')
ax.set_ylabel('纬度 (度)')
ax.set_title('地面站覆盖百分比分布')
ax.grid(True, alpha=0.3)
# 添加颜色条
plt.colorbar(scatter, ax=ax, label='覆盖百分比 (%)')
plt.suptitle('全球覆盖分析', fontsize=16)
plt.tight_layout()
return fig, axes
8.2 动画生成
class OrbitAnimator:
"""轨道动画生成器"""
def __init__(self, simulation_data):
self.simulation_data = simulation_data
self.fig = None
self.ax = None
def create_3d_animation(self, output_file='orbit_animation.mp4', fps=30, dpi=100):
"""创建 3D 轨道动画"""
from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt
# 设置图形
self.fig = plt.figure(figsize=(12, 10))
self.ax = self.fig.add_subplot(111, projection='3d')
# 绘制地球
self._plot_earth()
# 初始化卫星点
n_satellites = len(self.simulation_data['satellite_positions'])
satellites_points = []
satellites_lines = []
colors = plt.cm.rainbow(np.linspace(0, 1, n_satellites))
for i in range(n_satellites):
# 卫星当前位置点
point, = self.ax.plot([], [], [], 'o', color=colors[i], markersize=6, label=f'Sat {i+1}')
satellites_points.append(point)
# 卫星轨道线（最近 N 个点）
line, = self.ax.plot([], [], [], color=colors[i], linewidth=1.0, alpha=0.5)
satellites_lines.append(line)
# 设置图形属性
self.ax.set_xlabel('X (km)')
self.ax.set_ylabel('Y (km)')
self.ax.set_zlabel('Z (km)')
self.ax.set_title('卫星轨道动画')
# 设置坐标轴范围
all_positions = np.vstack(self.simulation_data['satellite_positions'])
max_range = np.max(np.abs(all_positions)) / 1000 * 1.1 # 转换为 km
self.ax.set_xlim(-max_range, max_range)
self.ax.set_ylim(-max_range, max_range)
self.ax.set_zlim(-max_range, max_range)
# 计算总帧数
n_frames = len(self.simulation_data['time'])
frame_step = max(1, n_frames // (fps * 60)) # 目标：60 秒动画
def update(frame):
"""更新动画帧"""
idx = frame * frame_step
if idx >= n_frames:
return satellites_points + satellites_lines
# 更新每颗卫星
for i in range(n_satellites):
positions = self.simulation_data['satellite_positions'][i]
# 当前点位置
x = positions[idx, 0] / 1000
y = positions[idx, 1] / 1000
z = positions[idx, 2] / 1000
satellites_points[i].set_data([x], [y])
satellites_points[i].set_3d_properties([z])
# 轨道线（最近 100 个点）
start_idx = max(0, idx - 100)
x_line = positions[start_idx:idx+1, 0] / 1000
y_line = positions[start_idx:idx+1, 1] / 1000
z_line = positions[start_idx:idx+1, 2] / 1000
satellites_lines[i].set_data(x_line, y_line)
satellites_lines[i].set_3d_properties(z_line)
# 更新时间标题
time_seconds = self.simulation_data['time'][idx]
hours = int(time_seconds // 3600)
minutes = int((time_seconds % 3600) // 60)
self.ax.set_title(f'卫星轨道动画 - 时间：{hours:02d}:{minutes:02d}')
return satellites_points + satellites_lines
# 创建动画
anim = FuncAnimation(self.fig, update, frames=n_frames//frame_step, interval=1000/fps, blit=False)
# 保存动画
print(f"正在生成动画，请稍候...")
anim.save(output_file, writer='ffmpeg', fps=fps, dpi=dpi)
print(f"动画已保存到：{output_file}")
plt.close(self.fig)
return anim
def _plot_earth(self):
"""绘制地球"""
R_e = 6378.137 # 地球半径 (km)
# 创建球体
u = np.linspace(0, 2 * np.pi, 50)
v = np.linspace(0, np.pi, 50)
x = R_e * np.outer(np.cos(u), np.sin(v))
y = R_e * np.outer(np.sin(u), np.sin(v))
z = R_e * np.outer(np.ones(np.size(u)), np.cos(v))
# 绘制地球
self.ax.plot_surface(x, y, z, color='lightblue', alpha=0.3, edgecolor='none')
9. 卫星可见性与覆盖分析
9.1 可见性预测算法
class VisibilityPredictor:
"""卫星可见性预测器"""
def __init__(self, constellation, min_elevation=10):
self.constellation = constellation
self.min_elevation = min_elevation
self.earth_radius = 6378137.0 # 地球半径 (m)
self.mu = 3.986004418e14 # 地球引力常数
def predict_visibility(self, ground_point, start_time, duration_hours=24, time_step=60):
"""预测卫星可见性
参数:
ground_point: 地面点 (lat, lon) (度)
start_time: 开始时间 (datetime 对象)
duration_hours: 预测时长 (小时)
time_step: 时间步长 (秒)
返回:
visibility_schedule: 可见性时间表
"""
lat, lon = ground_point
lat_rad = np.radians(lat)
lon_rad = np.radians(lon)
# 时间数组
total_time = duration_hours * 3600
time_array = np.arange(0, total_time, time_step)
# 初始化可见性数组
n_satellites = len(self.constellation.satellites)
visibility = np.zeros((n_satellites, len(time_array)))
elevations = np.zeros((n_satellites, len(time_array)))
# 计算每个时间步
for t_idx, t in enumerate(time_array):
# 地球自转角度
omega_e = 7.2921150e-5 # rad/s
theta = omega_e * t
# 地面点 ECI 坐标
x_ground = self.earth_radius * np.cos(lat_rad) * np.cos(lon_rad + theta)
y_ground = self.earth_radius * np.cos(lat_rad) * np.sin(lon_rad + theta)
z_ground = self.earth_radius * np.sin(lat_rad)
ground_pos = np.array([x_ground, y_ground, z_ground])
# 检查每颗卫星
for sat_idx, sat in enumerate(self.constellation.satellites):
# 获取卫星位置
r_eci, _ = self.constellation.get_satellite_position_velocity(sat_idx, t)
# 计算相对向量
r_rel = r_eci - ground_pos
r_rel_norm = np.linalg.norm(r_rel)
# 计算仰角
ground_norm = ground_pos / np.linalg.norm(ground_pos)
sin_elevation = np.dot(r_rel, ground_norm) / r_rel_norm
elevation = np.degrees(np.arcsin(sin_elevation))
elevations[sat_idx, t_idx] = elevation
# 检查可见性
if elevation >= self.min_elevation:
visibility[sat_idx, t_idx] = 1
# 构建可见性时间表
visibility_schedule = {
'time': time_array,
'time_dt': [start_time + timedelta(seconds=float(t)) for t in time_array],
'visibility': visibility,
'elevations': elevations,
'ground_point': ground_point,
'min_elevation': self.min_elevation
}
return visibility_schedule
def analyze_visibility_schedule(self, visibility_schedule):
"""分析可见性时间表"""
visibility = visibility_schedule['visibility']
time_array = visibility_schedule['time']
n_satellites, n_times = visibility.shape
analysis = {
'total_satellites': n_satellites,
'total_time_hours': time_array[-1] / 3600,
'time_step_seconds': time_array[1] - time_array[0] if len(time_array) > 1 else 0,
'satellite_visibility': {},
'summary': {}
}
# 每颗卫星的可见性统计
for sat_idx in range(n_satellites):
vis_times = visibility[sat_idx]
n_visible = np.sum(vis_times)
percentage = n_visible / n_times * 100
# 找到可见时间段
visible_periods = []
in_period = False
start_idx = 0
for t_idx, visible in enumerate(vis_times):
if visible and not in_period:
in_period = True
start_idx = t_idx
elif not visible and in_period:
in_period = False
visible_periods.append({
'start': time_array[start_idx],
'end': time_array[t_idx-1],
'duration': time_array[t_idx-1] - time_array[start_idx]
})
# 处理最后一个时间段
if in_period:
visible_periods.append({
'start': time_array[start_idx],
'end': time_array[-1],
'duration': time_array[-1] - time_array[start_idx]
})
analysis['satellite_visibility'][sat_idx] = {
'visible_count': int(n_visible),
'visible_percentage': percentage,
'visible_periods': visible_periods,
'average_period_duration': np.mean([p['duration'] for p in visible_periods]) if visible_periods else 0,
'max_period_duration': np.max([p['duration'] for p in visible_periods]) if visible_periods else 0
}
# 总体统计
# 任何卫星可见的时间
any_sat_visible = np.any(visibility, axis=0)
any_visible_percentage = np.sum(any_sat_visible) / n_times * 100
# 同时可见的卫星数量
simultaneous = np.sum(visibility, axis=0)
analysis['summary'] = {
'any_satellite_visible_percentage': any_visible_percentage,
'average_simultaneous_satellites': np.mean(simultaneous),
'max_simultaneous_satellites': np.max(simultaneous),
'min_simultaneous_satellites': np.min(simultaneous),
'coverage_gaps': self._find_coverage_gaps(any_sat_visible, time_array)
}
return analysis
def _find_coverage_gaps(self, any_sat_visible, time_array):
"""找出覆盖间隔"""
gaps = []
current_gap_start = None
for i, visible in enumerate(any_sat_visible):
if not visible and current_gap_start is None:
current_gap_start = time_array[i]
elif visible and current_gap_start is not None:
gap_duration = time_array[i] - current_gap_start
gaps.append({
'start': current_gap_start,
'end': time_array[i],
'duration': gap_duration
})
current_gap_start = None
# 处理最后一个间隔
if current_gap_start is not None:
gap_duration = time_array[-1] - current_gap_start
gaps.append({
'start': current_gap_start,
'end': time_array[-1],
'duration': gap_duration
})
return gaps
def get_next_pass(self, visibility_schedule, sat_idx, current_time_idx=0):
"""获取下一次卫星过境"""
visibility = visibility_schedule['visibility'][sat_idx]
time_array = visibility_schedule['time']
# 找到下一个可见时间段
for i in range(current_time_idx, len(visibility)):
if visibility[i]:
start_idx = i
# 找到结束时间
for j in range(start_idx, len(visibility)):
if not visibility[j]:
end_idx = j - 1
break
else:
end_idx = len(visibility) - 1
# 计算过境详情
elevations = visibility_schedule['elevations'][sat_idx]
max_elevation_idx = np.argmax(elevations[start_idx:end_idx+1]) + start_idx
return {
'satellite_index': sat_idx,
'start_time': time_array[start_idx],
'end_time': time_array[end_idx],
'duration': time_array[end_idx] - time_array[start_idx],
'max_elevation_time': time_array[max_elevation_idx],
'max_elevation': elevations[max_elevation_idx],
'start_elevation': elevations[start_idx],
'end_elevation': elevations[end_idx]
}
return None
9.2 链路预算分析
class LinkBudgetAnalyzer:
"""链路预算分析器"""
def __init__(self, frequency=10e9): # 默认 10GHz
self.frequency = frequency # Hz
self.c = 299792458 # 光速 (m/s)
self.k = 1.380649e-23 # 玻尔兹曼常数
def calculate_free_space_loss(self, distance):
"""计算自由空间损耗
参数:
distance: 距离 (m)
返回:
fs_loss: 自由空间损耗 (dB)
"""
wavelength = self.c / self.frequency
fs_loss = 20 * np.log10(distance) + 20 * np.log10(self.frequency) + 20 * np.log10(4 * np.pi / self.c)
return fs_loss
def calculate_link_margin(self, tx_power, tx_gain, rx_gain, distance, system_temperature=500, data_rate=100e6, additional_losses=3.0):
"""计算链路余量
参数:
tx_power: 发射功率 (dBW)
tx_gain: 发射天线增益 (dBi)
rx_gain: 接收天线增益 (dBi)
distance: 距离 (m)
system_temperature: 系统噪声温度 (K)
data_rate: 数据速率 (bps)
additional_losses: 附加损耗 (dB)
返回:
link_margin: 链路余量 (dB)
"""
# 自由空间损耗
fs_loss = self.calculate_free_space_loss(distance)
# 接收功率
rx_power = tx_power + tx_gain + rx_gain - fs_loss - additional_losses
# 噪声功率
noise_power = 10 * np.log10(self.k * system_temperature * data_rate)
# 载噪比
cnr = rx_power - noise_power
# 所需 Eb/N0 (假设 QPSK 调制，BER=1e-6)
required_ebn0 = 10.5 # dB
# 链路余量
link_margin = cnr - required_ebn0 - 10 * np.log10(data_rate)
return {
'received_power_dBw': rx_power,
'noise_power_dBw': noise_power,
'cnr_dB': cnr,
'link_margin_dB': link_margin,
'free_space_loss_dB': fs_loss
}
def analyze_constellation_link(self, constellation, ground_station, satellite_params, ground_station_params):
"""分析星座链路预算"""
results = {}
for sat_idx, sat in enumerate(constellation.satellites):
# 模拟 24 小时链路预算
time_step = 300 # 5 分钟
total_time = 24 * 3600
time_array = np.arange(0, total_time, time_step)
link_margins = []
distances = []
for t in time_array:
# 获取卫星位置
r_eci, _ = constellation.get_satellite_position_velocity(sat_idx, t)
# 计算距离
lat, lon = ground_station
lat_rad = np.radians(lat)
lon_rad = np.radians(lon)
# 考虑地球自转
omega_e = 7.2921150e-5
theta = omega_e * t
# 地面站 ECI 坐标
R_e = 6378137.0
x_ground = R_e * np.cos(lat_rad) * np.cos(lon_rad + theta)
y_ground = R_e * np.cos(lat_rad) * np.sin(lon_rad + theta)
z_ground = R_e * np.sin(lat_rad)
ground_pos = np.array([x_ground, y_ground, z_ground])
# 计算距离
distance = np.linalg.norm(r_eci - ground_pos)
distances.append(distance)
# 计算链路余量
link_result = self.calculate_link_margin(
tx_power=satellite_params['tx_power'],
tx_gain=satellite_params['tx_gain'],
rx_gain=ground_station_params['rx_gain'],
distance=distance,
system_temperature=ground_station_params['system_temperature'],
data_rate=satellite_params['data_rate']
)
link_margins.append(link_result['link_margin_dB'])
# 统计结果
results[sat_idx] = {
'min_link_margin': np.min(link_margins),
'max_link_margin': np.max(link_margins),
'avg_link_margin': np.mean(link_margins),
'min_distance': np.min(distances),
'max_distance': np.max(distances),
'avg_distance': np.mean(distances),
'link_available_percentage': np.sum(np.array(link_margins) > 0) / len(link_margins) * 100,
'time_series': {
'time': time_array,
'link_margin': link_margins,
'distance': distances
}
}
return results
10. 性能优化与扩展
10.1 并行计算优化
from concurrent.futures import ProcessPoolExecutor, ThreadPoolExecutor
import multiprocessing as mp
class ParallelOrbitSimulator:
"""并行轨道模拟器"""
def __init__(self, n_workers=None):
self.n_workers = n_workers or mp.cpu_count()
def simulate_constellation_parallel(self, constellation, duration_hours=24, time_step=60):
"""并行模拟星座"""
total_time = duration_hours * 3600
n_steps = int(total_time / time_step) + 1
n_satellites = len(constellation.satellites)
# 初始化结果数组
positions = np.zeros((n_satellites, n_steps, 3))
velocities = np.zeros((n_satellites, n_steps, 3))
# 并行计算每颗卫星的轨道
with ProcessPoolExecutor(max_workers=self.n_workers) as executor:
# 准备参数
args = [(sat_idx, constellation, total_time, time_step) for sat_idx in range(n_satellites)]
# 并行执行
results = list(executor.map(self._simulate_single_satellite, args))
# 整理结果
for sat_idx, (pos, vel) in enumerate(results):
positions[sat_idx] = pos
velocities[sat_idx] = vel
return {
'time': np.linspace(0, total_time, n_steps),
'satellite_positions': positions,
'satellite_velocities': velocities
}
@staticmethod
def _simulate_single_satellite(args):
"""模拟单颗卫星（工作进程函数）"""
sat_idx, constellation, total_time, time_step = args
n_steps = int(total_time / time_step) + 1
positions = np.zeros((n_steps, 3))
velocities = np.zeros((n_steps, 3))
for step in range(n_steps):
t = step * time_step
r_eci, v_eci = constellation.get_satellite_position_velocity(sat_idx, t)
positions[step] = r_eci
velocities[step] = v_eci
return positions, velocities
def calculate_coverage_parallel(self, constellation, ground_points, duration_hours=24, min_elevation=10):
"""并行计算覆盖"""
total_time = duration_hours * 3600
time_step = 60
n_steps = int(total_time / time_step)
n_points = len(ground_points)
# 初始化结果数组
visibility = np.zeros((n_points, n_steps))
# 并行计算每个地面点的覆盖
with ProcessPoolExecutor(max_workers=self.n_workers) as executor:
args = [(point_idx, constellation, ground_points, total_time, time_step, min_elevation) for point_idx in range(n_points)]
results = list(executor.map(self._calculate_single_point_coverage, args))
for point_idx, vis in enumerate(results):
visibility[point_idx] = vis
return {
'time': np.linspace(0, total_time, n_steps),
'visibility': visibility
}
@staticmethod
def _calculate_single_point_coverage(args):
"""计算单点覆盖（工作进程函数）"""
point_idx, constellation, ground_points, total_time, time_step, min_elevation = args
lat, lon = ground_points[point_idx]
n_steps = int(total_time / time_step)
visibility = np.zeros(n_steps)
for step in range(n_steps):
t = step * time_step
# 检查是否有卫星可见
visible = False
for sat_idx in range(len(constellation.satellites)):
r_eci, _ = constellation.get_satellite_position_velocity(sat_idx, t)
# 计算仰角
elevation = constellation.calculate_elevation(r_eci, lat, lon, t)
if elevation >= min_elevation:
visible = True
break
visibility[step] = 1 if visible else 0
return visibility
10.2 GPU 加速计算
try:
import cupy as cp
class GPUOrbitSimulator:
"""GPU 加速轨道模拟器"""
def __init__(self):
self.device = cp.cuda.Device(0)
self.mu = cp.float64(3.986004418e14)
def propagate_orbits_gpu(self, initial_states, delta_t):
"""GPU 加速轨道传播
参数:
initial_states: 初始状态数组 (n_satellites, 6)
delta_t: 时间增量 (标量或数组)
返回:
final_states: 最终状态数组
"""
# 将数据转移到 GPU
states_gpu = cp.array(initial_states, dtype=cp.float64)
if cp.isscalar(delta_t):
delta_t_gpu = cp.float64(delta_t)
else:
delta_t_gpu = cp.array(delta_t, dtype=cp.float64)
# 定义 GPU 核函数
kernel = cp.ElementwiseKernel(
in_params='T x, T y, T z, T vx, T vy, T vz, T dt',
out_params='T fx, T fy, T fz, T fvx, T fvy, T fvz',
''' // 位置和速度
T r = sqrt(x*x + y*y + z*z); T r3 = r * r * r;
// 二体问题加速度
T ax = -mu * x / r3; T ay = -mu * y / r3; T az = -mu * z / r3;
// RK4 积分
// k1
T k1_vx = ax; T k1_vy = ay; T k1_vz = az; T k1_x = vx; T k1_y = vy; T k1_z = vz;
// k2
T x2 = x + 0.5 * dt * k1_x; T y2 = y + 0.5 * dt * k1_y; T z2 = z + 0.5 * dt * k1_z;
T vx2 = vx + 0.5 * dt * k1_vx; T vy2 = vy + 0.5 * dt * k1_vy; T vz2 = vz + 0.5 * dt * k1_vz;
T r2 = sqrt(x2*x2 + y2*y2 + z2*z2); T r3_2 = r2 * r2 * r2;
T k2_vx = -mu * x2 / r3_2; T k2_vy = -mu * y2 / r3_2; T k2_vz = -mu * z2 / r3_2;
T k2_x = vx2; T k2_y = vy2; T k2_z = vz2;
// k3
T x3 = x + 0.5 * dt * k2_x; T y3 = y + 0.5 * dt * k2_y; T z3 = z + 0.5 * dt * k2_z;
T vx3 = vx + 0.5 * dt * k2_vx; T vy3 = vy + 0.5 * dt * k2_vy; T vz3 = vz + 0.5 * dt * k2_vz;
T r3 = sqrt(x3*x3 + y3*y3 + z3*z3); T r3_3 = r3 * r3 * r3;
T k3_vx = -mu * x3 / r3_3; T k3_vy = -mu * y3 / r3_3; T k3_vz = -mu * z3 / r3_3;
T k3_x = vx3; T k3_y = vy3; T k3_z = vz3;
// k4
T x4 = x + dt * k3_x; T y4 = y + dt * k3_y; T z4 = z + dt * k3_z;
T vx4 = vx + dt * k3_vx; T vy4 = vy + dt * k3_vy; T vz4 = vz + dt * k3_vz;
T r4 = sqrt(x4*x4 + y4*y4 + z4*z4); T r3_4 = r4 * r4 * r4;
T k4_vx = -mu * x4 / r3_4; T k4_vy = -mu * y4 / r3_4; T k4_vz = -mu * z4 / r3_4;
T k4_x = vx4; T k4_y = vy4; T k4_z = vz4;
// 最终状态
fx = x + dt/6.0 * (k1_x + 2.0*k2_x + 2.0*k3_x + k4_x);
fy = y + dt/6.0 * (k1_y + 2.0*k2_y + 2.0*k3_y + k4_y);
fz = z + dt/6.0 * (k1_z + 2.0*k2_z + 2.0*k3_z + k4_z);
fvx = vx + dt/6.0 * (k1_vx + 2.0*k2_vx + 2.0*k3_vx + k4_vx);
fvy = vy + dt/6.0 * (k1_vy + 2.0*k2_vy + 2.0*k3_vy + k4_vy);
fvz = vz + dt/6.0 * (k1_vz + 2.0*k2_vz + 2.0*k3_vz + k4_vz);
''',
name='propagate_orbit',
preamble='T mu = 3.986004418e14;'
)
# 执行核函数
x = states_gpu[:, 0]
y = states_gpu[:, 1]
z = states_gpu[:, 2]
vx = states_gpu[:, 3]
vy = states_gpu[:, 4]
vz = states_gpu[:, 5]
fx, fy, fz, fvx, fvy, fvz = kernel(x, y, z, vx, vy, vz, delta_t_gpu)
# 组合结果并传回 CPU
final_states = cp.stack([fx, fy, fz, fvx, fvy, fvz], axis=1)
return cp.asnumpy(final_states)
except ImportError:
print("CuPy 未安装，GPU 加速不可用")
11. 实际应用案例
11.1 Starlink 星座模拟案例
def simulate_starlink_constellation():
"""模拟 Starlink 星座"""
print("开始模拟 Starlink 星座...")
# 创建模拟器
simulator = SatelliteOrbitSimulator()
# Starlink Phase 1 参数（简化）
# 轨道壳 1: 550km, 53 度倾角，72 个轨道面，每面 22 颗卫星
constellation_1 = simulator.create_walker_constellation(
T=1584, # 72 * 22
P=72, F=1, altitude=550, inclination=53.0
)
# 轨道壳 2: 540km, 53.2 度倾角，72 个轨道面，每面 22 颗卫星
constellation_2 = simulator.create_walker_constellation(
T=1584, P=72, F=1, altitude=540, inclination=53.2
)
# 模拟运行（简化版，使用较小时长和步长）
print("模拟星座运行...")
simulation_data = simulator.simulate_constellation(
duration_hours=2, # 2 小时模拟
time_step=300 # 5 分钟步长
)
# 定义地面站
ground_stations = [
(40.7128, -74.0060, 'New York'),
(51.5074, -0.1278, 'London'),
(35.6762, 139.6503, 'Tokyo'),
(-33.8688, 151.2093, 'Sydney'),
(-23.5505, -46.6333, 'Sao Paulo'),
(28.6139, 77.2090, 'Delhi'),
(39.9042, 116.4074, 'Beijing')
]
# 计算覆盖
print("计算地面站覆盖...")
coverage_data = simulator.calculate_coverage(
ground_stations=ground_stations,
duration_hours=2,
min_elevation=25 # Starlink 最小仰角要求
)
# 分析结果
analyzer = OrbitAnalyzer()
analysis_results = analyzer.analyze_constellation_coverage(coverage_data)
# 可视化
visualizer = OrbitVisualizer()
# 绘制 3D 轨道
print("生成可视化...")
fig_3d, ax_3d = visualizer.plot_3d_orbit(
simulation_data, title="Starlink 星座轨道模拟"
)
plt.savefig('starlink_3d_orbit.png', dpi=300, bbox_inches='tight')
# 绘制地面站覆盖分析
fig_cov, axes_cov = visualizer.plot_coverage_analysis(coverage_data, analysis_results)
plt.savefig('starlink_coverage_analysis.png', dpi=300, bbox_inches='tight')
# 打印关键结果
print("\n=== Starlink 星座模拟结果 ===")
print(f"总卫星数：{constellation_1.T + constellation_2.T}")
print(f"轨道壳 1: {constellation_1.T}颗卫星，{constellation_1.P}个轨道面")
print(f"轨道壳 2: {constellation_2.T}颗卫星，{constellation_2.P}个轨道面")
print("\n地面站覆盖统计:")
for name, data in analysis_results['station_coverage'].items():
print(f" {name}: {data['coverage_percentage']:.1f}% 覆盖")
print("\n全局统计:")
global_metrics = analysis_results['global_metrics']
print(f" 全球平均覆盖：{global_metrics['global_coverage_percentage']:.1f}%")
print(f" 平均同时覆盖地面站数：{global_metrics['avg_simultaneous_coverage']:.1f}")
plt.show()
return {
'constellation_1': constellation_1,
'constellation_2': constellation_2,
'simulation_data': simulation_data,
'coverage_data': coverage_data,
'analysis_results': analysis_results
}
def analyze_link_budget_example():
"""链路预算分析示例"""
print("\n=== 链路预算分析示例 ===")
# 创建示例星座
simulator = SatelliteOrbitSimulator()
constellation = simulator.create_walker_constellation(
T=66, P=6, F=1, altitude=550, inclination=53.0
)
# 定义地面站
ground_station = (40.7128, -74.0060) # 纽约
# 卫星参数
satellite_params = {
'tx_power': 10.0, # dBW (10W)
'tx_gain': 30.0, # dBi
'data_rate': 100e6 # 100 Mbps
}
# 地面站参数
ground_station_params = {
'rx_gain': 45.0, # dBi
'system_temperature': 150.0 # K
}
# 创建链路预算分析器
link_analyzer = LinkBudgetAnalyzer(frequency=12e9) # 12 GHz
# 分析链路
link_results = link_analyzer.analyze_constellation_link(
constellation, ground_station, satellite_params, ground_station_params
)
# 打印结果
print("\n链路预算分析结果:")
for sat_idx, result in link_results.items():
print(f"\n卫星 {sat_idx}:")
print(f" 最小链路余量：{result['min_link_margin']:.2f} dB")
print(f" 最大链路余量：{result['max_link_margin']:.2f} dB")
print(f" 平均链路余量：{result['avg_link_margin']:.2f} dB")
print(f" 链路可用时间：{result['link_available_percentage']:.1f}%")
print(f" 最小距离：{result['min_distance']/1000:.1f} km")
print(f" 最大距离：{result['max_distance']/1000:.1f} km")
return link_results
def performance_benchmark():
"""性能基准测试"""
print("\n=== 性能基准测试 ===")
import time
# 测试不同规模的星座
constellation_sizes = [10, 50, 100, 500]
simulation_times = []
for size in constellation_sizes:
print(f"\n测试 {size} 颗卫星的星座...")
# 创建星座
simulator = SatelliteOrbitSimulator()
constellation = simulator.create_walker_constellation(
T=size, P=int(np.sqrt(size)), F=1, altitude=550, inclination=53.0
)
# 计时
start_time = time.time()
# 模拟 2 小时
simulation_data = simulator.simulate_constellation(
duration_hours=2, time_step=300
)
elapsed_time = time.time() - start_time
simulation_times.append(elapsed_time)
print(f" 模拟时间：{elapsed_time:.2f} 秒")
print(f" 每颗卫星平均时间：{elapsed_time/size:.4f} 秒")
# 绘制性能图表
plt.figure(figsize=(10, 6))
plt.plot(constellation_sizes, simulation_times, 'bo-', linewidth=2, markersize=8)
plt.xlabel('卫星数量')
plt.ylabel('模拟时间 (秒)')
plt.title('星座模拟性能基准测试')
plt.grid(True, alpha=0.3)
plt.xscale('log')
plt.yscale('log')
# 添加趋势线
x_log = np.log10(constellation_sizes)
y_log = np.log10(simulation_times)
coeffs = np.polyfit(x_log, y_log, 1)
trend_line = 10**coeffs[1] * np.array(constellation_sizes)**coeffs[0]
plt.plot(constellation_sizes, trend_line, 'r--', alpha=0.7, label=f'趋势：O(n^{coeffs[0]:.2f})')
plt.legend()
plt.savefig('performance_benchmark.png', dpi=300, bbox_inches='tight')
plt.show()
return constellation_sizes, simulation_times
11.2 完整示例：运行星座模拟
def main():
"""主函数：运行完整的星座模拟"""
print("=" * 60)
print("Python 卫星通信模拟：低轨星座的轨道力学计算")
print("=" * 60)
# 1. 模拟 Starlink 星座
starlink_results = simulate_starlink_constellation()
# 2. 链路预算分析
link_results = analyze_link_budget_example()
# 3. 性能基准测试
perf_results = performance_benchmark()
# 4. 创建动画（可选）
print("\n=== 生成轨道动画 ===")
animator = OrbitAnimator(starlink_results['simulation_data'])
try:
animator.create_3d_animation(
output_file='constellation_animation.mp4',
fps=30, dpi=100
)
print("动画生成完成!")
except Exception as e:
print(f"动画生成失败：{e}")
print("请确保已安装 ffmpeg")
print("\n" + "=" * 60)
print("模拟完成!")
print("=" * 60)
return {
'starlink': starlink_results,
'link_budget': link_results,
'performance': perf_results
}
if __name__ == "__main__":
# 运行主程序
results = main()
12. 结论与展望
12.1 主要成果
本文构建了一个完整的 Python 低轨卫星星座轨道力学计算框架，具有以下特点：
- 完整的轨道力学模型：从基本的二体问题到包含 J2 摄动、大气阻力、太阳辐射压力的精确模型
- 灵活的星座设计：支持 Walker 星座等多种构型
- 全面的分析工具：覆盖计算、可见性分析、链路预算等关键功能
- 高效的计算方法：支持并行计算和 GPU 加速
- 丰富的可视化：3D 轨道显示、星下点轨迹、覆盖分析图表
- 实用的应用案例：Starlink 星座模拟、性能基准测试
12.2 技术挑战与解决方案
在开发过程中遇到的主要挑战及解决方案：
- 计算效率问题：通过并行计算和 GPU 加速优化大规模星座模拟
- 数值稳定性：使用高阶数值积分方法和自适应步长控制
- 坐标系统转换：实现完整的 ECI、ECEF、轨道平面坐标转换
- 可视化复杂性：开发专门的 3D 可视化工具，处理地球曲面显示
12.3 未来发展方向
- 更精确的摄动模型：加入月球和太阳引力摄动、地球高阶重力场模型
- 星座优化算法：使用机器学习优化星座构型参数
- 实时轨道确定：集成 GPS 测量数据，实现实时轨道确定
- 碰撞规避：加入卫星碰撞风险评估和规避机动规划
- 网络模拟：扩展为完整的卫星通信网络模拟器
- 云平台部署：将系统部署为 Web 服务，提供在线模拟能力
12.4 实际应用价值
本框架可应用于：
- 星座设计与评估：帮助设计新的卫星星座系统
- 任务规划：支持卫星任务规划和地面站调度
- 教育培训：用于航空航天专业的教学和培训
- 科研分析：支持卫星通信相关科学研究
- 系统验证：验证商业星座系统的性能指标
12.5 代码获取与使用
完整代码可通过开源社区获取，安装和使用方法：
依赖库包括：
- NumPy, SciPy: 科学计算
- Matplotlib: 数据可视化
- Pandas: 数据处理
- Astropy: 天文计算（可选）
- CuPy: GPU 加速（可选）
附录
A. 轨道要素转换公式
详细的开普勒轨道要素与笛卡尔坐标转换公式...
B. 常用轨道参数表
典型低轨卫星星座参数参考...
C. 参考文献
- Vallado, D. A. (2013). Fundamentals of Astrodynamics and Applications.
- Wertz, J. R. (2011). Space Mission Engineering: The New SMAD.
- Lang, T. J. (2011). Satellite Communications.
D. 术语表
- LEO: 低地球轨道
- ECI: 地心惯性坐标系
- ECEF: 地固坐标系
- RAAN: 升交点赤经
- J2 摄动: 地球扁率引起的轨道摄动
通过本文介绍的 Python 卫星轨道模拟框架，研究人员和工程师可以快速构建、分析和优化低轨卫星星座系统。该框架结合了理论深度和工程实用性，为卫星通信系统的设计和分析提供了强大的工具支持。随着商业航天和卫星互联网的快速发展，这样的模拟工具将变得越来越重要。
