オンボード自律軌道決定 — 探査機自身が地上を頼らず「今どこにいるか」を知る

この回で学ぶこと

カルマンフィルタによる軌道決定の回では、地上局が受け取ったレンジング・ドップラー観測を、JPLのODPやMONTEのような大型計算機上でカルマンフィルタにかけ、逐次的に軌道を更新していく仕組みを学びました。あの議論には、あえて明示していなかった前提が1つあります。フィルタを回している計算機は地上にあり、演算能力にもメモリにもほぼ制約がないという前提です。倍精度(あるいはそれ以上)の浮動小数点演算、ギガバイト単位のメモリ、必要ならスレッドを並列に走らせられる余裕——地上の運用センターでは、こうした計算資源をほぼ無尽蔵に使えます。

しかし、探査機自身が「今、自分はどこにいて、どちらに向かって、どれくらいの速度で動いているか」をリアルタイムに知らなければならない場面があります。

• 小惑星近傍でのタッチダウン: 目標天体表面のごく近くまで自律的に降下し、危険な地形を避けながら着地点を選ぶ局面では、地球からの光行時間(小惑星リュウグウなら片道約20分前後)が判断のたびに乗ってしまう地上遠隔操縦はそもそも成立しません。
• 惑星への自律着陸(EDL: Entry, Descent, and Landing): 火星の大気突入から着陸まではわずか数分間で完結し、しかもその間は地球との通信そのものが困難・不可能なことが多い「セブン・ミニッツ・オブ・テラー」のフェーズです。
• 深宇宙での長期自律運用: 通信リンクが太陽合(ソーラーコンジャンクション、前々回近くの話題)などで長期間途絶する場合や、往復の光行時間があまりに長く、地上からのリアルタイム制御が原理的に不可能な将来の深宇宙ミッションでは、探査機は数日から数週間、自分自身の判断だけで航法・誘導・制御を継続する必要があります。

この回では、こうしたオンボード自律軌道決定(onboard autonomous orbit determination)を扱います。強調しておきたいのは、ここで使う数学的な道具立てそのもの——状態空間モデル、予測ステップと更新ステップ、カルマンゲイン——は前回学んだものとまったく同じだということです。ベクトル $\hat{\mathbf{x}}_k$ や行列 $P_k, K_k$ の意味も式変形も一切変わりません。変わるのは、そのフィルタをどこで、何を使って、どれだけの制約の中で回すかという工学的な文脈です。この回ではカルマンフィルタの基本式を再導出せず、探査機搭載コンピュータという特殊な実行環境がフィルタ設計に強いる3つの制約——(1) 状態次元と演算量のトレードオフ、(2) 搭載センサーの統合(センサーフュージョン)、(3) 単精度演算下での数値的安定性——を順に見ていきます。

直感的導入 — 同じ数式、まったく違う実行環境

前回のカルマンフィルタの5本の式をもう一度眺めてみましょう。

$$\hat{\mathbf{x}}_k^{-} = F_{k-1}\hat{\mathbf{x}}_{k-1}, \qquad P_k^{-} = F_{k-1}P_{k-1}F_{k-1}^\top + Q_{k-1}$$ $$K_k = P_k^{-}H_k^\top\big(H_k P_k^{-}H_k^\top + R_k\big)^{-1}$$ $$\hat{\mathbf{x}}_k = \hat{\mathbf{x}}_k^{-} + K_k\big(\mathbf{z}_k - H_k\hat{\mathbf{x}}_k^{-}\big), \qquad P_k = (I-K_k H_k)P_k^{-}$$

数式だけを見れば、これを実行するのに必要なのは「行列の掛け算と足し算、それに1回の逆行列計算」だけです。地上の計算機であれば、状態次元 $n$ が数百・数千に膨らんでも何の問題もありません。ところが、これを深宇宙探査機の搭載コンピュータ(オンボードコンピュータ、OBC)で動かそうとすると、事情はまったく違ってきます。

