呈貢區近30年植被覆蓋度時空變化特征與土地利用驅動

楊旭超, 張 軍, 李 杰, 勞潔英, 吳志娟

(1.云南大學 資源環境與地球科學學院, 昆明 650504;2.昆明市土地開發整理中心, 昆明 650504; 3.云南省地礦測繪研究院, 昆明 650504)

當前，受人類活動以及氣候變暖等因素影響，生態退化已發展成為全球環境變化的主要問題，受到各國政府及科學部門高度關注[1]。植被是陸地生態系統重要組成部分，可將土壤、大氣和水分等自然過程進行連接，并在物質循環、能量流動和信息傳遞等方面均發揮著重要作用[2]。同時，其還具有截留降雨、保土固土等功能，是土壤侵蝕與水土流失的主要監測因子，對于區域生態退化具有實時客觀的指示效果[3]。植被狀況的好壞，主要通過植被覆蓋度(FC)來表示。FC作為區域氣候模型、水土流失監測、土地沙漠化評價和分布式水文模型的重要因子，是反映地表植被分布特征和描述生態系統的重要參數，具有明顯的季候和年際變化特征，通過監測其空間分異及演化趨勢將對區域生態環境變化具有重要指示意義。

傳統的植被覆蓋度研究方法受測量條件限制，主要以野外實地測量為主，包括目估法、采樣法、儀器測量法及模型法等。而大區域的植被覆蓋監測更需建立在多光譜、多時相和多尺度的遙感數據基礎上。近年來，隨著對地觀測系統的發展，國內外諸多學者基于遙感手段，在植被監測領域進行了大量研究，探索出許多監測方法。如羅亞等[4]對模型法、植被指數法和混合像元分解模型法進行了具體研究。像元二分模型是一種簡單的混合像元分解模型，其不需依賴實測數據，模型計算簡單可靠、輸入參數通用易得，被廣泛應用于區域尺度上植被覆蓋動態變化及空間格局特征等研究中。如肖強[5]、廖清飛[6]、李杰[7]、李鈺溦[8]等學者分別以黃土高原、青海省東部、云南怒江、中國北方等地區為研究區，對其植被覆蓋時空演變遙感監測與分析。

縱觀國內外，目前在植被覆蓋度領域的研究多是采用大面積、低分辨率的遙感時序數據，如黃森旺等[9]基于GIMMS AVHRR-NDVI數據和氣象數據分析了三北防護林工程區25年的土地退化情況，Piao等[10]基于SPOT VGT-NDVI和MODIS-NDVI數據分析了歐亞大陸溫帶和北極地區植被生長趨勢的變化；對于中高分辨率的衛星數據的研究主要以Landsat系列遙感數據，但研究較少。上述研究為研究區域今后土地開發利用、環境保護治理等提供有用信息，但其均采用典型年份或等間距年份的數據，并未采用逐年的時序數據，這使得研究內容存在一定局限性[11]。

在云南省委、省政府建設“現代新昆明”的戰略背景下，昆明市呈貢縣自2011年7月8日撤縣設區，更名為呈貢新區。近年來，該區搶抓世界經濟形勢好轉和“一帶一路”建設新機遇，生態建設牢固樹立“綠水青山就是金山銀山”的發展理念，嚴守生態劃定紅線，確保森林覆蓋率達39%以上。因此，關注呈貢區植被時空演化特征、評估未來發展趨勢，對區域生態建設具有重要的實踐意義。賀晉云等[12]研究表明，近50多年我國西南地區的干旱災害發展具有面積增大和頻率增大的趨勢，研究區所處的滇東高原腹地正好位于該范圍內，氣候變化導致極端天氣和氣象災害的頻發，勢必會對呈貢區生態環境產生嚴重影響。鑒于此，本文以呈貢區1990—2018年Landsat TM/ETM/OLI等數據為研究對象，采用小波分析、趨勢分析、變異系數、灰度預測等方法，探討研究區植被覆蓋的時空格局、演化規律；此外，結合Sentinel-2A數據，進行土地覆蓋信息提取，以獲取研究區首尾時段的土地覆蓋類型空間分布圖，構建2000—2017年研究區土地利用轉移矩陣，在植被覆蓋變化大格局下對土地轉移以及植被變化信息進行分析。通過此以期為區域可持續發展、生態修復工程規劃以及態環境保護提供科學依據。

