劉乙奇 黃志鵬 于廣平 黃道平
(1.華南理工大學 自動化科學與工程學院,廣東 廣州 510640;2.廣州工業智能研究院,廣東 廣州 511458)
隨著城市化水平和工業化水平的不斷提高,我國城市污水中工業污水的含量顯著增加[1]。活性污泥法作為一種成本較為低廉、處理效果較佳的污水處理工藝,已經成為城市處理污水的主要工藝之一[2]。據不完全統計,活性污泥法是目前世界上90%的污水處理廠采用的處理方法[3]。活性污泥法的穩定運行離不開正常的泥水分離,然而,由絲狀菌過度增殖引起的污泥膨脹常常會阻礙正常的污泥分離和回流,從而使得出水水質惡化,大面積水體被污染,乃至整個污水處理系統癱瘓,污泥膨脹也因此被國內外專家一致認為是活性污泥工藝的癌癥[4]。
污泥膨脹的成因復雜,而且不同污水廠,甚至同一污水廠不同時期發生的污泥膨脹中的優勢菌種都可能不盡相同[5],因而難以建立準確的機理模型。基于多元統計分析的故障診斷方法可以避免復雜的建模過程,因此,面向污水處理過程的多元統計分析故障診斷方法近年來成為了研究的熱點[6]。它結合歷史觀測數據,用多元投影的方法將多變量樣本空間分解成主元子空間和殘差子空間,并分別在這兩個空間中構造反映空間變化的統計量,然后將在線觀測數據分別向兩個空間投影,并計算相應的統計量用于過程監控[7]。在多元統計分析方法中,典型相關分析(CCA)通過最大化兩組變量之間的相關性[8],在每組變量中適當構造若干個有代表性的綜合性指標,然后考察這些指標的相關性來反映兩組原始變量間的相關性,從而識別具有內協方差強相關的故障[9]。Chen等[10]提出了一種基于CCA的故障檢測方法,并在氧化鋁蒸發過程中進行了驗證。然而,在污泥膨脹的早期階段,整個污水處理系統的運行過程相對比較穩定,由上述方法構造的平方預測誤差(SPE,ESP)統計量所包含的故障信號不夠明顯,常常導致早期污泥膨脹檢測效果不夠穩定,甚至難以被檢測,最終導致故障的漏報問題。特征提取是解決這一問題的途徑之一。Stepanic等[11]利用時域故障特征參數和包絡分析后的故障特征頻率等統計量有效地識別了滾動軸承的故障。另外,由于基于CCA的過程監控是基于數據之間的相關性而非因果關系,所以它無法準確定位故障源[12]。其中一個解決辦法是利用貢獻圖[13],盡管貢獻圖有其優點,但非故障變量的貢獻值常常由于受到故障變量的影響而偏大,從而引起誤診。為此,Alcala等[14]提出了重構貢獻圖作為替代方案。此外,貢獻圖只能找出故障源,無法展示故障的傳播路徑。多元格蘭杰因果關系(MVGC)能夠利用觀測變量的時間序列量化各個觀測變量之間的因果關系,這可以作為分析故障傳播路徑的依據。胡瑾秋等[15]將格蘭杰因果關系檢驗應用于煉化系統故障診斷中。但如果不明顯的故障特征被送入MVGC分析,會導致分析結果有較大的偏差,甚至得出錯誤的結果。總的來說,每種方法都或多或少存在一定的局限,通過改進并整合它們,以此來實現優劣互補,是一個行之有效的解決方案。
基于上述討論,本文提出了一種新型的全生命周期故障診斷方法。首先利用離線數據構造出CCA模型,再結合基于絕對平均振幅值(AMAV)的特征提取導出新統計量SPEAMAV(ESP,AMAV)及其控制限,用于故障檢測。檢測出故障后立刻進入故障分離。為了克服貢獻圖的拖尾效應,發生故障后,結合污泥膨脹的特點,根據污泥體積指數(SVI,ISV)將觀測到的歷史數據重新排列,并重新建立CCA模型,用新的CCA模型導出的貢獻圖進行故障分離,貢獻值較大的變量被初步確定為可能的故障變量。為了進一步確定故障源,用AMAV處理每個可能的故障變量后送入MVGC分析,根據各變量間的因果關系值得出確切的故障變量。最后,為了得到整個系統的故障傳播路徑,用MVGC分析全部變量,并以故障變量為起點,得出完整的故障傳播路徑。
本文所提出新型全生命周期故障診斷方法的流程圖如圖1所示,具體說明如下:

