Smartphone Path Estimation and Localization

对论文 Smartphone Based Indoor Path Estimation and Localization Without Human Intervention (IEEE Transaction on Mobile Computing 2022) 的阅读整理。

Smartphone Path Estimation and Localization

论文情况

• 标题：Smartphone Based Indoor Path Estimation and Localization Without Human Intervention
• 作者：Junyoung Choi , Gyujin Lee, Sunghyun Choi , Fellow, Saewoong Bahk
• 会议：IEEE Transaction on Mobile Computing 2022
• 源码：未开源

1 Introduction

应用：分析人员密集情况，在办公楼里合理安排 Wi-Fi AP、在商场里布局商品摆放等。

为获取用户的准确位置，初始位置十分重要，因此认为干预不可避免。

本文提出了 PYLON 系统：

• 使用被 Wi-Fi 或 BLE 覆盖的建筑入口（如门）作为地标。
• PYLON 运行在手机上，使用 IMU 传感器、Wi-Fi 和 BLE 的 RSSI 来收集用户的移动数据（步数、步频、走路方向）。服务器使用这些数据经 floor plan 算法估计用户轨迹和路标修正（landmark correction）。

PYLON 的中心组成即为 landmark correction。PYLON 根据 RSSI 使用数据挖掘算法生成一个用户经过的 virtual room。PYLON 估计用户实际经过的空间，在执行 floor plan mapping 后，PYLON 利用真实世界的 floor plan 估计用户经过每道门（door）的时间。由 IMU 构建的粗糙轨迹会通过用户经过的 door 来进行修正。

与其他定位系统不同，PYLON 不需要 RSSI 指纹、Wi-Fi AP 的位置、BLE beacon、用户的初始位置。也不需要额外硬件和软件的修正。

主要贡献：

• 设计了 PYLON，利用了室内环境真实的 door，包括 floor plan mapping、经过 door 时的时间检测和估计，进行预测轨迹的修正。
• 不需要额外数据、硬件和软件的辅助。
• 在两栋办公楼里，使用 Android 手机进行了实验，获得了高精度的结果。

2 Background

2.1 Infrastructure-Based Localization

• 使用相机、声学设备和红外设备进行室内辅助定位。
• 需要放置基础设施和特定的硬件。
• 可以实现较高的精度。

2.2 Fingerprint-Based Localization

• 利用 Wi-Fi 或 BLE 的广播信号。
• 需要测量在每个 point of interest 的 RSSI，建立 fingerprint 数据库。
• 基于指纹的定位难以使用在动态环境（AP 节点改变），数据库也需要在不同点进行复制以维护系统。

2.3 Model-Based Localization

• 使用 radio propagation model 进行 RSSI 的估算。
• 室内障碍物会影响广播信号的传播，降低模型的计算精度。

2.4 Dead Reckoning

航位推算（Dead Reckoning）：

• 集合加速度计、陀螺仪、磁力计等惯性传感器。
• 传感器误差累积迅速，只能短时间运用。
• 只跟踪行走轨迹，而不是用户的精确位置。

3 System Overview

图-1 PYLON 框架

• data collection module：手机 Wi-Fi/BLE 的 RSSI，读取 IMU 传感器数据。
• data processing module：计算 RSSI 值的差异，根据 IMU 数据粗略计算用户的行走路径。

3.1 Notable RSSI Signature

调查 Wi-Fi/BLE 信号的空间和时间特征后，发现接收到的信号强度会在用户穿过 door 时显著增强（或减弱），如图-2。因此 RSSI 可以用于检测用户穿过 door 的准确时间，以及区分不同的房间。

图-2 用户经过门时明显的 RSSI 变化

3.2 Smartphone Operation

• data collection module：手机周期性触发 Wi-Fi 扫描来从多个 AP 接收数据并测量 RSSI。BLE beacon 周期性发送广播包，手机蓝牙通过接收到的广播包计算 RSSI。IMU 模块则从手机中收集加速度、角速度、地磁等信息。

• data processing module：原始 RSSI 数据被转化为 RSSI 栈差异（RSSI stacking differences）来表示当前 RSSI 与其他 RSSI 信息的累积差异。这些数据用于限制设备依赖性以及生成 virtual room。

• IMU 传感器：陀螺仪数据用于沿时间积分来获取用户的步行方向，并通过 low-pass 过滤器来过滤噪声。通过计算加速度计的数据峰值个数计算步数，每一个峰值对应一步。通过计算加速度峰值的时间间隙，可以得到步频来计算步长。结合这些行走数据，可以推导出用户的行走路径。

