基于Sen+Mann-Kendall的北京植被变化趋势分析

王佃来，刘文萍，黄心渊

1.北京林业大学 信息学院，北京 1000832.首钢工学院 计算机系，北京 100144

基于Sen+Mann-Kendall的北京植被变化趋势分析

王佃来1，2，刘文萍1，黄心渊1

1.北京林业大学 信息学院，北京 100083
2.首钢工学院 计算机系，北京 100144

1 引言

北京作为中国的首都和一个国际化的大都市，生态环境对其政治经济文化发展至关重要，因此全面直观准确地了解北京近14年的植被变化状况和趋势对以后生态环境评估，城市规划发展有着重要的意义，并且准确的植被变化数据可以为今后城市的整体规划与生态发展决策提供科学依据。

北京市的土地利用与土地覆盖状况近十几年来有了明显改变，许多学者近年来对北京市的植被覆盖变化做了研究。2001年陈云浩等对北京海淀区22年植被覆盖等级变化的空间过程和变化趋势进行了研究，得出海淀区植被覆盖整体呈增加趋势的结论[1]。2003年李晓琴对北京山区植被覆盖景观格局动态变化进行研究，结果表明昌平区、怀柔县、和密云县植被覆盖状况整体呈增加趋势；平谷县植被覆盖呈减少趋势；门头沟植被覆盖有转差的趋势[2]。2005年张学霞等研究了北京近五十年植被对全球变暖的响应及时效，结果表明北京市植被生长状况总体上在变好或生长季节在延长[3]。2008年张本昀等对北京山区的植被覆盖动态变化进行的研究结果表明山区覆盖度呈下降趋势[4]。2009年，郭锦等基于3S技术研究了北京市植被覆盖空间格局变化，结果表明北京市植被覆盖程度呈下降趋势[5]。祁燕等基于RS与GIS技术研究了北京市植被覆盖度变化，结果表明北京市整体植被覆盖度呈增长趋势[6]。张萍等对北京市延庆县植被覆盖变化进行了研究，研究表明延庆县植被变化明显，植被恢复远高于退化，植被总体上得到恢复[7]。总之，北京地区的植被覆盖情况近年有了明显改变，一方面通过植树造林、城市绿化、退耕还林和绿化荒山等手段，使植被活动向着增强的方面发展；另一方面，随着城市化和工业化的不断深化，植被退化趋势也比较明显，出现了沙尘暴，热岛效应等环境问题。

上述学者的研究主要通过计算各像元的覆盖度 fc，根据覆盖度取值范围将研究区域分为无植被覆盖、低植被覆盖、中植被覆盖、高植被覆盖和全植被覆盖4个等级，然后根据年际各等级的面积变化来分析植被覆盖变化情况，但对变化趋势的分析不够明确，不能直观呈现植被覆盖变化趋势。本文基于SPOT VEGETATION 1998年到2011年NDVI数据使用非参数检验法Mann-Kendall检验和Sen's Slope Estimator结合与参数检验方法一元线性回归变化斜率法对北京14年植被覆盖变化趋势进行分析，并结合GIS技术与北京市行政区图叠加，直观呈现了北京近14年植被覆盖变化时空分布格局。

2 研究区域概述

2.1 地理位置

北京市位于北纬39°56′，东经116°20′，雄踞于华北平原的北端，北以燕山山地与内蒙古高原接壤，西以太行山与山西高原毗连，东北与松辽大平原相通，东南距渤海约150 km，往南与黄淮海平原连片。全市总面积16 807.8 km2。

2.2 自然条件

北京市气候属于典型的暖温带半湿润大陆性季风气候，境内地貌复杂，山地高峰与平原之间相对高差悬殊，从而引起明显的气候垂直地带性。以海拔700～800 m为界，此界以下到平原，为暖温带半湿润季风气候；此界以上中山区为温带半湿润—半干旱季风气候；约在海拔1 600 m以上为寒温带半湿润—湿润季风气候。夏季炎热多雨，冬季寒冷干燥，春、秋短促。由于境内地形复杂，生态环境多样化，致使北京市植被种类组成丰富，植被类型多样，并且有明显垂直分布规律。自然条件对该区植被的发育是有利的，但由于北京市历史悠久，人类的生产活动对植被的结构和分布有着深刻的影响。目前平原地区主要是农田和城镇，广大山区占优势的植被是次生落叶阔叶灌丛和少量落叶阔叶林及温性针叶林[8]。

