軌道決定の不確かさの長期伝搬 — 観測が途切れた期間に誤差楕円体はどう育つか

この回で学ぶこと

カルマンフィルタの回で、私たちは観測が届くたびに推定値と誤差共分散を更新し続ける逐次推定の仕組みを学びました。あの議論の主役は「更新ステップ」でした。新しい観測が来るたびにカルマンゲインが働き、共分散 $P$ は縮んでいく——観測がある限り、フィルタは賢くなり続けます。

しかし現実の深宇宙運用では、観測は来ないのが普通です。DSN(深宇宙ネットワーク)のアンテナは数十のミッションが取り合う有限資源であり、1つの探査機に割り当てられる追跡パスは巡航フェーズなら週に数回、静穏期には週1回以下ということも珍しくありません。地球の自転による可視時間帯の制約、アンテナのメンテナンス、他ミッションのクリティカルイベント優先——こうした理由で、数日から数週間、観測がまったく更新されない期間が当たり前に発生します。

この期間、軌道推定はどうなっているのでしょうか。フィルタの言葉で言えば、更新ステップが一度も走らず、予測ステップ

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

だけが延々と繰り返される状態です。共分散は単調に育ち続けます。この回では、この「観測更新なしの長時間伝搬」を主役に据えて、次の問いに数式で答えていきます。

• 不確かさは時間とともにどの方向に、どんな速さで育つのか。
• 線形な共分散伝搬(ガウス分布の誤差楕円体)は、いつまで信用できるのか。なぜ長時間伝搬で誤差分布は「バナナ型」に歪むのか。
• 線形近似が破綻するとき、実務では何を使うのか(モンテカルロ法とそのコスト)。
• 育った不確かさは、次の追跡パスでのアンテナ捕捉や、フライバイ・軌道投入といった一発勝負のイベントの成否にどう跳ね返るのか。

カルマンフィルタが「観測がある間」の理論だったとすれば、今回は「観測がない間」の理論です。そして実務の航法チームが日々格闘しているのは、むしろこちらの時間帯なのです。

直感的な導入 — 観測の空白期間に何が起きるか

まず言葉でイメージを作りましょう。最後の追跡パスが終わった瞬間、フィルタの手元には推定値 $\hat{\mathbf{x}}_0$ と共分散 $P_0$ があります。位置誤差の標準偏差が数kmという、よく引き締まった状態だとします。ここから観測なしで1週間が経過すると、何が起きるか。

第一に、初期誤差そのものが力学によって引き伸ばされます。特に効くのが速度誤差です。速度が毎秒1mmずれていれば、1日(86400秒)後には単純計算でも86mの位置ずれになります。しかも軌道力学では話はもっと悪くて、後で見るように速度誤差は軌道エネルギー(半長径)の誤差を意味し、それが周回の速さそのものの誤差に化けるため、位置ずれは時間に比例してどこまでも育ちます。

第二に、モデル化しきれない力がじわじわ効きます。摂動の回で見たように、太陽輻射圧係数の誤差、微小なガスリーク、モデル外の加速度は、プロセスノイズ $Q$ として押し込まれています。観測がないあいだ、この $Q$ は毎ステップ足し込まれ続け、拭い去られることがありません。

第三に、これが今回の核心ですが、誤差が大きくなってくると誤差の「形」が変わり始めます。共分散行列で表せるのは楕円体(ガウス分布)だけですが、実際の誤差分布は軌道という曲がった道に沿って伸びるため、楕円体では表現できない曲がった、バナナのような形に歪んでいきます。共分散という道具そのものの限界に突き当たるのです。

順番に定式化していきましょう。

線形共分散伝搬の定式化

状態遷移行列による誤差の輸送

真の軌道力学は非線形な微分方程式 $\dot{\mathbf{x}} = f(\mathbf{x}, t)$ に従います($\mathbf{x}$ は位置・速度の6次元状態ベクトル)。参照軌道 $\mathbf{x}^*(t)$ のまわりの微小なずれ $\delta\mathbf{x}(t) = \mathbf{x}(t) - \mathbf{x}^*(t)$ は、1次のテイラー展開により線形な変分方程式

