以常见的某型液压挖掘机为研究对象, 进行土石挖掘试验。挖掘机的主要操作是挖掘循环和行走, 为简化分析, 将履带位置固定, 挖掘机仅重复挖掘作业。具体操作模式是将油门固定于最大开度, 进行全挖掘深度、全挖掘半径、全速回转并90° 卸料。利用数据采集仪获得零部件的载荷-时间历程信息, 以15个挖掘循环进行一次数据统计。试验现场、样机及相关细节见图1。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 挖掘机挖掘试验Fig.1 Excavation operation test of excavator

本文主要针对主泵出口压力载荷, 提取其挖掘段载荷-时间历程信息, 开展载荷谱编制的相关研究。首先提取试验测得的数据转折点, 然后进行雨流计数, 得到均值和幅值载荷分布直方图^[11]。根据载荷统计规律可知, 工程机械的载荷谱均值常服从正态分布, 幅值服从威布尔分布。图2为均、幅值载荷柱状图及其对应的拟合检验P-P图。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 均、幅值载荷柱状图及P-P图Fig.2 Mean and amplitude value histogram and P-P plot

由图2可知, 载荷的均、幅值数据分布在代表理论分布的对角线附近, 很难直接判断其拟合效果。为得到定量性结论, 采用卡方检验判断多种常用模型的拟合效果, 如表1所示。

表1

Table 1

表1(Table 1)

表1 卡方拟合检验结果 Table 1 Fitting results of Chi-square test 卡方检验 p 检验值 均值 正态 0.0015 19.54 威布尔 1.94× 10-11 89.29 对数正态 3.85× 10-10 79.85 幅值 正态 4.0169× 10-35 172.31 威布尔 1.4677× 10-102 386.66 对数正态 7.3319× 10-26 107.24

表1 卡方拟合检验结果 Table 1 Fitting results of Chi-square test

通常采用p=0.05作为检验是否通过的临界值, 观察表1, 不论何种分布, 二者均未通过卡方检验, 且幅值的拟合效果更差。因此, 基于篇幅的考虑, 本文只对效果更差的幅值载荷进行进一步分析。针对幅值载荷, 从表1可以知道, 上述3种常用分布都不能够满足幅值载荷分布函数的拟合要求, 因此, 可考虑采用混合分布函数对幅值载荷进行拟合。

当前工程机械载荷分布常用的几种函数分别是正态分布、对数正态分布以及威布尔分布, 概率密度函数分别如下:

式中: f1(x)为正态分布的概率密度函数; μ1为正态分布均值; σ1为正态分布的标准差; x为每级载荷的中值。

式中: f2(x)为对数正态分布的概率密度函数; μ2为对数正态分布均值; σ2为对数正态分布的标准差。

f3(x)=βαxαβ-1e-(x/α)β(3)

式中: f3(x)为威布尔分布的概率密度函数; α为威布尔分布尺度参数; β为威布尔分布的形状参数。

进行数据拟合时, 卡方检验是最常用的检验准则, 主要用于对比观察值与拟合值的偏离程度, 见式(4):

χ2=∑i=1kpi-N[F(xi+1)-F(xi)]2N[F(xi+1)-F(xi)](4)

式中: χ2为载荷分布拟合检验的总卡方值; k为组数; i为级数; pi为观察值的频数; N为载荷的总频数; F为累积分布函数; xi+1为第 i级载荷组的上限值; xi为第 i级载荷组的下限值。

在进行卡方拟合检验时, 采用Luise分组法^[9]来确定载荷直方图的组距, 见式(5):

Δ=R1+3.3lgN(5)

式中: R为载荷变化范围。

为了提高分布估计时的精度以及满足卡方检验的需求, 适当调整 Δ值, 以保证每组的最小频次不低于5。

利用分段函数组成的混合分布来描述载荷谱的分布情况, 其主要内容是将若干个不同参数的概率密度函数以截取的方式结合在一起组成新的概率密度函数。组成该分布函数的关键在于分段函数分段数的确定, 同时也是未知参数估计的前提。对于任一给定的数学模型, 有多种不同的确定函数个数的方法, 文献[10]已经对这些方法进行了阐述。其中AIC准则算法相较于其他方法更加稳定, 其结果也更加接近最优解。该准则是基于惩罚函数的思想提出的, 当纠正对数似然偏差降到最小时, 对应的个数就认为是函数的最优解, 即为本文所需模型的基本函数的个数。

AIC=2m-2lnL(m)(6)

式中: m为分段数; lnL(m)= ∑i=1nIn[f(x_i, θ )]是关于 m的极大对数似然函数值, f(xi, θ)是含有 m的概率密度函数。

