楊家鑫，唐圣金，*，李良，孫曉艷，祁帥，司小勝

1.火箭軍工程大學 導彈工程學院，西安 710025

2.火箭軍工程大學 作戰(zhàn)保障學院，西安 710025

3.火箭軍裝備部駐鄭州軍代表室， 鄭州 450000

在工程應用中，設備的性能會隨著時間累積呈現(xiàn)一種緩變的退化趨勢［1］，當其退化水平超過失效閾值時，可能會導致設備停機，甚至引發(fā)事故。剩余壽命是指設備在當前時刻距離失效時刻的時間長度［2］，準確的剩余壽命預測有助于操作人員在設備接近失效前做出正確的維護決策，從而有效提高設備的可靠性安全性，并降低設備的經(jīng)濟成本［3-5］。近年來，隨著航空航天、軍事等高可靠性領域的快速發(fā)展，統(tǒng)計數(shù)據(jù)驅動的剩余壽命預測方法已經(jīng)成為可靠性領域研究的熱門問題［6-7］。

統(tǒng)計數(shù)據(jù)驅動的剩余壽命預測可以充分利用同類設備的狀態(tài)監(jiān)測數(shù)據(jù)［8］。其基本原理：根據(jù)歷史退化數(shù)據(jù)或失效壽命數(shù)據(jù)估計表示退化過程共性特征的未知參數(shù)；然后根據(jù)實時檢測的現(xiàn)場退化數(shù)據(jù)在線更新漂移參數(shù)；最后利用剩余壽命的概率密度函數(shù)（Probability Density Function，PDF）預測剩余壽命。由于設備的退化過程具有隨機性，因此剩余壽命預測經(jīng)常使用隨機過程建模［9］，基于隨機過程的剩余壽命預測方法具有物理解釋清晰、數(shù)學性質(zhì)好等優(yōu)點［10］。Gebraeel 等［11］首次提出通過貝葉斯方法將現(xiàn)場退化信息與先驗退化信息合理融合，并利用得到的后驗信息進行剩余壽命預測。該方法開始主要應用于隨機系數(shù)回歸模型［11-13］，之后，被推廣到了伽馬過程［14］、逆高斯過程［15］、基本維納過程［16］等隨機過程退化模型。值得注意的是，因維納退化過程模型同時適用于單調(diào)和非單調(diào)的退化系統(tǒng)，近些年已被廣泛使用［17-18］。然而，在工程應用中，具有非線性特征的退化系統(tǒng)廣泛存在，基本線性退化模型難以有效跟蹤非線性退化過程［19-20］，可能會導致剩余壽命預測精度偏低。此外，在工程應用中，經(jīng)常會遇到不完美或缺失的先驗退化信息，這可能會導致先驗參數(shù)估計不準確或無法估計出先驗參數(shù)［13］，進而導致剩余壽命預測精度降低甚至失?。?，21］。針對先驗退化信息不完美或缺失的問題，現(xiàn)有的文獻主要有2 種解決方法。

第1 種方法是結合期望最大化算法（Expectation Maximization，EM）和貝葉斯方法。該方法最早由Wang 等［22］提出，用于擬合帶有漂移參數(shù)的自適應布朗運動模型。然后，Si 等［21］基于線性維納退化過程模型采用貝葉斯更新方法和EM算法對模型參數(shù)和剩余壽命分布進行更新，并推導了一個精確封閉的剩余壽命分布函數(shù)。然后，該機制被擴展到了隱含線性維納退化過程［23］、非線性維納退化過程［24-25］和隱含非線性維納退化過程［26］。聯(lián)合EM 算法和貝葉斯方法可以克服先驗信息不完美和缺失對剩余壽命預測結果的影響。Tang 等［27］通過假設模型漂移參數(shù)固定，推導了隱含線性維納退化過程模型的參數(shù)解析表達式，證明了該聯(lián)合算法與傳統(tǒng)極大似然估計（Maximum Likelihood Estimation， MLE）方 法具有相同估計結果的結論。然而，為什么該方法能夠獲得比傳統(tǒng)的貝葉斯方法更好的預測結果的機理仍有待進一步深入研究。

