基于MODIS和CLCD數據的黃土高原土地利用變化檢測及其對植被碳吸收模擬的影響

吳 巧 麗,張 鑫 陽,蔣 捷

(北京建筑大學測繪與城市空間信息學院,北京 100044;自然資源部城市空間信息重點實驗室,北京建筑大學,北京 102616)

0 引言

為保護和修復生態系統,我國先后開展了三北防護林體系建設、天然林保護工程、退耕還林還草工程和沙漠化綜合治理等重大生態工程[1]。黃土高原作為我國第一批退耕還林還草試點地區,近20年來為提高我國森林覆蓋度作出了巨大貢獻[2,3],1999—2019年黃土高原的植被覆蓋率由31.6%提高到65%,有效遏制了該地區的水土流失問題[4],同時,實施退耕還林還草工程對于黃土高原地區生態系統碳匯能力提升具有重要作用[5,6]。目前,大量研究聚焦在土地利用/土地覆被變化(Land-Use and Land-Cover Change,LUCC)對黃土高原氣候、生態環境和碳循環等的影響[7,8],缺乏對長時序土地利用土地覆被(Land Use and Land Cover,LULC)產品的應用準確性研究,制約著黃土高原退耕還林還草工程生態效益評估的準確性。

LULC產品是陸面生態過程模型進行植被—大氣碳通量模擬的核心驅動數據,其產品精度直接影響區域尺度到全球尺度總初級生產力(Gross Primary Production,GPP)的模擬精度。目前,全球已經公開發布了多種不同時空分辨率的LULC產品,如MODIS LULC、CLCD LULC、GlobeLand和World Cover等[9,10]。然而,不同LULC產品的分類算法、樣本選取和驗證方法各異,導致產品應用于局部區域的精度存在差異,增加了區域尺度GPP估算的不確定性[11]。既有研究結果表明,MODIS LULC(MCD12Q1)數據在黃土高原地區的分類精度不足60%[9],導致基于MCD12Q1數據估算的GPP難以準確分析LUCC對黃土高原陸地生態系統碳循環的影響。

綜上,本文引入時空分辨率和分類精度更高的CLCD LULC產品,對比分析CLCD LULC和MODIS LULC產品在表達2001—2019年黃土高原LUCC時空變化特征上的差異,并進一步基于長時序CLCD LULC產品、氣象數據和多源遙感數據驅動GPP遙感機理過程模型,模擬黃土高原地區在退耕還林還草工程實施近20年來植被GPP的時空變化特征,定量分析LULC產品不確定性對黃土高原GPP估算結果的差異,以期提高黃土高原地區的碳匯估算精度,為精確量化國家重大工程生態效益以及政府調整和執行生態修復政策提供科學依據。

1 研究方法與數據來源

1.1 研究區概況

黃土高原涉及陜西、山西、甘肅、寧夏、河南、內蒙古和青海共7個省區,面積約6.4×105km2(圖1),地勢西北高東南低,地形復雜、溝壑縱橫,主要地貌類型有黃土丘陵、黃土山地、黃土臺地和河谷平原4類[12];植被類型從東南到西北依次為森林植被帶、森林草原植被帶、典型草原植被帶、荒漠草原植被帶和草原荒漠帶[13]。該地區屬于溫帶大陸性季風氣候,年平均溫度介于4～14 ℃,年均降水量在200～750 mm之間,跨越了半濕潤、半干旱以及干旱區3個區域[14]。黃土高原惡化的生態環境嚴重制約當地社會經濟可持續發展。為減緩水土流失,提高植被覆蓋度,黃土高原實施了坡改梯、淤地壩建設和退耕還林還草等一系列重大生態建設工程,LULC格局發生顯著變化[15-17],緩解了水土流失嚴重的問題,為提升陸地生態系統固碳能力奠定了堅實基礎[18]。

圖1 黃土高原地貌空間分布

1.2 LUCC變化像元檢測和FGM模擬方案設計

為精確估算退耕還林還草工程實施近20年黃土高原地區陸地生態系統碳吸收的時空變化特征,本文首先基于不同源LULC數據,對比分析2001—2019年黃土高原LUCC時空變化特征,將研究期內發生土地利用變化的像元作為主要研究對象[19],若19年土地利用類別均一致,則識別為“穩定像元”,反之則為“變化像元”,提取變化像元制作土地利用轉移矩陣;其次,基于MCD12Q1數據和CLCD數據,結合多源遙感數據和氣象數據驅動Farquhar GPP Model(FGM)模型,生產2001—2019年黃土高原地區500 m空間分辨率、月尺度的GPP遙感產品,以CLCD GPP產品分析近20年黃土高原地區GPP的時空變化特征;最后,以CLCD LULC產品估算的GPP為基準,定量分析MCD12Q1分類誤差對GPP空間分布、各類植被GPP年際變化特征、GPP年總量和年際變化趨勢的影響。具體流程如圖2所示。