• 手机将 RSSI 栈差异、用户的移动和路径上传到服务器。

3.3 Server Operation

图-3 Landmark Correction 过程

• landmark correction：使用 door 为路标或者参考点，因为用户更换房间时必然要过门。图-2 为该过程的四个步骤：

  1. 在 RSSI 栈差异上使用数据挖掘方法，生成 virtual room；
  2. 使用 floor plan mapping 算法匹配 virtual room 和真实世界的 floor plan；
  3. 使用 floor plan mapping 计算用户在什么时候通过门，此处假设已知这些门在世界 floor plan 的位置；
  4. 如果判断用户穿过门，则使用门的位置信息来修正用户路径的误差。

• particle filter：经过修正后的路径仍然可能存在穿墙而过的问题，因此使用 particle filter 来确保用户不会“撞墙”。particle filter 为非参数形式的 Bayesian 估计。

4 Path Estimation

4.1 Step Detection

PYLON 通过加速度计读取来发现用户的步行方式。由于脚跟落地时，加速度计三个轴上加速度的大小会出现一次极大值，如图-4。

图-4 使用加速度计检测加速度峰值

从加速度计中读取加速度值并计算 magnitude：

$a = \sqrt{a_x^2 + a_y^2 + a_z^2}$

然后使用 low-pass filter 过滤高频噪声，通过在线的方式进行过滤再使用下面的加权平均式：

$a_{avg}^t = (1 - \alpha) \cdot a_i + \alpha \cdot a_{avg}^{t-1 }$

其中 $a_i$ 为第 $i$ 个加速度计的原始 magnitude，$a_{avg}^t$ 为经过 low-pass filter 后的平均值，$\alpha$ 的默认值为 $0.8$。

之后，使用滑动窗口寻找 $a_{avg}$ 的峰值。如果 $a_{avg}$ 为窗口（大小为 $t_w$）内极大值，且其 slope sign change 小于 0，就考虑 $a_{avg}$ 为一个峰值。由于人类的步频小于 3Hz，因此设置 $t_w = 3Hz$。

4.2 Step Length Estimation

人的步长受到身高、体重、时间、环境等影响。此处使用 [1] 中的步长计算算法：

• 使用遗传模型（generic model）并收集用户数据来训练模型。收集了来自 23 个有不同生理特征用户的 4000 多步数据进行训练。选择频率模型（frequency model）作为 generic step model 并发现步长与步频有明显的线性关系。步长 $L_s$ 与步频 $f_s$ 之间的模型表达为：

  $L_s = a \cdot f_s + b$

• $a$ 和 $b$ 需通过实验标定。

4.3 Walking Direction

通过对陀螺仪角速度的积分求解行走方向。转弯时，身体的转轴始终垂直指向地面（重力方向）。由于重力方向的轴有最大的加速度值，先通过加速度计的重力加速度确定手机的方向，然后利用单轴旋转求取旋转角度。

4.4 Location Update

基本上，用户每走一步，系统就要更新用户的位置。用户的每一步都视为沿着行走方向以估计都步长前进。由于没有用户初始位置的先验知识，因此只检测用户的行走轨迹。用户第 $k$ 步的位置为：

$\begin{aligned} x_k &= x_{k-1} + L_s \cdot \cos (\theta) \\ y_k &= y_{k-1} + L_s \cdot \sin (\theta) \end{aligned}$

$\theta$ 为行走方向。如图-5，由于未知用户的初始位置，无法确定用户的真实位置。因此使用 landmark correction 来进行修正。

图-5 位置初始位置时，用户的可能行走轨迹

5 Landmark Correction Part1: Virtual Room Generation

5.1 RSSI Stacking Difference

如果环境中存在 $n$ 个 Wi-Fi AP 和 BLE beacon，则表示 $n$ 个 RSSI 采样为：

$R = [r_1, r_2, \cdots, r_n]$

$r_i$ 为第 $i$ 个 Wi-Fi AP 或 BLE beacon 的 RSSI 数值。由于 RSSI 的数值不稳定，在不同设备或相同设备不同时间的测量值不同。因此，使用 RSSI 数值的差异关系来消除 RSSI 值对设备的依赖性。

使用 RSSI stacking difference 来表示一个 RSSI 与其他 RSSI 之间的累积差异。RSSI stacking difference 表示特定地点和时间的 RSSI 间隔关系，比原始数据更加稳定。第 $i$ 个 RSSI 的 stacking difference 表示如下：

$w_i = \sum_{j=1}^n I(r_i - r_j > 0) (r_i - r_j)$

