ARISE Kaggle部活動記 #4

ARISE analytics 渡邉です。ARISE Kaggle部の活動記へようこそ。

今回、Google Smartphone Decimeter Challenge 2022に参加し20位銀メダルを獲得できました。コンペのポイントと解法の紹介を行いたいと思います。

コンペ概要

目的

Andoridスマートフォンで受信したGNSSのデータから、その受信位置を予測することが目的。

データセット

スマートフォンを車内に設置し、ドライブする。数十分のドライブの間中、スマートフォンと人工衛星は通信を行っている。このとき、いくつの人工衛星と通信できるかはドライブの経路に依存するため、データ量は経路ごとに増減することとなる。

私の解法

3つの手法を試した。(どれも私が位置から考えたわけではなく、他のkaggle参加者が公開してくださった手法を基にしています。添付する実装例は、その公開されたコードをリンクとしています。)

 

1. 衛星通信データWLS
   • 衛星データから各時点でスマートフォンの位置を推定する。
     
2. カルマンフィルタ平滑化
   • 各時点の位置を、時系列全体の位置の情報を使って尤もらしくする。
     
3. 加加速度（躍度）最小化
   • ドライブであるから、急加速減速はしないはず。それをカルマン平滑化の結果とバランスを取りながら修正する。
     

解法の説明

衛星通信データWLS

【この手法の目的】
衛星データから受信デバイスの位置と速度を特定する

【問題のモデル化】
変数を次のように与える。

受信デバイスの位置\( x,y,z \)、受信デバイスの時刻\(c \) 、  \(i\)番目の衛星の位置\(x^i_s, y^i_s, z^i_s\)、 \(i\)番目の衛星の時刻 \(t^i_s\)、受信デバイスと\(i\)番目の衛星の距離 \(\rho^i\)

ここで、次の関係が成り立つ。

\(\rho^i=\sqrt{(x^i_s-x)^2+(y^i_s-y)^2+(z^i_s-z)^2}\)　・・・(1)

\(\rho^i\)は既知であるように見える。なぜなら、例えば12:00:00に発信された信号を12:01:00に受信したならば、60x[電波の秒速]で求められるからである。しかし、ここで問題が発生する。発信された時刻は、衛星が記録するため正確だが、受信した時刻は受信デバイスによる記録のため、不正確である。この時刻のずれを\(d\)とおく。ここまでの議論を次式で表す。

\(\rho^i=c(t+d-t_s)\)

\(=c(t-t_s)+cd\)

\(=P+cd\)・・・(2)

\(P\)を疑似距離とよぶ。それぞれの時刻から計算できるが、 \(cd\)分ズレているため、こう呼ばれているのだろう。

(1),(2)を合わせて考えると、

\(P+cd = \sqrt{(x^i_s-x)^2+(y^i_s-y)^2+(z^i_s-z)^2}\)　・・・(3)

となる。既知の変数と未知の変数に分けて考えると、

◆既知
・\(P\)　（発信した時刻と受信した時刻は分かるから）
・\(c\)　（物理定数）
・\(x^i_s, y^i_s, z^i_s\)　（受信した信号に記録されている）

◆未知
・\(x,y,z,d\)

よって、最低4機の衛星と通信できれば受信デバイスの位置が求められることが分かった！

しかし、4機以上の衛星と通信できることの方が多い。このとき、ノイズを多分に含むため、すべての方程式を満たす解は見つからない。このようなとき、できるだけそれぞれの方程式を満たさない量（正確な呼び方を筆者は知らないが、誤差と呼ぶことにする）を小さくする点を解とすることが考えられる。誤差を小さくするよう逐次的に計算する方法を逐次最小二乗法とよぶ。さらに工夫すると、衛星ごとにデータの信頼性は異なると考えられるため、衛星ごとに誤差に重みをつけることもできる。これはWLS（Weighted Least Squer）とよばれる。

