JP6509685B2 - 移動速度推定装置、方法、及びプログラム - Google Patents

Description

本発明は、移動速度推定装置、方法、及びプログラムに係り、特に、未観測地点又は未観測時点における移動物体の速度を推定する移動速度推定装置、方法、及びプログラムに関する。

２次元平面上を移動する多数の物体（人、車など）の緯度経度情報を一定時刻間隔である一定時間観測された多数の移動物体の時空間情報から、当該時間内における観測されていない場所や時刻での速度を推定する問題は、人流、交通流の解析のための要素技術として重要である。

この問題に対し、時空間ボクセルデータにおいて、観測された時空間ボクセルから観測されていない時空間ボクセルでの移動物体群の平均速度を推定する方法として、時空間ボクセルデータを行列もしくはテンソルデータとみなし、行列、テンソル分解に基づく手法が提案されている（例えば、非特許文献１、２）。

Wu-Jun and Yeung, Dit-Yan: Relation regularized matrix factorization, International Conference on Artificial Intelligence (IJCAI), 2009. Wang, Yilun and Zheng, Yu and Xue, Yexiang: Travel time estimation of a path using sparse trajectories, International Conrence of Kowledge Discovery and Data Mining (ACM SIGKDD), 2014.

しかし、上記の非特許文献１、２に記載の、行列、テンソル分解に基づく手法では、近傍のボクセルでの速度は類似しているという情報を用いて定式化しているが、速度場の時空間相関までは手法に反映させていないため推定精度に限界があった。さらに、時空間外挿入、すなわち、あるボクセルでの将来時刻での速度推定は原理的に推定不可能であった。

本発明は、上記問題点を鑑みてなされたものであり、未観測地点又は未観測時点における移動物体の速度を精度よく推定することができる移動速度推定装置、及びプログラムを提供することを目的とする。

上記目的を達成するために、本発明に係る移動速度推定装置は、一定時刻間隔で観測された移動物体の位置情報及び時刻情報を含む２次元位置情報データを入力として、前記２次元位置情報データが含まれる２次元平面をメッシュに分割し、かつ、時間方向に量子化したボクセルを生成する時空間量子化部と、

前記ボクセルの各々に対し、前記ボクセルに存在する移動物体の平均速度を算出し、時空間ボクセルデータを生成する平均速度算出部と、

前記時空間ボクセルデータに基づいて、前記ボクセルの各々の平均速度の重み付き線形結合により、未観測のボクセルにおける前記移動物体の速度を推定する速度推定部と、を含んで構成されている。

本発明に係る移動速度推定方法は、時空間量子化部、平均速度算出部、及び速度推定部を含む移動速度推定装置における移動速度推定方法であって、

前記時空間量子化部が、一定時刻間隔で観測された移動物体の位置情報及び時刻情報を含む２次元位置情報データを入力として、前記２次元位置情報データが含まれる２次元平面をメッシュに分割し、かつ、時間方向に量子化したボクセルを生成し、

前記平均速度算出部が、前記ボクセルの各々に対し、前記ボクセルに存在する移動物体の平均速度を算出し、時空間ボクセルデータを生成し、

前記速度推定部が、前記時空間ボクセルデータに基づいて、前記ボクセルの各々の平均速度の重み付き線形結合により、未観測のボクセルにおける前記移動物体の速度を推定する。

本発明に係るプログラムは、コンピュータを、上記の移動速度推定装置の各部として機能させるためのプログラムである。

本発明の移動速度推定装置、方法、及びプログラムによれば、２次元平面をメッシュに分割し、かつ、時間方向に量子化したボクセルを生成し、ボクセルの各々に対し、ボクセルに存在する移動物体の平均速度を算出して、時空間ボクセルデータを生成し、ボクセルの各々の平均速度の重み付き線形結合により、未観測のボクセルにおける移動物体の速度を推定することにより、未観測地点又は未観測時点における移動物体の速度を精度よく推定することができる、という効果が得られる。