3 研究数据及研究方法

3.1 数据来源

本文使用的数据包括两部分：（1）1998年4月到2011 年12月空间分辨率为1 km的SPOT-4 VGT S10 NDVI数据，该数据下载于http：//free.vgt.vito.be/，数据区域为东南亚（SE-Asia）。为保证数据的精确性，进行以下预处理：首先使用VGTExtract软件从下载的东南亚数据根据北京地区经纬度范围大致裁剪出北京地区数据，然后结合北京地区行政区划矢量图使用ENVI软件的掩模功能裁剪出北京地区的精确遥感数据。（2）北京行政区划矢量图.以上两部分数据投影格式为Geographic Lon/Lat。

SPOT-4 VGT数据由1998年3月升空的SPOT-4太阳同步近极地轨道卫星搭载的VEGETATION（VGT）传感器产生，数据接收日期为1998年4月，该数据主要用于植被覆盖观测。SPOTVGT提供的归一化植被指数NDVI是使用红光通道反射值（B2）和近红外通道反射值（B3）计算得到，即 NDVI=(B3-B2)/(B3+B2)，其值介于 -1和+1之间。VGT数据产品分为两大类，原型数据（Prototype）和合成数据（Synthesis）。

VGT.P（Prototype）数据产品主要为科研人员提供高质量的物理量原型数据，便于他们研建算法和应用模型。数据经过严格的系统误差订正并重采样为经纬投影，像元分辨率为1 km，像元亮度值是地物在大气顶层的反射率。除提供4个波段原始数据外，还根据用户需要提供相关辅助参数，如大气状况、系统信息（太阳的天底角、方位角、视场角、接收时间）和地形数据等。

VGT.S（Synthesis）产品是经过大气纠正的地表反射率数据，并运用多波段合成技术来获得，空间分辨率为1 km。VGT.S产品包括每天合成的4个波段的光谱反射率及NDVI数据集（S1），每10天合成的4个波段的光谱反射率及10天最大化NDVI数据集（S10）。采用最大值合成法（MVC）合成时间系列NDVI数据集，它可以最小化云和大气散射的影响。NDVI（S10）数据在发布之前，VITO所做的预处理包括大气校正、辐射校正和几何校正，并将NDVI图像上 -1到-0.1的值设置为 -0.1，再使用公式DN=(NDVI+0.1)/0.004转换到0～250的DN值。

3.2 研究方法

3.2.1 数据处理方法

本文使用的数据VGT NDVI S10数据是由10天数据采用最大值合成法（MVC）生成的，并且进行了大气校正、辐射校正和几何校正。本文在研究植被覆盖变化趋势时，数据选取策略是选取一年中植被最茂盛时期的NDVI值，但该时间点不容易确定。由于气候的影响，对于不同年份来说每年8月11日，8月21日或9月1日都有可能是一年中植被最茂盛的时间，因此选取一年中36旬的每个像元的NDVI最大值MAXNDVI作为研究数据。

其中，MAX(NDVI，i)为第i年的该像元最大化NDVI值，也是一年内该像元上植被最茂盛时期NDVI值；NDVI(i，j)为第i年第j旬的NDVI值。

3.2.2 一元线性回归变化斜率法

回归分析是研究多个变量之间统计联系的一种重要方法，是研究植被长时序变化趋势的重要方法[9]。对一组时间自变量x与NDVI因变量y数据，可以用如下的数学模型来描述：

式中，a，k是未知常数，ε是随机误差。利用观测值(xi，yi) (i=1，2，…，n)可以求出未知参数k：

其中

对于NDVI长时序数据，采用最小二乘法线性拟合后得到相应的线性方程，方程的斜率k说明像元NDVI值的多年度变化趋势（Weiss等，2001；Fuller，1998；Camberlin等，2007），k＞0，植被活动增强，k＜0，植被活动减弱。

3.2.3 Sen+Mann-Kendall趋势分析

线性回归法要求时间序列数据符合正态分布，并且易受噪声干扰。Sen趋势度是经过计算序列的中值，它可以很好地减少噪声的干扰，但其本身不能实现序列趋势显著性判断，而Mann-Kendall方法本身对序列分布无要求且对异常值不敏感，因此引入该方法可完成对序列趋势显著性检验。1998年到2011年SPOT VEGETATION NDVI时间序列数据由于受于大气和云层等因素的影响，很可能存在部分异常值，并且该序列的分布特征没有定论，所以采用上述两种方法结合可以增强方法的抗噪性，并在一定程度上提高检验结果的准确性。因此该方法在今后的研究当中应该受到足够的重视和利用[10]。

