基于MWSA的熱力系統單參數時序預測方法

肖鵬飛， 倪 何， 金家善

(海軍工程大學 a. 動力工程學院； b. 船舶與海洋工程學院，武漢 430033)

熱力系統具有部件眾多、結構復雜和工況多變等特點，設備故障的耦合性強且排故成本高[1].為了提高系統的全壽命健康管理水平，如何通過系統的歷史運行參數預測系統未來一段時間內的運行狀態，并將得到的預測結果作為系統運行管理策略和裝備維修的參考，已成為當前熱力系統科學管理的重要研究方向之一[2-3].

熱力系統的運行參數大多為非平穩時間序列.在非平穩時間序列趨勢項的提取方面，文獻[4]對風速時序數據進行預測并基于排列熵(PE)的經驗模態分解(EMD)趨勢項提取算法，結合人工神經網絡(ANN)的預測模型對提取的趨勢進行預測，通過計算預測數據與實際數據的相關系數來驗證預測精度，相比未引入排列熵算法的EMD結合ANN的EMD-ANN算法具有更好的預測效果；陳亮等[5]提出基于集合經驗模態分解(EEMD)和奇異值分解(SVD)的排列熵非平穩時序趨勢提取方法，并通過與文獻[6]和文獻[7]中的方法進行對比，驗證了該算法的優越性，但由于EEMD分解后的各個本征模態函數(IMF)含噪信息不一，且在分解后的噪聲分量和信號分量的分界問題方面選擇直接去掉最高頻分量的方法而未做量化分析，導致最高頻IMF分量中的信號信息丟失，其余分量未進行降噪處理，使后續趨勢提取誤差增加；劉海江等[8]提出利用EMD結合小波閾值去噪(WTD)的降噪方法，在保留原信號信息的基礎上對噪聲進行了更有效地消除，但因EMD算法中的“端點效應”問題導致重構信號誤差較大；黃禮敏[9]在EMD算法的基礎上提出了中值回歸經驗模態分解(MREMD)方法，有效解決了EMD分解過程中的“端點效應”問題，并且優化了包絡線生成方式，但同樣未對分解后的IMF分量進行定量分界，導致后續計算結果誤差較大.Behera等[10]利用整合滑動平均自回歸(ARIMA)模型，通過對公共部門每天接到的電話次數進行預測，并與實際數據對比分析，驗證了模型的有效性，其預測結果可作為部門管理層規劃員工日常工作量的參考.但由于熱力系統中的參數歷史數據中包含較多干擾，直接采用ARIMA模型預測誤差較大；仇琦等[11]在ARIMA預測模型的基礎上提出了基于EMD結合ARIMA的EMD-ARIMA算法預測模型，通過對光電功率信號經EMD分解后的預測驗證了模型的有效性，但這兩種方法都未對原始非平穩信號的趨勢進行提取，導致預測干擾較大.

基于上述研究提出一種熱力系統單參數時序預測模型，由MREMD、WTD、基于奇異值分解和排列熵的K-means聚類算法(SVD-PE-KCA)和ARIMA共4部分組成，簡稱MWSA.利用對原始信號的分解、降噪和趨勢提取，使預測結果更加精確，可為熱力系統的運行管理、控制優化和維修保障研究提供底層輸入，具有一定的理論創新和工程應用價值.

1 預測步驟和算法

基于MWSA的單參數時序預測流程如圖1所示.

首先，利用MREMD算法將歷史運行參數的時序數據信號分解為若干個IMF分量和1個殘余分量(Rn)；其次，計算各個IMF分量與原信號的互相關系數和均方根誤差，篩選出不符合閾值條件的分量進行WTD，并將去噪后的IMF分量與原本符合篩選條件的分量以及殘余分量重組成新的IMF分量；再次，利用SVD將重組后的新IMF分量分解，重構成按照奇異值升序排列的奇異值分量，并采用基于優化參數的PE算法計算各個IMF分量的熵值，利用K-means聚類算法針對各IMF分量熵值進行分類；最后，選取熵值較低的一類IMF分量重構為當前運行參數的變化趨勢項，并采用ARIMA模型對其進行預測.

