基于隨機(jī)森林回歸的玉米單產(chǎn)估測(cè)

王鵬新 齊 璇 李 俐 王 蕾 許連香

(1.中國(guó)農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院， 北京 100083; 2.農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)災(zāi)害遙感重點(diǎn)實(shí)驗(yàn)室， 北京 100083;3.中國(guó)農(nóng)業(yè)大學(xué)土地科學(xué)與技術(shù)學(xué)院， 北京 100083)

0 引言

作物長(zhǎng)勢(shì)的動(dòng)態(tài)監(jiān)測(cè)及產(chǎn)量的準(zhǔn)確估測(cè)，能夠?yàn)檗r(nóng)業(yè)經(jīng)營(yíng)者的田間管理和國(guó)家糧食政策的制定提供有效支撐[1-2]。近年來(lái)，隨著遙感技術(shù)的迅速發(fā)展，大范圍、多維空間的作物長(zhǎng)勢(shì)監(jiān)測(cè)和產(chǎn)量估測(cè)成為可能。目前，經(jīng)驗(yàn)回歸模型是作物產(chǎn)量估測(cè)的常用方法之一[3]。

經(jīng)驗(yàn)回歸模型通常選取與作物產(chǎn)量密切相關(guān)的特征參數(shù)進(jìn)行估產(chǎn)。在此類研究中，植被指數(shù)(Vegetation index, VI)應(yīng)用廣泛[4]。任建強(qiáng)等[5]以美國(guó)玉米為研究對(duì)象，以各州為估產(chǎn)區(qū)，通過(guò)篩選的歸一化植被指數(shù)(Normalized difference vegetation index, NDVI)與玉米單產(chǎn)間的最佳估產(chǎn)模型對(duì)2011年各州玉米單產(chǎn)進(jìn)行了估算，并推算全國(guó)玉米單產(chǎn)，結(jié)果表明，全國(guó)玉米單產(chǎn)的相對(duì)誤差僅為2.12%。王愷寧等[6]選取Landsat 8 OLI衛(wèi)星遙感數(shù)據(jù)，計(jì)算冬小麥灌漿期歸一化植被指數(shù)、比值植被指數(shù)(Ratio vegetation index, RVI)、綠度植被指數(shù)(Greenness vegetation index, GVI)和增強(qiáng)植被指數(shù)(Enhanced vegetation index, EVI)4種植被指數(shù)，并與冬小麥單產(chǎn)建立單植被指數(shù)和多植被指數(shù)的神經(jīng)網(wǎng)絡(luò)和SVM模型，結(jié)果表明，多植被指數(shù)SVM模型的估產(chǎn)精度高于神經(jīng)網(wǎng)絡(luò)模型。LIAQAT等[7]以巴基斯坦整個(gè)印度河流域?yàn)檠芯繀^(qū)域，通過(guò)多種植被指數(shù)，如土壤調(diào)整植被指數(shù)(Soil adjusted vegetation index, SAVI) 和改良土壤調(diào)整植被指數(shù)(Modified soil adjusted vegetation index, MSAVI)等，與小麥單產(chǎn)建立逐步回歸模型，結(jié)果表明SAVI與小麥單產(chǎn)的決定系數(shù)R2和Pearson相關(guān)系數(shù)分別為0.74和0.88。然而，作物單產(chǎn)除與植被指數(shù)相關(guān)外，還與土壤含水率和生長(zhǎng)狀態(tài)密切相關(guān)[8]。因此，可通過(guò)綜合作物生長(zhǎng)過(guò)程中的水分脅迫指標(biāo)和生長(zhǎng)狀態(tài)指標(biāo)提高作物單產(chǎn)估測(cè)精度。其中，條件植被溫度指數(shù)(Vegetation temperature condition index, VTCI)是基于歸一化植被指數(shù)和地表溫度(Land surface temperature，LST)的散點(diǎn)圖呈三角形的基礎(chǔ)上提出的[9]，可用于定量化地表征作物水分脅迫信息，并已成功應(yīng)用于陜西省關(guān)中地區(qū)干旱監(jiān)測(cè)及冬小麥單產(chǎn)估測(cè)預(yù)測(cè)等[10-12]。葉面積指數(shù)(Leaf area index，LAI)可表征植物的生長(zhǎng)狀態(tài)和光合作用能力，是作物長(zhǎng)勢(shì)監(jiān)測(cè)及單產(chǎn)估測(cè)的重要指標(biāo)[13]。此外，不同生育時(shí)期發(fā)生水分脅迫對(duì)作物單產(chǎn)的影響程度不同[14]，可通過(guò)賦予不同生育時(shí)期特征變量不同的權(quán)重，構(gòu)建綜合特征參數(shù)進(jìn)行作物單產(chǎn)估測(cè)以提高估測(cè)精度。王鵬新等[15]利用重采樣粒子濾波算法同化VTCI和LAI，并基于組合熵的方法構(gòu)建加權(quán)VTCI和LAI與冬小麥單產(chǎn)的線性回歸模型，結(jié)果表明不同管理模式下影響冬小麥單產(chǎn)的主要因子不同。