また、各時点の衛星の速度も得られるため、同様に各時点の受信デバイスの速度も求めることができる。

【このコンペにおける工夫点】
実装はこちら

位置と速度それぞれをscipy.optimize.least_squaresによって行った。工夫点は次の二つ。

1.信号の信頼性がデータに含まれていたため、それで重みづけした。
・RawPseudorangeUncertaintyMeters
・PseudorangeRateUncertaintyMetersPerSecond

2.WLSの損失関数にはコーシーロスを採用した。(\(loss(z) = ln(1+z)\) )
・デフォルトの線形ロス(\(loss(z) = z\) )よりもスコアが向上した。
・外れ値の影響を抑えられたからと考えている。外れ値により異常に大きな誤差が出てしまった場合でも、比較的小さな誤差をみなす。

 

 

カルマンフィルタ平滑化

【この手法の目的】
ノイズが載った時系列位置データから、ノイズを取り除く

【問題のモデル化】
WLSにより各時刻における位置データが得られたが、これらのデータは真の位置にノイズが載ったものと考える。また、ある時刻と次の時刻の間、車は等速直線運動に従うとする。これを状態空間モデルを用いて表現する。

\(\begin{bmatrix} x^{(t+1)} \\ y^{(t+1)} \\ z^{(t+1)} \end{bmatrix}\)=\(\begin{bmatrix} 1 &0 &0 \\ 0 &1 &0 \\ 0 &0 &1 \end{bmatrix}\) \(\begin{bmatrix} x^{(t)} \\ y^{(t)} \\ z^{(t)} \end{bmatrix}\)+\(\begin{bmatrix} w_x^{(t)} \\ w_y^{(t)} \\ w_z^{(t)} \end{bmatrix}\)

\(\begin{bmatrix} o_x^{(t+1)} \\ o_y^{(t+1)} \\ o_z^{(t+1)} \end{bmatrix}\)= \(\begin{bmatrix} 1 &0 &0 \\ 0 &1 &0 \\ 0 &0 &1 \end{bmatrix}\) \(\begin{bmatrix} x^{(t+1)} \\ y^{(t+1)} \\ z^{(t+1)} \end{bmatrix}\)+\(\begin{bmatrix} \mu_x^{(t+1)} \\ \mu_y^{(t+1)} \\ \mu_z^{(t+1)} \end{bmatrix}\)

\(W^{(t)} = \begin{bmatrix} w_x^{(t)} \\ w_y^{(t)} \\ w_x^{(t)} \end{bmatrix}\),\(W^{(t)} \sim N(V^{(t)}, \sigma_v^{(t)})\),\(V^{(t)}=\begin{bmatrix} v_x^{(t)} \\ v_y^{(t)} \\ v_z^{(t)} \end{bmatrix}\)

\(M^{(t)}=\begin{bmatrix} \mu_x^{t} \\ \mu_y^{t} \\ \mu_z^{t} \end{bmatrix}\),\(M^{(t)}\sim N(0,\sigma_x^{(t)})\)

状態空間モデルが決定すれば、あとはカルマンフィルタ平滑化をかけるだけである。詳細はこの本が詳しい。

【このコンペにおける工夫点】
実装はこちら

1.各時点のシステムノイズの共分散を推定した
・共分散はハイパーパラメータとして扱いカルマン平滑化を行う。今回はWLSより共分散を推定することができた。
http://ceres-solver.org/nnls_covariance.html

2.Xiaomiのシステムノイズの共分散を大きくした
・このコンペは各時点の速度も正解データとして与えられていた。その結果、Xiaomiのスマホだけ速度の推定精度が悪いことがわかっていた。これをカルマンフィルタ平滑化に反映するため、Xiaomiスマホの推定時だけ \(\sigma_v\)を10倍した。

 

加加速度（躍度）最小化

【この手法の目的】
予測誤差を最小化しつつ、加加速度を最小化する

