模型参数不确定性对排水管网功能状态评价的影响

吴 珊， 程玉林， 侯本伟， 李化雨， 李 俊

(1.北京工业大学建筑工程学院， 北京 100124； 2.北京市市政工程设计研究总院有限公司， 北京 100082)

排水管网水力模型广泛应用于城市排水系统的现状评估、改造规划以及新建城市排水系统的规划设计等方面[1]. 在实际建模与分析计算中，降雨荷载的时空分布[2]、模型参数的设置[3]、管道使用过程中的淤积[4]等都会导致模型参数取值存在不确定性，影响模型评价结果. 因此，有必要进行排水管网系统的不确定性分析，以探究参数不确定性对节点和管道功能状态评价的影响.

不确定性是指事件出现或发生的结果是不能准确确定的[5]，Tung等[6]认为水环境系统模型的不确定性来源是缺乏完善的信息对问题定义和求解中涉及的数据、现象和过程进行准确描述. Refsgaard等[7]将水环境系统模拟过程的不确定性分为输入数据、模型参数和模型结构的不确定性. Deletic等[8]将模型的不确定性分为输入数据、模型结构和模型校核的不确定性. 周玉文等[4]认为排水管网系统的不确定性来自于设计流量、设计过程和管网使用过程中的不确定性，并建议将不确定性的参数表示为随机变量，采用可靠度理论分析排水管网系统的风险.

在排水管网的不确定性分析方面，张子贤等[9]考虑管道直径、粗糙系数、水力坡度属性参数的不确定性，采用一次二阶可靠度法(first order reliability method，FORM)对雨水管道的设计过水能力进行了可靠度分析. 李芊等[10]以管道的管径、粗糙系数、铺设偏差3类参数为随机变量，采用Monte Carlo模拟方法，分析了雨水管网系统的可靠度. Thorndahl等[11]以暴雨深度、持续时间、峰值强度等降雨荷载特征参数为随机变量，采用FORM方法和Monte Carlo方法，分析了合流制排水管网系统的可靠度. Dotto等[12]以降雨时间步长、降雨随机空间分布、系统补偿为随机变量，采用马尔可夫链蒙特卡洛方法，分析了降雨数据不确定性对城市排水管网模型径流量计算结果的影响. Gouri等[13]考虑管渠的粗糙系数、坡度，尺寸属性参数的不确定性，采用FORM方法分析了管渠排水能力的可靠度. Tung[14]以径流系数、降雨的随机不确定性、降雨采样误差等参数为随机变量，采用概率点估计(probabilistic point estimation，PPE)方法研究了蓄水池设计的可靠性.

上述文献表明，排水管网模型在降雨输入以及产流、汇流、管网水动力计算过程中均存在不确定性的数据或参数，而已有的研究多关注其中某个过程中数据或参数不确定性的影响；对于全过程中各类不确定性参数对模型计算结果的影响关注较少. 本文建立的不确定性分析模型考虑了降雨输入和产流、汇流、管网水动力计算过程中多个不确定参数的影响，并比较了不确定性分析模型与确定性工况评价结果的差异性，分析了参数不确定性对节点和管道功能状态评价的影响.

1 模型不确定性参数

本文采用美国环境保护署开发的暴雨洪水管理模型(storm water management model，SWMM)软件进行区域排水系统水力模型的建模和计算，经过近50年的不断完善，当前版本的SWMM 5.1已经被广泛应用于环保、城市防洪、生态修复、规划设计等领域[15]. 降雨模型为SWMM荷载参数的主要输入信息，产流模型、汇流模型和管网水动力模型是其计算原理的3个主要组成部分[16].

