隨機相關性影響的多部件系統剩余壽命預測

石 慧， 康 輝， 任謙力， 曾建潮， 谷豐收

(1.太原科技大學 電子信息工程學院，太原 030024； 2.中北大學 大數據與視覺計算研究所，太原 030051； 3.哈德斯菲爾德大學 計算與工程學院， 英國 哈德斯菲爾德 HD1 3DH)

近年來，隨著信息和傳感技術的高速發展，設備健康監測變得越來越智能化。以機械設備中最重要的部件之一齒輪箱為例，其健康狀況直接影響到整個設備的安全可靠性，當齒輪箱中的部件發生故障時，會加速其他部件的退化，導致機械設備損壞，不僅會產生巨大的經濟損失，還會引發一系列安全隱患[1-2]。因此有必要針對系統的關鍵部件進行實時監測并預測其剩余壽命，以便對部件進行及時的更換或維修。

系統的剩余壽命預測是故障預測與健康管理(prognostics and health management，PHM)的核心?，F階段預測剩余壽命的方法有4種：基于知識的預測方法、基于物理模型的方法、基于數據驅動的方法以及混合預測的方法。基于知識的預測方法通常需要具備有關系統的專門知識[3]，同時在復雜機械設備中，故障的物理機理很難獲得，從而導致基于物理模型的方法適用性不強[4-5]。而基于數據驅動的方法不依賴物理模型，利用接收到的數據進行預測建模，具有較高的預測精度和較強的泛化能力[6-7]。現有的一些基于統計及概率知識的數據驅動剩余壽命預測方法需假定退化過程并進行參數估計，然而假設的模型可能與實際的模型不符，并且參數尋優過程容易陷入局部最優。核密度估計方法不附加任何假設，是一種非參數估計方法，它從數據本身出發來研究樣本的分布特征。同時與深度學習等方法相比，核密度估計方法能夠清楚地表征系統退化特征的變化。因此，在預測與健康管理領域具有較高的關注度[8]。

Hu等[9]提出非參數核密度估計方法來估計風速的概率密度函數，驗證了該方法的準確性和靈活性。Wang等[10]將核密度估計方法引入魯棒隨機配置網絡，以減少噪聲數據或離群點對最終學習模型的影響。Zhang等[11]提出基于核密度估計的非線性系統異常檢測方法來模糊正常樣品之間的差異，強化潛在異常樣品的異常。Sidibé等[12]利用核密度估計系統可靠性函數，研究了核平滑參數的可變性對成本模型的影響。雖然核密度估計方法已應用于系統的故障預測與健康管理，但大多數文獻都是針對系統中的單個部件進行可靠性建模研究。

現代的設備系統越來越復雜，系統中有許多重要的子系統和部件協同工作，部件之間退化過程中存在隨機相關性影響[13]。例如風電機組齒輪箱的齒輪和軸承退化相互影響。構建多部件系統的PHM系統就要充分考慮部件間的隨機相關性[14]。

前期研究多部件系統隨機相關性的文獻主要集中于由部件故障觸發的相互作用。一些文獻研究部件隨機相關性時僅考慮部件故障會對其他部件的退化產生影響，建立基于時間的維修決策模型[15-16]。而在實際運行過程中，獲取故障率需要大量破壞性試驗，研究具有隨機相關的多部件故障預測與健康管理問題更關注的是系統故障前各部件的連續隨機退化特征與實時剩余壽命的關聯性。Rasmekomen等[17]指出部件間的隨機相關性不僅可以由部件故障觸發，還可能由部件的退化觸發，并將此定義為退化相互作用，使用回歸模型來表征部件退化狀態對其他部件退化速率的影響。Bian等研究一個部件的退化程度如何影響系統中其他部件的退化率，并將這種由隨機相關引起的部件間相互作用稱為退化率相互作用(degradation-rate-interaction，DRI)，建立了離散DRI模型來描述系統中相互依賴的部件退化過程之間的相互作用。Do等[18]針對具有隨機相關性和經濟相關性的兩部件系統，提出了基于狀態的維修模型。Shahraki等[19]針對由隨機相關部件組成的復雜系統進行系統建模并制定選擇性維修策略。Xi等[20]將多部件系統中隱藏的退化相關性視為噪聲影響建立狀態空間模型，提出了在線剩余壽命預測方法。上述文獻大多假設部件退化模型，并進行參數估計。宋仁旺等[21]利用Copula函數對齒輪箱振動和噪聲之間的相關性進行建模，求得齒輪箱的剩余壽命。Li等[22]研究了考慮隨機退化的多部件系統的視情維修問題，采用Lévy Copula函數對系統的隨機相關性進行建模。楊志遠等[23]基于伽馬過程建立系統退化模型，使用Copula函數描述退化過程的相關性，進行系統可靠度分析和剩余壽命預測。Copula函數可以靈活構造多變量分布，但通過Copula函數構建多部件系統中部件間的聯合分布函數，不能明顯地反映部件間的隨機相關性特征，不同Copula函數的選擇也會對剩余壽命的預測結果產生較大偏差。

