CN103914586A - 基于glue多目标水质模型率定的水质模拟方法 - Google Patents

Abstract

本发明涉及一种基于GLUE多目标水质模型率定的水质模拟方法，其步骤包括：1）根据全局敏感性分析方法对水质模型进行参数的敏感性分析，并根据得到的敏感性指数筛选出多目标敏感参数；2）通过筛选出的多目标敏感参数的先验样本分布情况生成水质的随机参数组；3）对所述随机参数组中的水质指标进行模拟，并通过似然函数进行筛选，通过迭代得到各项水质指标的敏感参数交集，最终使参数收敛至更小区间内，实现对水质的高精度模拟。本发明对多项输出指标同时率定来确定最优参数组合，能够更准确、有效地模拟水质指标,为后续的水环境水质管理、预测、评估等工作提供更为准确的定量化理论支撑。

Description

基于GLUE多目标水质模型率定的水质模拟方法

技术领域

本发明涉及一种基于水质模型的参数自动率定的水质模拟方法，主要用于目前主流的多指标水质模型的率定、校准工作。

背景技术

目前，随着对水环境认知的深入，水环境模型的结构也变得日益复杂。常见的水质模型，例如WASP、QUAL、QUASAR、MIKE（ECOLAB）、BASINS（RCHRES）等，均是可同时模拟多项水质指标的复杂模型系统。尽管这些模型框架逐渐完备，能够更全面地反映人类对水环境改变的认识，但同时也使得模型变量、参数不断增多，导致率定（潘登,任理.2012.分布式水文模型在徒骇马颊河流域灌溉管理中的应用Ⅰ.参数率定和模拟验证[J].中国农业科学,(03):471-479.）难度加大。参数率定是水质模拟的重要环节，而参数率定结果的准确与否，是模型能否正确反映水质变化的规律的关键。

由于参数数量巨大，且许多参数无法通过实验来判断其取值，因此越来越多的模型参数率定工作都依赖于计算机算法。随着计算技术的发展，应用于复杂水环境模型参数识别的优化算法层出不穷。目前较为常见的参数率定方法有模拟退火算法（孟令群,郭建青.2009.利用混沌粒子群算法确定河流水质模型参数[J].地球科学与环境学报31(2):169-172）、遗传算法（杨斌斌,王文川.2010.多目标进化算法在新安江模型参数率定中的应用比较研究[J].水文.(03):38-42.）等确定性参数率定方法以及基于贝叶斯理论的GLUE方法（Beven,K,2012.Causalmodels as multiple working hypotheses about environmental processes.Comptes RendusGeoscience,344(2):77-88.Liu Y,Yang P J,Hu C,Guo H C.2008.Water quality modeling for loadreduction under uncertainty:A Bayesian approach[J].Water Research.42:3305–3314.Liu Y;FreerJ;Beven K,et al.2009.Towards a limits of acceptability approach to the calibration of hydrologicalmodels:Extending observation error.Journal of Hydrology,367(1–2):93-103.）。与其他确定性的参数率定方法不同，GLUE方法不将模型参数取一个确定值，而是将参数锁定在一定范围内，研究其分布情况，这种方式很好地体现了水环境复杂系统的不确定性特征。

不论采用何种算法率定模型参数，目前均是通过建立似然函数来对比模拟值与实测值，通过似然度判断模型参数是否有效。但当前水质模型复杂的非线性结构常常导致多组不同参数均能满足模型要求，从而无法从似然度来判断参数取值是否合理（Lindim C,Pinho JL,VieiraJ.M.P.2011.Analysis of spatial and temporal patterns in a large reservoir using water quality andhydrodynamic modeling[J].Ecological Modelling,222:2485-2494）。特别是对于包含一定误差或数据步长不足的的监测数据来说，即便似然度很高，也很难保证参数组合能够正确反映水质变化规律。

发明内容

为了克服现有方法中的单纯利用单一的似然函数进行水质模型率定时所出现的“异参同效”现象所导致的参数错误，本发明提供一种对多项输出指标同时率定来确定最优参数组合，进而实现水质模拟的方法，即通过综合考虑多项模拟指标，来进一步筛选满足条件的参数组合。该方法可以较好地避免因单个输出指标模拟而引起的参数偏差，能够更准确、有效地模拟水质指标。

