楊 濤,張鈺奇,付春健,趙華東

(1.鄭州大學 機械與動力工程學院,鄭州 450001; 2.河南省智能制造研究院,鄭州 450001)

0 引 言

水工閘門由于工作環(huán)境惡劣,在運行過程中會不可避免地出現(xiàn)性能退化,最終導致破壞。一旦閘門發(fā)生破壞,會造成巨大的經(jīng)濟損失。因此準確預測閘門剩余壽命(Remaining Useful Life,RUL)對于實時評估閘門健康狀態(tài),確定合適的維修策略以保證閘門運行安全具有重要意義。

目前研究主要認為銹蝕是影響閘門壽命的最大因素[1]。Kuz’Mitskii等[2]基于銹蝕分階段模型預測了閘門剩余壽命;李偉康等[3]考慮到銹蝕會影響材料性能,基于D-H曲線對閘門進行RUL預測;然而,目前的研究存在著預測指標單一、精度不高,難以給出剩余壽命的概率分布等問題。

目前,閘門上的銹蝕模型一般為確定性函數(shù),難以體現(xiàn)銹蝕過程的不確定性。而Gamma過程可以很好地模擬銹蝕演化過程[4],體現(xiàn)閘門銹蝕過程的不確定性。因此,考慮采用Gamma過程進行銹蝕模擬,通過數(shù)值模擬方法得到閘門退化數(shù)據(jù)。在獲取閘門的退化數(shù)據(jù)后,同時考慮到工程實際中設(shè)備運行工況復雜,反映設(shè)備退化的指標往往并不唯一。因此,如何將多種信息進行融合來共同反映設(shè)備退化狀態(tài),已成為剩余壽命預測的關(guān)鍵。目前,許多研究通過將多種特征融合成健康因子(Health Index,HI)來解決這一問題。任子強等[5]將多維數(shù)據(jù)融合成一個健康指標來表征發(fā)動機的退化過程,提高了RUL預測的精度;趙廣社等[6]基于歐式距離構(gòu)建了健康因子,并對發(fā)動機進行了RUL預測。

RUL預測的方法可以大致分為基于數(shù)據(jù)驅(qū)動和退化模型兩種[7]。其中,基于退化模型的方法不依賴于大量同類型設(shè)備的全壽命周期數(shù)據(jù),而是根據(jù)歷史數(shù)據(jù)建立退化模型來描述設(shè)備的退化趨勢,因此得到了廣泛的運用[8]。LI等[9]利用維納過程刻畫了渦輪發(fā)動機的退化過程;張英波等[10]基于Gamma過程建立了行星架RUL預測模型。

粒子濾波是基于退化模型實現(xiàn)RUL預測的有效手段,在非線性、非高斯系統(tǒng)中具有明顯優(yōu)勢?；诹Ｗ訛V波的RUL預測方法可以給出剩余壽命的概率分布,改善長期預測的不確定性問題,相較于單一的預測值,更具有參考價值。李玥鋅等[11]基于維納過程對電池退化過程進行建模,并基于粒子濾波算法得到了其剩余壽命概率分布;王璽等[12]基于粒子濾波算法對陀螺儀進行了RUL預測,取得較高的精度。

基于上述分析,本文提出了一種多特征信息融合預測閘門剩余壽命的方法,利用非線性維納過程對閘門進行退化建模,并基于粒子濾波算法實時跟蹤閘門退化曲線并進行RUL預測,同時給出了閘門剩余壽命的概率分布,可以更好地指導工程實際。

1 閘門銹蝕過程模擬

由于閘門運行工況復雜,影響其銹蝕的因素(如水流、干濕循環(huán)等)具有不確定性,因此閘門的銹蝕過程也伴隨著隨機性。此外,銹蝕過程還具有累積性。為了表征銹蝕隨時間發(fā)展的不確定性和累積性,考慮采用Gamma隨機過程模擬銹蝕深度的發(fā)展過程。

Gamma過程是一種具有非負增量的隨機過程,特別適合模擬具有微小增量的累計漸變過程,如銹蝕、疲勞、蠕變等[13]。假設(shè)t時刻的銹蝕量d(t)滿足Gamma分布Ga(α(t),β),其概率密度函數(shù)為

(1)

