楊志遠(yuǎn)， 趙建民， 程中華， 李俐瑩， 池闊

(1.陸軍工程大學(xué)石家莊校區(qū) 裝備指揮與管理系， 河北 石家莊 050003；2.河北科技大學(xué) 信息學(xué)院， 河北 石家莊 050000)

0 引言

隨著科學(xué)技術(shù)的進(jìn)步和生產(chǎn)工藝的改進(jìn)，現(xiàn)代軍用裝備和民用產(chǎn)品的結(jié)構(gòu)日益復(fù)雜，功能集成程度越來(lái)越高。復(fù)雜系統(tǒng)往往存在多個(gè)性能指標(biāo)，性能退化會(huì)導(dǎo)致系統(tǒng)故障，此外，內(nèi)部或外部沖擊也可能導(dǎo)致系統(tǒng)發(fā)生故障。由于系統(tǒng)部件受共同的運(yùn)行環(huán)境、應(yīng)力、結(jié)構(gòu)設(shè)計(jì)等因素影響，不同退化過程(即性能指標(biāo)退化過程)及沖擊對(duì)不同退化過程的影響間往往存在相關(guān)關(guān)系[1]。如現(xiàn)在應(yīng)用較廣泛的微機(jī)電系統(tǒng)(MEMS)中，通常存在摩擦磨損、疲勞裂紋、電介質(zhì)老化等多種退化，以及由于過載沖擊導(dǎo)致的部件損傷或失效[2-3]，而由于共同的環(huán)境應(yīng)力狀態(tài)，這些故障模式間往往存在相關(guān)關(guān)系。因此，在評(píng)估此類系統(tǒng)的可靠性時(shí)，需要考慮多種失效模式及退化相關(guān)性的影響。

傳統(tǒng)可靠性研究是基于壽命數(shù)據(jù)展開的，由于現(xiàn)代工業(yè)產(chǎn)品可靠性的提高，在短時(shí)間內(nèi)難以得到足夠的壽命數(shù)據(jù)。但是，傳感器技術(shù)的發(fā)展使得產(chǎn)品狀態(tài)監(jiān)測(cè)很容易實(shí)現(xiàn)，因此基于狀態(tài)監(jiān)測(cè)數(shù)據(jù)的可靠性和剩余壽命預(yù)測(cè)研究應(yīng)用越來(lái)越廣泛[4-5]?；谕瑯拥览恚跔顟B(tài)監(jiān)測(cè)數(shù)據(jù)對(duì)故障相關(guān)進(jìn)行研究在實(shí)際中更容易實(shí)現(xiàn)。性能退化模型是進(jìn)行系統(tǒng)可靠性分析的基礎(chǔ)，目前主要包括退化軌跡模型、退化量分布模型和基于隨機(jī)過程的退化模型等[6]，相對(duì)于前兩種模型，基于隨機(jī)過程的退化模型能夠描述性能退化在時(shí)間軸上的不確定性，更符合工程實(shí)際。對(duì)存在多個(gè)相關(guān)性能退化過程的系統(tǒng)可靠性分析時(shí)，首先需要建立退化相關(guān)性模型。退化相關(guān)性模型用于定量描述系統(tǒng)多個(gè)退化過程間的相關(guān)關(guān)系，主要包括多維退化模型[7-9]、退化率相關(guān)模型[10-12]和基于Copula函數(shù)的退化相關(guān)模型[13-16]。其中多維退化模型由于其函數(shù)結(jié)構(gòu)限制，應(yīng)用范圍較窄，且在涉及兩個(gè)以上部件退化時(shí)模型較為復(fù)雜。而退化率相關(guān)模型是考慮一個(gè)部件退化對(duì)其余部件退化速率的影響，在實(shí)際中二者之間的定量關(guān)系難以有效確定，限制了該方法的應(yīng)用研究。相對(duì)而言，Copula函數(shù)在定量描述隨機(jī)變量相關(guān)性方面具有很高的靈活性，并且提供了分析相關(guān)結(jié)構(gòu)對(duì)可靠性影響的參數(shù)化方法，在金融、能源、氣象等領(lǐng)域已得到廣泛應(yīng)用[17]。此外，在工程實(shí)際中，由于影響系統(tǒng)不同退化過程間相關(guān)性的因素十分復(fù)雜，很多情況下難以對(duì)其機(jī)理進(jìn)行全面分析，而Copula函數(shù)將相關(guān)性結(jié)構(gòu)與邊緣分布分開進(jìn)行考慮，通過分析退化數(shù)據(jù)即可對(duì)系統(tǒng)不同退化過程間的相關(guān)性進(jìn)行定量描述，避免了上述問題的出現(xiàn)。