代表的な耐放射線フライトコンピュータであるRAD750(BAE Systems製、火星探査車キュリオシティやパーサヴィアランス、多数の周回機・探査機に搭載)は、動作クロック約110〜200MHz、実効演算性能はおおむね数百MIPS(1秒あたり数億命令)程度です。これは、みなさんの手元のスマートフォンや、地上局の運用サーバーと比べて3桁から4桁近く遅い計算資源です。しかもこのCPUは軌道決定のためだけに存在するわけではなく、姿勢制御、通信、電源管理、故障診断、科学観測機器の制御など、探査機のあらゆるサブシステムのソフトウェアを同時にリアルタイムでこなさなければなりません。搭載メモリも、宇宙放射線によるビット反転(シングルイベントアップセット)対策のECC(誤り訂正符号)付きメモリを使うため実効容量がさらに削られ、地上のワークステーションとは比較にならないほど貴重な資源です。

つまり、オンボード自律軌道決定の設計とは、「地上でやっていることと数学的に同じことを、極端に貧しい計算資源の中で、かつ他のタスクと計算資源を奪い合いながら、それでも要求されるリアルタイム性を満たすようにやり遂げる」という、純粋にエンジニアリング上の最適化問題なのです。

数式的な議論

状態次元と演算量のトレードオフ

まず、カルマンフィルタ1サイクルの演算量を数え上げてみましょう。状態次元を $n$、観測次元を $m$ とします。

予測ステップの共分散伝播 $P_k^{-} = F_{k-1}P_{k-1}F_{k-1}^\top + Q_{k-1}$ は、$n\times n$ 行列同士の掛け算を2回行う演算であり、その計算量は

$$C_{\text{predict}} = O(n^3)$$

です。更新ステップでは、イノベーション共分散 $S_k = H_kP_k^{-}H_k^\top + R_k$ の計算に $O(n^2 m + nm^2)$、その逆行列 $S_k^{-1}$ の計算に $O(m^3)$、カルマンゲイン $K_k$ の計算とそれを使った状態・共分散の更新にさらに $O(n^2 m)$ 程度の演算が必要です。まとめると、1サイクルあたりの総演算量はおおむね

$$C_{\text{cycle}} = O(n^3) + O(n^2 m + nm^2 + m^3)$$

典型的な軌道決定フィルタでは観測次元 $m$(距離・視線方向などの観測数)は状態次元 $n$(位置・速度・バイアスなどの状態数)よりずっと小さいため、演算量は状態次元 $n$ の3乗、すなわち $P_k^{-}$ の伝播が支配的になります。これがオンボードフィルタ設計における最も重要な事実です。状態次元を1つ増やすたびに、演算コストは線形にではなく3乗で効いてくる。

具体的な数字で見てみましょう。位置3成分・速度3成分の基本的な状態(6次元、$n=6$)に、IMUのバイアス3成分を追加して $n=9$ にしたとします。演算コストの比は

$$\left(\frac{9}{6}\right)^3 = \left(\frac{3}{2}\right)^3 = 3.375$$

状態次元をたった1.5倍にしただけで、共分散伝播の演算コストは3.4倍近くに膨れ上がります。地上の計算機ならこの違いを気にする理由はほとんどありませんが、限られたCPU時間を姿勢制御ループや故障診断タスクと分け合っている搭載コンピュータでは、この3乗則が設計上の死活問題になります。オンボードフィルタの設計者は、したがって次のような判断を常に迫られます。

• そもそも推定すべき状態は何か(位置・速度だけで十分か、それともセンサーバイアスや姿勢誤差まで状態に含めて同時推定すべきか)を最小限に絞り込む。
• 更新ステップの周期 $\Delta t$ をどこまで長く取れるか(更新頻度 $f=1/\Delta t$ を下げれば、単位時間あたりの総演算負荷は下がるが、その分カルマンフィルタが予測に頼る時間が長くなり $P_k^{-}$ の不確かさの増大——前回見た $P_k^{-}=F P_{k-1}F^\top+Q$ の $Q$ の蓄積——を許容しなければならない)。
• 状態や共分散行列そのものを、複数のサブシステム(たとえば姿勢と軌道)に分割し、それぞれ小さい次元のフィルタを独立に回してから緩く結合する分割・階層化フィルタ(federated filter、分散カルマンフィルタとも呼ばれる設計思想)を使うか。

