CN110889201B

CN110889201B - 基于混杂动态系统的水质机理建模及健康风险评估方法

Info

Publication number: CN110889201B
Application number: CN201911072914.8A
Authority: CN
Inventors: 赵峙尧; 王小艺; 周宇琴; 王立; 许继平; 于家斌; 白玉廷
Original assignee: Beijing Technology and Business University
Current assignee: Beijing Technology and Business University
Priority date: 2019-11-05
Filing date: 2019-11-05
Publication date: 2020-07-03
Anticipated expiration: 2039-11-05
Also published as: CN110889201A

Abstract

本发明公开了一种基于混杂动态系统的水质机理建模及健康风险评估方法，属于环境工程技术领域。本发明包括：首先构建一般性水质机理模型，并为弥补一般性水质机理模型在有效建模水质演化过程的不足，利用混杂动态系统，构建基于混杂动态系统的水质机理模型；在已知水质测量数据的基础上，利用果蝇优化算法估计模型的未知参数；利用改进的交互多模型实现水质演化状态估计，结合健康度，实现水质健康风险评估。本发明建立了能有效描述水体不同状态和水质退化过程的水质机理模型，并能准确估计水质机理模型的未知参数，与已有方法相比，提高了水质监测的精度和速度，能有效实现水质评价，且考虑更加全面、准确，克服了单值评价带来的偶然性。

Description

基于混杂动态系统的水质机理建模及健康风险评估方法

技术领域

本发明属于环境工程技术领域，具体涉及一种基于混杂动态系统的水质机理建模及健康风险评估方法。

背景技术

水是地球上生命的基础，但这种宝贵的资源越来越受到威胁。在中国，据《2018年中国环境公报》可知，111个监测营养状态的湖泊(水库)中，I类水质的湖泊(水库)占6.3％；Ⅱ类水质的湖泊(水库)占30.6％；Ⅲ类水质的湖泊(水库)占29.7％；Ⅳ类水质的湖泊(水库)占17.1％；Ⅴ类水质的湖泊(水库)占8.1％；劣Ⅴ类水质的湖泊(水库)占8.1％。因此，水环境问题已迫不容缓。水质机理建模能有效描述水体内部环境变化，水质评价能有效评估水质情况，二者结合能有效描述水质健康状况，对水污染防治具有现实意义，这是解决水环境问题的重要思路。

目前，已有的水质建模方法主要分为两类：基于机理的水质建模方法和数据驱动的水质建模方法。(1)基于机理的水质建模方法。这类方法主要是通过分析水体内部的机理演化过程来建立水质模型。已有的经典水质机理模型包括：MIKE、WASP和EFDC等。(2)数据驱动的水质建模方法。这类方法利用数据分析模型和工具对水质指标的实测数据进行建模。常用的数据驱动模型包括：人工神经网络、时间序列模型和支持向量机等。两类方法各有优势，然而也各有缺点。首先，数据驱动的水质建模效果取决于数据量的数量和质量，且无法合理解释水质指标演化过程。其次，已有的经典水质机理模型通常包含大量未知参数，仅给定经验取值范围，而水质参数决定了水质模型的精确度，在此基础上预测水质变化或评价水质健康状态难以满足精度要求。最后，水质指标，例如叶绿素、总氮和总磷等，会随温度、光照等自然因素和过度耕种、过度使用化肥等人为因素的变化而变化，且自然因素和人为因素均同时具有规律性(自然行为和人为行为都存在周期性)和随机性，这导致湖泊(水库)水体演化过程也具有规律性和随机性，因此同一湖泊(水库)在不同时期内的水质模型参数不尽相同，即单一特性的水质演化机理模型(一般性水质机理模型)无法有效描述水体长时间内的演化过程。

发明内容

本发明针对已有的水质建模存在的问题，提供了一种基于混杂动态系统的水质机理建模及健康风险评估方法。首先，本发明为了弥补一般性水质机理模型在有效建模水质演化过程的不足，即无法有效描述水体不同状态和水质退化过程的问题，提出基于混杂动态系统的水质机理建模方法；其次，为了准确估计水质机理模型的模型参数，提出基于果蝇优化算法(Fruit fly optimization algorithm，FOA)的参数辨识方法，最后，为了弥补单值水质评价方法的偶然性，利用改进的交互多模型实现水质演化混杂状态估计，结合健康度，实现水质健康风险评估，为水质建模和评价提供一种新思路。

本发明的一种基于混杂动态系统的水质机理建模及健康风险评估方法，包括如下步骤：

