石金昊，朱衛(wèi)紅，田 樂，金日

（延邊大學(xué) 地理與海洋科學(xué)學(xué)院/長白山濕地生態(tài)系統(tǒng)功能與生態(tài)安全重點(diǎn)實(shí)驗(yàn)室，吉林 延吉133002）

0 引 言

【研究意義】據(jù)2010年《第一次全國污染源普査公報》數(shù)據(jù)顯示，面源污染已經(jīng)成為我國水體的主要污染源，是我國湖泊水質(zhì)惡化的主要原因。面源污染具有隨機(jī)性、復(fù)雜的遷移轉(zhuǎn)化機(jī)理、時空變化幅度大等特點(diǎn)[1]。面源污染是流域不同土地利用類型對降雨徑流的綜合響應(yīng)，其產(chǎn)生、輸出、遷移和轉(zhuǎn)化與土地利用類型及其不同組合與空間格局特征密切相關(guān)[2]，流域景觀的組成格局和空間格局變化也會影響面源污染物的產(chǎn)生、遷移和轉(zhuǎn)化。降雨對地表的沖刷通過陸地匯流將污染物帶入水體，土地利用和土地覆蓋方式影響了地表各類污染物的分布和負(fù)荷，用地類型的改變也影響了天然的水文過程。因此，運(yùn)用景觀生態(tài)學(xué)原理探究景觀格局變化與面源污染之間的關(guān)系成為面源污染控制的關(guān)鍵。

【研究進(jìn)展】隨著水文模型發(fā)展和應(yīng)用的不斷成熟，其中SWAT 模型被越來越多地用來評估流域尺度下的面源污染，分析其時空分布進(jìn)而識別關(guān)鍵污染區(qū)域和關(guān)鍵污染期，也用于分析和評價污染控制管理措施對水環(huán)境的影響[3]。Uriarte 等[4]研究了1977—2011年流域土地利用變化對水質(zhì)的影響，發(fā)現(xiàn)河流水質(zhì)惡化是因?yàn)槌鞘谢湍翀龅陌l(fā)展。Kibena 等[5]發(fā)現(xiàn)城市化進(jìn)程的加快和農(nóng)業(yè)用地的增加，使該地區(qū)的面源污染加劇。而李懷恩等[6]、宋林旭等[7]運(yùn)用SWAT 模型模擬了多年徑流、泥沙、氮磷等污染物濃度，表明不同的土地利用條件下產(chǎn)出的氮磷污染物濃度有較大差異，且植被覆蓋的增加可以很好地降低面源污染負(fù)荷。Wang 等[8]、孫麗娜等[9]利用SWAT 模型和CLUE-S 模型相結(jié)合的方法，分析了長時期土地利用變化下面源污染的輸出特征，并結(jié)合2 種模型模擬未來不同土地利用情景下的面源污染特征。【切入點(diǎn)】上述研究通常集中在流域內(nèi)的土地利用類型的構(gòu)成如何解釋面源污染的變化，很少有學(xué)者研究土地利用模式（景觀格局）的空間配置，例如斑塊的大小、形狀，邊緣配置或景觀內(nèi)土地利用類型的斑塊的空間連通性等。但土地利用類型指標(biāo)只能粗略預(yù)測水環(huán)境的變化，大多數(shù)研究沒有區(qū)分不同的景觀格局與面源污染的關(guān)系，所以難以全面反映土地利用變化對面源污染的影響?！緮M解決的關(guān)鍵問題】因此，基于SWAT模型分析布爾哈通河流域1986、1996、2006、2016年的面源污染時空分布特征，并結(jié)合實(shí)測的水文水質(zhì)數(shù)據(jù)進(jìn)行率定及驗(yàn)證?；谒鶚?gòu)建的布爾哈通河流域SWAT 模型，將重分類后的土地利用數(shù)據(jù)代入SWAT模型，分別模擬1986、1996、2006 和2016年4 期土地利用現(xiàn)狀下的面源污染時空分布，計(jì)算1986年和2016年布爾哈通河各子流域的景觀格局指數(shù)，分析子流域面源污染和景觀格局之間的關(guān)系，以期為該流域有效減少非點(diǎn)源污染提供科學(xué)依據(jù)。

