造林后表層土壤有機碳變化及影響因素分析

王 藝 杰，于 麗 君，張 穩，孫 文 娟*

(1.中國科學院植物研究所植被與環境變化國家重點實驗室，北京 100093;2.中國科學院大學，北京 100049;3.中國科學院大氣物理研究所大氣邊界層物理和大氣化學國家重點實驗室,北京100029)

0 引言

大氣CO2濃度升高[1]是全球變化的重要問題，樹木通過光合作用固碳是緩解全球CO2濃度上升的有效策略之一[2],因此，造林和再造林被認為是增加陸地碳固存、緩解大氣CO2濃度上升的有效管理方式[2]。陸地生態系統的大部分碳儲存在土壤中，土壤碳儲量大于植被和大氣碳儲量之和[3]，其碳輸入與輸出平衡對于大氣CO2濃度有重要影響[4],尤其是土壤有機碳比無機碳參與地球生物化學循環的活躍性更高[5]，因此，合理評估大范圍造林后土壤有機碳密度變化，可為國家實施森林管理、規劃減排計劃提供可靠的數據支持，且有助于國際社會明確人工林土壤有機碳庫在抵消碳排放中的作用。

由于缺乏可靠的造林后土壤有機碳密度估算方法，目前主要采用測定數據外推或基于土壤有機碳密度變化與環境要素關系的簡單假設進行估算。例如：方精云等基于土壤固碳占整個人工林生態系統固碳量約1/3的比例，結合森林清查數據折算出1981-2000年中國土壤碳匯[6]；Nave等將時間作為唯一自變量估算美國造林后0～10 cm的土壤有機碳密度變化量[7]。這些簡單估算忽略了影響造林后土壤變化的諸多因素及其相互作用的復雜性，導致結果不確定性較大。有研究指出，造林后土壤有機碳密度變化與造林年限存在典型的時間函數關系[8,9]。Nave等進一步提出，如果土壤有機碳密度與植被的結構、生產力以及氣候和土壤的物理特性等因素存在定量響應關系，則可以通過植被生產力等數據估計土壤有機碳密度[10]。因此，本文通過收集全球造林后土壤有機碳變化的測定數據，分析影響人工林表層(0～20 cm)土壤有機碳密度變化的因素，并建立人工林表層土壤有機碳指標與植被、氣候、土壤和土地利用因素關系的定量模型，以此估算我國不同區域造林后土壤有機碳變化。

1 數據及研究方法

1.1 數據來源與篩選

本文2009-2020年全球造林數據來源于Web of Science和中國知網等文獻數據庫，2009年前數據根據已發表的meta分析或綜述文獻[8,9,11-15]獲得，對初步獲取的340篇中英文文獻進一步篩選：文獻記述中包含造林前土地利用類型、土壤容重、造林樹種、林齡信息，土壤采樣深度大于20 cm，且必須包含對照組；若文獻中采用時空代替法進行長時序研究，則將造林年限小于3年的土壤作為對照組。經過篩選，共有104篇文獻滿足條件。

樣點數據在土壤、氣候環境、樹種和造林年限等方面的代表性，是建立區域適用的造林后土壤有機碳指標估算模型的必要前提。經上述篩選后的數據共包含123個研究樣點(圖1)，390組成對數據，造林年限跨度從6個月到91年，包含針葉樹23種(油松、馬尾松、杉木、赤松、黑松、輻射松、濕地松、云杉、冷杉等)、闊葉樹31種(楊屬、桉屬、合歡屬、櫟屬等)。樣點的土壤質地包括粉壤土(粉粒含量<60%)、黏土(黏粒含量>60%)、砂土(砂粒含量>90%)等，涵蓋從低緯度熱帶到中高緯度寒溫帶等主要陸地生態系統氣候帶，造林前土壤有機碳含量(有機碳重量百分比)范圍為0.03%～6.15%，整理和計算的土壤數據均為0～20 cm深度。

