張質(zhì)明,王曉燕,李明濤 (1.北京建筑大學(xué)北京應(yīng)對(duì)氣候變化和人才培養(yǎng)基地,北京 100044；.首都師范大學(xué)資源環(huán)境與旅游學(xué)院,北京 100048)

基于全局敏感性分析方法的WASP模型不確定性分析

張質(zhì)明1,2,王曉燕2*,李明濤2(1.北京建筑大學(xué)北京應(yīng)對(duì)氣候變化和人才培養(yǎng)基地,北京 100044；2.首都師范大學(xué)資源環(huán)境與旅游學(xué)院,北京 100048)

為了定量討論WASP(The water quality analysis simulation program,水質(zhì)分析模擬程序)模型中各參數(shù)對(duì)模擬結(jié)果的影響及其不確定性,本文在Simulink環(huán)境下對(duì)模型進(jìn)行重組,利用Sobol方法,對(duì)影響DO、CBOD、NH3-N、NO3--N等4項(xiàng)輸出的各參數(shù)敏感性進(jìn)行了研究,重點(diǎn)討論了輸入條件如時(shí)空變化對(duì)參數(shù)敏感性的影響.結(jié)果表明:Sobol全局敏感性分析方法能夠有效地篩選出 WASP模型中對(duì) DO、CBOD、NH3-N、NO變化模擬較為敏感的參數(shù),在實(shí)現(xiàn)模型參數(shù)“本土化”中具有較大的應(yīng)用潛力;對(duì)比總敏感性與一階敏感性指數(shù)之間的差異, WASP模型當(dāng)中的若干參數(shù)(如GP1、PNH3、anc等)耦合性較強(qiáng),不宜使用一次一個(gè)變量法對(duì)這些變量進(jìn)行敏感性分析;在參數(shù)敏感性變化方面,即使對(duì)于同一條河流,WASP的參數(shù)敏感性指數(shù)會(huì)因邊界條件的時(shí)空變化而發(fā)生改變.敏感性受空間變化的影響不大,但受時(shí)間(季節(jié))變化的影響顯著;輸入變量隨時(shí)間的變化,會(huì)引起模型不確定性在各季節(jié)上的較大的差異;WASP模型當(dāng)中存在明顯的“異參同效”作用,僅靠傳統(tǒng)方法進(jìn)行全局尋優(yōu)率定,并不能體現(xiàn)水體各組分之間相互轉(zhuǎn)化的過(guò)程,還需要結(jié)合實(shí)驗(yàn)對(duì)關(guān)鍵參數(shù)進(jìn)行率定.

不確定性分析；敏感性；Sobol方法；WASP模型；邊界條件

水質(zhì)模型中的不確定性研究已經(jīng)成為當(dāng)前水質(zhì)模擬領(lǐng)域中的重要問(wèn)題[1-7].由于在水質(zhì)模型的應(yīng)用中,眾多的可調(diào)參數(shù)會(huì)造成一定的不確定性,特別是當(dāng)“異參同效”[8-10]現(xiàn)象出現(xiàn)時(shí),多種參數(shù)組合令模型結(jié)果均可滿足模型的最優(yōu)條件,無(wú)法通過(guò)水質(zhì)模型準(zhǔn)確解析實(shí)際中的污染物遷移轉(zhuǎn)化規(guī)律.

為識(shí)別造成不確定性的主要參數(shù)[11-12],需要進(jìn)行敏感性分析[13-14],方法主要分為局部敏感性分析與全局敏感性分析兩種[15-16].局部敏感性分析的優(yōu)勢(shì)在于計(jì)算量小[17],但由于該方法建立在模型輸出與參數(shù)的一階偏導(dǎo)數(shù)關(guān)系上,對(duì)參數(shù)取值鄰域的可導(dǎo)性提出了要求.為了能夠更為準(zhǔn)確的評(píng)估參數(shù)敏感性,現(xiàn)在更多的研究?jī)A向于使用全局敏感性分析方法.目前常見(jiàn)的全局敏感性分析方法包括定性的Morris法、FAST法、GLUE法以及定量的Extend FAST法、Sobol法、基于ANN的權(quán)值分析法等[18].其中Sobol方法基于方差分解的原理,可用于非線性、非單調(diào)的數(shù)學(xué)模型、結(jié)果穩(wěn)健可靠、并且能夠給出對(duì)參數(shù)敏感性的定量評(píng)價(jià),已經(jīng)廣泛應(yīng)用于環(huán)境及其他領(lǐng)域大型非線性模型中[19-23].

