大九湖泥炭藓沼泽植被指数时空变化研究

张 唯，张振泽，刘福江*，梁天欣

（1. 中国地质大学（武汉）地理与信息工程学院，湖北 武汉 430000）

泥炭藓沼泽（sphagnum swamp）经过长时间的生长积累，可以将大气碳以平均每年每平方米12～23 g的速度固定[1]。尽管泥炭藓沼泽仅仅拥有相当于地球陆地表面的3%的覆盖面积[2]，但其储存了地球上约1/3的碳源[3]。而泥炭藓(Sphagnum)作为泥炭藓沼泽的优势物种，每年固定的碳[4]要远多于其他世界陆生植物积累的碳[5]。赵素婷[6]等于2011年中旬对大鄂西亚高山区泥炭藓沼泽开展了深入研究，并考察了该地区的高级植被群类型多样性，对大九湖泥碳藓沼泽在生态系统稳定领域的重要作用，以及主要植被群丛种类等做出了判断和研究。因此对泥炭藓沼泽植被生长趋势以及长时间序列变化监测是了解气候变化对泥炭藓沼泽影响、预估未来泥炭藓沼泽如何变化、检测生态系统是否稳定的有效工具［7］。

归一化植被指数（NDVI）作为应用最广泛的指数之一，能对较大空间范围内的植被生长状况进行监测。Boelman[8]等以由泥炭藓组成的苔原植物种群为主要研究目标，研究成果指出，NDVI在控制南北极苔原植物的生物量年际变动过程中具有优异的作用。同时，NDVI 是评价植被生长状况的重要指标，对气候变化极为敏感[9-10]。利用NDVI在时间和空间尺度上监测植被生长状况并研究其对气温和降水的响应已成为当前研究的热点之一。马蓉[11]等研究发现中国温带草原植被生长季NDVI 在1982—2015年总体呈增加趋势；生长季内，降水是最重要的影响因子。袁雷[12]等研究发现羌塘国家级自然保护区生长季NDVI 随时间变化呈变好趋势，气温与降水对NDVI 都有重要的正向影响。因此，有必要针对不同的研究区域，具体分析其NDVI与温度、降水等气象条件的相关性变化特征。

大九湖泥地区随着人类活动的加剧和气候变化，使得在1986—2005年该地区的生态环境受到了影响，在近十几年该地区政府出台了一系列的保护政策，使得该地区生态得到改善。在上述长时间的变化中，对泥炭藓沼泽植被生长状况发生的变化研究仍然匮乏；对该地区植被指数变化对气温、降水响应特征的认识不系统。基于该现状，本文利用Landsat 系列数据获取时序NDVI 特征，并结合气温、降水资料，研究了1986—2021年大九湖泥炭藓沼泽地区的时空变化特征，为指导后续对泥炭藓沼泽植被保护提供科学依据，为该区域生态环境发展提供理论支持。

1 研究区概况

大九湖湿地公园（109°56～110°11'E、31°34'～31°33'N）为亚高山湿地公园，位于神农架国家级自然保护区西部。由于研究区附近隆起地形的阻隔作用，以及长时间以来的植被发育，在研究区内形成了湖北省特有的、保存最为完整的、全球罕见的中纬度亚高山沼泽湿地，主要包括湿生泥炭沼泽、湿生草本沼泽以及湖泊等湿地类型[13]。其中，泥炭层深达2 m，尤其是研究区主要地带（110°0′ 09.91″ E，31°29′ 30.2″ N）存在着覆盖度接近百分百的泥炭藓。

当地政府在1986年至2001年对湿地内的天然河道进行了矫正，加大了疏通落水孔的力度，并且实施了开沟排水工程，开通了主要的南北运河。2002年以来，当地政府的招商引资导致外来商人将大九湖区域开发为各种经济作物的种植基地，使得湿地生态遭到了严重破坏，泥炭藓难以再度生长[14]，最终导致沼泽的退化。

2006年中旬，大九湖国家湿地公园正式设立。自2008年年末开始，包括“大九湖湿地保护/修复”在内的“大九湖湿地保护/修复/公园建设项目”等一系列政策开始实施。其中包括基础设施建设、开发旅游项目、生态移民搬迁等一期工程。第一期湿地恢复工程于2010年年末完成并通过国家审批。近些年，大九湖泥炭藓湿地研究区内的农业用地大部分被其他土地利用类型所取代[15]，生态系统逐步得到恢复。

2 数据来源与处理

