京津冀地區植被NDVI動態變化及其與氣候因子的關系

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

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

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

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

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

1 研究區概況、數據來源及研究方法

1.1 研究區概況

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

1.2 數據來源

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

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

1.3 研究方法

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

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指數R/S分析法(rescaled rang analysis method)，又稱之為重新標度極差分析法，最早是由Hurst在總結尼羅河的多年水文觀測資料時提出的一種分析方法[20]，后來經過Mandelbrot和Wallis進一步補充和完善，將其發展成一種研究時間序列的分型理論[21]，目前在水文學、經濟學、氣候學、地質學等領域有著廣泛的應用。Hurst指數的取值范圍為(0

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

2 結果與分析

2.1 植被覆蓋時空演變特征

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

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

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

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

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

表1 京津冀地區植被NDVI變化的顯著性統計

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

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

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

表2 京津冀地區植被NDVI未來變化趨勢統計

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

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

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

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

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

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

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

表3 京津冀地區植被NDVI與降水和氣溫最大相關系數顯著性檢驗統計

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

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

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

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

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

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

表4 京津冀地區植被NDVI對降水和氣溫變化最大響應滯后期統計

3 結 論

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

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

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

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

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