植被在地球上的能量交换和生态循环中起着重要作用^[1]，对于生态脆弱的草原矿区而言，植被受气候、采矿、放牧等人类活动的影响较大。利用多光谱遥感数据得到的植被覆盖度(fraction of vegetation coverage, FVC)作为衡量植被状况、指示植被动态变化的指标^[2-4]应用广泛。文献[5—9]分别对沙漠化地区、阴山北麓草原区、内蒙古大青山自然保护区、甘南草原区、淮南矿区等地的植被覆盖度动态变化进行了研究。从方法上看，利用线性混合理论对混合光谱进行应用和解释是较为普遍的方法，像元二分模型在构建时NDVIveg(NDVI for “pure” vegetation pixel)和NDVIsoil(NDVI for “pure” bare soil pixel)的取值也得到重视^[10-12]。在典型区域长时序植被覆盖度变化研究中传统的静态多时相比较研究仍有应用^[13-14]，但此种方法对于草原区FVC变化分析时易受到气温、降水等年际物候差异等的影响，难以科学揭示生态是否退化及其驱动因素等根本问题。已有的研究中对FVC进行了分级^[15-16]，这些基于分级的分析能够在一定程度上表现区域植被演化特征，但难以进一步提高定量表达的精度。此外，多数长时序植被覆盖度研究中只采用一种趋势分析方法，如回归分析、Sen趋势度和Manner-Kendall趋势检验相结合的方法等^[17-19]，多种分析方法的对比研究较少。

针对上述问题，本文以内蒙古呼伦贝尔市宝日希勒露天矿区为例，基于1985—2015年年度最大合成NDVI数据，采用一元线性回归法和Sen+Mann-Kendall法，进行FVC的长时间序列分析，并对比两种趋势分析方法的差异。研究结果有助于获取科学评价长时序煤炭开发活动对地表生态影响的基础数据，同时为类似区域长时序植被变化监测及FVC变化趋势分析提供方法参考。

1 研究区概况与数据获取 1.1 研究区概况

宝日希勒矿区位于内蒙古自治区呼伦贝尔市陈巴尔虎旗，地理坐标为49°10′N—49°40′N，119°10′E—120°30′E，矿区面积约为255 km²。

考虑到矿区开采可能影响周边区域，本文以自然地形界线和行政界线划定研究边界，即以莫尔格勒河和海拉尔河为研究区北、西、南的界线，东侧以陈巴尔虎旗和海拉尔区的乡镇行政边界为界线划定研究区，总面积约为1792 km²，地势由东北向西南逐渐降低，地跨森林草原与干旱草原两个地带。

1.2 数据来源

本文原始数据为1985—2015年的Landsat TM、ETM+、OLI遥感影像数据。为消除年际间物候差异等的影响，经中科院地理科学与资源研究所陆地表层格局与模拟实验室进行影像处理后，得到各监测年年度最大合成NDVI数据。由于1993年和2012年的图像质量不佳，因此研究中使用了29 a的年度最大合成数据。