圖1 新型全生命周期故障診斷方法的流程圖Fig.1 Flow chart of the new full life cycle fault diagnosis method
(1)故障檢測。先選定觀測變量和樣本用于訓練CCA模型,得出訓練集的SPE序列并由AMAV處理,計算出SPEAMAV控制限。隨后將來自測試集的SPE序列和訓練集的SPE序列合并,經AMAV處理后取出測試集部分的SPEAMAV用于故障檢測。一旦識別出早期故障,立即進入故障分離;否則,如果新樣本到達,更新測試集,重新返回故障檢測。
(2)故障分離。發現早期故障后,根據ISV值從小到大的順序重新排列歷史觀測樣本,然后重新建立CCA模型以提高模型精度。計算出每個變量的貢獻值,貢獻值較大的幾個變量是可能的故障變量,通過AMAV處理后送入MVGC分析,得出它們之間的因果關系值,確定出故障源。
(3)故障診斷。以故障變量作為故障路徑的起點,根據MVGC分析的結果,尋找因果關系值最大的方向,畫出故障傳播路徑。
在描述信號中的穩定分量時常用到絕對平均振幅值,它通常可以檢測信號中能量的波動[16]。絕對平均振幅值是時間序列絕對值的算術平均值,假設y(t)為時間序列,則時間區間[1,N]內y(t)的AMAV計算式為
(1)

經AMAV特征提取后的時間序列為
(2)
首先,在t=1時刻只能得到當前時刻的故障信息,無法與過去的信息進行比較,因此,對于時間序列的第一個值y(1),特征提取后僅僅相當于取其絕對值。在本文中用到的時間序列均為正值,所以特征提取后相當于保留了y(t)的原值。當t≥2時,依次計算出y(t)在對應時間段內的AMAV值,這樣可以很好地利用了之前的時間序列中所包含的故障信息。每一個元素都積累了之前所有元素所包含的故障信息,因此得出的新時間序列比原始的時間序列包含更為明顯的故障特征。在所提方法中,AMAV主要有兩個功能:①結合CCA得出新的統計量SPEAMAV及其控制限;②在進行MVGC分析之前處理原始觀測變量,從而有效地細化故障特征。
污泥膨脹是一類典型的漂移、微小故障,早期故障信號十分微弱,難以被檢測到。更為特殊的是,由于實際污水處理系統中的微生物存在自我調節的作用以及各種外部擾動的影響等原因,導致在污泥膨脹的初期階段,ISV值難以被及時、正確地檢測并真實反映出活性污泥的狀況,僅僅根據ISV值難以精準檢測出早期污泥膨脹的微小故障信號。為此,本文提出了基于AMAV-CCA的故障檢測方法。
基于CCA的故障檢測方法可以認為是基于主成分分析(PCA)和基于偏最小二乘(PLS)的故障檢測方法的有效拓展[17]。然而,由CCA導出的SPE統計量對微弱故障不夠敏感,常常導致早期故障被忽略。為了克服這一局限,本文結合基于AMAV的特征提取和CCA,推導出新的統計量SPEAMAV用于故障檢測。
假設原始數據集為X∈Rn×m,其中每行表示一個觀測樣本,每列表示一個觀測變量,X(1:n1,:)作為訓練集,X(n1+1:n,:)作為測試集。將X標準化后劃分為輸入變量X1∈Rn×m1和輸出變量X2∈Rn×m2,其中m1和m2分別為X1和X2所含的變量個數且m=m1+m2。
用訓練集估計相應的協方差矩陣:
(3)
(4)
(5)
ΣT=UΛVT
(6)
輸入空間和輸出空間的投影矩陣分別定義為
(7)
(8)
其中,k為主元個數,且k≤min{m1,m2}。
計算訓練集的殘差矩陣
etrain=X2(1:n1,:)LS-X1(1:n1,:)JSΛk
(9)
其中,Λk=diag(r1,r2,…,rk),r1,r2,…,rk為ΣT的奇異值且r1≥r2≥…≥rk。
計算訓練集的平方預測誤差時間序列
(10)
通過式(2)對ESP,train進行特征提取,得出ESP,train,AMAV,則SPEAMAV控制限δ2計算式如下:

