程小云 屈霞萍 张学宇 邓亚娟

（长安大学运输工程学院，陕西 西安 710064）

常发性拥堵多由容量不足、出行需求过大和信号控制不佳这3方面的问题所造成，具有相同形成周期与相似时空范围等典型特征［1］。随着机动车保有量增长，一些常发性拥堵点会随拥堵传播演变为区域性拥堵，最终导致难以疏导的系统性拥堵，严重影响居民出行安全和效率。因此，有关常发性拥堵点的识别、预测及改善措施等的研究相继涌现，为缓解交通拥堵措施与策略的提出提供了有力支撑。然而，已有研究更关注对“源头”的疏导，忽视了对扩散路径的控制，使得拥堵蔓延难以抑制，这也是导致缓堵效果减弱的重要原因之一。由于常发性拥堵具有较为固定的产生和传播模式，因此，通过深入挖掘常发性拥堵点的传播机理，系统分析传播时空路径，有助于从源头和路径上抑制拥堵蔓延，为精准有效缓堵措施的提出提供理论依据。

针对拥堵判别方法，贾若等［2］主要从传统交通流理论方法和机器学习进行总结，并指出了有待解决的研究方向。本研究将拥堵识别方法归纳为拥堵评价模型的构建［3-4］、基于轨迹分段的拥堵检测法［5-7］以及基于网格划分的拥堵区域识别方法［8-10］这三大类。在拥堵评价模型构建方面，王曦等［3］建立了基于高斯混合分布的交通拥堵评价模型，通过实例验证了模型对混合交通流拥堵评价的有效性；王忻［4］提出了一种改进的隐马尔可夫模型的交通拥堵态势识别机制，有效拟合了城市道路相邻交叉口交通拥堵状况。在基于轨迹分段的拥堵检测法方面，付子圣等［5］利用基于密度的空间聚类以及模糊C均值聚类对出租车GPS轨迹段进行二次聚类，精细化识别道路拥堵时空范围；郭雪婷等［6］基于出租车GPS数据，利用边界矩形法识别道路，依据出租车行程速度判别路段的拥堵状况；Chaurasia等［7］基于地图匹配，通过出租车轨迹分段识别拥堵路段。在基于网格划分的拥堵区域识别方面，一般根据具体问题将研究区域划分为不同尺度的网格单元，将车辆轨迹拆分到各网格中，利用速度或拥堵模式等指标判别拥堵区域，李勇［8］将时间帧大小取为30 min，以轨迹速度为指标识别拥堵区域；杨海强等［9-10］采用不同时间帧大小，以车辆行程速度与拥堵模式指标共同识别常发性拥堵区域。

在交通拥堵传播机理方面，研究方法主要分为3种：一是基于理论模型的仿真分析，主要有改进SIS模型［11］、SIR模型［12］、反应-扩散动力学模型［13］、灾害蔓延理论［14］和基于元胞自动机构建模型［15］。二是在考虑拥堵时空关联关系的基础上构建交通拥堵传播模型，韦伟［16］考虑交通流时空不平稳性，提出了一种融合偏最小二乘回归的时空自回归模型（PLS-STAR）刻画了交通拥堵传播过程；欧吉顺［17］利用图论和机器学习构建了交通拥堵时空关联关系推理模型，对交通拥堵传播模式进行挖掘。三是运用相关算法对交通拥堵的传播进行分析，王美红［18］提出了CTO-Apriori算法，用于挖掘拥堵状态在时空上的演化轨迹，从而发现拥堵的源头区域；Saeedmanesh等［19］通过动态聚类算法考虑拥塞区域的增长和收缩来模拟拥塞形成过程；文献［20-22］都通过拥堵路段或网格的时空关系得到拥堵传播树，其中文献［20］利用马尔科夫模型，而文献［22］利用动态贝叶斯网络分析确定拥堵传播概率。

国内外学者通过借鉴其他研究领域的模型算法，或对已有模型、算法进行改进优化，更加精准地识别了拥堵路段或区域，从不同角度揭示了拥堵传播机理，但由于数据源的特征及理论模型与现实情况的差异，使得该研究方向仍存在一些局限性。具体表现为，在常发性拥堵区域识别中，已有研究采用试验车或检测器获取数据，其有限的空间覆盖率与人为因素的干预，导致无法进行大范围区域研究或识别精度不高；在拥堵传播机理分析中，仿真模型很难反映路网拥堵时的真实状态，且大部分研究忽略了对拥堵传播多路径及其概率的研究。总之，目前针对交通拥堵问题预防和治理的研究多呈现出静态性和局部性。