本发明采用的技术方案如下：

一种基于GLUE多目标水质模型率定的水质模拟方法，其步骤包括：

1）根据全局敏感性分析方法对水质模型进行参数的敏感性分析，并根据得到的敏感性指数筛选出多目标敏感参数；

2）通过筛选出的多目标敏感参数的先验样本分布情况生成水质的随机参数组；

3）对所述随机参数组中的水质指标进行模拟，并通过似然函数进行筛选，通过迭代得到各项水质指标的敏感参数交集，最终使参数收敛至更小区间内，实现对水质的高精度模拟。

进一步地，步骤1）所述全局敏感性分析方法可以是Sobol方法、EFAST方法等。

进一步地，步骤3）设置信水平为0.5～1，进一步优选为0.75，用于判断单个参数值是否能够体现水质变化特征，并且通过对比不同水质指标的参数优化过程中每一次迭代的结果，来判断率定过程是否出现了过拟合现象。

进一步地，步骤3）每次筛选出模拟效果最好的样本，并依据最好的样本特征进行所述敏感参数交集的更新，进行若干次的迭代计算，最终使得敏感参数收敛至一个小区间内。

进一步地，步骤3）每次迭代生成样本数为5000个，每次根据似然函数值，筛选其中最优的15%的样本作为下一次迭代的样本。

进一步地，步骤3）所述似然函数优选采用NSE效率系数，也可以通过计算MSE均方根误差或者平均相对误差等方式来替代。

进一步地，步骤3）所述似然函数使用纳什系数（NSE）来衡量模拟值与观测值之间的拟合度，纳什系数表达式为：

$E_{\mathrm{ns}}=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(Q_{\mathrm{ob}}-Q_{\mathrm{sim}})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(Q_{\mathrm{ob}}-Q_{\mathrm{ob}\_\mathrm{average}})}^2},$

其中，Q_ob为观测值，Q_sim为模拟值，Q_{ob_average}为观测平均值，n为观测的次数。

优选地，所述敏感参数包括：硝化耗氧模块中的E12、K12、KNIT，沉积物耗氧模块中的SOD，水生植物生长模块中的GP1、PNH3，呼吸作用模块的E1R，沉淀模块的fD5，氮循环系统中的fon、anc。

进一步地，选取置信区间为（0.25，0.75）内的参数值构成交集的最大迭代次数内的交集。

与现有技术相比，本发明的有益效果如下：

1）该方法在率定的过程当中同时兼顾了多种水质指标的似然度，因此率定出来的结果可以同时满足多个水质指标的模拟效果，避免了因追求单个水质指标的似然度导致“过拟合”现象。从模型系统整体上提升了模拟的可信度；

2）该方法利用模型内部模块的耦合机制，通过多指标的模拟效果对模型参数取值进行验证，缩小参数的可行取值范围，从而能够大大减小模型参数所带来的模拟结果的不确定性。

3）由于该方法率定出来的参数是多项指标的模拟效果相互制约而得到的，因此可以过滤掉多种不满足条件的参数组合，为具有“异参同效”现象的复杂模型的参数率定工作提供一个更为可靠的方法。

采用本发明方法可以为后续的水环境水质管理、预测、评估等工作提供更为准确的定量化理论支撑。

附图说明

图1是本发明实施例的多目标率定算法流程图。

图2（a）—图2（d）是本发明实施例的水质模型的敏感参数图。

图3（a）—图3（j）是本发明实施例的依据不同指标进行率定的收敛过程对比图。

图4（a）—图4（d）是本发明实施例的水质模型率定效果图。

具体实施方式

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

本发明的多目标参数率定方法借鉴了GLUE方法的原理，同样是基于贝叶斯算法的思想（Zhang W T,George B.Arhonditsisa.2009.A Bayesian hierarchical framework for calibratingaquatic biogeochemical models[J].Ecological Modelling,(220):2142–2161），每次筛选出模拟效果最好的样本，之后依据最好的样本的特征进行参数集合的更新，进行若干次的迭代计算，最终使得参数收敛至一个小区间内，从而给出基于参数不确定性的响应范围。随着每一步迭代过程的进行，似然函数基本处于递增的状态（George B.Arhonditsis,Song S,Conrad L,et al.2007.Eutrophication risk assessment using Bayesian calibration of process-based models:Application to a mesotrophic lake[J].Ecological Modelling,(208):215–222）。这也就意味着，从总体上说，参数优化过程中迭代的后一步，会比上一步令参数取值更优。

