毛澤龍, 王治華,*, 吳 瓊, 劉成瑞

(1. 北京航空航天大學(xué)航空科學(xué)與工程學(xué)院, 北京 100083; 2. 北京空間飛行器總體設(shè)計(jì)部, 北京 100080;3. 北京控制工程研究所研發(fā)中心, 北京 100080)

0 引 言

目前,針對高可靠長壽命產(chǎn)品的可靠性分析,性能退化建模方法已廣泛應(yīng)用。為簡化模型,早期研究大多在產(chǎn)品退化過程僅具備單一退化規(guī)律的前提下開展。工程實(shí)際中,一方面,產(chǎn)品自身物理機(jī)理或運(yùn)行環(huán)境、負(fù)載情況等的動態(tài)變化,會導(dǎo)致退化過程呈現(xiàn)出兩階段甚至多階段的演變特點(diǎn)[1-3];另一方面,由于功能多樣性和組成復(fù)雜性,產(chǎn)品往往具有多個(gè)性能指標(biāo),且這些指標(biāo)間存在耦合關(guān)系[4-6]。此時(shí),單一指標(biāo)單一階段的性能退化模型往往不能準(zhǔn)確反映高可靠長壽命產(chǎn)品的真實(shí)退化規(guī)律。針對此類產(chǎn)品,如何建立合理退化模型準(zhǔn)確描述多階段、多指標(biāo)退化規(guī)律成為相關(guān)研究的熱點(diǎn)和難點(diǎn)。

針對單一指標(biāo)多階段性能退化問題,基于隨機(jī)過程的方法能較好描述退化過程中由產(chǎn)品內(nèi)部失效機(jī)理及外界環(huán)境因素導(dǎo)致的時(shí)變不確定性,已廣泛應(yīng)用[7]。其中,Wiener過程由于能描述非單調(diào)退化過程以及良好的數(shù)學(xué)性質(zhì)和物理解釋,應(yīng)用最為普遍[8-9]。文獻(xiàn)[10]針對液力耦合器退化過程呈現(xiàn)先慢后快的特征,以振動幅值為性能指標(biāo),基于Wiener過程建立了兩階段退化模型。文獻(xiàn)[11]針對液力耦合器第一階段為單調(diào)退化,第二階段為非單調(diào)退化的特點(diǎn),進(jìn)一步提出了基于Gamma過程和Wiener過程的兩階段退化模型。文獻(xiàn)[12]針對電容器容量退化特征,建立三階段Wiener過程評估其可靠性。針對產(chǎn)品退化過程中存在的非線性規(guī)律,文獻(xiàn)[13]基于Wiener過程建立了兩階段非線性模型;文獻(xiàn)[14]進(jìn)一步建立了發(fā)動機(jī)的多階段非線性退化模型。

針對多指標(biāo)聯(lián)合退化問題,指標(biāo)間相關(guān)性的合理描述至關(guān)重要。目前考慮指標(biāo)相關(guān)性的建模方法主要有以下幾類:基于多維隨機(jī)過程(如多維Wiener過程)建模[15-18];將多指標(biāo)融合為單一指標(biāo)后建模[19-20];基于Copula函數(shù)的建模[4,21-29]。其中,基于多維隨機(jī)過程的建模方法適用于各指標(biāo)退化規(guī)律相近的情況;對于第二種思路,融合原始指標(biāo)后得到的新指標(biāo),難以明確其物理意義及相應(yīng)失效閾值,限制了其應(yīng)用范圍;而基于Copula函數(shù)的建模方法能夠有效克服以上問題,因此得到了廣泛應(yīng)用。文獻(xiàn)[21]基于Wiener過程和Copula函數(shù)建立雙指標(biāo)退化模型,各指標(biāo)的邊緣分布通過Wiener過程描述,指標(biāo)間的相關(guān)關(guān)系由Copula函數(shù)表征。文獻(xiàn)[22]進(jìn)一步研究退化過程具有個(gè)體差異時(shí)的雙指標(biāo)建模問題。文獻(xiàn)[23]以某型繼電器為研究對象,提出了基于Copula函數(shù)和Wiener過程的雙指標(biāo)加速退化數(shù)據(jù)建模方法,并建立了基于貝葉斯馬爾可夫鏈蒙特卡羅的模型參數(shù)估計(jì)方法,擴(kuò)充了雙指標(biāo)退化模型的使用范圍。文獻(xiàn)[4,24-26]進(jìn)一步針對指標(biāo)數(shù)量更多的情況開展了研究?？紤]到某些產(chǎn)品各指標(biāo)退化過程是單調(diào)的,文獻(xiàn)[27-28]采用逆高斯過程和Copula函數(shù)建立雙指標(biāo)退化模型。

