卓 越，嚴(yán)海軍

(中國農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，北京 100083)

1 研究背景

化肥使用對于提高農(nóng)作物產(chǎn)量具有重要作用。但是長期以來，我國化肥投入結(jié)構(gòu)不合理，肥料利用率低，造成養(yǎng)分比例失調(diào)。過度施用化肥不僅增加了生產(chǎn)成本，還導(dǎo)致產(chǎn)量下降，同時對環(huán)境和農(nóng)作物造成嚴(yán)重污染[1-2]。實(shí)踐證明，精準(zhǔn)施肥可以節(jié)約肥料、增加糧食產(chǎn)量、均衡土壤養(yǎng)分、減少環(huán)境污染[3-4]，然而如何確定最佳施肥量是精準(zhǔn)施肥的難點(diǎn)。因此，探究作物產(chǎn)量與施肥量等因素之間的關(guān)系，合理構(gòu)建施肥模型、尋找最佳施肥量、實(shí)現(xiàn)精準(zhǔn)施肥是被關(guān)注的研究熱點(diǎn)之一。

施肥決策模型主要有目標(biāo)產(chǎn)量法、營養(yǎng)診斷法和肥料效應(yīng)函數(shù)法[5]。其中肥料效應(yīng)函數(shù)法是通過大量的田間試驗(yàn)獲得施肥量與產(chǎn)量的關(guān)系，進(jìn)而確定最優(yōu)施肥量，為目前廣泛應(yīng)用的一種施肥決策方法[6]。為了模擬施肥量等因素與產(chǎn)量之間的非線性關(guān)系，很多學(xué)者使用神經(jīng)網(wǎng)絡(luò)的方法構(gòu)建模型。馬成林等[7]先采用數(shù)據(jù)包絡(luò)分析法對數(shù)據(jù)進(jìn)行預(yù)處理，再通過BP神經(jīng)網(wǎng)絡(luò)建立模型，最終確定了最佳施肥方案。于合龍等[8]通過bagging算法生成多個BP神經(jīng)網(wǎng)絡(luò)并通過拉格朗日乘子法進(jìn)行集成從而建立模型，指出該方法優(yōu)于常規(guī)BP神經(jīng)網(wǎng)絡(luò)，并通過模型得到了最佳施肥方案。楊曉輝等[9]分別使用模擬退火算法和遺傳算法對BP神經(jīng)網(wǎng)絡(luò)進(jìn)行優(yōu)化，進(jìn)一步提高了施肥模型的擬合精度。王福林等[10]在模型輸入中加入玉米種植密度，使用BP神經(jīng)網(wǎng)絡(luò)建模，得到了最佳種植密度與最佳施肥量模型，并進(jìn)行了模型驗(yàn)證。Dong等[11]采用小波神經(jīng)網(wǎng)絡(luò)建立施肥模型，指出該方法的擬合精度優(yōu)于隨機(jī)森林回歸與支持向量回歸，并使用模型確定了最佳施肥量。

以上研究表明，使用BP神經(jīng)網(wǎng)絡(luò)或優(yōu)化BP神經(jīng)網(wǎng)絡(luò)可以描述施肥量等因素與作物產(chǎn)量之間的關(guān)系，從而制定施肥方案。然而BP神經(jīng)網(wǎng)絡(luò)存在兩大缺陷。第一，BP神經(jīng)網(wǎng)絡(luò)是一種個體學(xué)習(xí)器，有性能提升的瓶頸[12]。第二，BP神經(jīng)網(wǎng)絡(luò)是一種“黑箱”模型，無法評估每個特征對輸出結(jié)果的影響程度，模型的可解釋性差[10,13]。梯度提升樹是一種常用的集成學(xué)習(xí)算法，該算法訓(xùn)練多個個體學(xué)習(xí)器，再通過一定的策略結(jié)合，從而形成一個強(qiáng)學(xué)習(xí)器，具有擬合精度高、解釋性強(qiáng)等特點(diǎn)[14-17]。使用梯度提升樹算法建立的模型能有效地解決上述問題。本文以玉米“3414”試驗(yàn)數(shù)據(jù)作為訓(xùn)練樣本，通過插值算法對數(shù)據(jù)進(jìn)行擴(kuò)充，采用梯度提升樹算法建立施肥模型，并與常見的幾種建模方法進(jìn)行對比。通過分析施肥模型求解出最大產(chǎn)量與最優(yōu)施肥量，以期提供施肥指導(dǎo)。

