基于Copula函數(shù)的導(dǎo)彈多性能退化可靠性評(píng)估

戴宗亮， 李小兵， 王慧杰

(空軍工程大學(xué) 防空反導(dǎo)學(xué)院，陜西 西安 710051)

0 引 言

導(dǎo)彈在整個(gè)壽命周期中，絕大部分時(shí)間處于貯存狀態(tài)[1]。在長(zhǎng)期貯存過(guò)程中，導(dǎo)彈必然會(huì)受到周圍各種環(huán)境應(yīng)力的影響，難免引起設(shè)備性能參數(shù)的變化，導(dǎo)致整個(gè)系統(tǒng)功能異?；蚩煽慷冉档蚚2]，目前，常用的解決多性能退化參數(shù)可靠性評(píng)估的方法有2種: 1)假設(shè)多個(gè)性能參數(shù)之間相互獨(dú)立，可以按照串聯(lián)系統(tǒng)的方法進(jìn)行處理，具有簡(jiǎn)單快速的優(yōu)點(diǎn)，但分析不夠全面，評(píng)估結(jié)果誤差較大; 2)充分考慮性能參數(shù)之間的相關(guān)性，如聯(lián)合概率密度法和狀態(tài)空間法[3,4]，優(yōu)點(diǎn)是充分考慮了參數(shù)之間的相關(guān)性，但在特征參數(shù)較多時(shí)存在計(jì)算量大和建模困難。Copula函數(shù)[5,6]作為一種靈活、穩(wěn)健的相關(guān)性分析函數(shù)，在分析變量間相關(guān)性時(shí)具有諸多的優(yōu)點(diǎn)。

本文提出了基于Copula函數(shù)的多性能退化可靠性評(píng)估模型。

1 多性能參數(shù)退化數(shù)據(jù)可靠性模型

根據(jù)導(dǎo)彈的退化機(jī)理[7]，導(dǎo)彈性能的退化是m個(gè)退化參數(shù)直接影響的結(jié)果，其可靠度為

R(t)=P{Y1(t)≤D1,Y2(t)≤D2,…,Ym(t)≤Dm}

ym(t)dy1dy2…dym

(1)

式中f(y1(t),y2(t),…,ym(t))為y1(t),y2(t),…,ym(t)的聯(lián)合概率分布。由式(1)可知，只要知道聯(lián)合概率密度函數(shù)f(y1(t),y2(t),…,ym(t))，根據(jù)已知的失效閾值，即可求得多退化模型的可靠度R(t)。要知道這些性能參數(shù)在統(tǒng)計(jì)學(xué)上的相關(guān)性，通常可以通過(guò)性能參數(shù)間的協(xié)方差運(yùn)算獲得其相關(guān)性。其方差—協(xié)方差矩陣可表示為

(2)

當(dāng)cov(yi(t),yj(t))=0，表示第i和第j個(gè)性能參數(shù)獨(dú)立,i,j=1,2,…,m。當(dāng)V的非對(duì)角線元素均為零，則可認(rèn)為產(chǎn)品性能參數(shù)均相對(duì)獨(dú)立，可將其等效成串聯(lián)系統(tǒng)，可靠度為

R(t)=P{Y1(t)≤D1,Y2(t)≤D2,…，Ym(t)≤Dm}

=P{Y1(t)≤D1}·P1{Y2(t)≤D2}·…·

(3)

當(dāng)V的非對(duì)角線元素非全零，則可認(rèn)為產(chǎn)品性能參數(shù)不全相對(duì)獨(dú)立，則要求出其聯(lián)合概率密度函數(shù)f(y1(t),y2(t),…,ym(t))，以此求得可靠度。因此，建立多性能參數(shù)產(chǎn)品的可靠性模型的關(guān)鍵是建立求取各性能參數(shù)的聯(lián)合概率密度函數(shù)f(y1(t),y2(t),…,ym(t))。

如果性能參數(shù)退化量Y1(t),Y2(t),…,Ym(t)的邊緣分布均服從為正態(tài)分布，則多性能參數(shù)產(chǎn)品性能參數(shù)退化量的聯(lián)合分布密度函數(shù)為

(4)

式中Y=(Y1(t),Y2(t),…,Ym(t));μ和Σ分別為Y的均值和協(xié)方差矩陣。

綜上可得，該模型主要有兩個(gè)難點(diǎn):聯(lián)合概率密度f(wàn)(y1(t),y2(t),…,ym(t))的確定，對(duì)于參數(shù)邊緣分布不同的情況，理論推導(dǎo)其聯(lián)合分布十分繁瑣；在多退化性能相關(guān)的條件下，可靠度的計(jì)算，協(xié)方差矩陣V的大規(guī)模積分計(jì)算量較大。

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

2.1 Copula函數(shù)[8]

