SCE標定結(jié)合EnKF同化遙感和WOFOST模型模擬冬小麥時序LAI

許 偉，秦其明，張?zhí)碓?，龍澤?/p>

SCE標定結(jié)合EnKF同化遙感和WOFOST模型模擬冬小麥時序LAI

許 偉，秦其明※，張?zhí)碓?，龍澤?/p>

（北京大學遙感與地理信息系統(tǒng)研究所，北京 100871）

WOFOST（world food studies）模型可用于模擬冬小麥全生育期內(nèi)的時序葉面積指數(shù)（leaf area index, LAI），各器官生物量以及最終產(chǎn)量，對冬小麥的長勢監(jiān)測與產(chǎn)量預估有著重要意義。但將WOFOST模型用于中國具體區(qū)域的冬小麥生長模擬時，存在著參數(shù)定標困難、模擬結(jié)果不夠準確等嚴重問題。目前對該模型的定標大多依靠研究者的經(jīng)驗進行，雖已總結(jié)出了一套從標定到模擬應用的研究方法，但在區(qū)域模擬時仍然存在很多問題。為此，該文以較易獲取的LAI為參考指標，結(jié)合潛在生長水平模式下的WOFOST模型在衡水地區(qū)的應用，提出了一種“區(qū)域優(yōu)化標定，像元同化修正”的研究方法：首先在區(qū)域尺度上對WOFOST模型進行優(yōu)化標定，利用擴展傅里葉幅度靈敏度檢驗法（extend fourier amplitude sensitivity test, EFAST）分析模型各個參數(shù)的敏感性，在此基礎上選擇了可以迅速找到全局最優(yōu)解的SCE（shuffled complex evolution）算法對總敏感度最高的5個參數(shù)進行優(yōu)化，并將優(yōu)化前后的時序LAI曲線進行對比；其次運用第一步確定的模型最優(yōu)參數(shù)，在對區(qū)域內(nèi)每個像元進行模擬時，結(jié)合Sentinel-2衛(wèi)星數(shù)據(jù)反演所得的各個像元LAI，利用集合卡爾曼濾波（ensemble kalman filter, EnKF）在像元尺度上對LAI進行同化修正，并結(jié)合采樣點的2次實測LAI數(shù)據(jù)對同化所得結(jié)果進行驗證。試驗發(fā)現(xiàn)，優(yōu)化標定后的WOFOST模型模擬所得LAI曲線更接近所給的LAI真值，在此基礎上結(jié)合數(shù)據(jù)同化模擬得出的衡水地區(qū)每個像元LAI的2達到0.87，RMSE僅為0.62。因此，與原來只能通過經(jīng)驗進行定標的方法相比，該方法有效地解決了WOFOST模型在具體應用中亟待解決的復雜標定問題，并且結(jié)合同化修正有效地提高了模型在各個像元的模擬精度，2由0.70～0.83提升至了0.87，RMSE由0.89～1.36降低至了0.62。同時該文也提供了從模型標定到具體模擬整個過程中各個環(huán)節(jié)的思路與方法，有利于促進WOFOST模型在區(qū)域尺度上的應用。

模型；遙感；冬小麥；時序LAI；WOFOST；SCE

0 引 言

小麥作為中國三大糧食作物之一，在農(nóng)業(yè)生產(chǎn)中有著不可替代的地位，而葉面積指數(shù)（leaf area index，LAI）作為小麥生長過程中的一個重要參數(shù)，體現(xiàn)了小麥光合、呼吸和蒸騰作用等生物物理過程的能力，與小麥的生長動態(tài)以及生態(tài)系統(tǒng)生產(chǎn)力密切相關[1-2]。因此準確高效地獲取小麥LAI及其動態(tài)變化信息對小麥的長勢監(jiān)測與產(chǎn)量預估具有重要意義[3]。

獲取小麥LAI的方法一般有4種。第一種是利用儀器進行實地測量，如由美國LI-COR公司生產(chǎn)的LAI2000型植物冠層分析儀。實地測量的LAI數(shù)據(jù)準確度高，但是只能定點定時測量，難以獲取一定空間尺度和長時間序列的LAI數(shù)據(jù)。第二種方式是利用遙感數(shù)據(jù)，通過建立小麥光譜反射率與LAI的定量關系，反演得到小麥LAI，一般分為經(jīng)驗模型和物理模型。經(jīng)驗模型的本質(zhì)是通過建立光譜反射率與LAI的回歸關系來進行LAI的反演[4]；物理模型則以輻射傳輸模型為基礎，通過分析電磁輻射在植被冠層內(nèi)部的傳輸與交互作用，找到植物參量與冠層反射率之間的物理關系[5]，目前最常用的是PROSAIL模型[6]。利用遙感反演的手段可以快速高效地獲取大面積的小麥LAI，但是由于衛(wèi)星重返周期以及天氣的限制，難以獲取小麥全生育期的時序LAI。第三種方式是利用作物生長模型結(jié)合研究區(qū)的天氣數(shù)據(jù)，模擬小麥全生育期的生長狀況，進而可以得到小麥全生育期的LAI，常用的有WOFOST（world food studies）模型[7]。但一般模型參數(shù)眾多，在特定區(qū)域使用前必須進行本地化。第四種方式則是在作物生長模型模擬的基礎上，結(jié)合地面實測數(shù)據(jù)或者衛(wèi)星遙感數(shù)據(jù)進行數(shù)據(jù)同化。同化常用的方法有優(yōu)化法[如四維變分法（four-dimensional variational data assimilation, 4DVar）[8]]、濾波法[如卡爾曼濾波（kalman filter, KF）[9]]等，通過最小化觀測值與模擬值之間的差異，修正所得的時序LAI曲線，使得模擬結(jié)果更加準確。