降雨模型输入时，芝加哥雨型由于能较好反映降雨过程的平均特性，在缺乏当地降雨资料时使用也比较方便，成为工程应用最为广泛的一种设计雨型[17]；对于芝加哥雨型的分配计算，已知当地暴雨强度公式及参数时，雨峰位置系数是其分配计算过程中存在的主要不确定性参数. 对于产流、汇流和管网水动力计算方面，现有的模型参数灵敏度分析结果表明：产流模型的不渗透百分比、不渗透性洼地蓄水参数、渗透性洼地蓄水参数均对峰值流量、径流量较为敏感[18-19]. 对于透水地表，Horton模型是SWMM下渗计算的常用方法，在Horton模型参数对计算结果影响方面，最大下渗速率、最小下渗速率、衰减常数参数对峰值流量的影响较大[20]. 汇流模型的不渗透性粗糙系数、渗透性粗糙系数、子汇水区宽度对峰值流量、峰现时间以及径流量影响较大[21]. 管网水动力模型的管道粗糙系数、管径、进水偏移、出水偏移等参数对管道过流量和节点水深等结果影响较大[10,13].

根据上述文献的分析，本节总结了表1中所示的14类不确定性参数. 在现状排水管网模型分析中，对于不能通过实测资料直接或间接获取的参数，其取值范围主要参考室外排水设计规范[22]、SWMM用户手册[23]和相关的文献[17，24]. 对于模型参数的随机分布类型源自参考文献[25-28]. 其中，服从均匀分布的随机变量参数，可根据其参数取值范围确定分布参数；对于服从正态分布的随机变量参数，变量均值为参数取值范围的中值，对于变量标准差，假定该参数取值范围为其3σ界限，即随机变量有99.74%的可能性落在参数取值范围内.

2 参数不确定性分析

2.1 计算原理

根据结构可靠度理论[5]，将结构或系统的某个失效状态对应的功能函数表示为

Z(X)=R(X)-S(X)

(1)

式中：X=(X1,X2,…,Xn)为影响结构或系统状态的n个随机变量；R(X)表示结构或系统的抵抗能力；S(X)为外荷载作用下结构或系统的反应. 当Z(X)>0时，R>S，结构为可靠状态. 当Z(X)=0时，R=S，结构为安全临界状态. 当Z(X)<0时，R

则与功能函数Z(X)对应的结构或系统的失效概率为

(2)

式中f(X)为随机变量向量的联合概率密度函数.

对于本文关注的排水管网水力模型，分别定义节点溢流、管道满管为节点和管道的失效模式，对应的功能函数可以表示为

(3)

式中：X为影响模型(系统)状态的随机变量集；H0为模型节点对应的检查井井深；H(X)为模型计算得到的节点最大水深；D0为模型管道的直径；W(X)为模型计算得到的相对管底的管道最大水位.

由式(3)可以计算得到考虑不确定性工况下，节点、管道的失效概率为

(4)

同时根据公式

(5)

统计溢流节点i的溢流时间GTi(X)、溢流体积GVi(X)、管道j的充满度GFj(X) 的期望值.

2.2 Monte Carlo随机模拟

考虑模型参数的不确定性时，由于节点水深H(X)、管道最大水位W(X)是关于X的非线函数，且当X中的变量个数较多时，直接采用式(4)(5)求解节点和管道的可靠度和状态期望值是非常困难的，本节采用Monte Carlo随机模拟方法求解. Monte Carlo法是描述和估计函数统计量的一种常用近似解法，对于高度非线性功能函数，可以得到具有较高精度的近似解. 其所依据的原理为伯努利大数定理[29]：当试验次数N很大时，可以用事件发生的频率来代替事件发生的概率. 采用Monte Carlo模拟求解式(4)(5)时，相应值的估计值为

(6)

(7)

(8)

式中：N为模拟抽样的次数；x(k)为第k次模拟抽样得到的随机变量X的样本值；I(gnode(x(k))≤0)为{0,1}二元取值的示性函数，当gnode(x(k))≤0时，I=1，否则，I=0.

一般以式(9)所示估计值的变异系数作为衡量Monte Carlo计算准确度的依据：

(9)

在本文计算中，式(9)中G(X)不仅适用于GTi(X)、GVi(X)、GFj(X)，还适用于I(gnode(x(k)) ≤0)、I(gpipe(x(k)) ≤0).

