黃健熙 高欣然 黃 海 馬鴻元 蘇 偉 朱德海

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

0 引言

作物收獲期對作物產(chǎn)量與品質(zhì)影響顯著，收獲過早或過晚都不利于作物的豐產(chǎn)增收[1]。對于冬小麥而言，其主產(chǎn)區(qū)多采用冬小麥-夏玉米輪作的種植制度，收獲過晚將直接影響夏玉米的及時耕種。此外，隨著農(nóng)村土地制度改革的不斷推進，傳統(tǒng)的以家庭為單位的小農(nóng)經(jīng)濟正向大規(guī)模的機械化農(nóng)業(yè)轉(zhuǎn)型[2]，優(yōu)化收割次序、提前調(diào)度農(nóng)機，以實現(xiàn)作物的準時收獲是機械化農(nóng)業(yè)發(fā)展帶來的新問題，因此，在實際生產(chǎn)中迫切需要對作物成熟期進行大區(qū)域精準預報。

在農(nóng)學意義上來說，成熟期是指作物的營養(yǎng)存貯器官成長到可收獲的階段。影響成熟期變化的因素較為多樣，形成了不同的研究方法。文獻[3-4]著重分析氣象因素(如溫度、光周期、降水等)與作物成熟期的相關(guān)性，通過在單點尺度上建立氣象因子與成熟期的統(tǒng)計模型估測實驗年份的成熟期。段金省[5]利用多年的氣象數(shù)據(jù)和某年份作物的成熟期數(shù)據(jù)建立查找表，認為氣候條件越相近的年份生長的作物成熟期也越接近，并以此為理論基礎(chǔ)估測成熟期。基于氣象因子的估測方法易于操作，在使用氣象預報數(shù)據(jù)或區(qū)域氣候模型[6]的情況下能夠較快進行成熟期的預報，但是該類統(tǒng)計方法難以在實驗區(qū)域外推廣，且無法應對年際間種植品種和田間管理措施等因素的變化。

另一類方法是通過分析生育期內(nèi)作物生化參量的變化特征提取成熟期，這類方法多以遙感數(shù)據(jù)作為數(shù)據(jù)源，主要分為兩種技術(shù)路線：①分析生育期內(nèi)植被指數(shù)的時間序列曲線，計算指數(shù)時序曲線的曲率、波峰等特征，提取成熟期[7-8]，主要方法有：閾值法[9]、Logistic曲線擬合法[10]和滑動平均法[11]等。②監(jiān)測生育期內(nèi)作物水分[12]、葉綠素[13]、氮素[14]、干物質(zhì)[15]等指標的變化，使用回歸分析等統(tǒng)計學方法建立成熟期估測模型，進而實時監(jiān)測作物是否達到成熟所需指標。基于遙感數(shù)據(jù)監(jiān)測作生化參量提取成熟期的方法能夠在大區(qū)域上提取作物的成熟期，但是無法進行預報。此外，由于云、雨天氣的影響，作物生育期內(nèi)不可避免地會遇到遙感數(shù)據(jù)缺失或數(shù)據(jù)質(zhì)量較差的問題，將給時間序列曲線和生化指標監(jiān)測帶來誤差，影響成熟期提取精度，而且這類方法基于回歸分析構(gòu)建模型，對于品種調(diào)整和氣候變化缺乏適應性，難以在實驗區(qū)域外推廣應用。

作物成熟期受到氣候條件、作物品種特性和耕種時間的共同影響[16]，現(xiàn)有研究方法大多基于歷史數(shù)據(jù)對成熟期建模，不能很好地反映當年作物品種及耕種時間的變化。如果根據(jù)當年作物實際生長發(fā)育情況動態(tài)調(diào)整模型中與作物品種相關(guān)的參數(shù)，將能極大提高模型的適應性。因此，引入作物生長模型，并將遙感數(shù)據(jù)作為觀測數(shù)據(jù)，構(gòu)建代價函數(shù)優(yōu)化作物模型參數(shù)，不僅可以利用遙感數(shù)據(jù)大范圍、實時性強的優(yōu)點，而且又能以氣象預報數(shù)據(jù)驅(qū)動機理過程模型，實現(xiàn)區(qū)域尺度的成熟期預測。

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