本文選擇潛在生長水平模式的WOFOST模型結(jié)合氣象站提供的每日最高與最低氣溫數(shù)據(jù)和日照數(shù)據(jù)模擬衡水地區(qū)冬小麥的生長，以獲得小麥全生育期的時序LAI。如前文所述，本研究開展的關鍵是對模型的參數(shù)進行標定使其適用于衡水地區(qū)冬小麥生長的模擬，針對這方面的研究，國內(nèi)外學者已經(jīng)開展了很多工作。張素青等[10]在河南省夏玉米主產(chǎn)區(qū)對WOFOST模型進行了校準和驗證。黃健熙等[11]根據(jù)經(jīng)驗將模型參數(shù)劃分為敏感參數(shù)與非敏感參數(shù)，敏感參數(shù)通過農(nóng)氣站點記錄的作物與環(huán)境參數(shù)計算取值范圍，再通過優(yōu)化法確定，非敏感參數(shù)則采用默認值或查閱相關文獻得到。何亮等[12]利用擴展傅里葉幅度靈敏度檢驗法（extend fourier amplitude sensitivity test, EFAST）[13]分析了WOFOST模型各個參數(shù)的敏感性，再利用馬爾科夫蒙特卡洛方法（markov chain monte carlo，MCMC）對這些敏感參數(shù)進行優(yōu)化。可以看出，現(xiàn)有工作中對WOFOST模型參數(shù)標定大多數(shù)根據(jù)經(jīng)驗進行的，在優(yōu)化敏感參數(shù)的過程中也未采取科學高效的優(yōu)化算法，主觀性較大。何亮等采用的EFAST可以較好地得到模型各個參數(shù)的敏感性，但MCMC優(yōu)化算法依賴于初值的選取，難以收斂。

一般來說，當研究區(qū)域內(nèi)的農(nóng)作物屬于同一品種時，WOFOST模型的參數(shù)是唯一確定的，在此情況下每個像元若采用不同的參數(shù)反而是不合理的，也即研究區(qū)域內(nèi)的每個像元共用同一套標定后的模型參數(shù)。但在每個像元內(nèi)的作物生長可能由于土壤、水分、播種期等條件的不同而有一定的差異，僅利用同一套參數(shù)的潛在水平下的WOFOST模型結(jié)合氣象數(shù)據(jù)對作物的生長進行模擬在不同的像元可能存在著一定的偏差。針對此問題，范麗穎[14]在利用經(jīng)驗知識將WOFOST模型進行標定的基礎上，利用了集合卡爾曼濾波（ensemble kalman filter, EnKF）[15]方法對LAI的模擬進行同化，以此來提高LAI的模擬精度。但其WOFOST模型的標定僅僅是依靠人為的經(jīng)驗判定，可能會引入較大的誤差。

因此，綜合以上WOFOST模型在實際應用中的問題，本文提出一種“區(qū)域優(yōu)化標定，像元同化修正”的研究方法。首先，針對WOFOST模型過程中過于依賴人為經(jīng)驗的問題，本文先通過EFAST全局敏感性分析得到模型敏感性最大的5個參數(shù)，再利用了SCE（ShuffledComplexEvolution）[16]優(yōu)化算法結(jié)合研究區(qū)域單個站點的觀測數(shù)據(jù)對這5個參數(shù)在相關研究者們探索出的大致取值范圍內(nèi)進行了優(yōu)化，從而實現(xiàn)模型的標定，標定后的參數(shù)則會被應用于整個研究區(qū)域進行作物生長的模擬。其次，為了進一步提高模型模擬的精度，在模型標定的基礎上，本文結(jié)合Sentinel-2[17]衛(wèi)星數(shù)據(jù)反演得到幾個特定日期每個像元的LAI值，運用EnKF對模型在每個像元的模擬進行數(shù)據(jù)同化，以修正由于不同像元的環(huán)境差異對作物生長帶來的不同影響，由此得到衡水地區(qū)冬小麥全生育期的LAI。

1 研究區(qū)與數(shù)據(jù)來源

1.1 研究區(qū)概況

本文以河北省衡水市作為試驗區(qū)。衡水市位于河北省東南部（115°0′－116°34′E，37°03′－38°23′N），占地面積約8 815 km2。衡水市平均海拔約為21 m，西南部地勢較高，東北部地勢較低，屬于大陸季風氣候區(qū)。衡水市共包括2個市轄區(qū)（桃城區(qū)、冀州區(qū)），1個縣級市（深州市）與8個縣（棗強縣、武邑縣、武強縣、饒陽縣、安平縣、故城縣、景縣、阜城縣）。

