利用CLUE-S模型对四川省未来10年土地利用变化的预测与分析

樊文杰, 戴晓爱*, 谢一茹, 高怡凡

(1.成都理工大学地球科学学院， 成都 610059； 2.北京师范大学环境演变与自然灾害教育部重点实验室， 北京 100875)

土地是人类赖以生存和发展的物质基础，是人类存在、生活和生产的源泉[1-2]。土地利用是人类根据其自身发展需求，对土地加以改造开发利用的过程[3]。随着人口数量的不断增加以及社会的进步，人地矛盾日益尖锐，因此对土地利用进行科学的模拟和规划至关重要。土地利用变化模型通过分析以往的土地利用数据来理解土地利用系统的机理[4]，并评估土地政策对自然生态和经济发展的影响，为土地利用提供决策支持，以达到人类社会的可持续发展。

土地利用变化模型层出不穷，其中使用较为广泛的模型包括：元胞自动机模型、马尔科夫链模型、灰色系统预测模型、CLUE-S(conversion of land use and its effects at small regional extent)模型等。中国学者已经应用不同的土地利用变化模型对不同尺度下的研究区土地利用做了丰富的研究[5-12]。其中， CLUE-S模型使用广泛，得到了非常好的验证。此外，与其他模型相比，CLUE-S模型能够对不同土地利用方式进行同步模拟。

素有天府之国美誉的四川省，其GDP总量在中国位列前茅，生态系统服务多样[13]。然而四川省作为人口大省，伴随着近年来人口的持续增加，该地区的生态保护压力、耕地安全压力和经济发展压力日益增加，经济发展与粮食安全、生态保护方面的协调发展至关重要。土地利用变化模拟能够在政策目标约束下实现对四川省未来土地利用变化空间格局的预测，预测结果对当前四川省面临的问题将具有重要的参考意义。当前关于四川省的土地利用变化模拟研究较为丰富，但大部分研究区域仅针对市、县一级[14-16]，未针对四川省全域进行探索，且模拟情景单一，多为对现状的分析[17-18]。

因此，现使用土地利用变化模型CLUE-S，对四川省全域范围进行土地利用时空格局动态演变模拟，探究未来数年的土地利用变化格局，探寻缓和人地矛盾、调控区域生态安全的方法，以期为四川省实现高质量、可持续发展提供科学参考。

1 研究区域

1.1 研究区域概况

四川省(26°03′～34°19′N、97°21′～108°12′E)位于中国西南部，辖21个市(州)，总人口约8 375万人，总面积48.67万km2；位于青藏高原地区向长江中下流平原的过渡地带，整体地势呈西高东低。川西地区为高原山地气候，以高原和山地为主，森林草地分布范围广，是主要的生态系统服务供给区域；川东地区为亚热带季风气候，以丘陵和盆地为主，城市化水平高、经济发达，是经济开发重点区域。四川省2015年各土地利用类型面积从大到小依次为：草地(17.53万km2)、林地(16.56万km2)、耕地(11.93万km2)、未利用地(1.89万km2)、建筑用地(0.47万km2)、水域(0.29万km2)，空间分布如图1所示。

图1 四川省2015年土地利用分布Fig.1 Land use distribution in Sichuan Province in 2015

1.2 数据来源与处理