第2 種方法是融合多源信息方法。通常，用于退化模型參數(shù)估計的數(shù)據(jù)有3 種：同類設備的歷史退化數(shù)據(jù)（包括加速退化數(shù)據(jù)［28］）、同類設備的失效壽命數(shù)據(jù)和被評估設備的現(xiàn)場退化數(shù)據(jù)，如何利用這3 種數(shù)據(jù)所包含的多源退化信息進行剩余壽命預測稱為融合多源信息的剩余壽命預測方法。Tang 等［29］基于線性維納過程模型得到了參數(shù)估計結果與退化數(shù)據(jù)特征之間的關系，為合理融合多源信息提供了理論依據(jù)。然而，在實際應用中，特別是對于一些新設備，可能缺乏足夠的失效壽命數(shù)據(jù)或歷史退化數(shù)據(jù)，因此，有必要合理利用這2 種數(shù)據(jù)來計算模型中的未知參數(shù)。目前，已有學者采用融合多源信息的方法進行剩余壽命預測［30-32］。例如，由于加速退化測試數(shù)據(jù)與現(xiàn)場退化數(shù)據(jù)存在不同應力，Wang 等［30］提出了一種基于聯(lián)合似然函數(shù)的貝葉斯評估方法，將加速退化測試數(shù)據(jù)與失效壽命數(shù)據(jù)相結合。Pang 等［31］進一步研究了基于非線性維納退化過程的融合加速退化測試數(shù)據(jù)與失效壽命數(shù)據(jù)的方法。Zhang 等［32］提出了一種新的貝葉斯框架，用于融合二元退化數(shù)據(jù)和失效壽命數(shù)據(jù)進行可靠性分析。然而，這些文獻都沒有考慮隨機效應的影響。最近，Tang 等［33］證明，考慮隨機效應的剩余壽命預測在理論上會得到更好的剩余壽命預測結果，但文獻［33］在估計未知參數(shù)時只利用同類設備的歷史退化數(shù)據(jù)，沒有融合失效壽命數(shù)據(jù)。Tang 等［29］針對基本線性維納過程和非線性維納過程提出了一種融合失效壽命數(shù)據(jù)并考慮隨機效應的剩余壽命預測方法。之后，王鳳飛等［34］進一步提出了一種融合多源信息并考慮隨機效應的剩余壽命預測方法，并通過與僅使用退化數(shù)據(jù)或失效壽命數(shù)據(jù)的預測方法進行對比分析，驗證了其方法的優(yōu)越性。然而，由于測量儀器和環(huán)境等因素的影響，測量的退化數(shù)據(jù)可能存在測量誤差，不能直接反映真實的退化狀態(tài)［35］。Si 等［36］考慮了隨機效應、退化過程不確定性和測量誤差3 層不確定性對剩余壽命預測的影響，提出了基于Kalman 濾波算法的隨機參數(shù)在線更新方法，并得到了剩余壽命分布的解析表達式。對于存在測量誤差的退化設備，Tang 等［37］證明了如果不考慮測量誤差的影響可能會導致過早維護，從而增加維護成本。因此，在剩余壽命預測過程中需要考慮測量誤差的影響。最近，Yang等［38］基于隱含線性維納過程提出了一種考慮隨機效應的融合多源信息剩余壽命預測方法，并通過實驗驗證了其方法的有效性。然而，對于隱含非線性維納退化過程，文獻中尚未查閱有學者做過相關工作。

綜上所述，在先驗退化信息不完美和缺失情況下的剩余壽命預測方法還沒有得到充分研究。因此，本文主要基于隱含非線性維納退化過程模型進行融合多源信息的剩余壽命預測研究。首先，得出了隱含非線性維納過程模型參數(shù)估計的性質(zhì)，這個工作在文獻中尚未有學者報道過。其次，利用參數(shù)估計的性質(zhì)，基于隱含非線性維納過程模型，提出了融合多源信息的剩余壽命預測方法。當先驗退化信息不完美時，該方法使用現(xiàn)場退化數(shù)據(jù)估計模型固定參數(shù)，然后使用失效壽命數(shù)據(jù)估計漂移參數(shù)先驗信息。當先驗退化信息缺失時，該方法通過融合失效壽命數(shù)據(jù)和歷史退化數(shù)據(jù)估計模型先驗參數(shù)。然后，在上述2 種情況下，基于現(xiàn)場退化數(shù)據(jù)在線更新漂移參數(shù)和實際退化狀態(tài)。最后，通過仿真數(shù)據(jù)和2 個實例驗證了該方法的有效性。

1 隱含非線性維納退化過程模型參數(shù)估計的性質(zhì)

先驗參數(shù)描述設備的共同退化特征。準確估計模型先驗參數(shù)，有助于提高剩余壽命預測的精度。在工程應用中，經(jīng)常會遇到帶有測量誤差的非線性退化系統(tǒng)，線性退化模型可能會導致參數(shù)估計結果偏離真實值，進一步導致剩余壽命預測精度降低［19］。因此，本節(jié)基于隱含非線性維納退化過程模型分析參數(shù)估計的性質(zhì)。

令Y(t)表示設備在t 時刻測量的退化狀態(tài)，則退化過程可表示為

式中：x0為設備初始退化狀態(tài)；λ 為漂移參數(shù)表征設備退化速率；b 為非線性系數(shù)；σB為擴散系數(shù)；B(t)為標準布朗運動，表示退化過程的不確定性；ε 為測量誤差。不失一般性，令x0=0。

假設有n 個同類設備，每個設備在相同時刻t1，t2，…，tm檢測退化狀態(tài)，則第i 個設備在ti時刻的監(jiān)測值可表示為

令Δyi，j=yi(tj)-yi(tj-1)，其 中1 ≤i ≤n，1 ≤j ≤m，λi表示第i 個設備的退化速率，εi，j表示第i 個設備 在ti時 刻的測 量誤差。則Δyi，j可 表示為

通過MATLAB 函數(shù)“FMINSEARCH”最大化式（8）所示的未知參數(shù)的對數(shù)似然函數(shù)ln L，可以得出參數(shù)? 的估計值，接著將其代入式（5）～式（7），可 得和的 值，進 一 步 可 得 到參數(shù)σ2

ε。

目前已有學者給出了基本線性維納退化過程模型［29］和隱含線性維納退化過程模型［38］的無偏參數(shù)估計性質(zhì)。文獻［40］給出了隱含非線性維納退化過程模型無偏估計參數(shù)的期望，如引理1所示。引理1［33］隱含非線性維納過程無偏估計參數(shù)和的期望為

