JP2008064566A - 軌道推定方法とその装置 - Google Patents

Abstract

【課題】既設置の発信局から移動体に到達する電波を用いることにより新規に送受信機を必要としない軌道推定方法とその装置を提供する。
【解決手段】軌道推定方法は、複数の送信局から複数の移動体に電波を同時に送信し、２局の受信局において、上記移動体に搭載された中継器により中継送信された電波を受信し、上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出し、上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることのできる局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差から直接上記移動体の軌道を推定する。
【選択図】図３

Description

この発明は、移動体の軌道を推定する方法およびその装置に関するものである。

位置が未知の不法電波送信局などの位置を評定する方法では、２機の衛星を介して得られる電波の到来時間差（Ｔｉｍｅ Ｄｉｆｆｅｒｅｎｃｅ ｏｆ Ａｒｒｉｖａｌ：ＴＤＯＡ）とドップラー周波数差（Ｆｒｅｑｕｅｎｃｙ Ｄｉｆｆｅｒｅｎｃｅ ｏｆ Ａｒｒｉｖａｌ：ＦＤＯＡ）情報を用いて、位置を評定する方法が提案されている（例えば、非特許文献２参照）。この方法では、位置が未知の送信局から衛星に向けて送信される送信電波は、メインローブの他にサイドローブなど、空間的に幅をもっているので、２機の衛星により受信することができる。そして、衛星は中継器を搭載しており、電波を地上の受信局へ送り返す中継送信機能を持っているため、受信局では衛星を介して送り返された電波を受信することができる。２局の受信局の受信信号処理を行うことにより、信号の到来時間差およびドップラー周波数差を得ることができる。

この位置が未知な送信局の位置を評定するためには、衛星の位置および速度からなる軌道情報が必要となる。そこで、衛星の軌道を推定する従来の方法として、衛星の距離や方位角を観測し、その情報を用いて軌道推定を行っている。この方式では、光学望遠鏡により衛星からの光を観測したり、衛星との間で電波の送受信することにより、観測局から見た衛星の距離、仰角および方位角を測定したりしている（例えば、非特許文献１参照）。

有本好徳、他２名、「静止衛星測位システムにおける軌道決定精度について」、宇宙航行エレクトロニクス研究会資料、社団法人電子情報通信学会、１９９７年、ＳＡＮＥ８９−２６、ｐ．４５−５２ Ｄ．Ｐ．Ｈａｗｏｒｔｈ、他３名、「Ｉｎｔｅｒｆｅｒｅｎｃｅ Ｌｏｃａｌｉｚａｔｉｏｎ Ｆｏｒ ＥＵＴＥＬＳＡＴ Ｓａｔｅｌｌｉｔｅｓ−Ｔｈｅ Ｆｉｒｓｔ Ｅｕｒｏｐｅａｎ Ｔｒａｎｓｍｉｔｔｅｒ Ｌｏｃａｔｉｏｎ Ｓｙｓｔｅｍ」、Ｉｎｔｅｒｎａｔｉｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｓａｔｅｌｌｉｔｅ Ｃｏｍｍｕｎｉｃａｔｉｏｎｓ、Ｖｏｌ．１５、ｐ．１５５−１８３

しかし、衛星からの光を観測するためには、大型の光学望遠鏡を設置しなければならないし、衛星との間で電波の送受信を行う送受信機を別途必要になったりするので、装置全体の規模が大きくなるという問題がある。

この発明は、既設置の発信局から移動体に到達する電波を用いることにより新規に送受信機を必要としない軌道推定方法とその装置を提供することである。

この発明に係わる軌道推定方法は、複数の送信局から複数の移動体に電波を同時に送信し、２局の受信局において、上記移動体に搭載された中継器により中継送信された電波を受信し、上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出し、上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることのできる局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差から直接上記移動体の軌道を推定するものである。

この発明に係わる移動体軌道推定によれば、位置、仰角、方位角情報ではなく、到来時間差とドップラー周波数差情報を用いるので、放送局などの正規な送信局からの電波を用いることができ、軌道推定のために新たに電波を送信する必要がなく、その結果、軌道推定用に送信局を設置する必要がない。

