CN105930670A - 基于模型参数不确定性的河流突发污染事故动态预测方法 - Google Patents

Abstract

本发明公开了基于模型参数不确定性的河流突发污染事故动态预测方法。S1.利用参数率定、经验数据以及查阅文献的方法生成不确定性参数；S2.利用S1中不确定性参数值和初始断面污染物浓度计算事故点以后各断面各时刻的污染物浓度值；S3.选取似然函数计算S2中不同参数组所对应的似然值；S4.利用S3中的似然值可以估算出在一定置信度水平下模型预测结果的不确定范围，扩展到事故发生以后的整个事件维度可以获得模型预测结果的不确定性区间；S5.利用实测数据对模拟预测结果不断进行更新和校正，重复上述步骤，可求得最新更新的预测结果。该方法结合不确定性方法、动态更新理论、普适似然不确定性算法，实现了基于模型参数不确定性的河流突发污染事故水质预测。

Description

基于模型参数不确定性的河流突发污染事故动态预测方法

技术领域

本发明涉及污染物模拟领域，具体提出一种基于普适似然不确定性算法和动态更新的突发性水污染事故预警方法。

背景技术

河流、湖泊以及水库是人类重要的淡水生态资源，沿江临湖地区也通常是人类活动频繁和各种生物生息繁衍的重要区域。然而，随着经济的高速发展，人类的生产和生活活动对水资源安全造成了严重影响，这引发了众多经济问题和社会问题。目前，尚无有效手段彻底阻止人类活动对于水环境的恶劣影响，在今后相当长的时间范围内，水环境问题仍将是我们不得不面对的突出问题。

在众多水环境问题中，河流突发污染事故是一类发生频率高，危害性严重的事故。河流突发污染事故通常具有以下特征：突发、不易预知，不确定性高以及危害巨大等，这使得突发污染事故发生后的水质预测工作和应急处理工作充满困难。为使应急决策者更好地应对突发污染事故，对事故发生后污染物浓度时空变化规律的描述显得尤为重要。

污染事故的发生、发展以及演变具有很大的不确定性，主要表现为：(1)发生时间、空间的不确定性；(2)污染物的不确定性；(3)事故流域特性的不确定性；(4)污染方式的不确定性；(5)事故记录数据的不确定性等。这些不确定性的客观存在为污染物浓度的预测工作带来了困难和挑战。

确定性水质模型虽然理论系统完善、甚至可以准确细致地表现污染物迁移扩散过程，但由于河流环境本身的复杂性使得机理建模很难准确，而初始条件和水文数据又很难全面获知，再加上模型参数的率定工作很难进行，尤其是突发污染事故要求模型尽快给出预测数据，使得短期建立起优秀的确定性水质模型变得非常困难。在这样的背景下，随着不确定性理论研究的逐步成熟，不确定性水质模型逐渐成为了相关工作者的研究热点。

不确定性水质模型主要包括以统计分析为主要研究方法的随机理论模型，以时间序列分析、人工神经网络、遗传算法等为代表的主要基于数据处理的数学模型，以及以灰色系统、模糊数学为处理手段的灰色预测模型和模糊模型。

目前关于参数或者变量方面的不确定性分析，用的最多的是蒙特卡罗方法，该方法基于随机抽样进行模拟，其结果与参数分布情况的关系较大，该方法适用于各类复杂的非线性系统。

发明内容

为了解决现有技术的不足，本发明提供了一种基于模型参数不确定性的河流突发污染事故动态预测方法。

一种基于模型参数不确定性的河流突发污染事故动态预测方法，包括以下步骤：

S1.不确定性参数组的生成；

S2.模拟预测：利用S1中不确定性参数值和初始断面污染物浓度计算事故点以后各断面各时刻的污染物浓度值；

S3.似然值计算：选取似然函数计算S2中污染物浓度值所对应的似然值；

S4.概率密度函数计算：利用S3中的似然值估算出在模型预测结果的不确定范围，扩展到事故发生以后的整个事件维度获得模型预测结果的不确定性区间；

