基于临界滑水速度的湿滑路面全域滑水风险量化评价

摘要

基于修正LuGre摩擦模型提出临界滑水速度的快速计算方法，以临界滑水速度为滑水风险事件发生的判断指标；基于湿滑路段水膜厚度分布以及车速、轮迹等车辆行驶状况的概率模型，计算滑水风险事件发生的概率，从而对滑水风险进行量化，并根据风险概率将滑水风险划分为5个等级。结果表明：临界滑水速度主要受水膜厚度的影响，与路面的高程特征密切相关；当路面存在车辙、坑槽等表面病害时，病害区域的临界滑水速度显著下降。车辆滑水风险与临界滑水速度分布基本一致，但轮迹带区域的风险相对较高；降雨强度较大时，无病害或低病害路段的低风险区域超过96%，路段整体处于低风险等级，车辙等严重病害处出现大量中高风险区域，导致路段整体风险等级提高。

关键词

道路安全; 临界滑水速度; 滑水风险; 量化评价

路面湿滑是导致道路交通安全事故频发的重要原因。国内外统计数据显示^［1］，道路潮湿、积水等引发的交通事故基本占到交通事故总数的12%以上。在湿滑路面上，轮胎‒路面‒水膜的相互作用产生动水压力，当车辆速度与水膜厚度超过一定阈值后，轮胎在动水压力的作用下与路面发生分离^［2］，出现滑水现象，从而导致车辆行驶稳定性下降，极易引发交通安全事故。因此，实时、准确地评价湿滑路面上的滑水风险对于保障雨天道路交通安全至关重要。

现阶段主要以临界滑水速度为指标^［3-4］判断车辆是否发生滑水事故，从而评价车辆滑水风险。在一定的水膜厚度下，使轮胎和路面发生完全分离的最小速度即为临界滑水速度。精准计算临界滑水速度是评价滑水风险的前提，目前国内外学者对临界滑水速度的求解方法主要分为3类，即现场试验法、理论解析法和数值模拟法。在现场试验方面，美国国家航空航天局（National Aeronautics and Space Administration，NASA）^［5］在20世纪60年代开展了一系列滑水试验，提出了NASA公式；此后，Gallaway等^［6］根据现场试验数据回归，建立了考虑水膜厚度、轮胎特性等参数的临界滑水速度经验公式。在理论解析方面，学者们主要通过水力学理论分析动水压力与水膜厚度、车轮胎压等影响因素之间的关系，进而求出上述因素对滑水速度的影响规律^［7-8］。随着有限元等方法的逐步完善，学者们开始通过数值模拟手段分析临界滑水速度。Fwa等^［9］、Srirangam^［10］、Anupam等^［11］以流固耦合理论为核心建立轮胎‒路面‒水膜的相互作用模型，对滑水现象进行仿真分析；国内，董斌^［12］、朱晟泽^［13］、贾兴利等^［14］通过数值仿真手段分析了水膜厚度、动水压力、道路线形等因素对临界滑水速度的影响规律。整体而言，滑水试验获取的经验公式准确性高，但是成本高，测试条件苛刻，泛用性有限；滑水速度解析计算公式具有一定的物理意义，但是对于复杂湿滑状态下的滑水特征还需要进一步研究；数值模拟能够较为准确地对滑水现象进行仿真，但是模型精细化程度的提高带来的是计算时间的成倍增加，限制了模型的实际应用。仅以单一临界滑水速度为阈值，得出车辆滑水或不滑水的二元结论，难以量化表征真实条件下路段的滑水风险，因此还需要考虑路面自身特性、车辆特性和交通特性等因素，进一步加强湿滑路面上车辆滑水风险的量化研究。

