長江中下游地區耕地土壤有機碳密度變化率驅動因素

張海星 ，趙智強 ，王瑞鋒 ，劉 靜 ，楚天舒 ,4※

（1. 中國農業大學 資源與環境學院,北京 100193；2. 中國農業大學 工學院,北京 100083；3. 上海市建設用地和土地整理事務中心,上海 200003；4. 自然資源部大都市區國土空間生態修復工程技術創新中心,上海 200003）

0 引言

全球氣候變化與極端天氣頻發對國家安全和經濟社會可持續發展帶來了嚴重威脅，已成為全球性挑戰。2015 年，聯合國可持續發展目標和《巴黎氣候協定》共同倡議，各國努力建設可持續低碳經濟。2019 年中國溫室氣體排放的二氧化碳當量為140 億t，占到全球排放量的27%[1]。從生產部門分析，農業生產溫室氣體排放量占人為溫室氣體排放總量的23%，但耕地也是巨大的碳庫，其固碳量約占陸地生態系統總固碳量的21%[2]。隨著農業綠色發展的大力推進，探索提高耕地土壤有機碳含量的途徑，顯得至關重要。

近年來，國內外圍繞農田土壤有機碳展開了大量研究。HAMID 等[3]研究表明，在長期耕作條件下，種植年限的增加會降低農田總有機碳的可用性，其中2018 年的有機碳含量較2014 年降低了45.32%。而ZHOU 等[4]在中國黃土高原區研究發現，相較于小麥—小麥輪作，1984—2010 年苜蓿地表層土壤有機碳含量增加了24.4%。并且在退耕還草中，長期適度放牧也可以促進有機碳的積累[5]。董林林等[6]通過10 a 定位試驗發現，稻麥秸稈還田可增加太湖地區水田土壤有機碳含量0.18～0.46 g/(kg·a)并增強其穩定性。從中國科學院桃源農業生態實驗站25 a雙季稻定位試驗結果可知，氮磷鉀平衡施肥處理組的土壤有機碳含量較不施肥處理組高出7.7%[7]。相對于不施用氮肥，長期施用氮肥提高了農田土壤有機碳含量5.2%～12.0%，但其時間響應取決于土壤類型，即在有機碳含量較高的土壤類型中具有放大效應[8]。在西班牙穆爾西亞省果園中，長期采用免耕與綠肥的組合方式可抵消耕作對土壤團聚體破裂的影響，并促進新團聚體的形成，與僅免耕相比，免耕與綠肥組合可提高有機碳含量14.0%[9]。在東北黑土區，相對于傳統翻耕模式，保護性耕作能顯著增加土壤表層有機碳含量25.8%[10]，并能增加大團聚體比例及團聚體穩定性[11]。總的來說，多種生產措施可促進耕地土壤固碳增匯，而如何系統性評估土壤固碳的相關措施，找出主要影響因素，還有待進一步研究。

長江中下游地區作為中國傳統農業產區，如何實現耕地固碳增匯，是該地區農業綠色發展的關鍵問題。當前缺乏系統性探究該地區耕地土壤有機碳密度變化特征與驅動因素的分析。因此，本研究擬以長江中下游地區為研究區域，收集與整理耕地土壤有機碳密度長期觀測數據，量化不同耕地土壤有機碳密度變化率的差異特征，應用隨機森林模型探究土壤有機碳密度變化率的主要驅動因素，解析耕地土壤有機碳密度變化率對環境和管理因素的非線性響應關系，以期為耕地固碳增匯提供參考。

1 材料與方法

1.1 數據來源與篩選方法