更多的證明細節(jié)參見文獻［40］。然而，從引理1 并不能得知參數(shù)估計的精度與哪種因素有關，因此，本節(jié)進一步推導無偏估計參數(shù)的方差表達式，以獲取更多信息。具體結果在定理1 中給出。

定理1 隱含非線性維納退化過程模型無偏參數(shù)估計參數(shù)μ?λ、σ2B和σ2λ的方差為

上述已經(jīng)獲得了隱含非線性維納退化過程模 型 的 未 知 參 數(shù)和的 方 差。同 樣 地，也需要進行分析。令P=φΩ+F，其中φ=，那 么，則的表達式可表示為

具體推導過程不再重復。

基于定理1 和定理2，可以得出以下結論：

1） 根據(jù)式（10）、式（12）和式（13）可以得知，隱含非線性維納退化過程模型的參數(shù)μλ和σ2λ的估計精度主要受到單元數(shù)量和檢測時間長度影響。單元數(shù)量n 越多，檢測時間tm越長，得到的參數(shù)估計結果越接近真實值。因此，當歷史退化數(shù)據(jù)缺失甚至沒有時，將失效壽命數(shù)據(jù)融合到退化模型中能夠在一定程度上提高μλ和σ2λ的估計精度。

2） 根據(jù)式（11）和式（16）可以得知，表示模型退化不確定性的和估計精度與單元數(shù)量n和檢測次數(shù)(m-1)成正比。換句話說，在估計模型未知參數(shù)時，使用的數(shù)據(jù)量越多，估計的和σ2ε越接近真實值。

值得注意的是，上述2 點結論都是假設非線性系數(shù)b 為真實值。然而，非線性系數(shù)b 的估計是否影響上述結論以及非線性系數(shù)b 的估計精度與哪種因素有關，這里很難從理論上進行分析。本文在第4 節(jié)通過仿真實驗說明非線性系數(shù)b 的參數(shù)估計性質(zhì)，從仿真結果可以看出，對于隱含非線性維納退化過程，其參數(shù)估計與σ2B和σ2ε的參數(shù)估計具有相似的性質(zhì)。

2 隱含非線性維納退化過程模型剩余壽命預測

2. 1 基于Kalman 濾波算法參數(shù)在線更新

在維納退化過程模型參數(shù)中，漂移參數(shù)的隨機性描述同類設備之間的個體差異［41］。為了適應被評估設備的個體退化特征，需要使用現(xiàn)場退化數(shù)據(jù)對漂移參數(shù)在線更新。由于檢測的退化數(shù)據(jù)存在測量誤差，不能直接反映真實的退化狀態(tài)。因此，可以使用Kalman 濾波算法對漂移參數(shù)在線更新，并實時估計設備當前退化狀態(tài)［20］。

對于式（1）所示的隱含非線性維納退化過程，當檢測到被評估設備在tk時刻的退化狀態(tài)Y1：k時，tk時刻的狀態(tài)空間方程可表示為［20］

式 中：ψk=σBB(tk)-σBB(tk-1) 表 示tk時 刻 的σBB(tk)與σBB(tk-1)的差值；λk為漂移參數(shù)；xk表示tk時刻的真實退化狀態(tài)；yk表示tk時刻測量的退 化狀態(tài)；εk表示tk時刻的測量誤差；{ψk}k≥1和{εk}k≥1獨立同分布。

由于測量誤差的影響，真實退化狀態(tài)xk和漂移參數(shù)λk都是隱含狀態(tài)，為便于公式推導，進一步將式（17）狀態(tài)空間方程轉化為

式中：zk表示tk時刻的狀態(tài)最優(yōu)的估計值；Ak表示tk時刻的狀態(tài)轉移矩陣；ηk表示tk時刻狀態(tài)方程中噪聲矩陣；C 表示測量方程的觀測矩陣。

根據(jù)文獻［20］，可以通過Kalman 濾波算法更新上述狀態(tài)方程。令表示在k個退化數(shù)據(jù)條件下tk時刻的真實退化狀態(tài)估計值，表示在k 個退化數(shù)據(jù)條件下tk時刻設備的漂移系數(shù)，表示tk時刻真實退化狀態(tài)的方差，表示tk時 刻 漂 移 系 數(shù) 的 方 差，表示tk時刻真實退化值與漂移系數(shù)的協(xié)方差，則隱含非線性維納過程模型的漂移參數(shù)在線更新結果和估計的當前退化狀態(tài)可表示為

2. 2 剩余壽命預測

剩余壽命是指設備在當前時刻距離失效時刻的時間長度［2］。定義設備失效閾值為w，在tk時刻被評估設備的觀測數(shù)據(jù)為Y0：k={ y0，y1，…，yk}，真實退化狀態(tài)為X0：k={ x0，x1，…，xk}。根據(jù)首達時間的概念［29，42］，一旦獲得設備在tk時刻的真實退化量X0：k，則剩余壽命可以表示為

由于測量誤差的影響，真實退化量xk為隱含狀態(tài)，不能直接用于剩余壽命預測。針對該問題，基于Kalman 濾波算法參數(shù)更新結果，可以得出考慮狀態(tài)估計的剩余壽命PDF 表達式，具體推導過程參見文獻［20］。

3 隱含非線性維納過程先驗參數(shù)估計

