松濤水庫上游流域土地利用類型變化對氮負(fù)荷影響

雷湘齡 葉長青 朱麗蓉 薛楊

（1.海南省林業(yè)科學(xué)研究院 海南海口 571100；2.海南大學(xué)生態(tài)與環(huán)境學(xué)院 海南海口 570228；3.海南文昌森林生態(tài)系統(tǒng)國家定位觀測研究站 海南文昌 571300；4.海南大學(xué)旅游學(xué)院 海南海口 570228）

隨著各國對水污染問題的重視，點(diǎn)源污染已逐步得到有效的控制和治理，非點(diǎn)源污染成為水污染的主要來源[1]，北京密云水庫、安徽巢湖、天津于橋水庫、上海淀山湖、云南洱海等水域，非點(diǎn)源污染比率超過了點(diǎn)源污染[2]。土地利用類型的變化使土地下墊面發(fā)生改變，不適當(dāng)?shù)耐恋乩梅绞胶娃r(nóng)業(yè)管理措施在降雨-徑流驅(qū)動(dòng)因子的作用下加劇了土壤侵蝕，化肥中過量的氮、磷污染物附著于泥沙顆粒及可溶性氮磷污染物進(jìn)入水體，從而形成非點(diǎn)源污染[3-5]。

水文模型模擬是定量描述流域非點(diǎn)源污染過程的重要手段[6-7]。SWAT（Soil and Water Assessment Tool）模型被稱為是以農(nóng)業(yè)和森林為主的流域中較有前途的具有連續(xù)模擬能力的非點(diǎn)源污染模型。Anna等[8]使用SWAT模型對歐洲南部的多瑙河流域水質(zhì)進(jìn)行評估，發(fā)現(xiàn)流域內(nèi)主河道的總氮和總磷主要來源于農(nóng)作物種植施用的化肥；陳健等[9]和耿潤哲等[10]的研究結(jié)果也表明，農(nóng)業(yè)活動(dòng)所使用的化肥、農(nóng)藥等帶來的營養(yǎng)物流失是流域內(nèi)發(fā)生非點(diǎn)源污染的主要原因。大量的研究表明，土地利用變化會(huì)對非點(diǎn)源污染產(chǎn)生影響，Sushil等[11]在印度莫羅爾流域的研究說明，采用等高耕種和種植植被過濾帶能夠減少泥沙和土壤中營養(yǎng)物質(zhì)的流失；劉晶晶等[12]模擬湖北省清江流域宜都段2020年的土地利用情景，結(jié)果顯示，在耕地減少、林地增加、城鎮(zhèn)面積擴(kuò)大的情景下，總氮和總磷負(fù)荷分別減少23.77%和29.53%；秦耀民等[13]對陜西黑河流域不同土地利用情景進(jìn)行模擬，發(fā)現(xiàn)流域內(nèi)非點(diǎn)源污染負(fù)荷隨著林地面積的增加、耕地和未利用地的減少而減少；Lin等[14]在伊利湖西部流域的研究結(jié)果顯示，采用合理的農(nóng)業(yè)措施能夠降低流域內(nèi)地表水體中氮、磷等營養(yǎng)物的負(fù)荷。

近年來非點(diǎn)源污染成為當(dāng)下研究的熱點(diǎn)之一，國內(nèi)已對我國西北部、南部、東部典型流域的非點(diǎn)源污染開展研究，但鮮見關(guān)于熱帶地區(qū)河流的研究報(bào)道。海南省位于我國的最南端，屬于熱帶季風(fēng)氣候，松濤水庫是海南省最大的水庫。本研究以松濤水庫上游流域 1975—2014年的氣象數(shù)據(jù)和 1988、2002、2016年土地利用圖構(gòu)建SWAT模型，分析土地利用變化（砍伐熱帶雨林種植橡膠等經(jīng)濟(jì)林）對非點(diǎn)源污染氮負(fù)荷的影響，并識別造成非點(diǎn)源污染的“源-匯”區(qū)。

1 材料與方法

1.1 研究區(qū)概況

