WO2018180461A1

WO2018180461A1 - 測位システム、測位装置およびコンピュータプログラム

Info

Publication number: WO2018180461A1
Application number: PCT/JP2018/009702
Authority: WO
Inventors: 修一吉川
Original assignee: Kpnetworks; Kpnetworks Ltd; Nidec Corp
Current assignee: Kpnetworks; Kpnetworks Ltd; Nidec Corp
Priority date: 2017-03-28
Filing date: 2018-03-13
Publication date: 2018-10-04
Anticipated expiration: 2019-09-28
Also published as: JPWO2018180461A1

Abstract

測位システム（１）は複数の受信機（２０－Ａ～２０－Ｎ、Ｎ：４または５以上の整数）および演算回路（３１）を備える。演算回路は、各受信機からの入射信号波の時間波形を示すデータからＱ・ｐ＝ｓ（ｐ：基準位置からの波源位置の各座標成分を各成分として含むベクトル、ｓ：基準位置から各受信機までの距離の二乗と各受信機の距離パラメータの二乗との差を並べたベクトル）の関係式を満たす正則正方行列Ｑの逆行列を算出し、当該逆行列をベクトルｓに作用させてベクトルｐ中の波源の各座標成分を算出する。Ｑおよびｐは非線形項を線形項に置換した変数を含む。演算回路は、各受信機および基準受信機（２０－Ａ）の出力時間波形の相関を最大にする各受信機の時間シフト量と信号波の伝搬速度との積である各受信機の距離パラメータの項を、基準距離の項の成分と乗算される行または列の成分として含む行列Ｑの逆行列を求める。

Description

測位システム、測位装置およびコンピュータプログラム

　本開示は、測位システム、測位装置およびコンピュータプログラムに関する。

　特開２０１２－１９８０６６号公報および特開平１１－１６０４０９号公報は、対象物の位置を検出する位置検出装置を開示する。

　近年、Bluetooth（登録商標。以下同じ。）規格を利用した測位システムが普及しつつある。より具体的には、Bluetooth規格のうち、Bluetooth Low Energy（以下「ＢＬＥ」と略記する。）規格の信号を発生させる信号発生器を利用して、当該信号発生器の位置を推定するＢＬＥ測位システムが普及しつつある。

　ＢＬＥ測位システムの多くは３点測位アルゴリズムによる計測処理を採用している。３点測位アルゴリズムによる計測処理は、３つの受信機の各々の設置位置を示す座標情報、および、信号発生器と各受信機との各距離の情報を用いる。測位システムは、所定の演算を行い、３元連立非線形方程式の解として信号発生器の絶対位置座標を出力する。

特開２０１２－１９８０６６号公報特開平１１－１６０４０９号公報

　３点測位アルゴリズムでは３元連立非線形方程式の解を算出することが必要となる。しかしながら３元連立非線形方程式を解くことには困難が伴う。そのため、通常は偏微分行列の逆行列演算とその繰り返し演算を利用するニュートン法が用いられる。

　ニュートン法では適性値への収束判定が必要であり、期待しない別解へ収束することがあり得る。また、ハードウェアやソフトウェアの負担が大きい等の問題がある。

　本願の、限定的ではない例示的なある実施形態は、ニュートン法とは異なる測位アルゴリズムを利用した測位システムを提供する。

　本開示の例示的な実施形態によれば、測位システムは、５台以上の複数の受信機であって、各々が、波源から出力された信号波を含む入射信号波を受信して、受信した前記入射信号波の時間波形のデータを出力する、複数の受信機と、前記複数の受信機の各々から前記時間波形のデータを受け取って、行列Ｑ・ベクトルｐ＝ベクトルｓ（ただし、ベクトルｐ：基準位置からの前記波源の位置の各座標成分を各成分として含むベクトル、ベクトルｓ：前記基準位置から前記各受信機までの距離の二乗と、各受信機の距離パラメータの二乗との差を並べたベクトル）の関係式を満たす、少なくとも五次元の正則正方行列Ｑの逆行列Ｑ^-1を算出し、算出した前記逆行列Ｑ^-1を前記ベクトルｓに作用させて、前記ベクトルｐに含まれる前記波源の位置の各座標成分を算出する演算回路を備え、前記各受信機の距離パラメータを、前記波源の位置から前記複数の受信機のうちの基準受信機までの距離と、前記波源の位置から前記複数の受信機の各々までの距離との差として定義し、前記各受信機の位置の各座標成分と前記波源の位置の各座標成分との差の二乗の和が、前記各受信機と前記波源との距離の二乗に等しいという第１の関係式を、前記各受信機の位置の各座標成分と前記波源の位置の各座標成分との差の二乗の和が、前記波源と前記基準受信機との距離である基準距離と前記各受信機の距離パラメータとの和の二乗に等しいという第２の関係式に変形し、さらに、前記波源の位置の各座標成分を含む二次の項、一次の項、前記基準距離の二乗の項、および、前記基準距離と前記各受信機の距離パラメータとの積の項、の和が、前記基準位置から前記各受信機までの距離の二乗と前記各受信機の距離パラメータの二乗との差、に等しいという第３の関係式に変形し、前記ベクトルｐが、前記二次の項と前記基準距離の二乗の項との差を線形成分に置換した成分、および、前記基準距離の項の成分をさらに含むよう定義したときにおいて、前記演算回路は、受信機ごとに、前記各受信機から出力された前記時間波形のデータを時間軸に沿ってシフトさせながら、前記基準受信機から出力された前記時間波形のデータとの相関を演算し、前記相関を最大にする時間シフト量を決定し、決定した前記時間シフト量と前記波源から出力された信号波の伝搬速度との積を、前記各受信機の距離パラメータとして算出し、前記各受信機の距離パラメータの項を、前記基準距離の項の成分と乗算される行または列の成分として含む前記行列Ｑの逆行列Ｑ^-1を求める。

