基于元胞自动机的随机车流模拟与桥梁疲劳响应分析

周永兵，李 睿，刘海证,2，何永伟

(1.昆明理工大学 建筑工程学院，云南 昆明 650504；2.云南省建筑工程设计院，云南 昆明 650041 3.云南睿德道路桥梁工程设计有限公司，云南 昆明 650500)

0 引言

随着我国公路交通基础设施建设的快速发展，公路的服务水平和通行能力得到了极大的改善，然而由于交通流量的不断增加，在车辆荷载的反复作用下公路桥梁面临的疲劳问题越来越严重[1]。云南省处于高原山区，受地形地质条件的限制，截止到2017年低等级公路所占比例达到92.45%，同时全省货运量逐年增加，仅2017年公路运输所占比例达到90.31%[2]。另一方面，由于云南拥有丰富的水电及矿产资源，用于大型机组设备及矿物原材料运输的重载车辆较多[3]，致使云南山区低等级公路上超载超限现象异常突出，从而加剧了桥梁结构的疲劳损伤。采用准确、合理的疲劳车辆荷载模型用于桥梁抗疲劳分析和疲劳可靠性评估是确保桥梁安全使用的关键[4]，然而我国现行的《公路钢结构桥梁设计规范》(JTG D64—2015)[5]提出的疲劳荷载计算模型没有考虑到不同地域交通荷载的差异性，其疲劳荷载未必适合在不同地区使用。所以对于云南这样的高原山区省份，对山区低等级公路桥梁交通荷载调查和车辆荷载作用效应研究意义重大。

目前，采用动态称重(WIM)系统监测桥梁结构日常运营的车辆荷载，记录车辆通过监测断面的车重、车速、轴重、轴距、车时距等瞬时信息，利用实测的交通荷载参数进行随机车流模拟及结构响应分析得到了广泛发展[6-7]。普遍采用的蒙特卡洛(Monte Carlo)方法是基于实测交通荷载信息构造概率分布模型，进一步通过随机抽样得出交通特征参数来模拟随机车流[8]，其高效简单却无法模拟出车流的微观运动状态。元胞自动机(Cellular Automata)模型通过一定的演化规则可以全面地描述车辆的微观变化，该方法在交通荷载模拟中有较好的应用与发展前景[9]。

本研究根据云南某低等级公路桥梁上WIM系统的监测数据，对该区段的典型车辆进行分类统计和交通荷载参数研究，得到了各类车型比例、实时交通流量分布统计以及车重、车速、车长、轴距和轴重、车时距的概率分布模型及特征参数。在建立统计模型的基础上，采用多个元胞来模拟一辆车的精细化多元胞自动机模型，实现了对实际交通荷载状况的模拟和车辆作用下桥梁结构的荷载效应分析。

1 实测交通荷载参数研究

1.1 车型分类及实时交通流量分布统计

针对云南山区的低等级公路桥梁实际运营交通荷载的调查，在某三级公路上一新建单箱三室简支钢箱梁桥上安装WIM系统。通过对该路段进行交通监测，在2018年3月份连续29 d采集到49 199个有效车辆荷载数据，月平均日交通量为1 696 veh/d。参照《中国汽车车型手册》[10]的车型类别，利用系统聚类法对交通荷载数据进行统计分析得到云南省山区低等级公路车型分类及比例，如表1所示。结果表明：云南山区低等级公路交通荷载车型包括5类，第1类车型占95.17%，第2类和第3类车型各占1.8%，第4类车型占0.46%，第5类车型占0.78%，充分表明由于低等级公路行驶条件和通行能力限制，第1类车占绝大部分，但也存在不少大型车辆。

一天的交通流量在随时间不断改变，根据总车流量及各类车型的实时交通流量统计表明：总车流量分别在9:00左右和13:00—18:00时段内出现车流高峰，第1类车与总的车流量变化规律一致；其余4类车流量较小，在24 h内无明显变化规律。

表1 车型分类及比例Tab.1 Vehicle type classification and proportion

1.2 车重

车重一般服从对数正态、正态、威布尔和高斯混合分布等。车重分布统计首先假设其服从上述某种分布，通过最大期望(EM)算法编程进行参数极大似然估计[3]，其原理是一次随机试验有多个可能产生的结果X1，X2，…，Xn，假设结果Xi的发生概率与参数θ相关，即Xi发生概率记为P(Xi，θ)，则最大似然法便是寻找使结果Xi发生概率最大化的参数θ。最后通过拟合优度检验方法K-S法进行分布拟合检验，确定其合理的分布类型和概率密度函数[11]。

