植被類型對樹輪重建賀蘭山NDVI的影響

李明明，李剛

(1.中國科學院西北生態環境資源研究院，沙漠與沙漠化重點實驗室，甘肅 蘭州 730000；2，中國科學院大學，北京 100049；3.甘肅祁連山國家級自然保護管理局，甘肅 張掖 734000)

樹木作為組成森林中獨立的生命有機體，對氣候變化極為敏感且為減緩氣候變化帶來的影響發揮了重要作用，在其生長過程中樹木年輪記錄了大量的氣候、環境等方面的信息[1-2].然而，樹木的生長又分為初級生長和次級生長[3-4]，初級生長包括發芽、枝葉的生長及根系的衍生，次級生長包括樹干、枝及根系的生長.一方面樹木冠層的葉片通過光合作用合成有機質并通過樹木的韌皮部輸送到樹干以及樹根的各部位，另一方面，樹干將根系吸收的水和無機鹽等營養物質通過邊材運送到樹冠，并與根系一起支撐整棵樹木[5].因此，研究植被冠層與主干間的生長關系，有利于我們明確氣候變化背景下樹木生長的動態變化，更好地了解樹木的生理機制以及為應對氣候變化提供參考依據.

隨著遙感技術的進步與發展，利用遙感衛星監測探究植被的地表覆蓋度、生物量、凈初級生產力變化的研究日益增多[6-8]，也為研究氣候變化對植被動態的影響提供重要參考.歸一化植被指數(normalized difference vegetation index，NDVI)作為使用較普遍的遙感數據之一，它是植物生長狀態及植被空間分布密度的最佳指示因子，能反映冠層繁茂程度和植被覆蓋度的動態變化[9].但NDVI僅代表樹木冠層的變化情況，而生長過程中樹木冠層與主干之間存在差異，且由于遙感衛星數據記錄較短，限制了我們對過去植被變化歷史的認識與理解.利用樹輪資料對NDVI進行重建則可有效解決這一問題，因為樹輪資料具有定年準確、分辨率高、樣本量大等優點[10]，被廣泛用于重建過去不同時期的溫度、降水、干濕指數、徑流量等環境信息[11-15]，同樣可以利用樹輪資料重建NDVI，不僅可以彌補NDVI時間序列長度上的不足，且對揭示年代至百年尺度的植被生長變化規律、預估未來氣候情景下植被生長變化具有重要意義.

國內外學者已經開展了NDVI與樹輪資料的相關研究[16-22]，汪青春等[16]利用樹輪資料重建了青海布哈河流域的NDVI，王文志等[17]利用樹輪資料重建了祁連山地區的NDVI.此外，Kaufmann等[18]發現在北半球中高緯度地區兩者的相關性受不同緯度和樹種影響；D’Arrigo等[19]研究發現美國北部地區樹輪最大晚材密度與NDVI呈顯著相關，然而在北美高緯度地區、歐洲地區和加拿大的部分地區并沒有發現顯著相關性[20-22]，表明樹輪資料和NDVI之間的響應關系可能因研究區域和樹種或植被類型不同有所差異.以上研究在不同空間尺度上證實了NDVI與樹輪資料存在相關性，但結果存在差異.NDVI重建過程有諸多不確定性，可能是因為NDVI包含區域內所有植被的生長信息，不同植被類型的生理過程迥異，而現有的研究中幾乎均未考慮NDVI中植被類型對重建結果的影響[16-22]，因此直接將NDVI結果與樹輪資料結合進行重建，可能導致重建的NDVI序列出現偏差.如何提高重建結果的可靠性，使重建的NDVI更加準確反映植被變化是一個值得探討的問題.本研究選取賀蘭山為研究區域，嘗試將NDVI中林地與草地信息分離，并單獨使用林地NDVI與樹木徑向生長的關系對該地區NDVI進行重建，以期更準確地了解賀蘭山地區林地NDVI的變化歷史，為預測未來森林動態提供理論.

1 材料與方法

1.1 試驗區概況