実施の形態１．
図１は、この発明の実施の形態１に係わる移動体軌道推定装置のブロック図である。図２は、実施の形態１に係わる移動体軌道推定装置により２機の衛星の軌道を推定する様子を示す概念図である。
この発明の実施の形態１に係わる移動体軌道推定装置１は、図１、２に示すように、２機の軌道が未知の衛星２の軌道を推定するために、６個の到来時間の差τ_ｍ（ｍは１〜６の整数）と６個のドップラー周波数の差ｆ_ｍ（ｍは１〜６の整数）を用いているので、位置ｐ_ｒ，ｍが既知の６つの送信局３、１つの基準送信局７、各送信局３、７から各衛星２に到来し、各衛星２で中継された電波を受信する２つの受信局４、到来時間の差τ_ｍとドップラー周波数の差ｆ_ｍを算出するＴＤＯＡ／ＦＤＯＡ計算装置５、到来時間の差τ_ｍとドップラー周波数の差ｆ_ｍを用いて各衛星２の位置ｐ_ｓｉ（ｉは衛星を示し、１または２である）および衛星２の速度ｖ_ｓｉを推定する衛星軌道計算装置６を備える。この送信局３、７にはそれぞれ１台の発信機が設置され、受信局４にはそれぞれ１台の受信機が設置されている。なお、到来時間の差τ_ｍとドップラー周波数の差ｆ_ｍの算出は、一般的な方法であるので、説明は省略する。

図３は、実施の形態１に係わる衛星軌道計算装置６の機能ブロック図である。
実施の形態１に係わる衛星軌道計算装置６は、図３に示すように、収束条件、送信局３、７の位置、電波波長λ、電波伝搬速度ｃなど各種パラメータを入力する設定パラメータ入力部１１、収束演算のための軌道の初期値を設定する軌道初期値設定部１２、到来時間の差τ_ｍとドップラー周波数の差ｆ_ｍを用いて収束演算を行う収束演算部１３、収束演算により得た軌道が収束条件を満たしているか否かを判定する収束条件判定部１４を有する。

設定パラメータ入力部１１は、収束条件のための設定パラメータε_ｐ、ε_ｖ、収束演算のための６つの送信局３の位置ｐ_ｒ，ｍ（ｍは１〜６）、基準送信局７の位置ｐ_ｒ０、電波波長λ、電波伝搬速度ｃを入力する。
軌道初期値設定部１２は、２機の衛星２の初期軌道χ（ハット）を設定する。
収束演算部１３は、式（１）、式（２）に関する収束演算を行い、軌道χを算出する。このとき、Ｇ（χ（ハット）、ｍ）は、式（３）にχ（ハット）を代入して算出する。

また、式（３）の各要素は式（４）から式（９）により算出する。

収束条件判定部１４は、収束演算部１３で算出した軌道χに関するステップ値Δχ（ハット）が式（１０）、式（１１）の条件式を満足するか否かを判定し、満足するときその軌道を衛星２の軌道とし、満足していないときには再度収束演算部１３で収束演算を行う。

図４は、実施の形態１に係わる衛星軌道計算装置６での軌道推定手順を示すフローチャートである。
軌道推定手順を開始すると、図４に示すように、ステップＳ１０１で、収束条件、送信局３の位置、基準送信局７の位置など各種パラメータを読み込む。
ステップＳ１０２で、ＴＤＯＡ／ＦＤＯＡ計算装置５から到来時間の差とドップラー周波数の差に関する観測値を読み込む。
ステップＳ１０３で、収束演算のための軌道初期値を読み込む。
ステップＳ１０４で、収束演算を式（１）〜式（９）に基づいて行い、軌道χ（ハット）を推定する。
ステップＳ１０５で、算出した軌道χ（ハット）について式（１０）、式（１１）の収束条件式を満足しているか否かを判断し、満足しているとき軌道推定手順を終了し、満足していないときステップＳ１０４に戻る。

次に、式（１）、式（２）で表される軌道の推定値χ（ハット）およびステップ値Δχ（ハット）を収束演算し、式（１０）および式（１１）の収束判定条件式で判定することにより衛星軌道χを推定できる根拠を説明する。
位置が既知の６つの送信局３から２機の衛星２に同時に送信を行い、衛星２を経由して得られる信号を各受信局４で受信し、到来時間の差τ_ｍおよびドップラー周波数の差ｆ_ｍを計算することにより、第１送信局３から第７送信局３それぞれから衛星２への到来時間の差τ_ｍおよびドップラー周波数の差ｆ_ｍを観測することができる。

