JP7279801B2 - 電波発信源位置推定システム及び電波発信源位置推定方法 - Google Patents

Description

本開示は、電波の発信源の位置を推定する手法と、それを利用した電波発信源位置推定方法に関する。

近年、無線技術の進展と普及に伴って周波数資源がひっ迫しており、周波数を時間・空間・周波数領域で有効に活用する重要性が高まっている。そのため、一般に各国の電波監督機関が無線の利用者に周波数の割り当てを行い、許可を受けた周波数や電波強度の範囲内で電波を利用する方式がとられてきた。しかし、無線局免許を取得せずに電波利用する不正無線局が過大な出力で電波を送信したり、免許を受けた無線局が意図せず免許を受けていない他の帯域で電波を送信するなどして、電波干渉や通信障害を引き起こす問題が起こっている。これに対し、日本では総務省が非特許文献１に示される電波監視システム(DEURAS : Detect Unlicensed Radio Stations)を全国に配備し、電波の強度・到来方向を測定して不法無線局の位置を推定し、電波の適正な利用を図ってきた。これは例えば、推定した電波の発信源の位置を地図情報を参照して不法無線局か免許を取得した合法的な無線局かの判断をしたり、推定した位置の付近を可動式の電波センサで詳細探索したりする対処方法が考えられる。

電波発信源の位置推定方法として、例えば非特許文献２や特許文献１では電波センサで測定した電波強度Received Signal Strength Indicator (RSSI)を用いた方法が提案されている。一般に、電波の受信強度は電波センサが発信源から離れるにつれて、すなわち電波の伝搬距離が増大するにつれて減少する。非特許文献２では伝搬距離と受信強度との関係式を設定し、対象エリア内の各地点に発信源が存在する場合の受信電力の期待値を地点ごとに算出し、電波センサの受信電力の測定値が上述の期待値に最も適合する位置を発信源の位置と推定する手法が提案されている。

しかし、都市部においては建物や物体によって電波の反射・遮蔽が発生するため受信強度は距離の単純な関数とはならず、実際の測定値とその地点の期待値との間に差が生じ、位置推定精度が大きく劣化する。そこで特許文献１では電波センサから見た方位ごとに遮蔽物やブレークポイントの位置を設定し、領域ごとに伝搬距離と受信強度との関係式の減衰パラメータを設定して期待値を算出する方法が提案されている。尚、これらの手法は電波センサの測定物理量が到来時間Time of Arrival (ToA)や到来時間差Time Differential of Arrival (TDoA)や到来方向Angle of Arrival (AoA)である場合でも適用できる。

また、都市環境での建物や物体による反射・遮蔽の影響をより詳細に考慮する手法が特許文献２で提案されている。特許文献２は電波強度ではなくチャネルインパルス応答（CIR）を測定する点と、それに伴い伝搬距離と受信強度との関係式ではなく位置指紋を用いる点や、適合度合いを評価する方法が異なるが、対象エリアの任意の地点から電波が送信された場合の測定値の期待値を取得し、測定値が最も適合する位置を発信源の推定位置とする点で共通している。この先行例では、都市環境の複雑な電波伝搬状態を高確度で把握するため、車両に電波発信源を搭載してトレーニング信号を送信しながら対象となるエリア内を走行し、電波センサでこのトレーニング信号を受信する。そして、測定データをクリギングで空間補間して、対象エリアの任意の地点から電波が送信された場合の測定値の期待値を取得している。

このように、建物や物体による電波の反射・遮蔽が発生する環境において電波の発信源の位置を推定するにあたり、電波センサの測定値の期待値を対象エリア内の各点に設定して、電波センサの測定値に最も適合する位置を発信源の位置とする方式が提案されている。

特開2017-67529号公報 特許第6399512号公報

総務省 電波監視システム http://www.tele.soumu.go.jp/j/adm/monitoring/moni/type/deurasys/ 原晋介：位置推定における統計的推定理論，IEICE Fundamentals Review, 4-1, 32/38 (2010)

以下の分析は、本願発明者によって与えられたものである。特許文献１や特許文献２のように電波センサの測定値の期待値を都市の建物や物体やブレークポイントに対応して不連続に設定すると、同程度の期待値になる地点が不連続に分布する。特許文献2の方式では一見連続的に期待値が分布するが、同程度の期待値になる地点が大きく離れたところにも存在する。すると、電波センサの測定値に同程度に適合する位置が大きく離れた位置に同時に発生するが、その中で最も適合する位置を最終的な推定位置として出力する。一方、実際の測定では電波の発信源の送信出力や伝搬状態が頻繁に変動するため、電波センサの測定値も変動する。すると、その測定値に同程度に適合する複数の地点のなかで最も適合する地点が遷移することがあり、それらの地点が大きく離れているために推定位置が突然不連続に大きく変化する、という課題がある。この場合、発信源が実際に高速で大きく移動した場合と推定誤差が大きい場合との切り分けができないため、位置推定した発信源の分析が難しくなる。例えば、推定対象の発信源が移動体に搭載された発信源なのか移動していない発信源なのか、の判断が難しく、また、免許を取得した無線局なのか不法無線局なのかの判断が難しくなる。尚、位置推定の結果に対して時系列に移動平均やカルマンフィルターなどの時系列フィルタ処理をしても、複数の候補地点の中間の位置が推定結果として出力されて誤差が大きくなるため、時系列フィルタは効果的ではない。

