基于SWAT模型的赤水河流域面源污染研究

薛亞婷，孫文錦，鄒長(zhǎng)武，高姍，陳結(jié)平

(成都信息工程大學(xué) 資源環(huán)境學(xué)院，成都 610225)

0 引言

水資源是人類生產(chǎn)和生活所必需的重要自然資源。受人類活動(dòng)的影響，社會(huì)在不斷進(jìn)步的同時(shí)也為水環(huán)境帶來巨大壓力，其中面源污染逐漸成為水環(huán)境破壞的主要原因[1-2]。面源污染即非點(diǎn)源污染，主要來源包括施肥不當(dāng)、農(nóng)藥使用、畜牧養(yǎng)殖等[3]。面源污染是由降雨作用和徑流共同作用產(chǎn)生，成因復(fù)雜，過程繁復(fù)，不具有明顯規(guī)律，且范圍一般較廣。面源污染中農(nóng)業(yè)面源污染占主導(dǎo)位置。因此，通過調(diào)整農(nóng)業(yè)面源污染的問題進(jìn)而管理水環(huán)境的舉措具有重要意義和價(jià)值。赤水河流域由于工廠污水肆意排放、污水排放管網(wǎng)不合理等問題，存在著較為嚴(yán)重的面源污染問題。

面源污染研究方法主要包括野外監(jiān)測(cè)法、類型劃分法、綜合試驗(yàn)場(chǎng)法、人工模擬試驗(yàn)法等，但都存在或多或少的問題。如野外監(jiān)測(cè)法費(fèi)時(shí)費(fèi)力、效率低、周期長(zhǎng)，只能作為一種輔助手段；綜合試驗(yàn)場(chǎng)法雖然工作量小、花費(fèi)低，但是精度不高；類型劃分法雖然提高了精度，但是費(fèi)用及工作量也大大提升。而面源污染模型常常與GIS技術(shù)、實(shí)地監(jiān)測(cè)方法相結(jié)合，對(duì)流域中污染物負(fù)荷量進(jìn)行了量化模擬，逐漸成為了研究流域污染物時(shí)空分布、模擬污染物地表遷移的最常用方法。

中國(guó)正在努力通過將SWAT模型與其他應(yīng)用軟件集成進(jìn)行應(yīng)用，包括流量和沉積物產(chǎn)量、非點(diǎn)源污染、模型不確定性和模型組合的研究。SWAT模型被廣泛應(yīng)用于徑流模擬、面源污染計(jì)算、地表水質(zhì)模擬、泥沙及氮磷等營(yíng)養(yǎng)物質(zhì)的轉(zhuǎn)移方面[4-7]。目前，國(guó)內(nèi)外有學(xué)者利用SWAT模型對(duì)流域氮磷污染問題進(jìn)行了深入研究[7]，也有學(xué)者將SWAT模型與MODFLOW地下水模型進(jìn)行耦合，模擬潛水蒸發(fā)量對(duì)區(qū)域徑流的影響[8]，還有學(xué)者對(duì)面源污染負(fù)荷的定性及定量分析計(jì)算，進(jìn)而推算出與降水量的正比關(guān)系，或者通過將SWAT模型與多方模型進(jìn)行組合構(gòu)建大尺度的區(qū)域綜合模擬模型。SWAT模型可以很方便地提取土地利用類型、土壤類型、天氣、農(nóng)業(yè)管理等數(shù)據(jù)，以實(shí)現(xiàn)對(duì)流域的模擬，具有運(yùn)算效率高、省時(shí)省力、長(zhǎng)期模擬等優(yōu)點(diǎn)。本研究選擇沱江二級(jí)支流、絳溪河一級(jí)支流赤水河流域進(jìn)行面源污染研究分析，以SWAT模型為支撐，將其與GIS、RS等技術(shù)相結(jié)合，以DEM數(shù)據(jù)、土地利用、土壤類型及屬性、氣象數(shù)據(jù)等為基礎(chǔ)數(shù)據(jù)[9-10]，構(gòu)建合理的模型，分析污染成因，找尋污染預(yù)防和控制對(duì)策，為絳溪河流域乃至沱江流域水環(huán)境綜合管理提供依據(jù)。

1 材料與研究方法

1.1 研究區(qū)概況

圖 1 赤水河流域位置示意Figure 1 Location of Chishui River Basin

