排水管道排放污染物贝叶斯统计算法-SWMM耦合溯源模型

摘 要

城市排水管道中广泛存在着超标工业废水偷排漏排进入市政管道，导致污水处理厂运行效率下降甚至活性污泥大批死亡，造成处理水质不达标，排放水体水质恶化的严重污染事件。现有的物理溯源方法效率低下反馈时间长，不利于偷排漏排事故的快速有效识别。而由于排水管道系统拓扑结构及水力条件复杂多变，其环境反问题具有高度不确定性，成为常规数学模型溯源方法应用于排水管道中的一大难题。针对此现状，该研究以排水管道中偷排工业超标污染物质为研究对象，开展了基于在线水质监测数据的排水管道偷排漏排污染物质溯源模型的研究，构建了基于贝叶斯统计推理算法、SWMM水力水质模型以及Matlab编程的排水管网污染物质溯源数学模型，对偷排超标污染物质的排放位置、排放量和排放时段等三个未知的排放特征参数进行了统计反推，获得了排放节点、排放量和排放时段的概率分布；并重点研究了统计反演算法中游走步长的取值、采样监测的时间间隔及在线监测点的空间布置等实际限制条件对反演结果精度的影响；反演模型采用了先进的马尔科夫蒙特卡罗（MCMC）抽样算法，大大提高了溯源的效率，在有限的反馈时间内可提供更精确的溯源结果，从而提高对偷排漏排超标废水的管理水平以及排放事故的应急处理效率。主要研究内容和结论如下：

① 开展了SWMM水质模型、贝叶斯统计方法及MATLAB软件的耦合集成研究。结果表明，MATLAB强大的数据处理功能，能够实现MCMC抽样算法涉及的大量计算，提高反演效率；SWMM与贝叶斯耦合的数学模型可对排水管网未知污染源进行有效的识别追踪，并得出有实际指导意义的排放特征参数的范围，为物理溯源提供更精细的搜索路径和方向。

② 贝叶斯统计MCMC抽样算法中游走步长与反演效率和精度的关联关系研究。结果表明：游走步长过小可能使抽样无法在整个后验空间搜索，使抽样结果陷于局部解，无法准确反演污染源参数的后验概率密度；适当游走步长取值可以实现对污染物反演参数整个后验空间的抽样，避免使反演参数限于局部解，提高反演效率和反演精度；但是步长过大会使抽样样本跳出后验空间，只有很少的样本落在了后验空间，使抽样结果精度降低，若要在过大步长下提高精度就必须增加MCMC链长，这样就会导致抽样时间延长，抽样效率降低。

③ 采样监测的时间间隔及在线监测点的空间布置对反演算法有效性和精度的影响研究。结果表明：监测点的观测数据过少时，无法准确反演污染源参数的后验概率密度；当观测数据增多、采样监测时间间隔比较小时，观测数据没有反映监测点的浓度序列，此时也无法准确反演污染源参数的后验概率密度；只有适

I

重庆大学硕士学位论文

当的采样监测时间间隔才能使反演结果的精度和效率更高；水质监测点增加时可提高溯源效率。

④ 通过构建实际案例溯源模型，对节点瞬时排放和节点连续排放污染源源项单变量参数、双变量参数和三变量参数分别进行抽样反演，并将抽样结果和理论解对比来评估统计溯源模型的精度，结果表明：在污水管网污染事件中，只有一个未知反演参数时反演效果最佳，可以精确确定其发生的范围；两个和三个反演参数时都能获得合理的污染源源项范围，但是反演参数越多，污染物排放量、排放时刻和排放节点的后验概率密度越离散。

研究构建的溯源模型可模拟污水管道系统中保守及非保守污染物质的排放、迁移、转输全过程，并基于贝叶斯理论对排水管网中可能的偷排漏排污染点源进行统计意义上的反推并给出各个可能取值区间的概率，从而筛选出可能的污水排放节点、排放量和排放时间。通过贝叶斯方法、Matlab与SWMM的耦合集成，水质预测结果与实际管网的动态运行工况更接近，提高污染物传输模拟的精度，从而提高反演分析的精度。同时贝叶斯统计反演方法充分考虑了先验信息、观测噪声以及模型误差的影响，在泄漏源反算的不确定性分析中独具优势，能较好应对排水管网系统溯源反问题的不确定性。提供有实际指导意义的溯源结果。论文的研究成果可为水体污染控制、智慧排水监管系统等领域提供技术支持。

关键词：城市排水系统，污染物溯源，贝叶斯方法，MCMC抽样，SWMM模型

II

英文摘要

ABSTRACT

In the urban drainage system, accidents caused by the illicit discharge of industry waste water to the municipal sewage system frequently occur. They often have severe impacts on the activated sludge process in the sewage treatment plant, and even cause deterioration of effluent quality due to activated sludge poisoning. Hence it is important to identify illicit discharge location and discharging history. Existing source identification methods based on physical approaches have low efficiency and take a long feedback time, which is insufficient in the identification of the drain accidents. Meanwhile, limited source inference methods based on mathematic model have difficulties in the application in the drainage pipe system, in which source inference problem is highly uncertain due to the complex topology and hydraulic conditions. Hence, in this research, the illicitly discharged industrial pollutants in the drainage system was taken as the research object. A mathematical source identification model based on Bayesian inference algorithm, SWMM model and MATLAB programming in drainage networks was constructed. The three unknown discharge characteristic parameters including discharge location, discharge amount, and discharge time period, were statistically retrieved, and the probability distributions of the variables were obtained. The influence of inference parameters such as the value of the walking step size in the statistical inversion algorithm、the spatial arrangement of the on-line monitoring points and the time interval of the sampling monitoring, etc., were studied in detail. The model used the advanced Markov Monte Carlo. The MCMC sampling algorithm greatly improves the efficiency and provides more accurate results within a limited feedback time. Hence it could improve the management level of waste water and the emergent treatment efficiency in the discharge accidents. The main research contents and conclusions are as follows:

① The integration of SWMM, Bayesian method and MATLAB was carried out. This integration can output the time series of SWMM simulation results through MATLAB and match the corresponding monitoring data sequence to improve the analysis accuracy. Hence it could provide a dynamic feedback that is closer to the real condition of the pipe network. Results show that the coupled source identification model can provide a useful range of the unknown discharge source parameters and enhance the efficiency of a physical-based model.

III

重庆大学硕士学位论文

② Bayesian MCMC sampling algorithm was applied to study the relationship between walking step size and inversion efficiency. Results show that if the step size is too small, the sampling will be limited to the local solution while a large step size can move in a larger search space and avoid being confined to a local solution. Hence a larger step size can realize the sampling of the whole posterior space. However, a step size that is too large would increase the number of samples and hence an extended inversion time.

③ Influence of monitoring point spatial arrangement and monitoring time interval on the effectiveness and accuracy of inversion algorithm was investigated. Results show that when the observation data of the monitoring point is too small, the posterior probability density of the pollution source parameters cannot be accurately inverted. This is because the observation data does not reflect the concentration sequence of the monitoring point when observation data increase and monitoring time interval is relatively small. A more complete pollutant concentration curve means a longer monitoring time. Hence it can provide a more accurate Inversion results.

④A real drainage network was taken as a study case using the source-identification model developed. Three scenarios that reflects the availability of the source information are investigated including a single-variable, double-variable, and three-variable that represents one, two and three discharge characteristics of the discharge accident needs to be identified. Two discharging mode was studied including an instantaneous and a continuous discharge. Inferenced source information’s were compared with the theoretical solution to evaluate the accuracy of the statistical traceability model. Results show that the model provides accurate inversion results in cases that involves only one unknown inversion parameter. However in cases where two or three inversion parameters needs to be inferred, the model can provide a reasonable range of source information’s. But more inversion parameters, more discrete of the posterior probability density。

Since the purpose of the Bayesian inference method is to find out the possible values and give the probabilities of all the values, the basic idea is to narrow down the prior range and focus on the points where the probability distribution is large. Therefore, the result of SWMM-Bayesian inference model is reasonable and practical. The research results of the paper can provide technical support for water pollution control and drainage management .

IV

英文摘要

Keywords: Urban drainage system， Pollution source identification， Bayesian

inference， MCMC， SWMM

V

目 录

目 录

中文摘要 .......................................................................................................................................... I 英文摘要 ....................................................................................................................................... III 1 绪 论 ...................................................................................................................................... 1

1.1 选题背景及意义 ...................................................................................................................... 1 1.2 国内外研究现状 ...................................................................................................................... 3 1.2.1 基于物理搜索的污染物溯源方法 ................................................................................... 3 1.2.2 基于数学模型的污染物溯源方法 ................................................................................... 4 1.3 污染物统计法数学模型溯源 .................................................................................................. 9 1.4 研究目的、内容及技术路线 ................................................................................................ 10 1.4.1 研究目的 ......................................................................................................................... 10 1.4.2 研究内容 ......................................................................................................................... 11 1.4.3 技术路线 ......................................................................................................................... 12

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建 ........................... 13

2.1 单纯贝叶斯反演模型构建 .................................................................................................... 13 2.1.1 贝叶斯定理 ..................................................................................................................... 13 2.1.2 先验分布函数定义 ......................................................................................................... 14 2.1.3 似然函数构造 ................................................................................................................. 14 2.1.4 后验概率密度函数构造 ................................................................................................. 15 2.1.5 贝叶斯反演过程 ............................................................................................................. 15 2.2 MCMC抽样Matlab抽样 ..................................................................................................... 15 2.2.1 MCMC抽样算法 ............................................................................................................. 15 2.2.2 MCMC Matlab编程 ........................................................................................................ 16 2.3 贝叶斯MCMC反演模型构建 ............................................................................................. 18 2.3.1 贝叶斯-MCMC反演模型Matlab编程 ......................................................................... 18 2.3.2 贝叶斯-MCMC反演模型数值验证 .............................................................................. 21 2.3.3 游走步长及马氏链长对贝叶斯反演模型的影响 ......................................................... 21 2.4 本章小结 ................................................................................................................................ 26

3 SWMM、贝叶斯与MATLAB耦合的污染物溯源模型 ........................ 27

3.1 SWMM水质模型 ................................................................................................................... 27 3.2 SWMM模型集成技术难点 ................................................................................................... 27

VII

重庆大学硕士学位论文

3.3 SWMM、贝叶斯及MATLAB耦合污染物溯源模型 ........................................................ 29 3.4贝叶斯-SWMM污染物溯源模型的特点 ............................................................................. 30 3.5 本章小结 ................................................................................................................................ 31

4 排水管道污染物溯源模型的验证与应用 ......................................................... 33

4.1 仿真数值试验设计 ................................................................................................................ 33 4.1.1 节点瞬时排放事件数值试验设计 ................................................................................. 34 4.1.2 节点连续排放事件数值试验设计 ................................................................................. 36 4.2 节点瞬时排放污染事件溯源 ................................................................................................ 39 4.2.1 单个未知参数反演 ......................................................................................................... 40 4.2.3 两个未知参数反演 ......................................................................................................... 59 4.2.4 三个未知参数反演 ......................................................................................................... 67 4.3 节点连续排放污染事件溯源 ................................................................................................ 74 4.3.1 单个未知参数反演 ......................................................................................................... 75 4.3.2 两个未知参数反演 ......................................................................................................... 79 4.3.3 三个未知参数反演 ......................................................................................................... 85 4.4 溯源结果分析与讨论 ............................................................................................................ 88 4.5 本章小结 ................................................................................................................................ 89

5 结论与建议 ......................................................................................................................... 91

5.1 结论 ........................................................................................................................................ 91 5.2 建议 ........................................................................................................................................ 92

致 谢 ...................................................................................................................................... 95 参考文献 ...................................................................................................................................... 97 附 录 .................................................................................................................................... 101

A. 作者在攻读硕士学位期间发表的论文 ............................................................................... 101 B. 作者在攻读硕士学位期间参加的相关科研项目 ................................................................ 101 C. 作者在攻读硕士学位期间申报与获批的专利 ................................................................... 101

VIII

1 绪 论

1 绪 论

1.1 选题背景及意义

城市排水系统中频繁发生工厂偷排漏排超标废水对污水处理厂活性污泥法工艺造成严重冲击的污染事故，甚至导致活性污泥中毒死亡出水水质恶化，污染城市水环境。

重庆市长寿污水处理厂，采用氧化沟工艺,设计规模为8万立方米/日，处理规模为4万立方米/日，该污水处理厂在2004年期间常遭受到不明工业废水的“侵袭”，致污水处理系统处于崩溃边缘，如在2004年5月21日，由于被不明工业废水冲击，造成活性污泥全部中毒死亡，污水处理系统瘫痪；河北省沧州市运东污水处理厂采用微孔曝气氧化沟处理工艺，规模为10万立方米/日，日产中水3万立方米，在2015年07月22日由于不明化工企业偷排废水，进厂水质指数严重超标，导致污水处理系统处于崩溃边缘；陕西省咸阳市东郊污水处理厂，采用A2/O+纤维滤池处理工艺，其计规模为17万立方米/日，处理规模为10万立方米/日，其中工业废水占40%，由于不明工业企业向管网偷排超标污水，使在凌晨3～5点期间污水处理厂经常出现进水CODCr超过设计标准2～3倍的峰值现象，这种间歇式的冲击负荷在很大程度上影响了污水处理厂的正常运行，加大了运行成本。

图1.1 工厂偷排漏排 图1.2 污泥中毒上浮 Fig1.1Illcit discharge Fig1.2 Sludge poisoning

工厂排放的废水是污水管网中有毒物质的主要来源，有毒物质包含挥发酚、氰化物等非金属有毒污染物和重金属有毒污染物。这些难降解的有毒物质进入污水处理厂就会对活性污泥法工艺造成严重冲击。工业废水一般占城市污水的30%～70%，有的甚至高达80%。工业废水中常见的有毒有害物质如表1.1所示。

1

重庆大学硕士学位论文

表1.1工业废水有毒物质来源

Table1.1 The source of toxic substances of industrial wastewater

有害物质 酸 碱 汞及其化合物 镉及其化合物 六价铬及其化合物 砷及其化合物

酚 氰化物

铅及其化合物

油 硫化物

工业废水主要来源

化工、矿工、钢铁、有色金属冶炼、机械、电镀工业等 化纤、制碱、造纸、印染、皮革、电镀工业及石油炼厂等 氯碱、汞制剂农药、化工、仪表、电镀、汞精炼工业等 金属矿山、冶炼、电镀、金属处理、电池、特种玻璃工业等 矿山、冶炼、电镀、金属处理、电池、特种玻璃工业等 矿石处理、冶炼、化工、玻璃、涂料、农药、化肥工业等 焦化、煤气、煤油、合成树脂、化工、染料、制药工业等 焦化、煤气、电镀、金属清洗、有机玻璃、丙烯睛合成、煤油工

业及黄金工业等

冶炼、化工、农药、汽油防爆、含铅油漆、搪瓷工业等 煤油、机械、食品加工、油田、天然气加工工业等 化工、皮革、煤气、焦化、染色、黍占胶纤维、煤油、油田、天

然气加工工业等

游离氯 有机磷、有机氯

造纸、织物漂白、化工工业等

农药、化工工业等

工业废水的成分和性质具有复杂多变性，再加管理水平相对较低，极易导致携带有毒有害物质的工业废水排入污水厂，对污水处理厂造成严重的冲击负荷，主要有以下三种风险：① 超标排放。高浓度大量的工业企业废水排入污水厂，严重超过污水厂的设计负荷，导致工业废水降解率不足，出水指标高于设计排放标准，污水处理厂出水水质不达标；② 污泥中毒。污水处理厂的污泥在长期驯化已适应一定的水质要求，当工业企业废水进入污水处理厂时，活性污泥不能在短时间内适应，特别是当大量高浓度或含有毒有害物质的工业企业废水排入时，就会导致污泥活性降低、中毒、解体、甚至死亡的后果；③ 环境残留。重金属和其它难降解物质进入污水处理厂后，不能被活性污泥降解，污水处理厂的出水中仍含有该重金属等难降解物质。

目前市政管理部门缺乏快速有效的污染物质监管控制方法，通常当污水处理厂的活性污泥接触到有毒有害物质并造成损害时，才发现工业企业超标废水进入污水处理厂。甚至导致活性污泥中毒死亡而必须更换并重新驯化污泥。在此期间，废水的有效处理将无法进行，造成新的环境污染，经济及环境损失巨大。如果含有毒有害等难降解物质的工业废水在进入污水处理厂之前能够被预警检测出，可

2

1 绪 论

有效防止活性污泥工艺遭受损坏。尤其是对中小城市而言，工业企业经营管理水平和加工技术相对较低，污水排放事故频繁发生，甚至一些工业企业为了节省生产成本而偷排漏排超标废水，导致难以处理，极大影响了污水处理厂的正常运行，处理效果不稳定，出水水质受到影响。

面对此类突发的污染排放事故，实时准确的定位污染排放源成为事故应急救援中亟待解决的关键问题。在突发情况下，要求在污染源排放位置、强度及排放时间等信息均缺失的有限条件下，以最短的时间确定未知排放源停止排放，为应急决策提供依据。因此，研究非法排放事件下的污染物定位及其源强反算，不仅可及时进行处理，且可在水污染发生后及时对水污染企业进行责任追究，提高环保部门对水环境质量的监管水平，防治水污染事件发生，更好地保护城市水环境。

城市排水系统污染物质溯源，也称为污染源识别定位问题，或者称作事故重构，源反演，或源参数识别，主要是基于监测的浓度数据，追踪定位偷排漏排进入排水管道的污染物质的来源，寻找出污染源排放节点、排放时间、排放强度等关键信息。因此，在规范排水管网监控管理、快速处理突发性污染事故以及水环境质量保障方面具有重大的意义。

1.2 国内外研究现状

20世纪70年代，美国、德国等发达国家开始从对江河湖泊等地表水的自动监测向对市政排水管道在线监测发展，水质监测除常规指标水温、PH值、电导率、浊度等外也涉及到总磷、总氮、总有机碳、生物需氧量、化学需氧量、铅、汞等重金属指标。水质自动监控系统可以快速分析水样，并通过有线或无线的传输数据方式将快速分析得到的数据传输至检测中心，以便管理者及时掌握管网的运行状况，并根据监测数据优化调控管网。我国过去一直专注于对城市污水处理技术的研究，随着国家环保监管力度的加大和“水十条”的发布，我国在城市排水方面管理方面也逐步加强，也开始对排水管道的水质进行监测。

针对排水及自然水环境中的突发水质污染事故，国内外学术界及工程界提出了多种污染事故溯源方法。目前采用的污染物质溯源方法主要分为物理溯源和数学模型溯源两大类。

1.2.1 基于物理搜索的污染物溯源方法

物理溯源是通过物理工具现场监测，通过树枝状管网分叉处上、下游节点流量、水质、连通性进行监测分析，逐一排查，直接确定污染物质排放的方法。目前基于物理方法溯源按照其实现方式主要分为以下几种：

① 示踪法

这类方法主要包括同位素示踪法、微生物示踪法等，通过用示踪剂对研究对

3

重庆大学硕士学位论文

象进行标记，分析示踪剂分布和转移推演出研究对象的来源和转移规律[1]。方晶晶等人[2]使用多元素同位素及微生物示踪法，对河北平原地下水中污染物硝酸盐的来源和迁移转化规律进行追踪，揭示了地下水中硝酸盐的来源。

② 仪器搜索定位法

该方法通过使用红外线照相机、水面或水下自主机器人等对污染物进行追踪分析。例如Mathieu Lepot等[3]应用红外线照相机对非法连接和入渗的下水道进行了检测和侧向量化；Jay A．Fan'ell等[4]，将智能算法嵌入水下机器人，将追踪污染物的过程划分成：发现污染带、跟踪、偏离、重新发现和确认污染源等行为，通过对污染物浓度实时检测来判断当前所处的行为，然后通过对比来规划水下机器人前进方向和步长，通过行为之间的不断切换，最终实现污染物追踪溯源的目的。在加州圣克利门蒂岛附件海域通过使用此算法的水下机器人实现了追踪污染源的功能。

③ GIS地理信息技术溯源法

由于GIS具有强大的空间分析能力，如果排水管网信息以GIS电子网络的形式存放，则可以通过GIS的空间节点属性，基于管道节点流量平衡以及污染物质化学质量平衡方程，进行逆向拓扑分析，逐级溯源，推算可能的排放位置。但是由于可能排放的节点众多，单独根据GIS功能无法准确定位，通常需结合现场物理方法进行筛选识别。

1.2.2 基于数学模型的污染物溯源方法

污水管网分支繁多及污染事故发生的随机性和非连续性，物理方法寻找污染源的效率低下且信息反馈滞后，不利于应急处理。与此同时，一些排污单位为了躲避监测，采用接入清水稀释管道或者明暗两套排污口等迷惑性措施，使得事故性污染源往往无法准确被物理方法所辨识。随着水质模型在排水系统工程中的广泛应用，排水系统数学模型溯源也逐渐开始发展[5]~[10]。如郑伟等[11]将ARCGIS与一维水质模型进行集成，以重庆渝北区为研究对象，构建了城市污水管网溯源监控系统，该系统不仅具有传统污水管网信息管理的功能，而且可以检测污水管网的水质，当污水管网水质超标时能够及时预警，并对潜在的超标污水排放工业企业进行溯源分析。

水质模型可在较短时间内提供大量不同排放情景的模拟结果，为突发性水污染事故溯源提供数据支撑[12]~[17]。

① 数学模型溯源原理

数学模型溯源的原理是求解排水管道中污染物的对流传输方程的反问题，其核心是利用水质模型计算预测污染物质在水中的迁移转换规律。水质模型是水体水质变化规律的数学描述，排水管道通常可以视为一维流动，污染物质在其中的

