基于蒙特卡洛仿真的河湖水质预测及富营养化风险评估方法

技术领域

本发明设计一种基于蒙特卡洛仿真的河湖水质预测方法及富营养化风险评估方法，属于环境工程技术领域。

背景技术

随着现代社会的飞速发展，水体富营养化已成为一个全球性的重大水环境问题。目前，在世界范围内，受富营养化影响的沿海地区有500多个。此外，许多世界淡水湖泊、溪流和水库遭受富营养化。相关研究表明，我国当前已有80％的湖泊出现富营养化现象。在我国除一些人烟稀少的淡水湖区域以及还保持着原始状态的湖泊之外，其他湖泊营养盐水平的浓度已经超过了湖泊富营养化的标准浓度。尤其是我国的几大淡水湖的营养盐浓度都已经超过了磷、氮富营养化的发生浓度，总氮浓度甚至超过10倍之上。当下，我国的太湖以及巢湖都已经完全变为富营养化状态。因此，水环境问题已迫不容缓。水质预测能有效的掌握水质情况，预防大规模水污染事件的发生，对水污染防治具有现实意义。所以实现合理有效的水质预测，并在此基础上评估富营养化风险是解决水环境问题的重要思路。

已有的水质预测方法主要分为两类：基于机理的水质预测方法和数据驱动的水质预测方法。(1)基于机理的水质预测方法。这类方法主要是通过分析水体内部的机理演化过程来建立水质模型，再基于该水质模型对水质指标进行预测。常用的机理模型包括：QUAL2E、WASP和 EFDC等。QUAL2E是一个一维水质模型，适用于模拟完全混合的枝状河流水质；它假定存在主流输送机理，即假定平流与扩散混合都是沿着河流的主流向，而在河流的横向与垂向上水质组分是完全均匀混合的，允许河流沿程有多个污染源、取水口以及支流汇入；其基本方程是一维平流—扩散质量迁移方程。WASP(The water quality analysissimulation program) 是美国环保局(EPA)开发并推荐使用的用于各类水体水质模拟与分析的软件；它可用于模拟水动力学、河流一维不稳定流、湖泊和河口三维不稳定流、常规污染物和有毒污染物在水中的迁移和转化规律；其基本原理是平移—扩散方程。EFDC(TheEnvironmental Fluid Dynamics Code)模型是涵盖一维到三维地表水水质的数学模型；它可实现河流、湖泊、水库、湿地系统、河口和海洋等水体的水动力学、水质模拟、污染物迁移和泥沙输移等，是一个多参数有限差分模型。(2)数据驱动的水质预测方法。这类方法利用数据分析模型和工具对水质指标的实测数据进行建模，再通过此模型对水质指标进行预测。常用的数据驱动模型包括：人工神经网络、时间序列模型和支持向量机等，该类模型均是利用不同时刻的水质指标观测数据建立表征水质演化过程的数据驱动模型，并在此基础上实现水质过程预测。两类方法各具有一定优势，然而已有方法大都得到的是水质指标的取值变化预测(单值预测)，其结果通常具有一定的偶然性，精确度和可信度难以令人信服，在此基础上的富营养化评价也会存在偶然性和不确定性。

发明内容

本发明为了弥补单值预测方法带来的不足，结合水质动力学模型，提出了一种基于蒙特卡罗仿真的河湖水质预测及富营养化风险评估方法，为水质预测和富营养化评估的不精确性问题提供一种新思路和可行解决方案。

本发明基于蒙特卡洛仿真的河湖水质预测方法及富营养化风险评估方法，通过以下步骤实现：

步骤一：基于水质指标，DO、Phyt、BOD、OP、PO₄、ON、NH₃-N和NO₃-N，建立水质动力学模型；DO、Phyt、BOD、OP、PO₄、ON、NH₃-N和NO₃-N分别为溶解氧、浮游植物、碳生化需氧量、有机磷、正磷酸盐、有机氮、氨氮和硝酸盐氮。