送信局３の３次元位置ｐ_ｒ，ｍを式（１２）、２機の衛星２の３次元位置および速度を式（１３）で表す。ここで、式（１４）のｐ_ｓｉは衛星２の３次元位置で、式（１５）のｖ_ｓｉは衛星２の３次元速度である。
このとき、到来時間の差τ_ｍおよびドップラー周波数の差ｆ_ｍは、それぞれ式（１６）、式（１７）で表される。

ここで、｜｜・｜｜は、ベクトルのノルムを表し、物理的には距離を意味する。λは電波波長、ｃは電波伝播速度である。式（１６）は、送信局３と衛星２の距離の差を電波伝播速度ｃで除算した量の差、すなわち、２機の衛星２への電波の到来時間の差を意味する。
また、式（１７）は、送信局３と衛星２とを結ぶ線上への衛星速度の射影成分を電波波長λで除算した量の差、すなわち、２機の衛星２のドップラー周波数の差を意味する。

式（１６）および式（１７）に示すように、２機の衛星２の３次元位置および速度の未知変数は計１２個あり、１つの送信局３による送信で、到来時間の差およびドップラー周波数の差の２本の方程式を立てることができる。よって、６つ以上の送信局３を用いて同時に送信を行い、到来時間の差およびドップラー周波数の差を用いてそれぞれ６本以上、計１２本以上の方程式を立てれば、２機の衛星２の３次元位置ｐ_ｓｉおよび速度ｖ_ｓｉを求めることができる。

また、実際の到来時間の差およびドップラー周波数の差の観測値には受信局４の時計時間誤差、受信局４でのローカル周波数シフト、衛星内部でのローカル周波数シフトなどの誤差要因が存在するので、これらをキャンセルするために既知位置の基準送信局７（位置：ｐ_ｒ０）を用意し、基準送信局７による送信も同時に行う。そして、基準送信局７による到来時間の差およびドップラー周波数の差を基準とした到来時間の差およびドップラー周波数の差を用いることにする。すなわち、式（１６）および式（１７）の代わりに式（１８）および式（１９）を用いる。

このように、６つ以上の送信局３と基準送信局７の計７局以上で２機の衛星２の３次元位置ｐ_ｓｉおよび３次元速度ｖ_ｓｉを推定することができる。
次に、式（１８）および式（１９）の連立方程式を衛星２の位置ｐ_ｓｉおよび速度ｖ_ｓｉについて解く方法を説明する。
第ｍ送信局３による到来時間の差およびドップラー周波数の差の観測値をそれぞれτ_{ｏｂｓ，ｍ}、ｆ_{ｏｂｓ，ｍ}とおき、式（２０）のようにベクトル状に表記する。
また、θ_ｍ（χ）は、式（１８）および式（１９）で表す到来時間の差τ_ｍ（χ）およびドップラー周波数の差ｆ_ｍ（χ）をベクトル状に表したものであり、式（２１）のように定義できる。

式（２２）に表される観測誤差共分散行列Ｖは、受信機雑音などによる到来時間差観測雑音の分散値σ_τ ^２およびドップラー周波数差観測雑音の分散値σ_ｆ ^２を対角項にもつ。

このように定義された値を用いると、推定問題を式（２３）のように定式化できる。
すなわち、式（２３）を評価関数として、衛星の軌道χについて最小化する最小二乗問題となる。

次に、式（２３）を評価関数とする最小二乗法の解を収束演算で求める方法を説明する。今、評価関数式（２３）をＦ（χ）と呼ぶ。まず、式（２４）に従って、θ_ｍ（χ）の線形一次近似を行う。

ここで、Ｇ（χ，ｍ）は、式（２５）に示すようにθｍ（χ）のχについての偏微分であり、ベクトル表示すると式（２６）で表される。そして、Ｇ（χ，ｍ）を数式により計算する場合、Ｇ（χ，ｍ）は式（３）で表される。そして、▽_ｐｓ１τ_ｍ（χ）、▽_ｐｓ２τ_ｍ（χ）、▽_ｐｓ１ｆ_ｍ（χ）、▽_ｐｓ１ｆ_ｍ（χ）、▽_ｖｓ１ｆ_ｍ（χ）、▽_ｖｓ１ｆ_ｍ（χ）は、それぞれ式（４）、式（５）、式（６）、式（７）、式（８）、式（９）で求められる。

すると、評価関数Ｆ（χ）は式（２７）で表される。ここで、θ_{ｔｍｐ，ｍ}は式（２８）で表される。

