考慮環(huán)境改變的部件核密度剩余壽命預測

趙 彬, 李佳娟, 石 慧, 任謙力, 康 輝

(太原科技大學 電子信息工程學院,太原 030024)

風電機組在運行過程中受風速風向等環(huán)境因素影響載荷周期波動,葉片、主軸和齒輪箱等關鍵部件易發(fā)生故障。其關鍵部件的損壞會產生經濟損失并引發(fā)一系列安全隱患,因此對風電機組的關鍵部件進行實時監(jiān)測并預測其剩余壽命,可以為更合理地制定維修策略提供依據[1-2],并且提高系統(tǒng)的可用性和可靠性。

現有的一些數據驅動的剩余壽命預測方法需對模型假設和參數估計,但是風電機組等大型復雜裝備運行周期長且破壞性試驗造價昂貴,難以獲得足夠的失效數據,假設的模型往往與真實的物理模型有一定差距,且參數估計的最優(yōu)解不能保證全局最優(yōu)。因此不適用通過統(tǒng)計方法[3-4]進行剩余壽命預測,而深度學習剩余壽命預測方法能夠描述特征數據和剩余壽命的非線性關系,但無法得出模型剩余壽命的概率分布[5-6]。

核密度估計是一種非參數估計的方法,它從數據本身出發(fā),無需假設數據分布[7-8]。Lang等[9]使用核密度估計對半導體制造工藝中的異常點進行監(jiān)測。Hu等[10]采用核函數估計樣本的局部密度,在此基礎上提出了局部核回歸密度估計器和分層策略結合多尺度鄰域的信息細化樣本的異常因子進行數據聚類及異常檢測。宋仁旺等[11]考慮多退化量影響建立核密度估計和隨機濾波結合的模型的系統(tǒng)進行剩余壽命預測,僅根據故障數據采用核密度方法估計部件初始壽命分布。

風電機組在運行中可能遭遇雷電、雨雪、側風等惡劣天氣,這些環(huán)境變化對系統(tǒng)運行會產生強干擾,對風機葉片、齒輪等部件造成沖擊,影響系統(tǒng)的退化以及剩余壽命,故在建立環(huán)境變化的剩余壽命預測的模型中,環(huán)境的改變被建模為部件受到的沖擊。Sidibé等[12]考慮運行環(huán)境變化導致部件可靠度改變,采用核密度的方法估計系統(tǒng)可靠度函數并建立基于時間的維修決策模型。Wang等[13]以液壓滑閥和微機電系統(tǒng)為研究對象,建立了考慮沖擊影響的系統(tǒng)可靠度模型。

受到環(huán)境影響部件的隨機沖擊過程與自然退化過程之間具有相關性,這種相關性不可忽視,現有很多文獻都從失效的角度考慮相關性,如Wang等[14]認為累積沖擊會影響系統(tǒng)的退化率,并且在不同退化階段影響不同;而Gao等[15]則研究了不同種類的沖擊對系統(tǒng)退化量和退化率影響,并推導可靠度模型;Dong等[16]認為沖擊會影響系統(tǒng)發(fā)生的故障率,并對系統(tǒng)其進行維修決策;以上文獻僅描述了沖擊對退化單方面的影響,然而當退化到一定程度時,系統(tǒng)抵抗沖擊的能力變弱,會使得沖擊導致的失效更容易發(fā)生,由此可見,退化和沖擊之間并不是毫無關聯,而是相互依賴的。如Fan等[17]重點研究了退化對沖擊的影響,認為自然退化影響沖擊的到達率;Sun等[18]將沖擊分為非致命沖擊和致命沖擊,認為致命沖擊發(fā)生概率會隨著時間的變化而增加,從沖擊的角度解釋了系統(tǒng)抵抗沖擊的能力隨時間變化而減弱的情況。以上的文獻均是從失效的角度考慮退化和沖擊的相關性,并將此過程建模為競爭性失效[19],這種失效過程缺乏考慮由環(huán)境變化而導致的閾值變化。常春波等[20]考慮系統(tǒng)退化后突發(fā)失效閾值會降低,建立了基于變失效閾值的競爭失效模型。Hao等[21]研究了沖擊導致的硬失效與軟失效之間的相關性,且認為硬失效閾值受退化狀態(tài)的影響且是正比關系。Feng等[22]考慮了隨機沖擊的閾值和自然退化且是一種線性關系。

