基于小波分析的汛期月徑流量預測模型研究

周 婭 ，郭 萍 ，楊 柳 ，宋培培

（1.貴州省水利水電勘測設計研究院，貴州 貴陽 550002；2.中國農業大學水利與土木工程學院，北京 100083）

崗南水庫位于河北省石家莊市平山縣境內，位于滹沱河干流上，是一座以防洪為主，兼顧灌溉、發電、工業及城市生活用水等綜合利用的大（I）型水庫。水庫承擔著石家莊和下游鐵路、公路、華北油田及冀中平原的防洪保護任務，并保障區域農業、工業和城市用水安全[1]。小覺水文站是崗南水庫入庫徑流控制水文站，上游來水相對天然，且受水庫調節作用影響較小，長系列水文數據代表性較好。本文首先以主成分分析法提取出汛期氣象因子的主成分，然后通過A Trous小波提取出各序列的隨機項、趨勢項和周期項，分別對各項建立BP神經網絡預測模型。

1 研究區概況

滹沱河屬于海河流域的子牙河水系，全長587 km，流域面積24690 km2。河流流經山西省境內的繁峙等五縣市，于陽泉市盂縣閆家莊進入河北省，然后由石家莊市區北穿京廣鐵路，向東橫貫河北省，并從衡水市的安平縣一直向東流至獻縣臧橋，在此和滏陽河相匯流入大海[1]（見圖1）。流域年氣溫為4.0℃～8.8℃，年內降水主要集中在6月～9月，約占全年降水量的70%～80%，上游位于太行山的背風坡，雨量較小，徑流量也較低，年徑流深僅為20 mm～50 mm左右，而經太行山的東坡進入河北省時，降雨量增加，徑流量也隨之驟增，年徑流深在50 mm～150 mm之間，然后經黃壁莊水庫進入山麓沖積平原，年徑流再次減少，徑流深降至50 mm以下。

圖1 滹沱河流域范圍

2 方法簡介

2.1 主成分分析法

1933年，Hotelling首先提出主成分分析，該方法利用降維的思想，以信息損失量很少為前提，將多個指標轉換為幾個綜合指標，由原始變量的線性組合得到每一個綜合指標，每個綜合指標之間互不相關[2]。

主成分分析的數學模型為

可用矩陣表示為

經過標準化處理的矩陣就是X的相關系數矩陣R，如果主成分分析的一切計算不是從協方差矩陣出發，而是直接從樣本相關系數矩陣出發，就等于先對數據進行標準化，然后再利用協方差矩陣進行主成分分析[2]。主成分分析的基本步驟是：首先結合數據，判斷是否需要進行主成分分析；其次結合主成分的特征值和累計特征值貢獻率來確定抽取的主成分或因子的數據；最后，進行主成分分析，將抽取出來的主成分分析存為新變量，以便繼續分析[3]。

2.2 小波變換—A Trous

對水文時間序列f（t），t=1，2，…，N進行小波分解，令C0（t）=f（t），A Trous小波分解過程如下：

A Trous小波重構過程如下：

根據樣本序列長度，低通濾波器選用db4小波，選擇最大分解尺度P=3，采用MATLAB小波分解工具箱分解得到樣本的隨機項、趨勢項和周期項，采用BP神經網絡對各項分別進行預測，再重構得到預測模型的輸出值[4]。

3 汛期月徑流預測

汛期采用1969年～2009年的7月、8月、9月、10月的16個氣象因子：x1五臺山降水量、x2五臺山風速，x3五臺山氣壓、x4五臺山氣溫、x5五臺山相對濕度、x6五臺山日照時數、x7五臺山水汽壓、x8原平降水量、x9原平風速、x10原平氣壓、x11原平氣溫、x12原平相對濕度、x13原平日照時數、x14原平水汽壓、x15小覺站降水量、x16小覺站蒸發量，以及小覺站徑流量值作為分析基礎。

因子分析前，首先進行KMO檢驗和巴特利球體檢驗。KMO檢驗用于檢查變量間的相關性和偏相關性，取值在0～1之間。KMO統計量越接近于1，變量間的相關性越強，因子分析的效果越好。實際分析中，KMO統計量在0.7以上時效果比較好；當KMO統計量在0.5以下時不適合應用因子分析法。采用主成分分析方法分析以上數據，得出汛期實測數據對應的KMO值為0.819（大于0.8），意味著變量間的相關性較強，說明變量很適合作因子分析。巴特利特球形檢驗法是以相關系數矩陣為基礎的。它的零假設相關系數矩陣是一個單位陣，即相關系數矩陣對角線的所有元素均為1，所有非對角線上的元素均為0。巴特利特球形檢驗法的統計量Sig值是根據相關系數矩陣的行列式得到的。如果Sig值較大，且其對應的相伴概率值小于指定的顯著水平0.05時，拒絕零假設，表明相關系數矩陣不是單位陣，原有變量之間存在相關性，適合進行主成分分析；反之，數據不適合進行主成分分析。對本次采用的數據進行Bartlett的球形度檢驗，其Sig值小于顯著水平0.05，說明變量之間存在相關關系。所選取的影響因子的數據符合因子分析的條件，適合做主成分分析。進行主成分分析時，由表1可知前4個主成分貢獻率累計百分比為88.943%，對應的初始特征值分別是 8.364、2.975、2.190、0.702，根據初始特征值計算得到主成分表達式的系數矩陣（表2），進一步計算得到汛期四個主成分 Y1、Y2、Y3、Y4作為預測輸入。

