基于Google Earth Engine的典型峰叢洼地石漠化時空演變與驅動因子分析
——以西疇縣為例

徐紅楓, 王 妍, 劉云根, 肖羽芯, 黎 舟

(1.西南林業大學, 昆明 650224; 2.云南省山地農村生態環境演變與污染治理重點實驗室, 昆明 650224)

我國西南巖溶地區是石漠化分布最廣的區域[1]，其主要地貌特征是地表起伏大、多山地[2]、多碳酸鹽巖地區，土壤層淺薄，土地貧瘠[3]。根據《巖溶地區石漠化綜合治理規劃大綱(2006—2015)》區劃，峰叢洼地石漠化治理區是我國石漠化綜合治理八大類型區之一，該區出露碳酸鹽巖古老、堅硬、層厚質純，且連片分布，水文系統具有典型的二元結構，地表水發達而地下水發育，因此缺水、少土，耕地資源匱乏，人地矛盾的沖突尤為嚴重[4]。西疇縣地處云南省文山壯族苗族自治州中部，位于云貴高原向廣西盆地過渡的斜坡地帶，是典型的峰叢洼地地區，石漠化問題一直制約著西疇縣的經濟、社會、環境可持續發展[5]。30年來，西疇縣持續進行石漠化治理工作，人們通過退耕還林、坡改梯、經濟林營造等一列石漠化治理工程[6]，探索出了“六子登科”、“五法治水”等石漠化綜合治理模式，不僅取得良好的生態和社會效益，西疇人民在治理石漠化所產生的西疇精神也被列為云南脫貧攻堅精神的典型代表。

近些年來已經有不少學者針對石漠化信息提取及驅動因子進行了相關研究，在石漠化信息提取方面蘇旺德等以ETM，DEM和地質巖性數據構建了南汀河石漠化[7];羅旭玲等認為基于NPP，NDVI、地表反照率和坡度可以較為精準地反演石漠化[8];為了揭示石漠化的驅動因子，李陽兵等基于人機交互解譯，對黔中高原后寨河地區的石漠化演變驅動因子進行了分析，得出較多坡耕地的存在是石漠化的驅動因素[9]；王正雄等基于地理探測器對巖溶槽谷地區的石漠化驅動因素進行分析，得出巖性和土地利用、巖性和坡度的組合共同驅動巖溶槽谷地區的石漠化形成[10]。以往的研究多使用傳統的方法進行數據處理，處理效率不高，在大區域和長時間跨度的研究上局限性較大，基于云計算的可用性，本文借助Google Earth Engine云計算平臺，以典型峰叢洼地石漠化地區西疇縣作為研究區，旨在揭示其1990—2020年30年間的石漠化演變規律和驅動因子，為后續的石漠化演變及更大尺度生態恢復研究提供一定的參考和技術指導。

1 研究區和研究方法

1.1 研究區概況

西疇縣位于東經104°22′—104°58′，北緯23°06′—23°37′，其下轄有2個鎮、7個鄉：西灑鎮、興街鎮、蚌谷鄉、蓮花塘鄉、新馬街鄉、柏林鄉、法斗鄉、董馬鄉、雞街鄉，總面積1 506 km2。其中石漠山區面積達99.9%，巖溶面積達75.4%，境內巖溶峰叢溶蝕洼地，一般連片出現，山巒起伏，上層基巖裸露，地形封閉破碎，數百個水淹洼塘和大小不等的山間小盆地鑲于叢山巒障中。且洼地多為橢圓或條形，寬數十米至數百米；溶洞、豎并、天窗等巖溶現象較為普遍。

1.2 數據源與研究方法

1.2.1 Google Earth Engine簡介 Google Earth Engine(GEE)是谷歌旗下的云計算平臺，可以快速、批量處理大量數據[11]。GEE上提供了BP數量級的公開地球觀測數據，有超過500萬張影像，200多個公共數據集可供免費使用[12]。GEE的API接口支持JavaScript和python，可以通過編程調用數值、數組、矩陣、柵格數據和圖表等相關的多種模塊和函數，用于實現影像預處理、圖像分類、影像時間序列分析等操作。網址(https:∥earthengine.google.com/)。