本発明の実施の形態に係る移動速度推定装置の構成を示すブロック図である。 （Ａ）処理対象の２次元平面を示す図、及び（Ｂ）２次元平面を矩形メッシュに分割した様子を示す図である。 時間的量子化を説明するための図である。 （Ａ）タイムスタンプｔでの移動物体の観測データの例を示す図、（Ｂ）タイムスタンプｔでの北方向の平均速度を算出する方法を説明するための図、及び（Ｃ）時空間ボクセルデータの例を示す図である。 本発明の実施の形態に係る移動速度推定装置における移動速度推定処理ルーチンを示すフローチャートである。 本発明の実施の形態に係る手法の実験結果である。

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

＜システム構成＞
本発明の実施の形態に係る移動速度推定装置の構成について説明する。図１に示すように、本実施の形態に係る移動速度推定装置１００は、ＣＰＵと、ＲＡＭと、後述する移動速度推定処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この移動速度推定装置１００は、機能的には図１に示すように入力部１０と、演算部２０と、出力部４０とを備えている。

入力部１０は、一定時刻間隔で観測された各移動物体の位置情報及び時刻情報を含む２次元位置情報データを受け付ける。また、入力部１０は、推定対象となる未観測の地点及び時刻を受け付ける。

各移動物体はＧＰＳなどの位置センサにより、ある一定時刻間隔で位置情報が観測されるものとする。ただし、ＧＰＳなどのセンサで位置を計測する際、ある時刻で計測が失敗すると、その時刻の前後で移動軌跡は分断され得るので、ある移動物体の観測軌跡は必ずしも連続しているとは限らない。本実施の形態に係る移動速度推定装置１００はその様な観測データでも適用可能である。

演算部２０は、観測データ２２、時空間量子化部２４、平均速度算出部２６、時空間ボクセルデータ２８、速度推定部３０とを含んで構成されている。

観測データ２２は、入力部１０で受け付けた一定時刻間隔で観測された各移動物体の２次元位置情報データを記憶している。

＜時空間量子化部＞
時空間量子化部２４は、前記２次元位置情報データが含まれる２次元平面をメッシュに分割し、かつ、時間方向に量子化したボクセルを生成する。

本実施の形態では、位置センサの誤差や、人流、交通流での応用を勘案し、ある地点（２次元平面上の一点）、ある時点（ある時刻）での移動物体の速度推定ではなく、ある領域、ある時間帯での平均的な速度情報を推定する。そのために、観測データを空間（平面）的、時間的に量子化する必要がある。

２次元平面をメッシュに分割する空間的量子化については、図２（Ａ）に示すように、処理対象となる２次元平面を、経度方向および緯度方向にある適切な距離幅で分割することで、図２（Ｂ）に示すように、２次元平面を矩形メッシュに分割する。この距離幅は、二次元平面上でどれくらいの領域サイズで移動物体の速度を推定するかという実要求に依存して決定される。この空間的量子化を行う処理により、処理対象なる２次元平面がＩ×Ｊのメッシュに分割される。図２（Ｂ）に、Ｉ＝８、Ｊ＝１２の例を示す。

次いで、時間方向の量子化について説明する。今、図３に示すように、３つの移動物体が、ある時間内において一定時刻間隔で観測されているとする。一般にＧＰＳなどの位置センサで移動物体を観測する場合、移動物体毎のサンプリング時刻が必ずしも同期しているとは限らない。したがって、ある時刻で観測される移動物体の個数が少ないと、その時刻での未観測地点での移動物体の推定が困難となる。そこで、各移動物体の時刻情報をある一定時間幅で量子化する。図３では、３つの移動物体の時刻情報(ｔ1，ｔ2，ｔ3)が、タイムスタンプｔ＝｛１，２，３，４｝に量子化される。この時間的量子化により、例えば、タイムスタンプｔ＝1では、計６点の位置情報が観測されたことになる。

