軌道摂動論と追跡予測 — 二体問題から現実の軌道力学へ

この回で学ぶこと

前回、カルマンフィルタによる逐次的な軌道決定を学びました。あの議論の中で、状態遷移方程式

$$\mathbf{x}_k = F_{k-1}\, \mathbf{x}_{k-1} + \mathbf{w}_{k-1}$$

の $F_{k-1}$ について、「実際の軌道決定では二体問題や多体問題の運動方程式を数値積分することで得られる」と述べつつも、具体的な数式は等速直線運動 $F = \begin{bmatrix}1 & T\\ 0 & 1\end{bmatrix}$ という単純化した例で済ませ、$F$ の中身そのものには踏み込みませんでした。カルマンフィルタのアルゴリズム自体は、$F$ がどれほど精密であっても粗くても同じ5本の式で動きます。しかし、$F$ の精度が低ければ、フィルタがどれほど数学的に「最適」であっても、予測 $\hat{\mathbf{x}}_k^{-} = F_{k-1}\hat{\mathbf{x}}_{k-1}$ そのものが真の軌道からずれてしまうことは避けられません。

この回では、その $F$ の中身、すなわち探査機や地球周回衛星が実際にどのような力を受けて運動しているかを見ていきます。出発点は最も単純な二体問題(two-body problem)——中心天体と探査機だけがあり、両者の間に働く逆二乗則の重力だけを考える理想化されたケプラー運動——です。しかし現実には、地球の形が完全な球ではないこと、大気があること、太陽光が圧力を及ぼすこと、月や太陽が「第三体」として余分な重力を及ぼすことなど、二体問題からのさまざまなズレ=**摂動(perturbation)**が存在します。これらを $F$ にどこまで正確に組み込むかが、軌道予測の精度を、ひいては「次にいつどこにアンテナを向けるべきか」という追跡運用そのものの精度を左右します。

直感的な導入 — なぜ「ズレ」を無視できないのか

まず二体問題そのものを式で確認しておきましょう。中心天体(地球)の重力定数を $GM$、探査機の位置ベクトルを $\vec r$ とすると、二体問題の運動方程式は

$$\ddot{\vec r} = -\frac{GM}{r^3}\,\vec r, \qquad r \equiv |\vec r|$$

というきわめてシンプルな形をしています。この式の解は厳密に閉じた形(楕円・放物線・双曲線のケプラー軌道)で書け、軌道要素(半長径・離心率・軌道傾斜角・昇交点赤経・近地点引数・平均近点角の6つ)は時間が経っても変化しません。もし現実の探査機がこの式に厳密に従うなら、一度軌道要素を決定してしまえば、何年先の位置でも解析的に計算でき、追跡予測は非常に簡単な話で終わります。

しかし実際にはそうなりません。地球は完全な球ではなく赤道方向に膨らんだ回転楕円体に近い形をしており、その重力ポテンシャルは中心からの距離だけで決まらず緯度にも依存します。低軌道(LEO)では希薄とはいえ大気があり、探査機の運動エネルギーをじわじわと奪います。太陽光は光子の運動量として実際に力を及ぼしますし、月や太陽は「中心天体」ではないものの、無視できない重力を及ぼします。これらをまとめて摂動加速度 $\vec a_{pert}$ として、運動方程式を

$$\ddot{\vec r} = -\frac{GM}{r^3}\,\vec r + \vec a_{pert}(\vec r, \vec v, t)$$

と書き直したものが、現実の軌道力学モデルの出発点です。$\vec a_{pert}$ を無視して二体問題だけで軌道を予測すると、時間が経つにつれて予測位置は真の位置からどんどん外れていきます。前回の言葉で言えば、これは「モデル化しきれていない加速度成分」がプロセスノイズ $\mathbf{w}_{k-1}$ として押し込まれてしまうことを意味し、$Q_{k-1}$ を大きく見積もらざるを得なくなります。$Q$ が大きいと、予測共分散 $P_k^{-} = F_{k-1}P_{k-1}F_{k-1}^\top + Q_{k-1}$ の増大が速くなり、次の観測が来るまでの「今どこにいるか」の確信度が急速に失われます。つまり、摂動をどれだけ精密に $F$(あるいは非線形なら $f$)に組み込めるかは、カルマンフィルタが実用上どれだけ「賢く」振る舞えるかに直結するわけです。

