混凝土壩變形Wavelet-EGM-PE-ARIMA組合預(yù)測模型

(1.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，南京 210098； 2.河海大學(xué) 水資源高效利用與工程安全國家工程研究中心，南京 210098； 3.河海大學(xué) 水利水電學(xué)院，南京 210098; 4. 河北農(nóng)業(yè)大學(xué) 理工學(xué)院，河北保定 071066)

1 研究背景

自20世紀(jì)50年代以來，我國共修建約9.8萬座堤壩，取得了舉世矚目的成就，其中大部分是混凝土壩，它們在防洪、發(fā)電和灌溉等方面取得了巨大的經(jīng)濟(jì)效益[1]。變形是大壩壩體、壩基結(jié)構(gòu)性態(tài)變化的綜合反映，是衡量大壩結(jié)構(gòu)是否正常的重要標(biāo)志，因此，科學(xué)地開展混凝土壩變形性態(tài)安全監(jiān)控理論的研究，并對未來某時(shí)刻大壩變形進(jìn)行有效地預(yù)測，是壩工界的重要研究課題。

到目前為止，關(guān)于混凝土壩變形預(yù)測方法的研究已取得了一些成果。例如：陳久宇教授[2]使用統(tǒng)計(jì)回歸分析方法成功分離變形效應(yīng)量監(jiān)測數(shù)據(jù)中的分量；吳中如等[3-6]全面解讀了大壩變形預(yù)測統(tǒng)計(jì)模型、混合模型及確定性模型；顧沖時(shí)等[7-8]建立了大壩變形預(yù)測時(shí)空分布模型。此外，隨智能算法的興起，同時(shí)考慮到傳統(tǒng)模型相對智能算法的擬合精度較低，灰色理論、小波分析、神經(jīng)網(wǎng)絡(luò)等智能算法不斷被引進(jìn)到混凝土壩變形預(yù)測體系中來，例如：羅德河等[9]通過小波分解主值函數(shù)項(xiàng)、周期函數(shù)項(xiàng)和改進(jìn)后的平穩(wěn)時(shí)間序列，建立大壩變形的小波分析與ARMA預(yù)測模型；俞艷玲等[10]通過蛙跳算法優(yōu)化背景值和平滑系數(shù)、尋找最優(yōu)定解條件以及殘差優(yōu)化等方法，提出了改進(jìn)的非等間距GM(1,1)大壩位移預(yù)測模型;萬臣等[11]借助新維BP神經(jīng)網(wǎng)絡(luò)對時(shí)間序列的滾動(dòng)預(yù)測和馬爾科夫鏈模型對隨機(jī)擾動(dòng)誤差進(jìn)行修正的特點(diǎn)，構(gòu)建了一種新維BP神經(jīng)網(wǎng)絡(luò)-馬爾科夫鏈大壩沉降預(yù)測模型。這些方法已經(jīng)成功應(yīng)用于實(shí)際工程中，并取得了較大的社會、經(jīng)濟(jì)效應(yīng)，但也存在一些局限性。

大壩變形時(shí)間序列可視為一系列不同頻率信號的疊加，其中，水壓變形和溫度變形主要由庫水位的變化和溫度的變化引起，表現(xiàn)為高頻信號，具有明顯的周期性特征，可視為總變形中的周期性分量；而時(shí)效變形主要反映壩體混凝土和基巖的徐變、塑性變形以及基巖地質(zhì)結(jié)構(gòu)的壓縮變形，同時(shí)還包括壩體裂縫引起的不可逆位移以及自身體積變形[1]，表現(xiàn)為低頻信號，體現(xiàn)為總變形中的趨勢性分量[12]。因此，本文充分利用現(xiàn)代計(jì)算機(jī)技術(shù)，基于復(fù)合建模思想，分別針對不同頻率的信號建立模型，以實(shí)現(xiàn)對混凝土壩變形時(shí)間序列的有效預(yù)測。

2 組合模型建立理論

2.1 小波理論

小波分析具有良好的視頻局域化性質(zhì)，因此可利用小波理論多分辨分析的特點(diǎn)[13]對多頻信號進(jìn)行分解。大壩變形時(shí)間序列可以視為由不同頻率的信號序列組成，其中水壓和溫度變化具有一定的規(guī)律性，呈現(xiàn)一定的周期性，兩者引起的變形表現(xiàn)為高頻信號；而長期的時(shí)效變形則表現(xiàn)為低頻信號。因此可將高頻信號當(dāng)成由水壓和溫度變化引起的周期性分量F，將低頻信號當(dāng)成時(shí)效引起的趨勢性分量M?；诖耍疚母鶕?jù)小波多尺度分解出高、低頻信號，實(shí)現(xiàn)對周期性分量F和趨勢性分量M的提取。