式中:d為銹蝕量;α(t)和β分別為形狀參數(shù)和尺度參數(shù);Γ(α(t))為Gamma函數(shù)。

根據(jù)Gamma過程的性質(zhì),t時刻銹蝕量的均值可以表示為

E[d(t)] =α(t)β。

(2)

當采用Gamma過程描述退化過程時,一般認為尺度參數(shù)β為常數(shù),形狀參數(shù)α(t)與時間的冪成正比[14],即

α(t)=ctb;c>0,b>0 。

(3)

式中c和b均為常數(shù)。對于特定退化過程b一般為定值,文獻[15] 給出了不同退化類型下b的取值,閘門銹蝕與鋼筋銹蝕屬于同類,取b=1。

當確定參數(shù)b之后,模型中的參數(shù)c和β可以結(jié)合既有結(jié)構(gòu)的銹蝕觀測數(shù)據(jù)采用極大似然估計法或矩估計法來進行估計。然而,由于缺乏閘門長期銹蝕監(jiān)測數(shù)據(jù),同時考慮到t時刻的銹蝕深度與平均銹蝕速率有關(guān),因此利用已有銹蝕率的統(tǒng)計數(shù)據(jù)對參數(shù)c和β進行反向標定[16],即首先假定參數(shù)β為某個固定值,根據(jù)式(2)并結(jié)合銹蝕速率統(tǒng)計數(shù)據(jù)即可確定參數(shù)c。

當確定參數(shù)之后,閘門銹蝕深度發(fā)展過程便可以采用Gamma順序采樣法進行模擬:

(1)將總時間T等分為n個子區(qū)間。

(2)當t=ti時,生成服從Ga(α(ti)-α(ti-1),β)的隨機數(shù)Δdi,Δdi便為時段ti-ti-1的銹蝕增量。則ti時刻的銹蝕累積深度為di=di-1+Δdi。

(3)重復以上過程,直至達到總時間T。

2 健康因子的構(gòu)建

閘門的退化過程可由多種特征參數(shù)變化來表示,不同的特征對于閘門退化過程的敏感性不同,并且各個特征之間存在一定的相關(guān)性。若將多個特征指標進行單獨建模,不僅不符合多變量相互耦合實際情況,也會增加建模的復雜性。因此,對于多特征的情況,常常采用信息融合的方法,利用降維等手段將多維數(shù)據(jù)轉(zhuǎn)化為一維復合健康指標,這種方法不僅充分考慮了多變量相互關(guān)聯(lián)的情況,也更有利于對閘門退化過程進行建模。因此,將多種特征指標進行融合得到能夠表征閘門退化過程的健康因子是進行壽命預測的關(guān)鍵。

2.1 數(shù)據(jù)獲取

閘門銹蝕會使構(gòu)件截面面積變小,空間結(jié)構(gòu)整體效應發(fā)生改變,從而影響閘門的強度和剛度。同時,閘門銹蝕會引起自振頻率和振型等模態(tài)參數(shù)的改變。為此,本文綜合考慮了閘門靜力學和動力學特性,選取了應力、變形、主梁撓度,“干模態(tài)”下前10階自振頻率、位移模態(tài)、應變模態(tài),以及考慮流固耦合作用下的“濕模態(tài)”前10階自振頻率、位移模態(tài)、應變模態(tài),共63種特征指標來共同反映閘門退化過程。

2.2 特征指標篩選

特征篩選是指從特征集中選擇出最優(yōu)的子集來反映閘門退化過程。考慮到單調(diào)性強的特征能夠很好地表征閘門的退化趨勢,因此首先對特征進行單調(diào)性篩選。Spearman相關(guān)系數(shù)可用于評估2個變量之間相關(guān)關(guān)系的程度,考慮到時間序列是嚴格單調(diào)的,因此計算每一種特征時序與對應時間序列的Spearman作為該種特征的單調(diào)性得分。同時,考慮到變化程度大的特征能更好地反映退化過程。因此,綜合考慮單調(diào)性和離散性對特征指標進行篩選。

對于某一特征參數(shù)序列X=(x1,x2,…,xn),相應的時間序列T1=(t1,t2,…,tn),則:

(1)單調(diào)性指標為

(4)

