張永,楊自輝,郭樹江,王強強,詹科杰,張劍揮,魏懷東
(1.甘肅農業大學林學院,甘肅 蘭州730070;2.甘肅民勤荒漠草地生態系統國家野外科學觀測研究站,甘肅 民勤733300)
綠洲防護帶(農田防護林,防風固沙林)固定沙丘、防止流沙向綠洲蔓延、消減沙塵暴發展,是保護綠洲的生態屏障。因此,綠洲防護帶植被的動態監測與評價,是探討綠洲生態安全的重要手段,傳統的監測,主要通過人工實地考察,建立定量觀測樣地調查等,不適合較大范圍定量監測。近幾年,隨著國內外爭相發射資源衛星,如國產高分系列,資源衛星系列,歐洲航空局哨兵衛星系列,利用衛片和航拍進行荒漠化動態監測成為首選的方法之一。
植被覆蓋度(fractional vegetation cover,FVC)是指植被冠層的垂直投影面積占區域總面積的百分比,是描述防護帶(包括天然林和人工林)長勢與密度的重要指標。目前通過衛星數據進行植被覆蓋度的研究,主要從數據源的選擇、計算指數的選取、變化成因3個方面開展。傳統的數據源主要是中低分辨率數據,如Modis、Landsat系列數據[1-2]。近幾年,隨著高分數據[3]存量的增加,以及小型無人機[4]的快速發展,高分辨率數據作為數據源成為熱點。此外,利用ASD光譜儀,野外實測獲取光合、非光合植被、裸土、陰影的端元光譜信息,使用線性光譜模型計算覆蓋度[5]。計算指數的選擇關系到覆蓋度估算的精度。因此,基于NDVI[6](normalized difference vegetation index)歸一化植被指數、EVI[7](enhanced vegetation index)增強型植被指數、TAVI[8](terrain adjust vegetation index)地形調節植被指數等基于像元二分模型[9-10]計算植被覆蓋度被廣泛采用。
土地利用變化遙感監測,能夠定量反演研究區各土地利用類型在時間和空間上的動態變化。利用衛星數據進行土地類型分類,其精度成為關鍵。目前,主要從分類數據源、分類方法和分類器選擇等方面提高分類質量。對于大尺度區域來說Modis、Landsat[11]成為首選數據。中小尺度區域利用高分數據、無人機航片[12]比較適合。基于CART決策樹[13-15]、面向對象分類方法[16]、支持向量機法[17]等分類器可以提高分類精度。此外,利用紋理和地形[18]、歷史土地利用數據[19]作為輔助,也成為分類過程中的重要參考指標。
本研究以民勤綠洲區及綠洲邊緣防風固沙林為研究對象,通過植被覆蓋度和基于CART自動決策樹分類技術組合多源數據進行土地利用分類,開展防護林帶監測與評價,探討民勤綠洲防護帶植被變化趨勢,為綠洲防護帶建設提供基礎數據支撐。
民勤綠洲位于甘肅河西走廊東段,是石羊河延伸到沙漠腹地出現的非地帶性景觀,是指石羊河下游紅崖山水庫以下沖積成的狹長而平坦的綠洲帶。甘肅民勤縣東西北三面被巴丹吉林和騰格里沙漠包圍,近九成面積為荒漠化土地,是中國重要的沙塵暴策源地之一。民勤綠洲像插入沙漠中的楔子,阻擋兩大沙漠的合攏。20世紀50年代,民勤綠洲邊緣(這里指主要受西北風危害的綠洲西北邊緣一線)地下水位為1~3 m,主要植被為紅柳(Tamarixramosissima),俗稱柴灣,植被生長旺盛,綠洲有石羊河多條水系的地表徑流滋潤,生態環境良好。1958年建紅崖山水庫,流域終端湖青土湖來水斷流。而后,隨著流域農業發展,耕地面積逐步擴大,地表徑流不能滿足流域農業生產需求,綠洲內來水逐漸減少,靠打井提灌滿足農業生產,地下水位逐漸下降,時至今日,地下水位下降至25 m左右,植被無法利用,造成衰敗死亡。在這個過程中,為了綠洲不受風沙危害,20世紀60年代,開始不間斷在民勤綠洲邊緣大面積治沙造林,2007年石羊河流域重點治理工程實施,促進綠洲防護帶建設進展,幾十年的防護帶建設,綠洲植被如何變化,需要有一個量化的說明。本研究主要探討該區20年來植被變化。研究區范圍是結合衛星遙感圖,通過綠洲外圍分布的固沙林實地調查生長狀況,包括植被的蓋度、密度等來確定,見圖1。

