馮海林，李秀秀

(西安電子科技大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院， 陜西 西安 710126)

0 引 言

隨著科學(xué)技術(shù)的快速發(fā)展，航空航天、電子工業(yè)、機械等領(lǐng)域出現(xiàn)了大量具有高可靠性與長壽命的產(chǎn)品，即在常應(yīng)力下很難觀測到產(chǎn)品失效,在加速壽命試驗中也很少甚至無失效數(shù)據(jù). 這給傳統(tǒng)的基于壽命數(shù)據(jù)的可靠性分析帶來了極大挑戰(zhàn). 在此背景下，基于性能退化的建模思想，為這些長壽命、高可靠性產(chǎn)品的可靠性分析提供了新的途徑[1-2]. 從產(chǎn)品性能退化的角度分析其失效情況，一般需要選擇一個與產(chǎn)品壽命或可靠性高度相關(guān)且可以測量的或從測量數(shù)據(jù)中可以提取的變量，稱其為性能可靠性特征量[3-4].刻畫產(chǎn)品性能退化所提取的特征量稱為性能退化量，性能退化量受產(chǎn)品使用環(huán)境等因素的影響，且隨時間的變化而變化，具有一定的隨機性.因此，建立合理的性能退化隨機過程模型是準確預(yù)測產(chǎn)品剩余壽命的關(guān)鍵環(huán)節(jié)之一. 基于標準Wiener過程所衍生的性能退化過程模型是近年來廣受關(guān)注且可描述多種典型產(chǎn)品性能退化的一類隨機過程，其一般形式為

X(t)=vΛ(t)+σB(Λ(t)),

(1)

其中，{X(t)：t≥0}表示產(chǎn)品性能的退化量.v是漂移參數(shù)，σ是擴散參數(shù)，且v和σ均為常數(shù).B(·)是標準的Wiener過程，Λ(·)通常表示時間尺度的單調(diào)增函數(shù)，且Λ(0)=0和B(0)=0.

根據(jù)以上隨機效應(yīng)Wiener過程模型和當(dāng)前的性能退化量，建立剩余壽命的退化模型，獲得剩余壽命的概率密度函數(shù)和可靠度函數(shù). 通過分析退化模型參數(shù)，可知模型的參數(shù)含有隨機效應(yīng)，對數(shù)似然函數(shù)比較復(fù)雜，很難最優(yōu)化. 因此，選用EM算法進行參數(shù)估計，再融合同類產(chǎn)品的退化數(shù)據(jù)，對該Wiener退化模型中的未知參數(shù)進行估計，并與已有模型進行對比，以說明本模型的實用性和精確性.

1 含隨機效應(yīng)的Wiener退化模型

令X(t)為產(chǎn)品在t時刻的性能退化量的測量值，建立隨機效應(yīng)Wiener退化模型為

X(t)=vΛ(t)+(ζv+a)B(Λ(t)),

(2)

由于產(chǎn)品制造工藝、制造材料的差異以及受工作環(huán)境的動態(tài)影響，同類型不同產(chǎn)品的性能退化過程存在差異. 為描述這種差異性，在退化模型中常采用將相關(guān)參數(shù)隨機化的方法，即將Wiener退化模型(1)中的漂移參數(shù)設(shè)成正態(tài)分布.

模型(2)與模型(1)相比，有σ=ζv+a，表示v和σ呈線性關(guān)系. 通過加入誤差項a，縮小了性能退化測量值與真實值之間的偏差，從而提高了產(chǎn)品可靠性評估的精確度. 該Wiener模型不僅體現(xiàn)了產(chǎn)品個體在同時具有高退化率與高波動性的性質(zhì)，而且也解決了產(chǎn)品的退化性能在測量時會產(chǎn)生誤差的問題；同時，大大降低了計算難度，更符合實際產(chǎn)品的退化數(shù)據(jù)，從而能更準確地預(yù)測產(chǎn)品的剩余壽命.

根據(jù)Wiener退化模型(1)，考慮漂移參數(shù)的隨機效應(yīng)，可知模型(2)中退化量X(t)的概率密度函數(shù)(PDF)為

(3)

2 基于Wiener退化過程的剩余壽命預(yù)測

本文的主要目的是利用新隨機效應(yīng)Wiener退化模型(2)描述產(chǎn)品的退化軌跡，最終預(yù)測產(chǎn)品的剩余壽命. 下面將由首達時的概念推導(dǎo)隨機效應(yīng)Wiener退化模型(2)計算的剩余壽命的分布形式. 為此，首先給出Wiener過程首達時的分布.

