薛曉玉, 王曉云, 段含明,2, 頡耀文

(1.蘭州大學(xué) 資源環(huán)境學(xué)院, 蘭州 730000; 2.西華師范大學(xué) 國土資源學(xué)院, 四川 南充 637002)

近年來，全球范圍的氣候異常變化、人類活動不斷增強，這些問題正影響著陸地生態(tài)系統(tǒng)的平衡，導(dǎo)致區(qū)域生態(tài)環(huán)境的破壞，直接或間接地影響到植被NPP變化，使其變化驅(qū)動的研究顯得尤為重要[1-2]。NPP不僅可以直接表征陸地生態(tài)系統(tǒng)的質(zhì)量狀況，而且在保持生態(tài)系統(tǒng)平衡和調(diào)節(jié)區(qū)域生態(tài)環(huán)境中扮演著重要的角色[3]。凈初級生產(chǎn)力(Net Primary Productivity，NPP)是指植物在單位時間和面積上由光合作用產(chǎn)生的有機物總量減去自身呼吸作用消耗有機物的剩余量，這部分能量是其自身用于生長和繁殖的能量，也是生態(tài)系統(tǒng)中其他成員生存和繁衍的物質(zhì)基礎(chǔ)[4]。

植被NPP的估算方法很多，根據(jù)模型對輸入?yún)?shù)的選擇，可將現(xiàn)有模型大體分為3類：氣候相關(guān)模型、機理模型和光能利用率模型[5-7]。氣候相關(guān)模型是根據(jù)野外獲取的采樣點實測數(shù)據(jù)來建立NPP與氣候因子的簡單統(tǒng)計回歸模型，但該模型由于缺乏生態(tài)學(xué)理論作為依據(jù)或只能對潛在NPP進(jìn)行模擬，估算的準(zhǔn)確性不高[5]。機理模型是以氣候、土壤和植被等因子作為輸入?yún)?shù)，模擬生態(tài)系統(tǒng)中的光合作用、呼吸作用和土壤微生物的分解過程，從而模擬植被的能量流動與物質(zhì)轉(zhuǎn)換，涉及到的輸入?yún)?shù)較多，且其中部分參數(shù)不易獲取，單純的機理模型發(fā)展具有局限性[7]。光能利用率模型是以光合作用原理為依據(jù)，以NDVI數(shù)據(jù)、太陽輻射數(shù)據(jù)、氣溫和降水?dāng)?shù)據(jù)及植被類型數(shù)據(jù)為輸入因子，與其他模型相比，對地面實測數(shù)據(jù)的依賴較小，且模型的輸入?yún)?shù)相對較少，易于收集，在區(qū)域尺度和大尺度測算NPP的研究以及全球碳循環(huán)研究中被廣泛采用[8-13]。

NPP的產(chǎn)量不僅僅受到植被本身的影響，而且還會受到外部條件的限制[14]。目前，對于影響NPP的自然因素研究主要集中在氣溫、降水和地形等方面的研究[15-17]。隨著社會生產(chǎn)力的不斷發(fā)展，人類活動對生態(tài)環(huán)境產(chǎn)生的影響逐漸表現(xiàn)出全球性，其對植被生長的影響也日益顯著，因此，定量分析氣候變化和人類活動如何影響植被變化成為當(dāng)前的研究熱點[18-20]。北方農(nóng)牧交錯帶作為中國北方典型的農(nóng)、牧業(yè)系統(tǒng)過渡帶，地理位置突出，對外界的抵抗力較弱，且對氣候變化和人類活動的干擾敏感。20世紀(jì)70年代末期以來，國家先后實施了“三北防護(hù)林”、“退耕還林”和“京津沙塵暴源治理”等生態(tài)工程，出臺了草地生態(tài)保護(hù)補償獎勵政策，努力促使北方農(nóng)牧交錯帶的植被恢復(fù)和好轉(zhuǎn)[21]。深入研究這一區(qū)域NPP時空變化及其驅(qū)動因子，不僅可以揭示氣候變化、人類活動與植被變化的相互作用機制，還可為改善區(qū)域的生態(tài)環(huán)境提供參考。

1 研究區(qū)概況與數(shù)據(jù)

1.1 研究區(qū)概況

