水電站入庫流量預測–基於自定義損失函數的循環神經網路建模方法

• 2020 年 11 月 5 日
• 筆記
• LSTM, 損失函數, 深度學習

從志在必得到鎩羽而歸——記一次大數據競賽經歷

最近參加了一個比賽，在工業大數據產業創新平台上，是一個水電站入庫流量預測問題。簡單看了一下題目，嚯，這個方向以前有做過啊，不說了~開整。

賽題背景：對進入水電站水庫的入庫流量進行精準預測，能夠幫助水電站對防洪、發電計劃調度工作進行合理安排。入庫流量受到降水、蒸發、土壤（直接影響地表徑流、地下徑流）、上游來水等諸多因素的綜合影響（如下圖所示）。通過數據分析方法，基於大量的歷史數據和可獲取的監測數據，實現水電站入庫流量的精準預測，將對水電站生產帶來顯著的安全和經濟價值。

賽題描述：基於歷史數據和當前觀測資訊，對電站未來7日入庫流量進行預測。預測時間間隔為3個小時，即對未來56個時間節點的入庫流量數據進行預測。其中，歷史數據包括：入庫流量數據、環境觀測數據，遙測站降雨數據，降雨預報數據。

一.問題分析及思路簡述

這是一個典型的時序數據預測問題。首先對數據進行分析。本次題目的原始數據時間段為2013/1/1——2018/1/7，共計5年的歷史數據，目標對2018年3段數據(初賽)，2019年5段數據(複賽)進行預測，各類數據情況如下。

1.入庫流量數據：這是本題的預測目標，每個數據點之間時間間隔為3小時。通過觀察可以看出，主辦方已經將這些數據進行了歸一化，流量數據都分布於0~1之間。通過對數據文件的時間節點進行掃描，發現缺失2014-11-20 11:00:00到2014-12-31 23:00:00之間的數據。

2.環境表數據：環境表的數據間隔為天，即每隔24小時一個數據點。數據文件中共有三類數據，包括溫度，風速，風向。其中溫度和風速已經進行了歸一化操作，數據依然分布在0~1之間，而風向則是原始值，數值較大且沒有明顯區分度（數值處於999001~999016之間）。

3.遙測站降雨數據：這個數據文件較大，其時間粒度最細，每隔1個小時即有一個數據點。另外每個數據點包含39個測站，每一個數據即代表某個測站在某一小時內的降雨量。可以看出，數據也已經進行過歸一化，但其數值一般都比較小，普遍小於0.1。

4.降雨預報數據：降雨預報的數據以天為單位，即時間間隔為24小時。每一個時間節點上，包含5個數據點，代表未來5天的降雨情況。數據也都已經被歸一化到0~1範圍之內。經過掃描發現有部分日期的降雨數據缺失。如下表：

表1 降雨數據缺失情況

經過對數據的簡單分析，發現數據存在下述問題：①數據缺失，建立模型之前，需要對缺失數據進行處理，一般使用填充法或者是在拼接訓練樣本時跳過缺失的數據段。②時間尺度不統一，預測目標入庫流量數據的時間間隔為3小時，遙測站降雨數據時間間隔為1小時，而環境表和降雨預報數據的時間間隔都是24小時。在建模之前需要對各類數據的時間尺度進行統一。③特徵較為繁雜，這裡主要指的是降雨量預報數據和遙測站降雨數據；其中每個時間點的數據都包含較多特徵，在建立模型時是否將其納入，都是我們後續要完成的工作。

由於數據時序性突出且在每個時間點包含多種特徵（即數據有包含多個時間節點，每個時間節點又包含多個特徵數據），因此計劃使用RNN模型，將多個時間節點的多條特徵數據丟入模型一把梭。其建模過程大概分為以下幾步：①數據預處理及特徵選擇，②模型建立及訓練，③結果預測。其中需要額外注意的是：本次比賽的評分規則為納什效率係數（NSE），因此在建立模型的過程中，需要對損失函數（誤差判別函數）進行特殊的設計。評分規則截圖如下：

 二.數據預處理及特徵選擇

