基于歷史氣象資料和WOFOST模型的區域產量集合預報

馬鴻元 黃健熙,2 黃 海 張曉東,2 朱德海,2

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

0 引言

作物生長模型是采用數學模型方法描述作物光合、呼吸、蒸騰、營養等機理過程，以特定時間步長來動態模擬作物生長發育期間的生理生化參數的數學模型，能較為準確地描述光、溫、水、肥等因子以及田間栽培和管理措施對作物生長和發育的影響[1-2]。自20世紀60年代作物模型研究領域創始以來，作物模型逐漸成為農學、土壤學、植物生理學、氣象學等學科的重要研究工具，尤其在作物產量模擬方面，作物模型能夠動態、定量地模擬作物生長發育和產量形成，因此在產量預報[3]、農業氣象服務[4]、產量的氣候變化響應模式[5]等領域發揮了重要作用。如歐盟委員會聯合研究中心(JRC)很早就將作物模型應用于作物監測和產量預報，在農業遙感監測項目(MARS)中建立了專門的作物生長監測系統(CGMS)，定時發布基于作物模型的產量月報[6]；美國[7]、印度[8]、巴西[9]等農業大國也都很早開展了基于作物模型的產量預報，我國李振海[10]、帥細強等[11]、陳上等[12]以預報或歷史氣象資料為基礎初步開展了DSSAT、Oyrza2000及CERES-Maize模型在作物產量預報中的應用研究。

隨著應用領域的不斷深入，目前作物模型已從最早的田間尺度模擬擴展到了大區域尺度的模擬[13-14]，在氣候變化情景模擬和遙感數據同化方面應用尤為廣泛，但由于作物模型大多按天積分運算，需要輸入全生育期的完整氣象要素才能模擬得到最終的產量結果。基于作物模型的產量預報方法主要還是受制于預報期內的氣象資料的獲取，因此根據不同的研究需要，形成了各種氣象資料的獲取方法。

在月度、季度預報中，一些研究直接利用區域氣候模型生成的數據作為作物模型的輸入，進而預測農作物產量[15]，由于生成的氣象數據存在一定的誤差，還需要對氣候模型的輸出進行糾偏和訂正[16-17]，同時低分辨率的氣候模型模擬的逐日天氣數據與氣象站點實測數據有較大的尺度差異，直接使用氣候模型產生的逐日氣象數據時不可避免地存在時間和空間尺度不匹配的問題[18]。

為了實現氣候模型與作物模型的高效耦合，必須將氣候/天氣模型的數據進行降尺度，常用的方法一是統計降尺度[19]和隨機天氣發生器法等，其中天氣發生器在作物模型中應用較廣，多項研究通過天氣生成器對各氣候模型的月平均數據降尺度生成逐日氣象要素，作為作物模型的驅動數據[20-21]；二是通過動力學模型降尺度[22]，該方法利用復雜的動力學方程把氣候模型的模擬結果降尺度到更小的網格和時間尺度上。動力降尺度模型考慮了局地的地形因子，從而能更好地反映出局地氣候特征，對極端天氣狀況有較好的模擬效果，但此方法計算復雜費時，應用在區域化的作物模型中成本較高。

而在中短期預報中，數值天氣模式則可以和作物模型相結合，BAIGORRIA等[16]使用COAPS的區域模型生成了20個數值天氣預報集合成員，耦合陸面模型后輸入到CERES-Maize模型，其中缺失的輻射量要素使用了天氣發生器生成，最終使用主成分回歸分析進行了產量預報。LEE等[23]則開發了一種氣象數據服務端，利用站點觀測數據和數值預報數據直接生成ORYZA 2000作物模型的氣象驅動數據，大幅提高了作物模型的部署效率。DUMONT等[24-25]使用STICS作物模型進行了大量的預報氣象驅動數據比較工作，驗證了使用歷史平均資料和天氣發生器等不同氣象驅動策略對當年小麥產量預報的影響，結果顯示使用兩種策略沒有明顯差異，使用天氣發生器精度稍好，但使用歷史平均資料的計算效率則更高。