在可靠性模型中，系統(tǒng)同時(shí)受退化過程和外界沖擊過程的模型稱為退化閾值沖擊(DTS)模型。在DTS模型研究中，往往需要考慮退化過程和外界沖擊過程的相互關(guān)聯(lián)[18]。文獻(xiàn)[18-20]研究了沖擊過程對(duì)退化過程的影響，發(fā)現(xiàn)沖擊會(huì)導(dǎo)致系統(tǒng)退化量的增加；文獻(xiàn)[21]考慮了退化過程對(duì)沖擊過程到達(dá)率的影響，分析了系統(tǒng)可靠性；文獻(xiàn)[22]考慮了沖擊過程對(duì)退化過程和沖擊失效閾值的影響；文獻(xiàn)[23]在系統(tǒng)可靠性分析中考慮了退化過程和沖擊過程的相互影響，發(fā)現(xiàn)沖擊會(huì)導(dǎo)致系統(tǒng)退化量突增，系統(tǒng)退化過程會(huì)影響沖擊過程失效閾值；文獻(xiàn)[24]考慮了沖擊過程對(duì)沖擊失效閾值及退化速率的影響，建立了系統(tǒng)可靠性模型。在上述關(guān)于DTS模型的研究中，為計(jì)算方便，大部分是基于正態(tài)分布退化軌跡模型開展的，且均考慮了退化過程與沖擊過程的相互影響，并未分析不同退化過程間的相關(guān)性。如上文分析，復(fù)雜系統(tǒng)通常存在多個(gè)退化過程，且不同退化過程間往往是相關(guān)的。

本文考慮3種相關(guān)關(guān)系：1)沖擊過程會(huì)導(dǎo)致系統(tǒng)退化量的突增；2)不同退化過程間的相關(guān)關(guān)系；3)沖擊導(dǎo)致不同退化過程退化突增量之間的相關(guān)關(guān)系。由于每次外界沖擊都會(huì)導(dǎo)致系統(tǒng)退化量的突增，即系統(tǒng)不同退化過程的退化突增量受到同一沖擊過程的影響，因此考慮退化突增量之間的相關(guān)性是合理的。在此基礎(chǔ)上，本文基于Gamma過程和Copula函數(shù)建立系統(tǒng)可靠度模型，并提供相應(yīng)參數(shù)估計(jì)方法。

1 系統(tǒng)描述

系統(tǒng)包含m個(gè)退化過程，并受到外界沖擊的影響。外界沖擊會(huì)引起性能退化量的突增，系統(tǒng)實(shí)際退化由連續(xù)退化過程和沖擊引起的累積損傷組成。系統(tǒng)會(huì)發(fā)生退化故障和沖擊故障兩種故障模式。當(dāng)系統(tǒng)任意退化過程超過其特定的故障閾值后，即發(fā)生退化故障；當(dāng)外界沖擊幅值超過某一強(qiáng)度時(shí)，系統(tǒng)發(fā)生沖擊故障。

令Di(t)表示系統(tǒng)第i個(gè)退化過程在時(shí)刻t的總體退化量,i=1,2,…,m，m為退化過程數(shù)量。假定系統(tǒng)初始退化量為0，即Di(0)=0，相應(yīng)的故障閾值為L(zhǎng)i；Wj(j=1,2,3,4)為第j次沖擊的幅值，系統(tǒng)沖擊故障閾值為A. 基于此，系統(tǒng)故障過程如圖1所示。圖1中曲線表示系統(tǒng)退化過程，包括多個(gè)連續(xù)退化過程和外界沖擊導(dǎo)致的退化量突增部分，tj表示第j次沖擊到達(dá)時(shí)刻，Yij表示第j次沖擊對(duì)第i個(gè)性能指標(biāo)造成的退化增量。

圖1 系統(tǒng)故障過程Fig.1 System failure process

由圖1可以看出，系統(tǒng)的故障發(fā)生是由沖擊過程和m個(gè)退化過程競(jìng)爭(zhēng)失效造成的。如上所述，在本文中考慮系統(tǒng)退化過程相關(guān)性主要包括如下3個(gè)方面：

1)沖擊過程會(huì)導(dǎo)致系統(tǒng)退化量的突增；

2)不同退化過程間的相關(guān)關(guān)系；

3)沖擊導(dǎo)致不同退化過程退化突增量之間的相關(guān)關(guān)系。

1.1 退化過程模型

根據(jù)系統(tǒng)描述，退化過程Di(t)由兩部分組成：一是在時(shí)間t內(nèi)不考慮沖擊影響的退化量Xi(t)；二是沖擊引起的累積損傷量Si(t)。則Di(t)可表示為Di(t)=Xi(t)+Si(t)。

首選考慮退化過程Xi(t)。退化過程通常具有隨機(jī)性，而隨機(jī)過程是刻畫系統(tǒng)退化過程中不確定性的有效方法。其中，Gamma過程由于其良好的性質(zhì)，廣泛應(yīng)用于連續(xù)退化過程建模，包括材料的磨損、侵蝕、裂紋等[25]。本文采用Gamma過程來(lái)描述退化過程Xi(t)。因此，在任意時(shí)間間隔[u,t]t>u≥0內(nèi)，退化增量Xi(t)-Xi(u)均服從Gamma分布Γ(αi(t-u),βi)，其概率密度函數(shù)為

(1)

(2)

下面考慮沖擊引起的累積損傷模型，通常假設(shè)外界沖擊到達(dá)服從Poisson過程。令λ(t)表示t時(shí)刻Poisson過程到達(dá)率，則在t時(shí)間內(nèi)沖擊發(fā)生次數(shù)N(t)服從以下分布：

