任剛紅,杜坤, 和麗蓉,徐冰峰,杜雨

(1.昆明理工大學 建筑工程學院， 昆明 650500；2.中建二局第三建筑工程有限公司，武漢 430022)

管網(wǎng)水力模型被越來越多的水廠用于優(yōu)化供水調(diào)度、指導運營管理，如何使水力模型比較準確的反映管網(wǎng)實際運行狀態(tài)，保證決策結果的可靠性，是目前許多水廠面臨的難題。在管網(wǎng)水力模型中，相對于管道長度、管徑等參數(shù)，管道阻力系數(shù)具有較大不確定性，需要根據(jù)實測的節(jié)點水壓及管道流量進行識別，以保證管網(wǎng)水力模型精度。

與傳統(tǒng)水力平差計算正好相反，管網(wǎng)阻力系數(shù)識別以監(jiān)測的節(jié)點水壓及管道流量作為已知量反算模型中管道阻力系數(shù)，學者們通常將該反問題轉(zhuǎn)換為優(yōu)化問題進行求解。袁一星等[1]提出了CGA-DFP混合優(yōu)化算法進行管網(wǎng)阻力系數(shù)識別。王卓然[2]將SCEM-UA算法與EPANET水力計算模塊相結合識別管網(wǎng)阻力系數(shù)。詹書俊等[3]通過建立以模型計算值與監(jiān)測值差的多目標優(yōu)化問題，采用NSGA算法求解優(yōu)化問題實現(xiàn)管網(wǎng)阻力系數(shù)識別；Dini等[4]提出了基于蟻群算法的管網(wǎng)參數(shù)識別方法；劉永鑫等[5]利用遺傳算法求解管網(wǎng)連續(xù)性及能量方程識別管網(wǎng)阻力系數(shù)；信昆侖等[6]運用全局靈敏度法進行管道摩阻靈敏度分析，采用NSGA-II算法對靈敏度較大的管道進行參數(shù)識別。

由于管網(wǎng)中監(jiān)測點數(shù)量遠少于管道數(shù)(即已知量個數(shù)少于未知量個數(shù))，管網(wǎng)阻力系數(shù)識別是欠定的優(yōu)化問題。針對該問題的處理，Kang等[7]根據(jù)管道的管材及管齡對管道分組，并假設同組管道阻力系數(shù)相等，將欠定問題轉(zhuǎn)換為超定問題求解。Wu等[8]根據(jù)管道在管網(wǎng)中的位置對管道分組減少未知量個數(shù)，使管道阻力系數(shù)識別結果唯一、可靠。Mallick等[9]從識別結果穩(wěn)健性角度出發(fā)探討了管道分組問題，結果表明，管道分組這一參數(shù)化方法能有效降低識別結果方差，但由于實際管網(wǎng)拓撲結構不同，及監(jiān)測點數(shù)量、布置差異，不存在唯一準則適用于所有管網(wǎng)。

本文提出基于先驗信息的供水管網(wǎng)阻力系數(shù)識別算法，基本思路是根據(jù)管道的管材、管齡等先驗信息，對管道阻力系數(shù)進行估計，并將估計值作為觀測值引入目標函數(shù)，采用加權最小二乘法求解優(yōu)化問題，識別管道阻力系數(shù)。與現(xiàn)有方法相比，所提出算法無需對管道分組，通過利用先驗信息將欠定問題轉(zhuǎn)化為超定，克服了現(xiàn)有方法中管道分組不唯一的缺點；再者，通過權重系數(shù)權衡先驗信息與測量信息，避免了參數(shù)過擬合問題。此外，推導供水管網(wǎng)雅克比矩陣解析式用于構造搜索向量，提高了管道阻力系數(shù)識別計算效率。

1 基于先驗信息的供水管網(wǎng)阻力系數(shù)識別框架

實際工程中管道阻力系數(shù)除了與管道的管材及管齡相關外，還與管道內(nèi)壁涂料厚度、腐蝕程度及管網(wǎng)水力狀態(tài)等隨機因素相關。所提出識別算法一方面承認基于先驗信息的管道阻力系數(shù)包含一定有用信息，另一方面要求參數(shù)識別結果應盡量減少模型計算值與實測值間的差異。基于此，優(yōu)化問題的目標函數(shù)可構建為