衡水市地處河北沖積平原，是中國華北平原冬小麥的主產(chǎn)區(qū)之一。根據(jù)農(nóng)氣站物候資料，衡水地區(qū)冬小麥的播種時間一般約為10月中上旬，收獲時間一般為次年6月初。衡水市一般采用灌溉式方式種植小麥，水肥條件充足，受水分脅迫或養(yǎng)分脅迫等因素的影響較小。

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

1.2.1 地面實測數(shù)據(jù)

LAI地面實測數(shù)據(jù)利用LAI2000儀器測量得到，實測數(shù)據(jù)采樣點的分布如圖1所示。采樣點共有22個，每個縣各分布2個，編號如圖1所示。在每個采樣點附近隨機選取5個采樣位置并記錄經(jīng)緯度，每個位置進行5次LAI的測量，取平均值作為對應位置的LAI真值。

圖1 研究區(qū)域與采樣點分布

Fig.1 Studyareaandsamplingpoints distribution

在小麥的整個生育期內(nèi)，共進行了2次LAI的實地測量，測量日期分別為2017-03-29—2017-04-01以及2017-05-04—2017-05-06，分別位于冬小麥生長發(fā)育的拔節(jié)期與抽穗—灌漿期。這兩個生育期在冬小麥的生長過程中具有重要的意義，并且兩個時期的生長狀況不同，LAI相差較大，因此利用它們來對LAI的模擬進行同化修正可以有效的提升模擬精度。

兩次地面測量共得到了110組實測的LAI值。主要用于：建立Sentinel-2衛(wèi)星數(shù)據(jù)的LAI經(jīng)驗反演模型；SCE優(yōu)化算法中的LAI真值輸入；驗證同化之后WOFOST模型輸出的LAI精度。

1.2.2 氣象站數(shù)據(jù)

為進行WOFOST模型的標定，本文選擇位于衡水市地區(qū)的國家氣象站——南宮氣象站（臺站編號54705，位于115.38°E，37.22°N）提供的2016年9月－2017年7月的中國地面氣候資料日值數(shù)據(jù)（http://data.cma.cn）作為氣象數(shù)據(jù)驅(qū)動WOFOST模型。在潛在生長水平模式下，WOFOST模型需要輸入的驅(qū)動數(shù)據(jù)有日最高氣溫、日最低氣溫與日照時數(shù)，其余的如降水、肥力等數(shù)據(jù)則采用潛在水平下模型的默認數(shù)據(jù)。

為得到整個衡水地區(qū)每個像元內(nèi)的LAI模擬曲線，本文選取了21個河北省內(nèi)的國家氣象臺站（張北、蔚縣、石家莊、邢臺、豐寧、圍場、張家口、懷來、承德、遵化、青龍、秦皇島、霸州、唐山、唐海、樂亭、保定、饒陽、泊頭、黃驊、南宮）的氣象數(shù)據(jù)進行插值，并提取出衡水市內(nèi)每個像元上的天氣數(shù)據(jù)作為WOFOST模型的驅(qū)動數(shù)據(jù)。

1.2.3 Sentinel-2數(shù)據(jù)

Sentinel-2衛(wèi)星數(shù)據(jù)在可見光、近紅外、短波紅外波長范圍內(nèi)共提供13個波段。本文選擇空間分辨率為10 m的紅光波段與近紅外波段反射率數(shù)據(jù)，與兩次地面測量相對應，所選擇的Sentinel-2遙感影像的過境日期為2017-03-29與2017-05-04。這2景數(shù)據(jù)與地面實測數(shù)據(jù)一起構建LAI反演的經(jīng)驗模型，從而反演出衡水地區(qū)的LAI，為后續(xù)像元尺度的同化修正提供真值。

2 試驗方法

本文的試驗流程可概括如下：首先利用單點的氣象數(shù)據(jù)驅(qū)動潛在水平下的WOFOST模型，結(jié)合EFAST算法得到影響模型運行最敏感的5個參數(shù)；再利用對應的單點實測LAI數(shù)據(jù)結(jié)合SCE算法對這幾個敏感參數(shù)進行優(yōu)化標定；接著將標定后的參數(shù)運用于整個衡水地區(qū)的每個像元，利用插值得到的每個像元的氣象數(shù)據(jù)作為驅(qū)動，并將Sentinel-2反演得到的LAI數(shù)據(jù)作為真值對模型的模擬進行同化修正，最終得到每個像元的生長模擬結(jié)果，技術路線圖如圖2所示。

2.1 基于Sentinel-2的LAI反演

傳統(tǒng)的基于NDVI模型反演的LAI在濃密植被覆蓋地區(qū)往往出現(xiàn)飽和情況，為解決這一問題，Sun等[18]提出了反轉(zhuǎn)差值植被指數(shù)（inverted difference vegetation index，IDVI）用于高植被覆蓋地區(qū)的LAI反演。

式中red和nir分別為紅光和近紅外波段的地表反射率。IDVI在濃密植被地區(qū)與LAI有較好的線性關系，而在植被稀疏地區(qū)的相關性較差。

