2000年以來中國潛在植被凈初級生產(chǎn)力的時空分布模擬

畢 凡,潘竟虎

西北師范大學(xué)地理與環(huán)境科學(xué)學(xué)院, 蘭州 730070

陸地生態(tài)系統(tǒng)是一個巨大、龐雜的系統(tǒng),而碳元素循環(huán)在氣候變化和環(huán)境演化中起著舉足輕重的作用,也是陸地生態(tài)系統(tǒng)各圈層相互作用、相互影響、相互制約的關(guān)鍵因素[1]。作為陸地碳循環(huán)的重要組成部分[2—4],植被凈初級生產(chǎn)力(NPP)是綠色植被一定時間內(nèi)進行光合作用所產(chǎn)生的有機物量減去綠色植被自身的維持呼吸和生長呼吸消耗以后的殘余部分[5—6]。NPP是估算陸地生態(tài)系統(tǒng)碳源、碳匯以及描述碳循環(huán)的重要內(nèi)容。潛在自然植被的凈初級生產(chǎn)力(PNPP)則是指在沒有人類干擾的情況下自然植被的凈初級生產(chǎn)力,一般認(rèn)為是人類對于生態(tài)環(huán)境未造成影響時,自然狀態(tài)下可能產(chǎn)生的NPP值[7]。對于PNPP的定量評估,有助于理解人類活動對生態(tài)系統(tǒng)的影響程度,為實現(xiàn)“碳達峰”“碳中和”目標(biāo)提供科學(xué)依據(jù),為應(yīng)對氣候變化和指導(dǎo)生態(tài)系統(tǒng)恢復(fù)工作提供重要參考。

目前,國內(nèi)外已有許多學(xué)者對植被PNPP展開了研究,Klemm等[8]利用全球動態(tài)植被模型MC2研究了21世紀(jì)氣候變化對潛在植被類型的影響,并模擬了NPP的時空分布。任正超等[9]基于CSCS(Comprehensive sequential classification system)理論,采用改進的CASA(Carnegie-Ames-Stanford Approach)模型模擬中國潛在自然植被NPP及其時空分布格局,并給出了其對地形和氣候的響應(yīng)模式。潘竟虎[7]等利用分類回歸樹模型模擬出了中國潛在歸一化植被指數(shù)(NDVI),用改進的CASA模型模擬了中國PNPP的空間分布。趙軍等[10]以綜合順序分類系統(tǒng)為基礎(chǔ),用NPP分類指數(shù)模型模擬了內(nèi)蒙古自治區(qū)各植被類型的PNPP,并分析了NPP與氣候因子間的關(guān)系。唐正宇[1]采用改進的集成生物圈模型(IBIS)模擬了伊犁河流域植被的PNPP,并對人類活動的影響作出分析。總體而言,現(xiàn)今對于PNPP的研究還較為缺乏,在模擬方式和理論基礎(chǔ)上均不夠成熟,不同學(xué)者模擬出的PNPP及其分布格局也存在較大的差異[7]。

目前對PNPP的模擬多借助于實際NPP的估算模型,常見的方法有：(1)以綜合順序分類法模擬潛在植被類型,利用光能利用率模型對植被PNPP進行模擬[9—11]。(2)識別無人為影響點,利用機器學(xué)習(xí)的方法估算PNPP[12]。(3)利用氣候模型或過程模型來模擬PNPP[1,13]。(4)通過改進實際NPP估算模型中的主要參數(shù),達到計算PNPP的目的[7,14],此方法多用于光能利用率模型的改進。CASA模型是計算NPP及PNPP時應(yīng)用最為廣泛的模型,雖然該模型充分考慮到環(huán)境條件的差異和植被的特性,但在參數(shù)的求算過程中仍然存在以下不足：(1)CASA模型只能在估算光合有效輻射的吸收比例(FPAR)中的比值植被指數(shù)最大值(SRmax)時根據(jù)不同植被類型來取值,不能很好地反映植被類型與NPP的關(guān)系；(2)CASA模型中最大光能利用率(εmax)的取值固定為0.389 gC/MJ,實際上不同植被類型的光能利用率并不相同,最大光能利用率的取值自然也應(yīng)該有所差異；(3)在CASA模型中水分脅迫系數(shù)Wε(x,t)的計算用到了土壤水分子模型,求算過程涉及到大量參數(shù),數(shù)據(jù)的獲取較為困難,并且精度難以保證[15]。基于此,本文引入CSCS法的氣候指標(biāo),改進CASA模型中水分脅迫因子和最大光能利用率的計算取值,利用潛在葉面積指數(shù)(PLAI)估算潛在光合有效輻射的吸收比例(PFPAR),最終模擬出2000—2020年中國PNPP的空間分布,以期提出一種PNPP時空模擬的新視角,豐富人類活動對植被NPP影響的研究,為預(yù)測氣候變化和人類活動對生態(tài)環(huán)境的影響提供科學(xué)依據(jù),進而為制定合理的生態(tài)修復(fù)政策提供參考。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)來源及處理

