高寒草甸是在寒冷湿润的气候条件下, 以耐寒的多年生草本植物为建群种的高寒草地生态系统类型, 其广泛分布于青藏高原^[1]。高寒草甸是高寒地区重要的生态安全屏障和畜牧业基础^[2]。由于自身生境脆弱及长期人为活动干扰, 青藏高原上的高寒草甸发生了不同程度的退化^[3]。

多年的退化恢复实践表明, 需基于“分区-分类-分级-分段”的技术体系进行青藏高原高寒草地的退化治理^[4]。在这个技术体系中, 草地退化分级是一个重要的问题。通常, 随退化程度增加, 草甸植物群落盖度、高度、物种数、优良牧草比例逐渐下降^[5-6]。因此, 在尺度较小的研究中, 主要通过植物群落结构指标^[7-8](包括物种多样性、物种组成等)和群落功能指标(主要是植物群落生物量^[9-10])对草甸植物群落的退化程度进行界定。在空间尺度较大的研究中, 也常通过遥感影像计算草地的NDVI, 并反演草地植被盖度或生物量, 进而判断草地的退化程度^[11-16]。

如何定量、客观地识别退化阈值是草地退化阶段划分研究需要解决的关键问题^[17]。目前多根据经验和差异分析(如方差分析), 将退化草甸划分为重度退化、中度退化和轻度退化等级别^[18-21]。虽然这个分级体系没有明确统一的退化阈值识别方法, 但其在大范围的退化草甸修复治理中发挥了重要作用。草甸植物群落退化是一个非线性的过程, 且群落结构退化和功能退化可能也并不同步(如物种数下降, 但生物量增加)^[22]。因此, 在对退化草地进行小规模的精准修复时, 如果仅将草地退化划分为重度、中度和轻度等阶段, 无法准确处理处于各退化阶段之间或群落各指标退化程度矛盾的情形。本研究以滇西北典型高寒草甸为研究对象, 采用Mann-Kendall突变分析对野外调查获得的样方数据进行植物群落结构指标和功能指标退化阈值的识别, 并基于退化阈值的组合情况划分退化阶段, 以期解决上述关键问题, 进而完善草地退化阶段划分体系、为退化草地的精准治理提供依据。

1 材料与方法 1.1 研究区概况

研究区位于云南省迪庆藏族自治州香格里拉市区西南约7 km处的石卡雪山山麓(99°38′49.7″E, 27°48′03″N, 平均海拔3310 m)。香格里拉市属寒温带山地季风气候, 年均温6.9℃, 日温差大, 太阳辐射强烈; 年均降雨量619.5 mm, 干湿季分明, 全年71%以上的降水量集中在6—9月^[23]。由于气候暖化和旅游干扰影响^[23-24], 研究区内草甸表现出不同程度的退化, 其中, 旅游踩踏和车辆辗轧是草甸退化的主要原因^[25]。研究区中受干扰较小的草甸, 其植物群落以禾本科(Gramineae)和莎草科(Cyperaceae)植物为优势种, 主要伴生种有矮紫苞鸢尾(lris ruthenicar)、马蹄黄(Spenceria ramalana)和坚杆火绒草(Leontopodium franchetⅡ)等。

1.2 数据获取与分析 1.2.1 植物群落野外调查和室内分析

于2018年7月中旬, 即香格里拉市高寒草甸植物生长茂盛、受旅游践踏干扰影响最集中的阶段^{[23, 26]}, 在研究区内选取较少受旅游干扰(即人为活动少、没有道路分布)的区域作为“健康草甸”的参照, 随机设置18个1 m×1 m的调查样方; 在有旅游干扰的区域, 依据道路宽度判断旅游干扰强度, 在不同干扰强度区域随机布设42个1 m×1 m的样方进行植物群落调查。在每个1 m×1 m的样方内记录出现植物的名称、盖度和高度, 然后齐地刈割样方内植物的地上部分。回到实验室后, 将收集的植物样品放入烘箱烘干至恒重并称量, 得到植物群落的地上生物量。

1.2.2 植物群落数据地理编码

采用Mann-Kendall非参数检验法对植物群落调查数据进行突变检验。Mann-Kendall突变检验的数据需满足独立性和连续性要求^[27-28], 野外随机采样得到的植物样方数据是相互独立的, 但缺乏连续性。本研究通过地理编码的方式解决调查样方连续性问题, 地理编码过程如下：1)将研究区域划分成若干网格单元; 2)在网格内, 从出现第一个调查样方的位置沿某个方向开始连续编号。由于从不同的方向对调查样方进行编码得到的数据序列并不相同, 因此突变分析得到的突变点也可能不同。为验证从不同方向编码得到的数据序列具有类似的突变点, 本研究分别按由南至北、由北至南、由西至东、由东至西四个方向对调查样方进行地理编码。图 1所示的是从北至南的编码过程, 根据调查样方出现的位置依次设置编号为G-1、G-2、......。其他3个方向的编码过程与此类似。对调查样方完成编码后, 相应地, 每个样方的植物群落指标也获得对应的编码信息。

