静止衛星は地球の赤道上空約36,000 kmで、地球の自転と同期して常に同じ経度にとどまっているはずです。ところが実際には、何もしなければ静止衛星は数年のうちに軌道が大きく傾き、最終的には使い物にならなくなります。この「見えない力」の正体は何でしょうか？
答えは月と太陽の重力です。地球のまわりを周回する人工衛星は、地球の重力だけでなく、月や太陽からの重力も受けています。二体問題では地球と衛星の2つの天体しか考えませんが、現実には第三の天体（月・太陽・他の惑星）からの重力が衛星の軌道をじわじわと変化させます。これが第三体摂動です。
第三体摂動を理解すると、以下のような実際の宇宙ミッションの課題に答えられるようになります。
- 静止衛星の軌道維持: GEO衛星が月・太陽の摂動で傾斜角が年間約0.75°ずつ増加する問題と、それに必要なステーションキーピングの$\Delta V$
- 月遷移軌道の設計: 月の重力を積極的に利用した弱安定境界（WSB）軌道の概念
- GPS衛星の長期予測: MEO軌道（高度約20,200 km）の衛星が月・太陽摂動で受ける離心率変動の予測
本記事の内容
- 二体問題の限界と摂動論の必要性
- 第三体からの摂動力の導出
- 摂動ポテンシャルとルジャンドル展開
- 平均化法による長期変動の解析
- コズライの共鳴（離心率・傾斜角の長期振動）
- Pythonで月の重力摂動による軌道要素変化をシミュレーション
- GEO衛星への影響と軌道維持
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
二体問題の限界 — なぜ摂動を考えるのか
完全な二体問題の世界
二体問題では、衛星は地球の重力だけを受けて運動します。このとき衛星の軌道はケプラー楕円で完全に記述され、軌道要素（半長軸 $a$、離心率 $e$、傾斜角 $i$、昇交点赤経 $\Omega$、近地点引数 $\omega$、真近点角 $\nu$）のうち、最初の5つは永久に変化しません。真近点角だけが時間とともに変化し、衛星は同じ楕円上を何度も周回します。
しかし現実の衛星は「完璧な楕円」を描きません。地球は完全な球ではないし（J2摂動）、月や太陽の重力も作用し、大気抵抗も存在します。特に高軌道の衛星では、月・太陽の重力が最も支配的な摂動源となります。
摂動の大きさの目安
摂動がどの程度の影響を持つかを見積もるために、高度36,000 km（GEO軌道）にある衛星に対する各摂動加速度の大きさを比較してみましょう。
| 摂動源 | 加速度の大きさ (m/s²) | 地球重力に対する比率 |
|---|---|---|
| 地球重力（主項） | $\sim 0.22$ | 1 |
| J2（地球扁平率） | $\sim 3 \times 10^{-5}$ | $\sim 10^{-4}$ |
| 月の重力 | $\sim 5 \times 10^{-6}$ | $\sim 2 \times 10^{-5}$ |
| 太陽の重力 | $\sim 2 \times 10^{-6}$ | $\sim 10^{-5}$ |
| 太陽輻射圧 | $\sim 5 \times 10^{-8}$ | $\sim 2 \times 10^{-7}$ |
表を見ると、月と太陽の摂動加速度は地球重力の10万分の1程度と非常に小さいです。しかし、この小さな力が何年にもわたって蓄積されると、軌道要素は大きく変化します。たとえば、GEO衛星の傾斜角は年間約0.75°ずつ増加し、15年程度で最大約15°に達します。衛星通信のアンテナ追尾が困難になるレベルです。
このように、短期的には無視できる程度の力でも、長期的には軌道に大きな変化をもたらすのが第三体摂動の特徴です。では、この摂動力を数学的にどう扱うのかを見ていきましょう。
第三体からの摂動力の導出
問題設定
地球を中心天体（質量 $M_E$）、衛星を軌道天体（質量 $m$、$m \ll M_E$）、月（または太陽）を第三体（質量 $M_3$）とします。地球を原点とする座標系において、衛星の位置ベクトルを $\bm{r}$、第三体の位置ベクトルを $\bm{r}_3$ と表記します。
衛星に作用する力は、地球からの重力と第三体からの重力の合計です。
$$ m\ddot{\bm{r}} = -\frac{GM_E m}{r^3}\bm{r} + GM_3 m \left(\frac{\bm{r}_3 – \bm{r}}{|\bm{r}_3 – \bm{r}|^3} – \frac{\bm{r}_3}{r_3^3}\right) $$
ここで重要なのは、右辺第2項の構造です。衛星にかかる第三体の重力は2つの部分から成っています。
第1の部分 $GM_3(\bm{r}_3 – \bm{r})/|\bm{r}_3 – \bm{r}|^3$ は、第三体が衛星を直接引く力です。第三体から衛星の方を向いた重力そのものです。
第2の部分 $-GM_3\bm{r}_3/r_3^3$ は、第三体が地球を引く力です。地球を原点にした座標系を使っているため、地球自体が第三体に引かれて加速する効果を差し引く必要があります。この項は間接項と呼ばれます。
摂動加速度
$m$ で割って整理すると、衛星の運動方程式は次のようになります。
$$ \ddot{\bm{r}} = -\frac{\mu}{r^3}\bm{r} + \bm{a}_{\text{pert}} $$
ここで $\mu = GM_E$ は地球の重力パラメータ、$\bm{a}_{\text{pert}}$ は摂動加速度です。
$$ \bm{a}_{\text{pert}} = GM_3\left(\frac{\bm{r}_3 – \bm{r}}{|\bm{r}_3 – \bm{r}|^3} – \frac{\bm{r}_3}{r_3^3}\right) $$
この式の意味を直感的に理解しましょう。もし衛星が地球の中心にいたら（$\bm{r} = \bm{0}$）、2つの項は完全にキャンセルして $\bm{a}_{\text{pert}} = \bm{0}$ になります。つまり、摂動加速度は「衛星の位置での第三体の重力」と「地球の中心での第三体の重力」の差です。これを潮汐力（tidal force）と呼びます。地球と衛星の位置が離れているほど、この差が大きくなります。
摂動力の基本構造がわかったところで、次にこの力をポテンシャルとして表現し、ルジャンドル多項式を使って系統的に展開する方法を見ていきます。
摂動ポテンシャルとルジャンドル展開
摂動ポテンシャルの定義
摂動加速度はポテンシャルの勾配として表すことができます。摂動ポテンシャル $R$ を次のように定義します。
$$ R = GM_3\left(\frac{1}{|\bm{r}_3 – \bm{r}|} – \frac{\bm{r} \cdot \bm{r}_3}{r_3^3}\right) $$
すると $\bm{a}_{\text{pert}} = \nabla_{\bm{r}} R$ が成り立ちます。第1項は第三体による直接的な重力ポテンシャル、第2項は間接項のポテンシャルです。
ルジャンドル展開
人工衛星の軌道半径 $r$ は第三体までの距離 $r_3$ に比べて小さいです。月までの距離は約384,400 kmですが、GEO軌道半径は約42,164 kmなので、$r/r_3 \approx 0.11$ です。LEO（高度400 km）なら $r/r_3 \approx 0.018$ とさらに小さくなります。
この $r/r_3 \ll 1$ という条件を利用して、$1/|\bm{r}_3 – \bm{r}|$ をルジャンドル多項式で展開します。衛星と第三体の位置ベクトルのなす角を $\psi$ とすると
$$ \frac{1}{|\bm{r}_3 – \bm{r}|} = \frac{1}{r_3}\sum_{n=0}^{\infty}\left(\frac{r}{r_3}\right)^n P_n(\cos\psi) $$
ここで $P_n$ は $n$ 次のルジャンドル多項式です。最初のいくつかの項を書き出すと
$$ P_0(\cos\psi) = 1, \quad P_1(\cos\psi) = \cos\psi, \quad P_2(\cos\psi) = \frac{1}{2}(3\cos^2\psi – 1) $$
この展開を摂動ポテンシャルに代入してみましょう。$n=0$ の項は $1/r_3$ で、$\bm{r}$ に依存しないため摂動力を生みません（勾配がゼロ）。$n=1$ の項は $r\cos\psi/r_3^2 = \bm{r}\cdot\bm{r}_3/r_3^3$ で、間接項と正確にキャンセルします。
したがって、摂動ポテンシャルの実質的な寄与は $n \geq 2$ の項から始まります。
$$ R = \frac{GM_3}{r_3}\sum_{n=2}^{\infty}\left(\frac{r}{r_3}\right)^n P_n(\cos\psi) $$
主要項（$n=2$）の近似
$r/r_3 \ll 1$ なので、最低次の $n=2$ の項が最も大きな寄与を持ちます。$n=2$ のみを残すと
$$ R \approx \frac{GM_3}{r_3}\left(\frac{r}{r_3}\right)^2 P_2(\cos\psi) = \frac{GM_3 r^2}{2r_3^3}(3\cos^2\psi – 1) $$
この式が第三体摂動ポテンシャルの主要項です。物理的に見ると、この摂動ポテンシャルの大きさは $r^2$ に比例します。つまり、軌道半径が大きい衛星ほど第三体摂動の影響を強く受けます。これはGEO衛星がLEO衛星よりも月・太陽摂動を強く受ける理由を直観的に説明しています。
摂動ポテンシャルが得られたので、次はこのポテンシャルが軌道要素にどのような変化を引き起こすかを、ラグランジュの惑星方程式を使って調べます。
ラグランジュの惑星方程式と摂動の作用
ラグランジュの惑星方程式
摂動ポテンシャル $R$ が与えられたとき、軌道要素の時間変化率はラグランジュの惑星方程式で記述されます。ここでは主要な3つの軌道要素（離心率 $e$、傾斜角 $i$、昇交点赤経 $\Omega$）に注目します。
$$ \frac{de}{dt} = \frac{\sqrt{1-e^2}}{na^2 e}\left[\frac{\partial R}{\partial \omega} – \frac{1-e^2}{e}\cdot\frac{1}{\sqrt{1-e^2}}\frac{\partial R}{\partial M}\right] $$
$$ \frac{di}{dt} = \frac{1}{na^2\sqrt{1-e^2}\sin i}\left[\cos i \frac{\partial R}{\partial \omega} – \frac{\partial R}{\partial \Omega}\right] $$
$$ \frac{d\Omega}{dt} = \frac{1}{na^2\sqrt{1-e^2}\sin i}\frac{\partial R}{\partial i} $$
ここで $n = \sqrt{\mu/a^3}$ は平均運動（1周期あたりの角速度）、$M$ は平均近点角、$\omega$ は近地点引数です。
これらの方程式はかなり複雑ですが、ポイントは「摂動ポテンシャル $R$ の偏微分を通じて、各軌道要素がどう変化するかが決まる」ということです。$R$ の具体的な形さえわかれば、計算は機械的に進められます。
短周期変動と長周期変動
第三体の摂動ポテンシャルを軌道要素で表すと、衛星の真近点角 $\nu$（またはその代わりの平均近点角 $M$）と第三体の位置の両方に依存します。この依存性は2種類の変動を引き起こします。
短周期変動は衛星の公転周期（LEOなら約90分、GEOなら24時間）のスケールで振動する変動です。1周ごとに平均すると消えるため、長期的な軌道進化には直接影響しません。
長周期変動は月の公転周期（約27.3日）や太陽の見かけの公転周期（約365.25日）のスケール、あるいはそれよりさらに長い時間スケールで変化する成分です。これが軌道要素の永年的な変化（セキュラー変動）を引き起こします。
長期的な軌道進化を理解するためには、短周期変動を平均化して取り除く「平均化法」が有効です。次のセクションでこの手法を詳しく見ていきましょう。
平均化法による長期変動の解析
平均化法の考え方
平均化法（averaging method）は、摂動ポテンシャルを「速い変数」（衛星の平均近点角 $M$ や第三体の経度）について1周分平均することで、短周期変動を取り除く手法です。イメージとしては、振動するノイズを取り除いてトレンドだけを抽出するローパスフィルタのようなものです。
単平均化（single averaging） では、衛星の平均近点角 $M$ について1周分の平均をとります。
$$ \bar{R} = \frac{1}{2\pi}\int_0^{2\pi} R \, dM $$
二重平均化（double averaging） では、さらに第三体の軌道運動についても平均をとります。
$$ \bar{\bar{R}} = \frac{1}{2\pi}\int_0^{2\pi}\bar{R} \, dM_3 $$
二重平均化を行うと、摂動ポテンシャルは $a$、$e$、$i$、$\omega$ のみの関数になり、$M$ と $\Omega$ に依存しなくなります（太陽摂動の場合）。
二重平均化ポテンシャル
$n=2$ の主要項に対して二重平均化を行うと、次の結果が得られます。
$$ \bar{\bar{R}} = \frac{GM_3 a^2}{4r_3^3}\left[(1 + \tfrac{3}{2}e^2)(3\cos^2 i – 1) + \tfrac{15}{2}e^2\sin^2 i \cos 2\omega\right] $$
ここで $r_3$ は第三体の平均距離です。この式から2つの重要な特徴が読み取れます。
第1に、右辺第1項は $\omega$ に依存しません。これはセキュラー項と呼ばれ、$e$ と $i$ の長期的な変化率を決定します。
第2に、右辺第2項は $\cos 2\omega$ を含みます。近地点引数 $\omega$ が時間とともに変化すると、この項は周期的に振動します。この項がコズライの共鳴（Kozai-Lidov mechanism）の源です。
セキュラー変動
二重平均化ポテンシャルをラグランジュの惑星方程式に代入すると、$a$ の永年変動はゼロであることが示されます。これは重要な結果で、第三体摂動は軌道のエネルギー（したがって半長軸）を長期的には変化させないことを意味します。
一方、$\Omega$（昇交点赤経）には永年的な変化が生じます。太陽の摂動による $\Omega$ の永年変化率は
$$ \dot{\Omega}_{\odot} \approx -\frac{3n_{\odot}^2}{4n}\cos i $$
ここで $n_{\odot}$ は太陽の平均運動（$\approx 2\pi / 365.25$ rad/day）です。月の摂動も同様の形を持ちますが、月の質量と距離に応じた係数が異なります。
平均化法で長期的な傾向がわかりましたが、実はこの枠組みの中でさらに劇的な現象が隠れています。それが次に説明するコズライの共鳴です。
コズライの共鳴 — 離心率と傾斜角の結合振動
コズライの発見
1962年、日本の天文学者・古在由秀は、小惑星の軌道に対する木星の摂動を研究する中で驚くべき現象を発見しました。高い傾斜角を持つ天体の軌道では、離心率 $e$ と傾斜角 $i$ が数十年〜数百年の時間スケールで大きく振動するのです。同時期にリドフも同様の結果を独立に得たため、この現象はコズライ-リドフ機構（Kozai-Lidov mechanism）と呼ばれます。
日常のアナロジーでたとえるなら、ブランコの左右の揺れ（傾斜角の変動）と前後の揺れ（離心率の変動）がエネルギーをやりとりしながら、一方が大きくなると他方が小さくなるような関係です。
保存量
二重平均化ポテンシャルのもとでは、半長軸 $a$ に加えて、以下の量が保存されます。
$$ H = \sqrt{1-e^2}\cos i $$
この $H$ はしばしばコズライ積分と呼ばれます。物理的には、軌道角運動量の $z$ 成分（赤道面に垂直な成分）に比例する量です。
$H$ が保存されることの意味を考えましょう。$H = \sqrt{1-e^2}\cos i$ が一定なので、離心率 $e$ が増加すると $\sqrt{1-e^2}$ が減少し、$\cos i$ が増加する必要があります。つまり離心率が大きくなるとき、傾斜角は小さくなります。逆に離心率が小さくなるとき、傾斜角は大きくなります。これがコズライ共鳴の本質です。
臨界傾斜角
コズライ共鳴が起こるかどうかは初期傾斜角に依存します。臨界条件は次の式で与えられます。
初期離心率 $e_0 \approx 0$ の場合、コズライ振動が起こる条件は
$$ \cos^2 i_0 < \frac{3}{5} \quad \Longleftrightarrow \quad 39.2° < i_0 < 140.8° $$
この範囲の傾斜角を持つ軌道では、離心率が0から大きな値まで周期的に振動する可能性があります。極端な場合、ほぼ円軌道だった天体が非常に細長い楕円軌道になり、中心天体に衝突するほど近づくこともあります。
コズライ振動の周期
コズライ振動の特性的な時間スケールは次の式で見積もれます。
$$ T_{\text{Kozai}} \sim \frac{M_E}{M_3}\frac{P_3^2}{P}\left(1-e_3^2\right)^{3/2} $$
ここで $P$ は衛星の公転周期、$P_3$ は第三体の公転周期、$e_3$ は第三体の離心率です。月による摂動の場合、GEO衛星のコズライ周期は数十年のオーダーになります。
コズライ共鳴は天然の天体（不規則衛星、系外惑星のホットジュピター形成）でも重要な役割を果たしますが、人工衛星の場合はGEO衛星や中高度軌道衛星の軌道進化に大きな影響を与えます。ここからは、Pythonを使ってこれらの理論を数値的に確認しましょう。
Pythonで月の摂動による軌道要素変化をシミュレーション
シミュレーションの設計
ここでは、GEO軌道に近い衛星に対する月の重力摂動をシミュレーションします。衛星の運動方程式を地球の重力と月の重力の和として記述し、数値積分で軌道を追跡します。月は地球のまわりを円軌道で公転すると簡略化します。
まず、基本定数と初期条件を設定するコードを書きます。
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# --- 定数 ---
mu_earth = 3.986004418e14 # 地球の重力パラメータ [m^3/s^2]
mu_moon = 4.9028e12 # 月の重力パラメータ [m^3/s^2]
r_moon = 384400e3 # 月の平均軌道半径 [m]
n_moon = 2 * np.pi / (27.3217 * 86400) # 月の平均運動 [rad/s]
# --- 初期軌道（GEO付近） ---
a0 = 42164e3 # 半長軸 [m]（GEO）
e0 = 0.001 # 離心率（ほぼ円）
i0 = np.radians(23.0) # 傾斜角 [rad]（初期23°）
RAAN0 = np.radians(0.0) # 昇交点赤経 [rad]
omega0 = np.radians(0.0) # 近地点引数 [rad]
nu0 = 0.0 # 真近点角 [rad]
次に、軌道要素から状態ベクトル（位置・速度）への変換関数を準備します。
def kep_to_cart(a, e, i, RAAN, omega, nu, mu):
"""ケプラー要素 → 慣性系の位置・速度ベクトル"""
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 / p) * np.array([-np.sin(nu), e + np.cos(nu), 0.0])
# 回転行列（PQW → ECI）
cos_O, sin_O = np.cos(RAAN), np.sin(RAAN)
cos_i, sin_i = np.cos(i), np.sin(i)
cos_w, sin_w = np.cos(omega), np.sin(omega)
R = np.array([
[cos_O*cos_w - sin_O*sin_w*cos_i,
-cos_O*sin_w - sin_O*cos_w*cos_i,
sin_O*sin_i],
[sin_O*cos_w + cos_O*sin_w*cos_i,
-sin_O*sin_w + cos_O*cos_w*cos_i,
-cos_O*sin_i],
[sin_w*sin_i, cos_w*sin_i, cos_i]
])
r_eci = R @ r_pqw
v_eci = R @ v_pqw
return r_eci, v_eci
def cart_to_kep(r_vec, v_vec, mu):
"""慣性系の位置・速度 → ケプラー要素"""
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 - r_vec / r
e = np.linalg.norm(e_vec)
# 半長軸
energy = v**2 / 2 - mu / r
a = -mu / (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
そして、運動方程式を定義して数値積分を実行します。
def moon_position(t):
"""月の位置（赤道面上の円軌道と近似）"""
# 月の軌道傾斜角（赤道面に対して約28.5°の最大値）
i_moon = np.radians(23.44) # 黄道傾斜角で近似
theta = n_moon * t
x = r_moon * np.cos(theta)
y = r_moon * np.sin(theta) * np.cos(i_moon)
z = r_moon * np.sin(theta) * np.sin(i_moon)
return np.array([x, y, z])
def equations_of_motion(t, state):
"""地球重力 + 月の第三体摂動"""
r_vec = state[0:3]
v_vec = state[3:6]
r = np.linalg.norm(r_vec)
# 地球重力
a_earth = -mu_earth / r**3 * r_vec
# 月の摂動
r_moon_vec = moon_position(t)
r_sat_moon = r_moon_vec - r_vec
d = np.linalg.norm(r_sat_moon)
r3 = np.linalg.norm(r_moon_vec)
a_moon = mu_moon * (r_sat_moon / d**3 - r_moon_vec / r3**3)
a_total = a_earth + a_moon
return np.concatenate([v_vec, a_total])
# --- 初期状態ベクトル ---
r0, v0 = kep_to_cart(a0, e0, i0, RAAN0, omega0, nu0, mu_earth)
state0 = np.concatenate([r0, v0])
# --- 数値積分（5年間） ---
T_years = 5
t_span = (0, T_years * 365.25 * 86400)
t_eval = np.linspace(t_span[0], t_span[1], 5000)
print("数値積分を実行中...")
sol = solve_ivp(equations_of_motion, t_span, state0,
method='DOP853', t_eval=t_eval,
rtol=1e-12, atol=1e-14, max_step=3600)
print(f"完了: {sol.success}")
積分結果から各時刻の軌道要素を計算し、グラフ化します。
# --- 軌道要素の時間変化を計算 ---
N = len(sol.t)
a_hist = np.zeros(N)
e_hist = np.zeros(N)
i_hist = np.zeros(N)
RAAN_hist = np.zeros(N)
omega_hist = np.zeros(N)
for k in range(N):
r_k = sol.y[0:3, k]
v_k = sol.y[3:6, k]
a_hist[k], e_hist[k], i_hist[k], RAAN_hist[k], omega_hist[k], _ = \
cart_to_kep(r_k, v_k, mu_earth)
t_years = sol.t / (365.25 * 86400)
# --- 可視化 ---
fig, axes = plt.subplots(4, 1, figsize=(10, 12), sharex=True)
axes[0].plot(t_years, (a_hist - a0) / 1e3, color='steelblue', linewidth=0.5)
axes[0].set_ylabel('Δa [km]')
axes[0].set_title('第三体摂動（月）によるGEO衛星の軌道要素変化')
axes[0].grid(True, alpha=0.3)
axes[1].plot(t_years, e_hist, color='darkorange', linewidth=0.5)
axes[1].set_ylabel('離心率 e')
axes[1].grid(True, alpha=0.3)
axes[2].plot(t_years, np.degrees(i_hist), color='green', linewidth=0.5)
axes[2].set_ylabel('傾斜角 i [deg]')
axes[2].grid(True, alpha=0.3)
axes[3].plot(t_years, np.degrees(RAAN_hist), color='purple', linewidth=0.5)
axes[3].set_ylabel('RAAN Ω [deg]')
axes[3].set_xlabel('時間 [年]')
axes[3].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('third_body_perturbation_geo.png', dpi=150, bbox_inches='tight')
plt.show()
このシミュレーション結果から、いくつかの重要な特徴が読み取れます。
-
半長軸 $a$ はほぼ一定 — 数値的に確認すると、半長軸の変動は数km程度の短周期振動にとどまり、長期的なトレンドはありません。これは理論で予測した「第三体摂動は半長軸を永年的に変えない」という結果と一致しています。
-
離心率 $e$ は月の公転周期に対応する振動を示す — 約27日周期の短周期変動に加え、より長い周期の変動が重畳しています。振幅は $10^{-3}$ オーダーで、GEO衛星の運用にはそれほど影響しません。
-
傾斜角 $i$ が長期的に変化する — 初期値23°から数年で1°以上の変化が生じます。実際のGEO衛星では月と太陽の両方の摂動が組み合わさり、年間約0.75°の傾斜角増加が生じます。
-
昇交点赤経 $\Omega$ が永年的に回転する — 長期的なドリフトに短周期の振動が重畳した振る舞いを示しています。
これらの数値的な結果は、前節で導出した解析理論と定性的に良く一致しています。
コズライ共鳴のシミュレーション
コズライ共鳴をより明確に観察するため、高傾斜角（初期傾斜角60°）の軌道でシミュレーションを行います。
# --- コズライ共鳴のシミュレーション ---
# 高傾斜角の軌道（i=60°）で長期積分
i0_kozai = np.radians(60.0)
e0_kozai = 0.01
r0_k, v0_k = kep_to_cart(a0, e0_kozai, i0_kozai, RAAN0, omega0, nu0, mu_earth)
state0_k = np.concatenate([r0_k, v0_k])
T_years_k = 20
t_span_k = (0, T_years_k * 365.25 * 86400)
t_eval_k = np.linspace(t_span_k[0], t_span_k[1], 10000)
print("コズライ共鳴シミュレーション実行中...")
sol_k = solve_ivp(equations_of_motion, t_span_k, state0_k,
method='DOP853', t_eval=t_eval_k,
rtol=1e-12, atol=1e-14, max_step=3600)
print(f"完了: {sol_k.success}")
# 軌道要素の計算
N_k = len(sol_k.t)
e_kozai = np.zeros(N_k)
i_kozai = np.zeros(N_k)
omega_kozai = np.zeros(N_k)
for k in range(N_k):
_, e_kozai[k], i_kozai[k], _, omega_kozai[k], _ = \
cart_to_kep(sol_k.y[0:3, k], sol_k.y[3:6, k], mu_earth)
t_years_kozai = sol_k.t / (365.25 * 86400)
# コズライ積分の検証
H_kozai = np.sqrt(1 - e_kozai**2) * np.cos(i_kozai)
fig, axes = plt.subplots(3, 1, figsize=(10, 10), sharex=True)
axes[0].plot(t_years_kozai, e_kozai, color='darkorange', linewidth=0.5)
axes[0].set_ylabel('離心率 e')
axes[0].set_title('コズライ共鳴（初期傾斜角60°, GEO軌道）')
axes[0].grid(True, alpha=0.3)
axes[1].plot(t_years_kozai, np.degrees(i_kozai), color='green', linewidth=0.5)
axes[1].set_ylabel('傾斜角 i [deg]')
axes[1].grid(True, alpha=0.3)
axes[2].plot(t_years_kozai, H_kozai, color='crimson', linewidth=0.5)
axes[2].set_ylabel('コズライ積分 H')
axes[2].set_xlabel('時間 [年]')
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('kozai_resonance.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"Hの変動幅: {H_kozai.max() - H_kozai.min():.6f}")
print(f"Hの初期値: {H_kozai[0]:.6f}")
このグラフから、コズライ共鳴の特徴的な振る舞いが確認できます。
-
離心率と傾斜角が逆相関で振動する — 離心率が増加するとき傾斜角は減少し、その逆も成り立ちます。これはコズライ積分 $H = \sqrt{1-e^2}\cos i$ の保存から必然的に従う関係です。
-
コズライ積分 $H$ はほぼ一定 — 数値積分の誤差を除いて、$H$ の値は高い精度で保存されています。これは理論的な保存量が数値計算でも確認できることを示しており、シミュレーションの信頼性を裏付けます。
-
振動の時間スケールは数年〜数十年 — 実際のGEO衛星や廃棄軌道に投入された衛星にとって、このタイムスケールは無視できない長さです。
コズライ共鳴は人工衛星だけでなく、太陽系外惑星の軌道進化や連星系の進化にも重要な役割を果たしている普遍的な力学現象です。次に、これらの知見が実際のGEO衛星の運用にどう影響するかを見ていきましょう。
GEO衛星への影響と軌道維持
傾斜角の永年変化
GEO衛星にとって月・太陽の第三体摂動で最も問題になるのは、傾斜角の長期的な増加です。月と太陽の摂動を合わせると、GEO衛星の傾斜角は年間約 $0.75°$〜$0.95°$ の割合で増加します。
初期傾斜角がゼロのGEO衛星の場合、傾斜角の年間変化率は次のように近似できます。
$$ \frac{di}{dt} \approx \frac{3}{4}\frac{n_3^2}{n}(1 + \tfrac{3}{2}e_3^2)\sin 2\epsilon_3 $$
ここで $n_3$ は第三体の平均運動、$\epsilon_3$ は第三体の軌道面と赤道面のなす角です。月の場合 $\epsilon_{\text{moon}} \approx 18.3°$〜$28.6°$（18.6年周期で変動）、太陽の場合 $\epsilon_{\odot} \approx 23.44°$（黄道傾斜角）です。
南北ステーションキーピング
傾斜角が増加すると、衛星は静止軌道面から上下にずれて、地上のアンテナから見ると南北に振動して見えます。通信衛星のビーム幅は限られているため、傾斜角が許容範囲（通常0.05°〜0.1°）を超えないように定期的に軌道修正を行う必要があります。
この傾斜角修正のための南北ステーションキーピングに必要な年間 $\Delta V$ は
$$ \Delta V_{\text{NS}} \approx v_c \cdot \Delta i_{\text{年間}} \approx 3.07 \, \text{km/s} \times 0.013 \, \text{rad} \approx 40\sim50 \, \text{m/s/year} $$
これは東西ステーションキーピング（J2と月・太陽による経度方向ドリフト修正、年間約2 m/s）に比べて桁違いに大きく、GEO衛星の総 $\Delta V$ バジェットの大部分を占めます。15年のミッション寿命で600〜750 m/s もの $\Delta V$ が必要です。
廃棄軌道と長期軌道進化
ミッション終了後のGEO衛星は、通常300 km以上高い「墓場軌道」（graveyard orbit）に移されます。しかし、墓場軌道に投入された衛星も月・太陽摂動を受け続けます。コズライ共鳴によって離心率が増加すると、近地点高度がGEOベルトまで下がり、運用中の衛星と衝突するリスクが生じます。
このため、墓場軌道の設計では初期の $\omega$（近地点引数）を適切に選び、離心率の長期的な増大を抑制するデザインが研究されています。
実際のミッション設計ではこうした定性的な理解に加えて、数値的に正確な長期軌道伝播が必要になります。次の記事では、複数の摂動源を組み合わせた長期軌道伝播の手法を詳しく解説します。
まとめ
本記事では、月・太陽による第三体摂動の理論を解説しました。
- 二体問題の限界: 現実の衛星は地球の重力だけでなく、月や太陽からの重力（第三体摂動）も受ける。特に高軌道では第三体摂動が支配的な摂動源となる
- 摂動力の構造: 第三体摂動は「衛星位置での第三体の重力」と「地球中心での第三体の重力」の差（潮汐力）として表される
- ルジャンドル展開: 摂動ポテンシャルをルジャンドル多項式で展開し、$n=2$ の主要項が $r^2/r_3^3$ に比例することを示した。高軌道ほど影響が大きい
- 平均化法: 短周期変動を取り除き、長期的な軌道進化を解析する。半長軸は第三体摂動で永年変化しないが、傾斜角と昇交点赤経は変化する
- コズライ共鳴: 傾斜角が約39.2°を超える軌道では、離心率と傾斜角がエネルギーを交換しながら大きく振動する。コズライ積分 $H = \sqrt{1-e^2}\cos i$ が保存量となる
- GEO衛星への影響: 傾斜角が年間約0.75°増加し、南北ステーションキーピングに年間40〜50 m/sのΔVが必要
次のステップとして、以下の記事も参考にしてください。