圖1 民勤綠洲與固沙林區Fig.1 Minqin oasis and sand-fixing forest
Landsat TM數據分辨率為30 m,這種數據對于中尺度區域開展遙感定量監測,具有成本低,處理周期短,效率高等優點。因此,本研究選擇1996、2006和2016年Landsat TM、OLI影像。云覆蓋均小于10%,成像質量較好,時間范圍為7月末至8月初,每期3景。由于選用的數據屬于未經處理的原始數據,因此,先對9景數據進行輻射校正和大氣校正,并采用二次多項式進行幾何校正,將校正誤差控制在一個像元以內。
選用12景2016年8 m多光譜、2 m全色高分1號數據,時間范圍在6-9月,通過大氣校正、正射校正等處理最后融合為2 m多光譜影像。對比兩種分辨率下計算結果的誤差大小。所有數據經過鑲嵌和裁剪處理,獲得研究區影像。
1.3.1植被覆蓋度的計算 在原有沙地、裸地的基礎上栽植人工固沙林,以及在低覆蓋度土地栽植固沙林恢復植被,在衛星遙感圖上表現為植被覆蓋度的增加。因此利用遙感影像反演植被覆蓋度,監測和評價綠洲防護帶的植被變化。本研究采用基于NDVI的像元二分模型[9]計算植被覆蓋度。計算公式為:
(1)
式中:fc為植被覆蓋度,NDVI是通過計算得到的3期研究區NDVI實際值,NDVIveg和NDVIsoil分別對應研究區范圍內純植被的NDVI值和土壤的NDVI值。根據研究區NDVI值的分布范圍,以置信度5%分別獲取3期NDVI上下閾值,作為NDVIveg和NDVIsoil。最后利用Arcgis中的重分類工具,結合民勤綠洲植被類型實際分布和荒漠化地區土地利用類型劃分方法,將植被覆蓋度劃分為5級(包括低植被覆蓋度0~5%、中低植被覆蓋度5%~15%、中等植被覆蓋度15%~30%、中高植被覆蓋度30%~60%、高植被覆蓋度60%以上)。此種分類方法,較為精確的反映干旱荒漠地區各覆蓋度類型空間分布和面積變化。
1.3.2基于CART決策樹分類與面積變化統計 CART決策樹分類技術,是通過對測試變量和目標變量構成的訓練數據集的循環二分形成二叉樹形式的決策樹結構[20]。本研究結合野外實測數據選取訓練樣本和多源數據,基于ENVI 5.3的CART計算工具,對民勤綠洲進行土地利用類型分類。并對分類后結果進行面積統計,最后建立3期影像各地物類型面積轉移矩陣。多源數據[14]相應計算公式為:

圖2 1996年植被覆蓋分布Fig.2 Vegetation cover distribution in 1996

圖3 2006年植被覆蓋分布Fig.3 Vegetation cover distribution in 2006

圖4 2016年植被覆蓋分布Fig.4 Vegetation cover distribution in 2016