其中 $I$ 为特征函数：

$I(z > 0) = \begin{cases} 1, &\mathrm{if}\ z > 0 \\ 0, &\mathrm{else} \end{cases}$

5.2 Virtual Room Generation

PYLON 基于 RSSI stacking difference 生成 virtual room 并确定用户经过的房间数。这里使用数据挖掘和 K-means 聚类技术。K-means 中的类别数 $k$ 为用户所经过的房间数。使用肘部法则（elbow method）来确定类别数量。

图-6 用户行走轨迹示例

如图-6a，一个用户按字母顺序走过走廊和房间。手机在这个过程中收集 RSSI 数据并计算 RSSI stacking difference。服务器使用肘部法则确定聚类数并执行聚类算法，图-6 中 $k = 4$，即用户经过的房间（准确来说是空间）数。RSSI stacking difference 和路径估计信息都会被分为 4 类。

图-6b 即为 K-means 产生的 4 个 virtual room。每个 virtual room 的状态通过聚类随机分配。图中产生了 4 次状态迁移。

最终得到 virtual room 之间的关系，即 3 与 1、2、4 相连。

5.3 Virtual Graph Generation

基于 virtual room 和其关系，可以构建 virtual room 的逻辑 floor plan，使用无向图 $G = (N, E)$ 来进行表示。$n \in N$ 为 virtual room，$(u, v)$ 表示 virtual room $u$ 和 $v$ 互相连通，后面用 $G_v$ 来表示。

如图-6c 即为一个 $G_v$。

5.4 Physical Graph Generation

现在需要将 $G_v$ 与真实的 floor plan 进行匹配，在此之前，也需要一个真实 floor plan 的逻辑无向图 $G_p$。

图-7 真实环境的 floor plan

将房间和走廊都记为 $G_p$ 的一个节点，并如图-7a 将走廊分为多段。将走廊按垂直或水平分为多段，如 $n_1$ 和 $n_3$ 的分段。走廊垂直和水平的重合部分 $n_9$ 视为一个独立节点。因此，在生成 virtual room 时，$n_9$ 可能根据聚类结果而被视为一个独立区域。

根据真实 floor plan，可以得到 $G_p$ 如图-7b。

6 Landmark Correction Part2: From Floor Plan Mapping to Path Correction

有了 $G_v$ 和 $G_p$，下一步就是把它们进行比对合并。在此之前，先定义 candidate graph $G_c$（$G_p$ 的子集）并包含与 $G_v$ 相同数量的节点。

6.1 Candidate Graph Generation

一般上，$G_v$ 对节点数小于 $G_p$，除非用户经过了 floor plan 中的所有房间。因此存在多个 $G_p$ 的候选子图 $G_c$，但是不是所有 $G_p$ 的子图都可以是 $G_c$。如果子图中的一个节点与其他节点不相互可达，那么这个子图不能作为 $G_c$。

图-8 真实 floor plan 的子图

图-8a 为两个 $G_p$ 的子图，所有的节点相互可达，因此这两个图可以为 $G_c$，而图-8b 中的两个子图则不可以为 $G_c$。对于图-7b，如果 $G_v$ 包含 4 个节点，将会得到 39 个$G_c$，floor plan mapping 的目的就是选择一个 $G_c$ 作为 $G_v$。

6.2 Backbone Node Mapping

现在需要根据 $G_c$ 进行 floor plan mapping，首先需要进行 backbone node mapping。

使用中介中心度（betweenneww centrality），节点的中介中心度是对经过该节点的最短路径条数的度量。一般来说，叶节点的中介中心度为 0。节点 $v$ 的中介中心度可以以如下的形式定义：

$g(v) = \sum_{s \ne v \ne t} {\sigma_{st}(v) \over \sigma_{st}}$

其中 $\sigma_{st}$ 为从 $s$ 到 $t$ 的最短路径条数，$\sigma_{st}(v)$ 为经过 $v$ 的最短路径条数。将有非 0 中介中心度的节点称为 backbone node，然后在 $G_c$ 和 $G_v$ 之间执行 floor plan mapping。

图-9 Candidate Graph 和 Virtual Graph 的匹配

图-9 即为两个 $G_c$ 进行 backbone node mapping。图中左图为 $G_c$，右图为 $G_v$，将被匹配的点标记为红色。由于 $n_1$ 和 $n_3$ 中介中心度为 3 为 backbone node 且在 $G_c$ 中最高，因此最先参与匹配。对 39 个 $G_c$ 都进行 backbone node mapping，然后对剩下的点执行 dead-end node mapping。

