基于隨機(jī)森林模型的黃土高原草地凈初級(jí)生產(chǎn)力時(shí)空格局及未來演變趨勢(shì)模擬

劉歡歡， 陳 印， 劉 悅， 剛成誠(chéng)

（1.西北農(nóng)林科技大學(xué)草業(yè)與草原學(xué)院，陜西 楊凌 712100；2.西北農(nóng)林科技大學(xué)水土保持研究所，陜西 楊凌 712100；3.中國(guó)科學(xué)院水利部水土保持研究所，陜西 楊凌 712100）

草地是陸地生態(tài)系統(tǒng)主要組成部分之一，約占陸地面積30%［1］。草地生態(tài)系統(tǒng)在維持生態(tài)系統(tǒng)功能以及服務(wù)社會(huì)經(jīng)濟(jì)中均發(fā)揮重要功能，包括保護(hù)生物多樣性、水土保持以及減緩溫室效應(yīng)等［2］。植被凈初級(jí)生產(chǎn)力是衡量生態(tài)系統(tǒng)固碳能力的重要指標(biāo)之一［3］。草地凈初級(jí)生產(chǎn)力（Net Primary Pro?duction，NPP）不僅能夠直接反映草地生態(tài)系統(tǒng)生產(chǎn)力的高低，而且是評(píng)價(jià)草地生態(tài)系統(tǒng)碳收支及調(diào)節(jié)草地生態(tài)過程的主要因子［4］。因此，科學(xué)評(píng)估草地NPP及其時(shí)空動(dòng)態(tài)對(duì)草地生態(tài)系統(tǒng)管理及草地資源保護(hù)至關(guān)重要。

傳統(tǒng)野外采樣估算草地NPP 的方法雖然結(jié)果較為準(zhǔn)確，但由于其需要耗費(fèi)大量的人力、物力和財(cái)力，在大空間尺度草地NPP 的估算中較少使用。遙感技術(shù)和生態(tài)模型的發(fā)展為此提供了重要的方法和技術(shù)手段。遙感技術(shù)通過利用遙感指標(biāo)數(shù)據(jù)與草地生物量的相關(guān)關(guān)系來估算草地NPP，如馬文紅等［5］基于NDVI 以及地上生物量數(shù)據(jù)模擬了1982—2006 年中國(guó)北方草地生物量空間分布及動(dòng)態(tài)；王鶯等［6］利用EVI 與草地干物質(zhì)產(chǎn)量回歸模擬了2006—2008 年甘南草地NPP 變化。模型模擬同樣是大空間尺度草地NPP 估算的常用方法，如Guo等［7］利用CASA模型模擬了2001—2018年內(nèi)蒙古草地NPP；張赟鑫等［8］利用CASA模型和Miami 模型模擬了中亞草地NPP；Zhang 等［9］利用BIOME-BGC 模型模擬了呼倫貝爾草原NPP。近年來機(jī)器學(xué)習(xí)方法在植被生產(chǎn)力估算中亦有廣泛的應(yīng)用。Jung 等［10］利用機(jī)器學(xué)習(xí)中模型樹集成算法模擬了全球總初級(jí)生產(chǎn)力；Yao 等［11］利用模型樹集成算法模擬了1982—2015 年中國(guó)0.1°空間分辨率的總初級(jí)生產(chǎn)力。隨機(jī)森林（Random Forest，RF）是基于決策樹的集成機(jī)器學(xué)習(xí)算法，目前已被應(yīng)用于全球［12］、我國(guó)西北地區(qū)［13］及青藏高原［14］不同類型草地NPP的估算。

黃土高原是世界上最大的黃土覆蓋區(qū)，氣候干旱，水土流失嚴(yán)重，生態(tài)環(huán)境極其脆弱［15-16］。隨著退耕還林還草工程等生態(tài)工程的實(shí)施，黃土高原植被變化及其環(huán)境效應(yīng)備受關(guān)注［17］。目前，在黃土高原草地NPP 時(shí)空動(dòng)態(tài)及驅(qū)動(dòng)力分析方面已開展了一定的研究［18-19］。但由于黃土高原地區(qū)地形復(fù)雜，不同研究結(jié)果之間存在較大差異。在全球變暖背景下，黃土高原草地NPP 未來的演變趨勢(shì)仍不清楚。因此，本研究基于氣候、植被、土壤和地形因子，利用隨機(jī)森林模型模擬了2002—2020 年黃土高原草地NPP 的時(shí)空動(dòng)態(tài)，并基于共享社會(huì)經(jīng)濟(jì)路徑（Shared Socioeconomic Pathway，SSPs）4 個(gè)未來氣候情景數(shù)據(jù)，估算了黃土高原草地NPP的未來演變趨勢(shì)，以期為黃土高原草地生態(tài)系統(tǒng)管理及可持續(xù)發(fā)展提供科學(xué)依據(jù)。