步骤一：基于水质指标建立一般性水质机理模型；所述的水质指标包括溶解氧DO、浮游植物Phyt、碳生化需氧量BOD、有机磷OP、正磷酸盐PO₄、有机氮ON、氨氮NH₃-N和硝酸盐氮NO₃-N。

步骤二：对待研究的水体，从水体环境分析获得M个离散模态，建立基于混杂动态系统的水质机理模型，简称混杂水质机理模型，以描述水体不同状态和水质退化过程；

混杂水质机理模型表示为

其中，

表示离散模态空间，n表示给离散模态

分配连续过程变量空间，F表示连续动态方程集合，Π表示模态转移概率矩阵，Init表示混杂水质机理模型的初始条件，由初始时刻的混杂状态的分布来描述。混杂水质机理模型的混杂状态表示为s＝(q,x)，x为连续过程变量，是由水质指标组成的向量。

对混杂水质机理模型中的模型参数通过参数敏感性分析方法计算出敏感性参数和非敏感性参数。

步骤三：基于果蝇优化算法辨识水质参数；混杂水质机理模型中的敏感性参数是未知参数，根据待研究的水体的历史数据，通过果蝇优化算法求解所有模态的未知参数的最优值；果蝇优化算法中的味道浓度根据水质指标在各时刻的误差来设置，误差越小，味道浓度越佳，代表未知参数越优。

步骤四：利用步骤二输出的混杂水质机理模型，基于改进交互多模型算法对待研究的水体的水质健康风险进行评估。

本发明与现有技术相比，具有以下优势：

(1)本发明建立基于混杂动态系统的水质机理模型，弥补了一般性水质机理模型在有效建模水质演化过程的不足，更加符合实际水体水质的演化过程。

(2)本发明利用果蝇优化算法估计混杂水质机理模型的模型参数，克服了已有的经典水质机理模型通常包含大量未知参数，且仅给定经验取值范围的问题，有效提高了未知模型参数的辨识效率。

(3)本发明利用改进的交互多模型实现水质模型的状态估计，并结合健康度实现水质健康风险评估，克服了单值水质评价方法的偶然性，可实现对水体长时间的演化过程中的水质健康监测。

附图说明

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

图2是本发明基于混杂动态系统的水质机理模型图；

图3是本发明基于果蝇优化算法辨识水质参数过程的8个水质指标演化曲线和真实值曲线图；

图4是本发明基于果蝇优化算法辨识水质参数过程的水质演化误差减小曲线图；

图5是本发明基于果蝇优化算法辨识水质参数过程的离散模态切换结果图；

图6是本发明基于改进交互多模型算法的水质健康风险评估算法步骤图；

图7是本发明基于改进交互多模型算法的水质健康风险评估的8个水质指标估计值曲线和真实值曲线图；

图8是本发明基于改进交互多模型算法的水质健康风险评估过程的模态转换图；

图9是本发明基于改进交互多模型算法的水质健康风险评估过程的模态概率图；

图10是本发明基于改进交互多模型算法的水质健康风险评估结果。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图对本发明作进一步的详细和深入描述。

本发明的基于混杂动态系统的水质机理建模及健康风险评估方法，首先结合水质演化机理，构建水质演化机理模型，即一般性水质机理模型；其次，为弥补了一般性水质机理模型在有效建模水质演化过程的不足，利用混杂动态系统，构建基于混杂动态系统的水质机理模型；并在已知水质测量数据的基础上，利用果蝇优化算法估计模型的未知参数；最后，利用改进的交互多模型实现水质演化状态估计，结合健康度，实现水质健康风险评估。本发明建立能有效描述水体不同状态和水质退化过程的混杂水质机理模型，且果蝇优化算法能够准确估计混杂水质机理模型的未知参数；水质健康风险评估方法能有效估计混杂水质机理模型的混杂状态，实现水质评价，且考虑更加全面、准确，克服了单值评价带来的偶然性。

本发明是基于混杂动态系统的水质机理建模及健康风险评估方法，通过如下步骤实现：

步骤一：建立一般性水质机理模型。

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

其中，C_DO,C_Phyt,C_BOD,C_OP,

C_ON,

分别表示水质指标DO、Phyt、BOD、OP、PO4、ON、NH3-N和NO3-N的浓度；Phyt浓度的由叶绿素a的浓度表征；参数k₁,k₂,…,k₁₃为模型未知参数，其含义和单位如表1所示。

表1水质动力学模型参数

符号	含义	单位
			k<sub>1</sub>	20℃复氧速率	day<sup>-1</sup>
k<sub>2</sub>	20℃还原速率	day<sup>-1</sup>
			k<sub>3</sub>	20℃硝化速率	day<sup>-1</sup>
