CN108897964A - 一种污水管网超标排放工业废水的贝叶斯统计溯源方法 - Google Patents

Abstract

本发明公开了一种污水管网超标排放工业废水的贝叶斯统计溯源方法，它包括：1、在未知参数先验信息范围内随机产生初始点；2、模拟出当前参数对应监测点污染物浓度时间序列，与实际监测数据比较来获得未知参数的后验概率密度；3、根据建议分布产生候选参数，模拟出对应监测点污染物浓度时间序列，与实际监测数据比较似然程度来获得未知参数的后验概率密度，4、抽取一个随机数，判断该候选值是否接受，输出接受值及后验概率密度；5、重复步骤3和4至迭代完成。本发明的优点是：有效地缩小未知参数的取值范围，利用MCMC抽样方法特点，在保证抽样合理性的条件下，减少工作量及抽样时间，提高溯源效率。

Description

一种污水管网超标排放工业废水的贝叶斯统计溯源方法

技术领域

本发明属于市政工程信息技术，具体涉及一种集成贝叶斯统计推理算法、SWMM模型以及Matlab编程的污水管网中超标排放工业废水溯源的方法。

背景技术

在污水管网中，含有高浓度的重金属及其他有毒物质的偷排工业废水通常对污水处理厂活性污泥法工艺造成严重冲击，甚至导致活性污泥中毒死亡，出水水质不达标。而其中的溯源问题具有高度不确定性，常规数学模型溯源方法不适宜，通常溯源反推的结果为唯一值，一旦该反推结果与实际偏离，溯源结果便失去实际指导意义。

中国专利申请号为201810114905.X，名称为一种识别雨水管网污水直排污染源的方法,它包括：1、生成运行SWMM软件所需的初始输入*.inp文件；2、在Matlab中定义SWMM输入参数管网中污水直排节点编号、排放浓度、排放量三个参数的取值范围、取样步长和样本数；3、在Matlab中采用三层循环对所有样本组合进行取样，更新*.inp文件中相应位置的参数信息；4、在Matlab中用命令“！swmm5.exe<input file><report file>”调用Matlab外部的swmm5.exe程序，利用更新后的*.inp文件中的信息进行计算，并将运算结果储存在*.rpt文件中；5、计算第i次取样时的似然函数值；6、按似然函数值从大到小对参数组合进行排序。

1.该专利的推求未知参数的方法为枚举法，即对三个未知参数所有可能的取值逐一进行计算评估，参数的候选过程没有直接利用实际监测数据信息进行反馈，下一个参数候选值的选取完全是一随机过程，当前参数取值对下一候选值不产生影响。

2.该专利抽样计算工作量大，时间长。假设管网需反推的总节点数为NX，进行3个未知污染源信息参数的溯源计算，该专利中的枚举法需要正向求解NX³次污染物质对流传输方程，计算量大、反推时间长，溯源反馈慢。

发明内容

本发明所要解决的技术问题是提供一种污水管网超标排放工业废水的贝叶斯统计溯源方法。它利用贝叶斯理论MCMC抽样方法，能够在抽样过程中及时的反映样本的概率分布特点，选取后验分布概率大的区域进行集中抽样，对于较小区域进行较少抽样，减少在整个样本空间进行搜索的工作量，反映由于污水管网的复杂水力情况带来的溯源不确定性，有效地筛选出污水管网中超标废水偷排的位置、排放量、排放时间取值，辅助现场物理溯源，减少物理溯源的工作量，缩短污染事故反馈时间。

本发明所要解决的技术问题是通过这样的技术方案实现的，它包括以下步骤：

步骤1：在SWMM软件中构建污水管网模型，生成初始*.inp文件，然后在Matlab中打开初始*.inp文件，对*.inp文件进行读写，对所有可能偷排的节点、排放量和排放时间三个参数进行编辑，外部调用swmm5.exe程序，直接利用更新后的参数输入运行SWMM，将包括监测点的污染物浓度时间序列运行结果储存在*.rpt文件中；

步骤2：在Matlab中实现对*.rpt文件的读写，对排放节点、排放量和排放时间三个未知参数在先验信息范围内进行污染物溯源统计MCMC抽样，每次抽样后，需同步更新*.inp文件中相应的参数变量，外部运行SWMM，将更新的节点污染物浓度时间序列储存在*.rpt文件中；