我國大部地區為典型的季風氣候，年際變率大，且季節分布不均[26]，因此，基于氣象資料多年平均值作為預報期的驅動數據并不能很好地體現氣象要素的波動性，而且災害性的異常天氣會被平滑，這對模擬潛在的產量損失十分不利；而如果將歷史同期所有年份的氣象要素都作為集合成員進行驅動模擬，則可以重現本地歷史天氣情景下的產量[27]，通過集合預報的技術手段得出產量預報的不確定性，將以往單一數值的預報提升為概率預報[28-29]。

本文以多年氣象資料為依托，以Python平臺的WOFOST模型(Python crop simulation environment，PCSE)為基礎，構建使用歷史資料作為預報期氣象驅動數據的系統框架，區域化分布式模擬研究區內多年冬小麥產量，并系統分析產量集合的統計特征，結合農業氣象觀測資料與產量統計數據分析歷史集合預報方法的預報能力和精度，以期為區域化產量概率預報提供科學依據。

1 數據與方法

1.1 研究區概況

選取河北省保定市、衡水市冬小麥主產區為研究區(圖1)，該區地處黃淮海平原，屬暖溫帶半濕潤季風氣候，多年平均氣溫3～15℃，年降水量500～800 mm，四季分明，日照充足。冬小麥在本區糧食組成中占有較大比重。研究獲取了衡水市2009年29個觀測點的冬小麥田間觀測數據[30]，主要包括播種與收獲日期、物候期、葉面積指數(Leaf area index，LAI)和產量等，試驗點的分布充分考慮了該地區冬小麥的生長分布狀況，具有很好的代表性。

圖1 研究區示意圖Fig.1 Location of study region

1.2 氣象資料

所用的氣象驅動數據是中國科學院青藏高原研究所發布的中國區域高時空分辨率地面氣象要素驅動數據集(China meteorological forcing dataset)01.06.0014版本[31-32]，該驅動數據集包括近地面氣溫、地表氣壓、近地面空氣比濕、近地面全風速、降水率和向下長、短波輻射7個要素，空間分辨率為0.1°，時間分辨率為3 h。該數據集是以國際上現有的Princeton再分析資料、GLDAS資料、GEWEX-SRB輻射資料以及TRMM降水資料為背景場，融合了中國氣象局常規氣象觀測數據制作而成的。經過前期處理，將這些要素轉換為與作物模型量綱、格式一致的氣象輸入文件，包括日照輻射量、日最高溫、日最低溫、清晨水汽壓、風速和降水6個要素。

1.3 WOFOST模型

WOFOST(World food studies)模型是由荷蘭瓦赫寧大學和世界糧食研究中心共同開發的一種機理性作物生長模型[33]，該模型通過利用光截獲量和CO2同化量作為生長驅動過程，利用作物物候發育過程來描述作物生長，模擬的生長過程主要包括作物的碳同化、呼吸、蒸騰、干物質的分配等。模型主要有兩種應用模式，潛在模式下作物生長僅僅由溫度和太陽輻射決定，沒有考慮其他生長限制因素；水分限制模式下作物生長受水分限制。WOFOST模型是一種通用模型，通過定制不同的參數可以模擬不同種類的作物，具有很強的適用性，是大范圍環境條件下的各種作物模擬的有利工具。模型最初由Fortran語言開發，而今WOFOST研究組已經開發了適合Python平臺的版本PCSE，通過配置豐富的第三方庫，大大加強了模型應用的靈活性。

1.4 產量預報框架

圖2 基于歷史氣象資料的作物模型產量集合預報技術流程Fig.2 Flow chart of research technology route

傳統的作物模型運行過程只是輸入數據到輸出數據的單一映射，無法提供模型模擬的不確定性信息，但作物模型本身，以及輸入參數、氣象驅動存在著不確定性。如果以集合的思想將作物模型放入一個外部框架內，用不同年份的氣象驅動模擬的結果作為集合成員，可以使模擬結果的集合代表其概率分布，最終將單一數值的模擬輸出轉換為概率分布。基于上述考慮，本文引入多年(1979—2015年)歷史氣象資料，開發出以從某一預報時刻模擬該年作物生長和產量狀況的概率預報的集合預報框架。該框架以Python平臺下WOFOST(PCSE)模型為基礎，作物品種參數和區域化由本課題組歷史標定成果[34-39]確定，通過外部程序框架的搭建，完成多年氣象輸入文件生成、預報時刻氣象集合成員文件生成、分布式PCSE模擬、集合預報成員后處理等主要功能，從而實現區域化集合預報作物長勢和產量的主要功能(圖2)。利用該預報框架開展產量集合預報時，首先進行前期數據準備，預先收集作物品種參數、田間管理信息及播種至預報時刻的氣象資料，然后根據預報時刻生成模型的氣象輸入文件集合，該集合的每個成員文件由播種日至預報時刻的實際氣象和預報時刻至成熟日的歷史資料拼接而成，利用該輸入集合進行區域化的分布式模擬，最終得到模擬產量的集合并進行相應的后處理過程。

