Multi-Scenario Simulation of Land Use Change in Longxi River Basin Based on PLUS Model

龙溪河流域(29°38′-30°53′N，106°49′-108°05′E)发源于重庆市梁平区明达镇，自东北向西南进入长寿区，流域面积约3 112 km²，是长江北岸一级支流(图 1)。地貌以低山丘陵为主，海拔为103~1 120 m。丘陵地带地势起伏较小，多由低矮的山丘和缓坡组成，是山地和平原间的过渡地带。平原部分则地势平坦，土壤肥沃，是主要的农业区。气候属中亚热带湿润季风气候区，年均气温16.8 ℃，年均降水量1 192.4 mm。土壤主要有水稻土、黄壤、石灰土和黄褐土。截至2020年底，流域人口约206.7万人。流域水土流失面积约占流域面积的40%，且以水力侵蚀中的面蚀和沟蚀为主。流域干流年均水质达Ⅲ类考核标准。

主要包括空间信息数据、社会经济数据、空间规划数据(表 1)。将所有栅格和矢量数据统一坐标系为：CGCS2000_3_Degree_GK_CM_108E，以减少图形偏差。

DEM数据。利用ArcGIS 10.2中Spatial Analyst Tools工具中的提取分析功能，以龙溪河流域边界为基准，通过掩膜提取工具从重庆市DEM图中提取出龙溪河流域的DEM图。利用Spatial Analyst Tools工具中的表面分析和邻域分析功能，从龙溪河流域DEM数据中提取坡度、坡向、高程变异系数、地表粗糙度和地形起伏度等空间数据。

土地利用类型数据。包含2000、2005、2010、2015、2020年的土地利用数据，土地利用类型分为旱地、水田、林地、草地、水域和建设用地6类。

土壤类型数据。以龙溪河流域边界为基准，利用ArcGIS 10.2对重庆市土壤类型空间分布图进行掩膜裁剪，提取龙溪河流域土壤类型数据。参照中国科学院的土壤类型分类标准，将龙溪河流域土壤类型分为水稻土、紫色土、中性紫色土、黄壤、石灰岩等10类。

气象数据。收集2000-2020年龙溪河流域周边气象站点数据，采用克里金插值法对离散气象数据进行空间插值。利用栅格计算器计算流域年均气温和年均降水量。

社会经济数据。人口密度和GDP数据来源于中国科学院资源环境科学数据平台，利用ArcGIS 10.2对其进行掩膜裁剪，提取龙溪河流域人口分布和GDP数据；利用ArcGIS 10.2空间分析工具，采用欧式距离法，获得到点、线、面状的距离数据。

空间规划数据。以龙溪河流域边界数据为基准，利用Date Management Tools模块中的Raster Processing命令对重庆市“三区三线”数据进行裁剪，提取龙溪河流域对应数据。

自然因素是影响土地利用变化的主要因素，其与社会经济因素交互作用后，因子解释力得到增强。自然因子影响土地利用变化发生的概率，如高程越高、坡度越大，立地条件越差，土地受扰动的程度越低，土地利用发生变化的概率越小。土壤和水热条件越好，土地利用发生变化的概率越大。社会经济因子才是驱动土地利用变化根本动力，经济发展水平越高、人口密度越大、到河网或道路(居民点)的距离越近，土地利用的驱动力越强，土地利用发生变化的概率越高。基于龙溪河流域自然地理条件与社会经济状况的实际，考虑距离对人为活动的影响，研究选取5项自然环境因子(高程、坡度、土壤类型、年均气温、年均降水量)和5项社会经济因子(人口密度、人均GDP、到河网距离、到居民点距离、到主要交通道路距离)。在ArcGIS 10.2软件中，将10个驱动因子统一处理为30 m×30 m的栅格数据，对到主要交通道路、居民点和河网的距离，利用ArcGIS 10.2中的欧氏距离工具进行分析，其中到居民点距离不是主要影响因素，因此图中没有展示(图 2)。