赤水河是絳溪河的支流之一，絳溪河位于沱江的北側(cè)，絳溪河的發(fā)源地是仁壽縣的龍泉山脈，其流向?yàn)槲髂?東北流向，流經(jīng)鎮(zhèn)金、三岔、草池三鎮(zhèn)。赤水河自北匯入絳溪河，由龍泉流向簡(jiǎn)陽方向，流經(jīng)石盤鎮(zhèn)、簡(jiǎn)陽、太平橋鎮(zhèn)，上游蜿蜒于山際，下游流經(jīng)淺丘村鎮(zhèn)，河岸多臺(tái)地。全長(zhǎng)約25 km，流域面積約60 km2，平均流量0.68 m3·s-1。

赤水河流域遍及村莊、工廠、農(nóng)家樂、牲畜養(yǎng)殖場(chǎng)。由于工廠肆意排放污水、牲畜業(yè)管理不嚴(yán)、污水處理不到位等原因，該流域一直存在相當(dāng)嚴(yán)重的面源污染問題，故而本研究對(duì)該流域污染成因進(jìn)行分析，從而為該區(qū)域水質(zhì)改善和治理提供幫助。

1.2 基礎(chǔ)資料

表 1 空間數(shù)據(jù)信息Table 1 Spatial data information

表 2 屬性數(shù)據(jù)信息Table 2 Attribute data information

圖 2 土地利用類型數(shù)據(jù)Figure 2 Land use type data

圖 3 土壤類型數(shù)據(jù)Figure 3 Soil type data

SWAT模型運(yùn)行所需的數(shù)據(jù)為數(shù)字高程模型、土地利用數(shù)據(jù)、土壤類型數(shù)據(jù)、土壤性質(zhì)數(shù)據(jù)、氣象數(shù)據(jù)、實(shí)測(cè)水文水質(zhì)數(shù)據(jù)，如表1及表2所示。其中土地利用的數(shù)據(jù)源為2018年4月2日的Landsat 8遙感影像數(shù)據(jù)，來源于地理空間數(shù)據(jù)云。其他屬性數(shù)據(jù)及空間數(shù)據(jù)也為2018年的相關(guān)數(shù)據(jù)。水質(zhì)數(shù)據(jù)來源于水污染監(jiān)測(cè)站監(jiān)測(cè)建設(shè)項(xiàng)目，徑流數(shù)據(jù)由實(shí)際測(cè)量得到，并且選擇使用SWAT模型中國(guó)大氣同化驅(qū)動(dòng)數(shù)據(jù)集(CMADS V1.0)所提供的氣象數(shù)據(jù)作為參考。

利用DEM、土地利用數(shù)據(jù)、土壤類型數(shù)據(jù)構(gòu)建空間數(shù)據(jù)庫。參照文獻(xiàn)[11-12]將土壤利用數(shù)據(jù)、土壤類型數(shù)據(jù)進(jìn)行分類，如圖2及圖3所示。其中，林地占最大比例76.73%，剩余的依次為城市用地10.45%、耕地8.34%、水域8.35%、荒地0.63%。流域內(nèi)土壤類型為紫泥土和灰粽土，分別占比67.57%和32.43%。

利用氣象數(shù)據(jù)、土壤性質(zhì)數(shù)據(jù)構(gòu)建屬性數(shù)據(jù)庫。由于研究區(qū)域較小，沒有合適的氣象站點(diǎn)數(shù)據(jù)，故選擇使用SWAT模型中國(guó)大氣同化驅(qū)動(dòng)數(shù)據(jù)集(CMADS V1.0)所提供的氣象數(shù)據(jù)作為參考，其數(shù)據(jù)格式可直接用于SWAT模型運(yùn)行[13]。

1.3 研究方法

1.3.1 SWAT模型

SWAT模型是美國(guó)農(nóng)業(yè)部研究中心開發(fā)的物理流域模型，可對(duì)流域進(jìn)行模擬，建立徑流、土地利用、土壤類型、污染物之間的重要聯(lián)系[14]。SWAT模型具有運(yùn)算效率高、能長(zhǎng)期模擬等優(yōu)良特點(diǎn)[15]?；贒EM構(gòu)建研究區(qū)河網(wǎng)，將赤水河流域劃分為25個(gè)子流域，設(shè)置土地利用的分類閾值為10%，土壤數(shù)據(jù)的分類閾值20%，最終整個(gè)研究區(qū)域被劃分為101個(gè)水文響應(yīng)單元(HRU)[16-17]，進(jìn)而對(duì)徑流、水質(zhì)參數(shù)率定及檢驗(yàn)，并對(duì)赤水河流域面源污染進(jìn)行時(shí)間和空間特征分析。