松濤水庫上游流域位于海南省白沙縣內(nèi)（圖1），流域內(nèi)山地多，森林資源豐富，平均氣溫高，日照長、雨量豐沛，以臺風(fēng)帶來的降雨為主；多年平均降雨量為 1 859 mm，年最大降雨量可達(dá)2 637.8 mm，存在明顯的干濕兩季，11月至次年4月為干季，降雨量占全年的15.9%，5—10月為濕季，占全年降雨量的84.1%。

圖1 研究區(qū)域圖

以流域內(nèi)福才水文站點(diǎn)的水質(zhì)數(shù)據(jù)進(jìn)行模擬研究，控制的集水面積為 492.7 km2。根據(jù) 2012年海南省統(tǒng)計(jì)年鑒，白沙縣一年施用氮肥8 688 t，耕地和水旱田面積共有268.38 km2，平均氮肥施用320 kg/hm2（折純氮）。流域內(nèi)有2個(gè)排污口，分別用于南開鄉(xiāng)衛(wèi)生院排污和生活用水排污，氨氮負(fù)荷分別為0.03和0.04 t一年（表1）；而根據(jù)福才站的實(shí)測數(shù)據(jù)，氨氮負(fù)荷每年約有5 t；由此可見，在該流域內(nèi)水質(zhì)污染以非點(diǎn)源為主。

表1 研究區(qū)內(nèi)點(diǎn)源排放統(tǒng)計(jì)

1.2 方法

1.2.1 模型簡介SWAT模型中，水文循環(huán)陸地階段基于水量平衡的方程為：

式中：SWt表示土壤最終含水量（mm）；SWo表示第i天的土壤初始含水量（mm）；t表示時(shí)間（d）；Rday表示第i天的降水量（mm）；Qsurf表示地表徑流總量，單位為mm/hm2；Ea表示第i天的蒸散發(fā)量（mm）；Wseep表示第i天離開土壤剖面底部的滲透水流和旁通水流的水量（mm）；Qgw表示第i天回歸流的水量（mm）。

修正的通用土壤流失方程（USLE）如下：

式中：sed表示某天的產(chǎn)沙量（t）；Qsurf表示地表徑流總量（mm/hm2）；qpeak表示洪峰流量（m3/s）；areahru表示水文響應(yīng)單元HRU的面積（hm2）；KUSLE表示方程中的土壤可侵蝕因子0.013 t·m2·h/(m2·t·cm)；CUSLE表示土壤覆蓋與農(nóng)業(yè)管理措施因子；PUSLE表示水土措施保持因子；LSUSLE表示地形因子；CFRG表示粗糙度因子。

1.2.2 數(shù)據(jù)來源構(gòu)建SWAT模型需要建立流域的空間數(shù)據(jù)庫和屬性數(shù)據(jù)庫以保證模型模擬的精確性，構(gòu)建SWAT模型所需數(shù)據(jù)來源見表2。

表2 構(gòu)建模型所需的數(shù)據(jù)來源

1.2.3 模型的構(gòu)建

1.2.3 .1 數(shù)字高程模型數(shù)據(jù)（DEM）采用SWAT模型對 DEM 進(jìn)行填洼，計(jì)算河流流向，對流累積柵格進(jìn)行預(yù)處理，通過設(shè)置上游匯水面積的閾值劃分子流域，共劃分16個(gè)子流域。

1.2.3 .2 土地利用圖使用 ENVI軟件的監(jiān)督分類方法，結(jié)合谷歌地球影像圖，采用人機(jī)交互的方式對30 m×30 m遙感影像進(jìn)行解譯，將研究區(qū)的土地利用類型分為七類：漿紙林、水域、橡膠林、耕地、居民區(qū)、果園、天然林。

1.2.3 .3 子流域和水文響應(yīng)單元?jiǎng)澐諷WAT模型屬于分布式的水文模型，通過離散化來體現(xiàn)流域的空間分布[15]。在子流域上疊加土地利用圖、土壤類型圖和坡度圖，進(jìn)一步劃分水文響應(yīng)單元。

