南淝河城市河道水质模型的空间敏感性分析

谢如意，黄静水，周　琪

(同济大学环境科学与工程学院，上海　200092)

南淝河位于安徽省合肥市，是巢湖污染最严重的支流之一。合肥是安徽省的政治、经济和文化中心。2005年～2015年，合肥市人口增长了70%，达到了780万，国内生产总值增长了560%,达到了5 660亿元。作为城市主要的供水来源和收纳水体，城市的快速发展给南淝河带来了巨大的压力，水质修复成为需要解决的重要问题。

水质模型是进行河道水质研究的重要工具之一，广泛应用于情境管理和水质响应研究等方面，水质模型可以为水质管理提供决策支持。WASP、QUAL2K、EFDC、MIKE等都是现在常用的水质模型。但随着水质模型日益复杂，模型采用的参数也越来越多，而且大部分参数是无法精确测定的，为减少模型率定的困难，可以着重研究一部分对模型输出变化敏感的参数。

敏感性分析研究的是模型输入因子的变化如何影响模型输出的变化[1]。敏感性分析又分为局部敏感性和全局敏感性。全局敏感性对参数在整个取值空间的敏感性都进行分析，因此可以对参数的相互作用以及输入因子和模型输出的非线性关系进行研究。敏感性分析不仅为模型率定提供依据，也是对研究对象进行系统认识的有效工具之一。

在国内，敏感性分析方法在水质模型上已经有了一定的应用。张永祥等[2]将WASP模型应用于长河水体富营养化分析，采用的是单点扰动的局部敏感性分析方法，扰动幅度高达50%，未能很好地反映参数的真实敏感性。张质明等[3]以WASP水质模型的应用为例，通过Sobol方法确定模型的敏感参数，提出了一套基于GLUE法的多目标模型参数率定方法。Yi等[4]对建立的滇池水质模型的47个参数和7个外部扰动因子进行时间和空间的敏感性分析，阐明了敏感性分析对模型参数识别和不确定分析的重要性。本文的研究目的，是通过MATLAB的Simulink模拟仿真平台重建WASP水质模型的富营养化模块，结合SAFE工具包对多参数过程的复杂水质模型进行空间全局敏感性分析，探究水质模型参数的空间差异，进而揭示南淝河河道系统水质变化特征。

1　研究方法

1.1　WASP水质模型的选取

WASP是美国环保署推荐使用的水质模拟软件，是为分析池塘、湖泊、水库、河口和沿海水域的一系列水质问题而设计的多箱模型，可模拟常规污染物和有毒污染物在水中的迁移转化规律，被称为“万能水质模型”。

WASP的水质模拟分为两个模块，有毒物质模型TOXI和富营养化模型EUTRO。EUTRO模块用以模拟传统污染物的迁移转化规律，分为溶解氧、氮、磷和浮游植物四个子系统(图1)，包括8个状态变量和若干反应过程(表1)，具体的反应方程见WASP用户指导手册[5]。

图1　EUTRO模块各子系统及状态变量[5]

过程NH+4NO-3DIPPHYTCBODDODONDOP硝化作用√√√反硝化作用√√有机氮矿化√√浮游植物死亡√√√√√√浮游植物生长吸收√√√√√有机磷矿化√√氧化作用√√沉降√√√√大气复氧√底泥耗氧√内源呼吸√

1.2　Simulink环境下水质模型的构建

WASP自身并没有敏感性分析功能，只能借助其他的工具进行敏感性分析。同时，自WASP5后的版本不再提供源码，而且WASP本身是有界面的封装好的水质模拟软件，当需要多次批量修改参数运行调用时，操作十分不方便。最后，并没有专门针对WASP模型开发的敏感性分析工具，直接对WASP水质模型进行敏感性分析存在比较多的技术困难和操作不便性。

Simulink是MATLAB内部的可视化仿真工具，可用于实现动态系统的建模、仿真与分析。MATLAB/Simulink的主要优点有：(1)可视化的直观建模方式，简单明了，且各模块可以进行封装，大大增加了模型的可读性；(2)模型设计和模拟快速准确，还为用户提供了一个图形化的调试工具以辅助用户进行系统开发；(3)MATLAB/Simulink的仿真程序可以与MATLAB工具包实现良好的兼容，适于作为进一步研究和开发工艺控制程序的最佳平台。

