基于局部自适应带宽扩散核密度估计的载荷外推

Load extrapolation based on local adaptive bandwidth diffusion kernel density estimation

摘要

为实现由有限载荷数据外推得到全周期载荷谱，针对传统自适应带宽优化算法在参数选择上存在局限性，本文提出一种基于局部自适应带宽扩散核密度估计的载荷外推方法。该方法首先将二维雨流矩阵降维至一维等效幅值，然后基于局部积分均方误差优化局部带宽，利用局部最优带宽，通过扩散核密度估计构建概率密度分布模型，最后结合蒙特卡洛模拟方法外推目标频次载荷。对某特种车辆综合传动装置的预处理载荷数据进行对比验证，结果表明：与传统方法相比，本文所提方法得到的概率密度分布曲线和累计频次曲线更接近实际等效幅值，相关系数与决定系数均更趋近于1，其中相关系数分别为0.983 8与0.999 6，决定系数分别为0.967 9与0.999 1；均方根误差也更小，分别为5.05×10^－5与15.9。

Abstract

To achieve the extrapolation of the full-cycle load spectrum from limited load data， this paper proposes a load extrapolation method based on local adaptive bandwidth diffusion kernel density estimation， in view of the limitations of traditional adaptive bandwidth optimization algorithms in parameter selection. This method first reduces the two-dimensional rainflow matrix to one-dimensional equivalent amplitude， then optimizes the local bandwidth based on the local integral mean square error， and uses the local optimal bandwidth to construct a probability density distribution model through diffusion kernel density estimation. Finally， the target frequency load is extrapolated by combining the Monte Carlo simulation method. The preprocessed load data of a certain special vehicle's comprehensive transmission device are compared and verified. The results show that compared with the traditional method， the probability density distribution curve and cumulative frequency curve obtained by the method proposed in this paper are closer to the actual equivalent amplitude. The correlation coefficients and determination coefficients are both closer to 1， with the correlation coefficients being 0.983 8 and 0.999 6 respectively， and the determination coefficients being 0.967 9 and 0.999 1 respectively. The root mean square errors are also smaller， being 5.05×10^-5 and 15.9 respectively.

Graphical abstract

关键词

机械工程 / 特种车辆 / 综合传动装置 / 载荷谱 / 载荷外推 / 局部自适应带宽 / 扩散核密度估计

Key words

mechanical engineering / special vehicle / integrated transmission device / load spectrum / load extrapolation / local adaptive bandwidth / diffusion kernel density estimation

0 引 言

载荷谱是零部件抗疲劳优化设计、剩余疲劳寿命预测与可靠性设计试验的重要依据^［1-4］。特种车辆的传动装置与工作需求不同于一般车辆，其行驶载荷更复杂多变、随机性更强，且载荷数据获取较困难。在现役特种车辆上加装传感器采集全寿命周期内零部件的载荷数据，将会造成极大的资源浪费^［5］。典型工况下的载荷数据具有一定的统计规律，基于有限载荷数据的统计规律对未知载荷数据进行外推，能够有效推动零部件研发设计进程。因此，开展载荷外推方法的研究，具有重要的理论意义和工程应用价值。

随着数理统计技术的发展，载荷外推方法逐渐形成参数法载荷外推、非参数法载荷外推与时域载荷外推3种形式^［6，7］。参数法基于过往经验，通过假设变量属于某一已知分布模型计算总体的概率分布，常用的分布模型有威布尔分布与正态分布^［8-10］；非参数法利用非参数估计方法直接得到数据的概率分布^［11-13］；时域载荷外推利用一维分布拟合时域载荷数据，由具体的一维概率密度函数实现载荷外推^{［14，15］}。由于非参数法适用于随机性较强的载荷，且无需进行分布假设，越来越受学者们的青睐。

核密度估计是常用的非参数估计方法，其精度主要取决于带宽的选择。固定带宽对所有数据采用相同的取值，无法单独调节。由于载荷数据具有随机性，采用固定带宽计算得到的概率密度分布可能与实际样本分布相差较大，而能够根据数据密集程度变化的自适应带宽，可在一定程度上提高估计值的准确性。传统自适应带宽采用不同的优化算法修正固定带宽，而优化算法本身需要进行参数选择，若参数选择不当，会降低估计精度^［16］。

