赫迪，胡鵬程，王靖，王良，房世波

(1. 中國氣象科學(xué)研究院災(zāi)害天氣國家重點(diǎn)實(shí)驗(yàn)室，北京 100082；2. 澳大利亞昆士蘭大學(xué)農(nóng)業(yè)與食品科學(xué)學(xué)院，布里斯班 QLD4067；3. 中國農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院，北京 100193；4. 山東省農(nóng)業(yè)科學(xué)院，山東 濟(jì)南 250100)

食用油供給安全是我國糧食安全戰(zhàn)略的重要組成部分，當(dāng)前國內(nèi)食用油市場(chǎng)需求不斷增長(zhǎng)，然而我國油料作物的單產(chǎn)和種植面積卻一直徘徊不前[1-4]，食用油自給率低于30%[5]，食用油供給形勢(shì)嚴(yán)峻。 油菜是我國第二大油料作物，除了用于高品質(zhì)的食用油生產(chǎn)外，也正被廣泛應(yīng)用于生物能源等工業(yè)用油領(lǐng)域[6-9]，因此發(fā)展油菜生產(chǎn)對(duì)保障我國食用油供給安全有重要意義。

長(zhǎng)江流域是我國油菜的優(yōu)勢(shì)主產(chǎn)區(qū)，其冬油菜種植主要以水稻-油菜輪作形式利用冬季閑田進(jìn)行，不影響其他主糧作物的生產(chǎn)。 長(zhǎng)江流域可擴(kuò)大利用的油菜種植面積近399 萬hm2[10-11]，還有很大的發(fā)展空間。 定量評(píng)估長(zhǎng)江流域冬油菜產(chǎn)量潛力分布特征，明確油菜優(yōu)勢(shì)產(chǎn)區(qū)，對(duì)如何擴(kuò)大推廣冬油菜在長(zhǎng)江流域的種植有重要指導(dǎo)作用。

相比于傳統(tǒng)的光溫生產(chǎn)潛力和氣候生產(chǎn)潛力逐級(jí)訂正法，作物生長(zhǎng)模型APSIM 可以在逐日尺度上考慮作物生長(zhǎng)與各環(huán)境因子的關(guān)系，從而更準(zhǔn)確地獲得作物的潛在產(chǎn)量，現(xiàn)已被廣泛應(yīng)用于小麥、玉米、油菜、馬鈴薯等的產(chǎn)量評(píng)估[12-17]，但由于模型運(yùn)行需要的參數(shù)和輸入量較多，大范圍應(yīng)用受到限制。 機(jī)器學(xué)習(xí)近年來被廣泛應(yīng)用于農(nóng)業(yè)領(lǐng)域[18]，利用作物模型的輸入和輸出變量，構(gòu)建機(jī)器學(xué)習(xí)模型反演作物模型的模擬結(jié)果，可簡(jiǎn)化作物模型的輸入需求并提高計(jì)算效率。

本研究以長(zhǎng)江流域冬油菜產(chǎn)區(qū)為研究區(qū)，將作物生長(zhǎng)模型APSIM 和機(jī)器學(xué)習(xí)技術(shù)結(jié)合，建立油菜產(chǎn)量潛力預(yù)測(cè)模型，利用該模型評(píng)估長(zhǎng)江流域縣級(jí)冬油菜產(chǎn)量潛力，并對(duì)產(chǎn)量潛力的空間分布特征進(jìn)行分析，同時(shí)結(jié)合縣級(jí)產(chǎn)量的歷史統(tǒng)計(jì)數(shù)據(jù)定量當(dāng)前不同產(chǎn)量潛力區(qū)域的產(chǎn)量差，以期為擴(kuò)大長(zhǎng)江流域冬油菜生產(chǎn)提供科學(xué)指導(dǎo)。

1 研究區(qū)域與研究方法

1.1 研究區(qū)域和數(shù)據(jù)來源