(3)

在沖擊未造成系統(tǒng)故障時(shí)，對(duì)于同一退化過程，退化增量Yij為獨(dú)立同分布的非負(fù)隨機(jī)變量,j=1,2,…,N(t)?；诖?，在時(shí)間t內(nèi)對(duì)于第i個(gè)退化過程，由沖擊引起的累積退化量為

(4)

假設(shè)Yij的分布函數(shù)為FYi(si)，si為沖擊對(duì)第i個(gè)退化過程造成的退化增量值，則Si(t)的分布函數(shù)可表示為

(5)

(6)

式中：Φ(·)為標(biāo)準(zhǔn)正態(tài)分布函數(shù)。

1.2 相關(guān)性退化分析

1.2.1 Copula函數(shù)簡(jiǎn)介

Copula函數(shù)能夠全面描述多元隨機(jī)變量相關(guān)結(jié)構(gòu)，是研究隨機(jī)變量相關(guān)性的有力工具。Copula理論通過分別考慮邊緣分布和相關(guān)結(jié)構(gòu)，在Sklar定理基礎(chǔ)上構(gòu)造隨機(jī)變量聯(lián)合分布。本文中假設(shè)系統(tǒng)存在m個(gè)退化過程，因此這里對(duì)m維Copula函數(shù)進(jìn)行描述。令多元隨機(jī)變量(X1，X2，…，Xm)的邊緣分布函數(shù)分別為F1(x1)，F(xiàn)2(x2)，…，F(xiàn)m(xm)，則基于m維Copula函數(shù)C(u1，u2，…，um)的聯(lián)合分布函數(shù)可表示為

H(x1，x2，…，xm)=
C(F1(x1)，F(xiàn)2(x2)，…，F(xiàn)m(xm)|θ)=
C(Fi(xi)i=1,2,…,m|θ)，

(7)

式中：C(F1(x1)，F(xiàn)2(x2)，…，F(xiàn)m(xm)|θ)表示m維Copula函數(shù)，θ為Copula函數(shù)中的參數(shù)向量，其取值直接影響著隨機(jī)變量間相關(guān)關(guān)系強(qiáng)弱。當(dāng)多元隨機(jī)變量相互獨(dú)立時(shí)，C(u1，u2，…，ui,…,um)=u1u2…ui…um，ui為取值范圍為[0,1]的隨機(jī)變量。

上述聯(lián)合分布的概率密度函數(shù)可表示為

(8)

由于在描述隨機(jī)相關(guān)方面的靈活性和有效性，Copula函數(shù)被廣泛應(yīng)用于金融、社會(huì)、環(huán)境等許多領(lǐng)域。因此，在本文中采用Copula函數(shù)對(duì)退化過程相關(guān)性進(jìn)行建模分析。

1.2.2 相關(guān)性退化模型

如1.2.1節(jié)所述，本文中的退化過程相關(guān)性包括：1)退化過程[X1(t),X2(t),…,Xm(t)]之間的相關(guān)性；2)沖擊對(duì)不同退化過程造成的退化增量[S1(t),S2(t),…,Sm(t)]之間的相關(guān)性。在不同時(shí)刻對(duì)系統(tǒng)退化水平進(jìn)行測(cè)量，令Di(tl)表示第i個(gè)退化過程在第l次測(cè)量的退化量值。為便于實(shí)際應(yīng)用，這里假設(shè)在不同時(shí)間間隔內(nèi)，不同退化過程及沖擊引起的退化過程增量間的相關(guān)性可忽略。即若l≠l′，則Di(tl)-Di(tl-1)與Di′(tl′)-Di′(tl′-1)相互獨(dú)立，否則需要考慮其相關(guān)關(guān)系。其中，i′=1，2，…，m，且i≠i′. 在文獻(xiàn)[13]中也采用了類似假設(shè)，分析多元Wiener退化過程的相關(guān)性。在此假設(shè)下，基于Sklar定理，[X1(t),X2(t),…,Xm(t)]在相同時(shí)間間隔[u,t]t>u≥0內(nèi)增量的聯(lián)合分布函數(shù)可表示為

HX(t-u)(x1,x2,…,xm)=CX(FXi(t-u)(xi)i=1,2,…,m|θ1)，

(9)

[S1(t),S2(t),…,Sm(t)]在相同時(shí)間間隔[u,t]t>u≥0內(nèi)增量的聯(lián)合分布函數(shù)可表示為

HS(t-u)(s1,s2,…,sm)=CS(FSi(t-u)(si)i=1,2,…,m|θ2)，

(10)

式中：CX(·)與CS(·)表示相應(yīng)的Copula函數(shù)；θ1和θ2為對(duì)應(yīng)的Copula參數(shù)向量。在此基礎(chǔ)上，令u=0，即可分別得到t時(shí)刻[X1(t),X2(t),…,Xm(t)]和[S1(t),S2(t),…,Sm(t)]的聯(lián)合分布函數(shù)。

2 系統(tǒng)可靠度模型

