無跡卡爾曼濾波算法在污水處理BSM1過程中的應用研究

曾 靜 桑耀凱 李 元

(沈陽化工大學信息工程學院 遼寧 沈陽 110142)

0 引 言

污水處理是水循環的重要環節，涉及到復雜的生物和化學反應。污水處理過程(Wastewater Treatment Process,WWTP)是一個典型的非線性系統，其出水水質與環境的可持續性密切相關，通常由環境立法進行管理，在滿足環境法規的嚴格要求的同時，還要求確保工藝安全并將運行成本降至最低。污水入水流量和成分的顯著變化導致設計污水處理廠的先進控制和監測方案的復雜性增加。污水處理過程的典型仿真模型BSM1(Benchmark Simulation Model no.1)由活性污泥1號模型(ASMI)和雙指數沉淀速度模型組成，作為標準仿真平臺廣泛應用于污水過程的控制、優化、系統監視及故障診斷等研究中[1-2]。

關于污水處理過程的控制算法有許多研究成果。早期利用比例積分(Proportional-Integral，PI)控制進行污水處理廠營養物去除，PI控制可以達到運行中的穩定性要求，但不能有效處理約束[3]。近年來模型預測控制(Model Predictive Control，MPC)等先進控制方法也引起了人們的廣泛關注。英國蘭卡斯特的污水處理廠的應用證明了模型預測控制的有效性[4]。Francisco等[5]為BSM1基準仿真模型1號開發了MPC控制器，文獻[6-8]考慮了污水處理廠MPC的同步設計和控制問題。Zeng等[9]將經濟的MPC算法(Economic Model Predictive Control，EMPC)應用于BSM1，設計考慮終端代價(Terminal Cost)的經濟模型預測控制器，并將EMPC性能與常規MPC、PI進行了分析比較，在保證系統穩定性的基礎上顯著提高了出水質量，同時降低了全局能耗指標。

雖然關于污水處理廠控制系統設計的研究成果很多，但對污水處理廠狀態估計的關注相對較少。狀態估計是基于輸出測量和系統模型重構系統狀態的過程。狀態估計與污水處理廠的控制和監測密切相關，是污水處理廠運行中的一個重要課題。污水處理中的許多相關狀態是不可測量的或受到嚴重的噪聲污染。擴展卡爾曼濾波(Extended Kalman Filter，EKF)曾被應用于BSM1的狀態估計[10]。然而EKF濾波器線性化存在局限性,例如對于較強非線性模型需要計算的步長盡可能小。Julier等[11-13]提出了一種新的非線性濾波方案來替代EKF在目標跟蹤領域的應用，后來被命名為無跡卡爾曼濾波(Unscented Kalman Filter，UKF)。

基于上述考慮，本文針對污水處理的典型模型BSM1設計了無跡卡爾曼濾波器，并與擴展卡爾曼濾波的狀態估計結果進行比較。主要工作包括：(1) 詳細地描述EKF和UKF的算法理論及二者的優缺點，并描述無跡卡爾曼濾波的理論依據。(2) 將EKF與UKF在同一噪聲下應用于BSM1仿真，將仿真結果進行對比分析，得出二者方案各自的優缺點和最好的估計方案。(3) 在不同天氣條件不同噪聲大小狀況下，對EKF、UKF兩種方案在BSM1應用進行比較,分析比較結果。

1 EKF算法描述

擴展卡爾曼濾波算法理論在于對原系統方程和測量方程采取一階泰勒展開來近似逼近實際系統非線性模型。通過求偏導計算所需的雅可比矩陣確定出當前的系統矩陣F(k)和輸出矩陣H(k),為了不失一般性，我們設系統模型的狀態方程和輸出方程如下：

x(k)=f(x(k-1))+w(k-1)

(1)

z(k)=h(x(k))+v(k)

(2)

式中：f(·)為系統狀態方程、h(·)為系統輸出方程，二者為非線性函數；x為系統的狀態變量；w(k-1)為系統k-1時刻過程噪聲；v(k)為系統k時刻測量噪聲；z(k)為系統在k時刻測量值。其中噪聲均服從零均值的高斯分布。為了確定系統矩陣和輸出矩陣，對狀態方程和測量方程求偏導，得：

(3)

(4)

(5)

