海南島近20 年GPP 變化格局及驅動因素分析

雷濟舟，崔 嵬，朱濟帥，張潤卿，趙俊福，章 杰，張 翔，孫仲益,5

（1. 海南大學 生態與環境學院，海口 570228; 2. 國家林業和草原局發展研究中心，北京 100714;3. 海南長光衛星信息技術有限公司，海口 570311; 4. 海南省生態環境監測中心，海口 571126;5. 海南省農林環境過程與生態調控重點實驗室，海口 570228）

陸地生態系統（以下稱為陸生系統）總初級生產力（Gross Primary Production，GPP）是指大氣中的CO2通過植被光合作用進入陸生系統并轉化為有機碳的總量[1]，在一定程度上決定著碳匯時空分布格局[2]，影響著生態系統物質及能量流動[3]；為實現區域碳平衡，減緩氣候變化起到關鍵作用[4]，具有重要意義[5]。但受到生態系統及其與環境要素間相互作用的影響，GPP 的時空分布具有較大的變異性與異質性特征[6]，因此，對GPP 分布格局的驅動因素定量分析將會有助于提升區域碳循環的理解，也更有助于揭示區域生態系統對人類活動與區域氣候變化的響應與反饋。氣候變化以及人類活動作為影響熱帶地區GPP 分布格局的主要驅動因素，相關研究已取得豐富成果，對于前者主要集中于水分[7]、溫度[8-9]、CO2濃度[10]以及太陽輻射[11]等氣象因素對GPP 時空格局的影響[12]，已通過機理模型[13]、遙感模型[14]及統計方法[15]進行了量化分析與機理解釋。而后者，相關研究也證實了諸如城市建設[16]、退耕還林還草[17]、刀耕火種[18]等人類活動會顯著改變區域GPP 分布[19]，同時土地利用及覆被變化（Land Use and Cover Change，LUCC）作為人類活動的直接體現[20]，是最為經典的表征指標[21]。熱帶地區水熱條件充足，各種環境因素都會對GPP 產生影響，許多研究對此進行了討論，其中在亞馬遜熱帶地區通過機理模型揭示出輻射是生產力變化的主導因素之一[22]；熱帶地區水熱條件雖然相對充足，但時空間的不均勻分配使得年際降水對GPP 變化產生較為關鍵影響[23]，同時GPP 對溫度的變化也有著較高的敏感性[24]。在全球尺度的研究中已有結果表明，溫度、輻射以及水分對于GPP 的影響程度分別為13.07%、-7.24%和11.74%，具體表現為低緯度地區水分對于GPP 變化起到主導作用，中高緯度地區溫度起到主導作用[6]。隨著人類在熱帶地區的活動愈發頻繁[25]，越來越多的研究者發現開荒等農業生產行為致使的土地利用變化才是導致GPP 改變的主要驅動因素，這一觀點得到在亞馬遜雨林、東南亞熱帶雨林等地的相關研究支持[25-26]。綜上所述，熱帶地區植被結構復雜，環境因素多變，雖然科學界對GPP 時空分布格局的驅動因素及響應機理已基本達成共識，但主導因素及其相對貢獻大小依舊為爭論焦點。

海南島是我國第一大熱帶島嶼，作為相對獨立的地理單元是進行熱帶地區GPP 研究的理想切入點；并且近年來隨著海南自由貿易港等政策的落實[27-28]，海南島建設程度加深，人才引進以及城市擴張等一系列活動導致土地覆蓋變化劇烈，為探究人類活動與氣候變化對GPP 的影響提供了契機。因此，本研究選擇海南島作為研究靶區，對其近20 年GPP 動態變化趨勢進行歸因分析，具體目標為：探究2001—2019 年間海南島GPP 時空變化格局；量化氣象要素與人類活動對于年際間GPP 變化的貢獻率，揭示主導因素。旨在為解明熱帶地區生態系統響應區域氣候變化提供科學參考也為海南生態文明試驗區的政策制定提供理論依據。

1 材料與方法

