基于Box-Cox 變換與隨機系數回歸模型的非線性退化設備剩余壽命預測方法

楊保奎，張建勛，李慧琴，司小勝

火箭軍工程大學 智劍實驗室，西安 710025

隨著科學技術的不斷進步，各類設備呈現出集成化、智能化、多元化、復雜化的發展態勢（例如，航空發動機、導彈、核電站等）。與此同時，這些設備的造價也愈發高昂。受內在機制演化、外部環境改變等因素影響，設備往往會不可避免地發生退化，導致其使用性能降低及健康狀態劣化。如果單個設備模塊或者關鍵部件出現退化，則有可能造成整個設備發生故障，甚至引起失效并造成不可估量的人員財產損失。實踐表明，壽命預測與健康管理（Prognostics and Health Man‐agement，PHM）技術對于降低設備的故障率和提高設備的可靠性具有重要的理論研究價值與實際工程意義［1-4］。

剩余壽命（Remaining Useful Life，RUL）預測技術作為PHM 領域的核心，旨在通過設備的退化監測數據準確預測設備的剩余使用壽命，是降低設備計劃外失效風險、及時對退化設備進行維護、減少設備維護時間和費用的關鍵技術［5-6］。作為一類典型的數據驅動的設備剩余壽命預測方法，基于隨機系數回歸模型的剩余壽命預測研究開展較早。20 世紀90 年代初，Lu 和Meeker［7］首次提出了一般隨機系數回歸模型用于設備退化數據的建模及壽命預測。在該模型基礎上，出現了很多進一步的擴展研究，例如：Elwany 和Gebraeel［8-9］對隨機系數回歸模型參數估計的方法進行了改進并應用于剩余壽命預測中；Wang［10］對隨機系數回歸模型的建模思想和應用中存在的問題進行了總結；Fang 等［11］提出一種自適應預測缺失數據的回歸模型，利用自適應回歸函數對信號特征及其相應的失效時間進行建模。目前，隨機系數回歸模型作為一類能較好體現退化過程隨機性的線性模型，國內外學者不斷研究其在剩余壽命預測領域的應用與創新。Yan等［12］利用一種迭代更新回歸模型刻畫軸承健康指標的退化趨勢，并結合基于蒙塞爾變換的關鍵點檢測自適應方法，提高軸承剩余壽命的預測性能；萬昌豪等［13］克服工程實際中剩余壽命預測先驗信息不足的問題，基于非線性隨機系數模型進行退化建模；Tang 等［14］針對一類軟故障退化系統，提出一種具有時間效應的隨機系數自回歸模型描述退化過程，并利用Bayesian 理論周期性更新隨機系數，同時優化了動態條件下的維護策略；Ahsan 和Lemma［15］利用自回歸模型預測了燃氣渦輪發動機的剩余壽命。

然而，需要注意的是，隨機系數回歸模型較為依賴歷史退化監測數據的質量。在現有的基于隨機系數回歸模型的剩余壽命預測研究中，特別是對于一些非線性設備的退化數據，主要通過時間尺度變換［16］、對數變換［17］等變換技術進行線性化處理，將非線性退化數據變換為近似線性退化數據，然后建立線性退化模型用于剩余壽命預測。例如，Zhou 等［18］構建一種考慮退化速率與退化波動率依賴關系的廣義非線性Wiener 退化模型，該模型具有時變的均值方差比，并推導出近似剩余壽命的封閉解，進一步確定估計參數初始值和時間尺度函數的方法；Yu 等［19］利用年齡依賴狀態空間模型構建了一種考慮三源不確定性的非線性漂移Wiener 模型，并用渦扇發動機退化數據和合金疲勞裂紋數據進行實驗驗證；Zhang等［20］提出了一種非線性漂移分數布朗運動的鋰電池剩余壽命預測方法，該模型將非線性函數的所有參數定義為鋰離子電池退化模型的特定隱藏狀態變量；Xi 等［21］用與時間相關的非線性函數代替漂移項，用以描述系統退化過程復雜的非線性特征；Shahkar 和Khorasani［22］提出了一種非線性多變量Bayesian 模型，以解決單次測量無法得到足夠系統狀態信息的情況。但是，時間尺度變換、對數變換等技術的適用目標較為單一、變換函數的形式較為有限，限制了其在實際中的應用。同時，對于非線性退化數據，直接采用非線性隨機模型進行建模的方法較大地依賴于非線性函數形式的選擇，而對于工程實際中的非線性退化數據選擇何種非線性函數本身是一個難題。

