国際宇宙ステーション（ISS）は高度約400 kmの円軌道を周回していますが、何もしなければ数年のうちに大気圏に再突入して燃え尽きてしまいます。一方で、GPS衛星は高度約20,200 kmで数十年にわたって安定に運用されています。この差は何が原因でしょうか？ 「今から5年後、衛星はどこにいるか？」という問いに答えるためには、軌道の長期伝播（long-term propagation）が必要です。
短期予測（数日程度）であれば、衛星の運動は二体問題に近いケプラー軌道で十分に近似できます。しかし、数週間〜数年という時間スケールになると、大気抵抗・地球の非球対称性（J2）・月太陽の重力・太陽輻射圧といった摂動が蓄積し、軌道は大きく変化します。長期伝播は、これらの摂動をすべて考慮して軌道の進化を追跡する技術です。
軌道の長期伝播を理解すると、以下の応用が可能になります。
- 軌道寿命推定: 大気抵抗による軌道減衰を予測し、衛星がいつ再突入するかを見積もる
- コンステレーション維持: 数百機の衛星群の軌道進化を予測し、ステーションキーピング計画を立てる
- スペースデブリの長期追跡: 数十年後のデブリの分布を予測し、衝突リスクを評価する
本記事の内容
- 短期予測と長期予測の違い
- ガウスの摂動方程式（軌道要素の時間変化率）
- 主要な摂動源（J2、大気抵抗、第三体、太陽輻射圧）
- 数値積分法（RK4, RK78, DOP853）の選び方
- Pythonで低軌道衛星の数年間の軌道進化をシミュレーション
- 軌道寿命の推定
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
短期予測と長期予測の違い
短期予測（数時間〜数日）
短期予測では、衛星の位置を高い精度（数百m〜数km）で求めることが目的です。このためには、すべての摂動力を正確にモデル化し、直交座標系の運動方程式を小さな時間刻みで数値積分します。TLE/SGP4モデルは数日間の予測に適したこの種の手法です。
短期予測では、衛星の現在位置を正確に知ることが重要ですが、軌道の形状そのものの変化は小さいため、ケプラー要素はほぼ一定とみなせます。
長期予測（数週間〜数年）
長期予測では、軌道の「形状と向き」がどう変化するかに関心があります。半長軸がどう縮小するか、傾斜角がどう変動するか、離心率がどう振動するか — これらを何年にもわたって追跡します。
位置の絶対精度は二の次で、軌道要素のトレンド（傾向）が正しいことが重要です。そのため、短周期変動を平均化して計算効率を上げる手法（解析的半解析的手法）や、軌道要素に直接作用する摂動方程式を使う方法が好まれます。
二つのアプローチ
長期伝播には大きく分けて2つの方法があります。
直交座標法（Cowell法）: 運動方程式 $\ddot{\bm{r}} = -\mu\bm{r}/r^3 + \bm{a}_{\text{pert}}$ を直接数値積分します。実装が単純で、どんな摂動でも扱えますが、長期積分では計算コストが高く、数値誤差が蓄積しやすいです。
変分法（VOP法）: ラグランジュまたはガウスの摂動方程式を使い、軌道要素の時間変化率を直接積分します。摂動がないときは軌道要素が一定（右辺がゼロ）なので、長期積分でも数値的に安定です。
長期伝播では変分法が有利ですが、本記事ではまず両方の考え方を理解した上で、実装にはCowell法をベースとしたアプローチを使います。ガウスの摂動方程式は、結果を解釈するための理論的な枠組みとして活用します。
ガウスの摂動方程式
摂動加速度の分解
ガウスの摂動方程式では、摂動加速度 $\bm{a}_{\text{pert}}$ を衛星の瞬時軌道面内の3方向に分解します。
- $a_R$: 動径方向（地球から衛星に向かう方向、正が外向き）
- $a_S$: 接線方向（軌道面内で速度に垂直、正が軌道運動の向き）
- $a_W$: 法線方向（軌道面に垂直、正が角運動量ベクトルの向き）
この座標系は衛星の位置に応じて回転するRSW座標系と呼ばれます。イメージとしては、車のダッシュボード（R: 前方向/S: 横方向/W: 上方向）のように衛星と一緒に動く座標系です。
方程式の全体像
ガウスの摂動方程式の完全な形は次の通りです。$p = a(1-e^2)$ は半通径です。
半長軸の変化率:
$$ \frac{da}{dt} = \frac{2}{n\sqrt{1-e^2}}\left[e\sin\nu \cdot a_R + \frac{p}{r} \cdot a_S\right] $$
離心率の変化率:
$$ \frac{de}{dt} = \frac{\sqrt{1-e^2}}{na}\left[\sin\nu \cdot a_R + \left(\cos\nu + \frac{e + \cos\nu}{1 + e\cos\nu}\right) a_S\right] $$
傾斜角の変化率:
$$ \frac{di}{dt} = \frac{r\cos(\omega + \nu)}{na^2\sqrt{1-e^2}} a_W $$
昇交点赤経の変化率:
$$ \frac{d\Omega}{dt} = \frac{r\sin(\omega + \nu)}{na^2\sqrt{1-e^2}\sin i} a_W $$
近地点引数の変化率:
$$ \frac{d\omega}{dt} = \frac{\sqrt{1-e^2}}{nae}\left[-\cos\nu \cdot a_R + \sin\nu\left(1 + \frac{r}{p}\right) a_S\right] – \frac{r\sin(\omega+\nu)\cos i}{na^2\sqrt{1-e^2}\sin i} a_W $$
物理的解釈
これらの方程式から、各方向の摂動がどの軌道要素に影響するかがわかります。
動径方向 $a_R$ と接線方向 $a_S$ の摂動: 半長軸 $a$、離心率 $e$、近地点引数 $\omega$ を変化させます。特に $a_S$（接線方向の力）は半長軸の変化に直接寄与し、大気抵抗のような進行方向と逆向きの力は $a$ を減少させます。
法線方向 $a_W$ の摂動: 傾斜角 $i$ と昇交点赤経 $\Omega$ を変化させます。軌道面の向きを変える力です。第三体摂動の $z$ 成分がこの役割を果たします。
重要な洞察として、接線方向の力が最もエネルギー（したがって半長軸）を効率的に変化させることがわかります。ロケットのマヌーバが接線方向に行われるのは、この理由からです。
ガウスの摂動方程式は「どの摂動がどの軌道要素を変えるか」を理解するための道具です。では、実際にどのような摂動源が衛星に作用するのかを、一つずつ見ていきましょう。
主要な摂動源
J2摂動（地球の扁平率）
地球は完全な球ではなく、赤道方向にわずかに膨らんだ扁球体です。この非対称性による追加の重力場は、球面調和関数で展開されます。最大の項が $J_2$ 項です。
$$ a_{J2,R} = -\frac{3}{2}\frac{J_2 \mu R_E^2}{r^4}\left(1 – 3\sin^2 i \sin^2(\omega + \nu)\right) $$
$$ a_{J2,S} = -\frac{3}{2}\frac{J_2 \mu R_E^2}{r^4}\sin^2 i \sin 2(\omega + \nu) $$
$$ a_{J2,W} = -\frac{3}{2}\frac{J_2 \mu R_E^2}{r^4}\sin 2i \sin(\omega + \nu) $$
ここで $J_2 = 1.08263 \times 10^{-3}$、$R_E = 6378.137$ km です。J2摂動の永年的な効果は
$$ \dot{\Omega}_{J2} = -\frac{3}{2}nJ_2\left(\frac{R_E}{p}\right)^2 \cos i $$
$$ \dot{\omega}_{J2} = \frac{3}{4}nJ_2\left(\frac{R_E}{p}\right)^2 (5\cos^2 i – 1) $$
これらはよく知られたJ2による昇交点の回帰と近地点の移動です。太陽同期軌道は $\dot{\Omega}_{J2} = 360°/\text{year}$ を達成する傾斜角を選ぶことで実現されます。
大気抵抗
高度約1000 km以下の低軌道衛星では、希薄ながらも存在する大気が衛星に抗力を及ぼします。大気抵抗は常に衛星の速度と逆向きに作用するため、ガウスの摂動方程式では主に $a_S < 0$（逆接線方向）の成分を持ちます。
大気抵抗による加速度は
$$ \bm{a}_{\text{drag}} = -\frac{1}{2}\rho v_{\text{rel}}^2 \frac{C_D A}{m}\hat{\bm{v}}_{\text{rel}} $$
ここで $\rho$ は大気密度、$v_{\text{rel}}$ は大気に対する衛星の相対速度、$C_D$ は抗力係数（通常2.0〜2.5）、$A$ は衛星の断面積、$m$ は衛星の質量です。$C_D A / m$ の逆数を弾道係数 $B = m/(C_D A)$ と呼びます。
大気密度 $\rho$ は高度に対して指数関数的に減少します。
$$ \rho(h) \approx \rho_0 \exp\left(-\frac{h – h_0}{H}\right) $$
ここで $H$ はスケールハイト（高度約400 kmでは約60 km）、$\rho_0$ は参照高度 $h_0$ での密度です。太陽活動が活発な時期には大気が膨張し、同じ高度でも密度が数倍〜10倍に増加します。
大気抵抗は軌道からエネルギーを奪い続けるため、半長軸 $a$ が永年的に減少します。楕円軌道の場合、大気は近地点付近でのみ有意な密度を持つため、近地点でブレーキがかかります。これにより遠地点高度が下がり、軌道は徐々に円軌道に近づいた後、最終的にスパイラル状に降下します。
第三体摂動（月・太陽）
月と太陽による第三体摂動については前の記事で詳しく解説しました。ここでは長期伝播の文脈で重要なポイントだけをまとめます。
第三体摂動の加速度は
$$ \bm{a}_3 = GM_3\left(\frac{\bm{r}_3 – \bm{r}}{|\bm{r}_3 – \bm{r}|^3} – \frac{\bm{r}_3}{r_3^3}\right) $$
高軌道（GEO、MEO）では第三体摂動が最も重要な摂動源です。低軌道（LEO）では大気抵抗やJ2に比べて小さいですが、長期伝播では無視できません。
太陽輻射圧（SRP）
太陽から放射される光子が衛星に当たると、その運動量が衛星に伝わって力が作用します。これが太陽輻射圧（Solar Radiation Pressure, SRP）です。
$$ \bm{a}_{\text{SRP}} = -P_{\odot}\frac{C_R A_{\odot}}{m}\frac{AU^2}{r_{\odot}^2}\hat{\bm{r}}_{\odot} $$
ここで $P_{\odot} \approx 4.56 \times 10^{-6}$ N/m² は1 AUでの太陽輻射圧、$C_R$ は反射係数（1〜2）、$A_{\odot}$ は太陽に対する投影面積、$r_{\odot}$ は太陽までの距離です。
SRPは大きな太陽電池パドルを持つ衛星や、面積質量比が大きい軽量デブリに対して特に重要です。離心率に周期的な変動を引き起こし、GPS衛星の軌道決定では主要な誤差源の一つです。
各摂動源の特徴を理解したところで、次にこれらの摂動を含む運動方程式を数値的に解く方法を見ていきます。
数値積分法の選択
軌道力学で使われる数値積分法
軌道の長期伝播では、微分方程式を数値的に解く必要があります。使われる主な手法を比較します。
RK4（古典的ルンゲクッタ4次法）
$$ \bm{y}_{n+1} = \bm{y}_n + \frac{h}{6}(\bm{k}_1 + 2\bm{k}_2 + 2\bm{k}_3 + \bm{k}_4) $$
固定刻みの4次精度法です。実装が簡単で理解しやすいですが、刻み幅の自動制御ができないため、長期積分では効率が悪くなります。短期のテストやプロトタイピングには適しています。
RK7(8)（Dormand-Prince 8次法, DOP853）
高次の埋め込み型ルンゲクッタ法です。7次の解と8次の解を同時に計算し、その差から誤差を推定して刻み幅を自動的に調整します。SciPyの solve_ivp で method='DOP853' として利用でき、軌道力学の長期積分で広く使われます。
その他の手法
Adams-Bashforth-Moulton法（多段法）: 過去の複数のステップの情報を使うため、関数評価の回数が少なくて済みます。NASAのJPLなどで惑星探査機の軌道計算に使われてきました。
Gauss-Legendre法（陰的ルンゲクッタ法）: エネルギー保存性に優れたシンプレクティック積分法の一種です。保存系（大気抵抗がない場合）の長期積分に適しています。
長期積分での注意点
軌道の長期伝播で特に注意すべきは以下の点です。
数値誤差の蓄積: 何百万ステップもの積分を行うと、1ステップあたりの微小な誤差が蓄積します。相対許容誤差 rtol と絶対許容誤差 atol を十分小さく設定する必要があります（典型的には rtol=1e-12, atol=1e-14）。
エネルギー保存の監視: 保存系（大気抵抗なし）では、軌道エネルギー $\mathcal{E} = v^2/2 – \mu/r$ が保存されます。この量の変動をモニターすることで、数値積分の品質を確認できます。
刻み幅の適応: 近地点付近（速度が大きく、曲率が大きい）では細かい刻みが必要で、遠地点付近では粗い刻みで十分です。適応刻み幅制御を持つ積分法（DOP853など）が効率的です。
理論的な準備が整ったので、いよいよPythonで長期伝播のシミュレーションを実装しましょう。
Pythonで低軌道衛星の軌道進化をシミュレーション
モデルの設計
ここでは、高度600 kmの円軌道にある衛星の3年間の軌道進化をシミュレーションします。以下の摂動を考慮します。
- J2摂動（地球の扁平率）
- 大気抵抗（指数大気モデル）
- 月の第三体摂動（簡略化）
まず、定数と大気モデルを定義します。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# --- 定数 ---
mu = 3.986004418e14 # 地球の重力パラメータ [m^3/s^2]
R_E = 6378.137e3 # 地球の赤道半径 [m]
J2 = 1.08263e-3 # J2係数
omega_earth = 7.2921159e-5 # 地球の自転角速度 [rad/s]
# 大気モデル（指数大気の参照値テーブル）
# (基準高度 [km], 基準密度 [kg/m^3], スケールハイト [km])
atm_table = [
(200, 2.789e-10, 37.105),
(250, 7.248e-11, 45.546),
(300, 2.418e-11, 53.628),
(350, 9.518e-12, 53.298),
(400, 3.725e-12, 58.515),
(450, 1.585e-12, 60.828),
(500, 6.967e-13, 63.822),
(550, 1.454e-13, 60.828),
(600, 3.614e-14, 71.835),
(700, 5.245e-15, 88.667),
(800, 1.170e-15, 124.64),
(900, 3.614e-16, 181.05),
(1000, 1.454e-16, 268.00),
]
def atmospheric_density(h_km):
"""指数大気モデルによる密度推定 [kg/m^3]"""
if h_km > 1000:
return 0.0
if h_km < 200:
h_km = 200
# テーブルから適切な区間を選択
for k in range(len(atm_table) - 1):
if atm_table[k][0] <= h_km < atm_table[k + 1][0]:
h0, rho0, H = atm_table[k]
return rho0 * np.exp(-(h_km - h0) / H)
h0, rho0, H = atm_table[-1]
return rho0 * np.exp(-(h_km - h0) / H)
次に、各摂動の加速度を計算する関数を実装します。
# --- 衛星パラメータ ---
C_D = 2.2 # 抗力係数
A_sat = 10.0 # 断面積 [m^2]（大型衛星）
m_sat = 1000.0 # 質量 [kg]
B_coeff = m_sat / (C_D * A_sat) # 弾道係数 [kg/m^2]
# --- 月のパラメータ ---
mu_moon = 4.9028e12
r_moon_orbit = 384400e3
n_moon = 2 * np.pi / (27.3217 * 86400)
i_moon = np.radians(23.44)
def accel_j2(r_vec):
"""J2摂動加速度"""
r = np.linalg.norm(r_vec)
x, y, z = r_vec
factor = -1.5 * J2 * mu * R_E**2 / r**5
z2_r2 = (z / r)**2
ax = factor * x * (1 - 5 * z2_r2)
ay = factor * y * (1 - 5 * z2_r2)
az = factor * z * (3 - 5 * z2_r2)
return np.array([ax, ay, az])
def accel_drag(r_vec, v_vec):
"""大気抵抗加速度"""
r = np.linalg.norm(r_vec)
h_km = (r - R_E) / 1e3
if h_km > 1000 or h_km < 0:
return np.zeros(3)
rho = atmospheric_density(h_km)
if rho < 1e-20:
return np.zeros(3)
# 大気の共回転速度
v_atm = np.cross([0, 0, omega_earth], r_vec)
v_rel = v_vec - v_atm
v_rel_mag = np.linalg.norm(v_rel)
a_drag = -0.5 * rho * v_rel_mag / B_coeff * v_rel
return a_drag
def accel_moon(r_vec, t):
"""月の第三体摂動加速度"""
theta = n_moon * t
r_moon_vec = r_moon_orbit * np.array([
np.cos(theta),
np.sin(theta) * np.cos(i_moon),
np.sin(theta) * np.sin(i_moon)
])
r_sat_moon = r_moon_vec - r_vec
d = np.linalg.norm(r_sat_moon)
r3 = np.linalg.norm(r_moon_vec)
return mu_moon * (r_sat_moon / d**3 - r_moon_vec / r3**3)
運動方程式を統合して数値積分を実行します。
def eom_full(t, state):
"""全摂動を含む運動方程式"""
r_vec = state[0:3]
v_vec = state[3:6]
r = np.linalg.norm(r_vec)
# 地球の点質量重力
a_gravity = -mu / r**3 * r_vec
# 摂動
a_j2 = accel_j2(r_vec)
a_drag = accel_drag(r_vec, v_vec)
a_moon = accel_moon(r_vec, t)
a_total = a_gravity + a_j2 + a_drag + a_moon
return np.concatenate([v_vec, a_total])
# --- ケプラー要素 → 直交座標 ---
def kep_to_cart(a, e, i, RAAN, omega, nu, mu_val):
"""ケプラー要素から位置・速度ベクトルを計算"""
p = a * (1 - e**2)
r_mag = p / (1 + e * np.cos(nu))
r_pqw = r_mag * np.array([np.cos(nu), np.sin(nu), 0.0])
v_pqw = np.sqrt(mu_val / p) * np.array([-np.sin(nu), e + np.cos(nu), 0.0])
cO, sO = np.cos(RAAN), np.sin(RAAN)
ci, si = np.cos(i), np.sin(i)
cw, sw = np.cos(omega), np.sin(omega)
R = np.array([
[cO*cw - sO*sw*ci, -cO*sw - sO*cw*ci, sO*si],
[sO*cw + cO*sw*ci, -sO*sw + cO*cw*ci, -cO*si],
[sw*si, cw*si, ci]
])
return R @ r_pqw, R @ v_pqw
def cart_to_kep(r_vec, v_vec, mu_val):
"""位置・速度から軌道要素を計算"""
r = np.linalg.norm(r_vec)
v = np.linalg.norm(v_vec)
h_vec = np.cross(r_vec, v_vec)
h = np.linalg.norm(h_vec)
e_vec = np.cross(v_vec, h_vec) / mu_val - r_vec / r
e = np.linalg.norm(e_vec)
energy = v**2 / 2 - mu_val / r
if abs(energy) < 1e-20:
a = np.inf
else:
a = -mu_val / (2 * energy)
i = np.arccos(np.clip(h_vec[2] / h, -1, 1))
n_vec = np.cross([0, 0, 1], h_vec)
n = np.linalg.norm(n_vec)
if n > 1e-10:
RAAN = np.arccos(np.clip(n_vec[0] / n, -1, 1))
if n_vec[1] < 0:
RAAN = 2 * np.pi - RAAN
else:
RAAN = 0.0
if n > 1e-10 and e > 1e-10:
omega = np.arccos(np.clip(np.dot(n_vec, e_vec) / (n * e), -1, 1))
if e_vec[2] < 0:
omega = 2 * np.pi - omega
else:
omega = 0.0
if e > 1e-10:
nu = np.arccos(np.clip(np.dot(e_vec, r_vec) / (e * r), -1, 1))
if np.dot(r_vec, v_vec) < 0:
nu = 2 * np.pi - nu
else:
nu = 0.0
return a, e, i, RAAN, omega, nu
# --- 初期条件（高度600kmの円軌道, 傾斜角51.6°） ---
a0 = R_E + 600e3 # 半長軸
e0 = 0.001 # ほぼ円軌道
i0 = np.radians(51.6) # ISS軌道に近い傾斜角
RAAN0 = 0.0
omega0 = 0.0
nu0 = 0.0
r0, v0 = kep_to_cart(a0, e0, i0, RAAN0, omega0, nu0, mu)
state0 = np.concatenate([r0, v0])
# --- 3年間の積分 ---
T_years = 3.0
t_span = (0, T_years * 365.25 * 86400)
t_eval = np.linspace(t_span[0], t_span[1], 8000)
print("長期伝播シミュレーション開始...")
sol = solve_ivp(eom_full, t_span, state0,
method='DOP853', t_eval=t_eval,
rtol=1e-12, atol=1e-14,
max_step=60.0)
print(f"完了: {sol.success}, ステップ数: {sol.t.shape[0]}")
# --- 軌道要素の抽出 ---
N = len(sol.t)
a_arr = np.zeros(N)
e_arr = np.zeros(N)
i_arr = np.zeros(N)
RAAN_arr = np.zeros(N)
h_arr = np.zeros(N) # 高度
for k in range(N):
r_k = sol.y[0:3, k]
v_k = sol.y[3:6, k]
a_arr[k], e_arr[k], i_arr[k], RAAN_arr[k], _, _ = \
cart_to_kep(r_k, v_k, mu)
h_arr[k] = (np.linalg.norm(r_k) - R_E) / 1e3 # [km]
# 近地点・遠地点高度
hp_arr = a_arr * (1 - e_arr) - R_E # [m]
ha_arr = a_arr * (1 + e_arr) - R_E # [m]
t_years_arr = sol.t / (365.25 * 86400)
結果を可視化します。
fig, axes = plt.subplots(4, 1, figsize=(10, 14), sharex=True)
# 高度（近地点・遠地点）
axes[0].plot(t_years_arr, hp_arr / 1e3, color='steelblue',
linewidth=0.5, label='Perigee')
axes[0].plot(t_years_arr, ha_arr / 1e3, color='darkorange',
linewidth=0.5, label='Apogee')
axes[0].set_ylabel('Altitude [km]')
axes[0].set_title('低軌道衛星の軌道進化（高度600km, i=51.6°, 3年間）')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 離心率
axes[1].plot(t_years_arr, e_arr, color='crimson', linewidth=0.5)
axes[1].set_ylabel('Eccentricity e')
axes[1].grid(True, alpha=0.3)
# 傾斜角
axes[2].plot(t_years_arr, np.degrees(i_arr), color='green', linewidth=0.5)
axes[2].set_ylabel('Inclination [deg]')
axes[2].grid(True, alpha=0.3)
# 昇交点赤経
axes[3].plot(t_years_arr, np.degrees(RAAN_arr) % 360, color='purple',
linewidth=0.5)
axes[3].set_ylabel('RAAN [deg]')
axes[3].set_xlabel('Time [years]')
axes[3].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('long_term_propagation_leo.png', dpi=150, bbox_inches='tight')
plt.show()
このシミュレーション結果から、低軌道衛星の軌道進化の重要な特徴が読み取れます。
-
近地点高度・遠地点高度は大気抵抗により徐々に低下する — 大気抵抗は主に近地点でブレーキをかけるため、遠地点高度が先に低下し、軌道は次第に円形化します。高度が下がるほど大気密度が増すため、降下速度は加速的に増加します。
-
離心率は複雑な振動をしながら変化する — J2摂動が離心率に短周期の振動を引き起こし、大気抵抗が楕円→円の方向への長期トレンドを生み出しています。
-
傾斜角はJ2の効果でわずかに振動する — 低軌道では第三体摂動よりもJ2の効果が支配的で、傾斜角の変化は小さいです。長い目で見ると月・太陽摂動による微小な長期変動もあります。
-
昇交点赤経はJ2により永年的に回帰する — 傾斜角51.6°の場合、$\dot{\Omega}_{J2}$ は負の値を取り、昇交点は年間数十度の割合で西方に移動します。これがJ2摂動の最も顕著な長期効果です。
これらの結果を使って、次に衛星の軌道寿命を推定する方法を考えてみましょう。
軌道寿命の推定
数値的な軌道寿命
シミュレーション結果から、衛星の高度がある閾値（一般に120 km）を下回ったとき、衛星は再突入すると判定できます。上のシミュレーションでは3年では再突入に至りませんでしたが、高度降下のレートから外挿することで寿命を概算できます。
# --- 軌道寿命の概算 ---
# 半長軸の時間変化率を計算
a_km = a_arr / 1e3
da_dt = np.gradient(a_km, sol.t / (365.25 * 86400)) # [km/year]
# 最後の1年のデータから平均降下率を推定
mask_last = t_years_arr > (T_years - 1.0)
avg_descent_rate = np.mean(da_dt[mask_last]) # [km/year]
# 残り高度から寿命を推定
h_final = (a_arr[-1] - R_E) / 1e3 # 現在の平均高度 [km]
h_reentry = 120.0 # 再突入高度 [km]
if avg_descent_rate < 0:
remaining_years = (h_final - h_reentry) / abs(avg_descent_rate)
total_life = T_years + remaining_years
print(f"シミュレーション終了時の平均高度: {h_final:.1f} km")
print(f"平均降下率（最後の1年）: {avg_descent_rate:.2f} km/year")
print(f"推定残存寿命: {remaining_years:.1f} 年")
print(f"推定総軌道寿命: {total_life:.1f} 年")
else:
print(f"高度が上昇しています（異常）: rate = {avg_descent_rate:.4f} km/year")
この結果は、衛星の質量・面積比（弾道係数）と太陽活動の状態に大きく依存します。太陽活動極大期には大気が膨張するため、軌道寿命は極小期に比べて数分の1に短くなることもあります。
解析的な軌道寿命推定
大気抵抗だけを考慮した場合、円軌道の軌道寿命はKingの公式で近似できます。
$$ T_{\text{life}} \approx \frac{H}{v_c}\frac{B_{\text{coeff}}}{\rho_0}\exp\left(\frac{h_0 – h_{\text{initial}}}{H}\right) $$
ただしこの近似は大雑把なもので、実際にはシミュレーションによる数値的な推定がより信頼できます。
弾道係数の影響
弾道係数 $B = m/(C_D A)$ は軌道寿命を大きく左右する重要なパラメータです。以下の比較で、異なる弾道係数を持つ衛星の軌道進化を可視化します。
# --- 弾道係数の違いによる軌道寿命の比較 ---
B_values = [25.0, 50.0, 100.0, 200.0] # [kg/m^2]
colors = ['red', 'orange', 'blue', 'green']
fig, ax = plt.subplots(figsize=(10, 6))
for B_val, color in zip(B_values, colors):
# 簡易的な軌道減衰モデル（円軌道近似）
h_current = 600.0 # [km]
dt_days = 1.0
t_list = [0.0]
h_list = [h_current]
t_total = 0.0
max_time = 30 * 365.25 # 最大30年
while h_current > 120 and t_total < max_time:
r_current = (R_E + h_current * 1e3)
v_circ = np.sqrt(mu / r_current)
rho = atmospheric_density(h_current)
# 1軌道あたりの半長軸減少量（円軌道近似）
da_rev = -2 * np.pi * rho * r_current**2 / B_val
T_orbit = 2 * np.pi * r_current / v_circ # 軌道周期 [s]
# 1日あたりの変化
n_rev_day = 86400.0 / T_orbit
da_day = da_rev * n_rev_day # [m/day]
h_current += da_day / 1e3 # [km]
t_total += dt_days
t_list.append(t_total / 365.25)
h_list.append(h_current)
ax.plot(t_list, h_list, color=color, linewidth=1.5,
label=f'B = {B_val:.0f} kg/m²')
ax.set_xlabel('Time [years]')
ax.set_ylabel('Altitude [km]')
ax.set_title('弾道係数の違いによる軌道減衰の比較（初期高度600km）')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_ylim(100, 650)
plt.tight_layout()
plt.savefig('ballistic_coefficient_comparison.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、弾道係数が軌道寿命に与える影響が明確に読み取れます。
-
弾道係数が小さい（面積質量比が大きい）衛星ほど早く減衰する — $B = 25$ kg/m² の衛星（たとえば薄い大面積デブリ）は数年で再突入しますが、$B = 200$ kg/m² の密実な衛星は数十年持ちます。
-
降下は加速的 — 最初はゆっくり降下していますが、高度が下がるにつれて大気密度が急増するため、降下速度が加速します。最終段階では数日〜数週間で一気に高度が下がります。
-
「膝」の存在 — 高度約300 km付近でグラフの傾きが急変します。この高度以下では大気密度が急激に増し、軌道の崩壊が不可避的になります。
これらの知見は、衛星の設計（面積質量比の最適化）やデブリの長期追跡において実用的な重要性を持ちます。
まとめ
本記事では、摂動下における軌道の長期伝播の理論と実装を解説しました。
- 短期予測と長期予測の違い: 短期予測は位置精度を重視し、長期予測は軌道要素のトレンドを重視する。それぞれに適した手法が異なる
- ガウスの摂動方程式: 摂動加速度をRSW方向に分解し、各軌道要素の時間変化率を記述する。接線方向の力が最もエネルギーを効率的に変化させる
- 主要な摂動源: J2（地球扁平率）、大気抵抗、第三体摂動（月・太陽）、太陽輻射圧の4つが主要。低軌道では大気抵抗とJ2、高軌道では第三体摂動が支配的
- 数値積分法: DOP853（Dormand-Prince 8次法）が長期積分に適している。誤差許容値の設定とエネルギー保存の監視が重要
- 軌道寿命: 弾道係数と太陽活動が軌道寿命を大きく左右する。高度600 kmの衛星でも弾道係数次第で寿命は数年〜数十年まで変わる
次のステップとして、以下の記事も参考にしてください。