注：審圖號為GS (2016)1665號，底圖無修改。

1.2 自變量選取與缺失值補充

依據已有研究[10-12,16]，本文先驗性地選取植被類型、造林前土地利用類型、造林年限、喬木地上/地下有機碳密度、土壤初始有機碳密度、氣候濕潤度、土壤pH值、土壤黏粒+粉粒含量(clay+silt)作為自變量，對土壤有機碳變化進行分析。其中植被類型、造林前土地利用類型和造林年限信息在文獻的樣點描述中均有記述，個別文獻中未明確記述的自變量通過以下公式進行填補。

ABC=SAB×TD×0.5/1 000

(1)

UBC=SUB×TD×0.5/1 000

(2)

式中：ABC、UBC分別為喬木地上、地下有機碳密度(Mg C hm-2)；SAB、SUB分別為單株喬木地上、地下生物量(kg/株)；TD為喬木密度 (株/hm2)，利用樹木異速生長方程計算[17]；常數0.5為林木生物量含碳量[18]。

SOCDbf=BD×D×SOC

(3)

(4)

式中：SOCDbf為土壤初始有機碳密度(Mg C hm-2)；BD為土壤容重(g/cm3)，缺失值通過Adams方程計算[19]；D為土壤深度(cm)；SOC為土壤有機碳含量(%)。

IdM=P/(T+10)

(5)

式中：IdM為氣候濕潤度[20]；P為年降水量 (mm),T為年均氣溫(℃),二者缺失值采用AgMERRA樣點區的1980-2010年氣候平均值代替[21]。

本研究從GSDE(Global Soil Dataset for use in Earth System Models)數據庫[22]獲取樣點坐標對應的土壤pH值及其機械組成數據，為使GSDE數據庫中的土壤深度分層與本文關注的土壤深度一致，需對GSDE數據進行土壤深度轉化，計算公式為：

(6)

(7)

Material_content0～xi=Material_massxi/Soil_massxi

(8)

式中：Soil_mass為單位面積的土壤質量；Material_mass為單位面積的物質質量；Material_content為物質百分含量(%)；0～xi表示從0 cm到xi深度的值，i為土壤層序數。

1.3 土壤有機碳指標計算

本文表征土壤有機碳變化的指標包括土壤有機碳密度占生態系統有機碳密度的比例(PSOCD)和土壤有機碳密度變化量(ΔSOCD)。因不同土地利用方式下土壤容重不同，相同深度下土壤質量也不同[23]，因此需要先進行土壤等質量法校正(Equivalent Soil Mass，ESM)，使造林前后的土壤在相同質量下進行有機碳密度對比[9,12,23]，計算公式為：

SOCafn×100

(9)

ΔSOCD=SOCDcorr-SOCDbf

(10)

PSOCD=SOCDcorr/(SOCDcorr+ABC+UBC)

(11)

式中：SOCDcorr表示校正后的土壤有機碳密度(Mg C hm-2)；SOCDaf表示造林后土壤有機碳密度(Mg C hm-2)，其計算同式(3)；i為土壤層序數；n為最底層土壤序數。

1.4 模型建立與退耕還林后土壤固碳情景模擬

研究全球造林后土壤有機碳指標變化的文章以meta分析為主[8,11]，但該方法只能研究分類變量，當變量是數值型(特別是連續型)數值時，需先劃分為分類變量[8,11]，該過程會造成大量信息損失，且無法定量分析影響因素的總效應。本文根據前人研究結果，選取對土壤有機碳指標有顯著或潛在影響的自變量[8,10-13]，采用逐步線性回歸模型，選擇AIC(Akaike′s Information Criterion)最小的方程形式，分析其對土壤碳指標的影響。鑒于回歸模型易受數據極值影響，要求分析數據滿足方差齊性和正態性等[24]，本文同時使用隨機森林模型[25]分析土壤有機碳指標變化的影響因素。首先將單個決策樹的特征空間最大值設為3(即自變量總數的1/3)，根據遍歷調優的結果，將隨機森林建立子樹的數量設為1 100；然后使用交叉檢驗方法，以均方根誤差(RMSE)(式(12))檢驗模擬值和觀測值之間的偏差[26,27]，選擇相應變量代入隨機森林模型。