依据龙溪河流域土地利用变化规律，考虑政策对流域土地利用变化的影响，研究设定4种可能的土地利用变化情景：自然发展情景、耕地保护情景、生态保护情景和旱改水发展情景。

自然发展情景：土地利用变化主要由流域自然环境特点和社会经济发展状况所决定。本研究认为未来土地利用变化过程仍按历史演化趋势演进，展现出自然发展状态。

耕地保护情景：防止耕地“非农化”和“非粮化”，切实改进耕地占补平衡方式。本研究考虑加强地势平缓区耕地的保护利用，而陡坡区耕地则强调利用方式的多样性。

生态保护情景：严守生态保护红线，严禁建设占用生态保护用地和毁林开荒，加大生态修复力度，强化水源涵养与水土保持，禁止生态用地转换为耕地。

旱改水发展情景：聚焦占优补优，占水田补水田。本研究考虑地形坡度、交通区位、灌溉条件、产业发展等因素，实施旱地垦造水田整治工程，提高耕地设施配套程度与耕地质量。

将2010年和2020年土地利用数据输入PLUS模型的“Extract Land Expansion”模块中，将地类未发生变化的栅格单元value值设为0，发生变化的记为原地类代码，再将用地扩张数据和驱动因子数据输入LEAS模块，生成土地利用发展概率。本研究将随机森林数量设置为20，采样率为0.01，采样方式为随机采样，特征个数为10，线程为15(图 3)。

从图 4可以看出，影响旱地扩张的主要因素为坡度、年均气温和高程。坡度影响土壤侵蚀和地表水分的流动，坡度越陡越不利于水分的保持。气温影响水分蒸发，影响土壤水分含量。高程影响降水分布与农业生产布局。水田扩张受坡度和高程的影响，相比旱地，其年均降水量的贡献有所增加。年均气温、高程和坡度对林地的扩张影响较大，一定程度上决定林地的生长条件、分布和变化趋势。影响草地扩张的主要因素为人口密度、到河网距离和到主要交通道路距离，它们关联人类活动对土地类型的改变，影响水资源可得性和交通便利性。水域扩张的主要因素为土壤类型、坡度和年均气温，它们对水体的化学成分、水量和水循环产生重要影响。影响建设用地扩张的主要因素为人口密度。随着人口密度的增加，对住房、基础设施和经济活动的需求也会相应增长，进而导致建设用地需要进一步扩张。

1) 限制转换区设置。自然水体(河流、湖泊、水库等)自然状态下具有相对稳定性，不易发生改变，即水域的改变需要较大的投入。本研究利用ArcGIS 10.2，将历史土地利用数据重分类，把水域设置为0，表示不可转换，其他地类设置为1，表示可转换。

2) 邻域效应。反映周围环境对特定土地利用类型的影响，其强弱和方向决定变化的概率和方向。将30 m×30 m栅格数据的邻域设置为3。PLUS模型中递减阈值的衰减系数用于控制产生新斑块的难易，其值为0~1。衰减系数越大，越难以转换。本研究将衰减系数设置为0.9，以更好地模拟变化的实际情况^[26]。扩张系数用于调整模型生成新的土地利用斑块的能力，其值为0~1。数值越小，表示越难生成新的斑块，本研究将扩张系数设置为0.1。

3) 转换成本矩阵。该矩阵为衡量一种土地利用类型能否转化为另一种的重要指标。当允许转化时，成本设为1；不允许转化时，成本设为0。考虑到水域相对稳定，建设用地转化难度较大，本研究设定水域和建设用地为不可转化，其他用地类型为可相互转化(表 2)。

4) 邻域权重。取值范围为0~1。数值越大，邻域影响越大，即某种土地类型的扩张能力越强。根据各类用地以往变化规律确定邻域权重，计算公式为：

式中：W_i为第i类土地类型领域权重；TA_i为第i类土地利用扩张面积(hm²)；TA_min为各类土地利用扩张最小地类扩张面积(hm²)；TA_max为各类土地利用扩张最大地类扩张面积(hm²)^[27]。

