1982—2022 年黃土高原歸一化植被指數變化的時空特征

方德昆,閆小月,張學珍,羅海江

1.中國環境監測總站,國家環境保護環境監測質量控制重點實驗室,北京 100012

2.中國科學院地理科學與資源研究所,中國科學院陸地表層格局與模擬重點實驗室,北京 100101

植被覆蓋是地表生態環境的重要組成部分,是保障生態質量的基礎,也是實現“碳中和”的重要途徑。 植被具有涵養水源、保持水土、固碳釋氧和改善局地小氣候的功能[1-2];此外,植被結構的改善有利于野生動植物種群的增加,對生物多樣性的恢復起重要作用[3]。 因而,植被覆蓋狀況是區域生態環境質量評價的重要內容。

衛星遙感是區域植被覆蓋監測的主要技術手段。 基于衛星遙感技術發展了一系列植被指數用于監測植被生長狀況,如:歸一化植被指數(Normalized Vegetation Index,NDVI)、比值植被指數(Ratio Vegetation Index,RVI)、增強型植被指數(Enhanced Vegetation Index,EVI)和土壤調節植被指數(Soil-Ad justed Vegetation Index,SAVI)等。其中,NDVI 物理意義清晰,數學計算過程簡單,使用范圍最廣,成為我國生態環境質量評價的重要指標之一[4]。 與其他各指數相比,NDVI 具有產品完善、時間序列長等優點,尤其適合開展區域生態環境質量變化的長期研究。

黃土高原地處我國北方半干旱氣候區,自21世紀初國家實施“退耕還林(草)工程”以來,黃土高原土地利用結構和格局發生了深刻變化。 特別是十八大以來,在習近平生態文明思想的指引下,生態恢復和生態建設逐步融入區域發展戰略,山水林田湖草沙冰綜合治理措施逐步實施,區域生態環境顯著改善。 本文擬采用2 套衛星遙感的NDVI 數據產品(GIMMS 3g 和MODIS),通過數據融合的方法,建立一套長時間序列的NDVI 數據集,在此基礎上揭示1982—2022 年黃土高原NDVI 變化的時空特征。

1 研究方法

本研究采用了GIMMS NDVI 3g 和MODIS NDVI 數據集。 GIMMS NDVI 3g 是美國航空航天局(NASA)全球監測與模型研究(GIMMS)項目發布的第三代基于NOAA 系列衛星AVHRR 傳感器影像資料的NDVI 數據集,涵蓋1981 年7 月至2015 年12 月的全球NDVI 數據,空間分辨率為1/12°(約8 km),時間分辨率為1 d。 數據生產過程中已處理了校準損失、軌道漂移、火山爆發等不利因素造成的干擾。 目前該數據集已被廣泛應用于全球各地的長期植被覆蓋變化研究[5-7]。

MODIS NDVI 數據集是NASA 地球觀測系統(EOS)計劃發布的全球NDVI 數據集,由地球遙測衛星Terra 和Aqua 上的中分辨率成像光譜儀(MODIS) 影像資料合成。 本研究采用編號為MOD13A2 v6 的數據產品,該產品使用基于一定合成周期的最大值合成法消除云層和氣溶膠的干擾噪聲,涵蓋2000 年2 月至2023 年2 月的全球NDVI 數據,空間分辨率為1 km,時間分辨率為16 d。 目前該數據集已被廣泛應用于植被動力學研究[8-9]。

由于GIMMS NDVI 3g 和MODIS NDVI 數據集的時間跨度均不足以單獨覆蓋目標研究時段,故需融合2 套數據重構1982—2022 年的NDVI時間序列。 現有研究表明,GIMMS NDVI 3g 和MODIS NDVI 之間存在高度相關性[10-11],故本研究采用方差匹配法[式(1)和式(2)]重構覆蓋1982—2022 年且均一化的NDVI 時間序列。 方差匹配法的核心思想是通過調整其中一個數據集,使2 個時間序列重疊時段數據的平均值和標準差保持一致,從而確保重構序列前后2 段的方差可比。

式中:Mi表示i 年的MODIS NDVI;和分別表示重疊時段(即2003—2013 年)的MODIS NDVI和GIMMS NDVI 3g 的平均值; SM和SG分別表示重疊時段MODIS NDVI 和GIMMS NDVI 3g 的標準差;表示2003—2022 年MODIS NDVI 的平均值。