1.2.4 參數(shù)的率定與驗(yàn)證使用 SWAT-CUP的SUFI-2算法對SWAT模擬輸出的數(shù)據(jù)進(jìn)行率定和驗(yàn)證。SWAT 模型中涉及徑流與營養(yǎng)物負(fù)荷的參數(shù)眾多，本研究通過查閱相關(guān)文獻(xiàn)和 SWAT-CUP的多次迭代確定參數(shù)，利用 SUFI-2 的全局敏感性分析方法評價(jià)參數(shù)的敏感性，敏感性越高表示該參數(shù)對徑流或營養(yǎng)物負(fù)荷的貢獻(xiàn)率越大，選用的徑流參數(shù)及其率定與驗(yàn)證見表3、4。

表3 SWAT模型徑流敏感性參數(shù)及其率定值

以Nash-Sutcliffe效率系數(shù)Ens和確定性系數(shù)R2為評價(jià)指標(biāo)，R2越趨近于 1，表明實(shí)測值與模擬值的擬合度越高[16]。當(dāng)Ens>0.75時(shí)，模擬結(jié)果良好；當(dāng)0.36≤Ens≤0.75時(shí)，結(jié)果令人滿意，說明模型適用于該流域；當(dāng)Ens<0.36，模擬效果較差，表明模型并不適用[17]。

表4 徑流參數(shù)的率定與驗(yàn)證

1.2.4 .1 徑流的率定和驗(yàn)證根據(jù) Mann-Kendall法對年徑流系數(shù)突變檢驗(yàn)的結(jié)果，選取 1977—1996年為模型率定期（1975—1976年設(shè)置為模型預(yù)熱期），1997—2014年為驗(yàn)證期。徑流參數(shù)率定與驗(yàn)證的R2和Ens均達(dá)到0.90以上，說明SWAT模型在該流域內(nèi)對徑流具有較好的模擬效果。

1.2.4 .2 氨氮的率定和驗(yàn)證福才水質(zhì)監(jiān)測數(shù)據(jù)為2012—2014年，設(shè)置2012年為率定期，2013—2014年為驗(yàn)證期。氨氮率定與驗(yàn)證的Ens均在0.54以上，說明SWAT模型適用于該流域營養(yǎng)物負(fù)荷的模擬，見表5、6。

表5 SWAT模型氨氮敏感性參數(shù)及其率定值

1.2.4 .3 總氮的率定和驗(yàn)證根據(jù)福才水文站點(diǎn)的相關(guān)監(jiān)測數(shù)據(jù)得到 2013年的總氮數(shù)據(jù)，于2014年分別對下游的南豐水文站點(diǎn)氨氮和總氮濃度數(shù)據(jù)進(jìn)行監(jiān)測，計(jì)算南豐站點(diǎn)氨氮占總氮的百分比，約為25%；以此比例估算福才站點(diǎn)2014年的總氮濃度，并對估算出的總氮數(shù)據(jù)進(jìn)行率定與驗(yàn)證（表7、圖2~4），以2013年為率定期，2014年為驗(yàn)證期。

圖2 1988、2002、2016年土地利用下徑流的率定與驗(yàn)證

表6 氨氮的率定與驗(yàn)證

表7 總氮的率定與驗(yàn)證

圖3 1988、2002、2016年土地利用下氨氮的率定與驗(yàn)證

圖4 1988、2002、2016年土地利用下總氮（TN）的率定與驗(yàn)證

2 結(jié)果與分析

2.1 土地利用類型的變化

研究區(qū)內(nèi)土地利用類型以天然林和經(jīng)濟(jì)林為主，1988—2016年漿紙林、橡膠林和耕地的面積有所增加，橡膠林面積占比由1988年的9.74%增加至 2016年的 11.30%，漿紙林由 5.19%增加至9.49%；而天然林面積逐漸減少，由 80.15%降至64.63%（表8），全年總氨空間分布情況見圖5。

圖5 1988、2002、2016年土地利用下的全年總氮空間分布圖

表8 土地利用類型及其面積比

2.2 總氮負(fù)荷的時(shí)空分布