为此，本文提出一种基于局部自适应带宽的扩散核密度估计算法，并将其应用于载荷外推。该算法首先对二维雨流矩阵进行降维，然后以局部积分均方误差（Local mean integrated squared error，LMISE）作为局部损失函数，通过求解局部最小损失函数优化局部带宽。基于最优局部带宽进行扩散核密度估计，从而得到载荷的概率密度分布模型，最后结合蒙特卡洛模拟（Monte Carlo，MC）算法实现目标频次的载荷外推。以某特种车辆综合传动装置为研究对象，通过对比多种传统载荷外推方法，验证该方法的准确性与有效性。

1 非参数估计方法

1.1　核密度估计

核密度估计是一种由数据自身分布情况决定的非参数估计方法，该方法在未知数据先验知识的情况下，完全基于样本数据的分布规律推断概率密度函数，进而建立总体分布模型，突破了参数法中对样本母体分布的过度依赖，对随机载荷数据具有更强的适应性^［17-19］。

假设x₁，x₂，x₃，…，x_n 是随机变量x的样本数据，且其概率密度函数为f（x），则利用核密度估计得到关于x的概率密度函数估计量为：

F\left(x;h\right)=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}K\left(x,x_i;h\right)=\frac{\mathrm{1}}{nh}{\displaystyle \sum _{i=\mathrm{1}}^n}K\left(\frac{x-x_i}{h}\right)（1）

式中：h为带宽；n为样本量；K（·）为核函数。

核函数需满足归一性、对称性与非负性，表达式为：

\left\{\begin{array}{l}\int K\left(u\right)\mathrm{d}u=\mathrm{1}\\ \int uK\left(u\right)\mathrm{d}u=\mathrm{0}\\ \int u^{\mathrm{2}}K\left(u\right)\mathrm{d}u>\mathrm{0}\end{array}\right.（2）

核函数种类较多，其中高斯核函数较为常用，表达式为：

K\left(u\right)=\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }}}{\mathrm{e}}^{-\frac{u^{\mathrm{2}}}{\mathrm{2}}}（3）

核密度估计的主要参数选取为核函数与带宽，相较于核函数，带宽对分布模型的拟合效果影响更大。带宽的最佳值通常由积分均方误差（Mean integrated squared error，MISE）的最小值确定，表达式为：

\begin{array}{l}\mathrm{M}\mathrm{I}\mathrm{S}\mathrm{E}\left(h\right)=\int E{\left[F\left(x;h\right)-f\left(x\right)\right]}^{\mathrm{2}}\mathrm{d}x=\\ \frac{h^{\mathrm{4}}}{\mathrm{4}}\int {\left[f^{″}\left(x\right)\int u^{\mathrm{2}}K\left(u\right)\mathrm{d}u\right]}^{\mathrm{2}}\mathrm{d}x+\\ \frac{\mathrm{1}}{nh}\int K^{\mathrm{2}}\left(u\right)\mathrm{d}u+o\left(h^{\mathrm{4}}\right)+o\left(\frac{\mathrm{1}}{nh}\right)\end{array}（4）

式中：E（·）为求均值的函数。

忽略式（4）中的高阶无穷小，可得到渐进积分均方误差（Asymptotic mean integrated squared error，AMISE），表达式为：

\mathrm{A}\mathrm{M}\mathrm{I}\mathrm{S}\mathrm{E}\left(h\right)=\frac{\mathrm{1}}{nh}\int K^{\mathrm{2}}\left(u\right)\mathrm{d}u+\frac{h^{\mathrm{4}}}{\mathrm{4}}\int {\left[f^{″}\left(x\right)\int u^{\mathrm{2}}K\left(u\right)\mathrm{d}u\right]}^{\mathrm{2}}\mathrm{d}x（5）

令式（5）中带宽h的导数为零，可得：

h={\left\{\frac{\int {\left[K\left(u\right)\right]}^{\mathrm{2}}\mathrm{d}u}{\int {\left[f^{″}\left(u\right)\right]}^{\mathrm{2}}\mathrm{d}u\times {\left[\int u^{\mathrm{2}}K\left(u\right)\mathrm{d}u\right]}^{\mathrm{2}}}\right\}}^{\mathrm{0.2}}\times n^{-\mathrm{0.2}}（6）

若选用高斯核函数，则由经验公式得到固定带宽为：