實(shí)際上，大壩變形時(shí)間序列f(t)為離散數(shù)據(jù)序列，故采用離散小波變換對大壩變形時(shí)間序列進(jìn)行處理分析。小波變換可表示為[14]：

(1)

(2)

(3)

多尺度離散小波分解和重構(gòu)一般通過Mallat算法得以實(shí)現(xiàn)[15-16]。假設(shè)將信號分解為n層，其分解步驟如圖1所示[17]：首先將原始信號f(t)分為高頻信號d1和低頻信號a1，再將低頻信號a1分為高頻信號d2和低頻信號a2，以此類推直至分解得到高頻信號dn和低頻信號an。

圖1 小波分解步驟Fig.1 Wavelet decomposition steps

圖1中的關(guān)系可表示為:

f(t)=an+d1+d2+…+dn。

(4)

2.2 EGM(1,1)模型

本文選取的大壩變形序列時(shí)間較短，變形實(shí)測值數(shù)據(jù)不多，而灰色系統(tǒng)理論恰好主要對樣本量少且信息貧乏的不確定系統(tǒng)進(jìn)行研究挖掘，從而實(shí)現(xiàn)由部分到系統(tǒng)、由已知到未知的目標(biāo)，最終對系統(tǒng)的特征進(jìn)行合理闡述，且進(jìn)行預(yù)測預(yù)報(bào)。GM(1,1)模型是灰色系統(tǒng)的基本模型，具有多種形式，其中EGM模型最適合非指數(shù)增長[18-19]，故本文采用EGM(1,1)模型。設(shè)原始數(shù)據(jù)序列為X(0)=[x(0)(1),x(0)(2),…,x(0)(n)]，EGM(1,1)建模基本步驟如下。

Step 1：引入二階弱化算子D2，令

X(0)D=[x(0)(1)d1,x(0)(2)d2,…,x(0)(n)dn]，

(5)

令

X(0)D2=[x(0)(1)d12,x(0)(2)d22,…,x(0)(n)dn2]。

(6)

Step 2：對序列X(0)D2進(jìn)行累加，由此生成X(0)的1-GAO序列為

X(1)=[x(1)(1),x(1)(2),…,x(1)(n)] 。

(7)

(8)

其中

k=1，2，…，n-1。

k=1，2，…，n-1 。

(9)

式(9)即為預(yù)測方程，由于該式是對原始數(shù)據(jù)進(jìn)行累加計(jì)算，故需進(jìn)行一次累減處理，得到還原序列見式(10)，代入不同的k值可得到一組擬合數(shù)據(jù)。

k=1，2，…，n-1 。

(10)

2.3 周期外延模型

周期外延PE模型把一個(gè)復(fù)雜的時(shí)間序列看成由有限個(gè)具有不同周期的周期波相互疊加而成。而大壩變形時(shí)間序列存在由水壓和溫度造成的周期性項(xiàng)，故可以提前對大壩變形時(shí)間序列進(jìn)行小波分解提取周期性分量，再通過周期外延模型識別其周期波，再將識別的周期波分別外延，最后進(jìn)行疊加處理進(jìn)行預(yù)測。設(shè)周期性序列為F={f(k)}=[f(1),f(2),…,f(n)]，PE模型的建模步驟如下。

Step 1：計(jì)算均值生成函數(shù)。該函數(shù)計(jì)算公式[20]為

i=1，2，…，m， 1≤m≤M。

(11)

式中:nm表示為≤n/m的最大整數(shù);M表示為≤n/2的最大整數(shù)。由此均值生成函數(shù)可用矩陣式(12)表示。

(12)

k=1，2，…，n。

(13)

Step 2：提取優(yōu)勢周期。采用方差分析的方法，通過式(14)檢驗(yàn)序列{f(k)}是否存在長度為m的隱含的優(yōu)勢周期。

(14)

(15)

(16)

式中F(m)表示自由度為(m-1,n-m)的F分布。

本文中，設(shè)置信水平α=0.5。如果F(m)>Fα(m-1,n-m)，則可認(rèn)為x(k)隱含的優(yōu)勢周期長度為m。

Step 3：序列{f(k)}減去Step 1中周期為m的函數(shù)得到新序列{f′(k)}，然后對新序列{f′(k)}循環(huán)Step 1和Step 2，從而提取存在的其他優(yōu)勢周期。