1.3.2 SWAT-CUP

SWAT-CUP是被廣泛應(yīng)用在對(duì)SWAT進(jìn)行參數(shù)敏感性分析方面的工具，它包含了SUFI2、GLUE、Para So1、MCMC等敏感性分析方法，SWAT-CUP將其與SWAT模型鏈接，對(duì)目標(biāo)函數(shù)進(jìn)行多元回歸分析，以此達(dá)到對(duì)模型參數(shù)敏感度分析的目的[18]。本研究選用SWAT-CUP提供的SUFI2算法以進(jìn)行全局參數(shù)敏感性分析[19]，在眾多參數(shù)中找到最為敏感的參數(shù)進(jìn)行率定和調(diào)整。

3.故障原因按系統(tǒng)排列可使文章清晰，例如：進(jìn)氣系統(tǒng)、燃油系統(tǒng)、電控系統(tǒng)、機(jī)械部分。我將作者總結(jié)的7條重新排列：

1.4 模型評(píng)價(jià)指標(biāo)

評(píng)價(jià)指標(biāo)選用R2(相關(guān)性系數(shù))和Ens(納什系數(shù))，R2反映模擬結(jié)果與實(shí)測(cè)數(shù)值的符合程度，Ens說明模擬結(jié)果的準(zhǔn)確程度[6]?；谝陨蟽蓚€(gè)評(píng)價(jià)指標(biāo)來評(píng)價(jià)SWAT模型在赤水河流域的適用性，取相關(guān)性系數(shù)R2>0.6，納什系數(shù)Ens>0.5作為評(píng)價(jià)標(biāo)準(zhǔn)。其中，R2可利用線性回歸方程求得，當(dāng)R2=1時(shí)，數(shù)據(jù)完全相關(guān)；當(dāng)R2<1時(shí)，值越接近0表示吻合程度越低。Ens為1時(shí)表示模擬值與實(shí)測(cè)值二者達(dá)到完全一致。

2 結(jié)果與分析

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

本次試驗(yàn)選取若干相關(guān)參數(shù)，得到表3所示的運(yùn)行結(jié)果，其中，t檢驗(yàn)值反映參數(shù)的敏感性程度，絕對(duì)值越大越敏感；P值代表參數(shù)敏感性的顯著程度，越接近0即越顯著。一般選取t的絕對(duì)值>1且P<0.5的參數(shù)進(jìn)行校正。通過多次迭代運(yùn)算，得到對(duì)赤水河流域的徑流校準(zhǔn)結(jié)果。

表 3 敏感性分析結(jié)果Table 3 Results of sensitivity analysis

2.2 徑流參數(shù)率定及驗(yàn)證

表 4 年尺度徑流率定結(jié)果 Table 4 Results of annual runoff rate-setting

表 5 月尺度徑流率定結(jié)果 Table 5 Results of monthly runoff rate-setting

選取2013年為模型預(yù)熱期，2014—2015年為校準(zhǔn)期，2016年為驗(yàn)證期，進(jìn)行模型在赤水河流域的適用性驗(yàn)證。對(duì)徑流進(jìn)行率定的思路是依次從年、月、日3種尺度對(duì)參數(shù)進(jìn)行調(diào)整[20]。從年尺度來說，2014—2015年實(shí)測(cè)與模擬徑流結(jié)果的相對(duì)誤差Re值<20%，即年尺度徑流模擬效果良好，結(jié)果如表4所示。從月尺度來說，率定期內(nèi)月徑流R2>0.6，Ens>0.5，模擬效果良好，結(jié)果如表5所示，其對(duì)比效果如圖4所示。從日尺度來說，提取2014—2015年日徑流量模擬量，結(jié)果如圖5所示。

完成徑流參數(shù)率定后，對(duì)模型進(jìn)行驗(yàn)證，判斷精度是否達(dá)到要求[7]。選取2016年實(shí)測(cè)數(shù)據(jù)以驗(yàn)證模擬結(jié)果，得到如表6所示評(píng)價(jià)結(jié)果，驗(yàn)證期內(nèi)月徑流R2>0.6，Ens>0.5，模擬效果良好，對(duì)比效果如圖6所示。提取日徑流量模擬結(jié)果，2016年驗(yàn)證內(nèi)模擬效果如圖7所示。

