制动防抱死系统是汽车行驶过程中必不可少的主动安全系统之一，而如何精准预估路面附着系数是车辆实施理想防抱死控制的前提条件之一^[1]. 现有路面附着系数估计方法主要分为Cause-based识别和Effect-Based识别方法2大类^[2]. Cause-based识别方法分析与路面附着系数相关的各种物理因素，如利用光学传感器^[3]或者激光传感器^[4]，建立反映各个因素与路面附着系数关系的数学模型，通过测量相关因素，利用所设计的数学模型计算路面附着系数. 其成本较高，硬件复杂，商业推广价值不大，且识别精度过于依赖传感器信号的损失程度，在面对崎岖不平路面、行驶信号损耗严重的情况下，预估准确率不高^[5]. Effect-Based识别方法通过分析路面变化引起的车辆动力学响应来估算路面附着系数. 该类方法对工作环境要求不高，仅对车辆进行动力学分析即可进行路面附着系数辨识. 其包括多种技术途径，如卡尔曼滤波估计器^[6-7]、滑模观测器^[8]、神经网络^[9-10]等.

相比于卡尔曼滤波估计器和滑模观测器，神经网络处理非线性关系的能力更强，具有更加强大的容错性. Park等^[11]采用具有传感器特征选择的深度集成网络方法在极短时间内判别路面类型；Sun等^[12]提出遗传算法改进神经网络，利用路面纹理三维数据预估路面附着系数；林棻等^[13]设计改进Keras模型路面附着系数估计器，极大提高了路面附着系数预估精度. 上述方法都能实现基本的估计功能，但是均无法及时感知车辆在路面较大突变情况下的路面附着系数变化，若在2类路面交界处发生紧急制动，则可能因为路面附着系数预估不准确而造成防抱死控制效果变差，甚至引发车轮抱死危险. Elman神经网络被认为是典型的动态递归神经网络，其适应时变特性能力较强，能够高效处理非线性影响下的不确定性，该类方法在路面附着系数估计领域具有更大的发展潜力^[14]. 伍文广等 ^[15]引入Elman神经网络建立路面附着系数估计模型，与传统BP网络相比，提高了预估精度. 然而，传统Elman神经网络算法容易陷入局部极小点^[16]，导致路面附着系数预估不准确而引起车辆滑移率控制效果不佳，因此该方法在路面附着系数辨识应用中仍有较大提升空间.

在汽车行驶状态和路面工况时变条件下，为了克服采用传统Elman神经网络的路面附着系数估计方法可能引起路面附着系数预估不准确，进而导致车辆滑移率控制效果不佳的缺陷，设计基于变权重PSO-Elman神经网络的路面附着系数估计方法，提出将传统Elman神经网络的权值矩阵作为初始粒子群位置，以神经网络绝对误差和作为适应度值来调整粒子群位置，达到神经网络权值矩阵最优效果，以相比于传统Elman神经网络提高权值矩阵更新稳定性.

若不考虑悬架影响，构建如图1所示动力学模型^[17]，x、y、z轴的方向分别为行驶过程中汽车的纵向、横向和垂向.

汽车动力学方程^[18]如下 ：

(1)$ \begin{split} m\left( {{{\dot v}_x} - {v_y}\dot \gamma } \right) =& - \left( {{F_{y1}}+{F_{y2}}} \right)\sin \; \delta - {F_{\text{D}}}\cos \; \delta - \qquad\\ &\left( {{F_{x1}}+{F_{x2}}+{F_{{\mathrm{f}}1}}+{F_{{\text{f2}}}}} \right)\cos \; \delta - \\ &{F_{x3}} - {F_{x4}} - {F_{{\text{f3}}}} - {F_{{\text{f4}}}},\end{split} $

(2)$ \begin{split} m\left( {{{\dot v}_y}+{v_x}\dot \gamma } \right) =& \left( {{F_{y1}}+{F_{y2}}} \right)\cos \; \delta +{F_{y3}}+{F_{y4}} - \\&\left( {{F_{x1}} + {F_{x2}} + {F_{{\text{f}}1}} + {F_{{\text{f2}}}}} \right)\sin \; \delta - {F_{\text{D}}}\sin \; \delta ,\end{split} $

(3)$ \begin{split} {I_{{z}}}\frac{{{{\text{d}}^2}\gamma }}{{{\text{d}}{t^2}}} =&{l_{\text{a}}}\left[ \left( {{F_{y1}}+{F_{y2}}} \right)\cos \; \delta - \right. \\&\left. \left( {{F_{x1}}+{F_{x2}}+{F_{{\text{f1}}}}+{F_{{\text{f2}}}}} \right)\sin \; \delta \right] - \\ &{l_{\text{b}}}\left( {{F_{y3}}+{F_{y4}}} \right) - {t_{\text{r}}}\left( {{F_{x4}}+{F_{{\text{f4}}}} - {F_{x3}} - {F_{{\text{f3}}}}} \right)/2 - \\ &{t_{\text{f}}}\left[ \left( {{F_{x2}}\cos \; \delta +{F_{{\text{f2}}}}\cos \; \delta +{F_{y2}}\sin \; \delta } \right) -\right. \\& \left. \left( {{F_{x1}}\cos \; \delta +{F_{{\text{f1}}}}\cos \; \delta +{F_{y1}}\sin \; \delta } \right) \right]/2, \end{split} $

(4)$ {J_{{\text{W}}i}}{\dot \omega _i} = {F_{xi}}{R_{\text{ω }}} - {T_{{\text{b}}i}} - {T_{{\text{f}}i}}, $

(5)$ {F_{\text{D}}} = {{{C_{\text{D}}}Av}}/{{21.15}}. $