＜平均速度算出部＞
平均速度算出部２６は、ボクセルの各々に対し、当該ボクセルに存在する移動物体の８方向の平均速度を算出し、時空間ボクセルデータを生成する。

ここで、時空間量子化により変換されたデータから、各メッシュでの各タイムスタンプについて、８方向の平均速度を算出する方法について説明する。タイムスタンプ１，・・・・，Ｔに渡り、少なくとも一つ移動物体が存在するメッシュ領域（ｉ，ｊ）（ｉ∈｛１，・・・，Ｉ｝，ｊ∈｛１，・・・，Ｊ｝）の集合をＤ_sで表すものとする。便宜上、そのようなメッシュ領域がＭ個存在するとし、ｓ₁，ｓ₂，・・・，ｓ_M∈Ｄ_sと表すものとする。また、ｓ_mには合計Ｔ_m個のタイムスタンプで移動物体が存在していたとし、ｔ_mkをＴ_mにおける第ｋ番目のタイムスタンプとする。図３に示した様に、時間方向での移動物体の観測数は明らかにｓ_m毎に異なり得るので、ｔ_mkをｍ依存の表記としている。明らかに、ｔ_mk∈｛１，・・・，Ｔ｝である。

図４（Ａ）に、Ｉ＝Ｊ＝３で空間量子化された、３×３のメッシュにおけるタイムスタンプｔでの移動物体の観測データの例を示す。図４（Ａ）の２点間を結ぶ矢印は、その時間窓(タイムスタンプt)での移動ベクトルを表している。各移動物体の移動始点(矢印の起点)での方向を東(E)、西(W)、南(S)、北(N)、北東(NE)、北西(NW)、南西(SW)、南東(SE)の８方向に量子化する。具体的には、各移動ベクトルに対し、８方向の単位ベクトル（ｅ_E＝（１，０）’，ｅ_W＝（−１，０）’，ｅ_S＝（０，１）’，ｅ_N＝（０，１）’，ｅ_NE＝（√２／２，√２／２）’，ｅ_NW＝（−√２／２，√２／２）’，ｅ_SW＝（−√２／２，−√２／２）’，ｅ_SE＝（√２／２，−√２／２）’）との内積演算を行い、その値が最大となる方向を、分類方向とすれば良い。図４（Ａ）に移動方向の８方向量子化の例を示す。なお、あるメッシュ内でタイムスタンプｔでの移動物体の観測データが１つだけしか得られなかった場合には、当該移動物体について得られた次の観測データとの差分から移動ベクトルを求めればよい。

次いで、各メッシュ領域において以下の処理を行う。今、メッシュ領域ｓ_m、タイムスタンプｔに着目すると、図４（Ｂ）では北(Ｎ)方向に速度ｖ₁，ｖ₂，ｖ₃で移動する３つの移動物体が存在するので、メッシュ領域ｓ_m、タイムスタンプｔでの北方向の平均速度はｖ₁，ｖ₂，ｖ₃の算術平均として、

とする。この処理をすべてのＤ_s、ｔ_mk∈｛１，・・・、Ｔ_m｝に対し、各８方向独立に移動物体の平均速度を算出することで、

を得る。なお、あるメッシュ内でタイムスタンプｔでの移動物体の観測データが１つだけしか得られなかった場合には、当該移動物体について得られた次の観測データとの差分から速度ｖを求めればよい。以上の処理を実行することにより、図４（Ｃ）に示す様に、サイズＩ×Ｊ×Ｔの時空間ボクセルデータＸが生成できる。ただし、平均速度算出部２６では、時空間ボクセルデータのＩ×Ｊ×Ｔの全ボクセルにおいて平均速度が算出されるわけではなく、ｓ_m，ｔ_mk（ｋ＝１，・・・，Ｔ_m，ｍ＝１，・・・，Ｍ）に対応するボクセルでの平均速度が算出されたことになる。

時空間ボクセルデータ２８には、平均速度算出部２６で生成された時空間ボクセルデータが記憶される。