h=\mathrm{1.06}\sigma n^{-\mathrm{0.2}}（7）

式中：\sigma为样本数据的标准差。

若考虑样本数据分布的稀疏程度，基于固定带宽，计算自适应因子λ_i，可得到自适应带宽及其核密度估计。自适应因子λ_i 与自适应带宽核密度估计计算公式分别为：

{\lambda }_i={\left[\frac{F\left(x\right)}{{\left[{\displaystyle \prod _{i=\mathrm{1}}^n}F\left(x\right)\right]}^{-n}}\right]}^{-\alpha }（8）

F\left(x;h\right)=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}\frac{\mathrm{1}}{{\lambda }_{i}h}K\left(\frac{x-x_i}{{\lambda }_{i}h}\right)（9）

式中：\alpha为敏感性参数，其值取0.5时拟合效果较好^［20］。

1.2　扩散核密度估计

由于核密度估计未考虑样本数据的范围，在进行分布估计时会产生边界偏差^［21］。针对该问题，Botev等^［22］在传统核密度估计的基础上提出扩散核密度估计。该方法通过设置初始边界条件，使得边界外无概率分布，在保证分布估计值准确性的同时解决边界偏差问题。

联立式（1）与式（3）得到的高斯核密度估计为偏微分方程的唯一解，表达式为：

\frac{\partial }{\partial h}F\left(x;h\right)=\frac{\mathrm{1}}{\mathrm{2}}\times \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}F\left(x;h\right),x\in R,h>\mathrm{0}（10）

初始边界条件的定义为：

\left\{\begin{array}{l}F\left(x;\mathrm{0}\right)=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}\delta \left(x-x_i\right)\\ \frac{\partial }{\partial x}F\left(x;h\right)\left|{}_{x=x_i}\right.=\frac{\partial }{\partial x}F\left(x;h\right)\left|{}_{x=x_n}\right.=\mathrm{0}\end{array}\right.（11）

式中：δ（x-x_i ）为x_i 处的狄克拉函数。

联立式（10）与式（11），得到扩散核密度估计量为：

F\left(x;h\right)=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}k\left(x,x_i;h\right),x\in \left[\mathrm{0,1}\right]（12）

k\left(x,x_i;h\right)={\displaystyle \sum _{k=-\mathrm{\infty }}^{+\mathrm{\infty }}}K\left(x,\mathrm{2}k+x_i;h\right)+K\left(x,\mathrm{2}k-x_i;h\right)（13）

当带宽h趋近于无穷时，式（13）可表示为^［23］：

k\left(x,x_i;h\right)={\displaystyle \sum _{k=-\mathrm{\infty }}^{+\mathrm{\infty }}}{\mathrm{e}}^{-\frac{k^{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{2}}h}{\mathrm{2}}}\mathrm{c}\mathrm{o}\mathrm{s}\left(k\mathrm{\pi }x\right)\mathrm{c}\mathrm{o}\mathrm{s}\left(k\mathrm{\pi }{x}_i\right)（14）

扩散核密度估计量可表示为：

F\left(x;h\right)\approx {\displaystyle \sum _{k=\mathrm{0}}^{n-\mathrm{1}}}{a}_k{\mathrm{e}}^{-\frac{k^{\mathrm{2}}{\mathrm{\pi }}^{\mathrm{2}}h}{\mathrm{2}}}\mathrm{c}\mathrm{o}\mathrm{s}\left(k\mathrm{\pi }x\right)（15）

