簡(jiǎn) 彩，姚 芳，趙志領(lǐng)

(1.哈爾濱工業(yè)大學(xué) 市政環(huán)境工程學(xué)院，哈爾濱 150090；2.華僑大學(xué) 土木工程學(xué)院，福建 廈門(mén) 361021)

應(yīng)用Kalman濾波算法計(jì)算供水管網(wǎng)漏失量

簡(jiǎn) 彩1，姚 芳1，趙志領(lǐng)2

(1.哈爾濱工業(yè)大學(xué) 市政環(huán)境工程學(xué)院，哈爾濱 150090；2.華僑大學(xué) 土木工程學(xué)院，福建 廈門(mén) 361021)

目前，國(guó)際上對(duì)供水管網(wǎng)漏失量化理論研究較少，而常用的計(jì)算漏失量的方法也均存在一定不足.為此,基于信息控制領(lǐng)域常用的Kalman濾波分析算法，構(gòu)建管網(wǎng)系統(tǒng)漏失模型的狀態(tài)方程與觀測(cè)方程，并遞推計(jì)算出管網(wǎng)的漏失量.進(jìn)一步在實(shí)驗(yàn)室進(jìn)行漏失模擬實(shí)驗(yàn)，利用Kalman濾波算法建立漏失模型并分析實(shí)驗(yàn)數(shù)據(jù)，結(jié)果表明，應(yīng)用基于Kalman濾波算法的漏失模型計(jì)算供水管網(wǎng)漏失量是可行且可靠的，可為供水行業(yè)制定有針對(duì)性的供水管理及技術(shù)措施、降低漏失量、提高管網(wǎng)控漏水平提供依據(jù).

供水管網(wǎng) ；漏失量計(jì)算； Kalman濾波； 數(shù)學(xué)模型；漏失實(shí)驗(yàn)?zāi)M

城市供水管網(wǎng)漏失不僅會(huì)造成寶貴水資源的浪費(fèi)和嚴(yán)重的經(jīng)濟(jì)損失，還會(huì)引發(fā)嚴(yán)重的社會(huì)問(wèn)題，如道路、地基下沉[1]等，潛在危害極大.因此，研究供水管網(wǎng)漏失問(wèn)題并最終找到控制解決辦法是極其重要且迫切的任務(wù).管網(wǎng)漏失控制的一個(gè)重要前提是較準(zhǔn)確地計(jì)算出管網(wǎng)漏失量[2].當(dāng)前國(guó)際上常用的計(jì)算漏失量的方法包括夜間最小流量法、水平衡分析法和基于管網(wǎng)漏失量水力模型的分析法.夜間最小流量法[3]簡(jiǎn)單易行，但是該方法過(guò)分依賴主觀經(jīng)驗(yàn)而缺乏理論支撐，導(dǎo)致判斷結(jié)果的精確度與可信度不高.水平衡分析法[4]通常包括依賴統(tǒng)計(jì)數(shù)據(jù)推算而得的表觀漏失水量和免費(fèi)供水量，數(shù)據(jù)的準(zhǔn)確性不足.基于管網(wǎng)漏失量水力模型的分析法，具有較強(qiáng)的漏失辨識(shí)能力，但是模型中參數(shù)的確定多是依賴經(jīng)驗(yàn)取值，導(dǎo)致模型內(nèi)核存在一定缺陷.此外，李文博[5]將灰色關(guān)聯(lián)動(dòng)態(tài)分析理論應(yīng)用于供水管網(wǎng)的漏失量化研究中，該研究側(cè)重于管網(wǎng)漏失控制工作方面，并未提出量化漏失量的具體方法.高金良等[6]嘗試用盲源分離理論來(lái)研究管網(wǎng)漏失量，具有一定創(chuàng)新性，但也難以擺脫幅值不確定性的問(wèn)題.

本研究利用Kalman濾波算法，構(gòu)建并求解供水管網(wǎng)系統(tǒng)漏失數(shù)學(xué)模型，并在實(shí)驗(yàn)室漏失模擬平臺(tái)上進(jìn)行實(shí)驗(yàn)，對(duì)比分析模型模擬計(jì)算的漏失量和實(shí)測(cè)的真實(shí)漏失量數(shù)據(jù)，驗(yàn)證了該模型算法的可行性與可靠性，為控制漏失工作做好了準(zhǔn)備.

