王鵬新 胡亞京 李 俐 喬 琛

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

0 引言

葉面積指數(shù)(Leaf area index，LAI)[1]是評估籽粒潛在產(chǎn)量的重要指標，對作物長勢監(jiān)測與產(chǎn)量估測、預測均具有重要意義。然而，站點實測、遙感反演或作物模型模擬等單一方法都無法給出精度可靠且時空連續(xù)的LAI數(shù)據(jù)。數(shù)據(jù)同化算法[2-3]能夠?qū)⑦b感觀測和作物生長模型有效結(jié)合，是大面積作物產(chǎn)量估測、預測領(lǐng)域的研究熱點。文獻[4]將通過ENVISAT ASAR和MERIS數(shù)據(jù)提取得到LAI與CERES-Wheat模擬得到LAI進行同化，結(jié)果表明，同化后LAI有效提高了冬小麥產(chǎn)量預測精度。農(nóng)作物產(chǎn)量是多種環(huán)境參數(shù)共同作用的結(jié)果，如土壤含水率、氣象條件等，因此基于多參數(shù)的農(nóng)作物估產(chǎn)已經(jīng)成為研究趨勢。文獻[5]將CERES-Maize模擬的LAI、SM與遙感觀測的LAI、SM采用改進的EnKF算法進行同化，結(jié)果表明，基于同化LAI和SM構(gòu)建的雙參數(shù)產(chǎn)量估測模型精度更高。

干旱是我國最主要的農(nóng)業(yè)氣象災(zāi)害之一，影響著農(nóng)業(yè)的可持續(xù)發(fā)展和國家糧食安全，有效地進行干旱監(jiān)測對作物長勢和產(chǎn)量估測至關(guān)重要[6]。文獻[7]提出了基于條件植被溫度指數(shù)(Vegetation temperature condition index，VTCI)的干旱監(jiān)測方法，該方法被證明適用于陜西省關(guān)中平原和河北中部平原的干旱監(jiān)測和預測。文獻[8]采用EnKF算法對LAI和VTCI雙參數(shù)進行同化，以估測關(guān)中平原的冬小麥單產(chǎn)，研究表明，基于多參數(shù)構(gòu)建的估產(chǎn)模型為小麥產(chǎn)量估算提供了可靠的方法。此外，不同生育時期LAI、土壤水分等狀態(tài)參數(shù)對作物單產(chǎn)的影響程度不同，文獻[9]應(yīng)用熵值的組合預測方法計算冬小麥LAI、地上生物量和土壤含水率在不同生育時期的權(quán)重，進而對冬小麥單產(chǎn)進行估測，然而熵值法是基于指標數(shù)據(jù)的差異程度來計算權(quán)重，受異常數(shù)據(jù)的影響較大。隨機森林(Random forest, RF)[10]算法能夠有效處理非線性數(shù)據(jù)，對異常值不敏感，能夠降低異常值的干擾，從而獲得較合理的結(jié)果。文獻[11]利用隨機森林算法研究發(fā)現(xiàn)，降水、土壤肥料等變量是影響柳條稷產(chǎn)量的重要因子。文獻[12]以河南省冬小麥為研究對象，探討了隨機森林算法在冬小麥遲凍災(zāi)害氣候產(chǎn)量分類預測中的適用性，研究表明，隨機森林法可以有效應(yīng)用于區(qū)域冬小麥晚凍災(zāi)害氣候產(chǎn)量影響的預測研究。

