基于Kalman-ARIMA模型的大壩變形預(yù)測

薛 洋，楊 光，許 雷(.河海大學(xué)水利水電學(xué)院，南京 20098；2.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，南京 20098)

變形是水工建筑物運(yùn)行狀態(tài)的重要效應(yīng)量，也是最主要的監(jiān)測量[1,2]，變形能夠直觀可靠地反映大壩的結(jié)構(gòu)性態(tài)和安全狀況[3]，因此，變形預(yù)測是大壩安全監(jiān)測系統(tǒng)的關(guān)鍵組成部分，也是防汛調(diào)度的關(guān)鍵，對及時掌握大壩變形狀態(tài)、確保大壩安全運(yùn)行具有十分重要的作用[4,5]。在壩工實際問題中，影響大壩變形的因素極為復(fù)雜，除庫水壓力(水位)外，還有溫度、降雨、滲流、施工、地基、周圍環(huán)境和時效等[1,2]多種非線性因素，從而使得準(zhǔn)確有效地對大壩變形進(jìn)行預(yù)測具有一定的難度。目前，常用的大壩變形預(yù)測方法有多元線性回歸分析法、有限元方法、灰色GM模型、模糊聚類分析法、遺傳算法和時間序列分析法等。其中，基于三大經(jīng)典監(jiān)控模型(統(tǒng)計模型、確定性模型及混合模型)的變形預(yù)報與評價應(yīng)用得最為廣泛，但其存在參數(shù)確定困難、所建數(shù)學(xué)模型較為復(fù)雜、要求線性方程組非奇異與無病態(tài)且與大壩實際變形狀況相吻合[6]等缺點，難以滿足大壩變形預(yù)測的精度要求。

大壩變形監(jiān)測數(shù)據(jù)序列一般呈現(xiàn)出顯著的非平穩(wěn)性和周期性，并且在監(jiān)測過程中，易受自動化監(jiān)測資料噪聲污染。為提高大壩變形預(yù)測精度，及時準(zhǔn)確地反映大壩變形狀態(tài)，本文提出Kalman-ARIMA模型，即先利用卡爾曼濾波法對變形監(jiān)測數(shù)據(jù)序列進(jìn)行降噪處理，剔除隨機(jī)干擾噪聲(包括干擾噪聲和動態(tài)噪聲)，使得監(jiān)測數(shù)據(jù)能夠更準(zhǔn)確地反映大壩變形的實際狀況，然后利用自回歸差分滑動平均混合模型(ARIMA模型)對濾波后的數(shù)據(jù)進(jìn)行大壩變形預(yù)測，從而達(dá)到提高能大壩變形預(yù)測精度的目的。

1 卡爾曼濾波法

卡爾曼濾波是一種從與被提取信號有關(guān)的觀測量中通過算法估計出所需信號的濾波算法[7,8]。目前，該方法廣泛應(yīng)用于動態(tài)數(shù)據(jù)處理過程中，它具有最小無偏差性，并且能夠最大限度地剔除數(shù)據(jù)中包含的隨機(jī)干擾誤差即噪聲污染，得到真實值的最優(yōu)估計[9]。大壩變形監(jiān)測序列可以看作是一個隨機(jī)離散系統(tǒng)，隨機(jī)線性離散系統(tǒng)的卡爾曼濾波方程表達(dá)式為：

xk+1=Ak+1 | kxk+wk

(1)

zk+1=Ck+1xk+1+vk+1

(2)

式中：xk為tk時刻狀態(tài)向量；zk+1為tk+1的觀測向量；Ak+1|k為tk時刻至tk+1時刻的狀態(tài)轉(zhuǎn)移矩陣；Ck+1為tk+1時刻的觀測矩陣；wk為系統(tǒng)噪聲向量；vk+1為觀測噪聲向量。

系統(tǒng)噪聲向量wk和觀測噪聲向量vk的統(tǒng)計特性滿足以下條件：

(3)

式中：Q(k)為系統(tǒng)噪聲向量的方差矩陣；R(k)為觀測噪聲向量的方差矩陣；δkj為狄克拉函數(shù)，當(dāng)k=j時，δij=1，當(dāng)k≠j時，δkj=0。

(4)

Pk=[I-KkCk]Pk|k-1

(5)

其中：

(6)

Pk|k-1=Ak|k-1Pk-1ATk|k-1+Qk-1

(7)