本文所需氣象數(shù)據(jù)來自于中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/)提供的全國696個氣象站點的2000—2020年氣象觀測數(shù)據(jù),采用AMMRR(Analytic method based on multiple regression and residues)[16—17]的插值方法插值為1km分辨率的柵格數(shù)據(jù)。歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)數(shù)據(jù)來自于Level 1 and Atmosphere Archive and Distribution_System(LAADS)網(wǎng)站2000—2020年的MOD13A3數(shù)據(jù),空間分辨率為1 km,時間分辨率均為月。地表反射率數(shù)據(jù)來自于LAADS網(wǎng)站的MOD09A1數(shù)據(jù),其空間分辨率為1 km,時間分辨率為8天。數(shù)字高程模型(Digital Elevation Model,DEM)數(shù)據(jù)取自中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心(http://www.resdc.cn/Default.aspxDEM),空間分辨率為1 km。植被類型數(shù)據(jù)取自中國1∶100萬植被圖(https://www.resdc.cn/Default.aspx),空間分辨率為1 km。中國干濕分區(qū)矢量邊界是由地之圖網(wǎng)(http://map.ps123.net/china/11468.html)提供的中國干濕分區(qū)圖導(dǎo)入ArcGIS,經(jīng)過幾何校正后手動矢量化得到。

1.2 研究方法

1.2.1PNPP的估算方法

本文以CASA模型為基本計算式,在改進了模型中水分脅迫系數(shù)和最大光能利用率等參數(shù)的基礎(chǔ)上,通過毛德華[14]提出的計算潛在光合有效輻射吸收比例(PFPAR)的方法計算出PFPAR,來代替原CASA模型中的光合有效輻射吸收比例(FPAR),從而達到計算PNPP的目的。這種方法的優(yōu)點在于避免了使用氣象模型時過多強調(diào)氣象數(shù)據(jù)與植被間的線性關(guān)系,從而無法精確估算潛在條件下植被NPP的問題。此外,與基于遙感數(shù)據(jù)估算的實際NPP相比,解決了空間上數(shù)據(jù)一致性較低的問題[14]。PNPP的計算式為：

PNPP(x,t)=PFPRA(x,t)×SOL(x,t)×Tε1×Tε2×Wε(x,t)×εmax×0.5

(1)

式中,PFPAR (x,t)為潛在光合有效輻射吸收比例,SOL(x,t) 表示像元x在t月的太陽總輻射量,Tε1和Tε2表示溫度對光能利用率的影響,Wε為水分脅迫影響系數(shù),εmax為理想狀態(tài)下的最大光能利用率,常數(shù)0.5表示植被所能利用的太陽有效輻射占太陽總輻射的比例。式(1)中主要參數(shù)的計算如下：

(1)潛在光合有效輻射吸收比例(PFPAR)計算

本文中的PFPAR是計算PNPP的重要參數(shù),在以往對實際NPP的研究中,FPAR是用NDVI來代表植被生長的狀況,但NDVI由于遙感的現(xiàn)勢性不能反映出“潛在”的特性,因此本文在改進CASA模型的基礎(chǔ)上,采用毛德華使用的PFPAR計算方法[14],計算PFPAR來代替CASA模型中的FPAR,此方法被用來模擬2000—2010年東北地區(qū)沼澤濕地植被PNPP時,表現(xiàn)出較好的效果。由于此方法是通過整合VPM(Vegetation photosynthesis model)模型和TEM(Terrestrial ecosystem model)模型子模塊得到,而VPM模型和TEM模型都可用于區(qū)域和全球尺度的研究,因此本文將其用于中國PNPP的估算。PFPAR的計算式如下：

PFPRA=1-e-k×PLAI

(2)

式中,k取值為 0.5,PLAI為潛在葉面積指數(shù),參數(shù)的確定主要依據(jù)地表水平衡理論,參考區(qū)域大氣模型系統(tǒng)中陸面過程方案以及生物圈-大氣傳輸模式,使用水文平衡理論可估算平衡潛在蒸發(fā)和降水所需的葉面積,由此產(chǎn)生的葉面積規(guī)避了人類活動的影響[18]。PLAI計算式如下：

PLAI=LAImin+fsw×fst×(LAImax-LAImin)

(3)

其中,LAImax和LAImin分別是最大和最小葉面積指數(shù),其值由不同植被類型決定；fst和fsw分別表示土壤溫度對植被的影響和土壤水分對植被生長的影響。

土壤水分對植被光合作用的影響參照VPM模型中利用地表水分指數(shù)(Land surface water index,LSWI)來計算[14],公式如下：

(4)

式中,LSWImax是各像元LSWI在生長季中的最大值。

地表水分指數(shù)LSWI 的計算式如下：

(5)

式中,ρNIR代表近紅外波段地表反射率,ρMIR代表中紅外波段地表反射率。

溫度對植被光合作用的影響系數(shù)fst計算式如下：

(6)

式中,T是大氣溫度(℃)；Tmax,Tmin,Topt分別表示各單元生長季(4月—10月)植被光合速率所對應(yīng)的最高溫度、最低溫度和最適溫度。

(2)水分脅迫系數(shù)Wε(x,t) 計算

在CASA模型中水分脅迫系數(shù)Wε(x,t)的計算用到了土壤水分子模型,求算過程涉及到大量的參數(shù),數(shù)據(jù)的獲取較為困難,并且精度難以保證。張美玲等[15]在模擬2004—2008年中國草地NPP的研究中,改進了水分脅迫系數(shù),在CASA模型應(yīng)用中取得了較好的結(jié)果,因此,本文引入其對于水分脅迫系數(shù)的改進辦法,計算式如下：

(7)

L(K)=K+0.906K0.5+0.22

(8)

式中,EET為估算的實際蒸散，PET為潛在蒸散。∑θ為>0℃月積溫,K為濕潤度指標(biāo),計算式如下：

(9)

式中,P為月降水量。

(3)最大光能利用率εmax計算

最大光能利用率εmax的取值參考朱文泉等[19]在模擬中國典型植被光能利用率研究中的取值,如表1所示。

表1 最大光能利用率取值/(gC/MJ)

1.2.2時空變化分析方法

(1)年際變化率

研究時間段內(nèi)PNPP年際變化率的計算采用最小二乘法實現(xiàn),計算式如下[20—22]：

(10)

式中,Solpe是2000—2020年P(guān)NPPi的年際變化趨勢；i代表年份,n代表總年份(n=21),PNPPi代表第i年的PNPP。

(2)穩(wěn)定性分析

PNPP空間變化的穩(wěn)定性分析通過變異系數(shù)實現(xiàn)[23—24],計算式如下：

(11)

2 結(jié)果與分析

2.1 PNPP的空間分布

2000—2020年,中國年均PNPP的空間分布如圖1所示。由圖1可知,PNPP在空間分布上總體呈現(xiàn)自西北向東南遞增的趨勢。以400 mm等降水量線為界,位于其東側(cè)區(qū)域的PNPP值幾乎都大于600 gC/m2,而在其西側(cè)的值則幾乎都小于600 gC/m2。PNPP的高值區(qū)域(>800 gC/m2)主要集中在降水充沛、氣溫利于植被生長的大興安嶺、小興安嶺、長白山、華南沿海、東南沿海以及藏南和云南邊境地區(qū)；貴州、重慶和四川雖然也處于濕潤區(qū)內(nèi),但由于其特殊的地形和氣候的影響,夏季潮濕,冬季陰冷,陰天及云霧天氣較多,日照時間較短,導(dǎo)致太陽輻射量少,PNPP值不足800 gC/m2。PNPP次高值區(qū)主要分布在400 mm等降水量線附近,包括新疆南部、西藏中北部、青海中部、甘肅中南部、寧夏、山西、河北、遼寧、吉林北部、黑龍江南部和內(nèi)蒙古東北部地區(qū),PNPP值一般在400—800 gC/m2。PNPP低值區(qū)主要分布在西北地區(qū),這些地方多沙漠、戈壁,年降水量不足200 mm,植被稀少、氣候干燥,植被生產(chǎn)力較低,PNPP值多在0—400 gC/m2之間。新疆的伊犁河谷及阿爾泰山地區(qū)、甘肅的祁連山等地存在的特殊地形導(dǎo)致這些地區(qū)降水量相對較高,PNPP值也高出其周圍地區(qū)。按干濕區(qū)劃來看,干旱區(qū)的PNPP值最低,濕潤區(qū)的PNPP值最高,年均PNPP大小排序為：濕潤區(qū)(898.84 gC/m2)>半濕潤區(qū)(808.95 gC/m2)>半干旱區(qū)(588.399 gC/m2)>干旱區(qū)(335.33 gC/m2)。

圖1 2000—2020年中國年均PNPP空間分布Fig.1 Spatial distribution of annual mean PNPP in China during 2000—2020PNPP: 潛在凈初級生產(chǎn)力Potential net primary productivity

圖2 中國PNPP年均值年際變化 Fig.2 Interannual variation of annual mean value of PNPP in China

2.2 PNPP的年際變化

2000—2020年,中國PNPP年均值為663.62 gC/m2,PNPP年均最高值出現(xiàn)在2013年,為706.19 gC/m2,最低值出現(xiàn)在2020年,為628.41 gC/m2。總體來看,20年間PNPP以年均2.9 gC/m2的速度緩慢增加,但年際波動較大(圖2)。

根據(jù)模擬的PNPP空間分布,計算2000—2020年柵格尺度上PNPP的年際變化,結(jié)果如圖3所示。由圖3可知,20年間PNPP增加的區(qū)域與PNPP減小的區(qū)域面積基本相等。PNPP年際變化率在-2—6 gC/m2(基本不變和小幅減少)的面積就占到了全國陸地總面積的68.24%,這說明PNPP的年際變化雖存在一定的波動,但變化幅度總體而言較小。PNPP年際增加的區(qū)域占全國陸地面積的50.81%,PNPP顯著增加的區(qū)域主要分布在兩大片區(qū)：一是東北地區(qū),包括大小興安嶺、長白山、內(nèi)蒙古呼倫貝爾和遼寧中東部地區(qū)；二是西南地區(qū)的川鄂黔交界地區(qū)。PNPP年際變化減少的區(qū)域則占全國陸地面積的49.19%,PNPP顯著減少的區(qū)域集中在青藏高原,散布在西北干旱區(qū)以及海南島的南部。

圖3 2000—2020年中國PNPP的年際變化率Fig.3 Annual change rate of PNPP in China from 2000 to 2020

2.3 PNPP的空間穩(wěn)定性

為了探索PNPP在空間上的穩(wěn)定性,基于像元尺度計算了中國PNPP的變異系數(shù),借鑒前人對變異系數(shù)的分級標(biāo)準(zhǔn)[24],將PNPP波動程度劃分為五個等級(表2),分析全國植被PNPP的波動變化特征。2000—2020年,中國PNPP變異系數(shù)的平均值為0.088,表明中國PNPP的變化整體較小。圖4是PNPP變異系數(shù)的空間分布,由圖4和表2可知,中國大部分地區(qū)PNPP的變異系數(shù)都較小,變異系數(shù)<0.1的低波動地區(qū)面積達660.82萬km2,占全國總面積的68.87%,說明中國大部分地區(qū)的PNPP都處于較為穩(wěn)定的狀態(tài)。PNPP中等波動變化(0.1<變異系數(shù)≤0.15)和相對較高波動變化(0.15<變異系數(shù)≤0.2)的地區(qū)主要在青藏高原、長江中上游和華北地區(qū)中部,面積分別為245.74萬km2和33.93萬km2,占全國總面積的25.61%和3.54%。高波動變化(變異系數(shù)>0.2)的地區(qū)則主要集中在貴州,零散分布在西藏東部和新疆部分地區(qū),所占比例為1.98%。

圖4 中國PNPP的變異系數(shù)空間分布Fig.4 Spatial distribution of variation coefficient of PNPP in China

表2 中國PNPP的變異系數(shù)面積統(tǒng)計

3 討論

3.1 與相關(guān)研究成果比較

目前尚未有普遍認(rèn)可的PNPP估算方法,相關(guān)成果較少。將本文模擬的植被PNPP多年平均值與其他學(xué)者的研究結(jié)果進行對比(表3)：任正超等[9]利用CSCS方法和MODIS數(shù)據(jù)模擬得到中國1982—2012年自然植被的PNPP平均值為586.74 gC/m2。王玉濤[12]通過識別無人為影響點,估算了2001—2017年中國PNPP,計算出中國PNPP年均值為529.16 gC/m2。潘竟虎等[7]采用分類回歸樹模型模擬中國潛在NDVI,用改進的CASA和潛在NDVI計算得出中國PNPP多年均值為468.94 gC/m2。張美玲等[11]基于草原綜合順序分類法和CASA模型對中國2004—2008年的PNPP進行了計算,得到中國PNPP年均值為503.8 gC/m2。本文估算出的中國PNPP多年均值為663.62 gC/m2,高于前人的研究結(jié)果,其原因可能是本文使用朱文泉[19]模擬的中國植被最大光能利用率,其主要植被類型的最大光能利用率取值基本都大于原式中的固定取值(0.389 gC/MJ),而張美玲等人的研究中其PNPP模擬值是生物量實測值與氣象要素實測值之間進行回歸后用插值方法獲得,受插值方法的限制,其PNPP模擬結(jié)果較低[11]。潘竟虎等的研究中由于其數(shù)據(jù)限制,海南等地的氣象站點較少,氣溫空間化的結(jié)果并不理想,整體偏低[7],可能導(dǎo)致模擬出的PNPP值較低。任正超等[9]和王玉濤[12]模擬的結(jié)果雖與本文的結(jié)果較為接近,但由于采用的估算模型和時間區(qū)間不同,導(dǎo)致模擬的結(jié)果也不盡相同。

表3 不同學(xué)者估算中國PNPP的結(jié)果比較

從模擬結(jié)果的空間分布上來看,本文計算出的PNPP值海南省最高,為1132.8 gC/m2,其次是臺灣、云南和福建,PNPP值分別為1076.5 gC/m2、1030.6 gC/m2和961.1 gC/m2,PNPP值最小的省份是新疆,為323.1 gC/m2,而在任正超等[9]的研究中,全國PNPP值最高的省份是海南省(664.3 gC/m2),其次是福建(520.6 gC/m2)和云南(518.9 gC/m2),最低值同樣是新疆(54.6 gC/m2)。

從模擬的PNPP變化趨勢來看,本文發(fā)現(xiàn)中國自然植被PNPP在年際上總體表現(xiàn)出波動中增加的趨勢,這與任正超[9]和張美玲等[11]的研究結(jié)果類似,但與趙東升等[25]的研究結(jié)果相反。此外,本文模擬出的PNPP最高值在2013年,這與王玉濤[12]的模擬結(jié)果一致,這可能是由于2013年全國太陽輻射均值高于其他年份,2020年達到最低。本文發(fā)現(xiàn)青藏高原和四川盆地PNPP的波動較高,這是因為青藏高原地處亞洲腹地,區(qū)內(nèi)高寒生態(tài)系統(tǒng)極其脆弱,對氣候變化的響應(yīng)極為顯著[26]。四川盆地四周相接壤的高原阻隔了盆地內(nèi)外低層大氣的交換,且盆地暴雨頻繁,降水強度、持續(xù)時間和影響范圍年際變化較大[27],致使該地PNPP波動較大。

3.2 研究結(jié)果的不確定性

本文中的不確定性主要來自與數(shù)據(jù)有關(guān)的不確定性和模型相關(guān)的不確定性兩方面。在數(shù)據(jù)方面,由于植被類型數(shù)據(jù)更新時間長且可獲取的時間有限,因此本文研究的時間序列只選擇了2000—2020年,沒有進行更長時間序列的分析。受數(shù)據(jù)獲取的限制,中國氣象數(shù)據(jù)網(wǎng)中站點較多的年太陽輻射數(shù)據(jù)只提供到2016年,本文2017—2020年的太陽輻射數(shù)據(jù)選擇了太陽輻射日值數(shù)據(jù)集,但其在青藏高原的站點較少,一定程度上影響了太陽輻射的空間化,可能導(dǎo)致估算出的2017—2020年P(guān)NPP值整體偏低。在模型方面,CASA模型是依據(jù)北美地區(qū)植被所建立的,全球各地植被類型差異較大,模型中參數(shù)的修改較為困難,本文將CSCS法中的積溫和濕潤度引入CASA模型水分脅迫系數(shù)的計算中,簡化了計算過程,但同時也必然會在一定程度上犧牲CASA模型的過程機理優(yōu)勢。此外,本文最大光能利用率的取值借鑒朱文泉等[19]的研究,依據(jù)不同的植被類型取值,但其在當(dāng)時的研究對灌木等的實測數(shù)據(jù)較少,加之植被分類的精度也可能也會對模擬結(jié)果產(chǎn)生一定的影響。在廣泛獲取高精度數(shù)據(jù)的基礎(chǔ)上,改進最大光能利用率取值,開展更長時間序列的PNPP時空分布研究,深入分析城市化推進、重大生態(tài)工程的實施等人類活動對NPP的影響,將是下一步研究的重要方向。

4 結(jié)論

本文將CSCS模型中的積溫和濕潤度計算方法引入CASA模型,改進了模型中水分脅迫系數(shù)的計算,利用潛在葉面積指數(shù)估算出PFPAR,模擬了中國PNPP的時空分布。主要結(jié)論如下：

(1)基于CASA模型模擬PNPP,避免了使用氣象模型時過多強調(diào)氣象數(shù)據(jù)與植被間的線性關(guān)系,從而無法精確估算潛在條件下植被NPP的問題。與基于遙感數(shù)據(jù)估算的實際初級生產(chǎn)力相比,解決了空間上數(shù)據(jù)一致性較低的問題。

(2)中國PNPP空間分布差異顯著,整體呈現(xiàn)出東南沿海和東北大小興安嶺地區(qū)高,西北地區(qū)低的分布格局,以400 mm等降水量線為界,其東側(cè)區(qū)域的PNPP值大都>600 gC/m2,而其西側(cè)的值則幾乎都<600 gC/m2。

(3)2000—2020年,中國PNPP值總體呈現(xiàn)上升趨勢,年均PNPP最高值出現(xiàn)在2013年,為706.19 gC/m2,最低值出現(xiàn)在2020年,為628.41 gC/m2。PNPP年際變化小,大部分地區(qū)的PNPP空間變化穩(wěn)定,PNPP變異系數(shù)多年平均值為0.088,青藏高原,四川盆地和貴州等地的PNPP波動變化較大。