基于Copula函數(shù)的導(dǎo)彈部件非線性退化研究*

賀志遠(yuǎn)，呂衛(wèi)民，胡文林

(海軍航空大學(xué)， 山東煙臺 264001)

0 引言

導(dǎo)彈部件的許多故障是由于內(nèi)部材料的性能改變，比如：金屬材料的腐蝕，橡膠材料的老化以及電子元器件的失效所導(dǎo)致。研究導(dǎo)彈部件的性能退化過程，合理制定維修決策，對于裝備的保障工作具有重要價值。

隨著工業(yè)制造水平的提高，現(xiàn)代產(chǎn)品的失效機理變得越來越復(fù)雜。產(chǎn)品具有復(fù)雜內(nèi)部結(jié)構(gòu)和許多功能，退化過程往往有兩個或多個特征參量，并且它們可能是相關(guān)的。過去關(guān)于退化分析的大多數(shù)研究只涉及一個性能特征，當(dāng)這個關(guān)鍵性能特征參數(shù)值降低到失效閾值時，產(chǎn)品失效[1-3]。因此，正確處理各性能退化量之間的相關(guān)性，建立多性能退化模型是十分必要的。

目前，在性能多退化模型中，最簡單的情況是各特征量線性相關(guān)。文獻(xiàn)[4]采用主成分分析法，提出了一個PCA-CMAC模型來研究設(shè)備的性能退化。對于性能特征之間的復(fù)雜相關(guān)關(guān)系，常用函數(shù)來描述性能退化之間的多重相關(guān)性。文獻(xiàn)[5]采用逆時間尺度變換的逆高斯過程來構(gòu)建邊緣退化過程，通過Copula進(jìn)行多性能退化建模。文獻(xiàn)[6]認(rèn)為性能退化特性受Winner過程和適應(yīng)的Copula函數(shù)控制，以將特征的退化路徑結(jié)合在一起。文獻(xiàn)[7]建立了具有時間變換的Wiener過程模型評估休眠系統(tǒng)的可靠性，并通過多變量Copula函數(shù)描述退化特征的相關(guān)性。

某導(dǎo)彈部件具有典型的多性能退化特征，對其退化傳統(tǒng)研究方法較為單一。主要存在兩個方面的難題：1)該導(dǎo)彈部件屬于典型的機電一體化設(shè)備，內(nèi)部包含了電子部分和機械部分，各部分退化規(guī)律不同，導(dǎo)致整體退化軌跡呈現(xiàn)非線性。傳統(tǒng)的線性退化模型已經(jīng)無法有效的研究其退化過程。2)退化過程中，往往存在多個具有退化趨勢的特征參量，只研究某一退化特征，很可能降低退化分析的準(zhǔn)確性。針對以上問題，文中首先建立基于Wiener過程的非線性退化模型，再利用Copula函數(shù)進(jìn)行多性能參數(shù)的相關(guān)性建模，最后通過實例分析，對導(dǎo)彈部件的退化過程進(jìn)行研究。

1 基于Wiener過程的性能退化建模

1.1 線性退化

設(shè)t時刻產(chǎn)品性能參數(shù)的退化量為X(t)，且退化過程可用線性Wiener過程進(jìn)行表達(dá)，則退化模型為：

X(t)=μt+σB(t)

(1)

式中：μ為漂移系數(shù)，σ為擴散系數(shù)，B(t)為標(biāo)準(zhǔn)布朗運動。

假設(shè)該產(chǎn)品的失效閾值為ω，其性能退化軌跡由式(1)的Wiener過程描述，則產(chǎn)品的壽命T可定義為：

T={t:X(t)≥ω|X(0)<ω}

(2)

設(shè)μ和σ為固定未知參數(shù)，由文獻(xiàn)[8]可知，產(chǎn)品壽命T服從逆高斯分布，其概率密度函數(shù)為：

(3)

相對應(yīng)的壽命分布函數(shù)為：

(4)

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

1.2 非線性退化

導(dǎo)彈壽命周期任務(wù)剖面多變，內(nèi)部設(shè)備的退化機理復(fù)雜，退化軌跡也并非全部與時間線性相關(guān)，根據(jù)這種退化特征，本節(jié)建立基于Wiener隨機過程的非線性退化模型來描述導(dǎo)彈設(shè)備的退化。假設(shè)該設(shè)備退化過程只有一個關(guān)鍵性能參數(shù)，則非線性退化模型表示如下：

(5)

式中：X(0)為初始時刻的性能退化量(為方便研究，令X(0)=0)；μ(t;θ)為非線性函數(shù)；σB為擴散系數(shù)；B(t)為標(biāo)準(zhǔn)布朗運動，顯然，當(dāng)μ(t;θ)=μ時，式(5)轉(zhuǎn)化為1.1所述的線性Wiener過程。

