(57)【特許請求の範囲】
【請求項１】
少なくとも２機の人工衛星による、地上を飛翔する飛翔体を撮影し撮影画像データとして送信する段階と、
前記少なくとも２機の人工衛星による、各々の人工衛星の位置情報、速度情報、及び姿勢情報を送信する段階と、
受信アンテナを備えた基地局による、前記少なくとも２機の人工衛星より送信される少なくとも２枚の画像データと、前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とを受信する段階と、
前記基地局における位置解析装置による、前記アンテナにより受信した少なくとも２枚の画像データを解析して前記各々の人工衛星の頂角及び方位角を算出する段階と、
前記頂角及び方位角と、前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とを用いて前記飛翔体の位置を解析する段階と
を含む飛翔体探知方法であって、
前記各々の人工衛星の頂角及び方位角をそれぞれ
0005567805.tif 000135
としたとき、次式によって定義される観測ベクトルｙを生成し、
0005567805.tif 000136
一方で、前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とから次式に基づいて座標変換行列を生成し、
0005567805.tif 000137
一方で、飛翔体位置初期値をＸｒｅｆとしたとき、
前記飛翔体位置初期値と前記人工衛星の位置情報と前記座標変換行列（式４）、（式５）、（式７）とから次式で定められる係数を生成し、
0005567805.tif 000138
さらに、前記（式１１ａ）〜（式１１ｃ）により定められる係数と、前記観測ベクトルｙとから次式で定められるパラメータ推定処理
0005567805.tif 000139
を実行することによって飛翔体位置を特定する
ことを特徴とする方法。
【請求項２】
地上を飛翔する飛翔体を撮影可能な観測器と、撮影した画像データと各々の位置情報、速度情報、及び姿勢情報とを送信する送信器とを備えた少なくとも２機の人工衛星と、
前記少なくとも２機の人工衛星より送信される少なくとも２枚の画像データと前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とを受信するためのアンテナと、前記アンテナにより受信した前記少なくとも２枚の画像データと前記各々の位置情報、速度情報、及び姿勢情報とに基づいて前記飛翔体の位置を解析する位置解析装置とを備えた基地局と
で構成される飛翔体探知システムであって、
前記位置解析装置は、
前記少なくとも２枚の画像データを解析して前記各々の人工衛星の頂角及び方位角を算出する算出部と、
前記各々の人工衛星の頂角及び方位角をそれぞれ
0005567805.tif 000140
としたとき、次式によって定義される観測ベクトル
0005567805.tif 000141
を生成する観測ベクトル生成部と、
前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とから次式に基づいて座標変換行列
0005567805.tif 000142
を生成する座標変換行列生成部と、
飛翔体位置初期値をＸｒｅｆとしたとき、
前記飛翔体位置初期値と前記人工衛星の位置情報と前記座標変換行列（式４）、（式５）、（式７）とから次式で定められる係数
0005567805.tif 000143
を生成する係数生成部と、
前記（式１１ａ）〜（式１１ｃ）により定められる係数と、前記観測ベクトルｙとから次式で定められるパラメータ推定処理
0005567805.tif 000144
を実行するパラメータ推定処理部と
を備えたことを特徴とするシステム。
【請求項３】
地上を飛翔する飛翔体を撮影可能な観測器と、撮影した画像データと各々の位置情報、速度情報、及び姿勢情報とを送信する送信器を備えた少なくとも２機の人工衛星と、
前記少なくとも２機の人工衛星より送信される少なくとも２枚の画像データと前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とを受信するためのアンテナと、前記アンテナにより受信した少なくとも２枚の画像データと前記各々の位置情報、速度情報、及び姿勢情報とに基づいて前記飛翔体の位置を解析する位置解析装置とを備えた基地局と
で構成される飛翔体探知システムにおける位置解析装置であって、
前記少なくとも２枚の画像データを解析して前記各々の人工衛星の頂角及び方位角を算出する算出部と、
前記各々の人工衛星の頂角及び方位角をそれぞれ
0005567805.tif 000145
としたとき、次式によって定義される観測ベクトル
0005567805.tif 000146
を生成する観測ベクトル生成部と、
前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とから次式に基づいて座標変換行列
0005567805.tif 000147
を生成する座標変換行列生成部と、
飛翔体位置初期値をＸｒｅｆとしたとき、
前記飛翔体位置初期値と前記人工衛星の位置情報と前記座標変換行列（式４）、（式５）、（７）とから次式で定められる係数
0005567805.tif 000148
を生成する係数生成部と、
前記（式１１ａ）〜（式１１ｃ）により定められる係数と、前記観測ベクトルｙとから次式で定められるパラメータ推定処理
0005567805.tif 000149
を実行するパラメータ推定処理部と
を備えたことを特徴とする装置。
【請求項４】
地上を飛翔する飛翔体を撮影可能な観測器と、撮影した画像データと各々の位置情報、速度情報、及び姿勢情報とを送信する送信器とを備えた少なくとも２機の人工衛星と、
前記少なくとも２機の人工衛星より送信される少なくとも２枚の画像データと前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とを受信するためのアンテナと、前記アンテナにより受信した前記少なくとも２枚の画像データと前記各々の位置情報、速度情報、及び姿勢情報とに基づいて前記飛翔体の位置を解析する位置解析装置とを備えた基地局と
で構成される飛翔体探知システムの位置解析装置において実行されるコンピュータプログラムであって、
前記少なくとも２枚の画像データを解析して前記各々の人工衛星の頂角及び方位角を算出する算出ステップと、
前記各々の人工衛星の頂角及び方位角をそれぞれ
0005567805.tif 000150
としたとき、次式によって定義される観測ベクトル
0005567805.tif 000151
を生成する観測ベクトル生成ステップと、
前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とから次式に基づいて座標変換行列
0005567805.tif 000152
を生成する座標変換行列生成ステップと、
飛翔体位置初期値をＸｒｅｆとしたとき、
前記飛翔体位置初期値と前記人工衛星の位置情報と前記座標変換行列（式４）、（式５）、（式７）とから次式で定められる係数
0005567805.tif 000153
を生成する係数生成ステップと、
前記（式１１ａ）〜（式１１ｃ）により定められる係数と、前記観測ベクトルｙとから次式で定められるパラメータ推定処理
0005567805.tif 000154
を実行するパラメータ推定処理ステップと
を備えたことを特徴とするコンピュータプログラム。
【請求項５】
Δｔを人工衛星に搭載された観測器における観測時間間隔としたとき、（式１２）から求められる時刻ｔにおける飛翔体の位置ベクトル
0005567805.tif 000155
と時刻ｔ-Δｔにおける飛翔体の位置ベクトル
0005567805.tif 000156
とを用いて次式
0005567805.tif 000157
により前記飛翔体の速度ベクトルを算出する段階をさらに含む請求項１記載の方法。
【請求項６】
Δｔを人工衛星に搭載された観測器における観測時間間隔としたとき、（式１２）から求められる時刻ｔにおける飛翔体の位置ベクトル
0005567805.tif 000158
と時刻ｔ-Δｔにおける飛翔体の位置ベクトル
0005567805.tif 000159
とを用いて次式
0005567805.tif 000160
により前記飛翔体の速度ベクトルを算出する速度ベクトル算出部をさらに備えたことを特徴とする請求項２記載のシステム。
【請求項７】
Δｔを人工衛星に搭載された観測器における観測時間間隔としたとき、（式１２）から求められる時刻ｔにおける飛翔体の位置ベクトル
0005567805.tif 000161
と時刻ｔ-Δｔにおける飛翔体の位置ベクトル
0005567805.tif 000162
とを用いて次式
0005567805.tif 000163
により前記飛翔体の速度ベクトルを算出する速度ベクトル算出部をさらに備えたことを特徴とする請求項３記載の装置。
【請求項８】
Δｔを人工衛星に搭載された観測器における観測時間間隔としたとき、（式１２）から求められる時刻ｔにおける飛翔体の位置ベクトル
0005567805.tif 000164
と時刻ｔ-Δｔにおける飛翔体の位置ベクトル
0005567805.tif 000165
とを用いて次式
0005567805.tif 000166
により前記飛翔体の速度ベクトルを算出する速度ベクトル算出ステップをさらに備えたことを特徴とする請求項４記載のコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、飛翔するロケット、ミサイル等を探知する方法及びシステムに関し、より詳細には、数ｋｍ〜数十ｋｍの上空を飛翔する物体（以下、「飛翔体」という。）を探知し、飛翔方向や到達地点等を推測する方法及びシステムに関する。
【背景技術】
【０００２】
近年、ロケット、ミサイル等の飛翔体は、いずれも不慮の墜落や爆発といった事故を引き起こしてしまった場合には、墜落による飛翔体それ自体についての人的損害及び物的損害にとどまらず、墜落地点や着弾点に対する被害が極めて大きくなることが予想され、重大な損害を及ぼしかねない。
【０００３】
そのため、従来、レーダシステムを搭載した車両や船舶を地上あるいは海上に展開し、飛翔体が正常に飛行しているか、あるいは、飛翔体がどこからどこに向かって飛翔しているのか等について監視するシステムの構築が行われてきた。
【０００４】
しかしながら、地上に展開した車両ないし海上に展開した船舶の行動範囲は主として国内に限られてくるために、その監視領域は、そもそも全地球的規模から見れば限定されたものとなっている。したがって、国や地域をまたがるような遠方から発射された飛翔体の探知は非常に困難であった。また、レーダシステムを搭載した車両を地上に展開した場合、地形によっては山等の障害物がその監視領域を狭めることになるという課題もあった。
【０００５】
そこで、車両や船舶を使用した飛翔体監視システムを補完すべく、人工衛星を使ったシステムを構築する試みもなされ、例えば、非特許文献１などのような早期警戒衛星システムの研究開発が行われている。
【０００６】
非特許文献１に開示された技術は、長楕円軌道を有し、観測機を搭載した２機の人工衛星を使用するものである。一般に、長楕円軌道を有する人工衛星は高高度（例えば高度４万ｋｍ）を飛翔するので、地球上の多くの領域を監視下におくことができるものと推察される。
【０００７】
図１９に、２機の人工衛星を周回させた場合（Ｓａｔｅｌｌｉｔｅ１軌道、Ｓａｔｅｌｌｉｔｅ２軌道）に監視可能な領域を示す。図１９において、Ｓａｔｅｌｌｉｔｅ１及びＳａｔｅｌｌｉｔｅ２は、交互に北半球と南半球とを行ったり来たりする周回を行っている。具体的には、１日のうち１２時間は、Ｓａｔｅｌｌｉｔｅ１が北半球を飛翔（北半球を監視）し、Ｓａｔｅｌｌｉｔｅ２が南半球を飛翔（南半球を監視）する一方で、次の１２時間は、Ｓａｔｅｌｌｉｔｅ１が南半球を飛翔（南半球を監視）し、Ｓａｔｅｌｌｉｔｅ２が北半球を飛翔（北半球を監視）する。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】http://www.lockheedmartin.com/products/SpaceBasedInfraredSystemHigh/（ロッキード・マーチン社ホームページ）
【非特許文献２】http://en.wikipedia.org/wiki/Space-Based_Infrared_System（米国ウィキペディア「Space-Based Infrared System」の項）
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、この方法では、北半球及び南半球をそれぞれ１機の人工衛星が周回するために、２機が同時に監視（ステレオ視）できる領域は、図１９のハッチング部分のみとなり領域としては狭い。つまり、１機の人工衛星の監視によって飛翔体の存在は確認できるものの、２機の人工衛星によるステレオ視によって飛翔体の３次元位置を正確に突き止めることが可能なエリア、及び飛翔方向を正確に予測することが可能なエリアは限定的である。
【００１０】
また、図２０に示すよう複数の静止衛星（Ｓａｔｅｌｌｉｔｅ３、Ｓａｔｅｌｌｉｔｅ４）を使用して広い監視領域を確保しようとする試みもある。この場合、おおむね図２０において示したハッチング部分のエリアを確保できる。
【００１１】
しかしながら、静止衛星の位置は事前に容易に特定されることから、第三者がこの静止衛星の監視領域外での飛翔体の発射を行うことは容易となる。また、静止衛星は通常赤道上を飛翔していることから、監視可能エリアの内であっても地球上の高緯度の場所は雲を通しての画像となるため位置特定の信頼性低下してしまうという問題がある。また、ステレオ視による位置決め可能エリアは上述のとおり広範囲となるが、静止衛星は赤道上に位置するため高緯度になるほど位置特定精度が悪くなるという問題が生じる。
【００１２】
あるいは、複数の低高度の人工衛星で監視する方法も考えられるが、低高度であるが故にその監視領域は極めて狭く、せいぜい５００ｋｍ四方にとどまるものであって、広い監視領域からの飛翔体探知を迅速に行おうとすると、数十機の低高度人工衛星を必要とする。更に言えば、これら低高度人工衛星によりステレオ視で観測を行って、飛翔体の３次元位置を正確に推定するためには、精度上倍の数の低高度人工衛星が必要となるなど、システム構築を考えるうえで現実的でない。
【００１３】
そこで、本願発明の目的は、かかる従来技術における諸問題を解決することである。すなわち、本願発明は、複数の人工衛星から構成される飛翔体検出システムであって、観測器の出力であるデータ（例えば、撮像画像データ）をステレオ（立体）視することにより飛翔体の３次元位置を正確に推定ないし特定する飛翔体探知方法及びシステムを提供することを目的とする。
【課題を解決するための手段】
【００１４】
本発明は、少なくとも２機の人工衛星による、地上を飛翔する飛翔体を撮影し撮影画像データとして送信する段階と、前記少なくとも２機の人工衛星による、各々の人工衛星の位置情報、速度情報、及び姿勢情報を送信する段階と、受信アンテナを備えた基地局による、前記少なくとも２機の人工衛星より送信される少なくとも２枚の画像データと、前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とを受信する段階と、前記基地局における位置解析装置による、前記アンテナにより受信した少なくとも２枚の画像データを解析して前記各々の人工衛星の頂角及び方位角を算出する段階と、前記頂角及び方位角と、前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とを用いて前記飛翔体の位置を解析する段階とを含む飛翔体探知方法であって、前記各々の人工衛星の頂角及び方位角から観測ベクトルを生成する段階と、前記各々の人工衛星の位置情報、速度情報、及び姿勢情報に基づいて座標変換行列を生成する段階と、飛翔体位置初期値をＸrefとしたとき、前記飛翔体位置初期値と前記人工衛星の位置情報と前記座標変換行列とから係数を生成する段階と、前記係数と、前記観測ベクトルとから定められるパラメータ推定処理を実行することによって飛翔体位置を特定する段階とを含むことを特徴とする。
【００１５】
また、本発明は、地上を飛翔する飛翔体を撮影可能な観測器と、撮影した画像データと各々の位置情報、速度情報、及び姿勢情報とを送信する送信器とを備えた少なくとも２機の人工衛星と、前記少なくとも２機の人工衛星より送信される少なくとも２枚の画像データと前記各々の人工衛星の位置情報、速度情報、及び姿勢情報とを受信するためのアンテナと、前記アンテナにより受信した前記少なくとも２枚の画像データと前記各々の位置情報、速度情報、及び姿勢情報とに基づいて前記飛翔体の位置を解析する位置解析装置とを備えた基地局とで構成される飛翔体探知システムであって、前記位置解析装置は、前記少なくとも２枚の画像データを解析して前記各々の人工衛星の頂角及び方位角を算出する算出部と、前記各々の人工衛星の頂角及び方位角から観測ベクトルを生成する観測ベクトル生成部と、前記各々の人工衛星の位置情報、速度情報、及び姿勢情報に基づいて座標変換行列を生成する座標変換行列生成部と、飛翔体位置初期値をＸrefとしたとき、前記飛翔体位置初期値と前記人工衛星の位置情報と前記座標変換行列とから係数を生成する係数生部と、前記係数と、前記観測ベクトルとから定められるパラメータ推定処理を実行するパラメータ推定処理部とを備えたことを特徴とする。
【発明の効果】
【００１６】
本発明にかかる飛翔体探知方法及びシステム等によれば、複数の人工衛星から送信されてくる複数の撮像データに写しだされた飛翔体を立体視的に解析することにより、従来技術に比べて正確に飛翔体の３次元位置を特定することができる。また、かかる撮像データを所定の時間間隔（Δｔ）で連続して受信及び解析等処理することにより、飛翔体の３次元位置のみならず速度ベクトルも推定することもでき、飛翔体の不慮の墜落や爆発といった事故に際しても、その飛翔体の墜落地点や着弾点を正確に予測することが可能となり、墜落地点や着弾点における被害を未然に防ぐ、ないし最小限に抑止することができる。
【図面の簡単な説明】
【００１７】
【図１】本発明の一実施形態にかかる飛翔体探知システムの構成例を説明する説明図である。
【図２】本発明の一実施形態にかかる飛翔体探知システムにおける位置解析装置の構成例を説明する説明図である。
【図３】本発明の一実施形態にかかる飛翔体探知方法及びシステムの監視領域の例を説明する説明図である。
【図４】本発明の一実施形態にかかる飛翔体探知方法及びシステムにおける人工衛星から地球を見下ろした様子を説明する説明図である。
【図５】本発明の一実施形態にかかる飛翔体探知方法及びシステムにおける衛星機体座標系を説明する説明図である。
【図６】本発明の一実施形態にかかる飛翔体探知方法及びシステムにおける衛星軌道座標系（Ｏ）及び地心方向軌道座標系（ＥＤＯ）を説明する説明図である。
【図７】本発明の一実施形態にかかる飛翔体探知方法及びシステムにおいて撮像された撮像画面と地心方向軌道座標系（ＥＤＯ）との関係を説明する説明図である。
【図８】本発明の一実施形態にかかる飛翔体探知方法及びシステムにおける２機の人工衛星から撮像された撮像画像を説明する説明図である。
【図９】本発明の一実施形態にかかる飛翔体探知方法及びシステムにおける衛星機体座標系（Ｓ）と飛翔体ベクトルとの関係を説明する説明図である。
【図１０】本発明の一実施形態にかかる飛翔体探知方法及びシステムにおける地球中心固定座標系（ＥＣＥＦ）とＮＥＤ（Ｎｏｒｔｈ−Ｅａｓｔ−Ｄｏｗｎ）座標系との関係を説明する説明図である。
【図１１】本発明の一実施形態にかかる飛翔体探知方法及びシステムにおいて撮像された撮像画像とＮＥＤ（Ｎｏｒｔｈ−Ｅａｓｔ−Ｄｏｗｎ）座標系との関係を説明する説明図である。
【図１２】本発明の他の実施形態にかかる飛翔体探知方法及びシステムにおける３機の人工衛星から撮像された撮像画像を説明する説明図である。
【図１３】本発明の一実施形態（２衛星、衛星姿勢情報有の場合）にかかる飛翔体探知システムの位置解析装置における座標系変換理論を説明する説明図である。
【図１４】本発明の他の実施形態（２衛星、衛星姿勢情報無の場合）にかかる飛翔体探知システムの位置解析装置における座標系変換理論を説明する説明図である。
【図１５】本発明の他の実施形態（３衛星の場合）にかかる飛翔体探知システムの位置解析装置における座標系変換理論を説明する説明図である。
【図１６】本発明の一実施形態（２衛星、衛星姿勢情報有の場合）にかかる飛翔体探知システムの位置解析装置における演算処理フローを説明する説明図である。
【図１７】本発明の他の実施形態（２衛星、衛星姿勢情報無の場合）にかかる飛翔体探知システムの位置解析装置における演算処理フローを説明する説明図である。
【図１８】本発明の他の実施形態（３衛星の場合）にかかる飛翔体探知システムの位置解析装置における演算処理フローを説明する説明図である。
【図１９】２機の人工衛星を使用した従来の監視システムにおける監視領域の例を説明する説明図である。
【図２０】２機の静止衛星を使用した従来の監視システムにおける監視領域の例を説明する説明図である。
【発明を実施するための形態】
【００１８】
以下、本発明にかかる飛翔体探知方法及びシステム等を実施するための形態について詳述する。
【実施例１】
【００１９】
［システム構成例］
図１に、本発明の一実施形態にかかる飛翔体探知システムの構成例を示す。システム１００は、地球を周回する２機（ないし３機）の人工衛星１０１ａ〜人工衛星１０１ｂ（１０１ｃ）と、人工衛星１０１ａ〜人工衛星１０１ｂ（１０１ｃ）から送信される、目標とする飛翔体が撮影された画像データや各人工衛星自体の位置や姿勢等の各種データを受信するためのアンテナ１０３を備えた基地局１０２とからなる。基地局１０２内には、アンテナ１０３が受信した画像データ等を解析して目標とする飛翔体の位置を解析するための位置解析装置１０５が備えられており、アンテナ１０３と位置解析装置１０５とは有線信号線１０４で接続されている。なお、有線信号線１０４は、その一部又は全部を無線インタフェース（例えば、無線ＬＡＮを含む）としてもよい。また、基地局１０２は建屋として説明したが、航空機、船舶、車両等の移動体であってもよい。あるいは小型のアンテナ１０３を搭載した携帯用の位置解析装置１０５を携行する形態であってもよい（この場合、アンテナ１０３を備えた位置解析装置１０５それ自体が基地局となる）。
【００２０】
なお、人工衛星１０１ａ〜人工衛星１０１ｃには、地球表面を観測しその画像データを撮影する観測器（不図示）と、撮影した画像データと人工衛星自身の位置情報及び姿勢情報とを地球上の基地局に送信するための送信器（不図示）とが少なくとも備わっている。
【００２１】
図２に、位置解析装置１０５の本発明の一実施形態にかかる飛翔体探知システムにおける位置解析装置の構成例を示す。位置解析装置１０５は、地球を周回する２機ないし３機の人工衛星から送信される目標とする飛翔体が撮影された画像データと、人工衛星それ自体の位置及び姿勢等の各種データとを受信するためのアンテナ１０３及び受信部２０１と、有線若しくは無線インタフェース１０４を介して接続されている。
【００２２】
また、位置解析装置１０５は、アンテナ１０３及び受信部２０１が受信した信号をデータとして装置本体に送信するための通信インタフェース１０５１と、プログラムによって様々な数値計算や論理演算、機器制御などを行う中央処理装置（ＣＰＵ）１０５２と、ＲＡＭ及び／又はＲＯＭ等で構成されるメモリ１０５３と、磁気又は光学媒体等で構成される記憶ディスク１０５４とで構成されている。本発明の実施に必要なソフトウェア及び／又はプログラムは、通常、記憶ディスク１０５５にインストールないし記憶され、ソフトウェア及び／又はプログラムの実行時には、必要に応じてメモリ１０５３内のＲＡＭ等にその全部又は一部のソフトウェアモジュールが読み出され、ＣＰＵ１０５２において演算実行される。
【００２３】
さらに、必要に応じて、特定の行列演算等を高速に処理するよう設計されたディジタルシグナルプロセッサ（ＤＳＰ）１０５４、処理結果等を表示するためのディスプレイ装置等で構成される表示部１０５６、キーボード、マウス、ペン、ＯＣＲ、光学又は磁気等のリーダ等からなる入力部１０５７を含むように構成してもよい。
【００２４】
なお、図２においては図示していないが、位置解析装置１０５内の各ユニットはバス等の信号線で適宜接続されている。また、位置解析装置１０５は、後述する本発明にかかる計算処理ないし論理演算処理を実行可能なプログラムをインストール可能であり、アンテナ１０３及び受信部２０１と適切に接続可能なインタフェースを備えたものであれば、汎用コンピュータないし携帯情報端末（携帯電話等通信機能を備えた機器を含む）等を位置解析装置１０５として用いることもできる。
【００２５】
また、本明細書においては、いくつかの実施例に基づいて本発明にかかる飛翔体探知方法及びシステム等を説明するが、数式等で例示したものを含めた計算及び演算の各々は、主として位置解析装置１０５とコンピュータプログラムとが協調動作して実行される演算処理によるものである。以下、本発明にかかる飛翔体探知方法及びシステム等が前提とする理論的技術的背景（例えば、衛星の軌道について、種々の座標系について等）を踏まえた実装上の処理について詳述する。
【００２６】
図３に、本発明の一実施形態にかかる飛翔体探知方法及びシステムにおいて実施される人工衛星の構成（及び周回軌道例）を示す。地球表面を監視するための人工衛星は３機（Ｓａｔｅｌｌｉｔｅ１〜Ｓａｔｅｌｌｉｔｅ３）で構成されており、主として北半球を監視領域とした飛翔体探知システムとなっている。なお、図においては３機の人工衛星を用いているが、２機の人工衛星のみ（例えば、Ｓａｔｅｌｌｉｔｅ１及びＳａｔｅｌｌｉｔｅ２）を使用することも可能である。あるいは図３に示した３機の人工衛星のうちの２機を適宜切り替え選択して用いるように実施することも可能である。
【００２７】
Ｓａｔｅｌｌｉｔｅ１〜Ｓａｔｅｌｌｉｔｅ３までの３機それぞれの人工衛星の軌道は、離心率が０．０７５、軌道傾斜角が４５度、昇交点赤経は１２０度ずつ、均等に配置されるものである。また、近点離角はいずれも２７０度となっている。
【００２８】
また、北半球を監視領域とするためには、北半球上空では（南半球に対して）相対的にゆっくりとした速度で飛翔するような軌道とすることが望ましい。これを踏まえて一実施形態においては、離心率０．７５、近点離角２７０度にそれぞれ設定される。ここで、軌道傾斜角を４５度に設定すると、人工衛星は緯度の高い日本上空なども飛翔するようになり、図４に示すように日本付近上空を飛翔している（図中ちょうど原点付近に衛星があるものとする）場合には北極点を越えた領域まで監視することが可能になる。
【００２９】
さらに、３機の人工衛星を使用して北半球を監視するためには、３機が順次北半球上空に飛来するように、互いの軌道の関係を調整する必要がある。この場合、一実施形態において、昇交点赤径は１２０度ずつ均等に配置されることが望ましい。
【００３０】
以上の条件を備えた構成によれば、３機の人工衛星のうちの２機によってステレオ視可能な範囲は、図３のハッチング部に示すとおりとなり広範である。
【００３１】
［地球中心固定座標系定義］
次に、地球固定座標系（ＥＣＥＦ）について説明する。地球中心固定座標系（ＥＣＥＦ）とは、地球重心を原点とし、ｚ軸は地球の地殻に固定された慣用地軸極を通り、ｘ軸は赤道面と平均グリニッジ子午線との交点を通り、ｙ軸はこれらのｘ、ｚ軸と右手系をなすような座標系である。
【００３２】
［衛星機体座標系定義］
次に、図５に基づいて、衛星機体の座標系（Ｓ）について説明する。いま、地球固定座標系原点（図中に示したｘｙｚ座標の原点）から見た衛星の位置ベクトルを
0005567805.tif 000002
とし、地球中心から見た飛翔体の位置ベクトルを
0005567805.tif 000003
とし、飛翔体から衛星へのベクトルを
0005567805.tif 000004
とする。このとき、ｒSと、衛星の位置ベクトルＲSと、飛翔体の位置ベクトルＲ0との間には次式の関係がある。
【００３３】
【数１】
0005567805.tif 000005
・・・・・・（式１）
【００３４】
なお、衛星機体の太陽電池パネル方向の単位ベクトルは
0005567805.tif 000006
衛星機体のセンサー軸方向の単位ベクトルは
0005567805.tif 000007
ＹＳとＰＳと右手系を成す単位ベクトルは
0005567805.tif 000008
飛翔体頂角は
0005567805.tif 000009
で表されている。
【００３５】
［衛星軌道座標系定義］
次に、図６（ａ）に基づいて、衛星軌道座標系（Ｏ）について説明する。いま、図６（ａ）において示す位置に楕円軌道上を周回する衛星Ｓがあるとすると、衛星Ｓにおける軌道座標系（Ｏ）は、軌道面内で速度ベクトル方向の単位ベクトル
0005567805.tif 000010
をｘ軸、上記速度ベクトル方向とそれに直交する軌道面外方向ベクトルに直交する単位ベクトル
0005567805.tif 000011
をｙ軸、及びこれらのｘ、ｙ軸と右手系をなす軸をｚ軸とする座標系となる。
【００３６】
[観測方程式の導出]
次に、図８に、２機の人工衛星に搭載された観測機器が同時に観測を行った場合の撮影画像の様子を示す。図８（ａ）は一方の人工衛星が撮影した地球の画像であり、図８（ｂ）は他方の人工衛星が撮影した地球の画像である。図８（ａ）（ｂ）における
0005567805.tif 000012
は、飛翔体を示す。原理的にこの２枚の撮影画像データを解析することによって地球に対するステレオ視（立体視）が実現されている。ここで、図８（ａ）（ｂ）に示された撮像データから直接得られる観測量は、画像中心からの角度（頂角）である
0005567805.tif 000013
と、画像の縦横を基準にしたときの画像回転角（方位角）である
0005567805.tif 000014
とからなる４つのパラメータである。このとき、図９に示すように、衛星機体座標系（Ｓ）から見た飛翔体位置は、これら２つの角度と、飛翔体から衛星へのベクトルｒSとを用いて次式のとおり表される。
【００３７】
【数２】
0005567805.tif 000015
・・・・・・（式２）
【００３８】
ここで、ベクトルｒSに付されている｜Sは、そのベクトルの基底が衛星機体座標系（Ｓ）であることを意味する。
【００３９】
そして、最終的に求めるべきベクトルは、地球の中心を原点とした地球中心固定座標系（ＥＣＥＦ）における飛翔体の３次元位置ベクトル
0005567805.tif 000016
であるが、Ｒ0は、衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）への座標変換行列ＣS2ECEFと（式１）及び（式２）との関係より、各衛星の撮像画像における観測量
0005567805.tif 000017
に基づいて、次式で示される観測方程式を解くことで求められる。
【００４０】
【数３】
0005567805.tif 000018
・・・・・・（式３）
【００４１】
＜２衛星、衛星姿勢情報有の場合の座標系変換＞
ここで、衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）への座標変換行列ＣS2ECEFは、図１３に示されるように衛星軌道座標系（Ｏ）を介して行われる。
【００４２】
各衛星の既知のロール姿勢角
0005567805.tif 000019
ピッチ姿勢角
0005567805.tif 000020
ヨー姿勢角
0005567805.tif 000021
に基づく衛星機体座標系（Ｓ）から軌道座標系（Ｏ）への座標変換行列
0005567805.tif 000022
及び各衛星の既知の３次元ベクトルＲS（ｓ＝１，２）と速度ベクトルｖS（ｓ＝１，２）に基づく軌道座標系（Ｏ）から地球中心固定座標系（ＥＣＥＦ）への座標変換行列
0005567805.tif 000023
により衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）への座標変換行列ＣS2ECEFは次式のとおり決定される。
【００４３】
【数４】
0005567805.tif 000024
・・・・・・（式４）
【００４４】
【数５】
0005567805.tif 000025
・・・・・・（式５）
【００４５】
ここで、Ｃ1（θR）、Ｃ2（θP）、Ｃ3（θY）は、それぞれ、
【００４６】
【数６】
0005567805.tif 000026
・・・・・・（式６a）
0005567805.tif 000027
・・・・・・（式６ｂ）
0005567805.tif 000028
・・・・・・（式６ｃ）
で表される。
【００４７】
そして、ＣO2ECEFは、
【００４８】
【数７】
0005567805.tif 000029
・・・・・・（式７）
で表され、ｅsxO，ｅsyO，ｅszOについては、それぞれ、
【００４９】
【数８】
0005567805.tif 000030
・・・・・・（式８a）
0005567805.tif 000031
・・・・・・（式８ｂ）
0005567805.tif 000032
・・・・・・（式８ｃ）
【００５０】
で表される。演算子×は、ベクトル積を意味する。
【００５１】
＜２衛星、衛星姿勢情報有の場合の演算処理フロー＞
さらに、上記理論と計算機の特性とをより具体的に踏まえた演算処理例について詳説する。すわなち、観測方程式（式３）は、非線形の方程式となっているため、演算解を得るためには、変数ｘに関するｘref（飛翔体位置の初期値とする）周りの１次テイラー級数近似による線形化を行い、繰り返し演算処理を実行することで解を得ることが現実的である。
【００５２】
[一般的な非線形方程式の解法]
いま、（式９ａ）に示されるような一般的な非線形方程式が与えられた場合、一次テイラー級数展開により近似を行うと（式９ｂ）となり、（式９ｂ）を変数ｘについて解くと（式９ｃ）となる。変数ｘの具体的な計算は（式９ｄ）に示されるように初期値ｘrefを与えて、繰り返し計算により、所定範囲に収束するまで行うことで求められる。
【００５３】
【数９】
0005567805.tif 000033
・・・・・・（式９ａ）
0005567805.tif 000034
・・・・・・（式９ｂ）
0005567805.tif 000035
・・・・・・（式９ｃ）
0005567805.tif 000036
・・・・・・（式９ｄ）
ここで、
0005567805.tif 000037
【００５４】
［飛翔体の位置算出アルゴリズム］
したがって、(式３）で表わされる非線形な観測方程式の解法を本発明にかかる飛翔体探知方法ないしシステムに適用した場合には、図１６に示される処理フローとなる。以下、図１６に基づいて説明する。
【００５５】
まず、Ｓ１６０１において処理をスタートする。具体的には、衛星１０１ａ〜１０１ｃのうちの２機の衛星から得られる撮像画像から、観測量算出処理により飛翔体頂角及び方位角を求めるものである。すわなち、人工衛星より送信された撮像画像データは、アンテナ１０３及び受信部２０１により受信され、位置解析装置１０５内のＣＰＵ１０５２の制御により通信インタフェース１０５１を介してメモリ１０５３ないしディスク１０５５に記録され、撮像画像が解析処理されて飛翔体頂角及び方位角が求められる（Ｓ１６０２）。このとき、処理に必要なプログラムが、例えばディスク１０５５から適宜メモリ１０５３に読み込まれ実行される。併せて、各衛星の姿勢情報、速度情報、３次元位置情報が各衛星からそれぞれ送信されて、適宜アンテナ１０３及び受信部２０１により受信され位置解析装置１０５内に取り込まれる（Ｓ１６０５、Ｓ１６０６、Ｓ１６０７）。
【００５６】
次に、Ｓ１６０３において、飛翔体頂角及び方位角から、次式で示される観測量算出処理により、観測ベクトルyを算出生成する（Ｓ１６０４）。
【００５７】
【数１０】
0005567805.tif 000038
・・・・・・（式１０）
【００５８】
次に、Ｓ１６０８において、各衛星の位置、速度、姿勢角から式４，５，７に示した座標変換行列を算出する。
【００５９】
次に、飛翔体位置の初期値（上述の変数ｘに関するｘref参照）を決定し（Ｓ１６０９）、上記座標変換行列（式４，５，７）と衛星の位置と飛翔体の初期値とから次式で示される各係数を算出する（Ｓ１６１０）。
【００６０】
【数１１】
0005567805.tif 000039
・・・・・・（式１１ａ）
0005567805.tif 000040
・・・・・・（式１１ｂ）
0005567805.tif 000041
・・・・・・（式１１ｃ）
【００６１】
次に、求めた係数を用いて次式を繰り返し計算することにより最終的に求める飛翔体の位置ベクトルが算出される。
【００６２】
【数１２】
0005567805.tif 000042
・・・・・・（式１２）
【００６３】
すなわち、Ｓ１６１１において（式１２）に示されるパラメータ推定処理が実行されると、都度飛翔体位置が算出され（Ｓ１６１２）、この位置が予め定めた所定の範囲内かどうかが判断される（Ｓ１６１３）。所定の範囲内でなければ（Ｓ１６１３において、Ｎｏ）、Ｓ１６１０に戻り係数生成処理を行う。Ｓ１６１３において所定の範囲に入っていると判断されれば（Ｙｅｓ）、Ｓ１６１４に進み処理を終了する。
【００６４】
なお、（式１１ｃ）に示した係数は関数ｆ（・）のｘによる偏微分であり、次式により導出される。
【００６５】
【数１３】
0005567805.tif 000043
・・・・・・（式１３）
ここで、Ｌ及びＩは、それぞれ
【００６６】
【数１４】
0005567805.tif 000044
・・・・・・（式１４ａ）
0005567805.tif 000045
・・・・・・（式１４ｂ）
で表される。
【００６７】
［飛翔体の飛翔方向（速度）算出］
また、飛翔体の飛翔方向（速度）は、Δｔを人工衛星に搭載された観測器における観測時間間隔としたとき、（式１２）から求められる時刻tにおける飛翔体の位置ベクトル
0005567805.tif 000046
と時刻ｔ-Δtにおける飛翔体の位置ベクトル
0005567805.tif 000047
とを用いて次式より算出される。
【００６８】
【数１５】
0005567805.tif 000048
・・・・・・（式１５）
【００６９】
以上のように処理された結果、飛翔体が撮影された衛星画像（衛星写真）から飛翔体の位置及び速度方向（ベクトル）を特定することができる。
【実施例２】
【００７０】
実施例１においては、２機の人工衛星の観測時刻が「同時」であることを前提としたが、例えばカルマンフィルタのような動的な状態遷移を組み入れた推定アルゴリズムを用いると、次に示すように観測時刻は必ずしも「同時」である必要はない。以下に詳述する。
【００７１】
一般的にカルマンフィルタは、観測更新及び時間伝播の２つの処理構成をとり、それぞれ次式にて表されるものである。
＜観測更新処理＞
【００７２】
【数１６】
0005567805.tif 000049
・・・・・・（式１６）
【００７３】
ここで、ｆ（・）は観測関数、ｘは状態変数、ｙは観測量、Ｐは状態変数の共分散行列、Ｒは観測量の分散、Ｈは観測関数の状態変数による偏微分である。また、下付き添え字ｉはそれがｉ番目の時刻におけるものであることを示し、上付き添え字は、−が更新前、＋が更新後をそれぞれ意味する。
【００７４】
＜時間伝播処理＞
【数１７】
0005567805.tif 000050
・・・・・・（式１７）
【００７５】
ここで、ｇ（・）は時間伝播関数、Φiは状態変数ｘの状態遷移行列、Ｑは状態変数ｘのプロセスノイズである。
【００７６】
以上説明したカルマンフィルタを本発明にかかる飛翔体探知方法ないしシステムに適用した場合、観測関数ｆ（・）、状態変数ｘ、観測量ｙ等の各要素は次のとおりとすればよい。
【００７７】
【数１８】
0005567805.tif 000051
・・・・・・（式１８ａ）
0005567805.tif 000052
・・・・・・（式１８ｂ）
0005567805.tif 000053
・・・・・・（式１８ｃ）
0005567805.tif 000054
・・・・・・（式１８ｄ）
0005567805.tif 000055
・・・・・・（式１８e）
0005567805.tif 000056
・・・・・・（式１８ｆ）
0005567805.tif 000057
・・・・・・（式１８ｇ）
0005567805.tif 000058
・・・・・・（式１８ｈ）
【実施例３】
【００７８】
実施例１及び実施例２においては、人工衛星の姿勢角が既知であることを前提にし、衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）までの座標変換を軌道座標系（Ｏ）を介して算出したが、本発明はこれに限定されない。すわなち、姿勢角が既知でない場合であっても本発明の実施は可能であり、例えば、代替的な座標系（後述のＮＥＤ座標系など）を介して演算処理すればよい。
【００７９】
加えて、人工衛星の姿勢角が既知の場合であっても未知の場合であっても人工衛星の姿勢は精密に制御されている必要はなく、かかる人工衛星から撮影された画像に飛翔体が含まれていれば足りる。例えば、撮影された画像における地球は、必ずしもその中心が画面に収まっている必要はなく、例え地球の中心が撮影画像から大きくはみ出ていたとしても、画像処理によってその中心位置を推定演算可能である。このことは、人工衛星のロール角とピッチ角とが画像によって推定可能であることを意味している。
【００８０】
また、例えば撮影画像は、地球の画像の一部が撮影されていれば回転していても良い。この場合、観測画像データと地図における地形データとのパターンマッチング処理技術によって、地球の北方向及び東方向を推定可能である。このことは、人工衛星のヨー角に相当する角度が画像によって推定可能であることを意味している。より詳細には、画像処理により地球中心及び地球の北方向、東方向を特定し、各衛星のロール角、ピッチ角、ヨー角に相当する
0005567805.tif 000059
を推定できれば、衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）までの座標変換を、例えば図１０に示すＮＥＤ（Ｎｏｒｔｈ−ｅａｓｔ−ｄｏｗｎ）座標系を介して演算すればよい。
【００８１】
［ＮＥＤ座標系の定義及び地球中心固定座標系との関係］
図１０において、地球中心Ｏを原点とするｘｙｚ座標系（ＥＣＥＦ座標系）において人工衛星Ｓが図に示した位置（ｘ軸からｚ軸の周りに回転させた方位角λs，ｘｙ平面からの仰角θs）にあるとすると、ＮＥＤ座標系は衛星Ｓを原点とした３つの座標軸Ｎ，Ｅ，Ｄで表される。
【００８２】
［撮像画像とＮＥＤ座標系との関係］
次に、図１１を用いて撮像画面におけるＮＥＤ座標系と衛星ロール角、ピッチ角、ヨー角との関係を説明する。図中では、撮像画面の中心Ｏ1と地球の中心Ｏ2とは若干のずれを伴っている。そして、衛星機体座標系におけるピッチ軸
0005567805.tif 000060
がＯ2を原点とする横座標（Ｘ軸）と一致し、衛星機体座標系におけるロール軸
0005567805.tif 000061
がＯ2を原点とする縦座標（Ｙ軸）と一致するものとする。
【００８３】
また、ＮＥＤ座標系における東方向の軸をＥ、同座標系における北方向の軸をＮとすると、
衛星のピッチ角は
0005567805.tif 000062
衛星のロール角は
0005567805.tif 000063
衛星のヨー角は
0005567805.tif 000064
となる。図中右下に示した時計回りの矢印は、上記角度を定義するときの正の方向を示す。
【００８４】
なお、図１１においては全地球が画像に納まったものとして説明しているが、本発明はこれに限定されないことは上述したとおりである。地球の画像は一部分が納まっていればよく、地球の部分的な画像であっても地球の中心、北方向、東方向は推定可能である。
【００８５】
＜２衛星、衛星姿勢情報無の場合の座標系変換＞
衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）までの座標変換行列ＣS2ECEFは、図１４に示されるようにＮＥＤ座標系を介して行われる。
【００８６】
衛星機体座標系（Ｓ）からＮＥＤ座標系への座標変換行列
0005567805.tif 000065
と、各衛星の既知の３次元位置ベクトル
【数１９】
0005567805.tif 000066
・・・・・・（式１９）
に基づく、ＮＥＤ座標系から地球中心固定座標系（ＥＣＥＦ）までの座標変換行列ＣNED2ECEF（ＲS）により、衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）までの座標変換行列ＣS2ECEFは次のとおり決定される。
【００８７】
【数２０】
0005567805.tif 000067
・・・・・・（式２０ａ）
0005567805.tif 000068
・・・・・・（式２０ｂ）
ここで、
【００８８】
【数２１】
0005567805.tif 000069
・・・・・・（式２１ａ）
0005567805.tif 000070
・・・・・・（式２１ｂ）
0005567805.tif 000071
・・・・・・（式２１ｃ）
0005567805.tif 000072
・・・・・・（式２１ｄ）
である。なお、
【００８９】
【数２２】
0005567805.tif 000073
・・・・・・（式２２ａ）
0005567805.tif 000074
・・・・・・（式２２ｂ）
となる。
【００９０】
＜２衛星、衛星姿勢情報無の場合の演算処理フロー＞
本実施例における具体的演算処理例について説明する。
【００９１】
［飛翔体の位置算出アルゴリズム］
次に、(式３）で表わされる非線形な観測方程式を現実的に演算処理することを踏まえた処理フローは、図１７に示されるとおりである。以下、図１７に基づいて説明する。
【００９２】
まず、Ｓ１７０１において処理をスタートする。具体的には、衛星１０１ａ〜１０１ｃのうちの２機の衛星から得られる撮像画像から、観測量算出処理により飛翔体頂角及び方位角を求めるものである。すわなち、送信された画像は、アンテナ１０３及び受信部２０１により受信されて位置解析装置１０５内のＣＰＵ１０５２の制御により通信インタフェース１０５１を介してメモリ１０５３ないしディスク１０５５に記録され、撮像画像が解析処理されて飛翔体頂角及び方位角が求められる（Ｓ１７０２）。このとき、処理に必要なプログラムが、例えばディスク１０５５から適宜メモリ１０５３に読み込まれ実行される。
【００９３】
次に、Ｓ１７０３において、飛翔体頂角及び方位角から、次式で示される観測量算出処理により、観測ベクトルｙが算出生成される（Ｓ１７０４）。
【００９４】
【数２３】
0005567805.tif 000075
・・・・・・（式２３）
【００９５】
次に、Ｓ１７０５において、撮像画像における地球表面の地形形状と、他の地図データにおける地形形状とをパターンマッチング処理技術により比較し、人工衛星自身の姿勢差（ロール角、ピッチ角、ヨー角）を算出する。これによって、２衛星のロール角、ピッチ角、ヨー角
0005567805.tif 000076
が求められてメモリ１０５３に格納される（Ｓ１７０６）。
【００９６】
併せて、各衛星の速度情報及び３次元位置情報が各衛星からそれぞれ送信されて、適宜アンテナ１０３及び受信部２０１で受信して位置解析装置１０５内に取り込まれる（Ｓ１７０７、Ｓ１７０８）。
【００９７】
次に、Ｓ１７０９において、各衛星の位置、速度、姿勢角から（式２０a，２０ｂ，２１ｄ）に示した座標変換行列を算出する。
【００９８】
次に、飛翔体位置の初期値（上述の変数ｘに関するｘref参照）を決定し（Ｓ１７１０）、上記座標変換行列（式２０a，２０ｂ，２１ｄ）と衛星の位置と飛翔体の初期値とから次式で示される各係数を算出する（Ｓ１７１１）。
【００９９】
【数２４】
0005567805.tif 000077
・・・・・・（式２４a）
0005567805.tif 000078
・・・・・・（式２４ｂ）
0005567805.tif 000079
・・・・・・（式２４ｃ）
【０１００】
次に、求めた係数を用いて次式を繰り返し計算することにより最終的に求める飛翔体の位置ベクトルが算出される。
【０１０１】
【数２５】
0005567805.tif 000080
・・・・・・（式２５）
【０１０２】
すなわち、Ｓ１７１２において（式２５）に示されるパラメータ推定処理が実行されると、都度飛翔体位置が算出され（Ｓ１７１３）、この位置が予め定めた所定の範囲内かどうかが判断される（Ｓ１７１４）。所定の範囲内でなければ（Ｓ１７１４において、Ｎｏ）、Ｓ１７１１に戻り係数生成処理を行う。Ｓ１７１４において所定の範囲に入っていると判断されれば（Ｙｅｓ）、Ｓ１７１５に進み処理を終了する。
【０１０３】
［飛翔体の飛翔方向（速度）算出］
また、飛翔体の飛翔方向（速度）は、Δｔを人工衛星に搭載された観測器における観測時間間隔としたとき、（式２５）から求められる時刻ｔにおける飛翔体の位置ベクトル
0005567805.tif 000081
と時刻ｔ-Δtにおける飛翔体の位置ベクトル
0005567805.tif 000082
とを用いて次式より算出される。
【０１０４】
【数２６】
0005567805.tif 000083
・・・・・・（式３１）
【０１０５】
以上のように処理された結果、飛翔体が撮影された衛星画像（衛星写真）から飛翔体の位置及び速度方向（ベクトル）を特定することができる。
【実施例４】
【０１０６】
実施例１〜実施例３においては、２機の人工衛星によって観測された２枚の画像データに対して上述したとおりのステレオ視（立体視）アルゴリズムを用いて飛翔体の３次元位置及び飛翔方向（速度）を算出する例について説明した。実施例４では、３機の人工衛星によって観測された３枚の画像データを用いて飛翔体の３次元位置及び飛翔方向（速度）を算出する方法等について説明する。
【０１０７】
３機の人工衛星を使用した場合の特徴の１つは、人工衛星のヨー軸周りのヨー角姿勢
0005567805.tif 000084
も同時に推定できることである。そのため、人工衛星のロール角、ピッチ角、ヨー角姿勢が未知な場合、あるいは画像パターンマッチングを正確に行えずに北方向ないし東方向を正確に把握できない場合であっても、地球の中心さえ撮影画像から推定できれば対象とする飛翔体の３次元位置を特定することができる。図１２は、これら３機の人工衛星に搭載されている観測機器から観測された画像データを示している。
【０１０８】
この場合、直接的に得られる観測量は、撮影画像中心からの角度（頂角）である
0005567805.tif 000085
と、画像の縦横を基準にしたときの画像回転角（方位角）である
0005567805.tif 000086
とからなる６つのパラメータである。
【０１０９】
このとき、対象となる飛翔体の３次元位置
【数２７】
0005567805.tif 000087
・・・・・・（式３２）
は、２機の人工衛星の場合でも説明したように、衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）への座標変換行列ＣS2ECEFと前述（式１）及び（式２）との関係から、各衛星の撮像画面での観測量
0005567805.tif 000088
との間の下記の関係を導くことができる。
【０１１０】
【数２８】
0005567805.tif 000089
・・・・・・（式３３）
【０１１１】
そして、３機の人工衛星を使用し、人工衛星のヨー軸周りのヨー角姿勢も推定する場合は、衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）への座標変換は地心方向軌道座標系（ＥＤＯ）を経由して行う。
【０１１２】
［衛星の地心方向軌道座標系（ＥＤＯ）の定義］
まず、図６（ｂ）に基づいて、衛星軌道座標系（Ｏ）と衛星の地心方向軌道座標系（ＥＤＯ）との関係から地心方向軌道座標系（ＥＤＯ）について説明する。一般に、衛星軌道の離心率はゼロでは無い（すなわち衛星軌道は真円では無い）ので、衛星軌道座標系（Ｏ）と衛星の地心方向軌道座標系（ＥＤＯ）とは異なるものとなる。
【０１１３】
いま、図６（ｂ）において示す位置に楕円軌道上を周回する衛星Ｓがあるとすると、軌道座標系（Ｏ）における衛星Ｓは、軌道面内で速度ベクトル方向の単位ベクトル
0005567805.tif 000090
と、上記速度ベクトル方向とそれに直交する面外方向ベクトルに直交する単位ベクトル
0005567805.tif 000091
とで表される。一方で、地心方向軌道座標系（ＥＤＯ）における衛星Ｓは、地心方向の単位ベクトル
0005567805.tif 000092
と、地心方向に直交し、軌道面内にある単位ベクトル
0005567805.tif 000093
とで表される。
【０１１４】
［撮像画像と地心方向軌道座標系（ＥＤＯ）との関係］
次に、図７に撮像画像と地心方向軌道座標系（ＥＤＯ）との関係を示す。いま、図中の矩形外枠が撮像画面の外枠であったとすると、その撮像画面の中心Ｏ1と地球の中心Ｏ2とは若干のずれを伴っている。そして、衛星機体座標系におけるピッチ軸
0005567805.tif 000094
がＯ2を原点とする横座標（Ｘ軸）と一致し、衛星機体座標系におけるロール軸
0005567805.tif 000095
がＯ2を原点とする縦座標（Ｙ軸）と一致するものとし、地心方向に直交し、軌道面内にあるベクトル方向のＲＳ−ＰＳ平面への射影された単位ベクトルを
0005567805.tif 000096
とし、ｘＥＤＯ，ｚＥＤＯと右手系を成す単位ベクトル方向のＲＳ−ＰＳ平面への射影された単位ベクトルを
0005567805.tif 000097
とすれば、
衛星のピッチ角は
0005567805.tif 000098
衛星のロール角は
0005567805.tif 000099
衛星のヨー角は
0005567805.tif 000100
となる。なお、図中右下に示した時計回りの矢印は、上記角度を定義するときの正の方向を示す。
【０１１５】
＜３衛星の場合の座標系変換＞
衛星機体座標系（Ｓ）から地球中心固定座標系（ＥＣＥＦ）への座標変換行列ＣS2ECEFは、図１５に示されるように地心方向軌道座標系を介して行われる。
【０１１６】
衛星機体座標系（Ｓ）から地心方向軌道座標系（ＥＤＯ）への座標変換行列ＣS2EDOと、地心方向軌道座標系（ＥＤＯ）から地球中心固定座標系（ＥＣＥＦ）への座標変換行列ＣEDO2ECEFとから次式で表される。
【０１１７】
【数２９】
0005567805.tif 000101
・・・・・・（式３４a）
0005567805.tif 000102
・・・・・・（式３４ｂ）
ここで、
0005567805.tif 000103
・・・・・・（式３４ｃ）
0005567805.tif 000104
・・・・・・（式３４ｄ）
0005567805.tif 000105
・・・・・・（式３４ｅ）
である。
【０１１８】
また、衛星機体座標系（Ｓ）から地心方向軌道座標系（ＥＤＯ）への座標変換行列ＣS2EDOは、画像から地球中心を推定することで得られる各衛星のロール角、ピッチ角
0005567805.tif 000106
と、推定パラメータとして扱うヨー角
0005567805.tif 000107
とに基づき、各軸周りの回転行列を用いて次式で表される。
【０１１９】
【数３０】
0005567805.tif 000108
・・・・・・（式３５）
ここで、
【０１２０】
【数３１】
0005567805.tif 000109
・・・・・・（式３６a）
0005567805.tif 000110
・・・・・・（式３６ｂ）
0005567805.tif 000111
・・・・・・（式３６ｃ）
である。
【０１２１】
[観測方程式の導出]
以上より、３衛星のステレオ視（立体視）により、各衛星のヨー軸も推定する場合には、次式で示される観測方程式を解くことで
【数３２】
0005567805.tif 000112
・・・（式３７）
飛翔体の３次元位置
【数３３】
0005567805.tif 000113
・・・（式３８）
及び各衛星のヨー軸
0005567805.tif 000114
を算出することができる。
【０１２２】
＜３衛星の場合の演算処理フロー＞
本実施例における具体的演算処理例について説明する。
【０１２３】
［飛翔体の位置算出アルゴリズム］
次に、(式３７）で表わされる非線形な観測方程式を現実的に演算処理することを踏まえた処理フローは、図１８に示されるとおりである。以下、図１８に基づいて説明する。
【０１２４】
まず、Ｓ１８０１において処理をスタートするが、具体的には、衛星１０１ａ〜１０１ｃの３機の衛星から得られる撮像画像から、観測量算出処理により飛翔体頂角及び方位角が求められる。すわなち、送信された画像は、アンテナ１０３及び受信部２０１により受信され位置解析装置１０５内のＣＰＵ１０５２の制御により通信インタフェース１０５１を介してメモリ１０５３ないしディスク１０５５に記録され、撮像画像が解析処理されて飛翔体頂角及び方位角が求められる（Ｓ１８０２）。このとき、処理に必要なプログラムが、例えばディスク１０５５から適宜メモリ１０５３に読み込まれ実行される。
【０１２５】
次に、Ｓ１８０３において、飛翔体頂角及び方位角から、次式で示される観測量算出処理により、観測ベクトルｙが算出生成される（Ｓ１８０４）。
【数３４】
0005567805.tif 000115
・・・・・・（式３９）
【０１２６】
次に、Ｓ１８０５において、撮像画像における地球表面の地形形状と、他の地図データにおける地形形状とをパターンマッチング処理技術により比較し、人工衛星自身の姿勢差（ロール角、ピッチ角）を算出する。これによって、３衛星のロール角、ピッチ角
0005567805.tif 000116
が求められてメモリ１０５３に格納される（Ｓ１８０６）。
【０１２７】
併せて、各衛星の速度情報、３次元位置情報が各衛星からそれぞれ送信されて、適宜アンテナ１０３及び受信部２０１で受信して位置解析装置１０５内に取り込まれる（Ｓ１８０７、Ｓ１８０８）。
【０１２８】
次に、Ｓ１８０９において、各衛星の位置、速度、姿勢角から（式３４ｂ，３６ｂ，３６ｃ）に示した座標変換行列を算出する。
【０１２９】
次に、飛翔体位置の初期値（上述の変数ｘに関するｘref参照）を決定し（Ｓ１８１０）、上記座標変換行列（式３４ｂ，３６ｂ，３６ｃ）と衛星の位置と飛翔体の初期値とから次式で示される各係数を算出する（Ｓ１８１１）。
【数３５】
0005567805.tif 000117
・・・・・・（式４０ａ）
0005567805.tif 000118
・・（式４０ｂ）
0005567805.tif 000119
・・・・・・（式４０ｃ）
0005567805.tif 000120
・・・（式４０ｄ）
0005567805.tif 000121
・・（式４０ｅ）
【０１３０】
以上求めた係数を用いて次式を繰り返し計算することにより最終的に求める飛翔体の位置ベクトルが算出される。
【０１３１】
【数３６】
0005567805.tif 000122
・・・・・・（式４１）
【０１３２】
すなわち、Ｓ１８１２において（式４１）に示されるパラメータ推定処理が実行されると、都度飛翔体位置及び衛星姿勢が算出され（Ｓ１８１３）、この位置が予め定めた所定の範囲内かどうかが判断される（Ｓ１８１４）。所定の範囲内でなければ（Ｓ１８１４において、Ｎｏ）、Ｓ１８１１に戻り係数生成処理を行う。Ｓ１８１４において所定の範囲に入っていると判断されれば（Ｙｅｓ）、Ｓ１８１５に進み処理を終了する。
【０１３３】
なお、（式４０ｃ、４０ｄ、４０e）に示した係数は関数ｆ（・）の推定パラメータである飛翔体の３次元位置
【数３７】
0005567805.tif 000123
・・・・・・（式４２）
による偏微分、及び各衛星のヨー軸
0005567805.tif 000124
による偏微分であり、次式により導出される。
【０１３４】
【数３８】
0005567805.tif 000125
・・・・・・（式４３）
ここで、
【０１３５】
【数３９】
0005567805.tif 000126
0005567805.tif 000127
・・・（式４４）
であり、
【０１３６】
【数４０】
0005567805.tif 000128
・・・（式４５a）
0005567805.tif 000129
・・・（式４５ｂ）
0005567805.tif 000130
・・・（式４５ｃ）
である。また、
【０１３７】
【数４１】
0005567805.tif 000131
・・・（式４６）
である。
【０１３８】
［飛翔体の飛翔方向（速度）算出］
また、飛翔体の飛翔方向（速度）は、Δｔを人工衛星に搭載された観測器における観測時間間隔としたとき、（式４１）から求められる時刻tにおける飛翔体の位置ベクトル
0005567805.tif 000132
と時刻ｔ-Δtにおける飛翔体の位置ベクトル
0005567805.tif 000133
とを用いて次式より算出される。
【０１３９】
【数４２】
0005567805.tif 000134
・・・・・・（式４７）
【０１４０】
以上のとおり、２機の人工衛星より撮影された画像から飛翔体の位置及び速度ベクトルを特定可能であるとともに（実施例１〜実施例３）、３機の人工衛星より撮影された画像を使用すれば、さらにヨー角姿勢を推定することが可能となる（実施例４）。
【０１４１】
また、本発明は必ずしも図１に示したシステム全体として構成してはじめて完成し得るものではなく、２機ないし３機の人工衛星が撮影した画像（及び必要であれば撮影した人工衛星の位置情報等）データのみを使用し、対象とする飛翔体の３位次元位置及び速度ベクトルの特定を行うプログラムないしかかるプログラムが実行されるコンピュータという構成をとることもできる。
【０１４２】
なお、念のために申し添えておくが、特許請求の範囲、明細書、要約書、図面に記載された全ての技術的要素、ならびに、方法ないし処理ステップは、これら要素及び／又はステップの少なくとも一部が相互に排他的となる組み合わせを除く任意の組み合わせによって本発明にかかる方法及びシステムならびにプログラムの構成要素となりうる。
【０１４３】
なお、本発明は、上述した実施形態のいずれの個別具体的な詳細記載に制限されることはない。さらに、本発明の技術的範囲は、上述の説明のみによってではなく、特許請求の範囲の記載によってその外延が特定されるものであり、特許請求の範囲と均等となる置換ないし変更も本発明の技術的範囲となるものである。
【符号の説明】
【０１４４】
１００  飛翔体探知システム
１０１ａ〜ｃ  人工衛星
１０２  基地局
１０３  アンテナ
１０４  信号線
１０５  位置解析装置
２０１  受信部
１０５１  通信インタフェース
１０５２  ＣＰＵ
１０５３  メモリ
１０５４  ＤＳＰ
１０５５  ディスク
１０５６  表示部
１０５７  入力部
【図１】
0005567805.tif 000167
【図２】
0005567805.tif 000168
【図４】
0005567805.tif 000169
【図５】
0005567805.tif 000170
【図６】
0005567805.tif 000171
【図７】
0005567805.tif 000172
【図８】
0005567805.tif 000173
【図９】
0005567805.tif 000174
【図１０】
0005567805.tif 000175
【図１１】
0005567805.tif 000176
【図１２】
0005567805.tif 000177
【図１３】
0005567805.tif 000178
【図１４】
0005567805.tif 000179
【図１５】
0005567805.tif 000180
【図１６】
0005567805.tif 000181
【図１７】
0005567805.tif 000182
【図１８】
0005567805.tif 000183
【図３】
0005567805.tif 000184
【図１９】
0005567805.tif 000185
【図２０】
0005567805.tif 000186