喀斯特山区耕地时空演变与影响因子分析

张 扬，周忠发※，黄登红，朱 孟，吴 跃，孙建伟

（1.贵州师范大学地理与环境科学学院/喀斯特研究院，贵阳 550001；2.贵州省喀斯特山地生态环境国家重点实验室培育基地，贵阳 550001；3.国家喀斯特石漠化防治工程技术研究中心，贵阳 550001）

0 引 言

耕地是人类赖以生存和发展的物质基础[1]，具有农业生产、空间承载以及环境保护等多种功能[2]，关系国家安全和社会稳定，是保障粮食安全和国民经济高质量发展的宝贵资源[3]。伴随中国新型工业化、城镇化建设的深入推进，在经济建设活动驱动下大量优质耕地转为非农用地[4]。由城乡要素转移引发的耕地利用非粮化、粗放化、边际化等问题日益突出[5-7]，导致近年来耕地质量下降、空间破碎、生态问题等问题频发[8]，加之后备资源不断减少、激励约束机制尚不健全，耕地合理利用和占补平衡面临多重压力。为此，《中共中央国务院关于加强耕地保护和改进占补平衡的意见》和新修订的《中华人民共和国土地管理法》均明确提出推进耕地质量提升和保护、大力实施土地整治、严格控制耕地转为非耕地，并将切实保护耕地定为中国的基本国策；2020年中央一号文件亦明确提出粮食生产要稳面积，坚守耕地和永久基本农田保护红线。耕地数量及质量与自然环境及经济社会发展阶段密切相关，特别是在喀斯特山区，耕地数量少、质量差且分布破碎，往往是保障发展与保护安全的矛盾所在。因此，研究耕地的时空演变规律以及探测其影响因子，已成为喀斯特山区提升农业发展质量、促进耕地可持续利用和推进乡村产业振兴的重要途径[9]，对研究区域可持续发展和区域粮食安全具有重要指导意义[10]。

国内外关于耕地的研究成果丰富[11]，研究内容集中在时空演变与驱动机制[12-14]、耕地保护与管理方法[15-16]、功能评价与时空演变[17]、坡耕地时空动态变化[18-19]、耕地变化与粮食安全[20]等方面。其中，由于耕地与其他地类相互转换的频繁性和复杂性，耕地的时空演变一直是学术界关注的热点[21]，诸多学者对耕地的时空演变格局[22-23]、时空动态变化特征[24-26]、时空演变驱动力[27-28]及其主导成因[29-30]进行了积极的探索，对耕地时空演变规律和趋势的把握已成为重要研究内容。例如，张英男等[31]以县域为基本单元分析了黄淮海平原耕地功能演变的时空特征及其驱动机制，刘旭华等[27]基于国家尺度定量探索耕地利用变化的驱动力，Yao等[23]分析20世纪80年代以来全球耕地变化的基本特征及其空间格局，赵馨等[32]利用深度学习探索不同石漠化环境下耕地的发展，史小祺等[33]研究贵州省近40年来耕地功能转型特征及其演变差异。综上所述，现有成果大多从特征、机制、成因等方面对耕地时空演变及其影响因素进行了分析，较少考虑与社会经济发展的耦合和顾及影响因素效应水平的空间异质性，驱动机制的定量研究较为鲜见，且空间尺度大多集中于中东部平原地区，西部喀斯特山区耕地长时间尺度的演变分析及其影响因子探测则较少。

贵州省作为中国南方典型的喀斯特山区，喀斯特地貌分布广泛，人均耕地少、耕地破碎化、耕地质量低，成为社会经济高质量发展的重要限制因素。在经济转型发展的关键期，耕地保护与利用的矛盾必将愈加突显，探索喀斯特山区耕地时空演变规律，定量分析耕地分布特征和时空演变驱动因素对区域粮食安全、生态建设、乡村振兴、产业结构调整等都具有重要意义。鉴于此，本文基于2000—2018年贵州省耕地利用时空动态数据，结合GIS空间分析功能和探索性空间数据分析、地理探测器等数理模型，研究耕地分布与时空演变规律及其驱动因素，以期为喀斯特山区耕地保护、耕地资源可持续利用与区域发展战略制定提供科学依据和理论基础。

