榆林地区植被时空分异特征及其影响因素研究

2022-03-11聂桐董国涛蒋晓辉郭欣伟党素珍郑嘉昊李立缠王江

生态环境学报 2022年1期

聂桐，董国涛, ，蒋晓辉，郭欣伟，党素珍，郑嘉昊，李立缠，王江

1. 西北大学，陕西西安 710127；2. 黄河水利委员会黄河水利科学研究院，河南郑州 450003；3. 黑河水资源与生态保护研究中心，甘肃兰州 730030

植被是生态系统联系大气圈、土壤圈和水圈的重要部分，也是直接反映生态环境变化的敏感因子（Wei et al.，2018），且植被作为陆地生态系统的重要组成部分，提供生态服务并维持陆地生态系统的功能（Du et al.，2020）。植被对自然环境和人类活动的变化具有敏感性，反映气候变化和人类活动在短时间内的影响（Meng et al.，2020）。归一化植被指数（NDVI）不仅是应用于全球和区域植被生长的重要指标（Yuan et al.，2019），也是用于监测植被时空变化的重要指标。目前对植被覆盖变化驱动力的研究集中于自然因子和人类活动因子的耦合作用，已有许多学者对植被动态变化的影响因素开展研究。

国内学者利用遥感数据监测区域内植被覆盖空间分布特征及动态变化（贾坤等，2013），对黄土高原（孙锐等，2020）、华北地区（刘斌等，2015）、三江源（刘宪锋等，2013）等地区的植被覆盖变化进行探索，研究发现气候变化是影响植被覆盖时空分异性的主要影响因素。全国尺度上，植被 NDVI空间异质性明显（高江波等，2019），水分以及温度在不同空间呈现一定的限制性。在小区域尺度上，地形要素（夏安全等，2020）、地貌类型（庞国伟等，2021）、土壤类型（吕家欣等，2020）的不同影响植被形成和生长状态，导致各地区植被分布存在明显的差异性。因此，各类自然因素对植被NDVI影响存在时空上的叠加效应，导致植被时空分布变化出现不稳定性和复杂性。随着城市化进程的加快，人类活动对植被覆盖变化产生显著影响，并且对植被覆盖的影响不是单一的，而是双重的（赵维清等，2019）。国内采取退耕还林还草等工程对植被恢复产生正影响，但过快的城市化发展对植被覆盖也产生负影响（于璐等，2020）。人类活动对国内植被覆盖影响具有明显的空间分异性，从整体来看人为活动负影响呈现波动减少趋势，正影响表现出增大趋势（王建邦等，2019）。例如秦岭山区和三江源地区，人为因素对植被覆盖呈现正影响关系（李辉霞等，2011；邓晨晖等，2018），表明生态保护等措施取得了初步成效。

目前研究植被NDVI时空变化及影响因素多数是线性回归分析、趋势分析与相关分析方法（肖强等，2016；郑杰等，2016；孔冬冬等，2017），但在长时间序列中植被变化与影响因素之间少有严格标准的线性关系（Hein et al.，2011），定量分析植被时空特征的驱动力仍比较困难。国内外许多学者对半干旱半湿润区域的植被覆盖时空变化进行研究（Gholamnia et al.，2019；Liu et al.，2019），陕西省榆林市位于半湿润区向半干旱区过渡地带（栾金凯等，2018），北部是毛乌素沙漠南缘，南部是黄土高原的腹地，沟壑纵横决定了榆林市生态环境脆弱，因此对榆林市植被时空分异性及其驱动力研究对生态保护具有现实意义。地理探测器模型既可以探测数值型数据，也可以探测定性数据（王劲峰等，2017），同时可以探测两因子交互作用于因变量。本文通过地理探测器模型定量分析NDVI空间分异性，将自然因素与人类活动因素进行整体研究，并找出影响因子对植被的最佳适应范围，为促进榆林市植被保护与生态格局建设提供科学依据。

1 数据和方法

1.1 研究区概况

