CN104239694A - 一种城轨列车转向架的故障预测与视情维修方法 - Google Patents

Abstract

本发明提供一种城轨列车转向架的故障预测与视情维修方法，采用基于生存分析的方法确定转向架的构架、弹簧装置、连接装置、轮对和轴箱、驱动机构、基础制动装置的最佳寿命分布模型，并得出各子系统的可靠性特征函数，然后采用进化算法优化的神经网络模型对转向架各子系统的故障率进行计算。最后，将转向架各子系统的故障率和安全运营天数作为协变量进行比例风险模型建模，在成本优化的基础上，得出转向架系统视情维修的阈值和控制限。

Description

一种城轨列车转向架的故障预测与视情维修方法

技术领域

本发明涉及一种转向架的故障预测与视情维修方法，尤其涉及一种城轨列车转向架的故障预测与视情维修方法。 

背景技术

城市轨道列车转向架系统是城轨列车的重要部件，用来传递各种载荷，并利用轮轨间的黏着来保证牵引力的产生，其性能很大程度上决定了车辆运行的稳定性和安全性。因此，如何准确地评定城轨列车转向架系统的可靠性，预测转向架系统的故障率，对于城轨列车运行的稳定性和安全性具有重要意义。同时，转向架系统在使用过程中，需要及时进行检修、维护、保养，其产生的费用是城轨列车维护养护费用的主要组成部分，因此准确地评定城轨列车转向架系统的运行状态，制定合理的转向架系统维修计划，对于保证城轨列车安全经济地运行具有重要意义。 

研究系统可靠性的方法主要有统计模型法和机理模型法。统计模型法主要是通过分析系统失效时间的概率分布函数，并得出系统相关的可靠性指标，如平均故障间隔时间、失效率等。机理模型法主要是通过对系统的结构原理进行可靠性分析。在维修决策建模的研究，现有技术中提出综合考虑了运行条件、负载大小等因素对系统寿命带来的影响，建立旋转机械的维修模型并取得了较好的效果，表明了协变量在对系统寿命分布规律研究中的重要作用。 

以上研究在转向架系统可靠性分析和视情维修上都取得了一定的成果，但是，仍然一定的问题：(1)目前上述研究大多针对铁路客车、货车，而对城轨列车转向架可靠性的研究较少；(2)上述方法基于监控或检修数据，有时不能准确地掌握转向架的寿命时间，在时间的不确定问题上缺乏考虑；(3)对于传统的系统寿命或系统性能衰退数学描述通常是基于一维时间变量来衡量的，但是对于转向架等复杂系统而言，维修阈值的模糊性以及失效模式的随机性等因素使得仅依据一维时间变量无法客观地予以评估，必须同时考虑性能参数、失效类型等相关协变量对系统运行的影响，为制定维修优化决策提供更有力的依据。 

发明内容

本发明对转向架系统的各关键部件进行可靠性分析，通过对转向架系统关键部件历史检修数据的整理和分析，建立故障率预测模型，并针对整个转向架系统建立比例风险模型。 

具体采取如下技术方案： 

一种城轨列车转向架视情维修方法，转向架包括构架、弹簧装置、连接装置、轮对和轴箱、驱动机构、基础制动装置六个子系统，该方法依次包括如下步骤：(1)根据采集的历史故障数据进行删失处理，基于生存分析的方法确定转向架 各子系统的最佳寿命分布模型，并得出各子系统的可靠性特征函数，计算出各子系统的可靠度，将各子系统中可靠度最低的子系统确定为转向架最薄弱环节，进行对该最薄弱环节重点关注；(2)采用进化算法优化的神经网络模型计算转向架各子系统的故障率；(3)将转向架各子系统的安全运营天数和计算得到的故障率作为协变量进行比例风险模型建模，得出转向架系统视情维修的阈值和控制限，控制限上限是失效阈值,在运行过程中，一旦发现系统状态值超过上限，则系统此时为不可用状态，按规定必须进行修复性修理或更换后才重新投入使用；控制限下限是预防性维修或更换阈值，代表系统的潜在故障开始出现，如果系统的状态值超过了下限，则应考虑对部件作相应的排故或预防性维修工作，如果系统的状态值低于下限，则不用考虑维修。 

优选地，步骤(1)包括如下步骤： 

1)建立轮对和轴箱、弹簧装置、连接装置两参数威布尔分布模型， 

故障分布函数为： 

$F\left(t\right)=1-\exp \left\{{-\left[\frac{t}{\mathrm{\η}}\right]}^{\mathrm{\β}}\right\}$

可靠度函数为： 

$R\left(t\right)=\exp \left\{{-\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}}\right\}$

概率密度函数为： 

$f\left(t\right)=\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \left[{-\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}}\right]$

失效率函数为： 

$\mathrm{\λ}\left(t\right)=\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}$

其中，t≥0，β＞0且η＞0，t为故障间隔时间，β和η分别是分布的形状参数和尺度参数； 

2)建立构架、驱动机构、基础制动装置三参数威布尔分布模型， 

故障分布函数为： 

$F\left(t\right)=1-\exp \{-{\left[\frac{t-\mathrm{\γ}}{\mathrm{\η}}\right]}^{\mathrm{\β}}\}$

可靠度函数为： 

$R\left(t\right)=\exp \{-{\left(\frac{t-\mathrm{\γ}}{\mathrm{\η}}\right)}^{\mathrm{\β}}\}$

其概率密度函数为： 

$f\left(t\right)=\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t-\mathrm{\γ}}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \left[{-\left(\frac{t-\mathrm{\γ}}{\mathrm{\η}}\right)}^{\mathrm{\β}}\right]$

失效率函数为： 