a_k=\left\{\begin{array}{l}\mathrm{1},k=\mathrm{0}\\ \frac{\mathrm{2}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}\mathrm{c}\mathrm{o}\mathrm{s}\left(k\mathrm{\pi }{x}_i\right),k=\mathrm{1,2},\cdots ,n-\mathrm{1}\end{array}\right.（16）

2 基于局部自适应带宽扩散核密度估计的载荷外推

基于局部自适应带宽扩散核密度估计的载荷外推主要包括3个部分，分别为二维雨流矩阵数据降维、基于局部自适应带宽扩散核密度估计以及目标频次的载荷外推。其中，基于局部自适应带宽扩散核密度估计的载荷外推流程如图1所示。

2.1　二维雨流矩阵数据的降维

雨流计数法是对幅值与均值进行计数的双参数统计法，其计数原则与材料的应力应变迟滞环相同，能够反映零部件承受载荷本身的变化特征^［24］。考虑到均值的存在，二维均幅值谱无法应用于零件级的台架试验，需将二维谱转化为一维谱，因此本文方法在载荷外推前，需通过平均应力修正进行数据降维。

目前，广泛应用的平均应力修正模型有Gerber模型、Goodman模型、Soderberg模型与Smith模型。相较其他3种模型，Goodman模型更安全，且求解方法更简单，更易实现。因此，选用该模型进行数据降维，其表达式为：

T_{\mathrm{e}}=\frac{T_{\mathrm{A}}}{\mathrm{1}-T_{\mathrm{M}}/T_{\mathrm{b}}}（17）

式中：T_e为修正后的等效幅值；T_A为载荷幅值；T_M为载荷均值；T_b为扭转强度极限。

2.2　基于局部自适应带宽的扩散核密度估计

为了提高数据处理的速率与精度，将载荷数据映射至［0，1］，采用线性归一化方法进行数据标准化处理，表达式为：

x_i^{\mathrm{s}}=\frac{x_i-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}{x_{\mathrm{m}\mathrm{a}\mathrm{x}}-x_{\mathrm{m}\mathrm{i}\mathrm{n}}}（18）

式中：x_i^{\mathrm{s}}为标准化后的数据；x_max为数据x的最大值；x_min为数据x的最小值。

计算数据间的最小差值与最大差值，并根据分块数N进行均分，分别作为带宽h与间隔长度H，并选用高斯核函数。采用不同带宽计算整个样本间隔的估计值，表达式为：

F\left(v;h\right)=\int x_{\left(v~u_n\right)}^{\mathrm{s}}{k}_h\left(u_n\right)\mathrm{d}{u}_n（19）

式中：v为估计点；v~u_n 为估计点v相对于块内数据点数量u_n 的状态。

为了获取最佳估计值，通过间隔长度H调整核函数的形状，从而优化各局部间隔的固定带宽。为得到可优化局部间隔的间隔长度H，引入局部积分均方误差（Local mean integrated squared error，LMISE），其表达式为：

\mathrm{L}\mathrm{M}\mathrm{I}\mathrm{S}\mathrm{E}\left(h;H\right)=\int E{\left[F\left(v;h\right)-f\left(v\right)\right]}^{\mathrm{2}}\rho \left(v\sim x^{\mathrm{s}};H\right)\mathrm{d}v（20）

式中：v~x^s为估计点v相对于x^s的状态；ρ（v~x^s；H）为使数据点x^s的区间长度H的均方差集中的权函数，本文采用高斯权函数进行表达，表达式为：

\rho \left(v\sim x^{\mathrm{s}};H\right)=\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}\mathrm{\pi }}H}{\mathrm{e}}^{-\frac{{\left(v-x^{\mathrm{s}}\right)}^{\mathrm{2}}}{\mathrm{2}{H}^{\mathrm{2}}}}（21）

通过剔除固定带宽h的无关项，得到数据点x^s处的局部损失函数，表达式为：

C_N^{x^{\mathrm{s}}}\left(h;H\right)=\mathrm{L}\mathrm{M}\mathrm{I}\mathrm{S}\mathrm{E}\left(h;H\right)-\int f^{\mathrm{2}}\left(v\right)\rho \left(v~x^{\mathrm{s}};H\right)\mathrm{d}v（22）

为了实现带宽自适应调整，通过LMISE选择各数据点处的带宽h与间隔长度H，获取局部估计损失函数最小时的带宽h_c与间隔长度H_c。局部估计损失函数为：

\begin{array}{l}{\widehat{C}}_N^{x^{\mathrm{s}}}\left(h;H\right)=\int \left[k\left(v-x_i^{\mathrm{s}};h\right)\times k\left(v-x_j^{\mathrm{s}};h\right)\times \right.\left.\rho \left(v\sim x^{\mathrm{s}};H\right)\right]\times \mathrm{d}v-\frac{\mathrm{2}}{n^{\mathrm{2}}}\sum _{i\ne j}{k}_h\left(x_i^{\mathrm{s}}-x_j^{\mathrm{s}}\right)\times \\ \rho \left(v\sim x^{\mathrm{s}};H\right),i\mathrm{、}j=\mathrm{1,2},\cdots ,N\end{array}（23）