為準確反演出研究區(qū)LAI，本文基于NDVI與IDVI構建統(tǒng)計回歸模型，以解決2種指數(shù)各自在高植被覆蓋區(qū)和低植被覆蓋區(qū)的不足。

式中是Logistic形式的動態(tài)比例因子，NDVI為NDVI達到飽和時的閾值，為Logistic方程中的幅度系數(shù)。NDVI和在研究區(qū)域應用之前應當根據(jù)先驗知識進行更新和調(diào)整。在本文中的具體表達式為

本文所使用的數(shù)據(jù)以及LAI的反演過程與Sun等[18]一致，通過文獻[18]可知，利用IDVI反演得到的LAI的為0.91，RMSE為0.81.

圖2 技術路線圖

2.2 EFAST

EFAST是一種基于方差理論的定量全局敏感性分析方法，對樣本要求低，計算高效[13]，適用于WOFOST這種參數(shù)眾多且關系復雜的模型的敏感性分析。其算法的基本思想是分解模型參數(shù)對模型模擬結(jié)果的方差，把參數(shù)敏感性分為2種類型：各參數(shù)對結(jié)果的影響以及參數(shù)之間耦合對模型結(jié)果的影響。其中單個參數(shù)獨立作用是用主敏感度（main effect）衡量，而參數(shù)相互作用（interaction）則用總敏感度（total effect）和主敏感度之差來衡量。

模型的總方差分解為

式中為V參數(shù)x輸入變化單獨引起的模型方差，V為參數(shù)x通過參數(shù)x作用貢獻的耦合方差，V為參數(shù)x通過參數(shù)x、x作用貢獻的耦合方差，依此類推。

通過歸一化處理參數(shù)x的主敏感度S定義為

參數(shù)x的總敏感度為

式中?i為與參數(shù)x無關的所有其他參數(shù)方差之和（即去掉式（11）中所有含下標的方差項）。

2.3 SCE

SCE算法是段青云博士于1992年在Arizona大學水文與水資源系作研究時提出的全局優(yōu)化算法，最初應用于水文模型的優(yōu)化中[16]。它可以解決MCMC算法依賴于初值而難以收斂的問題，并且可以避免優(yōu)化過程中陷入局部極值區(qū)而無法收斂于全局最優(yōu)值。該算法繼承了下山單純形（downhill simplex）方法、控制隨機搜索（controlled random search）方法，競爭演化（competitive evolution）方法和洗牌（shuffling）算法等算法的特色，通過最小化代價函數(shù)，迅速找到全局最優(yōu)解。

3）樣本點排序：把個樣本點函數(shù)值（x,f）按升序排列，排序后仍記為(x,f)，=1,…,。其中,將排序后的點存在中；

4）劃分復合形群體：將劃分為個復合形1,···,A，每個復合形含有個點。

5）復合形演化：通過CCE方法[16]，對每個復合形分別演化。

6）復合形混合：把演化后的每個復合形的所有點組合，生成新的點集，再次按函數(shù)值f升序排列；

7）收斂判斷：如果滿足收斂條件則停止，否則返回步驟（4）。

其中的CCE算法是下山單純形方法的拓展，通過競爭機制決定復雜形參與演化的是哪些點。該算法根據(jù)每個點的目標函數(shù)值，采用特定方法來確定每個點被選擇的概率，從而保證目標函數(shù)值較高的點，參加演化的概率大于目標函數(shù)值較低的點，同時目標函數(shù)最低的點也有參加演化的機會。

2.4 EnKF

EnKF由Evensen提出，其核心思想是利用集合成員的統(tǒng)計來得到誤差協(xié)方陣，以避免復雜算子的計算[15]。EnKF在處理各種數(shù)據(jù)不確定性時比較靈活，算法易于實現(xiàn)和操作，因此被廣泛應用于地表數(shù)據(jù)同化研究中。

同化的過程分為預測和分析2個部分。

首先通過模型的運行得到LAI的一步預測值：

但一步預測值只是通過模型向前運行得到的，要得到時刻LAI的估計值，還需要將該時刻的衛(wèi)星反演值同化進去

根據(jù)以上處理，可以得到

由此可以方便地計算出濾波增益的值，進而得到時刻LAI的估計值。同理可以得到時刻LAI估計值的誤差P。

3 結(jié)果與分析

3.1 敏感性分析結(jié)果

WOFOST模型可輸出作物全生育期逐日的葉面積指數(shù)、地上總干物質(zhì)量以及最終產(chǎn)量，對于不同的目標輸出，模型輸入?yún)?shù)的影響程度也不同（例如可以模擬出最優(yōu)LAI的模型不一定能夠模擬出最優(yōu)的最終產(chǎn)量）。由于本文同化的參數(shù)為LAI，因此選取作物生育期最大葉面積指數(shù)（LAImax）作為分析對象，結(jié)合南宮氣象站提供的氣象數(shù)據(jù)作為驅(qū)動，在單點采用EFAST算法對WOFOST模型的14項重要輸入?yún)?shù)進行敏感性分析。參數(shù)的取值范圍與最終分析結(jié)果如表1所示。

