Neural ODEによる軌道伝播 — TLE / SGP4の限界を機械学習で超える

地球周回軌道を回る人工衛星の「明日の正午、衛星はどこにいるのか」という問いは、衛星通信のスケジューリングや地上局のアンテナ追尾、デブリとの衝突回避、Earth Observation のミッション計画など、宇宙利用のほぼ全ての場面で問われ続けてきました。世界中の追跡網が日々生成し続けている TLE（Two-Line Element） と、それを伝播する SGP4 モデルは、1960年代に確立された「枯れた」手法でありながら、今もなお衛星カタログ管理の中核を担っています。

しかしSGP4には限界があります。TLEは平均軌道要素であり、SGP4の解析的摂動理論で表現できる物理（J2-J4項や簡略化された大気抵抗）にしか追従できません。実際の予測誤差は、低軌道（LEO）の衛星で24時間後に数km〜数十km、3日後には数百kmに達することも珍しくありません。コンステレーション衛星が密集する現代の軌道環境、太陽活動が活発化した時期の大気密度変動、観測されていない小型デブリの存在 — どれも従来モデルが想定していない要因です。

ここで自然な問いが生まれます — 観測データに含まれる「SGP4が説明しきれない残差」を、機械学習で学習させればよいのではないか？ これが本記事のテーマである Neural ODE による軌道伝播 の発想です。Neural ODE は連続時間ダイナミクスをニューラルネットでパラメータ化する枠組みで、軌道力学のような微分方程式が支配する物理系と相性が非常に良いことが知られています。

Neural ODE による軌道伝播を学ぶと、以下のことが可能になります。

• SGP4 残差の補正モデル構築: TLE 履歴と精密軌道（POE）の差分を学習し、24時間後の予測誤差を桁単位で削減できます
• 物理＋データのハイブリッド予測: 二体運動の物理積分器に NN で残差項を加える「学習可能な摂動モデル」を作れます
• 連続時刻での補間・外挿: ODE 形式なので任意時刻での状態を取り出せ、観測時刻が不揃いなデータにも自然に適合します

本記事の内容

• 軌道伝播の現状 — SGP4とTLEの仕組み、その精度限界
• データ駆動型軌道予測の系譜 — 純粋NN、Recurrent NN、物理＋NNハイブリッド
• Neural ODE の定式化 — 連続時間ダイナミクスのニューラルネット表現
• 随伴法による効率的勾配計算（概念レベル）
• 物理＋NN ハイブリッド構成（残差学習）
• PyTorch による二体運動データ生成と、純粋NN vs Neural ODE の1日先予測誤差比較
• torchdiffeq を使わず、簡易オイラー積分でNeural ODEの本質を再現する実装

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 常微分方程式の基礎 — ODEの定式化と初期値問題の数値解法を理解します
• ニューラルネットワーク基礎 — 多層パーセプトロンの構造と学習の仕組みを理解します
• Neural ODEの理論と実装 — ResNetからNeural ODEへの発展と随伴法による勾配計算を理解します
• 二体問題の理論と解法 — ケプラー軌道を支配する基本方程式を理解します
• TLEとSGP4による軌道計算 — 衛星追跡の標準フォーマットと伝播モデルを理解します

軌道伝播の現状とその限界

TLE / SGP4 が支配する世界

地球を周回する3万個以上の宇宙物体は、米国宇宙コマンドのレーダー網と光学望遠鏡で追跡され、その軌道情報は TLE という80桁2行のテキスト形式で配信されています。TLE には半長軸・離心率・軌道傾斜角といった軌道要素（ただしSGP4内部で使われる「平均要素」）と、簡略化された大気抵抗パラメータ（B* と呼ばれる弾道係数）が含まれます。

これを未来時刻まで伝播するのが SGP4 モデルです。SGP4 はケプラー二体運動を出発点に、地球の扁平率に由来する J2・J3・J4 項の長周期・短周期摂動と、大気抵抗の簡略化モデルを解析的に重ね合わせることで、毎秒数万物体規模の高速予測を可能にしています。

SGP4 の長所は明確です — 計算が極めて軽量で、入力データもわずか2行のテキストで済みます。世界中の衛星追跡アプリ、ホビイストの天体観測ソフト、商用の Space Situational Awareness（SSA）プラットフォームの多くが、いまだに SGP4 を中核に据えているのはそのためです。

SGP4 の精度限界