1.1 研究區概況海南島（18°10′—20°10′N，108°37′—111°03′E），地處熱帶北緣，屬熱帶季風海洋性氣候，溫差較低，全年高溫，年平均氣溫22.5～25.6 ℃。年平均降雨約為1 640 mm，降水充足但分配不均，雨旱兩季分明。海南島地勢為四周低平，中間高聳，呈穹隆山地形，中部山區以五指山、鸚哥嶺為隆起核心，向外圍逐級下降，由山地、丘陵、臺地、平原構成環形層狀地貌，梯級結構明顯，土地覆蓋類型豐富，以森林、農田以及草原為主，其中，森林生態系統包括常綠闊葉林、常綠針葉林、落葉闊葉林以及混交林等。

1.2 研究數據

1.2.1 土地利用類型數據海南島土地利用類型數據集來源于MODIS（Moderate-resolution Imaging Spectroradiometer）Land Cover 產品（MCD12Q1），時間分辨率為1 a，空間分辨率為500 m，本研究時間范圍為2001—2019 年[29]。本研究采用了IGBP（International Geosphere Biosphere Programme）分類方案所確定的17 個土地覆蓋方案。結合相關研究[30]以及本研究實際情況，對土地覆蓋類型進行重分類，重分類情況見表1。

表1 根據IGBP 分類系統對土地利用類型的重分類

1.2.2 GPP 數據2001—2019 年海南島GPP 數據集為MODIS 產品MOD17A2HV006。該產品為4 級標準產品[31]，主要通過光能利用率模型進行計算，逐8 天合成空間分辨率為500 m 的GPP 產品，從區域尺度到全球尺度已經被廣泛應用[32]。本研究將GPP 數據求和處理成月值尺度后轉化為年尺度進行后續分析。

1.2.3 氣象數據本研究所使用2001—2019 年逐月Ta、PAR 以及相對濕度（Relative Humidity，RH）氣象數據，其中Ta、RH 數據均源自于國家地球系統科學數據中心（http://www.geodata.cn/），其空間分辨率均為1 km，是利用全國基本氣象站觀測數據進行空間插值并綜合地形數據所繪制的空間分布數據[33-34]。根據Ta 與RH 數據并結合Tetens 經驗公式[35]對海南島逐月VPD(飽和水汽壓差)的數據進行計算；PAR(光合有效輻射)的數據為Global Land Surface Satellite（GLASS）產品（http://www.glass.umd.edu/index.html），其空間分辨率為0.01°，時間分辨率為1 d 是基于多源遙感數據和地面實測數據進一步反演得到的長時間序列、高精度的全球地表遙感產品[36]。

1.3 研究方法

1.3.1 技術路線首先以海南島年際間用地類型未改變的柵格作為數據樣本，利用其多年VPD、Ta 及PAR 年距平值作為解釋變量，GPP 年距平均值作為目標變量，通過高斯過程回歸（Gaussian Process Regression， GPR）與隨機森林（Random Forest，RF）兩種機器學習算法構建預測模型，GPR 是常用的監督分類學習方法，是一種基于貝葉斯方法的非參數概率模型，回歸的目的是通過學習樣本，經過訓練得到輸入變量與輸出變量之間的函數關系，常常用于小樣本回歸分析。RF 是經典的基于分類和回歸樹的集成學習算法，可以解釋若干自變量對因變量的作用；其次，選擇較好的機器學習算法后，利用年際間VPD、PAR 以及Ta 的差異作為輸入數據，對其年際間氣象要素差異所引起的GPP 差異進行估算，完成氣象因素對GPP 影響的相對貢獻率的計算；最后，以年際間用地類型發生改變的柵格作為研究目標，計算相鄰年間LUCC 改變的柵格其GPP 的變化值，進而完成LUCC 對GPP 影響的相對貢獻率的計算；最終進行海南島GPP 變化主導因素的探討（圖1）。

圖1 研究流程圖