$\mathrm{\λ}\left(t\right)=\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t-\mathrm{\γ}}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}$

其中，t为故障间隔时间，β＞0且η＞0，β和η分别是分布的形状参数和尺度参数,γ是分布的位置参数。 

优选地，步骤(2)包括如下步骤：建立的基于PSO的BP神经网络模型，即采用PSO优化BP神经网络连接权值，并对BP神经网络进行训练，在建模过程中，需确定PSO的一些参数，设置如下： 

1)粒子数目 

粒子数目一般取值20～40，对于大多数问题来说，30个粒子就可以取得很好的结果，本节选取粒子数目为30； 

2)粒子最大速率 

粒子最大速率决定粒子在一次飞行中可以移动的最大距离，必须限制粒子最大速率，否则，粒子就可能跑出搜索空间，粒子的最大速率通常设定为粒子范围的宽度，本节选取为0.5； 

3)加速常数 

学习因子c₁和c₂表示粒子受社会知识和个体认知的影响程度，通常设为相同值以给两者同样的权重c₁＝c₂＝2，在找到这两个最优位置时，粒子根据如下的公式来更新自己的速度和位置， 

$v_i^{k+1}=w*v_i^k+c_1{r}_1(p_i^k-x_i^k)+c_2{r}_2(p_g^k-x_i^k)$

$x_i^{k+1}=x_i^k+v_i^{k+1}$

式中，表示第i个粒子在第k+1代时的飞行速度；表示第i个粒子在第k+1代时的位置；表示第i个粒子到第k代的最优位置；表示种群到第k代的最优位置；为个体认知；为种群认知；w为惯性权重；是粒子的速度；c₁和c₂为学习因子；r₁和r₂为均匀分布在[0,1]内的随机数；i＝1,2…N； 

4)适应度函数 

适应度函数用于区分种群中个体好坏的标准，也是进行选择的唯一依据，其值越大越好，根据本文所确定的网络评价指标，此处采用网络输出与期望输出间的均方根误差(RMSE)的倒数为评价标准，其适应度函数为 

f＝1/RMSE； 

5)初始化，将输入输出数据归一化到[-1,1]区间，初始化产生的粒子种群代表不同权重与阈值组合的神经网络，产生初始BP网络结构； 

6)评价，依据4)中公式计算其适应度； 

7)更新位置和速度，通过比较适应度不断更新粒子的位置，适应度最优的个体极值即为全局极值，对应的权值与阈值为粒子种群的当前最优解，更新速度； 

8)算法停止的判断，达到最大迭代次数或适应度达到期望误差，停止迭代； 

9)最优解的生存，迭代停止时，全局权值对应的权值与阈值为训练样本的最优解； 

10)将最优解带入BP神经网络模型中学习，用来对故障率进行预测。 

优选地，步骤(3)包括如下步骤： 

1)建立转向架的比例风险模型，其风险函数为： 

$h\left(t\right|X\left(t\right))=\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp (\mathrm{\γ}\·X\left(t\right))$

求得故障率阈值h^*： 

$h^{*}\≥\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp (\mathrm{\γ}\·X\left(t\right))$

其中，t为故障间隔时间，β＞0且η＞0，β和η分别是分布的形状参数和尺度参数,γ是分布的位置参数，X(t)为协变量； 

2)以工作时间t为横坐标、以的值为纵坐标作图，图形中上线为故障率阈值对应的控制限的上限，下线为故障率阈值对应的控制限的下限。 

本发明具有如下有益效果： 

(1)准确及时地掌握转向架系统关键部件的故障情况、可靠性现状以及可靠性的变化趋势，准确地掌握转向架系统的薄弱部件； 

(2)实现对转向架关键部件故障率进行科学、准确地预测； 

(3)通过转向架系统及其各部件的运行状态，合理地做出维修决策，并将维修成本因素考虑在内，对维修决策进行优化，避免维修成本的浪费，并结合预测的故障率对转向架系统进行预防维修。 

附图说明

图1是本发明主要流程图。 

图2是PSO-BP预测模型建模流程图。 

图3是可修系统失效过程的趋势变化图。 

图4是累计风险函数曲线图。 

图5是视情维修决策控制限示意图。 

图6是PHM模型输出维修决策图。 

具体实施方式

一、转向架各子系统可靠度计算 

如图1所示，首先，准确地掌握转向架系统各关键部件的可靠性趋势，根据采集的历史故障数据进行删失处理，基于生存分析理论对转向架系统进行可靠性分析，根据各关键部件故障数据得出最佳的分布模型，及时掌握各关键部件的可靠性趋势，并通过对比分析，得出转向架系统的薄弱环节；然后，对转向架系统各关键部件故障率进行预测模型研究，收集各关键部件的历史检修数据，进行整理分析，以可靠性较低部件的故障率采用粒子群算法优化的BP神经网络(PSO-BP)建立故障率预测模型；最后，根据检修数据对转向架系统视情维修决策。将转向架各部件的状态参数作为协变量，筛选出显著影响转向架运行安全的协变量，将成本因素考虑在内，对转向架视情维修的阈值和控制限进行优化，并结合预测的故障率对模型的有效性进行验证。 

删失数据统计：在转向架系统的故障数据统计中，主要包括两类删失数据：一类是区间删失数据，假定检修工作可靠，则故障发生在本次检修和上次检修之间，故障具体时间未知，例如序号2的时间为(35,36]，是指在第35天检修时未发现故障而第36天检修时发生故障，则表明故障发生在第35天和第36天之间，故障时间为一个区间值、不确定值，具体故障时间未知；一类是右删失数据，统计期的开始和结尾会有截尾数据，故障时间大于跟踪到的某个值，具体未知，例如序号59的时间为(43,]，表明在统计期结束时，轮对轴箱系统还未发生故障，故障时间大于43天而具体故障时间未知，轮对轴箱系统的故障时间删失数据如表3-2所示。 

