鄭 宇，程香菊，王兆禮,2，賴成光,2

(1.華南理工大學(xué)土木與交通學(xué)院，廣東 廣州 510641； 2.華南理工大學(xué)亞熱帶建筑科學(xué)國家重點實驗室，廣東 廣州 510641)

面源污染導(dǎo)致河流湖泊等水體出現(xiàn)富營養(yǎng)化和生物多樣性損失，并威脅人們的用水安全[1-2]。我國面源污染尤其是農(nóng)業(yè)面源污染已經(jīng)成為一個重要的水環(huán)境問題[3-4]。面源污染的產(chǎn)生和發(fā)展與地形、氣象、土壤、土地等諸多因素有關(guān)[5-6]，而流域景觀的組成格局和空間格局變化會改變這些因素，進而影響面源污染物的產(chǎn)生、遷移和轉(zhuǎn)化[7-9]。因此，以景觀生態(tài)學(xué)原理為指導(dǎo)，研究流域面源污染與景觀變化之間的關(guān)系，有助于對景觀的組成格局和空間格局進行有效的規(guī)劃和管理，從而改善流域的面源污染。近年來，流域水文水質(zhì)模型結(jié)合數(shù)理統(tǒng)計分析的方法被越來越多地應(yīng)用于研究面源污染與景觀格局的關(guān)系[10]，但大部分研究重點關(guān)注景觀的組成格局(例如土地利用面積比例)與面源污染的關(guān)系，而忽視了景觀的空間格局(例如景觀的斑塊密度、形狀指數(shù)等)與面源污染的關(guān)系，因而難以全面反映景觀對面源污染的影響[11]。

韓江流域是廣東省內(nèi)的第二大流域，從20世紀80年代開始，隨著經(jīng)濟的發(fā)展，流域內(nèi)的景觀發(fā)生了較大變化，對生態(tài)環(huán)境造成了深刻影響。韓江的上游梅江水質(zhì)日趨惡化，在某些河段存在富營養(yǎng)化現(xiàn)象，部分省界斷面水質(zhì)已經(jīng)嚴重惡化，威脅下游潮州、汕頭等市人民的用水安全。本文基于水土評價工具(soil and water assessment tool, SWAT)模型構(gòu)建面源污染模型，定量模擬和分析韓江流域面源污染的時空分布特征，識別關(guān)鍵污染地區(qū)，并應(yīng)用多元線性回歸、冗余分析等方法，探討流域面源污染與景觀格局之間的關(guān)系，以期為流域水環(huán)境管理提供依據(jù)。

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

1.1 研究區(qū)概況

以韓江流域潮安站以上部分為研究區(qū)，在SWAT模型中處理劃分為48個子流域(圖1)，總面積約29 077 km2。韓江流域擁有800多萬人口，流域內(nèi)的多年平均降雨量為1 600 mm，主要土地利用類型是林地。流域內(nèi)地形變化較大，高程在20～1 500 m之間，坡度在11.7°～20.0°之間。近年來，由于城市化的快速發(fā)展和水利工程的建設(shè)，流域的水質(zhì)逐漸惡化[12]。

1.2 數(shù)據(jù)來源

本文所用的主要數(shù)據(jù)見表1，其中，依照全國土地利用分類標(biāo)準(zhǔn)，土地利用類型數(shù)據(jù)被重新劃分為水田、旱地、草地、林地和城鎮(zhèn)5個主要類型；氣象數(shù)據(jù)包括日平均降雨、最高和最低氣溫、相對濕度、風(fēng)速及日照時長等；水文水質(zhì)數(shù)據(jù)包括徑流量、含沙量、氨氮質(zhì)量濃度等的監(jiān)測數(shù)據(jù)。除了表1所列數(shù)據(jù)，還根據(jù)研究區(qū)內(nèi)12個縣的統(tǒng)計年鑒等資料統(tǒng)計了各縣的人口、禽畜與水產(chǎn)養(yǎng)殖數(shù)量、耕地面積和化肥用量等數(shù)據(jù)，按照第一次全國污染源普查得到的排污系數(shù)成果，估算出研究區(qū)的生活污水、水產(chǎn)養(yǎng)殖、禽畜養(yǎng)殖和無機化肥施用產(chǎn)生的污染負荷并添加到模型中，以提高模擬的精度。模型中的農(nóng)業(yè)管理措施通過走訪和文獻調(diào)研獲得，例如研究區(qū)的主要作物被設(shè)定為雙季水稻和花生，早稻于每年的4月種植，7月收割；晚稻于每年的8月種植，11月收割；花生于每年的3月種植，7月收割。

