基于替代模型的非點(diǎn)源污染模擬不確定性分析

馬夢(mèng)蝶　李傳奇　崔佳偉　楊幸子　王德振

摘要：為探究非點(diǎn)源污染機(jī)理模型參數(shù)的不確定性對(duì)模型精度和模擬結(jié)果的影響，以廣利河流域?yàn)檠芯繀^(qū)，運(yùn)用小波神經(jīng)網(wǎng)絡(luò)建立研究區(qū)非點(diǎn)源污染模擬SWAT模型的替代模型，解決使用蒙特卡羅方法進(jìn)行不確定性分析時(shí)運(yùn)算量大的問(wèn)題，從而實(shí)現(xiàn)對(duì)非點(diǎn)源污染模擬模型參數(shù)的不確定性分析。研究結(jié)果表明：小波神經(jīng)網(wǎng)絡(luò)模型不僅能夠很好地代替SWAT非點(diǎn)源污染模擬模型，而且能大幅度減少蒙特卡羅方法的計(jì)算時(shí)間：研究區(qū)氨氮和總磷輸出結(jié)果在90%置信水平下的置信區(qū)間分別為4.02×10～ 25.85×10、9.01×10 -28.38×103 kg/a，相比總磷，氨氮輸出結(jié)果的離散程度更高，不確定性程度更大。

關(guān)鍵詞：不確定性分析;非點(diǎn)源污染模擬;蒙特卡羅方法;拉丁超立方抽樣方法;小波神經(jīng)網(wǎng)絡(luò);廣利河流域

中圖分類(lèi)號(hào)：X522

文獻(xiàn)標(biāo)志碼：A

doi： 10.3969/j.issn.1000-1379.2019.06.015

近年來(lái)，隨著點(diǎn)源污染問(wèn)題被有效控制，非點(diǎn)源污染引起的污染問(wèn)題備受關(guān)注。相比集中排放的點(diǎn)源污染，非點(diǎn)源污染具有隨機(jī)性、滯后性、模糊性和控制難度大等特點(diǎn)[1]。采用機(jī)理模型對(duì)非點(diǎn)源污染負(fù)荷進(jìn)行量化分析與評(píng)價(jià)是非點(diǎn)源污染治理與研究的重要手段[2]，但非點(diǎn)源污染模型本身就是對(duì)自然系統(tǒng)的概化，其模型結(jié)構(gòu)、輸人數(shù)據(jù)和參數(shù)在應(yīng)用時(shí)存在很多不確定性，影響模擬結(jié)果的準(zhǔn)確性。其中，模型參數(shù)的不確定性是非點(diǎn)源污染模擬不確定性的重要來(lái)源[3]。因此，有必要了解非點(diǎn)源污染模擬時(shí)參數(shù)的不確定性對(duì)模型精度和模擬結(jié)果的影響程度，客觀評(píng)價(jià)模型模擬結(jié)果的準(zhǔn)確性，從而為管理者制定有效的污染防治措施提供借鑒[4]。

蒙特卡羅不確定分析方法是目前使用較為廣泛的一種非點(diǎn)源污染模擬不確定性分析方法[5]。該方法不用考慮狀態(tài)函數(shù)是否線(xiàn)性，隨機(jī)變量是否符合正態(tài)分布等問(wèn)題，可運(yùn)用模擬結(jié)果的不確定性來(lái)反映模型參數(shù)的不確定性。Sohrabi等[6]在敏感性分析的基礎(chǔ)上，應(yīng)用拉丁超立方分層抽樣和蒙特卡羅模擬得到了MACRO模型參數(shù)的不確定性對(duì)徑流和非點(diǎn)源污染模擬結(jié)果的影響：余紅等[7]通過(guò)蒙特卡羅方法模擬評(píng)估了SWAT模型參數(shù)的不確定性對(duì)徑流、泥沙和營(yíng)養(yǎng)物輸出結(jié)果的影響：邢可霞[8]以滇池流域?yàn)槔褂妹商乜_方法模擬評(píng)估了HSPF模型非點(diǎn)源污染輸出的不確定性程度。在使用蒙特卡羅方法進(jìn)行不確定性分析時(shí)需要對(duì)模型進(jìn)行多次運(yùn)算，對(duì)于計(jì)算量較大的非點(diǎn)源污染模擬模型往往需要消耗大量時(shí)間。