表3-2故障时间的删失数据 

各子系统分布模型的确定：首先将轮对轴箱的故障数据导入Minitab的工作表，选用可靠性/生存统计工具，采用极大似然估计法分别对指数分布、对数正态分布、两参数Weibull分布和三参数Weibull分布进行了参数估计及拟合度检验，计算得到的参数值如表3-3所示。 

表3-3轮对轴箱系统寿命分布模型 

部分故障数据点落在两参数Weibull分布、对数正态分布、指数分布的置信区间外，拟合程度最差的为指数分布，三参数Weibull分布拟合效果较好，结合表3-3中的A-D统计量可知，三参数Weibull的A-D统计量最小，为0.464，拟合效果最好，所以轮对轴箱系统最佳分布为三参数Weibull分布。 

同理，可以确定构架、弹簧装置、连接装置、驱动机构和基础制动装置的故障分布模型。 

可靠性特征值是用来评价设备可靠性的指标，主要包括累计分布函数、可靠度函数、概率密度函数、故障率函数、平均故障间隔时间等。 

表3-5中给出了转向架各关键子系统的平均故障间隔时间MTBF。结果显示基础制动装置的平均故障间隔时间最高，约为71天，而弹簧装置和轮对轴箱较低，分别为36.38天、37.02天。 

表3-5各关键子系统MTBF 

(一)轮对轴箱系统可靠性特征值 

轮对轴箱服从两参数Weibull分布模型，其形状参数β为1.30274，尺度参数η为6.65490，阈值γ为30.8751，故可得故障时间概率密度函数f(t)、累计分布函数F(t)、可靠度函数R(t)和失效率函数λ(t)如下所示： 

故障时间概率密度函数为： 

$f\left(t\right)=\frac{1.30274}{{6.65490}^{1.30274}}{(t-30.8751)}^{0.30274}{e}^{-{\left(\frac{t-30.8751}{6.65490}\right)}^{1.30274}}$

累计分布函数为： 

$F\left(t\right)=1-e^{-{\left(\frac{t-30.8751}{6.65490}\right)}^{1.30274}}$

可靠度函数为： 

$R\left(t\right)=e^{-{\left(\frac{t-30.8751}{6.65490}\right)}^{1.30274}}$

失效率函数为： 

$\mathrm{\λ}\left(t\right)=\frac{1.30274}{6.65490}{\left(\frac{t-30.8751}{6.65490}\right)}^{0.30274}$

式中，t表示当前无故障工作的时间。 

(二)构架系统可靠性特征值 

构架系统服从三参数Weibull分布模型，其形状参数β为1.11033，尺度参数η为17.3204，阈值γ为30.5765，故可得故障时间概率密度函数f(t)、累计分布函数F(t)、可靠度函数R(t)和失效率函数λ(t)如下所示。 

故障时间概率密度函数为： 

$f\left(t\right)=\frac{1.11033}{{17.3204}^{1.11033}}{(t-30.5765)}^{0.11033}{e}^{-{\left(\frac{t-0.5765}{17.3204}\right)}^{1.11033}}$

累计分布函数为： 

$F\left(t\right)=1-e^{-{\left(\frac{t-30.5765}{17.3204}\right)}^{1.11033}}$

可靠度函数为： 

$R\left(t\right)=e^{-{\left(\frac{t-30.5765}{17.3204}\right)}^{1.11033}}$

失效率函数为： 

$\mathrm{\λ}\left(t\right)=\frac{1.11033}{17.3204}{\left(\frac{t-30.5765}{17.3204}\right)}^{0.11033}$

(三)弹簧装置系统可靠性特征值 

弹簧装置系统服从两参数Weibull分布模型，其形状参数β为7.83533，尺度参数η为38.6704，故可得故障时间概率密度函数f(t)、累计分布函数F(t)、可靠度函数R(t)和失效率函数λ(t)如下所示。 

故障时间概率密度函数为： 

$f\left(t\right)=\frac{7.83533}{{38.6704}^{7.83533}}{t}^{6.83533}{e}^{-{\left(\frac{t}{38.6704}\right)}^{7.83533}}$

累计分布函数为： 

$F\left(t\right)=1-e^{-{\left(\frac{t}{38.6704}\right)}^{7.83533}}$

可靠度函数为： 

$R\left(t\right)=e^{-{\left(\frac{t}{38.6704}\right)}^{7.83533}}$

失效率函数为： 

$\mathrm{\λ}\left(t\right)=\frac{7.83533}{38.6704}{\left(\frac{t}{38.6704}\right)}^{6.83533}$

(四)连接装置系统可靠性特征值 

连接装置系统服从两参数Weibull分布模型，其形状参数β为3.20939，尺度参数η为61.6529，故可得故障时间概率密度函数f(t)、累计分布函数F(t)、可靠度函数R(t)和失效率函数λ(t)如下所示。 

故障时间概率密度函数为： 

$f\left(t\right)=\frac{3.20939}{{61.6529}^{3.20939}}{t}^{2.20939}{e}^{-{\left(\frac{t}{61.6529}\right)}^{3.20939}}$

累计分布函数为： 

$F\left(t\right)=1-e^{{-\left(\frac{t}{61.6529}\right)}^{3.20939}}$

可靠度函数为： 

$R\left(t\right)=e^{-{\left(\frac{t}{61.6529}\right)}^{3.20939}}$

失效率函数为： 