本文以CERES-Maize模型模擬和MODIS遙感數(shù)據(jù)反演的LAI和VTCI為基礎(chǔ)數(shù)據(jù)，應(yīng)用粒子濾波算法(Particle filter，PF)對LAI和VTCI進行同化，采用隨機森林回歸算法計算同化LAI和VTCI的權(quán)重，并構(gòu)建最優(yōu)夏玉米產(chǎn)量估測模型，以期為單點和區(qū)域夏玉米產(chǎn)量估測提供重要依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)包括保定市、廊坊市、石家莊市、衡水市和滄州市的53個縣(區(qū))，覆蓋范圍為114°32′～117°36′E，36°57′～39°50′N，地處中國北方糧食生產(chǎn)基地黃淮海平原區(qū)，是典型的玉米、小麥種植區(qū)和重要的糧食生產(chǎn)基地[13]。該區(qū)位于河北省中部平原，海拔0～100 m，地勢西北高、東南低，覆蓋面積約為53萬km2。研究區(qū)屬于溫帶大陸性季風氣候，年平均氣溫4～13℃，年輻射總量4 390～5 180 MJ/m2，無霜期110～220 d，年降水量350～700 mm，降水主要集中在夏季，占全年的65%～70%。采用文獻[14]提出的作物分類方法提取夏玉米種植區(qū)域，其中，夏玉米識別的總體精度大于80%，因此可作為研究區(qū)域夏玉米單產(chǎn)估測和預測的可靠依據(jù)。通過疊加河北省行政區(qū)縣邊界圖，獲得研究區(qū)2013—2018年各縣(區(qū))的玉米種植分布圖，在研究區(qū)選取保定市雄縣雄州鎮(zhèn)南十里鋪村、保定市定州市劉早鎮(zhèn)北王莊村、保定市定興縣固城南太平莊村、衡水市冀州市南午鎮(zhèn)劉瓦窯村、衡水市景縣青蘭鄉(xiāng)南申莊村、衡水市饒陽縣城關(guān)鎮(zhèn)西直沃村、滄州市滄縣興濟鎮(zhèn)宋官屯村、滄州市泊頭縣四營鄉(xiāng)大趙屯村8個典型的夏玉米種植區(qū)作為研究樣點(圖1)。

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

1.2.1CERES-Maize模型

CERES-Maize[15]模型能夠模擬玉米生長季以天為步長的LAI、土壤剖面各層體積含水率、產(chǎn)量等，驅(qū)動模型運行需輸入包括氣象數(shù)據(jù)、土壤參數(shù)、田間管理參數(shù)和作物參數(shù)等數(shù)據(jù)。氣象數(shù)據(jù)來源于中國氣象科學數(shù)據(jù)共享服務(wù)網(wǎng)(http:∥cdc.cma.gov.cn/)，為2013—2018年河北省中部平原研究樣點的逐日最低和最高氣溫、降水量和日照時數(shù)等。模型所必需的逐日太陽輻射量Rs，則根據(jù)經(jīng)驗計算公式得[16]

(1)

式中Rmax——逐日天文輻射量，MJ/m2

as、bs——經(jīng)驗系數(shù)[16]，as取0.18，bs取0.55

w——逐日日照時數(shù)

W——逐日可照時數(shù)(最大時長)

土壤數(shù)據(jù)來源于文獻[17]，主要包括土壤類型、容重、田間持水量、飽和水含量、有機質(zhì)含量、有效氮含量、pH值等。田間管理數(shù)據(jù)主要包括夏玉米品種、播種方式、高度、植株密度、行距、施肥日期和施肥量、灌溉日期和灌溉量均來自文獻[18-19]。

作物遺傳特性參數(shù)包括出苗到幼苗期結(jié)束所需的積溫、光周期敏感系數(shù)，抽絲到收獲成熟所需的積溫、單株最大穗粒數(shù)、最大灌漿速率參數(shù)以及一片葉生長所需積溫。作物遺傳參數(shù)的標定采用模型自帶的調(diào)參程序包獲取[20]。

1.2.2旬LAI的生成