針對以上問題，本文提出基于Box-Cox 變換（Box-Cox Transformation，BCT）與隨機系數回歸模型的非線性退化設備剩余壽命預測方法。相比于時間尺度變換、對數變換等變換技術，BCT 的優點在于可以使變換后的回歸模型滿足線性、獨立性、齊方差性、正態性，同時BCT 的函數形式更具有一般性，可以將多種常見的變化包含在特例中。在采用BCT 和隨機系數回歸模型的基礎上，運用Bayesian 理論與蒙特卡羅-期望最大化（Monte Carlo Expected Maximization，MCEM）算法對退化模型的參數進行在線更新?？紤]到本文方法中，不僅需要在Bayesian 框架下進行退化模型參數的后驗估計，還需對模型其他未知參數進行估計，而常用的期望最大化（Ex‐pectation Maximization，EM）算法需要在其E 步求解完全數據對數似然函數的條件期望，但該條件期望顯示表達式的獲得較為困難，因此本文采用MCEM 算法［23］來近似計算E 步的條件期望。根據模型參數的在線估計，推導剩余壽命的概率密度函數及點估計值，并通過數值仿真和鋰電池實際退化數據驗證本文所提方法的有效性。

1 基于BCT 與隨機系數回歸模型的退化建模

1.1 退化建模

對于非線性退化設備的退化量X(t)，首先采用Box-Cox 變換對退化數據進行線性化處理，然后再采用線性隨機系數回歸模型對變換后的數據進行建模。Box-Cox 變換是在1964 年由Box和Cox［24-25］提出的一種應用非常廣泛的廣義冪變換方法，通過計算分析變換參數λ得到最優的BCT，用以改善原始退化數據的線性、獨立性、齊方差性以及正態性，其一般形式可表示為

式中：X(t)為原始退化量；λ為變換參數；Z(t，λ)為變換后的退化量；相應逆變換可表示為

針對變換后的退化變量Z(t，λ)，考慮到BCT的非線性數據線性化能力，采用式（3）線性隨機系數回歸模型對其時變演變過程進行建模

1.2 模型參數估計

1）離線估計

2）參數在線更新

根據式（3）的性質特點和Bayesian 公式，可知：tk時刻模型隨機參數φ、θ的后驗分布為兩變量正態分布，相應的后驗分布參數為

1.3 基于MCEM 算法的參數校準

由于隨機參數θ、φ無法直接測量，通過式（7）直接計算超參數Θ的極大似然估計值較為困難。現有研究中常用EM 算法解決此類參數估計問題，但對于EM 算法，難以求解其E 步積分的顯式表達式，且迭代優化過程較長。為解決這一問題，本文引入MCEM 算法來優化模型超參數估計過程，提高收斂速度，具體實現過程如下。

進一步，結合式（8）、式（10）可得

2 剩余壽命預測

在1.2、1.3 節求出隨機參數的后驗估計后，接下來是求解設備RUL 的分布。為此，本文假定退化過程Z(t，λ)到達給定的失效閾值時發生故障（失效閾值為w，且w>0），這樣，求解剩余壽命的分布就轉化為求解退化過程到達失效閾值的時間分布。

將退化過程中tk時刻至失效時刻之間的間隔表示為T，即為設備tk時刻的RUL，那么就有Z(λ)(tk+t)=w，則剩余壽命的條件累積分布函數為

