隠れマルコフモデルを用いたマップマッチングの基礎論文をまとめてみる

マップマッチングの論文を調べていたら、隠れマルコフモデルを用いてマップマッチングをする原初の論文が公開されていたので、せっかくだから勉強してプログラムを組んでみました。
この手法はST-Matchingと言われる手法です。
元論文はしたのになりますので、気になったら読んでみてください（英語なので読むの大変ですが、比較的平易に書いてくださってるのでその点では読みやすいかなと）

今回のキーワード

• 隠れマルコフモデル(HMM)
• ビタビアルゴリズム

そもそもマップマッチングとは？

GPSデータからその人が移動したであろう経路を推測するための操作
下の記事が参考になったのでお勧めです

今回はその中でも確率的マップマッチングに該当します。
正直、上の記事に大事なことは全部書いてありますが、今回はさらに詳しく書いていこうと思います

Step0：使用データとライブラリについての解説

今回マップマッチングはpythonを用いて実装していきますが、その前段階としていくつかの知識が必要になります

元データ

今回は以下のようなGPSデータを対象としています。
mymapで適当に作った仮想のGPSデータです。青のマーカーが1分毎(とする)GPS座標の位置、オレンジの線が想定する経路です。
このオレンジ戦を何とかして計算していきます。

取得する道路データの範囲

取得する範囲は元データを囲うように定義した、以下の範囲です。
この範囲のデータをOpenStreetMapという道路のWikiPediaみたいなサイトから取得します

使用するライブラリ

今回はpythonで実装するため、以下のライブラリを使って実装しています

ライブラリ名	説明
geopandas	データ分析などをするときに(ほぼ)必須のライブラリであるpandasを地図データに対応させたライブラリ
osmnx	OpenStreetMapのデータから道路ネットワークのグラフデータ(後述)を作成して、こねるためのライブラリ
scikit-learn	機械学習でよく使うライブラリ。今回はコサイン類似度だけ使います

道路ネットワークについて

osmnxで道路のデータを扱うときはグラフデータという形式で扱います。
グラフデータは大きく以下の2つの要素で構成されています。

1. ノード：いわゆる交差点などに該当する道路の端点
2. エッジ：ノード同士をつなぐもの。イメージする道路はこっち

今回の対象範囲の道路ネットワークのグラフデータは以下のようになっています
白丸がノードで線がエッジになります

以上が、最低限の基礎知識です。他にも最短経路を求めるためのダイクストラ法や、図形のデータを扱うライブラリであるshaplyなどがありますが、無くてもなんとかなるので気になったら調べてみてください。

隠れマルコフモデルを用いたマップマッチングの手法について

今回の論文の手法(ST-Mathing)では大きく分けて以下の2つの分析を行い、合計3つの指標を作成・計算しています

空間分析

空間分析では、以下の2つの指標を作成します。

1. 観測確率
   観測確率はある点$p_i$に対して、候補点$c_i^j$がそれぞれどれくらいあり得るのかの確率を求めた指標です。数式としては以下になります
   $N(c_i^j)=\dfrac{1}{\sqrt{2}\sigma}e^{-\frac{(x_i^j-\mu)^2}{2\sigma^2}}$
   論文だと$\mu=0$,$\sigma=20[m]$で定義されいています
   これはGPSの誤差は点$p_i$と$c_i^j$の2点間の距離、$dist(c_i^j,p_i)$についての正規分布$N(\mu,\sigma^2)$として表すことができるからです
   →雑にいうと近いほうが指標の値が高くなります。下の画像でいうと青点がGPSとした時に、赤点と黄点だと赤点の方があり得るよねってことです
   

