勾志竟，宮志宏，劉布春

(1.天津市氣象信息中心，天津 300074;2.天津市氣候中心，天津 300074;3.中國(guó)農(nóng)業(yè)科學(xué)院農(nóng)業(yè)環(huán)境與可持續(xù)發(fā)展研究所，北京 100081)

0 引 言

水稻是中國(guó)最主要的糧食作物之一，其產(chǎn)量的預(yù)測(cè)對(duì)于農(nóng)業(yè)政策調(diào)整和保障糧食安全問(wèn)題具有重要的意義[1]。由于其受多種因素影響，水稻產(chǎn)量的準(zhǔn)確預(yù)測(cè)較為困難。常規(guī)的產(chǎn)量預(yù)測(cè)主要包括基于統(tǒng)計(jì)回歸原理預(yù)測(cè)、灰色理論生長(zhǎng)模擬[2]、衛(wèi)星遙感預(yù)測(cè)[3]等方法，這些方法簡(jiǎn)便易行，可以預(yù)測(cè)產(chǎn)量。然而，水稻產(chǎn)量與多種因素之間存在著復(fù)雜的非線性關(guān)系，傳統(tǒng)的方法預(yù)測(cè)精度往往不高，隨著基于人工智能方法的不斷研究和優(yōu)化[4-5]，BP神經(jīng)網(wǎng)絡(luò)因其適用于分析產(chǎn)量與影響產(chǎn)量的眾多因素之間的復(fù)雜非線性關(guān)系，在產(chǎn)量預(yù)測(cè)方面得到較為廣泛的應(yīng)用[6]。

張成才等[7]采用播種面積、農(nóng)田有效灌溉面積、農(nóng)用機(jī)械總動(dòng)力、化肥施用量、農(nóng)藥施用量等指標(biāo)建立了BP神經(jīng)網(wǎng)絡(luò)糧食產(chǎn)量預(yù)測(cè)模型；姜新等[8]基于作物參數(shù)和農(nóng)機(jī)水平建立了神經(jīng)網(wǎng)絡(luò)產(chǎn)量預(yù)測(cè)模型；郭亞菲等[9]提出了基于主成分分析和粒子群優(yōu)化神經(jīng)網(wǎng)絡(luò)的糧食產(chǎn)量預(yù)測(cè)模型。但是這些模型對(duì)水稻產(chǎn)量的長(zhǎng)期變化趨勢(shì)預(yù)測(cè)能力偏弱，同時(shí)未能充分考慮生長(zhǎng)期的氣象要素對(duì)產(chǎn)量的影響。

而長(zhǎng)短期記憶網(wǎng)絡(luò)(long-short term memory，LSTM)作為一種改進(jìn)的時(shí)間循環(huán)神經(jīng)網(wǎng)絡(luò)，在長(zhǎng)時(shí)間序列預(yù)測(cè)問(wèn)題處理中有獨(dú)特的優(yōu)勢(shì)，在機(jī)器翻譯、語(yǔ)音識(shí)別等方面已經(jīng)有很多應(yīng)用[10-11]，但目前鮮有在水稻產(chǎn)量預(yù)測(cè)中的應(yīng)用。因此，該文擬基于TensorFlow深度學(xué)習(xí)框架，構(gòu)建LSTM水稻產(chǎn)量預(yù)測(cè)方法，以水稻不同生長(zhǎng)期的氣象要素作為預(yù)測(cè)因子對(duì)水稻產(chǎn)量進(jìn)行預(yù)測(cè)，并將預(yù)測(cè)結(jié)果與BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型進(jìn)行比較，由此探討LSTM在水稻產(chǎn)量長(zhǎng)期趨勢(shì)預(yù)測(cè)中的優(yōu)越性，以期可為水稻產(chǎn)量長(zhǎng)期趨勢(shì)預(yù)測(cè)提供一種新方法。

1 資料與方法

1.1 研究區(qū)域概況