設(shè)σB為固定未知參數(shù)，μ(t;θ)是關(guān)于t的可導(dǎo)函數(shù)，則產(chǎn)品壽命T的概率密度函數(shù)為：

(6)

根據(jù)文獻(xiàn)[9]，設(shè)

μ(t;θ)=λbtb-1

(7)

式中：λ為給定的漂移系數(shù)；b為固定參數(shù)。

將式(7)代入式(6)可得：

(8)

假設(shè)模型參數(shù)空間不具有隨機效應(yīng)，但不同項在退化中具有可變性，用漂移系數(shù)λ表示這種可變性，為了研究隨機效應(yīng)λ，引入以下定理[6]：

引理：若Z～N(μ,σ2)，且A，B，C，ω∈R，則

(9)

(10)

相應(yīng)產(chǎn)品壽命T的概率分布函數(shù)為：

(11)

可靠度函數(shù)可表示為：

R(t)=1-FT(t)

(12)

2 基于Copula函數(shù)多性能退化建模

2.1 Copula函數(shù)的概念和性質(zhì)

Copula函數(shù)的核心思想[10]可通過Sklar定理實現(xiàn)，若HX，Y(x,y)為聯(lián)合分布函數(shù)，則存在一個Copula函數(shù)C(u,v)，對于x,y∈(-,)，滿足

HX，Y(x,y)=C(FX(x),FY(y))

(13)

式中：FX(x),FY(y)分別為X和Y的邊緣分布。

Copula具有多種構(gòu)造形式，包括代數(shù)方法、幾何方法以及Archimedean族[11]生成元構(gòu)造法。文中主要對比研究Gumbel Copula、Frank Copula、Clayton Copula和Gaussian Copula。

Gumbel:

C(u,v)=exp{-[(-lnu)θ+(-lnv)θ]1/θ}

(14)

式中：θ∈[1,)。當(dāng)θ→時，u,v→完全相關(guān)。

Frank:

(15)

式中：θ∈(-,0)∪(0,+)。當(dāng)θ>0時，u,v→正相關(guān)；當(dāng)θ<0時，u,v→負(fù)相關(guān)。

Clayton:

C(u,v)=max[(u-θ+v-θ-1)-1/θ,0]

(16)

式中：θ∈(0,)。當(dāng)θ→時，u,v→完全相關(guān)。

Gaussian:

(17)

式中：θ∈(-1,1)。當(dāng)θ=±1時，u,v完全相關(guān)。

2.2 相關(guān)性建模

產(chǎn)品有m個性能退化參數(shù)，分別為(X1,X2,…,Xm)，分別反映了產(chǎn)品退化過程中的不同特征。本節(jié)利用Copula函數(shù)對退化過程的不同特征進(jìn)行相關(guān)性建模，為了方便研究，選取具有代表性的二元性能退化參數(shù)進(jìn)行建模分析，并假設(shè)它們的退化軌跡符合第1.2節(jié)所描述的Wiener隨機過程。

設(shè)產(chǎn)品退化量X1和X2不互相獨立，對應(yīng)的CDF分別為FT(x1)和FT(x2)，聯(lián)合分布函數(shù)為F(x1,x2)。根據(jù)Sklar定理，可得:

F(x1,x2)=C(FT(x1),FT(x2))

(18)

根據(jù)Copula函數(shù)的性質(zhì)，聯(lián)合分布函數(shù)的概率密度函數(shù)為：

f(x1,x2)=c(FT(x1),FT(x2);α)·fT(x1)·fT(x2)

(19)

式中：c(fT(x1),fT(x2);α)為C(FT(x1),FT(x2))的概率密度函數(shù)；α為反映Copula函數(shù)中相關(guān)性的系數(shù)。

為求得相關(guān)系數(shù)α，由式(19)可得，模型的對數(shù)似然函數(shù)為：

(20)

式中：γ包含模型中所有未知參數(shù)。

3 參數(shù)估計

3.1 TSML法

對于1.1節(jié)線性退化模型，將式(4)代入式(20)，可得：

(21)

式中：γ=(α,μ1,σ1,μ2,σ2)

對于式(21)，可以直接進(jìn)行極大化處理，求得未知參數(shù)，但由于γ維數(shù)較高，求解過程復(fù)雜。因此，采用兩步極大似然估計法(TSML法)對線性退化模型進(jìn)行參數(shù)估計。主要思路是，先極大化式(20)中不含未知參數(shù)α的部分：

(22)

式中：γ1=(μ1,σ1)；γ2=(μ2,σ2)。

根據(jù)式(4)，分別對fT(x1)和fT(x2)進(jìn)行極大似然估計，通過二維搜索法求得相應(yīng)參數(shù)(μ1,σ1)和(μ2,σ2)。再進(jìn)行第二步似然估計，將相應(yīng)參數(shù)值代回式(20)中求得α的估計值。