賀蘭山(N 38°27′～39°30′、E 105°20′～106°41′)位于寧夏回族自治區與內蒙古自治區交接地帶(圖1)，地處東亞夏季風的西北緣，是我國干旱與半干旱、草原與荒漠的過渡地帶，生態環境敏感、脆弱[23].賀蘭山地區為典型的溫帶大陸性氣候，優勢植被主要為針葉林，分布在海拔2 000 m以上，2 000 m以下主要為草地[24].根據氣象站觀測資料顯示(圖3)，賀蘭山全年干旱少雨，多年平均氣溫為9.2 ℃，最冷月(1月)平均氣溫為-8 ℃；最熱月(7月)平均氣溫為23.8 ℃；年降水量在250～438 mm，主要集中在5～9月，約占全年降水量的80%左右，年蒸發量在2 000 mm以上[25].由于年降水量與蒸發量相差較大，空氣干燥；日照強烈，年平均日照達3 000 h以上，無霜期短.

圖1 賀蘭山土地覆被及樹輪采樣點Figure 1 Land cover and tree-ring sampling site in Helan Mountains

1.2 樣本采集與年表建立

本研究油松樣本采集自賀蘭山蘇峪口國家森林公園，采樣點選取受人為干擾較小且遠離旅游區的原生林，海拔在2 100～2 300 m.共選取98棵油松樣樹，分別用生長錐從垂直和平行山坡的方向上在樹高1.3 m處采集2根樹芯.對采集的樣芯按照國際通用的試驗處理步驟進行樣品前處理[27]：將樣芯晾干、固定、打磨至顯微鏡下樹木年輪輪廓清晰可見，并剔除斷裂及破損的部分樣芯，共篩選得到183根樣芯用于本研究.使用LINTAB-6系統對樹木年輪寬度進行測量，精度為0.01 mm.之后利用TSAPwin程序完成交叉定年[28]，通過COFECHA程序對交叉定年結果進行檢驗.對完成定年后的樹木年輪寬度數據使用ARSTAN軟件來去除與樹木年齡有關的生長趨勢[29]，并基于負指數函數對數據進行擬合，得到樹輪寬度標準年表(圖2).選取1920～2018年(EPS>0.85)的賀蘭山油松樹輪寬度年表用于后續分析.

圖2 賀蘭山地區樹輪寬度標準年表及樹芯樣本量Figure 2 Tree-ring standard chronology for Helan Mountains and sample size expressed as number of cores

1.3 氣象數據