局部估计损失函数取决于带宽h与间隔长度H，间隔长度H改变时，最佳带宽h^* 也会改变，因此引入局部带宽优化参数λ_L表示两者间的关系，表达式为：

{\lambda }_{\mathrm{L}}=\frac{h^{\mathrm{*}}}{H}（24）

局部带宽优化参数λ_L调节用于局部优化的间隔长度H：在λ_L较小的情况下，即λ_L趋近于0时，固定带宽h在长间隔内得到优化；在λ_L较大的情况下，即λ_L趋近于1时，固定带宽h同样在长间隔内得到优化。对于各数据点x，可重复获取局部优化带宽h_λ，表达式为^［25］：

h_{\lambda }\left(x^{\mathrm{s}}\right)=\frac{\int {\rho }_{\lambda }\left[x^{\mathrm{s}}\sim u_n;H_{\lambda }\left(u_n\right)\right]{\overline{h}}_{\lambda }\left(x^{\mathrm{s}}\right)\mathrm{d}{u}_n}{\int {\rho }_{\lambda }\left[x^{\mathrm{s}}\sim u_n;H_{\lambda }\left(u_n\right)\right]\mathrm{d}{u}_n}（25）

式中：H_λ （u_n ）为局部带宽优化参数λ_L调节下的间隔长度；ρ_λ ［x^s~u_n；H_λ （u_n ）］为局部带宽优化参数λ_L调节下的权函数；{\overline{h}}_{{\lambda }_L}\left(x^{\mathrm{s}}\right)为重复得到局部优化带宽的均值。

对所得自适应带宽进行逆归一化处理，得到符合原始样本数据的自适应带宽，将其代入式（15）进行扩散核密度估计，即可获得样本数据的概率分布模型。

2.3　目标频次的载荷外推

目前，用于目标寿命周期载荷谱编制的载荷外推常采用MC算法。MC算法作为一种统计试验方法，基于统计抽样试验获取近似解。

基于所构建概率分布模型中变量的估计值，可从已知的母体载荷数据中进行目标次数的抽样，从而达到目标频次的载荷外推。概率分布模型可视为相关随机变量的概率分布，因此在实际操作中，首先生成一组服从均匀分布的随机数，再由概率分布模型将其转换为服从该模型的随机数。将这些随机数视为样本载荷数据的索引，经目标次数筛选，最终得到满足目标频次要求的载荷数据。

3 实例验证

本文通过采集某特种车辆综合传动装置的载荷数据，对所提方法进行分析及验证。

3.1　载荷数据采集及预处理

试验对象为某型履带车辆的综合传动装置，为测量其输出转矩，在该装置汇流行星排的输出轴上安装转矩传感器。转矩传感器放置位置示意图如图2所示。为了获取车辆的行驶状态，同步采集主要参数，如表1所示。

在实际测试过程中，由于测试系统、测试环境和人为因素的干扰，获取的综合传动装置的输出转矩数据经常存在杂波与异值点等问题，在获取实测载荷数据后，需进行滤波降噪与奇异点修正等预处理。

本文选用小波变换方法对实测载荷进行滤波降噪。经综合对比不同参数下的处理结果，最终选择Sym4小波基函数、10层分解层数、贝叶斯估计降噪方法、中值法阈值规则，同时假定噪声不随尺度变化，以此实现最优降噪效果。

在数据预处理环节，采用基于工况分割的奇异点修正方法。通过工况划分获取滤波降噪后各工况的载荷数据，根据不同工况的载荷规律，利用箱线图法对各工况的异常值进行识别并修正。预处理后的汇流排输出转矩数据如图3所示。

3.2　结果对比

本文采用四点雨流计数法对预处理后的汇流排输出转矩数据进行载荷循环统计处理，并合并统计的各工况雨流矩阵，得到合成工况的雨流矩阵，如图4所示。

