長江中下游地區(qū)耕地土壤有機(jī)碳密度變化率驅(qū)動因素

張海星 ，趙智強(qiáng) ，王瑞鋒 ，劉 靜 ，楚天舒 ,4※

（1. 中國農(nóng)業(yè)大學(xué) 資源與環(huán)境學(xué)院,北京 100193；2. 中國農(nóng)業(yè)大學(xué) 工學(xué)院,北京 100083；3. 上海市建設(shè)用地和土地整理事務(wù)中心,上海 200003；4. 自然資源部大都市區(qū)國土空間生態(tài)修復(fù)工程技術(shù)創(chuàng)新中心,上海 200003）

0 引言

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

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

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

1 材料與方法

1.1 數(shù)據(jù)來源與篩選方法

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

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

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

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

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

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

1.3 隨機(jī)森林算法

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

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

1.4 施肥策略優(yōu)化

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

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

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

2 結(jié)果與分析

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

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

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

2.2 隨機(jī)森林模型的構(gòu)建

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

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

2.3 驅(qū)動因素分析

2.3.1 驅(qū)動因素重要性排序

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

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

2.3.2 驅(qū)動因素的邊際依賴性分析

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

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

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

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

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

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

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

2.4 施肥策略優(yōu)化

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

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

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

綜合上述研究成果可知，有機(jī)無機(jī)肥配施是提升長江中下游地區(qū)耕地土壤有機(jī)碳密度的重要措施，并且當(dāng)?shù)負(fù)碛谐渥愕募S肥資源。因此，建議長江中下游地區(qū)大力推廣有機(jī)無機(jī)肥配施，并且有機(jī)肥氮量占氮肥總量比例控制在10%～30%，以助力化肥減量化行動，促進(jìn)農(nóng)田土壤固碳增匯，實現(xiàn)區(qū)域種養(yǎng)循環(huán)。

3 結(jié)論

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

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

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

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