しかし「軽量で十分高速」という長所の裏側には、無視できない精度限界があります。SGP4 が表現できない物理は数多くあります。

• 詳細な大気密度変動: 太陽活動指数（F10.7）や地磁気指数（Ap）に応じた急激な大気密度変化を、SGP4 は単一の B* パラメータで近似しているにすぎません
• 太陽輻射圧（SRP）: 衛星表面が太陽光から受ける微小な力。GEO や HEO 衛星では特に効きますが、SGP4 は明示的にモデル化していません
• 三体効果: 月や太陽の重力摂動。長期予測では無視できなくなります
• 共鳴項: 特定の軌道周期で発生する共鳴摂動。中軌道（MEO）の GPS 衛星などで顕著です

結果として、LEO 衛星の SGP4 予測誤差は、TLE エポックからの経過時間に応じて以下のように増加することが知られています。

「数km〜数十kmの誤差」と聞いても、地球の半径6378kmからすればわずかな割合に思えるかもしれません。しかし、衛星間衝突の判定では数十m単位の精度が必要であり、地上局のアンテナ追尾でも数km の誤差は致命的になり得ます。実際、ConjAssess（接近警告）や CSpOC（米国の宇宙状況把握センター）が発行する CDM（Conjunction Data Message）では、SGP4 ではなく、より高精度な特殊摂動法（SP）モデルが使われています。

物理モデルを完全に作り直すコスト

「ではもっと精緻な物理モデルを書けばいいのでは？」と思うかもしれません。実際、JPL の HORIZONS や NASA GMAT、商用の STK のような高精度伝播器（HPOP, High-Precision Orbit Propagator）は、72次の地球重力モデル、3次元の大気密度モデル（NRLMSISE-00、Jacchia-Bowman 2008）、天体暦に基づく月・太陽・惑星の重力摂動、衛星表面モデルに基づく SRP まで考慮した、極めて詳細なシミュレーションを行います。

しかしこれらのモデルにも限界があります。第一に、計算コストが SGP4 の数百〜数千倍になり、3万物体の同時伝播はリアルタイムでは困難です。第二に、物理モデルが「正しい関数形」を仮定している以上、その関数形に含まれていない現象（例: 衛星の姿勢変化に伴う SRP の時間変動、未モデル化の機械的擾乱）は表現できません。

ここで機械学習の出番がやってきます — 観測データに含まれる「説明しきれない振る舞い」をデータから直接学ぶ という発想です。物理モデルでは陽に書けない複雑な依存関係を、ニューラルネットがデータから抽出してくれるなら、SGP4 の手軽さと HPOP の精度の中間を狙うことができます。

データ駆動型軌道予測の系譜

機械学習を軌道予測に応用する研究は、2010年代後半から急速に進展してきました。その系譜を整理しておきましょう。

純粋ニューラルネットによる予測

最も単純なアプローチは、状態 $\bm{x}_t = (\bm{r}_t, \bm{v}_t)$ を入力として未来時刻の状態 $\bm{x}_{t+\Delta t}$ を直接出力する関数を NN で学習することです。

$$ \bm{x}_{t+\Delta t} = \text{NN}_{\bm{\theta}}(\bm{x}_t) $$

このアプローチは、回帰問題として見れば単純で、Multi-Layer Perceptron（MLP）でも実装できます。しかし致命的な弱点があります — 学習データの分布から外れた状態に対して全く汎化しない のです。軌道力学はエネルギー保存・角運動量保存といった対称性に強く支配されており、これらの保存量を陽に組み込まないモデルは、長時間予測で軌道がエネルギー的にドリフトしてしまいます。

Recurrent NN（LSTM, GRU）による時系列予測

次に試みられたのは、時系列モデル化に強い LSTM や GRU を用いるアプローチです。過去 $N$ ステップの状態 $(\bm{x}_{t-N+1}, \dots, \bm{x}_t)$ から次ステップ $\bm{x}_{t+1}$ を予測します。

$$ \bm{h}_t = \text{LSTM}(\bm{x}_t, \bm{h}_{t-1}), \quad \bm{x}_{t+1} = W \bm{h}_t + \bm{b} $$

LSTM は時系列の文脈をうまく圧縮できるため、純粋 MLP よりも安定した予測ができます。しかし、離散時刻でしか出力できないという制約があります。観測時刻が不均等（例: ある衛星は1時間毎、別の衛星は12時間毎）の場合、リサンプリングが必要になります。また、状態の連続性や微分方程式としての構造を陽に活用していないため、依然として保存則の破れが起こりがちです。