　本開示の例示的な実施形態によれば、波源の位置を推定するためのアルゴリズムにおいて、演算に含まれる非線形項を線形化する。これにより、ニュートン法を利用する場合と比較して、多重解を回避でき、繰り返し演算を行う必要がなくなる。よってハードウェアへの負荷やソフトウェアの処理を軽減することができる。

図１は、本開示の例示的な実施形態にかかる測位システム１の構成を模式的に示す図である。図２は、測位システム１が導入された環境の第一の例を示す図である。図３は、測位システム１が導入された環境の第二の例を示す図である。図４は、図２および図３における、棚２２および２３を含む４つの棚を示す図である。図５は、棚２５に載置された商品２４ａに貼り付けられた信号発生器１０を示す図である。図６は、棚２２および２３の間を走行する自走車両２４ｂに実装された信号発生器１０を示す図である。図７は、棚２２および２３の間を歩行する買い物客のスマートフォン２４ｃに実装された信号発生器１０を示す図である。図８は、信号発生器１０のハードウェアの構成を示す図である。図９は、測位装置３０のハードウェアの構成を示す図である。図１０は、波源位置Ｐと各受信機との位置関係を示す図である。図１１は、波源位置Ｐから出力された信号波が各受信機２０－Ｋ（Ｋ：Ａ，Ｂ，Ｃ，Ｄ，Ｅ）に到達した時刻ｔ_k（ｋ：ａ，ｂ，ｃ，ｄ，ｅ）を模式的に示す波形図である。図１２は、各受信機２０－Ｋからの出力された信号（出力信号）の時間波形の例を模式的に示す図である。図１３は、ＣＰＵ３１の処理ブロックの構成を示す図である。図１４は、ＣＰＵ３１の処理手順を示すフローチャートである。

　上述した通り、３点測位アルゴリズムにニュートン法を利用した測位システムでは、適正な解に収束しないことがあるほか、さらに、別の問題も存在し得る。すなわち、信号発生器と各受信機との各距離の情報の取得方法である。

　信号発生器と各受信機との各距離の情報は通常は既知ではない。直接、距離情報を得ることはできないため、ＢＬＥ測位システムでは、各受信機が受信した信号発生器からの信号の受信信号レベルを利用する。具体的には、事前に信号発生器と受信機とを１ｍ離れた距離に置き、その時の受信パワー値Ｐ＿ｏを予め取得しておく。その後、信号発生器から信号を受信すると、その時の受信信号レベルＰ＿ｉを利用して計測時の信号発生器と受信機との距離情報を得る。具体的には、計測時の距離情報（ｒ＿ｉ）²＝Ｐ＿ｏ／Ｐ＿ｉにより、求めることができる。

　しかしながら、上述の処理には複数の問題が内包されている。まず、受信パワー値Ｐ＿ｏを予め取得しておくことが現実的には困難である場合がある。たとえば、施設内を移動する買い物客の位置を推定するために測位システムが用いられる場合を想定する。買い物客が所有する電子機器（たとえばスマートフォン）を信号発生器とすると、すべての電子機器を受信機の１ｍの位置に設置して受信パワー値を事前に求めることは現実的ではない。

　利用される信号発生器の種類や数が予め決まっている施設（たとえば工場）内に設置された測位システムの場合には、上述した１ｍの位置の受信パワー値を事前に求めることが可能な場合はある。しかしながら、バッテリの残量等に起因して、事前計測時の信号発生器の送信出力と、測位時の信号発生器の送信出力とが異なる場合があり得る。受信機についても、事前計測時と測位時とで計測状況が異なる場合があり得る。

　また、信号発生器から送信された信号が各受信機に到達するまでの時間（Time Of Flight;TOF）を利用して、信号発生器と各受信機との各距離を求める方法も考えられる。しかしながら、ＴＯＦを求めるためには信号発生器が信号を送信した時刻を正確に知る必要がある。つまり、信号発生器が信号を送信した時刻を各受信機に正確に伝える仕組みが必要になる。さらに、信号発生器と各受信機とは完全に同期した時計で動作する必要があり、そのための仕組みも必要になる。

　上述した問題は、全て測位誤差の要因になる。

　本願発明者は、演算処理の負荷を軽減しつつ、測位誤差を抑え得る測位アルゴリズムについて検討を重ねた。その結果、新たな測位アルゴリズムを構築し、当該測位アルゴリズムを実装した測位システムを完成させるに至った。

　本開示による測位システムでは、５台以上の複数の受信機で、波源から出力された信号波をそれぞれ受信する。波源は、電磁波、音波等の信号波を出力する信号発生器である。

　以下、添付の図面を参照しながら、測位システム、測位装置およびコンピュータプログラムプログラムの構成例を説明する。ただし、必要以上に詳細な説明は省略する場合がある。たとえば、既によく知られた事項の詳細説明や実質的に同一の構成に対する重複説明を省略する場合がある。これは、以下の説明が不必要に冗長になるのを避け、当業者の理解を容易にするためである。なお、発明者は、当業者が本開示を十分に理解するために添付図面および以下の説明を提供するのであって、これらによって特許請求の範囲に記載の主題を限定することを意図するものではない。以下の説明においては、同一または類似する構成要素には、同一の参照符号を付している。

　図１は、本開示の例示的な実施形態にかかる測位システム１の構成を模式的に示す。　測位システム１は、信号発生器１０と、複数の受信機２０－Ａ、２０－Ｂ、・・・、２０－Ｎ（Ｎ：４または５以上の整数）と、測位装置３０とを備えている。

　測位システム１の測位装置３０は、信号発生器１０から出力された信号波を利用して、信号発生器１０の位置を推定する。

