一种改进的不确定性水质模型参数率定方法

张质明,王晓燕,潘润泽(1.首都师范大学资源环境与旅游学院,北京 100048；.北京建筑大学环境与能源工程学院,北京应对气候变化研究和人才培养基地,城市雨水系统与水环境教育部重点实验室,北京 100044)

一种改进的不确定性水质模型参数率定方法

张质明1,2,王晓燕1*,潘润泽2(1.首都师范大学资源环境与旅游学院,北京 100048；2.北京建筑大学环境与能源工程学院,北京应对气候变化研究和人才培养基地,城市雨水系统与水环境教育部重点实验室,北京 100044)

为了能够使得水质模型的参数率定更符合实际情况,减少水质模型当中常出现的异参同效的现象,提出一种能够结合实验与算法共同对模型参数进行识别的方法,即通过设置约束条件,将单纯的“与实测值进行对比”的一般率定过程转化为"能够对模型内部过程进行一定程度控制"的率定过程.以WASP模型在北运河的某河段中进行的水质模拟结果表明,通过引入约束条件可以从更为系统的角度有效地去除部分不符合污染物转化规律的参数组合,使得最终水质模型能够更为准确地反映实际的水质变化过程,同时也减少了模型的不确定性与模型参数率定的“异参同效”现象;在模型的非线性结构与约束条件的作用下,相同子模块中原本独立的参数会呈现一定的相关性;随着对区域的水质变化过程研究的深入,当引入新增约束条件或是降低容错比例时,约束条件的作用将会进一步增强.

水质模型；约束条件；参数率定；异参同效

当前水质模型成为水环境定量化管理中的重要工具,如WASP、QUAL2K、MIKE等复杂的机理模型在各类水体中得到了大量应用.为了使得机理模型能够更为详细地体现水体中污染物的迁移转化的过程,模型内部的结构复杂程度显著提升,参数的率定的难度也随之增加.

为了解决模型的参数率定问题,许多方法也随着计算机技术的发展而大量涌现,目前较为常见的若干参数率定方法包括:遗传算法、模拟退火算法、粒子群算法[1]、GLUE方法等以及各类方法所衍生出来的修正算法[2-4].

然而这类依靠似然函数进行率定的方法无法很好地解决复杂模型中的“异参同效”的问题.尽管机理模型具备理论概念的框架,但由于率定的似然函数仅仅针对观测值与最终模型的输出结果,忽视了对模型内部运算过程的验证,当多个参数组合的似然度近似,便无法判断参数取值的正确性.虽然这种现象可以通过大量的历史实测数据的验证来进行缓解,但由于水质监测的成本非常高[5-6],许多地方的监测投入不足或是由于历史数据缺失严重,很大程度上影响了模型在研究区域的可行性.

在模型的实际应用中,由于用于率定的实测数据本身具有不确定性,即使提供很长的时间序列以供模型的验证,也需要考虑这段时期是否具有代表性.当环境因外界原因(如污染源、水力条件等)产生较大程度变化时,用于描述机理的参数取值也不应是固定不变的.而如果忽视这个问题,进而通过机器算法不断调整参数,迫使模拟值不断逼近带有不确定性的实测值,就有可能会出现一些不合理的参数组合.实际上,水质模型参数取值的准确与否,不仅需要让最终模拟结果与实测值具有较高似然度,同时也应在模拟过程中体现水质变化的实际规律.

由于水环境模型参数的取值准确性非常关键,一些研究者直接采用实验方法对敏感关键参数进行测定,用于研究区的水质模拟[7-9];同时,利用实验室模拟的结果,也可以对模型的结构进行改进[10].当实验条件设计合理时,这类通过实验测定的参数的方法比通过计算机算法率定所得到的参数更准确,但这种针对特定参数的获取方法往往成本较高.此外,在具有复杂结构的模型中,能够直接测定的参数并不多.

为了在参数率定当中既能够较好兼顾实验室模拟的针对性与机器率定方法的效率,本文以WASP模型在北运河的某河段中模拟氮的转化为例,提出一种能够结合实验与算法共同对模型参数进行识别的方法.不同于直接测定特定参数,而是利用实验结果作为“约束条件”,既可以节约成本,也为引入相同研究区的其他研究结果提供了可能性,并且该方法避免了因“异参同效”所导致的问题,在很大程度上提高模型的率定的准确程度.

1 数据与方法

