知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2024086501
(43)【公開日】2024-06-27
(54)【発明の名称】相対位置検出システム及び相対位置検出方法
(51)【国際特許分類】
G01S 19/43 20100101AFI20240620BHJP
【ＦＩ】
G01S19/43
【審査請求】未請求
【請求項の数】6
【出願形態】ＯＬ
(21)【出願番号】P 2022201660
(22)【出願日】2022-12-16
(71)【出願人】
【識別番号】518212241
【氏名又は名称】公立大学法人公立諏訪東京理科大学
(74)【代理人】
【識別番号】110002697
【氏名又は名称】めぶき弁理士法人
(74)【代理人】
【識別番号】100104709
【弁理士】
【氏名又は名称】松尾 誠剛
(72)【発明者】
【氏名】小林 誠司
【テーマコード（参考）】
5J062
【Ｆターム（参考）】
5J062CC07
5J062GG02
(57)【要約】
【課題】ＧＮＳＳ衛星からの電波を使って、誤差数ｃｍ程度の高精度位置測位を可能とする、相対位置検出システム及び相対位置検出方法を提供する。
【解決手段】本発明は、ＧＮＳＳ衛星から送られる搬送波位相の瞬時値を用いてアンテナ間の相対位置ベクトルを３次元空間を探索して測位する相対位置検出システム及び相対位置検出方法である。第一段階では粗い刻みで探索を行ない、第二段階では細かい精度で探索することにより高精度測位を短時間で実現する。
【選択図】図１
【解決手段】本発明は、ＧＮＳＳ衛星から送られる搬送波位相の瞬時値を用いてアンテナ間の相対位置ベクトルを３次元空間を探索して測位する相対位置検出システム及び相対位置検出方法である。第一段階では粗い刻みで探索を行ない、第二段階では細かい精度で探索することにより高精度測位を短時間で実現する。
【選択図】図１
【特許請求の範囲】
【請求項1】
複数のアンテナと、
前記複数のアンテナのそれぞれに接続され複数のＧＮＳＳ衛星からの電波を受信して搬送波位相の瞬時値を検出する位相検出装置と、
前記搬送波位相の瞬時値を用いて前記複数のアンテナ間の相対位置ベクトルを測位する測位手段を備え、
前記測位手段は、
前記ＧＮＳＳ衛星の衛星位置を取得する衛星位置情報取得手段と、
前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置を用いて第１推定ベクトルを探索して求める第１の探索手段と、
前記搬送波位相の瞬時値、前記第１推定ベクトル及び前記ＧＮＳＳ衛星の衛星位置を用いて第２推定ベクトルを探索して求める第２の探索手段とを備えることを特徴とする相対位置検出システム
前記複数のアンテナのそれぞれに接続され複数のＧＮＳＳ衛星からの電波を受信して搬送波位相の瞬時値を検出する位相検出装置と、
前記搬送波位相の瞬時値を用いて前記複数のアンテナ間の相対位置ベクトルを測位する測位手段を備え、
前記測位手段は、
前記ＧＮＳＳ衛星の衛星位置を取得する衛星位置情報取得手段と、
前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置を用いて第１推定ベクトルを探索して求める第１の探索手段と、
前記搬送波位相の瞬時値、前記第１推定ベクトル及び前記ＧＮＳＳ衛星の衛星位置を用いて第２推定ベクトルを探索して求める第２の探索手段とを備えることを特徴とする相対位置検出システム
【請求項2】
前記第２の探索手段における第２推定ベクトルの変化刻み（Δ２）が、前記第１の探索手段における第１推定ベクトルの変化刻み（Δ１）よりも小さい値であることを特徴とする請求項１に記載の相対位置検出システム。
【請求項3】
前記第１の探索手段は、前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置に加えて前記相対位置ベクトルの初期値を用いて第１推定ベクトルを探索して求めることを特徴とする請求項１に記載の相対位置検出システム。
【請求項4】
複数のアンテナにそれぞれ接続した複数の位相検出装置を用いて複数のＧＮＳＳ衛星からの電波を受信して搬送波位相の瞬時値を検出し、前記搬送波位相の瞬時値を用いてアンテナ間の相対位置ベクトルを測位する相対位置検出方法であって、
前記測位は、前記複数のＧＮＳＳ衛星の衛星位置を取得する衛星位置取得工程と、
前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置を用いて第１推定ベクトルを探索して求める第１の探索工程と、
前記搬送波位相の瞬時値、前記第１推定ベクトル及び前記ＧＮＳＳ衛星の衛星位置を用いて第２推定ベクトルを探索して求める第２の探索工程とを備えることを特徴とする相対位置検出方法。
前記測位は、前記複数のＧＮＳＳ衛星の衛星位置を取得する衛星位置取得工程と、
前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置を用いて第１推定ベクトルを探索して求める第１の探索工程と、
前記搬送波位相の瞬時値、前記第１推定ベクトル及び前記ＧＮＳＳ衛星の衛星位置を用いて第２推定ベクトルを探索して求める第２の探索工程とを備えることを特徴とする相対位置検出方法。
【請求項5】
前記第２の探索工程における第２推定ベクトルの変化刻み（Δ２）が、前記第１の探索工程における第１推定ベクトルの変化刻み（Δ１）よりも小さい値であることを特徴とする請求項４に記載の相対位置検出方法。
【請求項6】
第１の探索工程は、前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置に加えて前記相対位置ベクトルの初期値を用いて第１推定ベクトルを探索して求めることを特徴とする請求項４に記載の相対位置検出方法。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、相対位置検出システム及び相対位置検出方法及に関する。
【背景技術】
【0002】
地滑り等の災害が発生しそうな危険地帯においては、地盤や構造物のわずかな変位を監視するためにＧＰＳ（Global Positioning System）、またはＧＮＳＳ（Global Navigation Satellite System）を用いた地盤監視システムが提案されている。ここでＧＰＳは米国が運営する人工衛星システムを表す用語であり、ＧＮＳＳは米国を含めて日本、ヨーロッパ、中国などの国々が打ち上げた人工衛星システムを表す用語である。本発明ではＧＮＳＳに統一して説明するが、本発明がＧＰＳに適用できることは言うまでもない。
【0003】
ＧＮＳＳの位置検出方法は、コード情報を検出するコード測位と、搬送波位相を検出する搬送波位相測位の２通りに大別できる。カーナビゲーションなど、広く一般的に用いられているのはコード情報を検出する方法である。ＧＮＳＳ衛星から送られるコード情報は約１ＭＨｚで変化するので、コード測位により数メートル～十メートル程度の精度で地球上の位置（緯度、経度及び高度）を知ることができる。
【0004】
コード測位により位置だけでなく全世界で統一された正確な時刻（ＧＰＳ時刻）も検出できる。ＧＰＳ時刻に閏秒補正を施せば全世界共通のＵＴＣ（Coordinated Universal Time）に変換できる。またＵＴＣに９時間を加えれば日本標準時（JST）に変換できる。
【0005】
コード測位よりも高い位置精度が必要な場合は、より高い周波数（約1.5GHz）の搬送波位相測位を用いることにより、誤差数センチメートル以下の高精度で測位できる。搬送波位相測位としては、ＲＴＫ（Real-Time-Kinematic）、あるいはキャリア位相相対測位の名称で呼ばれる方法がある。
【0006】
搬送波位相測位のシステムとして特許文献１には、複数のＧＮＳＳ衛星から送出される電波の搬送波位相の瞬時値を、複数のＧＮＳＳアンテナとそれらに接続された複数のＧＮＳＳ受信機で受信して所定のＧＰＳ時刻に搬送波位相の瞬時値を取得し、複数のＧＮＳＳアンテナ間の相対位置を検出する相対位置検出システムが開示されている。
【0007】
特許文献１に記載されているシステムにおいては搬送波位相の瞬時値のみを用いるので、計測装置から伝送する情報量を減らすことが可能となり、例えばＬＰＷＡ（Low Power Wide Area）のように、伝送速度が低い無線通信技術を適用して携帯電話回線が圏外となるような場所においても高精度で相対位置検出を可能となる。
【0008】
しかし特許文献１のシステムでは３次元空間を探索する為に多くの演算が必要で、測位時間が長いという問題があった。例えば８立方メートル（2mx2mx2m）の３次元空間を１mmの精度で探索した場合、探索点の総数は８０億点（=2000 x 2000 x 2000）と膨大な数になる。１点の探索にかかる時間が５０ナノ秒であったとしても、８０億点の探索には４００秒と長い時間が必要である。このため探索時間を短縮することが求められていた。
【先行技術文献】
【特許文献】
【0009】
【特許文献1】特開２０２２－１１２０２２号公報
【発明の概要】
【発明が解決しようとする課題】
【0010】
そこで本発明では、搬送波位相の瞬時値を用いて基線ベクトルを探索することにより、測位する際に探索にかかる時間を大幅に短縮することを可能とした相対位置検出システム及び相対位置検出方法を提供することを目的とする。
【課題を解決するための手段】
【0011】
［１］本発明の相対位置検出システムは、複数のアンテナ（２Ａ、２Ｂ）と、前記複数のアンテナのそれぞれに接続され複数のＧＮＳＳ衛星からの電波を受信して搬送波位相の瞬時値を検出する位相検出装置（３Ａ、３Ｂ）と、前記搬送波位相の瞬時値を用いて前記複数のアンテナ間の相対位置ベクトルを測位する測位手段（４）を備え、前記測位手段は、前記ＧＮＳＳ衛星の衛星位置を取得する衛星位置情報取得手段（４３）と、前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置を用いて第１推定ベクトルを探索して求める第１の探索手段（４１）と、前記搬送波位相の瞬時値、前記第１推定ベクトル及び前記ＧＮＳＳ衛星の衛星位置を用いて第２推定ベクトルを探索して求める第２の探索手段（４２）とを備えることを特徴とする相対位置検出システムである。
【0012】
［２］本発明の相対位置検出システムにおいて、前記第２の探索手段における第２推定ベクトルの変化刻み（Δ２）が、前記第１の探索手段における第１推定ベクトルの変化刻み（Δ１）よりも小さい値であることが好ましい。
【0013】
［３］本発明の相対位置検出システムにおいて、
前記第１の探索手段（４１）は、前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置に加えて前記相対位置ベクトルの初期値を用いて第１推定ベクトルを探索して求めることが好ましい。
前記第１の探索手段（４１）は、前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置に加えて前記相対位置ベクトルの初期値を用いて第１推定ベクトルを探索して求めることが好ましい。
【0014】
［４］本発明の相対位置検出方法は、複数のアンテナにそれぞれ接続した複数の位相検出装置を用いて複数のＧＮＳＳ衛星からの電波を受信して搬送波位相の瞬時値を検出し、前記搬送波位相の瞬時値を用いてアンテナ間の相対位置ベクトルを測位する相対位置検出方法であって、前記測位は、前記複数のＧＮＳＳ衛星の衛星位置を取得する衛星位置取得工程（Ｓ３）と、前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置を用いて第１推定ベクトルを探索して求める第１の探索工程（Ｓ４）と、前記搬送波位相の瞬時値、前記第１推定ベクトル及び前記ＧＮＳＳ衛星の衛星位置を用いて第２推定ベクトルを探索して求める第２の探索工程（Ｓ５）とを備えることを特徴とする相対位置検出方法である。
【0015】
［５］本発明の位相検出方法において、前記第２の探索工程における第２推定ベクトルの変化刻み（Δ２）が、前記第１の探索工程における第１推定ベクトルの変化刻み（Δ１）よりも小さい値であることが好ましい。
【0016】
［６］本発明の位相検出方法において、第１の探索工程（Ｓ４）は、前記搬送波位相の瞬時値及び前記ＧＮＳＳ衛星の衛星位置に加えて前記相対位置ベクトルの初期値を用いて第１推定ベクトルを探索して求めることが好ましい。
【発明の効果】
【0017】
本発明の相対位置検出システム及び相対位置検出方法は、刻み値の異なる二段階の探索を実施することにより、演算回数を激減させて短時間で測位を完了させることが可能な、相対位置検出システム及び相対位置検出方法を提供する。すなわち、本発明の相対位置検出システム及び相対位置検出方法は、第一段階では粗い刻みで探索を行ない、第二段階では細かい精度で探索することにより高精度測位を短時間で実現する。
【図面の簡単な説明】
【0018】
【図1】本発明の相対位置検出システム１の構成例を示す図である。
【図2】相対位置検出システム１を構成する位相検出装置３の構成例を示す図である。
【図3】位相検出装置３を構成する瞬時位相計算部３３の構成例を示す図である。
【図4】瞬時位相計算部３３の内部信号波形を模式的に表した図である。
【図5】測位手段４を構成する基線ベクトル探索手段４１の構成例を示す図である。
【図6】測位手段４を構成する基線ベクトル探索手段４２の構成例を示す図である。
【図7】基線ベクトル探索手段４１及び基線ベクトル探索手段４２を構成するEnergy演算手段９の構成例を示す図である。
【図8】実施形態２における相対位置検出方法のフローチャートである。
【図9】測位工程（Ｓ２）を構成する第１の探索工程（Ｓ４）のフローチャートである。
【図10】測位工程（Ｓ２）を構成する第２の探索工程（Ｓ５）のフローチャートである。
【図11】第１の探索工程及び第２の探索工程を構成するＥｎｅｒｇｙ演算工程のフローチャートである。
【図12】本発明を実施した実験結果の例である。
【発明を実施するための形態】
【0019】
本発明の実施形態について説明する。以降では衛星番号をカッコで囲んで表記する。例えば「衛星番号ｎ」を「（ｎ）」として表し、「基準衛星ｒ」を「（ｒ）」と表記する。また２台のアンテナ及び位相検出装置をＡとＢと表示する。例えば位相検出装置２Ａで検出されたＧＮＳＳ衛星Ｓａｔ１の搬送波位相の瞬時値φをφＡ（１）と表記する。以降では４つのＧＮＳＳ衛星（Ｓａｔ１～Ｓａｔ４）を受信できた場合を例として説明するが、より多くのＧＮＳＳ衛星からの電波を受信した場合に適用できることは言うまでもない。
【0020】
〔実施形態１〕
【0021】
図１は、本発明の実施形態１における搬送波位相を用いた相対位置検出システム１の構成例を示す図である。上空にある４つのＧＮＳＳ衛星（Ｓａｔ１～Ｓaｔ４）からの電波が２つのアンテナ（２Ａ及び２Ｂ）で受信され、電気信号に変換される。ここでアンテナ２Ａは基準位置に設置される。２つのアンテナにはそれぞれ位相検出装置（３Ａ及び３Ｂ）が接続される。２台の位相検出装置（３Ａ及び３Ｂ）はコード測位と搬送波位検出を実施し、得られた情報を測位手段４に伝送する。測位手段４は後述する搬送波位相測位を実施し、アンテナ２Ａからアンテナ２Ｂへの相対位置ベクトルＰ２を出力する。（以降では２つのアンテナ間の距離Ｒを「基線長」と呼ぶこともある。）
【0022】
２台の位相検出装置３Ａと３Ｂはそれぞれ、コード測位によりアンテナ２Ａと２Ｂの置かれた概略位置座標（緯度、経度及び高度）を求めて測位手段４に伝送する。コード測位では、アンテナ２Ａと２Ｂが置かれた概略位置座標を、数メートルから数十メートルの誤差で得ることが出来る。また２台の位相検出装置３Ａと３Ｂはそれぞれ、所定のＧＰＳ時刻で搬送波位相をサンプリングし、搬送波位相の瞬時値φＡとφＢを求め、測位手段４に伝送する。（以降ではこのようにして所定のＧＰＳ時刻でサンプリングした搬送波位相の瞬時値を「瞬時位相」と表記する。）
【0023】
［位相検出装置の構成］
図２は、相対位置検出システム１を構成する位相検出装置３の構成例を示す図である。複数の計測地点に設置されている位相検出装置３は、フロントエンド３０、衛星信号受信部３１、コード測位部３２、瞬時位相計算部３３及びタイミング発生手段３４で構成される。このうち衛星信号受信部３１、コード測位部３２及び瞬時位相計算部３３はいずれも衛星毎に複数個設ける。
図２は、相対位置検出システム１を構成する位相検出装置３の構成例を示す図である。複数の計測地点に設置されている位相検出装置３は、フロントエンド３０、衛星信号受信部３１、コード測位部３２、瞬時位相計算部３３及びタイミング発生手段３４で構成される。このうち衛星信号受信部３１、コード測位部３２及び瞬時位相計算部３３はいずれも衛星毎に複数個設ける。
【0024】
図２において受信アンテナ２により衛星電波を受信し、フロントエンド３０は微弱信号をフィルタで抽出し増幅してから低周波数信号に変換し、複数の衛星信号受信部３１に供給する。
【0025】
衛星信号受信部３１は、各衛星から検出される航法メッセージ、擬似距離Ｐｒ、搬送波位相Ω、搬送波の周波数オフセットＤｐを出力する。このうち周波数オフセットＤｐは一般に、「ドップラー周波数」と呼ばれる。コード測位部３２は航法メッセージと擬似距離を用いてコード測位を行い、各衛星の位置（Ｘ，Ｙ，Ｚ）と受信アンテナ２の概略位置座標（緯度・経度・高度）を求める。コード測位部３２はまた、各衛星で受信した電波の時刻誤差εを下式１で求めて出力する。
ここでＲａｎｇｅは各衛星から受信アンテナ２までの距離を表し、各衛星位置と受信アンテナ２の概略位置座標の距離差として求めることができる。Ｃｖは光速を表し、ＳＤは衛星時計の補正パラメータから計算される衛星時計の遅延時間である。衛星時計の補正パラメータは航法メッセージに含まれている。
【数1】
【0026】
タイミング発生手段３４は水晶発振器とカウンター等で構成され、位相検出装置３内部の動作クロックを供給する。動作クロック周波数が高くなると消費電力が増大するので、できるだけ低い周波数の動作クロックが望ましい。タイミング発生手段３４はまた、動作クロックをカウントして装置時刻Ｔｒを出力する。装置時刻Ｔｒには、動作クロックの周期（周波数の逆数）での時刻誤差が含まれる。（動作クロックが１ｋＨｚの場合、最大で１ｍｓの誤差が含まれる）。また装置時刻Ｔｒには、各衛星からの電波がアンテナ２まで到達するに要した時間により、さらに時刻誤差が含まれる。
【0027】
そこで複数の瞬時位相計算部３３は、装置時刻Ｔｒに含まれる時刻誤差の影響を補正し、所定時刻Ｔｓにおいてサンプリングされた搬送波位相（即ち瞬時位相）に補正して出力する。
【0028】
図３は、位相検出装置３を構成する瞬時位相計算部３３の構成例を示すブロック図である。補正時間検出手段３３１は、ＵＴＣに合致している所定時刻Ｔｓ、装置時刻Ｔｒと時刻誤差εから、下式２に従って位相補正時間δを算出して出力する。
【数2】
【0029】
以上により位相補正手段３３２は、下式３に従って瞬時位相φを計算することができる。
ここで関数ｍｏｄは２πを法とする剰余演算である。
【数3】
【0030】
図４は、瞬時位相計算部３３の内部信号波形を模式的に表した図である。いまＧＮＳＳ衛星から送出される搬送波が（Ａ）に示す波形であったとすると、衛星信号受信部３１内部の（図示しない）位相検出回路からは（Ｂ）に示す搬送波位相φが検出される。この搬送波位相は０～２πの範囲に畳み込まれているので、２πから０に変化する毎に２πを加えることにより、（Ｃ）に示す積算位相Ωが得られる。積算位相Ωを所定時刻Ｔｓでサンプリングすることができれば良いが、装置時刻Ｔｒには誤差があり、所定時刻Ｔｓには一致しない。
【0031】
そこで式３は、積算位相が傾きＤｐで直線変化するものと近似して積算位相Ωを補正し、所定時刻Ｔｓにおける瞬時位相φを求めるものである。ここで直線変化ではなく、高次多項式を用いて積算位相を補正することも可能である。
【0032】
図２に示した位相検出装置３は、以上説明した瞬時位相φを測位手段４に伝送する。またコード測位部３２により得られたアンテナ２の概略位置座標を測位手段４に伝送する。瞬時位相及び概略位置座標の伝送は、ＬＰＷＡ等の無線通信手段を用いることもできる。
【0033】
［測位手段の構成］
瞬時位相を使っての高精度測は、測位手段４で実施される。図１に示した測位手段４は基線ベクトル探索手段４１、基線ベクトル探索手段４２、衛星位置情報取得手段４３及びシステムコントローラ４５で構成される。
瞬時位相を使っての高精度測は、測位手段４で実施される。図１に示した測位手段４は基線ベクトル探索手段４１、基線ベクトル探索手段４２、衛星位置情報取得手段４３及びシステムコントローラ４５で構成される。
【0034】
システムコントローラ４５はＣＰＵなどで構成され、測位手段４の動作を司る。システムコントローラ４５は、２台の位置検出装置３Ａ及び３Ｂから得られた概略位置情報（緯度、経度及び高度）を地心直交座標系に変換し、２つのアンテナの相対位置を求めることにより、３次元探索の中心位置となる初期ベクトルＰ０を出力する。システムコントローラ４５はまた、位相検出装置３Ａから供給されるアンテナ２Ａの概略位置情報を地心直交座標系に変換し、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）として出力する。
【0035】
衛星位置情報取得手段４３は、航法メッセージを使ってケプラー方程式を解き、各衛星の位置（Ｘ，Ｙ，Ｚ）を求める。ここで航法メッセージは、各衛星から速度５０bpsで放送されている情報であり、市販ＧＮＳＳ受信機を設置して衛星電波を受信すれば入手できる。あるいは国土地理院等が公開しているインターネットサーバから航法メッセージをダウンロードすることもできる。
【0036】
基線ベクトル探索手段４１は、初期ベクトルＰ０を原点に３次元探索を実施し、第１推定ベクトルＰ１を出力する。この３次元探索では、探索の刻みΔ１を４ｃｍとして、おおざっぱに探索する。このようにして得られた第１推定ベクトルＰ１を中心にして、基線ベクトル探索手段４２は刻みΔ２＝１ｍｍと精度の高い第二段階の探索を行い、第２推定ベクトルＰ２を求め、測位結果として出力する。
【0037】
図５は、測位手段４を構成する基線ベクトル探索手段４１の構成例を示す図である。基線ベクトル探索手段４１は、Ｅｎｅｒｇｙ演算手段９，推定基線ベクトル探索手段４１１、探索刻みメモリ４１２、探索コントローラ４１３及び最小Ｅｎｅｒｇｙ検出手段４１５で構成される。
【0038】
探索刻みメモリ４１２は不揮発メモリなどで構成され、探索刻みΔ１として所定値（４ｃｍ）を出力する。この探索刻みΔ１を最小単位として後述する３次元空間の探索が行われる。探索コントローラ４１３はＣＰＵ等で構成され、３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）を所定範囲で順次変化させ、初期ベクトルＰ０を原点とする３次元空間を探索する。推定基線ベクトル探索手段４１１は、下式４の演算により推定基線ベクトルＰｓを出力する。
【数4】
【0039】
Ｅｎｅｒｇｙ演算手段９は、推定基線ベクトルの確からしさを示す指標としてＥｎｅｒｇｙを出力する。ここでＥｎｅｒｇｙが小さければ小さいほど、推定基線ベクトルが正しい値であると考えられる。最小Ｅｎｅｒｇｙ検出手段４１５は、Ｅｎｅｒｇｙをモニターしながら、Ｅｎｅｒｇｙを最小とした推定基線ベクトルを第１推定ベクトルＰ１として出力する。
【0040】
図６は、測位手段４を構成する基線ベクトル探索手段４２の構成例を示す図である。基線ベクトル探索手段４２は、第１推定ベクトルＰ１を中心として、推定基線ベクトルを刻みΔ２＝１ｍｍで変化させて３次元を探索し、第２推定ベクトルＰ２を求め、測位結果として出力する。基線ベクトル探索手段４２の構成は、探索の刻みΔ２＝が１ｍｍであることと、探索中心が第１推定ベクトルＰ１であることを除いて、基線ベクトル探索手段４１と同一である。
【0041】
図７は、基線ベクトル探索手段４１及び基線ベクトル探索手段４２を構成するEnergy演算手段９の構成例を示す図である。Ｅｎｅｒｇｙ演算手段９は、演算コントローラ９１，二重位相差演算手段９２，格子ベクトル演算手段９３，理論位相差演算手段９４，評価指標演算手段９５及び指標合成演算手段９６で構成される。
【0042】
演算コントローラ９１は、基準衛星（ｒ）と、基準衛星以外の衛星（ｎ）を順次抽出して衛星ペア（ｎ、ｒ）を形成する。例えば基準衛星（ｒ）としてＳａｔ１を選択した場合には、３つの衛星ペア（２，１）、（３，１）、（４，１）を順次抽出する。
【0043】
二重位相差演算手段９２は、基準衛星（ｒ）と、他の衛星（ｎ）から得られた瞬時位相に対して次式５の演算により二重位相差θを演算して出力する。二重位相差を用いることにより電離層における伝搬遅延の変化等、様々な外乱を除去することができる。
【数5】
【0044】
格子ベクトル演算手段９３は、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）から衛星（ｎ）への方向ベクトルｖ（ｎ）を求め、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）から基準衛星（ｒ）への方向ベクトルｖ（ｒ）を求め、２つの方向ベクトルの差分として下式６により格子ベクトルｋを求める。
【数6】
【0045】
理論位相差演算手段９４は、格子ベクトルｋと推定基線ベクトルＰｓとの内積演算（・）により、理論二重位相差Φを下式７で求める。
ここで関数ｍｏｄ（）は２πを法とする剰余演算である。
【数7】
【0046】
以上の処理により得られた理論二重位相差Φと、二重位相差演算手段９２により実測された二重位相差θは、何れも位相情報である。そこで評価指標演算手段９５は下式８のように２つの位相情報の差異を評価して評価指標Ｅを求める。
【数8】
【0047】
指標合成演算手段９６は、下式９に示すように全ての衛星ペアに関する評価指標EをＥｎｅｒｇｙに加算してＥｎｅｒｇｙを求める。
【数9】
【0048】
【0049】
［ステップＳ０］
最初にステップＳ０において、上空にあるＧＮＳＳ衛星（Ｓａｔ１～ＳＡｔ４）が送信する電波を２つのアンテナ（２Ａ及び２Ｂ）で受信し、２つのアンテナにそれぞれ位相検出装置（３Ａ及び３Ｂ）を接続し、それぞれ搬送波位相の瞬時値を検出する。
最初にステップＳ０において、上空にあるＧＮＳＳ衛星（Ｓａｔ１～ＳＡｔ４）が送信する電波を２つのアンテナ（２Ａ及び２Ｂ）で受信し、２つのアンテナにそれぞれ位相検出装置（３Ａ及び３Ｂ）を接続し、それぞれ搬送波位相の瞬時値を検出する。
【0050】
［ステップＳ１］
ステップＳ１において探索の初期ベクトルＰ０をセットする。地滑り等の災害検出に用いる場合は、アンテナ２Ａと２Ｂの設置位置がおよそ解っている。そこで２つのアンテナの位置座標を引き算することにより、初期ベクトルＰ０を求めることができる。あるいは位相検出装置（３Ａ及び３Ｂ）がコード測位を実施して得られる概略位置情報の差として、初期ベクトルＰ０を定めることもできる。
ステップＳ１において探索の初期ベクトルＰ０をセットする。地滑り等の災害検出に用いる場合は、アンテナ２Ａと２Ｂの設置位置がおよそ解っている。そこで２つのアンテナの位置座標を引き算することにより、初期ベクトルＰ０を求めることができる。あるいは位相検出装置（３Ａ及び３Ｂ）がコード測位を実施して得られる概略位置情報の差として、初期ベクトルＰ０を定めることもできる。
【0051】
ステップＳ２において、２つのアンテナから得られた瞬時位相を使って、高精度の搬送波位相測位を実施し、アンテナ２Ａからアンテナ２Ｂへの測位結果Ｐ２を出力する。このステップＳ２で実施される搬送波位相測位の詳細は、ステップＳ３からステップＳ５の３ステップで構成される。
【0052】
［ステップＳ３］
ステップＳ３において、航法メッセージを用いてケプラー方程式を解き、各衛星の位置情報（Ｘ，Ｙ，Ｚ）を取得する。
ステップＳ３において、航法メッセージを用いてケプラー方程式を解き、各衛星の位置情報（Ｘ，Ｙ，Ｚ）を取得する。
【0053】
［ステップＳ４］
ステップＳ４において、第１の探索工程を実施する。即ち初期ベクトルＰ０を中心にし、オフセットベクトルを粗い刻み（Δ１＝４ｃｍ）で変化させてＰ０に加え、３次元空間を探索し第１推定ベクトルＰ１を求める。
ステップＳ４において、第１の探索工程を実施する。即ち初期ベクトルＰ０を中心にし、オフセットベクトルを粗い刻み（Δ１＝４ｃｍ）で変化させてＰ０に加え、３次元空間を探索し第１推定ベクトルＰ１を求める。
【0054】
［ステップＳ５］
ステップＳ５において、第２の探索工程を実施する。即ち第１推定ベクトルＰ１を中心にして、オフセットベクトルを細かい刻み（Δ２＝１ｍｍ）で変化させてＰ１に加え、３次元空間を高精度探索し、第２推定ベクトルＰ２を求めて測位結果とする。
ステップＳ５において、第２の探索工程を実施する。即ち第１推定ベクトルＰ１を中心にして、オフセットベクトルを細かい刻み（Δ２＝１ｍｍ）で変化させてＰ１に加え、３次元空間を高精度探索し、第２推定ベクトルＰ２を求めて測位結果とする。
【0055】
［第１の探索工程］
第１の探索工程（ステップＳ４）のフローチャートを図９に示す。図９は、測位工程（Ｓ２）を構成する第１の探索工程（Ｓ４）のフローチャートである。以降では８ｍ立方の３次元空間を４ｃｍ刻みで３次元探索すること想定し、３つの係数（ＫＸ、ＫＹ、ＫＺ）をそれぞれ ―５０から＋５０まで１０１通りに変化させる例について説明する
第１の探索工程（ステップＳ４）のフローチャートを図９に示す。図９は、測位工程（Ｓ２）を構成する第１の探索工程（Ｓ４）のフローチャートである。以降では８ｍ立方の３次元空間を４ｃｍ刻みで３次元探索すること想定し、３つの係数（ＫＸ、ＫＹ、ＫＺ）をそれぞれ ―５０から＋５０まで１０１通りに変化させる例について説明する
【0056】
［ステップＳ４０］
ステップＳ４０において、３つのスキャン係数（ＫＸ、ＫＹ、ＫＺ）を全て「－５０」として初期化する。
ステップＳ４０において、３つのスキャン係数（ＫＸ、ＫＹ、ＫＺ）を全て「－５０」として初期化する。
【0057】
［ステップＳ４１］
ステップＳ４１において、初期ベクトルＰ０にオフセットベクトル（Δ1KX,Δ1KY,Δ1KZ）を加えて推定基線ベクトルＰｓとする。
ステップＳ４１において、初期ベクトルＰ０にオフセットベクトル（Δ1KX,Δ1KY,Δ1KZ）を加えて推定基線ベクトルＰｓとする。
【0058】
［ステップＳ４２］
ステップＳ４２では、このようにして得られた推定基線ベクトルに対して図１１のフローチャートに記載する工程を施し、Ｅｎｅｒｇｙを算出する。図１１は、第１の探索工程及び第２の探索工程を構成するＥｎｅｒｇｙ演算工程のフローチャートである。
ステップＳ４２では、このようにして得られた推定基線ベクトルに対して図１１のフローチャートに記載する工程を施し、Ｅｎｅｒｇｙを算出する。図１１は、第１の探索工程及び第２の探索工程を構成するＥｎｅｒｇｙ演算工程のフローチャートである。
【0059】
［ステップＳ４３］
ステップＳ４３では、Ｅｎｅｒｇｙがそれまでの最小値であったか否かを判断する。最小値であった場合に処理はステップＳ４４に遷移し、第１推定ベクトルＰ１＝Ｐｓとして第１推定ベクトルＰ１を更新する。
ステップＳ４３では、Ｅｎｅｒｇｙがそれまでの最小値であったか否かを判断する。最小値であった場合に処理はステップＳ４４に遷移し、第１推定ベクトルＰ１＝Ｐｓとして第１推定ベクトルＰ１を更新する。
【0060】
［ステップＳ４５］
ステップＳ４４では、３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）のそれぞれが、―５０から＋５０までの整数値を全て取り終えたかどうかを判断する。全ての値を取り終えていた場合、第１推定ベクトルＰ１を探索結果として処理は終了する。
ステップＳ４４では、３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）のそれぞれが、―５０から＋５０までの整数値を全て取り終えたかどうかを判断する。全ての値を取り終えていた場合、第１推定ベクトルＰ１を探索結果として処理は終了する。
【0061】
［ステップＳ４６］
３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）のそれぞれが―５０から＋５０までの整数値を全て取り終えていない場合、処理はステップＳ４６に遷移する。ステップＳ４６では、３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）の何れかに整数「１」を加え、あるいは初期値―５０に戻すアップデートを行った後に、ステップＳ４１に遷移して処理を続行する。
３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）のそれぞれが―５０から＋５０までの整数値を全て取り終えていない場合、処理はステップＳ４６に遷移する。ステップＳ４６では、３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）の何れかに整数「１」を加え、あるいは初期値―５０に戻すアップデートを行った後に、ステップＳ４１に遷移して処理を続行する。
【0062】
以上述べたように第１の探索工程（ステップＳ４）では、３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）のそれぞれを―５０から＋５０まで変化させ、Ｅｎｅｒｇｙ（Ｐｓ）を最も小さくする第１推定ベクトルＰ１を探索して求める。この探索では、ステップＳ４１～Ｓ４６までのループを約１０３万回（＝１０１の３乗回）評価することになる。
【0063】
［第２の探索工程］
第２の探索工程（ステップＳ５）の詳細構成例を図１０のフローチャートに示す。図１０は、測位工程（Ｓ２）を構成する第２の探索工程（Ｓ５）のフローチャートである。
ステップＳ４とステップＳ５の違いは、初期ベクトルと探索の変化刻みである。即ちステップＳ４においては、初期ベクトルＰ０を中心に探索していたが、ステップＳ４では第１推定ベクトルＰ１を中心とした探索となる。またステップＳ４では推定基線ベクトルＰｓの変化刻みを（Δ１＝４ｃｍ）と大きな値にして、８ｍ立方の広い範囲を粗く探索した。これに対してステップＳ５では推定基線ベクトルＰｓの変化刻みを（Δ２＝１ｍｍ）と小さな値に設定し、５ｃｍ立方の狭い範囲を高精度（１ｍｍ）で探索して第２推定ベクトルＰ２を求めるように構成されている。
第２の探索工程（ステップＳ５）の詳細構成例を図１０のフローチャートに示す。図１０は、測位工程（Ｓ２）を構成する第２の探索工程（Ｓ５）のフローチャートである。
ステップＳ４とステップＳ５の違いは、初期ベクトルと探索の変化刻みである。即ちステップＳ４においては、初期ベクトルＰ０を中心に探索していたが、ステップＳ４では第１推定ベクトルＰ１を中心とした探索となる。またステップＳ４では推定基線ベクトルＰｓの変化刻みを（Δ１＝４ｃｍ）と大きな値にして、８ｍ立方の広い範囲を粗く探索した。これに対してステップＳ５では推定基線ベクトルＰｓの変化刻みを（Δ２＝１ｍｍ）と小さな値に設定し、５ｃｍ立方の狭い範囲を高精度（１ｍｍ）で探索して第２推定ベクトルＰ２を求めるように構成されている。
【0064】
初期ベクトルと探索の変化刻みが異なることを除くと、第１の探索工程（ステップＳ４）と第２の探索工程（ステップＳ５）は同一であるので、各工程の動作説明を省略する。
【0065】
第２の探索工程（ステップＳ５）においても、３次元スキャン係数（ＫＸ、ＫＹ、ＫＺ）のそれぞれを―５０から＋５０まで変化させる。従って第２の探索工程において、ステップＳ５１、Ｓ４２，Ｓ４３，Ｓ４５、Ｓ４６と回るループを約１０３万回（＝１０１の３乗回）評価することになる。
【0066】
［ステップＳ４２の詳細］
ステップＳ４２としてＥｎｅｒｇｙを算出する工程は、図１１に記載したＳ４２１からＳ４２８の８工程で実現される。
ステップＳ４２としてＥｎｅｒｇｙを算出する工程は、図１１に記載したＳ４２１からＳ４２８の８工程で実現される。
【0067】
［ステップＳ４２１］
ステップＳ４２１において、位相検出装置３Ａから供給されるアンテナ２Ａの概略位置情報を地心直交座標系に変換し、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）を求める。
ステップＳ４２１において、位相検出装置３Ａから供給されるアンテナ２Ａの概略位置情報を地心直交座標系に変換し、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）を求める。
【0068】
［ステップＳ４２２］
ステップＳ４２２において、基準衛星（ｒ）と、基準衛星以外の衛星（ｎ）を順次抽出して衛星ペア（ｎ、ｒ）を形成する。例えば衛星数が４つで、基準衛星としてＳａｔ１を選択した場合、３つの衛星ペア（２，１）、（３，１）、（４，１）を順次抽出する。
ステップＳ４２２において、基準衛星（ｒ）と、基準衛星以外の衛星（ｎ）を順次抽出して衛星ペア（ｎ、ｒ）を形成する。例えば衛星数が４つで、基準衛星としてＳａｔ１を選択した場合、３つの衛星ペア（２，１）、（３，１）、（４，１）を順次抽出する。
【0069】
［ステップＳ４２３］
ステップＳ４２３において、基準衛星（ｒ）と、他の衛星（ｎ）から得られた瞬時位相に対して式５の演算を実施して、二重位相差θを求める。
ステップＳ４２３において、基準衛星（ｒ）と、他の衛星（ｎ）から得られた瞬時位相に対して式５の演算を実施して、二重位相差θを求める。
【0070】
［ステップＳ４２４］
ステップＳ４２４において、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）から衛星（ｎ）への方向ベクトルｖ（ｎ）を求め、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）から基準衛星（ｒ）への方向ベクトルｖ（ｒ）を求め、式６により格子ベクトルｋを求める。
ステップＳ４２４において、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）から衛星（ｎ）への方向ベクトルｖ（ｎ）を求め、アンテナ２Ａの３次元概略位置（ＡＸ，ＡＹ，ＡＺ）から基準衛星（ｒ）への方向ベクトルｖ（ｒ）を求め、式６により格子ベクトルｋを求める。
【0071】
［ステップＳ４２５］
ステップＳ４２５において、格子ベクトルｋと推定基線ベクトルＰｓとの内積演算により、理論二重位相差Φを式７の演算により求める。
ステップＳ４２５において、格子ベクトルｋと推定基線ベクトルＰｓとの内積演算により、理論二重位相差Φを式７の演算により求める。
【0072】
［ステップＳ４２６］
理論二重位相差Φと実測された二重位相差θは、何れも位相情報である。そこでステップＳ４２６において、式８に示すように２つの位相情報の差異を評価し、推定相対位置ベクトルＰｓの正しさを表す評価指標Ｅを求める。
理論二重位相差Φと実測された二重位相差θは、何れも位相情報である。そこでステップＳ４２６において、式８に示すように２つの位相情報の差異を評価し、推定相対位置ベクトルＰｓの正しさを表す評価指標Ｅを求める。
【0073】
［ステップＳ４２７］
ステップＳ４２７において、式９に示した演算により評価指標ＥをＥｎｅｒｇｙに加算する。
ステップＳ４２７において、式９に示した演算により評価指標ＥをＥｎｅｒｇｙに加算する。
【0074】
［ステップＳ４２８］
ステップＳ４２８において、可能性のある全ての測位ペア（ｎ，ｒ）についてステップＳ４２２～ステップＳ４２７の演算を行ったか否かを判断する。残りがあった場合に処理はステップＳ４２２に戻り、新たな測位ペアに関する評価指標演算が行われる。全ての測位ペアに関する処理が終わった場合、Ｅｎｅｒｇｙが得られる。
ステップＳ４２８において、可能性のある全ての測位ペア（ｎ，ｒ）についてステップＳ４２２～ステップＳ４２７の演算を行ったか否かを判断する。残りがあった場合に処理はステップＳ４２２に戻り、新たな測位ペアに関する評価指標演算が行われる。全ての測位ペアに関する処理が終わった場合、Ｅｎｅｒｇｙが得られる。
【0075】
以上説明したように本発明の相対位置検出システム及び相対位置検出方法は、３次元探索測位を、刻み値の異なる二段階の３次元探索として実施する。この結果、探索演算の繰り返し回数を激減させる。上述した例では、第一と第２の探索工程を合わせて約２０６万回の評価で測位解を得ることが出来る。これを従来法（一回の探索工程）で実施した場合は約８億回の評価が必要であった。本発明により絶大な演算時間削減効果が得られることが解る。
【0076】
［実験確認］
建物の屋上にアンテナ２Ａと２Ｂを９．８４ｍ離して設置し、本発明を適用して１秒間隔で５００回の測位を行った。この実験結果を図１２に示す。５００回の測位で基線長のバラつきは３．５ｃｍと小さく、正しい基線長が得られたことが解る。またこの実験で、２段階の３次元探索に要した演算時間は、市販のＣＰＵ（3GHz動作のIntel社 Core-i5）を用いた場合に約１００ミリ秒で完了した。この実験により、高速度演算が実現可能であることが確認された。
建物の屋上にアンテナ２Ａと２Ｂを９．８４ｍ離して設置し、本発明を適用して１秒間隔で５００回の測位を行った。この実験結果を図１２に示す。５００回の測位で基線長のバラつきは３．５ｃｍと小さく、正しい基線長が得られたことが解る。またこの実験で、２段階の３次元探索に要した演算時間は、市販のＣＰＵ（3GHz動作のIntel社 Core-i5）を用いた場合に約１００ミリ秒で完了した。この実験により、高速度演算が実現可能であることが確認された。
【0077】
以上に本発明の構成を述べたが、上述の実施形態に限られるものではなく、本発明の要旨を逸脱しない範囲で種々変形実施可能となるものである。例えば、下記に示すような変形実施も可能である。
【0078】
上述の実施形態においては、３次元空間を探索するものとして説明した。しかし例えばアンテナ２Ａが固定され、アンテナ２Ｂが上下方向に動かない場合には、２次元空間を探索するように簡略化して演算速度をさらに高めることも可能である。さらにアンテナ２Ｂの移動方向が予め解っている場合には、１次元空間を探索するように簡略化することもできる。
【0079】
上述の実施形態においては、瞬時位相の単位をラジアンとして説明した。しかし瞬時位相を定数２πで除することにより、位相の一回転を単位とすることもできる。このとき瞬時位相は０～１までの範囲で変化する。
【0080】
上述の実施形態において、基線ベクトル探索手段４１はＥｎｅｒｇｙを最も小さくする第１推定ベクトルＰ１を探索結果として出力するものとした。しかし、例えばＥｎｅｒｇｙ）を最も小さくした第１推定ベクトルＰｓ１と、Ｅｎｅｒｇｙを２番目に小さくした第１推定ベクトルＰｓ２とを求め、アンテナ２Ｂの移動状況など他の情報を使ってＰｓ１、Ｐｓ２のいずれかを選択して探索結果とすることもできる。
【0081】
上述の実施形態においては、衛星から送られる搬送波周波数が約1.5GHzとして説明したが、1.2GHzの搬送波を使う衛星を使うことも可能である。例えば1.5GHzと1.2GHzの２種類の搬送波に対して本発明を適用することも可能であり、さらに高精度の計測が可能になる。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】