表1 WOFOST模型參數(shù)對冬小麥生育期最大葉面積指數(shù)（LAImax）的影響

從表1中可以看出，對作物生育期最大葉面積指數(shù)（LAImax）總敏感度最大的5個模型輸入?yún)?shù)依次為：SPAN、TBASE、CVL、CVS、RGRLAI。其中SPAN、TBASE 和RGRLAI與作物葉片的存活時間相關，CVS、CVL則與呼吸作用積累同化物以及同化物轉(zhuǎn)化成干物質(zhì)密切相關。

3.2 SCE優(yōu)化算法結(jié)果

對于上述的5個總敏感度最高的參數(shù)，進一步在該點采用SCE算法對模型進行優(yōu)化標定，對于其他參數(shù)，采用相關文獻中的推薦值[11-12,19-28]。最終WOFOST模型輸入?yún)?shù)如表2所示。

3.3 冬小麥全生育期LAI反演結(jié)果與對比驗證

3.3.1 SCE優(yōu)化前后LAI結(jié)果對比

為展示SCE優(yōu)化方法的效果，分別利用SCE算法優(yōu)化前后的WOFOST模型模擬采樣點HS01所在區(qū)域內(nèi)的冬小麥生長狀況，如圖3所示。

從圖3中可以看出，優(yōu)化后的時序LAI曲線更加接近所給的2個日期的LAI真值，也即通過SCE優(yōu)化后的WOFOST模型更適合對衡水地區(qū)冬小麥的生長狀況進行模擬。

表2 SCE算法優(yōu)化WOFOST模型標定后部分參數(shù)取值

注：“/”左側(cè)為優(yōu)化前的參數(shù)值，“/”右側(cè)為優(yōu)化后的參數(shù)值。

Note: The values before optimization are on the left hand of ‘/’ and that after optimization are on the right hand of ‘/’.

圖3 HS01采樣點在WOFOST模型經(jīng)SCE優(yōu)化前后的冬小麥時序LAI曲線

3.3.2 EnKF同化結(jié)果與驗證

利用標定后的WOFOST模型模擬衡水地區(qū)每個像元內(nèi)的冬小麥生長情況，并結(jié)合Sentinel-2反演所得到的LAI值進行數(shù)據(jù)同化，并分別從時間和空間上對所得到的結(jié)果進行展示。

時間上選取采樣點HS01所得到的時序LAI曲線進行繪制，如圖4所示；空間上選取2017-03-30與2017-05-05這2個日期整個衡水市的LAI分布情況進行繪制，如圖5所示。

圖4 WOFOST模型經(jīng)EnKF同化后在HS01采樣點的冬小麥全生育期時序LAI曲線

圖5 SCE定標與EnKF同化后模型模擬的特定日期的衡水地區(qū)LAI分布

可以發(fā)現(xiàn)，在圓點處，原本光滑的LAI曲線會出現(xiàn)向點靠近的凹陷與凸起，這種現(xiàn)象說明了WOFOST模型的模擬與真實值有一定的偏差，而數(shù)據(jù)同化可以通過在某些時刻加入觀測值，將模擬曲線拉回到與真實值更加接近的水平，進而提高模型的模擬精度。

圖5中，模擬結(jié)果在衡水市的北部有一片空白區(qū)域，這是因為所選擇的兩景Sentinel-2影像沒有覆蓋到這些地區(qū)，因此依賴于觀測值的數(shù)據(jù)同化無法進行。為了直觀的展示出小麥種植區(qū)的LAI分布，在遙感影像覆蓋區(qū)域，本文利用2017衡水地區(qū)冬小麥的種植面積[29]將模擬結(jié)果進行掩膜提取。

為驗證本文研究方法所模擬得到的衡水地區(qū)LAI的精度，將圖5中2個日期的LAI模擬值分布圖中與采樣點相對應像元的值分別提取出來，并與地面實測數(shù)據(jù)進行對比驗證，如圖6所示。

圖6 WOFOST模擬所得LAI與實測數(shù)據(jù)比較

從圖6中可以看出，同化得到的LAI與實測LAI的散點基本分布在1:1線附近，當LAI值較低時，同化得到的LAI存在高估現(xiàn)象，當LAI值較高時，同化得到的LAI存在低估現(xiàn)象。同化結(jié)果與實測數(shù)據(jù)存在較好的相關性，2為0.87，RMSE為0.62，同化結(jié)果與實測數(shù)據(jù)之間的差異較小。而只采取查閱文獻與經(jīng)驗標定的方法對參數(shù)進行校準的WOFOST模型模擬得到的LAI與真實值之間的2在0.70～0.83之間，RMSE在0.89～1.36之間[11]。因此，利用本文提出的方法所模擬的LAI的2有所提高，RMSE有所降低，同化效果較好，所得LAI的精度較高，可以有效表征冬小麥的長勢情況。

4 討 論