k<sub>4</sub>	底泥需氧量	g/(m<sup>2</sup>·day)
			k<sub>5</sub>	浮游植物生长速率	day<sup>-1</sup>
k<sub>6</sub>	浮游植物死亡速率	day<sup>-1</sup>
			k<sub>7</sub>	20℃反硝化速率	day<sup>-1</sup>
k<sub>8</sub>	有机物沉降速率	m/day
			k<sub>9</sub>	溶解有机氮矿化速率	day<sup>-1</sup>
k<sub>10</sub>	有机磷溶解比例	None
			k<sub>11</sub>	有机氮溶解比例	None
k<sub>12</sub>	浮游植物呼吸温度系数	None
			k<sub>13</sub>	氨优先选择系数	None

表中，None表示没有单位。

令

θ＝(k₁,k₂,k₃,k₄,k₅,k₆,k₇,k₈,k₉,k₁₀,k₁₁,k₁₂,k₁₃)^Τ

其中，x表示水质指标组成的向量，θ表示模型未知参数组成的向量，且假设模型参数是慢变的，θ中任意参数的微分为0，则公式(1)可表示为：

上式中，f(x,θ)是水质机理模型函数，

和

分别为x和θ的一阶微分，0为数值0的向量。

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

其中，x(k)和θ(k)分别为k时刻的水质指标和模型参数；x(k-1)、θ(k-1)、ω_x,k-1和ω_θ,k-1分别为k-1时刻的水质指标、模型参数、水质指标过程噪声和模型参数过程噪声；h_d(h_d≤1day)是预测步长。水质指标的过程噪声

模型参数过程噪声

其中

是高斯分布，Σ_x、Σ_θ分别为水质指标和模型参数的过程噪声协方差阵。

则f通过8个水质指标和13个水质模型参数描述了8个水质指标的水质演化过程，

其中，

是实数集，

分别代表8×1维的实数矩阵空间和13×1维的实数矩阵空间。

步骤二：建立基于混杂动态系统的水质机理模型。

混杂系统是由“连续性子系统”和“离散性子系统”相互作用而构成的一类系统。当“连续性子系统”包含动态特性时，则称该类混杂系统为混杂动态系统(Hybrid dynamicalsystem,HDS)。在混杂动态系统中，连续性子系统的动态特性会基于系统动态模型随着时间发展不断的演化；离散性子系统则会受到外部事件和内部自发事件的影响产生状态的变化。二者相互作用，使系统在整体上呈现离散位置的变化，在局部上呈现连续位置的演化，表现出相比于连续动态系统和离散事件系统更加复杂的动态行为。目前，混杂动态系统已广泛应用于工业过程控制和运载器运动控制中。

本发明基于一般性水质机理模型，结合混杂动态系统，建立能描述水体不同状态和水质退化过程的水质机理模型，即混杂水质机理模型

水质指标x随不同的模型参数θ变化，将x作为连续过程变量，将不同模型参数对于那个不同离散模态，建立混杂水质机理模型

具体如下：

1)设本发明的水质监测中混杂动态系统的离散模态的个数为M，其中第i个离散模态为q_i，i＝1,2,…,M，M为正整数。不同离散模态对应不同的水体状态，每个离散模态对应一组模型参数θ。

2)定义离散模态空间

表示为

“：＝”表示定义的意思。

3)

表示给离散模态q分配连续过程变量空间，

离散模态q的连续过程变量空间表示为

定义在混杂状态空间上的混杂状态

表示混杂水质机理模型空间。

4)F＝[f₁,f₂,…,f_M]表示连续动态方程集合，f_i＝f(x,θⁱ)表示离线模态q_i的连续动态方程，x为连续过程变量，θⁱ为离散模态q_i的模型未知参数θ；考虑到每个模态下水质机理模型的模型未知参数θ个数为13，则具有M个模态的混杂水质机理模型

未知参数个数为M×13。这意味着模态的增加，提高了模型未知参数辨识的难度，影响了辨识的准确度。因此，通过参数敏感性分析方法计算水质机理模型f(x,θ)的模型敏感性参数(简称敏感参数)α和模型非敏感性参数β的取值，即θⁱ＝(αⁱ,β)，其中设敏感参数个数为η。对模态q_i，由式(3)，存在f_i(x,αⁱ)＝f(x,θⁱ)＝f(x,αⁱ,β)，即：

其中，对于q_i模态，

模型未知参数过程噪声

为模型敏感参数的过程噪声协方差阵；αⁱ(k-1)和

