基于修正q-威布尔分布的矿用卡车可靠性分析

Reliability analysis of mining trucks based on new alpha transformed q-Weibull distribution

摘要

为了更加准确地描述露天矿矿用卡车的失效规律，提高可靠性分析的准确性，构建了一种新的alpha变换。在此基础上，提出了一种四参数修正q-威布尔分布模型，并采用蜣螂优化算法与极大似然估计相结合的方式对模型的参数进行估计。通过实例对比验证了使用修正q-威布尔分布模型评估矿用卡车可靠性的合理性和有效性。数值试验结果表明，利用修正q-威布尔分布模型对矿用卡车故障间隔时间进行分析，制定相应的预防性维修周期能够更好地保障矿用卡车安全、稳定运行。

Abstract

In order to more accurately describe the failure pattern of mining trucks in open-pit mines and improve the accuracy of reliability analysis, a new alpha transformation is constructed. On this basis, a four-parameter new alpha transformed q-Weibull distribution model is proposed, and the parameters of the model are estimated using a combination of dung beetle optimization algorithm and maximum likelihood estimation. The rationality and effectiveness of using the new alpha transformed q-Weibull distribution model to evaluate the reliability of mining trucks are verified through comparison of examples. Numerical test results show that using the model to analyze the time between faults of mining trucks and formulating corresponding preventive maintenance cycles can better ensure the safe and stable operation of mining trucks.

Graphical abstract

关键词

矿用卡车 / 可靠性分析 / 修正q-威布尔分布 / 蜣螂优化算法 / 预防性维修周期 / 极大似然估计

Key words

mining trucks / reliability analysis / new alpha transformed q-Weibull distribution / dung beetle optimization algorithm / preventive maintenance cycle / maximum likelihood estimation

0 引言

随着信息技术、人工智能技术的飞速发展，矿产资源开发利用的智能化特征日益显著，例如，煤矿露天开采中，以设备智能化为代表的智慧露天矿建设飞速发展^[1]，矿用卡车作为电铲-卡车间断开采工艺和半连续开采工艺中最为重要的运输设备，在智慧露天矿建设进程中扮演着重要的角色。随着无人驾驶等技术的逐步推广和应用，矿用卡车能够越来越多地“感知”到自身的运行状态。矿用卡车安全、稳定运行对矿山生产至关重要。因此，及时有效地辨识出矿用卡车在充满噪声、浮尘和颠簸的恶劣环境中持续工作而引起的机械性能退化特征，并对机器失效行为和失效时间进行有效预判，实现智能设备的预防性维修，对保障智慧矿山建设有着重要的意义。

围绕矿用卡车的失效规律开展可靠性分析是实现预防性维修的重要基础性工作。在此方面，刘威等^[2]基于马尔科夫蒙特卡罗方法，分析了矿用卡车的失效时间与类别。RAHIMDL^[3]在可靠性框图的基础上，建立了矿用卡车的故障树。HE等^[4]和QARAHASANLOU等^[5]应用比例风险模型分析了矿用卡车的可靠性。概率分布在可靠性分析中占有重要地位，通过数学公式来描述数据规律，能够给出具体的数值指标，具有严密的理论基础，可以对系统可靠性提供准确的量化评估^[6]。

威布尔分布是描述机械系统及零部件寿命数据分布规律最常用的一种分布形式。TORAMAN^[6]通过威布尔分布、正态分布和对数正态分布对矿用卡车的可靠性进行分析，得出故障间隔时间数据多服从威布尔分布。YUAN等^[7]用二维威布尔分布对矿用卡车二维寿命数据进行建模。KUMAR等^[8]和MONIRI-MORAD等^[9]通过威布尔分布对矿用卡车进行可靠性分析，并计算了预防性维修的时间间隔。然而威布尔分布的失效率函数只能是递增、递减或常数^[10]，因此不适用于更复杂的生命周期数据。多年来，研究人员一直在开发威布尔分布的各种修正形式^[11]。

在非广延统计力学的背景下，PICOLI等^[12]将威布尔分布推广为q-威布尔分布（q-Weibull distribution，q-Weibull）。与只能描述单调失效率函数的威布尔分布相比，q-威布尔分布能模拟失效率的各种类型，包括单峰、浴盆型、单调型（单调递减、递增）和常数。此外，q-威布尔分布能够呈现短尾分布和长尾分布，充分拟合具有不同形状和尾部属性的数据集。樊立天等^[13]利用q-威布尔分布对医疗设备进行可靠性分析，验证了q-威布尔分布的有效性。MACHADO DE ASSIS等^[14]利用q-威布尔对水电设备的主要可靠性参数进行了计算，表明q-威布尔分布对故障现象的建模更接近实际情况，可以更好地规划维护和更换策略。MARCHETTI等^[15]利用q-威布尔分布对油气设施安全屏障老化问题进行了研究。然而以往的研究大多是基于原始q-威布尔分布模型进行计算，只可控制分布的尾部，没有针对模型的前端进行改进，并且仍未有研究利用q-威布尔分布对露天矿矿用卡车进行可靠性分析。