年平均气温和年降水数据来源于资源环境数据云平台(http://www.resdc.cn/)。

2 研究方法 2.1 植被覆盖度计算

采用像元二分模型计算FVC，计算方法为^[15]

式中，NDVI_soil为裸土或无植被覆盖区域的NDVI值；NDVI_veg为完全被植被所覆盖像元的NDVI值。本文中依据NDVI频率统计表，使用0.5%置信度确定NDVI值取值的置信区间。

2.2 植被覆盖度时间变化特征分析

研究区在1985—2015年第i年各像元FVC均值采用下式计算

式中，fc_i为研究区第i年的FVC的均值；fc_ij为研究区内第i年中第j个像元的FVC取值；i为年的序列号，取值为i=1，2，…，29；n为像元个数。

研究区FVC的时间变化特征分析采用转移矩阵法，其中FVC等级划分参考现有成果^[16]并结合研究区实际情况确定分级界限，见表 1。

2.3 植被覆盖度趋势变化对比分析

本文采用一元线性回归法^[20]和Sen+Mann-Kendall法^[21]进行FVC趋势变化分析。采用上述两种方法结合进行研究区的植被覆盖度趋势分析可以增强检验的抗噪性，并在一定程度上提高检验结果的准确性^[22]。

2.3.1 一元线性回归法

一元线性回归法计算如下^[4]

式中，n为监测时间。当slope>0时，说明FVC呈增加趋势，反之降低。

一元线性回归法的趋势检验采用t检验法^[23]，计算方法如下

式中，代表一元线性回归法得到的系数slope；i=1，2，3，…，29为年序号；n=29；，其中y_i为第i年的FVC的实际值，为第i年的FVC的回归值；x_i为时间自变量；α表示显著性水平。

2.3.2 Sen+Mann-Kendall法

Sen趋势度计算如下^[17]

式中, 1 < i < j < n；β代表斜率，β>0表示时间序列呈上升趋势，反之下降。

Mann-Kendall 趋势检验统计量S的计算如下^[19]

根据时间序列长度n值大小的不同，显著性检验统计量的选取有所不同。本文时间序列长度为29，采用检验统计量Z来进行趋势检验，检验中取显著水平α=0.05，Z_1-α/2=Z_0.975=1.96。

Z值计算方法如下

式中，，其中n为序列中数据个数，m为序列中结(重复出现的数据组)的个数，t_i为结的宽度(第i组重复数据组中的重复数据个数)。

3 结果与分析 3.1 植被覆盖度总体变化和多时相比较

研究区1985—2015年FVC均值的年际变化情况如图 1所示，FVC空间分布如图 2所示。

FVC的线性趋势不显著，fc整体呈波动微弱下降趋势，在1988、1991、1999和2008年出现峰值，分别为0.688、0.708、0.708、0.709；谷值出现在2007年，为0.434；均值为0.585，中位数为0.562。

图 3选取FVC均值接近的1985、1995、2005和2015年4期数据进行多时相FVC比较，并对1985年和2015年FVC作转移矩阵。结果表明：研究区FVC年际差异大，1985—2015年高覆盖的面积和中高覆盖的面积共减少251.03 km²，其中157.23 km²变为低值被覆盖区域，52.87 km²变为无植被覆盖区域，占总转移面积的83.70%，FVC呈明显下降趋势，FVC低值区域也在不断变化。通过气象数据可知，1985、1995、2005和2015年研究区年平均气温分别为-1.814、1.312、0.211、0.743℃，年降水量分别为472.138、349.816、410.352、301.259 mm。结合图 2和气象数据可得，FVC空间分布情况受气温和降水影响大，多时相的静态对比无法准确揭示草原的长时序变化趋势，据此判断草原是否退化等容易引起误判。相比之下，遥感时序分析方法将时序影像视为连续时态记录的完整档案，可以更为精确地获取地物在时序上的变化趋势和规律^[24-26]。

3.2 植被覆盖度空间演化趋势分析

1985—2015年FVC变化趋势程度如图 3所示。一元线性回归法和Sen+Mann-Kendall法趋势检验分析结果见表 2。

从图 3可得FVC下降最集中区域为研究区中部宝日希勒矿区所在位置，其他大部分区域的FVC变化斜率在0值左右，两种趋势分析方法得到的研究区FVC变化程度基本一致，但Sen+Mann-Kendall法在植被退化和植被改善区域面积有所增加。一元线性回归法分析结果表明，植被极显著增加的区域占1.95%，植被极显著减少的区域占3.35%，48.83%的地区植被变化不显著；Sen+Mann-Kendall法趋势检验分析结果表明，植被严重退化区域占6.50%，40.51%的地区植被变化不明显，植被明显改善区域占3.32%。

两种趋势分析法分析结果的空间分布差异如图 4所示。

由图 4可得，Sen趋势度法和一元线性回归法所模拟的植被覆盖度变化趋势基本上是一致的，两者之差在(-0.035 5，0.023 6)，绝大多数区域无差异。两种趋势分析法得出植被变化趋势相同的像元占比78.10%，Sen+Mann-Kendall法较一元线性回归法而言，植被退化面积占比和植被改善面积占比分别多了4.65%、3.67%；负值和正值反映出Sen趋势度法在植被退化和植被改善上要比一元线性回归法更为敏感，克服了一元线性回归法易受噪声影响的缺点。Sen+Mann-Kendall方法对于植被退化和植被改善的选择更为敏感，植被变化不显著的区域减少，结果更加精确。两种方法得出的结论基本一致，均为30年间研究区的植被覆盖度的改善面积大于植被覆盖度的退化面积；在草原矿区长时序植被覆盖演化趋势的分析中，Sen+Mann-Kendall方法能更准确地反映其变化趋势。

4 结论

本文以宝日希勒矿区及周边为研究区，应用1985—2015年的Landsat TM、ETM+、OLI遥感影像数据进行FVC的变化趋势和空间差异分析，主要结论如下：

(1) 研究区平均FVC的线性趋势并不十分显著，整体呈波动下降趋势。

(2) 采用两种趋势方法分析30年间FVC演化趋势，两者结果基本一致，均为研究区中部宝日希勒矿区所在位置FVC退化明显。

(3) 综合比较两种FVC趋势分析方法的结果，Sen+Mann-Kendall法较一元线性回归法而言，植被退化面积占比和植被改善面积占比分别多4.65%、3.67%，该方法在FVC改善和退化的反映上更为敏感，植被变化不显著的区域减少，能监测出轻微的植被变化趋势，趋势变化结果更加精确，更能准确地反映草原矿区长时序的FVC演化趋势。