汽车滑行试验道路阻力和风阻系数不确定度分析

摘要

对汽车滑行试验过程中滑行速度进行统计学分析，发现各个轮次的滑行车速偏差近似与平均滑行车速呈线性关系。基于滑行车速，分别利用v-F方法和v-t方法拟合得到汽车的滑行阻力和风阻系数，分析其对应优化问题的凸性和滑行阻力、风阻系数的不确定度；结果表明，v-F方法对应拟合问题是凸的，而v-t方法对应的拟合问题是非凸的。在拟合道路载荷时，需要选择泛化能力较强的迭代算法，同时需要根据参数的物理意义，设置合理的初始值和约束边界，以便得到合理的拟合参数值和道路阻力；而且相比v-F方法，v-t方法得到的风阻系数不确定度更小，置信区间也更窄。另外，基于车辆滑行速度的统计规律，进行蒙特卡洛模拟，分别应用v-F和v-t方法，得到道路阻力和风阻系数近似为正态分布；结果表明，v-F方法会放大随机因素的影响，v-t方法得到的道路阻力和风阻系数均更为稳定。

关键词

汽车滑行道路阻力; 风阻系数; 随机过程; 优化问题凸性; 不确定度; 蒙特卡洛模拟

汽车行驶过程中的行驶阻力主要包括空气阻力、滚动阻力、传动系阻力等。GB 18352.6—2016《轻型汽车污染物排放限值及测量方法（中国第六阶段）》^［1］及SAE J 2263—1996《Road Load Measurement Using Onboard Anemometry and Coastdown Techniques》^［2］均详细描述了滑行试验测试汽车阻力的方法，并对测试试验环境、车辆状态要求和数据处理方法有明确说明。按照标准规定：滑行试验时，首先将汽车设置为空挡，然后加速，直到比选定汽车速度高10 km/h，测量速度滑行过程中的速度和时间，多次重复试验，直到试验结果满足统计要求。

汽车空挡滑行过程时，其动力学方程可描述为：

$$F=m_{\mathrm{e}}a=D_{\mathrm{m}\mathrm{e}\mathrm{c}\mathrm{h}}+D_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}+D_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}}$$

（1）

式中： $m_{\mathrm{e}}$为包括旋转质量在内的试验车等效质量； a为汽车加速度； $D_{\mathrm{m}\mathrm{e}\mathrm{c}\mathrm{h}}$为机械阻力， $D_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}$为空气阻力，$D_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}}$为坡度阻力，可分别由下式表述：

$$D_{\mathrm{m}\mathrm{e}\mathrm{c}\mathrm{h}}=A_{\mathrm{m}}+B_{\mathrm{m}}v+C_{\mathrm{m}}{v}^{\mathrm{2}}$$

（2）

$$D_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}=\frac{\mathrm{1}}{\mathrm{2}}\rho A{v}_{\mathrm{r}}^{\mathrm{2}}{C}_{\mathrm{d}}$$

（3）

$$D_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{v}}=mg\mathrm{d}h/\mathrm{d}s$$

（4）

式中： $A_{\mathrm{m}}$、$B_{\mathrm{m}}$、$C_{\mathrm{m}}$为机械阻力系数； $v$为车速； $\rho$为空气密度； $C_{\mathrm{d}}$为风阻系数； $v_{\mathrm{r}}$为空气相对于汽车的速度； $m$为汽车质量； $g$为重力加速度； $\mathrm{d}h/\mathrm{d}s$为车辆前进方向试验道路坡度的正弦函数。

最终得到的动力学方程如下：

$$\begin{array}{l}F=m_{\mathrm{e}}a=A_{\mathrm{m}}+B_{\mathrm{m}}v+C_{\mathrm{m}}{v}^{\mathrm{2}}+\\ \frac{\mathrm{1}}{\mathrm{2}}\rho A{v}_r^{\mathrm{2}}{C}_{\mathrm{d}}+mg\mathrm{d}h/\mathrm{d}s\end{array}$$

（5）

上述标准对于阻力方程做了简化，使用速度的二次多项式描述道路阻力，即：

$$F=m_{\mathrm{e}}a=A+Bv+C{v}^{\mathrm{2}}$$

（6）