1 研究区域、方法与数据

1.1 研究区概况

贵州省地处云贵高原东部，地跨 103°36′～109°35′E、24°37′～29°13′N，南北长约 509 km，东西宽约 595 km，总面积为17.62万km2。北靠四川和重庆，南邻广西，西毗云南，东连湖南，地势西北高东南低，海拔为148～2 900 m，境内山脉众多、重峦叠峰、绵延纵横、山高谷深（图1）。全省喀斯特地貌出露面积为 10.91万 km2，占总面积的61.92%，部分地区石漠化严重，是中国南方喀斯特地貌发育最典型的地区之一。气候类型为亚热带湿润季风气候，多年平均气温为 15.5 ℃，年降水量在 800～1 600 mm。截至2018年底，全省森林覆盖率58.50%，是长江、珠江上游重要生态屏障；省域内92.5%的面积为山地和丘陵，耕地破碎化严重，耕地地块坡度大，耕地总量少质量低；全省辖9个地级行政区（6个地级市、3个自治州），88个县级行政区，总人口3 600万人，GDP为1.48万亿元。

图1 研究区地形分布示意图Fig.1 Terrain distribution of the study area

1.2 研究方法与数据来源

1.2.1 研究方法

1）核密度估计。核密度估计是一种用于估计概率密度函数的非参数方法，通过对各要素点建立平滑的圆形表面，基于要素点到参考位置的距离建立其峰值或核来创建平滑的连续表面[34]。耕地分布具有空间离散、点多面广特征，难以直观识别其时空的整体性和连续性，核密度分析以光滑的曲面渐进式传输中心强度，体现空间位置的差异性以及中心强度随距离衰减的特性[35]，能够有效测度山区耕地分布的空间差异性和连续性。与普通的行政区耕地面积密度不同，本文采用核密度识别耕地变化密度特征。其计算公式如下：

式中f(x,y)为位置(x,y)的耕地密度估计，空间上表现为单位面积上所承载的耕地面积， hm2/km2；n为观测数量；h为距离衰减阈值（即带宽）；k为空间权重函数；di为位置(x,y)距第i个观测位置的距离，km。

2）空间分布均衡性分析。Lorenz曲线既可用于研究单个地理单元的空间分布，也可对多个地理单元的空间分布进行对比分析[36]。采用Lorenz曲线分析耕地空间分布的均衡性，以贵州省县级行政区数量为横坐标，以贵州省耕地面积累计百分比为纵坐标，绘制贵州省耕地洛伦兹曲线。Lorenz曲线的弯曲程度代表耕地空间分布的均衡性大小，曲线与绝对平均线的靠近程度和耕地在区域内分布的均匀程度正相关。

3）探索性空间数据分析。探索性空间数据分析通过对事物或现象空间分布格局的描述与可视化[37]，空间自相关分析是数据分析在空间中的表现，反映某个属性或地理现象在某区域和邻近区的依赖程度。本文运用莫兰指数（Moran’sI）分析耕地面积的空间集聚特征，以此判断变量在空间中是否具有相关性及相关程度的大小。全局空间自相关描述研究区域各地理要素属性值之间的关联程度和空间特征，可衡量区域之间整体上的空间关联与空间差异程度[38]，全局莫兰指数是广泛应用的全局自相关统计量，计算式如下：

式中I是莫兰指数；n是研究区的数目；xi，xj是区域i，j的耕地面积，hm2，¯x是xi的平均值，hm2；wij是研究区域i，j的对称空间权重矩阵元素。Moran’sI反映的是研究区域之间属性值的相似程度，即研究区域之间的空间聚集分布现象明显与否。若Moran’sI=0，则不存在空间自相关。