3. 1 融合失效壽命數(shù)據(jù)的先驗參數(shù)估計

在實際剩余壽命預測中，經(jīng)常會遇到先驗退化信息不完美的情況，傳統(tǒng)的MLE 方法可能無法估計模型先驗參數(shù)或錯誤地估計模型先驗參數(shù)。為解決這個問題，Tang 等［29］提出了一種合理融合失效壽命數(shù)據(jù)和現(xiàn)場退化信息的剩余壽命預測方法。然而，Tang 等［29］沒有考慮測量誤差的影響，這可能會增大剩余壽命預測的不確定性，因此本節(jié)在文獻［29］的基礎上，進一步提出針對隱含非線性維納過程模型的融合失效壽命數(shù)據(jù)方法，其流程如圖1 所示。

圖1 融合失效壽命數(shù)據(jù)剩余壽命預測流程Fig.1 Flow of remaining useful life prediction based on failure time data fusion

根據(jù)第2 節(jié)得出的參數(shù)性質(zhì)可知，σ2B、b 和σ2ε的估計精度與現(xiàn)場退化數(shù)據(jù)的使用量成正比，因此，在該預測方法實施過程中，每當檢測到1 個新的現(xiàn)場退化數(shù)據(jù)就需要將圖1 的流程重新執(zhí)行1 次。

定義y0：m={ y0，y1，…，ym}為現(xiàn)場退化數(shù)據(jù)，T1：n′={T1，T2，…，Tn′}為n′個 設 備 的 失 效 壽 命 數(shù)據(jù)或是首次達到失效閾值w 的失效時間數(shù)據(jù)，則融合失效壽命數(shù)據(jù)剩余壽命預測方法的先驗參數(shù)估計分為2 個步驟。

第1 步 基于現(xiàn)場退化數(shù)據(jù)，使用MLE 獲取先驗固定參數(shù)、b 和σ2ε。

對于單個設備的退化過程建模，由于不存在個體之間的差異，因此在建立隱含非線性維納退化過程模型時，不考慮隨機效應的影響，視λ 為 固定系數(shù)。基 于式（1），令Δyj=yj-yj-1，，Δtj=tj-tj-1，1 ≤j ≤m，Δy=[Δy1，Δy2，…，Δym]T，Δv=[ Δv1，Δv2，…，Δvm]T，Δt=[ Δt1，Δt2，…，Δtm]T。則 有σε2F)，其中Ω=diag{ Δt1，Δt2，…，Δtm}，F(xiàn) 的元素見式（4），f1，1=1。令D=Ω+?F，其中，則單個設備的對數(shù)似然函數(shù)可表示為

將式（26）代入式（24），則可以得到b 和? 的剖面似然函數(shù)表達式為

通過MATLAB 的“FMINSEARCH”函 數(shù)最大化式（27）可以獲得b 和? 的估計值。將其代入式（25）和式（26）可以得到λ 和，進一步得到估計值。

第2 步 使 用 失 效 壽 命 數(shù) 據(jù) 求 取μλ和σ2λ的先驗分布。

根據(jù)Peng 和Tseng［42］的結論，隱含維納退化過程壽命分布不受測量誤差影響。則在給定和 失 效 壽 命 數(shù) 據(jù)T1：n′={T1，T2，…，Tn′}時，根 據(jù)非線性維納過程的性質(zhì)，第v 設備的故障時間Tv服 從 逆 高 斯 分 布［29］，其中1 ≤v ≤n′。Tv的 似 然函數(shù)L 為

根據(jù)式（28），失效壽命時間關于λv的對數(shù)似然函數(shù)可以表示為

給定非線性模型的失效壽命數(shù)據(jù)T1：n′={T1，T2，…，Tn′}，則 可以通 過最大化 式（29）得 到λ={ λ1，λ2，…，λn′}。接 下 來，先 驗 漂 移 參 數(shù)可 以通過式（30）計算。

3. 2 融合多源信息的隱含非線性維納過程先驗參數(shù)估計

當先驗退化信息缺失且同時存在先驗退化信息和失效壽命數(shù)據(jù)時，根據(jù)第2 節(jié)獲取的參數(shù)估計性質(zhì)可知，通過融合盡可能多的先驗退化數(shù)據(jù)可以提高模型未知參數(shù)估計的精度。最近，王鳳飛等［34］針對該問題提出了融合多源信息并考慮隨機效應的剩余壽命預測方法。然而，文獻［34］并沒有考慮測量誤差的影響，本節(jié)將該方法進一步推廣到隱含非線性維納退化過程，具體流程如圖2所示。

圖2 融合多源信息剩余壽命預測流程Fig.2 Flow of remaining useful life prediction based on multi-source information fusion

第2 步 基于歷史退化數(shù)據(jù)和失效壽命數(shù)據(jù)計算先驗漂移參數(shù)μλ和σ2λ。