根据公式(1)，研究区2010-2020年旱地、水田、林地、草地、水域和建设用地的领域权重分别为0、0.47、0.53、0.51、0.50、1.00。领域权重为0表示该土地利用类型扩张能力极弱，易被其他地类侵占，参考研究区2000-2020年土地利用变化趋势，将旱地的领域权重调整为0.10(表 3)。

5) 用地发展需求。利用2000年和2010年数据，通过马尔科夫链模块对研究区2020年各类土地利用发展需求进行预测。将预测数据输入CRAS模块中预测研究区2020年各类土地利用类型。从表 4可以看出，2020年预测的结果与实际数据较接近，这表明利用PLUS模型中的马尔科夫模型预测目标年份的土地利用栅格需求是相对可靠的。

利用Kappa系数评价PLUS模型土地利用模拟结果的精度，计算公式为：

式中：p_o为模拟与实际土地利用一致栅格数量占总体栅格数量比例；P_e为随机条件下期望的正确率；1代表理想情况下模拟正确的栅格比例。Kappa系数为-1~1。

从图 5可以看出，研究区2020年土地利用模拟结果和实际值较接近，两者一致性较高，Kappa系数为0.958，说明该模型能用于研究区未来土地利用变化模拟。使用PLUS模型及选取的10个驱动因子和参数，能较准确地模拟研究区土地利用的未来空间分布。

假设研究区2020-2030年间土地利用变化规律与2010-2020年保持一致，基于2000、2010、2020年土地利用现状资料，利用PLUS模型中的Markov Chain模块，对研究区自然发展情景下2030年土地利用变化趋势进行模拟(表 5和图 6a)。

从表 6可以看出，自然发展情景下建设用地快速扩张，其面积增加3 313.11 hm²，增长率达42.16%，而旱地和水田减少面积最多，分别为1 948.99、1017.26 hm²。

自然发展情景下2030年桂溪、黄沙和云台的旱地面积减少较多，分别为724.43、244.92、225.81 hm²，分别占旱地总减少面积的36.75%、12.42%、11.45%；水田面积减少主要集中于云集、长寿湖和新立，分别为113.76、100.49、77.80 hm²，分别占水田总减少面积的11.13%、9.84%、7.62%；林地面积减少以长寿湖、裴兴和云集为主，分别为41.95、34.77、26.94 hm²，分别占林地总减少面积的12.99%、10.76%、8.34%；草地面积减少主要集中在回龙、桂溪和三溪，减少面积分别为16.79、7.56、5.55 hm²。自然情景下不同区域水域面积变化差异较明显，但总体呈上升趋势，其中长寿湖水域面积增长最大，达11.12 hm²。桂溪、云台和黄沙的建设用地扩张显著，分别增加789.18、276.64、255.44 hm²，分别占总扩张面积的23.82%、8.35%、7.71%。

自然发展情景下，产生上述结果的原因在于城镇化进程进一步加速，大量农村人口逐渐向城镇集中，从而驱使建设用地的快速扩张，导致大量耕地(旱地和水田)、林草地和水域向建设用地转换，尤其是立地条件较好的区域。而且，在乡村振兴过程中，农村特色主导产业发展对建设用地的需求加剧，这也是导致研究区建设用地快速增长的主要原因之一。建设用地增加、旱地减少均主要发生在城镇建设用地迅速扩张区。水田减少、林草地减少主要发生在农业特色主导产业发展基础较为突出的区域。但是，自然发展情景仅是以往演化趋势延续的假设，是缺少政策作用下的理想演化状态，偏离政策作用后的实际情况，为此，研究区自然发展情景下的土地利用预测必须引入政策约束因素予以修正。

耕地保护情景中，基于粮食安全和基本农田保护政策考量，优先保护耕地，限制基本农田向其他土地利用类型转变。在考虑自然发展情景基础上，为限制耕地向非农用途的转换，参考龙溪河流域近年来耕地向非农建设用地转换的速率，将旱地和水田向建设用地的转移概率减少20%，得到2030年耕地保护情景下的模拟结果(表 7和图 6b)。

