基于Stacking模型的土壤綜合肥力評(píng)價(jià)
——以富川瑤族自治縣植煙區(qū)為例

鄒天祥,梁志鵬,龔佳林,周 萌,沈文杰,2,3,張介棠,范東升,盧燕回

(1.中山大學(xué)地球科學(xué)與工程學(xué)院,廣東 珠海 519000;2.廣東省地質(zhì)過程與礦產(chǎn)資源探查重點(diǎn)實(shí)驗(yàn)室,廣東 珠海 519000;3.廣東省地球動(dòng)力作用與地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室,廣東 珠海 519000;4.廣東微碳檢測科技有限公司,廣東 清遠(yuǎn) 511500;5.中國煙草總公司廣西壯族自治區(qū)公司,南寧 530022)

【研究意義】土壤是植物生長的基礎(chǔ),其營養(yǎng)元素和理化性質(zhì)決定了作物對(duì)養(yǎng)分的吸收和轉(zhuǎn)化效率;氣候?yàn)檗r(nóng)業(yè)發(fā)展提供了光、熱、水等物質(zhì)條件,并在很大程度上決定當(dāng)?shù)刈魑锓N植制度;地形地貌通過影響環(huán)境中的溫濕度、養(yǎng)分含量等,對(duì)作物種植和生長產(chǎn)生制約。因此,準(zhǔn)確評(píng)價(jià)土壤肥力及其影響因素至關(guān)重要,以便針對(duì)不同地區(qū)制定適宜的作物利用方案,實(shí)現(xiàn)作物優(yōu)質(zhì)生長和土地資源可持續(xù)發(fā)展?！厩叭搜芯窟M(jìn)展】土壤肥力是農(nóng)作物種植研究的重要課題,但研究者多側(cè)重于土壤營養(yǎng)物質(zhì)和礦質(zhì)成分的儲(chǔ)量和有效性研究,對(duì)土壤綜合肥力的探索不足,相關(guān)標(biāo)準(zhǔn)尚未統(tǒng)一。傳統(tǒng)綜合肥力評(píng)價(jià)方法(包括主成分分析法[1]、指數(shù)法、最小數(shù)據(jù)集法、內(nèi)梅羅指數(shù)法[2]、模糊綜合評(píng)價(jià)法[3]等)存在過程繁瑣、評(píng)價(jià)指標(biāo)多樣、主觀影響大和標(biāo)準(zhǔn)難以統(tǒng)一的缺點(diǎn)。近年來,隨著信息技術(shù)的快速發(fā)展和大數(shù)據(jù)行業(yè)的興起,機(jī)器學(xué)習(xí)在多個(gè)領(lǐng)域得到廣泛應(yīng)用。例如賴紅松等[4]采用遺傳算法-模擬退火(GASA)優(yōu)化支持向量機(jī)(SVM)參數(shù)算法實(shí)現(xiàn)對(duì)農(nóng)田地力的等級(jí)評(píng)價(jià);李小剛等[5]運(yùn)用BP神經(jīng)網(wǎng)絡(luò)結(jié)合模糊綜合評(píng)價(jià)實(shí)現(xiàn)對(duì)耕地土壤肥力的等級(jí)劃分。盡管已有關(guān)于土地肥力評(píng)價(jià)的機(jī)器學(xué)習(xí)和大數(shù)據(jù)案例和嘗試,但針對(duì)不同地區(qū)、不同作物土壤肥力的劃分標(biāo)準(zhǔn)無法通用,需要針對(duì)具體地區(qū)作物進(jìn)行相應(yīng)的土壤肥力評(píng)價(jià)?！颈狙芯壳腥朦c(diǎn)】煙草是我國重要經(jīng)濟(jì)作物,可塑性強(qiáng)、環(huán)境適應(yīng)性廣。富川瑤族自治縣是賀州市煙草的重要產(chǎn)區(qū),煙草產(chǎn)業(yè)是當(dāng)?shù)貙?shí)施精準(zhǔn)扶貧的重要窗口。然而當(dāng)?shù)蒯槍?duì)植煙區(qū)土壤綜合肥力方面的研究相對(duì)薄弱,尚未進(jìn)行系統(tǒng)的土壤肥力評(píng)價(jià)調(diào)查?！緮M解決的關(guān)鍵問題】為優(yōu)化富川縣植煙區(qū)烤煙的種植布局,合理利用植煙區(qū)生態(tài)資源,本文以富川植煙區(qū)為例,采用Stacking堆疊算法構(gòu)建植煙區(qū)的土壤綜合肥力評(píng)價(jià)模型,實(shí)現(xiàn)對(duì)土壤綜合肥力評(píng)價(jià)的標(biāo)準(zhǔn)化和快捷化。

