一种新的威布尔比例失效模型参数估计方法

张路朋，王 娟，王 攀，陈 亮

(1.国网河北省电力公司邯郸供电分公司，河北 邯郸 056000；2.国网河北省电力公司石家庄供电分公司，石家庄 050051)

2017-06-05

张路朋(1988-)，男，工程师，主要从事电力系统继电保护工作。

一种新的威布尔比例失效模型参数估计方法

张路朋1，王 娟1，王 攀1，陈 亮2

(1.国网河北省电力公司邯郸供电分公司，河北 邯郸 056000；2.国网河北省电力公司石家庄供电分公司，石家庄 050051)

针对比例失效模型在参数估计中存在计算速度慢、收敛效果差以及初始数据要求高等缺陷，提出了基于极大似然估计方法的分步估计算法，首先对形状参数和尺度参数进行极大似然估计，然后将计算结果代入比例失效模型，再对状态指标回归系数进行极大似然估计，经算例分析表明，该算法在计算速度、计算时间、数据样本选择方面具有较大的优势，证实了所提算法的有效性。

比例失效模型；参数估计；极大似然估计；分步估计算法

比例失效模型(Proportional Hazards Model，PHM)由Cox首先提出，最初用于医学研究，随着设备可靠性技术的不断发展，逐步成为设备状态评估与状态维修的典型模型，PHM能有效地将设备状态信息用于设备的状态维修，从而在考虑经济性和可靠性的前提下获得设备的最大可用寿命，并进行预知性维修或者更换决策。

目前应用最为广泛的是威布尔比例失效模型(Weibull Proportional Hazards Model, WPHM)，文献[1-2]将WPHM应用于风电机组的状态检修中，通过采集风机的温度和振动等数据，制定维修阈值，通过最小成本法计算出最佳维修时机，取得了较好的效果。

WPHM在应用过程中，参数的估算始终是一大难题，目前应用最多的方法是极大似然估计算法，文献[3-5]利用极大似然估计算法，将WPHM各参数进行整体估计，在参数估计过程中，采取的做法是利用历史故障数据，将形状参数、尺度参数、状态指标回归系数同时计算，这一做法大大增加了参数估计的计算量。但这种算法对初始数据的要求较高，同时可能会出现结果不收敛的情况。为了避免这种情况的发生，对WPHM的参数估计算法进行了改进，提出了分步估计算法。所谓分步估计算法，即将WPHM形状参数、尺度参数、状态指标回归系数分开估计，这种方法收敛效果更好，计算速度更快，且对初始数据的要求较低。

1 威布尔比例失效模型

威布尔比例失效模型兼顾设备自然老化和状态协变量冲击两方面的影响，在设备状态检修中得到广泛应用，如式(1)。

(1)

式中：h0(t)为基本失效函数，由运行时间t唯一决定，在威布尔比例失效模型中，用威布尔故障率函数来表示，反应变压器随时间的自然老化[6]；η是寿命参数；β>1是形状参数；φ(Z)为状态协变量函数，由状态协变量集合Z唯一决定，反应设备受状态协变量集合Z的冲击，设状态协变量集合Z包含p个状态协变量，分别用Z1～Zp表示，其回归系数分别为γ1～γp。

提出的分步估计算法就是通过极大似然估计算法，对WPHM的参数η，β与[γ1，γ2，…，γp]分开估计。

首先对式(1)所示的威布尔比例失效模型进行分析发现：

a. WPHM的基本失效函数h0(t)为威布尔分布，描述的是设备在理想条件下的自然劣化过程，其分布只取决与运行时间t。

b. 状态指标函数φ(Z)描述的是设备在实际工况下受运行环境影响的程度，其分布只取决于选取的状态指标集Z=[Z1，Z2，…，Zp]，即 φ(Z)受实际工况的影响。

c.h0(t)和φ(Z)在数学上没有联系。

以上结论进一步表明，可将基本失效函数h0(t)中涉及到的形状参数β和尺度参数η与状态指标函数φ(Z)中涉及到的γ=[γ1，γ2，…，γp]分开估计。采取先求取β和η，从而得到基本失效函数h0(t)，然后将其作为已知量代入式(1)，对γ=[γ1，γ2，…，γp]进行计算，最终得到威布尔比例失效模型。