(2)
(3)
(4)
式中:NDVI為歸一化植被指數;NDWI為歸一化差異水體指數(normalized difference water index);NDBI為歸一化裸露指數(normalized difference bare index)。bandgreen和bandswir1、bandnir、bandred分別對應TM影像中的綠波段、中紅外波段、近紅外波段和紅波段。NDVI可以消除地形和大氣的影響,對植被信息較為敏感。NDWI可以較好地消除土壤和植被背景的影響,從而突出水體信息;NDBI可以有效地分辨裸地信息。
根據植被覆蓋度計算方法,利用Arcgis進行密度分割和制圖,獲得3期植被覆蓋度圖(圖2~4)。
根據3期植被覆蓋度圖,利用Arcgis重分類工具,對各覆蓋度類型的面積計算統計。統計結果見表1。
根據圖2和表1可以看出:1996年民勤綠洲防護帶的低覆蓋度植被主要分布在綠洲外圍,占總面積的50.57%。中低覆蓋度植被主要分布在中部、西南部地區,占總面積的25.90%。中等覆蓋度植被主要分布在西南部、中西部、東南部地區,占總面積的6.48%。中高覆蓋度植被主要分布在農田周邊和西南部地區,占總面積的5.49%。高覆蓋度植被主要分布在中部農田附近,占總面積的11.56%。
根據圖3和表1可以看出:2006年民勤綠洲防護帶的低覆蓋度植被主要分布在綠洲外圍且面積較大,范圍較廣,占總面積的67%,中低覆蓋度植被主要分布在中部農田附近和東南部地區,所占比例8.73%,中等覆蓋度植被分布在東南部、西南部、中部地區,所占比例5.32%,中高覆蓋度植被分布在東南部、中部綠洲內部,所占比例5.35%,高覆蓋度植被主要呈條帶狀分布在中部農田附近,所占比例3.62%。
圖4和表1看出2016年各覆蓋度類型:低覆蓋度植被主要分布在西南部和中南部,所占比例49.39%,中低覆蓋度植被主要分布在東部、中部、北部、東南部,范圍較廣,所占比例26.68%,中等覆蓋度植被分布在東南部和中部,所占比例5.53%,中高覆蓋度植被分布在東南部、西南部、中部地區,所占比例4.06%,高覆蓋度植被主要分布于東南部和西南部農田附近以及中部地區,所占比例5.22%。

表1 1996、2006和2016年植被覆蓋度等級面積及所占比例Table 1 1996, 2006, 2016 vegetation coverage area and proportion

表2 1996、2006和2016年不同植被覆蓋度等級面積變化及所占比例Table 2 1996, 2006 and 2016 different vegetation coverage grade area changes and the proportion
結合圖1~4和表2可以發現,1996-2006年10年間,低覆蓋度和高覆蓋度面積有所增加,分別增加32.48%和27.65%。中低覆蓋度、中等覆蓋度、中高覆蓋度面積有所減少,分別減少66.28%、17.94%、2.61%,重點監測樣區中,劉家地和老虎口由原來的中低覆蓋度退化為低覆蓋度,南湖鄉部分區域由中低覆蓋度轉變為低覆蓋度。2006-2016年10年間,低覆蓋度和中高覆蓋度面積有所減少,分別減少26.28%和24.00%。中低覆蓋度、中等覆蓋度、高覆蓋度面積增加,分別增加了205.47%、4.12%、44.22%。重點監測樣區中,青土湖區由低覆蓋度改善為中低覆蓋度和中等覆蓋度,劉家地和老虎口由低覆蓋度轉變為中低覆蓋度,南湖鄉由低覆蓋度轉變為中低覆蓋度。總體來看1996-2016年20年間植被覆蓋度總體呈現上升趨勢,整個主綠洲區生態環境朝著良性發展。

圖5 1996年土地利用Fig.5 Land use in 1996
本研究主要監測綠洲防護帶的空間分布和面積變化,因此將民勤綠洲土地利用類型分為3類:農田防護林[新疆楊(Populusalbavar.pyramidalis)]、防風固沙林[梭梭(Haloxylonammodendron)、白刺(Nitrariatangutorum)、檉柳(Tamarixchinensis)]、其他地類(農田、水域、裸巖、鹽堿地、沙地)。此外,利用本研究采用的土地利用分類方法得到分類結果(圖5~7)。
對各土地利用類型面積進行計算統計,結果見表3和表4。
根據圖5和表3可知1996年各土地利用類型的空間分布和具體面積。其中農田防護林主要呈條帶狀零星分布在綠洲中部、中北部地區,占總面積的0.29%,防風固沙林主要分布在綠洲外圍,分布范圍較廣,占總面積的63.64%,其他地類主要呈狹長形分布于綠洲區外圍和綠洲中部地區,占總面積的36.07%。

圖6 2006年土地利用Fig.6 Land use in 2006