综上，本研究首先通过道路空间栅格化与时间切片构建时空立方体数据模型，通过立方体单元属性赋值，实现对出租车GPS轨迹点的时空映射；其次，综合栅格内轨迹行程速度与轨迹数两个指标进行各时间帧常发性拥堵区域的判别；最后，利用STC算法（Spatial-Temporal Congestions Algorithm）对常发性拥堵区域构建频繁拥堵传播树，并引入动态贝叶斯网络，通过贝叶斯估计进行参数学习，从而得到常发性拥堵传播路径及相应传播概率。该方法可从时空维度精细化辨识拥堵特征，动态分析网络拥堵传播方向及其可能性，为寻找及控制拥堵传播过程中的关键路径奠定基础。

1 研究区域与数据预处理

1.1 研究区域和数据概况

本研究所用数据为2017年4月10日至4月14日5天内的出租车行驶轨迹数据，共包含11 600辆运营的出租车，平均每天早高峰（07：00—09：00）约产生254万条数据，行驶轨迹覆盖了西安市绕城高速围合区域内的所有路段。由于早高峰时段较晚高峰更固定，更易发生由长距离通勤而产生的交通拥堵现象。因此，选取早高峰时间段07：00—09：00作为研究时段。

1.2 数据预处理方法

时空立方体模型是由二维空间坐标和一维时间坐标组成的一个三维立方体。其中，二维坐标系表示现实世界的平面空间，一维时间轴表示状态变化的时间序列。任意给定一个时间点对应的截面即为现实物理空间在该时间点的瞬时状态。该数据模型能够直观表达路网交通运行状态随时间变化的过程，为清晰展示拥堵传播的时空路径提供数据框架。

1.2.1 数据筛选

出租车GPS数据包括车辆ID、经度、纬度、GPS时间、速度、载客状态等字段。首先，按研究范围进行数据筛选；其次，按载客状态选取载客运营的车辆数据。出租车载客状态字段值为1表示载客，0表示空车。因出租车巡游揽客时存在车辆慢行寻客或停车待客的现象，无法反映所行驶道路的真实交通状态。因此，为提高数据可靠性，剔除GPS数据中载客状态字段值为0且速度很小的数据。

1.2.2 网格划分

常发性拥堵识别精度易受到网格大小的影响［23］。目前在交通拥堵研究中，采用的网格尺度主要有200 m×200 m［9］、250 m×250 m［10］、890 m×890 m［24］，或将研究区域划分为40行×40列的网格（每个网格尺寸约为2.7 km×0.926 km）［8］。网格尺度确定方法主要有经验取值法［8，10］和控制性因素取值法，如采用数据空间步长与失效网格数［9］，利用拥堵点空间聚类程度［24］作为影响拥堵识别精度的控制性因素来确定网格尺度。网格尺度过小时，网格内无法采集到足够多的有效轨迹点，影响拥堵识别精度；网格尺度过大时，则包含多条路段，拥堵识别结果难以精确到某一条路段上。因此，为使大部分网格只包含一条路段，根据西安市路网特征［25］以及已有研究［10］，采用250 m×250 m的网格划分研究空间单元，记为G（g1，g2），其中g1为行编号，g2为列编号。

1.2.3 时段划分

交通拥堵状况具有时变性，需将研究时段划分为多个时间帧以提高拥堵识别精度。现有拥堵研究中，时间帧大小取值有30 min［8］，10 min［9］，5 min［10］，甚至24 h［24］。然而，若时间帧粒度太小，易造成时空立方体单元内的车辆轨迹点过少，导致求取的平均速度以及轨迹数量难以表征网格内的拥堵状况；时间帧粒度太大，则严重影响拥堵识别的精度。因此，本研究综合考虑车辆轨迹点数与识别精度，取时间帧大小为10 min，则第i天的第j个时间帧编号记为Ti，（ji∈［1，5］，j∈［1，12］）。

1.2.4 轨迹拆分