从表 8可以看出，耕地保护情景下建设用地扩张速度减缓，比自然发展情景减少292.73 hm²。旱地、水田和林地均呈不同程度地减少，其中旱地和水田减少最多，分别为1856.17、815.20 hm²，但旱地和水田面积相较于自然发展情景有明显增加，分别增加92.82、202.05 hm²，其中旱地增加主要集中在云台、海棠和桂溪，分别为102.42、94.66、79.39 hm²，水田增加主要集中于新立、云集，分别为26.89、21.30 hm²。

耕地保护情景下2030年桂溪、黄沙和邻封旱地面积显著减少，分别为645.04、268.41、236.90 hm²，分别占旱地总减少面积的34.34%、14.29%、12.61%；水田减少面积主要集中于云集、长寿湖和桂溪，分别为92.46、88.88、65.43 hm²，分别占水田总减少面积的11.28%、10.84%、7.98%；林地减少面积主要集中在长寿湖、裴兴和云集，分别为41.81、34.44、27.07 hm²，分别占林地总减少面积的12.95%、10.67%、8.39%；水域面积总体呈上升趋势，其中长寿湖水域面积增长最多，达11.25 hm²，其次为龙溪河流域，增长4.59 hm²；桂溪、邻封和黄沙建设用地扩张较显著，分别为701.05、318.38、279.12 hm²，分别占建设用地总扩张面积的23.22%、10.55%、9.25%。

耕地保护情境下，研究区建设用地扩张速度较自然发展情景明显放缓，尽管旱地和水田面积仍在减少，但减少的速度得到有效控制。目前，区域经济正处于由传统增长模式向高质量发展转变，建设对土地的需求逐渐降低，此外，占优补优、占水田补水田政策的实施，使建设占用耕地的面积不仅显著减少，且占用耕地的质量也较以往有所降低，从而促使土地利用结构进一步优化。当然，促使耕地减少的主要原因仍是城镇建设用地扩张及农村现代农业产业发展。基本农田保护红线的划定是限制耕地向其他地类转换的主要因素之一，在做好耕地保护的同时应根据不同区域的资源环境条件，尽可能优化土地利用结构，适当增加农村产业发展用地，为乡村振兴过程中城乡要素资源的流动提供用地保障。

生态保护情景中，限制生态保护红线内的土地类型向其他用途转换。在优先考虑自然发展情景的基础上，为加快生态恢复和保护，考虑到耕地因撂荒向林草地转换的速度，将旱地和水田向林草地的转移概率增加20%，且因立地条件限制和现有政策文件对林草地和水域保护的要求，目前建设占用林草地的可能性和概率均不高，故将林草地和水域向建设用地转移的概率减少30%，得到2030年生态保护情景下的模拟结果(表 9和图 7a)。

从表 10可以看出，受生态保护红线约束，对比2030年自然发展情景，生态保护情景下林地面积增加131.97 hm²，主要集中在竹山、屏锦和回龙，分别为35.56、33.97、13.72 hm²；草地面积增加4.35 hm²；水域面积增加6.78 hm²；建设用地增长面积减少206.61 hm²；旱地和水田均呈减少趋势，其减幅分别为1 869.93、1 026.40 hm²。

生态保护情景下2030年桂溪、黄沙和云台旱地面积显著减少，分别为714.61、237.93、212.23 hm²，分别占旱地总减少面积的37.76%、12.57%、11.21%；水田面积减少主要在云集、长寿湖和新立，分别为110.57、100.82、77.73 hm²，分别占水田总减少面积的10.72%、9.78%、7.54%；林地面积减少以长寿湖、裴兴和杠家为主，分别为37.24、35.13、26.13 hm²，分别占林地总减少面积的19.42%、18.32%、13.63%；草地面积减少主要集中在回龙、桂溪和三溪，分别为18.11、6.93、5.55 hm²；水域面积增长最多的是长寿湖，达11.70 hm²。建设用地扩张最多的是桂溪、云台和黄沙，分别为774.73、260.44、250.11 hm²，分别占总扩张面积的24.95%、8.39%、8.05%。