榆林市位于陕西省最北部，处于 107°28′—111°15′E，36°57′—39°35′N 之间。东临黄河与山西相望，西连宁夏、甘肃，北邻内蒙古，南接本省延安市，地域总土地面积43578 km2（图1）。榆林市地貌北部为风沙草滩区，占总面积的42%，南部为黄土丘陵沟壑区，占总面积的58%。全市地势呈现西南高东南低，西南部平均海拔 1600—1800 m，其他各地平均海拔 1000—1200 m。榆林市是中国日照高值区之一，年平均日照时数2593.5—2914.4 h。2018年底，全市常住人口341.78万人（榆林市统计局，2019）。

图1 研究区地理位置Figure 1 The geographical position of the study area

1.2 数据来源与处理

本文选取植被NDVI、坡度、坡向、土壤类型、降水、气温、土地利用类型、人口密度和GDP等因子作为研究数据，用于探测榆林市的植被覆盖空间分异机制。NDVI采用的 MODIS数据（https://ladsweb.nascom.nasa.gov/）为时间分辨率16 d、空间分辨率为250 m的MOD13Q1 NDVI，并采用最大值合成。DEM原始数据空间分辨率为90 m，坡向、坡度由DEM数据计算获取。土壤和土地利用数据选自中国科学院资源环境科学数据中心（http://www.resdc.cn），空间分辨率为1 km。2000—2018年降水气温数据来源于中国气象数据网（http：//data.cma.cn），通过反距离加权平均法内插进行预处理生成2000、2005、2010、2015、2018年5期数据。2000—2018年人口密度GDP数据源自陕西省统计年鉴，通过空间插值对分县人口 GDP统计数据进行空间化处理生成5期数据（杨强等，2016）。坡度和土壤类型分别按照《土地利用现状调查技术规程》《1∶100万中华人民共和国土壤图》进行划分预处理。各因子数据经过投影转换、榆林市行政边界掩膜、重采样处理后，与250 m的NDVI数据像元大小保持一致。

1.3 研究方法

1.3.1 最大值合成法与等级划分

（1）最大值合成法（王正兴等，2003）（MVC）。选择具有 NDVI值的月合成数据取最大像元值生成，合成2000—2018年植被NDVI数据。（2）等间距法（彭文甫等，2016）。根据植被NDVI值大小，按照等间距方法进行等距划分，将植被 NDVI划分为低（0—0.2）、中低（0.2—0.4）、中（0.4—0.6）、中高（0.6—0.8）和高（>0.8）等5个植被覆盖度等级。

1.3.2 信息提取与自然因子分级

由于在实际操作时需要考虑精度与效率的平衡，最终选择在ArcGIS 10.4软件中借助渔网工具建立5 km×5 km格网，生成均匀分布的采样点共计1688个（图1）。将坡向、坡度、土壤类型、降水、气温、土地利用类型、人口密度和GDP数据采用自然间断点（刘宪锋等，2015）分级法分别划分为 9类、5类、9类、9类、9类、6类、9类和9类。将各自然、人文因子类型量数据以及榆林市NDVI值提取到点，再将各采样点因子数据和榆林市 NDVI数据导入地理探测器中进行运算。

1.3.3 地理探测器模型

地理探测器包括风险探测、因子探测、生态探测和交互探测4个组成部分，主要用来探测研究对象的空间分异性，同时揭示研究对象与驱动因子之间关系，是一种新的空间统计方法（王劲峰等，2017）。

（1）因子探测：主要用于探测NDVI的空间分异性，以及探测某自然人文因子（表 1）多大程度上解释了属性NDVI的空间分异，图2为各因子空间分布情况。具体计算：将NDVI图层与各因子分类图层叠加取值，进行探测分析，最终用q值表示（彭文甫等，2019），表达式为：

图2 研究区各因子的空间分布Figure 2 Spatial distribution of factors in the study area

表1 自然、人文因子指标Table 1 The indicators of natural and human factors

式中：

h=1, …；

