邓元杰，侯孟阳，张 晓，贾 磊，李园园，姚顺波*，龚直文，刘广全

(1.西北农林科技大学经济管理学院，西北农林科技大学资源经济与环境管理研究中心，陕西 杨凌 712100；2.中国水利水电科学研究院，北京 100038)

森林在维持陆地生态系统平衡、保障人类基本生存条件、减缓气候变化和实现社会经济可持续发展等方面起着重要的基础性作用[1]。而林地作为森林资源的载体，不仅是林业可持续发展的基础，同时也是维护森林生态服务和环境价值的前提[2]，因此辨析林地的时空变化特征，及其驱动力因素，对了解区域生态环境变化、合理利用资源和协调人地关系具有重要科学意义。

已有学者对林地变化的研究主要是描述研究区林地时空变化特征，并从研究区的自然社会等宏观背景出发，定性分析可能引起林地变化的原因[2-3]。但这些研究无法定量揭示引起林地变化的驱动机制，导致决策者在制定林地保护政策时缺乏针对性与有效性，进而不利于林地资源保护及区域的可持续发展。有学者在林地时空变化特征描述的基础上，综合自然、人文、社会、经济等因素，构建合适的计量模型探究引起林地变化的原因[4-6]，然而线性回归模型容易忽略引起土地利用变化的驱动力(生物、自然和社会经济)或制度(政策)在空间上存有的高度异质性，不能处理林地变化过程中的空间因素变量，导致难以客观准确地揭示林地变化的原因。而Logistic回归模型不仅能考虑引起土地变化驱动力因素的空间异质性，处理土地变化过程中的空间因素变量[7]，还可以在区域土地利用/覆被动态变化的基础上定量分析其与社会、经济、技术、政策、自然环境等驱动因子的相互关系，进而综合不同时空尺度区域的土地利用/覆被变化过程和驱动力影响[8]。Shahbazian等[9]基于Logistic回归建模，探究了城市空间扩张的驱动力因素，并进一步模拟了未来城市空间扩张的方向，而李云龙等[10]则聚焦于耕地，运用Logistic回归模型剖析了影响耕地变化的驱动力因素。此外，还有学者从研究区土地利用变化的整体格局出发，对每一种土地利用类型变化均进行了Logistic回归建模，识别出影响各土地利用类型变化的不同驱动因素[11]。

然而，目前较少有学者将Logistic模型运用在探究林地变化的驱动因素上，虽然有学者运用Logistic回归模型对林地变化的驱动因素做了一定的解释[12-13]，但他们在构建模型时忽略了林地增减变化与其所对应驱动因素的空间异质性，即在构建林地变化Logistic回归模型时，将林地增加和减少均纳入Logistic回归模型中，而未将林地增加和减少进行单独建模。使得在探究林地变化的驱动因素时，无法有效将林地增加和减少的驱动因素单独分离，进而不利于对林地变化规律的准确把握。基于此，本研究聚焦于陕西秦巴山区林地变化，在分析2000—2018年林地时空变化的基础上，对林地增加和林地减少分别构建Logistic模型，探究影响陕西秦巴山区林地变化的驱动因素，以期为该区林地保护政策制定和林地资源可持续利用提供理论依据。

1 材料与方法

1.1 研究区概况

秦巴山区是秦岭地区和大巴山地区的简称，位于我国中西部地区(105°30′～110°0′E,31°40′～34°40′N)，地跨陕西、甘肃、四川、重庆、湖北、河南5省1市，其主体位于陕西南部地区，北倚秦岭，南靠巴山，汉江穿流而过，形成了两山夹一川的地貌特征，地形起伏大，山地面积广，盆地小且东西长，南北短，辖宝鸡、西安、汉中、安康和商洛5市38个县，该区土地面积约为86 995 km2，占陕西省总土地面积的42.31%。

1.2 数据来源与预处理