通过对试验测试的速度、时间数据可得到汽车的加速度，利用式（6）进行二次多项式曲线拟合即可得到待确定参数，从而得到各个车速下的道路阻力。显然，由式（5）到式（6）做了大量的简化。标准对于滑行试验场的环境风速、坡度、汽车磨合状态等均做了详实规定，并对最终的二次曲线进行了环境修正，从而忽略坡度，环境条件等对于道路阻力的影响。分析式（6）中各个参数的物理意义，A、B可认为为机械阻力系数，尽管滚动阻力并不能由车速的线性关系表示，但在车速相对较低时（$v\le \mathrm{120}\mathrm{k}\mathrm{m}/\mathrm{h}$），使用线性关系描述滚动阻力系数是较为合理的^［3］。因此，机械阻力中的二次项系数相比风阻系数的影响可忽略不计，C可认为仅仅与风阻系数有关。尽管系数C并不能完全准确描述汽车的风阻系数，但可认为C与风阻系数是强线性相关的，滑行试验也是一种间接评价汽车风阻系数的方式^［4］。

然而，滑行阻力的测试结果存在很大的不确定性，测试得到的滑行曲线和拟合参数可能变化很大^［5］，这主要来源于以下几个方面：

（1） 物理环境和汽车状态的变化，即测试过程中环境中的环境风速、坡度、路面凹凸情况、温湿度、汽车状态的变化等。尽管标准对于试验精度和试验环境做了要求，但并不能完全消除环境因素的影响。而且，从汽车使用场景的角度考虑，汽车真实的驾驶环境本身就具备很大的随机性，复杂环境下的道路阻力和风阻系数反而更具备真实性，汽车阻力、机械阻力系数、风阻系数等对环境因素的依赖关系更值得关注，以便评估使用场景下的真实能耗。Saurabhh等^［6］利用滑行试验研究了温湿度对于道路阻力的影响，并得到了各个阻力系数的标准差。Andrew^［7］研究了不同油温下传动阻力的变化。

环境因素对于阻力的影响最明显体现在风阻系数的变化。首先，环境中的自然风速带来了偏航。Barden^［8］通过路试得到了不同路面环境下的偏航角概率分布，结果表明，偏航角的概率分布与路面环境存在明显的依赖关系，平均的偏航角度可达7.9。J1252^［9］给出了平均风速下的风阻系数，即利用风洞试验结果和统计的平均风速进行积分，得到考虑环境风的平均风阻系数。Howell^［10］计算了EPA工况下，不同车型考虑环境风速的平均风阻系数，并进行了线性拟合。统计结果显示，相比风洞中无偏航下的风阻系数，EPA工况下的平均风阻系数要高约4.4%。此外，环境风中的湍流成分也对风阻系数存在影响。Howell^［11］描述了一种理论估计的方法，即利用准静态假设，估计风阻系数的非定常量。Howell等^［12］利用该方法，计算了多个车型在不同路面环境下的风阻系数非定常量，结果表明，湍流在高车速下对于风阻系数的影响较小，而在低车速下，增大风阻系数超过0.1。

（2） 滑行试验过程中的采集误差。滑行试验过程中主要通过高精度GPS得到汽车车速，但测试过程中，可能存在信号丢失。此外，采集设备存在系统误差、随机误差以及信号噪声干扰等。尽管可以通过传感器标定、重复试验、滤波等方式部分消除，但这种方式也具备一定的不确定性，这主要源于测试环境并不具备高度的重复性。

（3） 数据处理的过程也具备不确定性。按照标准，需要计算各个基准车速下的加速度。这种加速度的计算与选取的速度间隔有关，而且这种计算方法也是一种加速度近似，具备一定误差。龚春忠等^［13］研究了速度间隔的选取大小对于拟合阻力曲线的影响。理想测试精度下，可直接根据前后时间步的车速得到加速度。但是，车速采集的误差使计算的加速度具备很大的震荡性。此外，非线性曲线拟合本身也具备不确定性。曲线拟合本质上是优化问题，得到的极值点依赖于优化问题的凸性，约束条件，初始值，选择的优化算法等。

综上，上述描述的多方面因素的不确定性带来了测试得到的阻力和风阻系数的不确定性，合理的评估方式是利用统计方法，分析环境因素的概率分布，得到对应道路阻力和风阻系数的统计规律等。本文基于某电动SUV滑行试验数据，首先分析滑行车速的统计学规律。而后，针对动力学方程，分析不同拟合方法对应优化问题的凸性、优化算法和不确定度等。最后，基于车速的统计规律，进行蒙特卡洛模拟，研究滑行道路阻力和风阻系数的不确定度和分布规律。

1 试验方法

本文试验车辆为岚图某款电动后驱SUV和四驱SUV，后驱配置车辆为样车阶段试制车，四驱配置车辆为小批量量产阶段车辆，滑行试验分别在襄樊达安汽车检测中心和中汽中心盐城汽车试验场进行。滑行试验按照标准进行，试验分三段进行，后驱的试验数据更为完整，本文选择后驱的滑行试验数据进行详细分析。此外，后驱车辆在中国汽车工程研究院风洞中心进行风洞试验，进行滑行试验和风洞试验的风阻系数对比。

