2013年2月15日、ロシアのチェリャビンスク上空で直径約17mの小惑星が大気圏に突入し、衝撃波で約1,500人が負傷しました。この小惑星は事前に検出されていませんでした。もし十分早い段階で観測され、軌道が決定されていれば、警報を出すことが可能だったかもしれません。
軌道決定（Orbit Determination, OD）とは、天体や人工衛星の観測データ（位置、角度、距離など）から、その軌道要素（軌道の形状・向き・位置を決める6つのパラメータ）を推定する問題です。軌道を知ることは、衛星の追跡と管制、デブリの衝突回避、小惑星の地球衝突リスク評価など、宇宙活動のあらゆる場面で不可欠です。
軌道決定を理解すると、以下のような応用が見えてきます。
- 宇宙状況認識（SSA）: 軌道上の衛星やデブリの位置を追跡し、衝突を予防する
- 惑星防衛: 地球近傍小天体（NEO）の軌道を早期に決定し、衝突リスクを評価する
- 深宇宙探査: 探査機の航行精度を向上させ、目標天体への正確な到達を実現する
本記事の内容
- 軌道決定問題の定式化
- ガウスの初期軌道決定法（3回の角度観測から軌道要素を推定）
- 最小二乗法による軌道改良
- 観測方程式の線形化とヤコビアン行列
- Pythonでのガウス法の実装とシミュレーション
- 誤差の評価と共分散解析
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
軌道決定問題の全体像
軌道決定とは何か
軌道決定の本質を日常のアナロジーで説明しましょう。夜空に一瞬光った流れ星を見たとき、その軌跡を延長して「どこから来たのか、どこへ行くのか」を推定する作業に似ています。ただし、1回の目撃情報では軌跡の方向しかわかりません。異なる時刻に3回以上観測すれば、3次元的な軌道を特定できます。
数学的に言えば、ケプラー軌道は6つの軌道要素（$a, e, i, \Omega, \omega, \nu_0$ または等価なパラメータ）で完全に決まります。したがって、最低6個の独立した観測量があれば、原理的に軌道を一意に決定できます。
軌道決定の2つの段階
実際の軌道決定は通常、2つの段階で行われます。
第1段階: 初期軌道決定（Initial Orbit Determination, IOD）
少数の観測データ（典型的には3回の角度観測）から、軌道の概略値を求めます。ガウスの方法、ラプラスの方法、ギブスの方法などが知られています。
第2段階: 軌道改良（Orbit Improvement / Differential Correction）
初期軌道を出発点として、多数の観測データを用いて最小二乗法で軌道要素を高精度に推定します。ここでは観測方程式の線形化とヤコビアン行列の計算が必要になります。
初期軌道決定は「ゼロから軌道を見つける」問題であり、軌道改良は「既知の軌道をデータに合わせて微調整する」問題です。前者は非線形で解析的な工夫が必要であり、後者は線形代数の枠組みで系統的に扱えます。
まずは古典的なガウスの方法から始めましょう。ガウスは1801年に小惑星ケレスの発見観測（41日間にわたるわずかな角度観測データ）から軌道を決定し、見失ったケレスの再発見に貢献したことで有名です。
ガウスの初期軌道決定法
問題設定
ガウスの方法は、3回の角度観測（赤経 $\alpha$ と赤緯 $\delta$ の組）から軌道を決定します。地上観測所から天体の方向（視線ベクトル）のみが観測され、距離は未知とします。
観測データ: – 時刻 $t_1, t_2, t_3$ における天体の方向（視線ベクトル）: $\hat{\bm{L}}_1, \hat{\bm{L}}_2, \hat{\bm{L}}_3$ – 同時刻の観測所の位置ベクトル: $\bm{R}_1, \bm{R}_2, \bm{R}_3$
求めるもの: – 天体の軌道要素6個（または等価的に、ある時刻の位置と速度）
視線ベクトルと位置の関係
天体の位置ベクトル $\bm{r}_i$ は、観測所の位置 $\bm{R}_i$ と天体の方向 $\hat{\bm{L}}_i$、そして未知の距離（スラントレンジ）$\rho_i$ を使って次のように書けます。
$$ \bm{r}_i = \bm{R}_i + \rho_i \hat{\bm{L}}_i, \quad (i = 1, 2, 3) $$
ここで $\rho_i > 0$ が3つの未知数です。これらが決まれば3つの位置ベクトルが求まり、位置と速度（ギブスの方法で計算可能）から軌道要素が得られます。
ラグランジュ係数の導入
二体問題の解において、時刻 $t$ の位置ベクトルは基準時刻 $t_2$ の位置ベクトルと速度ベクトルの線形結合で表されます。
$$ \bm{r}_1 = f_1 \bm{r}_2 + g_1 \dot{\bm{r}}_2 $$
$$ \bm{r}_3 = f_3 \bm{r}_2 + g_3 \dot{\bm{r}}_2 $$
ここで $f_i$ と $g_i$ はラグランジュ係数と呼ばれ、時間差 $\tau_i = t_i – t_2$ と軌道パラメータの関数です。この係数の正確な値は軌道がわからないと計算できませんが、短い時間間隔では級数展開で近似できます。
$$ f_1 \approx 1 – \frac{\mu}{2r_2^3}\tau_1^2, \quad g_1 \approx \tau_1 – \frac{\mu}{6r_2^3}\tau_1^3 $$
$$ f_3 \approx 1 – \frac{\mu}{2r_2^3}\tau_3^2, \quad g_3 \approx \tau_3 – \frac{\mu}{6r_2^3}\tau_3^3 $$
基本方程式の導出
$\bm{r}_1$ と $\bm{r}_3$ の式から $\dot{\bm{r}}_2$ を消去すると
$$ \bm{r}_2 = c_1 \bm{r}_1 + c_3 \bm{r}_3 $$
ここで
$$ c_1 = \frac{g_3}{f_1 g_3 – f_3 g_1}, \quad c_3 = -\frac{g_1}{f_1 g_3 – f_3 g_1} $$
$\bm{r}_i = \bm{R}_i + \rho_i \hat{\bm{L}}_i$ を代入すると
$$ \bm{R}_2 + \rho_2 \hat{\bm{L}}_2 = c_1(\bm{R}_1 + \rho_1 \hat{\bm{L}}_1) + c_3(\bm{R}_3 + \rho_3 \hat{\bm{L}}_3) $$
整理すると
$$ c_1 \rho_1 \hat{\bm{L}}_1 – \rho_2 \hat{\bm{L}}_2 + c_3 \rho_3 \hat{\bm{L}}_3 = \bm{R}_2 – c_1 \bm{R}_1 – c_3 \bm{R}_3 $$
右辺は既知量、左辺は3つの未知スラントレンジ $\rho_1, \rho_2, \rho_3$ を含みます。このベクトル方程式は3つのスカラー方程式を与えるので、3つの未知数に対して解くことができます。行列形式で書くと
$$ \begin{pmatrix} L_{1x} c_1 & -L_{2x} & L_{3x} c_3 \\ L_{1y} c_1 & -L_{2y} & L_{3y} c_3 \\ L_{1z} c_1 & -L_{2z} & L_{3z} c_3 \end{pmatrix} \begin{pmatrix} \rho_1 \\ \rho_2 \\ \rho_3 \end{pmatrix} = \bm{R}_2 – c_1 \bm{R}_1 – c_3 \bm{R}_3 $$
反復解法
上の式の $c_1, c_3$ はラグランジュ係数に依存し、ラグランジュ係数は $r_2$（まだ未知）に依存します。そこで、以下の反復手順で解きます。
- $r_2$ の初期推定値を仮定する（たとえば $r_2 = R_{\oplus} + 500 \, \text{km}$）
- ラグランジュ係数の近似値を計算する
- $c_1, c_3$ を計算する
- スラントレンジ $\rho_1, \rho_2, \rho_3$ を求める
- 位置ベクトル $\bm{r}_1, \bm{r}_2, \bm{r}_3$ を計算する
- 新しい $r_2 = |\bm{r}_2|$ で収束判定する
- 収束していなければ手順2に戻る
収束後、3つの位置ベクトルと時刻から軌道要素を決定します（ギブスの方法やヘリック-ギブスの方法を使用）。
ガウス法は初期軌道を与えてくれますが、その精度は観測データの量と質に依存します。より多くの観測データが得られたとき、軌道の精度を向上させるのが次に説明する最小二乗法による軌道改良です。
最小二乗法による軌道改良
問題設定
初期軌道決定で得られた軌道要素を出発点として、$N$ 個の観測データ（$N > 6$）に最も良くフィットするように軌道要素を修正します。
軌道要素ベクトルを $\bm{x} = (a, e, i, \Omega, \omega, M_0)^T$ とし、$i$ 番目の観測量を $z_i$ とします。観測量は軌道要素の非線形関数です。
$$ z_i = h_i(\bm{x}) + \epsilon_i $$
ここで $h_i(\bm{x})$ は観測方程式（軌道要素から観測量を予測する関数）、$\epsilon_i$ は観測誤差です。
線形化
$h_i(\bm{x})$ を現在の軌道推定値 $\bm{x}_0$ の周りでテイラー展開し、1次の項までで打ち切ります。
$$ h_i(\bm{x}_0 + \delta\bm{x}) \approx h_i(\bm{x}_0) + \sum_j \frac{\partial h_i}{\partial x_j}\bigg|_{\bm{x}_0} \delta x_j $$
残差（O-C: Observed minus Computed）を $\Delta z_i = z_i – h_i(\bm{x}_0)$ と定義すると
$$ \Delta z_i \approx \sum_j H_{ij} \delta x_j + \epsilon_i $$
行列形式で書くと
$$ \Delta\bm{z} = \bm{H}\delta\bm{x} + \bm{\epsilon} $$
ここで $\bm{H}$ はヤコビアン行列（偏微分行列）であり、各要素は
$$ H_{ij} = \frac{\partial h_i}{\partial x_j}\bigg|_{\bm{x}_0} $$
最小二乗解
観測誤差の共分散行列を $\bm{W}$ とすると、重み付き最小二乗問題の正規方程式は
$$ (\bm{H}^T \bm{W}^{-1} \bm{H})\delta\bm{x} = \bm{H}^T \bm{W}^{-1} \Delta\bm{z} $$
この解は
$$ \boxed{\delta\bm{x} = (\bm{H}^T \bm{W}^{-1} \bm{H})^{-1} \bm{H}^T \bm{W}^{-1} \Delta\bm{z}} $$
軌道要素を $\bm{x}_0 + \delta\bm{x}$ に更新し、残差が十分小さくなるまで反復します。
推定精度の評価
推定された軌道要素の共分散行列は
$$ \bm{P} = (\bm{H}^T \bm{W}^{-1} \bm{H})^{-1} $$
$\bm{P}$ の対角要素の平方根が各軌道要素の推定誤差の標準偏差を与えます。非対角要素はパラメータ間の相関を表します。
$$ \sigma_{x_j} = \sqrt{P_{jj}} $$
この共分散行列は「どの観測データがどのパラメータの推定に寄与しているか」を定量的に示してくれます。観測が少ない方向のパラメータは誤差が大きくなり、多くの観測がある方向は誤差が小さくなります。
理論の説明が終わったところで、次にPythonを使ってガウスの初期軌道決定法をシミュレーションしましょう。
Pythonによるガウス法のシミュレーション
問題設定
まず、既知の軌道から人工的な観測データを生成し、ガウス法でその軌道を復元できるか検証します。これにより、アルゴリズムの正しさと精度を確認できます。
import numpy as np
def kepler_to_cartesian(a, e, i, Omega, omega, nu, mu):
"""ケプラー軌道要素から直交座標（位置・速度）に変換"""
# 軌道面内の位置と速度
p = a * (1 - e**2)
r = p / (1 + e * np.cos(nu))
# 軌道面座標
r_pqw = np.array([r * np.cos(nu), r * np.sin(nu), 0])
v_pqw = np.sqrt(mu / p) * np.array([-np.sin(nu), e + np.cos(nu), 0])
# 回転行列（PQW → IJK）
cos_O, sin_O = np.cos(Omega), np.sin(Omega)
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_ijk = R @ r_pqw
v_ijk = R @ v_pqw
return r_ijk, v_ijk
def propagate_kepler(r0, v0, dt, mu):
"""二体問題の解析的軌道伝播（ユニバーサル変数法の簡易版）"""
from scipy.integrate import solve_ivp
def eom(t, state):
r_vec = state[:3]
v_vec = state[3:]
r_mag = np.linalg.norm(r_vec)
a_vec = -mu / r_mag**3 * r_vec
return np.concatenate([v_vec, a_vec])
state0 = np.concatenate([r0, v0])
sol = solve_ivp(eom, [0, dt], state0, method='RK45',
rtol=1e-12, atol=1e-14)
return sol.y[:3, -1], sol.y[3:, -1]
# テスト軌道の定義（LEO衛星）
mu_earth = 3.986e14 # m^3/s^2
R_earth = 6378e3 # m
# 軌道要素
a = R_earth + 700e3 # 半長軸: 高度700 km
e_orb = 0.01 # 離心率
inc = np.radians(45) # 傾斜角: 45°
Omega = np.radians(30) # 昇交点赤経: 30°
omega = np.radians(60) # 近地点引数: 60°
nu0 = np.radians(0) # 初期真近点角: 0°
# 初期位置・速度
r0, v0 = kepler_to_cartesian(a, e_orb, inc, Omega, omega, nu0, mu_earth)
print("=== テスト軌道（真値）===")
print(f"半長軸 a = {a/1000:.1f} km")
print(f"離心率 e = {e_orb:.4f}")
print(f"傾斜角 i = {np.degrees(inc):.1f}°")
print(f"初期位置 r0 = {r0/1000} km")
print(f"|r0| = {np.linalg.norm(r0)/1000:.1f} km")
print(f"初期速度 v0 = {v0/1000} km/s")
print(f"|v0| = {np.linalg.norm(v0)/1000:.3f} km/s")
このコードでは、既知の軌道要素からデカルト座標系での位置・速度を計算しています。高度700 km、傾斜角45°のLEO軌道を対象としています。結果として、軌道半径は約7,078 km、軌道速度は約7.5 km/sとなり、LEOの典型的な値と一致します。
観測データの生成とガウス法の実装
次に、この軌道上の3時刻での観測データを生成し、ガウス法で軌道を復元します。
import numpy as np
from scipy.integrate import solve_ivp
def kepler_to_cartesian(a, e, i, Omega, omega, nu, mu):
"""ケプラー軌道要素から直交座標に変換"""
p = a * (1 - e**2)
r = p / (1 + e * np.cos(nu))
r_pqw = np.array([r * np.cos(nu), r * np.sin(nu), 0])
v_pqw = np.sqrt(mu / p) * np.array([-np.sin(nu), e + np.cos(nu), 0])
cos_O, sin_O = np.cos(Omega), np.sin(Omega)
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]
])
return R @ r_pqw, R @ v_pqw
def propagate(r0, v0, dt, mu):
"""二体問題の数値積分"""
def eom(t, state):
rv = state[:3]
vv = state[3:]
rm = np.linalg.norm(rv)
return np.concatenate([vv, -mu / rm**3 * rv])
sol = solve_ivp(eom, [0, dt], np.concatenate([r0, v0]),
method='RK45', rtol=1e-12, atol=1e-14)
return sol.y[:3, -1], sol.y[3:, -1]
# 軌道パラメータ（真値）
mu_earth = 3.986e14
R_earth = 6378e3
a_true = R_earth + 700e3
e_true = 0.01
i_true = np.radians(45)
Om_true = np.radians(30)
w_true = np.radians(60)
nu0_true = np.radians(0)
r0, v0 = kepler_to_cartesian(a_true, e_true, i_true,
Om_true, w_true, nu0_true, mu_earth)
# 3時刻での位置を計算（観測間隔: 10分）
T_orb = 2 * np.pi * np.sqrt(a_true**3 / mu_earth)
dt = 600 # 10分間隔
r1_true = r0.copy()
r2_true, v2_true = propagate(r0, v0, dt, mu_earth)
r3_true, _ = propagate(r0, v0, 2 * dt, mu_earth)
# 観測所の位置（地球中心の慣性系で固定と仮定 — 簡略化）
R_obs = np.array([R_earth, 0, 0]) # 赤道上の観測所
# 視線ベクトル（方向のみ — 距離は未知と仮定）
rho1_true = r1_true - R_obs
rho2_true = r2_true - R_obs
rho3_true = r3_true - R_obs
L1 = rho1_true / np.linalg.norm(rho1_true)
L2 = rho2_true / np.linalg.norm(rho2_true)
L3 = rho3_true / np.linalg.norm(rho3_true)
print("=== 観測データ（真値から生成）===")
print(f"t1 = 0 s, t2 = {dt} s, t3 = {2*dt} s")
print(f"軌道周期 T = {T_orb:.1f} s ({T_orb/60:.1f} min)")
print(f"観測間隔 = {dt} s (周期の {dt/T_orb*100:.1f}%)")
# === ガウス法の実装 ===
tau1 = 0 - dt # t1 - t2
tau3 = 2*dt - dt # t3 - t2
# 初期推定: r2_mag を地球半径+高度の平均で仮定
r2_mag = R_earth + 700e3 # 初期推定値
max_iter = 20
tol = 1e-6
for iteration in range(max_iter):
# ラグランジュ係数の近似
f1 = 1 - mu_earth / (2 * r2_mag**3) * tau1**2
g1 = tau1 - mu_earth / (6 * r2_mag**3) * tau1**3
f3 = 1 - mu_earth / (2 * r2_mag**3) * tau3**2
g3 = tau3 - mu_earth / (6 * r2_mag**3) * tau3**3
d = f1 * g3 - f3 * g1
c1 = g3 / d
c3 = -g1 / d
# 右辺ベクトル
rhs = R_obs - c1 * R_obs - c3 * R_obs # 全時刻で同じ観測所
# 視線方向行列
M = np.column_stack([c1 * L1, -L2, c3 * L3])
# スラントレンジを解く
rho_vec = np.linalg.solve(M, rhs)
rho1, rho2, rho3 = rho_vec
# 位置ベクトル
r1_est = R_obs + rho1 * L1
r2_est = R_obs + rho2 * L2
r3_est = R_obs + rho3 * L3
r2_mag_new = np.linalg.norm(r2_est)
# 収束判定
rel_change = abs(r2_mag_new - r2_mag) / r2_mag
if iteration < 5 or rel_change < 0.01:
print(f"反復 {iteration+1}: |r2| = {r2_mag_new/1000:.2f} km, "
f"変化率 = {rel_change:.2e}")
if rel_change < tol:
print(f"\n収束（{iteration+1} 反復）")
break
r2_mag = r2_mag_new
# 速度の推定（ラグランジュ係数から）
v2_est = (r3_est - f3 * r2_est) / g3
# 結果の表示
print(f"\n=== ガウス法の結果 ===")
print(f"推定 r2 = {r2_est/1000} km")
print(f"真値 r2 = {r2_true/1000} km")
print(f"位置誤差 = {np.linalg.norm(r2_est - r2_true)/1000:.3f} km")
print(f"\n推定 v2 = {v2_est/1000} km/s")
print(f"真値 v2 = {v2_true/1000} km/s")
print(f"速度誤差 = {np.linalg.norm(v2_est - v2_true)*1000:.3f} m/s")
このシミュレーションの結果から、ガウス法の実用的な特性が理解できます。
- 反復は数回で収束します。$r_2$ の初期推定が合理的であれば（正しい高度帯を仮定すれば）、5〜10回の反復で $10^{-6}$ の相対精度に達します。
- 位置の推定精度はkm程度になります。これは観測間隔が軌道周期の数%であり、ラグランジュ係数の級数展開の精度に依存します。観測間隔が短いほど級数展開の精度は上がりますが、短すぎると方向変化が小さく、数値的に不安定になります。
- 速度の推定精度はm/s程度です。速度は2つの位置の差分から求めるため、位置の誤差が増幅されます。
観測誤差の影響評価
実際の観測には誤差があります。観測角度にガウスノイズを加えて、軌道決定の精度がどう変化するかを評価しましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def kepler_to_cartesian(a, e, i, Omega, omega, nu, mu):
"""ケプラー軌道要素から直交座標に変換"""
p = a * (1 - e**2)
r = p / (1 + e * np.cos(nu))
r_pqw = np.array([r * np.cos(nu), r * np.sin(nu), 0])
v_pqw = np.sqrt(mu / p) * np.array([-np.sin(nu), e + np.cos(nu), 0])
cos_O, sin_O = np.cos(Omega), np.sin(Omega)
cos_i, sin_i = np.cos(i), np.sin(i)
cos_w, sin_w = np.cos(omega), np.sin(omega)
R_mat = 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]
])
return R_mat @ r_pqw, R_mat @ v_pqw
def propagate(r0, v0, dt, mu):
"""二体問題の数値積分"""
def eom(t, state):
rv = state[:3]
vv = state[3:]
rm = np.linalg.norm(rv)
return np.concatenate([vv, -mu / rm**3 * rv])
sol = solve_ivp(eom, [0, dt], np.concatenate([r0, v0]),
method='RK45', rtol=1e-12, atol=1e-14)
return sol.y[:3, -1], sol.y[3:, -1]
def gauss_iod(L1, L2, L3, R_obs, tau1, tau3, mu, r2_init):
"""ガウスの初期軌道決定（簡略版）"""
r2_mag = r2_init
for _ in range(30):
f1 = 1 - mu / (2 * r2_mag**3) * tau1**2
g1 = tau1 - mu / (6 * r2_mag**3) * tau1**3
f3 = 1 - mu / (2 * r2_mag**3) * tau3**2
g3 = tau3 - mu / (6 * r2_mag**3) * tau3**3
d = f1 * g3 - f3 * g1
c1 = g3 / d
c3 = -g1 / d
rhs = R_obs - c1 * R_obs - c3 * R_obs
M = np.column_stack([c1 * L1, -L2, c3 * L3])
try:
rho_vec = np.linalg.solve(M, rhs)
except np.linalg.LinAlgError:
return None, None
r2_est = R_obs + rho_vec[1] * L2
r3_est = R_obs + rho_vec[2] * L3
r2_mag_new = np.linalg.norm(r2_est)
if abs(r2_mag_new - r2_mag) / r2_mag < 1e-8:
v2_est = (r3_est - f3 * r2_est) / g3
return r2_est, v2_est
r2_mag = r2_mag_new
return None, None
# 軌道パラメータ（真値）
mu_earth = 3.986e14
R_earth = 6378e3
a_true = R_earth + 700e3
e_true = 0.01
i_true = np.radians(45)
r0, v0 = kepler_to_cartesian(a_true, e_true, i_true,
np.radians(30), np.radians(60),
np.radians(0), mu_earth)
dt = 600 # 観測間隔 10分
r1_true = r0.copy()
r2_true, v2_true = propagate(r0, v0, dt, mu_earth)
r3_true, _ = propagate(r0, v0, 2*dt, mu_earth)
R_obs = np.array([R_earth, 0, 0])
# 真の視線ベクトル
L1_true = (r1_true - R_obs) / np.linalg.norm(r1_true - R_obs)
L2_true = (r2_true - R_obs) / np.linalg.norm(r2_true - R_obs)
L3_true = (r3_true - R_obs) / np.linalg.norm(r3_true - R_obs)
# モンテカルロシミュレーション
noise_levels = [0.1, 0.5, 1.0, 5.0, 10.0] # arcsec
n_trials = 200
results = {}
for noise_arcsec in noise_levels:
noise_rad = noise_arcsec * np.pi / (180 * 3600)
pos_errors = []
vel_errors = []
for trial in range(n_trials):
# 視線ベクトルにノイズを加える
L1_noisy = L1_true + np.random.randn(3) * noise_rad
L2_noisy = L2_true + np.random.randn(3) * noise_rad
L3_noisy = L3_true + np.random.randn(3) * noise_rad
L1_noisy /= np.linalg.norm(L1_noisy)
L2_noisy /= np.linalg.norm(L2_noisy)
L3_noisy /= np.linalg.norm(L3_noisy)
r2_est, v2_est = gauss_iod(L1_noisy, L2_noisy, L3_noisy,
R_obs, -dt, dt, mu_earth,
R_earth + 700e3)
if r2_est is not None:
pos_err = np.linalg.norm(r2_est - r2_true) / 1000 # km
vel_err = np.linalg.norm(v2_est - v2_true) # m/s
if pos_err < 1000: # 外れ値除外
pos_errors.append(pos_err)
vel_errors.append(vel_err)
results[noise_arcsec] = {
'pos_mean': np.mean(pos_errors),
'pos_std': np.std(pos_errors),
'vel_mean': np.mean(vel_errors),
'vel_std': np.std(vel_errors),
}
print(f"ノイズ {noise_arcsec:5.1f} arcsec: "
f"位置誤差 = {np.mean(pos_errors):.2f} ± {np.std(pos_errors):.2f} km, "
f"速度誤差 = {np.mean(vel_errors):.2f} ± {np.std(vel_errors):.2f} m/s")
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
noise_vals = list(results.keys())
pos_means = [results[n]['pos_mean'] for n in noise_vals]
pos_stds = [results[n]['pos_std'] for n in noise_vals]
vel_means = [results[n]['vel_mean'] for n in noise_vals]
vel_stds = [results[n]['vel_std'] for n in noise_vals]
axes[0].errorbar(noise_vals, pos_means, yerr=pos_stds,
fmt='o-', linewidth=2, color='#2196F3',
capsize=5, markersize=8)
axes[0].set_xlabel('観測ノイズ [arcsec]', fontsize=12)
axes[0].set_ylabel('位置誤差 [km]', fontsize=12)
axes[0].set_title('観測ノイズ vs 位置推定誤差', fontsize=13)
axes[0].grid(True, alpha=0.3)
axes[1].errorbar(noise_vals, vel_means, yerr=vel_stds,
fmt='s-', linewidth=2, color='#F44336',
capsize=5, markersize=8)
axes[1].set_xlabel('観測ノイズ [arcsec]', fontsize=12)
axes[1].set_ylabel('速度誤差 [m/s]', fontsize=12)
axes[1].set_title('観測ノイズ vs 速度推定誤差', fontsize=13)
axes[1].grid(True, alpha=0.3)
plt.suptitle('ガウス法の観測ノイズに対する感度分析（モンテカルロ, N=200）',
fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
このモンテカルロシミュレーションの結果から、軌道決定の精度に関する重要な知見が得られます。
- 位置誤差と速度誤差は観測ノイズにほぼ線形に比例します。これは問題が局所的に線形であることの反映であり、最小二乗法による軌道改良が有効であることを示唆しています。
- 1 arcsecの観測精度で、位置誤差はkm程度、速度誤差は数m/s程度です。現代の地上レーダーは角度精度0.01°（36 arcsec）程度、光学望遠鏡は0.1〜1 arcsec程度の精度を持つため、初期軌道決定としてはこの精度で十分な場合が多いです。
- 誤差の標準偏差（エラーバー）も観測ノイズとともに増加しており、ノイズが大きいほど推定の不確実性も大きくなることがわかります。これが、多数の観測データを用いた軌道改良（最小二乗法）が重要な理由です。
軌道決定の実務的考慮事項
観測の幾何学的配置
軌道決定の精度は、観測データの量だけでなく幾何学的配置にも大きく依存します。3回の角度観測が全て同じ方向から行われた場合、軌道の面外方向の情報がほとんど得られず、精度が著しく低下します。理想的には、衛星の軌道面に対して異なる方向から観測するのが望ましいです。
GDOP（Geometric Dilution of Precision）の概念は、GPSでおなじみですが、軌道決定にも同様に適用されます。ヤコビアン行列 $\bm{H}$ の特異値が軌道決定の「感度方向」を示しており、小さい特異値に対応する方向は推定精度が悪くなります。
観測データの種類
実際の軌道決定では、角度観測だけでなくさまざまな種類の観測データが使われます。
| 観測種類 | 精度 | 特徴 |
|---|---|---|
| 角度（赤経・赤緯） | 0.1〜10 arcsec | 光学望遠鏡、パッシブ観測 |
| レンジ（距離） | 1〜10 m | レーダー、レーザー測距 |
| レンジレート（距離変化率） | 0.01〜1 mm/s | ドップラー計測 |
| GPS | 1〜10 m | 衛星搭載GPS受信機 |
これらの観測データを組み合わせて使うことで、軌道決定の精度と信頼性が向上します。
摂動の考慮
本記事ではケプラー軌道（二体問題）を前提としましたが、実際の高精度軌道決定では以下の摂動を考慮する必要があります。
- 地球の非球対称性（$J_2$ 以下の重力調和項）
- 第三天体の引力（月、太陽）
- 大気抵抗（LEO）
- 太陽輻射圧
- 固体潮汐・海洋潮汐
これらの摂動は軌道伝播（観測方程式の計算）に含める必要があり、計算コストが大幅に増加します。現代の軌道決定ソフトウェア（GMAT, Orekit, ODTK など）はこれらの摂動を精密にモデル化しています。
まとめ
本記事では、軌道決定の理論と実装について解説しました。
- 軌道決定は「観測から軌道を推定する」逆問題であり、初期軌道決定と軌道改良の2段階で行う
- ガウスの方法は3回の角度観測から軌道を決定する古典的手法であり、ラグランジュ係数の級数展開と反復法で解く
- 最小二乗法による軌道改良は、観測方程式を線形化してヤコビアン行列を計算し、正規方程式を解く
- 推定精度は共分散行列 $\bm{P}$ で評価でき、観測の幾何学的配置と観測精度に依存する
- 観測ノイズの影響はモンテカルロシミュレーションで定量的に評価できる
軌道決定は天体力学の古典的な問題でありながら、宇宙デブリの増加に伴い、近年ますます重要性を増しています。特に、低軌道の混雑化により、高精度かつ迅速な軌道決定の需要が高まっています。
次のステップとして、以下の記事も参考にしてください。