天津市寧河區(qū)位于39°21′N，117°49′，屬于大陸性季風(fēng)氣候，暖溫帶半干旱半濕潤(rùn)風(fēng)帶，年均降水量為500 mm ～700 mm，降水量70%集中在6-8月份，地貌屬于海積、沖積平原，土壤含鹽分較高[12]。天津?qū)幒訁^(qū)“小站稻”既有食味特性，又有良好的耐鹽堿特性[13]，一般在4月上旬播種，5月底至6月上旬移栽，6月中下旬進(jìn)入返青期，8月上旬進(jìn)入孕穗期，9月上旬進(jìn)入乳熟期，9月下旬至10月中旬成熟。

1.2 數(shù)據(jù)來(lái)源

水稻產(chǎn)量數(shù)據(jù)和生育期地面氣象觀測(cè)資料來(lái)自天津市氣象信息中心，水稻產(chǎn)量數(shù)據(jù)為寧河區(qū)單位面積(每畝)水稻產(chǎn)量(kg)，生育期包括移栽-返青期(周期1)、分蘗期(周期2)、孕穗期(周期3)、抽穗期(周期4)、成熟期(周期5)。

地面氣象觀測(cè)資料來(lái)自天津市寧河區(qū)氣象觀測(cè)站，包括1989-2015年逐日平均風(fēng)速、最大風(fēng)速、日照時(shí)數(shù)、最高溫、最低溫、平均相對(duì)濕度、最小相對(duì)濕度及降水量。

1.3 研究方法

LSTM神經(jīng)網(wǎng)絡(luò)最早是由Hochreiter等[14]提出，是一種改進(jìn)的RNN算法，增加了忘記門(mén)，可以學(xué)習(xí)時(shí)間序列長(zhǎng)短依賴信息，能夠克服梯度消失的問(wèn)題，主要由忘記門(mén)、輸入門(mén)、候選門(mén)和輸出門(mén)組成，其基本結(jié)構(gòu)如圖1所示。

圖1 LSTM內(nèi)部結(jié)構(gòu)

忘記門(mén)(forget gate layer)控制著應(yīng)該忘記哪些信息，其計(jì)算公式如下所示：

ft=σ(Wf[ht-1,xt]+bf)

(1)

其中，Wf為忘記門(mén)權(quán)重，ht-1為上一個(gè)單元的輸出，xt為當(dāng)前單元的輸入，bf為偏置，σ為sigmod函數(shù)。

輸入門(mén)(input gate layer)通過(guò)激活函數(shù)決定著當(dāng)前時(shí)刻的輸入信息xt，有多少信息被加到信息流里面，其計(jì)算公式如下所示：

it=σ(Wi[ht-1,xt]+bi)

(2)

gt=σ(Wg[ht-1,xt]+bg)

(3)

其中，Wi、Wg為權(quán)重，bi、bg為偏置。

候選門(mén)(candidate layer)用來(lái)計(jì)算當(dāng)前的輸入和過(guò)去的記憶所具有的信息的總和，其計(jì)算公式如下所示：

Ot=σ(Wo[ht-1,xt]+bo)

(4)

ht=ot*tanh(St)

(5)

輸出門(mén)(output gate layer)控制著有多少記憶信息被用于下一階段的更新中，其計(jì)算公式如下所示：

St=ft*St-1+gt*it

(6)

本研究采用LSTM算法實(shí)現(xiàn)水稻產(chǎn)量預(yù)測(cè)模型，模型輸入為天津市寧河區(qū)水稻5個(gè)生長(zhǎng)周期的氣象要素觀測(cè)數(shù)據(jù)，即移栽-返青期平均風(fēng)速、移栽-返青期最大風(fēng)速、移栽-返青期平均風(fēng)速日照時(shí)數(shù)、移栽-返青期最高溫等，模型輸出為天津市寧河區(qū)單位面積的水稻產(chǎn)量。

1.4 模型評(píng)價(jià)指標(biāo)

本研究采用均方根誤差和平均絕對(duì)誤差來(lái)評(píng)價(jià)預(yù)測(cè)模型的優(yōu)劣程度。

(1)均方根誤差(root mean square error，RMSE)。

(7)

(2)平均絕對(duì)誤差(mean absolute error，MAE)。