選取2013—2018年玉米主要生育期(7—9月)MODIS葉面積指數(shù)產(chǎn)品(MCD15A3H)，利用MODIS數(shù)據(jù)處理工具MRT進行拼接、重采樣、投影轉(zhuǎn)換和裁剪等預處理，得到研究區(qū)LAI產(chǎn)品。采用上包絡(luò)線S-G(Savitzky-Golay)濾波[14]對原始LAI進行平滑處理，以消除云、大氣等因素引起LAI數(shù)據(jù)驟降現(xiàn)象。將玉米每旬所包含的多時相LAI的最大值作為每旬的LAI值，得到旬尺度LAI[21]。根據(jù)8個研究樣點的經(jīng)緯度坐標計算其在MODIS影像上的像素坐標，并以該像素為中心，3像素×3像素范圍內(nèi)所含像素的LAI均值作為玉米主要生育期單點同化的觀測數(shù)據(jù)[22]。

由于CERES-Maize作物生長模型運行得到的是以天為步長的LAI數(shù)據(jù)，因此取研究樣點玉米主要生育時期內(nèi)每旬LAI數(shù)據(jù)的均值作為樣點旬尺度模擬LAI數(shù)據(jù)。

1.2.3旬VTCI的生成

選擇2013—2018年7—9月的MODIS日地表反射率(MYD09GA)和日地表溫度產(chǎn)品(MYD11A1)，利用MRT對原始數(shù)據(jù)進行預處理生成每天的LST產(chǎn)品和NDVI產(chǎn)品。采用最大合成技術(shù)對每天的同一像素的LST和NDVI分別進行合成，生成旬LST和NDVI最大值合成產(chǎn)品；對2013—2018年7—9月旬LST和NDVI最大值合成產(chǎn)品的同一像素LST和NDVI再進行最大值合成，生成多年旬LST和NDVI最大值合成產(chǎn)品；對2013—2018年7—9月旬LST最大值合成產(chǎn)品的同一像素LST取最小值，生成多年旬LST最大-最小合成產(chǎn)品。VTCI的計算方法為

(2)

其中

Tmax(NDVIi)=a+bNDVIi

(3)

Tmin(NDVIi)=a′+b′NDVIi

(4)

式中NDVI——歸一化植被指數(shù)

VTCI——條件植被溫度指數(shù)

T(NDVIi)——某一像素的NDVI值為NDVIi時的地表溫度

Tmax(NDVIi)——當NDVIi值等于某一特定值時所有像素地表溫度的最大值，稱為熱邊界

Tmin(NDVIi)——當NDVIi值等于某一特定值時所有像素地表溫度的最小值，稱為冷邊界

a、b、a′、b′——待定系數(shù)，通過研究區(qū)域內(nèi)NDVI和T的散點圖近似獲得

研究所需模擬VTCI不能由CERES-Maize模型直接模擬得到，文獻[23]將VTCI與野外調(diào)研實測的土壤表層(0～10 cm和10～20 cm)含水率進行相關(guān)分析，表明VTCI與土壤表層含水率有較好的相關(guān)性，因此通過模型運行獲得樣點以天為步長的土壤(0～20 cm)含水率，取2013—2018年夏玉米生育期(7—9月)內(nèi)每旬水分數(shù)據(jù)均值與觀測VTCI進行線性回歸，獲取同化所需的模擬VTCI數(shù)據(jù)[24]。

1.2.4玉米區(qū)域調(diào)查統(tǒng)計數(shù)據(jù)

本研究中各縣的玉米播種面積、總產(chǎn)量數(shù)據(jù)來源于2013—2016年的《河北農(nóng)村統(tǒng)計年鑒》，玉米單產(chǎn)數(shù)據(jù)由總產(chǎn)量和播種面積計算得到。

1.2.5總體研究技術(shù)路線