对于包含n个随机变量参数X=(X1,X2,…,Xn)的排水管网模型，假定随机变量参数相互独立. 利用SWMM，并采用Monte Carlo模拟抽样计算的步骤如下：

步骤1建立研究区域SWMM. 根据已知条件、确定性参数、不确定性参数初始值，建立研究区域SWMM，存储为inp格式文本文件.

步骤2变量选择与计算设定. 选取对SWMM计算结果影响较大的n个重要参数作为模型的随机变量参数，并确定随机变量的分布类型(见表1)；设定Monte Carlo模拟次数N或模拟估计值变异系数的精度要求[δ].

步骤3变量随机抽样与模拟结果数值统计. 进行第k次模拟抽样，根据变量参数的随机分布类型分别产生n个随机数x(k)=(x1,x2, …,xn)，将这些随机数作为本次抽样的模型参数值写入模型inp文件中，调用swmm5.dll模型求解函数库进行模型计算，并读取rpt报告文件中节点溢流、管道充满度计算结果.

步骤4模拟估计值统计. 根据式(6)～(8)计算当前结果(前k次模拟后)对应的失效概率、溢流时间、体积、充满度的估计值，并根据式(9)计算模拟估计值的变异系数δk.

步骤5收敛终止判别. 当模拟次数k=N或δk≤[δ]时，终止抽样；否则，转步骤3.

模拟流程图如图1所示.

图1 Monte Carlo随机模拟的流程Fig.1 Flowchart of Monte Carlo simulation

3 案例分析

本文以北京市某研究区排水管网(见图2)为例，将参数不确定性分析模型应用到该区管网排水可靠性能评估. 该研究区总面积为3.1 km2，雨水管网设计重现期为3～5年一遇，共有检查井764个，管道766根(总长度22.36 km)，出水口7个. 河道防洪设计标准为50～100年一遇，出水口设置为自由出流. 采用表1所示的参数作为本案例模型的不确定性参数.

图2 研究区排水管网示意图Fig.2 Drainage network of the study area

3.1 参数灵敏度分析

在3年一遇设计降雨条件下，取表1中各参数的最小值为模型预设参数值，然后以参数取值区间值(最大值—最小值)的{0%，25%，50%，75%，100%}为预设参数值的5个增量值，研究参数不同增量值对模型计算结果(溢流节点个数、满管管道根数)的影响量，进而评价参数的灵敏度. 以管道粗糙系数C为例，单独取C值的5个增量值，其他参数同时取5个增量值时，模型计算结果中溢流节点个数、满管管道根数的变化情况如图3(a)(b)所示. 同理，分析表1中各参数对溢流节点个数、满管管道根数的变化情况. 图3(c)(d)为其他参数增量值为取值区间的50%时，各参数取值单独变化对溢流节点个数、满管管道根数的影响趋势图. 结果表明：C、Kw、r、Pi、O2、Dc对溢流节点个数、满管管道根数具有正相关关系，D、Rmin具有负相关关系，其他参数对结果的影响程度相对较小.

图3 参数对溢流节点、满管管道数量的影响Fig.3 Influence of parameters on the number of overflow nodes and full pipes

则采用以下方法将各参数进行灵敏度排序：在其他参数同时取5个增量值时，统计各参数在其取值范围内单独变化对应的溢流节点个数最大值与最小值之差，并将5个差值取绝对值相加；差值总和越大，其对溢流节点个数影响越大，其影响程度排序越靠前，由此得出各参数对溢流节点个数影响程度的排序. 同理可得各参数对满管管道根数影响程度排序. 将二者排序序号相加并按从小到大排序，排序结果如表2所示.

表2 3年一遇重现期设计降雨条件下各参数灵敏度排序Table 2 Order of sensitive parameters under the designed rainfall with a 3-year recurrence period

同理分析{3, 5, 10, 20}年重现期设计降雨条件下各参数灵敏度排序，排序结果见表3. 根据表3所示排序结果，本节以C、Kw、r、Pi、O2等5个影响较大参数为例，分析参数不确定性对模型计算结果的影响.