本開示は、上記の課題を鑑み、建物や物体による電波の反射・遮蔽が発生する都市環境において、推定した電波発信源の位置の変動が大きい場合でも正確に位置を把握する方法を提供することを目的とする。

上記の課題を解決するため、本開示の電波発信源位置推定装置は、既知位置の参照発信源から送信された電波の電波特徴量測定値に基づき該既知位置から電波が送信された場合の電波特徴量測定値の期待値を取得する学習データ生成部と、該既知位置の期待値を合成することにより対象エリアの任意の位置から電波が送信された場合の測定値の期待値を算出する空間分布合成部と、該期待値と電波センサで測定した受信電力との適合度合いを基に対象エリアの各地点に電波源が存在する尤度分布を計算する尤度算出部と、尤度算出部が算出した尤度分布をもとに尤度が最大値や極大値になる位置を算出して発信源の位置を推定する位置推定部と、対象エリアにおける尤度の空間分布と発信源の推定位置を表示する表示部と、を備える。

さらに、表示する極大値は、最大値との尤度の差分が設定した閾値の範囲内の極大値に限定する。さらに、それらの極大値のうち、指定した範囲内の期間において尤度が最大値になったことのある極大値を準発信源位置として表示し、それ以外の極大値を参考発信源位置として表示する。

また、指定した範囲内の期間において尤度が最大値になった場所を表示する。

また、本開示の電波発信源位置推定方法は、指定した複数の時刻において対象エリアの各地点における尤度を計算し、地点ごとに複数の時刻の尤度を統計処理し、その各地点の統計処理後の尤度の空間分布を表示する。

本開示によれば、建物や物体による電波の反射・遮蔽が発生する都市環境において、推定した電波発信源の位置が不連続に大きく変化する場合でも安定して位置を把握することができる。

本実施の形態の電波発信源位置推定方法の構成図である。 本実施の形態の電波発信源位置推定システムの構成図である。 本実施の形態の位置推定装置の概略構成を表す図である。 本実施の形態の位置推定装置の処理構成を表す図である。 本実施の形態の学習時のシステムの運用を説明するための図である。 本実施の形態の参照発信源から電波を送信していない場所の値を補完する方法を説明する図である。 本実施の形態の学習用データをもとに普通クリギングで空間補完して求めた対象領域のすべてのグリッドの期待値分布である。 本実施の形態の尤度算出部の処理構成を表す図である。 本実施の形態の尤度分布の例を表す図である。 本実施の形態の位置推定部の処理構成を表す図である。 本実施の形態の尤度分布を色の濃淡で表示するとともに、推定した発信源位置、準発信源位置、参考発信源位置をそれぞれ異なるマーカーで表示した図である。 本実施の形態の学習フェーズの処理を表す図である。 本実施の形態の位置推定フェーズの処理を表す図である。 本実施の形態の出力された尤度分布の一例を表す図である。 本実施の形態のセンサの配置を表す図である。 本実施の形態のセンサ毎に合成した期待値の空間分布を表す図である。 本実施の形態の尤度分布の一例を表す図である。 本実施の形態の結合尤度分布の一例を表す図である。 本実施の形態の推定された尤度分布の一例を表す図である。 本実施の形態の受信電力の測定値と、その参照発信源とセンサとの距離との関係をプロットした図である。 本実施の形態の期待値分布の一例を表す図である。 本実施の形態の各電波センサの尤度分布の例を表す図である。 本実施の形態の結合尤度分布の例を表す図である。

本開示の実施の形態について図面を参照して以下、詳細に説明する。

[第１の実施の形態]
はじめに、第1の実施の形態の電波発信源位置推定システムの構成について説明する。

図１は第１の実施の形態の電波発信源位置推定システムの運用を説明するための図である。位置推定の対象領域に電波センサが配置されており、電波特徴量を測定する。電波特徴量としては、受信した電波の強度RSSIや、チャネルインパルス応答（CIR）や、ToA、ＴＤｏＡ、ＡｏＡなどを用いることができる。この実施の形態では電波特徴量がＲＳＳＩである場合を想定して記述する。この場合、検出対象の周波数の電波を受信してその受信強度を記録する。各電波センサは時刻同期されており、その受信強度は時刻情報とともに解析サーバーに転送する。データの送信は、Wi-Fi等の無線LANやLTE等のモバイル網等の無線NWを介してもよく、またEthernet等の有線NWを介してもよい。電波センシング処理サーバーでは、電波センサから送信される受信強度が所定の閾値を上回った時に、以下に示す手法で電波発信源の位置推定を行う。