2 材料與方法

2.1 數(shù)據(jù)集描述與處理

“3414”試驗(yàn)是國內(nèi)普遍采用的研究肥料效應(yīng)的田間試驗(yàn)方案[18]，其中“3”指氮、磷、鉀3個因素，“4”指不施肥、最佳施肥量的0.5倍、最佳施肥量和最佳施肥量的1.5倍4種施肥水平，“14”指共有14種處理。

本文建模使用數(shù)據(jù)來自吉林省榆樹市10處玉米地的“3414”試驗(yàn)[19]。在試驗(yàn)區(qū)土壤類型、氣候等條件基本一致時，影響作物產(chǎn)量的因素為土壤氮、磷、鉀含量和氮、磷、鉀肥施用量這6個因素。依據(jù)當(dāng)?shù)氐膶?shí)際生產(chǎn)情況和專家經(jīng)驗(yàn)，得到氮、磷、鉀的最佳施用量分別為180、75和75 kg/hm2。表1列舉了10個試驗(yàn)區(qū)的土壤養(yǎng)分含量，表2列舉了試驗(yàn)區(qū)1的“3414”試驗(yàn)處理與對應(yīng)的產(chǎn)量。

表1 各試驗(yàn)區(qū)的土壤養(yǎng)分含量 mg/kg

表2 試驗(yàn)區(qū)1的14種施肥處理與產(chǎn)量 kg/hm2

由表2可以分析玉米施肥量與產(chǎn)量的關(guān)系。第2、3、6、11組數(shù)據(jù)中磷肥和鉀肥均處于最佳施用量，隨著氮肥施肥量的增加，產(chǎn)量呈先增加后減小的趨勢。觀察第4、5、6、7組數(shù)據(jù)，在氮肥和鉀肥處于最佳施用量時，隨著磷肥施肥量的增加，產(chǎn)量也呈現(xiàn)先增加后減小的趨勢。觀察第6、8、9、10組數(shù)據(jù)，對鉀肥也獲得相同結(jié)果。這種現(xiàn)象與實(shí)際情況相符，表明通過“3414”數(shù)據(jù)構(gòu)建施肥模型是合理的。

表2提供的施肥量與產(chǎn)量數(shù)據(jù)在構(gòu)建施肥模型時樣本數(shù)不夠，為此，可以通過插值算法在一定范圍內(nèi)推求出新的數(shù)據(jù)點(diǎn)，進(jìn)行數(shù)據(jù)集擴(kuò)充[11]。常用的插值算法有線性插值法、多項(xiàng)式插值法和樣條插值法。其中線性插值法快速簡單但精度較差，而且在插值點(diǎn)處斜率會發(fā)生變化；多項(xiàng)式插值法是線性插值法的推廣，精度有所提高，但是在區(qū)間邊緣容易出現(xiàn)振蕩現(xiàn)象；樣條插值法使用特殊分段多項(xiàng)式進(jìn)行插值，可以避免振蕩問題，精度較高。為了獲得更好的插值效果，本文采用2次樣條插值法。通過表2中的第2、3、6、11組數(shù)據(jù)可以建立試驗(yàn)區(qū)1的氮肥施肥量與產(chǎn)量的插值曲線，其結(jié)果如圖1所示。從插值曲線上兩個相鄰的原始數(shù)據(jù)點(diǎn)之間選擇4個插值點(diǎn)從而可以得到16組數(shù)據(jù)。使用同樣的方法對磷肥施肥量、鉀肥施肥量進(jìn)行插值。去掉重復(fù)的數(shù)據(jù)再加上原始數(shù)據(jù)，每個試驗(yàn)區(qū)可以擴(kuò)充為50組數(shù)據(jù)，對其余9個試驗(yàn)區(qū)的數(shù)據(jù)進(jìn)行相同的處理，可以將數(shù)據(jù)擴(kuò)充為500組。