6.3 Dead-end Node Mapping

dead-end node 为经过 backbone node mapping 后未被匹配的点。一般来说，dead-end node 即为叶节点。dead-end node 的节点仅依靠逻辑图 $G_c$ 和 $G_v$ 难以实现，因此，使用每个节点的相对位置进行匹配。

在 $G_c$ 当中，可以依靠真实 floor plan 确定每个节点的确切位置，将每个节点的位置定义为其所属房间（走廊）的中心。而 PYLON 根据 RSSI 生成 $G_v$ 时，也能由轨迹估计得到节点的相对位置。由此可以对 dead-end node 进行匹配。

余弦相似度（cosine similarity）通过内积测量两个非零向量的相似度，是对旋转的判断。向量 $A$ 和 $B$ 越接近，其余弦相似度越接近 1：

$cosine\ similarity = \cos (\theta) = {A \cdot B \over ||A|| \cdot ||B||}$

图-10 Dead-end Node 的匹配

此处以图-10 为例说明 dead-end 的匹配，dead-end node 的匹配被标记为绿色。图-10a 中 $G_v$ 和 $G_c$，计算 $n_1$ 到 $n_2, n_8, n_9$ 的三个向量（应该是通过坐标来建立），因此可以得到 9 个余弦相似度。然后根据相似度排序，最先匹配相似度最高的节点，本例中 $n_1$ 到 $n_8$ 的向量有最高的相似度，因此最先匹配 $n_8$。执行 dead-end node mapping 直到所有节点被匹配。

6.4 Final Candidate Graph Selectioon

图-11 Candidate Graph 的平均匹配相似度

如果 $G_c$ 精确匹配 $G_v$，则前面所计算的所有匹配 dead-end node 的余弦相似度为 1。因此选择平均余弦相似度最大的 $G_c$ 为 $G_v$。图-11 说明第 25 个 $G_c$ 的匹配度最好。

匹配完成后，就可以得知对应到真实环境下用户轨迹在走廊和房间之间的转移，从而根据 landmark 执行轨迹修正。

6.5 Door Passing Time Detection

由于 door 的位置固定，因此使用 door 作为 landmarks。如果知道用户穿过门的时间，则可以对用户的轨迹和旋转进行修正。使用图-6b 中的聚类结果为例说明如何得到穿过门的时间。

图-12 使用滑动窗口前后的聚类结果

然而仅根据聚类很难得到穿过门的时间，如图-12a，在状态发生转换时会产生乒乓效应（ping-pong effect），这意味着聚类结果不能很好反应穿过门的时间。

此处使用滑动窗口来使乒乓效应最小化，窗口大小定义为：

$W = {2 \over f_s} \cdot {1 \over I_{beacon}} \cdot S_N$

其中 $f_s$ 为步频，$I_{beacon}$ 为 beacon 间隔，$S_N$ 为手机可以接收到信号的 Wi-Fi AP 或 BLE beacon 到数量。由于 PYLON 生成 virtual room 时，手机对 RSSI 的采样则是一接收到信号就进行，因此设置一个滑动窗口单位为 beacon 的数量。

对原始聚类结果使用滑动窗口来生成新的状态：

$state_{new} (i + W) = mode(state_{old} (i+1, i+2, ..., i+W))$

$mode(\cdot)$ 表示出现最多的值。最终得到了图-12b 中的过滤结果。

6.6 Path Correction

图-13 对路径进行同步对齐

如果用户穿过门，则系统就能确定用户确切的位置和旋转：

• 首先，先得到用户穿过门的时间。
• 如图-12b，状态从 $n_3$ 变为 $n_6$ 十分重要，因为用户在这个时刻穿过门。
• 在图-13a 中，将穿过门的点标为蓝色，作为 anchor。
• 假设用户垂直的穿过门，PYLON 就能将轨迹进行旋转得到图-13b 的结果。
• 每一次用户穿过门，都要进行一次位置修正。

6.7 Particle Filter

图-14 Legacy、LC 和 PYLON 的路径估计

用户“穿墙而过”的情况仍然存在，如图-14。Legacy 为使用 dead reckoning 得到的路径，LC 为使用 landmark correction 得到的路径。

• 用户穿过门后，LC 表现出明显的优越。
• Legacy 需要提供初始的位置和旋转，仅使用 IMU 数据进行计算。
• LC 存在“穿墙而过”的现象。

采用粒子滤波解决“穿墙”的问题。

图-15 使用粒子滤波进行位置调整

