吳曉輝, 蔡忠義, 李全祥
(1.中國(guó)飛行試驗(yàn)研究院,陜西 西安 710089; 2.空軍工程大學(xué) 裝備管理與安全工程學(xué)院,陜西 西安 710051)
融合內(nèi)外場(chǎng)退化數(shù)據(jù)的可靠性評(píng)估方法
吳曉輝1, 蔡忠義2, 李全祥1
(1.中國(guó)飛行試驗(yàn)研究院,陜西 西安 710089; 2.空軍工程大學(xué) 裝備管理與安全工程學(xué)院,陜西 西安 710051)
針對(duì)性能退化過(guò)程服從Wiener過(guò)程的產(chǎn)品,文章運(yùn)用貝葉斯統(tǒng)計(jì)推斷法,提出了一種融合內(nèi)場(chǎng)加速退化試驗(yàn)(accelerated degradation test,ADT)與外場(chǎng)退化信息的可靠性評(píng)估方法。考慮到內(nèi)外場(chǎng)應(yīng)力環(huán)境差異,提出了基于修正系數(shù)的Wiener過(guò)程雙參數(shù)修正模型;建立了退化數(shù)據(jù)模型,得到各應(yīng)力下分布參數(shù)估計(jì)值;將各加速應(yīng)力下分布參數(shù)估計(jì)值折算到正常工作應(yīng)力下,構(gòu)成未知參數(shù)先驗(yàn)分布的數(shù)據(jù)樣本;構(gòu)建了外場(chǎng)退化數(shù)據(jù)下未知參數(shù)的后驗(yàn)分布函數(shù),采用馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)方法,得到參數(shù)的后驗(yàn)分布均值。實(shí)例分析結(jié)果驗(yàn)證了該文方法的正確性。
Wiener過(guò)程;退化數(shù)據(jù)融合;貝葉斯推斷;應(yīng)力環(huán)境差異
高可靠產(chǎn)品在外場(chǎng)正常工作應(yīng)力環(huán)境下很難獲取足夠多的失效數(shù)據(jù),導(dǎo)致傳統(tǒng)基于壽命數(shù)據(jù)的評(píng)估方法不能適用,但在外場(chǎng)通過(guò)測(cè)量產(chǎn)品關(guān)鍵性能參數(shù)隨時(shí)間的變化情況,收集外場(chǎng)退化信息,可以外推出產(chǎn)品壽命與可靠性水平[1]。同時(shí),為了豐富產(chǎn)品退化數(shù)據(jù)樣本、提高評(píng)估精度,可以借助于同類產(chǎn)品在加速退化試驗(yàn)(accelerated degradation test,ADT)中的性能退化數(shù)據(jù),但是從實(shí)驗(yàn)室收集到的數(shù)據(jù)與外場(chǎng)實(shí)測(cè)數(shù)據(jù)之間存在一定的差異。若將產(chǎn)品外場(chǎng)退化信息與同類產(chǎn)品ADT退化信息進(jìn)行融合評(píng)估,則需解決如下問(wèn)題:
(1) 考慮外場(chǎng)與實(shí)驗(yàn)室之間的應(yīng)力環(huán)境差異[2]。文獻(xiàn)[2]在產(chǎn)品退化過(guò)程服從Wiener過(guò)程的場(chǎng)合下,提出了基于廣義加速模型的退化過(guò)程參數(shù)修正模型來(lái)度量外場(chǎng)與實(shí)驗(yàn)室之間的應(yīng)力環(huán)境差異。
(2) ADT中同類產(chǎn)品個(gè)體之間的退化差異[3-8]。文獻(xiàn)[5]對(duì)Wiener過(guò)程的雙系數(shù)都進(jìn)行了隨機(jī)化處理,假定了擴(kuò)散系數(shù)服從逆伽瑪分布的前提下漂移系數(shù)服從正態(tài)分布,但假定的分布類型未經(jīng)分布假設(shè)檢驗(yàn),可能出現(xiàn)與實(shí)際退化數(shù)據(jù)擬合不一致的情況;文獻(xiàn)[7]認(rèn)為個(gè)體退化速率之間的差異主要在于漂移系數(shù),僅將Wiener過(guò)程的漂移系數(shù)看作服從正態(tài)分布的隨機(jī)變量,計(jì)算相對(duì)簡(jiǎn)便,但是忽略了擴(kuò)散系數(shù)對(duì)性能退化的影響。
(3) ADT與外場(chǎng)退化信息之間的融合方法。這類研究的主要解決方法是采用貝葉斯方法對(duì)多源信息進(jìn)行統(tǒng)計(jì)推斷,建立未知參數(shù)的估計(jì)模型,運(yùn)用最大期望法(expectation maximization,EM)[9-10]或馬爾科夫鏈蒙特卡洛(Markov chain Monte Carlo,MCMC)方法[11],得到未知參數(shù)估計(jì)值。如文獻(xiàn)[10]采用加速因子,將加速應(yīng)力下樣本數(shù)據(jù)折算到正常工作應(yīng)力下,并以此推斷出未知參數(shù)的先驗(yàn)分布并加以檢驗(yàn),使得未知參數(shù)先驗(yàn)分布具有良好的擬合適應(yīng)性。
因此,本文針對(duì)退化過(guò)程服從Wiener過(guò)程的產(chǎn)品,提出了融合ADT信息與外場(chǎng)退化信息的評(píng)估方法。首先,提出考慮外場(chǎng)差異的分布參數(shù)修正模型;其次,分別建立了步進(jìn)ADT數(shù)據(jù)與外場(chǎng)退化數(shù)據(jù)模型及分布參數(shù)估計(jì)函數(shù);然后,確定了未知參數(shù)的先驗(yàn)分布,建立了未知參數(shù)后驗(yàn)分布函數(shù);最后,通過(guò)實(shí)例分析驗(yàn)證了該方法的正確性。
假設(shè)產(chǎn)品在t時(shí)刻的性能退化量為X(t),X(t)為一元連續(xù)隨機(jī)過(guò)程且服從Wiener過(guò)程,可記為:
X(t)=μt+σB(t)+x0
(1)
其中,μ為漂移系數(shù);σ為擴(kuò)散系數(shù);B(t)為標(biāo)準(zhǔn)布朗運(yùn)動(dòng);x0為性能退化初值。
對(duì)于加速試驗(yàn),一般認(rèn)為Wiener過(guò)程的漂移系數(shù)μ與應(yīng)力S(本文應(yīng)力只討論大小)有關(guān),而擴(kuò)散系數(shù)σ與S無(wú)關(guān),則加速模型一般可以表示為:
μ=exp[a+bφ(S)]
(2)
其中,a、b為待定常數(shù);φ(S)為S的函數(shù)。
與內(nèi)場(chǎng)(實(shí)驗(yàn)室)應(yīng)力環(huán)境相比,產(chǎn)品在外場(chǎng)工作環(huán)境所發(fā)生的失效機(jī)理更加復(fù)雜,通常是由多方面應(yīng)力共同作用而出現(xiàn)的。因此,對(duì)于外場(chǎng)工作環(huán)境中產(chǎn)品分布模型的參數(shù)進(jìn)行一定的修正,則有:
μf=μAk1
(3)
其中,k1、k2為修正系數(shù);μf、μA為外場(chǎng)和試驗(yàn)應(yīng)力環(huán)境下的漂移系數(shù);σf、σA為外場(chǎng)和試驗(yàn)應(yīng)力環(huán)境下的擴(kuò)散系數(shù)。
經(jīng)推導(dǎo)可知,首次達(dá)到l時(shí)的壽命服從逆高斯分布,則外場(chǎng)應(yīng)力環(huán)境下產(chǎn)品的可靠度、概率密度函數(shù)分別為:
(1) 步進(jìn)加速退化數(shù)據(jù)?,F(xiàn)有m個(gè)試樣在n個(gè)加速應(yīng)力下進(jìn)行步進(jìn)ADT。由于每個(gè)新應(yīng)力下的退化量初值等于上一步應(yīng)力下的退化量末值,則步進(jìn)ADT中各應(yīng)力下的性能退化量可表示為:
X(t)=
(2) 外場(chǎng)退化數(shù)據(jù)。對(duì)于外場(chǎng)應(yīng)力環(huán)境下的某待估產(chǎn)品,通過(guò)測(cè)量得到該產(chǎn)品的性能退化數(shù)據(jù)Xf,其表達(dá)式如下:
Xf=[Xf(t1),Xf(t2),…,Xf(tK′)],
采用貝葉斯方法,以外場(chǎng)實(shí)測(cè)的性能退化數(shù)據(jù)Xf為評(píng)估對(duì)象,建立似然估計(jì)函數(shù);將ADT退化信息作為先驗(yàn)信息,推斷出參數(shù)集Θ的先驗(yàn)分布類型及其超參數(shù)估計(jì),表達(dá)式如下:
π(Θ|Xf)∝L(Xf|Θ)f(Θ)
(10)
其中,π(Θ|Xf)為待估參數(shù)集Θ的后驗(yàn)分布函數(shù);L(Xf|Θ)為似然估計(jì)函數(shù);f(Θ)為先驗(yàn)分布函數(shù)。
(1) 確定先驗(yàn)分布的類型。由文獻(xiàn)[11]推導(dǎo)出的加速應(yīng)力與Wiener過(guò)程參數(shù)之間的關(guān)系可知,加速因子Apq與漂移系數(shù)、擴(kuò)散系數(shù)之間滿足如下關(guān)系式:
將(2)式代入(11)式可知:
Apq=exp{b[φ(Sp)-φ(Sq)]}
(12)
(2) 確定先驗(yàn)分布的超參數(shù)估計(jì)。確定參數(shù)的先驗(yàn)分布類型后,可得各參數(shù)的先驗(yàn)分布密度函數(shù),建立似然估計(jì)函數(shù),估計(jì)出先驗(yàn)分布密度函數(shù)中的未知超參數(shù)。
確定所有待估參數(shù)的先驗(yàn)分布后,由(10)式可知,外場(chǎng)退化數(shù)據(jù)Xf下參數(shù)集Θ的后驗(yàn)分布形式為:
表1 ADT應(yīng)力下分布參數(shù)的估計(jì)值 10-4
表2 ADT應(yīng)力下分布參數(shù)估計(jì)值折算到正常應(yīng)力下的估計(jì)值 10-4
表3 外場(chǎng)修正系數(shù)k1、k2的數(shù)據(jù)樣本
表4 備選分布下參數(shù)的檢驗(yàn)統(tǒng)計(jì)量
k2~Γ(ak2,bk2)=Γ(49.15,24.95)。
(3) 參數(shù)后驗(yàn)分布確定。根據(jù)(14)式和求得的未知參數(shù)分布模型,采用WinBUGS軟件對(duì)上式進(jìn)行基于Gibbs抽樣的MCMC仿真計(jì)算,得到馬爾科夫鏈?zhǔn)諗繒r(shí)的參數(shù)集Θ后驗(yàn)均值。同時(shí),根據(jù)待估產(chǎn)品外場(chǎng)退化數(shù)據(jù)情況,估計(jì)出不同采樣次數(shù)t下參數(shù)后驗(yàn)均值并與文獻(xiàn)[2]中的參數(shù)真值進(jìn)行比較,結(jié)果見(jiàn)表5所列。
表5 不同采樣次數(shù)下參數(shù)后驗(yàn)均值
由表5可知,隨著外場(chǎng)采樣退化數(shù)據(jù)的不斷豐富,融合先驗(yàn)ADT與外場(chǎng)退化數(shù)據(jù)的參數(shù)后驗(yàn)均值越接近于真實(shí)值,估計(jì)精度越高。
(4) 評(píng)估結(jié)果分析。將本文采用貝葉斯方法來(lái)融合ADT與外場(chǎng)退化數(shù)據(jù)的評(píng)估方法記為M1,將只利用外場(chǎng)退化數(shù)據(jù)的評(píng)估方法記為M2。將兩者得到的參數(shù)Θ后驗(yàn)均值代入(9)式、(10)式,得到產(chǎn)品在外場(chǎng)工作環(huán)境下的分布參數(shù)估計(jì)值,見(jiàn)表6所列。
表6 外場(chǎng)環(huán)境下不同方法的分布參數(shù)估計(jì)
將表6中的外場(chǎng)應(yīng)力環(huán)境下的分布參數(shù)估計(jì)值代入(5)式、(6)式,得到M1的可靠度,并與M2值、真值進(jìn)行對(duì)比,結(jié)果如圖1所示。
圖1 不同方法的可靠度曲線
由圖1可以看出:① M1與真值之間的擬合程度要明顯優(yōu)于M2與真值之間的擬合程度;② 隨著外場(chǎng)采樣次數(shù)的不斷增多,M1擬合曲線逐漸向真值曲線接近。這是由于M1采用貝葉斯方法,對(duì)ADT與外場(chǎng)信息進(jìn)行了融合處理,將ADT信息作為參數(shù)估計(jì)的先驗(yàn)信息,使得參數(shù)的后驗(yàn)估計(jì)值與真值更加接近,且隨著外場(chǎng)信息的不斷充實(shí),參數(shù)后驗(yàn)估計(jì)值不斷接近真值,產(chǎn)品外場(chǎng)可靠性評(píng)估結(jié)果也不斷接近外場(chǎng)實(shí)際水平。
(1) 本文考慮到外場(chǎng)與內(nèi)場(chǎng)之間的應(yīng)力環(huán)境差異,提出了基于修正系數(shù)的Wiener過(guò)程雙參數(shù)修正模型,從總體上對(duì)同類產(chǎn)品ADT下分布參數(shù)進(jìn)行修正。
(2) 考慮到同類產(chǎn)品個(gè)體之間的退化差異,將Wiener過(guò)程參數(shù)進(jìn)行隨機(jī)化處理,利用加速因子,將各加速應(yīng)力下分布參數(shù)估計(jì)值折算到正常工作應(yīng)力下,構(gòu)成未知參數(shù)先驗(yàn)分布的數(shù)據(jù)樣本;采用CvM檢驗(yàn)法,從備選分布中確定未知參數(shù)最優(yōu)的先驗(yàn)分布類型。
(3) 采用MCMC法,仿真得到未知參數(shù)集的后驗(yàn)分布均值,并通過(guò)不斷豐富外場(chǎng)采樣數(shù)據(jù)來(lái)及時(shí)更新后驗(yàn)分布均值,實(shí)現(xiàn)產(chǎn)品外場(chǎng)可靠性的實(shí)時(shí)評(píng)估。
[1] 姜同敏.可靠性試驗(yàn)技術(shù)[M].北京:北京航空航天大學(xué)出版社,2012.
[2] 王立志,姜同敏,李曉陽(yáng),等.融合加速試驗(yàn)及外場(chǎng)使用信息的壽命評(píng)估方法[J].北京航空航天大學(xué)學(xué)報(bào),2013,39(7):947-951.
[3] WHITMORE G A.Normal-gamma mixtures of inverse Gaussian distribution [J].Scandinavian Journal of Statistics,1986,13(1):211-220.
[4] WANG X.Wiener processes with random effects for degradation data [J].Journal of Multivariate Analysis,2010,101(2):340-351.
[5] 王小林,郭波,程志君.融合多源信息的維納過(guò)程性能退化產(chǎn)品的可靠性評(píng)估[J].電子學(xué)報(bào),2012,40(5):977-982.
[6] PENG C Y,TSENG S T.Mis-specification analysis of linear degradation models [J].IEEE Transactions on Reliability,2009,58(3):444-455.
[7] 唐圣金,郭曉松,周召發(fā),等.步進(jìn)應(yīng)力加速退化試驗(yàn)的建模與剩余壽命估計(jì)[J].機(jī)械工程學(xué)報(bào),2014,50(16):33-40.
[8] TANG S J,YU C Q,WANG X,et al.Remaining useful life prediction of lithium-ion batteries based on the Wiener process with measurement error [J].Energies,2014,7(2):520-547.
[9] 彭寶華,周經(jīng)論,孫權(quán),等.基于退化與壽命數(shù)據(jù)融合的產(chǎn)品剩余壽命預(yù)測(cè)[J].系統(tǒng)工程與電子技術(shù),2011,33(5):1073-1078.
[10] 徐廷學(xué),王浩偉,張?chǎng)?EM算法在Wiener 過(guò)程隨機(jī)參數(shù)的超參數(shù)估計(jì)中的應(yīng)用[J].系統(tǒng)工程與電子技術(shù),2015,37(3):707-712.
[11] 王浩偉,徐廷學(xué),趙建忠.融合加速退化和現(xiàn)場(chǎng)實(shí)測(cè)退化數(shù)據(jù)的剩余壽命預(yù)測(cè)方法[J].航空學(xué)報(bào),2014,35(12):3350-3357.
Reliabilityassessmentmethodwithintegrationofinfieldandoutfielddegradationdata
WU Xiaohui1, CAI Zhongyi2, LI Quanxiang1
(1.Chinese Flight Test Establishment, Xi’an 710089, China; 2.College of Equipment Management and Safety Engineering, Air Force Engineering University, Xi’an 710051, China)
Aiming at the degradation process of products which obeys Wiener process, reliability assessment method with accelerated degradation test(ADT) data and field degradation data is put forward by using Bayesian statistical inference. Given the difference between field and laboratory stress environment, Wiener process double-parameters modified model with the correction factors is constructed. Degradation model is built to obtain estimated values of distribution parameters under each stress. Estimated values under accelerated stress are converted into those under regular stress, which constitute the data sample of prior distribution of unknown parameters. Posterior distribution function of unknown parameters under field degradation data is built and Markov chain Monte Carlo(MCMC) method is used to obtain mean values of posterior distribution. The accuracy and practical applicability of the presented method is verified by an example.
Wiener process; degradation data integration; Bayesian inference; stress environment difference
2016-06-07;
2016-10-17
國(guó)家自然科學(xué)基金資助項(xiàng)目(71601183);陜西省自然科學(xué)基金資助項(xiàng)目(2014JQ2-7045)
吳曉輝(1981-),男,陜西咸陽(yáng)人,中國(guó)飛行試驗(yàn)研究院高級(jí)工程師.
10.3969/j.issn.1003-5060.2017.12.002
TB114.3
A
1003-5060(2017)12-1589-05
文獻(xiàn)[2]中的數(shù)據(jù),根據(jù)其所提出的評(píng)估方法,融合先驗(yàn)ADT與外場(chǎng)退化信息,對(duì)待估產(chǎn)品外場(chǎng)壽命與可靠性進(jìn)行評(píng)估,具體分析如下:
(責(zé)任編輯胡亞敏)