評価関数Ｆ（χ）のχについての偏微分▽_χＦ（χ）は、式（２９）で表される。ここで▽_χは、χについての偏微分を意味する。

偏微分値▽_χＦ（χ）＝０_１２×１としたものをχについて解くと最適解が得られ、このとき、χ（ハット）は式（１）、収束演算のステップ値Δχ（ハット）は式（２）で表される。よって、式（１）および式（２）の収束演算により最適な衛星の軌道χを推定することができる。なお、収束判定では、収束演算のステップ値Δχ（ハット）から、位置および速度に相当する要素をそれぞれ選択してベクトル上に並べた値に対し、式（１０）、式（１１）の判定条件で判定する。

なお、実施の形態１に係わる移動体軌道推定装置１は、２機の衛星２の軌道を推定したが、衛星２の数はこれに限るものではない。例えば、衛星２の数をＮ機とし、ｄ次元の衛星の位置と速度とを推定する場合、必要な方程式数は２×Ｎ×ｄである。次に、送信局数をＭとし、Ｍ局の送信局３から同時送信することにより、２Ｍ（Ｎ−１）の方程式が得られる。これらから、方程式の数が未知数を上回るためには、２Ｍ（Ｎ−１）が２×Ｎ×ｄ以上を満たす送信局数Ｍを用意すればよい。収束演算式の導出は、実施の形態１で説明した軌道推定方法と同様である。なお、基準送信局を用いる場合には、加えて１局必要となる。

このように位置が既知の既に電波を出射している放送局などの送信局３を用いて衛星２の軌道を推定することができるので、軌道推定のためだけに別途衛星２に電波を照射する必要がなく、移動体軌道推定装置１の新規な送信局の設置は不要となる。
なお、実施の形態１に係わる移動体軌道推定装置１は、衛星２に限らず任意の移動体の軌道を推定することができる。

図５は、送信局位置評定装置に移動体軌道推定装置１を適用したときの構成図を示す。
次に、非特許文献２に記載されている位置が未知な送信局の位置を評定する送信局位置評定装置に実施の形態１に係わる移動体軌道推定装置１を適用することを説明する。
このときは、図５に示すように、軌道推定用受信局４およびＴＤＯＡ／ＦＤＯＡ計算装置５は、送信局位置評定装置１０の位置評定用受信局およびＴＤＯＡ／ＦＤＯＡ計算装置として兼用することができ、構成が簡単な送信局位置評定装置を提供することができる。

なお、実施の形態１ではＧ（χ，ｍ）を数式を用いて導出しているが、近似的に求めてもよい。近似的に求める場合、まず、Ｇ（χ，ｍ）を式（３０）のように列ベクトルに分解する。このとき、各列ベクトルは式（３１）で表され、Ｕ_ｊは、第ｊ要素だけが１であり、他の要素がすべて零であるような１２次元ベクトルである。ε_ｊはｊに応じて十分小さい値に設定すればよい。

実施の形態２．
図６は、この発明の実施の形態２に係わる移動体軌道推定装置のブロック図である。図７は、実施の形態２に係わる移動体軌道推定装置により２機の衛星の軌道を推定する様子を示す概念図である。
この発明の実施の形態２に係わる移動体軌道推定装置１Ｂは、実施の形態１に係わる移動体軌道推定装置１と送信局３ａ〜３ｃの局数が必ずしも６局以上ではなく、それに伴って衛星軌道計算装置６Ｂが異なっており、それ以外は同様であるので、同様な部分に同じ符号を付記して説明は省略する。
実施の形態１に係わる移動体軌道推定装置１では、２機の衛星２の１２個の未知変数の解を求めるために、６つ以上の位置が既知の送信局３から同時に送信されている電波を用いて、１２本以上の方程式を立てている。
しかし、実際に送信局３を多数設置することが難しい場合がある。そこで、実施の形態２に係わる移動体軌道推定装置１Ｂでは、観測回数を増やすことにより必要な数以上の方程式をたてる。この場合、衛星２の位置および速度は時々刻々と変化するため、衛星２の位置および速度について、直接連立方程式を解くことはできない。そこで、衛星２の位置および速度を未知変数とするのではなく、ケプラリアン軌道６要素を未知変数として推定する。すなわち、２機の衛星２のケプラリアン軌道６要素のための１２本の方程式を立てる。

