基于GA-RF模型土壤堅實度對黑土區(qū)大豆產量的影響

周修理，張萍萍，秦 娜，霍東旭，喬金友*

（1.東北農業(yè)大學電氣與信息學院，哈爾濱 150030；2.東北農業(yè)大學工程學院，哈爾濱 150030；3.黑龍江省主要農作物生產機械化材料化技術創(chuàng)新中心，哈爾濱 150030）

黑土性狀好、肥力高，是適宜大豆生長的優(yōu)質土壤。黑龍江省位于世界三大黑土區(qū)之一的東北平原，大豆種植面積與產量均居全國首位，在保障國家糧食安全中具有不可替代地位[1]。近年來，隨著我國農業(yè)現(xiàn)代化水平提高，黑龍江地區(qū)農業(yè)機械化程度快速發(fā)展，農業(yè)機械作業(yè)對黑土壓實風險也隨之上升。過度壓實導致土壤結構破壞、水滲透能力減弱、作物根系生長受阻，影響糧食綜合生產能力提升和農業(yè)可持續(xù)發(fā)展[2]。

農業(yè)機械作業(yè)造成土壤壓實已成為制約世界各國農業(yè)發(fā)展的重要問題，國內外學者就機械壓實對土壤物理特性、作物生長發(fā)育影響等展開大量研究。在機械壓實對土壤物理特性影響上，Kumar等通過研究不同耕作系統(tǒng)壓實狀況發(fā)現(xiàn)，機械壓實破壞土壤結構，導致土壤容重與土壤堅實度增加，耕作可在一定程度上緩解壓實影響，疏松土壤[3]；Sivarajan等研究輪式機械作業(yè)對土壤壓實影響得出，機械作業(yè)對測區(qū)土壤壓實作用顯著，壓實后表層土壤易于恢復而深層土壤恢復較慢[4]；喬金友等在使用JD-904型拖拉機進行壓實試驗基礎上，結合多因素方差分析，發(fā)現(xiàn)土壤深度與壓實次數(shù)為影響土壤堅實度最主要因素，隨壓實次數(shù)增加，土壤受壓實作用范圍及程度均增加[5]。關于機械壓實對作物生長發(fā)育影響，肖芬芳等通過不同壓強下水田土壤壓實試驗，結合單因素方差分析及線性回歸擬合得出，機械壓實對水稻根系發(fā)育及產量造成顯著影響，土壤受壓實強度與甘蔗成根數(shù)量、主根長度、單穴產量之間為負相關[6]；Alvarez等研究表明，不同作物受壓實影響狀況不同，機械壓實導致土壤密度增加進而玉米干物質含量降低，與大豆相比，土壤堅實度對玉米根系發(fā)育影響更明顯[7]；Tesfaye等研究發(fā)現(xiàn)，隨壓實次數(shù)增加，土壤體積密度與滲透阻力均呈非線性增長，壓實地塊甘蔗產量較未壓實地塊出現(xiàn)大幅下降[8]；齊紅志等通過研究機械壓實作業(yè)，對玉米植株生長發(fā)育影響發(fā)現(xiàn)，0～20 cm深度土壤受大型農業(yè)機械壓實作業(yè)影響最大，耕層土壤緊實度增加對玉米根系起到抑制作用，影響莖、葉生長，最終使玉米減產幅度達20%以上[9]。

各國學者針對機械壓實問題開展大量研究，但大多基于簡單統(tǒng)計分析與線性擬合，未建立機械壓實對作物產量影響預測模型，無法科學指導機械壓實問題消減研究。隨著人工智能興起，機器學習算法發(fā)展迅速。相較于傳統(tǒng)模型，機器學習算法可深入挖掘數(shù)據(jù)內部關系，具有更強數(shù)據(jù)擬合能力與預測精度[10]。結合農業(yè)生產機械化實際引入機器學習方法以探究機械作業(yè)造成土壤壓實對黑土區(qū)大豆產量影響具有重要理論意義。

