基于Shapley值組合預測的玉米單產估測

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

(1.中國農業大學信息與電氣工程學院， 北京 100083； 2.中國農業大學土地科學與技術學院， 北京 100083)

0 引言

糧食的生產與安全對提升我國農業的經營管理水平、完善農作物的種植結構和確保糧食安全等具有重要意義[1-3]，因此及時準確估測農作物的產量有利于維持國家穩定和促進經濟發展。

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

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

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

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

1 材料與方法

1.1 研究區概況

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

1.2 數據獲取與處理

1.2.1VTCI時間序列的生成

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

(1)

(2)

LSTmax(NDVIi)=a+bNDVIi

(3)

LSTmin(NDVIi)=a′+b′NDVIi

(4)

式中NDVI——歸一化植被指數

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

LST——地表溫度

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

LST(NDVIi)——研究區域內某一像素的NDVI值為NDVIi時的地表溫度

a、b、a′、b′——LST和NDVI散點圖近似得到的待定系數

1.2.2LAI時間序列的生成

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

1.2.3玉米生育時期VTCI和LAI的計算及異常點數據處理

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

1.3 研究方法

1.3.1極限梯度提升算法

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

(5)

xi——樣本所對應的特征變量VTCI與LAI

fk——第k個決策樹的預測函數

樹集成優化模型可以定義為

(6)

(7)

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

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

Ω(f)——正則化項

γ——復雜度參數

T——樹中葉子節點個數

λ——固定系數

ω——葉子節點量化權重向量

求解式(6)的優化問題，即求解CART樹的結構。通過保留訓練好的前t-1輪樹模型不變，在第t輪時添加一個新的預測函數，迭代計算得到最終預測結果[25]

(8)

ft(xi)——第t輪加入的新的預測函數

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

(9)

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

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

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

c——給定的臨界值

1.3.2隨機森林算法

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

(10)

式中y(x)——估測研究區域玉米的單產

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

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

1.3.3組合預測模型

XGBoost算法與RF算法均為統計理論的機器學習方法，但二者各有優勢。XGBoost算法通過引入正則化項在很大程度上避免了過擬合問題的出現，同時采用優化算法降低了問題計算復雜度，在處理大規模數據集方面具有明顯的優勢，已被廣泛應用于不同領域的研究[27-29]。相對于XGBoost，RF更適合處理高維數據，對噪聲數據的容忍性較高。通過結合XGBoost算法和RF算法的優點，參考經濟學中合作博弈論的Shapley值理論，通過計算確定組合預測模型總誤差在單一預測模型之間的分布，并在此基礎上確定各單一預測模型的權重，具體步驟如下：

(1)求取Shapley值

(11)

(12)

式中i——預測模型序號，即XGBoost與RF預測模型序號

E′i——第i個預測模型的Shapley值，即XGBoost與RF預測模型計算得到的均方誤差

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

|s|——預測模型個數

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

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

(2)計算權重

組合預測模型中若某單一預測模型所獲得的誤差分配值越大，表示預測精度越低，在組合預測模型中的權重就越小。基于此原則，預測模型i的權重λi定義為

(13)

最終的組合預測模型可描述為

(14)

式中Y——組合預測模型估測的玉米單產

Yi——模型XGBoost與模型RF估測的玉米單產

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

2 結果與分析

2.1 估產模型構建

構建XGBoost估產模型與RF估產模型，借鑒Shapley值理論，計算得到XGBoost與RF估產模型的Shapley值。根據Shapley值確定XGBoost估產模型與RF估產模型權重，進而完成對組合預測模型的構建。

2.1.1XGBoost與RF的估產模型構建

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

表1 玉米各生育時期的權重Tab.1 Weight results of each growth stage of maize

分別將XGBoost估產模型與RF估產模型得到的玉米各生育時期的特征權重進行加權VTCI和LAI計算，進而構建基于加權VTCI和LAI與玉米單產之間的回歸估測模型，將2010—2017年(除2012年)數據代入XGBoost估產模型與RF估產模型，進行可視化分析(圖2)。可以看出，XGBoost估產模型R2為0.31，均方根誤差為940.91 kg/km2，RF估產模型R2為0.30，均方根誤差為947.50 kg/km2，XGBoost估產模型與RF估產模型均通過顯著性檢驗(P<0.001)。結果表明，XGBoost估產模型與RF估產模型玉米單產估測精度均較高，可用于研究區域各縣(區)玉米的單產估測。

2.1.2組合預測模型的構建

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

Y=0.48YXGBoost+0.52YRF

(15)

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

2.2 玉米估測單產的精度評價

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

2.3 區域玉米單產的估測與時空變化分析

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

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

3 討論

經過對不同估產模型特點的研究分析發現，在玉米單產估測過程中，不同估產模型對特征因素的信息提取方式不完全一樣。同時以往單一估產模型均注重特征因素對產量的影響，而缺少考慮估產模型對特征因素的信息提取以及對估測產量的影響。本文在考慮VTCI和LAI作物長勢監測指標的基礎上，也充分考慮了估產模型對特征因素的信息提取以及對估測產量的影響，利用RF對噪聲數據的容忍性較高以及XGBoost較低的計算復雜度等優點，借鑒經濟學中Shapley值理論計算得到XGBoost和RF估產模型的權重，進而構建組合估產模型，實現了單一估產模型的綜合利用。在今后的研究中，可以在Shapley值理論的基礎上嘗試加入其他估產模型，以期進一步全面地提取特征因素的信息，避免單一估產模型在估產過程中造成有用信息的浪費，從而實現玉米單產估測精度的進一步提升。

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

4 結論

(1)通過借鑒組合預測模型思想，利用經濟學合作博弈論Shapley值利益分配理論確定了XGBoost估產模型與RF估產模型的權重，進而構建組合預測模型。結果表明，基于Shapley值的組合預測模型估測精度較高，且優于單一估產模型的精度。

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