式中：$ m $为整车质量；$ v $为车辆速度；$ {v_x} $、$ {v_y} $为车辆沿x、y轴的速度；$ {\dot v_x} $、$ {\dot v_y} $为车辆沿x、y轴的加速度；$ \gamma $、$ \dot \gamma $为车辆横摆角、横摆角速度；$ {I_{{z}}} $为转动惯量；$ {l_{\text{a}}} $、$ {l_{\text{b}}} $为车辆质心到前后轴距离；$ {t_{\text{f}}} $、$ {t_{\text{r}}} $为前后轴距；$ {J_{{\text{W}}i}} $为车轮惯性力矩；$ {\omega _i} $和$ {\dot \omega _i} $为角速度和角加速度；$ {F_{\text{D}}} $为空气阻力；$ {T_{{\text{b}}i}} $和$ {T_{{\text{f}}i}} $为车轮制动力矩和滚动阻力力矩；$ {R_{\text{ω }}} $为车轮半径；$ {F_{xi}} $、$ {F_{yi}} $为车轮纵向力、横向力；$ {C_{\text{D}}} $为空气阻力系数；$ A $为车辆迎风面积；$ {F_{{\text{f}}i}} $为滚动阻力；$ \delta $为前轮转向角；下标i可取1、2、3、4，分别表示汽车左前、右前、左后和右后轮.

各轮胎垂直载荷满足：

(6)$ \left. \begin{gathered} {F_{z1}} = \frac{{{l_{\text{b}}}mg}}{{2\left( {{l_{\text{a}}}+{l_{\text{b}}}} \right)}} - \frac{{m{a_x}{h_{\text{g}}}}}{{2\left( {{l_{\text{a}}}+{l_{\text{b}}}} \right)}} - \frac{{m{a_y}{h_{\text{g}}}{l_{\text{b}}}}}{{({l_{\text{a}}}+{l_{\text{b}}}){t_{\text{f}}}}}, \\ {F_{z2}} = \frac{{{l_{\text{b}}}mg}}{{2\left( {{l_{\text{a}}}+{l_{\text{b}}}} \right)}} - \frac{{m{a_x}{h_{\text{g}}}}}{{2\left( {{l_{\text{a}}}+{l_{\text{b}}}} \right)}}+\frac{{m{a_y}{h_{\text{g}}}{l_{\text{b}}}}}{{({l_{\text{a}}}+{l_{\text{b}}}){t_{\text{f}}}}}, \\ {F_{z3}} = \frac{{{l_{\text{a}}}mg}}{{2\left( {{l_{\text{a}}}+{l_{\text{b}}}} \right)}}+\frac{{m{a_x}{h_{\text{g}}}}}{{2\left( {{l_{\text{a}}}+{l_{\text{b}}}} \right)}} - \frac{{m{a_y}{h_{\text{g}}}{l_{\text{a}}}}}{{({l_{\text{a}}}+{l_{\text{b}}}){t_{\text{f}}}}}, \\ {F_{z4}} = \frac{{{l_{\text{a}}}mg}}{{2\left( {{l_{\text{a}}}+{l_{\text{b}}}} \right)}}+\frac{{m{a_x}{h_{\text{g}}}}}{{2\left( {{l_{\text{a}}}+{l_{\text{b}}}} \right)}}+\frac{{m{a_y}{h_{\text{g}}}{l_{\text{a}}}}}{{({l_{\text{a}}}+{l_{\text{b}}}){t_{\text{f}}}}}. \\ \end{gathered} \right\} $

(7)$ {a_x} = {{\dot v}_x} - {v_y}\dot \gamma ,\quad {a_y} = {{\dot v}_y}+{v_x}\dot \gamma . $

式中：$ g $为重力加速度，$ {h_{\text{g}}} $为汽车质心高度，$ {a_x} $、$ {a_y} $分别为质心加速度沿x、y轴的分量.

各车轮纵向滑移率^[19]如下：

(8)$ {\lambda _i} = \left({{{v_x} - {R_{\text{ω }}}{\omega _i}}}\right)/{{{v_x}}}. $

各车轮侧偏角计算公式如下：

(9)$ \left. \begin{gathered} {\alpha _1} = \delta - {\text{arctan}}\;\; \frac{{{v_y}+{l_{\text{a}}}\dot \gamma }}{{{v_x}+{t_{\text{f}}}\dot \gamma /2}}, \\ {\alpha _2} = \delta - {\text{arctan}}\;\;\frac{{{v_y}+{l_{\text{a}}}\dot \gamma }}{{{v_x} - {t_{\text{f}}}\dot \gamma /2}}, \\ {\alpha _3} = {\text{arctan}}\;\;\frac{{{l_{\text{b}}}\dot \gamma - {v_y}}}{{{v_x}+{t_{\text{r}}}\dot \gamma /2}}, \\ {\alpha _4} = {\text{arctan}}\;\;\frac{{{l_{\text{b}}}\dot \gamma - {v_y}}}{{{v_x} - {t_{\text{r}}}\dot \gamma /2}}. \\ \end{gathered} \right\} $

采用Dugoff轮胎模型描述轮胎摩擦力^[20]：

(10)$ \left. \begin{gathered} {F_{xi}} = {C_\sigma }\frac{{{\lambda _i}}}{{1+{\lambda _i}}}f\left( {{\varepsilon _i}} \right), \\ {F_{yi}} = {C_\alpha }\frac{{\tan \; {\lambda _i}}}{{1+{\lambda _i}}}f\left( {{\varepsilon _i}} \right). \\ \end{gathered} \right\} $

