车辆换道行为对交通流及交通安全有重要影响.车辆变更车道是最常见的驾驶行为之一, 也是最容易引发交通事故或交通拥堵的行为之一, 研究人员开展了广泛的车辆换道模型研究, 研究成果应用于车辆自适应巡航控制(adaptive cruise control, ACC)和计算机仿真(computer simulation, CS)^[1].一般的换道类型有3类, 强制换道(mandatory lane changing, MLC)、择机决策换道(discretionary lane changing, DLC)和随机换道(random lane changing, RLC).Li等^[2]认为现有对称两车道元胞自动机(symmetric two-lane cellular automata, STCA)模型没有考虑车辆换道时间的不均匀性, 提出了换道时间的离散性分段值, 并假设换道过程中, 换道车辆同时占据2条相关车道来近似模拟连续的换道过程, 结果表明考虑换道时间的元胞自动机模型的通行效率有所降低.Singh等^[3]也有类似的结论；Li等^[4]在此基础上提出换道的时空持续过程与当前车辆速度有关, 明确给出了换道时间函数, 并得出低速车辆的换道概率呈现二维分布的结论.但上述模型的理论假设并未有实测数据的支持.李玉洁等^[5]探讨了不同交通密度下施工作业区换道点的问题；李慧轩^[6]对车辆换道行为微观数据进行提取分析, 建立了换道持续时间与多因素的关系模型.

生存分析(survival analysis)是在综合考虑相关因素(内因和外因)的基础上, 对相关领域中与事件发生时间有关的问题提供统计规律的分析与推断方法^[7], 目前已广泛应用于相关领域.在交通运输研究领域, 生存分析方法己经逐渐应用于交通事故分析^[8-10]、行人或非机动车过街分析^[11-12]、实时停车预测^[13]、机动车寿命分析^[14]及机动车燃油消耗量建模^[15-16]、交通拥堵持续时间分析^[17-18]、交通事件持续时间分析^[19-22]、出行时间^[23]、居民出行行为分析^[24]及其他与持续时间相关的交通现象^[25]等方面, 目前生存分析方法被用来研究驾驶行为^[26-27].生存分析方法在研究交通事件持续时间及其影响因素方面具有一定的优势, 该方法考虑每个观测出现某一结局的时间长短, 可考虑删失数据(censored data), 且无需假定数据及误差服从某一分布.鉴此, 本文利用半参数生存分析方法, 建立相应的数学模型, 对高速公路施工作业区车辆车型和换道起点至合流点距离对换道持续时间规律影响进行探讨, 期望为交通管理与控制和未来智能网联车自主换道研究提供一定的理论基础.

1 模型和方法 1.1 基于风险的车辆换道持续时间

令X代表目标车辆换道持续时间, 为非负随机变量.它的累计分布函数为

$ F\left( t \right) = P\left( {X \le t} \right) = \int_0^t {f\left( x \right){\rm{d}}x}, \;\;\;\forall t \ge 0. $

(1)

式中:P(·)表示事件{·}发生的概率, f(x)为生存时间X的概率密度函数, t为任意给定时间.生存时间X的概率密度函数可表示为

$ f\left( t \right) = \frac{{{\rm{d}}F\left( t \right)}}{{{\rm{d}}t}} = \mathop {\lim }\limits_{\Delta t \to {0^ + }} \frac{{P\left( {t \le X \le t + \Delta t} \right)}}{{\Delta t}}, \;\;\;\;\forall t \ge 0. $

(2)

式中:Δt为时间t的瞬时增量.

目标车辆换道持续时间X的生存函数表示生存时间(即换道持续时间)大于t的概率, 具体的表达式为

$ \begin{array}{l} S\left( t \right) = P\left( {X > t} \right) = 1-F\left( t \right) = \int_t^\infty {f\left( x \right){\rm{d}}x}, \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\forall t \ge 0. \end{array} $

(3)

换道从车辆发生横向位移开始, 在换道时间X已经持续到t的情况下, 它在时间Δt内完成换道过程行驶至目标车道中心线的瞬时可能性用h(t)表示, 即为换道持续时间X的风险函数, 用公式表示为