(11)

至此,CCA模型已建立完畢。對于測試集X(n1+1:n,:),同樣計算出殘差和平方預測誤差時間序列
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時刻發生故障。由于AMAV描述的是信號中穩定分量的大小,而ESP,train作為由正常樣本推導出的統計量,含有的穩定分量較多,因此在故障檢測中將其利用起來能提高故障檢測的準確率。
當ESP,AMAV檢測到早期故障時,為了更好地決策支持,需要找出故障源。由于CCA模型是基于數據之間的相關性而不是用因果關系去定義具體的數據關系,所以故障往往會被掩蓋在相關關系中。為此,文中引入了貢獻圖用于故障分離。但傳統的貢獻圖存在故障拖尾問題。

(14)
然后通過U′(:,1:k)和V′(:,1:k)將e′分別投影至輸入和輸出空間,得到反映每個原始觀測變量的殘差矩陣
e=[e′(U′(:,1:k))T,e′(V′(:,1:k))T]
(15)
其中,e的第i列即表示第i個變量的殘差時間序列。
第i個變量的貢獻值計算式為
(16)
其中,ξi表示單位矩陣的第i列。
重排歷史觀測樣本有兩個目的:①保證訓練集的樣本都是最接近正常運行狀態的樣本,以提高CCA模型的精確度;②使測試集中的故障樣本是最接近真實故障狀態的樣本,計算出的殘差向量中包含更多的故障信息,以此得出的貢獻值更為準確,貢獻值最大或者較大的變量被認為是可能的故障變量。
在初步確定可能的故障變量后,可以利用MVGC分析得出它們之間的因果關系值,據此確定出最終的故障源。在MVGC分析之前,利用AMAV處理原始觀測變量,以獲取更為明顯的故障特征,從而提高MVGC分析的準確度。需要注意的是,由1.1節可知,經AMAV特征提取后的時間序列的第一個值本質上是原始值,故障信息從特征提取后的第二個值才開始積累。因此原始觀測變量經AMAV處理后需要將第一個值剔除,即從時間序列的第二個值開始分析。
兩個變量x1(t)、x2(t)之間的格蘭杰因果關系可以定義為:若在包含了變量x1(t)、x2(t)過去信息的條件下,對變量x2(t)的預測效果要優于只有x2(t)的過去信息的預測效果,則x1(t)是x2(t)的格蘭杰根原因[15]。為了量化x2(t)對x1(t)的格蘭杰因果關系,MVGC中使用了向量自回歸模型:
(17)
式中,p為模型階數,A11,j、A12,j、A21,j、A22,j為模型系數,ε1(t)、ε2(t)為時間序列的殘差。模型中x1(t)的回歸計算為
(18)
剔除變量x2(t)的影響后,得到
(19)
(20)
格蘭杰因果關系分析很容易推廣到多變量的情況,詳情參見MVGC工具箱[20]。
為了驗證所提方法的有效性,由某采用氧化溝法的污水處理廠提供的現場數據進行驗證。氧化溝又稱循環曝氣池,是活性污泥法的變種。世界上第一座氧化溝污水處理廠于1954年在荷蘭建立,歷經幾十年的發展,國外先后開發了多種多樣的氧化溝工藝[21]。氧化溝工藝流程圖如圖2所示。

圖2 氧化溝工藝流程圖Fig.2 Process flow diagram of oxidation ditch
該污水處理廠的主要參數如下:平均進水流量約為17萬 m3/d;平均水力停留時間(tHR)為16.5 h;生物固體停留時間(tSR)為15~22 d。所用的實驗數據均以一天為間隔進行采樣,共觀測了213組數據。觀測變量共15個,如表1所示。大約從第70天起,由于進水溫度較低,活性污泥微生物生態系統的平衡被打破,使得絲狀菌在競爭中處于優勢,絲狀菌污泥膨脹現象略有發生,持續時間約半年以上。

