近年来，我国汽车数量的持续增加引发了不少安全问题. 据统计^[1-2]可知，车车事故是交通事故中最常见且伤亡人数较高的情形. 研究汽车主动安全系统的先进辅助驾驶技术来提高汽车安全性已成为亟需解决的问题^[3].

目前，针对汽车主动安全系统的先进辅助驾驶技术大多数聚焦于紧急避撞. Cui等^[4]提出两阶段紧急转向避撞控制算法，该算法引导车辆通过环形路径进入相邻车道，使用转向控制来引导车辆沿中心线行驶. Zhu等^[5]基于五次多项式和模型预测控制算法(model predictive control,MPC)建立控制框架，综合考虑避碰效果、路径跟踪精度和车身稳定性要求. Pan等^[6]通过学习驾驶员数据，在紧急避撞过程中考虑驾驶员的反应特性，确保安全和稳定.

在不可避撞场景下，车辆无法完全避免碰撞. Vangi等^[7]研究车辆速度变化和碰撞位置对乘员损伤的影响，提出自适应调整纵向制动和转向的方法减小碰撞损伤，研究相关的碰撞模型^[8-9]. Parseh等^[10]结合车辆运动规划、车辆动力学和事故重建模型，讨论如何将车辆碰撞后运动纳入碰撞损伤最小化的研究中. 为了减小车辆在不可避撞场景下的碰撞损伤，Wang等^[11]提出自动驾驶车辆在不可避撞场景下的路径规划算法，以减小车辆碰撞损伤. 该算法采用模型预测控制，在代价函数中考虑碰撞损伤,碰撞损伤取决于碰撞角度、车辆速度和碰撞车辆的质量比. Parseh等^[12]提出以控制车辆极限的转向角率和制动率为基础创建离线轨迹库，利用基于事故数据的碰撞损伤评估模型，从轨迹库中选择最优轨迹并进行跟踪. Li等^[13]采用神经网络建立非线性碰撞损伤评估模型，在轨迹库对每一条轨迹进行潜在碰撞损伤的估计，得到碰撞损伤最小的轨迹.

以上研究对不可避撞场景车辆碰撞损伤最小化作出了积极探索，但仍具有一定的局限性. 1)仅考虑碰撞时或碰撞后的单一情况，未利用预测的碰后状态来优化碰撞前的车辆运动，以减小车辆碰撞损伤. 2)车辆碰撞模型由于碰撞位置的不同，建模公式有所差异，采用矢量方式将各种碰撞情形统一为一种建模方式. 3)轨迹库的生成仅基于前轮极限转向角和制动率，不能充分发挥轮胎纵向的力学性能，通过引入附加横摆力矩以提供车辆更大的潜力和灵活性. 4)在不可避撞场景下车辆易失稳，缺乏轨迹跟踪和横摆稳定性的协同控制.

参考Zhou等^[14]对车辆碰撞模型的研究，本文采用矢量方式将各种碰撞情形统一为一种建模方式，建立四自由度碰撞模型，提出分层结构的碰撞损伤最小化策略. 在决策层考虑车辆动力学约束，通过求解最优控制问题生成离线轨迹库，利用四自由度碰撞模型和碰撞前的车辆状态来估计碰后自车失稳风险，在线选择得到最优轨迹. 在运动控制层，基于模型预测控制建立轨迹跟踪与横摆稳定性协同控制器，保证轨迹的跟踪精度和车辆的稳定性.

考虑模型精度和计算量，在Simulink中建立包括车身和轮胎在内的非线性车辆模型. 如图1所示，忽略车辆悬架系统垂直运动、车身俯仰运动及侧倾运动，建立包括车辆纵向、侧向、横摆和4个车轮旋转的七自由度非线性车辆模型. 整车的动力学方程如下.

(1)$ \begin{split} m({{\dot v}_x} - {v_y}{\omega _z}) = \;&({F_{x{\mathrm{fl}}}}+{F_{x{\mathrm{fr}}}})\cos\; \delta +{F_{x{\mathrm{rl}}}}+{F_{x{\mathrm{rr}}}}- \\& ({F_{y{\mathrm{fl}}}}+{F_{y{\mathrm{fr}}}})\sin\; \delta . \end{split}$

(2)$ \begin{split} m({{\dot v}_y}+{v_x}{\omega _z}) =& ({F_{y{\mathrm{fl}}}}+{F_{y{\mathrm{fr}}}})\cos\; \delta +{F_{y{\mathrm{rl}}}}+{F_{y{\mathrm{rr}}}}+ \\& ({F_{x{\mathrm{fl}}}}+{F_{x{\mathrm{fr}}}})\sin\; \delta . \end{split}$

(3)$\begin{split} {I_{zz}}{{\dot \omega }_z} =\;& {l_{\rm{f}}}\{ ({F_{y{\mathrm{fl}}}} + {F_{y{\mathrm{fr}}}})\cos\; \delta + ({F_{x{\mathrm{fl}}}} + {F_{x{\mathrm{fr}}}})\sin\; \delta \} - {l_{\rm{r}}}({F_{y{\mathrm{rl}}}} + \\ &{F_{y{\mathrm{rr}}}}){\text+}0.5b({F_{x{\mathrm{rr}}}} - {F_{x{\mathrm{rl}}}}){\text+}\;0.5b\{ ({F_{x{\mathrm{fr}}}}\cos\; \delta - {F_{y{\mathrm{fr}}}}\times \\ &\sin\; \delta ) - ({F_{x{\mathrm{fl}}}}\cos\; \delta - {F_{y{\mathrm{fl}}}}\sin\; \delta )\} .\\[-5pt]\end{split} $

4个车轮的旋转动力学方程可以表示为

(4)$ {J_{{\mathrm{w}}ij}}{\dot \omega _{ij}} = {T_{{\mathrm{d}}ij}} - {T_{{\mathrm{b}}ij}} - {F_{xij}}R. $

式中：$ m $为整车质量，$ {v}_{x}、{v}_{y} $分别为车辆在车身坐标系下沿$ x $轴和$ y $轴方向的速度分量，$ {\omega }_{z} $为横摆角速度，$ \delta $为转向轮转角，$ {l}_{{\mathrm{f}}}、{l}_{{\mathrm{r}}} $分别为整车质心到前、后轴的距离，$ {I}_{zz} $为车身绕$ z $轴的转动惯量，$ b $为轮距，$ R $为车轮滚动半径，$ {F}_{xij} $和$ {F}_{yij}\mathrm{分}\mathrm{别} $为4个轮胎的纵向力和横向力，$ {T}_{\mathrm{d}ij} $和$ {T}_{{\mathrm{b}}ij} $分别为4个轮胎所受的驱动力和制动力，$ {J}_{{\mathrm{w}}ij} $和$ {\omega }_{ij}\mathrm{分}\mathrm{别} $为4个车轮的转动惯量和转动角速度，$ ij\in \{{\mathrm{fl}},{\mathrm{fr}},{\mathrm{rl}},{\mathrm{rr}}\} $分别指左前轮、右前轮、左后轮和右后轮.

车辆在高速行驶时，轮胎容易进入非线性区域，因此本研究的轮胎模型选用魔术公式. 在单一工况下，轮胎纵向力$ {F}_{x0} $和侧向力$ {F}_{y0} $的一般表达形式为

(5)$ \left. \begin{gathered} {F_{x0}} = {D_x}\sin\;\{ {C_x}\arctan\; [{B_x}\lambda - {E_x}({B_x}\lambda - \arctan \;({B_x}\lambda ))]\} , \\ {F_{y0}} = - \mu {D_y}\sin\; \{ (1.25 - 0.25\mu ){C_y}\arctan \;[(2 - \mu ){B_y}\alpha - \\ \;\;\;\;\;\;\;\;(2 - \mu ){E_y}((2 - \mu ){B_y}\alpha - \arctan\; ((2 - \mu ){B_y}\alpha ))]\} . \\ \end{gathered} \right\} $