実施の形態２においては、図７に示すように、３つの送信局３ａ〜３ｃが設置されているので、１２本の方程式を立てる、すなわち、６回分の観測を得るには、例えば、各送信局３ａ〜３ｃが２回ずつ送信すればよい。例えば、時系列に１局ずつ送信するとすれば、３ａ、３ｂ、３ｃ、３ａ、３ｂ、３ｃの順番で送信すればよい。なお、必ずしも１局ずつ時系列に送信する必要はなく、３つの送信局のうち複数を同時送信してもよいし、各送信機で均等に２回ずつにしなくてもよい。また、送信局は３である必要はなく、例えば、１つでもよく、その場合は、異なる時刻で６回観測することになる。ここで大切なのは、異なる時刻における到来時間の差及びドップラー周波数の差の観測値があること、すなわち、衛星の位置及び速度が時々刻々と変化するため、連立方程式を直接位置及び速度について解かず、ケプラリアン軌道６要素について解くことである。

例えば、時系列に１局ずつ３ａ、３ｂ、３ｃ、３ａ、３ｂ、３ｃの順番で送信する場合、ＴＤＯＡ／ＦＤＯＡ計算装置５は、送信局３ａからの１回目の送信によりτ_１とｆ_１、送信局３ｂからの２回目の送信によりτ_２とｆ_２、送信局３ｃからの３回目の送信によりτ_３とｆ_３、送信局３ａからの４回目の送信によりτ_４とｆ_４、送信局３ｂからの５回目の送信によりτ_５とｆ_５、送信局３ｃからの６回目の送信によりτ_６とｆ_６を、時系列に遅延時間の差とドップラー周波数の差としての観測値を得る。

図８は、この発明の実施の形態２に係わる衛星軌道計算装置の機能ブロック図である。
実施の形態２に係わる衛星軌道計算装置６Ｂは、図８に示すように、収束条件、送信局３ａ〜３ｃの位置、電波波長λ、電波伝搬速度ｃなど各種パラメータを入力する設定パラメータ入力部１１Ｂ、収束演算のためのケプラリアン軌道６要素の初期値を設定する軌道初期値設定部１２Ｂ、到来時間の差τ_ｍとドップラー周波数の差ｆ_ｍを用いて収束演算を行う収束演算部１３Ｂ、収束演算により得たケプラリアン軌道６要素が収束条件を満たしているか否かを判定する収束条件判定部１４Ｂを有する。

設定パラメータ入力部１１Ｂは、収束条件のための設定パラメータε_Ｍ０、ε_ｅ、ε_Ａ、ε_ω、ε_Ω０、ε_ｉ０、収束演算のための３つの送信局３ａ〜３ｂの位置ｐ_ｒ，ｍ（ｍは１〜３）、電波波長λ、電波伝搬速度ｃを入力する。
軌道初期値設定部１２Ｂは、２機の衛星２の初期ケプラリアン軌道６要素ξ（ハット）を設定する。
収束演算部１３Ｂは、式（３２）、式（３３）に関する収束演算を行い、ケプラリアン軌道６要素ξ（ハット）を算出する。このとき、Ｇ（ξ（ハット），ｔ_ｋ）は、式（３４）にξ（ハット）を代入して算出する。

また、式（３４）のＧ_{χ→θ，ｋ}、Ｇ_{ξ→χ，ｋ}は、それぞれ式（３５）と式（３６）から算出する。

式（３５）の▽_χτ（χ_ｋ，ｔ_ｋ）と▽_χｆ（χ_ｋ，ｔ_ｋ）は、それぞれ式（３７）と式（３８）から算出する。

式（３６）の位置に関する各要素は、式（３９）と式（４０）であり、式（３９）と式（４０）の各要素は、式（４１）から算出する。なお、Ω_ｓｉ，ｋは、式（４２）から算出する。但し、Ω_ｅ’は地球自転速度定数、ｔ_ｅｒは恒星時である。

また、式（３６）の速度に関する各要素は、式（４３）と式（４４）であり、式（４３）と式（４４）の各要素は、式（４５）から算出する。

また、式（４１）、式（４５）の▽_ξｓｉｘ_{ｏｒｂ，ｓｉ，ｋ}、▽_ξｓｉｙ_{ｏｒｂ，ｓｉ，ｋ}は、式（４６）と式（４７）の要素であり、式（４６）と式（４７）の各要素は、式（４８）から算出する。