本文以廣利河流域?yàn)檠芯繀^(qū)，運(yùn)用小波神經(jīng)網(wǎng)絡(luò)建立研究區(qū)非點(diǎn)源污染模擬模型的替代模型，在使用蒙特卡羅方法進(jìn)行不確定分析需要大量模型運(yùn)算時(shí)直接調(diào)用替代模型，減少計(jì)算負(fù)荷，實(shí)現(xiàn)對(duì)研究區(qū)非點(diǎn)源污染模擬模型參數(shù)的不確定性分析，并采用切比雪夫不等式和變異系數(shù)來(lái)描述模型輸出結(jié)果的不確定性程度。

1 資料與方法

1.1 研究區(qū)概況

廣利河位于山東省東營(yíng)市，全長(zhǎng)約60 km，流域面積510 km2。流域?qū)儆诘湫偷臏貛Т箨懶詺夂騾^(qū)，多年平均氣溫為12.8℃，年降水量為536.7 mm（多集中在6-9月），種植的主要作物為小麥、玉米、水稻和棉花，主要土地利用類(lèi)型為耕地、城鎮(zhèn)用地，主要土壤類(lèi)型為沖積土和鹽土。

1.2 數(shù)據(jù)獲取與處理

研究所需要的數(shù)據(jù)可以分為空間數(shù)據(jù)和屬性數(shù)據(jù)兩類(lèi)。空間數(shù)據(jù)中的數(shù)字高程數(shù)據(jù)為分別率為30 mX30 m的GDEM—DEM數(shù)據(jù)，該數(shù)據(jù)來(lái)源于地理空間數(shù)據(jù)云：土地利用數(shù)據(jù)選取寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心的中國(guó)地區(qū)土地覆蓋綜合數(shù)據(jù)集，數(shù)據(jù)格式為GRID.分別率為1 000 mxl 000 m;土壤類(lèi)型數(shù)據(jù)采用黑河計(jì)劃數(shù)據(jù)管理中心世界土壤庫(kù)（ HWSD）中的1：100萬(wàn)中國(guó)土壤數(shù)據(jù)集。使用ArcGIS工具在行政區(qū)域圖基礎(chǔ)上對(duì)三類(lèi)數(shù)據(jù)進(jìn)行裁剪、投影、掩膜處理，得到研究區(qū)的空間數(shù)據(jù)集。研究區(qū)屬性數(shù)據(jù)包括氣象、土地管理及水文水質(zhì)數(shù)據(jù)，其中：氣象數(shù)據(jù)采用CMADS數(shù)據(jù)集，時(shí)間為2008-2016年：土地管理數(shù)據(jù)主要包括農(nóng)作物種類(lèi)、耕作方式，灌溉、施肥等作物種植模式，該數(shù)據(jù)的獲取主要由區(qū)域調(diào)研與查閱研究區(qū)歷年統(tǒng)計(jì)年鑒獲取：研究區(qū)水文水質(zhì)數(shù)據(jù)由相關(guān)部門(mén)提供，包括2015-2016年廣利河明海閘斷面的實(shí)測(cè)月徑流量，以及氨氮和總磷含量。

1.3 研究方法

1.3.1 SWAT模型

SWAT模型是美國(guó)農(nóng)業(yè)部（USDA）農(nóng)業(yè)研究中心（ARS）歷經(jīng)30 a開(kāi)發(fā)的具有較強(qiáng)物理機(jī)制的半分布式水文模型，可用于模擬水文、泥沙、污染物等的遷移轉(zhuǎn)化過(guò)程，是目前為止國(guó)內(nèi)使用最多的機(jī)理模型[9-10]。本文選取SWAT模型對(duì)研究區(qū)內(nèi)非點(diǎn)源污染進(jìn)行模擬，運(yùn)用SWAT-CUP軟件對(duì)模型模擬結(jié)果進(jìn)行自動(dòng)率定、校準(zhǔn)和參數(shù)敏感性分析。其中：率定和校準(zhǔn)采用SWAT-CUP軟件中SUFI2算法，敏感性分析采用模型自帶的LH -OAT方法[11]，模擬效果采用確定系數(shù)R2和納什系數(shù)E來(lái)評(píng)價(jià)[12-14]。