式中：$ \mu $为路面附着系数，$ \lambda $为纵向滑移率，$ \alpha $为轮胎侧偏角，$ {B}_{x}、{B}_{y} $均为刚度因子，$ {C}_{x}、{C}_{y} $ 均为形状因子，$ {D}_{x}、{D}_{y} $均为峰值因子，$ {E}_{x}、{E}_{y} $均为曲率因子. 4个轮胎的$ \mathrm{\alpha }\mathrm{、}\mathrm{\lambda } $的计算式分别为

(6)$ \left. \begin{gathered} {\alpha _{{\mathrm{fl}}}},{\alpha _{{\mathrm{fr}}}} = \frac{{{v_y}+{l_{\mathrm{f}}}{\omega _z}}}{{{v_x} \mp 0.5b{\omega _z}}} - \delta , \\ {\alpha _{{\mathrm{rl}}}},{\alpha _{{\mathrm{rr}}}} = \frac{{{v_y} - {l_{\mathrm{r}}}{\omega _z}}}{{{v_x} \mp 0.5b{\omega _z}}}. \\ \end{gathered} \right\} $

(7)$ \left. \begin{gathered} {\lambda _{{\mathrm{fl}}}} = \frac{{R{\omega _{{\mathrm{fl}}}} - ({v_x} - 0.5b{\omega _z})}}{{\max \;\{{v_x} - 0.5b{\omega _z},R{\omega _{{\mathrm{fl}}}}\}}}, \\ {\lambda _{{\mathrm{fr}}}} = \frac{{R{\omega _{{\mathrm{fr}}}} - ({v_x}+0.5b{\omega _z})}}{{\max \;\{{v_x}+0.5b{\omega _z},R{\omega _{{\mathrm{fr}}}}\}}}, \\ {\lambda _{{\mathrm{rl}}}} = \frac{{R{\omega _{{\mathrm{rl}}}} - ({v_x} - 0.5b{\omega _z})}}{{\max \;\{{v_x} - 0.5b{\omega _z},R{\omega _{{\mathrm{rl}}}}\}}}, \\ {\lambda _{{\mathrm{rr}}}} = \frac{{R{\omega _{{\mathrm{rr}}}} - ({v_x}+0.5b{\omega _z})}}{{\max\; \{{v_x}+0.5b{\omega _z},R{\omega _{{\mathrm{rr}}}}\}}}. \\ \end{gathered} \right\} $

刚度因子、形状因子、峰值因子、曲率因子的表达式分别为

(8)$ \left.\begin{array}{l}C_x=b_0, \;B_x C_x D_x=\left(b_3 F_z^2 + b_4 F_z^2\right) / \exp\;({b_5 F_z}) , \\E_x=b_6 F_z^2+b_7 F_z+b_8 , \\D_y=a_1 F_z^2+a_2 F_z, \;C_y=a_0 , \\B_y C_y D_y=a_3 \sin\;\left(a_4 \arctan \;\left(a_5 F_z\right)\right) , \\E_y=a_6 F_z^2+a_7 F_z+a_8 .\end{array}\right\} $

式中：$ {a}_{0} \sim {a}_{8}、{b}_{0} \sim {b}_{8} $均为常数，数值如表1所示. 4个轮胎的垂直载荷计算式分别为

(9)$ {F_{z{\mathrm{fl}}}},{F_{z{\mathrm{fr}}}} = \frac{{mg}}{2}\left(\frac{{{l_{\mathrm{r}}}}}{L} - \frac{{{a_x}{h_{\mathrm{g}}}}}{{gL}} \mp \frac{{2{a_y}{h_{\mathrm{g}}}{l_{\mathrm{r}}}}}{{gbL}}\right), $

(10)$ {F_{z{\mathrm{rl}}}},{F_{z{\mathrm{fl}}}} = \frac{{mg}}{2}\left(\frac{{{l_{\mathrm{f}}}}}{L}+\frac{{{a_x}{h_{\mathrm{g}}}}}{{gL}} \mp \frac{{2{a_y}{h_{\mathrm{g}}}{l_{\mathrm{r}}}}}{{gbL}}\right). $

式中：$ {h}_{{\mathrm{g}}} $和$ L $分别为质心离地高度和轴距，$ {a}_{x} $和$ {a}_{y} $分别为车辆纵向加速度和侧向加速度.

当轮胎处于转向加制动联合工况时，轮胎的纵向力和侧向力应限制在附着椭圆上，考虑滑移率和侧偏角的轮胎模型表达式为

(11)$ \left. \begin{gathered} {\sigma _x} = \frac{\lambda }{{1+\lambda }},\;{\sigma _y} = \frac{{\tan \alpha }}{{1+\lambda }},\;\sigma {\text{ = }}\sqrt {\sigma _x^2+\sigma _y^2} . \\ {F_x} = \frac{{{\sigma _x}}}{\sigma }{F_{x0}}\left( {{\sigma _x}} \right),\;{F_y} = \frac{{{\sigma _y}}}{\sigma }{F_{y0}}\left( {{\sigma _y}} \right). \\ \end{gathered} \right\} $

为了预测两车碰撞后的自身车辆状态，提高碰撞模型精度，引入侧倾这一自由度，建立碰撞时的四自由度车辆动力学模型. 由于本文的研究场景主要针对交叉路口，两车的碰撞情形更复杂，通过参考Zhou等^[14]对追尾碰撞工况下车辆碰撞的研究，采用矢量方式将各种碰撞情形统一为一种建模方式. 如图2所示，假设碰撞力仅是水平的，碰撞位置由坐标$ ({d}_{x},{d}_{y},{d}_{z}) $确定. 车辆纵向、侧向、横摆和侧倾的动力学方程可以写为

(12)$ m({\dot v_x} - {v_y}{\omega _z}) = {F_x}, $

(13)$ m({\dot v_y}+{v_x}{\omega _z}) - {m_{\mathrm{R}}}h{\dot \omega _x} = {F_y}+{F_{y{\mathrm{f}}}}+{F_{y{\mathrm{r}}}}, $

(14)$ {I_{zz}}{\dot \omega _z}+{I_{xz}}{\dot \omega _x} = {d_x}{F_y} - {d_y}{F_x}+{l_{\mathrm{f}}}{F_{y{\mathrm{f}}}} - {l_{\mathrm{r}}}{F_{y{\mathrm{r}}}}, $

(15)$ \begin{split} &{I_{xx{\mathrm{s}}}}{{\dot \omega }_x}+{I_{xz}}{{\dot \omega }_z} - {m_{\mathrm{R}}}h({{\dot v}_y}+{v_x}{\omega _z}) = {F_y}({d_z} - h)+ \\ &({m_{\mathrm{R}}}gh - {K_{\mathrm{s}}})\psi - {D_{\mathrm{s}}}{\omega _x}. \end{split}$

式中：

(16)$ \left. \begin{gathered} {F_{y{\mathrm{f}}}} = {C_{\mathrm{f}}}\left(\delta - \frac{{{v_y}+{l_{\mathrm{f}}}{\omega _z}}}{{{v_x}}}\right),\left| {{F_{y{\mathrm{f}}}}} \right| < \mu mg\frac{{{l_{\mathrm{r}}}}}{L}; \\ {F_{y{\mathrm{r}}}} = {C_{\mathrm{r}}}\left(\frac{{ - {v_y}+{l_{\mathrm{r}}}{\omega _z}}}{{{v_x}}}\right),\left| {{F_{y{\mathrm{r}}}}} \right| < \mu mg\frac{{{l_{\mathrm{f}}}}}{L}. \\ \end{gathered} \right\} $