$$\delta\dot{\mathbf{x}}(t) = A(t)\,\delta\mathbf{x}(t), \qquad A(t) \equiv \left.\frac{\partial f}{\partial \mathbf{x}}\right|_{\mathbf{x}^*(t)}$$

に従います。この線形系の解は状態遷移行列(state transition matrix, STM) $\Phi(t, t_0)$ を使って

$$\delta\mathbf{x}(t) = \Phi(t, t_0)\,\delta\mathbf{x}(t_0)$$

と書け、$\Phi$ 自身は行列微分方程式

$$\frac{d}{dt}\Phi(t, t_0) = A(t)\,\Phi(t, t_0), \qquad \Phi(t_0, t_0) = I$$

を参照軌道と一緒に数値積分することで得られます。カルマンフィルタの回で $F_{k-1}$ と呼んでいた行列の正体はこの $\Phi(t_k, t_{k-1})$ です。

観測更新なしの共分散伝搬

初期誤差が平均ゼロ・共分散 $P(t_0)$ の確率変数だとすると、線形変換された $\delta\mathbf{x}(t)$ の共分散は

$$\boxed{P(t) = \Phi(t, t_0)\, P(t_0)\, \Phi(t, t_0)^\top + Q(t)}$$

となります。カルマンフィルタの予測ステップと同じ形ですが、意味合いが違います。フィルタでは $t_k - t_{k-1}$ は観測間隔(数秒〜数分)でしたが、いまは $t - t_0$ が数日から数週間です。1ステップの近似式ではなく、長い時間区間全体にわたる誤差の輸送を1本の式で表しています。

第2項の $Q(t)$ は、区間中に注入され続けるプロセスノイズの累積です。連続時間の白色雑音加速度(パワースペクトル密度 $q$、単位 $\text{m}^2/\text{s}^3$)としてモデル化すると、

$$Q(t) = \int_{t_0}^{t} \Phi(t, \tau)\, B\, q\, B^\top\, \Phi(t, \tau)^\top\, d\tau$$

という畳み込み積分になります($B$ は雑音が加速度として速度成分に入ることを表す行列)。この積分の構造を1次元(位置 $r$ と速度 $\dot r$)の単純な等速モデルで具体的に評価すると、よく知られた結果

$$Q(t) = q \begin{bmatrix} \dfrac{(t-t_0)^3}{3} & \dfrac{(t-t_0)^2}{2} \\[2mm] \dfrac{(t-t_0)^2}{2} & t-t_0 \end{bmatrix}$$

が得られます。ここで注目すべきはべきの違いです。速度分散は経過時間の1乗($q\,t$、ランダムウォーク)でしか育ちませんが、それが積分されて位置に届くため、位置分散は経過時間の3乗で育ちます。観測の空白期間が2倍になると、モデル化誤差起因の位置の不確かさ(標準偏差)は $2^{3/2} \approx 2.8$ 倍になる——空白期間の長さに対して不確かさの成長がいかに非線形かが、この $t^3$ に凝縮されています。

第1項も見ておきましょう。同じ1次元モデルで $\Phi = \begin{bmatrix} 1 & t-t_0 \\ 0 & 1 \end{bmatrix}$、初期共分散を対角 $\mathrm{diag}(\sigma_{r0}^2, \sigma_{v0}^2)$ とすると、位置分散は

$$\sigma_r^2(t) = \sigma_{r0}^2 + (t-t_0)^2\,\sigma_{v0}^2$$

となり、初期速度誤差が時間に比例して位置誤差に化けることが読み取れます。長時間伝搬では、初期位置誤差よりも初期速度誤差のほうがはるかに重い、というのが第一の教訓です。

誤差はどの方向に育つのか — 進行方向の永年成長

上の1次元モデルは本質を掴むには十分ですが、実際の軌道力学では誤差の成長には強い方向性があります。誤差を軌道に沿った局所座標系——動径方向(radial)、進行方向(along-track/transverse)、軌道面外方向(cross-track)の3軸、いわゆるRTN系——で見ると、圧倒的に速く育つのは進行方向です。その理由は軌道力学の基本関係から直接導けます。