【問題のモデル化】
ある時刻と次の時刻の間、車は緯度\(\phi(t)\) ・経度\(\psi\)それぞれにおいて、等加加速度運動に従うとする。

公式を参考に、緯度について書くと、\(\phi(t+\delta t) = \phi(t) + \delta t \dot{\phi(t)} + (1/2)\delta t^2 \ddot{\phi}(t) + (1/6)\delta t^3 \dddot{\phi}(t) \)となり、位置、速度、加速度を状態としたシステム方程式を考えると、次式のように表現できる。

\(\begin{bmatrix} \phi(t+\delta t) \\ \dot\phi(t+\delta t) \\ \ddot\phi(t+\delta t) \\ \psi(t+\delta t) \\ \dot\psi(t+\delta t) \\ \ddot\psi(t+\delta t) \end{bmatrix}\)=\(\begin{bmatrix} 1 & \delta t &(1/2)\delta t^2 & 0 & 0 & 0\\ 0 & 1 & \delta t & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & \delta t & (1/2)\delta t^2 \\ 0 & 0 & 0 & 0 & 1 & \delta t \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}\)\(\begin{bmatrix} \phi(t) \\ \dot\phi(t) \\ \ddot\phi(t) \\ \psi(t) \\ \dot\psi(t) \\ \ddot\psi(t) \end{bmatrix}\)+\(\begin{bmatrix} (1/6)\delta t^3 & 0 \\ (1/2)\delta t^2 & 0 \\ \delta t & 0 \\ 0 & (1/6)\delta t^3 \\ 0 & (1/2)\delta t^2 \\ 0 & \delta t \end{bmatrix}\)\(\begin{bmatrix} \dddot\phi(t) \\ \dddot\psi(t) \end{bmatrix}\)

上記は連続時間についての式であるため、データに即して離散時間で表す。ただし、サンプリング周期は十分小さいとする。

\(\begin{bmatrix} \phi[n+1] \\ \dot\phi[n+1] \\ \ddot\phi[n+1] \\ \psi[n+1] \\ \dot\psi[n+1] \\ \ddot\psi[n+1] \end{bmatrix}\)= \(\begin{bmatrix} 1 & \delta t &(1/2)\delta t^2 & 0 & 0 & 0\\ 0 & 1 & \delta t & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & \delta t & (1/2)\delta t^2 \\ 0 & 0 & 0 & 0 & 1 & \delta t \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}\)\(\begin{bmatrix} \phi[n] \\ \dot\phi[n] \\ \ddot\phi[n] \\ \psi[n] \\ \dot\psi[n] \\ \ddot\psi[n] \end{bmatrix}\)+\(\begin{bmatrix} (1/6)\delta t^3 & 0 \\ (1/2)\delta t^2 & 0 \\ \delta t & 0 \\ 0 & (1/6)\delta t^3 \\ 0 & (1/2)\delta t^2 \\ 0 & \delta t \end{bmatrix}\)\(\begin{bmatrix} \dddot\phi[n] \\ \dddot\psi[n] \end{bmatrix}\)

簡単のため、次のように置きなおす。

\(\Phi[n]\)=\(\begin{bmatrix} \phi[n] \\ \dot\phi[n] \\ \ddot\phi[n] \\ \psi[n] \\ \dot\psi[n] \\ \ddot\psi[n] \end{bmatrix}\), \(A=\begin{bmatrix} 1 & \delta t &(1/2)\delta t^2 & 0 & 0 & 0\\ 0 & 1 & \delta t & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & \delta t & (1/2)\delta t^2 \\ 0 & 0 & 0 & 0 & 1 & \delta t \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}\) , \(B=\begin{bmatrix} (1/6)\delta t^3 & 0 \\ (1/2)\delta t^2 & 0 \\ \delta t & 0 \\ 0 & (1/6)\delta t^3 \\ 0 & (1/2)\delta t^2 \\ 0 & \delta t \end{bmatrix}\),\(U[n]=\begin{bmatrix} \dddot\phi[n] \\ \dddot\psi[n] \end{bmatrix}\)