図２は第1の実施の形態の電波発信源位置推定システムの構成である。解析サーバー３０内に記憶装置４０と位置推定装置５０が内蔵されている。センサ２０と記憶装置４０が有線または無線で接続されており、センサ２０で測定された電波特徴量が保存される。また、学習時に参照発信源位置測定装置１０で測定された参照発信源の位置情報が記憶装置４０に保存される。参照発信源位置測定装置１０はＧＰＳなどのＧＮＳＳ受信機を用いることができ、測位衛星からの信号を受信して自己位置を推定し、装置内に測位情報を保存する。測定した測位情報は逐次記憶装置４０に転送してもよく、また学習が全て終了した後に一括で記憶装置４０に転送しても良い。測位情報の転送は、ＵＳＢ等の有線接続でもよく、またＷｉ－Ｆｉなどの無線ＬＡＮを用いても良い。そして、記憶装置４０に保存された参照発信源の位置とセンサ２０の測定値を位置推定装置５０が取り込んで、位置推定を行う。

図３は位置推定装置５０の概略構成を表す図である。位置推定装置５０には、前処理部５１と、学習データ生成部５３及び空間分布合成部５４から成る期待値生成部５２と、尤度算出部５５と、位置推定部５６と、表示部５７とが含まれる。

図４は位置推定装置５０の処理構成を表す図である。

まず前処理部５１が、前述の記憶装置４０から電波センサのモード情報と電波特徴量測定値や参照発信源位置情報を取り込む。モード情報はこれらのデータが学習用か位置推定用かのモードを指定する情報である。また、測定値が所定の閾値を下回った場合は電波が送信されていないとしてその時刻の測定データを除去したり、測定値の線形単位と対数単位の相互変換や正規化や標準化などの前処理や、隣接する時刻の測定データを時系列で移動平均したり、ノイズフィルタリングや異常値除去などのデータクレンジング処理を含んでも良い。そして、モード情報に従って、入力データを学習用か推定用かに分類する。

学習データ生成部５３は、時系列データである電波特徴量測定値と参照発信源位置情報を比較し、同時刻のデータを結合する。これにより、学習時における電波特徴量測定値とその時の参照発信源の位置が紐づけられる。また、発信源位置が所定の距離以内にある複数の測定データを平均などの統計処理によって集約するなどの前処理を含んでもよい。

空間分布合成部５４では、解析対象領域において学習データが存在しない地点に参照発信源が位置して電波を送信した場合の電波特徴量の測定値の期待値を、学習データから合成して算出する。これにより、解析対象領域を所定の間隔で分割した全地点について、その地点から参照発信源が電波を送信した場合の各センサの測定値の期待値を取得する。これを期待値の空間分布と呼ぶことにする。

図５から図７を用いて、期待値の空間分布の算出方法を説明する。図５は学習時のシステムの運用を説明するための図である。図１と同様に電波センサを配置し、電波特徴量を測定する準備を整える。そして、参照発信源から電波を送信させながら対象領域内を移動する。この時、参照発信源と共にＧＮＳＳ受信機を保持し、参照発信源の位置を測位する。この参照発信源から送信される電波を全センサで測定し、解析サーバーに送信して保存する。この測定値が学習用データであり、多くの位置について収集することが望ましい。そのため、実際の運用では参照発信源とＧＮＳＳ受信機を車両に搭載し、対象となるエリア内の道路を車両で移動しながら送信する。しかしながら、車両で走行できるのは車道のみであり、歩道や建物がある場所などは走行できないため、対象エリア内のすべての地点において参照発信源から送信することはできない。また、対象エリアが野原など何も障害物がない場所であったとしても、すべての地点から参照信号を発信することは現実的ではない。このように、参照発信源から電波を送信していない場所が発生するので、その地点の値を補完する必要がある。

図６は参照発信源から電波を送信していない場所の値を補完する方法を説明する図である。まず、対象領域を小さなグリッドに分割する。その後、センサ毎に以下の処理を行う。通常、それぞれのグリッド内の複数の地点から参照信号が送信されるので、分割したそれぞれのグリッド内から参照発信源が電波を送信した時のすべての測定値の代表統計量をそのグリッドの期待値として算出する。代表統計量は、平均値、最大値、中央値、７５パーセンタイル値、などを使用することができる。その値を図６では色の濃淡で表示する。

次に、図７に示すように、参照発信源が通過しておらず値が得られていないグリッドの期待値を、周囲の値が得られているグリッドの期待値で空間補完する。空間補間の方法としては、線形補完、立体補完、スプライン補完、クリギングなどの方法があるが、地理情報システムで広く利用される空間補間技術であるクリギングが最も適していると考えられる。クリギングの種類としては、測定値をそのまま補完する普通クリギングでもよく、また、数１で示される伝搬距離と受信強度との関係式をトレンドとして用いた不偏クリギング（または傾向付きクリギング）でもよい。不偏クリギングの場合は、学習用データが存在するグリッドの期待値と数式１の結果の差分について普通クリギングを施し、補完した結果に数式１の電力の和を期待値として出力することに相当する。ここで、(x,y)は電波発信源の位置座標、(x_n,y_n)は、電波センサｎの位置座標、d_n(x,y)は、電波センサｎと電波発信源との距離、(α,β)は伝搬定数を表す。図７は図６の学習用データをもとに普通クリギングで空間補完して求めた、対象領域のすべてのグリッドの期待値分布である。これが１つのセンサについての期待値分布であり、この処理をすべてのセンサについて行う。

