スタートラッカーによる姿勢決定 — 恒星パターンから探査機の「向き」を知る

この回で学ぶこと

前回扱った光学航法(OPNAV)は、探査機が搭載カメラで近傍の小惑星や彗星を撮像し、視差を利用して自分の位置を求める技術でした。しかしその議論の中で、実はもう一つの問いにも軽く触れていました。恒星は実質的に無限遠にあるため、恒星への視線方向は探査機の位置にはよらず、探査機の向きだけで決まる、という性質です。今回はこの性質を正面から扱い、恒星を撮像することで探査機の姿勢(attitude)、すなわち「自分がどちらを向いているか」を求める**スタートラッカー(Star Tracker, STT)**を掘り下げます。

「位置」と「向き」はまったく別の情報です。探査機が太陽系のどこにいるかを完璧に知っていても、機体がどちらを向いているか分からなければ、高利得アンテナ(HGA)を地球に正確に向けることも、太陽電池パネルを太陽に向けて発電することも、光学航法カメラを目標天体に向けることもできません。これまでのレッスンで暗黙のうちに前提としてきた「アンテナが正しい方向を向いている」「太陽電池が発電できている」という条件は、すべて姿勢決定・姿勢制御(Attitude Determination and Control, ADC)という土台の上に成り立っています。姿勢が分からなければ、探査機は文字通り宇宙の中で自分の向きを見失った、目隠しをされた状態になってしまいます。

具体的には、(1) 恒星パターンをカタログと照合して姿勢を特定するロスト・イン・スペース(lost-in-space)アルゴリズムの仕組み、(2) 回転行列とクォータニオン(四元数)による3次元姿勢表現、そしてなぜ宇宙工学ではオイラー角ではなくクォータニオンが好まれるのか(ジンバルロックの回避)、(3) 複数の恒星観測から最適な姿勢を推定する最小二乗問題であるWahbaの問題とその代表的な解法(TRIAD法・QUEST法)、そして(4) スタートラッカー単体の弱点をジャイロスコープで補うセンサーフュージョンの考え方を、順を追って数式で確認していきます。

直感的導入 — 「どこにいるか」ではなく「どちらを向いているか」

大航海時代の航海士は、六分儀で恒星の高度を測って自分の船の位置(緯度)を求めていました。スタートラッカーがやっていることは、これとは似て非なるものです。恒星は探査機からの距離によらず天球上のほぼ同じ位置に見えるため、恒星をいくら観測しても探査機が太陽系のどこにいるかは分かりません。その代わり、恒星が画像のどの位置にどんな配置で写っているかを見れば、「カメラが天球に対してどちらを向いているか」、つまり探査機の姿勢が一意に定まります。

その原理は、前回紹介した恒星ペアの角距離という「指紋」の話がそのまま使えます。夜空を見上げてオリオン座やカシオペア座を見分けられるのは、星々の相対的な位置関係(パターン)を知っているからです。スタートラッカーも同じ発想で、撮像した恒星どうしの角距離の組み合わせを、地上であらかじめ作っておいた恒星カタログのデータベースと照合し、「今画像に写っているのはどの恒星か」を、探査機の姿勢について何の事前情報もない状態から特定します。これがロスト・イン・スペースアルゴリズムです。

しかし「どの恒星がどこに写っているか分かった」だけでは終わりません。カメラ座標系(探査機本体に対して固定された座標系)で見た恒星の方向と、カタログが与える慣性座標系(天球基準系)での恒星の方向という、2つの異なる座標系での同じ恒星の記述を、矛盾なく結びつける回転を求める必要があります。この「複数のベクトル観測から最適な回転を求める」という問題こそが、姿勢決定の数学的な核心です。そしてその回転を計算機の中でどう表現するかによって、計算のしやすさも、扱いにくい落とし穴(ジンバルロック)の有無も大きく変わってきます。

数式による定式化

姿勢を表す数学(1): 回転行列

探査機の姿勢とは、数学的には、慣性座標系(天球基準系)で表されたベクトル $\vec r$ を、探査機本体に固定された座標系(ボディ座標系)で表したベクトル $\vec s$ に変換する線形写像のことです。この変換は $3\times3$ の回転行列(方向余弦行列) $A$ によって