本研究選取“有機碳”“有機質”“長期定位”“耕地”等關鍵詞，分別在中國知網（CNKI）和Web of Science文獻數據庫進行檢索，收集與整理了公開發表的有關長江中下游地區耕地土壤有機碳密度變化的研究論文。進而，按照如下標準進行文章篩選：1）研究區域為長江中下游地區，從省級行政區劃上涵蓋：上海、江蘇、浙江、安徽、江西、湖北、湖南；2）試驗方式為大田試驗，試驗年限不少于2 a，試驗地點、時間、基礎條件清楚，優選長期定位試驗；3）試驗土樣為表土（0～20 cm），土壤有機碳密度數據準確易獲取。此外，若文章中部分土壤基本信息缺失，則在世界土壤數據庫（Harmonized World Soil Database,version 1.2）中按經緯度查詢以獲取[12]。若文章中氣象基本信息缺失，如年降水量和年平均氣溫等，則在東安格利亞大學氣候研究中心氣象數據庫（Climate Research Unit）查詢以獲取[13]。

本研究按照上述方法篩選出長期試驗（表1），提取試驗年限、地點、耕地基本信息（耕地利用類型、試驗年限、土壤pH 值、黏粒含量、土壤有機碳含量）、作物基本信息（輪作模式、化肥氮量、有機肥氮量、秸稈氮量）、氣象基本信息（年平均降水量和年平均氣溫），建立長江中下游地區耕地表土有機碳密度變化率數據庫。其中，化肥氮量、有機肥氮量、秸稈氮量是由周年施用量乘以氮含量計算獲得。

表1 試驗基本信息Table 1 Basic information of experiment

1.2 耕地土壤有機碳密度變化率計算方法

若文獻未直接給出土壤有機碳密度變化率數據，則依據《2006 年IPCC 國家溫室氣體清單指南2019 修訂版》中方法計算獲取，具體如下：

式中ΔCSO為耕地土壤有機碳密度變化率，kg/(hm2·a)；CSO,t為第t年耕地土壤有機碳密度，kg/(hm2·a)；CSO,0為初始年耕地土壤有機碳密度，kg/(hm2·a)；t為試驗年限，a；c為耕地土壤有機碳含量，g/kg；b為耕地土壤容重，g/cm3；h為耕層厚度，cm；100 為單位轉化系數。

1.3 隨機森林算法

隨機森林算法作為典型機器學習算法之一，具有良好的處理非線性問題的能力，廣泛應用于農業預測與建模。在本研究中，將耕地基本信息、作物基本信息、氣象基本信息作為特征變量，而將表土有機碳密度變化率作為目標變量，采用隨機抽樣的方法確定訓練集與測試集。通過自助法采樣獲得構建N棵樹所需的n個子集，每次未抽中的數據稱為袋外數據，用于內部誤差估計和特征變量的重要性計算。對測試集的目標變量真值與預測值計算平均絕對誤差、均方誤差和均方根誤差，并進行誤差估計。選取對特征變量j隨機排列后預測值平均精確率減少值（mean decrease accuracy,MDA）表征特征變量的重要性并對其排序，篩選出影響耕地土壤有機碳密度變化率較大的特征變量作為驅動因素，基礎原理詳見式（2）。最終構建耕地土壤有機碳密度變化率的預測模型。算法采用Python 語言和scikit-learn 機器學習工具實現。

式中ij表示特征變量j的MDA 值；s表示原始數據集使用隨機森林回歸的模型分數；K表示對數據集中特征變量j數據進行重新隨機排列的總次數，sk,j表示在第k次對特征變量j隨機排列，原始數據集其他列不變的情況下，使用隨機森林回歸的模型分數；R2表示預測的決定系數；ytrue表示目標變量的實際值；ypred表示隨機森林回歸模型基于各特征變量得出的預測值；ytrue.mean表示目標變量實際值的均值。

1.4 施肥策略優化

2022 年11 月農業農村部發布《到2025 年化肥減量化行動方案》，將減少化肥用量、提高有機肥用量作為目標任務，促進農業綠色轉型。江蘇省作為國家農業綠色發展先行區，對于長江中下游地區農業綠色低碳發展具有引領示范作用。因此，本文以江蘇省為優化對象，選取典型水田和旱地輪作模式，將隨機森林模型與線性優化法相結合，以土壤有機碳密度變化率最大為目標函數，以氮肥總量和有機無機肥配施比例為約束條件，仿真分析出適宜的有機肥氮量和化肥氮量，以期協同實現耕地土壤固碳與化肥減量化，具體模型如下：