但是实测值往往存在一定的误差，且似然函数的选择存在一定的人为因素，因此如果按照单一水质指标进行迭代时，可能存在某些特殊的不符合常理的参数组合使得水质指标的模拟效果很高，而造成“过拟合”的现象，显然这种参数组合是难以顾及其他水质指标的模拟效果的。对于WASP这类具有多种输出的模型来说，当分别对各个水质指标的模拟进行参数率定时，最终结果未必相同（Andrew M,Sincock,Howard S,et al.2003,Calibration andsensitivity analysis of a river water quality model under unsteady flow conditions[J].Journal ofHydrology,277:214–229）。在一个完整的水质模型体系中，应当兼顾各项指标的模拟效果，需要选取可以满足各项水质指标的模拟精度的参数组合，也即它们的交集部分。本发明的多目标参数率定方法正是基于此而进行改进的，该方法的流程如图1所示。

由于传统的GLUE方法是通过似然度来求一组参数的分布情况，当参数的维度较大时，则几乎不可能出现完全相同的参数组。因此完全通过不同水质指标率定出来的参数组交集来判定是否出现过拟合现象显然是不现实的。在对污染物组分以及水力条件较为固定的一段水域的模拟中，每一个参数的取值都应当介于一定的固定范围内，以体现区域的水质变化特征。因此参数往往只要满足在这个合理的区间内，基本上就能够体现水质变化特征。本实施例中设置信水平为0.5～1，优选为0.75，用于判断单个参数值是否能够体现水质变化特征，并且通过对比不同水质指标的参数优化过程中每一次迭代的结果，来判断率定过程是否出现了过拟合现象。在本实施例中，每次迭代生成样本数为5000个，每次根据似然函数值，筛选其中最优的15%的样本作为下一次迭代的样本。似然函数使用纳什系数（Nash-Suttcliffe系数）来衡量模拟值与观测值之间的拟合度，纳什系数表达式为：

$E_{\mathrm{ns}}=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(Q_{\mathrm{ob}}-Q_{\mathrm{sim}})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(Q_{\mathrm{ob}}-Q_{\mathrm{ob}\_\mathrm{average}})}^2},$

式中:Q_ob为观测值，Q_sim为模拟值，Q_{ob_average}为观测平均值，n为观测的次数。

通过不同指标的同时模拟，一方面可以避免单个水质指标率定可能引起的“过拟合”现象，同时也可以降低参数选取的不确定性。这是因为具有复杂结构的模型在模块之间常常具有关联作用，这样就使得单个参数的取值能够同时影响多个模块的计算过程，从而导致不同水质指标的模拟效果发生变化。这样就有利于形成一种相互验证参数的机制：当一个参数组合不能满足其他水质指标的模拟精度时，即便它可以让某一项水质指标的模拟效果很好，也会被视为“过拟合”现象而不被采纳。能够同时满足多个水质指标的参数组个数远远少于满足单个指标的参数组个数。如此就会大大减小其取值的范围，最终降低模拟的不确定性。

下面结合具体案例说明本发明方法的具体实施过程。

由于本发明基于GLUE方法进行分析，需要反复修改模型参数，并进行大次数的模拟。为了方便起见，利用Matlab环境下的Simulink工具，将WASP6.0模型按照模型说明中的内容进行了重写，这样可以通过撰写Matlab的运行脚本程序来实现对该模型的分析。

WASP模型包括EUTRO与TOXI两个主要模块，其中EUTRO模块用于模拟常规污染物，而TOXI模块用于模拟如重金属之类的有毒物质。DO、CBOD、硝态氮、氨氮这四个指标仅通过EUTRO模块就可以进行模拟，为了简化工作，本发明在Simulink中仅实现EUTRO中的必要模块。其中主要涉及的机理过程为：与DO有关的氧化、硝化、沉淀、水生植物生长、大气复氧、呼吸等过程；与CBOD有关的氧化、死亡、反硝化、沉淀等过程；与硝态氮有关的反硝化作用、硝化作用、水生植物生长等过程；与氨氮有关的藻类死亡、矿化、硝化作用、水生植物生长等过程。在Matlab环境中将各个子过程分别写入各自的m-function并定义其输入、输出变量，最终集成在.mdl文件中。