1 研究區(qū)概況

布爾哈通河位于42°27′N—43°23′N，129°46′E—129°38′E 之間，是國際性河流圖們江的重要支流之一（圖1），源于安圖縣哈爾巴嶺山脈東南麓，橫穿延邊朝鮮族自治州中部平原盆地，最終匯入嘎呀河。河流全長約172 km，流域面積約7 064 km2，多年平均徑流量約63 198 萬m3[10]。全流域內(nèi)有耕地11.3 萬hm2，耕地率16%，其中水田2.85 萬hm2，占耕地面積的25.2%，且多集中在布爾哈通河的河谷盆地，是吉林省著名的水稻產(chǎn)區(qū)。由于開發(fā)歷史較早，墾殖率高達(dá)20%～30%，河谷內(nèi)的土地大多被開發(fā)成水田和建筑用地，河谷二側(cè)坡地開墾的旱田、果園非常普遍，原有的森林資源面積急劇減少，破碎化程度不斷加劇，水土流失面積占流域總面積的39%～49%[11]，加重流域面源污染風(fēng)險。

圖1 布爾哈通河流域DEMFig.1 DEM of Burhatong River Basin

圖2 布爾哈通河流域1986—2016年4 期土地利用Fig.2 The land use map for the fourth period of 1986—2016 in the Burhatong River Basin

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

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

①土地利用數(shù)據(jù)：本研究基于1986、1996、2006年和2016年布爾哈通河流域的Landsat系列衛(wèi)星遙感影像，先對遙感影像數(shù)據(jù)做處理工作，主要包括幾何校正、波段融合、波段合成、影像拼接與裁剪等。在ArcGIS10.1 和Cognition902 軟件的利用下，運(yùn)用面向?qū)ο蠓诸惙椒ńY(jié)合人工目視解譯方法進(jìn)行土地利用信息提取工作。從TM 圖像中把土地利用類型分為7 類：農(nóng)田、林地、草地、建設(shè)用地（住宅、商業(yè)和工業(yè)用地）、濕地、水域、裸地（圖2）；②為了劃分研究所需要的流域，利用NASA 提供的ASTER DEM（分辨率30 m）；③土壤數(shù)據(jù)來源于寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心，基于世界土壤數(shù)據(jù)庫中國土壤數(shù)據(jù)集（V1.1）；④氣象數(shù)據(jù)來自中科院水利水電科學(xué)研究院的SWAT 模型中國大氣同化驅(qū)動數(shù)據(jù)集[12]（CMADS，http://www.cmads.org/），時間尺度為2008—2017年，數(shù)據(jù)包含了SWAT 模型所需的降水、氣溫、風(fēng)速、太陽輻射、相對濕度；⑤水文數(shù)據(jù)來自榆樹川站點(diǎn)月徑流數(shù)據(jù)，時間尺度2010—2013年。月水質(zhì)監(jiān)測數(shù)據(jù)（總磷），時間尺度2011—2013年，該研究僅模擬流域總磷的時空分布。

2.2 分析方法

2.2.1 流域面源污染模擬

利用上述數(shù)據(jù)構(gòu)建SWAT 模型后，采用SWAT-CUP 進(jìn)行校準(zhǔn)和不確定性分析，以月徑流量、月總磷量為目標(biāo)對模型進(jìn)行校準(zhǔn)和驗(yàn)證，為判斷模型的模擬效果，分別選取決定系數(shù)R2和納什效率系數(shù)NSE評價模型精度[13]。校準(zhǔn)模型之后，將1986、1996、2006年的土地利用數(shù)據(jù)分別帶入模型，模擬4 期土地利用下的流域面源污染特征。

2.2.2 景觀格局分析