隨機(jī)森林(Random forest，RF)回歸模型是一種流行的機(jī)器學(xué)習(xí)模型，具有抗過(guò)擬合和預(yù)測(cè)精度高的特點(diǎn)[16-17]。應(yīng)用隨機(jī)森林回歸估測(cè)作物單產(chǎn)(尤其是通過(guò)綜合指數(shù)估測(cè)單產(chǎn))的研究相對(duì)較少。因此，本文以河北省中部平原地區(qū)為研究區(qū)域，選取條件植被溫度指數(shù)和葉面積指數(shù)為特征變量，通過(guò)隨機(jī)森林回歸算法獲取玉米主要生育時(shí)期各個(gè)特征變量的權(quán)重，進(jìn)而構(gòu)建加權(quán)特征變量與玉米單產(chǎn)間的回歸模型，以期為作物長(zhǎng)勢(shì)監(jiān)測(cè)及單產(chǎn)估測(cè)提供新思路。

1 材料與方法

1.1 研究區(qū)概況

河北省中部平原處于東經(jīng)114°32′～117°36′，北緯36°57′～39°50′之間(圖1)，包括石家莊市、保定市、廊坊市、衡水市和滄州市的部分或全部地區(qū)，包含53個(gè)縣(區(qū))。該區(qū)域?qū)倥瘻貛Т箨懶约撅L(fēng)氣候，四季分明，降水集中，是華北平原重要的農(nóng)業(yè)生產(chǎn)區(qū)之一。該地區(qū)年降水量在350～700 mm之間，且時(shí)空分布不均，降水主要集中在夏季，占全年的65%～70%，降水量由南向北逐漸減少。冬小麥-夏玉米輪作是該地區(qū)的主要耕作制度，該地區(qū)夏玉米出苗到拔節(jié)期一般在7月上旬至7月中旬、拔節(jié)到抽雄期在7月下旬至8月上旬、抽雄到乳熟期在8月中旬至9月上旬、乳熟到成熟期在9月中旬至9月下旬。通過(guò)王鵬新等[18]提出的基于時(shí)間序列葉面積指數(shù)傅里葉變換的作物種植區(qū)域提取方法提取了2010—2018年研究區(qū)域玉米種植區(qū)。

圖1 研究區(qū)域位置及玉米種植區(qū)(2010年)Fig.1 Location of study area and planting area of maize(2010)

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

1.2.1時(shí)間序列VTCI和LAI生成

選取2010—2018年每年7—9月Aqua-MODIS日地表溫度產(chǎn)品MYD11A1及日地表反射率產(chǎn)品MYD09GA，經(jīng)MRT預(yù)處理后獲得研究區(qū)域日LST和日NDVI產(chǎn)品，應(yīng)用最大值合成技術(shù)生成每年7—9月旬時(shí)間尺度的NDVI和LST最大值合成產(chǎn)品，基于多年某一旬的NDVI和LST最大值合成產(chǎn)品，運(yùn)用最大值合成技術(shù)分別生成多年的旬NDVI和LST最大值合成產(chǎn)品，基于每年7—9月旬LST最大值合成產(chǎn)品，運(yùn)用最小值合成技術(shù)生成多年的旬LST最大-最小值合成產(chǎn)品。VTCI取值范圍為0～1，其值越接近0，表明越干旱，作物受水分脅迫程度越重，其值越接近1，表明越濕潤(rùn)，作物受水分脅迫程度越輕或不受水分脅迫，VTCI計(jì)算公式為