根據(jù)第1節(jié)對(duì)系統(tǒng)故障過程的描述，可以看出系統(tǒng)有退化故障和沖擊失效兩種故障模式，在建立可靠度模型時(shí)，需要同時(shí)考慮兩種故障模式的存在以及退化間的相關(guān)性。

假設(shè)每次的沖擊幅值Wj是獨(dú)立同分布的隨機(jī)變量，分布函數(shù)為Wj～FWj(w),j=1,2,…,N(t)。這里假設(shè)沖擊幅值服從正態(tài)分布，則Wj的分布函數(shù)為

(11)

式中：μW和σW分別為沖擊幅值Wj的期望和標(biāo)準(zhǔn)差。

在此基礎(chǔ)上，系統(tǒng)的可靠度函數(shù)為

(12)

式中：φ(·)表示標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)；cS(·)為Copula密度函數(shù)。由(12)式可以看出，當(dāng)系統(tǒng)存在的退化過程較多時(shí)，由于涉及到高維積分的運(yùn)算，系統(tǒng)可靠性的估計(jì)存在一定困難，這時(shí)可采用數(shù)值積分方法得到系統(tǒng)可靠度隨時(shí)間變化的情況。事實(shí)上，在工程實(shí)踐中，對(duì)于某個(gè)系統(tǒng)而言，導(dǎo)致其最終故障的通常是少數(shù)退化速率較快的退化過程，因此在可靠性估計(jì)時(shí)分析這些主要的退化過程即可。

由(12)式還可以看出，系統(tǒng)可靠度計(jì)算存在一定困難，特別是當(dāng)涉及退化過程較多時(shí)。系統(tǒng)可靠度計(jì)算的難點(diǎn)主要在于兩類退化存在相關(guān)性。由上文分析可知，不同退化過程Xi(t)和不同累計(jì)損傷過程Si(t)間均存在相關(guān)關(guān)系，從而導(dǎo)致總退化過程[D1(t),D2(t),…,Dm(t)]間也存在相關(guān)性。采用適當(dāng)?shù)腃opula函數(shù)直接描述總退化過程間的相關(guān)關(guān)系結(jié)構(gòu)，可以簡(jiǎn)化可靠度函數(shù)形式，得到系統(tǒng)可靠度的近似值?；赟klar定理，由(7)式可得[D1(t),D2(t),…,Dm(t)]在相同時(shí)間間隔[u,t]t>u≥0內(nèi)增量的聯(lián)合分布函數(shù)為

HD(t-u)(d1,d2,…,dm)=
CD(FDi(t-u)(di)i=1,2,…,m|θ0)，

(13)

式中：CD(FDi(t-u)(di)i=1,2,…,m|θ0)表示相應(yīng)Copula函數(shù)，F(xiàn)Di(t-u)(di)表示第i個(gè)退化過程在時(shí)間區(qū)間[u,t]內(nèi)增量的分布函數(shù)，di為第i個(gè)退化過程總退化量值，θ0為以上Copula函數(shù)中的參數(shù)向量。

為得到總退化過程在t時(shí)刻退化量Di(t)的分布函數(shù)，以Δt為單位將時(shí)間區(qū)間[0,t]離散化，即t=JΔt，J為足夠大的正整數(shù)。令ΔXi,b=Xi(bΔt)-Xi((b-1)Δt)，b=1,2,…,J，則有Xi(bΔt)=ΔXi,1+ΔXi,2+…+ΔXi,b. 根據(jù)Gamma過程性質(zhì)，可知ΔXi,b相互獨(dú)立，且均服從參數(shù)為αiΔt和βi的Gamma分布。由Lyapunov中心極限定理可得Di(JΔt)的近似分布為

(14)

在此基礎(chǔ)上將時(shí)間連續(xù)化，即可得到Di(t)的近似分布為

(15)

基于此，系統(tǒng)可靠度函數(shù)可近似為

(16)

一般而言，(16)式得到的系統(tǒng)可靠度與(12)式不完全相同，但其計(jì)算量較小，特別是在系統(tǒng)涉及的退化過程較多時(shí)，采用(12)式很難求得精確的結(jié)果，且計(jì)算量很大。只要(16)式得到的近似值精度較高，就可以在工程實(shí)際中應(yīng)用。

3 參數(shù)估計(jì)與Copula函數(shù)選擇

從第2節(jié)的分析可以看出，退化和可靠性模型中需要估計(jì)的參數(shù)為Θ=(θ1,θ2,λ,μW,σW,(αi,βi,μYi,σYi)i=1，2，…，m).

為便于分析，將參數(shù)向量記為Θ=(θ1,θ2,Θ1,Θ2)，其中Θ1=(λ,μW,σW)，Θ2=(αi,βi,μYi,σYi)i=1，2，…，m。由各參數(shù)含義可知，Θ1表示外界沖擊過程的特征，而Θ2表示系統(tǒng)邊際退化特征，剩余的(θ1,θ2)則表示退化相關(guān)關(guān)系特征。在相關(guān)研究中，通常假設(shè)外界沖擊過程特征Θ1是已知的，其具體值可以通過監(jiān)測(cè)外部環(huán)境沖擊(如應(yīng)力、振幅等突變)的到達(dá)時(shí)間和強(qiáng)度統(tǒng)計(jì)得到。因此，下面主要對(duì)Θ2和(θ1,θ2)的估計(jì)進(jìn)行分析。