長(zhǎng)江流域指長(zhǎng)江干流和支流流經(jīng)的廣大區(qū)域，位于東經(jīng)90°33＇～122°25＇、北緯24°30＇～35°45＇之間，氣候主要受副熱帶季風(fēng)環(huán)流影響，以濕潤(rùn)的亞熱帶季風(fēng)氣候?yàn)橹鳌?長(zhǎng)江流域冬油菜的主要種植區(qū)分布于長(zhǎng)江中下游地區(qū)，以水稻-油菜輪作為主，冬油菜秋季播種，次年春季收獲[19]。本研究選取長(zhǎng)江中下游地區(qū)的玉溪、樂山、江口、長(zhǎng)沙、武漢和南京6 個(gè)代表站點(diǎn)(表1)，收集油菜生長(zhǎng)發(fā)育、田間管理、生物量和產(chǎn)量數(shù)據(jù)；其氣象數(shù)據(jù)來源于中國氣象局科學(xué)數(shù)據(jù)共享網(wǎng)，包括各個(gè)試驗(yàn)站點(diǎn)試驗(yàn)?zāi)攴菁?961—2010 年的逐日氣象資料。 氣象要素包括平均氣溫、最高氣溫、最低氣溫、降水量、日照時(shí)數(shù)、平均風(fēng)速，逐日太陽總輻射值由日照時(shí)數(shù)根據(jù)Angstrom 公式計(jì)算獲得。

表1 各站點(diǎn)的地理信息、氣候特征及油菜品種

縣域尺度的氣象資料(1985—2021 年)來自于美國航天航空局(NASA)的Prediction of Worldwide Energy Resource(POWER)數(shù)據(jù)庫(https:/ /power.larc.nasa.gov/)，氣象要素包括最高氣溫、最低氣溫、降水量、太陽總輻射。

1.2 APSIM-Canola 模型

作物生長(zhǎng)模型APSIM(Agricultural Production Systems sIMulator)可以逐日模擬溫度、光周期、輻射、土壤水分和氮素水平的變化對(duì)作物的發(fā)育過程、生長(zhǎng)過程以及籽粒產(chǎn)量的影響，目前已廣泛應(yīng)用于農(nóng)業(yè)生產(chǎn)系統(tǒng)管理、氣候變化影響評(píng)估和產(chǎn)量潛力評(píng)估等研究中。 APSIM 的油菜模擬模塊(APSIM-Canola)對(duì)油菜的發(fā)育期模擬主要由積溫驅(qū)動(dòng)，考慮了春化作用和光周期對(duì)開花的影響；生物量的積累主要根據(jù)冠層截獲的太陽輻射量和光能利用率計(jì)算，并考慮每天的溫度、水分和氮素對(duì)生物量積累的影響；籽粒產(chǎn)量是利用動(dòng)態(tài)收獲指數(shù)進(jìn)行模擬，在灌漿階段，收獲指數(shù)線性增加。APSIM-Canola 在中國的適用性和參數(shù)化過程已經(jīng)得到充分的驗(yàn)證和分析[14，20-21]，模型中影響油菜發(fā)育期、生物量和產(chǎn)量的重要參數(shù)包括春化敏感性、光周期敏感性、灌漿周期、光能利用效率和最大收獲指數(shù)。

模型參數(shù)優(yōu)化選擇了全局優(yōu)化的方法，即首先根據(jù)參數(shù)范圍，基于參數(shù)在范圍內(nèi)是均勻分布的假設(shè)生成全部的參數(shù)組合，即為目標(biāo)函數(shù)的可能解，然后從中選取使目標(biāo)函數(shù)最小的參數(shù)組合，即為目標(biāo)函數(shù)的最優(yōu)解。 本研究選取均方根誤差(RMSE)和標(biāo)準(zhǔn)均方根誤差(NRMSE)作為目標(biāo)函數(shù)對(duì)模型模擬結(jié)果進(jìn)行評(píng)價(jià)，計(jì)算方法如下:

其中，Oi和Si是相對(duì)應(yīng)的觀測(cè)值和模擬值是觀測(cè)值的平均值，n 是觀測(cè)樣本的數(shù)量。