分别为模态q_i在k-1时刻的模型未知参数和模型未知参数过程噪声；αⁱ(k)为模态q_i在k时刻的模型未知参数。

5)

为模态转移概率矩阵，转移概率π_ij为：

其中，

表示k-1时刻的模态q_i转移到k时刻的模态q_j的概率；且

根据马尔科夫性，系统在k时刻的状态仅取决于系统在k-1时刻的状态，而与先前时刻无关。设

表示混杂水质机理模型的离散模态概率向量，则有：

p(k)＝p(k-1)·Π (7)

其中，p₁为离散模态q₁的概率；p₂为离散模态q₂的概率；p_M为离散模态q_M的概率。

6)Init是系统的初始条件。设

表示概率密度函数。混杂水质机理模型的初始条件可由混杂状态s(0)＝(q(0),x(0))的分布来描述。一般地，第j个模态q_j的初始概率密度函数

为

其中，q_j(0)和x(0)分别为初始模态和初始状态，μ_j(0)为初始时刻x的均值，Σ_j(0)初始时刻x的方差。对于

满足p_j(0)≥0，且

步骤三：基于果蝇优化算法辨识水质参数。

一般水质机理模型包含有大量未知参数，仅给定经验取值范围，而水质参数决定了水质模型的精确度，在此基础上预测水质变化或评价水质健康状态难以满足精度要求。本发明建立的混杂水质机理模型

的复杂程度高于一般水质机理模型，这提高了水质参数的辨识难度。对此本发明基于果蝇优化算法辨识水质参数，从而实现对混杂水质机理模型

中参数的辨识。

果蝇优化算法是一种基于果蝇觅食行为来寻求全局优化的新方法，属于仿生智能优化算法。果蝇是一种感官知觉优于其他物种的昆虫，尤其是嗅觉和视觉。果蝇的嗅觉器官能很好的搜集漂浮在空气中的各种气味，甚至能够嗅到40公里以外的食物源。然后，飞到食物位置附近后亦可使用敏锐的视觉发现食物和同伴聚集的位置，并且向该方向飞去。果蝇优化算法在寻优过程中具有一定的随机性，为了确保果蝇群体向着正确的方向飞行，该算法引入了气味浓度判定值和气味浓度判定函数。果蝇优化算法具有计算过程简单、全局收敛性强、执行的时间短等优点，可应用于求解最优解，在求解数学函数极值、数据挖掘、广义回归神经网络参数优化等领域已有初步成果。

本发明利用果蝇优化算法，估计基于混杂动态系统的水质机理模型的未知参数。在基于果蝇优化算法的水质机理建模过程中，将寻优过程中的参数代入味道浓度函数中，获得该组参数的味道浓度Smell。在果蝇种群中，对比每个果蝇的味道浓度，并找出最优的味道浓度和该果蝇位置(最佳位置)，其余果蝇利用视觉飞往最佳位置。最后在不断迭代过程中获得最优参数解。

本发明基于果蝇优化算法的水质机理建模方法的一个实现流程如表2所示：

表2基于果蝇优化算法的水质机理建模方法

通过果蝇优化算法获得各个模态q_i的敏感参数的最优值

完成混杂水质机理的建模。

为验证基于果蝇优化算法的水质机理建模的有效性，采用误差最大值MaxE、误差均值MAE、误差方均根RMSE为验证指标，计算如下：

对第ic个水质指标x_ic∈x，基于果蝇优化算法的水质机理建模方法的MaxE为MaxE₁(ic)、MAE为MAE₁(ic)、RMSE为RMSE₁(ic)有

步骤四：基于改进交互多模型算法的水质健康风险评估。

交互多模型(Interacting Multiple Model,简称IMM)是一种基于滤波的递归估计器，能够有效估计混杂动态系统的混杂状态(离散的系统模态和连续的过程变量)分布，具有自适应的特点。交互多模型并行滤波器可采用滤波算法，例如粒子滤波算法、卡尔曼算法和扩展卡尔曼算法等。

本发明基于混杂动态系统的水质机理模型，采用交互多模型算法估计水质演化状态，并结合健康度，实现水质健康风险评估。本发明实施例中，系统模型的初始混杂状态分布如公式(8)，选用扩展卡尔曼滤波算法为交互多模型的并行滤波算法，则基于改进交互多模型的水质健康风险评估方法的一个实现流程如表3所示。

表3基于改进交互多模型算法的水质健康风险评估方法

设k时刻8个水质指标估计量为

8个水质指标观测量为y(k)，8个水质指标误差为Error₂(k)：