生态保护红线是限制生态用地向其他地类转换的关键因素之一，但是为满足城镇建设和乡村产业发展需要，建设用地扩张仍是导致旱地、水田、林草地等减少的主要原因。即便在生态保护情境下，仍要协调生态保护、耕地保护与经济发展之间的关系，做到“三线”动态协调共生。对立地条件相对较差的区域，林草地分布较多，应尽量减少其向建设用地转换，而对地势相对平缓的优质耕地集中分布区，应减少非农建设占用耕地。低丘缓坡区则是应成为非农建设占用耕地的区域，但应做好生态承载力评估，降低大面积开发对生态完整性的影响。

旱改水发展情景下，响应占优补优，占水田补水田，增加水田面积，优先考虑在自然发展情景基础上，基于立地条件和农户意愿，将旱地向水田的转移概率增加20%，水田向旱地转移的概率减少30%，得到2030年旱改水发展情景下的模拟结果(表 11和图 7b)。

从表 12可以看出，旱改水发展情景下，水田比自然发展情景增加164.94 hm²，主要集中在新立、云集和花桥，分别为19.21、17.70、16.64 hm²。旱地面积比自然发展情景减少121.83 hm²。林地面积减少282.95 hm²，草地面积减少69.56 hm²，水域面积增加6.44 hm²，三者变化趋势与自然发展情景差异不大。建设用地面积同自然发展情景基本一致。

旱改水发展情景下2030年桂溪、黄沙和云台旱地面积显著减少，分别为715.07、244.53、227.19 hm²，分别占旱地总减少面积的34.16%、11.68%、10.85%；水田减少面积主要在云集、长寿湖和桂溪，分别为96.06、90.38、66.95 hm²，分别占水田总减少面积的11.21%、10.55%、7.81%；林地面积减少以长寿湖、裴兴和云集为主，分别为42.40、35.41、27.01 hm²，分别占林地总减少面积的13.37%、11.16%、8.52%；草地面积减少主要集中在回龙、桂溪和三溪，分别为18.57、7.53、5.55 hm²；水域面积增加最大的是长寿湖，达11.15 hm²。桂溪、云台和黄沙的建设用地面积扩张明显，分别为778.19、275.05、255.51 hm²，分别占总扩张面积的23.81%、8.42%、7.82%。

研究区地势平缓，水田分布广泛，旱地与水田相嵌时有发生，且水资源丰沛，水利设施配套较齐全，做好耕地内部的优化调整，既能提高耕地质量，完善农田配套设施，又能提升耕地产能。为进一步提高耕地产能，满足占水田补水田的平衡要求，考虑立地条件和农户意愿，在研究区开展旱地改水田是极其合适的，不仅减轻占优补优的压力，更能提高水田田块的连片程度，为现代农业发展提供本底条件。但在旱改水实施过程中一定要考虑水源条件、设施配套条件、农户意愿、田块连片度等因素，不能一概而论。

考虑龙溪河流域的自然地理条件与社会经济状况，本研究选择了对流域斑块级土地利用作用较大的5项自然环境因子(高程、坡度、土壤类型、年均气温、年均降水量)和5项社会经济因子(人口密度、人均GDP、到河网距离、到居民点距离、到主要交通道路距离)。PLUS模型模拟结果发现，Kappa系数为0.958，模拟值与实际值一致性程度较高，表明因素选择合理且因素间的组合较好。这说明在开展流域土地利用变化模拟时，并不是因素选择得越多越好，选择最合适的关键影响因素才是最合理且符合实际的。