降雨是導(dǎo)致地表土壤養(yǎng)分流失的外部動(dòng)力條件[18]。雨滴落下沖刷地表，使得土壤中的氮、磷等營養(yǎng)元素流失[19]，在降雨-徑流驅(qū)動(dòng)因子的作用下進(jìn)入水體，從而產(chǎn)生了非點(diǎn)源污染。松濤水庫上游流域 5—10月為濕季，降雨量占全年的84.1%，其總氮輸出約占全年 85%；干季降雨量占全年的 15.9%，總氮輸出占比約為 15%。有研究表明，山東省小清河流域6—9月的降雨量占全年的70%，在該時(shí)間段內(nèi)非點(diǎn)源污染負(fù)荷的輸出占比為全年的一半以上[20]；清水河非點(diǎn)源污染物總氮、總磷與月均降雨量、產(chǎn)水量的變化趨勢基本一致，當(dāng)降雨量達(dá)到峰值時(shí)，總氮、總磷輸出量也到達(dá)最大值[21]。上述研究均與本文的研究結(jié)論相符，因此在該流域內(nèi)非點(diǎn)源污染負(fù)荷以濕季輸出為主。

降雨產(chǎn)生徑流，而徑流在遷移的過程同時(shí)還受到土地利用方式的影響，因此不同的土地利用方式下非點(diǎn)源污染物的輸出存在較大的差異[22]。采用SWAT模型分別模擬1988、2002、2016年三期土地利用下總氮負(fù)荷的變化，結(jié)果表明，1988—2016年流域全年氮負(fù)荷量表現(xiàn)為顯著的增加趨勢；與 1988年相比，2002年流域內(nèi)全年總氮負(fù)荷增加了114.6%，2016年增加了171.65%，在此期間流域內(nèi)天然林面積減少，耕地、漿紙林、橡膠林的面積增加。經(jīng)濟(jì)林面積的增加意味著施肥量的增加，但作物對氮肥當(dāng)季的利用率僅為30%~35%[23]，剩余的營養(yǎng)元素附著于土壤顆粒，當(dāng)發(fā)生降雨時(shí)，隨徑流遷移入水體，從而造成了流域內(nèi)非點(diǎn)源污染氮負(fù)荷量的增加。

由總氮負(fù)荷的空間分布（圖5）可知，1988年土地利用下的總氮負(fù)荷量最大的是7和16號子流域，2002年是 3、5、7、8、16號子流域的總氮輸出量最大；相比于 1988年，2002年屬于總氮負(fù)荷最大梯度等級的子流域增加了3個(gè)，其中3和8號子流域的耕地面積增加，5號子流域的橡膠面積增加；16個(gè)子流域中總氮輸出增長幅度最大的是1號子流域，漲幅為 223.58%，其耕地面積增加。2016年3、5、6、7、8、16號子流域的總氮輸出量最大，相比于1988年，總氮負(fù)荷最大梯度等級的子流域增加了4個(gè)，增長幅度最大的是8號子流域，漲幅為407.66%，其漿紙林和耕地面積大幅度增加。總氮輸出量高的子流域內(nèi)均以耕地和經(jīng)濟(jì)林為主，1988—2016年間天然林面積減少，由人工開墾種植的耕地、經(jīng)濟(jì)林面積增加；土地下墊面的改變，施肥量的增加，造成氮負(fù)荷量的增加。

2.3 土地利用與總氮負(fù)荷的關(guān)系

土地利用類型的單位面積氮負(fù)荷量從高到低依次為：耕地>果園>橡膠林>漿紙林>天然林。不同土地利用類型的氮負(fù)荷量與植被的覆蓋度相關(guān)，土壤中氮元素的流失是地表徑流與土壤中氮素相互作用的結(jié)果，當(dāng)植被覆蓋度增加，土壤顆粒與地表徑流充分作用，使水流速度減慢，同時(shí)使得隨著徑流遷移的粗顆粒沉淀；反之，泥沙中細(xì)顆粒含量增加，所以植被覆蓋度的增加有利于泥沙中細(xì)顆粒的富集，從而導(dǎo)致泥沙中全氮含量增加[24]。在研究區(qū)內(nèi)，耕地的單位面積氮負(fù)荷量最大，達(dá)到7.70 kg/hm2，在同等施肥條件下，相比于果園、橡膠林和漿紙林，耕地的植被覆蓋率最低；其次是果園，而天然林在當(dāng)?shù)厣L多年，植被覆蓋率高，無人工干預(yù)不施用化肥，因此單位面積氮負(fù)荷最低，為1.76 kg/hm2，其土壤中的氮營養(yǎng)物質(zhì)主要來源于生物固氮和微生物[25]。

