(57)【要約】【課題】送信電力の情報がなくても、送信元になる無線装置の位置を推定する技術を提供する。【解決手段】第１取得部１２ａはＴＯＡとＲＳＳを取得する。ＴＤＯＡサンプル測定部１４はＴＯＡからＴＤＯＡを導出する。第１因子グラフ処理部１６は、ＴＤＯＡを使用して、無線装置の位置座標の第１推定値を導出する。特定部１８は、第１推定値の近くに配置された少なくとも４つのセンサ１０を特定する。第２取得部１２ｂはトレーニング信号のＲＳＳを取得する。選択部２０は、第２取得部１２ｂにおいて取得したＲＳＳのうち、特定した少なくとも４つのセンサ１０に対するＲＳＳを選択する。係数計算部２２は、選択したＲＳＳをもとに係数を計算する。第２因子グラフ処理部２４は、第１取得部１２ａにおいて取得したＲＳＳからＤＲＳＳを導出し、導出したＤＲＳＳと係数とを使用しながら、無線装置の位置座標の第２推定値を導出する。【選択図】図１
2018146473.tif 000002
【特許請求の範囲】
【請求項１】
対象となる無線装置からの信号が複数のセンサのそれぞれにおいて受信されており、各センサでのサンプルのＴＯＡ（Ｔｉｍｅ  Ｏｆ  Ａｒｒｉｖａｌ）とＲＳＳ（Ｒｅｃｅｉｖｅ  Ｓｉｇｎａｌ  Ｓｔｒｅｎｇｔｈ）を取得する第１取得部と、
前記第１取得部において取得したサンプルのＴＯＡからセンサ間のサンプルのＴＤＯＡ（Ｔｉｍｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｏｆ  Ａｒｒｉｖａｌ）を導出するＴＤＯＡサンプル測定部と、
前記ＴＤＯＡサンプル測定部において導出したセンサ間のサンプルのＴＤＯＡを使用して、ピタゴラスのＴＤＯＡベースの因子グラフ処理を実行することによって、無線装置の位置座標の第１推定値を導出する第１因子グラフ処理部と、
前記第１因子グラフ処理部において導出した第１推定値の近くに配置された少なくとも４つのセンサを特定する特定部と、
既知の送信装置からのトレーニング信号が複数のセンサのそれぞれにおいて受信されており、各センサでのサンプルのＲＳＳを取得する第２取得部と、
前記第２取得部において取得したＲＳＳのうち、前記特定部において特定した少なくとも４つのセンサに対するＲＳＳを選択する選択部と、
前記選択部において選択したＲＳＳをもとに、係数を計算する係数計算部と、
前記第１取得部において取得したサンプルのＲＳＳからセンサ間のサンプルのＤＲＳＳ（Ｄｉｆｆｅｒｅｎｔｉａｌ  ＲＳＳ）を導出し、導出したＤＲＳＳと前記係数計算部において計算した係数とを使用しながら、ＤＲＳＳベースの因子グラフ処理を実行することによって、無線装置の位置座標の第２推定値を導出する第２因子グラフ処理部と、
を備えることを特徴とする推定装置。
【請求項２】
対象となる無線装置からの信号が複数のセンサのそれぞれにおいて受信されており、各センサでのサンプルのＴＯＡ（Ｔｉｍｅ  Ｏｆ  Ａｒｒｉｖａｌ）とＲＳＳ（Ｒｅｃｅｉｖｅ  Ｓｉｇｎａｌ  Ｓｔｒｅｎｇｔｈ）を取得するステップと、
取得したサンプルのＴＯＡからセンサ間のサンプルのＴＤＯＡ（Ｔｉｍｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｏｆ  Ａｒｒｉｖａｌ）を導出するステップと、
導出したセンサ間のサンプルのＴＤＯＡを使用して、ピタゴラスのＴＤＯＡベースの因子グラフ処理を実行することによって、無線装置の位置座標の第１推定値を導出するステップと、
導出した第１推定値の近くに配置された少なくとも４つのセンサを特定するステップと、
既知の送信装置からのトレーニング信号が複数のセンサのそれぞれにおいて受信されており、各センサでのトレーニング信号に対するサンプルのＲＳＳを取得するステップと、
トレーニング信号に対するサンプルのＲＳＳのうち、特定した少なくとも４つのセンサに対するＲＳＳを選択するステップと、
選択したＲＳＳをもとに、係数を計算するステップと、
対象となる無線装置からの信号に対するサンプルのＲＳＳからセンサ間のサンプルのＤＲＳＳ（Ｄｉｆｆｅｒｅｎｔｉａｌ  ＲＳＳ）を導出し、導出したＤＲＳＳと、計算した係数とを使用しながら、ＤＲＳＳベースの因子グラフ処理を実行することによって、無線装置の位置座標の第２推定値を導出するステップと、
を備えることを特徴とする推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、推定技術に関し、未知の発信源の位置を推定する推定方法およびそれを利用した推定装置に関する。
【背景技術】
【０００２】
無線の位置を推定するために、ＴＯＡ（Ｔｉｍｅ  Ｏｆ  Ａｒｒｉｖａｌ）、ＲＳＳ（Ｒｅｃｅｉｖｅ  Ｓｉｇｎａｌ  Ｓｔｒｅｎｇｈ）等が使用されているが、見通しがなく、かつマルチパス環境である屋内環境での使用には、ＲＳＳの方が適している。ＲＳＳによる処理を簡易にするために、超平面方程式、因子グラフも使用される（例えば、非特許文献１参照）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】Ｃｈｉｎ−Ｔｓｅｎｇ  Ｈｕａｎｇ，  Ｃｈｅｎｇ−Ｈｓｕａｎ  Ｗｕ，  Ｙａｏ−Ｎａｎ  Ｌｅｅ  ａｎｄ  Ｊｉｕｎｎ−Ｔｓａｉｒ  Ｃｈｅｎ、「Ａ  Ｎｏｖｅｌ  Ｉｎｄｏｏｒ  ＲＳＳ−Ｂａｓｅｄ  Ｐｏｓｉｔｉｏｎ  Ｌｏｃａｔｉｏｎ  Ａｌｇｏｒｉｔｈｍ  Ｕｓｉｎｇ  Ｆａｃｔｏｒ  Ｇｒａｐｈｓ」、ＩＥＥＥ  ＴＲＡＮＳＡＣＴＩＯＮＳ  ＯＮ  ＷＩＲＥＬＥＳＳ  ＣＯＭＭＵＮＩＣＡＴＩＯＮＳ、ＪＵＮＥ  ２００９、ＶＯＬ．８、ＮＯ．６、ｐ．３０５０−３０５８
【発明の概要】
【発明が解決しようとする課題】
【０００４】
信号の送信元になる無線装置の位置をＲＳＳをもとに推定する場合、無線装置における送信電力の情報が必要になる。しかしながら、推定対象となる無線装置が、不法無線局や未知発信源である場合、送信電力の情報を取得することが困難である。そこで、２つのセンサのそれぞれにおけるＲＳＳを減算することによって導出されたＤＲＳＳ（Ｄｉｆｆｅｒｅｎｔｉａｌ  ＲＳＳ）をもとに、推定がなされる。ＤＲＳＳでは、無線装置を囲むような少なくとも４つの観測スポットが必要になる。
【０００５】
本発明はこうした状況に鑑みてなされたものであり、その目的は、送信電力の情報がなくても、送信元になる無線装置の位置を推定する技術を提供することである。
【課題を解決するための手段】
【０００６】
上記課題を解決するために、本発明のある態様の推定装置は、対象となる無線装置からの信号が複数のセンサのそれぞれにおいて受信されており、各センサでのサンプルのＴＯＡ（Ｔｉｍｅ  Ｏｆ  Ａｒｒｉｖａｌ）とＲＳＳ（Ｒｅｃｅｉｖｅ  Ｓｉｇｎａｌ  Ｓｔｒｅｎｇｔｈ）を取得する第１取得部と、第１取得部において取得したサンプルのＴＯＡからセンサ間のサンプルのＴＤＯＡ（Ｔｉｍｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｏｆ  Ａｒｒｉｖａｌ）を導出するＴＤＯＡサンプル測定部と、ＴＤＯＡサンプル測定部において導出したセンサ間のサンプルのＴＤＯＡを使用して、ピタゴラスのＴＤＯＡベースの因子グラフ処理を実行することによって、無線装置の位置座標の第１推定値を導出する第１因子グラフ処理部と、第１因子グラフ処理部において導出した第１推定値の近くに配置された少なくとも４つのセンサを特定する特定部と、既知の送信装置からのトレーニング信号が複数のセンサのそれぞれにおいて受信されており、各センサでのサンプルのＲＳＳを取得する第２取得部と、第２取得部において取得したＲＳＳのうち、特定部において特定した少なくとも４つのセンサに対するＲＳＳを選択する選択部と、選択部において選択したＲＳＳをもとに、係数を計算する係数計算部と、第１取得部において取得したサンプルのＲＳＳからセンサ間のサンプルのＤＲＳＳ（Ｄｉｆｆｅｒｅｎｔｉａｌ  ＲＳＳ）を導出し、導出したＤＲＳＳと係数計算部において計算した係数とを使用しながら、ＤＲＳＳベースの因子グラフ処理を実行することによって、無線装置の位置座標の第２推定値を導出する第２因子グラフ処理部と、を備える。
【０００７】
本発明の別の態様は、推定方法である。この方法は、対象となる無線装置からの信号が複数のセンサのそれぞれにおいて受信されており、各センサでのサンプルのＴＯＡ（Ｔｉｍｅ  Ｏｆ  Ａｒｒｉｖａｌ）とＲＳＳ（Ｒｅｃｅｉｖｅ  Ｓｉｇｎａｌ  Ｓｔｒｅｎｇｔｈ）を取得するステップと、取得したサンプルのＴＯＡからセンサ間のサンプルのＴＤＯＡ（Ｔｉｍｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｏｆ  Ａｒｒｉｖａｌ）を導出するステップと、導出したセンサ間のサンプルのＴＤＯＡを使用して、ピタゴラスのＴＤＯＡベースの因子グラフ処理を実行することによって、無線装置の位置座標の第１推定値を導出するステップと、導出した第１推定値の近くに配置された少なくとも４つのセンサを特定するステップと、既知の送信装置からのトレーニング信号が複数のセンサのそれぞれにおいて受信されており、各センサでのトレーニング信号に対するサンプルのＲＳＳを取得するステップと、トレーニング信号に対するサンプルのＲＳＳのうち、特定した少なくとも４つのセンサに対するＲＳＳを選択するステップと、選択したＲＳＳをもとに、係数を計算するステップと、対象となる無線装置からの信号に対するサンプルのＲＳＳからセンサ間のサンプルのＤＲＳＳ（Ｄｉｆｆｅｒｅｎｔｉａｌ  ＲＳＳ）を導出し、導出したＤＲＳＳと、計算した係数とを使用しながら、ＤＲＳＳベースの因子グラフ処理を実行することによって、無線装置の位置座標の第２推定値を導出するステップと、を備える。
【０００８】
なお、以上の構成要素の任意の組合せ、本発明の表現を方法、装置、システム、記録媒体、コンピュータプログラムなどの間で変換したものもまた、本発明の態様として有効である。
【発明の効果】
【０００９】
本発明によれば、送信電力の情報がなくても、送信元になる無線装置の位置を推定できる。
【図面の簡単な説明】
【００１０】
【図１】本発明の実施例に係る推定装置の構成を示す図である。
【図２】図１の推定装置によるＰＴＤＯＡベースの因子グラフ処理の手順を示すフローチャートである。
【図３】図１の推定装置による下り処理の手順を示すフローチャートである。
【図４】図１の推定装置による上り処理の手順を示すフローチャートである。
【図５】図１の推定装置による計算処理の手順を示すフローチャートである。
【図６】図１の推定装置によるＤＲＳＳベースの因子グラフ処理の手順を示すフローチャートである。
【発明を実施するための形態】
【００１１】
本発明を具体的に説明する前に、まず概要を述べる。本発明の実施例は、不法無線局や未知発信源（以下、「無線装置」という）の位置を推定する推定装置に関する。ジオロケーションの近年のトレンドは、低価格、高精度、簡易性を確立することであるが、このジオロケーションは、未知の無線装置の位置を推定することにも使用されつつある。ここで、「未知」は、ターゲットである無線装置のＴＯＡと送信電力が分からないことを示す。無線でのジオロケーションは、電磁波の性質を観測することによって、無線装置の地理的座標を計算する処理として定義される。無線装置の地理的位置を推定するために、ＴＯＡ、ＲＳＳ、ＴＤＯＡ（Ｔｉｍｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｏｆ  Ａｒｒｉｖａｌ）、それらの組合せにおいて電磁波の性質が観測される。
【００１２】
これまで、ＲＳＳとＴＯＡの測定の組合せや、ＲＳＳとＤＯＡ（Ｄｉｒｅｃｔｉｏｎ  Ｏｆ  Ａｒｒｉｖａｌ）の測定の組合せが、ジオロケーションに使用されている。これらの組合せよって、ＲＳＳだけ、ＴＯＡだけ、ＤＯＡだけの場合よりも精度が向上する。しかしながら、方向角を得るために追加の構成が必要になる。また、ＲＳＳとＴＯＡの組合せでは、未知のターゲットに対する絶対的な送信電力とタイムスタンプの情報がなければ、未知の無線装置の位置を推定できない。
【００１３】
本実施例では、未知の無線装置の位置を推定するために、ＰＴＤＯＡ（Ｐｙｔｈａｇｏｒｅａｎ  Ｔｉｍｅ  Ｄｉｆｆｅｒｅｎｃｅ  Ｏｆ  Ａｒｒｉｖａｌ）因子グラフとＤＲＳＳ因子グラフとが組み合わされる。因子グラフのアルゴリズム構造では、グローバル関数が複数のローカル関数に分解される。ローカル関数として、因子グラフは２種類のノードを含む。１つ目は因子ノードであり、因子ノードは、更新したソフト情報（ＳＩ）メッセージを計算する。２つ目は変数ノードであり、変数ノードは、サムプロダクトアルゴリズムを実行し、因子ノードに値を渡す。これら因子ノードと変数ノードとの間では、測定データの平均値と分散値が交換される。
【００１４】
ジオロケーション技術において因子グラフは、ＴＯＡを使用してターゲットの位置を推定するためにサムプロダクトアルゴリズムを実行する。ガウス分布であると想定されるＴＯＡ測定結果の平均値と分散値だけがノードの間で交換される。各センサはＴＯＡを測定し、平均値と分散値が計算される。しかしながら、因子グラフをベースとしたＴＯＡは、ターゲットの完全な時間同期を必要とし、ターゲットが未知である場合にそれは不可能である。そのため、タイムスタンプの知識を不要とするために、因子グラフをベースとしたピタゴラスのＴＤＯＡが導入される。因子グラフをベースとしたピタゴラスのＴＤＯＡは、因子グラフをベースとしたＴＯＡを修正したものであり、因子グラフにノードが追加される。追加されたノードでは、２つのセンサ間のＴＯＡデータの減算が実行される。
【００１５】
因子グラフをベースとしたＲＳＳもジオロケーション技術に使用される。いくつかの観測スポットから送信信号系列が送信され、ＲＳＳがワット単位で取得される。ＲＳＳはｄＢ単位に変換されて因子グラフに入力される。トレーニングによって得られたＲＳＳ情報参照データは、内挿法によってＰＤＰ（Ｐｏｗｅｒ  Ｄｅｃａｙ  Ｐｒｏｆｉｌｅ）を生成するために使用される。ここで、局所的な線形性が、ターゲットの位置を推定するために利用される。これは、屋内において高精度である。
【００１６】
未知の無線装置における絶対的な送信電力の知識を不要とするために、因子グラフをベースとしたＤＲＳＳが提案される。ＤＲＳＳは、２つのセンサのそれぞれにおけるＲＳＳを減算することによって、ＲＳＳを修正した技術である。この減算は、因子グラフに新たなノードの追加を必要とする。しかしながら、ＤＲＳＳでは、少なくとも４つの適切な観測スポットからのＲＳＳ情報を必要とする。それは、ＤＲＳＳのアルゴリズムが差分領域で動作するので、４つかそれ以上の観測スポットがないと線型方程式が縮退してしまって位置が定まらないからである。なお、線型方程式が縮退することは、マトリクスのランクが一つ下がることともいえる。
【００１７】
この観測スポットは、ターゲットに近い位置を有し、ターゲットを囲むように配置されなければならない。ＤＲＳＳ因子グラフでは、そのような適切な観測スポットが既知であり、ターゲットが観測スポットのグリッドの中に存在すべきである。実際は、多くの観測スポットが配置され、それらの中から適切な観測スポットを選択するための別のアルゴリズムを必要とする。このような別のアルゴリズムとして、ＲＡＤＡＲアルゴリズムが使用されるが、ＲＡＤＡＲアルゴリズムはターゲットの送信電力の知識を必要とする。そのため、ここでは、ＲＡＤＡＲアルゴリズムの代わりにＰＴＤＯＡ因子グラフを使用する。ＰＴＤＯＡ因子グラフは、屋外環境におけるＤＲＳＳ因子グラフのための観測スポットを選択するために使用される。
【００１８】
図１は、推定装置１００の構成を示す。推定装置１００は、センサ１０と総称される第１センサ１０ａ、第２センサ１０ｂ、第３センサ１０ｃ、第Ｎセンサ１０ｎ、取得部１２、ＴＤＯＡサンプル測定部１４、第１因子グラフ処理部１６、特定部１８、選択部２０、係数計算部２２、第２因子グラフ処理部２４を含む。取得部１２は、第１取得部１２ａ、第２取得部１２ｂを含む。第１因子グラフ処理部１６は、ユークリッド因子ノード処理部４０、ユークリッド変数ノード処理部４２、ユークリッド距離推定因子ノード処理部４４、サムプロダクトアルゴリズム変数ノード処理部４６、ピタゴラス因子ノード処理部４８、ピタゴラス変数ノード処理部５０、変換因子ノード処理部５２、サムプロダクトアルゴリズム処理部５４、判定部５６、出力部５８を含む。第２因子グラフ処理部２４は、ＲＳＳ因子ノード処理部８０、ＲＳＳ変数ノード処理部８２、ＤＲＳＳ因子ノード処理部８４、ＤＲＳＳ変数ノード処理部８６、線形平面因子ノード処理部８８、出力部９０を含む。サムプロダクトアルゴリズム変数ノード処理部４６は、下り処理部６０、上り処理部６２を含み、ピタゴラス因子ノード処理部４８は、下り処理部６４、上り処理部６６を含み、ピタゴラス変数ノード処理部５０は、下り処理部６８、上り処理部７０を含み、変換因子ノード処理部５２は、下り処理部７２、上り処理部７４を含む。
【００１９】
Ｎ個のセンサ１０は、座標（Ｘｈ，Ｙｈ）に配置される。ここで、ｈはセンサ１０のインデックスを示す。これは、１０００×１０００ｍ２のエリアで、１００ｍのメッシュ構成グリッドを有する。Ｎ個のセンサ１０のうち、４つのセンサ１０が観測スポットとされる。４つの観測スポットはターゲットを囲むように配置され、観測スポットの座標は（ｘｍ，ｙｍ）と示される。ここで、ｍは観測スポットのインデックスを示す。ここで、センサ１０の位置は推定装置１００に知られている。また、ターゲットの座標は（ｘ，ｙ）と示される。ターゲットは、静的な未知の無線装置であり、対象となる無線装置である。ターゲットから送信された信号は、センサ１０に受信される。
【００２０】
第１取得部１２ａは、各センサ１０において受信された信号であって、かつターゲットから送信された信号をもとに、サンプルのＴＯＡを取得する。ＴＯＡは、送信装置から送信され、かつ受信装置で受信された信号間の時間差を示す。ここでは、送信装置がターゲットに相当し、受信装置がセンサ１０に相当する。ＴＯＡは次のように示される。
【数１】
2018146473.tif 000003
ここで、右辺の前者はセンサ１０が信号を受信した受信時間であり、後者はターゲットが信号を送信した送信時間である。未知のターゲットに対して送信時間は取得できないが、受信時間は取得できる。第１取得部１２ａは、サンプルのＴＯＡをＴＤＯＡサンプル測定部１４に出力する。
【００２１】
ＴＤＯＡサンプル測定部１４は、第１取得部１２ａからのサンプルのＴＯＡを受けつけ、サンプルのＴＯＡからセンサ１０間のサンプルのＴＤＯＡを導出する。これは、送信時間の知識を不要にするためであり、具体的には、次のように、センサ１０間のＴＯＡを減算することによって、ＴＤＯＡが導出される。
【数２】
2018146473.tif 000004
ＴＤＯＡの導出はすべてのセンサ１０に対してなされる。
【００２２】
第１因子グラフ処理部１６は、ＴＤＯＡサンプル測定部１４において導出したセンサ１０間のサンプルのＴＤＯＡを使用して、ＰＴＤＯＡベースの因子グラフ処理を実行することによって、無線装置の位置座標の第１推定値を導出する。以下では、第１因子グラフ処理部１６を詳細に説明する。
【００２３】
第１因子グラフ処理部１６におけるユークリッド因子ノード処理部４０は、ＴＤＯＡサンプル測定部１４からＴＤＯＡを受けつける。ユークリッド因子ノード処理部４０は、ＴＤＯＡに電磁波の伝搬速度を乗算することによって、サンプルのＴＤＯＡをユークリッド距離の差に変換する。特に、ユークリッド因子ノード処理部４０のＥｉ，ｊは、Ｋ個のＴＤＯＡをユークリッド距離の差に変換し、平均値と分散値を計算する。これらは、因子グラフにおいて交換されるメッセージである。
【００２４】
ユークリッド変数ノード処理部４２には、２つの因子ノードだけが直接接続される。ユークリッド変数ノード処理部４２は、ユークリッド因子ノード処理部４０のＥｉ，ｊからユークリッド距離推定因子ノード処理部４４のＤｉ，ｊへのメッセージを直接通過させる。このメッセージは、式（３）のような平均値と、式（４）のような分散値の形式で示される。
【数３】
2018146473.tif 000005
ここで、「ｍ」が「平均値」を示し、「σ２」が「分散値」を示し、矢印は、因子グラフのノード間におけるメッセージの進む方向を示す。
【００２５】
ユークリッド距離推定因子ノード処理部４４のＤｉ，ｊは、次のようにユークリッド距離の差を等価なユークリッド距離に変換する。
【数４】
2018146473.tif 000006
【００２６】
サムプロダクトアルゴリズム変数ノード処理部４６のｒｈは、サムプロダクトアルゴリズムによって、ユークリッド距離推定因子ノード処理部４４のＤｉ，ｊとピタゴラス因子ノード処理部４８のＣｈからのメッセージを更新する。添え字ｈは１つ目あるいは２つ目のセンサ１０のいずれかにかかわらないインデックスを示す。また、サムプロダクトアルゴリズム変数ノード処理部４６は、更新したメッセージを目的地ノードに進める。
【００２７】
サムプロダクトアルゴリズム変数ノード処理部４６の下り処理部６０におけるｒｈからピタゴラス因子ノード処理部４８のＣｈへの更新メッセージは式（９）と式（１０）のように示される。一方、サムプロダクトアルゴリズム変数ノード処理部４６の上り処理部６２におけるｒｈからユークリッド距離推定因子ノード処理部４４のＤｉ，ｊへの更新メッセージは式（１１）と式（１２）のように示される。しかしながら、このステップでの最初の繰り返しは、ピタゴラス因子ノード処理部４８のＣｈからのメッセージを示す平均値と分散値の初期値を必要とする。
【数５】
2018146473.tif 000007
ここで、ｈは、ｈ＝ｉであるかｈ＝ｊである。また、ｈ＜ｎに対してａ＝Ｄｈ，ｎ→ｒｈであり、ｈ＞ｎに対してａ＝Ｄｎ，ｈ→ｒｈである。ノードｒｈとＤｉ，ｊとの間のメッセージの交換は、タイムスタンプの知識の必要性を除去するための重要なステップである。
【００２８】
ピタゴラス因子ノード処理部４８のＣｈは、メッセージにピタゴラス演算を実行する。ピタゴラス因子ノード処理部４８の下り処理部６４は、式（１３）から式（１６）で示されるように、変数ノードΔｘｈとΔｙｈの相対的な距離のメッセージを導出する。ピタゴラス因子ノード処理部４８の上り処理部６６は、式（１７）と式（１８）で示されるようにｒｈのメッセージを置き換える。
【数６】
2018146473.tif 000008
ここで、ｍ２Δｘｈ→Ｃｈ、ｍ２Δｙｈ→Ｃｈ、δ２Δｘｈ→Ｃｈ、δ２Δｙｈ→Ｃｈの初期値は、最初の繰り返しにおいて「０」に設定される。
【００２９】
さらに、ピタゴラス因子ノード処理部４８のＣｈにおいて交換されたメッセージの安定性を維持するために、下り処理部６４は、次の制限を加える。
【数７】
2018146473.tif 000009
ここで、εは限定目的のパラメータであり、経験的に選択される。また、上り処理部６６は、次の２つの制限を加える。
【数８】
2018146473.tif 000010
【数９】
2018146473.tif 000011
ここで、δも限定目的のパラメータであり、経験的に選択される。
【００３０】
ピタゴラス変数ノード処理部５０は、ユークリッド変数ノード処理部４２と同様に、相対的な距離のノードΔｘｉとΔｙｉに直接つながる２つのノードだけが存在するので、メッセージはそれらの間で直接交換される。つまり、ピタゴラス因子ノード処理部４８のＣｈと変換因子ノード処理部５２のＡｈとＢｈとの間で交換される。ピタゴラス変数ノード処理部５０の下り処理部６８の処理は、式（１９）と式（２０）に示され、上り処理部７０の処理は、式（２１）と式（２２）に示される。
【数１０】
2018146473.tif 000012
ここで、平均値と分散値は、最初の繰り返しにおけて、初期値「０」にされる。
【００３１】
変換因子ノード処理部５２のＡｈとＢｈは、変数ノードｘｉｎｉｔとｙｉｎｉｔで計算されるターゲットの位置を推定するために、相対的な距離のメッセージを変換する。変換因子ノード処理部５２の下り処理部７２の処理は、式（２３）と式（２４）に示され、上り処理部７４の処理は、式（２５）と式（２６）に示される。
【数１１】
2018146473.tif 000013
【００３２】
サムプロダクトアルゴリズム処理部５４は、繰り返し処理において、変数ノードｘｉｎｉｔとｙｉｎｉｔに対して、式（２７）と式（２８）に示されるサムプロダクトアルゴリズムにしたがってメッセージを加算する。
【数１２】
2018146473.tif 000014
ｘｉｎｉｔに対する式（２７）と式（２８）は、ｙｉｎｉｔにも同様に適用される。
【００３３】
判定部５６は、変換因子ノード処理部５２の処理結果が収束しているかを判定する。ここでは、式（２７）の計算結果、式（２８）の計算結果と、それまでに導出されていた式（２７）の計算結果、式（２８）の計算結果との変化量がしきい値よりも小さくなった場合に、収束と判定される。なお、収束の判定はｙｉｎｉｔに対してもなされるとともに、平均値と分散値のいずれか一方をもとになされてもよい。収束していなければ、判定部５６は、サムプロダクトアルゴリズム処理部５４の処理結果を変換因子ノード処理部５２の上り処理部７４のＡｈとＢｈに戻す。
【００３４】
出力部５８では、判定部５６において繰り返しが収束したと判定された場合、変数ノードｘｉｎｉｔとｙｉｎｉｔが、式（２７）と式（２８）と同様にすべての入ってくるメッセージを加算する。
【数１３】
2018146473.tif 000015
ｙｉｎｉｔに対する演算は、式（２９）と式（３０）と同様である。（ｍｘｉｎｉｔ，ｍｙｉｎｉｔ）の最終的な値は、ＰＴＤＯＡ因子グラフで提供されるターゲットの荒い位置推定値を示し、前述の第１推定値に相当する。
【００３５】
特定部１８は、第１因子グラフ処理部１６において導出した第１推定値の近くに配置された４つのセンサ１０を観測スポットとして特定する。つまり、変数ノード（ｘｉｎｉｔ，ｙｉｎｉｔ）でのＰＴＤＯＡ因子グラフ推定アルゴリズムの結果は、ＤＲＳＳ因子グラフの初期位置として使用され、４つの最も適切な観測スポットを選択するためのターゲットの位置の粗い推定になる。
【００３６】
第２取得部１２ｂは、既知の送信装置からのトレーニング信号が複数のセンサ１０のそれぞれにおいて受信されており、各センサ１０でのサンプルのＲＳＳを取得する。既知の送信装置は、既知の位置座標に配置される。第２取得部１２ｂは、サンプルのＲＳＳに対して平均化処理を実行する。平均化の期間は、十分に長く、かつ十分なサンプルの量になるように設定される。これは、シャドウイングと瞬間的なフェージングの影響を低減することによって、パスロスフェージングだけが残り、エラーフリーにするためである。第２取得部１２ｂは、平均化したサンプルのＲＳＳを選択部２０に出力する。選択部２０は、第２取得部１２ｂにおいて取得したＲＳＳのうち、特定部１８において特定した４つのセンサ１０、つまり観測スポットに対するＲＳＳを選択する。
【００３７】
これまでのＲＴＤＯＡ因子グラフの処理に続いて、ＤＲＳＳ因子グラフの処理が実行される。ＤＲＳＳ因子グラフでは、ターゲット（ｘ，ｙ）の位置推定が、受信信号の電力とＰＤＰとの間の関係を介して得られる。ここで、受信信号の電力はセンサ１０および第１取得部１２ａで得られ、ＰＤＰは前述のトレーニング信号から得られる。また、第２取得部１２ｂにおいて、観測スポットにおいて受信したトレーニング信号には平均がなされる。トレーニング信号とターゲットからの信号のＲＳＳを含むすべてのＲＳＳサンプルはガウス的仮定にしたがう。そのため、これは、十分に長い平均だけによるパスロスモデルによってモデル化される。前述のごとく、十分に長い平均によって、瞬間のフェージング変動とシャドウイングが除外されるからである。その結果は、次のように示される。
【数１４】
2018146473.tif 000016
【００３８】
ここで、Ｐはパスロスのべき指数を示し、ｄＢを単位とする。ｄ、ｄ０、ｆｃ、ｆ、ｎは、ターゲットあるいは観測スポットからｈ番目のセンサ１０までのユークリッド距離（メートル単位）、パスロス指数モデルの参照距離（メートル単位）、キャリアの周波数（Ｈｚ単位）、光の速度（ｍ／ｓ）、パスロスのべき指数の係数をそれぞれ示す。
【００３９】
第１取得部１２ａは、各センサ１０において受信された信号であって、かつターゲットから送信された信号をもとに、サンプルのＲＳＳを取得する。第１取得部１２ａは、瞬間的なフェージング変動とシャドウイングを除去するために、サンプルのＲＳＳを十分長い期間で平均化する。そのため、サンプルのＲＳＳは、多くの独立した要因が集められることによって、測定エラーは、「０」を含まないガウスノイズでモデル化されるので、パスロス指数関数モデルで表現される。
【００４０】
第２因子グラフ処理部２４は、第１取得部１２ａにおいて取得したサンプルのＲＳＳからセンサ１０間のサンプルのＤＲＳＳを導出し、ＤＲＳＳベースの因子グラフ処理を実行することによって、無線装置の位置座標の第２推定値を導出する。なお、ＤＲＳＳベースの因子グラフ処理を実行する際に、ＤＲＳＳと、係数計算部２２において計算した係数とを使用するが、係数計算部２２において計算した係数については後述する。以下では、第２因子グラフ処理部２４を詳細に説明する。
【００４１】
第２因子グラフ処理部２４のＲＳＳ因子ノード処理部８０のＨｐｈは、ＲＳＳサンプルをワット単位からｄＢ単位に変換する。ＲＳＳ変数ノード処理部８２のＮｐｈは、ＲＳＳ因子ノード処理部８０のＨｐｈからＤＲＳＳ因子ノード処理部８４のＧｐｉ，ｊにメッセージを送る。
【数１５】
2018146473.tif 000017
【００４２】
ＤＲＳＳ因子ノード処理部８４のＧｐｉ，ｊは、２つのセンサ１０間においてＲＳＳサンプルを減算し、次のように、ＤＲＳＳサンプルを得る。
【数１６】
2018146473.tif 000018
これは、次のように示してもよい。
【数１７】
2018146473.tif 000019
ＤＲＳＳ因子ノード処理部８４は、ＤＲＳＳ変数ノード処理部８６のＦｐｉ，ｊでの繰り返し処理へのＤＲＳＳ入力メッセージとして使用するために次のような平均値と分散値を計算する。ここで、平均値と分散値は、測定値誤差がガウス分布にしたがうと想定されることによって導出される。
【数１８】
2018146473.tif 000020
【００４３】
ＤＲＳＳ変数ノード処理部８６のＮｐｉ，ｊは、次のように、ＤＲＳＳ因子ノード処理部８４のＧｐｉ，ｊからのメッセージを線形平面因子ノード処理部８８のＦｐｉ，ｊに直接渡す。
【数１９】
2018146473.tif 000021
【００４４】
係数計算部２２は、選択部２０において選択したＲＳＳをもとに、係数を計算する。具体的に説明すると、係数計算部２２は、次のように、４つの適切な観測スポットから送られたトレーニング信号のＲＳＳからセンサ１０間のサンプルのＤＲＳＳを導出する。
【数２０】
2018146473.tif 000022
係数計算部２２は、４つの適切な観測スポットから送られたトレーニング信号をもとに、線形平面最小２乗（ＬＳ）方程式を解くことによってＰＤＰを生成する。線形平面最小２乗（ＬＳ）方程式は次のように示される。
【数２１】
2018146473.tif 000023
ここで、ａｘｉ，ｊ、ａｙｉ，ｊ、ａＰｉ，ｊは、平面方程式の係数である。ｘとｙはターゲットの位置を示す。
【００４５】
また、係数のベクトルａ、観測スポットの座標とＤＲＳＳの行列Ｂ、定数の行列Ｃは、次のように示される。
【数２２】
2018146473.tif 000024
その結果、ＬＳ解は次のように示される。
【数２３】
2018146473.tif 000025
その結果、係数計算部２２において導出される係数は次の通りである。
【数２４】
2018146473.tif 000026
【００４６】
線形平面因子ノード処理部８８のＦｐｉ，ｊは、ターゲットの位置を近似するために線形平面を生成する。線形平面因子ノード処理部８８は、ＤＲＳＳ変数ノード処理部８６から平均値と分散値とを入力するとともに、係数計算部２２から係数も入力する。線形平面因子ノード処理部８８は、入力した平均値、分散値、係数をもとに、無線装置の位置座標ｘ、ｙの平均値と分散値とを次のように導出する。
【数２５】
2018146473.tif 000027
【００４７】
さらに、線形平面因子ノード処理部８８は、導出した平均値と分散値とに対してサムプロダクトアルゴリズムを実行する。
【数２６】
2018146473.tif 000028
式（４６）、式（４７）は、ｘに対するサムプロダクトアルゴリズムを示すが、ｙに対するサムプロダクトアルゴリズムも同様に示される。線形平面因子ノード処理部８８は、サムプロダクトアルゴリズムの実行結果のメッセージが収束するまで、処理を繰り返し実行する。収束についての条件は前述の通りであるので、ここでは説明を省略する。
【００４８】
出力部９０は、線形平面因子ノード処理部８８において収束した平均値と分散値とを入力する。出力部９０は、入力した平均値と分散値とに対してサムプロダクトアルゴリズムを実行する。
【数２７】
2018146473.tif 000029
式（４８）、式（４９）は、ｘに対するサムプロダクトアルゴリズムを示すが、ｙに対するサムプロダクトアルゴリズムも同様に示される。繰り返し回数が終了しなければ、線形平面因子ノード処理部８８に戻される。繰り返し回収が終了した場合、平均値が、ターゲットの位置（ｘ，ｙ）とされる。このターゲットの位置が第２推定値に相当する。
【００４９】
この構成は、ハードウエア的には、任意のコンピュータのＣＰＵ、メモリ、その他のＬＳＩで実現でき、ソフトウエア的にはメモリにロードされたプログラムなどによって実現されるが、ここではそれらの連携によって実現される機能ブロックを描いている。したがって、これらの機能ブロックがハードウエアのみ、ソフトウエアのみ、またはそれらの組合せによっていろいろな形で実現できることは、当業者には理解されるところである。
【００５０】
以上の構成による推定装置１００の動作を説明する。図２は、推定装置１００によるＰＴＤＯＡベースの因子グラフ処理の手順を示すフローチャートである。式（１）、式（２）が実行される（Ｓ１０）。式（３）、式（４）が実行される（Ｓ１２）。式（５）〜式（８）が実行される（Ｓ１４）。第１因子グラフ処理部１６は下り処理を実行する（Ｓ１６）。収束していれば（Ｓ１８のＹ）、式（２９）、式（３０）が実行される（Ｓ２０）。繰り返しが終了すれば（Ｓ２２のＹ）、第１推定値が出力される（Ｓ２６）。収束していない場合（Ｓ１８のＮ）、あるいは繰り返しが終了していない場合（Ｓ２２のＮ）、第１因子グラフ処理部１６は上り処理を実行する（Ｓ２４）。ステップ１４に戻る。
【００５１】
図３は、推定装置１００による下り処理の手順を示すフローチャートである。式（９）、式（１０）が実行される（Ｓ５０）。ｍｒｈ→Ｃｈ＜εでなければ（Ｓ５２のＮ）、式（１３）〜式（１６）が実行される（Ｓ５６）。ｍｒｈ→Ｃｈ＜εであれば（Ｓ５２のＹ）、ｍｒｈ→Ｃｈ＝εとされ（Ｓ５４）、式（１３）〜式（１６）が実行される（Ｓ５６）。式（１９）、式（２０）が実行される（Ｓ５８）。式（２３）、式（２４）が実行される（Ｓ６０）。式（２７）、式（２８）が実行される（Ｓ６２）。
【００５２】
図４は、推定装置１００による上り処理の手順を示すフローチャートである。式（２５）、式（２６）が実行される（Ｓ１００）。式（２１）、式（２２）が実行される（Ｓ１０２）。ｍ２ｒｈ→Ｃｈ−ｍΔｘｈ→Ｃｈ＜０である場合（Ｓ１０４のＹ）、ｍΔｘｈ→Ｃｈ＝ｓｉｇｎ（ｍΔｘｈ→Ｃｈ）×（ａｂｓ（ｍ２ｒｈ→Ｃｈ）−δ）とされる（Ｓ１０６）。これはｙに対してもなされる。ｍ２ｒｈ→Ｃｈ−ｍΔｘｈ→Ｃｈ＜０でない場合（Ｓ１０４のＮ）、ステップ１０６がスキップされる。式（１７）、式（１８）が実行される（Ｓ１０８）。σ２Ｃｈ→ｒｈ＞ｍａｘ（σ２Ｄｈ→ｒｈ）である場合（Ｓ１１０のＹ）、σ２Ｃｈ→ｒｈ＝ｍａｘ（σ２Ｄｈ→ｒｈ）とされる（Ｓ１１２）。σ２Ｃｈ→ｒｈ＞ｍａｘ（σ２Ｄｈ→ｒｈ）でない場合（Ｓ１１０のＮ）、ステップ１１２がスキップされる。式（１１）、式（１２）が実行される（Ｓ１１４）。観測スポットが特定される（Ｓ１１６）。
【００５３】
図５は、推定装置１００による計算処理の手順を示すフローチャートである。式（３７）が実行される（Ｓ１５０）。式（３８）、式（３９）が実行される（Ｓ１５２）。式（４０）が実行される（Ｓ１５４）。式（４１）が実行される（Ｓ１５６）。
【００５４】
図６は、推定装置１００によるＤＲＳＳベースの因子グラフ処理の手順を示すフローチャートである。式（３２）が実行される（Ｓ２００）。式（３３）〜式（３５）が実行される（Ｓ２０２）。式（３６）が実行される（Ｓ２０４）。式（４２）〜式（４５）が実行される（Ｓ２０６）。式（４６）、式（４７）が実行される（Ｓ２０８）。収束していれば（Ｓ２１０のＹ）、式（４８）、式（４９）が実行される（Ｓ２１２）。繰り返しが終了すれば（Ｓ２１４のＹ）、第２推定値が出力される（Ｓ２１６）。収束してない場合（Ｓ２１０のＮ）、あるいは繰り返しが終了していない場合（Ｓ２１４のＮ）、ステップ２０６に戻る。
【００５５】
本発明の実施例によれば、ＰＴＤＯＡを使用するので、タイムスタンプの知識を不要にできる。また、ＰＴＤＯＡを実行するので、ターゲットの位置を粗く推定できる。また、ターゲットの位置が粗く推定されるので、複数のセンサから、４つの観測スポットを選択できる。また、ターゲットを囲むような４つの観測スポットが選択されるので、ＤＲＳＳを実行できる。また、各センサでのサンプルのＲＳＳをもとに、センサ間のサンプルのＤＲＳＳを導出して処理を実行するので、相対的な値を使用できる。また、相対的な使用するので、絶対的な送信電力の情報を不要にできる。また、絶対的な送信電力の情報が不要になるので、送信電力の情報がなくても、ＲＳＳをベースにして送信元になるターゲットの位置を推定できる。また、送信元になるターゲットの位置が推定されるので、不法無線局や未知発信源を探索できる。
【００５６】
以上、本発明を実施例をもとに説明した。この実施例は例示であり、それらの各構成要素の組合せにいろいろな変形例が可能なこと、またそうした変形例も本発明の範囲にあることは当業者に理解されるところである。
【００５７】
本実施例において、特定部１８は４つのセンサ１０を観測スポットとして特定している。しかしながらこれに限らず例えば、特定部１８は、４つ以上のセンサ１０を観測スポットとして特定してもよい。本変形例によれば、推定精度を向上できる。
【符号の説明】
【００５８】
１０  センサ、  １２  取得部、  １４  ＴＤＯＡサンプル測定部、  １６  第１因子グラフ処理部、  １８  特定部、  ２０  選択部、  ２２  係数計算部、  ２４  第２因子グラフ処理部、  ４０  ユークリッド因子ノード処理部、  ４２  ユークリッド変数ノード処理部、  ４４  ユークリッド距離推定因子ノード処理部、  ４６  サムプロダクトアルゴリズム変数ノード処理部、  ４８  ピタゴラス因子ノード処理部、  ５０  ピタゴラス変数ノード処理部、  ５２  変換因子ノード処理部、  ５４  サムプロダクトアルゴリズム処理部、  ５６  判定部、  ５８  出力部、  ６０  下り処理部、  ６２  上り処理部、  ６４  下り処理部、  ６６  上り処理部、  ６８  下り処理部、  ７０  上り処理部、  ７２  下り処理部、  ７４  上り処理部、  ８０  ＲＳＳ因子ノード処理部、  ８２  ＲＳＳ変数ノード処理部、  ８４  ＤＲＳＳ因子ノード処理部、  ８６  ＤＲＳＳ変数ノード処理部、  ８８  線形平面因子ノード処理部、  ９０  出力部、  １００  推定装置。
【図１】
2018146473.tif 000030
【図２】
2018146473.tif 000031
【図３】
2018146473.tif 000032
【図４】
2018146473.tif 000033
【図５】
2018146473.tif 000034
【図６】
2018146473.tif 000035