式中:n為數(shù)據(jù)的個數(shù);分別將X,T1中的元素按照從小到大的順序編秩;Ri表示xi在X中的秩次;Qi表示ti在T1中的秩次。

(2)離散性指標為

為發(fā)現(xiàn)最優(yōu)特征指標,首先選擇出單調(diào)性較強的特征指標,然后再對離散性高的特征指標進行選擇。

2.3 特征融合及健康因子的構(gòu)建

多特征參數(shù)可以很好地表征設(shè)備退化狀態(tài),也會使設(shè)備建模和數(shù)據(jù)處理復雜化。因此,要將多維數(shù)據(jù)融合成一維數(shù)據(jù),即通過采用健康因子來表征閘門的退化過程。

主成分分析(Principal Component Analysis,PCA)是一種數(shù)據(jù)融合和特征提取的重要方法,主要思想是通過尋找一組正交基,將高維特征映射到低維空間,從而達到降維的目的。其主要步驟為:

(1)將性能退化指標構(gòu)成特征數(shù)據(jù)矩陣Y為

(7)

(3)計算P的特征值m1,m2,…,mp并從大到小進行排序,計算對應的特征向量e1,e2,…,ep。

(4)根據(jù)特征值mq,計算第q個主成分對應的貢獻率

(8)

一般取貢獻率ξ>80%的主成分可以較全面地反映出特征的信息[17],將其當做特征融合后的健康因子來反映設(shè)備退化過程。

3 閘門壽命預測

3.1 粒子濾波原理

粒子濾波是一種基于蒙特卡洛(Monte Carlo,MC)理論的貝葉斯估計算法,可用于處理非線性、非高斯系統(tǒng)的狀態(tài)估計問題。貝葉斯算法是粒子濾波的基礎(chǔ),通常采用狀態(tài)轉(zhuǎn)移方程和觀測方程來描述系統(tǒng)的狀態(tài)空間

(9)

式中:xk-1、xk分別表示k-1和k時刻系統(tǒng)的狀態(tài)值;fk表示系統(tǒng)的狀態(tài)傳遞函數(shù);εk-1表示系統(tǒng)的過程噪聲;yk表示k時刻系統(tǒng)的觀測值;hk表示觀測函數(shù);vk表示系統(tǒng)的測量噪聲。

貝葉斯濾波是對系統(tǒng)的后驗分布進行求解,首先根據(jù)狀態(tài)轉(zhuǎn)移方程,根據(jù)k-1時刻的觀測值,可預測得到k時刻系統(tǒng)狀態(tài)的先驗分布

p(xk|y1:k-1)=

(10)

式中:p(xk|xk-1)由狀態(tài)轉(zhuǎn)移方程定義;y1:k-1為1:(k-1)時刻的系統(tǒng)觀測值。

在此基礎(chǔ)上,根據(jù)k時刻觀測值,利用貝葉斯理論得到k時刻的后驗概率密度分布

式中p(yk|xk)由狀態(tài)測量方程定義。

由式(10)、式(11)可知,上述過程存在著復雜的積分運算,往往難以得到精確解。為了解決上述問題,粒子濾波算法借鑒蒙特卡羅方法將k時刻的后驗分布離散加權(quán)為

(12)

由式(12)可知,得到粒子對應的權(quán)值分布,就可以實現(xiàn)對狀態(tài)概率分布的近似估計,假設(shè)q(x0:k|y1:k)表示過去時刻狀態(tài)的后驗估計,同時狀態(tài)和觀測信息遵循一階馬爾科夫過程,則序貫重要性采樣粒子權(quán)重可轉(zhuǎn)換為

(13)

然而,隨著迭代次數(shù)的增加,權(quán)重有可能只存在少數(shù)粒子中,大部分粒子的權(quán)重可能小到忽略不計,造成估計性能的下降。因此,常常通過采用重采樣方法降低粒子退化問題。

3.2 退化模型的建立

閘門在運行過程中由于工作環(huán)境惡劣,其退化過程呈現(xiàn)出明顯的不確定性。Wiener過程可以很好地描述復雜設(shè)備退化過程中的不確定性,在RUL預測中得到了廣泛運用?？紤]到以上因素,利用非線性維納過程來刻畫閘門隨機退化過程,即