1.3 基于機(jī)器學(xué)習(xí)技術(shù)的產(chǎn)量潛力預(yù)測(cè)模型構(gòu)建

產(chǎn)量潛力預(yù)測(cè)模型的構(gòu)建主要分為兩步，首先利用APSIM 模型在代表站點(diǎn)進(jìn)行多年(1961—2010 年)無肥力脅迫的產(chǎn)量潛力模擬，并統(tǒng)計(jì)每年的關(guān)鍵氣象信息，包括生長(zhǎng)季(10—12 月、次年1—5 月)及逐月的平均太陽輻射、最低氣溫、最高氣溫、總降雨量，獲得長(zhǎng)期多點(diǎn)位的氣象和產(chǎn)量數(shù)據(jù)集；然后將全部數(shù)據(jù)的80%作為訓(xùn)練集、20%作為驗(yàn)證集，利用訓(xùn)練集的資料進(jìn)行機(jī)器學(xué)習(xí)模型的訓(xùn)練。 將數(shù)據(jù)集中的氣象資料作為輸入值、產(chǎn)量數(shù)據(jù)作為輸出值，利用機(jī)器學(xué)習(xí)算法訓(xùn)練機(jī)器學(xué)習(xí)模型，訓(xùn)練過程中用10-fold 交差驗(yàn)證的方法對(duì)模型進(jìn)行調(diào)參，最后用驗(yàn)證集的數(shù)據(jù)對(duì)機(jī)器學(xué)習(xí)模型進(jìn)行驗(yàn)證。 本研究通過比較3 種常用的機(jī)器學(xué)習(xí)算法——神經(jīng)網(wǎng)絡(luò)(kNN)、支持向量機(jī)(SVM)和隨機(jī)森林(Random Forest)，最終選用建模效果最好的隨機(jī)森林算法進(jìn)行建模。 機(jī)器學(xué)習(xí)模型的訓(xùn)練全部在R 語言程序中調(diào)用caret 程序包實(shí)現(xiàn)。

1.4 長(zhǎng)江流域冬油菜產(chǎn)量潛力預(yù)測(cè)和產(chǎn)量差計(jì)算

利用訓(xùn)練好的機(jī)器學(xué)習(xí)模型和縣域尺度的氣象資料(1985—2021 年)對(duì)長(zhǎng)江流域近37 年的油菜潛在產(chǎn)量進(jìn)行預(yù)測(cè)，并對(duì)其多年平均產(chǎn)量的區(qū)域特征進(jìn)行分析，分析由R 語言程序中的sf 程序包實(shí)現(xiàn)。 對(duì)有觀測(cè)資料記錄的縣級(jí)油菜產(chǎn)量(2008—2012 年)進(jìn)行提取，并與預(yù)測(cè)的潛在產(chǎn)量對(duì)比分析，明確當(dāng)前長(zhǎng)江流域油菜產(chǎn)量差的分布特征。

2 結(jié)果與分析

2.1 APSIM-Canola 模型驗(yàn)證

通過全局優(yōu)化的方法進(jìn)行模型參數(shù)調(diào)試后，APSIM-Canola 模型可以準(zhǔn)確模擬6 個(gè)代表性站點(diǎn)7 個(gè)油菜品種在不同播期和密度種植條件下的發(fā)育期、生物量和產(chǎn)量(圖1)。 其中，對(duì)開花期和成熟期的模擬誤差分別為2.4 d 和5.1 d，相對(duì)誤差分別為1.7%和2.4%；對(duì)生物量和產(chǎn)量的模擬誤差分別為2 168.6 kg/hm2和332 kg/hm2，相對(duì)誤差分別為20.6%和14.3%。 可見，生物量的模擬誤差略高，這主要由生物量的測(cè)量誤差較高及部分試驗(yàn)的生物量是由產(chǎn)量和收獲指數(shù)轉(zhuǎn)換獲得而造成。

圖1 APSIM-Canola 模型對(duì)6 個(gè)站點(diǎn)油菜發(fā)育期、生物量和產(chǎn)量的實(shí)測(cè)值與模擬值比較