步骤3：在Matlab中提取*.rpt文件中相应的模拟结果的时间序列，结合实际监测节点处的污染物质浓度时间序列，利用贝叶斯理论反演计算，考虑监测过程及模拟过程的误差，计算似然函数值和待反推的未知参数的后验概率分布；

步骤4，在Matlab中输出后验概率直方图，进行贝叶斯统计溯源结果的可视化分析，推求未知的偷排节点、排放量和排放时间的取值的概率分布。

本发明的技术效果是：

1、将污水管网系统中超标工业废水偷排的溯源问题转化为数学问题，采用贝叶斯MCMC统计反演方法，实现对大量可能的排放特征值进行筛选，得实际意义的排放信息的集合，利用了SWMM强大计算功能，实现管网复杂水力条件下污染物的传输模拟，与实际管网的动态运行工况更接近，提高了预测溯源的精度。

2、采用了贝叶斯方法进行统计反演，充分考虑了先验信息、测量误差及模型误差的影响，结果采用未知参数概率分布的形式，避免了传统求解采用固定唯一值来表征求解结果所导致的求解偏离结果便失去意义的极端情况，在高度不确定排水管网系统的溯源反问题求解中独具优势。

3、利用Matlab强大的数据处理功能，贝叶斯统计溯源根据当前参数取值结果与实际监测结果的吻合程度决定下一个参数的选取方向，能够进行有目的性的搜索，保证参数选取过程花更多的时间在最重要的区域，产生的样本更能够模仿目标分布样本，在保证抽样合理性的条件下，降低一般随机抽样所需要的工作量，减少抽样时间，提高溯源效率，从而快速筛选出超标工业废水偷排可能的节点、排放浓度和排放量。

本发明利用城市排水模型SWMM对污染物质在排水管网中的迁移和传输转化过程进行精确预测，基于贝叶斯理论对超标污染物排放节点、排放量及排放时间三个未知参数进行统计法反推，能够较好解决排水管网溯源问题的不确定性，面对复杂的多节点管网溯源问题，能够依据检测样本特点，降低抽样工作量，减少抽样计算时间，更加及时的应对超标废水偷排情况。

附图说明

本发明的附图说明如下：

图1：污染物统计溯源MCMC抽样流程图；

图2：实例的区域污水管网平面布置图；

图3：检测浓度序列图

图4：反演参数M(g)的后验概率直方图；

图5：反演参数T(min)的后验概率直方图；

图6：反演参数Jx的后验概率直方图；

具体实施方式

下面结合附图和实施例对本发明做进一步说明：

本发明的构思是：联合Matlab编程软件、SWMM排水模型及贝叶斯统计推理，对污水管网中可能的超标工业废水偷排污染源进行高效统计筛选。利用贝叶斯MCMC的特点，在保证反推结果精度的前提下，减少工作量，缩短溯源时间，为水环境保护提供强有力的帮助。

本发明利用城市排水模型SWMM对污染物质在排水管网中的迁移和传输转化过程进行精确预测，基于贝叶斯理论对超标污染物排放节点、排放量及排放时间三个未知参数进行统计法反推，包括以下步骤：

步骤1、在SWMM软件中构建污水管网模型，生成初始*.inp文件，后在Matlab中打开初始*.inp文件，对*.inp文件进行读写，对所有可能偷排的节点、排放量和排放时间三个参数进行编辑，外部调用swmm5.exe程序，直接利用更新后的参数输入运行SWMM，将包括监测点的污染物浓度时间序列运行结果储存在*.rpt文件中；

基础数据SWMM模型中直接添加并定义属性，或者外部文件输入的方式。先验概率是指根据以往经验和分析得到的概率。后验概率是指在得到“结果”的信息后重新修正的概率。在排水管网污染物质溯源问题中，需反推超标工业废水的排放节点、排放量和排放时间三个未知参数，这些参数具有随机性，但都分布在一个取值范围内，这个范围可以根据经验或主观判断进行确定，即是获得一个较窄的先验信息范围(可以根据经验或主观判断进行确定的未知参数范围)，缩小后续反演的搜索范围。同时认为排放量和排放时间的取值服从取值空间内的连续型均匀分布，因此它们的先验概率密度函数为：

排放节点X₃服从其取值范围的离散点匀分布，其概率函数为：

设某节点为排放节点情况下的先验总概率密度函数p(X)为：

