謝如意，黃靜水，周　琪

(同濟(jì)大學(xué)環(huán)境科學(xué)與工程學(xué)院，上海　200092)

南淝河位于安徽省合肥市，是巢湖污染最嚴(yán)重的支流之一。合肥是安徽省的政治、經(jīng)濟(jì)和文化中心。2005年～2015年，合肥市人口增長(zhǎng)了70%，達(dá)到了780萬，國(guó)內(nèi)生產(chǎn)總值增長(zhǎng)了560%,達(dá)到了5 660億元。作為城市主要的供水來源和收納水體，城市的快速發(fā)展給南淝河帶來了巨大的壓力，水質(zhì)修復(fù)成為需要解決的重要問題。

水質(zhì)模型是進(jìn)行河道水質(zhì)研究的重要工具之一，廣泛應(yīng)用于情境管理和水質(zhì)響應(yīng)研究等方面，水質(zhì)模型可以為水質(zhì)管理提供決策支持。WASP、QUAL2K、EFDC、MIKE等都是現(xiàn)在常用的水質(zhì)模型。但隨著水質(zhì)模型日益復(fù)雜，模型采用的參數(shù)也越來越多，而且大部分參數(shù)是無法精確測(cè)定的，為減少模型率定的困難，可以著重研究一部分對(duì)模型輸出變化敏感的參數(shù)。

敏感性分析研究的是模型輸入因子的變化如何影響模型輸出的變化[1]。敏感性分析又分為局部敏感性和全局敏感性。全局敏感性對(duì)參數(shù)在整個(gè)取值空間的敏感性都進(jìn)行分析，因此可以對(duì)參數(shù)的相互作用以及輸入因子和模型輸出的非線性關(guān)系進(jìn)行研究。敏感性分析不僅為模型率定提供依據(jù)，也是對(duì)研究對(duì)象進(jìn)行系統(tǒng)認(rèn)識(shí)的有效工具之一。

在國(guó)內(nèi)，敏感性分析方法在水質(zhì)模型上已經(jīng)有了一定的應(yīng)用。張永祥等[2]將WASP模型應(yīng)用于長(zhǎng)河水體富營(yíng)養(yǎng)化分析，采用的是單點(diǎn)擾動(dòng)的局部敏感性分析方法，擾動(dòng)幅度高達(dá)50%，未能很好地反映參數(shù)的真實(shí)敏感性。張質(zhì)明等[3]以WASP水質(zhì)模型的應(yīng)用為例，通過Sobol方法確定模型的敏感參數(shù)，提出了一套基于GLUE法的多目標(biāo)模型參數(shù)率定方法。Yi等[4]對(duì)建立的滇池水質(zhì)模型的47個(gè)參數(shù)和7個(gè)外部擾動(dòng)因子進(jìn)行時(shí)間和空間的敏感性分析，闡明了敏感性分析對(duì)模型參數(shù)識(shí)別和不確定分析的重要性。本文的研究目的，是通過MATLAB的Simulink模擬仿真平臺(tái)重建WASP水質(zhì)模型的富營(yíng)養(yǎng)化模塊，結(jié)合SAFE工具包對(duì)多參數(shù)過程的復(fù)雜水質(zhì)模型進(jìn)行空間全局敏感性分析，探究水質(zhì)模型參數(shù)的空間差異，進(jìn)而揭示南淝河河道系統(tǒng)水質(zhì)變化特征。

1　研究方法

1.1　WASP水質(zhì)模型的選取

WASP是美國(guó)環(huán)保署推薦使用的水質(zhì)模擬軟件，是為分析池塘、湖泊、水庫、河口和沿海水域的一系列水質(zhì)問題而設(shè)計(jì)的多箱模型，可模擬常規(guī)污染物和有毒污染物在水中的遷移轉(zhuǎn)化規(guī)律，被稱為“萬能水質(zhì)模型”。

