基于多重分形原理的大壩位移分析與預(yù)測研究

謝雨航，鄧念武，劉勇軍

(1.武漢大學(xué)水利水電學(xué)院，湖北 武漢 430072；2.中國長江三峽集團有限公司，湖北 宜昌 443133)

作為一門新興的理論體系，分形理論有獨特的魅力，目前已經(jīng)被證實該理論可以較好地應(yīng)用到自然科學(xué)以及社會科學(xué)中[1]。用分形理論分析雜亂無章的復(fù)雜系統(tǒng)，可發(fā)現(xiàn)它們?nèi)匀痪哂袛?shù)學(xué)規(guī)律，分形理論為時間序列分析供了新的分析途徑。大壩安全監(jiān)測數(shù)據(jù)屬于時間序列，監(jiān)測數(shù)據(jù)也具有隨機性、多尺度變化等復(fù)雜的特性，分形理論可以有效識別時間序列包含的內(nèi)在規(guī)律。本文將從多重分形去趨勢波動分析、非對稱多重分形去趨勢波動互相關(guān)分析及分形預(yù)測等方面對多重分形進行研究，以分析非平穩(wěn)時間序列，評估兩時間序列的相關(guān)性并預(yù)測大壩位移發(fā)展趨勢。

1 多重分形理論

1.1 多重分形去趨勢波動分析(MF-DFA)原理

Kantelhardt J W和Zschiegner S A等學(xué)者結(jié)合去趨勢波動分析(Detrended Fluctuation Analysis，簡稱DFA)方法與多重分形理論(Multifractal theory，簡稱MF)[2]，改進性地提出了MF-DFA方法。該方法統(tǒng)計時間序列在每個分割區(qū)間上波動的平均值，計算時間序列波動函數(shù)的冪律性和廣義Hurst指數(shù)，運用Hurst指數(shù)衡量時間序列結(jié)構(gòu)的平穩(wěn)和非平穩(wěn)性，并使用多重分形奇異譜衡量時間序列波動的奇異性，能準確分析非平穩(wěn)時間序列的長程相關(guān)性，并能夠避免對時間序列相關(guān)性的誤判。MF-DFA對非平穩(wěn)時間序列的多重分形特性分析更具優(yōu)越性。

1)構(gòu)造新的時間序列。假設(shè)有一個時間序列xi(i=1,2,3,…,n)，現(xiàn)構(gòu)造新的時間序列Xi：

(1)

2)分割時間序列。將時間序列Xi以相同的長度s(5≤s≤n/4)分割成不重疊的ns(ns=int(n/s))段數(shù)據(jù)，為充分利用數(shù)據(jù)的長度，再從數(shù)據(jù)的反方向用相同的長度分段，得到2ns段數(shù)據(jù)。

3)利用最小二乘擬合計算序列2ns個區(qū)間中每個區(qū)間的局部趨勢：

(2)

則整個序列的q階波動函數(shù)Fq(s)為

(3)

(4)

4)計算廣義Hurst指數(shù)h(q)、標度指數(shù)τ(q)、奇異譜f(α)與奇異指數(shù)α：

Fq(s)∝sh(q)

(5)

對ln[Fq(s)]做ln(s)的一次線性擬合關(guān)系圖，一次項系數(shù)即q為階廣義Hurst指數(shù)h(q)。

h(q)與標度指數(shù)τ(q)有關(guān)，其解析關(guān)系試如下：

τ(q)=qh(q)-1

(6)

對于描述多重分形序列的方法是奇異譜f(α)，可利用勒讓德(Legendre)轉(zhuǎn)換得到f(α)與τ(q)之間的關(guān)系：

α=τ′(q)=h(q)+qh′(q)

(7)

f(α)=q[α-h(q)]+1

(8)

奇異指數(shù)α用來描述序列的尖頂、峰值、波動和頻率等奇異性特征[3]。

1.2 非對稱多重分形去趨勢波動互相關(guān)(AMF-DCCA)分析

在復(fù)雜的系統(tǒng)中，常有效應(yīng)量序列受到自變量序列的影響而發(fā)生變化，而MF-DFA方法只能反映時間序列的自相關(guān)特性，無法體現(xiàn)效應(yīng)量序列與自變量序列之間的聯(lián)系。通過在效應(yīng)量序列中引入非對稱的思想，提出了非對稱多重分形去趨勢波動互相關(guān)分析(AMF-DCCA)方法[4]，評估兩時間序列是否存在非對稱相關(guān)性，算法如下：

1)假設(shè)有兩個長度均為n的原型觀測序列xi，yi(i=1,2,3,…,n)，構(gòu)造新的時間序列Xi，Yi：

(9)

(10)

以原型觀測序列xi為例，q階平均波動函數(shù)為：

(11)

(12)

3)對比MF-DFA方法，得出平均波動趨勢的波動函數(shù)為：

(13)

4)通過q階波動函數(shù)Fq(s)與時間尺度s之間存在的冪律關(guān)系可以計算廣義Hurst指數(shù)h1(q)：