2 试验数据处理与分析方法

滑行数据的处理涉及数字滤波、统计学分析、非线性拟合方法及对应优化算法、拟合参数不确定度分析方法等，下面进行描述。

2.1　滑行速度-时间数据分析

对于采集到的速度时间数据，本文使用Kalman滤波方法过滤采样误差，再剔除滑行过程中的无效数据，最终得到的各个轮次高速段、中速段、低速段的速度时间数据如下图1所示。

  

  

  

图1  滑行速度-时间

Fig.1  Coasting velocity versus time

分析各个速度区间的均值及标准差。定义速度均值和标准差为：

$$v_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}={\displaystyle \sum _{j=\mathrm{1}}^m}{v}_{ij}/m$$

（7）

$$v_{i-\mathrm{s}\mathrm{t}\mathrm{d}}=\sqrt[]{\frac{\sum (v_{ij}-v_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}{)}^{\mathrm{2}}}{m-\mathrm{1}}}$$

（8）

式中： $v_{ij}$为第j次滑行第i个速度采样点的滑行速度，$v_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$为第i个速度采样点的平均滑行速度； m为总滑行次数； n为每次滑行的速度采样点数。如图2所示，可以发现，在开始滑行时，$v_{i-\mathrm{s}\mathrm{t}\mathrm{d}}$的速度标准差较小，而随着速度均值的减小，其速度标准差近似线性增大。

图2  滑行车速平均值与标准差

Fig.2  Mean and standard error of coasting velocity

进一步，可定义速度偏差量为：

$$v_{i-\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}=v_{ij}-v_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$$

（9）

如果把汽车滑行过程作为一个随机过程来分析，则滑行车速可表述为：$v\left(t\right)=\mu \left(t\right)+\epsilon (t)$，其中$v\left(\mathrm{t}\right)$是在时间t上的随机变量，$\mu \left(\mathrm{t}\right)$是在时间t上的期望，即为式（9）中定义的$v_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$，$\epsilon (t)$是在时间t上的误差项或随机波动，即为$v_{i-\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}$。如图3为低速段滑行，速度偏差量随速度均值减小而近似线性地增大，只是线性关系的权重系数不同。这说明每次采集的速度偏差并不是独立的，而是相关的。对（$v_{i-\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}$，$v_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$）进行线性回归分析，即：

$$v_{i-\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}=k_j{v}_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}+b_j,j=\mathrm{1,2},\cdots ,m$$

（10）

图3  滑行速度偏差量

Fig.3  Coasting velocity deviation

定义回归分析的残差为：

$$s_j=\sqrt[]{\frac{\sum (v_{i-\mathrm{e}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}}-(k_j{v}_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}+b_j{))}^{\mathrm{2}}}{n-\mathrm{3}}}$$

（11）

式中： n—3代表回归分析的自由度。

如果假定，各个轮次的汽车滑行过程中，汽车状态和滑行时的环境因素是相对稳定的。（$k_j,b_j$）与各个轮次的滑行是一一对应的，可以将其理解为滑行时汽车和环境中相对稳定的因素的影响，表征每一次滑行过程中的平均风速、方向、温湿度、汽车状态等，而线性回归的残差可理解为滑行过程中随机因素的影响，表征湍流风速、风速方向瞬时变化、测试随机误差等。因此，滑行速度可表示为：

$$v_{ij}=k_j{v}_{i-\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}+b_j+{\epsilon }_{ij}$$

（12）

式中：${\epsilon }_{ij}~N(\mathrm{0},s_j^{\mathrm{2}})$，服从正态分布。而（$k_j,b_j$）的分布表征的是汽车状态、环境因素带来的阻力、速度等较为稳定的变化。显然，这取决于滑行过程的环境平稳性与汽车状态平稳性。按照上述标准规定，一般滑行过程的轮次至少为6次，本文中滑行次数超过10次，但样本量并不足以准确分析（$k_j,b_j$）分布律，这里假设（$k_j,b_j$）和$s_j$服从Weibull分布，其分布概率密度函数如下：