目前針對(duì)參數(shù)對(duì)最終結(jié)果的所造成的影響方面的研究較多,而對(duì)模型內(nèi)部各個(gè)子模塊所描述的污染物轉(zhuǎn)化過(guò)程所具有的不確定性方面的研究并不多見(jiàn),且關(guān)于參數(shù)敏感性隨時(shí)空條件變化會(huì)發(fā)生什么樣的改變的研究并不多.為討論WASP水質(zhì)模型在不同時(shí)空條件下的重點(diǎn)參數(shù)及其所引起的不確定性,本文以北運(yùn)河通州段為例,利用聚類分析方法確定出河道污染特征,作為空間上的劃分依據(jù);按照春、夏、秋、冬四季作為時(shí)間上的劃分依據(jù),分別對(duì) DO、CBOD、氨氮、硝態(tài)氮等 4個(gè)方面的模擬進(jìn)行敏感性分析,基于Sobol法識(shí)別出各條件下的重要敏感參數(shù),用來(lái)研究不同輸入條件下模型參數(shù)的敏感性變化.之后對(duì)WASP的分解模型使用蒙特卡羅模擬,確認(rèn)參數(shù)不確定性對(duì)模型中各個(gè)模塊所造成影響的大小,獲取模型中的敏感模塊.旨在為WASP模型的進(jìn)一步優(yōu)化、減少不確定性方面提供理論依據(jù).

1 數(shù)據(jù)與方法

1.1 分析方法

1.1.1 Sobol方法 Sobol方法是基于方差分解的思想對(duì)敏感性進(jìn)行評(píng)估的[24]:

Sobol證明過(guò)分解式的唯一性并且所有分解項(xiàng)均可以通過(guò)多重積分求得[25]總方差為

偏方差可以通過(guò)分解式計(jì)算得到

敏感度系數(shù)1,2,...,sii iS 為

Si為 xi的一階敏感性指數(shù),用于定量描述 xi在函數(shù)f(x)中所造成的影響;叫做因素的s階敏感性指數(shù),用于定量描述這s個(gè)參數(shù)共同作用對(duì)函數(shù) f(x)的影響.因此對(duì)于一個(gè) s個(gè)參數(shù)的模型來(lái)說(shuō),變量總敏感性指數(shù)可以表示為:

1.1.2 敏感性度量目標(biāo)函數(shù) 目前關(guān)于模型參數(shù)敏感性度量主要從兩個(gè)方面入手:一種是觀察參數(shù)變化對(duì)模型輸出值大小的直接影響,另外一種是通過(guò)參數(shù)變化對(duì)似然函數(shù)值造成的影響[26].本研究分別采用這 2種方法進(jìn)行計(jì)算,似然函數(shù)使用平均相對(duì)誤差以及Nash-Suttcliffe系數(shù)來(lái)衡量模擬值與觀測(cè)值之間的擬合度,納什系數(shù)表達(dá)式為:

式中:Qob為觀測(cè)值,Qsim為模擬值,Qob_average為觀測(cè)平均值,n為觀測(cè)的次數(shù).當(dāng) Qob=Qsim時(shí),Ens=1;若Ens＜0,則說(shuō)明模型模擬無(wú)效.

1.1.3 模型的結(jié)構(gòu)分解及參數(shù) Simulink是MATLAB的可視化仿真工具,它基于 MATLAB的框圖設(shè)計(jì)環(huán)境,可用于實(shí)現(xiàn)動(dòng)態(tài)系統(tǒng)建模、仿真與分析.由于WASP在5.0之后的版本不再提供源碼,為便于對(duì)模型的敏感性分析與參數(shù)率定,將 WASP6.0模型按照模型手冊(cè)中的內(nèi)容基于MATLAB/Simulink環(huán)境進(jìn)行了重寫(xiě).將 DO、CBOD、氨氮、硝氮的各個(gè)子過(guò)程(如:與 CBOD有關(guān)的氧化、死亡、反硝化、沉淀等作用機(jī)理)分別寫(xiě)入各自的m-function并定義其輸出變量,最終集成在.mdl文件內(nèi).在本研究中,通過(guò)蒙特卡羅模擬方法,循環(huán)利用MATLAB中的sim命令調(diào)用該模型,獲取自定義的輸出變量獲取各子過(guò)程的模擬結(jié)果,確定模型及其各子模塊的輸出范圍,用于對(duì)模型的不確定性來(lái)源進(jìn)行分析.

建模所涉及的機(jī)理過(guò)程、模型參數(shù)如表 1所示:

表1 模型所涉及的機(jī)理過(guò)程與參數(shù)Table 1 Mechanism processes and parameters referred to in the simulation of the four indicators

WASP應(yīng)用的范圍很廣,很多河流都有WASP研究的案例,參數(shù)的選取范圍也各不相同.特別是對(duì)于水動(dòng)力條件受閘壩人為干擾明顯的河流,參數(shù)的取值范圍可能會(huì)與一般河流不同:例如水閘截流時(shí)的情形,參數(shù)范圍就應(yīng)當(dāng)變化,甚至選取一些適用于湖泊水體的參數(shù)作為參考.因此本文模型參數(shù)的范圍參照各類水體中的研究案例[27-31].