　信号発生器１０は、内蔵電池または外部から供給される電力を利用して信号を生成し、信号波として出力する電子回路である。信号発生器１０は、信号波として電磁波または音波を放射する。本明細書では、信号発生器１０を「波源」と呼ぶこともある。本明細書では、信号発生器１０は、ＢＬＥ規格に準拠する２．４ＧＨｚ帯域の電磁波を出力するとして説明する。

　信号発生器１０は、測位対象が人である場合には、人が所有する電子機器に内蔵され得る。または、信号発生器１０は、測位対象が物である場合には当該物に貼り付けられてもよいし、当該物に内蔵されてもよい。

　測位システム１は、測位システム１が設置された空間内で、信号発生器１０の位置を推定する。「空間」とは、本実施形態では主として三次元空間を想定している。図には、Ｘ軸、Ｙ軸およびＺ軸が示されている。ただし、「空間」は二次元空間であってもよい。たとえば、床面を走行する自走車両を測位する場合には、測位システム１は、床面である二次元空間上において、信号発生器１０の位置を推定してもよい。

　複数の受信機２０－Ａ、２０－Ｂ、・・・、２０－Ｎは、各々が、アンテナ装置（図示せず）を内蔵し、当該アンテナ装置を用いて波源から出力された信号波を受信する。アンテナ装置は、上述のＢＬＥ規格に準拠する電磁波を受信可能である。各受信機は、受信した信号波の波形データを出力する。波形データは当該信号波の受信強度の大きさを示す。受信強度データを出力する。なお、アンテナ装置の構成は公知であるため、具体的な構成の説明は省略する。なお、信号発生器１０が音波を発生させる場合には、各受信機はマイクを有する装置であり得る。

　測位装置３０は演算回路（ＣＰＵ）３１を有している。演算回路３１は、複数の受信機２０－Ａ、２０－Ｂ、・・・、２０－Ｎの各々から信号波の波形データを受け取って所定の演算を行い、波源の位置を示す各座標成分を算出する。演算の詳細は後に詳述する。

　図２は、測位システム１が導入された環境の第一の例を示す。図２の例では、測位システム１は、物が載置される棚２２および２３等を有する工場、書店等に構築されている。複数の受信機の各々の位置が、「●」によって示されている。複数の受信機は、概ね分散して配置されている。なお、図２はＸＹ平面図であり、ｚ軸方向の関係は示されていないが、実際には各受信機の位置は高さ方向（ｚ軸方向）の座標も利用して表され得る。

　本実施形態では、複数の受信機は、予め定められた位置に設置されているとする。ある基準位置Ｏを原点とし、３つの軸（Ｘ軸、Ｙ軸およびＺ軸）を定めたとき、各受信機の位置は、Ｘ軸、Ｙ軸およびＺ軸の各値によって表すことができる。以下では、受信機２０－Ａの座標をＡ（ｘ_a，ｙ_a，ｚ_a）などと表記し、ｘ_a，ｙ_aおよびｚ_aをそれぞれ「ｘ座標」「ｙ座標」「ｚ座標」と呼ぶことがある。図面上に例示された受信機２０－Ｊのｘ座標はｘ_J、ｙ座標はｙ_J、ｚ座標はｚ_Jである。

　なお、後述の理由から、本実施形態では演算に利用する複数の受信機の各ｘ座標が全て同じ値になることはないとする。換言すると、測位装置３０は、ｘ座標が少なくとも１つ異なる複数の受信機を利用して演算を行う。ｙ座標、および、ｚ座標についても同様である。すなわち、複数の受信機の少なくとも１つのｙ座標は異なっており、複数の受信機の少なくとも１つのｚ座標も異なっている。

　図３は、測位システム１が導入された環境の第二の例を示す。図面から理解されるように、図３では、複数の受信機が環境内の一部の区画Ｒにまとまって設けられている。なお図３もまたＸＹ平面図であり、ｚ軸方向の関係は示されていないが、実際には各受信機の位置は高さ方向（ｚ軸方向）の座標も利用して表され得る。図３の構成は、上述した位置の条件を満たすように位置が調整されて一体化された複数の受信機のユニットを予め製造しておき、測位システム１の導入環境の天井部分に固定することによって容易に実現される。

　本実施形態にかかる測位処理は、図２および図３のいずれの態様でも利用可能である。

　次に、図４から図７を参照しながら、測位システム１における測位の対象を説明する。

　信号発生器１０は、固定された位置に設置された物に設けられてもよいし、位置が変化し得る物、たとえば自走車両、持ち運び可能な電子機器に設けられてもよい。持ち運び可能な電子機器は、たとえば、携帯電話、スマートフォン、電子タグ装置である。

　図４は、図２および図３における、棚２２および２３を含む４つの棚を示す。それぞれの棚は、物を載置するための複数の棚２５を有している。本実施形態では、複数の棚２５はＹ軸方向およびＺ軸方向に複数存在し得るとする。

　図５は、棚２５に載置された商品２４ａに貼り付けられた信号発生器１０を示す。信号発生器１０から出力された信号波は、５台以上の受信機によって受信される。図５では、受信機２０－Ｓおよび２０－Ｔが信号発生器１０から出力された信号波をそれぞれ受信している状態が示されている。後述の演算によれば、測位システム１は、商品２４ａが載置された位置を推定し得る。

　図６は、棚２２および２３の間を走行する自走車両２４ｂに実装された信号発生器１０を示す。本例でも、受信機２０－Ｓおよび２０－Ｔが、信号発生器１０から出力された信号波をそれぞれ受信している状態が示されている。測位システム１は、自走車両２４ｂの位置をリアルタイムで推定し得る。

　図７は、棚２２および２３の間を歩行する買い物客のスマートフォン２４ｃに実装された信号発生器１０を示す。本例でも、受信機２０－Ｓおよび２０－Ｔが、信号発生器１０から出力された信号波をそれぞれ受信している状態が示されている。測位システム１は、買い物客の位置をリアルタイムで推定し得る。