K-means算法是一种广泛应用的聚类算法, 采用划分的方法将数据集中的对象按照某一种方式聚类到合适的簇中, 保证每个对象只属于一个簇。本文选取数据集中所有数据对象的误差平方和最小作为聚类依据。但是, 簇的数目是需要决策者按照以往经验进行指定。而AIC准则恰好解决了簇的个数问题。根据AIC准则得到段数为 m,载荷聚类为 C1、C2、…、Cm。对每一簇的数据分别采用MLE法, 得到分段混合分布函数的表达式如下:

f(xm, θ)=λ1f1(x, θ1), x⊂C1λifi(x, θi), x⊂Ciλmfm(x, θm), x⊂Cm(7)

式中: x为载荷变量; f(xm, θ)为分段混合分布函数的概率密度函数; fi(x, θi)为各段对应的概率密度函数; λi为各分段函数的比例系数; θ为各分段函数的参数集合。

首先对处理后的幅值数据进行卡方检验, 选取拟合的分布函数分别是正态分布、对数正态分布和威布尔分布。为保证每个组内具有充分多的数据, 此处选取3.35 MPa作为幅值载荷直方图的组距, 得到如图4所示的一维幅值谱。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 一维幅值谱Fig.4 One dimensional amplitude spectrum

观察表1, 结合临界卡方检验值9.49可知, 该幅值载荷不满足任何一种单一分布的拟合要求, 同时可知对数正态分布拟合下的卡方检验值最小。因此, 选取对数正态分布作为各分段函数的拟合函数。段数个数分别为1, 2, …, 10时, AIC的计算结果见图5。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 分段数与AIC值Fig.5 Components number and AIC value

因此当 m=5时, 对应的AIC值最小。采用K-means法对数据进行聚类, 并将聚类后的各组数据分别进行对数正态拟合, 得到各组估计参数见表2。将表2中的参数结果代入式(7)中, 取显著性水平 α=0.05进行卡方检验, 得到对应的卡方检验值为6.5216( p=0.3674), 因此利用本文方法可以得到拟合效果较好的分布函数。

表2

Table 2

表2(Table 2)

表2 各组数据对数正态分布拟合参数估计值 Table 2 Estimated parameter values of each group by logarithmic normal distribution λ μ2 σ2 1 0.5685 1.6221 0.2279 2 0.1644 2.4127 0.1779 3 0.1130 2.9049 0.1135 4 0.0822 3.2352 0.0664 5 0.0719 3.4808 0.0823

表2 各组数据对数正态分布拟合参数估计值 Table 2 Estimated parameter values of each group by logarithmic normal distribution

为验证本文方法的稳定性, 将本文方法应用到已知分布参数的仿真数据中。已知3组各5000个数据分别服从N₁~(0, 1)(分布1), N₂~(3, 1)(分布2)和N₃~(6, 1)(分布3)。对数据进行截取, 在服从N₁~(0, 1)分布的数据中选取小于1的数据; 服从 N2~(3, 1)分布的数据中选取[1, 5]之内的数据; 而在服从 N3~(6, 1)分布的数据中选取大于5的数据, 重复1000次操作后可知各分布参数为 N'1~(-0.2876, 0.7934), N'2~(3.0002, 0.8796)和 N'3~(6.2872, 0.7936)。利用本文方法, 在数据聚类划分后, 采用极大似然估计法对参数进行估计, 结果见图6~图8。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 分布1参数估计值Fig.6 Estimated parameter values of distribution 1

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 分布2参数估计值Fig.7 Estimated parameter values of distribution 2

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 分布3的参数估计值Fig.8 Estimated parameter values of distribution 3

结合图6~图8可知, 对于1000组仿真数据, 利用聚类法得到的3个分布的参数估计值的分布接近正态分布, 采用上述直方图的均值代表各估计值, 与原均值的对比结果见表3。由于仿真数据的生成具有随机性, 截取后的位置参数值会发生偏移, 而分布2两边截取, 偏移程度相抵消。同时, 数据量的变化使形状参数发生改变。因此结合表1可知, 估计值在原值附近发生不同程度的偏移, 且偏移倍数在[0.82, 1.04]内, 造成的误差为0~18%。

表3

Table 3

表3(Table 3)

表3 参数估计结果分析 Table 3 Analysis of parameter estimation values 原均值a 估计均值b b/a |b-a|/a 分布1 位置参数 -0.2876 -0.2355 0.82 0.18 形状参数 0.7934 0.8269 1.04 0.04 分布2 位置参数 3.0002 3.0007 1.00 0.00 形状参数 0.8796 0.7809 0.89 0.11 分布3 位置参数 6.2872 6.2359 0.99 0.01 形状参数 0.7936 0.8265 1.04 0.04

表3 参数估计结果分析 Table 3 Analysis of parameter estimation values

随机选取一组仿真数据分析观察数据拟合效果, 数据总样本量为131 908, 取值范围为[-4.2486, 10.2170]。3个分布选取样本量分别为42 188、47 739以及41 981, 取值范围分别为[-4.2486, 0.9998]、[1.0001, 4.9997]、[5.0002, 10.2170]。