＜速度推定部＞
速度推定部３０は、時空間ボクセルデータ２８に基づいて、時空間クリギング手法を用いて、ボクセルの各々の平均速度の重み付き線形結合により、未観測のボクセルにおける移動物体の速度を推定する。

ここで、時空間ボクセルデータから時空間クリギング手法（非特許文献３）を用いてｓ_m∈Ｄ_sにおける未観測のボクセルでの速度を推定する方法について説明する。ただし、推定法は方向に依存しないため、以下ではある一方向での速度の推定法について説明することとし、実際には、以下の処理を各８方向について逐次実施すれば良い。

［非特許文献３］：Cressie, Noel and Wikle K. Christopher: Statistics for spatio-temporal data, Wiley, 2011.

Ｖ（ｓ₀；ｔ₀）を時空間ボクセルデータでの、未観測地点又は未観測時点に対応するボクセルでの平均速度を表すものとすると、クリギング法では次式に示す様に、予測値＾Ｖ（ｓ₀；ｔ₀）を観測データの重み付き線形結合として推定する。

ただし、ｗは推定すべき重みベクトルであり、未観測のボクセル毎に推定される。~ｖは値が既知の観測値ベクトルである。表記を簡単にするため、

とする。ただし、記法‘´’は転置演算を表す。Ｌ＝Ｔ₁＋・・・＋Ｔ_M とすると、ｗおよび~ｖの次元はＬとなる。さらに、~ｖ（ｉ）は~ｖの第ｉ要素を表すものとする。

重みベクトルは予め推定されているものとする。重みベクトルの推定では、２次定常性を仮定する。すなわち、以下の式を用いる。

ここで、μは未知の平均値（定数）とし、ｈ＝（ｈ_s，ｈ_t）はラグベクトルを表す。ｈ_sは空間ラグ（ｓ_mとｓ_m′間の距離）を、ｈ_tは時間ラグ(タイムスタンプｔとｔ′の差)を各々表す。式(6)はＶ（ｓ_m；ｔ）とＶ（ｓ_m′；ｔ′）の共分散がｈ_sおよびｈ_tのみに依存し、位置や時刻に依らないことを意味する。

式(4)および式(5)より、

を得る。ここで、

はＬ次元の１ベクトルとする。

Ｖ（ｓ₀；ｔ₀）をｓ₀およびｔでの真値(未知)とすると、期待自乗予測誤差は、

となる。式(7)の制約条件の下で期待自乗予測誤差をｗに関して最小化することにより、重みの最適解は次式のように求まる。

ここで、Ｃ^-1はコバリオグラム行列の逆行列を表し、コバリオグラム行列の第(ｉ，ｊ)要素は、

に相当する。また、ｃ₀はＬ次元の共分散ベクトルを表す。

Ｖ（ｓ₀；ｔ₀）は未知ゆえ、ｃ₀は直接計算できないが、Ｃ同様、二次定常性の性質を用いて推定可能である。上記のように、時空間ボクセルデータ２８に基づいて、未観測のボクセル毎に、重みベクトルｗが推定される。

出力部４０は、速度推定部３０によって推定された、未観測のボクセルにおける移動物体の速度を結果として出力する。

＜本発明の実施の形態に係る移動速度推定装置の作用＞
次に、本発明の実施の形態に係る移動速度推定装置１００の作用について説明する。移動速度推定装置１００は、入力部１０において一定時刻間隔で観測された各移動物体の２次元位置情報データを受け付けると、観測データ２２に格納する。そして、入力部１０において、推定対象となる未観測の地点及び時刻を受け付けると、移動速度推定装置１００は、図５に示す移動速度推定処理ルーチンを実行する。

まず、ステップＳ１００において、時空間量子化部２４は、観測データ２２から各移動物体の２次元位置情報データを読み込み、２次元位置情報データが含まれる２次元平面をメッシュに分割し、かつ、時間方向に量子化したボクセルを生成する。

