衛星の軌道減衰と再突入予測 — 大気抵抗モデルから寿命を推定する

2022年2月、SpaceXが打ち上げたばかりのStarlink衛星49機のうち40機が、打ち上げ直後に大気圏に再突入して失われました。原因は地磁気嵐による大気膨張です。太陽からの荷電粒子が地球の磁気圏に流入し、大気が加熱・膨張して高度200 km付近の大気密度が通常の約50%も増加しました。軌道投入直後の低い高度にいたStarlink衛星群は、想定を超える大気抵抗を受けて軌道を維持できなくなったのです。

この事例が示すように、「衛星はなぜ落ちるのか」「いつ落ちるのか」という問いは、宇宙ミッションの計画から宇宙デブリの管理まで、極めて実用的な問題です。衛星の軌道減衰と再突入予測は、大気抵抗のモデル化と密接に結びついています。

軌道減衰と再突入予測を理解すると、以下の応用が可能になります。

• 軌道寿命の設計: コンステレーション衛星の運用寿命と自然減衰による廃棄を計画する
• 再突入リスク評価: 制御不能となった大型衛星やロケット上段の落下範囲を予測する
• ミッション後廃棄の25年ルール: 国際ガイドラインに基づくデオービット計画の策定

本記事の内容

• 衛星はなぜ落ちるか — 大気抵抗の物理
• 大気密度モデル（指数大気、NRLMSISE-00）
• 弾道係数と軌道減衰の関係
• Kingの公式（軌道寿命の解析近似）
• 再突入高度と空力加熱
• Pythonで軌道減衰シミュレーション
• ISS/Starlinkの軌道維持戦略
• 太陽活動と大気膨張の影響

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 二体問題の理論と解法
• 軌道要素の物理的意味と計算
• 軌道の長期伝播

衛星はなぜ落ちるか — 大気抵抗の物理

低軌道にも大気はある

地球の大気は地表から宇宙に向かって連続的に薄くなっていきますが、明確な「境界」はありません。一般に高度100 km（カーマンライン）を宇宙との境界とする慣習がありますが、これは空力揚力だけでは飛行できなくなる高度に基づく定義であり、大気が完全になくなるわけではありません。

実際には、高度400 km（ISS軌道）でも大気密度は $\sim 10^{-12}$ kg/m³ 程度あります。地表の密度 $\sim 1.2$ kg/m³ に比べれば12桁も小さいですが、衛星が時速28,000 km（秒速7.7 km）で飛行していることを考えると、この希薄な大気でも無視できない抗力を生みます。

イメージとしては、空気抵抗を受けながら坂を下りるボールに似ています。最初は坂の傾斜が緩やかなのでゆっくり加速しますが、高度が下がるほど空気が濃くなり、ブレーキが強くなります。しかし同時に、速度を失った衛星はさらに高度が下がるため、大気はさらに濃くなる — この正のフィードバックによって、最終段階では急激に高度が低下します。

大気抵抗の基本式

衛星に作用する大気抵抗（ドラッグ）は、流体力学における抗力の式で表されます。

$$ \bm{F}_{\text{drag}} = -\frac{1}{2}\rho v_{\text{rel}}^2 C_D A \hat{\bm{v}}_{\text{rel}} $$

ここで

• $\rho$: 大気密度 [kg/m³]
• $v_{\text{rel}}$: 大気に対する衛星の相対速度 [m/s]
• $C_D$: 抗力係数（無次元、通常2.0〜2.5）
• $A$: 衛星の進行方向に対する投影断面積 [m²]
• $\hat{\bm{v}}_{\text{rel}}$: 相対速度の単位ベクトル

符号がマイナスなのは、抗力が常に相対速度と逆向き（衛星を減速させる方向）に作用するからです。

衛星の加速度にすると

$$ \bm{a}_{\text{drag}} = -\frac{1}{2}\rho v_{\text{rel}}^2 \frac{C_D A}{m}\hat{\bm{v}}_{\text{rel}} = -\frac{1}{2}\frac{\rho v_{\text{rel}}^2}{B}\hat{\bm{v}}_{\text{rel}} $$