1.1 MREMD分解

MREMD采用自回歸(AR)模型對信號端點延拓，并用優化包絡線擬合的方法改善了EMD算法中的“端點效應”問題，具體步驟和算法如下.

步驟一對于長度為N的初始時間序列x0(t)，采用如下式所示的p階AR模型對x(t)的左右端點進行預測延拓，即x(t)為延拓后的新序列，使左右端點處于延拓后時間序列的相鄰兩個極值點之間.

xt=φ0+φ1xt-1+φ2xt-2+…+φpxt-p+μt

(1)

式中：xt(t=p+1,p+2, …,N)為均值點；φ0，φ1，…，φp為p+1個實數；μt為零均值的白噪聲序列.

h1,0(t)=x0(t)-m1,0(t)

(2)

步驟三將計算得到的信號分量h1,0(t)作為原始信號，重復步驟一和步驟二進行迭代計算，直到經過l次迭代后的信號分量h1,l(t)滿足如下式所示的終止條件.

(3)

z=1, 2, …,k+1}≥68.5%

(4)

(5)

步驟四將h1,l(t)作為1階本征模態函數(即I1分量)輸出，x0(t)-I1得到分量R1，并將R1作為原始信號重復步驟一至步驟四，直到殘余信號成為單調函數或分離不出新的IMF分量為止，殘余信號提取過程表示為

(6)

式中：下標n為能夠分解出的IMF分量個數.

經過上述分解后，原始信號x0(t)可表示為所有IMF分量和最終殘余分量之和，表達式為

(7)

1.2 IMF分量選取

對于MREMD分解后的IMF分量，噪聲主要集中在高頻分量中，在工程實際應用時一般通過舍棄高頻分量的方法來達到降噪的目的，但這同樣也會丟棄部分的信號信息，使得分析結果產生誤差.此外，雖然低頻分量通常包含了較多的信息成分，但同樣也存在一些噪聲，僅靠舍棄高頻分量的方法并不能夠實現完全降噪的目的.為此，本文采用王彬蓉等[13]提出的IMF分量篩選方法對所有分量進行篩選，在降噪的同時盡可能多地保留原始信號的信息，具體步驟和算法如下.

步驟一計算各階IMF分量與原始信號的相關系數ci和均方根誤差RMSEi，表達式為

(8)

(i=1, 2, …,n)

步驟二計算IMF分量相關系數篩選閾值和均方差篩選閾值，表達式為

(9)

步驟三選取滿足如下條件的IMF分量，進行小波閾值降噪為

ci≤λIMF， RMSEi≥δIMF

(10)

1.3 WTD

WTD是利用小波變換將原始信號分解成多層近似系數和細節系數，由于經過小波變換后的噪聲信息主要集中在絕對值較小的細節系數中，所以可通過將絕對值小于規定閾值的細節系數設置為0，并將剩余小波系數(即分解得到的近似系數和保留的細節系數)通過小波逆變換重構回原始信號的方法，以達到去除噪聲的目的，具體步驟和算法如下.

步驟一閾值選取.傳統WTD的閾值選取方法有固定閾值法、自適應閾值法、啟發式閾值法和極大極小閾值法等[14].本文采用改進的固定閾值法，通過細節系數對每一層的噪聲進行估計，同時利用系數ln(j+1)逐層降低細節系數的閾值[15]，從而盡可能多地保留蘊含在高頻分量中的真實信號，其閾值選取準則為

(11)

式中：λj為第j層小波細節系數的噪聲閾值；d(j)為第j層小波細節系數；median為中間值函數，即將每層系數按照降序排列后，取其中間數的值(當系數個數為奇數時)或中間兩個數的均值(當系數個數為偶數時).

步驟二閾值處理.常用閾值處理函數有硬閾值函數、軟閾值函數和復合閾值函數[16]等，依次表示為

1.4 趨勢項提取

采用SVD、PE和K-means聚類算法的組合以實現信號的進一步降噪和趨勢項的提取，具體步驟和算法如下.