　図５～図７の例では、信号発生器１０は１つだけ示されているが、複数存在してもよい。ただし、後に詳述するように、個々の信号発生器１０の位置を推定するためには、各信号発生器１０から出力された信号波の波形を識別する必要がある。そのため、各信号発生器１０は、異なる波形を有する信号波を放射することを要する。測位システム１は、並列的に処理を行うことにより、受信した信号波に含まれる、各信号発生器１０を一意に特定する識別情報から各信号発生器１０の位置を推定することができる。

　次に、信号発生器１０および測位装置３０の構成を説明する。

　図８は、信号発生器１０のハードウェアの構成を示す。信号発生器１０は、高周波信号を生成するためのＩＣ１１と、記憶装置１２と、アンテナ１４とを有する。記憶装置１２はたとえばフラッシュＲＯＭであり、信号発生器１０毎に一意の識別情報１３が格納されている。ＩＣ１１は、アンテナ１４を利用して識別情報を周期的に送信する。

　図９は、測位装置３０のハードウェアの構成を示す。

　測位装置３０は、ＣＰＵ３１と、メモリ３２と、通信回路３３とを有しており、これらは内部バスで接続されている。ＣＰＵ３１は、後述の処理により、個々の信号発生器１０の位置を推定し、推定した位置を示す位置情報を生成する演算回路である。メモリ３２は、たとえばＤＲＡＭであり、ＣＰＵ３１の処理に関連して利用されるワークメモリである。

　通信回路３３は、たとえば、１または複数の通信コネクタを有する通信回路である。通信コネクタは、入力インタフェース３４ａ、および、波源の座標のデータを出力する出力インタフェース３４ｂを含む。入力インタフェース３４ａおよび出力インタフェース３４ｂは一体化され１つの通信コネクタとして実装されてもよい。

　入力インタフェース３４ａは、受信機２０－Ａ～２０－Ｎの各々から高周波電気信号を受信する入力端子である。高周波電気信号は、各受信機が受信した電磁波（信号波）を変換して生成された信号である。

　出力インタフェース３４ｂは、たとえばイーサネット（登録商標）規格の有線通信を行う通信端子であり、波源の座標のデータを出力する。波源の座標のデータに代えて、出力インタフェース３４ｂは、波源の座標を画像化した映像信号を出力してもよい。このとき、出力インタフェース３４ｂは、ＤＶＩ端子等の画像信号の出力端子であり得る。図９には、表示装置３５と接続された出力インタフェース３４ｂが示されている。

　次に、本願発明者が完成させた、本実施形態による測位アルゴリズムを説明する。

　本実施形態では５台の受信機を利用する。上述のとおり、受信機２０－Ａの位置Ａを（ｘ_a，ｙ_a，ｚ_a）などと表す。また、波源位置Ｐから受信機２０－Ａまでの距離をｒ_aなどと表す。測位すべき波源（信号発生器１０）の位置Ｐを（ｘ，ｙ，ｚ）と表す。

　図１０は、波源位置Ｐと各受信機との位置関係を示す。図１０には、波源位置Ｐ、各受信機の位置Ａ～Ｅが示されている。また、図１０には基準位置として原点Ｏが示されている。原点Ｏは任意に定め得る。各受信機の位置は、原点Ｏの位置（０，０，０）を基準として決定されている。

　波源位置Ｐと受信機２０－Ａの位置Ａとの間には以下の関係が成り立つ。

　同様に、波源位置Ｐと他の各受信機の位置との間にも以下の関係が成り立つ。

　ここでパラメータβｂ～βｅを導入する。パラメータβｂは、波源位置Ｐから受信機２０－Ｂまでの距離ｒ_bと波源位置Ｐから受信機２０－Ａまでの距離ｒ_aとの差を表す。よって以下のように表記される。なお本明細書では、パラメータβｂ～βｅを「距離パラメータ」と呼ぶことがある。
（数６）
　βｂ＝ｒ_b－ｒ_a

　同様に、βｃ、βｄおよびβｅも、波源位置Ｐから各受信機までの距離ｒ_iと波源位置Ｐから受信機２０－Ａまでの距離ｒ_aとの差を表しており、まとめて以下のように表記される。
（数７）
　βｃ＝ｒ_c－ｒ_a
　βｄ＝ｒ_d－ｒ_a
　βｅ＝ｒ_e－ｒ_a

　数６および７を利用して数２～５を変形すると、以下の数８が得られる。

　数１および数８の各式を展開して変形すると、以下の式が得られる。

　いま、数９に含まれる未知数を、以下のように置き換える。

　上述の変数の置き換えは、二次の項（非線形項）を一次の項（線形項）に変更する操作である。つまり非線形方程式を線形化することに相当する。

　すると、数９は行列を用いて以下のように表現することができる。

　数１１は行列Ｑ、列ベクトルｐおよびｓを用いて以下のように表現できる。
（数１２）
　Ｑ・ｐ＝ｓ

　なお、ベクトルｓは、基準位置である原点Ｏから各受信機までの距離の二乗から、パラメータβｂ～βｅの二乗を減算して並べて得られている。いま、上述したパラメータの定義に準じて、パラメータβａを、波源位置Ｐから受信機２０－Ａまでの距離ｒ_aと、波源位置Ｐから受信機２０－Ａまでの距離ｒ_aとの差を表すパラメータβａを定義する。すると、βａ＝０である。つまり、ベクトルｓの第１行は、基準位置である原点Ｏから受信機Ａまでの距離の二乗から、パラメータβａ²を減算して得られていると言える。ベクトルｓは、基準位置である原点Ｏから各受信機までの距離の二乗から、パラメータβａ～βｅの二乗を減算して並べて得られている、ということができる。

　上述のように定義したパラメータβａを用いると、行列Ｑの（１，５）成分は「２・βａ」であると考えることもできる。