L——变量植被 NDVI属性Y或自然、人文因子X的分层；

Nh和N——层h和全区的单元数；

和σ2——层h和全区的Y值的方差。

α和β——层内方差之和和全区总方差。q值越大说明某因子对榆林市植被NDVI的空间分异性解释力越强。

（2）交互作用探测：交互作用探测用于识别不同自然人文因子之间的交互作用，即评估两因子共同作用时是否会增加或减弱对因变量植被NDVI的解释力，或者两因子共同作用对NDVI的影响是否相互独立。首先，计算两类因子对植被NDVI的q值；其次，计算自然人文因子交互作用后的q值，对q(X1)、q(X2)与q(X1∩X2)进行比较（王劲峰等，2017）。

（3）风险区探测：风险区探测用于判断因子内部子区域之间的属性均值是否有显著的差别（Liu et al.，2019），用于判断适应植被生长自然、人文因子的最佳范围，用t统计量来检验：

式中：

Y——子区域植被NDVI；

——子区域h内植被NDVI属性均值；

nh——子区域h内样本数量；

φ——方差。

（4）生态探测：生态探测用于比较两类自然人文因子X1和X2对植被NDVI的空间分异性的影响是否有显著的差异，也可判断具体X1X2分两类因子共同作用对植被NDVI空间分布产生显著差异，以F统计量来衡量：

式中：

NX1及 NX2——两个自然因子的样本数量（Liu et al.，2019）；

αX1和βX2——由两个自然人文因子形成分层的层内方差之和（孔冬冬等，2017）；

L1和L2——变量X1和X2分层数目。

2 结果与分析

2.1 植被NDVI时空变化

2000—2018年榆林地区植被覆盖呈现明显增加趋势（图3a），增速为0.11/10 a。表2为2000—2018年植被NDVI动态变化数据，2000年中低等级植被覆盖区占榆林市面积的74.57%，到2018年中低等级植被覆盖区面积只占到10.59%。中等级以及中高等级植被覆盖区面积比 2000年增加了76.61%。2000—2010年间植被改善明显，植被覆盖面积呈现增加趋势，2010—2015年植被改善出现停滞，中低等级植被覆盖面积增加，气候变化引起的降水不稳定性也会导致区域NDVI变化具有较大的不确定性。2015—2018年中高等级植被覆盖面积达到近20年最大值。总之，2000—2018年低以及中低等级植被覆盖面积明显减少，而中等级以上植被覆盖面积明显增加，中高等级植被覆盖面积增速最快。为了进一步判断植被覆盖趋势变化的突变点，通过M-K检验（图3b）（Yue et al.，2002）发现，2000—2018年UF的值大于0，则表明序列呈上升趋势，其中2008年可能是突变点，通过滑动t检验（图3c）验证2008年确实是突变点，即2008年以后植被增长较为明显。

图3 （a）2000—2018年植被覆盖时间变化趋势（b）NDVI时间序列M-K检验（c）NDVI时间序列的滑动t检验Figure 3 (a) Temporal variation trend of vegetation cover from 2000 to 2018; (b) M-K test of NDVI time series;(c) Sliding T-test of NDVI time series

表2 2000—2018年植被覆盖面积Table 2 Area of vegetation cover from 2000 to 2018

2000—2018年榆林市植被 NDVI从中低等级向中以及中高等级植被覆盖区面积转变明显，说明榆林市植被恢复取得了重大进展。选择4期数据（图4）对比发现，2000—2018年植被覆盖空间分布变化明显。2000年榆林市植被NDVI中低等级植被覆盖区集中于榆林市西部，中等级集中于榆林市东北以及东南部；2018年榆林市植被NDVI中高等级集中于榆林市东部黄土丘陵区，中低等级集中于榆林市西北部的风沙区。2000—2010年，植被覆盖等级增加明显，东部地区由中低等级向中等级和中高等级转变。2010—2015年，植被覆盖等级由中高和中等级向中低等级转变，可能的原因与大范围干旱事件有关（刘彦随等，2017）。