步骤三：选取Chl_a、总磷TP和总氮TN三个水质指标，计算三个水质指标的营养状态指数，通过权值归一化得出综合营养状态指数TLI；对TLI做分级处理，并计算得出TLI的概率分布和处于不同营养程度的概率；上述Chl_a为Phyt中叶绿素的浓度；TP为OP与PO₄；TN为ON、NH₃-N>3-N。

步骤二：利用蒙特卡罗仿真预测水质指标的演化过程，获得各水质指标在不同时间点上取值的概率分布。蒙特卡罗仿真的基本原理是当问题或对象本身具有概率特征时，可以用计算机模拟的方法产生抽样结果，根据抽样计算统计量或者参数的值；随着模拟次数的增多，可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。在基于蒙特卡罗仿真的预测过程中，水质指标演化过程会被多次模拟。对每一次模拟，每一个水质指标在每一个预测时间点的数值均不相同，这些数值同时具有规律性和随机性。规律性体现在预测值是基于确定的水质动力学模型、水质指标和模型参数初值而产生的；随机性体现在每一个水质指标在每一个预测时间点受噪声影响而得到的预测值不尽相同。当模拟次数充分大时，结合概率统计规律，可以获得各水质指标在不同时间点上取值的概率分布。

步骤三：富营养化风险评估方法的构建

选取Chl_a、总磷TP和总氮TN三个水质指标，计算三个水质指标的营养状态指数，通过权值归一化得出综合营养状态指数TLI；对TLI做分级处理，并计算得出TLI的概率分布和处于不同营养程度的概率。上述Chl_a为Phyt中叶绿素的浓度；TP为OP与PO₄；TN为ON、NH₃-N>3-N。

本发明的优点在于：

1)本发明基于蒙特卡洛仿真的河湖水质预测及富营养化风险评估方法，利用蒙特卡罗仿真预测8个水质指标的演化过程，获得8个水质指标取值的概率分布，实现水质预测。弥补了单一定值的水质预测的不精确性，降低了预测的偶然性。

2)本发明基于蒙特卡洛仿真的河湖水质预测及富营养化风险评估方法，通过构造综合营养状态指数，结合水质指标预测结果，计算综合营养状态指数的概率分布和处于不同营养程度的概率，实现富营养化风险评估。提高了水体的富营养化风险评估的准确度。

附图说明

图1是本发明整体流程图；

图2是水体生态系统演化过程图；

图3a是本发明所涉及的水质指标中溶解氧DO演化曲线图；

图3b是本发明所涉及的水质指标中浮游植物Phyt演化曲线图；

图3c是本发明所涉及的水质指标中碳生化需氧量BOD演化曲线图；

图3d是本发明所涉及的水质指标中有机磷OP演化曲线图；

图3e是本发明所涉及的水质指标中正磷酸盐PO₄演化曲线图；

图3f是本发明所涉及的水质指标中有机氮ON演化曲线图；

图3g是本发明所涉及的水质指标中氨氮NH₃-N演化曲线图；

图3h是本发明所涉及的水质指标中硝酸盐氮NO₃-N演化曲线图；

图4a是本发明所涉及的水质指标中不同时刻下溶解氧DO的概率密度函数图；

图4b是本发明所涉及的水质指标中不同时刻下浮游植物Phyt的概率密度函数图；

图4c是本发明所涉及的水质指标中不同时刻下碳生化需氧量BOD的概率密度函数图；

图4d是本发明所涉及的水质指标中不同时刻下有机磷OP的概率密度函数图；

图4e是本发明所涉及的水质指标中不同时刻下正磷酸盐PO₄的概率密度函数图；

图4f是本发明所涉及的水质指标中不同时刻下有机氮ON的概率密度函数图；

图4g是本发明所涉及的水质指标中不同时刻下氨氮NH₃-N的概率密度函数图；

图4h是本发明所涉及的水质指标中不同时刻下硝酸盐氮NO₃-N的概率密度函数图；

图5是应用本发明方法获得的不同时刻下综合营养状态指数TLI概率密度函数图；

图6a是应用本发明方法，预测天数为100天的贫营养概率曲线图；