WASP的水質(zhì)模擬分為兩個(gè)模塊，有毒物質(zhì)模型TOXI和富營(yíng)養(yǎng)化模型EUTRO。EUTRO模塊用以模擬傳統(tǒng)污染物的遷移轉(zhuǎn)化規(guī)律，分為溶解氧、氮、磷和浮游植物四個(gè)子系統(tǒng)(圖1)，包括8個(gè)狀態(tài)變量和若干反應(yīng)過程(表1)，具體的反應(yīng)方程見WASP用戶指導(dǎo)手冊(cè)[5]。

圖1　EUTRO模塊各子系統(tǒng)及狀態(tài)變量[5]

過程N(yùn)H+4NO-3DIPPHYTCBODDODONDOP硝化作用√√√反硝化作用√√有機(jī)氮礦化√√浮游植物死亡√√√√√√浮游植物生長(zhǎng)吸收√√√√√有機(jī)磷礦化√√氧化作用√√沉降√√√√大氣復(fù)氧√底泥耗氧√內(nèi)源呼吸√

1.2　Simulink環(huán)境下水質(zhì)模型的構(gòu)建

WASP自身并沒有敏感性分析功能，只能借助其他的工具進(jìn)行敏感性分析。同時(shí)，自WASP5后的版本不再提供源碼，而且WASP本身是有界面的封裝好的水質(zhì)模擬軟件，當(dāng)需要多次批量修改參數(shù)運(yùn)行調(diào)用時(shí)，操作十分不方便。最后，并沒有專門針對(duì)WASP模型開發(fā)的敏感性分析工具，直接對(duì)WASP水質(zhì)模型進(jìn)行敏感性分析存在比較多的技術(shù)困難和操作不便性。

Simulink是MATLAB內(nèi)部的可視化仿真工具，可用于實(shí)現(xiàn)動(dòng)態(tài)系統(tǒng)的建模、仿真與分析。MATLAB/Simulink的主要優(yōu)點(diǎn)有：(1)可視化的直觀建模方式，簡(jiǎn)單明了，且各模塊可以進(jìn)行封裝，大大增加了模型的可讀性；(2)模型設(shè)計(jì)和模擬快速準(zhǔn)確，還為用戶提供了一個(gè)圖形化的調(diào)試工具以輔助用戶進(jìn)行系統(tǒng)開發(fā)；(3)MATLAB/Simulink的仿真程序可以與MATLAB工具包實(shí)現(xiàn)良好的兼容，適于作為進(jìn)一步研究和開發(fā)工藝控制程序的最佳平臺(tái)。

為便于對(duì)模型進(jìn)行敏感性分析，將WASP6.0模型按照模型手冊(cè)中的內(nèi)容在MATLAB/Simulink環(huán)境下重新建模。第一步，將EUTRO模塊中4個(gè)子系統(tǒng)8個(gè)狀態(tài)變量的所有子過程寫入Simulink中，如圖2a、圖2b所示。第二步，建立對(duì)流項(xiàng)、彌散項(xiàng)和源匯項(xiàng)的方程(圖2c)。因?yàn)楸狙芯筷P(guān)注水質(zhì)模型，所以水動(dòng)力部分在EPDRiv1模型中算出后，直接在本模型中使用，具體的結(jié)果見Huang等的文章[6]。第三步，將所有方程封裝為一個(gè)“塊”，使用者只能對(duì)每個(gè)“塊”的參數(shù)進(jìn)行修改，而無法直接修改內(nèi)部方程，這樣可以提供更為簡(jiǎn)潔的使用界面。最后，根據(jù)南淝河的實(shí)際情況，將研究河段劃分為45個(gè)小段，每段長(zhǎng)約200～500 m，采用一個(gè)“塊”來模擬，然后將45個(gè)“塊”連接起來，建立南淝河水質(zhì)模型(圖2e)。水質(zhì)模型構(gòu)建完成后，為檢查新模型的正確性，直接采用在WASP模型中手動(dòng)率定的參數(shù)值，檢驗(yàn)Simulink模型的正確性，部分參數(shù)值與WASP參數(shù)有區(qū)別。重建模型的模擬結(jié)果與WASP手動(dòng)率定結(jié)果相關(guān)性如表2所示，各個(gè)指標(biāo)的模擬結(jié)果相關(guān)性比較高，說明重構(gòu)的南淝河水質(zhì)模型結(jié)構(gòu)上是準(zhǔn)確的，可以在此基礎(chǔ)上對(duì)模型敏感性進(jìn)行分析。