ここで $B = m/(C_D A)$ は弾道係数です。弾道係数が大きい（重くて小さい）衛星は大気抵抗の影響を受けにくく、弾道係数が小さい（軽くて大きい）衛星は影響を強く受けます。

1軌道あたりのエネルギー損失

円軌道を仮定すると、衛星の力学的エネルギーは $\mathcal{E} = -\mu/(2a)$ です。大気抵抗による1軌道あたりのエネルギー損失は次のように近似できます。

$$ \Delta \mathcal{E} = \oint \bm{F}_{\text{drag}} \cdot d\bm{s} \approx -\pi \rho a \frac{C_D A}{m} v^2 a $$

ここで円軌道速度 $v = \sqrt{\mu/a}$ を使うと、1軌道あたりの半長軸の変化量が得られます。

$$ \frac{d\mathcal{E}}{da} = \frac{\mu}{2a^2} $$

したがって

$$ \Delta a \approx -2\pi \rho \frac{a^2}{B} \cdot v^2 \cdot \frac{a}{\mu} \cdot \frac{a^2}{\mu} \cdot \mu $$

正確に計算すると

$$ \Delta a_{\text{rev}} = -2\pi \rho \frac{a^2}{B} $$

この式は非常にシンプルで、1周あたりの半長軸の減少量は大気密度 $\rho$、軌道半径の2乗 $a^2$、弾道係数 $B$ の逆数で決まります。

大気抵抗の物理がわかったところで、次に大気密度をどのようにモデル化するかを見ていきましょう。

大気密度モデル

指数大気モデル

最も単純な大気モデルは、密度が高度に対して指数関数的に減少するという仮定です。

$$ \rho(h) = \rho_0 \exp\left(-\frac{h – h_0}{H}\right) $$

ここで $h_0$ は参照高度、$\rho_0$ はその高度での密度、$H$ はスケールハイト（密度が $1/e$ に減少する高度差）です。

スケールハイト $H$ は高度とともに増加します。これは、高高度ほど大気の組成が軽い元素（ヘリウム、水素）に変わり、温度も変化するためです。以下に代表的な値を示します。

指数大気モデルは単純ですが、大まかな軌道寿命推定には十分な精度を持ちます。特に概算や教育目的では広く使われます。

NRLMSISE-00モデル

実用的な軌道予測には、より精密な大気モデルが必要です。NRLMSISE-00（Naval Research Laboratory Mass Spectrometer and Incoherent Scatter radar Extended model）は、以下の要因を考慮した経験的大気モデルです。

• 高度: 地表から1000 km以上まで対応
• 緯度・経度: 昼夜や緯度による密度変動
• 太陽活動: F10.7太陽電波フラックス（10.7 cm波長の電波強度）で指標化
• 地磁気活動: Ap指数で磁気嵐の影響を考慮
• 季節変動: 夏半球の大気膨張など

太陽活動の影響は特に大きく、高度400 kmでの大気密度は太陽活動極大期に極小期の約10倍に達することがあります。このため、軌道寿命の予測では太陽活動の予報が重要な不確実性要因となります。

大気の共回転

もう一つ注意すべき点は、地球の大気は地球の自転と一緒に回転しているということです。赤道上の低軌道衛星の場合、大気の共回転速度は約465 m/sです。衛星の軌道速度（約7.7 km/s）に比べれば小さいですが、正確な抗力計算には考慮が必要です。

大気に対する衛星の相対速度は

$$ \bm{v}_{\text{rel}} = \bm{v}_{\text{sat}} – \bm{\omega}_E \times \bm{r}_{\text{sat}} $$

ここで $\bm{\omega}_E = (0, 0, 7.292 \times 10^{-5})$ rad/s は地球の自転角速度ベクトルです。

大気モデルが理解できたところで、次に弾道係数が軌道寿命にどう影響するかを定量的に見ていきましょう。

弾道係数と軌道寿命

弾道係数の定義と代表値