氣象數據下載自中國氣象科學數據共享服務網(http：//cdc.nmic.cn /)，所選臺站為銀川氣象站，其距離采樣點50 km，可用時段1951～2018年.選取該站的平均氣溫、平均最高氣溫、平均最低氣溫、降水量、相對濕度等要素用于本研究.

帕爾默干旱指數(palmer drought severity index，PDSI)用于分析溫度和降水的綜合影響，數據下載自荷蘭皇家氣象研究所(KNMI)(http://climexp.knmi.nl)，格點位置：N 38°42′，E 105°45′，時段1934～2018年.

1.4 土地覆被數據

本文選用清華大學研發的全球土地覆蓋產品[26](FROM-GLC10)，下載網址：http://data.ess.tsinghua.edu.cn/fromglc2017v1.html.該數據通過對覆蓋全球的陸地衛星專題制圖儀和增強型專題制圖儀傳感器(Landsat TM/ETM+)影像進行分類，獲得10 m分辨率的FROM-GLC10，該數據共分為10個不同土地覆蓋類型，包括：1 耕地，2 林地，3 草地，4 灌木叢，5 濕地，6 水體，7 苔原，8 建設用地，9 裸地，10 冰川與永久積雪.在ARCGIS10.2中將FROM-GLC 2017數據與NDVI影像疊加，將不同植被類型進行分類，裁剪出賀蘭山地區林地和草地的NDVI影像數據.

1.5 NDVI數據

NDVI數據下載自中國科學院資源科學數據中心(http://www.resdc.cn)，該數據是基于SPOT-VEGETATION NDVI衛星遙感數據合成產品，空間分辨率為1 km，時段為1998～2018年.該數據通過最大值合成法(MVC)得到每月1期的月值數據，進一步消除了云和衛星偏角等噪音影響[30].該數據有效反映了在空間和時間尺度上的植被覆蓋分布和變化狀況，對植被變化狀況監測、植被資源合理利用和其它生態環境相關領域的研究有十分重要的參考意義[31].將分類提取后林地、草地的NDVI影像數據轉化得到NDVI時間序列[30].

2 結果與分析

2.1 賀蘭山NDVI的變化特征

1998～2018年賀蘭山5～9月林地NDVI的平均值為0.277，高于草地(0.119)以及未分類情況(0.17)，且林地NDVI變化幅度略高于草地NDVI(圖3-A).對比21年NDVI月均值發現(圖3-B)，1～3月NDVI維持在全年較低水平，4月下旬植被緩慢生長NDVI開始上升，5月植被生長旺盛NDVI上升速度較快，5～8月持續增長，Gao等[32]2017年在賀蘭山地區微樹芯監測油松形成層細胞分生速率的結果顯示(圖3-B)，油松在6～7月形成層細胞分生速率最快.8月溫度開始下降(圖3-C)，但8月的降水量全年最高(41.22 mm)，此時水熱條件下植被仍可生長，NDVI在8月達到峰值，但不同植被類型的NDVI峰值差異較大，林地為0.30，草地為0.13.隨著9月降水開始下降，溫度進一步降低，植被冠層光合作用強度降低，生長速率開始減緩趨于停滯，植被逐漸停止生長(圖3-B).

圖3 賀蘭山地區不同植被NDVI年際(A)、年內(B)變化趨勢，形成層細胞分生速率及銀川氣象站溫度和降水變化(C)Figure 3 Variations of NDVI of different vegetation types on annual (A) and intra-annual (B) scales and cambium cell growth rate in Helan Mountains；distribution of monthly mean temperature and monthly total precipitation of Yinchuan meteorological station (C)

2.2 不同植被類型NDVI與樹輪寬度指數的關系

1998～2018年不同植被類型NDVI與樹輪寬度指數的相關性分析結果顯示(圖4)，上一年7月至當年3月，NDVI與樹輪寬度指數呈負相關且不顯著；在生長季，未分類NDVI與樹輪寬度指數僅在6月和7月呈顯著正相關，相關系數分別為0.484和0.528(n=21，P<0.05)；草地NDVI與樹輪年表在6月和7月也呈顯著正相關(r=0.561，r=0.603，n=21，P<0.01)，且草地的相關系數高于未分類.草地NDVI與樹輪寬度指數在6月和7月呈顯著正相關，相關系數分別為0.561和0.603(n= 21，P<0.01)，林地NDVI與樹輪寬度指數在5～7月均達到0.01顯著性水平，7月相關性最高，相關系數為0.685(n=21，P<0.01).林地NDVI與樹輪寬度指數的相關系數在6～7月高于草地和未分類的情況.不同月份組合的相關性分析結果也表明，林地NDVI(6～8月)與樹輪寬度的相關系數最大(r=0.645，n=21，P<0.01)，5～9月的月份組合兩者同樣具有較高的統計相關性(r=0.634，n=21，P<0.01).NDVI與樹輪寬度指數的顯著相關性說明兩者對限制因子的響應結果較一致，樹輪寬度指數可以反映生長季特別是夏季(6～8月)林地NDVI的變化.

圖4 1998～2018年賀蘭山地區不同植被類型NDVI與樹輪年表的相關系數Figure 4 Correlation coefficients between NDVI of different vegetation types and tree-ring chronology in Helan Mountains during 1998～2018

2.3 樹輪寬度與不同植被類型NDVI對氣候因子的響應

為進一步探究樹木徑向生長對氣候因子的響應，本研究分析了樹輪年表、林地NDVI、草地NDVI與溫度、降水以及PDSI的關系(圖5-A).結果表明，在1951～2018年樹木徑向生長與溫度呈負相關，與降水呈正相關，顯著相關的月份較少，兩者在6月相關性最強；與PDSI的相關性高于對溫度及降水的響應，幾乎在整個水文年均呈現顯著正相關，與5月和6月份相關系數最高，分別為0.56和0.566(n= 68，P<0.05).因為PDSI是降水、溫度及土壤水分等多種因素的綜合指標[33]，是水分供需累積效應的近似度量值，可以更好地指示土壤水分變化，表征土壤對于樹木生長的可供水量.

1998～2018年6～8月林地NDVI和溫度、降水的相關性結果顯示(圖5-B)，林地NDVI與降水主要是正相關關系，而與溫度主要為負相關關系，6月NDVI與降水的相關系數為0.452(n=21，P<0.05).林地NDVI同樣表現出對降水響應的滯后效應.草地NDVI對溫度、降水的響應結果(圖5-C)，逐月相關性不顯著，且相關系數低于林地NDVI與氣候因子的相關性.因此，林地NDVI與溫度、降水及PDSI的相關性高于草地，表明林地對氣候變化可能更敏感.

圖5 樹輪年表(A)與氣候因子、林地NDVI(B)和草地NDVI(C)與氣候因子的相關性分析Figure 5 Correlation coefficients between tree-ring chronology (A)，forestland NDVI (B) and grassland NDVI (C) and climate variables

2.4 賀蘭山夏季林地NDVI的重建

樹輪年表與賀蘭山地區夏季林地NDVI的相關系數為0.645(n=21，P<0.01).基于樹輪年表重建賀蘭山地區1920～2018年的林地NDVI，利用最小二乘法建立線性回歸方程：

NDVI=0.662×STD+0.221

(1)

公式(1)中NDVI表示林地NDVI在6～8月的重建值，STD表示賀蘭山樹輪標準化年表.根據擬合方程，本文重建了賀蘭山1920～2018年林地NDVI的變化(圖6-A).為檢驗重建的可靠性，利用逐一剔除法對擬合方程進行穩定性檢驗，誤差縮減值(RE)為0.357(P<0.01，n=99)，說明重建方程是穩定可靠的.

根據1920～2018年賀蘭山林地夏季NDVI重建結果(圖6-A)，林地生長狀況最好的5 a分別是：1950，1956，1963，1970和1979；最差的5 a分別是：1928，1947，1966，1973和1982年.對NDVI進行10 a低通濾波分析結果顯示，1950～1979年NDVI高于平均值(0.284)，森林生長狀況較好；1922～1933年低于平均值，森林生長狀況較差.將重建的林地NDVI結果(圖6-A)與銀川氣象站點6～8月的溫度、降水等器測數據進行對比(圖6-B～C)，發現1951～1981年該地區的溫度低于6～8月的多年平均溫度22.5 ℃，降水較平均值37.6 mm/月略高，而NDVI處于整個序列平均值之上，樹木生長狀況較好.根據重建結果，賀蘭山地區植被生長較差的時段是1922～1933年，通過旱澇指數發現(圖6-E)，造成植被生長較差的主要原因是持續的干旱.在1920～1930年左右，該地區發生過極端干旱事件，主要是降水的減少導致植被生長狀況較差，該階段NDVI處于低值水平.此外，1982年左右賀蘭山地區再次出現極端干旱事件，NDVI出現極端低值年，但NDVI對該極端干旱事件存在一定的滯后性.為進一步探究重建的林地NDVI與氣象因子的關系(表1)，對原始數據進行一階差分去除趨勢噪音的影響，發現重建的林地NDVI和降水、PDSI的相關性提高.重建NDVI序列與降水的相關系數為0.586(n=69，P<0.01)；與PDSI的相關系數為0.508(n=84，P<'0.01),可見重建結果是可靠的.

圖6 1920～2018年賀蘭山夏季林地NDVI變化歷史(A)與溫度(b)、降水(C)、PDSI(D)、旱澇指數(E)對比Figure 6 Comparison between reconstructed summer NDVI (A) and temperature (B)，precipitation (C)，PDSI (D)，drought and flood index (E) in the Helan Mountains

表1 1920～2018年NDVI與溫度、降水、PDSI、旱澇指數的Pearson相關系數

3 討論

通過NDVI變化特征可以了解賀蘭山地區植被的生長動態歷史.賀蘭山地區不同時段NDVI的增長速率相異，但不同植被類型的年內增長趨勢大體一致，與溫度、降水基本維持同步變化.因為樹木開始生長需要特定的條件，比如當積溫和降水到達一定程度后樹木形成層才開始活動[34-35].賀蘭山地區樹木徑向生長的開始時間大致為4月下旬或者5月初，局部受地形的影響可能在5月之前就已經開始生長，初期生長較為緩慢，這也是5月之前NDVI值較低的主要原因.將NDVI中的林地類型提取后發現其NDVI值高于草地和未分類，可能是由于賀蘭山地區建群種為常綠針葉樹種，一年四季并未完全枯黃落葉，而草地多為一年生草本植物，其經歷發芽、生長、枯黃等過程，導致其NDVI值不及林地.尤其是賀蘭山地區林地冠層的變化，即使該地區針葉林的冠層常綠，但隨著生長季前期溫度升高和降水增多，進入生長季樹木冠層綠度逐漸發生變化，這一變化遙感衛星可以較好的捕獲.另外，2017年在賀蘭山利用微樹芯對形成層細胞分生速率監測結果顯示[32]，細胞分生速率最快的時期在6～7月，而NDVI在8月達到最大值，這可能是因為在樹木生長初期，僅有一部分物質和能量用于樹木徑向生長，還有一部分用于冠層中新生枝葉的生長[36]，冠層密度和綠度逐漸增加，NDVI也逐漸上升.

此外，NDVI是植物葉面由紅光和近紅外2個波段合成計算得到的植被指數[7]，主要反映植物葉片的綠度，與光合作用密切相關.樹輪寬度反映樹木徑向生長狀況，主要取決于光合作用與呼吸作用所產生的部分碳水化合物[34]，另一部分碳水化合物用于根、莖、葉生長發育以及呼吸作用等.森林冠層生長與徑向生長之間存在聯系，例如當干旱發生時，樹木徑向生長與冠層生長減退，表現出較好的一致性.賀蘭山地區植被生長最差的時段是1922～1933年，造成植被生長下降的主要原因是降水的減少和持續的干旱，該極端干旱事件已被較多研究結果所證實[37-39].此外，1982年左右賀蘭山地區出現了極端干旱事件，但NDVI對干旱的響應存在滯后性，體現了樹木對極端干旱的“遺留效應”[40].陳峰等[33]在賀蘭山地區發現油松生長主要受水分的限制，我們通過對比分析多氣候序列指標，同樣發現該地區植被生長受干旱條件的控制，降水是影響樹木徑向生長的主要原因.

賀蘭山地區植被類型主要為林地和草地，生長季在5～9月，分類后林地NDVI值高于草地NDVI，且林地NDVI和樹輪年表的相關系數高于草地NDVI與樹輪年表的相關系數.兩者NDVI的差異可能與植被類型、葉片形狀、葉綠素和水分的含量等因素有關[41-42]，葉綠素對藍色(470 nm)和紅色(670 nm)波段最為敏感，葉綠素作為主要的吸光物質，直接影響到植物光合作用的光能利用效率，含量越低，藍、紅波段吸收減弱，可見光波段反射率升高近紅外反射率減弱，反之葉綠素含量越多，藍、紅波段吸收增強，可見光波段反射率降低，近紅外反射率增強[35]，且近乎所有的近紅外輻射由于被散射掉而很少被吸收.楊紅飛等[43]發現不同植被類型反射光譜曲線形態和特征不同，造成林地和草地冠層光譜反射差異較大.賀蘭山地區樹木在夏季冠層茂密，地面生物量大，相對葉綠素含量較高；該地區草地稀疏，地面生物量低，葉綠素含量和水分含量較低.李明明在賀蘭山地區通過提取NDVI中不同植被冠層物候與油松形成層物候進行分析后發現[44]，NDVI中林地、草地的物候參數存在差異，林地的植被結束生長時間受8月最低溫影響顯著，草地則不顯著.這也說明了不同植被類型對溫度的不同敏感度，也說明不同植被NDVI所包含的植被信息可能存在差異.因此不同植被的生理過程存在差異，例如林地和草地的根系深淺不同，且對環境因子的變化敏感度也存在差異[45].將林地NDVI單獨提取后，樹輪年表與林地的NDVI相關性有所上升，揭示了考慮不同植被類型對重建NDVI的必要性.

4 結論

根據賀蘭山1998～2018年間林地NDVI數據與樹輪年表的分析結果，驗證了植被類型對NDVI重建產生影響.對比分析草地、林地以及未分類的NDVI與樹輪年表的相關性結果表明，分類后的林地NDVI能更好的表示樹木生長信號.根據林地NDVI與樹輪年表重建1920～2018年6～8月林地NDVI的變化序列，顯示在過去的98 a樹木生長狀況波動較大，受20世紀20至30年代持續干旱事件的影響，該地區林地生長狀況較差；20世紀40至80年代降水增加樹木生長狀況較好；20世紀80年代至今生長基本保持穩定狀態.這種樹木生長狀況的變化受局地的水分狀況影響，極端干旱事件的“遺留效應”影響樹木徑向生長.本研究證實了不同植被類型會影響NDVI中樹木生長信號的反映，并為基于樹木年輪重建NDVI的研究提供參考依據.