(11)$ f\left({\varepsilon }_{i}\right)=\left\{\begin{array}{*{20}{l}}\left(2-{\varepsilon }_{i}\right){\varepsilon }_{i},& {\varepsilon }_{i} \leqslant 1;\\ 1,& {\varepsilon }_{i} > 1. \end{array}\right. $

(12)$ {\varepsilon _i} = \frac{{\mu {F_{zi}}\left( {1+{\alpha _i}} \right)}}{{2\sqrt {C_\sigma ^2 \lambda _i^{\text{2}}+C_\alpha ^2 {{\tan }^2}\;{\alpha _i}} }}. $

式中：$ {C_\sigma } $、$ {C_\alpha } $为轮胎纵向刚度、侧偏刚度；$ \mu $为路面附着系数；$ {\varepsilon }_{i}＞1 $表示各个车轮居于线性状态，$ {\varepsilon _i} \leqslant 1 $表示各个车轮居于非线性状态.

神经网络输入与输出应当满足一定相关性. 路面附着系数与汽车动力学模型中多种变量联系性较强，将式(10)~(12)代入式(1)~(5)，获取与路面附着系数密切相关的车辆运动状态参数. 表达式如下：

(13)$ \left. \begin{gathered} m{a_x} = \sum\limits_{i = 1}^4 {\left[ {\left( {{A_i} - {C_i}} \right)f\left( {{\varepsilon _i}} \right)} \right] - B} , \\ m{a_y} = \sum\limits_{i = 1}^4 {\left[ {\left( {{A_i} - {E_i}} \right)f\left( {{\varepsilon _i}} \right)} \right] - D.} \\ \end{gathered} \right\} $

(14)$ {A_i} = {C_{\sigma }}\frac{{{\lambda _i}}}{{1+{\lambda _i}}}. $

(15)$ B = \sum\limits_{i = 1}^4 {{F_{{\text{f}}i}}\cos \; \delta } +{F_{\text{D}}}\cos \; \delta . $

(16)$ {C_i} = {C_\alpha }\frac{{\tan \; {\lambda _i}}}{{1+{\lambda _i}}}{\varepsilon _i}\sin \; \delta . $

(17)$ D = \sum\limits_{i = 1}^4 {{F_{{\text{f}}i}}\sin \; \delta } +{F_{\text{D}}}\sin \; \delta . $

(18)$ {E_i} = {C_\alpha }\frac{{\tan \; {\lambda _i}}}{{1+{\lambda _i}}}{\varepsilon _i}\sin \; \delta . $

由于车轮非线性状态下的分析对结果的影响相对较小，在不考虑车轮非线性状态下，将式(11)、(12)代入式(13)：

(19)$ \left. \begin{gathered} {H_1}\left( \mu \right) = m{a_x} = \sum\limits_{i = 1}^4 {\left[ {\left( {{A_i} - {C_i}} \right)\left( {2 - \mu {G_i}} \right)\mu {G_i}} \right] - B,} \\ {H_2}\left( \mu \right) = m{a_y} = \sum\limits_{i = 1}^4 {\left[ {\left( {{A_i} - {E_i}} \right)\left( {2 - \mu {G_i}} \right)\mu {G_i}} \right] - D.} \\ \end{gathered} \right\} $

(20)$ {G_i} = \frac{{{F_{zi}}\left( {1+{\alpha _i}} \right)}}{{2\sqrt {C_\sigma ^2 \lambda _i^2+C_\alpha ^2 {{\tan }^2}\;{\alpha _i}} }}. $

可以看出，$ {H_1}\left( \mu \right) $、$ {H_2}\left( \mu \right) $是关于μ的二次函数. 与路面附着系数μ密切相关的车辆状态参数共有18个，分别是$ {\lambda _i} $、$ {\alpha _i} $、$ {v_x} $、$ {v_y} $、$ {a_x} $、$ {a_y} $、$ {\omega _i} $、$ \delta $、$ \dot \gamma $.

$ \mu $和上述18个车辆状态参数间存在非线性关系，其函数关系如下：

(21)$ \mu = f\left( {{\lambda _i},{\alpha _i},{v_x},{v_y},{a_x},{a_y},{\omega _i},\delta ,\dot \gamma } \right). $

根据式(21)确定神经网络模型的输入向量组为$ \left[ {{\lambda _i},{\alpha _i},{v_x},{v_y},{a_x},{a_y},{\omega _i},\delta ,\dot \gamma } \right] $，其中角速度$ {\omega _i} $可由轮速传感器测量得到；车辆的纵向速度$ {v_x} $、横向速度$ {v_y} $、纵向加速度$ {a_x} $、横向加速度$ {a_y} $以及车辆横摆角速度$ \dot \gamma $可以由三维运动姿态测量系统测量得到；车轮纵向滑移率$ {\lambda _i} $可由式(8)计算得到；车轮侧偏角$ {\alpha _i} $由式(9)计算得到；前轮转向角$ \delta $可由方向盘转角传感器测量并换算得到.

基于MATLAB/Simulink软件对汽车动态特性建模以获取数据集，车辆运动状态参数见表1^[11].

采集车辆制动过程中18个车辆状态参数，作为神经网络数据集. 设置路面附着系数取值范围为$ 0.1 \leqslant \mu \leqslant 1 $.0，步数为0.1，采样频率为200 Hz. 当路面附着系数≤0.5时，车辆制动的初始速度设为15 m/s；当路面附着系数＞0.5时，车辆制动的初始速度设为30 m/s. 车速一旦降为0，即停止数据采集. 综合统计后获得数据集样本空间大小为8.4万组（其中高附着路面有4.9万组，低附着路面有3.5万组）. 划分数据集为训练集、验证集与测试集3部分，以训练后续神经网络.