1 研究區概況

作為“現代新昆明”建設示范區，呈貢區于2011年7月8日撤縣設區，更名為“呈貢新區”，隸屬云南省昆明市，同時也是昆明市政府駐地[13]。呈貢區位于滇東高原腹地，滇池東岸，位于102°45′—103°00′E，24°42′—25°00′N，土地總面積為461 km2，下轄10個街道辦事處，65個社區居委會。地勢東高西低，海拔為1 743～2 768 m(圖1)。屬低緯度高原季風氣候，年平均日照時數2 200 h，年均溫度14.7℃，年平均降水量789.6 mm。

圖1 呈貢區位置示意圖

2 數據與研究方法

2.1 數據源與預處理

基于本研究時間尺度為29 a(1990—2018年)，影像數據采用美國航空航天局(NASA)提供的1990—2018年逐年的Landsat數據(表1)，空間分辨率為30 m，平均云量均低于5%，數據質量良好，根據數據的可獲得性，盡量挑選3—5月為研究時相，該時段為植被的主要生長季，使其具有可比性，以滿足研究需求；此外，鑒于Sentinel-2A數據在空間分辨率及光譜信息等方面均優于Landsat數據，故選用2018年Sentinel-2A數據(數據來源:https:∥scihub.copernicus.eu/dhus/#/home)進行研究區土地利用信息的提取。地形數據選用NASA(http:∥reverb.echo.nasa.gov/reverb/)提供的GDEMV2 30 m DEM數據，探究研究區植被覆蓋空間分布的高程差異。為提高研究的可靠性，根據2018年3月的實地考察，并結合Google Earth 1990年12月歷史數據，實現了2018年和1990年地類樣本點的數據采集。地類樣本點作為研究區多年FVC均值分級，以確定植被區及對應年份影像分類精度驗證的依據。

預處理方面，因Landsat及Sentinel-2A數據均經過幾何精校正，故只需再對其進行輻射定標及大氣校正，以消除大氣散射、吸收、反射引起的誤差，得到地表反射率即可。其中，對于Landsat數據，本研究采用基于MODTRAN4+輻射傳輸模型的FLAASH大氣校正方法，該法可有效去除水蒸氣/氣溶膠散射效應；而Sentinel-2A則采用歐洲航天局的Snap軟件進行輻射定標、大氣校正，其中大氣校正在Sen2Cor-2.4.0模塊下進行。最后依據呈貢區行政區劃矢量數據對預處理后的影像進行鑲嵌和裁剪。

表1 研究區遙感影像獲取情況

2.2 研究方法

2.2.1 植被覆蓋度估測

(1) 歸一化植被指數(NDVI)。歸一化植被指數是植被生長狀態及植被覆蓋度最佳指示因子[14]，其計算公式為：

NDVI=(NIR-R)/(NIR+R)

(1)

式中：NIR為近紅外波段的反射率；R為紅光波段的反射率。

(2) 基于像元二分模型估算植被覆蓋度。像元二分模型是將一個像元的地表組成,分為有植被覆蓋、無植被覆蓋兩部分，其中植被覆蓋度可以看作是植被的權重[15-16]。因此可以用NDVI計算植被覆蓋度：

FC=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(2)

式中：；NDVIsoil為完全由裸土覆蓋像元的NDVI值；NDVIveg則代表完全由植被覆蓋像元的NDVI值。本文根據各年影像數據，構建土地利用類型ROI，確定NDVIsoil和NDVIveg閾值，并提取森林類型內NDVI概率分布的95%下側分位數所對應的NDVI值為NDVIveg參數；NDVIsoil參數則統一采用裸地覆蓋地類內5%下側分位數。研究區FC時間序列平均值需進行重分類后，將植被覆蓋度劃分為4個等級；0～0.4為低植被覆蓋區；0.4～0.6為中等植被覆蓋區；0.6～0.8為中高植被覆蓋區；0.8～1為高植被覆蓋區。

2.2.2 小波分析 小波分析能清晰揭示出隱藏在時間序列中的多種變化周期及在不同時間尺度中的變化趨勢，兼具時域和頻域多分辨功能，能對系統未來發展趨勢進行定性估計。計算公式如下：

Morlet母小波形式為：

λ=π-1/4eicte-t2/2

(3)

子小波形式為：

(4)

離散小波變換形式為：

(5)

式中：t為時間；c為無量綱頻率；a為尺度因子，表示小波周期長度；b為時間因子，反映時間上的平移；Wf(a，b)為小波變換系數，可同時反映時間參數b和頻域參數a的特性；*表示共軛復數。小波方差FC時間序列的波動能量隨尺度(a)的分布情況確定演化過程中存在的主周期。計算公式為：

(6)

2.2.3 趨勢分析法 一元線性回歸趨勢分析能夠模擬出影像中每個柵格的變化趨勢，進而反映研究區FC空間變化趨勢特征，計算過程如下：

(7)

式中：Slope為1990—2018年回歸趨勢斜率；n為研究期總年數；i為年份；MFCi為第i年FC最大值。

2.2.4 變異系數法 本文通過變異系數(Cv)反映FC時間序列的波動情況，計算方法如下：

(8)

2.2.5 面向對象的最近鄰法 面向對象的最近鄰分類法根據分割結果選取訓練樣本對象及特征空間，通過計算待分類影像對象和各類別訓練樣本之間的距離，將待分影像對象歸類至距離最近的樣本對象所屬類別中。其表達式如下：

(9)

以Landsat4 TM(1990年)以及Sentinel-2A(2018年)兩期時相相近的遙感影像為土地利用信息提取數據源，基于eCognition軟件，采用面向對象的最鄰近分類法獲取研究區兩個時段的土地覆蓋類型空間分布圖，并基于樣本點構建混淆矩陣進行精度驗證，確保分類結果的可靠性。構建2000—2017年研究區土地利用轉移矩陣，在植被覆蓋變化大格局下對土地轉移以及植被變化信息進行分析。

3 結果與分析

3.1 時間變化特征

Landsat4 TM(1990年)以及Sentinel-2A(2018年)兩期時相相近的遙感影像為土地利用信息提取數據源，基于eCognition軟件，采用面向對象的最鄰近分類法進行土地利用信息提取，以獲取研究區兩個時段的土地覆蓋類型空間分布圖，并基于樣本點構建混淆矩陣進行精度驗證，確保分類結果的可靠性。通過此，構建2000—2017年研究區土地利用轉移矩陣，以在植被覆蓋變化大格局下對土地轉移以及植被變化信息進行分析。

如圖2所示，1990—2018年呈貢區植被覆蓋整體呈緩慢增加趨勢，NDVI值在1990—2018年由0.227 6增加到0.254 1，FC由0.624 4到0.660 3，NDVI及FC的線性傾向性分別為0.03%/10 a，0.07%/10 a，兩者波動變化具一定相似性，其中，最大值均出現在2016年，NDVI及FC分別為0.377 4，0.726 1，最小值出現2001年，分別為0.583 2，0.152 9。值得注意的是，研究區植被覆蓋在2004—2009年出現持續性下降趨勢，線性傾向性分別為-0.11/10 a(FC)，-0.06/10 a(FC)，其原因可能與近年來我國西南地區長時間尺度干旱事件有關；而自2015年至今，NDVI及FC的增速分別為0.28/10 a，0.10/10 a，植被覆蓋增長趨勢較為明顯，表明呈貢新區在森林保護、造林綠化以及石漠化、難造林地荒山荒坡生態修復等工程均取得一定成效。

根據時間因子、尺度因子構建小波方差圖，由圖3可知，呈貢區的FC主要存在高頻14 a變化周期，該周期也是研究區植被變化的主要周期，主要發生在1995—2010年；此外，小波方差圖中還存在一個較為明顯的峰值，其時間尺度為7 a，對應著第二峰值。主要由以上兩個周期的波動控制著植被覆蓋在整個時間域內的變化特征。

3.2 空間變化特征

3.2.1 空間分布特征 根據樣本點統計，1990—2018年，FC值分布具一定層次性，并隨時序呈波動變化。如圖4所示，統計得低植被覆蓋區面積為53.15 km2，占研究區總面積11.53%，大部分分布于呈貢區西部，東部少量分布，因其地貌類型為水域，故植被分布極少；FC值介于0.4～0.6的區域面積為76.3 km2，占16.55%，對比1990年、2018年土地利用圖，發現該區域主要為裸地及建設用地類型，植被覆蓋較少；FC值介于0.6～0.8的區域占比最多，為面積251.08 km2的中高植被覆蓋區，占比達54.46%，主要分布于呈貢區中部大部分區域，地類以耕地、裸地及建設用地為主；高植被覆蓋區域面積為80.47 km2，占比17.46%，分布于呈貢區的東北及東南向，該區域基本為林地，故植被覆蓋度最高。總體來說，除水域覆蓋區，研究區整體呈現外圍植被覆蓋較高，并向中部區域遞減的趨勢，植被覆蓋度以中等為主。

圖2 呈貢區NDVI，FC年際變化曲線

圖3 呈貢區FC小波變化情況

以研究區DEM重分類柵格圖構建ROI，經統計高程與植被覆蓋度的關系，如圖5所示，海拔為1 700～2 800 m，植被覆蓋度總體呈現急增—緩降—緩增趨勢。約88.88%的區域介于1 800～2 200 m，該范圍植被覆蓋度值為0.4～0.8，說明在該海拔范圍內，大部分區域為中高等植被覆蓋度；海拔1 700～1 800 m的區域主要是水域，植被覆蓋度低；而海拔2 200 m以上區域面積占研究區面積少，但由于該區域大多為林地，故植被覆蓋度較高。

3.2.2 空間趨勢特征 采用趨勢分析法模擬研究區影像中每個柵格單元29 a的變化趨勢，得到植被覆蓋空間變化趨勢。經統計，K值域為-0.045 2～0.037 3，表明植被變化趨勢存在空間差異。參照不同的FC變化性質將趨勢線斜率(K)劃分為5個等級，并繪制變化趨勢空間分布圖(圖6)。研究區約39%的地區植被覆蓋得到改善，其中輕微改善和明顯改善區域面積分別為129.21，50.68 km2，占比為28.03%，10.99%，主要分布于呈貢區東部大片區域的山地等地區；約117.67 km2的區域存在退化趨勢，占25.52%，分布于呈貢新區，這是由于呈貢新區的發展，建筑用地增多，故植被覆蓋度退化。植被基本沒有發生變化的區域面積約為163.44 km2，占35.45%，主要分布于水域，其他較均勻地遍布于呈貢區。整體而言，研究區1990—2018年植被改善區域大于退化區域，近半區域植被覆蓋得到了改善，植被發展前景較好。

圖4 呈貢區FC空間分布

圖5 呈貢區FC與高程對應關系

3.2.3 空間波動特征 研究區FC變異特征空間分布如圖7所示，經統計近29 a植被覆蓋度呈較低及低波動的變化態勢合占35.32%，低值區主要分布于呈貢區東部和西部兩大片區域，這部分地區地類多以林地、水域為主，變化較為穩定；中值區面積為170.6 km2，占比37.01%，分布于呈貢城區外圍；高值區占27.67%，主要分布于呈貢城區，這與近年來加速城鎮化建設密切相關，因其將耕地、林地及裸地轉為建設用地，故其空間波動較大。總的來說，由于呈貢新區的建設，導致該區中部大范圍區域空間波動較大，而外圍區域波動較小且更穩定。

3.2.4 土地利用變化特征 以上主要側重于對呈貢區近29 a植被覆蓋變化進行探析，而關于其土地利用變化尚不明晰，針對此，提取研究區1990年和2018年的土地利用信息(圖8)，在植被覆蓋變化大格局下對土地轉移以及植被變化信息進行進一步剖析。由表2可知，過去29 a間研究區林地、建設用地增加，耕地和裸地減少，而水域基本無變化。對比植被覆蓋度，發現植被覆蓋度高的區域為未發生地類改變的林地，呈貢區中部大部分為中等植被覆蓋度區域，地類以耕地、裸地和建設用地為主，水域覆蓋度最低；對比植被覆蓋變化趨勢空間分布圖，植被退化的區域主要為耕地、裸地轉變為建設用地的區域，植被改善的區域為呈貢區東部大部分的裸地轉變為林地區域；對比植被覆蓋變異特征空間分布圖，可知波動較高的為中部新增的建設用地區域，及東部的部分裸地轉為林地的區域，較穩定的為水域及東部部分林地。由轉移矩陣可知，轉移比例最大的為耕地和建設用地，其次為裸地和林地，比例最小的為水域。轉換為林地主要由裸地、耕地和建設用地構成，其中裸地轉換比例達25.94%；高達31.14%的建設用地及26.93%的耕地轉換為裸地，這與呈貢區大面積規劃建設及土地整改密切相關，并影響該區域植被覆蓋度時空特征。

圖6 FC變化趨勢空間分布 圖7 FC變異特征空間分布

圖8 1990年及2018年土地利用類型表2 1990年、2018年研究區土地利用類型面積轉移矩陣

4 結論與討論

本文采用像元二分法計算植被覆蓋度得到了呈貢區1990—2018年逐年的植被覆蓋度分布情況，將小波變換及趨勢分析等方法用于分析其時空變化特征，采用面向對象的高精度分類法得到首尾兩期土地利用變化。研究得出29 a年間呈貢區植被覆蓋整體呈緩慢增加趨勢，研究區在空間上呈現外圍植被覆蓋較高，并向中部區域遞減的趨勢，植被覆蓋度以中等為主，占研究區面積一半以上。變化趨勢表明呈貢區29 a年間植被改善區域大于退化區域。由于呈貢新區的建設，導致該區中部大范圍區域空間波動較大，而外圍區域波動較小且更穩定。呈貢區植被退化的區域主要為耕地、裸地轉變為建設用地的區域，植被改善的區域為呈貢區東部大部分的裸地轉變為林地區域；波動較高的為中部新增的建設用地區域；轉換為林地主要由裸地、耕地和建設用地構成，較高轉換比例的建設用地和耕地轉換為裸地。

本研究由于Landsat數據的可獲取性，每年選取的數據是3—5月期間的數據，故年際變化存在一定誤差。1990年的土地利用信息是基于谷歌影像進行精度驗證，部分樣本點可能不正確，故驗證效果有待改善。雖本文研究周期長，對于植被覆蓋的時間特征進行了回歸分析，但未對其未來發展趨勢進行預測，且未考慮氣候因子對植被覆蓋度的影響。今后研究重點應是在如今極端氣候事件頻現以及人類活動干預增強的背景下，如何定量區分氣候變化、人類活動以及其他驅動因子對植被時空變化的貢獻率。