假設有n 個同類設備，每個設備在相同時刻t1，t2，…，tm檢測退化狀態(tài)，基于式（1），令yi，j表示第i 個 設 備 在tj時 刻 測 量 的 退 化 狀 態(tài)，Δyi，j=yi(tj)-yi(tj-1)，，Δtj=tj-tj-1，Δyi=[Δyi，1，Δyi，2，…，Δyi，m]T，Δv=[ Δv1，Δv2，…，Δvm]T，Δt=[ Δt1，Δt2，…，Δtm]T和，其中，1 ≤i ≤n，1 ≤j ≤m，λi表示第i 個設備的退化 速 率。 則Δyi服 從 均 值 為μλΔv，方 差 為的多元正態(tài)分布。此外，假設存在n′組失效壽命數(shù)據(jù)T1：n′={T1，T2，…，Tn′}，Tv為第v 個單元的失效壽命數(shù)據(jù)，其中1 ≤v ≤n′。令D=Ω+?F，其中。受文獻［34］啟發(fā)，的剖面對數(shù)似然函數(shù)可以表示為

通過MATLAB 函數(shù)“FMINSEARCH”最大化式（31），可以得到μλ和的估計值。

融合多源信息是指融合歷史退化數(shù)據(jù)信息、失效壽命數(shù)據(jù)信息和被評估設備現(xiàn)場退化數(shù)據(jù)信息。然而，當先驗退化信息不完美時，如果基于傳統(tǒng)的極大似然方法直接使用該數(shù)據(jù)估計模型的未知參數(shù)，可能會導致參數(shù)的估計結果偏離真實值較遠，從而導致剩余壽命預測結果精度偏低。對于這一特殊情況，在使用融合多源信息剩余壽命預測方法時，可以將歷史退化數(shù)據(jù)轉換成失效壽命數(shù)據(jù)，僅融合失效壽命數(shù)據(jù)信息和被評估單元現(xiàn)場退化數(shù)據(jù)信息進行剩余壽命預測。

4 實 驗

本節(jié)通過仿真實驗和2 個實例來驗證文章提出方法的有效性和優(yōu)越性。首先，使用蒙特卡羅算法生成隱含非線性維納退化過程的仿真數(shù)據(jù)以驗證第2 節(jié)獲取的參數(shù)估計性質(zhì)。然后，使用鋰電池退化數(shù)據(jù)驗證在先驗退化信息不完美情況下的優(yōu)越性，使用疲勞裂紋數(shù)據(jù)驗證在先驗退化信息缺失情況下的優(yōu)越性。

4. 1 仿真實驗

Yang 等［38］分析了隱含線性維納退化過程模型的參數(shù)性質(zhì)，得到了參數(shù)估計精度與單元數(shù)量和檢測次數(shù)以及檢測時間成正比的結論。然而，對于隱含非線性維納退化過程模型，非線性系數(shù)b 的參數(shù)估計性質(zhì)難以通過理論推導得出，且當非線性系數(shù)b 為估計值時，第2 節(jié)獲取的性質(zhì)是否會受到影響有待驗證。本節(jié)試圖使用仿真實驗獲取影響非線性系數(shù)b 估計精度的因素，以及進一步驗證第2 節(jié)的參數(shù)估計性質(zhì)。

假 設 模 型 參 數(shù) 為μλ=1，=0.04，=0.25，b=1.5，=4，令單元數(shù)量N=80，每個單元在相同的時刻檢測退化數(shù)據(jù)，檢測次數(shù)m=40，檢測時間間隔Δt=1 h。仿真數(shù)據(jù)的部分單元退化路徑如圖3 所示。為了研究退化數(shù)據(jù)特征與參數(shù)估計精度之間的關系，分別使用5、20 和80組單元退化數(shù)據(jù)估計模型參數(shù)，來驗證單元數(shù)量對參數(shù)估計精度的影響。此外，每當生成1 個新的退化數(shù)據(jù)就將參數(shù)再估計1 次，以驗證檢測時間長度和檢測次數(shù)對參數(shù)估計精度的影響。由于仿真數(shù)據(jù)隨機生成，如果僅用1 次生成數(shù)據(jù)估計參數(shù)，其估計結果有可能難以完全體現(xiàn)退化數(shù)據(jù)與參數(shù)估計精度之間的關系。因此使用相同的方法生成20 次退化數(shù)據(jù)，每一次的參數(shù)估計方式相同。估計的參數(shù)變化曲線結果如圖4～圖8 所示。

圖3 部分仿真退化路徑Fig.3 Some degradation paths

圖4 μλ 隨單元數(shù)量N 變化的參數(shù)估計結果Fig.4 Estimation of μλ with change of number N of units

從圖4 和圖5 可以看出，參數(shù)估計變化曲線在使用單元數(shù)量N 相同時，退化數(shù)據(jù)檢測時間tm越長，得到的參數(shù)估計結果越接近真實值。這是因為檢測時間tm越長，式（10）、式（12）和式（13）中Δv′D-1Δv 數(shù)值越大，所以估計得到的μλ和不確定性越小，估計結果越準確。此外，當檢測時間達到一定數(shù)值時，參數(shù)估計結果隨著檢測時間長度的增加發(fā)生的變化不明顯，而通過增加單元數(shù)量N，參數(shù)的估計精度可以得到進一步提升，這是因為此時式（10）和式（12）中的和式（13）已經(jīng)很小，μλ和σ2λ的估計精度主要受到n的影響。因此，當先驗退化信息不完美時，僅使用失效壽命數(shù)據(jù)用來估計退化模型的先驗漂移參數(shù)可以確保其估計精度，其中失效壽命數(shù)據(jù)反映單元數(shù)量N。當同時存在失效壽命數(shù)據(jù)和退化數(shù)據(jù)時，合理地將2 種退化數(shù)據(jù)融合能夠進一步提高模型先驗漂移參數(shù)估計的精度。

