漓江流域林地时空变化及驱动力分析

付理宁, 王浩宇*, 冯宝坤, 王成,, 朱岚巍

(1.桂林理工大学测绘地理信息学院, 桂林 541000; 2.云南师范大学地理学部, 昆明 650000; 3.中国科学院空天信息创新研究院数字地球重点实验室, 北京 100094)

土地利用/土地覆被变化(land use/cover change,LUCC)可以直接表现人类活动对自然生态系统的影响,也是区域可持续发展的关键,长期以来一直作为全球气候、环境变化重点关注研究的内容[1],对引起土地利用/土地覆盖变化的驱动力分析即为总结和归纳土地利用类型和结构发生改变的因素。林地作为土地利用类型之一,不仅具有重要的生态价值、经济价值和社会价值,还是林业可持续发展的基础[2]。因此,了解林地时空变化及掌握其驱动力因素,对社会发展过程中出现的环境问题及实现可持续发展具有重要意义[3]。

目前,国内外学者在土地利用类型的空间动态变化方面的研究颇多[4-8],如刘根林等[9]采用土地利用转移矩阵、土地利用动态度等指标研究瑞兴于地区土地利用景观格局时空演变特征,并从地形地貌、社会经济及政策法规方面做出驱动力分析;艾敏等[10]以哈尔滨市呼兰区为研究区,从自然环境、社会经济、人口3个方面进行定性分析,采用监督分类、土地利用变化模型等方法分析研究区内土地利用变化特征。丁圆圆等[11]通过辐射传输法反演地表温度,分析永曲河流域地表温度与土地覆盖类型、地形因子及归一化植被指数之间的关系,其中土地覆盖类型的不同对地表温度产生的差异较明显。在林地时空动态变化研究方面,如唐源远等[12]基于土地利用变更数据采用转移矩阵、相关性分析等方法研究长株潭生态绿心林地景观格局动态变化,并从气象、人口及湿度方面进行驱动力分析;崔绿园等[13]基于土地利用遥感监测数据采用土地利用动态度、主成分分析等方法对重点区域林地时空变化特征进行分析,并从社会经济方面揭示了林地变化的影响因素。土地利用变化在研究方法上常采用监督分类等提取土地利用信息,根据土地利用转移矩阵及分类后动态变化等指标分析地类变化情况[14-17];从驱动力分析方面看,对于土地利用或林地动态变化的研究,在驱动力分析方面考虑了地形、气象、社会经济等影响因子[14, 18-19]。

漓江流域是典型的喀斯特地貌区,保水蓄水能力差,生态环境脆弱,容易导致石漠化和土地生产力下降等生态问题。漓江流域景观可持续发展监测与评估作为表征桂林市可持续发展进程与潜力综合案例的一部分,而林地变化又是表征区域景观变化的重要特征之一,因此,有必要对漓江流域林地资源动态变化及驱动力进行研究。基于上述分析,利用2013年、2017年、2020年共3期Landsat-8 OLI影像,采用随机森林分类方法进行土地利用分类并提取林地,借助土地利用转移矩阵及单一土地利用动态度等指标分析漓江流域林地时空变化特征,结合地理探测器[20]从地形、气象、经济及人口等方面进行驱动力分析,为漓江流域可持续发展的林地资源管理和保护方面提供相关参考信息。

1 数据与研究区

1.1 研究区

漓江流域位于广西壮族自治区桂林市境内,属于珠江流域西江水系,地处109°45′E～110°40′E、24°18′N～25°41′N,全长214 km,总面积达6 050 km2,流经兴安、灵川、阳朔、平乐及桂林市区,水质清澈,具有典型的喀斯特地貌特征,地理位置如图1所示。流域内旅游资源丰富,享有山青、水秀、洞奇、石美“四胜”之誉。属亚热带湿润季风气候,年平均气温约19 ℃,年均降雨量大于1 000 mm,呈现自西北向东南递减的趋势。流域内植被覆盖率达62%,包括了中国中亚热带所有的植被类型,其中高海拔地带的森林植被保存较好。

图1 漓江流域地理位置Fig.1 Geographical location of Lijiang River Basin

1.2 数据源