$\mathrm{\λ}\left(t\right)=\frac{3.20939}{61.6529}{\left(\frac{t}{61.6529}\right)}^{2.20939}$

(五)驱动机构可靠性特征值 

驱动机构服从三参数Weibull分布模型，其形状参数β为1.37215，尺度参数η为11.6081，阈值γ为30.6272，故可得故障时间概率密度函数f(t)、累计分布函数F(t)、可靠度函数R(t)和失效率函数λ(t)如下所示。 

故障时间概率密度函数为： 

$f\left(t\right)=\frac{1.37215}{{11.6081}^{1.37215}}{(t-30.6272)}^{0.37215}{e}^{-{\left(\frac{t-30.6272}{11.6081}\right)}^{1.37215}}$

累计分布函数为： 

$F\left(t\right)=1-e^{-{\left(\frac{t-30.6272}{11.6081}\right)}^{1.37215}}$

可靠度函数为： 

$R\left(t\right)=e^{-{\left(\frac{t-0.6272}{11.6081}\right)}^{1.37215}}$

失效率函数为： 

$\mathrm{\λ}\left(t\right)=\frac{1.37215}{11.6081}{\left(\frac{t-30.6272}{11.6081}\right)}^{0.37215}$

(六)基础制动装置可靠性特征值 

基础制动装置服从三参数Weibull分布模型，其形状参数β为1.11546，尺度参数η为36.7574，阈值γ为35.9950，故可得故障时间概率密度函数f(t)、累计分布函数F(t)、可靠度函数R(t)和失效率函数λ(t)如下所示。 

故障时间概率密度函数为： 

$f\left(t\right)=\frac{1.11546}{{36.7574}^{1.11546}}{(t-35.9950)}^{0.11546}{e}^{-{\left(\frac{t-35.9950}{36.7574}\right)}^{1.11546}}$

累计分布函数为： 

$F\left(t\right)=1-e^{-{\left(\frac{t-35.9950}{36.7574}\right)}^{1.11546}}$

可靠度函数为： 

$R\left(t\right)=e^{-{\left(\frac{t-35.9950}{36.7574}\right)}^{1.11546}}$

失效率函数为： 

$\mathrm{\λ}\left(t\right)=\frac{1.11546}{36.7574}{\left(\frac{t-35.9950}{36.7574}\right)}^{0.11546}$

二、采用进化算法优化的神经网络模型对转向架各子系统的故障率进行计算。 

PSO-BP预测模型的建立包括以下四步，建模流程图如图2所示。 

(1)获取实验数据样本，对数据样本进行筛选处理，确定训练数据样本和测试数据样本； 

(2)利用PSO算法优化确定BP相关参数，确定PSO的适应度函数，以及粒子群的规模、迭代代数等。 

(3)将PSO训练得到的最优参数以及训练数据样本作为输入进行BP模型训练。 

(4)将测试数据输入至训练得到的预测模型，检验模型的预测精度，最后得到PSO-BP预测模型。 

各关键子系统故障率预测是指利用系统的历史故障率预测将来的故障率。根据地铁车辆故障数据的特性，以“周”为时间单位，统计各关键子系统的故障次数，统计130条故障次数据，可通过下列公式计算得到各关键系统的故障率。 

${\mathrm{\λ}}_i=\frac{{\mathrm{\Σ}}_{i=1}^n{r}_i}{{\mathrm{\Σ}}_{i=1}^n{t}_i}=N_0/T$

n——抽样台数 

t_i——检测周期内第i台设备的实际工作时间，单位为小时 

r_i——检测周期内第i台设备出现的故障次数 

N₀——检测周期内系统累计故障次数 

T——检测周期内总工作时间。 

建立的基于PSO的BP神经网络模型，即采用PSO优化BP神经网络连接权值，并对BP神经网络进行训练。在建模过程中，需确定PSO的一些参数，设置如下： 

(1)粒子数目 

粒子数目一般取值20～40。对于大多数问题来说，30个粒子就可以取得很好的结果，本节选取粒子数目为30。 

(2)粒子最大速率 

粒子最大速率决定粒子在一次飞行中可以移动的最大距离。必须限制粒子最大速率，否则，粒子就可能跑出搜索空间。粒子的最大速率通常设定为粒子范围的宽度。本节选取为0.5。 

(3)加速常数 

学习因子c₁和c₂表示粒子受社会知识和个体认知的影响程度，通常设为相同值以给两者同样的权重c₁＝c₂＝2。在找到这两个最优位置时，粒子根据如下的公式(4.33)和(4.34)来更新自己的速度和位置。 

$v_i^{k+1}=w*v_i^k+c_1{r}_1(p_i^k-x_i^k)+c_2{r}_2(p_g^k-x_i^k)---(4-31)$

$x_i^{k+1}=x_i^k+v_i^{k+1}---(4-32)$

式中，表示第i个粒子在第k+1代时的飞行速度；表示第i个粒子在第k+1代时的位置；表示第i个粒子到第k代的最优位置；表示种群到第k代的最优位置；为个体认知；为种群认知；w为惯性权重；是粒子的速度；c₁和c₂为学习因子；r₁和r₂为均匀分布在[0,1]内的随机数；i＝1,2…N。 

(4)适应度函数 

适应度函数用于区分种群中个体好坏的标准，也是进行选择的唯一依据，其值越大越好。根据本文所确定的网络评价指标，此处采用网络输出与期望输出间的均方根误差(RMSE)的倒数为评价标准。其适应度函数为 

f＝1/RMSE          (4-33) 

(5)初始化。将输入输出数据归一化到[-1,1]区间，初始化产生的粒子种群代表不同权重与阈值组合的神经网络，产生初始BP网络结构。 

