知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2024-11-28
(45)【発行日】2024-12-06
(54)【発明の名称】位置標定装置
(51)【国際特許分類】
G01S 5/06 20060101AFI20241129BHJP
G01S 19/40 20100101ALI20241129BHJP
【ＦＩ】
G01S5/06
G01S19/40
【請求項の数】
13
(21)【出願番号】P 2020165988
(22)【出願日】2020-09-30
(65)【公開番号】P2022057630
(43)【公開日】2022-04-11
【審査請求日】2023-08-29
(73)【特許権者】
【識別番号】000006013
【氏名又は名称】三菱電機株式会社
(74)【代理人】
【識別番号】100095407
【弁理士】
【氏名又は名称】木村 満
(74)【代理人】
【識別番号】100131152
【弁理士】
【氏名又は名称】八島 耕司
(74)【代理人】
【識別番号】100147924
【弁理士】
【氏名又は名称】美恵 英樹
(74)【代理人】
【識別番号】100148149
【弁理士】
【氏名又は名称】渡邉 幸男
(74)【代理人】
【識別番号】100181618
【弁理士】
【氏名又は名称】宮脇 良平
(74)【代理人】
【識別番号】100174388
【弁理士】
【氏名又は名称】龍竹 史朗
(72)【発明者】
【氏名】福島 浩文
(72)【発明者】
【氏名】古橋 秀之
(72)【発明者】
【氏名】井上 透
【審査官】山下 雅人
(56)【参考文献】
【文献】特開２０１０－０６０３０３（ＪＰ，Ａ）
【文献】特開２００７－１０７９２９（ＪＰ，Ａ）
【文献】特表２０００－５１２３８０（ＪＰ，Ａ）
【文献】特開２００９－１９８４３５（ＪＰ，Ａ）
【文献】国際公開第２０１７／１０９９５１（ＷＯ，Ａ１）
【文献】特許第６６２５２９５（ＪＰ，Ｂ１）
【文献】特許第６２６１８３１（ＪＰ，Ｂ１）
【文献】特開２００５－２８３１８７（ＪＰ，Ａ）
【文献】米国特許出願公開第２００７／０１２４８２４（ＵＳ，Ａ１）
【文献】米国特許第０５５９４４５２（ＵＳ，Ａ）
【文献】特開２０１９－０９５２３４（ＪＰ，Ａ）
【文献】特開２０１８－００４６０１（ＪＰ，Ａ）
【文献】特開２０２０－０４１９５２（ＪＰ，Ａ）
【文献】特開２００６－３４９４７０（ＪＰ，Ａ）
【文献】国際公開第２０１９／１６７２６８（ＷＯ，Ａ１）
【文献】国際公開第２０１８／１４６７２３（ＷＯ，Ａ１）
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０１Ｓ ５／００－ ５／１４
Ｇ０１Ｓ１９／００－１９／５５
(57)【特許請求の範囲】
【請求項1】
複素数で表される第１の信号および第２の信号、前記第１の信号と前記第２の信号とのそれぞれ時間差および周波数差が入力されて、前記第１の信号に対して前記時間差を持たせた前記第２の信号の複素共役をとり、かつ前記周波数差で回転する絶対値が１である複素数を乗算して得られる信号と前記第１の信号とを乗算した値を決められた時間だけ積分して、前記第１の信号と前記第２の信号との間の２次元相関値を計算する相関関数計算部と、前記時間差と前記周波数差を変化させて前記２次元相関値を計算し、前記２次元相関値が最大のピークをとる前記時間差であるピーク時間差および前記２次元相関値が最大のピークをとる前記周波数差であるピーク周波数差を検出するピーク検出部と、
位置が未知である電波源が放射する電波が第１の衛星で中継されて第１の観測局で受信された第１の受信信号を前記第１の信号とし、前記電波源が放射する電波が前記第１の衛星とは異なる第２の衛星で中継されて第２の観測局で受信された第２の受信信号を前記第２の信号として計算された前記２次元相関値の前記ピーク時間差から、前記第１の衛星および前記第２の衛星の位置を使用して、前記電波源から前記第１の衛星までの電波の到達時間と、前記電波源から前記第２の衛星までの電波の到達時間の差である到達時間差を計算する到達時間差計算部と、
前記第１の受信信号を前記第１の信号とし、前記第２の受信信号を前記第２の信号として計算された前記２次元相関値の前記ピーク周波数差から、前記第１の衛星および前記第２の衛星の位置および速度ベクトルを使用して、前記電波源から前記第１の衛星までで発生する周波数偏移と、前記電波源から前記第２の衛星までで発生する周波数偏移との差である到達周波数差を計算する到達周波数差計算部と、
前記第１の衛星および前記第２の衛星の位置を使用して、前記到達時間差が得られる前記電波源の位置を表す線である等到達時間差線を計算する等到達時間差線計算部と、
前記第１の衛星および前記第２の衛星が静止衛星かどうかを判別する静止衛星判別部と、
前記静止衛星判別部において、前記第１の衛星および前記第２の衛星の両方が静止衛星であると判別されると、前記到達周波数差から前記電波源の移動を表す速度ベクトルである電波源速度を推定する電波源速度推定部と、
前記静止衛星判別部において、前記第１の衛星および前記第２の衛星の両方が静止衛星であると判別されると、前記電波源速度と、前記等到達時間差線と、前記第１の衛星および前記第２の衛星の位置および速度ベクトルとを使用して、前記電波源の位置に関する電波源位置推定情報を生成する第２の測位部とを備えた位置標定装置。
【請求項2】
前記電波源速度推定部は、前記電波源の位置が存在すると想定する範囲である想定位置範囲に含まれる前記等到達時間差線上の位置で、前記電波源速度が存在すると想定する範囲である想定速度範囲に存在する前記電波源速度を求める、請求項１に記載の位置標定装置。
【請求項3】
前記想定位置範囲を記憶する想定位置範囲記憶部と、前記想定速度範囲を記憶する想定速度範囲記憶部と、
前記想定位置範囲および前記想定速度範囲をユーザが入力する入力手段とをさらに、備えた請求項２に記載の位置標定装置。
【請求項4】
前記電波源速度推定部が、前記第１の衛星および前記第２の衛星の位置および速度ベクトルを使用して、前記等到達時間差線の上に設定された想定位置での、前記第１の衛星が移動することで発生する周波数偏移と前記第２の衛星が移動することで発生する周波数偏移との差である衛星周波数差と、前記第１の衛星までで発生する前記電波源速度による周波数偏移と前記第２の衛星までで発生する前記電波源速度による周波数偏移との差である電波源速度周波数差との和と、前記到達周波数差との差が決められた範囲になるように前記電波源速度を推定し、前記第２の測位部が、前記電波源速度と、前記第１の衛星および前記第２の衛星の位置および速度ベクトルを使用して、前記想定位置での前記衛星周波数差を前記到達周波数差から前記電波源速度周波数差を減算した値に決め、前記衛星周波数差が得られる前記電波源の位置である推定等周波数差線を計算する推定等周波数差線計算部を含み、
前記電波源位置推定情報が、前記等到達時間差線および前記推定等周波数差線である、請求項１から請求項３の何れか１項に記載の位置標定装置。
【請求項5】
前記電波源が存在すると想定する位置である想定位置に対して、前記第１の衛星および前記第２の衛星の位置および速度ベクトルを使用して、前記第１の衛星が移動することで発生する周波数偏移と、前記第２の衛星が移動することで発生する周波数偏移との差である衛星周波数差を計算する衛星周波数差計算部をさらに備え、前記電波源速度推定部が、前記第１の衛星および前記第２の衛星の位置および速度ベクトルを使用して、前記第１の衛星までで発生する前記電波源速度による周波数偏移と前記第２の衛星までで発生する前記電波源速度による周波数偏移との差である電波源周波数差が、前記想定位置範囲に含まれる前記等到達時間差線の上に設定した前記想定位置に対する前記衛星周波数差を前記到達周波数差から減算した周波数差になるように、前記想定速度範囲に含まれる前記電波源速度を前記想定位置に応じて推定し、
前記第２の測位部が、前記等到達時間差線および推定した前記電波源速度を使用して算出した推定等周波数差線を前記電波源位置推定情報とする、請求項２または請求項３に記載の位置標定装置。
【請求項6】
前記電波源速度推定部が推定する前記電波源速度の範囲の前記想定速度範囲に対する割合である可能速度割合が決められた下限値以上である前記想定位置である可能位置を、前記電波源位置推定情報とする、請求項５に記載の位置標定装置。
【請求項7】
前記第２の測位部が、前記想定位置範囲に含まれる前記等到達時間差線の上に設定された各前記想定位置で推定される前記電波源速度の前記可能速度割合に基づき、前記可能位置よりも前記電波源が存在する可能性が高いと推定される最尤位置を求める最尤位置推定部を含み、前記電波源位置推定情報が、前記最尤位置、または前記最尤位置および前記可能位置である、請求項６に記載の位置標定装置。
【請求項8】
前記最尤位置推定部が、各前記可能位置で推定される前記電波源速度の前記可能速度割合をすべての各前記可能位置で推定される前記電波源速度の前記可能速度割合の和である全体可能速度割合で除算した値を重み付け係数として、各前記可能位置の重み付け平均した位置に基づき前記最尤位置を求める、請求項７に記載の位置標定装置。
【請求項9】
前記最尤位置推定部が、前記可能速度割合が最大になる前記可能位置を前記最尤位置として求める、請求項７に記載の位置標定装置。
【請求項10】
前記第２の測位部が、前記第１の衛星および前記第２の衛星の位置および速度ベクトルを使用して、前記衛星周波数差が得られる前記電波源の位置である推定等周波数差線を計算する推定等周波数差線計算部を含み、前記電波源位置推定情報が、前記推定等周波数差線を含む、請求項５から請求項９の何れか１項に記載の位置標定装置。
【請求項11】
前記到達時間差計算部が、位置が既知の参照局が放射する電波が前記第１の衛星で中継されて前記第１の観測局で受信された第３の受信信号を前記第１の信号とし、前記参照局が放射する電波が前記第２の衛星で中継されて前記第２の観測局で受信された第４の受信信号を前記第２の信号として計算された前記２次元相関値の前記ピーク時間差である参照局ピーク時間差と、前記第１の受信信号を前記第１の信号とし、前記第２の受信信号を前記第２の信号として計算された前記２次元相関値の前記ピーク時間差である電波源ピーク時間差と、前記参照局から前記第１の衛星までの電波の到達時間と前記参照局から前記第２の衛星までの電波の到達時間との差である参照局時間差とに基づき、前記到達時間差を計算し、前記到達周波数差計算部が、前記第３の受信信号を前記第１の信号とし、前記第４の受信信号を前記第２の信号として計算された前記２次元相関値の前記ピーク周波数差である参照局ピーク周波数差と、前記第１の受信信号を前記第１の信号とし、前記第２の受信信号を前記第２の信号として計算された前記２次元相関値の前記ピーク周波数差である電波源ピーク周波数差と、前記参照局から前記第１の衛星までで発生する周波数偏移と、前記参照局から前記第２の衛星までで発生する周波数偏移との差である参照局周波数差とに基づき、前記到達周波数差を計算する、請求項１から請求項１０の何れか１項に記載の位置標定装置。
【請求項12】
前記到達時間差計算部が、前記第１の観測局、前記第２の観測局、前記第１の衛星および前記第２の衛星の位置を使用して、前記到達時間差を計算し、前記到達周波数差計算部が、前記第１の観測局および前記第２の観測局の位置、前記第１の衛星および前記第２の衛星の位置および速度ベクトルを使用して、前記到達周波数差を計算する、請求項１から請求項１０の何れか１項に記載の位置標定装置。
【請求項13】
前記静止衛星判別部において、前記第１の衛星および前記第２の衛星の少なくとも１つが静止衛星ではないと判別されると、前記第１の衛星および前記第２の衛星の位置および速度ベクトルを使用して、前記等到達時間差線の上で前記到達周波数差が得られる位置を、前記電波源の位置として求める第１の測位部をさらに備える、請求項１から請求項１２の何れか１項に記載の位置標定装置。
【発明の詳細な説明】
【技術分野】
【0001】
本開示は、移動体に搭載されて衛星通信をする通信装置である電波源の位置を決める位置標定装置に関する。
【背景技術】
【0002】
衛星通信をする通信装置（電波源）の位置を決める（位置標定）ために、位置標定する対象の通信装置が異なる２個の衛星でそれぞれ中継されて地上局で受信される２個の受信信号についての到達時間差(TDOA)と到達周波数差(FDOA)を使用する方法が知られている（例えば、非特許文献１参照）。
【0003】
通信に使用される衛星は、静止衛星、準天頂軌道衛星、低軌道衛星などの種類がある。２個の衛星がどちらも静止衛星の場合には、２個の衛星でそれぞれ中継されて受信される２個の受信信号の周波数差の中で、電波源の移動により発生するドップラ効果による周波数偏移が、衛星の移動による周波数偏移よりも大きくなる場合がある。そのような場合には、従来の位置標定装置では、検出された周波数差が得られるような電波源の位置を表す等FDOA線を決めることができず位置標定できない場合がある、等FDOA線が求められる場合でも、等FDOA線の誤差が大きく、標定位置の精度が大きく劣化する場合がある。
【0004】
従来の位置標定装置は、位置標定できない場合には、電波源の位置を表示しない。位置標定できた場合は、精度が劣化している場合もそのまま表示する。
【先行技術文献】
【非特許文献】
【0005】
【文献】Ｄ． Ｐ． Ｈａｗｏｒｔｈ， ｅｔ．ａｌ，“Ｉｎｔｅｒｆｅｒｅｎｃｅ Ｌｏｃａｌｉｚａｔｉｏｎ Ｆｏｒ Ｅｕｔｅｌｓａｔ Ｓａｔｅｌｌｉｔｅｓ －Ｔｈｅ Ｆｉｒｓｔ Ｅｕｒｏｐｅａｎ Ｔｒａｎｓｍｉｔｔｅｒ Ｌｏｃａｔｉｏｎ Ｓｙｓｔｅｍ，” Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｓａｔｅｌｌｉｔｅ Ｃｏｍｍｕｎｉｃａｔｉｏｎｓ， Ｖｏｌ．１５， ｐｐ．１５５－１８３， １９９７
【発明の概要】
【発明が解決しようとする課題】
【0006】
従来の位置標定装置は、中継する衛星が２個とも静止衛星である場合に、電波源の位置に関する情報を適切に提示し難いという課題がある。
【課題を解決するための手段】
【0007】
本開示に係る位置標定装置は、複素数で表される第１の信号および第２の信号、第１の信号と第２の信号とのそれぞれ時間差および周波数差が入力されて、第１の信号に対して時間差を持たせた第２の信号の複素共役をとり、かつ周波数差で回転する絶対値が１である複素数を乗算して得られる信号と第１の信号とを乗算した値を決められた時間だけ積分して、第１の信号と第２の信号との間の２次元相関値を計算する相関関数計算部と、時間差と周波数差を変化させて２次元相関値を計算し、２次元相関値が最大のピークをとる時間差であるピーク時間差および２次元相関値が最大のピークをとる周波数差であるピーク周波数差を検出するピーク検出部と、位置が未知である電波源が放射する電波が第１の衛星で中継されて第１の観測局で受信された第１の受信信号を第１の信号とし、電波源が放射する電波が第１の衛星とは異なる第２の衛星で中継されて第２の観測局で受信された第２の受信信号を第２の信号として計算された２次元相関値のピーク時間差から、第１の衛星および第２の衛星の位置を使用して、電波源から第１の衛星までの電波の到達時間と、電波源から第２の衛星までの電波の到達時間の差である到達時間差を計算する到達時間差計算部と、第１の受信信号を第１の信号とし、第２の受信信号を第２の信号として計算された２次元相関値のピーク周波数差から、第１の衛星および第２の衛星の位置および速度ベクトルを使用して、電波源から第１の衛星までで発生する周波数偏移と、電波源から第２の衛星までで発生する周波数偏移との差である到達周波数差を計算する到達周波数差計算部と、第１の衛星および第２の衛星の位置を使用して、到達時間差が得られる電波源の位置を表す線である等到達時間差線を計算する等到達時間差線計算部と、第１の衛星および第２の衛星が静止衛星かどうかを判別する静止衛星判別部と、静止衛星判別部において、第１の衛星および第２の衛星の両方が静止衛星であると判別されると、到達周波数差から電波源の移動を表す速度ベクトルである電波源速度を推定する電波源速度推定部と、静止衛星判別部において、第１の衛星および第２の衛星の両方が静止衛星であると判別されると、電波源速度と、等到達時間差線と、第１の衛星および第２の衛星の位置および速度ベクトルとを使用して、電波源の位置に関する電波源位置推定情報を生成する第２の測位部とを備えたものである。
【発明の効果】
【0008】
本開示によれば、中継する衛星が２個とも静止衛星である場合に、電波源の位置に関する情報を従来よりも適切に提示できる。
【図面の簡単な説明】
【0009】
【図1】従来の位置標定装置で位置が未知の電波源の位置を標定する処理を説明する概念図である。
【図2】実施の形態１に係る位置標定装置で、位置が未知の電波源の位置を標定する処理を説明する概念図である。
【図3】実施の形態１に係る位置標定装置の構成を示すブロック図である。
【図4】本開示に係る位置標定装置で使用するＸＹＺ座標系を説明する模式図である。
【図5】実施の形態１に係る位置標定装置の動作を説明するフローチャートである。
【図6】実施の形態１に係る位置標定装置が出力する電波源位置推定情報の例を説明する図である。
【図7】実施の形態２に係る位置標定装置で、位置が未知の電波源の位置を標定する処理を説明する概念図である。
【図8】実施の形態２に係る位置標定装置の構成を示すブロック図である。
【図9】実施の形態２に係る位置標定装置の動作を説明するフローチャートである。
【図10】実施の形態３に係る位置標定装置の構成を示すブロック図である。
【図11】実施の形態３に係る位置標定装置の動作を説明するフローチャートである。
【図12】実施の形態３に係る位置標定装置が出力する電波源位置推定情報の１例を説明する図である。
【図13】実施の形態３に係る位置標定装置が出力する電波源位置推定情報の別の１例を説明する図である。
【図14】実施の形態３に係る位置標定装置において電波源位置推定情報を示す２個の例での存在確度を示すグラフである。
【図15】実施の形態４に係る位置標定装置の構成を示すブロック図である。
【図16】実施の形態４に係る位置標定装置が使用する可能方向範囲、想定速度範囲および許容方向範囲の関係を例により説明する図である。
【図17】実施の形態４に係る位置標定装置の動作を説明するフローチャートである。
【発明を実施するための形態】
【0010】
実施の形態１．
図１は、従来の位置標定装置で位置が未知の電波源の位置を標定する処理を説明する概念図である。図１は、非特許文献１に開示されている方法に基づいている。図１の例では、位置が未知である電波源７０から電波が送信されている。電波源７０が放射するメインローブの電波８１が衛星(＃１)７１に向けて送信される。同時に、電波源７０はサイドローブの電波８２を放射する。衛星(＃１)７１は電波８１を中継して電波８３を送信する。衛星(＃１)７１からの電波８３は、監視局のアンテナ(＃１)７２で受信される。アンテナ(＃１)７２は、電波８３を受信して受信信号９１を生成する。電波源７０が放射するサイドローブの電波８２は、衛星(＃２)７３により中継され、監視局のアンテナ(＃２)７４で電波８４として受信される。アンテナ(＃２)７４は、電波８４を受信して受信信号９２を生成する。
図１は、従来の位置標定装置で位置が未知の電波源の位置を標定する処理を説明する概念図である。図１は、非特許文献１に開示されている方法に基づいている。図１の例では、位置が未知である電波源７０から電波が送信されている。電波源７０が放射するメインローブの電波８１が衛星(＃１)７１に向けて送信される。同時に、電波源７０はサイドローブの電波８２を放射する。衛星(＃１)７１は電波８１を中継して電波８３を送信する。衛星(＃１)７１からの電波８３は、監視局のアンテナ(＃１)７２で受信される。アンテナ(＃１)７２は、電波８３を受信して受信信号９１を生成する。電波源７０が放射するサイドローブの電波８２は、衛星(＃２)７３により中継され、監視局のアンテナ(＃２)７４で電波８４として受信される。アンテナ(＃２)７４は、電波８４を受信して受信信号９２を生成する。
【0011】
アンテナ(＃１)７２とアンテナ(＃２)７４は、別の監視局に設置されてもよい。アンテナ(＃１)７２が設置される監視局が第１の監視局であり、アンテナ(＃２)７４が設置される監視局が第２の監視局である。受信信号９１が、第１の衛星である衛星(＃１)７１で中継されて第１の観測局で受信された第１の受信信号である。受信信号９２が、第２の衛星である衛星(＃２)７３で中継されて第２の観測局で受信された第２の受信信号である。アンテナ(＃１)７２とアンテナ(＃２)７４が同じ監視局に設置される場合には、その監視局が第１の監視局であり、かつ第２の監視局である。
【0012】
受信信号９１，９２が、位置標定装置５０Ｘに入力される。位置標定装置５０Ｘは、従来の方式で位置標定する位置標定装置である。位置標定装置５０Ｘは、衛星(＃１)７１と衛星(＃２)７３からそれぞれ中継されて受信された受信信号９１，９２の相関を取ることで、到達時間差(ＴＤＯＡ：Time Difference Of Arrival)と到達周波数差(ＦＤＯＡ：Frequency Difference Of Arrival)を推定する。推定したＴＤＯＡが一定である面と地球の表面との交線である等ＴＤＯＡ線６１と、推定したＦＤＯＡが一定である面と地球の表面との交線である等ＦＤＯＡ線６２の交点を、位置標定装置５０Ｘは電波源７０の位置と推定する。推定された電波源７０の位置すなわち等ＴＤＯＡ線６１と等ＦＤＯＡ線６２の交点を、標定位置６３と呼ぶ。
【0013】
位置が未知の電波源７０から放射された電波８１がアンテナ(＃１)７２で電波８３として受信されるまでの経路と、電波８２がアンテナ(＃２)７４で電波８４として受信されるまでの経路とにおいて、衛星(＃１)７１以降の経路と衛星(＃２)７３以降の経路は分かっている。受信信号の相関を取ることで得られる時間差と周波数差から、衛星(＃１)７１以降および衛星(＃２)７３以降の経路で発生する時間差および周波数差を減算して、到達時間差(ＴＤＯＡ)と到達周波数差(ＦＤＯＡ)を計算する。ＴＤＯＡは、電波８１が衛星(＃１)７１に電波が到達する時間と、電波８２が衛星(＃２)７３に到達する時間の差である。ＦＤＯＡは、電波８１が衛星(＃１)７１に受信されるまでに発生する周波数偏移と、電波８２が衛星(＃２)７３に受信されるまでに発生する周波数偏移との差である。
【0014】
まず、ＴＤＯＡおよびＦＤＯＡから電波源７０の位置を標定する従来の方法を説明する。
【0015】
以下の変数を定義する。
ｒ：電波源７０の位置。地球の中心を原点とする３次元座標系で表す。
ｒ１：衛星(＃１)７１の位置。
ｒ２：衛星(＃２)７３の位置。
ｒｍ１：アンテナ(＃１)７２の位置。
ｒｍ２：アンテナ(＃２)７４の位置。
ｃ：光速。
ｖ１：衛星(＃１)７１の速度ベクトル。
ｖ２：衛星(＃２)７３の速度ベクトル。
ｆ：電波源７０からの電波の搬送波の周波数。
ｆ１：衛星(＃１)７１でのアップリンク周波数とダウンリンク周波数の差。
ｆ２：衛星(＃２)７３でのアップリンク周波数とダウンリンク周波数の差。
Δｔ：衛星(＃１)７１、衛星(＃２)７３で中継される電波の時間差。
Δｆ：衛星(＃１)７１、衛星(＃２)７３で中継される電波の周波数差。
Δｔｍ：衛星(＃１)７１または衛星(＃２)７３からの経路で発生する時間差。
Δｆｍ：衛星(＃１)７１または衛星(＃２)７３からの経路で発生する周波数差。
Δｔr：衛星(＃１)７１または衛星(＃２)７３までに発生するＴＤＯＡ。
Δｆr：衛星(＃１)７１または衛星(＃２)７３までに発生するＦＤＯＡ。
ｕ１：電波源７０から衛星(＃１)７１へ向かう単位ベクトル。
ｕ２：電波源７０から衛星(＃２)７３へ向かう単位ベクトル。
ｕ１ｍ：アンテナ(＃１)７２から衛星(＃１)７１へ向かう単位ベクトル。
ｕ２ｍ：アンテナ(＃２)７４から衛星(＃２)７３へ向かう単位ベクトル。
ｓ１(t)：アンテナ(＃１)７２で受信される信号９１。
ｓ２(t)：アンテナ(＃２)７４で受信される信号９２。
Ｒ：地球の半径。地球は真球と仮定する。
ｒ：電波源７０の位置。地球の中心を原点とする３次元座標系で表す。
ｒ１：衛星(＃１)７１の位置。
ｒ２：衛星(＃２)７３の位置。
ｒｍ１：アンテナ(＃１)７２の位置。
ｒｍ２：アンテナ(＃２)７４の位置。
ｃ：光速。
ｖ１：衛星(＃１)７１の速度ベクトル。
ｖ２：衛星(＃２)７３の速度ベクトル。
ｆ：電波源７０からの電波の搬送波の周波数。
ｆ１：衛星(＃１)７１でのアップリンク周波数とダウンリンク周波数の差。
ｆ２：衛星(＃２)７３でのアップリンク周波数とダウンリンク周波数の差。
Δｔ：衛星(＃１)７１、衛星(＃２)７３で中継される電波の時間差。
Δｆ：衛星(＃１)７１、衛星(＃２)７３で中継される電波の周波数差。
Δｔｍ：衛星(＃１)７１または衛星(＃２)７３からの経路で発生する時間差。
Δｆｍ：衛星(＃１)７１または衛星(＃２)７３からの経路で発生する周波数差。
Δｔr：衛星(＃１)７１または衛星(＃２)７３までに発生するＴＤＯＡ。
Δｆr：衛星(＃１)７１または衛星(＃２)７３までに発生するＦＤＯＡ。
ｕ１：電波源７０から衛星(＃１)７１へ向かう単位ベクトル。
ｕ２：電波源７０から衛星(＃２)７３へ向かう単位ベクトル。
ｕ１ｍ：アンテナ(＃１)７２から衛星(＃１)７１へ向かう単位ベクトル。
ｕ２ｍ：アンテナ(＃２)７４から衛星(＃２)７３へ向かう単位ベクトル。
ｓ１(t)：アンテナ(＃１)７２で受信される信号９１。
ｓ２(t)：アンテナ(＃２)７４で受信される信号９２。
Ｒ：地球の半径。地球は真球と仮定する。
【0016】
位置標定では、アンテナ(＃１)７２で受信される信号ｓ１(t)と、アンテナ(＃２)７４で受信される信号の共役ｓ２
*(t)の間の、式(1)で計算される相関関数(Cross ambiguity function :CAF)を計算する。信号ｓ１(t)が相関関数を計算する上での第１の信号であり、信号ｓ２(t)が第２の信号である。第１の受信信号を第１の信号とし、第２の受信信号を第２の信号として、相関関数CAFを計算する。
【0017】
【数1】
【0018】
式(1)で計算する相関関数ＣＡＦがピーク値を取るτが、時間差Δｔである。また、ピークを取るνが、周波数差Δｆである。Δｔは、信号ｓ１(t)と信号ｓ２(t)の経路の長さの差により発生するので、下に示す式で表せる。
【0019】
【数2】
【0020】
周波数差Δｆは、衛星(＃１)７１でアンテナ(＃１)７２で発生するドップラ周波数偏移と、衛星(＃２)７３およびアンテナ(＃２)７４で発生するドップラ周波数偏移の差なので、以下の式で計算する。従来は、電波源７０の速度を考慮しないで、周波数差を計算している。
【0021】
【数3】
【0022】
式(3)で使用する各単位ベクトルは、下に示す式で計算する。
【0023】
【数4】
【0024】
Δtの中で、電波８１が衛星(＃１)７１に到達した後の経路と、電波８２が衛星(＃２)７３に到達した後の経路との差による成分であるΔｔｍを、下に示す式(８)により計算する。式(９)に示すように、ΔｔからΔｔｍを減算して、Δｔｒを計算する。Δｔｒに関する電波源７０の位置ｒに対する条件式は、式(10)である。位置標定装置５０Ｘは、Δｔｒに対して、式(10)を満足する位置ｒを計算する。
【0025】
【数5】
【0026】
Δｆの中で、電波源７０の位置ｒおよび速度ベクトルｖのどちらにもよらないで発生する周波数差であるΔｆｍを、下に示す式(11)により計算する。式(12)に示すように、Δｆ
からΔｆｍを減算して、Δｆｒを計算する。Δｆｒに関する電波源の位置ｒに対する条件式は、式(13)である。位置標定装置５０Ｘは、Δｆｒに対して、式(13)を満足する位置ｒを計算する。
からΔｆｍを減算して、Δｆｒを計算する。Δｆｒに関する電波源の位置ｒに対する条件式は、式(13)である。位置標定装置５０Ｘは、Δｆｒに対して、式(13)を満足する位置ｒを計算する。
【0027】
【数6】
【0028】
電波源７０は地球の表面に存在するとして、電波源７０の位置ｒに関して以下の条件が成立するものとする。
【0029】
【数7】
【0030】
従来の位置標定装置５０Ｘは、式(10)および式(14)を満足する位置ｒを表す等ＴＤＯＡ線６１を求める。また、式(13)および式(14)を満足する位置ｒを表す等ＦＤＯＡ線６２を求める。等ＴＤＯＡ線６１と等ＦＤＯＡ線６２の交点を、電波源７０の標定位置６３とする。
【0031】
電波源７０が移動していない場合、あるいは衛星(＃１)７１および衛星(＃２)７３の両方または一方が静止衛星でない場合は、従来の位置標定装置５０Ｘでも電波源７０の位置を適切な精度で標定できる。電波源７０が移動していない場合は、式(13)で周波数差を正しく計算できる。衛星(＃１)７１および衛星(＃２)７３の両方または一方が静止衛星でない場合は、衛星の速度に対して電波源７０の速度が小さいので、電波源７０の速度を考慮しない式(13)でも周波数差が得られる電波源７０の位置ｒを適切な精度で計算できる。
【0032】
静止衛星でない衛星の移動速度は電波源７０の移動速度よりも５０倍以上大きいと考えられる。これは、静止衛星でない衛星は、数時間以下の周期で地球を周回しており、地球の表面に対して毎秒１km程度以上の速度で移動するからである。電波源７０は、速度の上限は秒速３０ｍ程度と考えられる。船舶の場合は時速３０ノット(時速５５．６ｋｍ、秒速１５．４ｍ)が速度の上限と言われている。自動車の場合は、例えば時速１００ｋｍ（秒速２７．８ｍ）程度が速度の上限と考えられる。
【0033】
静止衛星の地球の表面に対する移動速度は秒速１ｍ程度である。電波源７０がある程度速く移動する場合には、電波源７０の移動速度が静止衛星の移動速度よりも大きくなる。その結果、衛星(＃１)７１および衛星(＃２)７３の両方が静止衛星である場合には、周波数差は主に電波源７０が移動することにより発生する。そのため、衛星(＃１)７１および衛星(＃２)７３の両方が静止衛星であり、かつ電波源７０が移動している場合には、電波源７０が移動することにより発生する周波数差を考慮しない式(13)を充足する電波源７０の位置が求められない場合がある。式(13)を充足する電波源７０の位置が求められても、実際の電波源７０の位置から大きく離れた位置である場合がある。
【0034】
本開示に係る位置標定装置は、電波源７０が移動する速度ベクトルを推定して、衛星(＃１)７１および衛星(＃２)７３の両方が静止衛星であり、かつ電波源７０が移動している場合に、従来よりも適切な位置標定に関する情報をユーザに提示するものである。図２は、実施の形態１に係る位置標定装置で、位置が未知の電波源の位置を標定する処理を説明する概念図である。従来との最も大きな違いは、移動する電波源７０の速度ベクトルｖを推定することである。
【0035】
本開示に係る位置標定装置５０では、周波数差Δｆに電波源７０の速度ベクトルｖに起因して発生する周波数差も考慮する。位置標定装置５０の処理を説明するために、以下の変数を定義する。
ｖ：電波源７０の速度ベクトル。電波源速度とも呼ぶ。
еFDOA(ｖ,ｒ)：電波源７０の速度ベクトルｖにより発生する周波数差。電波源周波数差と呼ぶ。
Δｆｃ(ｒ)：衛星(＃１)７１および衛星(＃２)７３で発生する周波数差。
ｕｄ：ｕ１とｕ２の差を表すベクトル。ｕｄ＝ｕ１‐ｕ２。
ｖ：電波源７０の速度ベクトル。電波源速度とも呼ぶ。
еFDOA(ｖ,ｒ)：電波源７０の速度ベクトルｖにより発生する周波数差。電波源周波数差と呼ぶ。
Δｆｃ(ｒ)：衛星(＃１)７１および衛星(＃２)７３で発生する周波数差。
ｕｄ：ｕ１とｕ２の差を表すベクトル。ｕｄ＝ｕ１‐ｕ２。
【0036】
位置標定装置５０では、周波数差ΔｆおよびΔｆｒに関して、下に示す式を使用する。周波数差Δｆは、式(15)示す成分を有する。式(15)によるΔｆ対しては、式(12)で計算されるΔfrに対する電波源７０の位置ｒおよび速度ベクトルｖの関係式は、式(16)になる。式(16)における第３項が、電波源７０が移動することを考慮することで発生する周波数差である。
【0037】
【数8】
【0038】
電波源７０の速度ベクトルｖが推定できると、電波源周波数差であるеFDOA(ｖ,ｒ)、を下に示す式(17)で計算する。ｕ１(r)，ｕ２(r)は、ｖを推定する際に使用した電波源７０の位置ｒでのベクトルを使用する。さらに、式(18)に示すように、ΔｆｒからеFDOA(ｖ,ｒ)を減算してΔｆｃ(ｒ)を計算する。Δｆｃ(ｒ)は、ＦＤＯＡであるΔｆｒから電波源周波数差еFDOA(ｖ,ｒ)を減算した周波数差である。Δｆｃ(ｒ)は、衛星が移動することで発生する周波数差の推定値である。Δｆｃ(ｒ)を推定ＦＤＯＡと呼ぶ。Δｆｃ(ｒ)に対して、式(19)で計算される電波源７０の位置ｒを結ぶ線を、推定等ＦＤＯＡ線６４と呼ぶ。
【0039】
【数9】
【0040】
図３を参照して、実施の形態１に係る位置標定装置の構成を説明する。位置標定装置５０は、周波数変換部１Ａ，１Ｂ、ＡＤ変換部２Ａ，２Ｂ、直交検波部３Ａ，３Ｂ、相関関数計算部４、ピーク検出部５、衛星データ取得部６、ＴＤＯＡ計算部７、ＦＤＯＡ計算部８、等ＴＤＯＡ線計算部９、静止衛星判別部１０、第１の測位部１１、電波源速度推定部１２、第２の測位部１３、記憶部１４を有する。第１の測位部１１は、等ＦＤＯＡ線計算部１５を有する。第２の測位部１３は、ＦＤＯＡ修正部１６、推定等ＦＤＯＡ線計算部１７を有する。
【0041】
位置標定装置５０には、受信信号９１，９２が入力される。周波数変換部１Ａは、受信信号９１を中間周波数に変換する。ＡＤ変換部２Ａは、周波数変換部１Ａが出力する信号を決められた時間の刻み幅で決められたビット数のデジタル値に変換する。直交検波部３Ａは、ＡＤ変換部２Ａでデジタル化された信号を直交検波して、複素ベースバンド信号ｓ１(t)を出力する。周波数変換部１Ｂ、ＡＤ変換部２Ｂおよび直交検波部３Ｂは、同様な処理を受信信号９２に対して実施する。直交検波部３Ｂは、複素ベースバンド信号ｓ２(t)を出力する。複素ベースバンド信号は、複素数で表される信号である。
【0042】
相関関数計算部４は、決められた期間Ｔのｓ１(t)およびｓ２(t)に対して、式(1)により相関関数CAFを計算する。ピーク検出部５は、τとνを変化させて相関関数CAFを計算し、相関関数CAFが最大のピークを取るτおよびνであるΔｔおよびΔｆを検出する。Δｔは、相関関数CAFが最大のピークをとるピーク時間差である。Δｆは、相関関数CAFが最大のピークをとるピーク周波数差である。信号９１，９２から、中継する衛星(＃１)７１および衛星(＃２)７３を決める。衛星データ取得部６は、中継する衛星(＃１)７１および衛星(＃２)７３の位置および速度ベクトルなどを含む衛星データを外部から取得する。
【0043】
ＴＤＯＡ計算部７は、位置ｒ１，ｒ２，ｒｍ１，ｒｍ２を使用して、式(9)によりΔｔからΔｔｒを計算する。Δｔｒが、到達時間差(ＴＤＯＡ)である。ＦＤＯＡ計算部８は、位置ｒ１，ｒ２，ｒｍ１，ｒｍ２と速度ベクトルｖ１，ｖ２とを使用して、式(12)によりΔｆからΔｆｒを計算する。Δｆｒが、到達周波数差(ＦＤＯＡ)である。なお、式(12)には位置ｒ１，ｒ２，ｒｍ１，ｒｍ２が含まれていないが、単位ベクトルｕ１ｍは、式(6)により位置ｒ１，ｒｍ１から計算され、単位ベクトルｕ２ｍは、式(7)により位置ｒ２，ｒｍ２から計算される。Δｆｒは、間接的に位置ｒ１，ｒ２，ｒｍ１，ｒｍ２を使用して計算される。
【0044】
等ＴＤＯＡ線計算部９は、位置ｒ１，ｒ２を使用して、式(10)および式(14)を満足する位置ｒを表す等ＴＤＯＡ線６１を計算する。等ＦＤＯＡ線計算部１５は、位置ｒ１，ｒ２と速度ベクトルｖ１，ｖ２とを使用して、式(13)および式(14)を満足する位置ｒを表す等ＦＤＯＡ線６２を計算する。静止衛星判別部１０は、衛星(＃１)７１が静止衛星かどうかを判別し、衛星(＃２)７３が静止衛星かどうかを判別する。第１の測位部１１は、衛星(＃１)７１および衛星(＃２)７３の少なくとも一方が静止衛星でない場合に、等ＴＤＯＡ線６１と等ＦＤＯＡ線６２の交点を電波源７０の標定位置６３（位置ｒ）として決定する。
【0045】
等ＴＤＯＡ線６１は、到達時間差(ＴＤＯＡ)が得られる電波源７０の位置ｒを表す線である等到達時間差線である。等ＴＤＯＡ線計算部９は、等到達時間差線を計算する等到達時間差線計算部である。等ＦＤＯＡ線６２は、到達周波数差(ＦＤＯＡ)が得られる電波源７０の位置ｒを表す線である等到達周波数差線である。等ＦＤＯＡ線計算部１５は、到達周波数差を計算する等到達周波数差線計算部である。
【0046】
電波源速度推定部１２は、衛星(＃１)７１および衛星(＃２)７３の両方が静止衛星である場合に、電波源７０の速度ベクトルｖを推定する。第２の測位部１３は、等ＴＤＯＡ線６１と速度ベクトルｖを使用して算出した推定等ＦＤＯＡ線６４を、電波源位置推定情報６５として生成する。
【0047】
第２の測位部１３は、ＦＤＯＡ修正部１６、推定等ＦＤＯＡ線計算部１７を有する。ＦＤＯＡ修正部１６は、推定された速度ベクトルｖからＦＤＯＡ(＝Δｆr)を修正してΔｆｃ(ｒ)を推定する。推定されたΔｆｃ(ｒ)の中から決められた個数のΔｆｃ(ｒ)をほぼ等間隔になるように抽出する。抽出したΔｆｃ(ｒ)に対して、推定等ＦＤＯＡ線計算部１７は、式(14)および式(19)を満足する電波源７０の位置ｒである推定等ＦＯＤＡ線６４を計算する。
【0048】
Δｆｃ(ｒ)は、衛星(＃１)７１が移動することで発生する周波数偏移と、衛星(＃２)７３が移動することで発生する周波数偏移との差である衛星周波数差である。衛星周波数差は、Δｆｒから電波源周波数差еFDOA(ｖ,ｒ)を減算した周波数差でもある。ＦＤＯＡ修正部１６は、衛星周波数差を計算する衛星周波数差計算部である。推定等ＦＤＯＡ線６４は、Δｆｃ(ｒ)が得られる電波源７０の位置ｒである推定等周波数差線である。推定等ＦＤＯＡ線計算部１７は、推定等周波数差線を計算する推定等周波数差線計算部である。
【0049】
記憶部１４は、位置ｒ１，ｒ２，ｒｍ１，ｒｍ２、速度ベクトルｖ１，ｖ２、周波数ｆ，ｆ１，ｆ２などの処理に必要なデータを記憶する。記憶部１４は、ｕ１，ｕ２，ｕ１ｍ，ｕ２ｍ，Δｔｍ，Δｆｍなどの処理の途中で使用するデータも記憶する。記憶部１４は、想定位置範囲３１、想定速度範囲３２も記憶する。想定位置範囲３１は、電波源７０の位置ｒが存在すると想定する範囲である。想定速度範囲３２は、電波源７０の速度ベクトルｖが存在すると想定する範囲である。電波源速度推定部１２は、想定位置範囲３１、想定速度範囲３２を使用して電波源７０の速度ベクトルｖを推定する。電波源速度推定部１２は、想定位置範囲３１に含まれる等ＴＤＯＡ線６１上の位置で、想定速度範囲３２に存在するｖを推定する。
【0050】
電波源速度推定部１２が電波源７０の速度ベクトルｖを推定する方法を説明する。以下の変数を定義する。
ｒTDOAk：想定位置範囲３１に含まれる等ＴＤＯＡ線６１の上に分散して配置した位置。ｎを配置した位置の数として、ｋは、ｋ＝１，２，…，ｎである。均等に配置することが望ましい。
Ｑk(ｒTDOAk, ｖ)：速度ベクトルｖを計算する際に使用する誤差関数。
ε：誤差関数Ｑk(ｒTDOAk, ｖ)に対する閾値。
ｖmax：想定速度範囲３２で規定する|ｖ|の上限値。
ｖmin：想定速度範囲３２で規定する|ｖ|の下限値。
(ｖｘ, ｖｙ, ｖｚ)：速度ベクトルｖを直交座標系で表現したベクトル。
Δｖ：速度ベクトルｖを探索する際の、直交座標系の各成分を変化させる刻み幅。
ｒTDOAk：想定位置範囲３１に含まれる等ＴＤＯＡ線６１の上に分散して配置した位置。ｎを配置した位置の数として、ｋは、ｋ＝１，２，…，ｎである。均等に配置することが望ましい。
Ｑk(ｒTDOAk, ｖ)：速度ベクトルｖを計算する際に使用する誤差関数。
ε：誤差関数Ｑk(ｒTDOAk, ｖ)に対する閾値。
ｖmax：想定速度範囲３２で規定する|ｖ|の上限値。
ｖmin：想定速度範囲３２で規定する|ｖ|の下限値。
(ｖｘ, ｖｙ, ｖｚ)：速度ベクトルｖを直交座標系で表現したベクトル。
Δｖ：速度ベクトルｖを探索する際の、直交座標系の各成分を変化させる刻み幅。
【0051】
想定位置範囲３１は、例えば、緯度が東経１３０度～１４０度、経度が２５度～４５度などのように設定する。想定位置範囲３１は、緯度と経度が同時に変化するような境界線を有するものでもよい。想定速度範囲３２は、例えば、スカラー値である|ｖ|に対して、上限値ｖmaxと下限値ｖminとにより設定する。上限値ｖmaxは、例えば秒速１５ｍに設定する。下限値ｖminは、例えば秒速０ｍに設定する。想定速度範囲３２は、速度ベクトルｖの方向も制限するものでもよい。
【0052】
地球中心を原点とし地球の自転とともに回転する直交座標系であるＸＹＺ座標系を使用する。図４を参照して、ＸＹＺ座標系を説明する。図４は、本開示に係る位置標定装置で使用するＸＹＺ座標系を説明する模式図である。ここで、Ｒcは静止軌道の半径(Ｒc≒35786km)である。ＸＹＺ座標系では、地球６０が自転する軸にＺ軸を設定する。北向きが、Ｚ軸の正の向きである。赤道７７を通る平面をＸＹ平面とする。衛星(＃１)７１および衛星(＃２)７３がＹＺ平面に対して対称に配置されるように、Ｙ軸を設定する。衛星(＃１)７１および衛星(＃２)７３が存在する側がＹ軸の正の向きとする。つまり、衛星(＃１)７１および衛星(＃２)７３の中間の経度がゼロ度になり、衛星(＃１)７１および衛星(＃２)７３の中間を通る経線７８が、ＹＺ平面上に存在する。
【0053】
Ｙ軸およびＺ軸に直交するように、Ｘ軸を設定する。Ｙ軸の正の向きに対して、東側に９０度の方向がＸ軸の正の向きである。なお、静止衛星である衛星(＃１)７１および衛星(＃２)７３が、静止軌道７６上に存在する場合で計算する。静止衛星が静止軌道から微小にずれた位置に存在する場合があるが、静止軌道からのずれは、ＸＹＺ座標系を回転して補正することで対応できる。
【0054】
ｒTDOAkは、想定位置範囲３１に含まれる等ＴＤＯＡ線６１の上に、例えば緯度が０．５度刻みで設定する。ｒTDOAkは、緯度の小数点以下が０または０．５の値を取るように設定してもよいし、０または０．５以外の値を取るように設定してもよい。想定位置範囲３１の広さに応じて、適切な間隔で適切な個数のｒTDOAkを設定する。ｒTDOAkは、想定位置範囲３１に含まれる等ＴＤＯＡ線６１の上に設定した想定位置である。
【0055】
各ｒTDOAkで電波源７０の速度ベクトルｖを求める際に使用する誤差関数Ｑk(ｒTDOAk, ｖ)は、下に示す式で定義する。
【0056】
【数10】
【0057】
電波源速度推定部１２は、速度ベクトルｖをＱk(ｒTDOAk, ｖ)が十分小さくなるように、すなわち下記の式を満足するように決める。εは、速度ベクトルｖが必要な精度で求められるように決める。
【0058】
【数11】
【0059】
式(20)で計算する誤差関数Ｑk(ｒTDOAk, ｖ)が式(21)を満足するように速度ベクトルｖを決定することは、ｖにより発生する周波数差であるеFDOA(ｖ,ｒTDOAk)が、下に示す式(22)を満足するようにｖを決めることを意味する。еFDOA(ｖ,ｒ)は、衛星(＃１)７１までで発生する電波源７０が移動する速度ベクトルｖによる周波数偏移と、衛星(＃２)７３までで発生する速度ベクトルｖによる周波数偏移との差である電波源周波数差である。電波源速度推定部１２は、速度ベクトルｖを想定位置ｒTDOAkに応じて推定する。電波源速度推定部１２は、еFDOA(ｖ,ｒTDOAk)がΔｆｒからΔｆｃ(ｒTDOAk)を減算した周波数差になるように、速度ベクトルｖを推定する。
【0060】
【数12】
【0061】
装置のＨ／Ｗ的な制約などのために、十分な計算時間、メモリ等が確保できない場合も考えられる。式(23)および式(24)に示すように、Ｑk(ｒTDOAk, ｖ)をＴＤＯＡ線に沿って平均を取った値が、十分小さくなるようなＶを計算してもよい。
【0062】
【数13】
【0063】
式(23)および式(24)を使用することで算出される速度ベクトルの数が削減され、それに伴い算出するFDOA線、標定点の数を減少し、演算負荷及び装置のメモリ等の負荷を低減することが可能となる。
【0064】
電波源速度推定部１２は、速度ベクトルｖの大きさ|ｖ|を以下の式で計算する。
【0065】
【数14】
【0066】
想定速度範囲３２は、下に示す式(26)のように、速度ベクトルｖの大きさ|ｖ|を制限する。
【0067】
【数15】
【0068】
電波源速度推定部１２は、ｖｘ，ｖｙ，ｖｚのそれぞれを、Δｖずつ変化させて、式(21)および式(26)を満足するｖをすべて計算する。Δｖの大きさは、必要な精度でｖが計算でき、かつ計算時間が過度に長くならないように決める。
【0069】
ＦＤＯＡ修正部１６は、速度ベクトルｖを推定する際に使用したｒTDOAkでのｕ1(ｒTDOAk)，ｕ２(ｒTDOAk)と推定された速度ベクトルｖとから式(17)によりеFDOA(ｖ,ｒ)を計算する。さらに、式(18)によりΔｆｃ(ｒ)を計算する。Δｆｃ(ｒ)が修正されたＦＤＯＡである。静止衛星は地球の中心から約35786ｋｍの距離を有する静止軌道にある。電波源７０の位置ｒが数100ｋｍ程度変化しても、ｕ1(ｒ)，ｕ２(ｒ)の変化は小さい。ｒが変化した際のｕ1(ｒ)とｕ２(ｒ)との差ベクトルｕ２(ｒ)‐ｕ１(ｒ)の変化は、ｕ1(ｒ)，ｕ２(ｒ)の変化よりも小さい。そのため、位置ベクトルｒが正確に求まっていない場合でも、еFDOA(ｖ,ｒ)を推定できる。
【0070】
推定されたΔｆｃ(ｒ)の中から、想定位置範囲内３１の大きさに応じた個数のΔｆｃ(ｒ)を抽出する。Δｆｃ(ｒ)の個数は、想定位置範囲内３１で適切な間隔で等ＦＤＯＡ線６４が配置されるように決める。Δｆｃ(ｒ)は、位置ｒが等間隔に近くなるように抽出する。抽出した各Δｆｃ(ｒ)に対して、推定等ＦＤＯＡ線計算部１７は、式(14)および式(19)を満足する電波源７０の位置ｒである推定等ＦＯＤＡ線６４を計算する。第２の測位部１３は、想定位置範囲内３１内の等ＴＤＯＡ線６１および推定等ＦＤＯＡ線６４を電波源位置推定情報として出力する。
【0071】
動作を説明する。図５は、実施の形態１に係る位置標定装置の動作を説明するフローチャートである。ステップＳ０１で、アンテナ(＃１)７２が衛星(＃１)７１で中継された電波８３を受信して第１の受信信号ｓ１(ｔ)を生成する。同時に、アンテナ(＃２)７４が衛星(＃２)７３で中継された電波８４を受信して第２の受信信号ｓ２(ｔ)を生成する。ステップＳ０２で、ピーク検出部５が、ｓ１(ｔ)とｓ２(ｔ)の相関関数が最大のピークになるピーク時間差Δｔおよびピーク周波数差Δｆを検出する。ステップＳ０３で、ＴＤＯＡ計算部７がＴＤＯＡであるΔｔｒを計算する。ＦＤＯＡ計算部８がＦＤＯＡであるΔｆｒを計算する。ステップＳ０４で、等ＴＤＯＡ線計算部９がΔｔｒに対して等ＴＤＯＡ線６１を計算する。
【0072】
ステップＳ０５で、静止衛星判別部１０が衛星(＃１)７１および衛星(＃２)７３が静止衛星かどうか判別する。ステップＳ０６で、衛星(＃１)７１および衛星(＃２)７３の両方が静止衛星かどうかチェックする。Ｓ０６で、衛星(＃１)７１および衛星(＃２)７３の少なくと一方が静止衛星でない（Ｓ０６でＮＯ）場合は、ステップＳ０７で、等ＦＤＯＡ線計算部１５がΔｆｒに対して等ＦＤＯＡ線６２を計算する。ステップＳ０８で、第１の測位部１１が、等ＴＤＯＡ線６１と等ＦＤＯＡ線６２の交点を電波源７０の標定位置６３(位置ｒ)として決定する。
【0073】
Ｓ０６で、衛星(＃１)７１および衛星(＃２)７３の両方が静止衛星である（Ｓ０６でＹＥＳ）場合は、ステップＳ０９で、電波源速度推定部１２が電波源７０の速度ベクトルｖを推定する。ステップＳ１０で、ＦＤＯＡ修正部１６が速度ベクトルｖを考慮してΔｆｒを修正して、Δｆｃ(ｒ)を計算する。ステップＳ１１で、推定等ＦＤＯＡ線計算部１７がΔｆｃ(ｒ)に対して推定等ＦＤＯＡ線６４を計算する。ステップＳ１２で、第２の測位部１３は、等ＴＤＯＡ線６１および想定位置範囲内３１に含まれる推定等ＦＤＯＡ線６４を電波源位置推定情報として出力する。
【0074】
図６に、位置標定装置５０が位置標定した結果の例を示す。
位置標定する条件は、以下である。
衛星(＃１)７１の位置：経度が２度の静止軌道。
衛星(＃２)７３の位置：経度が－２度の静止軌道。
衛星(＃１)７１の速度ベクトルｖ１：
ｖ１＝(ｖx,ｖy,ｖz)＝(０．５，－０．３，－０．９)[ｍ/秒]。
衛星(＃２)７３の速度ベクトルｖ２：
ｖ２＝(ｖx,ｖy,ｖz)＝(－０．３，－０．２，０．９)[ｍ/秒]。
Δｔr＝－０．２[ミリ秒]。
Δｆr＝２４[Hz]。
想定位置範囲３１：緯度２０～３５度、経度－１４～－１度。
ｖmax：１５[ｍ/秒]
位置標定する条件は、以下である。
衛星(＃１)７１の位置：経度が２度の静止軌道。
衛星(＃２)７３の位置：経度が－２度の静止軌道。
衛星(＃１)７１の速度ベクトルｖ１：
ｖ１＝(ｖx,ｖy,ｖz)＝(０．５，－０．３，－０．９)[ｍ/秒]。
衛星(＃２)７３の速度ベクトルｖ２：
ｖ２＝(ｖx,ｖy,ｖz)＝(－０．３，－０．２，０．９)[ｍ/秒]。
Δｔr＝－０．２[ミリ秒]。
Δｆr＝２４[Hz]。
想定位置範囲３１：緯度２０～３５度、経度－１４～－１度。
ｖmax：１５[ｍ/秒]
【0075】
想定位置範囲３１内での推定等ＦＤＯＡ線６４と、等ＴＤＯＡ線６１を電波源位置推定情報６５として表示する。推定等ＦＤＯＡ線６４には、衛星(＃１)７１および衛星(＃２)７３の移動により発生するＦＤＯＡであるΔｆc(ｒ)の値が分かるように表示する。推定等ＦＤＯＡ線６４は、図を見やすくするために４本だけ表示している。実際には、適切な間隔で適切な個数の推定等ＦＤＯＡ線６４を表示する。等ＴＤＯＡ線６１を、想定位置範囲３１内でだけ表示してもよい。想定位置範囲３１は、枠線ではなく、透過するハッチングで表示してもよい。想定位置範囲３１内でだけ推定等ＦＤＯＡ線６４を表示する場合は、想定位置範囲３１を表示しなくてもよい。想定位置範囲３１の内部および外部で推定等ＦＤＯＡ線６４を表示し、想定位置範囲３１も表示してもよい。
【0076】
中継する衛星が２個とも静止衛星であり、位置標定できない場合でも、位置標定装置５０は、電波源７０の速度ベクトルｖを推定して、速度ベクトルｖに基づく電波源位置推定情報６５を提示できる。中継する衛星が２個とも静止衛星である場合に、電波源７０の位置に関する情報を従来よりも適切に提示できる。
【0077】
実際の測位では、監視局の受信機材や通信衛星の発振器の周波数ずれや、通信衛星の軌道（位置、速度）にオフセットずれ（定常的な偏差）が発生することは避けられない。そのため、受信信号ｓ１（ｔ），ｓ２（ｔ）には時間あるいは周波数のオフセットずれが生じる。時間あるいは周波数のオフセットずれが発生していると、位置標定の精度が悪くなる。
【0078】
位置標定装置は、等ＦＤＯＡ線計算部１５を備えず、第１の測位部１１が、位置ｒ１，ｒ２と速度ベクトルｖ１，ｖ２とを使用して、等ＴＯＤＡ線６１上で式(13)を満足する位置を求めて、標定位置６３としてもよい。式(13)を満足する位置は、ＦＤＯＡを満足する位置である。等ＦＤＯＡ線計算部１５を有するかどうかによらず、第１の測位部１１が、等ＴＯＤＡ線６１上でＦＤＯＡを満足する位置を求めて標定位置６３とする。
以上のことは、他の実施の形態でもあてはまる。
以上のことは、他の実施の形態でもあてはまる。
【0079】
実施の形態２．
実施の形態２に係る位置標定装置５０Ａは、位置が既知な参照局であるリファレンス局が放射する電波が２個の衛星でそれぞれ中継されて受信される第３の受信信号および第４の受信信号も使用して、電波源の位置を標定するように実施の形態１を変更した場合である。実施の形態２に係る位置標定装置では、第１の受信信号および第２の受信信号で検出するピーク時間差Δｔおよびピーク周波数差Δｆから、第３の受信信号および第４の受信信号検出するΔｔrefおよびΔｆrefを減算することで、第１の受信信号および第２の受信信号に存在する時間あるいは周波数のオフセットずれを除去する。
実施の形態２に係る位置標定装置５０Ａは、位置が既知な参照局であるリファレンス局が放射する電波が２個の衛星でそれぞれ中継されて受信される第３の受信信号および第４の受信信号も使用して、電波源の位置を標定するように実施の形態１を変更した場合である。実施の形態２に係る位置標定装置では、第１の受信信号および第２の受信信号で検出するピーク時間差Δｔおよびピーク周波数差Δｆから、第３の受信信号および第４の受信信号検出するΔｔrefおよびΔｆrefを減算することで、第１の受信信号および第２の受信信号に存在する時間あるいは周波数のオフセットずれを除去する。
【0080】
図７は、実施の形態２に係る位置標定装置で、位置が未知の電波源の位置を標定する処理を説明する概念図である。図７では、位置が既知であるリファレンス局７５からも電波が放射される。リファレンス局７５が放射するメインローブの電波８５が衛星(＃１)７１に向けて送信される。同時に、リファレンス局７５はサイドローブの電波８６を放射する。衛星(＃１)７１は電波８５を中継して電波８７を送信する。衛星(＃１)７１からの電波８７は、監視局のアンテナ（＃１）７２で受信される。アンテナ（＃１）７２は、電波８７を受信して受信信号９３を生成する。受信信号９３が第３の受信信号である。リファレンス局７５が放射するサイドローブの電波８６は、衛星(＃２)７３により中継され、監視局のアンテナ（＃２）７４で電波８８として受信される。アンテナ（＃２）７４は、電波８８を受信して受信信号９４を生成する。受信信号９４が第４の受信信号である。
【0081】
リファレンス局７５を使用する場合について、ＴＤＯＡおよび周波数差(ＦＤＯＡ)から電波源７０の位置を標定する方法を説明する。
【0082】
以下の変数を定義する。
ｒｃ：リファレンス局７５の位置。
Δｔref：リファレンス局７５が放射する電波が衛星(＃１)７１、衛星(＃２)７３で中継されることで発生する時間差。
Δｆref：リファレンス局７５が放射する電波が衛星(＃１)７１、衛星(＃２)７３で中継されることで発生する周波数差。
Δｔtar-ref：ΔｔとΔｔrefの差。
Δｆtar-ref：ΔｆとΔｆrefの差。
Δｔｓ：リファレンス局７５が放射する電波で、衛星(＃１)７１および衛星(＃２)７３までに発生する時間差。
Δｆｓ：リファレンス局７５が放射する電波で、衛星(＃１)７１および衛星(＃２)７３までに発生する周波数差。
ｕ１ｃ：リファレンス局７５から衛星(＃１)７１へ向かう単位ベクトル。
ｕ２ｃ：リファレンス局７５から衛星(＃２)７３へ向かう単位ベクトル。
ｓ３(t)：アンテナ(＃１)７２で受信される信号９３。
ｓ４(t)：アンテナ(＃２)７４で受信される信号９４。
ｒｃ：リファレンス局７５の位置。
Δｔref：リファレンス局７５が放射する電波が衛星(＃１)７１、衛星(＃２)７３で中継されることで発生する時間差。
Δｆref：リファレンス局７５が放射する電波が衛星(＃１)７１、衛星(＃２)７３で中継されることで発生する周波数差。
Δｔtar-ref：ΔｔとΔｔrefの差。
Δｆtar-ref：ΔｆとΔｆrefの差。
Δｔｓ：リファレンス局７５が放射する電波で、衛星(＃１)７１および衛星(＃２)７３までに発生する時間差。
Δｆｓ：リファレンス局７５が放射する電波で、衛星(＃１)７１および衛星(＃２)７３までに発生する周波数差。
ｕ１ｃ：リファレンス局７５から衛星(＃１)７１へ向かう単位ベクトル。
ｕ２ｃ：リファレンス局７５から衛星(＃２)７３へ向かう単位ベクトル。
ｓ３(t)：アンテナ(＃１)７２で受信される信号９３。
ｓ４(t)：アンテナ(＃２)７４で受信される信号９４。
【0083】
リファレンス局７５から放射されてアンテナ(＃１)７２で受信される信号ｓ３(t)と、アンテナ(＃２)７４で受信される信号の共役ｓ４
*(t)の間でも、式(1)と同様に計算される相関関数CAFを計算する。第３の受信信号ｓ３(t)を第１の信号とし、第４の受信信号ｓ４(t)を第２の信号として、相関関数CAFを計算する。
【0084】
リファレンス局７５から放射される電波の受信信号に関して、式(1)で計算する相関関数がピーク値を取るτが、時間差Δｔrefである。また、ピークを取るνが、周波数差Δｆrefである。Δｔrefは、信号ｓ３(t)と信号ｓ４(t)の経路の長さの差により発生するので、下に示す式で表せる。また、リファレンス局７５から衛星までで発生する時間差Δｔsは、下に示す式で表せる。
【0085】
【数16】
【0086】
Δｔsは、リファレンス局７５から衛星(＃１)７１までの電波の到達時間と、リファレンス局７５から衛星(＃３)７３までの電波の到達時間との差である参照局時間差である。
【0087】
周波数差Δｆrefおよびリファレンス局７５から衛星までで発生する周波数差Δｆsは、以下の式で計算する。
【0088】
【数17】
【0089】
Δｆｓは、リファレンス局７５から衛星(＃１)７１までで発生する周波数偏移と、リファレンス局７５から衛星(＃２)７３までで発生する周波数偏移との差である参照局周波数差である。
【0090】
式(29)および式(30)で使用する単位ベクトルｕ１ｃ(r)，ｕ２ｃ(r)は、下に示す式で計算する。
【0091】
【数18】
【0092】
実施の形態２では、電波源７０から放射される電波が受信されて生成された受信信号ｓ１(t)およびｓ２(t)のピーク時間差Δｔと、リファレンス局７５から放射される電波が受信されて生成された受信信号ｓ３(t)およびｓ４(t)のピーク時間差Δｔrefとの差であるΔｔtar-refを計算する。また、受信信号ｓ１(t)およびｓ２(t)のピーク周波数差Δｆと、受信信号ｓ３(t)およびｓ４(t)のピーク周波数差Δｔrefとの差であるΔｆtar-refを計算する。ピーク時間差Δｔを電波源ピーク時間差と呼び、ピーク周波数差Δｆを電波源ピーク周波数差と呼ぶ。ピーク時間差Δｔrefを参照局ピーク時間差と呼び、ピーク周波数差Δｆrefを参照局ピーク周波数差と呼ぶ。
【0093】
Δｔtar-refは、下に示す式(33)に示すような成分を有する。
【0094】
【数19】
【0095】
Δｆtar-refは、下に示す式(34)に示すような成分を有する。
【0096】
【数20】
【0097】
位置が既知であるリファレンス局７５、衛星(＃１)７１および衛星(＃２)７３の位置ｒｃ，ｒ１，ｒ２を使用して、Δｔｓを式(28)により計算する。Δｔｓは、リファレンス局７５から衛星(＃１)７１および衛星(＃２)７３までの時間差である。さらに式(35)に示すように、Δｔtar-refからΔｔｓを減算したものをΔｔｒとする。
【0098】
【数21】
【0099】
位置が既知であるリファレンス局７５、衛星(＃１)７１および衛星(＃２)７３の位置ｒｃ，ｒ１，ｒ２と、衛星(＃１)７１および衛星(＃２)７３の速度ベクトルｖ１，ｖ２とを使用して、Δｆｓを式(30)により計算する。Δｆｓは、リファレンス局７５から衛星(＃１)７１および衛星(＃２)７３までで発生する周波数差である。さらに式(36)に示すように、Δftar-refからΔfｓを減算した周波数差をΔfｒとする。
【0100】
【数22】
【0101】
式(35)で計算するΔｔｒと式(36)で計算するΔｆｒを使用する実施の形態２では、衛星(＃１)７１からアンテナ(＃１)７２までの経路と、衛星(＃２)７３からアンテナ(＃２)７４までの経路とにおいて時間差または周波数差を発生させる要因およびその他のずれを発生させる要因の影響を排除して、電波源７０の位置ｒを推定できる。
【0102】
図８を参照して、実施の形態２に係る位置標定装置の構成を説明する。位置標定装置５０Ａには、受信信号９１，９２に加えて受信信号９３，９４が入力される。そのため、位置標定装置５０Ａは、周波数変換部１Ａ，１Ｂ，１Ｃ，１Ｄ、ＡＤ変換部２Ａ，２Ｂ，２Ｃ，２Ｄ、直交検波部３Ａ，３Ｂ，３Ｃ，３Ｄを有する。相関関数計算部４Ａ、ピーク検出部５Ａ、ＴＤＯＡ計算部７Ａ、ＦＤＯＡ計算部８Ａおよび記憶部１４Ａを変更している。
【0103】
周波数変換部１Ｃは、受信信号９３を中間周波数に変換する。ＡＤ変換部３Ｃは、周波数変換部１Ｃが出力する信号をデジタル値に変換する。直交検波部３Ｃは、ＡＤ変換部２Ｃでデジタル化された信号を直交検波して、複素ベースバンド信号ｓ３(t)を出力する。周波数変換部１Ｄ、ＡＤ変換部２Ｄおよび直交検波部３Ｄは、同様な処理を受信信号９４に対して実施する。直交検波部３Ｄは、複素ベースバンド信号ｓ４(t)を出力する。
【0104】
相関関数計算部４Ａは、決められた期間Ｔのｓ１(t)およびｓ２(t)に対して、式(1)により相関関数CAFを計算する。ピーク検出部５Ａは、相関関数CAFが最大のピークを取るΔｔおよびΔｆを検出する。Δｔは、ｓ１(t)およびｓ２(t)の間での相関関数CAFがピークになるピーク時間差である。Δｆは、ｓ１(t)およびｓ２(t)の間での相関関数CAFがピークになるピーク周波数差である。
【0105】
相関関数計算部４Ａおよびピーク検出部５Ａは、決められた期間Ｔのｓ３(t)およびｓ４(t)に対しても同様に動作する。相関関数計算部４Ａは、決められた期間Ｔのｓ３(t)およびｓ４(t)に対して、式(1)により相関関数CAFを計算する。ピーク検出部５Ａは、相関関数CAFが最大のピークを取るΔｔrefおよびΔｆrefを検出する。Δｔrefは、ｓ３(t)およびｓ４(t)の間での相関関数CAFがピークになるピーク時間差である。Δｆrefは、ｓ３(t)およびｓ４(t)の間での相関関数CAFがピークになるピーク周波数差である。
【0106】
ＴＤＯＡ計算部７Ａは、式(33)に示すようにΔｔからΔｔrefを減算してΔｔtar-refを計算する。位置ｒ１，ｒ２，ｒｃを使用して、式(28)に示すようにΔｔｓを計算する。さらに式(35)に示すようにΔｔtar-refからΔｔｓを減算して、Δｔｒを計算する。ＴＤＯＡ計算部７Ａは、ΔｔとΔｔrefとΔｔｓとに基づき、Δｔｒを計算する。
【0107】
ＦＤＯＡ計算部８Ａは、式(34)に示すようにΔｆからΔｆrefを減算してΔｆtar-refを計算する。位置ｒ１，ｒ２，ｒｃおよび速度ベクトルｖ１，ｖ２を使用して、式(30)に示すようにΔｆｓを計算する。さらに式(36)に示すようにΔｆtar-refからΔｆｓを減算して、Δｆｒを計算する。ＦＤＯＡ計算部８Ａは、ΔｆとΔｆrefとΔｆｓとに基づき、Δｆｒを計算する。
【0108】
記憶部１４Ａは、記憶部１４と比較して、処理に必要なデータと、処理の途中で使用するデータに関して記憶するデータを変更している。記憶部１４Ａは、位置ｒｍ１，ｒｍ２、Δｔｍ，Δｆｍ、ベクトルｕ１ｍ，ｕ２ｍなどを記憶しない。記憶部１４Ａは、位置ｒｃ、Δｔｃ，Δｆｃ、ベクトルｕ１ｃ，ｕ２ｃなどを記憶する。
【0109】
動作を説明する。図９は、実施の形態２に係る位置標定装置の動作を説明するフローチャートである。図９について、実施の形態１の場合の図５と異なる点を説明する。ステップＳ０１Ａで、アンテナ(＃１)７２が電波８３を受信して第１の受信信号ｓ１(ｔ)を生成し、電波８７を受信して第３の受信信号ｓ３(ｔ)を生成する。同時に、アンテナ(＃２)７４が電波８４を受信して第２の受信信号ｓ２(ｔ)を生成し、電波８８を受信して第４の受信信号ｓ４(ｔ)をする。ステップＳ０２Ａで、ピーク検出部５が、ｓ１(ｔ)およびｓ２(ｔ)に関してピーク時間差Δｔおよびピーク周波数差Δｆを検出し、ｓ３(ｔ)およびｓ４(ｔ)に関してピーク時間差Δｔrefおよびピーク周波数差Δｆrefを検出する。ステップＳ０３Ａで、ＴＤＯＡ計算部７Ａが、ΔｔとΔｔrefからＴＤＯＡであるΔｔｒを計算する。ＦＤＯＡ計算部８Ａが、ΔｆとΔｆrefからＦＤＯＡであるΔｆｒを計算する。
【0110】
ステップＳ０４以降の動作は、位置標定装置５０と同じである。
【0111】
リファレンス局が送信する電波を受信して得られる第３の受信信号および第４の受信信号も使用することで、第１の受信信号および第２の受信信号に存在する時間あるいは周波数のオフセットずれを除去する。リファレンス局を使用することで、リファレンス局を使用しない場合よりも、電波源を位置標定する精度が高くなる。
【0112】
実施の形態３．
実施の形態３は、実施の形態１を基に、等ＴＤＯＡ線上の位置を電波源位置推定情報として出力するようにした場合である。図１０は、実施の形態３に係る位置標定装置の構成を示すブロック図である。図１０について、実施の形態１の場合の図３とは異なる点を説明する。実施の形態２あるいは他の実施の形態を基に、同様な変更をしてもよい。
実施の形態３は、実施の形態１を基に、等ＴＤＯＡ線上の位置を電波源位置推定情報として出力するようにした場合である。図１０は、実施の形態３に係る位置標定装置の構成を示すブロック図である。図１０について、実施の形態１の場合の図３とは異なる点を説明する。実施の形態２あるいは他の実施の形態を基に、同様な変更をしてもよい。
【0113】
位置標定装置５０Ｂは、想定情報入力部１８を追加している。想定情報入力部１８は、ユーザが使用して想定位置範囲３１および想定速度範囲３２Ｂを入力する。想定速度範囲３２Ｂは、速度ベクトルｖの大きさ|ｖ|の範囲だけでなく、ｖの方向も制限する範囲である。想定情報入力部１８は、想定位置範囲３１および想定速度範囲３２Ｂをユーザが入力する入力手段である。
【0114】
位置標定装置５０Ｂでは、電波源速度推定部１２Ｂ、第２の測位部１３Ｂ、記憶部１４Ｂを変更している。第２の測位部１３Ｂは、ＦＤＯＡ修正部１６Ｂを変更している。電波源速度推定部１２Ｂは、想定速度範囲３２Ｂの範囲内で速度ベクトルｖを探索する。電波源速度推定部１２Ｂは、起動されるとすぐに想定速度範囲３２Ｂを読み込む。
【0115】
第２の測位部１３Ｂでは、ＦＤＯＡ修正部１６Ｂを変更している。ＦＤＯＡ修正部１６Ｂは、想定位置範囲３１に含まれる等ＴＤＯＡ線上に配置した位置である各ｒTDOAkで、推定ＦＤＯＡであるΔｆｃ(ｒTDOAk)を計算する。ＦＤＯＡ修正部１６Ｂは、各ｒTDOAkに対して、式(37)によりΔｆｃ(ｒTDOAk)を計算する。
【0116】
【数23】
【0117】
第２の測位部１３Ｂは、存在確度計算部１９および最尤位置推定部２０も有する。存在確度計算部１９は、各ｒTDOAkで計算されたすべての速度ベクトルｖを含む錐体を求める。錐体はすべての速度ベクトルｖを包含する方向の範囲を規定する。錐体で規定される速度ベクトルｖの方向の範囲を許容方向範囲６６と呼ぶ。電波源速度推定部１２Ｂは、推定速度範囲３２Ｂ内の速度ベクトルｖを求めているので、許容方向範囲６６は推定速度範囲３２Ｂに含まれる。存在確度計算部１９は、許容方向範囲６６が想定速度範囲３２Ｂに占める割合（存在確度）を計算する。さらに、存在確度が決められた下限値（存在判断下限値３３）以上である場合に、その点に電波源７０が存在する可能性があると判断する。電波源７０が存在する可能性があるｒTDOAkを、可能位置６７と呼ぶ。すべての可能位置６７について存在確度を重み付け係数として求めた重心に至近の等ＴＤＯＡ線６１上の位置を最尤位置６８と呼ぶ。最尤位置推定部２０が、最尤位置６８を求める。最尤位置６８は、可能位置６７よりも電波源７０が存在する可能性が高いと考えられる位置である。可能位置６７および最尤位置６８は、電波源位置推定情報として出力する。
【0118】
存在判断下限値３３は、例えば１０％に設定する。存在判断下限値３３は、想定位置範囲３１内における存在確度の最大値に対する比率で規定してもよい。存在判断下限値３３は、固定値および存在確度の最大値に対する比率の両方で、ＡＮＤ条件またはＯＲ条件で使用するような複数の値でもよい。ユーザが存在判断下限値３３を入力する手段を備えてもよい。
【0119】
推定等ＦＤＯＡ線計算部１７は、可能位置６７と判定された位置ｒTDOAkに対して、式(37)で計算される推定ＦＤＯＡであるΔｆc(ｒTDOAk)が得られる位置ｒである推定等ＦＯＤＡ線６４を計算する。
【0120】
記憶部１４Ｂは、想定位置範囲３１、想定速度範囲３２Ｂおよび在判断下限値３３を記憶する。記憶部１４Ｂは、処理の途中で使用するデータとして、記憶部１４とは異なるデータを記憶する。
【0121】
想定速度範囲３２Ｂでは、速度ベクトルｖの方向を指定するために、地平座標系での仰角と方位角を使用する。地平座標系をＸＹＺ座標系に変換する方法を説明する。
以下の変数を定義する。
φ：電波源７０の位置ｒの緯度。
θ：電波源７０の位置ｒの経度。
ξ：電波源７０の速度ベクトルｖの地平座標系での仰角。
ψ：電波源７０の速度ベクトルｖの地平座標系での方位角。南向きでψ＝０度。東向きでψ＝９０度。
ξmax：電波源７０の速度ベクトルｖの地平座標系での仰角の上限値。
ξmin：電波源７０の速度ベクトルｖの地平座標系での仰角の下限値。
ψmax：電波源７０の速度ベクトルｖの地平座標系での方位角の上限値。
ψmin：電波源７０の速度ベクトルｖの地平座標系での方位角の下限値。
[Ｍ]:地平座標系からＸＹＺ座標系への変換行列。
以下の変数を定義する。
φ：電波源７０の位置ｒの緯度。
θ：電波源７０の位置ｒの経度。
ξ：電波源７０の速度ベクトルｖの地平座標系での仰角。
ψ：電波源７０の速度ベクトルｖの地平座標系での方位角。南向きでψ＝０度。東向きでψ＝９０度。
ξmax：電波源７０の速度ベクトルｖの地平座標系での仰角の上限値。
ξmin：電波源７０の速度ベクトルｖの地平座標系での仰角の下限値。
ψmax：電波源７０の速度ベクトルｖの地平座標系での方位角の上限値。
ψmin：電波源７０の速度ベクトルｖの地平座標系での方位角の下限値。
[Ｍ]:地平座標系からＸＹＺ座標系への変換行列。
【0122】
想定速度範囲３２Ｂでは、下に示す式で速度ベクトルｖの範囲を制限する。
【0123】
【数24】
【0124】
地平座標系で極座標で表す速度ベクトルｖは、変換行列[Ｍ]により、下に示す式(41)のようにＸＹＺ直交座標系での値に変換できる。
【0125】
【数25】
【0126】
変換行列[Ｍ]は、下に示す式(42)で表される。
【0127】
【数26】
【0128】
電波源速度推定部１２Ｂは、起動されるとすぐに想定速度範囲３２Ｂを読み込み、想定速度範囲３２Ｂを式(41)によりＸＹＺ座標系での範囲に変換する。電波源速度推定部１２Ｂは、速度ベクトルｖを探索する際には、想定速度範囲３２Ｂの範囲内で探索する。
【0129】
直交座標系での速度ベクトルＶは、変換行列[Ｍ－１]により、下に示す式(43)で表されるように地平座標系で極座標で表す速度ベクトルに変換できる。
【0130】
【数27】
【0131】
変換行列[Ｍ－１]は、下に示す式(44)で表される。
【0132】
【数28】
【0133】
存在確度計算部１９の処理を説明するため、以下の変数を定義する。
Ｇ(ｒTDOAk)：位置ｒTDOAkでの許容方向範囲６６の面積。
Ｇdet：想定速度範囲３２Ｂの面積。
γ(ｒTDOAk)：Ｇ(ｒTDOAk)のＧdetに対する割合。存在確度と呼ぶ。
γmin：ｒTDOAkに電波源７０が存在すると判断する存在確度γに対する下限値。存在判断下限値３３として、記憶部１４Ｂに記憶される。
ｒcent：可能位置６７であるｒTDOAkの重み付け重心。
ｒmost：最尤位置６８を表す。ｒcentに至近のＴＯＤＡ線６１上の位置。
Ｇ(ｒTDOAk)：位置ｒTDOAkでの許容方向範囲６６の面積。
Ｇdet：想定速度範囲３２Ｂの面積。
γ(ｒTDOAk)：Ｇ(ｒTDOAk)のＧdetに対する割合。存在確度と呼ぶ。
γmin：ｒTDOAkに電波源７０が存在すると判断する存在確度γに対する下限値。存在判断下限値３３として、記憶部１４Ｂに記憶される。
ｒcent：可能位置６７であるｒTDOAkの重み付け重心。
ｒmost：最尤位置６８を表す。ｒcentに至近のＴＯＤＡ線６１上の位置。
【0134】
存在確度計算部１９は、各ｒTDOAkに対して、そのｒTDOAkで計算されたすべての速度ベクトルｖを、式(43)により地平座標系での速度ベクトルに変換する。地平座標系での速度ベクトルｖをすべて含む範囲を求める。求めた範囲が、許容方向範囲６６である。式(45)に示すように、許容方向範囲６６の面積であるＧ(ｒTDOAk)を推定速度範囲３２Ｂの面積で除算して、存在確度γ(ｒTDOAk)を計算する。
【0135】
【数29】
【0136】
ここで、推定速度範囲３２Ｂの面積は、仰角ξおよび方位角ψを直交するように設定したξψ平面において、推定速度範囲３２Ｂが占める範囲の面積である。許容方向範囲６６も、ξψ平面における面積である。
【0137】
可能位置６７は、存在確度γ(ｒTDOAk)がγmin以上（γ(ｒTDOAk)≧γmin）であるｒTDOAkである。最尤位置推定部２０が、すべての可能位置６７であるｒTDOAkを重み付けして重心ｒcentを求める。
【0138】
【数30】
【0139】
存在確度γ(ｒTDOAk)は、電波源速度推定部１２Ｂが推定する速度ベクトルｖの範囲の想定方向範囲３２Ｂに対する割合である可能速度割合である。式(46)の右辺の分母は、すべての可能位置６７で推定される速度ベクトルｖの存在確度γ(ｒTDOAk)の和である全体可能速度割合である。各可能位置６７の存在確度γ(ｒTDOAk)を全体可能速度割合で除算した値を重み付け係数として、各ｒTDOAkを重み付け平均した位置、すなわち重み付け重心ｒcentを求める。最尤位置６８は、重み付け重心ｒcentに基づき決める。
【0140】
動作を説明する。図１１は、実施の形態３に係る位置標定装置の動作を説明するフローチャートである。図１１について、実施の形態１の場合の図５と異なる点を説明する。ステップＳ０９Ｂで、電波源速度推定部１２Ｂが、速度ベクトルｖの方向も制限する想定速度範囲３２Ｂを考慮して電波源７０の速度ベクトルｖを推定する。Ｓ０９Ｂの次にステップＳ１３で、各ｒTDOAkに対して計算されたすべての速度ベクトルｖを含む許容方向範囲６６を求める。ステップＳ１４で、許容方向範囲６６から各ｒTDOAkの存在確度γ(ｒTDOAk)を計算する。ステップＳ１５で、γ(ｒTDOAk)≧γminであるｒTDOAkを、可能位置６７とする。ステップＳ１６で、可能位置６７の重み付け重心ｒcentを計算する。ステップＳ１７で、ｒcentに至近のＴＯＤＡ線６１上の位置を求めて、最尤位置６８であるｒmostとする。ステップＳ１１Ｂで、推定等ＦＤＯＡ線計算部１７が可能位置６７に対して推定等ＦＤＯＡ線６４を計算する。ステップＳ１２Ｂで、第２の測位部１３Ｂは、想定位置範囲内３１に含まれる等ＴＤＯＡ線６１、推定等ＦＤＯＡ線６４、可能位置６７および最尤位置６８を電波源位置推定情報として出力する。電波源位置推定情報として、可能位置６７を通る推定等ＦＤＯＡ線６４あるいは可能位置６７を出力しなくてもよい。
【0141】
ステップＳ１８で、想定位置範囲３１または想定速度範囲３２Ｂが変更されたかチェックする。想定位置範囲３１または想定速度範囲３２Ｂが変更されている（Ｓ１８でＹＥＳ）場合は、Ｓ０９Ｂへ戻る。想定位置範囲３１または想定速度範囲３２Ｂが変更されていない（Ｓ１８でＮＯ）場合は、処理を終了する。
【0142】
図１２および図１３に、位置標定装置５０Ｂが位置標定した結果の例を示す。図１２が速度ベクトルｖの大きさ|ｖ|とｖの仰角ξに制約条件を設定し、ｖの方位角ψを制約しない場合である。図１３が速度ベクトルｖの大きさ|ｖ|、仰角ξおよび方位角ψを制約する場合である。
【0143】
位置標定装置５０Ｂは、電波源位置推定情報６５Ｂとして、可能位置６７と最尤位置６８も生成する。図では、可能位置６７は白丸で表示し、最尤位置６８は×印で表示する。
【0144】
設定した制約条件は、以下である。図１２に示す場合をCase１とし、図１３に示す場合をCase２とする。
γmin＝０．１
(Case１)
|ｖ|≦１５[ｍ/秒]
－１２[度]≦ξ≦１２[度]
－１８０[度]≦ψ≦１８０[度]
(Case２)
１３[ｍ/秒]≦|ｖ|≦１５[ｍ/秒]
－１２[度]≦ξ≦１２[度]
１３５[度]≦ψ≦１４５[度]
γmin＝０．１
(Case１)
|ｖ|≦１５[ｍ/秒]
－１２[度]≦ξ≦１２[度]
－１８０[度]≦ψ≦１８０[度]
(Case２)
１３[ｍ/秒]≦|ｖ|≦１５[ｍ/秒]
－１２[度]≦ξ≦１２[度]
１３５[度]≦ψ≦１４５[度]
【0145】
想定位置範囲３１に含まれる等ＴＤＯＡ線６１の上の各位置ｒTDOAkでの存在確度γは、図１４に示すようになる。図１４では横軸を位置ｒTDOAkの緯度φとしている。Case１、Case２のどちらでも、想定位置範囲３１に含まれる各位置ｒTDOAkで、存在確度γ＞γminである。そのため、想定位置範囲３１に含まれる等ＴＤＯＡ線６１の上の各位置ｒTDOAkは可能位置６７である。各可能位置６７について、その存在確度γをユーザが参照できる。
図を見やすくするため、可能位置６７を４個、表示している。可能位置６７は、適切な間隔で適切な個数を表示すればよい。
図を見やすくするため、可能位置６７を４個、表示している。可能位置６７は、適切な間隔で適切な個数を表示すればよい。
【0146】
Case１とCase２では、Case２の方が想定方向範囲３２Ｂを狭く設定している。そのため、位置ｒTDOAkの中で緯度φが大きいものの存在確度γが大きくなる。その結果、Case２での最尤位置６８は、Case１での最尤位置６８よりも緯度が大きい位置になる。衛星(＃１)７１および衛星(＃２)７３の速度ベクトルｖ１，ｖ２が異なる状況で、速度ベクトルｖの大きさ|ｖ|と方向の制約をより狭く設定すると、想定位置範囲３１に含まれる等ＴＤＯＡ線６１の上の各位置ｒTDOAkで、可能位置６７でないものが発生する場合がある。また、設定した想定方向範囲３２Ｂに応じて、最尤位置６８が変化する。
【0147】
位置標定装置５０Ｂでは、中継する衛星が２個とも静止衛星であり、位置標定できない場合でも、電波源７０の速度ベクトルｖを推定して、速度ベクトルｖに基づく電波源位置推定情報６５Ｂを提示できる。中継する衛星が２個とも静止衛星である場合に、電波源７０の位置に関する情報を従来よりも適切に提示できる。
【0148】
最尤位置６８に対応する推定等ＦＤＯＡ線６４を表示してもよい。可能位置６７に対応させて推定等ＦＤＯＡ線６４を表示しているが、推定等ＦＤＯＡ線６４は表示しなくてもよい。
【0149】
各可能位置６７の中で存在確度γが最大になる可能位置６７を、最尤位置６８としてもよい。存在確度γが最大になる可能位置６７は、例えば２分探索法により求めればよい。
【0150】
存在確度γ(ｒTDOAk)としては、直交座標系での許容方向範囲６６および想定速度範囲３２Ｂの立体角の比として計算してもよい。あるいは、重み付け和を利用して存在確度γ(ｒTDOAk)を計算してもよい。重み付け和とは、重み付け関数ｇ(ξ, ψ)を使用してＧ(ｒTDOAk)を計算することである。重み付け関数ｇ(ξ,ψ)は、ある方向(ξ,ψ)が許容方向範囲６６に含まれない場合に０を返し、含まれる場合に０でない値を返す関数である。重み付け関数g(ξ,ψ)は、０でない値を返す場合に、ξまたはψの少なくとも一方の値に応じて返す値を変化させる。そうすることで、ξまたはψの上限または下限に近い範囲での重みｇ(ξ, ψ)を小さくしてＧ(ｒTDOAk)を計算することができる。なお、許容方向範囲６６のξψ平面に面積であるＧ(ｒTDOAk)を計算することは、方向(ξ,ψ)が許容方向範囲６６に含まれない場合に０を返し、許容方向範囲６６に含まれる場合に１を返す関数（仮にｈ(ξ, ψ)とする）を、想定速度範囲３２Ｂに含まれるξおよびψで積分することと等価である。
以上のことは、他の実施の形態でもあてはまる。
以上のことは、他の実施の形態でもあてはまる。
【0151】
実施の形態４．
実施の形態４は、実施の形態３を基に、電波源の移動を表す速度ベクトルｖの推定方法を変更した場合である。実施の形態４に係る位置標定装置５０Ｃでは、各ｒTDOAkに対してＦＤＯＡであるΔｆｒと修正ＦＤＯＡであるΔｆｃ(ｒTDOAk)の差に等しい、電波源周波数差であるеFDOA(ｖ,ｒTDOAk)を計算する。
実施の形態４は、実施の形態３を基に、電波源の移動を表す速度ベクトルｖの推定方法を変更した場合である。実施の形態４に係る位置標定装置５０Ｃでは、各ｒTDOAkに対してＦＤＯＡであるΔｆｒと修正ＦＤＯＡであるΔｆｃ(ｒTDOAk)の差に等しい、電波源周波数差であるеFDOA(ｖ,ｒTDOAk)を計算する。
【0152】
図１５は、実施の形態４に係る位置標定装置５０Ｃの構成を示すブロック図である。図１５について、実施の形態３の場合の図１０とは異なる点を説明する。位置標定装置５０Ｃでは、電波源速度推定部１２Ｃおよび存在確度計算部１９Ｃを変更している。
【0153】
電波源速度推定部１２Ｃは、各ｒTDOAkに対して式(47)で計算されるеFDOA(ｖ,ｒTDOAk)に対して、式(48)を満足するような電波源７０の速度ベクトルｖを計算する。なお、Δｆｃ(ｒTDOAk)は、式(37)で計算する。
【0154】
【数31】
【0155】
式(48)を満足する速度ベクトルｖを計算する方法を説明するために、以下の変数を定義する。
ｕｄ(ｒTDOAk)：ｕ１(ｒTDOAk)とｕ２(ｒTDOAk)の差ベクトル。
(ｕｘ, ｕｙ, ｕｚ):ｕｄ(ｒTDOAk)を直交座標で表したベクトル。
φd：ｕｄ(ｒTDOAk)を極座標で表現した際のＸＹ平面となす角度。
θd：ｕｄ(ｒTDOAk)を極座標で表現した際のＹＺ平面となす角度。
ζ：速度ベクトルｖとｕｄ(ｒTDOAk)とがなす角度。
ρ：速度ベクトルｖのｕｄ(ｒTDOAk)の回りの回転角度。
ζmax：еFDOA(ｖ,ｒTDOAk)が得られる角度ζの上限値。
ζmax：еFDOA(ｖ,ｒTDOAk)が得られる角度ζの下限値。
ｕｄ(ｒTDOAk)：ｕ１(ｒTDOAk)とｕ２(ｒTDOAk)の差ベクトル。
(ｕｘ, ｕｙ, ｕｚ):ｕｄ(ｒTDOAk)を直交座標で表したベクトル。
φd：ｕｄ(ｒTDOAk)を極座標で表現した際のＸＹ平面となす角度。
θd：ｕｄ(ｒTDOAk)を極座標で表現した際のＹＺ平面となす角度。
ζ：速度ベクトルｖとｕｄ(ｒTDOAk)とがなす角度。
ρ：速度ベクトルｖのｕｄ(ｒTDOAk)の回りの回転角度。
ζmax：еFDOA(ｖ,ｒTDOAk)が得られる角度ζの上限値。
ζmax：еFDOA(ｖ,ｒTDOAk)が得られる角度ζの下限値。
【0156】
ｕｄ(ｒTDOAk)は、下に示す式(49)で計算する。ｕ１(ｒTDOAk)およびｕ２(ｒTDOAk)は、式(4)および式(5)においてｒ＝ｒTDOAkとして計算する。
【0157】
【数32】
【0158】
ｕｄ(ｒTDOAk)を極座標で表現する際の角度φdおよびθdは、以下の式で計算する。
【0159】
【数33】
【0160】
式(49)を、式(48)に代入して、下に示す式(53)が得られる。式(53)を変形して、式(54)となる。
【0161】
【数34】
【0162】
式(54)をcosζについて解くと、式(55)となる。
【0163】
【数35】
【0164】
式(55)を満足する角度ζを、速度ベクトルｖの大きさ|ｖ|の制約条件である式(38)を考慮して求める。余弦関数は１以下の値しか取らないので、式(55)を満足する角度ζが存在するためには、|ｖ|の上限値であるｖmaxは、下に示す式(56)を満足する必要がある。ｖmaxが式(56)を満足しない場合は、式(55)を満足する速度ベクトルｖは存在しない。
【0165】
【数36】
【0166】
ｖmaxが式(56)を満足する場合には、角度ζの上限値ζmaxは、式(57)で計算する。
【0167】
【数37】
【0168】
ｖmaxが式(56)を満足する場合には、角度ζの下限値ζminは、ｖminの値を場合分けして下に示す式で計算する。
【0169】
【数38】
【0170】
電波源速度推定部１２Ｃは、各ｒTDOAkに対して、式(56)が成立しない場合は、速度ベクトルｖは存在しないとする。式(56)が成立する場合は、式(57)で角度ζの上限値ζmaxを計算し、式(58-1)または式(58-2)で角度ζの下限値ζminを計算する。ζmaxおよびζminを使用して、式(49)と式(38)を満足する速度ベクトルｖとして、以下の範囲が得られる。
【0171】
【数39】
【0172】
速度ベクトルｖは、ベクトルｕd(ｒTDOAk)となす角度ζがζmax以下かつζmin以上であり、ｕd(ｒTDOAk)に平行な成分が(c/f)(еFDOA(ｖ,ｒ)/|ｕd(ｒTDOAk)|)であるベクトルである。式(59)～(61)を満足するすべての速度ベクトルｖの始点をｕd(ｒTDOAk)の始点に一致させると、ｖの先端は、ｕd(ｒTDOAk)に垂直な平面で円またはドーナツ状に存在する。
【0173】
速度ベクトルｖを、下に示すように変換行列[Ｔ]を使用してＸＹＺ座標系での値に変換できる。
【0174】
【数40】
【0175】
変換行列[Ｔ]は、下に示すように表される。
【0176】
【数41】
【0177】
式(57)から式(60)で規定される速度ベクトルｖの方向は、式(49)と式(38)を満足する速度ベクトルｖの方向である。式(49)と式(38)を満足する速度ベクトルｖの方向の範囲を可能方向範囲６９と呼ぶ。
【0178】
存在確度計算部１９Ｃは、ＸＹＺ座標系で表現された可能方向範囲６９を式(43)により地平座標系での方向の範囲に変換する。角度ζおよびρで表現される方向を、直接に地平座標系での方向に変換してもよい。その場合には、下に示す式で変換できる。
【0179】
【数42】
【0180】
存在確度計算部１９Ｃは、図１６に示すように、可能方向範囲６９と想定速度範囲３２Ｂとの重複部分を許容方向範囲６６として求める。図１６は、可能方向範囲６９、想定速度範囲３２Ｂおよび許容方向範囲６６の関係を例により説明する図である。図１６では、方位角ψが正の範囲だけを表示している。方位角ψは、３６０[度](２π[rad])の範囲で変化する。可能方向範囲６９に含まれる方向であり、かつ想定速度範囲３２Ｂにも含まれる方向が、許容方向範囲６６である。図１６では、想定速度範囲３２Ｂを、|ξ|≦１２[度]、７５[度]≦ψ≦１５０[度]に設定している。環状の領域である可能方向範囲６９と想定方向範囲３２との両方に含まれる方向の範囲が、許容方向範囲６６である。図１６では、許容方向範囲６６にハッチングを付して表示している。
【0181】
動作を説明する。図１７は、実施の形態４に係る位置標定装置の動作を説明するフローチャートである。図１７について、実施の形態３の場合の図１１と異なる点を説明する。ステップＳ０９Ｃで、電波源速度推定部１２Ｃが、速度ベクトルｖの大きさ|ｖ|を考慮し、ｖの方向を考慮しないで電波源７０の速度ベクトルｖを推定する。推定したすべての速度ベクトルｖの方向を含む範囲を、可能方向範囲６９とする。ステップＳ１３Ｃで、可能方向範囲６９と想定速度範囲３２Ｂとが重なり合う部分を、存在確度計算部１９Ｃが許容方向範囲６６として求める。
その他の点では、位置標定装置５０Ｂと同様に動作する。
その他の点では、位置標定装置５０Ｂと同様に動作する。
【0182】
位置標定装置５０Ｃでは、中継する衛星が２個とも静止衛星であり、位置標定できない場合でも、電波源７０の速度ベクトルｖを推定して、速度ベクトルｖに基づく電波源位置推定情報６５Ｂを提示できる。中継する衛星が２個とも静止衛星である場合に、電波源７０の位置に関する情報を従来よりも適切に提示できる。
【0183】
位置標定装置５０Ｃでは、電波源７０の速度ベクトルｖの可能方向範囲６９を解析的に求めている。そのため、位置標定装置５０Ｃは、位置標定装置５０Ｂよりも可能方向範囲６９を効率的に求めることができる。
【0184】
各実施の形態の自由な組み合わせ、あるいは各実施の形態の変形や一部の構成要素を省略すること、あるいは一部の構成要素の省略や変形をした各実施の形態の自由な組み合わせが可能である。
【符号の説明】
【0185】
５０、５０Ａ、５０Ｂ、５０Ｃ、５０Ｘ 位置標定装置、
１Ａ，１Ｂ，１Ｃ，１Ｄ 周波数変換部、
２Ａ，２Ｂ，２Ｃ，３Ｄ ＡＤ変換部、
３Ａ，３Ｂ，３Ｃ，３Ｄ 直交検波部、
４、４Ａ 相関関数計算部、
５、５Ａ ピーク検出部、
６ 衛星データ取得部、
７、７Ａ ＴＤＯＡ計算部、
８、８Ａ ＦＤＯＡ計算部、
９ 等ＴＤＯＡ線計算部、
１０ 静止衛星判別部、
１１ 第１の測位部、
１２、１２Ｂ、１２Ｃ 電波源速度推定部、
１３、１３Ｂ、１３Ｃ 第２の測位部、
１４、１４Ａ、１４Ｂ、１４Ｃ 記憶部、
１５ 等ＦＤＯＡ線計算部、
１６、１６Ｂ ＦＤＯＡ修正部、
１７ 推定等ＦＤＯＡ線計算部、
１８ 想定情報入力部、
１９、１９Ｃ 存在確度計算部、
２０ 最尤位置推定部、
３１ 想定位置範囲、
３２、３２Ｂ 想定速度範囲、
３３ 存在判断下限値、
６０ 地球、
６１ 等ＴＤＯＡ線（等到達時間差線）、
６２ 等ＦＤＯＡ線（等到達周波数差線）、
６３ 標定位置、
６４ 推定等ＦＤＯＡ線（推定等周波数差線）、
６５、６５Ｂ 電波源位置推定情報、
６６ 許容方向範囲、
６７ 可能位置、
６８ 最尤位置、
６９ 可能方向範囲、
７０ 電波源、
７１ 衛星(＃１)（第１の衛星）、
７２ アンテナ（＃１）、
７３ 衛星(＃２)（第２の衛星）、
７４ アンテナ（＃２）、
７５ リファレンス局（参照局）、
７６ 静止軌道、
７７ 赤道、
７８ 衛星(＃１)７１および衛星(＃２)７３の中間を通る経線、
８１ 電波、
８２ 電波、
８３ 電波、
８４ 電波、
８５ 電波、
８６ 電波、
８７ 電波、
８８ 電波、
９１ 受信信号（第１の受信信号）、
９２ 受信信号（第２の受信信号）、
９３ 受信信号（第３の受信信号）、
９４ 受信信号（第４の受信信号）。
１Ａ，１Ｂ，１Ｃ，１Ｄ 周波数変換部、
２Ａ，２Ｂ，２Ｃ，３Ｄ ＡＤ変換部、
３Ａ，３Ｂ，３Ｃ，３Ｄ 直交検波部、
４、４Ａ 相関関数計算部、
５、５Ａ ピーク検出部、
６ 衛星データ取得部、
７、７Ａ ＴＤＯＡ計算部、
８、８Ａ ＦＤＯＡ計算部、
９ 等ＴＤＯＡ線計算部、
１０ 静止衛星判別部、
１１ 第１の測位部、
１２、１２Ｂ、１２Ｃ 電波源速度推定部、
１３、１３Ｂ、１３Ｃ 第２の測位部、
１４、１４Ａ、１４Ｂ、１４Ｃ 記憶部、
１５ 等ＦＤＯＡ線計算部、
１６、１６Ｂ ＦＤＯＡ修正部、
１７ 推定等ＦＤＯＡ線計算部、
１８ 想定情報入力部、
１９、１９Ｃ 存在確度計算部、
２０ 最尤位置推定部、
３１ 想定位置範囲、
３２、３２Ｂ 想定速度範囲、
３３ 存在判断下限値、
６０ 地球、
６１ 等ＴＤＯＡ線（等到達時間差線）、
６２ 等ＦＤＯＡ線（等到達周波数差線）、
６３ 標定位置、
６４ 推定等ＦＤＯＡ線（推定等周波数差線）、
６５、６５Ｂ 電波源位置推定情報、
６６ 許容方向範囲、
６７ 可能位置、
６８ 最尤位置、
６９ 可能方向範囲、
７０ 電波源、
７１ 衛星(＃１)（第１の衛星）、
７２ アンテナ（＃１）、
７３ 衛星(＃２)（第２の衛星）、
７４ アンテナ（＃２）、
７５ リファレンス局（参照局）、
７６ 静止軌道、
７７ 赤道、
７８ 衛星(＃１)７１および衛星(＃２)７３の中間を通る経線、
８１ 電波、
８２ 電波、
８３ 電波、
８４ 電波、
８５ 電波、
８６ 電波、
８７ 電波、
８８ 電波、
９１ 受信信号（第１の受信信号）、
９２ 受信信号（第２の受信信号）、
９３ 受信信号（第３の受信信号）、
９４ 受信信号（第４の受信信号）。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】