对数据进行聚类后得到3组数据的取值分别为[-4.2486, 1.3848]、[1.3857, 4.6214]、[4.6220, 10.2170], 样本量分别为43 678、44 767以及43 463。通过MLE得到正态分布的参数估计值见表4。

表4

Table 4

表4(Table 4)

表4 各组数据正态分布拟合参数估计值 Table 4 Estimated parameter values of each group by normal distribution λ μ σ 1 0.3311 -0.2344 0.8242 2 0.3394 3.0045 0.7817 3 0.3295 6.2394 0.8270

表4 各组数据正态分布拟合参数估计值 Table 4 Estimated parameter values of each group by normal distribution

根据参数估计结果得到数据拟合效果见图9。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 频数直方图与混合正态分布密度函数Fig.9 Frequency histogram and mixed normal distribution density function

根据文献[9]提出的方法, 本文依旧选取上述3种分布函数进行函数分布估计。单一分布函数的拟合参数值见表5。

表5

Table 6

表5(Table 6)

表5 单一分布函数和分段函数拟合参数估计值 Table 6 Estimated parameter values of each group by single distribution function and piecewise function 函数 正态分布 对数正态分布 威布尔分布 μ1 σ1 μ2 σ2 α β 单一分布 11.33 8.88 2.16 0.70 12.57 1.40 分段 26.59 2.06 1.96 0.54 34.83 16.15

表5 单一分布函数和分段函数拟合参数估计值 Table 6 Estimated parameter values of each group by single distribution function and piecewise function

根据表5的参数值, 可以得到每一级的卡方值, 并将各级的最小卡方值标记出, 绘成图10。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 三种分布函数拟合检验的各级卡方值Fig.10 Chi-square value of each level obtained by three kinds of distribution function fitting test

从图10中可以看出, 对前5级数据, 采用对数正态分布拟合较好, 正态分布拟合第6、7级较好, 威布尔分布对尾部数据(即8、9级数据)拟合较好。将载荷数据分成3部分, 分别进行对数正态、正态以及威布尔分布的参数估计, 进而得到式(7)中的参数估计结果见表5。

根据参数估计值以及式(4)得到分段函数各级卡方值及总卡方值, 见表6。由表6可知, 分段函数的各级卡方检验值及总卡方检验值(ZK)在一定程度上明显小于单一分布时的卡方检验量, 但分段函数的总卡方值大于此分布下的临界卡方检验值(LK), 因此该分段函数没有通过卡方检验。

表6

Table 6

表6(Table 6)

表6 分段函数的卡方拟合检验值 Table 6 Chi-square test values of piecewise function 级数 卡方值 对数正态 1 37.55 2 4.25 3 0.02 4 15.22 5 31.21 正态 6 0.45 7 0.85 威布尔 8 2.13 9 0.02 ZK 91.71 LK 5.99

表6 分段函数的卡方拟合检验值 Table 6 Chi-square test values of piecewise function

为验证本文方法对拟合函数的依赖性, 各分段函数比例系数不变, 分别选取正态分布和威布尔分布作为本文方法的各段拟合函数。分别计算得到的参数估计值见表7。

表7

Table 7

表7(Table 7)

表7 各组数据正态分布和威布尔分布拟合参数 Table 7 Estimated parameter values of each group by normal and Weibull distribution function 正态分布 威布尔分布 μ1 σ1 α β 1 5.20 1.22 5.69 4.49 2 11.34 2.00 12.18 6.41 3 18.38 2.08 19.30 9.98 4 25.46 1.67 26.22 17.92 5 32.59 2.72 33.86 12.74

表7 各组数据正态分布和威布尔分布拟合参数 Table 7 Estimated parameter values of each group by normal and Weibull distribution function

同时, 得到卡方检验值见表8。

表8

Table 8

表8(Table 8)

表8 分段混合正态及混合威布尔分布的卡方检验值 Table 8 Chi-square test values of piecewise mixed normal and Weibull distribution 正态分布 威布尔分布 H 1 0 P 1.8458× 10-5 0.0892 χ2 31.7206 10.9750

表8 分段混合正态及混合威布尔分布的卡方检验值 Table 8 Chi-square test values of piecewise mixed normal and Weibull distribution

利用本文方法, 改变分段混合分布函数的拟合函数类型。根据单一分布时的拟合卡方检验值, 当选取卡方检验值略大于对数正态分布的威布尔分布为各段的拟合函数, 依旧可以得到符合拟合条件的分段混合分布函数。正态分布对应的卡方检验值是对数正态分布卡方值的2倍以上, 当拟合函数为正态分布时, 不能够通过卡方检验。因此可知, 当拟合函数发生改变时, 本文方法也能得到较好的结果。从而证明本文方法具有较广的拟合函数适用范围。