Step 4：優(yōu)勢周期的疊加。疊加同一時(shí)刻的不同周期可得

(18)

式(18)即為PE模型。

2.4 ARIMA模型

ARIMA模型將預(yù)測對象隨時(shí)間推移而形成的數(shù)據(jù)序列視為一個(gè)隨機(jī)序列，用一定的數(shù)學(xué)模型來近似描述這個(gè)序列，這個(gè)模型一旦被識別就可以通過時(shí)間序列的過去值及現(xiàn)在值來預(yù)測未來值。

在本文中，由EGM模型和PE模型殘差項(xiàng)組成的隨機(jī)性分量Z可視為非平穩(wěn)隨機(jī)時(shí)間序列，在此可采用ARIMA模型首先對隨機(jī)性分量Z進(jìn)行差分處理，將其轉(zhuǎn)化為平穩(wěn)時(shí)間序列，再進(jìn)行回歸建模處理。

該模型的數(shù)學(xué)表達(dá)式為

(19)

式中：φm(m=1,2,…，p)表示自回歸系數(shù)；θj(j=1,2,…，q)表示移動(dòng)平均系數(shù)；p表示自回歸部分的階數(shù)；q為移動(dòng)平均部分的階數(shù)；at為白噪音部分。

ARIMA模型的實(shí)現(xiàn)步驟[21]如下。

Step 1：計(jì)算差分階數(shù)d。本文中的隨機(jī)性分量Z需差分處理轉(zhuǎn)化為平穩(wěn)時(shí)間序列，差分階數(shù)d經(jīng)過相關(guān)圖檢驗(yàn)法確定為1。

Step 2：計(jì)算自回歸階數(shù)p和移動(dòng)平均階數(shù)q。通常采用最小信息準(zhǔn)則(AIC)和貝葉斯信息準(zhǔn)則(BIC)來確定p和q，本文考慮到BIC準(zhǔn)則收斂性相對較好，最終通過BIC準(zhǔn)則確定q=2和q=0。

Step 3：建立ARIMA模型。由Step 1 和Step 2確定的各系數(shù)，建立ARIMA(2,1,0)模型對隨機(jī)性分量Z進(jìn)行擬合預(yù)測。

3 組合模型建立步驟

綜上所述，實(shí)現(xiàn)混凝土壩變形Wavelet-EGM-PE-ARIMA組合預(yù)測模型的流程見圖2，建模具體步驟如下。

圖2 組合預(yù)測模型流程Fig.2 Flowchart of the combinatorial prediction model

Step 1：提取趨勢性分量M和周期性分量F。利用小波多分辨分析功能，對大壩變形時(shí)間序列f(t)進(jìn)行多尺度分解出高頻信號d1—d7和低頻信號a7，趨勢性分量M即為低頻信號a7，周期性分量F即為高頻信號的疊加。

Step 2：建立趨勢性分量M和周期性分量F擬合預(yù)測模型。運(yùn)用EGM(1,1)模型對Step 1中提取的趨勢性分量M進(jìn)行建模擬合預(yù)測，運(yùn)用周期外延模型對Step 1中提取的周期性分量F分別進(jìn)行建模擬合預(yù)測，并得到殘差項(xiàng)Δz1和Δz2組成隨機(jī)性分量Z。

Step 3：建立隨機(jī)性分量Z擬合預(yù)測模型,在Step 2基礎(chǔ)上，運(yùn)用ARIMA模型對隨機(jī)性分量Z進(jìn)行建模擬合預(yù)測。

Step 4：建立組合模型。通過小波重構(gòu)Step 2中的3個(gè)模型得到組合預(yù)測模型，實(shí)現(xiàn)對大壩變形的預(yù)報(bào)預(yù)測。

4 工程實(shí)例

某混凝土重力壩最大壩高181 m，主要由河床泄洪壩段、左右岸廠房壩段和左右岸非溢流壩段組成。其中右岸廠房壩段布置多條倒垂線，本文選取右岸廠房17#壩段IP01yc17測點(diǎn)的實(shí)測順河向位移，建立混凝土壩變形Wavelet-EGM-PE-ARIMA組合預(yù)測模型。選取2013年9月1日至2016年6月5日(測值間隔7 d)的觀測數(shù)據(jù)，共計(jì)141個(gè)，其中對前135個(gè)(2013年9月1日至2016年3月27日)觀測數(shù)據(jù)進(jìn)行建模分析，后6個(gè)(2016年4月3日至2016年5月8日)觀測數(shù)據(jù)進(jìn)行預(yù)測對比分析。