Sen趋势度计算公式为：

使用趋势度β来判断时间序列趋势的升降，当β＞0时，时间序列呈上升的趋势，反之呈下降的趋势。

Mann-Kendall趋势检验法过程如下：对于序列Xt=(x1，x2，…，xn)，先确定所有对偶值(xi，xj，j＞i)中 xi与xj的大小关系（设为S）。做如下假设：H0：序列中的数据随机排列，即无显著趋势，H1：序列存在上升或下降单调趋势。检验统计量S由公式（5）计算：

根据时间序列长度n值大小的不同，显著性检验统计量的选取有所不同：

当n＜10时，直接使用统计量S进行双边趋势检验。在给定显著性水平α下，如果|S|≥Sα/2则拒绝H0认为原序列存在显著趋势，否则接受H0认为序列趋势不显著。如果S＞0，则认为序列存在上升趋势，S=0，无趋势，S＜0认为序列存在下降趋势。

当n≥10时，统计量S近似服从标准正态分布，使用检验统计量Z进行趋势检验，Z值由公式（7）计算：

本文中时间序列长度为14（1998年-2011年），所以采用检验统计量Z来进行趋势检验，检验中取显著水平α=0.05，Z1-α/2=Z0.975=1.96。当 β＞0且|Z|＞1.96序列呈显著上升趋势，当 β＞0且|Z|≤Z1-α/2序列呈上升但不显著趋势，同理当β＜0且|Z|＞1.96时序列呈显著下降趋势，当β＜0且|Z|≤Z1-α/2序列呈下降但不显著趋势。

4 结果与分析

4.1 年度NDVI最大值的平均值变化趋势

年度NDVI最大值的平均值时间序列在一定程度上反映了北京地区的植被变化情况。从表1和图1中可以看出从1998年至2011年期间，虽然数据波动明显，但是整体上北京地区的植被变化趋势是上升的。将1998-2011时间段细分为三段：1998-2003，2004-2007，2008-2011，并对三个时间段的数据分别使用线性回归进行曲线拟合，得到表2数据，通过表2可以看出：1998-2003年期间植被变化呈下降趋势，2003-2007年植被变化呈直线上升趋势，2008-2011年植被变化呈缓慢上升趋势，这与北京申奥前与申奥成功后加强城市绿化结果是一致的。

表1 NDVI年度最大值的平均值

图1 1998-2011年NDVI年度最大值的平均值散点图

表2 分时间段线性回归结果

4.2 一元线性回归变化斜率法

利用回归法得到的北京地区植被变化趋势（如表3所示）反映了1998-2011年期间北京地区植被变化趋势是“整体改善、局部恶化”。约24%的地区植被呈退化趋势，其中严重退化区域约为8%，轻微退化区域约占16%。约76%的地区植被变化趋势上升，明显改善区域占26.33%，轻微改善区域占49.57%。图4由图2和图3叠加并添加比例尺和图例产生，该图清楚地表明了北京地区植被变化趋势的空间分布格局：植被变化增强的区域集中在北京市城区、怀柔区北部、延庆县东北部、密云县除城区附近区域、门头沟区部分地区和平谷的部分区域，并以北京市城区、怀柔区北部、延庆县东北部、密云县植被改善最为明显。这一现象说明近年来的城区绿化和郊区的退耕还林效果显著，特别是北京市城区的植被改善和北京市政府一系列绿化政策是分不开的。植被轻微退化和严重退化的区域呈现马蹄形包围着城区，大部分位于靠近北京市城区的北部、东部和南部地区，主要集中在昌平、顺义、通州和大兴四个区，上述四个区出现植被轻微甚至显著退化现象跟近几年来该地区房地产业大力发展有密切联系，人为破坏植被现象比较严重。

表3 回归分析法北京地区植被变化趋势

4.3 基于Sen+Mann-Kendall法趋势分析

表4表明北京地区植被活动整体上呈上升趋势，但是部分地区也存在着严重退化现象，这与回归分析结果相一致。7.46%的地区植被退化明显，14.25%的地区植被轻微退化，41.99%的地区植被轻微改善，36.30%的地区植被明显改善，植被变化存在上升趋势区域的面积达到78.29%，由以上数据可以看出北京地区植被变化趋势整体上是向着上升的趋势发展的。