摂動加速度の分類とオーダー比較

摂動加速度 $\vec a_{pert}$ を代表的な4つの要因に分解しましょう。

$$\vec a_{pert} = \vec a_{J_2\text{以下の非球対称重力}} + \vec a_{\text{大気抵抗}} + \vec a_{\text{太陽輻射圧}} + \vec a_{\text{第三体}}$$

非球対称重力(J2項)

以前、重力ポテンシャルを球面調和関数で展開する方法を学びました。地球のように自転による扁平がある天体では、展開の中で最も支配的な項は $n=2$ の帯調和項 $J_2$ で、その摂動加速度の大きさはおおむね

$$a_{J_2} \sim \frac{3}{2}\, J_2 \left(\frac{R_E}{r}\right)^2 \frac{GM}{r^2}$$

というオーダーでした(導出は前回のレッスンを参照してください)。地球の場合、$J_2 = 1.08263\times10^{-3}$、赤道半径 $R_E = 6378.137\ \text{km}$ です。$J_2$ が他の帯調和項($J_3, J_4,\ldots$、代表的には $10^{-6}$ オーダー)より3桁近く大きいため、地球周回衛星の軌道決定において、非球対称重力の摂動は事実上 $J_2$ 項が支配していると言ってよいほどです。

大気抵抗

低軌道(LEO)では、希薄な大気による抗力が無視できません。抗力は相対速度ベクトル $\vec v_{rel}$(大気の共回転を考慮した、探査機から見た大気の相対速度)と逆向きに働き、

$$\vec a_{\text{drag}} = -\frac{1}{2}\,\rho(h)\, v_{rel}\, \left(\frac{C_D A}{m}\right) \vec v_{rel}$$

という形で書けます。$\rho(h)$ は高度 $h$ における大気密度、$C_D$ は抗力係数(通常 2.0〜2.2程度)、$A$ は探査機の断面積、$m$ は質量で、$C_D A/m$ の逆数は**弾道係数(ballistic coefficient)**と呼ばれます。大気密度は高度とともに指数関数的に減少するモデル

$$\rho(h) \approx \rho_0 \exp\!\left(-\frac{h-h_0}{H}\right)$$

でおおまかに近似できます($H$ はスケールハイト、地球低軌道高度でおよそ50〜60 kmですが、太陽活動によって大きく変動します)。この指数的な減衰のため、大気抵抗の重要性は高度に極めて敏感です。国際宇宙ステーション(ISS)が飛ぶ高度400 km付近では無視できない摂動である一方、静止軌道(GEO、高度約35786 km)では実質的にゼロとみなせます。

太陽輻射圧(SRP)

太陽光子が探査機表面に当たって運動量を伝える力で、

$$\vec a_{\text{SRP}} = \frac{P_{sr}}{c}\, C_R \left(\frac{A}{m}\right) \hat u_{\odot}$$

と書けます。$P_{sr}/c \approx 4.56\times10^{-6}\ \text{N/m}^2$ は1天文単位における太陽輻射圧(光の運動量流束を光速 $c$ で割ったもの)、$C_R$ は反射特性を表す係数(完全吸収なら1、完全反射なら2、実際の衛星表面ではおよそ1.2〜1.5)、$\hat u_{\odot}$ は太陽から探査機への方向単位ベクトルです。太陽輻射圧は大気抵抗と違って高度に依存しないため、大気のないGEOや深宇宙巡航フェーズでは相対的に重要性が増します。

第三体重力(月・太陽)

地球以外の天体(月、太陽)も探査機に直接重力を及ぼします。ただしこれは天体そのものの重力ではなく、地球中心と探査機とで受ける第三体重力の「差」(潮汐的な効果)が効いてきます。

$$\vec a_{3\text{rd}} = GM_{3\text{rd}} \left[ \frac{\vec r_{3\text{rd}} - \vec r}{|\vec r_{3\text{rd}} - \vec r|^3} - \frac{\vec r_{3\text{rd}}}{|\vec r_{3\text{rd}}|^3} \right]$$

探査機と地球中心が第三体から見てほぼ同じ距離にある場合(地球中心からの距離 $r$ が第三体までの距離 $d$ に比べて十分小さい場合)、この差はテイラー展開により