S5.利用实测浓度数据比较模拟预测结果后作出判断，如果实测浓度数据与模拟预测结果之间具有误差，则进入S2，重复上述步骤，得到更新的预测结果；

如果实测浓度数据与模拟预测结果相符，则结束。

所述的S1步骤中，根据模型的不确定性参数，并确定其初始分布；对于选定的河流突发污染事故模型，模型的不确定性参数为E，u，k，其中E为纵向弥散系数，u为纵向水流流速，k为综合衰减速率系数，E，u，k的初始分布取均匀分布，根据历史数据和相关经验分别确定E，u，k的参数取值范围，然后按照蒙特卡罗方法随机获取E，u，k的组合θ，θ表示如下：

θ＝[Rand_E,Rand_u,Rand_k] (1)

式中，Rand代表对模型参数按照其分布进行随机取值。

所述的方法，将初始断面各时刻的污染物浓度值和参数组θ代入公式(2)，迭代算得事故点以后各断面各时刻的污染物浓度值，代入不同参数组θ_i(i＝1，2，…，N)，N为蒙特卡罗模拟次数，分别获得与其对应的各断面各时刻的污染物浓度值序列，再利用蒙特卡罗方法进行统计得到河流中任意位置、任意时刻的污染物浓度值分布，即概率密度函数；

${C_i}^{j+1}={C_i}^j+(E_i\frac{{C_{i+1}}^j+{C_{i-1}}^j}{{\mathrm{\Δx}}^2}+u_i\frac{{C_{i-1}}^j}{\mathrm{\Δ}x})\mathrm{\Δ}t-(E_i\frac{2{C_i}^j}{{\mathrm{\Δx}}^2}+u_i\frac{{C_i}^j}{\mathrm{\Δ}x}+k_i{C_i}^j)\mathrm{\Δ}t---\left(2\right),$

式中：C为平均污染物浓度，mg/L；u为纵向水流流速，单位为km²/h；E为纵向弥散系数，单位为km²/h；x为河水流动距离，单位为km；k为综合衰减速率系数，单位为h^-1；t为时间，单位为h；i表示第i段河段，j来表示第j个时间间隔；Δx和Δt分别表示空间和时间间隔大小。

所述的方法，过程如下：

设有统计独立的随机变量X_i(i＝1，2，…，k)，其对应的概率密度函数分别为f_x1，f_x2，…，f_xk，功能函数式为Y＝f(x₁，x₂，…，x_k)；

首先根据各随机变量X_i的分布，产生N组随机数x₁，x₂，…，x_k值，然后计算功能函数值Y_j＝f(x₁，x₂，…，x_k)(j＝1，2，…，N)，接下来根据计算结果，等间隔地选取m组间隔为s的Y_j的取值区间，m组区间互相连接、不重叠，并且包含全部Y_j取值，再分别统计落入每个取值区间的Y_j个数，即Quantity(Y_j)，并按式(3)统计变量密度，

$Density\left(Y_j\right)=\frac{Quantity\left(Y_j\right)}{s\·N}---\left(3\right)$

当N足够大，s足够小，m足够大时，根据伯努利大数定理和正态随机变量的特性以及频率分布直方图与概率密度函数之间的关系有：全部Density(Y_j)构成了Y_j的概率密度函数，而该概率密度函数即为被估计变量Y的实际取值分布。

似然函数表征某个参数组所预测出来的浓度值与实测浓度值之间的差异程度，采用Nash‐Sutcliffe确定性系数：

$L=R^2=1-\frac{{\mathrm{\&sigma;}}_e^2}{{\mathrm{\&sigma;}}_o^2};{\mathrm{\&sigma;}}_e^2<{\mathrm{\&sigma;}}_o^2---\left(4\right)$

式中：L代表似然值；R²代表确定性系数；代表预测序列的误差方差；代表实测序列的方差；计算似然值之前需要首先获取断面的实测污染物浓度值序列，结合预测值，代入公式(4)即可得到不同参数组θ_i所对应的似然值。