在可靠性模型中，进行参数估计所用到的数据主要有故障数据和结尾数据，在参数估计过程中，暂不考虑结尾数据，只采用故障数据。参数求取分为β，η的求取和γ=[γ1，γ2，…，γp]求取，采用的算法是极大似然估计算法。

2 参数估计方法

2.1 形状参数和尺度参数的求取

基本失效函数用威布尔失效率函数表示，表达式为

(2)

β，η的求取采取的是连续函数的极大似然估计方法，数理统计原理中强调，连续函数的极大似然估计采用概率密度函数，所以首先要求解威布尔分布的概率密度函数。根据可靠性知识[7-8]，可推出威布尔分布的概率密度函数为

(3)

基本失效函数h0(t)只取决于运行时间，所以对其估计时，只需要统计设备的历史故障时间，然后将其从小到大排列，得到历史故障时间集合T={t1，t2，…，tp}， 的极大似然估计算法如下。

为了便于计算，令m=ηβ，得到似然函数如式(4)所示：

(4)

对式(4)两端同时取对数，并分别求β,m的偏导，代入m=η-β，令其等于零可得：

(5)

采用牛顿法对上式进行收敛计算，可以求解出β，η的具体值。

2.2 状态协变量参数的估计

求取γ=[γ1，γ2，…，γp]，不但需要考虑故障时的运行工况，还需要考虑正常时的运行工况。有别于β，η的求取，γ=[γ1，γ2，…，γp]的求取采用的是散点极大似然估计，因此，不需要先求取概率密度函数。采用固定间隔时间采集的状态指标检测值来对γ=[γ1，γ2，…，γp]进行散点估计。设在ti时刻设备的状态监测值为Z(ti)，从开始投入运行到故障截止，一共采集了n个样本点数据，则似然函数可以写成如式(6)的形式：

(6)

(7)

对上式两边同时取对数可得

(8)

上式中含有积分项，对γ求偏导较难，考虑到在连续监测中，部件的状态协变量Z(ti)在较小的一段时间内变化幅度不大，可以把{Z(t)|t≥0}看作是右连续的阶跃过程，于是采用把积分项写成分段积分和的形式。假设0

(9)

分部积分可得

(10)

定义ev=exp[γ·Z(tv)](v=0,2,…,j)，则式(10)的积分项运用分部积分求解，可进一步推出

(11)

带入式(8)可得

(12)

式(12)对γ=[γ1，γ2，…，γp]中的γk求偏导并令其等于零，由于β，η此时都为常数，于是可得

(13)

令上式等于零，采用牛顿法对上式进行收敛计算，可以求解γ=[γ1，γ2，…，γp]的具体值。

各参数求出之后，代入式(1)可以得到威布尔比例失效模型。

3 算例分析

变压器长期运行于过负荷状态，内部会发生油气变化，此时需要时刻监测变压器的运行状态，以便制定合理的检修策略。利用油色谱分析技术可以得到变压器油中溶解的氢气(H2)、一氧化碳(CO)、二氧化碳(CO2)、甲烷(CH4)、乙烯(C2H4)、乙烷(C2H6)、乙炔(C2H2)等气体体积分数[9-10]。在进行算例分析时，利用变压器的油色谱分析数据来计算变压器的威布尔比例失效模型，可为制定变压器检修策略提供依据。

对某型号主变压器的历史故障数据和油色谱监测数据进行统计，表1为该型号主变压器的8次历史故障时间，表2为该型号第8次维修周期中的油色谱分析数据。运用MATLAB软件进行计算。利用表1的历史故障时间并结合2.1节方法，可对形状参数和尺度参数进行计算，计算结果如表3所示，利用表2的数据结合2.2节方法，可对该型号变压器的以及H2、CO、CO2、CH4、C2H4、C2H6、C2H2的回归系数进行计算，计算结果如表3所示。利用分步计算法进行计算时，计算参数γ=[γ1，γ2，…，γp]用时2 s，计算 用时14 s，整个计算过程用时16 s。

表1 某型号变压器的历史故障数据部分样本

主变压器序号12345678故障天数/d856106895514281222113113971646

表2 某型号变压器第8次维修周期内油色谱监测数据部分样本 μL/L

采样时间φ(H2)φ(CO)φ(CO2)φ(CH4)φ(C2H4)φ(C2H6)φ(C2H2)2016-01-1510881563.00.30.90.42016-07-16161002893.10.41.03.72016-10-15291192223.50.61.04.32016-11-18311283974.20.81.14.82016-12-16391372114.51.11.06.42017-01-15441383285.31.01.77.2