數據預處理一般包含野值剔除、空缺值填充、數據歸一化等操作。由於此次賽會提供的數據基本上都經過了歸一化，就不必再進行額外的歸一化操作。又由於降雨量數值的躍動性較強，極其不穩定（如圖1），因此常規的空缺值填充方法（如均值填充、相似特徵回歸等）所得到的結果可信度較差，在數據預處理階段，便不對空缺值進行填充，僅在拼接模型時將其跳過。同樣的，由於降雨量數據強烈的不穩定性（如圖1），這裡不對其進行野值剔除操作，常規的3σ準則等野值判別方式也無法有效工作。

圖1 部分降雨量統計示意圖

在特徵選擇方面，第一步工作是對遙測站降雨數據進行降維，即對39個測站的數據進行加和，形成一個單一的降雨量特徵。並將其時標和入庫流量數據進行統一，變為每3小時1個數據點。第二步的工作是驗證降雨量預報工作的可信度，具體方法是使用協方差函數，計算降雨量序列和降雨量預報序列的相關係數，其結果如下：

第一天預報： 0.6411514582529538

第二天預報： 0.556989653623267

第三天預報： 0.5372333462465179

第四天預報： 0.4850024133929662

第五天預報： 0.4228036934449081

可以看出，相關係數都大於0，即表明未來5天的降雨預報與真實的降雨量之間都是正相關的，可以說預測結果都較為準確，因此將這5個特徵全部納入循環神經網路模型。另外基於對數據的直觀理解，在環境表中，風速和溫度可以直接影響水域的蒸發速率，而風向關係不是很大，因此將風向這個特徵進行移除。

三.模型建立及訓練

1.訓練模型拼接經過上一步的特徵選擇工作，現在一共剩下9種特徵數據，也就是說在一個時間節點，共有[溫度、風速、D1降雨、D2降雨、D3降雨、D4降雨、D5降雨、實際降雨量（近3小時）、入庫流量]9條數據。由於題目最終要基於前一個月的數據對下個月前7天的入庫流量進行預測，這裡直接將時間窗口拉滿，使用前30天的數據作為模型輸入，後7天的數據作為模型輸出。數據的時間間隔為3小時，30天的數據共有240個時間節點，7天的數據共有56個時間節點。最終的一個訓練樣本由2部分組成，第一部分是維度為(240*9)的輸入，第二部分使維度為(56*1)的輸出，輸出中僅包含接下來7天內的入庫流量特徵數據。如圖2所示，單個樣本的輸入為綠色框中的部分，包括240個時間節點的全部特徵數據；樣本的輸出為紅色框內的部分，即後續56個時間節點的入庫流量值。

圖2 單個樣本輸入及輸出示意

 在進行拼接的過程中，採用窗口滑動的方法逐點獲取樣本，同時注意要將所有特徵的時間節點對齊，所有特徵數據都是時間連續且一致的情況下，再將其存入訓練樣本集合。最終得到11189條訓練樣本。

圖3 訓練樣本及其維度

 2.LSTM循環神經網路構建

①前向傳播計算結構

這裡採用LSTM循環神經網路模型對數據進行建模，可以同時將數據的時序性和多特徵性納入考量。在每一個時間節點，LSTM單元的輸入由2部分組成——當前時間節點的輸入，以及上一個時間節點的輸出。

圖4 LSTM循環神經網路模型

上圖是LSTM單元的示意圖，其公式不再具體描述。訓練樣本中每一個時間節點的數據，都作為輸入的一部分進入LSTM單元，經計算後得到輸出並再次循環進入該LSTM單元；直到最後一個時間節點，獲得一個特徵向量，維度為(1*神經元數)，將它乘以一個維度為(神經元數*56)的矩陣，得到一個長度為56的向量，使用Sigmoid激活函數計算後，向量中的所有值都位於(0,1)之間。該向量即為模型的最終輸出結果，代表未來56個時間節點的入庫流量預測。