X(t)=X(0)+λeθt+σBB(t) 。

(14)

式中:X(t)表示t時刻的系統(tǒng)狀態(tài);X(0)為系統(tǒng)初始狀態(tài);λ為漂移系數(shù);eθt為時間t和參數(shù)θ的非線性函數(shù),用于表示退化過程中的非線性;σB為擴散系數(shù);B(t)為標準維納過程。

因此,閘門的退化過程可以描述為

xk=xk-1+λk-1(eθk-eθ(k-1))+ηk。

(15)

式中ηk=σB(B(tk)-B(tk-1))服從于正態(tài)分布N(0,σ2Δt),Δt=tk-tk-1。

事實上,由于測量結(jié)果與真實狀態(tài)之間不可避免地存在誤差,因此,測量狀態(tài)和真實狀態(tài)之間的關(guān)系可以表示為

yk=xk+v。

(16)

式中v為測量噪聲,服從于正態(tài)分布N(0,γ2)。

同時,由于過程噪聲可以由其他參數(shù)的不確定性實現(xiàn)[18]。因此,根據(jù)退化模型可以得出其狀態(tài)空間模型,即

由式(17)可得,利用Yk=[xk,λk,θk,σk,γk] 可表示出系統(tǒng)的狀態(tài),并結(jié)合第3.1節(jié)所敘述的粒子濾波原理,利用在線數(shù)據(jù)實時更新未知參數(shù),便可以實現(xiàn)RUL預測。

3.3 剩余壽命預測

在得到狀態(tài)空間方程之后,根據(jù)第3.2節(jié)所述的粒子濾波算法對模型參數(shù)進行更新,可以得到t時刻的系統(tǒng)狀態(tài),根據(jù)退化模型進行遞推,可以對t+φ時刻的系統(tǒng)狀態(tài)進行預測[19]。退化模型遞推公式為

1≤ζ≤φ。

(18)

圖1 RUL預測示意圖Fig.1 Schematic diagram of RUL prediction

4 案例驗證

4.1 有限元模型建立

某溢洪道露頂式弧形鋼閘門尺寸為12 m×10.5 m(寬×高),主要材料為Q345B,彈性模量取2.06×105N/mm2,密度取7 850 kg/m3,泊松比取0.3,作用水頭取9.4 m,主梁計算跨度為7.6 m,現(xiàn)場原型如圖2所示。閘門為空間薄壁結(jié)構(gòu),采用Shell181殼單元模擬。閘門濕模態(tài)仿真時,上游水體長度取超過閘門高度10倍的長度,取150 m,采用Fliud30單元模擬,水流和閘門的相互作用表現(xiàn)在閘門迎水面上,將其設(shè)為流固耦合面。

圖2 閘門現(xiàn)場Fig.2 Site photo of the gate

4.2 閘門銹蝕深度演化過程模擬

根據(jù)文獻[20] 提出的銹蝕線性模型,閘門上的銹蝕深度變化可根據(jù)式(19)計算,即

(19)

式中:dL(t)為閘門t時刻的銹蝕量;t為閘門運行年限;v0為銹蝕速率;tI為閘門開始銹蝕的年限,一般取3～7 a,本文取5 a。

閘門上的銹蝕分為全面銹蝕和局部銹蝕,根據(jù)工程實際,全面銹蝕常常發(fā)生在面板、支臂、縱梁、肋板、次梁等結(jié)構(gòu)上,局部銹蝕常常發(fā)生在下主梁積水處、面板邊側(cè)和邊梁漏水處。根據(jù)統(tǒng)計我國閘門銹蝕速率的均值在(0.015～0.070)mm/a范圍內(nèi),考慮面板和局部銹蝕區(qū)域與水長期接觸,因此對面板和局部銹蝕區(qū)域銹蝕速率取v1=0.07 mm/a,對于支臂、縱梁等其他區(qū)域取v2=0.06 mm/a。

假定尺度參數(shù)β=0.1,結(jié)合式(2)、式(19)便可以得到形狀參數(shù)α(t)為

α(t)=c(t-5),t>5 。

(20)

式中:v1=0.06 mm/a時,c=0.6;v2=0.07 mm/a,c=0.7。