式中p(X_f)代表取值区间[a_f,b_f]上的均匀分布；f表示第f个参数；H代表需要反推的参数的个数；X_f代表第f种未知参数；п为连乘符号；K代表第K个排放节点；m代表共有m个排放节点；a、b分别为参数取值的上限和下限。

步骤2、在Matlab中实现对*.rpt文件的读写，对排放节点、排放量和排放时间三个未知参数在先验信息范围内进行污染物溯源统计MCMC抽样，每次抽样后，需同步更新*.inp文件中相应的参数变量，外部运行SWMM，将更新的节点污染物浓度时间序列储存在*.rpt文件中。

MCMC抽样是根据无明确数学表达式的后验概率分布进行抽样，从而缩小抽样空间，减少先验信息范围大而带来的工作量，产生大量服从后验概率分布的随机向量序列，生成一条候选值依赖于当前参数值(当前参数指第t次抽样样本X_t)的马尔科夫链，使得抽样过程集中在在后验概率分布大的区域，在这些区域密集抽样，而同时兼顾概率分布较小的区域，使得产生的样本更能够模仿目标分布样本而非绝对的唯一值。MCMC抽样样本为先验信息范围内的排放节点、排放量、排放时间的不同组合。

步骤3、在Matlab中提取*.rpt文件中相应的模拟结果的时间序列，结合实际监测节点处的污染物质浓度时间序列数据，利用贝叶斯理论反演计算，考虑监测过程及模拟过程的误差，计算似然函数值和待反推的未知参数的后验概率分布。

说明：步骤3嵌套在步骤2中的，为了强调步骤3这步的操作功能，把它单独列出来写，如果以步骤2、步骤3这个抽样的过程作为整个大流程来写，则步骤2会显得冗长。

上述步骤2和步骤3通过如图1所示的污染物溯源统计MCMC抽样流程图来实施，该流程开始于步骤S1；

在步骤S2，在未知参数的先验信息范围内随机抽取样本产生初始点X₀；

根据专家经验或主观判断，确定参数的取值范围，得到一个较窄的先验信息范围。如节点数不可能超过实际节点数。排放量不可能低于0，也不可能无限制的高。

在步骤S3，令t＝0，当前参数X_t＝X₀，利用SWMM一维污染物质对流传输模型计算出当前参数X_t对应的监测点处污染物浓度值，并与实际监测数据比较似然程度，计算出该参数对应的后验概率密度p(X_t|y)，输出X_t及p(X_t|y)。

利用贝叶斯理论反演计算，结合实际监测节点监测的污染物浓度时间序列数据，构造表征污染物扩散模型模拟的污染物浓度理论值和实测数据拟合程度的似然函数p(y|X_t)，充分实现利用已知信息进行未知参数反推：

式中：n为测量数据个数；排放强度为M；排放节点为Jx；排放时间为t；已知监测点观测到污染物浓度观测值Y＝{Y₁，Y₂，…，Y_i，…，Y_n}；污染物扩散模型模拟的污染物浓度理论值N＝{N₁，N₂，…，N_i，…，N_n}；为污染事件发生后在监测点第i个节点观测到的污染物浓度值Y_i；在未知参数X_t为M，Jx，T的条件下，通过污染物扩散模型模拟的监测点第i个节点处污染物浓度理论值为N_i(M，Jx，T|X_t)；δ是似然函数标准差。

根据贝叶斯理论，后验概率＝标准似然度*先验概率，可以计算出未知参数X_t的后验概率密度p(X_t|y)，代表了不同参数值的可能性，也即是未知参数X_t统计意义下的解：

式中，α为比例常数；其它参数意义同上。

在步骤S4，假定新的参数候选值X^*服从以当前参数X_t为均值，方差为σ的正态分布(建议分布)，根据建议分布在当前参数X_t状态下在先验信息范围内抽取X^*，利用水质模型计算出候选参数X^*对应的污染物浓度时间序列，与实际监测数据比较似然程度，按照式(5)计算得到候选参数X^*的后验概率密度P(X^*|y)；

从后验分布中直接抽取样本非常困难，MCMC抽样是通过建议分布来获得后验概率密度函数的带权子样，综合到误差分布通常满足正态分布和建议分布为对称分布，因此建议分布采取正态分布。