1.5 驗證方法與誤差分析

以1979—2015年氣象資料為歷史集合進行集合預報，以2009年29個觀測站的田間實測產量為觀測值，以2009年每月1日為預報起始時刻，同一日期下歷史集合模擬結果的平均值作為預報值，觀測值與預報值的差值即為預報誤差。

集合預報給出的預報是一個概率分布或者是一個包含多個成員的集合，若需要對集合預報進行確定性驗證，一般是對集合成員的平均值或集合中位數進行驗證[40]。集合成員的平均值計算公式為

(1)

Mi——集合成員值

n——集合成員數量

當有了成員集合后，還可進行某一產量水平的概率預報，假設集合中所有成員的權值相等，計算某一事件的成員數與集合成員總數的比例，即為該事件發生的概率，計算公式為

(2)

式中pi——模擬結果中某一事件(如產量高于7 000 kg/hm2)出現與否的標識，當該事件發生時pi=1，否則為0

平均差會由于差值的正負抵消，因此本文采用平均絕對誤差(Mean absolute error，MAE)表示誤差，另外還將均方根誤差(Root mean square error，RMSE)與皮爾遜相關系數(Pearson correlation coefficient，PCC)作為精度評價指標，相關計算式為

(3)

(4)

(5)

式中m——驗證樣本數量

i——驗證樣本序號

Fi——第i實測點模型預報值

Oi——第i實測點觀測值

其中，MAE越小說明預報越精確，PCC越接近于1說明預報模型的預報能力越強。通過對每個觀測點進行1979—2008年的集合模擬，對預報的誤差進行分布擬合，最后統計3種指標作為預報能力。

2 結果分析

2.1 單點預報分析

圖3為深州市2009年3月1日(圖3a)、4月1日(圖3b)、5月1日(圖3c)、6月1日(圖3d)的葉面積與籽粒干物質量預報結果，其中不同顏色的線條表示不同集合成員的全生育期模擬值，左上方藍色直方圖分別表示6月1日的LAI與最終產量的集合成員分布頻率。紅色虛線標記的是該點的實測產量(6 786 kg/hm2)，藍色豎線指示預報日期，該時刻之前的氣象驅動資料為2009年實際氣象數據，該天及之后的氣象驅動為1979—2008年的歷史集合成員。結果顯示預報開始后集合成員的模擬結果逐漸出現分化，LAI的集合分布在5.0～6.5之間，產量的集合分布在6 000～8 300 kg/hm2之間，不同年份的模擬結果有較大的差異，但隨著預報時刻的不斷更新，LAI的集合成員模擬值變化幅度減小，而且基于歷史氣象的模擬集合其變化區間能夠較好地涵蓋實測值，說明歷史氣象數據的集合化表達可以用來部分表達模型模擬的不確定性，在小麥拔節后期即可得到對小麥長勢、LAI較為準確的集合預報。

模型模擬的產量即存貯器官干物質質量(Total dry weight of storage organs，TWSO)的集合平均值與實測值相差始終在500 kg/hm2左右，而直到5月初產量模擬結果仍有很大變異性(圖4)，這是因為冬小麥生長前期是以營養生長為主，同化的干物質主要分配到葉、根、莖中；生育中期是小麥營養生長和生殖生長并進時期；在小麥生長后期同化物已停止向營養器官中分配，主要分配給貯藏器官；同時根、莖、葉等營養器官中的干物質會有部分轉移到穗部形成產量，因此在抽穗后小麥籽粒開始形成、積累大量干物質并進行器官間的養分轉移直至成熟，所以在冬小麥的產量形成中抽穗—灌漿的生殖生長階段十分關鍵，這一時期作物遭受氣象脅迫作用會十分明顯地影響產量，這也說明在后續實際業務預報中，在冬小麥生長早期集合成員模擬的變異性較大，此時進行預報效果較差，而最終產量預報結果對抽穗—灌漿階段的氣象預報最為敏感，從生育后期即5月開始進行華北地區產量預報有一定的預報價值。這段時間氣象數據的不確定性將決定產量預報結果的成敗，因此在生長早期進行業務產量預報時還需要結合中長期的數值天氣預報產品。

