使用由图 1所示，26 242个节点和43 134条边组成的高精度深圳路网数据。为了提高算法效率，将城市分割为83$ \times $43个1000 m$ \times $1000 m的地理子区，以地理子区为单位进行交通异常检测。

图  1  深圳路网和地理子区GPS轨迹

采用由13 584辆配备GPS接收器的出租车采集的GPS数据，平均数据采集频率为15 s/次，使用2014年10月所有工作日共916 851 503个GPS坐标记录。由于出租车空载寻客时一般车速较慢，路段速度不具有代表性，因此删除出租车在空载状态下的GPS数据。由于少量GPS采集设备故障，还需要对GPS数据进行清洗，删除出租车日出行次数(由载客状态转变为空载状态为一次出行)大于500次的GPS设备所产生数据。为保证单位时间内充足的GPS记录数量，只使用6:00-23:00时间段的数据。经过以上数据预处理，平均每天有13 446 868条出租车坐标点记录被保留下来。

在一个地理子区中，满足以下条件的GPS轨迹点组成轨迹$T = {p_1} \to {p_2} \to \cdots \to {p_n}$。

1) GPS轨迹点有可能匹配到与其临近的任意一条路段上，所以需要采用缓冲半径公式，以GPS轨迹点为中心确定搜索范围，选择范围内的路段作为候选匹配路段，有：

式中，参数$\alpha $是道路网定位误差；参数$w$是单向道路的宽度；参数$\beta $是车辆定位的最大误差；参数$m$是车辆宽度^[15]。经调研，深圳市最宽道路为深南大道，路面净宽80 m，除去16 m绿化带宽度，单向道路宽32 m，其他参数的取值参考文献[15]：$\alpha = 5$，$\beta = 10$，$m = 2$，得到$\rho = 30$。本文将搜索范围略微扩大，取值35 m，并删除35 m范围内没有路段的GPS轨迹点。

2) ${p_i}$对应的时间标签值递增。

3) ${p_i}$和${p_{i + 1}}$之间的时间间隔小于45 s。

如图 1所示，虚线框表示一个地理子区，图中${p_1} \to {p_2} \to \cdots \to {p_n}$为地理子区中一条满足条件的轨迹。

为了计算出租车的行驶路径和路段的通行速度，使用ST-matching地图匹配算法^[16]将GPS轨迹匹配到路段上。该算法以路径为匹配对象，结合了路网拓扑结构和地理信息，是一种全局的地图匹配算法，对低频的GPS数据也能保持较高的匹配精度，并能方便地获取路径和路径的旅行时间。算法步骤如下：

1) 获得候选路径。

对于一条GPS轨迹$T = {p_1} \to {p_2} \to \cdots \to {p_n}$，当${p_i}$与路段$e_i^j$的距离小于$d = 35{\text{ m}}$时，将$e_i^j$加入${p_i}$的候选匹配边集合中，${p_i}$在$e_i^j$上的投影$c_i^j$则为${p_i}$的第j个候选匹配点，将n个轨迹点相邻两点的候选匹配点两两组合，可得到该轨迹的若干条候选匹配路径。

2) 空间分析。

ST-matching算法中，路网的地理信息和拓扑结构都被用来评价候选匹配路径与GPS轨迹的匹配程度，地理信息评价值用观测可能性来衡量，拓扑结构评价值用传递可能性来衡量。

观测可能性用${p_i}$与$c_i^j$的距离来衡量${p_i}$匹配到$c_i^j$的概率。定义${p_i}$匹配到$c_i^j$的观测可能性$f(c_i^j)$为：

式中，$x_i^j$是${p_i}$与$c_i^j$的距离，${\text{dist}}(c_i^j, {p_i})$越小，${p_i}$点匹配到e_i^j的可能性越大。传递可能性$V(c_{i - 1}^t \to c_i^s)$基于候选匹配路径上相邻两点$c_{i - 1}^t$和$c_i^s$的最短路径长度计算得到，计算公式为：

式中，${d_{i - 1 \to i}} = {\text{dist}}({p_{i - 1}}, {p_i})$为${p_{i - 1}}$与${p_i}$之间的欧式距离；${w_{(i - 1, t) \to (i, s)}}$为$c_{i - 1}^t$和$c_i^s$的最短路径长度。

结合观测可能性和传递可能性，$c_{i - 1}^t \to c_i^s$的空间分析函数定义为：

3) 路径匹配。