本文在對部件間的退化隨機相關性特征及對退化狀態的影響進行系統分析的基礎上，研究多部件系統基于非參數核密度估計的實時剩余壽命預測方法。首先針對多部件系統中的復雜雙向隨機相關性進行建模；引入窗寬因子構建基于自適應窗寬的核密度估計退化模型以解決高密度區域過度平滑或低密度區域噪聲干擾的問題；然后通過首達時間的概念，求解部件的剩余壽命概率密度函數；最后在斜齒輪箱試驗臺上進行試驗，利用調制信號雙譜分析對傳感器數據進行特征提取，在識別非線性影響的同時考慮幅值和相位的影響，從而有效地抑制隨機噪聲的干擾，驗證了所提模型的準確性和適用性。

1 核密度估計模型的構建

1.1 隨機相關性分析

從多部件系統部件間存在的復雜隨機相關性的特征以及對部件連續退化狀態的不同影響角度出發，可將部件間的隨機相關性分為三類：單一單向隨機相關性、單一多向隨機相關性和雙向隨機相關性。隨機相關性示意圖，如圖1所示。

(a) 單一單向隨機相關性

(b) 單一多向隨機相關性

(c) 雙向隨機相關性圖1 隨機相關性示意圖Fig.1 Schematic diagram of stochastic dependence

(1) 單一單向隨機相關性：是指某個部件的退化只會對系統中的單個部件的退化產生單向影響；部件a自身的退化只對部件b的退化產生單向影響，部件a不受其他任何部件退化的影響。

(2) 單一多向隨機相關性：是指某個部件的退化會對系統中多個部件的退化產生單向影響或者是某個部件的退化會受到系統中其他多個部件退化的單向影響；部件a的退化會對部件b、部件c等多個部件的退化產生影響，其自身退化不受其他任何部件退化的影響。

(3) 雙向隨機相關性：是指部件的退化會受到系統中其他部件的影響，同時自身的退化也會對其他部件的退化產生影響。部件a自身的退化會對部件b和部件c等多個部件的退化產生影響，同時，這些受影響的部件退化到一定程度也會對部件a的退化產生影響，這種影響是雙向存在的。

本文考慮更為復雜的部件間具有雙向隨機相關性的多部件系統進行研究。假設一個系統中的部件i和部件j具有雙向隨機相關性，即一個部件的退化會加劇另一個部件的退化，反之亦然。隨著現代傳感器技術的發展，可以通過傳感器監測系統中各部件的歷史退化狀態數據及實時運行數據，利用這些樣本數據進行隨機相關性建模并建立剩余壽命預測模型。

1.2 雙向隨機相關性的模型構建

(1)

(2)

(3)

(4)

(5)

(6)

將式(5)代入式(6)，可以求出

(7)

式中，σk-1為k-1個已知樣本的標準差。

(8)

(9)

(10)

1.3 自適應窗寬下核密度估計退化分布模型的實時更新

(11)

部件i在tk時刻基于自適應窗寬的核密度估計為

(12)

2 部件退化分布的計算

(13)

(14)

tk-1+m時刻獲得m個樣本數據后，概率密度函數為

(15)

同理，在任意時刻獲得新的樣本數據后，都可獲得相應的累積退化概率密度函數。