$$f\left(x\right)=\left\{\begin{array}{c}\mathrm{0}\\ \frac{k}{\lambda }{\left(\frac{x}{\lambda }\right)}^{k-\mathrm{1}}{e}^{-{\left(\frac{x}{\lambda }\right)}^k}\end{array}\right.$$

（13）

式中： $k$为形状参数， $\lambda$为比例参数。如图4所示，为$k_j$和$s_j$的Weibull分布拟合结果。如上所述，由于滑行的样本量不足，并不能得到较好的拟合结果，但仍可看到使用Weibull分布是较为合理的。这里，可以通过改变Weibull分布的比例参数，调整其方差，即调整滑行过程的环境平稳性与汽车状态平稳性。

图4  $k_j$和$s_j$ Weibull分布拟合

Fig.4  Weibull distribution fitting of $k_j$ and $s_j$

2.2　汽车阻力拟合方法

如前所述，汽车滑行时的动力学方程如（6）所述，通过对阻力和滑行速度进行二次多项式拟合，即可得到阻力系数，这种方法定义为v-F方法。显然，加速度的计算与速度间隔的选取有关，即：

$$a=\frac{\left(v+∆v\right)-\left(v-∆v\right)}{∆t}=\frac{\mathrm{2}∆v}{∆t}$$

（14）

若考虑速度采样过程中的随机误差，即：

$$\begin{array}{l}a=\frac{\left(v+∆v+{\epsilon }_{\mathrm{1}}\right)-\left(v-∆v+{\epsilon }_{\mathrm{2}}\right)}{∆t}=\\ \frac{\mathrm{2}∆v}{∆t}+({\epsilon }_{\mathrm{1}}-{\epsilon }_{\mathrm{2}})/∆t\end{array}$$

（15）

式中： $∆v$为速度间隔； $∆t$为时间差； ${\epsilon }_{\mathrm{1}}$和${\epsilon }_{\mathrm{2}}$为速度采样的随机误差，尽管通过多次重复试验、滤波等可以部分消除。但是，$a$的准确性与$∆v$和$({\epsilon }_{\mathrm{1}}-{\epsilon }_{\mathrm{2}})$的相对大小有关。理想状况下，加速度应该通过前后时间步的速度得到，即$∆v$足够小，但因为采样误差等因素的存在，此方法得到的加速度震荡性很大，如图5所示。

图5  车辆滑行速度与加速度

Fig.5  Vehicle coasting velocity and acceleration

而如果增大$∆v$，计算的加速度则可能存在偏移，因为这种方式只是一种加速度的近似计算，可能放大误差，最终拟合的阻力和阻力系数也与$∆v$的大小有关。此外，式（6）是存在理论解的，即：

$$v=-\frac{B+\mathrm{t}\mathrm{a}\mathrm{n}\left(\left(\frac{t}{m_{\mathrm{e}}}\sqrt[]{\mathrm{4}AC-B^{\mathrm{2}}}-\mathrm{2}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{B+\mathrm{2}C{v}_{\mathrm{0}}}{\sqrt[]{\mathrm{4}AC-B^{\mathrm{2}}}}\right)\right)/\mathrm{2}\right)\sqrt[]{\mathrm{4}AC-B^{\mathrm{2}}}}{\mathrm{2}C}$$

（16）

式中： $v_{\mathrm{0}}$为初始速度。直接通过速度时间数据进行拟合，这种方法称为v-t方法。

2.3　非线性曲线拟合对应优化问题分析

曲线拟合问题可转化为优化问题，即在约束条件下，找到参数估计值，使对应的损失函数最小。理论上，非线性优化问题是否能找到最优解很大程度上取决于其是否凸性，对于凸优化问题，其局部最优解即为全局最优解，而对于非凸优化问题，并不一定能找到其全局最优解，其解很大程度取决于初始值，优化算法等，由此可能带来拟合参数和道路阻力的不确定性。判断一个优化问题的凸性有多种方法，本文计算损失函数的黑塞矩阵，并判断其在约束集内是否正定来评估非线性优化问题的凸性。

对于v-F方法，其对应的损失函数即为：

$$L=\sum (F_k-(A+B{v}_k+C{v_k}^{\mathrm{2}}{))}^{\mathrm{2}},k=\mathrm{1,2}\dots N$$

（17）

二次多项式拟合对应的问题，即在采样数据（$F_k,v_k$）下，找到（A，B，C）使$L$达到最小值。首先需要判断$L$是否为凸函数。其黑塞矩阵如下：

$$H=\left(\begin{array}{ccc}\mathrm{2}N& {\displaystyle \sum _{k=\mathrm{1}}^N}\mathrm{2}{v}_k& {\displaystyle \sum _{k=\mathrm{1}}^N}\mathrm{2}{v_k}^{\mathrm{2}}\\ {\displaystyle \sum _{k=\mathrm{1}}^N}\mathrm{2}{v}_k& {\displaystyle \sum _{k=\mathrm{1}}^N}\mathrm{2}{v_k}^{\mathrm{2}}& {\displaystyle \sum _{k=\mathrm{1}}^N}\mathrm{2}{v_k}^{\mathrm{3}}\\ {\displaystyle \sum _{k=\mathrm{1}}^N}\mathrm{2}{v_k}^{\mathrm{2}}& {\displaystyle \sum _{k=\mathrm{1}}^N}\mathrm{2}{v_k}^{\mathrm{3}}& {\displaystyle \sum _{k=\mathrm{1}}^N}\mathrm{2}{v_k}^{\mathrm{4}}\end{array}\right)$$

