夜空を見上げると、無数の星が独特のパターンを描いています。オリオン座のベルトの3つ星、北斗七星の柄杓の形 — 人間は古来、これらの星座パターンを認識して方角を知り、大海原を航海してきました。実は人工衛星も全く同じことをしています。宇宙空間に浮かぶカメラで恒星を撮影し、その配置パターンをデータベースと照合することで、自分がどちらを向いているかを驚くほど高い精度で決定するのです。この装置がスタートラッカー（Star Tracker）です。
スタートラッカーの精度は、数アーク秒（$1$ アーク秒 $= 1/3600°$）に達します。これは約 $20$ km先の1円玉を見分ける角度分解能に相当し、太陽センサや地磁気センサとは桁違いの精度です。地球観測衛星が地表のわずか数メートルの解像度で撮像できるのも、通信衛星がピンポイントでビームを地上局に向けられるのも、スタートラッカーによる高精度姿勢決定があってこそです。
スタートラッカーを理解すると、以下の技術領域への理解が深まります。
- 高精度地球観測: 光学センサの指向精度の根幹をなす姿勢決定技術
- 天文観測衛星: ハッブル宇宙望遠鏡やJWSTなど、サブアーク秒の姿勢安定性が要求されるミッション
- 画像処理とパターン認識: 星像の重心検出やパターンマッチングは、コンピュータビジョンの基本技術と共通する
- 姿勢決定アルゴリズム: スタートラッカーの出力は、TRIAD法やQUEST法などの姿勢決定アルゴリズムへの主要な入力となる
本記事の内容
- スタートラッカーの基本構成と動作フロー
- 星像の撮像とサブピクセル重心検出
- 星カタログとパターンマッチングの原理
- トライアングルマッチング法の詳細
- Lost-in-Space問題とその解法
- 運用上の制約（太陽回避角、月、SAA）
- Pythonでの星パターンマッチングの簡易実装と精度解析
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
スタートラッカーの基本構成
ハードウェア構成
スタートラッカーは、一言で言えば「宇宙用の高感度カメラ + 画像認識コンピュータ」です。街中の防犯カメラが人の顔を認識するのと本質的に同じ技術が、宇宙では恒星の配置パターンを認識するために使われています。
基本的なハードウェア構成は以下のとおりです。
光学系: 視野角（FOV）は一般に $8° \times 8°$ 〜 $20° \times 20°$ 程度の広角レンズです。焦点距離 $f$ は数十mmが一般的です。視野が広いほど一度に多くの星が写りますが、各星像の位置精度が低下するトレードオフがあります。
イメージセンサ: CCD（Charge-Coupled Device）またはCMOS（Active Pixel Sensor）が使用されます。近年はCMOSセンサが主流で、解像度は $512 \times 512$ 〜 $2048 \times 2048$ ピクセル程度です。宇宙環境での放射線耐性が重要な設計要件であり、放射線によるホットピクセルの蓄積が長期ミッションでの課題となります。
バッフル: センサへの迷光（太陽光や地球反射光）の入射を防ぐ遮光筒です。一般に、太陽回避角（太陽からの最小離角）は $20°$ 〜 $45°$ に設定されます。
プロセッサ: 星像の検出、重心計算、パターンマッチング、姿勢算出を行う搭載コンピュータです。リアルタイム処理が求められるため、DSPやFPGAが用いられることが多いです。
検出可能な等級
スタートラッカーが検出する恒星の明るさは、天文学の等級（magnitude）で表されます。等級が小さいほど明るく、1等級の差は明るさの約2.512倍（$10^{0.4}$ 倍）に相当します。
$$ \frac{F_1}{F_2} = 10^{0.4(m_2 – m_1)} $$
ここで $F$ はフラックス（光束）、$m$ は見かけの等級です。
典型的なスタートラッカーは、5.5等級〜7等級程度までの恒星を検出します。夜空で肉眼で見える星がおよそ6等級までなので、スタートラッカーは人間の目と同等かやや暗い星まで捉えることができます。全天で5.5等級までの星は約9,000個、6.5等級までだと約50,000個あります。視野角 $16° \times 16°$ のセンサの場合、一般にどの方向を向いても5〜30個程度の星が視野内に入ります。
星の数が多すぎても少なすぎても問題です。多すぎるとパターンマッチングの計算量が爆発し、少なすぎると同定に必要な情報が不足します。適切なしきい値等級の設定が重要な設計パラメータとなります。
スタートラッカーのハードウェア構成がわかったところで、次に星像の取得から姿勢算出に至る処理フローを詳しく見ていきます。
処理フロー: 撮像から姿勢算出まで
概要
スタートラッカーの処理は、大きく4つのステップに分かれます。これは「写真を撮って → 星を見つけて → どの星か特定して → 姿勢を計算する」という流れです。
ステップ1: 撮像（Image Acquisition） → ステップ2: 星像検出・重心計算（Centroiding） → ステップ3: 星同定（Star Identification） → ステップ4: 姿勢算出（Attitude Determination）
各ステップを順に解説していきます。
ステップ1: 撮像
イメージセンサが露光し、星像を取得します。露光時間は衛星の角速度に依存します。角速度が大きいと、露光中に星像がぶれて（smearing）点像から線像になり、重心精度が劣化します。
角速度 $\omega$ [rad/s]、露光時間 $t_{\text{exp}}$ [s]、焦点距離 $f$ [mm]、ピクセルサイズ $p$ [μm] のとき、星像のぶれ量（ピクセル数）は次のように見積もれます。
$$ \Delta_{\text{smear}} = \frac{\omega \cdot t_{\text{exp}} \cdot f}{p} $$
たとえば、$\omega = 0.1$ °/s $= 1.745 \times 10^{-3}$ rad/s、$t_{\text{exp}} = 0.1$ s、$f = 50$ mm、$p = 12$ μm のとき、
$$ \Delta_{\text{smear}} = \frac{1.745 \times 10^{-3} \times 0.1 \times 50}{0.012} \approx 0.73 \text{ pixel} $$
1ピクセル未満のぶれであれば許容範囲ですが、これを超えるとデスメアリング補正が必要になります。
ステップ2: 星像検出とサブピクセル重心計算
撮像画像から星像を検出し、その中心位置をサブピクセル精度で求める処理です。この精度がスタートラッカー全体の性能の土台となります。
星像の検出: まず、画像にしきい値処理を適用し、背景ノイズより有意に明るいピクセル群（星像候補）を抽出します。連結成分ラベリングにより、各星像候補をクラスタリングします。
サブピクセル重心計算: 星像の光学的な広がり（PSF: Point Spread Function）は通常ガウシアン形状に近似できるため、ピクセル値の重み付き平均により、ピクセルの整数座標よりも細かい精度で中心位置を求めることができます。
$N$ 個のピクセルからなる星像クラスタの重心 $(\bar{x}, \bar{y})$ は次のように計算されます。
$$ \begin{align} \bar{x} &= \frac{\sum_{k=1}^{N} (I_k – I_{\text{bg}}) \cdot x_k}{\sum_{k=1}^{N} (I_k – I_{\text{bg}})} \\ \bar{y} &= \frac{\sum_{k=1}^{N} (I_k – I_{\text{bg}}) \cdot y_k}{\sum_{k=1}^{N} (I_k – I_{\text{bg}})} \end{align} $$
ここで $I_k$ はピクセル $k$ の強度値、$I_{\text{bg}}$ は背景ノイズレベル、$(x_k, y_k)$ はピクセル $k$ の座標です。背景を差し引くことで、暗いピクセルの寄与を抑制し、重心精度を向上させます。
この重み付き重心法で達成可能な精度は、典型的に $0.05$ 〜 $0.2$ ピクセルです。さらに高い精度が必要な場合は、ガウシアンフィッティングにより $0.01$ ピクセル以下の精度も実現できます。
ステップ3からステップ4への流れ
重心計算で得られたイメージ座標 $(\bar{x}, \bar{y})$ は、光学系のパラメータ（焦点距離、歪曲収差）を用いて、機体座標系での方向ベクトルに変換されます。イメージセンサの中心を原点として、焦点距離 $f$ の位置にセンサがあるピンホールカメラモデルでは、星 $i$ の方向ベクトルは次のとおりです。
$$ \bm{b}_i = \frac{1}{\sqrt{\bar{x}_i^2 + \bar{y}_i^2 + f^2}} \begin{pmatrix} \bar{x}_i \\ \bar{y}_i \\ f \end{pmatrix} $$
ここで $\bm{b}_i$ は機体座標系での単位方向ベクトルです。この方向ベクトル群 $\{\bm{b}_1, \bm{b}_2, \ldots\}$ と星カタログから得られる慣性座標系での方向ベクトル群 $\{\bm{r}_1, \bm{r}_2, \ldots\}$ を対応付け、回転行列を求めるのがステップ3（星同定）とステップ4（姿勢算出）の仕事です。
ここまでで処理フローの全体像を把握しました。次に、最も計算集約的で知的なステップである「星同定」— すなわちパターンマッチングの手法を詳しく見ていきましょう。
星同定: パターンマッチングの原理
星同定問題とは
星同定の問題は次のように定式化されます。「イメージセンサ上で検出された $n$ 個の星像に対応する星カタログ中の恒星を特定せよ」。
これは本質的にパターン認識問題であり、いくつかの困難を伴います。
- 衛星の姿勢が未知: 当然ながら姿勢を求めるための処理なので、事前に姿勢がわからない
- 検出されない星がある: 感度限界以下の暗い星はカタログにあっても検出されない
- 偽星像（spurious）が存在する: 宇宙線のヒット、ホットピクセル、迷光などにより偽の星像が生じる
- 位置精度の限界: 重心計算のノイズにより、測定された方向ベクトルには誤差がある
これらの困難を乗り越えて頑健に星を同定するために、さまざまなアルゴリズムが開発されてきました。
星間角距離の不変性
星同定アルゴリズムの多くは、一つの重要な幾何学的性質に基づいています。それは「2つの星の間の角距離は、衛星の姿勢に依存しない」という事実です。
星 $i$ と星 $j$ の角距離 $\theta_{ij}$ は、内積で計算できます。
$$ \cos\theta_{ij} = \bm{b}_i \cdot \bm{b}_j $$
ここで $\bm{b}_i, \bm{b}_j$ は機体座標系での方向ベクトルです。重要なのは、慣性座標系での方向ベクトル $\bm{r}_i, \bm{r}_j$ からも同じ角距離が得られることです。
$$ \cos\theta_{ij} = \bm{r}_i \cdot \bm{r}_j $$
これは回転行列 $\bm{A}$ が内積を保存するため（$\bm{A}$ は直交行列）です。
$$ \bm{b}_i \cdot \bm{b}_j = (\bm{A}\bm{r}_i)^T(\bm{A}\bm{r}_j) = \bm{r}_i^T \bm{A}^T \bm{A} \bm{r}_j = \bm{r}_i \cdot \bm{r}_j $$
この不変量を利用すれば、姿勢が未知でも、検出された星間の角距離パターンをカタログの星間角距離と照合することで星を同定できるのです。
トライアングルマッチング法
トライアングルマッチング法は、星同定の古典的かつ直感的なアルゴリズムです。3つの星が作る三角形の「形」は姿勢に依存しないので、この形をカタログの三角形と照合します。
アルゴリズムの手順は以下のとおりです。
事前準備（オフライン）:
- 星カタログからセンサの検出限界等級以下の星を選び、カタログ星リストを作成
- カタログ星リストの全ペアについて角距離を計算し、角距離テーブルを構築
- 必要に応じて、全三つ組の三角形特徴量を計算し、三角形データベースを構築
星同定（オンライン）:
- 検出された星像から3つの星（三つ組）を選ぶ
- 3つの星間角距離 $(\theta_{12}, \theta_{13}, \theta_{23})$ を計算する
- これらの角距離をソートして正規化し、三角形の「特徴量」とする
- 三角形データベースを検索し、特徴量が一致する（許容誤差内の）カタログ三角形を見つける
- 一致が一意であれば、3つの星の同定が完了
- 残りの星は、同定済みの星との角距離を用いて順次同定する
三角形の特徴量としてよく使われるのは、3つの角距離をソートした組 $(d_1, d_2, d_3)$（$d_1 \leq d_2 \leq d_3$）です。あるいは、三角形の3つの内角 $(\alpha, \beta, \gamma)$ も姿勢不変な特徴量であり、角距離の絶対値に依存しないためスケール不変という利点があります。
Lost-in-Space問題
スタートラッカーの星同定には2つの動作モードがあります。
トラッキングモード（Tracking Mode）: 前の時刻の姿勢推定値がわかっている場合、次の時刻でどの星が視野内にいるかを予測できます。予測位置の近傍で星を探すだけなので、計算量は少なく高速です。通常運用ではこのモードが使われます。
ロストインスペースモード（Lost-in-Space Mode）: 衛星の姿勢に関する事前情報が全くない場合のモードです。打上げ直後、安全モードからの復帰、長時間の通信途絶後などに発生します。「宇宙で迷子になった」状態から、星像だけで姿勢を復元しなければなりません。
Lost-in-Space問題はトラッキングモードと比較して計算量が格段に大きくなります。カタログに $N$ 個の星がある場合、全ペアの角距離は $\binom{N}{2} \approx N^2/2$ 個、全三つ組は $\binom{N}{3} \approx N^3/6$ 個にもなるため、効率的な検索アルゴリズムとデータ構造（k-d木やハッシュテーブル）が不可欠です。
近年は、幾何学的ハッシング法やpyramid法など、Lost-in-Space問題をより高速に解くアルゴリズムが開発されています。これらは、角距離の組をハッシュ値に変換して $O(1)$ の検索を実現するものです。
星同定のアルゴリズムを理解したところで、次にスタートラッカーの精度と、その精度に影響を与える要因を分析します。
スタートラッカーの精度解析
精度の構成要素
スタートラッカーの姿勢決定精度は、複数の誤差源の寄与が合成されて決まります。主要な誤差源は以下のとおりです。
1. 重心検出誤差: ピクセルノイズ（読み出しノイズ、暗電流ノイズ、光子ショットノイズ）に起因する星像重心の不確定性です。これがスタートラッカーの精度の根本的な限界を決めます。
重心検出のS/N比を $\text{SNR}$ とすると、1つの星像の重心精度 $\sigma_c$（ピクセル単位）は近似的に次のように表されます。
$$ \sigma_c \approx \frac{\sigma_{\text{PSF}}}{\text{SNR}} $$
ここで $\sigma_{\text{PSF}}$ はPSF（点像分布関数）のガウシアン幅（ピクセル単位）です。
2. 光学歪曲収差: レンズの収差により、像面上の星の位置が理想的なピンホールモデルから系統的にずれます。この誤差は校正データにより補正可能ですが、残差が精度を制限します。
3. 星カタログの精度: カタログ中の恒星の位置精度は、現代のカタログ（Hipparcos: $\sim 1$ mas、Gaia: $\sim 0.02$ mas）では非常に高く、通常は支配的な誤差源にはなりません。
4. 固有運動: 恒星はわずかながら固有運動を持っています。数年スケールでは無視できる場合が多いですが、長期ミッションでは補正が必要です。
姿勢精度の見積もり
$n$ 個の星が検出され、各星の方向ベクトルの精度が $\sigma_{\text{star}}$ [rad] であるとき、3軸の姿勢精度は次のように見積もれます。
ボアサイト周りの回転（ロール）とボアサイトに垂直な2軸の回転（ピッチ、ヨー）では精度が異なります。ボアサイトに垂直な方向の精度は次のとおりです。
$$ \begin{equation} \sigma_{\text{cross}} \approx \frac{\sigma_{\text{star}}}{\sqrt{n}} \end{equation} $$
これは $n$ 個の独立な測定の統計的平均化によるものです。一方、ボアサイト周りの回転精度は次のようになります。
$$ \begin{equation} \sigma_{\text{roll}} \approx \frac{\sigma_{\text{star}}}{\sin\theta_{\text{rms}} \cdot \sqrt{n}} \end{equation} $$
ここで $\theta_{\text{rms}}$ は検出された星群のボアサイトからの平均離角です。$\sin\theta_{\text{rms}}$ がボアサイト周りの回転に対するレバーアームに相当し、視野の端に星が分布しているほどロール精度が向上します。逆に、すべての星がボアサイト近傍に集中していると $\sin\theta_{\text{rms}} \approx 0$ となり、ロール精度が極端に劣化します。
典型的な値として、$\sigma_{\text{star}} = 5$ arcsec、$n = 15$ 個、$\theta_{\text{rms}} = 5°$ の場合、
$$ \sigma_{\text{cross}} \approx \frac{5}{\sqrt{15}} \approx 1.3 \text{ arcsec} $$
$$ \sigma_{\text{roll}} \approx \frac{5}{0.087 \times \sqrt{15}} \approx 14.8 \text{ arcsec} $$
ボアサイト周りの精度がクロス軸に比べて約10倍悪いことがわかります。これはスタートラッカー固有の特性であり、高いロール精度が必要な場合は、ボアサイト方向が異なる2台以上のスタートラッカーを搭載する設計がとられます。
精度の理論的な枠組みを理解したところで、次にスタートラッカーの運用上の制約を見ていきましょう。
運用上の制約
太陽回避角
スタートラッカーにとって最大の「敵」は太陽光です。太陽は恒星と比較して圧倒的に明るく（見かけの等級: $-26.7$ 等）、直接光がセンサに入射するとCCD/CMOSが飽和し、場合によっては素子が損傷する可能性もあります。
そのため、スタートラッカーは太陽から一定の角度以内を避けて運用する必要があります。この角度を太陽回避角（Sun Exclusion Angle）と呼び、典型的には $20°$ 〜 $45°$ です。バッフルの設計により、この角度はある程度制御できます。
太陽回避角 $\theta_{\text{excl}}$ がもたらす影響は、全天のうちスタートラッカーが使用できない領域の立体角として定量化できます。
$$ \Omega_{\text{excl}} = 2\pi(1 – \cos\theta_{\text{excl}}) $$
全天の立体角は $4\pi$ なので、使用不可の割合は次のとおりです。
$$ f_{\text{excl}} = \frac{\Omega_{\text{excl}}}{4\pi} = \frac{1 – \cos\theta_{\text{excl}}}{2} $$
$\theta_{\text{excl}} = 30°$ の場合、$f_{\text{excl}} = (1 – \cos 30°)/2 \approx 6.7\%$ であり、全天の $93.3\%$ の方向でスタートラッカーを使用できます。
月と地球の影響
太陽だけでなく、月と地球もスタートラッカーに影響を与えます。
月: 満月時の見かけの等級は約 $-12.7$ 等です。太陽ほどではありませんが、月が視野内に入ると広範囲のピクセルが飽和し、星の検出が困難になります。月回避角は一般に $10°$ 〜 $20°$ に設定されます。
地球: 低軌道衛星の場合、地球は全天の大きな部分を占めます。地球の大気は太陽光を散乱し（地球反射光: アルベド）、地平線付近で明るい背景光を生じます。また、地球の縁（リム）からの散乱光はバッフルに入射して迷光の原因となります。地球回避角は $10°$ 〜 $30°$ 程度です。
南大西洋異常帯（SAA）
南大西洋異常帯（South Atlantic Anomaly）は、バンアレン帯の内側の放射線帯が地表に近づいている領域です。南米大陸の東方、大西洋上空にあり、低軌道衛星がこの領域を通過するとき、高エネルギー粒子のフラックスが通常の100倍〜1000倍になります。
これらの高エネルギー粒子がイメージセンサに衝突すると、局所的な高輝度ピクセル（宇宙線ヒット）が生じ、偽の星像として検出される可能性があります。SAA通過中はこのような偽星像の発生率が著しく増加するため、星同定アルゴリズムにはロバストな偽星像除去機能が必要です。
対策としては、以下の方法が一般的です。
- 時間フィルタリング: 宇宙線ヒットは1フレームのみに現れるが、恒星は複数フレームにわたって同じ位置に存在する
- 形状フィルタリング: 宇宙線ヒットはPSFと形状が異なる（シャープすぎるか、長い軌跡を残す）
- SAA領域での検出しきい値の引き上げ: 偽検出を減らす代わりに暗い星を犠牲にする
以上の運用制約を踏まえた上で、次にPythonでトライアングルマッチング法の簡易実装を行い、星同定と姿勢精度の関係を数値的に確認しましょう。
Pythonでの実装: トライアングルマッチング法
簡易星カタログの構築
まず、シミュレーション用の簡易星カタログと、星の観測をシミュレートする関数を実装します。実際の Hipparcos カタログの代わりに、球面上にランダムに分布する星を生成します。
import numpy as np
from itertools import combinations
def generate_star_catalog(n_stars=500, seed=42):
"""球面上にランダムな星カタログを生成（慣性座標系）"""
rng = np.random.default_rng(seed)
# 球面上の一様分布
phi = rng.uniform(0, 2 * np.pi, n_stars)
cos_theta = rng.uniform(-1, 1, n_stars)
sin_theta = np.sqrt(1 - cos_theta**2)
catalog = np.column_stack([
sin_theta * np.cos(phi),
sin_theta * np.sin(phi),
cos_theta
])
# 等級を割り当て（明るい星ほど少ない）
magnitudes = rng.exponential(1.5, n_stars) + 1.0
magnitudes = np.clip(magnitudes, 1.0, 7.0)
return catalog, magnitudes
def observe_stars(catalog, magnitudes, attitude_matrix, fov_half_angle_deg=10.0,
mag_limit=6.0, noise_sigma_arcsec=5.0, rng=None):
"""指定姿勢でのスタートラッカー観測をシミュレート"""
if rng is None:
rng = np.random.default_rng(0)
fov_half = np.radians(fov_half_angle_deg)
# 機体座標系に変換
body_vectors = (attitude_matrix @ catalog.T).T
# ボアサイト（z軸）からの角距離
cos_angles = body_vectors[:, 2]
# FOV内かつ検出限界以上の星を選択
in_fov = cos_angles > np.cos(fov_half)
bright_enough = magnitudes < mag_limit
detected = in_fov & bright_enough
# 検出された星のインデックスとベクトル
det_indices = np.where(detected)[0]
det_vectors = body_vectors[detected].copy()
# 測定ノイズを追加
noise_rad = noise_sigma_arcsec / 3600 * np.pi / 180
for i in range(len(det_vectors)):
# 法線方向の微小ランダム回転
noise = rng.normal(0, noise_rad, 3)
det_vectors[i] += noise
det_vectors[i] /= np.linalg.norm(det_vectors[i])
return det_indices, det_vectors
# カタログ生成
catalog, magnitudes = generate_star_catalog(n_stars=2000, seed=42)
print(f"カタログ星数: {len(catalog)}")
print(f"6等級以下の星数: {np.sum(magnitudes < 6.0)}")
このコードでは、2000個の恒星を球面上に一様分布させ、等級（明るさ）を指数分布で割り当てています。実際のカタログでは明るい星ほど数が少なく、暗い星が圧倒的に多いことを反映しています。
角距離テーブルの構築とトライアングルマッチング
次に、カタログ星間の角距離テーブルを構築し、トライアングルマッチング法で星同定を行います。
import numpy as np
from itertools import combinations
def build_angular_distance_table(catalog, magnitudes, mag_limit=6.0,
max_angle_deg=25.0):
"""カタログ星ペアの角距離テーブルを構築"""
bright = np.where(magnitudes < mag_limit)[0]
max_angle = np.radians(max_angle_deg)
pairs = []
distances = []
for i, j in combinations(bright, 2):
cos_d = np.dot(catalog[i], catalog[j])
cos_d = np.clip(cos_d, -1, 1)
d = np.arccos(cos_d)
if d < max_angle:
pairs.append((i, j))
distances.append(d)
return np.array(pairs), np.array(distances)
def triangle_match(obs_vectors, obs_indices_true, catalog, pairs, pair_distances,
tolerance_arcsec=30.0):
"""トライアングルマッチングによる星同定"""
tolerance = tolerance_arcsec / 3600 * np.pi / 180
n_obs = len(obs_vectors)
if n_obs < 3:
return {}
# 観測星の角距離を計算
obs_pairs = []
obs_dists = []
for i, j in combinations(range(n_obs), 2):
cos_d = np.clip(np.dot(obs_vectors[i], obs_vectors[j]), -1, 1)
obs_pairs.append((i, j))
obs_dists.append(np.arccos(cos_d))
# 最初の3つの星で三角形を構成
tri_obs = list(combinations(range(min(n_obs, 5)), 3))
best_match = {}
best_count = 0
for tri in tri_obs:
i, j, k = tri
# 3つの角距離を計算
d_ij = np.arccos(np.clip(np.dot(obs_vectors[i], obs_vectors[j]), -1, 1))
d_ik = np.arccos(np.clip(np.dot(obs_vectors[i], obs_vectors[k]), -1, 1))
d_jk = np.arccos(np.clip(np.dot(obs_vectors[j], obs_vectors[k]), -1, 1))
obs_triangle = sorted([d_ij, d_ik, d_jk])
# カタログの三角形を検索
# ペアテーブルからij角距離に一致するペアを探す
candidates_ij = np.where(np.abs(pair_distances - d_ij) < tolerance)[0]
for c_ij in candidates_ij:
cat_i, cat_j = pairs[c_ij]
# cat_i, cat_j が決まった状態で cat_k を探す
# cat_k は cat_i との角距離が d_ik に近く、
# かつ cat_j との角距離が d_jk に近い星
mask_ik = ((pairs[:, 0] == cat_i) | (pairs[:, 1] == cat_i))
candidates_ik = np.where(mask_ik & (np.abs(pair_distances - d_ik) < tolerance))[0]
for c_ik in candidates_ik:
p = pairs[c_ik]
cat_k = p[1] if p[0] == cat_i else p[0]
if cat_k == cat_j:
continue
# cat_j - cat_k の角距離を確認
cos_jk = np.clip(np.dot(catalog[cat_j], catalog[cat_k]), -1, 1)
d_jk_cat = np.arccos(cos_jk)
if np.abs(d_jk_cat - d_jk) < tolerance:
# 三角形が一致
match = {i: cat_i, j: cat_j, k: cat_k}
# 残りの星を順次同定
for m in range(n_obs):
if m in match:
continue
for ref_obs, ref_cat in match.items():
cos_d = np.clip(
np.dot(obs_vectors[m], obs_vectors[ref_obs]), -1, 1
)
d_obs = np.arccos(cos_d)
# カタログから距離が一致する星を探す
mask = ((pairs[:, 0] == ref_cat) | (pairs[:, 1] == ref_cat))
cands = np.where(
mask & (np.abs(pair_distances - d_obs) < tolerance)
)[0]
for c in cands:
p = pairs[c]
candidate = p[1] if p[0] == ref_cat else p[0]
if candidate not in match.values():
match[m] = candidate
break
if m in match:
break
if len(match) > best_count:
best_count = len(match)
best_match = match.copy()
return best_match
# テスト: ランダム姿勢での星同定
def random_rotation_matrix(rng):
"""ランダムな回転行列を生成"""
q = rng.normal(0, 1, 4)
q /= np.linalg.norm(q)
q0, q1, q2, q3 = q
return np.array([
[1-2*(q2**2+q3**2), 2*(q1*q2-q0*q3), 2*(q1*q3+q0*q2)],
[2*(q1*q2+q0*q3), 1-2*(q1**2+q3**2), 2*(q2*q3-q0*q1)],
[2*(q1*q3-q0*q2), 2*(q2*q3+q0*q1), 1-2*(q1**2+q2**2)]
])
print("角距離テーブルを構築中...")
pairs, pair_distances = build_angular_distance_table(
catalog, magnitudes, mag_limit=6.0, max_angle_deg=25.0
)
print(f"ペア数: {len(pairs)}")
rng = np.random.default_rng(123)
A_true = random_rotation_matrix(rng)
det_indices, det_vectors = observe_stars(
catalog, magnitudes, A_true,
fov_half_angle_deg=10.0, mag_limit=5.5,
noise_sigma_arcsec=5.0, rng=rng
)
print(f"\n検出された星の数: {len(det_indices)}")
print("トライアングルマッチング実行中...")
match_result = triangle_match(
det_vectors, det_indices, catalog, pairs, pair_distances,
tolerance_arcsec=30.0
)
# 結果の評価
correct = 0
for obs_idx, cat_idx in match_result.items():
if det_indices[obs_idx] == cat_idx:
correct += 1
print(f"同定された星: {len(match_result)} / {len(det_indices)}")
print(f"正しい同定: {correct} / {len(match_result)}")
if len(match_result) > 0:
print(f"正答率: {correct/len(match_result)*100:.1f}%")
このコードでは、まず全カタログ星ペアの角距離テーブルを構築し、次にランダムな姿勢で観測をシミュレートしてトライアングルマッチングを実行しています。結果として、検出された星の大半が正しく同定されることが確認できます。30 arcsec のマッチング許容誤差は、5 arcsec の測定ノイズに対して十分なマージンを持たせた値です。
姿勢精度のモンテカルロ解析
最後に、スタートラッカーの姿勢決定精度をモンテカルロシミュレーションで評価します。星が同定された後の姿勢算出には、ここでは簡易的にSVDベースの手法を使用します（詳細なTRIAD法とQUEST法は次の記事で解説します）。
import numpy as np
import matplotlib.pyplot as plt
def estimate_attitude_svd(body_vectors, ref_vectors):
"""SVDによる姿勢推定（Wahba問題の解）"""
n = len(body_vectors)
B = np.zeros((3, 3))
for i in range(n):
B += np.outer(body_vectors[i], ref_vectors[i])
U, S, Vt = np.linalg.svd(B)
det_sign = np.linalg.det(U) * np.linalg.det(Vt.T)
D = np.diag([1, 1, det_sign])
A_est = U @ D @ Vt
return A_est
def attitude_error_deg(A_true, A_est):
"""回転行列の誤差角度 [deg]"""
R_err = A_est @ A_true.T
cos_angle = np.clip((np.trace(R_err) - 1) / 2, -1, 1)
return np.degrees(np.arccos(cos_angle))
# --- モンテカルロシミュレーション ---
n_trials = 300
noise_levels = [1.0, 3.0, 5.0, 10.0, 20.0] # arcsec
results = {sigma: [] for sigma in noise_levels}
rng = np.random.default_rng(42)
for trial in range(n_trials):
A_true = random_rotation_matrix(rng)
for sigma in noise_levels:
det_indices, det_vectors = observe_stars(
catalog, magnitudes, A_true,
fov_half_angle_deg=10.0, mag_limit=5.5,
noise_sigma_arcsec=sigma, rng=rng
)
if len(det_indices) < 3:
continue
# 真のカタログベクトル（ノイズなし）
ref_vectors = catalog[det_indices]
# SVDで姿勢推定
A_est = estimate_attitude_svd(det_vectors, ref_vectors)
# 誤差計算
err = attitude_error_deg(A_true, A_est)
results[sigma].append(err)
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: ノイズレベルごとの姿勢誤差分布
colors = ['#00bcd4', '#4caf50', '#ff9800', '#e91e63', '#9c27b0']
for sigma, color in zip(noise_levels, colors):
errs = np.array(results[sigma])
errs_arcsec = errs * 3600
axes[0].hist(errs_arcsec, bins=30, alpha=0.5, color=color,
label=f'$\\sigma$ = {sigma}" (mean: {np.mean(errs_arcsec):.1f}")')
axes[0].set_xlabel('Attitude error [arcsec]', fontsize=12)
axes[0].set_ylabel('Count', fontsize=12)
axes[0].set_title('Attitude Determination Error Distribution', fontsize=13)
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)
axes[0].set_facecolor('#1a1a2e')
# 右: ノイズ vs 平均姿勢誤差
mean_errors = [np.mean(results[s]) * 3600 for s in noise_levels]
std_errors = [np.std(results[s]) * 3600 for s in noise_levels]
axes[1].errorbar(noise_levels, mean_errors, yerr=std_errors,
fmt='o-', color='#00bcd4', linewidth=2, markersize=8,
capsize=5, capthick=2, label='Monte Carlo result')
# 理論値: σ_att ≈ σ_star / sqrt(n), n ≈ 15
n_stars_avg = 15
theory_errors = [s / np.sqrt(n_stars_avg) for s in noise_levels]
axes[1].plot(noise_levels, theory_errors, '--', color='#ff9800',
linewidth=2, label=f'Theory: $\\sigma / \\sqrt{{{n_stars_avg}}}$')
axes[1].set_xlabel('Star measurement noise $\\sigma$ [arcsec]', fontsize=12)
axes[1].set_ylabel('Mean attitude error [arcsec]', fontsize=12)
axes[1].set_title('Attitude Error vs Measurement Noise', fontsize=13)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
axes[1].set_facecolor('#1a1a2e')
fig.patch.set_facecolor('#0f0f23')
for ax in axes:
ax.tick_params(colors='white')
ax.xaxis.label.set_color('white')
ax.yaxis.label.set_color('white')
ax.title.set_color('white')
for spine in ax.spines.values():
spine.set_color('white')
plt.tight_layout()
plt.savefig('star_tracker_accuracy.png', dpi=150, bbox_inches='tight',
facecolor='#0f0f23')
plt.show()
# 数値結果の表示
print("\n=== 姿勢決定精度 (Monte Carlo) ===")
print(f"{'ノイズσ [arcsec]':>20} {'平均誤差 [arcsec]':>20} {'3σ誤差 [arcsec]':>20}")
for sigma in noise_levels:
errs = np.array(results[sigma]) * 3600
print(f"{sigma:>20.1f} {np.mean(errs):>20.1f} {3*np.std(errs):>20.1f}")
左のヒストグラムから、測定ノイズ $\sigma$ が大きくなるにつれて姿勢決定誤差の分布が広がり、平均値も増大する傾向が明確に見て取れます。$\sigma = 5$ arcsec（典型的なスタートラッカーのノイズレベル）では、姿勢決定誤差の平均は約 $1$ 〜 $2$ arcsec 程度であり、これは $\sqrt{n}$ の統計的平均化効果と整合します。
右のグラフでは、モンテカルロ結果と理論値 $\sigma / \sqrt{n}$ の比較が示されています。両者がよく一致しており、SVDベースの姿勢推定が理論的な下限に近い性能を実現していることが確認できます。実際のスタートラッカーでは、光学歪曲や焦点距離の校正誤差なども加わるため、この理想的な性能からやや劣化しますが、数アーク秒の精度は十分に達成可能です。
まとめ
本記事では、スタートラッカーの動作原理を星像の取得から姿勢算出まで体系的に解説しました。
- 処理フロー: 撮像 → サブピクセル重心検出 → 星同定（パターンマッチング）→ 姿勢算出の4段階で動作する
- 星同定の鍵: 2つの星の角距離は衛星の姿勢に依存しないという不変性を利用してパターンマッチングを行う
- トライアングルマッチング法: 3つの星が作る三角形の角距離パターンをカタログと照合して星を同定する
- Lost-in-Space問題: 姿勢の事前情報がない場合の星同定は計算量が大きく、効率的なアルゴリズムが必要
- 精度特性: 姿勢精度はクロス軸方向で $\sigma_{\text{star}}/\sqrt{n}$、ロール軸方向ではレバーアーム効果により劣化する
- 運用制約: 太陽回避角、月/地球の影響、SAA通過時の偽星像が実運用上の課題
スタートラッカーは太陽センサと相補的な関係にあります。太陽センサは太陽方向の1本のベクトルを提供し、スタートラッカーは視野内の複数の星の方向ベクトルを提供します。では、これらの複数のベクトル観測から、3軸の姿勢回転行列をどのように最適に推定するのでしょうか？次の記事では、この問題を体系的に解くTRIAD法とQUEST法を解説します。
次のステップとして、以下の記事も参考にしてください。