通過MODIS產(chǎn)品MCD15A3H、MYD09GA和MYD11A1獲得觀測LAI與VTCI，運行 CERES-Maize模型得到的模擬LAI數(shù)據(jù)與土壤含水率數(shù)據(jù)，并根據(jù)MODIS數(shù)據(jù)獲得觀測VTCI與土壤含水率的數(shù)據(jù)相關(guān)性得到模擬VTCI。通過PF算法得到研究樣點的同化LAI和VTCI，利用隨機森林回歸算法計算玉米4個主要生育期的LAI與VTCI權(quán)重，進而得到加權(quán)LAI與VTCI，運用玉米統(tǒng)計單產(chǎn)與加權(quán)LAI與VTCI構(gòu)建的單產(chǎn)估測模型，并選擇最優(yōu)估產(chǎn)模型估測玉米單產(chǎn)以及精度評價，研究技術(shù)路線如圖2所示。

2 研究方法

2.1 粒子濾波算法

粒子濾波[25]是一種基于蒙特卡羅仿真的濾波算法，不受線性高斯模型約束。采用一組加權(quán)粒子近似表示狀態(tài)變量的后驗概率，主要包括預測和更新2個步驟。

2.1.1預測階段

(5)

(6)

(7)

式中uk——模型驅(qū)動參數(shù)

(8)

式中θk+1——k+1時刻的觀測噪聲

2.1.2更新階段

(9)

更新粒子權(quán)重

(10)

(11)

根據(jù)粒子濾波本身的特點，隨著無效采樣粒子數(shù)據(jù)的增加，使得大量的計算浪費在對后驗概率分布幾乎不起作用的粒子上，估計性能下降[27]，因此，本文在算法實施過程中采用殘差重采樣的方法，復制高權(quán)重粒子，去除低權(quán)重粒子。

2.2 隨機森林回歸算法

(1)通過Bootstrap重抽樣的方法獲得與原始樣本數(shù)據(jù)集元素數(shù)量相等的K個訓練樣本Q，并生成K棵回歸樹{h(x,θk)}，其中原始樣本數(shù)據(jù)集包括玉米產(chǎn)量影響因子集及玉米單產(chǎn)數(shù)據(jù)，本文將2013—2016年研究區(qū)域各縣(區(qū))4個生育時期(出苗-拔節(jié)期、拔節(jié)-抽雄期、抽雄-乳熟期、乳熟-成熟期)[21]的LAI、VTCI作為夏玉米產(chǎn)量影響因子集(共212組數(shù)據(jù))，表示為

(12)

式中x——研究區(qū)域各縣(區(qū))玉米產(chǎn)量影響因子及玉米單產(chǎn)數(shù)據(jù)

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

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

(2)在RF算法中，每個分裂節(jié)點隨機抽取所有變量中的M個變量作為當前節(jié)點的特征子集進行分裂。文獻[28]試驗認為，對于回歸問題m=M/3是較好的選擇，因此本文中m=1。

(3)Bootstrap重抽樣過程中，每次未被抽中的樣本概率接近于1/e，即袋外數(shù)據(jù)(Out of bag，OOB)[29]，可用來評估自變量對因變量的影響程度，即玉米4個生育時期LAI或VTCI對產(chǎn)量的影響程度。若將xi(i=1,2,3,4)為輸入變量，則變量xi在隨機森林中的重要性V(xi)計算式為

(13)

式中xi——玉米4個生育時期的LAI或VTCI

e′xi——每棵樹對應(yīng)的袋外數(shù)據(jù)誤差率

exi——隨機置換變量xi的序列后每棵樹對應(yīng)的袋外數(shù)據(jù)誤差率

3 結(jié)果與分析

3.1 LAI與VTCI同化結(jié)果與分析

