基于多源信息的隱含非線性維納退化過程剩余壽命預測

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

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

2.火箭軍工程大學 作戰保障學院，西安 710025

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

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

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

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

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

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

1 隱含非線性維納退化過程模型參數估計的性質

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

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

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

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

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

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

ε。

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

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

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

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

具體推導過程不再重復。

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

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

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

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

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

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

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

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

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

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

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

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

2. 2 剩余壽命預測

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

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

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

3. 1 融合失效壽命數據的先驗參數估計

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

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

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

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

第1 步 基于現場退化數據，使用MLE 獲取先驗固定參數、b 和σ2ε。

對于單個設備的退化過程建模，由于不存在個體之間的差異，因此在建立隱含非線性維納退化過程模型時，不考慮隨機效應的影響，視λ 為 固定系數。基 于式（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 的元素見式（4），f1，1=1。令D=Ω+?F，其中，則單個設備的對數似然函數可表示為

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

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

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

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

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

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

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

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

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

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

假設有n 個同類設備，每個設備在相同時刻t1，t2，…，tm檢測退化狀態，基于式（1），令yi，j表示第i 個 設 備 在tj時 刻 測 量 的 退 化 狀 態，Δ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，方 差 為的多元正態分布。此外，假設存在n′組失效壽命數據T1：n′={T1，T2，…，Tn′}，Tv為第v 個單元的失效壽命數據，其中1 ≤v ≤n′。令D=Ω+?F，其中。受文獻［34］啟發，的剖面對數似然函數可以表示為

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

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

4 實 驗

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

4. 1 仿真實驗

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖17 μλ 和在不同循環次數的后驗參數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

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

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

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

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

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

圖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 結 論

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

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

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