圖1 韓江流域的地理位置、水系、站點與子流域分布

1.3 流域面源污染模擬

SWAT是美國農(nóng)業(yè)部開發(fā)的流域尺度的時間上連續(xù)、空間上呈半分布式的模型，可用于模擬土地利用、土地管理措施、氣候變化等因素對流域水文水質(zhì)的影響。本文利用地形、氣象、土地利用、土壤、管理措施等數(shù)據(jù)構(gòu)建SWAT模型，用SWAT-CUP(SWAT calibration/uncertainty program)校準(zhǔn)和不確定性分析程序，以徑流量、含沙量和氨氮質(zhì)量濃度為目標(biāo)對所建模型進行校準(zhǔn)和驗證，并選用決定系數(shù)R2、納什系數(shù)ENS和相對誤差ER等指標(biāo)來評價模型的擬合效果[13]，選用fp因子和fr因子來評價模型的不確定性[14]。用校準(zhǔn)后的模型分別加載3期的土地利用數(shù)據(jù)，輸出各子流域出口處1990—2014年的總氮、總磷、有機氮、有機磷、氨氮、無機磷等污染物的多年平均質(zhì)量濃度，以進行進一步的分析。

表1 研究采用的主要數(shù)據(jù)

1.4 景觀格局分析

景觀格局可以分為組成格局和空間格局，并用相應(yīng)的景觀指標(biāo)來衡量[15-16]。本文子流域景觀格局的組成特征用某一土地利用類型面積百分比PLS和最大斑塊指數(shù)PLP來衡量，這兩個指標(biāo)是確定景觀中優(yōu)勢元素的重要依據(jù)，也是決定景觀中生物多樣性和優(yōu)勢物種的重要因素。景觀格局的空間特征用某一土地利用類型斑塊密度DP和形狀指數(shù)ILS來衡量，其中斑塊密度反映了景觀的破碎程度，形狀指數(shù)反映了景觀的形狀復(fù)雜度。這些指標(biāo)的計算公式如下：

(1)

(2)

(3)

(4)

式中：n為子流域中某種土地利用類型的斑塊數(shù)量；ai為斑塊i面積，m2；A為子流域面積，m2；ei為斑塊i的邊緣長度，m。

1.5 統(tǒng)計分析

在RStudio中，利用多元線性回歸和冗余分析研究景觀格局與面源污染的關(guān)系。多元線性回歸能夠反映單個因變量與多個自變量之間的關(guān)系，而冗余分析是多元線性回歸的擴展，能夠同時反映所有因變量與自變量間的關(guān)系。由于自變量數(shù)目較多，本文在多元線性回歸中使用全子集回歸對自變量集進行縮減，并根據(jù)回歸診斷的結(jié)果選擇最優(yōu)的回歸方程；在冗余分析中使用向前選擇對自變量集進行縮減。

2 結(jié)果與討論

2.1 流域的景觀格局特征