圖 4 2014—2015年月徑流量Figure 4 Monthly runoff from 2014 to 2015

圖 5 2014—2015年日徑流量Figure 5 Daily runoff from 2014 to 2015

圖 6 2016年月徑流量Figure 6 Monthly runoff in 2016

圖 7 2016年日徑流量Figure 7 Daily runoff in 2016

2.3 水質(zhì)參數(shù)率定及驗(yàn)證

表 6 月尺度徑流驗(yàn)證結(jié)果 Table 6 Monthly runoff verification results

表 7 月尺度污染負(fù)荷率定結(jié)果 Table 7 Results of monthly pollution load rate determination

表 8 月尺度污染負(fù)荷驗(yàn)證結(jié)果 Table 8 Monthly pollution load verification results

表 9 2014—2016年4—9月面源污染輸出總量 Table 8 Total output of non-point source pollution from April to September, 2014—2016

通過計(jì)算可得實(shí)測(cè)泥沙、總氮和總磷的月均負(fù)荷量。基于2015年實(shí)測(cè)水質(zhì)數(shù)據(jù)來率定模型。運(yùn)行得到2015年泥沙、總氮(TN)、總磷(TP)實(shí)測(cè)數(shù)據(jù)與年尺度模擬數(shù)值對(duì)比，相對(duì)誤差Re分別為2.8%、-10.8%、-4.5%，均小于20%，即年尺度污染負(fù)荷模擬效果良好，月均率定結(jié)果如表7所示。

輸出月均污染負(fù)荷模擬結(jié)果，得到率定期內(nèi)月均污染負(fù)荷的實(shí)測(cè)值與模擬值R2>0.6，模擬效果良好。完成參數(shù)率定后，基于2016年實(shí)測(cè)水質(zhì)數(shù)據(jù)來進(jìn)行水質(zhì)驗(yàn)證。運(yùn)行得到2016年泥沙、總氮(TN)、總磷(TP)實(shí)測(cè)數(shù)據(jù)與年尺度模擬數(shù)值對(duì)比，相對(duì)誤差Re均小于20%，即年尺度污染負(fù)荷模擬效果良好，月均驗(yàn)證結(jié)果良好(表8)。

2.4 面源污染負(fù)荷年內(nèi)分析

通過運(yùn)行校準(zhǔn)后的SWAT模型對(duì)赤水河流域2014—2016年的泥沙、總氮(TN)、總磷(TP)負(fù)荷進(jìn)行模擬，得到該流域2014—2016年的泥沙輸出總量、總氮輸出總量和總磷輸出總量[21]分別為：3 422.4 t、13 721 kg、11 383 kg。從年尺度的變化觀察得到泥沙、總氮、總磷與降水?dāng)?shù)值的年變化形式有很強(qiáng)的相似性，并且一年之中都是先增大再減小。降雨量集中在6—9月。2014年降水從一月起開始增加，7、8月達(dá)到最大降水量，分別為252.07 mm和242.09 mm；2015年8、9月達(dá)到最大降水量，分別為171.77 mm和217.56 mm；2016年6、7月達(dá)到最大降水量，分別為134.36 mm和192.95 mm。泥沙、總氮、總磷也有著年尺度的強(qiáng)相似性。總氮在2014年8月達(dá)到一整年內(nèi)最高9 026 kg,在2015年9月達(dá)到最高10 463 kg，在2016年8月達(dá)到最高14 478.31 kg?？偭自?014年8月達(dá)到一整年內(nèi)最高694 kg,在2015年9月達(dá)到最高1 631.022 kg，在2016年6月達(dá)到最高998 kg。泥沙在2014年8月達(dá)到一整年內(nèi)最高2 960.35 t,在2015年7月達(dá)到最高2 065.595 t，在2016年8月達(dá)到最高2 854.21 t，按月尺度輸出分布如圖8所示。

面源污染負(fù)荷大多在雨季輸出，即每年的4—9月，其污染輸出總量占全年比重如表9所示。

圖 8 2014、2015和2016年月尺度負(fù)荷輸出分布Figure 8 Monthly load output distribution in 2014, 2015 and 2016

2.5 面源污染負(fù)荷空間分析