本文的目的是使用WOFOST模型模擬衡水地區(qū)冬小麥的生長狀況，從而得到全生育期的時序LAI曲線。在模型的標定上，本文利用EFAST模型獲取最敏感的幾個參數(shù)，既避免了憑經(jīng)驗選取參數(shù)的主觀性，又減小了后續(xù)優(yōu)化過程的計算量；同時，本文利用的SCE算法可以使待優(yōu)化的參數(shù)快速收斂到全局最優(yōu)值，解決了MCMC算法依賴于初值選取的問題。在模型的區(qū)域尺度應用方面，本文通過EnKF結(jié)合遙感數(shù)據(jù)反演LAI對標定后的模型進行數(shù)據(jù)同化，使得模型在每個像元上模擬得到的LAI準確度更高?？偟膩碚f，本文所提出的“區(qū)域優(yōu)化標定，像元同化修正”的研究方法對WOFOST模型在中國區(qū)域尺度上的應用具有一定的參考意義。

在結(jié)果驗證方面，采取時序的LAI實測數(shù)據(jù)或已有的時序LAI產(chǎn)品與模型模擬所得到的LAI曲線進行對比更能反應本文研究工作的意義。但時序LAI的實測耗時耗力，現(xiàn)有的MODIS-LAI產(chǎn)品分辨率過大，且精度又達不到應用需求[26]，因此本文只選擇了小麥生長過程中兩個時間點的實測LAI數(shù)據(jù)在空間尺度上對WOFOST模型的模擬結(jié)果進行了驗證。從驗證的結(jié)果可以看出，模型在某些采樣點的模擬值并不夠準確，但WOFOST模型的目的在于獲取小麥生長過程中LAI變化趨勢，其單點模擬LAI的準確度并不能達到遙感反演結(jié)果的水平，與其他WOFOST模型的應用方法比較，本文的研究方法所模擬得到的LAI精度已經(jīng)達到了較高的水準。同時，這些不準確的點也反映了在區(qū)域內(nèi)每個像元上僅僅使用標定的WOFOST模型有很大偏差，在定標之后進行“像元同化修正”是必要的。

5 結(jié) 論

本文通過潛在生長水平下的WOFOST模型對衡水地區(qū)冬小麥的LAI進行模擬分析，提出了一種“區(qū)域優(yōu)化標定，像元同化修正”的研究方法。結(jié)果表明，利用EFAST算法對模型的參數(shù)敏感性進行分析，再結(jié)合地面實測數(shù)據(jù)，利用SCE算法對最敏感的幾個參數(shù)進行優(yōu)化后的模型更適用于該地區(qū)。為準確地模擬區(qū)域尺度上的冬小麥時序LAI，本文結(jié)合了Sentinel-2衛(wèi)星數(shù)據(jù)反演得到的LAI，利用EnKF算法對優(yōu)化后的WOFOST模型進行數(shù)據(jù)同化，所得到的LAI的2為0.87，RMSE為0.62，均達到較好的水平，模擬精度較高。本文的意義在于提供了從模型標定到具體模擬整個過程中各個環(huán)節(jié)的思路與方法，使得對WOFOST模型的應用不再僅僅依靠經(jīng)驗，這有利于WOFOST模型在中國區(qū)域尺度上的應用。

[1] Fan W J, Gai Y Y, Xu X R, et al. The spatial scaling effect of the discrete-canopy effective leaf area index retrieved by remote sensing[J]. Science China Earth Sciences, 2013, 56(9): 1548－1554.

[2] Liu Yang, Liu Ronggao, Chen Jingming, et al. Current status and perspectives of leaf area index retrieval from optical remote sensing data[J]. Journal of Geo-Information Science, 2013, 15(5):734－743.

[3] 劉良云．高光譜遙感在精準農(nóng)業(yè)中的應用研究[R]．中國科學院遙感應用研究所博士后出站報告，2002．

[4] Tillack A, Clasen A, Kleinschmit B, et al. Estimation of the seasonal leaf area index in an alluvial forest using high-resolution satellite-based vegetation indices[J]. Remote Sensing of Environment, 2014, 141(141): 52－63.

[5] Rasmus Houborg, Henrik Soegaard. Combining vegetation index and model inversion methods for the extraction of key vegetation biophysical parameters using Terra and Aqua MODIS reflectance data[J]. Remote Sensing of Environment, 2007, 106(1): 39－58.

[6] Baret F, Jacquemoud S, Guyot G, et al. Modeled analysis of the biophysical nature of spectral shifts and comparison with information content of broad bands[J]. Remote Sensing of Environment, 1992, 41(2/3): 133－142.

[7] Diepen C A, Wolf J, Keulen H, et al. WOFOST: A simulation model of crop production[J]. Soil Use & Management, 2010, 5(1): 16－24.

[8] Rabier F, Thépaut J N, Courtier P. Extended assimilation and forecast experiments with a four‐dimensional variational assimilation system[J]. Quarterly Journal of the Royal Meteorological Society, 1998, 124(550): 1861－1887.

[9] Kalman R E. A new approach to linear filtering and prediction problems[J]. Journal of Basic Engineering Transactions, 1960, 82: 35－45.

[10] 張素青，張建濤，李繼蕊，等. WOFOST模型在河南省夏玉米主產(chǎn)區(qū)的校準與驗證[J]. 河南農(nóng)業(yè)科學，2014，43(8)：152－156.