　ベクトルｐに含まれる要素のうち、ｘ、ｙおよびｚが、波源Ｐの位置を表す未知数である。行列Ｑの逆行列をＱ^-1と表すと、ベクトルｐは以下の式により求めることができる。
（数１３）
　ｐ＝Ｑ^-1・ｓ

　数１３は具体的には以下のように表現される。

　数１３および数１４が成立する条件として、行列Ｑの逆行列が存在すること、換言すれば、行列Ｑが正則正方行列であることが必要である。そのためには行列式が０でないことが条件となる。当該条件を満たすためには、（ｉ）行列Ｑのある行と他の行との間には整数倍の関係がないこと、および、（ｉｉ）行列Ｑのある列と他の列との間には整数倍の関係がないこと、が必要である。数１４に示されるように、行列Ｑの行に注目すると、上記関係（ｉ）は明らかに満たさない。上記関係を満たす場合には、５つの受信機のうちの少なくとも２つの位置が少なくとも同じになる必要があるが、そのような配置はありえないからである。

　そこで行列Ｑの列の関係を検討する。５つの受信機のすべてのｘ座標、ｙ座標またはｚ座標が一致する場合には、上記関係（ｉｉ）を満たす。そのため、５つの受信機を選択する際に、そのような条件を満たさない受信機の組を採用しなければよい。たとえば、ｚ座標については、高さが少なくとも１つ異なる５台の受信機を選択する。ｘ座標の関係およびｙ座標の関係も同様である。

　次に、行列Ｑの逆行列を求める際の、パラメータβｂ～βｅの取り扱いを説明する。数６および７に定義されるように、各パラメータは、波源位置Ｐから各受信機までの距離と波源位置Ｐから受信機２０－Ａまでの距離ｒ_aとの差である。当該距離の差は、波源位置Ｐから出力された信号波の到達時間の差に比例する。そこで以下、到達時間の差を求める方法を説明する。

　図１１は、波源位置Ｐから出力された信号波が各受信機２０－Ｋ（Ｋ：Ａ，Ｂ，Ｃ，Ｄ，Ｅ）に到達した時刻ｔ_k（ｋ：ａ，ｂ，ｃ，ｄ，ｅ）を模式的に示す波形図である。波源から全方向に放射された信号波は、伝搬距離に応じて減衰しながら、各受信機２０－Ｋに到達する。各受信機２０－Ｋは各時刻ｔ_k以降継続して信号波を受信する。

　本願発明者は、各受信機２０－Ｋが各時刻ｔ_k以降に受信した信号波の波形には相互に強い相関があることに着目した。たとえば、受信機２０－Ａが時刻ｔ_a以降に受信した信号波の時間波形Ａと、受信機２０－Ｂが時刻ｔ_b以降に受信した信号波の時間波形Ｂとは、振幅を考慮しなければ、概ね同じ形状または同じ位相を有していると言える。

　各受信機２０－Ｋが、図１１に示すような波形を出力することが可能であれば、各時刻ｔ_kを特定することは容易である。しかしながら実際には、各受信機２０－Ｋは、図１１に示すような態様で信号波の時間波形を出力することは事実上できない。その理由は、各受信機２０－Ｋは、信号発生器１０からの信号波以外にも種々の波源から放射された信号波（入射信号波）を絶えず受信し、当該入射信号波の時間波形データを出力しているからである。さらに、種々のノイズが、最終的に得られた時間波形データに重畳され得る。

　図１２は、各受信機２０－Ｋからの出力された信号（出力信号）の時間波形の例を模式的に示す。図示される各時間波形から、信号発生器１０から受信した信号波の時間波形を抽出することは困難である。

　そこで、本願発明者は、取得した１つの時間波形に対して、取得した他の時間波形をそれぞれ時間軸に沿ってシフトさせながら、両者の類似度が最も強くなる時間シフト量Ｕを決定することを考えた。以下、時間波形Ａに対して、時間波形Ｂを時間軸に沿ってシフトさせる例を説明する。

　本願発明者は、類似度を評価するために数１５に示す相関関数ｆ（Ｕ）を導入し、数１６に示すように相関関数ｆ_b（Ｕ）の値を最大化する時間シフト量Ｕを決定することに想到した。

数１５および１６において、ｇ_aは受信機２０－Ａからの出力信号を表す時間関数であり、ｇ_bは受信機２０－Ｂからの出力信号を表す時間関数である。

　決定されたシフト量Ｕは、受信機２０－Ａが信号波を受信した時刻ｔ_aと受信機２０－Ａが信号波を受信した時刻ｔ_bとの時間差（ｔ_b－ｔ_a）を表す。

　信号発生器１０から出力された信号波の伝搬速度をｖとすると、数１７に示す演算により、パラメータβｂ～βｅを得ることができる。
（数１７）
　βｂ＝ｖ・argmax f_b(Ｕ)
　βｃ＝ｖ・argmax f_c(Ｕ)
　βｄ＝ｖ・argmax f_d(Ｕ)
　βｅ＝ｖ・argmax f_e(Ｕ)

　ｆ_c（Ｕ）～ｆ_e（Ｕ）はそれぞれ、数１５において、ｇ_bを、受信機２０－Ｃからの出力信号を表す時間関数ｇ_c、受信機２０－Ｄからの出力信号を表す時間関数ｇ_d、受信機２０－Ｅからの出力信号を表す時間関数ｇ_eに置き換えて得られた相関関数である。

　伝搬速度ｖは、信号波が電磁波で１気圧の大気を伝搬する場合、約３０万キロメートル／秒である。または、伝搬速度ｖは、信号波が音波で、１気圧の乾燥空気を伝搬する場合、約３３１．５メートル／秒である。