1.1 研究區(qū)概況

選擇河南省冬小麥種植區(qū)作為研究區(qū)，其覆蓋范圍為31.38°～36.37°N, 110.35°～116.64°E，總面積約5.425×106hm2(中國國家統(tǒng)計局，2015年)，地處黃淮海平原，屬暖溫帶季風氣候區(qū)，年平均降水量約650 mm，年平均日照時數(shù)約2 300 h，四季分明，雨熱同期，光照充足。冬小麥是該地區(qū)種植的主要糧食作物之一，通常在10—11月播種，2—3月返青，5—6月成熟。由于氣象條件、播種時間和種植品種的差異，研究區(qū)內(nèi)冬小麥的成熟期在空間上和年際間也存在一定差異，適用于進行冬小麥成熟期的預測與驗證。本文所用冬小麥種植區(qū)通過隨機森林方法提取得到[17]，空間分辨率為500 m，在研究區(qū)內(nèi)包含35個國家級農(nóng)業(yè)氣象觀測站(圖1)。

圖1 研究區(qū)與農(nóng)業(yè)氣象站點分布Fig.1 Map of study area and agrometeorological sites

1.2 研究數(shù)據(jù)

1.2.1遙感數(shù)據(jù)

遙感數(shù)據(jù)來自2017—2018年MODIS LAI產(chǎn)品MCD15A3H (https:∥ladsweb.modaps.eosdis.nasa.gov/)，其時間分辨率為4 d，空間分辨率為500 m，通過使用ArcGIS工具對影像進行投影轉(zhuǎn)換、鑲嵌和裁剪等處理，輸出覆蓋整個研究區(qū)的MODIS LAI數(shù)據(jù)。

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

氣象數(shù)據(jù)來自2017—2018年歐洲中期天氣預報中心提供的氣象要素數(shù)據(jù)集以及集合預報格點數(shù)據(jù)集TIGGE(https:∥apps.ecmwf.int/datasets/data/tigge/levtype=sfc/type=cf/)，包括近地面氣溫、近地面氣壓、露點溫度、近地面全風速、太陽輻射、降水量等要素，其時間分辨率為6 h，空間分辨率為0.25°，其中預報數(shù)據(jù)的時限為16 d。對原始氣象數(shù)據(jù)進行預處理得到時間分辨率為1 d，空間分辨率為500 m的數(shù)據(jù)，包括：日最高溫、日最低溫、日照輻射量、水汽壓、平均風速與降水量6個要素，目前在多個領(lǐng)域得到廣泛應用[18-20]。超出預報數(shù)據(jù)時限的氣象數(shù)據(jù)參照本課題組之前的研究方法[21]，用歷史氣象數(shù)據(jù)代替。

1.2.3作物與土壤數(shù)據(jù)

作物和土壤數(shù)據(jù)分別來自農(nóng)業(yè)氣象站記錄的《作物生育狀況觀測記錄年報表》和《土壤水分觀測記錄年報表》，作物數(shù)據(jù)具體包括站點的區(qū)站號、臺站名稱、作物名稱、品種名稱、品種類型、作物關(guān)鍵生育期、臺站經(jīng)緯度、臺站海拔、株高、密度、田間管理措施等；土壤數(shù)據(jù)具體包括土壤質(zhì)量含水率、土壤相對濕度、土壤容重、田間持水量、凋萎濕度等。部分站點還包括關(guān)鍵生育期葉面積、植株干質(zhì)量等數(shù)據(jù)。

我國工程建設(shè)監(jiān)理制度的推行雖然走過了20多年，但有序競爭的市場環(huán)境還沒有形成。2010年，中國水利工程協(xié)會組織開展了首批水利建設(shè)市場主體信用評價工作，但采取的是工程施工和監(jiān)理市場主體自愿申請參加的方式，僅有137家水利建設(shè)市場主體自愿申報進行信用評級。其中監(jiān)理市場主體只有67家自愿申報進行信用評級，信用評級是在市場主體申報材料的基礎(chǔ)上，由全國性行業(yè)自律組織——中國水利工程協(xié)會進行初審、復審和組織評審委員會終審，信用評價缺乏強制性和約束性，信用評價結(jié)果缺乏公信力，實際的信用監(jiān)管效果還不夠理想，優(yōu)勝劣汰的市場準入和清出機制沒有形成，制約監(jiān)理活動良性發(fā)展的因素依然存在。