(12)

式中：Pi和Oi分別為樣本i的模擬值和觀測值；n為樣本量。

退耕還林后土壤有機碳指標變化的情景模擬使用全國尺度的土壤理化數據、年均溫和年降水數據，造林區域、造林密度和樹種數據來源于中國《造林技術規程》《主要樹種齡級與齡組劃分》[28,29]。中國《造林技術規程》主要以年積溫和年降水量為標準劃分造林氣候區[28]，得到9個造林氣候區，除極干旱區不適宜喬木生長外，其余造林氣候區分別為寒溫帶區、中溫帶區、暖溫帶區、亞熱帶區、熱帶區、半干旱區、干旱區和高寒區。

2 結果分析

2.1 土壤有機碳指標變化影響因素定量分析

逐步線性回歸模型分別定量描述了PSOCD和ΔSOCD與植被、氣候、土壤和土地利用因素的關系，如表1所示:對PSOCD有正向作用的自變量為土壤黏粒+粉粒含量(P<0.001)與土壤初始有機碳密度(P<0.001)，其回歸系數分別為0.15和0.36；有負向作用的為土壤pH值(P<0.01)、喬木地上(P<0.001)和地下(P<0.05)有機碳密度、造林年限(P<0.01)，其系數分別為-1.70、-0.28、-0.33和-0.15；造林前土地利用類型對PSOCD的影響排序為農田>草地>荒地>林地(P<0.001)。對ΔSOCD有正向作用的是土壤黏粒+粉粒含量(P<0.001)、喬木地下有機碳密度(P<0.01)、造林年限(P<0.01)和氣候濕潤度(P<0.05)，其回歸系數分別為0.26、0.36、0.16和0.07；有負向作用的為土壤初始有機碳密度(P<0.001)、土壤pH值(P<0.1)和喬木地上有機碳密度(P<0.05)，其系數分別為-0.15、-0.87和-0.07；造林前土地利用類型對ΔSOCD的影響排序為農田>荒地>草地>林地(P<0.05)，植被類型對其影響排序為針闊混交林>常綠闊葉林>常綠針葉林>落葉闊葉林>落葉針葉林(P<0.05)。

在隨機森林模型中，自變量與因變量的關系沒有解析表達式，僅能通過偏相關曲線進行判斷，因此將隨機森林模型的偏相關曲線與逐步線性回歸各自變量對因變量的作用方向相結合(圖2、圖3)，便于理解隨機森林模型的解析效果。由圖2可知，隨機森林模型與逐步線性回歸模型結果基本一致，即PSOCD隨土壤初始有機碳密度和土壤黏粒+粉粒含量增加而增加，隨喬木地上/地下有機碳密度、土壤pH值和造林年限增加而減少。其中，土壤黏粒+粉粒含量高有利于土壤形成團聚體，吸附并保護碳，因此有利于土壤有機碳密度的增加[11,12]。土壤pH值影響樹木生長和土壤碳過程，低pH值會降低樹木的生長速度，從而減少土壤的碳輸入，同時會降低土壤有機質的分解速率，有利于土壤有機碳積累[30]；反之，高pH值地區樹木生長速度更快，土壤有機碳分解也更快，綜合結果是堿性土壤的PSOCD隨著pH值增加迅速降低(圖2c)。PSOCD隨造林年限增加而下降，主要是由于土壤有機碳密度的增加速率小于喬木有機碳密度的增加速率。已有研究通常將氣候帶作為分類變量分析溫度和降水對土壤有機碳的影響[12,15]，認為從溫帶地區到亞熱帶地區，土壤有機碳密度有增大趨勢[11,12,19]，這與本研究結果一致(圖2d)。