3 剩余壽命預測模型的建立

(16)

由式(15)可知

(17)

將式(17)代入式(16)，可得

(18)

則tk-1時刻部件i的剩余壽命概率密度函數為

(19)

4 試驗驗證與分析

試驗中將監測到的部件振動信號數據利用調制信號雙譜(modulation signal bispectrum，MSB)變換進行特征提取，得到可以表征部件退化過程的退化狀態曲線，同時確定部件的閾值。將特征提取后的退化數據代入文中的剩余壽命預測模型，便可通過式(16)～式(19)求得部件的剩余壽命概率密度函數，求得部件的剩余壽命。剩余壽命預測試驗驗證流程如圖2所示。

圖2 剩余壽命預測試驗流程圖Fig.2 Experimental flow chart of remaining useful life prediction

4.1 數據采集

本文采用英國哈德斯菲爾德大學效率與性能工程中心的斜齒輪箱試驗臺進行試驗，試驗臺的基本結構，如圖3所示。兩個齒輪箱(GB1和GB2)采用背靠背的結構進行安裝，它們分別與驅動電機和負載電機相連。在GB1和GB2的箱體表面分別安裝傳感器，用于采集信號。采用16通道數據采集系統(data acquisition system，DAS)采集振動信號、聲學信號和油溫信號，電機額定轉速為1 500 r/min，振動信號的采樣頻率為96 kHz。

圖3 斜齒輪箱試驗臺Fig.3 Helical gearbox test bench

本文使用的兩個斜齒輪箱的幾何參數,如表1所示。

表1 齒輪箱規格Tab.1 Gearbox specification

圖4顯示了安裝在兩個齒輪箱中的齒輪分布情況，試驗采用齒輪正反交錯搭接的嚙合方式。

圖4 齒輪箱示意圖Fig.4 Schematic diagram of gearboxes

軸頻和嚙合頻率可以通過表2所示的公式計算。

表2 軸頻和嚙合頻率Tab.2 Shaft frequency and mesh frequency

為模擬工業設備在實際運行過程中環境和負載的隨機變化情況，在試驗期間，首先控制驅動電機以正弦曲線變化的轉速運轉30 min，然后在5種不同負載條件下進行勻速運轉(全速的50%)，最后在同樣的5種負載下以不同的恒定速度(全速的70%)運行。詳細的測試條件,如圖5所示。

圖5 循環運行條件Fig.5 Cycle operating conditions

試驗持續了838 h，監測到GB2的振動出現明顯峰值后停止試驗。打開齒輪箱后，觀察到GB2在低速階段發生磨損，齒根部位磨損較為嚴重，如圖6所示。

(a)

(b)圖6 838 h的齒輪磨損情況Fig.6 Gear wear of 838 hours

4.2 數據處理

試驗運行了838 h，齒輪在運行初期樣本特征波動較大，為齒輪嚙合階段，使用這個階段的數據來分析部件的磨損退化過程是不恰當的，因此選擇運行300 h之后的數據來分析，對振動數據進行MSB變換，進行特征提取并預測其剩余壽命。

傳統的頻域特征提取方法有頻譜分析、功率譜分析、常規雙譜分析等，頻譜和功率譜分析在對信號經過傅里葉變換后僅對幅值信息進行分析，忽略了相位信息。而常規的雙譜B(fx,fc)通過傅里葉變換后在頻域內可表示為

B(fx,fc)=E[X(fc)X(fx)X*(fc+fx)]

(20)

式中：X(f)為信號x(t)的離散傅立葉變換；fx為調制頻率；fc為載波頻率；X*為X的復共軛。

對于振動信號x(t)其MSB變換在頻域中表示為

BMS(fx,fc)=

E[X(fc+fx)X(fc-fx)X*(fc)X*(fc)]

(21)

MSB對常規雙譜改進后在進行退化特征提取時，同時考慮信號幅值和相位的變化，能夠更好地說明調制信號和載波信號之間的非線性關系，充分表示信號的調制特性，要比常規的雙譜分析方法更準確[26-27]。

以幅值和相位來表示，可將式(21)寫為