(1)

式中：nH為水壓監(jiān)測點數(shù)；n為管網(wǎng)中節(jié)點數(shù)；mq為管道流量監(jiān)測點數(shù)，m為管道數(shù)；Ck為管道阻力系數(shù)經(jīng)驗值；wH、wq、wC分別為節(jié)點水壓、管道流量及管道阻力系數(shù)權重系數(shù)(分別為水壓與流量監(jiān)測值誤差方差、管道阻力系數(shù)估計值方差的倒數(shù))；Ho為水壓監(jiān)測值、H(C)為對應的模型計算值；qo為管道流量監(jiān)測值、q(C)為對應的模型計算值。優(yōu)化問題的約束條件為供水管網(wǎng)質(zhì)量與能量守恒方程。為便于推導，將式(1)化為矩陣形式

(2)

f(Ck+ΔCk)=

(3)

式(3)的線性展開式為

f(Ck+ΔCk)≈

(4)

(5)

式中：JH(Ck)為nH×m矩陣;Jq(Ck)為mq×m矩陣;I為m×m的單位矩陣;[JH(Ck)Jq(Ck) I]T為(nH+mq+m)×m矩陣。因為nH+mq+m>m，即矩陣[JH(Ck)Jq(Ck)I]T的行向量個數(shù)大于列向量個數(shù)，管道阻力系數(shù)修正值為

(6)

圖1 供水管網(wǎng)阻力系數(shù)識別流程圖Fig.1 Water supply pipe network resistance coefficient identification flow char

2 供水管網(wǎng)雅克比矩陣計算

如圖1所示，在管網(wǎng)阻力系數(shù)識別過程中需要反復計算雅克比矩陣用于構造搜索向量，但目前大多數(shù)學者[10-13]采用有限差分法估算管網(wǎng)雅克比矩陣，需要逐個擾動參數(shù)反復進行管網(wǎng)水力平差計算，故計算量巨大，不利于大型管網(wǎng)阻力系數(shù)識別。鑒于此，本文推導了供水管網(wǎng)雅克比矩陣解析式，以提高參數(shù)識別計算效率。管網(wǎng)質(zhì)量與能量守恒方程為

(7)

式中：A為管網(wǎng)銜接矩陣，q為管道流量向量；Q為節(jié)點流量向量；H為節(jié)點水壓向量；h為管道水頭損失向量。式(7)的微分式為

(8)

(9)

(10)

根據(jù)式(10)還可得

(11)

(12)

根據(jù)式(10)、式(11)及式(12)，管道水頭損失對管道阻力系數(shù)的向量微分方程可寫為

Δh=ΔB-1SΔC

(13)

其中

當管網(wǎng)中存在水泵時，矩陣B中對應元素為(cb)-1|q|1-c，設水泵方程為hpump=a-bqc，a、b、c為水泵性能參數(shù)。根據(jù)式(7)，可得

Δh=-ATΔH

(14)

將式(14)帶入式(13)，可得

BATΔH=SΔC

(15)

根據(jù)式(15)，可得

ABATΔH=ASΔC

(16)

根據(jù)式(16)，可得

ΔH=(ABAT)-1ASΔC

(17)

同樣，管道流量的向量微分方程為

Δq=SΔC+BΔh

(18)

將式(14)帶入式(18)，可得

Δq=SΔCΔBATΔH

(19)

將式(17)帶入式(19)，可得

Δq=SΔCΔBAT(ABAT)-1ASΔC

(20)

根據(jù)式(17)、式(20)，節(jié)點水壓及管道流量對管道阻力系數(shù)的雅克比矩陣的解析式為

(21)

3 基于數(shù)值仿真的監(jiān)測值生成