Zhang Suqing, Zhang Jiantao, Li Jirui, et al. Calibration and validation of WOFOST in main maize-producing regions in Henan[J]. Journal of Henan Agricultural Sciences, 2014, 43(8): 152－156. (in Chinese with English abstract)

[11] 黃健熙，賈世靈，馬鴻元，等. 基于WOFOST模型的中國主產(chǎn)區(qū)冬小麥生長過程動態(tài)模擬[J]. 農(nóng)業(yè)工程學報，2017，33(10)：222－228.

Huang Jianxi, Jia Shiling, Ma Hongyuan, et al. Dynamic simulation of growth process of winter wheat in main production areas of China based on WOFOST model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(10): 222－228. (in Chinese with English abstract)

[12] 何亮，侯英雨，趙剛，等. 基于全局敏感性分析和貝葉斯方法的WOFOST作物模型參數(shù)優(yōu)化[J]. 農(nóng)業(yè)工程學報，2016，32(2)：169－179.

He Liang, Hou Yingyu, Zhao Gang, et al. Parameters optimization of WOFOST model by integration of global sensitivity analysis and Bayesian calibration method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(2): 169－179. (in Chinese with English abstract)

[13] Mcculloch A. Sensitivity analysis in practice: A guide to assessing scientific models[J]. Publications of the American Statistical Association, 2005, 101(473): 398－399.

[14] 范麗穎. 基于集合Kalman濾波的作物生長模型與遙感數(shù)據(jù)的同化研究[D]. 北京：中國農(nóng)業(yè)科學院，2012.

[15] Evensen G. The ensemble Kalman filter: Theoretical formulation and practical implementation[J]. Ocean Dynamics, 2003, 53(4): 343－367.

[16] Duan Q Y, Gupta V K, Sorooshian S. Shuffled complex evolution approach for effective and efficient global minimization[J]. Journal of Optimization Theory and Applications, 1993, 76(3): 501－521.

[17] Richter K, Atzberger C, Vuolo F, et al. Evaluation of Sentinel-2 spectral sampling for radiative transfer model based LAI estimation of wheat, sugar beet, and maize[J]. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 2011, 4(2): 458－464.