$ \begin{array}{l} h\left( t \right) = \mathop {\lim }\limits_{\Delta t \to 0} \frac{{P\left( {t < X < t + \Delta t|X \ge t} \right)}}{{\Delta t}} = \\ \;\;\;\;\;\;\;\;\;\;\mathop {\lim }\limits_{\Delta t \to 0} \frac{{S\left( t \right)-S\left( {t + \Delta t} \right)}}{{\Delta tS\left( t \right)}}, \;\;\;\forall t \ge 0. \end{array} $

(4)

1.2 换道持续时间Cox比例风险模型

为分析某类影响因素对换道持续时间的影响, 采用Cox比例风险模型分析换道持续时间数据.

1.2.1 Cox比例风险模型建立

Cox比例风险模型是一种乘法风险率半参数模型, 不需要生存数据服从某种分布的信息, 估计的性质不依赖于选取的生存时间分布, 但却可以通过一个模型来分析生存时间的分布规律以及协变量对生存时间的影响.

换道持续时间X>0, 以及与X相关的影响因素变量(协变量)构成的向量为Z=[Z₁, …, Z_q], 则考虑协变量的风险函数

$ h\left( {t/\mathit{\boldsymbol{Z}}} \right) = {h_0}\left( t \right)\exp \left( {\mathit{\boldsymbol{\beta Z}}} \right), \;\;\;\forall t \ge 0. $

(5)

式中：h₀(t)为基准风险函数, 它是全部协变量都为零或标准状态下的风险函数, 一般未知；Z为影响因素向量, 可能的影响因素包括交通条件属性、驾驶人属性、环境属性、车辆属性；β为对应影响因素的系数向量.

分析某2个影响因素向量Z₁、Z₂的相对风险性时, 由式(5)推导出两协变量的相对风险比率HR的计算式为

$ {\rm{HR = h}}\left( {t/{\mathit{\boldsymbol{Z}}_1}} \right)/h\left( {t/{\mathit{\boldsymbol{Z}}_2}} \right) = \exp \left( {\mathit{\boldsymbol{\hat \beta }}\left( {{\mathit{\boldsymbol{Z}}_1}-{\mathit{\boldsymbol{Z}}_2}} \right)} \right). $

(6)

式中：$ {\mathit{\boldsymbol{\hat \beta }}}$为β的估计值.

该比值是一个常数, 即风险率是成比例的, HR为Z₁、Z₂的函数.若HR>1, 则相较于Z₂, 目标换道车辆在Z₁的影响下以较短时间来完成换道的概率更大, 亦即Z₂对目标车辆换道持续时间的影响大于Z₁.

1.2.2 换道持续时间生存函数及累积分布函数

由生存函数S(t)和风险函数h(t)之间的关系, 并结合式(5), 可得考虑影响因素下车辆换道持续时间生存函数, 即带协变量的累积生存概率

$ S\left( t \right) = {\left\{ {\exp \left[{-{H_0}\left( t \right)} \right]} \right\}^{\exp \left( {\mathit{\boldsymbol{\beta Z}}} \right)}} = {\left[{{S_0}\left( t \right)} \right]^{\exp \left( {\mathit{\boldsymbol{\beta Z}}} \right)}}. $

(7)

式中：H₀(t)为累计基准的车辆换道持续时间风险函数, 见式(8)；S₀(t)为基准的车辆换道持续时间生存函数, 表示协变量都被忽略(Z=0)时的换道持续时间函数, 见式(9).

$ {H_0}\left( t \right) = \int_0^t {{h_0}\left( u \right){\rm{d}}u} . $

(8)

$ {S_0}\left( t \right) = \exp \left[{-{H_0}\left( t \right)} \right] = \exp \left[{-\int_0^t {{h_0}\left( u \right){\rm{d}}u} } \right]. $

(9)

于换道车辆而言, 式(7)中的累计生存概率即为给定时间车辆仍没有完成换道的概率.由此可得出, 换道持续时间为t时, 考虑影响因素的车辆换道累积持续时间分布函数(Accumulative Duration of MLC Distribution Function)

$ F\left( t \right) = P\left( {X \le t} \right) = 1- S\left( t \right) = 1- {\left[{{S_0}\left( t \right)} \right]^{\exp \left( {\beta Z} \right)}}. $