(6)评价。依据公式(4-33)计算其适应度。 

(7)更新位置和速度。通过比较适应度不断更新粒子的位置，适应度最优的个体极值即为全局极值，对应的权值与阈值为粒子种群的当前最优解。依据公式(4-31)和公式(4-32)更新速度。 

(8)算法停止的判断。达到最大迭代次数或适应度达到期望误差，停止迭代。 

(9)最优解的生存。迭代停止时，全局权值对应的权值与阈值为训练样本的最优解。 

(10)将最优解带入BP神经网络模型中学习，用来对故障率进行预测。 

三、基于成本优化比例风险模型的转向架视情维修决策。 

在比例风险模型中，假定系统的故障次数服从非齐次Possion分布(NHPP)，则在给定时间t内，系统故障的次数为N(t)，系统故障的风险函数表达式如式(5-1)所示： 

$\mathrm{\υ}\left(t\right)=\underset{\mathrm{\Δt}\→\∞}{\lim }\frac{P\{N(t+\mathrm{\Δt})-N\left(t\right)\≥1|Z\left(t\right)\}}{\mathrm{\Δt}}=\frac{\mathrm{dE}\left(N\left(t\right)\right)}{\mathrm{dt}}---(5-1)$

其中，Z(t)是系统运行到时间t时协变量的检测量。 

对于非齐次Possion分布，有 

$P\{N(t+\mathrm{\Δt})-N\left(t\right)\≥n|Z\left(t\right)\}=\frac{{\left[H_t\left(\mathrm{\Δt}\right)\right]}^n}{n!}\exp (-H_t\left(\mathrm{\Δt}\right))---(5-2)$

其中,H_t(Δt)是在时间[t,t+Δt]内系统的平均故障次数。 

NHPP的基准风险函数主要有指数分布、Gumbel分布、Weibull分布三种，如式(5-3)至式(5-5)所示。 

υ₀＝α           (5-3) 

υ₀＝αexp(γt)        (5-4) 

υ₀＝αγt^γ-1       (5-5) 

故障过程中PHM模型的基准风险函数可为以上三种分布中的任何一种，PHM模型为： 

υ(t|X(t))＝υ₀exp(γ^TX(t))      (5-6) 

而系统在t时刻之后的可靠性函数和t时刻之后失效的概率密度函数如式(5-7)和式(5-8)所示： 

$R\left(t\right|X\left(t\right))=P\{N(t+\mathrm{\Δt})-N\left(t\right)=0|X\left(t\right)\}=\exp \{-{\∫}_t^{t+\mathrm{\Δt}}\mathrm{\υ}\left(u\right|X\left(u\right))\mathrm{du}\}---(5-7)$

$f\left(t\right|X\left(t\right))=-R_t^{\′}\left(t\right)={\mathrm{\υ}}_0\exp \left({\mathrm{\γ}}^{T}X\left(t\right)\right)\exp \{-{\∫}_t^{t+\mathrm{\Δt}}\mathrm{\υ}\left(u\right|X\left(u\right))\mathrm{du}\}---(5-8)$

全参数PHM的求解步骤： 

全参数PHM是基准风险函数已知的情形(服从某种分布)，对基准函数的形状参数、尺度参数和协变量的回归系数进行估计，可以使这些参数的估计达到较高的精度。 

在工程实践中，Weibull分布函数被广泛运用于全参数PHM基准风险函数，简称Weibull PHM。Weibull分布由瑞典物理学家Waloddi Weibull在研究材料裂纹中发现，之后被学者们广泛地应用在寿命数据统计和可靠性研究上。 

在Weibull分布中，形状参数、尺度参数分别表征故障率的变化趋势、坐标尺度的放大或缩小。 

具体Weibull PHM的建模过程主要分为四个步骤，步骤如下所示^[62]。 

(1)-2LL检验 

-2LL检验是检验样本是否符合Weibull分布最常用的方法。 

Weibull分布的可靠度函数为： 

$F\left(t\right)=1-\exp (-{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}})---(5-9)$

故障分布函数为： 

$R\left(t\right)=1-F\left(t\right)=\exp (-{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}})---(5-10)$

通过对等式(5-10)两边取两次对数可得式(5-11)，如下所示： 

$\ln (-\ln \left(R\left(t\right)\right))=\ln \ln \left(\frac{1}{1-F\left(t\right)}\right)=\mathrm{\β}\ln t-\mathrm{\β}\ln \mathrm{\η}---(5-11)$

由式(5-11)可知，若样本服从Weibull分布，根据lnt与ln(-ln(R(t)))之间的函数拟合的直线斜率应大于1，直线的斜率是Weibull分布形状参数的估计值，-2LL检验法同时也是观察样本数据是否符合Weibull分布最直观的方法。 

(2)构造似然函数 

Weibull PHM为： 

$\mathrm{\υ}\left(t\right|X\left(t\right))=\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \left({\mathrm{\γ}}^{T}X\left(t\right)\right)---(5-12)$

系统在t时刻之后的可靠性函数和在t时刻系统故障的概率密度函数分别为 

$R\left(t\right|X\left(t\right))=P\{N(t+\mathrm{\Δt})-N\left(t\right)=0|X\left(t\right)\}=\exp \{-{\∫}_t^{t+\mathrm{\Δt}}\mathrm{\υ}\left(u\right|X\left(u\right))\mathrm{du}\}---(5-13)$

$f_t\left(t\right|X\left(t\right))=\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \left({\mathrm{\γ}}^{T}X\left(t\right)\right)\exp \{-{\∫}_t^{t+\mathrm{\Δt}}\mathrm{\υ}\left(u\right|X\left(u\right))\mathrm{du}\}---(5-14)$

