摘要：为认识南淝河的水质特性，本文在MATLAB/Simulink环境下构建南淝河水质模型，通过敏感性分析对城市河道的水质过程特征进行深入探究。以污水处理厂为分界点，用Morris敏感性分析方法对总共41个参数，进行上下游的敏感性分析。结果发现，由于污水处理厂大量尾水的排入，使南淝河水质系统从上游以浮游植物生长死亡为主导过程的偏自然系统，转换为下游以反硝化、硝化作用和SOD耗氧等过程为主导的受人类影响严重的系统。该结论为认识河道水质特性和选择恰当的水质修复方法提供理论依据，也为水质模型敏感性分析提供了新的分析思路，同时证明了敏感性分析方法是研究河流系统特征的有效方法之一。
南淝河位于安徽省合肥市，是巢湖污染最严重的支流之一。合肥是安徽省的政治、经济和文化中心。从2005年到2015年，合肥市人口增长了70%，达到了780万，国内生产总值增长了560%达到了5 660亿元。作为城市主要的供水来源和收纳水体，城市的快速发展给南淝河带来了巨大的压力，水质修复成为需要解决的重要问题。
水质模型是进行河道水质研究的重要工具之一，广泛应用于情境管理和水质响应研究等方面，水质模型可以为水质管理提供决策支持。WASP、QUAL2K、EFDC、MIKE等都是现在常用的水质模型。但随着水质模型日益复杂，模型采用的参数也越来越多，而且大部分参数是无法精确测定的，为减少模型率定的困难，可以着重研究一部分对模型输出变化敏感的参数。
敏感性分析研究的是模型输入因子的变化如何影响模型输出的变化。敏感性分析又分为局部敏感性和全局敏感性。全局敏感性对参数在整个取值空间的敏感性都进行分析，因此可以对参数的相互作用以及输入因子和模型输出的非线性关系进行研究。敏感性分析不仅为模型率定提供依据，也是对研究对象进行系统认识的有效工具之一。
在国内，敏感性分析方法在水质模型上已经有了一定的应用。张永祥等将WASP模型应用于长河水体富营养化分析，采用的是单点扰动的局部敏感性分析方法，扰动幅度高达50%，未能很好的反应参数的真实敏感性。张质明等以WASP水质模型的应用为例，通过Sobol方法确定模型的敏感参数，提出了一套基于GLUE法的多目标模型参数率定方法。Yi Xuan等对建立的滇池水质模型的47个参数和7个外部扰动因子进行时间和空间的敏感性分析，阐明了敏感性分析对模型参数识别和不确定分析的重要性。本文的研究目的，是通过MATLAB的Simulink模拟仿真平台重建WASP水质模型的富营养化模块，结合SAFE工具包对多参数过程的复杂水质模型进行空间全局敏感性分析，探究水质模型参数的空间差异，进而揭示南淝河河道系统水质变化特征。
1 研究方法
1.1 WASP水质模型的选取
WASP是美国环保署推荐使用的水质模拟软件，是为分析池塘、湖泊、水库、河口和沿海水域的一系列水质问题而设计的多箱模型，可模拟常规污染物和有毒污染物在水中的迁移转化规律，被称为“万能水质模型”。
WASP的水质模拟分为两个模块，有毒物质模型TOXI和富营养化模型EUTRO。EUTRO模块用以模拟传统污染物的迁移转化规律，分为溶解氧、氮、磷和浮游植物四个子系统(如图1所示)，包括8个状态变量和若干反应过程(如表2所示)，具体的反应方程见WASP用户指导手册。
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中，如图2-a，2-b。第二步，建立对流项、弥散项和源汇项的方程(图2-c)。因为本文关注于水质模型，所以水动力部分在EPDRiv1模型中算出后，直接在本模型中使用，具体的结果见黄静水的文章。第三步，将所有方程封装为一个“块”，使用者只能对每个“块”的参数进行修改，而无法直接修改内部方程，这样可以提供更为简洁的使用界面。最后，根据南淝河的实际情况，将研究河段划分为45个小段，每段长约200-500m，采用一个“块”来模拟，然后将45个“块”连接起来，建立南淝河水质模型(如图2-e)。水质模型构建完成后，为检查新模型的正确性，直接采用在WASP模型中手动率定的参数值，检验Simulink模型的正确性，部分参数值与WASP参数有区别。重建模型的模拟结果与WASP手动率定结果相关性见表2，各个指标的模拟结果相关性比较高，说明重构的南淝河水质模型结构上是准确的，可以在此基础上对模型敏感性进行分析。
图2 Simulink环境下建立南淝河水质模型
选用Simulink工具对WASP水质模型进行重建，有以下有几个优点：第一，模型运行和敏感性分析工具包集合在MATLAB统一平台上，可以直接实现敏感性分析;第二，可以有针对性地对水质模型结构、参数设置等进行调整，比如WASP中大多数参数是集总式的，无法对同一参数在不同位置设置不同的值。本研究为对模型参数进行空间敏感性分析，在重建模型时，将模型参数改进为分布式，即不同河段的参数值可以不同。
表3 敏感性分析的参数及其取值范围
SAFE(Sensitivity Analysis for Everybody)是一款包含众多敏感性分析方法及可视化函数的MATLAB工具包。其提供的敏感性分析方法包括：Morris法、RSA法、Sobol法、FAST法、动态识别分析法和基于密度的敏感性分析方法等。SAFE不仅提供了多种全局敏感性分析方法，还能够轻松的加入其它的敏感性分析方法。其提供的敏感性分析方法，都支持对敏感性指标进行鲁棒性和收敛性的评价。SAFE工具包中还提供了大量的可视化工具，可以对敏感性分析结果进行更好的研究和相互对比。这款工具包既能让非专业人员方便地操作，也允许有经验的使用者在此基础上开发新的功能。
本研究结合SAFE工具包中提供的LH-OAT采样方法及Morris敏感性分析方法，以望塘污水处理厂为分界点，将研究河段分为上下游，分别进行参数敏感性分析。整个模型参与敏感性分析的参数如表3所示，其参数取值范围来自于文献报道中的值。
1.3 Morris分析方法
Morris法是一种定性的全局敏感性分析方法，主要用于参数的筛选和排序。Morris法属于一次只改变一个变量的方法，即选取模型某一变量xi，其余参数值固定不变，在变量取值变化范围内随机改变xi的值，运行模型得到目标函数值，用模型输出对模型输入的变化率EEs来表示参数变化对输出值的影响程度，多次变化后的平均变化率作为全局敏感性的指标，公式如下：
其中：为比例因子，g代表模型输出，M为参与敏感性分析的参数个数，i代表第i个输入参数，j代表第j次模型运算，r为模型运行总次数，Si为全局敏感性参数。
EEs的平均值(μ)代表参数对输出结果总体敏感性大小，EEs的标准差(σ)代表某个输入参数与其他参数相互作用的程度。高标准差代表某个参数的敏感性在其整个取值空间内都发生变化。Morris法的敏感性指数u和σ都只反映参数敏感性的相对大小，只能用于判断敏感性的相对大小，其值并不反映参数的真实敏感性。
为了更好的实现采样，本研究采用LH-OAT采样方法，模型运行次数为r×(M+1)，极大的减少了模型的运算成本，适用于结构复杂、参数众多，对参数进行筛选和排序的情形。具体原理见A. van Griensven 。本研究的分析参数为M = 41个，采样大小选择r = 80，目标函数选择纳什系数(NSE)，分别考虑NH4+、NO3-、DIP、Chla、DO单个水质指标模拟效果的纳什系数，以及五个指标一起的模拟效果的纳什系数。
2 南淝河基本概况和模型建立
2.1 研究河段的基本概况
南淝河长约70 km，流域集水面积约1 446 km2，最终汇入巢湖。从董铺水库至橡胶坝一段流经合肥城区，长约17 km，如图3所示，本研究模拟河段为图中红色河段，长约12 km。自董铺水库下游约7 km处，河流两岸完全渠道化，为硬质河岸。城市河段内有两条支流，一条支流为四里河，其上游为大房郢水库;另一条支流为板桥河，沿岸主要为农业用地。董铺水库是合肥市主要的饮用水源地，拦截了南淝河上游的大部分清洁来水，导致南淝河城市河段基本没有清洁水来源。研究河段内有两个污水处理厂，望塘污水处理厂位于城市河段的上游，年均排放量约为2.21 m³/s;王小郢污水处理厂位于橡皮坝上游约1 km处，年均排放量约3.66 m³/s。污水处理厂尾水是南淝河城市河段水量的主要来源，旱季和雨季分别占河道流量的75%和53%。董铺水库下游3 km处，有一城中村，其未处理的生活污水直排进入南淝河。因此南淝河城区段是一段典型的以污水处理厂尾水为主要补给的城市化河段。
2.2 模型建立
2015年10月对南淝河城市河段14个采样点(R1-R14))进行采样，每个采样点的间距约为500~1 000 m，高空间密度的采样点布设可以更好地获得水质变化的空间信息。监测指标包括温度(T)、溶解氧(DO)、酸碱度(pH)、浊度、电导率(EC)、总溶解性磷(TDP)、硝氮(NO3-)、氨氮(NH4+)、总溶解性氮(TDN)、总氮(TN)、溶解性有机碳(DOC)等。具体的采样、实验过程、水质结果及模型建立，可参见课题组之前的研究。本文是在其基础上，利用MATLAB/Simulink平台和SAFE工具包，对已建立的水质模型进行参数敏感性的空间分析。
3 空间敏感性分析结果与讨论
以望塘污水处理厂为分界点，分别对上下游河段进行敏感性分析，Morris分析结果如图4所示，横轴u*代表参数的总体敏感性，纵轴σ代表参数间相互作用大小。图4最后一组NSE_T的目标函数为五个指标的总体纳什系数，代表上下游模型整体输出对参数的敏感性情况。表4表示的是对各个指标及总体模拟结果，敏感性排名前十的参数排序。
3.1 氨氮指标
污水处理厂上游，对NH4+来说，最敏感的四个参数是SOD、DP、f和K’e，除SOD外，都是和浮游植物系统相关的参数。同时，DP、SOD、Ia和f的σ值较高，说明这些参数与其他参数的相互作用比较强，属于高度敏感非线性的参数，比较难以识别。在敏感性排序前10位中，SOD反应底泥耗氧过程，KNIT反应的是硝化过程，f、DP、K’e、K1C、Ia、PNH3、aPC代表的是与浮游植物相关的过程，fD5是与CBOD沉降相关的参数。
污水处理厂下游，对NH4+来说，KmNc的u*和σ都远远高于其他参数，KmNc是氮循环过程中浮游植物限制的半饱和常数，反应的是有机氮矿化生成氨氮的过程。K12、KNIT、E12反应硝化过程，KmNc、K71反应的是有机氮矿化过程，PNH3、K’e和aNC反应的是浮游植物生长过程，DP、fon反应的浮游植物死亡过程。
污水处理厂上下段对比来看，对NH4+来说，上段是浮游植物相关参数占主导的系统，虽然硝化过程、底泥耗氧和CBOD沉降对NH4+浓度都有一定影响，但远不及浮游植物的影响。而且出现的SOD、fD5等参数都不是直接与氨氮反应过程相关的参数，很可能说明这些参数是通过影响浮游植物子系统和DO子系统，来间接对上游的氨氮浓度产生影响的。与上游不同，下游是多过程共同作用于NH4+变化的体系，与氨氮相关的四个直接过程的参数都在参数排序前十名中有体现。
3.2 硝酸盐指标
污水处理厂上游，NO3-对SOD最敏感，u和σ值都高于其他参数。说明上段的SOD耗氧影响了上游的DO浓度，进而影响水体中的NO3-浓度。K2D、KNO3代表反硝化作用，KNIT、K12代表硝化作用，f、DP、K’e、aNC代表浮游植物相关的过程，fD5与CBOD沉降有关，可能与SOD一样，通过影响DO进而对NO3-浓度间接产生影响。
污水处理厂下游，对NO3-来说，KNO3和K2D的敏感性和相互作用都远远高于其他参数，说明污水处理厂下游的NO3-主要受反硝化作用的影响，其他过程参数的影响基本可以忽略不计。
上下游对比来看，上游NO3-浓度主要受SOD的控制，其他过程共同发生作用。污水处理厂尾水排入河道，给河道带来了大量的NO3-负荷，直接改变了河道水质的连续性，下游的NO3-浓度被反硝化作用主导。上下游的主导过程已经发生了显著地变化，在进行水质模拟的时候要格外注意。尽管下游河道反硝化能力大幅度提升，但若污水处理厂排放负荷过大，超过河道的氮滞留能力，最终还是会对汇入湖泊的水质造成影响。
3.3 无机磷指标
污水处理厂上游，DIP浓度对浮游植物死亡速率最敏感(DP)，其次是与浮游植物生长相关的一系列参数。同时DIP在上游的参数敏感性排序与Chla的参数敏感性排序有很高的一致性，这主要是因为DIP与其他子系统的联系很少，只与浮游植物子系统发生密切联系，所以参数排序出现比较高的一致性，同样的结果在Yi Xuan对滇池水质模型的敏感性分析中也提到过。
污水处理厂下游，DIP除了对与浮游植物相关的一系列参数敏感，还对与有机磷矿化的部分参数敏感，如K83、E83和KmPc等。由图4可知，有机磷矿化的相关参数的σ值都比浮游植物相关参数的低，这也说明了浮游植物系统参数的相互作用非常大。
上下两段对比来看，磷系统自身的转化过程相对稳定，参数的相互作用程度低，主要是受浮游植物的生长与死亡的影响，上下游均有这个特点。污水处理厂对磷系统的影响不大。
3.4 浮游植物指标
污水处理厂上游，Chla浓度对DP最敏感，而且其参数间相互作用也最强。上游的敏感参数的u*和σ有很好的线性关系，与DIP的结果有一致性，这也说明了浮游植物相关参数的敏感性和相互作用可能存在内在的线性关系。
污水处理厂下游，Chla浓度仍是对K’e和DP最敏感，但K’e的交互作用远大于DP，光限制很可能是下游浮游植物的主要限制因素。
对比来看，上下游Chla浓度主要受浮游植物自身相关的一些参数的影响，说明在富营养化模块中，浮游植物系统是最核心最重要的子系统。其次，上下游敏感性排名前十的参数中都出现了KmP，有可能上下游浮游植物的营养限制都是磷限制。与浮游植物相关的各个过程：生长包括的三个限制项、死亡和沉降过程都在敏感性参数中有体现。这说明，浮游植物浓度是一个极其敏感的指标，与之相关的各个过程都直接地影响其浓度，更容易出现异参同效的情况，所以在参数识别和率定的时候要格外注意，模拟准确的难度也更大。
3.5 溶解氧指标
污水处理厂上游，对DO来说，敏感性排名前10的参数都有比较高的敏感性和相互作用。对DO影响最大的，依次是K’e、SOD、DP、K’e、K1C、PNH3、aPC、f、Ia、Is都是与浮游植物相关的参数。只有SOD代表底泥耗氧，K12代表硝化作用耗氧。说明上游的溶解氧平衡，主要受浮游植物的影响，但硝化作用和底泥耗氧都有参与其中。与其他参数相比而言，CBOD耗氧过程相关参数的敏感性偏低。
污水处理厂下游，敏感性参数已经与上段发生了很大变化。SOD的敏感性最大，K’e的交互作用最强，其他参数的敏感性和交互作用差不多。SOD、ES代表底泥耗氧; K2D、KNO3代表反硝化作用;KD、ED、fD5代表了CBOD耗氧;只有K1C、f和K’e代表浮游植物生长。由此可见，下游DO是多个过程共同作用的结果，不再是浮游植物占主导的系统了。
上下游对比来看，上游是以浮游植物光合作用产生氧气为主的自养型系统，下游已经成为耗氧过程占主导的系统。这个结果与黄静水等人直接通过水质模型计算DO平衡的结果相似。这种情况下，由于污水处理厂尾水的大量排入以及沿岸污染物的进入，会使下游水体中的DO持续降低，甚至出现黑臭现象。
3.6 总体分析
从总体纳什系数来看，上游最敏感的参数依次是K’e、DP、f，而且K’e的敏感性和交互作用远大于其他的参数，这说明上游河段的模拟结果对光最敏感。排名前十的参数中，除了SOD，其余9个参数都直接与浮游植物相关，充分反映了污水处理厂上游系统是一个以浮游植物生长死亡为主导的偏自然的河流系统。
下游最敏感的参数是K2D和KNO3，而且它们的σ值也排在前两位，说明反硝化作用对下游的模拟结果最敏感。除此之外，下段还主要受SOD和浮游植物生长死亡的影响，各个过程的影响在下游均有不同程度的体现。
以总体模拟的纳什系数为目标函数进行参数敏感性分析，发现污水处理厂上下游的参数敏感性发生了显著变化。上游是以浮游植物生长死亡为主导的自养型系统。下游由于河岸带变化、污水处理厂稀释等原因，浮游植物的浓度急剧减少，不再以浮游植物为主导。由于污水处理厂排放大量的硝酸盐进入河道，以及支流或者排口排入的大量有机物，导致下游河道溶解氧显著降低，反硝化作用强烈。这些变化反映了由于污水处理厂尾水的大量排入，打破了河流水质过程的连续性，在尾水排入的地方水质发生突变，河流上下游的主导过程发生了显著改变，一定程度可以揭示污水处理厂对河道系统带来的巨大影响。
4 结论
本研究选择在MATLAB/Simulink环境下重建南淝河水质模型，并结合SAFE工具包，实现对41个参数的空间敏感性分析发现参数对不同水质子系统敏感性存在差异，这些差异反映了水质系统的复杂性。同时，参数敏感性的空间差异，反映了水质过程的空间差异和研究对象的系统特征，揭示了污水处理厂尾水作为主要补给源对城市河道的巨大影响。污水处理厂上游是以浮游植物生长死亡为主的偏自然河道，河流保持在比较健康的状态;下游转变成以反硝化作用、硝化作用、SOD耗氧过程为主的受人为活动影响严重的河道，水质污染逐渐加重，甚至出现黑臭的情况。该方法很好地解决了WASP软件无法进行敏感性分析，以及无法设置分布式参数的弊端，是研究河道水质过程和系统特征的又一有效手段。
5 建议
根据分析结果，对南淝河治理提出整体性建议，针对上游河道，治理主要以截污为主，杜绝未处理生活污水直排的现象，增加河道上游清洁水来源，即可保持上游较好的水质状态。针对下游河道的治理，有以下几点建议：(1)进一步改进污水处理厂处理工艺，降低排放尾水中的硝氮和有机物负荷;(2)利用排放河道以及周边可以利用的条件尽可能减少尾水对河道水质的影响，比如尾水河道旁侧生态处理、生态生物处理模式等;(3)提升支流来水水质;(4)恢复下游河道的生态护岸，在护岸上种植水生植物和其他植物，增加下游河道的生物净化作用，同时营造柔美生态岸线，集防洪、生态、景观和自净等功能于一体。
原标题：净水技术|南淝河城市河道水质模型的空间敏感性分析