將1986年和2016年的土地利用數(shù)據(jù)由矢量數(shù)據(jù)轉(zhuǎn)化為柵格數(shù)據(jù)，然后按照SWAT 劃分出的子流域分別裁剪成44 份柵格數(shù)據(jù)，最后輸入Fragstats 4.2 軟件運(yùn)行，計(jì)算景觀格局指數(shù)。本研究在斑塊類型水平上選取PLAND（Percentage of Landscape）和PD（Patch Densit）來表征景觀格局的組成特征，這2 個指數(shù)是確定子流域優(yōu)勢斑塊和分析子流域不同類型斑塊破碎度的重要依據(jù)。在景觀水平上選取表征斑塊形狀復(fù)雜程度的LSI指數(shù)、表征斑塊聚集程度和連通性的CONTAG指數(shù)及COHESION指數(shù)、表征子流域景觀異質(zhì)性的SHDI指數(shù)，這幾個簡單且通用的景觀度量，能比較全面地反映子流域景觀空間特征。

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

冗余分析（RDA）現(xiàn)在被廣泛用于確定環(huán)境因素和生態(tài)指標(biāo)之間的關(guān)系，以44 個子流域?yàn)闃颖?，使用CANOCO5.0軟件分析面源污染負(fù)荷和景觀格局的關(guān)系，其中箭頭的長度代表該變量被排序圖解釋的程度，箭頭越長影響程度越高，當(dāng)各景觀指標(biāo)箭頭與水質(zhì)指標(biāo)箭頭之間夾角小于90°，二者關(guān)系為正相關(guān)，大于90°為負(fù)相關(guān)，等于90°則不存在相關(guān)性[14]。

3 結(jié)果與分析

3.1 流域景觀格局變化

林地作為布爾哈通河流域的主要土地利用類型，2016年的林地面積占流域面積的73%以上，并且從1986—2016年林地的面積有著小幅度的增加。由于開墾時間較早，農(nóng)業(yè)用地在流域內(nèi)的面積占比僅次于林地，在1986年達(dá)到21.45%，但到2016年縮減到了20.57%。受城市化進(jìn)程和經(jīng)濟(jì)發(fā)展的需要，建設(shè)用地從1986—2016年呈不斷擴(kuò)張的趨勢，其面積占比增長了0.92%，其余土地利用類型，僅有小幅變動。從1986—2016年，面積占比最大的林地有5 503.45 hm2被開墾成了農(nóng)田（表1），而大量的濕地被侵占成了農(nóng)田。因?yàn)?002 退耕還林工程在吉林省全面啟動，到了2016年，共有6 718.72 hm2的農(nóng)田轉(zhuǎn)換為林地，其次有4 577.50 hm2的農(nóng)田變?yōu)榱私ㄔO(shè)用地。1986—2016年，共有1 182.28 hm2的建設(shè)用地轉(zhuǎn)換成了農(nóng)田，農(nóng)田成為建設(shè)用地的主要轉(zhuǎn)變方向。

表1 1986—2016年土地利用轉(zhuǎn)移矩陣Table 1 1986—2016 Land use transfer matrix

本研究主要選取了林地、農(nóng)田、建設(shè)用地這3 種面積占比前3 的用地類型來分析其斑塊破碎化程度的改變，PD是表征流域斑塊破碎化程度的重要指標(biāo)，由表2 可知，1986—2016年，建設(shè)用地的PD平均值基本保持平衡，林地和農(nóng)田的PD平均值有著上升的趨勢，其中農(nóng)田的PD平均值上升幅度較大，破碎化程度較高。景觀水平上表征子流域斑塊聚集程度的CONTAG指數(shù)平均值呈下降趨勢，而斑塊連通性指數(shù)COHESION平均值基本維持不變。LSI和SHDI分別表征流域斑塊形狀復(fù)雜程度和景觀異質(zhì)性的指數(shù)，1986—2016年，二者的平均值都呈上升趨勢。

表2 1986—2016 布爾哈通河景觀格局變化總體趨勢Table 2 General trends of land use and landscape pattern changes in Burhatong River from 1986 to 2016

注 n 表示景觀中斑塊類型的數(shù)量。

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

3.2.1 模型校準(zhǔn)及驗(yàn)證