(1)

其中

Lmax(Ni)=a+bNi

(2)

(3)

式中L(Ni)——在研究區(qū)域內(nèi)，某一像素的NDVI值為Ni時(shí)的地表溫度

Lmax(Ni)、Lmin(Ni)——研究區(qū)域當(dāng)NDVI值為Ni時(shí)所有像素地表溫度最大值和最小值

a、b、a′、b′——待定系數(shù)，由研究區(qū)域LST和NDVI的散點(diǎn)圖近似得到

選取研究區(qū)域2010—2018年每年7—9月MODIS葉面積指數(shù)產(chǎn)品MCD15A3H，該產(chǎn)品是基于Terra和Aqua衛(wèi)星上的MODIS傳感器獲得的，與MOD15A2和MYD15A2產(chǎn)品相比，MCD15A3H產(chǎn)品既有較高的時(shí)間分辨率(4 d)又有較高的空間分辨率(500 m)，有利于作物長(zhǎng)勢(shì)監(jiān)測(cè)及產(chǎn)量估測(cè)。利用MRT對(duì)產(chǎn)品進(jìn)行預(yù)處理得到研究區(qū)域葉面積指數(shù)產(chǎn)品，原始葉面積指數(shù)產(chǎn)品由于云和大氣等因素的影響存在數(shù)據(jù)驟降的現(xiàn)象，因此通過(guò)上包絡(luò)線S-G(Savitzky-Golay)濾波對(duì)原始葉面積指數(shù)產(chǎn)品進(jìn)行平滑處理[18]，經(jīng)上包絡(luò)線S-G濾波平滑處理后的葉面積指數(shù)更加符合玉米生長(zhǎng)情況。為使LAI與VTCI具有相同的時(shí)間尺度，將玉米各旬所包含的多時(shí)相LAI的最大值作為各旬的LAI值，并對(duì)上包絡(luò)線S-G濾波后的LAI進(jìn)行歸一化處理，最大值為7，最小值為0。

1.2.2VTCI和LAI計(jì)算

依據(jù)玉米4個(gè)主要生育時(shí)期的劃分，將玉米各生育時(shí)期包含的多旬VTCI和LAI的平均值作為該生育時(shí)期的VTCI和LAI值，如將7月上旬至7月中旬VTCI的平均值作為出苗到拔節(jié)期的VTCI值。再疊加研究區(qū)域行政邊界圖，將各縣(區(qū))包含的所有像素的VTCI和LAI的平均值作為該縣(區(qū))的VTCI和LAI值。以此類推，計(jì)算得到研究區(qū)域2010—2018年各縣(區(qū))玉米各生育時(shí)期的VTCI和LAI值。

1.2.3玉米單產(chǎn)數(shù)據(jù)的來(lái)源及異常數(shù)據(jù)處理

通過(guò)查閱《河北農(nóng)村統(tǒng)計(jì)年鑒》得到研究區(qū)域各縣(區(qū))2010—2016年玉米播種面積和總產(chǎn)量數(shù)據(jù)，玉米單產(chǎn)由總產(chǎn)量和播種面積計(jì)算得到。

將VTCI和LAI與玉米單產(chǎn)進(jìn)行回歸分析的殘差的置信區(qū)間在[-4 000,4 000] kg/hm2以外的單產(chǎn)數(shù)據(jù)視為異常數(shù)據(jù)，在構(gòu)建模型時(shí)將其剔除。

1.3 研究方法

隨機(jī)森林回歸對(duì)噪聲數(shù)據(jù)集容忍度較高，對(duì)高維數(shù)據(jù)集具有良好的預(yù)測(cè)能力[19-20]。它是由一組沒(méi)有聯(lián)系的回歸決策樹(shù){h(x,θk)，k=1,2,…,K}構(gòu)成的K棵集成決策樹(shù)，表示為