1.2 研究區(qū)域與數(shù)據(jù)來(lái)源

研究區(qū)北運(yùn)河為唯一發(fā)源于北京境內(nèi)的水系,閘壩多,流速緩慢,天然水量較少,多來(lái)自污水處理廠處理排放的再生水.課題組在 2009年4～10月以及2010年的4月、7月10月、12月對(duì)北運(yùn)河通州段干流附近的共計(jì)14個(gè)樣點(diǎn)的流量、流速以及水質(zhì)(包括COD、CBOD、氨氮、硝態(tài)氮、DO、水溫等指標(biāo))進(jìn)行了同步采樣監(jiān)測(cè).其中采樣點(diǎn)布設(shè)及河段概化示意如圖1所示.

圖1 采樣點(diǎn)空間分布與河道概化Fig.1 Location of monitoring sites of North Canal in Tongzhou

如圖 1所示,監(jiān)測(cè)點(diǎn)的重要排污口以及匯入的支流均加設(shè)置了監(jiān)測(cè)點(diǎn).根據(jù)監(jiān)測(cè)數(shù)據(jù),北運(yùn)河通州段總磷濃度平均為0.70mg/L,超過(guò)地表水V類水標(biāo)準(zhǔn)[32].COD濃度為16.14～79.00mg/L,平均濃度為45.38mg/L,超過(guò)了地表水V類水標(biāo)準(zhǔn).氨氮濃度為 1.88～25.80mg/L,占總氮的 78.0%,是北運(yùn)河氮污染的主要指標(biāo).

2 結(jié)果與討論

2.1 水質(zhì)時(shí)空分布特征分析

普遍認(rèn)為,水質(zhì)變化的速率與水溫、水力條件、水質(zhì)狀況等因素有很大的關(guān)系.在 WASP模型中,水質(zhì)的變化速率的計(jì)算完全取決于這些狀態(tài)變量以及相關(guān)常值參數(shù).在模型計(jì)算過(guò)程中,狀態(tài)變量取值差異,將會(huì)導(dǎo)致相關(guān)參數(shù)敏感性的差異.

在主要的狀態(tài)變量中,水溫主要受到季節(jié)因素的影響;水力條件主要受到水系匯流狀況、季節(jié)因素以及閘壩控制的影響;河道水質(zhì)狀況(包括 DO、COD、氨氮、BOD、總磷等指標(biāo))除了受到季節(jié)因素的影響[33]以外,還與河段流經(jīng)區(qū)域的土地利用有關(guān)[34-36],因此狀態(tài)變量的變異可以分解為時(shí)間、空間兩方面分別進(jìn)行討論.

為了研究水質(zhì)污染狀況隨空間變化是否引起模型參數(shù)敏感性的變化,本研究采用聚類分析法來(lái)識(shí)別北運(yùn)河不同河段的污染物組成特征.根據(jù)對(duì)監(jiān)測(cè)水質(zhì)、流量數(shù)據(jù)的相關(guān)性分析,發(fā)現(xiàn)在諸多的狀態(tài)變量中存在一定的相關(guān)性.為了避免重復(fù)考慮這些因素,篩選出相對(duì)較為獨(dú)立的DO、氨氮、COD等3個(gè)指標(biāo)歸一化后進(jìn)行聚類分析,結(jié)果如圖2所示.

圖2 采樣點(diǎn)聚類分析結(jié)果Fig.2 Dendrogram showing clustering of sampling sites of the North Canal

如圖2所示,除8號(hào)采樣點(diǎn)以外,其余采樣點(diǎn)監(jiān)測(cè)到的水質(zhì)特征隨與空間的變化而變.采樣點(diǎn)1與采樣點(diǎn)2為代表的人口密集的居住區(qū)納污河段(類型 I),工業(yè)污染源較少,主要污染物來(lái)源為集中處理設(shè)施;3、4、5、6、7代表城鄉(xiāng)結(jié)合部的水質(zhì)特征(類型II),9、10、11、12、13、14代表流經(jīng)農(nóng)業(yè)用地、林地等人口較稀疏區(qū)域的河段(類型 III),主要污染來(lái)自面源.由于水溫隨季節(jié)變化明顯,且作為重要變量參與WASP模型中的多個(gè)模塊的計(jì)算,因此溫度水平按季節(jié)進(jìn)行劃分.由于北運(yùn)河水量受閘壩影響明顯,因此為保證流量的典型性,選用置信水平為 90%的監(jiān)測(cè)值來(lái)代表其正常范圍.根據(jù)此3類河段的各指標(biāo)監(jiān)測(cè)范圍,確定模型輸入變量的上下界,其中受到季節(jié)性影響較大的溫度及流量數(shù)據(jù)范圍見(jiàn)表2.

表2 研究區(qū)域季節(jié)性指標(biāo)數(shù)據(jù)范圍Table 2 Range of the seasonal indexes

2.2 模型不確定性分析

