この記事で学べること

・上左図のようなノイズを含む軌跡データから上右図のようなデータを得るための前処理手法とその流れ

キーワード: 異常値除去, ノイズフィルタリング, 滞留点検出

はじめに

こんにちは。NTTドコモの古山凌です。
業務では位置情報の分析を行っています。

本記事ではGNSS経路情報の簡単なEDA(探索的データ分析)を行いながら、そこに含まれるノイズを確認し、実際に前処理を行うまでの一例をご紹介させていただきます。

実際に行う操作としては異常値除去、ノイズフィルタリング、滞留点検出と、この記事を読んでいただいている皆様には既にご存知な部分も多いかと思いますが、

軌跡データをGoogle Earth Pro上で可視化することで、それぞれの前処理の効果を視覚的に楽しんでいただけるような記事にしたいと思います。

まずは軌跡データを可視化してみる

この記事ではKaggleのGoogle Smartphone Decimeter Challengeコンペのデータセットとして公開されているスマートフォンの軌跡データの一部を用います。

データの中身は以下のような位置情報の時系列データになっています。

自動車に乗せられたスマートフォンの位置情報のログ
約1000msおきに記録される
緯度(latDeg),経度(lngDeg)情報を持っている

この位置情報データをKMLファイルとして出力しGoogle Earth Proに読み込ませると、地図上に位置情報データをプロットした図を取得することができます。

赤い点は約１秒おきに記録された位置情報を示し、各点を黄色い直線で結ぶことによって、スマートフォンが通った(とされる)道を可視化しています。

前処理の方針を考える

さて、ここからは今回のデータセットにどのようなノイズが含まれているかを確認し、どのような前処理をどのような手順で行うかを考えていきましょう。

まずは、先ほど可視化した自動車の軌跡にどのようなノイズが含まれているかを考えてみます。
以下に考察をまとめましたのでご覧ください。

全体的に軌跡がジグザグしている
信号待ちなどの駐停車地点においてもポイントが散乱している
立体道路の下層部分を通る際に記録された位置情報が外れてしまっている

駐停車地点においてもポイントが複数の車線を跨って散乱していることから、今回のデータは少なくとも数m~10m程度の誤差があることがわかります。

そのため全体的に、歩道にはみでて走行しているように見えたり、何度も車線変更を行いながら走行しているように見えています。

また、立体道路の下層を走行する際に記録されたはずの位置情報は、
地形の影響を受け、正しい位置情報が記録できていないことがわかります。

さて、ここからは上に挙げた3点の問題を解決できるような前処理を行っていきたいと思います。

ノイズフィルタリング

まずは1点目に挙げた、「全体的に軌跡がジグザグしている」問題を解決します。
今回扱うような計測誤差を含む系列データからノイズを取り除くにはノイズフィルタリングが有効です。位置情報の系列データに対して用いられるノイズフィルタリングの一般的な手法として

Mean Filter (平均値フィルタ)
Median Filter (中央値フィルタ)
Kalman Filter (カルマンフィルタ)
Particle Filter (粒子フィルタ)

が挙げられます。今回のデータはサンプリングレートが1000msと低くなく、かつ外れ値の割合が少ないこと、また実装が簡単なことからMean Filter(平均値フィルタ)を採用します。

# Mean Filter (直近5ポイントの平均を算出)
df["latDeg_MeanFilter"] = df["latDeg"].rolling(5).mean()
df["lngDeg_MeanFilter"] = df["lngDeg"].rolling(5).mean()

Mean Filterを適用することで以下の結果が得られました。

Mean Filterを適用することで、全体的にジグザグだった軌跡がかなり滑らかになりました。走行中の車線を判別することが難しかった軌跡が、Mean Filter適用後では走行中の車線の判別が容易であることがわかります。

しかし一方で、Mean Filterを適用することにより立体道路下層部分の走行時に記録された外れ値が後方5step分のポイントに影響を及ぼしてしまい、計5箇所のポイントが道路外にプロットされてしまっていることがわかります。

よって、Mean Filterによって外れ値が周辺のポイントに影響を及ぼすこと避けるために、"異常値除去"→"Mean Filter"の順序で前処理を行うことにします。

なお、外れ値の影響を小さくするためにMean FilterではなくMedian Filterを選択することも有効です。しかしKalman Filterなどの別の選択肢を取る際にも、事前に異常値除去を行うことは有効ですので、以下、異常値検出についても参考にしていただけると良いと思います。

異常値検出

ここでは立体道路の下層部分を走行中に記録された外れ値の検出および除去を行いましょう。

まずは、今回の位置情報の系列データから各ポイントにおける移動速度[km/h]を算出し、さらに1秒前との移動速度の差[km/h]を算出しました。

1秒前との移動速度の差[km/h]を折れ線グラフにすると以下のようになります。また標準偏差を算出すると33.11になりました。このデータには信号待ちの時間などの実際には停車している時間のデータも含まれていることを考えると、数m~10m程度の測定誤差は速度に換算するとかなり大きいことがわかります。