初始时，粒子（采样点）生成在估计点附近，所有的粒子有相同的权重。根据用户的移动，每个粒子逐步位移，并被约束在 floor plan 下。第 $k$ 个位置的第 $i$ 个粒子的更新为：

$x_k^i = x_{k-1}^i + (L_s + \delta_k^i) \cdot \cos (\theta + \gamma_k^i) \\ y_k^i = y_{k-1}^i + (L_s + \delta_k^i) \cdot \sin (\theta + \gamma_k^i)$

其中 $L_s$ 和 $\theta$ 为估计步长和行走方向，$\delta_k^i$ 和 $\gamma_k^i$ 为 0 均值高斯噪声。

PYLON 检测是否每个粒子的移动破坏了室内约束（移动到了墙内），这些粒子则认为是 dead particle，并设权重为 0。剩下的 live particle 吸收 dead particle 的权重。然后经过重采样，在 live particle 周围产生新粒子，权重为 dead particle 的权重。这使得所有粒子的权重和超过 1，因此把所有粒子的权重进行归一化。最终取粒子的加权平均得到新的点的位置。

实验中去粒子数为 1000。

7 Performance Evaluation

• 实验路径：

图-16 实验路径

• 匹配精度：

图-17 匹配精度和余弦相似性

• 定位误差：

图-18 定位误差

• AP 或 beacon 数量的影响（初始放置 11 个设备）：

图-19 移走部分 AP 或 beacon 后的定位精度

• 行走距离和速度的影响：

图-20 行走距离和速度对定位精的影响

• 不同实验环境的影响：

图-21 两栋不同的办公楼

图-22 不同办公楼的定位精度

附件

1 Low-Pass Filter（低通滤波）

$y_n = \alpha x_n + (1 - \alpha) y_{n-1}$

其中 $x_n$ 为本次采样值，$y_n$ 为本次滤波输出值，$y_{n-1}$ 为上次滤波输出值。

2 Particle Filter（粒子滤波）

• https://blog.csdn.net/qq_36758914/article/details/123910958

3 Elbow Method（肘部法则）

附图-1 肘部法则

3.1 介绍

• 肘部法则用于指导 K-means 算法中 $k$ 的取值。
• 由附图-1，随着 $k$ 的增大，cost不断减少，当 cost 下降趋于缓慢时，取该值为 $k$。

3.2 算法

1. 对于 $n$ 个点的数据集，迭代计算 $k$ from 1 to n，每次聚类完成后计算每个点到其所属的簇中心的距离的平方和；
2. 平方和会逐渐变小，直到 $k$ 时平方和为 0，因为每个点都是它所在的簇中心本身。
3. 在这个平方和变化过程中，会出现一个拐点也即“肘”点，下降率突然变缓时即认为是佳的 $k$ 值。

3.3 代码

# clustering dataset
# determine k using elbow method

from sklearn.cluster import KMeans
from scipy.spatial.distance import cdist
import numpy as np
import matplotlib.pyplot as plt

x1 = np.array([3, 1, 1, 2, 1, 6, 6, 6, 5, 6, 7, 8, 9, 8, 9, 9, 8])
x2 = np.array([5, 4, 5, 6, 5, 8, 6, 7, 6, 7, 1, 2, 1, 2, 3, 2, 3])

plt.plot()
plt.xlim([0, 10])
plt.ylim([0, 10])
plt.title('Dataset')
plt.scatter(x1, x2)
plt.show()

# create new plot and data
plt.plot()
X = np.array(list(zip(x1, x2))).reshape(len(x1), 2)
colors = ['b', 'g', 'r']
markers = ['o', 'v', 's']

# k means determine k
distortions = []
K = range(1, 10)
for k in K:
    kmeanModel = KMeans(n_clusters=k).fit(X)
    distortions.append(sum(np.min(cdist(X, kmeanModel.cluster_centers_, 'euclidean'), axis=1)) / X.shape[0])

# Plot the elbow
plt.plot(K, distortions, 'bx-')
plt.xlabel('k')
plt.ylabel('Distortion')
plt.title('The Elbow Method showing the optimal k')
plt.show()

3.4 参考

• https://www.biaodianfu.com/k-means-choose-k.html
• https://blog.csdn.net/feichong621/article/details/109035295

参考

• [1] F. Li, C. Zhao, G. Ding, J. Gong, C. Liu, and F. Zhao, “A reliable and accurate indoor localization method using phone inertial sensors,” in Proc. ACM Conf. Ubiquitous Comput., 2012, pp. 421–430.