Kk=Pk|k-1CTk[CkPk|k-1Ck+Rk]-1

(8)

式(4)～(8)即為隨機(jī)線性離散系統(tǒng)卡爾曼濾波的基本方式。

2 ARIMA模型

ARIMA模型又稱博克斯-詹金斯模型，是一種時間序列預(yù)測方法。近些年，該模型廣泛應(yīng)用于大壩變形預(yù)測領(lǐng)域，它不僅可以考慮變形監(jiān)測數(shù)據(jù)在時間序列上的依存性，還可以考慮數(shù)據(jù)受到隨機(jī)干擾的波動性，從而能夠較為準(zhǔn)確地進(jìn)行大壩變形的預(yù)測。時間序列基本模型包括自回歸模型(AR)、移動平均模型(MA)和自回歸滑動平均模型(ARMA)[10]，這3種模型都是針對平穩(wěn)時間序列而言的，而大壩變形監(jiān)測數(shù)據(jù)一般都是非平穩(wěn)的，包含長期趨勢、季節(jié)波動、循環(huán)波動和不規(guī)則變動，故需要對原數(shù)據(jù)進(jìn)行有限次差分處理后產(chǎn)生一個平穩(wěn)的時間序列，再進(jìn)行建模。ARIMA模型正是基于此思想對非平穩(wěn)數(shù)據(jù)進(jìn)行處理與建模的，ARIMA(p,d,q)中，AR代表自回歸模型，MA代表滑動平均模型，I表示差分，p為自回歸階數(shù)，d為差分次數(shù)，q為滑動平均階數(shù)，其數(shù)學(xué)表達(dá)式為：

Φ(B)▽dxt=Θ(B)εt

(9)

其中：

Φ(B)=1-Φ1B-Φ2B2-…-ΦpBp

(10)

Θ(B)=1-Θ1B-Θ2B2-…-ΘpBp

(11)

▽=1-B

(12)

式中：xt(t=1,2,3,…)為時間序列；εt為均值是0，方差是σε的正態(tài)白噪聲過程；Φi(i=1,2,…,p)和Θj(j=1,2,…,q)為模型的待估系數(shù)；B為后移差分算子。

3 大壩變形預(yù)測的Kalman-ARIMA模型

3.1 卡爾曼濾波處理分析

大壩變形監(jiān)測數(shù)據(jù)屬于隨機(jī)線性離散系統(tǒng)，本文采用(α-β)卡爾曼濾波模型[11]，即將監(jiān)測點的變形和變形速度看成狀態(tài)參數(shù)，由于大壩變形量比較小，其加速度也很小，故將變形加速度看成動態(tài)噪聲向量，則狀態(tài)方程和觀測方程為：

(14)

3.2 時間序列分析與建模

時間序列模型常用的建模方法有Box-Jenkins法、Pandit-Wu S.M法(也稱為DDS法)以及長自回歸、白噪化法等，針對大壩變形監(jiān)測數(shù)據(jù)，可采用Box-Jenkins法建立ARIMA模型[12,13]，其步驟如下：

(1)平穩(wěn)化。運(yùn)用差分法或?qū)?shù)據(jù)先進(jìn)行對數(shù)變換再進(jìn)行差分運(yùn)算，提取大壩變形監(jiān)測數(shù)據(jù)序列中的確定性信息，使其自相關(guān)函數(shù)呈現(xiàn)指數(shù)衰減(即拖尾)的特征以達(dá)到平穩(wěn)的目的。

(2)識別并估計模型，包括模型識別、模型定階及參數(shù)估計。模型識別是在一定標(biāo)準(zhǔn)下(如序列是零均值過程)，根據(jù)樣本自相關(guān)函數(shù)(ACF)和樣本偏自相關(guān)函數(shù)(PACF)的形式，結(jié)合殘差方差圖、最佳準(zhǔn)則函數(shù)(AIC準(zhǔn)則、SBC準(zhǔn)則和BIC準(zhǔn)則)等定階法，確定ARMA模型的p、q值。然后利用矩估計、最小二乘估計或極大似然估計法對所建ARMA模型進(jìn)行參數(shù)估計。