圖1 試驗(yàn)區(qū)1氮肥施肥量與產(chǎn)量的插值曲線 圖2 不同學(xué)習(xí)率下模型迭代次數(shù)與決定系數(shù)的關(guān)系曲線

2.2 梯度提升樹

梯度提升樹(Gradient boosting decision tree，GBDT)是由Friedman[20]于2001年提出的一種集成學(xué)習(xí)算法。其主要思想是每次建立的新模型都基于上一個模型損失函數(shù)的負(fù)梯度，通過多個弱學(xué)習(xí)器合成為一個強(qiáng)學(xué)習(xí)器。當(dāng)弱學(xué)習(xí)器為回歸樹時，其計(jì)算過程如下[21]。

步驟1：輸入訓(xùn)練數(shù)據(jù)集D={(x1,y1),(x2,y2), … ,(xN,yN)}。模型的輸出為F(x)，損失函數(shù)為L(y,F(x))。損失函數(shù)L的種類很多，常見的有平方差損失函數(shù)、絕對損失函數(shù)、Huber損失函數(shù)等。在梯度提升樹算法中通常使用平方差損失函數(shù)。

L(y,F(x))=(y-F(x))2

(1)

步驟2：初始化模型F0(x)。

(2)

式中：γ為葉子結(jié)點(diǎn)輸出值。

步驟3：對m= 1, 2, … ,M進(jìn)行M次迭代，總共生成M個回歸樹。

(3)

(3)計(jì)算回歸樹Tm葉節(jié)點(diǎn)的最佳輸出值。

(4)

(4)更新模型

Fm(x)=Fm-1(x)+ν·γjm(x∈Rjm)

(5)

公式(5)中通常加入學(xué)習(xí)率ν來控制模型學(xué)習(xí)的速度，即每次更新Fm(x)之前，把葉子節(jié)點(diǎn)的輸出乘以學(xué)習(xí)率ν(0<ν≤1)，以小的步長逐漸逼近最佳結(jié)果。

步驟4：迭代結(jié)束，生成模型。

2.3 模型構(gòu)建

在構(gòu)建模型之前首先要劃分?jǐn)?shù)據(jù)集，一般選取2/3～4/5的樣本數(shù)據(jù)用于訓(xùn)練，剩余樣本用于測試[22]。梯度提升樹模型中很多參數(shù)需要調(diào)整，為了評估模型在不同參數(shù)下的效果，需要從訓(xùn)練集中選取一部分?jǐn)?shù)據(jù)作為驗(yàn)證集，用于模型調(diào)整模型參數(shù)[23]。由于本試驗(yàn)的數(shù)據(jù)量較少，使用單獨(dú)劃分出的驗(yàn)證集進(jìn)行參數(shù)調(diào)整不具有代表性，因此采用交叉驗(yàn)證的方式調(diào)整模型的參數(shù)。最終將500組數(shù)據(jù)隨機(jī)分為400組訓(xùn)練數(shù)據(jù)和100組測試數(shù)據(jù)，以土壤氮、磷、鉀含量和氮、磷、鉀肥施用量為輸入量，以玉米產(chǎn)量為輸出量。采用5折交叉驗(yàn)證和參數(shù)搜尋的方法調(diào)整模型參數(shù)。

在調(diào)節(jié)模型參數(shù)的過程中，通過score函數(shù)計(jì)算模型的決定系數(shù)R2來評判模型的優(yōu)劣。首先調(diào)節(jié)迭代次數(shù)M和學(xué)習(xí)率ν。這兩個參數(shù)均可控制模型的擬合程度，改變其中一個參數(shù)會影響另一個參數(shù)的最佳值。通過參數(shù)搜尋的方式尋找迭代次數(shù)M與學(xué)習(xí)率ν的最優(yōu)組合，不同學(xué)習(xí)率下迭代次數(shù)與模型決定系數(shù)的關(guān)系如圖2所示。由圖2可以看出，隨著迭代次數(shù)的增加，模型的精度越來越高；學(xué)習(xí)率較小時需要更大的迭代次數(shù)才能使模型達(dá)到較高的精度。迭代次數(shù)過少會導(dǎo)致模型的欠擬合，過多會導(dǎo)致模型的過擬合并且增加計(jì)算時間，綜合考慮選擇迭代次數(shù)M=100，學(xué)習(xí)率ν=0.1。同樣采用參數(shù)搜尋的方式可以確定回歸樹的最大深度為7，葉子節(jié)點(diǎn)最少樣本數(shù)為5，最大特征數(shù)為3。

