人工衛星や惑星の軌道を記述するには、6つの独立なパラメータが必要です。古典的軌道要素（Keplerian elements）は、軌道の形状・向き・衛星の位置を幾何学的に直感的に表す6つのパラメータです。
本記事の内容
- 古典的軌道6要素 $(a, e, i, \Omega, \omega, \nu)$ の定義
- 各要素の幾何学的意味
- 軌道要素から位置・速度ベクトルへの変換
- ECI座標系への回転行列
- Pythonでの実装と3D可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
古典的軌道6要素
軌道の形状を決める2要素
半長軸 $a$: 楕円軌道の大きさを決めます。円軌道では半径に等しくなります。
$$ a = \frac{r_p + r_a}{2} $$
ここで $r_p$ は近地点距離、$r_a$ は遠地点距離です。
離心率 $e$: 軌道の形状（扁平さ）を決めます。
$$ e = \frac{r_a – r_p}{r_a + r_p} $$
$e = 0$ で円、$0 < e < 1$ で楕円、$e = 1$ で放物線、$e > 1$ で双曲線です。
軌道面の向きを決める2要素
軌道傾斜角 $i$: 軌道面と赤道面のなす角度です。
- $i = 0°$: 赤道軌道
- $i = 90°$: 極軌道
- $i > 90°$: 逆行軌道
昇交点赤経 $\Omega$（RAAN）: 春分点方向から昇交点（衛星が南→北に赤道面を横切る点）までの角度です。
軌道面内の向きと衛星位置を決める2要素
近地点引数 $\omega$: 昇交点から近地点までの軌道面内の角度です。
真近点角 $\nu$（true anomaly）: 近地点から衛星までの軌道面内の角度です。時間とともに変化する唯一の要素です。
軌道面内の位置・速度
軌道方程式より、動径距離 $r$ は:
$$ r = \frac{a(1 – e^2)}{1 + e\cos\nu} $$
近地点座標系（PQW座標系）での位置ベクトル:
$$ \bm{r}_\text{PQW} = \begin{pmatrix} r\cos\nu \\ r\sin\nu \\ 0 \end{pmatrix} $$
速度ベクトルは、比角運動量 $h = \sqrt{\mu a(1-e^2)}$（$\mu = GM$）を用いて:
$$ \bm{v}_\text{PQW} = \frac{\mu}{h}\begin{pmatrix} -\sin\nu \\ e + \cos\nu \\ 0 \end{pmatrix} $$
ECI座標系への変換
PQW座標系からECI（Earth-Centered Inertial）座標系への変換は、3つの回転行列の積で表されます:
$$ \bm{R} = R_3(-\Omega) \cdot R_1(-i) \cdot R_3(-\omega) $$
ここで $R_1(\theta)$ は $x$ 軸周りの回転、$R_3(\theta)$ は $z$ 軸周りの回転です:
$$ R_3(\theta) = \begin{pmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad R_1(\theta) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos\theta & -\sin\theta \\ 0 & \sin\theta & \cos\theta \end{pmatrix} $$
変換行列の各成分を展開すると:
$$ \bm{R} = \begin{pmatrix} \cos\Omega\cos\omega – \sin\Omega\sin\omega\cos i & -\cos\Omega\sin\omega – \sin\Omega\cos\omega\cos i & \sin\Omega\sin i \\ \sin\Omega\cos\omega + \cos\Omega\sin\omega\cos i & -\sin\Omega\sin\omega + \cos\Omega\cos\omega\cos i & -\cos\Omega\sin i \\ \sin\omega\sin i & \cos\omega\sin i & \cos i \end{pmatrix} $$
ECI座標系での位置・速度:
$$ \bm{r}_\text{ECI} = \bm{R} \cdot \bm{r}_\text{PQW}, \quad \bm{v}_\text{ECI} = \bm{R} \cdot \bm{v}_\text{PQW} $$
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 地球の重力定数
mu = 3.986004418e14 # [m^3/s^2]
R_earth = 6371e3 # [m]
def orbital_elements_to_eci(a, e, i, Omega, omega, nu, mu):
"""軌道要素からECI座標系の位置・速度ベクトルを計算"""
# 半通径と比角運動量
p = a * (1 - e**2)
h = np.sqrt(mu * p)
# 動径距離
r = p / (1 + e * np.cos(nu))
# PQW座標系での位置・速度
r_pqw = np.array([r * np.cos(nu), r * np.sin(nu), 0.0])
v_pqw = (mu / h) * np.array([-np.sin(nu), e + np.cos(nu), 0.0])
# 回転行列
cO, sO = np.cos(Omega), np.sin(Omega)
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 ]
])
r_eci = R @ r_pqw
v_eci = R @ v_pqw
return r_eci, v_eci
# ISS風の軌道パラメータ
a_iss = R_earth + 408e3 # 半長軸 [m]
e_iss = 0.0007 # 離心率
i_iss = np.radians(51.6) # 軌道傾斜角
Omega_iss = np.radians(30) # RAAN
omega_iss = np.radians(60) # 近地点引数
# 真近点角を0~360°で変化させて軌道を描画
nu_array = np.linspace(0, 2*np.pi, 500)
positions = np.array([
orbital_elements_to_eci(a_iss, e_iss, i_iss, Omega_iss, omega_iss, nu, mu)[0]
for nu in nu_array
])
fig = plt.figure(figsize=(16, 6))
# (a) 3D軌道
ax1 = fig.add_subplot(131, projection='3d')
ax1.plot(positions[:, 0]/1e6, positions[:, 1]/1e6, positions[:, 2]/1e6,
'b-', linewidth=2, label='ISS-like orbit')
# 地球の描画
u = np.linspace(0, 2*np.pi, 50)
v = np.linspace(0, np.pi, 50)
xe = R_earth/1e6 * np.outer(np.cos(u), np.sin(v))
ye = R_earth/1e6 * np.outer(np.sin(u), np.sin(v))
ze = R_earth/1e6 * np.outer(np.ones_like(u), np.cos(v))
ax1.plot_surface(xe, ye, ze, alpha=0.3, color='cyan')
ax1.set_xlabel('X [×1000 km]')
ax1.set_ylabel('Y [×1000 km]')
ax1.set_zlabel('Z [×1000 km]')
ax1.set_title('ISS軌道（3D）')
ax1.legend()
# (b) 異なる軌道傾斜角
ax2 = fig.add_subplot(132, projection='3d')
inc_vals = [0, 28.5, 51.6, 90, 98.7]
colors = ['red', 'orange', 'blue', 'green', 'purple']
for inc, col in zip(inc_vals, colors):
pos = np.array([
orbital_elements_to_eci(a_iss, 0.0, np.radians(inc), 0, 0, nu, mu)[0]
for nu in nu_array
])
ax2.plot(pos[:, 0]/1e6, pos[:, 1]/1e6, pos[:, 2]/1e6,
color=col, linewidth=1.5, label=f'i={inc}°')
ax2.plot_surface(xe, ye, ze, alpha=0.2, color='cyan')
ax2.set_title('軌道傾斜角の比較')
ax2.legend(fontsize=8)
# (c) 異なる離心率
ax3 = fig.add_subplot(133)
e_vals = [0.0, 0.2, 0.5, 0.8]
for e_val in e_vals:
a_val = 2 * R_earth
r_vals = a_val * (1 - e_val**2) / (1 + e_val * np.cos(nu_array))
x_vals = r_vals * np.cos(nu_array) / R_earth
y_vals = r_vals * np.sin(nu_array) / R_earth
ax3.plot(x_vals, y_vals, linewidth=2, label=f'e = {e_val}')
theta_earth = np.linspace(0, 2*np.pi, 100)
ax3.plot(np.cos(theta_earth), np.sin(theta_earth), 'c-', linewidth=2, alpha=0.5)
ax3.fill(np.cos(theta_earth), np.sin(theta_earth), alpha=0.2, color='cyan')
ax3.plot(0, 0, 'ko', markersize=5)
ax3.set_xlabel('x / $R_E$')
ax3.set_ylabel('y / $R_E$')
ax3.set_title('離心率と軌道形状')
ax3.set_aspect('equal')
ax3.legend()
ax3.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('orbital_elements.png', dpi=150, bbox_inches='tight')
plt.show()
まとめ
本記事では、軌道の6要素と座標変換について解説しました。
- 6要素: $(a, e, i, \Omega, \omega, \nu)$ で軌道の形状・向き・衛星位置を完全に記述
- $a$（半長軸）, $e$（離心率）: 軌道の形状
- $i$（傾斜角）, $\Omega$（RAAN）: 軌道面の向き
- $\omega$（近地点引数）, $\nu$（真近点角）: 軌道面内の向きと位置
- PQW→ECI変換は3つの回転行列 $R_3(-\Omega) R_1(-i) R_3(-\omega)$ の積
次のステップとして、以下の記事も参考にしてください。