圖5 σ2λ 隨單元數(shù)量N 變化的參數(shù)估計結果Fig.5 Estimation of σ2λ with change of number N of units

從圖6 和圖7 可以看出，在整個仿真周期內(nèi)，當使用單元數(shù)量N 相同時，隨著檢測次數(shù)的增加，參數(shù)估計變化曲線逐漸向真實值收斂，當檢測次數(shù)m 相同時，使用的單元數(shù)量N 越多，參數(shù)估計結果越準確，這是因為式（11）和式（16）的數(shù)值與n(m-1)成反比，即和σ2ε估計精度與檢測次數(shù)和單元數(shù)量成正比。值得注意的是，從圖8 可以看出，非線性系數(shù)b 的參數(shù)估計變化曲線與圖6和圖7 具有相同特征，即b 的估計精度同樣與單元數(shù)量和檢測次數(shù)成正比。因此，在剩余壽命預測過程中，如果能夠獲取的現(xiàn)場退化數(shù)據(jù)很多，如4.2 節(jié)的鋰電池退化數(shù)據(jù)，則在先驗退化信息不完美的情況下，當現(xiàn)場退化數(shù)據(jù)檢測次數(shù)m很大時，可以不使用歷史退化數(shù)據(jù)估計模型固定參數(shù)。

圖6 σ2B 隨單元數(shù)量N 變化的參數(shù)估計結果Fig.6 Estimation of σ2B with change of number N of units

圖7 σ2 ε 隨單元數(shù)量N 變化的參數(shù)估計結果Fig.7 Estimation of σ2 ε with change of number N of units

圖8 b 隨單元數(shù)量N 變化的參數(shù)估計結果Fig.8 Estimation of b with change of number N of units

綜上分析，仿真實驗驗證了第2 節(jié)的參數(shù)估計性質(zhì)，得出了參數(shù)估計性質(zhì)不受非線性系數(shù)b影響的結論，并獲取了影響非線性系數(shù)b 估計精度的影響因素，為文章針對隱含非線性維納過程提出的融合多源信息方法進行剩余壽命預測奠定了理論基礎。

由圖4～圖8 可以看出，隱含非線性維納退化過程模型各個參數(shù)的收斂程度不同。這是因為參數(shù)估計的方差數(shù)值越大，參數(shù)估計結果的不確定性越大，因此收斂程度越差。根據(jù)第2 節(jié)推導的參數(shù)估計方差表達式可知，退化數(shù)據(jù)特征對不同的參數(shù)的影響不同，因此各個參數(shù)的收斂程度存在差異。

4. 2 融合失效壽命方法剩余壽命預測

本節(jié)使用NASA 發(fā)布的鋰電池容量退化數(shù)據(jù)，驗證文章提出的考慮測量誤差且融合多源信息的剩余壽命預測方法在先驗退化信息不完美的情況下剩余壽命預測的有效性。該數(shù)據(jù)由4 個同類型的鋰電池在室溫下以充電、放電和阻抗測量3 種工作模式反復循環(huán)的容量退化數(shù)據(jù)組成［43］。鋰電池長時間的休息可能會出現(xiàn)容量再生現(xiàn)象，又稱弛豫效應［44］。在實際使用過程中，由于不同的鋰電池休息時間不固定，可能會導致剩余壽命無法準確預測，因此需要將原始數(shù)據(jù)中的弛豫效應消除，消除弛豫效應后的鋰電池容量退化路徑如圖9 所示。如何消除弛豫效應的更多細節(jié)參見文獻［27，44］。通常，當電池容量退化到額定容量的70%～80% 時，電池被認為失效［27］。本節(jié)實驗中假設鋰電池失效閾值為額定容量的70%（1.4 A·h），并選擇B0007 電池為被評估單元，其他電池用于模型先驗參數(shù)估計。

圖9 鋰電池容量退化路徑Fig.9 Degradation paths of lithium battery capacity

從 圖9 可 以 看 出，B0005 和B0007 電 池 的 容量退化曲線具有非線性特征，而B0006 和B0018電池的容量退化曲線具有線性特征，如果直接使用歷史退化數(shù)據(jù)估計模型參數(shù)，會弱化非線性特征，從而降低剩余壽命預測精度［29］，這種情況被視為先驗退化信息不完美。實際上，這4 組退化數(shù)據(jù)不僅包含了退化數(shù)據(jù)信息，同時也包含了壽命數(shù)據(jù)信息，可以利用其壽命數(shù)據(jù)信息預測剩余壽命［29］，此時使用B0005、B0006、B0018 電池的壽命數(shù)據(jù)預測B0007 電池的剩余壽命就是特殊情況下的融合多源信息方法的剩余壽命預測問題。為簡單起見，將文章提出的考慮測量誤差且融合失效壽命數(shù)據(jù)的剩余壽命預測方法稱為M0，文獻［29］中提出的非線性融合失效壽命數(shù)據(jù)方法稱為M1，文獻［20］中基于隱含非線性維納退化過程模型，僅使用歷史退化數(shù)據(jù)估計模型未知參數(shù)的常用剩余壽命預測方法稱為M2。根據(jù)消除弛豫效應的數(shù)據(jù)，可以計算出B0005、B0006 和B00018 的失效時間分別為85.55、69.35 和51.02 次循環(huán)，這些時間被用來計算M0和M1的先驗漂移參數(shù)。