弾道係数 $B = m/(C_D A)$ は、衛星の大気抵抗に対する「頑丈さ」を表す指標です。単位はkg/m²で、大きいほど大気抵抗に強い（落ちにくい）ことを意味します。

典型的な衛星の弾道係数を比較します。

Starlink衛星は薄い平板状の設計のため弾道係数が非常に小さく、大気抵抗の影響を強く受けます。これは運用中は軌道維持のためにイオンスラスタで高度を維持する一方、ミッション後は自然減衰で速やかに再突入できるという利点にもなります。

軌道寿命の高度依存性

同じ弾道係数の衛星でも、初期軌道高度によって軌道寿命は桁違いに変わります。高度が100 km上がるごとに大気密度は約1桁下がるため、軌道寿命は指数関数的に増大します。

大まかな目安として、$B = 50$ kg/m² の衛星の場合

高度800 km以上のデブリは自然減衰では何世紀もかかるため、デブリ問題において特に深刻です。国際ガイドラインでは、LEO衛星のミッション後25年以内の再突入が推奨されています（近年は5年に短縮する動きがあります）。

軌道寿命の定量的な推定には、解析的な近似公式が便利です。次にKingの公式を紹介します。

Kingの公式 — 軌道寿命の解析近似

公式の導出の考え方

円軌道の衛星が大気抵抗のみで減衰する場合、半長軸の時間変化率は

$$ \frac{da}{dt} = -2\pi \rho(a) \frac{a^2}{B} \cdot \frac{1}{T_{\text{orb}}} = -\frac{\rho(a) a v_c}{B} $$

ここで、円軌道速度 $v_c = \sqrt{\mu/a}$ を使いました。これに指数大気 $\rho = \rho_0 \exp(-(a – a_0)/H)$ を代入すると

$$ \frac{da}{dt} = -\frac{\rho_0}{B}\exp\left(-\frac{a – a_0}{H}\right) a \sqrt{\frac{\mu}{a}} $$

この微分方程式を近似的に解くと、Kingの公式が得られます。

Kingの公式

円軌道の衛星の軌道寿命は次の近似式で推定できます。

$$ T_{\text{life}} \approx \frac{H B}{\rho_0 a_0 v_{c,0}} $$

ここで $a_0$、$v_{c,0}$、$\rho_0$ は初期軌道での値です。この式は次のように解釈できます。

• $H$ が大きい（スケールハイトが大きい）と、高度が下がっても密度の増加が緩やかなため、寿命が延びる
• $B$ が大きい（弾道係数が大きい）と、同じ密度でも受ける減速が小さいため、寿命が延びる
• $\rho_0$ が大きい（初期高度での密度が高い）と、寿命は短くなる

公式の精度と限界

Kingの公式は概算には便利ですが、以下の限界があります。

1. 円軌道を仮定 — 楕円軌道の場合、大気抵抗は主に近地点で作用するため、減衰パターンが異なる
2. 一定のスケールハイト — 実際にはスケールハイトは高度とともに変化する
3. 太陽活動の変動を無視 — 11年周期の太陽活動変動は軌道寿命に大きく影響する
4. J2摂動を無視 — J2による近地点高度の振動は、大気抵抗の効き方に影響する

実用的にはKingの公式でオーダー推定を行い、精密な寿命予測には数値シミュレーションを使うという使い分けが一般的です。

解析的な近似の限界を理解したところで、次にPythonで数値的な軌道減衰シミュレーションを実装します。

Pythonで軌道減衰シミュレーション

円軌道近似による高速シミュレーション

まず、円軌道を仮定した簡易モデルで様々な条件での軌道寿命を計算します。このモデルは1軌道ごとの平均的な減衰を計算するため、数値的に非常に効率的です。

import numpy as np
import matplotlib.pyplot as plt

# --- 定数 ---
mu = 3.986004418e14      # 地球の重力パラメータ [m^3/s^2]
R_E = 6378.137e3          # 地球半径 [m]