$$\vec s = A\,\vec r$$

と書けます。回転行列は長さと角度を保存する変換(剛体の回転)を表すため、その列ベクトル(あるいは行ベクトル)は互いに直交する単位ベクトルの組でなければならず、次の直交性の条件を満たします。

$$A^{\mathsf T} A = A A^{\mathsf T} = I, \qquad \det A = +1$$

($\det A = -1$ となる直交行列は回転ではなく鏡映を含む変換なので除外されます。)この条件を満たす行列全体の集合を特殊直交群 $SO(3)$ と呼び、探査機の姿勢は $SO(3)$ の元として表現されます。

回転行列は $3\times 3=9$ 個の成分を持ちますが、直交性の条件($A^{\mathsf T}A=I$)がもたらす束縛は6個(対称行列 $A^{\mathsf T}A - I$ の独立成分の数)あるため、実質的な自由度は $9-6=3$ しかありません。これは直感とも一致します。剛体の向きは「ロール・ピッチ・ヨー」のような3つの独立した角度で決まるはずだからです。この「9成分だが自由度3」という冗長性が、姿勢表現をめぐる工学的な議論の出発点になります。

姿勢を表す数学(2): オイラー角とジンバルロック

自由度が3なら、素直に3つの角度で姿勢を表現したくなります。これがオイラー角です。たとえば「ヨー角 $\psi$(機体をZ軸まわりに回転)→ピッチ角 $\theta$(新しいY軸まわりに回転)→ロール角 $\phi$(さらに新しいX軸まわりに回転)」という順番(3-2-1系列)で3回の基本回転を合成すれば、任意の姿勢を表現できます。回転行列はこの3つの基本回転行列の積 $A = A_x(\phi)A_y(\theta)A_z(\psi)$ として書けます。

この表現の直感的な分かりやすさは大きな利点ですが、致命的な弱点があります。角速度 $\vec\omega = (p,q,r)$(ボディ座標系での回転角速度)とオイラー角の時間微分の関係(運動学方程式)は、3-2-1系列の場合

$$\dot\phi = p + (q\sin\phi + r\cos\phi)\tan\theta, \qquad \dot\theta = q\cos\phi - r\sin\phi, \qquad \dot\psi = \frac{q\sin\phi + r\cos\phi}{\cos\theta}$$

となります。$\dot\phi$ と $\dot\psi$ の式に $\tan\theta$ および $1/\cos\theta$ が現れていることに注目してください。ピッチ角が $\theta \to \pm90^\circ$ に近づくと、この係数は発散します。このとき、ヨー軸まわりの回転とロール軸まわりの回転が実質的に同じ回転軸を指すようになってしまい(2つの回転軸が重なってしまい)、3自由度あったはずの回転が数学的に2自由度分しか表現できなくなります。これが**ジンバルロック(gimbal lock)**です。

ジンバルロックは物理的に姿勢制御ができなくなるわけではなく、あくまで「オイラー角というパラメータ化の選び方」がその姿勢の近傍で特異点を持つという、表現上の問題です。しかし探査機の姿勢制御則やオンボードのフィルタがこの運動学方程式を数値的に積分している場合、$\theta$ が $90^\circ$ に近づくと $1/\cos\theta$ が急激に増大し、数値誤差が爆発的に拡大するという実務上深刻な問題を引き起こします。姿勢がどの向きを取りうるか探査機の運用者が事前に完全に把握できるとは限らない以上、「ある特定の姿勢の近傍で計算が破綻する」表現方法を姿勢決定・姿勢制御の中核に据えるのは危険です。

姿勢を表す数学(3): クォータニオン(四元数)

この問題を回避するために宇宙工学で標準的に使われるのが**クォータニオン(四元数)**による姿勢表現です。クォータニオンはベクトル部 $\vec q_v = (q_1,q_2,q_3)$ とスカラー部 $q_4$ からなる4つの実数の組 $q = (q_1,q_2,q_3,q_4)$ で、単位ベクトル $\hat n$ のまわりに角度 $\Theta$ だけ回転する姿勢を

$$\vec q_v = \hat n \sin\!\left(\frac{\Theta}{2}\right), \qquad q_4 = \cos\!\left(\frac{\Theta}{2}\right)$$