また、式（４５）の▽_ξｓｉｘ（ドット）_{ｏｒｂ，ｓｉ，ｋ}、▽_ξｓｉｙ（ドット）_{ｏｒｂ，ｓｉ，ｋ}は、式（４９）と式（５０）の要素であり、式（４９）と式（５０）の各要素は、式（５１）から算出する。

また、式（４８）、式（５１）の▽_ξＥ_ｓ１，ｋ、▽_ξＥ_ｓ２，ｋ、は、それぞれ式（５２）、式（５３）であり、式（５２）と式（５３）の各要素は式（５４）から算出する。

また、式（５１）の▽_ξＥ（ドット）_ｓ１，ｋ、▽_ξＥ（ドット）_ｓ２，ｋは、それぞれ式（５５）、式（５６）であり、式（５５）と式（５６）の各要素は式（５７）から算出する。

収束条件判定部１４Ｂは、収束演算部１３Ｂで算出したケプラリアン軌道６要素ξに関するステップ値Δξ（ハット）が式（５８）、式（５９）、式（６０）、式（６１）、式（６２）、式（６３）の条件式を満足するか否かを判定し、満足するときそのケプラリアン軌道６要素を衛星２の軌道とし、満足していないときには再度収束演算部１３Ｂで収束演算を行う。

図９は、実施の形態２に係わる衛星軌道計算装置６Ｂでの軌道推定手順を示すフローチャートである。
軌道推定手順を開始すると、図９に示すように、ステップＳ２０１で、収束条件、送信局３の位置など各種パラメータを読み込む。
ステップＳ２０２で、ＴＤＯＡ／ＦＤＯＡ計算装置５から到来時間の差とドップラー周波数の差に関する観測値を読み込む。
ステップＳ２０３で、収束演算のためのケプラリアン軌道６要素の初期値を読み込む。
ステップＳ２０４で、収束演算を式（３２）〜式（５７）に基づいて行い、ケプラリアン軌道６要素ξ（ハット）を推定する。
ステップＳ２０５で、算出したケプラリアン軌道６要素ξ（ハット）について式（５８）から式（６３）の収束条件式を満足しているか否かを判断し、満足しているとき軌道推定手順を終了し、満足していないときステップＳ２０４に戻る。

次に、式（３２）と式（３３）で表されるケプラリアン軌道６要素の推定値ξ（ハット）およびステップ値Δξ（ハット）を収束演算し、式（５８）から式（６３）の収束判定条件式で判定することによりケプラリアン衛星６要素ξを推定できる根拠を説明する。
衛星２の軌道は、ケプラー則と呼ばれる軌道式に基づいて、軌道式の設定パラメータであるケプラリアン軌道６要素（Ｍ_０、ｅ、Ａ、ω、Ω_０、ｉ_０）を決定することにより決定される。ここで、Ｍ_０は平均近点離角、ｅは離心率、Ａは軌道長半径、ωは近地点引数、Ω_０は昇交点赤径、ｉ_０は軌道傾斜角である。

次に、ケプラリアン軌道６要素を推定するアルゴリズムの定式化を説明する。
ξを２機の衛星２のケプラリアン軌道６要素を含むベクトル値とし、式（６４）、式（６５）で定義する。

時刻ｔ_ｋにおける到来時間の差およびドップラー周波数の差の観測値をそれぞれτ_{ｏｂｓ，ｋ}、ｆ_{ｏｂｓ，ｋ}とおく。τ_{ｏｂｓ，ｋ}、ｆ_{ｏｂｓ，ｋ}は、既知の基準送信局７と、時刻ｔ_ｋに使用する送信局３（位置：ｐ_ｒ，ｋ）により得られたとする。また、時刻ｔ_ｋにおける２機の衛星２の位置および速度を式（６６）とおく。ここで、式（６７）で表されるｐ_ｓｉ，ｋは３次元位置、式（６８）で表されるｖ_ｓｉ，ｋは３次元速度である。

ケプラリアン軌道６要素ξと衛星の軌道χは、ケプラーの法則で対応付けられており、χ_ｋ＝χ_ｋ（ξ，ｔ_ｋ）と表すことができる。このとき、時刻ｔ_ｋにおける到来時間の差τ（χ_ｋ，ｔ_ｋ）およびドップラー周波数の差ｆ（χ_ｋ，ｔ_ｋ）の観測モデルはそれぞれ式（６９）、式（７０）で与えられる。