1）敏感参数筛选

利用Sobol方法对WASP模型进行参数的一阶敏感性分析及总敏感性分析，由于敏感性指数的精确度与采样数有很大关系，本实施例中尝试进行了多次运行，但结果差异不大。各参数的敏感性指数处于基本稳定的状态，其各项参数敏感性如图2（a）—图2（d）所示：其中，parameter为范围，ID为参数的编码，每个编码对于一个参数；图中S代表一阶敏感性系数，表示参数变动的直接影响；ST代表总敏感性系数，表示参数变动所产生的全部影响；两者之差ST-S代表参数之间的相互作用。

由图2（a）—图2（d）中的一阶敏感性指数与总敏感性的比例上可以看出：在DO与CBOD的模拟方面，参数的敏感性主要体现在单个参数变化的直接影响上，说明模型的结构耦合程度相对较弱；而硝态氮与氨氮的参数敏感性主要体现在参数之间的相互作用上，说明模型的耦合程度相对较高。对于耦合度较高的硝态氮与氨氮来说，参数率定难度更大，因为敏感参数的取值会对这两种指标同时产生显著影响；而DO与CBOD的参数取值相对较为独立，不会同时对这两个指标的取值同时产生过大的影响。

从图2（a）—图2（d）中敏感参数的分布情况来看，对于不同水质指标的模拟，其敏感参数不一样。将总敏感性指数大于0.1的参数作为需要进一步参与率定的敏感参数，包括：硝化耗氧模块中的E12、K12、KNIT，沉积物耗氧模块中的SOD，水生植物生长模块中的GP1、PNH3，呼吸作用模块的E1R，沉淀模块的fD5，氮循环系统中的fon、anc等共计10个敏感参数。这些敏感参数对于模型的模拟效果会产生极大的作用，而其他的参数则对模型模拟效果产生作用很小。因此为了简化计算，在率定过程中仅对这10个参数进行率定。

2）多目标参数率定

由于基于贝叶斯原理的率定结果提供的是参数的可行分布而非一个固定值，这就为参数的多目标寻优提供了可能。根据DO、CBOD、氨氮、硝态氮的模拟似然度分别进行迭代计算，每一步迭代所得到的参数取值分布如图3（a）—图3（j）所示。图中“1”（即灰色框）表示根据指标1即DO的率定迭代过程，“2”表示根据指标2即CBOD的率定迭代过程，“3”表示根据指标3即氨氮的率定迭代过程，“4”表示根据指标4即硝态氮的率定迭代过程。

随着迭代过程的进行，参数值范围一般来说是呈减小趋势的。按照某一项水质指标进行率定时，当参数范围缩减到一定程度，就不再与其他水质指标的率定参数范围形成交集，此时，模型系统对现实水环境的模拟的精度已经到达极限，如果再进行迭代，则会造成参数的“过拟合”。在本实施例中，根据图中所示的结果可以看出，各个参数会依据不同指标确定的目标函数得到的不一致的收敛区间。这说明如果按照单一指标进行率定，则会出现模型的“过拟合”，造成参数组无法适应该模型系统中其他水质指标的模拟。造成这样现象的原因一方面来自于模型机理的概化设计，另一方面则来自于实测值的误差，使得水质模拟不可能达到极高精度。作为一个体现水环境系统的动力模型，其参数率定必须从整体出发，令参数率定结果能够适应整个系统的模拟，必须兼顾各项指标的模拟效果，即便是单个水质指标的模拟结果似然度会因此而降低；若在率定过程中过度追求单一指标的模拟，则会造成对水质变化规律方面的误解。

如图3（a）—图3（j）所示，尽管最终的收敛区间不同，但可以从中发现一个规律：DO与CBOD的优化参数收敛结果较为接近，而氨氮与硝态氮的优化参数收敛结果较为接近。DO与CBOD之间的变化主要体现水体溶解氧被耗氧物质消耗的过程，氨氮与硝态氮之间的变化主要体现硝化与反硝化作用的动态过程，这两组过程应当同时在一个系统下进行模拟，而不是分别采取不同的参数组合来应对不同的水质指标的模拟。