2.2 基于機(jī)器學(xué)習(xí)技術(shù)的產(chǎn)量潛力預(yù)測(cè)模型驗(yàn)證

機(jī)器學(xué)習(xí)模擬可以準(zhǔn)確復(fù)現(xiàn)APSIM-Canola模型在6 個(gè)代表站點(diǎn)的長(zhǎng)期模擬結(jié)果。 驗(yàn)證數(shù)據(jù)集的驗(yàn)證結(jié)果(圖2)顯示，隨機(jī)森林模型可以解釋81% 的產(chǎn)量潛力變異， 模擬誤差為475 kg/hm2，相對(duì)誤差為16.4%。 該結(jié)果與APSIMCanola 模型對(duì)觀測(cè)數(shù)據(jù)的模擬精度接近，進(jìn)一步說明機(jī)器學(xué)習(xí)模型可以用于大區(qū)域產(chǎn)量潛力的預(yù)測(cè)。

圖2 隨機(jī)森林模型對(duì)APSIM-Canola 模型模擬的油菜產(chǎn)量潛力的復(fù)現(xiàn)和驗(yàn)證

機(jī)器學(xué)習(xí)建模過程中，通過對(duì)輸入因子的重要性進(jìn)行分析，最終選取6 個(gè)預(yù)測(cè)變量進(jìn)行建模，包括生長(zhǎng)季(10 月—次年5 月)的平均最高溫度、平均最低溫度、降雨量及7 月平均最低溫度、1 月平均最高溫度和2 月平均太陽輻射。 在縣域尺度上，利用上述6 個(gè)變量和隨機(jī)森林模型即可對(duì)每年的油菜產(chǎn)量潛力進(jìn)行預(yù)測(cè)。

2.3 長(zhǎng)江流域冬油菜的產(chǎn)量潛力和產(chǎn)量差分布特征

1985—2021 年長(zhǎng)江流域冬油菜平均產(chǎn)量潛力分布特征如圖3 所示。 冬油菜潛在產(chǎn)量分布大致呈西南低、東北高的趨勢(shì)，平均產(chǎn)量潛力變化范圍在1 000～4 000 kg/hm2之間。 云南地區(qū)的平均潛在產(chǎn)量?jī)H有1 000 ～1 500 kg/hm2，這可能是因?yàn)楸狙芯磕M的潛在產(chǎn)量為雨養(yǎng)潛在產(chǎn)量，該地區(qū)潛在產(chǎn)量較低可能與西南地區(qū)季節(jié)性干旱有關(guān)。 江蘇、安徽、湖北、重慶大部分地區(qū)及江西、湖南、貴州的北部地區(qū)為冬油菜的主要優(yōu)勢(shì)產(chǎn)區(qū)，多年平均潛在產(chǎn)量均可達(dá)到3 000 kg/hm2以上。

圖3 1985—2021 年長(zhǎng)江流域縣級(jí)冬油菜平均潛在產(chǎn)量空間分布特征

2008—2012 年平均的縣級(jí)統(tǒng)計(jì)資料顯示，統(tǒng)計(jì)產(chǎn)量的空間分布特征與潛在產(chǎn)量的分布特征基本一致，也呈北高南低、西低東高的趨勢(shì)，進(jìn)一步說明了機(jī)器學(xué)習(xí)模型預(yù)測(cè)的潛在產(chǎn)量分布的合理性。 產(chǎn)量差的分布特征則與統(tǒng)計(jì)產(chǎn)量和潛在產(chǎn)量不同，在江蘇大部分地區(qū)、安徽和湖北北部、四川大部分地區(qū)、云貴高原地區(qū)，油菜產(chǎn)量差大多低于1 000 kg/hm2，說明這些地區(qū)的品種選擇和農(nóng)田管理均接近產(chǎn)量最優(yōu)水平，進(jìn)一步的產(chǎn)量提高需要優(yōu)化品種選育。 而在江西、湖南大部分地區(qū)及安徽和湖北北部，產(chǎn)量差可達(dá)2 000 kg/hm2，縣級(jí)統(tǒng)計(jì)產(chǎn)量?jī)H有潛在產(chǎn)量的50%左右，有很大的產(chǎn)量提升空間(圖4)。