局部莫兰指数可以有效揭示影像分割在邻域对象间某一属性的相似性，局部莫兰指数的数值大小与邻域对象属性相似度呈正相关，与对象之间的可分性呈负相关，具体计算式如下：

式中Ii是第i个区域局部莫兰指数，其他变量的含义与莫兰指数公式中含义相同。当局部莫兰指数为正时，表示该区域与周边地区的空间差异小；当局部莫兰指数为负时，表示该区域与周边地区的空间差异大[39]。

4）耕地变化动态度。耕地变化动态度表示耕地资源分布的稳定情况，是反映区域耕地变化剧烈程度的重要指标[40]。利用耕地变化动态度来表达一定时间段内研究区耕地资源数量与分布的变化速度，定量评价耕地的变化速度。耕地动态度（K）表达式为[19]

式中K为研究时段内耕地动态度，%/a；Ua、Ub为期初和期末耕地的数量；T为研究时段长（5 a）。由于本文T的时段设定为年，故K值就是喀斯特山区耕地的年变化率。

5）土地利用转移矩阵。土地利用状态转移矩阵描述一段时间内系统中状态与状态的转移，可以全面又具体地分析出变化的结构特征，反映了各土地利用方式转化的流向与数量，能够揭示不同土地利用方式之间面积流动的时空演变特征[41]。其数学模型为

式中Sij为研究期初和期末土地利用状态，n为土地利用类型数。

6）地理探测器。影响耕地面积变化的因子错综复杂，归纳起来有社会经济和自然因素两方面，但从耕地资源的丰缺与相对于人类对耕地的需求而言，人类的社会经济活动对耕地面积变化起决定性作用[42]。地理探测器是探讨空间分异并揭示其背后驱动力的一组统计学方法，研究复杂地理因素的驱动作用机理，其原理是基于因变量与对其产生决定性影响的自变量在空间分布上具有显著一致性的假设[43]，被广泛用于探测土地利用变化[44]等相关领域。本文利用地理探测器的因子探测法定量评估影响研究区耕地分布及耕地时空演变的因素和影响强度。地理探测器中的因子解释力用q值度量，表达式为

式中L为变量Y或影响因子X的分层，即分类或分区；Nh为层h的单元数；N为研究区整体的单元数；σh2层h的Y值的方差；σ2是和研究区整体的Y值的方差。SSW和SST分别为层内方差之和与全区总方差。q为影响因素对耕地分布及耕地时空演变的影响程度，值域为[0，1]，值越大，说明影响因素对耕地分布及耕地时空演变的解释力越强。

1.2.2 数据来源

贵州省2000、2005、2010、2015、2018年5期土地利用数据以及道路数据、气象数据、降雨数据、行政区划数据均来源于中国科学院资源环境数据云平台（http://www.resdc.cn），其中土地利用数据空间分辨率为30 m。社会经济统计数据来源于《贵州省统计年鉴》（2000—2019）。不透水面数据（GAIA）来源于清华大学地球系统科学系科研资源数据库（http://data.ess.tsinghua.edu.cn/），空间分辨率为30 m。研究区DEM数据来源于地理空间数据云（http://www.gscloud.cn/），空间分辨率为30 m。

2 耕地时空分布与时空演变格局特征

2.1 时空分布格局特征

2.1.1 耕地面积密度分布特征

利用式（1）求得贵州省2000、2005、2010、2015、2018年耕地密度分别为 27.20、27.28、27.07、26.89、26.52 hm2/km2。总体上看，贵州省耕地面积在空间上呈现出“西北多、东南少”的典型特征（图2）。从图2可以发现，贵州省耕地面积密度变化较小，从2000年到2018年的核密度分布格局大体一致。耕地面积密度呈团状分布，具体特征表现为：耕地总量较高的区域主要位于贵州省西北部毕节市，北部遵义市以及中部安顺市，东部黔东南州黄平县等地带，而省会贵阳市域及周边县（区）、遵义市大娄山区域及赤水市、黔东南州大部、铜仁市梵净山区域和南部喀斯特地貌发育的册亨县、望谟县、罗甸县耕地总量较少。

