基于延安市不同植被類型NDVI對氣象因子影響及空間異質性

張家政，閔志強，王得軍，李清順，孫景梅，李宏韜，李崇貴*

(1.西安科技大學 測繪科學與技術學院，陜西 西安 710054；2.國家林業和草原局西北調查規劃設計院 旱區生態水文與災害防治國家林業局重點實驗室，陜西 西安 710048)

植被與氣候變化響應研究是目前生態與環境研究中的重要內容之一[1-3]，對土壤系統與大氣系統的能量傳輸具有至關重要的作用[4-5]。目前，研究方法更多是基于植被指數(normalized difference vegetation index，NDVI)進行不同尺度分析。如白旭陽等[6]通過采用對不同驅動因子進行貢獻率量化的方法分析1982-2015年氣候和人類活動對于新疆瑪納斯流域植被變化的影響，氣候變化對于NDVI變化的貢獻率高于人類活動對于NDVI變化的貢獻率，并且氣溫和降水對NDVI存在顯著性的影響；M.Lamchin等[7]通過分析1982-2014年亞洲地區氣象因子(溫度、降水和蒸散發等)與植被指數NDVI相關關系，結果表明溫度是導致亞洲地區植被指數NDVI變化最主要的因素；許玉鳳等[8]從年、季、月等不同尺度研究NDVI動態變化趨勢以及對氣候因子(日照時數、氣溫、降水、地表溫度等數據)響應，結果顯示貴州植被生長是受溫度調控，且氣溫升高、氣候變暖是促使植被生長季延長的氣候效應之一。很多研究將全區NDVI均值與氣象因子做相關性分析[9-10]，有些研究將NDVI在時間尺度作為一個連續變量，與氣象因子做空間相關性分析[7,11]。前者只考慮到了相關性關系，但是忽略了NDVI空間變化對于氣象因子響應；后者也沒有分析不同植被類型與氣象因子的響應。因此，為了解不同植被類型NDVI對于氣象因子響應，本研究基于長時序Landsat數據集與植被覆蓋數據，探索延安市不同植被類型時空變化以及與氣象因子(平均溫度、總降水、最高溫度和最低溫度)的響應，并采用最小二乘法+時空地理加權回歸模型[12-13]揭示NDVI與氣象因子擬合的空間異質性。

Google Earth Engine(GEE)平臺不僅提供存儲海量的遙感數據，而且廣泛用于變化監測[14-15]、土地利用動態監測[16]、農業應用[17]和災害監測與評估[18]等。GEE提供了JavaScript 和Python的API和基于Web的交互式開發環境，可以輕松地利用Google云的強大空間地理分析功能，使得長時序植被動態分析成為可能。

1 材料與方法

1.1 研究區概況

延安市(107°41′-110°31′E，35°21′-37°31′N)位于我國黃河流域中部，屬于黃土丘陵溝壑區，是黃土高原乃至全國水土流失最嚴重的地區之一。該市總面積為3.67萬km2，平均海拔1 200 m。延安地區夏天炎熱少雨，冬季干燥寒冷，年均氣溫為9～12℃，年均降水量在100～800 mm，氣溫年、日較差大。為遏制嚴重的水土流失，1999年國家開始實施生態修復工程。

1.2 研究方法