2 研究方法

2.1 WOFOST模型的本地化

WOFOST模型是由瓦赫寧根大學和瓦赫寧根研究中心共同開發(fā)的一種機理性作物生長模型，能夠模擬生育期內(nèi)作物對氣象和其他環(huán)境因子的反應，實現(xiàn)不同類型作物發(fā)育階段的量化、生物量和作物產(chǎn)量的模擬[22]，目前已經(jīng)廣泛應用于指導農(nóng)業(yè)生產(chǎn)和模擬作物長勢等多個方面[23-26]。為了使模型能夠準確地模擬研究區(qū)內(nèi)冬小麥的生長發(fā)育過程，首先需要對模型參數(shù)進行標定，實現(xiàn)WOFOST模型的本地化，使用鄭州農(nóng)業(yè)氣象試驗站的歷史觀測數(shù)據(jù)或查閱文獻[27-28]來確定模型的參數(shù)取值，部分參數(shù)校準值見表1。

2.2 MODIS LAI數(shù)據(jù)預處理

本文將研究區(qū)的MODIS LAI數(shù)據(jù)按照時間順序疊加，在各像元上生成對應的LAI時間序列曲線，使用Savitzky-Golay(S-G)濾波方法剔除MODIS LAI數(shù)據(jù)中由云、水汽、氣溶膠等帶來的噪聲。由于500 m分辨率影像受到混合像元因素的影響，且MODIS LAI數(shù)據(jù)的尺度低于WOFOST模型模擬LAI的尺度，為了減小系統(tǒng)誤差給同化造成的影響，采用歸一化的方法處理LAI數(shù)據(jù)，即把MODIS LAI和模型模擬LAI均映射到0～1內(nèi)，在保留MODIS LAI趨勢信息的情況下使兩種LAI數(shù)據(jù)源尺度相同，公式為

表1 WOFOST模型中主要參數(shù)校準值Tab.1 Calibrated values of main parameters in WOFOST

注：表中DVS表示作物的發(fā)育階段。

LAIn=(LAIo-LAImin)/(LAImax-LAImin)

(1)

式中LAIo、LAIn——歸一化前、后葉面積指數(shù)

LAImax、LAImin——該時間序列上葉面積指數(shù)的最大值和最小值

2.3 待優(yōu)化參數(shù)的選擇

由文獻[29-31]可得，對作物成熟期影響較大的作物參數(shù)有IDEM(出苗日期)、TBASEM(出苗下限溫度)、TSUMEM(從播種到出苗期有效積溫)、TSUM1(出苗期至開花期有效積溫)、TSUM2(開花期至成熟期有效積溫)等。由于TBASEM和TSUMEM均反映作物從播種期到出苗期的特性，而本文選擇LAI作為耦合變量，且遙感LAI數(shù)據(jù)只能反映作物出苗之后的生長狀況，因此無法有效地優(yōu)化播種期相關(guān)參數(shù)。分析歷史數(shù)據(jù)發(fā)現(xiàn)，在研究區(qū)TSUM1、TSUM2、IDEM 3個參數(shù)沒有表現(xiàn)出明顯的空間規(guī)律，且年際間變化較大，因此將三者設(shè)為待優(yōu)化參數(shù)，根據(jù)文獻[28-29]及歷史實測數(shù)據(jù)規(guī)定其初始值與上下限如表2所示。

2.4 SCE-UA算法與代價函數(shù)構(gòu)建

優(yōu)化算法采用復合形混合演化算法(Shuffled complex evolution-University of Arizona, SCE-UA)。該算法是DUAN等[32-33]在1993年提出的引入了復雜形分割、混合思想的全局優(yōu)化算法，在樣本空間的搜索效率、計算速率和全局搜索最優(yōu)的能力上表現(xiàn)突出且對待優(yōu)化參數(shù)初始值不敏感[34]，為遙感與作物生長模型同化的實際應用提供了可行的方法。

表2 待優(yōu)化參數(shù)初始值及取值范圍Tab.2 Initial value and range of parameters for optimization