BMS(fx,fc)=E〈|X(fc+fx)|X(fc-fx)||X*(fc)|·

|X*(fc)|exp[jφMS(fx,fc)]〉

(22)

式(22)中MSB的總相位可通過以下方式計算

φMS(fx,fc)=φ(fc+fx)+φ(fc-fx)-

φ(fc)-φ(fc)

(23)

φ(fc)和φ(fx)耦合時，相位關系可以表示為

φ(fc+fx)=φ(fc)+φ(fx)

(24)

φ(fc-fx)=φ(fc)-φ(fx)

(25)

將式(24)和式(25)代入式(23)，可以得到MSB的總相位為零，它的幅值由其4個分量幅值的乘積所確定，因此，如果(fc+fx)和(fc-fx)來源于fx和fc的非線性影響，在雙頻BMS(fx,fc)處會出現明顯的雙譜峰。這樣對MSB的表述將更為準確。與此相反，如果fx和fc無非線性作用，或者它們是隨機噪聲，其MSB的相位隨時間變化，最終平均結果為零，由此MSB具有很好地抑制噪聲作用。圖7給出了采用MSB對本文試驗臺數據特征提取的調制雙譜分析圖。

(a) MSB幅值譜

(b) MSB相干譜圖7 MSB雙譜分析Fig.7 MSB bispectrum analysis

如圖7所示，通過MSB雙譜分析先從幅值譜發現具有峰值的信號幅值，同時，該坐標點對應的相位譜中也有較大的相干系數，則說明該幅值可以反映部件的退化特征，若相位譜中對應著非常小的相干系數，說明該幅值不是通過調制作用生成，可能是由噪聲產生的，從而將其剔除。其在識別非線性影響的同時考慮幅值和相位的影響，從而有效地抑制隨機噪聲的干擾。提取能夠表征部件退化的振動信號特征數據出來，便可得到部件的退化狀態曲線圖。

圖4中GB2中的齒輪5和GB1中的齒輪3通過主軸進行連接，二者之間具有較強的相關性，故將齒輪5和齒輪3作為部件i和部件j進行研究。采用MSB特征提取方法從MSB頻譜中將能夠表征部件退化的特征幅值提取出來，便可得到部件i的退化狀態曲線如圖8所示，GB2在838 h監測到磨損嚴重，發生故障。其故障閾值為0.229 mm/s2。

圖8 部件i的退化狀態Fig.8 Degradation state of component i

4.3 結果分析

為了說明選擇不同的窗寬對核密度估計結果的影響，圖9給出了部件i在450 h自適應窗寬與固定窗寬下的核密度估計結果。

(a) 樣本分布情況

(b) 核密度估計結果圖9 兩種窗寬下核密度估計結果Fig.9 Results of kernel density estimation under two bandwidths

由圖9可知，基于自適應窗寬的核密度估計對樣本的擬合程度更高，與固定窗寬下的核密度估計結果相比，更能夠反映樣本的分布情況。

將MSB特征提取后的退化數據代入模型，可根據式(16)～式(19)求得部件i在當前監測時刻的剩余壽命概率密度函數。圖10給出本次試驗在7個監測時間下的剩余壽命預測結果?？梢钥闯?，隨著監測時間的增加，樣本數據增大，本文提出方法的剩余壽命預測結果越來越接近真實值。說明本文提出的基于自適應核窗寬的核密度估計剩余壽命預測方法可以很好的對部件i剩余壽命的概率密度函數進行估計。

圖10 核密度估計方法在不同監測時間下剩余壽命預測結果Fig.10 Remaining useful life prediction results of kernel density estimation method under different monitoring time

Do等考慮具有隨機相關性和經濟相關性的兩部件系統，使用參數估計進行建模計算。為了對本文所用的非參數方法進行進一步驗證，使用相同的樣本數據，基于參數估計的方法進行剩余壽命預測，與本文方法進行比較。