按照似然值的大小对不同参数组θ_i的预测结果进行统计加权，即按照式3统计时，Quantity(Yj)不再代表落入某个取值区间的结果数目，而是代表计算结果落入某个取值区间的似然值之和，相应的N代表全部似然值之和，进而获得各断面各时刻的预测结果分布函数；在此基础上，估算出在置信度水平下模型预测结果的不确定性范围，扩展到事故发生以后的整个时间维度获得模型预测结果的不确定性区间。

所述的S5是利用实测数据对模拟预测结果不断进行更新、校正的过程，其具体计算过程是：新的数据序列获取以后，先按照步骤5计算参数组θ_i的观测似然值，然后使用贝叶斯公式对前一次计算得到的似然值结果，即先验似然值，进行更新，得到后验似然值结果，再按照步骤6进行更新后的模拟预测，贝叶斯公式表示如下：

$L\left(Y\right|{\mathrm{\&theta;}}_i)=\frac{L\left({\mathrm{\&theta;}}_i\right|Y)\&CenterDot;L\left({\mathrm{\&theta;}}_i\right)}{C}---\left(5\right)$

式中：L(Y|θ_i)为后验似然值；L(θ_i|Y)为观测似然值；L(θ_i)为先验似然值；C为归一化因子；Y为预测因变量(即预测结果)；θ_i为第i组参数组；

当再有新的实测数据获取时，前一次更新过的后验似然值结果即变为当前计算的先验似然值结果，重复上述步骤，即可求得更新的预测结果。

本发明的有益效果：

给出了河流突发污染事故不确定性的河流突发污染事故水质预测一般框架，结合蒙特卡罗方法和GLUE算法实现前述框架，然后以此为基础，引入动态更新方法对次水质预测框架进行改进，并可利用事故发生后的污染物浓度实时数据对模型预测结果进行动态校正，使得模型预测结果随着污染事故的持续进行仍然保持良好的准确性和可靠性。

附图说明

图1为基于模型参数不确定性的河流突发污染事故动态预测方法的流程图；

图2为模型预测结果的动态更新流程图；

图3为河流Q以及事故发生点、监测点位置示意图；

图4为模拟效率X随着似然值临界值m选取的变化曲线；

图5为监测点A的实测值与不确定性预测区间；

图6为监测点B的实测值与不确定性预测区间；

图7为监测点C的实测值与不确定性预测区间；

图8为监测点B的污染物浓度实测值与不确定性预测区间对比。

具体实施方案

下面结合附图对本发明进一步描述。

基于模型参数不确定性的河流突发污染事故动态预测的基本步骤如图1所示，主要包含：参数组生成，模拟预测，似然函数选取与似然值计算以及概率密度函数计算。

模拟预测虽然实现了不确定性预测，但是尚未利用后续测得的污染物浓度值数据对预测结果进行动态校正，其预测准确性随着时间的推移而无法保证。所以，只有不断利用后续测得的浓度值来实时更新、校正预测结果才能使模型预测结果一直保持较高的准确性。因此，在图1的基础上，对预测结果动态更新，其过程如图2所示：当没有新的数据对似然函数结果进行动态更新时，直接利用这个似然函数值进行模拟预测，即图1的一次计算；当需要新数据对似然函数结果进行动态更新时，以这个函数值作为先验似然值结合观测似然值利用贝叶斯公式计算后验似然值从而达到动态更新的目的。

一种基于模型参数不确定性的河流突发污染事故动态预测方法，包括以下步骤：

S1.不确定性参数组的生成；

S2.利用S1中不确定性参数值和初始断面污染物浓度计算事故点以后各断面各时刻的污染物浓度值；

S3.选取似然函数计算S2中不同参数组所对应的似然值；

S4.利用S3中的似然值可以估算出在一定置信度水平下模型预测结果的不确定范围，扩展到事故发生以后的整个事件维度可以获得模型预测结果的不确定性区间；