在1990年和2010年土地利用背景下，整個流域內(nèi)各土地利用類型的面積占流域總面積的百分比如表2所示。和其他土地利用類型相比，盡管城鎮(zhèn)用地的面積比例是最小的，但是增長率卻是正值且是最大的，說明研究區(qū)正在經(jīng)歷城市化的初始階段，各種土地有轉(zhuǎn)換為城鎮(zhèn)用地的趨勢。在2010年土地利用背景下，部分典型子流域中各土地利用類型的面積百分比PLS、最大斑塊指數(shù)PLP、斑塊密度DP、形狀指數(shù)ILS如表3所示。可以發(fā)現(xiàn)，林地在每個子流域中都擁有最高的面積百分比和最大斑塊指數(shù)，說明研究區(qū)內(nèi)的土地利用類型以林地為主。人口密集的子流域，例如興寧盆地附近的33、36、37、46號子流域，農(nóng)業(yè)和城鎮(zhèn)用地的面積百分比更高。水田的面積百分比一般大于旱地，說明流域內(nèi)的主要農(nóng)業(yè)活動是水稻種植。上游子流域的草地面積百分比普遍大于下游子流域。盡管草地的面積百分比一般大于水田，水田的最大斑塊指數(shù)卻常常大于草地，說明流域內(nèi)草地的分布相比水田更為分散。城鎮(zhèn)的斑塊密度和形狀指數(shù)相比其他土地利用類型更低，表明在人類活動的影響下，城鎮(zhèn)斑塊的破碎程度更小，形狀更規(guī)則。

表2 各土地利用類型面積占全流域面積的百分比 %

2.2 流域的面源污染時空分布特征

在SWAT-CUP軟件中通過全局敏感性分析篩選出較為敏感的重要參數(shù)，例如地表徑流滯時、基流退水常數(shù)、曼寧系數(shù)、徑流曲線數(shù)、不同土地利用類型的流失系數(shù)等，再用內(nèi)置的SUFI-2(sequential uncertainty fitting version 2)算法進行校準(zhǔn)和驗證，結(jié)果如表4所示。在校準(zhǔn)期和驗證期，徑流量、含沙量、氨氮質(zhì)量濃度的決定系數(shù)R2和納什效率系數(shù)ENS基本都超過0.7，相對誤差ER基本都小于8.5%，說明模型的擬合效果良好。所有站點的fp因子都大于0.5，同時fr因子都比較小，表明模型具有較小的不確定性。因此，所建立的模型基本上能反映流域的面源污染規(guī)律。

應(yīng)用SWAT模型，分別加載1990年和2010年的土地利用數(shù)據(jù)進行模擬，并輸出各子流域出口處的總氮、總磷質(zhì)量濃度以及它們的增量，結(jié)果如圖2、圖3所示。從空間上看，各子流域的污染負荷相差較大，總氮、總磷質(zhì)量濃度呈現(xiàn)明顯的空間差異性。位于西部寧江流域地區(qū)(28號及其上游的子流域)的總氮質(zhì)量濃度較高，并在興寧市及附近地區(qū)(36號子流域)達到最大值3.91 mg/L，這主要是因為該地區(qū)位于農(nóng)業(yè)較為發(fā)達的寧江平原，農(nóng)田徑流、禽畜糞便、生活污水的產(chǎn)生量較大，產(chǎn)生了較大的總氮流失。流域東部的汀江流域地區(qū)的總磷質(zhì)量濃度普遍較高，并在龍巖市永定區(qū)及附近地區(qū)(21號子流域附近)達到最大值2.15 mg/L，這可能是因為這些地區(qū)的草地和林地面積比重普遍較大，而草地相比其他土地利用類型具有更高的磷流失風(fēng)險[17]。從1990—2010年，大部分子流域的總氮負荷呈現(xiàn)增加趨勢。雖然在大部分子流域，農(nóng)田、城鎮(zhèn)等人類用地的面積比自然土地面積小，但是在整個子流域的出口處(45號子流域)，各污染物的質(zhì)量濃度仍然較高且呈現(xiàn)上升趨勢，這反映出流域內(nèi)人類活動對水環(huán)境的影響大于自然過程。

表3 部分子流域2010年各土地利用類型的景觀格局指標(biāo)

表4 SWAT模型的擬合效果評價

(a) 1990年

(b) 2010年

(c) 增加量

圖2 各子流域的總氮質(zhì)量濃度模擬值

(a) 1990年

(b) 2010年

(c) 增加量

圖3 各子流域的總磷質(zhì)量濃度模擬值