1 研究區(qū)概況與方法

1.1 研究區(qū)概況

黃土高原地處我國(guó)黃河流域中部，地理位置位于100°54′～114°33′E，33°43′～41°16′N，總面積達(dá)63×104km2［20］（圖1），橫跨陜西、甘肅、寧夏、山西、青海、河南以及內(nèi)蒙古7 個(gè)省份，是我國(guó)四大高原之一，同時(shí)也是世界上最大的黃土覆蓋區(qū)［15］。黃土高原包括了干旱、半干旱以及半濕潤(rùn)氣候區(qū)，年降雨量在150～800 mm，海拔在200～3000 m，不同區(qū)域在氣候上有著明顯差異［21］。黃土高原的植被類型以森林、草地以及灌木地為主，草地面積約占總面積60%［22］。

圖1 研究區(qū)概況及采樣點(diǎn)分布Fig.1 Study area and sampling points

1.2 數(shù)據(jù)獲取及處理

通過文獻(xiàn)、國(guó)家地球系統(tǒng)科學(xué)數(shù)據(jù)中心數(shù)據(jù)以及GEDI生物量數(shù)據(jù)（https：//daac.ornl.gov/）獲取黃土高原草地地上生物量數(shù)據(jù)1788個(gè)，根據(jù)根冠比估算地下生物量，使用草地生物量乘以轉(zhuǎn)換系數(shù)（0.475）估算草地NPP（g C·m-2），作為實(shí)測(cè)NPP數(shù)據(jù)［23］。

為訓(xùn)練模型，選取19個(gè)特征參數(shù)包括土壤理化因子（S）、地形因子（T）、氣象因子（A）和植被因子（B）4 類。氣象因子中溫度和降水?dāng)?shù)據(jù)來源于國(guó)家地球系統(tǒng)科學(xué)數(shù)據(jù)中心-黃土高原分中心（http://loess.geodata.cn），根據(jù)月平均溫度數(shù)據(jù)合成年平均溫度（TEM）數(shù)據(jù)以及年生長(zhǎng)季平均溫度（TEM4-10）數(shù)據(jù)，根據(jù)月最低溫度數(shù)據(jù)合成為年最低溫度（TMN）數(shù)據(jù)，根據(jù)逐月降水量數(shù)據(jù)合成年降水量（PRE）數(shù)據(jù)以及年生長(zhǎng)季降水量（PRE4-10）數(shù)據(jù)［24］。數(shù)字高程模型（DEM）數(shù)據(jù)來自NASA 的SRTM 數(shù)據(jù)（http://srtm.csi.cgiar.org/），坡度（Slope）數(shù)據(jù)由DEM 數(shù)據(jù)通過ArcGIS 計(jì)算得來。土壤理化因子數(shù)據(jù)均來自世界土壤數(shù)據(jù)庫（https://data.isric.org/）。日光誘導(dǎo)葉綠素?zé)晒猓⊿IF）數(shù)據(jù)來自GOSIFv2 數(shù)據(jù)［25］（https://globalecology.unh.edu），植被有效光合輻射吸收比例（FPAR）數(shù)據(jù)、歸一化植被指數(shù)（NDVI）以及蒸散量（ET）來自MODIS 數(shù)據(jù)產(chǎn)品（https://modis.gsfc.nasa.gov/）。未來降雨量和溫度數(shù)據(jù)來源于WorldClim（https://www.worldclim.org/）。黃土高原空間范圍數(shù)據(jù)以及黃土高原草地分布數(shù)據(jù)來自中國(guó)科學(xué)院資源環(huán)境與數(shù)據(jù)中心（https://www.resdc.cn/）空間分辨率為1 km。以上數(shù)據(jù)分辨率均通過ArcGIS 重采樣為1 km，采用D_WGS_1984 投影。本文所有環(huán)境因子見表1。

表1 隨機(jī)森林（RF）中使用的環(huán)境因子Tab.1 Environmental factors in RF