本文采用美国Landsat 陆地系列卫星TM、ETM、OLI 数据作为数据源，均下载自USGS 官网(https://www.usgs.gov/)，数据空间分辨率为30 m。影像行列号为125，38，时间序列为1986—2021年，选取植被生长季数据。基于ENVI软件，计算得到最大NDVI 值[16]，最后经筛选得到36 景影像。其中，矢量数据来源于OSM 官网（https://master.apis.dev.openstreetmap.org/）。气温与降水数据来源于国家科技基础条件平台国家地球系统科学数据中心(http://www.geodata.cn/)月值数据，空间分辨率为1 km。

3 研究方法

3.1 NDVI信息提取

NDVI 是最常见的植被生长动态指标之一，取值范围在-1～1 之间。大于0 表示不同程度的植被覆盖，越接近1植被生长状况越良好。其计算公式如下：

式中，NIR 为近红外波段；RED 为红光波段；利用ENVI完成对NDVI数据的提取。

3.2 Theil-Sen median趋势分析

Theil-Sen median 趋势分析为非参数统计趋势分析方法。该方法在分析过程中不要求数据服从正态分布，且不受异常值影响[17-18]。对研究区泥炭藓沼泽进行逐像元分析，其计算公式为：

式中，n 为时间序列长度；Xj和Xi为第i年和第j年的NDVI 最大值。当β>0 时表示植被指数呈上升趋势，反之则为下降趋势。

3.3 Mann-Kendall显著性检验

Mann-Kendall 检验是一种非参数检验方法，该方法同样不需要样本服从特定的分布，不受异常值干扰[19]。将数据像元构造为一组时间序列，判断显著性差异[20-21]：

式中，ZMK为单个序列的统计检验量；S统计量由NDVI 的前后影像数据计算得到；var(S)为S统计量的方差，以上为化简后的求解公式；sgn 为逻辑判别函数符号；给定显著性水平α为0.05。

3.4 相关与偏相关分析

植被生长对气温和降水的响应存在一定的滞后性[22]，采用相关分析法分别对NDVI 与前1、2、3 个月的滑动窗口的气温、降水的响应进行分析；利用偏相关分析探究研究区NDVI 对气温、降水的响应[23-24]。

式中，i为样本数；xi为第i年的植被指数值；yi为不同时序各气象要素的值；-x为植被指数平均值；-y为各要素平均值。

式中，Rxy，z为将z影响剔除后x 与y 之间的偏相关系数。

4 研究结果分析

4.1 植被指数的年际变化

根据大九湖地区长时间的发展情况，该地区在2005年存在明显的拐点，因此将变化趋势分为两部分进行研究，分别为1986—2005年开发阶段，2006—2021年维护阶段。

4.1.1 泥炭藓沼泽NDVI的空间变化

采用Theil-Sen median趋势分析与M-K统计检验方法对NDVI逐像元进行分析，得到植被变化空间分布图。本文将研究区内植被指数划分为5 种变化类型，如表1。

表1 植被指数变化类型分类

表2 1986-2005年气温和降水与NDVI平均相关系数

表3 2006-2020年气温和降水与NDVI平均相关系数

其中，1986—2005年大九湖泥炭藓沼泽植被大部分地区呈退化趋势，在整个地区西南处以及中部呈现出显著退化趋势。显著退化的面积约1.63 km2，约占研究区总面积的15.96%；微显著退化的面积约4.53 km2，占比约44.42%；2006—2021年研究区内植被整体呈改善趋势。显著改善的面积约4.7 km2，占比约46.07%；微显著改善的面积约3.61 km2，占比约35.43%；1986—2021年整体来看，研究区内中部地区有显著退化的趋势，其余地区大部分呈改善趋势。显著退化的面积约1.21 km2，占比约为11.87%；微显著退化的面积约为1.89 km2，占比约为18.58%；微显著增长的面积约为3.04 km2，占比约29.83%；显著增长的面积约3.25 km2，占比约31.83%。

4.1.2 泥炭藓斑块NDVI的年际变化

上述研究已验证研究区NDVI 空间变化趋势，为更好地研究该地区的NDVI 变化对泥炭藓单一植被影响，选取该地区经实地考察泥炭藓覆盖度超过90%的泥炭藓斑块（110°0′ 09.91″ E，31°29′ 30.2″ N）为实验斑块，变化情况如下：

1）1986—2005年研究区斑块经Theil-Sen median趋势分析得整体为下降趋势，MK 显著性检验结果为显著下降。可见人为的开发行为对泥炭藓这一植被的生长产生了影响，使其线性下降率为-0.018/(10 a)，结果如图4。其中，在1987年以及1998年出现突然下降，在1988—1990年出现了缓慢上升的趋势。

2）2006—2021年该斑块趋势分析结果呈上升趋势，MK 显著性检验结果为微显著上升，在经治理后使其线性增长率为0.046/(10 a)，其中，在2009年出现突降，在2018年达到最大值。

3）整体分析1986—2021年该斑块呈一定的上升趋势，MK 显著性检验结果为不显著上升，线性增长率为0.021/(10 a)，结果如图1所示。

图1 泥炭藓斑块NDVI趋势变化

4.2 NDVI对气温和降水的响应

论文分析了大九湖地区泥炭藓沼泽地区NDVI 在时间和空间尺度上的变化。由于NDVI 除受人为因素的影响外，同时也受到气温和降水的制约，因此有必要对研究区内NDVI 与气温和降水的响应特征进行一步研究。

4.2.1 NDVI与气温和降水的滞后效应

选取1986—2020年最大值NDVI 数据与气温和降水数据经相关分析可得，该地区1986—2005年NDVI与气温和降水呈负相关关系；2006—2020年NDVI 与气温呈负相关关系与降水呈正相关关系。其中NDVI与前一个月气温相关系数最大，与降水前3 个月相关系数最大。因此，该地区NDVI滞后于气温1个月，滞后于降水3 个月。该结论与解晗[22]等在研究黄河流域植被指数对气温和降水的响应和魏榕[23]等雅碧江流域生长季植被对气象因子的响应分析中一致。

4.2.2 NDVI与气温和降水偏相关分析

根据上述气温和降水响应的滞后结论，选取前期1 月气温、前期3 月降水与NDVI 数据进行偏相关分析。得到结论如下：

1）1986—2005年，NDVI与降水偏相关系数以负相关为主，偏相关系数小于0地区约占72.7%。偏相关系数以负相关为主的区域主要集中在研究区中部以及北部；偏相关系数以正相关为主的地区主要分布于研究区的边缘地带，并且偏相关系数普遍较高。

与气温的偏相关系数同样以负相关为主，偏相关系数小于0的地区约占57.02%。其中以负相关为主的地区分布在研究区东部边缘地带以及中部地区；以正相关为主的地区分布在北部以及西部边缘地区，边缘地带偏相关系数普遍较高。

2）在2006—2020年，NDVI与降水偏相关性以正相关为主，偏相关系数小于0的地区约占77%。其中以正相关为主的在整个研究区均有分布，偏相关系数大于0.5的地区在东部和西部分布较为集中；负相关的地区在研究区东部和中部呈现出零散分布的状态。

与气温偏相关性以负相关为主，偏相关系数小于0 的区域约占65.47%。其中以负相关为主的区域在整个研究区均有分布，小于-0.5的区域集中于西部以及北部地区。以正相关为主的区域零散分布于整个研究区，大于0.5的地区在中部地区分布较多。

5 结语

论文利用Theil-Sen median趋势分析和Mann-Kendall 显著性检验方法，分析其时序变化，并结合研究区气象资料，分析近36 a来大九湖泥炭藓沼泽植被指数的时空演变规律，得到主要结论为：

1）1986—2005年研究区植被指数主要呈退化趋势；2006—2021年研究区植被指数主要呈改善趋势。经Theil-Sen median 趋势分析和M-K 统计检验方法对NDVI 逐像元进行分析，得到植被变化结果为1986—2005年显著退化的面积约占研究区总面积的15.96%，微显著退化的面积约44.42%；2006—2021年显著改善的面积约占46.07%，微显著改善的面约占35.43%。可见自1986年开始的一系列开发活动以及2002年政府的招商引资对该地区的植被指数产生了巨大的影响，使得整个研究区60.38%区域的植被呈现出退化趋势；直到2006年政府提出建立大九湖国家湿地公园，至2021年已有约81.5%的区域植被呈改善趋势，并且优于开发之前。

2）1986—2005研究区泥炭藓斑块NDVI显著下降；2006—2021年该斑块NDVI为微显著上升。其中，经一元回归分析可得1986—2005年其线性下降率为-0.018/（10a）；2006—2021年其线性增长率为0.046/（10 a）。在研究区整体植被指数变化的趋势下，泥炭藓这单—物种也受到了的影响，其中1986—2005年受该区域开发活动影响明显，使其植被指数显著下降；而自2006—2021年，尽管大九湖政府已建立大九湖湿地公园，使得研究区大面积地区植被指数得到了显著改善，而对泥炭藓这单一物种恢复程度不如研究区其他区域。

3）该地区NDVI 对气温、降水响应存在滞后效应，分别滞后于气温1个月，滞后于降水3个月。通过对研究区NDVI与前1、2、3、个月降水数据逐像元进行相关性分析可得1986—2005年以及2006—2020年NDVI与降水前3个月相关系数绝对值最大，与气温前一个月相关系数绝对值最大，该结论与现有研究基本一致。

4）相比气温因子，1986—2020年研究区的NDVI对降水的响应更强，与降水的相关性和偏相关性也较高，与现有的研究结论基本一致。而在研究区降水的相关性更强，这与泥炭藓沼泽以泥炭藓为主要物种有着密不可分的关系，作为低等植物，无法通过根部从土壤中吸收水分，泥炭藓从而进化出依靠叶片获取水分的能力，使其对降水更为依赖，导致该区域与降水有着更高的相关性。

6 展望

本文分析了近36 a来大九湖泥炭藓沼泽植被指数的时空演变规律，能够为今后大九湖泥炭藓沼泽地区维持沼泽地区的水量平衡，结合泥炭藓的生物特性开展恢复等生态工作提供可行方向。但本文的研究尚且存在一定的局限性。首先，由于缺乏研究区实测数据，无法对研究结果进行定量的精度评价；其次，气象资料与NDVI 数据在空间分辨率上存在差异，相关分析存在一定的不确定性；最后，由于时间的关系，本文只结合相关政策分析了NDVI年际变化的趋势，后续可以结合土地利用数据具体分析人类活动对该区域植被变化的影响方式、影响规律和驱动机制。