2.3 流域的景觀格局變化與面源污染變化的關(guān)系

以1990—2010年土地利用背景下各污染物質(zhì)量濃度的增量為因變量，各景觀格局指標(biāo)的增量為自變量，建立多元線性回歸方程如下：

yTN=0.536xDP5+0.093xPLP2

(5)

yON=-0.132xDP1+0.056xPLP2-0.003xILS5

(6)

yAN=-0.016xPLP5-0.015xILS3+
0.014xILS1-0.002xPLS3

(7)

yTP=-0.233xDP1+0.079xPLP2-0.005xILS5

(8)

yOP=-0.184xDP1+0.071xPLP2-0.004xILS5

(9)

yMP=-0.073xDP1+ 0.009xPLP2-0.002xPLS3

(10)

式中：yTN、yON、yAN、yTP、yOP、yMP分別為總氮、有機氮、氨氮、總磷、有機磷、無機磷質(zhì)量濃度的增量，mg/L；xPLS、xPLP、xDP、xILS分別為面積百分比、最大斑塊指數(shù)、斑塊密度、形狀指數(shù)的增量；數(shù)字下標(biāo)1、2、3、4、5分別表示水田、旱地、林地、草地和城鎮(zhèn)。

經(jīng)過全子集回歸的篩選，進入不同回歸方程的自變量不同，說明景觀格局對不同污染物的貢獻程度各異。旱地最大斑塊指數(shù)的增量與大部分污染物質(zhì)量濃度的增量呈現(xiàn)正相關(guān)關(guān)系，水田斑塊密度的增量與大部分污染物質(zhì)量濃度的增量呈現(xiàn)負相關(guān)關(guān)系，林地形狀指數(shù)和面積百分比的增量與氨氮質(zhì)量濃度的增量呈現(xiàn)輕微的負相關(guān)關(guān)系。城鎮(zhèn)斑塊密度的增量與總氮質(zhì)量濃度的增量呈現(xiàn)強正相關(guān)關(guān)系，表明城鎮(zhèn)斑塊可能是氮污染的主要貢獻者；相反城鎮(zhèn)的形狀指數(shù)增量與有機氮、總磷、有機磷質(zhì)量濃度增量呈現(xiàn)輕微的負相關(guān)關(guān)系；城鎮(zhèn)的面積百分比和最大斑塊指數(shù)沒有進入多數(shù)回歸方程，表明城鎮(zhèn)的組成格局對面源污染的影響沒有空間格局的影響大。

在多元線性回歸的基礎(chǔ)上，通過冗余分析來同時反映所有因變量與自變量間的關(guān)系。在冗余分析中，因變量的總方差被分解成為約束性和非約束性方差。由表5可見，前兩個軸變量已經(jīng)能解釋81.7%的總方差，以這兩個變量為橫軸和縱軸，將各自變量、因變量和子流域繪制成圖4所示的三序圖。從圖4可以看出，大部分子流域沿R2軸方向分布，而R2軸與城鎮(zhèn)的xPLP、xDP和xILS形成的夾角較小，說明流域內(nèi)各子流域之間的城鎮(zhèn)化程度差異較大。進入冗余分析模型中的景觀指標(biāo)主要是城鎮(zhèn)和旱地的景觀指標(biāo)，這說明污染物質(zhì)量濃度的變化對城鎮(zhèn)和農(nóng)田景觀指標(biāo)的變化更敏感。城鎮(zhèn)的最大斑塊指數(shù)增量、斑塊密度增量和形狀指數(shù)增量以及旱地的形狀指數(shù)增量、林地的面積百分比增量與各污染物質(zhì)量濃度的增量呈現(xiàn)負相關(guān)關(guān)系，而旱地和草地的面積百分比增量與各污染物質(zhì)量濃度的增量呈現(xiàn)正相關(guān)關(guān)系。

表5 冗余分析中各軸變量解釋的方差比例 %

2.4 討 論

