基于CLUE-S模型的煤矿城市土地利用变化模拟

赵明松 徐少杰 邓 良 刘斌寅 王世航 吴运金

(1.安徽理工大学空间信息与测绘工程学院， 淮南 232001；2.安徽理工大学矿山采动灾害空天地协同监测与预警安徽省教育厅重点实验室， 淮南 232001；3.安徽省地质测绘技术院， 合肥 230022； 4.生态环境部南京环境科学研究所， 南京 210042;5.国家环境保护土壤环境管理与污染控制重点实验室， 南京 210042)

0 引言

土地利用/土地覆盖变化(Land use and land cover change, LUCC)改变了区域的自然景观、物质循环、能量流动以及各种生态过程，影响着全球环境和区域生态经济系统的可持续发展，是当前土地变化学科和景观生态学研究的热点和前沿问题[1-3]。近年来，国内外学者对土地利用的研究主要集中在时空格局变化、驱动机制分析[4-5]、模拟预测[6-8]和生态环境效应[9-10]等方面。不同空间尺度上，模拟LUCC的模型较多，主要分为数量变化与空间变化的土地利用变化动态模拟模型及其混合模型。常见模型主要有Markov模型[11]、CLUE-S模型[11-14]、CA及改进模型[15-17]、系统动力学模型[18-19]和FLUS(Future land use simulation)模型[20]等。不同模型在LUCC模拟时具有各自的独特优势，其中，CLUE-S(Conversion of land use and its effects at small regional extent)模型因适用于小范围内的土地利用数量和空间位置变化的模拟与预测，且能够对土地利用变化与其社会、经济、技术、政策及自然环境等驱动因子相互关系进行定量分析，应用较广泛。

在CLUE-S模型中，运用Logistic逐步回归分析方法来计算每一种地类在区域内每个像元出现的概率[21]，分析各土地利用类型的空间分布与驱动因素之间的关系。传统的Logistic回归忽略了不同土地利用之间存在的空间依赖关系，即空间自相关[22]。Logistic回归模型假定土地利用格局的均质性，忽略了空间自相关对参与建模的各个变量之间相互作用的影响，BESAG[23]在传统Logistic回归模型基础上引入空间自相关空间权重，构建Autologistic回归模型。引入空间自相关因子，可以有效避免拟合残差的强自相关，提高模型精度。段增强等[24]对CLUE-S模型进行了改进，引入了动态计算的邻域分析因子，对北京市海淀区1991—2001年土地利用变化进行了多方案模拟。吴桂平等[25]通过将某一点某土地利用类型出现的概率定义为各驱动因子和一个空间自相关的虚拟变量的形式，解决了统计分析中空间自相关效应的影响。土地利用空间格局变化是引起土壤质量变化的主要原因，其中建设用地扩张不可避免会影响土壤质量尤其是耕地资源损失[26]。同时土壤质量也影响土地利用方式，如“永久基本农田”等高土壤质量的土地禁止转变为其他类型用地[27]。因此，在土地利用变化模拟中需要考虑土壤质量或土壤属性的影响。

两淮矿区是全国13个大型煤电基地之一，煤炭资源丰富，主要包括安徽省中北部淮河两岸的淮南和淮北两大矿区。淮南矿区主要包括淮南市区及辖县凤台县，矿产资源丰富。煤炭产业为区域经济的持续发展和人民生活水平的改善做出巨大贡献的同时，也带来一系列生态环境问题。因长期的采煤活动产生了大面积的地下采空区，引发土地沉陷，使得大量高产优质耕地、村庄等地表塌陷损毁[3, 28]。虽然已有研究利用GIS和RS技术分析淮南市土地利用程度变化过程的时空演化特征[29]，但是对研究区的土地利用变化模拟并未涉及。因此，本文选取淮南市为研究区，运用GIS和RS等技术分析1985—2016年研究区土地利用景观格局变化特征，并采用CLUE-S模型模拟未来土地利用空间格局，为区域国土空间规划管理、生态环境建设、塌陷区治理等提供数据基础。