表1 觀測變量Table 1 Observed variables
為了評估基于AMAV-CCA模型的故障檢測性能,本文使用誤報率(FAR,RFA)、漏報率(MAR,RMA)和故障診斷準確率(A)來評估故障檢測的效果。誤報率描述了正常樣本被錯誤地識別為故障樣本的比例;漏報率描述了故障樣本被錯誤地識別為正常樣本的比例;故障診斷準確率則是兩者的綜合指標,描述了樣本被正確識別的比例[22]。即
(21)
(22)
(23)
式中,NFP為被錯誤識別的正常樣本數,NFN為被錯誤識別的故障樣本數,NTN為被正確識別的正常樣本數,NTP為被正確識別的故障樣本數。
故障檢測方法應當保證盡可能低的RFA、RMA和盡可能高的準確率A。誤報率太高會導致污水處理廠的工作人員在運行正常時頻繁地去排查各個設備的故障,從而浪費人力物力;漏報率太高則情況更為嚴重,可能導致污泥膨脹已經發生卻未被察覺,使得處理不達標的污水排入河流。
故障檢測的目的是盡早發現早期故障,從而保證有足夠的時間進行后續的分析和維護。為了對比所提方法的改善效果,分別使用傳統的CCA故障檢測方法和AMAV-CCA故障檢測方法對實驗數據進行故障檢測。首先將15個觀測變量劃分為兩個部分,前8個觀測變量組成輸入變量X1∈R213×8,后7個觀測變量組成輸出變量X2∈R213×7。對于213個數據樣本,前50個作為訓練集,用于訓練CCA模型,后163個作為測試集,用于測試統計量的性能。將數據極差標準化后,通過1.2節介紹的步驟,顯著性水平選為0.01,可得出相應的故障檢測結果,如圖3所示。從圖3(a)中可以明顯看出,由CCA導出的ESP大約能在測試集的第80天(實際的第130天)穩定地檢測到污泥膨脹,而污泥膨脹在測試集的第20天(實際的第70天)開始逐漸發生,這比實際情況晚了約60 d,不利于污水處理廠的維護。從圖3(b)中可知,ESP,AMAV明顯有更為優越的性能,由于結合了AMAV,它能在測試集的第24天(實際的第74天)檢測出故障,比SPE早了約36 d,并且在這之后不會有遺漏。此外,隨著污泥膨脹程度的逐漸加深,ESP,AMAV的值在明顯增大,說明ESP,AMAV的大小能指示出污泥膨脹的嚴重程度,這能作為污水處理廠維護的一個參考依據。另外,在測試集的第36、40、42、63、148天左右,ESP的值突然增大,而ESP,AMAV的值則發生了較為明顯的跳變,這表明該污水處理廠的污泥膨脹大約在這幾天發生了質變。SPE能檢測到故障質變的時間,卻不能跟蹤故障的劣化情況,SPEAMAV能做到這一點,因此有更為穩定的檢測效果。SPE和SPEAMAV的性能評價指標如表2所示。

圖3 兩種方法的故障檢測結果Fig.3 Fault detection results of two methods

表2 兩種統計量的性能評價指標Table 2 Performance evaluation metrics of two statistics
傳統的貢獻圖具有“故障拖尾效應”,即由于故障變量的影響,會導致非故障變量的貢獻值變大,從而導致錯誤的結果[23]。為了解決這一問題,結合污泥膨脹的特點,本文通過根據ISV值的大小重排歷史觀測樣本的方式來提高診斷精度。從2.3節可知,SPEAMAV最早能在第74天檢測出污泥膨脹,為此,利用前74 d的歷史觀測數據集,根據第1.3節介紹的步驟,得出的貢獻圖如圖4所示。

圖4 重排樣本后得到的貢獻圖Fig.4 Contribution plots after sample rearrangement
由圖4可見,變量2(即θ)的貢獻值最高,變量10(即BOD5,o)和變量15(即ISV)的貢獻值也略高。由此可以初步確定θ是故障源,但考慮到貢獻圖的拖尾效應,再用前74 d的觀測數據進行MVGC分析,以確定θ、BOD5,o和ISV之間的因果關系。根據1.4節的步驟,在進行MVGC分析前先用AMAV處理每個原始變量,得到更為明顯的故障特征,結果如表3所示。