として表します(この4成分・スカラー最後という並びは、航法データメッセージで見たCCSDSの姿勢メッセージ規格 $(Q_1,Q_2,Q_3,Q_C)$ とも整合する慣例です)。定義から直ちに単位ノルム拘束 $q_1^2+q_2^2+q_3^2+q_4^2=1$ が成り立ちます。この4成分・1拘束という組み合わせで、回転行列と同じ「実質自由度3」を過不足なく表現しています。

クォータニオンから回転行列を構成する式は

$$A(q) = (q_4^2 - |\vec q_v|^2)\,I + 2\,\vec q_v \vec q_v^{\mathsf T} - 2 q_4\,[\vec q_v \times]$$

で与えられます。ここで $[\vec q_v\times]$ は $\vec q_v$ との外積を表す歪対称行列 $[\vec q_v\times]\vec x = \vec q_v \times \vec x$ です。三角関数を含まない多項式演算だけで回転行列に変換できる点も、オンボードコンピュータでの実装上の利点です。

姿勢の時間変化(運動学方程式)も、クォータニオンでは特異点のないシンプルな線形微分方程式で書けます。

$$\dot q = \frac{1}{2}\,\Omega(\vec\omega)\, q, \qquad \Omega(\vec\omega) = \begin{pmatrix} 0 & r & -q & p \\ -r & 0 & p & q \\ q & -p & 0 & r \\ -p & -q & -r & 0 \end{pmatrix}$$

(ここでの $p,q,r$ はオイラー角の節と同じくボディ角速度成分で、行列 $\Omega$ の添字 $q$ とは別物であることに注意してください。)この式にはオイラー角のときのような $\tan\theta$ や $1/\cos\theta$ の項が一切現れず、姿勢がどんな向きであっても、係数が発散することはありません。これが、宇宙機の姿勢決定・姿勢制御においてオイラー角ではなくクォータニオンが標準的に採用される最大の理由です。実装上は、数値積分の誤差でノルムが1からわずかにずれていくため、定期的に $q \leftarrow q/|q|$ と正規化する処理が加えられます。

Wahbaの問題: 複数の恒星観測から最適な姿勢を求める

スタートラッカーが $N$ 個の恒星を検出し、それぞれについて前回導入したカメラ座標系(ボディ座標系)での視線方向単位ベクトル $\hat s_i$ と、恒星カタログが与える慣性座標系での方向 $\hat r_i$ の対応が(ロスト・イン・スペースアルゴリズムによって)得られたとします。両者を結びつける回転行列 $A$、すなわち $\hat s_i \approx A\,\hat r_i$ を満たす $A\in SO(3)$ を求める問題がWahbaの問題です。実際の観測にはノイズが乗るため、厳密な等式は成り立たず、次の重み付き二乗誤差を最小化する問題として定式化されます。

$$J(A) = \frac{1}{2}\sum_{i=1}^{N} w_i \left\| \hat{s}_i - A\, \hat{r}_i \right\|^2 \ \longrightarrow\ \min_{A \in SO(3)}$$

重み $w_i$ は各恒星の観測精度(明るさ・センタリング精度)に応じて大きくすることが多く、通常 $\sum_i w_i = 1$ に規格化します。この最小化問題を解いてみましょう。ノルムを展開すると

$$J(A) = \frac{1}{2}\sum_i w_i\left(\hat s_i\!\cdot\!\hat s_i - 2\hat s_i\!\cdot\!(A\hat r_i) + (A\hat r_i)\!\cdot\!(A\hat r_i)\right)$$

となりますが、$\hat s_i,\hat r_i$ は単位ベクトルであり、かつ $A$ は直交行列なので $(A\hat r_i)\cdot(A\hat r_i) = \hat r_i^{\mathsf T}A^{\mathsf T}A\,\hat r_i = \hat r_i \cdot \hat r_i = 1$ です。したがって $\hat s_i \cdot \hat s_i$ と $(A\hat r_i)\cdot(A\hat r_i)$ はどちらも $A$ によらない定数(いずれも1)であり、$J(A)$ を最小化することは