然而在可實時監(jiān)測的系統(tǒng)運行過程中,用線性關系描述沖擊閾值和自然退化之間的關系并不準確,剩余壽命預測建模更需要關注失效前的隨機沖擊過程與連續(xù)退化過程之間的隨機相關性,系統(tǒng)受到沖擊影響后可能同時改變當前退化量以及之后的退化速率,進而影響其剩余壽命分布。

本文假設惡劣環(huán)境對部件的影響為部件受到的隨機沖擊,在對狀態(tài)監(jiān)測數據隨機統(tǒng)計分析基礎上,將自適應核密度估計引入實時狀態(tài)及剩余壽命預測中,提出了一種考慮環(huán)境沖擊影響的核估計剩余壽命預測模型。首先,當部件運行環(huán)境改變后,其受外界環(huán)境的隨機沖擊影響,分析了沖擊與自然退化的隨機相關性;其次,利用自適應核密度估計對部件的連續(xù)自然退化過程建模,并引入可變突發(fā)失效閾值,推導出環(huán)境改變后部件核密度估計的剩余壽命,最后采用風電機組退化數據和齒輪箱磨損實測數據進行分析驗證其有效性。

1 部件所受沖擊與自然退化過程隨機相關性分析

將部件運行的不同環(huán)境記為OEm(m=1,2,3,…),部件運行環(huán)境惡劣可能導致退化加速,其剩余壽命減少,可靠性降低,對剩余壽命預測影響較大,因此本文考慮部件環(huán)境OE1改變?yōu)榄h(huán)境OE2,即環(huán)境更加惡劣的情況。假設T1時刻部件運行環(huán)境改變,部件受到沖擊影響,部件所受沖擊與連續(xù)退化過程之間存在隨機相關性,隨機相關性影響可分為以下三種情形(如圖1所示):

圖1 環(huán)境改變的部件退化過程示意圖Fig.1 Schematic diagram of the process of component degradation as the environment changes

(1) 當tT1時,部件處于惡劣環(huán)境OE2,存在自然退化和沖擊,部件在退化過程中,當退化量小于L時受到外界幅值大于D0的沖擊則發(fā)生失效,定義為突發(fā)失效,部件突發(fā)失效閾值記為D0,如圖1(a)所示。

(2) 部件在環(huán)境改變受到沖擊后,部件的退化狀態(tài)發(fā)生改變,其中沖擊的幅值為Wn,n=0,1,2,…,∞,沖擊造成的退化增量為Yn,n=0,1,2,…,∞,經過t時間后退化量大于失效閾值H0,部件也會發(fā)生失效,定義為自然退化失效,TH為部件自然退化失效閾值的時刻,如圖1(b)所示。

