基于先验信息的供水管网阻力系数识别

任刚红,杜坤, 和丽蓉,徐冰峰,杜雨

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

管网水力模型被越来越多的水厂用于优化供水调度、指导运营管理，如何使水力模型比较准确的反映管网实际运行状态，保证决策结果的可靠性，是目前许多水厂面临的难题。在管网水力模型中，相对于管道长度、管径等参数，管道阻力系数具有较大不确定性，需要根据实测的节点水压及管道流量进行识别，以保证管网水力模型精度。

与传统水力平差计算正好相反，管网阻力系数识别以监测的节点水压及管道流量作为已知量反算模型中管道阻力系数，学者们通常将该反问题转换为优化问题进行求解。袁一星等[1]提出了CGA-DFP混合优化算法进行管网阻力系数识别。王卓然[2]将SCEM-UA算法与EPANET水力计算模块相结合识别管网阻力系数。詹书俊等[3]通过建立以模型计算值与监测值差的多目标优化问题，采用NSGA算法求解优化问题实现管网阻力系数识别；Dini等[4]提出了基于蚁群算法的管网参数识别方法；刘永鑫等[5]利用遗传算法求解管网连续性及能量方程识别管网阻力系数；信昆仑等[6]运用全局灵敏度法进行管道摩阻灵敏度分析，采用NSGA-II算法对灵敏度较大的管道进行参数识别。

由于管网中监测点数量远少于管道数(即已知量个数少于未知量个数)，管网阻力系数识别是欠定的优化问题。针对该问题的处理，Kang等[7]根据管道的管材及管龄对管道分组，并假设同组管道阻力系数相等，将欠定问题转换为超定问题求解。Wu等[8]根据管道在管网中的位置对管道分组减少未知量个数，使管道阻力系数识别结果唯一、可靠。Mallick等[9]从识别结果稳健性角度出发探讨了管道分组问题，结果表明，管道分组这一参数化方法能有效降低识别结果方差，但由于实际管网拓扑结构不同，及监测点数量、布置差异，不存在唯一准则适用于所有管网。

本文提出基于先验信息的供水管网阻力系数识别算法，基本思路是根据管道的管材、管龄等先验信息，对管道阻力系数进行估计，并将估计值作为观测值引入目标函数，采用加权最小二乘法求解优化问题，识别管道阻力系数。与现有方法相比，所提出算法无需对管道分组，通过利用先验信息将欠定问题转化为超定，克服了现有方法中管道分组不唯一的缺点；再者，通过权重系数权衡先验信息与测量信息，避免了参数过拟合问题。此外，推导供水管网雅克比矩阵解析式用于构造搜索向量，提高了管道阻力系数识别计算效率。

1 基于先验信息的供水管网阻力系数识别框架

实际工程中管道阻力系数除了与管道的管材及管龄相关外，还与管道内壁涂料厚度、腐蚀程度及管网水力状态等随机因素相关。所提出识别算法一方面承认基于先验信息的管道阻力系数包含一定有用信息，另一方面要求参数识别结果应尽量减少模型计算值与实测值间的差异。基于此，优化问题的目标函数可构建为

(1)