轨迹拆分即将连续的轨迹点映射到上述构建的时空立方体模型中，其本质是通过时空匹配为每个轨迹点赋予立方体单元相应的空间单元格编号G（g1，g2）及时间帧编号Ti，j。具体轨迹拆分过程如图1所示。

如图1所示，任意提取第i天的某车辆轨迹，记为tr1=｛p1，p2，…，p10｝，该轨迹跨越了多个时空立方体单元，经过步骤（1）与步骤（2），各轨迹点被拆分到时空立方体单元中，并被赋予网格编号G（g1，g2）与时间帧编号Ti，j，如轨迹点p1被赋予网格编号G（5，5）与时间帧编号Ti，1，轨迹点p3被赋予网格编号G（5，5）与时间帧编号Ti，2。轨迹拆分前后的字段如表1和表2所示，网格编号与时间帧编号为轨迹拆分后的新增字段。

图1 出租车轨迹拆分Fig.1 Splitting of taxi trajectory

表1 轨迹拆分前数据表主要字段Table 1 Main fields of data table before trajectory splitting

表2 轨迹拆分后数据表主要字段Table 2 Main fields of data table after trajectory splitting

2 常发性拥堵区域识别

2.1 拥堵识别指标及阈值的确定

文中在已有研究［10］的基础上，考虑行程速度波动性，采用加权平均行程速度与车辆轨迹数相结合的方法进行拥堵网格识别。轨迹行程速度与轨迹数存在以下4种关系［9］：行程速度降低时轨迹数增加（趋于拥堵）；行程速度降低时轨迹数减少（异常现象，约占7.3%）；行程速度提高时轨迹数增加（异常现象，约占11.3%）；行程速度提高时轨迹数减少（趋于畅通）。因此，同时采用平均行程速度与车辆轨迹数判别拥堵可有效避免数据异常造成的识别精度损失。当特定时间帧内某一网格轨迹数和轨迹平均行程速度与所有时间帧下两指标对应的平均水平值相差较大，且轨迹平均行程速度低于平均水平值而轨迹数大于平均水平值时，即可判定该网格在这一特定时间帧内处于拥堵状态［10］。

需要注意的是，某时间帧某网格内会存在多条车辆轨迹，若采用轨迹行程速度的算术平均值［10］，个别过低或过高的速度会使识别结果出现偏差。因此本研究考虑了行程速度的分布特性，依据速度分段统计轨迹数，以区间轨迹数在总轨迹数中的占比作为区间各轨迹行程速度的权重，表示其对行程速度平均值的影响程度，从而计算轨迹行程速度的加权平均值，减小识别结果偏差。再者，本研究采用某一时空立方体单元内的拥堵模式表征车辆轨迹数与加权平均行程速度的共同状态，可记为二维向量MTi，j=[NTi，，jVTi，j]，NTi，j指时间帧Ti，j内的车辆轨迹数，VTi，j指时间帧Ti，j内所有轨迹加权平均行程速度，采用拥堵模式的差异表征两个指标的变化情况。

2.1.1 拥堵识别指标

（1）加权平均行程车速

按10 km/h的速度间隔将时间帧Ti，j内某一网格的轨迹行程速度划分为r个区间（即［0，10 km/h），［10 km/h，20 km/h）等），第l个区间有xl条轨迹，则时间帧Ti，j内某一网格所有轨迹加权平均行程速度VTi，j的计算如式（1）所示：

式中，NTi，j是时间帧Ti，j内某一网格的总轨迹数，wl是第l个区间的权重值，V ly是第l个区间内第y条轨迹的行程速度。

其中，wl、V ly分别按下式计算：

式中，pm是车辆轨迹按GPS数据时间戳排列后的第m个轨迹点，tz，n、tz，1分别是第z辆车的第n个轨迹点与第1个轨迹点的时间戳，D( )pm+1-pm是两轨迹点的欧式距离。

（2）拥堵模式差异

利用欧式距离量化同一网格内在时间帧Ti，1及时间帧Ti，2的拥堵模式差异，公式如下：

2.1.2 判别方法

时间帧Ti，j内某一网格被识别为拥堵时其拥堵模式MTi，j需满足条件：

式中，Ti，1，…，Ti，j，…，Ti，12是某网格内的12个时间帧；δ是判定阈值；是某一网格12个时间帧拥堵模式的平均值，

