徐 勇, 黃雯婷, 靖娟利, 張占奕, 李明杰, 歐昱賢, 盧夢緣, 竇世卿

(桂林理工大學 測繪地理信息學院, 廣西 桂林 541006)

植被作為地理環(huán)境的重要組成部分，在陸地表層系統(tǒng)中所占的比例很高，是生物圈及其生態(tài)系統(tǒng)的核心和功能部分，其變化不僅對全球及區(qū)域能量循環(huán)、物質(zhì)的生物化學循環(huán)具有重要的影響，同時也對區(qū)域及全球氣候變化有重要的貢獻[1-2]。相較于其他數(shù)據(jù)，遙感數(shù)據(jù)具有時間和空間上的連續(xù)性，故常被用來監(jiān)測全球和區(qū)域植被時空分布及變化趨勢[3-6]，為研究植被的空間分布格局和預測未來發(fā)展趨勢奠定了堅實的基礎(chǔ)。

歸一化植被指數(shù)(normalized difference vegetation index， NDVI)對植被生物物理特征十分敏感，作為植被生長狀況的重要度量參數(shù)，被有效使用在監(jiān)測植被動態(tài)變化當中。國內(nèi)外學者利用NDVI數(shù)據(jù)在區(qū)域和全球尺度研究植被覆蓋時空演變特征和植被覆蓋對氣候變化的響應機制，并取得了豐碩的成果。de Jong等[7]利用GIMMS NDVI3 g數(shù)據(jù)反演了全球植被覆蓋變化趨勢，并在像元尺度上探測了植被覆蓋變化的轉(zhuǎn)折點。結(jié)果表明，全球超過26%地區(qū)的植被呈單調(diào)顯著變化。方精云等[8]發(fā)現(xiàn)生長季延長和植被生長速度加快是中國植被覆蓋狀況改善的重要原因，且溫度上升和夏季降水增加是導致植被活動增強的主要氣候因子。孟丹等[9]發(fā)現(xiàn)京津冀地區(qū)植被覆蓋呈改善趨勢，且整體上植被生長與降水相關(guān)關(guān)系為正，與氣溫相關(guān)關(guān)系為負。晏利斌等[10]利用GIMMS NDVI3 g分析了京津冀地區(qū)生長季植被覆蓋時空變化特征及其與氣候因子的相關(guān)關(guān)系。結(jié)果表明，京津冀地區(qū)植被覆蓋總體呈上升趨勢，且降水對植被生長的作用強于氣溫。

在全球氣候變暖和人類活動加強的背景下，研究植被覆蓋動態(tài)變化及其與氣候變化之間的相關(guān)關(guān)系和響應機制具有重要的現(xiàn)實意義。已有研究大多側(cè)重于探究植被覆蓋時空變化特征[11-14]，而對預測其未來變化趨勢鮮有研究，且在分析植被NDVI與降水和氣溫相關(guān)關(guān)系時，并未顧及地理位置信息給降水和氣溫插值精度帶來的影響，降低了研究結(jié)果的精度。因此，本文基于2001—2019年MODIS NDVI時間序列，利用Theil-Sen Median趨勢分析法、Mann-Kendall顯著性檢驗法、R/S分析法和Person相關(guān)分析等數(shù)學方法，結(jié)合ANUSPLIN氣象插值模型，反演2001—2019年京津冀地區(qū)植被覆蓋時空演變特征，并預測其未來變化趨勢，進一步分析植被NDVI對降水和氣溫變化的響應機制和時滯效應，研究結(jié)果可為區(qū)域乃至全球植被覆蓋時空演變特征及其與降水和氣溫響應關(guān)系研究提供理論支撐，對京津冀地區(qū)林業(yè)生態(tài)工程的實施效果評估和生態(tài)環(huán)境建設(shè)結(jié)果評估具有重要的現(xiàn)實意義。

1 研究區(qū)概況、數(shù)據(jù)來源及研究方法

1.1 研究區(qū)概況