根據(jù)模型運(yùn)行得到的輸出文件(.rch)，從中提取各個(gè)河段的面源污染輸出負(fù)荷量[22]，根據(jù)子流域空間分布結(jié)合對(duì)應(yīng)的面積，計(jì)算得到各子流域單位面積輸出負(fù)荷，空間分布如圖9所示。

3 結(jié)論與討論

本研究以赤水河流域?yàn)檠芯繉?duì)象，將SWAT模型與GIS、RS等技術(shù)相結(jié)合，構(gòu)建水質(zhì)評(píng)價(jià)模型，對(duì)赤水河流域內(nèi)的面源污染負(fù)荷進(jìn)行分析，揭示其時(shí)空分布特征，得出以下結(jié)論。

圖 9 2014—2016年泥沙、總氮(TN)、總磷(TP)負(fù)荷空間分布Figure 9 Spatial distribution of sediment, total nitrogen (TN) and total phosphorus (TP) loads from 2014 to 2016

(1)SWAT模型建立后的參數(shù)率定與校驗(yàn)中，率定期和驗(yàn)證期的徑流和泥沙、總氮(TN)、總磷(TP)模擬在年、月、日尺度的誤差控制在合理范圍內(nèi)，說明該模型在赤水河流域具有良好的適用性。

(2)由矯正后模型的模擬結(jié)果推知面源污染負(fù)荷大多在雨季輸出。流域上游的單位面積污染負(fù)荷較少，泥沙負(fù)荷大致在0.09到0.17 t·hm-2之間，總氮負(fù)荷大致在1.54到2.68 kg·hm-2之間，總磷負(fù)荷大致在0.09到0.13 kg·hm-2之間。中部穿過城鎮(zhèn)受人類活動(dòng)影響，污染負(fù)荷明顯高于其他地區(qū)。下游農(nóng)業(yè)生產(chǎn)活動(dòng)較多，污染負(fù)荷也相對(duì)較高。中下游污染嚴(yán)重處泥沙負(fù)荷達(dá)到1.46 t·hm-2以上，總氮負(fù)荷達(dá)到14.34 kg·hm-2以上，總磷負(fù)荷達(dá)到1.50 kg·hm-2以上。總氮、總磷的空間特征受農(nóng)業(yè)活動(dòng)影響明顯，耕地和城市用地所在地區(qū)污染比較嚴(yán)重，總氮負(fù)荷在6.32 kg·hm-2以上，總磷負(fù)荷在0.41 kg·hm-2以上。而泥沙受降雨、地形、土地利用和土壤類型等綜合因素的影響，在耕地所在地區(qū)，泥沙負(fù)荷多在0.42 t·hm-2以上，土壤類型多為灰棕沖積土。

(3)通過對(duì)結(jié)果進(jìn)行分析尋找面源污染成因，為尋找提高流域水環(huán)境質(zhì)量的方法提供思路，為其他地區(qū)農(nóng)業(yè)面源污染研究提供借鑒。

本研究對(duì)赤水河流域進(jìn)行了不同時(shí)間尺度的徑流、水質(zhì)模擬，結(jié)果表明SWAT模型在該小流域具有一定的適用性，徑流模擬效果及污染負(fù)荷模擬效果均良好。另外，由上述分析可得面源污染負(fù)荷大多在雨季輸出，即每年的4—9月，其污染輸出總量占全年比重如表9所示，從中推知全年約80%以上的污染負(fù)荷輸出集中在流域汛期。由此推斷降雨對(duì)面源污染輸出起到至關(guān)重要的作用，是其產(chǎn)生的直接推動(dòng)力[23-25]。最后，由于本研究的流域面積很小，所以為減少工作量有部分參數(shù)的賦值為主觀設(shè)定，基于臨近區(qū)域或更大尺度上的數(shù)據(jù)計(jì)算得到，對(duì)精度會(huì)有一定影響，若有更為準(zhǔn)確的數(shù)據(jù)可進(jìn)一步對(duì)內(nèi)容進(jìn)行完善。SWAT模型在運(yùn)行及校正過程中涉及到非常多參數(shù)，而采集到的數(shù)據(jù)僅為基礎(chǔ)數(shù)據(jù)，若想要得到更準(zhǔn)確的數(shù)據(jù)，可參考模型說明搜集其他相關(guān)數(shù)據(jù)進(jìn)行模擬及校正。