为了克服基于传统Elman神经网络的路面附着系数估计方法容易引起路面附着系数预估不准确，导致车辆滑移率控制效果不佳的问题，提出基于变权重PSO-Elman神经网络的路面附着系数估计系统.

Elman神经网络由输入层、隐含层、承接层和输出层组成，承接层与隐含层神经元个数相同，其拓扑结构如图2所示^[21].

Elman神经网表达式如下：

(22)$ \left. \begin{gathered} \mu \left( k \right) = g\left[ {{{\boldsymbol{w}}_3}{\boldsymbol{x}}\left( k \right)+{{{b}}_2}} \right], \\ {\boldsymbol{x}}\left( k \right) = f\left\{ {{{\boldsymbol{w}}_{\text{1}}}{{\boldsymbol{x}}_{\text{c}}}\left( k \right)+{{\boldsymbol{w}}_2} {{\boldsymbol{u}}\left( {k - 1} \right)} +{{\boldsymbol{b}}_1}} \right\}, \\ {{\boldsymbol{x}}_{\text{c}}}\left( k \right) = {\boldsymbol{x}}\left( {k - 1} \right). \\ \end{gathered} \right\} $

式中：$ \mu (k) $为$ k $时刻估计输出的路面附着系数预估值，$ {\boldsymbol{x}}(k) $与$ {{\boldsymbol{x}}_{\text{c}}}(k) $为隐含层和承接层在$ k $时刻的输出向量；$ {\boldsymbol{u}}(k - 1) $为$ k - 1 $时刻输入层的输入向量， $ {{\boldsymbol{w}}_3} $、$ {{\boldsymbol{w}}_2} $、$ {{\boldsymbol{w}}_1} $分别为隐含层到输出层、输入层到隐含层、承接层到隐含层权值矩阵，$ {{\boldsymbol{b}}_1} $为输入层阈值向量，$b_2 $为输出层阈值，$ f\left( \cdot \right) $、$ g\left( \cdot \right) $分别为隐含层神经元和输出神经元的传递函数.

Elman神经网络误差函数表达式如下：

(23)$ E\left( \mu \right) = {\left[ {{\mu _{\text{a}}} - \mu \left( k \right)} \right]^2}/2. $

式中：$ {\mu _{\mathrm{a}}} $为实际路面附着系数.

进行传统Elman神经网络模型训练，并确定初始权值矩阵为

(24)$ {\boldsymbol{p}} = \left[ {{{\boldsymbol{w}}_1},{{\boldsymbol{w}}_2},{{\boldsymbol{w}}_3}}^{\mathrm{T}} \right]. $

将初始权值矩阵$ {\boldsymbol{p}} $作为粒子群初始种群之一辅助粒子群初始化，在优化过程中为粒子提供更好的起点，增加搜索空间的多样性，加速收敛过程，提高优化算法的效率和准确性.

通过引入粒子群算法降低网络训练绝对误差和，实现网络权值优化，提高神经网络权值更新稳定性.

将基于Elman神经网络权值矩阵优化细分为4个步骤. 确定最大迭代次数和适应度函数；计算粒子群初始适应度函数值与初始权值矩阵；为了更有效平衡粒子全局、局部搜索能力，采用线性惯性权重递减策略对粒子群算法进行改进^[22]并输出改进粒子群权重；得到优化后的神经网络权值矩阵，通过训练输出路面附着系数估计值.

1）最大迭代次数和适应度函数确定.

过少的迭代次数可能导致网络未能充分学习，而过多的迭代次数会增加训练时间和计算成本. 为了兼顾考虑网络性能、训练时间和计算成本等因素，根据文献[22]，设定迭代次数$ {T_{{\text{max}}}} $=100，且在迭代过程中权值会不断逼近最优解. 在训练过程中，将迭代100次后得到的权值矩阵，作为神经网络的权值矩阵进行训练.

为了评价神经网络估计精度，引入神经网络绝对误差和：

(25)$ {E_{{\text{sumy}}}} = \sum {\left| {E\left( \mu \right)} \right|} . $

神经网络绝对误差和越小，网络精度越高.

在神经网络训练过程中，以$ {E_{{\text{sumy}}}} $作为适应度值来调整粒子群位置以调整相应神经网络参数，适应度值越低，表明输出值越准确.

2）初始最优适应度函数值与初始最优权值矩阵计算.

粒子群算法以群体智能为寻优原理，利用模型适应度值，使得输出值迅速向最优解靠近.

Elman神经网络误差函数对权值矩阵中的每一个值求偏导：

(26)$ \left.\begin{array}{l}\dfrac{\partial E}{\partial w_{3 j}}=-\delta_0 x_j(k) , \\\dfrac{\partial E}{\partial w_{2 j i}}=-\delta_0 \delta_1 w_{3 j} u_i(k-1) , \\\dfrac{\partial E}{\partial w_{1 j l}}=-\delta_0 \delta_2 w_{3 j} x_{{\mathrm{c}}, j}(k) , \\\delta_0=\mu_{\mathrm{a}}-\mu(k), \\\delta_1=\dfrac{\partial f}{\partial w_{2 j i}},\\\delta_2=\dfrac{\partial f}{\partial w_{1 j l}}.\end{array}\right\} $

式中：i=1, 2,···, q，q为输入层神经元的个数；j=1, 2, ···, n，l=1, 2, ···, n，n为隐含层/承接层神经元个数；w_3j表示w₃第j个值，$\partial {w_{2ji}} $表示w₂第j行i列的值，$\partial {w_{1ji}} $表示w₁第j行l列的值；${x_j}\left( k \right) $表示${{\boldsymbol{x}}}\left( k \right) $第j个值，${u_i}\left( {k - 1} \right) $表示${{\boldsymbol{u}}}\left( {k - 1} \right) $第i个值，${u_{{\mathrm{c}},j}}\left( k \right) $表示${{\boldsymbol{u}}_{{\mathrm{c}}}}\left( k \right) $第j个值.