由于之前的研究發(fā)現(xiàn)，與順序同化算法相比，SCE-UA算法更易受到觀測數(shù)據(jù)誤差的影響，當觀測LAI偏低時同化后的LAI偏低較嚴重[35]，因此在構(gòu)建代價函數(shù)時，先將遙感觀測LAI和WOFOST模擬LAI歸一化，得到代價函數(shù)

(2)

代價函數(shù)J為時間序列遙感觀測LAI與模型模擬LAI的誤差平方和。SCE-UA算法通過全局搜索，當滿足鄰近5個代價函數(shù)值均達到設(shè)定閾值或者迭代次數(shù)達到設(shè)定次數(shù)時停止運算并得到最優(yōu)參數(shù)集，將所得最優(yōu)參數(shù)集連同其他參數(shù)輸入WOFOST模型，實現(xiàn)模型模擬冬小麥生長發(fā)育過程和成熟期的預測。

3 結(jié)果與分析

3.1 WOFOST模型的標定

在本研究中，本地化的WOFOST模型是準確預測冬小麥成熟期的必要條件。耦合變量LAI是表征植被冠層結(jié)構(gòu)的基本參數(shù)，在WOFOST模型中是直接影響作物的光合、呼吸等重要過程的狀態(tài)變量[36]。因此，采用2016年鄭州和泛區(qū)農(nóng)業(yè)氣象站全生育期實測LAI數(shù)據(jù)對標定后的WOFOST模型模擬結(jié)果進行驗證。圖2是WOFOST模型模擬冬小麥LAI時間序列與實測LAI對比。

圖2 WOFOST模型模擬冬小麥LAI時間序列與實測LAI的對比Fig.2 Comparisons of WOFOST simulated time series LAI and field-measured LAI

冬小麥LAI在出苗期至越冬期前增長較緩慢，越冬期停止增長，返青期開始增長迅速，在開花期前達到最大值，開花期至成熟期逐漸降低，成熟期之后迅速降低。從圖2可以看出，標定后的WOFOST模型模擬的鄭州站和泛區(qū)站的LAI時間序列曲線的趨勢特征與實測值基本符合，經(jīng)計算，模擬LAI和實測LAI間的RMSE為0.92 m2/m2，說明模擬值與實測值具有良好的一致性，WOFOST模型模擬的冬小麥LAI精度較高。此外，部分點的模擬產(chǎn)生了一些偏差，可能的原因有：①年際間冬小麥生育期的變化導致WOFOST模型用2018年的作物參數(shù)無法模擬出較為準確的開花期和成熟期，影響了LAI的模擬精度。②模型標定是以鄭州站的觀測數(shù)據(jù)為準，因此不一定適用于其他地區(qū)冬小麥的模擬。

3.2 單點同化結(jié)果分析

開花期和成熟期分別是冬小麥營養(yǎng)生長和生殖生長階段的結(jié)束日期，是評價模型模擬發(fā)育期合理的關(guān)鍵指標。因此本文使用研究區(qū)內(nèi)農(nóng)業(yè)氣象站的開花期和成熟期觀測日期共同檢驗實驗結(jié)果的精度，圖3為2018年滑縣站同化前后模型模擬LAI和MODIS LAI的時間序列曲線。

圖3 滑縣站同化前后模型輸出LAI和MODIS LAI的時間序列曲線Fig.3 Changing curves of WOFOST simulated LAI, assimilated LAI and MODIS LAI in Huaxian station

從圖3可以看出，滑縣站同化后WOFOST模型模擬的出苗期、開花期與成熟期略晚于同化前。從LAI時序曲線的趨勢分析，同化后LAI曲線的變化趨勢更接近于MODIS LAI，表明WOFOST模型優(yōu)化TSUM1等參數(shù)能夠使其輸出的LAI曲線趨勢更接近遙感數(shù)據(jù)，而且并沒有受到遙感數(shù)據(jù)尺度較低的影響。相對應地，同化后LAI明顯高于同化前LAI和MODIS LAI。通過分析滑縣氣象數(shù)據(jù)和模型的逐日模擬結(jié)果，發(fā)現(xiàn)原因是同化前模型模擬的開花期在DOY110附近，由前文可知，冬小麥的快速生長時期發(fā)生在返青期至開花期之間，同化前返青期至開花期氣溫較低，光照不足，不利于冬小麥的快速生長；而同化后模擬的開花期較晚，為DOY122，該階段氣溫回升，使同化后的冬小麥在返青期至開花期期間長勢優(yōu)于同化前。由此看出，準確地模擬開花期不僅是成熟期模擬的基礎(chǔ)，更對WOFOST模型生物量、產(chǎn)量等多方面的模擬精度有很大的影響。

