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

(1.中國農業(yè)大學信息與電氣工程學院， 北京 100083; 2.農業(yè)農村部農業(yè)災害遙感重點實驗室， 北京 100083;3.中國農業(yè)大學土地科學與技術學院， 北京 100083)

0 引言

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

1.2 數據獲取與處理

1.2.1時間序列VTCI和LAI生成

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

(1)

其中

Lmax(Ni)=a+bNi

(2)

(3)

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

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

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

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

1.2.2VTCI和LAI計算

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

1.2.3玉米單產數據的來源及異常數據處理

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

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

1.3 研究方法

隨機森林回歸對噪聲數據集容忍度較高，對高維數據集具有良好的預測能力[19-20]。它是由一組沒有聯系的回歸決策樹{h(x,θk)，k=1,2,…,K}構成的K棵集成決策樹，表示為

(4)

式中x——各縣(區(qū))玉米4個生育時期VTCI或LAI值及玉米單產數據

K——決策樹的數量

θk——獨立同分布隨機向量

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

(1)袋裝法通過有放回地隨機抽樣，從原始樣本數據集中重復抽樣得到K個與原始樣本數據集相等的訓練樣本N，每個訓練樣本構成一棵決策樹。每次進行Bootstrap重抽樣時，未被抽中的樣本的概率為(1-1/N)N，當N趨向于無窮大時，未被抽中樣本的概率越接近1/e，約為0.368，即原始樣本中有36.8%的數據未被抽中，這些數據被稱為袋外數據(Out of bag, OOB)，因其未參與回歸樹的構建，故可用來估計預測袋外數據誤差(OOB誤差)及評估自變量對因變量的影響程度。另外，基于OOB預測誤差可以檢驗模型的泛化能力，不需再使用測試集檢驗模型的精度。通過袋裝法得到的K個訓練樣本都不相同，保證了回歸樹的差異性。

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

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

基于隨機森林回歸估測玉米單產的流程如圖2所示。

圖2 基于隨機森林回歸估測玉米單產的流程圖Fig.2 Flow chart for estimating maize yield based on random forest regression

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

(3)每棵回歸樹自上向下分裂生長，直到到達某個葉子節(jié)點輸出估測值，所有回歸樹構成隨機森林。將所有回歸樹輸出的玉米單產求平均值即可得到最終的玉米單產估測結果。

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

(5)

變量xj在整個隨機森林中的重要性得分為

(6)

式中NOOB——袋外數據樣本數

f(xn)——袋外數據中第n個樣本值

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

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

由于隨機性的引入，模型每次給出的變量重要性評分略有差異，故將10次運行結果的平均值進行歸一化處理，作為各個變量的權重。

1.4 估產模型構建

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

2 結果與分析

2.1 玉米各生育時期的權重

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

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

2.2 單產估測模型選擇

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

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

2.3 2012年各縣(區(qū))玉米單產估測的精度評價

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

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

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

2.4 玉米估測單產的時空變化規(guī)律

圖4 基于隨機森林回歸的玉米單產估測結果Fig.4 Estimate results of maize yields based on random forest regression

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

3 討論

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

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

4 結論

(1)通過隨機森林回歸確定玉米主要生育時期VTCI和LAI的權重，構建加權VTCI和LAI與玉米單產的單變量和雙變量估產模型。結果表明，基于隨機森林回歸的雙變量估產模型精度最高。

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