基于Shapley值組合預(yù)測(cè)的玉米單產(chǎn)估測(cè)

王鵬新 喬 琛 李 俐 周西嘉 許連香 胡亞京

(1.中國(guó)農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院， 北京 100083； 2.中國(guó)農(nóng)業(yè)大學(xué)土地科學(xué)與技術(shù)學(xué)院， 北京 100083)

0 引言

糧食的生產(chǎn)與安全對(duì)提升我國(guó)農(nóng)業(yè)的經(jīng)營(yíng)管理水平、完善農(nóng)作物的種植結(jié)構(gòu)和確保糧食安全等具有重要意義[1-3]，因此及時(shí)準(zhǔn)確估測(cè)農(nóng)作物的產(chǎn)量有利于維持國(guó)家穩(wěn)定和促進(jìn)經(jīng)濟(jì)發(fā)展。

干旱作為影響農(nóng)作物生長(zhǎng)發(fā)育及產(chǎn)量的重要因素已受到廣泛的關(guān)注[4-5]。王鵬新等[6]在歸一化植被指數(shù)(Normalized difference vegetation index，NDVI)和地表溫度(Land surface temperature，LST)的散點(diǎn)圖呈三角形區(qū)域分布的基礎(chǔ)上，提出了條件植被溫度指數(shù)(Vegetation temperature condition index，VTCI)的干旱監(jiān)測(cè)方法，可用于反映農(nóng)作物生長(zhǎng)過程中的水分虧缺信息。孫威等[7]對(duì)VTCI冷熱邊界的確定方法進(jìn)行了完善，并驗(yàn)證了利用 VTCI 進(jìn)行干旱監(jiān)測(cè)的可行性。農(nóng)作物產(chǎn)量不僅受到水分脅迫的影響，還與作物的生長(zhǎng)狀態(tài)密切相關(guān)[8]，葉面積指數(shù)(Leaf area index，LAI)通過對(duì)作物的光合速率和干物質(zhì)的積累量的反映能較好地表征農(nóng)作物的生長(zhǎng)狀態(tài)[9]。

遙感技術(shù)對(duì)比傳統(tǒng)統(tǒng)計(jì)調(diào)查方法，憑借其覆蓋范圍廣、重訪周期短等獨(dú)特優(yōu)勢(shì)被廣泛應(yīng)用于作物長(zhǎng)勢(shì)監(jiān)測(cè)及產(chǎn)量估測(cè)，同時(shí)對(duì)大規(guī)模的農(nóng)業(yè)生產(chǎn)調(diào)查、評(píng)價(jià)、監(jiān)測(cè)和管理具有獨(dú)特的作用[10]。隨著大數(shù)據(jù)技術(shù)的發(fā)展，遙感大數(shù)據(jù)將推動(dòng)農(nóng)業(yè)遙感估產(chǎn)的發(fā)展[11]。目前，遙感估產(chǎn)方法中常采用的統(tǒng)計(jì)模型、機(jī)理模型和半機(jī)理模型等均能夠較好地對(duì)作物進(jìn)行估產(chǎn)[12]。但由于在實(shí)際應(yīng)用中機(jī)理和半機(jī)理模型存在需要輸入較多的參數(shù)問題，因此機(jī)理和半機(jī)理模型存在一定的局限性。而統(tǒng)計(jì)模型其估產(chǎn)的精確度依賴于選取遙感影像的時(shí)相，對(duì)作物的生長(zhǎng)和產(chǎn)量形成的機(jī)理解釋性不強(qiáng)[12]，因此在實(shí)際應(yīng)用中同樣具有一定局限性。在作物生長(zhǎng)過程中，經(jīng)常受到各種因素影響，同時(shí)，這些因素在作物不同生育時(shí)期產(chǎn)生不同的影響，即使采用相同的估算方式也會(huì)得到不同的結(jié)果與估測(cè)精度。雖然其估產(chǎn)精度不同，但是這些模型可以提供不同的有用信息，如將一些方法丟棄，就會(huì)失去有用信息，從而影響估產(chǎn)的精度[13]。因此，可以將單個(gè)估測(cè)模型提供的有用信息進(jìn)行組合，以提高估產(chǎn)的精度。組合預(yù)測(cè)(Combination forecasting，CF)將不同預(yù)測(cè)模型進(jìn)行有效組合，可視為對(duì)無限逼近真實(shí)數(shù)據(jù)生成過程的有效補(bǔ)充[14]。