表3 不同重现期设计降雨条件下灵敏参数排序Table 3 Order of sensitive parameters under different return periods of designed rainfall

为验证本节选取的5个不确性参数的代表性，在{3, 5, 10, 20}年重现期条件下，选取5个和14个不确定性参数时，节点和管道失效概率等于1的个数统计如图4所示，由图中数据可知，选取5个不确性参数的分析结果与选取14个不确性参数相应结果的差距较小. 另外，模型计算结果中失效概率大于0.1的所有节点的溢流时间、溢流体积平均值表明：5个参数工况相对14个参数工况的变化率分别为{-10.29%, -10.73%, -15.84%, -0.20%}和{-9.79%, -8.09%, -9.13%, -0.27%}. 模型结果中所有管道充满度的平均值表明：5个参数工况相对14个参数工况的变化率分别为{-2.66%, -2.38%, -1.34%, 0. 0072%}. 因此，本节方法选择的5个不确定性参数具有较强的代表性，可以作为代表性参数对模型不确定性进行分析.

图4 节点和管道失效概率为1的数量统计Fig.4 Number of nodes andpipelines with failure probability of 1

3.2 不确定性模型计算结果分析

在不确定性分析中，模型考虑了C、Kw、r、Pi、O2共5个参数的不确定性，D取设计管径值，其他参数为表1中取值范围的中值. Monte Carlo模拟次数N设定为5 000次，为验证Monte Carlo模拟计算结果的收敛性，统计模型中所有节点或管道GT(X)、GVi(X)、GF(X)、Pf,node、Pf,pipe估计值的变异系数，最大值分别为{0.052,0.061, 0.012, 0.057, 0.047}(不计入溢流时间期望值小于0.01 min、溢流体积期望值小于0.01 m3的溢流节点)，说明Monte Carlo模拟结果具有较好的收敛性. 图5仅绘制了任选3个节点、3根管道在{3, 5, 10, 20}年降雨重现期下节点失效概率、管道充满度估计值随模拟次数的变化情况，并标注了相应的变异系数，可以发现曲线变化基本趋于平缓，且可认为模拟结果达到收敛.

图5 模型模拟值收敛性Fig.5 Convergence of model simulation values

1) 节点、管道可靠性计算结果

在{3, 5, 10, 20}年重现期条件下，分别统计节点、管道的失效概率如图6所示. 结果表明：①随着降雨强度(重现期)的增加，失效概率为0的节点和管道数量在逐渐减少，失效概率为1的节点和管道数量在逐渐增加. ②{3, 5, 10, 20}年重现期设计降雨条件下，分别有{92.15%, 87.04%, 63.87%, 44.24%}的节点和{79.68%, 67.17%, 35.90%, 30.88%}的管道长度的失效概率为0；这意味着，即便考虑参数不确定性时，该管网中的这些节点和管道仍不发生溢流或满管，为可靠节点或管道. ③{3, 5, 10, 20}年重现期设计降雨条件下，分别有{0.65%, 1.44%, 3.40%, 8.38%}的节点和{7.94%, 16.49%, 29.77%, 55.59%}的管道长度的失效概率为1；这意味着这些节点和管道肯定会发生溢流或满管，应采取措施提高这部分管道和节点的可靠度. ④对于失效概率介于0～1的节点和管道，存在着节点溢流、管道满管的可能性，应对失效概率较高的节点、管道给予关注.

图6 节点、管道失效概率分布Fig.6 Failure probability distribution of nodes and pipelines

2) 节点溢流时间、体积计算结果

统计全网中失效概率大于零的各个节点的溢流时间、体积均值见图7(a)(b). 结果表明：随着重现期的增加，节点溢流个数在逐渐增多，溢流时间、体积均值的分布范围在逐渐增大.