圖2 技術路線

1.3 FGM模型及其驅動數據

FGM模型[20]基于Farquhar模型進行葉片尺度的光合作用速率模擬,利用二葉模型進行葉片到冠層尺度的光合作用速率轉換,并結合陰陽葉LAI比例估算冠層尺度GPP。FGM模型解決了Fick定律、Farquhar模型和Ball-Woodrow-Berry氣孔導度模型聯立求解胞間CO2濃度時,復雜度高、計算效率低、參數化難度高的問題,在保留模型機理性強優勢的同時,提高了模型運算效率,更適合大尺度GPP模擬。目前,FGM模型已成功應用于中國[20]和歐洲[21]的GPP估算,因此,本文選用FGM模型定量分析LUCC對黃土高原GPP估算結果的影響,模型詳細參數及運算過程參見文獻[20-22]。

FGM模型的驅動數據包括遙感數據、氣象數據和大氣CO2濃度觀測數據(表1)。其中,遙感數據包括MODIS LULC[23]和CLCD LULC產品[24]、全球陸表特征參量(GLASS)的葉面積指數(LAI)數據[25]和下行短波輻射(DSR)數據[26];氣象數據來自國家地球系統科學數據共享服務平臺—黃土高原科學數據中心,是中國學者基于機器學習算法、中國氣象站點觀測數據和CRUNCEP氣象再分析數據集生產的一套空間分辨率更高的氣象參數產品,主要包括氣溫、相對濕度以及飽和水汽壓差(VPD)[27-29];由于缺乏時空連續的CO2濃度產品,本文采用美國夏威夷島MLO觀測站高精度和長時間跨度的大氣CO2濃度實測數據,該數據可以捕捉大氣CO2濃度的平均長期變化趨勢[30]。

表1 FGM模型的驅動數據集

由于多源遙感數據的時空分辨率、投影坐標系、空間覆蓋范圍和存儲格式均存在差異,為獲取配套的模型驅動參數,需進行數據預處理,使不同輸入數據的時空分辨率和投影一致,為GPP模擬奠定數據基礎。本文將遙感數據轉為tiff格式后,對數據進行拼接和裁減,重投影到WGS84地理坐標系,采用最鄰近插值法將數據分辨率重采樣到500 m;氣象數據采用最鄰近插值法將空間分辨率重采樣到500 m,由于氣象數據時間分辨率為月尺度,本研究將LAI和DSR產品重采樣到月尺度。最終,本文利用表1數據驅動FGM模型,估算出2001—2019年黃土高原地區500 m空間分辨率的GPP月尺度數據,并將GPP產品和LULC數據均轉投影為Albert_Conic_Equal_Area等面積投影坐標系,避免數據投影和面積變形對GPP年總量估算精度的影響。

1.4 土地利用數據及重分類方法

1)MODIS LULC數據包括MOD12Q1(1 km空間分辨率)和MCD12Q1(500 m空間分辨率)兩種數據,本文選擇較高空間分辨率的MCD12Q1產品。根據監督決策樹分類方法,MCD12Q1產品包含5種不同的土地覆蓋分類方案,本文按照國際地圈生物圈(IGBP)計劃定義分為17類[31],包括常綠針葉林(ENF)、常綠闊葉林(EBF)、落葉針葉林(DNF)、落葉闊葉林(DBF)、混交林(MF)、郁閉灌叢(CSH)、稀疏灌叢(OSH)、多樹草原(WSA)、稀樹草原(SAV)、草地(GRA)、濕地(WET)、農田(CRO)、農田與自然植被混合地類(CRV)、冰雪(PSI)、城市建成區(URB)、荒地(BRN)和水體(WAT)。