1.2.2 數據源與研究方法 Landsat陸地觀測衛星具有30 m良好空間分辨率，且具有長的時間覆蓋，所以基于Google Earth Engine調用長時間序列的Landast影像作為數據源，對西疇縣1990—2020年近30年的石漠化演變進行分析。其中植被覆蓋度和巖石裸露率使用了Landsat 5/8 Surface Reflectance Tier 1數據集，該數據集已經使用LaSRC進行了大氣校正，道路和建筑用地的提取使用了Landsat 5/8 Collection 1 Tier 1 TOA Reflectance數據集。坡度數據在GEE上基于SRTM Digital Elevation 30 m數據計算得到，年均降水量數據單位為0.1 mm，人口密度數據單位為人/km2，GDP數據單位為萬元/km2，以上柵格數據均下載自中科院資源環境科學與數據中心(網址為：http:∥www.resdc.cn)。

基于GEE的超強計算能力，以及其提供的BP級免費公開遙感數據，通過JavaScript API編程在線調用Landsat 5影像309景和Landsat8影像174景，以云分量20為閾值進行影像篩選、然后去云、進行影像集的自動拼接和裁剪、對水體進行掩膜，并基于像元二分模型計算得到研究區植被覆蓋度、巖石裸露率，通過計算指數和紋理特征實現道路和建筑用地的提取，最后將結果導出，下載到本地使用ENVI 5.3和ArcGIS 10.5進行后續處理，通過疊加制圖得到了西疇縣1990—2020年石漠化等級分布圖，整體流程見圖1。

圖1 石漠化制圖流程圖

2 石漠化指標

2.1 植被覆蓋度

(1) 歸一化植被指數(NDVI)。石漠化過程往往伴隨著植被的退化[13]。歸一化植被指數(NDVI)可以較好地反映地表植被狀況[14]，公式為：

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

(1)

式中:NDVI為像元歸一化植被指數值；NIR為近紅外波段；RED為紅光波段。

(2) 運用像元二分模型[15]。在年度合成NDVI頻率累積表上取頻率5%對應的NDVI值作為NDVImin，取累積頻率為95%的NDVI值作為NDVImax[16]，公式為：

FVC=(NDVI-NDVImin)/(NDVImax-NDVImin)

(2)

2.2 巖石裸露率

(1) 歸一化巖石指數(NDRI)。基于石漠化地區的地物特性，張曉倫提出用Landsat數據的短波紅外(SWIR)與近紅外(NIR)波段，計算歸一化巖石指數(NDRI)提取石漠化分布結果精度高[17]。公式為:

NDRI=(SWIR-NIR)/(SWIR+NIR)

(3)

(2) 巖石裸露率(Fr)。在年度合成NDRI頻率累積表上取頻率衛5%對應的NDRI值作為NDRImin，取累積頻率為95%的NDRI值作為NDRImax，計算公式為：

Fr=(NDRI-NDRImin)/(NDRImax-NDRImin)

(4)

2.3 建筑用地與道路提取

基于Google Earth Engine，以Landsat 5 Collection 1 Tier 1 TOA Reflectance(1990年、2000年和2010年)和Landsat 8 Collection 1 Tier 1 TOA Reflectance為數據源(2020年)提取西疇縣不透水面(建筑用地和道路)，計算歸一化差異水體指數(MNDWI)，來提取水體[18]，計算歸一化建筑指數(NDBI)來提取部分建筑特征[19]，并運用灰度共生矩陣(GLCM)計算其他紋理特征[20]，然后選取樣本點，取其中70%的樣本點作為訓練樣本，使用GEE上的隨機森林分類器進行訓練，另外的30%樣本點作為驗證樣本，生成混淆矩陣，并計算Kappa系數。樣本點分成兩類：第一類為不透水面，主要包括道路和建筑用地；第2類為透水面，主要包括植被、裸地、農用地水體等。樣本點的拾取都在GEE平臺上完成，2020年樣本的拾取基于Google高精度在線地圖，1990年、2000年、2010年樣本點的拾取則是基于Landsat5。