1.1 研究区域概况

选用北运河中一段水力条件与污染物汇入状况均较为简单的河段进行模拟,该河段总长约30km.数据采用课题组在2009年4月～10月(模型率定期)及2010年的4、7、10、12月(验证期)的非降雨时段对榆林庄至杨家洼闸段(图2)的河道流量及水质(包括COD、CBOD、氨氮、硝态氮、DO、盐度、叶绿素、pH、水温等指标)的同步监测成果.

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

图1 采样点空间分布与河道概化Fig.1 Location of monitoring sites of North Canal in Tongzhou

1.2 EUTRO模块

WASP模型中的EUTRO模块可以用来模拟常规污染物的迁移转化.该河段水质受氮污染的影响严重,因此本文采用其中的关于 CBOD、DO、氨氮、硝态氮等主要模块来进行该河段污染物降解的模拟.由于WASP模型的封装性,仅能观测到某项指标模拟的总体结果,因此本文基于MATLAB环境下利用Simulink对该模块进行了拆分,便于观测各部分的计算结果.

各项指标的计算过程如下,其中各项参数及变量的意义详见WASP模型的手册:

式中:第1项为浮游植物的死亡、第2项为CBOD的氧化、第3项为CBOD的沉淀、第4项为反硝化作用.

式中:第1项为大气复氧、第2项为因CBOD的氧化作用(碳化作用)、第3项为硝化作用、第4项为底泥耗氧、第5项为水生植物生长(光合作用)、第6项为水生植物的呼吸作用.

氨氮:

式中:第1项为水生植物死亡、第2项为矿化作用、第3项为水生植物生长、第4项氨氮硝化作用.

硝态氮:

式中:第1项为硝化作用、第2项为水生植物生长、第3项为反硝化作用.

1.2.1 率定方法 EUTRO当中涉及的水质指标转化过程比较多,不能通过实验方法一一进行测定.但由于参数所描述的水质变化过程应当遵循研究区域的水质变化的特征,与实验室模拟当中获取的结果不应差异过大.此外由于水环境作为一个开放系统,水体在自净过程中受到许多来自于人为、自然的外界干扰作用,因此在模型参数取值的选择上也应当考虑不确定性,大部分常值参数(除分子量比之类的具有确切值的参数以外)均应该考虑在一定范围内的变化所造成的影响.

在诸多率定方法中,GLUE算法将模型参数的分布状况替代了单一参数的取值[11-13],在模型不确定性的研究中得到了大量的应用[14-16].由于GLUE方法所借鉴的贝叶斯理论,是将前一次模拟得到的参数后验分布作为第 2次模拟时的先验分布.可以不断增加估计参数的信息,从而率定出最优的参数组合范围.每组随机生成的参数对应一个似然值,需要借助一个准则来判断是否接受这组参数.在本研究中,参数的范围根据先前的研究来确定,参数组合的抽样通过拉丁超立方的方式进行选取,经过模拟后,筛选出似然值最高的5%的参数组合,作为可被接受的参数集.

1.2.2 条件约束 由于北运河的氮污染物排放强度大,水体中的氮污染远远高于国外河流,因此其污染物的迁移转化也具有一定的特殊性,参数的本土化就显得非常重要.而在许多模型率定的研究中,参数的先验分布往往借鉴其他研究案例或者默认值来选取,其合理性有待商榷.为了能够对参数更好地进行本土化,本研究在率定之前,全部的参数需要经过约束条件进行初筛.在约束条件与模型两者的作用下,参数的初始分布会由均匀分布产生了改变,并且这些参数相互之间的独立性也可能因此而消失.

本研究利用课题组2008年5～8月在北运河沙河闸与杨洼闸采样,通过实验室模拟测定所获得的结果[17],作为条件来对模型参数的率定过程进行约束.这样可以保障参数取值的本地化.

如前所述,WASP模型当中的EUTRO模块将BOD的耗氧过程分为CBOD与NBOD两个部分,因此结合模型的结构特征,约束条件主要考虑以下3方面的内容.

(1) CBOD与NBOD比例

在采用罗威邦测定仪(德国),将水样加入到带有磁力搅拌器的棕色培养瓶中,在20℃条件下培养 10d,记录每天的值.通过对比加与不加丙烯基硫脲(ATU,用于抑制氨氧化细菌的活性)的 2组水样,其浓度差值就是NBOD,每组设置3个平行样,结果取其浓度的平均值.

