人工衛星や宇宙機の姿勢を表現する方法はいくつかあります。
代表的なのは、基準座標系に対して、宇宙機に固定された機体座標系の自由度3の回転で表現する方法です。姿勢表現として代表的なものには、回転行列、オイラー角、クォータニオンを用いたものがあります。
本記事の内容
- 人工衛星で使用される座標系の定義
- 各座標系間の変換
- 姿勢運動学方程式
- Pythonでの座標変換シミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
人工衛星で使用される座標系
ECI座標系（地球中心慣性座標系）
ECI（Earth-Centered Inertial）座標系は、地球の中心を原点とし、恒星に対して固定された座標系です。
- 原点: 地球中心
- $X$ 軸: 春分点方向
- $Z$ 軸: 地球の自転軸（北極方向）
- $Y$ 軸: $Z \times X$ で右手系を構成
ニュートンの運動法則が直接適用できる慣性座標系です。
ECEF座標系（地球中心地球固定座標系）
ECEF（Earth-Centered Earth-Fixed）座標系は、地球に固定された座標系です。
- 原点: 地球中心
- $X$ 軸: 経度0度（グリニッジ子午線）と赤道の交点方向
- $Z$ 軸: 地球の自転軸（北極方向）
- $Y$ 軸: 右手系
ECIからECEFへの変換は、地球の自転角 $\theta_G$ を用いた $Z$ 軸周りの回転です。
$$ \bm{R}_{\text{ECEF}} = \bm{R}_z(\theta_G) \bm{r}_{\text{ECI}} $$
$$ \theta_G = \theta_{G0} + \omega_E (t – t_0) $$
ここで $\omega_E = 7.2921 \times 10^{-5}$ rad/s は地球の自転角速度です。
LVLH座標系（軌道座標系）
LVLH（Local Vertical Local Horizontal）座標系は、衛星の重心を原点とし、軌道に追従する座標系です。
- $x$ 軸: 速度方向（進行方向）
- $y$ 軸: 軌道面法線方向（角運動量方向）
- $z$ 軸: 天底方向（地球中心向き）
機体座標系
衛星の構体に固定された座標系です。姿勢が理想状態のときLVLH座標系と一致するように定義されます。
座標変換
ECI → LVLH変換
衛星のECI座標系での位置 $\bm{r}$ と速度 $\bm{v}$ から、LVLH座標系の基底ベクトルを構成します。
$$ \hat{\bm{z}}_L = -\frac{\bm{r}}{|\bm{r}|}, \quad \hat{\bm{y}}_L = -\frac{\bm{r} \times \bm{v}}{|\bm{r} \times \bm{v}|}, \quad \hat{\bm{x}}_L = \hat{\bm{y}}_L \times \hat{\bm{z}}_L $$
ECI→LVLH変換の回転行列は、
$$ \bm{R}_{L/I} = \begin{pmatrix} \hat{\bm{x}}_L^T \\ \hat{\bm{y}}_L^T \\ \hat{\bm{z}}_L^T \end{pmatrix} $$
姿勢行列
LVLH座標系から機体座標系への変換行列が姿勢行列 $\bm{R}_{B/L}$ です。
$$ \bm{r}_{\text{body}} = \bm{R}_{B/L} \bm{r}_{\text{LVLH}} $$
姿勢運動学方程式
クォータニオンの時間微分
クォータニオン $\bm{q}$ の時間微分は、角速度 $\bm{\omega}$ を用いて次のように表されます。
$$ \dot{\bm{q}} = \frac{1}{2} \bm{\Omega}(\bm{\omega}) \bm{q} $$
ここで、
$$ \bm{\Omega}(\bm{\omega}) = \begin{pmatrix} 0 & -\omega_x & -\omega_y & -\omega_z \\ \omega_x & 0 & \omega_z & -\omega_y \\ \omega_y & -\omega_z & 0 & \omega_x \\ \omega_z & \omega_y & -\omega_x & 0 \end{pmatrix} $$
Pythonでの座標変換シミュレーション
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# 定数
mu = 3.986004418e14 # 地球の重力定数 [m^3/s^2]
Re = 6378137.0 # 地球の赤道半径 [m]
omega_e = 7.2921e-5 # 地球の自転角速度 [rad/s]
def eci_to_ecef(r_eci, t, theta_G0=0):
"""ECI座標からECEF座標への変換"""
theta_G = theta_G0 + omega_e * t
cos_t = np.cos(theta_G)
sin_t = np.sin(theta_G)
R = np.array([[cos_t, sin_t, 0],
[-sin_t, cos_t, 0],
[0, 0, 1]])
return R @ r_eci
def eci_to_lvlh(r_eci, v_eci):
"""ECI座標からLVLH座標系の基底ベクトルを計算"""
r_norm = r_eci / np.linalg.norm(r_eci)
h = np.cross(r_eci, v_eci)
h_norm = h / np.linalg.norm(h)
z_L = -r_norm # 天底方向
y_L = -h_norm # 軌道面法線の逆
x_L = np.cross(y_L, z_L) # 速度方向
return np.array([x_L, y_L, z_L])
def ecef_to_geodetic(r_ecef):
"""ECEF座標から測地座標（緯度・経度・高度）への変換"""
x, y, z = r_ecef
lon = np.arctan2(y, x)
p = np.sqrt(x**2 + y**2)
lat = np.arctan2(z, p)
alt = np.sqrt(x**2 + y**2 + z**2) - Re
return np.degrees(lat), np.degrees(lon), alt
# 軌道伝搬（2体問題）
def two_body(t, state):
"""2体問題の運動方程式"""
r = state[:3]
v = state[3:]
r_norm = np.linalg.norm(r)
a = -mu / r_norm**3 * r
return np.concatenate([v, a])
# 初期条件（高度500km、軌道傾斜角51.6度の円軌道）
h = 500e3 # 高度 [m]
a = Re + h
v_circ = np.sqrt(mu / a)
inc = np.radians(51.6)
r0 = np.array([a, 0, 0])
v0 = np.array([0, v_circ * np.cos(inc), v_circ * np.sin(inc)])
state0 = np.concatenate([r0, v0])
# 1軌道周期分のシミュレーション
T_orb = 2 * np.pi * np.sqrt(a**3 / mu)
t_span = (0, T_orb)
t_eval = np.linspace(0, T_orb, 5000)
sol = solve_ivp(two_body, t_span, state0, t_eval=t_eval, method='RK45',
rtol=1e-12, atol=1e-12)
# 座標変換と可視化
lats = []
lons = []
alts = []
for i in range(len(sol.t)):
r_eci = sol.y[:3, i]
r_ecef = eci_to_ecef(r_eci, sol.t[i])
lat, lon, alt = ecef_to_geodetic(r_ecef)
lats.append(lat)
lons.append(lon)
alts.append(alt / 1e3)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# ECI軌道
axes[0, 0].plot(sol.y[0] / Re, sol.y[1] / Re, 'b-', linewidth=0.5)
circle = plt.Circle((0, 0), 1, color='blue', alpha=0.3)
axes[0, 0].add_patch(circle)
axes[0, 0].set_xlabel('X [Earth Radii]')
axes[0, 0].set_ylabel('Y [Earth Radii]')
axes[0, 0].set_title('Orbit in ECI (XY plane)')
axes[0, 0].set_aspect('equal')
axes[0, 0].grid(True)
# 地上軌跡
axes[0, 1].scatter(lons, lats, c=sol.t / 60, cmap='viridis', s=1)
axes[0, 1].set_xlim(-180, 180)
axes[0, 1].set_ylim(-90, 90)
axes[0, 1].set_xlabel('Longitude [deg]')
axes[0, 1].set_ylabel('Latitude [deg]')
axes[0, 1].set_title('Ground Track')
axes[0, 1].grid(True, alpha=0.3)
# 高度の時間変化
axes[1, 0].plot(sol.t / 60, alts, 'b-')
axes[1, 0].set_xlabel('Time [min]')
axes[1, 0].set_ylabel('Altitude [km]')
axes[1, 0].set_title('Altitude')
axes[1, 0].grid(True)
# LVLH座標系の基底ベクトルの時間変化
lvlh_x = []
for i in range(len(sol.t)):
R_lvlh = eci_to_lvlh(sol.y[:3, i], sol.y[3:, i])
lvlh_x.append(R_lvlh[0])
lvlh_x = np.array(lvlh_x)
axes[1, 1].plot(sol.t / 60, lvlh_x[:, 0], label='x')
axes[1, 1].plot(sol.t / 60, lvlh_x[:, 1], label='y')
axes[1, 1].plot(sol.t / 60, lvlh_x[:, 2], label='z')
axes[1, 1].set_xlabel('Time [min]')
axes[1, 1].set_ylabel('LVLH x-axis in ECI')
axes[1, 1].set_title('LVLH Frame Orientation')
axes[1, 1].legend()
axes[1, 1].grid(True)
plt.tight_layout()
plt.show()
print(f"=== 軌道パラメータ ===")
print(f"軌道高度: {h/1e3:.0f} km")
print(f"軌道周期: {T_orb/60:.1f} min")
print(f"軌道速度: {v_circ:.1f} m/s")
まとめ
本記事では、人工衛星の姿勢表現と座標系について解説しました。
- ECI座標系は慣性座標系でニュートンの法則が直接適用できる
- ECEF座標系は地球に固定され、地上の位置と対応付けられる
- LVLH座標系は衛星に追従し、姿勢制御の基準となる
- クォータニオンの運動学方程式 $\dot{\bm{q}} = \frac{1}{2}\bm{\Omega}\bm{q}$ で姿勢の時間発展を記述する
- 各座標系間の変換は回転行列で表現される
次のステップとして、以下の記事も参考にしてください。