$$a_{3\text{rd}} \sim \frac{2\, GM_{3\text{rd}}\, r}{d^3}$$

というオーダーで近似できます。$d$ の3乗で急減するため、地球に近い低軌道ほど第三体の影響は相対的に小さく、逆に地球から遠いGEOや、まして地球を離れた深宇宙軌道では急速に重要性が増します。

オーダー比較

以上4つの摂動加速度を、二体問題の中心重力加速度 $a_{2body} = GM/r^2$ に対する比(相対的な大きさ)として、代表的な2つの軌道高度で比較してみます(いずれも典型的な衛星パラメータを仮定した概算値で、実際のミッションでは機体形状や太陽活動によって数倍〜1桁程度変動します)。

摂動要因	LEO (高度400 km)	GEO (高度約35786 km)
$J_2$(地球扁平)	$\sim 1.4\times10^{-3}$	$\sim 3.7\times10^{-5}$
大気抵抗	$\sim 7\times10^{-7}$(太陽活動で大きく変動)	ほぼゼロ(無視可)
太陽輻射圧	$\sim 1.5\times10^{-8}$	$\sim 5.7\times10^{-7}$
第三体(月・太陽)	$\sim 1.3\times10^{-7}$	$\sim 3.2\times10^{-5}$

この表から読み取れる実務上の教訓は明快です。LEOでは $J_2$ が圧倒的に支配的な摂動であり、次いで大気抵抗が(特に太陽活動が高い時期には)無視できない。一方GEOでは大気抵抗はほぼ無関係になり、代わりに $J_2$・第三体・太陽輻射圧が同程度のオーダーで拮抗するため、GEO衛星のステーションキーピング(南北・東西方向の軌道保持)は主に月・太陽の重力と太陽輻射圧への対策として設計されています(実際のGEO衛星の年間デルタV予算は、東西方向で数 m/s、南北方向で40〜50 m/s程度が典型的とされ、南北方向の大部分は月・太陽の重力による軌道傾斜角の変化を打ち消すために使われています)。

J2摂動による軌道要素の永年変化

$J_2$ が最も支配的な摂動であることが分かったところで、それが軌道要素にどのような影響を及ぼすかを見てみましょう。二体問題では軌道要素は時間に対して一定でしたが、$J_2$ 摂動を受けると、軌道面全体がゆっくりと**歳差運動(プレセッション)します。摂動論(平均要素法、あるいはLagrangeの惑星方程式)を使ってこれを解析的に導出すると、次のような永年変化(secular variation、軌道1周期にわたって平均すると消えずに積み残る変化)**の結果式が得られます(導出はここでは扱わず、結果とその意味に焦点を当てます)。

昇交点赤経(RAAN)の歳差運動:

$$\dot\Omega \approx -\frac{3}{2}\, n\, J_2 \left(\frac{R_E}{p}\right)^2 \cos i$$

近地点引数の変化:

$$\dot\omega \approx \frac{3}{4}\, n\, J_2 \left(\frac{R_E}{p}\right)^2 \big(5\cos^2 i - 1\big)$$

ここで $n = \sqrt{GM/a^3}$ は平均運動(軌道1周にかかる角速度)、$a$ は半長径、$p = a(1-e^2)$ は半直弦(半通径)、$i$ は軌道傾斜角です。半長径 $a$、離心率 $e$、軌道傾斜角 $i$ 自体には(1周期平均すれば)永年変化はありませんが、昇交点赤経 $\Omega$ と近地点引数 $\omega$ は時間とともに一方向にゆっくり回転し続けるという結果になっています。

この式が意味するところは、単なる数学的な好奇心ではなく、追跡運用にとって死活的に重要です。もし軌道決定・予測ソフトウェアが $J_2$ 摂動を軌道力学モデル $F$ に組み込んでいなければ、$\Omega$ や $\omega$ が実際には毎日少しずつ回転しているにもかかわらず、予測軌道はいつまでも初期時刻の軌道面に固定されたままになります。数日〜数週間先の追跡パス予測をこの誤った軌道面で計算すれば、探査機が実際に現れる方向とアンテナのビームが向く方向とのあいだに、無視できないポインティング誤差が蓄積することになります。