1 研究区概况与数据来源

1.1 研究区概况

选择淮南市为研究区(32°32′45″～33°0′24″N，116°21′21″～117°11′59″E，图1)，主要包括淮南市辖区和凤台县，总面积为2 580.02 km2。研究区属暖温带和亚热带的过渡地带，年平均气温16.6℃，全年降水量893.4 mm。以淮河为界，研究区包含两种不同的地貌类型，淮河以南属于江淮丘陵的一部分，海拔17～75 m；淮河以北地势平坦，为河间浅洼平原，海拔20～24 m。研究区煤矿远景储量444亿t，探明储量180亿t，占安徽省的70%，占华东地区的32%。截至2017年研究区采煤塌陷区面积为278.60 km2，约35%的面积塌陷积水形成巨大湖泊。

图1 研究区位置Fig.1 Location of study area

1.2 数据来源与处理

(1)土地利用数据。1985、1995、2005年研究区土地利用数据来源于长江三角洲科学数据共享平台(http:∥nnu.geodata.cn)。2016年土地利用图，利用Landsat8 OLI影像解译获得。影像来源于地理空间数据云(http:∥www.gscloud.cn)，轨道号分别为121/38、121/37和122/37，日期为2016年7月25日、7月25日和9月2日。利用ENVI 5.3软件对OLI数据进行预处理，采用最大似然方法进行监督分类，Kappa系数为0.95。

(2)土壤数据。安徽省1∶500 000土壤图，来源于全国第二次土壤普查；典型土壤类型理化性质来源于2010—2011年土壤调查获取的206个典型样点数据集[30]，包括表层(0～20 cm)土壤有机质含量(SOM)、全氮含量(TN)、全磷含量(TP)和全钾含量(TK)数据。

(3)数字高程模型数据(SRTM DEM)来源于地理空间数据云(http:∥www.gscloud.cn)，空间分辨率90 m。在ArcGIS 10.2软件中，对DEM进行裁剪，并提取坡度和坡向。

(4)中国GDP空间分布公里网格数据集(单位：元/km2)和中国人口空间分布公里网格数据集(单位：人/km2)，来自于中国科学院资源环境科学数据中心数据注册与出版系统(http:∥www.resdc.cn/DOI)。

(5)年均温(MAT)和年均降水量(MAP)，来自中国农业科学院农业资源与农业区划研究所中国生态环境背景层面建造项目完成的栅格数据，空间分辨率1 km。

(6)研究区统计年鉴(2005—2016年)电子版，来源于研究区统计局(http:∥tjj.huainan.gov.cn/content/column/15642315)。

土壤数据、年均温、年均降水量等均为整个安徽省范围，中国GDP空间分布公里网格数据集和中国人口空间分布公里网格数据集为全国范围，重采样后裁剪研究区数据。利用ArcGIS 10.2软件中ArcToolbox中的Euclidean Distance工具制作与铁路、高速公路、国道省道、河流、城镇的空间距离图。

2 研究方法

2.1 土地利用动态度

选择单一土地利用动态度和综合土地利用动态度分析研究区土地利用变化的速度。

(1)单一土地利用动态度研究一定时间范围内某种土地利用类型的数量变化情况，计算公式为[31]

(1)

式中K——研究时段内某一土地利用类型动态度

Ua——研究期初某一土地利用类型面积

Ub——研究期末某一土地利用类型面积

T——研究时时段T定义为年时，K就是该研究区域某种土地利用类型年变化率。

(2)综合土地利用动态度描述区域土地利用变化的总体速度，用于反映变化的剧烈程度，计算公式为[32]

(2)

式中Si——研究期初第i类土地利用类型总面积

ΔSi-j——研究开始至结束时段内，第i类土地利用类型转换为其他土地利用类型的面积总和