3 結(jié)果與分析

3.1 預(yù)測結(jié)果與模型解釋

確定模型的最佳參數(shù)組合后使用全部訓(xùn)練數(shù)據(jù)重新訓(xùn)練模型，訓(xùn)練結(jié)束后使用測試集檢驗(yàn)?zāi)Ｐ?。為了便于觀察，從100組測試數(shù)據(jù)中選取30組數(shù)據(jù)，將預(yù)測值與真實(shí)值進(jìn)行對比，結(jié)果如圖3所示。由圖3可以看到產(chǎn)量的預(yù)測值與實(shí)際值基本吻合，只有少部分測試數(shù)據(jù)出現(xiàn)了偏差，表明模型的預(yù)測效果較好，準(zhǔn)確度較高。

圖3 部分測試集模型預(yù)測值與真實(shí)值對比 圖4 模型各輸入變量的相對重要度

梯度提升樹模型通過每個變量在回歸樹中出現(xiàn)的次數(shù)來計(jì)算每個輸入變量的重要度[15]，從而使模型具有一定的解釋性。圖4展示了每個輸入變量的相對重要度。可以看出，相對重要度最高的為土壤的含磷量SP，其次分別為施氮量FN、土壤含氮量SN、土壤含鉀量SK、施鉀量FK、施磷量FP。從氮、磷、鉀各元素的角度來看，氮元素和磷元素對產(chǎn)量的影響較大，而鉀元素對產(chǎn)量的影響較小。

當(dāng)然，通過職業(yè)院校技能大賽，也反映出了我們在教學(xué)中的一些薄弱環(huán)節(jié)。如教學(xué)投入不足，教學(xué)實(shí)習(xí)和實(shí)訓(xùn)設(shè)備不夠完善，選手不能適應(yīng)競賽中采用的現(xiàn)代企業(yè)新設(shè)備、新技術(shù)、新流程，或在規(guī)定時間內(nèi)完成不了比賽任務(wù)；基礎(chǔ)理論課教學(xué)與專業(yè)技能訓(xùn)練沒有有機(jī)結(jié)合。應(yīng)大力推行教學(xué)做一體化模式，使車間與教室合二為一，理論與實(shí)踐有機(jī)融合，努力培養(yǎng)更多高素質(zhì)、技能型專業(yè)人才和實(shí)踐應(yīng)用型能工巧匠。

3.2 模型比較

除了梯度提升樹之外，BP神經(jīng)網(wǎng)絡(luò)(back propagation neural network, BPNN)、支持向量回歸(support vector regression, SVR)、隨機(jī)森林(random forest, RF)也是解決非線性問題的有力工具。為了進(jìn)一步驗(yàn)證梯度提升樹模型的效果，分別使用以上3種機(jī)器學(xué)習(xí)算法建立模型，并與梯度提升樹算法建立的模型進(jìn)行對比。為了便于比較，使用各方法建模時，均以土壤氮、磷、鉀含量和氮、磷、鉀肥施用量作為輸入變量，使用產(chǎn)量作為輸出變量，并使用相同的訓(xùn)練集進(jìn)行訓(xùn)練。

與梯度提升樹對比的3種算法中，BP神經(jīng)網(wǎng)絡(luò)使用3層前饋網(wǎng)絡(luò)，隱含層個數(shù)確定為11[24]，模型中加入L2正則化項(xiàng)防止過擬合，使用雙曲正切激活函數(shù)，用牛頓法進(jìn)行迭代，建模之前對數(shù)據(jù)進(jìn)行歸一化處理；支持向量回歸模型中引入RBF核函數(shù)來解決非線性問題，使用訓(xùn)練集交叉驗(yàn)證和網(wǎng)格搜索的方式最終確定懲罰系數(shù)C=115，核函數(shù)系數(shù)為0.1；隨機(jī)森林模型通過交叉驗(yàn)證逐步調(diào)整模型參數(shù)，最終得到回歸樹個數(shù)Mt=250，回歸樹的最大深度為19，最大特征數(shù)為5。