選取2015年7—9月玉米種植區(qū)的研究樣點進行粒子濾波算法的同化結(jié)果分析，以衡水市景縣青蘭鄉(xiāng)南申莊村和保定市定興縣固城南太平莊村2個樣點為例，總體觀察模擬的LAI、遙感觀測的LAI以及同化的LAI的變化趨勢如圖3所示?？梢钥闯觯珻ERES-Maize模型模擬的LAI總體高于MODIS反演的LAI，遙感反演的LAI存在明顯的偏低現(xiàn)象，這是因為混合像元等存在，且在有云污染的情況下，還存在數(shù)據(jù)質(zhì)量偏低甚至缺失的情況，但MODIS反演的LAI能夠在區(qū)域尺度上反映環(huán)境因子對作物生長綜合影響，能定量地監(jiān)測和描述作物在區(qū)域尺度上實際生長狀況。經(jīng)過PF算法的同化LAI變化趨勢更為平滑，與夏玉米實際生長變化情況更為相符。將2015年的8個研究樣點經(jīng)過PF算法的同化，VTCI和觀測VTCI分別與中國氣象科學數(shù)據(jù)共享服務(wù)網(wǎng)提供的有關(guān)降水量數(shù)據(jù)進行對比分析，同樣以上述2個樣點為例進行研究分析(圖4)?？梢钥闯?，相比于原觀測VTCI，總體上經(jīng)過PF算法調(diào)整后的旬尺度VTCI能更好地與旬累積降水量數(shù)據(jù)結(jié)合，同化VTCI對樣點水分脅迫信息反映更為敏感；同化VTCI在0.5左右，屬輕旱范圍。具體表現(xiàn)為：衡水市景縣青蘭鄉(xiāng)南申莊村樣點(圖4a)8月中旬降水量大于90 mm，而觀測VTCI在0.25左右，經(jīng)過同化后，VTCI調(diào)整為0.35左右；保定市定興縣固城南太平莊村樣點(圖4b)8月中旬和8月下旬，幾乎無明顯降水，同化VTCI分別從0.55和0.63調(diào)整為0.4和0.55左右，因此同化VTCI實驗結(jié)果對原觀測圖像整體和局部均能做出較為符合實際情況的調(diào)整。

3.2 同化結(jié)果精度驗證

為驗證單點數(shù)據(jù)同化結(jié)果與外部觀測數(shù)據(jù)的響應(yīng)關(guān)系，將2013—2018年研究區(qū)8個樣點的單點LAI和VTCI的PF同化結(jié)果與原觀測結(jié)果進行線性回歸分析(圖5)，可以看出，LAI和VTCI的同化結(jié)果和觀測數(shù)據(jù)結(jié)果較為吻合，其中，所有樣點(共432組數(shù)據(jù))同化LAI與觀測LAI的決定系數(shù)R2為0.914 2(P<0.001)，均方根誤差(RMSE)為0.467 2 m2/m2，歸一化均方根誤差(NRMSE)為6.67%。所有樣點(共432組數(shù)據(jù))同化VTCI和觀測VTCI間的決定系數(shù)R2為0.628 6(P<0.001)，均方根誤差(RMSE)為0.003 9，歸一化均方根誤差(NRMSE)為14.6%。由于所選取的8個樣點較為均勻地分布在河北省中部平原，故假設(shè)該8個研究樣點特性可以代表河北省中部平原，將同化試驗中的同化LAI和VTCI分別與遙感觀測LAI和VTCI建立線性回歸方程，由此可以得到每年玉米主要生育時期以旬為單位的區(qū)域LAI同化數(shù)據(jù)和VTCI同化數(shù)據(jù)，同化后的數(shù)據(jù)紋理更加清晰，降低了像素間驟升驟減的現(xiàn)象。

3.3 夏玉米估產(chǎn)模型構(gòu)建

