喀斯特山区耕地时空演变特征及驱动力分析：以毕节市为例

田 康， 蒙友波， 陈 洋*， 张 慧， 廖艳梅， 钟九生

(1．贵州省自然资源勘测规划研究院， 贵州 贵阳 550001； 2.贵州师范大学 喀斯特研究院， 贵州 贵阳 550001； 3.国家喀斯特石漠化防治工程技术研究中心， 贵州 贵阳 550001)

0 引言

【研究意义】耕地作为人类生存发展的重要基础条件，是宝贵的自然资源，承载着粮食安全重任，关系着国家安全和社会稳定[1-2]。近年来，随着城市化进程的不断加快，大量农村人口涌入城市，农村耕种人口减少，大量耕地撂荒，同时大量耕地转变为非农用地，耕地非粮化问题日益突出[3]。除此之外，耕地质量下降，耕地面积减小，破碎化程度加剧[4]。尤其是喀斯特地区，由于其地表破碎，石漠化问题突出，不合理的人类活动，使得人地关系矛盾突出，耕地数量和质量难以得到保障，发展受到制约。探究喀斯特地区耕地的时空演变特征及驱动因素，对区域耕地资源的保护、农业发展质量的提升以及粮食安全均具有重要意义。【前人研究进展】目前，国内外关于耕地研究内容集中在耕地质量与评价体系构建[5-7]、耕地时空演变与驱动因素探究[8-12]、耕地保护和管理方法[13-16]、功能评价与演变[17-18]、坡耕地时空动态变化[19-20]、粮食安全[21]等方面。其中，耕地的变化及驱动力分析一直是诸多学者的研究重点，主要方向为基于特定方法、模型或政策对驱动因子的模拟分析[22]，或在时间尺度进行耕地的变化及驱动因素对比分析或基于国家、经济圈、流域、省域、市域和县域等不同空间尺度的探究[23-28]。如罗晓虹等[29]利用景观分析软件，从时间尺度上探究1990—2015年间三峡库区重庆段水田的时空分布及演变特征，认为在未来的产业结构调整和高标准农田建设中，应稳定水田的数量和质量；蒲红梅等[30]从省域尺度上探明贵州省冬季耕地利用的空间分布特征，分析冬季耕地利用率、利用方式及其影响因素，认为应根据实际情况合理推进土地利用方式，提高农民收入；舒瑞等[31]则采用空间自相关、二元Logistic和核密度分析等方法探讨了宁夏沿黄生态经济带耕地的时空演变及驱动力，认为该区域的耕地保护基本实现动态平衡，耕地变化受经济、人口等社会经济因子综合作用。【研究切入点】综上研究对耕地从不同的时空尺度和地域尺度进行了探究，但研究多集中于我国地势较为平坦的中部、东部以及北部地区，对于地形复杂，具有独特的地表地下二元结构的喀斯特地区研究较少，且对喀斯特地区耕地自然和社会经济双重驱动因素的探索分析较为缺乏。贵州省作为典型的喀斯特山区，土地资源缺乏，耕地破碎且耕地质量低，严重制约了当地社会经济发展。为此，2022年1月国务院印发“关于支持贵州在新时代西部大开发上闯新路的意见”〔国发(20222)号〕中也明确指出，要严格落实全省耕地保护任务与责任，强化耕地数量保护和质量提升，调整优化耕地布局。因此，开展贵州典型喀斯特山区的耕地时空分布及驱动因素分析十分必要。【拟解决的关键问题】基于2009年和2018年贵州省毕节市耕地数据，运用GIS空间分析功能，探明耕地时空分布和演变规律以及其驱动因素，以期为喀斯特山区耕地保护、耕地资源可持续利用、优化布局与区域发展战略制定提供科学依据，助推乡村振兴更好更快发展。

1 数据来源与方法

1.1 研究区概况