2)CLCD LULC數據是由武漢大學楊杰等基于Google Earth Engine平臺1982—2021年多傳感器觀測數據,利用隨機森林分類器得到的LULC分類結果,空間分辨率為30 m[24]。在黃土高原地區CLCD LULC產品土地利用分類精度高于MODIS LULC產品[32]。針對MODIS和CLCD LULC產品分類體系的異同,參照全國土地利用遙感分類系統標準[33],本文結合黃土高原的地形地貌特點,將兩套LULC分類體系按一級大類重分類為草地、林地、耕地、灌叢和非植被5種地類[20](表2)。兩種LULC統一采用Albert_Conic_Equal_Area等面積投影坐標系,黃土高原空間范圍矢量數據來源于中國科學院地理科學與資源研究所(https://www.resdc.cn/,2023年3月9日獲得)。

表2 兩種土地利用與土地覆蓋(LULC)產品重分類前后的類別對應關系

2 結果分析

2.1 2001—2019年黃土高原LUCC時空分布特征

基于2001—2019年MODIS MCD12Q1和CLCD LULC數據得到的黃土高原各地類面積存在較大差異(圖3)。MCD12Q1數據(圖3a)顯示,2001—2019年黃土高原耕地、林地和灌叢面積分別增加20 688 km2、8 993 km2和2 269 km2,但草地面積減少量(約26 635 km2)超過耕地和灌叢增加的總面積,與退耕還林還草實際不符,可能是MODIS分類方法混淆了耕地和草地所致;CLCD LULC數據(圖3b)顯示,2001—2019年黃土高原林地(+14 025 km2)和草地(+7 262 km2)面積增加,耕地(-18 843 km2)、灌叢(-1 727 km2)和非植被(-717 km2)面積減少,更符合退耕還林還草實際。

圖3 基于MCD12Q1和CLCD LULC數據的2001—2019年黃土高原各地類面積變化

基于2001—2019年MCD12Q1和CLCD LULC數據得到研究區各地類變化像元空間分布(圖4),統計各地類轉換占比(圖5)可知:①MCD12Q1數據(圖5a)顯示,研究期間黃土高原變化像元地類主要為草地(占變化像元總面積的42.47%)、灌叢(占比16.63%)和耕地(占比13.86%)。其中,草地主要轉為耕地(25.97%)和灌叢(12.07%),灌叢主要轉為林地(6.69%)和耕地(6.56%),少部分轉為草地(3.09%),耕地主要轉為草地(6.84%)和灌叢(4.77%)。此外,還有約8.07%的變化像元發生在非植被區域,主要轉為草地(7.18%);林地像元變化最少(1.73%),主要轉為灌叢(1.47%)。②CLCD LULC數據(圖5b)顯示,黃土高原變化像元地類主要為耕地和草地,分別占變化像元總面積的31.82%和31.00%,耕地主要轉為草地(24.32%)和非植被(5.07%);草地主要轉為耕地(18.93%)和林地(7.79%);還有9.36%的變化像元發生在非植被區域,主要轉為草地(8.59%);林地像元變化較小,僅占變化像元總面積的1.26%,主要轉為耕地(0.81%)。

圖4 基于MCD12Q1和CLCD LULC數據檢測的2001—2019年黃土高原變化像元的空間分布

圖5 基于MCD12Q1和CLCD LULC數據檢測的2001—2019年黃土高原土地利用轉移矩陣

為更好對比兩種土地利用產品在黃土高原的分類精度,本文在甘肅省慶陽市(107°38′46″E,35°43′51″N)選取典型耕地和草地混合地貌區。通過兩種土地利用產品與Google Earth高分辨率影像對比,發現MCD12Q1確實對耕地和草地存在分類混淆現象,未能有效捕捉到Google Earth影像上的林地擴張現象;而CLCD LULC產品分類更精細且與Google Earth影像更接近(圖6)。

圖6 2001年和2019年黃土高原地區MCD12Q1產品、CLCD LULC產品和Google Earth影像對比

進一步根據分類精度更高的CLCD LULC數據計算2001—2019年黃土高原7個省區各地類變化面積(圖7)。可以看出,7個省區耕地面積均有所減少,陜西、甘肅和山西位居前三,耕地面積分別減少6 146 km2、4 650 km2和3 228 km2;7個省區的林地面積均有所增加,山西、陜西和甘肅位居前三,分別增加6 303 km2、4 143 km2和1 945 km2;草地面積增量較大的為內蒙古、甘肅、陜西和寧夏,分別增加5 879 km2,2 289 km2、2 001 km2和1 593 km2,而山西草地面積減少最多,為4 933 km2。由此說明,退耕還林還草主要發生在黃土高原中部,尤其在黃河中游段西側。

圖7 基于CLCD LULC的2001—2019年黃土高原7個省區各地類面積變化

2.2 基于FGM模型的黃土高原總初級生產力估算

由上述研究可知,基于MODIS MCD12Q1和CLCD LULC數據得到的各地類空間分布、年際變化特征均存在顯著差異,利用長時序CLCD LULC數據能更有效地捕捉退耕還林還草工程實施近20年來黃土高原地區LUCC的時空分布特征。因此,本文選取CLCD LULC產品估算2001—2019年黃土高原GPP,并以CLCD估算的GPP為標準,定量分析MCD12Q1數據分類誤差對黃土高原GPP估算結果的影響。

由圖8a可知,2001—2019年黃土高原GPP年均值為5 210.2 Tg C a-1,呈東南高、西北低的空間分布特征,最高值約為1 850 Tg C a-1,陜西、山西和河南GPP年均值明顯高于其他省區;由圖8b可知,林地GPP(年均值為104.2 Tg C a-1)對黃土高原GPP貢獻最大,其次為草地(年均值為85.1 Tg C a-1)和耕地(年均值為82.6 Tg C a-1),灌叢貢獻較小(年均值為2.4 Tg C a-1)。從年際變化趨勢看,黃土高原地區陸地生態系統的碳吸收能力顯著增強(年際趨勢為+6.4 Tg C a-2,P<0.001),其中草地(年際趨勢為+2.6 Tg C a-2,P<0.001)和林地(年際趨勢為+2.5 Tg C a-2,P<0.001)貢獻最大。

圖8 基于CLCD數據和FGM模型估算的2001—2019年黃土高原GPP空間分布和各類型植被GPP年總量的年際變化趨勢

2.3 LUCC不確定性對GPP估算結果的影響

圖9a顯示,基于MCD12Q1估算得到2001—2019年的GPP存在高值高估、低值低估的問題,即黃土高原南部地區(包括青海、甘肅、陜西省南部、山西、河南和寧夏南部)的MCD12Q1 GPP相對CLCD GPP偏大,北部地區(包括內蒙古、陜西北部和寧夏北部)MCD12Q1 GPP相對CLCD GPP偏小;估算結果顯示,MCD12Q1 GPP相對CLCD GPP在林地存在明顯低估現象(-38.3 Tg C a-1),在灌叢(+33.6 Tg C a-1)、草地(+28.9 Tg C a-1)和耕地(+2.3 Tg C a-1)存在高估現象;相對CLCD GPP估算結果,基于長時序MCD12Q1數據驅動FGM模型會給黃土高原地區GPP年總量模擬值帶來9.6%的高估;從圖9b看,2001—2019年MCD12Q1 GPP年總量約為5 712.8 Tg C a-1,年際增長趨勢顯著(+7.85 Tg C a-2,P<0.001),相比CLCD估算的GPP年際增長趨勢產生22.7%的高估。

3 結論與討論

本文基于MODIS MCD12Q1和CLCD兩類典型土地利用數據,定量分析了2001—2019年黃土高原LUCC時空變化特征,揭示了MCD12Q1土地利用分類精度不足對陸地生態系統碳匯估算精度的影響,主要研究結論如下:

1)由于數據來源、分類標準和空間分辨率不同,MODIS MCD12Q1和CLCD兩種LULC數據集顯示的黃土高原LUCC的時空變化特征存在較大差異。黃土高原是退耕還林還草等一系列生態修復工程實施的重點區域,CLCD LULC數據空間分辨率更高,分類更精細,精度更高[24,32],該數據顯示近20年來黃土高原耕地面積減少,林地和草地面積增加。在氣候變化背景下,黃土高原各類植被GPP呈現顯著增加趨勢,證明退耕還林還草工程產生了重要的生態效益。值得注意的是,雖然MCD12Q1數據被廣泛應用于全球GPP產品的生產,但其局地分類精度不一定高。例如,雖然MCD12Q1數據也顯示出黃土高原林地面積呈逐年上升趨勢,但表現出耕地面積增加和草地面積減少的特征,這與退耕還草的實際和谷歌高分影像目視解譯結果均不相符,證明長時序MCD12Q1數據對黃土高原地區耕轉林和耕轉草解譯精度較低,無法精確量化退耕還林還草工程的生態效益。

