石 慧，康 輝，趙**，董增壽
(太原科技大學 電子信息工程學院，山西 太原 030024)
剩余壽命(Remaining Useful Life, RUL)預測是制定最優維修決策的前提和保證。根據剩余壽命預測結果提前進行維修準備，可以有效減少故障的發生，進而提高系統的可用性和可靠性，降低或避免故障造成的重大損失，同時減少維修保障費用，具有重要的研究意義和實用價值[1-3]。
剩余壽命預測是系統故障預測和健康管理領域的核心[4]。AHMADZADEH等[5-7]對剩余壽命預測建模方面的研究進展進行了綜述，其方法主要包括物理模型的方法、基于專家知識的方法和數據驅動的方法等。對大多數系統而言，難以建立準確的物理模型表征復雜系統的退化及故障的機理，同時獲取有關系統完整的專家知識成本較高，而數據驅動方法僅依賴于觀測數據，因此它們更常用于RUL預測中[8-9]。
現有的許多數據驅動實時剩余壽命預測方法需要進行模型結構假設及參數估計，需要假設其作為判斷依據的隨機取值的樣本符合某種特定模型結構，這些模型結構的假設與實際的物理模型之間常常存在較大的差距，同時參數估計的最優化過程有可能收斂到局部最小而不能保證全局最優[10]，因而預測模型將不能保證最終漸近收斂于真實的樣本模型。基于深度學習的數據驅動預測方法有強大的特征自學習能力，不需要事先對未知的退化模型進行假設，但其內部網絡結構不能清楚地表征系統退化特征的變化[11]。核密度估計(Kernel Density Estimation, KDE)不需要對數據分布進行假設，完全從數據出發研究分布特征，是一種非參數估計方法[12]。BORTOLOTI等[13]提出一種基于k均值的監督學習聚類方法，稱為監督核密度估計k均值，將核密度估計用于數據分類。HU等[14]利用核密度估計方法來建立風速的概率密度函數模型，通過與參數估計模型的比較驗證了該方法的準確性。GUAN等[15]提出一種基于改進的核密度估計模型的檢測算法估計風電機組的可靠性。張衛貞等[10]針對單部件系統建立相應的核密度估計模型，并進行壽命預測。但其僅針對單部件的退化過程及其剩余壽命預測展開研究。
大型復雜系統往往是由相互作用的多個部件組成，部件之間的相互依賴關系是不可忽略的。多部件之間的相互依賴關系包括經濟相關性、結構相關性和隨機相關性[16]。隨機相關性影響著部件的退化狀態和剩余壽命分布，是維修相關性產生的原因[17]，因此考慮隨機相關性影響的部件剩余壽命預測可以更為準確地反映系統的退化規律和未來狀態[18]。在以往的研究中，許多描述多部件系統相關性的可靠性模型是通過研究部件失效對其他部件性能的影響來建立的[19]。但在實際系統運行過程中難以獲得故障率，部件發生故障之前的各部件的連續退化狀態與實時剩余壽命的相關性在研究多部件系統預測與健康管理問題中更值得關注。針對多部件系統部件間相關性的建模，常用的方法有多維退化模型[20-21]、退化率相關模型[22-23]和基于Copula函數的退化相關模型[24-28]。多維退化模型一般采用多維正態分布來描述系統中各部件的退化過程，退化率相關模型僅考慮部件間隨機相關性對退化率產生的影響，忽視了其對部件退化狀態造成的影響，Copula函數相較于前兩種方法，在具體的多變量模型未知的情況下，可以靈活地構造多變量分布，是一種有效而靈活的隨機變量間相關性建模方法，近年來被廣泛應用于可靠性和剩余壽命建模中。LI等[24]利用Lévy Copula函數建立了多部件系統在相同環境下考慮部件間隨機相關性的可靠性模型。蓋炳良等[25]和YANG等[26]利用Copula函數建立了的隨機相關的兩部件系統退化模型。楊志遠等[27]基于伽馬過程建立系統退化模型，使用Copula函數描述退化過程的相關性，進行系統可靠度分析和剩余壽命預測。WANG等[28]提出一種基于維納過程的產品RUL預測方法，其相關性用Frank Copula函數來表征。但對于不同的應用問題，尤其存在連續退化隨機相關多部件系統，很難建立一個剩余壽命預測模型完全包含所有復雜系統的情況。
綜上所述，本文依靠傳感器獲得的實時監測數據，考慮多部件系統連續退化過程中部件間隨機相關性影響，提出一種基于Copula理論的多部件系統非參數核密度剩余壽命預測方法。首先基于核密度估計求得部件的退化分布函數；
然后采用Copula函數進行部件隨機相關性建模，利用赤池信息準則選擇最優的Copula函數，避免了不同類型的Copula函數可能導致的預測結果不準確的問題；
最后建立連續退化過程中隨機相關性影響的多部件系統非參數核密度剩余壽命預測模型并通過齒輪箱試驗臺進行試驗，驗證了所提模型的有效性和準確性。
(1)
式中：hn為部件i選用的核函數窗寬，n為部件i的單位退化增量的樣本數，K(·)為部件i選用的核函數。
(2)
窗寬hn影響擬合效果，對概率密度函數估計的準確性具有較大影響。取值過大時，平均的影響突出，使密度函數的細節變得模糊，不能清楚地反映密度函數的特點；
若取值過小，隨機性的影響會增強，使密度函數變的不光滑，不能很好地反映密度函數的趨勢。因此，隨著樣本數據的變化可采用均方誤差最小準則來選取最優的窗寬，進而準確估計其概率密度函數。常用的準則包括積分均方誤差準則(MeanIntegratedSquaredError,MISE)和漸進積分均方誤差準則(AsymptoticMeanIntegratedSquaredError,AMISE)。AMISE可以更精確地平衡偏差項和方差項的階數，更適合應用于核密度估計中最優窗寬的計算[30]，因此采用極小化AMISE來求解窗寬的最優解。其表達式為：
(3)
hopt=argminhAMISE(h)=
(4)
(5)
(6)
相應的累積退化概率密度函數為：
(7)
(8)
(9)
(10)
2.1 退化之間的相關性分析
建立部件間退化相關性模型，首先需要分析系統各部件之間是否存在相關性。常用的相關性分析方法有Pearson系數、Kendall系數和Spearman系數。Pearson系數具有很好的捕獲線性相關性的能力，但在非線性相關噪聲方面表現較差;Kendall秩相關系數用于測量兩個有序數列之間對應程度的非參數統計量，對于非高斯噪聲，Kendall秩相關系數比Pearson相關系數更適合，并且表現出比Spearman秩相關系數更強的魯棒性和有效性[31-33]，因此選擇Kendall系數進行部件退化的相關性分析。Kendall秩相關系數τ是衡量系統中各部件之間相關性的重要指標，τ的取值范圍為-1～1，當τ為0時，表示兩個隨機變量相互獨立。
0 (11) 其中： 采用Kendall秩相關系數τ判斷其部件間存在相關性后可通過Copula函數表征其相關關系。 假設系統中M個部件存在隨機相關性，其剩余壽命聯合分布函數F(t1,t2,…,tM)為M維分布函數，可用Copula函數表示為： F(t1,t2,…,tM)=C(F1(t1),F2(t2),…,FM(tM);θ)= C(u1,u2,…,uM;θ)。 (12) 式中θ為相關參數，用來表征相關性程度。 根據SKlar定理可知，選擇合適的Copula函數C就可以通過Copula函數得到退化分布函數的聯合分布函數。將部件的剩余壽命的退化分布函數式(8)帶入式(12)中得： (13) (14) (15) 第i個部件在考慮隨機相關性下的條件分布函數為： (16) (17) (18) 表1中列舉了幾種常見的Copula函數形式。 表1 常見的Copula函數 在實際應用中，Copula函數的種類很多，如表1所示為3種常用的Copula函數形式，選擇不同的Copula函數進行研究，可能會導致不同的分析結果。赤池信息準則(Akaike Information Criterion, AIC)、貝葉斯信息準則(Bayesian Information Criterion, BIC)和交叉驗證法(Cross-Validation, CV)是統計學中模型選擇和評價的重要工具。WEI等[34]通過實驗驗證了AIC和BIC比CV準則的計算速度更快，AIC準則可用于衡量估計模型的復雜度與退化數據擬合程度[9,29]，且能更好地避免參數過擬合，因此采用AIC準則來檢驗函數擬合的優劣性，AIC值越小表明模型越適合。 通過AIC準則選擇最優的Copula函數，AIC表達式為： (19) 式中：n為樣本個數，P和Q分別表示Copula函數的似然函數值和參數個數。 對于Copula函數中的未知參數θ，本文采用極大似然函數進行估計。首先利用非參數核密度估計得到退化分布函數，然后利用極大似然估計方法估計Copula函數中的未知參數θ。 (20) 通過式(20)求得不同Copula函數的參數θ及其似然函數值，將其似然函數值代入式(19)分別計算對應的AIC。選擇AIC值最小的Copula函數代入式(13)，通過式(13)～式(18)即可求解考慮隨機相關性的部件剩余壽命概率密度函數。 如圖1所示為本次試驗所用的齒輪箱試驗臺，由主試齒輪箱和陪試齒輪箱相連接構成，其中心距為150cm，試驗采用機械杠桿加載，扭矩采用轉矩轉速傳感器測量，試驗過程中對主、陪試箱體的振動、加速度、溫度及噪聲進行監測，主試齒輪箱中為正反交錯搭接的一對齒輪，齒輪的斷齒等效為齒輪的失效，軸承的外圈磨損等效為軸承的失效。 試驗臺架上共布置13個傳感器，傳感器安裝情況如圖2所示，1#、2#、3#、4#加速度傳感器分別布置在主試箱內軸承座的徑向位，7#、8#加速度傳感器分別布置在主試箱內軸承座的軸向位，5#、6#加速度傳感器分別布置在陪試箱內軸承座的徑向位； 試驗中加載分八級載荷，八級載荷的大小分別為349.5扭矩，430.7扭矩，492.2扭矩，555.6扭矩，612.9扭矩，693.4扭矩，734扭矩，822.7扭矩，前七級載荷的運行時間均為10個小時，在第八級載荷時監測到各部件發生故障，剩余壽命預測選取從第八級加載開始后的加速度數據進行分析。由圖2可知，4#和1#加速度傳感器的位置最靠近主試箱的齒輪，8#加速度傳感器的位置最靠近主試箱的軸承，因此將4#、1#和8#加速度傳感器的振動信號分別作為齒輪1、齒輪2和軸承的退化信號。其采樣頻率為25.6kHz，每次采樣持續60s，每隔9min記錄一次采樣數據。 齒輪箱的退化軌跡整體呈上升趨勢，為了更好地反映各部件的退化趨勢，采用均方幅值方法對數據進行特征提取，可由下式計算： (21) 式中：Ni為采樣點數，yj為振動信號幅值。 通過均方幅值法對齒輪箱齒輪1、齒輪2和軸承的數據進行處理后，得到如圖4～圖6所示的退化曲線圖。可以看出，齒輪1在發生故障時的故障閾值為L1=76.375mm/s2，齒輪2在發生故障時的故障閾值為L2=79.580mm/s2，軸承在發生故障時的故障閾值為L3=85.637mm/s2。齒輪1的真實壽命為77.17h，齒輪2的真實壽命為78.83h，軸承的真實壽命為83.51h。 利用齒輪箱所提取的特征值來驗證不考慮隨機相關性的齒輪1剩余壽命預測KDE模型。由式(10)可求得齒輪1的剩余壽命概率密度函數，圖7給出了不同監測時刻齒輪1的剩余壽命預測結果。從圖中可以看出，從監測的10h～70h，剩余壽命概率密度圖的方差逐漸減小，這是由于在監測初期的數據量少，隨著數據的增多，預測結果也越準確。 如表2所示為齒輪1對KDE模型的預測結果誤差分析。 表2 齒輪1剩余壽命的預測結果與實際值的誤差分析 表2給出了齒輪1在不同監測時刻的預測值、實際值和模型誤差值。由表中的數據可知，隨著監測時間的增加，數據量增大，該模型所預測的結果與真實值的誤差逐漸減小，對齒輪1剩余壽命的預測也更加準確。 同理，用同樣的方法來分別驗證不考慮隨機相關性的齒輪2和軸承的剩余壽命預測KDE模型，可以得到不同監測時刻下齒輪2和軸承的剩余壽命預測結果，如圖8和圖9所示。 齒輪2和軸承對KDE模型的預測結果誤差分析分別如表3和表4所示。 表3 齒輪2剩余壽命的預測結果與實際值的誤差分析 表4 軸承剩余壽命的預測結果與實際值的誤差分析 齒輪2和軸承的預測結果同樣說明，隨著監測時間的增加，該模型的預測誤差也逐漸減小，說明了基于核密度估計的剩余壽命模型對齒輪2和軸承剩余壽命預測的準確性。 通過式(11)求得齒輪1和齒輪2，齒輪1和軸承，齒輪2和軸承的Kendall秩相關系數τ分別為0.413 8，0.713 5，0.327 3，說明齒輪1，齒輪2和軸承為存在隨機相關的多部件系統。 常見的3種Copula函數二維概率密度圖如圖10所示。 由圖10可知，GumbelCopula函數可描述上尾相關關系，ClaytonCopula函數可描述下尾相關關系，FrankCopula適合描述具有對稱相關結構的相關關系。由于部件的退化程度越大，對其他部件的影響也越大，部件之間的退化存在尾部相關性，且這種相關性是非對稱的，可見使用上尾相關性更適合描述不同部件間的隨機相關性。 通過極大似然估計Copula函數的參數，然后根據AIC準則選擇最優的Copula函數。由式(19)得到如表5所示的3種常用Copula函數的AIC值。 表5 常用Copula函數的AIC值 由表5可知,在相同的時刻下，AIC值最小的為Gumbel Copula函數，因此根據AIC準則可知最優的Copula函數應選擇Gumbel Copula函數，可以更好地描述部件之間的相關性，這與所分析的Gumbel Copula函數特性一致。 將部件的剩余壽命分布函數式(9)和表1中的Gumbel Copula函數表達式代入式(13)中，可求得3部件的聯合分布函數為： (22) (23) (24) (25) (26) (27) (28) 用齒輪箱提取的特征值來驗證考慮隨機相關性的剩余壽命預測模型。可分別求得齒輪1、齒輪2和軸承的剩余壽命預測結果，其在不同監測時刻的預測結果如圖11～圖13所示。 圖11～圖13顯示了在不同監測時刻下，兩種不同模型對齒輪1、齒輪2和軸承的剩余壽命預測結果。從圖中可以看出，隨著監測時間的增加，兩種模型的剩余壽命概率密度函數曲線方差逐漸減小，而且真實值都在兩種模型預測的范圍之內。同時考慮隨機相關的剩余壽命預測模型結果更接近真實值。表6～表8分別對3個部件考慮隨機相關性的剩余壽命預測模型和不考慮隨機相關性的KDE模型預測結果進行誤差分析。 表6 兩種模型下齒輪1剩余壽命預測結果誤差分析 續表6 表7 兩種模型下齒輪2剩余壽命預測結果誤差分析 表8 兩種模型下軸承剩余壽命預測結果誤差分析 對比表中數據可知，隨著監測時間的增加，兩種模型預測的準確度均在增加。而且在相同監測時刻，基于Copula理論考慮相關性剩余壽命預測結果比KDE剩余壽命預測結果有更小的絕對誤差。由此表明，考慮部件之間的隨機相關性后，預測齒輪的剩余壽命比不考慮隨機相關性所預測的結果更為準確。 為進一步驗證本文所提出的模型，使用相同的樣本數據采用Gamma分布模型分別對3個部件進行剩余壽命預測，圖14～圖16分別給出了在第40 h和第60 h的Gumbel Copula模型與Gamma分布模型的預測結果對比圖。通過比較可以得出，Gumbel Copula模型方法的預測結果比Gamma分布的預測結果更接近于剩余壽命的真實值。隨著監測時間的增加，樣本增多，兩種模型的精確性也更高，而Gumbel Copula模型的預測誤差更小，說明考慮部件隨機相關性后可以更準確地對部件的剩余壽命進行預測。 本文考慮連續退化多部件系統中部件間存在的隨機相關性，首先基于核密度估計的方法進行退化分布建模；2.2 考慮隨機相關性的多部件剩余壽命預測模型
3.1 數據采集
9#、10#聲音傳感器分別懸掛在主試箱和陪試箱正上方40cm處；
11#溫度傳感器布置在主試齒輪箱內部；
12#轉速傳感器布置在驅動電機的輸出端；
13#轉矩傳感器布置在主、陪試箱的聯接軸處。主試箱齒輪嚙合情況如圖3所示。
3.2 單部件的剩余壽命預測
3.3 考慮隨機相關性的部件剩余壽命預測
然后利用Kendall秩相關系數進行相關性分析；
接著采用Copula函數表征部件退化的相關性，建立相應的剩余壽命預測模型；
最后通過試驗驗證了基于Copula的考慮隨機相關性的多部件剩余壽命預測結果比不考慮隨機相關性的剩余壽命預測結果更接近剩余壽命的真實值，證實了該方法的有效性。未來將在多部件預測的基礎上，建立混合Copula的聯合剩余壽命預測方法，進一步提高預測結果的準確性。