東江湖流域未來土地利用變化對面源污染負荷的影響評估

鄭 昊，賀華翔，李海明，牛存穩，王嘉浩，李昕陽

（1.天津科技大學 濱海地下水利用與保護研究室，天津 300457；2.中國水利水電科學研究院 流域水循環與調控國家重點實驗室，北京 100078；3.華北水利水電大學，河南 鄭州 450046）

0 引言

流域面源污染是生態環境保護的重點關注內容，隨著新時期國土空間規劃體系的建立，以生態環境保護和高質量發展為前提確立的土地利用格局，形成了新的流域面源污染時空分布特征。由于流域面源污染的分散性、隨機性、累積性等特點，其對土地利用格局變化的響應具有空間上的變異性和時間上的不均勻性［1］。

SWAT模型是流域面源污染模擬和污染治理效果評估的常用方法［2］。李亞嬌［3］等將4種典型面源污染模型SWAT、GLWF、SPARROW、HSPF 進行對比分析，其研究認為SWAT 模型模擬精度高，在國內具有普遍適用性。SWAT 模型被廣泛應用于土地利用變化對面源污染的評估，秦耀民［4］等對黑河流域不同情景的土地利用格局進行模擬，說明流域內實行退耕還林政策有利于面源污染的控制，劉瑤［5］等定量分析土地利用變化下昌江流域的面源污染負荷，確定水田在該地區對面源污染產生有較大的影響。但是土地利用數據作為SWAT面源污染模型的輸入邊界，難以實現土地利用變化與面源污染的時空特征關系的動態解析。因此，學者們嘗試將面源污染模型與土地利用模型相耦合，ZHANG 等［6］和LIU［7］等將SWAT 模型和CLUE-S 模型耦合，分別對密云水庫和渾太河流域不同土地利用情景下的面源污染進行評估；但是CLUE-S 模型對空間分配規則有主觀性且容易忽視土地利用數據空間的相關性，對結果準確度產生一定的影響［8］。石金昊［9］等利用SWAT 模型和FLUS 模型耦合預測2036年布爾哈通河流域不同土地利用情景下的面源污染特征，但是FLUS模型無法對缺少歷史資料地區進行模擬［10］。

PLUS 模型除了解決上述兩個模型存在的問題，還增加了動態模擬多類土地斑塊的生成和演化的功能，因此，本文將SWAT 模型和PLUS 模型相耦合，以東江湖流域為例，基于歷史土地利用格局，分析流域面源污染分布規律，預測未來不同土地利用變化情景下面源污染負荷變化量，并對其產生影響進行評估，為東江湖流域面源污染治理提供參考依據。

1 材料與方法

1.1 研究區概況

東江湖屬于湘江水系耒水上游，流域范圍涉及資興市、汝城縣、桂東縣和宜章縣4個縣市，由浙水、滁水、漚江三條主要河流匯入，流域面積4 639 km2。流域內地勢中心及西北偏低，東部較高，地表起伏較大。流域地形及行政區如圖1所示。研究區屬亞熱帶季風性濕潤氣候，雨量充沛，多年平均降雨量為1 487 mm［11］，降水集中在4-7月份，占全年降水量的60%。

圖1 流域地形及行政區劃圖Fig.1 Topography and administrative divisions of the basin

東江湖流域內主要土地利用類型為林地，面積為3 576.83 km2，占流域總面積的77.1%；其次為耕地，面積為636.76 km2，占流域總面積的13.2%；其余為草地、建設用地、水域等。耕地圍繞河流兩側地勢低緩的地區分布，以水田為主，旱地次之；建設用地跟耕地均分布在地勢平緩地區，建設用地從耕地心中向四周擴散，呈放射狀分布。人口分布和GDP 增長也集中在耕地和建設用地上，近年來，東江湖流域城市迅速擴張，人口與經濟穩步增長。土地利用、人口經濟分布如圖2所示。

圖2 土地利用、人口與經濟分布圖Fig.2 Distribution map of land use，population and economy

目前，東江湖流域面源污染問題突出，其中農村生活垃圾、農業生產及工礦企業產生的污染影響范圍廣，每年約9 000 余t COD、1 700 余t氨氮、300 余t總磷進入東江湖［12］，面源污染問題較為突出。國內學者對于東江湖流域面源污染研究較少，尤其是缺乏面源污染模型的構建與土地利用格局對面源污染的響應。本研究致力于對流域面源污染時空響應模擬，預測東江湖流域未來土地利用格局，解析土地利用格局對于面源污染的影響。

1.2 研究方法

研究耦合PLUS 模型和SWAT 模型用于評估未來土地利用變化對面源污染負荷的影響。利用PLUS 模型預測流域不同情景的土地利用，利用SWAT 模型模擬不同情景流域面源污染特征，將二者相耦合分析不同土地利用情景下面源污染的變化情況。