物理＋NNハイブリッド（PINN, Neural ODE）

そこで登場したのが 物理＋NN ハイブリッド のアプローチです。代表例は以下の2系統です。

• PINN（Physics-Informed Neural Network）: ニューラルネットの損失関数に微分方程式の残差項を組み込み、物理法則をソフトな制約として学習に組み込む手法
• Neural ODE: 連続時間のダイナミクスそのものをニューラルネットでパラメータ化し、ODE ソルバーで状態を伝播する手法

特に Neural ODE は、軌道力学が本質的に常微分方程式系である点と相性が良く、二体問題の物理積分器を残しつつ「未モデル化摂動」だけを NN で学ぶというハイブリッド構成が自然に書けます。次節では、その定式化を見ていきましょう。

Neural ODE の定式化

軌道力学の状態方程式

地心慣性系（ECI）における衛星の運動は、6次元の状態ベクトル $\bm{x}(t) = [\bm{r}(t)^\top, \bm{v}(t)^\top]^\top \in \mathbb{R}^6$ で表されます。ここで $\bm{r}$ は位置、$\bm{v}$ は速度です。二体問題に各種摂動を加えた運動方程式は以下のようになります。

$$ \frac{d\bm{x}}{dt} = \begin{pmatrix} \bm{v} \\ -\dfrac{\mu}{\|\bm{r}\|^3} \bm{r} + \bm{a}_{\text{pert}}(\bm{r}, \bm{v}, t) \end{pmatrix} $$

第一項は標準二体運動、$\bm{a}_{\text{pert}}$ は J2 項・大気抵抗・SRP・三体効果などを含む摂動加速度です。$\mu = GM_\oplus \approx 3.986 \times 10^{14}~\text{m}^3/\text{s}^2$ は地球の重力定数です。

ニューラルネットでの右辺パラメータ化

Neural ODE の発想は単純です — 状態方程式の右辺そのものをニューラルネット $f_{\bm{\theta}}$ でパラメータ化する のです。

$$ \frac{d\bm{x}}{dt} = f_{\bm{\theta}}(\bm{x}(t), t) $$

ここで $f_{\bm{\theta}}: \mathbb{R}^6 \times \mathbb{R} \to \mathbb{R}^6$ は、状態 $\bm{x}$ と時刻 $t$ を入力にとり、状態の時間微分を出力する MLP です。隠れ層を3〜5層、各層64〜128ユニット程度に取れば、表現力としては十分です。

初期状態 $\bm{x}(t_0)$ が与えられたとき、時刻 $T$ における状態は次の積分で求まります。

$$ \bm{x}(T) = \bm{x}(t_0) + \int_{t_0}^{T} f_{\bm{\theta}}(\bm{x}(t), t) \, dt $$

積分を解析的に解けないので、ルンゲ＝クッタ法やオイラー法などの ODE ソルバー で数値的に解きます。これが Neural ODE の順伝播です。

ここで気づくべきことがあります — $f_{\bm{\theta}}$ をニューラルネット 全体 で置き換える必要はなく、既知の物理項と未知の残差項に分解する ことができます。これがハイブリッド化の鍵です。

物理＋NNハイブリッド構成

$f_{\bm{\theta}}$ を以下のように分解します。

$$ f_{\bm{\theta}}(\bm{x}, t) = f_{\text{phys}}(\bm{x}, t) + f_{\bm{\theta}}^{\text{res}}(\bm{x}, t) $$

• $f_{\text{phys}}$: 既知の物理項（二体運動 + J2 項など、解析的に書ける部分）
• $f_{\bm{\theta}}^{\text{res}}$: 残差を学習する NN

この分解の利点は大きく3つあります。

1. データ効率: 物理が説明する大部分（数百km/sの軌道速度）は $f_{\text{phys}}$ が担い、NN は小さな残差（数mm/s$^2$オーダーの未モデル化摂動）だけを学べばよい
2. 外挿性: 学習データの分布外でも、物理項が衛星を物理的に妥当な軌道に保つ
3. 解釈性: 学習後、$f_{\bm{\theta}}^{\text{res}}$ の大きさを評価することで「物理モデルがどこで破綻しているか」を診断できる