在道路智能化发展趋势下，利用道路信息的智能感知手段对道路安全风险进行精细化评价，成为保障道路安全的新思路。鉴于此，本研究基于临界滑水速度对湿滑路面上的车辆滑水风险进行了量化评价。首先，分析了车辆滑水风险的成因，以临界滑水速度作为风险事件发生的判断指标，并基于修正LuGre摩擦模型提出了临界滑水速度的快速计算方法；然后，利用路段全域水膜厚度分布和车速、轮迹等车辆行驶状况的概率分布模型，计算了路段各点滑水风险事件发生的概率，对滑水风险进行量化，并提出了滑水风险等级划分方法；最后，基于3个典型路段阐明了湿滑路段车辆滑水风险的量化评价过程。

1 车辆滑水风险的量化评价流程

1.1　车辆滑水风险的失效准则与判断方法

车辆行驶在湿滑路面上，轮胎与路面间存在一定的水膜，水膜的润滑作用和动水压力作用导致轮胎与路面的接触面积减小、附着系数下降；特定的水膜厚度下，当车速超过某一阈值时，轮胎与路面完全分离，轮胎‒路面附着系数下降到零，产生完全滑水现象，此阈值为临界滑水速度。因此，选取车速作为滑水风险的判别指标，以临界滑水速度作为判断阈值，当车速大于等于临界滑水速度时则认为车辆出现滑水风险。

1.2　车辆滑水风险的量化方法

目前对于风险的量化评价主要从事故严重程度以及发生的概率两方面进行。然而，滑水事故的严重程度受众多因素影响，只能在发生后进行量化，难以直接根据路面湿滑状态预估。因此，对于湿滑路段，利用车速、轮迹等车辆行驶状况的概率分布，计算路段面域各点实际车速超过该点临界滑水速度的概率，以此作为车辆发生滑水风险事件的概率，从而对滑水风险进行量化。

1.3　风险等级的划分方法

在对滑水风险进行评价时，需要通过一定标准对量化指标值进行划分，建立风险等级。参考 ISO/TR 24971^［15］中对于风险概率等级的划分方法，将滑水风险等级划分为5个等级，各等级对应的概率如表1所示。

表1  风险等级的划分准则

Tab.1  Rules for risk level

风险等级	概率
A（高风险）	≥10-3
B（中高风险）	［10-4，10-3）
C（中风险）	［10-5，10-4）
D（中低风险）	［10-6，10-5）
E（低风险）	<10-6

1.4　车辆滑水风险的评价流程

湿滑路面全域滑水风险的量化评价主要分为两部分：一是获取发生滑水风险事件的判别指标，即计算特定路段湿滑状态下的临界滑水速度分布；二是量化湿滑路段全域滑水风险，即计算当车辆行至路面某点时滑水风险事件的发生概率，并依据风险等级划分方法对滑水风险等级进行评价。具体的评价流程如图1所示。

图1  湿滑路面车辆滑水风险评价流程

Fig.1  Flowchart of vehicle hydroplaning risk evaluation on wet pavement

基于修正LuGre摩擦模型对不同湿滑状态下车辆的临界滑水速度进行快速求解。然后，在已知路表全域高程信息以及降雨强度的基础上，通过二维浅水方程获取路段全域水膜厚度分布，以表征路段湿滑状态。根据临界滑水速度计算方法，路段各点的临界滑水速度可以通过该处的水膜厚度进行计算，因此获取湿滑路段全域的临界滑水速度分布。

在此基础上，考虑到车辆行驶状态的随机性，基于蒙特卡罗法，利用车速、轮迹等车辆行驶状况的概率分布模型，模拟车辆在湿滑路面上行驶的随机位置及随机速度，计算车辆行驶在路段各点时发生滑水风险事件的概率。最后，通过滑水风险等级划分方法获取湿滑路段全域的滑水风险等级分布，并对路段整体的滑水风险进行评价。

2 车辆临界滑水速度的计算方法

2.1　修正LuGre摩擦模型

LuGre摩擦模型是由Wit等^［16］提出的半理论摩擦模型。根据LuGre摩擦模型，可将轮胎与路面的接触界面等效为若干个微小、无质量的刷毛，其中路面刷毛刚性且固定，胎面刷毛为弹性，各接触点所受到的摩擦力可通过刷毛的瞬时变形以及对时间的导数计算。