2. 伝達確率
   伝達確率はとある候補点$c_i^j$から$c_{i-1}^j$への経路がどれだけそれっぽいかを示す指標です。数式としては以下になります。
   $V(c_{i-1}^t→c_i^s) = \dfrac{dist(p_i,p_{i-1})}{w_{(i-1,t)→(i,s)}}$
   ここで、$w_{(i-1,t)→(i,s)}$は候補点$c_{i-1}^t$から$c_i^s$への最短経路、$dist(p_i,p_{i-1})$はGPSの点$p_{i-1}$,$p_i$のユークリッド距離です。
   つまり、GPSの2点間の直線距離と、候補点の最短経路の長さを比較する指標です。値が大きく慣ればなるほど良いです。ただし、最短経路が直線距離よりも短くなることはないため、最大値は1になるはずです。
   イメージ図としては下図の感じで、黒線(ユークリッド距離)と青線(最短経路の長さ)を比較する感じです。
   

時間分析

空間分析を行うことで、各GPSの座標から求められる候補点がそれぞれどれくらいあり得そうかという指標を計算することができます。ただ、それだけだと車の速度などの指標を使っていないため、高速道路とその下の道の区別ができないといった欠点があります。それをなんとかするために、論文では時間分析を行なっています。
時間分析は道路の各セグメントに設定されている制限速度のデータと、ある候補点$c_{i-1}$と$c_{i}$間の平均速度のコサイン類似度を計算しています。論文ではコサイン距離と記載されていますが、数式的にはコサイン類似度です。

コサイン類似度とは

コサイン類似度とは2つのベクトル、a,bがどれくらい同じ方向に向いているかの指標です。-1から1の間の値を取り、1に近づくほど類似度が高くなります。
コサイン類似度はscikit-learnのcosine_similarityで計算することが出来ます

参考ページ
https://mathlandscape.com/cos-similar/

各道路セグメントの制限速度は以下のような行列として、定義されます。eは道路セグメントのエッジが該当します。osmnxであればmaxspeedという名前であるはずです。ただし、あまりデータが無いのでnullのことも多いです。

\begin{bmatrix}
e_{1max} & e_{2max} & ・・・ & e_{nmax}\\
\end{bmatrix}

候補点$c_{i-1}$と$c_{i}$間の平均速度$v_{(i-1,t)→(i,s)}$は、下図の青線の長さを$p_{i-1}$,$p_{i}$の時間差分で割ることによって計算します。

\begin{bmatrix}
v_{(i-1,t)→(i,s)} & v_{(i-1,t)→(i,s)} & ・・・ & v_{(i-1,t)→(i,s)}\\
\end{bmatrix}

指標を用いた最適経路の決定

上の3つの指標の値を隠れマルコフモデルの記号出力確率分布として扱って解くことで最適な経路を計算していきます。解くためのアルゴリズムとしては、動的計画法を用いたビタビアルゴリズムを用いて計算しています。
ビタビアルゴリズムの詳細は以下のリンクをご参考に。

実際の例を用いたマップマッチング

上の式などで理解できた場合は不要ですが、私は初見で理解できなかったので実際の処理の流れを図示してみようと思います。

今回の対象データ

下図のGPSデータを用います。一番上の画像から一部抜粋したものです
時系列的には、1分ごとに左から右に移動しているとします。それぞれ、$GPS_1,GPS_2,GPS_3$と呼称します

注意
今回の例では細かい計算を一切行いません。
各確率はわかりやすいようにキリの良い数値を適当に当てはめている点に留意して下さい。

Step1候補点の算出

各GPSのデータから、一定範囲内(青円)に存在する道路に対して候補点(赤点)を算出します。
便宜上それぞれの候補点に番号を振っています。
候補点は$C_1^1,C_1^2,C_1^3・・・C_3^3$と呼称します。図の1,2,3・・・と番号を振っている物に対応します。

Step2観測確率の計算

GPSと候補点から観測確率を計算します。

$GPS_1$に対して、($C_1^1,C_1^2,C_1^3,C_1^4$)の観測確率をそれぞれ(0.4,0.2,0.1,0.3)とする
$GPS_2$に対して、($C_2^1$)の観測確率を(0.5)とする
$GPS_3$に対して、($C_3^1,C_3^2,C_3^3$)の観測確率をそれぞれ(0.1,0.3,0.6)とする

Step3遷移確率の計算

各候補点どうしでの遷移確率は以下とします