BATES等[15]首次提出了組合預(yù)測(cè)理論，這一預(yù)測(cè)理論在提高精度的同時(shí)更充分利用了預(yù)測(cè)樣本所表達(dá)的信息，受到了國(guó)內(nèi)外廣大學(xué)者重視。組合預(yù)測(cè)根據(jù)各單一模型的信息貢獻(xiàn)度，進(jìn)而基于計(jì)算得到的單一模型權(quán)重構(gòu)建組合預(yù)測(cè)模型，從而實(shí)現(xiàn)減少預(yù)測(cè)誤差、提高預(yù)測(cè)精度的目標(biāo)[16]。特別是在時(shí)間序列中分析真實(shí)數(shù)據(jù)生成過程中，通常具有區(qū)制轉(zhuǎn)換(Regime shift)或參數(shù)漂移等特性。組合預(yù)測(cè)方法的引入，可減少由參數(shù)錯(cuò)誤或模型錯(cuò)誤帶來的預(yù)測(cè)誤差[17]，甚至在單一預(yù)測(cè)結(jié)果存在有偏性的情況下，通過組合能產(chǎn)生具有無偏性的預(yù)測(cè)結(jié)果[18]。因此，組合模型具有普遍性的優(yōu)點(diǎn)，最終的預(yù)測(cè)結(jié)果更接近實(shí)際數(shù)據(jù)。根據(jù)組合預(yù)測(cè)方法綜合手段的不同，可分為權(quán)重綜合法和區(qū)域綜合法兩種。與區(qū)域綜合法相比，權(quán)重綜合法得到了廣泛的應(yīng)用。目前，組合預(yù)測(cè)現(xiàn)在已經(jīng)在多個(gè)領(lǐng)域內(nèi)得到了應(yīng)用[19-20]，但是在農(nóng)業(yè)領(lǐng)域還少有報(bào)道。

本文以河北中部平原為研究區(qū)域，選取與玉米長(zhǎng)勢(shì)和產(chǎn)量密切相關(guān)的VTCI和LAI為特征變量，采用極限梯度提升樹(Extreme gradient boosting，XGBoost)和隨機(jī)森林(Random forest，RF)兩種機(jī)器學(xué)習(xí)算法模型分別估測(cè)研究區(qū)域的玉米單產(chǎn)，并借鑒經(jīng)濟(jì)學(xué)中合作博弈論Shapley值利益分配理論，通過組合預(yù)測(cè)模型中的權(quán)重合成思想，以單一預(yù)測(cè)模型的均方誤差為基礎(chǔ)確定單一模型權(quán)重得到組合預(yù)測(cè)模型，以期為玉米長(zhǎng)勢(shì)監(jiān)測(cè)和產(chǎn)量估測(cè)提供新方法。

1 材料與方法

1.1 研究區(qū)概況

河北中部平原屬黃淮海平原組成部分，是我國(guó)重要的玉米種植和生產(chǎn)基地[21]，位于114°32′～117°36′E與36°57′～39°50′N之間，包括石家莊等5個(gè)市的53個(gè)縣(區(qū))，其土地面積約為5.30×104km2。河北中部平原地處典型溫帶大陸性季風(fēng)氣候區(qū)，四季分明、雨熱同期。河北省中部平原年降水量350～800 mm。降水時(shí)空分布不均，南方的降水量比北方高，夏季降水量高于冬季，豐水年降水量與枯水年降水量相差較大。該地區(qū)的耕種管理制度為一年兩熟，結(jié)合該地區(qū)玉米實(shí)際生長(zhǎng)狀況，本文將河北中部平原玉米生長(zhǎng)劃分為4個(gè)生育時(shí)期：7月上旬至7月中旬的出苗-拔節(jié)期、7月下旬至8月上旬的拔節(jié)-抽穗期、8月中旬至9月上旬的抽穗-乳熟期和9月中旬至9月下旬的乳熟-成熟期。根據(jù)王鵬新等[22]提出的作物分類方法，進(jìn)而獲得研究區(qū)域2010—2018年玉米種植區(qū)分布圖，其中河北中部平原2017年玉米種植區(qū)分布圖如圖1所示。

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

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

