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

謝雨航，鄧念武，劉勇軍

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

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

1 多重分形理論

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

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

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

(1)

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

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

(2)

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

(3)

(4)

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

Fq(s)∝sh(q)

(5)

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

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

τ(q)=qh(q)-1

(6)

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

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

(7)

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

(8)

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

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

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

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

(9)

(10)

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

(11)

(12)

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

(13)

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

(14)

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

1.3 分形預測理論

根據分形應用中的數學基礎與方法[5-6]中的定義，分形維數可近似地表示為

N=C·r-D

(15)

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

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

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

(16)

C=Ni·rD

(17)

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

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

(18)

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

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

(19)

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

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

2 實例分析

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

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

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

圖1 大壩水平位移監測點位布置圖

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

圖2 環境量時間序列關系曲線圖

圖3 各測點時間序列關系曲線圖

圖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，且明顯大于其他測點時間序列奇異譜高度，說明該兩測點位移時間序列中大位移出現的概率較高。由圖1可知，觀測點C和D位于拱壩的中間位置，當環境量變化時位移變化會較大，拱壩中間部位位移量較大，符合拱壩的變形規律。

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

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

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

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

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

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

時間序列表現出多重分形特性的原因主要來源于兩個方面：①時間序列的非正態分布；②時間序列的長程相關性。對原始時間序列進行隨機重排得到重排序列，相位隨機化處理原始時間序列得到替代序列，根據原始序列、重排序列和替代序列可以有效區分上述兩種原因對多重分形的貢獻程度。

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

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

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

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

圖7 q～ln(s)～ln(Fq(s))三維函數關系圖

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

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

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

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

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

表1 測點A、C實測值、預測值和殘差統計表 mm

3 結 語

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