ほぼ円軌道(半長径 $a$)を考えます。平均運動は $n = \sqrt{GM/a^3}$ なので、半長径に誤差 $\delta a$ があると平均運動に

$$\delta n = \frac{\partial n}{\partial a}\,\delta a = -\frac{3}{2}\,\frac{n}{a}\,\delta a$$

の誤差が生じます。周回の角速度そのものがずれているので、進行方向の位置誤差(弧長)は

$$\delta s(t) \approx a\,\delta n\,(t - t_0) = -\frac{3}{2}\, n\, \delta a\,(t - t_0)$$

と、**時間に比例して際限なく成長(永年成長)**します。動径方向の誤差 $\delta a$ 自体は振動しつつも有界にとどまるのに対し、それが引き起こす進行方向の誤差だけが線形に発散していく——これが軌道誤差伝搬の最も重要な非対称性です。

さらに、進行方向の速度誤差 $\delta v$ は vis-viva の関係(円軌道なら $v = na$)を通じて半長径誤差

$$\delta a = \frac{2}{n}\,\delta v$$

に直結するので、上の式と組み合わせると

$$\delta s(t) \approx -3\,\delta v\,(t - t_0)$$

というきわめて簡潔な結果になります。つまり進行方向の速度を $\delta v$ 測り損ねると、進行方向の位置誤差は毎秒 $3\,\delta v$ のペースで育つ。冒頭の素朴な見積もり($\delta v \times t$)の3倍速です。しかも符号が負であることに意味があります。速く動くように見誤った(実際より $a$ が小さい)探査機は実際には先に進んでいるのではなく、周期が短くなって先回りします。エネルギーを与えると遅れ、奪うと進むという、軌道ランデブーでもおなじみの直感に反する振る舞いです。

この進行方向の永年成長は、誤差楕円体の言葉で言えば「楕円体が進行方向にどんどん細長く引き伸ばされていく」ことを意味します。数日〜数週間の伝搬後の誤差楕円体は、進行方向の主軸が他の2軸より1〜2桁長い、極端に扁平な葉巻型になるのが典型です。そして**進行方向の位置誤差とは、要するに「探査機がその場所を通過する時刻のずれ」**です。長時間伝搬後の軌道予測誤差の実体は、「どこを通るか」より「いつそこを通るか」の誤差なのだ、という見方は実務でも非常に有効です。

非線形性の壁 — 誤差分布はバナナ型に歪む

ここまでの議論はすべて、参照軌道まわりの1次テイラー展開(線形化)の上に成り立っています。線形変換はガウス分布をガウス分布に写すので、誤差分布は常に楕円体で表せる——これが共分散伝搬の暗黙の前提でした。しかし誤差が育ってくると、この前提そのものが崩れます。

理由は幾何学的に明快です。進行方向に大きく引き伸ばされた誤差分布を考えましょう。分布の各点(それぞれが「あり得る探査機の位置」)は、直線上ではなく曲がった軌道に沿って分布します。進行方向に $\pm1000$ kmの広がりを持つ分布は、実際には軌道の曲率に従って弧を描きます。ガウス分布の等確率面は直線的な楕円体ですから、これを楕円体で近似すると、弧の中央部では分布を取りこぼし、弧の外側(実際には探査機がほぼ存在し得ない領域)を過大に含んでしまいます。数値実験で多数のサンプル点を非線形伝搬させて散布図を描くと、点群が三日月状に曲がって分布する様子が見られ、この形状から**バナナ型分布(banana-shaped distribution)**と呼ばれます。

もう少し形式的に言うと、線形化の誤差はテイラー展開の2次以降の項

$$\delta\mathbf{x}(t) = \Phi(t,t_0)\,\delta\mathbf{x}_0 + \frac{1}{2}\,\delta\mathbf{x}_0^\top\, \frac{\partial^2 \boldsymbol{\phi}}{\partial \mathbf{x}_0^2}\,\delta\mathbf{x}_0 + \cdots$$

