基于SPOT-VEGETATION數據的神農架林區1998—2013年植被覆蓋度格局變化

劉家琰,謝宗強,申國珍,樊大勇,熊高明,趙常明,周友兵,徐文婷,*

1 中國科學院植物研究所植被與環境變化國家重點實驗室,北京　100093 2 中國科學院大學,北京　100049

植被是覆蓋地表的森林、灌叢、草地和農作物等的總稱,在陸地生態系統的碳水循環和能量交換方面有著重要的意義[1]。植被覆蓋度(fractional vegetation cover, FVC)是指單位面積內植被(包括葉、莖、枝)的垂直投影所占百分比[2],是描述植被群落及生態系統的重要參數[3- 6]。

全球氣候變化、人類活動都會對植被覆蓋產生重要的影響,使植被-土壤系統的蒸散量、植被表面的反射率和粗糙度等特性發生改變,從而進一步影響生態系統的能量平衡、碳循環過程以及生產力[7]。對植被覆蓋度進行不同尺度的監測十分必要。例如,對于干旱半干旱地區由于畜牧業和沙漠化造成的土地退化來說,植被覆蓋度是評估該區草地狀況的重要指標[8]。植被覆蓋度及其變化對水文、生態、全球變化等都具有重要意義,在植被變化[9]、生態環境調查、水土保持研究、蒸散量研究以及土地退化[10]、鹽漬化[11]和沙漠化[12]等研究領域都有廣泛的應用。

植被覆蓋度的監測方法可分為地面測量和遙感監測[13],而遙感由于其大范圍的數據獲取和連續觀測,已成為植被覆蓋度研究的主要手段[14]。GIMMS (Global Inventory Modeling and Mapping Studies)、MODIS (Moderate Resolution Imaging Spectroradiometer)、SPOT-VEGETATION等數據產品為研究長時間序列的植被覆蓋變化提供了可能。研究表明,通過建立遙感數據的植被指數或多波段光譜組合與地面實測植被覆蓋度間的回歸關系,可以較好的估算小范圍的植被覆蓋度。Shoshany等[15]利用TM數據建立多元線性回歸模型估算植被覆蓋度,相關系數達到0.88；Xiao等[16]利用ETM NDVI數據的線性回歸方程對美國新墨西哥州中部區域的植被覆蓋度進行擬合,相關系數達0.89。Graetz等[17]利用MSS數據與實測數據估算了澳大利亞南部半干旱地區稀疏草地的植被覆蓋度,有效地監測到了牧區植被覆蓋度在經受干旱、放牧脅迫前后的變化,為該牧區植被退化程度的監測提供了技術支撐；Boyd等[18]利用非線性回歸模型估算了美國西北部的針葉林蓋度,相關系數為0.56。像元分解模型可成功用于估算區域尺度的植被覆蓋度。孫久虎等[19]利用像元二分模型估算了北運河地區1994—2004年間植被覆蓋變化情況,結果表明研究區內植被退化嚴重,植被覆蓋度下降了9.99%；陳云浩等[13]利用植被覆蓋度計算模型對北京市海淀區的植被覆蓋進行動態遙感監測研究,結果顯示海淀區1975—1997年間西北部山區植被覆蓋度呈持續增長趨勢,而植被覆蓋度為25%—50%的區域則表現為由東向西的大面積減少,一定程度上反映了城市化過程對植被覆蓋的影響；李苗苗等[20]通過NDVI像元二分法估算密云水庫上游植被覆蓋度,溝谷、平地等地勢較低的區域植被覆蓋度受人為干擾較大,而地勢較高的山區植被覆蓋度變化較小,且經過驗證該結果具有較高精度(相關系數0.89,協方差0.02)；張喜旺等[21]對伊洛河流域的植被覆蓋度空間分異進行了研究,Detsch等[22]利用GIMM數據對植被覆蓋動態及季節性進行了研究。