図４に戻る。尤度算出部５５では、位置推定用の電波特徴量測定値と期待値の空間分布とを比較して、その地点に電波発信源が存在する尤度を解析領域の全グリッドについて算出する。この尤度は、期待値と測定値の適合度合いを表す指標の一種であり、より一般化する場合には尤度は適合度に含まれる。尤度の一例として、卓越した直接波が存在せずにたくさんの散乱波だけが受信されるマルチパスフェージングであるレイリーフェージングの場合、その受信強度を２乗した物理量の確率密度分布は指数関数となり、その尤度はこの指数関数を用いて得られる。図８は尤度算出部５５の処理構成を表す図である。センサnが測定した電波特徴量P_nと、センサｎの全グリッドの期待値の空間分布

とすると、対象領域の任意の位置(x,y)に発信源が存在する場合にセンサｎがP_nの受信電力の信号を受信する尤度

は数２で求められる。これを対象領域の全グリッドについて算出すれば、電波センサｎの発信源位置の尤度分布が得られる。全電波センサを考慮した結合尤度分布は各電波センサの尤度分布を乗算することによって得られるが、真数として扱うと最小値と最大値の差が大きく、結合尤度の値を地図上に色の濃淡で表すと最大値だけが強調されて他の極大値の位置を把握することが難しくなるため、数３のようにまず各センサの尤度を対数に変換し、全センサを考慮した尤度として対数尤度の和である結合対数尤度Lを算出する。これを対象領域の全グリッドについて算出すれば、全センサを考慮した結合対数尤度分布L(x,y)が得られる。これより後では、この結合対数尤度分布を尤度分布と呼ぶ。尤度分布の例を図９に示す。

ここで注意すべきこととして、数２の尤度と数３や図９の尤度分布は、発信源の電力が学習フェーズで使用した参照発信源の電力と等しい場合における尤度である。実際に位置推定の対象となる未知の発信源の送信電力は参照発信源の電力とは異なる場合がほとんどであるため、正確に位置を推定するためには位置推定の対象の発信源の送信電力も推定する必要がある。ここで、参照発信源の送信電力P_TX0と推定対象の発信源の送信電力P_TXの差分P_TX-P_TX0をΔPとする。すると、参照発信源との送信電力の差分がΔPの発信源が対象領域の任意の位置(x,y)に存在する場合にセンサｎがP_nの電力を受信する尤度

は数４で求められ、全センサを考慮した結合対数尤度分布L(x,y,ΔP)は数５で求められる。

尤度算出部５５では、ΔPを所定の範囲・間隔で掃引し、各ΔPにおける尤度分布を算出し、これら複数の尤度分布をまとめて電力・尤度分布として位置推定部５６に出力する。図９に示した尤度分布は単一のΔPについての尤度分布であるので、掃引したΔPの数だけそれぞれのΔPにおける尤度分布が存在し、それらをまとめて電力・尤度分布として出力する。

図４に戻り、位置推定部５６では、算出した電力・尤度分布をもとに発信源の位置を推定する。図１０は位置推定部５６の処理構成を表す図である。位置推定部５６は、時系列統計処理部６１と、電力推定部６２と、発信源位置推定部６３と、極大値算出部６５と、記憶部６４と、準発信源位置抽出部６６と、を有する。

時系列統計処理部６１は、位置推定する対象時刻の前後の時刻を含む尤度分布の時系列データを統計処理する。具体的には、各グリッドの尤度の値について設定した時間幅に渡って平均や重みづけ平均やカルマンフィルターなどのフィルタリング処理を行い、その結果をその時刻の尤度として出力する。これらの時系列フィルタリング処理は、それぞれの電力の差分値ΔPにおける尤度分布それぞれに対して行う。また、この処理は必ずしも行う必要はなく、推定対象の時刻の尤度の値をそのまま出力しても良い。

電力推定部６２は、推定対象の発信源の送信電力を推定し、その送信電力における尤度の空間分布を出力する。まず、時系列統計処理部６１が出力した電力・尤度分布のうち、全電力・グリッドにおける尤度の中での最大値を検索する。そして、それを与えるΔPを求め、そのΔPにおける全グリッドの尤度分布を出力する。

発信源位置推定部６３は、電力推定部６２が出力した尤度分布において、全グリッドの中で尤度が最大になるグリッドを抽出し、その位置を発信源位置として出力する。同時に、その位置と、尤度の計算に用いた測定値の測定時刻とを記憶部６４に保存する。

極大値算出部６５は、時系列統計処理部６１が出力した尤度分布において、全グリッドの尤度を隣接するグリッドの尤度と比較し、隣接するグリッドの尤度よりも大きい尤度を有するグリッドとその尤度を抽出する。そして、これらの極大値から最大値を判定し、最大値以外の極大値のうち最大値との差分が所定の閾値以内の極大値を準発信源位置候補として、その位置を出力する。図９の例では、グリッド１とグリッド２とグリッド３が極大値になっており、このうちグリッド２が最大値なので、グリッド１とグリッド３を極大値の位置として出力する。そして、グリッド１と３の極大値の尤度と、グリッド２の最大値の尤度との差分を計算し、所定の閾値以内であった場合に準発信源位置候補として出力する。尚、最大値以外の極大値が無い場合や、最大値との差が閾値以内の極大値が無い場合は、何も出力しない。