圖4 開花期、成熟期的預測值與實測值的回歸分析Fig.4 Regression analysis diagrams of predicted and measured anthesis and maturity date

在表3中，同化前開花期、成熟期取值范圍分別為DOY97～108、DOY135～152，同化后兩者取值范圍變?yōu)镈OY108～120、DOY141～156?？梢钥闯觯鄶?shù)點同化前的開花期和成熟期均早于同化后，且同化后開花期、成熟期的時間變異性大于同化前。對比各站點同化前后參數(shù)值，發(fā)現(xiàn)部分開花、成熟較晚的站點同化前開花期、成熟期均提前，說明標定的作物參數(shù)并不符合這些站點的實際情況，導致樣本內(nèi)部差異較小。同化前開花期、成熟期誤差較大也是因為在標定WOFOST模型時研究區(qū)內(nèi)的TSUM1、TSUM2等作物參數(shù)是基于歷史數(shù)據(jù)計算取值，沒有考慮到年際間品種與播種時間等因素的變化，同化后誤差明顯減小，說明同化MODIS LAI和WOFOST模型的方法可以提高研究區(qū)冬小麥開花期、成熟期的預測精度。

表3 部分站點同化前后開花期、成熟期日期對比Tab.3 Comparison of simulated and assimilated anthesis and maturity date DOY

開花期、成熟期的預測值與實測值的回歸分析見圖4。從圖4可以看出，同化后的開花期、成熟期與對應觀測日期基本符合，所成散點的趨勢線斜率均在1左右，說明同化后的開花期、成熟期與觀測日期具有良好的一致性，且RMSE分別為2.10 d和2.48 d，預測精度較高。誤差的主要來源有：①遙感數(shù)據(jù)空間分辨率較低，同化混合像元可能導致開花期模擬出現(xiàn)誤差，隨著模型的運行誤差逐漸累積，進而影響成熟期的預測精度。②由于研究區(qū)內(nèi)冬小麥品種和田間管理等方面存在差異，而大量參數(shù)在空間上的分布無法獲得，只能采用標定值，因此模型運行時其他參數(shù)取值相同也會帶來一定的誤差；使用SCE-UA算法最小化代價函數(shù)的過程中，可能很難準確地找到最接近實際情況的全局最優(yōu)參數(shù)集，因此優(yōu)化后的參數(shù)集也可能存在一定誤差。

3.3 開花期與成熟期空間分布

圖5 研究區(qū)同化前后開花期空間分布對比Fig.5 Comparison of anthesis date distribution before and after assimilation in study area

圖5、6分別為同化前后開花期、成熟期的空間分布圖。從圖5可見，由于只受氣象輸入數(shù)據(jù)的影響，同化前開花期模擬值呈較為明顯的區(qū)塊狀分布，相鄰的像元開花期時間十分接近。西北部和東部地區(qū)由于冬小麥播種時間較早且氣候適宜，因此開花期最早，南部、中部次之，而河南省西部海拔較高，北部緯度較高，因此氣溫相對較低，開花期最晚。同化后的開花期總體趨勢為從西南部向東北部逐漸推后，西部高海拔地區(qū)仍然最晚。同化前后開花期的空間分布變化主要表現(xiàn)在：①同化后研究區(qū)西南部開花期整體提前、研究區(qū)西北部和東部略有延后，原因是實際的IDEM、TSUM1數(shù)值與標定值相比發(fā)生了較大變化，導致同化后相應區(qū)域開花期普遍延后或提前。②同化后的開花期表現(xiàn)出了更詳細的空間變異性，在表現(xiàn)出開花期整體趨勢的前提下，可以看出相鄰冬小麥像元的開花期存在一定差異，說明用MODIS LAI優(yōu)化WOFOST模型的參數(shù)后可以在一定程度上表現(xiàn)出較細空間尺度上種植品種等方面的差異。

