2001-2020年天水市植被覆盖时空变化及其驱动机制研究

2022-11-29靳惠安唐志红姚晓军张大弘黄登冕张亚雄

西北林学院学报 2022年6期

靳惠安，唐志红，姚晓军，张大弘，黄登冕，张亚雄

(1.甘肃林业职业技术学院，甘肃天水 741020；2.西北师范大学地理与环境科学学院，甘肃兰州 730070；3.西北大学城市与环境学院，陕西西安 710127)

植被作为陆地生态系统的重要组成，对气候变化极为敏感[1]。在全球气候变暖和人类活动不断增强背景下，植被生长环境的变化对植被覆盖度产生重要影响[2]。植被覆盖度指植被(包括叶、茎、枝)在地面的垂直投影面积占统计区域总面积的百分比[3]。植被覆盖与区域生态系统的变化存在密切联系，如气候变化、能量和水分循环、微生物及人类活动等，是评价生态系统的基本指标，对探究地球系统不同圈层间相互作用和不同尺度的生态环境变化具有重要意义[4-7]。

卫星遥感影像因具有观测范围大、重访周期短、空间分辨率多样等优点，已成为全球及区域植被覆盖动态监测的主要手段[8]。基于多源遥感影像结合不同植被指数和模型来估算植被覆盖是当前植被遥感反演与评价的重要研究内容[9]。归一化植被指数(Normalized Difference Vegetation Index,NDVI)由于能消除部分太阳高度角、大气云层与地形等因素产生的辐射干扰[10]，已广泛应用于大范围地表植被动态监测[11-14]。植被的演变是自然条件与社会经济因素综合作用的结果[15]，目前，有关植被覆盖与自然因素相关性的研究，主要集中植被覆盖对降水和气温变化的响应。Liu等[16]对全球植被覆盖变化的研究认为植被改善与气候变暖密切相关，尤其在干旱地区，降水与植被变化的相关性更明显。李登科等[17]对退耕还林后陕西省植被覆盖度变化及其对气候响应的研究表明，退耕还林等生态工程的实施在一定程度上降低了植被对气候因子的敏感性。张珍珍等[10]对梵净山植被覆盖时空变化驱动因素的分析表明，人类活动是导致地表区域植被变化的主要因素。张学玲等[11]对武功山植被覆盖的空间分异特征的研究发现植被覆盖与地形因子存在较高的相关性。以上研究表明，受自然条件和区位的影响，不同区域植被覆盖变化的驱动因素存在差异。

近年随着西部大开发战略、丝绸之路经济带倡议和黄河流域高质量发展重要战略的实施，西部地区生态文明建设成为热点话题。天水市因地处黄土高原与秦岭山地交错的干旱半干旱过渡地带，是典型的农林生态复合区，且纵跨黄河、长江两大流域，其生态环境对该市及周边地区的经济、文化建设尤为重要。本文利用谷歌地球引擎(google earth engine,GEE)获取天水市2001-2020年Landsat系列遥感影像，构建NDVI时间序列数据，采用Theil-Sen Median趋势分析、Mann-Kendall检验、Hurst指数以及地理探测器方法，探究近20 a天水市植被覆盖动态变化与驱动因素，明晰植被覆盖对区域气候、自然条件及人类活动的响应特征，从而为区域生态环境保护与治理提供科学依据。

1 材料与方法

1.1 研究区概况

天水市(34°05′-35°10′N，104°35′-106°44′E)地处甘肃省东南部，位于我国黄河流域与长江流域的交接部位，也是丝绸之路经济带的重要节点城市(图1)。全市现辖两区(秦州、麦积)五县(武山、甘谷、秦安、清水与张家川回族自治县)，总面积约1.43×104km2，总人口370.40万。境内地势西北高、东南低，在地质构造上处于中祁连造山带和秦岭造山带的接合部位。渭河横贯中部，沿河接纳多条支流；南部的西汉水和永宁河是嘉陵江上游主要支流。天水市年均降水量500～700 mm，年均温8～12 ℃，属温带大陆性气候，是典型的农林牧交错复合区、气候变化敏感区和生态环境脆弱区[18]。因地处华北、华中、蒙新和喜马拉雅植物交汇处，境内森林资源丰富，主要分布在关山、小陇山和西秦岭林区。植被类型主要包括温带针叶、亚热带针叶、热带山地针叶、温带落叶阔叶、温带落叶灌丛，草甸为温带丛生禾草典型草原，旱地主要耕种耐寒经济作物、落叶果树。

图1 研究区地理位置

1.2 基础数据处理