式中：nH为水压监测点数；n为管网中节点数；mq为管道流量监测点数，m为管道数；Ck为管道阻力系数经验值；wH、wq、wC分别为节点水压、管道流量及管道阻力系数权重系数(分别为水压与流量监测值误差方差、管道阻力系数估计值方差的倒数)；Ho为水压监测值、H(C)为对应的模型计算值；qo为管道流量监测值、q(C)为对应的模型计算值。优化问题的约束条件为供水管网质量与能量守恒方程。为便于推导，将式(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的行向量个数大于列向量个数，管道阻力系数修正值为

(6)

图1 供水管网阻力系数识别流程图Fig.1 Water supply pipe network resistance coefficient identification flow char

2 供水管网雅克比矩阵计算

如图1所示，在管网阻力系数识别过程中需要反复计算雅克比矩阵用于构造搜索向量，但目前大多数学者[10-13]采用有限差分法估算管网雅克比矩阵，需要逐个扰动参数反复进行管网水力平差计算，故计算量巨大，不利于大型管网阻力系数识别。鉴于此，本文推导了供水管网雅克比矩阵解析式，以提高参数识别计算效率。管网质量与能量守恒方程为

(7)

式中：A为管网衔接矩阵，q为管道流量向量；Q为节点流量向量；H为节点水压向量；h为管道水头损失向量。式(7)的微分式为

(8)

(9)

(10)

根据式(10)还可得

(11)

(12)

根据式(10)、式(11)及式(12)，管道水头损失对管道阻力系数的向量微分方程可写为

Δh=ΔB-1SΔC

(13)

其中

当管网中存在水泵时，矩阵B中对应元素为(cb)-1|q|1-c，设水泵方程为hpump=a-bqc，a、b、c为水泵性能参数。根据式(7)，可得

Δh=-ATΔH

(14)

将式(14)带入式(13)，可得

BATΔH=SΔC

(15)

根据式(15)，可得

ABATΔH=ASΔC

(16)

根据式(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)

根据式(17)、式(20)，节点水压及管道流量对管道阻力系数的雅克比矩阵的解析式为

(21)

3 基于数值仿真的监测值生成

供水管网参数识别存在补偿误差问题，例如：当调整管道阻力系数或节点流量都能使模型计算值与监测值相符时，则无法分辨模型误差源于节点流量或管道阻力系数，即二者间存在补偿误差。如何获得有用监测值是在利用优化算法进行参数识别前首先应回答的问题。针对管网阻力系数识别，Ostfeld等[14]研究表明，应通过消火栓放水并记录放水量以减小节点流量补偿误差、获得有用监测值。Ormsbee等[15]指出消火栓放水至少应保证管网供水压力下降大于3.5 m，使管网处于高负荷水力“紧绷”状态，以加大监测值对管道阻力系数敏感度，否则，收集的监测值无用。

从算法验证角度来说，通过开展实地消火栓放水获得监测值，需要投入大量人力、财力，且影响管网正常运行，代价过高。再者，由于实际管网阻力系数未知，而监测值又存在误差，导致参数识别结果及模型准确性都失去参照，不利于算法验证。鉴于上述原因，为便于算法验证，本文参考文献[11]所采用的数值仿真法产生监测值。

1)管网水力模型构建。EPANET是目前使用最广泛的管网水力计算引擎，故在EPANET中构建管网水力模型开展相关研究。

2)管道阻力系数“真值”生成。考虑到实际管道阻力系数不仅与管龄及管材相关，还与管道内壁涂料厚度、腐蚀程度，及管网水力状态等随机因素相关。为准确反映实际情况，采用随机抽样法生成管道阻力系数；一般情况下，假定管道阻力系数真值服从N(Co,σ2)的正态分布，其中Co为根据先验信息估计的管道阻力系数。

3)监测值生成。应用EPANET进行管网水力计算，采用随机抽样法产生随机误差添加到计算的管道流量及节点水压中作为“真实”监测值。这里添加的随机误差可包括水压、流量本身的监测误差，同时还可包括节点流量的补偿误差。

4)管道阻力系数识别及结果评判。根据产生的“真实”监测值，应用所提出算法识别管网中各管道阻力系数。在评判识别结果时，一方面可将识别结果与管道阻力系数“真值”进行比较，另一方面可观察模型计算精度的改善情况。

4 案例分析

4.1 案例1

图2 举例管网1Fig.2 Example pipe network

表1给出了节点水压对管道阻力系数的雅克比矩阵(ABAT)-1(AS)，表2给出了管道流量对管道阻力系数的雅克比矩阵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)

根据表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分别为监测误差方差及管道阻力系数估计值方差的倒数。一般情况下，认为监测值误差服从正态分布，其中,水压监测值均方差σH=0.3 m，流量监测值均方差σq=2 L/s，则权重矩阵W为