$GPS_1からGPS_2$への遷移確率

from/to	$c_2^1$
$c_1^1$	0.9
$c_1^2$	0.7
$c_1^3$	0.7
$c_1^4$	0.6

また、$GPS_2からGPS_3$への遷移確率

from/to	$c_3^1$	$c_3^2$	$c_3^3$
$c_2^1$	0.6	0.7	0.9

Step4時間確率の計算

各候補点どうしでの時間確率は以下とします
ぶっちゃけOpenStreetMapのデータだと最高速の値があまり無いので、時間確率であまり差は出ません。

$GPS_1からGPS_2$への時間確率

from/to	$c_2^1$
$c_1^1$	0.9
$c_1^2$	0.9
$c_1^3$	0.9
$c_1^4$	0.7

また、$GPS_2からGPS_3$への時間確率

from/to	$c_3^1$	$c_3^2$	$c_3^3$
$c_2^1$	0.8	0.8	0.9

求めた指標を隠れマルコフモデルとして扱い、ビタビアルゴリズムで解く

ここまでの結果を図にすると以下のような形になります

ビタビアルゴリズムを用いて解く場合は、各移動($GPS_1→GPS_2$,$GPS_2→GPS_3$)に分けてそれぞれで最適な水位を見つけていくことになります。

Step1初期値を決める

ビタビアルゴリズムでは下図のような形で計算していきます。
今回の初期値は$c_1^1〜c_1^4$の値になります。

Step2 1回目の計算

候補点$c_1^1〜c_1^4$から$c_2^1$に向かうスコアが一番高いもの(一番それっぽいもの)を求めていきます。
スコアの計算式は、これまでに求めた指標をかけて出発点の確率を足したものです。

例)$c_1^1$から$c_2^1$に向かうスコアの計算
$スコア=c_2^1の観測確率 * (c_1^1からc_2^1への遷移確率)* (c_1^1からc_2^1への時間確率) + c_1^1のスコア(初期値)$
→
$スコア=0.5 * 0.9 * 0.9 + 0.4 = 0.805$

同じように$c_1^2〜c_1^4$までスコアを計算すると以下のようになります。

from/to	$c_2^1$
$c_1^1$	0.805
$c_1^2$	0.515
$c_1^3$	0.415
$c_1^4$	0.51

この結果から、スコアの一番高いルート$c_1^1 → c_2^1$が$c_2^1$への確率が一番高い経路となります。今回の例では$c_2^1$しかありませんでしたが、複数の候補点がある場合はその候補点までのそれぞれの一番確率が高い経路とスコアを保持します。
GPSデータの最終1つ手前まではこの計算を繰り返します。

Step3 2回目の計算

次に$c_2^1$から$c_3^1 〜 c_3^3$に向かうスコアを計算します。
それぞれのスコアは以下のようになります。

$スコア=c_3^1の観測確率 * (c_2^1からc_3^1への遷移確率)* (c_2^1からc_3^1への時間確率) + c_2^1のスコア$
→
$スコア=0.1 * 0.6 * 0.8 + 0.805 = 0.853$

from/to	$c_3^1$	$c_3^2$	$c_3^3$
$c_2^1$	0.853	0.937	1.291

この点で最終なので、各終点からのスコアを比較して一番高いものが最終的な経路となります。今回は$c_3^3$が一番高いため、$c_2^1 → c_3^3$の経路を採択します。
step2の結果から$c_2^1$への経路としては、$c_1^1 →c_2^1$なので、$c_1^1 →c_2^1 → c_3^3$が今回のマッチング結果となります。

イメージとしてはこんな感じです。本当は道路の形状情報を保持しているので綺麗に総形になります。

終わりに

上のような形で隠れマルコフモデルを用いてマップマッチングを行うことが出来ます。
読むのは大変でしたが、一回理解して仕舞えばそこまで難しいアルゴリズムでは無いので皆さんも一回実装してみて下さい。
全然整理できていませんが、私が作成したコードを以下のリポジトリにおいてますので、参考までに。