(14)

式中，h1(q)為廣義Hurst指數(shù)，但h1(q)表示的是兩個相關(guān)序列之間的互相關(guān)性。

1.3 分形預(yù)測理論

根據(jù)分形應(yīng)用中的數(shù)學(xué)基礎(chǔ)與方法[5-6]中的定義，分形維數(shù)可近似地表示為

N=C·r-D

(15)

式中：r為特征線度，如長度和時間等；N為r對應(yīng)的量，如大壩安全監(jiān)測中的位移、溫度和水位等；C為待定常數(shù)；D為分維數(shù)。

分維數(shù)D為常量時的分形分布稱為常維分形，(ln(Ni),ln(ri))呈線性關(guān)系，根據(jù)該直線上的任意兩個數(shù)據(jù)點可以確定該段直線的分形參數(shù)分維數(shù)D和常數(shù)C：

D=ln(Ni/NJ)/ln(rj/ri)i≠j

(16)

C=Ni·rD

(17)

設(shè)N與r之間的函數(shù)關(guān)系為N=f(r)，對于常維分形分布有f(r)=C/rD，求解出r得到任一函數(shù)關(guān)系的常維分形形式：

r=[C/f(r)]1/D

(18)

在實際應(yīng)用中函數(shù)關(guān)系通常是未知的，只存在一系列的數(shù)據(jù)點(Ni,ri)，因此需要尋找恰當(dāng)?shù)淖儞Q方法，以求出變換的具體形式。運用施行累計和的系列變換r和N，將數(shù)據(jù)進行一系列的變換，從中選出一種變換，讓變換后的數(shù)據(jù)能夠用分形分布來處理。變換方式如下：

由于ln(Ni)，ln(ri)(i=1,2,…,n)在一般情況下不能與分形分布模型符合良好，于是運用公式(1)將Ni構(gòu)造一階累計和序列S1i：

(19)

同樣可以構(gòu)造二階、三階及多階累計和序列S(J)i(J為階數(shù))。

建立各階累計和序列的變維分形模型，通過公式(16)求得各階的分形參數(shù)分維數(shù)D；比較各階分段變維分形模型，并選擇效果最好的變換，運用外插的方法開展預(yù)測計算工作。

2 實例分析

福建某雙曲拱壩，壩頂高程346.5 m，最大壩高81.5 m，壩頂中心線弧長153.676 m，壩頂設(shè)6個水平位移監(jiān)測點，如圖1所示。本文采用2007年1月至2015年12月共9年大壩徑向位移和環(huán)境量數(shù)據(jù)進行MF-DFA分析、AMF-DCCA分析以及分形預(yù)測。

2.1 大壩位移時間序列多重分形去趨勢波動分析

根據(jù)大壩變形監(jiān)測資料的周期性波動趨勢，結(jié)合時間尺度的合理取值，取時間尺度s=6,9,12,15,18,21,24，單位為月，根據(jù)MF-DFA算法得到各時間序列的h(q)～q關(guān)系圖和f(α)～α關(guān)系圖。

圖1 大壩水平位移監(jiān)測點位布置圖

圖2與圖3分別為溫度、水位和各測點徑向位移廣義Hurst指數(shù)h(q)～q關(guān)系。當(dāng)q=2時，水位和溫度時間序列的h(2)值分別為0.69和0.61，各觀測點時間序列的h(2)值均大于0.5，說明水位與溫度和各測點時間序列存在長程相關(guān)性。整體上看，溫度序列的大或小波動標度變化都較水位的大或小波動標度變化顯著；由于大壩處于穩(wěn)定的運行期，可以看出各測點序列的小波動標度單調(diào)遞減較各測點的大波動變化明顯，說明各測點出現(xiàn)小位移的概率較高。

圖2 環(huán)境量時間序列關(guān)系曲線圖

圖3 各測點時間序列關(guān)系曲線圖

圖4為各測點時間序列多重分形奇異譜曲線。位移觀測點A、B、C、D、E和F的分形奇異譜寬度Δα分別為1.279、0.523、1.149、0.573、1.150和1.026，奇異譜高度Δf(α)分別為-0.035、0.524、0.558、1.010、0.449和0.074。觀測點C和D時間序列奇異譜高度Δf(α)>0，且明顯大于其他測點時間序列奇異譜高度，說明該兩測點位移時間序列中大位移出現(xiàn)的概率較高。由圖1可知，觀測點C和D位于拱壩的中間位置，當(dāng)環(huán)境量變化時位移變化會較大，拱壩中間部位位移量較大，符合拱壩的變形規(guī)律。

圖4 各測點時間序列多重分形奇異譜曲線圖

2.2 大壩位移時間序列多重分形去趨勢波動互相關(guān)分析