其中$ {m}_{{\mathrm{R}}} $为簧上质量，$ {m}_{{\mathrm{NR}}} $为簧下质量，$ {\omega }_{x} $为侧倾角速度，$ \varPsi $为侧倾角，$ {F}_{y\mathrm{f}}、{F}_{y{\mathrm{r}}} $分别为前、后轴的总侧向力，$ {F}_{x} $为沿车辆纵向方向的碰撞力分量，$ {F}_{y} $为沿车辆侧向方向的碰撞力分量，$ {C}_{\mathrm{f}}、{C}_{{\mathrm{r}}} $分别为前、后轴的总轮胎侧偏刚度，$ {K}_{{\mathrm{s}}} $为总悬架侧倾角刚度，$ {D}_{{\mathrm{s}}} $为总悬架侧倾阻尼，$ h $为簧上质量与簧下质量质心间的垂向距离，$ {d}_{x}、{d}_{y} $和$ {d}_{z} $分别为碰撞点在车身坐标系下的坐标，$ {I}_{xz} $为簧上质量绕车身坐标系$ x $轴和$ z $轴的惯性积，$ {I}_{xx{\mathrm{s}}} $为车辆绕车身坐标系$ x $轴的转动惯量.

根据碰撞场景定义了另外2个局部坐标系，如图3所示，坐标系$ x-y $固定在被撞车辆上，坐标系$ {x}'-{y}' $固定在施撞车辆上，分别对应各自的纵轴和横轴. 碰撞的典型持续时间为0.1~0.2 s，这种量级的持续时间可以采用交叉项积分的梯形对式(12)~(15)进行整理，可得如下公式.

(17)$ m({V_{1x}} - {v_{1x}}) - \frac{{\Delta t}}{2}m({V_{1y}}{\varOmega _{1z}}+{v_{1y}}{\omega _{1z}}) = {P_x}, $

(18)$\begin{split} & m({V_{1y}} - {v_{1y}})+\frac{{\Delta t}}{2}m({V_{1x}}{\varOmega _{1z}}+{v_{1x}}{\omega _{1z}}) - {m_{\rm{R}}}h({\varOmega _{1x}} - \\& {\omega _{1x}}) = {P_y} - \frac{{\Delta t}}{2}{C_{\rm{f}}}\left(\frac{{{V_{1y}} + {l_{\rm{f}}}{\varOmega _{1z}}}}{{{V_{1x}}}}+\frac{{{v_{1y}} + {l_{\rm{f}}}{\omega _{1z}}}}{{{v_{1x}}}}\right) - \frac{{\Delta t}}{2} \times \\ &{C_{\rm{r}}}\left(\frac{{{V_{1y}} - {l_{\rm{r}}}{\varOmega _{1z}}}}{{{V_{1x}}}}+\frac{{{v_{1y}} - {l_{\rm{r}}}{\omega _{1z}}}}{{{v_{1x}}}}\right),\end{split}$

(19)$ \begin{split} &{I_{zz}}({\varOmega _{1z}} - {\omega _{1z}})+{I_{xz}}({\varOmega _{1x}} - {\omega _{1x}}){\text{ = }}{d_x}{P_y} - {d_y}{P_x} - \frac{{\Delta t}}{2} \times \\ &{C_{\rm{f}}}{l_{\rm{f}}}\left(\frac{{{V_{1y}}+{l_{\rm{f}}}{\varOmega _{1z}}}}{{{V_{1x}}}}+\frac{{{v_{1y}}+{l_{\rm{f}}}{\omega _{1z}}}}{{{v_{1x}}}}\right)+\frac{{\Delta t}}{2}{C_{\rm{r}}}{l_{\rm{r}}}\left(\frac{{{V_{1y}} - {l_{\rm{r}}}{\varOmega _{1z}}}}{{{V_{1x}}}}+ \right. \\&\left.\frac{{{v_{1y}} - {l_{\rm{r}}}{\omega _{1z}}}}{{{v_{1x}}}}\right), \\[-1pt]\end{split} $

(20)$ \begin{gathered} {I_{xx{\mathrm{s}}}}({\varOmega _{1x}} - {\omega _{1x}})+{I_{xz}}({\varOmega _{1z}} - {\omega _{1z}}) - {m_{\rm{R}}}h({V_{1y}} - {v_{1y}}) - \frac{{\Delta t}}{2} \times \\ {m_{\rm{R}}}h({V_{1x}}{\varOmega _{1z}}+{v_{1x}}{\omega _{1z}}) = {P_y}({d_z} - h) - \frac{{\Delta t}}{2}{D_{\mathrm{s}}}({\varOmega _{1x}}+{\omega _{1x}}). \\ \end{gathered} $

式中：$ \Delta t $为两车接触时间，$ {P}_{x}、{P}_{y} $分别为被撞车辆所受冲量在车身坐标系上的分量，$ {V}_{ix}、{V}_{iy} $分别为碰撞后车辆在车身坐标系下沿$ x $轴和$ y $轴方向的速度分量，$ {\varOmega }_{iz} $为碰撞后的横摆角速度，$ {\varOmega }_{ix} $为碰撞后的侧倾角速度，$ i=\mathrm{1,2} $分别表示被撞车辆和施撞车辆. 利用同样的方法，可得施撞车辆相应的4个方程，这里省略.

仅由式(17)~(20)不足以求解得到碰撞瞬间车辆所受的冲量，需要引入恢复系数$ e $和切向相互作用系数$ {\varepsilon }_{\mathrm{p}} $.

$ e $被定义为碰撞前、后两车在接触面法向方向上的速度差比值的相反数，如下所示：

(21)$ e = - \frac{{{V_{{\mathrm{n}}2}} - {V_{{\mathrm{n}}1}}}}{{{v_{{\mathrm{n}}2}} - {v_{{\mathrm{n}}1}}}}. $

式中：$ {V}_{{\mathrm{n}}1}、{V}_{{\mathrm{n}}2} $分别为碰撞后被撞车辆和施撞车辆在接触面法向方向上的速度分量，$ {v}_{{\mathrm{n}}1}、{v}_{{\mathrm{n}}2} $分别为碰撞前被撞车辆和施撞车辆在接触面法向方向上的速度分量.

(22)$\begin{split} {V_{{\rm{n}}1}} = &({V_{1x}}\cos\; {\theta _1} - {V_{1y}}\sin\; {\theta _1}+{d_y}{\varOmega _{1z}})\cos\; \varGamma + \\ &({V_{1x}}\sin\; {\theta _1}+{V_{1y}}\cos\; {\theta _1} - {d_x}{\varOmega _{1z}})\sin\; \varGamma .\end{split} $

(23)$ \begin{split} {V_{{\rm{n}}2}} = ({V_{2x}}\cos\; {\theta _2} - {V_{2y}}\sin\; {\theta _2}+{d_{y'}}{\varOmega _{2z}})\cos\; \varGamma+ \\ ({V_{2x}}\sin\; {\theta _2}+{V_{2y}}\cos\; {\theta _2} - {d_{x'}}{\varOmega _{2z}})\sin\; \varGamma .\end{split} $

(24)$ \begin{split} {v_{{\mathrm{n}}1}} =& ({v_{1x}}\cos\; {\theta _1} - {v_{1y}}\sin\; {\theta _1}+{d_y}{\omega _{1z}})\cos\; \varGamma+ \\ &({v_{1x}}\sin\; {\theta _1}+{v_{1y}}\cos\; {\theta _1} - {d_x}{\omega _{1z}})\sin\; \varGamma . \end{split}$