地上で軌道決定を行う場合、フィルタ設計者は「精度をどこまで追求できるか」を主に考えます。搭載コンピュータの場合は、精度の追求と演算量・メモリ・電力という有限な予算の中でリアルタイム性を確保することが同格の設計目標になる、というのがオンボード自律軌道決定の最初の、そして最も本質的な違いです。

観測源の違いとセンサーフュージョン

地上局での軌道決定では、観測 $\mathbf{z}_k$ はほぼレンジング(測距)とドップラー(レンジレート)、それにDDORのような角度情報でした。探査機搭載のフィルタでは、地球との電波リンクが常に使えるとは限らず(むしろ使えない局面のためにこそオンボードOD が必要なのでした)、代わりに探査機自身が携えているセンサーからの観測を統合しなければなりません。主なものを整理すると次のようになります。

光学航法カメラ: 前に学んだように、恒星や近傍天体を撮像し、視線方向の単位ベクトル $\hat{s}_{\text{cam}}$ を観測量として得ます。これは状態(位置)に対して非線形な観測関数であり、

$$\mathbf{z}_k^{\text{cam}} = h_{\text{cam}}(\mathbf{x}_k) + \mathbf{v}_k^{\text{cam}}, \qquad h_{\text{cam}}(\mathbf{x}) = \frac{\mathbf{r}_{\text{target}} - \mathbf{r}(\mathbf{x})}{\lVert \mathbf{r}_{\text{target}} - \mathbf{r}(\mathbf{x}) \rVert}$$

のように、目標天体の位置ベクトル $\mathbf{r}_{\text{target}}$ と自分の位置 $\mathbf{r}(\mathbf{x})$ の差を正規化した方向ベクトルとして表せます。これは距離情報を直接には含まない「方位のみ(bearing-only)」の観測であり、単独では距離が決まらないため、探査機自身の運動(視差)や複数天体の同時観測と組み合わせて使う必要があります。

慣性計測装置(IMU): ジャイロスコープと加速度計を組み合わせたIMUは、他の2つのセンサーとは性質が異なります。IMUは「今どこにいるか」を直接観測するのではなく、角速度と比力(重力以外の加速度)を測るセンサーです。オンボードフィルタでは、IMU出力は観測 $\mathbf{z}_k$ としてではなく、状態遷移方程式の**制御入力(強制項)**として使われるのが一般的です。

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

ここで $\mathbf{u}_{k-1}$ がIMUから得た角速度・加速度(ストラップダウン慣性航法における積分入力)です。これによって、光学航法カメラの観測が得られない(たとえば降下中の一瞬、カメラが照明条件やブラーで使えない)期間も、IMUの積分だけで状態を高頻度に伝播させ続けることができます。ただしIMUにはバイアス(温度や振動で緩やかに変化するオフセット誤差)やランダムウォーク誤差が必ず伴い、これを補正しないまま積分し続けると誤差が時間とともに際限なく蓄積します。この対策として、IMUバイアス $\mathbf{b}$ を状態ベクトルに追加し(状態拡大、state augmentation)、

$$\tilde{\mathbf{x}}_k = \begin{bmatrix} \mathbf{x}_k \\ \mathbf{b}_k \end{bmatrix}, \qquad \mathbf{b}_k = \mathbf{b}_{k-1} + \boldsymbol{\eta}_{k-1}, \quad \boldsymbol{\eta}_{k-1}\sim\mathcal{N}(\mathbf{0}, Q_b)$$

のようにバイアス自体をゆっくり変動する確率過程(ランダムウォーク)としてモデル化し、他の観測(カメラなど)が得られたタイミングでバイアスもろとも推定・補正するのが標準的な手法です。これがまさに前節で述べた「状態次元を増やすと演算コストが3乗で効いてくる」典型例であり、IMUバイアスをどこまで細かく状態化するか(軸ごとに独立か、温度依存性まで含めるかなど)は、演算予算と要求精度のせめぎ合いで決まります。

GNSS受信機: 地球周回衛星や月・火星近傍での運用など、GPSやGLONASS、あるいは将来の月測位網(LunaNet/Moonlightなど)の電波が届く範囲では、GNSS受信機は位置・速度をほぼ直接観測できる強力なセンサーです。ただし真の意味での深宇宙(地球や他の測位衛星群から遠く離れた惑星間空間や小惑星近傍)では、そもそもGNSS信号が届かないため使えません。GNSSが使える局面では観測モデルは線形に近く扱いやすいのに対し、使えない局面(多くの深宇宙ミッション)では光学航法とIMUの統合が自律航法の主力にならざるを得ない、という住み分けがある点は覚えておく価値があります。