大壩位移監(jiān)測序列的AMF-DCCA分析旨在分析大壩位移時間序列與環(huán)境量時間序列之間的相關(guān)程度，環(huán)境量的變化直接影響到大壩位移的變化，由于篇幅有限，選取A、C測點的位移時間序列為主要研究對象，驗證和分析大壩位移時間序列與環(huán)境量時間序列的互相關(guān)特性。由圖5各測點的h(q)～q函數(shù)關(guān)系圖可知，隨著q值的增大，h(q)單調(diào)遞減，說明各測點位移時間序列與環(huán)境量之間存在互相關(guān)性，大壩位移時間序列的波動受到環(huán)境量時間序列波動的影響，符合大壩的變形規(guī)律。由圖中可以看出，兩測點h(2)均大于0.5且小于1，說明大壩位移時間序列與環(huán)境量時間序列存在長程相關(guān)性，即環(huán)境量時間序列的改變會對大壩位移時間序列變化產(chǎn)生影響。當(dāng)q<0時，各測點環(huán)境量的h(q)變幅較q>0的大，說明環(huán)境量小幅度波動對大壩位移時間序列的互相關(guān)特性影響較大幅度波動的大，尤其溫度的小幅度波動最為顯著，即溫度時間序列的小幅度波動對大壩位移時間序列的波動影響較大。

圖5 測點位移時間序列與環(huán)境量間的h(q)～q關(guān)系曲線圖

2.3 大壩位移時間序列多重分形原因判定

為了運用分形預(yù)測理論對大壩變形開展預(yù)測工作，首先對影響大壩位移監(jiān)測時間序列多重分形特性形成的原因進行分析。由于篇幅有限，僅選取測點A和C進行分析。

時間序列表現(xiàn)出多重分形特性的原因主要來源于兩個方面：①時間序列的非正態(tài)分布；②時間序列的長程相關(guān)性。對原始時間序列進行隨機重排得到重排序列，相位隨機化處理原始時間序列得到替代序列，根據(jù)原始序列、重排序列和替代序列可以有效區(qū)分上述兩種原因?qū)Χ嘀胤中蔚呢暙I程度。

現(xiàn)對原始時間序列和重排后序列進行MF-DFA分析。根據(jù)MF-DFA算法取時間尺度s=[6,9,12,15,18,21,24]，取波動函數(shù)階數(shù)q=[-10,-8,-6,-4,-2,0,2,4,6,8,10]，最小二乘擬合階數(shù)m=1，得到各位移測點的原始序列、重排序列和替代序列的h(q)～q關(guān)系如圖6所示。

由圖6的兩測點位移時間序列h(q)～q關(guān)系曲線圖可知，重排序列的非線性特征明顯弱于原始序列，說明原始序列和替代序列的多重分形特性較強。即大壩位移時間序列的多重分形特性主要是長程相關(guān)性引起的。

圖6 測點位移時間序列h(q)～q關(guān)系曲線圖

圖7為測點A、C的q～ln(s)～ln(Fq(s))三維函數(shù)關(guān)系圖。當(dāng)q值由-10至10之間變化時，在時間尺度s較小時，波動函數(shù)值之間的差異明顯，而在時間尺度逐漸變大時，波動函數(shù)值之間的差異逐漸變小，說明小的時間尺度能夠識別局部時間內(nèi)的大小波動，大的時間尺度在一定程度上均化了大小波動的幅度。因此在做大壩變形預(yù)測時，選擇時間尺度s=12，即以年周期進行預(yù)測效果較好。

圖7 q～ln(s)～ln(Fq(s))三維函數(shù)關(guān)系圖

2.4 大壩位移時間序列分形理論預(yù)測分析

以測點A、C從2007年1月21日至2015年12月18日的測值序列為例，按照上述方法，計算1～3階的分段維數(shù)值，并作出兩測點的分段維數(shù)圖。

由圖8與圖9可看出，一階分段維數(shù)圖波動范圍最小，尾部也較平直，可以基于一階分段維數(shù)運用2007年1月21日至2014年12月23日期間序列建立模型，預(yù)測2015年1月23日和2015年12月18日期間的測值。測點A和C的實測值、預(yù)測值和殘差如表1所示。

圖8 A測點平移處理時間序列分段維數(shù)圖

圖9 C測點平移處理時間序列分段維數(shù)圖

表1 測點A、C實測值、預(yù)測值和殘差統(tǒng)計表 mm

3 結(jié) 語

本文利用多重分形原理對大壩位移進行分析研究，并進行了實例分析。MF-DFA分析表明，該大壩壩頂各測點出現(xiàn)小位移的概率較高，且拱冠出現(xiàn)大位移的概率大；AMF-DCCA分析表明，溫度變化對大壩位移時間序列的波動影響較大；運用分形預(yù)測理論，對各測點徑向位移時間序列進行了擬合和預(yù)測，結(jié)果表明分形預(yù)測理論可以很好地應(yīng)用到大壩位移分析中。相信隨著分形預(yù)測理論的不斷發(fā)展，分形預(yù)測理論對大壩安全監(jiān)測資料的分析和研究將會更加深入。