采用误差最大值、误差均值、误差方均根为验证指标，验证本发明基于改进交互多模型算法的水质健康风险评估方法的有效性。对第ic个水质指标x_ic∈x，误差Error₂(ic)∈Error₂，且Error₂(ic,k)∈Error₂(ic)，误差最大值MaxE₂(ic)、误差均值MAE₂(ic)、误差方均根RMSE如下：

实施例

步骤一：建立一般性水质机理模型；

基于水质指标演化机理，建立一般性水质机理模型，模型方程如公式(1)所示。

步骤二：建立基于混杂动态系统的水质机理模型；

基于一般性水质机理模型，结合混杂动态系统，建立基于混杂动态系统的水质机理模型

通过研究水质状态演化机理设定离散模态个数M＝3，即

基于混杂动态系统的水质机理模型如图2所示；设定过程变量为8个水质指标。利用参数敏感性分析确定单一模态下水质机理模型的非敏感性模型参数为(k₁,k₃,k₄,k₆,k₉,k₁₁)，如表4所示，敏感性模型参数个数Ms＝7，敏感性模型参数为α^Ms＝(k₂,k₅,k₇,k₈,k₁₀,k₁₂,k₁₃)^T。

表4非敏感性模型参数取值

设定模态转移概率矩阵Π：

步骤三：基于果蝇优化算法辨识水质参数；

本实施例中，水质指标数据来源于苏州白洋湾金墅站点2009年6月21日至2011年6月21日730天的水质数据，采样频率为每天1次，即h_d(h_d＝4hour＝1/6day)，共N_D＝730组数据。

设置果蝇优化算法中的参数：

最大迭代次数MaxGeneration＝500，

果蝇种群个数MaxGroup＝700，

味道浓度阈值ObjSmell＝10；

按照表4参数及上述参数计算8个水质指标演化结果，结果如图3所示，误差减小曲线如图4所示，水质指标演化误差为4.102，离散模态切换结果如图5所示，其优化参数如表5所示，8个水质指标中的误差最大值MaxE、误差均值MAE、误差方均根RMSE的结果如表5所示。

表5果蝇优化算法的优化参数

参数	模态q<sub>1</sub>	模态q<sub>2</sub>	模态q<sub>3</sub>
				k<sub>2</sub>	0.346887	2.061147	0.322873
k<sub>5</sub>	0.184434	0.303521	0.386834
				k<sub>7</sub>	0.404857	0.492667	0.352939
k<sub>8</sub>	0.594912	0.492328	0.550249
				k<sub>10</sub>	1.193282	0.542891	0.592165
k<sub>12</sub>	0.591601	1.191702	0.281865
				k<sub>13</sub>	1.283096	0.544441	0.651737

混杂水质机理模型建模结果如图3所示。在图3中，水质指标的真实值和混杂水质机理模型的模拟值被分别描述。为进一步反映优化结果的准确性，混杂水质机理模型的误差减少曲线如图4所示；模态切换如图5所示。同时，通过FFOA获得混杂水质机理模型的最优敏感参数如表5所示。最后，选用指标MAXE、MAE、RMSE为混杂水质机理模型的建模过程评价指标，计算结果如表6所示。

表6各水质指标的误差最大值、误差均值和误差方均根

指标	MAXE	MAE	RMSE
				DO	0.902163	0.05756	0.101157
Phyt	0.146091	0.006429	0.016872
				BOD	0.603974	0.034923	0.068712
OP	0.003825	0.00038	0.000531
				PO<sub>4</sub>	0.005953	0.000546	0.00076
ON	0.024295	0.004711	0.005911
				NH<sub>3</sub>-N	0.026744	0.006986	0.008285
NO<sub>3</sub>-N	0.040819	0.004755	0.005818

由图3-图4可知，水体长期演化过程中存在不同水质状态，即存在不同特征的水体特性，采用混杂水质机理模型能有正确描述水体长期演化过程中的具有不同特征的水体特性和不同水体特性下的水质演化过程。这一点也反映在图5中，混杂水质机理模型的模态在实验时间内均发生了变化，这说明水体内部演化过程具有不同特性。最后，8个水质指标的MAXE、MAE、RMSE指标计算值和指标综合计算值如表6所示。在表6中，混杂水质机理模型的三个评价指标计算值均较小，这进一步反应了混杂水质机理模型能有正确描述水体长期演化过程中的具有不同特征的水体特性和不同水体特性下的水质演化过程。

步骤四：基于改进交互多模型算法的水质健康风险评估；

