极值服从广义Pareto分布的扭转载荷外推方法研究

郑国峰，陈柏先，隗寒冰，杨昊民，严璐瑶

（1.重庆交通大学机电与车辆工程学院，重庆 400074；2.中国汽车工程研究院股份有限公司，重庆 401122）

0 引言

全寿命周期内的服役载荷谱，对零部件疲劳耐久性的准确开发具有重要意义。但对零部件全寿命周期下部件的服役载荷进行采集，需要耗费大量成本。考虑到零部件在典型工况下的服役载荷具有相应的统计规律，可通过掌握部分随机载荷统计规律，再基于统计规律进行外推，实现对未采集载荷谱的掌握。载荷谱外推的本质就是通过采集载荷谱的统计规律预测相同工况下未采集载荷谱。

时域外推法认为载荷谱时间序列中超过阈值的极大值和极小值服从某种分布，基于概率密度分布函数可以实现对时间序列规律的掌握和预测。杨子涵等［1］提出基于经验模态分解－峰值超阈值模型的载荷谱外推方法，通过对外推负荷筛选、阈值限定和极值排序，解决了传统时域外推方法在载荷外推时结果失真的问题。He等［2］对数控机床时域载荷外推方法中阈值的求解问题，提出了一种基于灰色关联分析的最优阈值确定模型，进一步提高数控机床载荷谱精度。郑国峰等［3］对广义Pareto分布模型形状参数的2种情况进行研究，发现形状参数为零时的载荷外推结果比形状参数不为零时的结果更为保守。Wang等［4］针对轨道车辆动态应力数据进行研究，得出一种轨道车辆长时间服役条件下疲劳寿命的评估方法，即通过极值理论对动态应力进行外推，发现当采集载荷过小时，会造成一定的误差。陈晋市等［5］提出将主成分分析与多种拟合检验准则相结合选取最优阈值的载荷谱外推方法，用于挖掘机主泵载荷外推研究，但受主成分分析方法影响较大。Yang等［6］提出一种改进的阈值选择方法，即分析极值样本的平均值来获取最佳阈值对拖拉机时域载荷进行外推，其结果具有较高的精度，同时保留了载荷谱的时频特性。基于以上问题，本文以均方误差为阈值选择标准，对商用车驾驶室稳定杆扭转载荷进行时域外推，该方法受采集载荷样本的影响较小，同时能够保留载荷谱的统计特性。对于时域外推法中载荷谱的阈值选取，常用的选取方法有经验法、灰色关联度法［2，7－8］、均方误差法［9］等。经验法需要依靠工程师的经验进行选取，具有很大的局限性。灰色关联度法是根据选择不同阈值所得超出量拟合曲线与样本的分布曲线的接近程度的研究方法，具有较高的精确性；均方误差法是运用统计学抽样来模拟样本中的极值，通过计算均方误差进行选择的方法。本文采用均方误差法进行阈值选取。

本文中提出了一种基于广义Pareto分布（generalized pareto distribution，GPD）函数的扭转载荷时域外推方法。本文介绍广义Pareto模型（Pareto模型为一种数学分布模型，即极值样本服从Pareto分布）的建立过程，同时以商用车驾驶室稳定杆的扭转随机载荷谱为外推对象，对比分析了外推前后载荷谱的穿级计数、功率谱密度、雨流图和潜在伪损伤量。结果表明：所构建的载荷外推算法能够准确地实现扭转载荷的外推。

1 基于GPD的扭转载荷外推理论

扭转载荷是部件发生疲劳失效的一种重要载荷形式。汽车结构中的典型部件（如稳定杆、传动轴等）在服役过程中，在扭转载荷作用下，其内部将呈现出交变切向应力，在较大幅值的应力长时间作用后，结构发生疲劳失效。与受水平载荷的典型部件不同，受纯扭转载荷的部件更多出现的是剪切疲劳失效，而非拉压疲劳失效。因此，部件全寿命周期内的扭转随机疲劳载荷谱是研究结构疲劳性能的关键。可根据采集的扭转随机载荷的统计规律，对未知载荷谱进行预测，实现目标里程下载荷谱的获取，完成部件耐久性评估。

1.1 GPD函数及其性质

对部件所受到的随机载荷谱进行外推前，需要先进行载荷谱的预处理，如载荷的峰谷值提取。对于载荷谱中峰谷值超过某一阈值的部分，定义为载荷谱的极值（也叫超出量），见图1。研究表明：随机载荷谱的极值服从GPD分布［3］。广义Pareto分布函数表达式为［10］：

式中：u为位置参数（实为所设置的阈值）；σ为尺度参数；ξ为形状参数。其概率密度函数为