圖10為被評估單元部分時刻的剩余壽命分布，圖11 為被評估單元在第100 次循環(huán)的剩余壽命分 布。從 圖10 和 圖11 可以看出，M0、M1和M2的剩余壽命分布都能夠覆蓋真實的剩余壽命值，且M0的圖形比M1和M2的圖形更集中于真實的剩余壽命值，這說明M0具有更高的預測精度。M0和M1的圖形對比說明，考慮測量誤差能夠提高剩余壽命預測精度。此外，通過對比M1和M2可知，雖然M2考慮了測量誤差的影響，但其剩余壽命的預測精度仍然比M1差，說明先驗退化信息不完美可能會導致剩余壽命預測結果出現(xiàn)更大的偏差。

圖10 基于鋰電池數(shù)據(jù)不同循環(huán)次數(shù)估計的剩余壽命分布Fig.10 Estimated remaining useful life at different cycle times based on lithium battery data

接下來從參數(shù)實時變化曲線來分析3 種預測方法的剩余壽命預測效果，如圖12 所示。圖12中，M2估計的始終比M0大，這是因為B0006和B0018 電池退化過程具有線性退化特征，導致M2估計的非線性系數(shù)b 偏小如圖12（c）所示，這進一步導致M2無法有效地跟蹤被評估電池的退化特征，于是被評估電池的部分非線性退化特征體現(xiàn)到了上。除此之外，M1估計的始終大于M0，這是因為M1沒有考慮測量誤差的影響，B0007 電池退化的不確定性都由體現(xiàn)。不同的是，對于M0，B0007 電池退化的不確定性由和σε2分擔。通常，越小，得到的分布圖形越尖銳，剩余壽命預測結果的不確定性越小。因此，這進一步體現(xiàn)了M0方法的優(yōu)越性。

圖12 M0、M1 和M2 在不同時刻的后驗參數(shù)Fig.12 Posterior parameters of M0， M1 and M2 at different cycle times

為了能夠更直觀地區(qū)分3 種預測方法的預測效果，計算各模型剩余壽命的均方誤差（Mean Squared Errors，MSEs）。第k 次循環(huán)的MSEs 計算表達式為

圖13 展示了3 種預測方法剩余壽命的MSEs，圖中M0的MSEs 始終小于M1和M2，這說明M0的預測精度高于M1和M2?？偟膩碚f，在先驗退化信息不完美的情況下，融合失效壽命數(shù)據(jù)方法能夠提高模型的參數(shù)估計精度，進一步提高剩余壽命預測的精度。此外，對于帶測量誤差的退化設備，預測剩余壽命時考慮測量誤差影響是必要的。

圖13 M0、M1 和M2 方法估計的剩余壽命均方誤差Fig.13 MSEs of remaining useful life of M0， M1 and M2

4. 3 融合多源信息方法剩余壽命預測

本節(jié)使用鋁合金疲勞裂紋退化數(shù)據(jù)驗證本文提出的考慮測量誤差且融合多源信息的剩余壽命預測方法在先驗退化信息缺失的情況下的有效 性。原 始 數(shù) 據(jù) 取 自Lu 和Meeker［45］，該 實 驗測試了21 個產(chǎn)品的裂紋長度，測試頻率為1 萬次循環(huán)。如果產(chǎn)品的裂紋長度超過40.64 mm，則被視為失效。在實驗驗證中使用終止時間為9 萬次循環(huán)的21 組退化數(shù)據(jù)，退化路徑如圖14 所示，詳細的數(shù)據(jù)介紹及使用方法參見Wang 等［46］。為了便于驗證融合多源信息方法的有效性，假設21 組退化數(shù)據(jù)中1 組退化數(shù)據(jù)為被評估單元的現(xiàn)場退化數(shù)據(jù)，6 組退化數(shù)據(jù)為歷史監(jiān)測數(shù)據(jù)，14 組退化數(shù)據(jù)僅能獲取其失效壽命［34］。為簡單起見，定義本文提出的考慮測量誤差的融合多源信息剩余壽命預測方法為M3。

圖14 疲勞裂紋退化路徑Fig.14 Degradation paths of fatigue crack

本節(jié)實驗以圖14 中的紅色實線為被評估單元退化數(shù)據(jù)，藍色虛線為歷史退化數(shù)據(jù)，用于估計M2的模型未知參數(shù)，黑色點劃線為失效壽命數(shù)據(jù)，根據(jù)失效閾值w=17.78 mm［34］計算其失效時間，用作M0和M3的失效壽命數(shù)據(jù)。M0估計的先驗參數(shù)如表1 所示。M2估計的先驗參數(shù)為，。M3估計的先驗參數(shù)為μλ=9.503，，，b=1.295，。

表1 M0在不同時刻的先驗參數(shù)Table 1 Prior parameters of M0 at different times

圖15 展示了被評估單元在部分時刻的剩余壽命分布，圖16 為被評估單元在9 萬次循環(huán)的剩余壽命分布。從圖15 和16 可以看出，M0和M2的剩余壽命分布都偏離真實剩余壽命值較遠，僅M3能夠較好地覆蓋真實剩余壽命值，這說明M3的預測效果優(yōu)于M0和M2。

