基于雙參數和粒子濾波同化算法的夏玉米單產估測

王鵬新 胡亞京 李 俐 喬 琛

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

0 引言

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

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

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

1 材料與方法

1.1 研究區概況

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

1.2 數據獲取與預處理

1.2.1CERES-Maize模型

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

(1)

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

as、bs——經驗系數[16]，as取0.18，bs取0.55

w——逐日日照時數

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

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

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

1.2.2旬LAI的生成

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

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

1.2.3旬VTCI的生成

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

(2)

其中

Tmax(NDVIi)=a+bNDVIi

(3)

Tmin(NDVIi)=a′+b′NDVIi

(4)

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

VTCI——條件植被溫度指數

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

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

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

a、b、a′、b′——待定系數，通過研究區域內NDVI和T的散點圖近似獲得

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

1.2.4玉米區域調查統計數據

本研究中各縣的玉米播種面積、總產量數據來源于2013—2016年的《河北農村統計年鑒》，玉米單產數據由總產量和播種面積計算得到。

1.2.5總體研究技術路線

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

2 研究方法

2.1 粒子濾波算法

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

2.1.1預測階段

(5)

(6)

(7)

式中uk——模型驅動參數

(8)

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

2.1.2更新階段

(9)

更新粒子權重

(10)

(11)

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

2.2 隨機森林回歸算法

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

(12)

式中x——研究區域各縣(區)玉米產量影響因子及玉米單產數據

K——決策樹的數量

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

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

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

(13)

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

e′xi——每棵樹對應的袋外數據誤差率

exi——隨機置換變量xi的序列后每棵樹對應的袋外數據誤差率

3 結果與分析

3.1 LAI與VTCI同化結果與分析

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

3.2 同化結果精度驗證

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

3.3 夏玉米估產模型構建

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

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

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

表2 加權LAI和加權VTCI的玉米單產估測模型Tab.2 Yield estimation models of maize between weighted LAI and VTCI

3.4 估產結果驗證

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

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

4 討論

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

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

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

5 結論

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

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