準発信源位置抽出部６６は、まず、位置推定の対象時刻から設定した時間を過去に遡った期間の発信源位置の情報を記憶部６４から取り込む。そして、準発信源位置候補が極大値算出部６５から出力された場合に、それら準発信源位置候補のうち過去の発信源位置と同じグリッドの準発信源位置候補を、準発信源位置として出力する。このとき、過去の発信源位置と厳密に同じ位置でなくともよく、過去の発信源位置から所定の距離以内にある準発信源位置候補を準発信源位置に含めても良い。そして、その他の準発信源位置候補を参考発信源位置として出力する。これら出力する発信源位置の確度は、発信源位置が最も高く、次に準発信源位置が高く、参考発信源位置が最も確度の低い発信源位置の候補である。図９の例では、グリッド１が設定した時間の範囲内の過去において尤度が最大になり発信源位置として推定した位置であり、準発信源位置と判定されたとする。グリッド２が尤度が最大であるので発信源位置として判定され、グリッド３は設定した時間の範囲内の過去においては尤度が最大になったことがないので参考発信源位置と判定されたとする。

図４に戻り、表示部５７は、入力された尤度分布と推定位置を、ディスプレイやプロジェクターなどの表示装置に表示する。ここで言う推定位置には、上記の図１０の発信源位置や準発信源位置や参考発信源位置を含む。図１１はその一例で、図９に示した尤度分布を色の濃淡で表示するとともに、推定した発信源位置、準発信源位置、参考発信源位置をそれぞれ異なるマーカーで表示している。また、ここでは対象時刻の発信源推定位置だけを表示しているが、設定した期間内の過去の推定位置も同時に表示してもよい。その場合、対象時刻から過去に遡った推定位置ほどマーカーの色を薄くする、などの処理を施しても良い。

次に、第1の実施の形態の電波発信源位置推定システムの処理について説明する。まず処理の概略を説明すると、本システムは、学習フェーズと推定フェーズに分けられる。学習フェーズでは、まず参照発信源から電波を送信しかつその位置を測位しながら対象領域内を移動し、その信号をセンサで受信する。そしてその測定値と測位情報をもとに対象領域内の全地点について、それぞれの地点から電波が送信された場合の各センサの測定値の期待値を合成し、期待値の空間分布を出力する。推定フェーズでは、その期待値の空間分布と、センサの測定値をもとに推定対象の電波発信源の位置を推定する。

図１２に、学習フェーズの処理を示す。まず、電波センサが測定した参照発信源の電波特徴量の測定値と、参照発信源の位置情報を取り込む（S１１）。次に、同時刻の測定値と参照発信源の位置情報を結合することにより、学習用データを生成する（S１２）。次に、空間分布の合成においては、まず解析領域をグリッドに分割し、参照発信源が通過したグリッドについては代表統計量を算出する。次に、参照発信源が通過していないグリッドについては、クリギングなどの空間補間手法を利用して空間分布を合成する（S１３）。これにより、対象領域の全グリッドについて、それぞれのグリッドから参照発信源が信号を送信した場合の電波センサの測定値の期待値を取得する。この処理を各電波センサについて行うことで、期待値の空間分布を算出し、出力する（S１４）。

次に、図１３を利用して位置推定フェーズの処理を説明する。まず、学習フェーズで作成した各センサの期待値の空間分布を取得する（S２１）。次に、電波センサが測定した電波特徴量測定値を取得する（S２２）。これは、リアルタイムで測定している場合には順次電波センサの測定値を取り込んで位置推定してもよく、そうでない場合は過去の測定値を時系列で順番に取り込んで位置推定しても良い。そして、同時刻に測定された複数センサの測定値に対して、時間ごとに以下の位置推定処理を行う。まず期待値の空間分布と電波センサの測定値から、解析の対象領域の全グリッドについて尤度を計算する（S２３）。そして、その尤度分布を基に、発信源位置を推定する（S２４）。この発信源位置には、対象時刻における尤度が最大値となる発信源位置と、尤度が極大値になる位置のうち尤度の最大値から設定した閾値の範囲内の尤度で、かつ設定した時間の範囲内の過去において尤度が最大値になった準発信源位置と、それ以外の尤度の最大値から設定した閾値の範囲内の尤度を取る尤度の極大値の位置である参考発信源位置を含む。そして、算出した尤度分布と、推定した発信源位置を表示装置に表示する（S２５）。この時、推定する対象時刻の推定位置だけでなく、過去の推定位置も同時に表示しても良い。