2.1 Wiener過程首達時的分布

給定產(chǎn)品的失效閾值ξ,產(chǎn)品的壽命T為性能退化量X(t)首次達到失效閾值ξ的時間，則

T=inf(t|X(t)≥ξ,t>0).

(4)

考慮模型(2)中的漂移參數(shù)μ和a是隨機變量，則產(chǎn)品的壽命T的根率密度函數(shù)(PDF)為

(5)

其相應(yīng)的概率分布函數(shù)為

(6)

因此，產(chǎn)品的可靠度函數(shù)可表示為

(7)

2.2 產(chǎn)品的剩余壽命預(yù)測

為預(yù)測某產(chǎn)品的剩余壽命，首先須建立此產(chǎn)品的剩余壽命預(yù)測模型. 假定特定產(chǎn)品在時刻tk和t+tk的性能退化量為X(tk)和X(t+tk). 按照Wiener過程獨立增量的性質(zhì)，可得

X(t+tk)=X(tk)+vΛ(t)+(ζv+a)B(Λ(t)).

(8)

產(chǎn)品的剩余壽命是指從當(dāng)前時刻到產(chǎn)品失效時刻的時間間隔. 如果已知當(dāng)前的性能退化量為X(tk)，根據(jù)壽命的定義可知產(chǎn)品在時刻tk的剩余壽命t為

M=inf(t|X(t+tk)≥ξ,t>0).

(9)

獲得產(chǎn)品剩余壽命的關(guān)鍵是找到剩余壽命的概率密度函數(shù)(PDF).由Wiener過程獨立增量性質(zhì)知,

M= inf(t|X(t+tk)≥ξ,t>0)=

inf(t|X(t+tk)-X(tk)≥ξ-X(tk),t>0)=

inf(t|X(t)≥ξ-X(tk),t>0),

(10)

則產(chǎn)品的剩余壽命t的PDF為

(11)

3 Wiener退化模型的參數(shù)估計

由退化數(shù)據(jù)估計參數(shù)

(12)

其中β為Λ(t;β)中未知參數(shù)的集合. 當(dāng)一個產(chǎn)品的性能退化量可觀察時，能夠更新v的分布. 假定可獲得N個實驗產(chǎn)品的退化數(shù)據(jù)，當(dāng)時間為(t1,t2,…,tn)′時觀察的退化量為x=(x1,x2,…,xn)′，令Δx=(Δx1,Δx2,…,Δxn)′，其中Δx1=x1和Δxj=xj-xj-1. 類似地，令λ1=Λ(t1)和λj=Λ(tj)-Λ(tj-1),1

(13)

其中，

令vi為第i個產(chǎn)品已實現(xiàn)的隨機效應(yīng). 以vi為條件，第i個產(chǎn)品的退化量{Xi(t);vi,t>0}有獨立的增量. ΔXi=(ΔXi1,ΔXi2,…,ΔXini)為退化量增量，其中ΔXij=X(tij)-X(tij-1),j=1,2,…,ni，ΔXi1=X(ti1). Δxi是增量的觀察值，N個產(chǎn)品的觀察數(shù)集D={Δxi,i=1,2,…,N}. 由對數(shù)似然函數(shù)得到極大似然估計，式(13)可改寫為式(14).

(14)

因式中含有隨機效應(yīng)，對數(shù)似然函數(shù)較復(fù)雜，很難最優(yōu)化. 因此，可用EM算法估計新隨機效應(yīng)Wiener退化模型中的未知參數(shù). 給定參數(shù)估計θ(m)，經(jīng)EM算法反復(fù)迭代Eθ(m)[l(θ)]后，用θ更新θ(m).

EM算法的求解分2步：

E-step給定觀測數(shù)據(jù)集D={Δxi,i=1,2,…,N}和初始參數(shù)值θ(0)，計算對數(shù)似然函數(shù)l(θ)關(guān)于未知參數(shù)θ的期望. 定義對數(shù)似然函數(shù)的期望Q(θ|θ(m))=E[l(θ;D)|θ(m),D],θ(m)表示m次迭代的估計參數(shù) .vi關(guān)于Δxi和θ的條件分布為

(15)

由式 (15)可改寫為式(14)，有

(16)

第2部分為

(17)

利用Q2(θ|θ(m))對ζ2求導(dǎo)，且令導(dǎo)數(shù)為0，式(18)改寫為式(17)：

(18)

將式(18)代入式(17)，得到

(19)