通过上述步骤得到各类车型的车重分布结果如表2所示：第1类车服从对数正态分布，其余4类车服从三峰正态分布，在第2类车重中，某峰值达到52 t左右，第3类达到57 t左右，第4类达到47 t左右，第5类达到64 t左右。在《汽车、挂车及汽车列车外廓尺寸、轴荷及质量限值》(GB1589—2016)[12]中，三轴车最大允许总质量限值25 t，四轴车36 t，五轴车43 t，六轴车49 t，可见低等级公路由于监管不力，车辆超载现象严重，桥梁结构存在极大的安全隐患。

表2 车重分布结果(单位：t)Tab.2 Result of vehicle weight distribution(unit：t)

1.3 车速

通过对车速的分布拟合和K-S检验，5类车的车速均服从正态分布，各类车型的车速分布结果如表3所示，可见第1类车速均值为37 km/h，第2类为30 km/h，第3，4，5类车速在28 km/h左右，整体车速较低，符合低等级公路的行驶条件。

表3 车速分布结果(单位：km/h)Tab.3 Result of vehicle speed distribution(unit：km/h)

1.4 车长

车长是交通流中车辆位置确定的重要参数，通过对车长的分布拟合和K-S检验，5类车的车长均服从正态分布，各类车型的车长分布结果如表4所示。

表4 车长分布结果(单位：m)Tab.4 Result of vehicle length distribution(unit：m)

1.5 轴距和轴重

在精细化随机车流模拟中，车长和车重参数不能准确反映实际交通荷载状态[13]，车辆荷载作用下桥梁结构经受的应力循环次数和幅值取决于轴距和轴重[4]，因此轴距和轴重是交通荷载参数分析中的重点。本研究将车辆各轴按车头到车尾的顺序编号为轴1～6，通过回归分析得出车辆轴距分布以及轴重的比例分配，各类车型的轴距均值和标准差如表5所示，轴重比例分配如表6所示。

表5 轴距分布结果(单位：m)Tab.5 Result of wheelbase distribution(unit：m)

表6 轴重比例分配Tab.6 Proportion of axle load distribution

1.6 车时距

车时距是车辆依次通过监测断面的时间间隔，是随机车流模拟中的重要参数。根据WIM系统的监测数据得到监测时间段内平均时间间隔为51.2 s，车辆较为稀疏。同样对时间间隔分布拟合和K-S检验，得到车时距服从对数正态分布，其中μ=3.59，σ=1.42。

2 精细化的元胞自动机模型

元胞自动机作为一种高效的动态模型，由元胞空间、元胞、演化规则、相邻元胞4个主要组成部分[14]。首先通过一定的规则、尺寸将模拟车道离散为许多个网格，即为元胞，每个元胞都有各自的状态参数取值；由元胞和明确的边界条件共同组成了元胞空间；演化规则则是对下一时刻的元胞状态参数进行更新的具体方法；相邻元胞[8]是指能对某特定元胞的演化规则产生影响，为该规则影响范围内的元胞集合。因此，元胞自动机就是具有各自状态参数的元胞，根据自身状态与其相邻元胞的状态，通过预先设置的演化规则，以一定的更新步长在模拟空间中自行演变，从而能够实现对交通流的准确模拟。

2.1 元胞长度和更新步长

在元胞自动机模型中，元胞长度是指元胞网格划分的尺寸大小，在交通状况的模拟和加载精度两方面，元胞长度起着决定性作用。首先，元胞长度的取值决定了模型中车速的连续程度和车流行驶的稳定性[14]。元胞长度的取值越大，则车速取值就越离散，致使车辆的加速度偏大，这样会导致较大误差。另一方面，元胞长度也直接决定了影响线加载的精度。采用多个元胞模拟一辆车时，车辆的轴距和轴重分布受控于元胞长度，合理的长度取值可使模型中车辆的轴距和轴重和实际情况相对应。更新步长的设置与驾驶员的反应时间有关。敬明[15]等人在对元胞自动机模型更新步长和元胞尺寸的研究表明：小的更新步长可以较好体现车辆的相对运动，车辆可以及时针对周围环境的变化做出反应，减少了行驶过程中为防止碰撞所产生的减速，车辆行驶平稳。

因此，本研究改变传统元胞自动机模型一般用一个长5 m或7.5 m的元胞代表一辆车的方法，采用多个元胞模拟一辆车，将元胞长度和更新步长设置为可调节的变量，根据具体情况来模拟实际交通状况：对于中小跨径的桥梁，可减小元胞长度提高模拟精度，反之增大元胞长度，同时在相应位置的元胞一端用集中荷载表示轴重。

2.2 换道规则