基于Zhao等的研究^［17］，可以通过纹理影响系数$\theta$、水膜修正系数Y_R和Y_F来表征路面纹理及水膜厚度对轮胎‒路面附着系数的影响。$\theta$为标准路面纹理下轮胎‒路面真实接触面积与当前纹理下的比值；Y_R为水膜作用下轮胎‒路面接触区域长度与干燥情况下的比值，Y_F为向上的动水压力与干燥情况下轮胎竖向力的比值。因此，修正LuGre摩擦模型可以表示为

$$\frac{\partial \overline{z}}{\partial t}=v_{\mathrm{r}}-\theta Y_{\mathrm{R}}\left[\frac{{\sigma }_{\mathrm{0}}\left|v_{\mathrm{r}}\right|}{g(v_{\mathrm{r}})}\overline{z}(t)+K(t)\left|wr\right|\overline{z}(t)\right]$$

（1）

$$\mu =\left(\theta Y_{\mathrm{R}}-Y_{\mathrm{F}}\right)\left[\theta Y_{\mathrm{R}}{\sigma }_{\mathrm{0}}\overline{z}(t)+{\sigma }_{\mathrm{1}}\frac{\partial \overline{z}(t)}{\partial t}+{\sigma }_{\mathrm{2}}{v}_{\mathrm{r}}\right]$$

（2）

$$g(v_{\mathrm{r}})={\mu }_{\mathrm{c}}+({\mu }_{\mathrm{s}}-{\mu }_{\mathrm{c}}){\mathrm{e}}^{-{\left|v_{\mathrm{r}}/v_{\mathrm{s}}\right|}^{\mathrm{0.5}}}$$

（3）

$$\theta =\frac{A_{\mathrm{c}\mathrm{0}}}{A_{\mathrm{c}}}$$

（4）

$$v_{\mathrm{s}}={\beta }_{\mathrm{1}}{\mathrm{e}}^{\mathrm{1}\mathrm{000}{\beta }_{\mathrm{2}}h+{\beta }_{\mathrm{3}}}+{\beta }_{\mathrm{4}}$$

（5）

$$Y_{\mathrm{R}}=\mathrm{1}-\frac{\mathrm{12}v\gamma r_{\mathrm{0}}^{\mathrm{2}}}{a\mathrm{\pi }pL}\left(\frac{\mathrm{1}}{h_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{2}}+\mathrm{2}\epsilon h_{\mathrm{m}\mathrm{i}\mathrm{n}}}-\frac{\mathrm{1}}{h_{\mathrm{0}}^{\mathrm{2}}+\mathrm{2}\epsilon h_{\mathrm{0}}}\right)$$

（6）

$$Y_{\mathrm{F}}=\frac{\rho w_{\mathrm{t}}r{v}^{\mathrm{2}}}{\mathrm{3}{F}_{\mathrm{n}}}\left(\sqrt[]{{\left(\frac{L}{\mathrm{2}r}\right)}^{\mathrm{2}}-{\left(\frac{h}{r}\right)}^{\mathrm{2}}+\frac{\mathrm{2}h}{r}\sqrt[]{\mathrm{1}-{\left(\frac{L}{\mathrm{2}r}\right)}^{\mathrm{2}}}}-\frac{L}{\mathrm{2}r}\right)$$

（7）