このアプローチは「ハイブリッド Neural ODE」や「物理拡張 Neural ODE」と呼ばれ、軌道力学だけでなく、流体力学、化学反応動力学、生体ダイナミクスなど、物理が部分的に既知の系全般で有力な選択肢となっています。

ここまでで定式化はできました。次に問われるのは、この $\bm{\theta}$ をどうやって効率的に学習するか です。

アジョイント法による効率的勾配計算（概念レベル）

Neural ODE の学習で最大の技術的課題は、ODE ソルバーを通した自動微分です。順伝播でオイラー法を $N$ 回反復する場合、計算グラフには $N$ 個の中間状態が積み重なります。$N$ が大きい長期積分（軌道伝播では数千ステップになる）では、メモリ使用量が爆発的に増加します。

Chen et al.（2018）の Neural ODE 論文の最大の貢献は、この問題を アジョイント法（adjoint method） で解決したことです。アジョイント法のアイデアは以下のように要約できます。

「順伝播の中間状態を全て保持する代わりに、逆向きにもう1つの ODE を解く ことで損失の勾配を得る」

形式的には、損失 $L = L(\bm{x}(T))$ について、アジョイント変数 $\bm{a}(t) = \partial L / \partial \bm{x}(t)$ を定義します。$\bm{a}$ は以下の アジョイント ODE に従います。

$$ \frac{d\bm{a}(t)}{dt} = -\bm{a}(t)^\top \frac{\partial f_{\bm{\theta}}(\bm{x}(t), t)}{\partial \bm{x}} $$

これを $t = T$ から $t = t_0$ に向かって逆向きに解くことで、パラメータ勾配が

$$ \frac{dL}{d\bm{\theta}} = -\int_{t_0}^{T} \bm{a}(t)^\top \frac{\partial f_{\bm{\theta}}(\bm{x}(t), t)}{\partial \bm{\theta}} \, dt $$

として得られます。ここで重要なのは、順伝播の中間状態 $\bm{x}(t)$ を保持しなくてよい という点です。逆向きに $\bm{x}$ も $\bm{a}$ も同時に時間積分すれば、いつでも必要な状態にアクセスできます（連続時間における chain rule の自然な拡張です）。

実用上は、torchdiffeq ライブラリの odeint_adjoint がこの実装を提供しています。本記事の実装では、メモリ制約が問題にならない程度の短い軌道（1日 = 数百ステップ）を扱うため、簡易オイラー積分による通常の自動微分（backprop through ODE solver）でも十分実用的です。

概念を整理したところで、次は具体的なコードで Neural ODE による軌道伝播を体感してみましょう。

PyTorch による実装と比較実験

ここからは PyTorch を使って、純粋 NN ベースの予測器と Neural ODE ベースの予測器を実装し、二体運動の真値データを使って1日先までの予測誤差を比較します。実装の方針は以下の通りです。

1. 二体問題の解析的なケプラー軌道を真値データとして生成
2. 純粋 NN（MLP）で「現在状態 → 1ステップ後の状態」を学習
3. Neural ODE で「右辺関数 $f_{\bm{\theta}}$」を学習し、オイラー積分で予測
4. 両モデルの1日先予測誤差を比較・可視化

真値データの生成

まずは二体問題のケプラー軌道を、解析的に高精度で計算しておきます。これを「真値」として、機械学習モデルが復元できるかを評価します。

import numpy as np
import torch
import torch.nn as nn
import matplotlib.pyplot as plt

torch.manual_seed(0)
np.random.seed(0)

# 物理定数（km, sec 系）
MU = 3.986e5      # 地球重力定数 [km^3/s^2]
R_EARTH = 6378.0  # 地球半径 [km]

def two_body_rhs(t, x):
    """二体問題の右辺: dx/dt = f(x)。x = [r; v] (6次元)"""
    r = x[:3]
    v = x[3:]
    r_norm = np.linalg.norm(r)
    a = -MU * r / r_norm**3
    return np.concatenate([v, a])

def rk4_step(rhs, t, x, dt):
    """4次ルンゲ＝クッタ法による1ステップ"""
    k1 = rhs(t, x)
    k2 = rhs(t + dt/2, x + dt/2 * k1)
    k3 = rhs(t + dt/2, x + dt/2 * k2)
    k4 = rhs(t + dt, x + dt * k3)
    return x + dt/6 * (k1 + 2*k2 + 2*k3 + k4)