似然函数为： 

$L(\mathrm{\β},\mathrm{\η},\mathrm{\γ})=\underset{i=1}{\overset{n}{\mathrm{\Π}}}f{\left(t_i\right|\mathrm{\β},\mathrm{\η},\mathrm{\γ})}^{{\mathrm{\δ}}_i}R{\left(t_i\right|\mathrm{\β},\mathrm{\η},\mathrm{\γ})}^{1-{\mathrm{\δ}}_i}---(5-15)$

Weibull PHM的似然函数为： 

$L(\mathrm{\β},\mathrm{\η},\mathrm{\γ})=\left\{\underset{i\∈D}{\mathrm{\Π}}\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t_i}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp \left({\mathrm{\γ}}_k{X}_k\left(t_i\right)\right)\right\}\left\{\underset{j=1}{\overset{n}{\mathrm{\Π}}}\exp \right[-H\left(t_j\right)\left]\right\}---(5-16)$

其中： 

$H\left(t_j\right)={\∫}_0^{T_j}\exp \left({\mathrm{\γ}}^{T}X\left(t\right)\right)d{\left(\frac{t_j}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}---(5-17)$

对式(5-16)两边取对数，得到Weibull PHM的对数似然函数，如式(5-18)所示。 

$\ln L(\mathrm{\β},\mathrm{\η},\mathrm{\γ})=\ln {\left(\frac{t}{\mathrm{\η}}\right)}^d+\underset{i=1}{\overset{d}{\mathrm{\Σ}}}\ln {\left(\frac{t_i}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}+\underset{i=1}{\overset{d}{\mathrm{\Σ}}}\ln \left({\mathrm{\γ}}_k{X}_k\left(t_i\right)\right)-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}H\left(t_j\right)---(5-18)$

(3)对各参数极大似然估计的求解 

极大似然估计求解过程如下。 

设未知参数所组成的向量为： 

θ＝[θ₁,θ₂,...,θ_p+3]＝[β,η,γ₁,γ₂,...,γ_p]     (5-19) 

对各参数的一阶偏导数组成的向量记为G(θ)＝g_i(θ)，其中： 

$g_i\left(\mathrm{\θ}\right)=\frac{\∂l}{\∂{\mathrm{\θ}}_i},i=\mathrm{1,2},...,p+3---(5-20)$

对数似然函数分别对各待估参数求一阶偏导数,并令其等于0，得： 

$\frac{\∂l}{\∂\mathrm{\β}}=\frac{d}{\mathrm{\β}}+\underset{i=1}{\overset{d}{\mathrm{\Σ}}}\ln \left(\frac{t_i}{\mathrm{\η}}\right)-\ln \left(\frac{t_i}{\mathrm{\η}}\right)\underset{j=1}{\overset{d}{\mathrm{\Σ}}}H\left(t_j\right)=0---(5-21)$

$\frac{\∂l}{\∂\mathrm{\η}}=-\frac{d}{\mathrm{\η}}-\frac{(\mathrm{\β}-1)d}{\mathrm{\η}}-\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\frac{\∂H\left(t_j\right)}{\∂\mathrm{\η}}=-\frac{\mathrm{d\β}}{\mathrm{\η}}-\frac{\mathrm{\β}}{\mathrm{\η}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}H\left(t_j\right)=0---(5-22)$

$\frac{\∂l}{\∂{\mathrm{\γ}}_k}=x_k\left(t_i\right)-x_k\left(t_j\right)\underset{j=1}{\overset{n}{\mathrm{\Σ}}}H\left(t_j\right)=0,k=\mathrm{1,2},...,p+1---(5-23)$

通过联立式(5-21)至式(5-23)求得各未知参数的估计值。 

(4)得出极大似然估计结果 

最终求解得到各待估参数的估计值从而得到： 

通过最大似然那估计求得Weibull分布的形状参数、尺度参数和协变量的回归系数，可以得出基准风险函数，从而进一步求得转向架系统的可靠度函数、累计故障分布函数和故障概率密度函数等。 

最小维修成本法的计算过程如下： 

单位时间平均维修成本为： 

$C\left(t\right)=\frac{C_r[(C_p/C_r)+E\left(N\left(t\right|X\left(t\right))\right)]}{t}---(5-28)$

其中，C_r是故障后更换送修的平均成本，C_p是预防维修的平均成本，C_p/C_r是预防维修成本与更换送修成本的比率。 

E(N(t|X(t)))是检修间隔期内预期的故障次数，且有： 

$E\left(N\left(t\right|X\left(t\right))\right)=H\left(t\right|X\left(t\right))={\∫}_0^t\exp (\mathrm{\γ}\·X\left(t\right))d{(u/\mathrm{\η})}^{\mathrm{\β}}---(5-29)$

当成本最小时，故障率函数斜率应最小。按式(5-29)绘出H(t)-t曲线，过点(0，-C_p/C_r)作函数的切线，切点处对应的是成本最小的优化维修间隔时间。 

视情维修决策控制限的确定，根据所求得的优化时间，从故障数据中找出与此最接近的失效时间与伴随变量数据，由此求得故障率阈值h^*，即： 

$h^{*}\≥\frac{\mathrm{\β}}{\mathrm{\η}}{\left(\frac{t}{\mathrm{\η}}\right)}^{\mathrm{\β}-1}\exp (\mathrm{\γ}\·X\left(t\right))---(5-30)$

对式(5-30)两边取对数，得： 

$\mathrm{\γ}\·X\left(t\right)\≤\ln \left(\frac{{\mathrm{\η}}^{\mathrm{\β}}{h}^{*}}{\mathrm{\β}}\right)-(\mathrm{\β}-1)\ln t---(5-31)$

以工作时间t为横轴、以为纵轴作图，所得结果即为维修阈值对应的控制限。然后根据当前运行时刻监控到的相关伴随变量计算出相应γ·X(t)，并在图中画点。若该点处于曲线下方，则处于正常工作区，可继续工作；若处于曲线上方，则处于性能衰退区，需要采取维修措施。若恰好落在上下曲线之间的阴影区，则处于临界区，应对于该状态加强监控跟踪，考虑视情维修。 

分析可修系统时，随着运行时间的累计，可修系统的状态逐渐劣化，如图3所示。设检测时刻为{T_i；i＝0,1,…,k},两相邻检测时刻的状态量的差值记为ΔX_(k,k+1)。对于复杂系统的监控及视情维修决策过程中，通常至少会设定状态阈值的上下限，其中上限L是失效阈值。在运行过程中，一旦发现系统状态值超过L，则系统此时为不可用状态，按规定必须进行修复性修理或更换后才重新投入使用；下限M是预防性维修或更换阈值，代表系统的潜在故障开始出现。如果系统的状态值超过了M，则应考虑对部件作相应的排故或预防性维修工作。 

不同类型的转向架性能指标是不同的，如地铁转向架使用同一类型的转向架，则可假定该类型转向架的故障数据符合同一分布，作为建模时的样本数据。以表5-1的32例维修时监控的故障数据为例，采用比例风险模型进行视情维修决策建模，其中模型协变量选取的是转向架各子系统的故障率和安全运营天数，模型的输出为转向架系统故障间隔时间。 

表5-1转向架历史故障数据的部分样本 

注：GJF、LDF、QDF、THF、ZDF、LJF分别是构架、轮对和轴箱、驱动机构、弹簧装置、基础制动装置、连接装置的故障率，GJT、LDT、QDT、THT、ZDT、LJT分别是构架、轮对和轴箱、驱动机构、弹簧装置、基础制动装置、连接装置在转向架系统维修时的安全运营天数。 

(一)协变量的筛选 

影响转向架系统性能劣化的协变量有很多，在对转向架的建模时，以转向架的故障间隔时间为自变量，主要考虑的协变量包括构架、轮对和轴箱、驱动机构、弹簧装置、基础制动装置、连接装置在转向架系统维修时的故障率和安全运营天数等。这些协变量经相关性分析，从中筛选出最显著的协变量纳入所建立的比例风险模型。 

(二)比例风险模型的建立 

首先检验寿命周期变量是否服从Weibull分布。针对故障间隔时间样本，利用上述检验方法进行计算，利用Minitab得出计算结果，通过-2LL检验法可知，故障间隔时间t与失效分布函数F(t)之间的对数曲线图近似拟合为一条直线，故障间隔时间符合Weibull分布。从调整的A-D检验量和P值来看，A-D统计量为0.367，P值大于0.250，表明故障间隔时间服从Weibull分布规律。并同时得出，Weibull分布的形状参数β为1.75488，尺度参数η为47.76002。 

在转向架的比例风险模型中，相应的风险函数为： 

$h\left(t\right|X\left(t\right))=\frac{1.75488}{47.76002}{\left(\frac{t}{47.76002}\right)}^{0.75488}\exp (\mathrm{\γ}\·X\left(t\right))---(5-32)$

其中，γ为与协变量向量对应的各回归系数组成的向量，γ和β、η一样都是待估计的参数。 

按照前述对比例风险模型的建模步骤，根据转向架失效数据样本，通过极大似然估计法求解得到各个待估参数。 

基于极大似然估计原理，首先运用SPSS软件对各协变量间的相关系数进行比较，各协变量间的相关系数见表。从表中可知，相关系数较大则说明变量间明显相关，通常在建模时，相关系数较大的协变量须剔除，相关系数较小的协变量被保留，从而使各变量之间具有较高的独立性。 

表5-2协变量之间的相关系数 

通过逐一比较，剔除明显相关的协变量后，实际用于建模的协变量主要包括：构架故障率、驱动机构故障率、连接装置的故障率和驱动、弹簧装置、基础制动装置的安全运营天数。以此作为初步建模时待纳入的协变量，然后采用协变量筛选的方法进行显著性协变量的筛选，最终确定出能够反映系统失效的最显著协变量。 

比例风险模型显著性协变量筛选结果显示，轮对和轴箱、弹簧装置、连接装置的故障率和构架、轮对和轴箱、驱动机构、弹簧装置、连接装置的安全运营天数等因无显著作用而被剔除，最终构架故障率、驱动机构故障率、连接装置的故障率和基础制动装置的安全运营天数等四个协变量被筛选进入模型，记为X(t)＝{x₁(t),x₂(t),x₃(t),x₄(t)}。 

Weibull分布的形状参数β为1.75488，尺度参数η为47.76002，可得相应的系统累积风险函数曲线，如图4所示。 

$F\left(t\right)=1-e^{-{\left(\frac{t}{47.76002}\right)}^{1.75488}---(5-33)}$

因没有搜集到确定的相关数据，先按一般情况假定预防和送修两类维修成本 的比率C_p/C_r为0.2^[66]，在图4中过(0，-0.2)向曲线画切线，求得的切点位置对应的故障时间为27天左右，即为成本最小的优化维修间隔时间。 

经极大似然估计得到Weibull分布的形状参数和尺度参数，并结合协变量筛选得出的回归系数，最终得到系统的风险函数为： 

$v\left(t\right|X\left(t\right))=\frac{1.75488}{47.76002}{\left(\frac{t}{47.76002}\right)}^{0.75488}\exp [-2.889x_1\left(t\right)-1.091x_2\left(t\right)+3.483x_3\left(t\right)-0.028x_4\left(t\right)]---(5-34)$

其中，x₁，x₂，x₃，x₄分别对应构架故障率、驱动机构故障率、连接装置的故障率、基础制动装置的安全运营天数。 

计算求得成本最小的优化维修间隔时间为27天，从故障数据中找出与此最接近的失效时间，分别为24天、29天，由此求得故障率阈值h^*分别为0.248、0.340，通过式(5-31)计算出视情维修决策的上下限，如图5所示，图中横坐标为故障时间，纵坐标为综合状态量，综合状态量可通过γ·X(t)计算。 

实际决策中，图5中的上线为视情维修的上限，如综合状态量超过上线，则应该立即视情维修；下线为视情维修的下限，如综合状态量处于上下限之间，则系统处于临近失效的过渡过程，需要加强对系统的监控或者对系统进行预防性维修；如综合状态量处于下限以下，一般不用考虑视情维修。 

为了验证模型的有效性，将后15组数据分别预测的构架、驱动机构、连接装置故障率、基础制动装置安全运营天数作为输入，PHM模型分析得出的维修建议与实际维修情况对比，如表5-3所示。 

表5-3PHM输出与实际维修情况对比 

其中，x₁，x₂，x₃，x₄分别是PSO-BP预测的构架故障率、驱动机构故障率、 连接装置故障率、基础制动装置的安全运营天数。 

结合图6可知，在这15组数据中，在时间间隔为56天处，综合状态量值为5.19，图中显示，超过维修决策的上限，则说明此时需要对转向架系统维修，时间间隔为68天时，综合状态量为1.96，在维修决策上限和下限之间，则说明此时需要对转向架系统进行监控。 

Claims (4)

1.一种城轨列车转向架的故障预测与视情维修方法，所述转向架包括构架、弹簧装置、连接装置、轮对和轴箱、驱动机构、基础制动装置六个子系统，其特征在于，该方法依次包括如下步骤：

（1）根据采集的历史故障数据进行删失处理，基于生存分析的方法确定转向架各子系统的分布模型，并得出各子系统的可靠性特征函数，计算出各子系统的可靠度，将各子系统中可靠度最低的子系统确定为转向架最薄弱环节；

（2）采用进化算法优化的神经网络模型计算转向架各子系统的故障率；

（3）将转向架各子系统的安全运营天数和计算得到的故障率作为协变量进行比例风险模型建模，得出转向架系统视情维修的阈值和控制限，控制限上限是失效阈值,在运行过程中，一旦发现系统状态值超过上限，则系统此时为不可用状态，按规定必须进行修复性修理或更换后才重新投入使用；控制限下限是预防性维修或更换阈值，代表系统的潜在故障开始出现，如果系统的状态值超过了下限，应对部件作相应的排故或预防性维修工作，如果系统的状态值低于下限，则不用考虑维修。

2.根据如权利要求1所述的城轨列车转向架的故障预测与视情维修方法，其特征在于，

步骤（1）包括如下步骤：

1）建立轮对和轴箱、弹簧装置、连接装置两参数威布尔分布模型，

故障分布函数为：

可靠度函数为：

                           

概率密度函数为：

                   

失效率函数为：

其中，，且，t为故障间隔时间，和分别是分布的形状参数和尺度参数；

2）建立构架、驱动机构、基础制动装置三参数威布尔分布模型，

故障分布函数为：

                           

可靠度函数为：

                             

其概率密度函数为：

                    

失效率函数为：

其中，t为故障间隔时间，且，和分别是分布的形状参数和尺度参数,是分布的位置参数。

3.根据前述权利要求所述的城轨列车转向架的故障预测与视情维修方法，其特征在于，步骤（2）包括如下步骤：建立的基于PSO的BP神经网络模型，采用PSO优化BP神经网络连接权值，并对BP神经网络进行训练，在建模过程中，确定PSO的参数，设置如下：

1)选取粒子数目为30；