これら性質の異なる観測を1つのカルマンフィルタに統合する考え方をセンサーフュージョンと呼びます。数式の上では、観測ベクトル $\mathbf{z}_k$ とその共分散 $R_k$ をセンサーごとにブロック化して1本の更新式にまとめる(バッチ的な統合更新)ことも、あるいは各センサーの観測が届いた時刻ごとに1つずつ逐次的に更新ステップを繰り返す(逐次統合更新)こともできます。

$$\mathbf{z}_k = \begin{bmatrix} \mathbf{z}_k^{\text{cam}} \\ \mathbf{z}_k^{\text{GNSS}} \\ \vdots \end{bmatrix}, \qquad H_k = \begin{bmatrix} H_k^{\text{cam}} \\ H_k^{\text{GNSS}} \\ \vdots \end{bmatrix}, \qquad R_k = \begin{bmatrix} R_k^{\text{cam}} & & \\ & R_k^{\text{GNSS}} & \\ & & \ddots \end{bmatrix}$$

搭載コンピュータでは、観測共分散 $R_k$ が対角行列(異なるセンサー間の誤差が無相関)であるとき、$m\times m$ の逆行列計算 $S_k^{-1}$ を1回のまとめた行列演算として行うのではなく、観測1成分ずつをスカラーの逐次更新として処理する実装がしばしば好まれます。スカラー更新であれば逆行列は単なる除算になり、$O(m^3)$ の行列逆演算そのものを避けられるためです。また、各センサーのサンプリングレートは大きく異なります(IMUは数百Hz、光学航法カメラは数秒に1回、GNSSは1Hz程度など)。このため実際のオンボードフィルタは、一定周期で全観測を同時に処理する単純な形ではなく、観測が届いた時刻ごとに、その時刻まで状態を予測してから、その1つの観測だけで更新するという非同期・可変レートの実装(マルチレートカルマンフィルタ)になるのが普通です。

数値的安定性 — 単精度演算がもたらす特有の困難

最後に、地上の軌道決定ではほとんど意識する必要のない問題を扱います。搭載コンピュータでは、消費電力・処理速度・耐放射線性の制約から、単精度浮動小数点(32ビット)、場合によってはさらに精度の低い固定小数点演算しか使えないことがあります。倍精度(64ビット)で悠々と回せる地上のフィルタとは異なり、有効桁数がおよそ7桁程度しかない単精度演算では、カルマンフィルタの繰り返し計算のなかで丸め誤差が着実に蓄積していきます。

とりわけ深刻なのが共分散更新式です。前回導出した基本形

$$P_k = (I - K_k H_k)\, P_k^{-}$$

は数学的には正しいのですが、これを浮動小数点演算でそのまま計算すると、丸め誤差の影響で更新後の $P_k$ が対称行列でなくなったり、固有値が負になったりすることがあります。共分散行列は本来、定義上つねに半正定値(すべての固有値が0以上)でなければなりません。固有値が負になるということは、「ある状態成分の不確かさの2乗である分散がマイナスになる」という、物理的にありえない状態にフィルタが陥ったことを意味し、そこから先の推定は発散してしまいます。この現象は演算精度が高い地上の計算機でも原理上は起こり得ますが、単精度演算下ではるかに顕在化しやすくなります。

これに対する古典的な対策の1つが、更新式をより数値的に頑健な形——Joseph form(ジョセフ形)——で計算することです。

$$P_k = (I - K_k H_k)\, P_k^{-}\, (I - K_k H_k)^\top + K_k R_k K_k^\top$$

これは前回、カルマンゲインを導出する途中で経由した一般形そのものです。単純化された形 $P_k=(I-K_kH_k)P_k^{-}$ はカルマンゲイン $K_k$ が理論上正確な値であるという前提のもとで初めて成立する簡略形であり、丸め誤差で $K_k$ がわずかにずれると対称性・半正定性が崩れやすいのに対し、Joseph formは($K_k$ の値によらず)2次形式の和として構成されているため丸め誤差があっても半正定性が壊れにくいという性質を持ちます。演算量は単純形より増えますが、限られた精度でも安定に動作させたいオンボードフィルタでは、この余分な演算コストを払う価値があります。