采用的遥感影像为2013年、2017年及2020年共3期Landsat-8 OLI数据,以及桂林市社会经济统计年鉴(2013—2020年),其他数据包括数字高程模型(digital elevation modal,DEM,空间分辨率30 m)、漓江流域边界矢量数据、行政区划以及英国国家大气科学中心(National Centre for Atmospheric Science,NCAS)制作的CRU TS(Climatic Research Unit gridded Time Series)温度和降水量数据。因研究区是典型的喀斯特地区,多云多雾,在Landsat-8数据选择方面,尽可能选择同时相或时间接近的影像,利用ENVI 5.3软件对影像进行了辐射定标、大气校正、多景镶嵌以及范围裁剪等,对DEM及气象数据用ArcGIS 10.6软件进行处理,得到研究区的坡度、坡向及年均降水量、年均气温数据等数据集。

1.3 样本数据

对预处理后的Landsat-8 OLI影像选取训练样本,参考2017版《土地利用分类》(GB/T 21010—2017)选择波段654(RGB)组合,参考相近时期影像、Google Earth高分影像,并结合漓江流域土地利用的现状及研究需要,将土地利用类型分为建筑用地、耕地、未利用土地、水体、草地和林地等6种类型。

2 研究方法

采用随机森林分类方法对研究区的土地利用类型进行分类,对林地的空间变化信息采用土地利用转移矩阵及单一土地利用动态度进行空间分析,最后利用地理探测器软件提取2013年、2017年和2020年影响漓江流域林地空间分布的共同驱动因子,继而对林地进行驱动力分析。

2.1 影像分类与精度评价

随机森林算法采用Bagging策略,强调随机有放回的抽取训练样本,根据基尼指数(Gini)最小原则进行分裂生成多棵决策树[21]。决策树作为一个整体来运行,需要识别两个参数来进行准确的类预测,即预测变量的数量和分类树的数量,最终由众多决策树以投票的方式形成具有N棵决策树的随机森林。此分类方法预测精度较高,异常值和噪声容忍度较好,训练速度快,在防止过拟合方面有明显效果,具有更准确的分类结果、用户自定义参数更少等优点[22]。利用验证样本集产生的混淆矩阵开展精度评价,混淆矩阵中包括总体精度、Kappa系数、用户精度以及生产者精度。

2.2 林地空间分析

2.2.1 土地利用转移矩阵

土地利用/覆被变化可以通过计算土地利用转移矩阵进行定量描述,可用式(1)表达不同土地利用类型之间的相互转化情况,表达式为

(1)

式(1)中:S为n×n矩阵,n为土地类型数;Sij为面积,i和j分别为研究期初和研究期末土地类型。

2.2.2 单一土地利用动态度

该指标表示研究区某一土地利用类型的面积在单位时间内的变化情况,某种土地利用类型的变化速度可用式(2)计算。

(2)

式(2)中:K为研究时间段内某一土地利用类型变化的剧烈程度;Ua、Ub分别为研究初期与研究结束时某种土地利用类型的面积;T为研究时长,当T的时间段设置为年时,K为该研究区某种土地利用类型的年变化率[9]。

2.3 驱动力分析

(1)因子探测。因子探测是对某因子多大程度上解释某属性空间分异方面的探测。将林地图层与各种影响因子图层进行叠加,探测所选择影响因子(具体影响因子如表1所示)间的重要性,用q度量,表达式为

表1 2013—2020年漓江流域影响因子指标Table 1 Indicators of impact factors in Lijiang River Basin from 2013 to 2020

(3)

(2)交互探测。交互探测用于探测识别不同的影响因子A和B共同作用后是否会增强或减弱对因变量的解释力。

3 结果与分析

3.1 分类结果

影像分类结果如图2所示,对分类结果进行精度验证,3期影像分类总体精度均在90%以上,且Kappa系数达0.8以上,满足数据处理和分析要求。

3.2 林地的时空变化

3.2.1 林地的变化特征

基于2013—2020年漓江流域地区土地利用转移矩阵,可以清晰看出研究区各土地利用类型之间的转移情况。由表2可知,2013—2020年漓江流域林地转出量较大,且与其他5种土地利用类型进行了不同程度的转化。其中,林地向耕地的转出量最大,转出面积达706.07 km2;未利用土地和建筑用地次之,转出面积分别为83.21 km2和57.37 km2;向林地转入的土地利用类型中,耕地向林地转入的面积最大,达到110.05 km2,林地的转出量多于转入量,面积呈现减少的趋势。

表2 漓江流域2013—2020年土地利用转移矩阵Table 2 Land use transfer matrix of Lijiang River Basin from 2013 to 2020

3.2.2 林地变化的动态度分析

2013—2017年、2017—2020年及2013—2020年3个时间段内林地的变化动态度如表3所示,可以看出,2013—2020年漓江流域的林地面积不断减少,且2017—2020年林地动态度变化最大。