土地利用類型對全年總氮輸出的貢獻(xiàn)率由大到小依次為：天然林>橡膠林>漿紙林>耕地>果園。天然林的單位面積氮負(fù)荷量最小，但其面積占比最大，因此輸出氮負(fù)荷量的貢獻(xiàn)率最大。1988—2016年，不同土地利用類型的氮負(fù)荷貢獻(xiàn)率與其面積變化相關(guān)，耕地、橡膠林、漿紙林的面積呈增加的趨勢，氮負(fù)荷貢獻(xiàn)率隨之增加，如1988年的耕地面積為25.56 km2，2016年為46.75 km2，其貢獻(xiàn)率從3.37%漲到了15.25%（表9）；1988—2016年漿紙林面積增加了21.19 km2，氮負(fù)荷貢獻(xiàn)率增加了一倍；相反，天然林面積減少了76.43 km2，貢獻(xiàn)率從1988年的58.4%下降至2016年的34.33%。

對1988、2002、2016年研究區(qū)內(nèi)各子流域不同土地利用類型的面積與其輸出的總氮量進(jìn)行Pearson相關(guān)性分析，見表8。1988年天然林與總氮負(fù)荷呈極顯著相關(guān)；2002年果園與總氮負(fù)荷顯著相關(guān)，漿紙林、天然林與總氮負(fù)荷呈極顯著相關(guān)；2016年耕地、橡膠林與總氮負(fù)荷顯著相關(guān)，漿紙林、天然林與總氮負(fù)荷呈極顯著相關(guān)。天然林在流域內(nèi)面積占比最大，因此在三期土地利用下均與總氮負(fù)荷呈極顯著相關(guān)，但隨著面積的減小，其相關(guān)系數(shù)逐年降低；相反，耕地、橡膠林、漿紙林的面積隨年限有所增加，氮負(fù)荷量隨之增加，其相關(guān)系數(shù)逐漸增大，甚至在2016年達(dá)到顯著或極顯著相關(guān)。

將景觀生態(tài)學(xué)的“源-匯”理論應(yīng)用于非點(diǎn)源污染中，“源”是土壤養(yǎng)分流失的源頭，“匯”則是能夠接納遷移物質(zhì)，并對非點(diǎn)源污染具有滯緩作用的場所[26-27]。結(jié)合不同土地利用類型氮負(fù)荷的貢獻(xiàn)率（表9）和Pearson相關(guān)分析（表10），研究區(qū)內(nèi)耕地在1988和2002年與總氮負(fù)荷呈負(fù)相關(guān)，但隨著面積的增加，2016年達(dá)到顯著相關(guān)，其貢獻(xiàn)率也逐年增加；橡膠林、漿紙林也隨著面積的增加，相關(guān)系數(shù)和貢獻(xiàn)率增加，因此耕地、橡膠林和漿紙林對非點(diǎn)源污染氮負(fù)荷起到“源”的作用。天然林由于在流域內(nèi)面積占比最大，其貢獻(xiàn)率最高且與總氮負(fù)荷呈極顯著相關(guān)，而天然林具有涵養(yǎng)水源、保持水土等生態(tài)功能，因此起到“匯”的作用，滯緩非點(diǎn)源污染的發(fā)生。

表9 不同土地利用類型的總氮負(fù)荷

表10 土地利用類型的面積與總氮負(fù)荷的Pearson相關(guān)系數(shù)

3 討論與結(jié)論

3.1 討論

降雨是非點(diǎn)源污染產(chǎn)生的直接驅(qū)動(dòng)力[28]，研究區(qū)內(nèi)存在明顯的干濕兩季，濕季降雨量占全年的 84.1%，其總氮輸出量占全年的 85%，總氮負(fù)荷與月降雨量的變化趨勢基本一致，隨著降雨量的增加，流域內(nèi)總氮輸出升高。