采用相對誤差(RE)、均方根誤差(RMSE)和平均絕對誤差(MAE)作為模型的評價(jià)指標(biāo)。圖5分別展示了4種模型在相同測試集上的相對誤差。由圖5(a)中可以看，出梯度提升樹模型的RE主要在0～1%的范圍內(nèi)波動，有少部分測試數(shù)據(jù)的RE在1%～2%，平均相對誤差為0.46%。圖5(b)中BP神經(jīng)網(wǎng)絡(luò)模型的RE主要在0～1%的范圍內(nèi)波動，有少部分測試數(shù)據(jù)的RE在1%～2%，有極少測試數(shù)據(jù)的RE為2%～3%，平均相對誤差為0.54%。圖5(c)中支持向量回歸模型的RE波動較大，波動范圍主要在1%～6%，平均相對誤差為3.19%。圖5(d)中隨機(jī)森林模型的RE主要在0～2%的范圍內(nèi)波動，少部分測試數(shù)據(jù)的RE位于2%～6%，平均相對誤差為1.00%。可以看出4種模型中，支持向量回歸模型的相對誤差最大，在數(shù)值上明顯高于其他3種模型；而梯度提升樹模型的相對誤差最小。

圖5 4種模型在相同測試集上的相對誤差

圖6展示了4種模型(GBDT、BPNN、SVR和RF)的均方根誤差與平均絕對誤差。其中梯度提升樹模型的RMSE和MAE分別為62.2和48.7kg/hm2，BP神經(jīng)網(wǎng)絡(luò)模型的RMSE和MAE分別為78.5和56.5 kg/hm2，支持向量回歸模型的RMSE和MAE分別為371.6和337.5 kg/hm2，隨機(jī)森林模型的RMSE和MAE分別為133.8和104.6 kg/hm2。支持向量回歸模型的RMSE和MAE明顯高于其他3種模型，梯度提升樹、BP神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林3種模型的RMSE和MAE較小，其中梯度提升樹模型的RMSE和MAE最小。

圖6 4種模型的均方根誤差與平均絕對誤差

以上結(jié)果表明梯度提升樹算法建模效果最優(yōu)，BP神經(jīng)網(wǎng)絡(luò)和隨機(jī)森林次之，支持向量回歸最差。

3.3 確定最優(yōu)施肥量

依據(jù)梯度提升樹算法構(gòu)建的施肥模型可表示為：

Y=F(SN,SP,SK,FN,FP,FK)

(6)

式中：Y為產(chǎn)量，kg/hm2；SN、SP、SK分別為土壤中氮、磷、鉀的含量，mg/kg；FN、FP、FK分別為氮、磷、鉀肥施用量，kg/hm2；F為產(chǎn)量與土壤氮、磷、鉀含量和氮、磷、鉀肥施用量的函數(shù)關(guān)系。

因此在給定土壤養(yǎng)分含量的情況下，通過求解非線性規(guī)劃問題可以計(jì)算最大產(chǎn)量以及相應(yīng)的施肥量。以試驗(yàn)區(qū)1為例，已知：

(1)Y=F(SN,SP,SK,FN,FP,FK)

(2)SN= 102 mg/kg，SP= 45 mg/kg，SK= 156 mg/kg

(3)0

經(jīng)過計(jì)算，可得試驗(yàn)區(qū)1的最大產(chǎn)量與相應(yīng)的最佳施肥量，即試驗(yàn)區(qū)1的最佳施氮量為193 kg/hm2，施磷量為80 kg/hm2，施鉀量為73 kg/hm2，此時最大產(chǎn)量為10 161 kg/hm2。

使用相同的方法可以計(jì)算其他9個試驗(yàn)區(qū)的最佳施肥量與產(chǎn)量。表3列舉了各試驗(yàn)區(qū)的最優(yōu)施肥量與產(chǎn)量。