同时，加权平均行程速度与车辆轨迹数需足：

当某时空立方体单元内的交通运行状态同时满足式（6）、（8）与（9）时，即可认为时间帧Ti，j内该网格拥堵。对于判定阈值δ，值越小，算法准确率越高，漏检率相应越低［10］，因此文中取δ为0.4［10］。

2.2 常发性拥堵区域识别方法

已有研究［1，10，26］将某一网格或路段在一天内或多天内拥堵总时长是否超过某一阈值作为常发性拥堵区域的评判方法。本研究针对常发性拥堵产生的时空稳定性，在已有研究基础上精细化识别某一网格常发性拥堵发生的具体时段，为更细粒度的拥堵传播机理分析提供数据基础，具体过程如下。

若网格g在第i天的时间帧Ti，j内发生拥堵，则记为Di，g(Ti，j)=1，否则记为Di，g(Ti，j)=0，则网格g在时间帧Ti，j内发生拥堵的总天数表示为

若满足

则认为网格g在时间帧Ti，j常发性拥堵，其中ε为拥堵天数阈值。常发拥堵路段分为周、月和年常发拥堵路段［27］，本文主要研究周常发拥堵路段，依据文献［27］，ε取值为4。

3 交通拥堵传播机理分析

常发性拥堵具有发生、传播、消散过程的时空规律性。在时空立方体模型的框架下，基于上述交通拥堵识别算法，进一步分析拥堵传播机理，以确定拥堵传播网络中的关键路段，为提出针对性的阻断及预防拥堵传播的措施提供参考。首先采用动态识别传播关系的STC算法依据路网实际交通状况构建拥堵传播树，其次利用频繁项集挖掘算法重构频繁拥堵传播子树，最后引入动态贝叶斯网络建立拥堵传播模型。

3.1 构建拥堵传播树

引用动态识别传播关系的STC算法［22］构建拥堵传播树，其原理是从最后一个时间帧开始寻找该时间帧每一拥堵网格在前一时间帧内的父节点，即若网格G1在拥堵产生时间上位于网格G2之前，且空间上相邻，则G1为G2的父节点。该算法每次只挖掘两个相邻时间帧的拥堵依赖关系，无需跨越多个时间帧，因此算法的复杂度只取决于每个时间帧的拥堵网格数。该算法的伪代码如下所示。

给定常发性拥堵网格A（时间帧1）、B（时间帧2）、C（时间帧3），则上述算法中的t为3。首先筛选出时间帧3内的拥堵网格C，由于上一时间帧无拥堵网格，因此i=3的循环中得到的拥堵传播树为空；i=2的循环中，网格C为上一时间帧的拥堵网格，网格B被筛选出来，与网格C进行时间帧与空间上的相邻性判别，若同时相邻，则认为网格B向网格C产生拥堵传播，为两网格构建拥堵传播关系；否则不构成传播关系。同理可判断网格A与网格B。若最终判定网格B与C、网格A与B之间具有传播关系，则产生网格“A→B→C”的拥堵传播树。

3.2 频繁拥堵传播子树构建

本研究采用频繁项集挖掘算法寻找具有频繁传播关系的网格对，以重构频繁拥堵传播子树。设项集I=［A1，A2，…，An］（n=1，2，…，Q）是研究时段内全路网中所有具有传播关系网格对的集合，其中An指第n组具有传播关系的网格对；设事务数据库D=［B1，B2，…，Bm］（m=1，2，…，M）是所有拥堵传播树的集合，其中事务Bm指第m个由若干个项组成的拥堵传播树。

寻找具有频繁传播关系的网格对即挖掘传播树中的频繁1项集［An］。由于交通拥堵的传播具有动态性，在考虑拥堵传播关系的复现性时，存在事务Bm（拥堵传播树）同时包含多个项An（拥堵传播网格对）的现象。因此本研究提出一种以出现频率加权的频繁项集支持度计算方法，以考虑拥堵传播网格对在传播树中多次出现的情况，如式（12）所示，其中，事务权重，即出现频率，为该项集在传播树（事务）中出现的次数占总传播次数的比值。由于拥堵传播关系网格对An具有时空相邻特性，即可将相邻的两个时间帧视为一个时间帧组合，而每一组具有传播关系的网格对对应一个时间帧组合，则本研究中总传播次数为11个时间帧组合（共12个时间帧）。