表3 参数计算值(分步计算法)

参数名称βηγ1γ2γ3γ4γ5γ6γ7参数值9.0372853.30.13420.09260.07430.16320.21530.17200.1682计算用时/s 2 14

利用整体估计方法计算结果如表4所示，计算用时为37 s。

表4 参数计算值(整体计算法)

参数名称βηγ1γ2γ3γ4γ5γ6γ7参数值9.0423853.50.13380.09270.07490.16520.21430.17260.1692计算用时/s37

表4与表3对比发现，利用提出的分步估计算法求解结果与整体估计算法结果相差较小，可满足工程误差要求，但是分步计算法总体用时16 s，而整体计算方法用时37 s，是分步计算方法的两倍多，可见分步估计算法在计算速度、计算时间方面具有较大的优势，而这种优势随着数据量的增加会逐步扩大。同时在数据样本选择方面，分步计算法对数据要求不高，计算更容易收敛，而整体计算方法对数据要求较高，计算结果容易发散。

将表3计算结果代入式(1)中，可以得到变压器的威布尔比例失效函数，可根据该失效函数制定相应的状态检修策略，为主变运行维护人员提供检修建议。

4 结束语

对传统比例失效模型参数估计方法进行了改进，提出一种基于极大似然估计算法的新的参数估计方法，即分步估计算法，该算法认为威布尔比例失效模型的参数与在数学上互不相关，可以利用极大似然估计算法进行分开计算。算例分析表明，分步估计算法在计算速度、计算时间、数据样本选择方面具有较大的优势，证实了所提方法的有效性。

[1] Ardine，Banjevicd．Optimization a mine haul truck wheel motors' condition monitoring program use for proportional hazards modeling [J]．Journa l of Quality in Maintenance Engineering，2001，7(4)：286 - 301．

[2] Xiang Wu，Sarah M. Ryan．Optimal Replacement in the Proportional Hazards Model With Semi-Markovian Covariate Process and Continuous Monitoring[J]． IEEE Tanscctions on Reliability，2011，60(3)：580-589．

[3] 张继全，赵洪山．基于比利失效模型的风机齿轮箱轴承检修[J]．电源技术，2011，35(5)： 570-573．

[4] 陈寒雨．基于设备故障率评估的视情维修研究[D]．成都：电子科技大学，2011．

[5] 马维青，于瑶章．基于比例失效模型的配电网设备状态检测[J]．信息技术，2015，(15)：206-209．

[6] 赵洪山，张路朋．基于可靠度的风电机组预防性机会维修策略[J]．中国电机工程学报，2014，34(22)：3777-3783．

[7] 张路朋．风电机组的状态机会维修策略研究[D]．保定：华北电力大学，2015．

[8] 张路朋，赵洪山．基于时间延迟的风机齿轮箱状态优化维修[J]．中国电力，2014，47(11)：108-111．

[9] DL/T 722-2014．变压器油中溶解气体分析和判断导则[S]．

[10] 李守学，司昌健，张春丰，等．某500 kV变压器油色谱异常监测及缺陷诊断[J]．吉林电力，2017，45(1)：44-46．

A New Method of Parameter Estimation for Proportional Hazards Model

Zhang Lupeng1，Wang Juan1,Wang Pan1,Chen Liang2

(1.State Grid Hebei Electric Power Corporation,Handan Power Supply Branch,Handan 056000,China；2.State Grid Hebei Electric Power Company Shijiazhuang Power Supply Branch,Shijiazhuang 050051,China)

Aiming at the slow computation speed,poor convergence effect and high initial data requirements in parameter estimation of proportional hazards model,this paper proposes step estimation algorithm based on maximum likelihood estimation.Firstly,this paper has shape parameter and scale parameter estimated by maximum likelihood estimation,then it substitutes the result into proportional hazards model,and then it has regression coefficient of condition index estimated by maximum likelihood estimation.Through the calculation example analysis,step estimation algorithm is more effective in computation speed,computation time and data sample selection.

proportional hazards model;parameter estimation;maximum likelihood estimation;step estimation algorithm

TM41

B

1001-9898(2017)05-0018-03

本文责任编辑：靳书海