注：正負號表示變量在逐步線性回歸模型中的符號方向，未出現則表示變量在逐步線性回歸中被剔除， *、**、***分別表示P<0.05、P<0.01、P<0.001，下同。

圖3 不同自變量與 ΔSOCD的隨機森林模型偏相關曲線及其在逐步線性回歸模型中的作用方向Fig.3 Partial correlation curves of ΔSOCD and independent variables based on random forest model and the action direction of independent variables in stepwise linear regression model

由圖3可知，除土壤pH值、氣候濕潤度和喬木地上有機碳密度外，其他自變量在隨機森林模型中的偏相關趨勢與在逐步線性回歸中的作用方向一致，說明ΔSOCD隨土壤黏粒+粉粒含量(小于80%時)、喬木地下有機碳密度和造林年限(小于60年時)增加而增加，隨土壤初始有機碳密度增加而減少，而與土壤pH值的關系更接近二次曲線(圖3c)。Rowley等指出酸性土壤通過Fe3+和Al3+的吸附作用保護土壤有機碳，堿性土壤主要通過Ca2+保護土壤有機碳，只有中性土壤是土壤微生物利用土壤有機碳的“機會窗口”[31]，該結論與本研究曲線形狀表現的意義一致且偏相關曲線的拐點也與“機會窗口”的pH值范圍基本一致。ΔSOCD與氣候濕潤度的關系也表現為二次曲線，當溫度與降水協同更有利于土壤呼吸增加時，土壤有機碳密度也隨之降低，表現為ΔSOCD處于低值，從“三基點”的角度解釋，達到最大土壤呼吸的氣候濕潤度小于達到最大土壤有機碳輸入的氣候濕潤度，ΔSOCD在該點之前隨著氣候濕潤度增加而減少，之后隨著氣候濕潤度增加而增加，因此，ΔSOCD與氣候濕潤度的曲線趨勢與文獻[15]中氣候帶對土壤有機碳密度變化速率的影響趨勢相似。

喬木地上有機碳密度的作用在兩個模型中存在差異(圖3e)，喬木地上、地下有機碳密度在逐步線性回歸中的非標準化系數分別為-0.07和0.36(表1)。實際情況下喬木地上、地下生長同步，因此，從兩者系數的比例看，只有喬木地下與地上有機碳密度比例小于1∶5.14的樹種，土壤有機碳密度才會隨著喬木有機碳密度增加而減少，可近似認為根冠比小于0.19的樹種造林后土壤有機碳密度會下降。由于針葉樹和闊葉樹的生物量分配策略不同，闊葉樹的根系發達[32]，使得地下生物量更高的闊葉樹提高了根系對土壤有機碳的輸入，故樹種的根冠比越大，越有利于土壤有機碳積累。有研究顯示，根據文獻數據計算出華北落葉松的根冠比是0.26，而野外采樣獲得的華北落葉松的根冠比是0.17[33]，不同研究中由于樣點針葉樹種的根冠比存在差異，導致土壤有機碳密度變化的計算結果出現土壤有機碳損失[13,34]、累積或者保持不變[35]。