隨機森林（Random forest，RF）是基于Bagging集成學習思想分類器，具有訓練樣本數(shù)量需求少、分類精度高等優(yōu)點，在作物識別、地塊分類、機械故障診斷等領域應用廣泛[11-12]，但壓實對作物產量影響預測尚缺乏研究。本文利用隨機森林實現(xiàn)不同深度土壤堅實度對大豆產量影響預測，針對隨機森林超參數(shù)選擇主觀性強、泛化能力差問題，使用遺傳優(yōu)化算法（Genetic algorithm，GA）對隨機森林進行改進，構建基于GA-RF土壤堅實度對大豆產量影響預測模型，可有效識別因受機械壓實影響而導致大豆減產嚴重地塊，對提高大豆品質，保護土壤資源，促進農業(yè)機械化向全程全面高質高效發(fā)展具有重要意義。

1 數(shù)據(jù)獲取與處理

1.1 試驗區(qū)概況

試驗地點位于黑龍江省哈爾濱市東北農業(yè)大學向陽試驗示范基地（44°04′N，125°42′E）。試驗區(qū)處于第一積溫帶，氣候為中溫帶大陸季風性氣候，年平均溫度3.5～4.5℃，大于10℃有效積溫為2 600～2 700℃，年均降雨量400～600 mm，集中在7～8月，無霜期為135～145 d，土壤類型為典型黑土，土壤質地為粉壤土（砂粒5.3%，粉粒68.5%，黏粒28.9%）。

1.2 試驗設計

針對當前黑土區(qū)大豆生產作業(yè)中多種經營規(guī)模并存、不同規(guī)模農業(yè)機械同時使用情況，設計基于不同類型拖拉機、不同壓實次數(shù)的梯度壓實試驗。在春適期播種大豆，選用凱斯210大型拖拉機、約翰迪爾904中型拖拉機、約翰迪爾280小型拖拉機分別進行2、4、6、8、10、12次壓實試驗，以模擬不同機械化程度土壤受壓實情況。為盡可能消除地塊差異對試驗結果的影響，每組壓實處理設置3次重復，設置未經壓實地塊為對照組（CK），試驗過程中控制拖拉機行駛速度保持勻速且一致。

土壤因機械作業(yè)導致壓實情況可通過土壤堅實度、土壤含水率、土壤容重等進行表征[13]，其中土壤堅實度又稱土壤滲透阻力，是基于土壤硬度合成指標，由柱塞壓入土壤時所受土壤摩擦力、壓縮力和抗剪力等構成，對作物種子發(fā)芽、幼苗出土、根系生長影響較大；相較于其他土壤物理特性指標，土壤堅實度通過土壤堅實度儀獲取，簡單易行，故選用不同深度土壤堅實度衡量土壤受機械壓實作業(yè)影響情況。壓實完成，在每個試驗區(qū)拖拉機壓實輪轍中線處隨機選取3個測點，使用由荷蘭Eijkelkamp公司生產PV6.08土壤堅實度測定儀測定0～80 cm土層深度土壤堅實度（MPa），該儀器每1 cm自動記錄一個土壤堅實度數(shù)值。土壤堅實度測量時，雙手握持PV6.08型貫穿阻力測試儀勻速插入土壤，以避免因受力不均對測量結果產生干擾。

大豆成熟后，采集壓實輪轍兩側2 m2大豆樣本，脫粒，測定籽粒含水率，計算標準含水率[14]大豆產量（kg·hm-2）。

土壤堅實度測量與大豆樣本采集現(xiàn)場，分別如圖1和2所示。

圖1 土壤堅實度測量Fig.1 Soil penetration resistance measurement

圖2 大豆樣本采集Fig.2 Soybean sample collection

1.3 數(shù)據(jù)處理