表4 Sen+Mann-Kendall法趋势分析结果表

图2 北京市行政区划图

图3 回归分析法北京地区植被变化趋势图

图4 图2与图3叠加图

图5 Sen+Mann-Kendall法北京地区植被变化趋势图

图6 图2与图5叠加图

图6是图2与图5叠加并添加比例尺和图例生成的，从图6中可以直观地看出植被变化存在上升趋势的区域主要集中在以下地区：北京市城区、延庆县、怀柔区、密云县、门头沟和房山区，这与回归分析时的结果基本一致。植被变化呈下降趋势的地区主要集中在靠近北京市城区的北部、东部和南部，呈马蹄状包围着北京城区，只有西部地区没有合拢，这与回归分析法得到的结果一致。植被变化呈下降趋势的区域主要出现昌平区、顺义区、通州区、大兴区和房山区的东部。这也与回归分析法结果基本上一致。

分析表3可得到如下结论：采用回归分析法，当k≤0时，植被变化退化区域占总面积的24.10%，当k＞0时，植被变化上升地区占总面积的75.90%。分析表4可以得出：采用Sen+Mann-Kendall法，当β＜0时，植被退化区域为21.71%，β≥0时，植被上升区域为78.21%。综合分析上面两组数据得出两种方法在分析植被变化退化或上升区域的差异仅为2.39%，由此可知这两种方法针对植被变化趋势结果具有良好的一致性，Sen+Mann-Kendall法有很高的准确性。

5 结论

（1）1998年-2011年间，北京地区植被活动整体呈增强趋势，植被变化趋势上升的区域占北京市总面积的75%以上。北京市城区、延庆县、怀柔区和密云县的大部分地区植被活动趋势显著上升。但是，局部地区植被活动趋势显著下降，且下降区域呈马蹄状包围着北京城区，主要集中在北京市城区的北部、东部和南部，以昌平、顺义、通州和大兴靠近城区的部分退化最为严重。

（2）Sen趋势度是序列的中值，有良好的抗噪性，Mann-Kendall检验对数据的分布无要求，其对时序显著性判断有坚实的统计学理论基础，并且方法本身有一定的抗噪性，在遥感数据存在一定的噪声情况下，两种方法结合能较好的克服噪声对分析结果的影响；此外该方法分析结果与一元线性回归法在植被变化下降（上升）区域的差异仅为2.39%，说明该方法有较高的准确性。该方法与GIS技术结合可以直观呈现北京地区植被变化趋势的时空分布格局，并且生成的植被变化趋势图与回归法有良好的空间一致性。基于以上优点，该方法可以广泛应用到其他区域植被变化趋势分析中。

[1]陈云浩，李晓兵，史培军.基于遥感的植被覆盖变化景观分析[J].生态学报，2002，22（10）：1581-1586.

[2]李晓琴，孙丹峰，张凤荣.基于遥感的北京山区植被覆盖景观格局动态分析[J].山地学报，2003，21（3）：272-280.

[3]张学霞，葛全胜，郑景云.近50年北京植被对全球变暖的响应及其时效—基于遥感数据和物候资料的分析[J].生态学杂志，2005，24（2）：123-130.

[4]张本昀，喻铮铮，刘良云，等.北京山区植被覆盖动态变化遥感监测研究[J].地域研究与开发，2008，27（1）：108-112.

[5]郭锦，张晓丽，赵丽琼，等.基于3S技术的北京市植被覆盖空间格局变化研究[J].安徽农业科学，2009，37（17）：8264-8266.

[6]祁燕，王秀兰，冯仲科，等.基于RS与GIS的北京市植被覆盖度变化研究[J].林业调查规划，2009，34（2）：1-4.

[7]张萍，彭道黎，李万德，等.延庆县植被覆盖动态变化监测[J].湖北农业科学，2009，48（5）：1132-1136.

[8]霍亚贞.北京自然地理[M].北京：北京师范学院出版社，1989.

[9]韩秀珍，李三妹，罗敬宁，等.近20年中国植被时空变化研究[J].干旱区研究，2008，25（6）：753-759.

[10]蔡博峰，于嵘.基于遥感的植被长时序趋势特征研究进展评价[J].遥感学报，2009，13（6）：1177-1186.