($\boldsymbol{\phi}$ は非線形な解の流れ)を無視したことに由来します。2次項の寄与は $\delta\mathbf{x}_0$ の大きさの2乗に比例するため、誤差が小さいうちは線形近似は優秀だが、誤差が育つにつれて急速に劣化する。長時間伝搬では、(1) 初期誤差が力学で引き伸ばされて大きくなること、(2) 伝搬時間が長いほど流れの曲率の影響が蓄積することの両方が効いて、いつか必ず線形近似が破綻する時刻が来ます。ケプラー力学は角度変数(どこを回っているか)に対して強い非線形性を持つため、破綻はまさに進行方向——最も誤差が大きい方向——で最初に起きます。

非線形性がもたらす実務上の問題は2つあります。第一に、平均そのものがずれます。ガウス分布の初期誤差を非線形に伝搬させると、伝搬後の分布の平均は「初期平均を非線形伝搬した点」(参照軌道上の点)から系統的にずれます。バナナの重心は弧の内側に落ちるからです。第二に、共分散が分布の広がりを正しく表さなくなります。楕円体の3σ面の内側に本来の99.7%が入らなくなり、衝突確率や捕捉成功確率のような「分布の裾」に依存する計算が大きく狂います。

モンテカルロ法 — 力ずくの非線形伝搬とそのコスト

線形近似が信用できないとき、最も素直で最も信頼できる代替はモンテカルロ法です。初期分布 $\mathcal{N}(\hat{\mathbf{x}}_0, P_0)$ から $N$ 個のサンプル $\{\mathbf{x}_0^{(i)}\}_{i=1}^N$ を引き、それぞれを完全な非線形運動方程式で数値積分し、時刻 $t$ でのサンプル群 $\{\mathbf{x}^{(i)}(t)\}$ から標本平均・標本共分散、あるいは分布そのもの(ヒストグラム、確率の直接計数)を作ります。

$$\hat{\boldsymbol{\mu}}(t) = \frac{1}{N}\sum_{i=1}^{N} \mathbf{x}^{(i)}(t), \qquad \hat{P}(t) = \frac{1}{N-1}\sum_{i=1}^{N} \big(\mathbf{x}^{(i)}(t) - \hat{\boldsymbol{\mu}}\big)\big(\mathbf{x}^{(i)}(t) - \hat{\boldsymbol{\mu}}\big)^\top$$

線形化を一切使わないので、バナナ型の歪みも平均のずれも、力学モデルが正しい限り忠実に再現されます。問題は計算コストです。モンテカルロ推定の標準誤差は $1/\sqrt{N}$ でしか縮まないため、精度を1桁上げるにはサンプル数が100倍必要です。特に「3σ相当の稀な事象の確率」(たとえば捕捉失敗確率や衝突確率)を直接計数で求めたい場合、生起確率 $p$ の事象を相対誤差 $\epsilon$ で推定するのに必要なサンプル数はおよそ

$$N \approx \frac{1-p}{p\,\epsilon^2}$$

で、$p = 10^{-3}$ を10%の精度で見積もるだけで $N \sim 10^5$ 本の高精度軌道伝搬が必要になります。1本の伝搬が高次重力場・精密摂動モデル込みで数秒かかるとすれば、単純計算で数日分のCPU時間です(幸いモンテカルロは完全に並列化できるため、実務ではクラスタ計算で現実的な時間に収めます)。

このコストのため、実務では両極端の間を埋める中間的な手法も使われます。シグマポイントと呼ばれる少数(状態次元 $n$ に対して $2n+1$ 個)の代表点だけを非線形伝搬して平均・共分散を再構成するアンセンテッド変換(UKFの伝搬部と同じ考え方)、分布を複数のガウス成分の混合で表現し各成分は線形伝搬する**ガウス混合モデル(GMM)**などです。使い分けの相場観としては、日常的な追跡予測は線形共分散伝搬で十分、フライバイ設計や衝突確率評価のようなクリティカルな解析はモンテカルロで検証、その間をUT/GMMで補う、という階層になっています。

実務での使われ方

次パスの捕捉不確かさ予算 — 育った共分散をアンテナビームで受け止める