次に、ステップＳ１０２において、平均速度算出部２６は、ボクセルの各々に対し、当該ボクセルに存在する移動物体の８方向の各々の平均速度を算出し、時空間ボクセルデータを生成し、時空間ボクセルデータ２８に格納する。

ステップＳ１０４において、速度推定部３０は、上記ステップＳ１０２で生成された時空間ボクセルデータに基づいて、時空間クリギング手法を用いて、ボクセルの各々の平均速度の重み付き線形結合により、推定対象となる未観測の地点及び時刻に対応する未観測のボクセルにおける移動物体の速度を推定し、出力部４０により出力し、移動速度推定処理ルーチンを終了する。ここで、推定対象となる未観測の地点及び時刻に対応する未観測のボクセルの各々について、ボクセルの各々に対する重みは、観測データ２２から読み込んだ各移動物体の２次元位置情報データに基づいて求められる。
また、未観測のボクセルとして、内挿により推定されるものと、外挿により推定されるものとがある場合には、内挿により推定される未観測のボクセルにおける移動物体の速度を推定した後に、外挿により推定される未観測のボクセルにおける移動物体の速度を推定する。

＜実施例＞
本発明の実施の形態に係る手法の有効性を確認すべく、図６に、ローマ市内でのタクシーの移動軌跡データを用いた実験結果を示す。本データは、２０１４年２月４日の午前１０時から１０時４０分までのローマにおける３１２台のタクシーのＧＰＳにより計測された移動軌跡のデータである。観測領域は、南北に２．８４ｋｍ、東西に５．０３ｋｍで、メッシュサイズは１０ｍ四方とした。本データに対し、１０時から１０時３０分までのタイムスタンプのデータを用いて、前述した時空間ボクセルデータを作成した。そして、本発明の実施の形態に係る手法、および、従来法（非特許文献１）を適用して、推定誤差比較を行った。実験評価では、10-fold cross-validation で評価し、誤差の評価尺度として、平均相対誤差：

を用いた。ここで、^ｖ_iは第ｉテストデータに対する推定値で、ｖ_iは真値を表す。Ｎはテストデータ総数である。

図６の内挿とは、１０時から１０時３０分までの観測データから生成した時空間ボクセルデータを用いて10-fold cross-validation により評価した結果に相当する。そして、図６の外挿とは、１０時から１０時４０分までの観測データから生成した時空間ボクセルデータに対し、１０時から１０時３０分までに相当する時空間ボクセルデータのみから、１０時３０分から１０時４０分に相当する時空間ボクセルの速度を推定する問題である。明らかに、外挿は将来予測に相当し、内挿に比べより困難な問題である。行列・テンソル分解に基づく方法では、原理的に外挿問題は適用できないため、外挿問題に対しては、時系列解析手法で著名な従来法（ベクトル自己回帰モデル：VAR）と比較した。

図６の結果から明らかなように、本発明により内挿に関しては顕著に精度の高い推定が実現できている。これは、従来法が、時空間的に近いボクセル同士は類似した値を有するという仮定に基づいているのに対し、本発明の実施の形態に係る手法では、時空間相関というより正確な統計量を用いて推定していることに起因する。外挿問題では、内挿に比べ誤差が大きいことが確認できる。これは前述した通り、外挿問題が内挿問題に比べ問題の難易度が高いことを実証している。外挿問題では本発明の実施の形態に係る手法はＶＡＲとほぼ同等な結果である。ただし、ＶＡＲが将来予測のための手法であるのに対し、本発明の実施の形態に係る手法では、内挿と外挿を同じ枠組みで実現できるという特長を有する。以上の結果より、本発明の実施の形態に係る手法が移動物体の速度の時空間内挿、外挿が適切に推定でき、特に内挿に関しては従来法に対する顕著な優位性を確認した。