一般的な手法では尤度が最大値となる位置を推定値として出力していたが、本開示では同時に他の極大値の位置を算出し、過去に尤度が最大となった位置の情報を参照して、指定した期間内に尤度が最大となった極大値を準発信源位置として出力する。さらに、その他の、極大値の位置を参考発信源位置として出力する。これが効果をもたらす事例を説明する。ある時刻ｔ₁において推定対象の発信源が図１１のグリッド２から電波を送信し、その電波を測定して推定された尤度分布が図１１に示される尤度分布であるとする。この時、グリッド２の尤度が全推定領域で最大であるため位置推定システムはグリッド２を発信源位置と推定し、正しく推定できている。その次に位置推定をする時刻ｔ_２において、推定対象の発信源は同じ位置にとどまっていて電波の送信を継続しているが、発信源やセンサの周囲の車両や人などの反射物・遮蔽物が移動して伝搬環境が変動したため、電波センサの測定値が微小に変動し、図１１に示す尤度分布が微小に変動して図１４に示す尤度分布が出力されたとする。この時、全体の尤度分布は図１１とほぼ同じであるが、グリッド１の尤度がグリッド２の尤度よりも大きくなって最大となり、グリッド１が発信源位置と推定されている。一般的な手法であれば尤度や適合度が最大になる位置だけを出力するので、発信源位置がグリッド２からグリッド１に大きく移動したと推定される。この場合、本当に推定対象の発信源が移動したのか、それとも誤差が大きいか、の判断が難しい。また、グリッド２に免許を取得した既知の無線局がある場合に、図１４の尤度分布を出力する時刻ｔ_２の場合にはグリッド１に別の無線局があると判断されてしまう。このような事例は、先行文献に示される技術のうち、非特許文献２に示される連続的な期待値分布を使用する場合よりも、特許文献１や特許文献２に示される電波センサの測定値の期待値を不連続に設定する場合の方が多く発生する現象である。一方、本開示の手法であれば、グリッド２にも発信源がいる可能性があることが分析できる。そして、発信源位置と準発信源位置が大きく離れていることから、発信源が高速で移動した可能性よりも位置推定誤差が大きい事例である可能性が高いことが分析でき、グリッド１とグリッド２のそれぞれに発信源が存在する場合を想定して無線局の分析をすることができる。

[実施例１]
実際の例を用いて、第１の実施の形態の電波発信源位置推定システムの実施例を説明する。

まず、図１５に示すようにセンサ１～４を配置する。センサの(x,y)座標は、単位をｍとし、センサ１（-5270, 950）、センサ２（-4000, 450）、センサ３（-5400, -800）、センサ４（-3800, -900）である。各センサではアンテナで受信した電波の電力の測定を行い、測定値をサーバーにLTEで転送する。

次に学習フェーズとして、参照発信源を車両に搭載して、送信電力を30ｄBに設定して送信を開始し、図１５に示した経路を走行し、その時の位置情報をGPSで測定する。そして、その時の電波の受信電力を全センサで測定する。次に、センサの測定値と参照発信源の位置情報を、時刻が同じデータ同士を結合することにより、学習データを生成する。

次に、空間分布の合成を行う。まず、対象領域とグリッドを設定する。全センサの中でのｘ、ｙ座標それぞれの最小値と最大値を算出し、対象領域の最小座標地点と最大座標地点を特定する。本実施例の場合、最小値は(-5400, -900)、最大値は（-3800, 950）である。そして、その地点の外側に余白として2000mをとって対象領域とする。本実施例の場合、対象領域の最小座標地点は(-7400, -2900)、最大座標地点は（-800, 2950）である。そして、この対象領域を20m長のグリッドで分割する。

次に、センサ毎に、学習データを有する全グリッドについて、そのグリッド内のデータの平均値を算出してそのグリッドの代表統計量を算出する。これが、そのグリッドから参照発信源が送信した信号をそのセンサが測定した時の、測定値の期待値である。そして、この期待値を用いて学習データが存在しないグリッドの期待値を普通クリギングで空間補間する。このようにしてセンサ毎に合成した期待値の空間分布を図１６に示す。期待値の大小を色の濃淡で表示しており、明るい色がその地点から送信された電波の受信強度が大きいことを表し、暗い色がその地点から送信された電波の受信強度が小さいことを表す。ここで、センサ３の期待値はセンサ３を中心に等方的に分布している。これは、対象領域がセンサ３から見通しが良い場合や、遮蔽物が等方的に分布する場合などによく見られる現象である。一方、センサ４の分布は、異方的であり、かつセンサ４に近い領域が強度が弱く、少し離れた領域の強度が最も強い。これは、センサから特定の方向にだけ電波の遮蔽物がある場合や、センサがビルの屋上などにあってアンテナの指向性がビルの足元の方向が弱い場合などによく見られる現象である。このようにして期待値の空間分布を算出し、学習フェーズを終了する。

次に、位置推定フェーズに入る。まず、全電波センサで位置推定対象からの電波の受信強度を測定する。ある時刻ｔ₁において、センサ１～４のセンサの受信強度が、p₁=-89.9dB, p₂=-77.0dB, p₃=-85.4dB, p₄=-82.9dBであった。

この測定値を用いて、電力・尤度分布を算出する。ΔPを-20～20の範囲で5dB刻みに設定し、センサ毎に数４を用いて各ΔPにおける全グリッドの尤度