(4)

式中x——各縣(區(qū))玉米4個(gè)生育時(shí)期VTCI或LAI值及玉米單產(chǎn)數(shù)據(jù)

K——決策樹(shù)的數(shù)量

θk——獨(dú)立同分布隨機(jī)向量

為了提高模型的預(yù)測(cè)精度并防止出現(xiàn)過(guò)擬合情況，以隨機(jī)森林回歸算法結(jié)合袋裝法得到訓(xùn)練樣本子集，并結(jié)合隨機(jī)子空間法得到節(jié)點(diǎn)分裂特征[21]。

(1)袋裝法通過(guò)有放回地隨機(jī)抽樣，從原始樣本數(shù)據(jù)集中重復(fù)抽樣得到K個(gè)與原始樣本數(shù)據(jù)集相等的訓(xùn)練樣本N，每個(gè)訓(xùn)練樣本構(gòu)成一棵決策樹(shù)。每次進(jìn)行Bootstrap重抽樣時(shí)，未被抽中的樣本的概率為(1-1/N)N，當(dāng)N趨向于無(wú)窮大時(shí)，未被抽中樣本的概率越接近1/e，約為0.368，即原始樣本中有36.8%的數(shù)據(jù)未被抽中，這些數(shù)據(jù)被稱為袋外數(shù)據(jù)(Out of bag, OOB)，因其未參與回歸樹(shù)的構(gòu)建，故可用來(lái)估計(jì)預(yù)測(cè)袋外數(shù)據(jù)誤差(OOB誤差)及評(píng)估自變量對(duì)因變量的影響程度。另外，基于OOB預(yù)測(cè)誤差可以檢驗(yàn)?zāi)Ｐ偷姆夯芰Γ恍柙偈褂脺y(cè)試集檢驗(yàn)?zāi)Ｐ偷木?。通過(guò)袋裝法得到的K個(gè)訓(xùn)練樣本都不相同，保證了回歸樹(shù)的差異性。

(2)隨機(jī)子空間法通過(guò)袋裝法得到K棵回歸樹(shù)后，每個(gè)分裂節(jié)點(diǎn)隨機(jī)抽取所有變量(特征)中的Mtry個(gè)變量(特征)作為當(dāng)前節(jié)點(diǎn)分裂的特征子集，根據(jù)分類回歸樹(shù)(Classification and regression tree，CART)方法在特征子集中選擇最優(yōu)分裂方式進(jìn)行分裂。通過(guò)隨機(jī)子空間法得到的回歸樹(shù)具有隨機(jī)性和獨(dú)立性。在隨機(jī)森林回歸中，樹(shù)的數(shù)量K和隨機(jī)選擇的節(jié)點(diǎn)分裂變量(特征)Mtry決定著模型的預(yù)測(cè)能力。

圖3 OOB誤差隨回歸樹(shù)數(shù)量的變化曲線Fig.3 Changing curves of OOB errors with number of regression trees

基于隨機(jī)森林回歸估測(cè)玉米單產(chǎn)的流程如圖2所示。

圖2 基于隨機(jī)森林回歸估測(cè)玉米單產(chǎn)的流程圖Fig.2 Flow chart for estimating maize yield based on random forest regression

(1)將研究區(qū)域各縣(區(qū))2010—2016年玉米4個(gè)生育時(shí)期的VTCI或LAI值及玉米單產(chǎn)數(shù)據(jù)作為原始樣本(共357組數(shù)據(jù))輸入模型，通過(guò)Bootstrap重抽樣得到K個(gè)訓(xùn)練樣本子集并生成K棵回歸樹(shù)。VTCI和LAI估測(cè)玉米單產(chǎn)的OOB誤差隨樹(shù)的數(shù)量K變化曲線如圖3所示，可以看出，當(dāng)K為500時(shí)，OOB誤差趨于平穩(wěn)，故將K設(shè)為500。

(3)每棵回歸樹(shù)自上向下分裂生長(zhǎng)，直到到達(dá)某個(gè)葉子節(jié)點(diǎn)輸出估測(cè)值，所有回歸樹(shù)構(gòu)成隨機(jī)森林。將所有回歸樹(shù)輸出的玉米單產(chǎn)求平均值即可得到最終的玉米單產(chǎn)估測(cè)結(jié)果。