在步骤S5，由于污染物溯源统计MCMC抽样是根据后验概率分布情况抽样的，则可能会出现在一个局部空间最优的循环抽样(例如整个后验概率分布有两个以上的波峰，MCMC抽样可能会在其中不是最优的一个波峰循环抽样，容易造成局部解)，设置随机数，使得能够出现跳出局部空间的抽样，从均匀分布闭区间上U[0，1]中抽取一个随机数u，判断该候选值是否接受，以扩大随机搜索的范围避免搜索限于局部解。从当前位置X_t获得候选值X^*，马尔科夫链从X_t位置移动到X^*的接受概率：

如果

候选值被接受，则在第t+1抽样的样本值为X^*，即X_t+1＝X^*，X_t＝X_t+1，t＝t+1，p(X_t|y)＝P(X^*|y)，输出当前参数X_t及对应的后验概率分布p(X_t|y)；否则在t+1次抽样的样本链停留在原来的状态X_t，即X_t+1＝X_t，X_t＝X_t+1，t＝t+1。

在步骤S6，根据实际情况，设置迭代最终次数，重复步骤S4，步骤S5直至迭代完成。

在步骤S7，流程结束。

步骤4，在Matlab中输出后验概率直方图，进行贝叶斯统计溯源结果的可视化分析，推求未知的偷排节点、排放量和排放时间的取值的概率分布。

实施例

图2为区域污水管网平面布置图，管网中水流流向如箭头所示，该污水管网共有管道64段，65个节点，管径为400～800mm，其中PFK1为该污水管网的总排放口，在管网下游主干管上节点J1处布设水质监测点。该污水管网流态恒定，污染物质在污水管道的迁移传输服从一维水质模型，降解过程为一级衰减，衰减系数为0.25。假定某日1:00时，在节点28处的某工厂瞬时排入重量为1000kg的污染物BX，在下游监测点J1处观测到BX浓度变化曲线。现利用监测节点J1被扰动后的带有噪声的浓度场分布反演上游BX排放源的排放强度、排放位置和排放时间。

本发明所要解决的问题就是在未知超标工业废水偷排污染源任何信息的情况下，通过在下游监测点J1连续监测水样中污染物浓度，反推出所有可能偷排的节点、排放量及排放时间。在监测点1处首次观测时间为1:10，后每隔10min观测一次，共6个数据，观测浓度序列见图3。

基于经验和实际情况，确定先验信息取值的上下限，待反演未知参数先验分布为该区间内的均匀分布，对应的先验概率密度函数分别为：

按照本发明所述步骤在SWMM5.0构建排水模型并在Matlab软件中编写贝叶斯MCMC抽样程序，迭代次数为4000。在Matlab里面，利用unifrnd等统计函数实现了抽样中初始值的随机抽取，利用normrnd、linspace等统计函数生成建议分布随机数，并且计算出后验概率密度函数，并且输出当前参数及对应的后验概率密度，最后利用hist等统计函数生成后验概率密度直方图，得到概率密度直方图，图4、图5和图6，从三个参数的后验概率直方图可以看出污染物排放量、排放时间、排放节点的后验概率密度较集中，可以根据后验概率直方图中的未知参数概率大小情况推断污染事件的未知参数排放量M∈[900000g，1260000g]，排放时间在T∈[35min，72min]之间，排放节点在Jx∈[20，45]之间。

表1所示，通过将污染事件的反演参数排放量M(g)、排放时间T(min)的抽样结果与真实值对比，可知污染物排放时间和排放量中值误差较小。从本实例可以看出，本发明方法是有效可靠的。

表1：MCMC抽样结果和实际值对比表。

本发明通过贝叶斯的反演算法，三个未知污染物排放参数的取值范围大大缩小，可能的偷排节点数由65个减少至26个，节点数减少了近70％，数学模型溯源的结果极大的缩小了未知参数的取值范围，将减少现场人工物理溯源的搜索范围及工作量，提高溯源效率，且抽样高效合理，抽样次数大大减少，克服以往数学模型溯源工作量大且耗时的缺点。

Claims (4)

1.一种污水管网超标排放工业废水的贝叶斯统计溯源方法，其特征是，包括以下步骤：

步骤1：在SWMM软件中构建污水管网模型，生成初始*.inp文件，然后在Matlab中打开初始*.inp文件，对*.inp文件进行读写，对所有可能偷排的节点、排放量和排放时间三个参数进行编辑，外部调用swmm5.exe程序，直接利用更新后的参数输入运行SWMM，将包括监测点的污染物浓度时间序列运行结果储存在*.rpt文件中；

