(57)【要約】【課題】既設置の発信局から移動体に到達する電波を用いることにより新規に送受信機を必要としない軌道推定方法とその装置を提供する。【解決手段】軌道推定方法は、複数の送信局から複数の移動体に電波を同時に送信し、２局の受信局において、上記移動体に搭載された中継器により中継送信された電波を受信し、上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出し、上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることのできる局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差から直接上記移動体の軌道を推定する。【選択図】図３
2008064566.tif 000002
【特許請求の範囲】
【請求項１】
複数の送信局から複数の移動体に電波を同時に送信し、
２局の受信局において、上記移動体に搭載された中継器により中継送信された電波を受信し、
上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出し、
上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることのできる局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差から直接上記移動体の軌道を推定することを特徴とする軌道推定方法。
【請求項２】
少なくとも１つの送信局から複数の移動体に電波を時系列に送信し、
２局の受信局において、上記移動体に搭載された中継器により中継送信された電波を受信し、
上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出して、
上記送信局の同時の送信では上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることができない局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差からケプラリアン軌道６要素を推定することを特徴とする軌道推定方法。
【請求項３】
複数の移動体に電波を同時に送信する複数の送信局と、
上記移動体に搭載された中継器により中継送信された電波を受信する複数の受信局と、
上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出するＴＤＯＡ／ＦＤＯＡ計算装置と、
上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることのできる局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差から直接上記移動体の軌道を推定する移動体軌道計算装置と、
を備えることを特徴とする移動体軌道推定装置。
【請求項４】
複数の移動体に電波を時系列に送信する少なくとも１つの送信局と、
上記移動体に搭載された中継器により中継送信された電波を受信する複数の受信局と、
上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出するＴＤＯＡ／ＦＤＯＡ計算装置と、
上記送信局の同時の送信では上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることができない局数の送信局が存在するときにも、得られた到来時間の差およびドップラー周波数の差からケプラリアン軌道６要素を推定し、上記推定されたケプラリアン軌道６要素から任意の時刻の上記移動体の３次元の位置および３次元の速度を推定することを特徴とする移動体軌道推定装置。
【請求項５】
上記移動体の機数は２であり、上記送信局の局数は６以上であり、
上記移動体軌道計算装置は、同時に送信された電波により得られた到来時間の差およびドップラー周波数の差より、上記移動体の３次元の位置および３次元の速度を推定することを特徴とする請求項３に記載する移動体軌道推定装置。
【請求項６】
上記移動体の機数は２であり、上記送信局の局数は１局以上であり、
上記移動体軌道計算装置は、時系列に送信された電波により得られた到来時間の差およびドップラー周波数の差より、上記移動体のケプラリアン軌道６要素を推定し、上記推定されたケプラリアン軌道６要素から任意の時刻の上記移動体の３次元の位置および３次元の速度を推定することを特徴とする請求項４に記載する移動体軌道推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、移動体の軌道を推定する方法およびその装置に関するものである。
【背景技術】
【０００２】
位置が未知の不法電波送信局などの位置を評定する方法では、２機の衛星を介して得られる電波の到来時間差（Ｔｉｍｅ  Ｄｉｆｆｅｒｅｎｃｅ  ｏｆ  Ａｒｒｉｖａｌ：ＴＤＯＡ）とドップラー周波数差（Ｆｒｅｑｕｅｎｃｙ  Ｄｉｆｆｅｒｅｎｃｅ  ｏｆ  Ａｒｒｉｖａｌ：ＦＤＯＡ）情報を用いて、位置を評定する方法が提案されている（例えば、非特許文献２参照）。この方法では、位置が未知の送信局から衛星に向けて送信される送信電波は、メインローブの他にサイドローブなど、空間的に幅をもっているので、２機の衛星により受信することができる。そして、衛星は中継器を搭載しており、電波を地上の受信局へ送り返す中継送信機能を持っているため、受信局では衛星を介して送り返された電波を受信することができる。２局の受信局の受信信号処理を行うことにより、信号の到来時間差およびドップラー周波数差を得ることができる。
【０００３】
この位置が未知な送信局の位置を評定するためには、衛星の位置および速度からなる軌道情報が必要となる。そこで、衛星の軌道を推定する従来の方法として、衛星の距離や方位角を観測し、その情報を用いて軌道推定を行っている。この方式では、光学望遠鏡により衛星からの光を観測したり、衛星との間で電波の送受信することにより、観測局から見た衛星の距離、仰角および方位角を測定したりしている（例えば、非特許文献１参照）。
【０００４】
【非特許文献１】有本好徳、他２名、「静止衛星測位システムにおける軌道決定精度について」、宇宙航行エレクトロニクス研究会資料、社団法人電子情報通信学会、１９９７年、ＳＡＮＥ８９−２６、ｐ．４５−５２
【非特許文献２】Ｄ．Ｐ．Ｈａｗｏｒｔｈ、他３名、「Ｉｎｔｅｒｆｅｒｅｎｃｅ  Ｌｏｃａｌｉｚａｔｉｏｎ  Ｆｏｒ  ＥＵＴＥＬＳＡＴ  Ｓａｔｅｌｌｉｔｅｓ−Ｔｈｅ  Ｆｉｒｓｔ  Ｅｕｒｏｐｅａｎ  Ｔｒａｎｓｍｉｔｔｅｒ  Ｌｏｃａｔｉｏｎ  Ｓｙｓｔｅｍ」、Ｉｎｔｅｒｎａｔｉｎａｌ  Ｊｏｕｒｎａｌ  ｏｆ  Ｓａｔｅｌｌｉｔｅ  Ｃｏｍｍｕｎｉｃａｔｉｏｎｓ、Ｖｏｌ．１５、ｐ．１５５−１８３
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかし、衛星からの光を観測するためには、大型の光学望遠鏡を設置しなければならないし、衛星との間で電波の送受信を行う送受信機を別途必要になったりするので、装置全体の規模が大きくなるという問題がある。
【０００６】
この発明は、既設置の発信局から移動体に到達する電波を用いることにより新規に送受信機を必要としない軌道推定方法とその装置を提供することである。
【課題を解決するための手段】
【０００７】
この発明に係わる軌道推定方法は、複数の送信局から複数の移動体に電波を同時に送信し、２局の受信局において、上記移動体に搭載された中継器により中継送信された電波を受信し、上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出し、上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることのできる局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差から直接上記移動体の軌道を推定するものである。
【発明の効果】
【０００８】
この発明に係わる移動体軌道推定によれば、位置、仰角、方位角情報ではなく、到来時間差とドップラー周波数差情報を用いるので、放送局などの正規な送信局からの電波を用いることができ、軌道推定のために新たに電波を送信する必要がなく、その結果、軌道推定用に送信局を設置する必要がない。
【発明を実施するための最良の形態】
【０００９】
実施の形態１．
図１は、この発明の実施の形態１に係わる移動体軌道推定装置のブロック図である。図２は、実施の形態１に係わる移動体軌道推定装置により２機の衛星の軌道を推定する様子を示す概念図である。
この発明の実施の形態１に係わる移動体軌道推定装置１は、図１、２に示すように、２機の軌道が未知の衛星２の軌道を推定するために、６個の到来時間の差τｍ（ｍは１〜６の整数）と６個のドップラー周波数の差ｆｍ（ｍは１〜６の整数）を用いているので、位置ｐｒ，ｍが既知の６つの送信局３、１つの基準送信局７、各送信局３、７から各衛星２に到来し、各衛星２で中継された電波を受信する２つの受信局４、到来時間の差τｍとドップラー周波数の差ｆｍを算出するＴＤＯＡ／ＦＤＯＡ計算装置５、到来時間の差τｍとドップラー周波数の差ｆｍを用いて各衛星２の位置ｐｓｉ（ｉは衛星を示し、１または２である）および衛星２の速度ｖｓｉを推定する衛星軌道計算装置６を備える。この送信局３、７にはそれぞれ１台の発信機が設置され、受信局４にはそれぞれ１台の受信機が設置されている。なお、到来時間の差τｍとドップラー周波数の差ｆｍの算出は、一般的な方法であるので、説明は省略する。
【００１０】
図３は、実施の形態１に係わる衛星軌道計算装置６の機能ブロック図である。
実施の形態１に係わる衛星軌道計算装置６は、図３に示すように、収束条件、送信局３、７の位置、電波波長λ、電波伝搬速度ｃなど各種パラメータを入力する設定パラメータ入力部１１、収束演算のための軌道の初期値を設定する軌道初期値設定部１２、到来時間の差τｍとドップラー周波数の差ｆｍを用いて収束演算を行う収束演算部１３、収束演算により得た軌道が収束条件を満たしているか否かを判定する収束条件判定部１４を有する。
【００１１】
設定パラメータ入力部１１は、収束条件のための設定パラメータεｐ、εｖ、収束演算のための６つの送信局３の位置ｐｒ，ｍ（ｍは１〜６）、基準送信局７の位置ｐｒ０、電波波長λ、電波伝搬速度ｃを入力する。
軌道初期値設定部１２は、２機の衛星２の初期軌道χ（ハット）を設定する。
収束演算部１３は、式（１）、式（２）に関する収束演算を行い、軌道χを算出する。このとき、Ｇ（χ（ハット）、ｍ）は、式（３）にχ（ハット）を代入して算出する。
【００１２】
【数１】
2008064566.tif 000003
【００１３】
また、式（３）の各要素は式（４）から式（９）により算出する。
【００１４】
【数２】
2008064566.tif 000004
【００１５】
収束条件判定部１４は、収束演算部１３で算出した軌道χに関するステップ値Δχ（ハット）が式（１０）、式（１１）の条件式を満足するか否かを判定し、満足するときその軌道を衛星２の軌道とし、満足していないときには再度収束演算部１３で収束演算を行う。
【００１６】
【数３】
2008064566.tif 000005
【００１７】
図４は、実施の形態１に係わる衛星軌道計算装置６での軌道推定手順を示すフローチャートである。
軌道推定手順を開始すると、図４に示すように、ステップＳ１０１で、収束条件、送信局３の位置、基準送信局７の位置など各種パラメータを読み込む。
ステップＳ１０２で、ＴＤＯＡ／ＦＤＯＡ計算装置５から到来時間の差とドップラー周波数の差に関する観測値を読み込む。
ステップＳ１０３で、収束演算のための軌道初期値を読み込む。
ステップＳ１０４で、収束演算を式（１）〜式（９）に基づいて行い、軌道χ（ハット）を推定する。
ステップＳ１０５で、算出した軌道χ（ハット）について式（１０）、式（１１）の収束条件式を満足しているか否かを判断し、満足しているとき軌道推定手順を終了し、満足していないときステップＳ１０４に戻る。
【００１８】
次に、式（１）、式（２）で表される軌道の推定値χ（ハット）およびステップ値Δχ（ハット）を収束演算し、式（１０）および式（１１）の収束判定条件式で判定することにより衛星軌道χを推定できる根拠を説明する。
位置が既知の６つの送信局３から２機の衛星２に同時に送信を行い、衛星２を経由して得られる信号を各受信局４で受信し、到来時間の差τｍおよびドップラー周波数の差ｆｍを計算することにより、第１送信局３から第７送信局３それぞれから衛星２への到来時間の差τｍおよびドップラー周波数の差ｆｍを観測することができる。
【００１９】
送信局３の３次元位置ｐｒ，ｍを式（１２）、２機の衛星２の３次元位置および速度を式（１３）で表す。ここで、式（１４）のｐｓｉは衛星２の３次元位置で、式（１５）のｖｓｉは衛星２の３次元速度である。
このとき、到来時間の差τｍおよびドップラー周波数の差ｆｍは、それぞれ式（１６）、式（１７）で表される。
【００２０】
【数４】
2008064566.tif 000006
【００２１】
ここで、｜｜・｜｜は、ベクトルのノルムを表し、物理的には距離を意味する。λは電波波長、ｃは電波伝播速度である。式（１６）は、送信局３と衛星２の距離の差を電波伝播速度ｃで除算した量の差、すなわち、２機の衛星２への電波の到来時間の差を意味する。
また、式（１７）は、送信局３と衛星２とを結ぶ線上への衛星速度の射影成分を電波波長λで除算した量の差、すなわち、２機の衛星２のドップラー周波数の差を意味する。
【００２２】
式（１６）および式（１７）に示すように、２機の衛星２の３次元位置および速度の未知変数は計１２個あり、１つの送信局３による送信で、到来時間の差およびドップラー周波数の差の２本の方程式を立てることができる。よって、６つ以上の送信局３を用いて同時に送信を行い、到来時間の差およびドップラー周波数の差を用いてそれぞれ６本以上、計１２本以上の方程式を立てれば、２機の衛星２の３次元位置ｐｓｉおよび速度ｖｓｉを求めることができる。
【００２３】
また、実際の到来時間の差およびドップラー周波数の差の観測値には受信局４の時計時間誤差、受信局４でのローカル周波数シフト、衛星内部でのローカル周波数シフトなどの誤差要因が存在するので、これらをキャンセルするために既知位置の基準送信局７（位置：ｐｒ０）を用意し、基準送信局７による送信も同時に行う。そして、基準送信局７による到来時間の差およびドップラー周波数の差を基準とした到来時間の差およびドップラー周波数の差を用いることにする。すなわち、式（１６）および式（１７）の代わりに式（１８）および式（１９）を用いる。
【００２４】
【数５】
2008064566.tif 000007
【００２５】
このように、６つ以上の送信局３と基準送信局７の計７局以上で２機の衛星２の３次元位置ｐｓｉおよび３次元速度ｖｓｉを推定することができる。
次に、式（１８）および式（１９）の連立方程式を衛星２の位置ｐｓｉおよび速度ｖｓｉについて解く方法を説明する。
第ｍ送信局３による到来時間の差およびドップラー周波数の差の観測値をそれぞれτｏｂｓ，ｍ、ｆｏｂｓ，ｍとおき、式（２０）のようにベクトル状に表記する。
また、θｍ（χ）は、式（１８）および式（１９）で表す到来時間の差τｍ（χ）およびドップラー周波数の差ｆｍ（χ）をベクトル状に表したものであり、式（２１）のように定義できる。
【００２６】
【数６】
2008064566.tif 000008
【００２７】
式（２２）に表される観測誤差共分散行列Ｖは、受信機雑音などによる到来時間差観測雑音の分散値στ２およびドップラー周波数差観測雑音の分散値σｆ２を対角項にもつ。
【００２８】
【数７】
2008064566.tif 000009
【００２９】
このように定義された値を用いると、推定問題を式（２３）のように定式化できる。
すなわち、式（２３）を評価関数として、衛星の軌道χについて最小化する最小二乗問題となる。
【００３０】
【数８】
2008064566.tif 000010
【００３１】
次に、式（２３）を評価関数とする最小二乗法の解を収束演算で求める方法を説明する。今、評価関数式（２３）をＦ（χ）と呼ぶ。まず、式（２４）に従って、θｍ（χ）の線形一次近似を行う。
【００３２】
【数９】
2008064566.tif 000011
【００３３】
ここで、Ｇ（χ，ｍ）は、式（２５）に示すようにθｍ（χ）のχについての偏微分であり、ベクトル表示すると式（２６）で表される。そして、Ｇ（χ，ｍ）を数式により計算する場合、Ｇ（χ，ｍ）は式（３）で表される。そして、▽ｐｓ１τｍ（χ）、▽ｐｓ２τｍ（χ）、▽ｐｓ１ｆｍ（χ）、▽ｐｓ１ｆｍ（χ）、▽ｖｓ１ｆｍ（χ）、▽ｖｓ１ｆｍ（χ）は、それぞれ式（４）、式（５）、式（６）、式（７）、式（８）、式（９）で求められる。
【００３４】
【数１０】
2008064566.tif 000012
【００３５】
すると、評価関数Ｆ（χ）は式（２７）で表される。ここで、θｔｍｐ，ｍは式（２８）で表される。
【００３６】
【数１１】
2008064566.tif 000013
【００３７】
評価関数Ｆ（χ）のχについての偏微分▽χＦ（χ）は、式（２９）で表される。ここで▽χは、χについての偏微分を意味する。
【００３８】
【数１２】
2008064566.tif 000014
【００３９】
偏微分値▽χＦ（χ）＝０１２×１としたものをχについて解くと最適解が得られ、このとき、χ（ハット）は式（１）、収束演算のステップ値Δχ（ハット）は式（２）で表される。よって、式（１）および式（２）の収束演算により最適な衛星の軌道χを推定することができる。なお、収束判定では、収束演算のステップ値Δχ（ハット）から、位置および速度に相当する要素をそれぞれ選択してベクトル上に並べた値に対し、式（１０）、式（１１）の判定条件で判定する。
【００４０】
なお、実施の形態１に係わる移動体軌道推定装置１は、２機の衛星２の軌道を推定したが、衛星２の数はこれに限るものではない。例えば、衛星２の数をＮ機とし、ｄ次元の衛星の位置と速度とを推定する場合、必要な方程式数は２×Ｎ×ｄである。次に、送信局数をＭとし、Ｍ局の送信局３から同時送信することにより、２Ｍ（Ｎ−１）の方程式が得られる。これらから、方程式の数が未知数を上回るためには、２Ｍ（Ｎ−１）が２×Ｎ×ｄ以上を満たす送信局数Ｍを用意すればよい。収束演算式の導出は、実施の形態１で説明した軌道推定方法と同様である。なお、基準送信局を用いる場合には、加えて１局必要となる。
【００４１】
このように位置が既知の既に電波を出射している放送局などの送信局３を用いて衛星２の軌道を推定することができるので、軌道推定のためだけに別途衛星２に電波を照射する必要がなく、移動体軌道推定装置１の新規な送信局の設置は不要となる。
なお、実施の形態１に係わる移動体軌道推定装置１は、衛星２に限らず任意の移動体の軌道を推定することができる。
【００４２】
図５は、送信局位置評定装置に移動体軌道推定装置１を適用したときの構成図を示す。
次に、非特許文献２に記載されている位置が未知な送信局の位置を評定する送信局位置評定装置に実施の形態１に係わる移動体軌道推定装置１を適用することを説明する。
このときは、図５に示すように、軌道推定用受信局４およびＴＤＯＡ／ＦＤＯＡ計算装置５は、送信局位置評定装置１０の位置評定用受信局およびＴＤＯＡ／ＦＤＯＡ計算装置として兼用することができ、構成が簡単な送信局位置評定装置を提供することができる。
【００４３】
なお、実施の形態１ではＧ（χ，ｍ）を数式を用いて導出しているが、近似的に求めてもよい。近似的に求める場合、まず、Ｇ（χ，ｍ）を式（３０）のように列ベクトルに分解する。このとき、各列ベクトルは式（３１）で表され、Ｕｊは、第ｊ要素だけが１であり、他の要素がすべて零であるような１２次元ベクトルである。εｊはｊに応じて十分小さい値に設定すればよい。
【００４４】
【数１３】
2008064566.tif 000015
【００４５】
実施の形態２．
図６は、この発明の実施の形態２に係わる移動体軌道推定装置のブロック図である。図７は、実施の形態２に係わる移動体軌道推定装置により２機の衛星の軌道を推定する様子を示す概念図である。
この発明の実施の形態２に係わる移動体軌道推定装置１Ｂは、実施の形態１に係わる移動体軌道推定装置１と送信局３ａ〜３ｃの局数が必ずしも６局以上ではなく、それに伴って衛星軌道計算装置６Ｂが異なっており、それ以外は同様であるので、同様な部分に同じ符号を付記して説明は省略する。
実施の形態１に係わる移動体軌道推定装置１では、２機の衛星２の１２個の未知変数の解を求めるために、６つ以上の位置が既知の送信局３から同時に送信されている電波を用いて、１２本以上の方程式を立てている。
しかし、実際に送信局３を多数設置することが難しい場合がある。そこで、実施の形態２に係わる移動体軌道推定装置１Ｂでは、観測回数を増やすことにより必要な数以上の方程式をたてる。この場合、衛星２の位置および速度は時々刻々と変化するため、衛星２の位置および速度について、直接連立方程式を解くことはできない。そこで、衛星２の位置および速度を未知変数とするのではなく、ケプラリアン軌道６要素を未知変数として推定する。すなわち、２機の衛星２のケプラリアン軌道６要素のための１２本の方程式を立てる。
【００４６】
実施の形態２においては、図７に示すように、３つの送信局３ａ〜３ｃが設置されているので、１２本の方程式を立てる、すなわち、６回分の観測を得るには、例えば、各送信局３ａ〜３ｃが２回ずつ送信すればよい。例えば、時系列に１局ずつ送信するとすれば、３ａ、３ｂ、３ｃ、３ａ、３ｂ、３ｃの順番で送信すればよい。なお、必ずしも１局ずつ時系列に送信する必要はなく、３つの送信局のうち複数を同時送信してもよいし、各送信機で均等に２回ずつにしなくてもよい。また、送信局は３である必要はなく、例えば、１つでもよく、その場合は、異なる時刻で６回観測することになる。ここで大切なのは、異なる時刻における到来時間の差及びドップラー周波数の差の観測値があること、すなわち、衛星の位置及び速度が時々刻々と変化するため、連立方程式を直接位置及び速度について解かず、ケプラリアン軌道６要素について解くことである。
【００４７】
例えば、時系列に１局ずつ３ａ、３ｂ、３ｃ、３ａ、３ｂ、３ｃの順番で送信する場合、ＴＤＯＡ／ＦＤＯＡ計算装置５は、送信局３ａからの１回目の送信によりτ１とｆ１、送信局３ｂからの２回目の送信によりτ２とｆ２、送信局３ｃからの３回目の送信によりτ３とｆ３、送信局３ａからの４回目の送信によりτ４とｆ４、送信局３ｂからの５回目の送信によりτ５とｆ５、送信局３ｃからの６回目の送信によりτ６とｆ６を、時系列に遅延時間の差とドップラー周波数の差としての観測値を得る。
【００４８】
図８は、この発明の実施の形態２に係わる衛星軌道計算装置の機能ブロック図である。
実施の形態２に係わる衛星軌道計算装置６Ｂは、図８に示すように、収束条件、送信局３ａ〜３ｃの位置、電波波長λ、電波伝搬速度ｃなど各種パラメータを入力する設定パラメータ入力部１１Ｂ、収束演算のためのケプラリアン軌道６要素の初期値を設定する軌道初期値設定部１２Ｂ、到来時間の差τｍとドップラー周波数の差ｆｍを用いて収束演算を行う収束演算部１３Ｂ、収束演算により得たケプラリアン軌道６要素が収束条件を満たしているか否かを判定する収束条件判定部１４Ｂを有する。
【００４９】
設定パラメータ入力部１１Ｂは、収束条件のための設定パラメータεＭ０、εｅ、εＡ、εω、εΩ０、εｉ０、収束演算のための３つの送信局３ａ〜３ｂの位置ｐｒ，ｍ（ｍは１〜３）、電波波長λ、電波伝搬速度ｃを入力する。
軌道初期値設定部１２Ｂは、２機の衛星２の初期ケプラリアン軌道６要素ξ（ハット）を設定する。
収束演算部１３Ｂは、式（３２）、式（３３）に関する収束演算を行い、ケプラリアン軌道６要素ξ（ハット）を算出する。このとき、Ｇ（ξ（ハット），ｔｋ）は、式（３４）にξ（ハット）を代入して算出する。
【００５０】
【数１４】
2008064566.tif 000016
【００５１】
また、式（３４）のＧχ→θ，ｋ、Ｇξ→χ，ｋは、それぞれ式（３５）と式（３６）から算出する。
【００５２】
【数１５】
2008064566.tif 000017
【００５３】
式（３５）の▽χτ（χｋ，ｔｋ）と▽χｆ（χｋ，ｔｋ）は、それぞれ式（３７）と式（３８）から算出する。
【００５４】
【数１６】
2008064566.tif 000018
【００５５】
式（３６）の位置に関する各要素は、式（３９）と式（４０）であり、式（３９）と式（４０）の各要素は、式（４１）から算出する。なお、Ωｓｉ，ｋは、式（４２）から算出する。但し、Ωｅ’は地球自転速度定数、ｔｅｒは恒星時である。
【００５６】
【数１７】
2008064566.tif 000019
【００５７】
また、式（３６）の速度に関する各要素は、式（４３）と式（４４）であり、式（４３）と式（４４）の各要素は、式（４５）から算出する。
【００５８】
【数１８】
2008064566.tif 000020
【００５９】
また、式（４１）、式（４５）の▽ξｓｉｘｏｒｂ，ｓｉ，ｋ、▽ξｓｉｙｏｒｂ，ｓｉ，ｋは、式（４６）と式（４７）の要素であり、式（４６）と式（４７）の各要素は、式（４８）から算出する。
【００６０】
【数１９】
2008064566.tif 000021
【００６１】
また、式（４５）の▽ξｓｉｘ（ドット）ｏｒｂ，ｓｉ，ｋ、▽ξｓｉｙ（ドット）ｏｒｂ，ｓｉ，ｋは、式（４９）と式（５０）の要素であり、式（４９）と式（５０）の各要素は、式（５１）から算出する。
【００６２】
【数２０】
2008064566.tif 000022
【００６３】
また、式（４８）、式（５１）の▽ξＥｓ１，ｋ、▽ξＥｓ２，ｋ、は、それぞれ式（５２）、式（５３）であり、式（５２）と式（５３）の各要素は式（５４）から算出する。
【００６４】
【数２１】
2008064566.tif 000023
【００６５】
また、式（５１）の▽ξＥ（ドット）ｓ１，ｋ、▽ξＥ（ドット）ｓ２，ｋは、それぞれ式（５５）、式（５６）であり、式（５５）と式（５６）の各要素は式（５７）から算出する。
【００６６】
【数２２】
2008064566.tif 000024
【００６７】
収束条件判定部１４Ｂは、収束演算部１３Ｂで算出したケプラリアン軌道６要素ξに関するステップ値Δξ（ハット）が式（５８）、式（５９）、式（６０）、式（６１）、式（６２）、式（６３）の条件式を満足するか否かを判定し、満足するときそのケプラリアン軌道６要素を衛星２の軌道とし、満足していないときには再度収束演算部１３Ｂで収束演算を行う。
【００６８】
【数２３】
2008064566.tif 000025
【００６９】
図９は、実施の形態２に係わる衛星軌道計算装置６Ｂでの軌道推定手順を示すフローチャートである。
軌道推定手順を開始すると、図９に示すように、ステップＳ２０１で、収束条件、送信局３の位置など各種パラメータを読み込む。
ステップＳ２０２で、ＴＤＯＡ／ＦＤＯＡ計算装置５から到来時間の差とドップラー周波数の差に関する観測値を読み込む。
ステップＳ２０３で、収束演算のためのケプラリアン軌道６要素の初期値を読み込む。
ステップＳ２０４で、収束演算を式（３２）〜式（５７）に基づいて行い、ケプラリアン軌道６要素ξ（ハット）を推定する。
ステップＳ２０５で、算出したケプラリアン軌道６要素ξ（ハット）について式（５８）から式（６３）の収束条件式を満足しているか否かを判断し、満足しているとき軌道推定手順を終了し、満足していないときステップＳ２０４に戻る。
【００７０】
次に、式（３２）と式（３３）で表されるケプラリアン軌道６要素の推定値ξ（ハット）およびステップ値Δξ（ハット）を収束演算し、式（５８）から式（６３）の収束判定条件式で判定することによりケプラリアン衛星６要素ξを推定できる根拠を説明する。
衛星２の軌道は、ケプラー則と呼ばれる軌道式に基づいて、軌道式の設定パラメータであるケプラリアン軌道６要素（Ｍ０、ｅ、Ａ、ω、Ω０、ｉ０）を決定することにより決定される。ここで、Ｍ０は平均近点離角、ｅは離心率、Ａは軌道長半径、ωは近地点引数、Ω０は昇交点赤径、ｉ０は軌道傾斜角である。
【００７１】
次に、ケプラリアン軌道６要素を推定するアルゴリズムの定式化を説明する。
ξを２機の衛星２のケプラリアン軌道６要素を含むベクトル値とし、式（６４）、式（６５）で定義する。
【００７２】
【数２４】
2008064566.tif 000026
【００７３】
時刻ｔｋにおける到来時間の差およびドップラー周波数の差の観測値をそれぞれτｏｂｓ，ｋ、ｆｏｂｓ，ｋとおく。τｏｂｓ，ｋ、ｆｏｂｓ，ｋは、既知の基準送信局７と、時刻ｔｋに使用する送信局３（位置：ｐｒ，ｋ）により得られたとする。また、時刻ｔｋにおける２機の衛星２の位置および速度を式（６６）とおく。ここで、式（６７）で表されるｐｓｉ，ｋは３次元位置、式（６８）で表されるｖｓｉ，ｋは３次元速度である。
【００７４】
【数２５】
2008064566.tif 000027
【００７５】
ケプラリアン軌道６要素ξと衛星の軌道χは、ケプラーの法則で対応付けられており、χｋ＝χｋ（ξ，ｔｋ）と表すことができる。このとき、時刻ｔｋにおける到来時間の差τ（χｋ，ｔｋ）およびドップラー周波数の差ｆ（χｋ，ｔｋ）の観測モデルはそれぞれ式（６９）、式（７０）で与えられる。
【００７６】
【数２６】
2008064566.tif 000028
【００７７】
ここで、χｋはξの関数、すなわちχｋ＝χｋ（ξ，ｔｋ）であるため、τ（χｋ，ｔｋ）およびｆ（χｋ，ｔｋ）をそれぞれτ（ξ，ｔｋ）およびｆ（ξ，ｔｋ）と表すことができる。
次に、θｏｂｓ，ｋは時刻ｔｋで得られた到来時間の差の観測値τｏｂｓ，ｋおよびドップラー周波数の差の観測値ｆｏｂｓ，ｋを格納したベクトル値で、式（７１）で定義する。
また、ｋを観測回数、θ（ξ，ｔｋ）を、式（６９）および式（７０）に示した数学モデルをベクトル状に表したものとし、式（７２）で定義する。
【００７８】
【数２７】
2008064566.tif 000029
【００７９】
以上の定義から、推定問題を式（７３）のように定式化することができる。
【００８０】
【数２８】
2008064566.tif 000030
【００８１】
式（７３）は、上記評価関数を２機の衛星２のケプラリアン軌道６要素ξについて最小化する最小二乗問題である。
次に、式（７３）の最小二乗問題の解を収束演算で求める方法を説明する。
評価関数式（７３）をＦ（ξ）と呼び、θ（ξ，ｔｋ）の線形一次近似を式（７４）に基づいて行う。ここで、Ｇ（ξ，ｔｋ）は、式（３４）であり、計算式について式（３５）から式（５７）である。
【００８２】
【数２９】
2008064566.tif 000031
【００８３】
すると、式（７３）は、式（７５）となる。ここで、θｔｍｐ，ｋは、式（７６）である。
【００８４】
【数３０】
2008064566.tif 000032
【００８５】
Ｆ（ξ）の偏微分▽Ｆ（ξ）は、式（７７）で表される。
【００８６】
【数３１】
2008064566.tif 000033
【００８７】
そして、式（７７）を０１２×１に等しいとしたものをξについて解くと、式（３２）で表されるケプラリアン軌道６要素の最適解ξ（ハット）と式（３３）で表されるステップ値Δξ（ハット）が得られる。よって、式（３２）および式（３３）の収束演算により最適なケプラリアン軌道６要素を推定することができる。なお、収束判定では、収束演算のステップ値Δξ（ハット）から、同じ種類のケプラリアン軌道６要素を選択してベクトル上に並べた値に対し、式（５８）から式（６３）に基づいて判定する。
【００８８】
このように実施の形態２に係わる移動体軌道推定装置１Ｂは、利用できる送信局３の局数が衛星２の軌道に関する全ての未知変数の解を求めるために必要な局数より少なくても、電波を送信局３から時系列的に送信し、得られた到来時間の差およびドップラー周波数の差を用いてケプラリアン軌道６要素を推定できるので、その得られたケプラリアン軌道６要素の推定値からケプラー則により任意の時刻の衛星軌道を推定することができる。
【００８９】
また、実施の形態２に係わる移動体軌道推定装置１Ｂは、既に放送している放送局を送信局３として利用することができるので、軌道推定のために別途衛星２に電波を照射する必要がなく、構成が簡素化できる。
【００９０】
なお、実施の形態２に係わる移動体軌道推定装置１Ｂは、衛星２に限らず任意の移動体の軌道を推定することができる。
また、実施の形態１と同様に、軌道を推定できる衛星の機数は２機に限るものではなく、任意の機数の衛星の任意の次元の軌道推定も可能である。例えば、移動体の数をＮ機とし、ｄ次元の位置と速度を推定するとき、ケプラーの法則のケプラリアン軌道６要素のように、任意時刻におけるｄ次元の軌道位置と速度を表すことができる数式法則があると仮定する。この法則におけるパラメータ要素が２ｄ個であるとすると、必要な方程式数は２Ｎｄである。次に、送信局数をＭとし、観測タイミング数をＫとする。例えば、各タイミングにおいてＭ局の送信局による同時送信を行うとすると、２Ｍ（Ｎ−１）Ｋ個の方程式が得られる。これらの定義から、方程式の数が未知数を上回るためには、２Ｍ（Ｎ−１）Ｋ≧２Ｎｄを満たす送信局数Ｍおよび観測タイミング数Ｋを用意すればよい。なお、基準送信局を用いる場合には、加えて１局必要となる。収束演算式の導出は、上記の導出に沿って容易に拡張できるため、ここでは説明を省略する。
【００９１】
また、実施の形態２ではＧ（ξ，ｔｋ）を数式を用いて導出しているが、近似的に求めてもよい。近似的に求める場合、まず、Ｇ（ξ，ｔｋ）を式（７８）のように列ベクトルに分解する。このとき、各列ベクトルは式（７９）で表され、Ｕｊは、第ｊ要素だけが１であり、他の要素がすべて零であるような１２次元ベクトルである。εｊはｊに応じて十分小さい値に設定すればよい。
【００９２】
【数３２】
2008064566.tif 000034
【図面の簡単な説明】
【００９３】
【図１】この発明の実施の形態１に係わる移動体軌道推定装置のブロック図である。
【図２】実施の形態１に係わる移動体軌道推定装置により２機の衛星の軌道を推定する様子を示す概念図である。
【図３】実施の形態１に係わる衛星軌道計算装置の機能ブロック図である。
【図４】実施の形態１に係わる衛星軌道計算装置での軌道推定手順を示すフローチャートである。
【図５】実施の形態１に係わる移動体軌道推定装置を未知位置送信局位置評定装置に適用したときの全体のブロック図である。
【図６】この発明の実施の形態２に係わる移動体軌道推定装置のブロック図である。
【図７】実施の形態２に係わる移動体軌道推定装置により２機の衛星の軌道を推定する様子を示す概念図である。
【図８】実施の形態２に係わる衛星軌道計算装置の機能ブロック図である。
【図９】実施の形態２に係わる衛星軌道計算装置での軌道推定手順を示すフローチャートである。
【符号の説明】
【００９４】
１、１Ｂ  移動体軌道推定装置、２  衛星、３、３ａ〜３ｃ  送信局、４  受信局、５  ＴＤＯＡ／ＦＤＯＡ計算装置、６、６Ｂ  衛星軌道計算装置、７  基準送信局、１０  送信局位置評定装置、１１、１１Ｂ  設定パラメータ入力部、１２、１２Ｂ  軌道初期値設定部、１３、１３Ｂ  収束演算部、１４、１４Ｂ  収束条件判定部。
【図１】
2008064566.tif 000035
【図２】
2008064566.tif 000036
【図３】
2008064566.tif 000037
【図４】
2008064566.tif 000038
【図５】
2008064566.tif 000039
【図６】
2008064566.tif 000040
【図７】
2008064566.tif 000041
【図８】
2008064566.tif 000042
【図９】
2008064566.tif 000043