图6b是应用本发明方法，预测天数为100天的中营养概率曲线图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是基于蒙特卡洛仿真的河湖水质预测方法及富营养化风险评估方法，如图1所示，具体通过如下步骤实现：

步骤一：建立水质动力学模型。

水体生态系统演化过程，如图2所示，水体生态系统主要分为以下四个过程：溶解氧平衡过程、浮游植物动力学过程、磷循环和氮循环。水体生态系统主要由八个水质指标描述其演化过程：溶解氧(DO)、浮游植物(Phyt)、碳生化需氧量(BOD)、有机磷(OP)、正磷酸盐(PO₄)、有机氮(ON)、氨氮(NH₃-N)和硝酸盐氮(NO₃-N)。基于上述水质指标演化机理，针对水体生态系统的富营养化过程建立水质动力学模型，其数学模型为上述8个水质指标方程组成的一组非线性连续微分方程组。

其中，{C_DO,C_Phyt,C_BOD,C_OP,C_ON,}分别表示水质指标DO、Phyt、BOD、OP、 PO₄、ON、NH₃-N和NO₃-N的浓度。Phyt的浓度以叶绿素a(Chl_a)的浓度为准，参数{k₁,k₂,…,k₁₃}>

表1水质动力学模型参数

令θ＝(k₁，k₂，k₃，k₄，k₅，k₆，k₇，k₈，k₉，k₁₀，k₁₁，k₁₂，k₁₃)^T，且假设模型参数是慢变的，则式(1)表示为：

上式中，f(x,θ)为水质动力学模型函数。

步骤二：基于蒙特卡洛仿真的河湖水质预测方法的构建。

将式(2)写成离散化形式如下：

上式中，x∈R^8×1,θ∈R^13×1,f:R^8×1×R^13×1→R^8×1；水质指标的过程噪声模型参数过程噪声其中，Σ_x、Σ_θ为相应的噪声协方差阵；h是预测步长，h≤1day，t为时间点，t＝0,h,2h,…,T，T为预测天数。

在基于蒙特卡罗仿真的预测过程中，水质指标演化过程会被多次模拟。对每一次模拟，每一个水质指标在每一个预测时间点的数值均不相同，这些数值同时具有规律性和随机性。规律性体现在预测值是基于确定的水质动力学模型、水质指标和模型参数初值而产生的；随机性体现在每一个水质指标在每一个预测时间点受噪声影响而得到的预测值不尽相同。当模拟次数充分大时，结合概率统计规律，可以获得各水质指标在不同时间点上取值的概率分布。基于蒙特卡罗仿真的水质预测的具体流程步骤为：

A、令水质指标的初值为x₀，模型参数的初值为θ₀，蒙特卡罗仿真粒子数为N，预测步长为h，初始时刻t＝0。

B、根据水质指标和模型参数的初值(x₀,θ₀)，对N个粒子赋初值，即i为第i个粒子，令

x₀⁽ⁱ⁾＝x₀,θ₀⁽ⁱ⁾＝θ₀>

C、对于时刻t＝h,2h,…,T,基于t-1时刻的水质指标和模型参数值{x_t-1⁽ⁱ⁾,θ_t-1⁽ⁱ⁾}，利用式(3)>t⁽ⁱ⁾,θ_t⁽ⁱ⁾}。

D：若t＜T，t←t+h，返回进行步骤C，否则进行步骤E。

E：对计算水质指标和模型参数的概率密度函数；

上式中，P(x_t|(x₀,θ₀))为基于(x₀,θ₀)的水质指标概率密度函数，P(θ_t|(x₀,θ₀))为基于(x₀,θ₀)的模型参数概率密度函数，δ为狄拉克函数。

步骤三：富营养化风险评估方法的构建。

由于水体生态系统主要由八个水质指标描述其演化过程：溶解氧(DO)、浮游植物(Phyt)、碳生化需氧量(BOD)、有机磷(OP)、正磷酸盐(PO₄)、有机氮(ON)、氨氮(NH₃-N)和硝酸盐氮>3-N)。选取Chl_a、总磷(TP,TP＝OP+PO₄)和总氮(TN,TN＝ON+NH₃-N+NH₃-N)三个水质指标，计算三个水质指标的营养状态指数，通过权值归一化得出综合营养状态指数TLI；对TLI做分级处理，并计算得出TLI的概率分布和处于不同营养程度的概率。具体流程步骤为：