供水管網(wǎng)參數(shù)識別存在補償誤差問題，例如：當調(diào)整管道阻力系數(shù)或節(jié)點流量都能使模型計算值與監(jiān)測值相符時，則無法分辨模型誤差源于節(jié)點流量或管道阻力系數(shù)，即二者間存在補償誤差。如何獲得有用監(jiān)測值是在利用優(yōu)化算法進行參數(shù)識別前首先應回答的問題。針對管網(wǎng)阻力系數(shù)識別，Ostfeld等[14]研究表明，應通過消火栓放水并記錄放水量以減小節(jié)點流量補償誤差、獲得有用監(jiān)測值。Ormsbee等[15]指出消火栓放水至少應保證管網(wǎng)供水壓力下降大于3.5 m，使管網(wǎng)處于高負荷水力“緊繃”狀態(tài)，以加大監(jiān)測值對管道阻力系數(shù)敏感度，否則，收集的監(jiān)測值無用。

從算法驗證角度來說，通過開展實地消火栓放水獲得監(jiān)測值，需要投入大量人力、財力，且影響管網(wǎng)正常運行，代價過高。再者，由于實際管網(wǎng)阻力系數(shù)未知，而監(jiān)測值又存在誤差，導致參數(shù)識別結果及模型準確性都失去參照，不利于算法驗證。鑒于上述原因，為便于算法驗證，本文參考文獻[11]所采用的數(shù)值仿真法產(chǎn)生監(jiān)測值。

1)管網(wǎng)水力模型構建。EPANET是目前使用最廣泛的管網(wǎng)水力計算引擎，故在EPANET中構建管網(wǎng)水力模型開展相關研究。

2)管道阻力系數(shù)“真值”生成?？紤]到實際管道阻力系數(shù)不僅與管齡及管材相關，還與管道內(nèi)壁涂料厚度、腐蝕程度，及管網(wǎng)水力狀態(tài)等隨機因素相關。為準確反映實際情況，采用隨機抽樣法生成管道阻力系數(shù)；一般情況下，假定管道阻力系數(shù)真值服從N(Co,σ2)的正態(tài)分布，其中Co為根據(jù)先驗信息估計的管道阻力系數(shù)。

3)監(jiān)測值生成。應用EPANET進行管網(wǎng)水力計算，采用隨機抽樣法產(chǎn)生隨機誤差添加到計算的管道流量及節(jié)點水壓中作為“真實”監(jiān)測值。這里添加的隨機誤差可包括水壓、流量本身的監(jiān)測誤差，同時還可包括節(jié)點流量的補償誤差。

4)管道阻力系數(shù)識別及結果評判。根據(jù)產(chǎn)生的“真實”監(jiān)測值，應用所提出算法識別管網(wǎng)中各管道阻力系數(shù)。在評判識別結果時，一方面可將識別結果與管道阻力系數(shù)“真值”進行比較，另一方面可觀察模型計算精度的改善情況。

4 案例分析

4.1 案例1

圖2 舉例管網(wǎng)1Fig.2 Example pipe network

表1給出了節(jié)點水壓對管道阻力系數(shù)的雅克比矩陣(ABAT)-1(AS)，表2給出了管道流量對管道阻力系數(shù)的雅克比矩陣S-BAT(ABAT)-1(AS)。

表1 雅克比矩陣(ABAT)-1(AS)Table 1 Jacobian matrix(ABAT)-1(AS)

表2 雅克比矩陣S-BAT(ABAT)-1(AS)Table 2 Jacobian matrixS-BAT(ABAT)-1(AS)

根據(jù)表1及表2給出的雅克比矩陣，所提出算法第一次迭代時的搜索向量能構建,如表3所示。

表3 第一次迭代時搜索向量[JH(C1) Jq(C1) I]TTable 3 Search for vectors for the first iteration[JH(C1) Jq(C1) I]T

其中，JH(C1)為矩陣(ABAT)-1(AS)的第一行(詳表1)；Jq(C1)為矩陣S-BAT(ABAT)-1(AS)的最后一行(詳表2)。如前述，權重矩陣W中，元素wH、wq、wC分別為監(jiān)測誤差方差及管道阻力系數(shù)估計值方差的倒數(shù)。一般情況下，認為監(jiān)測值誤差服從正態(tài)分布，其中,水壓監(jiān)測值均方差σH=0.3 m，流量監(jiān)測值均方差σq=2 L/s，則權重矩陣W為