1 self.w = tf.Variable(tf.truncated_normal([self.num_nodes, self.output_dim], -0.1, 0.1),name='w')   # output_dim 是指輸出維度，即一個樣本的Y包含幾個值
2 self.b = tf.Variable(tf.zeros([self.output_dim]),name='b')
3 self.batch_in = tf.placeholder(tf.float32, [None, self.train_data.shape[1], self.input_dim],name='batch_in')
4 self.batch_out = tf.placeholder(tf.float32, [None, self.output_dim],name='batch_out')
5 lstm_cell = tf.nn.rnn_cell.BasicLSTMCell(self.num_nodes,forget_bias=1.0,state_is_tuple=True)
6 output, final_state = tf.nn.dynamic_rnn(lstm_cell, self.batch_in, time_major=False, dtype=tf.float32)
7 self.y_pre = tf.nn.sigmoid(tf.matmul(final_state[1], self.w) + self.b,name="y_pre")

上述程式碼中，w，b為最後輸出結果的權重參數和偏置，用來將RNN最後一個時間節點輸出的特徵向量轉換為一個長度為56的最終輸出向量。batch_in，batch_out則是輸入的佔位符，代表訓練樣本的輸入(X)和輸出(Y)，其第一維設置為None，使其可以匹配任意數量的樣本，方便成批開展訓練。在第5行中的num_nodes代表神經元數目，同時也是特徵向量的長度。y_pre則代表前向傳播的最終結果，也是最終得到的預測值。

②定義損失函數及反向傳播計算

之後需要定義損失函數及反向傳播計算。這裡重點對損失函數的設計進行詳細說明，本次題目的評分規則比較特殊，不同於一般的交叉熵函數、均方根誤差等評價方式，這裡使用納什效率係數(NSE)作為評價指標。

在進行NSE計算的過程中，將數據分成兩段，前2天的誤差結果給予0.65的權重，後5天的結果給予0.35的權重。根據題目要求，NSE越大，評分越高。也就是說需要公式右側的兩個減數：

 儘可能小。因此，在設計損失函數時，將nse’設置為損失函數，並在訓練模型的過程中使其朝著損失變小的方向進行訓練。

1 self.batch_out_mean = tf.reshape(tf.reduce_mean(self.batch_out,axis=1),[-1,1])
2 self.nse_left = 0.65*tf.reduce_sum(tf.square(self.batch_out[:,:16]-self.y_pre[:,:16]),axis=1)/tf.reduce_sum(tf.square(self.batch_out[:,:16]-self.batch_out_mean),axis=1)
3 self.nse_right = 0.35*tf.reduce_sum(tf.square(self.batch_out[:,16:]-self.y_pre[:,16:]),axis=1)/tf.reduce_sum(tf.square(self.batch_out[:,16:]-self.batch_out_mean),axis=1)
4 self.nse = tf.reduce_mean(self.nse_left+self.nse_right,name="nse")
5 self.train_op = tf.train.RMSPropOptimizer(self.lr).minimize(self.nse)
6 self.saver = tf.train.Saver()

上述程式碼中，self.nse代指的是上方公式中的nse’變數。首先計算這一批真實值中每一個樣本向量的均值，再分別計算nse’的兩個部分，賦予其0.65及0.35的權重，之後再加和。之後在第5行，設置優化器，使模型在反向傳播計算的過程中，沿著nse’減小的方向優化。反向傳播的過程構造完畢。

3.模型訓練及參數調優

①梯度爆炸問題的解決