1.3.2 蒙特卡羅方法

蒙特卡羅方法是基于隨機(jī)數(shù)的一種不確定性分析方法（蒙特卡羅模擬），也是生成隨機(jī)數(shù)的一種抽樣方法（蒙特卡羅抽樣）。應(yīng)用蒙特卡羅方法進(jìn)行模型參數(shù)的不確定性分析一般分為三個(gè)步驟：假定參數(shù)符合一定的概率分布，通過(guò)抽樣獲取模型參數(shù)樣本：將參數(shù)樣本代人模型依次計(jì)算：對(duì)模型輸出結(jié)果進(jìn)行統(tǒng)計(jì)分析，計(jì)算統(tǒng)計(jì)特征值以描述模型的不確定性程度[6]。在利用蒙特卡羅方法模擬時(shí)，合適的參數(shù)樣本可有效減少模型模擬次數(shù)，因此在參數(shù)抽樣時(shí)，分別采用蒙特卡羅隨機(jī)抽樣和拉丁超立方分層抽樣[15]兩種方法，并對(duì)比分析兩種抽樣結(jié)果的優(yōu)劣，以選取合適的參數(shù)樣本。

1.3.3 小波神經(jīng)網(wǎng)絡(luò)模型

小波神經(jīng)網(wǎng)絡(luò)（ WA-ANN）是基于小波變換思想構(gòu)成的一種新的神經(jīng)網(wǎng)絡(luò)模型，由輸入層、隱含層、輸出層構(gòu)成。它的原理是用小波函數(shù)代替?zhèn)鹘y(tǒng)神經(jīng)網(wǎng)絡(luò)隱含層節(jié)點(diǎn)中的激發(fā)函數(shù)（如Sigmod函數(shù)），用小波函數(shù)中伸縮因子和平移因子來(lái)確定傳統(tǒng)神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)中的權(quán)值及閾值16]。小波神經(jīng)網(wǎng)絡(luò)不僅充分繼承了小波變換良好的局部化性質(zhì)和神經(jīng)網(wǎng)絡(luò)的自學(xué)習(xí)功能，而且避免了神經(jīng)網(wǎng)絡(luò)等在結(jié)構(gòu)設(shè)計(jì)上的盲目性，因此具有更強(qiáng)的學(xué)習(xí)能力、更快的收斂速度和更高的精度17]。只要參數(shù)設(shè)置及訓(xùn)練合理，小波神經(jīng)網(wǎng)絡(luò)可以以一定的精度逼近任意非線(xiàn)性函數(shù)。小波神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)見(jiàn)圖1（其中：x、x2、…、x為輸入樣本;y1、y2、…、y為輸出樣本;φ1、φ2、…、φ為小波函數(shù)）。

在學(xué)習(xí)階段，通過(guò)樣本訓(xùn)練不斷調(diào)整神經(jīng)網(wǎng)絡(luò)中的參數(shù)和權(quán)值以及隱含層節(jié)點(diǎn)個(gè)數(shù)，直到誤差迭代至容許誤差范圍之內(nèi)，從而使小波神經(jīng)網(wǎng)絡(luò)模型具有識(shí)別能力，得到相應(yīng)的輸出數(shù)據(jù)集。基于小波神經(jīng)網(wǎng)絡(luò)的記憶識(shí)別功能，建立研究區(qū)非點(diǎn)源污染模擬模型的替代模型。

1.3.4 統(tǒng)計(jì)分析

將切比雪夫不等式和變異系數(shù)C作為統(tǒng)計(jì)指標(biāo)來(lái)反映模型輸出結(jié)果的不確定性程度，切比雪夫不等式可用于計(jì)算不同概率下污染負(fù)荷輸出結(jié)果的置信區(qū)間，變異系數(shù)C可反映模型污染負(fù)荷輸出結(jié)果的離散程度。切比雪夫不等式：