式中：$\overline{z}$为轮印中刷毛的平均变形量；v为轮胎运动速度；r为轮胎半径；w为轮胎转动角速度；v_r为轮胎的相对速度，v_r=wr-v；v_s为Stribeck速度；μ为轮胎与路面间的附着系数；μ_c为库仑摩擦系数；μ_s为最大静摩擦系数；F_n为轮胎受到的竖向力；L为轮印长度； K（t）为竖向力分布影响的函数，可近似简化为 K（t）=7/6L；A_c为当前纹理下轮胎‒路面真实接触面积，A_c0为标准路面纹理下轮胎‒路面真实接触面积；γ为水的黏度；ρ为水的密度；r₀为轮印等效半径；p为胎压；h为水膜厚度；h_min为轮胎与路面间的最小水膜厚度；h₀为初始水膜厚度（为简化计算，取h₀=0.01h）；$\epsilon$为沥青路面平均断面深度；w_t为轮胎宽度；a为轮印系数，对于矩形轮印，取18.1；σ₀、σ₁、σ₂、β₁、β₂、β₃、β₄为模型的标定参数。

采用GripTester摩擦测试车的实测数据测试不同车速和水膜厚度下的摩擦系数，通过非线性回归方法对模型参数进行标定。最终得到的标定参数以及其他汽车计算参数如表2所示。

表2  LuGre摩擦模型的计算参数

Tab.2  Parameters of LuGre friction model

标定参数	取值	汽车参数	取值
σ0	481.69	Fn/N	3 469
σ2	0	p/Pa	200 000
μc	0.63	wt/m	0.205
μs	1.24	L/m	0.164
β1	1.46	b/m	0.125
β2	-0.69	r/m	0.316
β3	0.57	r0/m	0.103
β4	5.24

注：  b为轮印宽度。

2.2　轮胎临界滑水速度计算方法

修正LuGre摩擦模型能够较为准确地建立车速、水膜厚度与附着系数之间的映射关系，如图2所示。对于特定水膜厚度，当车速增加至某一临界值时，附着系数为零，轮胎与路面分离，将此时的速度临界值视为当前条件下的临界滑水速度。

图2  临界滑水速度计算

Fig.2  Calculation of critical hydroplaning velocity

单个轮胎临界滑水速度计算方法如下：在已知水膜厚度及其他参数下，按照0.1 km·h^-1的计算间隔，计算车速在10 km·h^-1至150 km·h^-1范围内轮 胎‒路面间的附着系数，附着系数为零的最小车速即为当前条件下的轮胎临界滑水速度。

计算不同水膜厚度下轮胎的临界滑水速度，并将计算结果与李强等的临界滑水速度计算模型^［8］进行对比，如图3所示。结果表明，本研究的临界滑水速度计算模型与李强等的临界滑水速度计算结果^［8］具有极高的一致性，两者相对误差小于5%。因此，可认为通过修正LuGre摩擦模型得到的轮胎临界滑水速度是相对可靠的。

图3  不同水膜厚度下轮胎的临界滑水速度

Fig.3  Critical hydroplaning velocity under different water film thicknesses

2.3　车辆临界滑水速度计算方法

前文所提的临界滑水速度模型是针对单个轮胎进行计算，考虑到车辆实际行驶过程中水膜厚度的不均匀分布，4个轮胎可能对应着不同的临界滑水速度。因此，将汽车重心作为控制点，当车辆重心行驶至某点时，获取4个轮胎所在位置的水膜厚度，分别计算对应的临界滑水速度，取最小值作为汽车重心到达该位置时的临界滑水速度，如图4所示。同时，为保持汽车4个轮胎均在路域范围内，汽车的重心位置应控制在图中的白色区域，此时4个轮胎均位于路域范围内；当汽车重心位于花纹区域时，轮胎超出路域范围，因此不考虑汽车重心位于花纹区域的情况。

图4  车辆临界滑水速度计算方法

Fig.4  Calculation method of vehicle critical hydroplaning velocity

3 典型路段滑水风险量化评价

3.1　典型路段临界滑水速度分布

为说明所提出的滑水风险量化评价方法的合理性，选取3条典型路段进行滑水风险分析。利用车载LiDAR移动测绘系统获取路段全域的高程模型。各路段的高程特征如图5所示。路段1为具有一定横坡和纵坡且无明显病害的普通道路，横坡、纵坡的坡度均约为1%；路段2为具有局部坑槽的拱形路面，横坡坡度约为2.5%，纵坡则在纵断面方向约40 m处存在变坡点，变坡点前后的纵坡坡度分别约为0和0.5%；路段3为存在严重车辙的道路，横坡坡度约为1%，纵坡坡度约为3%，最大车辙深度为 15 mm。