基于參數估計的剩余壽命預測方法結果，如圖11所示。圖11中：虛線為基于參數估計方法計算的剩余壽命預測結果；實線代表本文利用自適應核窗寬核密度估計方法得到的剩余壽命預測結果；“*”號為部件i的真實剩余壽命結果。從圖11可以看出，隨著監測時間的增加，樣本數據量增大，采用非參數核密度估計方法進行剩余壽命預測更加接近真實的剩余壽命，說明本文提出的非參數核密度估計方法與參數估計相比，有一定的優越性。

圖11 參數估計與核密度估計方法比較結果Fig.11 Comparison between parameter estimation and kernel density estimation

表3給出了參數估計和核密度估計方法對部件剩余壽命預測結果的誤差分析。從表3可以看出，隨著監測時間的推移，獲得的樣本數據越多，兩種方法對部件剩余壽命的預測就越精確。在同一監測時間下，本文提出的核密度估計的方法要比參數估計方法更精確，進一步驗證了非參數核密度估計方法的有效性和精確性。

表3 參數估計和核密度估計剩余壽命預測誤差分析Tab.3 Error analysis of residual life prediction of parameter estimation and kernel estimation

圖12用三維網格曲面圖和二維等高線圖繪制了部件i的平均剩余壽命和部件i與部件j退化的關系圖，以此來說明部件j對部件i的隨機相關性影響。圖12(a)中，在部件i自身退化量一定的情況下，部件i的平均剩余壽命隨著部件j退化量的增加而減少；在部件j退化量一定的情況下，部件i的剩余壽命隨著部件i自身退化量的增加而減少。圖12(b)中，從淺色區域到深色區域，部件i的平均剩余壽命逐漸減小，同樣反映出部件i的平均剩余壽命不僅與自身退化有關，還會受到部件j退化量的影響。說明部件j對部件i具有隨機相關性影響。

(a) 三維網格圖

(b) 二維等高線圖圖12 部件j對部件i的隨機相關性影響Fig.12 The stochastic dependence effect of component j on component i

同理，圖13用三維網格曲面圖和二維等高線圖來說明部件i對部件j的隨機相關性影響。圖13(a)中，在部件j自身退化量一定的情況下，部件j的平均剩余壽命隨著部件i退化量的增加而減少；在部件i退化量一定的情況下，部件j的剩余壽命隨著部件j自身退化量的增加而減少。圖13(b)中，從淺色區域到深色區域，部件j的平均剩余壽命逐漸減小，同樣反映出部件j的平均剩余壽命不僅與自身退化有關，還會受到部件i退化量的影響。說明部件i對部件j具有隨機相關性影響。

(a) 三維網格圖

(b) 二維等高線圖圖13 部件i對部件j的隨機相關性影響Fig.13 The stochastic dependence effect of component i on component j

綜上可知，部件i和部件j具有雙向隨機相關性影響，進一步說明了本文考慮多部件系統部件間的雙向隨機相關進行建模的必要性。

為了更好的說明所提模型的通用性，采用Xi等考慮多部件隨機相關性試驗所用的數據，再次對文中模型進行驗證分析。

部件i的退化狀態,如圖14所示。

圖14 部件i的退化狀態Fig.14 Degradation state of component i

將退化數據代入本文核密度估計模型，并與參數估計模型進行比較，得到兩種方法下部件i在不同監測時刻的剩余壽命預測結果,如圖15所示。

圖15 參數估計與核密度估計方法比較結果Fig.15 Comparison between parameter estimation and kernel density estimation

從圖15可以看出，隨著監測時間的增加，兩種方法的預測結果越來越接近真實值，而本文提出的核密度方法比參數估計的方法更加準確，再一次證明了文中提出的自適應核密度估計模型的適用性。

5 結 論

針對多部件系統部件連續退化過程中存在的雙向隨機相關性，基于自適應窗寬的核密度估計方法建立了實時剩余壽命預測模型。首先針對部件間存在的雙向隨機相關性影響，引入條件核密度估計來進行建模，然后采用自適應窗寬的核密度估計方法求得相應的概率密度函數，最后建立實時剩余壽命預測模型，得到部件的剩余壽命。利用斜齒輪箱試驗臺進行驗證分析，結果表明該方法可以很好的進行壽命預測，通過和基于參數估計的剩余壽命預測方法進行比較，進一步驗證了該方法的準確性和有效性。