近16年祁連山植被NPP時空格局及其對氣候變化的響應

蘭云飛， 李傳華

(1.北京師范大學地理科學學部， 北京 100875； 2.北京師范大學災害風險科學研究院， 北京 100875； 3.西北師范大學地理與環境科學學院， 甘肅 蘭州 730070； 4.甘肅省綠洲資源環境與可持續發展重點實驗室， 甘肅 蘭州 730070)

植被凈初級生產力(Net primary productivity，NPP)是指綠色植物在單位時間和單位面積上，由光合作用所產生的有機干物質總量中扣除自養呼吸后的剩余部分[1-2]。自工業革命以來，人們對于環境的過度干擾與破壞已經給整個生態系統造成了巨大的壓力。大量化石燃料的燃燒使得大氣中CO2等其他溫室氣體濃度增加，并觸發了一系列環境問題，如全球溫度升高、干旱問題加劇、植被地域分布改變等[3]。IPCC第五次氣候變化評估報告指出，以溫度升高為特點的氣候變化是過去半個世紀以來全球變化的主要內容[4]。其中，北半球中緯度地區增溫最快，預估未來氣候將持續變暖[5]。植被在陸地生態系統中占據重要位置，對氣候變化響應也最為敏感[6]，NPP作為植被生態系統中物質與能量轉換和傳遞的基礎，在氣候變化研究中可充當“指示劑”的作用[7]。因此，有必要為植被NPP的時空格局及其對氣候變化的響應進行長期序列的監測與分析。

隨著國內外碳循環研究的不斷深入，模擬估算NPP的諸多模型被提出，如氣候生產力模型、生理生態過程模型、光能利用率模型等[8-9]，遙感技術的發展為NPP的模擬估算提供了強有力的手段，并取得了豐富的成果[10-11]，其中，基于遙感數據的CASA模型在當前NPP的研究中應用最為廣泛，已在我國西部山區乃至全國碳循環的研究中得到了廣泛認可[12-15]。美國國家航空航天局(NASA)于2016年發布了EOL/MODIS(對地觀測系統/中分辨率成像光譜儀)遙感數據(MOD17A3)，該數據提供了現成的NPP遙感產品，大大簡化了NPP與碳循環的實驗工作，越來越多的被學者[16-18]接受，在此之后，MOD17A3HGF 6.0版的NPP年度數據產品相繼被開發出來，相較于5.5版的MOD17A3數據集，該數據集采用了新的生物屬性調查表(BPLUT) 和新版的全球模型與融合室(GMAO)的日值氣象數據對NPP數值進行模擬，提高了NPP的估算精度[19]。但對于我國西部山區NPP與碳循環的研究，利用MOD17A3產品的還較少，MOD17A3HGF的更是鮮見。

在研究植被與氣候變化的關系中，逐像元相關分析是一種最為主要的研究手段，該方法往往需要對地面氣象站點數據進行插值處理，此時氣象數據插值的精度則是需要認真思考的問題。已有研究[20]表明，提高氣象要素插值的精度有利于NPP的估算研究。目前在西部地區開展的研究以運用常規的插值方法(反距離權重、普通克呂格法等)較多[21-23]，由澳大利亞科學家Hutchinson編制的ANUSPLIN軟件目前在氣象數據的插值處理中已得到較為廣泛的應用[24-25]，但在我國西部地區的應用還較少，也有學者如張仁平等[26]利用ANUSPLIN軟件對新疆67個氣象臺站的氣溫、降水和日照時數插值處理，分析了新疆草地NPP對氣候變化的響應，但鑒于新疆國土面積廣闊，研究所用的站點還是略少。ANUSPLIN軟件已被證明比常規插值法的精度更高[27-28]，在復雜地表且氣象要素空間差異性較大的區域也非常適用[29]，但我國氣象站點東密西疏，要在西部山區應用ANUSPLIN軟件提高氣象數據插值精度，同時兼顧站點稀少帶來的誤差則顯得力不從心。祁連山是我國西北地區重要的生態保障區，生態環境戰略地位突出，而祁連山及周邊地區站點稀少且布局不均，通過離散的站點插值模擬還是會產生一定的誤差。