さらに進んだ対策として、共分散行列 $P$ そのものを直接更新するのではなく、$P = SS^\top$ を満たす平方根行列 $S$(あるいは $P=UDU^\top$ という形の上三角行列 $U$ と対角行列 $D$)を状態として持ち歩き、$S$ や $U, D$ を更新する 平方根フィルタ(square-root filtering) や UDフィルタ(UD分解フィルタ) という手法があります。これらは、$P$ の固有値の平方根に相当する量を直接扱うため、条件数(最大固有値と最小固有値の比)がおよそ平方根倍に改善され、同じ精度の演算でもより深い桁落ちに耐えられるという利点があります。歴史的にも、これらの手法はまさに演算資源の乏しかった1960〜70年代の搭載コンピュータ——アポロ計画の誘導コンピュータや初期の航法衛星——のために開発されたものであり(平方根フィルタはJ. E. Potterによるアポロ計画向けの考案が起源とされ、UDフィルタはG. J. BiermanとC. L. Thorntonによって発展させられました)、オンボード自律軌道決定という文脈は、その誕生の経緯そのものに直結しています。単精度演算という制約の中でカルマンフィルタを破綻させずに回し続けるという課題は、宇宙工学の歴史を通じて繰り返し現れてきた、根の深い問題なのです。

実務での使われ方

オンボード自律航法・軌道決定を実運用で最初期に確立したのはNASA/JPLのAutoNavシステムです。1998年に打ち上げられた技術実証機ディープ・スペース1号(Deep Space 1)に初めて搭載され、探査機自身が搭載カメラで小惑星・彗星を撮像し、恒星背景に対する視差から自律的に軌道を推定して、地上からのコマンドなしに接近フライバイを行うことに成功しました。AutoNavはその後、彗星探査機スターダスト(Stardust)やディープ・インパクト(Deep Impact)、そしてイオンエンジンで小惑星ベスタ・準惑星ケレスを周回した**ドーン(Dawn)**にも受け継がれています。ドーンでは、イオンエンジンによる連続的な低推力軌道変更を行うため、地上でのバッチ処理による軌道決定・軌道設計のサイクルでは追いつかず、搭載コンピュータ上でほぼリアルタイムに軌道を推定し続ける自律航法が運用上不可欠でした。

日本の小惑星探査機はやぶさ2の小惑星リュウグウへのタッチダウン運用も、オンボード自律航法の代表例です。地球とリュウグウの間の光行時間は往復で数十分に達し、地上からリアルタイムに降下を制御することは不可能でした。そのため探査機は、あらかじめ表面に投下しておいた反射性のターゲットマーカーをカメラで追跡し、レーザー距離計(LRF)による高度情報と組み合わせて、自律的に降下速度と姿勢を制御しながらタッチダウンを行いました。地上の運用チームは、降下シーケンスそのものをあらかじめ設計・アップロードしておき、実際の降下中は探査機の自律判断に委ねるという運用形態が取られています。

火星着陸機のEDL(突入・降下・着陸)フェーズも、オンボード自律航法が必須の局面です。NASAの火星探査車**パーサヴィアランス(Mars 2020)**は、**LVS(Lander Vision System)と呼ばれる画像航法系を搭載し、降下中にカメラで撮影した地表の画像を、あらかじめ搭載しておいた基準地図とパターンマッチングすることで自己位置を特定する地形照合航法(TRN: Terrain Relative Navigation)**を実現しました。LVSでは、TRNによる画像ベースの位置推定とIMUによる高頻度な状態伝播をカルマンフィルタで統合し、危険な岩やクレーターを避けた着陸地点選定(Hazard Avoidance)にリアルタイムで反映させています。これはこの回で述べた「光学観測(低頻度・非線形)とIMU(高頻度・積分による伝播)をカルマンフィルタで統合する」という設計そのものの、現代における到達点と言えるでしょう。

演習問題