圖7 2016年土地利用Fig.7 Land use in 2016
結合圖6和表3可以看出,2006年各土地利用類型:農田防護林主要分布于中部和中北部地區,分布范圍較小,占總面積的0.72%,防風固沙林主要分布于綠洲西北部、東部邊緣地帶,東南部有小范圍集中分布,所占比例74.10%,其他地類主要分布在綠洲中部、北部、東部邊緣地帶,東南部有少量分布,所占比例25.18%。
結合圖7和表3可以看出,2016年各土地利用類型中,農田防護林主要分布于綠洲中部、中北部,占總面積的0.47%。防風固沙林主要分布于中部農田外圍大部分地區,分布范圍較廣,面積較大,占總面積的76.73%。其他地類小范圍分布于綠洲中部、北部、東南部、西部地區,占總面積的22.79%。
根據圖1、5~7和表4看出,1996-2006年農田防護林面積增加了25.75 km2,增幅148.68%,說明在這10年間,阻止沙漠侵吞農田,開展農田防護林建設,取得一定成果。防風固沙林分布范圍明顯向綠洲邊緣擴大,增加面積630.12 km2,增幅達16.47%。此外,老虎口和南湖鄉,防風固沙林面積增加明顯,青土湖和劉家地以沙地等其他地類為主,沒有較大變化。2006-2016年的10年間,農田防護林面積出現減少,分布范圍發生萎縮,實際減少14.64 km2,減少比例33.99%。防風固沙林面積分布范圍進一步擴大,實際增加面積有158.38 km2,增幅3.55%。此外,青土湖、劉家地、老虎口,由裸地改善為防風固沙林地,南湖鄉也有部分裸地轉變為防風固沙林。總體上看,1996-2016年的20年間,農田防護林出現先快速增長后小幅減少趨勢。防風固沙林分布范圍和面積出現先快速增長后緩慢增長的趨勢,說明民勤綠洲,20年來人工造林工程成果顯著,有效遏制沙漠向綠洲推進。

表3 1996-2016年土地利用面積統計Table 3 Land use area statistics from 1996 to 2016

表4 1996-2016土地利用面積變化Table 4 Change of land use area from 1996 to 2016
通過表5可以發現,2006年較1996年,農田防護林面積有52.53%保持不變,有近一半的面積發生轉移,其中大部分面積轉化為防風固沙林和其他地類,所占比例分別為26.57%和20.90%;防風固沙林有88.97%的面積保持不變,其余主要轉變為其他地類,所占比例10.36%。2016年較2006年,農田防護林只有10.87%的面積保持不變,主要轉變為防風固沙林和其他地類,所占比例65.77%和23.38%;防風固沙林有82.33%面積保持不變,其余主要轉變為其他地類,所占比例17.42%。
高分1號融合后的2 m多光譜數據,較中低分辨率數據,具有地物紋理和輪廓更加清晰的優點,可以較為準確的判斷地物類別以及更加準確的獲取純植被和純土壤的NDVI值。因此,適合基于像元二分模型求取植被覆蓋度結果和中低分辨率土地利用分類結果的精度驗證。本研究將處理后的高分1號數據利用ENVI進行波段計算,獲取研究區NDVI值的范圍,計算植被覆蓋度用于對Landsat TM/OLI數據處理結果的精度驗證。同時在高分1號影像均勻選取每類地物的樣點(可以彌補中尺度區域選點困難,時間和人力成本高的缺點),結合野外實測數據建立驗證樣本,對基于CART決策樹分類方法得到的土地利用分類結果進行精度驗證。
根據表6植被覆蓋度精度驗證,可以看出,低覆蓋度、中低覆蓋度和中等覆蓋度分別有83.14%、87.83%和84.47%被正確分類。中高和高覆蓋度分類精度較低,有79.17%和73.42%被正確分類。總體分類精度為81.62%,Kappa系數為0.78。說明基于二分像元模型計算植被覆蓋度,總體分類精度良好,但是由于TM數據分辨率較低,混合像元的存在對于提取純植被和純裸地NDVI值有一定的干擾。
根據表7土地利用分類精度可以發現,防風固沙林分類精度最高,有91.58%被正確分類。農田防護林分類精度最低,只有65.21%被正確分類,這是因為農田防護林為零星分布,一般寬度為1~3 m,而TM數據單位像元的邊長為30 m,所以很難有純凈的農田防護林像元存在,只有范圍較大的林區和經濟林田可以較好地被分離出來。總體分類精度為74.95%,Kappa系數為0.71.說明總體分類精度較好。

表5 綠洲防護帶面積轉移矩陣Table 5 Oasis protection zone area transfer matrix (%)