S——与时段T对应的研究区土地利用变化速率,即综合土地利用动态度

2.2 土地利用变化重要性指数

土地利用变化重要性指数(Ci)反映区域主要土地利用变化的方向，筛选出土地利用变化的主要类型[33]。计算公式为

Ci=Ai/A×100%

(3)

(4)

式中Ci——第i类土地利用变化重要性指数，取值0～100%

Ai——第i类土地变化面积，km2

A——该区域各类土地变化面积之和，km2

Ci越大，说明第i类土地利用变化越占主导。利用ArcGIS软件计算1985—1995年、1995—2005年和2005—2016年的土地利用转移矩阵，计算土地利用变化重要性指数，识别主要的土地利用变化类型。

2.3 CLUE-S模型

CLUE-S模型是VERBURG等[21]在CLUE模型的基础上开发的，是一种在较小尺度上模拟土地利用变化及其环境效应的模型，包括空间分析模块和非空间分析模块。非空间分析模块独立于模型外，另需运用Markov模型进行土地利用需求预测，然后作为参数代入模型[34]。CLUE-S模型输入主要包括：土地利用转移规则、限制区域、土地利用需求量、空间分析。

2.3.1空间分析

在CLUE-S模型中，运用Logistic逐步回归分析方法来计算每一种地类在区域内每个像元出现的概率[21]，分析各土地利用类型的空间分布与驱动因素之间的关系。Logistic回归模型忽略了空间自相关对参与建模的各个变量之间相互作用的影响。本文选用BESAG[23]在传统Logistic回归模型基础之上引入空间自相关性构建的Autologistic回归模型，一般表达式为

(5)

式中Pi——每个空间单元为土地利用类型i的概率

β0——常数项β——向量X的系数

X——由一系列驱动因素构成的向量

yi——事件状态，为二值变量

wij——空间权重

r——事件状态yi的系数

本研究借助邻域因子构建空间权重矩阵，基本理念是每一个土地单元的状态是否发生变化，不仅取决于其自身所受到的影响，也同时受到一定周边范围内的用地状态的影响。空间权重设定为[35]

(6)

式中Sjn——栅格j周围每个邻域栅格与j的相似性，用地类型相同为1，否则为0

N——邻域栅格j的邻域总数

参照文献[36-37]，本研究选取16个土地利用变化的驱动因子，包括：DEM、坡度(Slope)、坡向(Aspect)、与铁路距离(Railway_d)、与高速公路距离(Freeway_d)、与国道省道距离(Road_d)、与河流距离(River_d)、与城镇距离(Town_d)、土壤有机质含量(SOM)、全氮含量(TN)、全磷含量(TP)、全钾含量(TK)、GDP、人口数量(POP)、年均温(MAT)、年均降水量(MAP)，并对驱动因子进行归一化处理(图2)。

2.3.2空间分配

空间分配是根据总概率对土地利用需求进行多次迭代分配的过程。具体公式[38-39]为

Ti,u=Pi,u+Eu+Iu

(7)

式中Ti,u——地类u在i栅格单元上的总概率

Pi,u——栅格单元i对于地类u的适宜性概率(通过Logistic回归方程求得)

Eu——设定的地类u的转移弹性

Iu——土地利用类型u的迭代变量值(所有地类均相同)

2.3.3模型结果检验

利用Kappa系数对CLUE-S的模拟结果进行精度检验。

利用SPSS 22软件进行Logistic回归分析，根据方程的系数β分析各驱动因子对不同土地利用类型的影响。ROC(Relative operating characteristics)方法被用于检验回归方程的解释能力，ROC值为ROC曲线以下的面积，其值介于0.5～1之间，越小解释能力越弱[40-41]。一般认为，当ROC值大于0.7时，计算地类的概率分布与真实地类分布有较好的一致性，驱动因子解释能力较好[42]。CLUE-S模拟在IDRISI 17.0软件中进行。