表3 θ、BOD5,o和ISV的因果關系數值表Table 3 Causal relationship numerical table of θ,BOD5,o and ISV
由表3可見,Fθ→ISV>FISV→θ,Fθ→BOD5,o>FBOD5,o→θ,再結合貢獻圖中θ的貢獻值最大,可以確定故障變量是θ。ISV常常作為污泥膨脹的檢測性指標,可以認為其他變量都是通過直接或間接地影響ISV的大小來引發污泥膨脹。因此,可將ISV作為故障傳播路徑的終點。θ、BOD5,o和ISV變量之間存在兩條故障路徑,一條是主要路徑θ→ISV,另一條是次要路徑θ→BOD5,o→ISV。
在第74天檢測出污泥膨脹并確定出故障變量θ后,為了更好地支持決策,還需要分析出故障傳播路徑。如果用前74 d的歷史觀測數據做MVGC分析,得出來的故障傳播路徑不能完全收斂至ISV,這是因為前74 d所蘊含的故障信息只夠確定出故障變量,不足以得到完整的故障路徑。為了解決這一問題,可通過適當收集并分析檢測出早期故障以后的數據,從而積累足夠的故障信息。根據以上分析,文中用MVGC分析了前85 d的樣本,相應的因果關系數值表如表4所示。按照表中因果關系值最大的路徑,首先確定故障變量θ作為故障傳播路徑的起點,ISV作為故障傳播路徑的終點,故障傳播路徑如圖5所示,其中實線表示主要故障路徑,虛線表示次要故障路徑。圖5表明,θ作為故障變量,與ISV之間有著最為緊密的聯系,導致故障的概率也最高,能通過影響活性污泥中細菌的生長狀況來直接影響ISV,或者通過影響其他變量來間接影響ISV,如ρMLSS。溫度的變化可以導致部分細菌硝化率降低甚至死亡,從而直接影響了ρMLSS的變化,最終又間接影響二沉池的沉降率和ISV。此外,根據實際經驗,進水的水質成分、污泥負荷等,而與之相關的BOD5,i、CODi、ρTKN,i、ρTP,i等變量中,只有ρTP,i與ISV有一定的關聯,且關聯性不如θ與ISV的關聯性高,因此,通過所提方法可得,θ是引起污泥膨脹的最主要因素,而ρTP,i是次要因素。

圖5 故障傳播路徑Fig.5 Fault propagation route

表4 所有變量的因果關系數值表Table 4 Causal relationship numerical table of all variables
從以上分析可以看出,文中所提方法能夠在第74天檢測出早期污泥膨脹并確定故障源θ,在第85天確定出故障傳播路徑,有助于污水處理廠及時對污泥膨脹進行維護。
本文提出了一種新型的全生命周期故障診斷方法,涵蓋了故障檢測、故障分離、故障傳播路徑分析,并將其應用于絲狀菌污泥膨脹檢測。在該方法中,基于AMAV的特征提取充分改善了CCA對早期故障不敏感的問題,新統計量SPEAMAV不僅有更好的檢測效果,而且還能作為判斷污泥膨脹嚴重程度的一個參考。此外,貢獻圖可充分彌補CCA無法實現故障分離的局限。通過重排歷史觀測數據集的方式,在一定程度上克服了貢獻圖的拖尾效應。當確定了故障發生的時間和故障變量后,利用AMAV處理原始變量,可以提取更明顯的故障特征,再通過MVGC分析得出因果矩陣,由此得出故障傳播路徑。現場數據的實驗結果表明,文中所提方法在對污泥膨脹的系統性檢測和分析中有著明顯的優勢。可以不失一般性地認為,該方法能拓展到其他漂移、累積型故障的診斷,如材料腐蝕、齒輪劣化等。然而,為了能及時、有效地控制住污泥膨脹,在不同的情況下需要采取不同的措施。今后將研究如何根據故障診斷的結論制定出正確的維護策略。