本研究參考趙松喬[22]、趙哈林等[23]的研究成果，將北方農(nóng)牧交錯帶界定為多年平均年降水量300～450 mm、降水年變率15%～30%和干燥度指數(shù)0.2～0.5范圍內(nèi)的區(qū)域[21]。為了保持空間的連貫性，并考慮不同時期的北方農(nóng)牧交錯帶邊界變化的特點，構(gòu)建了寬50 km的緩沖區(qū)，由此形成了本研究所采用的研究區(qū)[21]，其介于34°48′—47°19′N，101°43′—126°00′E，面積為83.5萬km2，海拔高度為100～4 700 m。該區(qū)位于東部季風(fēng)區(qū)與西北大陸性氣候區(qū)的過渡帶，年平均氣溫為0～9℃，其西北地區(qū)干旱少雨，植被稀少；東南地區(qū)降水豐沛，氣候濕潤。由于特殊的地理位置，北方農(nóng)牧交錯帶在地貌和氣候的雙重作用下，生態(tài)環(huán)境的承載力較小。

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

NDVI數(shù)據(jù)來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心(http:∥www.resdc.cn/)。該產(chǎn)品經(jīng)過了輻射校正、大氣校正、幾何校正和多波段合成技術(shù)等預(yù)處理，采用最大值合成法生成1998年以來的中國月度(1—12月)、季度(春、夏、秋、冬)、年度的植被指數(shù)數(shù)據(jù)集，數(shù)據(jù)的空間分辨率都為1 km。氣象數(shù)據(jù)為1998—2016年北方農(nóng)牧交錯帶及其周圍166個站點的日平均氣溫、日降水量和30個站點的逐月輻射數(shù)據(jù)，來自于中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/site/index.html)。對該數(shù)據(jù)進(jìn)行累加統(tǒng)計得到相應(yīng)站點的月平均氣溫和月降水量數(shù)據(jù)，與輻射月值數(shù)據(jù)均采用ANUSPLIN方法[24]進(jìn)行空間插值獲得與NDVI數(shù)據(jù)有相同投影和空間分辨率的柵格數(shù)據(jù)。植被類型數(shù)據(jù)來源于國家青藏高原科學(xué)數(shù)據(jù)中心(http:∥data.tpdc.ac.cn/zh-hans/)。該數(shù)據(jù)是基于IGBP分類系統(tǒng)的植被類型數(shù)據(jù)，生產(chǎn)年份為2000年，空間分辨率為1 km，中國地區(qū)的總體分類精度為71%，高于其他數(shù)據(jù)分類產(chǎn)品，可以作為光能利用率模型的輸入數(shù)據(jù)。本文對以上數(shù)據(jù)進(jìn)行了裁剪、投影變換、鑲嵌等預(yù)處理，作為CASA模型的主要輸入數(shù)據(jù)，最終的分析區(qū)域是將研究區(qū)內(nèi)非植被區(qū)域進(jìn)行掩膜處理的植被覆蓋區(qū)。

2 研究方法

2.1 基于CASA模型的NPP估算及驗證

CASA模型[25]是以NDVI數(shù)據(jù)、氣象數(shù)據(jù)(氣溫、降水、太陽輻射)以及植被類型數(shù)據(jù)為輸入信息來驅(qū)動模型。在模型中，NPP的值由植物吸收的光合有效輻射(APAR)和光利用率(ε)2個因子來表示，其估算公式如下：

NPP(x,t)=APAR(x,t)×ε(x,t)

(1)

式中:t為時間；x為空間位置；APAR(x,t)為像元x在t月份吸收的光合有效輻射[MJ/(m2·月)]；ε(x,t)為像元x在t月份的實際光利用率(g C/MJ)。具體解算過程參見文獻(xiàn)[25]。