采用梯度下降训练方法提高训练效率，各参数修正量与负梯度成正比：

(27)$ \Delta w=-\eta \frac{\partial E}{\partial w}. $

式中：$ \eta $为比例系数.

为了更直观地表示神经网络权值矩阵，将式(26)代入式(27)，得到

(28)$ \left.\begin{array}{l}\Delta w_{1 i j}=\eta \delta_0 x_j(k), \\\Delta w_{2 j i}=\eta \delta_0 \delta_1 w_{3 j} u_i(k-1) , \\\Delta w_{3 j l}=\eta \delta_0 \delta_2 w_{3 j} x_{{\mathrm{c}}, j}(k) .\end{array}\right\} $

式(28)即为Elman神经网络各层之间的权值矩阵，PSO算法中每一个粒子代表权值矩阵中的一个值.

粒子群优化算法的搜索空间维度为S维，粒子的种群规模为M，其组成的粒子群为

(29)$ {\boldsymbol{x}} = [{x_1},{x_2}, \cdots ,{x_M}]. $

第d个粒子的位置信息、速度信息为

(30)$ {{\boldsymbol{X}}_d} = [{x_{d1}},{x_{d2}}, \cdots ,{x_{dS}}], $

(31)$ {{\boldsymbol{V}}_d} = [{v_{d1}},{v_{d2}}, \cdots ,{v_{dS}}]. $

式中：$ d = 1,2, \cdots ,M $.

当前粒子搜索的个体最优位置信息为

(32)$ {{\boldsymbol{p}}_d} = [{p_{d1}},{p_{d2}}, \cdots ,{p_{dS}}]. $

种群中所有粒子的最优位置信息为

(33)$ {{\boldsymbol{p}}_{\text{g}}} = [{p_1},{p_2}, \cdots ,{p_M}]. $

初始化粒子群得到最优适应度函数值E_sumy-0与最优权值矩阵 $ {{\boldsymbol{p}}_{{\text{g-}}0}} = \left[ {{{\boldsymbol{w}}_{1 \text{-} 0}},{{\boldsymbol{w}}_{2 \text{-} 0}},{{\boldsymbol{w}}^{\mathrm{T}}_{3 \text{-} 0}}} \right] $.

3）变权重PSO权值矩阵优化.

通过变权重粒子群优化算法进行神经网络权值矩阵优化，流程图详见图3.

粒子以适应度值为依据调整速度与位置，自发靠近最优解的粒子位置. $ {E_{{\text{sumy}}}}\left( \cdot \right) $为最小优化问题的求解目标函数，第d个粒子当前的个体最优位置为