神農架被認為是我國華中地區的“生物基因庫”,具有豐富的生物多樣性和森林資源,是我國保存較為完整的亞熱帶亞高山森林生態系統區域[23-24]和華中地區唯一現存的原始林分布區。19世紀70、80年代,由于人口增長、土地利用/土地覆蓋劇烈變化,神農架林區覆蓋度大幅度下降[25-26]。1986年,神農架林區建立國家級森林和野生動物自然保護區,2000年全面停止對天然林的砍伐,實行“天保工程”和退耕還林工程。目前有關神農架林區植被覆蓋變化的研究較少,其中黃靖等[27]利用2003—2012年的MODIS EVI (enhanced vegetation index)數據分析了神農架林區植被動態變化及其與氣候因子的關系,發現這10年來林區植被覆蓋整體呈增加趨勢且東部增幅大于西部,氣溫是影響植被覆蓋的主要因子；姜哲等[28]利用Landsat數據對神農架林區1987年、2000年和2013年3個時期的土地覆蓋進行解釋,發現兩個時段內森林凈增長面積分別為14.70 km2和207.49 km2。但以上研究并未指出神農架林區植被覆蓋的變化程度, 且較少涉及其與社會經濟因素的關系。

本文利用1998—2013年的SPOT-VEGETATION歸一化植被指數數據,應用像元二分法估算了神農架林區的植被覆蓋度及其變化趨勢,分析了神農架國家級自然保護區內植被覆蓋度變化與保護區外的差異,為進一步有效地保護神農架林區原始天然林提供依據。

1　材料與方法

1.1　研究區概況

神農架林區(109°56′—110°58′ E, 31°15′—31°75′ N)位于鄂西邊陲,是全國唯一以林區命名的縣級行政單位,總面積3250 km2,擁有全球中緯度地區唯一一塊保存完好的原始林,森林覆蓋度為88.7%[28]；神農架國家級自然保護區位于神農架林區的西南部,面積為769.5 km2。神農架林區為秦嶺山系大巴山脈東段,地勢由西向東,由南向北逐漸降低,海拔3105.4 m的最高峰神農頂是華中第一峰。該區地處北亞熱帶季風區,氣溫偏低且多雨,年平均氣溫12℃,年平均降水量1185 mm[27]。區內植被以亞熱帶成分為主,兼有溫帶和熱帶成分,并具有明顯的垂直地帶性,海拔400—1000 m為亞熱帶常綠闊葉林帶,海拔1000—1700 m為北亞熱帶常綠落葉闊葉混交林帶,海拔1700—2200 m為暖溫帶落葉闊葉林帶,海拔2200—2600 m為溫帶針闊混交林帶,海拔2600—3000 m為寒溫性針葉林帶,3000 m以上為亞高山灌叢草甸帶。

1.2　數據源及處理

1.2.1遙感數據