毕节市位于贵州省西北部，乌蒙山腹地，地处川、滇、黔三省结合部，东靠贵阳市、遵义市，南连安顺市、六盘水市，西邻云南省昭通市、曲靖市，北接四川省泸州市，总面积2.69万km2。毕节市是乌江、北盘江、赤水河发源地，长江和珠江的生态屏障，拥有10 km以上长度河流共计193条，水资源丰富。季风气候比较明显，降雨量较为充沛，立体气候突出。境内地形地貌复杂多样，地势西高东低，交通不便，经济发展较为落后，是全国唯一以“开发扶贫、生态建设”为主题的试验区。

1.2 数据来源与处理

研究区土地利用数据以2009年和2018年的Landsat TM影像为数据源，运用ENVI5.3进行预处理，采用监督分类与人机交互解译的方法进行解译，借助Google Earth对解译后的图斑进行检查修正，最终得到2009年和2018年的土地利用数据。DEM数据是由地理空间数据云平台下载(https://www.gscloud.cn/search)，空间分辨率30 m，高程、坡度数据均由DEM数据计算所得，土壤侵蚀和地貌数据来自中国科学院资源环境数据云平台(http://www.resdc.cn)。人口和GDP数据来自各县的统计年鉴。降水和气温数据来自国家气象信息中心(http://data.cma.cn/)。

1.3 研究方法

1.3.1 格网分析法 为充分探讨该地区耕地分布特征，同时减小因行政边界和县域面积不同影响分析带结果，利用ArcGIS中的渔网功能，将研究区分为5 km×5 km的网格，共计1 219个，其中完整网格有927个，不规则网格292个。将耕地与格网数据进行叠加分析，统计每个格网中的耕地面积，并以此作为空间自相关分析的数据源，使结果更加具有对比性。

1.3.2 核密度分析 核密度分析是通过样本数据计算和估计数据聚集情况，并通过既定的距离衰减函数度量研究要素密度的变化情况，以此探索空间区域中的热点分布和变化特征[32-33]，实现耕地分布整体性和连续性的空间表达。核密度计算公式：

(1)

式中，f(x)为耕地地块分布核密度估计值(即单位面积上的耕地面积)，k(x)称为核函数，h为带宽，x-xi表示估计点到样本点x处的距离。

1.3.3空间自相关分析 空间自相关性是反映一个区域单元上某种地理现象，或某一属性值与邻近区域单元上同一现象或属性值相关性大小的重要指标，是空间域中值聚集程度的一种度量。常用的空间自相关性检验方法是用莫兰指数(Moran’s I)量化这种聚集属性，有全局空间自相关和局部空间自相关之分。全局Moran’s I的计算公式[34-35]：

(2)

式中，n为研究区空间观测对象的数量；xi和xj别为空间位置上第i和第j个观测对象的值(x所有对象的平均观测值)；wij为空间权重矩阵，表示空间位置上第i个和第j个观测对象的邻接关系。当对象i和j相邻时，wij=1；当对象i和j不相邻时，wij=0。

Moran’s I指数的取值范围[-1，1]。若I>0则表示空间正相关，说明研究对象在空间上有聚集性；当I<0时表示空间负相关，说明研究对象在空间上离散；若I=1，则表示计算结果未通过显著性检验，研究对象随机分布。

局部空间自相关(LISA, local indicators of spatial association)关注的是特定单元与相邻单元之间的空间关联性。局部Moran's I的计算公式：

(3)

式中，各变量的意义与式(2)相同。

1.3.4 耕地变化动态度 耕地变化动态度大小可以反映耕地变化的剧烈程度，计算公式[20]：

(4)

式中，K为耕地动态度，Ua为期初耕地面积，Ub为期末耕地面积，T为研究时段长度。

1.3.5 地理探测 综合考虑自然因子和社会因子对耕地的响应特征，运用地理探测器分别选取地貌(a1)、土壤侵蚀(a2)、高程(a3)、坡度(a4)、GDP(a5)、人口(a6)、降水(a7)和气温(a8)8个指标探测耕地时空分布特征的驱动因子。

2 结果与分析

2.1 耕地分布格局特征

2.1.1 时空动态度 毕节市各县域耕地面积及变化动态度见表1。从耕地类型看，在毕节市8个县区中，旱地在耕地中占据主导地位，这与喀斯特地区特殊的地形地貌密切相关，该地区地形复杂破碎，加之喀斯特地区特殊的地表地下二元结构，造成地表保水能力差，不利于水田的开垦建设。从其变化面积看，全市耕地面积处于减少状态，其中水田的变化最为明显，约减少1/2。从水田和旱地的分布县域看，旱地在威宁县的分布面积最广，七星关区和大方县次之，赫章县、织金县和纳雍县分布适中，金沙县的分布最少；而水田的分布则与旱地的分布情况相反，其主要分布于金沙县、织金县和黔西县，大方县和纳雍县次之，赫章县相对较少，威宁县分布最少，仅3.67 km2，为金沙县水田面积的1.89%。水田和旱地在一定程度上与各县域国土面积的大小、地形和海拔有关。从耕地的变化情况看，9年间耕地的综合变化动态度为-2.40%，整体上呈现减少的态势。其中，水田的变化最为明显，其变化动态度为-6.21%，是旱地变化的3倍左右。从各县域看，水田在威宁县和纳雍县的变化最大，其变化动态度分别为-10.60%和-10.20%，在黔西县的变化动态度最小，为-2.52%，最大约是最小变化动态度的5倍；旱地变化在金沙县和七星关区最明显，在威宁县变化最小，其变化动态度分别为-3.74%、-3.45%和-1.31%，其变化相差约3倍。

表1 毕节市各县域耕地面积及变化动态度

如图1所示，从单个网格的耕地面积占比看，在2009年时，各格网单元内耕地面积占46.20%，相当于100 hm2面积内有46.20 hm2耕地。其中，最大面积占86.95%，位于大方县；最小面积占1.25%，位于威宁县；除此之外，存在6个空网格，表明其中不存在耕地，存在与市的边界处。在空间分布上呈现“高低高低高”分布，分布较低的多位于各县交界处和地形复杂区，如威宁县的西南部、赫章县的东南部以及大方县和金沙县的交界处。在2018年时其分布特征与2009年时存在明显差异。其整体面积占比均呈下降态势，呈“高低高低”趋势分布。各网格单元耕地面积占36.09%。最大值为73.81%，位于纳雍县；最小值仅0.83%，位于七星关区；此外有8个网格为空值，存在与市的边界处。

图1 毕节市耕地占格网面积的比例

如图2所示，从耕地面积占比变化量看，全域的耕地面积占比大体表现为下降态势，最低下降62.81%，位于七星关区。部分网格单元呈上升态势，毕节市西部的赫章县和威宁县分布较集中，纳雍县、大方县、织金县和黔西县零星分布。从耕地综合动态度看，其与耕地的面积占比变化量空间分布大体相似，多表现为下降态势，毕节市的东部下降较明显，在毕节市西部存在部分上升，即存在部分单元格的耕地面积增加。

图2 毕节市耕地面积占比变化及综合变化动态度

综上来看，毕节市的耕地以旱地为主，主要分布于毕节市西部，由西到东大体呈下降态势，而水田则主要分布于毕节的东部，其分布区态势与旱地相反，由西到东呈增加态势，这可能受到海拔和地形的影响。县域单元内耕地的变化动态度均为下降态势，其中水田的变化最为明显，是旱地变化动态度的近3倍。在网格单元内，其综合变化动态势也呈下降趋势，其下降区多集中于毕节市的东部，而上升区多位与毕节市的西部，即威宁县和赫章县一带。

2.1.2 核密度 2009年和2018年毕节市耕地核密度如图3所示。毕节市耕地聚集性具有显著的差异性，在2009年时，毕节市耕地核密度为0.76～62.26个/km2，拥有3个极高密度区和8个高密度区；中高密度区、高密度区和极高密度区主要位于毕节市东部及东南部的金沙县、黔西县和织金县，其中黔西县分布最为密集，表明该地区耕地斑块较为破碎且集中；而低密度区、中低密度区和中密度区主要分布在毕节市的西部和中部，尤其是赫章县和纳雍县，表明这些地区耕地较为分散或斑块较为完整。在2018年时，其耕地核密度为2.42～62.70个/km2，数值上相对于2009年时变化较小，但从空间分布上看，变化较为明显。对比可知，在2018年时，毕节市耕地聚集性明显增加，虽极高密度区由原来的3个降至2个，但其局部的范围有所扩散，主要分布在黔西县；而高密度区由原来的8个增至19个，主要分布在除纳雍县和大方县以外的6个县区。低密度区和中低密度区向中密度区和中高密度区进行了大面积扩散，这主要是由于退耕还林还草、生态移民和易地扶贫搬迁等相关政策的实施，部分因坡度过大和符合退耕条件以及城市扩展等因素的耕地转换为其他用地，导致耕地不连片，破碎程度增加，密度增大。

图3 2009年和2018年毕节市耕地的核密度

2.1.3 空间自相关性 在2009年时，毕节市耕地面积的全局Moran’s I为0.452，Z值为21.35。在2018年时，耕地面积的全局Moran’s I为0.51，Z值为24.16。表明，毕节市耕地面积不是随机分布，而且具有明显的空间聚集性，且从全局Moran’s I的变化看，空间自相关性呈增大趋势，说明耕地聚集分布越来越强烈。由于其全局表现出明显的空间聚集性，因此研究进一步对其进行局部自相关分析，将其结果分为7类，分别为极冷点、冷点、较冷点、不显著、较热点、热点和极热点，其中极冷点、冷点和较冷点表示小面积耕地斑块的聚集区，其聚集程度逐渐减小；较热点、热点和极热点表示大面积耕地斑块的聚集区，其聚集程度逐渐增大；不显著表示中等面积耕地斑块且分布较为适中的区域(图4)。从局部自相关看，毕节市耕地在大范围内表现为不显著的相关特征，但同时也存在相关性明显的区域。在2009年时，毕节市大面积耕地斑块的聚集区主要分布在威宁县的中部，赫章县的西北部，七星关区、大方县和黔西县的西南部，小面积耕地斑块的聚集区则在各县区均有零星分布。在2018年时，毕节市耕地的局部自相关发生了明显变化，从大面积耕地斑块的聚集区分布看，威宁县的分布范围扩大近1倍，且部分向赫章县的东北部延伸，与赫章县内的热点分布区连接，而七星关区大面积耕地斑块的聚集区向不显著区转变。小面积耕地斑块的聚集区变化主要发生在赫章县的东南部，表现为小面积耕地斑块的聚集区向中等面积耕地斑块且分布较为适中的区域转变，其他区域变化较小。

图4 2009年和2018年毕节市耕地面积热点

2.2 耕地不同环境因素下的分布特征

2.2.1 不同地形 为去除县域面积不同对耕地分布造成的影响，同时探讨不同坡度梯度下的耕地分布特征，将坡度分为5个等级，并统计各个等级下耕地分布情况，以及不同等级下的耕地变化动态度。由图5可知，在两期耕地分布情况中，旱地的分布面积随着坡度的增大而增大，而水田则与之相反。一定程度上说明，旱地具有高坡度聚集性，水田则更倾向于低坡度区。究其原因：一方面在于水田对于水源的需求比旱地大，在坡度较低地区，水源更容易获取；另一方面地势平缓地区的保水能力更好。从耕地的变化动态度看，随着坡度的增加，其变化动态度也逐渐增大，旱地和水田的均在>25°内变化最强烈。这与退耕还林还草等政策有关。

图5 不同坡度梯度下的毕节市耕地分布情况

从耕地在不同高程海拔梯度下的分布情况(图6)看，旱地主要分布在低于2 400 m以下地区，大于2 400 m以上则分布较少，水田则全分布在2 400 m以下地区。可见，高程在一定程度上限制了水田的分布。从耕地的变化动态度看，旱地的动态变化度在低海拔地区较高，这主要是由于低海拔区易受到人类活动的影响，如城镇扩展、还林还草政策等多方面的因素。水田整体变化动态度大于旱地，在不同高程梯度上，随着海拔的增加，其变化动态也逐渐增大。综合看，海拔对水田的影响更强，受水田自身种植水稻的影响，而水稻对于气候的要求较高，适宜在气候条件相对稳定，温度适中的区域；当海拔过高时，气温相对于低海拔区要低，不利于水稻生长。

图6 不同高程梯度下的耕地分布情况

2.2.2 不同地貌 从不同地貌条件来看(图7)，2009年旱地主要分布区为侵蚀-剥蚀中山、溶蚀-侵蚀丘陵、溶蚀-侵蚀中山、溶蚀为主丘陵和溶蚀为主中山等地貌区，而水田主要分布在溶蚀-侵蚀中山、侵蚀-剥蚀中山和溶蚀为主丘陵，在溶蚀为主低山区则没有分布。2018年水田和旱地的主要分布地貌区，未发生明显变化，但在分布面积上均有所减小。从变化动态度看，旱地在侵蚀-剥蚀盆地、溶蚀-侵蚀中山、溶蚀为主低山和溶蚀为主盆地地貌区变化均较为剧烈。其中，侵蚀-剥蚀盆地主要为威宁草海位置所在地，为保护草海的生态系统，退耕周边部分耕地。溶蚀-侵蚀中山在毕节市的分布面积最广，基数相对较大，进而其变化动态度较强。溶蚀为主盆地所在位置多为各县区城区所在地，由于城镇化进程不断加快，部分耕地被占用或转变为其他用地，受人类活动干预最强，所以其变化动态度最强。水田在各地貌区的变化均较强烈，最低为溶蚀为主丘陵区，变化动态度为-4.49；最高为溶蚀-侵蚀盆地区，变化动态度为-8.89。

图7 不同地貌下的毕节市耕地分布

2.2.3 不同土壤侵蚀强度 从土壤侵蚀强度等级看(图8)，随着土壤侵蚀强度的加强，水田和旱地，均表现出递减趋势，土壤侵蚀强度越强，其分布面积越小。在2009年时，旱地在微度土壤侵蚀强度下的分布面积是极强度下的15倍，而水田在微度土壤侵蚀强度下的分布面积与极强度下的面积比也高达11倍。在2018年时，旱地在微度土壤侵蚀强度下的分布面积和极强度下的面积比有所增加，接近17倍，而水田在微度土壤侵蚀强度下的分布面积与极强度下的面积是2009年比例的3倍，高达30倍。从变化动态度看，水田和旱地的高变化动态度区均为极强度侵蚀区。可见，土壤侵蚀度也是影响耕地分布的重要因素。

图8 不同土壤侵蚀强度下的毕节市耕地分布

2.3 耕地时空演变及分布特征驱动因子

从表2看出，2009年各因子的贡献率(q值)为a5> a6> a3> a7> a8> a1> a2> a4，表明人口(a6)和GDP(a5)2个社会因子相比自然因子对耕地的驱动力较大，在自然因子中高程(a3)和降水(a7)的驱动力较为突出。2018年各因子的贡献率为a6> a1> a5> a3> a4> a7> a8> a2，自然因子整体依然小于社会因子的贡献率，但整体而言各指标的q值均较小，表明此时单个驱动因子对于该地区的耕地分布影响不大。从两期贡献率的变化情况看，自然因子中地貌(a1)贡献率未发生明显变化，表明基于此时的地貌分区对耕地的驱动影响相对稳定；其余指标均发生较明显变化，但社会因子依旧占据主导地位，可见人类活动对耕地时空演变及分布影响较大。综上，耕地的时空演变和分布特征受自然和社会多驱动因子的共同影响，其中社会因子对其影响相对较大，在该研究中GDP(a5)和人口(a6)对耕地的时空分布为主导因子；高程(a3)、坡度(a4)、降水(a7)和气温(a8)次之；土壤侵蚀(a2)影响最小；地貌(a1)相对稳定。

表2 耕地时空分布特征各因子的贡献值(q值)

从表3可知，两个因子间的交互作用相对于单个因子，其对耕地时空演变和分布特征的影响大。2009年a1∩a5、a1∩a6、a2∩a5、a3∩a5、a3∩a6、a3∩a7、a4∩a5、a5∩a6、a5∩a7、a5∩a8、a6∩a7和a6∩a8的贡献率q值均大于0.5，而a1∩a2的q值最小，仅0.142 9；a1∩a5的q值最大，为0.686 3。可见，在该时间点上，地貌、降水、气温等自然因子和GDP(a5)、人口(a6)社会因子的相互作用对耕地的影响较大。在2018年时，各因子间的交互作用贡献率q值均小于0.5，其中最大的为a5∩a6，其q值为0.423 8，最小的为a2∩a7，q值为0.066 6。此时社会因子交互作用的贡献率q值均大于自然因子交互作用的贡献率。综上，9年间耕地的分布主要是受到人类活动的影响，通过人为干预，使得耕地在一定程度上发生改变，这也与社会的发展趋势相符合。2009年社会经济不发达，地区主要以农耕等第一产业为主，其主要受到地形地貌以及气候等自然因素的制约，2018年已开始由第一产业向其他二、三产业转移，人类活动对土地的干预加强，如退耕还林还草、生态移民、城镇化建设和生态建设等，使得耕地发生转移，进而产生变化。

表3 耕地时空演变及分布特征驱动因子间的交互作用

3 讨论

前人对喀斯特山区的耕地研究基本位于贵州省经济较为发达地区，农业占经济发展的比重较小，无法为喀斯特贫困山区的农业发展提供较有效的技术支持，除此之外，其探究的影像因素相对单一[36]。本研究以位于乌蒙山区的毕节市为例，县域和格网2个单元对耕地的分布特征进行详细分析，在区域选择上具有较好的代表性，同时选取自然因子和社会因子，并运用地理探测手段进行影像因子分析，较全面地探究了耕地分布的影响因素。代仁丽等[37]通过对贵州省铜仁市的耕地分布特征及驱动因子进行探究发现，在耕地驱动因子中，坡度占主导作用。研究发现，毕节市耕地的分布主要受到社会因子的干扰较大，在耕地属性类别(旱地与水田)的分布时才主要受到坡度影响。从王雨枫等[38]对京津冀地区耕地利用时空分异及影响机制分析研究看，随着人民需求不断改变，社会经济及政策因素对耕地利用的影响力呈上升趋势，这与本研究结果相似，社会因子在一定程度上影响着耕地的分布，且影响程度大于自然因子。可见，耕地的面积和空间分布差异随着社会经济的发展产生改变，耕地的属性类别则主要受到自然因子的影响。本研究也存在不足，一方面在因素分析中仅选取自然和社会因子中的部分因子分析，因子的选取尚不够全面，无法全面揭示耕地时空演变的规律；另一方面，对于研究区域自然基地的影响分析尚不彻底，如土壤类型和岩性，在一定程度上可能也影响耕地的分布和变化，相对于平原地区，喀斯特山区居民点分布较为分散，但也会影响耕地的分布，除此之外路网密度等社会因素是否对耕地的分布和变化产生影响，下一步将加强该部分的研讨。

4 结论

根据2009年和2018年2期的耕地数据，结合相关的地形、地貌和土壤侵蚀数据分析表明，毕节市的耕地以旱地为主，其主要分布于毕节市西部，由西到东大体上呈下降态势，水田则主要分布于毕节的东部，其分布区态势与旱地相反，由西到东呈增加态势。研究区内耕地的变化动态度均为下降态势，其中水田的变化最明显，是旱地变化动态度的近3倍。毕节市耕地斑块的聚集性具有明显差异，2009年主要高核密度区集中在黔西县和金沙县，2018年各县区的中、高核密度聚集区均有所扩展。毕节市耕地的空间分布表现出空间自相关特征，其热点和极热点区主要分布在威宁县及大方县、黔西县和织金县的交界处。毕节市耕地的分布受到自然因子和社会因子的共同影响，其中社会因子对耕地的影响较为强烈。