確定形狀參數(shù)α(t)和尺度參數(shù)β后,便可以利用Gamma順序采樣法模擬閘門銹蝕發(fā)展過程。值得說明的是由于銹蝕過程存在著不確定性,利用Gamma順序采樣法會產(chǎn)生多條銹蝕曲線,通過實際檢測的銹蝕數(shù)據(jù)的增加,利用參數(shù)估計方法可減少這種不確定性。為了數(shù)據(jù)的可重復性和合理性,本文采用了設(shè)置隨機種子的方法,找到與線性模型相近的曲線,既能保證其銹蝕過程的不確定性,又保證其合理性,得到的銹蝕發(fā)展過程如圖3所示。圖3中線性增長模型為文獻[20] 提出的銹蝕模型,其無法體現(xiàn)銹蝕過程中的不確定性;考慮銹蝕深度發(fā)展過程的不確定性,利用Gamma過程模擬銹蝕深度發(fā)展過程時,其是由一系列不同長度的微小增長段構(gòu)成。

圖3 銹蝕深度演化過程Fig.3 Evolution process of corrosion depth

根據(jù)文獻[21] ,當最大應力發(fā)生處厚度≤16 mm時,取[σ] =225 MPa。通過對不同運行年限下的閘門進行仿真,結(jié)果表明在29 a時閘門上最大應力σmax>0.95[σ] =213.75 MPa達到失效的條件,圖4為閘門失效時刻應力云圖。

圖4 閘門應力云圖Fig.4 Stress contours of gate

4.3 閘門特征指標篩選

利用第2.2節(jié)所敘述的單調(diào)性和離散性指標對特征進行篩選,特征指標得分如圖5所示。圖5中橫坐標1～3分別為應力、變形、主梁撓度,4～33分別為干模態(tài)下前10階自振頻率、位移模態(tài)、應變模態(tài),34～63分別為濕模態(tài)下前10階自振頻率、位移模態(tài)、應變模態(tài)。

圖5 特征指標得分Fig.5 Scores of characteristic indices

綜合考慮,先取單調(diào)性得分>0.9的特征指標,然后再從單調(diào)性篩選之后的指標中進行離散性篩選?？紤]到10階之后的離散性得分趨于平穩(wěn)且數(shù)值較小,對健康因子的構(gòu)建影響不大,因此取離散性得分前10特征用于融合健康因子,共篩選出應力、變形干模態(tài)下3、4、5、8階位移模態(tài),濕模態(tài)下4、5、8階位移模態(tài),5階應變模態(tài)在內(nèi)的10種特征,10種特征的退化曲線如圖6所示。閘門所篩選出的變形、位移模態(tài)、應變模態(tài)等特征的初始狀態(tài)云圖如圖7所示。(由于應力圖4已給出,本處不再列出)。

圖6 特征指標退化曲線Fig.6 Degradation curves of feature indices

圖7 特征指標云圖Fig.7 Contours of feature indices

4.4 閘門健康因子的構(gòu)建

根據(jù)2.3節(jié)所論述的PCA方法對篩選后的特征進行融合以構(gòu)建健康因子,利用PCA方法得到的各主成分分量得分如圖8所示。

圖8 各主成分貢獻率Fig.8 Contribution rate of each principal component

其中,第一主成分貢獻率ξ=90.960 1%,滿足ξ>80%的要求,因此,取第一主成分構(gòu)建健康因子,健康因子曲線,如圖9所示。

圖9 健康因子曲線Fig.9 Health index curve

4.5 閘門剩余壽命預測

在得出健康因子退化曲線之后,便可以利用粒子濾波方法實時跟蹤退化曲線,并進行RUL預測。在利用粒子濾波方法進行壽命預測前,需先確定粒子的初始分布。然而事實上,對于一個特定的系統(tǒng)很難進行足夠多的壽命測試以獲取可靠的參數(shù)先驗分布[22]。因此,參數(shù)的先驗分布一般是根據(jù)經(jīng)驗得出,通過對閘門腐蝕退化過程進行反復試驗并結(jié)合先驗知識,發(fā)現(xiàn)初始分布為均勻分布更加符合閘門退化過程,具體參數(shù)如表1所示。

表1 初始參數(shù)分布Table 1 Initial parameter distribution