图5  典型路段的高程特征

Fig.5  Elevation characteristics of typical road sections

基于前期研究^［18］，采用二维浅水方程对水膜厚度面域分布进行预估，采用该模型可计算任意路面高程下、不同降雨强度时道路全域的水膜厚度时空分布，其控制方程如下所示：

$$\frac{\partial h}{\partial t}+\frac{\partial (u_{x}h)}{\partial x}+\frac{\partial (u_{y}h)}{\partial y}=q_{\mathrm{r}}(t)$$

（8）

$$\frac{\partial (u_{x}h)}{\partial t}+\frac{\partial \left(u_x^{\mathrm{2}}h+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\right)}{\partial x}+\frac{\partial (u_x{u}_{y}h)}{\partial y}=$$

$$gh{S}_{\mathrm{0},x}-gh{S}_{\mathrm{f},x}$$

（9）

$$\frac{\partial (u_{y}h)}{\partial t}+\frac{\partial (u_x{u}_{y}h)}{\partial x}+\frac{\partial \left(u_y^{\mathrm{2}}h+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\right)}{\partial y}=$$

$$gh{S}_{\mathrm{0},y}-gh{S}_{\mathrm{f},y}$$

（10）

式中：u_x和u_y分别为水流速度沿x和y方向的分量； q_r（t）为降雨强度；S_f为摩擦源项，S₀为底坡源项，可根据路面高程数据计算；g为重力加速度。

以降雨强度3 mm·min^-1、持续降雨180 s时的湿滑路面作为滑水风险分析对象，典型路段的水膜厚度分布如图6所示。3个路段的水膜厚度分布呈现道路中心线向路肩方向逐渐增加的趋势。路段2圆圈部分表示局部坑槽导致路面径流累积，此处的水膜厚度高达9.87 mm。从图6c可以看出，路段3大量路面径流在车辙凹陷处汇集，从而在轮迹带处形成长条的积水区域，水膜厚度最高可达22.77 mm，远超无病害路面。

图6  典型路段的水膜厚度分布

Fig.6  Water film thickness distribution of typical road sections

利用车辆临界滑水速度计算方法，在典型路段水膜厚度分布的基础上，计算得到典型路段的临界滑水速度分布，如图7所示。3个典型路段临界滑水速度的特征值如表3所示。

图7  典型路段的临界滑水速度分布

Fig.7  Critical hydroplaning velocity distribution of typical road sections

表3  典型路段的临界滑水速度特征值

Tab.3  Characteristic values of critical hydroplaning velocity for typical road sections

路段	临界滑水速度/（km·h-1）
	最小值	最大值	平均值
路段1	93.0	121.8	106.7
路段2	85.4	137.2	112.2
路段3	70.3	133.9	97.9

对于如路段1所示的普通道路，其临界滑水速度从道路中心线向两边逐渐降低；在路面坡度的排水作用下，即便是降雨强度为3 mm·min^-1的特大暴雨场景，路段全域内临界滑水速度仍保持在较高水平，最小值为93.0 km·h^-1，比路段限速80 km·h^-1高出16.3%。对于如路段2所示的路拱路段，其临界滑水速度呈现道路中心线高、两边低的特点；由于其表面存在局部坑槽，导致路段局部区域临界滑水速度降低，临界滑水速度的最小值为85.4 km·h^-1，仅比路段限速80 km·h^-1高出6.6%；除去坑槽区域后，临界滑水速度的最小值提高至93.4 km·h^-1，与路段1的结果相近。对于如路段3所示的严重车辙路段，其临界滑水速度显著降低，当汽车轮胎处在车辙位置时，临界滑水速度的最小值可低至70.3 km·h^-1，仅为路段限速80 km·h^-1的87.9%；车辆在该路段行驶时若不进行合理的车速控制，极易出现滑水风险。