由于GIMMS NDVI 3g 的時間跨度比MODIS NDVI 更長,本研究以GIMMS NDVI 3g 為參考標準,調整MODIS NDVI 的平均值和標準差,使其在重疊時段與GIMMS 3g 的平均值和標準差一致。由于數據可用性原因,本研究采用2 個數據集的重疊時段2003—2013 年作為參考時段進行數據調整,之后通過組合1982—2002 年的GIMMS NDVI 3g 數據和調整后的2003—2022 年MODIS NDVI 數據,構建1982—2022 年均一化的NDVI時間序列。 此外,由于2 個數據集的空間分辨率不同,故在應用方差匹配法進行調整之前,需對MODIS NDVI 進行像元合并。 本研究采用面積平均法進行此升尺度操作,將MODIS NDVI 的空間分辨率由1 km 升至8 km,確保與GIMMS NDVI 3g 一致。

利用均一化的1982—2022 年NDVI 時間序列,首先分析逐像元的NDVI 變化趨勢,為突出“退耕還林(草)工程”以及十八大以來NDVI 變化的空間特征,針對1982—2019 年、1982—2002年、2003—2022 年、2012—2022 年4 個時段分別進行分析,變化趨勢采用最小二乘法進行擬合。采用土地利用類型變化情況解釋黃土高原NDVI變化的原因,其中土地利用數據空間分辨率為1 km,來自中國多時期土地利用遙感監測數據集[12],包括1980、2000、2010、2020 年4 個時間點。 最后,以5 a 為一個時間段,對比分析了各時段之間各等級NDVI 的面積變化及主要NDVI 等值線的空間位置變化。

2 結果與討論

2.1 區域NDV I 變化的時空特征及與土地利用變化的關系

1982—2022 年黃土高原暖季(5—9 月)NDVI時序變化如圖1 所示。 黃土高原暖季NDVI 呈現出明顯的增加趨勢,增長率約為0.03/(10 a),暖季NDVI 在2018 年達到最大值(0.55),比最小值(0.40,出現在1982 年)高出約38%。 黃土高原暖季NDVI 具有明顯的年際波動趨勢,這與年際的降水不均勻、旱澇、氣候變化有關聯[13-14]。 其增加趨勢表現出“先慢后快”的特征,具體可劃分為2 個明顯的時間段:1982—2002 年,暖季NDVI表現為緩慢的上升趨勢,速率為0.01/(10 a),增幅僅為0.06;2003—2022 年,暖季NDVI 快速上升,速率高達0.06/(10 a), 增幅約為0.12。1982—2002 年,暖季NDVI 均值僅為0.44,但由于暖季NDVI 持續上升,均值在2012—2022 年達到較高水平,為0.52。 暖季NDVI 的增加反應了植被覆蓋率的提高,表明近40 年來黃土高原地表持續綠化,2003—2022 年綠化效果更為明顯,又以2012—2022 年尤為顯著。

圖1 1982—2022 年黃土高原暖季(5—9 月)NDVI 時序變化Fig.1 NDVI’s temporal tendency in the warm season (May-September) of the Loess Plateau from 1982 to 2022

黃土高原各時間段暖季NDVI 變化趨勢的空間分布如圖2 所示。 總體來看,近40 年來黃土高原86.4%區域的暖季NDVI 呈上升趨勢,其余區域略有下降[圖2(a)]。 其中黃土高原中部,即山西西部、陜西北部和寧夏中部暖季NDVI 上升速率較快,約為0.03/(10 a) ;這些地區主要為沙漠(毛烏素沙漠)和山區(呂梁山脈)。 暖季NDVI 呈下降趨勢的區域面積較小且分布零散,最大下降速率約為-0.02/(10 a) 。 通過與同期土地利用類型的空間分布[圖3(a)、3(d)]比較,可以發現暖季NDVI 上升速率較快的區域與同期耕地轉換成林地、草地的空間范圍基本一致;而在暖季NDVI 下降趨勢的區域中,主要的土地利用變化則是林、草地轉變為耕地及建設用地。

圖2 黃土高原暖季NDVI 變化趨勢空間分布Fig.2 Spatial distribution of NDV I’s temporal tendency in the warm season of the Loess Plateau

圖3 黃土高原土地利用類型空間分布Fig.3 Spatial distribution of land use types in the Loess Plateau