第一次迭代修正計算值為

表4給出了迭代過程中ΔC值及管道阻力系數(shù)。

最終識別結果，表5給出了參數(shù)識別前后模型誤差。

表4 代過程中ΔC值及最終識別結果Table 4 ΔC value in the process and final identification results

表5 參數(shù)識別前后模型計算值與真實值間的差Table 5 The difference between the value of the model and the real value before and after the parameter identification

續(xù)表5

由表4可知，管道1的阻力系數(shù)被較準確識別，其它管道阻力系數(shù)值基本不變。這是由于管道1為高位水池出水管，其管道流量遠大于其它管道，處于高負荷水力狀態(tài)，導致監(jiān)測值對管道1阻力系數(shù)敏感度遠大于其它管道。上述結論可通過表1、表2的雅克比矩陣進行說明。雅克比矩陣又稱靈敏度矩陣，反映了監(jiān)測值對參數(shù)的靈敏程度，其中元素值越大表明對應參數(shù)對模型計算精度影響越大且越容易被識別，反之亦然。例如：表1中的第一列代表了各節(jié)點水壓對管道1的阻力系數(shù)靈敏度，其中各值比其它各列的值均大了一個數(shù)量級以上，表明監(jiān)測值對管道1阻力系數(shù)靈敏度遠大于其它管道，即管網(wǎng)水力模型精度主要取決于管道1的阻力系數(shù)，且管道1阻力系數(shù)更容易識別。

由表5可知，參數(shù)識別前，節(jié)點水壓平均絕對誤差為0.48 m，管道流量平均絕對誤差為0.91 L/s；參數(shù)識別后，節(jié)點水壓平均絕對誤差為0.1 m，管道流量平均絕對誤差為0.41 L/s，模型計算誤差整體上明顯減小，這表明利用所提出算法識別管網(wǎng)阻力系數(shù)能提高模型計算精度。此外，節(jié)點1的水壓與真實值的差異僅為0.01 m，遠小于監(jiān)測誤差值0.11 m，這表明應用所提出的加權方法能有效防止參數(shù)過擬合。再者，根據(jù)表4可知，整個參數(shù)識別過程僅需3次迭代，表明所提出算法計算效率高。

4.2 案例2

圖3 舉例管網(wǎng)2Fig.3 Example pipe network

應用所提出算法識別管網(wǎng)阻力系數(shù)，限于篇幅原因，不對結果進行詳細列舉。總體而言，與案例1類似，流量較大的主供水管道4、6、11、42、47、51的阻力系數(shù)被較準確識別，與靈敏度分析結果一致。經(jīng)管道阻力系數(shù)識別，模型節(jié)點水壓平均絕對誤差由0.76 m降低到0.11 m，最大節(jié)點水壓計算誤差由1.5 m降低到0.48 m，模型計算精度有較大改善，這表明所提出算法可用于實際大型管網(wǎng)參數(shù)識別。此外，整個參數(shù)識別過程僅需3次迭代，且6、36節(jié)點水壓計算誤差小于監(jiān)測值隨機誤差，表明所提出算法計算效率高，能有效避免參數(shù)過擬合問題。

5 結論

結果表明，所提出算法通過3次迭代就能獲得最終識別結果，計算效率高且能避免參數(shù)過擬合問題。通過分析還發(fā)現(xiàn)，管網(wǎng)水力模型計算精度主要取決于管網(wǎng)中供水主管管道阻力系數(shù)，通過識別這些管道阻力系數(shù)，保持其他管道阻力系數(shù)不變，不失為一種可行的識別方法。值得說明的是，參數(shù)識別結果及識別后模型計算精度與監(jiān)測點數(shù)量、布置位置及數(shù)據(jù)采集時管網(wǎng)運行狀態(tài)密切相關，通常應用高負荷運行狀態(tài)下的水壓監(jiān)測數(shù)據(jù)能得到更準確的管道阻力系數(shù)校核值。通過優(yōu)化監(jiān)測點布置，采集消火栓放水試驗時監(jiān)測值能改善識別結果、提高模型計算精度。鑒于監(jiān)測點布置本身是一個復雜的優(yōu)化問題，超出了本文的研究范圍，在此不進行深入探討。在工程實踐中，Walski[16]建議可將監(jiān)測點布置在用水量較大的節(jié)點及管網(wǎng)外圍(遠離水源)的節(jié)點。