S5.利用实测数据对模拟预测结果不断进行更新和校正，重复上述步骤，可求得最新更新的预测结果。

所述的S1步骤中，参数组生成部分需要首先明确模型有哪些不确定性参数，并确定其初始分布。对于选定的河流突发污染事故模型，模型的不确定性参数为E，u，k。通常由于缺少足够的数据支持，E，u，k的初始分布可取均匀分布。所以实际的计算过程是根据历史数据和相关经验分别确定E，u，k的参数取值范围，然后按照蒙特卡罗方法随机获取E，u，k的组合θ，θ可表示如下：

θ＝[Rand_E,Rand_u,Rand_k] (1)

式中，Rand代表对模型参数按照其分布进行随机取值。

将初始断面各时刻的污染物浓度值和参数组θ代入公式(2)，可以迭代算得事故点以后各断面各时刻的污染物浓度值。代入不同参数组θ_i(i＝1，2，…，N)(N为蒙特卡罗模拟次数)可以分别获得与其对应的各断面各时刻的污染物浓度值序列，再利用蒙特卡罗方法进行统计即可得到河流中任意位置、任意时刻的污染物浓度值分布(概率密度函数)。通常突发污染事故发生后，对于整条河流而言只有部分断面设置了采样点，即只有部分河流位置有实测数据，所以模拟预测部分更加关注这些有采样点的河道位置的污染物浓度值分布。

${C_i}^{j+1}={C_i}^j+(E_i\frac{{C_{i+1}}^j+{C_{i-1}}^j}{{\mathrm{\&Delta;x}}^2}+u_i\frac{{C_{i-1}}^j}{\mathrm{\&Delta;}x})\mathrm{\&Delta;}t-(E_i\frac{2{C_i}^j}{{\mathrm{\&Delta;x}}^2}+u_i\frac{{C_i}^j}{\mathrm{\&Delta;}x}+k_i{C_i}^j)\mathrm{\&Delta;}t---\left(2\right)$

所使用的蒙特卡罗方法其基本原理如下：

由概率定义知，某事件发生的概率可以用大量试验中该事件发生的频率来估算，当试验足够多时，可以认为该事件发生的频率即为其发生的概率。因此，可以先对影响某一事件是否发生的随机变量进行大量随机抽样，然后把这些抽样值分别代入功能函数式，并统计计算结果，进而求得某一事件的发生概率。蒙特卡罗方法正是基于此思路进行分析的，其计算过程如下：

设有统计独立的随机变量X_i(i＝1，2，…，k)，其对应的概率密度函数分别为f_x1，f_x2，…，f_xk，功能函数式为Y＝f(x₁，x₂，…，x_k)。

首先根据各随机变量X_i的分布，产生N组随机数x₁，x₂，…，x_k值，然后计算功能函数值Y_j＝f(x₁，x₂，…，x_k)(j＝1，2，…，N)，接下来根据计算结果，等间隔地选取m组间隔为s的Y_j的取值区间(m组区间互相连接、不重叠，并且包含全部Y_j取值)，再分别统计落入每个取值区间的Y_j个数(Quantity(Y_j))，并按式(3)统计变量密度，

$Density\left(Y_j\right)=\frac{Quantity\left(Y_j\right)}{s\&CenterDot;N}---\left(3\right)$

当N足够大、s足够小、m足够大时，根据伯努利大数定理和正态随机变量的特性以及频率分布直方图与概率密度函数之间的关系可有：全部Density(Y_j)构成了Y_j的概率密度函数，而该概率密度函数即为被估计变量Y的实际取值分布。

蒙特卡罗方法主要用于解决两类问题，即纯数学求解问题和随机性问题，其对于解决具有随机性或者不确定性的实际问题有很强的适用能力。本申请所研究的污染物浓度模拟预测问题属于随机性问题，通过大量的随机试验从模型参数的先验分布中抽取参数代入河流突发污染事故预测模型进行运算，然后统计得到污染物浓度在各时空维度下的概率密度函数，这是对结果进行进一步不确定性分析的基础。