图4 榆林市植被NDVI 2000—2018年时空变化Figure 4 Spatio-temporal changes of vegetation NDVI from 2000 to 2018 in Yulin

2.2 探测因子影响力

因子探测器用于表征各自然、人文因子对植被NDVI的影响力，通过计算得出的q值大小表示对NDVI的解释力。由表3可知，2018年各自然、人文因子的q值大小排序为：人口密度 (0.224)>气温(0.208)>GDP (0.183)>土地利用类型 (0.162)>坡度(0.141)>降水 (0.033)>土壤类型 (0.015)>坡向(0.008)。根据彭文甫（2019）对主要因子的选取，在8类自然、人文因子中，本文认为人口密度和气温两类因子对植被NDVI的解释力最大，解释力在20%以上，因此是影响榆林市植被覆盖空间分异性的主要因素。GDP、土地利用类型、坡度q值分别达到0.183、0.162和0.141，解释力超过10%，为次级影响因素。降水、土壤类型和坡向对NDVI的解释力未超过5%，表明这3类因子对榆林市植被NDVI分布存在间接影响。

表3 自然因子的q值Table 3 q value of natural factors

2.3 因子时间变化

2000—2018年5期数据结果显示（图5），土壤类型、降水对植被的影响力减小，坡向、坡度、气温、土地利用类型、人口密度和GDP的q值总体呈现增加趋势。其中，2000—2005年坡向、坡度、降水、气温、土地利用类型GDP的q值增加外，其余因子q值呈减少趋势；2005—2010年，除土壤类型q值减少外，其余因子q值都以增加趋势为主；2010—2015年除土壤类型q值减少外，所有因子q值呈现增加趋势；2015—2018年土壤类型、降水q值减少，其余因子q值呈现增加趋势。

图5 2000—2018年自然人文因子q值变化Figure 5 Changes in natural and human factors q value from 2000 to 2018

2.4 探测因子生态作用分析

生态因子探测用于比较两因子对植被NDVI空间分布的影响是否有显著性差异。本文中，行因子与列因子有显著性差异标记为“Y”，无显著性差异标记为“N”。由表4可知，坡向与坡度、气温、土地利用类型、人口密度和GDP对榆林市植被NDVI空间分布的影响存在显著性差异；降水与气温、土地利用类型、人口密度和 GDP之间存在显著性差异；土壤类型与坡向、降水之间不存在显著性差异。表明坡向、土壤类型和降水因子对植被空间分布没有显著影响，与其他因子共同作用对植被分布的影响存在显著性差异。进一步说明坡向、土壤类型、降水对植被影响力小，气温、人口密度和GDP对植被变化影响力大。

表4 生态探测结果Table 4 Ecological detection results

2.5 探测因子交互作用分析

多数因子交互作用的q值大于单一因子的q值，因子交互作用类型为非线性增强和相互增强关系，不存在独立关系（表5）。如X5∩X7(0.361)>X2∩X5(0.324)>X5∩X6(0.307)>X5∩X8(0.285)>X5∩X4(0.254)>X3∩X5(0.239)，结果表明，气温与坡向、坡度、降水之间交互作用呈现非线性增强效应，与土地利用类型、人口密度和GDP呈现相互增强关系；人文因子与其他因子之间对植被 NDVI的交互作用，如X5∩X7(0.361)>X2∩X7(0.352)>X6∩X7(0.321)>X4∩X7(0.322)>X1∩X7(0.257)>X3∩X7(0.244)，人口密度与坡向、土壤类型、降水之间交互作用呈现非线性增强效应，与坡度、气温、土地利用类型和GDP呈现相互增强效应。总之，8类自然、人文因子的交互作用比单因子对植被NDVI的影响更加显著，多因子交互作用不是独立的，而是相互增强和非线性增强显著关系。

表5 交互探测因子结果Table 5 Interaction detection factors results

3 讨论

3.1 人口密度