频繁1项集［An］支持度按下式计算：

式中：P([An])为项集[An]发生的概率；C为拥堵传播树的总数；C([An])为包含1项集[An]的所有拥堵传播树以出现频率加权的支持度计数，

其中：αm为传播树Bm中1项集[An]的权重，即其所占时间帧组合数与总组合数11的比值；N(Bm)为决策变量，若[An]⊂Bm，则N(Bm)=1，否则N(Bm)=0。

如事务数据库中有3个拥堵传播树：B1=｛A→B→C→B→A→B｝，B2=｛D→A→B→E｝，B3=｛C→F→E→B｝，若按传统方法，1项集［A→B］的支持度计数为2，支持度为2/3；按本研究提出的加权频繁项集支持度计算方法，对于该1项集支持度计数的各传播树（事务）的权值分别为α1=2/11、α2=1/11。项集I=［A→B］的以时间帧出现频率加权的频繁项集支持度计数为C([An])=3/11，传播树（事务）总数C为3，则该网格对的支持度为1/11，表示拥堵由网格A向网格B传播的概率为1/11。若项集的支持度高于某一阈值，则称该项集为频繁项集。

本研究以3.1节中STC算法构建的拥堵传播树作为事务，采用加权的频繁项集支持度计算方法，取最大支持度的一半作为支持度阈值ε，挖掘具有频繁传播关系的拥堵网格的集合，通过重构得到频繁拥堵传播子树。

3.3 基于动态贝叶斯网络的拥堵传播模型

3.3.1 动态贝叶斯网络模型

动态贝叶斯网络是以概率网络为基础，将静态网络结构与时间信息结合，形成具有处理时序数据的新随机模型［23］。网络节点间具有相互独立关系，可大大减少建网过程中的计算量，并根据该独立关系进行贝叶斯网络的学习。因此，本研究利用动态贝叶斯网络建立拥堵传播模型，并进行参数学习，得到拥堵传播概率。

一个动态贝叶斯网络可以定义为（B0，B→），其中B0表示贝叶斯网络的先验分布，即动态贝叶斯的初始网络，B→表示由两个以上时间片段的贝叶斯网络组成的图形。具体的动态贝叶斯网络示意图如图2所示。

图2 动态贝叶斯网络［22-23］Fig.2 Dynamic Bayesian network

图2（a）是动态贝叶斯网络的先验网络；图2（b）为转移网络，各时间帧具有与初始网络相同的网络结构，相邻时间帧间会产生拥堵传播。假设有A、B、C三个网格，位置关系如图2所示，网格B的拥堵可能由同时间帧网格A及前一时间帧网格A与B引起，因此t-1时间帧内拥堵网格B的父节点包括三个：t-1时间帧内网格A与t-2时间帧内网格A和B。

3.3.2 参数学习及传播概率计算

极大似然估计适用于数据量大的情况，估计的参数能够较好地反映实际情况；而当采用的样本数据有限时，贝叶斯估计误差更小。由于本研究将轨迹数据拆分到了不同时空立方体单元内，并进行了常发性拥堵网格的识别，致使数据量减少，因此文中选择贝叶斯估计对构建的动态贝叶斯网络以及相关的样本数据进行参数学习来估计拥堵传播概率。

拥堵传播树中，时间帧t-1到时间帧t之间拥堵传播的概率为两时间帧之间的各网格间拥堵传播概率的乘积［22］，数学表达式为

其中，表示时间帧t内拥堵的网格i；Pa()是的父节点；是时间帧t内网格i的父节点拥堵的条件下其拥堵的概率。

一个包含多个时间帧的频繁拥堵传播树的传播概率计算需考虑所有相邻时间帧之间的传播概率，计算公式为

4 实例分析

本研究选取西安市南二环路东段区域为研究案例，采用上述方法识别常发性拥堵，分析其拥堵传播过程。如图3所示，所选区域主要由南二环路东段、翠华路、文艺南路、草场坡等构成。其中，南二环路为西安市快速路，承担着大量东西向的早晚高峰通勤出行，且与其相交的翠华路沿线分布了3所中小学，汇聚了大量的通学出行。因此，该区域道路早高峰通勤通学双重叠加、交通管理困难导致该区在早晚高峰呈现常态化拥堵。