土地利用数据来源于地理国情监测云平台(http://www.dsac.cn/)，共包括2000年和2018年两期(图1)，其空间分辨率为30 m。该数据集的生产是以2000年Landsat 7 ETM+和2018年Landsat 8 OLS_TIRS遥感影像为数据源，通过人工目视解译生成，后期经过实地验证，误差修正后的土地利用以及类型综合评价精度达到93%以上[14-15]。同时，利用Google Earth高分辨率卫星历史影像与同期Landsat系列卫星影像对数据中存在的异常区域进行了再次验证与修正，确保了土地利用数据的客观性与准确性。土地利用数据分类系统根据“国家资源环境遥感宏观调查与动态研究”项目中的分类标准[16]，将陕西秦巴地区土地利用类型重新分为耕地、林地、草地、水域、建设用地和未利用地6类，其中林地由有林地、灌木林、疏林地和其他林地4种林地类型构成。

本研究基于ArcGIS平台，对两期土地利用数据进行空间分析，将2018年相对于2000年林地增加栅格，即2000—2018年其他地类转变为林地的栅格赋值为1，林地保持不变的栅格赋值为0，将其作为林地增加驱动力回归模型的因变量(图2A)；将2018年相对于2000年林地减少栅格，即2000—2018年林地转变为其他地类的栅格赋值为1，林地保持不变的栅格赋值为0，将其作为林地减少驱动力回归模型的因变量(图2B)。

参考类似研究成果[12-13]，并基于秦巴山区实际特征和数据代表性、可获取性，从自然、社会经济、地理区位3个方面构建了林地变化驱动力体系(表1)。各具体指标数据来源及预处理方法如下：①自然驱动力，主要包括海拔、坡度两个地形因子，年降雨量倾向率、年均气温倾向率两个气象因子以及土壤有机质含量。其中，地形因子数据来源于美国国家航空航天局在2020年最新发布的NASA DEM(https://lpdaac.usgs.gov);坡度数据为基于DEM数据，运用ArcGIS10.5软件中的坡度函数计算获得；气象因子的倾向率用以表征气象因子变化趋势，正值表示趋向增加或提高，负值表示趋向减少或降低，倾向率的计算方法详见文献[17]，气象因子数据来源于中国气象数据网(https://data.cma.cn/)提供的秦巴山区及周边区域12个气象站点2000—2018年逐年气候资料，计算得到降雨和气温的倾向率；土壤有机质含量提取自联合国粮农组织(FAO)所构建的世界土壤数据库HWSD V1.2(http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HWSD_Data/)。②社会经济驱动力，主要包括陕西秦巴山区2000—2018年各县总人口密度变化和年地均GDP变化。各县总人口和GDP数据均来源于陕西省统计年鉴(2001年，2019年)。③地理区位驱动力，主要包括到道路的距离、到铁路的距离、到城镇的距离、到农村居民点的距离和到河流的距离。其中，道路、铁路和河流的数据源于国家基础地理信息中心公布的1∶100万全国基础地理数据库(https://webmap.cn)，且在此基础上结合OpenStreetMap和最新的Google Earth影像进行修正。农村居民点和城镇数据则从秦巴山区2018年的土地利用数据中提取得到。为保证空间数据的一致性，以上所有数据空间分辨率均为30 m，空间坐标系均为WGS_1984_UTM_Zone_49N，所有数据均使用ArcGIS 10.5软件实现了空间化表达。

表1 2000—2018年林地变化驱动力指标

1.3 Logistic回归模型

二分类Logistic回归模型可以对二分类因变量和自变量(或连续变量或混合变量)进行回归建模，设第i个案例的自变量xi=(x1、x2……xn)，二分因变量Yi取值0或1(Yi=0表示主体事件不发生，Yi=1表示主体事件发生)。即林地增加时Y=1，否则为0；林地减少时Y=1，否则为0。该模型方法是基于抽样数据，为各自变量产生回归系数，并通过这些系数来讨论模型中因变量与自变量的关系。设P为事件发生概率，取值范围为0～1，则1-P为该事件不发生的概率，这种概率可以用Logistic函数计算，其表达式为[18]：

P=(Y=1|x1,x2……xn)=

(1)

通常对公式(1)进行Logit变换，以转变为线性公式：

(2)

式中:x1、x2……xn为自变量；参数β0、β1……βn为待求的回归系数；P/(1-P)为1组事件的发生比，发生比率expβ是衡量自变量对因变量的重要指标。expβ表示自变量每增加1个单位，事件发生比的变化倍数(expβ<1，发生比减少；expβ=1，发生比不变；expβ>1，发生比增加)[18]。Waldχ2用来评价每个自变量对因变量的解释程度[18]。对Logistic回归方程拟合度的检验选用Homsmer-Lemeshow 指标(HL)，当HL指标统计显著表示模型拟合不好；相反，当HL指标统计不显著表示模型拟合程度较好。除此之外，对Logistic回归模型的解释能力检验，通常采用受试者工作特征曲线(ROC曲线)进行验证。该曲线常被用于度量 Logistic 回归模型的吻合度, ROC曲线下面积的值为 0.5～1.0, 其值越高, 表明Logistic回归分析模型的吻合度越好, 反之越差[19]。

2 结果与分析

2.1 2000—2018年秦巴山区林地变化特征

基于秦巴山区2000年和2018年两期土地利用数据，借助ArcGIS 10.5软件对其林地变化进行了空间分析，再运用Excel软件对分析结果进行统计可得，2000—2018年秦巴山区林地面积由614.84万hm2增长至618.88万hm2，林地变化以转入增加为主，18年间共转入23.22万hm2，转出19.23万hm2，净增加39 951.72 hm2。由林地变化对地形因子的响应可得(图3)，林地变化主要集中分布在海拔≥500～1 000 m、≥1 000～1 600 m，坡度≥15°～25°、≥25°范围。其中，在海拔<500 m和坡度<6°范围内，林地面积减少大于增加；而在海拔≥500 m和坡度≥15°范围内，林地面积则为增加大于减少。

林地变化结构见图4，陕西秦巴山区耕地转为林地是林地面积增加最主要的原因，18年间共有21.83万hm2耕地转为林地，占林地转入总面积的94%。此外也有少部分的草地(9 824.94 hm2)、未利用地(1 880.73 hm2)、建设用地(1 405.71 hm2)和水域(845.91 hm2)转为林地。在林地转出方面，林地转为耕地是导致林地面积减少最主要的原因，18年间共有16.95万hm2林地转为耕地，占林地转出总面积的88.17%。此外分别有11 303.91、4 829.04、3 553.20 和3 049.83 hm2的林地转为草地、未利用地、水域和建设用地。从林地变化的空间分布来看，林地增加和减少在秦巴山区各处均有分布。但林地增加主要分布在秦巴山区南部和西南部(图2A)；而林地减少主要位于秦巴山区东北部(图2B)。

2.2 林地变化的Logistic回归分析

在构建林地变化驱动力回归模型之前，对各自变量进行共线性问题诊断，结果自变量的容忍度(tolerance)为0.50～0.96，方差膨胀因子(variance inflation factor, VIF)均小于 10，表明自变量间共线性不严重，均可以纳入Logistic回归模型。

2.2.1 林地增加的Logistic回归分析

在林地增加的Logistic回归模型中，坡度用4个虚拟变量分别代表坡度Ⅱ(≥2°～6°)、坡度Ⅲ(≥6°～15°)、坡度Ⅳ(≥15°～25°)和坡度Ⅴ(≥25°)，坡度Ⅰ(<2°)作为它们的参照对象。在HL检验中，Sig.值为0.355，大于 0.05，统计不显著，即模型的拟合效果很好; 并由回归结果可知，似然比卡方检验的观测值为122.51(P<0.001)，且ROC曲线下的面积(AUC)大于0.7，该模型通过了各项检验，可用于林地增加的驱动力分析。

显著性水平和Waldχ2统计量结果见表2，2000—2018年，秦巴山区林地增加较为重要的自变量依次为海拔、到城镇的距离、年均气温倾向率、坡度、到农村居民点的距离、土壤有机质含量、到道路的距离和到铁路的距离。其中，海拔、到城镇的距离、年均温倾向率、到道路的距离和到铁路的距离对林地增加呈显著负向影响。具体而言，海拔每升高1 m，林地面积增加的概率将降低0.19%；到城镇、道路和铁路的距离每增加1 m，林地面积增加的概率将分别降低0.002%、0.010%和0.001%；年均气温升高将不利于林地面积增加。坡度、到农村居民点的距离以及土壤有机质含量对林地增加呈显著正向影响。具体而言，坡度Ⅱ(≥2°～6°)、坡度Ⅲ(≥6°～15°)、坡度Ⅳ(≥15°～25°)和坡度Ⅴ(≥25°)这4个坡度范围对林地增加的影响均为显著，其产生的正向影响均大于坡度Ⅰ(<2°)。且根据表3中Waldχ2统计量可得，坡度Ⅴ(≥25°)的Waldχ2统计量最大，表明其对林地面积增加的正向影响相比于其他坡度范围更大。即以坡度Ⅰ(<2°)作为参照对象，林地面积增加概率在≥25°的坡度范围内将会升高37.5%。到农村居民点的距离每增加1m，林地面积增加的概率将会升高0.004%；土壤有机质含量每增加1%，林地面积增加的概率将会升高4.64%。

表2 2000—2018年秦巴山区林地增加驱动力因素分析

2.2.2 林地减少的Logistic回归分析

在林地减少的Logistic回归模型中，坡度用4个虚拟变量分别代表坡度Ⅱ、坡度Ⅲ、坡度Ⅳ、坡度Ⅴ，坡度Ⅰ作为它们的参照对象。HL检验中，sig.值为0.148，大于 0.05，统计不显著，即模型的拟合效果很好; 由回归结果可知，似然比卡方检验的观测值为103.45(P<0.001)，且ROC曲线下的面积(AUC)大于0.7，该模型通过了各项检验，可用于林地减少的驱动力分析。

根据显著性水平和Waldχ2统计量(表3)可得，2000—2018年，秦巴山区林地减少较为重要的自变量依次为坡度、海拔、到城镇的距离、到道路的距离、到农村居民点的距离和年均气温倾向率。其中，海拔、到城镇的距离、到道路的距离及年均温倾向率对林地面积减少呈显著负向影响。具体而言，海拔每升高1m，林地面积减少概率将会降低0.12%；到城镇、道路的距离每增加1m，林地面积减少的概率将会分别降低0.010%和0.002%；年均温升高将有利于抑制林地面积减少。坡度和到农村居民点的距离对林地面积减少呈显著正向影响。具体而言，坡度Ⅱ(≥2°～6°)、坡度Ⅲ(≥6°～15°)、坡度Ⅳ(≥15°～25°)和坡度(≥25°)这4个坡度范围对林地面积减少的影响均为显著的，其产生的正向影响均大于坡度Ⅰ(<2°)。且根据表3中Waldχ2统计量可得，坡度Ⅳ(≥15°～25°)的Waldχ2统计量最大，表明其对林地面积减少的正向影响也最大。即以坡度Ⅰ(<2°)作为参照对象，林地面积减少概率将会增加153.25%。到农村居民点的距离每增加1 m，林地面积减少的概率将会升高0.003%

表3 2000—2018年秦巴山区林地减少驱动力因素分析

3 讨 论

1)2000—2018年，陕西秦巴山区林地面积增减共存，但总体呈增加态势，且主要集中分布在中高海拔和坡度>15°范围内，这与张碧桃等[20]研究结果一致。究其原因，陕西秦巴山区受退耕还林、长江天然林保护等生态修复工程的影响，坡耕地成为退耕还林的主要对象，再加之中国经济的快速发展，山区农民进城务工趋势明显，耕地撂荒严重，农用地成为林地增长的主要来源[21]。与此同时，道路修建、居民点建设和非法采矿等人类社会经济活动的开展是导致林地减少的重要原因[22]，但从总体来看，林地增加面积高于林地减少面积。此外，林地面积在低海拔和坡度<6°范围内减少大于增加，这与城市扩张、交通设施建设等活动侵占了少量林地密切相关[23]。

2)Logistic回归模型结果显示，海拔和坡度是影响林地变化的主导驱动因素，这与前人的研究结论相同[12-13]，然而具体到秦巴山区，随着海拔的增加，林地增加和减少的概率均呈逐步降低态势，这一结果与本研究相吻合。林地增加概率随海拔升高而降低，这与生态修复工程的开展多集中于开发更容易、成本较低的中低海拔地区有关，因为该区域水热条件充沛，工程实施可使区域林地在短时间内得到增加，实现生态修复目标[13]。同时，社会经济活动开发难度及其成本随海拔的升高而增加，高海拔地区林地受人为活动扰动小，因此林地减少概率随海拔上升而下降。林地增加概率随坡度上升而增大，此结论与Xie等[24]的研究结论一致。中国自20世纪末在全国开展的以退耕还林为主的一系列生态修复工程是引起林地随坡度上升而增加的重要原因[25]；相反，林地随坡度增加而减少，这可能与道路修建、森林砍伐、居民点建设和非法采矿密切相关[24]。其他两个自然因子气温倾向率和土壤有机质对林地变化也产生了较显著的影响。其中，从年均气温倾向率对林地变化影响来看，2000—2018年，研究区内年均气温整体呈升高态势。林地增加和减少的概率均随着年均气温倾向率的增加而降低，虽然温度适度升高有助于树木生长[26]，并一定程度上阻缓林地退化，但秦巴山区的植被生长对气温变化敏感度较高，年均气温升高反而不利于植被的生长与林地面积的扩大[4,27]，因为温度升高所诱导的干旱可能会导致研究区域树木生长缓慢甚至死亡[28]。土壤有机质对林地增加的影响为正，对林地减少的影响却不显著，这是由于土壤有机质含量较高的林地土壤肥力综合指数较大，秦巴山区人工林工程有效地改善了土壤养分，维持了土壤的肥力，提高了土壤持水的能力[29]，从而使得林地增加的可能性随着土壤有机质含量的增加而升高。