式中ΔCSO,pred(x,y)為隨機森林模型預測不同情況下耕地土壤有機碳密度變化率，kg/(hm2·a)；x為化肥氮量，kg/hm2；y為有機肥氮量，kg/hm2；a為氮肥總量，kg/hm2，參考農業農村部發布的作物科學施肥指導意見與當地推薦施肥量確定；r為有機肥氮量占氮肥總量比例，%；rmin為有機肥氮量占氮肥總量比例的最小值，%；rmax為有機肥氮量占氮肥總量比例的最大值，%。根據大田作物常規施肥策略與作物養分需求規律，若有機肥施用量過高，大田作物容易出現“貪青”現象，作物生殖生長期延后，直接影響作物產量。并且，兼顧長江中下游地區農業面源污染防控需要，參考《NY/T 4 163—2022 稻田氮磷流失綜合防控技術指南》，本研究將有機肥氮量占氮肥總量比例的最大值設置為30%。

根據策略優化結果，結合《NY/T 3 877—2021 畜禽糞便土地承載力測算方法》，定量評估江蘇省畜禽養殖是否能提供充足的有機糞肥資源支持優化后的施肥策略。

2 結果與分析

2.1 耕地土壤有機碳密度變化率特征分析

通過對59 處290 組長期定位試驗數據進行統計與分析可知，水田和旱地土壤有機碳密度變化率范圍分 別 為-1 548.15～3 577.10 kg/(hm2·a)和-261.89～3 245.01 kg/(hm2·a)。從土地利用類型來看，水田與旱地的土壤有機碳密度變化率差異不顯著（P=0.85）。從輪作模式來看（表2），針對水田，水稻—小麥、水稻—油菜的土壤有機碳密度變化率較高，但數據離散程度大；單作水稻的土壤有機碳密度變化率相對較低。針對旱地，蔬菜—蔬菜和油菜—棉花的土壤有機碳密度變化率較高，兩者間無顯著性差異；而玉米—玉米的土壤有機碳密度變化率最低。

表2 耕地土壤有機碳密度變化率Table 2 Cultivated land soil organic carbon density change rate (kg·hm-2·a-1)

2.2 隨機森林模型的構建

本研究以長期定位試驗數據為基礎，運用隨機森林算法，調節模型生成的決策樹數目超參數N=10 000，其余取默認值，獲得水田和旱地的土壤有機碳密度變化率模型（圖1）。其決定系數（R2）分別為0.63 和0.83，預測值與真實值高度相關[68]，表明隨機森林算法在耕地土壤有機碳密度變化率模型的構建上取得了較好的效果。

圖1 耕地土壤有機碳密度變化率（ΔCSO）的真實值與預測值散點圖Fig.1 Scatter diagrams of the true and predictive value of cultivated land soil organic carbon density change rate(ΔCSO)

2.3 驅動因素分析

2.3.1 驅動因素重要性排序

結合耕地土壤有機碳密度變化率的驅動因素重要性排序與系統聚類結果（表3）可知，對于水田，有機肥氮量、土壤pH 值、化肥氮量、秸稈氮量對有機碳密度變化率的影響較大，而試驗年限、黏粒含量、年平均降水量、土壤有機碳含量、年平均氣溫、水旱輪作、復種指數不是主要影響因素。對于旱地，有機肥氮量對有機碳密度變化率的影響最大，化肥氮量次之，而土壤有機碳含量、年平均降水量、試驗年限、土壤pH 值、年平均氣溫、黏粒含量、復種指數、秸稈氮量的影響較小。總的來說，有機肥氮量和化肥氮量是長江中下游地區耕地土壤有機碳密度變化率的主要驅動因素。