图2 各驱动因子栅格图Fig.2 Raster map of driving factors

3 结果与分析

3.1 土地利用动态度

图3 研究区1985—2016年土地利用类型Fig.3 Land use in Huainan from 1985 to 2016

研究区土地利用以耕地为主(图3)，4个时期均占总面积64%以上，最高达75.71%，其他地类面积比例由高到低依次为建设用地、水体、林地和草地。1985—2016年间，耕地面积减少299.93 km2，面积百分比减少11.62个百分点；建设用地和水体面积分别增加205.76、110.6 km2，面积百分比增加7.98个百分点和4.29个百分点；林地和草地面积分别减少2.41、14.02 km2，面积百分比减少不足1个百分点。水体增加主要发生在淮南市区的潘集区和凤台县。

从3个时段土地利用动态度来看(表1)，3个时段中，2005—2016年间各种土地利用变化速率最快，各土地利用动态度和综合土地利用动态度最高。3个时段研究区的综合土地利用动态度分别为0.05%、0.79%和13.46%。1985—1995年和1995—2005年两个时段土地利用变化缓慢，仅存在耕地面积的减少和建设用地面积的增加，其他土地利用类型变化不明显。2005—2016年，耕地面积减小270.71 km2，水体和建设用地面积分别增加92.17、195.41 km2，其中建设用地土地利用变化较为活跃，年变化率为5.19%。30年间，研究区耕地面积呈不断减少趋势，2005—2016年耕地年变化率为-1.28%，分别是1985—1995年和1995—2005年的25.6、12.8倍。水体从1995年面积开始增加，2005—2016年动态度达到3.4%。

表1 1985—2016年研究区土地利用动态度Tab.1 Land use dynamic index in Huainan from 1985 to 2016

3.2 土地利用类型空间变化

图4为研究区1985—2016年土地利用转移桑基(Sankey)图，表示不同时段各土地利用转入和转出情况。1985—1995年，土地利用变化不明显，仅存在少量耕地向建设用地转移。1995—2005年，土地转移主要发生在耕地、水体和建设用地之间。该时期，耕地转移为建设用地和水体面积为19.22、18.37 km2，分别占耕地转出面积的50.42%和48.19%。建设用地的转入面积和转出面积相差不大，主要转出为耕地，面积为18.14 km2，少量转移为水体，面积为0.39 km2。2005—2016年，建设用地类型的转移较为频繁，均有不同程度的转入和转出，未变化面积为235.46 km2，转入面积为302.18 km2，转出面积为106.78 km2。建设用地由耕地转入面积为277.82 km2，是建设用地未变动面积的1.18倍；水体转入面积为134.73 km2，其中93.44%由耕地转移而来，主要分布在研究区西部和北部的采煤塌陷区(图5a)，5.63%由建设用地转移而来。耕地面积的增加量主要来自于建设用地和水体，转入面积分别为98.90、35.04 km2。综合3个时段，水体转入面积依次为0、18.76、134.73 km2，且由耕地转入占0、97.92%、93.44%；水体转出面积依次为0、0.35、42.56 km2，其中2005—2016年间有82.33%流入耕地。

图4 1985—2016年土地利用类型转移Sankey示意图Fig.4 Sankey diagram of land use type conversion from 1985 to 2016

2005—2016年间研究区土地变化类型较多，变化速率较快，根据该时段土地利用变化重要性指数(Ci)，识别出该时段主要的土地利用变化类型。结果表明：2005—2016年间，共有20种土地变化类型，其中耕地→建设用地类型的Ci为47.66%，耕地→水体类型的Ci为21.60%，故耕地向建设用地和水体的转换为主要土地利用变化类型。