本研究將模型的不確定性分為整體與子模塊兩方面來(lái)進(jìn)行討論.模型的整體輸出不確定性主要體現(xiàn)參數(shù)變化給模型輸出指標(biāo)帶來(lái)的影響;子模型的不確定性可以體現(xiàn)整體輸出不確定性在各子模塊的分量.

2.2.1 模型整體輸出的不確定性 蒙特卡羅模擬中,模型參數(shù)可根據(jù)實(shí)際數(shù)據(jù)所展現(xiàn)的特征為依據(jù)來(lái)確定服從何種分布,但對(duì)于無(wú)法進(jìn)行參數(shù)界定概率分布的時(shí)候,一般情況可假設(shè)為均勻分布[37].本研究將這些參數(shù)設(shè)定為均勻分布.模擬所形成的條帶狀軌跡如圖 5所示.其中CBOD、NH3-N、NON等指標(biāo)擬合趨勢(shì)較好,而 DO 擬合效果不佳,這可能與模型輸入數(shù)據(jù)步長(zhǎng)有關(guān),由于DO在一天之內(nèi)的變化幅度也較其他指標(biāo)大,以步長(zhǎng)單位為 1d的輸入數(shù)據(jù)無(wú)法滿足其模擬精度.

圖3 2009年4月～2010年12月該河段內(nèi)年蒙特卡羅模擬下的模型的整體輸出Fig.3 Monte Carol simulation of the different indicators of WASP from April, 2009to December, 2010at Yulin Zhuang Segment觀測(cè)值為除枯水期以外的月監(jiān)測(cè)數(shù)據(jù)

以受閘壩干擾較小的榆林莊橋斷面為例,DO的模擬中,秋、冬的不確定性小于春、夏兩季;CBOD模擬中,夏季的不確定性最低;NON的模擬中冬季不確定性最低,其余季節(jié)不確定性相差不大;NH3-N模擬的不確定性四季相差不多.可見(jiàn),對(duì)于同一個(gè)參數(shù)范圍,由于輸入條件的隨時(shí)間的變化,參數(shù)對(duì)模型所造成的不確定性也可能隨變量水平的不同有著較大差異.

圖4 蒙特卡羅模擬下的模型的各子模塊的輸出Fig.4 Sub-models’ respond under Monte Carol simulation of WASP

2.2.2 子模塊不確定性 由于WASP模型的最終輸出是由若干個(gè)子模塊共同決定的,其中DO的模擬包括因水中CBOD的氧化所消耗的DO、藻類呼吸所消耗的 DO、因氨氮的硝化作用所消耗的DO、沉積物氧化分解所消耗的 DO、浮游植物生長(zhǎng)過(guò)程中因光合作用產(chǎn)生的DO、大氣復(fù)氧所增加的DO;CBOD的模擬包括CBOD氧化所消耗的量、反硝化菌因消耗碳源所減少的 CBOD、因沉淀作用所減少的CBOD、藻類死亡所產(chǎn)生的CBOD;硝態(tài)氮的模擬包括因反硝化作用所減少的硝態(tài)氮、浮游植物生長(zhǎng)所消耗的硝態(tài)氮,因氨氮硝化所增加的硝態(tài)氮;氨氮的模擬包括因礦化作用所減少的氨氮、因硝化作用所減少的氨氮、浮游植物生長(zhǎng)所消耗的氨氮、藻類死亡所產(chǎn)生的氨氮.根據(jù)蒙特卡羅模擬的結(jié)果,各個(gè)模塊的不確定性如圖4所示.

可以看出,對(duì)于DO的模擬,不確定性主要來(lái)自于沉積物、硝化作用等子模塊;CBOD的模擬方面,不確定性主要來(lái)自于沉淀作用;NO3-N模擬方面,不確定性主要是硝化作用與浮游植物生長(zhǎng); NH3-N的模擬方面,不確定性主要來(lái)自于浮游植物的死亡、生長(zhǎng)及硝化作用.然而,盡管NO3-N、NH3-N的模擬方面均有多個(gè)模塊存在較大的不確定性,但是最終疊加的結(jié)果的不確定性卻并不大.這說(shuō)明WASP模型當(dāng)中存在明顯的“異參同效”作用,各個(gè)子模塊此消彼長(zhǎng),僅靠傳統(tǒng)方法進(jìn)行全局尋優(yōu)率定,并不能體現(xiàn)水體各組分之間相互轉(zhuǎn)化的過(guò)程,還需要結(jié)合其他手段對(duì)關(guān)鍵參數(shù)進(jìn)行率定.

另外,由圖4可以看出,同一個(gè)子模塊,在不同的時(shí)期的不確定性也會(huì)發(fā)生很大的變化,這是由于輸入條件的改變所造成的.

2.3 模型參數(shù)敏感性分析