似然函数表征某个参数组所预测出来的浓度值与实测浓度值之间的差异程度，有多种可供选择，而如何选取带有一定的主观性，通常可采用Nash‐Sutcliffe确定性系数：

$L=R^2=1-\frac{{\mathrm{\&sigma;}}_e^2}{{\mathrm{\&sigma;}}_o^2};{\mathrm{\&sigma;}}_e^2<{\mathrm{\&sigma;}}_o^2---\left(4\right)$

式中：L代表似然值；R²代表确定性系数；代表预测序列的误差方差；代表实测序列的方差。

计算似然值之前需要首先获取某断面的实测污染物浓度值序列，结合前面计算得到的预测值，代入公式(4)即可算得不同参数组θ_i所对应的似然值，然后按照一定标准(似然值大于某一阈值，阈值的选取带有一定主观性)剔除明显无法模拟污染物扩散情况的参数组，以备后续使用。

按照似然值的大小对不同参数组θ_i的预测结果进行统计加权(即按照式3统计时，Quantity(Yj)不再代表落入某个取值区间的结果数目，而是代表计算结果落入某个取值区间的似然值之和，相应的N代表全部似然值之和)，进而获得各断面各时刻的预测结果分布函数。在此基础上，可以估算出在一定置信度(如95％)水平下模型预测结果的不确定性范围，扩展到事故发生以后的整个时间维度可以获得模型预测结果的不确定性区间。

动态更新过程是利用实测数据对模拟预测结果不断进行更新、校正的过程，其具体计算过程是：新的数据序列获取以后，先按照步骤5计算参数组θ_i的观测似然值，然后使用贝叶斯公式对前一次计算得到的似然值结果(先验似然值)进行更新，得到后验似然值结果，再按照步骤6进行更新后的模拟预测。贝叶斯公式在本应用中可表示如下：

$L\left(Y\right|{\mathrm{\&theta;}}_i)=\frac{L\left({\mathrm{\&theta;}}_i\right|Y)\&CenterDot;L\left({\mathrm{\&theta;}}_i\right)}{C}---\left(5\right)$

式中：L(Y|θ_i)为后验似然值；L(θ_i|Y)为观测似然值；L(θ_i)为先验似然值；C为归一化因子；Y为预测因变量(即预测结果)；θ_i为第i组参数组。

当再有新的实测数据获取时，前一次更新过的后验似然值结果即变为当前计算的先验似然值结果，重复上述步骤，即可求得最新更新的预测结果。

实施例

以河流Q上游某处发生的化学品泄漏污染事故为实例进行模拟计算与结果分析，各项数据源于文献。事故发生点以后的河段全长约为160km，将河道均分为320段(每段500米，即Δx＝500m)进行模拟计算(Δt＝60s)。按照计算流程，首先需要明确不确定性分析中的参数分布和范围。根据大量历史数据和现场测定结果可以确定事发河流流速u的取值范围为：0.1‐1.0(m/s)；通过参考相关文献，确定河流纵向弥散系数E为：100‐300(m²/s)；综合降解系数k为：0.58×10^‐6‐1.74×10^‐6(s^‐1)。模型参数E，u，k的初始分布取为均匀分布。选取污染事故发生后，监测数据记录相对完整的3个河流断面(分别记为监测点A，B，C)作为研究对象，河道与采样断面示意图如附图3所示。

1.有效参数组合的临界值选取

计算完所有参数组所对应的似然值后需要根据似然值的大小对参数组进行取舍，取舍标准带有一定的主观性。可以采用模拟效率(用符号X表示)对取舍标准进行分析，模拟效率表示被选择用于进行不确定性预测的参数组占全部由蒙特卡罗方法生成的参数组的百分比，其计算公式如下：

$X=\frac{N_{R^2>m}}{N_{all}}---\left(6\right)$

式中，X表示模拟效率；m表示似然值选取临界值；分子表示似然值大于临界值m的参数组数目；分母表示GLUE算法中由蒙特卡罗方法所生成的全部参数组数目。