3)地理区位因子也是影响秦巴山区林地变化的重要因素，到城镇的距离以及到道路的距离，均对林地变化产生负向影响，表明林地增加和减少的概率随着到道路和城镇的距离增加而下降，随着秦巴山区社会经济和城镇化的快速发展，一方面区域人口增长和经济发展所产生的大量建设用地扩张需求导致林地被占用[23]；另一方面，陕西秦巴山区地方政府在生态文明建设等一系列生态保护政策的引导下，越来越重视生态环境的保护和治理，追求环境与经济的协调发展[30]，通过在城市、道路周边植树造林，增加了城市周围林地覆盖率，改善了城市的生态环境。此外，模型结果显示，林地减少概率随着到农村居民点距离增加而增大，这可能与农户作为退耕还林等生态修复政策的主要参与者，其退耕地离农村居民点较近有关[31]；同时，离农村居民点越近的林地，受人类活动干扰也越大，如道路建设、毁林开荒等[24]，因此距离农村居民点越近，林地增加可能性越小。

4)通过Logistic回归分析，本研究虽较好地揭示了引起陕西秦巴山区林地转化的驱动机理，但林地转化的驱动机制极其复杂，与研究区的社会、经济、人口、自然因素、政府政策与价值观念等密切相关[8]，且受限于数据量化方法的可操作性，导致有些因素被遗漏或难以量化而未纳入模型。首先，诸如政府政策、自然保护区影响程度等因素在特定区域与时间内对林地转化驱动巨大，但难以精确合理地量化[10]，例如促使林地增加的退耕还林等生态修复政策；其次，受困于数据获取，导致本研究所使用的社会经济数据只能统计到区县一级，如能使用更精确、细化的社会经济数据，分析结果将更加精准。