2)長時序CLCD LULC數據表明,2001—2019年黃土高原7個省區的土地利用變化規律基本一致,其中,山西和陜西退耕還林效果顯著,而內蒙古以退耕還草為主,這與已有研究發現的草地退化規律一致[34,35],草地有可能進一步轉化為林地[36]。綜上,耕草轉換是黃土高原的主要LUCC變化特征,部分地區存在耕轉草、草轉林的自然演替現象。

3)長時序CLCD數據能較準確捕捉黃土高原地區退耕還林還草現象,2001—2019年黃土高原GPP呈顯著增長趨勢,其中,林地和草地GPP增長趨勢最顯著,表明退耕還林還草工程對陸地生態系統碳吸收的促進作用。但植被光合作用速率受氣候變化、大氣CO2濃度增加、植被變綠、土地利用變化等眾多因素共同作用的影響,需進一步分析植被碳吸收增長的具體原因。

4)研究假設CLCD LULC數據分類準確,定量分析了MCD12Q1數據分類誤差對黃土高原GPP估算精度的影響。對比CLCD GPP,MCD12Q1 GPP在黃土高原南部和北部地區分別存在高估和低估現象,基于MCD12Q1數據估算的GPP年總量均值高估9.6%,GPP年際變化趨勢則高估近22.7%;從各地類看,MCD12Q1 GPP低估了林地的GPP而高估了耕地、草地和灌叢的GPP。

雖然CLCD LULC數據在中國的分類精度更高[9,32],但目前無法完全避免土地利用分類誤差,這也給GPP估算帶來了一定的不確定性,未來有必要從土地利用分類誤差向GPP模擬誤差傳遞的角度開展進一步研究。