轨迹T的候选匹配路径组成候选匹配图${G'_T}({V'_T}, {E'_T})$，${V'_T}$是GPS轨迹点的候选点集合，${E'_T}$是相邻两点之间最短路径的集合，一条候选匹配路径${P_c}$可表示为${P_c}:c_1^{{s_1}} \to c_2^{{s_2}} \to \cdots \to c_n^{{s_n}}$，那么匹配到路径${P_c}$的可能性为：

匹配可能性最高的路径为$T$的最终匹配路径为：

图 2表示一条待匹配轨迹$T = {p_1} \to {p_2} \to \cdots \to $ ${p_n}$，三角形标识的点都是待选的候选匹配点，${p_5}$和${p_6}$各有两个候选匹配点，实线和虚线分别为两条候选匹配路径，经计算，经过$c_{\text{5}}^1$和$c_6^1$的路径匹配度更高。该路径的旅行时间为起点与终点，即${p_1}$与${p_6}$的时间标签的差值，${p_{i - 1}}$到${p_i}$的行驶速度为：

图  2  地图匹配算法示例

式中，${w_{(i - 1, t) \to (i, s)}}$是${p_{i - 1}}$到${p_i}$的最终匹配点之间的最短路径；${t_{(i - 1, t) \to (i, s)}}$是它们之间的时间差。如果同时间段一条路段上有多个速度记录，可计算得到路段的平均通行速度。

由于所使用的路网精度较高，可能存在因起点或终点不在同一条路段的细微差别，而将两条轨迹判定为经过不同路径的情况。为了保证路径判定有一定的容差率，将地理子区进一步划分为100 m$ \times $ 100 m的格子，如果两条轨迹途经的细分网格相同，则认为它们经过同一条路径。

本文采用基于密度的DBSCAN聚类算法^[17]对通过同一条路径所有轨迹的旅行时间进行聚类。该算法需要确定两个参数：EPS邻域和最小样本数MinPts，参数的选择对聚类的结果影响较大，较小的EPS和较大的MinPts会导致算法对路径旅行时间异常的判断不够敏感，反之，会过于敏感。使用欧式距离计算两点之间的距离为：

30 min是目前交通状态短时预测常用时间段之一^[18]，且能保证一条路径有足够数据量来进行聚类分析，所以，使用30 min作为时间窗，定义每个时间窗的数据为一个聚类数据集，并将聚类数据集中的每个实例除以数据集中实例的最大值，实现所有实例的归一化。根据文献[17]对参数EPS值和MinPts值选取方法的介绍，先确定MinPts取值为4，然后找到4-距离的突增位置来确定参数EPS的值，方法如下：首先计算每两个实例之间的欧式距离，找到与每个实例最邻近的第4个实例，二者之间的距离称为4-距离。然后对一个聚类数据集中实例的4-距离进行概率统计，发现绝大多数数据集的4-距离概率呈指数分布。最后对概率分布进行指数拟合，选取斜率为-1处对应的4-距离作为EPS值。图 3a是对某个聚类数据集的4-距离进行归一化后的概率分布图，图 3b是将图 3a中局部区域放大后，做拟合曲线斜率为-1处的切线，对应的横坐标为EPS取值。

图  3  DBSCAN算法参数选择

经过聚类后，聚类数据集中的实例被分成m个团簇，定义拥有最多实例的团簇为中心团簇，中心团簇中最大值为初始异常阈值，初步认为数据集中、所有小于初始异常阈值的对象是初始正常值，初始正常值的平均值$\left\langle t \right\rangle $为路径在所测时间段的正常旅行时间。由于存在少量聚类数据集，使实例的4-距离不满足指数分布，且拟合的效果不稳定，为减少DBSCAN参数对阈值计算结果的影响，根据常用异常检测方法三倍标准差准则^[19]，计算初始正常值的标准差$\sigma $。定义$\left\langle t \right\rangle + 3 \times \sigma $为路径在该时间段内的异常阈值。如图 4所示，图 4a为某条路径10月份工作日10:20-12:30所有旅行时间，虚线为每30 min路径旅行时间的异常阈值。每10 min计算一次待测路径的平均旅行时间，图 4b为某条待测路径在10:20-12:30不同时间窗进入地理子区时，路径的平均旅行时间，与对应时间窗的异常阈值对比，虚线空心点所表示的旅行时间值超出了异常阈值，由图可知，异常持续了50 min，即11:00-11:50。

图  4  路径旅行时间异常判断

交通事故或交通拥堵等交通异常事件的显著特点是造成路段通行速度降低，因此，大量的研究基于路段的速度、流量等属性进行路段状态识别。然而，这种异常识别方法也存在一些不足。一些路段的速度表现出不规则的波动给交通异常判断带来困难，如图 5a和图 5b所示。因此，选择路径的旅行时间作为研究对象，能综合考虑多个邻接路段的通行情况，减小异常检测误差，如图 4b所示。

图  5  路段速度计算

使用上一节阐述的方法计算各路径的异常阈值，图 6为某地理子区中检测出的4条异常路径旅行时间随时间的变化，4条路径的旅行时间均在第6个时间窗开始突增并超出异常阈值，于第8个时间窗附近达到相对高峰之后逐步恢复正常。由图 6a可知，该条路径正常通行时的旅行时间约50 s，而旅行时间异常值最高达280 s。

图  6  路径旅行时间随时间变化

为衡量异常路径覆盖路段受影响程度，引入R值进行计算：计算路径在各个时间窗的受影响程度=平均旅行时间/$\left\langle t \right\rangle $，取受影响程度最大的时间窗进行各路段$R$值的计算：

式中，${v_r}$为一个月中该时间窗的平均通行速度；${v_a}$为受影响程度最大时间窗的平均通行速度。图 7a所示为路径1在第8个时间窗依次经过路段${v_r}$和${v_a}$的对比。

图  7  路段速度分析

图 7b中加粗路段为异常路径覆盖路段，由$R$值分布可看出各路段受影响的程度。

为验证该交通异常检测结果，从新浪微博中搜索信息，发现如表 1所示事件。

该事件初始定位在图 7b所示的地理图标位置，事件描述与所测交通异常基本相符，微博所反映的事件发生在11:28，为图 6中箭头所示位置。而算法在11:00就检测出交通异常，比微博发出提前28分钟，微博发出时间基本位于路径旅行时间最大的时间窗内，从微博内容可看出此时交警还在赶往途中。图 6b展示了交通异常的产生和消散过程，图 7b展示了交通异常的集中区域。

将地理子区内所有交通异常路径经过的节点和边组成的网络定义为“基于路径的交通异常子网”，如图 8a中的粗线所示。建立基于路段的交通异常子网用于比较。定义路段在检测时间窗的车辆行驶平均速度${\mu _c}$，路段在工作日相同时段的车辆行驶平均速度为${\mu _m}$，标准差为${\sigma _e}$，当${\mu _c} < {\mu _m} - 3 \times {\sigma _e}$时认为路段的车辆行驶速度发生异常。提取车辆行驶速度发生异常的边组成的网络，得到“基于路段的交通异常子网”，如图 8b、图 8d所示。

图  8  异常网络规模

在发生交通事故的时间窗：基于路径的交通子网连成一个大团簇，如图 8a所示；而基于路段的交通异常子网分布较为杂乱，有几个零散的团簇，如图 8b所示。在交通运行正常的时间窗：基于路径的交通子网根本不存在，如图 8c所示；而基于路段的交通异常子网还是有一些零散的小团簇存在，这是由于个别路段偶发速度波动引起的，如图 8d所示。可以看出，基于路径的交通异常子网在交通正常与异常状态下的区别更为明显，有助于提高交通异常判断的准确性，上述结果也体现了基于路径的交通异常检测方法的优势：有利于减小偶发情况对整体检测结果的影响。

近年来，随着数据分析挖掘技术的发展，路段旅行时间和拥堵的判别方法已经日益成熟，出行者在出行时可以借助各种软件，查看到各个路段的实时路况信息，比如拥堵情况、运行速度、限速情况等。但是，局部的、可以自行消散的路段异常出现频繁，使交管部门较难准确识别出真正存在严重问题的区域。

本文针对有一定规模的道路交通异常突发情况，提出了基于路径旅行时间的交通异常检测算法。该算法不同于常用的基于路段通行能力和状态进行异常检测的方法，可以减少由单一的路段信息和车辆偶然因素导致的异常误判，检测出具有规模的、自行消散困难的、异于历史通行状态的交通异常网络，为交通管理部门提供可信度更高的交通异常监测结果。算法在现实路网和真实数据环境下的应用证明其检测的交通异常时空信息及时有效，所检测出的带有时间和空间信息的交通异常，能为交通管理者疏导交通提供重要管控和决策信息。

本文研究得到霍英东青年教师基金基础研究课题(141075)及中南大学创新驱动计划(2016CSX014)的资助，在此表示感谢。