本研究通過SWAT-CUP 2012 進(jìn)行參數(shù)的全局敏感性分析，最后篩選出敏感程度較大的23 個參數(shù)，利用SWAT-CUP 程序里的SUFI-2 參數(shù)估計(jì)最優(yōu)化方法對榆樹川站點(diǎn)2010—2013年的月徑流量和月總磷量進(jìn)行校準(zhǔn)及驗(yàn)證。徑流在校準(zhǔn)期的R2和NSE 分別達(dá)0.86、0.81，在驗(yàn)證期的R2和NES 分別達(dá)0.72、0.72。總磷在校準(zhǔn)期的R2和NSE 分別達(dá)0.62、0.57，在驗(yàn)證期的R2和NSE 分別達(dá)0.71、0.69。一般情況下，R2＞0.6，NSE＞0.5 時，模型能有較好的模擬效果。綜合來看，本研究構(gòu)建的布爾哈通河流域SWAT模型達(dá)到了模擬精度的要求，可以很好地模擬流域水文水質(zhì)狀況。

圖3 1986—2016年4 期土地利用下子流域面源污染分布特征Fig.3 Distribution characteristics of non-point source pollution in each sub-basin under land use from 1986 to 2016

3.2.2 不同土地利用的面源污染時空變化特征

時間尺度上看，流域總磷流失負(fù)荷峰值主要發(fā)生在6—8月，降雨對總磷負(fù)荷的作用方式主要通過對徑流和輸沙量的影響產(chǎn)生，一般情況下降雨量和總磷流失量呈正相關(guān)?？偭棕?fù)荷與徑流量大致一致，表現(xiàn)為總磷負(fù)荷隨徑流量增加而增加。已有研究[15]表明，在夏季降雨量多，總磷負(fù)荷量也隨之增多，汛期面源污染負(fù)荷量通常較高。空間尺度上看（圖3），總磷的流失主要發(fā)生在流域北部和東部地區(qū)，集中在9、10、32、33、34 號子流域，在2016年，這5 個子流域的總磷流失負(fù)荷占總負(fù)荷的30%。1986、1996、2006、2016年的總磷負(fù)荷分別為74.04、73.78、82.50、128.31 t。1986—1996年，總磷流失量保持在相對穩(wěn)定的狀態(tài)，但從2006年以后，總磷流失負(fù)荷開始呈急劇上升趨勢，這一差異的原因可能是進(jìn)入快速城市化發(fā)展階段后，城市的快速擴(kuò)張和人類生活生產(chǎn)污水的排放直接加速了總磷流失。分析結(jié)果表明，建設(shè)用地和農(nóng)田都與總磷流失負(fù)荷正相關(guān)（圖4（b））。

3.2.3 景觀格局與面源污染的關(guān)系

通過RDA 排序分析探討1986年和2016年景觀格局與非點(diǎn)源污染過程之間的關(guān)系。由圖4 所知，林地與TP負(fù)相關(guān)，而建設(shè)用地和農(nóng)田與TP正相關(guān)。1986年，景觀指標(biāo)對總磷污染的解釋率高達(dá)86.43%（表3），總磷流失負(fù)荷貢獻(xiàn)度較高的景觀指標(biāo)主要是農(nóng)田、PD（林地）及PD（農(nóng)田）（圖4（a））。而在2016年，景觀指標(biāo)對總磷污染的解釋率為77.60%（表3），總磷流失負(fù)荷貢獻(xiàn)度排前二的景觀指標(biāo)變成了建設(shè)用地和PD（林地）（圖4（b））。1986—2016年，農(nóng)田的PD值呈上升趨勢，并且農(nóng)田的PD值與總磷負(fù)相關(guān)，表明農(nóng)田斑塊的破碎化在一定程度上減緩了面源污染的流失。林地在布爾哈通流域的占比達(dá)到了73%，盡管林地對改善面源污染有著積極作用，但林地的PD值與總磷正相關(guān)，因?yàn)榱值匕邏K的破碎化會造成水土流失，從而加大面源污染流失風(fēng)險。CONTAG與TP負(fù)相關(guān)，這表明破碎化嚴(yán)重且景觀斑塊較為分散的區(qū)域面源污染負(fù)荷更高。COHESIO與TP負(fù)相關(guān)，表明水質(zhì)惡化更可能與流域內(nèi)分散分布的土地利用有關(guān)，所以提高流域內(nèi)土地利用斑塊之間的連通性，減少斑塊破碎化對養(yǎng)分流失改善有積極作用。SHDI 和LSI都與總磷正相關(guān)，表明流域內(nèi)景觀異質(zhì)性的增加和形狀的復(fù)雜化，能夠?qū)е挛廴矩?fù)荷輸出風(fēng)險水平的增加。