(8)

2 結(jié)果與分析

2.1 1989-2015年水稻單位面積產(chǎn)量趨勢(shì)

寧河區(qū)1989-2015年的水稻單產(chǎn)趨勢(shì)如圖2所示，可以看出寧河區(qū)1993年水稻單產(chǎn)明顯低于往年(按缺測(cè)處理)，2010年之前整體波動(dòng)不大，2012年水稻單產(chǎn)最高，2014年和2015年相比也出現(xiàn)了較大程度的減產(chǎn)。

圖2 寧河區(qū)1989-2015年水稻單產(chǎn)變化趨勢(shì)

2.2 1989-2015年水稻各生長(zhǎng)期氣象要素變化特征

為了消除原始?xì)庀笠財(cái)?shù)據(jù)量綱不同所帶來(lái)的不利，通常的做法是歸一化處理，將各氣象要素?cái)?shù)據(jù)轉(zhuǎn)換到[0,1]之間歸一化公式如下：

(9)

其中，xi為輸入數(shù)據(jù)，i為數(shù)據(jù)序號(hào)，xmax、xmin為x中的最大值和最小值。

本研究主要考慮天津市寧河區(qū)水稻5個(gè)生長(zhǎng)周期的氣象要素對(duì)水稻產(chǎn)量的影響，通過(guò)對(duì)1989-2015年寧河區(qū)氣象要素觀測(cè)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，得到如下結(jié)論：

經(jīng)分析，1989-2015年寧河區(qū)水稻移栽-返青期平均風(fēng)速、最大風(fēng)速、平均氣溫、最高氣溫幾乎無(wú)變化，日照總時(shí)數(shù)減少，最低氣溫降低，最小相對(duì)濕度和平均相對(duì)濕度變小，降雨量減少。

分蘗期平均風(fēng)速變大，最大風(fēng)速減小，日照總時(shí)數(shù)增多，平均氣溫、最低氣溫、降雨量幾乎無(wú)變化，最高氣溫升高，最小相對(duì)濕度、平均相對(duì)濕度變小。

孕穗期平均風(fēng)速、最大風(fēng)速變大，日照總時(shí)數(shù)減少，平均氣溫、最高氣溫、最低氣溫升高，最小相對(duì)濕度幾乎無(wú)變化，平均相對(duì)濕度變小。

抽穗期平均風(fēng)速變大，最大風(fēng)速幾乎無(wú)變化，日照總時(shí)數(shù)減少，平均氣溫、最高氣溫、最低氣溫降低，最小相對(duì)濕度、平均相對(duì)濕度變小，降雨量減少。

成熟期平均風(fēng)速、最大風(fēng)速變大，日照總時(shí)數(shù)減少，平均氣溫、最高氣溫、最低氣溫升高，最小相對(duì)濕度、平均相對(duì)濕度變小，降雨量幾乎無(wú)變化。

2.3 預(yù)測(cè)因子的選擇

預(yù)測(cè)因子的選擇很大程度上影響了預(yù)測(cè)的結(jié)果，該文利用Pearson相關(guān)系數(shù)法[15]來(lái)考察各氣象要素對(duì)水稻產(chǎn)量的影響力，公式如下：

(10)

其中，x與y分別為2個(gè)變量的觀測(cè)值。

若r>0，表示2個(gè)變量是正相關(guān)的；若r<0，則表示2個(gè)變量是負(fù)相關(guān)的，n為樣本數(shù)量(文中為26)。

各生長(zhǎng)周期各要素與水稻產(chǎn)量相關(guān)性系數(shù)如表1所示，其中，X1、X2、X3、X4、X5、X6、X7、X8、X9分別表示平均風(fēng)速、最大風(fēng)速、日照總時(shí)數(shù)、平均溫度、最高溫、最低溫、平均相對(duì)濕度、最小相對(duì)濕度及總降水量。

表1 各生長(zhǎng)周期各要素與水稻產(chǎn)量相關(guān)性系數(shù)