依據(jù)文獻[20]計算2013—2018年研究區(qū)域各縣玉米4個生育時期的LAI和VTCI。由于本文旬尺度LAI最大值為7 m2/m2，最小值為0 m2/m2，而旬尺度VTCI范圍為0～1，為使LAI與VTCI具有相同的取值范圍，對LAI數(shù)據(jù)進行歸一化處理。根據(jù)隨機森林回歸算法獲取夏玉米各生育時期LAI和VTCI的權(quán)重(表1)?？梢钥闯觯瑹o論觀測LAI和同化LAI在各生育時期的權(quán)重均表現(xiàn)為：抽雄-乳熟期與拔節(jié)-抽雄期權(quán)重較大，乳熟-成熟期權(quán)重次之，出苗-拔節(jié)期權(quán)重較小。依據(jù)農(nóng)學知識，出苗-拔節(jié)期為玉米的營養(yǎng)生長階段，光合作用的產(chǎn)物用來進行根系和葉片為中心的營養(yǎng)生長，生長前期LAI與玉米產(chǎn)量的相關(guān)性不大，因此該時期LAI的權(quán)重最小。抽雄-乳熟期是夏玉米產(chǎn)量形成的關(guān)鍵階段，是營養(yǎng)階段向生殖生長階段過渡的核心時期，拔節(jié)-抽雄期主要是玉米營養(yǎng)生長和生殖生長并進階段，因此，LAI對玉米單產(chǎn)的以抽雄-乳熟期與拔節(jié)-抽雄期影響較大；乳熟-成熟期是光合產(chǎn)物的運輸、轉(zhuǎn)移階段，對夏玉米產(chǎn)量起著重要作用，權(quán)重略低于拔節(jié)-抽雄期。觀測VTCI和同化VTCI在各生育時期的權(quán)重均為拔節(jié)-抽雄期和抽雄-乳熟期較大，出苗-拔節(jié)期次之，乳熟-成熟期最小。拔節(jié)-抽雄期是玉米營養(yǎng)生長和生殖生長并進階段，抽雄-乳熟期是營養(yǎng)階段轉(zhuǎn)向生殖生長階段，在受到干旱脅迫時，導致雄穗雌穗花期不遇或授粉不良，進而導致出現(xiàn)空稈或結(jié)實粒數(shù)減少的情況。由于生殖生長階段受阻后不能恢復，最終導致產(chǎn)量下降，生殖生長期的干旱對作物產(chǎn)量的影響最大[30]；出苗-拔節(jié)期玉米植株受旱，但復水后能夠快速恢復生長，因此出苗-拔節(jié)期干旱對產(chǎn)量影響相對較?。蝗槭?成熟期穗粒已經(jīng)形成，水分脅迫對玉米產(chǎn)量的影響最小。綜上，隨機森林回歸算法確定的觀測LAI、同化LAI、觀測VTCI和同化VTCI的權(quán)重較為合理。

表1 夏玉米各生育時期LAI與VTCI的權(quán)重Tab.1 Weights of LAI and VTCI at key growth stages of summer maize

基于玉米各主要生育期的權(quán)重，結(jié)合LAI觀測值、LAI同化值、VTCI觀測值和VTCI同化值分別計算各縣(區(qū))的加權(quán)LAI與加權(quán)VTCI，構(gòu)建夏玉米生育期加權(quán)LAI和加權(quán)VTCI與實際單產(chǎn)之間的線性回歸方程作為夏玉米單產(chǎn)模型(表2)。結(jié)果表明，無論是單參數(shù)還是雙參數(shù)單產(chǎn)估測模型，同化參數(shù)與玉米單產(chǎn)間的擬合效果好于觀測參數(shù)與玉米單產(chǎn)間的擬合效果，雙參數(shù)單產(chǎn)估測模型的決定系數(shù)均高于單參數(shù)單產(chǎn)估測模型的決定系數(shù)，基于觀測LAI與VTCI和同化LAI與VTCI的單產(chǎn)估測模型的R2分別為0.503 3和0.615 6，且均達到極顯著水平，而同化后的雙參數(shù)單產(chǎn)估測模型R2比觀測雙參數(shù)單產(chǎn)估測模型高0.112 3。綜上，基于同化LAI與VTCI的單產(chǎn)估測模型精度更高，可作為較優(yōu)的單產(chǎn)估測模型。

表2 加權(quán)LAI和加權(quán)VTCI的玉米單產(chǎn)估測模型Tab.2 Yield estimation models of maize between weighted LAI and VTCI

3.4 估產(chǎn)結(jié)果驗證