（18）

与参数$A,B,C$无关，计算黑塞矩阵的特征值。对于本文中的滑行试验数据，其特征值为［1.4×10¹¹， 10.7， 8.6×10⁵］，均大于0，由此可以判断其正定性，即损失函数为凸函数。若给定$A,B,C\ge \mathrm{0}$约束条件，该约束条件是凸集，对应的优化问题是凸优化问题，理论上，其全局最优解即为局部最优解，即v-F方法容易找到全局最优解的^［14］。

而对于v-t方法，其拟合过程虽然更加直接，但其对应的损失函数复杂，即：

$$L=\sum {\left(v_k-\frac{B+\mathrm{t}\mathrm{a}\mathrm{n}\left(\left(\frac{t}{m_{\mathrm{e}}}\sqrt[]{\mathrm{4}AC-B^{\mathrm{2}}}-\mathrm{2}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{B+\mathrm{2}C{v}_{\mathrm{0}}}{\sqrt[]{\mathrm{4}AC-B^{\mathrm{2}}}}\right)\right)/\mathrm{2}\right)\sqrt[]{\mathrm{4}AC-B^{\mathrm{2}}}}{\mathrm{2}C}\right)}^{\mathrm{2}},k=\mathrm{1,2}\dots N$$

（19）

其黑塞矩阵除了与采样点（$v_k,t_k$）有关，也与$A,B,C$的取值范围有关，数学分析判断比较困难。本文利用数值方法计算黑塞矩阵，对本文中的滑行试验速度时间数据，并根据$A,B,C$对应物理意义的约束范围判断损失函数的凸性，在约束条件：

$$\mathrm{100}\le A\le \mathrm{200},\mathrm{0}\le B\le \mathrm{1},\mathrm{0.03}\le C\le \mathrm{0.04}$$

（20）

对$A,B,C$在约束区间内进行均匀采样，并计算其黑塞矩阵的特征值。结果表明，在上述约束区间，损失函数的黑塞矩阵的特征值并不是恒大于0的，因此，该优化问题是非凸的，在拟合过程中，是可能陷入局部解的。

2.4　非线性曲线拟合的不确定度评价方法

一般来说，判断非线性拟合参数的不确定度是困难的^［15］。一种近似的方法是^［16］，基于线性近似，利用误差传递理论初步分析参数的不确定度，即对于函数：

$$X=f\left(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_l\right)$$

（21）

因为自变量$(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_l)$的误差带来$X$的误差可表示为：

$$V\left(X\right)={\left(\frac{\partial X}{\partial x_{\mathrm{1}}}\right)}^{\mathrm{2}}V\left(x_{\mathrm{1}}\right)+{\left(\frac{\partial X}{\partial x_{\mathrm{2}}}\right)}^{\mathrm{2}}V\left(x_{\mathrm{2}}\right)+\dots +\mathrm{2}\frac{\partial X}{\partial x_{\mathrm{1}}}\frac{\partial X}{\partial x_{\mathrm{2}}}C\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)+\dots$$

（22）

式中： $V\left(x_{\mathrm{1}}\right),V\left(x_{\mathrm{2}}\right),\dots$为自变量的方差，$C\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)$为$x_{\mathrm{1}},x_{\mathrm{2}}$的协方差。这种方式得到的不确定度要低于实际非线性拟合参数的不确定度，但也是一种较为合理的评价方式。此外，可以通过损失函数的协方差矩阵或者黑塞矩阵的逆矩阵进行分析，本文采用黑塞矩阵的逆矩阵进行参数的不确定度分析^［15］。同时，某些情况下，损失函数的黑塞矩阵或者协方差矩阵方法可能失效，同时使用F检验法，对参数的置信区间进行分析。

3 试验结果讨论

基于以上描述的方法，对本文中的滑行试验数据进行分析，详细结果如下。

3.1　道路阻力和风阻系数拟合算法