根据图7(a)(b)的均值区间划分，统计各区间的节点在5 000次抽样模拟中，溢流时间、溢流体积的变异系数最大值，见图7(c)(d)；结果表明：对于溢流时间均值大于{10 min, 1 min, 0.1 min}的节点，溢流时间变异系数最高可达{0.092, 1.069, 2.149}；对于溢流体积均值大于{1 000 m3, 100 m3, 10 m3}的节点，溢流体积变异系数最高可达{0.217, 1.162, 2.292}. 随机选取在{3, 5, 10, 20}年重现期设计降雨条件下均出现溢流的部分节点，分别将各节点在5 000次Monte Carlo抽样中得出的溢流时间、体积进行从大到小排序，首尾各去除总个数的2.5%，统计在置信度为95%的条件下，各节点溢流时间、体积波动区间见图7(e)(f)，发现：随着降雨强度的增加，节点溢流时间、溢流体积模拟值波动区间逐渐增大，说明参数不确定性对模拟值的影响越来越大.

以溢流时间、体积均值最大的节点为例，分析其在5 000次Monte Carlo抽样中，溢流时间、溢流体积的频度分布情况见图7(g)(h)，其中图例项括号内标注了模拟值的变异系数. 结果表明：随着降雨强度的增加，频度分布由集中趋于分散；而变异系数却由于均值的增加而逐渐减小.

图7 节点溢流时间、溢流体积统计Fig.7 Statistics of node overflow time and volume

结合以上分析发现：①溢流时间、体积均值越小，变异性越大. ②重现期越高，不确定性分析对节点溢流时间、溢流体积值的影响越大. ③降雨强度越大，溢流时间、体积频度分布越趋于分散. 因此，应综合衡量模型要求的计算精度与可以接受的变异系数，适当考虑变异性的影响.

3) 管道充满度计算结果

不同降雨强度下，全网中管道充满度模拟均值的频数分布见图8(a)，结果表明：随着降雨强度增加，低充满度管道根数减少，满管管道根数增加，而管道充满度模拟均值的分布范围逐渐减小. 在{3, 5, 10, 20}年重现期条件下，统计各管道充满度值在5 000次抽样模拟中的变异系数见图8(b)，结果表明：在{3, 5, 10, 20}年重现期条件下，变异系数超过0.05的比例分别为{9.14%, 9.40%, 11.23%, 5.74%}，超过0.1的比例仅为{1.57%, 2.48%, 2.61%, 1.04%}. 可认为管道充满度模拟值离散性较小.

统计各管道在5 000次Monte Carlo抽样中得出的充满度最大、最小值，绘制各管道充满度波动区间如图8(c)所示，结果表明：随着降雨强度的增加，管道充满度模拟值波动区间先增大后减小. 随机选取某一管道，分析其在5 000次Monte Carlo抽样中，充满度的频度分布情况如图8(d)所示，结果表明：随着降雨强度的增加，频度分布在{3, 5, 10}年重现期下由集中趋于分散，而20年重现期下由于管线充满度值最大为1又再次集中，这也是导致8(c)中模拟值波动区间先增大后减小的原因.

图8 管道充满度统计Fig.8 Statistics of pipeline filling

结合以上分析发现：①管道充满度值的离散性较小，说明变化区间相对均值的变化较小，建议采用均值统计管道充满度结果. ②随着降雨强度的增加，管道充满度波动区间先增大后减小，频度分布先分散后集中.

3.3 不确定性模型与确定性工况计算结果比较

将本文3.2节中考虑不确定性的模型计算结果与确定性参数工况结果进行比较. 其中，确定性工况中，C、Kw、r、Pi、O2等5个参数取表1中取值范围的中值，其他参数设置与不确定性模型相同. 计算{3, 5, 10, 20}年重现期下节点溢流、管道满管结果，根据计算结果，将节点、管道进行如下分类：A1为确定性工况中溢流的节点；A2为确定性工况未溢流且在不确定性模型分析时失效概率大于0的节点；B1为确定性工况中满管的管道；B2为确定性工况未满管且不确定性模型中失效概率大于0的管道.