$\dot\Omega$ の式は、実は工学的に積極的に利用されることでも有名です。$\dot\Omega$ を、地球が太陽の周りを1年で1周する角速度(約 $0.9856\ \text{deg/day}$)にちょうど一致させるように高度と軌道傾斜角を選ぶと、軌道面が常に太陽に対して同じ向きを保つ**太陽同期軌道(sun-synchronous orbit, SSO)**が実現できます。地球観測衛星の多くがこの軌道を採用しているのは、$\cos i$ の符号を選ぶことで $\dot\Omega$ の向きを制御できるという、$J_2$ 摂動の永年項をむしろ設計上の自由度として使っているからです(演習問題2で実際に確認します)。

摂動モデルの精度と追跡予測のトレードオフ

ここまでの内容を、カルマンフィルタの枠組みに戻して整理しましょう。予測ステップの共分散伝播式

$$P_k^{-} = F_{k-1}\, P_{k-1}\, F_{k-1}^\top + Q_{k-1}$$

において、$F_{k-1}$(あるいはEKFでの非線形伝播関数 $f$)に $J_2$、大気抵抗、太陽輻射圧、第三体重力をどこまで精密に組み込むかは、2つの経路で追跡予測の質に影響します。

1. バイアス(モデル化誤差そのもの): $F$ に含まれない力が実際に働いていれば、予測値 $\hat{\mathbf{x}}_k^{-} = F_{k-1}\hat{\mathbf{x}}_{k-1}$ そのものが系統的に真の状態からずれます。これはカルマンフィルタが仮定する「プロセスノイズは平均ゼロ」という前提を破ることになり、フィルタが自信満々に間違った位置を報告する(共分散 $P$ は小さいのに実際の誤差は大きい、いわゆるフィルタの発散や過信)原因になります。
2. 不確かさの増大速度: $F$ に含めきれない力を仕方なくプロセスノイズ $Q$ で吸収する場合(あるいは大気抵抗係数 $C_D$ や太陽輻射圧係数 $C_R$ を「考慮パラメータ(consider parameter)」として不確かさに含める場合)、$Q$ が大きいほど $P_k^{-}$ は観測のない期間に急速に膨らみます。

摂動モデルが精密であるほど、この両方が改善されます。バイアスが減り、$Q$ を小さく保てるため、次の観測(次の追跡パス)が来るまでの予測軌道の精度が高く保たれるわけです。これは実務上、DSN運用スケジューリングやアンテナ自動追尾と直接つながっています。

• 追跡パスの合間、探査機は数時間から数日にわたって地上局の視界の外にいます。この間、地上局は次のパスの開始時刻に、精密な事前軌道予測(プリディクト、predict)に基づいてアンテナを目的の方向へあらかじめ向けておく必要があります。摂動モデルが粗く予測誤差(角度誤差)が大きいと、アンテナのビーム幅からの外れが大きくなり、アンテナ自動追尾で学んだモノパルス方式などの捕捉(アクイジション)手続きに、より広い探索範囲とより長い時間が必要になります。
• 逆に摂動モデルが高精度であれば、初期捕捉の不確実性(サーチコーン)を狭くでき、捕捉時間を短縮できます。これは1回の追跡パスの中で科学データ取得に使える時間を実質的に増やすことに直結し、DSN運用スケジューリングで見た「有限のアンテナ資源を数十のミッションが取り合う」制約充足問題において、各ミッションが要求する必要パス時間・許容誤差そのものを左右する要因になります。
• さらに、追跡パスの間隔自体も摂動モデルの精度に依存して決められます。高精度な摂動モデルがあれば、観測データが得られない期間が長くても予測共分散 $P_k^{-}$ の増大が緩やかなため、追跡パスの頻度を減らしても許容誤差内に軌道推定を保てます。これは、限られたDSNアンテナ資源をより多くのミッションに割り当てる余地を生む、運用上の重要なトレードオフです。

まとめると、摂動モデルの精度と、必要な追跡データ更新頻度は、片方を上げればもう片方を下げられるという意味でトレードオフの関係にある、と言えます。

実務での使われ方

実際の軌道決定・軌道伝播ソフトウェアでは、摂動モデルの複雑さにいくつかの明確なレベルがあり、ミッションフェーズや用途に応じて使い分けられています。