1 Kalman濾波算法簡(jiǎn)介

Kalman濾波[7]由匈牙利數(shù)學(xué)家卡爾曼(R.E.Kalman)于1960年正式提出，隨后被廣泛應(yīng)用于各個(gè)領(lǐng)域，這是因?yàn)樗哂衅渌麨V波無(wú)法比擬的優(yōu)勢(shì)：采用的算法具有實(shí)時(shí)遞推性，即最優(yōu)化自回歸數(shù)據(jù)處理算法[8]；在分析隨機(jī)過(guò)程估計(jì)問(wèn)題時(shí)，卡爾曼將狀態(tài)空間方法靈活地融入進(jìn)時(shí)域中，提高了解決問(wèn)題的效率.

Kalman濾波處理的對(duì)象是隨機(jī)信號(hào)[9].首先利用一定的先驗(yàn)信息，以系統(tǒng)的觀測(cè)信號(hào)為輸入，以所求信號(hào)的估計(jì)值為輸出，將基本問(wèn)題描述為一個(gè)函數(shù)，即觀測(cè)方程

x=f(s,v).

(1)

式中：x表示觀測(cè)信號(hào)向量；s表示源信號(hào)向量，也是擬估計(jì)的狀態(tài)參數(shù)；v表示系統(tǒng)過(guò)程噪聲.

然后以觀測(cè)方程為約束，以合適的估計(jì)準(zhǔn)則為目標(biāo)函數(shù)，求出待求狀態(tài)參數(shù)的最優(yōu)值[10]，這就是Kalman濾波算法的應(yīng)用思路.通常來(lái)講，我們所了解的先驗(yàn)信息，是系統(tǒng)的內(nèi)部特性和系統(tǒng)狀態(tài)量的變化規(guī)律，一般表達(dá)為系統(tǒng)的狀態(tài)方程.所以，建立合理的狀態(tài)方程和觀測(cè)方程是正確應(yīng)用Kalman濾波算法的前提.Kalman濾波算法的工作原理如圖1所示.

圖1 Kalman濾波算法工作原理

2 應(yīng)用Kalman濾波算法分析計(jì)算管網(wǎng)漏失量

在管網(wǎng)總供水量和用戶總用水量均已知時(shí)，兩者之差即為所求的管網(wǎng)漏失量.管網(wǎng)總供水量一般較容易計(jì)算，本研究利用Kalman濾波算法計(jì)算實(shí)際操作中難以準(zhǔn)確測(cè)定的用戶總用水量.供水管網(wǎng)系統(tǒng)實(shí)際上是一個(gè)連續(xù)的隨機(jī)動(dòng)態(tài)系統(tǒng)，必須先將其離散化才能使用Kalman濾波算法，本研究以穩(wěn)定的運(yùn)行狀態(tài)為判斷依據(jù)按照時(shí)間段將其分割離散.然后建立關(guān)于所求狀態(tài)參數(shù)的狀態(tài)方程與觀測(cè)方程，最后依據(jù)估計(jì)準(zhǔn)則建立關(guān)于最優(yōu)解的遞推公式，從而求出待求狀態(tài)參數(shù)的最優(yōu)估計(jì)值，同時(shí)實(shí)現(xiàn)算法的遞推運(yùn)算和自我收斂.

2.1 建立狀態(tài)方程

在供水管網(wǎng)系統(tǒng)中提取關(guān)于所求狀態(tài)參數(shù)的先驗(yàn)信息，建立狀態(tài)方程：以所求的狀態(tài)量即當(dāng)前時(shí)刻的用戶用水量估計(jì)值為因變量，以前一時(shí)刻的用戶用水量估計(jì)值為自變量，并引入系統(tǒng)的控制輸入量以及系統(tǒng)過(guò)程噪聲，建立系統(tǒng)的狀態(tài)方程

(2)