人类活动对生态环境的影响不可忽视，本文将人口密度因子划分为9个分区（图2），分别用A1—A9表示（表 6）。榆林市植被 NDVI均值在 A2、A6、A7、A8、A9分区中超过0.6，说明在这4个分区范围植被生长较好。统计检验表明，A9与 A1、A2、A3、A4、A5、A6、A7区植被 NDVI之间存在显著性差异，在 100—104 person·km-2范围时，植被NDVI值较大。A9区主要位于榆林市东南部，以绥德县、吴宝贤、佳县南部为主，综合水热条件分析，该地区降水气温适宜，为植被生长提供有利的生存条件。人口密度过大对植被生长具有限制性，人口增长使得城市建设占用大量绿地，在坡地大面积开荒耕种，造成水土流失，破坏并减少了植被覆盖，因此在城市建设中，需要增加绿地建设。人口密度过小区域位于榆林市定边县（图2g），属于资源性严重缺水县，年平均降雨量325 mm，定边县的自然环境条件约束此地区植被生长。同时，随着国民生态意识的提高，人类社会经济活动加大对生态修复方面的投入。

表6 人口密度每2个分区的植被NDVI均值及其统计显著性（置信水平95%）Table 6 Vegetation NDVI mean and statistical significance for every 2 divisions of population density(95% confidence level)

3.2 气温

气候因子对植被的生长具有限制性。本文将气温划分为9类，分别用B1—B9表示。表7表明，B8区植被NDVI均值最大，表明这一范围促进植被生长。统计检验表明，B8区与 B1、B2、B3、B4、B5、B6、B9区植被NDVI均值存在显著差异，与B7区无显著差异。因此在10.67—10.91 ℃范围，适合植被的生长。气温对植被的生长具有明显的适宜范围，适度的增温对植被活动产生积极作用。半干旱半湿润地区气温对植被生长影响大于降水，气温升高会促进光合作用，提高净生产力；然而，温度过高导致蒸发量增加以及土壤湿度下降（吕渡等，2018）影响植被变化。降水与气温的交互作用显著增强了气温对植被的NDVI的影响，水热条件协同作用有助于植被生长。

表7 气温每2个分区的植被NDVI均值及其统计显著性（置信水平95%）Table 7 Vegetation NDVI mean and statistical significance for every 2 divisions of temperature (95% confidence level)

3.3 其他因子

由于篇幅的原因，只选择q值大于20%的主要因子做详细讨论。根据风险因子探测结果总结出适合植被生长的最佳类型或范围（表8）。在95%置信水平下通过统计检验表明，植被 NDVI均值越大，处于这一均值的范围更适合植被生长。不同自然、人文因子或因子内部分区的植被NDVI均值存在差异。

表8 自然、人文因子适应范围Table 8 Adaptation range of natural and human factors

随着各自然、人文因子的变化，植被NDVI均值呈现波动变化。水热条件是决定植被生产力、空间布局以及动态变化的主要因素。降水和气温分别在 520.48—541.71 mm 和 10.67—10.91 ℃植被NDVI均值达到最大，分别为0.593和0.650，在这一范围内，暖湿的水热条件更利于植被活动。当降水因子影响力较低时（图 5），人类活动因子影响力较大，这时NDVI主要受人类活动影响。可能的原因是，降水不足以满足区域内植被生长水分需求时，人类活动一定程度上减缓水分限制，通过增加水分供给和转变土地利用方式。

坡向和坡度分别在 337.5°—22.5°和 6°—25°植被NDVI均值最大，分别为0.595、0.613，地形因子在这一范围内适合于植被生长。研究发现，在陕西省地势较低区域（<6°），是人类活动主要区域，植被覆盖相对较低；随着坡度的不断增加（6°—15°），人类活动减少，植被覆盖度随之增加；在15°—25°之间，通过退耕还林还草等生态工程的实施，耕地向林草地转变逐渐占据了该地形的主导（赵亮等，2019），人工种树和林草自然恢复使该区域植被覆盖度有所增加；在陕北>25°地区，与关中和陕南陡坡植被生长不同，该区域对陕北植被生长有一定限制性，较少有天然林生长和人工林种植。随着坡向的变化，植被NDVI均值变化差别不大，南坡接收太阳辐射多与北坡，北坡相对湿润且蒸发较少适合植被生长。