選取2010—2018年每年7—9月空間分辨率和時(shí)間分辨率分別為1 000 m、1 d的MODIS日地表溫度產(chǎn)品(MYD11A1)及日地表反射率產(chǎn)品(MYD09GA)。利用MRT對(duì)日地表溫度和日地表反射率產(chǎn)品進(jìn)行預(yù)處理之后，得到研究區(qū)域的日LST和日NDVI產(chǎn)品，運(yùn)用最大值合成法生成旬LST與NDVI的最大值合成產(chǎn)品；基于生成的多年某一旬的NDVI和LST最大值合成產(chǎn)品再次使用最大值合成技術(shù)生成多年的旬NDVI和LST的最大值合成產(chǎn)品；基于生成的多年某一旬的LST最大值合成產(chǎn)品利用最小值合成技術(shù)逐像素提取最小值，計(jì)算得到多年旬LST的最大—最小值合成產(chǎn)品，并以此通過計(jì)算生成河北中部平原旬VTCI時(shí)間序列數(shù)據(jù)。VTCI的計(jì)算公式為[6,23]

(1)

(2)

LSTmax(NDVIi)=a+bNDVIi

(3)

LSTmin(NDVIi)=a′+b′NDVIi

(4)

式中NDVI——?dú)w一化植被指數(shù)

ρNIR、ρred——近紅外、紅光波段的反射率

LST——地表溫度

LSTmax(NDVIi)、LSTmin(NDVIi)——當(dāng)NDVIi為某一特定值時(shí)所有像素地表溫度最大值和最小值，即熱邊界和冷邊界

LST(NDVIi)——研究區(qū)域內(nèi)某一像素的NDVI值為NDVIi時(shí)的地表溫度

a、b、a′、b′——LST和NDVI散點(diǎn)圖近似得到的待定系數(shù)

1.2.2LAI時(shí)間序列的生成

選取2010—2018年7—9月的MODIS葉面積指數(shù)產(chǎn)品(MCD15A3H)，其空間分辨率和時(shí)間分辨率分別為500 m、4 d。利用上包絡(luò)線Savitzky-Golay濾波對(duì)經(jīng)MRT處理后得到的原始LAI產(chǎn)品進(jìn)行平滑處理以消除云層覆蓋、大氣溶膠等因素引起的數(shù)據(jù)驟降的現(xiàn)象[23]，濾波處理后的葉面積指數(shù)變化趨于平穩(wěn)且更加符合玉米的生長(zhǎng)物侯特征，解決了原始數(shù)據(jù)存在的質(zhì)量問題。由于LAI與VTCI時(shí)間分辨率和取值范圍不同，因此對(duì)濾波后的葉面積指數(shù)進(jìn)行歸一化處理，并通過觀察多時(shí)相MODIS的MCD15A3H原始數(shù)據(jù)相元統(tǒng)計(jì)直方圖將取值范圍設(shè)置為0～7，進(jìn)而計(jì)算得到研究區(qū)域2010—2018年玉米主要生育時(shí)期的LAI時(shí)間序列數(shù)據(jù)。

1.2.3玉米生育時(shí)期VTCI和LAI的計(jì)算及異常點(diǎn)數(shù)據(jù)處理

基于研究區(qū)玉米生育時(shí)期的劃分結(jié)果，將生育時(shí)期內(nèi)所包含的多旬VTCI和LAI平均值作為研究區(qū)域玉米該生育時(shí)期逐像素的VTCI和LAI值；通過疊加河北中部平原各縣(區(qū))行政邊界圖，將各縣(區(qū))所包含的玉米生育時(shí)期逐像素VTCI和LAI值的平均值作為生育時(shí)期縣(區(qū))尺度的VTCI和LAI值。同時(shí)在構(gòu)建回歸模型時(shí)，剔除加權(quán)VTCI與LAI殘差置信區(qū)間在[-4 000,4 000] kg/km2以外的異常點(diǎn)數(shù)據(jù)。

1.3 研究方法

1.3.1極限梯度提升算法

極限梯度提升算法(XGBoost)是一種基于梯度提升的決策樹(Gradient boosting decision tree)集成算法[24]，其基分類器主要包含分類和回歸樹(Classification and regression tree，CART)。本文玉米單產(chǎn)估測(cè)是回歸問題，所以其基模型選擇為回歸樹。含有K棵決策樹XGBoost的樹集成模型定義為

(5)

xi——樣本所對(duì)應(yīng)的特征變量VTCI與LAI