圖2　Simulink環(huán)境下建立南淝河水質(zhì)模型

指標(biāo)NH+4NO-3DIPChl?aDO模擬結(jié)果相關(guān)性098099074077099

選用Simulink工具對(duì)WASP水質(zhì)模型進(jìn)行重建，有以下有幾個(gè)優(yōu)點(diǎn)：第一，模型運(yùn)行和敏感性分析工具包集合在MATLAB統(tǒng)一平臺(tái)上，可以直接實(shí)現(xiàn)敏感性分析；第二，可以有針對(duì)性地對(duì)水質(zhì)模型結(jié)構(gòu)、參數(shù)設(shè)置等進(jìn)行調(diào)整，比如WASP中大多數(shù)參數(shù)是集總式的，無法對(duì)同一參數(shù)在不同位置設(shè)置不同的值。本研究為對(duì)模型參數(shù)進(jìn)行空間敏感性分析，在重建模型時(shí)，將模型參數(shù)改進(jìn)為分布式，即不同河段的參數(shù)值可以不同。

表3　敏感性分析的參數(shù)及其取值范圍

注：(a)表示參數(shù)范圍來自Rates,Constants,and Kinetics Formulations in Surface Water Quality Modeling[7](第二版)一書；(b)表示參數(shù)范圍來自《WASP6用戶手冊(cè)》[5]推薦值

SAFE(sensitivity analysis for everybody)是一款包含眾多敏感性分析方法及可視化函數(shù)的MATLAB工具包[8]。其提供的敏感性分析方法包括：Morris法、RSA法、Sobol法、FAST法、動(dòng)態(tài)識(shí)別分析法和基于密度的敏感性分析方法等。SAFE不僅提供了多種全局敏感性分析方法，還能夠輕松地加入其他的敏感性分析方法。其提供的敏感性分析方法，都支持對(duì)敏感性指標(biāo)進(jìn)行魯棒性和收斂性的評(píng)價(jià)。SAFE工具包中還提供了大量的可視化工具，可以對(duì)敏感性分析結(jié)果進(jìn)行更好的研究和相互對(duì)比。這款工具包既能讓非專業(yè)人員方便地操作，也允許有經(jīng)驗(yàn)的使用者在此基礎(chǔ)上開發(fā)新的功能。

本研究結(jié)合SAFE工具包中提供的LH-OAT采樣方法及Morris敏感性分析方法，以望塘污水處理廠為分界點(diǎn)，將研究河段分為上下游，分別進(jìn)行參數(shù)敏感性分析。整個(gè)模型參與敏感性分析的參數(shù)如表3所示，其參數(shù)取值范圍來自于文獻(xiàn)報(bào)道中的值。

1.3　Morris分析方法

Morris法[9]是一種定性的全局敏感性分析方法，主要用于參數(shù)的篩選和排序。Morris法屬于一次只改變一個(gè)變量的方法，即選取模型某一變量xi，其余參數(shù)值固定不變，在變量取值變化范圍內(nèi)隨機(jī)改變xi的值，運(yùn)行模型得到目標(biāo)函數(shù)值，用模型輸出對(duì)模型輸入的變化率EEs來表示參數(shù)變化對(duì)輸出值的影響程度，多次變化后的平均變化率作為全局敏感性的指標(biāo)如式(1)。

(1)

其中：Ci-比例因子;

g-模型輸出;

M-參與敏感性分析的參數(shù)個(gè)數(shù);

i-第i個(gè)輸入?yún)?shù);

j-第j次模型運(yùn)算;