(3) 設L是使部件突發(fā)失效閾值改變的退化量閾值,環(huán)境改變后,部件運行在退化量小于L時,表明部件性能較好,受到幅值大于D0的沖擊才會突發(fā)失效,當部件運行在退化量大于L小于H0階段,表明部件性能變差,抵御外界環(huán)境變化能力降低,受到幅值較小的沖擊即會發(fā)生突發(fā)失效。改變后的突發(fā)失效閾值為D1(D1

1.1 部件自然退化過程

設tk為當前監(jiān)測時間,假設在時間[0,tk]部件的退化增量ΔX1,ΔX2,…,ΔXk是獨立同分布隨機變量,f(Δx)為退化增量的概率密度函數,f(Δx)的核密度估計為

(1)

式中:k為時間[0,tk]的退化增量Xi樣本個數;h為窗寬;K為核函數。

核函數的種類很多,但是不同核函數對積分均方誤差的影響非常小,選擇常用的高斯核函數作為核密度估計的模型

(2)

采用積分均方誤差(mean integrated squared error,MISE)對核密度估計與退化增量的概率密度函數之間的誤差進行測量

(3)

(4)

代入高斯核函數K(Δxk)得到最優(yōu)窗寬hk

(5)

式中,σk為k個初始樣本特征退化增量ΔX1,ΔX2,…,ΔXk的方差,則自然退化過程X(t)在不同時刻的概率密度函數g(x)可通過退化增量的概率密度函數卷積求得

(6)

部件在t時刻累積分布函數記為G(x,t),當處于環(huán)境OE1中,此時不考慮沖擊影響,只考慮自然退化,可求得部件的可靠度函數為

(7)

1.2 環(huán)境改變沖擊對自然退化的隨機相關性影響

假設部件在環(huán)境OE2受到的沖擊遵循速率為λ的泊松過程{N(t),t>0},到t時刻時沖擊發(fā)生次數的概率為

(8)

則環(huán)境改變后t時刻由沖擊造成的退化量的總量為

(9)

若沖擊幅值都小于突發(fā)失效閾值,部件不會突發(fā)失效,第n次沖擊后部件不會突發(fā)失效的概率是

(10)

式中:n=1,2,…,∞,d=0,1;FW(w,Dd)為Wn的累積分布函數;Φ(·)為標準正態(tài)分布函數。

1.3 可靠度函數建立

部件在時間T1從OE1進入OE2后,可能會受到多次沖擊,設RS(t)為部件在時間T1環(huán)境改變后的可靠度函數;環(huán)境OE2中部件的退化不僅包含自然退化,也包括沖擊引起的退化,因此其退化過程建模為

XS(t)=X(t)+S(t)

(11)

(1) 部件運行環(huán)境剛剛改變還沒有受到沖擊,此時,考慮部件不發(fā)生自然退化失效的概率

RS1[t|N(t)=0]=P[X(t)R1(t,H0)

(12)

(2) 改變環(huán)境后到t時刻時沖擊次數不為0,若部件的退化量XS(t)

(13)

(3) 改變環(huán)境后時刻受到多次沖擊,部件退化量H0>XS(t)>L,突發(fā)失效閾值發(fā)生改變,將沖擊分為兩個階段,假設第j-1個沖擊時,部件的退化量XS(tj-1)

(14)

綜上所述,部件在t時刻的可靠度為

(15)

將式(12)、式(13)和式(14)代入式(15)可得系統(tǒng)可靠度函數表達式為

(16)

式中,tj-1與tj分別為第j-1次沖擊到達的時刻和第j次沖擊到達的時刻,均為隨機變量,計算過程中可用其均值替代。其推導過程為:

設{N(t),t≥0}是到達率為λ的泊松分布,若已知在[0,t)時間內有n個沖擊到達,則這n個沖擊到達的時間是相互獨立的隨機變量,服從[0,t)上的均勻分布。

(17)

設tj的密度函數為p(tj),分布函數為F(tj),t1,t2,…,tn,則tj的概率密度函數為

(18)

由于tj服從[0,t)的均勻分布,可得tj的概率密度函數

(19)

則tj的期望為

(20)

其中

(21)

求得tj的期望為

(22)

假設沖擊幅值的均值與造成退化量的均值間是比例關系,即μY=aμW,a為常數;某時刻的自然退化量為x2,受到第j-1次沖擊時的退化量為xj-1,受到第j次沖擊時的退化量為xj。由式(7)可以得到環(huán)境OE2中時間的自然退化量x2小于失效閾值H0的可靠度函數

(23)

同理可得時間t的自然退化量x2小于退化量L的可靠度函數

(24)

將式(23)和式(24)代入式(16),則部件在時間t的可靠度函數可寫為

(25)

其中,假設tk為當前監(jiān)測時間,改變環(huán)境后的部件剩余壽命概率密度函數為

(26)

因此將式(25)代入式(26),可求得環(huán)境改變后部件的剩余壽命概率密度函數。

2 數值試驗分析

本文以風力發(fā)電機組的關鍵部件葉片的退化為例進行考慮環(huán)境變化的核密度剩余壽命預測建模分析。由文獻[23]可知,風力發(fā)電機組葉片的自然退化增量Δx1,Δx2,…,Δxk服從形狀參數是α=0.542,尺度參數是β=1.147的Gamma分布,采用其退化增量數據作為本文部件的監(jiān)測數據,退化失效閾值為20 cm,通過蒙特卡洛分析方法對模型進行驗證。

根據Shafiee[24]對風力發(fā)電機的研究,由蒙特卡洛方法可得沖擊到達的速率λ=0.797 7。主要參數如表1所示。

表1 可靠度函數的參數Tab.1 Parameters of the environment reliability function

2.1 可靠度函數

圖2為部件在環(huán)境OE1的可靠度函數時間T1=0,100,200,300,400,500,600,700,800時由環(huán)境OE1改變?yōu)榄h(huán)境OE2時部件的可靠度函數;部件在OE2環(huán)境中比在OE1環(huán)境中可靠度下降得更快,越早改變環(huán)境,部件退化加快,越早失效。

圖2 在不同監(jiān)測時間T1環(huán)境改變以及不改變環(huán)境的部件的可靠度函數對比Fig.2 Comparison of reliability functions of components with and without environment changes at different monitoring time T1

2.2 剩余壽命預測

設時間T1=100,由Gamma分布的數據增量可得當TH=1 024時,部件失效,圖3為部件在時間T1=100時改變環(huán)境的剩余壽命概率密度函數;表2得出不同監(jiān)測時間的剩余壽命預測均值,隨著部件的不斷運行,監(jiān)測數據會隨之增加,剩余壽命概率密度函數左移的越來越快,這是因為部件受到沖擊之后,退化加速,平均剩余壽命減小。

表2 T1=100時環(huán)境改變的剩余壽命預測值

圖3 T1=100時環(huán)境改變剩余壽命概率密度函數Fig.3 Probability density function of remaining useful life of environmental change at time T1=100

2.3 敏感度分析

本文考慮環(huán)境改變后的競爭性失效模型中,退化和沖擊存在隨機相關性,所以有必要研究降低后的突發(fā)失效閾值D1、沖擊到達的速率λ、沖擊幅值的均值μW這些參數對部件可靠度函數和剩余壽命概率密度函數的影響。

如圖4(a)所示,考慮不同D1對部件退化的影響,得到D1=1.00,1.50,2.00,2.50,3.00,3.54的部件可靠性曲線和剩余壽命概率密度函數;圖4(a)中隨著D1的增大,部件的退化曲線右移,退化變慢;圖4(b)中則顯示為峰值降低,曲線右移,平均剩余壽命增加;在相同的退化條件下,平均剩余壽命降低,表明當部件退化到一定程度,抵抗沖擊的能力降低,很小的沖擊就會導致部件失效。

圖4 不同D1部件的結果圖對比Fig.4 Comparison of different results in different D1

如圖5所示,考慮不同對部件退化的影響,得到λ=0.797 7,2.393 1,3.988 5的部件可靠性曲線和剩余壽命概率密度函數,隨著的增大,環(huán)境OE2逐漸惡劣,沖擊到達的頻率更高,部件退化更快,部件平均失效時間前移。越大的λ值將導致部件在環(huán)境OE2中有越高的突發(fā)失效概率,導致部件更早失效。

圖5 部件在不同λ的結果圖對比Fig.5 Comparison of different results in different λ

為進一步對本文剩余壽命預測方法的效果進行評估,表3給出了在不同沖擊頻率條件下,T1=100時環(huán)境改變的部件不同監(jiān)測時間剩余壽命預測值的比較。由表3可知,相同監(jiān)測時間,隨著λ的增大,部件的剩余壽命預測值減小。

表3 不同λ的剩余壽命預測值Tab.3 Remaining useful life for different λ

如圖6所示,考慮沖擊幅值的均值μW對部件退化的影響,得到μW=0.25,0.50,1.00的部件可靠性曲線和剩余壽命概率密度函數,隨著μW的增大,沖擊的幅值越大,環(huán)境OE2逐漸惡劣,部件退化更快,部件平均失效時間前移。說明μW對環(huán)境OE2中部件可靠性的影響是不可忽略的,越大的μW值意味著每次受到沖擊幅值增大,部件運行環(huán)境更惡劣,將導致部件在環(huán)境OE2中有越高的突發(fā)失效概率,導致部件更早失效。

圖6 部件在不同μW的結果圖Fig.6 Comparison of different results in different μW

表4給出了在不同沖擊幅值的均值μW條件下,T1=100時環(huán)境改變的部件剩余壽命預測值的比較。由表4可知,相同監(jiān)測時間,隨著μW的增大,部件的剩余壽命預測值減小,隨著監(jiān)測時間的增加,可預測不同μW的部件剩余壽命值。即隨著時間的增加,可以預測部件處于不同環(huán)境的平均剩余壽命,驗證了本模型的有效性。

表4 不同μW的剩余壽命預測值Tab.4 Remaining useful life for different μW

如圖7所示,考慮沖擊幅值與造成退化量的均值間的比例關系a對部件退化的影響,由圖可得比例a越大,圖7(a)中可靠度左移,圖7(b)中峰值左移,這是由于沖擊對部件造成的退化量增大,部件失效越快。

圖7 部件在不同a的結果圖比較Fig.7 Comparison of different results in different a

3 實例分析

本文使用模擬風力發(fā)電機齒輪箱的試驗臺進行驗證,如圖8所示,該試驗臺由兩部分組成,分別是主試箱和陪試箱,試驗臺架的中心距為150 mm,試驗采用機械杠桿加載,扭矩采用轉矩轉速傳感器進行測量。試驗過程中對主試箱和陪試箱中的傳感器數據進行監(jiān)測。

圖8 齒輪試驗臺架Fig.8 Gear test bench

主試箱齒輪安裝位置如圖9所示,本試驗臺共安裝了11個傳感器,其中8個傳感器是加速傳感器,由于4#傳感器的安裝位置距離主試箱的軸承近,當傳感器接收到齒輪箱發(fā)生故障的信號時,此時信號衰減最小,可以很好地表征齒輪箱的狀態(tài),故本文選用4#傳感器。其余傳感如9#和10#為聲音傳感器,用于監(jiān)測齒輪箱的聲音信號,而11#為溫度傳感器用于監(jiān)測齒輪箱的溫度。為模擬風機在實際操作過程中環(huán)境和負載的隨機變化情況,試驗加載了八級載荷,從第八級載荷開始預測,在t∈[0,10]h時,齒輪處于嚙齒狀態(tài);t∈[10,68]h時,特征值不斷增大,齒輪進入正常磨損;t∈[68,77]h時,齒輪磨損加劇,在77.17 h發(fā)生斷齒,此時齒輪的故障閾值為y=77.325 mm/s2。由于齒輪在68 h之后磨損突然加劇,建模為隨機沖擊導致的退化加快。在進行數據采樣時,每隔9 min記錄一次采樣數據,每次采樣時長為60 s,采樣頻率為25.6 kHz。

圖9 主箱齒輪安裝位置Fig.9 Location of gear in the main box

在數據采集過程中,由于實際應用中存在著大量未知干擾,獲取的實際數據中也不可避免疊加了很多噪聲,本文使用均方幅值的方法,可以對原始振動數據進行去噪處理的同時并進行特征提取,從而能更好地表征齒輪退化狀態(tài),如圖10所示。采樣信號的均方幅值計算為

圖10 濾除噪聲特征提取后齒輪的退化示意圖Fig.10 Schematic diagram of gear degradation after noise filtering and feature extraction

(27)

式中:N為采樣點數;yM為第M個振動信號的幅值。

3.1 剩余壽命預測

如圖11所示,*為齒輪的真實剩余壽命,°為本文的模型。由圖11可知,當部件不斷運行時,剩余壽命概率密度函數左移的速度越來越快,這表明當部件受到沖擊之后,退化加速,平均剩余壽命減小,由此驗證了本文所提出模型的合理性與有效性。

圖11 環(huán)境改變后齒輪得剩余壽命概率密度函數Fig.11 Probability density function of remaining Useful Life of gears after environmental changes

表5中對不同時刻的齒輪剩余壽命的進行比較,在tk=68 h后,沖擊出現時,由于監(jiān)測數據集不夠,均方根誤差(root mean square error,RMSE)會偏大,隨著數據的增多,本模型的準確性也逐漸提高。

表5 不同監(jiān)測時間下的齒輪箱剩余壽命預測值

4 結 論

本文考慮部件在實際運行中,外界環(huán)境會發(fā)生改變,將惡劣環(huán)境對部件的影響建模為部件受到的隨機沖擊,考慮窗寬對核密度估計的影響,從數據本身出發(fā)使用自適應核密度估計得到自然退化模型;由于環(huán)境改變后部件受到沖擊產生的退化與連續(xù)自然退化過程具有隨機相關性,當退化到一定程度后會降低突發(fā)失效閾值,因此在環(huán)境改變后考慮變失效閾值,引入核密度估計建立考慮沖擊與連續(xù)自然退化過程之間隨機相關的模型,建立了考慮環(huán)境隨機沖擊影響的剩余壽命預測模型,預測部件的平均剩余壽命。最后以風電機組為研究對象,驗證了沖擊的到達率和使突發(fā)失效閾值改變的退化量均影響剩余壽命預測值,即環(huán)境改變后,環(huán)境的惡劣程度對部件的剩余壽命影響無法忽視。

本文研究僅考慮單個部件在環(huán)境改變后的剩余壽命預測,沒有考慮到更為復雜的多部件的情況,多部件系統(tǒng)中各部件之間的隨機相關性也會影響剩余壽命的變化,因此考慮多部件系統(tǒng)存在的隨機相關性,在運行環(huán)境改變后的剩余壽命預測值得進一步研究。