$$g(A) \equiv \sum_i w_i\, \hat s_i \cdot (A\hat r_i) = \mathrm{tr}\!\left(A\sum_i w_i\, \hat r_i \hat s_i^{\mathsf T}\right) = \mathrm{tr}(A B^{\mathsf T})$$

という利得関数(gain function) $g(A)$ を最大化することと等価です。ここで $B \equiv \sum_i w_i\, \hat s_i \hat r_i^{\mathsf T}$ を姿勢プロファイル行列と呼びます。つまりWahbaの問題は「$3\times3$行列 $B$ が与えられたとき、$\mathrm{tr}(AB^{\mathsf T})$ を最大にする回転行列 $A$ を見つけよ」という、より扱いやすい形に整理されました。

TRIAD法とQUEST法 — Wahbaの問題を実際に解く

TRIAD法(2ベクトル法) は、$N=2$ の場合に成り立つ最も単純な解法です。2つの恒星の観測ベクトル $\hat s_1,\hat s_2$(ボディ座標系)とカタログベクトル $\hat r_1,\hat r_2$(慣性座標系)から、それぞれの系で正規直交基底(トライアド)

$$\vec t_{1}^{\,b} = \hat s_1,\quad \vec t_{2}^{\,b} = \frac{\hat s_1 \times \hat s_2}{\lvert \hat s_1 \times \hat s_2\rvert},\quad \vec t_{3}^{\,b} = \vec t_1^{\,b} \times \vec t_2^{\,b}$$ $$\vec t_{1}^{\,r} = \hat r_1,\quad \vec t_{2}^{\,r} = \frac{\hat r_1 \times \hat r_2}{\lvert \hat r_1 \times \hat r_2\rvert},\quad \vec t_{3}^{\,r} = \vec t_1^{\,r} \times \vec t_2^{\,r}$$

を構成すると、姿勢行列は $A = \begin{bmatrix}\vec t_1^{\,b} & \vec t_2^{\,b} & \vec t_3^{\,b}\end{bmatrix}\begin{bmatrix}\vec t_1^{\,r} & \vec t_2^{\,r} & \vec t_3^{\,r}\end{bmatrix}^{\mathsf T}$ という単純な行列積で直接求まります。計算コストが極めて低い反面、3個以上の恒星情報を統計的に均すことができず、また第1ベクトル $\hat s_1$ を優先的に正確とみなす非対称な扱いになる(2番目のベクトルの誤差は $\vec t_2$ 以降にしか影響しないが1番目のベクトルの誤差は結果に直接乗る)という弱点があります。

3個以上の恒星をバランスよく使う厳密解が、Davenportのq法です。姿勢プロファイル行列 $B$ から、対称なトレースレス行列を含む $4\times4$ 行列

$$K = \begin{pmatrix} S - \sigma I & \vec z \\ \vec z^{\mathsf T} & \sigma \end{pmatrix}, \qquad S = B + B^{\mathsf T},\quad \sigma = \mathrm{tr}(B),\quad \vec z = \begin{pmatrix} B_{23}-B_{32} \\ B_{31}-B_{13} \\ B_{12}-B_{21} \end{pmatrix}$$

を作ると、$g(A)$ を最大化する最適クォータニオンは、この $K$ の最大固有値に対応する固有ベクトルとして得られることが示せます。すなわちWahbaの問題は、4次対称行列の固有値問題に帰着します。実運用の高速なオンボード計算では、この固有値問題を毎回律儀に解く代わりに、最大固有値が理想条件下( $\sum_i w_i \approx 1$ )で $\lambda_{\max}\approx 1$ に近いことを利用し、ニュートン・ラフソン法で固有値方程式(特性方程式)の根を数回の反復で高速に求めるQUESTアルゴリズム(QUaternion ESTimator)が広く使われています。QUESTはDavenportのq法と数学的に等価な解を、はるかに少ない計算量で与えるため、限られた計算資源しか持たない探査機の姿勢制御コンピュータでもリアルタイムに実行できます。

なお、観測される恒星の数が2個ぎりぎりの場合、視線方向(ボアサイト)を軸とした回転(ロール)についての姿勢決定精度は、ボアサイトに直交する2軸(チルト)方向の精度に比べて本質的に悪くなります。これは2つの恒星が張る平面の情報だけではボアサイトまわりの回転を弱くしか拘束できないためで、視野内でできるだけ角度的に離れた恒星を多数(実運用では数個〜十数個程度)使うことが、姿勢決定精度、とりわけロール方向の精度を向上させる鍵になります。