在對部分訓練樣本進行訓練的過程中，發現某些樣本在經歷幾輪訓練後損失值急劇攀升，成為某個較大的固定值後便不再有任何變動。初步分析是損失函數nse’中，有部分Q_o^t值與觀測序列均值接近，導致分母部分接近於0，以致損失函數值激增，導致梯度爆炸。為解決這個問題，分階段對訓練樣本進行測試，篩除會導致梯度爆炸的樣本，最終篩除的樣本為：6200~6460，7650~8100共計710個樣本，約佔樣本總數的6%。

②訓練樣本切分

經過樣本篩除之後，剩餘樣本共10479個，將這些訓練樣本劃分為3個集合，即訓練集、驗證集、測試集。其中訓練集佔比90%，共9431條樣本；驗證集和測試集佔比都為5%，分別含有524條訓練樣本。

③學習率調節及訓練終止條件設置

在每一輪次的訓練中，僅使用訓練集中的數據對模型參數進行優化調節，並在完成一輪訓練之後，使用當前的模型參數，對驗證集中的訓練樣本進行預測，並計算其nse’。如果該模型對訓練集的預測效果優於上一輪，則繼續訓練；如果模型連續3輪在驗證集中沒有取得更好的預測效果，則將學習率降低50%，以防止模型在最優解附近震蕩；如果模型連續4輪沒有在驗證集中取得更好的預測效果，則結束訓練，保存模型參數。訓練結束以後，使用得到的模型對測試集數據進行預測，計算其nse’結果，與驗證集差距不大的話，則說明模型的預測結果可信。

 1 last_loss = 9999.9
 2 last_4_converge = [1,1,1,1]
 3 last_3_converge = [1,1,1]
 4 for i in range(epochs):
 5     for j in range(int(len(self.train_data)/self.batch_size)):
 6         batch_i = self.train_data[j*self.batch_size:(j+1)*self.batch_size]
 7         batch_o = self.train_label[j*self.batch_size:(j+1)*self.batch_size]
 8         self.sess.run(self.train_op, feed_dict={self.batch_in:batch_i, self.batch_out:batch_o.reshape(self.batch_size,self.output_dim)})
 9     if (i+1)%1==0:
10         this_loss = self.sess.run(self.nse,feed_dict={self.batch_in:self.validate_X, self.batch_out:self.validate_Y.reshape(len(self.validate_Y),self.output_dim)})
11         if(this_loss>=last_loss):
12             last_3_converge.append(0)
13             last_3_converge.pop(0)
14             last_4_converge.append(0)
15             last_4_converge.pop(0)
16         else:
17             last_3_converge.append(1)
18             last_3_converge.pop(0)
19             last_4_converge.append(1)
20             last_4_converge.pop(0)
21         last_loss = this_loss
22         train_loss = self.sess.run(self.nse,feed_dict={self.batch_in:batch_i,self.batch_out:batch_o.reshape(self.batch_size,self.output_dim)})
23         print('epoch:%d train:'%(i+1),train_loss,'  validate:',this_loss,' ',last_4_converge)
24         if((1 not in last_4_converge) or (self.lr<0.000002)):  #連續4輪未收斂或學習率過小時，結束訓練
25             break
26         elif(1 not in last_3_converge):  #連續3輪未收斂，降低學習率
27             self.lr = self.lr*0.5

④模型訓練

根據上述的方法描述，構建LSTM深度學習模型，開展訓練，經過多輪調整測試，發現超參數在設置為下表中的結果時，預測效果較好。

表2 模型超參數設置建議

通過共計1000輪次的訓練(兩輪500次的訓練)，最終得到預測模型。其對訓練集上最後一個批次的數據預測偏差(nse’)為0.151，對驗證集中所有數據的預測偏差(nse’)為0.19552，對測試集中所有樣本的預測偏差(nse’)為0.19509。

 圖5 訓練結果