在隨機森林模型中，通常根據IncMSE指標(去除自變量后隨機森林預測誤差增加的百分數)明確自變量對因變量的影響程度，數值越大說明自變量越重要。基于此，對逐步線性回歸模型和隨機森林模型的自變量進行重要性排序(圖4)。由圖4a可知，隨機森林模型中，PSOCD自變量的重要性從大到小為土壤初始有機碳密度(76.59%)、喬木地上有機碳密度(73.73%)、喬木地下有機碳密度(42.97%)、氣候濕潤度(33.22%)、土壤pH值(28.39%)、土壤黏粒+粉粒含量(26.98%)、造林年限(26.12%)、植被類型(18.33%)和土地利用類型(17.01%)，而逐步線性回歸模型中自變量排序與之略有不同，但主導因素一致，均為土壤初始有機碳密度、喬木地上有機碳密度、喬木地下有機碳密度，而土壤pH值、土壤黏粒+粉粒含量和造林年限是次要因素。由圖4b可知，ΔSOCD的自變量在隨機森林模型中的重要性從大到小為土壤黏粒+粉粒含量(52.86%)、土壤初始有機碳密度(52.71%)、造林前土地利用類型(49.94%)、土壤pH值(41.89%)、氣候濕潤度(36.79%)、植被類型(28.24%)、喬木地下有機碳密度(21.88%)、喬木地上有機碳密度(21.75%)和造林年限(21.39%)。其中，喬木地上、地下有機碳密度和造林年限的回歸方程標準化系數排序靠前，但重要性排序靠后；而氣候濕潤度和土壤pH值的重要性排序靠前，但標準化系數排序靠后。兩種模型結果出現差異的原因可能是統計理論方法不同，逐步線性回歸模型是對全部數據進行方差分解，而隨機森林模型采用有放回抽樣形成不同子集，分析每個自變量對因變量預測結果的影響。但兩模型中土壤黏粒+粉粒含量和土壤初始有機碳密度均位居前列，而逐步線性回歸方程中造林前土地利用類型也達到顯著水平。綜合考慮可認為主導ΔSOCD的因素為土壤黏粒+粉粒含量、土壤初始有機碳密度和造林前土地利用類型，其他因素可歸為次要因素。

注：縱軸旁數字表示逐步線性回歸中標準化系數絕對值的排序。

2.2 退耕還林后土壤有機碳指標變化情景模擬

本文將10折交叉檢驗的均方根誤差(RMSE)作為自變量篩選依據，RMSE越小說明預測值與觀測值越接近。重復5次檢驗后得到土壤有機碳指標的RMSE與自變量數量的關系(圖5)，可以看出，PSOCD、ΔSOCD的RMSE在自變量數量大于4后，基本不變。結合圖4的分析結果，選擇圖4中IncMSE排序前四的自變量(PSOCD的模擬變量為喬木地上/地下有機碳密度、土壤初始有機碳密度和氣候濕潤度，ΔSOCD的模擬變量為土壤黏粒+粉粒含量、土壤pH值、土壤初始有機碳密度和造林前土地利用類型)作為隨機森林模型的驅動變量，進而模擬我國不同區域退耕還林后土壤有機碳時空變化(表2)。

圖5 基于隨機森林模型的0～20 cm土壤有機碳指標交叉檢驗Fig.5 Cross-validation of SOC0～20 estimation based on random forest model

表2 基于隨機森林模型估算我國退耕還林后0～20 cm PSOCD 和ΔSOCD變化Table 2 Estimation of PSOCD and ΔSOCD changes in 0～20 cm soil layer under the "Grain-for-Green" Program based on random forest model

由表2可知，我國各氣候區針葉樹的PSOCD區域均值從造林5年后的80.15%～92.87%下降到造林40年后的27.35%～61.09%，5～40年間降幅從寒溫帶(52.71%)到熱帶(61.61%)逐漸增大，從干燥氣候區(干旱和半干旱區，29.57%～35.37%)到濕潤氣候區(52.71%～61.61%)降幅也逐漸增大。闊葉樹PSOCD均值從造林5年后的61.21%～85.29%下降到造林40年后的25.67%～53.94%，期間降幅從寒溫帶(31.40%)到亞熱帶(41.92%)逐漸增大，從干燥氣候區(31.35%～33.73%)到熱帶和亞熱帶氣候區(34.21%～41.92%)降幅也逐漸增大。在內陸(干旱、半干旱和高寒)地區，闊葉樹造林40年后PSOCD均值小于針葉樹，而其他氣候區闊葉樹造林40年后PSOCD均值大于針葉樹(表2)，造林40年后暖濕地區的PSOCD均值小于冷干地區。退耕還林后ΔSOCD隨造林時間有不同程度增加，在內陸地區針葉樹土壤有機碳密度的增量和增幅均小于闊葉樹，而其他氣候區針葉樹土壤有機碳密度的增量和增幅均大于闊葉樹。造林后5～40年針葉樹土壤有機碳密度增幅最大的是中溫帶(5.43 Mg C hm-2)，增幅最小的是干旱區(0.23 Mg C hm-2)；闊葉樹土壤有機碳密度增幅最大的是高寒區(8.87 Mg C hm-2)，增幅最小的是半干旱區(1.40 Mg C hm-2)，高寒氣候區針葉樹和闊葉樹對土壤有機碳密度增量(造林40年后)和增幅(造林后5～40年間)的影響差異很大。