r-模型運(yùn)行總次數(shù);

Si-全局敏感性參數(shù)。

EEs的平均值(μ)代表參數(shù)對(duì)輸出結(jié)果總體敏感性大小，EEs的標(biāo)準(zhǔn)差(σ)代表某個(gè)輸入?yún)?shù)與其他參數(shù)相互作用的程度。高標(biāo)準(zhǔn)差代表某個(gè)參數(shù)的敏感性在其整個(gè)取值空間內(nèi)都發(fā)生變化。Morris法的敏感性指數(shù)u和σ都只反映參數(shù)敏感性的相對(duì)大小，只能用于判斷敏感性的相對(duì)大小，其值并不反映參數(shù)的真實(shí)敏感性[1]。

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

2.1　研究河段的基本概況

南淝河長(zhǎng)約70 km，流域集水面積約1 446 km2，最終匯入巢湖。董鋪水庫至橡膠壩一段流經(jīng)合肥城區(qū)，長(zhǎng)約17 km，如圖3所示，本研究模擬河段為圖中紅色河段，長(zhǎng)約12 km。

圖3　研究區(qū)域圖(黑色河段為本研究模擬河段)

董鋪水庫下游約7 km處，河流兩岸完全渠道化，為硬質(zhì)河岸。城市河段內(nèi)有兩條支流：一條支流為四里河，其上游為大房郢水庫；另一條支流為板橋河，沿岸主要為農(nóng)業(yè)用地。董鋪水庫是合肥市主要的飲用水源地，攔截了南淝河上游的大部分清潔來水，導(dǎo)致南淝河城市河段基本沒有清潔水來源。研究河段內(nèi)有兩個(gè)污水處理廠，望塘污水處理廠位于城市河段的上游，年均排放量約為2.21 m3/s；王小郢污水處理廠位于橡皮壩上游約1 km處，年均排放量約3.66 m3/s。污水處理廠尾水是南淝河城市河段水量的主要來源，旱季和雨季分別占河道流量的75%和53%。董鋪水庫下游3 km處，有一城中村，其未處理的生活污水直排進(jìn)入南淝河。因此南淝河城區(qū)段是一段典型的以污水處理廠尾水為主要補(bǔ)給的城市化河段。

2.2　模型建立

3　空間敏感性分析結(jié)果與討論

以望塘污水處理廠為分界點(diǎn)，分別對(duì)上下游河段進(jìn)行敏感性分析，Morris分析結(jié)果如圖4所示，橫軸u*代表參數(shù)的總體敏感性，縱軸σ代表參數(shù)間相互作用大小。圖4最后一組NSE_T的目標(biāo)函數(shù)為五個(gè)指標(biāo)的總體納什系數(shù)，代表上下游模型整體輸出對(duì)參數(shù)的敏感性情況。表4表示的是對(duì)各個(gè)指標(biāo)及總體模擬結(jié)果中敏感性排名前10的參數(shù)排序。

圖4　單指標(biāo)及總體模擬效果納什系數(shù)Morris敏感性分析結(jié)果圖

排序NH+4NO-3DIPChl?aDO總體USDSUSDSUSDSUSDSUSDSUSDS1SODKmNcSODK2DDpK＇eDpK＇eK＇eSODK＇eK2D2DpK71K2DKNO3fDpfDpSODK＇eDpKNO33fKNITKNITSODK＇efK＇efDpK2DfDp4K＇eK12DpE2DK1caPCK1cK1caPCKNO3SODK＇e5KNITDpfK＇eIaKmPcIaIafEDK1cSOD6K1cK＇eK＇eESaPCk83vs4IsPNH3KDPNH3f7IaE12fD5EDIsfopaPCE1CIafIaK1c8aPCfONaNCKDvs4IaIsvs4K12ESaPCKmNc9PNH3PNH3KNO3fD5E1CIsE1CKmPK1cfD5vs4Ia10fD5aNCK12vs3fopK1cKmPEcIsK1cIsIs