1.3.2 機器學習模型構建及評價本研究通過GPR 及RF 機器學習算法構建模型以分析氣象因素對GPP 的相對貢獻率。如表2 所示，海南島近20 年間各用地類型均有不同程度改變，裸地與濕地類型未變柵格所占比例較低，用于訓練模型的樣本數量較少，因此選擇GPR 算法；而林地、草地、耕地等其他類型用地樣本數量充足，利用RF 能夠取得較好效果。所有用地利用的模型訓練均采用五折交叉驗證進行精度的評估，以決定系數（Coefficient Of Determination，R2）、均方根誤差（Root Mean Square Error，RMSE，記作SE）作為模型精度評價指標。

表2 模型評價精度表

1.3.3 歸因分析方法為探究LUCC 以及氣象因素對于海南島GPP 的影響程度，本研究中LUCC以及氣象因素相對貢獻度計算公式如下：

式中，Con(LUCC)(i,i+1)為第i年至第i+1 年間LUCC對于GPP 變化的相對貢獻率，GPPi為i年海南島GPP 總值，GPP(LUCC)(i,i+1)為i至i+1 年間LUCC所驅動的GPP 變化值；Con(Climate)(i,i+1)為i年至第i+1 年間氣象因素對GPP 變化的相對貢獻率，f為1.3.1 所介紹的機器學習模型。VPD(i,i+1)、Ta(i,i+1)以及PAR(i,i+1)分別為第i年至第i+1 年間VPD、Ta 以及PAR 的變化量，其中，i為年份；n為各土地利用類型，n=1，2……7。

2 研究結果

2.1 海南島GPP 變化情況海南島時間上，近20 年期間GPP（以C 表示，下同）時間變化如圖2-a所示，GPP 年際變化明顯，總體上研究時期GPP 呈現上升趨勢（0.44 Tg·a-1）。2005 年海南島GPP 值最低僅為56.46 Tg C，2017 年后GPP 值持續增加，在2019 年達到峰值74.12 Tg C。其GPP 主要由3 種土地利用類型構成，林地為GPP 值最高的土地利用類型，約占到全年GPP 總值的85%，19 年期間GPP 值以0.46 Tg·a-1增加；耕地作為GPP 值僅次于林地的土地利用類型，其年GPP 值最大占比達到10.7%，耕地19 年期間變化幅度不明顯，其年GPP 值以0.07 Tg·a-1增加；草地GPP 值在19 年期間持續降低（-0.11 Tg·a-1）；空間上，近20 年間海南島87.8%的面積呈現出GPP 增加趨勢(圖2-b)；其中，海南島北部以及東北部是GPP 極顯著增加的集中區，中部山區部分地區GPP 無明顯變化，而海口市與三亞市周邊小部分區域的GPP 表現為減少趨勢，且碎片化程度較高，占比為8.9%。

圖2 海南島GPP 時空間變化趨勢

2.2 海南島LUCC 變化情況對近20 年海南島土地利用變化數據進行分析（表3），研究時期內轉移面積小于1 km2的分支不進行展示。結果顯示，2001—2005 年，海南島土地利用轉移總面積為2 736.94 km2。其中，林地的轉出面積最大，凈轉出面積為309.13 km2，主要轉化為耕地；2005—2010 年土地轉移總面積為4 164.04 km2，草地轉化成林地以及耕地轉化為林地為最主要的轉出方式；2010—2015 年期間，海南島的土地覆蓋變化主要發生在林地、草地以及耕地的地區面積分別為1 188.10 km2、1 074.56 km2以及1 070.85 km2；2015—2019 年最主要的轉移方式是林地轉化為耕地和草地轉化為林地，其轉化凈面積分別為541.25 km2和214.65 km2。

表3 海南島近20 年間土地利用每5 年的轉移情況