綜上可見,以上有關(guān)階段性特征的研究均圍繞單一指標(biāo)的退化過程展開,無法描述多指標(biāo)聯(lián)合退化時(shí)的情況;有關(guān)多指標(biāo)建模研究均未考慮指標(biāo)退化的階段性特征,可能導(dǎo)致可靠度評估結(jié)果出現(xiàn)偏差。對于產(chǎn)品退化過程中耦合性與階段性規(guī)律并存時(shí)可靠性分析問題有待進(jìn)一步研究。文獻(xiàn)[29]以飛機(jī)艙門鎖為研究對象,采用指標(biāo)間相關(guān)性描述退化過程的耦合性,對退化過程同時(shí)存在耦合性與階段性的建模問題進(jìn)行了探究。

本文通過Copula函數(shù)描述指標(biāo)間相關(guān)性,引入變點(diǎn)描述退化階段性,考慮不同階段指標(biāo)間耦合性規(guī)律變化,建立了雙指標(biāo)階段性退化模型,并給出了可靠性評估方法。另外,針對退化過程包含兩個(gè)變點(diǎn)且不在同一位置,而難以估計(jì)的問題,提出了一種考慮雙指標(biāo)變點(diǎn)的模型參數(shù)整體估計(jì)方法。最后,通過實(shí)例分析,驗(yàn)證了本文方法的建模合理性與可靠性分析有效性。

1 退化規(guī)律分析與模型引出

工程實(shí)際中,一些產(chǎn)品退化過程呈現(xiàn)多指標(biāo)和階段性規(guī)律并存的特點(diǎn)。例如,飛機(jī)艙門鎖[29]關(guān)鍵部分為連桿機(jī)構(gòu),連桿之間、連桿與基座間均通過軸和軸套連接,構(gòu)成了多個(gè)連接點(diǎn)。在使用中,連接點(diǎn)軸套會逐漸磨損,當(dāng)任一連接點(diǎn)軸套磨損量超出對應(yīng)的失效閾值,則艙門鎖失效。通過對艙門鎖的失效機(jī)理和退化數(shù)據(jù)分析,發(fā)現(xiàn)其中兩處連接點(diǎn)軸套的磨損最為嚴(yán)重,將這兩處軸套的磨損量作為艙門鎖的兩個(gè)性能指標(biāo),記為指標(biāo)1和指標(biāo)2,其性能退化測試數(shù)據(jù)如圖1所示。

圖1 艙門鎖性能退化數(shù)據(jù)Fig.1 Degradation data of cabin door lock mechanism

從圖1中可見,指標(biāo)1和指標(biāo)2的退化呈現(xiàn)明顯的兩階段特征,這是由于前期潤滑劑充足,軸與軸套間潤滑良好,軸套磨損速率較慢;隨著潤滑劑消耗,當(dāng)消耗到某一水平,軸套磨損速率迅速增大。同時(shí)需注意,指標(biāo)1變點(diǎn)在0.2×105次運(yùn)動循環(huán)附近,指標(biāo)2變點(diǎn)在0.3×105次運(yùn)動循環(huán)附近,兩指標(biāo)變點(diǎn)不在同一位置。此外,由于艙門鎖的結(jié)構(gòu)特點(diǎn),兩連接點(diǎn)的磨損也相互影響,即指標(biāo)1和指標(biāo)2之間還存在耦合關(guān)系。