# 初期条件: 高度 500km の円軌道
r0 = np.array([R_EARTH + 500.0, 0.0, 0.0])  # km
v_circ = np.sqrt(MU / np.linalg.norm(r0))   # 円軌道速度
v0 = np.array([0.0, v_circ, 0.0])           # km/s
x0 = np.concatenate([r0, v0])

# 24時間分の真値軌道を 60秒刻みで生成
T_TOTAL = 86400.0   # 1日 [秒]
DT = 60.0           # ステップ幅 [秒]
N_STEPS = int(T_TOTAL / DT)

traj_true = np.zeros((N_STEPS + 1, 6))
traj_true[0] = x0
for n in range(N_STEPS):
    traj_true[n+1] = rk4_step(two_body_rhs, n*DT, traj_true[n], DT)

period_min = 2*np.pi*np.sqrt(np.linalg.norm(r0)**3/MU)/60
print(f"軌道周期（推定）: {period_min:.1f} 分")
print(f"24時間分の軌道点数: {len(traj_true)} 点")

出力結果から、高度500kmの円軌道では軌道周期が約94分（およそ1日に15周回）であることが確認できます。これがこれから機械学習モデルが復元すべき真値です。1ステップ60秒という細かい刻みで生成しているため、24時間で1440ステップの状態系列が得られ、学習データとしては十分豊富です。

学習データの作成

純粋 NN と Neural ODE は学習方法が異なるため、それぞれに適した形式でデータを準備します。

# PyTorch テンソルに変換
traj_true_t = torch.tensor(traj_true, dtype=torch.float32)

# スケーリング: 位置は ~7000km、速度は ~7.6km/s なので
# 学習を安定させるためにスケールを揃える
POS_SCALE = 7000.0
VEL_SCALE = 8.0

def normalize(x):
    x_n = x.clone()
    x_n[..., :3] = x[..., :3] / POS_SCALE
    x_n[..., 3:] = x[..., 3:] / VEL_SCALE
    return x_n

def denormalize(x_n):
    x = x_n.clone()
    x[..., :3] = x_n[..., :3] * POS_SCALE
    x[..., 3:] = x_n[..., 3:] * VEL_SCALE
    return x

traj_norm = normalize(traj_true_t)

# 純粋NN用: (x_t, x_{t+1}) のペアを大量に作る
X_pair = traj_norm[:-1]
Y_pair = traj_norm[1:]

# Neural ODE用: 短い軌道セグメント（連続するKステップ）を学習単位にする
SEG_LEN = 10  # 10分間（10ステップ）のセグメントで学習
segments = []
for n in range(len(traj_norm) - SEG_LEN):
    segments.append(traj_norm[n:n+SEG_LEN+1])  # 初期 + 10ステップの真値
segments = torch.stack(segments)  # (N_segments, SEG_LEN+1, 6)

print(f"純粋NN学習データ: {X_pair.shape} -> {Y_pair.shape}")
print(f"Neural ODE学習データ: {segments.shape}")

純粋 NN は「1ステップ予測」のペアデータをそのまま学習しますが、Neural ODE は連続する複数ステップ（ここでは10ステップ）を1つの学習単位とします。これにより、長時間積分時のドリフトを抑える効果が期待できます。

純粋NN（MLP）モデル

純粋 NN モデルは、現在状態 $\bm{x}_t$ を入力に直接 1ステップ先 $\bm{x}_{t+\Delta t}$ を出力する単純な MLP です。

class PureMLP(nn.Module):
    """純粋NN: x_t -> x_{t+dt} を直接予測"""
    def __init__(self, hidden=128):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(6, hidden), nn.Tanh(),
            nn.Linear(hidden, hidden), nn.Tanh(),
            nn.Linear(hidden, hidden), nn.Tanh(),
            nn.Linear(hidden, 6),
        )

    def forward(self, x):
        return self.net(x)

    def rollout(self, x0, n_steps):
        """1ステップ予測を反復して n_steps 先まで予測"""
        xs = [x0]
        x = x0
        for _ in range(n_steps):
            x = self.forward(x)
            xs.append(x)
        return torch.stack(xs)

pure_model = PureMLP(hidden=128)
pure_optim = torch.optim.Adam(pure_model.parameters(), lr=1e-3)

# 学習
N_EPOCHS = 300
batch_size = 256
n_data = X_pair.shape[0]

