劉乙奇 黃志鵬 于廣平 黃道平
(1.華南理工大學(xué) 自動化科學(xué)與工程學(xué)院,廣東 廣州 510640;2.廣州工業(yè)智能研究院,廣東 廣州 511458)
隨著城市化水平和工業(yè)化水平的不斷提高,我國城市污水中工業(yè)污水的含量顯著增加[1]?;钚晕勰喾ㄗ鳛橐环N成本較為低廉、處理效果較佳的污水處理工藝,已經(jīng)成為城市處理污水的主要工藝之一[2]。據(jù)不完全統(tǒng)計,活性污泥法是目前世界上90%的污水處理廠采用的處理方法[3]?;钚晕勰喾ǖ姆€(wěn)定運行離不開正常的泥水分離,然而,由絲狀菌過度增殖引起的污泥膨脹常常會阻礙正常的污泥分離和回流,從而使得出水水質(zhì)惡化,大面積水體被污染,乃至整個污水處理系統(tǒng)癱瘓,污泥膨脹也因此被國內(nèi)外專家一致認(rèn)為是活性污泥工藝的癌癥[4]。
污泥膨脹的成因復(fù)雜,而且不同污水廠,甚至同一污水廠不同時期發(fā)生的污泥膨脹中的優(yōu)勢菌種都可能不盡相同[5],因而難以建立準(zhǔn)確的機理模型?;诙嘣y(tǒng)計分析的故障診斷方法可以避免復(fù)雜的建模過程,因此,面向污水處理過程的多元統(tǒng)計分析故障診斷方法近年來成為了研究的熱點[6]。它結(jié)合歷史觀測數(shù)據(jù),用多元投影的方法將多變量樣本空間分解成主元子空間和殘差子空間,并分別在這兩個空間中構(gòu)造反映空間變化的統(tǒng)計量,然后將在線觀測數(shù)據(jù)分別向兩個空間投影,并計算相應(yīng)的統(tǒng)計量用于過程監(jiān)控[7]。在多元統(tǒng)計分析方法中,典型相關(guān)分析(CCA)通過最大化兩組變量之間的相關(guān)性[8],在每組變量中適當(dāng)構(gòu)造若干個有代表性的綜合性指標(biāo),然后考察這些指標(biāo)的相關(guān)性來反映兩組原始變量間的相關(guān)性,從而識別具有內(nèi)協(xié)方差強相關(guān)的故障[9]。Chen等[10]提出了一種基于CCA的故障檢測方法,并在氧化鋁蒸發(fā)過程中進(jìn)行了驗證。然而,在污泥膨脹的早期階段,整個污水處理系統(tǒng)的運行過程相對比較穩(wěn)定,由上述方法構(gòu)造的平方預(yù)測誤差(SPE,ESP)統(tǒng)計量所包含的故障信號不夠明顯,常常導(dǎo)致早期污泥膨脹檢測效果不夠穩(wěn)定,甚至難以被檢測,最終導(dǎo)致故障的漏報問題。特征提取是解決這一問題的途徑之一。Stepanic等[11]利用時域故障特征參數(shù)和包絡(luò)分析后的故障特征頻率等統(tǒng)計量有效地識別了滾動軸承的故障。另外,由于基于CCA的過程監(jiān)控是基于數(shù)據(jù)之間的相關(guān)性而非因果關(guān)系,所以它無法準(zhǔn)確定位故障源[12]。其中一個解決辦法是利用貢獻(xiàn)圖[13],盡管貢獻(xiàn)圖有其優(yōu)點,但非故障變量的貢獻(xiàn)值常常由于受到故障變量的影響而偏大,從而引起誤診。為此,Alcala等[14]提出了重構(gòu)貢獻(xiàn)圖作為替代方案。此外,貢獻(xiàn)圖只能找出故障源,無法展示故障的傳播路徑。多元格蘭杰因果關(guān)系(MVGC)能夠利用觀測變量的時間序列量化各個觀測變量之間的因果關(guān)系,這可以作為分析故障傳播路徑的依據(jù)。胡瑾秋等[15]將格蘭杰因果關(guān)系檢驗應(yīng)用于煉化系統(tǒng)故障診斷中。但如果不明顯的故障特征被送入MVGC分析,會導(dǎo)致分析結(jié)果有較大的偏差,甚至得出錯誤的結(jié)果。總的來說,每種方法都或多或少存在一定的局限,通過改進(jìn)并整合它們,以此來實現(xiàn)優(yōu)劣互補,是一個行之有效的解決方案。
基于上述討論,本文提出了一種新型的全生命周期故障診斷方法。首先利用離線數(shù)據(jù)構(gòu)造出CCA模型,再結(jié)合基于絕對平均振幅值(AMAV)的特征提取導(dǎo)出新統(tǒng)計量SPEAMAV(ESP,AMAV)及其控制限,用于故障檢測。檢測出故障后立刻進(jìn)入故障分離。為了克服貢獻(xiàn)圖的拖尾效應(yīng),發(fā)生故障后,結(jié)合污泥膨脹的特點,根據(jù)污泥體積指數(shù)(SVI,ISV)將觀測到的歷史數(shù)據(jù)重新排列,并重新建立CCA模型,用新的CCA模型導(dǎo)出的貢獻(xiàn)圖進(jìn)行故障分離,貢獻(xiàn)值較大的變量被初步確定為可能的故障變量。為了進(jìn)一步確定故障源,用AMAV處理每個可能的故障變量后送入MVGC分析,根據(jù)各變量間的因果關(guān)系值得出確切的故障變量。最后,為了得到整個系統(tǒng)的故障傳播路徑,用MVGC分析全部變量,并以故障變量為起點,得出完整的故障傳播路徑。
本文所提出新型全生命周期故障診斷方法的流程圖如圖1所示,具體說明如下:
圖1 新型全生命周期故障診斷方法的流程圖Fig.1 Flow chart of the new full life cycle fault diagnosis method
(1)故障檢測。先選定觀測變量和樣本用于訓(xùn)練CCA模型,得出訓(xùn)練集的SPE序列并由AMAV處理,計算出SPEAMAV控制限。隨后將來自測試集的SPE序列和訓(xùn)練集的SPE序列合并,經(jīng)AMAV處理后取出測試集部分的SPEAMAV用于故障檢測。一旦識別出早期故障,立即進(jìn)入故障分離;否則,如果新樣本到達(dá),更新測試集,重新返回故障檢測。
(2)故障分離。發(fā)現(xiàn)早期故障后,根據(jù)ISV值從小到大的順序重新排列歷史觀測樣本,然后重新建立CCA模型以提高模型精度。計算出每個變量的貢獻(xiàn)值,貢獻(xiàn)值較大的幾個變量是可能的故障變量,通過AMAV處理后送入MVGC分析,得出它們之間的因果關(guān)系值,確定出故障源。
(3)故障診斷。以故障變量作為故障路徑的起點,根據(jù)MVGC分析的結(jié)果,尋找因果關(guān)系值最大的方向,畫出故障傳播路徑。
在描述信號中的穩(wěn)定分量時常用到絕對平均振幅值,它通??梢詸z測信號中能量的波動[16]。絕對平均振幅值是時間序列絕對值的算術(shù)平均值,假設(shè)y(t)為時間序列,則時間區(qū)間[1,N]內(nèi)y(t)的AMAV計算式為
(1)
經(jīng)AMAV特征提取后的時間序列為
(2)
首先,在t=1時刻只能得到當(dāng)前時刻的故障信息,無法與過去的信息進(jìn)行比較,因此,對于時間序列的第一個值y(1),特征提取后僅僅相當(dāng)于取其絕對值。在本文中用到的時間序列均為正值,所以特征提取后相當(dāng)于保留了y(t)的原值。當(dāng)t≥2時,依次計算出y(t)在對應(yīng)時間段內(nèi)的AMAV值,這樣可以很好地利用了之前的時間序列中所包含的故障信息。每一個元素都積累了之前所有元素所包含的故障信息,因此得出的新時間序列比原始的時間序列包含更為明顯的故障特征。在所提方法中,AMAV主要有兩個功能:①結(jié)合CCA得出新的統(tǒng)計量SPEAMAV及其控制限;②在進(jìn)行MVGC分析之前處理原始觀測變量,從而有效地細(xì)化故障特征。
污泥膨脹是一類典型的漂移、微小故障,早期故障信號十分微弱,難以被檢測到。更為特殊的是,由于實際污水處理系統(tǒng)中的微生物存在自我調(diào)節(jié)的作用以及各種外部擾動的影響等原因,導(dǎo)致在污泥膨脹的初期階段,ISV值難以被及時、正確地檢測并真實反映出活性污泥的狀況,僅僅根據(jù)ISV值難以精準(zhǔn)檢測出早期污泥膨脹的微小故障信號。為此,本文提出了基于AMAV-CCA的故障檢測方法。
基于CCA的故障檢測方法可以認(rèn)為是基于主成分分析(PCA)和基于偏最小二乘(PLS)的故障檢測方法的有效拓展[17]。然而,由CCA導(dǎo)出的SPE統(tǒng)計量對微弱故障不夠敏感,常常導(dǎo)致早期故障被忽略。為了克服這一局限,本文結(jié)合基于AMAV的特征提取和CCA,推導(dǎo)出新的統(tǒng)計量SPEAMAV用于故障檢測。
假設(shè)原始數(shù)據(jù)集為X∈Rn×m,其中每行表示一個觀測樣本,每列表示一個觀測變量,X(1:n1,:)作為訓(xùn)練集,X(n1+1:n,:)作為測試集。將X標(biāo)準(zhǔn)化后劃分為輸入變量X1∈Rn×m1和輸出變量X2∈Rn×m2,其中m1和m2分別為X1和X2所含的變量個數(shù)且m=m1+m2。
用訓(xùn)練集估計相應(yīng)的協(xié)方差矩陣:
(3)
(4)
(5)
ΣT=UΛVT
(6)
輸入空間和輸出空間的投影矩陣分別定義為
(7)
(8)
其中,k為主元個數(shù),且k≤min{m1,m2}。
計算訓(xùn)練集的殘差矩陣
etrain=X2(1:n1,:)LS-X1(1:n1,:)JSΛk
(9)
其中,Λk=diag(r1,r2,…,rk),r1,r2,…,rk為ΣT的奇異值且r1≥r2≥…≥rk。
計算訓(xùn)練集的平方預(yù)測誤差時間序列
(10)
通過式(2)對ESP,train進(jìn)行特征提取,得出ESP,train,AMAV,則SPEAMAV控制限δ2計算式如下:
(11)
至此,CCA模型已建立完畢。對于測試集X(n1+1:n,:),同樣計算出殘差和平方預(yù)測誤差時間序列
etest=X2(n1+1:n,:)LS-X1(n1+1:n,:)JSΛk
(12)
(13)
最后將ESP,train和ESP,test合并得到[ESP,train,ESP,test],再通過式(2)計算得到[ESP,train,AMAV,ESP,test,AMAV],一旦ESP,test,AMAV(i)≥δ2,則表明在第i時刻發(fā)生故障。由于AMAV描述的是信號中穩(wěn)定分量的大小,而ESP,train作為由正常樣本推導(dǎo)出的統(tǒng)計量,含有的穩(wěn)定分量較多,因此在故障檢測中將其利用起來能提高故障檢測的準(zhǔn)確率。
當(dāng)ESP,AMAV檢測到早期故障時,為了更好地決策支持,需要找出故障源。由于CCA模型是基于數(shù)據(jù)之間的相關(guān)性而不是用因果關(guān)系去定義具體的數(shù)據(jù)關(guān)系,所以故障往往會被掩蓋在相關(guān)關(guān)系中。為此,文中引入了貢獻(xiàn)圖用于故障分離。但傳統(tǒng)的貢獻(xiàn)圖存在故障拖尾問題。
(14)
然后通過U′(:,1:k)和V′(:,1:k)將e′分別投影至輸入和輸出空間,得到反映每個原始觀測變量的殘差矩陣
e=[e′(U′(:,1:k))T,e′(V′(:,1:k))T]
(15)
其中,e的第i列即表示第i個變量的殘差時間序列。
第i個變量的貢獻(xiàn)值計算式為
(16)
其中,ξi表示單位矩陣的第i列。
重排歷史觀測樣本有兩個目的:①保證訓(xùn)練集的樣本都是最接近正常運行狀態(tài)的樣本,以提高CCA模型的精確度;②使測試集中的故障樣本是最接近真實故障狀態(tài)的樣本,計算出的殘差向量中包含更多的故障信息,以此得出的貢獻(xiàn)值更為準(zhǔn)確,貢獻(xiàn)值最大或者較大的變量被認(rèn)為是可能的故障變量。
在初步確定可能的故障變量后,可以利用MVGC分析得出它們之間的因果關(guān)系值,據(jù)此確定出最終的故障源。在MVGC分析之前,利用AMAV處理原始觀測變量,以獲取更為明顯的故障特征,從而提高M(jìn)VGC分析的準(zhǔn)確度。需要注意的是,由1.1節(jié)可知,經(jīng)AMAV特征提取后的時間序列的第一個值本質(zhì)上是原始值,故障信息從特征提取后的第二個值才開始積累。因此原始觀測變量經(jīng)AMAV處理后需要將第一個值剔除,即從時間序列的第二個值開始分析。
兩個變量x1(t)、x2(t)之間的格蘭杰因果關(guān)系可以定義為:若在包含了變量x1(t)、x2(t)過去信息的條件下,對變量x2(t)的預(yù)測效果要優(yōu)于只有x2(t)的過去信息的預(yù)測效果,則x1(t)是x2(t)的格蘭杰根原因[15]。為了量化x2(t)對x1(t)的格蘭杰因果關(guān)系,MVGC中使用了向量自回歸模型:
(17)
式中,p為模型階數(shù),A11,j、A12,j、A21,j、A22,j為模型系數(shù),ε1(t)、ε2(t)為時間序列的殘差。模型中x1(t)的回歸計算為
(18)
剔除變量x2(t)的影響后,得到
(19)
(20)
格蘭杰因果關(guān)系分析很容易推廣到多變量的情況,詳情參見MVGC工具箱[20]。
為了驗證所提方法的有效性,由某采用氧化溝法的污水處理廠提供的現(xiàn)場數(shù)據(jù)進(jìn)行驗證。氧化溝又稱循環(huán)曝氣池,是活性污泥法的變種。世界上第一座氧化溝污水處理廠于1954年在荷蘭建立,歷經(jīng)幾十年的發(fā)展,國外先后開發(fā)了多種多樣的氧化溝工藝[21]。氧化溝工藝流程圖如圖2所示。
圖2 氧化溝工藝流程圖Fig.2 Process flow diagram of oxidation ditch
該污水處理廠的主要參數(shù)如下:平均進(jìn)水流量約為17萬 m3/d;平均水力停留時間(tHR)為16.5 h;生物固體停留時間(tSR)為15~22 d。所用的實驗數(shù)據(jù)均以一天為間隔進(jìn)行采樣,共觀測了213組數(shù)據(jù)。觀測變量共15個,如表1所示。大約從第70天起,由于進(jìn)水溫度較低,活性污泥微生物生態(tài)系統(tǒng)的平衡被打破,使得絲狀菌在競爭中處于優(yōu)勢,絲狀菌污泥膨脹現(xiàn)象略有發(fā)生,持續(xù)時間約半年以上。
表1 觀測變量Table 1 Observed variables
為了評估基于AMAV-CCA模型的故障檢測性能,本文使用誤報率(FAR,RFA)、漏報率(MAR,RMA)和故障診斷準(zhǔn)確率(A)來評估故障檢測的效果。誤報率描述了正常樣本被錯誤地識別為故障樣本的比例;漏報率描述了故障樣本被錯誤地識別為正常樣本的比例;故障診斷準(zhǔn)確率則是兩者的綜合指標(biāo),描述了樣本被正確識別的比例[22]。即
(21)
(22)
(23)
式中,NFP為被錯誤識別的正常樣本數(shù),NFN為被錯誤識別的故障樣本數(shù),NTN為被正確識別的正常樣本數(shù),NTP為被正確識別的故障樣本數(shù)。
故障檢測方法應(yīng)當(dāng)保證盡可能低的RFA、RMA和盡可能高的準(zhǔn)確率A。誤報率太高會導(dǎo)致污水處理廠的工作人員在運行正常時頻繁地去排查各個設(shè)備的故障,從而浪費人力物力;漏報率太高則情況更為嚴(yán)重,可能導(dǎo)致污泥膨脹已經(jīng)發(fā)生卻未被察覺,使得處理不達(dá)標(biāo)的污水排入河流。
故障檢測的目的是盡早發(fā)現(xiàn)早期故障,從而保證有足夠的時間進(jìn)行后續(xù)的分析和維護(hù)。為了對比所提方法的改善效果,分別使用傳統(tǒng)的CCA故障檢測方法和AMAV-CCA故障檢測方法對實驗數(shù)據(jù)進(jìn)行故障檢測。首先將15個觀測變量劃分為兩個部分,前8個觀測變量組成輸入變量X1∈R213×8,后7個觀測變量組成輸出變量X2∈R213×7。對于213個數(shù)據(jù)樣本,前50個作為訓(xùn)練集,用于訓(xùn)練CCA模型,后163個作為測試集,用于測試統(tǒng)計量的性能。將數(shù)據(jù)極差標(biāo)準(zhǔn)化后,通過1.2節(jié)介紹的步驟,顯著性水平選為0.01,可得出相應(yīng)的故障檢測結(jié)果,如圖3所示。從圖3(a)中可以明顯看出,由CCA導(dǎo)出的ESP大約能在測試集的第80天(實際的第130天)穩(wěn)定地檢測到污泥膨脹,而污泥膨脹在測試集的第20天(實際的第70天)開始逐漸發(fā)生,這比實際情況晚了約60 d,不利于污水處理廠的維護(hù)。從圖3(b)中可知,ESP,AMAV明顯有更為優(yōu)越的性能,由于結(jié)合了AMAV,它能在測試集的第24天(實際的第74天)檢測出故障,比SPE早了約36 d,并且在這之后不會有遺漏。此外,隨著污泥膨脹程度的逐漸加深,ESP,AMAV的值在明顯增大,說明ESP,AMAV的大小能指示出污泥膨脹的嚴(yán)重程度,這能作為污水處理廠維護(hù)的一個參考依據(jù)。另外,在測試集的第36、40、42、63、148天左右,ESP的值突然增大,而ESP,AMAV的值則發(fā)生了較為明顯的跳變,這表明該污水處理廠的污泥膨脹大約在這幾天發(fā)生了質(zhì)變。SPE能檢測到故障質(zhì)變的時間,卻不能跟蹤故障的劣化情況,SPEAMAV能做到這一點,因此有更為穩(wěn)定的檢測效果。SPE和SPEAMAV的性能評價指標(biāo)如表2所示。
圖3 兩種方法的故障檢測結(jié)果Fig.3 Fault detection results of two methods
表2 兩種統(tǒng)計量的性能評價指標(biāo)Table 2 Performance evaluation metrics of two statistics
傳統(tǒng)的貢獻(xiàn)圖具有“故障拖尾效應(yīng)”,即由于故障變量的影響,會導(dǎo)致非故障變量的貢獻(xiàn)值變大,從而導(dǎo)致錯誤的結(jié)果[23]。為了解決這一問題,結(jié)合污泥膨脹的特點,本文通過根據(jù)ISV值的大小重排歷史觀測樣本的方式來提高診斷精度。從2.3節(jié)可知,SPEAMAV最早能在第74天檢測出污泥膨脹,為此,利用前74 d的歷史觀測數(shù)據(jù)集,根據(jù)第1.3節(jié)介紹的步驟,得出的貢獻(xiàn)圖如圖4所示。
圖4 重排樣本后得到的貢獻(xiàn)圖Fig.4 Contribution plots after sample rearrangement
由圖4可見,變量2(即θ)的貢獻(xiàn)值最高,變量10(即BOD5,o)和變量15(即ISV)的貢獻(xiàn)值也略高。由此可以初步確定θ是故障源,但考慮到貢獻(xiàn)圖的拖尾效應(yīng),再用前74 d的觀測數(shù)據(jù)進(jìn)行MVGC分析,以確定θ、BOD5,o和ISV之間的因果關(guān)系。根據(jù)1.4節(jié)的步驟,在進(jìn)行MVGC分析前先用AMAV處理每個原始變量,得到更為明顯的故障特征,結(jié)果如表3所示。
表3 θ、BOD5,o和ISV的因果關(guān)系數(shù)值表Table 3 Causal relationship numerical table of θ,BOD5,o and ISV
由表3可見,F(xiàn)θ→ISV>FISV→θ,F(xiàn)θ→BOD5,o>FBOD5,o→θ,再結(jié)合貢獻(xiàn)圖中θ的貢獻(xiàn)值最大,可以確定故障變量是θ。ISV常常作為污泥膨脹的檢測性指標(biāo),可以認(rèn)為其他變量都是通過直接或間接地影響ISV的大小來引發(fā)污泥膨脹。因此,可將ISV作為故障傳播路徑的終點。θ、BOD5,o和ISV變量之間存在兩條故障路徑,一條是主要路徑θ→ISV,另一條是次要路徑θ→BOD5,o→ISV。
在第74天檢測出污泥膨脹并確定出故障變量θ后,為了更好地支持決策,還需要分析出故障傳播路徑。如果用前74 d的歷史觀測數(shù)據(jù)做MVGC分析,得出來的故障傳播路徑不能完全收斂至ISV,這是因為前74 d所蘊含的故障信息只夠確定出故障變量,不足以得到完整的故障路徑。為了解決這一問題,可通過適當(dāng)收集并分析檢測出早期故障以后的數(shù)據(jù),從而積累足夠的故障信息。根據(jù)以上分析,文中用MVGC分析了前85 d的樣本,相應(yīng)的因果關(guān)系數(shù)值表如表4所示。按照表中因果關(guān)系值最大的路徑,首先確定故障變量θ作為故障傳播路徑的起點,ISV作為故障傳播路徑的終點,故障傳播路徑如圖5所示,其中實線表示主要故障路徑,虛線表示次要故障路徑。圖5表明,θ作為故障變量,與ISV之間有著最為緊密的聯(lián)系,導(dǎo)致故障的概率也最高,能通過影響活性污泥中細(xì)菌的生長狀況來直接影響ISV,或者通過影響其他變量來間接影響ISV,如ρMLSS。溫度的變化可以導(dǎo)致部分細(xì)菌硝化率降低甚至死亡,從而直接影響了ρMLSS的變化,最終又間接影響二沉池的沉降率和ISV。此外,根據(jù)實際經(jīng)驗,進(jìn)水的水質(zhì)成分、污泥負(fù)荷等,而與之相關(guān)的BOD5,i、CODi、ρTKN,i、ρTP,i等變量中,只有ρTP,i與ISV有一定的關(guān)聯(lián),且關(guān)聯(lián)性不如θ與ISV的關(guān)聯(lián)性高,因此,通過所提方法可得,θ是引起污泥膨脹的最主要因素,而ρTP,i是次要因素。
圖5 故障傳播路徑Fig.5 Fault propagation route
表4 所有變量的因果關(guān)系數(shù)值表Table 4 Causal relationship numerical table of all variables
從以上分析可以看出,文中所提方法能夠在第74天檢測出早期污泥膨脹并確定故障源θ,在第85天確定出故障傳播路徑,有助于污水處理廠及時對污泥膨脹進(jìn)行維護(hù)。
本文提出了一種新型的全生命周期故障診斷方法,涵蓋了故障檢測、故障分離、故障傳播路徑分析,并將其應(yīng)用于絲狀菌污泥膨脹檢測。在該方法中,基于AMAV的特征提取充分改善了CCA對早期故障不敏感的問題,新統(tǒng)計量SPEAMAV不僅有更好的檢測效果,而且還能作為判斷污泥膨脹嚴(yán)重程度的一個參考。此外,貢獻(xiàn)圖可充分彌補CCA無法實現(xiàn)故障分離的局限。通過重排歷史觀測數(shù)據(jù)集的方式,在一定程度上克服了貢獻(xiàn)圖的拖尾效應(yīng)。當(dāng)確定了故障發(fā)生的時間和故障變量后,利用AMAV處理原始變量,可以提取更明顯的故障特征,再通過MVGC分析得出因果矩陣,由此得出故障傳播路徑。現(xiàn)場數(shù)據(jù)的實驗結(jié)果表明,文中所提方法在對污泥膨脹的系統(tǒng)性檢測和分析中有著明顯的優(yōu)勢??梢圆皇б话阈缘卣J(rèn)為,該方法能拓展到其他漂移、累積型故障的診斷,如材料腐蝕、齒輪劣化等。然而,為了能及時、有效地控制住污泥膨脹,在不同的情況下需要采取不同的措施。今后將研究如何根據(jù)故障診斷的結(jié)論制定出正確的維護(hù)策略。