根據剩余壽命的定義其取值為非負的，即T≥0，所以需要對增加這一約束條件，那么剩余壽命的分布函數可以進一步表示為

基于以上結果，可以通過微分求得剩余壽命的條件概率密度函數為

式中：?(·)為標準正態隨機變量的概率密度函數。

直接對式（20）求解剩余壽命的點估計較為困難，而上文已經得到Z(λ)(tk+t)=w，將式（17）中t+tk時刻預測的均值代替Z(λ)(tk+t)，即w=，則可以求得tk時刻RUL 的近似點估計為

式中：RULk為本文方法對應RUL 的點估計。

3 數值分析與實例驗證

3.1 數值分析

為驗證本文模型的有效性，首先通過數值例子進行仿真。基于數值仿真數據，對模型參數進行初始化：μφ0=0.02，=0.2，μθ0=0.05，=0.01，=0.004，λ=0.41，失效閾值給定為2.5，采樣間隔為5。為進一步體現本文方法的優越性，引入文獻［13］進行對比實驗，其模型為

式中：x0為初始狀態，且x0=0；Λ(t；θ)為非線性函數，一般傳統的非線性函數為exp(θt)?1 或tθ，本文令Λ(t；θ)=exp(θt)?1；

進行數值仿真時，令φ～N(0.02，0.000 2)，θ～N(0.005，1.0×10?6)，σ2=0.000 1，文獻［13］參數設置與本文模型保持一致，仿真出的軌跡如圖1 所示。可以看出經過BCT 后的退化軌跡明顯比原始退化仿真軌跡更趨向于線性軌跡；由于文獻［13］采用典型的非線性隨機系數模型，擬合更一般的退化數據效果明顯不如本文方法。此外，變換后的退化數據共有100 個采樣點，超過失效閾值的時刻大概為496。

圖1 數值仿真與估計的退化軌跡Fig.1 Degradation trajectory of numerical simulation and estimation

通過選取21 組退化仿真數據估計，得到模型參數估計值，具體通過表1 對比真實值可以看出，數據量越大，參數估計越準確。

表1 數據累積得到的參數估計結果Table 1 Parameter estimation results obtained by data accumulation

為驗證本文采用的BCT 對非線性退化數據線性化的效果，并進一步量化本文模型擬合退化數據的精確度，引入平均預測得分、均方誤差、決定系數、相關系數作為性能評價指標［27］。通過分析模型預測的性能，以佐證本文RUL 預測方法的有效性和性能提升程度。

1）平均預測得分（Score）計算公式為

式中：τi=Ti?；Ti為ti時刻非線性退化設備真實的RUL；為ti時刻非線性退化設備RUL 的預測值，本文用退化模型RUL 的點估計表示；n為設備退化總時間。該指標越小，預測結果越準確。

2）均方誤差（Mean Squared Error，MSE）計算公式為

式中：Z(ti，λ)為非線性退化數據經過BCT 得到的退化數據；Z′(ti，λ)為經過本文模型得到的退化數據。

MSE 描述了本文模型預測誤差平方的均值，該值越小誤差越小。

3）決定系數（也稱擬合優度，R2）計算公式為

決定系數刻畫了模型預測值與真實RUL 的密切程度，該指標取值一般在[0，1]，取值越接近1 表明擬合程度越好。

4）相關系數（也稱為關聯系數，r）是統計學中重要的評價指標之一，描述2 個變量之間的相互關系和密切程度。本文采用Pearson 相關系數刻畫BCT 前的退化數據和變換后退化數據的線性相關程度，其中變換前的相關系數為

變換后的相關系數可表示為

相關系數的絕對值介于0～1，通常情況下，該值越接近于1，表明2 個變量之間的線性相關程度越密切，越接近于0，2 個變量之間的線性相關程度越弱。

5）線性度（又稱非線性誤差，δ）計算公式為

式中：δ為線性度；xmax、xmin分別為退化數據的最大值和最小值；F為退化數據的擬合函數，可用最小二乘法擬合退化數據的回歸方程表示，即