在降雨條件下，非點(diǎn)源污染的負(fù)荷強(qiáng)度還受到土地利用方式的影響。不同土地利用類型的單位面積氮負(fù)荷量從高到低依次為：耕地>果園>橡膠林>漿紙林>天然林，其與植被覆蓋度和化肥施用量相關(guān)，因此耕地的單位面積氮負(fù)荷量最大，相反，天然林最小。同理，流域內(nèi)的氮負(fù)荷變化和空間分布與土地利用的分布相關(guān)，1988—2016年天然林面積減少，耕地、橡膠林、漿紙林和果園種植面積增加，在此期間總氮負(fù)荷呈增加的趨勢；相比于 1988年，在相同的劃分條件下 2002年總氮負(fù)荷最大梯度等級的子流域增加了3個(gè)，2016年增加了4個(gè)，總氮增加的子流域均是由于耕地或其他經(jīng)濟(jì)林的種植面積增加了。不同土地利用總氮輸出的貢獻(xiàn)率從高到低依次為：天然林>橡膠林>漿紙林>耕地>果園，貢獻(xiàn)率的大小與面積有關(guān)；Pearson相關(guān)性分析與貢獻(xiàn)率分析的結(jié)果基本一致，天然林在流域內(nèi)所占面積最大，氮負(fù)荷貢獻(xiàn)率最高，與總氮負(fù)荷的關(guān)系均表現(xiàn)為極顯著相關(guān)，但隨著面積的減小，其相關(guān)系數(shù)和貢獻(xiàn)率逐漸降低；耕地、橡膠林、漿紙林隨著面積的增加，貢獻(xiàn)率升高，與總氮負(fù)荷的相關(guān)系數(shù)逐漸增大，到了2016年達(dá)到顯著或極顯著相關(guān)。

當(dāng)前非點(diǎn)源污染已成為水體污染的一大主要污染來源，而由農(nóng)業(yè)造成的非點(diǎn)源貢獻(xiàn)率已超過50%。松濤水庫是海南省重點(diǎn)飲用水源保護(hù)區(qū)，其上游流域山多地廣，森林資源豐富，以農(nóng)業(yè)為主導(dǎo)產(chǎn)業(yè)，因此如何有效地防控非點(diǎn)源的發(fā)生已成為不容小覷的問題。本研究成果可為當(dāng)?shù)胤乐无r(nóng)業(yè)非點(diǎn)源污染提供參考。

3.2 結(jié)論

研究使用SWAT模型對松濤水庫上游流域進(jìn)行了非點(diǎn)源污染負(fù)荷的模擬，分析了土地利用變化對氮負(fù)荷的影響，結(jié)果如下：

（1）對以1988、2002和2016年土地利用數(shù)據(jù)為基礎(chǔ)構(gòu)建的SWAT模型進(jìn)行率定與驗(yàn)證，對月徑流率定和驗(yàn)證的納什效率系數(shù)Ens均為0.90以上，氨氮和總氮的Ens分別為0.54和0.46以上，說明SWAT模型適用于該流域的非點(diǎn)源污染模擬。

（2）降雨是造成非點(diǎn)源污染的外部動(dòng)力條件，濕季總氮負(fù)荷全年輸出約占85%，干季總氮負(fù)荷全年約占15%。1988—2016年流域內(nèi)氮負(fù)荷量呈增加的趨勢，氮負(fù)荷的空間分布與土地利用類型相關(guān)，總氮輸出量高的子流域內(nèi)均以耕地和經(jīng)濟(jì)林為主。

（3）單位面積氮負(fù)荷量從高到低依次為：耕地>果園>橡膠林>漿紙林>天然林；土地利用類型對全年總氮輸出的貢獻(xiàn)率由大到小依次為：天然林>橡膠林>漿紙林>耕地>果園；天然林與氮負(fù)荷的關(guān)系為極顯著相關(guān)，隨著面積的減小，其Pearson相關(guān)系數(shù)逐漸減小；隨著面積的增加，耕地、橡膠林、漿紙林與總氮負(fù)荷的Pearson相關(guān)系數(shù)逐漸增大，在2016年達(dá)到顯著或極顯著相關(guān)。