定義1m維Copula函數(shù)是一個(gè)定義域?yàn)閇0,1]m，值域?yàn)閇0,1]的函數(shù)，即函數(shù)C：[0,1]m→[0,1]，且滿足條件

1)C對(duì)于每個(gè)變量遞增；

2)對(duì)[0,1]m中任意uk=0,k≤m，滿足C(u1,…,uk,…,um)=0；

3)對(duì)[0,1]m中任意vk∈[0,1],k≤m，滿足C(1,1,…vk,…,1,1)=vk；

4)對(duì)[0,1]m中所有的ak,bk∈[0,1]m,k=1,2,…,m，且ak≤bk,k=1,2,…,m，有

(5)

對(duì)于二元Copula函數(shù)，任意u1,u2,v1,v2∈[0,1]，令u1≤u2,v1≤v2，則

C(u2,v2)-C(u1,v2)-C(u2,v1)+C(u1,v1)≥0

(6)

定義2 若?(u1,u2,…,um)∈Im,C1(u1,u2,…，um)≤C2(u1,u2,…,um)，則稱Copula函數(shù)C1

m維Copula函數(shù)的Freche-hoeffding上、下界分別為

(7)

Sklar定理：一個(gè)聯(lián)合分布函數(shù)H(y1,y2,…,ym)可分解為m個(gè)邊緣分布函數(shù)F1(y1),F2(y2),…,Fm(ym)和一個(gè)Copula函數(shù)，即

H(y1,y2,…,ym)=Cm(F1(y1),F2(y2),…,Fm(ym))

(8)

式中Cm(·)為一個(gè)m維的Copula 函數(shù)。如果F1(y1),F2(y2),…,Fm(ym)是連續(xù)的，則Copula函數(shù)唯一確定；否則要根據(jù)RandF1×RandF2×…×RandFm確定。相反，如果F1(y1),F2(y2),…,Fm(ym)為分布函數(shù)，那么函數(shù)H(y1,y2,…,ym)為一個(gè)聯(lián)合分布函數(shù)，其邊緣分布為F1(y1),F2(y2),…,Fm(ym)，通過(guò)Copula函數(shù)的密度函數(shù)cm(·)和邊緣分布函數(shù)F1(y1),F2(y2),…,Fm(ym)，可以方便地求出m維分布函數(shù)的密度函數(shù)如下

h(y1,y2,…,ym)=cm(F1(y1),F2(y2),…,

(9)

2.2 可靠性評(píng)估模型

根據(jù)Copula函數(shù)的理論基礎(chǔ)，將Copula函數(shù)引入到多性能退化量的可靠性評(píng)估模型，不難得出多性能參數(shù)產(chǎn)品的聯(lián)合密度函數(shù)為

f(y1,y2,…,ym)=cm(F1(y1),F2(y2),…,

(10)

聯(lián)合分布函數(shù)為

H(y1,y2,…,ym)=Cm(F1(y1),F2(y2),…,Fm(ym))

(11)

假設(shè)Ω=(v1,v2,…,vm,θ)為Copula 函數(shù)的參數(shù)向量，其中v1,v2,…,vm為各個(gè)性能參數(shù)的邊緣分布函數(shù)F1(y1),F2(y2),…,Fm(ym)的參數(shù)向量，θ為Copula函數(shù)的相關(guān)系數(shù)。根據(jù)精確極大似然估計(jì)[9]得似然函數(shù)為

(12)

將L(Ω)關(guān)于參數(shù)向量Ω極大化，得極大似然估計(jì)為

(13)

式(13)可以同時(shí)求得邊緣分布函數(shù)F1(y1),F2(y2),…,Fm(ym)的參數(shù)v1,v2,…,vm和θ，但當(dāng)維數(shù)較高時(shí)計(jì)算十分復(fù)雜，為了簡(jiǎn)化問(wèn)題，采用邊際分布的函數(shù)推斷方法(function inference method,FIM)對(duì)上式進(jìn)行參數(shù)估計(jì)，即二步極大似然估計(jì)法。首先利用最大似然估計(jì)法對(duì)各邊緣分布的參數(shù)v1,v2,…,vm進(jìn)行估計(jì)，并對(duì)參數(shù)θ進(jìn)行估計(jì)，其步驟如下：

1)對(duì)v1,v2,…,vm進(jìn)行估計(jì)得

(14)

Fm(ymj;vm);θ)

(15)

則可由Copula函數(shù)推出產(chǎn)品可靠度為

R(t)=Cm(F1(D1;v1),F2(D2;v2),…,Fm(Dm;vm);θ)

(16)

綜上，基于Copula函數(shù)的多性能退化模型的步驟為:

1)獲取或求取m組數(shù)據(jù)的邊緣分布函數(shù)；