步驟一SVD分解.將經過小波閾值降噪處理后的n個IMF分量I1～In和殘余分量Rn去均值，然后以列向量的形式組成矩陣AN×(n+1)，并計算其協方差矩陣B(n+1)×(n+1)：

AN×(n+1)=

(12)

(13)

將矩陣U和A中0特征值所對應的向量去掉，重構后的K個奇異值分量表示為

IN×K=AN×KUK×K=(P1P2…PK)

(14)

式中：PK為第K個奇異值分量.

步驟二PE計算.利用延遲相空間重構法[17]對IN×K中各個奇異值分量進行重構，如對分量P1重構表示為

(15)

式中：l=N-(m-1)τ；p1, i(i=1, 2, …,N)為第1個奇異值分量的第i個元素；m為嵌入維數；τ為延遲時間.

將重構矩陣的每行元素按升序排列，得到每行元素的索引值，然后將每行元素的索引值組成索引向量并構成索引矩陣.索引矩陣中每行元素的排列方式共有m!種可能，假設第b種排列方式在索引矩陣中的出現次數為ub，統計其出現的頻率為

fb=ub/m!，b∈(1,K)

(16)

由式(16)可計算得到各奇異值分量的排列熵為

(17)

步驟三K-means聚類.根據計算得到的排列熵值對IMF分量進行聚類，選取熵值較低的分量重構原始信號的趨勢項.首先，任意選取其中1個點ρ1作為1個類別中心，并選擇與ρ1距離最遠的ρ2點為另1個類別中心，計算熵值序列中其余對象點x與ρ1和ρ2的距離為

d(ρ,x)=|x-ρ|

(18)

然后，把其余各點劃分到與兩個類別中心距離近的一類，將兩類點的均值作為新的類別中心進行分類，重復上述步驟，直至兩個類別中心不再發生變化為止.最后，將兩類熵值點中平均熵值較低一類所對應的奇異值分量篩選出來[18]，重構為原始信號的趨勢項為

(19)

1.5 ARIMA預測模型

ARIMA常用于差分平穩時間序列的擬合和預測，具體步驟和算法如下.

步驟一對式(19)計算得到的趨勢項TN×1進行平穩性(Phillips-Perron)檢驗[19]，然后對檢驗結果為非平穩的時間序列作差分平穩化處理，表示為

(20)

(21)

步驟三經過AIC檢驗得到模型的自回歸系數多項式階數s和移動平均系數多項式階數q，并擬合ARIMA(s,d,q)模型表達為

yt=c+a1yt-1+…+asyt-s+εt-

b1εt-1-…-bqεt-q

(22)

步驟四對ARIMA模型的殘差序列進行LB(Ljung-Box)統計量檢驗[21]，表示為

(23)

?k′>0

2 預測案例

以某型船舶蒸汽動力系統的除氧器水位為案例，對上述基于MWSA的熱力系統單參數時序預測方法進行驗證.按1.1節的步驟和算法，對一段時間內監控系統采集的實際水位數據進行MREMD分解，得到7個IMF分量和1個殘余分量，如圖2所示.

圖2 分解得到的IMF分量和RCFig.2 Decomposed IMF component and RC

對于分解得到的分量，根據式(8)計算其與原始信號的相關系數和均方根誤差，如表1所示.

表1 各個IMF分量與原信號的相關系數和均方根誤差Tab.1 Correlation coefficient and root mean square error of each IMF component relative to the original signal

根據式(9)計算得到IMF分量相關系數的篩選閾值λIMF=0.301 7，均方差的篩選閾值δIMF=7.314 1，符合條件的分量有I4～I7.按1.3節的具體步驟和算法對I1～I3進行小波閾值降噪，經反復測試，選取Sym4小波基，分解層數為3層.采用固定閾值選取方法，根據式(11)計算各層的噪聲閾值.采用不同閾值函數處理后的各IMF分量的信噪比及均方根誤差，如表2所示.

表2 去噪后各IMF分量的信噪比、均方根誤差Tab.2 Signal-to-noise ratio and root mean square error of each IMF component after denoising