在NPP估算過程中，會由于一些不可控因素引入誤差，因此需要對估計結(jié)果進(jìn)行精度評價。本文采取與其他學(xué)者的成果進(jìn)行對比的方法，從而驗證CASA模型在本研究區(qū)的適用性。發(fā)現(xiàn)模擬結(jié)果與其他學(xué)者的總體相似(表1)，與李萌萌[26]基于改進(jìn)的CASA模型估算的陜西省植被NPP結(jié)果相比較，各類植被的NPP值基本一致，其中差別較大的是落葉闊葉林，可能的原因是研究區(qū)內(nèi)闊葉林種類不同所致；其與CEVSA模型[27]和GLO-PEM[28]模型估算的結(jié)果存在一定的差異，但其模擬結(jié)果均表現(xiàn)為林地NPP偏高，草地、農(nóng)田偏低，誤差的存在是必然的，但是其值的大小均在可以接受的合理變動范圍之內(nèi)，說明本研究利用CASA模型估算的NPP是具有較高準(zhǔn)確性的，并且具有一定價值，可作為區(qū)域NPP研究的基本數(shù)據(jù)。

表1 本文植被NPP模擬結(jié)果與其他模型研究對比

2.2 Theil-Sen Median趨勢分析

Theil-Sen Median趨勢分析[29]是計算樣本數(shù)據(jù)的中值，計算公式為：

(2)

式中:SNPP為Theil-Sen Median趨勢；NPPi，NPPj分別為第i年和第j年像元的NPP值。該計算方法計算的是n(n-1)/2個數(shù)據(jù)組合斜率的中位數(shù)。當(dāng)SNPP>0時；NPP呈上升趨勢,反之則呈下降趨勢。

2.3 偏相關(guān)方法

本研究的NPP是基于氣象因子作為輸入數(shù)據(jù)利用CASA模型模擬得到的，其與降水、氣溫具有一定的關(guān)系，因此通過相關(guān)分析方法[30]來分析NPP與降水量或氣溫的相關(guān)性，相關(guān)性系數(shù)公式如下：

(3)

在分析多個因素對NPP的影響時，暫不考慮其他因素的影響而單獨研究某一個因素與NPP的相關(guān)程度即為NPP的偏相關(guān)分析，偏相關(guān)系數(shù)計算公式如下：

(4)

式中:rxyz為變量z固定后變量x與變量y的偏相關(guān)系數(shù);rxy，rxz，ryz分別為變量x與變量y、變量x與變量z、變量y與變量z的相關(guān)系數(shù)。

本文采用T檢驗方法對以上的偏相關(guān)系數(shù)進(jìn)行顯著性檢驗，計算公式如下：

(5)

2.4 基于Thornthwaite模型估算植被潛在生產(chǎn)力

本研究選用Thornthwaite紀(jì)念模型[1]來估算潛在NPP，模型計算公式如下：

NPP=3000×[1-e-0.0009695(E-20)]

(6)

式中:NPP為潛在凈初級生產(chǎn)力[g C/(m2·a)]；E為年實際蒸散量(mm)。

另外，E的具體計算式如下：

(7)

L=2000+25T+0.05T3

(8)

式中:r為年總降水量(mm)；L為年最大蒸散量(mm)；T為年平均氣溫(℃)。

2.5 氣候變化和人類活動對植被NPP變化相對影響的評估方法

本文參考郭繼凱[31]、周妍妍[20]等有關(guān)NPP變化的情景設(shè)定方案，通過每年不同類型NPP的變化量來評估氣候變化和人類活動對植被生產(chǎn)力的相對貢獻(xiàn)。其中，第1種是實際NPP的變化量(NA)；第2種是潛在NPP的變化量(Np)；第3種為人類活動對NPP的影響(NH)，通過計算潛在NPP變化量跟實際NPP變化量的差值來獲得(NH=NP-NA)。

通過計算NA,NP,NH的斜率KA，KP，KH來區(qū)分氣候變化和人類活動對NPP的改善和退化作用。KA>0代表實際NPP趨于增加狀態(tài)，KA<0表示實際NPP逐漸減少。KP>0代表氣候變化有利于植被恢復(fù)，KP<0表明氣候變化導(dǎo)致植被退化。KH>0代表人類活動導(dǎo)致植被退化；KH<0表示人類活動促進(jìn)植被生長。

計算NA,NP和NH的變化趨勢，采用的方法是Theil-Sen Median方法。具體情景設(shè)定方案見表2。

表2 各種可能情景下氣候變化與人類活動在植被凈初級生產(chǎn)力變化中相對作用的評價方法

3 結(jié)果與分析

3.1 北方農(nóng)牧交錯帶植被NPP空間分布及變化趨勢