を算出する。次に、数５を用いて全センサの尤度を結合した結合尤度分布L(x,y,ΔP)を算出する。これが、電力・尤度分布である。

次に、位置推定部５６において送信電力の推定と尤度分布の推定を行う。まず、L(x,y,ΔP)の中で尤度が最大になるΔPを検索すると、ΔP=5のときに尤度が最大となった。ΔP=5のときの各センサの尤度分布

を図１７に、結合尤度分布L(x,y,ΔP=5)を図１８に示す。この尤度分布において、尤度が最大であるグリッドＡをこの時刻ｔ₁における発信源位置として推定する。次に、この尤度分布における尤度の極大値を算出し、数５で得られる結合対数尤度の値で、最大値との差分が設定した閾値である-10の範囲内の極大値を抽出すると、グリッドＢのみが抽出された。このグリッドＢを参考発信源位置として出力する。このようにして、本開示の電波発信源位置推定システムにより、発信源の位置を推定できる。尚、この時の実際の発信源位置は図１８に示されるようにグリッドＡの近傍であったとする。

次に、時刻ｔ₁から1秒後の時刻ｔ_２において、同様にして位置推定を行う。発信源の位置は図１９に示すようにほぼ変わらずグリッドAの近傍であるが、周囲の車両などが移動したため、時刻ｔ_２におけるセンサ１から４の受信強度が、p₁=-88.5dB, p₂=-77.2dB, p₃=-86.3dB, p₄=-83.3dBと若干変化した。この測定値を用いて推定した尤度分布を図１９に示す。図１８の尤度分布と類似しているが、各グリッドの尤度の値は若干変化している。まず、尤度が最大になるグリッドはグリッドBに変化し、尤度の極大値は数か所で算出され、そのうち尤度の最大値から-10以内の尤度Lを有する極大値のグリッドはグリッドAだけであった。この場合、まず尤度が最大になったグリッドBを発信源位置として出力する。そして、位置推定の対象時刻ｔ_２から設定した時間である10秒過去に遡った期間の発信源位置の情報を記憶部６４から取り込んで検出した極大値のグリッドと比較すると、グリッドAが設定した期間内のｔ₁において尤度が最大値になり発信源位置として推定されているので、グリッドAを準発信源位置として出力する。

一般的な手法ではグリッドBだけを発信源位置として出力するので発信源が時刻ｔ₁から時刻ｔ_２にかけて大きく移動したように見えるが、本開示では実際の発信源位置に近いグリッドAも準発信源位置として出力するので、グリッドAに発信源が存在する可能性があると判断でき、グリッドAとグリッドBの両方を発信源位置の候補として未知無線局の探索を行うことができる。尚、推定結果に対して移動平均のような時系列フィルタ処理を施すと、推定位置はグリッドAとグリッドBの中間に推定されて実際の位置からの誤差が大きく判断されるため、有効ではない。

このように、本開示は、推定した電力における尤度分布を表示し、かつ設定した過去の期間内に推定した発信源位置の情報を参照して準発信源位置を表示する。この手法が、特許文献１や特許文献２のように電波センサの測定値の期待値を都市の建物や物体やブレークポイントに対応して不連続に設定する場合に特に有効であることを、非特許文献２の方法に適用した場合を図示することによって説明する。

図２０は、図１５の学習フェーズで得られた受信電力の測定値と、その参照発信源とセンサとの距離との関係をプロットした図である。最小二乗法で数１にフィッティングすることにより、数１の伝搬定数(α,β)が得られる。この伝搬定数を用いて図１６に相当する期待値分布を算出すると、図２１のようにセンサを中心に等法的・連続的に減少する分布が得られる。これは、図１６の不連続な分布とは対照的である。そして、本実施例の時刻ｔ₁における測定値と数５を用いて図１７に相当する各電波センサの尤度分布を算出すると、図２２に示す分布が得られる。そして、結合尤度分布は図２３のようになり、尤度が最大になったグリッドCを発信源位置と推定する。

この手法では、電波を遮蔽する都市の建物やブレークポイントを考慮しない為、実際の発信源位置と離れた位置が発信源位置と推定されており、位置推定誤差が大きい。その一方で、尤度の算出に用いた期待値の分布が図２１に示すように連続的なので、最終的に算出する図２３の結合尤度分布も連続的になり、測定値が変動しても推定する発信源位置が不連続に大きく変動することはない。このように、本開示の手法は、測定値の期待値を都市の建物や物体やブレークポイントに対応して不連続に設定する場合に特に発生する推定位置の大きな変動に対して、推定対象の発信源位置の分析を可能にする推定結果を出力する。

以上、実施の形態を参照して本願発明を説明したが、本願発明は上記によって限定されるものではない。本願発明の構成や詳細には、発明のスコープ内で当業者が理解し得る様々な変更をすることができる。

本開示は、違法電波源の位置推定システムや、自動運転の位置推定システム、遭難者の位置推定システム、といった用途に適用できる。