1.2.1 數據來源 Landsat影像 本研究Landsat遙感影像來自于GEE平臺(https://earthengine.google.com/)上landsat surface reflectance data(陸地衛星地表反射率數據)。該數據集已經經過大氣輻射校正處理，其空間分辨率為30 m，時間分辨率為16 d。首先在GEE平臺上利用空間和時間獲取1993-2017年延安市Landsat數據集[1993-2011年Landsat 5 TM(‘LANDSAT/LT05/C01/T1_SR’)、2012年Landsat 7 ETM+(‘LANDSAT/LE07/C01/T1_SR’)和2013-2017年Landsat 8 OLI(‘LANDSAT/LT08/C01/T1_SR’)]；然后批量對每一景影像進行去云掩膜處理和裁剪處理，并計算歸一化植被指數(NDVI)；最后采用最大值合成法(maximum value composition,MVC)對每年NDVI影像進行最大值合成，以此消除大氣和云等干擾誤差，得到的年最大NDVImax代表相應各年的NDVI。

氣象數據 降水(PRE)、平均氣溫(Tmean)、最高溫度(Tmax)和最低溫度(Tmin)數據來源于國家科技基礎條件平臺——國家地球系統科學數據中心(http://www.geodata.cn)，其空間分辨率為1 km，選取時間序列從1993年1月至2017年12月逐月數據。因數據格式為.nc，所以先利用Python進行數據預處理，使之變成可以使用的柵格圖像，再將其合并成年尺度數據[19]。

土地覆蓋數據 土地覆蓋數據來源于歐洲航天局的1993-2017年共25期土地利用/覆蓋(ESA CCI-LC)數據集(http://maps.elie.ucl.ac.be/CCI/viewer.)，空間分辨率300 m。延安市植被覆被類型分為針葉林(coniferous forest)、闊葉林(broad-leaved forest)、灌木/草地(shrub/grassland)、鑲嵌林地(mosaic tree)(喬木分布>50%)、鑲嵌草地(mosaic grassland)(草本分布>50%)、人類占用地(human occupied land)和水體(water bodies)共7類(圖1)。

圖1 研究區植被覆蓋情況

1.2.2 數據處理方法 線性分析研究表明[20-21]，一元線性回歸模型(Slope)能真實有效地反映植被的年際變化。因此，本研究采用一元線性回歸對延安市每個像元1993-2017年植被NDVImax進行趨勢分析，計算公式為：

(1)

Mann-Kendall(M-K)檢驗 與以往研究者所采用的基于最小二乘法的時間序列線性回歸模型的檢驗方法相比較，M-K屬于非參數檢驗方法，該方法不需要樣本遵循一定的分布，也不受少數異常值的干擾[22]。M-K檢驗已經在氣象、水文等領域的時間序列分析中得到廣泛的應用。

重心遷移模型 為深入研究各植被類型時空集聚和遷移特征，本研究采用重心遷移模型進行植被時空演化的定量表達，以此揭示延安市不同區域植被類型活動的演變軌跡。其計算公式如下：

(2)

最小二乘回歸模型+時空地理加權回歸模型(OLS+GTWR)最小二乘法(ordinary least squares，OLS)常常在遙感變化監測中應用，也是常用的求解無約束最優化問題的回歸分析方法之一。具體表達式為：

(3)

式中,y代表因變量，本研究指的是NDVImax；β0為常數項，βj為第j個自變量的回歸系數；xij代表第i個樣本的第j個自變量；ζ為符合正態分布的隨機誤差項。

地理加權回歸模型(GWR)是空間異質性研究的經典模型，而時空地理加權回歸模型(GTWR)則是將時間屬性聯系到GWR的空間屬性中，以更好地反映研究區時空變化信息，使得模型估算結果更加有效[23-24]。具體表達式為：

(4)

式中,yi代表第i個樣本點的被解釋變量；ui，vi分別為第i個樣本點的經、緯度坐標；ti為第i個樣本點的時間坐標；β0(ui,vi,ti)為第i個樣本點的常數項；βk(ui,vi,ti)為第k個解釋變量在第i個樣本點的回歸系數。

1.2.2.3 獨立生活能力訓練:每周1次，每次40 min。向患者講解良好的生活習慣的重要性，對患者日常生活能力進行基礎評估，與患者共同制定個體化日常行為規范，督促其按時睡覺、按時起床、整理病床、洗漱、自動排隊進食服藥及定時參加適當的體育鍛煉，組織患者看電視、報刊和雜志。幫助患者建立良好的生活習慣，培養患者有規律的生活和自理能力。

本研究通過借助ArcGIS10.3平臺Zonal Statistics工具提取每種植被類型的年NDVImax、PRE、Tmean、Tmax和Tmin等數據，再采用OLS+GTWR模型的方法建立不同土地覆蓋類型年NDVImax與各種氣象因子的回歸模型。

2 結果與分析

2.1 NDVImax和PRE、Tmean、Tmax和Tmin年際變化特征

利用一元線性回歸方法分析延安市NDVImax和PRE、Tmean、Tmax和Tmin的年際變化趨勢(圖2)。NDVImax整體呈現波動增長趨勢，NDVImax均值由1993年的0.239增加到2017年的0.386，并在2013年達到最大值0.511，之后有所降低。NDVImax變化率為0.008/a(P<0.01)，植被生長整體呈現較好的生長趨勢；PRE、Tmax和Tmin總體呈現波動增長趨勢，其增長速率分別為4.619 mm/a(P<0.05)，0.035℃/a(P<0.01)和0.013℃/a，其中波動最大的是Tmax，由1993年的14.01℃增長到2017年15.49℃；Tmean總體則呈現平穩增長趨勢，其增長率為0.023℃/a。

圖2 1993-2017年延安市NDVImax和PRE、Tmean、Tmax和Tmin年際變化趨勢

2.2 土地覆蓋面積及土地覆蓋變化對植被指數NDVImax的影響

除水體在25 a基本保持不變(20.146 km2)以外，其他土地覆蓋類型面積均發生了不同程度的變化(圖3a)。針葉林和闊葉林面積均呈現增加趨勢，其變化率分別為0.185 km2/a(P<0.01)和0.007 km2/a(P<0.01)。特別的，闊葉林面積變化最大(增加3.85×103km2)；灌木/草地的面積呈現先減少后增加趨勢，1993-2007年草地/灌木呈現減少趨勢，變化率為-5.279/a(P<0.01)，而在2008年急劇增加為12.345×103km2；鑲嵌林地和鑲嵌草地面積均呈現減少趨勢，其變化率分別為-0.014 km2/a(P<0.01)和-7.223 km2/a(P<0.01)；人類占用地呈現總體增加趨勢[0.001 km2/a(P<0.01)]。特別的，在2006年人類占用地面積達到最大值(0.989×103km2)，而在2008年面積由2007的0.989×103km2急劇減少為0.986×103km2，其后保持穩定增加狀態。不同土地覆蓋類型NDVImax變化介于0.1～0.7，其NDVImax大小分別為針葉林>闊葉林>鑲嵌林地>鑲嵌草地>灌木/草地(圖3b)，且逐年呈現增加趨勢。其中針葉林增加幅度最大，分別為0.009 1/a；闊葉林增加幅度最小，為0.007 1/a；灌木/草地、鑲嵌林地和鑲嵌草地增加幅度一致，為0.008/a。

圖3 不同植被覆蓋的面積以及NDVImax年變化

2.3 不同植被覆蓋類型重心演變軌跡

重心分布變化能夠從空間上反映不同植被覆蓋類型的時空演變過程。為進一步探討延安市植被土地覆蓋類型重心的空間變化規律，本研究分析了1993-2017年植被覆蓋類型重心演變情況(圖4)。針葉林重心主要集中在洛川縣北部地區，其重心總體向西南部偏移2 516.03 m，其中，2004-2005年重心轉移距離最大(776.55 m)；闊葉林重心主要集中在富縣西部地區，其重心變化幅度最小，重心總體向東偏移675.96 m，其中，1994-1995年重心轉移距離最大(186.40 m)；灌木/草地重心主要集中在安塞縣中南部地區，其重心總體向東南偏移1 186.70 m，其中，1994-1995年重心轉移距離最大(331.29 m)；鑲嵌林地重心主要集中在富縣西部地區，其重心變化幅度最大、最復雜，重心總體向南偏移4 157.74 m，其中，2002-2003年重心轉移距離最大(1 240.62 m)；鑲嵌草地重心主要集中在甘泉縣東部地區，其重心總體向西南偏移1 896.36 m，其中2007-2008年重心轉移距離最大(1 235.66 m)。

圖4 重心轉移

2.4 植被指數NDVImax空間變化特征

基于Slope+M-K趨勢檢驗的分析方法得到NDVImax變化趨勢的空間分布(圖5a)。NDVImax以增加為主，面積占研究區總面積的80.40%[(顯著增加(Slope≥0.02)和增加(0.002≤Slope<0.02)的總和)]，遠大于NDVImax減少趨勢面積[(顯著減少(Slope≤-0.02)和減少(-0.02≤Slope<-0.002)的總和)]，表明植被覆蓋狀況呈現不斷改善趨勢。NDVImax呈現減少趨勢地區主要集中在延安市中部地區(寶塔區南部、甘泉縣北部和富縣中北部等)。該地區主要以灌木/草地為主，隨著城鎮化建設的加快推進，部分地區灌木/草地轉化為人類占用地加劇了植被退化。顯著增加區域主要分布在黃龍縣東部、宜川縣中北部和延長縣南部等區域。

通過對NDVImax進行顯著性檢驗(圖5b)，1993-2017年，除1993年外，NDVImax在其余時間段的UF曲線均>0，呈上升趨勢，且2004-2017年UF曲線超過95%顯著水平線，說明此時間段延安市NDVImax增長趨勢顯著，植被改善明顯。

圖5 年NDVImax變化趨勢及M-K檢驗

2.5 全區植被指數NDVI對氣象因子的響應

本研究基于像元對植被年NDVImax與各氣象要素進行相關性分析(圖6)。圖6a的均值為0.262，從像元尺度來看，研究區NDVImax與PRE呈正相關的面積89.53%，主要分布在安塞縣、延川縣、延長縣以及宜昌縣北部地區，尤其在延川縣延川鎮和關莊鎮、安塞縣張坪鄉等地NDVImax與PRE之間相關性較高，呈極顯著正相關(P<0.01)，PRE是這些地區植被生長的主要限制因子，其他地區植被與降水相關性不明顯；圖6b的均值為0.085，從像元尺度來看，研究區NDVImax與Tmean呈正相關的面積71.38%，主要分布在延川縣中部、宜川縣南部、黃龍縣中部以及洛川縣中部等地區。研究區NDVImax與Tmean呈顯著負相關(P<0.01)的地區占總面積的6.02%，主要分布在志丹縣南部橋鎮和永寧鎮等地區，Tmean是影響該地區植被生長的主要限制因子；圖6c的均值為0.119，研究區NDVImax與Tmin呈正相關的面積75.57%，主要分布在子長縣中南部、延川縣西部和延長縣中北部等地區，尤其在延川縣賈家坪鎮、子長縣王家灣鎮等地NDVImax與Tmin之間相關性較高，呈極顯著正相關(P<0.01)，占研究區總面積的1.61%；圖6d的均值為0.033，研究區NDVImax與Tmax呈正相關的面積57.02%，主要分布在洛川縣中部、富縣中南部以及黃龍縣。研究區NDVImax與Tmax呈極顯著負相關(P<0.01)的地區占總面積的3.02%，主要分布在志丹縣南部等地區。

圖6 1993-2017年NDVImax與PRE(a)、Tmean(b)、Tmin(c)和Tmax(d)相關系數空間分布

有研究表明，水熱條件的增加會促進植被生長，且氣候的變化對植被生長影響存在一定的滯后性[23]，這些因素都有可能造成氣候變化對植被無顯著性影響，也有可能是不同植被覆蓋對氣候變化存在差異影響，從而掩飾了氣候變化對于植被的作用。

2.6 不同植被類型NDVImax對氣象因子的響應機制

在不同植被覆蓋類型區氣象因子對各種植被覆蓋類型的作用有明顯差異。由表1知，針葉林與Tmean和Tmin顯著性相關，年NDVImax與Tmean和Tmin呈正相關，相關系數分別為0.431(P<0.05)和0.505(P<0.05)；闊葉林僅與Tmin顯著性相關，年NDVImax與Tmin呈正相關，相關系數分別為0.396(P<0.05)；灌木/草地和人類占用地與PRE和Tmin顯著性相關，年NDVImax與PRE和Tmin呈正相關，相關系數分別為0.402(P<0.05)、0.409(P<0.05)和0.397(P<0.05)、0.408(P<0.05)；鑲嵌林地和鑲嵌草地僅與Tmin顯著性相關，年NDVImax與Tmin呈正相關，相關系數分別為0.433(P<0.05)和0.414(P<0.05)；水體與各氣象因子均無顯著性關系，Tmax與各種覆蓋類型年NDVImax均無顯著相關性，說明Tmax在年尺度上對不同類型植被作用均不顯著。

表1 不同植被類型NDVImax與氣象因子相關系數

對不同植被覆蓋類型年NDVImax與氣象因子建立OLS模型，剔除方差膨脹因子(Variance Inflation Factor，VIF)>7.5的氣象因子，然后采用不同植被類型NDVImax的變化量作為因變量，剔除完VIF>7.5的因子作為自變量，對NDVImax的變化進行GTWR分析(表2)。相比OLS模型，GTWR模型對于各植被類型與氣象因子的擬合效果顯著提高，不同植被類型的GTWR模型擬合平均R2(0.419)遠大于OLS模型的擬合平均R2(0.100)，與OLS模型的AICc之差(ΔAICc=5254.78)遠>3，二者擬合差異明顯，GTWR模型擬合度高于OLS模型。

表2 GTWR與OLS模型結果比較

對GTWR模型各因子擬合系數采用箱線圖描述(圖7)。針葉林箱線圖中，PRE、Tmean、Tmax和Tmin對NDVImax平均回歸系數分別為0.040、-0.017、0.030和0.075，表明PRE、Tmax和Tmin對針葉林NDVImax呈正向影響，Tmean對針葉林NDVImax呈負向影響；上下四分位數范圍分別為-0.053～0.126、-0.152～0.081、-0.076～0.135和0.024～0.136，表明各年份PRE和Tmin密度系數變化幅度相對較大，而Tmean和Tmax密度系數變化差異不顯著；闊葉林箱線圖中，PRE、Tmean、Tmax和Tmin對NDVImax平均回歸系數分別為-0.037，-0.008，-0.127和0.112，表明PRE、Tmean和Tmax對闊葉林NDVImax呈正向影響，Tmin對闊葉林NDVImax呈負向影響；上下四分位數范圍分別為-0.093～0.077、-0.067～0.038、-0.185～-0.076和0.065～0.158，表明各年份PRE和Tmin密度系數變化幅度相對較大，而Tmean和Tmax密度系數變化差異不顯著；灌木/草地箱線圖中，PRE、Tmax和Tmin對NDVImax平均回歸系數分別為0.087、0.011和0.030，中位數分別為0.075、0.013和0.033，表明PRE、Tmax和Tmin對灌木/草地NDVImax呈正向影響；上下四分位數范圍分別為0.06～0.168、-0.074～0.094和-0.038～0.092，表明各年份PRE密度系數變化幅度相對較大，而Tmax和Tmin密度系數變化差異不顯著；鑲嵌林地箱線圖中，PRE、Tmean、Tmax和Tmin對NDVImax平均回歸系數分別為-0.017、-0.034、0.050和0.112，表明PRE和Tmean對鑲嵌林地NDVImax呈正向影響，Tmax和Tmin對鑲嵌林地NDVImax呈負向影響；上下四分位數范圍分別為-0.035～0.063、-0.038～-0.030、0.038～0.06和0.111～0.113，表明各年份PRE、Tmean、Tmax和Tmin密度系數變化差異不顯著；鑲嵌草地箱線圖中，PRE、Tmax和Tmin對NDVImax平均回歸系數分別為0.061、-0.094和0.097，表明PRE和Tmin對鑲嵌草地NDVImax呈正向影響，Tmax對鑲嵌草地NDVImax呈負向影響；上下四分位數范圍分別為-0.008～0.117、-0.216～-0.004和0.025～0.167，表明各年份Tmax密度系數變化幅度相對較大，而PRE和Tmin密度系數變化差異不顯著。

圖7 GTWR擬合回歸系數箱型

3 結論與討論

隨著國家對延安市實行大規模退耕還林等生態修復工作以來，延安市植被動態變化監測已成為目前研究的熱點[25-26]，但目前很少研究不同植被覆蓋類型變化及其與氣象因子關系，由于單一植被生態系統在水土保持、調節氣候等方面與多植被類型存在很大差距，所以研究不同植被覆蓋類型變化具有較深的意義。對此，本研究定量探究了不同植被覆蓋類型變化及其與氣象因子關系，揭示了不同植被類型對氣象因子的響應規律。與MODIS數據相比，Landsat衛星具有時間跨度長、空間分辨率高(30 m)的優勢，能夠為長時序植被動態監測提供更好的數據支撐。GEE平臺可以提供大量的數據，這樣為今后長時序、大幅度研究提供了便捷。

本研究基于ESA CCI-LC植被覆蓋分類數據以及Landsat數據集，使用Slope+Manna-Kendall非參數檢驗方法，分析全區和不同植被覆蓋類型區內NDVImax與氣象因子時間變化特征與變化趨勢，并應用重心遷移模型分析其空間變化特征，最后應用最小二乘法和時空地理加權回歸模型進一步揭示不同植被覆蓋類型區內NDVImax分異特征及其驅動力，得到以下結論。

1)1993-2017年，延安市全區植被年NDVImax整體呈現波動增長趨勢，NDVImax均值由1993年的0.239增加到2017年的0.386，針葉林、闊葉林、灌木/草地面積呈現增加趨勢，而鑲嵌林地和鑲嵌草地面積則呈現減少趨勢。不同植被覆蓋類型年NDVImax也呈現增長趨勢，其中針葉林增加幅度最大，為0.009 1/a。這與朱會利等[25]、何立恒等[26]的研究結果相一致，這表明自延安市生態修復工程以來，植被覆蓋狀況極大改善。

2)延安市全區NDVI與氣象因子(Tmean、PRE、Tmax和Tmin)之間沒有顯著相關性，但在不同植被覆蓋類型情況下，氣象因子對NDVImax存在顯著性作用。這說明對植被類型進行分類更加有利于揭示氣象因子對于植被的作用機制。其中，針葉林年NDVImax與Tmean和Tmin呈現顯著性正相關，相關系數分別為0.431(P<0.05)和0.505(P<0.05)；闊葉林年NDVImax僅與Tmin呈現顯著性正相關，相關系數分別為0.396(P<0.05)；灌木/草地和人類占用地年NDVImax與PRE和Tmin呈現顯著性正相關，相關系數分別為0.402(P<0.05)、0.409(P<0.05)和0.397(P<0.05)、0.408(P<0.05)；鑲嵌林地和鑲嵌草地年NDVImax僅與Tmin呈現顯著性相關，相關系數分別為0.433(P<0.05)和0.414(P<0.05)。

3)針葉林主要分布在洛川縣，遷移速率分別為100.64 m/a；灌木/草地重心主要分布在安塞縣，遷移速率為47.46 m/a；鑲嵌草地主要分布在甘泉縣，遷移速率為75.85 m/a；鑲嵌林地和闊葉林主要分布在富縣，遷移速率為分別為166.31 m/a和27.04 m/a。

4)Tmean對闊葉林和鑲嵌林地NDVImax呈現正向影響，對針葉林NDVImax呈現負向影響；PRE對各種植被類型NDVImax均呈現正向影響；Tmax對針葉林、闊葉林和灌木/草地NDVImax呈現正向影響，對鑲嵌草地和鑲嵌林地NDVImax呈現負向影響；Tmin對針葉林、灌木/草地和鑲嵌草地NDVImax呈現正向影響，對闊葉林和鑲嵌林地NDVImax呈現負向影響。