方程式を整理すると次になる。

\(0=\begin{bmatrix} A & -I \end{bmatrix}\)\(\begin{bmatrix} \Phi[n] \\ \Phi[n+1] \end{bmatrix}\)+\(BU[n]\)

全ての時点についてまとめて表す。

\(0\)=\(\begin{bmatrix} A & -I & 0 & \cdots & & 0 \\ \vdots & & & & & \\ 0 & \cdots & A & -I & \cdots & 0 \\ \vdots & & & & & \\ 0 & \cdots & 0 & 0 & A & -I \end{bmatrix}\)\(\begin{bmatrix} \Phi[1] \\ \vdots \\ \Phi[n] \\ \vdots \\ \Phi[N] \\ \end{bmatrix}\)+\(\begin{bmatrix} B & 0 & \cdots & & 0 \\ 0 & B & & & \\ \vdots & & \ddots & & \\ & & & & \\ 0 & & & & B \end{bmatrix}\)\(\begin{bmatrix} U[1] \\ \vdots \\ U[n] \\ \vdots \\ U[N] \end{bmatrix}\)

これが制約条件となる。簡単のために文字を置きなおす。

\(0 = DX+EJ\)

また、予測誤差を定義する。

\(\begin{bmatrix} \hat\Phi[1] \\ \vdots \\ \hat\Phi[n] \\ \vdots \\ \hat\Phi[N] \\ \end{bmatrix}\)-\(\begin{bmatrix} 1 & \cdots & 1 & 0 & \cdots & 0 & 0 & \cdots & & \\ 0 & \cdots & 0 & 1 & \cdots &1 & 0 & \cdots & & \\ & & & & & & \ddots & & & \\ 0 & \cdots & 0 & & & & & 1 & \cdots & 1 \end{bmatrix}\)\(\begin{bmatrix} \Phi[1] \\ \vdots \\ \Phi[n] \\ \vdots \\ \Phi[N] \\ \end{bmatrix}\)

これも簡単のために置きなおしておく。

\(Y\)-\(C\)\(X\)

ここまでで、最適化は次式のように定義できる。

◆決定変数
\(X=\begin{bmatrix} \Phi[1] \\ \vdots \\ \Phi[n] \\ \vdots \\ \Phi[N] \\ \end{bmatrix}\),\(ｍJ=\begin{bmatrix} U[1] \\ \vdots \\ U[n] \\ \vdots \\ U[N] \end{bmatrix}\)

◆目的変数
\(U^TRU+(Y-CX)^TL(Y-CX)\)
\(R\)と\(L\)は、加加速度最小化か予測誤差最小化どちらを重視するか決定する適当な重み

◆制約
\(0 = DX+EJ\)

【このコンペにおける工夫点】
実装はこちら

コースごとの\(R\)と\(L\)の調整
・カルマン平滑化により予測できていないことが推察された場合、Lを小さくし、予測誤差の重みを小さくした。
・例：両側に木々が生い茂っている見通しが悪いコース

終わりに

GNSSという歴史ある技術ではあるため、キャッチアップの量が膨大でしたが、ドメイン知識をつける楽しさが味わえたコンペでした。また、時系列データ平滑化や最適化についても新しい手法に触れられました。

ARISE analyticsでは、kaggle等の分析コンペティションで上位成績を残すと、インセンティブとして報奨がもらえるARISE Tech Master制度があるので引き続き金メダルを目指し、努力していきたいです。

分析はチームワークが大事です。役割を決めお互いの強みを持ち合うことで、一人で行うより何倍も効率的に分析課題に取り組めます。我々と一緒にKaggle部を盛り上げてくださる方はこちらのページからご連絡お待ちしております!!!!