劉彩花,景 浩,成一雄,牛 月
(太原理工大學(xué)水利科學(xué)與工程學(xué)院,山西 太原 030024)
國內(nèi)外統(tǒng)計資料顯示,地下滲漏引起大壩發(fā)生安全事故的概率為20%,失事原因為滲流量增大、滲透坡降加大等[1]。我國1981年的調(diào)查結(jié)果表明,31.7%的大壩安全事故是由滲漏管涌引起的[2]。因此,對大壩進行滲流監(jiān)測,通過建立恰當模型來監(jiān)控大壩運行狀態(tài)、預(yù)測大壩未來運行性狀意義重大。
由于滲流序列體現(xiàn)出復(fù)雜的非線性、非平穩(wěn)特征,可能包含周期性、趨勢性及隨機波動等成分,因此聯(lián)合運用多種理論分析方法建立耦合模型取代單一模型預(yù)測成為研究趨勢。本文利用經(jīng)驗?zāi)B(tài)分解法 (Empirical Mode Decomposition,EMD) 具有秉留原數(shù)據(jù)屬性的特點及時間序列ARI(Auto Regressive Integrated)模型具有預(yù)測精度較高與先進統(tǒng)計性的優(yōu)勢,提出了EMD-ARI滲流耦合模型,將此模型應(yīng)用于汾河水庫左壩岸滲流量的預(yù)報中,并將預(yù)報結(jié)果與ARI單一時序模型的預(yù)報結(jié)果進行比較和分析,以期為大壩安全運行狀態(tài)監(jiān)控及模型構(gòu)建提供參考。
EMD[3]可將原始序列分解為若干階頻率不同的固有模態(tài)函數(shù) (Intrinsic Mode Function,IMF)和一階趨勢項,分解步驟為:
(1)找出原始序列x(t)全部極大值點,進行3次樣條函數(shù)法插值,得各時刻上、下包絡(luò)線值xmax(t)、xmin(t) 及包絡(luò)線均值 m(t)
(2) 求類距平值函數(shù) h(t)
判斷h(t)是否為固有模態(tài)函數(shù)的依據(jù)為:①h(t)函數(shù)的極值點個數(shù)與跨零點個數(shù)相等或僅相差1;②m(t)恒等于零。滿足則h(t)為固有模態(tài)函數(shù),否則將 h(t)作為原始序列重復(fù)步驟 (1)、(2),直至同時滿足①、②為止。得到的一階固有模態(tài)函數(shù)I1(t)及剩余值序列為
(3)以r1(t)為新序列,重復(fù)以上步驟,直至分解出所有階IMF分量及一階趨勢項r(t)。
ARI[4,5](p,d)模型是d次差分運算與AR(p)模型的組合,表達式為:
式中,Xt為原序列;φ1,φ2,…,φp為自回歸系數(shù);p 為自回歸模型階數(shù);B為后移算子;B的次數(shù)表示后移期數(shù);d為差分運算次數(shù);Φ(B)為自回歸系數(shù)多項式;at為零均值白噪聲序列。
汾河水庫建于1958年,為大 (二)型水利樞紐,具有發(fā)電、防洪、灌溉等綜合效益。水庫于1990年建立了大壩滲流自動觀測系統(tǒng),監(jiān)測內(nèi)容包括左右壩基滲流、左壩岸滲流、古河床壩體繞滲、浸潤線等,其中,左壩岸滲流采用三角堰觀測,觀測頻率為每周一次,觀測設(shè)備完好。本文采用2010年7月2日~2012年11月17日汾河水庫左壩岸共125組滲流監(jiān)測數(shù)據(jù)進行模型構(gòu)建;采用2012年11月24日~2012年12月29日汾河水庫左壩岸共6組滲流監(jiān)測數(shù)據(jù)檢驗?zāi)P蛿M合結(jié)果與預(yù)報精度,部分滲流監(jiān)測數(shù)據(jù)如表1所示。
表1 汾河水庫左壩岸部分滲流量監(jiān)測數(shù)據(jù)L·s-1
基于MATLAB平臺對汾河水庫左壩岸滲流序列進行EMD分解,分解結(jié)果如圖1所示。分解得4階IMF分量與趨勢項,經(jīng)Hilbert變換及倒數(shù)運算后求得 4階 IMF分量的平均周期分別為 23.8、43.6、148、222.7天,表明滲流序列可能存在23.8、43.6、148、222.7天這4個波動周期;1階趨勢項總體呈上升趨勢,表明汾河水庫左壩岸滲流量有隨時間推移增加的趨勢。
圖1 汾河水庫左壩岸滲流序列EMD分解結(jié)果
應(yīng)用Fourier級數(shù)對4階IMF分量進行周期擬合[6],不考慮均值項的擬合結(jié)果為
應(yīng)用線性方程對趨勢項進行擬合,擬合結(jié)果為
式中,t為滲流量序列時間序號(1,2,…,125)。
水庫滲流序列經(jīng)提取可能周期與趨勢項后得剩余值序列,對其進行時序分析建立ARI模型的步驟[7]為:
(1)剩余值序列平穩(wěn)性檢驗。ADF單位根檢驗結(jié)果表明剩余值序列顯著非平穩(wěn),一階差分后平穩(wěn)。
(2)差分序列的純隨機性檢驗。采用LB檢驗統(tǒng)計量對差分后的剩余值序列進行純隨機性檢驗,檢驗結(jié)果為LB統(tǒng)計量的P(Pr>ChiSq)值均小于0.05,說明該序列為非白噪聲序列,滿足時序建模條件。
(3)差分序列擬合AR模型。模型階數(shù)由BIC準則確定,BIC值最小時相應(yīng)模型階數(shù)最優(yōu),其表達式如下
式中,BIC為信息準則函數(shù);N為序列長度;p為模型階數(shù);C為常數(shù);為模型殘差方差的極大似然估計值。
根據(jù)BIC準則確定最優(yōu)模型為AR(2),模型參數(shù)經(jīng)檢驗極顯著,還原差分后得ARI(2,1)模型表達式為
(2) 模型殘差檢驗。 ARI(2,1)模型殘差的檢驗結(jié)果為LB檢驗統(tǒng)計量的P(Pr>ChiSq)值均大于0.05,表明ARI(2,1)模型適合滲流序列。
汾河水庫EMD-ARI模型的結(jié)構(gòu)為
式中,yt為汾河水庫左壩岸滲流量序列;aj、bj為相應(yīng)各階IMF分量的傅氏系數(shù);ωj為相應(yīng)各階IMF分量的角頻率;ci為趨勢項各項系數(shù);ut為汾河水庫左壩岸滲流量剩余值序列;φk為時間序列ARI模型參數(shù);ut-k為t-k時刻的剩余值序列;εt為模型殘差項。
基于SAS軟件平臺運用Gauss-Newton[9](高斯-牛頓)法對全部參數(shù)進行重構(gòu),所得EMD-ARI模型參數(shù)重構(gòu)的最小二乘估計值如表2、3、4所示。
表2 IMF分量傅氏系數(shù)
表3 趨勢項系數(shù)
表4 ARI模型參數(shù)
汾河水庫EMD-ARI模型的回歸方程為
水庫左壩岸滲流序列只進行時序分析,經(jīng)平穩(wěn)性與隨機性檢驗后確定最優(yōu)模型為ARI(2,1),對模型參數(shù)進行條件最小二乘檢驗,檢驗結(jié)果如表5所示。
表5 條件最小二乘結(jié)果
由表5可知,t檢驗統(tǒng)計量的P(Pr>|t|)值全部小于檢驗水平0.01,模型參數(shù)均極顯著,將一階差分還原后得ARI(2,1)模型的表達式為
式中,yt為汾河水庫左壩岸滲流量序列;yt-1為t-1時刻滲流量序列;yt-2為 t-2時刻滲流量序列;yt-3為t-3時刻滲流量序列;εt為模型殘差項。
經(jīng)參數(shù)檢驗后對ARI(2,1)模型進行殘差檢驗,檢驗結(jié)果如表6所示。
從表6可以看出,ARI(2,1)模型殘差的LB檢驗統(tǒng)計量的P(Pr>ChiSq)值均大于0.05,殘差序列為白噪聲,表明滲流序列可建立ARI(2,1)模型。
表6 ARI(2,1)模型殘差自相關(guān)檢驗
表7 滲流模型預(yù)報誤差分析
圖2 EMD-ARI滲流模型擬合
圖3 ARI滲流模型擬合
(1)EMD-ARI模型與ARI模型的擬合結(jié)果分別如圖2、3所示。從圖2和圖3可以看出,兩個模型的預(yù)報值與實測值的吻合度均較高、擬合效果均較好。
(2)EMD-ARI模型與ARI模型的預(yù)報結(jié)果如表7所示。由表7可知,①EMD-ARI模型其預(yù)報值與實測值的相對誤差E總體偏??;②EMD-ARI模型第6步預(yù)測的相對誤差E為7.56%,小于ARI模型的11.65%;④E<10時EMD-ARI模型的預(yù)報合格率為100%,高于ARI模型的83%;④E<15,E<20時兩模型的預(yù)報合格率均為100%,滿足相對誤差E小于20%的要求。
(1)汾河水庫左壩岸滲流序列經(jīng)EMD分解,得4階IMF分量與1階趨勢分量,4階IMF分量表征滲流序列可能具有的4個周期成分,1階趨勢分量表明滲流量隨時間推移增加。
(2)EMD-ARI模型與ARI模型的相對誤差呈先增后減再突變增趨勢,即相對誤差表現(xiàn)為跳躍特征,且第6步預(yù)測的相對誤差的突變增表明隨預(yù)測步長的增大模型預(yù)測精度降低。
(3)本文1步預(yù)測的時間間隔為7天,ARI模型第6步預(yù)測的相對誤差為11.65%,誤差較大,因此適用于進行5步預(yù)測,即預(yù)測汾河水庫未來35天以內(nèi)的左壩岸滲流量;EMD-ARI模型第6步預(yù)測的相對誤差為7.56%,預(yù)測精度較高,可對汾河水庫左壩岸滲流量進行6步預(yù)測。
(4)EMD-ARI滲流耦合模型的預(yù)報精度高于ARI滲流單一模型的預(yù)報精度,該法切實可行,可用于滲流預(yù)測。
[1]陳文燕,朱林,王文韜.大壩安全監(jiān)測的現(xiàn)狀與發(fā)展趨勢[J].電力環(huán)境保護,2009,25(6):38-42.
[2]毛昶熙,段祥寶,李祖貽.滲流數(shù)值計算分析與程序應(yīng)用[M].南京:河海大學(xué)出版社,1999.
[3]趙雪花,安莉莉,袁旭琦.基于HHT和R/S分析的黃河上游年徑流序列演變模式分析[J].水電能源科學(xué),2013,31(7): 9-12.
[4]王燕.應(yīng)用時間序列分析[M].北京:中國人民大學(xué)出版社,2005.
[5]王振龍.應(yīng)用時間序列分析(第二版)[M].北京:中國統(tǒng)計出版社,2010.
[6]BOWERMAN B L,O’CONNELL R T.Forecasting and Time Seriers[M].3rd ed.Brooks/Cole.Duxbury,United States America,1993.
[7]肖枝洪,郭明月.時間序列分析與SAS應(yīng)用 [M].武漢:武漢大學(xué)出版社,2009.
[8]楊金芳,翟永杰,王東風(fēng),等.基于支持向量回歸的時間序列預(yù)測[J].中國電機工程學(xué)報,2005,25(17): 110-114.
[9]高惠璇,等.SAS系統(tǒng) SAS/STAT軟件使用手冊 [M].北京:中國統(tǒng)計出版社,1997.