fk——第k個(gè)決策樹的預(yù)測(cè)函數(shù)

樹集成優(yōu)化模型可以定義為

(6)

(7)

l(yi,i)=(yi-i)2

式中l(wèi)(yi,i)——損失函數(shù)，即均方誤差

Ω(f)——正則化項(xiàng)

γ——復(fù)雜度參數(shù)

T——樹中葉子節(jié)點(diǎn)個(gè)數(shù)

λ——固定系數(shù)

ω——葉子節(jié)點(diǎn)量化權(quán)重向量

求解式(6)的優(yōu)化問題，即求解CART樹的結(jié)構(gòu)。通過保留訓(xùn)練好的前t-1輪樹模型不變，在第t輪時(shí)添加一個(gè)新的預(yù)測(cè)函數(shù)，迭代計(jì)算得到最終預(yù)測(cè)結(jié)果[25]

(8)

ft(xi)——第t輪加入的新的預(yù)測(cè)函數(shù)

XGBoost模型可以通過輸出玉米4個(gè)生育時(shí)期VTCI或LAI的特征變量重要性來評(píng)估不同特征變量對(duì)玉米產(chǎn)量的影響程度。在XGBoost中常用基于增益、覆蓋度、頻率的特征重要性指標(biāo)進(jìn)行特征重要性評(píng)價(jià)，其中基于增益(gain)表示各特征變量對(duì)XGBoost模型中每棵樹采取每個(gè)特征變量的貢獻(xiàn)而計(jì)算的模型相對(duì)貢獻(xiàn)度為

(9)

式中G——VTCI與LAI基于增益獲取的特征重要性得分值

GL、GR——左、右子樹中的梯度

HL、HR——左、右子樹中的二階梯度

c——給定的臨界值

1.3.2隨機(jī)森林算法

隨機(jī)森林(RF)是一種基于袋裝法(bagging)理論實(shí)現(xiàn)的集成學(xué)習(xí)算法，其基模型是決策樹模型。隨機(jī)森林具有極好的準(zhǔn)確率，同時(shí)可以評(píng)估各特征在回歸問題上的重要性。隨機(jī)森林將若干沒有聯(lián)系的回歸決策樹{y(x,θk),k=1,2，…,K}構(gòu)成K棵集成決策樹[26]，可以定義為

(10)

式中y(x)——估測(cè)研究區(qū)域玉米的單產(chǎn)

x——特征變量的輸入，即VTCI與LAI

文獻(xiàn)[26]基于RF算法確定了河北中部平原玉米各個(gè)生育時(shí)期的權(quán)重，進(jìn)而構(gòu)建了加權(quán)VTCI和LAI與玉米單產(chǎn)之間的雙參數(shù)估產(chǎn)模型，結(jié)果顯示，雙參數(shù)估產(chǎn)模型精度相對(duì)較高，達(dá)到極顯著的水平(P<0.001)，且基于隨機(jī)森林算法確定的玉米各生育時(shí)期權(quán)重均較為合理。因此，本文引用其基于隨機(jī)森林算法確定的河北中部平原玉米各生育時(shí)期的VTCI和LAI權(quán)重進(jìn)行估產(chǎn)模型的構(gòu)建。

1.3.3組合預(yù)測(cè)模型

XGBoost算法與RF算法均為統(tǒng)計(jì)理論的機(jī)器學(xué)習(xí)方法，但二者各有優(yōu)勢(shì)。XGBoost算法通過引入正則化項(xiàng)在很大程度上避免了過擬合問題的出現(xiàn)，同時(shí)采用優(yōu)化算法降低了問題計(jì)算復(fù)雜度，在處理大規(guī)模數(shù)據(jù)集方面具有明顯的優(yōu)勢(shì)，已被廣泛應(yīng)用于不同領(lǐng)域的研究[27-29]。相對(duì)于XGBoost，RF更適合處理高維數(shù)據(jù)，對(duì)噪聲數(shù)據(jù)的容忍性較高。通過結(jié)合XGBoost算法和RF算法的優(yōu)點(diǎn)，參考經(jīng)濟(jì)學(xué)中合作博弈論的Shapley值理論，通過計(jì)算確定組合預(yù)測(cè)模型總誤差在單一預(yù)測(cè)模型之間的分布，并在此基礎(chǔ)上確定各單一預(yù)測(cè)模型的權(quán)重，具體步驟如下：