隨機(jī)森林回歸模型不但能精確地估測(cè)玉米單產(chǎn)，而且還可給出各個(gè)變量的重要性評(píng)分，即玉米4個(gè)生育時(shí)期VTCI或LAI對(duì)玉米單產(chǎn)的影響程度?；诨嵯禂?shù)和基于OOB誤差是常用的變量重要性評(píng)分的統(tǒng)計(jì)量，本研究中基于OOB估測(cè)誤差得到各變量的重要性。若xj(j=1,2,3,4)為輸入變量，則在第k棵樹(shù)上的重要性Ik為隨機(jī)置換變量前后袋外數(shù)據(jù)估測(cè)誤差的差值[22]。其計(jì)算公式為

(5)

變量xj在整個(gè)隨機(jī)森林中的重要性得分為

(6)

式中NOOB——袋外數(shù)據(jù)樣本數(shù)

f(xn)——袋外數(shù)據(jù)中第n個(gè)樣本值

fk(xn)、fk(x′n)——隨機(jī)置換變量前后第k棵樹(shù)上的袋外數(shù)據(jù)第n個(gè)樣本的估測(cè)值

I(·)——判別函數(shù)，當(dāng)f(xn)=fk(xn)或f(xn)=fk(x′n)時(shí)，取值為1，否則為0

由于隨機(jī)性的引入，模型每次給出的變量重要性評(píng)分略有差異，故將10次運(yùn)行結(jié)果的平均值進(jìn)行歸一化處理，作為各個(gè)變量的權(quán)重。

1.4 估產(chǎn)模型構(gòu)建

通過(guò)隨機(jī)森林回歸方法確定玉米主要生育時(shí)期VTCI和LAI的權(quán)重，計(jì)算2010—2018年各縣(區(qū))加權(quán)VTCI和LAI。對(duì)2010—2016年(除2012年，用來(lái)進(jìn)行精度驗(yàn)證)加權(quán)VTCI和LAI與玉米單產(chǎn)進(jìn)行回歸分析，選取擬合程度最優(yōu)的回歸模型對(duì)2012年各縣(區(qū))的玉米單產(chǎn)進(jìn)行估測(cè)及精度驗(yàn)證，并基于該模型逐像素估測(cè)2010—2018年研究區(qū)域的玉米單產(chǎn)。

2 結(jié)果與分析

2.1 玉米各生育時(shí)期的權(quán)重

基于隨機(jī)森林回歸模型運(yùn)行10次輸出的各變量重要性的平均值進(jìn)行歸一化處理，得到玉米各生育時(shí)期VTCI和LAI的權(quán)重(表1)?？梢钥闯?，玉米拔節(jié)—抽雄期和抽雄—乳熟期的VTCI權(quán)重較大，說(shuō)明受水分脅迫時(shí)對(duì)玉米單產(chǎn)的影響程度相對(duì)較大，主要是因?yàn)檫@兩個(gè)時(shí)期對(duì)水分脅迫較敏感，抽雄期前后發(fā)生水分脅迫會(huì)導(dǎo)致幼穗發(fā)育不良，果穗偏小，雄穗在抽出2～3 d后失去散粉能力，甚至有的雄穗不能抽出，或抽穗時(shí)間延遲，導(dǎo)致禿尖增長(zhǎng)，造成不同程度的玉米產(chǎn)量下降，水分脅迫較重的會(huì)造成雌穗部分不育甚至空稈。苗期—拔節(jié)期和乳熟—成熟期的VTCI權(quán)重相對(duì)較小，說(shuō)明發(fā)生水分脅迫對(duì)玉米單產(chǎn)的影響較小，主要是苗期發(fā)生一定程度的水分脅迫會(huì)使根向下生長(zhǎng)，有利于玉米植株后期的生長(zhǎng)發(fā)育，且后期有充足水分時(shí)能夠彌補(bǔ)之前減少的生長(zhǎng)量，乳熟期之后穗粒已經(jīng)形成，受水分影響不大[23]。LAI對(duì)玉米單產(chǎn)的影響以抽雄—乳熟期和乳熟—成熟期較大，苗期—拔節(jié)期和拔節(jié)—抽雄期較小，表明生長(zhǎng)前期LAI與玉米產(chǎn)量的相關(guān)性不大，主要是因?yàn)楣夂献饔玫漠a(chǎn)物用來(lái)進(jìn)行以根系和葉片為中心的營(yíng)養(yǎng)生長(zhǎng)，抽雄期時(shí)LAI達(dá)到最大，玉米進(jìn)入以果穗為中心的生殖生長(zhǎng)階段，LAI與產(chǎn)量的相關(guān)性開(kāi)始增大，這與姚小英等[24]的研究結(jié)果較一致。