表1 解釋的總方差（汛期）

表2 主成分表達式的系數矩陣（汛期）

由于小波分析具有將時間序列分解為多層的能力，可分離季節性趨勢時間序列中的趨勢項、周期項和隨機項，分離后各項所含頻率成分單一，有利于建模預測。本研究將主成分分析提取得到的四個主成分Y1、Y2、Y3、Y4和小覺站月徑流值Z進行A Trous小波變換分解，提取得到各數據序列的隨機項（Y1D、Y2D、Y3D、Y4D、ZD），趨勢項（Y1A、Y2A、Y3A、Y4A、ZA），周期項（Y1P、Y2P、Y3P、Y4P、ZP）。在此基礎上，采用MATLAB自帶BP神經網絡訓練工具箱（ntstool）建立各項預測模型，將以上數據作為輸入信號從輸入層進入網絡，經隱含層逐層處理后到達輸出層。每層的神經元狀態只影響下一層的神經元狀態。如果輸出層得不到期望輸出，則誤差反向傳播，根據誤差調整網絡權值和閾值，從而使BP神經網絡預測輸出不斷逼近目標值，從而建立小覺站汛期月徑流預測模型，模型采用的是訓練較快的trainscg函數。

3.1 隨機項

隨機項預測模型網絡輸入層、隱含層和輸出層神經元個數分別為5個、2個和1個。

圖2 汛期隨機項各數據集觀測值和預測值回歸分析圖

圖3 汛期隨機項觀測值和預測值對照圖

圖2 訓練結果顯示，訓練數據集、測試數據集和檢驗數據集中預測值與觀測值的相關系數R分別為0.9297、0.9553、0.9300，觀測值與預測值的相關系數為0.9332，實測值和預測值吻合精度很高，同時模型對應數據集的均方誤差（MSE）分別為 2.6×10-2、1.3×10-2、1.9×10-2，說明此模型也可以較好的模擬隨機項數據集，從預測值和實測值對比圖（圖3）可以看出，模型對隨機項的模擬精度很高，所選取的氣象因子能很好的反映小覺站汛期徑流量變化的隨機性。

3.2 趨勢項

趨勢項預測模型網絡結構輸入層、隱含層和輸出層神經元個數分別為5個、2個和1個。

圖4 汛期趨勢項各數據集觀測值和預測值回歸分析圖

圖5 汛期趨勢項觀測值和預測值對照圖

圖4訓練結果顯示，訓練數據集、測試數據集和檢驗數據集中預測值與觀測值的相關系數R分別為0.9970、0.9965、0.9986，觀測值與預測值的相關系數為0.9972，實測值和預測值吻合精度很高，同時模型對應數據集的均方誤差（MSE）分別為 1.4×10-3、1.2×10-3、1.0×10-3，說明此模型也可以較好的模擬趨勢項數據集，從預測值和實測值對比圖（圖5）可以看出，模型對趨勢項的模擬精度很高，可見選取的氣象因子能很好的反映小覺站汛期徑流量變化的趨勢。

3.3 周期項

周期項預測模型網絡結構輸入層、隱含層和輸出層神經元個數分別為5個、4個和1個。

圖6 汛期周期項各數據集觀測值和預測值回歸分析圖

圖7 汛期周期項觀測值和預測值對照圖

圖6 訓練結果顯示，訓練數據集、測試數據集和檢驗數據集中預測值與觀測值的相關系數R分別為0.8614、0.8631、0.9289，觀測值與預測值的相關系數為0.8737，實測值與預測值吻合精度較高，同時模型對應數據集的均方誤差（MSE）分別為 3.4×10-2、2.2×10-2、2.5×10-2，說明此模型也可以較好的模擬周期項數據集，從預測值和實測值對比圖（圖7）可以看出選取的氣象因子能較好好的反映小覺站汛期徑流量變化的周期性。

4 小結

圖8 汛期觀測值和預測值對照圖

本文采用主成分分析法提取出汛期氣象因子的主成分，然后采用A Trous小波提取出序列的隨機項、趨勢項和周期項，分別對各項建立神經網絡預測模型，得到較好的預測結果，模型能很好的反映隨機項、趨勢項的變化，對周期項變化的響應較隨機項和趨勢項的差，但也具有較高的精度，說明所建立的模型能很好的反映徑流的隨機性和趨勢性。采用A Trous小波重構法得到徑流預測值，得到汛期重構的預測值和小覺水文站實測值的對比圖（圖8），從預測值與實測值對比來看，模型能較好的模擬小覺站月徑流數據序列中的極值。