式中：X(k)表示k時(shí)刻管網(wǎng)中的用水量；U(k)表示k時(shí)刻對(duì)管網(wǎng)的控制輸入；Φ為狀態(tài)轉(zhuǎn)移矩陣，反映管網(wǎng)運(yùn)行狀態(tài)的變化情況；B表示與控制輸入相關(guān)的矩陣；W(k)表示過(guò)程噪聲.

2.2 建立觀測(cè)方程

確定Kalman濾波算法的約束條件：以與用戶用水量相關(guān)的觀測(cè)信號(hào)為輸入，以所求用戶用水量的估計(jì)值為輸出，建立供水管網(wǎng)系統(tǒng)的觀測(cè)方程

(3)

式中：Z(k)表示k時(shí)刻與用戶用水量相關(guān)的觀測(cè)信號(hào)；X(k)表示k時(shí)刻管網(wǎng)中的用水量；H為測(cè)量矩陣，表示用戶用水量與觀測(cè)信號(hào)的內(nèi)在關(guān)系；V(k)表示測(cè)量噪聲.

2.3 目標(biāo)函數(shù)及最優(yōu)解的計(jì)算

本研究擬采用管網(wǎng)用戶用水量的線性最小方差估計(jì)為估計(jì)準(zhǔn)則，也是求解最優(yōu)解的目標(biāo)函數(shù)，確定Kalman濾波算法的求解過(guò)程如下

1)在管網(wǎng)上一運(yùn)行狀態(tài)K-1的基礎(chǔ)上，結(jié)合管網(wǎng)系統(tǒng)的過(guò)程模型，對(duì)管網(wǎng)下一狀態(tài)K的用水量進(jìn)行預(yù)測(cè)，計(jì)算公式為

(4)

(5)

).

(6)

4)估計(jì)準(zhǔn)則為用戶用水量滿足線性最小方差估計(jì)，卡爾曼最優(yōu)濾波增益計(jì)算公式如下

(7)

式中：Rk表示V(k)的正定方差矩陣，Hk表示k時(shí)刻的測(cè)量矩陣，HkT表示測(cè)量矩陣的轉(zhuǎn)置矩陣.

(8)

3 管網(wǎng)漏失模擬實(shí)驗(yàn)驗(yàn)證Kalman濾波算法

由于在現(xiàn)場(chǎng)難以準(zhǔn)確測(cè)定不同工況下的用戶用水量和管網(wǎng)漏失量，采用在實(shí)驗(yàn)室中進(jìn)行模擬實(shí)驗(yàn)，用不同的閥門(mén)開(kāi)啟度和水表計(jì)量的方法來(lái)模擬管網(wǎng)用戶用水量和管網(wǎng)漏失量.通過(guò)水表計(jì)量的方式進(jìn)行模擬，一部分水表模擬管網(wǎng)中用戶用水量，另一部分水表模擬管網(wǎng)中的漏失量.用人工調(diào)節(jié)水泵運(yùn)行頻率的方法模擬管網(wǎng)的不同運(yùn)行工況.在管網(wǎng)系統(tǒng)某一穩(wěn)定運(yùn)行的工況下，實(shí)驗(yàn)?zāi)M管網(wǎng)發(fā)生漏失時(shí)的情況，記錄各項(xiàng)參數(shù)，包括水泵出口流量、管網(wǎng)入口流量、管網(wǎng)總供水量、用戶總用水量、管網(wǎng)漏失量等，然后改變運(yùn)行工況與實(shí)驗(yàn)條件，記錄參數(shù)的變化情況，建立管網(wǎng)系統(tǒng)的狀態(tài)方程與參數(shù)方程，采用Kalman濾波算法并結(jié)合相關(guān)實(shí)驗(yàn)數(shù)據(jù)，計(jì)算出管網(wǎng)中用戶用水量的變化情況，最后利用差值法求出管網(wǎng)漏失量.實(shí)驗(yàn)平臺(tái)如圖2所示.

圖2 管網(wǎng)漏失模擬系統(tǒng)

3.1 基于Kalman濾波算法建立管網(wǎng)系統(tǒng)漏失數(shù)學(xué)模型

3.1.1 建立管網(wǎng)系統(tǒng)漏失數(shù)學(xué)模型前的相關(guān)假設(shè)