綜上可見,此類產(chǎn)品退化特點(diǎn)可總結(jié)如下:考慮產(chǎn)品具有兩個(gè)性能指標(biāo),令Xi(t)表示第i個(gè)指標(biāo)的退化過程,其中i=1,2;當(dāng)任一指標(biāo)首次達(dá)到對應(yīng)失效閾值時(shí),則產(chǎn)品發(fā)生失效;兩指標(biāo)退化過程均呈現(xiàn)階段性,且在很多情況下,不同階段內(nèi)指標(biāo)相關(guān)性強(qiáng)弱不一致(該結(jié)論將在實(shí)例分析中說明)。根據(jù)兩指標(biāo)變點(diǎn)τ1、τ2間位置關(guān)系,問題可分為兩類:考慮最一般的情況,當(dāng)τ1≠τ2時(shí),不妨令τ1<τ2,變點(diǎn)τ1、τ2將產(chǎn)品退化過程分為3段：(0,τ1]、(τ1,τ2]、(τ2,+∞),各段內(nèi)指標(biāo)退化規(guī)律與耦合規(guī)律均可以不相同;當(dāng)τ1=τ2時(shí),產(chǎn)品退化過程為兩階段,可視為前述一般情況的簡化。

此類產(chǎn)品退化過程中耦合性和階段性規(guī)律并存,帶來了以下挑戰(zhàn):一是如何建立合理的退化模型表征退化過程中不同階段雙指標(biāo)間耦合性規(guī)律;二是由于產(chǎn)品兩指標(biāo)退化過程均包含變點(diǎn),且兩變點(diǎn)不在同一位置,如何建立有效的模型參數(shù)估計(jì)方法;最后,由于退化過程復(fù)雜,可靠度解析形式難以獲得,如何基于退化模型有效評估產(chǎn)品可靠性。

2 雙指標(biāo)階段性退化建模

2.1 邊緣退化建模

對于兩個(gè)性能指標(biāo)均呈現(xiàn)兩階段退化規(guī)律的情況,對指標(biāo)i,基于Wiener過程建立如下兩階段邊緣退化模型:

Xi(t)=[μi1Λi1(t)+σi1B(Λi1(t))]I(0,τi](t)+
[Xi(τi)+μi2Λi2(t-τi)+σi2B(Λi2(t-τi))]I(τi,∞)(t)

(1)

式中:τi為指標(biāo)i的變點(diǎn)(即在時(shí)刻τi,指標(biāo)i退化過程由第一階段進(jìn)入第二階段);μi1和μi2分別為指標(biāo)i第一和第二階段的漂移系數(shù),表征各階段的退化速率;σi1和σi2分別為指標(biāo)i第一和第二階段的擴(kuò)散系數(shù),表征各階段分散性;I(0,τi](t)、I(τi,∞)(t)為示性函數(shù);B(·)為標(biāo)準(zhǔn)Wiener過程;Λi1(t)=tγi1、Λi2(t-τi)=(t-τi)γi2為第一和第二階段的時(shí)間尺度變換函數(shù),分別表征指標(biāo)i在各階段退化的非線性特征;指標(biāo)i的邊緣退化模型參數(shù)可記為Θi=(μi1,μi2,σi1,σi2,γi1,γi2,τi),i=1,2。

(2)

則退化增量ΔXij在變點(diǎn)τi前后服從不同參數(shù)的正態(tài)分布可表示為

(3)

退化增量ΔXij的累積分布函數(shù)和概率密度函數(shù)分別為

(4)

(5)

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

2.2 基于Copula的雙指標(biāo)階段性退化建模

根據(jù)常用假設(shè)[21],本文通過Copula函數(shù)描述指標(biāo)間的相關(guān)性,在不同時(shí)間間隔內(nèi),指標(biāo)1與指標(biāo)2退化增量間的相關(guān)性可忽略,即僅考慮同一時(shí)間間隔內(nèi)兩指標(biāo)退化增量間的相關(guān)性;由此在邊緣退化模型基礎(chǔ)上,構(gòu)建雙指標(biāo)階段性退化模型。根據(jù)Sklar’s理論,兩指標(biāo)退化增量ΔX1j、ΔX2j的聯(lián)合分布函數(shù)F(ΔX1j,ΔX2j)可由各自邊緣分布函數(shù)F1(ΔX1j)、F2(ΔX2j)和Copula函數(shù)表示,即