京津冀地區(qū)(113°04′—119°53′E，36°01′—42°37′N)位于華北北端，地表形態(tài)復雜、地形起伏較大，主要表現(xiàn)為東南低、西北高。京津冀地區(qū)由北京和天津兩個直轄市，以及河北省組成，總面積約為2.17×105km2。受半濕潤半干旱大陸性氣候的影響，四季分明，冬季低溫干燥，夏季高溫濕潤。年平均氣溫10.4～11.9 ℃，年累計平均降水量375.5～684.7 mm。西北部土地利用類型主要以林地和草地為主，東南部主要以農(nóng)用地為主。京津冀地區(qū)位于環(huán)渤海經(jīng)濟圈核心區(qū)，經(jīng)濟發(fā)展水平、城市化速度和人口密度均高于其他地區(qū)。人口總數(shù)占全國總?cè)丝诘?.07%(中國統(tǒng)計年鑒2019)，而面積僅占全國總面積的2.25%。受氣候變化和人類活動的影響，該地區(qū)生態(tài)環(huán)境較為脆弱，植被生長對氣候變化的響應敏感。

1.2 數(shù)據(jù)來源

1.2.1 MODIS NDVI遙感數(shù)據(jù) MODIS NDVI數(shù)據(jù)來源于美國國家航空航天局發(fā)布的MOD13Q1C6數(shù)據(jù)集，時間跨度為2001年1月至2019年12月，時間分辨率為16 d，空間分辨率為250 m，每年23個時相。首先，利用MRT(MODIS Re-Projection Tools)，對原始格式為HDF的MOD13Q1C6數(shù)據(jù)集進行NDVI波段提取、批量鑲嵌、重采樣和投影轉(zhuǎn)換，輸出數(shù)據(jù)格式為Geo-tiff的影像文件，投影為Alberts Equal Area。然后，基于最大化合成法(maximum value composition, MVC)得到月最大NDVI時間序列，該處理可以減少大氣中云、顆粒以及太陽高度角的影響。最后，利用均值法計算得到年平均NDVI時間序列，并裁剪出2001—2019年覆蓋研究區(qū)的NDVI時間序列。

1.2.2 DEM數(shù)據(jù)和氣象數(shù)據(jù) 本文所使用的DEM數(shù)據(jù)來源于美國國家航空航天局發(fā)布的SRTM DEM數(shù)據(jù)，空間分辨率30 m。為了與NDVI時間序列保持相同空間分辨率，將DEM重采樣為250 m空間分辨率。氣象數(shù)據(jù)來源于國家氣象科學技術(shù)中心提供的《中國地面氣候資料月值數(shù)據(jù)集》(http:∥data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_MON.html)，主要包括2000—2019年降水和氣溫數(shù)據(jù)，時間分辨率為月，氣溫的精度為0.1 ℃，降水量精度為0.1 mm，數(shù)據(jù)經(jīng)過嚴格的精度控制，質(zhì)量良好。為了提高研究結(jié)果的精度，本文選取覆蓋研究區(qū)及其周邊地區(qū)的65個氣象站點數(shù)據(jù)，以DEM作為協(xié)變量，基于ANUSPLIN插值模型生成降水和氣溫時間序列。

1.3 研究方法

1.3.1 Theil-Sen Median趨勢分析和Mann-Kendall顯著性檢驗 趨勢分析是一種非參數(shù)統(tǒng)計的斜率估計方法。由于該方法不需要數(shù)據(jù)服從特定的分布特征，且受異常值的影響較小，故常用于估算水文、氣象、植被指數(shù)等長時間序列數(shù)據(jù)的變化趨勢[15-18]。Mann-Kendall顯著性檢驗是一種科學的數(shù)學統(tǒng)計方法，可用于檢驗時間序列數(shù)據(jù)變化趨勢的顯著程度[19]。本文將兩種方法結(jié)合起來，用于估算植被NDVI變化趨勢，并根據(jù)變化斜率判斷其變化趨勢的顯著程度。

Theil-Sen Median趨勢分析用來計算植被NDVI的變化趨勢，用slope表示NDVI變化趨勢，當slope>0時，植被NDVI呈上升趨勢；當slope=0時，植被NDVI基本保持不變；當slope<0時，植被NDVI呈下降趨勢。Mann-Kendall顯著性檢驗法被用來檢測NDVI時間序列變化趨勢的顯著性，在給定顯著性水平α下，當|Z|>Z1-α/2，表示時間序列在α水平上變化顯著，反之，則變化不顯著。本文定義變化趨勢在α=0.05下顯著時，為顯著變化；在α=0.01下顯著時，為極顯著變化。