Kappa系數通常用于一致性檢驗，也可以用于衡量分類的精度，Kappa系數越大，分類精度越高。公式如下：

(5)

式中:po是每一類正確分類的樣本數量之和除以總樣本數,也就是總體分類精度；pe代表偶然一致性,通常Kappa系數在0～0.200代表極低的一致性,0.210～0.400代表一般的一致性,0.400～0.600代表中等的一致性,0.610～0.800代表較高的一致性,0.810～1代表幾乎完全一致。

不透水面提取結果如下：1990—2020年四期不透水面提取的總體精度分別為0.816，0.860，0.865，0.982，Kappa系數分別為0.627，0.707，0.730，0.963。1990年、2000年、2010年三期提取精度在0.62～0.80之間，2020年的提取精度在0.963四期提取結果一致性較高，基本可以滿足研究需要。

2.4 水體掩膜

歸一化植被指數范圍在-1至1之間，通常地表云、雪以及水體覆蓋時，近紅外波段的反射率低于紅波段的反射率，所以NDVI值為負數，所以在GEE上計算植被覆蓋度和巖石裸露率時對NDVI小于0的部分進行掩膜，從而達到剔除水體的目的。

在參考前人標準[21]的基礎上，將石漠化劃分為無、潛在、輕度、中度、重度和極重度石漠化6個等級，并剔除了水體、道路和建設用地，得到了1990年、2000年、2010年、2020年總共四期的石漠化等級分布圖(圖2)。在Google Earth Pro上驗證解譯精度，然后制作混淆矩陣，并計算Kappa系數，2020年使用的驗證影像來源為Google Earth高精度影像，1990年、2000年、2010年使用的影像數據源為Landsat 5。四期石漠化等級分布圖的總體精度分別為0.715，0.735，0.778，0.895，Kappa系數分別為0.681，0.671，0.745，0.721。1990—2020年四期石漠化等級分布圖精度均在0.610～0.800，一致性較高。

圖2 1990-2020年石漠化等級分布

3 結果與分析

3.1 石漠化時間變化特征

由圖3可知，第一階段1990—2000年的10年間各級石漠化面積都有一定減少，無石漠化土地增加明顯，增加了19.349%；在第二階段2000—2010年，2010年石漠化分布情況相較2000年又有所惡化，無石漠化土地面積減少明顯，減少了22.874%，重度石漠化面積增加了8.940%；在第三階段，2010—2020年期間，西疇縣的無石漠化地面積上升明顯，較2010年面積增加了46.015%，潛在、輕度、中度、重度石漠化土地面積較2010年顯著減少。

3.2 石漠化轉移矩陣

轉移概率矩陣是指客觀事物從一個狀態轉移到另一個狀態，可以清晰表示不同時期同一區域石漠化類型相互轉換的數量化關系[22]。為分析不同等級石漠化相互之間的轉化關系，利用ArcGIS 10.5構建西疇縣石漠化轉移矩陣，在Origin 2021作圖得到石漠化轉移圖(圖4)。

圖3 1990-2020年各級石漠化面積比重

從圖4可以看出：(1) 1990—2000年西疇縣的潛在石漠化轉向無石漠化轉變最顯著，面積160.809 km2；其次，變化較大的是輕度石漠化轉換成無石漠化面積84.767 km2；其余類型之間的轉換面積相對較小，總體來說石漠化都在呈好轉趨勢。(2) 2000—2010年，無石漠化轉變為潛在石漠化變化最顯著，面積144.746 km2；無石漠化轉變為輕度石漠化的面積次之，有115.608 km2；除此之外，潛在和輕度石漠化向中度石漠化轉化面積較大分別為72.322 km2，75.303 km2，同時也有較大面積的中度石漠化向重度石漠化轉化，有79.321 km2。(3) 2010—2020年主要趨勢為潛在、輕度和中度石漠化轉變為無石漠化石，面積分別為186.49 km2，195.916 km2，186.985 km2；同時也有66.168 km2重度漠化轉化為無石漠化；(4) 1990—2020年潛在、輕度、中度、重度、極重度石漠化均有向無石漠化轉變，無石漠化面積增加總計648.476 km2，截至2020年西疇縣無石漠化土地面積共計940.854 km2；修復難度較大的極重度石漠化也得到了很好的修復，向重度、中度、輕度及潛在石漠化轉變，其中有20.208 km2轉變為無石漠化。