相關(guān)研究表明，土地利用類型可以被劃分成“源”景觀和“匯”景觀，前者對流域內(nèi)的泥沙流失、污染物流失具有較大貢獻[18]，而后者對泥沙流失和營養(yǎng)物流失的貢獻較小甚至有削弱作用[19]。林地通常被認為是一種“匯”景觀格局，因為它具有保持水土、凈化水質(zhì)的功能[20-21]。圖4亦表明，林地的面積百分比與大部分污染物的質(zhì)量濃度負相關(guān)。

由于削減徑流、反硝化等作用，草地也常被認為是一種“匯”景觀格局[22-23]，但也有研究指出草地增加了水環(huán)境污染的風(fēng)險[24-25]。圖4表明，草地面積百分比與污染物質(zhì)量濃度呈正相關(guān)關(guān)系，說明研究區(qū)內(nèi)的草地更像一種“源”景觀格局。這主要是由以下原因?qū)е碌模孩傺芯繀^(qū)的坡度較大，平均坡度達到16.6°，增加了營養(yǎng)物隨土壤流失的風(fēng)險[26-28]；②研究區(qū)受亞熱帶氣候影響，降雨量較為充沛，多年平均降雨量超過1 400 mm，且80%集中在雨季，增加了土壤顆粒隨著地表徑流流失的風(fēng)險[29]；③從1990—2010年，研究區(qū)有197.3 km2的林地被轉(zhuǎn)化成草地，占了草地來源的78.5%，而草地對污染物的削減作用比林地弱，從而導(dǎo)致了徑流中營養(yǎng)物通量的增加[30]。因此，盡管草地景觀具有調(diào)節(jié)徑流和削減污染的作用，依然不足以抗衡大坡度和強降雨的影響。本文得到的結(jié)果表明流域管理中考慮坡度的重要性，例如在坡度較大的地區(qū)退耕退草還林可能是降低流域面源污染的有效措施。為論證這一結(jié)論，設(shè)置了將10°以上地區(qū)的草地轉(zhuǎn)換為林地的管理情景并進行模擬，結(jié)果如表6所示，該情景下流域出口處的多年平均總氮通量削減了1.82%，總磷通量削減了1.16%。

由于化肥的施用以及土壤的松動，農(nóng)業(yè)用地是造成水環(huán)境惡化的重要因素之一[31]。式(5)～(10)和圖4表明農(nóng)業(yè)用地最大斑塊指數(shù)的增量與污染物質(zhì)量濃度的增量呈現(xiàn)正相關(guān)關(guān)系，但是，與水田相比，進入回歸模型和冗余分析模型中的旱地指標(biāo)更多，說明旱地與污染物質(zhì)量濃度的相關(guān)性比水田與污染物質(zhì)量濃度的相關(guān)性更強。這可能是由以下原因造成的：①農(nóng)民傾向于在旱地施用更多肥料，導(dǎo)致農(nóng)田徑流中的氮磷等營養(yǎng)物質(zhì)量濃度更高[32]；②旱地上的作物可能會產(chǎn)生更大的滲漏淋失[33]，導(dǎo)致旱地土壤中的有機質(zhì)和總氮質(zhì)量濃度比水田中的更高[34]；③水田具有更長的水力保留時間和缺氧環(huán)境，促進了反硝化作用，從而減少了總氮的產(chǎn)生。因此，調(diào)整農(nóng)業(yè)結(jié)構(gòu)，例如將旱地改為水田可能有助于減少肥料用量和土壤流失。為了證明這一結(jié)論，設(shè)置了將興寧市的旱地改為水田的管理情景并進行模擬，結(jié)果如表6所示，該情景下流域出口處的多年平均總氮通量削減了4.63%，總磷通量削減了7.31%。需要注意的是，本文提出的管理措施主要是為了佐證通過數(shù)理統(tǒng)計方法得到的面源污染與景觀指標(biāo)之間的關(guān)系，因此在設(shè)定模擬情景時保持了模型數(shù)據(jù)庫中各種土地利用類型的參數(shù)(如污染物去除率、城鎮(zhèn)不透水面積等)不變，沒有考慮設(shè)置過濾帶、采取低影響開發(fā)措施、改善污水處理工藝、設(shè)置植草水道等具體工況的效果，導(dǎo)致模擬結(jié)果中對污染物的削減量有限；另外，雖然考慮到措施的可行性，只對局部地區(qū)施行了改造措施，但是沒有進行詳細的經(jīng)濟效益評估。如何結(jié)合本文關(guān)于面源污染與景觀格局關(guān)系的結(jié)論，提出更加立體、有效、經(jīng)濟的流域管理方案，是一個有待進一步研究的課題。