1.2.3 植物群落数据突变分析

(1) 植物群落结构和功能指标构建

将物种丰富度(R)和优良牧草综合优势度(P_SDR)作为群落结构的指标。R通过统计调查样方内出现的植物种数得到。P_SDR的计算包括两个步骤：1)计算植物种在群落中的综合优势度^[29-30], 计算公式为：SDR_i=(RC_i+RH_i)/2, 式中SDR_i为物种i的综合优势度, RC_i为相对盖度, 其等于物种i的盖度与群落总盖度(所有物种盖度相加)的比值; RH_i为相对高度, 其等于物种i的高度与群落总高度(所有物种高度相加)的比值; 2)将禾本科植物和莎草科植物的SDR相加, 得到P_SDR。

将植物群落地上生物量(AGB)作为衡量植物群落功能的指标。

(2) Mann-Kendall突变分析

对地理编码得到的植物群落数据序列进行Mann-Kendall突变分析。对某个群落指标x, 其Mann-Kendall突变分析过程如下：

设样本序列为x, n为序列的样本长度(即样方数量), 令m_i为第i个样方x_i大于第j个样方x_j 1≤j≤i的累计数, 假定序列无变化趋势, 可以定义统计量：

假定序列为平稳正态相互独立随机变量序列, 则可得d_k的均值和方差分别为：

将d_k标准化：

得到一条UF的曲线, 随后按样本序列x的逆序x_n, x_n-1, …, x₁再重复上述过程, 得到UB, 给定显著性水平α =0.05, 若UF和UB在置信线之间出现交点, 那么交点便是突变的开始^[31]。

本研究采用“嵌套”Mann-Kendall分析的方式得到一个数据序列的多个突变点, 即：当数据序列检测到第一个突变点之后, 以该突变点为分界, 将数据序列分为小于该值和大于该值的两部分, 对这两部分数据分别进行突变分析, 此过程重复进行, 直至无新突变点出现。

1.2.4 草甸植物群落退化程度划分

将植物群落结构指标和功能指标的突变点整合起来划分草甸植物群落的退化程度、绘制草甸植物群落退化程度查询表。草甸植物群落退化程度查询表的绘制过程中, 按V_ij/V_maxi(式中V_ij为指标i的第j个突变阈值, V_maxi为指标i的最大突变阈值)的方式对坐标轴数值进行标准化。在Origin 9.1中完成作图。

2 结果与分析 2.1 植物物种丰富度的变化阈值

Mann-Kendall突变检测发现, 四组植物群落地理编码数据序列的物种丰富度均存在多个突变点。从南至北编码方向上, 突变点为10、13、14、15、16、17和20种; 从北至南编码方向上, 突变点为10、13、14、15、16和17种; 从西至东编码方向上, 突变点为10、11、13和17种; 从东至西编码方向上, 突变点为10、11、13和17种(图 2)。四组地理编码序列共有的突变点为10、13和17种, 将这3个数值作为植物物种丰富度的变化阈值。

2.2 优良牧草综合优势度的变化阈值

四组地理编码序列的优良牧草综合优势度均存在多个突变点。从南至北编码方向上, 突变点为0.18、0.25、0.30、0.32、0.36、0.39、0.45和0.59;从北至南编码方向上, 突变点为0.18、0.25、0.30、0.32、0.39、0.46和0.54;从西至东编码方向上, 突变点为0.18、0.25、0.30、0.34、0.39、0.45和0.60;从东至西编码方向上, 突变点有0.22、0.25、0.30、0.39、0.46和0.60(图 3)。四组地理编码序列共有的突变点为：0.25、0.30和0.39。将这3个数值作为优良牧草综合优势度的变化阈值。

2.3 植物群落地上生物量的变化阈值

四组地理编码序列的植物群落地上生物量均存在多个突变点。从南至北编码方向上, 地上生物量的突变点为31.98、63.77、74.42、91.68、95.00、111.08、124.32、132.36 g/m²; 从北至南编码方向上, 突变点为25.61、31.98、38.56、57.21、63.77、74.42、91.68、95.00、120.89 g/m²; 从西至东编码方向上, 突变点为25.61、31.98、86.64、91.68、95.00、106.97、124.32 g/m²; 从东至西编码方向上, 突变点为31.98、83.78、91.68、106.97、116.5 g/m²(图 4)。四组地理编码序列共同的突变点为31.98 g/m²和91.68 g/m², 将这两个数值作为植物群落地上生物量的变化阈值。

2.4 草甸植物群落退化程度的划分结果

根据植物物种丰富度(R)、优良牧草综合优势度(P_SDR)和地上生物量(AGB)的变化阈值可识别出草甸植物群落退化的阈值空间(图 5), 根据该阈值空间可将草甸退化定义为两种类型：“基本退化类型”和“过渡退化类型”(图 6)。

本研究将植物群落结构指标和功能指标的数值均处在最大变化阈值以上的草甸定义为健康草甸, 即：R≥17种、P_SDR≥0.39、AGB≥91.68 g/m²。