F(ΔX1j,ΔX2j)=C(F1(ΔX1j),F2(ΔX2j);θ)

(6)

式中:C(·)為Copula函數(shù),表征指標(biāo)退化增量間的相關(guān)關(guān)系;θ為Copula函數(shù)參數(shù),與指標(biāo)退化增量間的相關(guān)性大小有關(guān)。

兩指標(biāo)的變點(diǎn)τ1、τ2將測試時(shí)間段分成了3段,考慮到各段內(nèi)指標(biāo)間的相關(guān)性大小可能不一致,θ在各段內(nèi)取值也不相同,分別為θ1、θ2、θ3,即

θ=θ1I(0,τ1](tj)+θ2I(τ1,τ2](tj)+θ3I(τ2,∞)(tj)

(7)

針對兩個(gè)指標(biāo)的情況,目前應(yīng)用較廣的Copula函數(shù)有Gaussian Copula、 Gumbel Copula、Clayton Copula、Frank Copula等,這幾種Copula函數(shù)具有不同的耦合性描述效果。本文采用在退化領(lǐng)域已廣泛使用的Frank Copula函數(shù)描述指標(biāo)增量間的耦合性。Frank Copula函數(shù)及其密度函數(shù)具體形式見文獻(xiàn)[30]。綜合考慮指標(biāo)退化的階段性與指標(biāo)間的相關(guān)性,建立雙指標(biāo)階段性退化模型如下:

(8)

模型參數(shù)可記為Θ=(Θ1,Θ2,θ1,θ2,θ3)。

3 可靠性分析

確定模型參數(shù)估計(jì)值是利用雙指標(biāo)階段性模型進(jìn)行可靠性分析的前提。模型待估參數(shù)為Θ=(Θ1,Θ2,θ1,θ2,θ3)。由于參數(shù)較多,難以通過直接優(yōu)化對數(shù)似然函數(shù)獲得模型參數(shù)估計(jì)值,且模型參數(shù)中包含兩個(gè)變點(diǎn),進(jìn)一步給參數(shù)估計(jì)帶來困難。針對以上問題,本文提出了一種考慮雙指標(biāo)變點(diǎn)的模型參數(shù)整體估計(jì)方法。