(25)$ \begin{split} {v_{{\mathrm{n}}2}} = &({v_{2x}}\cos\; {\theta _2} - {v_{2y}}\sin\; {\theta _2}+{d_{y'}}{\omega _{2z}})\cos\; \varGamma + \\ &({v_{2x}}\sin\; {\theta _2}+{v_{2y}}\cos\; {\theta _2} - {d_{x'}}{\omega _{2z}})\sin\; \varGamma . \end{split} $

(26)$ \left. \begin{gathered} {d_x} = {d_1}\cos\; {\xi _1},\;{d_y} = {d_1}\sin\; {\xi _1}, \\ {d_{x'}} = {d_2}\cos\; {\xi _2},\;{d_{y'}} = {d_2}\sin\; {\xi _2}. \\ \end{gathered} \right\} $

式中：$ {\theta }_{1}、{\theta }_{2} $分别为两车的航向角，$ {\xi }_{1}、{\xi }_{2} $分别为车身坐标系中$ x $轴逆时针旋转到碰撞位置的夹角，$ \varGamma $为接触面坐标系与大地坐标系的夹角，$ {d}_{1}、{d}_{2} $分别为两车质心到接触点的距离.

切向相互作用系数$ {\varepsilon }_{{\mathrm{p}}} $被定义为接触面上切向冲量与法向冲量之比，

(27)$ {\varepsilon _{\mathrm{p}}} = {{{P_{\mathrm{t}}}}}/{{{P_{\mathrm{n}}}}}. $

式中：

(28)$ \left. \begin{gathered} {P_{\mathrm{t}}} = {P_y}\cos\; \varGamma - {P_x}\sin\; \varGamma , \\ {P_{\mathrm{n}}} = {P_y}\sin\; \varGamma +{P_x}\cos\; \varGamma . \\ \end{gathered} \right\} $

碰撞冲量从车身坐标系投影到大地坐标系，如下所示：

(29)$ \left. \begin{gathered} {P_x}\cos\; {\theta _1} - {P_y}\sin\; {\theta _1} = - [{P_{x'}}\cos\; {\theta _2} - {P_{y'}}\sin\; {\theta _2}], \\ {P_x}\sin\; {\theta _1}+{P_y}\cos\; {\theta _1} = - [{P_{x'}}\sin\; {\theta _2}+{P_{y'}}\cos\; {\theta _2}]. \\ \end{gathered} \right\} $

根据式(17)~(29)，可得

(30)$ {\boldsymbol{A}}{\boldsymbol{X}}=\left[\begin{array}{*{20}{l}}{{\boldsymbol{A}}}_{11}& {{{0}}}&{{\boldsymbol{A}}}_{13}\\{{{0}}}&{{\boldsymbol{A}}}_{22}&{{\boldsymbol{A}}}_{23}\\ {{\boldsymbol{A}}}_{31}&{{\boldsymbol{A}}}_{32}&{{\boldsymbol{A}}}_{33}\end{array}\right] {\boldsymbol{X}}={\boldsymbol{B}}. $

式(30)将碰撞后的车辆状态与碰撞前的车辆状态联系起来，通过碰撞前的车辆状态$ ({v}_{1x},{v}_{1y}, {\omega }_{1z},{\omega }_{1x}, {v}_{2x},{v}_{2y},{\omega }_{2z},{\omega }_{2x}) $，可以求解得到碰撞后的状态$ ({V}_{1x}, {V}_{1y},{\varOmega }_{1z},{\varOmega }_{1x},{V}_{2x},{V}_{2y},{\varOmega }_{2z},{\varOmega }_{2x}) $. 式中：

$\begin{split} {{\boldsymbol{A}}_{11}} = &\left[ {\begin{array}{*{20}{c}} m&0&{ - \dfrac{{\Delta t}}{2}m{V_{1y}}}&0 \\ 0&{m+\dfrac{{\Delta t}}{2}\dfrac{{{C_{\rm{f}}}+{C_{\rm{r}}}}}{{{V_{1x}}}}}&{\dfrac{{\Delta t}}{2}m{V_{1x}}+\dfrac{{\Delta t}}{2}\dfrac{{{C_{\rm{f}}}{l_{\rm{f}}} - {C_{\rm{r}}}{l_{\rm{r}}}}}{{{V_{1x}}}}}&{ - {m_{\rm{R}}}h} \\ 0&{\dfrac{{\Delta t}}{2}\dfrac{{{C_{\rm{f}}}{l_{\rm{f}}} - {C_{\rm{r}}}{l_{\rm{r}}}}}{{{V_{1x}}}}}&{\dfrac{{\Delta t}}{2}\dfrac{{{C_{\rm{f}}}{l_{\rm{f}}}^2+{C_{\rm{r}}}{l_{\rm{r}}}^2}}{{{V_{1x}}}}+{I_{zz}}}&{ - {I_{xz}}} \\ 0&{ - {m_{\rm{R}}}h}&{{I_{xz}} - \dfrac{{\Delta t}}{2}{m_{\rm{R}}}h{V_{1x}}}&{{I_{xx{\mathrm{s}}}}+\;\dfrac{{\Delta t}}{2}{D_{\mathrm{s}}}} \end{array}} \right],\\{{\boldsymbol{A}}_{22}} =& \left[ {\begin{array}{*{20}{c}} m&0&{ - \dfrac{{\Delta t}}{2}m{V_{2y}}}&0 \\ 0&{m+\dfrac{{\Delta t}}{2}\dfrac{{{C_{\rm{f}}}+{C_{\rm{r}}}}}{{{V_{2x}}}}}&{\dfrac{{\Delta t}}{2}m{V_{2x}}+\dfrac{{\Delta t}}{2}\dfrac{{{C_{\rm{f}}}{l_{\rm{f}}} - {C_{\rm{r}}}{l_{\rm{r}}}}}{{{V_{2x}}}}}&{ - {m_{\rm{R}}}h} \\ 0&{\dfrac{{\Delta t}}{2}\dfrac{{{C_{\rm{f}}}{l_{\rm{f}}} - {C_{\rm{r}}}{l_{\rm{r}}}}}{{{V_{2x}}}}}&{\dfrac{{\Delta t}}{2}\dfrac{{{C_{\rm{f}}}{l_{\rm{f}}}^2+{C_{\rm{r}}}{l_{\rm{r}}}^2}}{{{V_{2x}}}}+{I_{zz}}}&{ - {I_{xz}}} \\ 0&{ - {m_{\rm{R}}}h}&{{I_{xz}} - \dfrac{{\Delta t}}{2}{m_{\rm{R}}}h{V_{2x}}}&{{I_{xx{\mathrm{s}}}}+\;\dfrac{{\Delta t}}{2}{D_{\mathrm{s}}}} \end{array}} \right],\end{split} $