土壤类型划分为9类（图2），统计检验表明，1区植被 NDVI均值最高达到 0.6，钙层土（栗钙土、灰褐土、淡灰褐土、粘化黑垆土）较适合植被生长。相关研究表示，在降水较少地区，土壤水分是限制植被生长的重要因子。降水与土壤的交互作用增强了土壤类型对植被NDVI空间分布的影响，降水量的输入降低土壤水分的变异性（陈晓莹等，2020；彭文甫等，2019）。地貌类型影响土壤水分的分布（Hawley et al.，1983），榆林市土壤类型受地貌类型影响，西部多为风沙区以风沙土类型为主，表层土壤含水量极低，不利于植被生长。土地利用类型中耕地和林地土地类型适宜植被生长，NDVI均值达到最大0.614，林地土地类型为植被生长提供良好的地形条件。实施退耕还林还草工程以来，耕地转林地区域植被覆盖度正在持续增加，土地利用类型的转变对植被分布有直接影响。

GDP反映一个地区经济发展情况，经济发展过快对自然环境发展有所影响。将GDP划分为9分区，在1区植被NDVI均值最大，表明这一范围植被生长较好。统计检验表明，在266×104—420×104yuan·km-2范围内，植被覆盖较好。经济发展与生态环境相互影响、相互作用，两者互相协调影响区域长期可持续发展。社会经济发展从生态环境中获取自然资源，导致植被破坏、生态质量下降等不利影响，但是植被恢复工程也需要强大经济基础的支撑。根据已有研究（唐晓灵等，2020），榆林市生态环境治理投入的力度跟不上经济发展的速度，导致生态环境水平低于经济发展水平，处于不协调状态，生态环境问题影响到经济的可持续发展。为了经济与环境之间的平衡发展，采取退耕还林工程、封山育林、水土保持等措施，提高植被覆盖度，促进生态环境恢复。

人口密度以及经济发展对植被生长在一定范围内起到限制作用，退耕还林还草或者毁林开荒等人类活动都对榆林市植被带来或正或负影响。已有研究表明，黄土丘陵沟壑区退耕还林还草工程取得巨大成就，人类活动对陕西省植被恢复影响较大（黄麟等，2020；金凯等，2020），这与本文结论基本一致。建议协调好人类活动与自然环境的关系，对实现植被恢复和生态保护具有重要意义。

榆林地区处于中国退耕还林生态工程重点区域，自1999年实施退耕还林还草工程以来，植被得到明显改善，人类活动对生态环境产生不可忽视的积极影响。本文对2000—2018年NDVI的时空分异性进行分析，探究植被分布与自然因素和人类活动因素的影响机制。植被空间分布受到自然因素和人类活动的共同影响，与自然因素相比，近些年人类活动对植被恢复的影响更加显著。本文对年际NDVI进行探索，可以考虑生长季NDVI变化研究，此外，通过最大值合成法得到的NDVI数据较难与实测数据完全相同（彭文甫等，2016），深入研究可以考虑采用高空间分辨率遥感数据和详尽的实地调查，提高制图和分类的准确性。

以往对于榆林地区植被的研究大多采用线性回归和相关分析法（王涛等，2017；栾金凯等，2018），本文采用的地理探测器模型，能够有效的探测空间分异性，分析各种现象的驱动力以及多因子交互作用，在分析驱动力方面得到广泛运用（王劲峰等，2017）。本文运用地理探测器探测榆林地区植被空间分异性驱动力，得出的结果比较符合已有研究。值得注意的是，降水作用在陕北地区较小，可能与选择的尺度问题有关，在以后的研究中，可以选择更大的地理单元进行研究。