根据图1，Xj（j＝1，2，3，…，N）表示提取峰谷值后的采样点，N为载荷谱峰谷值的采样点总数；设定阈值为u，其中采样点中超过阈值的采样点定义为极值样本，记作Yk＝Xj－u（j＝1，2，3，…，N；k为样本序列号，j为采集点序列号）。将载荷谱中的极值样本分为上极值样本和下极值样本，即在载荷谱中提取出载荷谱的峰谷值，然后设定上阈值umax和下阈值umin，峰谷值中超过上阈值的部分就是载荷谱的上极值样本，超过下阈值的部分就是载荷谱的下极值样本，上下极值样本可分别表示为

式中：为上极值样本：为下极值样本。

1.2 GPD函数阈值

基于广义Pareto分布的载荷谱外推，阈值的选择决定了广义Pareto分布函数的参数，进而影响载荷谱外推结果。载荷谱阈值的确定首先需要基于极值样本的均值超出函数确定一个取值区间，然后以形状参数最小的均方误差为目标，确定最佳阈值［9］。

1.2.1 经验均值的超出函数

通过经验均值超出函数初步确定一个阈值区间范围，以形状参数最小的均方误差为目标，在该区间内通过增量步迭代，实现最优阈值的确定。经验均值超出函数为［4］

式中：Xi是极值样本；Nu是样本中的数量。根据GPD函数性质，超出量均值函数为

阈值区间的确定，需要绘制经验均值超出函数en（u）的散点图，随后根据式（6）可知，当形状参数和尺度参数一定时，经验均值超出函数和阈值具有线性关系。故在经验均值超出函数图中找到呈线性关系的区间，则该区间为阈值范围区间。

1.2.2 GPD函数最优阈值的确定

基于GPD函数的极值样本的矩分析，可以得到Pareto分布函数中形状参数与矩之间的关系表达式，如下：

由式（7）可知，GPD函数的形状参数ξ与所选定的阈值u相关。当阈值u偏大时，极值样本的数量将会过少，导致样本分布不稳定；当阈值u偏小时，极值样本的数量过多，将导致所得结果极大偏离估计值和真实值。因此，采用均方误差MSE（Mean square error）描述阈值u的选择情况。均方误差综合考虑阈值偏大或者偏小所带来的误差。均方误差通过如下公式获取［4］。

式中：Bia2（ξ）是形状参数偏差；Var（ξ）是形状参数方差。

形状参数偏差用于描述估计值与实际值之前的巨大偏差，通过如下公式计算。

式中：b为抽样次数。形状参数的方差用于描述样本分布的稳定性，通过如下公式计算。

采用均方误差估计阈值u时，使用自抽样方法对极值样本进行多次重复抽样［11］，得到自抽样法样本组{B}i，进而得到多个相对应的形状参数估计值。其中，最小形状参数均方误差估计值所对应的阈值为最佳阈值。

综上，GPD函数最优阈值的确定流程可总结如下（见图2）：

图2 阈值确定流程

1）在经验超出量函数确定的阈值区间中选择一个阈值作为初始阈值u0，得到极值样本，并记录极值样本中的数量Nj。

2）计算出形状参数的初始估计值，记为ξ0。

3）使用自抽样法对极值样本随机可重复选择Nj个样本数据，得到一组自抽样样本Bi。

4）计算出自抽样样本Bi对应的形状参数估计值ξi。

5）重复步骤3）和4），得到b组自抽样样本所对应的形状参数估计值组｛ξi｝。

6）计算自抽样样本组对应的Bia（ξ）和Var（ξ）的估计量，然后运用式（8）计算MSE（ξi）。

7）重置初始阈值uj＝uj＋j·Δu（j＝1，2，…N），重复步骤1）—步骤6），计算得到阈值样本组｛uj｝和对应的均方误差估计值组｛MSE（ξi）｝。

8）均方误差估计值组｛MSE（ξi）｝中的最小值对应的阈值uj为最佳阈值。

步骤5）中的抽样次数b对阈值的分析结果有显著影响。当抽样次数太少，计算结果将不稳定；而抽样次数超过一定量时，计算结果将趋于稳定。因此，要保证阈值的精确度，需要达到一定抽样次数。

1.3 GPD函数形状参数和尺度参数估计

GPD函数的形状参数和尺度参数的常用方法有：矩估计法［12］、极大似然估计法［13］、概率加权矩法［14］等。其中，基于极大似然的估计法，首先对含有未知参数的分布函数建立对数似然函数，对似然函数求导并等于零，求解方程得到估计参数的估计量。具体而言，对极值样本建立对数似然函数：

基于式（11）对形状参数和尺度参数求偏导，得到关于形状参数和尺度参数的方程组。

求解以上方程组，便可获取形状参数和尺度参数的估计值。

1.4 基于GPD的载荷谱外推流程

基于GPD函数的载荷谱外推的步骤总结如下：