(1)求取Shapley值

(11)

(12)

式中i——預(yù)測(cè)模型序號(hào)，即XGBoost與RF預(yù)測(cè)模型序號(hào)

E′i——第i個(gè)預(yù)測(cè)模型的Shapley值，即XGBoost與RF預(yù)測(cè)模型計(jì)算得到的均方誤差

s——包含模型XGBoost與RF的所有集合

|s|——預(yù)測(cè)模型個(gè)數(shù)

E(s)——集合s的收益

E(s-{i})——集合s中去除成員i后的收益

(2)計(jì)算權(quán)重

組合預(yù)測(cè)模型中若某單一預(yù)測(cè)模型所獲得的誤差分配值越大，表示預(yù)測(cè)精度越低，在組合預(yù)測(cè)模型中的權(quán)重就越小?；诖嗽瓌t，預(yù)測(cè)模型i的權(quán)重λi定義為

(13)

最終的組合預(yù)測(cè)模型可描述為

(14)

式中Y——組合預(yù)測(cè)模型估測(cè)的玉米單產(chǎn)

Yi——模型XGBoost與模型RF估測(cè)的玉米單產(chǎn)

以通過計(jì)算得到的2010—2017年玉米4個(gè)生育時(shí)期的VTCI和LAI數(shù)據(jù)作為特征變量，相對(duì)應(yīng)的玉米單產(chǎn)數(shù)據(jù)作為目標(biāo)變量(每年53組數(shù)據(jù)，共424組數(shù)據(jù))，以2010—2017年(除2012年)數(shù)據(jù)作為訓(xùn)練集合，2012年數(shù)據(jù)作為測(cè)試集合。由于特征變量與目標(biāo)變量之間差值較大，特征變量值的范圍在0～1之間，而目標(biāo)變量值在3 000～9 000 kg/km2之間。由于特征變量和目標(biāo)變量處于不同區(qū)間的值域差異可能對(duì)結(jié)果造成不同的影響，取對(duì)數(shù)后不會(huì)改變數(shù)據(jù)的性質(zhì)和相對(duì)關(guān)系，對(duì)訓(xùn)練數(shù)據(jù)集和測(cè)試數(shù)據(jù)集中的目標(biāo)變量進(jìn)行取對(duì)數(shù)處理，并通過對(duì)估產(chǎn)模型輸出的估測(cè)數(shù)據(jù)進(jìn)行取對(duì)數(shù)還原，得到最終的估測(cè)單產(chǎn)數(shù)據(jù)。采用均方根誤差(RMSE)、平均相對(duì)誤差(MRE)和決定系數(shù)(R2)等指標(biāo)對(duì)模型的精度進(jìn)行評(píng)價(jià)。

2 結(jié)果與分析

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

構(gòu)建XGBoost估產(chǎn)模型與RF估產(chǎn)模型，借鑒Shapley值理論，計(jì)算得到XGBoost與RF估產(chǎn)模型的Shapley值。根據(jù)Shapley值確定XGBoost估產(chǎn)模型與RF估產(chǎn)模型權(quán)重，進(jìn)而完成對(duì)組合預(yù)測(cè)模型的構(gòu)建。

2.1.1XGBoost與RF的估產(chǎn)模型構(gòu)建