由图4可知，耕地→建设用地土地变化通过占用耕地，建设用地得到扩张。多数情况下建设用地的扩张与地区经济规模密切相关，研究区GDP从2005年的267.15亿元增加至2016年的821.4亿元，区域建设用地扩张除了和经济相关，可能因煤炭开采村庄搬迁，使得新建设区产生。1995—2016年研究区煤矿产业发展迅速，尤其在2005年以来煤矿开采的工业增加值和产量持续增长(表2)。2005—2016年研究区工业增加值由126.87亿元增加到347.77亿元，其中煤矿开采行业的工业增加值在总增加值的占比达60%以上，2013年后有所下降。2010—2013年每年新开工采煤沉陷区村庄搬迁安置点10个以上。2015年煤沉陷区搬迁安置房开工建设4 900套，搬迁入住4 800户；2016年采煤沉陷区搬迁安置房开工建设4 500套，搬迁入住4 400户。

图5 研究区2005—2016年土地利用变化图谱Fig.5 Land use change of Huainan during 2005—2016

表2 1995—2016年研究区煤矿开采行业生产概况Tab.2 Production of coal mining industry in Huainan from 1995 to 2016

耕地→水体土地变化主要发生在凤台县和市区的潘集区(图5a)，总面积达125.89 km2。由于煤矿开采导致了大面积地表沉陷，截至2017年，研究区采煤沉陷区总面积达278.60 km2，且矿区潜水集聚形成大片水体。新增水体集中分布在潘谢矿区内(图5b)，潘谢矿区为淮南的主要矿区，是典型的高潜水位矿区，塌陷后积水严重。其中，凤台县的岗河和西淝河一带，塌陷水体在顾桥矿、张集矿和新集矿内均有分布，丁集矿内也存在；潘集区的泥河一带，采煤塌陷地大面积出现(图5a)，塌陷水体在潘一矿、潘二矿、潘三矿和潘北矿均分布。

3.3 土地利用变化模拟与预测

3.3.1最佳模拟尺度选择

选用50 m×50 m、100 m×100 m、150 m×150 m、200 m×200 m和250 m×250 m共5种模拟尺度，计算土地利用类型和驱动因子的Logistic回归方程。表3为不同模拟尺度下的ROC，结果表明采用150 m×150 m网格的各土地利用类型回归方程的ROC最大，因此确定该尺度为最佳模拟尺度。

表3 不同尺度各地类ROCTab.3 ROC of different scales

3.3.2Logistic回归分析

表4为研究区各土地利用类型Logistic回归分析结果，L和LS代表未加入和加入土壤因子的普通Logistic回归模型，AL和ALS代表未加入和加入土壤因子的Autologistic回归模型。结果表明，土壤质量因子和空间自相关性因子，对不同土地利用类型Logistics回归建模的贡献率不同。

当普通Logistics回归建模(L)中加入土壤质量因子后(LS)，除建设用地外，其余土地利用类型建模的ROC有不同程度的增加。如，耕地的ROC由0.617增加到0.680；林地的ROC由0.807增加到0.961；建设用地的ROC没有明显变化，为0.708和0.709。普通Logistics回归建模(L)中加入空间自相关性因子后(AL)，仅有耕地和建设用地建模的ROC有明显增加。耕地的L和AL模型的ROC由0.617增加到0.722；建设用地的L和AL模型的ROC由0.708增加到0.769。这说明在分析土地利用变化的空间驱动时不可忽略土地利用类型空间自相关关系的影响。普通Logistics回归建模(L)中同时加入土壤质量因子和空间自相关性因子时(ALS)，模型的ROC增加最多，具有最高的解释精度。如，耕地和建设用地的ROC分别增加0.201和0.133，明显大于传统Logistic回归分析的ROC。