レベル1: 解析的簡易モデル(SGP4/SDP4)。 NORAD(現・米宇宙軍18th SDS)が公開するTLE(Two-Line Element set)データと組み合わせて使われる SGP4/SDP4(Simplified General Perturbations) は、$J_2$ を中心とした主要な摂動項を解析的な補正式に押し込み、数値積分なしで高速に軌道を伝播できるモデルです。精度はTLEの鮮度にもよりますが数日で数km程度の誤差が生じ得る粗いモデルで、詳細な大気モデルや第三体重力は簡略化されています。数万個規模のカタログ天体(デブリを含む)を扱う宇宙状況把握(SSA)、接近解析(コンジャンクション・スクリーニング)、大まかな可視パス予測など、速度と網羅性が精度より優先される用途で広く使われています。

レベル2: 中精度の数値積分モデル。 Cowellの方法(摂動を含めた運動方程式をそのまま数値積分する方法)に、$J_2$〜$J_6$程度の帯調和項、簡易な指数大気モデルまたはHarris-Priesterモデル、平板1枚のキャノンボールモデルによる太陽輻射圧を組み合わせたもので、多くのミッション運用における日常的な軌道決定・追跡予測に使われる標準的な精度クラスです。

レベル3: 高精度モデル。 精密科学ミッションや、月・惑星探査機の確定軌道決定には、重力場を数十〜数百次の球面調和関数(たとえばEGM2008の低次項を切り出したもの、あるいは重力科学ミッションで直接推定された天体固有の重力場モデル)まで展開し、大気モデルもNRLMSISE-00やJacchia-Bowman(JB2008)のような太陽活動・地磁気指数を反映した精密モデルを使い、太陽輻射圧も探査機の実形状(パネルごとの反射率・向き)を考慮したマクロモデルで扱い、第三体重力もJPLのDE(Development Ephemeris)シリーズのような高精度暦を用いる、というレベルまで精緻化されます。相対論的補正(シャピロ遅延、測地歳差などの一般相対性理論の効果)を含めることもあります。JPLのMONTEやODP、NASAのGMAT (General Mission Analysis Tool)、商用のSTK/Astrogatorといった軌道決定・解析ソフトウェアが、このレベルの精密モデルを実装しています。

ミッションフェーズによる使い分け。 深宇宙巡航フェーズでは、大気抵抗は無関係になり、支配的な摂動は太陽輻射圧と(接近時の)惑星の重力場になるため、モデルはむしろ単純化されます。一方、惑星周回フェーズ、特に低高度で科学観測を行うフェーズ(たとえば火星や木星の重力場マッピングミッション)では、重力科学の回で見たように、摂動そのものが観測対象になるほど精密なモデルが要求されます。また、着陸・ランデブーのような近接運用フェーズでは、そもそも中心天体からの重力・摂動論の枠組みよりも、相対運動の力学(Clohessy-Wiltshire方程式など)が前面に出てきます。これらの精度レベルの選択自体が、ミッション設計における重要なトレードオフの1つです。軌道データの交換にはCCSDS 502.0-B (Orbit Data Messages) のような標準フォーマットが使われ、状態ベクトルだけでなく、どの摂動モデルでどの座標系・時刻系に基づいて生成されたデータかというメタ情報も明示的に記述されるようになっています。

演習問題

1. 高度400 km(地球中心からの距離 $r \approx 6778\ \text{km}$)を飛ぶ衛星について、本文の近似式 $a_{J_2} \sim \frac{3}{2}J_2(R_E/r)^2(GM/r^2)$ を用いて $J_2$ 摂動加速度の大きさ(絶対値、$\text{m/s}^2$)を計算してください($GM_{\text{地球}} \approx 3.986\times10^{14}\ \text{m}^3/\text{s}^2$、$R_E \approx 6378\ \text{km}$、$J_2 \approx 1.08263\times10^{-3}$)。また、これを本文の表にある大気抵抗の相対的な大きさ($\sim 7\times10^{-7}$、中心重力加速度に対する比)と比較し、なぜ低軌道の軌道決定において $J_2$ を最優先でモデルに組み込む必要があるのか説明してください。