1.2.1 PLUS模型

PLUS 模型是一個以元胞自動機構建的新模型［13］，包括基于擴張分析策略（Land Expansion Analysis Strategy，LEAS）的規則挖掘框架以及多類隨機斑塊種子（CA-based on multiple random patch seeds，CARS）生成機制，可以有效地探索土地利用變化的原因，具有較高的模擬精度［14］。首先，把提取的土地利用擴張圖導入LEAS 模塊，疊加多個驅動因素，并設置隨機采樣點，采用隨機森林算法，獲取各類土地利用類型擴張規則，并生成相應的土地發展潛力圖。其次，CARS 模塊是由情景模擬驅動的土地利用預測模型，模擬過程中綜合了土地利用需求、鄰域權重和轉換成本矩陣3 個部分，使土地利用量達到未來用地需求［15］。東江湖流域土地利用模擬過程如下。

土地利用變化受自然因素和人文因素的影響［16］，驅動因子是PLUS 模型構建的重要數據，是推動土地利用類型發生變化的基礎，為構建土地利用模擬模型，本研究將影響土地利用變化的主要驅動因子分為自然氣候因子和社會經濟因子兩類，選取了八個驅動因子和一個限制因子（如表1所示），將所有驅動因子設置成與土地利用數據相同的行列數，并選取水域作為限制因子。

表1 土地利用變化驅動因素與限制因素Tab.1 Driving factors and limiting factors of land use change

本文采用PLUS 模型自帶的Markov Chain 模塊對研究區未來的各個類型的土地利用需求進行預測，通過模型計算2010年和2020年各類用地轉換概率，預測2035年東江湖流域的各類土地的需求。研究區內林地和耕地占較大比例，根據各個地類歷史發展趨勢并結合當地的生態保護政策和耕地保有紅線，盡量控制耕地和林地向建設用地轉換，因此對轉移矩陣進行設定如表2所示。其中0表示該地類不可轉化，1表示可以轉化。

表2 土地利用轉移矩陣Tab.2 Land use transition matrix

鄰域權重參數即該地類的擴張強度，范圍從0～1，越接近于1 表明該種土地利用類型擴張能力越強，鄰域權重值的確定一般采取兩種方式，無量綱量處理法和經驗法［17］，本研究對兩種方法進行對比后選取無量綱量處理法，最終鄰域權重參數值如表3所示。

表3 鄰域權重參數Tab.3 Neighborhood weight parameters

1.2.2 SWAT模型

SWAT 面源污染負荷模型主要分為水文過程子模型、土壤侵蝕子模型和污染負荷子模型［18］。

（1）水文子模型構建。將研究區內的16 個雨量站和2 個國家氣象站提供氣象數據，寨前（三）和東江2 個水文站提供水文數據（如圖1所示）輸入到模型中，根據DEM 的坡度、河長計算流積和流向，將最小子流域的面積設定為15 000 hm2，模型自動將研究區劃分為16個子流域（如圖2所示）。

（2）土壤侵蝕子模型和營養物質循環過程構建。將研究區內土壤類型概化為五種土壤類型，分別為堆積土、低活性強酸土、雛形土、高活性強酸土和高活性淋溶土；為了適應污染負荷模型的建立，將研究區內土地利用類型劃分為耕地、草地、水域、林地、建設用地和未利用地6 種常見一級土地利用類型，坡度劃分為0～15°、15～45°、45～90° 3類，再與坡度結合起來生成水文響應單元。流域內以雙季稻為主，種植時間為4月初和8月初，收割時間為7月下旬和11月上旬，耕作面積設定為639 km2，施用強度設置為氮肥53.82 kg/hm2，磷肥39.10 kg/hm2，復合肥為28.41 kg/hm2。

收集東江湖流域SWAT 模型構建的主要輸入數據如表4所示，包括高程數據、土壤數據、土地利用數據、氣象觀測數據。采用寨前（三）水文站的逐月徑流量和漚江水庫的月氨氮和總磷濃度（mg/L）的實測數據對模型進行率定驗證，本研究主要模擬指標確定為氨氮和總磷。

表4 SWAT模型構建數據表Tab.4 SWAT model construction data table

1.3 模型精度驗證

1.3.1 PLUS模型精度驗證

采用Kappa系數對PLUS 模型精度進行檢驗。Kappa系數計算公式如下：

式中：Pa為正確模擬的比例；Pc為隨機情況期望模擬比例；Pp為理想情況下正確模擬比例［19］。