鑒于此，本研究采用以往應用匱乏的MOD17A3HGF NPP數據產品，氣溫與降水來自于祁連山及周邊多站點(圖1)的ANUSPLIN插值的柵格數據集，通過一元線性回歸分析、逐像元相關分析等方法定量研究2000—2015年祁連山植被NPP的時空變化及其影響因素，旨為祁連山森林碳匯評估、地表植被估產以及生態環境修復提供參考價值。

圖1 祁連山及周邊地面氣象站點分布Fig.1 Distribution of ground meteorological stations in Qilian Mountains and surrounding areas

1 材料與方法

1.1 研究區概況

祁連山位于西北內陸干旱半干旱地區，是黃土、蒙新、青藏三大高原的分界線[30]。東起烏鞘嶺，西至當金山口與阿爾金山相連，由一系列西北—東南走向的山脈與谷地組成，地勢自東向西傾斜，大部分海拔在3 500～5 000 m左右。祁連山兼具大陸性與高原性的氣候特征，自然條件復雜，水熱格局差異明顯。年均溫為0.6℃，降水量表現出從東向西遞減的變化趨勢，山區降水量為400～700 mm，而年蒸發量則表現為相反的趨勢[31]。植被分布受大氣環流與地形的共同影響海拔由低到高呈現荒漠草原、山地草原、山地森林草原、高山灌叢草甸、高寒草甸和高寒稀疏草甸的景觀格局。

1.2 數據及預處理

