鄭 昊，賀華翔，李海明，牛存穩(wěn)，王嘉浩，李昕陽

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

0 引言

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

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

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

1.2 研究方法

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

1.2.1 PLUS模型

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

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

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

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

表2 土地利用轉(zhuǎn)移矩陣Tab.2 Land use transition matrix

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

表3 鄰域權(quán)重參數(shù)Tab.3 Neighborhood weight parameters

1.2.2 SWAT模型

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

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

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

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

表4 SWAT模型構(gòu)建數(shù)據(jù)表Tab.4 SWAT model construction data table

1.3 模型精度驗證

1.3.1 PLUS模型精度驗證

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

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

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

1.3.2 SWAT模型率定驗證結(jié)果

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

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

表5 水量水質(zhì)率定驗證結(jié)果Tab.5 Verification results of water quantity and quality calibration

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

1.4 模型耦合及情景設計

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

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

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

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

（1）創(chuàng)新性實現(xiàn)了PLUS 模型和SWAT 模型耦合模擬，并應用于不同土地利用格局下面源污染動態(tài)模擬。兩者結(jié)合實現(xiàn)了面源污染負荷和土地利用格局的動態(tài)解析，體現(xiàn)了土地利用格局優(yōu)化對面源污染防治的實踐意義。

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

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

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