表1 玉米各生育時(shí)期的權(quán)重結(jié)果Tab.1 Weight results of each growth stage of maize

2.2 單產(chǎn)估測(cè)模型選擇

將隨機(jī)森林回歸方法計(jì)算得到的2010—2016年(除2012年)加權(quán)VTCI和LAI與玉米單產(chǎn)基于縣域尺度進(jìn)行線性回歸分析，建立不同變量的單產(chǎn)估測(cè)模型(表2)。結(jié)果表明，基于隨機(jī)森林回歸的加權(quán)VTCI和玉米單產(chǎn)的相關(guān)性最低(R2=0.001)，且沒(méi)有通過(guò)顯著性檢驗(yàn)；加權(quán)LAI與玉米單產(chǎn)的相關(guān)性次之(R2=0.296)；加權(quán)VTCI和LAI與玉米單產(chǎn)的相關(guān)性最高(R2=0.303)，模型達(dá)極顯著水平(P<0.001)，表明VTCI和LAI與玉米單產(chǎn)呈顯著的正相關(guān)關(guān)系。因此，基于雙變量估產(chǎn)模型的精度高于單變量模型的精度?；陔S機(jī)森林回歸雙變量估產(chǎn)模型估測(cè)玉米單產(chǎn)時(shí)，玉米單產(chǎn)受LAI影響較大，VTCI影響較小，原因可能是研究區(qū)域受人為因素的影響較大，當(dāng)發(fā)生水分脅迫時(shí)，通過(guò)及時(shí)灌溉緩解了當(dāng)?shù)睾登?，致使玉米單產(chǎn)對(duì)VTCI不敏感。綜上所述，基于隨機(jī)森林回歸的雙變量估產(chǎn)模型精度最高，可用于估測(cè)研究區(qū)域2012年各縣(區(qū))的玉米單產(chǎn)。

表2 加權(quán)VTCI和LAI與玉米單產(chǎn)間的線性回歸分析Tab.2 Linear regression analysis between weighted VTCI and LAI and maize yields

2.3 2012年各縣(區(qū))玉米單產(chǎn)估測(cè)的精度評(píng)價(jià)

基于隨機(jī)森林回歸雙變量估產(chǎn)模型及2012年加權(quán)VTCI和LAI對(duì)各縣(區(qū))玉米單產(chǎn)進(jìn)行估測(cè)(表3)。玉米估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的相對(duì)誤差以清苑區(qū)最低，為0.35%，以海興縣最高，為37.10%。其中，31個(gè)縣(區(qū))玉米估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的相對(duì)誤差在10%以下，7個(gè)縣(區(qū))在10%～15%，15個(gè)縣(區(qū))在15%以上，53個(gè)縣(區(qū))的平均相對(duì)誤差為9.85%，均方根誤差為824.77 kg/hm2。個(gè)別縣(區(qū))如海興縣、鹽山縣的相對(duì)誤差較大，原因可能是海興縣、鹽山縣瀕臨渤海，土壤鹽漬化嚴(yán)重，農(nóng)業(yè)生產(chǎn)條件較差，農(nóng)田水利設(shè)施建設(shè)和機(jī)械化水平較低，不適宜種植經(jīng)濟(jì)作物，種植冬小麥和夏玉米是僅有的選擇。近年來(lái)當(dāng)?shù)匾巡扇「脑熘佧}堿地的相關(guān)措施使玉米單產(chǎn)有所提高，但是玉米生產(chǎn)仍處于較低水平，玉米單產(chǎn)被高估，從而使估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的相對(duì)誤差較大。個(gè)別縣(區(qū))如正定縣、藁城區(qū)和新樂(lè)市實(shí)際玉米單產(chǎn)較高，估測(cè)單產(chǎn)偏低，玉米單產(chǎn)被低估，原因可能是這幾個(gè)縣(區(qū))是國(guó)家糧食豐產(chǎn)科技工程河北省項(xiàng)目區(qū)的核心區(qū)，田間管理及時(shí)，玉米單產(chǎn)受人為因素影響較大。