　数１４の右辺に含まれる、ｘ_k，ｙ_k，ｚ_k（ｋ：ａ，ｂ，ｃ，ｄ，ｅ）は、各受信機２０－Ｋ（Ｋ：Ａ，Ｂ，Ｃ，Ｄ，Ｅ）の位置の各座標値であり、既知である。よって数１４の右辺の逆行列Ｑ^-1およびベクトルｓは全て演算によって求めることができる。これにより、数１４の左辺の成分ｘ、ｙ、ｚ、つまり波源の位置（ｘ，ｙ，ｚ）を推定することができる。

　なお、本明細書では、各受信機が、信号発生器１０からの信号を受信して信号を出力するまでに要する時間（処理時間）は実質的に等しいとする。等しくない場合であっても、他の各受信機の処理時間から、受信機Ａの処理時間を減算した差の情報を予め取得し、得られた値Ｕに加えればよい。これにより、信号発生器１０からの信号が到達した時間差を正しく求めることができる。

　測位装置３０のＣＰＵ３１（図９）は、上述の原理に従った演算を実行する。以下、図１３および図１４を参照しながら、ＣＰＵ３１の動作を説明する。

　図１３は、ＣＰＵ３１の処理ブロックの構成を示す。また図１４は、ＣＰＵ３１の処理手順を示すフローチャートである。

　ＣＰＵ３１は、パラメータ演算器５１と、逆行列演算器５２と、ベクトル乗算器５３として機能する。図１３では、あたかも３つの構成要素が存在するように示されているが、実際には処理の単位を意味している。本実施形態では、ＣＰＵ３１は、図１４に示すフローチャートによる処理を行うコンピュータプログラムにしたがって動作する。ＣＰＵ３１は、コンピュータプログラムの命令に従い、時間によって、パラメータ演算器５１として動作し、逆行列演算器５２として動作し、ベクトル乗算器５３として動作する。処理ブロックから処理ブロックへの矢印は、データが次の演算に利用されることを意味している。

　ただし、パラメータ演算器５１、逆行列演算器５２およびベクトル乗算器５３の少なくとも１つをハードウェアによって実現してもよい。たとえばＦＰＧＡ（Field-programmable gate array）を用いることにより、パラメータ演算器５１、逆行列演算器５２およびベクトル乗算器５３を１つの集積回路上に実装することが可能である。

　以下では説明の便宜のため、パラメータ演算器５１、逆行列演算器５２およびベクトル乗算器５３が独立に設けられた構成要素であるとして説明する。

　パラメータ演算器５１は、受信機２０－Ａ～２０－Ｅの各々から出力信号のデータを受け取り、数１７に示す演算を行うことにより、パラメータβｂ～βｅを算出する（ステップＳ１）。算出した各パラメータは逆行列演算器５２に送られる。

　逆行列演算器５２は、パラメータ演算器５１からパラメータβｂ～βｅを受け取る。また逆行列演算器５２は、メモリ３２に格納されている各受信機の位置座標を示すデータを読み出す（ステップＳ２）。そして逆行列演算器５２は、パラメータβｂ～βｅおよび各受信機の位置座標を示すデータを用いて行列Ｑの逆行列を計算する（ステップＳ３）。逆行列の演算方法は公知であり、たとえば掃き出し法(sweep out method)を用いて計算することができる。なお、行列Ｑを入力すると逆行列Ｑ^-1を出力するコンピュータプログラムのライブラリもまた公知であり容易に入手可能である。そのようなライブラリプログラムを用いてＣＰＵ３１を逆行列演算器５２として動作させてもよい。

　ベクトル乗算器５３は、パラメータ演算器５１からパラメータβｂ～βｅを受け取る。また、ベクトル乗算器５３は、逆行列演算器５２によって計算された逆行列Ｑ^-1およびメモリ３２に格納されている各受信機の位置座標を示すデータを受け取る。後者は、原点位置から各受信機までの距離の二乗を並べたベクトルｓを求めるために利用される。受け取ったデータに基づいて、ベクトル乗算器５３は、演算Ｑ^-1・ｓを行う（ステップＳ４）。これにより、数１４左辺のベクトルｐを求めることができる。ベクトル乗算器５３は、求めたベクトルｐの成分であるｘ，ｙ，ｚを、波源の位置として出力する（ステップＳ５）。

　以上の処理により、波源である信号発生器１０の位置を推定することができる。

　上述の数１０に示す操作により、非線形項が線形化され、未知数が増えた。さらに、ベクトルｐの要素に未知数ｒ_aを導入した。これらにより、波源の位置を示すｘ、ｙ、ｚに加えて未知数が２つ増加した。５つの未知数を求めるためには、５つの連立方程式が必要となる。これが、５台の受信機からの信号を利用して、５行５列の行列Ｑが必要になった理由である。逆行列を求めるために、本願発明者は受信機の数を増やして行列のランクの減少を補った。

　アルゴリズムを線形化することにより、ニュートン法を利用する場合と比較して、多重解を回避でき、繰り返し演算を行う必要もないため、ハードウェアへの負荷やソフトウェアの処理を軽減することができる。ＢＬＥ測位システムの導入は比較的手軽に行うことができる。そのような性格上、導入のコストを抑えることが求められる。線形計測処理は、入力に対して出力が一意的に決定するためテーブル化も可能でありシステムの低コスト化に貢献する。

　「テーブル化」とは、予め測定されたβｂ、βｃ、βｄ、βｅの組と、それらから推定されるｘ、ｙ、ｚの組とを対応付けたテーブルを用意することを言う。実際に出力信号の時間波形データが得られ、βｂ、βｃ、βｄ、βｅの組が決定されたとき、テーブルを参照して一致する組を検索する。一致する組がテーブルに存在した場合には、対応付けられたｘ、ｙ、ｚの組を読み出して出力する。行列演算等が不要になるため、非常に高速に結果を出力することができる。さらにＣＰＵ３１の負荷を大きく抑制できる。テーブルのエントリ数が増えるほど、実測したβｂ、βｃ、βｄ、βｅの組に一致する組がテーブル上に存在する可能性が高くなる。