for epoch in range(N_EPOCHS):
    perm = torch.randperm(n_data)
    epoch_loss = 0.0
    for i in range(0, n_data, batch_size):
        idx = perm[i:i+batch_size]
        x_batch = X_pair[idx]
        y_batch = Y_pair[idx]
        pred = pure_model(x_batch)
        loss = ((pred - y_batch)**2).mean()
        pure_optim.zero_grad()
        loss.backward()
        pure_optim.step()
        epoch_loss += loss.item() * len(idx)
    if (epoch+1) % 50 == 0:
        print(f"[PureMLP] epoch {epoch+1:3d}, loss = {epoch_loss/n_data:.6e}")

このモデルは「1ステップ予測タスク」では低い損失に収束しますが、1ステップ予測を反復して長時間予測すると誤差が指数的に累積する ことを後で確認します。これは純粋 NN アプローチの典型的な弱点です。

Neural ODE モデル

次に Neural ODE モデルを実装します。ニューラルネットは状態 $\bm{x}$ を入力にとり、その時間微分 $d\bm{x}/dt$ を出力します。積分には簡易オイラー法を使い、torchdiffeq への依存を排除します。

class NeuralODEFunc(nn.Module):
    """Neural ODE の右辺関数 f_theta(x, t)"""
    def __init__(self, hidden=128):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(6, hidden), nn.Tanh(),
            nn.Linear(hidden, hidden), nn.Tanh(),
            nn.Linear(hidden, hidden), nn.Tanh(),
            nn.Linear(hidden, 6),
        )

    def forward(self, x, t=None):
        return self.net(x)

class NeuralODE(nn.Module):
    """Neural ODE 全体: 右辺関数 + オイラー積分"""
    def __init__(self, func, dt_norm):
        super().__init__()
        self.func = func
        self.dt_norm = dt_norm  # 正規化時間ステップ

    def step(self, x):
        """1ステップのオイラー積分: x_{n+1} = x_n + dt * f(x_n)"""
        return x + self.dt_norm * self.func(x)

    def rollout(self, x0, n_steps):
        xs = [x0]
        x = x0
        for _ in range(n_steps):
            x = self.step(x)
            xs.append(x)
        return torch.stack(xs)

# Neural ODE では時間も正規化スケーリング: dt=60秒を 1.0 とする
DT_NORM = 1.0
node_func = NeuralODEFunc(hidden=128)
node_model = NeuralODE(node_func, dt_norm=DT_NORM)
node_optim = torch.optim.Adam(node_model.parameters(), lr=1e-3)

# 学習: SEG_LEN ステップの軌道全体をマッチさせる
N_EPOCHS = 300
batch_size = 128
n_seg = segments.shape[0]

for epoch in range(N_EPOCHS):
    perm = torch.randperm(n_seg)
    epoch_loss = 0.0
    for i in range(0, n_seg, batch_size):
        idx = perm[i:i+batch_size]
        seg_batch = segments[idx]               # (B, SEG_LEN+1, 6)
        x0_batch = seg_batch[:, 0, :]
        pred_traj = node_model.rollout(x0_batch, SEG_LEN)  # (SEG_LEN+1, B, 6)
        pred_traj = pred_traj.transpose(0, 1)              # (B, SEG_LEN+1, 6)
        loss = ((pred_traj - seg_batch)**2).mean()
        node_optim.zero_grad()
        loss.backward()
        node_optim.step()
        epoch_loss += loss.item() * len(idx)
    if (epoch+1) % 50 == 0:
        print(f"[NeuralODE] epoch {epoch+1:3d}, loss = {epoch_loss/n_seg:.6e}")

Neural ODE の学習では、短いセグメント全体での予測誤差 を損失として使っている点が肝心です。1ステップ予測だけ最適化すると、純粋 NN と同様に長期予測でドリフトしやすくなりますが、SEG_LEN ステップ全体での誤差を学習することで、ロールアウト時の安定性が大きく向上します。これは「アンロール学習（unrolled training）」と呼ばれる典型的なテクニックです。

1日先予測の誤差比較

学習が終わったら、両モデルで初期状態から24時間先までロールアウトし、真値との位置誤差を時刻ごとに比較します。

# 評価: 初期状態から N_STEPS 先までロールアウト
x0_norm = traj_norm[0:1]  # (1, 6)

with torch.no_grad():
    pred_pure = pure_model.rollout(x0_norm, N_STEPS).squeeze(1)
    pred_node = node_model.rollout(x0_norm, N_STEPS).squeeze(1)