在WASP模型中，CBOD模块与DO模块结合得较为紧密，而氨氮与硝态氮的结合更为紧密，结合紧密的模块之间存在更多的共用变量与共用参数。参数组合的变化对紧密关联的模块所产生的作用也会有更大的可能性产生共变。这使得寻优过程中，结合较为紧密的模块就可能存在模拟效果的同优或同劣。

本实施例在多个目标参数优选过程中，选取置信区间为（0.25，0.75）内的参数值构成迭代过程中所产生的最后交集。得到的参数范围如表1所示。

表1.参数范围

表1中参数的取值范围区间比较小，说明在该方法的率定下，模型的参数不确定性能够显著降低。需要指明的是，本发明实施例选择的参数置信区间具有一定的随意性，可以根据研究的需要进行调整。该范围内的参数取值可以满足各项模型输出的精度要求，比单目标率定方法要更加可靠。最终的模拟效果如图5（a）—图5（d）所示。

可以看出，WASP在模拟硝化作用与反硝化作用方面的效果要远好于模拟水体溶解氧于CBOD。对于CBOD与DO的模拟效果欠佳，除了模型结构与数据精度的原因，还可能是由于一些概化的常量参数的介入导致的。一些易受时间因素影响的参数在模拟时段内会发生变化，致使原来固定参数组不能够适应整体的模拟时段，通过误差累积，有可能会造成无法解释的现象（如本发明实施例中的溶解氧与CBOD的模拟）。对于这类问题，应尝试将数据分段，利用多段不同的常值组合来描述整个模拟时段内的水质变化过程，但这无疑会加大模型的数据需求。由此可见，当数据量不是非常充足的时候，仅通过单一目标函数来率定参数组，则难以反映整体系统的变化过程。

以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求所述为准。

Claims (10)

1.一种基于GLUE多目标水质模型率定的水质模拟方法，其步骤包括： 

1）根据全局敏感性分析方法对水质模型进行参数的敏感性分析，并根据得到的敏感性指数筛选出多目标敏感参数； 

2）通过筛选出的多目标敏感参数的先验样本分布情况生成水质的随机参数组； 

3）对所述随机参数组中的水质指标进行模拟，并通过似然函数进行筛选，通过迭代得到各项水质指标的敏感参数交集，最终使参数收敛至更小区间内，实现对水质的高精度模拟。 

2.如权利要求1所述的方法，其特征在于：步骤1）所述全局敏感性分析方法为下列方法中的一种：Sobol方法、EFAST方法。 

3.如权利要求1所述的方法，其特征在于，步骤1）所述多目标敏感参数包括：硝化耗氧模块中的E12、K12、KNIT，沉积物耗氧模块中的SOD，水生植物生长模块中的GP1、PNH3，呼吸作用模块的E1R，沉淀模块的fD5，以及氮循环系统中的fon、anc。 

4.如权利要求1所述的方法，其特征在于：步骤3）设置信水平为0.5～1，用于判断单个参数值是否能够体现水质变化特征，并且通过对比不同水质指标的参数优化过程中每一次迭代的结果，来判断率定过程是否出现了过拟合现象。 

5.如权利要求4所述的方法，其特征在于：步骤3）设置所述置信水平为0.75。 

6.如权利要求1所述的方法，其特征在于：步骤3）每次筛选出模拟效果最好的样本，并依据最好的样本特征进行所述敏感参数交集的更新，进行若干次的迭代计算，最终使得敏感参数收敛至一定区间内。 

7.如权利要求1所述的方法，其特征在于：步骤3）所述似然函数优选采用NSE效率系数，也可以通过计算MSE均方根误差或者平均相对误差等方式来替代。 

8.如权利要求1所述的方法，其特征在于：步骤3）所述似然函数使用纳什系数来衡量模拟值与观测值之间的拟合度，纳什系数的表达式为： 

其中，Q_ob为观测值，Q_sim为模拟值，Q_{ob_average}为观测平均值，n为观测的次数。 

9.如权利要求1所述的方法，其特征在于：步骤3）每次迭代生成样本数为5000个，每次根据似然函数值，筛选其中最优的15%的样本作为下一次迭代的样本。 

10.如权利要求1所述的方法，其特征在于：步骤3）选取置信区间为（0.25，0.75）内的参数值构成迭代过程中所产生的最后交集。