經對比，本次構建的PLUS模型具有較好精度，以2020年土地利用數據為例，模擬與實際數據相比，Kappa系數可達0.90，高于Kappa精度達到0.8 以上的要求，說明PLUS 模式對研究區具有良好的適用性。

1.3.2 SWAT模型率定驗證結果

研究采用SWAT-CUP 軟件，選取SUFI-2 方法［18］進行月尺度流量及污染負荷模擬結果的率定驗證。

流量模擬方面，率定期為2009-2014年，驗證期為2015-2019年；污染負荷模擬方面，率定期為2019年1-6月，驗證期為2019年7-12月，均以納什效率系數NSE和線性回歸系數R2為評價指標。率定驗證結果情況如表5所示，擬合情況如圖3所示，結果表明模型模擬效果良好。

表5 水量水質率定驗證結果Tab.5 Verification results of water quantity and quality calibration

圖3 水量水質率定、驗證期擬合圖Fig.3 Water quantity and quality calibration and verification period fitting diagram

1.4 模型耦合及情景設計

將PLUS 模型和SWAT 模型進行精度驗證后，以PLUS 模型預測輸出的土地利用格局，嵌入SWAT模型空間數據庫中，替代原有土地利用模塊，與模型原有的土壤類型和坡度結合起來生成新的水文響應單元，重新運行模型，輸出未來不同情景下的面源污染負荷。耦合過程如圖4所示，在SWAT 其他參數數據保持不變情況下，預測不同情景土地利用格局下的面源污染空間分布特征。

圖4 模型耦合過程圖Fig.4 Model coupling process diagram

研究以2035年為典型年，設置兩種未來土地利用變化情景，具體情況如下：

（1）歷史趨勢情景：根據2010-2020年東江湖流域的土地利用類型變化趨勢，設置未來土地利用轉移矩陣，以此預測2035年土地利用布局，變化趨勢表現為耕地、草地和林地均向建設用地轉化。

（2）國土空間規劃情景：參考《郴州市國民經濟和社會發展第十四個五年規劃和2035年遠景目標綱要》的要求，郴州市土地利用發展趨勢遵循控制建設用地總量、確保林地數量有所增加和耕地保有底線，主要變化趨勢是在耕地保有紅線前提下，將耕地轉化為建設用地。

2 結果與分析

2.1 東江湖流域歷史土地利用演變分析

研究對比分析2010年、2020年兩個時期的土地利用數據，土地利用轉移矩陣如表6所示，2010-2020年土地利用轉換情況如圖5所示。經分析，與2010年相比，2020年建設用地增加最多，增加18.95 km2，其次是水域面積增加0.34 km2，其他的用地類型面積均減少，向建設用地轉出。體現了城市向外擴張侵占耕地、砍伐林地、吞并草地的過程，符合城市擴張規律。

表6 2010-2020年土地利用轉移矩陣km2Tab.6 Land use transition matrix from 2010 to 2020

圖5 2010年到2020年間土地利用轉換圖Fig.5 Land use conversion map from 2010 to 2020

2.2 東江湖流域面源污染時空演變特征

研究對比分析了2009年、2016年和2019年東江湖流域氨氮和總磷的面源污染特征，如圖6所示。

從時間尺度上看，氨氮和總磷輸出負荷集中在降水較大的4-7月份，流域面源污染受降水、徑流等自然過程影響，氨氮與總磷負荷與徑流量相關，表現為隨徑流量增加而增加，尤其是總磷，豐水期有明顯增加。從2009年到2016年，污染物負荷呈現逐漸增加趨勢，這除了與降水密切相關外還與城市迅速擴張和人口快速增長有關；2016年到2019年污染負荷產出相對減少，說明面源污染負荷受降水以及退耕還林等生態保護政策影響。

從空間尺度上看，將研究區按照最小子流域閾值為15 000 hm2劃分為16 個子流域，氨氮負荷較高子流域集中在河流西北部和中部徑流量較大的4 號、6 號、7 號和10 號子流域，其中4 號和6號子流域的氨氮負荷約占整個流域的50%。總磷負荷較高子流域集中在4 號、6 號、7 號、10 號和9 號子流域，約占整個流域總磷負荷的70%。主要原因如下：①因4 號子流域建設用地面積較大、人口較為密集，因此氨氮、總磷負荷較大；②4 號、6號、7號、10號子流域耕地沿河分布，耕地分布較為分散，易產生土壤侵蝕，污染物隨徑流進入東江湖；③子流域內因建設用地和水域擴張，造成對耕地、林地、草地的侵占，土壤侵蝕造成的面源污染負荷較高。

結合HRU輸出結果中各個地類總磷負荷，分析各用地類型對總磷污染負荷的貢獻，由表7 可以看出總磷在耕地上負荷最高，在未利用地上負荷最低，貢獻程度為耕地＞建設用地＞草地＞林地＞未利用地。由此可見，土地利用變化對污染負荷產生較大影響。