首先，为了验证上述分析的准确性，利用不同的优化算法进行v-t方法和v-F方法的阻力曲线拟合。本文利用Python SciPy优化包，应用Levenberg-Marquardt和BFGS两种常用局部优化算法，以及Differential Evolution和Basinhopping全局优化算法。Levenberg-Marquardt算法是一种用于非线性最小二乘问题的优化算法，结合了梯度下降和高斯-牛顿法的思想，常用于数据拟合和参数估计。BFGS是一种基于拟牛顿法，通过逐步改进目标函数的二阶近似矩阵来迭代求解的局部优化算法。Differential Evolution是一种全局优化算法，常用于求解连续型的非线性优化问题。该算法基于种群的演化过程，在搜索空间中寻找最优解。Basinhopping（盆地跳跃）也是一种全局优化算法，常用于解决具有多个局部极小点的非线性优化问题。该算法结合了局部搜索和随机搜索的策略，通过跳出当前局部极小点所在的“盆地”，探索搜索空间中的其他区域，以寻找全局最优解。此外，还考虑Brute这一种暴力算法，是一种简单直接的解决问题的方法，它通过尝试所有可能的解决方案来求解问题。

无论选择何种迭代算法，都需要根据参数的物理意义设置合理的约束边界和初始值，以便加快收敛速度，同时避免陷入局部解。理论上，v-F方法对应的优化问题是凸的，得到的解不依赖于初始值的选择，但是，如果选择的优化算法泛化能力不足，某些初始值可能会导致算法陷入局部解。如图6为在A₀=0.01，B₀=0，C₀=0.01时，v-F方法得到的道路阻力。可以看到，对于某些初始值，迭代算法是容易陷入局部解的。此时，拟合的残差可能很大。如图7为得到的拟合参数C，跳跃性很大，此时，拟合的道路阻力和风阻系数可能都与真实情况不符。而对于v-t方法，由于其对应优化问题是非凸的，除了初始值外，还需要对参数的约束边界进行合理设置，否则迭代算法甚至无法收敛。而且，由于v-t方法对应的损失函数复杂得多，其迭代收敛速度相比v-F方法较慢。

图6  不同优化算法的道路阻力拟合

Fig.6  Road load fitting by different optimization algorithms

图7  不同优化算法拟合参数值

Fig.7  Fitting parameter values of different optimization algorithms

综上，在处理数据时，优先选择泛化能力较高的优化算法，同时需要根据参数的物理意义设置合理的约束边界和拟合初始值，以便高效准确的找到合理解。实际处理滑行数据时，可以将v-t和v-F方法结合起来，将v-F方法得到的解作为v-t方法的初始值，得到所需的全局解。

3.2　速度间隔对于拟合阻力和阻力系数的影响分析

如上文所示，速度间隔$∆v$可能会对拟合的滑行阻力和拟合参数值存在影响。这里分别选取$∆v=\mathrm{0.2,0.5,1.0,2.0,3.0,4.0,5.0}$，使用v-F方法进行滑行阻力拟合。图8为得到的拟合参数值以及车速为120 km/h时的阻力。同时根据标准的一般过程，只取多个基准速度的阻力和车速进行拟合，得到拟合参数值和车速为120 km/h时的阻力。可以看到，速度间隔的大小对于得到的道路阻力影响不大，但是，拟合参数值对于速度间隔较为敏感，下文中，速度间隔$∆v=\mathrm{2}\mathrm{k}\mathrm{m}/\mathrm{h}$。

图8  不同速度间隔得到的道路阻力和拟合参数值

Fig.8  Road load and fitting parameter values for different velocity spacing

3.3　v-F方法和v-t方法的不确定度分析

这里统一选择Levenberg-Marquardt优化算法，同时设置参数初始值为：A₀=120，B₀=0.5，C₀=0.03，参数的约束范围分别为：

$$\mathrm{100}\le A\le \mathrm{200},\mathrm{0}\le B\le \mathrm{1},\mathrm{0.03}\le C\le \mathrm{0.04}$$

（23）

使用v-F方法和v-t方法进行拟合参数的不确定度分析，对于v-F方法，速度间隔$∆v=\mathrm{2}\mathrm{k}\mathrm{m}/\mathrm{h}$。如表1所示，为两种方法得到的拟合参数值以及标准差。可以发现，两种方法得到的拟合参数值比较接近，这也说明使用v-t方法进行滑行数据分析是可行的。而且，v-t方法得到的参数不确定度明显更小，如上文分析，v-F方法涉及加速度的计算，可能放大测试过程中随机因素的影响，因此，拟合过程中损失函数的残差也更大，相应的，其随机误差带来的拟合参数不确定度也越大。

表1  拟合参数值与标准差

Tab.1  Fitting parameters and standard deviation

方法	A	σA	B	σB	C	σC
v-F	157.093	4.323	0.797	0.107	0.036 2	0.000 6
v-t	156.377	1.430	0.806	0.035	0.036 2	0.000 2

使用F检验法，进行拟合参数的置信区间分析，如图9所示。可以看出，v-t方法的置信区间更小，即v-t方法的稳定性更好。v-F和v-t方法，95%的置信水平，C的置信区间分别为［0.035 0，0.037 4］和［0.035 8，0.036 6］。