本文中使用了中国科学院地理科学与资源研究所(http://www.resdc.cn)的2015年中国省级行政边界数据、2015年中国GDP空间分布公里网格数据集、全国1 km的数字高程模型和4期(2000年、2005年、2010年、2015年)土地利用数据；国家地球系统科学数据中心(http://www.geodata.cn)的2015年中国1 km分辨率平均气温数据集和平均降水量数据集；其他空间数据主要包括：Global Human Settlement(https://ghsl.jrc.ec.europa.eu/ghs_pop2019.php)的全球人口空间分布数据；Weiss等[19]提供的城市可达性数据；OpenStreetMap(http://download.geofabrik.de)的河流、湖泊和公路网数据；国际土壤参考资料中心(https://data.isric.org)的土壤有机碳密度和土壤质地数据。社会经济数据来自中国国家统计局和四川省十四五规划相关文件。研究中使用到的坡度数据是使用ArcGIS中Slope工具将高程数据作为输入生成的；可达性栅格数据是使用ArcGIS中Euclidean Distance将对应的矢量文件作为输入得到的。所有使用到的栅格数据经重采样处理后统一成1 km×1 km的分辨率。

2 研究方法

2.1 模型

研究中采用由荷兰瓦赫宁根大学的Verburg等[20]开发的CLUE-S模型。该模型的前身是CLUE模型，相比于CLUE模型而言更适合省级区域范围的研究。CLUE-S模型是以空间位置适宜性分析为基础，结合土地系统时空动态竞争和相互作用，来模拟土地利用变化。

CLUE-S模型可以分为需求量模块和空间分配模块[20-22]。需求量模块根据需求情景在整体层面上计算不同土地利用类型的总需求量，该模块通常是通线性规划或者马尔可夫过程等方法处理后作为CLUE-S模型的输入。空间分配模块将上述土地需求量反映为空间上某一位置土地利用类型的变化。CLUE-S模型需要土地政策及限制、土地利用类型转换、土地利用类型需求量和土地位置属性四类信息。土地政策及限制表明根据相关法律政策哪些区域的土地利用类型不允许转换，比如永久基本农田、国家级自然保护区、国家森林公园等。土地利用类型转换有转换规则和转换弹性两项，前者表明不同的各种土地利用类型能否相互转换，后者反映土地利用类型维持自身状态的能力。土地利用需求量是第一个模块得到的结果。土地位置属性信息用于将一系列相关因子作为输入使用逻辑回归计算出某个位置最佳的土地利用类型。

2.2 发展情景设计

基于研究区域自身特点和其发展战略，综合考虑生态效益、经济效益，设计了三种情景。第1种是生态效益型情景，该情景在满足政府预期经济增长速度的前提下，能够达到生态效益最大化。第2种是经济效益型情景，该情景以不减少生态效益为基础，最大化经济效益。第3种是综合效益性情景，该情景为情景1和情景2的折中实现，该情景首先要求生态效益增量为情景1下生态效益增加量的50%，然后最大化经济效益。

选取线性规划来计算土地利用需求量。根据模拟情景，可以确定目标函数[式(1)]、面积约束条件[式(2)～式(4)]、价值量约束条件[式(5)]为

X=argmaxX(KX)

(1)

s.t.Xi≥Xmin

(2)

Xj≤Xmax

(3)

(4)

(5)

式中：X为通过线性规划求解得到的各土地利用类型面积；Xi表示第i种土地利用类型的所占面积,km2；Xmin是第i种土地利用类型面积的最小值,km2；Xmax是第i种土地利用类型面积的最大值,km2；Varea是四川省的总面积,km2。

在情景1中，式(1)中K表示各土地利用类型的生态效益系数，式(5)中r表示到模拟终止年份GDP的增长率，Ci表示各土地利用类型的经济效益系数。在情景2和情景3中，式(1)中K表示经济效益系数，Ci表示生态效益系数。在情景2中，式(5)中r=0；而在情景3中，式(5)中r的值等于情景1中生态效益增长率的50%。

生态效益系数法、经济效益系数法[11,23]用来度量土地的生态效益和经济效益。经济效益系数是根据国家统计局公布的研究区域农业总产值、林业总产值、牧业总产值和渔业总产值以及研究区域全年GDP除以每种土地利用类型相应的面积计算得到的。此外，四川省内的未利用地主要是由裸岩石质地的土地组成，所以将未利用地的经济效益系数设置为0。生态效益系数根据四川省生态系统服务价值量[24]得到。计算得到的系数见表1。

表1 土地利用类型的生态效益系数和经济效益系数Table 1 Weight of ecosystem services value and economic benefit coefficients by land use types

根据四川省“十四五”规划文件，本文中设定在所有的模拟情景中耕地面积不再减少，草地的面积不少于当前草地面积的95%，未利用地的面积不少于当前未利用地面积的90%。水域的面积较短时期内往往难以发生骤变，因此采用GM(1,1)模型预测其面积。根据上述方程使用MATLAB软件中的线性规划相关函数求解得到2030年各土地利用类型的面积预测值。由于CLUE-S模型需要输入逐年的面积预测值，故依据2015年和2030年的数据使用线性插值法得到中间各年面积的预测值。

2.3 模型参数设定

2.3.1 Logistic回归分析

在CLUE-S模型中土地利用转换往往会发生在“适宜性”最佳的位置，这种适宜性依赖于研究区域的自然地理环境和社会经济发展状况。为了量化这种适宜性，根据研究区特性，选取了人均GDP、人口密度、气温、降水、高程、坡度等13种驱动因子，使用Logistic回归分析得到某一栅格单元转变为某种土地利用类型适宜性值，即

(6)

式(6)中：Pi是当栅格单元是第i种土地利用类型时的适宜性，X是驱动因子在该栅格单元的值。逐步回归可以使得与第i种土地利用类型相关性较小的驱动因子从回归方程种剔除出去，也就是将驱动因子对应的β系数设置为0。由于ROC(receiver operation characteristic)方法具有不受样本不均衡影响等一系列优点，故使用ROC方法对回归结果的好坏做出评价。ROC曲线下面积称为AUC，AUC的值介于0.5～1。AUC的值越大，说明拟合结果越好，驱动因子对栅格单元土地利用类型的解释性就越强。

由于要对每种土地利用类型进行回归分析，故把土地利用类型根据中国科学院土地利用分类系统将土地利用类型重分类成六大类[25]，再次使用重分类工具制作对应的二值图。CLUE-S模型没有内置Logistic回归分析功能，所以先使用ArcGIS中Raster to ASCII工具把13种驱动因子对应的栅格数据以及土地利用二值图像转换成ASCII格式，再使用模型自带的FileConvert工具把ASCII格式的数据转换成统计分析软件SPSS接受的格式，最后运行SPSS中的二元Logistic回归分析功能，得到的结果如表2所示。表2表明，所选驱动因子对耕地、草地和未利用地解释性非常好，其AUC值分别为0.900 512、0.872 227和0.845 733；对建设用地和林地的解释性较好，其AUC值分别为0.789 206和0.705 876。

表2 Logistic回归的分析结果与检验结果Table 2 Results of logistic regression and ROC test

2.3.2 转换规则和转换弹性

通过使用2000年、2005年、2010年、2015年四期四川省土地利用栅格数据进行统计分析，可以得到近些年来土地利用类型转换情况，然后使用上述数据作为土地利用转换矩阵，如表3所示。土地利用转换矩阵是一个二维数字矩阵，矩阵中的取值有0和1两种。如果第i行第j列的值为0，则表示的是不允许第i行所代表的土地利用类型向第j列代表的土地利用类型发生转换，否则可以发生转换。

表3 土地利用转换矩阵Table 3 Matrix of land use conversion

同样使用上述四期数据进行统计分析，可以得到近些年来各土地利用类型向其他土地利用类型转出的比率。比率越小，说明该土地利用类型维持类型不变的能力越强，弹性也就越大。结合转出比率数据以及专家经验，得到转出弹性系数设置结果，如表4所示。

表4 土地利用类型转换弹性Table 4 Conversion elasticity of land use types

2.4 模型验证

为了验证CLUE-S模型运行结果的准确性，在正式模拟之前，首先以四川省2010年土地利用图为起点，模拟四川省2015年土地利用情况。将模拟结果与真实结果作对比。如果相同或高度相似，则说明模拟效果好，反之亦然。在评价相似性时，采用Kappa系数法。具体而言，Kappa系数法需要输入两幅影像，综合考虑像元数量相似性以及位置相似性，来判断这一对影像是否具有空间分布一致性。本文中共采用以下三种Kappa系数：标准Kappa、反映数量相似性的KHis、反映位置相似性的KLoc，这些系数的计算公式为

(7)

(8)

(9)

式中：po为观测图和模拟图中像元值相同的像元个数与总像元数的比值，即总体精度；pe为服从观测图数据分布下的一致性期望概率；pmax为服从观测图数据分布下的一致性最大概率。标准Kappa的值等于KHis与KLoc的乘积。通常情况下Kappa值介于0～1，0表示观测图与模拟图之间具有极低的一致性，1表示观测图与模拟图之间具有高度一致性。

本文中使用了根据每一类别像元总数计算得到的KHis、根据每一类别像元空间分布计算得到的KLoc和标准Kappa。将上述两幅影像作为 Map Comparison Kit 软件的输入，可以得到KHis系数的值为0.999 77、KLoc系数的值为0.896 76、Kappa系数的值为0.896 55，表明CLUE-S模型参数的设定能较好地反映和预测研究区域土地利用情况。

3 结果与分析

3.1 不同情景下土地需求结果分析

通过线性规划计算出3种模拟情景下的土地利用需求面积，如表5所示。在3种情景下耕地的需求面积与2015年耕地面积基本相同。林地面积在生态效益情景下最大，经济效益情景下最小，综合效益型居于两者之间，这也体现出林地具有较高的生态效益，并且三种情景中林地面积与现状值相比都略有增长。从表5种可以看出草地、水域和未利用地的面积在3种情景中分别相同，这也符合本文的需求设定。建设用地的面积在经济效益型情景中最大，在生态效益型情景中最小，综合效益型居于两种情景之间，这也体现出建设用地具有较高的经济效益。

3.2 不同情景下土地模拟结果分析

运行CLUE-S模型可以得到不同模拟情景下的四川省2030年不同土地利用类型面积的大小及其空间分布，如表6和图2所示。通过比较表5和表6，发现两者对应数值的大小几乎保持一致，说明模型结果能够较为理想地反映设计的情景需求。

表5 不同模拟情景预测的土地利用需求面积Table 5 Area of land use demand predicted by different simulation scenarios

表6 不同模拟情景模拟的土地利用需求面积Table 6 Area of land use demand simulated by different simulation scenarios

图2 不同情景下土地利用分布Fig.2 Land use distribution under different scenarios

根据图2，从整体上看，耕地主要分布在成都平原、四川省的东北部地区、四川省的南部地区、凉山州和攀枝花市。林地分布在除了成都平原及其周边地区外的地区。草地和未利用地集中分布在甘孜州和阿坝州。水域分布在四川省的东北部以及凉山州和攀枝花市。建设用地分布在成都周边的城市、四川省的东北部以及攀西地区。这与当前的土地利用分布现状基本相同。

生态效益型情景模拟结果如图2(a)所示。该情景下的生态效益为108 405.5亿元，经济效益为71 913.6亿元，与其他两种情景相比其生态效益达到最大值。从地域分布来看：攀西地区的耕地面积减少，川东北地区北部耕地面积也减少，川西平原西南部地区的耕地面积增加；川西平原西南部地区的林地减少，攀西地区的林地会有所增加；四川盆地北部边沿和西部边沿的草地会减少，攀西地区西部的草地也会减少，川西北地区草地有所增加；四川盆地内水域面积会减少，攀西地区和川东北地区水域面积会增加；新的建设用地主要集中于攀西地区和川东北地区；未利用地主要在川西北地区内发生变迁。经济效益型情景模拟结果如图2(b)所示。该情景下的生态效益为107 857.9亿元，经济效益为78 346.3亿元，与其他两种情景相比其经济效益达到最大值。综合效益性情景模拟结果如图2-c所示。该情景下的生态效益为108 068.2亿元，经济效益为75 480.7亿元，其生态效益和经济效益都处于中间水平。在情景2和情景3中，各土地利用类型的动态变化与情景1基本保持一致。

三种情景中土地利用分布在成都平原经济区、川东北经济区和攀西经济区差异最为显著。与情景2相比，情景1在成都平原经济区有更多的林地(+484 km2)，在川东北经济区和攀西经济区建设用地面积会减少(分别为-371 km2和-612 km2),体现出了该情景对于生态效益增加的需求；与情景3相比，情景2在成都平原经济区林地面积会减少(-220 km2)，在川东北经济区和攀西经济区有更多的建设用地(分别为+151 km2和+276 km2),体现出了该情景对于经济效益增加的需求；与情景1相比，情景3在成都平原经济区的林地也会减少(-264 km2)，在川东北经济区和攀西经济区建设用地面积会增加(分别为+220 km2和+336 km2)。

4 结论

根据研究区的发展设定了3种发展情景，使用CLUE-S模型对四川省2030年土地利用变化进行了模拟研究。

(1)研究结果表明生态效益型情景能够最大化生态系统服务价值，生态系统服务价值总量为108 405.5亿元；经济效益型情景能够维持生态系统服务水平现状，还可以最大程度地满足经济发展对土地的需求，经济效益为78 346.3亿元；综合效益型情景综合考虑生态和经济发展，实现了两者的平衡，生态系统服务价值总量和经济效益分别为108 068.2亿元和75 480.7亿元。

(2)综合考虑3种情景，从整体的发展趋势上来看，与现状相比到2030年时，林地增加部分主要分布在攀西地区；水域增加部分主要分布在攀西地区；建筑用地增加部分主要分布在攀西地区和四川盆地。

(3)研究结果可以为四川省土地利用规划和管理提供参考，本文的方法也可以应用于我国的其他省份土地利用研究。在未来的研究中，可以通过收集更加丰富的土地利用数据，完善模型中的转换规则；使用不同的方法计算预测年份的各土地利用类型面积大小；将模型中使用的逻辑回归替换为随机森林来完成“适宜性”的计算。