假設(shè)兩個(gè)性能指標(biāo)的變點(diǎn)τ1、τ2分別位于測試時(shí)刻tw1、tw2(tw1

f(ΔX1j,ΔX2j)=c(F1(ΔX1j),F2(ΔX2j);θ)·
f1(ΔX1j)·f2(ΔX2j)

(9)

式中:c(·)為Frank Copula函數(shù)對應(yīng)的密度函數(shù)。

對于(t1,t2,…，tw1)時(shí)刻的測試數(shù)據(jù),由式(9)得到同一時(shí)刻各指標(biāo)退化增量的聯(lián)合密度函數(shù)為

(10)

進(jìn)一步由不同時(shí)刻退化增量獨(dú)立,建立(t1,t2,…，tw1)時(shí)刻的退化數(shù)據(jù)的似然函數(shù)為

(11)

相應(yīng)的對數(shù)似然函數(shù)為

(12)

類似地,(tw1+1,tw1+2,…,tw2)、(tw2+1,tw2+2,…,tn)時(shí)刻的測試數(shù)據(jù)的對數(shù)似然函數(shù)分別為

(13)

(14)

綜上,由式(12)～式(14),m個(gè)試樣在時(shí)間段[t1,tn]的對數(shù)似然函數(shù)可表示為

lnL(Θ)=lnL1+lnL2+lnL3

(15)

建立對數(shù)似然函數(shù)后,根據(jù)以下步驟估計(jì)模型參數(shù)。

步驟 1定義如下關(guān)于tw1、tw2的函數(shù)。

log-LF(tw1,tw2)=lnL(Θ)max|τ1=tw1,τ2=tw2

(16)

式中：log-LF(tw1,tw2)為假定變點(diǎn)τ1、τ2位于測試時(shí)刻tw1、tw2時(shí),對數(shù)似然函數(shù)lnL(Θ)的最大值。其中對lnL(Θ)的優(yōu)化基于遺傳算法實(shí)現(xiàn),本文通過調(diào)用Matlab工具箱的ga函數(shù)優(yōu)化lnL(Θ)。

步驟 2確定tw1、tw2可能的取值組合。由于tw1、tw2滿足t1

步驟 3計(jì)算函數(shù)log-LF(tw1,tw2)在上述取值組合對應(yīng)的值。取log-LF(tw1,tw2)最大時(shí)對應(yīng)的取值組合作為變點(diǎn)估計(jì)值,對應(yīng)的其他模型參數(shù)取值作為模型參數(shù)最優(yōu)估計(jì)。

得到模型參數(shù)估計(jì)值后,可結(jié)合退化模型及失效閾值計(jì)算產(chǎn)品可靠度。對于兩指標(biāo)聯(lián)合退化的產(chǎn)品,當(dāng)其任一指標(biāo)退化量首次達(dá)到失效閾值時(shí),即認(rèn)為其發(fā)生失效。具體地,設(shè)產(chǎn)品兩指標(biāo)對應(yīng)失效閾值分別為D1、D2,產(chǎn)品可靠度可定義如下:

(17)

由于產(chǎn)品退化過程復(fù)雜,可靠度的解析形式難以獲得,采用蒙特卡羅方法計(jì)算可靠度。蒙特卡羅方法的基本思想是應(yīng)用隨機(jī)抽樣技術(shù)得到符合待求解分布的樣本,通過對樣本值的觀察與統(tǒng)計(jì),實(shí)現(xiàn)對未知分布的分析。

根據(jù)雙指標(biāo)階段性退化模型及模型參數(shù)估計(jì)值,開展蒙特卡羅仿真,生成樣本容量為N的仿真退化數(shù)據(jù);統(tǒng)計(jì)各樣本失效時(shí)間Th,其中h=1,2,…,N;統(tǒng)計(jì)t時(shí)刻未失效樣本個(gè)數(shù)N′(t),即滿足Th>t的樣本個(gè)數(shù),則t時(shí)刻可靠度為R(t)=N′(t)/N。其中,生成仿真數(shù)據(jù)具體包括以下步驟。

步驟 1設(shè)定仿真時(shí)間和步長等參數(shù),為了確保仿真樣本在仿真期間發(fā)生失效,仿真時(shí)間應(yīng)盡量長一些。設(shè)仿真時(shí)間步長為Δt,仿真時(shí)刻數(shù)為M,則總仿真時(shí)長為MΔt。

步驟 2從退化增量聯(lián)合分布C(F1(ΔX1j),F2(ΔX2j);θ)生成仿真時(shí)間段[(s-1)Δt,sΔt]內(nèi)指標(biāo)退化增量[ΔX1(sΔt),ΔX2(sΔt)]T,其中s=1,2,…,M。

步驟 3重復(fù)步驟2,得到仿真時(shí)長內(nèi),各指標(biāo)仿真退化增量,記為

步驟 4對退化增量求和得到仿真樣本在各仿真時(shí)刻指標(biāo)退化量:

步驟 5重復(fù)步驟2～步驟4即可得到樣本容量為N的仿真退化數(shù)據(jù)。

4 實(shí)例分析

為了驗(yàn)證本文方法的有效性,采用該方法針對艙門鎖實(shí)例開展分析。將本文模型記為M0,如式(8)所示。為體現(xiàn)本文方法適用性,與如下兩種參考方法進(jìn)行對比:① 現(xiàn)有僅考慮指標(biāo)間耦合性的雙指標(biāo)模型,記為M1,模型具體形式見文獻(xiàn)[21],以說明在雙指標(biāo)退化中考慮階段性的意義;② 實(shí)例數(shù)據(jù)來源于文獻(xiàn)[29],該文獻(xiàn)針對艙門鎖的退化特點(diǎn)建立了一種分析方法,將文獻(xiàn)[29]方法記為M2。

采用本文模型對艙門鎖性能退化數(shù)據(jù)進(jìn)行退化建模和參數(shù)估計(jì),得到對數(shù)似然函數(shù)log-LF(tw1,tw2)隨變點(diǎn)tw1、tw2變化關(guān)系如圖2所示。

圖2 log-LF(tw1,tw2)隨變點(diǎn)變化圖Fig.2 log-LF(tw1, tw2) variation with change points tw1 and tw2

表1 本文方法M0參數(shù)估計(jì)結(jié)果

由表1可見,不論對于指標(biāo)1還是指標(biāo)2,第一、二階段模型參數(shù)都有明顯差異,說明同一指標(biāo)不同階段退化速率及分散性規(guī)律發(fā)生變化;此外,不同階段Copula參數(shù)的不同說明不同階段指標(biāo)間相關(guān)性規(guī)律存在差異。

為了對比各模型擬合效果,表2列出了3種模型的對數(shù)似然函數(shù)值和AIC(akaike information criterion)值。其中AIC值定義如下:

AIC=-2(log-LF)+2n′

(18)

式中:n′為未知模型參數(shù)個(gè)數(shù)。AIC值越小,模型擬合能力更優(yōu)。

表2 log-LF值和AIC值對比結(jié)果

由表2可見,本文模型M0的log-LF和AIC值均優(yōu)于參考模型,這是由于本文方法同時(shí)考慮產(chǎn)品退化過程中的階段性與相關(guān)性規(guī)律,能夠合理描述退化過程。

為了更為直觀說明3種模型擬合效果及參數(shù)估計(jì)準(zhǔn)確性,圖3給出了平均退化曲線對比結(jié)果。從圖3中可知,對于指標(biāo)1和指標(biāo)2,本文方法M0給出的平均退化曲線能夠有效描述樣本均值的變化規(guī)律。

圖3 平均退化曲線對比結(jié)果Fig.3 Comparision results of mean degradation curves

進(jìn)一步,基于M0、M1和M2對艙門鎖進(jìn)行可靠性評估。指標(biāo)1和指標(biāo)2失效閾值分別為D1=1.6、D2=0.9。圖4展示了基于各模型得到的產(chǎn)品可靠度曲線,可見相較于本文模型M0,僅考慮兩指標(biāo)耦合性的模型M1前期結(jié)果偏保守,后期結(jié)果偏危險(xiǎn);模型M2得到的可靠性評估結(jié)果最為保守。

圖4 可靠度曲線對比結(jié)果Fig.4 Comparative results of reliability curves

通過以上對比說明本實(shí)例中,當(dāng)產(chǎn)品具有雙指標(biāo)聯(lián)合、階段性退化特點(diǎn)時(shí),須綜合考慮兩方面特點(diǎn)進(jìn)行合理建模,方能開展有效的可靠性分析。

5 結(jié) 論

針對產(chǎn)品同時(shí)具有雙指標(biāo)耦合性與階段性退化規(guī)律,本文提出了一種雙指標(biāo)階段性退化模型,并建立了可靠性分析方法。針對兩指標(biāo)變點(diǎn)不在同一位置的一般情況,本文建立了一種考慮雙指標(biāo)變點(diǎn)的模型參數(shù)整體估計(jì)方法。通過實(shí)例分析,驗(yàn)證本文方法的有效性。結(jié)論如下:

(1) 雙指標(biāo)耦合退化且存在階段性規(guī)律時(shí),不同階段指標(biāo)間的耦合性強(qiáng)弱可能不一致,建立退化模型時(shí)應(yīng)予以考慮;

(2) 本文方法能有效表征退化過程中兩指標(biāo)的耦合性與階段性規(guī)律,可靠性評估結(jié)果更為準(zhǔn)確可信,符合工程實(shí)際需求。