下載 (http://www.vgt.vito.be/) 1998—2013年法國SPOT-VEGETATION衛星的的歸一化植被指數數據,數據空間分辨率為1 km。該數據為每10 d最大化合成(maximum value composite, MVC)的NDVI數據,將云和大氣散射的影響降到最低[29]。將數據換算,轉化為NDVI原始值。

1.2.2氣象數據

利用中國氣象數據網(http://data.cma.cn/)上的中國地面國際交換站氣候資料日值數據集(V 3.0),得到1998—2013年間中國各氣象臺站的氣溫、降水量數據,在ArcGIS 10.1中利用Kriging插值法進行插值,得到神農架林區1 km分辨率的平均氣溫、最低氣溫和降水量分布圖。

1.2.3DEM數據

從CGIAR CSI(http://www.cgiar-csi.org/)網站下載研究區域90 m分辨率的STRM DEM數據。

1.2.4物種數數據

根據Zhao等[30]發表的論文,提取神農架林區不同海拔段的群落物種數和喬木物種數據。該研究中,海拔每升高100 m設置一個樣方進行調查,或選取特殊的群落類型進行樣方調查,共設喬木樣方166個。對喬木進行每木調查,并記錄胸徑、樹高、冠幅等信息,以及樣地的海拔、坡度、坡向、經緯度等環境因素。

1.3　研究方法

1.3.1植被覆蓋度的計算

遙感圖像中的一個像元往往包含多個組分(如植被、裸土),各組分均對遙感器所觀測到的像元信息有所貢獻；我們常用像元分解模型將像元中的遙感信息(如波段、植被指數)分解,并進一步估算植被覆蓋度。像元二分模型是最簡單的一種像元分解模型,其只由植被覆蓋地表、無植被覆蓋地表兩部分構成。研究表明,NDVI對土壤背景變化敏感,且其大小取決于植被覆蓋度和葉面積指數,因此可以用NDVI反演植被覆蓋度,計算公式如下[31]：

(1)

式中,NDVIsoil表示裸地的歸一化植被指數,NDVIveg表示完全植被覆蓋的植被指數。

植被覆蓋度小于15%時,NDVI可將土壤與植被較好區分開；植被覆蓋度在25%—80%范圍內,NDVI隨植被覆蓋度線性增加；植被覆蓋度大于80%時,其檢測能力有所下降。盡管NDVI存在自身的局限性且許多其他的植被指數因考慮了土壤、大氣等因素得到發展, 但NDVI因其對植被檢測的高靈敏度、對植被覆蓋度較寬的檢測范圍、對地形及群落結構造成陰影干擾的消除和對太陽高度角及大氣帶來噪聲的削弱等特點被廣泛應用[32]。Leprieur等[33]比較了不同植被指數估算植被覆蓋度的能力,結果表明盡管全球環境監測指數(global environmental monitoring index, GEMI)考慮大氣的影響,修正型土壤植被指數(modified soil-adjusted vegetation index, MSAVI)考慮土壤背景的影響,但在實際應用中歸一化植被指數NDVI更具優勢。

求解公式(1)的關鍵在于NDVIsoil和NDVIveg的取值,而兩者的值在不同圖像中應是不一樣的,因此不可取定值,而應由圖像本身決定,以削弱干擾因素的影響[34]。通常NDVIsoil是NDVI圖像中的最小值,而NDVIveg是圖像中的最大值,Gutman和Ignatoo[5]提出以NDVI∞替代NDVIveg,以NDVI最小值替代NDVIsoil,然而該基于像元分解的密度分解的方法在實際中難以應用。李苗苗[34]將該模型進一步簡化,為排除異常值,取NDVI圖像累積頻率為95%、5%處的像元值分別作為NDVIveg和NDVIsoil,也就得到了以下公式：

(2)

該研究利用以上公式計算得到密云水庫上游植被覆蓋度,經精度驗證,實測值與估計值相關系數達99%,協方差0.06,表明該模型估計植被覆蓋度具有很高的準確性。李曉錦[35]利用該模型估算黃土高原植被覆蓋度,基于Landsant/TM數據提取的植被覆蓋度精度為95.02%, 同樣具有較高的精度。

圖1　NDVI在置信區間上的取值Fig.1　NDVI value in confidential interval

根據公式(2)計算時間序列上每一景(每一時相)影像中NDVI的頻率累積值, 取頻率為5%處的NDVI值作為NDVIsoil,取95%處的NDVI值作為NDVIveg(圖1),得到1998—2013年每旬的逐像元植被覆蓋度,按旬平均得到全年平均植被覆蓋度(下文中如無特殊說明,簡稱植被覆蓋度)。為進一步分析研究區內植被覆蓋度在時間、空間上的分異與變化,在所得逐像元植被覆蓋度的基礎上,計算多年平均植被覆蓋度(逐像元及研究區平均)、年平均最大植被覆蓋度(即年最大植被覆蓋度的多年平均值,同樣按逐像元及研究區內平均計算)。

1.3.2植被覆蓋度變化率的計算

為了更好地分析1998—2013年間神農架林區植被覆蓋度的變化程度以及變化在空間上的分布,需消除年季節動態導致的植被覆蓋度波動影響,可采用如下公式[36]計算每年生長季(5—9月) FVC的變化率：

FVCα= 直線斜率 / 均值 × 16 × 100

(3)

式中,FVCα為植被覆蓋度變化率,直線斜率為1998—2013年生長季平均FVC對年份回歸的直線斜率,均值為多年生長季平均的FVC值。所有空間分析及空間動態分析均在ArcGIS 10.0中進行。

1.3.3植被覆蓋度動態變化的影響因素分析

利用R語言[37]分析神農架林區植被覆蓋度的動態變化與距離居民地最近距離、距離公路最近距離、平均溫度、降水量、年最低溫度、海拔等因素的相關性。

1.3.4植被覆蓋度海拔梯度格局與物種數間關系分析

利用R語言分析神農架林區植被覆蓋度的動態變化與海拔梯度、不同海拔梯度上的物種數的相關關系。

2　結果

2.1　神農架林區平均植被覆蓋度

神農架林區1998—2013年間的多年平均植被覆蓋度為66.8%,年際間變動范圍64.3%—71.4%,其中2012年植被覆蓋度最低(圖2)。年平均最大植被覆蓋度為93.8%,年際間變動范圍92.5%—94.5%,2005年最大植被覆蓋度最低。保護區內和保護區外的多年平均覆蓋度分別為66.5%和66.8%；保護區內和保護區外的多年平均最大植被覆蓋度分別為94.3%和93.7%,保護區內植被覆蓋度顯著高于保護區外(成對t檢驗,t= 3.221,P< 0.01)。

圖2　神農架國家級保護區內外平均植被覆蓋度及最大覆蓋度(平均值±標準差)Fig.2　Mean and maximum fractional vegetation cover inside and outside the Shennongjia National Nature Reserve (average ± standard deviation)

1998—2013年間神農架林區多年平均覆蓋度在空間分布上存在差異(圖3),空間上最大植被覆蓋度為82.4%,最小為27.7%。空間聚類分析發現,高覆蓋度聚集區主要分布在陰峪河流域、九沖河流域、里叉河流域、溫水河流域以及宋洛河流域,而低覆蓋度聚集區主要在大九湖、松柏、老君山以及保護區內旅游公路兩側區域。

圖3　神農架林區1998—2013年平均植被覆蓋度空間分布Fig.3　The mean fractional vegetation cover pattern of Shennongjia Forest District

植被覆蓋度海拔格局分析結果表明,海拔1500—2000 m區域植被覆蓋度最高(R2= 0.185,P< 0.01)；植被覆蓋度在海拔1700 m以下隨著海拔的升高而增加,而在1700 m以上,隨著海拔的升高而降低(圖4)。根據Zhao等[30]有關喬木物種數和群落總物種數的調查數據,結合本研究估算的植被覆蓋度,覆蓋度與喬木物種數和總物種數隨海拔的變化趨勢相似(圖4),且覆蓋度與喬木物種數(r= 0.953,P< 0.01)、總物種數(r= 0.807,P< 0.01)之間有著極高的相關性；覆蓋度和喬木的物種豐富度均在海拔1500—2000 m之間最高,而這一海拔段也是典型地帶性植被常綠落葉闊葉混交林分布區域。

圖4　神農架林區植被覆蓋度海拔梯度格局及與物種數間的關系Fig.4　The mean fractional vegetation cover altitude pattern of Shennongjia Forest District and relationship with species richness

2.2　1998—2013年間神農架林區植被覆蓋的變化

神農架林區在1998—2013年間植被覆蓋度的平均變化率為1.45%,其中保護區內變化率為2.26%,保護區外的變化率為1.23%。單因素方差分析方法(one-way ANOVA)結果表明保護區內外植被覆蓋度的平均變化率有顯著差異(P< 0.01)。老君山頂、神農頂、青樹村、東溪村以及保護區西南部的相思嶺村等地植被蓋度均有明顯增加,保護區旅游公路兩側植被覆蓋度也有少量增加；植被覆蓋率降低的區域主要集中在大草坪和宋洛村(圖5)。

2.3　神農架地區植被覆蓋動態變化的影響因素分析

對植被覆蓋度變化率及其影響因素進行相關性分析,結果表明林區1998—2013年間植被覆蓋度的變化率受溫度、降水、年最低溫、距居民地距離、距公路距離等因素影響(溫度:r= 0.12,P< 0.01；降雨:r= 0.18,P< 0.01； 年最低溫:r= 0.04,P< 0.01；距居民地距離:r= -0.09,P< 0.01；距公路距離:r= -0.16,P< 0.01)。海拔與植被覆蓋度變化率無顯著相關性(表1)。

3　討論

神農架林區是全球中緯度地帶唯一現存的原始林分布區,保存有較好的自然植被[23]。其植被覆蓋度在空間上的分布存在明顯的差異。原始林主要分布在陰峪河流域、金猴嶺,溫水河流域以及宋洛河流域也分布有少量原始森林。因而這些區域植被覆蓋度較高。林區西部大九湖區域具有大量農田、濕地和少量裸巖,導致其植被覆蓋度較低；老君山以及保護區內旅游公路的西北段,由于海拔較高,植被類型主要是高山草甸和灌叢,因而植被覆蓋度較低。

植被覆蓋度多年平均變化率表明,1998—2013年神農架林區植被覆蓋在空間上有增有減, 且差異明顯,但整體變化率為1.45%,植被覆蓋呈現較好的增長趨勢,這與前人研究的結果一致,如姜哲等[28]利用Landsat TM數據對3期土地覆蓋進行分類,發現神農架林區森林面積凈增加207.49 km2；黃靖和夏智宏[27]利用Modis EVI數據研究2003—2012年神農架林區植被指數的變化,也發現植被指數呈現增長的趨勢。

表1　植被覆蓋度變化率與各因素相關性

*P< 0.05;**P< 0.01

保護區內植被覆蓋度空間分布于1998—2013年間有增有減,但總體保護效果優于保護區外。保護區內植被覆蓋度的增加幅度約為保護區外的2倍,表明國家級自然保護區發揮了良好的保護效果。除植被自然恢復外,2000年移栽人工日本落葉松林(Larixkaempferi(Lamb.) Carr.),使保護區東部(老君山頂)植被覆蓋度顯著增加；保護區北部(青樹村、東溪村)于2001年實施移民搬遷,植被人為干擾較少,處于恢復狀態,因而植被覆蓋度有所增加。保護區內同樣存在植被覆蓋度下降的區域,集中分布在陰峪河流域、落陽河流域、黃柏阡和興隆寺村、坪阡村等區域。保護區內植被覆蓋度降低,可能是由于受到2004年底、2005年初極端低溫的影響(其間最低溫度達-13.2℃),以及2008年初極端低溫(低達-17.7℃)的影響,尤其部分地區(陰峪河流域、落陽河流域等)因海拔較低,常綠植物占優勢的植被遭受到的冰凍雪災更為嚴重；除此之外,保護區內部分地區(如落陽河流域、黃柏阡、興隆寺村等)均有相當面積的華山松分布(1300—2700 m)[38-39],而神農架林區華山松大小蠹自2012年以來呈爆發趨勢,華山松林受災嚴重[40],導致其覆蓋度大幅度下降。再則,2012年保護區內水庫的修建(坪阡村附近),也導致植被覆蓋度有所下降。

保護區外植被覆蓋度亦同時存在增加和下降的趨勢。神農架機場和滑雪場等設施的修建導致保護區外部分地區(大草坪、宋洛村附近)植被減少,植被覆蓋度下降。而保護區外的公路對植被覆蓋度的影響較為復雜。林區內的兩條主要公路為307省道和209國道,307省道兩側植被覆蓋度有所減少,而209國道兩側植被覆蓋度呈明顯增加趨勢；其主要原因是209國道由于鄂西生態文化旅游圈項目的實施,于2008年完成道路兩側共15.8 km2的綠化工作,這與Pauleit等[41]關于城鎮化表現為人工綠地的增加和自然植被的減少的研究結果一致。

Ge等[42]對神農架長期監測樣地的分析表明,極端低溫事件極有可能威脅亞熱帶地區以常綠樹種為主的闊葉林的存在。這與本文神農架保護區內中低海拔以常綠物種為優勢的區域植被覆蓋度的降低較為明顯的研究結果一致。此外,人類干擾、自然災害等均會對植物生長產生影響,進而影響植被覆蓋度。因此,未來對保護區的保護和管理不僅僅應關注人類干擾,對自然干擾也應予以重視。