3.1 退化過程模型參數(shù)估計(jì)

如果能夠得到系統(tǒng)各退化過程Xi(t)以及每次沖擊對(duì)不同退化過程造成的增量Yij數(shù)據(jù)，則可用邊際推斷法(IFM)分兩步得到參數(shù)估計(jì)值2和(1,2)，這種方法在文獻(xiàn)[26]中有詳細(xì)步驟和應(yīng)用說(shuō)明，不再贅述。但在工程實(shí)際中，由于在系統(tǒng)運(yùn)行過程中難以準(zhǔn)確識(shí)別沖擊到達(dá)時(shí)刻并區(qū)分連續(xù)退化過程和累積損傷過程引起的退化量值，Xi(t)和Yij的數(shù)據(jù)通常難以獲得。因此，上述參數(shù)估計(jì)方法的實(shí)用性不強(qiáng)。

假設(shè)取N個(gè)相同系統(tǒng)進(jìn)行退化實(shí)驗(yàn)，在系統(tǒng)故障前定期對(duì)其退化狀態(tài)進(jìn)行檢測(cè)。第q個(gè)系統(tǒng)進(jìn)行故障前檢測(cè)的次數(shù)記為Kq，q=1,2,…,N，第i個(gè)退化過程第k次檢測(cè)的退化量記為Di(tq,k)，樣本值為di(tq,k)，k=1,2,…,Kq. 基于這些數(shù)據(jù)對(duì)Θ2和(θ1,θ2)進(jìn)行估計(jì)。由于模型涉及的參數(shù)較多，為降低計(jì)算難度，借鑒IFM，首先得到Θ2的估計(jì)值，在已知Θ2的基礎(chǔ)上再對(duì)(θ1,θ2)進(jìn)行估計(jì)。為方便描述，定義以下符號(hào)：ΔDi,q,k=Di(tq,k)-Di(tq,k-1),Δdi,q,k=di(tq,k)-di(tq,k-1),Δtq,k=tq,k-tq,k-1，tq,0=0.

根據(jù)Gamma過程和齊次Poisson過程特性，可得ΔDi,q,k的分布函數(shù)為

(17)

對(duì)應(yīng)的概率密度函數(shù)為

(18)

假設(shè)在時(shí)間區(qū)間Δtq,k內(nèi)的沖擊次數(shù)N(Δtq,k)已知，如果在實(shí)際中不能準(zhǔn)確記錄沖擊到達(dá)次數(shù)，則可以用E[N(Δtq,k)]=λΔtq,k代替。

根據(jù)Gamma過程和Poisson過程的獨(dú)立增量特性，可得退化數(shù)據(jù)樣本的對(duì)數(shù)似然函數(shù)為

(19)

式中：di表示實(shí)驗(yàn)中系統(tǒng)第i個(gè)退化過程樣本數(shù)據(jù)矩陣。

在上述大似然函數(shù)中，fΔDi,q,k(x)的表達(dá)式可能難以得到，在實(shí)際應(yīng)用中可采用數(shù)值微分的方法求解。此外，根據(jù)(14)式，在實(shí)驗(yàn)中取Δtq,k為較大的值時(shí)，可以近似認(rèn)為ΔDi,q,k服從正態(tài)分布。因此上述樣本最大似然函數(shù)可轉(zhuǎn)化為

(20)

在i=1,2,…,m的情況下最大化對(duì)數(shù)似然函數(shù)，即可得到Θ2的估計(jì)值2.

3.2 Copula相關(guān)性參數(shù)估計(jì)

在3.1節(jié)基礎(chǔ)上，對(duì)(θ1,θ2)的值進(jìn)行估計(jì)。由相關(guān)性退化模型可得[ΔD1,q,k,ΔD2,q,k,…,ΔDm,q,k]的聯(lián)合分布函數(shù)為

(21)

其對(duì)應(yīng)的概率密度函數(shù)為

(22)

由此可得退化數(shù)據(jù)樣本的輪廓對(duì)數(shù)似然函數(shù)為

lnL(θ1,θ2|d,2)=

(23)

式中：d表示全部退化樣本數(shù)據(jù)矩陣。

通過最大化(23)式即可得到(θ1,θ2)的最大似然估計(jì)值(1,2)。然而，由于樣本似然函數(shù)涉及到對(duì)多維積分求微分，且積分函數(shù)本身較為復(fù)雜，很難得到明確的似然函數(shù)解析表達(dá)形式，直接最大化(23)式得到參數(shù)估計(jì)值的思路很難實(shí)現(xiàn)。為此，本文提出一種基于貝葉斯方法的間接參數(shù)估計(jì)方法。

在貝葉斯方法中，將模型中的參數(shù)作為隨機(jī)變量處理。假設(shè)參數(shù)(θ1,θ2)的先驗(yàn)分布密度函數(shù)為π(θ1,θ2)，根據(jù)貝葉斯公式可得