P(k|k-1)=F(k-1)P(k-1|k-1)FT(k-1)+Q

(6)

式中：P(k-1|k-1)為k-1時刻的誤差的協方差矩陣；P(k|k-1)為k-1時刻預測k時刻的誤差協方差矩陣；Q為過程噪聲的協方差矩陣。測量更新如下：

K(k)=P(k|k-1)HT(k)[HT(k)P(k|k-1)H(k)+R]-1

(7)

(8)

P(k|k)=P(k|k-1)-K(k)H(k)P(k|k-1)

(9)

2 UKF濾波算法描述

(1) 利用初始估計值，生成初始的L個Sigma點：

(10)

(2) 利用狀態方程計算Sigma點的一步預測：

χ(i)(k|k-1)=f(χ(i)(k-1|k-1))

(11)

(3) 對經過非線性傳遞之后Sigma點加權計算：

(12)

(4) 計算協方差矩陣的一步預測：

(13)

(5) 利用測量方程計算一步Sigma點的測量值：

Z(i)(k|k-1)=h(χ(i)(k|k-1))

(14)

(6) 計算測量值的預測加權求和：

(15)

(7) 計算新息的協方差：

(16)

(8) 計算協方差：

(17)

(9) 計算Kalman增益：

(18)

(10) 更新誤差協方差矩陣：

PX(k|k)=PX(k|k-1)-K(k)PzzK(k)T

(19)

(11) 利用新息更新狀態，z(k)為當前實測輸出：

(20)

3 仿真實驗