式中：p、q常數。

線性度描述的是在規定變化范圍內的實際軌跡偏離擬合特性直線的程度，本文用來描述設備退化數據偏離擬合模型的程度，線性度越接近于0，說明偏離程度越小，線性程度就越好。

給定10、30、60、100 個采樣點的仿真退化數據，然后分別對BCT 后數據進行擬合，結果如圖2所示。結合圖1，發現每增加一定的數據量，擬合程度越高。觀察模型整體仿真退化軌跡，如圖2（c）所示，可以看出大約在15 個采樣點后，模型退化軌跡與變換后的退化軌跡重合；另外，模型仿真數據與BCT 后的數據之間的MSE 為0.001 4，表明本文模型能夠較為準確地刻畫出變換后的退化數據，體現了本文方法的有效性。

圖2 不同采樣點下得到的仿真退化軌跡Fig.2 Simulation degradation trajectories obtained at different sampling points

圖3 給出了應用本文提出的MCEM 算法更新模型超參數的過程。從圖3 可以得到，本文模型的隨機參數的均值μφ、μθ分別經過大約21、15 個采樣點后與真實的參數值保持相等，收斂速度較快，可見本文采用的MCEM 算法與Bayes‐ian 思想結合的方法估計參數效果十分明顯。

圖3 模型參數估計Fig.3 Model parameter estimation

圖4為本文模型與文獻［13］中的模型在第49個采樣點時的RUL 估計及其對應的概率密度函數，從圖中可以看出，本文方法預測的結果幾乎與真實的RUL 重合，且相較于文獻［13］隨機性更小，可見本文模型能較為準確地預測出相應時刻的RUL。

圖4 第49 個采樣點RUL 的概率密度函數Fig.4 Probability density function of RUL at the 49th sampling point

圖5 給出的是通過時間累積得到的RUL 預測結果，圖6 為本文模型預測的RUL 及相應的概率密度函數，明顯看出在預測初期，模型需要經過約17 個采樣點時間后可以得到接近真實的RUL，且后期預測效果精度較高，預測結果較為穩定。

圖5 RUL 預測Fig.5 RUL prediction

圖6 本文方法得到RUL 的概率密度函數Fig.6 Probability density function of RUL obtained by this method

圖7 為文獻［13］方法預測得到的RUL，預測結果雖能較為迅速地收斂于真實RUL，但預測精度略低于本文模型。圖8 為原始數值仿真數據未經過BCT、直接通過本文模型預測得到的RUL，與圖6 對比，可以看出雖然預測結果最終仍能收斂于真實RUL，但需要經過約50 個采樣點才能較為精準地估計出真實RUL，且隨機性較大?？梢姺蔷€性退化數據經過BCT 后得到近似線性退化數據再進行RUL 預測，效果要優于未經變換的退化數據。

圖7 文獻［13］方法預測的結果Fig.7 Results predicted by Ref.［13］

圖8 BCT 前數據估計RUL 的概率密度函數Fig.8 Probability density function of RUL estimated by data before BCT

3.2 實例驗證

鋰電池具有輸出電壓高、能量密度高、使用壽命長、環境污染低等特點，作為能源儲存介質，廣泛應用于航空航天、新能源汽車、新能源發電等領域［28-30］。鋰電池作為各種設備重要零部件之一，在循環使用時，由于其內部復雜的化學反應、材料的損耗，以及腐蝕、高溫環境影響或者使用不當，就可能導致鋰電池性能降低和健康狀態劣化，直至壽命終止。2019 年4 月19 日，位于美國亞利桑那州的公共服務公用事業公司（APS）電池系統儲能電站發生爆炸，直接原因是電芯熱失控后爆炸氣體累積，和氧氣發生化學反應后燃燒；2021 年4 月16 日，北京市豐臺區某儲能電站由于一座電池間鋰電池發生內短路故障，電池模組過熱失控引發起火爆炸，事故造成了多人遇難，財產損失高達上千萬［31］。因此，為了有效降低成本并維護設備及人員的安全，準確估計出鋰電池的剩余壽命，預測健康工作時間就凸顯出重要意義。