1 材料與方法

1.1 研究區(qū)概況

廣西省賀州市富川瑤族自治縣位于廣西東北部,總面積157 200 hm2,下轄9鎮(zhèn)3鄉(xiāng)?？h域地勢呈北高南低,四周是丘陵山脈和巖溶峰林,中部是寬闊的溶蝕平原[6]。研究區(qū)屬亞熱帶季風(fēng)氣候,陽光資源豐富,晝夜溫差明顯,年降雨量達(dá)1700 mm,多年平均氣溫19.1 ℃,歷年總積溫7008 ℃,全年≥10 ℃的活動(dòng)積溫為6099 ℃,平均日照時(shí)數(shù)為1573.6 h,區(qū)內(nèi)以輕壤土、中壤土為主要土壤類型,部分地區(qū)為砂土、黏土和重壤土,具備生產(chǎn)優(yōu)質(zhì)煙葉的生長環(huán)境和氣候條件。

1.2 評(píng)價(jià)因子選取

土壤肥力受眾多因素影響,可歸納為土壤條件、氣候環(huán)境和地形地貌3大類。僅土壤條件就包含土壤物理指標(biāo)、化學(xué)指標(biāo)、生物指標(biāo)和養(yǎng)分指標(biāo)等多種緯度, 涉及指標(biāo)繁多且部分指標(biāo)的測定和應(yīng)用存在困難[5]。本研究在同類研究[7-11]的基礎(chǔ)上,篩選出15種易于測定且重現(xiàn)性好的指標(biāo)作為影響因子,指標(biāo)的選取和簡要解釋見表1[12-13]。

表1 肥力評(píng)價(jià)指標(biāo)體系

圖1 采樣點(diǎn)空間分布

表2 土壤養(yǎng)分的測定指標(biāo)及測定方法

1.3 數(shù)據(jù)來源

1.3.1 土壤樣品采集和測試 根據(jù)植煙土壤類型、耕作模式等具體狀況,將采樣區(qū)域劃分為若干單元,盡量保證每個(gè)單元的土壤性狀均勻一致。具體采樣方法如下:對(duì)于連片的煙田,每3.33 hm2作為一個(gè)單元,采用“S”形布點(diǎn)采樣;不連片或狹小的煙田,則單獨(dú)作為一個(gè)采樣單元,采用“梅花”形布點(diǎn)采樣。土樣采集為耕作層0～20 cm表層土壤,采用多點(diǎn)混合采樣法,在每個(gè)采樣單元采集5～8個(gè)點(diǎn)混合而成混合樣。在本研究中,共采集1360個(gè)樣點(diǎn),其中廣西省河池市、百色市和賀州市鐘山縣地區(qū)總計(jì)1038個(gè)樣點(diǎn)用作訓(xùn)練和測試數(shù)據(jù),富川縣322個(gè)樣點(diǎn)用于研究區(qū)數(shù)據(jù),研究區(qū)采樣點(diǎn)空間分布見圖1。樣品經(jīng)自然風(fēng)干、去雜、研磨、過篩后,進(jìn)行養(yǎng)分指標(biāo)的測定分析(表2)。

1.3.2 氣候環(huán)境 氣候數(shù)據(jù)來源于中國地面氣候資料日值數(shù)據(jù)集(V3.0)。本研究選擇3—6月大田期數(shù)據(jù),統(tǒng)計(jì)廣西全域及周圍共計(jì)42個(gè)氣象站站點(diǎn)2010—2019年日值數(shù)據(jù),計(jì)算出各站點(diǎn)旺長期降雨量、大田期日均氣溫、大田期日均相對(duì)濕度、大田期活動(dòng)積溫和大田期日照時(shí)數(shù),采用合適的插值方法[14]獲取不同影響因子的氣候信息數(shù)據(jù)。

表3 各影響因子隸屬函數(shù)、閾值及其權(quán)重