使用自然断点法对各年度耕地总量进行 9级分类，2000—2018年贵州省耕地总量的时空分布呈现如下特征：①耕地面积总量大体呈现下降趋势，2000、2005、2010、2015和2018年耕地总量分别为497、498、494、490、484万hm²，2005年为拐点，2005年后耕地数量递减；②各年度耕地总量总体稳定，分级变化均发生在零星几个格网，且级别的变化不大，其余绝大部分格网未发生分级变化（图3）。这与近20年来贵州省实施的退耕还林还草、防护林生态工程、石漠化综合治理、水源地涵养林等生态建设以及工业化、城镇化和产业结构调整导致的耕地减少趋势相吻合。

2.1.2 空间分布均衡性分析

本文采用Lorenz曲线分析耕地空间分布的均衡性，根据区位熵计算结果，升序排列贵州省各县（区）的区位熵、各县（区）耕地面积占总耕地面积的百分比，并计算累计百分比。以贵州省县级行政区数量为横坐标，以贵州省耕地面积累计百分比为纵坐标，绘制耕地资源空间分布洛伦兹曲线（图4）。

通过Lorenz曲线进行判断，Lorenz曲线距离绝对平均线较远，说明耕地空间分布的均衡性不强。为了直观地分析耕地空间分布的空间差异性，采用基尼系数进行定量描述，在洛伦兹曲线的基础上，按照几何图形分块近似逼近的方法，求得基尼系数为0.326 3，分布较为不平均。结合Lorenz曲线及基尼系数进行判断，贵州省2000—2018年间的耕地空间分布的均衡性不显著，耕地的空间分布处于一种各区域分配较有差距的状态。贵州省境内喀斯特地貌分布广，部分地区石漠化严重，加之山脉众多、山高谷深，各区域由于自然条件的影响，导致耕地资源的空间分配不均衡。

图2 2000—2018年贵州省耕地核密度分布Fig.2 Kernel density distribution of cultivated land in Guizhou Province from 2000 to 2018

图3 2000—2018年贵州省耕地面积空间分布Fig.3 Spatial distribution of cultivated land in Guizhou Province from 2000 to 2018

图4 贵州省耕地资源空间分布洛伦兹曲线Fig.4 Lorenz curve of matching cultivated and land resources in Guizhou Province

2.1.3 空间分布关联性分析

使用Geoda空间分布关联性分析功能，全局自相关分析结果显示，2018年贵州省耕地面积的Moran’sI为0.528，经显著性检验Z值为42.754 2，P=0.001，呈现空间正相关。说明2018年贵州省耕地面积的空间分布并不是完全随机的，而是呈现部分空间位置上相邻的区域单元具有相同的属性值。局部空间自相关分析结果显示（图5），除不显著区域外，耕地面积的高-高聚类数量最多，其次为低-低聚类，而低-高聚类、高-低聚类极少。高-高聚类主要分布在毕节市、遵义市、安顺市、黔西南州、黔东南州、铜仁市的部分区域，如黔西县、金沙县、平坝区、织金县、务川县、道真县，这些县（区）与邻域自然条件比较相似，耕地面积均较高；低-低聚类分布最主要的区域是位于黔西南布依族苗族自治州，其余的零星分布于遵义市、黔东南市、铜仁市等部分区域。这部分区域大部分属于典型的喀斯特地貌或其它生态敏感脆弱区，耕地数量少、密度低；低-高聚类、高-低聚类的格网数量极少，基本都零星分布在贵州北面靠近与重庆市、四川省接壤的区域。

