ケプラー軌道は二体問題の解ですが、実際の人工衛星軌道は地球の重力場の非球対称性によって摂動を受けます。最も支配的な摂動は地球の扁平性（oblateness）に起因する $J_2$ 摂動です。この摂動により、軌道面が歳差運動し、近地点が回転します。
本記事の内容
- 地球の扁平性と $J_2$ 係数
- 摂動ポテンシャル
- RAANの歳差 $\dot{\Omega}$
- 近地点引数の回転 $\dot{\omega}$
- 太陽同期軌道の条件
- Pythonでの計算と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
地球の扁平性と $J_2$ 係数
地球は完全な球ではなく、自転による遠心力で赤道方向に膨らんだ回転楕円体です。赤道半径と極半径の差は約 21 km です。
地球の重力ポテンシャルは球面調和関数で展開できます:
$$ U = -\frac{\mu}{r}\left[1 – \sum_{n=2}^{\infty} J_n \left(\frac{R_E}{r}\right)^n P_n(\sin\phi)\right] $$
ここで $P_n$ はルジャンドル多項式、$\phi$ は地心緯度、$R_E$ は地球の赤道半径です。
支配的な項は $n = 2$ の帯状調和 $J_2$ です:
$$ \boxed{J_2 = 1.08263 \times 10^{-3}} $$
$J_2$ は $J_3$ の約1000倍であり、非球形重力の主要な摂動源です。
$J_2$ 摂動ポテンシャル
$J_2$ 項のみを取り出すと:
$$ U_{J_2} = \frac{\mu J_2 R_E^2}{2r^3}(3\sin^2\phi – 1) $$
$\sin\phi = \sin i \sin(\omega + \nu)$ を用いて:
$$ U_{J_2} = \frac{\mu J_2 R_E^2}{2r^3}\left[3\sin^2 i \sin^2(\omega + \nu) – 1\right] $$
軌道平均による長期変化
摂動の効果を1軌道周期で平均化すると、長期的な（secular）軌道要素の変化率が得られます。ラグランジュの惑星方程式を用いた結果:
半長軸と離心率
$$ \dot{a}_\text{sec} = 0, \quad \dot{e}_\text{sec} = 0 $$
軌道の形状は $J_2$ の長期摂動では変化しません。
RAAN（昇交点赤経）の歳差
$$ \boxed{\dot{\Omega} = -\frac{3}{2} n J_2 \left(\frac{R_E}{a}\right)^2 \frac{\cos i}{(1-e^2)^2}} $$
ここで $n = \sqrt{\mu/a^3}$ は平均運動（平均角速度）です。
$\dot{\Omega}$ の特徴:
- $i < 90°$（順行軌道）: $\dot{\Omega} < 0$（西向きに歳差）
- $i > 90°$（逆行軌道）: $\dot{\Omega} > 0$（東向きに歳差）
- $i = 90°$（極軌道）: $\dot{\Omega} = 0$
近地点引数の回転
$$ \boxed{\dot{\omega} = \frac{3}{2} n J_2 \left(\frac{R_E}{a}\right)^2 \frac{4 – 5\sin^2 i}{2(1-e^2)^2}} $$
$\dot{\omega} = 0$ となる条件:
$$ \sin^2 i = \frac{4}{5} \implies i = 63.4° \text{ or } 116.6° $$
この傾斜角をクリティカル傾斜角といい、モルニヤ軌道やツンドラ軌道で利用されます。
平均運動の補正
$$ \dot{M}_\text{sec} = n\left[1 + \frac{3}{2} J_2 \left(\frac{R_E}{a}\right)^2 \frac{1 – \frac{3}{2}\sin^2 i}{(1-e^2)^{3/2}}\right] $$
太陽同期軌道
太陽同期軌道（SSO）は、RAANの歳差率が地球の公転角速度と一致する軌道です:
$$ \dot{\Omega} = \dot{\Omega}_\text{sun} = \frac{2\pi}{365.25 \times 86400} \approx 1.991 \times 10^{-7} \text{ rad/s} $$
条件:
$$ -\frac{3}{2} n J_2 \left(\frac{R_E}{a}\right)^2 \frac{\cos i}{(1-e^2)^2} = \dot{\Omega}_\text{sun} $$
円軌道（$e = 0$）の場合:
$$ \boxed{\cos i = -\frac{2\dot{\Omega}_\text{sun}}{3 J_2 R_E^2} \cdot \frac{a^{7/2}}{\sqrt{\mu}}} $$
典型的なSSO:
- 高度 600 km → $i \approx 97.8°$
- 高度 800 km → $i \approx 98.6°$
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
# 定数
mu = 3.986004418e14 # [m^3/s^2]
R_E = 6378.137e3 # 赤道半径 [m]
J2 = 1.08263e-3
omega_sun = 2*np.pi / (365.25 * 86400) # 地球公転角速度 [rad/s]
def j2_secular_rates(a, e, i, mu, R_E, J2):
"""J2摂動による長期変化率を計算"""
n = np.sqrt(mu / a**3)
p2 = (1 - e**2)**2
factor = 1.5 * n * J2 * (R_E / a)**2
Omega_dot = -factor * np.cos(i) / p2
omega_dot = factor * (2 - 2.5 * np.sin(i)**2) / p2
return Omega_dot, omega_dot
# LEO (高度400km, i=51.6°)の例
a_leo = R_E + 400e3
i_leo = np.radians(51.6)
Od, wd = j2_secular_rates(a_leo, 0.0, i_leo, mu, R_E, J2)
print("=== LEO (400km, i=51.6°) ===")
print(f"RAAN歳差 Ω̇ = {np.degrees(Od):.4f} °/day = {np.degrees(Od)*365.25:.2f} °/year")
print(f"近地点回転 ω̇ = {np.degrees(wd):.4f} °/day")
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (a) RAAN歳差率 vs 軌道傾斜角
i_range = np.linspace(0, np.pi, 500)
altitudes = [400, 600, 800, 1200]
for alt in altitudes:
a = R_E + alt * 1e3
Od_vals, _ = j2_secular_rates(a, 0.0, i_range, mu, R_E, J2)
Od_deg_day = np.degrees(Od_vals) * 86400
axes[0, 0].plot(np.degrees(i_range), Od_deg_day, linewidth=2,
label=f'h = {alt} km')
axes[0, 0].axhline(np.degrees(omega_sun)*86400, color='red', linestyle='--',
linewidth=1.5, label='太陽同期条件')
axes[0, 0].axhline(0, color='black', linewidth=0.5)
axes[0, 0].set_xlabel('軌道傾斜角 $i$ [deg]')
axes[0, 0].set_ylabel('$\\dot{\\Omega}$ [deg/day]')
axes[0, 0].set_title('RAAN歳差率')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)
# (b) 近地点引数の回転率 vs 傾斜角
for alt in altitudes:
a = R_E + alt * 1e3
_, wd_vals = j2_secular_rates(a, 0.0, i_range, mu, R_E, J2)
wd_deg_day = np.degrees(wd_vals) * 86400
axes[0, 1].plot(np.degrees(i_range), wd_deg_day, linewidth=2,
label=f'h = {alt} km')
axes[0, 1].axhline(0, color='black', linewidth=0.5)
axes[0, 1].axvline(63.4, color='red', linestyle='--', alpha=0.7, label='63.4° (critical)')
axes[0, 1].axvline(116.6, color='red', linestyle='--', alpha=0.7)
axes[0, 1].set_xlabel('軌道傾斜角 $i$ [deg]')
axes[0, 1].set_ylabel('$\\dot{\\omega}$ [deg/day]')
axes[0, 1].set_title('近地点引数の回転率')
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)
# (c) 太陽同期軌道: 高度 vs 傾斜角
alt_range = np.linspace(200, 1500, 500)
a_range = R_E + alt_range * 1e3
n_range = np.sqrt(mu / a_range**3)
# cosφ = -2*omega_sun / (3*J2*R_E^2) * a^(7/2) / sqrt(mu)
cos_i_sso = -2 * omega_sun * a_range**(3.5) / (3 * J2 * R_E**2 * np.sqrt(mu))
valid = np.abs(cos_i_sso) <= 1
i_sso = np.degrees(np.arccos(cos_i_sso[valid]))
axes[1, 0].plot(alt_range[valid], i_sso, 'b-', linewidth=2)
axes[1, 0].set_xlabel('軌道高度 [km]')
axes[1, 0].set_ylabel('軌道傾斜角 $i$ [deg]')
axes[1, 0].set_title('太陽同期軌道の条件')
axes[1, 0].grid(True, alpha=0.3)
# 代表的な高度をマーク
for h_mark, label in [(500, '500km'), (600, '600km'), (800, '800km')]:
a_m = R_E + h_mark * 1e3
cos_i_m = -2 * omega_sun * a_m**3.5 / (3 * J2 * R_E**2 * np.sqrt(mu))
if abs(cos_i_m) <= 1:
i_m = np.degrees(np.arccos(cos_i_m))
axes[1, 0].plot(h_mark, i_m, 'ro', markersize=8)
axes[1, 0].annotate(f'{label}: {i_m:.1f}°', (h_mark, i_m),
textcoords="offset points", xytext=(10, -10), fontsize=10)
# (d) RAAN歳差の時間発展
t_days = np.linspace(0, 365, 1000)
orbits = [
(R_E+400e3, 0.0, np.radians(51.6), 'ISS-like (51.6°)'),
(R_E+800e3, 0.0, np.radians(98.6), 'SSO (98.6°)'),
(R_E+400e3, 0.0, np.radians(90), '極軌道 (90°)'),
]
for a, e, inc, name in orbits:
Od_val, _ = j2_secular_rates(a, e, inc, mu, R_E, J2)
Omega_t = np.degrees(Od_val * t_days * 86400)
axes[1, 1].plot(t_days, Omega_t, linewidth=2, label=name)
axes[1, 1].plot(t_days, np.degrees(omega_sun * t_days * 86400), 'k--',
linewidth=1.5, label='地球公転')
axes[1, 1].set_xlabel('時間 [日]')
axes[1, 1].set_ylabel('$\\Delta\\Omega$ [deg]')
axes[1, 1].set_title('RAAN歳差の時間発展')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('j2_perturbation.png', dpi=150, bbox_inches='tight')
plt.show()
まとめ
本記事では、$J_2$ 摂動と軌道要素の長期変化について解説しました。
- $J_2 = 1.08263 \times 10^{-3}$: 地球の扁平性を表す最も重要な摂動項
- RAAN歳差: $\dot{\Omega} = -\frac{3}{2}nJ_2(R_E/a)^2 \cos i/(1-e^2)^2$
- 近地点回転: $\dot{\omega} = \frac{3}{2}nJ_2(R_E/a)^2(2-2.5\sin^2 i)/(1-e^2)^2$
- クリティカル傾斜角: $i = 63.4°$ で $\dot{\omega} = 0$（モルニヤ軌道）
- 太陽同期軌道: $\dot{\Omega} = \dot{\Omega}_\text{sun}$ で軌道面が太陽と一定の角度を保つ
- 高度600 kmでSSO傾斜角は約97.8°
次のステップとして、以下の記事も参考にしてください。