2 模擬與分析

2.1 模型率定與校核

運(yùn)用收集到的廣利河明海閘斷面2015-2016年實(shí)測(cè)月徑流、氨氮、總磷數(shù)據(jù)對(duì)SWAT模型進(jìn)行率定和校核，其中2015年為率定期，2016年為校核期。模擬值與實(shí)測(cè)值對(duì)比見(jiàn)圖2（其中徑流用流量表示），可以看出擬合效果良好。率定期和校核期徑流、氨氮、總磷的確定系數(shù)與納什系數(shù)均在0.6以上，符合SWAT模型模擬的要求，具體見(jiàn)表1。根據(jù)LH -OAT方法敏感性分析的結(jié)果，選取對(duì)氨氮和總磷輸出結(jié)果較為敏感的6個(gè)參數(shù)用于后續(xù)不確定性分析，各參數(shù)物理意義及取值范圍見(jiàn)表2。

2.2 抽樣方法對(duì)比

為了選取合適的抽樣樣本，在Matlab平臺(tái)上分別編寫(xiě)蒙特卡羅方法和拉丁超立方方法抽樣程序，以SCS徑流曲線(xiàn)系數(shù)為例，抽取20組樣本進(jìn)行對(duì)比，結(jié)果見(jiàn)圖3。由圖3可以看出，拉丁超立方分層抽樣樣本在數(shù)軸上分布更加均勻，不存在樣本重合的問(wèn)題，而蒙特卡羅抽樣樣本則具有一定隨機(jī)性，在數(shù)軸上的分布沒(méi)有規(guī)律，樣本值之間的距離也忽近忽遠(yuǎn)。因此可以得出拉丁超立方抽樣樣本均勻性更好，更具代表性，在后續(xù)不同的抽樣組合中均采用此方法。

2.3 參數(shù)不確定性分析

利用拉丁超立方分層抽樣對(duì)氨氮和總磷輸出結(jié)果較為敏感的6個(gè)參數(shù)進(jìn)行100、10、1 000次抽樣。抽樣時(shí)假設(shè)各個(gè)參數(shù)相互獨(dú)立且符合均勻分布，其中100次和10次的抽樣參數(shù)輸入SWAT模型中計(jì)算得到污染負(fù)荷輸出結(jié)果，分別用于訓(xùn)練和檢測(cè)小波神經(jīng)網(wǎng)絡(luò)模型：1000次抽樣參數(shù)輸入訓(xùn)練好的小波神經(jīng)網(wǎng)絡(luò)模型中計(jì)算得到污染負(fù)荷輸出結(jié)果，用于統(tǒng)計(jì)分析。

根據(jù)小波神經(jīng)網(wǎng)絡(luò)的原理，以BP神經(jīng)網(wǎng)絡(luò)為基本框架，在Matlab平臺(tái)編寫(xiě)小波神經(jīng)網(wǎng)絡(luò)程序，模型輸入層神經(jīng)元個(gè)數(shù)為6，輸出層神經(jīng)元個(gè)數(shù)為2。經(jīng)過(guò)樣本訓(xùn)練，隱含層神經(jīng)元個(gè)數(shù)為8，即模型的最終結(jié)構(gòu)為6：8：2。采用SWAT模型的10組參數(shù)對(duì)小波神經(jīng)網(wǎng)絡(luò)模型的精度進(jìn)行檢測(cè)，結(jié)果見(jiàn)表3。

由表3可知，小波神經(jīng)網(wǎng)絡(luò)模型和研究區(qū)非點(diǎn)源污染模擬SWAT模型氨氮輸出結(jié)果擬合的平均相對(duì)誤差為1.91%.總磷輸出結(jié)果擬合的平均相對(duì)誤差為1 .23%，誤差較小，說(shuō)明小波神經(jīng)網(wǎng)絡(luò)模型可以代替研究區(qū)的SWAT模型。