　なお、テーブルのエントリ数を増加させることは極めて容易である。測位システム１が導入された環境で、波源となる信号発生器１０の位置を変えながら、都度、各受信機で信号波を受信してβｂ、βｃ、βｄ、βｅの組を取得し、波源の位置ｘ、ｙ、ｚの組を推定しておけばよいからである。検索のヒット率が十分高いテーブルを設けておけば、比較的処理能力が低いＣＰＵを採用しても、十分高速に波源の位置を示す情報を出力することができる。

　なお、上述の説明では、信号発生器１０の位置を推定するためには５台の受信機からの出力信号の時間波形データを利用すると説明したが、６台以上の受信機からの出力信号の時間波形データを利用してもよい。

　上述した実施形態では、ベクトルｐおよびｓを列ベクトルとしたときの行列Ｑを説明した。しかしながら、ベクトルｐおよびｓを行ベクトルとして行列Ｑを構築してもよい。転置の関係にあれば数学的には同じであることは明らかである。

　上述の処理によれば、信号発生器が受信機に信号発生時刻を伝えるシステムが不要となり、ＴＯＦ方式の欠点が解消される。

　本開示の誘導システムは、屋内または屋外を移動する移動体の位置の推定に用いられ得る。また測位した結果を利用する移動体の位置の制御に用いられ得る。

　１　測位システム、　１０　信号発生器、　２０－Ａ、２０－Ｂ、・・・、２０－Ｎ　受信機、　３０　測位装置、　３１　ＣＰＵ（演算回路）、　３２　メモリ、　３３　通信回路、　３４ａ　入力インタフェース、　３４ｂ　出力インタフェース、　５１　パラメータ演算器、　５２　逆行列演算器、　５３　ベクトル乗算器

Claims

