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

謝如意，黃靜水，周　琪

(同濟大學環境科學與工程學院，上海　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.