(10)

1.3 换道持续时间Cox比例风险模型半参数估计

换道持续时间Cox风险模型中需要估计的有参数β和h₀(t)；由h₀(t)和S₀(t)的关系, 可将问题转化为估计S₀(t).估计的基本思路是先给出β的不依赖于h₀(t)的偏似然函数, 然后极大化偏似然函数, 给出β的估计, 再沿用导出乘积限估计的类似方法估计S₀(t).

1.3.1 协变量参数β的估计

容量为n的样本中, m个目标车辆的换道持续时间是非删失数据, 从小到大排列后依次为t₁, t₂, …, t_m, 即t₁ < t₂ < … < t_m, 对应的影响因素的协变量矢量依次为Z₁, Z₂, …, Z_m, 另外n-m个目标车辆的换道持续时间删失.在时间t_i处的风险集用R_i=R(t_i)表示, 即所有换道持续时间不短于t_i的车辆样本的集合.当多个目标车辆的换道持续时间值相等即数据中存在“结”的情况下, 目标车辆换道持续时间的Cox风险模型的偏似然函数为

$ L\left( \mathit{\boldsymbol{\beta }} \right) = \prod\nolimits_{i = 1}^m {\frac{{\exp \left( {\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{S}}_i}} \right)}}{{{{\left[{\sum\nolimits_{j \in R\left( {{t_i}} \right)} {\exp \left( {\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{Z}}_j}} \right)} } \right]}^{{d_i}}}}}} . $

(11)

式中：d_i为目标车辆中换道持续时间为t_i的个数；S_i为d_i个样本的协变量的和, 若m_i为t_i时刻换道完成的全部目标车辆的集合, 则${\mathit{\boldsymbol{S}}_i} = \sum\nolimits_{j \in {m_i}} {{\mathit{\boldsymbol{Z}}_j}} . $

该函数与h₀(t)无关, 并能极大化而给出β的估计值$\mathit{\boldsymbol{\hat \beta }} $；L(β)是偏似然函数, 不能由样本观测值的发生概率推导出来.文献[28]给出了以推断为目的时L(β)可被当作通常的似然函数处理的合理性.

对式(11)两边取对数, 得到对数似然函数LL(β)为

$ \begin{array}{l} LL\left( \mathit{\boldsymbol{\beta }} \right) = \\ \sum\nolimits_{{\rm{i = 1}}}^{\rm{m}} {\left( {\mathit{\boldsymbol{\beta }}{S_i}} \right)}- \sum\nolimits_{i = 1}^m {{d_i}\ln \left[{\sum\nolimits_{j \in R\left( {{t_i}} \right)} {\exp \left( {\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{Z}}_j}} \right)} } \right]} . \end{array} $

(12)

求式(12)的最大值, 就能求得偏似然函数的最大似然估计量.再对式(12)求关于β的偏导数得

$ \begin{array}{l} \frac{{\partial LL\left( \mathit{\boldsymbol{\beta }} \right)}}{{\partial {\mathit{\boldsymbol{\beta }}_h}}} = \\ \sum\nolimits_{i{\rm{ = 1}}}^m {\left[{{S_{ih}}-{d_i}\frac{{\sum\nolimits_{j \in {\bf{R}}\left( {{t_i}} \right)} {{Z_{jh}}\exp \left( {\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{Z}}_j}} \right)} }}{{\sum\nolimits_{j \in {\bf{R}}\left( {{t_i}} \right)} {\exp \left( {\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{Z}}_j}} \right)} }}} \right]}, h = 1, 2, \cdots, q. \end{array} $

(13)

式中：S_ih为S_i=[S_i1, S_i2, …, S_iq]中的第h个元素.

解q个非线性方程组∂LL(β)/∂β_h=0, h=1, 2, …, q可求出β的最大似然估计量$\mathit{\boldsymbol{\hat \beta }} $, 采用Newton-Raphson迭代方法求得数值解.

1.3.2 换道持续时间函数的估计

协变量Z下目标车辆换道持续时间函数为

$ S\left( t \right) = {\left[{{S_0}\left( t \right)} \right]^{\exp \left( {\mathit{\boldsymbol{\beta Z}}} \right)}}. $