為探究不同深度層次土壤堅實度對大豆產量影響，將采集的0～80 cm共80個土壤堅實度數(shù)據(jù)按層次進行劃分，計算各試驗小區(qū)0～10，11～20，…，71～80 cm土壤堅實度平均值作為輸入特征xi，i=1，…，N，N為樣本個數(shù)，本試驗共54個試驗小區(qū)，每個試驗小區(qū)設置3處土壤堅實度測定點，N值為162。經計算得CK組大豆產量為2 529.95 kg·hm-2，試驗組大豆產量為1 816.48～2 518.85 kg·hm-2，較CK組產量下降0.44%～28.20%，即不同機械壓實作業(yè)組合均對大豆造成減產影響。依據(jù)式（1）計算各試驗組同CK組大豆產量變化率ci（%），將其按照表1離散化為正常、輕度影響、較嚴重影響、嚴重影響4類狀況作為待預測類別yi，以分類評估不同深度土壤堅實度對大豆產量影響狀況[15]。

表1 類別劃分Table 1 Category division

式（1）中，mi為地塊i大豆產量；m′為CK組大豆產量。

為消除由地塊差異及人為因素帶來的誤差，將原始數(shù)據(jù)基于3σ法則作異常值處理。當土壤堅實度數(shù)據(jù)滿足式（2）時視為正常值，否則視為異常值，予以剔除，式中σ為標準差，μ為均值。