(24)

式中：π(θ1,θ2|d)為(θ1,θ2)的后驗(yàn)分布密度函數(shù)；L(d|θ1,θ2)為似然函數(shù)，即L(θ1,θ2|d,2)。在得到π(θ1,θ2|d)的后驗(yàn)樣本后，即可將樣本均值作為參數(shù)的估計(jì)值。由于L(d|θ1,θ2)的復(fù)雜性，要得到π(θ1,θ2|d)很困難。為了避免多維積分運(yùn)算，令S=(siqk)m×N×Kq為模型參數(shù)，定義新的似然函數(shù)：

(25)

同時(shí)，構(gòu)造新的先驗(yàn)分布為

(26)

根據(jù)貝葉斯公式可知

π(θ1,θ2,S|d)∝π(θ1,θ2,S)L(d|θ1,θ2,S).

(27)

基于此，通過抽樣π(θ1,θ2,S|d)得到(θ1,θ2,S)的后驗(yàn)樣本，然后可得相關(guān)參數(shù)的估計(jì)值。關(guān)于后驗(yàn)樣本的抽樣，可采用馬爾可夫鏈蒙特卡洛(MCMC)方法來(lái)解決，MCMC方法在高維總體或復(fù)雜總體取樣問題中應(yīng)用已趨于成熟，也已經(jīng)有貝葉斯統(tǒng)計(jì)專業(yè)軟件能夠簡(jiǎn)便地實(shí)現(xiàn)MCMC方法，如WinBUGS軟件。MCMC方法具體原理和用法可參考文獻(xiàn)[27]，不再贅述。

下面對(duì)近似系統(tǒng)可靠度函數(shù)(16)式中的θ0進(jìn)行估計(jì)。在(14)式基礎(chǔ)上，容易得到相應(yīng)的對(duì)數(shù)似然函數(shù)為

lnL(θ0|d,2)=

(28)

可以看出，在似然函數(shù)中有一部分與θ0的計(jì)算無(wú)關(guān)，因此可以忽略。(28)式可簡(jiǎn)化為

lnL(θ0|d,2)=

(29)

通過最大化(29)式，即可得到θ0的估計(jì)值0. 當(dāng)然，同時(shí)以Θ2和θ0為未知參數(shù)直接最大化(28)式，即可一次性得到對(duì)應(yīng)估計(jì)值。

顯然，通過上述分兩步估計(jì)方法得到的估計(jì)值并不能保證樣本的完全似然函數(shù)達(dá)到最大，從而使得估計(jì)值可能不是全局最優(yōu)結(jié)果，但這種方法卻大大減少了計(jì)算量，特別是對(duì)于系統(tǒng)退化過程較多的情形。而且，這種方法得到的結(jié)果與直接采用極大似然估計(jì)得到的結(jié)果很接近，可以滿足工程應(yīng)用的需求[26]。

3.3 Copula函數(shù)選擇方法

關(guān)于Copula函數(shù)形式的選擇，在參數(shù)估計(jì)結(jié)果基礎(chǔ)上，可采用赤池信息準(zhǔn)則(AIC)進(jìn)行判斷。AIC表達(dá)式為

ζ=-2{max lnL(Θ)}+2p，

(30)

式中：ζ為AIC值；p為數(shù)學(xué)模型中未知參數(shù)的個(gè)數(shù)；L(Θ)為似然函數(shù)。根據(jù)AIC判斷準(zhǔn)則，具有最小AIC值的Copula相關(guān)結(jié)構(gòu)對(duì)數(shù)據(jù)的擬合效果最好。

4 算例分析

本文以美國(guó)Sandia國(guó)家實(shí)驗(yàn)室研發(fā)的一款微型發(fā)動(dòng)機(jī)為例，對(duì)競(jìng)爭(zhēng)失效系統(tǒng)可靠度模型進(jìn)行驗(yàn)證。該微型發(fā)動(dòng)機(jī)在運(yùn)行過程中同時(shí)存在磨損退化與沖擊過程，由部件間摩擦造成的磨損是微型發(fā)動(dòng)機(jī)的主要失效模式，同時(shí)沖擊會(huì)造成一定量的磨損碎片，當(dāng)沖擊過大時(shí)會(huì)導(dǎo)致其失效，關(guān)于該微型發(fā)動(dòng)機(jī)的具體結(jié)構(gòu)可參考文獻(xiàn)[28]。由微型發(fā)動(dòng)機(jī)失效機(jī)理可以看出，其失效過程同時(shí)包含了退化過程和沖擊過程。為了對(duì)本文的模型加以驗(yàn)證，在參考文獻(xiàn)[28]原有數(shù)據(jù)基礎(chǔ)上，假設(shè)系統(tǒng)存在兩個(gè)磨損退化過程。在原始失效過程實(shí)驗(yàn)中，采用退化軌跡模型對(duì)磨損過程進(jìn)行描述，系統(tǒng)在某一時(shí)刻的磨損量服從正態(tài)分布。由于磨損過程是非減的，采用Gamma過程對(duì)其進(jìn)行描述更為符合工程實(shí)際。考慮該系統(tǒng)存在兩個(gè)相同退化過程和一個(gè)沖擊過程，相關(guān)參數(shù)如表1所示。