采用切比雪夫不等式對(duì)小波神經(jīng)網(wǎng)絡(luò)模型輸出的1 000組氨氮和總磷的計(jì)算結(jié)果進(jìn)行區(qū)間估計(jì)，得到不同概率下氨氮和總磷輸出結(jié)果的置信區(qū)間。通過(guò)計(jì)算變異系數(shù)C，反映各污染負(fù)荷輸出結(jié)果的離散程度，結(jié)果見(jiàn)表4、表5。由表5可知，氨氮輸出結(jié)果的C值為23.11%，總磷輸出結(jié)果的C值為16.36%，表明相比總磷，研究區(qū)氨氮輸出結(jié)果的離散程度更大，不確定性也更大。

3 結(jié)果及討論

在采用蒙特卡羅方法研究非點(diǎn)源污染機(jī)理模型參數(shù)的不確定性中引入替代模型，分析了研究區(qū)模型參數(shù)不確定性對(duì)污染物輸出結(jié)果的影響，并解決了應(yīng)用蒙特卡羅方法進(jìn)行不確定性分析時(shí)計(jì)算量大，耗時(shí)長(zhǎng)的問(wèn)題。將小波神經(jīng)網(wǎng)絡(luò)作為非點(diǎn)源污染機(jī)理模型的替代模型，與以往研究中采用的替代模型相比[19-20]，精度更高、模擬誤差更小，即小波神經(jīng)網(wǎng)絡(luò)更適合作為非點(diǎn)源污染模擬的替代模型。但是小波神經(jīng)網(wǎng)絡(luò)作為黑箱模型，其輸出結(jié)果不具有時(shí)間性和空間性，在應(yīng)用方面存在局限。一般來(lái)說(shuō)，替代模型多應(yīng)用于不同領(lǐng)域模型的優(yōu)化迭代求解過(guò)程或者不確定性研究，而不將其直接作為數(shù)值模擬模型。

從氨氮和總磷輸出結(jié)果的置信區(qū)間看，不同置信水平氨氮和總磷輸出結(jié)果的區(qū)間變化范圍比較大，說(shuō)明參數(shù)的不確定性對(duì)兩種污染物輸出結(jié)果的影響程度都比較大，如果忽略模型輸入?yún)?shù)不確定性對(duì)輸出結(jié)果的影響，將會(huì)增大污染治理的防控風(fēng)險(xiǎn)，因此在評(píng)價(jià)與制定污染治理措施時(shí)，應(yīng)加強(qiáng)模型的不確定性分析以提高研究成果的可靠性。從氨氮和總磷變異系數(shù)的計(jì)算結(jié)果看，研究區(qū)氨氮輸出結(jié)果的離散程度更高，不確定性更大，主要原因：農(nóng)業(yè)施肥是非點(diǎn)源污染的主要來(lái)源，而農(nóng)業(yè)氮肥的施用量遠(yuǎn)高于磷肥，在徑流產(chǎn)污的過(guò)程中不確定性程度更高：氨氮化學(xué)性質(zhì)不穩(wěn)定，在遷移過(guò)程中容易發(fā)生物理或化學(xué)變化，物質(zhì)的存在形式不確定性較大。

4 結(jié)語(yǔ)

基于SWAT模型和小波神經(jīng)網(wǎng)絡(luò)替代模型對(duì)廣利河流域進(jìn)行了非點(diǎn)源污染模擬及不確定性分析，結(jié)果表明：運(yùn)用SWAT模型對(duì)廣利河流域進(jìn)行非點(diǎn)源污染模擬，模型率定和校核時(shí)期的R2和納什系數(shù)均在0.6以上，表明SWAT模型可以應(yīng)用于廣利河流域：小波神經(jīng)網(wǎng)絡(luò)模型作為一種黑箱模型，不僅能夠很好地代替SWAT非點(diǎn)源污染模擬模型，而且能大幅度減少蒙特卡羅方法模擬的計(jì)算時(shí)間，提高運(yùn)算效率：不同置信水平氨氮和總磷輸出結(jié)果的區(qū)間變化范圍比較大，說(shuō)明SWAT模型參數(shù)的不確定性對(duì)污染輸出的結(jié)果影響較大，且從變異系數(shù)的計(jì)算結(jié)果來(lái)看，氨氮輸出結(jié)果的不確定性程度更高。

參考文獻(xiàn)：