1) 不确定性模型失效概率与确定性工况计算结果差异性

统计上述4类节点、管道在不确定性模型中计算出的失效概率结果见表4. 从表中数据表明：①对于A1类节点和B1类管道，其中部分节点和管道的失效概率较低，小于0.5的节点和管道个数占对应类别总个数的比例最高分别为14.85%和4.86%；可见，相对于确定性工况，不确定性模型能够避免部分低风险节点和管道计算结果为溢流或满管. ②对于A2类节点和B2类管道，其中部分节点和管道失效概率较高，大于0.5的节点和管道个数占对应类别总个数的比例最高分别为8.57%和8.44%；可见，相对于确定性工况，不确定性模型能够识别潜在高风险节点和管道. 综合以上两点可以发现：参数取值的不确定性直接导致模型模拟结果的差异性，有必要对模型参数进行不确定性分析；不确定性模型对高风险概率的节点和管道具有较好的识别效果.

表4 不同可靠度范围节点个数占对应类别节点总个数比例Table 4 Number of nodes in different reliability ranges accounts for the proportion of the total number of nodes in the corresponding category %

2) 节点溢流时间、体积计算结果差异性

对于A1类节点，不确定性模型得到的各个节点溢流时间、体积模拟均值与确定性工况下相应结果的统计特征箱形图见图9(a)(b). 由图9(a)可知：在各重现期条件下，相对于确定性工况计算结果，不确定性模型得到的溢流时间均值的分布范围较大(离散性大)，均值、中值较小. 比较中值线的位置可以发现，不确定性模型结果中多于50%的节点溢流时间均值小于确定性工况节点溢流时间的最小值. 以上分析说明，采用确定工况评价A1类节点得到的溢流时间是偏大的(保守的). 该结论同样适用于图9(b)节点溢流体积的比较.

对于A2类节点，不确定性模型得到的各个节点溢流时间、体积均值统计见图9(c)(d). 发现：在{3, 5, 10, 20}年重现期条件下，溢流时间均值的最大值分别为{0.24 min, 0.38 min, 0.45 min, 0.58 min}，溢流体积均值的最大值分别为{1.63 m3, 5.21 m3, 1.62 m3, 1.28 m3}，相对确定工况全部溢流节点溢流时间、体积平均值的变化率分别为{11.90%, 16.33%, 22.83%, 29.13%}和{5.91%, 10.41%, 3.98%, 2.91%}.

图9 节点溢流时间、溢流体积比较Fig.9 Node overflow time, overflow volume comparison

3) 管道充满度计算结果差异性

不确定性模型得到的管道充满度均值相对确定工况管道充满度比较结果如图10所示. 结果表明：①对于B1类管道，在各重现期条件下，充满度均值基本处于在0.89以上，说明2类工况评价结果是总体相符的，而确定性工况计算的管道充满度是偏保守的. ②对于B2类管道，分别统计每根管道在不确定性模型中的均值相对确定性工况的变化值和变化率，如图9(b)所示，发现与确定性工况相比，不确定性模型模拟结果的相对变化值最高可达0.15，相对变化率最高可达25%.

图10 管道充满度比较Fig.10 Pipe filling degree comparison

4 结论

1) 降雨输入模型的雨峰位置系数、产流模型的不渗透百分比、汇流模型的子汇水区宽度、管网水动力模型的管道粗糙系数、出水偏移等参数对节点溢流、管道满管计算结果影响较大，建议考虑模型全过程参数的不确定性.

2) 对于不确定性分析结果中节点溢流时间、体积结果的统计与表征，建议在考虑模拟均值(估计)大小的同时，考虑模拟值的变异性.

3) 不确定性分析对高风险概率的节点和管道具有较好的识别效果；对于溢流时间、体积计算结果，仅不确定性分析识别出的风险节点均值相对确定性全网均值的变化率分别达到了29.13%、10.41%，对于管道充满度值计算结果，相对确定性工况的变化率达到了22.31%. 建议对排水管网进行可靠性分析.