$ {{\boldsymbol{A}}_{23}} = \left[ {\begin{array}{*{20}{c}} { - 1}&0&0&0 \\ 0&{ - 1}&0&0 \\ {{d_{y'}}}&{ - {d_{x'}}}&0&0 \\ { - ({d_{z'}} - h)}&0&0&0 \end{array}} \right]{{\boldsymbol{A}}_{13}} = \left[ {\begin{array}{*{20}{c}} { - 1}&0&0&0 \\ 0&{ - 1}&0&0 \\ {{d_y}}&{ - {d_x}}&0&0 \\ { - ({d_z} - h)}&0&0&0 \end{array}} \right], $

$\begin{split} {{\boldsymbol{A}}_{31}} = &\left[ {\begin{array}{*{20}{c}} 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \\ { - \cos\; {\theta _1}\cos\; \varGamma - \sin\; {\theta _1}\sin\; \varGamma }&{\sin\; {\theta _1}\cos\; \varGamma - \cos\; {\theta _1}\sin\; \varGamma }&{ - {d_y}\cos\; \varGamma +{d_x}\sin\; \varGamma }&0 \end{array}} \right], \\{{\boldsymbol{A}}_{32}} = &\left[ {\begin{array}{*{20}{c}} 0&0&0&0 \\ 0&0&0&0 \\ 0&0&0&0 \\ {\cos\; {\theta _2}\cos\; \varGamma +\sin\; {\theta _2}\sin\; \varGamma }&{ - \sin\; {\theta _2}\cos\; \varGamma +\cos\; {\theta _2}\sin\; \varGamma }&{{d_{y'}}\cos\; \varGamma - {d_{x'}}\sin\; \varGamma }&0 \end{array}} \right] ,\end{split}$

$ {{\boldsymbol{A}}_{33}} = \left[ {\begin{array}{*{20}{c}} { - \sin\; \varGamma - {\varepsilon _{\mathrm{p}}}\cos\; \varGamma }&{\cos\; \varGamma - \;{\varepsilon _{\mathrm{p}}}\sin\; \varGamma }&0&0 \\ {\cos\; {\theta _1}}&{ - \sin\; {\theta _1}}&{\cos\; {\theta _2}}&{ - \sin\; {\theta _2}} \\ {\sin\; {\theta _1}}&{\cos\; {\theta _1}}&{\sin\; {\theta _2}}&{\cos\; {\theta _2}} \\ 0&0&0&0 \end{array}} \right], $

$ \boldsymbol{X}=\left[V_{1 x} ,\quad V_{1 y} ,\quad \varOmega_{1 z} ,\quad \varOmega_{1 x} ,\quad V_{2 x},\quad V_{2 y} ,\quad \varOmega_{2 z} ,\quad \varOmega_{2 x} ,\quad P_x ,\quad P_y ,\quad P_{x^{\prime}} ,\quad P_{y^{\prime}}\right]^{\mathrm{T}}. $

$ {\boldsymbol{B}} = {\left[ {\mathop {mv}\nolimits_{1x} ,\;0,\;0,\;0,\;\mathop {mv}\nolimits_{2x} ,\;0,\;0,\;0,\;0,\;0,\;0,\;e\left( {\mathop v\nolimits_{1x} \left( {{\mathrm{sin}}\;\theta_1 \;{\mathrm{sin}}\;\varGamma +\; \cos \;\theta_1 \;\cos \varGamma } \right) - \mathop v\nolimits_{2x} \left( {{\mathrm{sin}}\;\theta_2 \;{\mathrm{sin}}\;\varGamma + \cos \;\theta_2 \;\cos\; \varGamma } \right)} \right)} \right]^{\mathrm{T}}}. $

Carsim的计算结果为验证四自由度碰撞模型的精度提供了基准. 为了验证本文中碰撞模型的准确性，将基于Matlab/Simulink平台搭建的四自由度碰撞模型与Carsim的计算结果进行对比. 仿真中，施撞车辆和被撞车辆为相同型号的车辆，车辆参数来自于Carsim B-Class车辆模型，如表2所示.

仿真碰撞场景如图4所示，施撞车辆和被撞车辆的车速均为15 m/s，当发生碰撞时，施撞车辆和被撞车辆均具有一定的航向角$ {\theta }_{1} $和$ {\theta }_{2} $. 假设两车碰撞前的初始侧向速度、横摆角速度和侧倾角速度均为零，接触时间$ \mathrm{\Delta }t $为0.15 s，碰撞过程开始于1 s，结束于1.15 s，仿真结果如图5所示.

从验证结果可知，采用四自由度碰撞模型计算得到的碰后状态与Carsim计算的结果吻合较好，预测的横摆角速度存在一定的偏差，但在误差允许范围内.

上层决策层由轨迹库、轨迹预测和碰撞损伤评估模型组成. 其中轨迹库基于数值优化，在考虑车辆动力学约束的情况下，离线生成车辆特定速度下的可行驶轨迹库. 轨迹预测用来预测对方车辆未来短时间内的行驶轨迹. 碰撞损伤评估模型在线评估轨迹库中轨迹的碰撞时风险和碰后自车失稳风险,在决策时以极短时间从轨迹库中确定出最优轨迹，并将其传给下层运动控制层，碰撞损伤最小化策略的结构图如图6所示.

轨迹库是离线生成的，在碰撞损伤最小化系统运行前就已创建好. 当碰撞损伤最小化系统运行时，根据碰撞损伤评估模型在轨迹库中选取最优的一条轨迹，并传至下层控制器. 轨迹库中的轨迹均是通过数值求解最优控制问题生成的^[12]. 最优控制问题可以描述为

(31)$ \left.\begin{gathered} \mathop {\min }\limits_{{\boldsymbol{x}},{\boldsymbol{u}}} \;\;\int_0^{{t_{\mathrm{f}}}} {{{[{{({{\boldsymbol{X}}_{\mathrm{g}}} - {\boldsymbol{X}}(t))}^2}+{{({{\boldsymbol{Y}}_{\mathrm{g}}} - {\boldsymbol{Y}}(t))}^2}]}^2}{\mathrm{d}}t+{v_x}({t_{\mathrm{f}}})} ; \\ {\mathrm{s}}.{\mathrm{t}}.\;\;\;\;\dot {\boldsymbol{x}}(t) = f({\boldsymbol{x}}(t),{\boldsymbol{u}}(t)), \\ \;\;\;\;\;\;\;\;G({\boldsymbol{x}}(t),{\boldsymbol{u}}(t)) \leqslant 0, \\ \;\;\;\;\;\;\;\;{\boldsymbol{x}}(0) = {{\boldsymbol{x}}_0}, \\ \;\;\;\;\;\;\;\;{\boldsymbol{x}}({t_{\mathrm{f}}}) = {{\boldsymbol{x}}_{\mathrm{T}}}. \\ \end{gathered} \right\}$

式中：$ {\boldsymbol{x}} $为状态量，$ {\boldsymbol{u}} $为控制量，$ f({\boldsymbol{x}},{\boldsymbol{u}}) $为车辆模型，下标0和T分别表示车辆的初始状态和最终状态，$ {{\boldsymbol{X}}}_{\mathrm{g}}、{{\boldsymbol{Y}}}_{\mathrm{g}} $分别为整车在全局坐标系下的目标位置，$ G\left({\boldsymbol{x}}\right(t),{\boldsymbol{u}}(t\left)\right) $为约束条件，$ {t}_{\mathrm{f}} $为最终时间.

轨迹库的生成采用三自由度车辆模型来描述车辆运动，整车的动力学模型描述如下：

(32)$ \left. \begin{gathered} m({{\dot v}_x} - {v_y}{\omega _z}) = {F_x}, \\ m({{\dot v}_y}+{v_x}{\omega _z}) = {F_{y{\mathrm{f}}}}+{F_{y{\mathrm{r}}}}, \\ {I_{zz}}{{\dot \omega }_z} = {l_{\rm{f}}}{F_{y{\mathrm{f}}}} - {l_{\rm{r}}}{F_{y{\mathrm{r}}}}+\Delta M, \\ \dot S_1 = {v_x}\cos\; \theta - {v_y}\sin\; \theta , \\ \dot S_2 = {v_x}\sin\; \theta +{v_y}\cos\; \theta . \\ \end{gathered} \right\} $

式中：$ S_1、S_2 $为车辆的全局坐标；$ \mathrm{\Delta }M $为依赖于车轮纵向力的附加横摆力矩；$ \theta $为车辆横摆角，认为其近似等于航向角，$ \dot{\theta }={\omega }_{z} $.

结合式(32)，整理可得$ f({\boldsymbol{x}},{\boldsymbol{u}}) $. 其中状态量$ {{\boldsymbol{x}}}={[{v}_{x},{v}_{y},\omega ,S_1,S_2]}^{{\mathrm{T}}} $，控制量$ {{\boldsymbol{u}}}={[\delta ,{F}_{x},\Delta M]}^{{\mathrm{T}}} $.

对于每一个目标位置，求解式(31)定义的最优控制问题. 在给定约束条件和成本的情况下，在车辆的初始位置和目标位置之间创建最优轨迹，如图7所示. 在不同的初始速度下，离线创建包含最优轨迹的集合，如图8所示. 当碰撞损伤最小化系统运行时，选择与当初车辆速度相对应的轨迹子集.

为了估计两车的碰撞损伤，对对方车辆进行轨迹预测. 自身车辆检测到即将发生碰撞到碰撞这一段时间常为1.0~2.0 s，针对这种短时间的轨迹预测，可以采用运动学模型. 常见的汽车运动学模型有恒速度模型、恒加速度模型、恒角速度模型、恒转向角速度模型和恒角速度加速度模型(constant turn rate and acceleration, CTRA)等. 根据本文的研究场景，在交叉路口选择恒角速度加速度模型进行轨迹预测^[15]. 车辆的位置、速度、加速度和航向角数学模型如下所示：

(33)$ \left[ \begin{gathered} x(t+T) \\ y(t+T) \\ {v_x}(t+T) \\ \theta (t+T) \\ a(t+T) \\ {\omega _z}(t+T) \\ \end{gathered} \right] = \left[ \begin{gathered} x(t) \\ y(t) \\ {v_x}(t) \\ \theta (t) \\ \;\;a \\ \;\;{\omega _z} \\ \end{gathered} \right]+\left[ \begin{gathered} \Delta x(T) \\ \Delta y(T) \\ \;\;\;Ta \\ \;\;\;T{\omega _z} \\ \;\;\;\;0 \\ \;\;\;\;0 \\ \end{gathered} \right]. $

式中：

(34)$ \begin{split} \Delta x(T) = &[({v_x}(t){\omega _z}+a\omega T)\sin\; (\theta (t)+{\omega _z}T)+a \;\times \\ & \cos\; (\theta (t)+{\omega _z}T) - {v_x}(t){\omega _z}\sin\; \theta (t) - a\;\times\\ &\cos\; \theta (t)]/{\omega _z}^2,\end{split} $

(35)$ \begin{split} \Delta y(T) =& [( - {v_x}(t){\omega _z} - a{\omega _z}T)\cos\; (\theta (t)+{\omega _z}T)+a\;\times \\&\sin\; (\theta (t)+{\omega _z}T)+{v_x}(t){\omega _z}\cos\; \theta (t) - a\;\times\\&\sin\; \theta (t)]/{\omega _z}^2. \end{split} $

其中$ x、y $分别为汽车全局坐标系下的$ X $轴和$ Y $轴位置，$ \theta $为汽车的航向角，$ {v}_{x} $为车体坐标系下的纵向速度，$ \omega $为汽车的横摆角速度，$ a $为汽车纵向的最大加速度.

车辆碰撞损伤主要取决于碰撞位置、碰撞时两车相对速度及碰后自车失稳风险，碰撞损伤评估模型综合考虑这三方面，以实现碰撞损伤的最小化.

由于发生碰撞的车辆位置对车辆碰撞和乘客伤害严重程度起着至关重要的作用，根据文献[16]的结果，建立车身碰撞位置与碰撞损伤之间的关系. 将代表被撞车辆车身的矩形划分为10个多边形，如图9所示. 其中，多边形10和5分别表示车辆的前部和后部. 车辆的左、右两侧分别分为4个多边形(1~4和6~9). 此外，一个或多个多边形组合得到碰撞位置，例如单个多边形2与碰撞位置P1相关联，多边形2和3一起组合得到碰撞位置P0.

定义多边形10，即车辆前方，作为碰撞多边形. 若对方车辆身上有任何点属于我方车辆的多边形10，则我方车辆的前部参与碰撞，即我方车辆撞击另一辆车，如图10所示. 通过分离轴定理，可以确定我方车辆的多边形10与对方车辆的哪些多边形发生碰撞，通过表3确定碰撞位置及碰撞损伤.

速度相关的指标常被用于评估碰撞损伤. 根据相关的碰撞数据分析^[17]可知，碰撞时刻两车相对速度$ \Delta V $作为可接受的车辆碰撞损伤衡量标准，被广泛采用. $ \Delta V $的计算公式为

(36)$ \Delta V = \sqrt {{{({V_{{e_X}}} - {V_{{o_X}}})}^2}+{{({V_{{e_Y}}} - {V_{{o_Y}}})}^2}} . $

式中：$ {V}_{{e}_{X}}、{V}_{{e}_{Y}}、{V}_{{o}_{X}}、{V}_{{o}_{Y}} $分别为自身车辆和对方车辆在全局坐标系下沿$ X $轴和$ Y $轴方向的速度分量.

车辆在遭受碰撞后往往会同时产生横摆运动，较大的横摆角速度极易影响碰后车辆的稳定性^[10]，导致二次碰撞甚至多次碰撞. 将影响碰后车辆稳定性的潜在风险定义为

(37)$ {\mathrm{{\rm{CSI}}}}(\omega ) = {\omega ^2}. $

式中：$ \omega $为预测出的碰后自车横摆角速度.

总碰撞损伤为

(38)$ {\rm{CSI}} = {k_1}{\rm{CSI}}(P)+{k_2}\Delta V+{k_3}{\rm{CSI}}(\omega ). $

式中：$ {k}_{1}、{k}_{2} $和$ {k}_{3} $分别为与碰撞位置、碰撞时两车相对速度和碰后自车失稳风险相关的碰撞损伤权重系数.

从上层决策层基于碰撞损伤评估模型得到最优轨迹，输出给下层运动控制层. 为了保证车辆在极限工况下轨迹跟踪的精度和稳定性，控制器采用基于分层控制架构的轨迹跟踪与横摆稳定性协同控制方法^[18].

(39)$ \left. \begin{gathered} m({{\dot v}_y}+{v_x}{\omega _z}) = {F_{y{\mathrm{f}}}}+{F_{y{\mathrm{r}}}}, \\ {I_{zz}}{{\dot \omega }_z} = {l_{\rm{f}}}{F_{y{\mathrm{f}}}} - {l_{\rm{r}}}{F_{y{\mathrm{r}}}}+\Delta M, \\ \dot S_1 = {v_x}\cos\; \theta - {v_y}\sin\; \theta , \\ \dot S_2 = {v_x}\sin\; \theta +{v_y}\cos\; \theta . \\ \end{gathered} \right\} $

基于模型预测控制算法，设计轨迹跟踪控制器. 根据式(16)、(39)，将车辆模型写为状态空间方程的形式:

(40)$ {\dot {\boldsymbol{x}}_{{\mathrm{t}}1}} = {{\boldsymbol{A}}_{\mathrm{c}}}{{\boldsymbol{x}}_{{\mathrm{t}}1}}+{{\boldsymbol{B}}_{\mathrm{c}}}{{\boldsymbol{u}}_{{\mathrm{t}}1}}+{{\boldsymbol{w}}_{\mathrm{c}}}. $

式中：${\boldsymbol{x}}_{{\mathrm{t}}1}={[{v}_{y},\omega ,X,Y]}^{{\mathrm{T}}} \text{，} $

$ {{\boldsymbol{A}}_{\mathrm{c}}} = \left[ \begin{gathered} \;\;\;\frac{{{C_{\rm{f}}}+{C_{\rm{r}}}}}{{m{v_x}}}\;\;\;\;\;\;\frac{{{l_{\rm{f}}}{C_{\rm{f}}} - {l_{\rm{r}}}{C_{\rm{r}}}}}{{m{v_x}}} - {v_x}\;\;\;\;\;0\;\;\;\;\;0 \\ \frac{{{l_{\rm{f}}}{C_{\rm{f}}} - {l_{\rm{r}}}{C_{\rm{r}}}}}{{{I_{zz}}{v_x}}}\;\;\;\;\;\;\;\frac{{l_{\rm{f}}^2{C_{\rm{f}}}+l_{\rm{r}}^2{C_{\rm{r}}}}}{{{I_{zz}}{v_x}}}\;\;\;\;\;\;\;\;0\;\;\;\;\;0 \\ \;\;\; - \sin\; \theta \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;0 \\ \;\;\;\;\;\cos\; \theta \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\;\;\;\;\;0 \\ \end{gathered} \right] , $

$ {{\boldsymbol{B}}_{\mathrm{c}}} = \left[ \begin{gathered} \; - \frac{{{C_{\rm{f}}}}}{m}\;\;\;\;\;0 \\ - \frac{{{l_{\rm{f}}}{C_{\rm{f}}}}}{{{I_{zz}}}}\;\;\;\frac{1}{{{I_{zz}}}} \\ \;\;\;\;0\;\;\;\;\;\;\;0 \\ \;\;\;\;0\;\;\;\;\;\;\;0 \\ \end{gathered} \right] \text{，} {{\boldsymbol{u}}_{{\mathrm{t}}1}} = \left[ \begin{gathered} \;\;\delta \\ \Delta {M_z} \\ \end{gathered} \right] \text{，} {{\boldsymbol{w}}_{\mathrm{c}}} = \left[ \begin{gathered} \;\;\;\;\;0 \\ \;\;\;\;\;0 \\ {v_x}\cos\; \theta \\ {v_x}\sin\; \theta \\ \end{gathered} \right]. $

使用欧拉方法，对上述状态空间方程进行离散化处理，忽略高阶项，得到离散状态空间方程：

(41)$ \left. \begin{gathered} {{\boldsymbol{x}}_{{\mathrm{t}}1}}\left( {k+1} \right) = {{\boldsymbol{A}}_{{\mathrm{t}}1}}{{\boldsymbol{x}}_{{\mathrm{t}}1}}(k)+{{\boldsymbol{B}}_{{\mathrm{t}}1}}{{\boldsymbol{u}}_{{\mathrm{t}}1}}(k)+{{\boldsymbol{w}}_{{\mathrm{t}}1}}, \\ {{\boldsymbol{y}}_{{\mathrm{t}}1}}(k) = {{\boldsymbol{C}}_{{\mathrm{t}}1}}{{\boldsymbol{x}}_{{\mathrm{t}}1}}(k). \\ \end{gathered} \right\} $

式中：${\boldsymbol{A}}_{{\mathrm{t}}1}=\boldsymbol{I}+{\boldsymbol{A}}_{{\mathrm{c}}}{T}_{{\mathrm{s}}} $；$ {\boldsymbol{B}}_{{\mathrm{t}}1}={\boldsymbol{B}}_{{\mathrm{c}}}{T}_{{\mathrm{s}}} $；$ {\boldsymbol{w}}_{{\mathrm{t}}1}={\boldsymbol{w}}_{{\mathrm{c}}}{T}_{{\mathrm{s}}} $；$ {T}_{{\mathrm{s}}} $为仿真步长，$ {T}_{{\mathrm{s}}}=0.01 $；

$ {{\boldsymbol{C}}}_{\mathrm{t} 1}=\left[\begin{array}{llll}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 0 & 1\end{array}\right] ; $

输出向量$ {\boldsymbol{y}}_{{\mathrm{t}}1}\left(k\right)={[{v}_{y},{\omega }_{z},Y]}^{{\mathrm{T}}} $.

定义模型预测控制的目标函数为

(42)$ \begin{split} \min \;\;&\sum\limits_{j = 1}^N {[{{({{\boldsymbol{y}}_{{\mathrm{des}}}}(k+j) - {{\boldsymbol{y}}_{{\mathrm{t}}1}}(k+j))}^{\mathrm{T}}}{{\boldsymbol{Q}}_{{\mathrm{t}}1}}({{\boldsymbol{y}}_{{\mathrm{des}}}}(k+j) - } \\ &{{\boldsymbol{y}}_{{\mathrm{t}}1}}(k+j))]+\sum\limits_{j=0}^{N - 1} {[{{\boldsymbol{u}}_{{\mathrm{t}}1}}{{(k+j)}^{\mathrm{T}}}{{\boldsymbol{R}}_{{\mathrm{t}}1}}{{\boldsymbol{u}}_{{\mathrm{t}}1}}(k+j)]} . \end{split} $

式中：状态向量的权重$ {\boldsymbol{Q}}_{{\mathrm{t}}1}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\;[{q}_{{v}_{y}},{q}_{{\omega }_{z}},{q}_{Y}] $，输出的权重$ {\boldsymbol{R}}_{{\mathrm{t}}1}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\;[{r}_{\delta },{r}_{\Delta M}] $，期望向量$ {\boldsymbol{y}}_{\mathrm{d}\mathrm{e}\mathrm{s}}\left(k\right)= [{v}_{y\_\mathrm{d}\mathrm{e}\mathrm{s}}, {\omega }_{z\_{\mathrm{des}}},{Y}_{{\mathrm{des}}}]^{{\mathrm{T}}} $.

转矩优化分配控制器是在多约束条件下求解最优轮胎纵向力分配，以满足附加横摆力矩$ \Delta M $和总纵向力$ {F}_{x} $的要求. 根据轮胎附着圆，考虑轮胎利用率最低，建立目标函数：

(43)$ J = \sum\limits_{ij = {\mathrm{fl,fr,rl,rr}}} {\frac{{F_{xij}^2+F_{yij}^2}}{{{{(\mu {F_{zij}})}^2}}}} . $

考虑到$ {F}_{yij} $对于优化问题而言是常数，将优化目标函数简化为

(44)$ J = \sum\limits_{ij = {\mathrm{fl,fr,rl,rr}}} {\frac{{F_{xij}^2}}{{{{(\mu {F_{zij}})}^2}}}} . $

由于优化问题的有效解要满足车辆纵向车速与直接横摆力矩的需求，可得

(45)$ \left. \begin{split} &{F_{x{\mathrm{fl}}}}\cos\; \delta +{F_{x{\mathrm{fr}}}}\cos\; \delta +{F_{x{\mathrm{rl}}}}+{F_{x{\mathrm{rr}}}} = {F_x}, \\& 0.5b({F_{x{\mathrm{fr}}}} - {F_{x{\mathrm{fl}}}})\cos\; \delta +0.5b({F_{x{\mathrm{rr}}}} - {F_{x{\mathrm{rl}}}})+ \\&\qquad {l_{\rm{f}}}({F_{x{\mathrm{fl}}}} +{F_{x{\mathrm{fr}}}})\sin\; \delta {\text{ = }}\Delta M. \end{split} \right\} $

式中：$ {F}_{x} $为保证当前车辆纵向速度的总纵向力，可以根据比例控制得到$ {F}_{x}={K}_{{\mathrm{p}}}({v}_{x\_{\mathrm{des}}}-{v}_{x}) $.

驱动电机的最大输出转矩和路面附着条件决定了车轮的最大纵向力，有

(46)$ \left| {{F_{xij}}} \right| \leqslant \min \;\{\left| {\mu {F_{zij}}} \right|,\left| {{T_{\max }}/R} \right|\}. $

综合式(43)~(45)，得到下层转矩优化分配控制器：

(47)$ \begin{split} & \min _{\boldsymbol{u}_{\min } \leqslant \boldsymbol{u}_{\mathrm{t} 2} \leqslant \boldsymbol{u}_{\max }}\left\|{{\boldsymbol{Q}}}_{\mathrm{u}}\left(\boldsymbol{u}_{\mathrm{t} 2} - \boldsymbol{u}_{\mathrm{d}}\right)\right\|^2 + k_{\mathrm{t} 2}\left\|\boldsymbol{R}_{\mathrm{v}}\left(\boldsymbol{B}_{\mathrm{t} 2} \boldsymbol{u}_{\mathrm{t} 2} - \boldsymbol{v}_{\mathrm{t} 2}\right)\right\|^2 = \\& \min _{\boldsymbol{u}_{\min } \leqslant \boldsymbol{u}_{\mathrm{t} 2} \leqslant \boldsymbol{u}_{\max }}\left\|\underbrace{\binom{\sqrt{k_{\mathrm{t} 2}} \boldsymbol{R}_{\mathrm{v}} \boldsymbol{B}}{\boldsymbol{Q}_{\mathrm{u}}}}_{\boldsymbol{A}_{\mathrm{t} 2}} \boldsymbol{u}_{\mathrm{t} 2}-\underbrace{\binom{\sqrt{k_{\mathrm{t} 2}} \boldsymbol{R}_{\mathrm{v}} \boldsymbol{v}_{\mathrm{t} 2}}{\boldsymbol{Q}_{\mathrm{u}} \boldsymbol{u}_{\mathrm{d}}}}_{\boldsymbol{b}_{\mathrm{t} 2}}\right\|^2= \\& \min _{\boldsymbol{u}_{\min } \leqslant \boldsymbol{u}_{\mathrm{l} 2} \leqslant \boldsymbol{u}_{\max }}\left\|\boldsymbol{A}_{\mathrm{t} 2} \boldsymbol{u}_{\mathrm{t} 2}-\boldsymbol{b}_{\mathrm{t} 2}\right\|^2 .\\[-1pt]\end{split} $

式中：

$ {\boldsymbol{u}}_{\mathrm{t}2}={[{{F}}_{{x}\mathrm{f}\mathrm{l}},{{F}}_{{x}\mathrm{f}\mathrm{r}},{{F}}_{{x}\mathrm{r}\mathrm{l}},{{F}}_{{x}\mathrm{r}\mathrm{r}}]}^{\mathrm{T}}\text{；}{{u}}_{\mathrm{d}}=0 \text{；} $

$ {\boldsymbol{v}}_{\mathrm{t}2}={[{{F}}_{{x}},\mathrm{\Delta }{M}]}^{\mathrm{T}} \text{，} {\boldsymbol{R}}_{\mathrm{v}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\;[{{r}}_{\mathrm{v}}{{F}}_{{x}},{{r}}_{\mathrm{v}}\mathrm{\Delta }{M}] \text{；} $

$\begin{array}{l}{{\boldsymbol{B}}_{{\rm{t}}2}} = \left[ {\begin{array}{*{20}{c}}{\cos {\delta _{\rm{f}}}} & {\cos {\delta _{\rm{f}}}} & 1 & 1\\{ - \dfrac{{b\cos {\delta _{\rm{f}}}}}{2} + {l_{\rm{f}}}\sin {\delta _{\rm{f}}}} & {\dfrac{{b\cos {\delta _{\rm{f}}}}}{2} + {l_{\rm{f}}}\sin {\delta _{\rm{f}}}} & { - \dfrac{b}{2}} & {\dfrac{b}{2}}\end{array}} \right], \\{{\boldsymbol{Q}}_{\rm{u}}} = {\mathop{\rm diag}\nolimits} \left[ {\dfrac{{{c_{{\rm{fl}}}}}}{{\mu {F_{z{\rm{fl}}}}}},\;\dfrac{{{c_{{\rm{fr}}}}}}{{\mu {F_{z{\rm{fr}}}}}},\;\dfrac{{{c_{{\rm{rl}}}}}}{{\mu {F_{z{\rm{rl}}}}}},\;\dfrac{{{c_{{\rm{rr}}}}}}{{\mu {F_{z{\rm{rr}}}}}}} \right].\end{array}$

$ {{\boldsymbol{R}}}_{\mathrm{v}}、{{\boldsymbol{Q}}}_{\mathrm{u}} $为加权矩阵；$ {c}_{ij} $分别为4个车轮纵向力在优化目标函数中的权重，$ ij\in \{{\mathrm{fl}},{\mathrm{fr}},{\mathrm{rl}},{\mathrm{rr}}\} $分别指左前轮、右前轮、左后轮和右后轮；$ {{r}_{{\mathrm{v}}}\mathrm{、}k}_{{\mathrm{t}}2} $为权重系数.

为了验证提出的碰撞损伤最小化算法，基于Matlab/Simulink平台搭建车辆模型和控制器模型，参考文献[16]模拟和分析2个场景. 第1个场景是T型路口，涉及2辆车. 随着车辆数的增加，碰撞损伤评估模型将变得更加复杂，所以第2个场景涉及多辆汽车.

该场景具体如图11所示，场景中自车以20 m/s的速度直行通过路口，它车以20 m/s的速度行驶进入T型路口. 自车通过CTRA模型预测出它车即将左转进入内侧车道的轨迹. 当自车检测到20 m/s对应的轨迹库均与所预测的它车轨迹发生碰撞，自车触发碰撞损伤最小化系统.

图12给出自车选择不同轨迹与它车碰撞的预测结果. 利用碰撞损伤评估模型，计算各条轨迹的碰撞损伤. 如表4所示为部分代表性轨迹，其中主要碰撞位置有5个，综合考虑碰撞位置、碰撞时两车相对速度和碰后自车失稳风险，确定出总碰撞损伤最小的为136号轨迹，仿真结果如图13所示. 图中，M为制动力矩.

多车场景具体如图14所示，场景中自车以15 m/s的速度直行通过路口，自南向北左转的它车1以15 m/s的速度行驶进入十字交叉路口，自北向南左转的它车2以15 m/s的速度行驶进入十字交叉路口. 自车通过CTRA模型预测出它车1、2的未来行驶轨迹. 当自车检测到15 m/s对应的离线轨迹库均与所预测的它车轨迹发生碰撞时，自车触发碰撞损伤最小化系统.

图15(a)、(b)分别给出当自车选择不同轨迹时，与它车1、2碰撞的预测结果. 其中自车与它车1的主要碰撞位置为B0和P2，与它车2的主要碰撞位置为P1、P2和F0. 通过碰撞损伤评估模型计算各条轨迹的碰撞损伤，如表5所示为部分代表性轨迹. 综合考虑碰撞位置、碰撞时两车相对速度和碰后自车失稳风险，确定碰撞损伤最小的为78号轨迹，仿真结果如图16所示.

(1)采用矢量方式，将各种碰撞情形统一为一种建模方式，建立四自由度碰撞模型. 通过四自由度碰撞模型和碰撞前的车辆状态，估计碰后自车失稳风险. 碰撞损伤评估模型综合考虑碰撞位置、碰撞时两车相对速度和碰后自车失稳风险，以尽可能减小车辆碰撞损伤.

(2)利用附加横摆力矩能够改变车辆横摆角速度和航向角的特点，在轨迹库的生成中引入附加横摆力矩以提供更大的潜力和灵活性，充分发挥出轮胎纵向的力学性能，而不仅仅基于前轮极限转向角和制动率.

(3)为了保证最终的轨迹跟踪精度和车辆稳定性，基于模型预测控制建立轨迹跟踪与横摆稳定性协同控制器，实现多目标综合控制.