基于XGBoost與RF算法輸出的特征重要性值，通過歸一化計(jì)算得到研究區(qū)域玉米各生育時(shí)期VTCI與LAI的權(quán)重(表1)。可以看出，玉米生育后期的LAI權(quán)重高于玉米生育前期LAI的權(quán)重，表明玉米生育前期的LAI對(duì)玉米產(chǎn)量影響程度較弱，玉米生育后期的LAI對(duì)玉米產(chǎn)量影響程度較強(qiáng)。原因可能是玉米LAI變化規(guī)律呈前期增長(zhǎng)緩慢、中期增長(zhǎng)快速、后期下降緩慢的偏峰曲線，其中玉米在出苗-拔節(jié)期和拔節(jié)-抽雄期主要以分化莖葉的營(yíng)養(yǎng)生長(zhǎng)為主，此階段玉米生長(zhǎng)迅速，葉片迅速增多增大；以抽雄-乳熟期為界，玉米進(jìn)入以生殖生長(zhǎng)為主的生育后期，光合產(chǎn)物的分配模式主要以果穗為中心，是玉米產(chǎn)量形成的重要時(shí)期。因此玉米生育后期的LAI對(duì)玉米產(chǎn)量的影響程度強(qiáng)于玉米生育前期。拔節(jié)-乳熟期的VTCI權(quán)重高于出苗-拔節(jié)期和乳熟-成熟期的權(quán)重，表明拔節(jié)-乳熟期的VTCI對(duì)玉米產(chǎn)量的影響程度較強(qiáng)，出苗-拔節(jié)期和乳熟-成熟期的VTCI對(duì)玉米產(chǎn)量的影響程度較弱。原因可能是在拔節(jié)-抽雄期、抽雄-乳熟期玉米進(jìn)入營(yíng)養(yǎng)生長(zhǎng)的高峰期，此階段玉米生長(zhǎng)迅速，對(duì)土壤中的水分吸收也最為急迫，若發(fā)生水分脅迫將減緩玉米根莖葉的生長(zhǎng)發(fā)育，降低光合作用對(duì)玉米干物質(zhì)積累速率，減少蛋白質(zhì)和有機(jī)質(zhì)的合成，造成玉米粒重明顯降低，最終影響玉米產(chǎn)量，因此拔節(jié)-抽雄期、抽雄-乳熟期的VTCI對(duì)玉米產(chǎn)量的影響程度較強(qiáng)。而在出苗-拔節(jié)期，玉米由于植株矮小，對(duì)水分的需求量較少，在乳熟-成熟期，玉米處于生育后期，生長(zhǎng)變緩，對(duì)一定的水分脅迫表現(xiàn)出一定的忍受力，因此出苗-拔節(jié)期、乳熟-成熟期VTCI對(duì)玉米產(chǎn)量的影響程度較弱。綜上所述，基于XGBoost算法與RF算法確定的研究區(qū)域玉米各生育時(shí)期的權(quán)重均較為合理。

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

分別將XGBoost估產(chǎn)模型與RF估產(chǎn)模型得到的玉米各生育時(shí)期的特征權(quán)重進(jìn)行加權(quán)VTCI和LAI計(jì)算，進(jìn)而構(gòu)建基于加權(quán)VTCI和LAI與玉米單產(chǎn)之間的回歸估測(cè)模型，將2010—2017年(除2012年)數(shù)據(jù)代入XGBoost估產(chǎn)模型與RF估產(chǎn)模型，進(jìn)行可視化分析(圖2)?？梢钥闯?，XGBoost估產(chǎn)模型R2為0.31，均方根誤差為940.91 kg/km2，RF估產(chǎn)模型R2為0.30，均方根誤差為947.50 kg/km2，XGBoost估產(chǎn)模型與RF估產(chǎn)模型均通過顯著性檢驗(yàn)(P<0.001)。結(jié)果表明，XGBoost估產(chǎn)模型與RF估產(chǎn)模型玉米單產(chǎn)估測(cè)精度均較高，可用于研究區(qū)域各縣(區(qū))玉米的單產(chǎn)估測(cè)。

2.1.2組合預(yù)測(cè)模型的構(gòu)建

為進(jìn)一步提高玉米單產(chǎn)估測(cè)精度，進(jìn)行組合預(yù)測(cè)模型的構(gòu)建。需要確定XGBoost與RF估產(chǎn)模型的權(quán)重，基于2010—2017年(除2012年)由單一估產(chǎn)模型輸出的各縣(區(qū))的估測(cè)單產(chǎn)數(shù)據(jù)(共371組)，分別計(jì)算XGBoost估產(chǎn)模型、RF估產(chǎn)模型和組合預(yù)測(cè)模型的均方誤差，將各模型的均方誤差代入式(11)，將得到的XGBoost與RF估產(chǎn)模型的Shapley值代入式(13)，則可以得到XGBoost與RF估產(chǎn)模型的權(quán)重分別為0.48與0.52，則組合預(yù)測(cè)模型為

Y=0.48YXGBoost+0.52YRF

(15)

通過計(jì)算，組合預(yù)測(cè)模型的R2為0.32，模型通過顯著性檢驗(yàn)且達(dá)顯著水平(P<0.001)。同時(shí)，組合預(yù)測(cè)模型均方根誤差為939.81 kg/km2，與單一XGBoost估產(chǎn)模型與RF估產(chǎn)模型相比，組合預(yù)測(cè)模型決定系數(shù)與均方根誤差均得到提升。結(jié)果表明，組合模型的玉米單產(chǎn)估測(cè)精度優(yōu)于單一估產(chǎn)模型，可用于河北中部平原的玉米單產(chǎn)估測(cè)。