为便于对模型进行敏感性分析，将WASP6.0模型按照模型手册中的内容在MATLAB/Simulink环境下重新建模。第一步，将EUTRO模块中4个子系统8个状态变量的所有子过程写入Simulink中，如图2a、图2b所示。第二步，建立对流项、弥散项和源汇项的方程(图2c)。因为本研究关注水质模型，所以水动力部分在EPDRiv1模型中算出后，直接在本模型中使用，具体的结果见Huang等的文章[6]。第三步，将所有方程封装为一个“块”，使用者只能对每个“块”的参数进行修改，而无法直接修改内部方程，这样可以提供更为简洁的使用界面。最后，根据南淝河的实际情况，将研究河段划分为45个小段，每段长约200～500 m，采用一个“块”来模拟，然后将45个“块”连接起来，建立南淝河水质模型(图2e)。水质模型构建完成后，为检查新模型的正确性，直接采用在WASP模型中手动率定的参数值，检验Simulink模型的正确性，部分参数值与WASP参数有区别。重建模型的模拟结果与WASP手动率定结果相关性如表2所示，各个指标的模拟结果相关性比较高，说明重构的南淝河水质模型结构上是准确的，可以在此基础上对模型敏感性进行分析。

图2　Simulink环境下建立南淝河水质模型

指标NH+4NO-3DIPChl⁃aDO模拟结果相关性098099074077099

选用Simulink工具对WASP水质模型进行重建，有以下有几个优点：第一，模型运行和敏感性分析工具包集合在MATLAB统一平台上，可以直接实现敏感性分析；第二，可以有针对性地对水质模型结构、参数设置等进行调整，比如WASP中大多数参数是集总式的，无法对同一参数在不同位置设置不同的值。本研究为对模型参数进行空间敏感性分析，在重建模型时，将模型参数改进为分布式，即不同河段的参数值可以不同。

表3　敏感性分析的参数及其取值范围

注：(a)表示参数范围来自Rates,Constants,and Kinetics Formulations in Surface Water Quality Modeling[7](第二版)一书；(b)表示参数范围来自《WASP6用户手册》[5]推荐值

SAFE(sensitivity analysis for everybody)是一款包含众多敏感性分析方法及可视化函数的MATLAB工具包[8]。其提供的敏感性分析方法包括：Morris法、RSA法、Sobol法、FAST法、动态识别分析法和基于密度的敏感性分析方法等。SAFE不仅提供了多种全局敏感性分析方法，还能够轻松地加入其他的敏感性分析方法。其提供的敏感性分析方法，都支持对敏感性指标进行鲁棒性和收敛性的评价。SAFE工具包中还提供了大量的可视化工具，可以对敏感性分析结果进行更好的研究和相互对比。这款工具包既能让非专业人员方便地操作，也允许有经验的使用者在此基础上开发新的功能。

本研究结合SAFE工具包中提供的LH-OAT采样方法及Morris敏感性分析方法，以望塘污水处理厂为分界点，将研究河段分为上下游，分别进行参数敏感性分析。整个模型参与敏感性分析的参数如表3所示，其参数取值范围来自于文献报道中的值。

1.3　Morris分析方法

Morris法[9]是一种定性的全局敏感性分析方法，主要用于参数的筛选和排序。Morris法属于一次只改变一个变量的方法，即选取模型某一变量xi，其余参数值固定不变，在变量取值变化范围内随机改变xi的值，运行模型得到目标函数值，用模型输出对模型输入的变化率EEs来表示参数变化对输出值的影响程度，多次变化后的平均变化率作为全局敏感性的指标如式(1)。

(1)

其中：Ci-比例因子;

g-模型输出;

M-参与敏感性分析的参数个数;

i-第i个输入参数;

j-第j次模型运算;

r-模型运行总次数;

Si-全局敏感性参数。

EEs的平均值(μ)代表参数对输出结果总体敏感性大小，EEs的标准差(σ)代表某个输入参数与其他参数相互作用的程度。高标准差代表某个参数的敏感性在其整个取值空间内都发生变化。Morris法的敏感性指数u和σ都只反映参数敏感性的相对大小，只能用于判断敏感性的相对大小，其值并不反映参数的真实敏感性[1]。

2　南淝河基本概况和模型建立

2.1　研究河段的基本概况

南淝河长约70 km，流域集水面积约1 446 km2，最终汇入巢湖。董铺水库至橡胶坝一段流经合肥城区，长约17 km，如图3所示，本研究模拟河段为图中红色河段，长约12 km。

图3　研究区域图(黑色河段为本研究模拟河段)

