方德昆,閆小月,張學(xué)珍,羅海江

1.中國(guó)環(huán)境監(jiān)測(cè)總站,國(guó)家環(huán)境保護(hù)環(huán)境監(jiān)測(cè)質(zhì)量控制重點(diǎn)實(shí)驗(yàn)室,北京 100012

2.中國(guó)科學(xué)院地理科學(xué)與資源研究所,中國(guó)科學(xué)院陸地表層格局與模擬重點(diǎn)實(shí)驗(yàn)室,北京 100101

植被覆蓋是地表生態(tài)環(huán)境的重要組成部分,是保障生態(tài)質(zhì)量的基礎(chǔ),也是實(shí)現(xiàn)“碳中和”的重要途徑。 植被具有涵養(yǎng)水源、保持水土、固碳釋氧和改善局地小氣候的功能[1-2];此外,植被結(jié)構(gòu)的改善有利于野生動(dòng)植物種群的增加,對(duì)生物多樣性的恢復(fù)起重要作用[3]。 因而,植被覆蓋狀況是區(qū)域生態(tài)環(huán)境質(zhì)量評(píng)價(jià)的重要內(nèi)容。

衛(wèi)星遙感是區(qū)域植被覆蓋監(jiān)測(cè)的主要技術(shù)手段。 基于衛(wèi)星遙感技術(shù)發(fā)展了一系列植被指數(shù)用于監(jiān)測(cè)植被生長(zhǎng)狀況,如:歸一化植被指數(shù)(Normalized Vegetation Index,NDVI)、比值植被指數(shù)(Ratio Vegetation Index,RVI)、增強(qiáng)型植被指數(shù)(Enhanced Vegetation Index,EVI)和土壤調(diào)節(jié)植被指數(shù)(Soil-Ad justed Vegetation Index,SAVI)等。其中,NDVI 物理意義清晰,數(shù)學(xué)計(jì)算過(guò)程簡(jiǎn)單,使用范圍最廣,成為我國(guó)生態(tài)環(huán)境質(zhì)量評(píng)價(jià)的重要指標(biāo)之一[4]。 與其他各指數(shù)相比,NDVI 具有產(chǎn)品完善、時(shí)間序列長(zhǎng)等優(yōu)點(diǎn),尤其適合開(kāi)展區(qū)域生態(tài)環(huán)境質(zhì)量變化的長(zhǎng)期研究。

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

1 研究方法

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

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

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

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

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

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

2 結(jié)果與討論

2.1 區(qū)域NDV I 變化的時(shí)空特征及與土地利用變化的關(guān)系

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

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

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

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

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

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

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

表1 2003—2022 年、2012—2022 年黃土高原暖季NDVI 增速大于0.05/(10 a) 區(qū)域的土地利用類型變化情況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) 的區(qū)域及其土地利用類型變化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 等值線變化的主要特征

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

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

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

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

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

3 結(jié)論

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