步骤2：在Matlab中实现对*.rpt文件的读写，对排放节点、排放量和排放时间三个未知参数在先验信息范围内进行污染物溯源统计MCMC抽样，每次抽样后，同步更新*.inp文件中相应的参数变量，外部运行SWMM，将更新的节点污染物浓度时间序列储存在*.rpt文件中；

步骤3：在Matlab中提取*.rpt文件中相应的模拟结果的时间序列，结合实际监测节点处的污染物质浓度时间序列，利用贝叶斯理论反演计算，考虑监测过程及模拟过程的误差，计算似然函数值和待反推的未知参数的后验概率分布；

步骤4，在Matlab中输出后验概率直方图，进行贝叶斯统计溯源结果的可视化分析，推求未知的偷排节点、排放量和排放时间的取值的概率分布。

2.根据权利要求1所述的一种污水管网超标排放工业废水的贝叶斯统计溯源方法，其特征是：

在步骤1中，排放量和排放时间的取值服从取值空间内的连续型均匀分布，它们的先验概率密度函数为：

排放节点X₃服从其取值范围的离散均匀分布，其概率函数为：

在某节点为排放节点情况下的先验总概率密度函数p(X)

式中p(X_f)代表取值区间[a_f,b_f]上的均匀分布；f表示第f个参数；H代表需要反推的参数的个数；X_f代表第f种未知参数；п为连乘符号；K代表第K个排放节点；m代表共有m个排放节点；a、b分别为参数取值的上限和下限。

3.根据权利要求2所述的一种污水管网中超标排放工业废水贝叶斯溯源方法，其特征是，在步骤2中，所述污染物溯源统计MCMC抽样的步骤包括：

步骤21在未知参数的先验信息范围内随机产生初始点X₀；

步骤22令t＝0，当前参数X_t＝X₀，利用SWMM一维污染物质对流传输模型计算出对应的污染物浓度值，并与实际监测数据比较计算吻合程度，通过贝叶斯理论计算出对应的后验概率密度p(X_t|y)，输出X_t及p(X_t|y)；

步骤23假定新的参数候选值X^*服从以当前参数X_t为均值，方差σ的正态分布为建议分布，根据建议分布在当前参数X_t状态下抽取X^*，利用水质模型计算出X^*对应的污染物浓度值，与实际监测数据比较似然程度，计算出未知参数的后验概率密度P(X^*|y)；

步骤24从均匀分布闭区间上U[0，1]中抽取一个随机数u，判断该候选值是否接受，以扩大随机搜索的范围避免搜索限于局部解，从当前位置X_t获得候选值X^*，马尔科夫链从X_t位置移动到X^*的接受概率：

如果

候选值被接受，则在第t+1抽样的样本值为X^*，即X_t+1＝X^*，X_t＝X_t+1，t＝t+1，p(X_t|y)＝P(X^*|y)，输出当前参数X_t及对应的后验概率分布p(X_t|y)，否则在t+1次抽样的样本链停留在原来的状态X_t，即X_t+1＝X_t，X_t＝X_t+1，t＝t+1；

步骤25，设置迭代次数，重复步骤23，步骤24直至迭代完成。

4.根据权利要求3所述的一种污水管网超标排放工业废水的贝叶斯统计溯源方法，其特征是：

在步骤3中，在利用贝叶斯理论反演计算中，构造表征污染物扩散模型模拟的污染物浓度理论值和实测数据拟合程度的似然函数p(y|X)：

式中：n为测量数据个数；排放强度为M；排放节点为Jx；排放时间为t；已知监测点观测到污染物浓度观测值Y＝{Y₁，Y₂，…，Y_i，…，Y_n}；污染物扩散模型模拟的污染物浓度理论值N＝{N₁，N₂，…，N_i，…，N_n}；污染事件发生后在监测点第i个节点观测到的污染物浓度值为Y_i；在未知参数X为M，Jx，T的条件下，通过污染物扩散模型模拟的监测点第i个节点处污染物浓度理论值为N_i(M，Jx，T|X)；δ是似然函数标准差；

根据贝叶斯理论，计算出未知参数X的后验概率密度：

式中，X为未知参数；p(X|y)为未知参数X的后验概率密度函数；α为比例常数。