3.1.1 NPP的空間分布 本文以CASA模型為基礎(chǔ)，制作了1998—2016年北方農(nóng)牧交錯帶的植被NPP數(shù)據(jù)，空間分辨率為1 km。植被NPP多年平均值空間分布見圖1，空間上呈西南、東北高，中間低的分布特征，其值為13.87～1 251.06 g C/(m2·a)，主要集中分布在200～500 g C/(m2·a)，總體平均值為364.93 g C/(m2·a)。植被多年平均NPP的高值區(qū)主要位于西南部的祁連山山麓地帶、中部的燕山山脈、呂梁山、恒山等地、東部的大興安嶺地區(qū)，而低值區(qū)主要分布在中西部植被稀少的干旱地區(qū)，植被NPP多年平均值為13.87～200，200～500，500～1 251.06 g C/(m2·a)時，分別占總面積的11.58%，74.43%，13.99%。

3.1.2 NPP的時空變化特征 1998—2016年北方農(nóng)牧交錯帶植被NPP年均值變化見圖2，近19 a間，植被NPP年均值變化范圍在302.70～459.76 g C/(m2·a)，呈顯著的上升趨勢，2013年的植被NPP年均值為最高，達(dá)到了459.76 g C/(m2·a)，2001年的植被NPP年均值最低，每年的植被NPP的平均增加量約為3.95 g C/(m2·a)。

圖1 1998-2016年北方農(nóng)牧交錯帶植被NPP多年平均值空間分布

圖2 1998-2016年北方農(nóng)牧交錯帶NPP年均值變化

采用Theil-Median趨勢分析方法模擬了每個像元的凈初級生產(chǎn)力變化趨勢，反映了19 a間研究區(qū)植被NPP變化的空間特征(圖3)。研究區(qū)內(nèi)植被NPP增加的區(qū)域明顯大于植被NPP減少的區(qū)域，說明植被生長狀況較好。顯著增加的區(qū)域主要分布在研究區(qū)中、西部的偏南地區(qū)及東北平原的部分零散區(qū)域；顯著減少的區(qū)域較少，主要分散分布在內(nèi)蒙古自治區(qū)東北部。其中，植被NPP增加速率大于3 g C/(m2·a)的區(qū)域面積占總面積的56.95%；植被NPP減小速率大于3 g C/(m2·a)的區(qū)域面積占總面積的5.42%；植被NPP增加速率和減小速率小于3 g C/(m2·a)的區(qū)域面積分別占總面積的13.18%，24.55%。

3.2 氣候因子的空間分布及變化趨勢

3.2.1 氣候因子的時空變化特征 根據(jù)降水和氣溫的線性擬合結(jié)果可以看出(圖4)，表明1998—2016年北方農(nóng)牧交錯帶的降水量呈上升趨勢，速率為3.945 4 mm/a(圖4A)。年平均氣溫呈上升趨勢，速率較小，線性回歸增加速率是0.007 3℃/a(圖4B)。近19 a間，2000年的總降水量最低，為316.29 mm，2012年降水量最高，年平均氣溫最低，分別為509.31 mm和4.74℃，2015年的年平均氣溫最高，為6.94℃。

圖3 1998-2016年北方農(nóng)牧交錯帶NPP的Theil-Sen Median趨勢分析結(jié)果

根據(jù)北方農(nóng)牧交錯帶及其周邊166個氣象站降水資料的插值數(shù)據(jù)進(jìn)行趨勢分析，1998—2016年的年降水量空間分布如圖5A所示，可以看出研究區(qū)內(nèi)的降水量具有明顯的空間差異，總體上東南高，西北低。研究區(qū)的最高年均降水量為715.01 mm，最低年均降水量為194.74 mm。受海陸位置和地形等地理環(huán)境的影響，研究區(qū)內(nèi)年均降水量具有明顯的區(qū)域性差異，自東南向西北逐漸減小，大興安嶺、陰山山脈、黃土高原北緣一線的西北地區(qū)年均降水量多在200～400 mm，東南季風(fēng)區(qū)距離海洋較近，雨水充足，降水量在400～600 mm。由研究區(qū)的年平均氣溫空間分布圖(圖5B)可以看出，該地區(qū)的氣溫也有明顯的空間差異，呈現(xiàn)南高北低，東高西低的趨勢。全年年均氣溫小于5℃的區(qū)域主要位于陰山山脈以北和大興安嶺以西以及研究區(qū)西南部青海省和甘肅省的祁連山地區(qū)，其中部分地區(qū)小于0℃。年平均氣溫在10℃以上主要分布在無定河上游的小部分區(qū)域，其余地區(qū)年平均氣溫為5～10℃。