本文主要通過MatLab和SPSS兩個(gè)軟件平臺進(jìn)行擬合處理及預(yù)測實(shí)現(xiàn)。由圖2中的流程可知：首先進(jìn)行小波多尺度分解提取高、低信號，通過試算，選擇sym8小波基對大壩變形時(shí)間序列進(jìn)行7層分解(見圖3)，然后分別提取趨勢性分量M和周期性分量F見圖4。

圖3 大壩變形時(shí)間序列的7層小波分解Fig.3 Seven levels of wavelet decomposition>of dam deformation time series

圖4 小波分解提取的趨勢性分量M和周期性分量FFig.4 Trend component M and periodic component Fextracted by wavelet decomposition

對于小波提取的趨勢性分量M，在MatLab平臺上通過編程建立EGM(1,1)模型進(jìn)行建模預(yù)測，并得到殘差Δz1。趨勢性分量M的擬合曲線及殘差Δz1見圖5，圖中趨勢性分量M的實(shí)測過程線和擬合過程線位于左縱坐標(biāo)軸，殘差Δz1位于右縱坐標(biāo)軸。

圖5 均值EGM(1,1)模型擬合曲線Fig.5 Fitted curve of EGM(1,1) model

對于小波提取的周期性分量F，在MatLab平臺通過編程建立PE模型進(jìn)行建模預(yù)測，并得到殘差Δz2。周期性分量F的擬合曲線及殘差Δz2見圖6。圖中趨勢性分量F的實(shí)測過程線和擬合過程線位于左縱坐標(biāo)軸，殘差Δz2位于右縱坐標(biāo)軸。

圖6 周期外延模型擬合曲線Fig.6 Fitted curve of PE model

對于由前2次擬合得到的殘差項(xiàng)組成的隨機(jī)性分量Z，本文通過在SPSS平臺上建立ARIMA模型對其建模預(yù)測。在建模過程中，通過相關(guān)圖檢驗(yàn)法和BIC準(zhǔn)則，最終確定ARIMA(2,1,0)。隨機(jī)性分量Z擬合曲線見圖7。

圖7 ARIMA模型擬合曲線Fig.7 Fitted curve of ARIMA model

最后通過小波重構(gòu)建立組合模型，繪制傳統(tǒng)模型和組合模型的擬合曲線見圖8，并采用復(fù)相關(guān)系數(shù)R和均方差s來比較2種模型的擬合精度,見表1。結(jié)果表明：相對傳統(tǒng)模型，該組合模型擬合精度較高，具有更好的擬合精度。

圖8 2種模型擬合曲線Fig.8 Fitted curves of two models

模型Rs組合模型0.9990.011統(tǒng)計(jì)模型0.9930.066

以2016年4月3日至2016年5月8日的6次測值作為預(yù)測對比依據(jù)，分別求得統(tǒng)計(jì)模型和該組合模型預(yù)測值，繪制2種模型預(yù)測值比較折線,見圖9，同時(shí)采用相對誤差來比較2種模型的預(yù)測精度,見表2。結(jié)果表明本文的組合模型預(yù)測精度相對統(tǒng)計(jì)模型較高。

圖9 2種模型預(yù)測值比較Fig.9 Comparison of prediction result betweentwo models

日期實(shí)測值/mm組合模型統(tǒng)計(jì)模型預(yù)測值/mm相對誤差/%預(yù)測值/mm相對誤差/%2016-04-033.123.171.503.081.212016-04-103.083.090.603.080.012016-04-173.153.242.953.053.232016-04-243.203.344.283.025.542016-05-013.303.251.412.9610.262016-05-083.253.191.862.8911.01

5 結(jié) 語

針對大壩水壓變形、溫度變形和時(shí)效變形的特點(diǎn)，基于小波多辨分析方法、均值EGM模型、周期外延模型以及ARIMA模型，提出一種混凝土壩變形Wavelet-EGM-PE-ARIMA組合預(yù)測模型。結(jié)合工程實(shí)例可知：

(1)該模型針對大壩變形時(shí)間序列的特點(diǎn)進(jìn)一步劃分出趨勢性項(xiàng)、周期性項(xiàng)，有助于加深對大壩變形時(shí)間序列變化規(guī)律的理解。

(2)該模型針對大壩變形監(jiān)測數(shù)據(jù)中的趨勢性項(xiàng)、周期性項(xiàng)以及隨機(jī)性項(xiàng)不同特點(diǎn)，分別進(jìn)行建模分析，能夠有效地提高模型擬合和預(yù)測精度，具有一定的應(yīng)用價(jià)值。