2.3 LUCC 對GPP 變化的影響研究時期，海南島LUCC 所驅動的GPP 變化主要發生在林地、草地以及耕地地區（圖3）LUCC 所驅動的GPP 年際變化明顯。其中，林地地區LUCC 所驅動的GPP呈增加趨勢，在2018—2019 年林地地區LUCC 所驅動的GPP 達到最大值8.41×10-2Tg C，占比57%；草地地區19 年期間呈顯著微弱增加趨勢，2005—2006 年草地地區達到LUCC 所驅動的GPP峰值7.56×10-2Tg C，2004—2005 年草地地區LUCC 所驅動的GPP 僅為-5.99×10-2Tg C；耕地地區在近20 年內也是呈現出微弱增加趨勢，在2008—2009 年LUCC 所驅動的GPP 值最大。

圖3 各土地利用類型在LUCC 所驅動GPP 變化

2.4 相對貢獻率變化LUCC 以及氣象因素對于海南島GPP 相對貢獻率的年際變化如圖4 所示，年際間LUCC 與氣象要素的相對貢獻率同頻率較高，但程度相差較大，氣象因素主導海南島年際間GPP 變化。LUCC 的相對貢獻率近20 年期間波動較為劇烈，其中2008—2009 年期間相對貢獻率達0.31%，而在2004—2005 年LUCC 對于GPP 的相對貢獻率最小為-0.24%；氣象因素對年際間GPP 變化的相對貢獻率逐年差異明顯，2005—2006 年其正向相對貢獻率達到最大值10.79%，而在2010—2011 年氣象要素對于GPP 的影響程度為歷年來最低值-14.99%。

圖4 相對貢獻率變化圖

3 討 論

3.1 LUCC 對GPP 的影響LUCC 一直是陸生系統GPP 的主要驅動因素。林地、草地以及耕地地區是海南島LUCC 發生的主要區域。研究時期LUCC 共變化15 528.40 km2。其中，植樹造林6 749.33 km2，累積GPP 達到0.23 Tg C，這與相關研究結論[37]基本一致即植樹造林對于GPP 增長有著較為良好的促進作用。這一現象的原因可能是椰林[38]、橡膠林[39]等林種具有較高的經濟價值。此外，近52%的新林地資源主要來源于草地，改善了海南島植被覆蓋，強化了林地生態系統的防線保護作用；耕地也在此過程中凈增加679.41 km2，GPP 在此過程中增加1.63×10-2Tg C。耕地的快速擴張一定程度改變原有的生態環境，由單一經濟作物所構成的生態系統抵抗力穩定性降低可能會在未來引發水土流失，所以如何合理規劃耕地資源是要深入考慮的。

LUCC 對GPP 年際間變化的相對貢獻率遠低于氣象因素，其一在于，海南島林地面積高（79.1%），并可占全島GPP 的84.8%，并且主要集中在中部山區，由于天然林保護等政策[40-41]，中部山區人類活動強度較低，所以在近20 年來，林地并未發生較大的面積變化（近20 年變化不足2%）；二是由于MCD12Q1 產品的將海南島土地利用與覆蓋分為14 類，而為方便分析本研究將其重分類為7 類，不同分類系統所造成的相對貢獻率影響(圖5-a)最高可相差3 倍（2004—2005 年）。尺度效應也是重要因素，如圖5-b 所示，將海南島近20 年逐年逐市縣進行拆解，并使用IGBP 原始分類系統，LUCC 在各市縣中對年際間GPP 變化的相對貢獻率普遍高于全島平均值，LUCC 相較于氣象因素而言作用范圍小，因此在進行全島LUCC 相對貢獻率計算時不同地區間差異相互抵消，因此也造成LUCC 的相對貢獻率較低。LUCC 相較于氣象因素，對于小尺度范圍的GPP年際間變化起到主要作用。LUCC 與氣象因素對于年際間GPP 變化的相對貢獻率具有較高的同頻性，其原因在于相對貢獻率的算法以前一年GPP 為基準，而非年際間GPP 差異為基準，好處在于避免出現年際間GPP 相差極小而導致的極端貢獻；但負效應變為同頻效應以2005—2006 年為例，極端事件第18 號熱帶風暴“達維”對全島GPP 產生顯著負面影響[42]，GPP 為近20 年最低值，因此2005—2006 年會產生極高的LUCC 與氣象因素的正向相對貢獻，各市縣貢獻率受整體變化所影響。