圖15 基于疲勞裂紋數(shù)據(jù)不同循環(huán)次數(shù)估計的剩余壽命分布Fig.15 Estimated remaining useful life at different cycle times based on fatigue crack data

圖16 在9 萬次循環(huán)M0、M2和M3估計的剩余壽命分布Fig.16 Estimated remaining useful life of M0，M2 and M3 at 0.09 million cycles

圖17為3 種預測方法的后驗漂移參數(shù)實時變化曲線。從3 種方法的模型先驗參數(shù)和后驗漂移參數(shù)的估計結果可知，M2估計的模型固定參數(shù)與M3相同。然而，M0估計的b 和均小于M2和M3的估計值。根據(jù)第2 節(jié)參數(shù)估計性質(zhì)和4.1 節(jié)仿真實驗可以得知，和b 的估計精度與單元的數(shù)量和檢測次數(shù)成正比。出現(xiàn)和b 偏小的原因是被評估設備現(xiàn)場退化數(shù)據(jù)的不確定性小，且檢測次數(shù)m 較少（即式（11）中的m 較小），從而導致和b 的估計結果出現(xiàn)較大偏差。值得注意的是，由于M0估計的b 和偏小，在融合失效壽命數(shù)據(jù)時，導致了M0的漂移參數(shù)在貝葉斯框架下受到影響，進一步導致剩余壽命預測出現(xiàn)失敗。這說明對于現(xiàn)場退化數(shù)據(jù)較少的設備，融合失效壽命數(shù)據(jù)剩余壽命預測方法可能無法準確預測剩余壽命。圖18 展示了3 種預測方法在不同時刻的剩余壽命的MSE。圖中，M3的MSE 整體上小于M0和M2，這說明M3具有更好的預測效果?？偟膩碚f，將盡可能多的先驗退化信息合理融合能夠提高模型參數(shù)估計精度，進而提高剩余壽命預測的精度。

圖17 μλ 和在不同循環(huán)次數(shù)的后驗參數(shù)Fig.17 Posterior parameters of μλ and at different cycle times

圖18 M0、M2 和M3 方法估計的剩余壽命均方誤差Fig.18 MSE of remaining useful life of M0， M2 and M3

為了驗證文章提出的考慮測量誤差且融合多源信息剩余壽命預測方法對帶有測量誤差的退化系統(tǒng)能夠提高剩余壽命預測的精度，定義文獻［34］提出的非線性融合多源信息的剩余壽命預測方法為M4。圖19 為被評估單元在部分時刻的剩余壽命分布，圖20 為被評估單元在9 萬次循環(huán)的剩余壽命分布。從圖19 和圖20 可以看出，M3和M4的剩余壽命分布都能很好地覆蓋真實的剩余壽命值，但由于疲勞裂紋退化數(shù)據(jù)的測量誤差較小，M3和M4的分布圖形十分接近。

圖19 不同循環(huán)次數(shù)M3和M4估計的剩余壽命分布Fig.19 Estimated remaining useful life of M3 and M4 at different cycle times

圖20 在9 萬次循環(huán)M3和M4估計的剩余壽命分布Fig.20 Estimated remaining useful life of M3 and M4 at 0.09 million cycles

為了顯著比較考慮測量誤差的效果，這里在原始疲勞裂紋退化數(shù)據(jù)的基礎上適當增加1 個方差為0.003 的正態(tài)分布誤差。

與圖13 的數(shù)據(jù)設 定相同，圖21 為M3和M4在部分時刻的剩余壽命分布，從圖21 可以看出，M3的剩余壽命分布圖形比M4更窄，且M3的分布圖形更集中于真實剩余壽命值，這一結果表明，在退化數(shù)據(jù)存在測量誤差的情況下，使用M3能夠得到比M4更好的預測效果。接下來，計算不同時間點的MSE，如圖22 所示?？梢钥闯鯩3的MSE 始終遠小于M4，這進一步體現(xiàn)了文章提出的剩余壽命預測方法的優(yōu)越性?？偟膩碚f，對于帶測量誤差的退化系統(tǒng)，考慮測量誤差影響是必要的。

圖21 增加測量誤差情況下M3 和M4 估計的剩余壽命分布Fig.21 Estimated remaining useful life of M3 and M4 with additional measurement errors

圖22 M3 和M4 方法估計的剩余壽命均方誤差Fig.22 MSE of remaining useful life of M3 and M4

5 結 論

準確的參數(shù)估計是保證剩余壽命預測精度的前提。本文首先分析了隱含非線性維納退化過程模型參數(shù)估計的性質(zhì)。然后，基于這些性質(zhì)，針對隱含非線性維納退化過程，提出了一種融合多源信息的剩余壽命壽命預測方法。最后，通過仿真實驗和2 個實例驗證了參數(shù)估計性質(zhì)和所提剩余壽命預測方法的有效性。綜上所述，本文的主要貢獻可以總結如下：

1） 基于隱含非線性維納退化過程的無偏參數(shù)估計結果，推導得到了模型參數(shù)的方差，分析得到了隱含非線性維納退化過程的參數(shù)估計精度與單元數(shù)量、檢測時間長度和檢測次數(shù)之間的關系，為合理融合多源信息進行剩余壽命預測提供了理論依據(jù)。

2） 針對帶有測量誤差的隨機退化系統(tǒng)，為應對先驗退化信息不完美或缺失等問題，提出了一種合理融合多源信息的剩余壽命預測方法。