Deng等研究我國退耕還林土壤有機碳密度變化時發現，在針葉樹造林后31～40年內，0～20 cm土壤固碳速率大于闊葉樹[36]，人工針葉林0～20 cm土壤的平均固碳速率為0.03 Mg C hm-2a-1，人工闊葉林為0.17 Mg C hm-2a-1[36]。與之相比，本研究預測我國造林40年后針葉林和闊葉林0～20 cm土壤平均固碳速率分別為0.21 Mg C hm-2a-1和0.22 Mg C hm-2a-1(氣候區內樣本量加權平均)，模型預測的退耕還林土壤固碳作用較高，尤其是針葉林土壤固碳速率。但陳越的研究顯示，青海、東北和新疆地區針葉樹造林后0～20 cm土壤固碳速率為0.58 Mg C hm-2a-1，大于相同區域中闊葉樹造林后土壤平均固碳速率(0.30 Mg C hm-2a-1)[27]，與本研究模擬的造林后5～40年間針葉林ΔSOCD增幅大于闊葉林ΔSOCD增幅結果一致。因此考察針葉樹造林后的ΔSOCD時，需考慮更長時間尺度的土壤有機碳積累。

3 結論與展望

本文基于收集的全球造林后土壤有機碳指標變化文獻數據，通過逐步線性回歸模型和隨機森林模型分析影響土壤有機碳變化的因素，并利用統計分析和隨機森林模型模擬中國造林后土壤有機碳指標的變化。結論如下：1)土壤有機碳密度占生態系統有機碳密度的比例(PSOCD)主要受土壤初始有機碳密度和喬木地上、地下有機碳密度影響；土壤有機碳密度變化量(ΔSOCD)主要受土壤質地和土壤初始有機碳密度影響。土壤黏粒+粉粒含量越高，越有利于土壤固碳；但土壤初始有機碳密度越高的地區，造林后更易出現土壤有機碳密度下降，因此在農田或荒地上造林更有利于表層土壤有機碳密度增加。土壤有機碳輸入隨著喬木生長而增加，樹種根冠比越大，越有利于土壤有機碳密度增加。但對整個生態系統而言，比土壤有機碳密度增加速度更快的森林碳增加速度使土壤固碳對生態系統固碳的貢獻率隨時間下降，因此考察人工林生態系統碳匯效應時，需考慮長時間尺度的土壤有機碳積累。2)隨機森林模型模擬結果說明，我國造林40年后ΔSOCD自北向南隨年均溫增大而增大，從內陸向沿海隨年降水量增大而增大，但造林后5～40年間，寒溫帶地區ΔSOCD的增幅更大，說明暖濕地區造林后植被生長速度快，土壤固碳效果更好，而寒冷地區的固碳潛力更大；造林40年后PSOCD自北向南、自內陸向沿海下降且幅度逐漸增大。研究結果可為未來基于大量數據利用隨機森林模型預測區域或國家尺度造林后土壤有機碳指標變化提供支持，同時為在管理措施和造林技術方面保護土壤碳庫提供指導。