在建立準(zhǔn)確的管網(wǎng)系統(tǒng)漏失數(shù)學(xué)模型的過(guò)程中，對(duì)于管網(wǎng)系統(tǒng)過(guò)程噪聲和觀測(cè)噪聲相關(guān)統(tǒng)計(jì)特性的收集往往存在較大困難，故需作如下假設(shè)：供水管網(wǎng)系統(tǒng)中的過(guò)程噪聲序列W(k)和觀測(cè)噪聲序列V(k)均為高斯白噪聲隨機(jī)序列，且均值為零；Qk為非負(fù)定方差矩陣；Rk為正定方差矩陣；W(k)和V(k)二者相互獨(dú)立，且均與管網(wǎng)初始運(yùn)行狀態(tài)無(wú)關(guān).

3.1.2 狀態(tài)方程的建立

狀態(tài)量為每一穩(wěn)定工況下模擬的管網(wǎng)中用戶用水量及其隨時(shí)間變化的變化量，即

(9)

式中：Qk表示k時(shí)刻管網(wǎng)中用戶用水量，L/s；aQk表示管網(wǎng)中用水量隨時(shí)間的變化率，L/s2；n表示不同工況的數(shù)目.

(10)

若aQk為正值，表示管網(wǎng)中用水量增加；若aQk為負(fù)值，表示管網(wǎng)中用水量減少.在本研究中，管網(wǎng)系統(tǒng)沒(méi)有其他控制輸入量，即U(k)為零，則供水管網(wǎng)的狀態(tài)方程為

(11)

(12)

式中δkj為Kronecker-δ函數(shù).

綜上知

(13)

3.1.3 觀測(cè)方程的建立

供水管網(wǎng)系統(tǒng)的觀測(cè)信號(hào)Zk=(Q1,Q2,Q3,...,Qk)，其中Qi(i=1,2,...,k)表示每一運(yùn)行工況下的用水總量，計(jì)算公式為

Qi=Qzi-α0Hiβ0，i=1,2,3,...,n.

(14)

式中：Qzi表示i工況時(shí)的管網(wǎng)總水量(L/s)，可在實(shí)驗(yàn)中測(cè)算出；Hi表示i工況時(shí)的管網(wǎng)供水壓力(m)，可用水泵出口壓力、最不利控制點(diǎn)壓力和管網(wǎng)中間某節(jié)點(diǎn)壓力三者的平均值進(jìn)行計(jì)算；α0表示漏損系數(shù)；β0表示漏損指數(shù)；利用管網(wǎng)供水壓力最高值和最低值兩組數(shù)據(jù)，可計(jì)算出漏損系數(shù)α0和漏損指數(shù)β0,計(jì)算公式為

Qz=Qk+α0Hkβ0,k=1,2,3,...,n.

(15)

式中：Qz表示k時(shí)刻管網(wǎng)總水量(L/s),Qk表示k時(shí)刻管網(wǎng)總用水量(L/s),Hk表示k時(shí)刻管網(wǎng)供水壓力(m),α0表示漏損系數(shù),β0表示漏損指數(shù).

則管網(wǎng)系統(tǒng)的觀測(cè)方程為

(16)

3.2 漏失模型的求解

3.3 實(shí)驗(yàn)數(shù)據(jù)分析整理

以單水源單漏點(diǎn)環(huán)狀管網(wǎng)為例，進(jìn)行24種不同運(yùn)行工況的模擬，利用SCADA在線數(shù)據(jù)監(jiān)測(cè)設(shè)備記錄實(shí)驗(yàn)數(shù)據(jù)(實(shí)際應(yīng)用中也可以使用相同測(cè)量方法)，實(shí)測(cè)得到管網(wǎng)中6 426組流量與壓力觀測(cè)信號(hào)，并對(duì)其進(jìn)行了預(yù)處理，結(jié)果如圖3、4所示.

圖3 預(yù)處理后壓力與流量觀測(cè)信號(hào)

圖4 每一工況下壓力與流量變化

Fig.4 Flow and pressure change trend chart under each working condition

表1 基于Kalman濾波算法計(jì)算的漏失量與實(shí)驗(yàn)測(cè)得漏失量對(duì)比