(34)$ {{\boldsymbol{p}}_d}\left( {T + 1} \right) = \left\{ \begin{gathered} {{\boldsymbol{p}}_d}\left( T \right),\;{E_{{\text{sumy}}}} \left[ {{{\boldsymbol{X}}_d}\left( {T + 1} \right)} \right] \geqslant {E_{{\text{sumy}}}} \left[ {{{\boldsymbol{p}}_d}\left( T \right)} \right] ; \\ {{\boldsymbol{X}}_d}\left( {T + 1} \right),\;{E_{{\text{sumy}}}} \left[ {{{\boldsymbol{X}}_d}\left( {T + 1} \right)} \right] \leqslant {E_{{\text{sumy}}}} \left[ {{{\boldsymbol{p}}_d} \left( T \right)} \right] . \\ \end{gathered} \right. $

最优解为粒子当前的全局最优位置：

(35)$ \begin{split} {E_{{\text{sumy}}}}&\left[ {{{\boldsymbol{p}}_{\text{g}}}\left( T \right)} \right] = \min \;\left\{ {E_{{\text{sumy}}}}\left[ {{{\boldsymbol{p}}_0}\left( T \right)} \right],\right. \\ &\left.{E_{{\text{sumy}}}}\left[ {{{\boldsymbol{p}}_1}\left( T \right)} \right], \cdots, {E_{{\text{sumy}}}}\left[ {{{\boldsymbol{p}}_M}\left( T \right)} \right]\right\}. \end{split} $

式中：$ {{\boldsymbol{p}}_{\text{g}}}\left( T \right) $为此轮迭代过程中的最优解；$ {E_{{\text{sumy}}}} \left[ {{{\boldsymbol{p}}_{\text{g}}}\left( T \right)} \right] $为此轮迭代过程中的最优适应度值；$ {{\boldsymbol{p}}_{{\text{g}} {\text{-}} T}} = {{\boldsymbol{p}}_{\text{g}}}\left( T \right) $，$ {E_{{\text{sumy-}}T}} = {E_{{\text{sumy}}}}\left[ {{{\boldsymbol{p}}_{\text{g}}}\left( T \right)} \right] $.

在粒子群算法的迭代寻优过程中，每个粒子分别得到在s维中的个体最优值$ {{{p}}_{ds}} $和全局最优值$ {{{p}}_{{\text{g}s}}} $，然后分别更新自身的速度和位置：

(36)$ \begin{split} {v_{ds}}\left( {T+1} \right) =& {w_{\text{g}}}{v_{ds}}\left( T \right)+{c_1}{r_1}\left[ {{{{p}}_{ds}} - {x_{ds}}\left( T \right)} \right]+ \\ &{c_2}{r_2}\left[ {{{{p}}_{{\text{g}}s}} - {x_{ds}}\left( T \right)} \right], \end{split}$

(37)$ {x_{ds}}\left( {T+1} \right) = {x_{ds}}\left( T \right)+{v_{ds}}\left( {T+1} \right). $

式中：$ s = 1,2,\cdots,S $；$ {w_{\text{g}}} $为权重；c₁、c₂为学习因子，值为2； r₁、r₂为 [0,1.0]的随机数；$ {v_{ds}}(T+1) $为下次迭代时第d个粒子在s维的速度，变化范围为[−0.5,0.5]，$ {x_{ds}}(T+1) $为该粒子对应的位置，变化范围为[−1.0,1.0].

由于标准粒子群算法中权重固定，粒子全局搜索能力和局部搜索能力比重不能得到有效的平衡，使用线性权重递减策略对标准粒子群权重$ {w_{\text{g}}} $进行改进，改进粒子权重为

(38)$ {w_{{\text{g}} {\text{-}} T}} = {w_{{\text{max}}}} - \frac{{\left( {{w_{{\text{max}}}} - {w_{{\text{min}}}}} \right) T}}{{{T_{{\text{max}}}}}}. $

式中：$ {w_{{\text{max}}}} $和$ {w_{{\text{min}}}} $分别为权重的最大值与最小值，分别取0.9、0.4.

随着每一次迭代过程，权值矩阵都会进一步被优化，直至达到最大迭代次数后得到优化后的权值矩阵.

4）神经网络模型训练.

若训练达到最大迭代次数100则训练结束，输出的最优解$ {{\boldsymbol{p}}_{{\text{g-10}}0}} = \left[ {{{\boldsymbol{w}}_{1 {\text{-}} 100}},{{\boldsymbol{w}}_{2 {\text{-}} 100}},{{\boldsymbol{w}}^{\mathrm{T}}_{3 {\text{-}} 100}}} \right] $作为Elman神经网络各层之间连接的权值矩阵. 此时神经网络模型进行最终训练并输出路面附着系数估计值.

准确识别路面附着系数是获得车辆理想防抱死控制效果的前提. 现有标准路面种类较少，无法实时匹配任一路面^[23]，因此，在滑移率控制过程中亟须拟合出路面附着系数特性曲线，用于实施匹配不同的路面工况. 理想滑移率估计单元如图4所示.

理想滑移率估计单元具体设计步骤如下.

1）基于路面附着系数特性曲线得到实际路面附着系数$ {\mu _{\text{a}}} $和与之相对应的理想滑移率$ {\lambda _{\text{p}}} $.

标准路面附着系数根据$ {\mu _\phi } {\text{-}} \lambda $曲线函数计算得到^[20]

(39)$ {\mu _\phi }\left( \lambda \right) = {c_1}(1 - {{\mathrm{e}}^{ - {c_2}\lambda }}) - {c_3}\lambda . $

式中：$ \lambda $为路面滑移率；c₁、c₂、c₃为路面状态参数；$ {\mu _\phi } $为可利用附着系数，即路面附着率，最大路面附着率即为路面附着系数.

根据式(39)确定实际路面附着系数和与之相对应的理想滑移率：

(40)$ \left. \begin{gathered} {\mu _{\text{a}}} = {c_1} - \dfrac{{{c_3}}}{{{c_2}}}\left(1+\lg \;\dfrac{{{c_1}{c_2}}}{{{c_3}}}\right), \\ {\lambda _{\text{p}}} = \dfrac{1}{{{c_2}}}\lg \;\dfrac{{{c_1}{c_2}}}{{{c_3}}}. \\ \end{gathered} \right\} $

常见的典型路面主要有覆冰、覆雪、潮湿鹅卵石、重潮湿沥青、典型潮湿沥青、轻潮湿沥青、干水泥路面和干沥青路面这8类.

基于式(40)，典型路面相关参数值见表2^[24]. 通过MATLAB软件绘制出典型路面附着系数特性曲线如图5所示.

2）$ {\mu _{\text{a}}} $-$ {\lambda _{\text{p}}} $曲线拟合得到理想滑移率.

将8类标准路面的$ {\lambda _{\text{p}}} $和$ {\mu _{\text{a}}} $通过函数拟合得到拟合曲线，分别使用幂函数拟合^[25]、高斯函数拟合^[26]、傅里叶函数拟合^[27]3种方式进行特性曲线拟合，如图6所示.

得到的拟合公式分别如下：

(41)$ {\hat \mu _{\text{p}}}\left( \lambda \right) = 24.19{\lambda ^{1.69}}, $

(42)$ \begin{split} {{\hat \mu }_{\text{f}}}\left( \lambda \right) =& - 0.618\;4\cos \; \left( {15.98\lambda } \right) - \\ &0.221\;8\sin \; \left( {15.98\lambda } \right)+0.700\;1, \end{split} $

(43)$ {\hat \mu _{\text{g}}}\left( \lambda \right) = 1.217{{\text{exp}}\;\left[{ - {{\left( {\dfrac{{\lambda - 1.6}}{{1.922}}} \right)}^2}}\right]}. $

式中：$ {\hat \mu _{\text{p}}} $、$ {\hat \mu _{\text{f}}} $、$ {\hat \mu _{\text{g}}} $分别为幂函数、傅里叶函数、高斯函数拟合曲线.

为了进一步量化评价拟合方法，采用相关系数（R²）评价上述3种方法的拟合精度：