为了解决上述问题，提出一种新的alpha变换方法，建立修正q-威布尔分布模型。采用蜣螂优化算法与极大似然估计相结合对模型的参数进行估计。在矿用卡车故障数据集上对其性能进行仿真研究，验证该模型的有效性，并给出矿用卡车的预防性维护建议。

1 修正q-威布尔分布模型

1.1 q-威布尔分布模型

非广延统计力学在数学中引入了指数和对数函数的推广^[16]，即\mathrm{e}\mathrm{x}{\mathrm{p}}_q\left(x\right)和\mathrm{l}{\mathrm{n}}_q\left(x\right)，分别称为q-指数函数和q-对数函数。这两个函数都有一个控制参数q，函数定义分别为

\mathrm{e}\mathrm{x}{\mathrm{p}}_q\left(x\right)=\left\{\begin{array}{l}{\left[\mathrm{1}+\left(\mathrm{1}-q\right)x\right]}^{\frac{\mathrm{1}}{\mathrm{1}-q}} \left[\mathrm{1}+\left(\mathrm{1}-q\right)x\right]>\mathrm{0}\\ \mathrm{0} \mathrm{其他}\end{array}\right.，

\mathrm{l}{\mathrm{n}}_q\left(x\right)=\frac{x^{\mathrm{1}-q}-\mathrm{1}}{\mathrm{1}-q}\left(x>\mathrm{0},q\ne \mathrm{1}\right)，

式中，x,q\in \mathbf{R}。

将q-指数函数扩展到威布尔分布中，得出q-威布尔分布的概率密度函数（probability density function，PDF）为

f_q\left(x\right)=\left(\mathrm{2}-q\right)\lambda \theta x^{\theta -\mathrm{1}}\mathrm{e}\mathrm{x}{\mathrm{p}}_q\left(-\lambda x^{\theta }\right)=

\begin{array}{c}\end{array}\left(\mathrm{2}-q\right)\lambda \theta x^{\theta -\mathrm{1}}{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{1}}{\mathrm{1}-q}}。

通过\left(\mathrm{2}-q\right)确保，使PDF归一化。并且

，

式中，x_{\mathrm{m}\mathrm{a}\mathrm{x}}={\left[\lambda \left(\mathrm{1}-q\right)\right]}^{-\mathrm{1}/\theta }，以保证PDF在时的概率性质。当q\to \mathrm{1}时，f_q\left(x\right)退化为威布尔分布的概率密度函数。

q-威布尔分布的累积分布函数（cumulative distribution function，CDF）、可靠度函数（reliability function）和失效率函数（hazard rate function）分别为

F_q\left(x\right)=\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}，

R_q\left(x\right)=\mathrm{1}-F_q\left(t\right)={\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}，

h_q\left(x\right)=\frac{f_q\left(x\right)}{R_q\left(x\right)}=\frac{\left(\mathrm{2}-q\right)\lambda \theta x^{\theta -\mathrm{1}}}{\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }}。

1.2 修正q-威布尔分布模型

由于q-威布尔分布的前端缺乏灵活性，本文提出一种新的alpha变换（new alpha transformation，NAT）方法来修正q-威布尔分布。

设F\left(x\right)为连续随机变量\mathit{X}的累积分布函数，\overline{F}\left(x\right)为连续随机变量\mathit{X}的可靠度函数，引入参数\alpha，将累积分布函数变换为幂次形式，增加灵活性。NAT变换扩展分布的CDF为

G\left(x;\alpha \right)=\frac{F^{\alpha }\left(x\right)}{\overline{F}\left(x\right)+F^{\alpha }\left(x\right)}=\mathrm{1}-\frac{\overline{F}\left(x\right)}{\overline{F}\left(x\right)+F^{\alpha }\left(x\right)}。

对CDF进行求导，得到变量\mathit{X}的概率密度函数为

g\left(x;\alpha \right)=\frac{F^{\alpha -\mathrm{1}}\left(x\right)f\left(x\right)\left(\alpha \overline{F}\left(x\right)+F\left(x\right)\right)}{{\left(\overline{F}\left(x\right)+F^{\alpha }\left(x\right)\right)}^{\mathrm{2}}}。

式中，f\left(x\right)为连续随机变量\mathit{X}的概率密度函数。

当\alpha =\mathrm{1}时，G\left(x;\alpha \right)退化为原始累积分布函数F\left(x\right)，g\left(x;\alpha \right)退化为原始概率密度函数f\left(x\right)。这一结论说明了NAT变换方法将原始分布拓展得到广义分布。

将q-威布尔的累积分布函数F_q\left(x\right)和概率密度函数f_q\left(x\right)代入NAT变换方法中，得出一个服从NAT变换的四参数扩展分布，称为修正q-威布尔分布（new alpha transformed q-Weibull distribution，NATQ），记为\mathit{X}\sim \mathrm{N}\mathrm{A}\mathrm{T}\mathrm{Q}\left(q,\theta ,\lambda ,\alpha \right)。累积分布函数和概率密度函数分别为