表3 2012年各縣(區(qū))玉米估測(cè)單產(chǎn)Tab.3 Estimated yields of maize in each county (district) in 2012

為了進(jìn)一步驗(yàn)證隨機(jī)森林回歸雙變量估產(chǎn)模型的精度，基于2012年各縣(區(qū))玉米實(shí)際單產(chǎn)與估測(cè)單產(chǎn)進(jìn)行線性回歸分析。結(jié)果表明，估測(cè)單產(chǎn)與實(shí)際單產(chǎn)間呈顯著的正相關(guān)關(guān)系(P<0.001)，R2達(dá)到0.540；估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的均方根誤差為631.64 kg/hm2，進(jìn)一步說(shuō)明基于隨機(jī)森林回歸雙變量估產(chǎn)模型的精度較高，可用于研究區(qū)域玉米單產(chǎn)估測(cè)。

2.4 玉米估測(cè)單產(chǎn)的時(shí)空變化規(guī)律

圖4 基于隨機(jī)森林回歸的玉米單產(chǎn)估測(cè)結(jié)果Fig.4 Estimate results of maize yields based on random forest regression

基于隨機(jī)森林回歸雙變量估產(chǎn)模型逐像素估測(cè)2010—2018年研究區(qū)域玉米單產(chǎn)(圖4)，并逐像素統(tǒng)計(jì)玉米估測(cè)單產(chǎn)。結(jié)果表明，2010、2012、2013年玉米估測(cè)單產(chǎn)相差不大，西部地區(qū)(包括石家莊市和保定市)玉米估測(cè)單產(chǎn)在6 600 kg/hm2左右，東部地區(qū)(包括滄州市)在6 100 kg/hm2左右，南部地區(qū)(包括衡水市)在6 800 kg/hm2左右，北部地區(qū)(包括廊坊市)在6 200 kg/hm2左右；2011年玉米估測(cè)單產(chǎn)略高于2010年；2014年玉米估測(cè)單產(chǎn)略低于2013年；2015、2016、2017年玉米估測(cè)單產(chǎn)略高于2014年，西部地區(qū)在6 900 kg/hm2左右，東部地區(qū)在6 100 kg/hm2左右，南部地區(qū)在7 000 kg/hm2左右，北部地區(qū)在6 100 kg/hm2左右，2018年西部地區(qū)和南部地區(qū)玉米單產(chǎn)在7 000 kg/hm2左右，東部地區(qū)和北部地區(qū)在6 500 kg/hm2左右。以2017年為例，西部地區(qū)玉米估測(cè)單產(chǎn)為6 868 kg/hm2，東部地區(qū)為6 051 kg/hm2，南部地區(qū)為6 833 kg/hm2，北部地區(qū)為6 045 kg/hm2。研究年份間2011年玉米單產(chǎn)最高，2014年玉米單產(chǎn)較低，原因可能是2011年降水量充沛，玉米單產(chǎn)高于常年，2014年玉米生育期內(nèi)發(fā)生階段性干旱且局部地區(qū)旱情較重，玉米單產(chǎn)下降。

3 討論

