董麗麗,徐淑琴,劉 楊,王云鶴

(東北農(nóng)業(yè)大學(xué) 水利與建筑學(xué)院,哈爾濱 150030)

大氣降水是一種復(fù)雜多變的自然現(xiàn)象,目前降水預(yù)報(bào)方法有成因分析法、水文統(tǒng)計(jì)法、模糊分析法、人工神經(jīng)網(wǎng)絡(luò)法等[1-2]。而降雨量序列受噪聲的影響,如果采用含有噪聲的水文序列進(jìn)行分析計(jì)算,將影響數(shù)據(jù)分析成果的精度。本文將采用小波消噪隨機(jī)耦合方法建立降雨量預(yù)報(bào)模型,可以用來(lái)進(jìn)行降水預(yù)報(bào),為水文工作者提供參考,以及指導(dǎo)農(nóng)業(yè)生產(chǎn)。

1 小波消噪隨機(jī)耦合模型的基本原理

實(shí)測(cè)水文序列由于各種因素的干擾,不可避免地含有噪聲。噪聲淹沒(méi)了水文序列的真實(shí)變化規(guī)律,若采用含有噪聲的水文序列進(jìn)行分析,將影響數(shù)據(jù)的可靠性和數(shù)據(jù)分析成果的精度。因此在數(shù)據(jù)建模之前,應(yīng)首先對(duì)序列進(jìn)行消噪處理。目前,小波消噪的主要方法是采用20世紀(jì)90年代Donoho等人提出的閾值消噪方法[3-4]。首先利用小波分解和重構(gòu)消除降雨量序列中的噪聲,然后對(duì)重構(gòu)的降雨量序列應(yīng)用小波隨機(jī)耦合方法進(jìn)行預(yù)測(cè)。

小波隨機(jī)耦合方法的基本思想是采用快速小波變換A Trous算法進(jìn)行小波分解,得到某尺度下的小波變換序列,然后對(duì)各小波變換序列進(jìn)行成分識(shí)別、互相關(guān)分析后,建立各小波變換序列相應(yīng)的數(shù)學(xué)模型,最后采用小波變換重構(gòu)算法得到所研究水文時(shí)間序列的小波隨機(jī)耦合模型[5-6]。

A Trous算法的分解過(guò)程如下[5]:

式中Cj(t)、Wj(t)(j=1,2,…,P)分別為在尺度j下的背景信號(hào)和細(xì)節(jié)信號(hào);P為尺度數(shù); h(k)為離散低通濾波器,h(k)=(1/16,1/4, 3/8,1/4,1/16)。稱 {W1(t),W2(t),…, WP(t),CP(t)}為在尺度P下的小波變換序列。A Trous算法的重構(gòu)過(guò)程如下[3]:

2 小波消噪隨機(jī)耦合模型實(shí)例應(yīng)用

查哈陽(yáng)農(nóng)場(chǎng)地處諾敏河下游、嫩江左岸。為了研究查哈陽(yáng)農(nóng)場(chǎng)的降雨量變化規(guī)律,對(duì)農(nóng)場(chǎng)合理開發(fā)利用水資源,制定合理的灌溉制度具有重要的指導(dǎo)意義。因此,收集了查哈陽(yáng)農(nóng)場(chǎng)1956～2008年每年作物生育期 (5～9月)降雨量觀測(cè)資料,對(duì)該地區(qū)降雨量變化進(jìn)行模擬,并預(yù)測(cè)未來(lái)發(fā)展趨勢(shì)。查哈陽(yáng)農(nóng)場(chǎng)1956～2008年 (5～9月)降雨量變化情況見圖1。

圖1 查哈陽(yáng)農(nóng)場(chǎng)1956～2008月降雨量變化圖Fig.1 Variation curve of month precipitation from 1956 to 2008 in Chahayang Farm

2.1 對(duì)實(shí)測(cè)降雨量序列進(jìn)行消噪平穩(wěn)化處理

由圖1可見,實(shí)測(cè)降雨量序列為非平穩(wěn)時(shí)間序列,需要進(jìn)行一次小波消噪和差分標(biāo)準(zhǔn)化處理,利用式[6](4)、式(5)將其轉(zhuǎn)化為平穩(wěn)序列,見圖2。

圖2 查哈陽(yáng)農(nóng)場(chǎng)月降雨消噪差分后標(biāo)準(zhǔn)化序列變化曲線Fig.2 Curve of month precipitation series after the treatment of denoising,difference and standardization in Chahayang Farm

2.2 對(duì)實(shí)測(cè)降雨量消噪序列進(jìn)行小波分解與重構(gòu)