圖3 深州市觀測點4個預報時期的LAI和TWSO集合模擬值與實測值對比Fig.3 Comparison of simulated and measured values of LAI and TWSO in four forecast periods of Shenzhou City

圖4 2009年深州市實測點不同起報時刻的冬小麥產量滾動集合預報及其箱形圖Fig.4 Rolling ensemble forecasts of winter wheat yields in 2009 and their box plots at different forecasting times in Shenzhou City

2.2 區域預報分析

通過分布式區域化模擬預報，得到研究區內區域化的產量集合預報結果。為了對產量預報結果進行驗證，本文以預報的集合中值產量和實測點產量進行對比驗證，可以看出隨著預報日期的推進，研究區南部和北部的產量差異加大，南部預報的產量有所提高，北部產量逐漸降低，說明2009年該區域南部整體天氣條件較好，實際產量優于歷史條件下的模擬(圖5)，而北部區域天氣條件差于歷史同期。從該研究區產量大于7 500 kg/hm2的區域化概率預報(圖6)，可以看出在2月初開始預報時，根據歷史同期氣象資料，大部分地區尤其是保定地區的產量高于7 500 kg/hm2的概率低于50%，然而隨著預報時刻的推進，絕大多數地區產量高于7 500 kg/hm2的概率大于70%，說明2009年氣象條件適宜冬小麥生長，天氣條件優于歷史同期。

2.3 預報性能定量分析

圖5 2009年8個預報時期下區域化集合產量預報(集合中位數)Fig.5 Regionalized yield ensemble forecast (median of ensemble) for eight forecast periods in 2009

圖6 2009年8個預報時期下產量高于7 500 kg/hm2的概率預報Fig.6 Probability forecast for event of production above 7 500 kg/hm2 on eight forecast dates in 2009

為了進一步評估集合預報的性能和誤差，本文對所有田間實測點的產量集合預報進行了統計，根據實測點的實測產量計算16個預報時期集合中位數的驗證指標MAE、RMSE、PCC。16個預報時期的MAE為458～578 kg/hm2, RMSE在598～758 kg/hm2之間，PCC在0.503～0.563之間(圖7)，集合平均數的情況與之類似，但由于篇幅所限未列出。總體而言基于歷史氣象的集合預報對于本年的產量有一定的預報能力，相對誤差在6.6%左右，集合中位數與實測值相關性也較好，說明歷史氣象集合的變異性很好地表達了氣象輸入參數的不確定性，有助于概率預報和數據等應用。但隨著預報期的推進集合中位數的誤差指標收斂較慢，直到5月中旬左右預報的性能達到最大值，說明氣象數據的不確定性仍然不足以完全表現實際產量的不確定性，尤其是潛在模式下模型對于水分脅迫和土壤特性的影響沒有得到表達，而且受制于氣象驅動數據空間分辨率的限制，對于精細尺度的產量無法僅憑氣象數據進行預報。

圖7 衡水市實測點不同起報時刻下產量滾動預報的集合中位數驗證結果Fig.7 Verification results of rolling ensemble forecast median of 29 observation points in Hengshui City at different forecasting times

3 討論

氣象數據是作物模型的主要驅動數據，利用作物模型模擬作物生長需要完整描述整個生育期的氣象數據，缺乏預報期氣象驅動數據是制約作物模型應用于產量預報的主要瓶頸，因而目前作物模型的應用主要局限于歷史產量估測、單點模擬驗證等，而作物模型機理明確、可擴展性強的優勢沒有得到發揮。氣象歷史資料對于本地的氣候模式和變異性有很好的描述能力，因此可以使用歷史資料作為預報期的輸入資料。本文采用潛在模式下的WOFOST模型作為預報模型，對2009年衡水保定地區的冬小麥產量進行了模擬預報，并使用集合預報的方法進行了后處理，結果發現基于歷史氣象的集合預報可以進行本年產量的預報，相對誤差在6.6%左右，預報精度較好，表明該方法可以作為缺乏預報期氣象資料時的補充手段，而且集合成員的分布可以表現產量的不確定性，可以根據集合成員提供產量水平的概率預報。