表3 耕地土壤有機碳密度變化率的驅動因素重要性排序Table 3 Sorting of the driving factors of cultivated land soil organic carbon density change rate

2.3.2 驅動因素的邊際依賴性分析

就水田而言，隨著有機肥氮量的增加，土壤有機碳密度變化率呈現“波動增長—下降”趨勢，最大值約為987 kg/(hm2·a)（圖2a）。隨著土壤pH 值從弱酸性至中性再至弱堿性，土壤有機碳密度變化率呈現出階梯式增長趨勢；當土壤pH 值為7.8 時，土壤有機碳密度變化率最高可達941 kg/(hm2·a)（圖2b）。當不施用化學氮肥時，土壤有機碳密度變化率極低；隨著化肥氮量增加，土壤有機碳密度變化率快速增加，而后在616～850 kg/(hm2·a)范圍間呈波動變化（圖2c）。另外，隨著秸稈氮量的增加，土壤有機碳密度變化率呈現波動增長趨勢，最大值約為904 kg/(hm2·a)（圖2d）。

圖2 水田驅動因素的部分依賴圖Fig.2 Partial dependence plots for the paddy field driving factors

科學調控土壤pH 值和秸稈還田量有利于耕地土壤固碳（圖2b 和圖2d）。在土壤pH 方面，雖然弱堿性土壤有利于固碳，但土壤pH 值是影響農作物產量的重要因素之一，土壤堿化可使農作物減產50%以上[69]。因此，本研究建議長江中下游地區合理調控水田土壤pH 值至中性，可兼顧實現耕地固碳與作物高產目標。在秸稈還田量方面，秸稈雖是重要碳源，但隨著秸稈還田量的增加，稻田甲烷排放量卻呈現線性增加趨勢，其增加的溫室效應會抵消土壤固碳的減排效益[70]。因此，未來研究還需要進行大量的田間試驗與模型仿真分析，探究適宜的秸稈還田量，協同實現土壤固碳與溫室氣體減排，以促進農田生態系統綠色低碳發展。

對于旱地而言，當不施用有機氮肥時，土壤有機碳密度變化率極低；隨著有機肥氮量增加，土壤有機碳密度變化率呈現出“增長—平穩”趨勢，最高值約為1 425 kg/(hm2·a)（圖3a）。對于一年一季的輪作模式來說，隨著化肥氮量的提升，土壤有機碳密度變化率也隨之增長，最大值在659 kg/(hm2·a)左右；而對于一年兩季的輪作模式來說，化肥氮量呈現“增長—波動下降”趨勢，當化肥氮量超過300 kg/hm2后，土壤有機碳密度變化率迅速提高，最高可達824 kg/(hm2·a)（圖3b）。

圖3 旱地驅動因素的部分依賴圖Fig.3 Partial dependence plots for the dry land driving factors

從上述分析可知，有機肥氮量和化肥氮量是長江中下游地區耕地土壤有機碳密度變化率的重要驅動因素。因此，進一步對有機肥氮量和化肥氮量進行雙驅動因素邊際依賴性分析，評估兩者是否協同促進耕地土壤有機碳密度變化率。從雙驅動因素部分依賴圖結果（圖4）可知，當有機肥氮量或化肥氮量過低時，耕地土壤有機碳密度變化率很低，不利于耕地土壤固碳；當有機肥氮量或化肥氮量過高時，耕地土壤有機碳密度變化率并非最大值，這說明有機肥或化肥過量施用，也不是土壤固碳的最佳途徑；相較于單驅動因素分析結果（圖2 和圖3），有機肥氮量和化肥氮量的雙驅動因素分析中，土壤有機碳密度變化率的最大值更大；這說明有機肥和化肥可協同提升耕地土壤有機碳密度。在水田中，有機肥氮量和化肥氮量分別為346～397 kg/hm2和422～533 kg/hm2，土壤有機碳密度變化率達到最大值，且較為穩定，約為1 156 kg/(hm2·a)。而在旱地中，有機肥氮量和化肥氮量分別在219～284 kg/hm2和338～394 kg/hm2范圍內，土壤有機碳密度變化率穩定在最大值，約為1 654 kg/(hm2·a)。由此可見，有機無機肥配施將有效促進耕地土壤固碳。