為驗證估測模型的精度，分別利用觀測和同化雙參數(shù)單產(chǎn)估測模型對河北省中部平原2015年各縣(區(qū))的夏玉米單產(chǎn)進行估測(表3)，結(jié)果表明，基于同化LAI和VTCI雙參數(shù)單產(chǎn)估測模型的所有縣(區(qū))估測單產(chǎn)的平均相對誤差8.43%，RMSE為631 kg/hm2，NRMSE為10.50%，而基于觀測LAI和VTCI雙參數(shù)單產(chǎn)估測模型的所有縣(區(qū))估測單產(chǎn)的平均相對誤差為12.57%，RMSE為763 kg/hm2，NRMSE為12.71%?；谕疞AI和VTCI雙參數(shù)構(gòu)建的單產(chǎn)估測模型分析53個縣(區(qū))估測單產(chǎn)與統(tǒng)計單產(chǎn)相對誤差，結(jié)果表明，欒城區(qū)、高碑店市、雄縣等26個縣(區(qū))相對誤差低于10%，定州市、河間市、文安縣等12個縣(區(qū))的相對誤差均低于15%，欒城區(qū)、孟村縣等10個縣(區(qū))相對誤差在15%～25%范圍內(nèi)，青縣、獻縣、黃驊市、鹽山縣和安平縣5個縣相對誤差超過25%。表3表明基于同化LAI和VTCI的夏玉米單產(chǎn)估測模型的估產(chǎn)精度較高，作為較優(yōu)估產(chǎn)模型可以為玉米產(chǎn)量的預估提供重要信息。

表3 2015年研究區(qū)域雙參數(shù)最優(yōu)估產(chǎn)模型的估測產(chǎn)量與統(tǒng)計單產(chǎn)的絕對誤差和相對誤差Tab.3 Absolute and relative errors of estimated yield and statistical yield in bi-variables optimal yield model of study area in 2015

將2013—2018年各縣(區(qū))的夏玉米的同化LAI和VTCI值代入較優(yōu)單產(chǎn)估測模型(表2)，得到河北省中部平原玉米區(qū)域單產(chǎn)估測結(jié)果，結(jié)果表明，研究年份中2017年玉米單產(chǎn)最高，2014年玉米單產(chǎn)較低，原因可能是2017年降水量充沛，玉米單產(chǎn)高于常年，2014年玉米生育期內(nèi)發(fā)生階段性干旱且局部地區(qū)旱情較重，玉米單產(chǎn)下降[14]。從單產(chǎn)的空間分布看，西部地區(qū)(保定市、石家莊市)最高、南部地區(qū)(衡水市)和北部地區(qū)(廊坊市)次之、東部地區(qū)(滄州市)最低(圖6)。具體表現(xiàn)為2013年夏玉米估測產(chǎn)量為6 296 kg/hm2，統(tǒng)計產(chǎn)量為6 773 kg/hm2，相對誤差為7.05%；2014年夏玉米估測產(chǎn)量為5 927 kg/hm2，統(tǒng)計產(chǎn)量為6 166 kg/hm2，相對誤差為3.89%；2015年夏玉米估測產(chǎn)量為6 010 kg/hm2，統(tǒng)計產(chǎn)量為6 627 kg/hm2，相對誤差為10.27%；2016年夏玉米估測產(chǎn)量為6 185 kg/hm2，統(tǒng)計產(chǎn)量為6 286 kg/hm2，相對誤差為1.64%；2017年夏玉米估測產(chǎn)量為6 448 kg/hm2；2018年夏玉米估測產(chǎn)量為6 439 kg/hm2；以2016年為例，西部地區(qū)夏玉米估測產(chǎn)量為6 625 kg/hm2，統(tǒng)計產(chǎn)量為7 183 kg/hm2，相對誤差為7.76%；東部地區(qū)夏玉米估測產(chǎn)量為5 458 kg/hm2，統(tǒng)計產(chǎn)量為5 038 kg/hm2，相對誤差為8.32%；南部地區(qū)夏玉米估測產(chǎn)量為6 586 kg/hm2，統(tǒng)計產(chǎn)量為5 699 kg/hm2，相對誤差為15.56%；北部地區(qū)夏玉米估測產(chǎn)量為5 789 kg/hm2，統(tǒng)計產(chǎn)量為5 266 kg/hm2，相對誤差為9.94%。