表3 冗余分析各排序軸方差解釋率Table 3 Total variance explained by the ordination axis

圖4 1986 和2016年污染物流失量與景觀指標(biāo)的冗余分析Fig.4 Redundant analysis of pollutant loss and landscape pattern index in 1986 and 2016

4 討論

通過構(gòu)建布爾哈通河流域SWAT 模型發(fā)現(xiàn)，流域降雨和徑流過程對面源污染流失有著重要影響，總磷流失量和徑流量正相關(guān)，流域豐水期多出現(xiàn)在6—9月，總磷產(chǎn)生量峰值也主要在6—9月，所以該時期是非點(diǎn)源污染的關(guān)鍵期。農(nóng)村和城市是非點(diǎn)源污染的主要來源，但在不同的歷史發(fā)展階段，二者對面源污染的影響也不同?；? 期的流域土地利用數(shù)據(jù)，通過構(gòu)建SWAT 模型分析因土地利用變化而導(dǎo)致的面源污染時空分布特征，發(fā)現(xiàn)受地形、農(nóng)業(yè)生產(chǎn)及城市發(fā)展影響的9、10、32、33 和34 號子流域總磷污染負(fù)荷較高，屬于流域污染的關(guān)鍵控制區(qū)。農(nóng)田和建設(shè)用地是面源污染的主要來源，而森林對面源污染改善有著積極影響，這與前人研究的結(jié)果基本一致[16]。具有污染輸出效應(yīng)的建設(shè)用地和農(nóng)田一直被視為面源污染的“源”景觀，其中建設(shè)用地對水質(zhì)的影響主要取決于流域內(nèi)城市化水平的高低。大量的生活垃圾、工業(yè)污水的排放及農(nóng)業(yè)生產(chǎn)使用的化肥、農(nóng)藥等會經(jīng)過城市不透水表面進(jìn)入河流。布爾哈通河中游到下游流域二岸分布城市、村屯、耕地等，來自城區(qū)的生活污水、農(nóng)業(yè)污水及工業(yè)廢水對流域面源污染有著負(fù)面影響。該流域傳統(tǒng)的農(nóng)業(yè)活動（即耕作或過度施肥）增加了土壤侵蝕的風(fēng)險，并且過量的養(yǎng)分應(yīng)用可能導(dǎo)致非點(diǎn)源污染經(jīng)地表徑流后流入到河流。此外，通過對流域的景觀格局的進(jìn)一步分析，發(fā)現(xiàn)了農(nóng)田對面源污染貢獻(xiàn)度降低的可能原因。有研究表明，農(nóng)田斑塊的破碎化程度越低，斑塊廊道聯(lián)通性越好，農(nóng)田污染物容易形成面源污染[17]。盡管33 號子流域的建設(shè)用地面積最大，但33 號子流域的污染輸出負(fù)荷比32 號和34 號小。已有研究表明，流域“源”和“匯”景觀的空間位置不同，其對水環(huán)境的影響也不同[18]。33號子流域的建設(shè)用地和農(nóng)田基本位于上游區(qū)域，下游區(qū)域則是大面積的森林，這對該流域的面源污染起到了更好的凈化作用。