2)粒子最大速率选取为0.5；

3)粒子根据如下的公式来更新自己的速度和位置，

          

                              

式中，表示第i个粒子在第k+1代时的飞行速度；表示第i个粒子在第k+1代时的位置；表示第i个粒子到第k代的最优位置；表示种群到第k代的最优位置；为个体认知；为种群认知；w为惯性权重；是粒子的速度；和为学习因子；和为均匀分布在[0,1]内的随机数；i=1, 2…N；

4)适应度函数为

                              ；

5)初始化，将输入输出数据归一化到[-1,1]区间，初始化产生的粒子种群代表不同权重与阈值组合的神经网络，产生初始BP网络结构；

6)评价，依据4）中公式计算其适应度；

7)更新位置和速度，通过比较适应度不断更新粒子的位置，适应度最优的个体极值即为全局极值，对应的权值与阈值为粒子种群的当前最优解，更新速度；

8)算法停止的判断，达到最大迭代次数或适应度达到期望误差，停止迭代；

9)最优解的生存，迭代停止时，全局权值对应的权值与阈值为训练样本的最优解；

10)将最优解带入BP神经网络模型中学习，用来对故障率进行预测。

4.根据前述权利要求所述的城轨列车转向架的故障预测与视情维修方法，其特征在于，步骤（3）包括如下步骤：

1)建立转向架的比例风险模型，其风险函数为：

求得故障率阈值：

                      

其中，t为故障间隔时间，且，和分别是分布的形状参数和尺度参数,是分布的位置参数，X(t)为协变量；

2) 以工作时间t为横坐标、以的值为纵坐标作图，图形中上线为故障率阈值对应的控制限的上限，下线为故障率阈值对应的控制限的下限。