1. 状態次元を $n=6$(位置・速度)から、IMUのバイアス3成分と時計誤差1成分を追加して $n=10$ に拡張したとき、共分散伝播ステップ $P_k^{-}=F_{k-1}P_{k-1}F_{k-1}^\top+Q_{k-1}$ の演算コストはおよそ何倍になるか、$O(n^3)$ の比を使って見積もってください。またこの結果から、状態拡大を安易に行うことがオンボードフィルタにとってなぜ危険なのかを説明してください。

2. 光学航法カメラの観測関数 $h_{\text{cam}}(\mathbf{x}) = (\mathbf{r}_{\text{target}}-\mathbf{r}(\mathbf{x}))/\lVert \mathbf{r}_{\text{target}}-\mathbf{r}(\mathbf{x})\rVert$ が状態 $\mathbf{x}$(位置成分を含む)について非線形である理由を、正規化(ノルムで割る操作)に着目して説明し、この観測をオンボードのカルマンフィルタに組み込むには(前回学んだ拡張カルマンフィルタの考え方を踏まえて)何を計算する必要があるか述べてください。

3. 単純化された共分散更新式 $P_k=(I-K_kH_k)P_k^{-}$ とJoseph form $P_k=(I-K_kH_k)P_k^{-}(I-K_kH_k)^\top + K_kR_kK_k^\top$ について、$K_k$ が理論上の正確な値からわずかに $\delta K$ だけずれた($K_k \to K_k+\delta K$)場合を考え、どちらの式のほうが結果として得られる $P_k$ の対称性・半正定性が保たれやすいか、直感的な理由とともに論じてください(2次形式の和として書けているかどうかに着目)。

4. IMUのバイアスをなぜ観測ノイズとしてではなく、状態ベクトルに追加してランダムウォークとしてモデル化するのか、「積分によって誤差が時間とともに蓄積する」というIMU固有の性質に立ち返って説明してください。またこのようにバイアスを状態化することで、他のセンサー(カメラなど)からの観測が得られたときに何が同時に補正されることになるか述べてください。

まとめと次回予告

オンボード自律軌道決定は、地上局で行うカルマンフィルタによる軌道決定と数学的にはまったく同じ枠組みを使いながら、探査機搭載コンピュータという極端に制約された実行環境の中でそれを動かすという、独自のエンジニアリング上の課題を抱えています。状態次元の増加が演算コストを3乗で押し上げるというトレードオフ、光学航法カメラ・IMU・(使える局面では)GNSS受信機という性質の異なるセンサーをどう統合するかというセンサーフュージョンの設計、そして単精度演算下で共分散行列の半正定性を保ち続けるための数値的な工夫(Joseph form、平方根フィルタ、UDフィルタ)——これらはいずれも、地上の潤沢な計算資源の下では表面化しない、「実際にハードウェアの上でフィルタを動かす」ことに固有の問題でした。ディープ・スペース1号のAutoNavから、はやぶさ2のタッチダウン、そして火星2020のLVSに至るまで、この分野は深宇宙探査の自律性を支える基盤技術として発展を続けています。

これまでは、1機の探査機が単独で自分自身の軌道を推定する状況を前提にしてきました。しかし将来のミッションでは、複数の探査機が互いに近接して隊列を組む編隊飛行(フォーメーションフライング)が計画されており、そこでは「自分がどこにいるか」だけでなく「自分が相手に対してどこにいるか」という相対航法が中心的な課題になります。次回は、この編隊飛行における相対航法の考え方に軽く触れ、これまで見てきた絶対軌道決定の枠組みがどう拡張されるのかを概観します。

参考文献

• G. J. Bierman, Factorization Methods for Discrete Sequential Estimation, Academic Press, 1977
• M. S. Grewal, A. P. Andrews, Kalman Filtering: Theory and Practice Using MATLAB, Wiley
• S. Bhaskaran, “Autonomous Navigation for Deep Space Missions,” SpaceOps Conference, 2012
• A. Johnson et al., “Implementation of a Map Relative Localization System for Precision Landing on Mars 2020,” AAS Guidance, Navigation and Control Conference, 2020
• 「はやぶさ2」プロジェクトチーム, タッチダウン運用に関する技術報告(JAXA)
• B. D. Tapley, B. E. Schutz, G. H. Born, Statistical Orbit Determination, Elsevier Academic Press