圖5 各指標(biāo)模擬參數(shù)敏感性指數(shù)Fig.5 Sobol’s index of the parameters in the simulation of each indicatorS為一階敏感性指數(shù),ST為總敏感性指數(shù),ST-S即為由參數(shù)之間交互、協(xié)同作用產(chǎn)生的影響

2.3.1 WASP模型主要敏感參數(shù)的確定 基于Sobol方法,利用輸出變化量以及納什系數(shù)(NSE)、平均相對(duì)誤差(MSE)兩個(gè)統(tǒng)計(jì)量作為似然函數(shù)進(jìn)行敏感性的度量.結(jié)果如圖5所示.

圖 5展示了以不同度量方法所得到參數(shù)總敏感性指數(shù)與一階敏感性指數(shù).總體上來(lái)講 3種方法結(jié)果具有一致性.這是因?yàn)閷?duì)于模型輸出有著顯著影響的參數(shù),其對(duì)模型的精度所產(chǎn)生的影響也相對(duì)較大.但是在進(jìn)行似然函數(shù)計(jì)算時(shí),各個(gè)參數(shù)相當(dāng)于在模型計(jì)算的基礎(chǔ)上又進(jìn)行了一次函數(shù)作用,因此參數(shù)相互之間的耦合作用就會(huì)進(jìn)一步增強(qiáng).如圖2所示,當(dāng)對(duì)于敏感性的度量選用似然函數(shù)時(shí),參數(shù)總敏感性一般要比一階敏感性大得多;而選用模型輸出變化幅度時(shí),參數(shù)的總敏感性與一階敏感性之間的差異較小.這就說(shuō)明,盡可能選用形式較為簡(jiǎn)單似然函數(shù),有助于在研究參數(shù)敏感性問(wèn)題中,減少因人為在使用似然函數(shù)時(shí)所造成的干擾.

綜合以上結(jié)果,得到WASP模型中對(duì)結(jié)果產(chǎn)生主要作用的參數(shù)(表3).

表3 WASP模型主要參數(shù)Table 3 The main sensitive parameters of WASP simulation at Tongzhou part of the North Canal

由表3可見(jiàn),除SOD、fD5、fon外,其他對(duì)最終結(jié)果產(chǎn)生重要影響的參數(shù),均參與多個(gè)模塊計(jì)算.盡管這些參數(shù)在各個(gè)模塊中的重要性不一致,但仍可能會(huì)與其他參數(shù)產(chǎn)生較強(qiáng)相互作用,如參數(shù)6、7、30,都在其相應(yīng)的模塊中出現(xiàn)了總敏感性指數(shù)遠(yuǎn)大于一階敏感性指數(shù)的現(xiàn)象.由于這類參數(shù)敏感性會(huì)因其他參數(shù)的取值不同而發(fā)生較大的變化.通過(guò)一次變化一個(gè)參數(shù)而固定其他參數(shù)進(jìn)行的局部敏感性分析就無(wú)法得到準(zhǔn)確的結(jié)果.因此“一次一個(gè)變量”的敏感性評(píng)估方法對(duì)于這類參數(shù)來(lái)說(shuō)并不適用.

對(duì)比利用蒙特卡羅模擬結(jié)果可以看出,敏感參數(shù)的所屬子模塊,其不確定性也較高:DO模塊中,參數(shù)SOD決定著沉積物耗氧量, E12決定著硝化作用的強(qiáng)弱,在DO的模擬中,這兩個(gè)參數(shù)的取值非常關(guān)鍵;CBOD模塊中,參數(shù)fD5能夠很大程度地決定模擬過(guò)程中沉淀量的多少,是模擬CBOD濃度的關(guān)鍵;硝態(tài)氮中的多個(gè)敏感參數(shù)也都集中在浮游植物生長(zhǎng)與硝化過(guò)程的子模塊中;氨氮模擬中的敏感參數(shù)也均集中在浮游植物的生長(zhǎng)、死亡、硝化過(guò)程等子模塊中.

但是由于邊界條件的改變會(huì)造成參數(shù)敏感性也隨之發(fā)生變化,所以盡管敏感參數(shù)在整個(gè)的模擬過(guò)程中呈主導(dǎo)地位,但在一定條件下的特殊時(shí)期,可能會(huì)發(fā)生變化.

2.3.2 輸入條件變化對(duì)參數(shù)敏感性的影響 為討論季節(jié)變化與空間變化所引起輸入條件的改變對(duì)參數(shù)敏感性的影響,本研究以四季的水質(zhì)水量均值數(shù)據(jù)作為模型的輸入數(shù)據(jù),用于體現(xiàn)四季的變化;以空間聚類分析結(jié)果的各類區(qū)域均值作為模型輸入,用于體現(xiàn)空間變化.