2.2 时空演变格局特征

2.2.1 土地利用转移矩阵

土地利用状态转移矩阵描述一段时间内系统中状态与状态的转移，为探明整个研究期的土地利用类型的转移形式，本文以研究期初（2000年）和期末（2018年）贵州省土地利用空间信息为基础进行空间叠置，得到近20年来贵州省土地利用转移矩阵，分析耕地空间变化的规律（表1）。

从表1可以看出，贵州省土地利用结构以林地和耕地为主，分别约占总面积的 53%和 28%，草地、建设用地、水域和未利用地的面积较小，约占总面积的20%。2000—2018年贵州省土地利用结构改变较为明显。从各类土地利用类型转换来看，2000—2018年的土地利用面积建设用地变化最为剧烈，其次为耕地，耕地变化以地类转出为主，面积净减少12.19万hm2，主要流向为林地、草地和建设用地。总体而言，2000—2018年贵州省耕地面积变化的主要趋势是大规模减少，而耕地的流向主要是林地、草地及建设用地，这和近年来大举实施的“退耕还林还草”以及社会经济的快速发展密不可分。

图5 耕地资源LISA聚集图Fig.5 LISA aggregation map of cultivated land resources

表1 2000—2018年土地利用转移矩阵Table 1 Transfer matrix of land use evolution during 2000—2018 hm2

2.2.2 耕地变化动态度分析

耕地变化动态度反映区域之间耕地变化的差异情况，以县域作为分析单元，进行动态度变化的分析。由图6可知，2000年至2005年为初期，2005年至2010年为中期，由图6可知，2010年至2015年为后期，2015年至2018年作为末期（下文同）。从空间分布上来说，耕地变化动态度的空间差异性不强，但从不同时期来看，各时期差异性较为明显。耕地变化动态度从初期的局部变化，过渡到中期变化较为强烈，后期的动态变化最为强烈，表明进入频繁的变化期，末期的变化动态度稍有回落。

2.2.3 耕地面积变化类型及其幅度分析

将 2000—2018年土地利用数据进行空间叠置，获取研究区内耕地面积总量及分布，以格网作为统计单元，将增加、减少两类型进行变化幅度的划分，由此得到研究区连续4个时间段内的耕地面积变化类型及其幅度（图7）。

①从空间分布上来看，贵州省耕地增减都较为均衡，没有明显的聚类现象；②从时间变化来看，初期耕地增减较多，中期、后期耕地以减少为主，零散分布增加类型格网，末期增加类型格网数量回升，但增加程度不大。结合具体的数据来看，耕地面积的减少量稳定，且在末期达到了较高的值，耕地面积的增加则由高值骤降至低值，并在低值处趋于稳定。从中期开始，耕地面积的减少值大大超过了增加值（图8）。2005年后，一方面由于退耕还林、城镇化、长江/珠江流域防护林、石漠化综合治理、水源地涵养林、城市周边防护林工程建设等，另一方面由于社会经济的发展，农村劳动力转移就业快速增加，耕地的数量逐期递减，由初期的耕地增加多进入到中、后期的耕地减少多。但在末期，随着“脱贫攻坚战”、“乡村振兴战略”的提出，农村劳动力部分回流，同时期耕地面积正增长的格网数量增加。

图6 2000—2018年贵州省耕地变化动态度空间格局Fig.6 Spatial pattern of cultivated land dynamic evolution in Guizhou Province from 2000 to 2018

图7 2000—2018年耕地面积变化类型及其幅度空间分布图Fig.7 Spatial distribution diagram of variation types and amplitude of cultivated land area from 2000 to 2018

图8 2000—2018年耕地面积变化Fig.8 Variation of cultivated land area from 2000 to 2018

2.3 耕地分布及其时空演变影响因子探测