実務での使われ方

ハードウェア構成。 現代のスタートラッカーは、視野角およそ$10^\circ\times10^\circ$〜$20^\circ\times20^\circ$程度のレンズと、CCDまたはCMOSイメージセンサ(近年はCMOSが低消費電力・高耐放射線性の観点から主流)を組み合わせた光学ヘッドと、恒星カタログの照合・Wahbaの問題の求解・追尾処理を行う専用の処理ユニットから構成されます。オンボードに保持する恒星カタログは、Hipparcos/Tychoカタログのような数万〜十数万個規模の全天恒星カタログから、十分明るく(実用上は5〜6等級程度まで)、かつ位置精度と固有運動の情報が確かな恒星を数千個程度に絞り込んだサブセットが使われるのが一般的です。

精度と更新レート。 代表的な宇宙用スタートラッカー(たとえばSodern社のHYDRAシリーズ、Terma社のHE-5AS、Ball Aerospace社のCT-602、Blue Canyon TechnologiesのNSTなど)は、ボアサイトに直交する2軸方向でおよそ1〜10秒角(arcsec)、ボアサイトまわりのロール方向でその数倍〜10倍程度の精度を、1〜10Hz程度の更新レートで継続的に出力します。これは前回扱ったスタートラッカー単体の測角精度(サブピクセルのセントロイド精度)を、複数の恒星の統計的な平均化によってさらに高めた結果です。

ジャイロスコープとのセンサーフュージョン。 スタートラッカーには弱点もあります。第一に、機体が高速に姿勢変化(スルー)している最中は、露光時間中に星像が画像上を流れてしまい(ブラー)、恒星の検出・パターンマッチングが困難になります。第二に、太陽・月・地球からの迷光が視野に入ると、恒星が眩しさに埋もれて検出できなくなる「ブラインド」状態に陥ります。そこで実務では、スタートラッカーを単独で使うのではなく、慣性計測装置(IMU)に含まれるジャイロスコープと組み合わせたセンサーフュージョンが標準です。ジャイロは角速度 $\vec\omega$ を高いレート(数百Hz)で出力でき、これを前節のクォータニオン運動学方程式 $\dot q = \tfrac12\Omega(\vec\omega)q$ に沿って積分すれば、短時間であれば高精度に姿勢を伝播できます。しかしジャイロには避けられないバイアス誤差(ゼロ点ドリフト)があり、積分を続けるほど姿勢誤差が時間とともに蓄積していきます。そこで、姿勢変化が緩やかな定常フェーズではスタートラッカーの絶対姿勢測定値でジャイロのバイアスドリフトを定期的に補正し(カルマンフィルタによる姿勢決定・ジャイロバイアス推定の同時フィルタがよく用いられます)、急激な姿勢変化やブラインド期間中はジャイロの積分値だけで姿勢を維持する、という役割分担が組まれます。この構成は木星・土星探査機のような長期ミッションから、はやぶさ2のような小天体探査機まで、ほぼすべての現代の深宇宙探査機の姿勢決定・姿勢制御系(ADCS)に共通する基本アーキテクチャです。

演習問題

1. 3-2-1系列のオイラー角運動学方程式 $\dot\psi = (q\sin\phi+r\cos\phi)/\cos\theta$ において、ピッチ角 $\theta = 89^\circ$ のときと $\theta=89.99^\circ$ のときとで、$1/\cos\theta$ の値がそれぞれ何倍に増大するかを計算し、この係数がジンバルロックに向かってどのように発散していくかを数値で確認せよ。またクォータニオンの運動学方程式 $\dot q = \tfrac12\Omega(\vec\omega)q$ がこの種の発散を持たない理由を、この回で学んだ内容をもとに説明せよ。