\begin{array}{l}{F}_{\mathrm{N}\mathrm{A}\mathrm{T}\mathrm{Q}}\left(x\right)=\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\times \\ {\left\{{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}+{\left(\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\right)}^{\alpha }\right\}}^{-\mathrm{1}},\end{array}\begin{array}{l}{f}_{\mathrm{N}\mathrm{A}\mathrm{T}\mathrm{Q}}\left(x\right)=(\mathrm{2}-q)\lambda \theta x^{\theta -\mathrm{1}}{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{1}}{\mathrm{1}-q}}\times \\ {\left(\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\right)}^{\alpha -\mathrm{1}}\left(\alpha {\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}+\right.\\ \left.\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\right)\left({\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}+\right.\\ {\left.{\left(\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\right)}^{\alpha }\right)}^{-\mathrm{2}}\mathrm{。}\end{array}

根据可靠度函数以及失效率函数的定义，可以得到NATQ分布的可靠度函数以及失效率函数分别为

\begin{array}{l}{R}_{\mathrm{N}\mathrm{A}\mathrm{T}\mathrm{Q}}\left(x\right)={\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\times \\ {\left\{{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}+{\left(\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\right)}^{\alpha }\right\}}^{-\mathrm{1}},\end{array}\begin{array}{l}{h}_{\mathrm{N}\mathrm{A}\mathrm{T}\mathrm{Q}}\left(x\right)=(\mathrm{2}-q)\lambda \theta x^{\theta -\mathrm{1}}{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{1}}{\mathrm{1}-q}}\times \\ {\left(\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\right)}^{\alpha -\mathrm{1}}{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{q-\mathrm{2}}{\mathrm{1}-q}}\times \\ \left(\alpha {\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}+\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\right)\times \\ {\left({\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}+{\left(\mathrm{1}-{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda x^{\theta }\right]}^{\frac{\mathrm{2}-q}{\mathrm{1}-q}}\right)}^{\alpha }\right)}^{-\mathrm{1}}\mathrm{。}\end{array}

NATQ分布的PDF和CDF关于参数\alpha的变化趋势见图1。假设q=\mathrm{1.5}、\theta =\mathrm{3.5}、\lambda =\mathrm{0.1}，以\alpha =\mathrm{1}为对照，可知参数\alpha对曲线的左侧产生明显的影响。由图1（a）可知，随着参数\alpha减小，PDF的峰值变小且左移，曲线转折点明显左移，分布前端曲线饱满。由图1（b）可知，随着参数\alpha减小，CDF的曲线转折点明显左移。当x\to \mathrm{\infty }时，极限为1，满足累积分布函数要求。

通过对两个人工数据集进行拟合（采用本文第2章所提参数估计方法进行拟合），NATQ分布与q-威布尔分布对比结果见图2。图2（a）和图2（b）分别表示分布前端概率较低与较高的情况，NATQ分布均有较好的拟合效果，且拟合曲线不再呈单一的变化趋势，具有较高的灵活性。

NATQ分布的失效率和可靠度函数见图3。给定不同的参数值，可以得出不同类型的失效率函数，如递增-递减-递增、浴盆型、单峰型、单调递增、单调递减和常数型。

NATQ分布的PDF关于参数q、\theta和\lambda的变化趋势见图4。

q和\theta为形状参数，分别控制分布的尾部和峰度。当q增大时，描述拖尾现象。当\theta增大时，峰值变大。\lambda为尺度参数，控制分布的尺度范围，当\lambda增大时，数据范围缩小。

2 修正q-威布尔分布的参数估计

2.1 DBO-MLE参数估计

设\mathit{X}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n\right)是总体服从NATQ分布且样本量为n的随机变量，参数为q、\theta、\lambda和\alpha。假设\mathit{\epsilon }=\left(q,\theta ,\lambda ,\alpha \right)，则极大似然估计（maximum likelihood estimation，MLE）中似然函数^[17]可以表示为

L_n\left(\mathit{\epsilon }\right)=\stackrel{n}{\prod _{i=\mathrm{1}}}{f}_{\mathrm{N}\mathrm{A}\mathrm{T}\mathrm{Q}}\left(x_i;\mathit{\epsilon }\right)。

假设对数似然函数为{\mathit{𝓁}}_n=\mathrm{l}\mathrm{o}\mathrm{g}\left(L_n\left(\mathit{\epsilon }\right)\right)，则

\begin{array}{l}{\mathit{𝓁}}_n={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{\infty }}}\left\{\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{2}-q\right)+\mathrm{l}\mathrm{o}\mathrm{g}\lambda +\left(\theta -\mathrm{1}\right)\mathrm{l}\mathrm{o}\mathrm{g}\left(x_i\right)+\right.\\ \begin{array}{c}\end{array}\mathrm{l}\mathrm{o}\mathrm{g}\theta +\frac{\mathrm{l}\mathrm{o}\mathrm{g}{T}_i}{\mathrm{2}-q}+\left(\alpha -\mathrm{1}\right)\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{1}-T_i\right)+\end{array}