（a）  v-F方法

(b)  v-t方法

图9  拟合参数置信区间

Fig.9  Confidence intervals of fitting parameters

如果对动力学方程不做简化，即：

$$\begin{array}{l}F=m_{\mathrm{e}}a=A_{\mathrm{m}}+B_{\mathrm{m}}v+C_{\mathrm{m}}{v}^{\mathrm{2}}+\\ \frac{\mathrm{1}}{\mathrm{2}}\rho A{v}_{\mathrm{r}}^{\mathrm{2}}{C}_{\mathrm{d}}+mg\mathrm{d}h/\mathrm{d}s\end{array}$$

（24）

$$v_{\mathrm{r}}^{\mathrm{2}}=v^{\mathrm{2}}+v_{\mathrm{w}}^{\mathrm{2}}+\mathrm{2}v{v}_{\mathrm{w}}\mathrm{c}\mathrm{o}\mathrm{s}\theta$$

（25）

$$C_{\mathrm{d}}=a_{\mathrm{0}}+a_{\mathrm{1}}Y+a_{\mathrm{2}}{Y}^{\mathrm{2}}+a_{\mathrm{3}}{Y}^{\mathrm{3}}+a_{\mathrm{4}}{Y}^{\mathrm{4}}$$

（26）

$$Y=\mathrm{t}\mathrm{a}{\mathrm{n}}^{-\mathrm{1}}\left(\frac{v_{\mathrm{w}}\mathrm{s}\mathrm{i}\mathrm{n}\theta }{v_{\mathrm{w}}+v\mathrm{c}\mathrm{o}\mathrm{s}\theta }\right)$$

（27）

式中： $v_{\mathrm{w}}$为环境风速； $\theta$为环境风速角；$Y$为偏航角度。这其中涉及的未知拟合参数有8个。进行拟合时，并根据风洞试验结果和转鼓试验，设置各个拟合参数的初始值。但是，损失函数具有强烈的非线性，参数之间也存在高度相关性，图10为得到的拟合参数相关系数矩阵。可以发现，（A_m，B_m）、（A_m，a₁）、（A_m，a₂）、（A_m，a₄）、（B_m，a₃）、（B_m，a₁）、（C_m，a₀）、（a₁，a₃）、（a₂，a₄）的相关系数都超过了0.95。显然，从数值分析的角度，拟合参数是存在严重冗余的^［17］，得到的拟合参数值不确定度很大，很容易导致过拟合，得到不符合实际的道路阻力和风阻系数。对于非线性曲线拟合可能的参数冗余，可以使用正则化技术来限制参数的自由度，提高模型的泛化能力。同时，需要更为准确的真实道路环境下的风阻系数描述，上述标准提出使用偏航角的四次多项式描述滑行过程中的风阻系数，但实际处理过程中还是使用二次多项式拟合滑行阻力。而且这种风阻系数的描述方法显然也不够合理，因为其假定每个滑行车速下的风阻系数与偏航角度是一一对应的，这显然忽略了流场的滞后效应，这里仅初步探讨，不进行深入研究。

图10  拟合参数相关系数矩阵

Fig.10  Correlation coefficient matrix of fitting parameters

3.4　蒙特卡洛方法分析滑行阻力和风阻系数的分布律

为了进一步分析滑行试验得到的道路阻力的分布规律，根据以上关于滑行速度的统计假设，通过蒙特卡洛模拟，研究v-t方法和v-F方法得到的道路载荷和风阻系数的统计规律。通过上述v-t方法拟合得到的滑行阻力系数，构建理想中的速度时间曲线，而后根据上述关于滑行速度的统计假设，加入速度偏差量。这里，模拟了真实汽车滑行试验过程，每一个蒙特卡洛采样对应多轮次的滑行试验，正反方向次数相同，总次数为Coasting-nums，而后进行滑行数据平均，使用上述方法得到对应的道路阻力和风阻系数。蒙特卡洛采样次数为3000。

首先设置Coasting-nums=6，标准中规定的最小滑行次数。如图所示，11为前300次样本得到的道路阻力和风阻系数。可以看到，即便进行了数据平均，阻力的随机性是很明显的，偏差可超过60 N。这里定义C_{d WT}为风洞试验得到的无偏航角度下的风阻系数，C_d-C_{d WT}为风阻系数偏差量。可以发现，两种方法得到的风阻系数期望值很接近。与风洞试验结果对比，期望值与5°偏航下的风阻系数更为接近，这也与参考文献中［4］的结果吻合。但是，v-t方法得到的风阻系数偏差量随机性更小，使用正态分布，对数据分布进行拟合。如图12所示，其中，v-t方法得到的风阻系数偏差量近似呈N（0.017 1， 0.005 2²），v-F方法得到的风阻系数偏差量呈N（0.017 6，0.007 7²）。此外，两种方法得到在车速为130 km/h时的阻力偏差较小，其中两种方法得到的阻力分别近似呈N（872.8， 14.3²）和N（872.7，14.4²）。