2)選擇適當(dāng)?shù)腃opula函數(shù)形式并求取參數(shù)估計(jì)值；

3)根據(jù)產(chǎn)品可靠度公式求得產(chǎn)品可靠度。

3 實(shí)例應(yīng)用

以某型地空導(dǎo)彈為例，通過(guò)MATLAB仿真軟件進(jìn)行仿真計(jì)算。已知該型導(dǎo)彈可靠性評(píng)定主要由參數(shù)X1及參數(shù)X2共同決定，在自然條件下連續(xù)進(jìn)行20次測(cè)量的數(shù)據(jù)如表1。

表1 各參數(shù)的連續(xù)測(cè)試數(shù)據(jù)

其中參數(shù)X1的失效閾值為2.5 Ω，參數(shù)X2的失效閾值為10 Ω。為確定參數(shù)X1,X2的分布類型，首先給出其頻率直方圖如圖1。

從圖1難以準(zhǔn)確看出參數(shù)X1,X2的分布類型，但參數(shù)均類似于正態(tài)分布，對(duì)參數(shù)進(jìn)行正態(tài)檢驗(yàn)，引入正態(tài)分析圖對(duì)參數(shù)進(jìn)行檢驗(yàn)的圖2所示。

圖1 參數(shù)頻率直方圖

從圖2可以看出:參數(shù)X1及X2服從正態(tài)分布，且求得X1～N(2.2,0.004)，X2～N(42,8)。參數(shù)X1及參數(shù)X2的二元頻率直方圖，如圖3所示。

圖3 X1及X2的二元頻率直方圖

根據(jù)圖3及參數(shù)X1,X2之間的相關(guān)性，選取二元正態(tài)Copula[10]作為其連接函數(shù)，有

C(u1,u2;ρ)=Φρ(Φ-1(u1),Φ-1(u2))

(17)

C(u1,u2;ρ)=Φρ(Φ-1(u1),Φ-1(u2))=

(18)

將ρl=0.952 3代入式(18)可得各時(shí)刻的導(dǎo)彈可靠性數(shù)據(jù)如圖4所示。其隨著時(shí)間而降低，符合實(shí)際規(guī)律，同時(shí)表明模型的有效性。

圖4 導(dǎo)彈各時(shí)刻的可靠度

4 結(jié) 論

本文提出了基于Copula函數(shù)的多性能退化可靠性評(píng)估模型，以較高的準(zhǔn)確度描述了導(dǎo)彈多性能指標(biāo)間的關(guān)系及系統(tǒng)的可靠性退化特征。模型在確保評(píng)估結(jié)果準(zhǔn)確的基礎(chǔ)上，大幅減小了計(jì)算量，為導(dǎo)彈多性能退化可靠性模型研究提供了一種新方法，具有較高的推廣性和實(shí)用性。

[1] QJ3153—2002.導(dǎo)彈貯存可靠性設(shè)計(jì)技術(shù)指南[S].

[2] 段鵬程,王學(xué)奇,霍建成,等.基于多傳感器的機(jī)載導(dǎo)彈環(huán)境應(yīng)力監(jiān)測(cè)裝置[J].傳感器與微系統(tǒng),2015,34(9):111-114.

[3] 曾 健,陳俊芳.Copula函數(shù)在風(fēng)險(xiǎn)管理中的應(yīng)用研究—以上證A股與B股的相關(guān)結(jié)構(gòu)分析為例[J].當(dāng)代財(cái)經(jīng),2005(2):34-38.

[4] Kim I T,Itoh Y.Accelerated exposure tests as evaluation tool for estimating life of organic coatings on steel bridges[J].Corrosion Engineering,Science and Technology,2007,42(3):242-252.

[5] 劉景森，金 勇.一種基于吉布斯抽樣的MUSIC多維參數(shù)聯(lián)合估計(jì)算法[J].傳感器與微系統(tǒng),2008,27(6):62-65.

[6] 戴邵武，王克紅，錢儉學(xué).基于AKPSO算法的加速度計(jì)快速標(biāo)定方法[J].傳感器與微系統(tǒng),2015,34(2):69-72.

[7] 王 凱.導(dǎo)導(dǎo)彈武器系統(tǒng)貯存環(huán)境監(jiān)測(cè)及貯存可靠性評(píng)定方法研究[D].哈爾濱:哈爾濱理工大學(xué),2012.

[8] Nelson R B.An introduction to Copulas[M].Berlin Heidelberg:Springer,2006.

[9] 王 琦,馮新喜,杜欽峰.三維傳感器極大似然配準(zhǔn)算法[J]. 傳感器與微系統(tǒng),2007,26(11):54-56.

[10] 張世英.Copula理論及其在金融分析上的應(yīng)用[M].北京：清華大學(xué)出版社,2008.