图3 拥堵传播研究区域Fig.3 Study area of congestion propagation

基于上述方法构建频繁拥堵传播子树，并在时空立方体中对其传播过程进行可视化，如图4所示。该研究范围内常发性拥堵网格为（72，26）、（73，26）、（74，26）、（75，26）、（75，25）。网格（72，26）、（73，26）、（74，26）、（75，26）之间均由南二环路东段连通，因此网格之间的拥堵传播即为南二环路东段的拥堵传播，具体的传播过程如图5所示，早高峰阶段，该传播过程从07：20开始持续了100min。

图4 频繁拥堵传播子树时空立方图Fig.4 Space-time trajectory cube set of frequent congestion propagation subtrees

通过对动态贝叶斯网络进行参数学习，得到频繁拥堵传播子树的具体传播概率，如图5所示。07：20—07：30时间段内网格（72，26）与（73，26）开始拥堵，07：40—07：50时段内，网格（72，26）继续拥堵的概率为80%，通过南二环路东段，其拥堵传播至网格（73，26）内的概率为93%；网格（73，26）继续拥堵的概率为85%，传播至网格（74，26）的概率为78%，同理可得其他时段传播过程。由图5可知，所选区域拥堵常发，传播较频繁且传播概率均较大。

进一步计算各时间帧之间网格拥堵传播的概率，具体计算过程：由式（14）计算可得，时间帧07：20—07：30到时间帧07：30—07：40之间拥堵的传播概率为49%（≈80%×93%×85%×78%），其中既包括由其他网格传播的拥堵（93%、78%），也包括本网格内部传播的拥堵（80%、85%），同理可得其他时间帧之间网格拥堵传播概率。

为验证本研究所提方法对拥堵传播过程分析的合理性，利用百度地图中该区域早高峰路况时变过程进行比对，如图6所示。

由图6，拥堵从长安路立交处开始沿南二环路东段传播，与拥堵传播图5中的传播过程一致，研究结果与实际情况基本吻合，验证了方法的有效性。

图5 频繁拥堵传播子树图Fig.5 Frequent congestion propagation subtree

图6 百度地图路况时变过程图Fig.6 Road conditions of Baidu Maps

拥堵治理措施可分为短效与长效两方面，短效缓堵措施如识别高峰时段内巡逻车疏导的重点路段；合理组织交通诱导屏或广播等其他媒体进行车流诱导的最佳路线［28］；长效缓堵措施如交叉口信号优化；规划关键路段的平行分流路段；采取单向交通组织，对关键路段进行结构优化［29］等。由上述分析结果可知，网格（72，26），（73，26），（74，26），（75，26）频繁产生拥堵。经实地调查，长安路立交设有无信号灯的行人过街通道，易对机动车行驶产生极大干扰，因此可在高峰时段增加巡逻车进行疏导或重新设计行人过街通道的位置及类型；其次，利用交通诱导屏、广播等在车辆进入立体交叉口前进行路径诱导，以减少南二环路东段交通量；最后，根据各网格在各时间帧之间的传播概率，确定拥堵传播关键时段的关键区域，增加巡逻车进行疏导。

5 结语

本研究旨在从时空维度辨识拥堵特征，分析拥堵传播路径及其概率，为缓堵策略的实施提供寻找关键路径的理论方法。基于时空立方体模型，采用车辆轨迹数与行程速度指标共同识别各时间帧内常发性拥堵网格，为拥堵传播机理的准确分析奠定基础；利用STC算法构建拥堵传播树，提出了考虑拥堵传播在时间上复现性特征的频繁项集挖掘方法，以重构频繁拥堵传播子树，明确常发性拥堵传播路径；基于动态贝叶斯网络量化分析了网格间拥堵传播可能性，为抑制拥堵蔓延的缓堵措施提供实施的关键路段和时段。

文中所提出的方法可为缓堵措施的精准实施提供一种基于实时交通数据的局部拥堵传播趋势评价方法，实现对“源头”疏导的同时对传播路径采取阻断措施抑制拥堵蔓延，为助力科学治堵提供理论支撑。后续研究中，将对拥堵传播关键路径精细化识别以及拥堵阻断策略进行更深入研究。