表1 微型發(fā)動(dòng)機(jī)可靠性模型參數(shù)值

在Copula函數(shù)的選擇上，由于阿基米德Copula族具有較為簡(jiǎn)單明晰的解析表達(dá)式，且結(jié)構(gòu)簡(jiǎn)單、具有良好的計(jì)算性質(zhì)，被廣泛應(yīng)用于金融、社會(huì)、工業(yè)等許多領(lǐng)域。這里以阿基米德Copula為例，對(duì)上述可靠度模型進(jìn)行說(shuō)明驗(yàn)證。阿基米德Copula族中的Gumbel Copula函數(shù)對(duì)變量上尾反映較為敏感，用于表示退化過程間的相關(guān)性；Clayton Copula函數(shù)對(duì)變量下尾變化反映靈敏，用于表示沖擊導(dǎo)致退化增量間的相關(guān)性。當(dāng)然，在工程實(shí)際中，可根據(jù)系統(tǒng)特點(diǎn)和退化特征對(duì)具體Copula相關(guān)結(jié)構(gòu)函數(shù)進(jìn)行選擇。兩種Copula函數(shù)的表達(dá)式如下：

Gumbel Copula:

CGum(u1,u2,…,un)=

Clayton Copula:

其中，τ∈[0,1]為Kendallτ指標(biāo)，能夠反映變量間的相關(guān)程度。由Kendallτ指標(biāo)與Copula參數(shù)θ之間的關(guān)系可以看出，對(duì)于Gumbel Copula函數(shù)，當(dāng)θ=1時(shí)CGum描述的隨機(jī)變量間的關(guān)系相互獨(dú)立，當(dāng)θ→∞時(shí)CGum描述的隨機(jī)變量間的關(guān)系趨于完全正相關(guān)。對(duì)于Clayton Copula函數(shù)，當(dāng)θ→0時(shí)CCla描述的隨機(jī)變量間的關(guān)系相互獨(dú)立，當(dāng)θ→∞時(shí)CGum描述的隨機(jī)變量間的關(guān)系趨于完全正相關(guān)。

4.1 系統(tǒng)可靠性分析

在以上參數(shù)和Copula函數(shù)基礎(chǔ)上，對(duì)微型發(fā)動(dòng)機(jī)系統(tǒng)可靠性和故障密度函數(shù)進(jìn)行計(jì)算，結(jié)果如圖2所示。

圖2 不同相關(guān)度下系統(tǒng)可靠度及故障概率密度Fig.2 System reliability and failure probability density under different dependency degrees

圖2中分別給出了τGum=0、τCla=0，τGum=0.5、τCla=0.8以及τGum=0.9、τCla=0.9時(shí)的系統(tǒng)可靠度和故障密度函數(shù)隨時(shí)間的變化情況。由圖2可以看出，退化相關(guān)性對(duì)系統(tǒng)可靠度和故障密度函數(shù)有一定影響，隨著系統(tǒng)不同退化過程以及沖擊導(dǎo)致退化增量間相關(guān)性的提高，系統(tǒng)可靠度也呈增高的趨勢(shì)。在這類存在退化相關(guān)性的系統(tǒng)中，獨(dú)立性假設(shè)導(dǎo)致系統(tǒng)可靠度估計(jì)過于保守，可能會(huì)降低系統(tǒng)使用和維修的經(jīng)濟(jì)性。此外，在系統(tǒng)運(yùn)行40 000 r之前，不同相關(guān)性程度下的系統(tǒng)可靠度基本相同，這是因?yàn)橄到y(tǒng)前期故障主要是外界沖擊引起的。

4.2 可靠度估算模型精確性分析

在下面的分析中令相關(guān)性程度為τGum=0.5、τCla=0.8，為驗(yàn)證可靠度估算模型(16)式的精確性，首先需要對(duì)系統(tǒng)退化過程進(jìn)行采樣。以時(shí)間間隔10 000 r對(duì)系統(tǒng)退化過程以及沖擊過程進(jìn)行模擬，得到1 000個(gè)采樣點(diǎn)數(shù)據(jù)如圖3所示。

圖3 相關(guān)性退化過程采樣結(jié)果Fig.3 Sampling results of dependent degradation processes

由圖3可以看出，樣本點(diǎn)數(shù)據(jù)綜合了Gumbel Copula和Clayton Copula的特點(diǎn)，在上尾和下尾部分均展示出較高的相關(guān)性。在此基礎(chǔ)上，對(duì)可靠度估算模型(16)式中的Copula函數(shù)進(jìn)行選擇并估計(jì)相關(guān)性參數(shù)θ0，結(jié)果如表2所示。由表2可以看出，t-Copula的AIC值最小，因此在可靠度估算模型(16)式中選取t-Copula描述退化間的相關(guān)關(guān)系，相應(yīng)的參數(shù)值在表2中已給出。在此基礎(chǔ)上，對(duì)系統(tǒng)可靠度的實(shí)際值和估計(jì)值進(jìn)行比較，結(jié)果如圖4所示。