基于改进交互多模型算法的水质健康风险评估步骤如图6所示。水质指标数据来源于苏州白洋湾金墅站点20011年6月21日至2012年6月21日365天的水质数据，采样频率为每天1次，即h(h＝4hour＝1/6day)，共N_T＝365组数据。

对任意一个模态q_i(i＝1,2,3)：

过程噪声系数矩阵Γ_w,i＝diag{1；1；1；1；1；1；1；1；0；0；0；0；0；0；0；}；

过程噪声方差

观测噪声矩阵Γ_v,i＝I₈，I₈为8×8单位矩阵；

观测噪声方差Q_v,i＝diag{0.0006；0.001；0.003；0.00001；0.00001；0.0001；0.00003；0.00001；}；

依据国标GBT14848-2017《地下水质量标准》，结合8个水质指标，设定水质健康标准如表7所示。。

表7水质健康标准

按照表5和上述参数计算8个水质指标演化结果，结果如图7所示，模态交互结果如图8所示，模态概率如图9所示。8个水质指标中的误差最大值MaxE、误差均值MAE、误差方均根RMSE，结果如表8所示。基于上述研究结果和表7的水质健康标准，水质健康风险评估结果如图10所示

表8 IMM的误差最大值、误差均值和误差方均根

指标	MAXE	MAE	RMSE
				DO	0.140945	0.013926	0.02085
Phyt	11.70148	0.663541	1.12886
				BOD	0.777577	0.12303	0.165695
OP	0.023556	0.002062	0.003348
				PO<sub>4</sub>	0.021172	0.001846	0.002986
ON	0.651383	0.041808	0.070607
				NH<sub>3</sub>-N	0.099915	0.008572	0.013417
NO<sub>3</sub>-N	0.001303	0.000117	0.000175

考虑到混杂水质机理模型在建模过程中的优势，本发明选用混杂水质机理模型和改进的交互多模型算法估计水质混杂状态的概率分布，获得水质监测值，包含不同模态的水质指标估计值、混杂水质机理模型的水质指标估计值、不同模态概率和水质综合模态。混杂水质机理模型的水质指标估计值如图7所示。为使结果更清晰，依据不同模态概率计算水质综合模态，结果如图8所示。在水质监测过程中，获得水体状态处于不同模态的概率，结果如图依据图9。同时，选用MAXE、MAE和RMSE为评估过程的评价指标，评价结果如表8所示。最后，基于不同模态的水质指标估计值、不同模态概率和水质健康标准，计算水质健康风险概率，结果如图10所示。

由图7可知，混杂水质机理模型的水质指标估计值和真实值非常接近，这说明该模型能准确估计水质混杂状态的概率分布。这一点也表现在图8-图9中。在图8-图9中，不同模态概率和水质综合模态均表明水质长期演化过程具有不同特性的特点。在表8中，8个水质指标的MAXE、MAE和RMSE表明IMM-EKF算法在水质监测过程的有效性和准确性。最后，结合不同模态的水质指标估计值、不同模态概率和WQHRAS，计算水质健康风险概率，结果如图10所示。在图10中，水质评估结果表示该阶段水体较为健康，水质指标处于健康空间的概率较大。

通过上述的本发明的实施例及具体实施步骤的说明，清楚地说明了采用本发明如何对水体进行监测和水质健康评价的实现过程。对实际中其他水体的监测和水质健康评价，都可以采用本发明方法来实现。在本说明书中的描述中，可能对可在上述参考文献或其它现有技术资料中获得的部分技术细节进行了省略、简化、变通等处理，这对于本领域的技术人员来说是可以理解的，并且这不会影响本说明书的公开充分性。同时，本领域的技术人员能够理解，对本发明的上述实施例能够做出各种修改、变型、以及替换，其均落入如所附权利要求限定的本发明的保护范围内。

Claims

1.一种基于混杂动态系统的水质机理建模及健康风险评估方法，其特征在于，包括：

步骤一：建立水质机理模型，其中，由水质指标组成的向量表示为x，由水质机理模型的未知参数组成的向量表示为θ，水质机理模型函数表示为f(x,θ)；

步骤二：对待研究的水体，将水质机理模型和混杂动态系统结合，建立描述水体不同状态和水质退化过程的混杂水质机理模型；

所述的混杂水质机理模型表示

其中：

(1)对水体环境分析获得M个离散模态，不同离散模态表示水体的不同状态，对应水质机理模型的不同模型参数；将第i个离散模态表示为q_i，离散模态空间

(2)n表示给离散模态q分配连续过程变量空间，

离散模态q的连续过程变量空间表示为

将在混杂水质机理模型的混杂状态空间上的混杂状态表示为s＝(q,x)，x为连续过程变量；