以2002 年“退耕還林(草)工程”政策付諸實施為時間節點,之前(1982—2002 年)黃土高原暖季NDVI 變化趨勢并不明顯(圖1)。 雖然70.9%的區域暖季NDVI 呈上升趨勢,但上升速率較小,僅為0.01/(10 a) 。 暖季NDVI 呈下降趨勢的區域面積較小,主要分布在山西中西部、河南與山西交界處以及陜西中部,且下降速率也較小,大多不超過-0.01/(10 a) 。 同期土地利用變化甚微[圖3(a)、圖3(b)]。 “退耕還林(草)工程”政策實施之后(2003—2022 年),黃土高原93.8%區域的暖季NDVI 呈上升趨勢[圖2(c)]。 其中,自黃土高原東北至西南方向出現一條醒目的“綠化帶”,橫穿內蒙古南部、山西北部、陜西中部、甘肅東南部和寧夏南部,暖季NDVI 上升速率平均高達0.04/(10 a) 。 同時期,黃土高原發生了大面積的“耕地范圍縮減、林草范圍擴增”的土地利用變化[圖3(b)、圖3(d)]。 可見,“退耕還林(草)工程”是近40 年黃土高原暖季NDVI 增加的重要原因。

十八大提出“加強生態文明建設”的總體要求以來,黃土高原暖季NDVI 的增加趨勢愈加顯著(圖1)。 其中黃土高原東北部、西南部,即山西北部、陜西北部和甘肅東部上升速率較大,約為0.07/(10 a) ,其空間分布與2003—2022 年暖季NDVI 上升速率高值區基本吻合。 因此,2012—2022 年黃土高原暖季 NDVI 的顯著上升是

2003—2022 年黃土高原地表綠化的主要成因。同時期,暖季NDVI 呈下降趨勢的區域主要包括山西南部的運城盆地和洛陽盆地,平均下降速率約為-0.05/(10 a) 。 通過與同期土地利用變化的空間分布對比[圖3(c)、圖3(d)],可以發現暖季NDVI 增加區域與耕地向林、草地轉換區域范圍基本一致,暖季NDVI 降低區域則與城市、耕地擴張區域基本一致。 可見,耕地繼續向林地、草地轉換可能是這一時期暖季NDVI 呈上升趨勢的主要成因;而城建用地和耕地范圍的擴增,則可能是黃土高原東南部暖季NDVI 呈下降趨勢的主要成因。

由上述結果可知,黃土高原暖季NDVI 變化與同期土地利用變化密切相關。 由于2003—2022 年暖季NDVI 增速較大,2012—2022 年增速尤為顯著,因此選取黃土高原2000、2010、2020 年土地利用類型,分別對2003—2022 年、2012—2022 年暖季NDVI 增速大于0.05/(10 a) 區域的土地利用變化情況進行分析,旨在進一步認識暖季NDVI 快速增加背后的土地利用變化因素。 其中,2000—2020 年土地利用變化對應于2003—2022 年暖季NDVI 變化,2010—2020 年土地利用變化對應于2012—2022 年暖季NDVI 變化。 將土地利用類型轉換面積由大至小排列,取前6 位,如圖4、表1 所示。 總體來看,在暖季NDVI 增速大于0.05/(10 a)區域內,土地利用類型變化主要表現為耕地向草地、林地的轉換。 2000—2020年,耕地轉換為草地、耕地轉換為林地、草地轉換為林地的面積共計6.70×104km2,占同期土地利用變化總面積的43.9%,其中耕地轉換為草地的面積百分比最大(為29.3%)。 2010—2020 年土地利用變化特征與2000—2020 年基本一致,均表現為“耕地面積減少、林草面積增加”。 可見,這種變化趨勢是2 個時期暖季NDVI 快速增加的主要原因。 進一步觀察2010—2020 年土地利用類型變化的空間分布情況可知[圖4(b)],黃土高原東部暖季NDVI 增加的主要原因可能是草地向林地的轉換,而其他區域則是耕地向草地的轉換。

表1 2003—2022 年、2012—2022 年黃土高原暖季NDVI 增速大于0.05/(10 a) 區域的土地利用類型變化情況Table 1 Areas of each type of land use changes for regions with NDVI’s grow th rate greater than 0.05/(10 a) in the warm season of the Loess Plateau from 2003 to 2022 and from 2012 to 2022,respectively

圖4 2003—2022 年、2012—2022 年黃土高原暖季NDV I 增速大于0.05/(10 a) 的區域及其土地利用類型變化Fig.4 Areas with NDVI’s growth rate greater than 0.05/(10 a) and their land use changes in the warm season of the Loess Plateau from 2003 to 2022 and from 2012 to 2022,respectively

2.2 典型暖季NDVI 等值線變化的主要特征