2. 半長径 $a = 7071\ \text{km}$(高度約700 kmのほぼ円軌道、$p\approx a$ と近似してよい)の太陽同期軌道を設計します。太陽同期条件は、昇交点赤経の歳差運動率が地球の公転角速度 $\dot\Omega = 360^\circ/365.2422\ \text{day} \approx 0.9856\ \text{deg/day}$ に一致することです。平均運動 $n = \sqrt{GM/a^3}$ を計算し(単位はrad/sからdeg/dayに変換すること)、本文の式 $\dot\Omega \approx -\frac{3}{2}nJ_2(R_E/p)^2\cos i$ を使って、この条件を満たす軌道傾斜角 $i$ を求めてください。$\cos i$ の符号がどうなるか(順行軌道か逆行軌道か)にも注目し、その物理的な意味を説明してください。

3. 静止軌道(GEO、地球中心からの距離 $r = 42164\ \text{km}$)を回る通信衛星を考えます。本文の第三体重力の近似式 $a_{3\text{rd}} \sim 2GM_{3\text{rd}}r/d^3$ を使い、月による摂動加速度のオーダーを見積もってください($GM_{\text{月}} \approx 4.905\times10^{12}\ \text{m}^3/\text{s}^2$、地球〜月の平均距離 $d \approx 3.844\times10^{8}\ \text{m}$)。この値を、GEOにおける中心重力加速度 $GM_{\text{地球}}/r^2$ と比較した相対的な大きさとして表し、本文の表の値($\sim3.2\times10^{-5}$)とおおむね一致することを確認してください。GEO衛星の運用で、南北方向のステーションキーピングに毎年比較的大きなデルタVが必要になる理由を、この計算結果から説明してください。

4. 摂動モデルの精度を落とし、$F$(あるいは非線形伝播関数 $f$)に $J_2$ を含めない粗い二体問題モデルだけでカルマンフィルタによる軌道決定を行ったとします。このとき、(a) 予測共分散 $P_k^{-}$ の増大の仕方、(b) 観測残差(イノベーション)の統計的な振る舞い、(c) 次回の追跡パスにおけるアンテナのビーム捕捉(アクイジション)にどのような影響が出るか、それぞれ前回とアンテナ自動追尾の回で学んだ内容を踏まえて説明してください。

まとめと次回予告

今回は、カルマンフィルタの状態遷移行列 $F$ の中身、すなわち軌道力学モデルそのものを精緻化しました。二体問題の理想的なケプラー運動から出発し、地球の非球対称重力($J_2$項)、大気抵抗、太陽輻射圧、第三体重力という4つの主要な摂動源を、それぞれの数式とオーダーで比較し、軌道高度によって支配的な摂動要因が大きく異なることを見ました。特に $J_2$ 摂動が昇交点赤経・近地点引数の永年的な歳差運動を引き起こすこと、そしてこれが太陽同期軌道のような実用的な軌道設計に積極的に利用されていることも確認しました。最後に、摂動モデルの精度が、追跡パス間の予測精度、ひいてはアンテナ捕捉やDSNスケジューリングといった運用そのものに直結するトレードオフの関係にあることを見ました。

ここまでのシリーズで、測距・ドップラーといった電波観測量から軌道を求める推定アルゴリズム(最小二乗法、カルマンフィルタ)と、その推定を支える力学モデル(二体問題とその摂動)の両方を手に入れました。しかし、軌道決定の精度を語る上でもう1つ欠かせない観測技術があります。次回は、複数の地上局で同時に受信した信号の到達時刻差から、驚くほど高い角度分解能で探査機の方向を決定する VLBI(超長基線電波干渉法)測地学を、その基線ベクトルや大気遅延補正の考え方とともに詳しく見ていきます。

参考文献

• D. A. Vallado, Fundamentals of Astrodynamics and Applications, 4th ed., Microcosm Press/Springer
• O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods and Applications, Springer
• B. D. Tapley, B. E. Schutz, G. H. Born, Statistical Orbit Determination, Elsevier Academic Press
• CCSDS 502.0-B, Orbit Data Messages, Blue Book
• J. H. Yuen (ed.), Deep Space Telecommunications Systems Engineering, JPL Publication 82-76
• D. Vallado and P. Crawford, “SGP4 Orbit Determination,” AIAA 2008-6770