(44)$ \left. \begin{gathered} {{R_1^2}} = 1 - {{\sum\limits_{i=1}^8 {{{\left( {{{\hat \mu }_{{\text{p}}i=1}} - {\mu _{{\text{a}}i}}} \right)}^2}} }}\bigg/{{\sum\limits_{i=1}^8 {{{\left( {{{\bar \mu }_{{\text{p}}i}} - {\mu _{{\text{a}}i}}} \right)}^2}} }}, \\ {{R_2^2}} = 1 - {{\sum\limits_{i=1}^8 {{{\left( {{{\hat \mu }_{{\text{f}}i}} - {\mu _{{\text{a}}i}}} \right)}^2}} }}\bigg/{{\sum\limits_{i=1}^8 {{{\left( {{{\bar \mu }_{{\text{f}}i}} - {\mu _{{\text{a}}i}}} \right)}^2}} }}, \\ {{R_3^2}} = 1 - {{\sum\limits_{i=1}^8 {{{\left( {{{\hat \mu }_{{\text{g}}i}} - {\mu _{{\text{a}}i}}} \right)}^2}} }}\bigg/{{\sum\limits_{i=1}^8 {{{\left( {{{\bar \mu }_{{\text{g}}i}} - {\mu _{{\text{a}}i}}} \right)}^2}} }}. \\ \end{gathered} \right\} $

式中：$ {\mu _{{\text{a}}i}} $为8类不同路面附着系数；$ {\overline{\mu }}_{{\text{p}}{i}}、{\overline{\mu }}_{{\text{f}}{i}}、{\overline{\mu }}_{{\text{g}}{i}} $为幂函数、高斯函数、傅里叶函数拟合平均值，$ {\widehat{\mu }}_{\text{p}i}、{\widehat{\mu }}_{\text{f}i}、{\widehat{\mu }}_{\text{g}i} $为与之对应的预测值，$ {R_1^2}、 {{R_2^2}}、{{R_3^2}} $为与之对应的相关系数.

使用MATLAB计算出拟合相关系数分别为0.9941、0.9989、0.9996. 傅里叶函数拟合相关系数更接近1.0，拟合效果最好. 因此选用傅里叶函数拟合曲线$ {\hat \mu _{\text{f}}}\left( \lambda \right) $作为理想滑移率-路面附着系数参考曲线，得到的理想滑移率为$ {\lambda _{\text{d}}}\left( {{\mu _{\text{f}}}} \right) $.

采用MATLAB/Simulink软件，对所提出的基于变权重PSO-Elman神经网络的路面附着系数估计方法的有效性进行仿真分析.

引入路面附着系数均方根误差$ {E_{{\text{rmse\_}}\mu }} $和滑移率均方误差$ {E_{{\text{mse\_}}\lambda }} $作为评价路面附着估计效果的直接评价指标和辅助评价指标，分别定量评估路面附着系数估计效果和滑移率控制效果：

(45)$ {E_{{\text{rmse\_}}\mu }} = \left[ {\frac{1}{N}\sum\limits_{k = 1}^N {{{\left[ {{\mu _{\text{a}}} - \mu \left( k \right)} \right]}^2}} }\right]^{1/2} , $

(46)$ {E_{{\text{mse\_}}\lambda }} = \frac{1}{M_{\mathrm{t}}}{\sum\limits_{k = 1}^{M_{\mathrm{t}}} {\left| {\lambda \left( k \right) - {\lambda _{\text{d}}}\left( k \right)} \right|} ^2}. $

式中：N为采样频率200 Hz情况下的测试样本数量，$ \lambda \left( k \right) $为k时刻车轮滑移率仿真值，$ {\lambda _{\text{d}}}\left( k \right) $为k时刻车轮理想滑移率，M_t为测试样本数量.

在仿真过程中须先对隐含层/承接层神经元个数$ n $进行设定. 采用“试算法”对$ n $逐步放大仿真，考虑神经网络输入量个数，根据文献[28]中所提出关于隐含层/承接层神经元个数的确定原则，将研究中所设计的隐含层/承接层神经元个数试算范围确定为$ {\text{4}} \leqslant n \leqslant {\text{13}} $，并依据路面附着系数均方根误差$ {E_{{\text{rmse\_}}\mu }} $大小来判断所设定的$ n $是否理想，$ {E_{{\text{rmse\_}}\mu }} $越小说明所设定的$ n $越理想. 仿真结果如表3所示. 可以看出，当n=8时，$ {E_{{\text{rmse\_}}\mu }} $最小，神经网络训练能够达到最好的效果，因此设本研究中隐含层/承接层神经元个数$ n $=8.

为了验证所提出的基于变权重PSO-Elman神经网络的路面附着系数估计方法的有效性，将参考文献[15]中Elman神经网络作为仿真对照组，更全面地评估所提出方法在路面附着系数估计方面的预估准确性和稳定性.

车轮滑移率除受到路面附着系数的影响，还受到车辆状态参数（比如轮胎状态、传感器测量精度）、防抱死控制算法等因素影响. 使用控制变量法，在不改变其他因素前提下，均采用文献[29]提出的区间二型模糊防抱死控制算法对车辆进行防抱死控制，只改变路面附着系数，以其对滑移率控制的影响评判变基于变权重PSO-Elman神经网络的路面附着系数估计方法的有效性.

模拟防抱死应用场景下较为常见、具有代表性且相对危险的4类工况，即覆冰、覆雪、潮湿鹅卵石和典型潮湿沥青路面工况，记为工况1~4，其中车辆在工况1~3中行驶的初始速度设为15 m/s，在工况4中行驶的初始速度为30 m/s.

代表性路面参数如表4所示^[20]. 如图7所示分别表示在工况1~4下路面附着系数预估值随时间变化曲线，路面附着系数均方根误差统计结果如表5所示. 其中，“实际值”、“Elman”、“变权重PSO-Elman”分别表示实际路面附着系数、使用基于传统Elman和变权重PSO-Elman神经网络的路面附着系数估计方法得到的路面附着系数预估值，p为基于变权重PSO-Elman神经网络的路面附着系数估计方法较基于传统Elman神经网络的路面附着系数估计方法预估均方根误差的降低比率.