4

1 绪 论

对流传输服从一维对流传输方程，其控制方程如公式1.1所示，表示污染物质沿着管道长度方向的对流、传输及降解全过程。

∂(AC)∂t

+

∂(QC)∂x

=∂x(AEx∂x)−AK1C+S(x,t,M） （1.1）

∂∂C

其中：A------过水断面面积，（m2）； Q------流量，（m3/h）；

C------污水管网中污染物的浓度分布（mg/L），是时间和位置的函数

C(x,t)；

t-------时间，（s）； x-------流程，（m）；

Ex------纵向离散系数，（m2/s）； K1------污染物一级降解系数，（s-1）；

S（x,t,M）------ 污染点源排放随时间变化的函数。

泄漏源定位问题，也称作源项反演或源参数识别等，是在方程1.1中浓度C(x,t)分布已知的前提条件下，反推其中的源项S(x,t,M）。C(x,t）来源于监测点所记录的污染物质浓度C。

方程1.1是一个偏微分方程，通常采用数值的方法进行求解。求解该方程的前提是流场已知，也即是流量或者流速等参数已经确定。排水管道水质模型通常综合水文、水力的计算来获得管道沿程的水量和流速分布，然后再数值离散求解方程1.1，是一个复杂耦合的过程。

标志着城市排水管道数学模型开始的是1971年由美国EPA开发的SWMM模型。SWMM包含水文、水力及水质三大计算模块，既可以对天然降雨进行模拟，也可以对节点定义的流量输入进行模拟。SWMM内嵌的一维水质模型可以在给定化学组分降解系数的前提下，将化学物质定义在节点排放，在管道系统沿程求解一维对流扩散方程，从而模拟化学物质在整个管道系统内部稀释、传输、降解全过程。其它基于SWMM 开发的一系列的城市排水软件例如XP-SWMM、Infowork ICM，Mike Mouse等软件对水质的模拟过程类似，被大量用于城市污水管道的流量计算及管道设计计算，以及合流制的溢流污染风险评估以及管控方案规划中[18]。杨雪[19]基于SWMM水质模型建立了合流制溢流污染规律研究模型，分析了溢流污水初期效应、溢流动态过程，以及不同降雨重现期、不同截流倍数溢流污水水质的变化及管道沉积物的影响；黄兵等[20]利用SWMM水质模型建立了昆明市主城区船房河流域合流制排水管网模型，通过对昆明市两场典型降雨进行模拟，分析了其在污水入流、雨水入流、时间峰值三者叠加共同作用下对合流制排水管网的影响；孙全民等[21]应用SWMM模型对柳州市竹鹅溪片区截流式合流制管网溢流进行了水质水量模拟。

5

重庆大学硕士学位论文

② 数学模型溯源技术难点

利用数学模型溯源实质上是逆向求解一维污染物对流传输方程(公式1.1)，也即是是求解对流传输过程的反问题。反问题是相对于正问题而言的，数学上通常把研究得相对比较充分或完备的问题称为正问题，与此相对应的另一个问题称为反问题[22]。在本研究中，正问题就是求解一维污染物对流传输方程1.1所描述的污染物在水体中混合迁移转化规律的问题，具体来说就是指在确定边界条件（上游及下游边界处浓度值）和初始条件下（也即是污染物的排放特征函数），求解污染物随时间和河流长度方向的变化规律。相应的反问题则是根据已知污染物时空分布的部分信息逆向反推初始条件，也即是污染物的排放特征函数S（x,t,M）。因此，数学溯源模型就是基于污染事故发生后所监测获取的污染物浓度分布C（x,t）对污染源函数S（x,t,M）中的三个未知参数x，t，M进行反向推导，进而计算出污染源的排放节点x、排放时间t和排放强度M等参数，从而实现识别查找未知的污染源的目的，这个过程也称为参数反演或者源项参数的反演。

对于污染物质的对流传输方程而言，正问题通常是适定的（well-posed），也即是根据一组提供的边界条件和初始条件，可以唯一获得污染物质沿河流随时间的变化分布。然而其反问题，也即是源项反演溯源问题是不适定的（ill-posed），宏观表现为一组给定的监测数据可能存在着多个污染源函数。

数学模型溯源问题不适定的根本原因在于在污染物质溯源过程中，通常监测节点数少于总节点数，所以控制方程的总数少于未知量数目，这是由有限的监测结果反推未知的排放源项参数问题本身所具有的特性所决定的，其结果就是使得数学模型的逆向求解在数学上是不适定的，表现为解的非唯一性或者不确定性。这成为困扰数学模型溯源的一大难题。

③ 数学模型溯源主要算法

由于污染物质的源项反演问题的不适定性，污染源识别问题成为国内外环境流体力学研究的难点和热点，目前的研究主要集中在给水管网中非法入侵污染物、地下水及天然河道突发污染事件识别，所采用的主要求解方法有直接法、优化方法和统计方法。

1）直接法

直接法求解一般是利用Tikhonov正则变换，将不适定性反问题转化为适定的问题来进行求。Tikhonov正则化方法研究成果主要集中在数学物理反问题领域

[23][24]

，在水环境模型参数估计也有一些初步的研究成果。闵涛等[25]应用Tikhonov

正则化方法，基于一维河流水质模型以及实际监测水质指标，对河流水质纵向弥散系数进行了反演估计。与此类似，栗苏文等[26]在基于一维河流水质模型 Dobbins BOD-DO 耦合模型的基础上，应用导数-正则化参数辩识方法，对Dobbins BOD-DO

6

1 绪 论

耦合模型中的耗氧系数、复氧系数、沉淀再悬浮系数以及光合作用产氧/生物呼吸耗氧速率等4个水质模型参数进行了反推。

Tikhonov正则化方法从理论的角度求解反问题，其求解过程复杂，通常适用于理论解已经知道的环境反问题，在实际工程中应用相对困难。

2）最优化方法。

最优化方法是基于系统优化原理，将源项反问题转化为寻求正问题参数的最优解来求解。最优化方法的出发点是将目标函数设为监测数据与模拟数据之差，应用最优化方法寻优算法使得目标函数最小化。在反问题求解过程中，目前主要采用的寻优算法包括脉冲普法、遗传算法(GA) [27][28]、模糊推理系统、贪心算法[29]、关系树—线性规划法[30]、模拟-优化法[31]、人工神经网络[32][33]等方法。

脉冲谱法是一种理论比较严密的半解析方法，其基本思想是：在时间域中给出输入信息，在复频域中确定未知量。李兰[34][35]采用一维水质方程模拟了河流中可降解有机污染物质在水体中的输移、转化、累积过程，将紊动扩散、分子扩散、断面内流速分布不均引起的污染物质沿纵向的浓度变化，采用纵向混合系数进行描述，并采用了脉冲谱法对纵向混合系数进行了逆向的推理反演。

遗传算法(GA)是基于生物进化理论的一种寻优算法，在寻找最优解的过程中，模仿生物的过程，通过搜索方向复制、交叉、变异等操作，不断调整其寻优区域，经过多次迭代循环，最终实现，优胜劣汰，找到最优解或接近最优解。以目标函数作为搜索信息是遗传算法最大的优点，因为在寻优的过程不需求解目标函数的导数值这些信息。对于一些，工程实践过程中常见的高度非线性的反问题，显示出了更高的灵活性性。基于遗传算法为主[36]的污染源溯源已经被用于大气中的污染物反演识别研究。Ng等[37]运用遗传算法，求解一维河流水质模型，率定了河流水质模型参数。Banik[38]采用排水软件SWMM构建了城市污水管网模型，在给定的污水排放模式下，求管道中的流速及流量及污染物指标的变化规律。然后基于监测数据，利用遗传算法求解多目标优化问题，研究了以监控非法污染物入侵为目的的污水水质监测点的最优布置。

模糊推理系统对于复杂系统或不精确系统的模拟非常有用，可以将水质模型参数提炼成一个优化模型，但当优化的目标函数比较模糊时，就需要用模糊优化的手段。Subhankar等[39]基于灰色模糊优化模型建立了河流系统水质管理模型。冯耀龙等[40]基于河流水质系统中各目标冲突而又模糊的特点，建立了河流水质管理的模糊优化模型，以“满意度”作为其目标函数，在河流水体满足一定水质要求的前提下，为系统各方的竞争与协调提供了更大的方便和灵活性，因此更符合实际情况。

贪心算法是用局部解构造全局解，即从问题的某一个初始解出发，根据具体

7

重庆大学硕士学位论文

问题的优化准则，在每一步的选择时都取当前问题的最优解，并且每一步只选择一个监测点，在每一次选择时都满足当前条件下的最优条件，重复此过程，直至所有监测点选取完或达到算法终止条件[41]。陶涛等[42]针对城市供水系统中突发水污染事件，利用混合整数规划模型表达了多目标水质监测点布置，通过模拟管网发生的水质污染事件，计算管网发生污染事件时对用户带来的影响值，利用贪心算法对模型求解，并结合ZJ市管网，验证了该算法对模型的求解效率，提出了该市基于突发污染事件的管网水质监测点优化布置方案。

关系树—线性规划法是在数据训练的基础上，用关系树(RT)描述污水管网监测节点与其上游所有节点之间线性规律。若已知某节点的污染物浓度，根据节点关系树的线性规律可求出上游所以节点的最大污染浓度和其发生的时间，依次求解下去，直到求出所有相关节点的最大污染物浓度和发生时间，其中污染物浓度最大的节点即是污染物注入点。

模拟-优化法是将污染源排放点和排放时间归结为一个最优化问题，并结合梯度下降法进行优化计算。在突发性污染事件中，污染源的位置、污染物的强度以及污染事件发生时间等参数的确定可以转化成数学上的最优化问题。

人工神经网络（ANN）是一种通过数学方法对人和动物神经网络的某些结构和功能进行抽象模拟的方法，是一种模仿结构及其功能的线性信息处理系统。此方法也被应用于市政给排水系统中，如Zou等[43]等建立了考虑输入不确定性下水质模型的人工神经网络；信昆仑等[44]利用人工神经网络实现了对供水管网污染源的追踪定位。Kim[45]在丹麦排水管道软件CHIMOUSE求解流场的基础上，利用BP-人工神经网络优化方法，展开了寻溯污水系统中致病菌菌源的算法研究。

最优化方法基于单纯的逆问题求解，侧重于单点寻优的最佳预测，对排水管网溯源问题所固有的高度随机性考虑较少，寻优的结果为唯一的取值，一旦寻优路径偏离，结果就不具有实际指导意义。

3）蒙特卡洛随机抽样统计方法

随机抽样的统计方法是基于概率论对所有可能的取值进行随机抽样，遍历所有可能参数的取值，然后把取值作为边界条件和初始条件带入方程1.1求解污染物质的浓度分布，将浓度分布与监测数据进行对比。杨一帆等[46]通过结合污染物对流传输模型，分析了大气应急污染监测的采样数据，采用蒙特卡洛模拟(Monte Carlosimulation)随机生成污染物源项参数，然后通过求解正问题（公式1.1）计算污染物浓度的空间分布，并与实际测量结果进行对比分析，由此对污染源进行定位识别。

随机抽样统计方法概念及计算过程简单，但是计算量极大，而且计算量随着管网节点数的增多呈指数增加，耗费时间长，效率低下反馈慢。通常仅适用于局

8

1 绪 论

部小型管网的溯源。

4）基于后验分布概率的统计方法

由于随机抽样方法效率低下耗时长，在此基础上产生了基于后验分布概率的统计方法抽样方法。其基本的溯源思路与随机抽样类似，也即是对可能的源项参数取值进行随机抽样，然后把取值作为边界条件和初始条件带入方程1.1，求解正问题获得污染物质的浓度分布，将浓度分布与监测数据进行对比。

后验分布统计方法与随机抽样不同的地方在于：1. 以概率密度函数的形式而非一个固定的值来表示求解结果，对抽样得到的取值实际发生的概率进行估计，也即是获得需要反演的参数的后验概率分布，而非单一解，如计算某一排放事故在指定节点发生的概率；2. 统计抽样的过程可融入误差等不确定因素，量化解的不确定性。因此后验分布统计方法对于不确定的反问题显示出较大的优越性。通常后验概率分布不存在明显的函数表达式，而通过贝叶斯理论，可以获得后验分布的取值。

贝叶斯反演方法不是遍历所有的参数，而是对概率较大的区域进行集中抽样，而对于概率较小的区间抽样较少，因此可以灵活高效地实现整个区域内的抽样检查，而且抽取的样本能够代表整体的概率分布。基于贝叶斯统计的概率统计方法的原理是：根据先验知识和正向模型，构建污染源参数的后验概率分布，然后对其进行估计，从而实现对污染源参数的反演。

由于贝叶斯统计反演算法的高效性，尤其是结合马尔科夫蒙特卡罗(MCMC)抽样方法，使得其成为各类溯源问题的有效方法[47]。Zuxin Xu等[48]提出了一种基于标记物种和化学质量平衡（CMB）模型的贝叶斯-蒙特卡罗统计模拟方法，对非暴雨进水情况进行初步评估；曹小群等[49]利用贝叶斯-蒙特卡罗方法求解对流-扩散方程的污染源项识别问题；陈海洋等[50]采用了贝叶斯-蒙特卡罗方法研究水体污染源项识别问题；姜继平等[51]利用贝叶斯理论，结合浓度时间序列方差假定和Adaptive Metropolis MCMC后验采样，建立了用于突发河流水污染应急溯源的贝叶斯推理方法；王晓波等[52]提出了基于事故树（Fault Tree，FT）与动态贝叶斯网络（Dynamic Bayes Network，DBN）的油气管道事故快速辨识及深度溯源分析方法；冯民权等[53]应用贝叶斯网络的逆向推理功能对供水工程水质风险进行计算，分析推理出供水工程水质风险的主要因素及其风险概率。

1.3 污染物统计法数学模型溯源

数学模型的溯源问题本身具有不适定的特性，也即是解的非唯一性，排水管道中，这种不确定性尤其明显。排水管道中，多个排放点浓度存在着空间差异，加上水质特征指标浓度的测量误差，源浓度测量和末端排放口浓度观测存在的非

9

重庆大学硕士学位论文

同步性等特征，更加增加了其中进行污染物质溯源的不确定性。

直接求解和最优化方法反演算法基于单纯的逆问题求解，侧重于单点寻优的最佳预测，对反问题所固有的随机性考虑较少，通常只能获得参数的“点估计”。对于环境水力学反问题中的非唯一问题，往往需要寻求的是所有的可行性解而不是“最优解”，此时最优化方法往往难以实现。目前，学术界普遍认为，对于求解不确定性较强的反演问题，获取一定置信度下的参数分布区间，比找出最优单点更具有实际物理意义。

贝叶斯方法的基本理念和正则化、最优化方法截然不同，其本质不是寻找某一个使得后验概率最高的“最优点”，而是在整个解空间内搜索，得出所有值的概率分布。因此，贝叶斯反演统计法能够有效解决反问题的不确定性这一难题，成为国内外环境流体力学研究的热点，集中在给水管网、地下水及天然河道污染事件溯源中，而类似的反演方法在多节点水力条件复杂的排水管道系统中的研究国内外尚未见系统地开展。因此将贝叶斯统计法推广应用到排水管网污染物溯源求解中，构建一个适用于城市排水管道的源项反演数值模型，研究多节点复杂水力条件下的可靠性、收敛性和稳定性，揭示相关的影响因子，具有重要的意义。

因此，本研究采基于EPA SWMM 水质模型，模拟整个污水管道系统中偷排漏排污染物质的迁移转换过程，求解对流传输方程正问题，然后通过SWMM水质模型与MATLAB软件的耦合，构建污染物质溯源模型，根据水质模型结果与在线监测污染物浓度值的相似程度，应用贝叶斯统计法对排水管网中可能的污染物排放源项进行MCMC抽样，从而有效地获得概率较高的污水排放节点、排放量和排放时间的集合，并开展数值试验，来检验基于贝叶斯统计MCMC抽样反演溯源模型的有效性和准确性。

1.4 研究目的、内容及技术路线

1.4.1 研究目的

针对现有的城市排水系统突发水污染事故溯源技术仍以物理方式人工现场搜索为主，面临着效率较低、响应速度慢及监测点数据利用不充分等局限性，本研究基于EPA SWMM的水力水质模块，集成MATLAB编程，耦合贝叶斯统计反演算法，并采用先进的MCMC抽样进行贝叶斯推理计算，构建排水管道超标污染物质的溯源模型，模拟污水管道系统中保守及非保守污染物质的排放、迁移、转输全过程，对排水管网中可能的偷排漏排污染点源进行统计意义上的反推，从而有效地筛选出可能的污水排放节点、排放量和排放时间，并给出各个可能取值区间的概率;重点研究统计反演算法中游走步长的取值、在线监测点的空间布置及采样监测的时间间隔等限制条件对反演结果精度的影响，确定最佳条件以提高溯源算

10

1 绪 论

法对污水管网溯源结果的有效性和准确度，以期为水体污染控制、智慧排水监管系统等领域提供技术支持。

1.4.2 研究内容

本课题的研究内容主要包括以下四个部分：

① SWMM模型、MATLAB软件及贝叶斯方法的耦合集成

通过MATLAB编程，调用统计函数实现贝叶斯推理计算，并将贝叶斯推理的结果输出为可读文件供SWMM调用，同时通过调用相关函数及指定参数，在MATLAB内部直接运行SWMM，并将SWMM模拟结果的时间序列输出反馈回至贝叶斯分析，实现在不启动SWMM界面前提下，将SWMM模拟结果的动态显示并反馈。结合贝叶斯统计MCMC抽样分析，从而有效地筛选出可能的污染物排放节点、排放量和排放时间，减少大量繁琐的物理溯源所需的人力、物力、财力。

② 贝叶斯统计MCMC抽样算法中游走步长对反演效率和精度的关联关系研究

MCMC抽样算法游走步长的选取直接影响着贝叶斯方法的效率，最佳取值与问题本身的特征有关。研究考察采用不同的游走步长对对反演精度和效率的影响，并确定在本研究中各反演参数抽样时最佳的游走步长。

③ 采样监测的时间间隔及在线监测点的空间布置对反演算法有效性及精度的影响研究

在最佳游走步长确定的基础上，考察瞬时监测数据（仅一个时间节点）以及连续取样监测（多个时间节点）两种监测模式的数据资料对贝叶斯反演算法的精度及效率的影响，并在此基础上，考察不同监测取样的时间间隔对反演精度和效率的影响，并确定在本研究中反演算法最佳的采样监测的时间间隔；增加水质监测点，考察在线监测点的空间布置对反演精度和效率的影响。

④ 反演参数的个数对贝叶斯算法可靠性的影响研究

考察在城市排水管道中常见的两种污染物质排放模式：1）节点瞬时排放；2）节点连续排放，在此基础上进行源项参数识别研究。

需要反演的污染源源项S(x,t,M)包括排放节点x、排放量M和排放时间t三个特征排放参数。这3个特征参数在实际情况中可能部分参数已知，或者3个参数均未知。为考察贝叶斯反演算法对于信息缺失度的敏感程度，在每种排放模式下，考察3种反演的工况：单个未知参数，两个未知参数和三个未知参数，分别代表排放节点、排放量、排放时间三个排放特征参数中有单个、两个、三个未知参数需要反推，其他参数已知的情形，比较分析反演参数的个数对贝叶斯算法有效性的影响。

11

重庆大学硕士学位论文

1.4.3 技术路线

本课题的技术路线如图1.3所示： 阅读相关文献

SWMM水力水质模型构建 Matlab贝叶斯MCMC抽样编程 SWMM模型与MATLAB耦合集成 不同污染物质排放模式情形下溯源 污染物瞬时排放 污染物连续排放 不同步长对反监测取样时间间隔及监测反演信息缺演算法精度的点的空间布置对反演算法失程度对反影响 精度的影响 演精度的影响 污染物SWMM-贝叶斯统计溯源模型 图1.3 课题技术路线图 Fig1.3 Schematic of Technical Solution

12

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建

排水管道系统由于其拓扑结构及水力条件复杂多变，其中的污染物溯源问题具有高度不确定性，而贝叶斯方法是一种随机性反演算法，其中已知信息和反推信息皆用概率分布表示，通过引入测量误差和算法误差直接量化解的不确定性，结果提供所有节点可能为污染源的概率分布，能够较好地应对溯源问题的高度不确定性。

贝叶斯方法将参数定义为随机变量，该变量的分布特征通过经验值取值范围或者先验分布可以确定，将先验分布和传统的似然函数结合得到参数的后验分布，基于此进行统计推理。贝叶斯方法融合了先验信息（在获得监测数据之前具有的对未知源项参数分布的信息，如：根据排放历史某些工业点具有更高的偷排概率）、样本（模拟数据）和总体（监测数据）三种信息，其本质就是“由果求因”，接近人类利用新信息不断调整反推未知的思维过程。

因此，本研究从统计反演的角度出发，采用贝叶斯推理构建了污染物质溯源数学模型。

2.1 单纯贝叶斯反演模型构建

2.1.1 贝叶斯理论

贝叶斯推理的原理是基于条件概率，将未知参数的先验信息转化为后验概率分布的方法。先验信息通过概率密度函数p(X)的形式表示，表示未获得观测值y前参数X的分布规律。这个概率密度函数通常可以通过参数可能的取值范围区间估计得到，代表对未知参数的主观判断。后验分布p(X|y)表示在获得了观测值y之后，参数X的分布规律。贝叶斯的反推过程是基于公式2.1，计算似然函数，获得参数的后验概率分布。

p(X|y)=

其中：X------模型参数；

y------观测数据；

p(X|y)------参数的后验概率密度函数； p(X)------参数的先验概率密度函数； p(y|X)------条件概率密度； p(y)------积分常数。

p(y|X)通常也称为似然函数，表示模型参数与观测数据的相似程度，似然函数值越大表示二者拟合程度越好，反之越差。

13

p(X)p(y|X)

p(y)

∝p(X)p(y|X) （2.1）

重庆大学硕士学位论文

2.1.2 先验分布函数定义

先验分布表述研究者在对数据进行分析之前已有的信息，是指用来描述随机变量位置状况的概率分布，先验信息是与检测结果无关的信息，代表参数的取值范围，通常来源于专家经验和主观判断。通常得到的先验信息是粗糙的，为了使原始的先验信息被贝叶斯统计方法所用，必须将其表达为概率密度函数的形式[54]。 先验分布的形式与取值不需要有严谨的推断过程，通常凭借主观判断或者经验参数取值范围来设定[55]~[58]。常见用于贝叶斯推理反算中的先验概率分布函数有均匀分布、正态分布、对数正态分布等。

在排水管网污染物质溯源问题中，源项的未知参数都分布在一个特定范围内，通常认为参数的取值概率满足这个取值区间之内的均匀分布。因此，当排放源参数的上下界可知时，可以利用区间上的均匀分布来表示。排水管道污染物溯源中，待反演的排放节点、排放量和排放时间三个参数的先验分布函数按照公式2.2 所示均匀分布进行定义。N代表所有需要反推的参数的个数，先验分布的总概率密度函数p(X) 为单个概率密度函数p(Xi)的积。

1

p(X)=

∏ni=1p(Xi) ， p(Xi)={

bi−ai

, Xi∈[ai,bi];

0, 其它。

（2.2）

式中p(Xi)代表取值区间[a，b]上的均匀分布。a、b分别为参数取值的上限和下限。i表示第i个参数。在溯源时，根据实际情况对未知参数的上下限进行设置。

2.1.3 似然函数构造

在基于马尔科夫蒙特卡罗（MCMC）方法的污水管网污染物溯源反演算法中，似然函数表示模型参数和观测数据的拟合程度，因此，似然函数的构造在贝叶斯反演计算中至关重要，它直接影响到污染源溯源统计反演结果的稳定性、可靠性和计算效率。

在排水管道溯源问题中，似然函数p(y|X)通过计算污染物扩散模型模拟的理论值T={T1，T2，⋯，Ti，⋯，Tn}与观测数据Y={Y1，Y2，⋯，Yi，⋯，Yn}的之差的模数来表征模拟值与实测值之间的拟合程度。其中，Ti为通过扩散模型模拟得到的第i个节点的污染物浓度理论值，Yi为污染事件发生后在监测点第i个节点观测到的污染物浓度值。

贝叶斯推理过程考虑了监测过程以及模拟过程的误差因素等数据的不确定性，并将这些不确定性综合考虑，用概率分布的形式表达。常见的造成误差的原因有水质扩散模型的预测偏差和传感器的测量误差。由于任何电子设备都会产生噪声，噪声会干扰甚至淹没信号，特别是当信号十分微弱的情况下[59]~[61]。因此，在污染物浓度检测过程中噪声误差是不能被忽略的。如果反推结果获得排放强度为M，排放节点为Jx，排放时间为t的情况下的误差为 εi=Ti(M，Jx，t|X)−Yi。

14

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建

根据误差理论，测量误差以及模型误差通常服从均值为零、标准偏差为σ的正态分布。因此，综合了误差概率分布的的似然函数可定义为：

p(y|X)=

(Ni(M，Jx，T|X)−Yi)n

∑exp([−i=1(2πσ2)n/22σ2

1

2

]) （2.3）

式中：n为测量数据个数。

2.1.4 后验概率密度函数构造

将公式（2.2）、（2.3）代入（2.1），可以将溯源问题转换成求解未知参数X的后验概率密度函数的问题：

p(X|y)=

∏ni=1p(Xi)

×

(Yi−Ni(M，Jx，T|X))n

∑exp([−i=1(2πσ2)n/22σ2

1

2

]) （2.4）

从（2.4）可以看出，若已知污染源未知参数X的先验分布p(X)，且已知监测点的观测数据Y={Y1，Y2，⋯，Yi，⋯，Yn}，就可以求得未知参数X的后验概率密度，这个概率密度代表了不同X的取值的可能性，也即是未知参数X统计意义下的解。

2.1.5 贝叶斯反演过程

贝叶斯反演实现的具体过程如图2.1所示。

气象数据 X1 参数的先验分布P（X） 扩散模型 先验抽样 更新 观测浓度C 计算浓度C’ 后验抽样（MCMC） 收敛 似然函数P（C|X1） 稳定的后验分布P（X|C）

图2.1 贝叶斯反演结构图

Fig2.1 Schematic of Beyesian Inference

2.2 MCMC抽样Matlab编程

2.2.1 MCMC抽样算法

在贝叶斯统计理论中，后验概率密度函数就是反问题的解，但直接求解后验概率密度函数公式2.4比较复杂。另外，随着未知变量个数n的增加，公式2.4的计算量呈指数级增加。因此，对公式2.4 的求解通常通过统计抽样的方法来实现。

15

重庆大学硕士学位论文

贝叶斯后验概率密度函数的统计抽样，其基本思想是把一个复杂的后验密度函数求解问题转化为一系列简单的抽样问题。抽样过程从任意一个污染物的排放量、排放节点和排放时间初始解开始，通过循环迭代不断调节抽样的范围，使得计算结果逼近实际监测数据值，最终获得未知参数在不同取值下的概率分布，也即是参数的后验分布，来估计污染物质的排放量、排放节点和排放时间三未知个参数。

马尔科夫链蒙特卡罗法（MCMC），是近期发展起来、可很好地用来求解贝叶斯后验概率密度的统计方法．其基本思想是对无明确数学表达式的概率分布P进行抽样，产生大量服从分布P的随机向量序列，通过建立合适的抽样算法使得该序列满足马尔科夫链性质，以保证新候选参数仅依赖于目前的参数，而与之前的候选参数无关。可见，MCMC方法利用马尔科夫链的特性，对所有的解进行搜索，能够保证马尔科夫链花更多的时间在最重要的区域，产生的样本更能够模仿目标分布样本。

在贝叶斯溯源算法中，概率分布P即是根据似然函数构造获得的后验分布。马尔科夫蒙特卡罗法（MCMC）对污染源未知参数X在先验信息范围内根据先验分布进行抽样，抽样时将污染物扩散模型模拟的污染物浓度理论值N={N1，N2，⋯，Ni，⋯，Nn}，接近监测点观测到的污染物浓度观测值Y={Y1，Y2，⋯，Yi，⋯，Yn}的未知参数组N(M，Jx，T|X)保留下来，偏差过大的去除，经过多次迭代得到的抽样结果就是未知的污染源参数信息。

2.2.2 MCMC Matlab编程

对未知参数进行遴选以及对遴选参数值进行取舍的过程也称为向候选值的转移。学者们研究了不同的转移算法来提高MCMC抽样的收敛速度，Metropolis-Hastings 算法是其中较为先进的一种。Metropolis等（1953）最先提出基于马尔科夫链抽样的Metropolis算法，后来Hastings（1970）给出了一般形式的MCMC方法即Metropolis-Hastings算法。

Metropolis-Hastings方法确定污染物的排放量、排放节点和排放时间三个未知参数，由第t次抽样的样本Xt转化到第t+1次抽样的样本Xt+1，最终生成一个马尔科夫链{Xt|t=0，1，2，⋯}，它的平稳分布即为未知参数的后验概率密度函数分布。在所有的Metropolis-Hastings抽样方法中，由建议分布q(∙|Xt)生成候选值X∗，候选值X∗即是污染源未知参数在第t+1次抽样中可能被接受的样本。因为从后验分布中直接抽取样本非常困难，通常抽样是通过选取建议分布来获得后验概率密度函数的带权子样。建议分布的选取直接影响着候选值X∗是否优于Xt、是否能被接受等性能，建议分布的定义域必须覆盖所有的后验分布这是建议分布选择的准则。如果候选值被接受，则在第t+1抽样的样本值为X∗，即Xt+1=X∗；否则在t+1次抽样的样本链停留在原来的状态Xt，即Xt+1=Xt。从当前位置Xt获得候选值X∗，

16

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建

Markov链从Xt位置移动到X∗的接受概率为[62]：

A(Xt，X)=

∗

p(X∗)q(Xt|X∗)

min{1，p(X)q(X∗|X)}tt

=

p(X∗)

min{1，p(X)}

t

（2.5）

式中：p(X∗)，p(Xt)为目标概率密度函数，在贝叶斯推理中为条件概率密度。 MCMC算法的反演过程分以下五步，通过Matlab编程实现： ① 在模型参数先验信息范围内随机参生模型参数初始点Xt；

② 利用水质模型计算出模型参数初始点对应的污染物浓度值，并计算出该 模型参数对应的条件概率密度p(Xt)；

③ 根据建议分布在当前参数Xt状态下抽取新的测试参数X∗，利用水质模型计算出测试参数X∗对应的污染物浓度值及条件概率密度p(X∗)； ④ 从均匀分布U[0，1]中抽取一个随机数u，如果u<

p(X∗)

min{1，p(X)}，则接受

t

该测试参数并设定为当前模型参数，即Xt+1=X∗，否则不接受该测试参数，即Xt+1=Xt；

⑤ 重复步骤③，④直至迭代完成。

通过利用该方法对后验概率密度函数抽样，抽样的结果就是污染源未知参数的估计值，即是所求污染物溯源的结果。MCMC抽样算法的流程如图2.2所示：

17

重庆大学硕士学位论文

开始 模型参数初始点Xt的生成 参数初始点的污染物浓度值及条件概率密度p(Xt)的计算 测试参数X∗对应的污染物浓度值及条件概率密度p(X∗)的计算 接受概率a的计算 a>u 否 是 Xt+1=X∗ Xt+1=Xt 否 迭代完成 是 结束

图2.2 MCMC抽样编程流程图

Fig2.2 Flow chart of MCMC sampling algorithm

2.3 贝叶斯- MCMC反演模型构建

2.3.1 贝叶斯-MCMC反演模型Matlab编程

根据贝叶斯理论及MCMC抽样算法步骤，在Matlab编程环境下进行编程，

18

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建

采用了Matlab里面的unifrnd、normrnd、linspace和hist等主要统计函数实现了抽样中初始值的随机抽取、建议分布随机数的生成以及后验概率密度直方图的生成。各个函数的统计分析功能如下：

X=unifrnd(a，b)：生成上下端点为[a，b]的连续均匀分布的随机数组X； X＝normrnd(μ，sigma)：生成服从正态分布(μ参数代表均值，sigma参数代表方差)的随机数；

X=logspace(a，b，n)：生成一个(1xn)数组，数据的第一个元素值为a，最后一个元素为b，n是总采样点数；

Counts=hist(a，abins)中abins：间隔数，也就是参数a统计的间隔数。 假设目标分布为p(x)，采用MCMC算法对其进行抽样。当proposal分布取正态分布q(x∗|x(i))=N（x(i)，σ2），在给定链长下，部分算法如图2.3所示：

19

重庆大学硕士学位论文

clcclear% Initialize the Metropolis samplerT = 50000; % Set the maximum number of iterationssigma =200; % Set standard deviation of normal proposal densityXmin = -10; Xmax = 20; % define a range for starting valuesX = zeros( 1 , T ); % Init storage space for our samplesseed=1; rand( 'state' , seed ); randn('state',seed ); % set the random seedX(1) = unifrnd( Xmin , Xmax ) % Generate start value% Start samplingt = 1;while t < T % Iterate until we have T samplest = t + 1;% Propose a new value for theta using a normal proposal densityX_star = normrnd(X(t), sigma )Y_Star=0.3*exp(-0.2*X_star*X_star)+0.7*exp(-0.2*(X_star-10.0)*(X_star-10.0));Y_n=0.3*exp(-0.2*X(t-1)*X(t-1))+0.7*exp(-0.2*(X(t-1)-10.0)*(X(t-1)-10.0));alpha = min([1 Y_Star/Y_n]);% Draw a uniform deviate from [ 0 1 ]u = rand;% Do we accept this proposal?if u < alphaX(t) = X_star; % If so, proposal becomes new stateelse X(t) = X(t-1);% If not, copy old stateendend% Display histogram of our samplesfigure( 1 ); clf;subplot( 3,1,1 );nbins = 200;Xbins = linspace( Xmin , Xmax , nbins );counts = hist( X , Xbins );bar( Xbins , counts/sum(counts) , 'k' );xlim( [ Xmin Xmax ] );ylim( [ 0 0.04 ] );xlabel( 'X' ); ylabel( 'p(X)' );% Overlay the theoretical densityy = 0.3*exp(-0.2*Xbins.*Xbins)+0.7*exp(-0.2*(Xbins-10.0).*(Xbins-10.0));hold on;plot( Xbins , y/sum(y) , 'r--' , 'LineWidth' , 3 ); 图2.3 MCMC主要算法代码

Fig2.3 MCMC code

其中：T代表抽样次数即Markov链长，sigma（σ）代表抽样的游走步长。

20

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建

2.3.2 贝叶斯-MCMC反演模型数值验证

在Matlab编程实现MCMC M-H算法的基础上，应用实例案例来验证MCMC M-H算法，并研究不同参数取值对MCMC抽样收敛性的影响。 ① 数值试验设计

实例案例具有目标分布：p(x)∝0.5exp(−0.3x2)+exp（−0.4(x−10)2，进行MCMC算法抽样。在此案例中，这个问题的真实解的概率分布也即是等于0.5exp(−0.3x2)+exp（−0.4(x−10)2。通过Matlab生成的抽样结果概率分布的直方图与实际的概率分布函数进行对比，可以看出反演模型是否真实表达了贝叶斯的反演过程。

在此算例中，建议分布取正态分布q(x∗|x(i))=N（x(i)，σ2）。Markov链长分别取500、5000、50000，σ即游走步长取20进行贝叶斯MCMC反演．部分代码如下所示。

② 数值试验验证结果

抽样结果如图2.3-2.15所示。图中光滑曲线为理论的概率分布函数曲线p（x），条形图代表贝叶斯反演结果的概率分布曲线。两者越接近，则说明反演的精度越高。

由图2.15可以看出，随着循环次数的增加，抽样结果逐渐逼近目标概率分布，也即是说明贝叶斯反推的结果与真实值吻合。根据数值试验理论，说明贝叶斯MCMC M-H 反演模型能够精确反演出真实的结果。贝叶斯MCMC反演模型得到验证。

2.3.3 游走步长及马氏链长对模型性能的影响

为了研究链长及游走步长对贝叶斯反演性能的影响，分别取500、5000、50000三组Markov链长，以及σ=1、5、10、15、20、200共6组游走步长，进行贝叶斯MCMC反演计算。

① 不同工况反演结果

21

重庆大学硕士学位论文

图2.4 MCMC算法抽样直方图（σ=1，链长500）

Fig2.4 Sampling histogram of MCMC algorithm (sigma =1, chain length 500)

图2.5 MCMC算法抽样直方图（σ=1，链长5000）

Fig2.5 Sampling histogram of MCMC algorithm (sigma =1, chain length 5000)

从图2.4与2.5可以看出，σ=1，MCMC链长500时，搜索一直在局部解附近游走，搜索结果陷于局部解，不能正确反演目标概率分布；当MCMC链长增加到5000时，搜索仍然陷于局部解附近，抽样不能正确反演目标概率分布。出现这种情况的原因可能是σ过小，使搜索无法在整个后验空间上抽样，陷于局部解，无法反演真实的后验分布。所以，在此基础上加大游走步长σ，应该会避免局部解的情况。

图2.6 MCMC算法抽样直方图（σ=5，链长500）

Fig2.6 Sampling histogram of MCMC algorithm (sigma =5, chain length 500)

22

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建

图2.7 MCMC算法抽样直方图（σ=5，链长5000）

Fig2.7 Sampling histogram of MCMC algorithm (sigma =5, chain length 5000)

从图2.6与2.7中可以看出，当σ取值增大到5，MCMC链长500 时，搜索范围增大，但抽样的结果显示在0附近的概率较大，仍然未能正确反演目标概率分布；当MCMC链长增加到5000时，搜索结果表明抽样的频率仍然在0附近非常密集，这表明在0附近的概率较大，在10右边的取值范围几乎没有抽样发生，抽样不能正确反演目标概率分布。虽然当σ=5时抽样结果仍然没有正确反演目标概率分布，但比σ=1时更加接近真实的后验分布。所以应该继续加大游走步长σ，应该提高后验分布的反演精度。

图2.8 MCMC算法抽样直方图（σ=10，链长500）

Fig2.8 Sampling histogram of MCMC algorithm (sigma =10, chain length 500)

图2.9 MCMC算法抽样直方图（σ=10，链长5000）

Fig2.9 Sampling histogram of MCMC algorithm (sigma =10, chain length 5000)

23

重庆大学硕士学位论文

图2.10 MCMC算法抽样直方图（σ=15，链长500）

Fig2.10 Sampling histogram of MCMC algorithm (sigma =15, chain length 500)

图2.11 MCMC算法抽样直方图（σ=15，链长5000）

Fig2.11 Sampling histogram of MCMC algorithm (sigma =15, chain length 5000)

图2.12 MCMC算法抽样直方图（σ=20，链长500）

Fig2.12 Sampling histogram of MCMC algorithm (sigma =20, chain length 500)

图2.13 MCMC算法抽样直方图（σ=20，链长5000）

Fig2.13 Sampling histogram of MCMC algorithm (sigma =20, chain length 5000)

从图2.8-2.13可以看出，σ取值增大到10、15、20，MCMC链长500时，MCMC搜索范围继续增大，抽样的结果逼近目标概率分布；当MCMC链长增加到5000时，抽样的结果不仅逼近目标概率分布，而且精度更高。从图2.11可以看出，σ=15，

24

2 基于MATLAB编程的贝叶斯-MCMC反演模型构建

MCMC链长增加到5000时，抽样结果与理论值基本接近。

图2.14 MCMC算法抽样直方图（σ=200，链长500）

Fig2.14 Sampling histogram of MCMC algorithm (sigma =200, chain length 500)

图2.15 MCMC算法抽样直方图（σ=200，链长5000）

Fig2.15 Sampling histogram of MCMC algorithm (sigma =200, chain length 5000)

图2.16 MCMC算法抽样直方图（σ=200，链长50000）

Fig2.16 Sampling histogram of MCMC algorithm (sigma =200, chain length 50000)

从图2.14-2.16可以看出，当σ = 200、MCMC链长为500或5000时，抽样结果精度太低；当MCMC链长增加到50000时，抽样结果才逼近目标概率分布。

② 结果讨论

通过上述分析得到以下结论：

1）建议分布选取正态分布N（x(i)，σ2）进行MCMC抽样时，游走步长σ和MCMC链条的长度对抽样结果影响较大。游走步长过小可能使抽样无法在整个后验空间搜索，使抽样结果陷于局部解，无法准确反演污染源参数的后验概率密度；适当游走步长取值可以实现对污染物反演参数整个后验空间的抽样，避免使反演

25

重庆大学硕士学位论文

参数限于局部解，提高反演效率和反演精度；但是步长过大会使抽样样本跳出后验空间，只有很少的样本落在了后验空间，使抽样结果精度降低，若要在过大步长下提高精度就必须增加MCMC链长，这样就会导致抽样时间延长，抽样效率降低。

2）单纯增大MCMC链条的长度可以增加反演精度，但是对局部解的影响作用不明显，仍然会使抽样结果陷于局部解，而且会使抽样时间延长，从而降低抽样效率。

3）设计不合理的马氏链可能使抽样结果限于局部最小解，且降低反演效率和反演精度。所以在进行MCMC抽样时要根据先验范围选择适当的游走步长，确定合适的MCMC链长长度，便可提高后验概率部分的反演精度和反演效率。

4）采用正态分布作为MCMC抽样的建议分布，可以实现对瞬时排放的稳态污染物质排放地点较精确的反演；

2.4 本章小结

本章基于贝叶斯理论，定义了先验分布函数，构造了似然函数和后验概率密度函数，在对候选参数的抽样过程中，选取正态分布 N（x(i)，σ2）作为建议分布，通过Matlab编程构建了贝叶斯MCMC M-H反演模型。然后通过数学模型试验，在理论解已知的条件下，对构建的贝叶斯反演模型进行了验证。

在此基础上，为了考察不同的数值参数对贝叶斯反演模型效率的影响， 采用了5组游走步长和3组循环链长，对贝叶斯反演模型进行了收敛性分析。结论如下：

① 构建的贝叶斯反演模型很好地实现了在给定的后验分布概率密度函数下的抽样，随着循环次数的增加，反演结果逐渐向理论值逼近，最后收敛到理论解附近。

② 步长σ和MCMC链条的长度对抽样结果影响较大。大的步长取值可以游走更大的搜索空间，避免限于局部解，可以实现对整个后验空间的抽样，但是方差过大会使抽样次数增加，反演时间延长，效率降低。

③ 增大MCMC链条的长度可以增加精度，但是对局部解的影响作用不明显；设计不合理的马氏链可能使抽样结果限于局部最小解导致结果不能收敛到真实值。

贝叶斯反演利用了先验信息，反演过程考虑了观测数据误差的影响，结果采用概率分布的表达形式，能够较好求解排水管网溯源这一高度不确定环境反问题。

26

3 SWMM、贝叶斯与MATLAB耦合的污染物溯源模型

3 SWMM、贝叶斯与MATLAB耦合的污染物溯源模型

SWMM水质模块求解一维污染物对流扩散模型，能够较精确地动态模拟排水管道系统中保守物质及非保守物质的迁移与转化过程。且SWMM可将水质模块与动态波水力计算模块耦合，实现包括压力流、逆流、溢流、泵站提升等多种复杂水力条件下的污染物传输模拟，与实际管网的动态运行工况更接近。MATLAB具有强大的统计分析计算功能，能够实现贝叶斯反演算法涉及的大量抽样计算，提高了统计反演法的效率。

因此，本研究耦合集成具有强大的数据处理功能Matlab编程软件和SWMM水质模型，构建了基于贝叶斯理论的SWMM溯源模型。通过MATLAB编程，调用统计函数实现贝叶斯推理计算，并将贝叶斯推理的结果输出为可读文件供SWMM调用；然后通过调用相关函数及指定参数，在MATLAB内部直接运行SWMM，并将SWMM模拟结果的时间序列输出反馈回至贝叶斯分析，从而有效地筛选出可能的污水直排节点、排放量和排放时间。

由于贝叶斯反演过程融入观测数据误差的影响，结果采用概率分布的表达形式，在分析反算数据的不确定性上具有独特优势，在实际应用中提供更多的未知污染源的信息，为物理搜索提供更明确高效的搜索路径， 不仅减少大量繁琐的人工识别所需的工作量，且提高溯源效率，实现更快速的反馈，起到保护水体环境的作用。

3.1 SWMM水质模型

本研究基于SWMM模型，在排水管网系统中构建一维水质模拟模型，求解一维污染物对流扩散方程3.1：

∂(AC)∂t

+

∂(QC)∂x

=∂x(AEx∂x)−AK1C+S(x,t,M） （3.1）

∂∂C

其中：A------过水断面面积，（m2）； Q------流量，（m3/h）；

C------污水管网中某种污染物的浓度，（mg/L）； t-------时间，（s）； x-------流程，（m）；

Ex------纵向离散系数，（m2/s）； K1------污染物一级降解系数，（s-1）；

S（x,t,M）------ 污染点源排放随时间变化的函数。

方程3.1是一个偏微分方程，通常采用数值的方法进行求解。求解该方程的前

27

重庆大学硕士学位论文

提是流场已知，也即是流量或者流速等参数已经确定。排水管道水质模型通常综合水文、水力的计算来获得管道沿程的水量和流速分布，然后再数值离散求解方程3.1，是一个复杂耦合的过程。SWMM包含水文、水力、水质3大计算模块，可以实现之间的耦合。SWMM模型架构如图3.1所示。

数据输入 全局对象数据库.INP 水文、水力、水质模拟 .OUT 路径/径流解算 图形用户界面 报告输出 图3.1 SWMM模型架构 Fig3.1 SWMM structure

.RPT

3.2 SWMM模型集成技术难点

SWMM引擎对数据的输入和输出提了严格的要求，SWMM的输入文件.inp为一文本文件，包含所有SWMM模型运行所需的各个参数的取值以及建模的空间数据。SWMM运行后自动生成分析报告文件*.rpt和运行结果文件*.out。

分析报告文件*.rpt是对整个模拟过程的总结，不包括模拟过程中各个对象各个时刻的信息记录，如降雨、管道流量、节点水位等，以文本形式保存，可以被直接提取。然而，SWMM模型在整个模拟过程中各个对象的时间序列，也即是整个监测过程的浓度变化曲线，以二进制方式保存在运行结果文件*.out中。该文件中的信息无法被直接提取，因此集成所面临的问题就是将SWMM模拟过程中的时间序列输出，能够直接被调用，用于进一步分析。

对SWMM模型进行外部调用时，通常需要读取到SWMM 模型运行输出的结果文件中的内容，传统的方法有两种：1.修改源代码，将模拟过程的时间序列输出，重新编译修改后的SWMM源代码生成可执行文件，然后才能实现模拟结果的动态显示；2.运行 SWMM可执行文件后，采用SWMM核心代码中提供的基于VB、C/C++和Delphi Pascal 语言平台的接口函数，逐层读取结构复杂的二进制模拟结果文件。

本研究通过直接在Matlab下对SWMM可执行文件采用命令，在不启动

28

3 SWMM、贝叶斯与MATLAB耦合的污染物溯源模型

SWMM界面前提下，将参数传输给SWMM水质模型并运行再输出结果，实现将SWMM模拟结果的动态显示及反馈，模拟时间序列数据可以进一步被调用。本集成方式操作简单，克服了以往实现SWMM模拟结果动态显示的工作量大且耗时的缺点，该技术已经申请专利“一种识别雨水管网污水直排污染源的方法”（201810114905.X）。

3.3 SWMM、贝叶斯及MATLAB耦合污染物溯源模型

以MATLAB为主要界面，将SWMM水质模型、MATLAB及贝叶斯模型进行集成。集成方法步骤如下：

步骤1：在SWMM软件中输入构建SWMM模型所需的基础数据，生成运行SWMM软件所需的初始输入*.inp文件；

步骤2：在Matlab中打开初始*.inp文件，对*.inp文件进行读写，然后对所以可能污水直排的节点、排放量和排放时间三个参数进行定位和修改；

步骤3：在Matlab调用swmm5.exe程序，利用更新后的*.inp文件中的信息进行计算，并将运算结果储存在*.rpt文件中；

步骤4：在Matlab中实现对*.rpt文件的读写，然后提取*.rpt文件中相应的模拟结果的时间序列数据样本，结合实际测量时间序列数据样本，进行MCMC抽样分析；

步骤5：在Matlab中用MCMC算法对三个参数在先验范围内进行抽样分析，每次抽样后，均需更新*.inp文件中相应位置的参数信息，并将更新的运算结果储存在*.rpt文件中，通过步骤4提取的数据进行计算，直至迭代完成。

水质模型与MATLAB的集成流程图如下：

29

重庆大学硕士学位论文

开始 构建SWMM模型，生成初始输入*.inp文件 Matlab中打开初始*.inp文件 读写*.inp文件 反演参数的定位和修改 调用Matlab外部的swmm5.exe程序运行更新的*.inp文件并储存结果于*.rpt文件 Matlab中*.rpt文件的读写 提取时间序列数据，进行MCMC抽样 直至迭代完成更新*.inp文件中相应位置的参数信息，并将更新的运算结果储存在*.rpt文件中 结束 图3.3 SWMM水质模型、MATLAB及贝叶斯的集成流程图 Fig3.3 Integrated flow chart of SWMM、MATLAB and Bayesian

3.4 贝叶斯-SWMM污染物溯源模型的特点

耦合集成SWMM模型、MATLAB、贝叶斯理论的溯源模型具有以下几个特点： ① 将污水管网系统中污水偷排漏排污染源识别问题转化为数学问题，对所有可能节点、排放量和排放时间三个参数进行水质模型预测，利用

30

3 SWMM、贝叶斯与MATLAB耦合的污染物溯源模型

Matlab软件不断读取不同源项取值下SWMM更新后的污染物浓度随时间变化的过程线，可快速预测不同排放情景，提高物理溯源效率。

② 利用了SWMM强大的水文、水力、水质计算功能，能够较精确地反映排水管道系统中保守物质及非保守物质的迁移与转化过程，且可将水质模块与动态波水力计算模块耦合，实现包括压力流、逆流、溢流、泵站提升等多种复杂水力条件下的污染物传输模拟，与实际管网的动态运行工况更接近，提高了分析的精度。

③ 采用了贝叶斯方法进行统计反演，充分考虑了先验信息、观测噪声以及模型误差的影响，结果采用概率密度的分布形式，能够获得实际意义的结果，在高度不确定排水管网系统的溯源反问题求解中独具优势。

④ 利用Matlab强大的数据处理功能，能够实现MCMC抽样法涉及的大量计算，快速、准确分析出各参数组合的似然度值，从而快速确定出污水直排可能的节点、排放浓度和排放量。

3.5 本章小结

本章开展了SWMM水质模型、贝叶斯统计推理及MATLAB软件的耦合集成研究，并基于SWMM的水力水质模块，集成MATLAB编程，耦合贝叶斯统计反演算法，并采用先进的MCMC抽样进行贝叶斯推理计算，构建了排水管道超标污染物质的溯源模型。溯源模型可模拟污水管道系统中保守及非保守污染物质的排放、迁移、转输全过程，并基于贝叶斯理论对排水管网中可能的偷排漏排污染点源进行统计意义上的反推并给出各个可能取值区间的概率，从而筛选出可能的污水排放节点、排放量和排放时间。

通过贝叶斯方法、Matlab与SWMM的耦合集成，水质预测结果与实际管网的动态运行工况更接近，提高污染物传输模拟的精度，从而提高反演分析的精度，同时贝叶斯统计反演，充分考虑了先验信息、观测噪声以及模型误差的影响，能较好应对排水管网系统溯源反问题的不确定性。

31

重庆大学硕士学位论文

32

4 排水管道污染物溯源模型的验证与应用

4 排水管道污染物溯源模型的验证与应用

本章中采用实际管网对构建的贝叶斯-SWMM溯源模型进行验证。由于污水管道物理试验受到水力相似、几何相似、和摩阻相似多个方面的约束，构建复杂规模较大，容易受到缩放因子过大不能反映真实的水力条件的影响，且反应的工况单一，而数值试验作为理论、物理、数值等三种最重要的实验方法之一，能够反复进行不同工况模型，方便地控制和调整参数。因此，本研究采用数值试验的方法对构建的污染物质溯源模型进行验证。

遵循数值试验方法常用的数据生成方法，提供真实的污染物质排放特征函数作为初始条件及边界条件，求解污染物质的对流传输方程正问题，以获得监测点处的污染物质浓度时间序列。将正问题求解获得的监测点浓度值进行扰动，作为一个“实际污染物浓度监测数据，将此数据作为反演溯源模型的监测数据。然后将真实的源项参数与反演算法获得的源项参数进行对比，可以评价算法的精确性。

本章主要通过数值试验生成理论解，并进行扰动作为反演的基础数据，研究对污水管网节点瞬时排放污染物及节点连续排放污染物两种排放模式。

需要反演的污染源源项S(x,t,M) 包括排放节点x、排放量M和排放时间t三个特征排放参数。这3个特征参数在实际情况中可能部分参数已知，或者3个参数均未知。为考察贝叶斯反演算法对于信息缺失度的敏感程度，在每种排放模式下，考察3种反演的工况：单个未知参数，两个未知参数和三个未知参数，分别代表排放节点、排放量、排放时间三个排放特征参数中有单个、两个、三个未知参数需要反推，其他参数已知的情形，比较分析反演参数的个数对贝叶斯算法有效性的影响。

4.1 仿真数值试验设计

利用SWMM构建排水管道水质模型，已知排放节点进行正向模拟，生成仿真试验数据的监测数据，对反演模型的有效性及可行性进行验证。

图4.1为某一小型污水管网的平面布置图，管网中水流流向如箭头所示，该污水管网共有管道64段，节点65个，管径为400～800mm，管道形状为圆形，粗糙系数设为0.01，PFK1为该污水管网的总排放口。其中布设废水超标排放点，设置不同的排放函数S(x,t,M)，利用SWMM求解污染物的对流传输过程，模拟污染物的转输过程，并在下游某节点处布置水质监测点，记录污染物从不同节点排放时，在监测点处观测到的污染物浓度时间序列。

将此时间序列进行扰动，作为数值试验中的观测数据，输入溯源模型。

33

重庆大学硕士学位论文

4.1.1 节点瞬时排放事件数值试验设计

① 污染物排放初始条件设置

假定某日1:00时，分别在节点53、40、50处的某些工厂瞬时排入重量为1吨的某可降解污染物质BX，在下游监测点1处观测此污染物。节点53、40、50为恒定进流，流量分别为10m3/s、1m3/s、1m3/s。污染物瞬时排放的时间序列如图4.2，污水管网SWMM模型模拟参数如表4.1。

为了模拟污染物在管道中的转输降解过程，假设BX在管道中的衰减系数为0.25 d-1。

表4.1 污水管网SWMM模型模拟参数

Table4.1 Simulation parameters of SWMM model for sewer network

参数 数值

管径（mm） 400～800

管段形状 CIRCULAR

管段粗糙系数

0.01

污染物衰减系数 模拟时间（h）

0.25

节点40恒定进流量（m3/s）

1

3

节点50恒定进流量（m3/s）

1

参数

BX瞬时排入量

（吨）

数值

1

节点53恒定进流量（m3/s）

10

图4.1 污水管网平面布置图

Fig4.1 Layout of the studied sewage network

34

4 排水管道污染物溯源模型的验证与应用

图4.2 污染物瞬时排放的时间序列

Fig4.2 Time series of instantaneous release of pollutants

② 瞬时污染物排放模拟结果

分别在节点53、40、50瞬时排入污染物BX时运行SWMM，在下游检测点1处观测到的BX浓度序列分别如图4.3、图4.4、图4.5。

图4.3 节点53为瞬时排放节点，检测点1处观测的BX浓度序列

Fig4.3 Node 53 is the instantaneous discharge node, and the concentration sequence at the

detection point 1

图4.4 节点40为瞬时排放节点，检测点1处观测的污染物浓度序列

Fig4.4 Node 40 is the instantaneous discharge node, and the concentration sequence at the

detection point 1

35

重庆大学硕士学位论文

图4.5 节点50为瞬时排放节点，检测点1处观测的污染物浓度序列

Fig4.5 Node 50 is the instantaneous discharge node, and the concentration sequence at the

detection point 1

从图4.3可以看出，节点53为BX瞬时排放节点时，在监测点1处1:05时首次观测到BX，BX浓度最高值发生在1.82h时，达0.275mg/L，到3h时降解到0.01mg/L；从图4.4可以看出，节点40为BX瞬时排放节点时，在监测点1处1:03时首次观测的BX，BX浓度最高值发生在1.6h时，达3.2mg/L，到3h时降解到0.025mg/L；从图4.5可以看出，节点50为BX瞬时排放节点时，在监测点1处观测到浓度最高值发生在1.3h时，达5.5mg/L，到3h时已完全降解。

③ 结果讨论

从模拟结果可以看出，从不同节点瞬时排入相同量的污染物，在下游同一监测点首次观测到污染物的时间不同，污染物浓度达到最高值的时间不同，污染物浓度最高值不同，模拟结束时的浓度也不同。原因是不同节点的流量不同，对污染物的稀释降解程度不同；不同节点距监测点的距离不同，排放点越靠近监测点越早观测到污染物出现。

总之，从不同节点瞬时排入相同量的污染物，通过运行SWMM模型模拟污染物在管道中的转输过程，在下游同一监测点观测到的污染物浓度序列不同。这为本研究通过采用贝叶斯统计MCMC抽样算法反演污染源源项研究方法和研究内容提供理论依据。

4.1.2 节点连续排放事件数值试验设计

① 污染物排放初始条件设置

图4.1为某一小型污水管网的平面布置图，管网中水流流向如箭头所示，该污水管网共有管道64段，节点65个，管径为400～800mm，管道形状选CIRCULAR，粗糙系数设为0.01，其中PFK1为该污水管网的总排放口。

假定某日1:00时，分别在节点53、40、50处的某些工厂以100000g/min的速度连续排入污染物BX，在下游监测点1处观测此污染物。节点53、40、50为恒定进流，分别为10m3/s、1m3/s、1m3/s，为利用SWMM模型模拟污染物在管道

36

4 排水管道污染物溯源模型的验证与应用

中的的转输过程，设置BX在管道中的衰减系数设为0.25，污染物瞬时排放的时间序列如图4.6，污水管网SWMM模型模拟参数如表4.2：

表4.2 污水管网SWMM模型模拟参数

Table4.2 Simulation parameters of SWMM model for sewer network

参数 数值

管径（mm） 400～800

管段形状 CIRCULAR

管段粗糙系数

0.01

污染物衰减系数 模拟时间（h）

0.25

3

节点50恒定进流量（m3/s）

1

参数 BX连续排入速度

（g/min）

节点53恒定进流量

（m3/s） 10

节点40恒定进流量（m3/s）

1

数值 100000

图4.6 污染物连续排放的时间序列

Fig4.6 Time series of continuous discharge of pollutants

② 污染物连续排放模拟结果

分别在节点53、40、50连续排入污染物BX时运行SWMM，在下游检测点1处观测到的BX浓度序列分别如图4.7、图4.8、图4.9。

图4.7 节点53为连续排放节点，检测点1处观测的BX浓度序列

Fig4.7 Node 53 is a continuous discharge node, and the concentration sequence at the detection

point 1

37

重庆大学硕士学位论文

图4.8 节点40为连续排放节点，检测点1处观测的BX浓度序列

Fig4.8 Node 40 is a continuous discharge node, and the concentration sequence at the detection

point 1

图4.9 节点50为连续排放节点，检测点1处观测的BX浓度序列

Fig4.9 Node 50 is a continuous discharge node, and the concentration sequence at the detection

point 1

从图4.7可以看出，节点53为BX连续排放点时，在监测点1处1:06时首次观测的BX，BX浓度在3h时为1.366mg/L；从图4.8可以看出，节点40为BX连续排放节点时，在监测点1处1:03时首次观测的BX，BX浓度在3h时为13.885mg/L；从图4.9可以看出，节点50为BX瞬时排放节点时，在监测点1处1:02时首次观测的BX，BX浓度在3h时为17.346mg/L。

③ 结果讨论

从模拟结果可以看出，从不同节点以相同速度连续排入污染物，在下游同一监测点首次观测到污染物的时间不同，在模拟结束时污染物浓度上升率也不同。

总之，从不同节点以相同速度连续排入污染物，通过运行SWMM模型模拟污染物在管道中的转输过程，在下游同一监测点观测到的污染物浓度序列不同。这为本研究通过应用贝叶斯统计MCMC抽样法反演污染源源项的研究方法和研究内容提供理论依据。

38

4 排水管道污染物溯源模型的验证与应用

4.2 节点瞬时排放污染事件溯源

图4.10为某小型污水管网平面布置图，管网中水流流向如箭头所示，该污水管网共有管道64段，65个节点，管径为400～800mm，管道形状选CIRCULAR，其中PFK1为该污水管网的总排放口，在管网下游主干管上某节点监测处布设水质监测点，假定该污水管网水文条件稳定，满足一维水质模型应用条件。

假定某日1:00时，在节点28处的某工厂瞬时排入重量为1吨的污染物BX，在下游监测点1处1:05时观测到BX污染带前锋，现利用监测节点1处噪声的浓度场部分分布数据反演上游BX排放源的排放强度、排放位置和排放时间，为了利用SWMM模拟污染物的转输过程，设置BX在管道中的衰减系数为0.25，考虑误差因素，贝叶斯反演计算假定监测点1处的污染物浓度序列近似正态分布，从图4.13可看出最高值发生在模拟开始后1.35h，污染物浓度达到0.13mg/L，模拟分析结束时污染物已降低到0.006mg/L。相关参数及先验信息见表4.3。

图4.10 污水管网平面布置图

Fig4.10 Layout of the studied sewage network

39

重庆大学硕士学位论文

表4.3 污水管SWMM模型参数

Table4.3 Parameters of SWMM model for sewage pipe

参数 类型 数值 参数 类型 数值

管径（mm）

已知 400～800 模拟时间（h）

已知 3

管段形状 已知 CIRCULAR

管段粗糙系数

已知 0.01

排放时间T（min）

待反演 60

污染物衰减系数

已知 0.25 排放节点Jx 待反演 J28

排放量 M（g）

待反演 1000000

4.2.1 单个未知参数反演

在三个污染物质排放特征参数（排放节点、排放量、排放时刻）中，假定只有一个参数未知需要通过溯源模型获得信息，而其余两个排放参数已知。在监测点实测污染物浓度数据的基础上，利用溯源模型进行未知参数的反演推理，获得该参数的概率分布。因此包括三个情景的反演计算，分别为排放节点、排放量、排放时刻的反演。

本节在进行贝叶斯MCMC抽样反演分析时采用不同的游走步长，分析不同的游走步长对反演精度和效率的影响，并确定在本研究中各反演参数抽样时最佳的游走步长。

① 未知参数排放节点反演分析 1）排放节点的反演分析

假定某日在1:00时，在节点Jx处的某工厂瞬时排入重量为1吨的污染物BX。在监测点1处首次观测时间为1:10，后每隔10min观测一次，共6个数据，检测时间间隔为10min的观测浓度序列见图4.11。

根据先验信息，待反演模型参数先验分布为均匀分布，其对应的先验概率密度函数分别为：

P(Jx)={

0.0156， 1≪x≪65

0， 其它

（4.1）

由于一般情况观测噪声为白噪声，所以假定观测误差服从标准差为σ=0.01的正态分布N(0，σ2)，则污染物溯源的似然函数可表示为：

p(Jx)=

(Ci(Jx)−Ci)6

∑exp([−i=1(2πσ2)6/22σ2

1

2

2

]) （4.2）

根据贝叶斯定理，模型参数的后验概率密度函数为：

p(Jx)=

式中：α为比例常数。

40

(Ci(Jx)−Ci)

αexp(∑6[−i=12σ2

]) （4.3）

4 排水管道污染物溯源模型的验证与应用

上式即为污水管网污染事件溯源的未知参数Jx的后验概率密度函数，通过MCMC抽样算法对模型进行抽样，就可以得到未知参数Jx的后验概率密度函数分布。

2）排放节点的反演结果分析

因为节点是离散点，所以算法中利用Proposal分布提取新的测试值时，分别以游走步长为1、3等概率随机在当前值得基础上左右抽取，取值范围为先验信息范围Jxϵ[1 65]，按照MCMC反演思路迭代2000次。步长为1、3时反演参数Jx的后验概率直方图分别见图4.12、4.13，参数Jx落在各区域的频数统计见表4.4。

图4.11 检测时间间隔为10min的观测浓度序列

Fig4.11 Detection concentration sequence with a time interval of 10 min

图4.12 步长=1时，参数Jx的后验概率直方图

Fig4.12 Step size=1,the posterior probability histogram of parameter Jx

表4.4 步长=1，参数Jx落在各区域的频数统计

Table4.4 Step size=1,frequency statistics of parameter Jx falling in each region

区域（x）

Jx

频数

1500

500

20～40

40～60

41

重庆大学硕士学位论文

图4.13 步长=3，参数Jx的后验概率直方图

Fig4.13 Step size=3,the posterior probability histogram of parameter Jx

从图4.12可以看出，当步长=1时，污水管道污染物瞬时排放节点Jx的后验分布较离散，这是由于不同节点和管段污染物离散和降解过程不同，所以在某些节点上抽样结果是空白的，但还是可以看出污染源排放节点的后验分布主要集中在Jx=[20,52]，主要分布在预先设定节点28的周围；从图4.13可以看出，当步长=3时，污染物瞬时排放节点Jx的后验分布比较分散，没有在真实解附近抽样，原因是由于游走步长为整数且过大，可能跳过真实解，而使后验分布不收敛。

3）结果讨论

通过上述分析，污水管道污染物瞬时排放的反演参数排放节点Jx，在贝叶斯MCMC抽样算法中游走步长为1时，可以提高反演精度和反演效率。这是因为污水管网污染物排放节点为整数离散点，在MCMC抽样反演中如果以大于1的步长游走搜索，则无法在整个后验空间内搜索，可能跳过真实值，使抽样结果无法反演真实的后验分布。

当游走步长=1时，污染源Jx参数的后验统计结果如图4.14所示：

图4.14 污染源参数Jx的后验统计结果

Fig4.14 The posterior statistical results of pollution source parameters Jx

42

4 排水管道污染物溯源模型的验证与应用

② 未知参数排放量反演分析 1）排放量的反演分析

假定某日在1:00时，在节点28处的某工厂瞬时排入重量为M（g）的污染物BX。在监测点1处首次观测时间为1:10，后每隔10min观测一次，共6个数据。

根据先验信息，待反演模型参数先验分布为均匀分布，其对应的先验概率密度函数分别为：

0.1×10−5， 500000≪M≪1500000

P(M)={ （4.4）

0， 其它

由于一般情况观测噪声为白噪声，所以假定观测误差服从标准差为σ=0.01的正态分布N(0，σ2)，则污染物溯源的似然函数可表示为：

p(M)=(2πσ2)

1

∑66/2exp(i=1[−

(Ci(M)−Ci)2

2σ2

]) （4.5）

根据贝叶斯定理，模型参数的后验概率密度函数为：

p(M)=βexp(∑6i=1[−

式中：β为比例常数。

上式即为污水管网污染事件溯源的未知参数M（g）的后验概率密度函数，通过MCMC抽样算法对模型进行抽样，就可以得到未知参数M（g）的后验概率密度函数分布。

2）排放量的反演结果分析

算法中Proposal为正态分布N(μ，σ^2)，其中σ即游走步长分别取10，10000，10000000，按照MCMC反演思路迭代2000次。σ分别取10，10000，10000000时反演参数M（g）的后验概率直方图分别见图4.15、4.16、4.17，参数M（g）落在各区域的频数统计见表4.5。非最佳游走步长的统计结果不做频数统计。

(Ci(M)−Ci)2

2σ2

]) （4.6）

图4.15 σ=10，参数M（g）的后验概率直方图

Fig4.15 sigma =10，the posterior probability histogram of parameter M（g）

43

重庆大学硕士学位论文

图4.16 σ=10000，参数M（g）的后验概率直方图

Fig4.16 sigma =10000,the posterior probability histogram of parameter M（g）

表4.5 σ=10000，参数M（g）落在各区域的频数统计

Table4.5 sigma =10000,frequency statistics of parameter M（g） falling in each region

区域（×106g）

M（g）

频数

400

1550

50

0.8～0.9

0.9～1.1

1.1～1.2

图4.17 σ=10000000，参数M（g）的后验概率直方图

Fig4.17 sigma =10000000，the posterior probability histogram of parameter M（g）

从图4.15可以看出，当σ=10时，抽样没有在整个后验空间内搜索，抽样结果限于局部解，若增加迭代次数则反演时间延长，反演效率降低；从图4.16可以看出，σ=10000时，参数排放量M（g）的后验概率密度分布非常集中，绝大多数集中在实际值1000000（g）附近，基本符合预先设定值，从图可以推断污染事件的排放量在M=[85000g,105000g]之间；从图4.17可以看出，σ=10000000时，由于游走步长过大使极少部分抽样样本落在后验空间范围内，无法反演参数M(g)的后验概率分布，若要提高抽样结果精度，需增加迭代次数，则反演时间延长，反演效率降低。

3）结果讨论

通过上述分析，污水管道污染物瞬时排放的未知反演参数排放量先验范围为M∈[500000 1500000]时，在贝叶斯MCMC抽样算法中游走步长为10000，可以提高反演精度和反演效率。如果游走步长相对于先验范围过小，则抽样结果会陷于局部解，无法反演真实的后验分布；如果游走步长先对于先验范围过大，则大

44

4 排水管道污染物溯源模型的验证与应用

多数的样本不会落后验空间内，使反演精度和反演效率降低。

如表4.6所示，当游走步长=10000时，通过将污染物排放量的抽样结果与真实值对比，可知污染物排放量的数学期望和中值误差均较小。污染源参数M（g）的反演结果误差柱状图见图4.18。

表4.6 污染源参数M（g）的后验统计结果

Table4.6 The posterior statistical results of pollution source parameters M（g）

参数 M（g）

实际值 1000000

中值 987454

中值误差/%

1.25

均值 975237.34

均值误差/%

2.48

M（g）420中值误差/% 均值误差/%2.481.25 图4.18 污染源参数M（g）的反演结果误差

Fig4.18 Inversion error of pollution source parameter M(g)

③ 未知参数排放时刻反演分析 1）排放时刻的反演分析

假定某日在T（min）时，在节点28处的某工厂瞬时排入重量为1吨的污染物BX。在监测点1处首次观测时间为1时10分，后每隔10min观测一次，共6个数据。

根据先验信息，待反演模型参数先验分布为均匀分布，其对应的先验概率密度函数分别为：

P(T)={

0.0083， 0≪T≪1200， 其它

（4.7）

由于一般情况观测噪声为白噪声，所以假定观测误差服从标准差为σ=0.01的正态分布N(0，σ2)，则污染物溯源的似然函数可表示为：

p(T)=(2πσ2)6/2exp(∑6i=1[−p(T)=γexp(∑6i=1[−

式中：γ为比例常数。

1

(Ci(T)−Ci)2

2σ2

]) （4.8）

根据贝叶斯定理，模型参数的后验概率密度函数为：

(Ci(T)−Ci)2

2σ2

]) （4.9）

45

重庆大学硕士学位论文

上式即为排水管网污染事件溯源的未知参数T（min）的后验概率密度函数，通过MCMC抽样算法对模型进行抽样，就可以得到未知参数T（min）的后验概率密度函数分布。

2）排放时刻的反演结果分析

由于在模拟中时间序列是以分钟为最小单位且为整数，所以算法中利用Proposal分布提取新的测试值时，分别以步长为1min、3min等概率随机在当前值得基础上左右提取，取值范围为先验信息范围Tϵ[0 120]，按照MCMC反演思路迭代2000次。步长为1min、3min时参数T（min）的后验概率直方图分别见图4.19、4.20，参数T（min）落在各区域的频数统计见表4.7。

图4.19 步长=1min，参数T（min）的后验概率直方图

Fig4.19 Step size=1min,the posterior probability histogram of parameter T（min）

表4.7 步长=1min，参数T（min）落在各区域的频数统计

Table4.7 Step size =1min,frequency statistics of parameter T（min）falling in each region

区域（min）

T（min）

频数

50

1900

50

40～50

50～70

70～80

图4.20 步长=3min，参数T（min）的后验概率直方图

Fig4.20 Step size=3min,the posterior probability histogram of parameter T（min）

从图4.19可以看出，当步长=1min时，污染物瞬时排放时刻T（min）的后验概率密度分布非常集中，污染源排放时间T在60min附近的频率最大，基本符合预先设定值，从图可以推断污染事件的排放时间大概是T=[53min,71min]之间；从

46

4 排水管道污染物溯源模型的验证与应用

图4.20可以看出，当步长=3min时，由于步长过大且为整数的原因使抽样没有在60min左右游走，使抽样结果精度降低，这样抽样结果不符合预先设定值。

通过上述分析，污水管道污染物瞬时排放的未知反演参数排放节点Jx，在贝叶斯MCMC抽样算法中游走步长为1min时，可以提高反演精度和效率。这是因为在模拟中时间序列是以分钟为最小单位且为整数，在MCMC抽样反演中如果以大于1min的步长游走搜索，则无法在整个后验空间内搜索，可能跳过真实值，使抽样结果无法反演真实的后验分布。

3）结果讨论

通过上述分析，污水管道污染物瞬时排放的反演参数T（min），在贝叶斯MCMC抽样算法中游走步长为1min时，反演精度高。

如表4.8所示，当游走步长=1min时，通过将污染物排放时刻T（min）的抽样结果与真实值对比，可知污染物排放时刻的数学期望和中值误差均较小。污染源参数T（min）的反演结果误差柱状图见图4.21。

表4.8 污染源参数T（min）的抽样结果和实际值比较

Table4.8 Comparison of sampling results and actual values of pollution source parameters T

(min)

参数 T（min）

实际值 60

中值 58

中值误差/%

3.33

均值 58.3945

均值误差/%

2.68

T(min) 420中值误差/% 均值误差/%3.332.68 图4.21 污染源参数T（min）的反演结果误差

Fig4.21 Inversion error of pollution source parameter T（min）

④ 瞬时监测与连续监测及采样监测的时间间隔对反演效果的影响

在上述先验信息和排放量步长=10000、排放时间步长=1min、排放节点步长=1的基础上，对污染源单个未知参数进行MCMC抽样，比较分析瞬时监测与连续监测及采样监测的时间间隔对反演算法有效性及精度的影响。

1）瞬时监测与连续监测对反演效果的影响

上节中在一个监测点观测六个数据，通过对单个未知参数进行贝叶斯MCMC

47

重庆大学硕士学位论文

抽样基本可以反演其后验概率分布，本节采用瞬时监测数据对单变量参数进行贝叶斯MCMC抽样，比较分析瞬时监测与连续监测对反演算法有效性及精度影响。

a）反演结果分析

节点1为监测点，在1.5h时观测1个数据（见图4.22），按照MCMC反演思路迭代2000次。排放时间T（min）为反演参数的后验概率直方图见图4.23，排放量M（g）为反演参数的后验概率直方图见图4.24，排放节点Jx为反演参数的后验概率直方图见图4.25，参数T（min）落在各区域的频数统计见表4.9，参数M（g）落在各区域的频数统计见表4.10，参数Jx落在各区域的频数统计见表4.11。

图4.22 监测点J1处观测一个数据的浓度序列

Table4.22 The concentration sequence of one observation data at the monitoring point J1

图4.23 参数T（min）的后验概率直方图

Fig4.23 The posterior probability histogram of parameter T（min）

表4.9 参数T（min）落在各区域的频数统计

Table4.9 Frequency statistics of parameterT（min）falling in each region

区域（min）

T（min）

频数

800

750

450

50～70

70～90

90～110

48

4 排水管道污染物溯源模型的验证与应用

图4.24 参数M（g）的后验概率直方图

Fig4.24 The posterior probability histogram of parameter M（g）

表4.10 参数M（g）落在各区域的频数统计

Table4.10 Frequency statistics of parameter M（g）falling in each region

区域（×106g）

M（g）

频数

1200

600

200

0.7～0.9

0.9～1.1

1.1～1.3

图4.25 参数Jx的后验概率直方图

Fig4.25 The posterior probability histogram of parameter Jx

表4.11 参数Jx落在各区域的频数统计

Table4.11 Frequency statistics of parameter Jx falling in each region

区域（x）

Jx

频数

1050

950

0～20

20～40

从图4.23、4.24、4.25可以看出，当只有一个观测数据时，单个未知反演参数排放时刻T（min）、排放量M（g）和排放节点Jx的后验概率密度均较离散，比较难推断污染源的相关信息，排放时刻T（min）为反演参数时只能推断出污染事件大概发生在T =[55min，100min]，排放量M（g）反演参数时只能推断出污染物排放量大概发生在M=[790000g，1200000g]之间，排放节点Jx为反演参数时只能推断出污染物排放节点大概发生在Jx=[2，33]之间。

如表4.12所示，当只有一个观测数据时，通过将污染物排放时刻T（min）、

49

重庆大学硕士学位论文

排放量M（g）的抽样结果与真实值对比，可知污染物排放量的数学期望和中值误差均很多、大。污染源参数T（min）、M（g）的反演结果误差柱状图见图4.26，污染源参数Jx的后验统计结果见图4.27。

表4.12 一个观测数据时污染源参数T（min）、M（g）的抽样结果和实际值对比 Table4.12 Comparison of sampling results and actual values of pollution source parameters T

(min) and M (g) for one observation data

参数 T（min） M（g/min）

实际值 60 100000

403020100T（min）中值误差/% M（g）均值误差/%中值 80 887880

中值误差/% 33.33 11.21

均值 77.3525 934730

均值误差/% 28.92 6.53

33.3328.9211.216.53 图4.26 一个观测数据时污染源参数T（min）、M（g）的反演结果误差 Fig4.26 Error of inversion result of pollution source parameters T (min) and M (g) in one

observation data

图4.27 一个观测数据时污染源参数Jx的后验统计结果

Fig4.27 The posteriori statistical result of pollution source parameter Jx for an observation data

50

4 排水管道污染物溯源模型的验证与应用

b）结果讨论

1个观测数据和时间间隔为10min的6个观测数据，单参数M（g）、T（min）的反演误差比较见表4.13和图4.28，单参数Jx的反演结果对照见图4.29。很明显可以看出，6个观测数据的MCMC抽样反演误差远小于1个观测数据的反演误差，可见，相对于1个观测数据而言，6个数据观据所反演的后验概率密度更有效，精度更高。导致这种结果的原因是1个观测数据即一个点无法唯一确定一条浓度序列，而真实值在监测点的浓度序列是唯一的，或者接近真实值的浓度序列也是相对真实浓度序列发生微小变化，这样的后验值是可以接受的，但是经过一个点的浓度序列是千差万别，相应的反演的参数也波动很大，误差很大，无法反演真实的后验概率分布。时间间隔为10min的6个观测数据则会减少反演结果的不确定性，误差相对较小。

表4.13 1个观测数据和6个观测数据，参数M（g）、T（min）反演误差对照

Table4.13 comparison of the inversion results of parameters M (g) and T (min) with one observation

data and six observation data

污染物排放量M（g）

项目 中值误差/% 均值误差/%

1个观测数据

11.21 6.53

6个观测数据

1.25 2.48

污染物排放时间T（min） 1个观测数据

33.33 28.92

33.3328.926个观测数据

3.33 2.68

3530252015105011.216.531.252.48瞬时监测连续监测瞬时监测3.332.68连续监测污染物排放量M（g）中值误差/% 污染物排放时间T（min）均值误差/% 图4.28 1个观测数据和6个观测数据，参数M（g）、T（min）的反演误差对照

Fig4.28 comparison of the inversion results of parameters M (g) and T (min) with one observation

data and six observation data

51

重庆大学硕士学位论文

图4.29 1个观测数据和6个观测数据，参数Jx的反演结果对照

Fig4.29 comparison of the inversion results of parameters Jx with one observation data and six

observation data

2）采样监测的时间间隔对反演效果的影响

本节在监测节点1处6个观测数据的基础上，研究观测数据时间间隔对反演算法精度和效率的影响。

a）采样监测的时间间隔为5min反演结果分析

节点1为监测点，检测数据为6个，时间间隔为5min，观测浓度序列如图4.30，按照MCMC反演思路迭代2000次。排放节点Jx为反演参数的后验概率直方图见图4.31,排放时间T（min）为反演参数的后验概率直方图见图4.32，排放量M（g）为反演参数的后验概率直方图见图4.33，参数Jx落在各区域的频数统计见表4.14，参数M（g）落在各区域的频数统计见表4.15，参数T（min）落在各区域的频数统计见表4.16。

图4.30 检测时间间隔为5min的观测浓度序列

Fig4.30 Detection concentration sequence with a time interval of 5 min

52

4 排水管道污染物溯源模型的验证与应用

图4.31 时间间隔5min，参数Jx的后验概率直方图

Fig4.31 Time Interval= 5 min, the posterior probability histogram of the parameter Jx

表4.14 时间间隔5min，参数Jx落在各区域的频数统计

Table4.14 Interval 5 min, Frequency statistics of parameter Jx falling in each region

区域（x）

Jx

频数

1300

700

0～20

20～40

图4.32 时间间隔5min，参数M（g）的后验概率直方图

Fig4.32 Time Interval = 5 min, the posterior probability histogram of the parameter M（g）

表4.15 时间间隔5min，参数M（g）落在各区域的频数统计

Table4.15 Interval = 5 min, Frequency statistics of parameter M（g）falling in each region

区域（×106g）

M（g）

频数

50

1800

150

0.8～0.9

0.9～1.1

1.1～1.3

图4.33 时间间隔5min，参数T（min）的后验概率直方图

Fig4.33 Interval = 5 min, the posterior probability histogram of the parameter T（min）

53

重庆大学硕士学位论文

表4.16 时间间隔5min，参数T（min）落在各区域的频数统计

Table4.16 Interval = 5 min,Frequency statistics of parameter T（min）falling in each region

区域（min）

T（min）

频数

1900

100

50～70

70～90

从图4.31、4.32、4.33可以看出，当观测数据时间间隔为5min时，单个未知反演参数排放节点Jx和排放量M（g）的后验概率密度均离散，比较难推断污染源的相关信息，排放节点Jx为未知反演参数时只能推断大概发生在Jx=[1，35]，排放量M（g）为未知反演参数时只能推断大概发生在M=[880000g，1100000g]之间，排放时刻T（min）的后验概率分布比较集中在真实值60min附近，可以推断出污染事件大概发生在T =[55min，76min]之间。

如表4.17所示，当观测数据时间间隔为5min时，通过将污染物排放时刻T（min）、排放量M（g）、排放节点Jx的抽样结果与真实值对比，可知污染物排放时间和排放量的数学期望和中值误差均较小，但是排放时间和排放量的后验概率分布并不集中在实际值附近，排放节点Jx的后验统计结果分布比较广泛。时间间隔5min的反演结果误差见图4.34，污染源参数Jx的后验统计结果见图4.35。

表4.17 时间间隔5min的抽样结果和实际值对比

Table4.17 Comparison of sampling results and actual values of time interval 5min 参数 T（min） M（g/min）

实际值 60 100000

420T（min）中值误差/% M（g）均值误差/%中值 62 972870

中值误差/%

3.33 2.71

均值 61.9654 991090

均值误差/%

3.28 0.89

3.333.282.710.89 图4.34 时间间隔5min的反演结果误差 Fig4.34 Inversion result error of time interval 5min

54

4 排水管道污染物溯源模型的验证与应用

图4.35 污染源参数Jx的后验统计结果

Fig4.35 The posterior statistical results of pollution source parameters Jx

b）采样监测的时间间隔为15min反演结果分析

节点J1为监测点，检测数据为6个，时间间隔为15min，观测浓度序列如图4.36，按照MCMC反演思路迭代2000次。排放时间T（min）为反演参数的后验概率直方图见图4.37，排放量M（g）为反演参数的后验概率直方图见图4.38，排放节点Jx为反演参数的后验概率直方图见图4.39，时间间隔15min参数T（min）落在各区域的频数统计见表4.18，时间间隔15min参数M（g）落在各区域的频数统计见表4.19，时间间隔15min参数Jx落在各区域的频数统计见表4.20。

图4.36 检测时间间隔为15min的观测浓度序列

Fig4.36 Detection concentration sequence with a time interval of 15min

55

重庆大学硕士学位论文

图4.37 时间间隔15min，参数T（min）的后验概率直方图

Fig4.37 Interval = 5 min, the posterior probability histogram of the parameter T（min）

表4.18 时间间隔15min，参数T（min）落在各区域的频数统计

Table4.18 Interval = 5 min, Frequency statistics of parameter T（min）falling in each region

区域（min）

T（min）

频数

2000 50～70

图4.38 时间间隔15min，参数M（g）的后验概率直方图

Fig4.38 Interval = 15min, the posterior probability histogram of the parameter M（g）

表4.19 时间间隔15min，参数M（g）落在各区域的频数统计

Table 4.19 Interval = 15min, Frequency statistics of parameter M（g）falling in each region

区域（×106g）

M（g）

频数

300

1400

300

0.7～0.9

0.9～1.1

1.1～1.3

图4.39 时间间隔15min，参数Jx的后验概率直方图

Fig4.39 Interval = 15min, the posterior probability histogram of the parameterJx

56

4 排水管道污染物溯源模型的验证与应用

表4.20 时间间隔15min，参数Jx落在各区域的频数统计

Table 4.20 Interval = 15min, frequency statistics of parameter Jx falling in each region

区域（x）

Jx

频数

300

1700

10～20

20～40

从图4.37、4.38、4.39可以看出，当观测数据时间间隔为15min时，单个未知反演参数排放量M（g）和排放时刻T（min）的后验概率密度均集中分布在实际值附近，且污染物的排放时刻T（min）呈现比较明显的正态分布。从图中可以推断，当排放量M（g）为未知反演参数时发生在M=[850000g，1090000g]之间，当排放节点Jx为未知反演参数时发生在Jx=[16，32]之间，当排放时间T（min）为未知反演参数时发生在T=[54min，66min]之间。

如表4.21所示，当观测数据时间间隔为15min时，通过将污染物排放时刻T（min）、排放量M（g）、排放节点Jx的抽样结果与真实值对比，可知污染物泄时间和排放量的数学期望和中值误差均很小，排放节点Jx的后验统计结果集中分布在真实排放点的周围。时T（min）、M（g）的抽样结果和实际值对比见图4.40，时间间隔15min参数Jx的后验统计结果见图4.41。

表4.21 时间间隔15min的抽样结果和实际值对比

Table4.21 Comparison of sampling results and actual values of time interval 15min 参数 T（min） M（g/min）

实际值 60 100000

21.510.50T（min）中值误差/% M（g）均值误差/%中值 61 985920

中值误差/%

1.67 1.41

均值 60.3240 999020

均值误差/%

0.54 0.10

1.670.540.11.41 图4.40 时间间隔15min的反演结果误差 Fig4.40 Inversion result error of time interval 15min

57

重庆大学硕士学位论文

图4.41 时间间隔15min，参数Jx的后验统计结果

Fig4.41 Interval = 15min,the posterior statistical results of parameters Jx

c）结果讨论

污染物监测节点J1处观测6个数据，时间间隔分别为5min、10min、15min的单参数M（g）、T（min）的反演误差比较见表4.22和图4.42，很明显可以看出，时间间隔15min的反演误差小于时间间隔10min的反演误差，时间间隔10min的反演误差小于时间间隔5min的反演误差。可见，相对时间间隔5min和10min而言，时间间隔15min所反演的后验概率密度更有效，精度更高。时间间隔5min、10min、15min 参数Jx的后验统计结果见图4.43，明显可以看出，时间间隔为15min时所反演的可能排放节点更集中分布在真实值附近。

通过以上分析得出，6个数据的时间间隔为15min时能更好的反映节点1处观测到的污染物浓度序列，所以，通过贝叶斯MCMC抽样算法识别污水管网突发污染事件污染源源项时，要尽可能使监测点观测到的数据能够反映此监测点污染物浓度序列。

表4.22 时间间隔5min、10min、15min，M（g）、T（min）的反演误差对照 Table4.22 Interval = 5min、10min and 15min,comparison of the inversion results of parameters

M (g) and T (min)

污染物排放量M（g）

项目

5min

中值误差/% 均值误差/%

33.33 28.29

10min 1.25 2.48

15min 1.41 0.1

5min 3.33 3.28

10min 3.33 2.68

15min 1.67 0.54

污染物排放时间T（min）

58

4 排水管道污染物溯源模型的验证与应用

353025201510505min10min2.481.251.410.115min3.333.283.332.6810min1.670.5415min33.3328.295min污染物排放量M（g）中值误差/% 污染物排放时间T（min）均值误差/% 图4.42 时间间隔5min、10min、15min，M（g）、T（min）的反演误差对照 Fig4.42 Interval = 5min、10min and 15min,comparison of the inversion results of parameters M

(g) and T (min)

图4.43 时间间隔5min、10min、15min,参数Jx的后验统计结果

Fig4.43 Interval = 5min、10min and 15min,the posterior statistical results of parameters Jx

4.2.3 两个未知参数反演

在三个污染物质排放特征参数（排放节点、排放量、排放时刻）中，假定有两个参数未知需要通过溯源模型获得信息，而其余一个排放参数已知。在监测点实测污染物浓度数据的基础上，利用溯源模型进行未知参数的反演推理，获得该

59

重庆大学硕士学位论文

参数的概率分布。因此包括三个情景的反演计算，分别为：污染物排放量和排放时刻同时反演、排放节点和排放量同时反演、排放时刻和排放节点同时反演。

在4.2.1节分析结果的基础上，对污染源两个未知参数进行贝叶斯MCMC抽样。

① 未知参数排放量和排放时间反演分析 1）排放时间和排放量的反演分析

假定某日在T（min）时，在节点28处的某工厂瞬时排入重量为M（g）的污染物BX。

根据先验信息，待反演模型参数先验分布为均匀分布，其对应的先验概率密度函数分别为：

0.1×10−5， 500000≪M≪1500000

P(M)={ （4.10）

0， 其它0.0083， 0≪T≪120

()PT={ （4.11）

0， 其它

由于一般情况观测噪声为白噪声，所以假定观测误差服从标准差为σ=0.01的正态分布N(0，σ2)，则污染物溯源的似然函数可表示为：

p(M，T)=(2πσ2)p((M，T)=

式中：α为比例常数。

上式即为排水管网污染事件溯源的未知参数T（min）、M（g）的后验概率密度函数，通过MCMC抽样算法对模型进行抽样，就可以得到未知参数T（min）、M（g）的后验概率密度函数分布。

2）排放时间和排放量的反演结果分析

按照MCMC反演思路迭代2000次。反演参数M（g）、T（min）的后验概率直方图分别见图4.44、图4.45，参数M（g）落在各区域的频数统计见表4.23，参数T（min）落在各区域的频数统计见表4.24。

1

∑66/2exp(i=1[−

(Ci(M，T)−Ci)2

2σ2

]) （4.12）

根据贝叶斯定理，模型参数的后验概率密度函数为：

(Ci(M，T)−Ci)

αexp(∑6[−i=12σ2

2

]) （4.13）

图4.44 参数M（g）的后验概率直方图

Fig4.44 The posterior probability histogram of parameter M（g）

60

4 排水管道污染物溯源模型的验证与应用

表4.23 参数M（g）落在各区域的频数统计

Table 4.23 Frequency statistics of parameter M（g）falling in each region

区域（×106g）

M（g）

频数

1600

400

0.9～1.1

1.1～1.2

图4.45 参数T（min）的后验概率直方图

Fig4.45 The posterior probability histogram of parameter T（min）

表4.24 参数T（min）落在各区域的频数统计

Table 4.24 Frequency statistics of parameter T（min）falling in each region

区域（min）

T（min）

频数

700

1280

20

40～50

50～70

70～80

从图4.44、4.45可以看出，反演参数排放量M（g）和排放时刻T（min）的后验概率密度分布均较离散，排放量在大于真实值的1100000g附近有较多分布，排放时间在小于真实值的50min附近有较多分布，从图中可以推断，污染事件的反演参数排放量发生在M=[920000g，1130000g]之间，排放时刻发生在T=[42min,73min]之间。根据Ct=C0e−kt可知，当排放节点一定时，排放量增大、排放时间提前，在同一监测点同一时间观测到的污染物浓度是有可能一样或者相近，所以这样的反演结果是存在的也是合理的。

如表4.25所示，通过将污染事件的反演参数排放量M（g）和排放时刻T（min）的抽样结果与真实值对比，可知污染物泄时间和排放节点的数学期望和中值误差均比较小，在10%之内。参数M（g）、T（min）反演结果误差见图4.46。

表4.25 参数M（g）、T（min）抽样结果和实际值对比

Tab4.25 Comparison of sampling results and actual values of parameters M (g) and T (min) 参数 M（g） T（min）

实际值 1000000 60

中值 1056340 59

中值误差/%

5.63 1.67

均值 1043670 57.0870

均值误差/%

4.37 4.86

61

重庆大学硕士学位论文

64204.861.675.634.37T（min）中值误差/% M（g）均值误差/% 图4.46 参数M（g）、T（min）反演结果误差 Fig4.46 Inversion result error of parameters M (g) and T (min)

② 未知参数排放时间和排放节点反演分析 1）排放时间和排放节点的反演分析

假定某日在T（min）时，在节点Jx处的某工厂瞬时排入重量为1000000g的污染物BX。

根据先验信息，待反演模型参数先验分布为均匀分布，其对应的先验概率密度函数分别为：

P(T)={

0.0083， 0≪T≪1200， 其它

0.0156， 1≪x≪65

0， 其它

（4.14）

P(Jx)={ （4.15）

由于一般情况观测噪声为白噪声，所以假定观测误差服从标准差为σ=0.01的正态分布N(0，σ2)，则污染物溯源的似然函数可表示为：

p(T，Jx)=

(Ci(T，Jx)−Ci)6

∑exp([−i=1(2πσ2)6/22σ2

1

2

]) （4.16）

根据贝叶斯定理，模型参数的后验概率密度函数为：

p((T，Jx)=

式中：β为比例常数。

上式即为排水管网污染事件溯源的未知参数T，Jx的后验概率密度函数，通过MCMC抽样算法对模型进行抽样，就可以得到未知参数T，Jx的后验概率密度函数分布。

2）排放时间和排放节点的反演结果分析

按照MCMC反演思路迭代2000次。T、Jx参数的后验概率直方图分别见图4.47、图4.48，参数T（min）落在各区域的频数统计见表4.26，参数Jx落在各区域的频数统计见表4.27。

(Ci(T，Jx)−Ci)

βexp(∑6[−i=12σ2

2

]) （4.17）

62

4 排水管道污染物溯源模型的验证与应用

图4.47 反演参数T（min）的后验概率直方图

Fig4.47 The posterior probability histogram of parameter T（min）

表4.26 参数T（min）落在各区域的频数统计

Table 4.26 Frequency statistics of parameter T（min）falling in each region

区域（min）

T（min）

频数

300

1700

40～50

50～70

图4.48 反演参数Jx的后验概率直方图

Fig4.48 The posterior probability histogram of parameter Jx

表4.27 参数Jx落在各区域的频数统计

Table 4.27 Frequency statistics of parameter Jx falling in each region

区域（x）

Jx

频数

1700

300

20～40

40～50

从图5.47、4.48可以看出，反演参数排放时刻T（min）的后验概率密度集中分布在小于实际值的57min附近，排放节点Jx的后验概率密度分布较离散，从图可推断污染事件排放时间大概发生在T=[47min，65min]之间，排放节点大概发生在Jx=[26，43]之间。根据Ct=C0e

−k

xu

可知，当排放量一定时，由于管段的长度和

水流速度等因素，污染物从不同节点提前排放，在同一监测点同一时间观测到的污染物浓度是有可能一样或者相近，所以这样的反演结果是存在的也是合理的。

如表4.28所示，通过将污染事件的反演参数排放时间T（min）的抽样结果与

63

重庆大学硕士学位论文

真实值对比，可知污染物泄时间的数学期望和中值误差均小，在10%内，参数T（min）的抽样结果误差见图4.49。如图4.50所示，污染源参数Jx的后验统计结果集中分布在真实排放点的周围，在很大程度上缩小了监测范围，可能排放节点为整个管网的四分之一。

表4.28 参数T（min）的抽样结果和实际值对比

Table4.28 Comparison of sampling results and actual values of parameters T (min) 参数 T（min）

实际值 60

中值 57

中值误差/%

5.00

均值 56.1175

均值误差/%

6.47

T（min）1050中值误差/% 均值误差/%56.47 图4.49 参数T（min）的反演结果误差 Fig4.49 Inversion result error of parameters T (min)

图4.50 污染源参数Jx的后验统计结果

Fig4.50 The posterior statistical results of pollution source parameters Jx

64

4 排水管道污染物溯源模型的验证与应用

③ 未知参数排放量和排放节点反演分析 1）排放量和排放节点的反演分析

假定某日在1:00时，在节点Jx处的某工厂瞬时排入重量为M（g）的污染物BX。

根据先验信息，待反演模型参数先验分布为均匀分布，其对应的先验概率密度函数分别为：

0.1×10−5， 500000≪M≪1500000

P(M)={ （4.18）

0， 其它0.0156， 1≪x≪65

P(Jx)={ （4.19）

0， 其它

由于一般情况观测噪声为白噪声，所以假定观测误差服从标准差为σ=0.01的正态分布N(0，σ2)，则污染物溯源的似然函数可表示为：

p(M，Jx)=(2πσ2)

p((M，Jx)=

式中：γ为比例常数。

上式即为排水管网污染事件溯源的未知参数M，Jx的后验概率密度函数，通过MCMC抽样算法对模型进行抽样，就可以得到未知参数M，Jx的后验概率密度函数分布。

2）排放量和排放节点的反演结果分析

按照MCMC反演思路迭代2000次。M、Jx参数的后验概率直方图分别见图4.51、图4.52，参数M（g）落在各区域的频数统计见表4.29，参数Jx落在各区域的频数统计见表4.30。

1

∑66/2exp(i=1[−

(Ci(M，Jx)−Ci)2

2σ2

]) （4.20）

根据贝叶斯定理，模型参数的后验概率密度函数为：

(Ci(M，Jx)−Ci)

γexp(∑6[−i=12σ2

2

]) （4.21）

图4.51 反演参数M（g）的后验概率直方图

Fig4.51 The posterior probability histogram of parameter M（g）

65

重庆大学硕士学位论文

表4.29 参数M（g）落在各区域的频数统计

Table 4.29 Frequency statistics of parameter M（g）falling in each region

区域（×106g）

M（g）

频数

800

1200

0.9～1.1

1.1～1.3

图4.52 反演参数Jx的后验概率直方图

Fig4.52 The posterior probability histogram of parameter Jx

表4.30 参数Jx落在各区域的频数统计

Table 4.30 Frequency statistics of parameter Jx falling in each region

区域（x）

Jx

频数

100

1400

500

10～20

20～40

40～60

从图4.51、4.52可以看出，反演参数排放量M（g）和排放节点Jx的后验概率密度分布较离散，且均偏离真实值，只能大概推测排放量在M=[890000g，1300000g]之间，排放节点在Jx=[15，52]之间。根据Ct=C0e

−k

xu

可知，当排放时间

一定时，由于管段的长度和水流速度等因素，不同量的污染物从不同节点排放，在同一监测点同一时间观测到的污染物浓度是有可能一样或者相近，所以这样的反演结果是存在的也是合理的。

如表4.31所示，通过将污染事件的反演参数排放量M（g）的抽样结果与真实值对比，可知污染物泄量的均值误差和中值误差均较大，参数Jx的抽样结果误差见图4.53。如图4.54所示，污染源Jx参数的后验统计结果集中分布在真实排放点的周围，但监测范围比较广。

表4.31 参数M(g)抽样结果和实际值对比

Table4.31 Comparison of sampling results and actual values of parameters M(g) 参数 M(g)

实际值 1000000

中值 1000000

66

中值误差/% 14.37

均值 1150300

均值误差/% 15.03

4 排水管道污染物溯源模型的验证与应用

M(g) 1614.3714中值误差/% 15.03均值误差/% 图4.53 参数M(g)的反演结果误差 Fig4.53 Inversion result error of parameters M(g)

图4.54 污染源参数Jx的后验统计结果

Fig4.54 The posterior statistical results of pollution source parameters Jx

4.2.4 三个未知参数反演

假定三个排放特征参数（排放节点、排放量、排放时刻）均未知，需要经过溯源模型提供信息。污染物从节点瞬时排入管网，排放量为M，排放时刻为模拟开始后T，排放节点参数Jx。

在4.2.1节分析结果的基础上，对污染源三个未知参数进行贝叶斯MCMC抽样。

① 三个未知参数反演分析

1） 三个未知参数反演分析

假定某日在T（min）时，在节点Jx处的某工厂瞬时排入重量为M（g）的污染物BX。

根据先验信息，待反演模型参数先验分布为均匀分布，其对应的先验概率密度函数分别为：

67

重庆大学硕士学位论文

1×10−5， 500000≪M≪1500000

P(M)={ （4.22）

0， 其它0.0083， 0≪T≪120

P(T)={ （4.23）

0， 其它0.0156， 1≪x≪65

P(Jx)={ （4.24）

0， 其它

由于一般情况观测噪声为白噪声，所以假定观测误差服从标准差为σ=0.01的正态分布N(0，σ2)，则污染物溯源的似然函数可表示为：

p(M，T，Jx)=

(Ci(M，T，Jx)−Ci)6

∑exp([−i=1(2πσ2)6/22σ2

1

2

]) （4.25）

根据贝叶斯定理，模型参数的后验概率密度函数为：

p((M，T，Jx)=αexp(∑6i=1[−

式中：α为比例常数。

上式即为排水管网污染事件溯源的未知参数M，T，Jx的后验概率密度函数，通过MCMC抽样算法对模型进行抽样，就可以得到未知参数M，T，Jx的后验概率密度函数分布。

2） 三个未知参数反演结果分析

按照MCMC反演思路迭代4000次。参数M（g）、T（min）、Jx的后验概率直方图见图4.55、4.56、4.57，参数M（g）落在各区域的频数统计见表4.32，参数T（min）落在各区域的频数统计见表4.33，参数Jx落在各区域的频数统计见表4.34。

(Ci(M，T，Jx)−Ci)2

2σ2

]) （4.26）

图4.55 反演参数M（g）的后验概率直方图

Fig4.55 The posterior probability histogram of parameter M（g）

表4.32 参数M（g）落在各区域的频数统计

Table4.32 Frequency statistics of parameter M（g）that falling in each region

区域（×106g）

M（g）

频数

2200

1800

0.9～1.1

1.1～1.2

68

4 排水管道污染物溯源模型的验证与应用

图4.56 反演参数T（min）的后验概率直方图

Fig4.56 The posterior probability histogram of parameter T（min）

表4.33 参数T（min）落在各区域的频数统计

Table4.33 Frequency statistics of parameter T（min）that falling in each region

区域（min）

T（min）

频数

1800

2000

200

30～50

50～70

70～80

图4.57 反演参数Jx的后验概率直方图

Fig4.57 The posterior probability histogram of parameter Jx

表4.34 参数Jx落在各区域的频数统计

Table4.34 Frequency statistics of parameter Jx that falling in each region

区域（x）

Jx

频数

100

3400

500

10～20

20～40

40～50

从图4.55、4.56、4.57可以看出，污染物排放量、排放时间、排放节点的后验概率密度较离散，只能大概推断污染事件的反演参数排放量在M=[900000g，1260000g]，排放时间在T=[35min，72min]之间，排放节点在Jx=[20，45]之间。根据Ct=C0e

−k

xu

可知，对于三个参数的反演，反演结果的可组合性更多，贝叶斯统

计MCMC抽样的目的就是找出可能的值，缩小先验范围，着重检查概率分布大的点，但也不能忽略概率小的点。

如表4.35所示，通过将污染事件的反演参数排放量M（g）、排放时间T（min）

69

重庆大学硕士学位论文

的抽样结果与真实值对比，可知污染物排放时间和排放量中值误差和数学期望都比较小，参数M（g）、T（min）反演结果误差见图4.58。如图4.59所示，污染源参数Jx的后验统计结果集中分布在真实排放点的周围，可能排放量节点范围为整个管网的三分之一。

表4.35 反演参数抽样结果和实际值对比

Table4.35 Comparison of sampling results and actual values of inversion parameters 参数 M（g） T（min）

实际值 1000000 60

151050T（min）中值误差/% M（g）均值误差/%中值 1085200 52

中值误差/%

8.52 13.33

均值 1093100 52.6320

均值误差/%

9.31 12.28

13.3312.288.529.31图4.58 反演参数M（g）、T（min）反演结果误差 Fig4.58 Inversion result error of parameters M (g) and T (min)

图4.59 污染源参数Jx的后验统计结果

Fig4.59 The posterior statistical results of pollution source parameters Jx

70

4 排水管道污染物溯源模型的验证与应用

② 监测点的空间布置对溯源效果的影响

增加水质监测节，比较分析监测点的空间布置对溯源效果的影响。图4.60为某小型污水管网平面布置图，节点21为增设的水质监测点。

按照MCMC反演思路迭代4000次。参数M（g）、T（min）、Jx的后验概率直方图见图4.61、4.62、4.63，参数M（g）落在各区域的频数统计见表4.36，参数T（min）落在各区域的频数统计见表4.37，参数Jx落在各区域的频数统计见表4.38。

图4.60 污水管网平面布置图 Fig4.60 Layout drawing of sewage networks

图4.61 反演参数M（g）的后验概率直方图

Fig4.61 The posterior probability histogram of parameter M（g）

71

重庆大学硕士学位论文

表4.36 参数M（g）落在各区域的频数统计

Table4.36 Frequency statistics of parameter M（g）that falling in each region

区域（×106g）

M（g）

频数

4000 0.9～1.1

图4.62 反演参数T（min）的后验概率直方图

Fig4.62 The posterior probability histogram of parameter T（min）

表4.37 参数T（min）落在各区域的频数统计

Table4.37 Frequency statistics of parameter T（min）that falling in each region

区域（min）

T（min）

频数

2100

1900

50～70

70～80

图4.63 反演参数Jx的后验概率直方图

Fig4.63 The posterior probability histogram of parameter Jx

表4.38 参数Jx落在各区域的频数统计

Table4.38 Frequency statistics of parameter Jx that falling in each region

区域（x）

Jx

频数

4000 20～40

从图4.61、4.62、4.63可以看出，污染物排放节点的后验概率密度分布集中，

72

4 排水管道污染物溯源模型的验证与应用

可推断污染事故反演参数排放节点在Jx=[28，35]之间；污染物排放量、排放时间的后验概率密度较离散，只能大概推断反演参数排放量在M=[910000g，1090000g]之间，排放时间在T=[57min，80min]之间。根据Ct=C0e

−k

x

u

可知，对于三个参数

的反演，反演结果的可组合性更多，贝叶斯统计MCMC抽样的目的就是找出可能的值，缩小先验范围，着重检查概率分布大的点，但也不能忽略概率小的点。

如表4.39所示，通过将污染事件的反演参数排放量M（g）、排放时间T（min）的抽样结果与真实值对比，可知污染物排放量中值误差和数学期望都比较小，排放时间的中值误差和数学期望都比较大，参数M（g）、T（min）反演结果误差见图4.64。如图4.65所示，污染源参数Jx的后验统计结果集中分布在真实排放点的周围，可能排放量节点范围为整个管网的三分之一。

从反演结果可以得出，增加水质监测点可以提高溯源精度，尤其是对排放节点的溯源更加精确。

表4.39 反演参数抽样结果和实际值对比

Table4.39 Comparison of sampling results and actual values of inversion parameters 参数 M（g） T（min）

实际值 1000000 60

20151050T（min）中值误差/% M（g）均值误差/%2.852.76中值 971520 69

中值误差/%

2.85 15.00

均值 972420 69.3560

均值误差/%

2.76 15.59

1515.59 图4.64 反演参数M（g）、T（min）反演结果误差 Fig4.64 Inversion result error of parameters M (g) and T (min)

73

重庆大学硕士学位论文

图4.65 污染源参数Jx的后验统计结果

Fig4.65 The posterior statistical results of pollution source parameters Jx

4.3 节点连续排放污染事件溯源

图4.10为某小型污水管网平面布置图，管网中水流流向如箭头所示，该污水管网共有管道64段，65个节点，管径为400～800mm，管道形状选CIRCULAR，其中PFK1为该污水管网的总排放口，在管网下游主干管上某节点监测处布设水质监测点，假定该污水管网水文条件稳定，满足一维水质模型应用条件。考察贝叶斯反演算法对于信息缺失度的敏感程度，在每种排放模式下，考察3种反演的工况：单个未知参数，两个未知参数和三个未知参数，分别代表排放节点、排放量、排放时间三个排放特征参数中有单个、两个、三个未知参数需要反推，其他参数已知的情形，比较分析反演参数的个数对贝叶斯算法有效性的影响。

假定某日1:00时，在节点28处的某工厂以100000g/min的速度连续排放污染物BX，在下游监测点1处1:05时观测到BX污染带前锋，现利用检测节点1处噪声的浓度场部分分布数据反演上游BX排放源的排放强度、排放位置和排放时间，为了利用SWMM模拟污染物的转输过程，设置BX在管道中的衰减系数为0.25。 监测点J1处的污染物浓度在模拟时间内不断上升，浓度最高达0.31mg/L，相关参数及先验信息见表4.40。

74

4 排水管道污染物溯源模型的验证与应用

表4.40 污水管SWMM模型水文参数

Table4.40 Parameters of SWMM model for sewage pipe

参数 类型 数值 参数 类型 数值

管径（mm）

已知 400～800

管段形状 已知 CIRCULAR

管段粗糙系数

已知 0.01

污染物衰减系数

已知 0.25 排放节点Jx 待反演 J28

模拟时间（h） 排放量 M（g/min） 排放时间T（min）

已知 3

待反演 100000

待反演 60

根据4.2.1节分析结果，污染物连续排放时，在监测点1处首次观测时间为1:15，后每隔15min观测一次，共8个数据，排放量M（g/min）游走步长为1000g/min，排放时间步长为1min，排放节点步长为1时，反演精度高，所有以此条件进行下面的贝叶斯MCMC统计反演。从图4.66监测点J1处的污染物浓度在模拟时间内不断上升，浓度最高达0.31mg/L。

图4.66 检测时间间隔为15min的观测浓度序列

Fig4.66 Detection concentration sequence with a time interval of 15min

4.3.1 单个未知参数反演

在三个污染物质排放特征参数（排放节点、排放量、排放时刻）中，假定只有一个参数未知需要通过溯源模型获得信息，而其余两个排放参数已知。

① 未知参数排放节点反演结果分析

假定某日在1:00时，在节点Jx处的某工厂以100000g/min的速度连续排放污染物BX。

按照MCMC反演思路迭代2000次。反演参数Jx的后验概率直方图分别见图4.67，参数Jx落在各区域的频数统计见表4.41。

75

重庆大学硕士学位论文

图4.67 连续排放反演参数Jx的后验概率直方图 Fig4.67 The posterior probability histogram of parameter Jx

表4.41 参数Jx落在各区域的频数统计

Table4.41 Frequency statistics of parameter Jx that falling in each region

区域（x）

Jx

频数

20

1980

10～20

20～40

从图4.67可以看出，污染物连续排放的反演参数排放节点的后验概率密度较离散，只能推测出污染事件排放点大概发生在Jx=[20，35]之间。

如图4.68所示，污染源Jx参数的后验统计结果集中分布在真实排放点的周围，在很大程度上缩小了监测范围，可能排放节点为整个管网的八分之一。

图4.68 污染源参数Jx的后验统计结果

Fig4.68 The posterior statistical results of pollution source parameters Jx

② 未知参数排放量反演结果分析

76

4 排水管道污染物溯源模型的验证与应用

假定某日在1:00时，在节点28处的某工厂以M（g/min）的速度连续排放污染物BX。

按照MCMC反演思路迭代2000次。对后验概率分布进行MCMC抽样，反演参数M（g/min）的后验概率直方图分别见图4.69，参数M（g/min）落在各区域的频数统计见表4.42。

图4.69 连续排放反演参数M（g/min）的后验概率直方图 Fig4.69 The posterior probability histogram of parameter M（g/min）

表4.42 参数M（g/min）落在各区域的频数统计

Table4.42 Frequency statistics of parameter M（g/min）that falling in each region

区域（×106g）

M（g/min）

频数

1900

100

0.9～1.1

1.1～1.4

从图4.69可以看出，污染物连续排放污染事件的反演参数排放量M（g/min）的后验概率密度分布非常集中，绝大多数集中在实际值100000（g/min）附近，且污染物排放强度呈现比较明显的正态分布。从图可以推断污染事件的排放量最可能发生在M=[96000g/min,104000g/min]之间。

如表4.43所示，通过将连续排放污染物参数排放量M（g/min）的抽样结果与真实值对比，可知污染物排放量的数学期望和中值误差均很小，污染源参数M（g/min）反演结果误差见图4.70。

表4.43 污染源参数M（g/min）的抽样结果与实际值对比

Table4.43 Comparison of sampling results and actual values of inversion parameters M（g/min） 参数 M（g/min）

实际值 100000

中值 99683

中值误差/%

0.32

均值 100360

均值误差/%

0.36

77

重庆大学硕士学位论文

M（g/min）0.40.350.3中值误差/% 均值误差/%0.360.32 图4.70 污染源参数M（g/min）反演结果误差 Fig4.70 Inversion result error of parameters M（g/min）

③ 未知参数排放时刻反演结果分析

假定某日在T时，在节点28处的某工厂以100000g/min的速度连续排放污染物BX。

按照MCMC反演思路迭代2000次。反演参数T（min）的后验概率直方图分别见图4.71，参数T（min）落在各区域的频数统计见表4.44。

图4.71 连续排放反演参数T（min）的后验概率直方图 Fig4.71 The posterior probability histogram of parameter T（min）

表4.44 参数T（min）落在各区域的频数统计

Table4.44 Frequency statistics of parameter T（min）that falling in each region

区域（min）

T（min）

频数

100

1800

100

40～50

50～70

70～80

从图4.71可以看出，污染物连续排放的反演参数初始排放时间T（min）的后验概率密度集中分布在55附近，可推测出污染事件初始排放时间发生在T=[49min，70min]之间。

如表4.45所示，通过将污染事件的反演参数初始排放时间T（min）的抽样结果与真实值对比，可知污染物排放节点Jx的数学期望和中值误差均较小，污染源参数T（min）反演结果误差见图4.72。

78

4 排水管道污染物溯源模型的验证与应用

表4.45 污染源参数T（min）的抽样结果与实际值对比

Table4.45 Comparison of sampling results and actual values of inversion parameters T（min） 参数 T（min）

实际值 60

中值 56

中值误差/%

6.67

均值 57.1060

均值误差/%

4.82

T（min）1050中值误差/% 均值误差/%6.674.82 图4.72 污染源参数T（min）反演结果误差 Fig4.72 Inversion result error of parameters T（min）

4.3.2 两个未知参数反演

在三个污染物质排放特征参数（排放节点、排放量、排放时刻）中，假定有两个参数未知需要通过溯源模型获得信息，而其余一个排放参数已知。在监测点实测污染物浓度数据的基础上，利用溯源模型进行未知参数的反演推理，获得该参数的概率分布。因此包括三个情景的反演计算，分别为：污染物排放量和排放时刻同时反演、排放节点和排放量同时反演、排放时刻和排放节点同时反演。

① 未知参数排放量和排放时间反演结果分析

假定某日在T时，在节点28处的某工厂以M（g/min）的速度连续排放污染物BX。

按照MCMC反演思路迭代2000次。反演参数M（g/min）、T（min）的后验概率直方图分别见图4.73、图4.74，参数M（g/min）落在各区域的频数统计见表4.46，参数T（min）落在各区域的频数统计见表4.47。

图4.73 连续排放反演参数M（g/min）的后验概率直方图 Fig4.73 The posterior probability histogram of parameter M（g/min）

79

重庆大学硕士学位论文

表4.46 参数M（g/min）落在各区域的频数统计

Table4.46 Frequency statistics of parameter M（g/min）that falling in each region

区域（×106g）

M（g/min）

频数

20

1980

0.8～0.9

0.9～1.1

图4.74 连续排放反演参数T（min）的后验概率直方图 Fig4.74 The posterior probability histogram of parameter T（min）

表4.47 参数T（min）落在各区域的频数统计

Table4.47 Frequency statistics of parameter T（min）that falling in each region

区域（min）

T（min）

频数

1200

800

30～50

50～70

从图4.73可以看出，污染物连续排放污染事件的反演参数排放量M（g/min）的后验概率密度分布非常集中，绝大多数集中在实际上100000（g/min）附近，且污染物排放强度呈现比较明显的正态分布，可推断排放量最可能发生在M=[94000g/min,105000g/min]之间；从图4.74可以看出，污染物连续排放的反演参数初始排放时间T（min）的后验概率密度较离散，只能推测出污染事件初始排放时间大概发生在T=[40min，62min]之间。

如表4.48所示，通过将污染事件的反演参数排放量M（g/min）、排放时间T（min）的抽样结果与真实值对比，可知污染物排放量M（g/min）的数学期望和中值误差均很小，排放时间T（min）的数学期望和中值误差均较大，污染源参数M（g/min）、T（min）反演结果误差见图4.75。

80

4 排水管道污染物溯源模型的验证与应用

表4.48 污染源参数M（g/min）、T（min）的抽样结果与实际值对比

Table4.48 Comparison of sampling results and actual values of inversion parameters M（g/min）

and T（min）

参数 M（g/min） T（min）

实际值 100000 60

2520151050T（min）中值误差/% 中值 100540 47

21.67中值误差/%

0.54 21.67

均值 100130 49.0815

均值误差/%

0.13 18.20

18.20.540.13M（g）均值误差/% 图4.75 污染源参数M（g/min）、T（min）反演结果误差 Fig4.75 Inversion result error of parameters M（g/min）and T（min）

② 未知参数排放时间和排放节点反演结果分析

假定某日在T时，在节点Jx处的某工厂以100000（g/min）的速度连续排放污染物BX。

按照MCMC反演思路迭代2000次。反演参数T（min）、Jx的后验概率直方图分别见图4.76、图4.77，参数T（min）落在各区域的频数统计见表4.49，参数Jx落在各区域的频数统计见表4.50。

图4.76 连续排放反演参数T（min）的后验概率直方图 Fig4.76 The posterior probability histogram of parameter T（min）

表4.49 参数T（min）落在各区域的频数统计

Table4.49 Frequency statistics of parameter T（min）that falling in each region

区域（min）

T（min）

频数

900

1100

40～80

80～1200

81

重庆大学硕士学位论文

图4.77 连续排放反演参数Jx的后验概率直方图 Fig4.77 The posterior probability histogram of parameter Jx

表4.50 参数Jx落在各区域的频数统计

Table4.50 Frequency statistics of parameter Jx that falling in each region

区域（x）

Jx

频数

550

1000

450

0～20

20～40

40～60

从图4.76、图4.77可以看出，污染物连续排放的反演参数初始排放时间T（min）、排放节点Jx的后验概率密度均较离散，只能大概推测出污染事件初始排放时间大概发生在T=[45min，110min]之间，排放节点大概在Jx=[13，50]之间。

如表4.51 所示，通过将连续污染事件的反演参数排放时间T（min）的抽样结果与真实值对比，排放时间T（min）的数学期望和中值误差均很大，污染源参数T（min）的反演结果误差见图4.78。如图4.79所示，污染源Jx参数的后验统计结果集中分布在真实排放点的周围，但范围比较广，可能排放节点为整个管网的二分之一。

表4.51 污染源参数M（g/min）的抽样结果与实际值对比

Table4.51 Comparison of sampling results and actual values of inversion parameters M（g/min） 参数 T（min）

实际值 60

中值 82

中值误差/% 36.67

均值 78.0670

均值误差/% 30.11

T（min）40200中值误差/% 均值误差/%36.6730.11 图4.78 污染源参数T（min）的反演结果误差 Fig4.78 Inversion result error of parameters T（min）

82

4 排水管道污染物溯源模型的验证与应用

图4.79 污染源参数Jx的后验统计结果

Fig4.79 The posterior statistical results of pollution source parameters Jx

③ 未知参数排放量和排放节点反演结果分析

假定某日在1:00时，在节点Jx处的某工厂以M（g/min）的速度连续排放污染物BX。

按照MCMC反演思路迭代2000次。反演参数M（g/min）、Jx的后验概率直方图分别见图4.80、图4.81，参数M（g/min）落在各区域的频数统计见表4.52，参数Jx落在各区域的频数统计见表4.53。

图4.80 参数M（g/min）的后验概率直方图

Fig4.80 The posterior probability histogram of parameter M（g/min）

83

重庆大学硕士学位论文

表4.52 参数M（g/min）落在各区域的频数统计

Table4.52 Frequency statistics of parameter M（g/min） that falling in each region

区域（×106g）

M（g/min）

频数

1950

50

0.9～1.1

1.1～1.2

图4.81 参数Jx的后验概率直方图

Fig4.81 The posterior probability histogram of parameter Jx

表4.53 参数Jx落在各区域的频数统计

Table4.53 Frequency statistics of parameter Jx that falling in each region

区域（x）

Jx

频数

500

1500

20～40

40～65

从图4.80可以看出，污染物连续排放污染事件的反演参数排放量M（g/min）的后验概率密度分布非常集中，绝大多数集中在实际上100000（g/min）附近，且污染物排放强度呈现比较明显的正态分布，可推断排放量最可能发生在M=[94000g/min,105000g/min]之间；从图4.81可以看出，污染物连续排放的反演参数初始排放时间Jx的后验概率密度较离散，只能大概推测出污染事件排放节点大概在Jx=[27，65]之间。

如表4.54所示，通过将污染事件的反演参数排放量M（g/min）的抽样结果与真实值对比，可知污染物排放量M（g/min）的数学期望和中值误差均很小，污染源参数M（g/min）的反演结果误差见图4.82。如图4.83所示，污染源Jx参数的后验统计结果分布不集中，分布范围比较广。

表4.54 污染源参数M（g/min）的抽样结果与实际值对比

Table4.54 Comparison of sampling results and actual values of inversion parameters M（g/min） 参数 M（g/min）

实际值 100000

中值 100230

84

中值误差/%

0.23

均值 100330

均值误差/%

0.33

4 排水管道污染物溯源模型的验证与应用

M（g）0.40.330.230.20中值误差/% 均值误差/% 图4.82 污染源参数M（g/min）的反演结果误差 Fig4.82 Inversion result error of parameters M（g/min）

图4.83 污染源参数Jx的后验统计结果

Fig4.83 The posterior statistical results of pollution source parameters Jx

4.3.3 三个未知参数反演

假定三个排放特征参数（排放节点、排放量、排放时刻）均未知，需要经过溯源模型提供信息。污染物从节点瞬时排入管网，排放量为M，排放时刻为模拟开始后T，排放节点参数Jx。假定某日在T时，在节点Jx处的某工厂以M（g/min）的速度连续排放污染物BX。

按照MCMC反演思路迭代4000次。反演参数M（g/min）、Jx和T（min）的后验概率直方图分别见图4.84、图4.85、图4.86，参数M（g/min）落在各区域的频数统计见表4.55，参数T（min）落在各区域的频数统计见表4.56，参数Jx落在各区域的频数统计见表4.57。

85

重庆大学硕士学位论文

图4.84 参数M（g/min）的后验概率直方图

Fig4.84 The posterior probability histogram of parameter M（g/min）

表4.55 参数M（g/min）落在各区域的频数统计

Table4.55 Frequency statistics of parameter M（g/min） that falling in each region

区域（×106g）

M（g/min）

频数

800

1700

1500

0.7～0.9

0.9～1.1

1.1～1.4

图4.85 参数T（min）的后验概率直方图

Fig4.85 The posterior probability histogram of parameter T（min）

表4.56 参数T（min）落在各区域的频数统计

Table4.56 Frequency statistics of parameter T（min）that falling in each region

区域（min）

T（min）

频数

1900

2100

0～40

40～80

图4.86 参数Jx的后验概率直方图

Fig4.86 The posterior probability histogram of parameter Jx

86

4 排水管道污染物溯源模型的验证与应用

表4.57 参数Jx落在各区域的频数统计

Table4.57 Frequency statistics of parameter Jx that falling in each region

区域（x）

Jx

频数

1400

2000

600

0～20

20～40

40～60

从图4.84、4.85、4.86可以看出，污染物连续排放的反演参数排放量M（g/min）、初始排放时间T（min）、排放节点Jx的后验概率密度均较离散，只能大概推测出污染事件排放量在M=[80000g/min，137000g/min]之间，初始排放时间大概发生在T=[0min，70min]之间，排放节点大概在Jx=[5，55]之间。

如表4.58所示，通过将连续污染事件的反演参数排放量M（g/min）、初始排放时间T（min）的抽样结果与真实值对比，可知污染物排放的初始排放时间T（min）的数学期望和中值误差均很大，排放量M（g/min）的均值误差和中值误差都比较小，但其后验概率密度比较分散， 污染源参数T（min）、M（g/min）的反演结果误差见图4.87。如图4.88所示，污染源Jx参数的后验统计结果分布比较广泛。这也表明如果两个潜在污染源到达水质监测点的距离相当，那么模型结果会受到干扰，这种情形下，需要通过更精确的先验信息来提高溯源的精确度，也即是现场监测数据来缩小先验信息的取值区间。

表4.58 污染源参数T（min）、M（g/min）的抽样结果与实际值对比

Table4.58 Comparison of sampling results and actual values of inversion parameters T（min）and

M（g/min）

参数 M（g/min） T（min）

实际值 100000 60

403020100T（min）中值误差/% M（g/min）均值误差/%中值 109543 41

中值误差/%

9.54 31.67

均值 108765 40.9722

均值误差/%

8.76 31.71

31.6731.719.548.76 表4.87 污染源参数T（min）、M（g/min）的反演结果误差 Fig4.87 Inversion result error of parameters T（min）and M（g/min）

87

重庆大学硕士学位论文

图4.88 污染源参数Jx的后验统计结果

Fig4.88 The posterior statistical results of pollution source parameters Jx

4.4 溯源结果分析与讨论

从上述分别对污水管网污染物节点瞬时排放及连续排放，污染源源项单个未知参数、两个未知参数以及三未知参数的贝叶斯MCMC抽样反演结果中，可以得到以下结论：

① 基于贝叶斯统计马尔科夫蒙特卡罗(MCMC)方法的城市污水管网污染物排放，污染源源项溯源反演算法可以有效地反演未知的污染源的排放特征信息。

② 通过对污染物瞬时排放单变量反演参数不同游走步长的比较，分析得出在本算法中污染物瞬时排放量M（g）游走步长为10000，排放节点Jx步长为1，排放时间T（min）步长为1时，可以提高反演精度和反演效率。

③ 污染物瞬时排放：污染源单个未知参数反演中，贝叶斯MCMC抽样算法可以很准确的识别污染源，具有反演精度高、误差小的特点，尤其在对排放时间反演更加准确。两个未知参数反演中，对排放量M（g）和排放时间T（min）识别比较准确，反演精度较高，误差较小，但是对排放节点Jx反演存在较大的误差。三个未知参数反演中，对排放量M（g）和排放时间（min）的识别比较准确，反演精度较高，误差在10%左右，对排放节点Jx的反演误差在25%左右，无法准确识别只能大概推断其范围。

④ 污染物连续排放：污染源单个未知参数反演中，排放量M（g/min）的反演精度高、误差在0.35%左右，排放时间T（min）和排放节点Jx的误差也较小在10%内。两个未知参数的反演中，排放量M（g）反演误差较小，排放时间T（min）

88

4 排水管道污染物溯源模型的验证与应用

和排放节点Jx的反演误差都比较大。三个未知参数反演中，排放量M（g）和排放节点反演误差在10%内，但是分布范围都比较广，排放时间T（min）误差较大，所以三个未知参数的反演结果只能推断其大概发生的范围。

根据污染物质的迁移转换控制方程Ct=C0e

−k

x

u

可知，导致双个未知参数和三

未知参数反演误差增大的原因是反演参数的多组合性，由于污水管网各管段长度及水流流速不同等原因，使不是真实值的参数组合到一起也可以得到同样或相近的结果。而贝叶斯统计溯源法的目的就是找出可能的值，并给出所有值的概率，其基本思想是缩小先验范围，着重检查概率分布大的点，但也不能忽略概率小的点。因此，贝叶斯方法溯源的结果是合理的，且有实际指导意义。这也是统计法溯源不同于确定性的优化法溯源的一大特色。

4.5 本章小结

本章在实际排水管网的基础上，通过数值试验生成理论解，并进行扰动后作为反演的基础数据，将贝叶斯反演结果与实际监测水质数据进行比较，检验了基于贝叶斯统计法-SWMM溯源模型的有效性和准确性，并得出主要贝叶斯反演参数的取值。

研究考察了污水管网节点瞬时排放污染物及节点连续排放污染物两种排放模式。同时为考察贝叶斯反演算法对于信息缺失度的敏感程度，在每种排放模式下，考察3种反演的工况：单个未知参数，两个未知参数和三个未知参数，分别代表排放节点、排放量、排放时间三个排放特征参数中有单个、两个、三个未知参数需要反推，其他参数已知的情形，比较分析反演参数的个数对贝叶斯算法有效性的影响。

开展了不同流态下反演参数的个数对贝叶斯算法可靠性的影响研究，研究结果表明：

① 通过对污染物瞬时排放时单参数选取不同游走步长反演结果的比较，分析得出在本算法中污染物瞬时排放量M（g）游走步长为10000，连续排放量M（g/min）步长为1000，排放节点Jx步长为1以及排放时间T（min）步长为1时，可以提高反演精度和反演效率。

② 污染物瞬时排放：污染源单个未知参数反演中，贝叶斯MCMC抽样算法可以很准确的识别污染源，具有反演精度高、误差小的特点，尤其在对排放时间反演更加准确；两个未知参数反演中，对排放量M（g）和排放时间T（min）的识别比较准确，反演精度较高，误差较小，对排放节点Jx的反演可将其范围缩小到四分之一到三分之一；三个未知参数反演中，对排放量M（g）和排放时间（min）的识别比较准确，反演精度较高，误差在10%左右，对排放节点Jx的反演可将范

89

重庆大学硕士学位论文

围缩小到三分之一。

③ 污染物连续排放：污染源单个未知参数反演中，排放量M（g/min）的反演精度高、误差在0.35%左右，排放时间T（min）反演的误差也较小在10%内，排放节点Jx范围缩小到五分之一；两个未知参数的反演中，对排放量M（g）反演误差较小，排放时间T（min）演误差都比较大，排放节点Jx的反演范围分布比较广；三个未知参数的反演中，对排放量M（g）和排放节点反演误差在10%内，但是分布范围都比较广，对排放时间T（min）误差较大，排放节点Jx的反演范围分布比较广，所以三参数的反演结果只能推断其大概发生的范围。根据𝐶𝑡=𝐶0𝑒

−𝑘

𝑥

𝑢

可知，导致两个参数和三个参数反演误差增大的原因是反演参数的多组合性，由于污水管网各管段长度及水流流速不同等原因，使不是真实值的参数组合到一起也可以得到同样或相近的结果。

利用贝叶斯推理方法进行污染源的反演能够充分利用参数的先验信息，综合考虑观测误差和模型误差，反演结果利用后验概率分布表示，能够给出具有实际指导意义的污染源范围及排放特征。

90

5 结论与建议

5 结论与建议

5.1 结论

本研究针对城市排水系统中频繁发生的工厂偷排漏排超标废水对污水处理厂活性污泥法工艺造成严重冲击的污染事件，开展了基于在线水质监测数据的排水管道偷排漏排污染物质溯源模型的研究，构建了基于贝叶斯统计推理算法、SWMM水力水质模型以及Matlab编程的排水管网污染物质溯源数学模型，对偷排超标污染物质的排放位置、排放量和排放时段等三个未知的排放特征参数进行了统计反推，获得了排放节点、排放量和排放时段的概率分布；并重点研究了统计反演算法中游走步长的取值、采样监测的时间间隔及在线监测点的空间布置等实际限制条件对反演结果精度的影响；反演模型采用了先进的马尔科夫蒙特卡罗（MCMC）抽样算法，大大提高了溯源的效率，在有限的反馈时间内可提供更精确的溯源结果，从而提高对偷排漏排超标废水的管理水平以及排放事故的应急处理效率。

主要研究内容和结论如下：

① 开展了SWMM水质模型、贝叶斯统计方法及MATLAB软件的耦合集成研究。该集成通过MATLAB将SWMM模拟结果的时间序列输出，与相应的实际监测数据序列进行匹配提高了分析精度，实现了在线监测动态反馈，与实际管网非运行工况更接近。MATLAB强大的数据处理功能，能够实现MCMC抽样算法涉及的大量计算，提高了MCMC反演效率，从而快速确定出污水偷排漏排可能的排放节点、排放时间和排放量。

② 贝叶斯统计MCMC抽样算法中游走步长与反演效率和精度的关联关系研究。结果表明：步长过小会使抽样限于局部解，大的步长可以游走更大的搜索空间，避免限于局部解，有效反演污染源参数的后验概率密度；较大的步长可以实现对整个后验空间的抽样，但是步长过大会使抽样次数增加，反演时间延长，效率降低。本研究提出的基于贝叶斯统计马尔科夫蒙特卡罗(MCMC)方法的城市污水管网污染物排放，污染源源项溯源反演算法，通过对污染物瞬时排放单变量反演参数不同游走步长的比较，分析得出在本算法中污染物瞬时排放量M（g）游走步长为10000，排放节点Jx步长为1，排放时间T（min）步长为1时，可以提高反演精度和反演效率。

③ 采样监测的时间间隔及在线监测点的空间布置对反演算法有效性和精度的影响研究。结果表明：监测点的观测数据过少时，无法准确反演污染源参数的后验概率密度；当观测数据增多、采样监测时间间隔比较小时，观测数据没有反映监测点的浓度序列，此时也无法准确反演污染源参数的后验概率密度；只有适

91

重庆大学硕士学位论文

当的采样监测时间间隔才能使反演结果的精度和效率更高；水质监测点增加时可提高溯源效率。

④ 通过实际管网案例开展了数值试验，将贝叶斯反演结果与实际监测水质数据进行比较，检验了基于贝叶斯统计法-SWMM溯源模型的有效性和准确性，并得出主要贝叶斯反演参数的取值。研究了两种污染物质的排放模式：

1）污染物瞬时排放：污染源单个未知参数反演中，贝叶斯MCMC抽样算法可以很准确的识别污染源，具有反演精度高、误差小的特点，尤其在对排放时间反演更加准确。两个未知参数反演中，对排放量M（g）和排放时间T（min）识别比较准确，反演精度较高，误差较小，但是对排放节点Jx反演存在较大的误差。三个未知参数反演中，对排放量M（g）和排放时间（min）的识别比较准确，反演精度较高，误差在10%左右，对排放节点Jx的反演误差在25%左右，无法准确识别只能大概推断其范围。

2）污染物连续排放：污染源单个未知参数反演中，排放量M（g/min）的反演精度高、误差在0.35%左右，排放时间T（min）和排放节点Jx的误差也较小在10%内。两个未知参数的反演中，排放量M（g）反演误差较小，排放时间T（min）和排放节点Jx的反演误差都比较大。三个未知参数反演中，排放量M（g）和排放节点反演误差在10%内，但是分布范围都比较广，排放时间T（min）误差较大，所以三个未知参数的反演结果只能推断其大概发生的范围。

总的来说，构建的溯源模型可模拟污水管道系统中保守及非保守污染物质的排放、迁移、转输全过程，并基于贝叶斯理论对排水管网中可能的偷排漏排污染点源进行统计意义上的反推并给出各个可能取值区间的概率，从而筛选出可能的污水排放节点、排放量和排放时间。通过贝叶斯方法、Matlab与SWMM的耦合集成，水质预测结果与实际管网的动态运行工况更接近，提高污染物传输模拟的精度，从而提高反演分析的精度。同时贝叶斯统计反演方法充分考虑了先验信息、观测噪声以及模型误差的影响，在泄漏源反算的不确定性分析中独具优势，能较好应对排水管网系统溯源反问题的不确定性，提供有实际指导意义的溯源结果。论文的研究成果可为水体污染控制、智慧排水监管系统等领域提供技术支持。

5.2 建议

本研究针对污水管网突发污染事件溯源，构建了贝叶斯统计MCMC算法溯源模型，虽然取得了一些成果，但仍然存在许多不足之处，还需要针对以下问题继续研究：

① 目前溯源模型还处于比较初级阶段，适用范围比较窄，只能用于水文条件稳定的污水管道，在应对水文条件复杂的污水管道溯源结果准确度无法保证。为

92

5 结论与建议

了能够在真实的污水管网突发污染事件中也能有效的反演出每个污染源参数，可以引入污水管道分段建模的思想，根据水文条件的不同进行分段建模，进一步提升模型的适用性。

② 本研究研究的污水管道污染物溯源模型还未应用到实际的水质预警管理系统中，需要不断改进模型建模和模型优化，将算法编写成易于调用的接口，集成到实际工程系统中，为水体污染控制、智慧排水监管系统等领域提供技术支持。

93

重庆大学硕士学位论文

94

致 谢

致 谢

青春如梦，岁月如花，流水似年，稍纵即逝。转眼间，三年的研究生生涯已接近尾声，重大带给我诸多成长与收获，这将是我一生宝贵的财富，她的春夏秋冬、一草一木将永远铭记在我的记忆里。

本论文从论文选题、论文开题、文献研读、系统开发、论文初稿、到定稿历时18个月完成。在此过程中，离不开我的导师向平副教授、柴宏祥教授和邵知宇老师的言传身教和悉心指导。导师渊博的知识，丰富的经历，严谨的治学态度，高屋建瓴的科学思维和勇于开拓的精神使我受益匪浅；此外，感谢柴老师和邵老师对我生活的种种关心，是您们让我明白了更多的为人之道，是您们让我懂得了研究生生活的意义。在论文付梓之际，谨向恩师致以最崇高的敬意和衷心的感谢！

同样，我要感谢重庆大学城市建设与环境工程学院何强老师课题组的各位老师们，感谢他们对我的指导、关心和帮助。感谢与我一同走过三年研究生时光的向钰和裔士刚，感谢她们在这一路上的相互陪伴和帮助，让我们共同进步和成长。感谢我的同门以及师弟师妹在我在学习上、生活上、工作上给予的支持与帮助，谢谢你们让我的研究生生活充满精彩；特别感谢李霜、张晓媛在课题上给予我的无私帮助，课题能顺利进行和完成你们的功劳不可或缺。

特别需要感谢我的家人，感谢你们的养育之恩，感谢你们在我成长路上的默默支持，感谢你们对我的信任，感谢你们带给我的奋斗的决心和勇气，我会不忘初心，勇往直前。

最后，衷心地感谢在百忙之中评阅论文以及参加答辩的各位专家、教授和答辩秘书，感谢你们对论文提出宝贵的意见！

郑卓乐

二O一八年五月

95

重庆大学硕士学位论文

96

参考文献

参考文献

[1] 张风宝，杨明义，赵晓光，等．磁性示踪在土壤侵蚀研究中的应用进展叨．地球科学进

展．2005(07)：75 1-756．

[2] 方晶晶. 河北平原邯郸地区地下水硝酸盐污染来源及迁移转化过程的多元素同位素及微

生物（E.coli）示踪[D].中国地质大学,2014.

[3] Lepot M, Makris K F, Fhlr C. Detection and quantification of lateral, illicit connections and

infiltration in sewers with Infra-Red camera: Conclusions after a wide experimental plan[J]. Water Research, 2017, 122:678.

[4] Farrell J A，Pang S，LiW．Chemical PlumeTracingvia all Autonomous Underwater

Vehicle[J]．IEEE Joumal of Oceanic Engineering．2005，30(2)：428442． [5] 郑伟. 城市污水管网有毒物质溯源监控技术研究[D]. 重庆大学, 2011.

[6] U.S.EPA. Exposure Models Library and Integrated model evaluation system,

EPA/600/C-92/002 [R].1996.

[7] Jobson H E. Predicting travel time and dispersion in river and streams [J]. Journal of Hydraulic

Engineering, 1997, 123(11):971-977.

[8] Cammara A S, Randall C W.The qual II model[J]. Journal of Environmental Engineering,

ASCE, 1984, 110(5):993.

of predictive models for surface water quality and [9] Thomann R V. The future ̓Golden Age

ecosystem management[J].Journal of Environmental Engineering, 1998, 124(2):94-103. [10] 徐祖信, 廖振良. 水质数学模型研究的发展阶段和空间层次[J]. 上海环境科学， 2003，

22(2):79-85.

[11] 李本刚,陶澎,曹军.水环境模型与水环境模型库管理[J]. 水科学进展, 2002, 13(1):14-12. [12] Ashley, R. M. And Crabtree, R.W. Sediment origins deposition and build-up in combined sewer

systems [J]. Water Science and Technology, 1992, 33(8):1-12.

[13] Matos, I.S., Aires, C.M.. Mathematical modeling of sulphides and hydrogen sulphide gas

build-up in the Costa do Estoril sewerage system [J]. Water Science and Technology, 1995, 31(7):255-261.

[14] Vollertsen J., Hvitved-Jacobsen T.. Sewer quality modeling-a dry weather approach [J]. Urban

Water, 2000, 2(4):295-303.

[15] 冯良记. 城市排水管网水动力及水质转化模型研究[D]. 大连：大连理工大学，2009. [16] 白子易. 下水道系统生化动力模式建立之研究[D]. 中坜: 国立中央大学, 2001. [17] 朱学兵. 城市排水管道数学模型及模拟[D]. 西安: 西安建筑科技大学,2004.

97

重庆大学硕士学位论文

[18] 周永潮.水力模型在城市排水系统设计与管理中的应用研究[D].同济大学;同济大学环境科

学与工程学院,2008.

[19] 杨雪.基于SWMM模型的合流制溢流污染研究及控制[D].北京建筑工程学院,2008. [20] 黄兵,朱晓敏,王树东,郑金龙,姚波,陈素. 基于SWMM的昆明市船房片区合流制排水系统

模拟[J]. 中国给水排水,2012,28(19):15-18+23.

[21] 孙全民,胡湛波,李志华,彭艺艺,陈立华,魏群.基于SWMM截流式合流制管网溢流水质水量

模拟[J].给水排水,2010,46(07):175-179.

[22] 肖庭延, 于慎根, 王彦飞.反问题的数值解法[ M] .北京:科学出版社, 2003:1-2.

[23] VICTOR I, STEFAN K.Identification of the diffusion coefficient in a one-dimensional

parabolic equation[J] .Inverse Problems, 2000 , 16(3):665-680.

[24] BARBAR K T.Regularization by projection with a posterior discretization level choice for level

choice for linear and nonlinear ill-posed prob-lem[ J] .Inverse Problem s, 2000 , 16(6):1523-1539.

[25] 闵涛, 周孝德.河流水质纵向弥散系数反问题的迭代算法[ J] .水动力学研究与进展:A 辑,

2003 , 18(5):547-552.

[26] 栗苏文, 李兰, 李志永.Dobbins 模型参数识别反问题的导数-正则化方法[J] .水利学报,

2004(7):104-108.

[27] 王佳鹤，金忠青,等.对流－扩散方程反问题的控制论求解方法[J].水动力学研究与进展A

辑.1996(2)

[28] 闵涛,周孝德,张世梅,等.对流-扩散方程源项识别反问题的遗传算法[J].水动力学研究与进

展,2004.19(4):520-524

[29] Preis A, Ostfeld A. Genetic algorithm for contaminant source characterization using imperfect

sensors[J]. Civil Engineering and Environmental Systems.2008, 25(1): 29-39

[30] 陶涛, 芦颖军, 信昆仑，等.基于模式识别的供水管网污染源追踪分析方法[J].天津大学学

报：自然科学与工程技术版》2013年 第5期

[31] Preis A, Ostfeld A. Contamination source identification in water systems: a hybrid model

trees-linear programming scheme[J]. Journal of Water Resources Planning and Management - ASCE, 2006, 132(4): 263-273

[32] 吕谋,王梦琳,孙贤忠，供水管网突发污染试验模拟及污染源定位研究[J]，青岛理工大学学

报. 2009.30（6）：1-17

[33] KimM., Choi C.Y., Gerba C.P., Source tracking of microbialintrusion in water systems using

artificialneural networks[J].Water Research. 2008.42(4):1308-1314

[34] 李兰.水质反问题模型的时域频域算法[ J] .水科学进展, 1998 , 9(3):218-223. [35] 李兰.BOD5-DO 参数反问题偶合模型的研究[ J] .水科学进展, 2000 , 11(3):255-259.

98

参考文献

[36] 宁莎莎,蒯琳萍.混合遗传算法在核事故源项反演中的应用[J].原子能科学技

术,2012,46(S1):469-472.

[37] Ng, A. W.M., Perera, B.J.C. Selection of genetic algorithm operators for river water quality

model calibration[J]. Engineeering Applications of Artificial Intelligen,16(5-6),529-541. [38] Banik，Optimal placement of water quality monitoring stations in sewer[J]，Procedia

Engineering，2015，（119）：1308–1317

[39] Subhankar Karmakar，P.P. Mujumdar. Grey fuzzy optimization model for water quality

management of a river system. Advances in Water Resources,29,1088-1105.

[40] 冯耀龙,韩文秀.河流系统水质管理模糊优化模型[J].天津大学学报（自然科学与工程技术

版）,2001,(5):698-701. DOI:10.3969/j.issn.0493-2137.2001.05.031.

[41] 杨洋 .基于贪心策略的自适应关键帧提取算法研究:[河北工业大学硕士学位论文]. 天津:

河北工业大学,2011.

[42] 陶涛, 吕存阵, 信昆仑,等. 基于突发污染事件的管网水质监测点优化布置[J]. 同济大学学

报(自然科学版), 2010, 38(11):1621-1625.

[43] Zou，R，Lung，W.S，Guo，H. Neural network embedded Monte Carlo approach for water quality

modelling under input information uncertainty.Journal of Computing in Civil Engineering,2002,16（2），135-142.

[44] 信昆仑,项宁银,陶涛,尹兆龙.基于用户水质投诉信息的供水管网污染源追踪定位[J].同济大

学学报(自然科学版),2014,42(02):283-286+291.

[45] Kim M,Suwon G,Choi C Y.; Gerba, C P. Development and evaluation of a

decision-supporting model for identifying the source location of microbial intrusions in real gravity sewer systems[J]，Water Research. 2013,47(13):4630-4638

[46] 杨一帆,张凯山.突发型大气污染源位置识别反演问题的数值模拟[J].环境科学学

报,2013,33(09):2388-2394.

[47] Andrieu C，Freitas ND，Doucet A，et a1．An Introduction to MCMC for Machine

Learning[J]．Machine Leaming．2003，50(1-2)：5-43．

[48] Xu Z, Wang L, Yin H, et al. Source apportionment of non-storm water entries into storm drains

using marker species: Modeling approach and verification[J]. Ecological Indicators, 2016, 61:546-557.

[49] 曹小群，宋军强，张卫民，等．对流-扩散方程源项识别反问题的 MCMC方［J］．水动

力学研究与进展: A 辑，2010，25(2):127-136．

[50] 陈海洋，腾彦果， 王金生， 等． 基于 Bayesian-MCMC 方法的水体污染识别反问题

［J］．湖南大学学报:自然科学版， 2012，39( 6) : 74-78．

[51] 姜继平,董芙嘉,刘仁涛,袁一星.基于河流示踪实验的Bayes污染溯源:算法参数、影响因素

99

重庆大学硕士学位论文

及频率法对比[J].中国环境科学,2017,37(10):3813-3825.

[52] 王晓波,李诗赞,代雪云,张楠.基于事故树与贝叶斯网络的管道泄漏事故溯源方法[J].油气储

运,2017,36(09):1013-1018.

[53] 冯民权,张园园.基于MATLAB的贝叶斯网络供水工程水质风险分析[J].黑龙江大学工程学

报,2015,6(01):5-11.

[54] 朱嵩. 基于贝叶斯推理的环境水力学反问题研究[D]. 浙江大学，2008.

[55] Voutilainen A，Kaipio J P. Sequential Monte Carlo estimation of aerosol size distributions[J].

Computational Statistics & Data Analysis the Official Journal of the International Association for Statistical Computing. 2004，48(4): 887-908.

[56] 李明亮. 基于贝叶斯统计的水文模型不确定性研究[D].清华大学,2012. [57] 李勇,易文德. 贝叶斯分析中先验分布优选的方法[J].自然科学版,2005(04):5-7. [58] 阙慧敏. 基于先验信息的统计预测方法及其应用研究[D].华北电力大学（北京）,2011. [59] 赵艳菊. 强噪声背景下机械设备微弱信号的提取与检测技术研究[D].天津大学,2009. [60] 陈曦晖. 强噪声干扰下行星轮系振动信号分析及其故障诊断技术研究[D].中国矿业大

学,2017.

[61] 电子部第四十一所.AV3981 职能噪声系数分析仪说明书[K].2000.

[62] Andrieu C，Freitas N D，Doucet A.et al.An introduction to MCMC for machine learning [J].

Machine Learning，2003，50：5—43.

100

附 录

附 录

A. 作者在攻读硕士学位期间发表的论文

[1] 邵知宇，郑卓乐，黄文钟，李霜，张晓媛，柴宏祥. 多功能雨水塘水位-流量曲线的构建

及水力学原理研究,中国给水排水, 2017,33(9).

B. 作者在攻读硕士学位期间参加的相关科研项目

[1] 国家“水体污染控制与治理”科技重大专项:绿色建筑与小区低影响开发雨水系统研究与示

范(2010ZX07320-001).

C. 作者在攻读硕士学位期间申报与获批的专利

[1] 邵知宇，柴宏祥，张晓媛，郑卓乐，何强. 发明专利：一种识别雨水管网污水直排污染源

的方法.专利号：201810114905.X.

101