结合表4和上述分析结果，本研究采用同时加入土壤质量因子和空间自相关因子的回归建模结果用于CLUE-S模拟。研究区内年均降水量(MAP)是影响耕地分布的主要驱动因子，β系数为-0.618。除此，耕地还受与铁路距离(Railway_d)、GDP和全钾含量(TK)因素影响。除水体外，土壤养分含量对于土地利用分布概率有不同程度的影响。土壤有机质含量(SOM)和TK对耕地转化有较为明显的促进作用，全氮含量(TN)对草地转化起促进作用。全磷含量(TP)对林地转化有一定的促进作用，对其他类型用地转化作用不明显。林地和草地的土壤全钾含量β系数均为负值，代表土壤全钾含量与林地、草地分布概率呈负相关。TK对建设用地扩展起促进作用但其作用程度小于耕地。建设用地分布概率受GDP、高程(DEM)、坡度(Slope)和人口数量(POP)影响较大，主要驱动因子为GDP，β系数为0.561，代表GDP对建设用地扩展起促进作用。

表4 研究区各土地利用类型传统Logistic回归和Autologistic回归的β系数Tab.4 β value of Logistic regression and Autologistic regression for different land use types in Huainan

3.3.3CLUE-S模型检验与预测

以2010、2016年土地利用数据作为模拟的首期、末期，通过研究区2010—2016年间土地利用转移矩阵来计算转移概率，得到面积转移矩阵。测算2016—2040年研究区土地利用需求面积(表5)。

表5 研究区2016—2040年土地利用需求面积Tab.5 Land use demand in Huainan from 2016 to 2040 km2

图6为研究区2016年的模拟土地利用图与土地利用现状图，总栅格数为114 663个，模拟正确栅格数共90 832个。模拟结果的Kappa系数为0.74，模拟效果理想。

图7 研究区2028、2034、2040年土地利用预测结果Fig.7 Prediction result of land use in Huainan in 2028, 2034 and 2040

以研究区2016年为起点，预测2028、2034、2040年研究区各个土地利用类型的空间分布。结果(图7)表明：2016—2040年期间地类变化趋势，主要表现为耕地、林地和草地面积减少，水体和建设用地面积增加。2016—2028年间，耕地资源大量流失，预计由1 653.46 km2减少至1 426.61 km2，水体和建设用地增量均超过100 km2；后期2028—2040 年间，预计各地类面积变化明显减缓，建设用地面积变化不足1 km2，集中分布在淮河一带，耕地资源仍大量流失。水体呈现“填充式”和向外蔓延式增加，分布在泥河和高塘湖以及谢家集区的鲁村湾塌陷区和钱家湖塌陷区。未来实施耕地保护的同时也要加强塌陷水体治理措施。

永久基本农田保护、生态红线以及城镇开发边界、湿地保护等国家刚性政策规定的区域，在土地管理中有着极其严格的用途管制，限制转换为其他土地利用类型。本研究在土地利用变化模拟中，对于这类限制转换区域，并没有实际划分出来，默认所有区域均可发生地类转换。因此在后期的研究中，需要充分考虑该类因素，使得模拟结果更具有可靠性。

4 结论

(1)研究区土地利用类型以耕地为主，1985—2016年，耕地、林地和草地呈减少趋势，其中耕地面积持续减少且流失速率加剧；建设用地和水体面积呈增加趋势。2005—2016年，各种土地利用类型变化速率最快，建设用地最剧烈。土地利用转移主要发生在耕地、水体和建设用地之间，其中，1995—2016年，水体转入面积的93.44%来自于耕地。建设用地扩张也是通过占用耕地实现。

(2)空间自相关因子和土壤质量因子对土地利用类型转移的驱动作用不可忽略，且这两种因子对各土地利用类型的Logistic回归模型具有不同程度的贡献作用。其中土壤有机质含量、土壤全钾含量、土壤全氮含量、土壤全磷含量等土壤质量因子对耕地、林地、草地类型转移均有显著影响；空间自相关性因子对耕地和建设用地分布模拟有显著作用。

(3)研究区土地利用模拟的Kappa系数为0.74。模拟未来土地利用变化发现，耕地、林地和草地面积减少，水体和建设用地面积增加。2016—2028年，土地利用变化较快，到后期2028—2040年土地利用变化将减缓。