# --- 指数大気モデル ---
# (基準高度[km], 密度[kg/m^3], スケールハイト[km])
atm_table = [
    (150, 2.070e-9, 22.523),
    (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, f107_factor=1.0):
    """
    指数大気モデルによる密度推定
    f107_factor: 太陽活動補正（1.0=標準, >1で活発期）
    """
    if h_km > 1200 or h_km < 100:
        return 0.0
    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]
            rho = rho0 * np.exp(-(h_km - h0) / H)
            return rho * f107_factor
    h0, rho0, H = atm_table[-1]
    return rho0 * np.exp(-(h_km - h0) / H) * f107_factor

def simulate_decay_circular(h_init_km, B_coeff, f107_factor=1.0,
                            max_years=100, h_reentry=120.0):
    """
    円軌道近似による軌道減衰シミュレーション

    Parameters:
        h_init_km: 初期高度 [km]
        B_coeff: 弾道係数 [kg/m^2]
        f107_factor: 太陽活動補正係数
        max_years: 最大計算年数
        h_reentry: 再突入判定高度 [km]

    Returns:
        t_years: 時間配列 [年]
        h_km: 高度配列 [km]
    """
    h = h_init_km  # [km]
    dt = 0.01      # 時間刻み [日]
    t = 0.0
    max_days = max_years * 365.25

    t_list = [0.0]
    h_list = [h]

    while h > h_reentry and t < max_days:
        r = R_E + h * 1e3          # 軌道半径 [m]
        v_c = np.sqrt(mu / r)       # 円軌道速度 [m/s]
        T_orb = 2 * np.pi * r / v_c # 軌道周期 [s]
        rho = atmospheric_density(h, f107_factor)

        if rho < 1e-20:
            # 大気が事実上なければ大きなステップで飛ばす
            t += 30.0
            t_list.append(t / 365.25)
            h_list.append(h)
            continue

        # 1軌道あたりの半長軸減少量 [m]
        da_rev = -2 * np.pi * rho * r**2 / B_coeff

        # 1日あたりの周回数
        n_rev = 86400.0 / T_orb

        # 1日あたりの高度変化 [km]
        dh_day = da_rev * n_rev / 1e3

        # 適応的な時間刻み（急降下時は細かく）
        if abs(dh_day) > 1.0:
            dt = 0.001
        elif abs(dh_day) > 0.01:
            dt = 0.1
        else:
            dt = 1.0

        h += dh_day * dt
        t += dt

        if t / 365.25 - t_list[-1] > 0.001 or h < h_reentry + 50:
            t_list.append(t / 365.25)
            h_list.append(max(h, 0))

    return np.array(t_list), np.array(h_list)

初期高度による軌道寿命の比較

まず、異なる初期高度での軌道減衰を比較します。

# --- 初期高度の比較 ---
fig, ax = plt.subplots(figsize=(10, 6))

B = 50.0  # 弾道係数 [kg/m^2]（典型的な衛星）
heights = [300, 400, 500, 600, 700, 800]
colors = plt.cm.viridis(np.linspace(0.2, 0.9, len(heights)))

for h0, color in zip(heights, colors):
    t, h = simulate_decay_circular(h0, B, max_years=50)
    label = f'h₀ = {h0} km'
    if h[-1] <= 120:
        # 再突入した場合、寿命を表示
        life = t[-1]
        label += f' (寿命 ≈ {life:.1f} yr)'
    ax.plot(t, h, color=color, linewidth=1.5, label=label)