(1) 季節(jié)變化對(duì)參數(shù)敏感性的影響:我國(guó)北方城市的河流季節(jié)性較強(qiáng),往往在不同的季節(jié)水量不同,并同時(shí)呈現(xiàn)不同的水質(zhì)特征.季節(jié)影響著輸入變量的同時(shí),也會(huì)影響參數(shù)在一定時(shí)間范圍內(nèi)所呈現(xiàn)的敏感性.以榆林莊橋斷面的監(jiān)測(cè)數(shù)據(jù)為例,通過(guò)對(duì)四季的數(shù)據(jù)變化,分別進(jìn)行參數(shù)的敏感性計(jì)算,結(jié)果如圖6所示.

由圖 6可知,即使對(duì)于相同的河段來(lái)說(shuō),輸入條件在季節(jié)上的變化也會(huì)引起同一個(gè)參數(shù)敏感性的顯著改變.其中春、夏的敏感性分布情況較為接近,而秋、冬敏感性分布較為接近.這就意味著在參數(shù)率定方面,如果依照季節(jié)的不同,對(duì)參數(shù)進(jìn)行有針對(duì)性的率定,可以顯著改善率定工作的效率.

圖6 不同季節(jié)條件下的各參數(shù)敏感性指數(shù)Fig.6 Sobol indices of the parameters in different seasons

對(duì)比圖5與圖6可以發(fā)現(xiàn),在CBOD的模擬中,盡管參數(shù) fD5(16號(hào)參數(shù))在整體上是不確定性的主要來(lái)源(尤其是春、夏兩季),但是在秋、冬兩季的模擬中,其他參數(shù)對(duì) CBOD的模擬影響程度有所加大.因此在確定該參數(shù) fD5之后,在秋冬兩季的模擬中繼續(xù)對(duì)其他參數(shù)的率定,可以作為對(duì)該時(shí)間段模擬輸出值的微調(diào),以提高模擬精度.

區(qū)分各個(gè)時(shí)間段內(nèi)影響水質(zhì)模擬的主要參數(shù),可以為改善局部的模擬效果提供一定得理論依據(jù).圖 3中所示,在該河段內(nèi),春、夏時(shí)段內(nèi)的CBOD模擬,夏、秋、冬時(shí)段內(nèi)的DO模擬,秋、冬時(shí)段內(nèi)的NH3、NO3模擬都只有一個(gè)主要參數(shù),通過(guò)率定相應(yīng)的主要參數(shù),就可以顯著縮小模擬的不確定性.

(2) 空間變化對(duì)參數(shù)敏感性的影響 根據(jù)聚類分析得到的 3種不同類型的河段進(jìn)行參數(shù)的敏感性計(jì)算(圖7).

圖7 不同區(qū)域類型條件下的各參數(shù)敏感性指數(shù)Fig.7 Sobol indicesof the parameters in different area of the North Canal

由圖7可見(jiàn),在水質(zhì)特征不同的河段上,參數(shù)敏感性大小存在一些變化:其中變化最大的地方存在于NO3-N模擬中,參數(shù)E12、KNIT(2、3號(hào)參數(shù))的敏感性在第3類區(qū)域中明顯比其他兩類區(qū)域要低,這是因?yàn)榈?類區(qū)域當(dāng)中NO3-N明顯比其他兩個(gè)區(qū)域的含量低所造成的.

但從整體上看,北運(yùn)河通州段內(nèi)的空間變化一般并不影響參數(shù)的主導(dǎo)地位.這有可能是由于北運(yùn)河整體水質(zhì)較差,雖然根據(jù)聚類分析能夠依據(jù)水質(zhì)特征區(qū)分出不同類型的河段,但各河段均屬重污染、缺水型的納污河流,對(duì)于這類水體的模擬來(lái)說(shuō),其參數(shù)的敏感性基本保持統(tǒng)一,也就是說(shuō)在參數(shù)率定過(guò)程中,其率定的優(yōu)先級(jí)隨著空間的改變不會(huì)發(fā)生變化.

3 結(jié)論

3.1 對(duì)比總敏感性與一階敏感性指數(shù)之間的差異,發(fā)現(xiàn)WASP模型當(dāng)中的若干參數(shù)耦合性較強(qiáng),不宜使用一次一個(gè)變量法對(duì)這些變量進(jìn)行敏感性分析.

3.2 WASP的參數(shù)敏感性指數(shù)會(huì)隨模擬對(duì)象的時(shí)、空變化而發(fā)生改變;其中時(shí)間(季節(jié))變化的影響甚至可能會(huì)導(dǎo)致敏感參數(shù)數(shù)量上的變化,而空間變化的影響力不大.

3.3 WASP模型的不確定性在各個(gè)季節(jié)上也有一定差異:DO的模擬中,秋、冬的不確定性小于春、夏兩季;CBOD模擬中,夏季的不確定性最低;NON的模擬中冬季不確定性最低,其余季節(jié)不確定性相差不大;NH3-N模擬的不確定性四季相差不多.

