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

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

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

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

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

模型；遙感；冬小麥；時(shí)序LAI；WOFOST；SCE

0 引 言

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

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

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

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

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

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

1.1 研究區(qū)概況

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

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

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

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

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

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

Fig.1 Studyareaandsamplingpoints distribution

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

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

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

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

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

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

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

2 試驗(yàn)方法

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

2.1 基于Sentinel-2的LAI反演

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

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

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

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

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

圖2 技術(shù)路線(xiàn)圖

2.2 EFAST

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

模型的總方差分解為

式中為V參數(shù)x輸入變化單獨(dú)引起的模型方差，V為參數(shù)x通過(guò)參數(shù)x作用貢獻(xiàn)的耦合方差，V為參數(shù)x通過(guò)參數(shù)x、x作用貢獻(xiàn)的耦合方差，依此類(lèi)推。

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

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

式中?i為與參數(shù)x無(wú)關(guān)的所有其他參數(shù)方差之和（即去掉式（11）中所有含下標(biāo)的方差項(xiàng)）。

2.3 SCE

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

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

4）劃分復(fù)合形群體：將劃分為個(gè)復(fù)合形1,···,A，每個(gè)復(fù)合形含有個(gè)點(diǎn)。

5）復(fù)合形演化：通過(guò)CCE方法[16]，對(duì)每個(gè)復(fù)合形分別演化。

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

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

其中的CCE算法是下山單純形方法的拓展，通過(guò)競(jìng)爭(zhēng)機(jī)制決定復(fù)雜形參與演化的是哪些點(diǎn)。該算法根據(jù)每個(gè)點(diǎn)的目標(biāo)函數(shù)值，采用特定方法來(lái)確定每個(gè)點(diǎn)被選擇的概率，從而保證目標(biāo)函數(shù)值較高的點(diǎn)，參加演化的概率大于目標(biāo)函數(shù)值較低的點(diǎn)，同時(shí)目標(biāo)函數(shù)最低的點(diǎn)也有參加演化的機(jī)會(huì)。

2.4 EnKF

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

同化的過(guò)程分為預(yù)測(cè)和分析2個(gè)部分。

首先通過(guò)模型的運(yùn)行得到LAI的一步預(yù)測(cè)值：

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

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

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

3 結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

表2 SCE算法優(yōu)化WOFOST模型標(biāo)定后部分參數(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采樣點(diǎn)在WOFOST模型經(jīng)SCE優(yōu)化前后的冬小麥時(shí)序LAI曲線(xiàn)

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

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

時(shí)間上選取采樣點(diǎn)HS01所得到的時(shí)序LAI曲線(xiàn)進(jìn)行繪制，如圖4所示；空間上選取2017-03-30與2017-05-05這2個(gè)日期整個(gè)衡水市的LAI分布情況進(jìn)行繪制，如圖5所示。

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

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

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

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

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

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

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

4 討 論

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

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

5 結(jié) 論

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

[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] 劉良云．高光譜遙感在精準(zhǔn)農(nóng)業(yè)中的應(yīng)用研究[R]．中國(guó)科學(xué)院遙感應(yīng)用研究所博士后出站報(bào)告，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ū)的校準(zhǔn)與驗(yàn)證[J]. 河南農(nóng)業(yè)科學(xué)，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模型的中國(guó)主產(chǎn)區(qū)冬小麥生長(zhǎng)過(guò)程動(dòng)態(tài)模擬[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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è)工程學(xué)報(bào)，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濾波的作物生長(zhǎng)模型與遙感數(shù)據(jù)的同化研究[D]. 北京：中國(guó)農(nóng)業(yè)科學(xué)院，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] 黃健熙，馬鴻元，田麗燕，等. 基于時(shí)間序列LAI和ET同化的冬小麥遙感估產(chǎn)方法比較[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 黃健熙，武思杰，劉興權(quán)，等. 基于遙感信息與作物模型集合卡爾曼濾波同化的區(qū)域冬小麥產(chǎn)量預(yù)測(cè)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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)量動(dòng)態(tài)集成預(yù)報(bào)方法[J]. 應(yīng)用氣象學(xué)報(bào)，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

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

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

秦其明，江蘇徐州人，教授，博士，博士生導(dǎo)師，主要研究方向?yàn)槎窟b感與地理信息系統(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標(biāo)定結(jié)合EnKF同化遙感和WOFOST模型模擬冬小麥時(shí)序LAI[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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