１０ 参照発信源位置測定装置
２０ センサ
３０ 解析サーバー
４０ 記憶装置
５０ 位置推定装置
５１ 前処理部
５２ 期待値生成部
５３ 学習データ生成部
５４ 空間分布合成部
５５ 尤度算出部
５６ 位置推定部
５７ 表示部
６１ 時系列統計処理部
６２ 電力推定部
６３ 発信源位置推定部
６４ 記憶部
６５ 極大値算出部
６６ 準発信源位置抽出部

Claims (14)

1. 既知位置の参照発信源から送信された電波の電波特徴量測定値に基づき該既知位置から電波が送信された場合の電波特徴量測定値の期待値を取得する学習データ生成部と、
   該既知位置の期待値を合成することにより対象エリアの任意の位置から電波が送信された場合の測定値の期待値を算出する空間分布合成部と、
   該期待値と電波センサで測定した受信電力との適合度合いを基に対象エリアの各地点に電波源が存在する適合度の分布を計算する適合度算出部と、
   適合度算出部が算出した適合度分布をもとに適合度が最大値や極大値になる位置を算出して発信源の位置を推定する位置推定部と、
   対象エリアにおける適合度の空間分布と発信源の推定位置を表示する表示部と、
   を備え、
   前記位置推定部において、指定した範囲内の期間において発信源位置と推定された位置の近傍の極大値を準発信源位置として出力し、それ以外の極大値を参考発信源位置として出力し、
   前記表示部はそれら発信源位置と準発信源位置と参考発信源位置を区別して表示することを特徴とする、電波発信源位置推定システム。
2. 前記電波特徴量は、受信電力であることを特徴とする、請求項１に記載の電波発信源位置推定システム。
3. 前記適合度は、受信強度の確率密度関数を基にした尤度であり、前記適合度の分布は尤度分布であり、前記適合度算出部は尤度算出部であることを特徴とする、請求項１または２に記載の電波発信源位置推定システム。
4. 前記位置推定部は、尤度が最大になる送信電力を推定し、その送信電力における尤度分布を算出し、その尤度分布において尤度が最大値と極大値になる位置を算出し、
   前記表示部は、推定した送信電力における尤度分布を表示することを特徴とする、請求項１から３のいずれか１項に記載の電波発信源位置推定システム。
5. 前記表示部は、算出した尤度の最大値の位置と極大値の位置を発信源位置候補として表示することを特徴とする、請求項１から４のいずれか１項に記載の電波発信源位置推定システム。
6. 前記位置推定部において、算出した尤度の極大値と尤度の最大値との差を計算し、その差が設定した閾値の範囲内の極大値だけを出力することを特徴とする、請求項５に記載の電波発信源位置推定システム。
7. 前記表示部は、指定した範囲内の期間において尤度が最大値になった場所を表示することを特徴とする、請求項１から６のいずれか１項に記載の電波発信源位置推定システム。
8. 既知位置の参照発信源から送信された電波の電波特徴量測定値に基づき該既知位置から電波が送信された場合の電波特徴量測定値の期待値を取得し、
   該既知位置の期待値を合成することにより対象エリアの任意の位置から電波が送信された場合の測定値の期待値を算出し、
   該期待値と電波センサで測定した受信電力との適合度合いを基に対象エリアの各地点に電波源が存在する適合度の分布を計算し、
   適合度算出部が算出した適合度分布をもとに適合度が最大値や極大値になる位置を算出して発信源の位置を推定し、
   対象エリアにおける適合度の空間分布と発信源の推定位置を表示し、
   前記発信源の位置を推定する際に、指定した範囲内の期間において発信源位置と推定された位置の近傍の極大値を準発信源位置として出力し、それ以外の極大値を参考発信源位置として出力し、
   前記空間分布と前記発信源の推定位置を表示する際に、それら発信源位置と準発信源位置と参考発信源位置を区別して表示することを特徴とする、
   電波発信源位置推定方法。
9. 前記電波特徴量は、受信電力であることを特徴とする、請求項８に記載の電波発信源位置推定方法。
10. 前記適合度は、受信強度の確率密度関数を基にした尤度であり、前記適合度の分布は尤度分布であることを特徴とする、請求項８または９に記載の電波発信源位置推定方法。
11. 前記発信源の位置を推定する際に、尤度が最大になる送信電力を推定し、その送信電力における尤度分布を算出し、その尤度分布において尤度が最大値と極大値になる位置を算出し、
    前記空間分布と前記発信源の推定位置を表示する際に、推定した送信電力における尤度分布を表示することを特徴とする、請求項８から１０のいずれか１項に記載の電波発信源位置推定方法。
12. 前記空間分布と前記発信源の推定位置を表示する際には、算出した尤度の最大値の位置と極大値の位置を発信源位置候補として表示することを特徴とする、請求項８から１１のいずれか１項に記載の電波発信源位置推定方法。
13. 前記発信源の位置を推定する際に、算出した尤度の極大値と尤度の最大値との差を計算し、その差が設定した閾値の範囲内の極大値だけを出力することを特徴とする、請求項１２に記載の電波発信源位置推定方法。
14. 前記空間分布と前記発信源の推定位置を表示する際に、指定した範囲内の期間において尤度が最大値になった場所を表示することを特徴とする、請求項８から１３のいずれか１項に記載の電波発信源位置推定方法。

Images