SSPs 情景在典型濃度路徑（RCPs）情景基礎(chǔ)上發(fā)展而來，用于定量描述氣候變化與社會(huì)經(jīng)濟(jì)發(fā)展路徑之間的關(guān)系，反映未來社會(huì)面臨的氣候變化適應(yīng)和減緩挑戰(zhàn)［26］。本研究分別選擇SSP126、SSP245、SSP370 和SSP585 4 個(gè) 情 景 下，ACCESSESM1-5、BCC-CSM2-MR、CanESM5、CMCC-ESM2和CNRM-CM6-1，共計(jì)5 個(gè)氣候系統(tǒng)模式2021—2040年（2030s）、2041—2060 年（2050s）、2061—2080 年（2070s）、2081—2100 年（2090s）4 個(gè)時(shí)間段的降水量和溫度作為輸入數(shù)據(jù)，利用RF 模型預(yù)測(cè)未來不同氣候情景下黃土高原草地NPP 的演變趨勢(shì)。RF模型中的其他環(huán)境因子仍使用2020年數(shù)據(jù)。

1.3 研究方法

1.3.1 建模方法 本研究采用RF 進(jìn)行建模，RF 是基于分類和回歸樹（Classification and Regression Tree，CART）而開發(fā)的集成算法，可用于分類問題以及回歸問題［27-28］。RF 的主要原理為通過結(jié)合多個(gè)決策樹，并平均其結(jié)果使得決策樹泛化誤差收斂從而產(chǎn)生更好的預(yù)測(cè)結(jié)果［29］。RF 可廣泛用于預(yù)測(cè)問題，并且應(yīng)用簡(jiǎn)單，無需復(fù)雜地調(diào)整參數(shù)，通過調(diào)整RF 中回歸樹的數(shù)量以及決策樹每處節(jié)點(diǎn)預(yù)測(cè)變量的數(shù)量來優(yōu)化模型，提高預(yù)測(cè)精度［29］。此外，RF 具有處理高維特征空間數(shù)據(jù)的能力，在遙感應(yīng)用中已有較多研究［30］。RF 模型性能通過決策系數(shù)（R2）以及均方根誤差（RMSE）評(píng)價(jià)其擬合能力以及預(yù)測(cè)性能。本研究使用R 4.0.5進(jìn)行隨機(jī)森林建模、參數(shù)優(yōu)化、預(yù)測(cè)以及空間可視化。

1.3.2 Theil-Sen Median趨勢(shì)分析和Mann-Kendall檢驗(yàn) Theil-Sen Median 趨勢(shì)分析是一種穩(wěn)健的非參數(shù)統(tǒng)計(jì)趨勢(shì)計(jì)算方法［31］，計(jì)算公式如下：

式中：xj和xi為樣本數(shù)據(jù)值。當(dāng)Senslope＆gt;0 時(shí)，表明該數(shù)據(jù)在此時(shí)間序列呈增長(zhǎng)趨勢(shì)，反之則表示呈下降趨勢(shì)。

Mann-Kendall檢驗(yàn)方法是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)方法，用于判斷趨勢(shì)的顯著性［30］。計(jì)算公式如下：

式中：S為Mann-Kendall 檢驗(yàn)需要構(gòu)建的統(tǒng)計(jì)量，ZMK是標(biāo)準(zhǔn)正態(tài)變量。本文在置信度水平α=0.05上來檢驗(yàn)NPP變化趨勢(shì)的顯著性。當(dāng)P小于0.05時(shí)表示NPP變化顯著，反之則表示不顯著。

2 結(jié)果與分析

2.1 NPP與環(huán)境因子相關(guān)性分析

在植被因子中，NPP 與SIF 相關(guān)性最強(qiáng)（r=0.82），其次為NDVI 和FPAR（圖2a）。在氣象因子中，NPP與ET相關(guān)性最高（r=0.91），與PRE4-10、PRE、TMN 相關(guān)性也較強(qiáng)，而與TEM 以及TEM4-10相關(guān)性較弱。NPP 與土壤因子整體相關(guān)性較弱，其中NPP與SAND_S 相關(guān)性最強(qiáng)（r=-0.27），與SOC_T 相關(guān)性最弱（r=0.11）。整體而言，NPP 與氣象因子、植被因子相關(guān)性較強(qiáng)，而與氣象因子中的TEM 和TEM4-10、土壤因子以及地形因子相關(guān)性弱。而在隨機(jī)森林模型因子重要性中（圖2b），對(duì)模型估算NPP影響最大的是ET。PRE、PRE4-10以及植被因子中的SIF 對(duì)NPP估算影響也較大。結(jié)合相關(guān)性分析和重要性排序，影響黃土高原草地NPP的環(huán)境因子主要為植被因子以及PRE、PRE4-10和ET。

圖2 NPP與環(huán)境因子相關(guān)性（a）及隨機(jī)森林因子重要性（b）Fig.2 The correlations between NPP and environmental factors(a),and the importance of random forest factors(b)