2.2 玉米估測(cè)單產(chǎn)的精度評(píng)價(jià)

將2012年數(shù)據(jù)分別代入XGBoost估產(chǎn)模型、RF估產(chǎn)模型與組合預(yù)測(cè)模型，將玉米的的估測(cè)單產(chǎn)與實(shí)際單產(chǎn)數(shù)據(jù)進(jìn)行數(shù)據(jù)可視化分析(圖3)。對(duì)比分析，可以看出XGBoost估產(chǎn)模型、RF估產(chǎn)模型和組合預(yù)測(cè)模型的估測(cè)單產(chǎn)與實(shí)際單產(chǎn)均呈顯著的正相關(guān)關(guān)系。同時(shí)，R2均不小于0.50，可以反映實(shí)際單產(chǎn)的波動(dòng)均有50%以上能被估測(cè)單產(chǎn)的波動(dòng)所描述，即玉米的實(shí)際單產(chǎn)與估測(cè)單產(chǎn)之間的誤差較小。其中，組合預(yù)測(cè)模型R2達(dá)到0.52，均方根誤差為831.14 kg/km2，平均相對(duì)誤差為9.86%，對(duì)比單一XGBoost估產(chǎn)模型與RF估產(chǎn)模型，組合預(yù)測(cè)模型的決定系數(shù)、均方根誤差與平均相對(duì)誤差均得到提升。進(jìn)一步表明，基于Shapley值的組合預(yù)測(cè)模型的估產(chǎn)精度優(yōu)于單一估產(chǎn)模型。

2.3 區(qū)域玉米單產(chǎn)的估測(cè)與時(shí)空變化分析

將組合估產(chǎn)模型應(yīng)用于河北中部平原2010—2018年逐像素玉米單產(chǎn)估測(cè)(圖4)。結(jié)果表明，研究區(qū)域玉米的估測(cè)單產(chǎn)隨年際變化呈現(xiàn)先減少后增加的波動(dòng)變化趨勢(shì)。其中，2010—2014年玉米估測(cè)單產(chǎn)整體呈現(xiàn)逐年下降的趨勢(shì)，并在2014年玉米平均單產(chǎn)達(dá)到最低，平均單產(chǎn)在6 000 kg/km2左右，原因可能是2014年河北中部平原降水較少且發(fā)生階段性干旱，導(dǎo)致玉米單產(chǎn)減少嚴(yán)重；2010、2011、2013年平均估測(cè)單產(chǎn)相差不大，均在6 350 kg/km2左右。對(duì)比2010—2014年玉米估測(cè)單產(chǎn)年際變化趨勢(shì)，2015—2018年的玉米單產(chǎn)整體呈逐年上升的趨勢(shì)，其中2018年玉米平均估測(cè)單產(chǎn)最高，約為6 500 kg/km2，2016年和2017年玉米平均估測(cè)單產(chǎn)次之，約為6 400 kg/km2，2015年玉米平均估測(cè)單產(chǎn)達(dá)到最低，約為6 200 kg/km2。

研究區(qū)域玉米的單產(chǎn)估測(cè)空間上呈現(xiàn)西部地區(qū)玉米估測(cè)單產(chǎn)最高，南部和北部地區(qū)玉米單產(chǎn)次之，東部地區(qū)玉米單產(chǎn)最低的分布特征。河北中部平原的北部地區(qū)中，2012年玉米的估測(cè)單產(chǎn)最高，2014年玉米估測(cè)單產(chǎn)最低，分別約為6 000 kg/km2和4 500 kg/km2，2015、2016、2017年玉米估測(cè)單產(chǎn)相差不大，均約為5 000 kg/km2，2011年和2018年玉米的估測(cè)單產(chǎn)略低于2012年；東部地區(qū)中，2016年玉米的估測(cè)單產(chǎn)最高，2014年玉米的估測(cè)單產(chǎn)最低，分別約5 000 kg/km2和4 500 kg/km2，2016—2018年東部地區(qū)的玉米估測(cè)單產(chǎn)相差不大，均約為4 700 kg/km2。南部地區(qū)玉米的估測(cè)單產(chǎn)分別在2012年和2014年達(dá)到最高和最低，分別約6 500、4 500 kg/km2，其余年份玉米估測(cè)單產(chǎn)相差不大。西部地區(qū)中，2012、2015年玉米的估測(cè)單產(chǎn)較高且相差不大，均約為6 700 kg/km2，2010、2014年玉米估測(cè)單產(chǎn)較低，均約為4 500 kg/km2。