圖4 1998-2016年北方農(nóng)牧交錯帶降水氣溫變化

圖5 1998-2016年北方農(nóng)牧交錯帶年降水量和年平均氣溫空間分布

3.2.2 氣候因子的空間變化特征 近19 a間，研究區(qū)大部分區(qū)域年降水量呈增加趨勢(圖6A)，呈增加、減少趨勢的面積分別占總面積的96.72%，3.28%，增加較為明顯的區(qū)域分布于陰山山脈以南的陜西省東北部和山西省西北部，西遼河平原吉林省鐵嶺市北部、松嫩平原等地。年降水量呈下降趨勢的區(qū)域主要位于陰山山脈以東的內(nèi)蒙古自治區(qū)中東部和六盤山以西的甘肅省東南部等地。在全球變暖的大趨勢下，1998—2016年研究區(qū)內(nèi)70.60%的區(qū)域年平均氣溫呈上升趨勢(圖6B)。增溫明顯的區(qū)域主要分布在大興安嶺的東北部、燕山北部、呂梁山北部、黃土高原西南部、祁連山東北部山麓地帶。年平均氣溫呈下降趨勢的區(qū)域占研究區(qū)總面積的29.40%，主要位于松嫩平原、渾善達(dá)克沙地東南部、鄂爾多斯高原等地。

3.3 NPP變化與影響因子的關(guān)系

3.3.1 NPP變化與氣候變化的關(guān)系 從植被NPP與降水的相關(guān)分布及相關(guān)顯著性分布(圖7A、圖8A)可以看出，研究區(qū)植被NPP與降水量之間主要呈正相關(guān)，植被NPP與年降水量的總體平均偏相關(guān)系數(shù)為0.369 2，研究區(qū)內(nèi)97.70%的像元NPP值和降水量正相關(guān)。NPP和降水量的正相關(guān)的區(qū)域主要分布在黃土高原、鄂爾多斯高原、河套平原、大青山、渾善達(dá)克沙地、大興安嶺、松嫩平原、西遼河平原等大部分地區(qū)，降水的增加可促進(jìn)植被生長狀況的改善，但這些區(qū)域的正相關(guān)關(guān)系并不顯著；而NPP和降水呈負(fù)相關(guān)關(guān)系極顯著的區(qū)域主要位于研究區(qū)西南部祁連山東北部山麓、呂梁山和太行山北部、燕山山脈周邊等地，這些區(qū)域海拔較高或者緯度較高，氣溫較低，說明降水的增多并不能改善該地的植被生長。圖7B和圖8B為研究區(qū)植被NPP與相應(yīng)時期平均氣溫的相關(guān)關(guān)系及相關(guān)顯著性分布，NPP與平均氣溫的整體平均偏相關(guān)系數(shù)為0.004 8，呈負(fù)相關(guān)的像元占比為50.49%，且負(fù)相關(guān)關(guān)系是極顯著的；NPP與平均氣溫呈正相關(guān)關(guān)系不顯著的區(qū)域較多，集中分布在青海省東北部、甘肅省東南部及松嫩平原南部等地區(qū)，研究區(qū)內(nèi)中部和東北部的大部分區(qū)域海拔較高、緯度較高，冬季長，熱量累積較少，從而表明熱量不足會抑制植被生長。

圖6 1998-2016年北方農(nóng)牧交錯帶年降水量和年平均氣溫變化趨勢空間分布

圖7 1998-2016年北方農(nóng)牧交錯帶植被NPP與氣象因子的空間相關(guān)性分布

圖8 1998-2016年北方農(nóng)牧交錯帶植被NPP與氣象因子的相關(guān)顯著性分布

