2001—2018年中國總初級生產力時空變化的遙感研究

張心竹,王鶴松,延 昊,艾金龍

1 北京林業大學林學院, 北京 100083 2 北京林業大學生態與自然保護區學院, 北京 100083 3 國家氣象中心, 北京 100081 4 益陽職業技術學院,生物與信息工程系, 益陽 413049

植被總初級生產力GPP(Gross Primary Production)是綠色植物通過光合作用吸收大氣CO2制造的有機物總量,驅動著大氣CO2含量的季節和年際變化,體現了陸地生態系統在自然條件下的生產能力[1],在全球碳循環中扮演著重要角色[2]。陸地生態系統在全球碳循環中作為一個主要的碳匯,能抵消大量的人為碳排放,所以準確量化GPP對生態系統功能評估和碳平衡研究至關重要,這已成為地球系統科學領域的研究熱點[3]。

量化GPP的主要方法有地面觀測和模型模擬兩種。地面觀測法多采用渦度相關技術,通過儀器直接獲得相關源區的碳通量,經處理后得到GPP[4]。由于通量站點數量少,分布稀疏,所以地面觀測數據常用于GPP模型的驗證分析[5]。在區域乃至全球尺度上,生態系統過程模型和光能利用率模型是模擬GPP的有效方法[6]。過程模型考慮植物生理生態過程,理論意義強,并且能夠模擬未來生產力變化。然而過程模型擁有不足,其結構復雜,包含大量理想化假設,所需參數較多,尤其是在大尺度應用上,模型精度受到一定影響[7]。

由于衛星遙感技術可以穩定、持續地獲取植被的大尺度空間動態變化信息[8],多種遙感植被數據被廣泛用于植被生產力估算。基于一層大葉模型假設的光能利用率(LUE)理論,建立了大量遙感光能利用率模型,如MODIS GPP標準產品,VPM模型、EC-LUE模型等[9- 11]。但不同模型對植被光能利用率的參數化方案存在較大差異。如MODIS GPP存在一定程度的低估,尤其是在農田地區,與其算法中沒有區分C3、C4植物的最大光能利用率參數有關[5,12]。為此,Yan[13]提出了區分C3、C4植物光能利用率的TEC模型。由于植物冠層對散射輻射較直射輻射具有更高的光能利用率,He[14]提出了區分陰生葉和陽生葉的雙葉光能利用率模型。Yan[15]進一步提出改進的雙葉光能利用率模型(DTEC),考慮了直接輻射和散射輻射及C3和C4植物差異的影響。

日光誘導葉綠素熒光(solar-induced chlorophyll fluorescence,SIF)是綠色植物在自然光照下吸收太陽能進行光合作用時重新釋放的紅光和近紅外光子,是植被瞬時光合作用活動的理想探針,在植物生長監測中有著巨大的應用潛力[16]。大量研究表明SIF與通量站點GPP顯著正相關[17-18]。近年來,GOSAT、OCO- 2、GOME- 2等衛星均能反演植被SIF信息,為大尺度植被GPP研究帶來了新的思路和方法[19-20]。但SIF衛星數據分辨率較粗,對GPP的大尺度模擬和驗證造成了困難。為解決這一問題,Li和Xiao[21]基于OCO- 2 SIF和氣象數據使用立方體集成法(Cubist),開發了時間分辨率為8 d,空間分辨率為0.05°的全球SIF數據集(GOSIF)。

應用各種大葉遙感光能利用率模型在中國已開展了大量研究,在GPP年際變化方面存在較大差異[6,22-26],而采用雙葉遙感GPP模型和基于日光誘導葉綠素熒光SIF的GPP數據集的中國區域研究還很少。為了更好地準確量化和分析中國區域植被GPP時空動態,本文使用GOSIF-GPP和DTEC-GPP對中國2001至2018年GPP進行模擬,著重研究：(1)中國陸地生態系統總初級生產力的時空格局及年際變化；(2)采用中國通量網的6個通量站點開展GPP模型驗證實驗；(3)比較兩種模型GPP的差異。

1 材料與方法

1.1 GOSIF-GPP產品