3 討論

經(jīng)過對(duì)不同估產(chǎn)模型特點(diǎn)的研究分析發(fā)現(xiàn)，在玉米單產(chǎn)估測(cè)過程中，不同估產(chǎn)模型對(duì)特征因素的信息提取方式不完全一樣。同時(shí)以往單一估產(chǎn)模型均注重特征因素對(duì)產(chǎn)量的影響，而缺少考慮估產(chǎn)模型對(duì)特征因素的信息提取以及對(duì)估測(cè)產(chǎn)量的影響。本文在考慮VTCI和LAI作物長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)的基礎(chǔ)上，也充分考慮了估產(chǎn)模型對(duì)特征因素的信息提取以及對(duì)估測(cè)產(chǎn)量的影響，利用RF對(duì)噪聲數(shù)據(jù)的容忍性較高以及XGBoost較低的計(jì)算復(fù)雜度等優(yōu)點(diǎn)，借鑒經(jīng)濟(jì)學(xué)中Shapley值理論計(jì)算得到XGBoost和RF估產(chǎn)模型的權(quán)重，進(jìn)而構(gòu)建組合估產(chǎn)模型，實(shí)現(xiàn)了單一估產(chǎn)模型的綜合利用。在今后的研究中，可以在Shapley值理論的基礎(chǔ)上嘗試加入其他估產(chǎn)模型，以期進(jìn)一步全面地提取特征因素的信息，避免單一估產(chǎn)模型在估產(chǎn)過程中造成有用信息的浪費(fèi)，從而實(shí)現(xiàn)玉米單產(chǎn)估測(cè)精度的進(jìn)一步提升。

本文雖然選取了與玉米長(zhǎng)勢(shì)和單產(chǎn)密切相關(guān)的VTCI和LAI作為特征變量，但玉米的生長(zhǎng)發(fā)育除了受到水分脅迫和生長(zhǎng)狀態(tài)的影響之外，還受到其他自然因素和人為因素的影響，例如溫度、土壤肥力、田間管理水平等因素。因此在基于VTCI與LAI作物生長(zhǎng)指標(biāo)的基礎(chǔ)上，未來研究應(yīng)進(jìn)一步綜合考慮與玉米單產(chǎn)相關(guān)性較大的其他因素。此外，基于Shapley值構(gòu)建的組合預(yù)測(cè)模型對(duì)河北中部平原各縣(區(qū))玉米單產(chǎn)估測(cè)精度雖較高，但缺少對(duì)農(nóng)學(xué)先驗(yàn)知識(shí)的考慮，今后研究中可以通過主觀賦權(quán)法進(jìn)一步修正XGBoost估產(chǎn)模型和RF估產(chǎn)模型權(quán)重，使模型權(quán)重更優(yōu)。

4 結(jié)論

(1)通過借鑒組合預(yù)測(cè)模型思想，利用經(jīng)濟(jì)學(xué)合作博弈論Shapley值利益分配理論確定了XGBoost估產(chǎn)模型與RF估產(chǎn)模型的權(quán)重，進(jìn)而構(gòu)建組合預(yù)測(cè)模型。結(jié)果表明，基于Shapley值的組合預(yù)測(cè)模型估測(cè)精度較高，且優(yōu)于單一估產(chǎn)模型的精度。

(2)將組合預(yù)測(cè)模型應(yīng)用于研究區(qū)域2010—2018年逐像素的玉米單產(chǎn)估測(cè)。結(jié)果表明，從研究區(qū)域玉米單產(chǎn)估測(cè)的年際變化看，河北中部平原2010—2018年各縣(區(qū))的玉米估測(cè)單產(chǎn)波動(dòng)變化，但總體呈先減少后增加的趨勢(shì)；從玉米單產(chǎn)估測(cè)的空間分布看，整體呈西部地區(qū)最高，南部和北部地區(qū)次之，東部地區(qū)最低的特征。