ここで、χ_ｋはξの関数、すなわちχ_ｋ＝χ_ｋ（ξ，ｔ_ｋ）であるため、τ（χ_ｋ，ｔ_ｋ）およびｆ（χ_ｋ，ｔ_ｋ）をそれぞれτ（ξ，ｔ_ｋ）およびｆ（ξ，ｔ_ｋ）と表すことができる。
次に、θ_{ｏｂｓ，ｋ}は時刻ｔ_ｋで得られた到来時間の差の観測値τ_{ｏｂｓ，ｋ}およびドップラー周波数の差の観測値ｆ_{ｏｂｓ，ｋ}を格納したベクトル値で、式（７１）で定義する。
また、ｋを観測回数、θ（ξ，ｔ_ｋ）を、式（６９）および式（７０）に示した数学モデルをベクトル状に表したものとし、式（７２）で定義する。

以上の定義から、推定問題を式（７３）のように定式化することができる。

式（７３）は、上記評価関数を２機の衛星２のケプラリアン軌道６要素ξについて最小化する最小二乗問題である。
次に、式（７３）の最小二乗問題の解を収束演算で求める方法を説明する。
評価関数式（７３）をＦ（ξ）と呼び、θ（ξ，ｔ_ｋ）の線形一次近似を式（７４）に基づいて行う。ここで、Ｇ（ξ，ｔ_ｋ）は、式（３４）であり、計算式について式（３５）から式（５７）である。

すると、式（７３）は、式（７５）となる。ここで、θ_{ｔｍｐ，ｋ}は、式（７６）である。

Ｆ（ξ）の偏微分▽Ｆ（ξ）は、式（７７）で表される。

そして、式（７７）を０_１２×１に等しいとしたものをξについて解くと、式（３２）で表されるケプラリアン軌道６要素の最適解ξ（ハット）と式（３３）で表されるステップ値Δξ（ハット）が得られる。よって、式（３２）および式（３３）の収束演算により最適なケプラリアン軌道６要素を推定することができる。なお、収束判定では、収束演算のステップ値Δξ（ハット）から、同じ種類のケプラリアン軌道６要素を選択してベクトル上に並べた値に対し、式（５８）から式（６３）に基づいて判定する。

このように実施の形態２に係わる移動体軌道推定装置１Ｂは、利用できる送信局３の局数が衛星２の軌道に関する全ての未知変数の解を求めるために必要な局数より少なくても、電波を送信局３から時系列的に送信し、得られた到来時間の差およびドップラー周波数の差を用いてケプラリアン軌道６要素を推定できるので、その得られたケプラリアン軌道６要素の推定値からケプラー則により任意の時刻の衛星軌道を推定することができる。

また、実施の形態２に係わる移動体軌道推定装置１Ｂは、既に放送している放送局を送信局３として利用することができるので、軌道推定のために別途衛星２に電波を照射する必要がなく、構成が簡素化できる。

なお、実施の形態２に係わる移動体軌道推定装置１Ｂは、衛星２に限らず任意の移動体の軌道を推定することができる。
また、実施の形態１と同様に、軌道を推定できる衛星の機数は２機に限るものではなく、任意の機数の衛星の任意の次元の軌道推定も可能である。例えば、移動体の数をＮ機とし、ｄ次元の位置と速度を推定するとき、ケプラーの法則のケプラリアン軌道６要素のように、任意時刻におけるｄ次元の軌道位置と速度を表すことができる数式法則があると仮定する。この法則におけるパラメータ要素が２ｄ個であるとすると、必要な方程式数は２Ｎｄである。次に、送信局数をＭとし、観測タイミング数をＫとする。例えば、各タイミングにおいてＭ局の送信局による同時送信を行うとすると、２Ｍ（Ｎ−１）Ｋ個の方程式が得られる。これらの定義から、方程式の数が未知数を上回るためには、２Ｍ（Ｎ−１）Ｋ≧２Ｎｄを満たす送信局数Ｍおよび観測タイミング数Ｋを用意すればよい。なお、基準送信局を用いる場合には、加えて１局必要となる。収束演算式の導出は、上記の導出に沿って容易に拡張できるため、ここでは説明を省略する。