表2 不同Copula函數(shù)參數(shù)估計(jì)值和AIC值

圖4 系統(tǒng)可靠度實(shí)際值與估計(jì)值對(duì)比Fig.4 Comparison of the real and estimated values of system reliability

圖4中分別給出了系統(tǒng)可靠度實(shí)際值和估計(jì)值的變化以及二者間差值的變化情況。由圖4可以看出，系統(tǒng)可靠度實(shí)際值和估計(jì)值基本相同，可靠度之間的差值最大不超過0.002. 基于此，可靠度估算模型(16)式的精確度較高，可以在一定范圍內(nèi)代替可靠度模型(12)式，從而簡(jiǎn)化計(jì)算。

4.3 退化相關(guān)性影響分析

沖擊到達(dá)的頻率對(duì)系統(tǒng)可靠性有直接影響，下面分析沖擊過程到達(dá)率參數(shù)λ對(duì)系統(tǒng)可靠度的影響。分別令λ=5×10-4、λ=5×10-5和λ=5×10-6，在此基礎(chǔ)上，得到在退化相關(guān)和不相關(guān)情況下的系統(tǒng)可靠度，如圖5所示。

圖5 參數(shù)λ對(duì)系統(tǒng)可靠度影響Fig.5 Influence of λ on system reliability

圖5(a)給出了不同情況下系統(tǒng)的可靠度變化情況。由圖5(a)可知，隨著λ的減小，即外界沖擊次數(shù)的減少，系統(tǒng)可靠度明顯提升；對(duì)于不同的λ值，退化相關(guān)與不相關(guān)情況下系統(tǒng)可靠度值的差異也有所不同。由圖5(b)具體分析其差異，可以看出：隨著λ的減小，在退化相關(guān)與不相關(guān)情況下系統(tǒng)可靠度的差值不斷增大。這是因?yàn)殡S著λ的減小，系統(tǒng)的主要故障模式由沖擊故障轉(zhuǎn)化為退化故障，所以此時(shí)退化相關(guān)性對(duì)系統(tǒng)可靠度的影響也會(huì)變大。此外，通過對(duì)比相關(guān)與不相關(guān)情況下系統(tǒng)可靠度差值變化，可以看出：在系統(tǒng)運(yùn)行前期，考慮和不考慮退化相關(guān)性對(duì)系統(tǒng)可靠度計(jì)算結(jié)果的影響不大。這是因?yàn)樵缙谙到y(tǒng)可靠度下降主要是由于沖擊故障概率提高所致，此時(shí)發(fā)生退化故障的概率很小，所以相應(yīng)的退化相關(guān)性對(duì)系統(tǒng)可靠度的影響可以忽略。特別地，如果這里不考慮沖擊故障，即沖擊只引起相應(yīng)退化增量，而不會(huì)導(dǎo)致系統(tǒng)故障，則當(dāng)λ=5×10-5時(shí)，在退化相關(guān)和不相關(guān)情況下的系統(tǒng)可靠度如圖6所示。

圖6 不考慮沖擊故障下系統(tǒng)可靠度Fig.6 System reliability without considering shock failure

由圖6可以看出，如果只考慮系統(tǒng)退化故障，則在退化相關(guān)和不相關(guān)情況下，系統(tǒng)可靠度有明顯區(qū)別?；谝陨戏治觯诟?jìng)爭(zhēng)失效系統(tǒng)可靠度評(píng)估中，如果系統(tǒng)主要失效模式為沖擊失效，則為方便計(jì)算，可以忽略系統(tǒng)不同退化過程間的相關(guān)性；否則必須充分分析退化相關(guān)性對(duì)系統(tǒng)可靠度的影響。

5 結(jié)論

本文針對(duì)同時(shí)存在多個(gè)退化過程和外界沖擊過程的競(jìng)爭(zhēng)失效系統(tǒng)，在考慮沖擊過程會(huì)導(dǎo)致系統(tǒng)退化量突增基礎(chǔ)上，考慮不同退化過程間的相關(guān)關(guān)系和沖擊導(dǎo)致不同退化過程退化突增量之間的相關(guān)關(guān)系?；贕amma過程和Copula函數(shù)得到了系統(tǒng)可靠性模型和解析表達(dá)式，為簡(jiǎn)化模型計(jì)算，給出了可靠性模型的近似結(jié)果?；谪惾~斯理論給出了模型參數(shù)估計(jì)方法，并給出了Copula函數(shù)選擇方法。通過實(shí)例分析，驗(yàn)證了可靠性近似模型的精確性，分析了相關(guān)性參數(shù)對(duì)系統(tǒng)可靠度的影響。結(jié)果表明：對(duì)于存在退化相關(guān)性且主要失效模式為退化失效的系統(tǒng)，獨(dú)立性假設(shè)導(dǎo)致系統(tǒng)可靠度估計(jì)過于保守，這可能會(huì)降低系統(tǒng)使用和維修的經(jīng)濟(jì)性。