3.2　典型路段车辆滑水风险等级评价

在路段临界滑水速度分布的基础上，利用车速、轮迹的概率分布模型，通过蒙特卡罗法模拟车辆在湿滑路段上行驶的随机位置及随机速度，从而计算路段各点实际车速超过该点临界滑水速度的概率，对各点的车辆滑水风险进行量化。

本研究采用Vadeby等^［19］基于实测数据建立的车速分布模型。该模型表明，当路段限速为80 km·h^-1时路段内行驶车辆的速度分布满足均值为75.3 km·h^-1、标准差为5.7 km·h^-1的正态分布。

对于车道内车辆的轮迹分布，基于董忠红等^［20］研究中的实测数据，采用如图8所示的驼峰型轮迹分布进行计算。利用上述车速、轮迹的概率分布，采用蒙特卡罗法随机生成200 000组工况，计算路段各点实际车速超过临界滑水速度的概率，得到典型路段的风险等级分布，如图9所示，分布特征如表4所示。

图8  道路轮迹的横向分布

Fig.8  Lateral distribution of wheel-path

图9  典型路段的滑水风险等级分布

Fig.9  Hydroplaning risk level distribution of typical road sections

表4  典型路段的滑水风险等级分布特征

Tab.4  Characteristics of hydroplaning risk level distribution for typical road sections

路段	等级A占比/%	等级B占比/%	等级C占比/%	等级D占比/%	等级E占比/%	平均风险等级
路段1	0	0	0.08	3.45	96.47	D
路段2	0	0.11	0.32	0.59	98.98	C
路段3	12.47	13.63	9.88	9.13	54.89	A

3个路段滑水风险等级分布特征与其临界滑水速度分布特征基本一致，但轮迹带处的滑水风险等级相对提高。对比路段1和2可知，两者处于风险等级C、D、E区域的占比接近，但由于路段2中存在零星的路面坑槽，导致路段2局部的风险等级可达到B级，零星高危区域的存在使路段2的平均风险等级达到C级，高于路段1。由路段3的结果可知，车辙的存在显著增加高风险等级的占比，车辙深度越深，车辆行驶越容易面临滑水风险。路段3中风险等级A、B的占比高达26.10%，主要位于车辙区域，致使全域的平均风险等级达到了A级。总体而言，降雨强度较大时，无病害或低病害路段全域基本处于低滑水风险等级，只有在坑槽等少量病害区域出现了中风险；对于车辙等严重的路面病害，一旦出现便会提高全域的滑水风险等级，若不及时进行路面状况修复以及合理的车速控制，路段内极易因车轮滑水频发而出现交通安全事故。

4 结论

（1）提出了湿滑路面车辆滑水风险的量化评价流程。以临界滑水速度为风险事件发生的判断指标，考虑路面湿滑状态以及车速、轮迹等车辆行驶状况，计算风险事件发生的概率，从而对滑水风险进行量化与等级评价。

（2）建立了车辆临界滑水速度的快速计算方法。在LuGre摩擦模型基础上，引入路面纹理及水膜厚度的影响，建立了不同车辆运动参数及水膜厚度与附着系数的映射关系，以附着系数为零时的最小车速为临界滑水速度。

（3）临界滑水速度主要受水膜厚度的影响，与路面的高程特征密切相关。车辆滑水风险与临界滑水速度分布基本一致，且轮迹带区域的滑水风险等级相对较高。当路面存在车辙、坑槽等表面病害时，病害区域的滑水风险显著提高。

利用道路信息感知手段，考虑路面湿滑状态，以及车速、轮迹等车辆行驶状况的概率分布，实现了湿滑路面上车辆滑水风险量化评价。

作者贡献声明

赵鸿铎：设计研究方法，统筹研究工作。

赵兰若：数据分析，论文写作。

蔡爵威：模型构建，算法设计。

马鲁宽：参与方法研究，修改论文。