采用馬里蘭大學的CS2-38 鋰電池的退化數據進行實例研究來驗證本文模型的有效性［32-33］。該型號鋰電池是在標準的恒定電流/恒定電壓協議下，對所有的鋰電池進行相同充電過程，在電壓達到4.2 V 之前，電流速率保持在0.5 C，之后電壓恒定在4.2 V，直至充電電流降到0.05 A 以下。放電過程的截止電壓為2.7 V，重復進行充放電實驗后，每組數據將會被隨機編號。目前，該數據集已被廣泛應用于各類鋰電池退化與剩余壽命相關實驗驗證中。

結合2.2 節，對CS2-38 鋰電池實際退化數據應用本文提出的BCT 方法，其對應的變換參數的估計值為λ=3.35。計算得到變換前退化數據的相關系數和線性度分別為?0.890 7、22.17%，變換后的分別為?0.970 7、12.96%，并且通過圖9 可以看出，鋰電池退化數據經過BCT 后線性度明顯高于變換前；此外，基于本文方法的擬合軌跡整體上與變換后的退化軌跡重合，但是本文模型由于考慮了退化過程樣本數據的差異性，導致整體退化過程的隨機性較高于參數φ=0 時的模型，二者MSE 分別為0.045 1、0.005 0，驗證了本文所提方法的合理性。

圖9 BCT 后及估計的退化軌跡Fig.9 Degradation trajectories after BCT and estimated

圖10、圖11 分別給出了鋰電池退化數據BCT 前、后經過本文模型估計得到的RUL 及其概率密度函數。作為非線性退化數據，雖然通過本文模型也能較為準確地估計出非線性設備的RUL，但隨著退化數據的累積，其預測精度和隨機性都大于經過BCT 后的預測結果。對比圖10、圖11，明顯看出通過本文方法估計RUL 的隨機性逐漸減小，預測精確性也較高，證明了本文基于BCT 方法的優越性。

圖10 本文模型RUL 的概率密度函數Fig.10 Probability density function for RUL of this model

圖11 BCT 前數據估計RUL 的概率密度函數Fig.11 Probability density function of RUL estimated by data before BCT

表2 給出了通過仿真退化數據與鋰電池退化數據在BCT 前后進行剩余壽命預測時，Score、MSE、R2、r·，t、δ這5 個性能指標對比結果。可以發現，仿真退化數據與鋰電池退化退化數據經過BCT 后用于剩余壽命預測，對比采用文獻［13］非線性退化模型剩余壽命預測，對應的模型性能指標Score、MSE 均小于變換前的退化數據，而擬合度大于變換前的退化數據；此外，相關系數較BCT 前的退化數據分別有10.05%、8.98%的提升，線性度較變換前的退化數據分別有70.77%、41.54%的提升。綜合以上對比結果表明，本文提出的方法在非線性退化設備剩余壽命預測中具有較好的適用性，與已有方法相比在預測性能上擁有一定的優勢。

表2 本文模型預測性能對比Table 2 Comparison of model prediction performance

4 結論

針對非線性服役設備剩余壽命預測問題，提出了基于BCT 與隨機系數回歸模型的非線性設備剩余壽命預測方法，采用BCT 對退化數據進行線性化處理，在此基礎上通過隨機系數回歸模型構建退化模型，運用Bayesian 理論與MCEM 算法在線更新模型參數，推導出了剩余壽命的分布函數以及其點估計值。最后通過數值分析與鋰電池非線性退化數據進行驗證研究，結果表明：針對非線性退化監測數據，經過BCT 后的退化數據進行建模進而預測剩余壽命，模型擬合度和預測準確性均有顯著提升，從而驗證了本文方法的有效性和潛在的工程應用價值。