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

（1.貴州省水利水電勘測(cè)設(shè)計(jì)研究院，貴州 貴陽(yáng) 550002；2.中國(guó)農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，北京 100083）

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

1 研究區(qū)概況

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

圖1 滹沱河流域范圍

2 方法簡(jiǎn)介

2.1 主成分分析法

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

主成分分析的數(shù)學(xué)模型為

可用矩陣表示為

經(jīng)過(guò)標(biāo)準(zhǔn)化處理的矩陣就是X的相關(guān)系數(shù)矩陣R，如果主成分分析的一切計(jì)算不是從協(xié)方差矩陣出發(fā)，而是直接從樣本相關(guān)系數(shù)矩陣出發(fā)，就等于先對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，然后再利用協(xié)方差矩陣進(jìn)行主成分分析[2]。主成分分析的基本步驟是：首先結(jié)合數(shù)據(jù)，判斷是否需要進(jìn)行主成分分析；其次結(jié)合主成分的特征值和累計(jì)特征值貢獻(xiàn)率來(lái)確定抽取的主成分或因子的數(shù)據(jù)；最后，進(jìn)行主成分分析，將抽取出來(lái)的主成分分析存為新變量，以便繼續(xù)分析[3]。

2.2 小波變換—A Trous

對(duì)水文時(shí)間序列f（t），t=1，2，…，N進(jìn)行小波分解，令C0（t）=f（t），A Trous小波分解過(guò)程如下：

A Trous小波重構(gòu)過(guò)程如下：

根據(jù)樣本序列長(zhǎng)度，低通濾波器選用db4小波，選擇最大分解尺度P=3，采用MATLAB小波分解工具箱分解得到樣本的隨機(jī)項(xiàng)、趨勢(shì)項(xiàng)和周期項(xiàng)，采用BP神經(jīng)網(wǎng)絡(luò)對(duì)各項(xiàng)分別進(jìn)行預(yù)測(cè)，再重構(gòu)得到預(yù)測(cè)模型的輸出值[4]。

3 汛期月徑流預(yù)測(cè)

汛期采用1969年～2009年的7月、8月、9月、10月的16個(gè)氣象因子：x1五臺(tái)山降水量、x2五臺(tái)山風(fēng)速，x3五臺(tái)山氣壓、x4五臺(tái)山氣溫、x5五臺(tái)山相對(duì)濕度、x6五臺(tái)山日照時(shí)數(shù)、x7五臺(tái)山水汽壓、x8原平降水量、x9原平風(fēng)速、x10原平氣壓、x11原平氣溫、x12原平相對(duì)濕度、x13原平日照時(shí)數(shù)、x14原平水汽壓、x15小覺(jué)站降水量、x16小覺(jué)站蒸發(fā)量，以及小覺(jué)站徑流量值作為分析基礎(chǔ)。

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

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

表2 主成分表達(dá)式的系數(shù)矩陣（汛期）

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

3.1 隨機(jī)項(xiàng)

隨機(jī)項(xiàng)預(yù)測(cè)模型網(wǎng)絡(luò)輸入層、隱含層和輸出層神經(jīng)元個(gè)數(shù)分別為5個(gè)、2個(gè)和1個(gè)。

圖2 汛期隨機(jī)項(xiàng)各數(shù)據(jù)集觀測(cè)值和預(yù)測(cè)值回歸分析圖

圖3 汛期隨機(jī)項(xiàng)觀測(cè)值和預(yù)測(cè)值對(duì)照?qǐng)D

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

3.2 趨勢(shì)項(xiàng)

趨勢(shì)項(xiàng)預(yù)測(cè)模型網(wǎng)絡(luò)結(jié)構(gòu)輸入層、隱含層和輸出層神經(jīng)元個(gè)數(shù)分別為5個(gè)、2個(gè)和1個(gè)。

圖4 汛期趨勢(shì)項(xiàng)各數(shù)據(jù)集觀測(cè)值和預(yù)測(cè)值回歸分析圖

圖5 汛期趨勢(shì)項(xiàng)觀測(cè)值和預(yù)測(cè)值對(duì)照?qǐng)D

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

3.3 周期項(xiàng)

周期項(xiàng)預(yù)測(cè)模型網(wǎng)絡(luò)結(jié)構(gòu)輸入層、隱含層和輸出層神經(jīng)元個(gè)數(shù)分別為5個(gè)、4個(gè)和1個(gè)。

圖6 汛期周期項(xiàng)各數(shù)據(jù)集觀測(cè)值和預(yù)測(cè)值回歸分析圖

圖7 汛期周期項(xiàng)觀測(cè)值和預(yù)測(cè)值對(duì)照?qǐng)D

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

4 小結(jié)

圖8 汛期觀測(cè)值和預(yù)測(cè)值對(duì)照?qǐng)D

本文采用主成分分析法提取出汛期氣象因子的主成分，然后采用A Trous小波提取出序列的隨機(jī)項(xiàng)、趨勢(shì)項(xiàng)和周期項(xiàng)，分別對(duì)各項(xiàng)建立神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型，得到較好的預(yù)測(cè)結(jié)果，模型能很好的反映隨機(jī)項(xiàng)、趨勢(shì)項(xiàng)的變化，對(duì)周期項(xiàng)變化的響應(yīng)較隨機(jī)項(xiàng)和趨勢(shì)項(xiàng)的差，但也具有較高的精度，說(shuō)明所建立的模型能很好的反映徑流的隨機(jī)性和趨勢(shì)性。采用A Trous小波重構(gòu)法得到徑流預(yù)測(cè)值，得到汛期重構(gòu)的預(yù)測(cè)值和小覺(jué)水文站實(shí)測(cè)值的對(duì)比圖（圖8），從預(yù)測(cè)值與實(shí)測(cè)值對(duì)比來(lái)看，模型能較好的模擬小覺(jué)站月徑流數(shù)據(jù)序列中的極值。