利用第3.1節(jié)所敘述的粒子濾波原理對閘門進行RUL預測,由于粒子濾波中的粒子是從初始分布中隨機產(chǎn)生的,導致每次循環(huán)預測結(jié)果在微小范圍內(nèi)波動。為了更好地說明粒子濾波進行RUL預測的原理和步驟,圖10—圖12是粒子濾波單次循環(huán)預測的結(jié)果。粒子濾波實時跟蹤退化曲線如圖10所示,由圖10可以看出所提方法能很好地跟蹤退化曲線,并能有效地預測閘門RUL。

圖10 實時跟蹤退化曲線Fig.10 Tracking of degradation process

狀態(tài)跟蹤集(即不同的預測時刻)和粒子數(shù)是衡量粒子濾波算法有效性的關(guān)鍵參數(shù)。表2為算法在不同的組合方案下的仿真,仿真次數(shù)為100次,從表2可看出在不同組合方案下,該算法的預測精度均較高,表明該算法具有較強的穩(wěn)定魯棒性、泛化適用性。

表2 不同組合下RUL預測結(jié)果對比Table 2 Comparison of RUL prediction results under different combinations

綜合考慮預測精度和時效性,將粒子數(shù)設(shè)為1 000,則在每一個預測時刻可以得到1 000個RUL預測值,從而得出RUL分布直方圖,部分預測時刻的RUL分布頻率直方圖如圖11所示。將粒子濾波預測中值當做預測值,由圖11可知,在95%的置信區(qū)間內(nèi),預測值與真實值吻合程度較高。隨著閘門運行年限的增加,在不同預測時刻可以得到預測的剩余壽命概率密度函數(shù),如圖12所示。從圖12可以看出,隨著閘門運行年限的增加,剩余壽命的概率分布區(qū)間范圍在逐漸減小,說明預測的不確定性在逐漸降低。

圖11 剩余壽命分布直方圖Fig.11 Histogram of RUL distribution

圖12 剩余壽命概率密度函數(shù)Fig.12 Probability density function of RUL

4.6 RUL預測結(jié)果及分析

為了驗證所提算法的穩(wěn)定性以及預測性能的優(yōu)劣水平,對每一預測時刻進行100次重復試驗。同時為了更好地指導工程實際,給出了各個預測剩余壽命區(qū)間出現(xiàn)的概率,部分預測時刻結(jié)果如表3所示,表3中預測值是指100次循環(huán)下預測RUL的均值,Var是指100次循環(huán)下預測RUL的方差,AE代表預測值與實際RUL的絕對誤差。

表3 100次重復試驗結(jié)果Table 3 Results of 100 replicate tests

從閘門運行第10年開始預測,按照上述方法可得到閘門不同運行年限下的RUL預測結(jié)果,如圖13所示。由圖13可以看出,RUL預測結(jié)果與實際RUL結(jié)果相近,且真實RUL均在95%的預測RUL置信區(qū)間內(nèi)。

圖13 剩余壽命預測結(jié)果Fig.13 RUL prediction results

為了客觀衡量預測的精度,引入均方根誤差(Root Mean Square Error,RMSE)、平均絕對誤差(Mean Absolute Error,MAE)和方差絕對誤差(Variance Absolute Error,VAE)來進行評價(式(21)),評價指標得分越小說明預測的精度越高,結(jié)果如表4所示。從表4可以看出3個評價指標得分均比較小,說明所提方法預測閘門RUL具有較高的精度。

表4 評價指標得分Table 4 Scores of evaluation indices

5 結(jié) 論

(1)利用應力、自振頻率、位移模態(tài)、應變模態(tài)等多特征參數(shù)相較于單一指標更能全面、準確地反映閘門退化過程。

(2)利用非線性維納過程對閘門退化過程進行建模,并基于粒子濾波對其進行剩余壽命預測,可以得到剩余壽命的概率分布以描述預測結(jié)果的不確定性,從而更好地指導閘門日常管理。同時,由評價指標RMSE為1.395 5,說明所提方法具有較高的預測精度。

(3)未來根據(jù)剩余壽命預測結(jié)果,綜合考慮經(jīng)濟性、安全性等因素以確定合理的維修策略是下一步研究的關(guān)鍵。