(3)診斷，包括模型的適應(yīng)性檢驗和參數(shù)的顯著性檢驗。模型的適應(yīng)性檢驗的實質(zhì)為檢驗被估模型殘差的獨(dú)立性。若殘差序列為白噪聲序列，說明擬合模型準(zhǔn)確刻畫了時間序列過程，顯著有效；若殘差序列不是白噪聲序列，說明擬合模型不顯著，相關(guān)信息有殘留，則返回步驟(2)，重新估計p、q值并建模，再進(jìn)行參數(shù)檢驗。參數(shù)的顯著性檢驗，即檢驗?zāi)Ｐ椭忻恳粋€未知參數(shù)是否顯著為零。當(dāng)變量的參數(shù)不顯著時，將該變量從模型中剔除，重新進(jìn)行(2)。

(4)預(yù)測。通過診斷檢驗的模型即可用于大壩變形預(yù)測。

診斷這一步驟極為必要，樣本的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)只對時間序列ARMA(p，q)過程的識別起指導(dǎo)作用，尤其當(dāng)p，q都不為0時，識別難度較大，還需借助系統(tǒng)的物理背景和先驗知識及實驗技巧和經(jīng)驗。

4 工程實例分析

某水電站位于皖南長江支流青弋江上游，是一座綜合性中型水利水電樞紐工程。樞紐包括大壩、溢洪道、泄洪中孔、泄水底孔、發(fā)電廠房和筏道。大壩是一座同心圓變半徑的混凝土重力拱壩，壩頂高程為126.3 m，最大壩高為76.3 m，壩頂弧長419 m，壩頂寬8 m，最大壩底寬53.5 m。經(jīng)過一次大壩安全監(jiān)測系統(tǒng)更新改造后，該工程具備了一套以變形觀測為主的大壩觀測系統(tǒng)。大壩共布置了16個倒垂線測點用來監(jiān)測壩基和F11、F32大斷層及兩壩肩巖體的變形。本文選取處于8號壩段的一個測點從2013年1月1日到2013年5月27日的位移觀測值進(jìn)行分析，共計147個位移值。其中將前137個實測值作為已知數(shù)據(jù)，后10個實測值作為未來實際位移值。首先對已知數(shù)據(jù)進(jìn)行卡爾曼濾波處理，再對濾波后的數(shù)據(jù)建立ARIMA模型，并進(jìn)行未來10個位移的預(yù)測，然后把位移預(yù)測結(jié)果與單純用ARIMA模型對原數(shù)據(jù)進(jìn)行預(yù)測的結(jié)果作對比，從而驗證本文提出的Kalman-ARIMA模型在大壩變形預(yù)測上的合理性和適用性。

利用3.1節(jié)所述的針對大壩變形數(shù)據(jù)的卡爾曼濾波法對已知位移數(shù)據(jù)進(jìn)行濾波處理，結(jié)果如圖1所示，并統(tǒng)計已知位移值與濾波值之間的殘差，發(fā)現(xiàn)絕大多數(shù)殘差值在±0.01 mm以內(nèi)，而且濾波值反應(yīng)的變化規(guī)律與實測位移值基本一致，說明該卡爾曼濾波過程是對原始觀測位移數(shù)據(jù)進(jìn)行了一個很好的平滑處理，剔除噪聲的效果較好[14]。

圖1 位移實測值與卡爾曼濾波值對比Fig.1 Comparison of the observed and kalman filter values

根據(jù)3.2節(jié)所述ARIMA模型的建立方法對位移實測值進(jìn)行建模，模型參數(shù)p,d,q分別為0,1,1。對卡爾曼濾波后的位移數(shù)據(jù)進(jìn)行建模，模型參數(shù)p,d,q分別為0,1,0。利用ARIMA(0,1,1)模型對實測數(shù)據(jù)進(jìn)行預(yù)測，利用ARIMA(0,1,0)模型對卡爾曼濾波后數(shù)據(jù)進(jìn)行預(yù)測，結(jié)果見表1和圖2?？梢钥闯?，對位移實測值進(jìn)行卡爾曼濾波處理，然后再進(jìn)行ARIMA模型的建立并預(yù)測，可以有效地減小預(yù)測值與實際原位移之間的殘差絕對值和殘差標(biāo)準(zhǔn)差，從而提高大壩位移的預(yù)測精度，產(chǎn)生這個效果的主要原因可以分析為：卡爾曼濾波過程可以剔除觀測過程中隨機(jī)干擾噪聲對監(jiān)測數(shù)據(jù)的影響，使得濾波后的數(shù)據(jù)更加逼近真實值，對大壩較為真實的位移值進(jìn)行建?？梢暂^為準(zhǔn)確地反應(yīng)壩體力學(xué)機(jī)理對位移的影響，從而可以更加精確地進(jìn)行變形預(yù)測。由此說明，Kalman-ARIMA模型應(yīng)用在大壩變形預(yù)測上是合理有效的。