以监测点B的实测数据计算不同似然值临界值m选取下的模拟效率X，可以得到如附图4所示的曲线。

从图中可见，模拟效率X随着临界值m的增大而逐渐下降。若m值选取过大，则X偏小，即可接受的参数组数目过少，此时可理解为寻找“最优解”，由前文分析可知这种情况会忽视各种不确定性的影响，其预测结果反而不准确。若m值选取过小，则X偏大，此时过多考虑了各种不确定性因素以至于让一些无法良好模拟当前河流污染物扩散状况的参数组也对最终预测结果有所贡献，很明显这种情况也是不合理的。当m取值0.5时，模拟效率近似为60％，此时既让模型模拟效率保证在50％以上，又使计算结果包含了曲线切线附近区域所对应的参数组(该区域即使m值变化较小，X值也会发生较大变化，若不包含该区域则会丢掉很多有价值的参数组)，所以选定似然值临界值为0.5。

2.在一定置信度下的预测与实际结果对比

进行5000次蒙特卡罗模拟，分别对监测点A，B，C取事故发生后15h，40h，50h以内的实测数据对预测结果动态更新，再分别统计监测点A，B，C的预测结果，可以获得监测点A，B，C在各时间维度下的模型预测浓度值分布(概率密度函数)。对于得到的概率密度函数，可以计算出在95％置信度下的最优置信区间，即预测结果的不确定性范围，其与实测值的对比图分别如附图5、6、7所示。

为评价不确定性预测区间的优劣，引入覆盖率(CR)的概念。覆盖率是指预测区间覆盖实测数据的比率，它是最常用的预测区间评价指标。CR越大，则表示预测区间越能包含全部实测数据。

从图中可以看出，大部分的实测值(CR_A＝97％；CR_B＝92％；CR_C＝90％)都落在置信度为95％的不确定性预测区间以内，但也有少部分的实测值落于不确定性预测区间外部。实测值落于预测区间外部的原因有很多：采样测量误差，历史数据信息偏差，模型本身准确性偏差等。

而若把置信度设定为90％，则有相对较少的实测值落于不确定性预测区间内，但由于不确定性预测区间变小，其预测不确定性反而降低，有助于模型结果的预测工作。总结可有，置信度取值越低，则不确定性预测范围越小，预测区间的可靠性越低，但有利于进行相对准确的结果预测；反之，置信度取值越高，则不确定性预测范围越大，预测区间的可靠性越高，但较难对结果进行准确预测。所以，在具体选定预测置信度时，需要综合考虑各种方面的因素。

3.实际应用过程

3中计算得到的不确定性预测结果是以某监测点自身测得的浓度数据对该监测点的不确定性预测结果进行动态更新，这样可以得到相对准确的预测结果(以监测点B的不确定性预测结果为例，其覆盖率高达92％)。然而，河流突发污染事故发生后，倘若等到监测点有可用实测数据对不确定性预测结果进行更新时往往时间过久(尤其对于河流中、下游的监测点而言时间更久，以监测点B为例需要等到污染事故发生后约30h才有可用实测数据获取)，而应急决策者却希望更早地得知污染物在整条河流中的扩散规律，以便从容准备应急方案，避免污染事故造成更大的经济损失和社会负面影响。

为解决上述矛盾，在实际应用时，可采取如下改进做法：

(1)先采用靠近源头的监测点实测数据计算模型参数组的似然值分布，并粗略认为该似然值分布可以描述河流其他断面的污染物扩散规律，然后即可据此分布分别计算河流各关注断面(通常被关注的河流断面都会设置监测点)的不确定性预测区间。

(2)当有更靠近某关注断面的监测点获取到实测数据时，则以该监测点实测数据对应的参数组似然值分布替换之前的似然值分布，并计算不确定性预测区间。

(3)最后，当该关注断面自身获取到实测数据后，则可如3一样计算较为准确的不确定性预测结果。

