随着人们对水环境认知的深入，水环境模型的结构也变得日益复杂.常见的水质模型，如WASP、QUAL、QUASAR、MIKE(ECOLAB)、BASINS(RCHRES)等均是可同时模拟多项水质指标的复杂模型系统.尽管这些模型的框架逐渐完备，能够更全面地反映人类对水环境改变的认识，但同时也使得模型变量、参数不断增多，导致率定难度加大(潘登等，2012).而参数取值的准确与否，是模型能否正确反映水质变化规律的关键.

由于参数数量巨大，且许多参数无法通过实验来判断其取值，因此，越来越多的模型参数率定工作都依赖于计算机算法.随着计算技术的发展，应用于复杂水环境模型参数识别的优化算法层出不穷.目前，较为常见的参数率定方法有模拟退火算法(孟令群，2009)、遗传算法(杨斌斌等，2010)等确定性参数率定方法，以及基于贝叶斯理论的GLUE方法(Liu， 2008，2009; Beven，2012).与其他确定性的参数率定方法不同，GLUE方法不将模型参数取一个确定值，而是将参数锁定在一定范围内，研究其分布情况，这种方式很好地体现了水环境复杂系统的不确定性特征.

不论采用何种算法率定模型参数，目前均是通过建立似然函数来对比模拟值与实测值，通过似然度判断模型参数是否有效.但当前水质模型复杂的非线性结构常常导致多组不同参数均能满足模型要求，从而无法从似然度来判断参数取值是否合理(Lindim，2011).特别是对于包含一定误差或数据步长不足的监测数据来说，即便似然度很高，也很难保证参数组合能够正确反映水质变化规律.

为了降低错选参数的风险，本文以WASP(Water Quality Analysis Simulation Program)水质模型为例，提出一个对多项输出指标同时率定来确定最优参数组合的方法，即通过综合考虑多项模拟指标，来进一步筛选满足条件的参数组合.该方法可以较好地避免因仅仅评价单个输出指标的模拟效果而引起的参数偏差.

WASP模型是美国EPA推荐的水质模拟软件，在对河流(杨家宽等，2005;Lin et al., 2011)、湖泊(Arhonditsis et al., 2005; Chao et al., 2010;史铁锤等，2010)、河口(Arioli et al., 2005)等不同特征的水体的模拟预测中已有了大量的研究案例，被称为万能的水质模型.

为了对WASP模型进行大次数模拟并便于修改参数，本研究在MATLAB环境下，利用其内部的可视化仿真工具Simulink对WASP中的溶解氧(DO)、碳生化需氧量(CBOD)、氨氮、硝态氮等4个指标的模拟模块(EUTRO模块)进行分解，各个模块设立监视器，监控各个子模块的运算结果.在MATLAB环境下通过矩阵存储参数取值，按列依次读取参数组合，并完成模型蒙特卡洛模拟脚本的撰写工作.

为避免率定过程中因参数维度提高导致过大的计算量，需要对参数敏感性进行讨论.对比局部敏感性分析方法，全局敏感性分析方法的优势主要体现在：①可以获取参数之间的相互作用；②得到的结果更加接近于实际情况；③可以用于非线性的模型.由于全局敏感性分析方法的这些优势，可用于在各指标与脆弱性之间的响应关系中定量评估各项指标的直接影响与间接影响.目前，常见的全局敏感性分析方法包括Morris法、FAST法、Sobol法、Extend FAST法、GLUE法等.

由于复杂水环境模型的非线性结构，本文以通过Sobol敏感性分析方法确定的敏感参数作为本方法的研究对象.Sobol方法是目前最常见的全局敏感性分析方法之一(Nossent et al., 2011; Yang，2011)，其思想是利用参数变化所引起的方差来衡量参数的敏感程度，特别适合于非线性的模型参数敏感性分析.

借鉴GLUE方法的原理，多目标参数率定方法同样基于贝叶斯算法的思想(Zhang et al., 2009)，每次筛选出模拟效果最好的样本，之后依据最好的样本的特征进行参数集合的更新，进行若干次的迭代计算，最终使得参数收敛至一个小区间内，从而给出基于参数不确定性的响应范围.随着每一步迭代过程的进行，似然函数基本处于递增的状态(Arhonditsis et al., 2007; Missaghi et al., 2010).这也就意味着从总体上说，参数优化过程中迭代的后一步，会比上一步令参数取值更优.

但实测值往往存在一定的误差，且似然函数的选择存在一定的人为因素.因此，如果按照单一水质指标进行迭代时，可能存在某些特殊的不符合常理的参数组合使得水质指标的模拟效果很高，而造成“过拟合”的现象，显然这种参数组合是难以顾及其他水质指标的模拟效果的.对于WASP这类具有多种输出的模型来说，当分别对各个水质指标的模拟进行参数率定时，最终结果未必相同(Sincock，2003).在一个完整的水质模型体系中，应当兼顾各项指标的模拟效果，需要选取可以满足各项水质指标的模拟精度的参数组合，也即它们的交集部分.具体算法的流程如图 1所示.

图 1(Fig. 1)

图 1 多目标率定算法流程图 Fig. 1 Multi-object calibration algorithm flow chart

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

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

研究所用数据来自于课题组承担的国家重大水专项课题成果.研究区北运河为唯一发源于北京境内的水系，闸坝多、流速缓慢、天然水量少，多来自污水处理厂处理排放的再生水.课题组分别在2009年4月至10月(模型率定期)及2010年的4、7、10、12月(验证期)的非降雨时段对榆林庄至杨家洼闸段(图 2)的河道流量及水质(包括COD、CBOD、氨氮、硝态氮、DO、盐度、叶绿素、pH、水温等指标)进行了同步采样监测.