さらにこのデータ(1秒前との移動速度の差)をヒストグラムにすると以下のような正規分布が得られました。

よって今回は平均値 ± 標準偏差 * 3を異常値の境界として、境界外の値を異常値として除去します。

以上の異常値除去の操作を行なったのち、軌跡データをGoogle Earth Pro上で可視化をすると以下のような結果が得られ、異常値がしっかりと除去できていることが確認できました。

さらに異常値除去を行った後、Median Filterを適用すると以下のような結果が得られました。これでこれまで目標としていた「全体的なジグザグの平滑化」および「立体道路下層の異常値除去」の2つが達成できていることがわかります。

滞留点検出

軌跡データの中から、しばらくその場所に滞在したことを示す点群を検出することを滞留点検出(Stay Point Detection)と呼びます。

閾値設定の工夫が必要になりますが、滞留点検出によってレストランやスーパーマーケットなどの任意の施設への訪問を検出することできますし、今回のような短時間の信号待ちを検出することも可能です。

滞留点検出のアルゴリズムはこちらの論文で提案されている、時間的な閾値と距離的な閾値を手動で決め、「何秒間以上、何メートル以上の移動をしなかった場合、それを滞留(信号待ち)とみなす」方法を実装します。

なお、同じ地点に再度訪れた場合にそれを滞留と見なさないように今回は「10秒以上かつ20秒以下の間において30m以上の移動がない点群を信号待ちとみなす」ように工夫を加えました。

from geopy.distance import geodesic

def Distance(lngDeg_i, latDeg_i, lngDeg_j, latDeg_j):
    return geodesic((latDeg_i, lngDeg_i), (latDeg_j, lngDeg_j))

def StayPointDetection(lngDeg=lngDeg, latDeg=latDeg, seq_len=df.shape[0], distThreh=0.03, maxTimeThreh=20, minTimeThreh=10):
    StayPointList = []
    start = 0
    end = seq_len - 1
    while end >= 0 :
        while end > start :
            dist = Distance(lngDeg[start], latDeg[start], lngDeg[end], latDeg[end])
            delta_t = 1 * (end - start)
            if dist <= distThreh and (delta_t <= maxTimeThreh and delta_t >= minTimeThreh):
                StayPointList.append(( start,  end))
                end = start - 1
                start = 0
            start += 1
        start = 0
        end -= 1
    return StayPointList

lngDeg = df["lngDeg_MeanFilter"].to_numpy()
latDeg = df["latDeg_MeanFilter"].to_numpy()
StayPointList = StayPointDetection()

さらに検出した滞留点群の重心を駐停車位置とみなし、滞留点群をその1点で置き換える操作をします。

for sp in StayPointList:
    sp_len = sp[1] - sp[0]
    df.loc[sp[0]:sp[1], "lngDeg_StayPoint"] = float(df.loc[sp[0]:sp[1], "lngDeg_MeanFilter"].sum() / sp_len)
    df.loc[sp[0]:sp[1], "latDeg_StayPoint"] = float(df.loc[sp[0]:sp[1], "latDeg_MeanFilter"].sum() / sp_len)

これによって、以下のような軌跡を得ることができました。信号待ち箇所およびスタート場所の駐車地点が比較的綺麗な点群で表現されたことがわかります。

まとめ

今回は、まずはじめに自動車内のスマートフォンのGNSS経路情報を可視化し、以下の問題を確認しました。

立体道路による計測誤差(異常値)
位置情報精度(数m~10m程度)の誤差
駐停車地点の位置情報のばらつき

さらにそれぞれの問題に対して

異常値検知
ノイズフィルタリング
滞留点検知

を行うこと綺麗な軌跡データを作成することができました。

なお今回はKaggleのGoogle Smartphone Decimeter Challengeコンペのデータセットの内、ほんの一部を抽出してうまくいったことを確認したに過ぎません。

例えば、今回適用した滞留点検知では「10秒間、30m以上の移動がなければ信号待ちとみなす」といった操作を行いましたが、この方法をデータの中身をよく確認せずに巨大なデータセットに適用した場合、「渋滞にハマっている自動車の軌跡を全て滞留とみなしてしまう」ことや、「実は歩行者のデータセットも含まれていて、その内のほとんどを滞留とみなしてしまう」なんてこともあると思います。

それぞれのデータセットに適した前処理を行うことが重要であることを伝えることも本記事の趣旨の一つですが、この記事を読んでいただいた皆様には、本記事を前処理の流れの一例として捉えていただき、日々の分析に役立てていただけますと幸いです。

最後までご覧いただき、ありがとうございました。

参考文献

滞留点検出アルゴリズム : Li, Quannan, et al. "Mining user similarity based on location history." Proceedings of the 16th ACM SIGSPATIAL international conference on Advances in geographic information systems. 2008.
軌跡データ全般 : Zheng, Yu. "Trajectory data mining: an overview." ACM Transactions on Intelligent Systems and Technology (TIST) 6.3 (2015): 1-41.

GNSS経路情報を可視化しながら前処理入門