變量移栽-返青期分蘗期孕穗期抽穗期成熟期X1-0.139-0.109-0.338-0.341-0.439X2-0.239-0.266-0.191-0.388-0.496X30.1040.0940.034-0.1130.29X4-0.237-0.0690.2780.0920.103X5-0.046-0.1170.2410.0380.052X6-0.2720.1320.350.1160.148X7-0.0970.2050.320.0610.087X8-0.2250.2630.2230.0180.032X9-0.2460.4440.1870.0440.241

(11)

最后經(jīng)過(guò)顯著性檢驗(yàn)，如公式(11)所示，并結(jié)合農(nóng)業(yè)氣象專家的建議，選取移栽-返青期最大風(fēng)速、移栽-返青期平均溫度、移栽-返青期最低溫、移栽-返青期最小相對(duì)濕度、移栽-返青期總降水量、分蘗期最大風(fēng)速、分蘗期平均相對(duì)濕度、分蘗期最小相對(duì)濕度、分蘗期總降水量、孕穗期平均風(fēng)速、孕穗期平均溫度、孕穗期最高溫、孕穗期最低溫、孕穗期平均相對(duì)濕度、孕穗期最小相對(duì)濕度、抽穗期平均風(fēng)速、抽穗期最大風(fēng)速、成熟期平均風(fēng)速、成熟期最大風(fēng)速、成熟期日照總時(shí)數(shù)、成熟期總降水量這21個(gè)要素作為水稻產(chǎn)量的預(yù)測(cè)因子。

2.4 TensorFlow框架

TensorFlow是谷歌的開(kāi)源機(jī)器學(xué)習(xí)算法框架[16]，其應(yīng)用場(chǎng)景非常廣泛，包括語(yǔ)音識(shí)別、自然語(yǔ)言處理、計(jì)算機(jī)視覺(jué)、數(shù)據(jù)分析及預(yù)測(cè)等[17]，它支持Python、C++等多種編程語(yǔ)言，計(jì)算模型采用有向圖(directed graph)，其中每個(gè)運(yùn)算操作(operation)作為一個(gè)節(jié)點(diǎn)(node)，節(jié)點(diǎn)間的連接為邊(edge)。

本研究采用Anaconda在Windows操作系統(tǒng)下安裝CPU版本的TensorFlow，Anaconda提供了Python的包管理和環(huán)境管理功能。

TensorFlow的具體安裝步驟如下：

(1)到Anaconda官網(wǎng)上下載電腦系統(tǒng)對(duì)應(yīng)版本的Anaconda并安裝。

(2)創(chuàng)建運(yùn)行環(huán)境，安裝Python3.6.10，命令如下：

conda create --name tensorflow python=3.6.10

(3)激活TensorFlow環(huán)境，命令如下：

activate tensorflow

(4)安裝TensorFlow，命令如下：

pip install tensorflow

(5)啟動(dòng)Anaconda，并安裝Spyder。

2.5 基于TensorFlow的LSTM預(yù)測(cè)模型

實(shí)驗(yàn)采用Python語(yǔ)言和TensorFlow庫(kù)實(shí)現(xiàn)LSTM預(yù)測(cè)模型，如圖3所示。

圖3 LSTM水稻產(chǎn)量預(yù)測(cè)模型

(1)輸入層：以天津市寧河區(qū)水稻生長(zhǎng)周期的21個(gè)氣象要素作為訓(xùn)練數(shù)據(jù)，水稻單產(chǎn)數(shù)據(jù)作為目標(biāo)數(shù)據(jù)，對(duì)原始?xì)庀笠睾退井a(chǎn)量數(shù)據(jù)進(jìn)行預(yù)處理，得到逐年的時(shí)間序列數(shù)據(jù)集，并將數(shù)據(jù)集劃分為訓(xùn)練集Tran_data和測(cè)試集Test_data，采用公式(9)對(duì)訓(xùn)練集和測(cè)試集進(jìn)行標(biāo)準(zhǔn)化處理。

訓(xùn)練數(shù)據(jù)和目標(biāo)數(shù)據(jù)實(shí)現(xiàn)如下：

X=tf.placeholder(tf.float32,[None,time_step,input_size])