圖4 有機肥與化肥氮量雙驅動因素的部分依賴圖Fig.4 Partial dependence plots of driving factors for joint effects of organic fertilizer nitrogen and chemical fertilizer nitrogen

2.4 施肥策略優化

在氮肥總量控制、有機無機肥配施比例約束情況下，采用隨機森林模型與線性規劃仿真分析可知：隨著有機肥氮量占比增加，江蘇省水田土壤有機碳密度變化率呈現“增加—波動”趨勢；當有機肥氮量占氮肥總量比例超過10%，水田土壤有機碳密度變化率在1 143～1 335 kg/(hm2·a)間波動（圖5a）。旱地土壤有機碳密度變化率隨有機肥占氮肥總量比例增加而增加，近似冪函數關系（R2=0.89）；當有機肥氮量占比超過10%，水田土壤有機碳密度變化率增長趨緩，最高可達1 039 kg/(hm2·a)（圖5b）。總的來說，有機肥氮量占氮肥總量比例控制在10%～30%，可穩定地促進耕地土壤固碳，同時助力化肥減量化行動。

圖5 耕地土壤有機碳密度變化率的仿真分析結果Fig.5 Simulation results of cultivated land soil organic carbon density change rate

與此同時，當有機肥替代比例為30%時，2020 年江蘇省耕地需要有機肥氮量為38.31 萬t。核算可知，2020年江蘇省畜禽養殖可提供有機肥氮量為42.20 萬t，糞肥資源充足。

綜合上述研究成果可知，有機無機肥配施是提升長江中下游地區耕地土壤有機碳密度的重要措施，并且當地擁有充足的糞肥資源。因此，建議長江中下游地區大力推廣有機無機肥配施，并且有機肥氮量占氮肥總量比例控制在10%～30%，以助力化肥減量化行動，促進農田土壤固碳增匯，實現區域種養循環。

3 結論

本研究應用隨機森林模型與線性規劃探究長江中下游地區耕地土壤有機碳密度變化率的驅動因素與提升途徑，得到以下結論：

1）本研究整理與分析59 處長期定位試驗數據發現，長江中下游地區水田和旱地土壤有機碳密度變化率范圍分別為-1 548.15～3 577.10 kg/(hm2·a)和-261.89～3 245.01 kg/(hm2·a)。從土地利用類型來看，水田與旱地的土壤有機碳密度變化率差異不顯著（P=0.85）。

2）對于水田而言，有機肥氮量、土壤pH 值、化肥氮量、秸稈氮量對土壤有機碳密度變化率的影響較大，平均精確率減少值（mean decrease accuracy,MDA）分別為0.49、0.32、0.23、0.15，而試驗年限、黏粒含量、年平均降水量、土壤有機碳含量、年平均氣溫、水旱輪作不是主要影響因素。對于旱地而言，有機肥氮量對土壤有機碳密度變化率的影響最大、化肥氮量次之，MDA 值分別為0.98、0.10，土壤有機碳含量、年平均降水量、試驗年限、土壤pH 值、年平均氣溫、黏粒含量、復種指數、秸稈氮量的影響較小。有機肥氮量和化肥氮量是長江中下游地區耕地土壤有機碳密度變化率的主要驅動因素。

3）有機無機肥配施、調控土壤pH 值至中性、秸稈還田有利于水田土壤有機碳密度變化率增加；有機無機肥配施顯著提高旱地土壤有機碳密度變化率；耕地有機肥氮量占氮肥總量比例10%～30%為宜。