采用上述的A Trous算法,取P=2,將查哈陽(yáng)農(nóng)場(chǎng)1956～2008年作物生育期月降雨量差分序列分解為W1(t)、W2(t)、C2(t)的變化曲線。將各小波分解序列進(jìn)行疊加,得到重構(gòu)序列,見圖3中(a)、(b)、(c)、(d)。

2.3 小波變換序列成分識(shí)別和互相關(guān)分析

圖3 查哈陽(yáng)農(nóng)場(chǎng)月降雨序列小波消噪分解與重構(gòu)Fig.3 Wavelet decomposition and reconstruction of the month precipitation series after the treatment of denoising in Chahayang Farm

根據(jù)式(6)、式(7)、式(8)計(jì)算各小波變換序列的自相關(guān)系數(shù)、方差譜密度和互相關(guān)系數(shù),并繪制自相關(guān)圖、方差譜密度圖和互相關(guān)圖,見圖4、圖5和圖6。根據(jù)圖3、圖4、圖5分析各小波變換序列的變化特性,可以近似認(rèn)為W1(t)代表原序列的隨機(jī)項(xiàng),W2(t)代表原序列的周期項(xiàng), C2(t)代表原序列的趨勢(shì)項(xiàng)。

由圖6可以看出,各小波變換序列互相關(guān)系數(shù)基本上落在95%容許限范圍以內(nèi),而且趨近于0。所以,可以認(rèn)為小波變換序列W1(t)、W2(t)和C2(t)互相關(guān)性較小,兩兩獨(dú)立。因此分析和處理月降雨量差分序列就轉(zhuǎn)化為對(duì)W1(t)、W2(t)和C2(t)的處理。

圖6 各小波變化序列互相關(guān)圖Fig.6 Cross-correlation map of each wavelet transform series

2.4 各小波變換序列模型的構(gòu)建

2.4.1 細(xì)節(jié)序列W1(t)隨機(jī)模型的建立

通過(guò)MAT LAB編程計(jì)算,細(xì)節(jié)序列W1(t)的均值為0.000 428 14≈0,方差為1.531 9,偏態(tài)系數(shù)為-0.268 6,所以認(rèn)為序列W1(t)需要進(jìn)行正態(tài)性轉(zhuǎn)化[7],經(jīng)正態(tài)性轉(zhuǎn)換后,W1(t)序列將服從ε～(0,0.007 9)的正態(tài)分布。

對(duì)序列W1(t)分別進(jìn)行自相關(guān)分析和偏相關(guān)分析,自相關(guān)圖具有拖尾性,而偏相關(guān)圖具有截尾性,所以初步判定模型形式為AR(p)模型。參考文獻(xiàn) [8-9],初步判定模型階數(shù)為6,屬于AR (p)模型。采用AIC準(zhǔn)則對(duì)AR(p)模型的階數(shù)進(jìn)行進(jìn)一步識(shí)別,AIC(6)=-364.363 7,AIC (5)=-19.535 3,AIC(7)=-100.785 7,因此確定模型的階數(shù)為6階,選定k=1∶6。對(duì)AR (6)模型參數(shù)進(jìn)行MAT LAB編程計(jì)算,建立如下自回歸模型:

2.4.2 細(xì)節(jié)序列W2(t)模型的建立

通過(guò)MAT LAB編程計(jì)算,細(xì)節(jié)序列W2(t)的均值為-0.000 063 22≈0,方差為0.230 7,偏態(tài)系數(shù)為-0.005≈0,所以可以認(rèn)為細(xì)節(jié)序列W2(t)是正態(tài)分布的。經(jīng)過(guò)對(duì)自相關(guān)圖和互相關(guān)圖的判斷,確定細(xì)節(jié)序列W2(t)模型屬于AR模型,模型的階數(shù)是3,對(duì)于序列W2(t),選定k=1, 2,3(偏相關(guān)系數(shù)超出95%容許限范圍)。建立自回歸模型為:

2.4.3 背景序列C2(t)模型的建立

背景序列C2(t)的均值為-0.000 491 36≈0,方差為0.057,偏態(tài)系數(shù)為-0.017 4≈0,所以可以認(rèn)為細(xì)節(jié)序列C2(t)是正態(tài)分布的。經(jīng)過(guò)對(duì)自相關(guān)圖和互相關(guān)圖的判斷,確定細(xì)節(jié)序列 C2(t)模型屬于AR模型,模型的階數(shù)是2,建立的背景序列C2(t)模型為:

2.5 小波消噪隨機(jī)耦合模型組合及擬合

將上述各小波變換序列模型疊加,就可以得到查哈陽(yáng)農(nóng)場(chǎng)月降雨量小波消噪隨機(jī)耦合模型,即:

采用建立的小波消噪隨機(jī)耦合模型對(duì)查哈陽(yáng)農(nóng)場(chǎng)1956～2008年的 (5～9月)降雨量進(jìn)行擬合,見圖7。

圖7 查哈陽(yáng)農(nóng)場(chǎng)月降雨量擬合曲線(1956～2008)Fig.7 Fitting curve of month precipitation from 1956 to 2008 in Chahayang Farm

2.6 小波消噪隨機(jī)耦合模型精度檢驗(yàn)及預(yù)測(cè)

對(duì)所建小波消噪隨機(jī)耦合模型的擬合數(shù)據(jù)進(jìn)行擬合效果檢驗(yàn)[10],擬合效果評(píng)價(jià)指標(biāo)后驗(yàn)差比值C=0.04,小誤差頻率p=1。經(jīng)過(guò)計(jì)算,平均相對(duì)誤差e=2.35%,相對(duì)最大誤差為2.77%,擬合準(zhǔn)確率1。對(duì)2005～2008年作物生育期的降雨量進(jìn)行預(yù)測(cè)檢驗(yàn),見表1,預(yù)測(cè)檢驗(yàn)合格率為100% (相對(duì)誤差e＜20%為合格),達(dá)到1級(jí)標(biāo)準(zhǔn)。運(yùn)用該模型對(duì)2009～2015年每年 (5～9月)降雨量進(jìn)行預(yù)測(cè),預(yù)測(cè)曲線見圖8。

表1 查哈陽(yáng)農(nóng)場(chǎng)作物生育期降雨量預(yù)測(cè)檢驗(yàn)Table 1 Pprediction test of the precipitation data of crop growth period in Chahayang Farm

圖8 查哈陽(yáng)農(nóng)場(chǎng)月降雨量預(yù)測(cè)曲線(2009～2015)Fig.8 month prediction forecasting curve from 2009 to2015 in Chahayang Farm

3 結(jié) 論

選擇53 a的查哈陽(yáng)農(nóng)場(chǎng)月降雨量觀測(cè)資料進(jìn)行分析,通過(guò)對(duì)各小波變換序列進(jìn)行互相關(guān)分析,用小波重構(gòu)算法得到降雨預(yù)測(cè)模型。由圖7可見資料的序列越長(zhǎng),建立的模型擬合度越好,因此預(yù)報(bào)精度也越高,該模型能夠反映該地區(qū)降雨量變化規(guī)律,是一種實(shí)用的預(yù)報(bào)模型,為防洪抗旱,合理制定灌溉制度,指導(dǎo)生產(chǎn),都是有益處的。同時(shí)也為當(dāng)?shù)毓芾聿块T管理和規(guī)劃水資源提供可靠依據(jù)。

[1]楊 旭,欒繼虹,馮國(guó)章.中長(zhǎng)期水文預(yù)報(bào)研究評(píng)述與展望 [J].西北農(nóng)業(yè)大學(xué)學(xué)報(bào),2000,28(6): 203-207.

[2]戴長(zhǎng)雷,遲寶明,李治軍.基于回歸分析與時(shí)序分析降水預(yù)報(bào)迭合模型的構(gòu)建與實(shí)現(xiàn)[J].河南師范大學(xué)學(xué)報(bào):自然科學(xué)版,2006,34(1):15-18.

[3]王文圣,丁 晶,李躍清.水文小波分析 [M].北京:化學(xué)工業(yè)出版社,2005.

[4]D L Donoho.Denoising by Soft-thresholding[J]. IEEE Transactions on Information Theory,1995,41 (3):613-627.

[5]徐淑琴,付 強(qiáng),孟凡香,等.基于小波隨機(jī)耦合模型的地下水動(dòng)態(tài)變化規(guī)律研究 [J].中國(guó)農(nóng)村水利水電,2008,(4):33-36.

[6]劉 東,付 強(qiáng).小波隨機(jī)耦合模型在三江平原井灌區(qū)地下水位預(yù)測(cè)中的應(yīng)用[J].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí), 2008,38(16):54-62.

[7]付 強(qiáng).數(shù)據(jù)處理方法及其應(yīng)用 [M].北京:科學(xué)出版社,2006:49-61.

[8]丁 晶,劉權(quán)授.隨機(jī)水文學(xué)[M].北京:中國(guó)水利水電出版社,1997:54-59.

[9]付 強(qiáng).農(nóng)業(yè)水土資源系統(tǒng)分析與綜合評(píng)價(jià) [M].北京:中國(guó)水利水電出版社,2005:228-235.

[10]陳南祥.地下水動(dòng)態(tài)預(yù)報(bào)模型的精度評(píng)價(jià) [J].工程勘察,1999,(3):35-38.