本研究采用的植被NPP數據來源于美國NASA EOS/MODIS的2000—2015年的MOD17A3HGF 6.0版本的全球NPP數據集(http://www.ntsg.umt.edu/)，空間分辨率為500 m。利用MRT工具及Python語言代碼，將下載的三景影像批量拼接，統一坐標后，裁剪生成2000—2015年祁連山地區NPP數據，對各年份NPP數據乘以轉換因子并剔除無效值(置為Nodata)，無效值在本研究中不參與數值的分析與計算。

氣溫和降水來自中國科學院資源環境科學與數據中心(http://www.resdc.cn/)的逐年年降水量和年平均氣溫的空間插值數據集，該數據集是基于全國2 400多個氣象站點日觀測數據，采用ANUSPLIN軟件，以DEM為協變量插值生成，空間分辨率為1 km×1 km。值得一提的是，本研究插值的精度主要通過平均絕對誤差(MAE)、平均相對誤差(MRE)和均方根誤差(RMSE)3個參數評估，計算公式見文獻[32]，相較于其它插值方法，ANUSPLIN插值精確性明顯提高。

1.3 研究方法

1.3.1趨勢分析 線性趨勢法是通過建立氣候變量與時間之間的一元線性回歸，從而得到1條直線來表示該氣候變量與時間的關系[33]?；谙裨囊辉€性回歸分析法，能夠通過一定時間段內單個像元的數值變化，計算每個柵格植被NPP的時空變化趨勢，單個像元多年線性回歸方程的趨勢線斜率即為年際變化速率[34]。本研究利用該方法逐像元計算2000—2015年祁連山植被NPP的變化速率，具體公式如下：

(1)

式(1)中，i為年時序，i=1,2,3,……,16，NPPi為第i年的NPP值；n為年數，本研究n=16；θslope表示NPP的變化趨勢，當θslope>0，監測時段內NPP是增加的，反之減少。該公式已廣泛用于植被指數和NPP的時間序列分析，具有很好的穩定性和置信度[35]。為了檢查回歸模型的有效性，顯著性的檢驗是很有必要的，結合NPP變化趨勢，使用Matlab語言進行顯著性檢驗，其趨勢劃分為5類：極顯著增加(Slope>0，P<0.01)，顯著增加(Slope>0，0.01≤P<0.05)，無明顯變化(P≥0.05)，顯著減少(Slope<0，0.01≤P<0.05)，極顯著減少(Slope<0，P<0.01)。

1.3.2相關性分析 相關性分析是用于衡量兩個變量間線性關系的強弱程度。本文采用Pearson相關系數法研究植被年NPP與氣溫和降水的關系，由于氣溫和降水柵格影像與NPP的分辨率不同，將逐年NPP影像重采樣統一分辨率為1 km×1 km，相關分析計算公式如下：

(2)

2 結果和分析

2.1 區域多年平均氣溫、年降水量和年NPP的時空分布

2.1.1氣溫、降水和NPP的時間變化 圖2a,2b分別為祁連山近16年植被NPP與年平均氣溫和年降水量的時間變化，二者都隨時間表現出明顯的年際波動變化，多年平均氣溫為-0.4℃，有弱的下降趨勢，多年平均降水量約405 mm，微弱增加，在氣溫和降水的共同作用下，植被NPP呈上升趨勢，增加幅度為13.2%。近16年NPP的變化介于174.8～224.3 gC·m-2·a-1之間，多年平均值為201.9 gC·m-2·a-1，為更深入的了解植被NPP對氣溫和降水的響應，在各年植被NPP、氣溫和降水的基礎上逐像元計算實現研究時段內祁連山的年平均氣溫、年降水和植被NPP的空間化表達(圖3a,3b,3c)。

圖2 2000—2015年祁連山植被NPP與氣溫和降水的時間變化Fig.2 Temporal changes of vegetation NPP and temperature and precipitation in Qilian Mountains from 2000 to 2015

2.1.2氣溫、降水和NPP的空間格局 祁連山多年平均NPP的空間格局差異明顯，表現出自東向西逐級遞減的分布規律，與祁連山的水(圖3a)熱(圖3b)空間分布的趨勢基本一致。90%以上的區域植被NPP值介于0～400 gC·m-2·a-1，最高值在700 gC·m-2·a-1以上，NPP的高值區面積較小，零星分布于祁連山的東南部。湟水流域河谷地帶氣溫較高、降水充沛，水熱條件較好，多年NPP平均值可達600 gC·m-2·a-1以上；中西部地區，地表景觀以荒漠、裸巖和冰川為主，自然環境惡劣，植被稀少，多年NPP均值接近于0。另外，冷龍嶺、托來山和青海南山等幾條西北—東南走向的高大山脈NPP也呈現與氣候環境相適應的低值分布狀態。

圖3 2000—2015年祁連山氣溫，降水和NPP的年均空間分布Fig.3 Annual average spatial distribution of temperature,precipitation and NPP in Qilian Mountains from 2000 to 2015

2.2 NPP的空間變化

由圖4可知，祁連山大部分地區植被NPP是增加的，占NPP分布區的70%，增加的地區主要集中在祁連山的中東部，湟水谷地、拉脊山以及冷龍嶺北麓NPP的增加值可達120 gC·m-2·a-1以上。NPP的空間變化主要以增減值為0～40 gC·m-2·a-1和-40～0 gC·m-2·a-1為主，分別占總面積的43%和29%，集中分布在青海湖西北部的山脈群，NPP減少值小于-40 gC·m-2·a-1僅有不到1%的比例，零星分布于大通河流域河谷地帶。

圖4 2000—2015年祁連山NPP的逐年累積變化Fig.4 Annual cumulative changes of NPP in Qilian Mountains from 2000 to 2015

2.3 NPP的總體趨勢

一元線性回歸分析逐像元計算植被NPP時空變化的結果(圖5a)表明，近16年來祁連山NPP的總體趨勢以增加為主，中東部如湟水谷地、冷龍嶺、青海南山等祁連山的周邊地區NPP增長較快，最大增量可達34 gC·m-2·a-1；NPP呈減少的區域零星分布于中部地區，最大減少量為-17 gC·m-2·a-1，NPP減少表明祁連山近年來局部地區有植被退化的趨勢。研究時段內NPP變化的顯著性檢驗結果(圖5b)與NPP總體趨勢的分布基本吻合，在NPP的分布區中，近57%的區域NPP的變化趨勢并不顯著，顯著增加和極顯著增加的面積分別占23%和20%，而顯著減少和極顯著減少的區域不到1%，說明研究時段內祁連山的生態健康水平總體趨于改善。

圖5 2000—2015年祁連山年均NPP空間變化趨勢及顯著性檢驗Fig.5 The spatial variation trend and significance test of the annual average NPP in Qilian Mountains from 2000 to 2015

2.4 NPP與氣溫和降水相關分析

逐像元計算了2000—2015年祁連山植被NPP與氣溫的相關系數，結果如圖6a所示。NPP與氣溫的相關系數介于-0.88～0.91，平均值為0.35，整體上呈現正相關為主，正負相關并存的空間格局。烏鞘嶺以東、西北部的大雪山以及青海南山的南坡NPP與氣溫負相關較為明顯，負相關的面積約占NPP分布區的12%。NPP與氣溫為正相關的分布范圍廣泛且集中連片，青海湖附近區域相關系數普遍可達0.60以上。

利用同樣的方法計算了NPP與降水的相關系數，由圖6b可知，NPP與降水的相關系數介于-0.76～0.87，基于像元計算的NPP與降水的空間關系同樣是正相關為主，正負相關并存，但NPP和降水為負相關的區域僅有10%，遠小于與氣溫為負相關的區域，負相關與關系不明顯(-0.10

圖6 NPP與氣溫和降水的相關性分布Fig.6 Correlation distribution of NPP with temperature and precipitation

總體上看，NPP與氣溫的相關性由中心(正相關為主)向四周(負相關為主)逐漸減弱，NPP與降水的相關性由中心(負相關為主)向四周(正相關為主)逐漸增強，NPP與氣溫和降水呈現正負相關交替互補的分布態勢。

3 討論

目前，利用MODIS產品估算模擬整個祁連山地區NPP的研究還很少，但在其它地區已有先驗的成果[36-37]可供參考，劉旻霞等[38]利用作物產量估算的NPP值與其MOD17A3產品進行相關分析，結果發現2種數據值呈顯著正相關(R2=0.77，P<0.01)，均值標準誤差為3.95。李玲等[39]利用祁連縣的凈初級生產力實測數據，得出模擬值和實測值之間的相關性達到顯著水平(R2=0.71，P<0.01)，驗證了MOD17A3數據應用于環青海湖地區估算NPP的可行性。本文使用的MOD17A3HGF產品是MOD17A3的更新版，精度質量與空間分辨率更高，其結果也具有可比性。

從NPP的時間變化上看，本研究與張禹舜[40]、孫力煒等[41]均得出近年來祁連山植被NPP是逐漸增加的一致結論，與祁連山植被覆蓋變化的研究結果[42-43]是相符合的，就NPP的空間分布而言，與上述學者[40-41]得出的結論也是一致的，且NPP的分布也較好的反映了水熱分布的信息。值得一提的是，本研究利用MOD17A3HGF產品估算的多年NPP均值與前述學者的結果相比有一定的偏差，這可能與研究時限、不同模型選取的參數以及估算方法的不同有關。

隨著衛星遙感技術的成熟，模擬預估植被NPP的模型逐步得以應用，有學者認為綜合多種算法的集合預估方法能夠有效地提高模型模擬精度，其模擬結果明顯好于單一算法，這將會是成為今后關于植被NPP研究的一個重要方向[44]。

總體上看，本研究利用了在山區插值精度較為理想的ANUSPLIN軟件插值的氣象數據集，較好的克服了以往我國西部山區因站點較少、插值方法不理想而導致的氣象數據插值精度不高的問題，為評估分析植被NPP與氣候變化的關系提供了更加可靠的依據，但是影響NPP的氣象因子較多，僅選擇氣溫和降水兩個因子并不能很好的表征氣候變化的復雜性，同時植被NPP還受區域地貌格局、人類活動的影響，更多影響植被NPP的因素還有待今后的深入研究。

4 結論

本研究利用MOD17A3HGF產品，基于像元尺度計算和分析了2000—2015年祁連山植被NPP的時空格局及其對氣候變化的響應，結果表明：近16年來祁連山NPP的平均值為201.9 gC·m-2·a-1，NPP均值空間分布自東向西逐漸減少，與水熱分布的趨勢一致。湟水流域河谷地帶氣溫較高、降水充沛，水熱條件較好，多年NPP均值可達600 gC·m-2·a-1以上；中西部地區地表景觀以荒漠、裸巖和冰川為主，自然環境惡劣，多年NPP均值接近于0。受氣溫和降水的共同影響，NPP呈波動增加的趨勢，湟水谷地、冷龍嶺等地區NPP值增加明顯，顯著增加和極顯著增加的區域可達23%和20%，顯著減少和極顯著減少的區域不到1%。NPP與氣溫和降水的關系呈現以正相關為主、正負相關互補并存的分布態勢，與氣溫的相關性由中心(正相關為主)向四周(負相關為主)逐漸減弱，與降水的相關性則由中心(負相關為主)向四周(正相關為主)逐漸增強。