1.3.2R/S分析及Hurst指數(shù)R/S分析法(rescaled rang analysis method)，又稱之為重新標度極差分析法，最早是由Hurst在總結(jié)尼羅河的多年水文觀測資料時提出的一種分析方法[20]，后來經(jīng)過Mandelbrot和Wallis進一步補充和完善，將其發(fā)展成一種研究時間序列的分型理論[21]，目前在水文學、經(jīng)濟學、氣候?qū)W、地質(zhì)學等領(lǐng)域有著廣泛的應用。Hurst指數(shù)的取值范圍為(0

1.3.3 Pearson相關(guān)分析 Pearson相關(guān)分析可用來分析兩列數(shù)據(jù)之間是否存在相關(guān)關(guān)系及相關(guān)程度[22]，相關(guān)系數(shù)越大，代表兩個因子間的相關(guān)性越好，反之，則代表相關(guān)性越差。以往研究表明，植被NDVI對降水和氣溫變化具有時滯效應，且滯后期往往為0～3個月[6,23-24]，故本文基于MATLAB在像元尺度上計算植被NDVI與同期及前期1～3個月降水和氣溫的相關(guān)系數(shù)，并根據(jù)其最大相關(guān)系數(shù)，得到植被NDVI對降水和氣候變化響應的滯后期。

2 結(jié)果與分析

2.1 植被覆蓋時空演變特征

2.1.1 植被NDVI時間變化特征 如圖1所示，2001—2019年京津冀地區(qū)植被NDVI整體呈波動上升趨勢，上升速率為0.002 2/a。京津冀地區(qū)植被NDVI最大值為0.432，最小值為0.381，分別出現(xiàn)在2018年和2001年。北京、天津和河北植被NDVI變化趨勢表現(xiàn)出較大的差異，北京和河北植被NDVI呈上升趨勢，且北京植被NDVI上升趨勢明顯高于河北，而天津植被NDVI呈波動下降趨勢。以上結(jié)果表明，2001—2019年京津冀地區(qū)植被覆蓋整體呈現(xiàn)改善態(tài)勢，尤以北京改善程度最為顯著，但天津植被覆蓋呈現(xiàn)退化趨勢。

圖1 京津冀地區(qū)植被NDVI時間變化特征

2.1.2 植被NDVI空間變化特征 如圖2a和表1所示，2001—2019年京津冀地區(qū)植被NDVI變化趨勢呈現(xiàn)出明顯的空間異質(zhì)性。植被NDVI上升區(qū)域占總面積的76.19%，其中極顯著上升(p<0.01)約占總面積的49.31%，從京津冀地區(qū)東南部呈帶狀延伸至西北部，主要包括北京、保定和石家莊西北地區(qū)、張家口和承德南部地區(qū)、秦皇島北部以及滄州、邢臺和衡水東南部分地區(qū)。以上地區(qū)土地利用類型以林地和草地為主，地勢較高，受人類負向擾動較小。自1998年以來，得益于退耕還林還草工程、三北防護林工程和環(huán)北京和天津地區(qū)防沙治沙工程等一系列生態(tài)林業(yè)工程的實施，以上地區(qū)在植被群落恢復和生態(tài)系統(tǒng)建設(shè)等方面取得了顯著的成果，植被覆蓋呈明顯上升趨勢，荒漠化程度得到了有效的控制，荒漠化面積大大縮減[25-28]。

圖2 2001—2019年京津冀地區(qū)植被NDVI變化趨勢(a)及其未來變化趨勢(b)

植被NDVI下降區(qū)域占總面積的23.81%，其中，極顯著下降(p<0.01)和顯著下降(p<0.05)區(qū)域約占總面積的10.52%，主要分布在京津冀地區(qū)東南部各城市中心及其周邊地區(qū)，尤以唐山、天津、廊坊、石家莊和邯鄲為著。以上地區(qū)主要位于地勢平坦的華北平原，地勢較低，受社會經(jīng)濟發(fā)展、人類活動增強、人口密度增加等的負向影響，各城市中心及其周邊地區(qū)植被活動減弱、植被覆蓋呈下降趨勢。

表1 京津冀地區(qū)植被NDVI變化的顯著性統(tǒng)計

2.1.3 植被NDVI未來變化趨勢R/S分析常用來探究時間序列數(shù)據(jù)變化趨勢的持續(xù)性與反持續(xù)性。京津冀地區(qū)植被NDVI Hurst指數(shù)均值為0.507，整體呈弱持續(xù)性變化。植被NDVI呈持續(xù)性變化的區(qū)域約占京津冀地區(qū)總面積的50.18%，略大于呈反持續(xù)性變化的區(qū)域(49.82%)。

為揭示京津冀地區(qū)植被NDVI未來變化的持續(xù)性特征，本文將趨勢分析和R/S分析結(jié)果進行重分類并進行疊置分析，得到各像元未來植被NDVI變化趨勢，并將結(jié)果分成4個類別(詳見圖2b和表2)：持續(xù)下降、下降、上升、持續(xù)上升。統(tǒng)計分析表明，未來植被NDVI下降和持續(xù)下降區(qū)域占總面積的一半以上(53.11%)，其中下降的比例占總面積的39.83%，均勻分布于京津冀地區(qū)；呈持續(xù)下降的區(qū)域占總面積的13.28%，主要分布在環(huán)渤海灣地區(qū)的天津、廊坊和滄州及其周邊地區(qū)。未來植被NDVI呈持續(xù)上升的區(qū)域分布較為集中，主要分布張家口和承德大部分地區(qū)，其他城市零星分布，約占京津冀地區(qū)總面積的36.90%。

由以上結(jié)論可以推斷，由于多項林業(yè)生態(tài)工程的實施，構(gòu)建起了一張喬木、灌木和草地相結(jié)合的防護網(wǎng)，不僅起到涵養(yǎng)水源、保持水土的作用[29-30]，也為張家口和承德及其周邊地區(qū)植被生長和繁衍提供了良好的環(huán)境，使得該地區(qū)未來植被NDVI繼續(xù)以上升為主；而天津、廊坊和滄州及其周邊地區(qū)，由于受城市人口密度增加和城區(qū)面積擴張的影響，加速了耕地、林地和草地向城市建設(shè)用地的轉(zhuǎn)換[31]，使得未來植被覆蓋繼續(xù)保持退化態(tài)勢，影響區(qū)域生態(tài)環(huán)境保護和生態(tài)文明建設(shè)。

表2 京津冀地區(qū)植被NDVI未來變化趨勢統(tǒng)計

2.2 植被NDVI對降水和氣溫的響應

2.2.1 植被NDVI與降水相關(guān)關(guān)系 為探究京津冀地區(qū)植被生長對氣候變化的響應機制，計算2001—2019年京津冀地區(qū)植被NDVI與同期及前1～3個月降水和氣溫之間的相關(guān)系數(shù)，根據(jù)最大相關(guān)系數(shù)揭示京津冀地區(qū)植被NDVI對降水和氣溫變化最大響應的空間分布特征。在以上研究基礎(chǔ)上，利用T檢驗法對植被NDVI與降水和氣溫最大相關(guān)關(guān)系的顯著性進行檢驗，結(jié)果如圖3所示。

如圖3a所示，京津冀地區(qū)植被NDVI與降水的最大相關(guān)系數(shù)為-0.832至0.945，相關(guān)系數(shù)均值為0.319。植被NDVI對降水變化的響應強度呈現(xiàn)明顯的地域差異，主要表現(xiàn)為西北高，東南低。由表3可知，京津冀地區(qū)植被NDVI主要與降水呈正相關(guān)，相關(guān)系數(shù)為正的區(qū)域約占京津冀地區(qū)總面積的87.35%。

如圖3b所示，保定西北部、北京西部和北部、張家口東南部、滄州東部、以及承德西南部地區(qū)，植被NDVI與降水呈極顯著和顯著正相關(guān)，分別占總面積的17.80%和16.44%。以上地區(qū)由于受生態(tài)林業(yè)工程的影響，植被覆蓋程度呈逐年遞增的趨勢，且該地區(qū)屬半濕潤半干旱地區(qū)，水分是限制植被生長的主要氣象因子[10]。一方面，區(qū)域植被覆蓋狀況的改善能夠起到防止水土流失、涵養(yǎng)水源的作用。另一方面，水分條件的改善能夠促進區(qū)域植被的生長和繁衍進程。植被NDVI與降水呈極顯著和顯著負相關(guān)的區(qū)域約占總面積的0.63%，主要集中在廊坊和邯鄲各城市中心及其周邊地區(qū)。以上地區(qū)由于受高速經(jīng)濟發(fā)展和快速人口增長的影響，城區(qū)面積由中心向邊緣擴展，導致植被覆蓋呈下降趨勢，而該地區(qū)研究時段內(nèi)降水呈上升趨勢，故植被NDVI與降水呈負相關(guān)。

圖3 京津冀地區(qū)植被NDVI與降水相關(guān)系數(shù)(a)及顯著性(b)空間分布

2.2.2 植被NDVI與氣溫相關(guān)關(guān)系 如圖4a所示，京津冀地區(qū)植被NDVI與氣溫的最大相關(guān)系數(shù)為-0.804至0.906。植被NDVI對氣溫變化的最大響應強度由東南到西北依次呈現(xiàn)“弱—強—弱”的空間分布格局。由表3可知，京津冀地區(qū)氣溫對植被生長的正向促進作用遠小于降水，植被NDVI與氣溫的最大相關(guān)關(guān)系均值僅為0.132，由此可以推斷降水是京津冀地區(qū)限制植被生長的主要氣候因子，這與晏利斌等[10]研究結(jié)果一致，即京津冀地區(qū)植被NDVI與降水和氣溫變化均呈正相關(guān)，但與降水的相關(guān)性更加顯著。

如圖4b和表3所示，植被NDVI與氣溫最大相關(guān)系數(shù)為正的區(qū)域約占研究區(qū)總面積的69.14%。其中，植被NDVI與氣溫呈極顯著正相關(guān)和顯著正相關(guān)的區(qū)域約占總面積的9.90%，呈帶狀零星分布于京津冀地區(qū)腹地。植被NDVI與氣溫呈負相關(guān)的面積約占京津冀地區(qū)總面積的30.86%，主要分布在京津冀地區(qū)兩翼。其中，呈顯著負相關(guān)和極顯著負相關(guān)的面積約占總面積的3.23%，主要分布在京津冀地區(qū)東南沿線各城市中心及其周邊地區(qū)。由于受城市擴張和人類活動的影響，以上地區(qū)植被覆蓋狀況較差且抵御外界擾動的能力較弱，且城市周邊地區(qū)土地利用類型以耕地為主，受人類農(nóng)業(yè)管理活動的影響較大，京津冀地區(qū)部分區(qū)域植被NDVI對氣溫變化呈不顯著相關(guān)甚至負相關(guān)[19,32]，表明氣溫對植被生長具有較強的抑制作用。

表3 京津冀地區(qū)植被NDVI與降水和氣溫最大相關(guān)系數(shù)顯著性檢驗統(tǒng)計

圖4 京津冀地區(qū)植被NDVI與氣溫相關(guān)系數(shù)(a)及顯著性(b)空間分布

2.2.3 植被NDVI對降水和氣溫變化最大響應的滯后效應 考慮到植被生長對降水和氣溫變化具有時滯效應，本文通過植被NDVI與降水和氣溫的最大相關(guān)系數(shù)，得到京津冀地區(qū)植被NDVI對降水和氣候變化最大響應滯后期，揭示植被生長對降水和氣溫變化時滯效應的空間分布特征。

如圖5a和表4所示，京津冀地區(qū)植被生長對降水變化存在明顯滯后效應，且滯后期表現(xiàn)出明顯的地域分異格局。植被NDVI對前2月和前3月降水變化最大響應的區(qū)域約占總面積的58.74%；其中，植被NDVI對前3月降水最大響應的面積最大，約為44.22%，主要分布在張家口西北部、承德東北部、唐山東部、滄州、邯鄲、石家莊西部和保定中部；植被NDVI對前1月降水最大響應的區(qū)域約占京津冀地區(qū)總面積的27.33%，僅次于前3月，呈片狀分布在張家口東南部、承德西南部和北京。

如圖5b和表4所示，相較于降水，京津冀地區(qū)植被生長對氣溫變化的滯后期較短。植被NDVI對當月和前1月氣溫變化最大響應的區(qū)域約占總面積的68.41%；其中，植被NDVI對當月氣溫變化最大響應的區(qū)域約占總面積的32.81%，呈片狀分布在張家口及其周邊地區(qū)；植被NDVI對前1月氣溫變化最大響應的面積最大，約為35.60%，呈帶狀從京津冀地區(qū)東北部延伸至西南部。植被NDVI對前3月氣溫變化最大響應區(qū)域僅占總面積的15.23%，主要集中分布在保定、石家莊和秦皇島西北部。

圖5 京津冀地區(qū)植被NDVI對降水(a)和氣溫(b)變化最大響應滯后期空間分布

由以上可知，京津冀地區(qū)植被NDVI對降水和氣溫變化均存在明顯的時滯效應，植被NDVI對降水變化的滯后期略長于氣溫，這一結(jié)論與前人的研究結(jié)果一致[23-24]。降水過程中的部分降水經(jīng)地表徑流流出，不能通過滲透作用進入土壤；而另一部分可被土壤吸收，補償土壤水分，改善土壤墑情。只有保留在土壤部分的降水，可以被認為是植被生長的“有效降水”。降水從地表經(jīng)過滲透作用被土壤吸收，而后通過根部吸收轉(zhuǎn)移到植物的各個組織，最后成為植物生長可用的水資源需要花費一定的時間[33-34]。因此，植物生長往往對降水變化具有一定的滯后效應。相較于降水，溫度的變化可以直接影響植物的呼吸作用強度和光合作用效率，對植被生長具有更加直接的作用，故植被生長對氣溫變化的滯后期更短。

表4 京津冀地區(qū)植被NDVI對降水和氣溫變化最大響應滯后期統(tǒng)計

3 結(jié) 論

本文利用多源數(shù)據(jù)和多數(shù)學分析模型，在不同時空尺度上定量反演了近19 a京津冀地區(qū)植被覆蓋時空演變特征，并預測其未來變化趨勢，進一步探究了植被NDVI對降水和氣溫變化的時空響應特征及時滯效應。主要研究結(jié)論如下：

(1) 得益于林業(yè)生態(tài)工程的實施，京津冀地區(qū)植被覆蓋整體呈現(xiàn)改善態(tài)勢。近19 a來京津冀地區(qū)植被NDVI上升區(qū)域約占總面積的76.19%。其中，極顯著上升約占總面積的49.31%，主要集中在京津冀地區(qū)西北部。

(2) 京津冀地區(qū)未來植被呈下降趨勢的面積略大于上升趨勢的面積，其中呈持續(xù)下降的區(qū)域約占總面積的13.28%，主要分布在天津、廊坊和滄州及其周邊地區(qū)，以上地區(qū)受社會經(jīng)濟發(fā)展、城市擴張和人口密度增加的負向影響較大。

(3) 植被NDVI對降水和氣溫變化的響應強度呈現(xiàn)明顯的地域差異，植被NDVI與降水的相關(guān)性主要表現(xiàn)為西北高，東南低，而植被NDVI與氣溫的相關(guān)性由東南到西北依次呈現(xiàn)“弱—強—弱”的空間分布格局。植被NDVI與降水的相關(guān)關(guān)系略強于氣溫。因此，降水被認為是京津冀地區(qū)控制植被生長的主要氣候因子。

(4) 京津冀地區(qū)植被生長對降水的氣溫的變化存在明顯的時滯效應。植被NDVI對前3月降水最大響應的面積最大，約占總面積的44.22%；植被NDVI對前1月氣溫變化最大響應的面積最大，約占總面積的35.60%。