表3 各試驗(yàn)區(qū)的最優(yōu)施肥量與產(chǎn)量 kg/hm2

由于梯度提升樹算法基于回歸樹，因此在一定土壤含量范圍內(nèi)可能得出相同的最佳施肥量。然而各個試驗(yàn)區(qū)的土壤養(yǎng)分含量不同，因此其最大產(chǎn)量也有所不同。由表3可知，試驗(yàn)區(qū)7得到的產(chǎn)量最大，達(dá)到13 242 kg/hm2。

4 討 論

合理的施肥方案既可以提高作物的產(chǎn)量，同時也能減少環(huán)境污染。通過建立施肥模型尋找施肥量等因素與產(chǎn)量之間的關(guān)系，從而指導(dǎo)施肥是實(shí)現(xiàn)精準(zhǔn)農(nóng)業(yè)的關(guān)鍵。本文提出了一種基于梯度提升樹算法的施肥模型，得到了適合當(dāng)?shù)氐淖罴咽┓柿糠桨浮?/p>

在4種建模方法的對比中，梯度提升樹算法建模效果最優(yōu)，BP神經(jīng)網(wǎng)絡(luò)和隨機(jī)森林其次，最后是支持向量回歸，這和Dong等[11]的研究結(jié)果十分接近。梯度提升樹算法和隨機(jī)森林算法均為基于決策樹的集成算法，建模精度較高。兩種算法的不同之處在于梯度提升樹的基學(xué)習(xí)器是串行生成，即每次生成的新學(xué)習(xí)器都是依據(jù)上一次的建模結(jié)果；隨機(jī)森林算法中的基學(xué)習(xí)器是并行生成，基學(xué)習(xí)器之間是獨(dú)立的。BP神經(jīng)網(wǎng)絡(luò)模型雖然也有較高的建模精度，但是相比其他3種方法，在建模過程中需要更多參數(shù)調(diào)整，并且不同的參數(shù)組合對建模結(jié)果有較大的影響，相比兩種基于決策樹的建模方法缺少可解釋性。支持向量回歸雖然可以通過引入核函數(shù)的方式解決非線性的問題，但是在擬合的精度上不如其他3種方法。通過模型計(jì)算得到各試驗(yàn)區(qū)的最佳施肥量盡管與當(dāng)?shù)氐耐扑]量稍有不同，但在產(chǎn)量上有所提高，能夠給當(dāng)?shù)厥┓侍峁┘夹g(shù)指導(dǎo)。

本文提出的施肥模型考慮了土壤養(yǎng)分含量和總施肥量，因此只適用于土壤類型、氣候等條件基本一致的地區(qū)，存在一定的局限性。為了擴(kuò)大施肥模型的適用范圍，使其具有更好的泛用性，今后應(yīng)該在更廣泛的尺度上收集數(shù)據(jù)，并且在建模時考慮更多的影響因素。

5 結(jié) 論

(1)使用梯度提升樹算法建立的施肥模型可以反映土壤養(yǎng)分含量、施肥量與作物產(chǎn)量之間的關(guān)系，并對產(chǎn)量進(jìn)行較高精度預(yù)測。

(2)對比4種建模方法發(fā)現(xiàn)，梯度提升樹模型最優(yōu)，BP神經(jīng)網(wǎng)絡(luò)和隨機(jī)森林模型次之，支持向量回歸模型最差，在今后實(shí)際應(yīng)用中可以優(yōu)先采用梯度提升樹算法進(jìn)行建模。

(3)相比BP神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林和支持向量回歸模型，梯度提升樹算法建立的施肥模型具有更好的解釋能力。通過分析本文建立的施肥模型發(fā)現(xiàn)，影響產(chǎn)量較大的因素是土壤含磷量與施氮量，鉀元素對產(chǎn)量的影響較小。

(4)在已知土壤氮、磷、鉀養(yǎng)分含量情況下，由施肥模型可以得到最優(yōu)施肥量方案和最大產(chǎn)量，從而有效指導(dǎo)施肥。