將上述介紹的UKF算法應用于BSM1模型，并與EKF在不同天氣和不同大小噪聲條件下的估計性能進行比較分析。BSM1的模型描述詳見文獻[10，14]。三種不同天氣條件數據文件分別為干燥天氣數據文件、降雨天的數據文件和暴風雨天的數據文件，詳見國際水協會的官網(http://www.benchmarkwwtp.org)。干燥天氣包含了兩周相同的動態干燥條件數據；雨天數據文件中包含一周的動態干燥條件數據和發生在第二周的一次長降雨事件；暴風雨天的數據文件中包含了一周的動態干燥條件和發生于第二周期間的兩次暴風雨事件。生物反應器和理想分離器的初始值通過模擬工廠100天的輸入，并在干燥條件下模擬14天計算得到的。由于三種不同天氣狀況的輸入數據前一周是完全相同的，因此仿真直接從第二周開始。在仿真中，假設對整個過程測量進行采樣是同步和周期性的，采樣周期取15 min，系統的狀態變量包括5個反應器及二次沉降池變量(SALK、XBH、XI、XS、XND、XP、XBA、SI、SS、SO、SNO、SNH、SND)見表1，反應器及沉降池可測量變量(SO，SALK，SNH，SNO，COD，CODf，BOD，TSS)見表2。在每組的仿真中，將相同的測量噪聲、過程噪聲序列和相同的初始值分別應用到EKF、UKF方案上，過程噪聲和測量噪聲均為有界高斯白噪聲。

表1 反應器中狀態變量

表2 測量輸出量

首先比較了UKF、EFK方案在干燥天氣下的性能，在本組仿真中78個狀態變量分別產生零均值，標準差為0.2x0的高斯白噪聲，并有界于[-0.2x0,0.2x0],測量傳感器自帶噪聲為零均值，標準差為0.2y0的高斯白噪聲，并有界于[-0.2y0,0.2y0]，y0=Cx0，其中：x0為初始狀態值；y0為初始的測量值；C為輸出測量矩陣。濾波方案所采用的初始估計值xinit均設置為1.02x0，根據協方差計算公式得出過程噪聲和測量噪聲的協方差矩陣分別為Q=diag((0.2x0)2)和R=diag((0.2y0)2)，初始的誤差方差矩陣設為P(0)=Q，設過程噪聲和測量噪聲以及狀態誤差之間均不相關。圖1為反應器3中實際過程狀態和UKF和EKF的濾波軌跡。其中：虛線為EKF軌跡；點線為UKF軌跡；實線為實際過程狀態軌跡。其他反應器和沉淀池具有與其類似的特點，可以看出EKF、UKF都能很好實現跟蹤，大多狀態跟蹤效果都挺好，只有個別狀態估計偏差較大如XI、XP、XBH。為了更加細致地比較EKF、UKF的性能優劣，我們需要定量計算一下二者在每個時間點的誤差均值并繪成圖像。比較方案將每組誤差統一標準化之后采用RMS(Root Mean Square)。將整個時間段[t0,ts]兩種方案每個狀態的最大估計誤差值作分母來標準化每個狀態誤差，t0、ts分別為初始時間和終止時間，狀態誤差如式(21)所示。

(21)

式中：ei(tk)為一個標準化的狀態變量的誤差序列，k為仿真時刻。計算78個狀態在時間段[t0,ts]內每一個時刻的整體誤差如式(22)所示。

(22)

采用誤差均值及最大值來比較兩種方案估計性能的優劣，計算公式如下：

(23)

Max|e|=Max(e(tk))

(24)

式中：K為[t0,ts]時間段計算的點數。根據式(21)-式(24)，通過定量計算EKF、UKF的誤差均值分別為2.94、2.39，最大標準化誤差分別為5.04、4.15，可以看出EKF給出的估計性能沒有UKF好。干燥天氣數據下UKF對比EKF性能大概提升了12.7%，圖2為二者的誤差軌跡，實線為EKF的誤差軌跡，虛線為UKF誤差軌跡，可以看到實線始終高于虛線，一開始二者比較相近，隨著輸入數據和仿真時長的增加，二者的曲線也在不斷波動變化。

圖1 干燥天氣下反應器3中EKF軌跡、UKF軌跡、 實際過程狀態軌跡

圖2 EKF誤差和UKF誤差軌跡

在另一組的仿真中，在不同天氣狀況分別對兩種估計方案的性能進行了測試。除了系統輸入數據外，所有的條件都與上一組仿真完全相同。過程噪聲取標準差為0.2x0的高斯白噪聲，測量噪聲取標準差為0.2y0的高斯白噪聲，圖3給出了雨天狀況下反應器3中的13種成分的實際狀態和兩種估計方案的軌跡。其中：虛線為EKF軌跡；點線為UKF軌跡；實線為實際過程狀態軌跡。其他反應器和沉淀池具有與反應器3類似的特點。此天氣狀況和當前噪聲環境下，EKF、UKF的估計誤差均值分別為2.6、2.26。UKF對比EKF估計精度提升了13.08%。暴風雨天狀況下，EKF、UKF的估計誤差均值分別為2.9、2.43，UKF較EKF估計精度提升16.2%。這表明了盡管惡劣天氣下輸入數據波動比較大，UKF算法仍然比EKF濾波算法的估計性能好。

圖3 雨天狀況下反應器3中EKF軌跡、UKF軌跡、 實際過程狀態軌跡

在最后一組仿真中，通過改變過程噪聲和測量噪聲的大小來測試EKF、UKF性能的優劣，初始條件xinit仍設為1.02x0，通過參數wQ和wR來調整噪聲的幅度，假設過程噪聲大小為wQx0，并有界于[-wQx0,wQx0]，則Q=diag((wQx0)2)，測量噪聲的大小為wRy0,并有界于[-wRy0,wRy0]，均為有界的高斯白噪聲，有y0=Cx0，R=diag((wRy0)2)。初始誤差方差陣P(0)=Q，不同大小過程噪聲和測量噪聲下EKF和UKF狀態估計的誤差均值如表3所示?？梢钥闯鲈诟稍锾?、雨天和暴雨天情況下，無論過程噪聲及測量噪聲大小如何改變，UKF估計性能始終優于EKF估計效果。

表3 不同噪聲、天氣條件下估計誤差均值比較

4 結 語

無跡卡爾曼濾波算法通過對非線性函數的概率密度分布近似來完成估計過程，與擴展卡爾曼算法相比不需要對Jacobin矩陣進行求導，也沒有忽略高階項，因此對非線性系統的狀態估計有較高的計算精度。本文將無跡卡爾曼濾波算法引入到污水處理過程帶的狀態估計問題中，并與擴展卡爾曼濾波算法相比較，結果表明無論天氣狀況如何，在相同模型噪聲與測量噪聲下的UKF估計精度始終高于EKF的估計精度，并且魯棒性好，UKF較EKF的估計精度提升了10%～20%左右，顯著提高了對復雜污水處理過程的狀態估計效果。