本文從3月開始進行滾動預報，3月初模型即有一定預報能力(PCC為0.505, MAE 為504 kg/hm2)，說明在本地的歷史多年氣象資料集合可以代表當年的氣象要素變率；而且預報誤差在后期逐漸減小，說明隨著氣象輸入數據不斷被更新為實測值，氣象數據的不確定在下降；但產量集合的收斂過程較慢，直到灌漿中后期達到最高預報精度，這可能與模型模擬的空間尺度較粗而且模型其他參數未能區域化有關，另外在潛在模式下作物生長的模擬未考慮水分等脅迫作用的影響，區域上的產量差異僅僅來自于氣象輸入的差異，因而使得在空間上產量的分布梯度較不明顯。

由于作物品種、土壤特性、管理措施對于產量的變異性有很強的影響，各種輸入參數的變異性共同構成了產量的變異性，為了進一步提高產量區域化預報的精度，需要對敏感參數進行區域化，然而目前作物模型區域化調參受制于數據缺乏因而難度較大，將來借助農業氣象觀測、土壤數據庫等數據進行參數區域化是進行模型產量預報的重要前期過程。另外衛星遙感數據是柵格類型的觀測數據源，具有覆蓋面積廣、空間連續、周期觀測、相對成本低的優勢，通過遙感觀測和作物模型的數據同化可以獲得大區域范圍的模型參數分布[41]，因此數據同化將會在模型產量預報方面有很大的潛力。

作物生長模型是一個非線性的復雜系統，通過模擬可以給出產量和其他生長參數的確定性模擬值，而集合預報的技術手段可以得到預報值不確定性信息，本文采用集合預報的思想進行了作物模型模擬，集合成員的生成方式是歷史多年的氣象輸入，而作物模型的不確定性來源除了氣象輸入外，還有模型本身、輸入參數的不確定性，因此多模式的作物模型、作物參數擾動都可以用來生成預報集合[42-43]，這樣的超級集合更有助于表達整個系統的不確定性，為提高模擬精度、概率預報和數據同化[44]提供了一條很好的途徑。和數值天氣預報的技術發展歷史相似，從單值的模擬拓展為概率預報是產量預報的發展方向，因此使用集合預報的思想進行作物模型模擬具有很好的應用前景。

本文初步探討集合預報手段在作物模型模擬中的技術流程和應用前景，因此僅通過較小規模研究區一年的模擬預報對技術方案的可行性和誤差進行了實驗，而本文的技術方案還需要更多驗證，多年大區域的驗證實驗還有待進行，而且解決預報期天氣輸入的根本方法在于結合數值天氣預報，對于集合數值預報和作物模型的結合需要更深一步的研究。此外，在氣候變化的大背景下，近年來的極端、異常天氣事件增多，需要采用其他方法得到預報期更高精度的氣象輸入集合，如歷史相似年份替代、天氣發生器等技術手段。

4 結論

(1)使用多年氣象資料作為驅動建立預報集合，集合的平均值與實測產量PCC最高為0.563，MAE最低為458 kg/hm2,表明歷史氣象集合對于本年度的作物產量具有一定的預報能力，能夠作為預報期的氣象數據補充手段，而集合的平均值或中位數能否代表當年實際氣候變異下的模擬產量，取決于在歷史氣象集合中當年氣象條件的一般性。

(2)抽穗至灌漿時期是作物產量形成的關鍵期，冬小麥產量對于這一階段的氣象條件十分敏感，在抽穗后開始產量預報有最高的預報精度。因此為進一步提高基于作物模型的產量預報精度，需要準確的抽穗至灌漿時期的氣象參數，作物模型與中長期數值預報結合是提高預報時效的必然要求。

(3)現階段研究在潛在模式下模擬預報，產量的變異性只與氣象數據中的溫度和輻射量有關，因此驅動數據的低空間分辨率和輸入參數無區域差異是制約預報精度的關鍵因素，使得模擬產量的空間梯度較小，預報精度難以進一步提高，為了解決輸入參數的估計問題，與遙感觀測進行數據同化是可能的解決途徑。