為了進一步探究近40 年來黃土高原暖季NDVI 的動態變化情況,以5(a)為一個時間段,對比分析了各時段各等級暖季NDVI 的面積變化及暖季NDVI 等值線的空間位置變化(圖5)。 總體來看,以暖季NDVI 等值線0.4 為界,黃土高原西北區域等值線數值較小且分布稀疏,這是由于西北區域分布有沙漠(烏蘭布和沙漠、庫布齊沙漠和毛烏素沙漠)及高原(鄂爾多斯高原),土地利用類型以草地和荒漠為主,植被覆蓋度較低。 而東南區域等值線數值較大且分布密集,這是由于該區域分布著山脈(呂梁山脈、子午嶺山脈和六盤山脈)及盆地(臨汾盆地、運城盆地和洛陽盆地),栽培植被、闊葉林廣布,植被覆蓋度較高。

圖5 黃土高原暖季NDVI 等值線分布Fig.5 Distribution of NDVI contours in the warm season of the Loess Plateau

以5(a)為一個時間段,即1985—1989 年、1995—1999 年、2005—2009 年、2015—2019 年,通過前、后2 個時段的對比來看,暖季NDVI 等值線的空間位置發生了明顯變化,主要表現為:低值等值線(0～0.4)由中部向西北方向推進(圖5),導致暖季NDVI 小于0.4 的面積減少了11.81×104km2(表2);而高值(0.6 ～0.1)等值線高值由東南向中部推進,導致暖季NDVI 大于0.6 的面積增加了6.81×104km2。 謝佩君等[15-16]通過分析1987—2020 年的“植被綠度線”位置的變化發現,植被綠度線明顯向西北方向移動,尤其是2005—2010 年北移幅度最大,在空間格局上,這一變化主要集中在黃土高原中部區域的黃土丘陵區,即陜西省北部榆林和延安一帶。 這表明黃土高原中部地區暖季NDVI 持續增加,植被恢復效果最為明顯。 此外,暖季NDVI 等值線空間位置變化、各值域面積變化在各時間段并不相同,即在“退耕還林(草)”前(1985—1989 年、1995—1999 年),暖季NDVI 等值線空間位置變化不明顯、各值域面積變化較小;“退耕還林(草)”后(2005—2009 年、2015—2019 年),暖季NDVI 等值線空間位置變化較為明顯,各值域面積變化較大,與前文分析結果一致。

表2 黃土高原暖季NDVI 各值域區間的面積Table 2 Areas for each range of NDVI in the warm season of the Loess Plateau

黃土高原植被覆蓋度增加致使生態質量有所恢復轉好,植被覆蓋度和植被凈初級生產力呈增加趨勢,其中水源涵養服務持續轉好的面積占比為15.46%[17]。 2000—2008 年,黃土高原生態系統凈固碳能力從0.011 Pg 增至0.108 Pg,植被固碳量增加速率為每年9.4 g/m2[18]。 值得注意的是,暖季NDVI 在黃土高原南部、東部和青海東部一帶呈緩慢下降趨勢。 這可能是由于:黃土高原南部、東部人口較多,1990—2015 年人口增幅為29.83%[19];區域城鎮化發展迅速,耕地分布較為廣泛,受人類活動的影響,植被覆蓋發生一定程度的退化[20]。 青海東部是黃土高原向青藏高原的過渡地帶,干旱缺水的自然條件使農業生產較難達到可持續發展的必要條件[21]。 局部植被覆蓋度減少對生態系統服務功能構成了一定威脅[18]。

3 結論

采用方差匹配方法對GIMMS NDVI 3g 和MODIS NDVI 進行了融合,重構形成了1982—2022 年黃土高原的NDVI 時間序列。 研究發現其間黃土高原暖季NDVI 呈顯著上升趨勢,其中1982—2002 年黃土高原暖季NDVI 出現小幅度的增加,2003—2019 年(尤其是十八大以來)黃土高原暖季 NDVI 上升趨勢明顯, 且增速較快。2003—2019 年的綠化區域主要分布在黃土高原中部,并向東北、西南方向延展,這些區域是開展“退耕還林(草)”的重點區域;同時,暖季NDVI等值線低值向西北方向明顯推進。 說明近20 年的生態建設工程取得了明顯成效,地表植被覆蓋度大幅提升。 值得注意的是,最近20 年間黃土高原南部、東部和青海東部一帶,由于城市擴展、氣候變化等因素,NDVI 在部分區域稍有下降。