對比圖5a和圖6a可以看出，研究區(qū)中部、東部和西南部地區(qū)成熟期相對提前(位于研究區(qū)內(nèi)的先后順序)，說明這些地區(qū)當年度夏季氣溫較高，光照充足，因此生育期提早結(jié)束。從圖6可以看出，同化后成熟期的空間分布除了增加空間細節(jié)之外，還可以看出，中北部地區(qū)成熟期大幅延后，研究區(qū)內(nèi)成熟期整體趨勢為從西南部向東北部逐漸延遲，主要原因是優(yōu)化后的TSUM2(開花期至成熟期有效積溫)大于標定值，更符合該地區(qū)的實際情況，提高了模型的預測精度。

圖6 研究區(qū)同化前后成熟期空間分布對比Fig.6 Comparison of maturity date distribution before and after assimilation in study area

4 討論

通過數(shù)據(jù)同化框架耦合遙感數(shù)據(jù)和作物生長模型，并結(jié)合氣象預報數(shù)據(jù)實現(xiàn)區(qū)域尺度作物成熟期的預報，雖然總體精度較高，但仍存在一定的誤差，誤差產(chǎn)生的原因有：①預報延伸期氣象驅(qū)動數(shù)據(jù)存在誤差，由于5月16日(DOY136)后的氣象數(shù)據(jù)是使用歷史相似年份的同期氣象數(shù)據(jù)替代的，因此與實際氣象數(shù)據(jù)有一定偏差，以后可以嘗試天氣發(fā)生器[37-38]、區(qū)域氣候模型[6]等多種中長期氣象預報方法。②本方法僅優(yōu)化了WOFOST模型的3個輸入?yún)?shù)，其他參數(shù)取值在研究區(qū)內(nèi)相同，在一定程度上也限制了模型的預測精度，可以考慮增加對LAI敏感性較高、在生育期不同階段波動較大的參數(shù)作為待優(yōu)化參數(shù)，如AMAXTB(葉片最大CO2同化速率)、SLATB(比葉面積)等。

冬小麥開花期的模擬精度會直接影響其成熟期、LAI模擬的精度，誤差過大時甚至會影響LAI模擬的尺度，因此通過數(shù)據(jù)同化方法優(yōu)化WOFOST模型參數(shù)進而提高開花期模擬精度對于成熟期預測具有重要意義。而本文使用的MODIS LAI數(shù)據(jù)存在混合像元效應且容易受到云的影響，使用濾波的方法平滑MODIS LAI也會導致一定的誤差，因此為進一步提高成熟期預測精度，可以考慮增加高分辨率遙感數(shù)據(jù)或非光學遙感數(shù)據(jù)(如SAR數(shù)據(jù))進行多遙感數(shù)據(jù)同化。此外，可以考慮將量化開花期模擬的不確定性量化，在開花期模擬集合的基礎(chǔ)上進行同化，也可以加入作物參數(shù)擾動、多來源氣象預報數(shù)據(jù)分別構(gòu)建預報集合[39-40]，這些集合更有助于表征整個系統(tǒng)的不確定性，進而為決策者提供多維度的信息。

5 結(jié)論

(1) 與觀測數(shù)據(jù)相比，同化后WOFOST模型模擬的開花期RMSE為2.10 d，預測的成熟期RMSE為2.48 d，預測精度較高，說明基于MODIS LAI和WOFOST模型同化的方法在大區(qū)域作物成熟期預測方面可行。

(2) 基于歸一化思想構(gòu)建LAI代價函數(shù)，能夠優(yōu)化WOFOST模型的IDEM、TSUM1、TSUM2參數(shù)，優(yōu)化后的LAI時間序列曲線尺度不受MODIS LAI數(shù)據(jù)尺度較低的影響。

(3) 冬小麥開花期的模擬誤差會傳遞到成熟期的模擬中，因此開花期的準確模擬對于成熟期的精準預測具有重要意義，為了提高預測精度，使用高分辨率遙感數(shù)據(jù)進行同化或使用SAR數(shù)據(jù)等進行多目標數(shù)據(jù)同化是可能的解決途徑。