2.2 RF模型精度評(píng)價(jià)

2.3 草地NPP空間動(dòng)態(tài)

黃土高原草地NPP 的空間分布呈現(xiàn)較強(qiáng)的異質(zhì)性（圖4）。2002—2020 年黃土高原草地年平均NPP 為276.55 g C·m-2。黃土高原草地NPP 的分布呈現(xiàn)從南向北逐漸遞減的趨勢(shì)，其中黃土高原西北地區(qū)草地NPP 最低，主要位于陜西北部、寧夏以及甘肅蘭州地區(qū)附近，最低值為76.52 g C·m-2；草地NPP 最高的區(qū)域主要分布在黃土高原南部地區(qū)，如陜西關(guān)中地區(qū)，最高值為634.86 g C·m-2；而在黃土高原東部地區(qū)草地分布較少，草地NPP 約為450 g C·m-2。

圖4 2002—2020年黃土高原草地NPP空間分布Fig.4 Spatial pattern of grassland NPP in the Loess Plateau during 2002-2020

2002—2020年，黃土高原草地NPP呈現(xiàn)總體增加的趨勢(shì)（圖5），其中16.59%的地區(qū)呈顯著增加趨勢(shì)，38.52%的地區(qū)呈輕微增加趨勢(shì)，主要分布在黃土高原東部以及南部；1.89%的地區(qū)草地NPP 變化不顯著；35.57%的地區(qū)呈輕微減少趨勢(shì)，7.43%的地區(qū)呈顯著減少趨勢(shì)，主要分布在黃土高原西北地區(qū)。

圖5 2002—2020年黃土高原草地NPP空間變化Fig.5 Spatial dynamic of grassland NPP in the Loess Plateau during 2002-2020

圖3 隨機(jī)森林模型精度驗(yàn)證Fig.3 Accuracy test of Random Forest model

2.4 未來不同氣候情景下黃土高原地區(qū)草地NPP演變趨勢(shì)

到本世紀(jì)末，4 個(gè)未來氣候情景下黃土高原草地NPP 均呈現(xiàn)增加的趨勢(shì)，其中SSP126 情景草地NPP增加幅度最少，為1.20 g C·m-2·（10a）-1；SSP585情景增幅最大，為8.10 g C·m-2·（10a）-1。SSP245 和SSP370 情景，草地NPP 的增速分別為3.90 g C·m-2·（10a）-1和4.00 g C·m-2·（10a）-1（圖6）。草地NPP 增加的區(qū)域主要集中在黃土高原南部地區(qū)，SSP585中草地NPP 在黃土高原南部地區(qū)增加最多。而草地NPP 下降的區(qū)域主要集中在黃土高原西部地區(qū)，其中在SSP126中草地NPP將下降最多。

圖6 不同氣候情景下2090s黃土高原草地NPP空間分布Fig.6 Spatial pattern of grassland NPP in the Loess Plateau in 2090s under different climate scenarios

3 討論

3.1 隨機(jī)森林模型估算草地NPP適用性

相比其他NPP 估算模型，隨機(jī)森林的優(yōu)勢(shì)在于可以更準(zhǔn)確地描述NPP與其他環(huán)境因子之間的非線性復(fù)雜關(guān)系，且隨機(jī)森林有較好抗過擬合能力［32-33］。本研究利用隨機(jī)森林模擬得到黃土高原草地的平均NPP（276.55 g C·m-2）。劉洋洋等基于CASA模型模擬的2000—2015 年黃土高原草地年均NPP 為202.93 g C·m-2，原因可能與使用草地分布數(shù)據(jù)不同有關(guān)［19］。魏建洲基于MOD17A3 NPP 數(shù)據(jù)模擬的黃土高原草地NPP 在208.04～262.66 g C·m-2之間，其高值與本研究結(jié)果較為接近［34］。Li等利用CASA模型分別模擬了石羊河流域及我國(guó)西北干旱半干旱區(qū)草地年均NPP，其結(jié)果分別為236.99 g C·m-2和277.97 g C·m-2［35-37］。相比于其他研究結(jié)果，本研究模擬草地NPP結(jié)果偏大，一方面在于不同研究中所使用的草地分布數(shù)據(jù)不同，如劉洋洋使用的草地分布數(shù)據(jù)來自Global Land Cover 2000 數(shù)據(jù)產(chǎn)品［19］，本研究所使用草地分布數(shù)據(jù)來源于中科院土地利用數(shù)據(jù)產(chǎn)品；另一方面基于隨機(jī)森林模型模擬草地NPP，受訓(xùn)練數(shù)據(jù)影響很大，本研究用于訓(xùn)練模型的實(shí)測(cè)草地NPP數(shù)據(jù)主要集中在黃土高原西部，東部地區(qū)的樣本點(diǎn)數(shù)量相對(duì)較少。但總體而言，本研究利用隨機(jī)森林模擬的黃土高原草地NPP 結(jié)果較為可信。