Tab.1 Measured water leakage contrast to the estimate with Kalman filter algorithm

實(shí)驗(yàn)工況算法求得漏失量/(L·s-1)實(shí)測(cè)漏失量/(L·s-1)漏失量差值/(L·s-1)相對(duì)誤差/%11.33841.21310.12539.3621.68771.83280.14518.6032.06311.99820.06493.1542.80502.80010.00490.1851.99042.01030.01991.0061.47871.45720.02151.4571.99211.99580.00370.1981.54471.56730.02261.4791.85231.82990.02241.21102.02642.01870.00770.38112.31442.29890.01550.67122.06372.01170.05202.52131.36041.33270.02772.04140.89240.93740.04505.04151.28621.28450.00170.13161.89371.88590.00780.41171.64271.63770.00500.31181.56691.57790.01100.70191.58581.61450.02871.81201.50161.48570.01591.06211.88051.82010.06043.21222.20572.21020.00450.20231.86471.86250.00220.12241.46101.47070.00970.66

3.4 實(shí)驗(yàn)結(jié)果分析

由實(shí)驗(yàn)結(jié)果可知，基于Kalman濾波分析模型算法模擬計(jì)算出的管網(wǎng)漏失量相對(duì)于實(shí)測(cè)的真實(shí)漏失量的誤差范圍為0.12%～9.36%，采用相關(guān)系數(shù)作為結(jié)果可靠性程度的評(píng)價(jià)指標(biāo)，相關(guān)系數(shù)是衡量研究變量之間線性相關(guān)程度的量.兩個(gè)變量越滿足線性相關(guān)關(guān)系，相關(guān)系數(shù)越接近1.計(jì)算得出模擬計(jì)算的漏失量與實(shí)測(cè)的真實(shí)漏失量二者的相關(guān)系數(shù)為0.985，這表明模擬計(jì)算的漏失量與實(shí)測(cè)的漏失量具有非常好的線性相關(guān)關(guān)系，同時(shí)考慮到兩者之間的相對(duì)誤差范圍，說(shuō)明利用Kalman濾波算法模型計(jì)算管網(wǎng)系統(tǒng)漏失量是可行且可靠的.

分析圖5可知，在管網(wǎng)系統(tǒng)實(shí)驗(yàn)?zāi)M的前幾個(gè)運(yùn)行工況中，模擬計(jì)算的漏失量與實(shí)測(cè)漏失量二者的數(shù)據(jù)擬合程度較差，以前5個(gè)工況為例，模擬計(jì)算的漏失量與實(shí)測(cè)的真實(shí)漏失量相關(guān)系數(shù)僅為0.968，原因是在建立數(shù)學(xué)模型時(shí)，將管網(wǎng)中第一個(gè)觀測(cè)信號(hào)設(shè)定為管網(wǎng)的狀態(tài)隨機(jī)變量Xk的初值，但由于Kalman濾波具有實(shí)時(shí)遞推、自我收斂的特性，隨著算法的遞推進(jìn)行，狀態(tài)量會(huì)逐漸收斂到最優(yōu)值，并保持穩(wěn)定狀態(tài).以后10個(gè)工況為例，模擬計(jì)算的漏失量與實(shí)測(cè)的真實(shí)漏失量相關(guān)系數(shù)為0.995.這體現(xiàn)了Kalman濾波算法的實(shí)時(shí)遞推性和自回歸最優(yōu)化的特點(diǎn)，也是Kalman濾波算法能夠用來(lái)計(jì)算供水管網(wǎng)漏失量的關(guān)鍵因素.

圖5 Kalman濾波算法計(jì)算的漏失量與實(shí)驗(yàn)測(cè)得漏失量對(duì)比

Fig.5 Estimate water leakage with Kalman filter algorithm contrast to the measured water leakage

4 結(jié) 語(yǔ)