(14)

只要估计出S₀(t), 即可得到协变量Z下目标车辆换道持续时间S(t)的估计.得到β的估计值后, 采用Kaplan-Meier乘积限估计方法得出S₀(t)的估计为

$ \begin{array}{l} {{\hat S}_0}\left( t \right) = \\ \left\{ \begin{array}{l} 1, t < {t_1};\\ \prod\nolimits_{t \ge {t_1}} {\exp \left[{-{d_i}\sum\nolimits_{j \in R\left( {{t_i}} \right)} {\exp \left( {\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{Z}}_j}} \right)} } \right]}, \;\;t \ge {t_1}. \end{array} \right. \end{array} $

(15)

由此可得到任影响因素Z下车辆换道持续时间生存函数的估计为(仅考虑t≥t₁)

$ \hat S\left( t \right) = {\left[{{{\hat S}_0}\left( t \right)} \right]^{\exp \left( {\mathit{\boldsymbol{\beta Z}}} \right)}} = {\left[{{{\hat S}_0}\left( t \right)} \right]^{\exp \left( {\sum\nolimits_{k = 1}^q {{{\mathit{\boldsymbol{\hat \beta }}}_k}{\mathit{\boldsymbol{Z}}_k}} } \right)}}. $

(16)

2 数据与实证分析 2.1 基础数据和换道持续时间样本提取

交通大数据分析中无人机(unmanned aerial vehicle, UAV)获取的交通数据较传统监控系统更全面, 不仅可以获取平均速度、密度、流率等宏观、中观交通数据, 而且可获取道路上包括车辆轨迹、换道和跟车数据在内车辆水平上的微观交通数据^[29].于2016年10月13日、14日12:00:00-18:00:00, 用大疆Phantom 4无人机对包茂高速陕西境内某路段施工作业区进行调研.对施工作业区进行高空监控摄像, 采集了自然驾驶状态下施工作业区的交通量数据及车辆的运动轨迹数据.分析车型Z₁和至合流点距离Z₂这2个因素对换道持续时间的影响.

样本为所调查的高速公路施工作业区路段的换道车辆, 总共收集了269个换道车辆的样本.如表 1所示中为样本示例, 每一行为一个换道车辆样本数据, 表示一次换道过程中换道车辆的各项指标, 包括样本车辆的编号(ID)、车辆类型(V_T)、目标车辆车速(v)及其开始换道的位置至施工作业区合流点的距离(l), 此外还包括了数据是否删失的信息.V_T为名义变量, 分为小型车(0)、中型车(1)和大型车(2)；l为数值型变量, 描述车辆换道起点距离合流点的距离；是否删失为分类变量, 数据删失编码为0, 数据非删失编码为1.

2.2 强制换道持续时间的总体分布及Cox模型估计结果

调研中发现, 在施工作业区, 大型车基本都靠右行驶, 极少量大型车有换道行为.图 1给出了车型和至合流点距离交互影响下车辆换道持续时间分布.所有样本的平均换道持续时间是9.87 s, 标准差为5.67 s.

分类变量“车型”以第1类作为参照水平, 得到的参数估计结果如表 2所示, 第6列给出了各协变量参数估计的显著水平；可以得出考虑车型和至合流点距离的施工作业区车辆强制换道耗时的Cox比例风险模型为

$ h\left( {t/\mathit{\boldsymbol{Z}}} \right)/{h_0}\left( t \right) = \exp \left( {-0.226{\mathit{\boldsymbol{Z}}_1}-0.004{\mathit{\boldsymbol{Z}}_2}} \right). $

(17)

车型(V_T)系数估计值为-0.226, 相对风险为0.798, 表明中型车的强制换道结束的风险略高于小型车, 是小型车换道耗时的0.798倍.但该系数的标准误差SE相对较大, 为0.157, 表明该系数的点估计量不精确而应该计算区间的估计量.该系数的95%置信区间为-0.226±1.96×0.157, 即(0.082, 0.534)；对应的相对风险的95%置信区间为(0.586, 1.085).风险率的置信区间包含1, 说明V_T在Cox模型中的重要性应该小心地解释, 亦即V_T对车辆在养护施工作业区换道持续时间的影响并不显著(检验p值为0.150＞0.05).换道起点至合流点距离(l)系数估计值为-0.004, 表明在施工作业区随着车辆至合流点距离的减小, 强制换道结束的风险逐渐增加, 系数的检验p值近似等于0, 说明l对车辆在施工作业区换道持续时间有显著影响.