通过式（17），将合成工况的二维均幅值雨流矩阵进行降维，得到一维等效幅值。为验证本文提出的基于局部自适应带宽扩散核密度估计方法的准确性，将其与以下方法进行对比分析：基于式（7）固定带宽扩散核密度估计的载荷外推法（Diffusion kernel density estimation-fixed bandwidth-experience，DKDE-FB-EXP）、基于改进Sheather-Jones法（Improved Sheather-Jones，ISJ）^［22］固定带宽扩散核密度估计的载荷外推法（Diffusion kernel density estimation-fixed bandwidth-improved Sheather-Jones，DKDE-FB-ISJ）、结合上述2种固定带宽与式（8）的自适应带宽扩散核密度估计载荷外推法（（Diffusion kernel density estimation-adaptive bandwidth-experience，DKDE-AB-EXP）与（Diffusion kernel density estimation-adaptive bandwidth-improved Sheather-Jones，DKDE-AB-ISJ））、基于局部自适应带宽核密度估计的载荷外推法（Kernel density estimation-local adaptive bandwidth，KDE-LAB），并以等效幅值数据作为参照。所有方法均选择高斯核函数，见式（3）。为提高计算效率，设定128个估计点。对于KDE-LAB方法与本文提出的方法，通过归一化等效幅值的最大差值、最小差值及分块数，设置相等的初始带宽h₀与初始间隔长度H₀，表达式为：

\begin{array}{l}{B}_i=\mathrm{\Delta }{x}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{s}}+\frac{\left(\mathrm{\Delta }{x}_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{s}}-\mathrm{\Delta }{x}_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\mathrm{s}}\right)\left(i-\mathrm{1}\right)}{N},\\ i=\mathrm{1,2},\mathrm{3},\cdots ,N\end{array}（26）

式中：B_i 为相等的初始带宽h₀与初始间隔长度H₀；∆x^s_min为标准化等效幅值x^s的最小差值；∆x^s_max为标准化等效幅值x^s的最大差值。

等效幅值与不同载荷外推方法的密度分布如图5所示。

由图5可知，DKDE-FB-EXP方法与DKDE-AB-EXP方法得到的密度分布曲线重合，DKDE-FB-ISJ方法与DKDE-AB-ISJ方法得到的密度分布曲线同样重合，表明2种带宽求取方法的固定带宽与自适应带宽估计结果一致。由于自适应带宽为固定带宽与自适应系数的乘积，因此估计结果相同时，采用固定带宽较为便捷。与其他外推方法相比，本文所提方法更贴近等效幅值的密度分布情况。为了进一步对比密度分布拟合效果，计算等效幅值与上述外推方法间的相关系数R、决定系数R²及均方根误差（Root mean square error，RMSE）。其中，相关系数R与决定系数R²越接近1、RMSE越小，说明拟合效果越好。等效幅值分布的拟合结果如表2所示。

结合图5与表2可知，DKDE-FB-ISJ方法与DKDE-AB-ISJ方法得到的估计结果与等效幅值的相关系数R和决定系数R²均高于DKDE-FB-EXP方法与DKDE-AB-EXP方法，且RMSE更小，表明基于ISJ方法的固定带宽扩散核密度估计与自适应带宽扩散核密度估计的效果较传统经验公式更优；KDE-LAB方法得到的相关系数R与决定系数R²也高于传统扩散核密度估计方法，且RMSE更小，说明基于局部自适应带宽的核密度估计方法处理效果更佳。由本文方法得到的相关系数R与决定系数R²更接于1，分别为0.983 8与0.967 9，RMSE仅为5.05×10^－5。这是由于在局部自适应带宽的核密度估计基础上，进一步引用扩散核密度估计，消除了范围外的概率密度分布，使得本文方法所得等效幅值分布与实际等效载荷分布具有更高的相似度，因此分布拟合效果更佳，更贴近等效载荷的实际分布规律。

载荷循环的累积频次曲线是编制载荷谱与试验谱的数据基础，决定着编制的准确性。因此，本文将等效幅值数据与上述多种载荷外推方法（包括本文方法）所得累积频次曲线进行对比验证。由于MC算法为随机采样，本文选择多次外推结果的平均值作为对比数据。等效幅值与不同载荷外推方法的累计频次曲线如图6所示。

由图6可知，通过本文方法外推得到的等效幅值累计频次曲线与实际等效幅值累计频次曲线更为接近，拟合效果更好。为了进一步评估各方法对等效幅值累积频次曲线的拟合效果（见表3），沿用上述相同的检验指标进行对比。

由表3可知，本文方法得到的相关系数R与决定系数R²更接近1，分别为0.999 6与0.999 1，且RMSE仅为15.9。这表明与传统载荷外推方法相比，本文方法所得等效幅值累积频次曲线与实际等效幅值累积频次曲线间的拟合度与相似度更高，效果更佳，可为载荷谱与试验谱编制提供更准确的数据基础。