2.3 東江湖流域多情景土地利用預測模擬結果

分析2010年、2015年和2020年三期土地利用變化轉移矩陣和變化趨勢可知，耕地、林地和草地面積均減少，建設用地面積增加，耕地、林地和草地均向建設用地轉化。2035年歷史趨勢情景下，建設用地增加44.61 km2；隨著建設用地擴張，林地減少面積最大，減少35.74 km2，其次是耕地、草地，流域由建設用地、水域向林地擴張的趨勢仍然明顯。

2035年國土空間規劃情景是在歷史趨勢情景下的基礎上，結合《郴州市國民經濟和社會發展第十四個五年規劃和2035年遠景目標綱要》，控制其他地類向林地擴張，林地面積增加2.01 km2，水域面積略有擴張，草地面積略有減少，耕地面積減少幅度稍大，減少19.60 km2，建設用地面積增加至20.07 km2，說明流域城市面積繼續擴張，擴張趨勢主要是耕地轉化為建設用地，建設用地向周圍耕地擴張，耕地面積減少主要發生在汝城縣。根據流域地形概況，該地區地形起伏稍大，建設用地擴張主要分布在地勢較平緩的耕地。兩種情景的土地利用格局以及變化差異如圖7所示。

圖7 兩種情景下的土地利用差異圖Fig.7 Land use difference map under two scenarios

2.4 基于土地利用預測下流域污染負荷變化量

經模擬，歷史趨勢情景下氨氮負荷為355.51 t，總磷負荷為3 918.08 t；國土空間規劃情景下氨氮負荷為353.40 t，總磷負荷為3 863.47 t。兩種情景下污染空間分布特征基本一致，如圖8所示，其中國土空間規劃情景比歷史趨勢情景氨氮減少2.12 t，總磷減少54.6 t。

圖8 2035年兩種情景下面源污染特征Fig.8 Characteristics of source pollution under the two scenarios in 2035

從各子流域的氨氮負荷變化情況看，負荷減少較多的子流域集中在流域中部和西部；從各子流域總磷負荷變化情況看，負荷減少較多的子流域集中在4號、6號、7號和10號。

國土空間規劃情景與歷史趨勢情景相比，4 號、5 號、6 號子流域建設用地和耕地均有減少，4 號氨氮、總磷分別減少0.99、14.26 t，5號氨氮、總磷分別減少0.27、3.99 t，6號氨氮、總磷分別減少0.46、9.80 t；7 號、10 號子流域主要表現為林地面積增加，沒有向耕地的轉化，7 號氨氮、總磷分別減少0.16、7.15 t，10 號氨氮、總磷分別減少0.35、8.25 t。總體上看，污染負荷增加是由建設用地的快速擴張和林地的減少導致的；流域上游污染較輕，上游丘陵地勢坡度較高，林地覆蓋面積大，建設用地較少，因此產生的污染物較少；流域下游地勢低洼，坡度較小，林地覆蓋面積較小，耕地和建設用地沿湖分布，易產生面源污染。因此對于林地和耕地的限制轉換對于流域的污染負荷削減具有重要作用。具體情況，如圖9所示。

圖9 兩種情景下污染物變化量Fig.9 Changes in pollutants under two scenarios

3 結論

（1）創新性實現了PLUS 模型和SWAT 模型耦合模擬，并應用于不同土地利用格局下面源污染動態模擬。兩者結合實現了面源污染負荷和土地利用格局的動態解析，體現了土地利用格局優化對面源污染防治的實踐意義。

（2）根據東江湖流域面源污染長系列模擬結果，解析了其時空演變特征和變化趨勢。從時間尺度上看，污染輸出負荷集中在降水量較大的豐水期，年際變化表現為面源污染先增加后減少的趨勢；從空間尺度上看，氨氮、總磷負荷較高集中在流域西北部和中部徑流量較大以及農田較為分散的子流域。

（3）東江湖流域2010-2020年間土地利用變化趨勢表現為建設用地向周圍林地、耕地、草地擴張。若2035年繼續擴張建設用地，國土空間規劃情景比歷史趨勢情景下建設用地少增加18.03 km2，林地相對增加26.82 km2，耕地相對增加4.84 km2，整體而言國土空間規劃情景下城市發展步伐放緩，有利于生態環境保護。

（4）經SWAT面源污染模型模擬顯示，國土空間規劃情景下比歷史趨勢情景下氨氮負荷和總磷負荷分別減少2.12 t 和54.6 t，通過對建設用地發展的合理限制以及對林地保護的布局優化，證實調整土地利用格局對減少面源污染具有明顯的正向效益。