圖4 2008—2012 年長(zhǎng)江流域縣級(jí)冬油菜平均統(tǒng)計(jì)產(chǎn)量、潛在產(chǎn)量和產(chǎn)量差空間分布特征

3 討論

目前我國長(zhǎng)江流域適于油菜種植的冬閑田約有778 萬hm2[11]，明確長(zhǎng)江流域的油菜優(yōu)勢(shì)產(chǎn)區(qū)及產(chǎn)量差的分布特征，可以科學(xué)指導(dǎo)和推廣該區(qū)域的油菜種植。 本研究結(jié)果表明長(zhǎng)江中下游的冬油菜潛在產(chǎn)量有很大的區(qū)域性差異，長(zhǎng)江中下游北部區(qū)域有更高的產(chǎn)量潛力，多年平均潛在產(chǎn)量可達(dá)3 000 ～4 000 kg/hm2。 張智[22]的研究表明長(zhǎng)江流域冬油菜的實(shí)際產(chǎn)量、試驗(yàn)產(chǎn)量、可獲得產(chǎn)量、潛在產(chǎn)量分別為1 810、2 437、3 702、3 996 kg/hm2，與本研究結(jié)果基本一致。 而Tian 等[11]利用改進(jìn)后的農(nóng)業(yè)生態(tài)區(qū)域模型(AEZ)分析的長(zhǎng)江流域油菜潛在產(chǎn)量分布則與本研究結(jié)果有所差異，其研究指出該流域油菜的潛在產(chǎn)量高值區(qū)主要分布在四川、云南及江西和湖南南部，與實(shí)際的縣級(jí)統(tǒng)計(jì)產(chǎn)量分布有所差異。 叢日環(huán)等[23]對(duì)長(zhǎng)江流域與油菜產(chǎn)量相關(guān)的氣候因子的分析研究發(fā)現(xiàn)，AEZ 模型獲得的潛在產(chǎn)量分布特征與長(zhǎng)江流域大于等于0 ℃積溫的分布模式十分接近；進(jìn)一步研究表明，長(zhǎng)江上游和下游地區(qū)大于等于0 ℃積溫與油菜產(chǎn)量并無顯著相關(guān)關(guān)系，但在中游地區(qū)呈負(fù)相關(guān)關(guān)系。 這可能是由于長(zhǎng)江流域日照時(shí)數(shù)和太陽輻射相對(duì)充足，太陽輻射不是油菜生長(zhǎng)的限制因子。 AEZ 模型是根據(jù)聯(lián)合國農(nóng)糧組織(FAO)和國際應(yīng)用系統(tǒng)分析研究所(IIASA)基于中國1961 年以來的統(tǒng)計(jì)資料共同開發(fā)的大尺度統(tǒng)計(jì)模型，其模擬結(jié)果可能更接近傳統(tǒng)的光照、溫度、水分逐級(jí)訂正的生產(chǎn)潛力計(jì)算結(jié)果。 而本研究預(yù)測(cè)的長(zhǎng)江流域油菜潛在產(chǎn)量分布可能與實(shí)際情況更為接近。

本研究根據(jù)隨機(jī)森林算法篩選出6 個(gè)用于估測(cè)油菜潛在產(chǎn)量的因子，篩選過程中發(fā)現(xiàn)生長(zhǎng)季最高溫度、1 月份最高溫度和2 月份太陽輻射均與潛在產(chǎn)量呈負(fù)相關(guān)關(guān)系，而生長(zhǎng)季降水量和7月份最低氣溫呈正相關(guān)關(guān)系。 這與前人研究結(jié)果基本一致，如張皓等[24]研究發(fā)現(xiàn)長(zhǎng)江流域冬油菜產(chǎn)量與蕾薹期和花期的溫度呈顯著負(fù)相關(guān)關(guān)系。值得注意的是，油菜在長(zhǎng)江流域的生長(zhǎng)季節(jié)為10月到次年5 月，而7 月份的平均最低氣溫在不同的機(jī)器學(xué)習(xí)算法中均被篩選為關(guān)鍵因子，相關(guān)分析結(jié)果也顯示7 月份的平均最低氣溫和油菜潛在產(chǎn)量的相關(guān)系數(shù)高達(dá)0.85。 進(jìn)一步分析發(fā)現(xiàn)，7月的平均最低氣溫與其他關(guān)鍵因子如生長(zhǎng)季降雨、最高氣溫、1 月最高氣溫和2 月的太陽輻射均有顯著的相關(guān)關(guān)系，相關(guān)系數(shù)分別為0.54、-0.80、-0.74 和-0.74，不同時(shí)間段氣象因子間的相互關(guān)系，可能與天氣系統(tǒng)的相互作用有關(guān)。