課題組在陜西關(guān)中平原的冬小麥干旱監(jiān)測(cè)及單產(chǎn)估測(cè)中采用客觀賦權(quán)法如熵值法確定VTCI的權(quán)重[25]，構(gòu)建的加權(quán)VTCI和冬小麥單產(chǎn)的回歸模型精度較高，但熵值法基于指標(biāo)的差異程度確定指標(biāo)權(quán)重，異常數(shù)據(jù)對(duì)權(quán)重影響較大，且可能使權(quán)重與實(shí)際相背，因此確定冬小麥主要生育時(shí)期VTCI的權(quán)重與實(shí)際水分脅迫對(duì)冬小麥單產(chǎn)的影響程度不符。在河北省中部平原地區(qū)應(yīng)用隨機(jī)森林回歸確定玉米主要生育時(shí)期VTCI和LAI的權(quán)重，結(jié)果表明隨機(jī)森林回歸確定的VTCI權(quán)重以拔節(jié)—抽雄期、抽雄—乳熟期的權(quán)重較大，根據(jù)實(shí)際水分脅迫對(duì)玉米單產(chǎn)的影響程度[26]可以看出,基于隨機(jī)森林回歸的權(quán)重結(jié)果更加符合實(shí)際情況。主要因?yàn)楦珊祵?duì)玉米單產(chǎn)的影響具有非線性的特征，隨機(jī)森林回歸模型對(duì)于非平衡數(shù)據(jù)比較穩(wěn)健，不易受到異常值的干擾，能有效處理非線性問(wèn)題。雖然基于隨機(jī)森林回歸確定的玉米主要生育時(shí)期VTCI和LAI的權(quán)重較合理，但是未考慮農(nóng)學(xué)先驗(yàn)知識(shí)，可通過(guò)結(jié)合主觀賦權(quán)法如改進(jìn)的層次分析法進(jìn)一步修正隨機(jī)森林回歸得到的權(quán)重，使權(quán)重更加符合實(shí)際情況。另外水分脅迫也會(huì)影響玉米的生長(zhǎng)狀態(tài)，即VTCI和LAI之間可能存在多元共線性的問(wèn)題，而隨機(jī)森林回歸模型對(duì)多元共線性不敏感，可以很好地預(yù)測(cè)多個(gè)變量的作用，因此隨機(jī)森林回歸模型的精度較高。

影響玉米單產(chǎn)的因素有很多，除了受到水分脅迫和生長(zhǎng)狀態(tài)的影響外，還受到其他因素如溫度、洪澇災(zāi)害、田間管理、玉米品種等的影響。楊笛[27]通過(guò)模擬氣候變化、肥料、種植面積、灌溉和品種5個(gè)驅(qū)動(dòng)因子對(duì)黃淮海夏玉米區(qū)玉米單產(chǎn)的影響，表明肥料和品種在玉米增產(chǎn)中的作用和地位隨時(shí)間在提高，種植面積的增長(zhǎng)及灌溉系數(shù)的減少不利于玉米增產(chǎn)。通過(guò)查閱《河北農(nóng)村統(tǒng)計(jì)年鑒》可以看出,研究年份間灌溉和肥料的使用較多，這可能與研究區(qū)域玉米高產(chǎn)有一定的聯(lián)系。另外，個(gè)別年份發(fā)生災(zāi)害如2016年研究區(qū)域部分縣(區(qū))玉米苗期發(fā)生雹災(zāi)，影響玉米出苗，7月又出現(xiàn)澇災(zāi)和病蟲(chóng)害使玉米單產(chǎn)略有下降。這些因素對(duì)玉米單產(chǎn)的影響不容忽視，綜合考慮與玉米單產(chǎn)相關(guān)性較大的因素是今后研究的重點(diǎn)。

4 結(jié)論

(1)通過(guò)隨機(jī)森林回歸確定玉米主要生育時(shí)期VTCI和LAI的權(quán)重，構(gòu)建加權(quán)VTCI和LAI與玉米單產(chǎn)的單變量和雙變量估產(chǎn)模型。結(jié)果表明，基于隨機(jī)森林回歸的雙變量估產(chǎn)模型精度最高。

(2)基于隨機(jī)森林回歸雙變量估產(chǎn)模型估測(cè)2010—2018年研究區(qū)域玉米單產(chǎn)，結(jié)果表明，玉米估測(cè)單產(chǎn)在空間上的分布特征為西部地區(qū)最高、北部和南部次之、東部最低，年際間的分布特征為在波動(dòng)中呈先減少后增加的趨勢(shì)。