(3)F＝[f₁,f₂,…,f_M]，表示连续动态方程集合，f_i＝f(x,θⁱ)表示离线模态q_i的连续动态方程，由水质机理模型函数得到，θⁱ为离散模态q_i的模型参数向量；

(4)

表示模态转移概率矩阵，转移概率π_ij为：

其中，

表示k-1时刻的模态q_i转移到k时刻的模态q_j的概率；且转移概率满足条件：

表示M×M维的实数矩阵空间；设p＝[p₁,p₂,…,p_M]表示混杂水质机理模型的离散模态概率向量，其中，第i个离散模态q_i的概率表示为p_i；根据马尔科夫性，有p(k)＝p(k-1)·Π，其中，p(k)、p(k-1)分别表示在k时刻和k-1时刻的离散模态概率向量；

(5)Init表示混杂水质机理模型的初始条件，由初始时刻的混杂状态的分布来描述；

对混杂水质机理模型中的模型参数通过参数敏感性分析方法计算出敏感性参数，设θⁱ＝(αⁱ,β),i＝1,2,…,M，其中αⁱ为模型敏感性参数向量，β为模型非敏感性参数向量，设敏感性参数个数为η；

步骤三：基于果蝇优化算法辨识混杂水质机理模型中的水质参数；

混杂水质机理模型中的敏感性参数是未知参数，根据待研究的水体的历史数据，通过果蝇优化算法求解所有模态的未知参数的最优值；果蝇优化算法中的味道浓度根据水质指标在各时刻的误差来设置，误差越小，味道浓度越佳，代表未知参数越优；

步骤四：利用步骤三输出的混杂水质机理模型对待研究的水体的水质健康风险进行评估。

2.根据权利要求1所述的方法，其特征在于，所述的步骤二中，初始时刻的混杂状态表示为s(0)＝(q(0),x(0))，第j个模态q_j的初始概率密度函数表示为

如下：

其中，x(0)表示初始时刻的连续过程变量向量，q_j(0)表示初始时刻的第j个模态，

是高斯分布，μ_j(0)为初始时刻x的均值，Σ_j(0)初始时刻x的方差，

表示第j个模态q_j的初始概率，模态q_i的初始概率p_j(0)满足条件：p_j(0)≥0，且

3.根据权利要求1或2所述的方法，其特征在于，所述的步骤三中，果蝇优化算法求解所有模态的未知水质参数的过程包括：

(a1)初始化，包括：

1.1)设置离散模态个数M、非敏感性模型参数β和单一模态下模型未知参数个数η；

1.2)设置模态转移概率矩阵

1.3)随机设定各离散模态的未知参数初值αⁱ(0)，i＝1,2,…,M；每个模态有η个未知参数；

1.4)设定数据长度为N_D；

1.5)计算目标函数初值，如下：

1.5.1)设定初始模态为q(0)＝q₁；

1.5.2)在k(k＝1,2,…,N_D)时刻，根据p(k-1)、q_i(k-1)和Π，计算p(k)＝p(k-1)·Π，并确定q_j(k)，计算水质指标估计值，对第ic个水质指标，设在k时刻的水质指标实际值为x(ic,k)，计算的估计值为

则误差

得到初始时刻的味道浓度

1.6)保存味道浓度Smell(0)为最佳味道浓度BestSmell；

1.7)设定果蝇种群个体数为MaxGroup，每个个体携带一组未知参数的值；

1.8)随机设定果蝇种群每个个体的初始位置[X_axis,Y_axis]；

1.9)设定最大迭代次数MaxGeneration和味道浓度阈值ObjSmell；

(a2)进行迭代，对第gen(gen＝1,2,…,MaxGeneration)次迭代，执行下面2.1～2.3；

2.1)种群寻优过程，计算每个果蝇的味道浓度，包括2.1.1～2.1.3：

2.1.1)种群中每个果蝇利用嗅觉搜寻食物的随机距离与方向；对第i个模态的第

个敏感参数

更新其在第gro个果蝇的位置：

X(gro,(i-1)×η+im)＝X_axis(gro,(i-1)×η+im)+Random value

Y(gro,(i-1)×η+im)＝Y_axis(gro,(i-1)×η+im)+Random value

Random value表示随机数值；[X_axis(gro,(i-1)×η+im),Y_axis(gro,(i-1)×η+im)]为第i个模态的第im个敏感参数的当前位置；

2.1.2)对第i个模态q_i，计算第gro个果蝇与零点的距离

计算第gro个果蝇的味道浓度判定值S(gro,(i-1)×η+im)＝1/Dist(gro,(i-1)×η+im)；