第一次迭代修正计算值为

表4给出了迭代过程中ΔC值及管道阻力系数。

最终识别结果，表5给出了参数识别前后模型误差。

表4 代过程中ΔC值及最终识别结果Table 4 ΔC value in the process and final identification results

表5 参数识别前后模型计算值与真实值间的差Table 5 The difference between the value of the model and the real value before and after the parameter identification

续表5

由表4可知，管道1的阻力系数被较准确识别，其它管道阻力系数值基本不变。这是由于管道1为高位水池出水管，其管道流量远大于其它管道，处于高负荷水力状态，导致监测值对管道1阻力系数敏感度远大于其它管道。上述结论可通过表1、表2的雅克比矩阵进行说明。雅克比矩阵又称灵敏度矩阵，反映了监测值对参数的灵敏程度，其中元素值越大表明对应参数对模型计算精度影响越大且越容易被识别，反之亦然。例如：表1中的第一列代表了各节点水压对管道1的阻力系数灵敏度，其中各值比其它各列的值均大了一个数量级以上，表明监测值对管道1阻力系数灵敏度远大于其它管道，即管网水力模型精度主要取决于管道1的阻力系数，且管道1阻力系数更容易识别。

由表5可知，参数识别前，节点水压平均绝对误差为0.48 m，管道流量平均绝对误差为0.91 L/s；参数识别后，节点水压平均绝对误差为0.1 m，管道流量平均绝对误差为0.41 L/s，模型计算误差整体上明显减小，这表明利用所提出算法识别管网阻力系数能提高模型计算精度。此外，节点1的水压与真实值的差异仅为0.01 m，远小于监测误差值0.11 m，这表明应用所提出的加权方法能有效防止参数过拟合。再者，根据表4可知，整个参数识别过程仅需3次迭代，表明所提出算法计算效率高。

4.2 案例2

图3 举例管网2Fig.3 Example pipe network

应用所提出算法识别管网阻力系数，限于篇幅原因，不对结果进行详细列举。总体而言，与案例1类似，流量较大的主供水管道4、6、11、42、47、51的阻力系数被较准确识别，与灵敏度分析结果一致。经管道阻力系数识别，模型节点水压平均绝对误差由0.76 m降低到0.11 m，最大节点水压计算误差由1.5 m降低到0.48 m，模型计算精度有较大改善，这表明所提出算法可用于实际大型管网参数识别。此外，整个参数识别过程仅需3次迭代，且6、36节点水压计算误差小于监测值随机误差，表明所提出算法计算效率高，能有效避免参数过拟合问题。

5 结论

结果表明，所提出算法通过3次迭代就能获得最终识别结果，计算效率高且能避免参数过拟合问题。通过分析还发现，管网水力模型计算精度主要取决于管网中供水主管管道阻力系数，通过识别这些管道阻力系数，保持其他管道阻力系数不变，不失为一种可行的识别方法。值得说明的是，参数识别结果及识别后模型计算精度与监测点数量、布置位置及数据采集时管网运行状态密切相关，通常应用高负荷运行状态下的水压监测数据能得到更准确的管道阻力系数校核值。通过优化监测点布置，采集消火栓放水试验时监测值能改善识别结果、提高模型计算精度。鉴于监测点布置本身是一个复杂的优化问题，超出了本文的研究范围，在此不进行深入探讨。在工程实践中，Walski[16]建议可将监测点布置在用水量较大的节点及管网外围(远离水源)的节点。

[1] 袁一星, 张志军. 供水管网校核模型参数估计与求解方法的研究[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] 王卓然. 给水管网管道阻力系数校正试验研究[D]. 哈尔滨：哈尔滨工业大学, 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的供水管网模型校核[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] 刘永鑫, 邹平华, 马月璇. 基于遗传算法的供水管网阻力系数辨识[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] 信昆仑, 詹书俊, 陶涛,等. 基于灵敏度分析的供水管网模型多目标校核[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.