图 2给出了根据生存分析的半参数方法考虑车型和至合流点距离两协变量计算得到的生存函数曲线.该图能反映施工作业区车辆换道实施过程中的换道在某时刻完成的概率.

由图形可知随着换道持续时间的增加生存函数的一般变化趋势；强制换道持续过程的生存时间大于5 s的概率约为80%, 生存时间大于10 s的概率为28%, 一半以上车辆的换道持续过程的生存时间处在5 ~10 s之间, 亦即20%的施工作业区车辆强制换道持续时间在5 s以内, 72%的车辆强制换道持续时间小于10 s, 一半以上(52%)车辆的强制换道持续时间处在5 ~10 s之间.

如图 3所示为施工作业区车辆换道持续时间的危险函数，由图 3可知施工作业区车辆换道持续时间的危险率是个增函数, 即某车辆换道持续到的时刻t越大, 其在时刻完成换道的可能性增加.

2.3 换道起点对强制换道持续时间的影响

为进一步深入分析换道起点至合流点距离(Dist)对换道持续时间的影响, 将换道起点至合流点距离分别划分为4组和7组, 分别得到不同的换道起点至合流点距离的生存函数和危险函数, 如图 4和5所示.

从图 4(a)可知, 在相同的生存概率下, 至合流点距离500 m内, 换道起点至合流点距离越近, 车辆强制换道持续时间越短；500 m以外, 车辆换道持续时间又开始降低；从图 5(a)的风险函数可知, 车辆换道时间持续到t, 在至合流点0~100 m范围内开始换道的车辆在t+Δt时刻完成换道的可能性最大.

为了得到分析结果, 对车辆换道起点至合流点距离进行分组, 得到结果如图 4(b)所示；至合流点距离400 m内, 随着换道起点至合流点距离越来越近, 车辆强制换道持续时间先减少后增加, 其中, 在至合流点距离为50~100 m范围内, 车辆的强制换道持续时间达到最短, 在至合流点距离为300~400 m范围内车辆的换道持续时间达到最大；400 m以外, 车辆换道持续时间降低, 后又有缓慢增大；从图 5(b)中的风险函数可知, 车辆换道时间持续到t, 在至合流点50~100 m范围内开始换道的车辆在t+Δt时刻完成换道的可能性最大.

3 结论

采用生存分析非参数方法, 构建了施工作业区车辆换道持续时间的风险分析模型, 以包茂高速陕西省境内某路段典型施工作业区使用UAV采集的车辆换道数据为例, 利用样本数据对模型参数进行估计, 定量分析了车型和换道车辆的换道起点至合流点距离等因素对车辆换道持续时间的影响.得出的主要结论如下：

(1) 强制换道持续过程的生存时间大于5 s的概率约为80%, 大于10 s的概率为28%, 一半以上处在5 ~10 s之间.

(2) 在至合流点距离为50~100 m的范围内, 车辆的强制换道持续时间达到最短；距离为300~400 m范围内达到最大.

(3) 至合流点距离400 m内, 随着换道起点至合流点距离越来越近, 车辆强制换道持续时间先减少后增加；400 m以外, 车辆换道持续时间降低, 后又有缓慢增大.车辆换道时间持续到t时刻, 在至合流点50~100 m范围内开始换道的车辆, 在t+Δt时刻完成换道的可能性最大.

研究成果研究可为交通管理与控制和未来智能网联车自主换道研究奠定一定的理论基础.在复杂多变的行车环境中, 多个因素动态协同影响车辆换道, 本文重点对影响车辆换道行为的车型和换道车辆换道起点至合流点距离2个因素进行了分析, 未来可针对目标车辆速度、不同交通时段、与目标车道前/后车的车头间距或时距等对车辆强制换道持续时间的影响进行研究.在实验路段选取方面, 在未来的研究中应选择多个路段多个施工作业区进行进一步深入对比研究.