異常值剔除導致數(shù)據(jù)集中部分數(shù)據(jù)出現(xiàn)缺失，對此使用拉格朗日插值法進行插補[16]。拉格朗日插值法可依據(jù)現(xiàn)有數(shù)據(jù)點分布構造插值函數(shù)L(x)，實現(xiàn)對未知樣本點預測。對有n+1個取值點且分布互異土壤堅實度數(shù)據(jù)D={(x0,y0),(x1,y1),…,(xn,yn)而言，設x為土層深度，y為其對應土壤堅實度值，對于其中存在缺失值求解過程如下：

①根據(jù)已知土壤堅實度土層數(shù)據(jù)求解拉格朗日基函數(shù)lj(x)：

②根據(jù)式（3）中拉格朗日基函數(shù)構造拉格朗日插值多項式：

③將缺失函數(shù)值對應點帶入式（4）求解得到對未知點土壤堅實度值L(x)估計值，完成缺失值填補工作。

原始數(shù)據(jù)經異常值與缺失值處理后，依照大豆受影響狀況類別整理試驗數(shù)據(jù)，結果如表2所示。

表2 試驗數(shù)據(jù)Table 2 Experimental data

2 研究方法

2.1 皮爾遜相關系數(shù)

皮爾遜相關系數(shù)（Pearson correlation coefficient）又稱皮爾遜積矩相關系數(shù)，通過待評估的兩組數(shù)據(jù)與各自離差相乘反映兩組數(shù)據(jù)之間相關程度，其計算公式如式（5）所示。

式（5）中，X、Y為待確定相關程度的兩組數(shù)據(jù)，分別有n個元素；和分別為兩組數(shù)據(jù)平均值。皮爾遜相關性系數(shù)取值介于-1與1之間，取值為正時代表X與Y為正相關，值為負時代表負相關。

使用皮爾遜相關系數(shù)法衡量不同深度土壤堅實度與大豆產量之間相關性，篩選不同大豆產量關聯(lián)度高的特征作為分類模型輸入。

2.2 改進隨機森林（GA-RF）模型

改進隨機森林通過引入遺傳算法對隨機森林模型優(yōu)化。遺傳算法也稱進化算法，通過不斷選擇、交叉、變異等操作，演化產生近似解，具有全局搜索能力[17]。利用遺傳算法對隨機森林中決策樹個數(shù)、決策樹最大深度、葉子節(jié)點含有最小樣本數(shù)、節(jié)點可分最小樣本數(shù)4個關鍵參數(shù)進行尋優(yōu)[18]，解決隨機森林中超參數(shù)設置盲目性、泛化性能差問題，實現(xiàn)土壤堅實度對大豆產量影響預測。

圖3為GA-RF模型執(zhí)行過程，步驟如下：

圖3 GA-RF模型流程Fig.3 GA-RF model flow

①將采集到的土壤堅實度與大豆產量數(shù)據(jù)進行預處理，按照3∶1隨機劃分訓練集和驗證集，使用訓練集數(shù)據(jù)用于模型建立；

②初始化遺傳算法種群參數(shù)，設置隨機森林中待優(yōu)化超參數(shù)取值范圍，采用實數(shù)編碼方式隨機產生一組初始參數(shù)編碼，用于模型初始化；

③以模型分類誤差為適應度函數(shù)，計算種群中個體適應度；

④根據(jù)適應度選擇個體，淘汰適應度低個體，保留適應度高個體作交叉、變異等操作，產生下一代；

⑤重復步驟③～④，直至滿足終止條件（迭代次數(shù)或停滯代數(shù)達到設定值），輸出隨機森林中決策樹個數(shù)、決策樹最大深度、葉子節(jié)點含有最小樣本數(shù)、節(jié)點可分的最小樣本數(shù)作為最優(yōu)超參數(shù)組合，建立GA-RF模型；

⑥使用構建的GA-RF模型對驗證集進行驗證，輸出預測結果。

2.3 模型評價

選擇準確率（Accuracy）以及宏平均（Macro average）后查準率（Precision）、召回率（Recall）、F1值（F-Measure）作為評價標準[19]，衡量模型分類效果，計算公式分別如式（6）～式（9）所示。

式（6）～（9）中，n為類別數(shù)，N為數(shù)據(jù)集樣本數(shù)，TP為正確分類為正例的正樣本數(shù)，F(xiàn)P為錯誤分類為正例的負樣本數(shù)，F(xiàn)N為錯誤分類為負例的正樣本數(shù)，TN為正確分類為負例的負樣本數(shù)。

準確率為正確分類樣本占總樣本比例；查準率為所有預測結果為正例中真正為正例比例；召回率為所有為正例中被正確預測的樣本比例；F1值為查準率和召回率的加權調和平均，綜合評價模型性能；宏平均即在多分類問題中計算每一類樣本性能指標的算術平均值作為最終評價結果，上述評價指標接近于1代表模型分類性能越好。

3 結果與分析

3.1 不同深度土壤堅實度與黑土區(qū)大豆產量相關性分析

試驗共采集162組數(shù)據(jù)，依據(jù)3σ法則與拉格朗日插值法對數(shù)據(jù)異常值及缺失值進行處理，用皮爾遜相關系數(shù)法計算不同深度土壤堅實度與黑土區(qū)大豆產量相關系數(shù)并進行相關性排序，結果如表3所示。

由表3可知，不同深度土壤堅實度與大豆產量之間均為負相關，其中負相關性最強為11～20 cm土壤堅實度，相關系數(shù)為-0.935；0～10 cm與21～30 cm土壤堅實度相關系數(shù)分別為-0.909和-0.872，不同大豆產量之間存在極強負相關；31～40 cm、41～50 cm與51～60 cm土壤堅實度與大豆產量相關系數(shù)絕對值為0.6～0.8，為強相關；61～70 cm與71～80 cm等深層次土壤堅實度與大豆產量相關性較低。選取與黑土區(qū)大豆產量之間存在強相關性土壤堅實度特征，作為GA-RF模型輸入?yún)⒘?，以降低模型復雜度，提高計算效率。

表3 不同深度土壤堅實度與大豆產量相關性Table 3 Correlation between soil penetration resistance at different depths and soybean yield

3.2 基于GA-RF模型土壤堅實度對黑土區(qū)大豆產量影響預測

3.2.1 基于GA的RF參數(shù)尋優(yōu)