3.2 黃土高原草地NPP的動(dòng)態(tài)變化

本研究發(fā)現(xiàn)2002—2020年黃土高原草地NPP呈增加趨勢(shì)，這與前人在該區(qū)域研究結(jié)果一致［19,38-40］。由于黃土高原氣候干旱，降水是主要的植被生長(zhǎng)限制因子［18］。2002—2020 年黃土高原降水量呈增加趨勢(shì)（圖7），促進(jìn)了黃土高原草地NPP 的增加。劉洋洋等［19］也認(rèn)為黃土高原降水增加導(dǎo)致的濕潤(rùn)環(huán)境為草地NPP 增加提供了有利條件。Gang 等［40］的研究結(jié)果同樣表明氣候變化對(duì)黃土高原草地NPP增加有著促進(jìn)作用。從空間變化上來看，黃土高原西部和南部草地NPP增加，與劉錚等［38］的研究結(jié)果相似。這些地區(qū)位于半干旱以及半濕潤(rùn)地區(qū)，降水量的增加導(dǎo)致了NPP的快速積累［41］。此外，自1999年以來，退耕還林還草工程通過圍欄封育、封山禁牧等措施使草地生態(tài)環(huán)境得到了較大改善，是導(dǎo)致黃土高原草地NPP增加的重要因素［40-42］。到本世紀(jì)末，4個(gè)情景中黃土高原草地NPP仍呈增加趨勢(shì)，原因在于未來黃土高原地區(qū)降水量和最低溫度將繼續(xù)增加，暖濕化的氣候條件有利于草地NPP 的積累，尤其是在SSP585情景中。

圖7 2002—2020年黃土高原年降水量變化趨勢(shì)Fig.7 Temporal trend of mean annual precipitation in the Loess Plateau from 2002 to 2020

3.3 RF估算結(jié)果的不確定性分析

本研究使用RF 模型模擬黃土高原草地NPP 精度較高，但仍存在一定的不確定性。首先，RF 模型性能受數(shù)據(jù)集數(shù)量及質(zhì)量影響較大［43］，黃土高原區(qū)域地形復(fù)雜，本文使用1 km的空間分辨率數(shù)據(jù)模擬草地NPP 導(dǎo)致部分區(qū)域草地NPP 模擬結(jié)果精度較低。其次，RF 模型作為經(jīng)驗(yàn)?zāi)Ｐ?，受?xùn)練數(shù)據(jù)質(zhì)量以及數(shù)量影響較大，本研究草地NPP實(shí)測(cè)數(shù)據(jù)空間分布不均勻，在黃土高原西北地區(qū)以及西南地區(qū)較多，而其余地方較少，未來應(yīng)考慮使用更多的實(shí)測(cè)數(shù)據(jù)進(jìn)行模擬。此外，RF 模型屬于黑箱模型，通過多次調(diào)試來確定合理模擬結(jié)果，但無法控制模型內(nèi)部運(yùn)行［44］。

4 結(jié)論

本研究基于野外實(shí)測(cè)數(shù)據(jù)，以植被、氣候、土壤和地形因子作為環(huán)境因子，利用隨機(jī)森林模型模擬了黃土高原草地NPP的時(shí)空動(dòng)態(tài)，并估算了未來不同氣候情景下，草地NPP 的演變趨勢(shì)，得出了以下主要結(jié)論：

（1）2002—2020 年，黃土高原草地平均NPP 為276.55 g C·m-2·a，空間上呈“東南高西北低”的分布特征。近19 a 黃土高原草地NPP 整體呈現(xiàn)上升趨勢(shì)，主要得益于20世紀(jì)以來黃土高原地區(qū)降雨量的增加及“退耕還林還草”等一系列生態(tài)工程。

（2）到本世紀(jì)末，在暖濕化的氣候背景下，黃土高原草地NPP 均將繼續(xù)保持增加的趨勢(shì)，其中在SSP585情景中將增加最多。

（3）基于氣候、環(huán)境、土壤和地形因子構(gòu)建了黃土高原草地NPP隨機(jī)森林模型，尚未包含人為活動(dòng)指標(biāo)因子，構(gòu)建綜合氣候、環(huán)境、土壤和人為因素等多因子的要素體系，將會(huì)進(jìn)一步提高草地植被生產(chǎn)力的估算精度。