该河段地势平缓，流速缓慢，河段狭长，没有支流及取排水口对河道水量的干扰，且沿河水力状况变化较小，流经区域主要为农地、林草地，主要污染来自面源，在没有降雨的时段，基本可以视为没有外源汇入，边界条件相对简单.

图 2(Fig. 2)

图 2 研究案例河段示意图 Fig. 2 Schematic diagram of case study segment in North Canal

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

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

图 3(Fig. 3)

图 3 Simulink环境下的WASP模型结构 Fig. 3 WASP model built by Simulink in Matlab

表1(Table.1)

表1 模型所涉及的机理过程与参数 Table.1 Mechanism processes and parameters referred to in the simulation of the four indicators

表1 模型所涉及的机理过程与参数 Table.1 Mechanism processes and parameters referred to in the simulation of the four indicators 模拟指标 子过程 涉及参数1 涉及参数2 涉及参数3 编号 名称 意义 编号 名称 意义 编号 名称 意义 DO 氧化作用 10 KD 20℃耗氧系数 18 ED 氧化温度系数 12 KBOD 半饱和系数   硝化作用 1 K12 20℃硝化速率 2 E12 硝化温度系数 3 KNIT 半饱和系数   沉积物耗氧 4 SOD 沉积物耗氧速率 5 Es 沉积物温度系数   浮游植物生长 6 GP1 浮游植物生长速率 7 PNH3 选择项系数   大气复氧 9 E1R 温度系数   呼吸 8 K1R 20℃呼吸速率 9 E1R 呼吸温度系数 CBOD 氧化作用 22 KDS 有机碳降解速率 23 EDS 有机碳降解温度系数 12 KBOD 半饱和系数   死亡 20 aoc 氧碳比 21 K1D 死亡率   反硝化作用 17 K2D 20℃反硝化速率 18 E2D 反硝化温度系数 19 KNO3 半饱和系数   沉淀作用 16 fD5 溶解态CBOD占比 硝态氮 反硝化作用 17 K2D 20℃反硝化速率 2 E12 硝化温度系数 19 KNO3 半饱和系数   硝化作用 1 K12 20℃硝化速率 2 E12 硝化温度系数 3 KNIT 半饱和系数   浮游植物生长 6 GP1 浮游植物生长速率 30 anc 浮游植物氮氧比 7 PNH3 选择项系数 氨氮 死亡 28 DP1 藻类死亡速率 30 anc 浮游植物氮氧比 29 fon 有机氮的比例   矿化作用 33 K71 20℃有机氮矿化速率 31 E71 矿化温度系数 32 KmPc 磷限制因子半饱和常数   硝化作用 1 K12 20℃硝化速率 2 E12 硝化温度系数 3 KNIT 半饱和系数   浮游植物生长 6 GP1 浮游植物生长速率 30 anc 硝化温度系数 7 PNH3 选择项系数

利用Sobol方法对WASP模型进行参数的一阶敏感性分析及总敏感性分析，由于敏感性指数的精确度与采样数有很大关系，本研究尝试进行了多次运行，但结果差异不大.各参数的敏感性指数处于基本稳定的状态，其各项参数敏感性如图 4所示.

图 4(Fig. 4)

图 4 Sobol方法参数敏感性分析结果(图中S代表一阶敏感性系数，表示参数变动的直接影响；St代表总敏感性系数，表示参数变动所产生的全部影响；两者之差(St-S)代表参数之间的相互作用) Fig. 4 Sobol first order sensitivity and total sensitivity index of WASP parameter included in the simulation of the four indicators

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

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

如前所述，基于贝叶斯原理的率定结果提供的是参数的可行分布而非一个固定值，这就为参数的多目标寻优提供了可能.根据DO、CBOD、氨氮、硝态氮的模拟似然度分别进行迭代计算，每一步迭代所得到的参数取值分布如图 5所示，其中，DO的参数率定迭代过程通过绿色箱型图表示，CBOD通过红色箱型图表示，氨氮通过黑色箱型图表示，硝态氮通过蓝色箱型图表示.

图 5(Fig. 5)

图 5 依据不同指标进行率定的收敛过程对比(图中横坐标的iterations代表GLUE算法的迭代次数) Fig. 5 Variation of the updating parameter by GLUE method with different water quality indexes

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

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

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

本研究在多目标参数优选中，对比各个指标寻优过程中每一步迭代所产生的参数集合(本研究选取的置信区间为(0.25，0.75))，确定能够构成交集的最大迭代次数，选取该交集部分作为多目标确定的参数范围，得到结果见表 2.

表2(Table.2)

表2 多目标确定下的参数范围 Table.2 The parameter range confirmed by multi-objective optimization

表2 多目标确定下的参数范围 Table.2 The parameter range confirmed by multi-objective optimization 参数 取值 参数 取值 K12 0.103~0.104 PNH3 2.40~3.95 E12 1.11~1.13 E1R 1.021~1.025 KNIT 0.23~0.35 fD5 0.35~0.68 SOD 1.51~2.98 fon 0. 368~0.700 GP1 0.291~0.312 anc 0.193~0.415

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

图 6(Fig. 6)

图 6 模型率定效果 Fig. 6 Monte Carlo simulation result with reduced uncertainty

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

1)本文以WASP水质模型的应用为案例，提出了一套基于GLUE方法结合多目标验证的水质模型参数率定方法，研究结果表明，该方法在率定的过程当中同时兼顾了多种水质指标的似然度，因此，率定出来的结果可以同时满足多个水质指标的模拟效果，避免了因追求单个水质指标的似然度导致“过拟合”现象，从模型系统整体上提升了模拟的可信度.

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

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