3.3 石漠化空間演變分析

為探究西疇縣近30年石漠化空間演變格局，將石漠化等級降低的區域視為好轉，石漠化等級提升的區域視為惡化，石漠化等級未發生變化的區域視為不變(圖5)。總的來看，1990—2020年西疇縣的石漠化有51.715 km2發生惡化，有850.183 km2石漠化土地發生好轉，在空間上整體都呈現好轉，也有少部分土地呈現惡化，主要分布在在蓮花塘鄉西北部、新馬街鄉的東北部，雞街鄉和西灑鎮的交界處也存在部分零星的惡化土地。其中1990—2000年西疇縣的東北部石漠化改善較大，包括雞街鄉、董馬鄉和法斗鄉，西南部存在惡化，尤其是興街鎮和新馬街鄉，還有蚌谷鄉和蓮花塘鄉交界處；2000—2010年北部和西部大面積石漠化土地都存在不同程度的惡化，尤其是董馬鄉、西灑鎮、蚌谷鄉以及法斗鄉北部；2010—2020年西疇縣石漠化總體是發生改善，尤其是西灑鎮、蚌谷鄉和興街鎮，蚌谷長箐片區、興街江龍片區、羅漢沖片區、三光片區等石漠化綜合治理示范區的設立治理效果十分顯著。

圖5 1990-2020年石漠化變化

3.4 貢獻因子獲取與分析

地理探測器(GeoDetector)是由王勁峰團隊開發的地理探測器軟件[23]。其理論核心是通過空間異質性來探測因變量與自變量之空間分布格局的一致性，據此度量因變量對自變量的解釋度即q值。q的值域為[0，1]，值越大，說明這個因子對Y值的解釋力越強。其模型：

(6)

為分析西疇縣石漠化驅動因子，分別選取了西疇縣2010年坡度、年均降水量、人口密度、GDP作為自變量[10]。在ArcGIS 10.5上將4個變量柵格采用自然間斷法重分類成6類，然后用500 m網格點對石漠化類型分布柵格與變量柵格進行采樣提取。

從單因子探測結果來看GDP(0.626)>年均降水(0.577)>人口密度(0.056)>坡度(0.022)，多因子交互探測如表1所示，坡度和GDP的交互作用下，對石漠化的影響最大(q=0.645)，GDP和人口密度對石漠化形成的影響則次之(q=0.639)。

表1 多因子交互探測結果

4 討論與結論

4.1 討 論

Google Earth Engine(GEE)存儲有長時間序列的多源遙感數據，且計算能力強大，隨著近幾年云計算的興起，越來越多的學者開始使用GEE云計算平臺開展相關研究。張滔等[11]基于GEE對京津冀地區2001—2015年植被變化與城鎮擴張的關系進行研究實現了數據的快速處理；裴杰等[24]利用GEE計算長時間序列的Landsat-NDVI對南洞地下河流域1988—2016年29年的植被覆蓋度變化進行了監測研究；陳黔等[25]基于GEE與機器學習對毛烏素沙地灌木覆蓋度進行估算，有效、快速地區分了灌木并對其植被覆蓋度進行估算。在巖溶地區的大尺度石漠化監測方面，GEE可以解決數據收集困難、數據預處理復雜和效率低下的問題，與傳統的石漠化遙感信息提取相比，利用GEE云計算平臺提取石漠化表征因子，只需要通過代碼的編寫，就能進行大尺度遙感數據的批量處理與計算，使得在相同數量級數據處理下效率得到大幅提升，GEE勢必成為實現大尺度石漠化監測的有效工具。