对于同向双车道交通，车辆在行驶过程中可在本车道继续行驶，也可以根据周围车辆位置的变化和自身行驶状态[6]，在满足换道条件下可选择另一个车道行驶，各车道的换道规则如下[16]。

图1 行车道向超车道换道示意图Fig.1 Schematic diagram of lane changing from carriageway to overtaking lane

(1)

同时，车辆满足上述换道条件后还需根据驾驶员的主观意愿决定最终是否换道，故引入行车道向超车道换道概率Pchange21来表征换道的随机性，rand

图2为从超车道向行车道换道示意，同理若超车道上车辆与前车间距gap1小于车辆的最大行驶速度Vmax和安全行车距离dsafe-len的最大值，且超车道上须换道车辆与行车道前方车辆间距gap21+大于gap1，与行车道后方车辆间距gap21-大于Vmax和dsafe-len的最大值，则超车道车辆满足换道条件，会考虑向行车道换道行驶。

图2 超车道向行车道换道示意图Fig.2 Schematic diagram of lane changing from overtaking lane to carriageway

其次，若超车道上须换道车辆与行车道前方车辆间距gap21+大于Vmax和dsafe-len的最大值，超车道上车辆也可换向行车道行驶。

(2)

同上，引入超车道向行车道换道概率Pchange12，rand

2.3 行驶规则

(3)匀速规则：车辆保持原来速度继续前行，车速更新为：Vi(t+dt)=Vi(t)；

(4)随机慢化：为体现车辆行驶的随机性，引入随机慢化概率：rand

2.4 边界条件

在元胞自动机模型中，车辆进出模拟车道需设置边界条件，包括周期性和开放性的边界条件。周期性的边界条件指车辆驶出后马上又会进入到模型中，这样会保持车辆数和车流的密度不变，对于实际交通状况的模拟有较大误差[18]。开放性的边界条件则是用一定概率α来产生车辆，和一定概率β来离开[19]，更能体现实际交通状况。本研究用开放性的边界条件，车辆产生概率α根据实测的车时距取值，对于车辆离开概率β设置为全开放边界。

3 随机车流模拟及荷载效应分析

3.1 随机车流建模过程

在实测交通荷载参数和元胞自动机模型的基础上，对某60 m跨径的单箱三室简支钢箱梁桥进行随机车流模拟。该桥位于云南某低等级公路，单向双车道，主梁断面如图3所示。主梁梁高3 m，顶板宽10.5 m，厚16 mm，设置8 mm厚的闭口加劲肋和14 mm 厚的板肋；底板宽7.5 m，厚14 mm，设置190 mm×16 mm的板肋；腹板厚16 mm，设置160 mm×14 mm 的板肋。

图3 钢箱梁横截面(单位：mm)Fig.3 Cross-section of steel box girder (unit:mm)

在元胞自动机模型中，车辆的轴距和轴重等参数可准确确定，而随着车辆位置的不断更新，桥梁所产生的荷载效应也在不断变化。影响线是获取桥梁结构在移动荷载作用下结构响应的基本工具，因此，根据随机车流和结构的影响线，利用叠加原理便可得出桥梁的荷载效应，其原理如下所示[20]：

(3)

建模流程如图4所示：首先，根据模拟精度和车轴信息将桥梁车道离散得到元胞空间，并且将WIM系统的监测数据分别拟合为相应的数学概率分布模型，得到车辆相关特征参数，包括车辆类型、车重、车速、车长、轴重和轴距、车时距。其次，基于实测的交通特征参数生成车辆序列，车辆按照设定的边界条件进入模拟车道，并根据实测的车速分布赋予车速值，通过换道规则和行驶规则在元胞空间随时间自行演化。最后，根据影响线加载方式，采用叠加原理计算得到模拟车流作用下桥梁结构的荷载响应，直至模拟总时长完成，得出该桥梁微观交通流及其随机荷载应力历程。

图4 随机车流模拟流程Fig.4 Process of random traffic flow simulation

根据具体交通荷载参数和桥梁信息，设定如表7所示的模型参数。首先结合桥梁跨径和模拟精度，元胞长度取0.1 m可以准确模拟出车辆的轴距分布；另外通过该桥梁平均车速为36.6 km/h，即10.2 m/s，则更新步长取0.2 s时，相当于车辆前进2 m统计一次影响线加载的效应，这样可达到模拟精度的要求。元胞自动机模型会根据车辆最大速度和周围环境自动得出每辆车的行驶速度[16]，因此模型需要运行一段时间后，交通流可达到稳定状态，以使结果准确可靠。双车道元胞自动机模型达到稳定的状态所需时间为[17]：

t=模拟长度×车道数量×10/元胞长度=12 000 s。