本研究結(jié)果顯示，長(zhǎng)江中下游的冬油菜優(yōu)勢(shì)產(chǎn)區(qū)主要分布于江蘇、安徽、湖北、重慶大部分地區(qū)及江西、湖南、貴州的北部地區(qū)，油菜多年平均潛在產(chǎn)量可達(dá)3 000 kg/hm2以上。 目前長(zhǎng)江中游二熟區(qū)和三熟區(qū)油菜種植面積最為廣泛，共計(jì)375.6 萬hm2，而長(zhǎng)江下游油菜種植面積相對(duì)較低，為66.13 萬hm2[23]。 長(zhǎng)江中下游冬閑田的遙感監(jiān)測(cè)顯示，當(dāng)前長(zhǎng)江中下游農(nóng)業(yè)區(qū)冬閑田面積仍占總耕地面積的45%左右，冬閑田主要集中于長(zhǎng)江中下游北部，其中湖北、湖南和江蘇是冬閑田分布面積最多的省份[10]。 結(jié)合本研究結(jié)果，建議在這些地區(qū)推廣冬油菜種植，充分利用冬閑田，提高油菜總產(chǎn)量。 同時(shí)，在湖南、江西、湖北和安徽的北部，當(dāng)前產(chǎn)量與潛在產(chǎn)量的產(chǎn)量差仍在2 000 kg/hm2，縣級(jí)統(tǒng)計(jì)產(chǎn)量?jī)H為潛在產(chǎn)量的50%左右，建議在這些地區(qū)進(jìn)一步開展油菜種植技術(shù)的推廣，如水氮管理、播期調(diào)控和病蟲害防控等[24-27]，以提高油菜產(chǎn)量。

4 結(jié)論

長(zhǎng)江流域是我國冬油菜種植的優(yōu)勢(shì)產(chǎn)區(qū)，油菜潛在產(chǎn)量變化范圍在1 000～4 000 kg/hm2之間，分布呈北高南低、西低東高的趨勢(shì)。 江蘇、安徽、湖北、重慶大部分地區(qū)及江西、湖南、貴州的北部地區(qū)，油菜多年平均潛在產(chǎn)量可達(dá)3 000 kg/hm2以上，是長(zhǎng)江流域的優(yōu)勢(shì)產(chǎn)區(qū)，應(yīng)重點(diǎn)推廣冬油菜種植，充分利用這些地區(qū)的冬季閑田。 與實(shí)際統(tǒng)計(jì)產(chǎn)量比較，油菜產(chǎn)量差在江西、湖南大部分地區(qū)及安徽和湖北北部可達(dá)2 000 kg/hm2，縣級(jí)統(tǒng)計(jì)產(chǎn)量?jī)H有潛在產(chǎn)量的50%左右，仍有很大的產(chǎn)量提升空間，在這些地區(qū)應(yīng)進(jìn)一步推廣高效的油菜種植技術(shù)，提高油菜的實(shí)際產(chǎn)量。 而在產(chǎn)量差較低的江蘇大部分地區(qū)、安徽和湖北北部、四川大部分地區(qū)、云貴高原地區(qū)，應(yīng)考慮通過選育適宜的高產(chǎn)品種進(jìn)一步提升產(chǎn)量潛力。