Y=tf.placeholder(tf.float32,[None,time_step,output_size])

其中，X為訓(xùn)練數(shù)據(jù)，Y為目標(biāo)數(shù)據(jù)。

(2)隱藏層：將輸入層處理后的數(shù)據(jù)傳入隱藏層，計(jì)算相對(duì)損失，對(duì)LSTM算法隱藏層參數(shù)layer_num、rnn_unit、weights、biases進(jìn)行優(yōu)化，將最終結(jié)果輸出。

其中，layer_num為隱藏層數(shù)量，rnn_unit為隱藏層單元數(shù)。輸入層和輸出層的權(quán)重為weights，偏置為biases。

輸入層和輸出層權(quán)重、偏置及誤差實(shí)現(xiàn)如下：

weights={'in':tf.Variable(tf.random.normal([input_size,rnn_unit])),'out':tf.Variable(tf.random.normal([rnn_unit,1]}

biases={'in':tf.Variable(tf.constant(0.1,shape=[rnn_unit,])),'out':tf.Variable(tf.constant(0.1,shape=[1,]))}

loss=tf.reduce_mean(tf.square(tf.reshape(pred,[-1])-tf.reshape(Y, [-1])))

(3)輸出層：通過(guò)隱藏層的輸出計(jì)算得到均方根誤差和平均絕對(duì)誤差，本研究網(wǎng)絡(luò)單元采用tf..nn.rnn_cell. BasicLSTMCell，默認(rèn)激活函數(shù)為tanh。

(4)優(yōu)化：神經(jīng)網(wǎng)絡(luò)最常用的優(yōu)化方法多為基于梯度的，本研究針對(duì)傳統(tǒng)神經(jīng)網(wǎng)絡(luò)梯度下降問(wèn)題采用自適應(yīng)距估計(jì)優(yōu)化器(Adam)，該算法將Momentum和RMSprop融合在一起，能夠根據(jù)參數(shù)調(diào)節(jié)相應(yīng)的學(xué)習(xí)率，并且內(nèi)存占用較少。同時(shí)，本研究采用Dropout方法防止網(wǎng)絡(luò)在訓(xùn)練過(guò)程中出現(xiàn)過(guò)擬合問(wèn)題。

Adam優(yōu)化器和Dropout方法實(shí)現(xiàn)如下：

train_op=tf.AdamOptimizer(lr).minimize(loss)

cell=tf.nn.rnn_cell.DropoutWrapper(cell,output_keep_prob=keep_prob)

(5)預(yù)測(cè)：通過(guò)saver.restore(sess,module_file)將保存的訓(xùn)練好的最優(yōu)模型參數(shù)恢復(fù)，輸入新的預(yù)測(cè)因子，即可預(yù)測(cè)水稻產(chǎn)量。

2.6 實(shí)驗(yàn)結(jié)果

文中基于TensorFlow構(gòu)建了LSTM長(zhǎng)短時(shí)記憶網(wǎng)絡(luò)，并通過(guò)與傳統(tǒng)的BP神經(jīng)網(wǎng)絡(luò)進(jìn)行對(duì)比實(shí)驗(yàn)，驗(yàn)證基于TensorFlow構(gòu)建的深度學(xué)習(xí)模型的優(yōu)越性。網(wǎng)絡(luò)參數(shù)及調(diào)優(yōu)方法參考文獻(xiàn)[18-19]，通過(guò)多次實(shí)驗(yàn)對(duì)模型進(jìn)行調(diào)優(yōu)得到BP神經(jīng)網(wǎng)絡(luò)和LSTM長(zhǎng)短時(shí)記憶網(wǎng)絡(luò)的參數(shù)，如表2所示。

表2 BP神經(jīng)網(wǎng)絡(luò)和LSTM參數(shù)

從表2可以看出，BP神經(jīng)網(wǎng)絡(luò)的網(wǎng)絡(luò)結(jié)構(gòu)為21-q-l，其中21為輸入層的節(jié)點(diǎn)數(shù)目，q為隱含層數(shù)目，采用經(jīng)驗(yàn)公式q=2*d+1，d為輸入層節(jié)點(diǎn)數(shù)目，l為輸出層節(jié)點(diǎn)數(shù)目，即水稻單位面積產(chǎn)量，學(xué)習(xí)率設(shè)為0.005，迭代次數(shù)設(shè)為1 000，誤差精度設(shè)為0.005 1。LSTM長(zhǎng)短期記憶網(wǎng)絡(luò)的輸入層節(jié)點(diǎn)數(shù)為21，輸出層節(jié)點(diǎn)數(shù)為1，隱藏層數(shù)設(shè)為3、隱藏層節(jié)點(diǎn)數(shù)為100，學(xué)習(xí)率為0.000 6，迭代次數(shù)為300。

實(shí)驗(yàn)結(jié)果如圖4所示?？梢钥闯?，2014年相比2013年產(chǎn)量有減小的趨勢(shì)，BP神經(jīng)網(wǎng)絡(luò)沒(méi)有擬合出來(lái)。相比BP神經(jīng)網(wǎng)絡(luò)，LSTM長(zhǎng)短期記憶網(wǎng)絡(luò)通過(guò)門(mén)單元(gate)可以更好地學(xué)習(xí)到長(zhǎng)時(shí)間序列數(shù)據(jù)的細(xì)節(jié)，進(jìn)而更好地預(yù)測(cè)出未來(lái)的單位面積水稻產(chǎn)量的發(fā)展趨勢(shì)。

圖4 BP神經(jīng)網(wǎng)絡(luò)與LSTM預(yù)測(cè)結(jié)果對(duì)比

兩種方法的預(yù)測(cè)結(jié)果如表3所示。

表3 BP神經(jīng)網(wǎng)絡(luò)和LSTM產(chǎn)量預(yù)測(cè)結(jié)果

從表3可以看出，除了2014年，LSTM長(zhǎng)短時(shí)記憶網(wǎng)絡(luò)對(duì)2010-2015年的寧河區(qū)小站稻畝產(chǎn)量預(yù)測(cè)精度都比BP神經(jīng)網(wǎng)絡(luò)要高。

其中，BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)結(jié)果的RMSE為75.12，MAE為65.64。LSTM長(zhǎng)短期記憶網(wǎng)絡(luò)預(yù)測(cè)結(jié)果的RMSE為34.77，MAE為33.37，相比BP神經(jīng)網(wǎng)絡(luò)RMSE和MAE都有明顯的下降，RMSE降低了53.7%，MAE降低了49.2%，預(yù)測(cè)精度更高。

3 結(jié)束語(yǔ)

該文利用水稻產(chǎn)量的時(shí)間序列數(shù)據(jù)進(jìn)行預(yù)測(cè)，基于TensorFlow深度學(xué)習(xí)框架構(gòu)建了LSTM水稻產(chǎn)量預(yù)測(cè)方法，采用1989-2015年天津市寧河區(qū)地面氣象觀測(cè)資料與產(chǎn)量數(shù)據(jù)，選取水稻不同生長(zhǎng)期的氣象要素作為預(yù)測(cè)因子，對(duì)水稻單產(chǎn)進(jìn)行預(yù)測(cè)。結(jié)果表明，BP神經(jīng)網(wǎng)絡(luò)由于過(guò)擬合問(wèn)題，預(yù)測(cè)效果不佳，LSTM長(zhǎng)短期記憶網(wǎng)絡(luò)的預(yù)測(cè)精度較高，可以更好地預(yù)測(cè)未來(lái)的糧食產(chǎn)量的長(zhǎng)期發(fā)展趨勢(shì)，可為水稻產(chǎn)量預(yù)測(cè)提供新的參考思路。水稻產(chǎn)量受經(jīng)濟(jì)和自然條件等多種因素影響，準(zhǔn)確的預(yù)測(cè)往往十分困難，文中的方法僅考慮了氣象要素對(duì)產(chǎn)量的影響，未來(lái)可以綜合考慮社會(huì)經(jīng)濟(jì)因素及生產(chǎn)投入、品種等因素做進(jìn)一步的產(chǎn)量預(yù)測(cè)研究。