注：US代表污水處理廠上游，DS代表污水處理廠下游

3.1　氨氮指標(biāo)

3.2　硝酸鹽指標(biāo)

3.3　無機(jī)磷指標(biāo)

污水處理廠上游，DIP濃度對(duì)浮游植物死亡速率最敏感(DP)，其次是與浮游植物生長(zhǎng)相關(guān)的一系列參數(shù)。同時(shí)DIP在上游的參數(shù)敏感性排序與Chl-a的參數(shù)敏感性排序有很高的一致性，這主要是因?yàn)镈IP與其他子系統(tǒng)的聯(lián)系很少，只與浮游植物子系統(tǒng)發(fā)生密切聯(lián)系，所以參數(shù)排序出現(xiàn)比較高的一致性，同樣的結(jié)果在Yi等滇池水質(zhì)模型的敏感性分析中也提到過[4]。

污水處理廠下游，DIP除了對(duì)與浮游植物相關(guān)的一系列參數(shù)敏感，還對(duì)與有機(jī)磷礦化的部分參數(shù)敏感，如K83、E83和KmPc等。由圖4可知，有機(jī)磷礦化相關(guān)參數(shù)的σ值都比浮游植物相關(guān)參數(shù)的低，這也說明了浮游植物系統(tǒng)參數(shù)的相互作用非常大。

上下兩段對(duì)比來看，磷系統(tǒng)自身的轉(zhuǎn)化過程相對(duì)穩(wěn)定，參數(shù)的相互作用程度低，主要是受浮游植物生長(zhǎng)與死亡的影響，上下游均有這個(gè)特點(diǎn)。污水處理廠對(duì)磷系統(tǒng)的影響不大。

3.4　浮游植物指標(biāo)

污水處理廠上游，Chl-a濃度對(duì)Dp最敏感，而且其參數(shù)間相互作用也最強(qiáng)。上游的敏感參數(shù)的u*和σ有很好的線性關(guān)系，與DIP的結(jié)果有一致性，這也說明了浮游植物相關(guān)參數(shù)的敏感性和相互作用可能存在內(nèi)在的線性關(guān)系。

對(duì)比來看，上下游Chl-a濃度主要受浮游植物自身相關(guān)的一些參數(shù)的影響，說明在富營(yíng)養(yǎng)化模塊中，浮游植物系統(tǒng)是最核心最重要的子系統(tǒng)。其次，上下游敏感性排名前10的參數(shù)中都出現(xiàn)了KmP，有可能上下游浮游植物的營(yíng)養(yǎng)限制都是磷限制。與浮游植物相關(guān)的各個(gè)過程：生長(zhǎng)包括的三個(gè)限制項(xiàng)、死亡和沉降過程都在敏感性參數(shù)中有體現(xiàn)。這說明，浮游植物濃度是一個(gè)極其敏感的指標(biāo)，與之相關(guān)的各個(gè)過程都直接地影響其濃度，更容易出現(xiàn)異參同效的情況，所以在參數(shù)識(shí)別和率定的時(shí)候要格外注意，模擬準(zhǔn)確的難度也更大。

3.5　溶解氧指標(biāo)

上下游對(duì)比來看，上游是以浮游植物光合作用產(chǎn)生氧氣為主的自養(yǎng)型系統(tǒng)，下游已經(jīng)成為耗氧過程占主導(dǎo)的系統(tǒng)。這個(gè)結(jié)果與Huang等直接通過水質(zhì)模型計(jì)算DO平衡的結(jié)果相似[6]。這種情況下，由于污水處理廠尾水的大量排入以及沿岸污染物的進(jìn)入，會(huì)使下游水體中的DO持續(xù)降低，甚至出現(xiàn)黑臭現(xiàn)象。

3.6　小結(jié)

下游最敏感的參數(shù)是K2D和KNO3，而且它們的σ值也排在前兩位，說明反硝化作用對(duì)下游的模擬結(jié)果最敏感。除此之外，下段還主要受SOD和浮游植物生長(zhǎng)死亡的影響，各個(gè)過程的影響在下游均有不同程度的體現(xiàn)。