[18] Sun Y, Ren H, Zhang T, et al. Crop leaf area index retrieval based on inverted difference vegetation index and NDVI[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15(11): 1662－1666.

[19] Huang J, Sedano F, Huang Y, et al. Assimilating a synthetic Kalman filter leaf area index series into the WOFOST model to improve regional winter wheat yield estimation[J]. Agricultural and Forest Meteorology, 2016, 216: 188－202.

[20] Huang J, Tian L, Liang S, et al. Improving winter wheat yield estimation by assimilation of the leaf area index from Landsat TM and MODIS data into the WOFOST model[J]. Agricultural and Forest Meteorology, 2015, 204: 106－121.

[21] Ma G, Huang J, Wu W, et al. Assimilation of MODIS-LAI into the WOFOST model for forecasting regional winter wheat yield[J]. Mathematical & Computer Modelling, 2013, 58(3/4): 634－643.

[22] Ma Hongyuan, Huang Jianxi, Zhu Dehai, et al. Estimating regional winter wheat yield by assimilation of time series of HJ-1 CCD NDVI into WOFOST-ACRM model with Ensemble Kalman filter[J]. Mathematical & Computer Modelling, 2013, 58(3/4): 759－770.

[23] Huang Jianxi, Ma Hongyuan, Su Wei, et al. Jointly assimilating MODIS LAI and ET products into the SWAP model to estimate winter wheat yield[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(8): 4060－4071.

[24] Tian Liyan, Li Zhongxia, Huang Jianxi, et al. Comparison of two optimization algorithms for estimating regional winter wheat yield by integrating MODIS leaf area index and world food studies model[J]. Sensor Letters, 2013, 11(6/7): 1261－1268.

[25] 黃健熙，馬鴻元，田麗燕，等. 基于時間序列LAI和ET同化的冬小麥遙感估產(chǎn)方法比較[J]. 農(nóng)業(yè)工程學報，2015，31(4)：197－203.

Huang Jianxi, Ma Hongyuan, Tian Liyan, et al. Comparison of remote sensing yield estimation methods for winter wheat based on assimilating time-sequence LAI and ET[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(4): 197－203.(in Chinese with English abstract)

[26] 黃健熙，武思杰，劉興權，等. 基于遙感信息與作物模型集合卡爾曼濾波同化的區(qū)域冬小麥產(chǎn)量預測[J]. 農(nóng)業(yè)工程學報，2012，28(4)：142－148.

Huang Jianxi, Wu Sijie, Liu Xingquan, et al. Regional winter wheat yield forecasting based on assimilation of remote sensing data and crop growth model with Ensemble Kalman method[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(4): 142－148. (in Chinese with English abstract)

[27] 邱美娟，宋迎波，王建林，等. 山東省冬小麥產(chǎn)量動態(tài)集成預報方法[J]. 應用氣象學報，2016，27(2)：191－200.

Qiu Meijuan, Song Yingbo, Wang Jianlin, et al. Integrated technology of yield dynamic prediction of winter wheat in Shandong province[J]. Journal of applied meteorological science, 2016, 27(2): 191－200.(in Chinese with English abstract)

[28] Wu D, Yu Q, Lu C, et al. Quantifying production potentials of winter wheat in the North China Plain[J]. European Journal of Agronomy, 2006, 24(3): 226－235.

[29] Sui J, Qin Q, Ren H, et al. Winter wheat production estimation based on environmental stress factors from satellite observations[J]. Remote Sensing, 2018, 10(6): 962.

Time-series LAI simulation of winter wheat based on WOFOST model calibrated by SCE and assimilated by EnKF

Xu Wei, Qin Qiming※, Zhang Tianyuan, Long Zehao

(100871)

WOFOST (world food studies) model can be used to simulate time-series LAI (leaf area index), the organs’ biomass, and the yields of winter wheat. Therefore, it is meaningful for the growth monitoring and production prediction of winter wheat. So far, the calibration of WOFOST usually relies on researches’ experience, which brings many problems while using the model in a specific area. As a result, we focus on the calibration problem and try to improve the accuracy of the simulated results in this paper. The potential production WOFOST was analyzed and LAI was chosen as the measure index because it was easy to obtain. In this study, we selected Hengshui as the study area, and two field experiments were carried out in this area during two different periods. One period was from 2017-03-29 to 2017-04-01 and the other was from 2017-05-04 to 2017-05-06. It was divided into 11 sampling areas and 5 sampling points in every area were obtained to measure the LAI, so we got approximately 110 measured data totally. A method called ‘Calibrating in area by optimization and correcting at pixel by assimilation’ was presented in this paper. Firstly, calibrating WOFOST model in local area: The weather data including sunshine duration data and the maximum and minimum air temperature data every day were used to run the WOFOST model. The data were from Nangong National Weather Station and can be downloaded in National Meteorological Information Center. Then the sensitivity of model parameters can be analyzed with EFAST (extend fourier amplitude sensitivity test) and the 5 most sensitive parameters were selected to optimize the model. It was worthwhile to note that there were different indices to evaluate the sensitivity of every parameter, such as main effect, interaction, and total effect, and the total effect was considered as the most important index in this study. As for the optimization, the SCE (shuffled complex evolution) algorithm was used which could find the global optimal solution fastly. It can solve the initial value dependence problem and local convergence problem which might exist in other optimization algorithms such as MCMC (Markov Chain Monte Carlo). In order to proof that the optimization was valid, the time-series LAI curves simulated were compared by WOFOST before and after optimization with SCE with measured values. It turned out the model after optimization was much more appropriate to simulate the growth of winter wheat in study area. Secondly, assimilating the model in every pixel in the study area: We interpolated weather data from 21 National Weather Stations in Hebei Province in order to run WOFOST in every pixel. Based on this, EnKF (Ensemble Kalman Filter) was used to assimilate LAI in every pixel with the remote sensing data from Sentinel-2. As a result, we could get the time-series LAI curve at every pixel. The LAI curve at point HS01 was illustrated and it was obvious that assimilation made a difference in the simulation. Additionally, the simulated LAI distribution maps were illustrated in Hengshui at date of 2017-03-30 and 2017-05-05. And the simulated LAI values of the pixels according to the sampling points were extracted. By comparing the simulated LAI with measured LAI, we found that2was increased from 0.70-0.83 to 0.87 and RMSE was decreased from 0.89-1.36 to 0.62. Therefore, the method proposed in this study solved the calibration problem and improved the accuracy of time-series LAI simulated compared with other studies. In addition, we provided specific theories and methods in every stage from calibration to application. It contributed to the application of WOFOST in our country.

models; remote sensing; winter wheat; time-series LAI; WOFOST; SCE

2018-08-10

2019-06-12

國家高分重大專項“GF-7衛(wèi)星高精度農(nóng)作物信息提取技術”（11-Y20A16-9001-17/18）；國家重點研發(fā)計劃課題“作物生長與生產(chǎn)力衛(wèi)星遙感監(jiān)測預測”（2016YFD0300603）

許 偉，安徽蕪湖人，博士生，主要研究方向為農(nóng)業(yè)定量遙感。Email：xuwei1995@pku.edu.cn

秦其明，江蘇徐州人，教授，博士，博士生導師，主要研究方向為定量遙感與地理信息系統(tǒng)。Email：qmqin@pku.edu.cn

10.11975/j.issn.1002-6819.2019.14.021

S127; S512.1; TP79

A

1002-6819(2019)-14-0166-08

許 偉，秦其明，張?zhí)碓?，龍澤? SCE標定結(jié)合EnKF同化遙感和WOFOST模型模擬冬小麥時序LAI[J]. 農(nóng)業(yè)工程學報，2019，35(14)：166－173. doi：10.11975/j.issn.1002-6819.2019.14.021 http://www.tcsae.org

Xu Wei, Qin Qiming, Zhang Tianyuan, Long Zehao. Time-series LAI simulation of winter wheat based on WOFOST model calibrated by SCE and assimilated by EnKF[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(14): 166－173. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2019.14.021 http://www.tcsae.org