基于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