　５台以上の複数の受信機であって、各々が、波源から出力された信号波を含む入射信号波を受信して、受信した前記入射信号波の時間波形のデータを出力する、複数の受信機と、
　前記複数の受信機の各々から前記時間波形のデータを受け取って、
　行列Ｑ・ベクトルｐ＝ベクトルｓ
（ただし、ベクトルｐ：基準位置からの前記波源の位置の各座標成分を各成分として含むベクトル、ベクトルｓ：前記基準位置から前記各受信機までの距離の二乗と、各受信機の距離パラメータの二乗との差を並べたベクトル）
の関係式を満たす、少なくとも五次元の正則正方行列Ｑの逆行列Ｑ^-1を算出し、算出した前記逆行列Ｑ^-1を前記ベクトルｓに作用させて、前記ベクトルｐに含まれる前記波源の位置の各座標成分を算出する演算回路を備え、
　　前記各受信機の距離パラメータを、前記波源の位置から前記複数の受信機のうちの基準受信機までの距離と、前記波源の位置から前記複数の受信機の各々までの距離との差として定義し、
　　前記各受信機の位置の各座標成分と前記波源の位置の各座標成分との差の二乗の和が、前記各受信機と前記波源との距離の二乗に等しいという第１の関係式を、
　　前記各受信機の位置の各座標成分と前記波源の位置の各座標成分との差の二乗の和が、前記波源と前記基準受信機との距離である基準距離と前記各受信機の距離パラメータとの和の二乗に等しいという第２の関係式に変形し、さらに、
　　前記波源の位置の各座標成分を含む二次の項、一次の項、前記基準距離の二乗の項、および、前記基準距離と前記各受信機の距離パラメータとの積の項、の和が、前記基準位置から前記各受信機までの距離の二乗と前記各受信機の距離パラメータの二乗との差、に等しいという第３の関係式に変形し、
　　前記ベクトルｐが、前記二次の項と前記基準距離の二乗の項との差を線形成分に置換した成分、および、前記基準距離の項の成分をさらに含むよう定義したときにおいて、
　前記演算回路は、受信機ごとに、
　　前記各受信機から出力された前記時間波形のデータを時間軸に沿ってシフトさせながら、前記基準受信機から出力された前記時間波形のデータとの相関を演算し、
　　前記相関を最大にする時間シフト量を決定し、
　　決定した前記時間シフト量と前記波源から出力された信号波の伝搬速度との積を、前記各受信機の距離パラメータとして算出し、
　　前記各受信機の距離パラメータの項を、前記基準距離の項の成分と乗算される行または列の成分として含む前記行列Ｑの逆行列Ｑ^-1を求める、測位システム。
　前記行列Ｑが五次元の正則正方行列であるときにおいて、前記基準位置を（０，０，０）、前記波源の位置を（ｘ，ｙ，ｚ）、前記基準距離をｒ_a、前記各受信機の位置を（ｘ_i，ｙ_i，ｚ_i）、前記各受信機の距離パラメータをβｉとすると（ただしｉ＝ａ，ｂ，ｃ，ｄ，ｅ）、
　前記行列Ｑ、前記ベクトルｐ、および前記ベクトルｓは数１のように表現される、請求項１に記載の測位システム。
　前記ベクトルｐに含まれる前記ωは、前記基準位置から前記波源までの距離の二乗と前記基準距離の二乗との差を表し、
　前記行列Ｑおよび前記ベクトルｓに含まれる前記β_b、β_c、β_dおよびβ_eは、受信機ごとに決定した、前記時間シフト量と前記波源から出力された信号波の伝搬速度との積である、請求項２に記載の測位システム。
　前記波源は、前記信号波として電磁波または音波を放射する、請求項１から３のいずれかに記載の測位システム。
　前記演算回路が算出した前記波源の座標を出力するインタフェースをさらに備えた、請求項１から４のいずれかに記載の測位システム。
　測位システムに用いられる測位装置であって、
　前記測位システムは、５台以上の複数の受信機であって、各々が、波源から出力された信号波を含む入射信号波を受信して、受信した前記入射信号波の時間波形のデータを出力する、複数の受信機を有し、
　前記測位装置は、
　前記複数の受信機の各々から前記時間波形のデータを受け取る入力端子と、
　行列Ｑ・ベクトルｐ＝ベクトルｓ
（ただし、ベクトルｐ：基準位置からの前記波源の位置の各座標成分を各成分として含むベクトル、ベクトルｓ：前記基準位置から前記各受信機までの距離の二乗と、各受信機の距離パラメータの二乗との差を並べたベクトル）
の関係式を満たす、少なくとも五次元の正則正方行列Ｑの逆行列Ｑ^-1を算出し、算出した前記逆行列Ｑ^-1を前記ベクトルｓに作用させて、前記ベクトルｐに含まれる前記波源の位置の各座標成分を算出する演算回路と、
　前記演算回路が算出した前記波源の位置の各座標成分のデータを出力する出力端子と
を備え、
　　前記各受信機の距離パラメータを、前記波源の位置から前記複数の受信機のうちの基準受信機までの距離と、前記波源の位置から前記複数の受信機の各々までの距離との差として定義し、
　　前記各受信機の位置の各座標成分と前記波源の位置の各座標成分との差の二乗の和が、前記各受信機と前記波源との距離の二乗に等しいという第１の関係式を、
　　前記各受信機の位置の各座標成分と前記波源の位置の各座標成分との差の二乗の和が、前記波源と前記基準受信機との距離である基準距離と前記各受信機の距離パラメータとの和の二乗に等しいという第２の関係式に変形し、さらに、
　　前記波源の位置の各座標成分を含む二次の項、一次の項、前記基準距離の二乗の項、および、前記基準距離と前記各受信機の距離パラメータとの積の項、の和が、前記基準位置から前記各受信機までの距離の二乗と前記各受信機の距離パラメータの二乗との差、に等しいという第３の関係式に変形し、
　　前記ベクトルｐが、前記二次の項と前記基準距離の二乗の項との差を線形成分に置換した成分、および、前記基準距離の項の成分をさらに含むよう定義したときにおいて、
　前記演算回路は、受信機ごとに、
　　前記各受信機から出力された前記時間波形のデータを時間軸に沿ってシフトさせながら、前記基準受信機から出力された前記時間波形のデータとの相関を演算し、
　　前記相関を最大にする時間シフト量を決定し、
　　決定した前記時間シフト量と前記波源から出力された信号波の伝搬速度との積を、前記各受信機の距離パラメータとして算出し、
　　前記各受信機の距離パラメータの項を、前記基準距離の項の成分と乗算される行または列の成分として含む前記行列Ｑの逆行列Ｑ^-1を求める、測位装置。
　測位システムの測位装置に設けられたコンピュータによって実行されるコンピュータプログラムであって、
　前記測位システムは、５台以上の複数の受信機であって、各々が、波源から出力された信号波を含む入射信号波を受信して、受信した前記入射信号波の時間波形のデータを出力する、複数の受信機を有し、
　前記コンピュータプログラムは、前記コンピュータに対し、
　前記複数の受信機の各々から取得された前記時間波形のデータを示すデータを読み出すステップと、
　行列Ｑ・ベクトルｐ＝ベクトルｓ
（ただし、ベクトルｐ：基準位置からの前記波源の位置の各座標成分を各成分として含むベクトル、ベクトルｓ：前記基準位置から前記各受信機までの距離の二乗と、各受信機の距離パラメータの二乗との差を並べたベクトル）
の関係式を満たす、少なくとも五次元の正則正方行列Ｑの逆行列Ｑ^-1を算出するステップと、
　算出した前記逆行列Ｑ^-1を前記ベクトルｓに作用させて、前記ベクトルｐに含まれる前記波源の位置の各座標成分を算出するステップと、
　前記演算回路が算出した前記波源の位置の各座標成分のデータを出力するステップと
　を実行させ、
　　前記各受信機の距離パラメータを、前記波源の位置から前記複数の受信機のうちの基準受信機までの距離と、前記波源の位置から前記複数の受信機の各々までの距離との差として定義し、
　　前記各受信機の位置の各座標成分と前記波源の位置の各座標成分との差の二乗の和が、前記各受信機と前記波源との距離の二乗に等しいという第１の関係式を、
　　前記各受信機の位置の各座標成分と前記波源の位置の各座標成分との差の二乗の和が、前記波源と前記基準受信機との距離である基準距離と前記各受信機の距離パラメータとの和の二乗に等しいという第２の関係式に変形し、さらに、
　　前記波源の位置の各座標成分を含む二次の項、一次の項、前記基準距離の二乗の項、および、前記基準距離と前記各受信機の距離パラメータとの積の項、の和が、前記基準位置から前記各受信機までの距離の二乗と前記各受信機の距離パラメータの二乗との差、に等しいという第３の関係式に変形し、
　　前記ベクトルｐが、前記二次の項と前記基準距離の二乗の項との差を線形成分に置換した成分、および、前記基準距離の項の成分をさらに含むよう定義したときにおいて、
　前記逆行列Ｑ^-1を算出するステップは、
受信機ごとに、
　　前記各受信機から出力された前記時間波形のデータを時間軸に沿ってシフトさせながら、前記基準受信機から出力された前記時間波形のデータとの相関を演算し、
　　前記相関を最大にする時間シフト量を決定し、
　　決定した前記時間シフト量と前記波源から出力された信号波の伝搬速度との積を、前記各受信機の距離パラメータとして算出し、
　　前記各受信機の距離パラメータの項を、前記基準距離の項の成分と乗算される行または列の成分として含む前記行列Ｑの逆行列Ｑ^-1を求める、コンピュータプログラム。