本研究建立了以Kalman濾波算法為基礎(chǔ)的供水管網(wǎng)系統(tǒng)漏失數(shù)學(xué)模型，充分發(fā)揮Kalman濾波算法自身的優(yōu)勢(shì)，成功模擬計(jì)算出管網(wǎng)漏失量，并采用相對(duì)誤差和相關(guān)系數(shù)兩個(gè)指標(biāo)來(lái)分析模擬計(jì)算出的數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的吻合程度，結(jié)果證明，模擬計(jì)算出的數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的吻合程度較高.這也證明了利用Kalman濾波算法研究計(jì)算管網(wǎng)漏失水量是完全可行可靠的，為下階段開(kāi)展降漏工作提供了有力的數(shù)據(jù)支撐.

[1] 董深，呂謀，盛澤斌,等.基于遺傳算法的供水管網(wǎng)反問(wèn)題流失定位[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2013, 45(2):106-128.

DONG S, LU M, SHENG Z B, et al. Inverse transient leakage location of water supply network based on genetic algorithm[J]. Journal of Harbin Institute of Technology,2013, 45(2):106-128.

[2] 伍悅濱，劉天順.基于瞬變反問(wèn)題分析的給水管網(wǎng)漏失數(shù)值模擬[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2005, 37(11):1483-1485.

WU Y B, LIU T X. Inverse transient analysis based numerical simulation of leakage in water distribution system[J]. Journal of Harbin Institute of Technology, 2005, 37(11):1483-1485.

[4] ALMANDOZ J. Leakage assessment through water distribution network simulation[J]. Journal of Water Resources Planning and Management, 2005,131(6):458-466.

[5] 李文博. 城市供水管網(wǎng)漏損率量化及控制研究[D]. 哈爾濱:哈爾濱工業(yè)大學(xué), 2010.

LI W B. Study on quantifying leakage percentage and controlling of urban water supply distributionsystem[D]. Harbin:Harbin Institute of Technology,2010.

[6] 高金良，李娟娟.區(qū)域供水管網(wǎng)盲源分離漏失量研究[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 2015， 47(6):33-37.

GAO J L, LI J J. Study on leakage of regional distribution network using blind source seperation[J]. Journal of Harbin Institute of Technology,2015， 47(6):33-37.

[7] KALMAN R E. A new approach to linear filtering and prediction problems[J]. Journal of basic Engineering, 1960,82(1):35-45.

[8] KEPPENNE C L. An ensemble recentering Kalman filter with an application to Argo temperature data assimilation into the NASA GEOS-5 coupled model[J]. Ocean Modelling, 2014,77:50-55.

[9] DI L， KAR S， ALSAADI F E, et al. Distributed Kalman filtering with quantized sensing state, signal processing[J]. IEEE Transactions on, 2015,63(19): 5180-5193 .

[10]KULIKOVA M V. Square-root adaptive wave filtering for marine vessels[C]//Control Conference (ECC). European, 2015: 3143-3148.

(編輯 劉 彤)

Applying Kalman filtering algorithm to calculate leakage of water distribution system

JIAN Cai1, YAO Fang1， ZHAO Zhiling2

(1.School of Municipal and Environmental Engineering，Harbin Institute of Technology, Harbin 150090,China； 2.College of Civil Engineering, Huaqiao University, Xiamen 361021, Fujian, China)

Considering the insufficiency of the present international study on the quantification of water pipeline leakage, and the deficiencies of the existing common methods of leakage calculation, this study combines Kalman filtering algorithm commonly used in the area of information control to build state equation and observation equation of pipeline leakage model. The simulation experiment of pipeline leakage calculation based on Kalman filtering algorithm conducted in the laboratory proved the feasibility of the model. Our study provides a powerful basis for policy-makings of water supply network management, leakage reduction, and improvement of leakage control level.

water distribution system；leakage calculation； Kalman filtering；mathematical model； leakage experimental simulation

10.11918/j.issn.0367-6234.201605097

2016-05-20

國(guó)家自然科學(xué)基金(51278148;51278206)；國(guó)家水體污染控制與治理科技重大專項(xiàng)(2014ZX07405002);廣東省教育部產(chǎn)學(xué)研結(jié)合項(xiàng)目(2011A090200040)

簡(jiǎn) 彩(1993—)，女，碩士研究生

簡(jiǎn) 彩，jiancheer@163.com

TU991

A

0367-6234(2017)08-0060-05