観測空白期間の終わり、次の追跡パスの開始時点で、地上局は伝搬された共分散 $P(t)$ が示す不確かさの中から探査機を「見つけ出す」必要があります。位置の不確かさ $\sigma_{pos}$ は距離 $\rho$ で割ることで角度の不確かさ $\sigma_{ang} \approx \sigma_{pos}/\rho$ に換算され、これがアンテナの**半値ビーム幅(HPBW)**と比較されます。DSNの34 m局のXバンド(8.4 GHz)のHPBWは約0.063°(およそ $1.1\times10^{-3}$ rad)です。深宇宙距離ではこの角度予算は通常余裕がありますが(演習問題3で確認します)、周波数軸の予算はずっとタイトです。進行方向の誤差は前述の通り「通過時刻のずれ」すなわち視線速度予測の誤差でもあり、ドップラー予測のずれが受信機のロックループの引き込み範囲を超えると、周波数掃引による探索が必要になって捕捉時間が延びます。つまり伝搬された共分散は、角度(ポインティング)と周波数(ドップラー予測)の両方の捕捉予算に変換されて次パスの運用計画に直接入力されるのです。要求される捕捉予算から逆算して「共分散がここまで育つ前に次のパスを入れる」という形で追跡スケジュール要求が決まる、という向きの使われ方もします。摂動の回で見た「モデル精度と追跡頻度のトレードオフ」の、共分散を介した定量版です。

フライバイ・軌道投入前の「最終OD更新から実行まで」の管理

不確かさ伝搬の管理が最もシビアになるのが、惑星フライバイや軌道投入(MOI/JOIなど)の直前です。軌道決定には観測データの取得・処理・解の収束確認・コマンド作成・アップリンクという一連の時間(データカットオフからマヌーバ実行まで、典型的には数時間から数日)が必ず必要で、最後のOD解が確定した瞬間から実行の瞬間まで、軌道推定は観測更新なしの伝搬に委ねられます。フライバイの照準精度は標的天体まわりのBプレーン(照準平面)上の誤差楕円で管理されますが、この楕円の最終的な大きさは「最終ODの共分散を実行時刻まで伝搬したもの+マヌーバ実行誤差」で決まります。データカットオフを遅らせれば伝搬時間が短くなり照準は正確になりますが、地上での処理・検証時間が圧迫されリスクが増える——このトレードオフの定量的な裏付けとして共分散伝搬解析が使われ、最終軌道修正マヌーバ(TCM)の実施時期やキャンセル判断(伝搬誤差がすでに許容楕円に収まっていればTCMを取りやめる)が下されます。カッシーニのタイタンフライバイ群や、はやぶさ2の小惑星近傍運用のような繰り返しの近接イベントでは、この「最終OD→伝搬→実行」のサイクル設計そのものが航法運用の骨格でした。

共分散の共有と「共分散リアリズム」

伝搬された共分散は機関の外にも流通します。CCSDSの軌道データ標準(CCSDS 502.0-B, Orbit Data Messages)では状態ベクトルとともに6×6共分散行列を交換でき、衝突リスク評価では**CDM(Conjunction Data Message, CCSDS 508.0-B)に両物体の伝搬共分散が載って衝突確率計算の入力になります。ここで深刻なのが共分散リアリズム(covariance realism)**の問題です。線形伝搬の限界、$Q$ の過小設定、モデル化誤差の見落としにより、報告される共分散が実際の誤差分布より小さすぎる(過信)ケースが実運用で繰り返し指摘されており、過信した共分散は衝突確率を桁で誤らせます。米宇宙軍やNASA CARA(Conjunction Assessment Risk Analysis)の運用では、実測残差と共分散の整合性を統計検定で監視し、必要ならスケーリングを施すという地道な検証が続けられています。「共分散は計算できる。しかしそれが正しいかは別問題」——長時間伝搬を扱う実務者の共通認識です。

演習問題

1. ほぼ円軌道(半長径 $a = 7000$ km、$GM = 3.986\times10^{14}\ \text{m}^3/\text{s}^2$)を回る衛星の進行方向速度に $\delta v = 1$ mm/s の推定誤差が残っているとします。(a) 平均運動 $n$ を計算してください。(b) 関係式 $\delta a = 2\delta v/n$ からこの速度誤差が意味する半長径誤差 $\delta a$ を求めてください。(c) 近似式 $\delta s \approx -3\,\delta v\,(t-t_0)$ を使い、観測空白が1日および7日続いたときの進行方向位置誤差をそれぞれ求めてください。動径方向の誤差(オーダーとして $\sim\delta a$ にとどまる)と比較し、誤差楕円体がどんな形に育つか述べてください。