表3 漓江流域2013—2020年林地变化动态度分析Table 3 Dynamic analysis of forest land change in Lijiang River Basin from 2013 to 2020

3.3 驱动力分析

3.3.1 因子探测

利用地理探测器软件分析各种因子对林地的影响,对选中的影响因子分别用不同的编码来标记,结果如图3所示。可以看出,2013—2020年对林地面积影响度在0.5以上的共同影响因子有X3、X6、X7(分别对应:第一产业产值、高程、坡度),其中共同影响因子中2017年和2020年的X6(高程)的影响力最大,q为0.857,解释力达85%以上,因此高程是2013—2020年林地共同影响因子中主要的影响因子。

图3 2013年、2017年、2020年影响因子探测Fig.3 Detection of impact factors in 2013, 2017 and 2020

3.3.2 交互探测

利用地理探测器探测林地与各驱动因子之间的相互关系如表4所示。研究发现各种影响因子对林地的交互作用类型存在双线性增强和非线性增强,无减弱和相互独立的情况出现;共同驱动因子间两两交互作用大于单一因子的影响,同时也反映了林地动态变化是自然和经济等多因素共同作用的结果。

表4 2013年、2017年、2020年漓江流域各影响因子交互作用情况Table 4 Interaction of various influencing factors in Lijiang River Basin in 2013, 2017 and 2020

3.4 结果讨论

运用ArcGIS 10.6软件将2013—2020年漓江流域流经各县区的林地与高程、坡度及第一产业产值3种共同影响因子进行叠加分析,进一步分析比较不同年份林地与3种共同影响因子之间的变化情况,如图4所示。图4(a)表明,不同等级的高程对各区县的林地面积存在一定的差异,其中2013—2020年兴安县、阳朔县、灵川县的高程变化值均在 1 450 m 以上,对应的林地面积在840 km2以上,以上区县的林地面积约为高程变化值较小的五城区的4倍。总体上,漓江流域流经各区县的高程变化越大,对应的林地面积相应较大。

图4 2013年、2017年、2020年漓江流域林地与共同驱动因子关系图Fig.4 Relationship between forest land and common driving factors in Lijiang River Basin in 2013, 2017 and 2020

不同等级的坡度对林地有不同程度的影响,林地随坡度的空间和时间分布情况与随高程的时空分布情况相似。从图4(b)可以看出:2013—2020年兴安县、阳朔县、灵川县的坡度均在16.5°以上,对应的林地面积均大于840 km2,而坡度在10°～14°的五城区、临桂区及永福县,其林地面积处在560 km2以内。整体上来看,2013—2020年漓江流域流经各区县的坡度越高,对应的林地面积越大。

第一产业产值的变化对林地的空间分布有不同程度的影响。从图4(c)中可以看出:2013—2017年、2017—2020年、2013—2020年3个阶段漓江流域流经各区县的第一产业产值均在不断增加,对应时间段的林地面积在不断减少。整体来看,2013—2020年漓江流域各区县的第一产业产值与林地面积呈现负相关的关系。

4 结论

研究从土地利用转移矩阵及动态度分析了漓江流域土地利用变化特征,从中提取出林地并分析与其他土地利用类型的转化情况,采用地理探测器分析了林地变化的驱动因素,主要结论如下。

(1)耕地是漓江流域内林地的主要转化类型,2013—2020年由林地转为耕地的面积为706.07 km2,转为未利用土地和建筑用地次之,转出面积分别为83.21 km2和57.37 km2;耕地向林地转入的面积最大,达到110.05 km2,且林地的转出量多于转入量,林地面积呈现减少的趋势。

(2)2013—2017年、2017—2020年以及2013—2020年时间段,漓江流域的林地面积从2013年的 4 689.37 km2不断减少到2020年的4 000.95 km2,2013—2020年间林地动态度整体上呈增大的趋势,其中2017—2020年林地动态度变化最大。

(3)高程、坡度和第一产业产值是2013—2020年影响漓江流域林地时空变化的共同影响因子,整体来看,2013—2020年漓江流域流经各区县的高程变化越大,对应的林地面积变化相应较大;坡度与林地的影响关系同高程与林地的影响关系类似;2013—2020年漓江流域流经各区县的第一产业产值与林地面积的关系呈现负相关。

对漓江流域林地变化以及产生变化的驱动因子进行分析,可以为有关部门管理林地方面提供相关依据。但是本文未考虑漓江流域边界涉及的全部县区(部分县区涉及漓江流域面积较少),在后续开展林地动态变化及驱动力分析的研究中会进行进一步的优化。