3.4 WASP模型當(dāng)中存在明顯的“異參同效”作用,各個(gè)子模塊此消彼長(zhǎng),僅靠傳統(tǒng)方法進(jìn)行全局尋優(yōu)率定,并不能體現(xiàn)水體各組分之間相互轉(zhuǎn)化的過(guò)程,還需要結(jié)合實(shí)驗(yàn)對(duì)關(guān)鍵參數(shù)進(jìn)行率定.3.5 對(duì)于在一定時(shí)空條件下敏感性會(huì)突增的參數(shù),應(yīng)當(dāng)在該條件下另行率定.

[1] Vivian P, Roberto J C. Qual2E model for the Corumbata′? River[J]. Ecological Modelling, 2006,198:269—275.

[2] Mehmet Y, Erdal K, Ridvan B. Simulation of river streams:Comparison of a new technique with QUAL2E [J]. Mathematical and Computer Modelling, 2007,46:292—305.

[3] Fan C H, Ko C H, Wang W S. An innovative modeling approach using Qual2K and HEC-RAS integration to assess the impact of tidal effect on River Water quality simulation [J]. Journal of Environmental Management, 2009,90:1824—1832.

[4] Lin C E, Chen C T, Kao C M, et al. Development of the sediment and water quality management strategies for the Salt-water River,Taiwan [J]. Marine Pollution Bulletin, 2011,63:528—534.

[5] Thorsen M, Refsgaard J C, Hansen S, et al. Assessment of uncertainty in simulation of nitrate leaching to aquifers at catchment scale [J]. Journal of Hydrology, 2001,242:210-227.

[6] Karl-Erich L, Katrin F, Martina B. Structural uncertainty in a river water quality modeling system [J]. Ecological Modelling,2007,204:289—300.

[7] Fabrizio B, Benedicte L, Eva L G, et al. Estimation of sampling uncertainty in lake-water monitoringin a collaborative field trial[J]. Trends in Analytical Chemistry, 2012,36:176-184.

[8] Keith B, Jim F. Equinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology [J]. Journal of Hydrology,2001,249:11-29.

[9] Reichert P, Omlin M. On the usefulness of over parameterized ecological models [J]. Ecological Modelling, 1997,95:289-299.

[10] Lindenschmidt K E. The effect of complexity on parameter sensitivity and model uncertainty in river water quality modeling[J]. Ecological Modelling, 2006,190:72-86.

[11] 王建平,程聲通,賈海峰.基于MCMC法的水質(zhì)模型參數(shù)不確定性研究 [J]. 環(huán)境科學(xué), 2006,27(1):24-30.

[12] Lindim C, Pinho J L, Vieira J M P. Analysis of spatial and temporal patterns in a large reservoir using water quality and hydrodynamic modeling [J]. Ecological Modelling, 2011,222:2485-2494.

[13] Cea L, Bermúdez M, Puertas J. Uncertainty and sensitivity analysis of a depth-averaged water quality model [J].Environmental Modelling and Software, 2011,26:1526-1539.

[14] Saltelli A, Annonis P. How to avoid a perfunctory sensitivity analysis [J]. Environmental Modelling and Software, 2010,25:1508-1517.

[15] Saltelli A, Tarantola S, Chan K P S. A quantitative model independent method for global sensitivity analysis of model output [J]. Technometrics, 1999,41(1):39-56.

[16] 張 巍,鄭 一,王學(xué)軍.水環(huán)境非點(diǎn)源污染的不確定性及分析方法 [J]. 農(nóng)業(yè)環(huán)境科學(xué)學(xué)報(bào), 2008,27(4):1290-1296.

[17] Griensven A, Meixner T, Grunwald S, et al. A global sensitivity analysis tool for the parameters of multi-variable catchment models [J]. Journal of Hydrology, 2006,324:10-23.

[18] 蔡 毅,邢 巖,胡 丹.敏感性分析綜述 [J]. 北京師范大學(xué)學(xué)報(bào)(自然科學(xué)版), 2008,44(1):9-16.

[19] Sobol' I M. Sensitivity estimates for non-linear mathematical models [J]. Mathematical Modelling and Computational Experiment, 1993,4(1):407-414.

[20] Florian P, Keith B, Marco R, et al. Multi-method global sensitivity analysis of flood inundation models [J]. Advances in Water Resources , 2008,31:1-14.

[21] Dimov I, Georgieva R, IvansovskaS, et al. Studying the sensitivity of pollutants’ concentrations caused by variations of chemical rates [J]. Journal of Computational and Applied Mathematics, 2010,235:391-402.

[22] Yang J. Convergence and uncertainty analyses in Monte-Carlo based sensitivity analysis [J]. Environmetal Modelling and Software, 2011,26:444-457.

[23] Thierry A Mar, Stefano T. Variance-based sensitivity indices for models with dependent inputs [J]. Reliability Engineering and System Safety, 2012,107:115—121.

[24] Nossent J P Elsen. Sobol’ sensitivity analysis of a complex environmental model [J]. Environmental Modelling and Software,2011,26(12):1515-1525.

