蘭云飛， 李傳華

(1.北京師范大學(xué)地理科學(xué)學(xué)部， 北京 100875； 2.北京師范大學(xué)災(zāi)害風(fēng)險(xiǎn)科學(xué)研究院， 北京 100875； 3.西北師范大學(xué)地理與環(huán)境科學(xué)學(xué)院， 甘肅 蘭州 730070； 4.甘肅省綠洲資源環(huán)境與可持續(xù)發(fā)展重點(diǎn)實(shí)驗(yàn)室， 甘肅 蘭州 730070)

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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)及預(yù)處理

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

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

1.3 研究方法

1.3.1趨勢分析 線性趨勢法是通過建立氣候變量與時(shí)間之間的一元線性回歸，從而得到1條直線來表示該氣候變量與時(shí)間的關(guān)系[33]。基于像元的一元線性回歸分析法，能夠通過一定時(shí)間段內(nèi)單個(gè)像元的數(shù)值變化，計(jì)算每個(gè)柵格植被NPP的時(shí)空變化趨勢，單個(gè)像元多年線性回歸方程的趨勢線斜率即為年際變化速率[34]。本研究利用該方法逐像元計(jì)算2000—2015年祁連山植被NPP的變化速率，具體公式如下：

(1)

式(1)中，i為年時(shí)序，i=1,2,3,……,16，NPPi為第i年的NPP值；n為年數(shù)，本研究n=16；θslope表示NPP的變化趨勢，當(dāng)θslope>0，監(jiān)測時(shí)段內(nèi)NPP是增加的，反之減少。該公式已廣泛用于植被指數(shù)和NPP的時(shí)間序列分析，具有很好的穩(wěn)定性和置信度[35]。為了檢查回歸模型的有效性，顯著性的檢驗(yàn)是很有必要的，結(jié)合NPP變化趨勢，使用Matlab語言進(jìn)行顯著性檢驗(yàn)，其趨勢劃分為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相關(guān)性分析 相關(guān)性分析是用于衡量兩個(gè)變量間線性關(guān)系的強(qiáng)弱程度。本文采用Pearson相關(guān)系數(shù)法研究植被年NPP與氣溫和降水的關(guān)系，由于氣溫和降水柵格影像與NPP的分辨率不同，將逐年NPP影像重采樣統(tǒng)一分辨率為1 km×1 km，相關(guān)分析計(jì)算公式如下：

(2)

2 結(jié)果和分析

2.1 區(qū)域多年平均氣溫、年降水量和年NPP的時(shí)空分布

2.1.1氣溫、降水和NPP的時(shí)間變化 圖2a,2b分別為祁連山近16年植被NPP與年平均氣溫和年降水量的時(shí)間變化，二者都隨時(shí)間表現(xiàn)出明顯的年際波動(dòng)變化，多年平均氣溫為-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對(duì)氣溫和降水的響應(yīng)，在各年植被NPP、氣溫和降水的基礎(chǔ)上逐像元計(jì)算實(shí)現(xiàn)研究時(shí)段內(nèi)祁連山的年平均氣溫、年降水和植被NPP的空間化表達(dá)(圖3a,3b,3c)。

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

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

圖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可知，祁連山大部分地區(qū)植被NPP是增加的，占NPP分布區(qū)的70%，增加的地區(qū)主要集中在祁連山的中東部，湟水谷地、拉脊山以及冷龍嶺北麓NPP的增加值可達(dá)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的總體趨勢

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

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

2.4 NPP與氣溫和降水相關(guān)分析

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

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

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

總體上看，NPP與氣溫的相關(guān)性由中心(正相關(guān)為主)向四周(負(fù)相關(guān)為主)逐漸減弱，NPP與降水的相關(guān)性由中心(負(fù)相關(guān)為主)向四周(正相關(guān)為主)逐漸增強(qiáng)，NPP與氣溫和降水呈現(xiàn)正負(fù)相關(guān)交替互補(bǔ)的分布態(tài)勢。

3 討論

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

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

隨著衛(wèi)星遙感技術(shù)的成熟，模擬預(yù)估植被NPP的模型逐步得以應(yīng)用，有學(xué)者認(rèn)為綜合多種算法的集合預(yù)估方法能夠有效地提高模型模擬精度，其模擬結(jié)果明顯好于單一算法，這將會(huì)是成為今后關(guān)于植被NPP研究的一個(gè)重要方向[44]。

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

4 結(jié)論

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