地理探测器是通过研究要素之间的空间分层异质性解释其背后驱动因子的方法[45]，综合考虑自然因素和社会经济因素对贵州省耕地空间异质性的影响，以格网为基本单元选取人口（a1）、GDP（a2）、路网密度（a3）、开发强度（a4）、降水（a5）、气温（a6）、平均坡度（a7）、地面起伏度（a8）、共计8个驱动因子。经地理探测器计算后，各因子对耕地空间分布及耕地时空演变的贡献率如表2所示。从表中可以看出，3期数据中贡献率排序差别不大，2000年、2018年各因子对耕地空间分布贡献率从大到小依次为a1、a2、a3、a7、a4、a6、a5、a8；2010年顺序为a1、a2、a3、a7、a4、a6、a8、a5。

表2 耕地空间分布与时空演变地理探测器变量和指标Table 2 Variables and indicators of spatial distribution and spatio-temporal evolution of cultivated land in Geo Detector

地理探测器分析结果显示研究区的耕地空间分布和时空演变受多种驱动因子的影响，对不同时期的影响因子进行探测，得出：研究区耕地空间上的差异受自然因素和社会经济因素共同作用的影响，但社会经济因素影响力大于自然因素。其中人口（a1）、GDP（a2）对对耕地空间分布起主导作用，人口、社会经济差异是贵州省耕地空间分异最主要原因；路网密度（a3）、平均坡度（a7）的q统计量超过了0.15，是影响其空间分异的重要因素；开发强度（a4）、气温（a6）在0.10～0.15之间，也是影响其空间分异较为重要的因素；降水（a5）、地面起伏度（a8）q值小于0.1，说明这些驱动因子一定程度上影响着贵州省耕地的空间分异，但是影响程度较小。

3 结 论

本文以贵州省2000—2018年耕地为研究对象，借助GIS分析功能、Geoda空间分布关联性分析功能以及地理探测器等技术方法，探明了喀斯特山区耕地空间分布特征及其时空演变规律，并探测其分布和时空演变的主要驱动因子。主要研究结论如下：

1）贵州省耕地面积在空间上呈现出“西北多、东南少”的典型特征，但整体上耕地资源的空间分配均衡性不显著；耕地面积密度呈团状分布，省会贵阳市域、生态脆弱敏感区及喀斯特地貌发育区域耕地总量较少；聚类特点与耕地面积密度特征高度吻合，高-高聚类区域自然条件比较相似，耕地面积均较高；低-低聚类分布区域大部分属于典型的喀斯特地貌或其他生态敏感脆弱区，耕地数量少、密度低。

2）贵州省耕地面积总量与耕地密度均在 2005出现拐点，2005年后耕地数量、耕地密度均递减；从各类土地利用类型转换来看，土地利用结构改变较为明显，耕地变化以地类转出为主，主要流向为林地、草地和建设用地；耕地动态变化频繁，从时间序列上来看，初期为局部变化，变化程度随时间推移越加强烈，末期的变化稍有回落。从格网尺度进行分析，贵州省耕地增减都较为均衡，增减格网数量初期较多，中期至末期以减量格网数量为主；各期耕地面积的减少量稳定，且在末期达到了较高的值，各期耕地面积的增加量则由高值骤降至低值，并在低值处趋于稳定。

3）2000、2010 、2018年各期数据中对耕地空间分布贡献率排序极为相似，影响研究区耕地空间分布的主要因子为人口、GDP，次要因子为路网密度、平均坡度，社会经济因素影响力大于自然因素。

考虑到耕地与农村外出务工人员数据、贫困人口数据等的相互作用及影响，在后续研究中将多维贫困数据、外出务工人员数据及迁徙数据等与本文中的耕地空间分布和时空演变特征相结合，挖掘更多关于喀斯特山区耕地地理空间数据与多维贫困数据的潜在影响机理，可为区域粮食安全、生态建设、乡村振兴、产业结构调整提供有效理论支撑。