\begin{array}{c}\end{array}\left.\mathrm{l}\mathrm{o}\mathrm{g}\left[\left(\alpha -\mathrm{1}\right)T_i+\mathrm{1}\right]-\mathrm{2}\mathrm{l}\mathrm{o}\mathrm{g}\left[T_i+{\left(\mathrm{1}-T_i\right)}^{\alpha }\right]\right\} \mathrm{。}（1）

式中，T_i={\left[\mathrm{1}-\left(\mathrm{1}-q\right)\lambda {x_i}^{\theta }\right]}^{\left(\mathrm{2}-q\right)/\left(\mathrm{1}-q\right)}。

如果{\mathit{𝓁}}_n的一阶偏导数存在，则q、\theta、\lambda和\alpha的极大似然估计可以通过解下列方程得到。

\frac{\partial {\mathit{𝓁}}_n}{\partial q}=\mathrm{0},\frac{\partial {\mathit{𝓁}}_n}{\partial \theta }=\mathrm{0},\frac{\partial {\mathit{𝓁}}_n}{\partial \lambda }=\mathrm{0},\frac{\partial {\mathit{𝓁}}_n}{\partial \alpha }=\mathrm{0}。

由于对数似然函数的一阶导数是非线性的，利用经典的数值优化方法来求解其参数复杂且低效，启发式优化方法可以用来解决这个问题。故本文采用蜣螂优化算法来求解似然方程的估计值。

蜣螂优化算法（dung beetle optimisation，DBO）是一种基于昆虫蜣螂行为的优化算法，与传统的优化算法不同，DBO能够快速收敛，具有高精度、高稳定性的特点。该算法主要通过蜣螂的位置变化，利用滚球、跳舞、繁殖、觅食和偷窃5种行为获得全局最优权值和阈值^[18]。本文将蜣螂优化算法与极大似然估计相结合，构建分布的参数估计方法，记作DBO-MLE参数估计模型。