[1] 袁一星, 張志軍. 供水管網(wǎng)校核模型參數(shù)估計與求解方法的研究[J]. 給水排水, 2005, 31(9):105-111.

YUAN Y X, ZHANG Z J. Study on estimation and approach of parameters of calibration model for water distribution network model [J]. Water & Waste Water, 2005, 31(9):105-111.(in Chinese)

[2] 王卓然. 給水管網(wǎng)管道阻力系數(shù)校正試驗研究[D]. 哈爾濱：哈爾濱工業(yè)大學, 2012.

WANG Z R. experimental study on calibration of pipeline roughness in water distribution networks[D]. Harbin：Harbin Institute of Technology, 2012.(in Chinese)

[3] 詹書俊, 陶濤. 基于NSGA—II的供水管網(wǎng)模型校核[J]. 給水排水, 2013, 39(3):158-160.

ZHAN S J, TAO R. Model checking of water supply network based on NSGA - II[J]. Water & Waste Water, 2013, 39(3):158-160.(in Chinese)

[4] DINI M, TABESH M. A new method for simultaneous calibration of demand pattern and hazen-williams coefficients in water distribution systems[J]. Water Resources Management, 2014, 28(7):2021-2034.

[5] 劉永鑫, 鄒平華, 馬月璇. 基于遺傳算法的供水管網(wǎng)阻力系數(shù)辨識[J]. 中國給水排水, 2014(23):113-116.

LIU Y X, ZOU P H, MA Y X. Identification of resistance coefficient of water supply network based on genetic algorithm[J]. China Water & Wastewater, 2014(23):113-116.(in Chinese)

[6] 信昆侖, 詹書俊, 陶濤,等. 基于靈敏度分析的供水管網(wǎng)模型多目標校核[J]. 同濟大學學報(自然科學版), 2014, 42(5):736-739.

XIN K L, ZHAN S J, TAO T，et al. Multi-objeetive Calibration of hydraulic model in water distribution network based on sensitivity analysis[J]. Journal of Tongji University(Natural Science), 2014, 42(5):736-739.(in Chinese)

[7] KANG D, LANSEY K, KANG D, et al. Demand and roughness estimation in water distribution systems [J]. Journal of Water Resources Planning & Management, 2010, 137(1):20-30.

[8] WU Z Y, CLARK C. Evolving effective hydraulic model for municipal water systems[J]. Water Resources Management, 2009, 23(1):117-136.

[9] MALLICK K N, AHMED I, TICKLE K S, et al. Determining pipe groupings for water distribution networks [J]. Journal of Water Resources Planning & Management, 2002, 128:130-139.

[10] LANSEY K E, El-SHORBAGY W, AHMED I,et al. Calibration assessment and data collection for water distributionnetworks[J].Journal of Hydraulic Engineering, 2001, 127(4): 270-279.

[11] KANG D, LANSEY K. Real-time demand estimation and confidence limit analysis for water distribution systems [J]. Journal of Hydraulic Engineering, 2009, 135(10): 825-837.

[12] PEREZ R, PUIG V, PASCUL J, et al. Pressure sensor distribution for leak detection in Barcelona water distribution network[J]. Water Science and Technology: Water Supply, 2009, 9(6): 715.

[13] MéNDEZ M, ARAYA J A, SNCHEZ L D. Automated parameter optimization of a water distribution system[J]. Journal of Hydroinformatics, 2013, 15(1):71-85.

[14] OSTFELD A, SALOMONS E, ORMSBEE L, et al. Battle of the water calibration networks [J]. Journal of Water Resources Planning & Management, 2012, 138(5):523-532.

[15] ORMSBEE L E, LINGIREDDY S. Calibration of hydraulic network models [J]. Journal AWWA,1997,89(2): 42-50.

[16] WALSKI T M. Technique for calibrating network models [J]. Journal of Water Resources Planning & Management,1983,109(4):360-372.