表1 實測數(shù)據(jù)預(yù)測值與卡爾曼濾波后數(shù)據(jù)預(yù)測值對比Tab.1 Comparison of the predictive values between observed and kalman filter data

圖2 兩種方法的殘差絕對值對比Fig.2 Comparison of the absolute residuals from the two methods

5 結(jié) 語

卡爾曼濾波是一種對動態(tài)觀測數(shù)據(jù)進(jìn)行實時處理的有效辦法，可以剔除觀測過程中隨機(jī)干擾噪聲對數(shù)據(jù)精度的影響，從而可以獲得更加接近真實值的觀測數(shù)據(jù)。本文提出卡爾曼濾波與ARIMA模型相結(jié)合的預(yù)測方法，即利用卡爾曼濾波對大壩原始變形數(shù)據(jù)進(jìn)行濾波處理，再根據(jù)對濾波后的數(shù)據(jù)進(jìn)行建模計算并預(yù)測。結(jié)合具體大壩位移觀測數(shù)據(jù)分析，結(jié)果表明，Kalman-ARIMA模型能夠有效地剔除大 壩位移數(shù)據(jù)中的隨機(jī)干擾噪聲，并得到比較精確的預(yù)測結(jié)果，因此，該模型在大壩變形預(yù)測方面具有推廣價值。

□

[1] 吳中如．水工建筑物安全監(jiān)控理論及其應(yīng)用[M]．南京：河海大學(xué)出版社，2003．

[2] 顧沖時，吳中如．大壩與壩基安全監(jiān)控理論和方法及其應(yīng)用[M]．南京：河海大學(xué)出版社，2006．

[3] 何金平．大壩安全監(jiān)測理論與應(yīng)用[M]．北京：中國水利水電出版社，2010．

[4] 沈細(xì)中，張文鴿，馮夏庭．大壩變形預(yù)測的ANFIS 模型[J]．巖土力學(xué)，2006，27(增)：119-122．

[5] 張 豪，許四法．基于經(jīng)驗?zāi)B(tài)分解和遺傳支持向量機(jī)的多尺度大壩變形預(yù)測[J]．巖石力學(xué)與工程學(xué)報，2011，30(增2)：3 681-3 688．

[6] 李守巨，劉迎曦，劉玉靜．基于進(jìn)化神經(jīng)網(wǎng)絡(luò)混凝土大壩變形預(yù)測[J]．巖土力學(xué)，2003，24(4)：635-638．

[7] 張志鵬，王偉平，鄭海超．卡爾曼濾波及其在時間序列預(yù)測中的應(yīng)用[J]．儀表技術(shù)，2010,(7)：37-39．

[8] 顏可珍，張 鄒，林 峰．卡爾曼時間序列模型預(yù)測路面平整度指數(shù)[J]．路基工程，2012,(4)：1-3．

[9] Sinopoli B, Schenato L, Franceschetti M, et al. Kalman filtering with intermittent observations[J]. Automatic Control, IEEE Transactions on, 2004,49(9):1 453-1 464.

[10] 潘迪夫，劉 輝，李燕飛．基于時間序列分析和卡爾曼濾波算法的風(fēng)電場風(fēng)速預(yù)測優(yōu)化模型[J]．電網(wǎng)技術(shù)，2008，32(7)： 82-86．

[11] 王 利，李亞紅，劉萬林．卡爾曼濾波在大壩動態(tài)變形監(jiān)測數(shù)據(jù)處理中的應(yīng)用[J]．西安科技大學(xué)學(xué)報, 2006,26(3)：353-357．

[12] 張樹京，齊立心．時間序列分析簡明教程[M]．北京：清華大學(xué)出版社，2003．

[13] 肖枝洪，郭明月．時間序列分析與SAS應(yīng)用[M]．武漢：武漢大學(xué)出版社，2012．

[14] 王井利，張春哲．自適應(yīng)卡爾曼濾波在地鐵監(jiān)測中的應(yīng)用[J]．沈陽建筑大學(xué)學(xué)報：自然科學(xué)版，2014,(2)：263-268．