1）输入采集载荷谱并提取载荷谱峰谷值Xj。

2）选择合适的上阈值umax和下阈值umin。

3）提取上阈值超出量Ymaxk和下阈值超出量Ymink。

4）采用GPD函数对超出量Yk进行拟合，获得GPD函数的参数。

5）将拟合得到的参数还原到GPD函数中，然后通过该函数随机生成数量相等的超出量。

6）生成的超出量替换原始载荷谱中的超出部分，从而获得外推信号。

2 载荷谱采集与外推

2.1 载荷谱采集

商用车驾驶室在转向等受到侧向力的工况时，驾驶室的稳定杆能够防止车身发生过大横向侧倾，使驾驶室保持平衡。故驾驶室稳定杆的本质是一个横置的扭杆弹簧，在功能上可以看成是一种特殊的弹性元件，其服役过程中主要受到扭转载荷的作用。以商用车驾驶室稳定杆为对象，首先介绍其扭转载荷采集的方法，然后进行扭转载荷外推研究与验证。

2.2 稳定杆载荷谱采集

在进行驾驶室稳定杆扭转载荷采集之前，需要分别进行贴片、组桥和标定工作，见图3。

图3 驾驶室稳定杆扭转载荷采集应变片粘贴、组桥与标定

将贴片、标定好后的驾驶室稳定杆装车，在典型的驾驶工况上进行载荷谱采集，得到100 km的载荷谱数据，见图4。

图4 用户典型工况100 km采集载荷谱数据

为了对外推算法进行验证，将采集的100 km用户城市工况载荷谱数据分为两部分，一部分载荷谱的里程长度为50 km，用于基于广义Pareto分布的扭转载荷外推；另一部分载荷谱的里程长度为50 km，用于与外推扭转载荷进行对比，验证载荷外推的效果。根据图4，驾驶室稳定杆前50 km载荷数据中最大载荷为960.4 N·m，最小载荷为－1 378 N·m；驾驶室稳定杆后50 km载荷数据中最大载荷为1 378 N·m，最小载荷为－1 192 N·m。

2.3 GPD函数参数拟合

2.3.1 阈值确定

根据式（8），采用均方误差MSE来确定阈值u，但首先需要经验均值的超出函数初步确定一个阈值区间范围，见图5。

图5 上下阈值的超出量分布图

根据试验载荷样本，初步确定上阈值的区间范围为350～400 N·m，下阈值的区间范围为－330～－380 N·m，见图6。初步确定上下阈值区间后，使用MSE来描述上下阈值区间内的阈值选择情况，上下阈值对应的MSE分别见图6（a）与图6（b）。

图6 上下阈值及其对应均方误差曲线

根据图6，在指定的阈值区间内，每个阈值所对应的MSE随阈值的增大而整体呈现波动上升的趋势，但分别在351 N·m和－331 N·m时取得最小值。因此，最终确定上阈值为351 N·m，下阈值为－331 N·m。

2.3.2 参数拟合及评价

基于确定的上、下阈值，从驾驶室稳定杆扭转载荷数据样本中，分别提取上、下阈值极值样本并进行拟合，见图7。

图7 峰谷值极值样本GPD函数拟合

根据图7，GPD函数与峰谷值极值样本拟合时，较小出现偏离参考线的情况，说明GPD函数拟合情况较为精确，具有较小的拟合参数偏差。为了对拟合结果进行量化评价，采用拟合优度对拟合结果进行评判［15］。

拟合优度是指回归曲线对离散数据的拟合程度，判定系数R2度量回归曲线与离散数据之间的吻合程度。R2的取值范围为［0，1］，当R2值越接近1，说明回归曲线对离散数据的拟合程度越好；反之，说明回归曲线对离散数据的拟合程度越差。R2值为0.7～1，相关性强；R2值为0.5～0.7，相关性一般；R2值为0～0.5，相关性弱［16］。

表1 基于峰谷值超出量的GPD函数拟合参数

根据表1，峰值极值样本的相关系数R2值为0.875 2，谷值极值样本的相关系数R2值为0.793，说明通过拟合得到峰谷值超出量的形状参数和尺度参数，与极值样本具有强相关性，能够保留极值样本的统计学特性。

3 载荷谱外推结果与分析

3.1 载荷谱外推

根据载荷谱外推步骤，通过GPD函数模型随机生成数量相等的峰值和谷值超出量，并替换原始载荷谱中的相应峰谷值，获得外推载荷谱。对驾驶室稳定杆前50 km扭转载荷进行外推，得到的外推扭转载荷见图8。

图8 基于GPD函数的外推载荷谱

根据图8，基于GPD函数外推得到的载荷谱，其极值的绝对值与原始载荷谱相比变大。极大值从原始数值960.4 N·m增大为1 939.6 N·m，极小值从原始数值－1 378 N·m增大为－1 409.4 N·m。