結果出來一看，哎，這模型厲害了！對於沒有訓練過的的測試集來說，其預測結果的平均nse’是0.19509，按照賽會給出的計分公式: score=(1-nse’)*100，理論上講，對於新的輸入樣本，預測結果的nse平均應該是0.805左右，咋地也能得到80來分的成績。後面又把測試集的預測結果影像畫出來一看，嘿，還真不錯。

 上面是測試集中找了幾個比較有代表性的數據樣本預測結果，不論是哪種形狀的數據，遞增的、大段固定值的、以天為單位強周期的，都能取得不錯的效果，尤其是對入庫流量的峰谷趨勢把握的相當到位。預測最差的也就是最後一張圖片，nse’的值約為1.064。看到這些，瞬間覺得信心爆棚，勝利似乎近在眼前。

但現實總是無情，在使用這個模型對目標數據進行預測時，效果很不理想。上傳預測結果後，基本上是負分，對於複賽中2019年的5段數據更是超出了下界，即意味著nse’的值超過了2.0。最後終究沒有衝進決賽圈（前30），折戟賽場，鎩羽而歸。

四.原因分析及個人感悟

技術層面：

對於已有的2013年初至2017年底的數據，本文所構建的LSTM循環神經網路表現良好，能夠對未來7天的入庫流量數據進行精準的預測。但對主辦方給出的3段2018年數據及5段2019年數據，則效果欠佳，損失函數nse’值擴大了5~20倍不等。由於主辦方並未公開2018年及2019年數據，無法做到按例分析。猜測出現該問題的原因可能是以下幾個方面：

①模型對已有數據(2013.1——2017.12)產生了過擬合。理論上來講，這種現象是不會發生的，因為在進行模型訓練的過程中，模型參數調整僅依賴於訓練集，而模型訓練程度則由驗證集進行判斷，測試集始終沒有參與其中。以前在別的時序數據預測項目中出現過模型在測試集上表現糟糕，而在訓練集和驗證集上表現良好的現象，可以推測出模型對訓練集合驗證集擬合過度。本次比賽中這種對訓練集、驗證集、測試集都表現良好而在新的樣本上表現差勁的現象還是第一次遇到。經過思考，為避免模型對已有數據產生過擬合，可以在設置終止條件時，控制學習率的大小，為其規定一個最小值，當學習率調整小於該值時，停止訓練。

②在解決梯度爆炸問題的過程中，剔除了710個樣本，這710個樣本中可能存在某種規律沒有被模型所掌握。關於這個問題有兩種解決方案：一是在損失函數中，在計算分母時判斷其值的大小，為其規定一個最小值，如0.01，避免梯度的爆炸；二是使用常規損失函數如rmse(均方根誤差)來訓練出第二個模型，預測結果由兩個模型進行均值融合後獲得。

③時間窗口設置過長。這次為了節約建模時間，把時間窗口設置到最大的30天，這樣便增加了模型的複雜度，在進行訓練時，更多的時間窗口便意味著特徵向量在LSTM結構里循環了更多的次數，也就更容易產生梯度爆炸/消失問題。另外，直觀上來講，30天前的降雨對當前的河流水量產生不了太大影響，因此設定一個較短的、符合實際情況的時間窗口更加合適。

心理層面：

通過這次比賽有以下收穫：一是不能鑽牛角尖，驕兵必敗。在初賽階段，我已經發現上傳的預測結果無法達到模型的理論得分，當時經過了500輪的訓練，模型的理論得分在70分左右，但上傳的2018年3條數據預測結果只得到了-20分的成績。由於對於自己的模型信仰堅定，我判斷是訓練還不夠不充分所導致。便又深入訓練了500輪，達到了理論得分80分的水平，在複賽中對2019年5條數據的預測結果直接崩盤，達到-100分開外……二是將所有特徵放進模型裡面一把梭，會使其參數量增大，模型複雜度也會跟著暴漲，根據奧卡姆剃刀法則，越複雜的模型也就越難以信任。未來的建模過程中還應投入更多的人工分析工作，不能一味信任複雜模型。

歡迎各位大佬進行指導，提出改進意見。