对于土地利用的情景设置，本研究重点考虑政策对龙溪河流域土地利用的作用，特别是与土地利用密切相关、对土地利用作用强度大的耕地保护、生态文明建设、经济高质量发展等方面，研究设置自然发展、耕地保护、生态保护和旱改水发展4大情景。PLUS模型模拟结果发现，不同情景间龙溪河流域主要土地利用类型差距较大，但并未出现数量级的差异。4种情景的横向对比发现，自然发展情景下旱地、水田、林地分别减少1 948.99、1 017.26、289.00 hm²，建设用地增加3 313.10 hm²；耕地保护情况下旱地、水田、林地分别减少1 856.17、815.20、288.47 hm²，建设用地增加3 020.38 hm²；生态保护情景下旱地、水田、林地分别减少1 869.93、1 026.40、157.03 hm²，建设用地增加3 106.49 hm²；旱改水发展情景下旱地、水田、林地分别减少2 071.11、852.24、282.95 hm²，建设用地增加3 269.41 hm²。综上，耕地保护情景下旱地和水田分别增加92.82 hm²和202.05 hm²，生态保护情景下林地分别增加131.97 hm²，旱改水发展情景下水田增加164.94 hm²。这说明模拟结果是非常适度的，符合政策的作用规律，即政策对流域土地利用的影响较大，但不会改变利于土地利用的根本格局，不会影响流域高质量发展。

但是，本研究选择的驱动因素和情景设置是否能用于类似其他流域的土地利用变化多情景模拟，仍需要进一步论证。因为不同流域的自然地理条件与社会经济状况有一定的差异，流域土地利用现状存在差异，未来土地利用的演化趋势受自然因素的影响与社会经济的驱动也有很大不同，基于PLUS模型的土地利用多情景模拟结果与实际情况间的符合程度会存在一定的差异。未来将进一步探索PLUS模型在其他流域的适用程度。

1) 使用构建的PLUS模型，对研究区2020年土地利用格局进行模拟，经检验，PLUS模型Kappa系数为0.958，模拟结果和实际值较接近，两者一致性较高。使用PLUS模型以及选取的10项驱动因子和参数，能较为准确地模拟研究区土地利用空间格局。

2) 自然发展情景下，流域土地利用演化仍按照以往的方式和速度进行，旱地、水田、林地分别减少1 948.99、1 017.26、289.00 hm²，建设用地增加3 313.10 hm²；耕地保护情景下，考虑到建设占用优质耕地的面积会减少，但并不以损害经济发展为代价，旱地、水田、林地分别减少1 856.17、815.20、288.47 hm²，建设用地增加3 020.38 hm²；生态保护情景下，基于陡坡耕地撂荒和对水域严格保护的实际，受林草地面积增加、耕地面积减少的驱动，旱地、水田、林地分别减少1 869.93、1 026.40、157.03 hm²，建设用地增加3 106.49 hm²；旱改水发展情景下，根据占优补优和耕地地力提升的需要，耕地内部利用方式的调整将成为未来土地利用变化方式之一，旱地、水田、林地分别减少2 071.11、852.24、282.95 hm²，建设用地增加3 269.41 hm²。

3) 与自然发展情景比较，耕地保护情景下旱地和水田分别增加92.82 hm²和202.05 hm²，生态保护情景下林地增加131.97 hm²，旱改水发展情景下水田增加164.94 hm²。耕地保护情景下建设用地增加量最小，旱地、水田减少量最小，林地减少量仅较生态保护情景多。考虑到耕地保护、经济发展和生态维持的协同需要，既要满足经济发展的需求，又不能造成耕地和林草地的大幅减少，做到“一要吃饭，二要建设，三要保护”，研究区未来土地利用格局调整推荐耕地保护情景下的模拟结果。

4) PLUS模型较好地剔除了流域土地利用驱动因素间的多重共线性与空间自相关问题，将不同政策作用下的非线性土地需求定量化，显著提高了模拟结果与实际情况的接近程度，适合用于流域土地利用变化的多情景模拟，能为流域国土空间格局优化提供科学依据。