4 结 论

（1）本文提出一种基于局部自适应带宽扩散核密度的载荷外推估计方法。该方法首先对二维均幅值雨流矩阵进行降维，通过带宽与间隔长度计算损伤函数，并利用局部最小损失函数选取带宽优化参数；然后使用最优带宽对一维等效幅值进行扩散核密度估计，构建概率密度分布模型；最后结合MC算法，最终实现目标频次的载荷外推。

（2）以某特种车辆综合传动装置为研究对象，采集某型履带车辆综合传动装置的汇流排输出转矩等数据。在数据预处理阶段，采用小波变换方法进行滤波降噪，并基于工况分割的奇异点修正策略，有效提高数据质量。

（3）将多种传统载荷外推方法与本文方法进行对比验证，结果表明：本文方法所得概率密度分布模型及外推累计频次曲线与实际等效幅值的相关系数分别为0.983 8与0.999 6，决定系数分别为0.967 9与0.999 1，均更接近1，且RMSE仅为5.05×10^－5与15.9，拟合效果更佳，可为后续载荷谱与试验谱编制、零件级台架试验等提供理论参考。

参考文献

[1] 徐保荣, 王涛, 梁梓. 基于操纵动作预测的履带车辆载荷谱编制方法与流程[J]. 兵工学报, 2022, 43(2): 252-259.

[2] Xu Bao-rong, Wang Tao, Liang Zi. Generating method and process of load spectra for tracked vehicles based on prediction of manipulating action[J]. Acta Armamentarii, 2022, 43(2): 252-259.

[3] 骆清国, 王旭东, 张更云, 等. 装甲车辆发动机载荷谱仿真与试验研究[J]. 兵工学报, 2012, 33(12): 1442-1447.

[4] Luo Qing-guo, Wang Xu-dong, Zhang Geng-yun, et al. Loading spectrum simulation and experiment study of the armored vehicle engine[J]. Acta Armamentarii, 2012, 33(12): 1442-1447.

[5] 牛文铁, 才福友, 付景静. 基于自适应带宽核密度估计的载荷外推方法研究[J]. 农业机械学报, 2021, 52(1): 375-384.

[6] Niu Wen-tie, Cai Fu-you, Fu Jing-jing. Load extrapolation method based on adaptive bandwidth kernel density estimation[J].Transactions of the Chinese Society for Agricultural Machinery, 2021, 52(1): 375-384.

[7] 于立娟, 安阳, 何佳龙, 等.机电装备载荷谱外推技术研究进展及发展趋势[J]. 吉林大学学报: 工学版, 2023, 53(4): 941-953.

[8] Yu Li-juan, An Yang, He Jia-long, et al. Research progress and development trend of extrapolation method in electromechanical equipment load spectrum[J].Journal of Jilin University (Engineering and Technology Edition), 2021, 52(1): 375-384.

[9] 张文胜. 某履带车辆传动系统的传动轴载荷谱编制方法研究[D]. 北京: 北京理工大学机械与车辆学院, 2016.

[10] Zhang Wen-sheng. Research on load spectrum compiling for transmission shaft in transmission system of tracked vehicle[D]. Beijing: School of Mechanical Engineering, Beijing Institute of Technology, 2016.

[11] 田维波, 宋绪丁, 万一品. 常用载荷谱外推算法概述[J]. 南方农机, 2019, 50(3): 103, 108.

[12] Tian Wei-bo, Song Xu-ding, Wan Yi-ping. Overview of commonly used load spectrum extrapolation algorithms[J].South Agricultural Machinery, 2019, 50(3): 103, 108.

[13] 刘岩, 张喜逢, 王振雨, 等. 载荷谱外推方法的对比[J]. 现代制造工程, 2011(11): 8-11.

[14] Liu Yan, Zhang Xi-feng, Wang Zhen-yu, et al. Contrast of extrapolations in compiling load spectrum[J]. Modern Manufacturing Engineering, 2011(11): 8-11.

[15] Chen C H, Yang Z J, He J L, et al. Load spectrum generation of machining center based on rainflow counting method[J]. Journal of Vibroengineering, 2017, 19(8): 5767-5779.