由于河流环境本身所具有的不确定性，不同河流断面所对应的模型参数组似然值分布本应不同，但是为了更早地进行不确定性预测，可以先忽略这些差异的影响而进行较为粗略的不确定性预测，以便让应急决策者尽早得到数据支持，并进行前期应急决策；当某关注断面自身获取到实测数据后，即可按照本专利方法进行更为准确地不确定性预测，此时的预测结果可以帮助决策者校正之前的相关决策，进而更有针对性地应对污染事故。

以监测点B为例，倘若要进行准确地不确定性预测，需要等到事故发生后40h才有足够的数据支持，所需时间明显过久。实际应用时，可以先用监测点A获取的实测浓度数据计算模型参数组的似然值分布，然后以该分布计算监测点B的不确定性预测区间，其不确定性预测结果如附图8

按照这种办法，在事故发生后15h即可对监测点B的污染物浓度变化规律进行不确定性预测，虽然预测结果准确性较差(覆盖率不到43％)，但却足以为应急决策者在对污染物在河流中的扩散规律丝毫不清楚时提供数据参考，便于其准备充足的方案应对污染事故。事故发生后30h，监测点B自身获取到实测数据后，即可据此计算得到3所述的更为准确的不确定性预测结果，决策者可以据此采取更进一步的污染事故应对策略。

以上阐述的是本发明给出的实例，仿真结果体现了本发明所提出的技术方案对于A河流突发污染事故的动态预测结果。需要指出的是，本发明不只限于上述实施例，对于其他突发污染事故，采用本发明的技术方案也能给出很好的预警效果。

Claims (7)

1.一种基于模型参数不确定性的河流突发污染事故动态预测方法，其特征在于，包括以下步骤：S1.不确定性参数组的生成；

S2.模拟预测：利用S1中不确定性参数值和初始断面污染物浓度计算事故点以后各断面各时刻的污染物浓度值；

S3.似然值计算：选取似然函数计算S2中污染物浓度值所对应的似然值；

S4.概率密度函数计算：利用S3中的似然值估算出在模型预测结果的不确定范围，扩展到事故发生以后的整个事件维度获得模型预测结果的不确定性区间；

S5.利用实测浓度数据比较模拟预测结果后作出判断，如果实测浓度数据与模拟预测结果之间具有误差，则进入S2，重复上述步骤，得到更新的预测结果；如果实测浓度数据与模拟预测结果相符，则结束。

2.根据权利要求1所述的方法，其特征在于，所述的S1步骤中，根据模型的不确定性参数，并确定其初始分布；对于选定的河流突发污染事故模型，模型的不确定性参数为E，u，k，其中E为纵向弥散系数，u为纵向水流流速，k为综合衰减速率系数，E，u，k的初始分布取均匀分布，根据历史数据和相关经验分别确定E，u，k的参数取值范围，然后按照蒙特卡罗方法随机获取E，u，k的组合θ，θ表示如下：

θ＝[Rand_E,Rand_u,Rand_k] (1)

式中，Rand代表对模型参数按照其分布进行随机取值。

3.根据权利要求2所述的方法，其特征在于，将初始断面各时刻的污染物浓度值和参数组θ代入公式(2)，迭代算得事故点以后各断面各时刻的污染物浓度值，代入不同参数组θ_i(i＝1，2，…，N)，N为蒙特卡罗模拟次数，分别获得与其对应的各断面各时刻的污染物浓度值序列，再利用蒙特卡罗方法进行统计得到河流中任意位置、任意时刻的污染物浓度值分布，即概率密度函数；

${C_i}^{j+1}={C_i}^j+\left(E_i\frac{{C_{i+1}}^j+{C_{i-1}}^j}{{\mathrm{\&Delta;x}}^2}+u_i\frac{{C_{i-1}}^j}{\mathrm{\&Delta;}x}\right)\mathrm{\&Delta;}t-\left(E_i\frac{2{C_i}^j}{{\mathrm{\&Delta;x}}^2}+u_i\frac{{C_i}^j}{\mathrm{\&Delta;}x}+k_i{C_i}^j\right)\mathrm{\&Delta;}t---\left(2\right),$