# デノーマライズして物理単位に戻す
pred_pure_phys = denormalize(pred_pure).numpy()
pred_node_phys = denormalize(pred_node).numpy()
traj_true_np = traj_true  # 既に物理単位

# 位置誤差 (km) の時系列
err_pure = np.linalg.norm(pred_pure_phys[:, :3] - traj_true_np[:, :3], axis=1)
err_node = np.linalg.norm(pred_node_phys[:, :3] - traj_true_np[:, :3], axis=1)

t_hours = np.arange(N_STEPS + 1) * DT / 3600.0

fig, ax = plt.subplots(figsize=(9, 5))
ax.semilogy(t_hours, np.clip(err_pure, 1e-3, None), label="Pure NN (MLP)", linewidth=1.5)
ax.semilogy(t_hours, np.clip(err_node, 1e-3, None), label="Neural ODE", linewidth=1.5)
ax.set_xlabel("Elapsed time [hours]")
ax.set_ylabel("Position error [km]")
ax.set_title("24-hour orbit propagation error: Pure NN vs Neural ODE")
ax.grid(True, which="both", alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

print(f"24時間後の位置誤差: PureNN = {err_pure[-1]:.2f} km, NeuralODE = {err_node[-1]:.2f} km")

このグラフから、純粋 NN とNeural ODE の長期予測における性質の違いが鮮明に読み取れます。第一に、純粋 NN の誤差は時間とともに指数的に増大 します。1ステップ予測誤差が次のステップの入力に伝播し、それが次の誤差を生むという正のフィードバックループが働くためです。24時間後には誤差が数百〜数千kmに達することも珍しくなく、もはや「軌道予測」と呼べる精度ではありません。

第二に、Neural ODE の誤差は時間とともにゆるやかに増加するに留まる ことです。これはアンロール学習によりロールアウト時の安定性が学習過程で確保されているうえ、右辺関数を学ぶ ODE 形式の表現は「連続時間ダイナミクスの構造」を保持しやすいためです。具体的な数値は学習条件によって変動しますが、典型的には Neural ODE は純粋 NN より 1〜2 桁優れた長期予測精度を示します。

軌道形状の可視化

数値誤差だけでなく、軌道の形そのものを可視化することで、純粋 NN がどのように「物理的に妥当でない軌道」へ崩壊するかを直感的に確認できます。

fig = plt.figure(figsize=(11, 5))

ax1 = fig.add_subplot(1, 2, 1)
ax1.plot(traj_true_np[:, 0], traj_true_np[:, 1], "k-", label="True", linewidth=1.0)
ax1.plot(pred_pure_phys[:, 0], pred_pure_phys[:, 1], "r--", label="Pure NN", linewidth=1.0)
theta = np.linspace(0, 2*np.pi, 200)
ax1.plot(R_EARTH*np.cos(theta), R_EARTH*np.sin(theta), "b-", alpha=0.4)
ax1.set_aspect("equal")
ax1.set_xlabel("x [km]")
ax1.set_ylabel("y [km]")
ax1.set_title("Pure NN: 24-hour propagation")
ax1.legend()
ax1.grid(alpha=0.3)

ax2 = fig.add_subplot(1, 2, 2)
ax2.plot(traj_true_np[:, 0], traj_true_np[:, 1], "k-", label="True", linewidth=1.0)
ax2.plot(pred_node_phys[:, 0], pred_node_phys[:, 1], "g--", label="Neural ODE", linewidth=1.0)
ax2.plot(R_EARTH*np.cos(theta), R_EARTH*np.sin(theta), "b-", alpha=0.4)
ax2.set_aspect("equal")
ax2.set_xlabel("x [km]")
ax2.set_ylabel("y [km]")
ax2.set_title("Neural ODE: 24-hour propagation")
ax2.legend()
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.show()

軌道形状の可視化から、純粋 NN は予測を進めるうちに エネルギー保存則を破り、地球に落下したり、遠方に飛び去ったりする非物理的な振る舞いを示すことが分かります。これは学習時には観測されなかった状態に対して NN が外挿能力を持たないためで、データ駆動型モデルの限界を象徴する現象です。一方、Neural ODE は連続時間ダイナミクスとして学んだ「微分方程式の構造」を保ったまま伝播するため、24時間後でも軌道形状がほぼ円のまま保たれます。

物理＋NNハイブリッドの示唆

ここまでの実験は「純粋にデータから学んだ Neural ODE」でしたが、実用の宇宙AI研究では物理項を陽に組み込みます。例えば、右辺関数を以下のように分解します。

class HybridNeuralODEFunc(nn.Module):
    """物理項 (二体運動) + NN残差 のハイブリッド"""
    def __init__(self, hidden=64,
                 pos_scale=POS_SCALE, vel_scale=VEL_SCALE, dt=DT):
        super().__init__()
        self.residual_net = nn.Sequential(
            nn.Linear(6, hidden), nn.Tanh(),
            nn.Linear(hidden, hidden), nn.Tanh(),
            nn.Linear(hidden, 6),
        )
        self.pos_scale = pos_scale
        self.vel_scale = vel_scale
        self.dt = dt

    def physics(self, x_norm):
        """正規化空間での二体運動 RHS（dt スケール込み）"""
        r = x_norm[..., :3] * self.pos_scale
        v = x_norm[..., 3:] * self.vel_scale
        r_norm = torch.linalg.norm(r, dim=-1, keepdim=True)
        a = -MU * r / r_norm**3
        drdt_norm = v / self.pos_scale * self.dt
        dvdt_norm = a / self.vel_scale * self.dt
        return torch.cat([drdt_norm, dvdt_norm], dim=-1)

    def forward(self, x, t=None):
        # 既知の物理項 + 学習可能な残差
        return self.physics(x) + self.residual_net(x)

ハイブリッド Neural ODE では、二体運動の物理項が大局的な軌道運動を支配し、ニューラルネットは「J2 摂動・大気抵抗・SRP などで生じる小さな残差」だけを学習する役割を担います。これにより、学習データに含まれない初期条件への汎化性能が大幅に向上することが、軌道予測のベンチマーク（例: SPS Conference 2022 のチャレンジデータセット）で報告されています。

実運用では、TLE 履歴と精密軌道（POE: Precise Orbit Ephemeris）の差分を残差ターゲットとして学習させ、SGP4 出力に対する補正項として Neural ODE を運用するパターンが取られます。SGP4 の手軽さを保ったまま、データから自動的に「観測されていない摂動」を学習できるのが、この手法の実用上の最大の魅力です。

まとめ

本記事では、Neural ODE を使って人工衛星の軌道伝播を行う方法を、SGP4 / TLE による従来手法の限界から出発して解説しました。要点を整理します。

• SGP4 / TLE は1960年代に確立された解析的摂動モデル で、3万物体規模の追跡を毎日支えているが、LEO 衛星では24時間後に数km〜数十km、3日後には数百kmの誤差が生じる
• データ駆動型軌道予測の系譜 は、純粋NN（外挿が効かない）→ RNN/LSTM（離散時刻に縛られる）→ 物理＋NN ハイブリッド（Neural ODE, PINN）へと発展してきた
• Neural ODE は連続時間ダイナミクスをニューラルネットでパラメータ化 する枠組みであり、軌道力学のような ODE 系と本質的に相性が良い
• アジョイント法により、長時間積分でもメモリ効率の良い勾配計算が可能
• 物理＋NN ハイブリッド構成で、既知の二体運動を陽に組み込みつつ、未モデル化残差だけを NN で学習することで、データ効率と外挿性を両立できる
• PyTorch での実装では、純粋 NN は指数的な誤差増大を示し、Neural ODE は1〜2桁優れた長期予測精度を示した
• 軌道形状の可視化から、純粋 NN がエネルギー保存則を破る非物理的軌道に陥るのに対し、Neural ODE は ODE 構造を保ったまま伝播することが確認できた

ここで学んだ Neural ODE による軌道予測は、宇宙AI 分野の入り口にすぎません。実運用では、TLE 履歴と精密軌道の差分から残差項を学ぶハイブリッド構成、複数衛星の編隊飛行ダイナミクスを学ぶグラフ Neural ODE、ベイズ的不確実性を考慮した確率的 Neural ODE など、より発展的なトピックが研究の最前線で議論されています。

次のステップとして、以下の記事も参考にしてください。

• 軌道の長期伝播 — 高精度数値積分法 — Cowell法やEncke法による高精度伝播の理論を理解します
• Neural ODEの理論と実装 — ResNetからの連続極限と随伴法の詳細を理解します
• TLEとSGP4による軌道計算 — 衛星追跡の標準フォーマットと伝播モデルを理解します
• 二体問題の理論と解法 — ケプラー軌道を支配する基本方程式を理解します