また、実施の形態２ではＧ（ξ，ｔ_ｋ）を数式を用いて導出しているが、近似的に求めてもよい。近似的に求める場合、まず、Ｇ（ξ，ｔ_ｋ）を式（７８）のように列ベクトルに分解する。このとき、各列ベクトルは式（７９）で表され、Ｕ_ｊは、第ｊ要素だけが１であり、他の要素がすべて零であるような１２次元ベクトルである。ε_ｊはｊに応じて十分小さい値に設定すればよい。

この発明の実施の形態１に係わる移動体軌道推定装置のブロック図である。 実施の形態１に係わる移動体軌道推定装置により２機の衛星の軌道を推定する様子を示す概念図である。 実施の形態１に係わる衛星軌道計算装置の機能ブロック図である。 実施の形態１に係わる衛星軌道計算装置での軌道推定手順を示すフローチャートである。 実施の形態１に係わる移動体軌道推定装置を未知位置送信局位置評定装置に適用したときの全体のブロック図である。 この発明の実施の形態２に係わる移動体軌道推定装置のブロック図である。 実施の形態２に係わる移動体軌道推定装置により２機の衛星の軌道を推定する様子を示す概念図である。 実施の形態２に係わる衛星軌道計算装置の機能ブロック図である。 実施の形態２に係わる衛星軌道計算装置での軌道推定手順を示すフローチャートである。

符号の説明

１、１Ｂ 移動体軌道推定装置、２ 衛星、３、３ａ〜３ｃ 送信局、４ 受信局、５ ＴＤＯＡ／ＦＤＯＡ計算装置、６、６Ｂ 衛星軌道計算装置、７ 基準送信局、１０ 送信局位置評定装置、１１、１１Ｂ 設定パラメータ入力部、１２、１２Ｂ 軌道初期値設定部、１３、１３Ｂ 収束演算部、１４、１４Ｂ 収束条件判定部。

Claims (6)

1. 複数の送信局から複数の移動体に電波を同時に送信し、
   ２局の受信局において、上記移動体に搭載された中継器により中継送信された電波を受信し、
   上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出し、
   上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることのできる局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差から直接上記移動体の軌道を推定することを特徴とする軌道推定方法。
2. 少なくとも１つの送信局から複数の移動体に電波を時系列に送信し、
   ２局の受信局において、上記移動体に搭載された中継器により中継送信された電波を受信し、
   上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出して、
   上記送信局の同時の送信では上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることができない局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差からケプラリアン軌道６要素を推定することを特徴とする軌道推定方法。
3. 複数の移動体に電波を同時に送信する複数の送信局と、
   上記移動体に搭載された中継器により中継送信された電波を受信する複数の受信局と、
   上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出するＴＤＯＡ／ＦＤＯＡ計算装置と、
   上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることのできる局数の送信局が存在するとき、得られた到来時間の差およびドップラー周波数の差から直接上記移動体の軌道を推定する移動体軌道計算装置と、
   を備えることを特徴とする移動体軌道推定装置。
4. 複数の移動体に電波を時系列に送信する少なくとも１つの送信局と、
   上記移動体に搭載された中継器により中継送信された電波を受信する複数の受信局と、
   上記受信した電波から上記移動体に係わる到来時間の差およびドップラー周波数の差を検出するＴＤＯＡ／ＦＤＯＡ計算装置と、
   上記送信局の同時の送信では上記複数の移動体の軌道を示す未知変数以上の到来時間の差およびドップラー周波数の差を得ることができない局数の送信局が存在するときにも、得られた到来時間の差およびドップラー周波数の差からケプラリアン軌道６要素を推定し、上記推定されたケプラリアン軌道６要素から任意の時刻の上記移動体の３次元の位置および３次元の速度を推定することを特徴とする移動体軌道推定装置。
5. 上記移動体の機数は２であり、上記送信局の局数は６以上であり、
   上記移動体軌道計算装置は、同時に送信された電波により得られた到来時間の差およびドップラー周波数の差より、上記移動体の３次元の位置および３次元の速度を推定することを特徴とする請求項３に記載する移動体軌道推定装置。
6. 上記移動体の機数は２であり、上記送信局の局数は１局以上であり、
   上記移動体軌道計算装置は、時系列に送信された電波により得られた到来時間の差およびドップラー周波数の差より、上記移動体のケプラリアン軌道６要素を推定し、上記推定されたケプラリアン軌道６要素から任意の時刻の上記移動体の３次元の位置および３次元の速度を推定することを特徴とする請求項４に記載する移動体軌道推定装置。

Images