以上説明したように、本実施の形態に係る移動速度推定装置によれば、２次元平面をメッシュに分割し、かつ、時間方向に量子化したボクセルを生成し、ボクセルの各々に対し、ボクセルに存在する移動物体の平均速度を算出して、時空間ボクセルデータを生成し、時空間クリギング手法を用いて、ボクセルの各々の平均速度の重み付き線形結合により、未観測のボクセルにおける移動物体の速度を推定することにより、未観測地点又は未観測時点における移動物体の速度を精度よく推定することができる。また、未観測地点又は未観測時点での移動物体の速度が推定できれば、渋滞予測、さらには予測に基づく誘導対策などに活用できる。

なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、上記実施の形態では、移動物体の移動方向を、東、西、南、北、北東、北西、南西、及び南東の８方向の何れかに分類する場合を例に説明したが、これに限定されるものではなく、上、下、右、左、右上、左上、左下、及び右下の８方向の何れかに分類するようにしてもよい。

また、上述の移動速度推定装置は、内部にコンピュータシステムを有しているが、コンピュータシステムは、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。

１０ 入力部
２０ 演算部
２２ 観測データ
２４ 時空間量子化部
２６ 平均速度算出部
２８ 時空間ボクセルデータ
３０ 速度推定部
４０ 出力部
１００ 移動速度推定装置

Claims (6)

1. 移動物体の二次元平面上の位置情報と時刻情報からなる時系列データを観測データとして、
   前記観測データが含まれる２次元平面をメッシュに分割し、かつ、時間方向に量子化したボクセルを生成する時空間量子化部と、
   前記観測データに基づいて、前記ボクセルの各々に対し、当該ボクセル内に存在する移動物体の各々について前記観測データから求められる、前記移動物体の各々の移動速度の平均値を、当該ボクセルに存在する移動物体の平均速度として算出し、時空間ボクセルデータを生成する平均速度算出部と、
   前記時空間ボクセルデータに基づいて、前記ボクセルの各々の平均速度の重み付き線形結合により、未観測のボクセルにおける前記移動物体の速度を推定する速度推定部と、
   を含む移動速度推定装置。
2. 前記平均速度算出部は、前記観測データから、各移動物体の各タイムスタンプにおける移動ベクトルを求め、前記移動ベクトルを用いて前記ボクセル内に存在する移動物体の平均速度を算出することを特徴とする請求項１記載の移動速度推定装置。
3. 前記平均速度算出部は、移動物体の移動方向を、東、西、南、北、北東、北西、南西、及び南東の８方向の何れか、又は上、下、右、左、右上、左上、左下、及び右下の８方向の何れかに分類し、各分類方向ごとに、前記ボクセルに存在する移動物体の平均速度を算出することを特徴とする請求項１又は２記載の移動速度推定装置。
4. 前記速度推定部は、
   前記観測データに基づいて求められた前記ボクセルの各々に対する重みを用いた、前記ボクセルの各々の平均速度の重み付き線形結合により、未観測のボクセルにおける前記移動物体の速度を推定する請求項１〜３のいずれか１項記載の移動速度推定装置。
5. 時空間量子化部、平均速度算出部、及び速度推定部を含む移動速度推定装置における移動速度推定方法であって、
   移動物体の二次元平面上の位置情報と時刻情報からなる時系列データを観測データとして、
   前記時空間量子化部が、前記観測データが含まれる２次元平面をメッシュに分割し、かつ、時間方向に量子化したボクセルを生成し、
   前記平均速度算出部が、前記観測データに基づいて、前記ボクセルの各々に対し、当該ボクセル内に存在する移動物体の各々について前記観測データから求められる、前記移動物体の各々の移動速度の平均値を、当該ボクセルに存在する移動物体の平均速度として算出し、時空間ボクセルデータを生成し、
   前記速度推定部が、前記時空間ボクセルデータに基づいて、前記ボクセルの各々の平均速度の重み付き線形結合により、未観測のボクセルにおける前記移動物体の速度を推定する
   移動速度推定方法。
6. コンピュータを、請求項１〜請求項４のいずれか１項に記載の移動速度推定装置として機能させるためのプログラム。