基于多特征融合的水工閘門剩余壽命預測

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

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

0 引 言

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

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

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

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

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

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

1 閘門銹蝕過程模擬

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

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

(1)

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

根據Gamma過程的性質,t時刻銹蝕量的均值可以表示為

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

(2)

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

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

(3)

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

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

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

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

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

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

2 健康因子的構建

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

2.1 數據獲取

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

2.2 特征指標篩選

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

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

(1)單調性指標為

(4)

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

(2)離散性指標為

為發現最優特征指標,首先選擇出單調性較強的特征指標,然后再對離散性高的特征指標進行選擇。

2.3 特征融合及健康因子的構建

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

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

(1)將性能退化指標構成特征數據矩陣Y為

(7)

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

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

(8)

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

3 閘門壽命預測

3.1 粒子濾波原理

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

(9)

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

貝葉斯濾波是對系統的后驗分布進行求解,首先根據狀態轉移方程,根據k-1時刻的觀測值,可預測得到k時刻系統狀態的先驗分布

p(xk|y1:k-1)=

(10)

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

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

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

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

(12)

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

(13)

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

3.2 退化模型的建立

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

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

(14)

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

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

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

(15)

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

事實上,由于測量結果與真實狀態之間不可避免地存在誤差,因此,測量狀態和真實狀態之間的關系可以表示為

yk=xk+v。

(16)

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

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

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

3.3 剩余壽命預測

在得到狀態空間方程之后,根據第3.2節所述的粒子濾波算法對模型參數進行更新,可以得到t時刻的系統狀態,根據退化模型進行遞推,可以對t+φ時刻的系統狀態進行預測[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,現場原型如圖2所示。閘門為空間薄壁結構,采用Shell181殼單元模擬。閘門濕模態仿真時,上游水體長度取超過閘門高度10倍的長度,取150 m,采用Fliud30單元模擬,水流和閘門的相互作用表現在閘門迎水面上,將其設為流固耦合面。

圖2 閘門現場Fig.2 Site photo of the gate

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

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

(19)

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

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

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

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

(20)

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

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

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

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

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

4.3 閘門特征指標篩選

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

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

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

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

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

4.4 閘門健康因子的構建

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

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

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

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

4.5 閘門剩余壽命預測

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

表1 初始參數分布Table 1 Initial parameter distribution

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

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

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

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

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

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

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

4.6 RUL預測結果及分析

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

表3 100次重復試驗結果Table 3 Results of 100 replicate tests

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

圖13 剩余壽命預測結果Fig.13 RUL prediction results

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

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

5 結 論

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

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

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