從β開始最大化式(18)，得到β(m+1)，且將結(jié)果代入式(17)，得到(ζ2)(m+1).

4 數(shù)據(jù)分析

基于本文的Wiener退化模型數(shù)據(jù)分析，以文獻[8]中2024-T351鋁合金樣本的疲勞裂紋擴展數(shù)據(jù)為例(在同一恒定振幅下獲得數(shù)據(jù))，所有樣品的初始裂紋長度相同；每個樣品的裂紋長度增量是在10 000個時間周期后測量的，每間隔5 000個時間周期測量1次，直到40 000個時間周期止. 實驗數(shù)據(jù)見表1[8]，得到的裂紋增長路徑圖見圖1.

表1 5種Wiener過程模型的MLE和AIC

圖1 30個樣本的裂紋增長路徑Fig.1 Fatigue crack growth paths of 30 testing units

WU等[8]通過縱向數(shù)據(jù)分析技術(shù)擬合數(shù)據(jù)，本文使用Wiener過程分析此組數(shù)據(jù). 研究中，簡單模型式(1)記為SM，其中參數(shù)v和σ為固定常數(shù)；由YE等[7]提出的隨機效應(yīng)模型記為RDV，其漂移參數(shù)v～N(μ，w2)，擴散參數(shù)σ=ζv；由PENG 等[5]提出的隨機效應(yīng)模型記為RD，其v～N(μ，w2)；由WANG[6]提出的隨機效應(yīng)模型記為RDRV，其漂移參數(shù)[v|σ]～N(μ,θσ2)，擴散參數(shù)σ-2～Ga(r,δ)；本文的隨機效應(yīng)模型記為NRDV,其漂移參數(shù)v～N(μ，w2)，擴散參數(shù)σ=ζv+a；對這幾種模型進行參數(shù)估計和AIC值計算.

通常認為裂紋增長遵循冪律，表示為Λ(t)=tb. 表1為SM(固定效應(yīng))、RDV(隨機效應(yīng)模型)、RD(單一效應(yīng)模型)、RDRV(混合效應(yīng)模型)和NRDV(新隨機效應(yīng)模型)5種模型參數(shù)的極大似然估計及對應(yīng)的AIC(Akaike information criterion).

AIC=2k-2 lnl，

其中，k為參數(shù)的數(shù)量，l為模型似然函數(shù)的最大值. AIC值越小，表示模型越符合實際. 由表1 知，“NRDV”最符合實際數(shù)據(jù)，因此，由該模型預(yù)測的剩余壽命更精確.

圖2 參數(shù)和之間的散點圖Fig.2 Scatter plot of the estimated drift parameter and estimated volatility parameter .

從圖2中可以看出，估計的漂移參數(shù)和擴散參數(shù)呈線性分布，且一次擬合系數(shù)為 0.583 6. 說明當(dāng)產(chǎn)品有更高的退化率時，其退化路徑會隨退化率的變化而波動. 本文提出的模型可以更好地解釋這一現(xiàn)象.

圖3 參數(shù)vi的標準正態(tài)Q-Q圖Fig.3 Standard normal Q-Q plot for parameter

5 結(jié) 論

高退化率產(chǎn)品也具有高波動性，這是一種常見現(xiàn)象. 由于均值和方差是2個互相獨立的參數(shù)，以前研究的Wiener退化過程不易描繪上述問題. 為了克服這一缺陷，本文提出了新隨機效應(yīng)的Wiener過程模型，即在簡單的Wiener退化模型上，引入隨機效應(yīng)，對v施加一個正態(tài)分布，并將擴散參數(shù)與漂移參數(shù)設(shè)定成線性關(guān)系， 此模型較現(xiàn)有隨機效應(yīng)維納過程模型(擴散參數(shù)和漂移參數(shù)都相互獨立)更符合實際情況. 對新的隨機效應(yīng)模型進行統(tǒng)計推斷，并用EM算法估計模型參數(shù)，驗證了本文模型的有效性和準確性.

本模型可用于退化試驗，以評估產(chǎn)品的可靠性，也可預(yù)測已知當(dāng)前退化量產(chǎn)品的剩余壽命. 結(jié)果表明，本模型更接近實際應(yīng)用. 因此，利用本模型估計的可靠性和預(yù)測的剩余使用壽命更準確. 同時，本模型還具有擴展性，而現(xiàn)有隨機效應(yīng)模型中的假設(shè)是有所限制的. 需要指出的是： 除了正態(tài)分布外，截斷正態(tài)分布也常用來表示隨機效應(yīng)分布.