当植物群落的结构和功能指标数值低于健康草甸时, 认为草甸发生了退化。根据变化阈值, “基本退化类型”可划分为三类：Ⅰ级退化, 即13种≤R < 17种、0.30≤P_SDR < 0.39、31.98 g/m²≤AGB < 91.68 g/m²; Ⅱ级退化, 即10种≤R < 13种、0.25≤P_SDR < 0.30、31.98 g/m²≤AGB < 91.68 g/m²; Ⅲ级退化, 即R < 10种、P_SDR < 0.25、AGB < 31.98 g/m²(表 1, 图 6)。

依据本研究采样点的数据情况, 通过植物群落退化程度查询表得到介于“基本退化类型”之间的4个“过渡退化类型”分别为：Ⅰ级正向过渡, 即Ⅰ级退化向健康草甸过渡阶段; Ⅰ级反向过渡阶段, 即Ⅰ级退化向Ⅱ级退化过渡阶段; Ⅱ级反向过渡阶段, 即Ⅱ级退化向Ⅲ级退化过渡阶段; Ⅲ级正向过渡阶段, 即Ⅲ级退化向Ⅱ级退化过渡阶段(表 1)。

3 讨论与结论 3.1 基于地理编码和突变分析识别植物群落退化阈值的可行性

识别退化阈值才能准确划分草地退化阶段^[17]。目前主要通过模型模拟和控制放牧实验^[32]的方式分析植物群落的变化阈值。模型模拟方面, 常用草原群落退化指数(SCDI)^[33]、CENTURY模型^[34]、MaxEnt模型和ROC曲线^[35-36]等方法, 通过模拟环境因子变化来研究生态系统的变化阈值, 但这些模型的构建及阈值计算过程复杂、模型运行结果易受生态系统中不确定因素的影响, 且复杂模型存在过度拟合问题, 会对模拟结果产生不良影响^[37-39]。控制放牧实验方法确定的退化阈值较可靠^[32], 但对实验时间和实验场地要求高, 且不适用于非放牧干扰导致的退化情况。相较于控制放牧实验, 野外样方调查数据的获取成本更低、效率更高。本研究通过野外样方调查数据识别植物群落结构指标和功能指标的退化阈值, 并基于退化阈值的组合特征划分研究区草甸退化类型, 丰富了草地退化阶段划分研究的方法。

从阈值识别结果来看, 本研究从4个方向对植物群落数据进行地理编码, 经突变分析所得的变化阈值较接近, 说明基于该方法进行阈值识别在技术层面是可行的, 得到的结果是可靠的。针对草甸植物群落结构退化过程和功能退化过程不一致的现象, 本研究通过建立植物群落结构指标和功能指标退化阈值之间的对应关系, 在识别“基本退化类型”的基础上, 识别出了介于3个基本退化类型之间的“过渡退化类型”, 解决了传统分类方法不能统一各单项指标划分结果相互矛盾的问题。

需要指出的是, 由于生态系统的层次性特征, 各生态要素特征会在不同时空尺度间存在差异^[40-42]。因此, 对草甸植物群落退化阈值的识别以及退化程度的划分也应在特定的时空范围内进行。本研究建立的草甸植物群落退化阈值识别方法只适用于某块草甸或同种草甸类型(即：草甸异质性不高)的监测, 不能用于异质性高的不同草甸地块或多种草甸类型混合的情形。其次, 作为参照的“健康草甸”的定义标准决定着草甸植物群落退化程度的划分结果, 在实际监测过程中, 需要根据当地情况确定“健康草甸”的内涵。

3.2 识别植物群落退化的“基本退化类型”和“过渡退化类型”的意义

草甸退化是一个复杂的生态过程, 但在干扰条件稳定且持续存在的条件下, 随退化程度加重, 植物群落表现出结构简单化、群落功能弱化的趋势：群落物种丰富度、优良牧草比例、地上生物量等表征指标逐渐降低, 且在不同退化程度间存在明显差异^{[3, 43]}。因此, 可将植物群落结构和功能指标变化的突变点(即退化阈值)组合起来作为“基本退化类型”的划分依据, 这种退化类型明确指出草甸所处的退化阶段, 从而为退化植物群落的恢复措施制定、恢复方向确定提供依据。

在自然条件下, 气候条件等因素时而有利于草甸恢复, 时而加剧草甸退化进程。因此, 退化草甸植物群落的变化也有两个方向：加剧退化或逐渐恢复——这两个变化过程可能是交替出现的。本研究将介于“基本退化类型”之间的植物群落定义为“过渡退化类型”, 归为该类型的植物群落存在两种可能的变化趋势：由较低的退化程度向较高的退化程度过渡、由较高的退化程度向较低的退化程度过渡。在建立退化程度查询表的基础上, 如果对某块草甸进行长期监测, 就可判断出退化区域草甸植物群落的“过渡趋势”, 进而为恢复措施的适时调整以及恢复效果评估提供更精准的依据。