3.3.2 氣候變化和人類活動對NPP的相對影響 從氣候變化和人類活動影響植被NPP變化的空間分布可以看出，氣候變化和人類活動對植被影響的積極作用明顯大于其對植被的消極貢獻(xiàn)(表3，圖9)。氣候變化主導(dǎo)植被恢復(fù)的面積占總面積的47.64%，主要分布在研究區(qū)東部的西遼河平原和松嫩平原，主要土地利用類型為農(nóng)田，氣候變化主導(dǎo)植被恢復(fù)的情況在研究區(qū)中西部的農(nóng)牧交錯帶呈零星分布。人類活動使植被NPP增加的面積占總面積的2.16%，集中分布在甘肅省蘭州市周邊、內(nèi)蒙古自治區(qū)烏蘭察布市和呼和浩特市邊緣交界處。二者共同作用促進(jìn)植被改善的面積占總面積的30.56%，主要分布研究區(qū)中西部的農(nóng)牧交錯地區(qū)。

表3 各影響因素主導(dǎo)植被NPP變化的比例分布

圖9 氣候變化和人類活動主導(dǎo)植被NPP變化的空間分布

氣候變化導(dǎo)致植被NPP退化的面積占總面積的0.56%，主要分散分布在甘肅省白銀市、內(nèi)蒙古自治區(qū)烏蘭察布市等地。人類活動主導(dǎo)植被NPP減少的面積占總面積的17.67%，集中分布在內(nèi)蒙古自治區(qū)中、東部等地，主要土地利用類型為草地和林地。二者共同作用導(dǎo)致植被退化的面積占總面積的1.41%，主要離散分布在甘肅省白銀市西部、內(nèi)蒙古自治區(qū)呼和浩特市東部、烏蘭察布市西部等地。由上述分析可知，在研究區(qū)內(nèi)的甘肅省蘭州市和白銀市交界處、內(nèi)蒙古自治區(qū)呼和浩特市和烏蘭察布市交界處，引起植被變化的影響因素交叉分布。

4 討論與結(jié)論

4.1 討 論

(1) 在分析NPP與影響因素的關(guān)系時，僅考慮了降水、氣溫和人類活動的影響；由于研究區(qū)范圍較大，在分析自然因素影響時，沒有考慮到地形地貌、土壤質(zhì)地等方面，對人類活動主控因素的確定沒有具體化。需要進(jìn)一步對研究區(qū)的NPP變化成因進(jìn)行分段分析，并結(jié)合多方面的影響因子來說明變化成因。

(2) 限于研究區(qū)的范圍較廣、研究尺度較大等原因，沒有進(jìn)行實地植被NPP的野外調(diào)查，缺少實地驗證，在以后的工作中，應(yīng)加強設(shè)計合理的野外觀測試驗，從而進(jìn)一步提高模型的模擬精度。

4.2 結(jié) 論

(1) 1998—2016年，研究區(qū)的NPP年均值空間上呈現(xiàn)西南、東北高，中間低的特點，研究區(qū)植被的NPP年均值主要集中分布在200～500 g C/(m2·a)。植被NPP增加的區(qū)域明顯大于植被NPP減少的區(qū)域，顯著增加的區(qū)域主要分布在研究區(qū)中、西部的偏南地區(qū)及東北平原的零星地區(qū)；顯著減少的區(qū)域分散分布于研究區(qū)的東北邊緣。

(2) 研究區(qū)內(nèi)氣象要素(降水和氣溫)的分布具有明顯的時空差異，年降水量呈東南高、西北低分布，而年平均氣溫分布呈現(xiàn)南高北低、東高西低，在研究時段內(nèi)，年降水量整體呈增加趨勢，大部分地區(qū)的年平均氣溫也呈上升趨勢。

(3) 北方農(nóng)牧交錯帶的植被NPP與降水量之間存在正相關(guān)關(guān)系，而與年平均氣溫具有明顯的負(fù)相關(guān)關(guān)系。研究時段內(nèi)，氣候變化和人類活動對植被影響的積極貢獻(xiàn)明顯大于對植被變化的消極作用，氣候變化、人類活動以及氣候變化和人類活動共同作用造成植被恢復(fù)的面積分別占恢復(fù)總面積的59.28%，2.68%，38.04%；其三者造成植被退化的面積，分別占退化總面積的2.87%，89.96%，7.17%。