3.2 载荷谱外推结果分析

为进一步对比扭转载荷的外推效果，从穿级计数、雨流计数、功率谱密度以及伪损伤分析等方面，对外推前后载荷谱进行对比，验证所提出的外推算法。

3.2.1 穿级计数分析

对前50 km采集到的原始载荷谱、GPD函数外推得到的载荷谱，以及后50 km采集到的原始载荷谱进行穿级计数对比分析，观察其幅值和频次的变化情况，见图9。

图9 原始载荷谱与外推载荷谱穿级计数曲线

根据图9，前50 km采集的原始载荷谱在基于GPD函数外推后，所得到的外推载荷谱在幅值较小处的频次基本吻合，在幅值较大处的频次有所增加。主要原因：原始载荷谱外推后其极值的绝对值大于原始载荷谱极值的绝对值，由于外推载荷极值大小生成具有随机性，是生成少量较大极值导致。

基于GPD函数外推得到的载荷谱，在负值扭转载荷中与后50 km采集的原始载荷谱更加吻合。但在正值扭转载荷中，外推载荷谱与实际采集的载荷谱之间具有一定的偏差，主要体现在：外推载荷谱出现了少量较大幅值的循环，而采集的后50 km原始载荷谱中间幅值的循环较多。

3.2.2 雨流分析

对驾驶室稳定杆前50 km 原始载荷谱、前50 km外推载荷谱和后50 km原始载荷谱进行雨流计数，绘制出From-to雨流图，见图10。

图10 原始载荷谱与外推载荷谱雨流计数

根据图10（a）和（b），前50 km采集的原始载荷谱相比后50 km采集的原始载荷谱，大幅值信号明显较少。图10（c）的外推载荷谱雨流图，有一个大幅值循环出现，表明基于GPD函数的外推得到的外推载荷谱有较大幅值的出现。总体来看，前50 km GPD函数外推得到的外推载荷谱与后50 km采集到的原始载荷谱的效果更加一致。

3.2.3 功率谱密度分析

将前50 km外推载荷谱与后50 km原始载荷谱进行能量变化分析，将其功率谱密度进行对比。

根据图11，在1 Hz处出现峰值，由振动理论可知，该频率为驾驶室稳定杆的固有频率。在20 Hz之前，前50 km外推载荷谱能量与后50 km原始载荷谱能量基本吻合；但20 Hz之后，前50 km外推载荷谱能量与后50 km原始载荷谱能量有一定差异，主要原因在于外推后的载荷谱中增加了部分大幅值信号，而低频对应小幅值，高频对应大幅值，因此使得外推后的载荷谱功率谱密度的高频部分高于原始载荷谱功率谱密度。但外推前后载荷谱的能量主要集中在前20 Hz，因此外推后的载荷谱能量与原始载荷谱基本一致。

图11 外推前后功率谱密度变化曲线

3.2.4 伪损伤分析

为了进一步验证载荷外推方法的准确性，对前50 km原始载荷谱、前50 km外推载荷谱和后50 km原始载荷谱进行伪损伤分析。运用ncode软件进行伪损伤分析，商用车扭杆的材料相关参数见表2。

表2 商用车扭杆的材料属性

分析结果如表3所示，前50 km外推载荷谱潜在伪损伤量略大于前50 km原始载荷谱，原因是外推后载荷谱的载荷循环次数略大于原始载荷谱；后50 km原始载荷谱潜在伪损伤量小于前50 km外推载荷谱，但循环次数大于前50 km原始载荷谱，其原因为后50 km原始载荷循环在小应力幅处的循环次数高于前50 km外推载荷，在大应力幅值处的循环次数低于前50 km外推载荷。此外，前50 km、后50 km原始载荷谱的循环次数与路面状况的随机性存在一定关系。后50 km原始载荷谱与外推载荷谱潜在伪损伤量差异性较小，表明GPD函数载荷外推在误差允许范围内，具有较好的结果。

表3 3种载荷潜在伪损伤量计算结果

4 结论

通过所构建的载荷外推算法进行外推，其外推载荷谱极值的绝对值均大于原始载荷谱极值的绝对值。外推后的载荷谱中增加了部分大幅值信号，而低频对应小幅值，高频对应大幅值，因此使得外推后的载荷谱功率谱密度的高频部分高于原始载荷谱功率谱密度，外推载荷潜在伪损伤量与原始载荷潜在伪损伤量差异较小。故基于GPD函数外推算法，所得外推结果与原始载荷谱基本一致，GPD函数能够准确实现商用车驾驶室稳定杆扭转载荷外推。因此，本文外推算法对商用车作业时驾驶室稳定杆承受扭转载荷具有较好适应性。