图11  蒙特卡洛试验得到的道路阻力和风阻系数

Fig.11  Road load and air drag coefficient from Monte Carlo simulations

  

  

图12  道路载荷概率分布

Fig.12  Probability distribution of road load

需要说明的是，上述蒙特卡洛试验，都根据标准的规范过程，每一次模拟都对应多个轮次的滑行试验，而后进行数据平均，得到对应的道路阻力和风阻系数，这相当于对环境因素的影响进行了一定程度的过滤，否则，得到的道路阻力和风阻系数的不确定度更大。Coasting-nums=1时，即不再进行数据平均，直接进行拟合，这是真实道路环境汽车阻力的近似模拟，如图13所示。此时，阻力和风阻系数的随机性更大，存在“双峰”，这是不同滑行方向的体现，在实际行驶过程中，存在顺风和逆风，上坡与下坡，空气阻力和坡度阻力的变化带来整车阻力的明显变化。

图13  道路载荷和风阻系数概率分布（Coasting-nums=1）

Fig.13  Probability distribution of road load and air drag coefficient (Coasting-nums=1)

此外，研究滑行次数对于得到的阻力和风阻系数的影响，如图14所示。统计上，提高滑行次数可降低道路阻力和风阻系数的不确定度。在Coasting-nums≥10时，滑行次数对于阻力和风阻系数的标准差影响很小，但整体上，v-t方法得到的风阻系数更符合稳定性要求。当Coasting-nums=20时，v-t方法得到的风阻系数标准差为0.002 6，具备较好的重复性。同时，对比四驱车辆和两驱车辆在不同试验场的阻力和风阻系数标准差，可以发现，四驱车辆得到的阻力和风阻系数不确定度更大，这也与真实情况相符，因为四驱车辆试验场的风场环境是要更不稳定。

图14  不同滑行轮次下的道路载荷期望和标准差

Fig.14  Expectation and standard deviation of road load and air drag coefficient for different coasting rounds

最后，改变Weibull分布的比例参数，研究v-t方法和v-F方法得到的风阻系数标准差，此时滑行次数均为6次。如上文所述，$s_j$表征的是测试过程中的随机性影响，增大$s_j$概率密度函数的比例参数，即增大滑行过程中的随机性。而$k_j$表征的是测试过程中的整体稳定性，增大$k_j$概率密度函数的比例参数，即增大滑行过程中各个轮次滑行试验之间总体汽车状态和环境状态的差异性。图15为不同比例参数下的道路阻力和风阻系数标准差。总体上，v-t方法得到的风阻系数不确定度更小，因为v-F方法计算加速度的过程会放大随机误差的影响。在$s_j$对应的比例参$\lambda =\mathrm{4}{\lambda }_{\mathrm{0}}$数，$k_j$对应的比例参数$\lambda =\mathrm{0.25}{\lambda }_{\mathrm{0}}$时，两种方法得到的风阻系数标准差相差约3倍。而$k_j$比例参数的变化，风阻系数的标准差区别较小。因此，v-t方法比v-F方法更具备稳定性，可用此方法，得到更为稳定的汽车滑行阻力和风阻系数。

图15  不同比例参数下的风阻系数标准差

Fig.15  Standard deviation of air drag coefficient for different scale parameter

（a） v-t方法 (b) v-F方法

4 结论

本文基于某电动SUV的滑行试验数据，对滑行数据得到道路载荷和风阻系数的不确定度进行了分析，可得到以下结论：

（1） 首先对得到的车辆滑行速度进行统计学分析，发现滑行车速的期望值与标准差近似呈线性关系。

（2） v-F方法对应拟合问题是凸的，而v-t方法对应的拟合问题是非凸的。优先选择泛化能力较强的迭代算法，同时需要根据参数的物理意义，设置合理的初始值和约束边界，以便避免陷入局部解，加快收敛速度，得到合理的拟合参数值和道路阻力。

（3） 使用v-F方法时，得到的拟合参数值与选取的速度间隔大小有关，但速度间隔对得到的道路阻力影响较弱； v-t方法得到的拟合参数更为稳定，其参数置信区间也更窄。

（4） 基于车辆滑行速度的统计规律，进行Monte Carlo模拟，分别应用v-t和v-F方法，得到道路阻力和风阻系数近似为正态分布，v-F方法会放大随机因素的影响，v-t方法得到的道路阻力和风阻系数均更为稳定。