[11]朴世龙，方精云.最近18年来中国植被覆盖的动态变化[J].第四纪研究，2001，21（4）：294-302.

[12]方精云，朴世龙，贺金生，等.近20年来中国植被活动在增强[J].中国科学：C辑，2003，33（6）：554-565.

[13]邱海军，曹明明.基于SPOT VEGETATION数据的中国植被覆盖时空变化分析[J].资源科学，2011，33（2）：335-340.

[14]Cleveland R B，Cleveland W S，McRae J E，et al.STL：a seasonal-trend decomposition procedure based on loess[J]. Journal of Official Statistics，1990，6（1）：3-73.

[15]Shabanov N V，Zhou Liming，Knyazikhin Y，et al.Analysis of interannual changes in northern vegetation activity observed in AVHRR data from 1981 to 1994[J].IEEE Transactions on Geoscience and Remote Sensing，2002，40（1）：115-130.

[16]Stow D A，Hope A，McGuire D，et al.Remote sensing of vegetation and land-cover change in arctic tundra ecosystems[J]. Remote Sensing of Environment，2004，89：281-308.

[17]Herrmann S M，Anyamba A，Tucker C J.Recent trends in vegetation dynamics in the African sahel and their relationship to climate[J].Global Environmental Change，2005，15：394-404.

[18]Zhou Liming，Tucker C J，Kaufmann R K，et al.Variations in northern vegetation activity inferred from satellite data ofvegetation index during 1981 to 1999[J].Journalof Geophysical Research，2011，106：20069-20083.

WANG Dianlai1，2,LIU Wenping1,HUANG Xinyuan1

1.College of Information,Beijing Forestry University,Beijing 100083,China
2.Department of Computer Science,Shougang Institute of Technology,Beijing 100144,China

The spatial distribution and dynamic change of vegetation cover in Beijing are analyzed based on SPOT VEGETATION's NDVI data from 1998 to 2011 using the method of slope of linear regression and Sen+Mann-Kendallanalysis.Experimental results show that the vegetation's change is significantly increased in the following areas：the urban of Beijing city,Yanqing county, Huairou and Pinggu district.On the contrary,the vegetation's decreasing areas locate at the north,east and south of the urban district in Beijing and surround like a Horseshoe-shaped.The two methods'sresults have good spatial consistency in increasing （or decreasing）vegetation coverregions.Sen+Mann-Kendall analysis can be widely applied in other areas for detecting the trends of vegetation's change thanks to its noise immunity and without requirements of data distribution.

SPOT VEGETATION;trend analysis in vegetation cover;linear regression;Sen's slope estimator;Mann-Kendall test

基于1998年到2011年长时序SPOT VEGETATION归一化植被指数数据，采用一元线性回归斜率变化法和Sen+ Mann-Kendall法对北京地区的植被变化趋势做了时空分析。实验结果表明：在1998年到2011年期间，北京市城区、延庆县、怀柔区和平谷区的植被变化趋势显著上升；而植被恶化区则集中在北京市城区北部、东部和南部，并以马蹄形包围北京市区。两种方法实验结果在植被上升（下降）区域具有一致性。同时，Sen+Mann-Kendall法以其良好的抗噪性和对数据分布无要求性可广泛应用到其他区域的植被变化趋势分析中。

SPOT VEGETATION；植被变化趋势分析；一元线性回归；Sen趋势度估计法；Mann-Kendall检验

A

TP399

10.3778/j.issn.1002-8331.1206-0282

WANG Dianlai,LIU Wenping,HUANG Xinyuan.Trend analysis in vegetation cover in Beijing based on Sen+Mann-Kendall method.Computer Engineering and Applications,2013,49（5）：13-17.

国家重点基础研究发展规划（973）（No.2009CB421105）；中央高校基本科研业务费专项资金资助（No.YX2011-28）。

王佃来（1972—），男，博士，高级工程师，研究领域为数字图像处理、软件设计与开发；刘文萍（1971—），通讯作者，女，博士，教授，研究领域为数字图像处理及模式识别；黄心渊（1965—），男，博士，教授，博导，研究领域为计算机动画、计算机可视化等。

E-mail：wangdl12345@126.com

2012-06-18

2012-09-20

1002-8331（2013）05-0013-05

CNKI出版日期：2012-11-06 http://www.cnki.net/kcms/detail/11.2127.TP.20121106.1607.003.html