根据测定结果,杨洼闸的CBOD与NBOD的平均耗氧量之比为18.25/23.79.然而,一方面该取值在自然环境下会受到外界因素干扰,应该只在一定程度上具有参考性;另一方面,作为约束条件,应当适当放宽范围,以便参数抽样的结果能够落在区间内,而非一个确切的点上.由于上游沙河闸的CBOD与NBOD的平均耗氧量之比为16.65/ 20.33,同时考虑到沿途排污口汇入时的干扰,变化范围约为 10%～15%,因此本研究设定相应容许范围为±15%,得到约束条件:

式中:CBOD为 BOD中碳化部分消耗的氧量, NBOD为硝化部分消耗的氧量.

有机物质的生物化学氧化反应一般分为两个阶段,第一阶段为碳氢化合物氧化为二氧化碳和水,称为碳化阶段;第二阶段氨被氧化为亚硝酸盐及硝酸盐,称为硝化阶段.在EUTRO模块当中,这两部分的耗氧分别对应溶氧模拟中的碳化与硝化两部分.其中,在碳化部分的计算,模型采用以 CBOD、溶解氧、温度的函数来实现,而关于硝化部分的耗氧量则通过氨氮、溶解氧、温度的函数来实现.

(2) 氨氮的降解速率

氨氮的降解系数可以通过模型模拟出来的氨氮变化速率来得到.在模型中,影响氨氮浓度的环节包括浮游植物的死亡/内源呼吸、矿化作用、浮游植物的生长、以及硝化作用4个环节的共同作用.作为氨氮变化其中的一个过程,硝化活性的衡量可通过 EUTRO中硝化过程对于氨氮衰减量的贡献来进行模拟.

根据对北运河氨氮降解系数的测定实验结果,该榆林庄桥至杨洼闸断面的氨氮降解系数为0.0089～0.1012/d.该实验是在自然条件下对水中氨氮的浓度进行连续监测得到的范围,因此在模型模拟时,该模块降解的量不应低于原来浓度的0.89%,且不应高于10.12%.根据北运河氨氮的浓度现状(北运河氨氮浓度为 10.1～23.3mg/L),限定其变化范围为:

其中:NH3_Degradation代表氨氮在单位时间内的降解量.

(3) 硝化活性

硝化过程是河流系统中氮循环的一个重要环节.通过硝化过程,氨氮转化为硝态氮与亚硝态氮,其速率反映了河流对氨氮降解过程的效率[18].经过对氨氧化过程和抑制氨氧化过程中硝化产物NO2

-和NO3-浓度变化情况可计算得到研究区水体的硝化活性范围为在0.24～4.29mg/(L⋅d).

在 EUTRO模型中,计算氨氮变化量的第四项为氨氮向的氮氧化物转化过程的模块.该模块的计算结果反映了模拟对象水体的氨氮浓度变化速度.

2 结果与讨论

2.1 可行参数约束结果

根据实验结果对参数集可行域范围的划定,可以发现一些参数的先验范围发生了变化.由于约束条件的存在,原本相互独立的参数取值呈现了一些相关性.例如图2中所示,参数ED的取值与参数 KD的取值就明显受到了约束条件的影响.这是由于ED与KD 2个参数同属氧化过程的模拟环节,该环节对本研究涉及到的约束条件都有显著的机理关系,因此约束的结果较为明显.在参数抽样环节中,对于存在强相关性的参数,不应逐一取值,而是应当考虑在筛选的参数组集合当中进行选取,以免出现错误的组合.

事实上,在水质模型中,污染物迁移转化规律的描述常常需要通过多个参数才能够完成.因此参数组合的确定十分重要,如果对先验参数不进行筛选,而直接进行独立取值,许多不合理的参数组合将可能出现,显然这些参数组合会歪曲真实的模拟过程.而这种错误在模型的应用中有时很难被发现:一方面,用于评估模型效果的似然度往往通过人为设定的似然函数进行计算而得到,而复杂模型当中“过拟合”或“同效”的现象可能会使得一些似然度看起来很高,但事实上是由不合理的参数组合所形成的“最优解”;另一方面,似然函数的计算往往只针对最终的模拟结果与实测值的对比,而非基于过程进行监控,对于里面子模块的响应正确与否无从验证.例如,本研究当中的氨氮模拟分为四项:水生植物死亡、矿化作用、水生植物生长、氨氮硝化,这4项的增减都会影响到最终氨氮浓度的大小.显然,如果不对子模块的响应进行验证,每个子模块响应的不确定性则会较大,在“同效”现象存在时,就很难识别出氨氮增减的真正原因.