[25] 李 睿.Sobol靈敏度分析方法在結(jié)構(gòu)動(dòng)態(tài)特性分析中的應(yīng)用研究 [D]. 長(zhǎng)沙:湖南大學(xué), 2003.

[26] 王綱勝,夏 軍,陳軍鋒.模型多參數(shù)靈敏度與不確定性分析 [J].地理研究, 2010,29(2):263-270.

[27] Tim A Wool, Robert B Ambrose, James L Martin, et al. Water quality analysis simulation program (WASP) Draft: User’s manual [M]. Atlanta: US Environmental Protection Agency, MS.Tetre. Tech., 2001.

[28] 王旭東,劉素玲,張樹(shù)深,等.白洋淀水域 WASP富營(yíng)養(yǎng)化模型改進(jìn)研究 [J]. 環(huán)境科學(xué)與技術(shù), 2009,32(10):19-24.

[29] Arhonditsis G B, Brett M T. Eutrophication model for Lake Washington (USA) Part I. Model description and sensitivity analysis [J]. Ecological Modelling, 2005,187:140-178.

[30] 路成剛.基于WASP7.3的南四湖水質(zhì)模擬分析研究 [D]. 青島:青島理工大學(xué), 2010.

[31] 史鐵錘,王飛兒,方曉波.基于 WASP的湖州市環(huán)太湖河網(wǎng)區(qū)水質(zhì)管理模式 [J]. 環(huán)境科學(xué)學(xué)報(bào), 2010,30(3):631-640.

[32] 于 洋.北運(yùn)河水體中氨氮的氧化過(guò)程及微生物響應(yīng)特征 [D].北京:首都師范大學(xué), 2012.

[33] 方曉波,駱林平,李 松,等.錢(qián)塘江蘭溪段地表水質(zhì)季節(jié)變化特征及源解析 [J]. 環(huán)境科學(xué)學(xué)報(bào), 2013,33(7):1980-1988.

[34] Whitehead P G. Steady state and dynamic modeling of nitrogen in the River Kennet: impacts of Land use change since the 1930s [J].The Science of the Total Environment, 2002,282:417-434.

[35] Vaze J, FrancisH S Chiew. Experimental study of pollutant accumulation on an urban road surface [J]. Urban Water,2002,4:379-389.

[36] 孫金華,曹曉峰,黃 藝.滇池流域土地利用對(duì)入湖河流水質(zhì)的影響 [J]. 中國(guó)環(huán)境科學(xué), 2011,31(12):2052-2057.

[37] 鄒 銳,朱 翔,賀 彬,等.基于非線性響應(yīng)函數(shù)和蒙特卡洛模擬的滇池流域污染負(fù)荷削減情景分析 [J]. 環(huán)境科學(xué)學(xué)報(bào),2011,31(10):2312-2318.

Uncertainty analysis of WASP based on global sensitivity analysis method.

ZHANG Zhi-ming1,2, WANG Xiao-yan2*,LI Ming-tao2
(1.Beijing Climate Change Response Research and Education Center, Beijing University of Civil Engineering and Architecture, Beijing 100044, China；2.College of Resources, Environment and Tourism, Capital Normal University, Beijing 100048, China). China Environmental Science, 2014,34(5)：1336～1346

To quantitatively evaluate parameters influence and uncertainty of model, WASP model was reorganized by Simulink in this paper. Sensitivity of the parameters related to model output of DO, CBOD, NH3-N and NO3--N was studied based on Sobol method. In particular, sensitivity changes with temporal and spatial variations of input were discussed. Global sensitivity analysis of Sobol method can identify the most sensitive parameters of the process simulation; The significant difference between the total sensitivity and first-order sensitivity index showed that some parameters (e.g. GP1, PNH3, anc, etc.) of the WASP model were strongly coupling, and“One variable at a time”method was inappropriate to evaluate the sensitivity of these parameters; Sensitivity index of WASP parameter changed with the temporal and spatial variation of boundary conditions, even for the same river; WASP model showed an obvious equifinality for different parameters, so the traditional methods such as the global optimization calibration, failed to simulate the mechanism process. Experiments were required as a verification for calibration of key parameters.

uncertainty analysis；sensitivity；Sobol method；WASP model；boundary conditions

X703

A

1000-6923(2014)05-1336-11

2013-09-10

國(guó)家自然科學(xué)基金項(xiàng)目(40971258,41271495);高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金聯(lián)合資助項(xiàng)目(20121108110006);北京市教委北京市應(yīng)對(duì)氣候變化研究基地(2014年)(市級(jí))專項(xiàng)(PXM2014_014210_000037)

* 責(zé)任作者, 教授, cxnwxy@sohu.com

張質(zhì)明(1984-),男,北京人,博士,主要從事水質(zhì)模擬、模型不確定性方面的研究.發(fā)表論文10余篇.