[16] He J L, Wang S X, Li G F,et al. Compilation of NC lathe dynamic cutting force spectrum based on two-dimensional mixture models[J]. The International Journal of Advanced Manufacturing Technology, 2018, 98(1/4): 251-262.

[17] Zhang J, Lu J W, Han W D, et al. Program load spectrum compilation for accelerated life test of parabolic leaf spring[J]. International Journal of Automotive Technology, 2019, 20(2): 337-347.

[18] 于佳伟, 郑松林, 赵礼辉, 等. 基于非参数核密度估计法的车辆大数据服役载荷外推方法[J]. 中国机械工程, 2021, 32(11): 1307-1314.

[19] Yu Jia-wei, Zheng Song-lin, Zhao Li-hui, et al. Extrapolation method of vehicle big-data service-loads based on non-parametric KDE method[J]. China Mechanical Engineering, 2021, 32(11): 1307-1314.

[20] Wang J X, Liu Y, Zeng X H, et al. Selection method for kernel function in nonparametric extrapolation based on multicriteria decision-making technology[J]. Mathematical Problems in Engineering, 2013, 2013(8): 1-11.

[21] 高天宇, 李焕良, 郑铮, 等. 基于雨流频次外推的ZLK50型装载机液压缸载荷谱编制[J]. 机械强度, 2017, 39(4): 951-956.

[22] Gao Tian-yu, Li Huan-liang, Zheng Zheng, et al.Compiling method of load spectrum of ZLK50 type loader hydraulic cylinder based on rain flow frequency extrapolation[J]. Journal of Mechanical Strength, 2017, 39(4): 951-956.

[23] Wang S X, Li G F, Yang H J, et al.Research on extreme value distribution of cutting force in CNC lathe based on generalized extreme value distribution[C]∥2019 4th International Conference on System Reliability and Safety(ICSRS). Piscataway, N J: IEEE, 2019:198-202.

[24] Fogle J, Agarwal P, Manuel L. Towards an improved understanding of statistical extrapolation for wind turbine extreme loads[J]. Wind Energy, 2008, 11(6): 613-635.

[25] Yuan K, Chen X Q, Gui Z P, et al.A quad-tree-based fast and adaptive kernel density estimation algorithm for heat-map generation[J]. International Journal of Geographical Information Science, 2019, 33(12): 2455-2476.

[26] Rosenblatt M. Remarks on some nonparametric estimates of a density function[J]. The Annals of Mathematical Statistics, 1956, 27(3): 832-837.

[27] Parzen E. On estimation of probability density function and mode[J]. The Annals of Mathematical Statistics, 1962, 33(3): 1065-1076.

[28] 李研. 轮式装载机载荷谱编制中非参数外推方法研究[D]. 长春: 吉林大学机械科学与工程学院, 2016.

[29] Li Yan. Research on non-parametric extrapolation method in compiling load spectra of wheel loader[D]. Changchun:School of Mechanical Science and Engineering, Jilin University, 2016.

[30] Abramson I S. On bandwidth variation in kernel estimates-a square root law[J]. The Annals of Statistics,1982, 10(4): 1217-1223.

[31] 张龙升. 基于扩散核密度估计的车辆后副车架疲劳损伤分析[D]. 秦皇岛: 燕山大学车辆与能源学院, 2021.

[32] Zhang Long-sheng. Fatigue damage analysis of automobile rear subframe based on kernel density estimation via diffusion[D]. Qinhuangdao:School of Vehicles and Energy, Yanshan University, 2021.

[33] Botev Z I, Grotowski J F, Kroese D P.Kernel density estimation via diffusion[J].The Annals of Stats,2010, 38(5): 2916-2957.

[34] Xu X Y, Yan Z, Xu S L. Estimating wind speed probability distribution by diffusion based kernel density method[J]. Electric Power Systems Research,2015, 121: 28-37.

[35] 阎楚良, 王公权. 雨流计数法及其统计处理程序研究[J]. 农业机械学报, 1982(4): 90-103.

[36] Yan Chu-liang, Wang Gong-quan. Investigation of rain flow counting method and its program of statistical treatment[J]. Transactions of the Chinese Society for Agriculture Machinery, 1982(4): 90-103.

[37] Smith A C, Brown E N. Estimating a state-space model from point process observations[J]. Neural Computation, 2003, 15(5): 965-991.

基金资助

国家自然科学基金项目(52175074)