[l]賀纏生，傅伯杰，陳利頂.非點(diǎn)源污染的管理及控制[J].環(huán)境科學(xué)，1998，19（5）：88-92，97.

[2] 夏軍，翟曉燕，張永勇.水環(huán)境非點(diǎn)源污染模型研究進(jìn)展[J].地理科學(xué)進(jìn)展，2012.31（7）：941-952.

[3]廖謙，沈珍瑤，農(nóng)業(yè)非點(diǎn)源污染模擬不確定性研究進(jìn)展[J].生態(tài)學(xué)雜志，2011，30（7）：1542-1550.

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

[5]鄭一，流域非點(diǎn)源污染模擬和管理的不確定性研究[J].中國(guó)農(nóng)學(xué)通報(bào)，2008，24（增刊）：101-105.

[6]

SOHRABI T M， SHIRMOHAMMADI A，MONTAS H.Un-certainty in Nonpoint Source Pollution Models andAssociated Risks[J]. Environmental Forensics， 2002，3（2）：179-189.

[7]余紅，沈珍瑤.大寧河流域非點(diǎn)源污染不確定性分析[J].北京師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2008，44（1）：86-91.

[8]邢可霞.流域非點(diǎn)源污染模擬及不確定性研究[D].北京：北京大學(xué)，2005：129.

[9] 劉莊，晁建穎，張麗，等.中國(guó)非點(diǎn)源污染負(fù)荷計(jì)算研究現(xiàn)狀與存在問(wèn)題[J].水科學(xué)進(jìn)展，2015，26（3）：432-442.

[10] 張楠，張占龐，周祖昊，等.基于SWAT的天津市非點(diǎn)源污染分布特性分析[J].人民黃河，2008，30（4）：60-61.

[11]左德鵬，徐宗學(xué).基于SWAT模型和SUFI-2算法的渭河流域月徑流分布式模擬[J].北京師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2012，48（5）：490-496.

[12] PANAGOPOULOSY，MAKROPOULOSC，BALTASE，et a1.SWAT Parameterization for the Identifcation of Critical DiffusePollution Source Areas Under Data Lmitations[J].EcologicalModelling， 2011， 222（19）： 3500-3512.

[13] MORIASI D N， ARNOLD J G，LIEW M W V， et al.Model Evaluation Cuidelines for Systematic Quantificationof Accuracy in Watershed Simulations[J].Transactions ofthe Asabe， 2007， 50（3）： 885-900.

[14] 張培培，李瓊，闞紅濤，等.基于SWAT模型的植草河道對(duì)非點(diǎn)源污染控制效果的模擬研究[J].農(nóng)業(yè)環(huán)境科學(xué)學(xué)報(bào)，2014，33（6）：1204-1209.

[15] 雷鵬，陳曉偉，張貴金，等.基于LHS-MC的堤防滲透破壞風(fēng)險(xiǎn)分析[J].人民黃河，2014，36（ 10）：45-47.

[16]侯澤宇，盧文喜，陳社明，基于小波神經(jīng)網(wǎng)絡(luò)方法的降水量預(yù)測(cè)研究[J].節(jié)水灌溉，2013（3）：31-34.

[17]蔡念，胡匡祜，李淑宇，等.小波神經(jīng)網(wǎng)絡(luò)及其應(yīng)用[J].中國(guó)體視學(xué)與圖像分析，2001，6（4）：239-245.

[18]

ZHANC Z，LU W X， ZHAO Y， et al.Development Tend-ency Analysis and Evaluation of the Water Ecological Car-rying Capacity in the Siping Area of Jilin Province in ChinaBased on System Dynamics and Analytic Hierarchy Process[J]. Ecological Modelling， 2014， 275： 9-21.

[19] 宋文博.伊通河流域非點(diǎn)源污染模擬及不確定分析[D].長(zhǎng)春：吉林大學(xué)，2016：71.

[20] 閆雪嫂，盧文喜，歐陽(yáng)琦.基于替代模型的非點(diǎn)源污染模擬不確定性分析：以石頭口門(mén)水庫(kù)匯水流域?yàn)槔齕J].中國(guó)環(huán)境科學(xué)，2017，37（8）：3011-3018.