2. 単位ベクトル $\hat n = (0,0,1)$(Z軸)のまわりに角度 $\Theta = 60^\circ$ だけ回転する姿勢を表すクォータニオン $q=(q_1,q_2,q_3,q_4)$ を求めよ。さらに本文の式 $A(q) = (q_4^2-|\vec q_v|^2)I + 2\vec q_v\vec q_v^{\mathsf T} - 2q_4[\vec q_v\times]$ を用いて対応する回転行列 $A$ を計算し、これが慣性系からボディ系への座標変換(本文の $\vec s = A\vec r$ の意味での姿勢行列)

$$A_z(\Theta) = \begin{pmatrix} \cos\Theta & \sin\Theta & 0 \\ -\sin\Theta & \cos\Theta & 0 \\ 0 & 0 & 1 \end{pmatrix}$$

と一致することを確認せよ。また $AA^{\mathsf T}=I$、$\det A = +1$ が成り立つことも確かめ、これが正しい回転行列であることを検証せよ。

3. あるスタートラッカーが2つの恒星を検出し、ボディ座標系での視線方向単位ベクトルが $\hat s_1 = (1,0,0)$、$\hat s_2 = (0, \tfrac{\sqrt2}{2}, \tfrac{\sqrt2}{2})$、対応する恒星カタログ上の慣性座標系での方向が $\hat r_1 = (0,1,0)$、$\hat r_2 = (-\tfrac{\sqrt2}{2},0,\tfrac{\sqrt2}{2})$ であったとする。TRIAD法の手順に従ってボディ座標系・慣性座標系それぞれのトライアド $\{\vec t_1,\vec t_2,\vec t_3\}$ を構成し、姿勢行列 $A$ を求めよ。

4. 姿勢決定にスタートラッカーとジャイロを併用しているある探査機で、太陽光の迷光により $\Delta t = 5$ 分間スタートラッカーがブラインド状態になり、その間はジャイロの積分値のみで姿勢を維持したとする。ジャイロのバイアスドリフト率が $\dot\epsilon = 0.05\ \text{deg/hr}$ であるとき、このブラインド期間中に蓄積するおおよその姿勢誤差を秒角(arcsec)単位で見積もれ。またこの誤差の大きさを、本文で挙げたスタートラッカー単体の代表的な姿勢決定精度(1〜10秒角)と比較し、なぜ長時間のブラインド状態が姿勢制御システムにとって望ましくないのかを論じよ。

まとめと次回予告

スタートラッカーは、前回の光学航法が近傍天体の視差から探査機の「位置」を求めたのとは対をなす技術で、恒星という無限遠の基準点への視線方向から探査機の「向き(姿勢)」を求めます。回転行列は姿勢を厳密に表現できますが冗長であり、オイラー角は直感的な反面ジンバルロックという特異点を抱え、クォータニオンはその両方の弱点を回避する形で宇宙工学の標準的な姿勢表現として定着しています。そして複数の恒星観測から最適な姿勢を求めるWahbaの問題は、TRIAD法による高速な近似解から、Davenportのq法・QUESTアルゴリズムによる統計的に最適な解まで、実務で使われる具体的なアルゴリズムへと発展してきました。この姿勢決定という土台があって初めて、探査機はアンテナを地球に、太陽電池を太陽に、そしてカメラを目標天体に正確に向けることができるのです。

次回は、位置と姿勢のどちらでもなく、地上局や天体への方向(角度)だけを手がかりに探査機の軌道を推定する**角度のみによる航法(Bearings-only navigation)**に軽く触れます。距離の情報が得られない、あるいは信頼できない状況で、方向の変化だけからどこまで軌道の情報を引き出せるのかを見ていきます。

参考文献

• G. Wahba, “A Least Squares Estimate of Satellite Attitude,” SIAM Review, Vol. 7, No. 3, 1965
• M. D. Shuster, S. D. Oh, “Three-Axis Attitude Determination from Vector Observations,” Journal of Guidance and Control, Vol. 4, No. 1, 1981
• M. D. Shuster, “Approximate Algorithms for Fast Optimal Attitude Computation,” AIAA Guidance and Control Conference, 1978
• F. L. Markley, J. L. Crassidis, Fundamentals of Spacecraft Attitude Determination and Control, Springer, 2014
• J. R. Wertz (ed.), Spacecraft Attitude Determination and Control, Kluwer Academic Publishers
• C. C. Liebe, “Accuracy Performance of Star Trackers — A Tutorial,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 38, No. 2