1.3.3 地形地貌 DEM數(shù)據(jù)來源于BIGMAP(http://www.bigemap.com),空間分辨率17.5 m。使用ArcGIS表面分析工具獲取海拔和坡度數(shù)據(jù)。

1.4 模糊數(shù)學(xué)模型

訓(xùn)練集數(shù)據(jù)的建立使用模糊數(shù)學(xué)模型法[15],通過隸屬函數(shù)將無規(guī)則分布的定量或定性描述數(shù)據(jù)轉(zhuǎn)換為0～1的隸屬值,其隸屬函數(shù)公式包括拋物線型、S型和反S型[16],公式如下:

(1)

(2)

(3)

式中,x1、x2、x3、x4分別表示隸屬函數(shù)的下限、下優(yōu)、上優(yōu)和上限。參評(píng)指標(biāo)的因子權(quán)重采用AHP層次分析法[8]計(jì)算,且判斷矩陣均通過一致性檢驗(yàn)。具體影響因子的隸屬函數(shù)閾值[8, 17-19]和權(quán)重如表3所示。其中,pH、有機(jī)質(zhì)、全氮、容重、孔隙度、旺長期降雨量、大田期日均氣溫、大田期相對(duì)濕度、大田期日照時(shí)數(shù)、海拔為拋物線形隸屬函數(shù),有效磷、速效鉀、大田期活動(dòng)積溫為S形隸屬函數(shù),坡度為反S形隸屬函數(shù)。土壤質(zhì)地采用直接賦值法,砂壤土、輕壤土、中壤土、黏土、砂土和重壤土分別賦值1.0、1.0、0.9、0.8、0.8和0.8。

得到對(duì)應(yīng)樣點(diǎn)的綜合肥力指數(shù)IFI(Integrated fertility index),計(jì)算公式如下:

(4)

式中,n為影響因子總數(shù),Pj為第j個(gè)影響因子權(quán)重,Xij表示第i個(gè)樣本第j個(gè)影響因子的隸屬度。IFI的取值范圍為0.1～1.0,IFI值越高表明該樣點(diǎn)的綜合肥力評(píng)價(jià)越好。

1.5 Stacking模型

Stacking 是一種混合學(xué)習(xí)方法[20-21],其核心思想在于通過整合多個(gè)差異性較大的模型,從而實(shí)現(xiàn)相較單一模型更佳的預(yù)測表現(xiàn)。通常將各獨(dú)立模型稱為初級(jí)學(xué)習(xí)器,而負(fù)責(zé)整合各初級(jí)學(xué)習(xí)器的模型被稱作次級(jí)學(xué)習(xí)器或者元學(xué)習(xí)器[22]。Stacking模型首先針對(duì)初始數(shù)據(jù)集進(jìn)行初級(jí)學(xué)習(xí)器的訓(xùn)練,并生成用于元學(xué)習(xí)器訓(xùn)練的次級(jí)訓(xùn)練集數(shù)據(jù),在此基礎(chǔ)上得到最終預(yù)測結(jié)果(圖2)。

由于數(shù)據(jù)存在噪聲,各模型在不同特征方面的表現(xiàn)存在差異。在理想情況下,采用Stacking訓(xùn)練方法可以提煉各模型在特征提取優(yōu)秀的部分,同時(shí)摒棄其在某些方面不佳的表現(xiàn)。然而,直接使用初級(jí)學(xué)習(xí)器的訓(xùn)練集來產(chǎn)生次級(jí)訓(xùn)練集存在過擬合風(fēng)險(xiǎn),需要采用交叉驗(yàn)證方法來形成次級(jí)訓(xùn)練集[23]。如圖3所示,以5折交叉驗(yàn)證為例,在訓(xùn)練初級(jí)學(xué)習(xí)器時(shí),將訓(xùn)練數(shù)據(jù)隨機(jī)劃分5份,每次用其中4份訓(xùn)練并預(yù)測剩下的1份數(shù)據(jù),重復(fù)5次即可獲得完整訓(xùn)練集的預(yù)測結(jié)果。同時(shí),每個(gè)初級(jí)學(xué)習(xí)器都需要對(duì)測試數(shù)據(jù)集進(jìn)行預(yù)測,從而得到5個(gè)測試集預(yù)測結(jié)果的平均值,作為次級(jí)學(xué)習(xí)器的一個(gè)特征緯度。

1.6 模型搭建

為了增強(qiáng)Stacking模型的泛化性能,需要考慮以下兩個(gè)方面:其一,選擇具有不同特征的學(xué)習(xí)器作為初級(jí)學(xué)習(xí)器以提高模型的多樣性,減少過擬合現(xiàn)象;其二,對(duì)每個(gè)選定的學(xué)習(xí)器進(jìn)行參數(shù)優(yōu)化,以提高模型的精度和泛化性能。通過模型的優(yōu)化組合和調(diào)參,最終選擇隨機(jī)森林[24](rf)、極端樹[25](et)、梯度提升樹(GBDT)、XGBoost[26]和Catboost作為初級(jí)學(xué)習(xí)器,選擇線性回歸(lr)作為元學(xué)習(xí)器,所有模型基于python搭建,建立過程如下:① 通過模糊數(shù)學(xué)模型獲取初始訓(xùn)練數(shù)據(jù)集,采取8∶2劃分,80%數(shù)據(jù)用作訓(xùn)練集,20%數(shù)據(jù)作為測試集。② 選擇隨機(jī)森林(rf)、極端樹(et)、梯度提升樹(GBDT)、XGBoost和Catboost作為初級(jí)學(xué)習(xí)器,使用網(wǎng)格交叉搜索確定模型參數(shù)[27]。③ 構(gòu)建一個(gè)1×830、1×208和5×208的矩陣;使用5折交叉訓(xùn)練,每一折訓(xùn)練166個(gè)訓(xùn)練集樣本,同時(shí)預(yù)測每一折的208個(gè)測試樣本。每一折預(yù)測產(chǎn)生166個(gè)預(yù)測值,堆疊成1×830的次級(jí)訓(xùn)練集數(shù)據(jù)矩陣;每一折預(yù)測形成1×208測試集預(yù)測值,填充后形成5×208的預(yù)測值矩陣,取平均得到次級(jí)數(shù)據(jù)集的預(yù)測數(shù)據(jù)。④ 基于線性回歸的元學(xué)習(xí)器對(duì)次級(jí)數(shù)據(jù)集進(jìn)行訓(xùn)練,從而構(gòu)建融合5種初級(jí)學(xué)習(xí)器的Stacking模型。

圖2 Stacking模型簡要流程

圖3 k折交叉訓(xùn)練流程

表4 Stacking01模型中各學(xué)習(xí)器性能指標(biāo)

2 結(jié)果與分析

2.1 模型性能度量

為有效度量模型性能,選取平均絕對(duì)誤差(MAE)、均方根誤差(RMSE)和決定系數(shù)(R2)為性能度量指標(biāo)。

(5)

(6)

(7)

Stacking01模型的泛化性能顯著優(yōu)于所有初級(jí)學(xué)習(xí)器,R2系數(shù)達(dá)0.9858,相較于初級(jí)學(xué)習(xí)器Rf01、Et01、GBDT01、Xgboost01和Catboost01模型有0.2066、0.2168、0.0432、0.0322和0.0063的提升,平均提升幅度為0.1010,表現(xiàn)出顯著優(yōu)勢。調(diào)用特征重要性指令表明,各學(xué)習(xí)器特征重要度前6的特征均為有機(jī)質(zhì)、全氮、穩(wěn)定積溫、速效鉀、pH和有效磷。為簡化分析,將模型特征精簡為重要度前6的關(guān)鍵因子并構(gòu)建Stacking02模型,在表5中展示各初級(jí)學(xué)習(xí)器和Stacking02模型的超參數(shù)集及性能度量。

Stacking02模型的MAE、RMSE和R2分別為0.0066、0.0084和0.9529。圖4-a展示了模型的線性擬合結(jié)果,擬合斜率為0.9438,呈現(xiàn)出優(yōu)秀的擬合結(jié)果。圖4-b和4-c分別描述了常規(guī)殘差分布及其區(qū)間占比,所有測試樣點(diǎn)的殘差均分布在-0.025～0.025,且集中于-0.01～0.01,占總數(shù)的77.88%。

2.2 富川縣植煙區(qū)肥力綜合評(píng)價(jià)預(yù)測

使用訓(xùn)練完成的Stacking02模型對(duì)研究區(qū)肥力進(jìn)行綜合評(píng)價(jià)預(yù)測,其預(yù)測結(jié)果見表6,研究區(qū)IFI的最小值、平均值、25%位點(diǎn)、50%位點(diǎn)和75%位點(diǎn)均高于廣西其它植煙區(qū)。圖5展示了預(yù)測結(jié)果的空間分布,根據(jù)四分位點(diǎn)把IFI劃分為4個(gè)等級(jí),分別為Ⅳ級(jí)(0.720～0.821)、Ⅲ級(jí)(0.821～0.851)、Ⅱ級(jí)(0.851～0.875)、Ⅰ級(jí)(0.875～0.930)。研究區(qū)北部綜合肥力水平優(yōu)于南部,位于北部的石家鄉(xiāng)、富陽鎮(zhèn)、葛坡鎮(zhèn)和朝東鎮(zhèn)Ⅰ級(jí)和Ⅱ級(jí)產(chǎn)區(qū)占比分別達(dá)66.67%、66.67%、61.29%和59.26%,而位于研究區(qū)南部的新華鄉(xiāng)、蓮山鎮(zhèn)、柳家鄉(xiāng)的Ⅰ級(jí)和Ⅱ級(jí)產(chǎn)區(qū)占比較低,分別為15.79%、11.11%和0。由表7可知,相較于IFI后25%的樣點(diǎn),IFI前25%的樣點(diǎn)平均土壤pH低0.36,平均全氮含量高0.30 g/kg,平均速效鉀含量高139.01 g/kg,低IFI區(qū)域存在土壤pH偏高而全氮和速效鉀含量偏低的問題[28-29]。

表5 Satcking02模型中各學(xué)習(xí)器性能度量和超參數(shù)集

圖4 Stacking02模型測試集預(yù)測結(jié)果可視化

圖5 富川縣植煙區(qū)肥力綜合評(píng)價(jià)空間分布

表6 富川植煙區(qū)預(yù)測結(jié)果與廣西其它地區(qū)對(duì)比

3 討 論

Stacking算法是通過組合多個(gè)模型以生成新模型的集成方法,其第一層學(xué)習(xí)器的目標(biāo)是對(duì)原始數(shù)據(jù)進(jìn)行特征抽取,并利用異質(zhì)性來表示不同特征。該模型的有效性取決于第一層學(xué)習(xí)器特征抽取的質(zhì)量,而第二層是基于對(duì)第一層數(shù)據(jù)的學(xué)習(xí),使得模型存在過擬合風(fēng)險(xiǎn)。因此,在訓(xùn)練初級(jí)學(xué)習(xí)器時(shí)需要使用交叉驗(yàn)證,并選擇簡單的模型用于結(jié)合初級(jí)學(xué)習(xí)器。

基于6種關(guān)鍵因子構(gòu)建的Stacking02模型對(duì)研究區(qū)的預(yù)測結(jié)果表明,研究區(qū)綜合肥力水平較好,但低IFI區(qū)域存在土壤pH偏高,且全氮和速效鉀含量偏低。研究表明土壤pH在5.5～6.5的烤煙質(zhì)量較好[30],而氮肥和有機(jī)肥的施用可致土壤pH下降[31-33],因此需要因地制宜地補(bǔ)充氮肥、有機(jī)肥和鉀肥。此外,研究區(qū)西南紅壤分布廣泛,其抗蝕能力差且長期連作使土壤全氮、全鉀、有效磷和有機(jī)質(zhì)降低,進(jìn)而導(dǎo)致土壤理化性質(zhì)惡化和煙葉品質(zhì)降低[34],建議在擴(kuò)大南部煙區(qū)種植面積的同時(shí),采用稻煙輪作、增加復(fù)種指數(shù)等管理措施。

4 結(jié) 論

相較于各初級(jí)學(xué)習(xí)器,Stacking模型在MAE、RMSE和R2性能度量上均有提升。通過對(duì)關(guān)鍵特征因子的優(yōu)選,Stacking02模型在強(qiáng)大的泛化性能和計(jì)算效率上取得了平衡,其R2決定系數(shù)達(dá)0.9529,表明該模型可擴(kuò)展至其它地區(qū)及作物的土壤綜合肥力評(píng)價(jià),進(jìn)而制定針對(duì)性的作物利用方案;采用Stacking02模型對(duì)研究區(qū)進(jìn)行肥力綜合評(píng)價(jià),結(jié)果表明研究區(qū)北部綜合肥力水平優(yōu)于南部。低IFI區(qū)域普遍存在pH偏高而全氮速效鉀含量偏低的問題,建議合理施用肥料、采用稻煙輪作等管理措施,在增加氮素鉀素的同時(shí)降低土壤pH。