(57)【要約】【課題】衛星の撮像装置において撮像する方向の誤差を推定する精度を向上する。【解決手段】観測行列算出部３２は、撮像画像と基準画像とのずれを、推定パラメータに観測行列を乗じた積と観測ノイズの和で表す場合の観測行列を計算する。観測特性行列算出部３３は、観測行列の転置行列、観測ノイズの共分散行列の逆行列および観測行列を乗算した観測特性行列を算出する。固有値算出部３４は、観測特性行列の固有値対角行列および固有ベクトル行列を計算する。選択部３５は、固有値対角行列の固有値の大きいものから順に定めた基準で限定固有値を選択し、限定部３６は、限定固有値に対応する限定固有ベクトル行列を生成する。推定部３７は、限定固有値に対応して推定パラメータの数を低減させた主成分空間に変換したパラメータを、最小二乗法を用いて推定し、変換部３８は、パラメータの推定結果を、元空間の推定パラメータに対する推定結果に変換する。【選択図】図５
2015141159.tif 000002
【特許請求の範囲】
【請求項１】
天体物を周回する衛星に搭載される撮像装置の撮像方向を決定する画像ナビゲーション装置であって、
前記撮像装置における画像センサ受光面の画素に対する、前記撮像装置のアラインメント誤差を含む前記撮像装置の視線方向が、前記撮像装置から観測点への前記撮像装置の位置誤差を含む視線方向に一致する関係から規定される、前記撮像装置の誤差の推定パラメータと前記画素の位置との線形化関係式を用いて、前記撮像装置で前記天体物上のランドマークを撮像した画像と基準ランドマーク画像とのずれである観測量を、前記推定パラメータに観測行列を乗じた積と観測ノイズの和で表す場合の前記観測行列を計算する観測行列算出部と、
前記観測行列の転置行列と、前記観測ノイズの共分散行列の逆行列と、前記観測行列とを、この順序で乗算した観測特性行列を算出する観測特性行列算出部と、
前記観測特性行列の固有値対角行列および固有ベクトル行列を計算する固有値算出部と、
前記固有値対角行列に含まれる固有値の大きいものから順に定めた基準で限定固有値を選択する選択部と、
前記固有ベクトル行列から、前記限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列を生成する限定部と、
前記限定固有ベクトル行列を用いて、前記限定固有値に対応して前記推定パラメータの数を低減させた主成分空間に前記線形化関係式を変換した観測方程式に含まれるパラメータを、最小二乗法を用いて推定する推定部と、
前記主成分空間に変換した観測方程式に含まれるパラメータの推定結果を、変換前の元空間における前記推定パラメータに対する推定結果に変換する変換部と、
を備える画像ナビゲーション装置。
【請求項２】
前記変換部で変換した推定結果を、前記天体物表面上の観測地点と前記画像センサ受光面の画素位置との関係を表す観測モデルに代入し、その観測モデルを満たす前記撮像装置の制御量を解くことで、与えられた観測地点を指向するための補正後制御量を計算する、補正部を備える請求項１に記載の画像ナビゲーション装置。
【請求項３】
前記推定部は、最小二乗法で計算する推定誤差共分散に、変動の不確定性の大きさを表す共分散行列を加算して、次回の推定の事前情報として利用する、請求項１または２に記載の画像ナビゲーション装置。
【請求項４】
前記天体物上のランドマークを撮像するごとに、
前記観測行列算出部は、前記観測行列を算出し、
前記観測特性行列算出部は、前記観測特性行列を算出し、
前記固有値算出部は、前記固有値対角行列および前記固有ベクトル行列を算出し、
前記選択部は、閾値以上の限定固有値を選択し、
前記限定部は、前記閾値以上の限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列を生成する、
請求項１から３のいずれか１項に記載の画像ナビゲーション装置。
【請求項５】
前記推定パラメータとして、
前記画像センサ受光面の位置アライメント誤差、及び、回転アライメント誤差、
集光レンズの位置アライメント誤差、及び、回転アライメント誤差、
ミラー回転軸の位置アライメント誤差、及び、回転アライメント誤差、
ミラー回転軸に対するミラー平面の位置アライメント誤差、及び、回転アライメント誤差、
ミラー角エンコーダのゼロ点誤差、
前記衛星に対する前記撮像装置の位置アライメント誤差、及び、回転アライメント誤差、
前記衛星の姿勢決定誤差、もしくは、
前記衛星の軌道決定誤差、
または、それらの組み合わせを含む、請求項１から４のいずれか１項に記載の画像ナビゲーション装置。
【請求項６】
天体物を周回する衛星であって、
撮像装置と、
前記撮像装置の撮像方向を決定する、請求項１から５のいずれか１項に記載の画像ナビゲーション装置と、
を備える衛星。
【請求項７】
天体物を周回する衛星に搭載される撮像装置の撮像方向を決定する画像ナビゲーション装置が行う画像ナビゲーション方法であって、
前記撮像装置における画像センサ受光面の画素に対する、前記撮像装置のアラインメント誤差を含む前記撮像装置の視線方向が、前記撮像装置から観測点への前記撮像装置の位置誤差を含む視線方向に一致する関係から規定される、前記撮像装置の誤差の推定パラメータと前記画素の位置との線形化関係式を用いて、前記撮像装置で前記天体物上のランドマークを撮像した画像と基準ランドマーク画像とのずれである観測量を、前記推定パラメータに観測行列を乗じた積と観測ノイズの和で表す場合の前記観測行列を計算する観測行列算出ステップと、
前記観測行列の転置行列と、前記観測ノイズの共分散行列の逆行列と、前記観測行列とを、この順序で乗算した観測特性行列を算出する観測特性行列算出ステップと、
前記観測特性行列の固有値対角行列および固有ベクトル行列を計算する固有値算出ステップと、
前記固有値対角行列に含まれる固有値の大きいものから順に定めた基準で限定固有値を選択する選択ステップと、
前記固有ベクトル行列から、前記限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列を生成する限定ステップと、
前記限定固有ベクトル行列を用いて、前記限定固有値に対応して前記推定パラメータの数を低減させた主成分空間に前記線形化関係式を変換した観測方程式に含まれるパラメータを、最小二乗法を用いて推定する推定ステップと、
前記主成分空間に変換した観測方程式に含まれるパラメータの推定結果を、変換前の元空間における前記推定パラメータに対する推定結果に変換する変換ステップと、
を備える画像ナビゲーション方法。
【請求項８】
前記変換ステップで変換した推定結果を、前記天体物表面上の観測地点と前記画像センサ受光面の画素位置との関係を表す観測モデルに代入し、その観測モデルを満たす前記撮像装置の制御量を解くことで、与えられた観測地点を指向するための補正後制御量を計算する、補正ステップを備える請求項７に記載の画像ナビゲーション方法。
【請求項９】
前記推定ステップでは、最小二乗法で計算する推定誤差共分散に、変動の不確定性の大きさを表す共分散行列を加算して、次回の推定の事前情報として利用する、請求項７または８に記載の画像ナビゲーション方法。
【請求項１０】
前記天体物上のランドマークを撮像するごとに、
前記観測行列算出部ステップで、前記観測行列を算出し、
前記観測特性行列算出ステップで、前記観測特性行列を算出し、
前記固有値算出ステップで、前記固有値対角行列および前記固有ベクトル行列を算出し、
前記選択ステップでは、閾値以上の限定固有値を選択し、
前記限定ステップでは、前記閾値以上の限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列を生成する、
請求項７から９のいずれか１項に記載の画像ナビゲーション方法。
【請求項１１】
前記推定パラメータとして、
前記画像センサ受光面の位置アライメント誤差、及び、回転アライメント誤差、
集光レンズの位置アライメント誤差、及び、回転アライメント誤差、
ミラー回転軸の位置アライメント誤差、及び、回転アライメント誤差、
ミラー回転軸に対するミラー平面の位置アライメント誤差、及び、回転アライメント誤差、
ミラー角エンコーダのゼロ点誤差、
前記衛星に対する前記撮像装置の位置アライメント誤差、及び、回転アライメント誤差、
前記衛星の姿勢決定誤差、もしくは、
前記衛星の軌道決定誤差、
または、それらの組み合わせを含む、請求項７から１０のいずれか１項に記載の画像ナビゲーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、人工衛星から画像を撮像し、その画像から抽出したランドマークを用いて、目標の指向方向と実際の指向方向との差を求め、最小二乗法などを適用して、撮像装置の誤差を推定する画像ナビゲーション装置、衛星および画像ナビゲーション方法に関する。
【背景技術】
【０００２】
画像ナビゲーション方法では、人工衛星から画像を撮影し、その画像から抽出した地球のエッジまたはランドマークを用いて、画像光学軸と衛星本体軸との差を求め、最小二乗法などを適用して、姿勢誤差およびアライメント誤差を推定する。誤差を推定するパラメータとして考えられているものは、５つの姿勢角およびアライメント誤差などである。詳しくは、例えば、衛星の姿勢をオイラ角で表した、ロール角、ピッチ角およびヨー角の３つの角度、ならびに、衛星本体に対する画像センサのアライメント誤差を表した、ロールアライメント偏差およびピッチアライメント誤差の２つの角度である。
【０００３】
画像ナビゲーション方法では、例えば、ニュートン法を用いた反復計算による非線形最小二乗推定を使用して、これらのパラメータを推定する。さらに、２４時間分の地球エッジ、あるいは、ランドマークの情報を使用して、推定したパラメータの時間変動を、フーリエ級数、および、多項式で近似する。誤差を推定するのは、これらのフーリエ級数、および、多項式の各係数を求めて、各時刻のパラメータを推測することで、次回の地球観測の指向精度を向上させようとするものである（例えば、特許文献１参照）。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特表２００２−５１３４６４号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
実際の衛星では、各種センサのアライメント誤差、各種センサのノイズ成分、および、バイアス成分、レンズやミラーなどの各種光学機器のアライメント誤差など、さまざまなずれ、あるいは、誤差が、位置単位では数〜数百μｍのオーダーで、角度単位では、数〜数百μradのオーダーで発生している。例えば、地球観測を行う静止衛星の場合では、１００μrad指向方向がずれると、地球表面上では約３.５ｋｍ観測地点がずれてしまう。地球表面上で１００ｍオーダーの位置決定精度で地球観測を行う場合、これらのずれまたは誤差は、観測精度に大きな影響をもたらす。
【０００６】
実際の衛星では、前述のように、さまざまなずれ／誤差が発生しているが、従来の画像ナビゲーション方法では、推定するパラメータとして、高々５つ程度のずれ／誤差のみを選んで、これらのずれ／誤差のみを推定する。観測条件が変化しても常に同じパラメータを推定し続けるため、観測条件によって高い精度で観測対象を指向することが難しい、といった課題がある。
【０００７】
本発明は、上述の事情に鑑みてなされたもので、衛星の撮像装置において撮像する方向の誤差を推定する精度を向上することを目的とする。
【課題を解決するための手段】
【０００８】
上述の目的を達成するため、この発明に係る画像ナビゲーション装置は、天体物を周回する衛星に搭載される撮像装置の撮像方向を決定する画像ナビゲーション装置であって、観測行列算出部、観測特性行列算出部、固有値算出部、選択部、限定部、推定部および変換部を備える。観測行列算出部は、撮像装置の誤差の推定パラメータと画素の位置との線形化関係式を用いて、撮像装置で天体物上のランドマークを撮像した画像と基準ランドマーク画像とのずれである観測量を、推定パラメータに観測行列を乗じた積と観測ノイズの和で表す場合の観測行列を計算する。推定パラメータと画素の位置との線形化関係式は、撮像装置における画像センサ受光面の画素に対する、撮像装置のアラインメント誤差を含む撮像装置の視線方向が、撮像装置から観測点への撮像装置の位置誤差を含む視線方向に一致する関係から規定される。観測特性行列算出部は、観測行列の転置行列と、観測ノイズの共分散行列の逆行列と、観測行列とを、この順序で乗算した観測特性行列を算出する。固有値算出部は、観測特性行列の固有値対角行列および固有ベクトル行列を計算する。選択部は、固有値対角行列に含まれる固有値の大きいものから順に定めた基準で限定固有値を選択する。限定部は、固有ベクトル行列から、限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列を生成する。推定部は、限定固有ベクトル行列を用いて、限定固有値に対応して推定パラメータの数を低減させた主成分空間に線形化関係式を変換した観測方程式に含まれるパラメータを、最小二乗法を用いて推定する。変換部は、主成分空間に変換した観測方程式に含まれるパラメータの推定結果を、変換前の元空間における推定パラメータに対する推定結果に変換する。
【発明の効果】
【０００９】
本発明によれば、観測特性行列の固有値対角行列に含まれる固有値の大きいものから順に定めた基準で限定固有値を選択し、限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列を用いて、限定固有値に対応して推定パラメータの数を低減させ、可観測性の大きいパラメータに限定してパラメータを推定するので、衛星の撮像装置において撮像する方向の誤差を推定する精度を向上することができる。
【図面の簡単な説明】
【００１０】
【図１】実施の形態に係る画像ナビゲーションの概念を示す図である。
【図２】画像ナビゲーションのパラメータ推定に関わる部分を示す図である。
【図３】撮像装置の構成配置及び座標系の定義を表す図である。
【図４】軌道固定座標系と地球との幾何学的な関係を表す図である。
【図５】本発明の実施の形態１に係る画像ナビゲーション装置の構成例を示すブロック図である。
【図６】実施の形態１に係る推定パラメータの推定の一例を示すフローチャートである。
【図７】本発明の実施の形態２に係る画像ナビゲーション装置の構成例を示すブロック図である。
【図８】実施の形態２に係る補正の動作の一例を示すフローチャートである。
【図９】本発明の実施の形態３に係る推定パラメータの推定の一例を示すフローチャートである。
【図１０】本発明の実施の形態４に係る補正の動作の一例を示すフローチャートである。
【発明を実施するための形態】
【００１１】
地球を周回する衛星が、カメラなどの撮像装置を用いて、地球を観測する場合を対象に考える。衛星には、撮像装置のアライメント誤差、姿勢決定誤差および軌道決定誤差があり、これらの影響で、撮像装置の視線方向を観測すべき地点に指向させても、場合によって地球表面上で数百ｍから数ｋｍ以上ずれた地点を指向することがある。画像ナビゲーションは、観測した画像を用いて、衛星が有するアライメント誤差、姿勢決定誤差および軌道決定誤差を推定し、推定した結果を用いて、次回の観測すべき地点をより精度よく指向するための手法である。以下、地球を周回する人工衛星に搭載する撮像装置を想定して説明するが、実施の形態は地球以外の天体物、例えば、他の惑星または衛星（月）を周回する人工衛星にも適用できる。
【００１２】
図１は、画像ナビゲーションの概念を示す図である。地球を周回する衛星４１に撮像装置４２を搭載し、地球観測撮像領域５３およびランドマーク撮像領域５４を撮像する。撮像装置４２は、１次鏡４４、２次鏡４３などのミラーを回転させることで、観測すべき地点およびランドマークを指向する。ランドマークは、地球上の島、半島、その他特徴のある地形の中からあらかじめ選定された数百の地点である。衛星には、ランドマーク周囲の地形情報を抽出した基準ランドマーク画像４８がデータベースに記憶されている。ランドマークは、画像ナビゲーションの推定処理において、実際に指向している視線方向４６を決定するために用いられる。画像ナビゲーションでは、複数のランドマークを観測し、観測したランドマーク画像４７と、データベースに保存された対応する基準ランドマーク画像４８とを比較し、双方の画像の間でどれぐらいの移動ずれ量があるかを、相関部４９にて計算する。撮像装置４２の視線方向をランドマークに指向させて、実際に観測してみると、観測したランドマーク画像４７に、基準ランドマーク画像４８からの移動ずれ量が検出される。
【００１３】
この移動ずれ量を推定部５０に渡し、アライメント誤差、姿勢決定誤差、軌道決定誤差を推定パラメータとして推定処理を行う。推定されたアライメント誤差、姿勢決定誤差および軌道決定誤差を、補正部５１に渡し、これらにより更新された光学系の数学モデルを用いて、観測すべき地点を指向するためのミラー回転角を計算する。このミラー回転角を用いてミラーを回転させ、補正観測部５２にてより精度よく観測すべき地点を指向した観測を行う。
【００１４】
図２は、画像ナビゲーションのパラメータ推定に関わる部分を示すブロック図である。画像ナビゲーション開始６１から、通常の地球観測を行うモード６３と、画像ナビゲーションにより推定処理を行う画像ナビゲーションモード６２の２種類のモードがある。地球観測モード６３を実行する場合は、補正部５１にて、その時に保持しているアライメント誤差、姿勢決定誤差および軌道決定誤差の値を用いた光学系の数学モデルに基づき、観測すべき地点を指向するためのミラー回転角を計算する。そして、補正観測部５２にて地球観測を行う。画像ナビゲーションによる推定処理を行う前の状態では、アライメント誤差、姿勢決定誤差および軌道決定誤差の値は、すべてゼロに設定され、それ以降で画像ナビゲーションによる推定処理を行えば、その都度、これらの値が更新される。
【００１５】
画像ナビゲーションモード６２を実行する場合は、ランドマーク（ＬＭ）の観測部６４にて、複数のランドマークを観測し、複数のランドマーク画像６５を取得する。基準ランドマーク画像の参照部６６にて対応する基準ランドマーク画像６７をデータベースから取得し、相関部４９にて、双方の画像の間でどれぐらいの移動ずれ量６９があるかを計算する。目的のランドマークを指向したはずが、実際に観測してみると、上記で計算される移動ずれ量６９があったことが検出される。この移動ずれ量６９を推定部５０に渡し、アライメント誤差、姿勢決定誤差および軌道決定誤差を推定パラメータ７１として推定処理を行う。
【００１６】
推定されたアライメント誤差、姿勢決定誤差および軌道決定誤差を、補正部５１に渡し、これらにより更新された光学系の数学モデルを用いて、観測すべき地点を指向するためのミラー回転角、および、その補正量７３を計算する。これを用いてミラーを回転させ、補正観測部５２にてより精度よく観測すべき地点を指向した観測を行う。
【００１７】
（１）座標系の定義
図３は、撮像装置の構成配置及び座標系の定義を表す図である。ここで、撮像装置の光学系における座標系の定義を説明し、撮像装置の誤差を推定するパラメータ（推定パラメータ）の取り方の例を示す。撮像装置の構成配置として、例えば、図３に示すものを考える。合わせて、衛星固定、あるいは、各光学機器固定の座標系を定義する。ただし、本実施の形態は、衛星に搭載する撮像装置がどのような構成配置であっても適用可能である。
【００１８】
撮像装置は、例えば、画像センサ２、集光レンズ３、２次鏡４、１次鏡６および開口部８から構成される。撮像装置の光学系における光軸１は、画像センサ２の受光面から集光レンズ３を通って、２次鏡４および１次鏡６で折り返されて、開口部８に向かう。
【００１９】
観測すべき平行光は、開口部８より入射し、１次鏡６により反射して２次鏡４に向かい、さらに２次鏡４により再度反射して、集光レンズ３により画像センサ２の受光面上に集光される。１次鏡６および２次鏡４は、１次鏡回転軸７および２次鏡回転軸５の周りにそれぞれ回転する機構を有する。
【００２０】
画像センサ２の受光面に固定された座標系を画像センサ受光面固定座標９の方向に取り、座標系｛Ｉ｝と記す。２次鏡４のミラー平面に固定された座標系を２次鏡ミラー平面固定座標１０の方向に取り、座標系｛Ｍ２｝と記す。１次鏡ミラー平面に固定された座標系を１次鏡ミラー平面固定座標１１の方向に取り、座標系｛Ｍ１｝と記す。衛星に固定された座標系を衛星固定座標１２の方向に取り、座標系｛Ｕ｝と記す。衛星の周回軌道に固定された座標系を、衛星の姿勢回転がゼロ位置である際の衛星固定座標１２と同じ方向に取り、座標系｛Ｙ｝と記す。
【００２１】
（２）推定パラメータの定義
撮像装置４２の光学系の構成配置として、例えば、図３に示すものを考えた場合、画像ナビゲーションの性能に影響を与える可能性のある、アライメント誤差、姿勢決定誤差および軌道決定誤差を列挙すると、以下に挙げる誤差項目が考えられる。それぞれの光学機器がずれを持つ可能性があるので、すべての光学機器に対して、位置アライメント誤差、および、回転アライメント誤差を考えればよい。
【００２２】
画像センサ２の受光面の位置アライメント誤差、および、回転アライメント誤差
集光レンズ３の位置アライメント誤差、および、回転アライメント誤差
２次鏡回転軸５の位置アライメント誤差、および、回転アライメント誤差
２次鏡回転軸５に対する２次鏡平面の位置アライメント誤差、および、回転アライメント誤差
１次鏡回転軸７の位置アライメント誤差、および、回転アライメント誤差
１次鏡回転軸７に対する１次鏡平面の位置アライメント誤差、および、回転アライメント誤差
１次鏡６および２次鏡４のミラー角エンコーダのゼロ点誤差
衛星４１に撮像装置４２を据え付ける際の位置アライメント誤差、および、回転アライメント誤差
衛星４１の姿勢決定誤差
衛星４１の軌道決定誤差
【００２３】
位置アライメント誤差は、機器を設置した時にずれる誤差のうち、直交する３軸方向（ｘ、ｙ、ｚ）に平行移動してずれる位置ずれを表し、３つの変数で定義できる。回転アライメント誤差は、このうち３軸方向（ｘ、ｙ、ｚ）それぞれを中心に回転してずれる回転ずれを表し、３つの変数で定義できる。
【００２４】
図３に示す構成配置を考えた場合、上記に挙げた誤差項目のうち、考慮しなくてよいもの、または、集約可能なものがいくつか存在する。集光レンズ３の中心を光学機器全体の測定基準とすると、集光レンズの位置アライメント誤差は、衛星の軌道決定誤差と同じになり、回転アライメント誤差は、衛星の姿勢決定誤差と同じになり、それぞれ集約可能である。ミラーの位置が変わっても、視線方向は変わらないので、ミラーの位置アライメント誤差は、画像ナビゲーションの性能に影響を与えない。従って、ミラー面垂直方向の回転アライメント誤差は、推定パラメータとして考慮しなくてよい。ミラー面垂直方向に対する回転は、鏡面が平行に回転するのみで、視線方向を変化させるものではなく、推定パラメータとして考慮しなくてよい。ミラー回転軸に対するミラー平面の回転アライメント誤差のうち、ミラー回転軸と同じ方向の回転誤差は、ミラー回転軸の回転アライメント誤差と同じになり、集約可能である。衛星に搭載光学機器を据え付ける際の位置アライメント誤差、および、回転アライメント誤差は、衛星の姿勢決定誤差と全く同じ変数として表現されるため、姿勢決定誤差に集約可能である。
【００２５】
また、上記の誤差項目のうち、画像ナビゲーションへの影響が十分小さいため、無視可能なものがいくつか存在する。通常、位置アライメント誤差は、１００μｍ以下、回転アライメント誤差は、１００μrad以下になる。この範囲における影響を考えると、通常の搭載光学機器構成において受光面が同方向に１００μｍ移動しても、受光面で観測される画像は最も大きくて０.０５pixel以下しか動かないことが、幾何学的な関係から導出できる。画像センサ受光面の位置アライメント誤差のうち、受光面垂直方向の位置誤差は影響が十分小さいため、この位置アライメント誤差は無視できる。画像センサ受光面が受光面垂直方向に１００μrad回転しても、受光面で観測される画像は、最も大きくて０.０５pixel以下しか動かないことが、幾何学的な関係から導出できる。画像センサ受光面の回転アライメント誤差のうち、受光面垂直方向の回転誤差は、影響が十分小さいため、この回転アライメント誤差は無視できる。画像センサ受光面が受光面水平方向に１００μrad回転しても、受光面で観測される画像は、最も大きくて０.０００５pixel以下しか動かないことが、幾何学的な関係から導出できる。画像センサ受光面の回転アライメント誤差のうち、受光面水平方向の回転誤差は、影響が十分小さいため、この回転アライメント誤差は無視できる。
【００２６】
以上をまとめると、画像ナビゲーションの性能に影響を与える、アライメント誤差、姿勢決定誤差、軌道決定誤差は、以下に挙げる１６個の項目を考えればよい。性能に影響を与えるこれらの各種誤差を、推定パラメータとして推定することで、画像ナビゲーションの性能を向上させることができる。
【００２７】
画像センサ受光面の位置アライメント誤差のうち、受光面水平方向の位置２成分：Δsx、Δsyと置く。
２次鏡回転軸の回転アライメント誤差のうち、回転軸の傾きを表す回転２成分：ΔM2rx、ΔM2rzと置く。
２次鏡回転軸に対する２次鏡平面の回転アライメント誤差のうち、ミラー平面の傾きを表す回転１成分：ΔM2txと置く。
１次鏡回転軸の回転アライメント誤差のうち、回転軸の傾きを表す回転２成分：ΔM1rx、ΔM1rzと置く。
１次鏡回転軸に対する１次鏡平面の回転アライメント誤差のうち、ミラー平面の傾きを表す回転１成分：ΔM1txと置く。
１次鏡、および、２次鏡のミラー角エンコーダのゼロ点誤差の２成分：ΔNS、ΔEWと置く。
衛星の姿勢決定誤差の回転３成分：ΔAx、ΔAy、ΔAzと置く。
衛星の軌道決定誤差の位置３成分：ΔLx、ΔLy、ΔLzと置く。
【００２８】
なお、搭載光学機器の構成配置が図３に示すものと異なる場合は、適切な推定パラメータを定義し直すことで、以降、同様の手順で推定を実施することができる。
【００２９】
（３）観測モデルおよびその線形化
以降では、前述した複数の推定パラメータをどのように数式モデルとして定式化し、本モデルを用いて、本発明の推定手法をどのように実施するかについて説明する。
【００３０】
アライメント誤差、姿勢決定誤差および軌道決定誤差を推定パラメータとして推定を行うためには、観測量と推定パラメータとの間の関係を数式モデルとして定式化する必要がある。地球の観測画像を取得するので、この場合の観測量は、地球表面上の任意の地点を観測した際、画像に映るその地点の画素位置である。アライメント誤差、姿勢決定誤差および軌道決定誤差のいずれも全くない理想状態では、任意の地点を観測した場合に、画像に映るその地点の画素位置は常にゼロ、すなわち、中心となるように構成される。したがって、観測量は、ランドマークを撮像した画像と基準ランドマーク画像とのずれである。
【００３１】
数式モデルを構築するには、画像センサ受光面の任意の画素から集光レンズ中心方向に向かって出力した視線方向ベクトルを出発点として、各光学機器の中に視線方向ベクトルを入射させ、アライメント誤差による視線方向のずれを含めながら光学機器から出力される視線方向ベクトルを代数計算する。図３の構成配置に従って視線方向を順番に光学機器を通過させて、最終的に搭載光学機器の開口部における視線方向ベクトルを代数計算する。これに、姿勢決定誤差および軌道決定誤差の影響を加えて、開口部から地球表面上に到達する視線方向ベクトルの位置を代数計算する。これによって、地球表面上の任意の地点と画像センサ受光面の画素との関係が、推定パラメータを含んだ形式で得られる。地球表面上の任意の地点を、観測すべき地点に固定すれば、最終的に、観測量である画素位置と推定パラメータとの間の関係式が得られる。
【００３２】
（３．１）集光レンズを通過した後の視線方向
例えば、図３の構成配置における、観測量と推定パラメータとの間の関係式を導出する。はじめに、画像センサ受光面の任意の画素位置ｐ＝［ｐx、ｐy］T（右肩のTは転置を表す）から集光レンズ３の中心方向に向かう、衛星固定座標系｛Ｕ｝で表される視線方向ベクトルＵｚinは、
Ｕｚin＝normalize（Ｕ［−ｆｃ  −ｐx−Δsx  ｐy＋Δsy］T）    （１）
で表される。ここでは紙面の便宜上、列ベクトルを行ベクトルの転置Tで記載している。ベクトルの左肩の文字は、ベクトルを表現する座標系を示す。ｆｃは、集光レンズ３の焦点距離である。視線方向ベクトルを単位ベクトルとするため、normalize（・）関数を付加してある。normalize（・）は、ベクトルｘをその大きさ｜ｘ｜で除して単位長さのベクトルｘ／｜ｘ｜とする操作を表す。
【００３３】
（３．２）２次鏡を通過した後の視線方向
次に、２次鏡４を通過した後の視線方向ベクトルUｚ2を代数計算する。そのために、衛星固定座標系｛Ｕ｝から２次鏡ミラー平面固定座標系｛Ｍ２｝へ方向ベクトルを変換する座標変換行列Ｍ２ＣＵを代数計算する。変換行列表記の右添え字は変換前の座標系、左肩の文字は変換後の座標系を示す。まず、座標系｛Ｕ｝から、２次鏡４の回転軸に沿った中間的な座標系｛Ａ２｝に変換し、その座標系上で回転アライメント誤差ΔM2rx、ΔM2rzによるずれを含めた座標系｛Ｂ２｝に変換する。そして、座標系｛Ｂ２｝上で２次鏡回転軸周りにミラーを回転させた座標系｛Ｃ２｝に変換して、座標系｛Ｃ２｝上で回転アライメント誤差ΔM2txによるずれを含めて最終的なミラー平面固定座標系｛Ｍ２｝に変換する。すなわち、
Ｍ２ＣＵ＝Ｍ２ＣＣ２・Ｃ２ＣＢ２・Ｂ２ＣＡ２・Ａ２ＣＵ    （２）
である。それぞれの座標変換行列は、以下のとおり。
Ａ２ＣＵ＝Ｃx（−π／２）Ｃz（−π／４）    （３）
Ｂ２ＣＡ２＝Ｃz（ΔM2rz）Ｃx（ΔM2rx）    （４）
Ｃ２ＣＢ２＝Ｃy（θEW＋ΔEW）    （５）
Ｍ２ＣＣ２＝Ｃx（ΔM2tx）    （６）
【００３４】
ここでθEWは、２次鏡４のミラー回転角である。Ｃx（・）、Ｃy（・）、Ｃz（・）は、それぞれｘ軸、ｙ軸、ｚ軸周りに座標を回転させる座標変換行列である。
【数１】
2015141159.tif 000003
【数２】
2015141159.tif 000004
【数３】
2015141159.tif 000005
【００３５】
代数計算を行うと、式（２）の座標変換行列は以下のように表される。
【数４】
2015141159.tif 000006
【数５】
2015141159.tif 000007
【数６】
2015141159.tif 000008
【数７】
2015141159.tif 000009
【００３６】
ミラー平面の法線方向がｚ軸となるミラーを用いた場合の、方向ベクトルを反射させる行列Ｍは、以下の式（１４）で表せる。
【数８】
2015141159.tif 000010
２次鏡４を通過した後の視線方向ベクトルUｚ2は、入射する視線方向ベクトルＵｚinを、衛星固定座標系｛Ｕ｝から２次鏡ミラー平面固定座標系｛Ｍ２｝へ変換し、座標系｛Ｍ２｝上でミラーによる反射行列Ｍを作用させ、再び２次鏡ミラー平面固定座標系｛Ｍ２｝から衛星固定座標系｛Ｕ｝に変換することで得られる。すなわち、
Ｕｚ2＝ＵＣＭ２・Ｍ・Ｍ２ＣＵ・Ｕｚin    （１５）
である。ＵＣＭ２は、Ｍ２ＣＵの逆変換なのでＭ２ＣＵの逆行列になるが、ここでの座標変換行列はすべて直交行列なので、逆行列は転置行列と同じになり、ＵＣＭ２＝Ｍ２ＣＵTにより計算可能である。
【００３７】
（３．３）１次鏡を通過した後の視線方向
次に、１次鏡６を通過した後の視線方向ベクトルＵｚ1を代数計算する。２次鏡４の時と同様に、まず、衛星固定座標系｛Ｕ｝から１次鏡ミラー平面固定座標系｛Ｍ１｝へ方向ベクトルを変換する座標変換行列Ｍ１ＣＵを代数計算する。
Ｍ１ＣＵ＝Ｍ１ＣＣ１・Ｃ１ＣＢ１・Ｂ１ＣＡ１・Ａ１ＣＵ    （１６）
それぞれの座標変換行列は、以下のとおり。
Ａ１ＣＵ＝Ｃz（−π／２）Ｃx（π／４）    （１７）
Ｂ１ＣＡ１＝Ｃz（ΔM1rz）Ｃx（ΔM1rx）    （１８）
Ｃ１ＣＢ１＝Ｃy（θNS＋ΔNS）    （１９）
Ｍ１ＣＣ１＝Ｃx（ΔM1tx）    （２０）
【００３８】
θNSは、１次鏡のミラー回転角である。代数計算を行うと、式（１６）の座標変換行列は以下のように表される。
【数９】
2015141159.tif 000011
【数１０】
2015141159.tif 000012
【数１１】
2015141159.tif 000013
【数１２】
2015141159.tif 000014
【００３９】
１次鏡６を通過した後の視線方向ベクトルＵｚ1は、入射する視線方向ベクトルUｚ2に対して、衛星固定座標系｛Ｕ｝から１次鏡ミラー平面固定座標系｛Ｍ１｝へ変換し、座標系｛Ｍ１｝上でミラーによる反射行列Ｍを作用させ、再び１次鏡ミラー平面固定座標系｛Ｍ１｝から衛星固定座標系｛Ｕ｝に変換することで得られる。すなわち、
Ｕｚ1＝ＵＣＭ１・Ｍ・Ｍ１ＣＵ・Ｕｚ2    （２５）
である。ＵＣＭ１は、Ｍ１ＣＵの逆変換なので、Ｍ１ＣＵの逆行列になるが、ここでの座標変換行列はすべて直交行列なので、逆行列は転置行列と同じになり、ＵＣＭ１＝Ｍ１ＣＵTにより計算可能である。
【００４０】
（３．４）衛星から出力される視線方向
次に、衛星から出力される視線方向ベクトルＹｚoutを代数計算する。入射する視線方向ベクトルＵｚ1に対して、姿勢決定誤差を表す座標変換行列を作用させて、衛星固定座標系｛Ｕ｝から軌道固定座標系｛Ｙ｝へ変換することで代数計算する。すなわち、
Ｙｚout＝Ｃz（ΔAz）Ｃy（ΔAy）Ｃx（ΔAx）ＹＣＵ・Ｕｚ1    （２６）
である。通常、地球を観測する際は、衛星固定座標系｛Ｕ｝が軌道固定座標系｛Ｙ｝に一致する姿勢で観測を行うので、ここでは、ＹＣＵ＝Ｉ（：単位行列）とする。
【００４１】
これまでの計算結果をまとめて、最終的に得られる視線方向ベクトルＹｚoutは、式（１）、式（１５）、式（２５）および式（２６）から、次のようになる。
Ｙｚout＝Ｃz（ΔAz）Ｃy（ΔAy）Ｃx（ΔAx）ＹＣＵ
・ＵＣＭ１・Ｍ・Ｍ１ＣＵ・ＵＣＭ・２Ｍ・Ｍ２ＣＵ
・normalize（Ｕ［−ｆｃ  −ｐx−Δsx  ｐy＋Δsy］T）    （２７）
式（２７）で得られた、衛星から出力される視線方向ベクトルＹｚoutを表す関係式を、次のように置く。
Ｙｚout＝ｈ（θEW，θNS，Δｘ，ｐ）    （２８）
ただし、変数Δｘ、ｐは、次のように置く。
【数１３】
2015141159.tif 000015
【００４２】
（３．５）視線方向が地球表面上に到達する位置
次に、衛星から出力される視線方向ベクトルＹｚoutと、Ｙｚoutが地球表面上に到達する位置との間の関係式を代数計算する。軌道固定座標系｛Ｙ｝の上で視線方向ベクトルを考える。撮像装置の開口部８の中心を軌道固定座標系｛Ｙ｝の原点に採り、地球中心方向をｚ軸、衛星の進行方向をｘ軸とする。
【００４３】
図４は、軌道固定座標系と地球との幾何学的な関係を表す図である。便宜上、仮想的な地球表面２７を想定する。軌道固定座標系｛Ｙ｝のｚ軸は、地球２１の中心２２に向かっている。ｚ軸と直交して地球表面に接する平面を仮想的な地球表面２７と定義し、地球２１の表面を平面として扱えるようにする。仮想的な地球表面２７の上に、観測する任意の地点２４をＬＭ＝（ＬＭx，ＬＭy）と取る。仮想的な地球表面２７は、軌道固定座標系｛Ｙ｝の原点２３から地球２１のエッジを見込む円錐２６を境界とする。仮想的な地球表面２７上の点と、本来の曲面の地球表面上の点との間は、別途、幾何学的な性質を用いて相互に変換することができるため、仮想的な地球表面２７を用いても何ら最終的な機能を失うものではない。
【００４４】
このとき、観測する任意の地点２４の座標ＬＭを、軌道固定座標系｛Ｙ｝上で、
ＬＭ＝Ｙ［ＬＭx  ＬＭy  Ｒse］T    （３０）
と表せる。式（３０）では、軌道固定座標系｛Ｙ｝の原点２３と仮想的な地球表面２７との距離２８をＲseとしている。衛星から出力された視線方向２５のベクトルＹｚoutが観測地点ＬＭの方向を向くとき、Ｙｚoutは次のように表される。
Ｙｚout＝normalize（Ｙ［ＬＭx−ΔLx  ＬＭｙ−ΔLy  Ｒse−ΔLz］T）  （３１）
式（３１）は、衛星の位置Ｙ［０，０，０］Tに、軌道決定誤差Ｙ［ΔLx，ΔLy，ΔLz］Tを加えて、観測地点ＬＭとの差を取ることで得られる。式（３１）で得られた、視線方向ベクトルＹｚoutと観測地点ＬＭとの間の関係を次以下の式（３２）で定義する。
Ｙｚout＝ｌ（ＬＭ，ΔLx，ΔLy，ΔLz）  （３２）
【００４５】
（３．６）画素位置と推定パラメータとの間の関係式
式（２８）のｈ（θEW、θNS、Δｘ、ｐ）と式（３２）のｌ（ＬＭ、ΔLx、ΔLy、ΔLz）は、衛星から出力される同じ視線方向ベクトルＹｚoutを表す。式（２８）では、アライメント誤差、および、姿勢決定誤差とＹｚoutとの関係を、式（３２）では、軌道決定誤差とＹｚoutとの関係を表す。地球表面上の任意の地点を観測した場合、式（２８）と式（３２）のＹｚoutは同じ値になるので、次の関係式が成り立つ。
Ｙｚout＝ｈ（θEW，θNS，Δｘ，ｐ）＝ｌ（ＬＭ，ΔLx，ΔLy，ΔLz）  （３３）
【００４６】
式（３３）により、地球表面上の任意の地点ＬＭと画像センサ受光面の画素位置ｐとの関係が、推定パラメータΔｘ、ΔLx、ΔLy、ΔLzを含んだ形式で得られる。すなわち、地球表面上の任意の地点を観測すべき地点に固定すれば、最終的に、観測量である画素位置ｐと推定パラメータΔｘ、および、ΔLx、ΔLy、ΔLzとの間の関係式が得られる。式（３３）は、非線形の場合の完全な観測モデル式である。次に、この関係式を線形化する。後述の推定では、この関係式を線形化した式を使用する。
【００４７】
（３．７）衛星から出力される視線方向の関係式を線形化
式（２８）を線形化する。通常、位置アライメント誤差は１００μｍ以下、回転アライメント誤差は１００μrad以下、姿勢決定誤差は１００μrad以下といった非常に小さな値となる。そこで、アライメント誤差についてテーラー展開を施し、線形化を行って数式モデルを簡単化する。これにより、２次以上の微小項は無視されるが、２次の微小項は、（１００μ）２ｍ＝（１０−４）２ｍ＝１０−８ｍ、すなわち１０ｎｍ以下、あるいは、（１００μ）２rad＝（１０−４）２rad＝１０−８rad、すなわち１０ｎrad以下のオーダーになると見積もられる。画像ナビゲーションでは、１０〜１００μｍオーダー、あるいは、１０〜１００μradオーダーの誤差を対象にするので、これらに比べると、１０−４倍小さな１０ｎｍ以下、あるいは、１０ｎrad以下の近似誤差しか生じ得ず、線形近似を行っても最終的な画像ナビゲーションの性能にはほとんど影響がないと言える。
【００４８】
上記範囲で最大のアライメント誤差、および、通常生じ得る最大の姿勢決定誤差と軌道決定誤差があった場合、観測量ｐx、ｐyは、１０μｍ〜１ｍｍオーダーの値になる。Ｕｚinは正規化されているから、式の中でｐx、ｐyは常にｐx／ｆｃ、ｐy／ｆｃの形で現れる。ｆｃは１０ｃｍ〜１ｍであるから、ｐx、ｐyとｆｃの大きさを比べると、通常、１０２〜３倍のオーダーでｆｃの方が大きい。ｐx／ｆｃ、ｐy／ｆｃは、０.０１以下の小さな値になるため、ｐx、ｐyを線形化できる。０．０１程度では、線形近似による誤差が大きくなる場合も考えられるが、通常、画像ナビゲーションでは１回の推定だけでは終了せず、衛星の運用中は、常時、観測と推定を繰り返す。これにより、アライメント誤差、姿勢決定誤差および軌道決定誤差の推定値が、次第に真値に近づいていくため、観測量ｐx、ｐyは１０μｍ程度に、ｐx／ｆｃ、ｐy／ｆｃは、１０−４程度に小さくなっていく。そのため、線形近似を行っても、何度か観測と推定を繰り返しある程度時間が経過した後の最終的な画像ナビゲーションの性能にはほとんど影響がないと言える。
【００４９】
アライメント誤差および姿勢決定誤差、ならびに、観測量は十分小さいので、式（２８）で得られた視線方向ベクトルＹｚoutにおいて、アライメント誤差Δsx、Δsy、ΔM2rz、ΔM2rx、ΔEW、ΔM2tx、ΔM1rz、ΔM1rx、ΔNS、ΔM1tx、姿勢決定誤差ΔAx、ΔAy、ΔAz、および、観測量ｐx、ｐyについて線形化し、２次以上の微小項を無視する。これにより、衛星から出力される視線方向ベクトルＹｚoutは、最終的に以下の式（３４）で表せる。
Ｙｚout＝Ｃh,0＋Ｃh,px・ｐx＋Ｃh,py・ｐy
＋Ｃh,Sx・ΔSx＋Ｃh,Sy・ΔSy
＋Ｃh,M2rx・ΔM2rx＋Ｃh,M2rz・ΔM2rz＋Ｃh,EW・ΔEW＋Ｃh,M2tx・ΔM2tx
＋Ｃh,M1rx・ΔM1rx＋Ｃh,M1rz・ΔM1rz＋Ｃh,NS・ΔNS＋Ｃh,M1tx・ΔM1tx
＋Ｃh,Ax・ΔAx＋Ｃh,Ay・ΔAy＋Ｃh,Az・ΔAz        （３４）
【００５０】
このとき、各係数は、以下のような構成となる。
【数１４】
2015141159.tif 000016
【数１５】
2015141159.tif 000017
【数１６】
2015141159.tif 000018
【数１７】
2015141159.tif 000019
【数１８】
2015141159.tif 000020
【数１９】
2015141159.tif 000021
【数２０】
2015141159.tif 000022
【数２１】
2015141159.tif 000023
【数２２】
2015141159.tif 000024
【数２３】
2015141159.tif 000025
【数２４】
2015141159.tif 000026
【数２５】
2015141159.tif 000027
【数２６】
2015141159.tif 000028
【数２７】
2015141159.tif 000029
【数２８】
2015141159.tif 000030
【数２９】
2015141159.tif 000031
【００５１】
（３．８）視線方向が地球表面上に到達する位置の関係式を線形化
次に、式（３１）を線形化する。式（３１）において、normalize（・）内部のベクトル全成分をＲseで割って、以下の式（５１）に変形する。
【数３０】
2015141159.tif 000032
【００５２】
通常、軌道決定誤差は３ｋｍ以下の値である。一方、Ｒseは３５７８６kmであり、ΔLx／Ｒse、ΔLy／Ｒse、ΔLz／Ｒseは１０−４以下の小さな値になる。そのため、Ｒseで割った軌道決定誤差ΔLx／Ｒse、ΔLy／Ｒse、ΔLz／Ｒseについてテーラー展開を施し、線形化を行って数式モデルを簡単化する。これにより、２次以上の微小項は無視されるが、２次の微小項は、１０−８オーダーの近似誤差しか生じ得ず、線形近似を行っても最終的な画像ナビゲーションの性能にはほとんど影響がないと言える。ここで、Ｒseで正規化した軌道決定誤差を次のように置く。
Δ'Lx＝ΔLｘ／Ｒse、Δ'Ly＝ΔLy／Ｒse、Δ'Lz＝ΔLz／Ｒse    （５２）
ΔL＝［Δ'Lx，Δ'Ly，Δ'Lz］T                              （５３）
これにより、衛星から出力される視線方向ベクトルＹｚoutは、最終的に次のような式になる。
Ｙｚout＝Ｃl,0＋Ｃl,Lx・Δ'Lx＋Ｃl,Ly・Δ'Ly＋Ｃl,Lz・Δ'Lz  （５４）
【００５３】
このとき、各係数は、以下のような構成となる。
【数３１】
2015141159.tif 000033
【数３２】
2015141159.tif 000034
【数３３】
2015141159.tif 000035
【数３４】
2015141159.tif 000036
【００５４】
ここで、式（３２）で定義した関数ｌ（ＬＭ、ΔLx、ΔLy、ΔLz）を、ΔLを用いて、以下の式（５９）で定義することもできる。
【数３５】
2015141159.tif 000037
【００５５】
（３．９）画素位置と推定パラメータとの間の線形化関係式
式（３４）と式（５４）は、衛星から出力される同じ視線方向ベクトルＹｚoutを表す。式（３４）では、アライメント誤差および姿勢決定誤差と、Ｙｚoutとの関係を、式（５４）では、軌道決定誤差とＹｚoutとの関係を表す。地球表面上の任意の地点を観測した場合、式（３４）と式（５４）のＹｚoutは同じ値になるので、
Ｙｚout＝Ｃh,0＋Ｃh,px・ｐx＋Ｃh,py・ｐy
＋Ｃh,Sx・ΔSx＋Ｃh,Sy・ΔSy
＋Ｃh,M2rx・ΔM2rx＋Ｃh,M2rz・ΔM2rz＋Ｃh,EW・ΔEW＋Ｃh,M2tx・ΔM2tx
＋Ｃh,M1rx・ΔM1rx＋Ｃh,M1rz・ΔM1rz＋Ｃh,NS・ΔNS＋Ｃh,M1tx・ΔM1tx
＋Ｃh,Ax・ΔAx＋Ｃh,Ay・ΔAy＋Ｃh,Az・ΔAz
＝Ｃl,0＋Ｃl,Lx・Δ'Lx＋Ｃl,Ly・Δ'Ly＋Ｃl,Lz・Δ'Lz            （６０）
ここで、Ｃh,0とＣl,0は同じ値となり、式から消去される。Ｃl,0のＬＭx、ＬＭyにおいて、幾何学的な性質からθEW、θNSを用いた表現に直すと、Ｃl,0の表現はＣh,0と完全に一致する。そこで、式（６０）を次のように行列の形にまとめる。
Ｈp・ｐ＋Ｈx・Δｘ＝ＨL・ΔL  （６１）
各記号は次のように定義される。
Ｈp＝［Ｃh,px  Ｃh,py］  （６２）
【数３６】
2015141159.tif 000038
ＨL＝［Ｃl,Lx  Ｃl,Ly  Ｃl,Lz］  （６４）
【００５６】
式（６１）によって、地球表面上の任意の地点ＬＭと画像センサ受光面の画素位置ｐとの関係が、推定パラメータΔｘ、ΔLを含んだ形式で得られる。そして、地球表面上の任意の地点を観測すべき地点に固定すれば、最終的に、観測量である画素位置ｐ（ランドマークを撮像した画像と基準ランドマーク画像とのずれ）と、推定パラメータΔｘおよびΔLとの間の関係式が得られる。式（６１）を観測モデル式という。
【００５７】
（４）理想状態での観測地点に対するミラー回転角
いかなる誤差もない理想状態の場合、すなわち、アライメント誤差、姿勢決定誤差および軌道決定誤差がゼロの場合を考える。この場合において、任意の地点ＬＭを決めたとき、その地点を指向するためのミラー回転角θEW、θNSを代数計算する。式（３３）に対して、理想状態Δｘ＝０、ｐ＝０、ΔL＝０を代入して、
ｈ（θEW，θNS，０，０）＝ｌ（ＬＭ，０）  （６５）
このとき、式（６５）は以下の式（６６）のような構成となる。
【数３７】
2015141159.tif 000039
【００５８】
これを、ミラー回転角θEW、θNSについて解いて、式（６７）を得る。
θEW＝（１／２）sin-1（ａ［１］）
θNS＝（１／２）sin-1（−ａ［２］／cos（２θEW））    （６７）
ただし、ａは、開口部の中心から地表面上の観測地点への方向ベクトル
ａ＝normalize（Ｙ［ＬＭx  ＬＭy  Ｒse］T）    （６８）
である。ａ[１]は、ベクトルａの第１成分値、ａ[２]は、ベクトルａの第２成分値を表す。式（６７）は、誤差のない理想状態における、任意の地点ＬＭを指向するためのミラー回転角θEW、θNSである。
【００５９】
（５）推定
画素位置と推定パラメータとの間の線形化関係で得られた観測モデルの式（６１）を、観測方程式の一般形へ変形し、観測方程式について推定パラメータを推定する。以下、推定の方法を説明する。
【００６０】
（５．１）観測方程式の一般形
得られた画像ナビゲーションのための観測モデル式（６１）は、ベクトルに３成分を持つ３つの方程式になる。
Ｙｚout＝Ｈp・ｐ＋Ｈx・Δｘ＝ＨL・ΔL    （６９）
Ｙｚoutは、長さ１の方向ベクトルなので、Ｙｚoutの３成分のうち１成分は冗長である。すなわち、他の２成分により残りの１成分は決定される。従って、Ｙｚoutのｘ方向、ｙ方向、ｚ方向のうち、ｘ方向、ｙ方向成分の方程式のみを使用する。Ｍxyを次のように定義する。
【数３８】
2015141159.tif 000040
【００６１】
Ｍxyを用いて観測モデル式を次のように置く。
Ｍxy・Ｈpｐ＋Ｍxy・Ｈx・Δｘ＝Ｍxy・ＨL・ΔL    （７１）
式（７１）を次のように観測方程式の一般形へ変形する。
ｐ＝（Ｍxy・Ｈp）-1（−Ｍxy・Ｈx・Δｘ＋Ｍxy・ＨL・ΔL）    （７２）
実際の観測量ｙにはノイズが含まれるので、このノイズ成分をνと置き、ｐ＋νを新たな変数ｙとする。そして、式（７２）を次のように変形する。
ｙ＝Ｈ・ｘ＋ν    （７３）
各記号は以下の式（７４）で定義される。
【数３９】
2015141159.tif 000041
【００６２】
Ｎ個のランドマーク地点を観測した場合に、観測した地点をｐ1、．．．、ｐNと置き、その他の記号もそれぞれの地点で定義して、以下の式（７５）を得る。
ｙ1＝ｐ1＋ν1、  Ｈ1＝（ＭxyＨp,1）-1［−ＭxyＨx,1  ＭxyＨL,1］
ｙ2＝ｐ2＋ν2、  Ｈ2＝（ＭxyＨp,2）-1［−ＭxyＨx,2  ＭxyＨL,2］
・・・
ｙN＝ｐN＋νN、  ＨN＝（ＭxyＨp,N）-1［−ＭxyＨx,N  ＭxyＨL,N］
（７５）
【００６３】
さらに、これをまとめて以下の式（７６）および式（７７）のように記述する。
【数４０】
2015141159.tif 000042
ｂy＝ｂＨ・ｘ＋ｂν    （７７）
以下、式（７７）として定式化した観測方程式の一般形を対象として、推定則を説明する。
【００６４】
（５．２）従来の最小二乗法における課題
画像ナビゲーションでは、推定すべきパラメータの数が極めて多くなる性質がある。図３に示す搭載光学機器の構成配置においては、画像ナビゲーションの性能に影響を与えるものとして、式（６３）のΔｘと式（６４）のΔLに示される１６個のパラメータが存在する。しかし、一般にこれら１６個のパラメータすべてを、精度良く推定することは不可能である。式（７７）のような線形の観測方程式で、観測ノイズがガウス分布、且つ、ホワイトノイズである場合は、通常、最適推定アルゴリズムである最小二乗法を用いて、パラメータを推定する。式（７７）に対して最小二乗法を適用し、これら１６個のパラメータを推定すると、推定結果ｘ’は、真値ｘと比べて大幅に異なる値が求まったり、最小二乗法の中にある逆行列計算が失敗することが多い。これは、１６個のパラメータのうち可観測性が極めて低いパラメータが含まれていることに起因する。ここで言う可観測性とは、各推定パラメータの変化が観測量に現れる度合い、つまり、各推定パラメータの推定のし易さのことである。推定パラメータの可観測性を低くさせる要因には、次の２つの場合がある。
・Ｓ／Ｎ比が小さい場合
・推定パラメータの間で従属性が強い場合
【００６５】
Ｓ／Ｎ比が小さい場合は、観測量の真値の大きさに対し、ノイズの大きさが大きくなりすぎ、真値がノイズに埋もれて推定がしに難くなる状態である。式（７７）で表せば、観測量の真値ｂＨ・ｘの大きさに対し、観測ノイズｂνの大きさが大きい成分が存在することである。これは、観測行列ｂＨの大きさ、真値ｘの大きさ、観測ノイズｂνの大きさに依存し、観測方程式の特徴により決定される。
【００６６】
推定パラメータの間で従属性が強い場合は、多数ある推定パラメータのうちの、いくつか複数の推定パラメータの間において、ある推定パラメータが他の推定パラメータの線形結合形式でほぼ表現できる状態である。このとき、推定パラメータがこの線形結合の方程式を満たす値であれば、どんな値であっても、ほぼ観測方程式（７７）を満たすことになる。観測方程式（７７）を満たしても、推定パラメータに自由度が余る状態であり、式（７７）だけでは推定パラメータを一意に決定し難い状態である。
【００６７】
（５．３）従来の推定則における課題
単純に式（７７）に対して最小二乗法を適用するだけでは、上述の課題がある。従来の画像ナビゲーションにおけるパラメータ推定則では、前述の最小二乗法が持つ課題を回避するために、多数ある推定パラメータのうち、５つ前後を選択してパラメータ数を減らし、選択したパラメータに対して推定を行う。この際、パラメータの選択方法において、体系的に統一した手順はなく、多くのシミュレーション反復作業による試行錯誤で、最終的な性能が良くなる推定パラメータを選択している。そのため、設計時において、パラメータ選択に手間が掛かる。さらに、パラメータを試行錯誤で選択するので、選択したパラメータの組み合わせが、最適なものかどうかが理論的に分らない。また、設計時に一度選択したら、画像ナビゲーションでは、その組み合わせを固定して使用し続け、観測条件が変化しても不変のままである。
【００６８】
画像ナビゲーションでは通常、ランドマーク地点周辺の日陰状況、あるいは、雲の分布状況に応じて、毎回異なる複数のランドマーク地点を指向しなくてはならない。また、地球全域に一様ではなく、ある領域に偏った複数のランドマーク地点を指向する可能性もある。そのため、観測地点の数や観測地点の分布は、毎回大きく変化する可能性がある（例えば、１０〜１００地点の間）。このように、画像ナビゲーションでは、その都度、観測方程式の特性が変化し、推定パラメータの可観測性も変化する。そのため、選択すべきパラメータの最適な組み合わせは、これに合わせて変化しなくてはならない。以下、上述のような課題を解決する方法を説明する。
【００６９】
（５．４）多数の推定パラメータに対する推定則
式（７７）の観測方程式に対して、ｂｙを観測した場合に推定パラメータｘを推定する具体的な計算方法を、これより以降で提供する。推定パラメータｘは、その個数、および、可観測性に何ら制限はなく、パラメータの個数が１以上のいくつであっても、また、各パラメータの可観測性が高い、低い、あるいは、不可観測である条件に関わらず、ここで提供する計算手順は実施可能である。
【００７０】
ｂｙを観測した場合、それに対応する観測行列ｂＨが決定される。画像ナビゲーションにおいて、ｂＨには、式（６２）〜（６４）、式（７５）および式（７６）から、２次鏡のミラー回転角θEW、１次鏡のミラー回転角θNS、観測地点ＬＭ、および、その他に物性値を示す定数パラメータが含まれている。
【００７１】
ノイズｂνは、ガウス分布、且つ、ホワイトノイズであると仮定し、Ｅ[ｂνi・ｂνTｉ＋Δｉ]＝ｂＲδ（Δi）であるとする。ここで、Ｅ[・]は期待値を計算する演算子、δ（・）は離散のデルタ関数である。すなわち、ｂνの共分散をｂＲとする。なお、実際にノイズｂνが、厳密にガウス分布、且つ、ホワイトノイズではない場合においても、特性が似ているノイズであれば、それによる推定性能の低下は小さいことが分かっている。
【００７２】
観測行列ｂＨとノイズｂνの共分散ｂＲが、推定パラメータｘの可観測性を決定する。ここで、次の式（７８）で行列ｂＨＨを数値的に計算する。
ｂＨＨ＝ｂＨT・ｂＲ-1・ｂＨ  （７８）
推定パラメータｘの次元数をｎxとすると、ｂＨＨはＲｎx×ｎxの正方対称行列である。ｂＨＨは、推定パラメータｘの可観測性を表す情報が含まれている。ここでは、行列ｂＨＨを観測特性行列という。
【００７３】
次に、観測特性行列ｂＨＨの固有値対角行列ｂＤＨＨ、および、固有ベクトル行列ｂＶＨＨを数値的に計算する。あるいは、ｂＨの特異値対角行列を２乗した行列、および、特異ベクトル行列を計算しても同じ結果を得る。固有値対角行列ｂＤＨＨの対角成分は、値の大きい順にソートされているものとする。固有ベクトル行列ｂＶＨＨは、直交行列である。従って、その逆行列はｂＶＨＨTになる。固有値分解では、一般に、以下の性質が成立する。一般に固有値は、Ａν＝λνを満たすλ、νで複数存在するので、Ａｖ＝λiｖi、ｉ＝１〜ｋと記述できる。λiｖiを集めて１つの行列等式にすれば、Ａ・［...ｖi...］＝［...ｖi...］・diag（...λi...）（diagは対角行列を表す。）←→ＡＶ＝ＶＤ。よってＡ＝ＶＤＶTと表せる。これを、観測特性行列ｂＨＨに適用すれば、
ｂＨＨ＝ｂＶＨＨ・ｂＤＨＨ・ｂＶＨＨT    （７９）
である。固有値対角行列ｂＤＨＨの対角成分には、観測特性行列ｂＨＨの固有値が大きいものから順に並んでいる。
【数４１】
2015141159.tif 000043
【００７４】
式（７９）の固有値分解によって、パラメータ空間が、推定パラメータｘの空間から、パラメータ推定方向の主成分空間σｘに変換される。そこでは、観測特性行列ｂＨＨの固有値と推定パラメータｘが１対１で対応しているのではなく、推定パラメータｘの線形結合が観測特性行列ｂＨＨの固有値に対応している。その対応関係は、固有ベクトル行列ｂＶＨＨを用いて、
σｘ＝ｂＶＨＨT・ｘ    （８１）
で表される。これは、主成分空間σｘのｉ番目成分は、ｉ番目固有ベクトルを係数とする、推定パラメータｘの線形結合で表される、ということを示す。すなわち、主成分空間σｘのｉ番目成分に対する推定パラメータｘの配合比率は、ｉ番目固有ベクトルの各成分値となる。この変換により、観測方程式は、以下の式（８２）で、主成分空間の観測行列σＨと推定パラメータσｘに変換される。
ｂｙ＝ｂＨ・ｘ＋ｂν＝ｂＨ・ｂＶＨＨ・ｂＶＨＨT・ｘ＋ｂν
＝σＨ・σｘ＋ｂν                            （８２）
ただし、σＨ＝ｂＨ・ｂＶＨＨである。
【００７５】
ここで、観測特性行列ｂＨＨのｎx個ある固有値の大きさ、すなわち、固有値対角行列ｂＤＨＨが、推定パラメータσｘの可観測性の度合いを表している。ｂＤＨＨにあるｉ番目の固有値は、推定パラメータσｘのｉ番目の成分値における可観測性に対応している。それぞれの固有値が大きければ大きいほど、それに対応するσｘの成分値の可観測性が高いことを意味する。すなわち、それに対応するσｘの成分値を精度良く推定可能であることを示す。値の大きな固有値が多く含まれれば、それだけ多くのパラメータを精度良く推定できることになる。逆に、値の小さな固有値が含まれれば、そのパラメータの推定精度は低くなり、推定結果は大きく上下に変化する可能性がある。固有値がゼロの場合は、不可観測の状態で、そのパラメータは全く推定できないことを示す。このように、ｂＤＨＨを見ることで、何個のパラメータは推定精度が高く、何個のパラメータは推定精度が低く、何個のパラメータは不可観測の状態なのかが、非常に簡単に分かる。
【００７６】
ゼロ、あるいは、相対的に値の小さな固有値が含まれていると、推定結果を劣化させる。そのため、これらの固有値に対応する推定パラメータを削除して、推定するパラメータの数を低減させる。ｂＤＨＨの対角上には、固有値が大きいものから順に並んでいるので、上位の何個目までをパラメータとして選定するかを決定すればよい。ここでは、例えば、ｎp個を選定すると仮定する。選定すべき固有値の大きさは、観測方程式の特徴によって変化する。例えば、実施の形態の画像ナビゲーションにおける観測方程式の場合、一般的な条件によるシミュレーションの結果から、値が４・１０８以上である固有値を選定すればよいことが分かっている。これを基に固有値を選定すると、例えば、ランドマークを２０地点観測した場合では、通常、ｎp＝５、すなわち、５個の固有値、１００地点観測した場合ではｎp＝７で、７個の固有値、１０００地点観測した場合ではｎp＝９で、９個の固有値を選定するのが最適となる。
【００７７】
なお、条件が異なる場合であっても、選定する固有値の個数をｎp＝１から順番に１個ずつ推定パラメータｘの次元数ｎxまで増やしてシミュレーションを実行し、最終的な推定性能が最も良くなるｎpの値を調べることが可能である。つまり、調べるケース数は、最大、推定パラメータｘの次元数ｎxまでである。実施の形態の場合、ｎx＝１６であるから、調べるケース数は、ｎp＝１から１６までの最大１６ケースである。本実施の形態では、このように少ないケース数を調べるだけで、パラメータの選定を決定することができる。
【００７８】
一方、総当たりでパラメータを選定する場合は、パラメータを５個選定すると決めたとき、１６個あるパラメータから５個の組み合わせを選ぶので、最大１６Ｃ５＝４３６８ケースを調べる必要がある。パラメータを６個選定すると決めたとき、１６Ｃ６＝８００８ケース、７個選定すると決めたとき、１６Ｃ７＝１１４４０ケースを調べる必要がある。何個のパラメータを選定すればよいかが分らないので、すべてを調べる必要があり、総数は、最大１６Ｃ１＋１６Ｃ２＋...＋１６Ｃ１６＝６５５３５となり、非常に膨大なケース数となる。それに対して、本実施の形態に従えば、最大１６ケースを調べるだけで選定可能である。
【００７９】
値の大きなものから順にｎp個の固有値を選定した後、次に、対応する固有ベクトル行列ｂＶＨＨを決定する。選定した固有値に対応する固有ベクトルのみを持つ固有ベクトル行列ｂ＊ＶＨＨを新たに生成する。上位のｎp個目までの固有値を選定するので、固有ベクトル行列ｂＶＨＨの１列目からｎp列目までを抜き取り、限定固有ベクトル行列ｂ＊ＶＨＨとする。式で表せば、
ｂ＊ＶＨＨ＝ｂＶＨＨ・＊Ｍ      （８３）
＊Ｍは列の抜き取り操作を行うＲｎx×ｎpの行列で、次の式（８４）である。
【数４２】
2015141159.tif 000044
【００８０】
対応する固有値行列もｂ＊ＤＨＨ＝＊ＭT・ｂＤＨＨ・＊Ｍと置く。パラメータの数を低減させた場合の推定パラメータσ＊ｘは、
σ＊ｘ＝ｂ＊ＶＨＨT・ｘ    （８５）
である。すると、観測方程式（８２）は、次の観測方程式（８６）となり、パラメータの数を低減させた場合の主成分空間の観測行列σ＊Ｈと推定パラメータσ＊ｘに変換される。
ｂｙ＝σ＊Ｈ・σ＊ｘ＋ｂν    （８６）
ただし、σ＊Ｈ＝ｂＨ・ｂ＊ＶＨＨである。
【００８１】
次に、観測方程式（８６）に基づいて推定則を構築する。元々の観測方程式（７７）ｂｙ＝ｂＨ・ｘ＋ｂνにおける最適推定アルゴリズムは、最小二乗法である。パラメータの数を低減させた場合の観測方程式（８６）に対しても、同様に最小二乗法を用いて、パラメータを推定すると、次の式（８７）を得る。
σ＊ｘ'＝σ＊Ｑ・σ＊ＨT・ｂＲ-1（ｂｙ−ｂＨ・ｘ'0）    （８７）
σ＊ｘ'は、パラメータの数を低減させた場合の主成分空間のパラメータσ＊ｘに対する推定結果である。ｘ'0は、前回の推定により求められた、元空間のパラメータｘに対する推定結果である。推定結果σ＊ｘの推定誤差共分散σ＊Ｑは、次の式（８８）で表される。
σ＊Ｑ＝（σ＊Ｑ0-1＋σ＊ＨT・ｂＲ-1・σ＊Ｈ）-1      （８８）
σ＊Ｑ0は、パラメータσ＊ｘの事前情報である。
【００８２】
式（８８）の推定則において、ｂ＊ＤＨＨを用いた表現に変形すると、式（８９）で表せる。
σ＊Ｑ＝(σ＊Ｑ0-1＋ｂ＊ＶＨＨT・ｂＨT・ｂＲ-1・ｂＨ・ｂ＊ＶＨＨ)-1
＝(σ＊Ｑ0-1＋ｂ＊ＶＨＨT・ｂ＊ＶＨＨ・ｂ＊ＤＨＨ・ｂ＊ＶＨＨT・ｂ＊ＶＨＨ)-1
＝(σ＊Ｑ０-1＋ｂ＊ＤＨＨ)-1                                    （８９）
ｂ＊ＤＨＨがσ＊Ｑ0-1に加算され、その逆行列が推定誤差共分散σ＊Ｑになることから、ｂ＊ＤＨＨの値の大きさが、直接、推定に寄与することが分かる。ｂＤＨＨのうち、値が小さい固有値を削除してｂ＊ＤＨＨとしたが、値が小さいため、削除前のσ＊Ｑ0-1と削除後のσ＊Ｑ0-1＋ｂ＊ＤＨＨとの値の大きさはほとんど変化せず、従って、推定誤差共分散σ＊Ｑにはほとんど影響を与えない。すなわち、値が小さい固有値を削除することでパラメータの数を低減しても、推定性能にはほとんど影響を与えないことが示される。なお、数値計算には、計算量の少ない、式（８９）の表現を用いることができる。
【００８３】
次に、主成分空間のパラメータσ＊ｘに対する推定結果σ＊ｘ'を、元空間の推定パラメータｘに対する推定結果ｘ'に変換する。変換は、次の関係式（９０）により行う。
ｘ'＝ｂ＊ＶＨＨ（ｂ＊ＶＨＨT・ｂ＊ＶＨＨ）-1・σ＊ｘ'      （９０）
式（９０）によって、式（８５）の関係を満たすｘ'のうち、ｘ'の２乗ノルムを最小にするものが求められる。これにより、パラメータｘの推定結果ｘ'を得る。元の推定パラメータｘの空間に対する推定誤差共分散Ｑは、以下の式（９１）で与えられる。
Ｑ＝ｂ＊ＶＨＨ（ｂ＊ＶＨＨT・ｂ＊ＶＨＨ）-1・σ＊Ｑ（ｂ＊ＶＨＨT・ｂ＊ＶＨＨ）-1・ｂ＊ＶＨＨT    （９１）
得られた推定誤差共分散Ｑは、次回の推定の際に、事前情報ｘ0'＝ｘ'、Ｑ0＝Ｑとして利用することができる。ただし、Ｑ0からσ＊Ｑ0への変換は、式（８５）から、次の関係式（９２）により行う。
σ＊Ｑ0＝ｂ＊ＶＨＨT・Ｑ0・ｂ＊ＶＨＨ    （９２）
【００８４】
（６）補正
以上で推定された誤差の推定パラメータを用いて、撮像装置の視線方向を補正することができる。パラメータの推定から従属的に導出される補正方法である。これは、画像ナビゲーションを行うに当たって、パラメータを推定した後に、次の観測地点をより精度よく指向するための処理である。
【００８５】
仮に、いかなる誤差もない理想状態の場合、すなわち、アライメント誤差、姿勢決定誤差、軌道決定誤差がゼロの場合は、目標の観測地点を決定した場合に、式（６７）で求まるミラー回転角θEW、θNSを用いて、撮像装置の１次鏡、および、２次鏡を回転させれば、完全に観測地点を指向することができる。そして、観測地点を中心とした正確な画像を撮像することができる。しかし、現実には必ず誤差が存在するので、求めたミラー回転角θEW、θNSを用いて観測地点を指向すると、目標の観測地点からずれた地点を指向することになる。そこで、前述の推定則により、アライメント誤差、姿勢決定誤差および軌道決定誤差をパラメータｘとして推定を行い、前述の推定結果ｘ'を得る。以下に、この推定結果ｘ'を考慮に入れて、目標の観測地点を指向するためのミラー回転角θEW,c、θNS,cを計算する方法を記載する。
【００８６】
（７）補正後ミラー回転角
推定結果ｘ'を考慮に入れてミラー回転角を計算すると、式（６７）で求められる誤差のない理想状態のミラー回転角θEW、θNSから値が変化する。この変化量をミラー回転角の補正量ΔEW、ΔNSと置き、補正後ミラー回転角θEW,c、θNS,cを、次の式（９３）のように表現する。
θEW,c＝θEW＋ΔEW
θNS,c＝θNS＋ΔNS    （９３）
【００８７】
目標の観測地点ＬＭを決定した場合における、補正後ミラー回転角θEW,c、θNS,cを計算するには、以下の方法を用いることができる。非線形の場合の観測モデル式（３３）において、式（９３）の関係、および、観測地点ＬＭ、推定により求めたパラメータｘ'、また、目標の観測地点が画像の中心に位置するためにｐ＝０、のそれぞれを代入して、式（９４）を得る。
ｈ（θEW＋ΔEW，θNS＋ΔNS，Δｘ'，０）＝ｌ（ＬＭ，ΔL'）    （９４）
この関係を満たすΔEW、ΔNSが、ミラー回転角の補正量となる。ここで、Δｘ'とΔL'は、定義式（７４）より、推定結果ｘ'の対応する部分である。
【数４３】
2015141159.tif 000045
ミラー回転角の補正量ΔEW、ΔNSは、式（９６）で計算できる。
ΔEW,ΔNS＝solveΔEW,ΔNS｛ｈ（θEW＋ΔEW，θNS＋ΔNS，Δｘ'，０）
＝ｌ（ＬＭ，ΔL'）                                    （９６）
式（９６）において、非線形方程式からΔEW、ΔNSを解く操作solveΔEWΔNS｛・｝には、一般的なニュートン法を用いて、実時間で高速に、数値的に解ΔEW、ΔNSを計算することができる。解ΔEW、ΔNSは小さな値となるので、ニュートン法の初期推定解ΔEW,0、ΔNS,0の値は、例えばゼロとすればよい。補正後ミラー回転角θEW,c、θNS,cは、式（９６）よりΔEW、ΔNSを求めることで、式（９３）より直ちに計算できる。
【００８８】
ミラー回転角の補正量ΔEW、ΔNSを用いて、目標の観測地点ＬＭを指向するために、撮像装置の１次鏡、および、２次鏡を、補正後ミラー回転角θEW,c、θNS,cの値に回転させることで、より精度よく観測地点を指向することができる。
【００８９】
実施の形態１．
図５は、本発明の実施の形態１に係る画像ナビゲーション装置の構成例を示すブロック図である。実施の形態１では、上述の画像ナビゲーションのうち、推定までを行う。画像ナビゲーション装置３０は、比較部３１、観測行列算出部３２、観測特性行列算出部３３、固有値算出部３４、選択部３５、限定部３６、推定部３７および変換部３８から構成される。撮像装置４２、記憶部４０および画像ナビゲーション装置３０は、衛星４１に搭載される（図１参照）。撮像装置４２の構成は決まっているので、例えば図３に基づいて、画像ナビゲーション装置３０に対して推定パラメータと座標変換行列は設定されている。
【００９０】
比較部３１は、地球表面のランドマークを撮像した画像と記憶部４０から読み込んだ基準ランドマーク画像とを比較し、それらのずれである観測量を観測行列算出部３２に渡す。観測行列算出部３２は、式（７７）の観測行列ｂＨを計算する。観測特性行列算出部３３は、観測行列ｂＨから、式（７８）の観測特性行列ｂＨＨを算出する。固有値算出部３４は、式（７９）を満たす観測特性行列ｂＨＨの固有値対角行列ｂＤＨＨ、および、固有ベクトル行列ｂＶＨＨを算出する。
【００９１】
選択部３５は、固有値対角行列ｂＤＨＨに含まれる固有値の大きいものから順に、定めた基準でｎp個の限定固有値を選択する。限定部３６は、固有ベクトル行列ｂＶＨＨから、限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列ｂ＊ＶＨＨを生成する。また、限定固有ベクトル行列ｂ＊ＶＨＨに対応する固有値行列ｂ＊ＤＨＨを生成する。
【００９２】
推定部３７は、限定固有ベクトル行列ｂ＊ＶＨＨを用いて、限定固有値に対応して推定パラメータの数を低減させた主成分空間に線形化関係式を変換した観測方程式（８２）に含まれるパラメータを、最小二乗法を用いて推定する。そして、式（８７）から、主成分空間のパラメータσ＊ｘに対する推定結果σ＊ｘ'を計算する。変換部３８は、関係式（９０）により、主成分空間のパラメータσ＊ｘに対する推定結果σ＊ｘ'を、元空間の推定パラメータｘに対する推定結果ｘ'に変換する。
【００９３】
なお、式（９０）において、主成分空間のパラメータσ＊ｘに対する推定結果σ＊ｘ'を、元空間のパラメータｘに対する推定結果ｘ'に変換する際、次の関係式により行ってもよい。
ｘ'＝ｂ＊ＶＨＨ・σ＊ｘ'    （９７）
これは、式（８５）の関係を満たすｘ'のうち、σｘ'からσ＊ｘ'へ低減した推定パラメータをゼロと置いた場合のｘ'が求まる。すなわち、次の関係式と等価である。
ｘ'＝ｂＶＨＨ・σｘ'，  σｘ'＝［σ＊ｘ'  ０ ... ０］T    （９８）
これにより、パラメータｘの推定結果ｘ'を得る。元の推定パラメータｘの空間に対する推定誤差共分散Ｑは、
Ｑ＝ｂ＊ＶＨＨ・σ＊Ｑ・ｂ＊ＶＨＨT    （９９）
で得られる。得られた推定誤差共分散Ｑは、次回の推定の際に、式（８７）および式（８８）の事前情報ｘ'0＝ｘ'、Ｑ0＝Ｑとして利用することができる。ただし、Ｑ0からσ＊Ｑ0への変換は、式（９２）と同様に、次の関係式（１００）により行う。
σ＊Ｑ0＝ｂ＊ＶＨＨT・Ｑ0・ｂ＊ＶＨＨ    （１００）
【００９４】
図６は、実施の形態１に係る推定パラメータの推定の一例を示すフローチャートである。観測行列算出部３２は、比較部３１から観測量を読み込んで（ステップＳ１１）、観測行列を計算する（ステップＳ１２）。観測特性行列算出部３３は、観測行列から観測特性行列を算出する（ステップＳ１３）。そして、固有値算出部３４は、観測特性行列の固有値対角行列と固有ベクトル行列を算出する（ステップＳ１４）。
【００９５】
選択部３５は、固有値の大きいものから順に、定めた基準で限定固有値を選択する（ステップＳ１５）。限定部３６は、固有ベクトル行列から、限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列ｂ＊ＶＨＨを生成する（ステップＳ１６）。
【００９６】
推定部３７は、限定固有値に対応して推定パラメータの数を低減させた主成分空間に線形化関係式を変換した観測方程式（８２）に含まれるパラメータを、最小二乗法を用いて推定する（ステップＳ１７）。そして、主成分空間のパラメータに対する推定結果を計算する。変換部３８は、主成分空間のパラメータに対する推定結果を、変換前の元空間の推定パラメータに対する推定結果ｘ'に変換する（ステップＳ１８）。
【００９７】
画像ナビゲーション装置３０は、コンピュータ、ＦＰＧＡ（Field Programmable Gate Array）もしくはＣＰＬＤ（Complex Programmable Logic Device）などのＰＬＤ（Programmable Logic Device）、または、行列演算などの演算素子を含む論理演算素子を用いた論理回路、などで構成することができる。
【００９８】
以上説明したように、本実施の形態１の画像ナビゲーション装置３０によれば、以下のような効果がある。第１に、衛星に対する画像ナビゲーションにおいて、アライメント誤差、姿勢決定誤差、軌道決定誤差をパラメータとして推定を行う際に、推定するパラメータの個数、および、その可観測性（各推定パラメータの推定のし易さ）に依らず、数学的な性能限界に近い高精度な推定結果を得られる。衛星に対する画像ナビゲーションの場合、衛星に存在するすべてのアライメント誤差、姿勢決定誤差、軌道決定誤差を考えると、例えば、１６個のパラメータが存在し、そのすべてを計算に取り込んで推定することで、従来より高精度な推定結果を得られる。
【００９９】
第２に、推定するパラメータの個数およびその可観測性に依らず、常に数値的に安定した計算方法および推定結果を得られる。画像ナビゲーションのような可観測性の低いパラメータが多数ある場合、単純に最小二乗法を適用すると、推定結果が激しく上下したり、真値と比べて大幅に異なる値が求まったり、最小二乗法の中にある逆行列計算が失敗することが多い。本実施の形態１の画像ナビゲーション装置３０では、これらの現象が発生せず、常に安定した計算方法および推定結果を提供する。
【０１００】
第３に、推定するパラメータの個数およびその可観測性に依らず、常に体系的に統一した推定手順を提供できる。どのような条件であっても、アルゴリズムとして本実施の形態の計算手順を同様に行うことができる。これにより、設計時におけるパラメータの選択に掛かる手間を削減できる。さらに、衛星による地球の観測条件に応じて、その都度リアルタイムで自動的にパラメータの選択を行うことができる。
【０１０１】
第４に、推定するパラメータの個数およびその可観測性に依らず、常に最適なパラメータの組み合わせによる推定手順を提供できる。従来の画像ナビゲーションでは、可観測性の低いパラメータが多数ある場合、そのうち影響の大きなパラメータ数個を試行錯誤で選択するので、選択したパラメータの組み合わせが、本当に最適なものかどうかが理論的に分らない。本実施の形態１の画像ナビゲーション装置３０では、ある行列の固有値を計算し、値の大きな固有値を選択するため、最適な選択であることが数値的に示される。
【０１０２】
以上の効果により、常に高精度な推定結果が得られることから、所定の観測地点を高精度で指向することができ、最終的な画像ナビゲーションの性能が向上する、といった効果がある。また、例えば、画像ナビゲーションで得られた推定結果ｘ'によって、本来得られるべき画像がどのようにずれたかが分かる。そこで、撮影した画像を補正することができる。
【０１０３】
実施の形態２．
図７は、本発明の実施の形態２に係る画像ナビゲーション装置の構成例を示すブロック図である。実施の形態２では、推定パラメータの推定結果に基づいて、撮像装置の制御量を補正する。実施の形態２の画像ナビゲーション装置３０は、実施の形態１の構成に、補正部３９が追加されている。その他の構成は、実施の形態１と同様である。
【０１０４】
補正部３９は、推定パラメータの推定結果ｘ'を考慮に入れて、例えば式（９６）により、撮像装置４２の制御量であるミラー回転角の補正量ΔEW、ΔNSを計算する。補正部３９は、ミラー回転角の補正量ΔEW、ΔNSから、式（９３）で補正後ミラー回転角θEW,c、θNS,cを算出し、撮像装置に入力する。撮像装置では、補正後ミラー回転角θEW,c、θNS,cの値に回転させることで、より精度よく観測地点を指向することができる。
【０１０５】
実施の形態２では、撮像装置４２として、図３の構成を想定している。その場合、撮像装置４２の制御量は、１次鏡６と２次鏡４のミラー回転角である。撮像装置４２の制御量は、ミラー回転角に限らない。例えば、２次鏡４を用いずに、画像センサ２と集光レンズ３の１次鏡６に対する角度を変化できるように構成した場合、制御量は、２次鏡４のミラー回転角および画像センサ２と集光レンズ３の回転角である。撮像装置４２全体の衛星に対する姿勢を変えられるような構成では、撮像装置の角度が制御量である。
【０１０６】
図８は、実施の形態２に係る補正の動作の一例を示すフローチャートである。図８のフローチャートは、実施の形態１の動作にステップが追加されている。ステップＳ１１〜ステップＳ１８は、実施の形態１の図６のフローチャートと同様である。変換部３８が、主成分空間のパラメータに対する推定結果を、変換前の元空間の推定パラメータに対する推定結果ｘ'に変換すると（ステップＳ１８）、補正部３９は、推定結果ｘ'を考慮して、ミラー回転角の補正量を計算し、補正後ミラー回転角を計算して、撮像装置に入力する（ステップＳ１９）。
【０１０７】
本実施の形態２の画像ナビゲーション装置３０によれば、実施の形態１の効果に加えて、以下の効果を奏する。実施の形態１の効果により、常に高精度な推定結果が得られることから、所定の観測地点を高精度で指向することができ、最終的な画像ナビゲーションの性能が向上する。また、前述の推定則により推定したパラメータの結果を用いて、所定の観測地点を指向するための補正後ミラー回転角を精度良く計算することで、所定の観測地点を高精度で指向することができる。
【０１０８】
実施の形態３．
実施の形態３では、推定で得られた推定誤差共分散を、次回の推定の際に、事前情報として利用する。画像ナビゲーションにおいて、今回の推定と次回の推定との間の経過時間が大きいなどの要因で、推定パラメータが変動する場合は、式（８９）において得られた推定誤差共分散Ｑに、変動の不確定性の大きさを表す共分散行列Ｗdを加算して、次回の推定の事前情報Ｑ0＝Ｑ＋Ｗdとして利用することができる。
【０１０９】
初回の推定、あるいは、事前情報がない場合は、式（８７）と式（８９）において、ｘ'0＝０、Ｑ0＝∞、すなわち、Ｑ0-1＝０とすることができる。パラメータが取り得る値の範囲を持つ場合は、Ｑ0にその範囲情報をガウス分布の分散として与えることができる。
【０１１０】
図９は、本発明の実施の形態３に係る推定パラメータの推定の一例を示すフローチャートである。図９のステップＳ１１〜ステップＳ１６は、実施の形態１の図６のフローチャートと同様である。選択部３５が、固有値の大きいものから順に定めた基準で限定固有値を選択し（ステップＳ１５）、限定部３６が、固有ベクトル行列から、限定固有値に対応する固有ベクトルのみを有する限定固有ベクトル行列を生成する（ステップＳ１６）と、推定部３７は、前回の推定で記憶しておいた事前情報である前回の推定誤差共分散を読み込む（ステップＳ２１）。前回の結果がない場合は、事前情報を０とする。
【０１１１】
推定部３７は、限定固有値に対応して推定パラメータの数を低減させた主成分空間に線形化関係式を変換した観測方程式（８２）に含まれるパラメータを、最小二乗法を用いて推定し（ステップＳ１７）、主成分空間のパラメータに対する推定結果を計算する。そのとき得られた推定誤差共分散に変動の不確定性の大きさを表す共分散行列Ｗdを加算して、次回の推定の事前情報として推定誤差共分散を記憶する（ステップＳ２２）。変換部３８は、主成分空間のパラメータに対する推定結果を、変換前の元空間の推定パラメータに対する推定結果ｘ'に変換する（ステップＳ１８）。
【０１１２】
実施の形態３の画像ナビゲーション装置３０によれば、推定パラメータが変動する場合、あるいは、事前情報がない場合において、より正確な推定を行うことができ、従って、より高精度な推定結果が得られる、といった効果がある。
【０１１３】
実施の形態４．
実施の形態４では、ランドマークを撮像するごとに、観測行列、観測特性行列、固有値対角行列および固有ベクトル行列を算出し、そのつど限定固有値を選択して、推定パラメータを推定する。
【０１１４】
画像ナビゲーションでは一般に、ランドマーク地点周辺の日陰状況、あるいは、雲の分布状況に応じて、毎回異なる複数のランドマーク地点を指向しなくてはならない。また、地球全域に一様ではなく、ある領域に偏った複数のランドマーク地点を指向する可能性もある。そのため、観測地点の数や観測地点の分布は、毎回大きく変化する可能性がある（例えば、１０〜１００地点の間）。このように、画像ナビゲーションでは、その都度、観測方程式の特性が変化し、推定パラメータの可観測性も変化するため、低減するパラメータの最適な組み合わせは、これに合わせて変化しなくてはならない。
【０１１５】
実施の形態１において、値の小さな固有値に対応する推定パラメータを削除して、推定するパラメータの数を低減させる方法を記載した。これによれば、事前に設計されたある閾値以上の大きさの固有値を選定すればよいので、自動的にパラメータを選定することができる。従って、画像ナビゲーションにおいて、毎回異なる複数のランドマーク地点を指向したその都度、実施の形態１に記載した推定則を実施し、その中で、固有値対角行列、および、固有ベクトル行列を計算し、ある閾値以上の大きさの固有値を選定することで、最適なパラメータの数を毎回決定し、パラメータの数を変化させてもよい。
【０１１６】
例えば、ある観測において、ランドマークを２０地点観測できた場合、ｎp＝５、すなわち、５個の固有値を選定し、次の観測において、今度は１００地点観測できた場合、ｎp＝７、７個の固有値を選定、その次の観測では、再び、２０地点しか観測できなかった場合、ｎp＝５、５個の固有値を選定する、といったように、観測条件が変化すればそれに合わせて、自動的にパラメータの数を変化させることができる。合わせて、固有ベクトル行列も毎回変化する。すなわち、主成分空間が変化する。これにより、毎回、高精度な推定結果が得られる。
【０１１７】
図１０は、本発明の実施の形態４に係る補正の動作の一例を示すフローチャートである。図１０のフローチャートは、図８のステップＳ１１〜ステップＳ１９を繰り返す構成になっている。また、選択部３５は、閾値以上の限定固有値を選択する（ステップＳ３５）。そして、限定部３６は、閾値以上の限定固有値に対応する限定固有ベクトル行列を生成する（ステップＳ３６）。その他のステップは、実施の形態２と同様である。ただし、変換部３８が補正後ミラー回転角を計算して、撮像装置に入力すると、ステップＳ１１に戻って、比較部３１から新たな観測量を読み込むことから繰り返す。
【０１１８】
本実施の形態４の画像ナビゲーション装置３０によれば、毎回異なる複数のランドマーク地点を指向したり、観測地点数が大きく変化するような、観測条件が大きく変化する場合であっても、常に高精度な推定結果が得られる、といった効果がある。また、自動的にパラメータを選定するため、観測条件が変化しても、手間が掛からない、といった効果がある。
【符号の説明】
【０１１９】
１  光軸、２  画像センサ、３  集光レンズ、４  ２次鏡、５  ２次鏡回転軸、６  １次鏡、７  １次鏡回転軸、８  開口部、２１  地球、３０  画像ナビゲーション装置、３１  比較部、３２  観測行列算出部、３３  観測特性行列算出部、３４  固有値算出部、３５  選択部、３６  限定部、３７  推定部、３８  変換部、３９  補正部、４０  記憶部、４１  衛星、４２  撮像装置、４３  ２次鏡、４４  １次鏡、４６  視線方向、４７  ランドマーク画像、４８  基準ランドマーク画像、４９  相関部、５０  推定部、５１  補正部、５２  補正観測部、６４  観測部、６５  ランドマーク画像、６６  参照部、６７  基準ランドマーク画像、７５  地球。
【図１】
2015141159.tif 000046
【図２】
2015141159.tif 000047
【図３】
2015141159.tif 000048
【図４】
2015141159.tif 000049
【図５】
2015141159.tif 000050
【図６】
2015141159.tif 000051
【図７】
2015141159.tif 000052
【図８】
2015141159.tif 000053
【図９】
2015141159.tif 000054
【図１０】
2015141159.tif 000055