ax.set_xlabel('Time [years]')
ax.set_ylabel('Altitude [km]')
ax.set_title(f'初期高度による軌道減衰の比較 (B = {B:.0f} kg/m²)')
ax.legend(loc='upper right', fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_ylim(100, 850)
plt.tight_layout()
plt.savefig('orbital_decay_heights.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、軌道減衰の特徴的なパターンが明確に見えます。

1. 高度300 kmでは数ヶ月で再突入する — 大気密度が高いため、急速に高度が低下します。低軌道ミッションでは定期的な軌道維持噴射が不可欠です。

2. 高度が上がるにつれて寿命は指数関数的に延びる — 高度600 kmなら数年、800 kmなら数十年以上です。

3. 最終段階の急降下 — いずれの高度でも、再突入直前に高度が急激に下がる「膝」が見られます。これは高度低下→密度増加→さらなる高度低下という正のフィードバックによるものです。

太陽活動の影響

太陽活動は大気密度に大きな影響を与えます。太陽活動が活発な時期（F10.7指数が高い時期）には、太陽からの紫外線放射が増加し、上層大気が加熱されて膨張します。

# --- 太陽活動の影響 ---
fig, ax = plt.subplots(figsize=(10, 6))

h0 = 500  # 初期高度 [km]
B = 50.0

# 太陽活動レベル
solar_levels = [
    (0.3, '太陽極小期 (F10.7≈70)', 'blue'),
    (1.0, '標準 (F10.7≈150)', 'green'),
    (3.0, '太陽極大期 (F10.7≈200)', 'orange'),
    (8.0, '太陽嵐時 (極端)', 'red'),
]

for f107, label, color in solar_levels:
    t, h = simulate_decay_circular(h0, B, f107_factor=f107, max_years=30)
    life_str = f' (寿命 ≈ {t[-1]:.1f} yr)' if h[-1] <= 120 else ''
    ax.plot(t, h, color=color, linewidth=1.5, label=label + life_str)

ax.set_xlabel('Time [years]')
ax.set_ylabel('Altitude [km]')
ax.set_title(f'太陽活動レベルによる軌道寿命の変化 (h₀={h0} km, B={B:.0f} kg/m²)')
ax.legend(loc='upper right', fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_ylim(100, 550)
plt.tight_layout()
plt.savefig('solar_activity_effect.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、太陽活動が軌道寿命に与える劇的な影響がわかります。

1. 太陽極大期と極小期で寿命が数倍〜10倍異なる — 同じ衛星でも、どの太陽活動周期に運用されるかで寿命が大きく変わります。コンステレーション計画では太陽活動の予測が不可欠です。

2. 太陽嵐のような極端なイベントは急激な軌道降下を引き起こす — 2022年のStarlink事例のように、短期間に大量のデブリが生成されたり、低軌道の衛星が失われたりするリスクがあります。

3. 不確実性が大きい — 太陽活動の長期予測は困難で、軌道寿命の予測にも大きな不確実性が伴います。ミッション設計では保守的な仮定（太陽極大期を想定）が推奨されます。

太陽活動の影響の次に、衛星が実際に再突入する際の物理を見てみましょう。

再突入高度と空力加熱

再突入のフェーズ

衛星の再突入は、一般に以下の3つのフェーズに分けられます。

フェーズ1: 自由分子流領域（高度120 km以上） — 大気の分子の平均自由行程が衛星の大きさより十分大きい領域です。個々の空気分子が衛星表面に衝突して運動量を受け渡します。クヌーセン数 $Kn = \lambda/L \gg 1$（$\lambda$: 平均自由行程、$L$: 衛星の代表長さ）で特徴づけられます。

フェーズ2: 遷移領域（高度80〜120 km） — 自由分子流から連続体流れへの遷移領域です。空力加熱が急速に増大し、衛星の表面温度が上昇します。

フェーズ3: 連続体流れ領域（高度80 km以下） — 通常の空気力学が成り立つ領域です。衝撃波が形成され、衛星前面の空気は数千〜数万度に加熱されます。

空力加熱

再突入時の最大加熱率は、衛星の速度と大気密度に依存します。ストークス点（衝撃波の直後、衛星の前面中心）での加熱率は、次の近似式で表されます。

$$ \dot{q} \approx C \sqrt{\frac{\rho}{R_n}} v^3 $$

ここで $C$ は定数（$\approx 1.83 \times 10^{-4}$ W·s³/(kg^{1/2}·m))、$R_n$ は衛星先端の曲率半径、$v$ は衛星速度です。

軌道速度 $v \approx 7.7$ km/s での加熱は極めて過酷で、多くの衛星構造材料の融点を超えます。制御されない再突入では、衛星は高度70〜80 km付近で分解し始め、大部分は燃え尽きますが、チタン合金やステンレス鋼などの高融点材料は地表に到達することがあります。

再突入の予測精度

再突入時期の予測精度は、再突入が近づくにつれて向上します。典型的な精度は以下の通りです。

この大きな不確実性は、主に大気密度の予測困難さに起因します。太陽活動や地磁気嵐のタイミングによって大気密度が急変するため、再突入の正確な時刻と場所の予測は本質的に困難です。

再突入の物理を理解したところで、次に実際の衛星（ISSとStarlink）がどのように軌道維持を行っているかを見てみましょう。

ISS/Starlinkの軌道維持戦略

ISSの軌道維持（リブースト）

ISSは高度約400 kmを周回しており、大気抵抗により年間約2 km/yearの高度低下が生じます（太陽活動に依存）。これを補うため、定期的にリブースト（軌道引き上げ噴射）が行われます。

ISSのリブースト戦略は以下の通りです。

• 頻度: 約2ヶ月に1回（太陽極大期にはより頻繁）
• ΔV: 1回あたり約1〜2 m/s
• 手段: プログレス補給船やCygnus補給船のスラスタ、またはZvezdaモジュールのエンジン
• 年間ΔV: 約10〜20 m/s（太陽活動依存）

ISSの大きな太陽電池パネルは弾道係数を小さくするため、大気抵抗の影響が大きいです。太陽活動極大期には月に1回以上のリブーストが必要になることもあります。

Starlinkの軌道維持と廃棄

SpaceXのStarlink衛星は、運用軌道（約550 km）での軌道維持にイオン推進（クリプトンイオンスラスタ）を使用しています。

Starlink衛星の軌道戦略は以下のようなものです。

1. 初期軌道投入: ファルコン9で高度約300 km付近に投入
2. 軌道上昇: イオンスラスタで約550 kmまで自力で上昇（数週間〜数ヶ月）
3. 運用軌道維持: 高度550 kmで軌道維持
4. ミッション後廃棄: イオンスラスタで高度を下げ、大気抵抗で自然減衰

Starlinkの弾道係数が小さい（約10 kg/m²）ことは、故障した衛星が高度550 kmから数年で自然に再突入することを意味しています。これは意図的な設計で、デブリの長期蓄積を防ぐためです。

# --- ISS と Starlink の軌道減衰比較 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# ISS（リブーストなしの場合）
t_iss, h_iss = simulate_decay_circular(420, B_coeff=76, max_years=5)
axes[0].plot(t_iss, h_iss, 'b-', linewidth=1.5)
axes[0].set_xlabel('Time [years]')
axes[0].set_ylabel('Altitude [km]')
axes[0].set_title('ISS（リブーストなしの場合）')
axes[0].set_ylim(100, 450)
axes[0].grid(True, alpha=0.3)
axes[0].axhline(y=120, color='red', linestyle='--', alpha=0.5,
                label='Re-entry altitude')
axes[0].legend()

# Starlink（推進なしの場合）
t_sl300, h_sl300 = simulate_decay_circular(300, B_coeff=10, max_years=3)
t_sl550, h_sl550 = simulate_decay_circular(550, B_coeff=10, max_years=10)

axes[1].plot(t_sl300, h_sl300, 'r-', linewidth=1.5, label='h₀ = 300 km（投入直後）')
axes[1].plot(t_sl550, h_sl550, 'b-', linewidth=1.5, label='h₀ = 550 km（運用軌道）')
axes[1].set_xlabel('Time [years]')
axes[1].set_ylabel('Altitude [km]')
axes[1].set_title('Starlink（推進なしの場合）')
axes[1].set_ylim(100, 600)
axes[1].grid(True, alpha=0.3)
axes[1].axhline(y=120, color='red', linestyle='--', alpha=0.5,
                label='Re-entry altitude')
axes[1].legend()

plt.tight_layout()
plt.savefig('iss_starlink_decay.png', dpi=150, bbox_inches='tight')
plt.show()

この比較グラフから、2つの衛星の軌道特性の違いが明確にわかります。

1. ISS（左図）: リブーストなしでは約1〜2年で再突入します。弾道係数は約76 kg/m²と決して小さくありませんが、高度400 kmの大気密度でこの大きな構造体には相当な抗力が作用します。ISSの運用にはリブーストが不可欠です。

2. Starlink（右図）: 投入高度300 kmでは数ヶ月で再突入します（2022年の事例）。運用軌道550 kmでも数年で再突入するため、故障した衛星もデブリとして長期間残留しません。これはStarlinkの小さな弾道係数による「設計された使い捨て性」です。

軌道寿命マップの作成

最後に、初期高度と弾道係数をパラメータとした軌道寿命マップを作成します。

# --- 軌道寿命マップ ---
h_range = np.arange(200, 901, 25)  # [km]
B_range = np.array([10, 20, 50, 100, 200])  # [kg/m^2]

lifetime_map = np.zeros((len(B_range), len(h_range)))

for i, B in enumerate(B_range):
    for j, h0 in enumerate(h_range):
        t, h = simulate_decay_circular(h0, B, max_years=200)
        if h[-1] <= 120:
            lifetime_map[i, j] = t[-1]
        else:
            lifetime_map[i, j] = 200  # 200年以上

fig, ax = plt.subplots(figsize=(10, 6))

for i, B in enumerate(B_range):
    ax.semilogy(h_range, lifetime_map[i, :], linewidth=2,
                label=f'B = {B} kg/m²')

ax.set_xlabel('Initial Altitude [km]')
ax.set_ylabel('Orbital Lifetime [years]')
ax.set_title('軌道寿命マップ（円軌道, 標準大気）')
ax.legend()
ax.grid(True, alpha=0.3, which='both')
ax.axhline(y=25, color='red', linestyle='--', alpha=0.5,
           label='25-year guideline')
ax.set_ylim(0.01, 300)
plt.tight_layout()
plt.savefig('orbital_lifetime_map.png', dpi=150, bbox_inches='tight')
plt.show()

この軌道寿命マップから、ミッション設計に直結する重要な情報が読み取れます。

1. 25年ルールの閾値高度 — 赤い破線（25年）と各曲線の交点が、自然減衰で25年以内に再突入する最大高度を示しています。$B = 50$ kg/m² の場合、約650 km以下であれば25年ルールを自然に満たせます。

2. 弾道係数の影響は高高度ほど大きい — 低軌道（300 km以下）ではどの弾道係数でもすぐに再突入しますが、600 km以上では弾道係数による差が数桁にまで広がります。

3. 800 km以上は深刻 — ほとんどの衛星/デブリにとって、高度800 km以上からの自然減衰は100年以上かかります。この高度帯のデブリ除去には能動的なデオービット（ドラッグセイル等）が必要です。

まとめ

本記事では、衛星の軌道減衰メカニズムと再突入予測の理論・実装を解説しました。

• 大気抵抗の物理: 低軌道にも希薄ながら大気が存在し、秒速7.7 kmで飛行する衛星に対して無視できない抗力を生む。1軌道あたりの半長軸減少は $\Delta a \approx -2\pi\rho a^2/B$ で近似される
• 大気密度モデル: 指数大気モデルが基本。実用的にはNRLMSISE-00が使われ、太陽活動（F10.7）と地磁気活動（Ap）を考慮する
• 弾道係数: $B = m/(C_DA)$ が軌道寿命を決定する鍵パラメータ。軽くて大きい衛星ほど早く減衰する
• Kingの公式: 円軌道の軌道寿命を $T \approx HB/(\rho_0 a_0 v_{c,0})$ で概算できる。精密な予測には数値シミュレーションが必要
• 太陽活動の影響: 太陽極大期と極小期で軌道寿命が数倍〜10倍異なる。11年周期の太陽活動予測がミッション設計の不確実性要因
• ISS/Starlink: ISSは定期的リブーストで軌道維持。Starlinkは低弾道係数を活かした自然減衰による廃棄設計

次のステップとして、以下の記事も参考にしてください。

• 軌道の長期伝播
• 衛星衝突回避の理論
• TLEとSGP4による軌道計算