首先，假设蜣螂种群数量为N，并将种群划分为4个部分，即N_{\mathrm{1}}、N_{\mathrm{2}}、N_{\mathrm{3}}和N_{\mathrm{4}}。其次，当蜣螂个体i\in N_{\mathrm{1}}时，给定随机数\delta \in \left(\mathrm{0,1}\right)，当时，蜣螂进行直线滚球行为；当\delta \ge \mathrm{0.9}时，通过跳舞重新定位。获得局部最优解{\mathit{x}}_{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t。当蜣螂个体i\in N_{\mathrm{2}}时，根据{\mathit{x}}_{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t形成繁殖区域并进行繁殖。当蜣螂个体i\in N_{\mathrm{3}}时，根据全局最优解{\mathit{x}}_{\mathrm{l}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t形成觅食区域并进行觅食。当蜣螂个体i\in N_{\mathrm{4}}时，根据{\mathit{x}}_{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t和{\mathit{x}}_{\mathrm{l}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t形成竞争食物的最优位置。最后，判断是否为可行解，计算并比较所有蜣螂的适应度值，更新全局最优解。

（1） 当蜣螂滚球（无障碍向前）时，位置更新公式为

{\mathit{x}}_i^{t+\mathrm{1}}={\mathit{x}}_i^t+ak{\mathit{x}}_i^{t-\mathrm{1}}+b|{\mathit{x}}_i^t-{\mathit{x}}_{\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{t}}^t|，（2）

式中：{\mathit{x}}_i^t=\left(q_i^t,{\theta }_i^t,{\lambda }_i^t,{\alpha }_i^t\right)为在NATQ分布参数形成的四维空间内第i个蜣螂在第t次迭代时的位置信息；a为自然因素，可以使蜣螂偏离原来的方向，其中，a=\mathrm{1}表示无偏差，a=-\mathrm{1}表示偏离原方向；k为滚球蜣螂的偏转系数，k\in (\mathrm{0,0.2}]；b为随机数，b\in (\mathrm{0,1})；{\mathit{x}}_{\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{t}}^t为全局最差位置；|{\mathit{x}}_i^t-{\mathit{x}}_{\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{t}}^t|为一个向量，表示向量相减后取绝对值，模拟光强变化。

（2） 当蜣螂遇到障碍无法前进时，需要通过跳舞来重新定位，位置更新公式为

{\mathit{x}}_i^{t+\mathrm{1}}={\mathit{x}}_i^t+\mathrm{t}\mathrm{a}\mathrm{n}\sigma \left|{\mathit{x}}_i^t-{\mathit{x}}_i^{t-\mathrm{1}}\right|，（3）

式中：\sigma为偏转角，\sigma \in \left[\mathrm{0},\mathrm{\pi }\right]，当\sigma为0或\mathrm{\pi }/\mathrm{2}时，正切函数无意义，表示蜣螂的位置不发生改变；|{\mathit{x}}_i^t-{\mathit{x}}_i^{t-\mathrm{1}}|为一个向量，表示第i只蜣螂在第t次迭代时与第t-\mathrm{1}次迭代时位置之间的偏移量。

（3） 采用边界选择策略来模拟蜣螂繁殖区域，位置更新公式为

{\mathit{x}}_i^{t+\mathrm{1}}={\mathit{x}}_{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t+{\mathit{b}}_{\mathrm{1}}\left({\mathit{x}}_i^t-{\mathit{L}}_{{\mathrm{b}}^{\mathrm{*}}}\right)+{\mathit{b}}_{\mathrm{2}}\left({\mathit{x}}_i^t-{\mathit{U}}_{{\mathrm{b}}^{\mathrm{*}}}\right)，（4）

式中：{\mathit{b}}_{\mathrm{1}}和{\mathit{b}}_{\mathrm{2}}为两个\mathrm{1}\times d维的相互独立的随机向量；d为优化问题的维数；{\mathit{L}}_{{\mathrm{b}}^{\mathrm{*}}}和{\mathit{U}}_{{\mathrm{b}}^{\mathrm{*}}}为繁殖区域的下界和上界。

（4） 当蜣螂觅食时，位置更新公式为

{\mathit{x}}_i^{t+\mathrm{1}}={\mathit{x}}_i^t+C_{\mathrm{1}}\left({\mathit{x}}_i^t-{\mathit{L}}_{{\mathrm{b}}^{\mathrm{l}}}\right)+{\mathit{C}}_{\mathrm{2}}\left({\mathit{x}}_i^t-{\mathit{U}}_{{\mathrm{b}}^{\mathrm{l}}}\right)，（5）

式中：{\mathit{L}}_{{\mathrm{b}}^{\mathrm{l}}}和{\mathit{U}}_{{\mathrm{b}}^{\mathrm{l}}}为觅食区域的下界和上界；C_{\mathrm{1}}为服从正态分布的随机数；{\mathit{C}}_{\mathrm{2}}为\mathrm{1}\times d维的\left(\mathrm{0,1}\right)上的随机向量。

（5） 蜣螂在偷窃时，位置更新公式为

{\mathit{x}}_i^{t+\mathrm{1}}={\mathit{x}}_{\mathrm{l}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t+Q\mathit{g}\left(\left|{\mathit{x}}_i^t-{\mathit{x}}_{\mathrm{g}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t\right|+\left|{\mathit{x}}_i^t-{\mathit{x}}_{\mathrm{l}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t\right|\right)，（6）

式中：{\mathit{x}}_{\mathrm{l}\mathrm{b}\mathrm{e}\mathrm{s}\mathrm{t}}^t为最佳食物来源，即全局最优解；Q为常数，其取值范围为\{\mathrm{0.1,0.5,1},\mathrm{1.5,2}\}；\mathit{g}为\mathrm{1}\times d维的符合正态分布的随机向量。

为解决对数似然函数相关的非线性微分方程组求解复杂以及解析解难以得到的问题，采用蜣螂优化算法来求解MLE问题，完成分布参数估计。

对于NATQ分布，DBO算法中使用的适应度函数为NATQ分布的对数似然函数，即式（1）。通过选择最优的参数估计值\widehat{q}、\widehat{\theta }、\widehat{\lambda }和\widehat{\alpha }，使适应度函数值最大化。适应度函数需要满足当时，，以及约束条件\lambda >\mathrm{0}、\theta >\mathrm{0}、x_i>\mathrm{0}、、\left(\alpha -\mathrm{1}\right)T_i+\mathrm{1}>\mathrm{0}、T_i+{\left(\mathrm{1}-T_i\right)}^{\alpha }>\mathrm{0}、\left(\mathrm{2}-q\right)>\mathrm{0}。

在NATQ分布MLE问题的求解过程中，不允许种群中存在不可行解，如果生成的解为不可行解，则认为它不如当前的解，并将它丢弃，而保留当前解。算法流程见图5。

2.2 蒙特卡罗模拟

对所提出的DBO-MLE参数估计进行编码，并通过仿真实验来评价其性能。针对可靠性工程中的6类典型应用场景，设计不同的失效率类型。q、\theta、\lambda和\alpha的值以及相应的失效率函数类型见表1。

假设样本大小n分别为100、300和500。由逆变换法生成服从NATQ分布的随机数。蒙特卡罗模拟500次，模拟步骤如下。

步骤1 选择数值N，即蒙特卡罗模拟次数。

步骤2 选择参数\mathit{\epsilon }=\left(q,\theta ,\lambda ,\alpha \right)的值。

步骤3 生成服从NATQ分布且样本量为n的随机数。

步骤4 用DBO-MLE方法进行参数估计得到估计值\widehat{\mathit{\epsilon }}。

步骤5 重复N次步骤3和步骤4。

步骤6 计算均值（MEAN）、偏差（Bias）和均方误差（MSE），表达式分别为

V_{\mathrm{M}\mathrm{E}\mathrm{A}\mathrm{N}}\left(\widehat{\mathit{\epsilon }}\right)=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}{\widehat{\mathit{\epsilon }}}_n^i，V_{\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}}\left(\widehat{\mathit{\epsilon }}\right)=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}\left({\widehat{\mathit{\epsilon }}}_n^i-\mathit{\epsilon }\right)，

E_{\mathrm{M}\mathrm{S}\mathrm{E}}\left(\widehat{\mathit{\epsilon }}\right)=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}{\left({\widehat{\mathit{\epsilon }}}_n^i-\mathit{\epsilon }\right)}^{\mathrm{2}}，

式中：\mathit{\epsilon }表示待估计参数值；{\widehat{\mathit{\epsilon }}}_n^i表示样本量为n，且第i次蒙特卡罗模拟的估计结果。

蒙特卡罗模拟结果见表2。

由表2可以看出，采用蜣螂优化算法与极大似然估计相结合的方法进行参数估计，所得到的分布参数估计值与真实值接近，随着样本数量n的增加，参数估计的Bias及MSE均逐渐减小且接近于0。说明DBO-MLE参数估计方法能够有效地捕捉数据的特征，并准确估计出参数值。

3 实例分析

通过3个真实的露天矿矿用卡车故障数据集，将NATQ分布与q-威布尔分布、q-指数威布尔分布（q-exponential-Weibull distribution，qE-W）、托普-莱昂威布尔分布（Topp Leone Weibull distribution，TLWD）和三参数扩展奇威布尔指数分布（extended odd Weibull exponential distribution，EOWEx）对比。

q-指数威布尔分布^[19]的概率密度函数为

f(x)=\left(\mathrm{2}-q\right)\alpha \lambda \gamma {\mathrm{e}}^{\lambda x^{\gamma }}{x}^{\gamma -\mathrm{1}}{\left[\mathrm{1}-\left(\mathrm{1}-q\right)\alpha \left({\mathrm{e}}^{\lambda x^{\gamma }}-\mathrm{1}\right)\right]}^{\frac{\mathrm{1}}{\mathrm{1}-q}}。

托普-莱昂威布尔分布^[20]的概率密度函数为

f(x)=\mathrm{2}\lambda \alpha \theta x^{\alpha -\mathrm{1}}{\mathrm{e}}^{-\mathrm{2}\theta x^{\alpha }}{\left(\mathrm{1}-{\mathrm{e}}^{-\mathrm{2}\theta x^{\alpha }}\right)}^{\lambda -\mathrm{1}}。

三参数扩展奇威布尔指数分布^[21]的概率密度函数为

f(x)=\alpha \lambda {\mathrm{e}}^{\alpha \lambda x}{\left(\mathrm{1}-{\mathrm{e}}^{-\lambda x}\right)}^{\alpha -\mathrm{1}}{\left[\mathrm{1}+\beta {\left({\mathrm{e}}^{\lambda x}-\mathrm{1}\right)}^{\alpha }\right]}^{\frac{-\mathrm{1}}{\beta }-\mathrm{1}}。

为了衡量分布拟合优良性，结合分布拟合度和参数数量，采用负对数似然（negative log-likelihod，NLL）、改正的赤池信息准则（consistent Akaike information criterion，CAIC）和Hannan-Quin信息准则（Hannan-Quin information criterion，HQIC）评估模型拟合效果，表达式分别为

V_{\mathrm{N}\mathrm{L}\mathrm{L}}=-L，

V_{\mathrm{C}\mathrm{A}\mathrm{I}\mathrm{C}}=\mathrm{2}{V}_{\mathrm{N}\mathrm{L}\mathrm{L}}+\frac{\mathrm{2}kn}{n-k-\mathrm{1}}，

V_{\mathrm{H}\mathrm{Q}\mathrm{I}\mathrm{C}}=\mathrm{2}{V}_{\mathrm{N}\mathrm{L}\mathrm{L}}+\mathrm{2}q\mathrm{l}\mathrm{o}\mathrm{g}\left(\mathrm{l}\mathrm{o}\mathrm{g}n\right)，

式中：L为最大化的对数似然函数值；n为样本量，k为参数的个数。

3个信息准则的值越小，模型的拟合效果越好。

数据集1为设备编号030103^#的630E型矿用卡车2018年1月起至2020年10月的机械故障间隔时间数据。数字特征见表3。

数据集2为设备编号030415^#的SF33900型矿用卡车2019年3月起至2020年11月的液压故障间隔时间数据。数字特征见表4。

数据集3为设备编号030418^#的SF33900型矿用卡车2018年3月起至2020年11月的焊修故障间隔时间数据。数字特征见表5。

3个数据集中每个模型的参数估计值、NLL、CAIC和HQIC结果见表6～表8，个数据集中每个模型拟合的PDF和CDF见图6。由上述结果可知，NATQ分布的3种信息准则数值均最小，NATQ分布拟合效果更好。分别将3个数据集中NATQ分布参数估计值代入可靠度函数，得到

R_{\mathrm{1}}\left(x\right)=\left\{{\left[\mathrm{1}-{\left(\mathrm{1}-\mathrm{0.655}\mathrm{0}{x}^{\mathrm{0.103}\mathrm{4}}\right)}^{\mathrm{1.154}\mathrm{7}}\right]}^{\mathrm{1.395}\mathrm{9}}+\right.

{\left.{\left(\mathrm{1}-\mathrm{0.655}\mathrm{0}{x}^{\mathrm{0.103}\mathrm{4}}\right)}^{\mathrm{1.154}\mathrm{7}}\right\}}^{-\mathrm{1}}{\left(\mathrm{1}-\mathrm{0.655}\mathrm{0}{x}^{\mathrm{0.103}\mathrm{4}}\right)}^{\mathrm{1.154}\mathrm{7}}，（7）

\begin{array}{l}{R}_{\mathrm{2}}\left(x\right)={\left(\mathrm{1}+\mathrm{1.751}\mathrm{0}\times {\mathrm{10}}^{-\mathrm{8}}{x}^{\mathrm{5.368}\mathrm{9}}\right)}^{-\mathrm{5.035}\mathrm{0}}\times \\ \left\{{\left[\mathrm{1}-{\left(\mathrm{1}+\mathrm{1.751}\mathrm{0}\times {\mathrm{10}}^{-\mathrm{8}}{x}^{\mathrm{5.368}\mathrm{9}}\right)}^{-\mathrm{5.035}\mathrm{0}}\right]}^{\mathrm{0.361}\mathrm{9}}+\right.\end{array}

{\left.{\left(\mathrm{1}+\mathrm{1.751}\mathrm{0}\times {\mathrm{10}}^{-\mathrm{8}}{x}^{\mathrm{5.368}\mathrm{9}}\right)}^{-\mathrm{5.035}\mathrm{0}}\right\}}^{-\mathrm{1}}，（8）

R_{\mathrm{3}}\left(x\right)=\left\{{\left[\mathrm{1}-{\left(\mathrm{1}-\mathrm{0.444}\mathrm{6}{x}^{\mathrm{0.186}\mathrm{8}}\right)}^{\mathrm{1.036}\mathrm{2}}\right]}^{\mathrm{17.580}\mathrm{3}}+\right.

{\left.{\left(\mathrm{1}-\mathrm{0.444}\mathrm{6}{x}^{\mathrm{0.186}\mathrm{8}}\right)}^{\mathrm{1.036}\mathrm{2}}\right\}}^{-\mathrm{1}}{\left(\mathrm{1}-\mathrm{0.444}\mathrm{6}{x}^{\mathrm{0.1868}}\right)}^{\mathrm{1.036}\mathrm{2}}。（9）

可靠性评价指标取故障间隔时间（mean time between failures，MTBF），其计算式为

V_{\mathrm{M}\mathrm{T}\mathrm{B}\mathrm{F}}={\int }_{\mathrm{0}}^{\mathrm{\infty }}x{f}_{\mathrm{N}\mathrm{A}\mathrm{T}\mathrm{Q}}\left(x\right)\mathrm{d}x。

计算得到3个数据集的MTBF分别为10.858 6 d、17.034 2 d、39.275 1 d。取不同的可靠度R\left(x\right)代入式（7）、式（8）和式（9），得到预防性维修周期，见表9。

由表9可知，矿用卡车的预防性维修周期随着R\left(x\right)的减小而增大。当数据集1、2、3在可靠度分别为0.4、0.6、0.5时，其预防性维修周期均小于对应的故障间隔时间。因此为保证矿用卡车的可靠性，分别取可靠度R_{\mathrm{1}}\left(x\right)=\mathrm{0.4}、R_{\mathrm{2}}\left(x\right)=\mathrm{0.6}和R_{\mathrm{3}}\left(x\right)=\mathrm{0.5}为预防性维修周期。

4 结论

为了更加准确地研究露天矿矿用卡车故障数据的失效率规律，提高其可靠性分析的准确性，提出了一种基于修正q-威布尔分布的可靠性分析模型。通过推理与验证，得出以下结论。

（1）修正q-威布尔分布能有效控制分布的前端变化，允许分布的前端和尾端同时具有更大的灵活性。

（2）通过蒙特卡罗模拟验证了蜣螂优化算法与极大似然估计相结合的参数估计方法的有效性。

（3）实例对比分析表明，修正q-威布尔分布在露天矿卡车数据集上的拟合效果最优。

（4）基于修正q-威布尔分布模型，计算得出矿用卡车3个数据集的故障间隔时间以及预防性维修周期。

在后续研究过程中，可以考虑将矿用卡车分为多个子系统，将故障数据类型更加细分，针对故障多发的子系统分别进行可靠性分析，为设备稳定高效运行提供科学有效的基础决策数据。

参考文献

[1] 付恩三,刘光伟,赵浩,智慧露天矿山总体框架及关键技术研究[J].工矿自动化,2021,47(8):27-32.

[2] FU Ensan, LIU Guangwei, ZHAO Hao,et al.Research on the overall framework and key technologies of intelligent open-pit mines[J].Industry and Mine Automation,2021,47(8):27-32.

[3] 刘威,郭直清,黄敏,基于MCMC算法的露天矿典型设备故障预测方法研究[J].煤炭科学技术,2019,47(10):51-57.

[4] LIU Wei, GUO Zhiqing, HUANG Min,et al.Research on fault prediction method of typical equipment in open-pit mine based on MCMC algorithm[J].Coal Science and Technology,2019,47(10):51-57.

[5] RAHIMDEL M J.Bayesian network approach for reliability analysis of mining trucks[J].Scientific Reports,2024,14(1):3415.

[6] HE Y S, KUSIAK A, OUYANG T H,et al.Data-driven modeling of truck engine exhaust valve failures: a case study[J].Journal of Mechanical Science and Technology,2017,31(6):2747-2757.

[7] QARAHASANLOU A N, ATAEI M, KHAOLUKAKAIE R,et al.Maintainability measure based on operating environment,a case study: Sungun copper mine[J].Journal of Mining and Environment,2018,8(3):511-521.

[8] TORAMAN S.System reliability analysis of large capacity electric mining trucks used in coal mining[J].Journal of Reliability and Statistical Studies,2023:81-98.

[9] YUAN F Q, BARABADI A, LU J M.Reliability modelling on two-dimensional life data using bivariate Weibull distribution:with case study of truck in mines[J].Eksploatacja i Niezawodnosc-Maintenance and Reliability,2017,19(4):650-659.

[10] KUMAR D, GUPTA S, YADAV P K.Reliability,availability and maintainability (RAM) analysis of a dragline[J].Journal of Mines,Metals and Fuels,2020,68(2):68-77.

[11] MONIRI-MORAD A, SATTARVAND J.A comparative study between the system reliability evaluation methods: case study of mining dump trucks[J]. Journal of Engineering and Applied Science,2023,70(1):103.

[12] ASSIS E M, BORGES E P, VIEIRA DE MELO S A B.Generalized q-Weibull model and the bathtub curve[J].International Journal of Quality & Reliability Management,2013,30(7):720-736.

[13] AHMAD Z, ALMASPOOR Z, KHAN F,et al.On predictive modeling using a new flexible Weibull distribution and machine learning approach:analyzing the COVID-19 data[J].Mathematics,2022,10(11):1792.

[14] PICOLI S, MENDES R S, MALACARNE L C. q-exponential,Weibull,and q-Weibull distributions:an empirical analysis[J].Physica A:Statistical Mechanics and Its Applications,2003,324(3/4):678-688.

[15] 樊立天,王浩文,凌庆庆,基于q-Weibull分布的医疗设备可靠性分析[J].机械强度,2023,45(2):392-398.

[16] FAN Litian, WANG Haowen, LING Qingqing,et al.Reliability analysis of medical equipment based on q-Weibull distribution[J].Journal of Mechanical Strength,2023,45(2):392-398.

[17] MACHADO DE ASSIS E, LIMA G A C, PRESTES A,et al. q-Weibull applied to Brazilian hydropower equipment[J].IEEE Transactions on Reliability,2019,68(1):122-132.

[18] MARCHETTI S, DI MAIO F,ZIO E,et al.Key performance indicators of aging safety barriers in oil and gas facilities[C]//Proceeding of the 33rd European Safety and Reliability Conference.September 3-7,2023.Research Publishing Services,2023:2230-2235.

[19] TSALLIS C.Possible generalization of Boltzmann-Gibbs statistics[J].Journal of Statistical Physics,1988,52(1):479-487.

[20] ALMONGY H M, ALMETWALLY E M, ALJOHANI H M,et al.A new extended Rayleigh distribution with applications of COVID-19 data[J].Results in Physics,2021,23:104012.

[21] XUE J K, SHEN B.Dung beetle optimizer:a new meta-heuristic algorithm for global optimization[J].The Journal of Supercomputing, 2023,79(7):7305-7336.

[22] SUNDARAM N, JAYAKODI G.A study on statistical properties of a new class of q-exponential-Weibull distribution with application to real-life failure time data[J].Reliability: Theory and Applications,2023,18(3):582-595.

[23] TUOYO D O, OPONE F C, EKHOSUEHI N.The Topp-Leone Weibull distribution:its properties and application[J].Earthline Journal of Mathematical Sciences,2021,7(2):381-401.

[24] AFIFY A Z, MOHAMED O A.A new three-parameter exponential distribution with variable shapes for the hazard rate:estimation and applications[J].Mathematics,2020,8(1):135.

基金资助

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

辽宁省教育厅基本科研项目(LJ212410147013)

辽宁省“兴辽英才计划”项目(XLYC2211085)