(4)

因此本次模拟结果需以12 000 s为起点计算，即模拟时长为稳定时间加上桥梁实际监测时间。最后由于低等级公路条件限制，故换道概率及加速度等取值较小。

3.2 模拟结果及校核

为了验证随机车流模拟的准确性，利用WIM系统实测数据对模型进行校核：根据实测的基本统计参数和元胞自动机模拟结果的统计参数进行对比，从车型比例、车速、车长、车头时距几个方面进行模型校核，通过实测与模拟结果对比显示，基于元胞自动机的随机车流模拟具有较高的精度。

表7 模型参数Tab.7 Model parameters

各类车型比例的实测与模拟结果对比如表8所示，各类车型的模拟结果最大误差不超过0.07%，模拟精度较高。

表8 车型比例的实测与模拟结果对比(单位：%)Tab.8 Comparison of vehicle type proportion between measured and simulated results(unit:%)

各类车型的车速实测与模拟结果对比如表9所示，各车速均值的误差均可控制在1.5%以内，表明车速的模拟状况与实际较为符合。

表9 车速实测与模拟结果对比(单位：km/h)Tab.9 Comparison of vehicle speed between measured and simulated results(unit：km/h)

各类车型的车长实测与模拟结果对比如表10所示，第2类车长的均值误差较大，达到9.64%，第1类误差达到3.54%，其余可控制在1.4%以内，可较好地模拟实际交通状况。

表10 车长实测与模拟结果对比(单位：m)Tab.10 Comparison of vehicle length between measured and simulated results(unit：m)

图5为车头时距的分布拟合图，由实测和模拟结果对比，二者规律一致。

图5 车头时距分布拟合实测与模拟结果对比Fig.5 Comparison of time headway distribution fitting between measured and simulated results

3.3 模拟车流荷载响应

将模拟的随机车流加载到钢箱梁跨中下缘的应力影响线上，车流周期取本次监测的时间29 d，得到了模拟车流作用的荷载应力历程如图6所示。

图6 模拟车流作用应力历程Fig.6 Stress history under simulated vehicle flow

由于随机车流作用下得到的荷载应力历程是无规律的，而雨流计数法[3]可以对类似复杂无规律的波动情况进行分析，并且对出现相同幅值的次数进行累计，进而可将复杂的随机情况变成概率统计问题，所以通过雨流计数法从上述应力历程中整理得到相应的应力循环次数和幅值，最终得出应力幅分布如图7所示。

在模拟车流作用，得出了桥梁结构在29 d的车流周期下的应力幅分布情况，通过对比我国的《公路钢结构桥梁设计规范》(JTG D64—2015)[5]、英国的BS5400规范[21]和美国的AASHTO规范[22]提出用于钢结构桥梁疲劳荷载计算模型作用下的应力幅值，结果如表11所示。结果表明：我国规范计算值>英国BS5400规范计算值>美国AASHTO规范计算值，我国规范取值较为安全，然而对于云南山区低等级公路桥梁实际荷载作用下的应力幅值，规范计算结果不能准确客观地描述各个地域复杂的交通荷载状况，并不能涵盖实际情况下存在的高应力幅值，即实际状况下单次作用荷载更大，因此桥梁的疲劳研究需要根据不同的交通特性提出具体适用的疲劳荷载，则通过随机车流模拟得到该地区的疲劳荷载谱具有重要意义。

图7 应力幅分布Fig.7 Stress amplitude distribution

表11 规范计算值Tab.11 Calculation values stipulated in specifications

4 结论

(1)通过WIM系统对云南山区低等级公路交通荷载的实时监测，探明了高原山区公路桥梁实际运营期间的交通荷载特性，研究分析了各类车型的车重、车速、车长、轴距、轴重和车时距的特征参数分布规律。

(2)基于传统的元胞自动机模型，对元胞的长度和更新步长作细化研究，提出了适用于山区公路交通流模型的换道、行驶的具体规则与边界条件的设置，以及模型参数的具体取值，完成了精细化的元胞自动机模型对随机车流的模拟。最后根据WIM系统实测数据校核了元胞自动机模型，校核结果表明：元胞自动机模拟车流得到的车流参数与实测参数较好地吻合，保证了模型的可靠性，同时精细化的元胞自动机模拟结果有较高的精度，在随机车流模拟中具有一定的优势。

(3)通过随机车流加载到钢箱梁跨中下缘的应力影响线得到了相应的荷载应力历程，并采用雨流计数法统计了模拟车辆荷载作用下的应力幅值和循环次数，可为云南省这样的山区低等级公路桥梁的疲劳寿命分析和抗疲劳设计提供参考依据。