Li和Xiao[27]基于GOSIF數據和通量站觀測GPP建立起八種不同植被類型的由SIF反演GPP的統計模型,進而由GOSIF數據反演得到2001—2018年全球GPP數據集。本研究使用的是2001—2018年中國區域的GOSIF GPP月尺度產品(http://globalecology.unh.edu)。

1.2 DTEC-GPP模型

DTEC模型是氣象和遙感數據驅動的雙葉光能利用率模型,區分了陰生葉、陽生葉對GPP的貢獻和直射輻射、散射輻射影響[14],并采用了一個改進的水分脅迫系數[13],模型主要結構如下：

GPP=(εmsu×APARsun+εmsh×APARshd)×Wε×Tε

(1)

式中,εmsu和εmsh分別是陽生葉和陰生葉的最大光能利用率,APARsun和APARshd分別是被陽生葉和陰生葉吸收的光合有效輻射[14],Wε和Tε分別是水分和溫度脅迫因子。對于C3植物：εmsh=3.78×Df1.8,εmsu=1.67 gC/MJ。對于C4植物：εmsh=5.78×Df1.8,εmsu=2.56 gC/MJ。

(2)

(3)

(4)

LAIshd=LAI-LAIsun

(5)

PARdif=PAR×Df

(6)

式中,PARdir和PARdif分別是直射光合有效輻射和散射光合有效輻射,LAIsun和LAIshd分別是陽生葉和陰生葉的葉面積指數,PAR是光合有效輻射,Df是散射比[15],其余參數的含義及計算方法見文獻[15]。

水分脅迫因子Wε的定義來自于TEC模型[13]：

Wε=E/EPT

(7)

式中,EPT是潛在蒸發散[28],E是實際蒸散[29]。

1.3 氣象數據

驅動DTEC-GPP模型的氣象因子包括空氣溫度、水汽壓、空氣相對濕度、降水量、風速和日照時數。本研究使用了來自于國家氣象信息中心(http://data.cma.cn)的中國區域2000個高密度臺站2001—2018年的月氣象數據,采用反距離加權(IDW)插值方法將2000個臺站的氣象數據插值到0.05 ×0.05°的空間分辨率。并根據聯合國糧食及農業組織(FAO)的方法[30]由日照時數等氣象數據計算出月地表總輻射和凈輻射。

1.4 遙感數據

驅動DTEC-GPP模型的LAI遙感數據來自于MOD15A2H LAI/FPAR產品(https://lpdaac.usgs.gov/products/mod15a2hv006/),這是Terra衛星上搭載的中分辨率成像光譜儀(MODIS)數據反演的植被參數產品,時間分辨率是8d,空間分辨率是500 m。本文將通過質量控制后的LAI數據重采樣到月尺度5 km空間分辨率,與月尺度的氣象數據結合驅動DTEC GPP模型。

1.5 渦度通量數據

渦度協方差(EC)方法可以測量陸地生態系統和大氣之間的CO2、水和能量通量,因此EC觀測結果常被用于評估各種GPP模型[4]。本文用來驗證GPP的通量數據來自于中國通量網(http://chinaflux.org)提供的六個通量站的GPP數據,包括森林、草地、農田三種生態系統類型,具體站點信息見表1[31]。

表1 中國通量網六個站點的簡要信息

Table 1 The brief description of six ux sites of Chinaflux

表1 中國通量網六個站點的簡要信息

站點名稱Site name縮寫Abbreviation數據年份Year經度/(E)Latitude緯度/(N)Longitude海拔/mAltitude植被類型Biome type優勢種Dominant species 長白山站CBS2003—2008128°05'45″42°24'9″738溫帶針闊混交林紅松、椴樹千煙州站QYZ2003—2008115°03'29.2″26°44'29.1″100亞熱帶常綠針葉林馬尾松、濕地松海北站HB2003—2008101°19'52″37°39'55″3358高寒金露梅灌叢草甸金露梅、矮蓑草當雄站DX2004—200891°03'58.90″30°29'50.21″4333草原化高寒草甸高山蓑草內蒙古站NMG2004—2008116°24'14.4″43°19'31.8″1200溫帶羊草草原羊草、冰草禹城站YC2003—2008116°34'12.72″36°49'44.4″28農田冬小麥、夏玉米

1.6 研究方法

對2001—2018年的DTEC和GOSIF GPP年數據,進行時空變化分析。采用一元線性回歸法在全國尺度和象元尺度分析18年間GPP值的趨勢傾向率,具體采用線性擬合方程(y=k×x+b)的斜率k表征18年間GPP的變化趨勢和幅度。k>0表示GPP為增長趨勢；k<0表示GPP為減小趨勢。

采用六個通量站觀測的GPP數據對DTEC和GOSIF估計的GPP進行了月尺度驗證,使用決定系數(R2)、均方根誤差(RMSE)和偏差(Bias)進行精度評估。

2 結果

2.1 GPP模擬驗證

圖1顯示DTEC模型和GOSIF模型都能較好模擬通量站點GPP的季節動態。在長白山站點,DTEC和GOSIF模型表現十分優異,有與站點GPP幾乎相同的季節動態。在千煙州站,二者在冬季均有一定程度的低估,這可能是由于千煙州站地面異質性大,空間代表性不佳。在當雄和海北草地站點上,兩種GPP模型在夏季生長季都存在高估,尤其是GOSIF GPP,Xiao和 Li[27]也指出GOSIF GPP在生產力偏低的地區存在一定程度的高估。禹城農田站由于種植冬小麥和夏玉米,一年有兩個生長高峰,DTEC和GOSIF模型均能捕捉到兩個生長高峰但同時也均存在低估。值得注意的是,遙感象元數據與地面通量塔觀測數據存在尺度差異,而且地面的高異質性也會使通量站數據的代表性較差,這都會影響驗證效果。

圖1 月尺度上模擬的GOSIF GPP、DTEC GPP和通量觀測GPP的季節動態Fig.1 Seasonal variations of monthly estimated GOSIF GPP, DTEC GPP and flux tower-observed GPP at six flux tower sitesGPP: 總初級生產車Gross primary productivity

對兩種GPP模擬結果的驗證統計見表2,GPP模擬精度隨著通量站點和生態系統類型變化。在森林站點的驗證結果最好,R2為0.88—0.95。在草地站點中,海北站、當雄站的精度較高,R2為0.83—0.91,而內蒙古站R2較低,小于0.6。在農田站點,GOSIF的R2明顯高于DTEC,但兩者的RMSE和 Bias均相對較大,存在明顯低估。GOSIF GPP在3個草地站均存在高估,在森林站和農田存在低估,而DTEC GPP僅在當雄站高估,在其余5個站點均為低估。此外,在各個站點GOSIF GPP精度略高于DTEC GPP。

所有站點的總體驗證結果(圖2)表明,GOSIF GPP的精度(R2=0.78; RMSE=54.67 gC m-2月-1; Bias=-5.15 gC m-2月-1)略好于DTEC模型(R2=0.74; RMSE=66.07 gC m-2月-1;bias=-30.35 gC m-2月-1),DTEC低估程度更大。王克清[32]在與本文相同的六個通量站點處對兩種GPP模型進行驗證,R2分別在0.72—0.89與0.67—0.94之間。Jia等人[33]在8個草地站點對8個模型進行驗證,R2在0.64—0.89之間。Li[24]在中國32個站點對EC-LUE模型進行驗證,C3站點和C4站點的R2分別為0.79和0.62。與上述研究相比,本文的兩種模型精度與之相近。

2.2 中國主要生態系統類型的GPP年際變化特征

2001—2018年主要生態系統類型GPP年值變化曲線見圖3。GOSIF GPP和DTEC GPP模擬的不同生態系統年均值大小順序相同,常綠闊葉林的GPP年均值最高,分別為2140和2014 gC m-2a-1,其次是常綠針葉林和混交林,分別是1713、1717 gC m-2a-1和1464、1452 gC m-2a-1,然后是落葉闊葉林、農田、灌叢、落葉針葉林、草地。草地生態系統平均生產力最低,兩種模型模擬結果均不到300 gC m-2a-1。此外,GOSIF GPP除在常綠針葉林生態系統稍低于DTEC GPP,在其余生態系統類型上均高于DTEC GPP。

從變化趨勢上看,兩個模型模擬的主要生態系統GPP年值從2001至2018年都呈顯著上升趨勢(P<0.01),僅GOSIF GPP在落葉針葉林生態系統中的上升趨勢不顯著。兩種模型的波動幅度有差異,GOSIF GPP波動較小,呈平穩上升,而DTEC GPP波動幅度較大。

表2 模擬的月尺度GOSIF GPP、DTEC GPP在6個通量站點的驗證統計

圖2 GOSIF GPP, DTEC GPP與通量站點GPP月數據的比較Fig.2 Scatter plot of GOSIF GPP and DTEC GPP vs. observed GPP at six flux tower sites on a monthly scale

圖3 2001—2018年中國不同植被類型 GPP 年際變化Fig.3 Interannual variations of GPP of different vegetation types from 2001 to 2018

2.3 中國陸地生態系統GPP時空變化特征

2.3.1中國陸地生態系統GPP多年空間分布格局

從2001—2018年GPP多年均值空間分布(圖4)來看,GOSIF GPP和DTEC GPP的空間分布格局基本相同,都展示了巨大的空間差異性,GPP高值主要分布在東南地區,低值主要分布在西北地區,總體上呈現由東南向西北、從沿海向內陸遞減的分布趨勢,這與已有研究結果[6,23-24,26]基本一致。

其中華南地區、華東地區、云南中南部及青藏高原南部GPP多年均值高于2000 gC m-2a-1,西南地區、長江中下游地區GPP多年均值為1500 gC m-2a-1左右。南方地區生產力高主要是由于氣候溫暖濕潤、光照充足和植被覆蓋度高。此外,東北森林地區也具有較高的GPP年均值,約為1000 gC m-2a-1。而西北的沙漠和荒漠地區、青藏高原中北部以及內蒙古中西部的草原區,年均GPP不足500 gC m-2a-1,這些地區由于水分脅迫和低溫制約,植物生長季偏短,導致了生產力低下。內蒙東部及青藏高原東南部為過渡地帶,GPP多在500 gC m-2a-1左右,與我國半干旱半濕潤分界線基本一致,說明降水是影響我國植被生產力分布的一個重要因素。華北、黃淮大部分地區屬于農業區,GPP多年均值在800—1500 gC m-2a-1之間。此外,新疆天山山脈附近GPP能達到500 gC m-2a-1左右,這是由于來自大西洋的盛行西風和來自北冰洋的氣流,碰到山坡抬升,產生較強降水,促使天山西部及西北部具有較高的植被覆蓋度[34],同時天山的高山融雪也會促使山前地帶綠洲的形成。

圖4 2001—2018年中國年均GPP空間分布格局Fig.4 Spatial pattern of annual GPP in China during 2001 to 2018

圖5 2001—2018年中國植被 GPP 年際變化Fig.5 Interannual variations of GPP in China from 2001 to 2018

2.3.22001—2018年中國陸地生態系統GPP總量年際變化

2001—2018年GOSIF GPP和DTEC GPP的年際變化如圖5,總體均呈顯著增加趨勢,但年際間有波動,其中DTEC GPP變化波動更劇烈。使用GOSIF數據得到的GPP年總量變化范圍為6.35—8.08 PgC/a,均值為7.23 PgC/a,年增幅約為0.094 PgC/a(P<0.01)。 DTEC模型得到的GPP年總量變化范圍為6.20—7.79 PgC/a,均值為6.93 PgC/a,年增幅約為0.073 PgC/a (P<0.01)。無論是GPP總量還是增長幅度,GOSIF GPP均大于DTEC GPP。DTEC GPP的年值曲線在2003、2005、2009、2014年分別出現過較大的波動,這與劉剛等[35]模擬的2001—2014年中國植被凈第一性生產力的結果類似,在2003、2005、2009、2014年都出現GPP低值,主要是這些年份遭受到不同程度的氣象災害。本研究期間,2001—2018年,中國GPP顯著上升,這是由適宜的氣候變化、CO2的施肥效應、氮沉降、人類活動(如植樹造林、農業灌溉)等多方面因素共同作用的結果[1,36-38]。

2.3.3中國陸地生態系統GPP多年時空變化特征

2001—2018年我國GPP時空變化見圖6,近18年間GOSIF GPP和DTEC GPP年變化趨勢基本相同,上升趨勢和下降趨勢并存,但呈上升趨勢的面積更大。變化趨勢率k值總體上呈東南地區大于西北地區、沿海地區大于內陸地區的分布。年GPP呈上升趨勢最高的地區包括華南、華中、華東地區、西南地區東南部以及海南島,趨勢率k能夠達到20—30 gC m-2a-1,其次為陜西、山西、河北、東北大部地區,上升幅度約15 gC m-2a-1左右。此外,GOSIF GPP數據結果顯示新疆天山山脈附近區域GPP也有明顯的上升趨勢,個別地區上升幅度甚至能達到20 gC m-2a-1。上述地區的GPP增長趨勢與已有的研究基本一致[22,35,39]。西北地區、西藏、內蒙古中東部變化幅度較小,基本在-5—5 gC m-2a-1之間。GPP下降趨勢比較嚴重的地區則分布在青藏高原南部、云南中北部和上海、天津等城市區域,下降幅度約為-10 gC m-2a-1。此外,DTEC GPP顯示河南南部地區有較明顯的下降趨勢,大約在-20 — -5 gC m-2a-1之間。青藏高原南部、云南中北部和河南南部地區主要是由于降雨偏少造成干旱事件引發GPP下降[25,40]。而上海、天津等城市GPP的下降趨勢是城市建設用地擴張造成植被生產力下降。總體上,我國長江以南大部分地區GOSIF GPP的年變化趨勢率k高于DTEC GPP。

圖6 2001—2018年中國年GPP變化趨勢Fig.6 Trend of annual GPP in China during 2001 to 2018

3 討論

3.1 中國GPP年際變化

在已有的中國GPP研究中,光能利用率模型、過程模型和以機器學習法為主的統計模型均有廣泛應用。由于時間段選取、模型選擇、參數確定、數據集來源不同以及各種誤差,GPP模擬結果也有較大差異(表3)。其中多年平均GPP最大值為使用支持向量回歸方法[41]得到的2000—2015年間的7.81 PgC/a,最小值為使用EC-LUE模型[42]得到的的2000—2009年間的5.38 PgC/a。本研究使用GOSIF和DTEC雙葉模型得到的GPP年總量分別為7.23PgC/a和 6.93 PgC/a,均處在已有結果的范圍之內。采用多種類型的GPP模型從多角度開展大尺度GPP研究,可以減小陸地生態系統碳循環研究的不確定性。

從年際變化趨勢上看,本研究發現中國植被生產力從2001至2018年呈上升趨勢,與已有基于衛星觀測和生態系統模型的研究結論基本一致[46]。但本研究中GOSIF和DTEC GPP上升趨勢分別為0.094 PgC/a和0.073 PgC/a,GOSIF GPP上升幅度更大。而已有研究得到的GPP上升趨勢為0.02 PgC/a至0.057 PgC/a(表3),低估了中國GPP增長的趨勢。有研究表明,目前對GPP年際上升趨勢的估算存在明顯低估的現象[6,43],與本研究結論一致。本文選取的GOSIF算法和DTEC模型是目前最新的研究方法,在機理上都具有深刻意義,以中國地區為研究對象也是首次。

表3 不同模型估計的中國年均GPP

3.2 GOSIF算法與DTEC模型對比

GOSIF算法與DTEC模型是兩類不同的模型,在不同方面擁有優勢。從機理上看,日光誘導葉綠素熒光SIF包含大量光合信息,能準確地反映植物光合作用,與植物生理過程高度耦合[16],在葉片、冠層、植株、生態系統等不同尺度上均與GPP顯著相關[17-18]。GOSIF是目前時間和空間分辨率最高的SIF全球產品,GOSIF GPP根據植物葉綠素進行光合作用時釋放的SIF輻射直接估計GPP,已經包含各種環境脅迫和人為因素的影響[21,27,47],能夠更加客觀的反應植物的真實生長狀態,捕捉更多空間分布上的細節。但是GOSIF GPP在生產力較低地區出現高估,精度受統計模型和訓練數據的影響。DTEC GPP根據溫度和水分脅迫下植物對光合有效輻射的利用率進行GPP估計,考慮了散射輻射的作用,采用改進后的水分脅迫參數,區分了C3、C4植物參數化,在冠層復雜的森林和干旱少雨的地區有良好表現[15],較傳統大葉LUE模型有較大改進。但由于模型輸入參數較多,并且沒有考慮人為灌溉和西部地區高山融雪的影響,可能會引入誤差。例如DTEC GPP主要是通過降雨量驅動土壤水分平衡方程估計水分脅迫,這導致DTEC GPP在人為經營程度高的地區,如河南和河北的農田、新疆天山附近的綠洲和一些集約化經營管理的林地會低估GPP。而GOSIF GPP包含了自然降雨和人工灌溉的影響,能較客觀的反映出自然和人為影響下的植被生產力。因此,使用這兩類GPP模型有助于理解自然和人為因素對植被生產力的不同影響。

4 結論

DTEC和GOSIF模型都揭示出2001至2018年中國陸地植被GPP呈顯著增加趨勢,平均每年分別增加0.073Pg C和0.094Pg C,在植被稠密的東南地區GPP增幅高于植被稀疏的西北地區。中國DTEC和GOSIF GPP多年平均值分別為6.93 Pg C和7.23 Pg C,GPP空間格局呈現出明顯異質性,從南到北、從東到西GPP年均值逐漸減小。在6個通量站的驗證表明DTEC和GOSIF模型表現良好,具有可用性。

比較而言,GOSIF GPP算法能客觀地反映植被生產力狀況,包括人為影響嚴重的灌溉農田和綠洲地區；而DTEC模型則更適合自然條件下植被生產力的模擬,研究光、溫、水等氣象因子對植被生產力的影響。未來仍需要改進GPP估計模型和方法,采用多種GPP模型研究方法,探索氣候和人為因素對植被生產力的影響以及植被對氣候變化的反饋。