除了組成格局，景觀的空間格局也與水質(zhì)相關(guān)。例如，水田的斑塊密度與有機氮和各種形態(tài)磷的質(zhì)量濃度負相關(guān)，旱地的形狀指數(shù)與大部分污染物的質(zhì)量濃度負相關(guān)。本文和前人的研究結(jié)果都表明，城鎮(zhèn)景觀的斑塊密度和形狀指數(shù)與各種形態(tài)的氮磷污染物質(zhì)量濃度負相關(guān)[35]，這主要是因為城鎮(zhèn)景觀受人類活動影響較大，因而擁有更密集的空間分布和更規(guī)則的形狀[36]。

式(5)～(10)和圖4表明林地的空間格局指標(biāo)與污染物質(zhì)量濃度的相關(guān)性沒有組成格局指標(biāo)與污染物質(zhì)量濃度的相關(guān)性強，而城鎮(zhèn)景觀則相反?？紤]到林地和城鎮(zhèn)分別是研究區(qū)每個子流域中面積最大和最小的土地利用類型，推測景觀的組成格局和空間格局對水質(zhì)的貢獻程度與景觀的面積有關(guān)。面積大的景觀如林地和草地，其空間格局與水質(zhì)的相關(guān)性更強；面積小的景觀如農(nóng)田和城鎮(zhèn)，其組成格局與水質(zhì)的相關(guān)性更強。

表6 不同管理情景下流域出口處污染物的多年平均通量和質(zhì)量濃度

3 結(jié) 論

a. 研究區(qū)的土地利用以林地、草地為主，水田、旱地次之。流域處于城鎮(zhèn)化的初始階段，盡管城鎮(zhèn)面積比例較小，但增長率最高。下游子流域的城鎮(zhèn)化程度比上游高，呈現(xiàn)出向上游蔓延的趨勢。在以后的研究中模擬未來土地利用情況或者進行流域管理規(guī)劃時，應(yīng)當(dāng)重視城鎮(zhèn)的影響。

b. 在當(dāng)前土地利用策略和人類活動的影響下，研究區(qū)內(nèi)的水質(zhì)呈現(xiàn)惡化趨勢。研究區(qū)西部的寧江流域地區(qū)是總氮污染的關(guān)鍵區(qū)，總氮質(zhì)量濃度最大3.91 mg/L；東部龍巖市永定區(qū)的總磷質(zhì)量濃度較高，最大值為2.15 mg/L。

c. 流域的面源污染與景觀格局存在密切的關(guān)系。景觀指標(biāo)能解釋81.7%的面源污染變化。景觀組成格局指標(biāo)與空間格局指標(biāo)的作用與景觀面積有關(guān)，景觀面積越大，其組成格局可能比空間格局對水質(zhì)的影響更大。林地的面積百分比與面源污染呈現(xiàn)負相關(guān)關(guān)系，其他景觀則相反。城鎮(zhèn)或農(nóng)田等人類景觀的斑塊密度、形狀指數(shù)等空間格局指標(biāo)與污染物質(zhì)量濃度具有較強的負相關(guān)性。旱地相比水田對研究區(qū)面源污染的貢獻更大。

d. 在坡度較大的地區(qū)退草還林，可以削減1.82%的總氮和1.16%的總磷；將農(nóng)業(yè)結(jié)構(gòu)由旱地為主調(diào)整到以水田為主，可以削減4.63%的總氮和7.31%的總磷。