根據(jù)皮爾遜相關系數(shù)分析結果，選擇0～10 cm、11～20 cm、21～30 cm、31～40 cm、41～50 cm、51～60 cm 6個與黑土區(qū)大豆產量相關性強的土壤堅實度特征，作為預測模型輸入?yún)⒘?，預測大豆產量受影響狀況。將預處理及特征選擇數(shù)據(jù)按3∶1隨機作訓練集與驗證集劃分，輸入訓練集數(shù)據(jù)到GARF模型中用于模型建立。使用GA對RF超參數(shù)尋優(yōu)時，GA種群規(guī)模設置為20，交叉變異概率設置為0.1，終止條件為最大迭代次數(shù)為150或停滯代數(shù)為10；RF中待確定的4個超參數(shù)尋優(yōu)范圍設置與最終尋優(yōu)結果見表4。

表4 GA-RF尋優(yōu)結果Table 4 GA-RF optimization results

根據(jù)表4尋優(yōu)結果可知，當RF中決策樹個數(shù)設置為11、決策樹最大深度為7，葉節(jié)點含有最小樣本數(shù)為2、節(jié)點可分最小樣本數(shù)為5時，GA-RF模型性能最優(yōu)。

3.2.2 GA-RF與RF對黑土區(qū)大豆產量影響預測性能比較

將訓練集數(shù)據(jù)分別應用于GA-RF與RF模型中進行訓練，使用驗證集數(shù)據(jù)驗證土壤堅實度對黑土區(qū)大豆產量影響狀況預測模型性能，繪制RF與GA-RF模型訓練集與驗證集預測結果與真實類別對比曲線，如圖4和5所示，圖中，X軸為樣本編號，Y軸為樣本類別，實線為大豆產量受影響狀況真實類別，虛線為預測類別。

圖4 RF預測結果Fig.4 RF prediction results

計算GA-RF與RF模型在訓練集與驗證集上準確率及宏平均后查準率、召回率、F1值等評價指標并作對比，結果如表5所示。

圖5 GA-RF預測結果Fig.5 GA-RF prediction results

表5 GA-RF與RF模型性能對比Table 5 Comparison between GA-RF and RF model

分析圖4、5及表5可知，RF與GA-RF模型在訓練集中均有較好預測性能，對大豆產量受影響類別預測準確率分別為95.04%、97.52%；但在應用驗證集時，RF模型準確率較訓練集出現(xiàn)大幅下降，GA-RF模型準確率下降幅度較小，具有更好的泛化性與魯棒性。GA-RF模型在驗證集中宏平均后準確率、查準率、召回率和F1值，較RF模型分別提高7.31%、8.47%、7.50%和8.15%，在實現(xiàn)土壤堅實度對大豆產量影響狀況分類預測時，表現(xiàn)更佳。

在實際生產應用中，更重要的是對受壓實作業(yè)影響導致產量大幅下降樣本的識別，比較圖4和5兩種模型預測結果可知，RF對正常類別樣本預測性能較好，但針對其他受影響類別樣本識別能力弱，對于嚴重影響類別其識別準確率僅為75%，這將對實際生產造成嚴重誤判，不利于后續(xù)壓實緩解方案實施。結合改進GA算法的RF模型對不同類別樣本區(qū)分度更好，其在驗證集中對嚴重影響類別預測準確率達到100%，可有效識別因壓實作業(yè)使大豆產量受嚴重影響的樣本，可為機械壓實消減提供正確指導。

3.2.3 GA-RF與其他機器學習模型對黑土區(qū)大豆產量影響預測性能比較

為進一步評估GA-RF模型在土壤堅實度對大豆產量影響預測性能上優(yōu)越性，分別使用支持向量機（Support vector machines，SVM）、K近鄰（Knearest neighbors，KNN）、決策樹（Decision tree，DT）和梯度提升樹（Gradient boosting decision tree，GBDT）等常用機器學習算法進行建模并同GA-RF驗證集結果進行對比分析，結果見圖6。

圖6 不同模型驗證集結果對比Fig.6 Comparison of results of different model validation sets