表6 植被覆蓋度精度驗證Table 6 Verification of vegetation coverage accuracy (%)
注:總體分類精度81.62%,Kappa系數:0.78。
Note:Overall accuracy 81.62%; Kappa coefficient 0.78.

表7 土地利用分類精度驗證Table 7 Land use classification accuracy verification (%)
注:總體分類精度74.95%,Kappa系數:0.71。
Note:Overall accuracy 74.95%; Kappa coefficient 0.71.
甘肅民勤地區處于石羊河下游,三面受巴丹吉林和騰格里沙漠包圍,是中國四大沙塵源區之一。特殊的地理環境,迫使民勤人民長期開展綠洲生態環境治理,2007年以來石羊河流域重點治理規劃實施,治理效果、生態環境變化趨勢等是人們關注的熱點問題,本研究利用Landsat TM結合高分1號數據,探討20年來防護林帶面積及其植被覆蓋度變化,研究結果與滑永春等[21]、韓濤等[22]、王紅霞等[23]的研究結論基本一致。
利用中等分辨率遙感數據,基于光譜特征進行監督分類,在建立分類樣本時容易出現“同物異譜”、“同譜異物”現象,導致各類別可分離度低,最終影響分類結果的精度。尤其是在荒漠地區,植被稀疏,覆蓋度低,獲取的植被NDVI值會低于一般理論值。因此,單一數據輔助分類效果不佳。CART決策樹分類,通過組合多源數據和目視解譯建立分類樣本進行土地利用分類,較傳統分類方法,精度容易把握,分類結果質量較高。
1996-2016 年間,民勤綠洲及其外圍防護林帶各等級植被覆蓋度中,低覆蓋度植被的面積先增后減,整體稍有減少,實際減少面積比例為2.33%,這一類型為流動沙丘、稀疏白刺群落以及一些沙地草本植物。中低覆蓋度(增加面積比例3.00%)、高覆蓋度(增加面積比例84.10%)植被的面積整體處于增加趨勢;中等覆蓋度(減少的面積比例14.56%)、中高覆蓋度(減少的面積比例25.98%)植被面積整體處于減少趨勢。民勤綠洲地下水位持續降低,1996、2006、2016年地下水位分別為14、21、23 m左右,已處于民勤生態地下水位臨界范圍(7~12 m)[24],靠地下水生長的植物(如檉柳、人工梭梭林)衰退死亡,一些種群(白刺)逐步適應依靠降水生長。低覆蓋度植被生長主要受降雨的影響,20年來變化不大;受地下水位下降的影響,原先營造的梭梭林退化,靠地下水位生長的檉柳林死亡,是中等覆蓋度和中高覆蓋度植被面積減少的原因;大面積營造固沙林、封沙育林育草、退耕還林等工程措施促進了低覆蓋度和高覆蓋度植被面積增加。
總體來看,20年來,民勤綠洲農田防護林面積先增加后減少,整體處于增加趨勢,增幅64.09%,面積增加了11.10 km2;防風固沙林面積整體在增加,先快速增長,后緩慢增長,增加了20.58%,增加面積為788.50 km2。隨著石羊河流域綜合治理規劃的實施,節水灌溉工程減少了綠洲防護林的生態用水量,農田防護林衰敗死亡,但大面積的生態經濟林建設,綠洲防護林面積增加。民勤生態環境建設,大面積固沙造林,綠洲外圍形成了近10 km寬的防風固沙林帶,防風固沙林面積增加。
在綠洲防護帶面積轉移變化中:1996-2006年10年間,有47.47%農田防護林轉變為防風固沙林和其他地類,所占比例分別為26.57%和20.90%。2006-2016年的10年,有89.13%農田防護林轉變為防風固沙林和其他地類,所占比例為65.77%和23.38%。1996-2006年的前10年,防風固沙林主要轉變為其他地類,所占比例為10.36%;2006-2016年的后10年,防風固沙林有17.42%轉變為其他地類。
利用同期高分1號數據計算植被覆蓋度,對Landsat TM數據結果進行精度驗證。植被覆蓋度精度驗證: 總體分類精度為81.62%,Kappa系數為0.78。說明基于二分像元模型計算植被覆蓋度,總體分類精度良好。土地利用分類精度驗證:防風固沙林達91.58%,農田防護林只有65.21%。