4 討論

粒子濾波不受線性高斯模型的約束，適用于在任意非線性非高斯動態(tài)系統(tǒng)，但濾波器迭代次數(shù)的增加會導致粒子退化現(xiàn)象。本文通過殘差重采樣方法在一定程度上解決了粒子退化問題，但隨著重采樣次數(shù)的增加，采樣后的粒子集會逐漸降低其多樣性。因此，既要保持粒子多樣性，又要解決粒子退化問題算法是研究的重點。此外，本研究將粒子濾波同化算法應(yīng)用于作物單產(chǎn)估測，在今后的研究中，可以將數(shù)據(jù)同化算法進一步擴展為作物長勢監(jiān)測等方面的研究。

本文研究LAI和VTCI雙參數(shù)在夏玉米主要生育期與產(chǎn)量的關(guān)系，并沒有將研究區(qū)的研究樣點進行分類，即旱作樣點和灌溉樣點，未來研究可以基于不同類型樣點和不同類型區(qū)域同化最優(yōu)變量獲取的估產(chǎn)模型進行研究區(qū)的單產(chǎn)估測，以提高夏玉米單產(chǎn)估測的精度。

本文選用的MODIS數(shù)據(jù)空間分辨率較低，今后可將高分辨率影像作為作物估產(chǎn)的數(shù)據(jù)源，且隨著農(nóng)業(yè)監(jiān)測手段越來越豐富，不同時間、不同類型、不同維度的傳感器為農(nóng)業(yè)估產(chǎn)提供數(shù)據(jù)，最終形成農(nóng)業(yè)大數(shù)據(jù)，因此利用大數(shù)據(jù)建立作物估產(chǎn)模型滿足不同類型、不同尺度的估產(chǎn)是需要研究的內(nèi)容。

5 結(jié)論

(1)基于RF算法對2013—2018年河北省中部平原夏玉米8個研究樣點進行單點同化。同化LAI和VTCI能較好地綜合模型模擬值和遙感觀測值的優(yōu)點，同化LAI可以減緩CERES-Maize模型模擬LAI的劇烈變化，改善了MODIS觀測值因受空間分辨率和云影響而導致的偏低現(xiàn)象。同化VTCI與旬累積降水量數(shù)據(jù)對比分析表明，同化VTCI的變化趨勢明顯改善，其效果優(yōu)于遙感反演的VTCI，更適合區(qū)域旱情的監(jiān)測。

(2)基于同化雙參數(shù)構(gòu)建的估產(chǎn)模型的決定系數(shù)R2最大為0.615 6，且達極顯著水平(P<0.001)，因此作為較優(yōu)估產(chǎn)模型。利用2015年53個縣(區(qū))夏玉米統(tǒng)計單產(chǎn)對較優(yōu)估產(chǎn)模型進行了精度驗證，結(jié)果表明，所有縣(區(qū))估測單產(chǎn)的平均相對誤差為8.43%，RMSE為631 kg/hm2，NRMSE為10.50%。用較優(yōu)單產(chǎn)估測模型對2013—2018年各縣(區(qū))的夏玉米產(chǎn)量進行估測，得到河北省中部平原玉米區(qū)域單產(chǎn)估測結(jié)果。結(jié)果表明，區(qū)域估測結(jié)果的時空分布與河北省中部平原夏玉米的實際情況相符，應(yīng)用RF算法同化LAI和VTCI的單產(chǎn)估測精度較高，適用于區(qū)域夏玉米單產(chǎn)的估算。