遥感影像通过GEE平台(https://code.earthengine.google.com)获取，数据选用2001-2020年Landsat系列卫星地表反射率(surface reflectance,SR)，空间分辨率30 m，时间分辨率16 d。为避免Landsat 7卫星数据产生条带影响，2001-2012年和2013-2020年分别使用Landsat 5和Landsat 8卫星数据。Landsat数据获取与预处理主要包括时间与云量控制、去云处理、指数计算与最大值合成。考虑季节因素对植被覆盖特征的影响，获取研究区各年6-10月影像，按云量<10%筛选，后利用最大值合成法获取各年的NDVI最大值。

气象站点数据来源于中国气象数据网(http://www.nmic.cn)，获取麦积台站2001-2020年逐日气温与降水，用于分析近20 a气温与降水变化。气象格点数据采用欧洲中期天气预报中心(https://www.ecmwf.int)的ERA5 0.125°×0.125°逐日格点再分析资料，经python语言程序批量处理后重采样生成1 km栅格数据，用于获得年积温与年降水。数字高程模型(DEM)来自于SRTM V4.1数据集，空间分辨率经重采样生成1 km栅格数据，用于获取研究区高程、坡度、坡向信息；地貌、土壤、土地利用(LUCC)数据来自于中国科学院资源环境科学数据中心(https://www.resdc.cn)1∶1 000 000地貌类型、土壤类型、土地利用空间分布数据集，空间分辨率均为1 km。以上数据经过批量裁剪与投影变换，利用地理探测器进行驱动因素分析。

1.3 研究方法

1.3.1 植被指数最大合成 NDVI是反映植被生长状态及空间分布的最佳指示因子[19]，研究植被指数的年际变化，采用国际通用的最大值合成方法(maximum value compositing,MVC)[20]，第i年的最大植被指数为

(1)

根据《土壤侵蚀分类分级标准》[21]，将植被覆盖分为：低覆盖(NDVI≤0.1)、较低(0.10.7)5级。

1.3.2 植被覆盖变化趋势 Theil-Sen Median趋势分析和Mann-Kendall检验方法结合可模拟单个像元的变化趋势，是表征长时间序列区域格局演变规律的重要方法[22]。该方法对数据误差具有较强的抵抗能力，对显著性水平的检验具有统计学理论基础，能客观表示长时间序列植被覆盖的演化趋势，计算公式

(2)

式中:SNDVI是一元线性拟合方程的斜率；NDVIi,j分别为第i和j年的最大NDVI值。当SNDVI>0时，表示植被覆盖呈改善趋势，反之，呈退化趋势。

Mann-Kendall检验，用来判断趋势的显著性。计算公式如下

设定{NDVIi}，i=2001，2002，…，2020

(3)

(4)

式中:NDVIi与NDVIj分别表示像元中第i和j年的最大NDVI值，n表示时间序列的长度，sgn为符号函数，统计量Z的取值范围为(-∞,+∞)。在给定显著性水平α下，当|Z|>u1-α/2时，表示研究序列在α水平上存在显著的变化。本文判断在α=0.05的水平上NDVI时间序列变化趋势的显著性，根据SNDVI实际情况，将SNDVI介于-0.05～0.05的区域划分为基本不变，SNDVI>0.05的区域表示植被改善，SNDVI<-0.05的区域表示植被退化。Mann-Kendall检验在95%置信度上的显著性结果分为显著变化(|Z|>1.96)和不显著变化(|Z|<1.96)。

1.3.3 Hurst指数 Hurst指数是定量描述时间序列上的持续性的有效方法，在生态学与气候学得到广泛应用。文中使用R/S分析法计算Hurst指数，以分析NDVI的持续性特征。对于时间序列{NDVI(t)}，t=1，2，…，n，定义均值序列，公式如下[23]

(5)

(6)

极差R(τ)=maxX(t,τ)-minX(t,τ)τ=1,2,…,n

(7)

τ=1,2,…,n

(8)

对于比值R(τ)/S(τ)R/S，如果存在R/S∝τH的关系，表明研究的时间序列内存在Hurst现象，H为Hurst指数，根据H判断NDVI序列是随机变化或存在持续性。Hurst指数取值范围为0～1，若0.5

1.3.4 地理探测器地理探测器是探测空间分异性，揭示驱动因素的统计学方法[25]。为定量分析天水市NDVI时空变化驱动机制，选取8个(年积温、年降水、土地利用类型、植被类型、土壤类型、高程、坡度、坡向)影响因子，利用因子探测、交互作用探测和生态探测分析各因子的贡献率。

因子探测可对NDVI空间分异性进行探测，计算公式如下

(9)

式中:q表示驱动因子(X)对NDVI(Y)的解释力，h为Y或X的分层，Nh和N为分层h和全区Y的单元数，σh2和σ2为二者的方差。

交互作用探测可在单个因子基础上，对气象、人类活动及地表等因素间的交互作用进行分析，用于判断因子间共同作用(减弱或增强)和相互独立作用对NDVI的解释力。

生态探测用于比较2个因子X1和X2对属性Y的空间分布的影响是否存在显著差异，以F统计量来衡量：

(10)

(11)

式中:Nx1，Nx2分别为X1和X2的样本量，L1，L2为二者的分层数目，SSWx1和SSWx2为X1和X2分层的层内方差之和。

2 结果与分析

2.1 植被覆盖时间变化特征

2001-2020年天水市NDVI总体呈显著波动增长趋势(图2)，在呈现波动上升趋势，增长速率为0.97%/a(P<0.05)。多年植被覆盖介于0.46～0.65，均值为0.55，其中NDVI最大值出现在2019、2020年，最小值出现在2003年。植被覆盖在2011年之前呈现缓慢波动上升态势，其原因是我国天然林保护工程及退耕还林工程的相继实施。2011-2013年植被覆盖表现出突增现象，可能原因是天然林保护工程和退耕还林工程中的人造林不断生长，使得植被覆盖有效提升。2013年后植被覆盖的波动变化，可能与气候变化以及局部地区的复垦存在关系。近20 a降水和积温均呈现上升趋势，表明天水市植被覆盖在良好的水热条件下呈现增长趋势。

图2 2001-2020年天水市NDVI、降水及积温年际变化

2.2 植被覆盖空间分布特征

2001-2020年天水市植被覆盖空间分布特征存在明显差异，如图3a所示，呈现西北低、东南高的分布特征。西北部属黄土高原丘陵沟壑区，植被覆盖较低；东南部属秦岭山地余脉，林地与农作物分布广泛，植被覆盖较高。低覆盖和较低覆盖区集中位于城镇、交通(铁路、公路)沿线及西北部渭河与散渡河河谷交汇地带；较高覆盖与高覆盖区主要分布在秦州区、麦积区东南部的小陇山林区，清水、张川县东部的关山林、草区，以及武山、甘谷县南部的西秦岭林区；中度覆盖区广泛分布于低、高覆盖区域的中间过渡地带。近20 a NDVI均值的分级统计表明，低覆盖区占全市总面积的1.5%，较低覆盖区占2.6%，中度覆盖区占42.9%，较高覆盖占28.9%，高覆盖区占24.1%。

图3 2001-2020年天水市多年平均NDVI及变化趋势空间分布

将Theil-Sen Median趋势分析和Mann-Kendall检验结果进行叠加分析，得到2001-2020年天水市NDVI变化趋势类型(图3b)。植被显著改善情况在各区县均存在，其中秦州区、麦积区西北部及清水县分布较为集中，植被的改善格局呈现由林区向其周边区域扩张的态势。植被退化区域主要集中于中部河谷地带，这与城镇规模的不断扩大存在关系，该区人类活动(如农业耕种、基础设施建设)较为频繁，导致地表植被覆盖降低，其他区域基本不变，主要以市内的高覆盖的林区为主。对植被覆盖变化类型进行统计(表1)，植被改善区域面积为9 356.49 km2，占全市总面积的65.43%；基本不变的区域为4 660.37 km2(32.59%)；植被退化区域为283.14 km2(1.98%)。

表1 NDVI变化趋势统计

2.3 植被覆盖变化持续性特征

天水市植被覆盖的Hurst指数均值为0.68，如图4所示，Hurst指数<0.5的区域占总面积的11%，其中强反持续性区域仅为0.1%，弱反持续性为10.9%；>0.5的区域占89%，其中，强持续性为33.9%，弱持续性为55.1%。以上结果表明天水市植被覆盖的持续性特征为以正向为主，小陇山林区、关山林区及西秦岭林区东部均表现为强持续性，表明在未来依然具有良好态势；小陇山、关山林区的外环表现出一定的反持续性，虽然该区域植被覆盖在过去明显好转，但反持续性表明在未来可能存在退化风险，其他区域普遍表现为弱持续性。

图4 2001-2020年天水市NDVI变化持续性

2.4 植被覆盖空间分异的驱动因素

2.4.1 驱动因子探测分析驱动因子探测结果显示(表2)，不同因子对NDVI空间分异的解释力为降水>土壤类型>积温>土地利用类型>地貌类型>高程>坡度>坡向。所有因子中降水、土壤类型、温度和土地利用类型对NDVI的解释力较强，q均在30%以上，其中，降水的解释力为56.8%，表明处于干旱半干旱交汇区的天水市，降水是影响NDVI的主导因素。土壤类型、积温和土地利用类型的解释力分别为38.5%、37%和35.4%，为次要的驱动因子；地貌类型、高程、坡度以及坡向的解释力较低，表明其对植被覆盖的解释力不足。

表2 因子探测结果

2.4.2 交互作用探测分析如表3所示，交互作用解释力最强的是降水和土地利用，2个因子交互作用的q为0.702，表明土地利用类型与降水量相互影响，并作为因变量解释NDVI的变化较强。其次，降水与土壤类型交互作用的解释力也较强，其q为0.681，降水与积温交互作用次之；此外，降水与多数因子交互产生的q均较高，最低值是降水和坡向，其q也达到了0.57。解释力最差的是高程与坡向组合，其q仅为0.222，2个因子交互作用时解释力不足。

表3 交互作用探测结果

2.4.3 土地利用转移矩阵利用2000和2020年两期土地利用数据(LUCC)叠加分析得到土地利用类型转移矩阵(表4)，结果显示，近20 a土地利用变化类型主要为耕地→森林、耕地→草地及草地→森林，面积变化分别占全市总面积的2.45%、3.6%、2.45%。由此可见，近20 a天然林保护、退耕还林还草工程的实施，对植被覆盖的增长具有贡献。

表4 土地利用转移矩阵

3 结论与讨论

3.1 讨论

本研究结果客观反映了天水市植被覆盖状况，多年均值为0.55，空间分布格局西北低、东南高。基于Theil-Sen Median趋势分析和Mann-Kendall检验的结果表明植被覆盖整体向好，地理探测器分析结果显示，气候因素、地表因素和人类活动共同作用影响了天水市植被变化，其中降水的单因子解释力最强，且降水与其余因子交互作用解释力均较高，由此可见，降水是植被变化的主要驱动因素，与前人的研究结果一致[15,26-28]。土壤类型是NDVI变化的重要因素，在降水为主要驱动因素的区域，土壤类型对植被生长和雨水再利用效率影响显著[29]，天水市内主要的半淋溶土分布区NDVI最高，初育土分布的黄土丘陵区域植被覆盖度较低，该区域也是水土保持与生态修复的重点区域，因此NDVI增速较快。气温作为另一重要的气候因素，与降水共同作用为植被生长提供适宜的水热条件，是区域植被的分异和变化的重要因子[24,30]。近20 a随着天水市降水增加为植被生长提供充足水分的同时，气温上升有效促进光合作用和根系对水肥的吸收、传输，利于植物碳水化合物的形成。此外，地貌类型的差异也对植被的分布存在影响，由于天水市地处黄河流域和长江流域的交错地带，地貌类型复杂多样，市内东部与南部地层褶皱隆起形成的山地地貌[31]，东北部的关山为陕、甘两省重要的天然林区和草场是黄土高原地区重要的生态屏障，南部的小陇山与西秦岭作为秦岭余脉是黄河、长江两大流域的重要分水岭；市区西北部受地质沉陷和黄土层沉积形成黄土丘陵地貌[31]，土质瘠薄，易造成水土流失。

除自然因素外，人类活动不断改造着地表环境，是影响植被覆盖变化的重要驱动因素。人类活动既对植被覆盖增长起促进作用，也会破坏植被覆盖使其退化[32]。自20世纪90年代，受气候暖干、人为因素的影响，天水地区作物受干旱威胁导致农业产量出现较大波动，天然林区亦遭受严重影响[33]。1998年国家林业局试点启动12个省区的天然林保护工程，1999年实施退耕还林工程，以遏制生态恶化、保护生物多样性，促进天然林资源保护。近年植被覆盖的改善与天然林保护、退耕还林还草等生态工程实施及各级部门的监管紧密相关，据资料显示，1999-2014年，天水市退耕还林面积23.4万hm2，随着2014年新一轮退耕还林工程的启动，退耕还林面积有所提升[34]。此外，天水市南北两山绿化提升改造工程对秦州—麦积南北两山植被覆盖的提升效果显著。因此，黄土高原地区植被覆盖的增加，也是人类活动积极作用的结果，但在人类活动剧烈影响区(如城镇扩张、市政工程与交通线建设)，植被覆盖亦存在退化[35]。Hurst指数结果表明，天水市植被覆盖在未来以正向持续性为主，部分反持续性区域(黄土丘陵区、林区外围部分区域)由于缺少了林地的保护作用，不能有效防水固土，极易造成水土流失，可能导致植被退化，需保持关注。