当将约束条件介入参数的选取时,不符合条件的参数组合会被舍弃.因此,即便这部分参数当中某些参数组合的似然度较高,也仍然不会被选择,从而减少模拟错误的可能性.随着约束条件的增加,可被筛选的参数组合的数量也会降低,因“同效”现象导致的错误也会随之减少.

在本研究中,通过约束条件来确定先验分布,会改变初始参数组合,在筛除一部分参数组合的同时,也可以在一定程度上减少因"异参同效"而引起的对污染物转化机理模拟的偏差.除了对参数初始范围的影响外,约束条件可能还会在更高维度的空间对参数的分布进行改变,但这难以通过图示对比看到,为了简单示意,图中仅展示了部分参数的二维分布,可以看到通过依据实测结果的筛选,这些参数的可行范围以及参数的分布情况发生了一定的变化.如图所示:落在浅色点覆盖范围内的任意一组组合都满足参数约束条件中所归定的规则,可视为符合北运河实测结果;而落在浅色点覆盖区域之外的部分,尽管其符合单个参数的选取范围(黑色区域),但与北运河的研究结果反映出来的水质转化规律相矛盾,将会在之后的率定过程中舍弃.在本研究中,经过约束条件筛选出来的参数大约占原样本的 1/5,这与本研究所采用的容错范围,约束条件的种类都有很大关系.毫无疑问,随着约束条件的增加,满足约束条件的参数样本会进一步减少.

2.2 模型率定结果

将符合约束条件的参数组合作为先验参数组,利用GLUE方法,对模型的参数取值进行率定.根据 WASP模型的敏感性分析研究筛选出来的10个最为敏感的参数[19]的率定值如图 3所示.

图3 GLUE算法对模型主要敏感参数的筛选结果Fig.3 Calibration of main sensitive parameters combination by GLUE algorithm

可以看出,经过模型参数的约束,一些模型参数在似然度较高的部分出现了较为明显的集中.尽管有一些参数仍然不能很好地通过 GLUE方法进行识别,这是由于实测数据数量较少或者模型本身概化所致,说明在这些参数当中"异参同效"的现象比较明显.

一般来说,GLUE方法的一个重要思想是强调参数组合对模拟结果的影响.但是目前 GLUE方法也是基于对模拟结果与实测值之间的似然度来进行参数组合的筛选,当实测数据不足,或是其中一部分数据缺乏代表性,就会在很大程度上干扰到最终所选取的参数组合.而本研究将部分已知的研究区的水质转化规律作为约束条件,把原本主观的先验分布转变为经过筛选的参数组合,从另一个角度确保参数组合在描述水质转化过程中的可靠性.这样率定出来的结果既可以保证对最终结果的拟合效果,也能够同时兼顾对水质转化过程的控制.

3 结论

3.1 传统的通过抽样方法获取的参数组合往往存在大量不符合污染物转化规律的情况,通过引入约束条件可以从更为系统的角度有效地去除部分不符合污染物转化规律的参数组合,使得最终水质模型能够更为准确地反映实际的水质变化过程,同时也因此减少了模型的不确定性与模型参数率定的"异参同效"现象;

3.2 本研究总结了CBOD与NBOD比例、氨氮的降解速率与硝化活性等三个在研究区河段进行研究的成果,作为筛选合理参数的约束条件.在模型的非线性结构与约束条件的作用下,相同子模块中原本独立的参数会呈现一定的相关性;

3.3 在本研究中,符合约束条件的参数集合约为原样本量的 1/5.随着对区域的水质变化过程研究的深入,当引入新增约束条件或是降低容错比例时,这一比例还会进一步降低.

[1] 刘苏宁,甘 泓,魏国孝.粒子群算法在新安江模型参数率定中的应用 [J]. 水利学报, 2010,(5):537-544.

[2] 郭建青,李 彦,王洪胜,等.应用单纯形-模拟退火混合算法估计河流水质参数 [J]. 水科学进展, 2004,(6):765-769.

[3] 王纲胜,夏 军,陈军锋.模型多参数灵敏度与不确定性分析 [J].地理研究, 2010,(2):263-270.

