东江湖流域未来土地利用变化对面源污染负荷的影响评估

郑 昊，贺华翔，李海明，牛存稳，王嘉浩，李昕阳

（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，通过对建设用地发展的合理限制以及对林地保护的布局优化，证实调整土地利用格局对减少面源污染具有明显的正向效益。