以總體模擬的納什系數(shù)為目標(biāo)函數(shù)進(jìn)行參數(shù)敏感性分析，發(fā)現(xiàn)污水處理廠上下游的參數(shù)敏感性發(fā)生了顯著變化。上游是以浮游植物生長(zhǎng)死亡為主導(dǎo)的自養(yǎng)型系統(tǒng)。下游由于河岸帶變化、污水處理廠稀釋等原因，浮游植物的濃度急劇減少，不再以浮游植物為主導(dǎo)。污水處理廠排放大量的硝酸鹽進(jìn)入河道，以及支流或者排口排入的大量有機(jī)物，導(dǎo)致下游河道溶解氧顯著降低，反硝化作用強(qiáng)烈。這些變化反映了污水處理廠尾水的大量排入，打破了河流水質(zhì)過程的連續(xù)性，在尾水排入的地方水質(zhì)發(fā)生突變，河流上下游的主導(dǎo)過程發(fā)生了顯著改變，一定程度可以揭示污水處理廠對(duì)河道系統(tǒng)帶來的巨大影響。

4　結(jié)論

本研究選擇在MATLAB/Simulink環(huán)境下重建南淝河水質(zhì)模型，并結(jié)合SAFE工具包，實(shí)現(xiàn)對(duì)41個(gè)參數(shù)的空間敏感性分析發(fā)現(xiàn)參數(shù)對(duì)不同水質(zhì)子系統(tǒng)敏感性存在差異，這些差異反映了水質(zhì)系統(tǒng)的復(fù)雜性。同時(shí)，參數(shù)敏感性的空間差異，反映了水質(zhì)過程的空間差異和研究對(duì)象的系統(tǒng)特征，揭示了污水處理廠尾水作為主要補(bǔ)給源對(duì)城市河道的巨大影響。污水處理廠上游是以浮游植物生長(zhǎng)死亡為主的偏自然河道，河流保持在比較健康的狀態(tài)；下游轉(zhuǎn)變成以反硝化作用、硝化作用、SOD耗氧過程為主的受人為活動(dòng)影響嚴(yán)重的河道，水質(zhì)污染逐漸加重，甚至出現(xiàn)黑臭的情況。該方法很好地解決了WASP軟件無法進(jìn)行敏感性分析，以及無法設(shè)置分布式參數(shù)的弊端，是研究河道水質(zhì)過程和系統(tǒng)特征的又一有效手段。

5　建議

根據(jù)分析結(jié)果，對(duì)南淝河治理提出整體性建議，針對(duì)上游河道，治理主要以截污為主，杜絕未處理生活污水直排的現(xiàn)象，增加河道上游清潔水來源，即可保持上游較好的水質(zhì)狀態(tài)。針對(duì)下游河道的治理，有以下幾點(diǎn)建議：(1)進(jìn)一步改進(jìn)污水處理廠處理工藝，降低排放尾水中的硝氮和有機(jī)物負(fù)荷；(2)利用排放河道以及周邊可以利用的條件盡可能減少尾水對(duì)河道水質(zhì)的影響，比如尾水河道旁側(cè)生態(tài)處理、生態(tài)生物處理模式等；(3)提升支流來水水質(zhì)；(4)恢復(fù)下游河道的生態(tài)護(hù)岸，在護(hù)岸上種植水生植物和其他植物，增加下游河道的生物凈化作用，同時(shí)營(yíng)造柔美生態(tài)岸線，集防洪、生態(tài)、景觀和自凈等功能于一體。

[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模型參數(shù)率定與敏感性分析[J].水資源與水工程學(xué)報(bào),2009,20(5):28-30.

[3]張質(zhì)明,王曉燕,于洋,等.基于GLUE法的多指標(biāo)水質(zhì)模型參數(shù)率定方法[J].環(huán)境科學(xué)學(xué)報(bào),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.