[4] 刘 毅,陈吉宁,杜鹏飞.环境模型参数优化方法的比较 [J]. 环境科学, 2002,(2):1-6.

[5] Mannina G, Viviani G. Water quality modelling for ephemeral rivers: Model development and parameter assessment [J]. Journal of Hydrology, 2010,393(3/4):186-196.

[6] Naddeo V. River water quality assessment: Implementation of non-parametric tests for sampling frequency optimization [J]. Land Use Policy, 2013,30(1):197-205.

[7] Raimonet M, L Vilmin, N Flipo, et al. Modelling the fate of nitrite in an urbanized river using experimentally obtained nitrifier growth parameters [J]. Water Research, 2015,73:373–387.

[8] 逄 勇,丁 玲,高 光.基于生态槽实验的藻类生长参数确定[J]. 环境科学, 2005,(3):78-82.

[9] 欧阳威,蔡冠清,黄浩波,等.基于原位土壤观测的 SWAT关键参数及模拟优化分析 [J]. 农业环境科学学报, 2014,(8):1601-1608.

[10] 王飞儿,杨 佳,李亚男,等.基于沉积物磷释放的 WASP水质模型改进研究 [J]. 环境科学学报, 2013,(12):3301-3308.

[11] Arhonditsis G B, Qian S S, Stow C A, et al. Eutrophication risk assessment using Bayesian calibration of process-based models: Application to a mesotrophic lake [J]. Ecological Modelling, 2007,208(2–4):215-229.

[12] Missaghi S, Hondzo M. Evaluation and application of a three-dimensional water quality model in a shallow lake with complex morphometry [J]. Ecological Modelling, 2010,221(11): 1512–1525.

[13] Zhang W, Arhonditsis G B. A Bayesian hierarchical framework for calibrating aquatic biogeochemical models [J]. Ecological Modelling, 2009,220(220):2142-2161.

[14] Beven K, Freer J. Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology [J]. Journal of Hydrology, 2001,249(1–4):11-29.

[15] Liu Y, P Yang, C Hu, et al. Water quality modeling for load reduction under uncertainty: A Bayesian approach [J]. Water Research, 2008,42(13):3305-3314.

[16] Ratto M, Tarantola S, Saltelli A. Sensitivity analysis in model calibration: GSA-GLUE approach [J]. Computer Physics Communications, 2001,136(3):212-224.

[17] 于 洋.北运河水体中氨氮的氧化过程及微生物响应特征 [D]. 2012.首都师范大学.

[18] 王 超,单保庆,赵 钰.滏阳河水系沉积物硝化速率分布及溶解氧的限制效应 [J]. 环境科学学报, 2015,35(6):1735-1740.

[19] 张质明,王晓燕,李明涛.基于全局敏感性分析方法的 WASP模型不确定性分析 [J]. 中国环境科学, 2014,34(5):1336-1346.

A parameter calibration method with constraint based on laboratory experimental result.

ZHANG Zhi-ming1,2,WANG Xiao-yan1*, PAN Run-ze2(1.College of Resources, Environment and Tourism, Capital Normal University, Beijing 100048, China；2.Beijing Climate Change Response Research and Education Center, Beijing University of Civil Engineering and Architecture, Beijing 100044, China). China Environmental Science, 2017,37(3)：956～962

In order to make the parameters of the water quality model more realistic and reduce the equivalent phenomenon of different parameters in water quality model, this paper proposed a method to identify the model parameters by setting laboratory experimental constraints, which could be used to make a certain degree of internal process control. The example of a water quality simulation by WASP model in a section of the North Canal showed that by this method, the model could be more accurate to reflect the actual water quality change process by reducing uncertainty and equivalent from combination of the parameters. Through the action of the nonlinear model structure and constraint conditions, the originally independent parameters under the same sub module began to demonstrate certain correlation. With the further study of the water quality change process, the constraint conditions would be further enhanced when introduced new constraint conditions or fault decrease tolerance.

water quality model；constraint；parameter calibration；equivalent

X703

A

1000-6923(2017)03-0956-07

张质明(1984-),北京人,讲师,博士,主要研究方向为水环境模拟.发表论文10余篇.

2016-07-01

北京市自然科学基金青年项目(8154044);国家自然科学基金项目(41271495);高等学校博士学科点专项科研基金联合资助项目(20121108110006)

* 责任作者, 教授, cnuwxy@sohu.com