董铺水库下游约7 km处，河流两岸完全渠道化，为硬质河岸。城市河段内有两条支流：一条支流为四里河，其上游为大房郢水库；另一条支流为板桥河，沿岸主要为农业用地。董铺水库是合肥市主要的饮用水源地，拦截了南淝河上游的大部分清洁来水，导致南淝河城市河段基本没有清洁水来源。研究河段内有两个污水处理厂，望塘污水处理厂位于城市河段的上游，年均排放量约为2.21 m3/s；王小郢污水处理厂位于橡皮坝上游约1 km处，年均排放量约3.66 m3/s。污水处理厂尾水是南淝河城市河段水量的主要来源，旱季和雨季分别占河道流量的75%和53%。董铺水库下游3 km处，有一城中村，其未处理的生活污水直排进入南淝河。因此南淝河城区段是一段典型的以污水处理厂尾水为主要补给的城市化河段。

2.2　模型建立

3　空间敏感性分析结果与讨论

以望塘污水处理厂为分界点，分别对上下游河段进行敏感性分析，Morris分析结果如图4所示，横轴u*代表参数的总体敏感性，纵轴σ代表参数间相互作用大小。图4最后一组NSE_T的目标函数为五个指标的总体纳什系数，代表上下游模型整体输出对参数的敏感性情况。表4表示的是对各个指标及总体模拟结果中敏感性排名前10的参数排序。

图4　单指标及总体模拟效果纳什系数Morris敏感性分析结果图

排序NH+4NO-3DIPChl⁃aDO总体USDSUSDSUSDSUSDSUSDSUSDS1SODKmNcSODK2DDpK＇eDpK＇eK＇eSODK＇eK2D2DpK71K2DKNO3fDpfDpSODK＇eDpKNO33fKNITKNITSODK＇efK＇efDpK2DfDp4K＇eK12DpE2DK1caPCK1cK1caPCKNO3SODK＇e5KNITDpfK＇eIaKmPcIaIafEDK1cSOD6K1cK＇eK＇eESaPCk83vs4IsPNH3KDPNH3f7IaE12fD5EDIsfopaPCE1CIafIaK1c8aPCfONaNCKDvs4IaIsvs4K12ESaPCKmNc9PNH3PNH3KNO3fD5E1CIsE1CKmPK1cfD5vs4Ia10fD5aNCK12vs3fopK1cKmPEcIsK1cIsIs

注：US代表污水处理厂上游，DS代表污水处理厂下游

3.1　氨氮指标

3.2　硝酸盐指标

3.3　无机磷指标

污水处理厂上游，DIP浓度对浮游植物死亡速率最敏感(DP)，其次是与浮游植物生长相关的一系列参数。同时DIP在上游的参数敏感性排序与Chl-a的参数敏感性排序有很高的一致性，这主要是因为DIP与其他子系统的联系很少，只与浮游植物子系统发生密切联系，所以参数排序出现比较高的一致性，同样的结果在Yi等滇池水质模型的敏感性分析中也提到过[4]。

污水处理厂下游，DIP除了对与浮游植物相关的一系列参数敏感，还对与有机磷矿化的部分参数敏感，如K83、E83和KmPc等。由图4可知，有机磷矿化相关参数的σ值都比浮游植物相关参数的低，这也说明了浮游植物系统参数的相互作用非常大。

上下两段对比来看，磷系统自身的转化过程相对稳定，参数的相互作用程度低，主要是受浮游植物生长与死亡的影响，上下游均有这个特点。污水处理厂对磷系统的影响不大。

3.4　浮游植物指标

污水处理厂上游，Chl-a浓度对Dp最敏感，而且其参数间相互作用也最强。上游的敏感参数的u*和σ有很好的线性关系，与DIP的结果有一致性，这也说明了浮游植物相关参数的敏感性和相互作用可能存在内在的线性关系。

对比来看，上下游Chl-a浓度主要受浮游植物自身相关的一些参数的影响，说明在富营养化模块中，浮游植物系统是最核心最重要的子系统。其次，上下游敏感性排名前10的参数中都出现了KmP，有可能上下游浮游植物的营养限制都是磷限制。与浮游植物相关的各个过程：生长包括的三个限制项、死亡和沉降过程都在敏感性参数中有体现。这说明，浮游植物浓度是一个极其敏感的指标，与之相关的各个过程都直接地影响其浓度，更容易出现异参同效的情况，所以在参数识别和率定的时候要格外注意，模拟准确的难度也更大。

3.5　溶解氧指标

