火星から地球に「Touchdown confirmed」という1行のテレメトリが届くまでに、約14分の片道光遅延を挟んで、機体側ではたった7分間の出来事が完結しています。秒速 $5.4 \text{ km/s}$（マッハ約24）で薄い大気に突入した1トン級の探査機を、地表でぴたりと静止させる。地球からの介入は不可能、すべては事前にプログラムされた自律制御で進む。NASA/JPLが「Seven Minutes of Terror（7分間の恐怖）」と名付けたこの一連のシーケンスが、火星探査の最大関門 EDL（Entry, Descent, Landing） です。
なぜそれほど難しいのか。火星大気は地球の約 $1/100$ という極端な薄さでありながら、超音速領域ではそれなりの空力抵抗を生み、しかも音速の壁を抜けるあたりで衝撃波と乱流が密接に絡み合います。「大気が薄すぎてパラシュートだけでは止まりきれない」、しかし同時に「大気はあるので空力加熱が無視できない」——この中途半端さこそが、地球用の弾道再突入カプセルや月用のロケット直下降ではそのまま流用できない理由です。歴代の火星着陸はエアバッグ、レトロロケット、そしてスカイクレーン（Sky-Crane）という、奇抜にも見える三つの戦略を独自に編み出してきました。
本記事では、火星EDLの物理を「何が一番怖いのか」を軸に段階別に解きほぐします。応用先として想定するのは、(1) 将来の有人火星ミッションで必要になるであろう数トン級の大型ペイロード輸送（既存のスカイクレーンは1トン級が上限）、(2) 月・火星・タイタンなど 大気がある他天体への着陸計画の比較設計、(3) 軌道力学コースで学ぶ突入軌道方程式の数値解法（scipy.integrate.solve_ivp の典型例）です。Mars 2020 のPerseveranceがどう成功し、ExoMars 2016のSchiaparelliがどう失敗したのか、その差はどこにあったのか。Pythonで突入軌道を積分しながら追体験していきましょう。
本記事の内容
- 火星大気の特殊性（地球の1%、しかし熱流束は意外と大きい理由）
- 突入軌道方程式（球座標版）と Sutton-Graves による加熱率の見積もり
- ピーク加熱とピーク減速のタイミング、そして両者がほぼ同時に来る理由
- 超音速パラシュート展開条件（Disk-Gap-Band、$M < 2.2$、$\bar{q} < 850 \text{ Pa}$）
- スカイクレーン Powered Descent と Terrain Relative Navigation（TRN）
- エアバッグ（Pathfinder/MER）、レトロロケット（Viking/Phoenix/InSight）、スカイクレーン（MSL/Mars 2020）の使い分け
- Schiaparelli 2016 がなぜ墜落したか — IMUサチュレーションから連鎖した論理エラー
- Pythonで突入軌道、ピーク加熱タイミング、突入角sensitivity、MarsGRAMライクな大気モデル可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
火星大気はなぜ「中途半端」なのか
地表気圧と密度
地球の地表気圧は約 $101 \text{ kPa}$、対する火星はわずか $600 \text{ Pa}$ 前後（季節と地形高度で $400$–$900 \text{ Pa}$ 程度変動）。比率にして $0.6\%$、ざっくり地球の100分の1です。「これなら抵抗が小さくて楽そう」と思いきや、突入時の機体は秒速数 km という超高速で、しかも質量は数百 kg から1トン超。空力抵抗の絶対量は、薄い大気の中でも莫大になります。
具体的に、抗力の力学的等価指標である動圧 $\bar{q} = \frac{1}{2}\rho V^2$ を考えると、地球低軌道からの再突入時のピーク動圧が $\sim 50 \text{ kPa}$ なのに対し、火星突入では $\sim 5\text{–}10 \text{ kPa}$。1桁小さいですが、ゼロにはほど遠い。むしろ「減速のためにこの動圧を最大限利用しなければならない」のが火星です。地球なら下層大気で確実に止まりますが、火星では地表に達するまでに使える「抵抗予算」が圧倒的に少ない。これが第一の困難です。
大気の薄さが招く本当の問題
ここで火星EDLの本質的な難しさが浮かび上がります。突入時の運動エネルギーは
$$ E_K = \frac{1}{2} m V^2 = \frac{1}{2} \times 1000 \times (5400)^2 \approx 1.46 \times 10^{10} \text{ J} $$
これは TNT 換算で約3.5トン分。これを「7分間」で「薄い大気」を使ってゼロまで持っていく必要があります。地球の場合、低軌道からの再突入では大気圏内を約 $30$ 分かけて滞空し、ピーク動圧 $50 \text{ kPa}$ の濃い空気で減速できます。火星では時間は1/4、大気は1/10。差し引きで「減速能力が圧倒的に足りない」。
このため、火星EDLでは超音速パラシュートを含めても $\sim 90\%$ までしか減速できず、最後の数百 m/s は必ずロケットエンジンによるアクティブ減速（Powered Descent）に頼る必要があります。地球の弾道カプセル（ソユーズ等）はパラシュートだけで地表まで安全に届きますが、火星では物理的にそれが不可能。これが「火星EDLが地球用カプセルの流用では絶対に解けない」根本理由です。
しかし加熱は地球並みかそれ以上
「大気が薄いのだから加熱も少ないはず」と思うのは早計です。突入速度が高ければ加熱はその3乗に比例して増えます（Sutton-Graves式、後述）。Perseveranceの場合、突入速度 $V_e \approx 5.4 \text{ km/s}$、ピーク加熱率は約 $130\text{–}160 \text{ W/cm}^2$。これは地球低軌道再突入（$\sim 100\text{–}200 \text{ W/cm}^2$）と同レベルです。「薄い」と「加熱しない」は別問題。これが第二の困難で、火星専用の耐熱シールド（PICA等のアブレータ）が常に必要になります。
それでは、この「薄いのに熱い、しかし減速予算が足りない」大気を、機体がどう乗りこなすのか。突入軌道の運動方程式から組み立てていきましょう。
突入軌道方程式
座標系と仮定
火星を半径 $R_M \approx 3389.5 \text{ km}$ の球体、自転は無視（突入時間が7分と短いため影響は小さい）、機体は質点として扱います。位置を球座標 $(r, \theta)$、速度を地表に対する大きさ $V$ と飛行経路角 $\gamma$（水平を $0$、上向きを正）で表します。ここで $\gamma$ が正なら上昇、負なら降下です。突入は典型的に $\gamma_e \approx -12^\circ$ から $-15^\circ$ で開始します。
機体に働く力は (1) 重力 $mg$、(2) 抗力 $D = \frac{1}{2}\rho V^2 C_D A$、(3) 揚力 $L = \frac{1}{2}\rho V^2 C_L A$ の3つ。揚力は機体姿勢（迎角）に依存しますが、Mars 2020のような中程度の揚抗比 $L/D \approx 0.24$ のカプセルでは、揚力の符号制御（bank angle modulation）でレンジを微調整します。
運動方程式の導出
地表から見た半径方向と接線方向の運動方程式を立てます。位置・速度の関係は
$$ \dot{r} = V \sin\gamma, \quad r\dot{\theta} = V \cos\gamma $$
次に、機体の運動量変化を考えます。慣性系での加速度のうち、接線方向（速度ベクトル方向）と法線方向（速度ベクトルに垂直、上向きが正）に分解すると、
$$ \dot{V} = -\frac{D}{m} – g\sin\gamma $$
$$ V\dot{\gamma} = \frac{L\cos\sigma}{m} – \left(g – \frac{V^2}{r}\right)\cos\gamma $$
ここで $\sigma$ はバンク角（揚力ベクトルを進行方向まわりに回転させる角度、$\sigma=0$ で揚力上向き）、$g = \mu_M/r^2$ は局所重力です。$\mu_M = 4.282837 \times 10^{13} \text{ m}^3/\text{s}^2$ は火星の重力定数。
第二式の右辺第二項 $-(g – V^2/r)\cos\gamma$ は遠心力補正項で、$V^2/r$ が重力と釣り合うと（$V \approx 3.5 \text{ km/s}$ 程度）軌道が水平に維持される効果を表します。突入初期 $V \approx 5.4 \text{ km/s}$ ではこの項が負（下向きに引き込む）、終盤ではほぼゼロになります。
大気密度モデル
火星大気は地球の ISA のような綺麗な単純指数則ではありませんが、まずは1次近似として等温大気の指数モデル
$$ \rho(h) = \rho_0 \exp\left(-\frac{h}{H}\right) $$
を使います。$\rho_0 \approx 0.020 \text{ kg/m}^3$（地表）、スケールハイト $H \approx 11.1 \text{ km}$。地球の $H \approx 8.5 \text{ km}$ より厚いのは、火星の重力が地球の38%しかなく、大気が「広がりやすい」ためです。
実運用ではこれに季節変動・ダストストーム・場所依存性を加えた Mars-GRAM モデルが使われますが、定性的な議論には指数モデルで十分です。
これらを連立すれば、初期条件 $(r_0, \theta_0, V_0, \gamma_0)$ から地表 $r = R_M$ までの軌道を数値積分できます。Pythonでは scipy.integrate.solve_ivp の標準的な題材です。次節でまず加熱率の見積もり式を導入してから、実装に進みます。
Sutton-Gravesによるよどみ点加熱率
よどみ点とは
突入機の前面、流れが完全に止まる点（よどみ点, stagnation point）は最も加熱が激しい部位です。ここでの加熱率を $\dot{q}_s$ （単位 $\text{W/cm}^2$）と書きます。流体力学・気体力学の詳細な解析からは複雑な形になりますが、Suttonと Graves が1971年に提案した実用的な簡略式が、半世紀後の今もミッション初期検討で使われています：
$$ \dot{q}_s = K \sqrt{\frac{\rho}{R_N}}\, V^3 $$
ここで $R_N$ はよどみ点の曲率半径 [m]、$\rho$ は局所大気密度 [kg/m³]、$V$ は機体速度 [m/s]。係数 $K$ は大気組成に依存し、火星（$\text{CO}_2$ 主体）では
$$ K_{\mathrm{Mars}} \approx 1.898 \times 10^{-4} \quad [\text{W} \cdot \text{m}^{-1/2} \cdot \text{kg}^{-1/2} \cdot \text{s}^2] $$
（最終的に $\dot{q}_s$ が $\text{W/cm}^2$ になる単位系で扱う場合、$\rho$ を kg/m³、$V$ を m/s、$R_N$ を m として上記係数を使うと $\text{W/m}^2$ になり、$10^{-4}$ をかけて $\text{W/cm}^2$ に変換します）。
地球（$N_2/O_2$）の係数 $K_{\mathrm{Earth}} \approx 1.7415 \times 10^{-4}$ と比べると、$\text{CO}_2$ 大気の火星の方が約9%加熱しやすい。これは $\text{CO}_2$ のモル質量が高く、衝撃波背後の温度が同じ速度でも上がりやすいためです。
なぜ $V^3$ なのか
この $V^3$ 依存性が、火星EDL設計を決定づけます。フライバイ軌道から直接突入する場合（典型 $V_e \approx 6\text{–}7 \text{ km/s}$）と、軌道投入してから降下する場合（$V_e \approx 5.4 \text{ km/s}$）では、たった $1.5$ 倍の速度比でも加熱率は $1.5^3 \approx 3.4$ 倍違います。だからこそPerseveranceは直接突入を選んだのに対し、より大型の機体ではいったん火星周回軌道に乗ってから降下する案が検討されます。
導出を簡単に追いましょう。よどみ点では大気が静止し、運動エネルギー $\frac{1}{2}\rho V^3 \cdot$ 単位面積・単位時間（質量流束 $\rho V$ に運動エネルギー $\frac{1}{2}V^2$ をかけたもの）が衝撃波を通じて熱として流入します。すなわち入射エネルギー流束は $\rho V^3$ に比例。これが対流伝熱で機体表面に渡るとき、境界層厚さがよどみ点曲率半径 $R_N$ に依存することから $1/\sqrt{R_N}$ が現れ、最終的に密度の平方根 $\sqrt{\rho}$ になる（厳密な気体力学的境界層解析の帰結）。結果として
$$ \dot{q}_s \propto \sqrt{\frac{\rho}{R_N}} V^3 $$
という形に収まります。$R_N$ が大きい（鈍い形状）ほど加熱が緩和されるのは、衝撃波が機体から離れて衝撃層が厚くなり、熱が拡散するため。これがアポロカプセル以来「鈍頭体（blunt body）」が突入機の標準形状である理由です。
ピーク加熱とピーク減速の関係
突入時、$\rho$ は高度低下とともに指数関数的に増え、$V$ は減速で減ります。両者の積 $\sqrt{\rho}V^3$ が最大になる点がピーク加熱、$\rho V^2$ が最大になる点がピーク動圧（≈ピーク減速）。簡単な解析（Allen-Eggers 解、定数 $C_D$、定数 $\gamma$ の弾道突入を仮定）では、両者はほぼ同じ高度で発生し、典型的に $30$–$40 \text{ km}$ あたり、突入から約 $80\text{–}100$ 秒後です。Perseveranceの実データもこの予測とよく一致しています。
ここまでで「加熱はどう計算するか」「いつピークが来るか」が見えてきました。では実際に運動方程式を数値積分して、軌道全体を再現してみましょう。
Pythonで突入軌道を積分する
設定とソルバ呼び出し
scipy.integrate.solve_ivp で積分します。状態ベクトルを $\bm{s} = [r, \theta, V, \gamma]^\top$ とし、上述の運動方程式を関数として渡します。
import numpy as np
from scipy.integrate import solve_ivp
# 火星定数
R_M = 3389.5e3 # 火星半径 [m]
MU_M = 4.282837e13 # 重力定数 [m^3/s^2]
# 大気モデル（指数則）
RHO0 = 0.020 # 地表密度 [kg/m^3]
H_SCALE = 11.1e3 # スケールハイト [m]
def rho_mars(h):
"""高度 h [m] における火星大気密度 [kg/m^3]"""
return RHO0 * np.exp(-h / H_SCALE)
# Mars 2020 Perseverance 相当のカプセル諸元
MASS = 3257.0 # 突入質量 [kg]
AREA = np.pi * (4.5 / 2)**2 # 基準面積（直径4.5m）[m^2]
CD = 1.45 # 抗力係数（鈍頭体, 超音速領域）
CL = 0.348 # 揚力係数（L/D≈0.24 に相当）
RN = 1.125 # よどみ点曲率半径 [m]
K_MARS = 1.898e-4 # Sutton-Graves係数（W·m^-1/2·kg^-1/2·s^2）
def dynamics(t, s, bank_sigma=0.0):
"""突入軌道のODE。bank_sigma はバンク角 [rad]（揚力ベクトルの傾き）"""
r, theta, V, gamma = s
h = r - R_M
rho = rho_mars(h)
g = MU_M / r**2
q_bar = 0.5 * rho * V**2 # 動圧
D = q_bar * CD * AREA # 抗力
L = q_bar * CL * AREA # 揚力
dr = V * np.sin(gamma)
dtheta = V * np.cos(gamma) / r
dV = -D / MASS - g * np.sin(gamma)
dgamma = (L * np.cos(bank_sigma) / MASS) / V \
- (g - V**2 / r) * np.cos(gamma) / V
return [dr, dtheta, dV, dgamma]
# 突入インターフェース（高度125km）からスタート
h_entry = 125e3
V_entry = 5400.0 # m/s
gamma_entry = np.deg2rad(-12.0)
s0 = [R_M + h_entry, 0.0, V_entry, gamma_entry]
# 地表到達でイベント終了
def hit_surface(t, s, **_):
return s[0] - R_M
hit_surface.terminal = True
hit_surface.direction = -1
sol = solve_ivp(dynamics, [0, 600], s0, events=hit_surface,
max_step=0.5, rtol=1e-9, atol=1e-9, dense_output=True)
print(f"突入から地表まで: {sol.t[-1]:.1f} 秒")
print(f"地表での速度: {sol.y[2][-1]:.1f} m/s")
print(f"地表での経路角: {np.rad2deg(sol.y[3][-1]):.1f} deg")
このコードを実行すると、おおよそ「突入から $250$ 秒前後で地表」「最終速度数百 m/s」という結果になります。実際の Mars 2020 では超音速パラシュート展開が突入から $240$ 秒で行われたので、純粋な空力減速だけでパラシュート展開条件まで到達するタイミングがリアルに再現されていることがわかります。最終速度が数百 m/s 残るのが「パラシュート + Powered Descentで残りを処理しなければならない」典型例で、火星EDLの厳しさを数値で実感できます。
軌道の可視化と各物理量の時間履歴
積分結果から、高度・速度・動圧・加熱率を時系列でプロットして、ピーク加熱とピーク減速のタイミングを見ます。
import matplotlib.pyplot as plt
t = sol.t
r, theta, V, gamma = sol.y
h = r - R_M
rho = rho_mars(h)
q_bar = 0.5 * rho * V**2 # 動圧 [Pa]
decel = q_bar * CD * AREA / MASS / 9.81 # 減速度 [g]
# Sutton-Graves加熱率：W/m^2 → W/cm^2 で扱う
q_dot = K_MARS * np.sqrt(rho / RN) * V**3 * 1e-4
fig, axes = plt.subplots(2, 2, figsize=(11, 8))
axes[0, 0].plot(t, h / 1000, 'b-')
axes[0, 0].set_xlabel('Time [s]'); axes[0, 0].set_ylabel('Altitude [km]')
axes[0, 0].set_title('Altitude vs time'); axes[0, 0].grid(alpha=0.3)
axes[0, 1].plot(t, V, 'g-')
axes[0, 1].axhline(450, color='r', ls='--', alpha=0.5, label='Mach 2.2 ≈ 450 m/s')
axes[0, 1].set_xlabel('Time [s]'); axes[0, 1].set_ylabel('Velocity [m/s]')
axes[0, 1].set_title('Velocity vs time'); axes[0, 1].legend(); axes[0, 1].grid(alpha=0.3)
axes[1, 0].plot(t, decel, 'r-')
peak_g = decel.max(); peak_g_t = t[decel.argmax()]
axes[1, 0].axvline(peak_g_t, color='k', ls=':', alpha=0.5)
axes[1, 0].set_xlabel('Time [s]'); axes[1, 0].set_ylabel('Deceleration [g]')
axes[1, 0].set_title(f'Peak deceleration: {peak_g:.1f} g @ t={peak_g_t:.0f}s')
axes[1, 0].grid(alpha=0.3)
axes[1, 1].plot(t, q_dot, 'm-')
peak_q = q_dot.max(); peak_q_t = t[q_dot.argmax()]
axes[1, 1].axvline(peak_q_t, color='k', ls=':', alpha=0.5)
axes[1, 1].set_xlabel('Time [s]'); axes[1, 1].set_ylabel(r'$\dot{q}_s$ [W/cm$^2$]')
axes[1, 1].set_title(f'Peak heating: {peak_q:.0f} W/cm² @ t={peak_q_t:.0f}s')
axes[1, 1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('mars_edl_profile.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから三つの重要な特徴が読み取れます。第一に、高度プロファイル（左上）は突入後 $60$–$80$ 秒で急速に密な大気に突入し、その後 $120\text{–}200$ 秒で「ほぼ水平飛行」のような緩やかな降下に転じます。これは揚力と遠心力の効果で、機体が大気上層を「滑空」している状態。Perseveranceの実フライトでも同じ挙動が報告されました。第二に、ピーク減速（左下）とピーク加熱（右下）はほぼ同じ時刻（突入後 $80\text{–}100$ 秒）に発生していること。これは前節で示した解析的予測通りで、$\rho V^2$ と $\sqrt{\rho}V^3$ の最大点が密度勾配のもとで一致するという定理を数値的に確認したことになります。第三に、終端速度（右上）はマッハ2.2相当の $\sim 450 \text{ m/s}$ を約 $200$ 秒後に下回り、ここがちょうど超音速パラシュート展開タイミングになります。
突入角sensitivity（コリドー幅）
次に、突入角 $\gamma_e$ を $-10^\circ$ から $-15^\circ$ までスキャンして、最終位置（着陸点ダウンレンジ）とピーク減速がどう変わるかを見ます。これが突入コリドー設計の基本です。
gammas = np.deg2rad(np.linspace(-15.0, -10.0, 11))
downranges = []
peak_decels = []
peak_heatings = []
for ge in gammas:
s0_g = [R_M + h_entry, 0.0, V_entry, ge]
sol_g = solve_ivp(dynamics, [0, 800], s0_g, events=hit_surface,
max_step=0.5, rtol=1e-9, atol=1e-9)
if sol_g.t_events[0].size == 0:
# 地表に到達せず大気圏外に跳ね返った（skip out）
downranges.append(np.nan); peak_decels.append(np.nan); peak_heatings.append(np.nan)
continue
# ダウンレンジ：θに半径をかけて距離換算
downranges.append(sol_g.y[1][-1] * R_M / 1000)
r_g, _, V_g, _ = sol_g.y
h_g = r_g - R_M
rho_g = rho_mars(h_g)
peak_decels.append((0.5 * rho_g * V_g**2 * CD * AREA / MASS / 9.81).max())
peak_heatings.append((K_MARS * np.sqrt(rho_g / RN) * V_g**3 * 1e-4).max())
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
gam_deg = np.rad2deg(gammas)
axes[0].plot(gam_deg, downranges, 'o-')
axes[0].set_xlabel('Entry flight path angle [deg]')
axes[0].set_ylabel('Downrange [km]')
axes[0].set_title('Downrange vs entry angle'); axes[0].grid(alpha=0.3)
axes[1].plot(gam_deg, peak_decels, 's-', color='red')
axes[1].axhline(15, color='k', ls='--', alpha=0.5, label='15g (crew limit)')
axes[1].set_xlabel('Entry flight path angle [deg]')
axes[1].set_ylabel('Peak deceleration [g]')
axes[1].set_title('Peak g-load'); axes[1].legend(); axes[1].grid(alpha=0.3)
axes[2].plot(gam_deg, peak_heatings, '^-', color='magenta')
axes[2].set_xlabel('Entry flight path angle [deg]')
axes[2].set_ylabel(r'Peak $\dot{q}_s$ [W/cm$^2$]')
axes[2].set_title('Peak heating'); axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('mars_edl_corridor.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、突入角が決定的に重要な三つのトレードオフが見えます。第一に、浅い突入角（$-10^\circ$側）ではダウンレンジが長く伸び、ピーク減速・ピーク加熱は緩やか。だが浅すぎると大気を貫いて「跳ね返って」（skip out）宇宙に戻ってしまうリスクがあり、ある臨界角以下では $\gamma > 0$ になって地表に到達しません。第二に、急な突入角（$-15^\circ$側）ではダウンレンジは短く、しかしピーク減速とピーク加熱が急増します。$-15^\circ$ では減速度が $\sim 12\text{–}13 g$ に達し、機器の構造限界に迫る。第三に、コリドー（許容突入角範囲）は概ね $-12^\circ$ ± $1\text{–}2^\circ$ で、実機の航法精度（$0.1^\circ$ 程度）を考えれば十分マージンがあります。Perseveranceの公称突入角は $\gamma_e = -11.99^\circ$、実測も $-12.04^\circ$ と1パーセントの精度に収まりました。
ここまでで「突入から地表 $\sim 11 \text{ km}$ あたりまで」の物理は把握できました。しかし機体は依然マッハ2程度。残る大気層では、パラシュートとロケットエンジンが主役を引き継ぎます。
Disk-Gap-Band 超音速パラシュート
なぜ「Disk-Gap-Band」形状か
火星のような薄い大気で超音速（マッハ 1.5–2.2）でパラシュートを展開するには、地球の一般的なリングスロット型では制御不能な振動（superSONIC parachute aerodynamic flapping）に陥ります。NASA Viking計画（1976）で開発された Disk-Gap-Band（DGB） 形状は、頂上の円板（Disk）、中央の隙間（Gap）、円筒状のバンド（Band）の3層構造で、超音速衝撃波背後の不安定流れを安定化します。
直径は Viking → Pathfinder → MER → MSL → Mars 2020 と段階的に拡大し、Mars 2020 では直径 21.5 m（プロジェクトADVANCEで超音速展開試験を新規実施）。展開時の制約は2つ：マッハ数 $M < 2.2$ と動圧 $\bar{q} < 850 \text{ Pa}$。これを越えるとキャノピーが裂けるか、安定展開できません。
展開トリガロジック
実機では (a) ナビゲーション解からマッハ数を推定、(b) 加速度計から動圧を推定、(c) 両条件を満たした最初の時刻に展開、という論理。Mars 2020はさらにRange Triggerを採用し、目標着陸点までの予測残距離に基づいて展開時刻を動的に調整します。これにより着陸誤差楕円が従来の $25 \text{ km}$ から $7 \text{ km}$ まで縮小しました。
パラシュート減速とバックシェル分離
展開直後の動圧 $\bar{q} \approx 800 \text{ Pa}$、抗力係数 $C_{D,p} \approx 0.6$、直径 21.5 m として抗力は
$$ D_p = \frac{1}{2} \times 0.014 \times 450^2 \times 0.6 \times \pi \times 10.75^2 \approx 1.5 \times 10^5 \text{ N} $$
機体質量 $\sim 3000 \text{ kg}$ で減速度 $\sim 5g$。これでマッハ 2 → マッハ 0.5（音速以下）まで急減速します。亜音速に入ったところでヒートシールドを分離（Heat Shield Separation, 突入から約260秒）し、続いてレーダーで地表までの絶対高度を取得、そしてバックシェル分離（パラシュートとカプセル後部を切り離す）。
ここで Powered Descent Vehicle（PDV、ローバ＋ジェットパック）が単独で落下を始めます。残る速度は秒速 $\sim 80\text{–}100 \text{ m/s}$。これを Powered Descent でゼロにします。
Powered Descent と スカイクレーン
MARS Landing Engine（MLE）
Mars 2020 の Descent Stage（ジェットパック部）は 8基の MLE スラスタ（各 3,000 N、燃料：ヒドラジン）を搭載。質量制御推進剤搭載量は約 $390 \text{ kg}$、総推力 $\sim 24 \text{ kN}$ で重量 $\sim 1900 \text{ kg}$ の PDV を支えます（火星重力 $g_M = 3.71 \text{ m/s}^2$ で重量は $\sim 7000 \text{ N}$、ホバリング推力比 $\sim 3.4$）。
推進剤消費はツィオルコフスキー則で見積もれます。目標 $\Delta V = 80 \text{ m/s}$、比推力 $I_{sp} = 220 \text{ s}$（ヒドラジン典型値）、初期質量 $m_0 = 1900 \text{ kg}$ なら
$$ m_p = m_0 \left(1 – \exp\left(-\frac{\Delta V}{I_{sp} g_0}\right)\right) = 1900 \left(1 – e^{-80/(220 \times 9.81)}\right) \approx 70 \text{ kg} $$
実機の搭載量 $390 \text{ kg}$ に対して使用量 $\sim 70 \text{ kg}$。残りは姿勢制御、Divert Maneuver（緊急横移動）、そしてホバリング時の長時間燃焼に充てられます。
スカイクレーン：1トン級ローバの「吊り下げ着陸」
ここからが Mars Science Laboratory（2012, Curiosity）が編み出した革新です。1トン級のローバを直接エアバッグで着地させるのは不可能（衝撃が大きすぎる、ローバの突起構造が壊れる）、ロケット脚で着地させるのもダスト巻き上げで光学・センサが汚染される。そこで考案されたのがスカイクレーン：
- PDV が高度 $\sim 20 \text{ m}$ で 「Constant Velocity Descent」（一定速度降下、$0.75 \text{ m/s}$）に移行
- ローバを4本のケーブル + 1本のアンビリカルで $7.5 \text{ m}$ 吊り下げ展開
- ローバが地表に「触れた」瞬間（推力急減で検知）にケーブルを瞬時に火薬で切断
- ジェットパックは自律的に横方向飛行 → 墜落してローバから離れる
これにより、(a) ローバの直下がロケット排気で汚染されず、(b) 着陸衝撃が車輪のサスペンションで吸収、(c) ローバが展開完了状態で着陸できる、という三重の利点が得られました。Curiosity と Perseverance はこの方式で、半世紀の火星探査史で最も精密かつ大型の着陸を達成しています。
Terrain Relative Navigation（TRN）
Mars 2020で初投入された新技術が TRN。パラシュート展開後の数十秒間、機体下面のカメラが地表を撮影し、事前にロードされたJezero クレータの軌道画像と特徴点マッチング、自機の絶対位置を推定します。これにより、たとえ航法誤差で危険な岩場や崖の真上に降下していても、Powered Descent中に横移動（divert）して安全地帯に着陸できます。
実際の Perseveranceは TRN によって、危険判定された一次着陸点から約 200 m横にdivertして着陸しました。これがなければ岩場に落ちて転倒していた可能性があり、TRNは Mars 2020 の最大の保険装置の一つです。
ここまでで「成功する火星EDLの構造」が見えてきました。逆に「失敗するとどうなるか」を、Schiaparelliの例で見ておきましょう。
Schiaparelli 2016：論理エラーが招いた墜落
概要
ESA-Roscosmos の ExoMars 2016 ミッションに搭載された Schiaparelli は、ESA独自の火星着陸技術実証機。2016年10月19日、突入は順調、パラシュート展開も成功。しかしPowered Descent 開始の判断ロジックでエラーが連鎖し、結果として高度 $3.7 \text{ km}$ で逆噴射を停止、機体は秒速 $300 \text{ m/s}$ で地表に衝突して大破しました。
失敗の連鎖
ESA の事故調査報告によると、原因は以下のような連鎖：
- パラシュート展開直後の急激な振動で IMU（慣性計測装置）が短時間サチュレーション（出力上限に張り付き）
- ナビゲーションフィルタは IMU データを積分して姿勢推定 → 飽和データを誤って積分し続け、「機体が地表下にいる」という非物理的な高度推定を出力
- システムは「高度ゼロ以下 → すでに着地している」と誤判断し、レトロロケット点火直後にエンジンカットオフ
- 実際にはまだ高度 $3.7 \text{ km}$ で秒速 $80 \text{ m/s}$、自由落下で加速し $\sim 300 \text{ m/s}$ で激突
教訓
この失敗が示すのは、火星EDLが「個々のサブシステムが正しく動くだけでは不十分」で、(a) センサ飽和時のフェイルセーフ、(b) ナビゲーションフィルタの異常値検出（sanity check）、(c) 複数センサのクロスチェック（レーダー高度計とIMU推定の整合性確認）といったシステム統合レベルの冗長性が不可欠だということ。Mars 2020 では IMU、レーダー高度計、ドップラーレーダー（速度計測）、TRNカメラの四重冗長で同様の単点故障を排除しています。
ここまでで物理と運用論理の両面が出揃いました。最後に、火星大気モデルそのものを少し深掘りして、運用上の不確実性がどこにあるかを見ておきます。
MarsGRAMライクな大気モデルと不確実性
単純指数モデルの限界
実際の火星大気は、高度プロファイルが単純な指数則からずれます。理由は (1) 温度が高度ごとに変動（成層圏で逆転層がある）、(2) ダスト含有量が密度に影響、(3) 季節（特に南半球の冬の二酸化炭素霜の蒸発）で全体気圧が $\pm 30\%$ 変動。
ここでは、温度プロファイル付きの指数則（密度が高度の関数として温度補正されたもの）と、ダスト含有量バリエーションを加えた簡易モデルを比較してみます。
import numpy as np
import matplotlib.pyplot as plt
# 高度配列
h = np.linspace(0, 130e3, 500)
# 1. 等温指数モデル
rho_iso = RHO0 * np.exp(-h / H_SCALE)
# 2. 温度プロファイル（簡易版：地表210K、上層140K に線形補間）
T0 = 210.0
T_top = 140.0
h_thermo = 100e3
T = T0 - (T0 - T_top) * np.clip(h / h_thermo, 0, 1)
# 圧力は静水圧平衡から積分（理想気体）
R_gas = 188.92 # CO2 比気体定数 [J/(kg·K)]
g_M = 3.71
# dp/dh = -rho g = -p/(RT) g → ln(p/p0) = -∫ g/(RT) dh
integ = np.concatenate([[0], np.cumsum(0.5 * (g_M/(R_gas*T[1:]) + g_M/(R_gas*T[:-1])) * np.diff(h))])
p = (RHO0 * R_gas * T0) * np.exp(-integ) # 圧力 [Pa]
rho_temp = p / (R_gas * T)
# 3. ダストストーム想定（密度が下層で +30%、上層で +60% 増加）
dust_factor = 1.0 + 0.3 * np.exp(-h / 30e3) + 0.3 * np.exp(-(h-50e3)**2 / (20e3)**2)
rho_dust = rho_temp * dust_factor
# 4. ダスト稀薄期（密度が -20%）
rho_clear = rho_temp * 0.8
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].semilogx(rho_iso, h/1000, 'k--', label='Isothermal exp')
axes[0].semilogx(rho_temp, h/1000, 'b-', label='Temperature-corrected')
axes[0].semilogx(rho_dust, h/1000, 'r-', label='Dust storm (high)')
axes[0].semilogx(rho_clear, h/1000, 'g-', label='Clear season (low)')
axes[0].set_xlabel(r'Density $\rho$ [kg/m$^3$]')
axes[0].set_ylabel('Altitude [km]')
axes[0].set_title('Mars atmospheric density profiles')
axes[0].legend(); axes[0].grid(alpha=0.3, which='both')
axes[1].plot(T, h/1000, 'b-')
axes[1].set_xlabel('Temperature [K]')
axes[1].set_ylabel('Altitude [km]')
axes[1].set_title('Temperature profile (simplified MarsGRAM)')
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.savefig('mars_atmosphere_models.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから二つの重要な読み取りができます。第一に、単純指数モデル（黒破線）と温度補正モデル（青実線）は地表近くではほぼ一致しますが、高度 $50\text{–}100 \text{ km}$ のピーク加熱領域で最大2倍の差が生じます。EDL のクリティカルな高度がまさにここなので、設計マージンを取るときには温度補正が不可欠。第二に、ダストストーム期（赤実線）と稀薄期（緑実線）の密度差はほぼ2倍で、これがそのままピーク動圧・ピーク加熱・パラシュート展開高度に伝搬します。Mars 2020 のEDLシステムは、Monte Carloシミュレーションで $\pm 30\%$ の密度変動を仮定して、全パラメータが許容範囲内に収まることを事前に検証しました。
Monte Carlo によるロバスト性検証の概念
実運用では、密度、突入角、突入速度、CD、CL の各パラメータをガウシアンでランダムサンプリングし、数千ケースの軌道シミュレーションを実行、最終着陸誤差や最大荷重の確率分布を作ります。簡易版を回してみましょう。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
n_runs = 200
results = {'final_V': [], 'peak_g': [], 'peak_q': [], 'downrange': []}
for _ in range(n_runs):
# パラメータ摂動
ge = np.deg2rad(np.random.normal(-12.0, 0.3))
Ve = np.random.normal(5400.0, 30.0)
cd_pert = CD * np.random.normal(1.0, 0.05)
rho_scale = np.random.normal(1.0, 0.15) # ±15% 密度変動
def dyn_pert(t, s):
r, theta, V, gamma = s
h_l = r - R_M
rho_l = RHO0 * np.exp(-h_l / H_SCALE) * rho_scale
g_l = MU_M / r**2
q_bar_l = 0.5 * rho_l * V**2
D_l = q_bar_l * cd_pert * AREA
L_l = q_bar_l * CL * AREA
return [V*np.sin(gamma), V*np.cos(gamma)/r,
-D_l/MASS - g_l*np.sin(gamma),
L_l/MASS/V - (g_l - V**2/r)*np.cos(gamma)/V]
s0p = [R_M + h_entry, 0.0, Ve, ge]
solp = solve_ivp(dyn_pert, [0, 800], s0p, events=hit_surface,
max_step=1.0, rtol=1e-7, atol=1e-7)
if solp.t_events[0].size == 0:
continue
r_p, _, V_p, _ = solp.y
h_p = r_p - R_M
rho_p = RHO0 * np.exp(-h_p / H_SCALE) * rho_scale
results['final_V'].append(V_p[-1])
results['peak_g'].append((0.5*rho_p*V_p**2*cd_pert*AREA/MASS/9.81).max())
results['peak_q'].append((K_MARS*np.sqrt(rho_p/RN)*V_p**3*1e-4).max())
results['downrange'].append(solp.y[1][-1] * R_M / 1000)
fig, axes = plt.subplots(2, 2, figsize=(11, 7))
axes[0,0].hist(results['final_V'], bins=20, color='green', alpha=0.7)
axes[0,0].set_xlabel('Final velocity [m/s]'); axes[0,0].set_title('Pre-parachute V')
axes[0,1].hist(results['peak_g'], bins=20, color='red', alpha=0.7)
axes[0,1].set_xlabel('Peak deceleration [g]'); axes[0,1].set_title('Peak g')
axes[1,0].hist(results['peak_q'], bins=20, color='magenta', alpha=0.7)
axes[1,0].set_xlabel(r'Peak heating [W/cm$^2$]'); axes[1,0].set_title('Peak heating')
axes[1,1].hist(results['downrange'], bins=20, color='blue', alpha=0.7)
axes[1,1].set_xlabel('Downrange [km]'); axes[1,1].set_title('Downrange distribution')
for ax in axes.flatten(): ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('mars_edl_montecarlo.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"Final V: mean={np.mean(results['final_V']):.0f}, std={np.std(results['final_V']):.0f} m/s")
print(f"Peak g : mean={np.mean(results['peak_g']):.1f}, std={np.std(results['peak_g']):.1f}")
print(f"Downrange: mean={np.mean(results['downrange']):.0f}, std={np.std(results['downrange']):.0f} km")
この Monte Carlo の結果が示すのは、火星 EDL におけるパラメータ不確実性が最も強く伝搬するのはダウンレンジで、その標準偏差は数十 km に達します。これが「着陸楕円（landing ellipse）」の正体で、$3\sigma$ で $\sim 25 \text{ km}$（MSL）から $\sim 7 \text{ km}$（Mars 2020、Range Triggerあり）まで縮小されてきた経緯がよくわかります。一方、ピーク減速とピーク加熱は分布の中央値が物理的な許容範囲（$15g$ 以下、$200 \text{ W/cm}^2$ 以下）に十分余裕があり、システム設計の安全マージンが効いていることを定量的に示しています。最終速度（パラシュート展開時）はおおよそ $400\text{–}500 \text{ m/s}$ にうまく収まっており、Disk-Gap-Band の展開条件（$M < 2.2 \Leftrightarrow V < 450\text{--}500 \text{ m/s}$）と整合的です。
エアバッグ・レトロロケット・スカイクレーンの使い分け
歴代の火星着陸方式を簡単に整理します。
| 方式 | ミッション | 着陸質量 | 利点 | 欠点 |
|---|---|---|---|---|
| エアバッグ | Pathfinder (1997), MER (2004) | $\sim 350 \text{ kg}$ | シンプル、安価、地形耐性 | 質量制限が厳しい、跳ね返り着陸点不確実 |
| レトロロケット | Viking (1976), Phoenix (2008), InSight (2018) | $\sim 350\text{–}600 \text{ kg}$ | 中質量対応、精密着陸 | ダスト巻き上げ、機体直下汚染 |
| スカイクレーン | MSL (2012), Mars 2020 (2021) | $\sim 900\text{–}1000 \text{ kg}$ | 大型ローバ可、汚染なし、展開状態着陸 | 複雑、有人スケールに拡張困難 |
将来の有人ミッション（数十トン級）では、スカイクレーンの単純な拡張は構造的に成立せず、Supersonic Retro-Propulsion（超音速のうちにロケットで逆噴射する）や、ハイパーソニック・インフレータブル空力減速装置（HIAD）といった次世代技術の研究が進んでいます。SpaceX Starship のような大型機の火星着陸計画は、これらの未実証技術が前提です。
まとめ
本記事では、火星EDLの「7分間の恐怖」を物理から運用までの全層で解説しました。
- 火星大気は地球の1%だが減速予算は致命的に足りない — 突入運動エネルギーが大きすぎ、パラシュートだけでは止まらない。Powered Descent が必須。
- 加熱は地球並み — Sutton-Graves 式 $\dot{q}_s = K\sqrt{\rho/R_N}\,V^3$ より、$V^3$ 依存性が支配的。$\text{CO}_2$ 大気のため係数 $K$ は地球より9%大きい。
- 突入軌道は球座標の4変数 ODE — $\dot{r}, \dot{\theta}, \dot{V}, \dot{\gamma}$ の連立系を
scipy.integrate.solve_ivpで積分すれば、ピーク加熱・減速・着陸時刻を全て予測できる。 - ピーク加熱とピーク減速はほぼ同時刻 — 突入後 $80\text{–}100$ 秒、高度 $30\text{–}40 \text{ km}$ で発生。
- 超音速パラシュート Disk-Gap-Band — マッハ $< 2.2$、動圧 $< 850 \text{ Pa}$ で展開。Mars 2020 の Range Trigger で着陸誤差を $7 \text{ km}$ まで縮小。
- スカイクレーン — 1トン級ローバの精密着陸を可能にし、汚染と着陸衝撃を同時解決する独創的方式。
- Terrain Relative Navigation — 地表画像と軌道画像のマッチングで絶対位置を取得、危険地形を自律 divert。
- Schiaparelli の教訓 — IMU 飽和から「高度マイナス」の非物理出力が発生、エンジン誤切断。システム統合レベルの冗長性が不可欠。
- Monte Carlo によるロバスト性検証 — 密度 $\pm 15\%$、突入角 $\pm 0.3^\circ$ 等の摂動下で全パラメータが許容範囲内に収まることを事前検証する。これが実運用の保証手段。
火星 EDL は、突入軌道力学、超音速空力、熱流体、構造設計、ナビゲーション、ロケット工学、システム工学が全て交差する宇宙工学の総合競技です。次のステップとして、突入機の熱・構造設計の詳細、惑星間遷移軌道の入り口、そして再突入予測の数学的フレームワークを扱う以下の記事を参考にしてください。