圖5 土地分類對相對貢獻率的影響及各市(縣)相對貢獻率變化

3.2 氣象因素對GPP 的影響本研究僅選取Ta、VPD 以及PAR[43]作為影響GPP 年際間變化的主要氣象驅動因素，未考慮如大氣CO2濃度、氮沉降等環境因素，這為本研究帶來一定不確定性，但根據Sun 等人研究表明[6]，環境要素中氣象因素是環境因素中主導GPP 年際間變化的關鍵因素，因此不確定性能夠最小程度降低。此外，本研究選取VPD 而非降水作為水分條件，主要由于像海南島這類低緯度濕潤地區，大氣的濕潤程度是限制植物生長的主要因素[44]。其次，對于氣象因素間的交互作用，本研究所采用機器學習算法（GPR 與RF），屬于數據驅動模式，能夠減弱自變量間的交互關系[15]。

氣象因素相對貢獻較高，分別選取正負相對貢獻率最高的2008—2009 年（2005—2006 年為最高正相對貢獻，但主要受極端事件“達維”影響，因此選取2008—2009 年）與2010—2011 年為例；2008—2009 年，PAR、Ta 及VPD 均呈現超距平的轉好態勢；而2010—2011 年PAR 與Ta 呈現出遠超平均年際間變化的降幅，同時水分制約也略有所增加（圖6）；可見本研究所提出的基于機器學習的數據驅動算法能夠捕捉氣象要素對GPP 的影響。除氣象要素與LUCC 外，人類的管理與經營活動，如耕地農業管理及政策支持[45-46]、橡膠林下經營[47]、天然林的保護[48]、極端事件以及病蟲害等異常事件也會對海南島GPP 產生影響[49]。

圖6 海南島各氣象要素變化圖

3.3 不確定性分析首先本研究所使用的GPP與LUCC 數據來源于MODIS 產品，MODIS 作為光學遙感產品應用于熱帶地區時不可避免精度會降低[50]，這是傳統光學遙感產品共同面臨問題。MODIS 的LUCC 產品在全球尺度的精度為73.6%[51]，本研究將可能出現異物同譜的相似類進行了重分類，共7 大類，這能夠大大降低其分類精度所帶來的不確定性(圖5-a)，這里不能夠忽視。其次，MODIS 的GPP 產品基于光能利用率算法，是現階段主流大區域尺度的遙感GPP 算法，包括EC-LUE[52]、CASA 以及VPM[6]等多種形式，其中可吸收光合有效輻射比率（fAPAR，fraction of Absorbed PAR）是模型核心部分，無論何種形式模型均將fAPAR 考慮為關于NDVI 等植被指數的函數；長時間序列植被指數主要源于MODIS 產品，因此本研究中直接使用了MODIS 的GPP 產品，避免二次計算時增加不確定性。此外關于計算方案引入的不確定性主要為，當年際間LUCC發生改變時，對GPP 的影響則全部歸于LUCC 的貢獻；但實則氣象要素依舊對該區域產生影響，可此部分的影響并未考慮到本研究中，低估了氣象因素造成的總體影響。綜上，鑒于本研究的時間范圍與尺度，結合各遙感產品的精度，最終確定MCD12Q1 為研究數據，該數據產品在國內的研究已廣泛應用于國內相關研究[53]，并取得較為可靠的成果。考慮到海南島全年較高的云覆蓋，產品精度有限。未來要進一步研究LUCC 對GPP 影響等問題時，更高精確度土地利用數據產品是不可或缺的；本研究使用機器學習模型雖然可以在一定程度上將氣象因素與人類活動進行解耦并單獨分析其對GPP 的影響，實際環境中GPP 的影響是多元驅動的，無法對影響因素的共同作用進行驅動分析。本研究中假設LUCC 發生改變時GPP的變化歸因為LUCC 所引起，將氣象因素作用歸于LUCC 貢獻中，這也帶來了一定不確定性。