經過近30年的石漠化治理，2020年西疇縣的無石漠化土地面積已經達到總面積的62.474%，石漠化治理取得了很大成效。1990—2000年的10年間，石漠化得到了一定的治理但是成效并不明顯，可能的原因是石漠化治理雖然已經開展，但是由于缺乏統一的規劃和資金，所以成效甚微。2000—2010年，石漠化情況有所惡化，此時西疇縣進入了經濟發展的新階段，峰叢洼地大量耕地轉化成建筑用地[26]，人們為謀求更多的生存空間，任意開采開荒，破壞林草植被，加劇了石漠化形成；后10年時間(2010—2020年)，西疇石漠化治理成效顯著，這和政府的規劃有很大的關系，西疇縣一期石漠化綜合治理工程(2011—2014年)，對磨合、三光、長箐3個小流域實施了封山育林、植樹造林、坡改梯等措施，2016年制定《西疇縣全面深化生態文明體制改革實施方案》，并相繼實施了蚌谷長箐片區、興街江龍、羅漢沖、三光片區等石漠化綜合治理工程，探索出了“六子登科”、“五法治水”的石漠化綜合治理模式。鑒于植物措施是鞏固土壤、減少水土流失的有效措施[27]，后續應當繼續推進石漠化治理工程，并加大對植樹造林與經濟林營造的扶持力度，在治理石漠化的同時增加人民的收入。

在驅動因子方面GDP(q=0.626)和年均降水(q=0.577)對西疇石漠化的形成貢獻較大，GDP的貢獻度最高，說明經濟的發展對自然資源的影響較大，如果不能處理好經濟發展與生態保護的關系，極易加劇石漠化的進程；其次是降水，其決定了植被生長的條件，尤其是西疇縣屬于典型的峰叢洼地，極度缺水，降水對植被的生長尤為重要。坡度對西疇縣石漠化形成的貢獻并不大(q=0.022)，但坡度與GDP交互作用下對石漠化的影響最大(q=0.645)，通常情況下，在適合人類活動的坡度范圍內，坡度越陡，人類的活動越容易造成水土流失[28]；人口密度和GDP的交互作用對石漠化的形成影響則次之(q=0.639)，一般人口密度越大、GDP越高，經濟活動越頻繁，越容易造成植被破壞與水土流失，形成石漠化。

4.2 結 論

(1) 從時間上來看1990—2020年西疇縣石漠化演變的總體趨勢是先小幅減少，再擴張，然后大幅減少。1990—2000年石漠化轉移的主要方向是潛在石漠化和輕度石漠化向無石漠化轉移；2000—2010年，石漠轉移方向是無石漠化向潛在和輕度石漠化轉移；2010—2020年潛在、輕度和中度石漠化向無石漠化石轉移最顯著。從空間上來看西疇縣石漠化分布北重南輕，在石漠化治理成效上，經過近30年的治理西疇石漠化整體都呈現好轉，石漠化改善面積為850.183 km2，尤其是西灑鎮、蚌谷鄉和興街鎮，蚌谷長箐片區、興街江龍片區、羅漢沖片區、三光片區等石漠化綜合治理示范區治理效果十分顯著。在蓮花塘鄉西北部、新馬街鄉的東北部、雞街鄉和西灑鎮的交界處也存在零星惡化，惡化面積為51.715 km2。

(2) 截至2020年，西疇縣無石漠化土地面積已經達940.854 km2，占西疇縣總面積的62.474%，無石漠化面積增加總計648.476 km2。說明西疇縣的石漠化治理工作較為成功，西疇縣人民因地制宜探索出的“六子登科”、“五法治水”等石漠化綜合治理模式效果明顯。

(3) 地理探測器的結果顯示年均降水量(q=0.626)和GDP(q=0.577)是驅動西疇縣石漠化形成的主要驅動因子，多因子交互探測結果顯示坡度和GDP交互作用(q=0.645)，人口密度和GDP交互作用(q=0.639)這兩對組合共同驅動石漠化形成。

(4) 鑒于Google Earth Engine(GEE)云計算平臺提供的BP級海量遙感數據和強大的計算能力，本研究利用GEE進行石漠化表征因子的快速提取，分析西疇縣1990—2020年石漠化的動態變化及驅動因子，旨在利用GEE快速地實現大尺度地學問題分析，為后續大尺度的石漠化遙感制圖提供有力的技術支持。