对第im个模型敏感参数赋值，

2.1.3)计算目标函数：

设定初始模态为q(0)＝q₁；

在k(k＝1,2,…,N_D)时刻，根据p(k-1)、qi(k-1)和Π，计算p(k)＝p(k-1)·Π，并确定q_j(k)，计算水质指标估计值和误差值；设第ic个水质指标的误差值为e₁(ic,k)；得到第gro个果蝇的味道浓度

2.2)对比寻优过程，包括2.2.1～2.2.2：

2.2.1)从第gen次迭代中的所有MaxGroup个果蝇中找到味道浓度最小的果蝇，标记该果蝇的味道浓度为min(Smell)，将该果蝇位置记为最优位置(BestX,BestY)，将该果蝇的误差记为最优误差；

2.2.2)如果min(Smell)＜BestSmell，则更新最佳味道浓度BestSmell＝min(Smell)，重新定位果蝇群体位置：X_axis＝BestX，Y_axis＝BestY；

2.3)判断是否达到最优，包括：

如果BestSmell＜ObjSmell，寻优结束，输出当前最优果蝇；

对第i个模态q_i中第im个模型敏感参数，计算如下：

最优位置到零点的距离

最优位置味道浓度判定值BestS((i-1)×η+im)＝1/BestDist((i-1)×η+im)，

计算最优参数

获得模态q_i的敏感参数的最优值

4.根据权利要求1或2所述的方法，其特征在于，所述的步骤四中，利用步骤三输出的混杂水质机理模型，采用交互多模型算法估计水质演化状态，对水质的健康风险进行评估，包括如下步骤：

步骤1)进行初始化，包括：

1.1)设定数据长度为N_T,k＝1,2,…,N_T，k-1＝0表示初始时刻；

1.2)确定k时刻的水质指标的观测量Y^k＝{y(0),y(1),…,y(k)}；

1.3)根据步骤三输出的最优值设置各模态q_j(j＝1,2,…,M)的模型敏感参数初始值α^j(0)；

1.4)设定模态q_j(j＝1,2,…,M)的混杂水质机理模型的过程噪声系数矩阵Γ_w,j和过程噪声方差Q_w,j；

1.5)设定模态q_j(j＝1,2,…,M)的混杂水质机理模型的观测噪声矩阵Γ_v,j和观测噪声方差Q_v,j；

1.6)设定模态q_j(j＝1,2,…,M)观测量的非线性函数h_j；

步骤2)在时间N_T内的每个时刻k，执行下面过程2.1～2.6：

2.1)对每个模态q_j(j＝1,2,…M)进行估计量交互，包括：

(1)预测模态概率p_j：

Y^k-1为k-1时刻的水质指标的观测量，p_j(k-1)为模态q_j在k-1时刻的概率；

(2)计算交互模态概率

(3)计算交互系统过程变量

E表示求期望；

(4)更新交互协方差矩阵

2.2)对于各模态q_j(j＝1,2,…M)，进行并行滤波，包括：

(1)预测模态q_j在k时刻的状态

(2)f_j表示模态q_j的连续动态方程，计算f_j的Jacobian矩阵A_j和h_j的Jacobian矩阵C_j：

(3)预测模态q_j在k时刻的协方差矩阵：

(4)计算模态q_j在k时刻的测量残差r_j：

(5)计算模态q_j在k时刻的测量残差协方差阵S_j：

(6)计算卡尔曼增益

(7)更新过程变量：

(8)更新协方差矩阵：

2.3)对于各模态q_j(j＝1,2,…M)，更新模态概率和模态识别，如下：

(1)计算似然函数

(2)更新模态概率

(3)设定概率阈值p_T，进行模态识别：

2.4)估计量融合，包括：

(1)计算融合的过程变量

(2)计算融合的协方差矩阵：

2.5)更新转移概率矩阵，如下：

如果k-1时刻处于q_i模态，在k时刻处于q_j模态，且q_i≠q_j，则：

A)定义初等矩阵Ξ如下：

其中，row表示矩阵的每一行，←ith表示第i行元素变化过程，←jth第j行元素变化过程；

B)更新k时刻的转移概率矩阵Π(k)＝Ξ·Π(k-1)·Ξ；

2.6)计算待研究的水体的健康度，包括如下：

首先，计算模态q_j(j＝1,2,…,M)健康度：

A)计算过程变量x_j的概率密度函数

其中，

是高斯分布；

B)设定模型健康空间

和不健康空间

C)计算模态q_j在k时刻的健康度

然后，计算水体在k时刻的健康度