3.2 MCMC算法

非線性退化模型復(fù)雜，未知參數(shù)較多，進(jìn)行相關(guān)性建模后將更加復(fù)雜，參數(shù)估計時，傳統(tǒng)的似然估計法已不再適用。本節(jié)采用一種基于貝葉斯理論的Gibbs算法[12]進(jìn)行參數(shù)估計。

設(shè)樣本γ為一個n維變量，γ=(γ1,γ2,…,γn)，其先驗分布為p(γ1,γ2,…,γn)。

Gibbs采樣的步驟[13]如下：

當(dāng)m足夠大時，γ(m)可以看作γ的真值，從而求得樣本γ的后驗分布q(γ1,γ2,…,γn)，進(jìn)而估計未知參數(shù)值。

4 實例分析

以2個月為周期對某導(dǎo)彈進(jìn)行彈上部件測試，獲取多個部件的測試參數(shù)，選取光纖陀螺儀作為研究對象。對于光纖陀螺儀主要采集了兩個退化特征參量，分別為零偏退化數(shù)據(jù)和標(biāo)度因數(shù)退化數(shù)據(jù)，退化數(shù)據(jù)如圖1所示。將零偏退化數(shù)據(jù)與標(biāo)度因數(shù)退化數(shù)據(jù)分別記為X1(t)和X2(t)，根據(jù)經(jīng)驗及工程實際，將失效閾值設(shè)定為ω1=0.5、ω2=1.5。

圖1 光纖陀螺儀退化數(shù)據(jù)

首先對數(shù)據(jù)進(jìn)行退化建模，分別建立基于Wiener過程的線性退化模型和非線性退化模型；其次，基于Copula函數(shù)的線性模型和非線性模型進(jìn)行相關(guān)性建模；最后，對不同模型分別進(jìn)行參數(shù)估計。

由于Copula函數(shù)表達(dá)形式多樣，不同的Copula函數(shù)會產(chǎn)生不同的退化概率分布，如果不能選擇合適的Copula函數(shù)，可能導(dǎo)致結(jié)果不正確。因此，對4種不同形式的Copula函數(shù)(Gumbel Copula，F(xiàn)rank Copula，Clayton Copula和Gaussian Copula)分別進(jìn)行建模以及參數(shù)估計，通過AIC準(zhǔn)則進(jìn)行檢驗，為相關(guān)性模型選擇合適的Copula函數(shù)，AIC值越小說明函數(shù)擬合效果越好。線性相關(guān)模型和非線性相關(guān)模型的AIC檢驗結(jié)果如表1所示。

通過AIC檢驗結(jié)果，線性相關(guān)模型采用Gumbel Copula函數(shù)擬合效果最好，而非線性相關(guān)模型采用Gaussian Copula函數(shù)擬合效果最好。采用相應(yīng)的函數(shù)，對產(chǎn)品的可靠度進(jìn)行分析，如圖2所示。曲線R1和R4為文中建立的非線性和線性相關(guān)模型評估的可靠度；曲線R2為傳統(tǒng)的單一性能模型評估的可靠度，選取零偏退化數(shù)據(jù)作為退化參量；曲線R3為根據(jù)光纖陀螺儀真實故障數(shù)據(jù)評估的可靠度。

表1 兩種相關(guān)模型的AIC檢驗結(jié)果

圖2 不同模型的可靠度曲線

由圖2可以看出曲線R1與真實可靠度曲線R3最為接近，證明文中建立的基于Copula函數(shù)的非線性退化模型能夠更好的描述光纖陀螺儀的退化軌跡。

5 結(jié)論

基于Copula函數(shù)建立了非線性退化模型，評估了某導(dǎo)彈部件的可靠度，通過與真實故障數(shù)據(jù)的評估結(jié)果對比，證明了模型的合理性。同時，得出以下結(jié)論：

1)許多導(dǎo)彈部件的某一性能退化軌跡與時間并非線性相關(guān)，傳統(tǒng)基于線性隨機過程的模型已經(jīng)無法滿足產(chǎn)品的可靠性評估的要求。文中研究的非線性過程能夠有效提高評估結(jié)果的準(zhǔn)確性。

2)多數(shù)具有機電一體化特征的導(dǎo)彈部件內(nèi)部結(jié)構(gòu)復(fù)雜，退化過程往往具有多個特征參量，利用Copula函數(shù)建立各參量之間的相關(guān)性模型，對于導(dǎo)彈部件的可靠性評估具有很好的適用性。

3)文中僅研究了具有退化軌跡的失效模式，實際上還存在其他失效模式，比如，競爭失效。各種失效模式之間的關(guān)系值得進(jìn)一步研究。