如图8所示分别表示在工况1~4下车轮滑移率随时间的变化曲线，滑移率均方误差见表6. 其中，“理想值”、“Elman”、“变权重PSO-Elman”分别表示理想滑移率、使用基于传统Elman和变权重PSO-Elman神经网络的路面附着系数估计方法得到的滑移率仿真值，q为基于变权重PSO-Elman神经网络的路面附着系数估计方法较基于传统Elman神经网络的路面附着系数估计方法在滑移率均方误差上的降低比率.

由图7可知，基于变权重PSO-Elman神经网络的路面附着系数估计方法所估计的路面附着系数更加接近真实值，在4类工况下路面附着系数均方根误差较传统Elman神经网络分别降低37.89%、47.42%、30.70%和26.47%.

如图8所示，在4类工况中使用基于变权重PSO-Elman神经网络的路面附着系数估计方法的车辆在防抱死过程中滑移率更加接近理想滑移率，均方误差较传统Elman神经网络分别降低6.64%、7.32%、6.63%和20.74%. 仿真结果表明使用基于变权重PSO-Elman神经网络的路面附着系数估计方法的车辆在防抱死过程中路面附着系数识别更为准确，滑移率控制更为精准，波动程度更小，制动安全性更高.

选取覆雪路面对接潮湿鹅卵石路面作为工况5进行变附着路面仿真，路面附着系数在1 s时刻由0.196突变为0.383，车辆初始速度为15 m/s.

如图9、10所示分别表示变附着系数路面行驶工况下路面附着系数、车轮滑移率随时间变化的曲线.

如图9所示，在工况5下路面附着系数突变瞬间（由0.196上升到0.383），使用基于Elman神经网络、变权重PSO-Elman神经网络的路面附着系数估计方法得到的路面附着系数预估值随时间变化曲线的收敛时间几乎相同，而使用基于变权重PSO-Elman神经网络路面的附着系数估计方法估计的路面附着系数更加接近真实值，相比传统Elman神经网络估值方法，均方根误降低了19.2%，表明基于变权重PSO-Elman神经网络的路面附着系数估计方法估计准确度更高.

如图10所示，在工况5下使用基于变权重PSO-Elman神经网络的路面附着系数估计方法的车辆在防抱死过程中滑移率更加接近理想滑移率，均方误差较传统Elman神经网络前后轮分别降低7.39%、6.44%，说明使用基于变权重PSO-Elman神经网络的路面附着系数估计方法的车辆在防抱死过程中路面附着系数识别更为准确，滑移率能够以更小波动保持在理想滑移率附近，制动稳定性更高.

仿真模拟定附着与变附着2类行驶工况下，所提出的基于变权重PSO-Elman神经网络的路面附着系数估计方法的粒子群权重变化的有效性. 其中，定附着系数选用覆冰路面工况，变附着系数选取覆雪路面对接潮湿鹅卵石路面工况，初始速度皆设为15 m/s.

如图11所示为定附着工况、变附着工况下路面附着系数预估值随时间的变化曲线，路面附着系数均方根误差如表7所示. 其中，“PSO-Elman”表示基于PSO-Elman神经网络的路面附着系数估计方法（未变权重）得到的路面附着系数预估值，p₁为基于PSO-Elman神经网络的路面附着系数估计方法较基于传统Elman神经网络的路面附着系数预估方法在均方根误差上的降低比率，p₂为基于变权重PSO-Elman神经网络的路面附着系数估计方法较基于传统Elman神经网络的路面附着系数预估方法在均方根误差上的降低比率，p₃为基于变权重PSO-Elman神经网络的路面附着系数估计方法较基于PSO-Elman神经网络的路面附着系数预估方法在均方根误差上的降低比率. 可以看出，在定附、变附工况下使用基于变权重PSO-Elman神经网络的路面附着系数估计方法估计的路面附着系数皆比基于PSO-Elman神经网络的路面附着系数估计方法所得结果更加接近真实值. 由表7可知，变权重PSO-Elman神经网络较PSO-Elman神经网络的路面附着系数均方根误差分别降低30.59%、10.02%，说明本研究针对粒子群实施“变权重”，有助于提升PSO-Elman神经网络关于路面附着系数识别的估计精度.

（1）设计基于变权重PSO-Elman神经网络的路面附着系数估计方法，利用变权重粒子群算法降低网络训练绝对误差和实现神经网络权值矩阵优化，克服传统Elman神经网络在汽车制动过程可能陷入局部最优的难题，提升路面附着系数的估计精度和控制系统鲁棒性.

（2）在定附着路面制动工况下，采用基于变权重PSO-Elman神经网络的路面附着系数估计方法相较于传统Elman神经网络方法，车轮滑移率控制更加接近理想滑移率，结果表明变权重PSO-Elman神经网络对于路面附着系数具有更好的估计效果.

（3）在变附着系数制动工况下，采用基于变权重PSO-Elman神经网络的路面附着系数估计方法相较于传统Elman神经网络方法，车轮滑移率控制更加接近理想滑移率，结果表明基于变权重PSO-Elman神经网络的路面附着系数估计方法鲁棒性更强.

（4）下一步，将针对所提出估计方法筹划进行路面附着系数估计实验，构建变参数动态仿真及实验对比验证与校正平台，开展面向紧急制动/非紧急制动、驱动、转向等工况的路面附着系数估计的有效性验证，循环优化相关参数，为理论研究提供更加可靠的支持和验证.