2. プロセスノイズだけによる位置分散の成長 $\sigma_r^2(t) = q\,(t-t_0)^3/3$ を考えます。モデル化しきれない加速度をPSD $q = 1\times10^{-12}\ \text{m}^2/\text{s}^3$ の白色雑音とするとき、観測空白が1日、3日、7日の場合の位置標準偏差 $\sigma_r$ をそれぞれ計算し、空白期間が7倍になると $\sigma_r$ が何倍になるかを確認してください($t^{3/2}$ 則)。

3. 地球から $5\times10^{8}$ km の距離にいる深宇宙探査機の、伝搬後の位置不確かさ(1σ)が進行方向に3000 kmまで育ったとします。(a) これを地球から見た角度の不確かさ(rad)に換算し、DSN 34 m局Xバンドの半値ビーム幅 約 $1.1\times10^{-3}$ rad と比較してください。(b) 一方、地球周回衛星を距離2000 kmから追尾する地上局を考えます。進行方向の位置不確かさがわずか10 kmでも、角度の不確かさは何radになりますか。(a)(b)の比較から、「角度の捕捉予算が問題になるのは深宇宙と地球周回のどちらの運用か」「深宇宙で代わりにタイトになるのはどの予算か」を論じてください。

4. 発生確率 $p = 3\times10^{-3}$ 程度と見込まれる捕捉失敗イベントの確率を、モンテカルロ法の直接計数で相対誤差 $\epsilon = 0.1$(10%)で推定したいとします。本文の式 $N \approx (1-p)/(p\,\epsilon^2)$ から必要サンプル数 $N$ を見積もり、1本の軌道伝搬に2秒かかるとして単一CPUコアでの所要時間を概算してください。また、同じ計算を線形共分散伝搬(ガウス分布の裾の解析計算)で済ませた場合に潜むリスクを、バナナ型歪みの観点から説明してください。

まとめと次回予告

今回は、カルマンフィルタの「観測がない側」の理論として、軌道決定の不確かさの長期伝搬を扱いました。観測更新なしの共分散は $P(t) = \Phi P_0 \Phi^\top + Q(t)$ に従って単調に育ち、初期速度誤差は時間に比例して位置誤差に化け($t^2$ の分散成長)、プロセスノイズは $t^3$ の分散成長をもたらします。方向で見れば、半長径誤差が平均運動誤差に変換されることで進行方向の誤差だけが永年成長し($\delta s \approx -3\,\delta v\, t$)、誤差楕円体は進行方向に細長い葉巻型に、さらに誤差が育つと線形近似の及ばないバナナ型に歪みます。その先はモンテカルロ法(と、UT・GMMのような中間手法)の領分であり、精度とCPU時間のトレードオフでした。そしてこの伝搬された不確かさこそが、次パスの捕捉予算、フライバイ前の最終OD運用、共分散リアリズムといった航法実務の意思決定を駆動している——これが今回の全体像です。

次回は打って変わって信号の世界に戻り、測距符号の設計を扱います。測距やレンジングの回では擬似雑音(PN)符号を「相関の鋭い便利な符号」として使ってきましたが、そもそも相関特性の良い符号はどうやって作るのか。複数の探査機・衛星が同じ帯域を共有しても互いに混信しない符号族を系統的に構成できるGold符号を、その代数的な構成法と相関特性の理論限界から見ていきます。

参考文献

• B. D. Tapley, B. E. Schutz, G. H. Born, Statistical Orbit Determination, Elsevier Academic Press
• 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
• Y. Luo and Z. Yang, “A review of uncertainty propagation in orbital mechanics,” Progress in Aerospace Sciences, Vol. 89, 2017
• CCSDS 502.0-B, Orbit Data Messages; CCSDS 508.0-B, Conjunction Data Message
• DSN Telecommunications Link Design Handbook, DSN No. 810-005