由表2可知，對于高頻的I1分量，采用復合閾值降噪的信噪比更高，均方差更小.但對于較低頻的I2和I3分量，采用硬閾值降噪的效果更佳.按1.4節的步驟和算法，將降噪后的I1～I3與I4～I7及殘余分量Rn一起進行SVD分解，得到8個奇異值分量即P1～P8，如圖3所示.

圖3 SVD分解后的奇異值分量Fig.3 Singular value components after SVD decomposition

根據式(15)～(17)計算各奇異值分量的排列熵值.在相空間重構時，運用互信息法[22]選取最佳延遲時間，在1～100 s的延遲時間范圍內，各奇異值分量的互信息值Mi1～Mi8變化如圖4所示.采用偽近鄰法[23]計算最佳嵌入維數，在1～10的嵌入維數變化范圍內，各奇異值分量的偽近鄰率FNNP1～FNNP8隨嵌入維數m的變化如圖5所示.

圖4 各分量互信息值隨延遲時間的變化曲線Fig.4 Variation curve of mutual information value with delay time

圖5 各分量偽近鄰率隨嵌入維數的變化曲線Fig.5 Variation of pseudo nearest neighbor rate with embedding dimension

取圖4中互信息值隨時間變化曲線的第一個極值點為各奇異值分量的最佳延遲時間τopt；取圖5中偽近鄰率小于5%時的嵌入維數為各奇異值分量的最佳嵌入維數mopt；然后將其代入式(15)，并通過式(16)和式(17)計算得到各分量的排列熵值，如表3所示.

表3 各奇異值分量的最佳延遲時間、最佳嵌入維數和排列熵值Tab.3 Optimal delay time, embedding dimension, and permutation entropy of each singular value component

對各分量的排列熵值進行K-means聚類，選取P2、P3和P4共3個分量進行趨勢項重構，結果如圖6所示.

根據1.5節的步驟和算法，對除氧器水位進行擬合和預測.首先，對如圖6所示的趨勢項進行Phillips-Perron檢驗，結果為非平穩序列，再經過4次差分計算后轉換為平穩序列；然后，對差分平穩化后的序列進行AIC，得到自回歸系數多項式階數s以及移動平均系數多項式階數q分別為1和3；在經過4次反差分計算和LB統計量檢驗后，擬合得到如式(22)所示的ARIMA(1, 4, 3)模型.

圖6 除氧器水位的原始項與趨勢項Fig.6 Original item and trend item of deaerator water level

取實際運行數據的前90%為訓練數據，取后10%為檢驗數據，預測除氧器水位在未來94 s的變化趨勢，結果如圖7所示.而在相同數據輸入下，采用文獻[23]中提出的無降噪處理的算法(MSOP)得到的預測結果如圖8所示.

由圖7和圖8 可知，隨著預測時長的增加，兩種方法的預測誤差都會增大.經計算，在采用MWSA方法時，預測結果與實際數據(檢驗數據)的均方差為1.388 7；而文獻[23]中提出的方法其預測結果與實際數據(檢驗數據)的均方差為2.133 4.可知，本文預測方法的準確度更高.

圖7 基于MOSP的除氧器水位趨勢項與預測項Fig.7 Trend and prediction of deaerator water level based on MSOP

圖8 基于MWSA的除氧器水位的趨勢項與預測項Fig.8 Trend and prediction of deaerator water level based on MWSA

3 結語

提出一種基于MWSA的熱力系統單參數時序預測方法，并以某型號船舶蒸汽動力系統的除氧器水位為例，根據實際運行數據預測其在未來一段時間內的變化趨勢.通過與基于MSOP的時序預測方法對比，證明所提方法的預測誤差較小，預測結果較為準確.該算法完全基于基本物理規律，利用純數學推導得出，在推導過程中未加入任何約束條件或有關具體設備的輔助方程，在數學上具有本質的普適性.利用滑動數據采集窗口可以對其他時段的數據進行分析驗證，亦可通過添加其他設備的數據接口，對其他熱力參數進行分析和預測.研究成果可用于熱力系統的運行管理、控制優化和維修保障研究，達到優化運行管理和控制策略，合理安排維修保障資源、減少設備故障和延長系統使用壽命的目的.