式中：C为平均污染物浓度，mg/L；u为纵向水流流速，单位为km²/h；E为纵向弥散系数，单位为km²/h；x为河水流动距离，单位为km；k为综合衰减速率系数，单位为h^-1；t为时间，单位为h；i表示第i段河段，j来表示第j个时间间隔；Δx和Δt分别表示空间和时间间隔大小。

4.根据权利要求3所述的方法，其特征在于，过程如下：设有统计独立的随机变量X_i(i＝1，2，…，k)，其对应的概率密度函数分别为f_x1，f_x2，…，f_xk，功能函数式为Y＝f(x₁，x₂，…，x_k)；首先根据各随机变量X_i的分布，产生N组随机数x₁，x₂，…，x_k值，然后计算功能函数值Y_j＝f(x₁，x₂，…，x_k)(j＝1，2，…，N)，接下来根据计算结果，等间隔地选取m组间隔为s的Y_j的取值区间，m组区间互相连接、不重叠，并且包含全部Y_j取值，再分别统计落入每个取值区间的Y_j个数，即Quantity(Y_j)，并按式(3)统计变量密度，

$Density\left(Y_j\right)=\frac{Quantity\left(Y_j\right)}{s\&CenterDot;N}---\left(3\right)$

当N足够大，s足够小，m足够大时，根据伯努利大数定理和正态随机变量的特性以及频率分布直方图与概率密度函数之间的关系有：全部Density(Y_j)构成了Y_j的概率密度函数，而该概率密度函数即为被估计变量Y的实际取值分布。

5.根据权利要求4所述的方法，其特征在于，似然函数表征某个参数组所预测出来的浓度值与实测浓度值之间的差异程度，采用Nash‐Sutcliffe确定性系数：

$L=R^2=1-\frac{{\mathrm{\&sigma;}}_e^2}{{\mathrm{\&sigma;}}_o^2};{\mathrm{\&sigma;}}_e^2<{\mathrm{\&sigma;}}_o^2---\left(4\right)$

式中：L代表似然值；R²代表确定性系数；代表预测序列的误差方差；代表实测序列的方差；计算似然值之前需要首先获取断面的实测污染物浓度值序列，结合预测值，代入公式(4)即可得到不同参数组θ_i所对应的似然值。

6.根据权利要求5所述的方法，其特征在于，按照似然值的大小对不同参数组θ_i的预测结果进行统计加权，即按照式3统计时，Quantity(Yj)不再代表落入某个取值区间的结果数目，而是代表计算结果落入某个取值区间的似然值之和，相应的N代表全部似然值之和，进而获得各断面各时刻的预测结果分布函数；在此基础上，估算出在置信度水平下模型预测结果的不确定性范围，扩展到事故发生以后的整个时间维度获得模型预测结果的不确定性区间。

7.根据权利要求6所述的方法，其特征在于，所述的S5是利用实测数据对模拟预测结果不断进行更新、校正的过程，其具体计算过程是：新的数据序列获取以后，先按照步骤5计算参数组θ_i的观测似然值，然后使用贝叶斯公式对前一次计算得到的似然值结果，即先验似然值，进行更新，得到后验似然值结果，再按照步骤6进行更新后的模拟预测，贝叶斯公式表示如下：

$L\left(Y\right|{\mathrm{\&theta;}}_i)=\frac{L\left({\mathrm{\&theta;}}_i\right|Y)\&CenterDot;L\left({\mathrm{\&theta;}}_i\right)}{C}---\left(5\right)$

式中：L(Y|θ_i)为后验似然值；L(θ_i|Y)为观测似然值；L(θ_i)为先验似然值；C为归一化因子；Y为预测因变量(即预测结果)；θ_i为第i组参数组；

当再有新的实测数据获取时，前一次更新过的后验似然值结果即变为当前计算的先验似然值结果，重复上述步骤，即可求得更新的预测结果。