由圖6可見，在驗證集中，KNN模型表現(xiàn)最差，其準確率、宏平均查準率、召回率、F1值分別為80.49%、77.99%、75.42%、76.52%；SVM與DT準確率均為85.37%，但在查準率與F1值上SVM略優(yōu)于DT；GBDT在準確率上和傳統(tǒng)RF模型接近，但召回率僅為81.67%，整體性能不如隨機森林模型。本文提出改進后隨機森林模型性能較SVM、KNN、DT、GBDT更為優(yōu)異，其準確率及宏平均查準率、召回率和F1值分別為95.12%、95.97%、96.25%、95.96%，較其他模型分別提高7.31%、5.63%、9.38%和10.91%，在對大豆產量受土壤堅實度影響狀況預測上具有更好識別能力與魯棒性。

4 討論

目前針對機械壓實對作物產量影響研究多基于簡單線性擬合，缺乏相關預測模型構建，通過GA-RF模型將機器學習方法引入機械壓實領域，可有效評估壓實風險，為機械化調整方案提供指導。GA-RF模型在土壤堅實度對黑土區(qū)大豆產量影響預測上準確率達95.12%，較傳統(tǒng)RF模型提高7.31%，有效克服傳統(tǒng)隨機森林模型參數(shù)選擇主觀性、泛化能力弱的問題，可有效識別受壓實影響地塊，評估壓實風險，為機械壓實對作物影響預測提供新思路。

皮爾遜相關分析結果表明，表層（0～30 cm）土壤堅實度增加對黑土區(qū)大豆帶來減產效應最為明顯，與Bottinelli等研究結果一致[20]。在實際生產作業(yè)中，農機壓實效應直接作用于表層土壤，使表層土壤堅實度增幅較大，阻礙作物根系生長發(fā)育，導致作物減產。因此，在黑土區(qū)大豆機械化生產過程中，應做好表層土壤疏松工作，避免在土壤含水量較高時進行機械作業(yè)[21]，同時通過增施有機肥等方式改善土壤結構，提高土壤抗壓能力。

試驗過程中樣本數(shù)據(jù)量較少，試驗結果存在一定局限性，導致分析結果可能存在偏差；構建土壤堅實度對大豆產量影響預測模型，僅針對單一年份數(shù)據(jù)進行訓練及驗證，對不同年份機械壓實作業(yè)對大豆產量影響以及不同受影響程度樣本壓實緩解措施需進一步研究。

5 結論

本文以黑土區(qū)大豆為研究對象進行不同梯度土壤壓實試驗，通過皮爾遜相關系數(shù)評估不同深度土壤堅實度與大豆產量之間相關性，結合遺傳算法對隨機森林超參數(shù)作優(yōu)化，構建基于GA-RF模型土壤堅實度對大豆產量影響預測模型，結果表明：

a.所測8個土層深度土壤堅實度與大豆產量均呈負相關；不同深度層次土壤堅實度對大豆產量影響效應不同。表層0～30 cm土壤堅實度對大豆產量影響最大，其相關性超過0.8；71～80 cm土壤堅實度與大豆產量相關性最小。

b.GA-RF模型引入遺傳算法對傳統(tǒng)RF模型超參數(shù)作優(yōu)化，克服傳統(tǒng)RF模型中魯棒性差，泛化能力低的問題；采用GA-RF模型實現(xiàn)土壤堅實度對黑土區(qū)大豆產量影響預測準確率達95.12%，較傳統(tǒng)RF模型提高7.31%，可更好預測土壤堅實度對黑土區(qū)大豆產量影響狀況。

c.GA-RF模型實現(xiàn)土壤堅實度對大豆產量影響預測效果較好，其宏平均后查準率、召回率和F1值分別為95.97%、96.25%和95.96%，較支持向量機、K近鄰、決策樹和梯度提升樹等常用機器學習算法中最優(yōu)模型，分別提高5.63%、9.38%和10.91%，分類性能更佳。