森林一直被認(rèn)為與農(nóng)業(yè)面源污染有著直接的聯(lián)系，其植被覆蓋可以減少土壤侵蝕，具有削減暴雨徑流、減少水土流失、吸附污染物的作用，可以有效減少地表徑流攜帶營養(yǎng)鹽進(jìn)入河流，對水質(zhì)惡化在一定程度上有削減作用[19]，但往往局限于景觀組成格局的角度分析得出該結(jié)論。以具有污染物消減效應(yīng)的森林為優(yōu)勢斑塊的流域景觀格局中，其對景觀破碎程度最為敏感[20]。林地景觀越破碎，斑塊越分散，越不利于抑制土壤侵蝕、徑流等過程，從而造成污染元素的流失。本研究也通過冗余分析發(fā)現(xiàn)林地景觀斑塊破碎化的增加，會增大流域面源污染產(chǎn)生的風(fēng)險，而林地作為布爾哈通河流域的主要用地類型，其景觀空間格局變化引起的生態(tài)變化，更值得注意和警惕。

區(qū)域水質(zhì)是多尺度環(huán)境因子的綜合反映，不單單受土地利用的影響，景觀格局的變化也被認(rèn)為是影響區(qū)域水質(zhì)的主要原因[21]，景觀格局與流域非點(diǎn)源污染研究是探究景觀格局與地表環(huán)境變化過程的典型模式。景觀形狀指數(shù)（LSI）、香濃多樣性指數(shù)（SHDI）與總磷污染呈正相關(guān)，主要因?yàn)檫@些指標(biāo)反映了子流域中斑塊的多樣性，這些多樣性的增大又間接地體現(xiàn)了景觀的異質(zhì)性與人類活動干擾的增強(qiáng)，因此也增加了面源污染的風(fēng)險，其他相關(guān)研究結(jié)果也表明這些指標(biāo)與水質(zhì)的惡化有密切聯(lián)系[22]。斑塊聚集度指數(shù)（CONTAG）反映的是景觀的分離與散布程度，其高值反映了高度聚集。例如，當(dāng)不同的土地利用最大限度地聚集時，CONTAG接近100，而當(dāng)土地利用最大限度地分散時接近與0。斑塊連通性指數(shù)（COHESIO）表示斑塊內(nèi)細(xì)胞的空間連通性和流域內(nèi)土地利用的物理連通性，高蔓延度值和高連通性說明景觀中的某種優(yōu)勢斑塊類型高達(dá)集聚并且形成了良好的連接性，景觀的破碎化程度較低，所以面源污染風(fēng)險也會隨之減小。有相關(guān)研究表明坡度主要通過影響產(chǎn)流量和產(chǎn)沙量進(jìn)而影響磷素的流失總量[23]，而在布爾哈通河流域坡度較高區(qū)域，在雨水充沛且多暴雨出現(xiàn)的條件下，會加重水土流失，致使大量未被植物吸收利用的氮磷營養(yǎng)元素隨水土流失進(jìn)入河流，是造成流域面源污染的重要原因。

5 結(jié)論

1）SWAT 模型能較好地模擬布爾哈通河流域總磷污染時空分布，反映研究區(qū)的總磷污染流失負(fù)荷。

2）布爾哈通河流域面源污染防治的關(guān)鍵控制區(qū)域在坡度較高的北部區(qū)域和農(nóng)田和城鎮(zhèn)集中的東部區(qū)域。

3）通過冗余分析發(fā)現(xiàn)林地與總磷負(fù)荷負(fù)相關(guān)，農(nóng)田和建設(shè)用地與總磷負(fù)荷正相關(guān)；林地斑塊破碎化越高，流域面源污染流失越嚴(yán)重，農(nóng)田斑塊的破碎化在一定程度上減緩了面源污染的擴(kuò)散。流域斑塊之間的高度流通性和聚集性，對養(yǎng)分流失改善有積極作用。流域內(nèi)景觀異質(zhì)性的增加和形狀的復(fù)雜化，能夠?qū)е挛廴矩?fù)荷輸出風(fēng)險的增加。對該區(qū)域進(jìn)行土地利用規(guī)劃與景觀格局優(yōu)化時，應(yīng)從景觀組成格局和空間格局二方面綜合考慮，才能更好地改善流域面源污染。