a、对于时刻t＝h,2h,…,T，粒子x⁽ⁱ⁾，i∈[1,N]：

(1)对于Chl_a，C_Phyt,t∈x_t⁽ⁱ⁾，有

TLI_t(Chl_a)＝10(2.500+1.086ln(C_Phyt,t))>

对于TP，C_OP,t、有

TLI_t(TP)＝10(9.436+1.624ln(C_OP,t+C_PO4,t))>

对于TN，C_ON,t、有

(2)以Chl_a为基准参数，则Chl_a、TP、TN的权重分别为：

上式中，则w₁、r₁₁分别是Chl_a的权重和与自身的相关系数，w₂、r₁₂分别是TP的权重和与基准参数Chl_a的相关系数，w₃、r₁₃分别是TN的权重与基准参数Chl_a的相关系数。

(3)结合各水质指标的营养状态指数和权重，计算TLI：

TLI_t＝w₁·TLI_t(Chl_a)+w₂·TLI_t(TP)+w₃·TLI_t(TN)>

b、：对TLI做分级处理。TLI的值越大营养程度越重，以0—100分级，TLI<30为贫营养 (Oligotrophic)，30<TLI<50为中营养(Mesotrophic)，TLI>50为富营养(Eutropher)。

c、计算TLI的概率密度函数：

d、根据综合营养状态指数TLI的概率分布计算处于不同营养程度的概率

上式中，为t时刻预测结果的TLI_t属于区间[E_m,E_m+1]的概率，B[E_m，E_m+1]为TLI_t属于区间[E_m,E_m+1]的个数，区间[E_m,E_m+1]为TLI_t的分级区间，E_m与E_m+1分别为指TLI_t的第m个分级区间的上限和下限，即当属于第m个分级区间时，TLI_t>m＜TLI_t＜E_m+1。

实施例1：

步骤一：建立水质动力学模型

基于水质指标演化机理，建立水质动力学模型，模型方程如式(1)所示。

步骤二：基于蒙特卡洛仿真的河湖水质预测方法的构建

水质指标初值：x₀＝[6.5,0.05019,3,0.07,0.06,0.35,0.36,0.1]^T，

水质指标过程噪声协方差阵：Σ_x＝diag{0.1,0.001,0.01,0.001,0.001,0.01,0.01,0.01}，

模型参数初值：θ₀＝[0.1,1.208,0.01,0.01,0.15,0.13,0.0009,0.001,0.025,0.075,0.9,0.35,1.02,1.3]^T，

模型参数过程噪声协方差阵：

Σ_θ＝diag{0.1,0.01,0.001,0.1,0.001,0.001,0.0001,0.0001,0.0001,0.01,0.001,0.01,0.1}。

预测步长h＝1/24(单位：day)、粒子数N＝1000和仿真时长T＝100(单位：day)。

按照式(1)～(4)和上述参数数值计算8个水质指标演化结果，结果如图3a～3h所示。按照式(5)计算不同时刻下的8个水质指标的概率密度函数，结果如图4a～4h所示。

步骤三：富营养化风险评估方法的构建

基于上述水质指标演化结果，构造综合营养状态指数，计算综合营养状态指数的概率分布和处于不同营养程度的概率，实现富营养化风险评估。选取Chl_a(即Phyt)、总磷(TP) 和总氮(TN)三个水质指标，计算三个水质指标的营养状态指数，通过权值归一化得出TLI；对TLI做分级处理，并计算得出TLI的概率分布和处于不同营养程度的概率。

中国湖泊(水库)的基准参数Chl_a与其他参数之间的相关系数r_1j如下表2所示。

表2中国湖泊(水库)的基准参数Chl_a与其他参数之间的相关关系r_ij

按照式(6)～(11)和表2计算综合营养状态指数的概率密度函数，结果如图5所示，按照式(12)计算水体处于不同营养状态的概率，结果如图6a、6b所示。