上下游对比来看，上游是以浮游植物光合作用产生氧气为主的自养型系统，下游已经成为耗氧过程占主导的系统。这个结果与Huang等直接通过水质模型计算DO平衡的结果相似[6]。这种情况下，由于污水处理厂尾水的大量排入以及沿岸污染物的进入，会使下游水体中的DO持续降低，甚至出现黑臭现象。

3.6　小结

下游最敏感的参数是K2D和KNO3，而且它们的σ值也排在前两位，说明反硝化作用对下游的模拟结果最敏感。除此之外，下段还主要受SOD和浮游植物生长死亡的影响，各个过程的影响在下游均有不同程度的体现。

以总体模拟的纳什系数为目标函数进行参数敏感性分析，发现污水处理厂上下游的参数敏感性发生了显著变化。上游是以浮游植物生长死亡为主导的自养型系统。下游由于河岸带变化、污水处理厂稀释等原因，浮游植物的浓度急剧减少，不再以浮游植物为主导。污水处理厂排放大量的硝酸盐进入河道，以及支流或者排口排入的大量有机物，导致下游河道溶解氧显著降低，反硝化作用强烈。这些变化反映了污水处理厂尾水的大量排入，打破了河流水质过程的连续性，在尾水排入的地方水质发生突变，河流上下游的主导过程发生了显著改变，一定程度可以揭示污水处理厂对河道系统带来的巨大影响。

4　结论

本研究选择在MATLAB/Simulink环境下重建南淝河水质模型，并结合SAFE工具包，实现对41个参数的空间敏感性分析发现参数对不同水质子系统敏感性存在差异，这些差异反映了水质系统的复杂性。同时，参数敏感性的空间差异，反映了水质过程的空间差异和研究对象的系统特征，揭示了污水处理厂尾水作为主要补给源对城市河道的巨大影响。污水处理厂上游是以浮游植物生长死亡为主的偏自然河道，河流保持在比较健康的状态；下游转变成以反硝化作用、硝化作用、SOD耗氧过程为主的受人为活动影响严重的河道，水质污染逐渐加重，甚至出现黑臭的情况。该方法很好地解决了WASP软件无法进行敏感性分析，以及无法设置分布式参数的弊端，是研究河道水质过程和系统特征的又一有效手段。

5　建议

根据分析结果，对南淝河治理提出整体性建议，针对上游河道，治理主要以截污为主，杜绝未处理生活污水直排的现象，增加河道上游清洁水来源，即可保持上游较好的水质状态。针对下游河道的治理，有以下几点建议：(1)进一步改进污水处理厂处理工艺，降低排放尾水中的硝氮和有机物负荷；(2)利用排放河道以及周边可以利用的条件尽可能减少尾水对河道水质的影响，比如尾水河道旁侧生态处理、生态生物处理模式等；(3)提升支流来水水质；(4)恢复下游河道的生态护岸，在护岸上种植水生植物和其他植物，增加下游河道的生物净化作用，同时营造柔美生态岸线，集防洪、生态、景观和自净等功能于一体。

[1]Pianosi F,Beven K,Freer J,etal.Sensitivity analysis of environmental models:A systematic review with practical workflow[J].Environmental Modelling & Software,2016(79):214-232.

[2]张永祥,王磊,姚伟涛,等.WASP模型参数率定与敏感性分析[J].水资源与水工程学报,2009,20(5):28-30.

[3]张质明,王晓燕,于洋,等.基于GLUE法的多指标水质模型参数率定方法[J].环境科学学报,2014,34(7):1853-1861.

[4]Yi X,Zou R,Guo H C.Global sensitivity analysis of a three-dimensional nutrients-algae dynamic model for a large shallow lake[J].Ecological Modelling,2016,327:74-84.

[5]Wool T A,Ambrose R B,Martin J L,etal.Water quality analysis simulation program (WASP) version 6.0 draft:User′s manual[Z].2013.

[6]Huang J,Yin H,Chapra S C,etal.Modelling dissolved oxygen depression in an urban river in China[J].Water,2017,9(7):520-538.

[7]Bowie G L,Mills W B,Porcella D B,etal.Rates constants and kinetics formulations in surface water quality modeling[Z].1985.

[8]Pianosi F,Sarrazin F,Wagener T.A Matlab toolbox for global sensitivity analysis[J].Environmental Modelling & Software,2015(70):80-85.

[9]Morris M D.Factorial sampling plans for preliminary computational experiments[J].Technometrics,1991,33(2):161-174.

[10]Griensven A Van ,Meixner T,Grunwald S,etal.A global sensitivity analysis tool for the parameters of multi-variable catchment models[J].Journal of Hydrology,2006,324(1-4):10-23.