謝雨航，鄧念武，劉勇軍

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

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

1 多重分形理論

1.1 多重分形去趨勢(shì)波動(dòng)分析(MF-DFA)原理

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

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

(1)

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

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

(2)

則整個(gè)序列的q階波動(dòng)函數(shù)Fq(s)為

(3)

(4)

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

Fq(s)∝sh(q)

(5)

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

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

τ(q)=qh(q)-1

(6)

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

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

(7)

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

(8)

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

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

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

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

(9)

(10)

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

(11)

(12)

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

(13)

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

(14)

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

1.3 分形預(yù)測(cè)理論

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

N=C·r-D

(15)

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

分維數(shù)D為常量時(shí)的分形分布稱為常維分形，(ln(Ni),ln(ri))呈線性關(guān)系，根據(jù)該直線上的任意兩個(gè)數(shù)據(jù)點(diǎn)可以確定該段直線的分形參數(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)，對(duì)于常維分形分布有f(r)=C/rD，求解出r得到任一函數(shù)關(guān)系的常維分形形式：

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

(18)

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

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

(19)

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

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

2 實(shí)例分析

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

2.1 大壩位移時(shí)間序列多重分形去趨勢(shì)波動(dòng)分析

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

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

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

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

圖3 各測(cè)點(diǎn)時(shí)間序列關(guān)系曲線圖

圖4為各測(cè)點(diǎn)時(shí)間序列多重分形奇異譜曲線。位移觀測(cè)點(diǎn)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è)點(diǎn)C和D時(shí)間序列奇異譜高度Δf(α)>0，且明顯大于其他測(cè)點(diǎn)時(shí)間序列奇異譜高度，說(shuō)明該兩測(cè)點(diǎn)位移時(shí)間序列中大位移出現(xiàn)的概率較高。由圖1可知，觀測(cè)點(diǎn)C和D位于拱壩的中間位置，當(dāng)環(huán)境量變化時(shí)位移變化會(huì)較大，拱壩中間部位位移量較大，符合拱壩的變形規(guī)律。

圖4 各測(cè)點(diǎn)時(shí)間序列多重分形奇異譜曲線圖

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

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

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

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

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

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

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

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

圖6 測(cè)點(diǎn)位移時(shí)間序列h(q)～q關(guān)系曲線圖

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

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

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

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

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

圖8 A測(cè)點(diǎn)平移處理時(shí)間序列分段維數(shù)圖

圖9 C測(cè)點(diǎn)平移處理時(shí)間序列分段維數(shù)圖

表1 測(cè)點(diǎn)A、C實(shí)測(cè)值、預(yù)測(cè)值和殘差統(tǒng)計(jì)表 mm

3 結(jié) 語(yǔ)

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