基于多重分形的混凝土重力壩水平位移波動分析

2022-02-12 02:34:14周蘭庭柳志坤龔云柱

河海大學學報(自然科學版) 2022年1期

周蘭庭，柳志坤，龔云柱

(1.河海大學水利水電學院, 江蘇南京 210098; 2.青島市經(jīng)濟發(fā)展研究院, 山東青島 266000)

時間序列的研究對于揭示事物發(fā)展變化規(guī)律具有重要的意義，但由于事物的發(fā)展常常受多種因素的影響，導致實際時間序列往往是非平穩(wěn)的，這給其分析研究帶來了困難。為了有效挖掘時間序列的內在規(guī)律，人們引入分形理論來處理非平穩(wěn)時間序列，并在環(huán)境保護、金融市場、微震監(jiān)測等[1-3]領域廣泛應用。推廣到水利領域，混凝土重力壩的水平位移受庫水位漲落、溫度變化等因素[4]的影響，監(jiān)測序列表現(xiàn)出高度的非線性和分形特征，分形理論最初在該領域的應用主要集中在大壩變形性態(tài)的單分形分析[5-6]，通過單分形來表征位移的長程相關性，描述時間序列波動的宏觀特征，但其局限性在于無法對位移過程中某一時刻的局部特征或不同波動層次的復雜精細特征進行刻畫，故引入多重分形來完善單分形理論。

在多重分形領域，Su等[7]應用多重分形去趨勢波動分析法(multifractal detrended fluctuation analysis，MF-DFA)對某混凝土重力壩現(xiàn)場觀測數(shù)據(jù)進行分析，確定現(xiàn)有大壩時間序列的多重分形標度行為，刻畫了大壩的長期行為和結構演變規(guī)律；胡江等[8]利用MF-DFA對某混凝土重力壩監(jiān)測序列的波動特征進行分析，從局部和整體兩個波動層次評價了大壩的工作性態(tài)及其演變規(guī)律；楊景文等[9]基于MF-DFA研究了某混凝土重力壩多個裂縫序列的多重分形特征，并通過序列替換發(fā)現(xiàn)裂縫序列的多重分形特征主要由序列自身的長程相關性造成；韋武昌等[10]利用滑動平均法改進MF-DFA，并將其應用于某混凝土壩的變形分析中，探討了多重分形強度與壩段運行機理的關系，對大壩的工作性態(tài)進行了安全評價；周斌等[11]將MF-DFA延伸到某混凝土重力壩的滲漏量監(jiān)測分析中，通過計算滲漏量的定標指數(shù)確定滲漏監(jiān)測序列的時程相關性，對多重分形在大壩變形分析中的應用有借鑒意義;周蘭庭等[12]利用MF-DFA、MV-MFDFA和A-MFDFA解析了大壩變形性態(tài)的多重分形特征及其對稱性。

前人的研究主要采用傳統(tǒng)方法探討混凝土重力壩監(jiān)測序列自身的多重分形特征，未考慮環(huán)境量對監(jiān)測序列多重分形波動特征的影響，胡江等[8]雖然在研究中引入了水位、溫度等環(huán)境量因素，但并未將其與位移的多重分形特征建立聯(lián)系。本文引入多重分形去趨勢相關性分析法(multifractal detrended cross -correlation analysis，MF-DCCA)和多重分形非對稱去趨勢相關性分析法(multifractal asymmetric detrended cross-correlation analysis，MF-A-DCCA)，在考慮位移序列自身波動特征的同時，兼顧位移與環(huán)境量間的相關性，研究位移在環(huán)境量影響下整體波動與局部趨勢的聯(lián)合多重分形特征。

1 計算方法

1.1 多重分形去趨勢相關性分析法

MF-DCCA[13-14]從多重分形視角考察兩個時序的相關性，在金融、環(huán)境等眾多實證研究中得到廣泛應用。將該方法用于混凝土重力壩水平位移與環(huán)境量的分析中，算法解析如下。

步驟1記大壩實測位移時間序列為x(1)(t)，某一環(huán)境量實測序列為x(2)(t)，序列長度均為N。按式(1)計算相應序列的累計離差和y(i)(t)(i=1，2)：

(1)

步驟2給定時間尺度s，對y(i)(t)進行分割，得到h個長度為s的等長度不重疊連續(xù)子區(qū)間，其中h=int(N/s)。由于N未必能夠整除s，會出現(xiàn)數(shù)據(jù)冗余問題，故采用逆序處理法回收數(shù)據(jù)信息，即在正序劃分的基礎上，從序列的末端開始重復前述操作，總計得到2h個子區(qū)間。

步驟3通過多項式擬合去除各子區(qū)間的趨勢項，按式(2)計算去趨勢后序列間的相關系數(shù)F2(s,v)：

(2)

式中：v為劃分的子區(qū)間；pv(i)(k)為第v個子區(qū)間的趨勢量，i=1，2。

步驟4計算兩序列的q階相關系數(shù)Fq(s)：

(3)

式中：q為波動階數(shù)。當q>0時，表明大波動的影響占優(yōu)勢；當q<0時，表明小波動的影響占優(yōu)勢；特別地，當q=0時，可得極限形式：

(4)

步驟5若序列間存在多重分形相關關系，則有如下冪律關系成立：

Fq(s)∝sH12(q)

(5)

式中H12(q)為序列間的廣義相關Hurst指數(shù)。對式(5)兩邊取對數(shù)得到式(6)，對式(6)最小二乘擬合，其斜率值即為H12(q)。

lgFq(s)=H12(q)lgs+lgC

(6)

式中C為常數(shù)。

觀察擬合得到的H12(q)，若其不隨q的變化而變化，說明序列間的相關關系是單分形的；若其隨q的變化呈非線性變化，說明序列間存在多重分形相關關系。當H12(q)處于不同的區(qū)間范圍時，可得到兩序列間相關關系的不同結論：若H12(q)>0.5，兩個序列間呈現(xiàn)正持續(xù)的相關性；H12(q)<0.5，兩個序列間呈現(xiàn)反持續(xù)的相關性；H12(q)=0.5，則兩個序列間的關系是隨機、不確定的。

步驟6求解兩序列間的相關性多重分形譜：

τ(q)=qH12(q)-1

(7)

α=τ′(q)

(8)

f(α)=qα-τ(q)

(9)

式中：τ(q)為標度函數(shù)；α為奇異性指數(shù)；f(α)為多重分形譜，其表征序列具有相同奇異性指數(shù)子區(qū)間的維數(shù)值。此外，常用分形譜寬度Δα和分形強度ΔH來定量描述位移序列的多重分形特性。Δα主要反映大壩變形的奇異性和性態(tài)演變的空間差異性，ΔH主要反映序列的分形強度，并與參數(shù)值成正比:

Δα=max(α(q))-min(α(q))

(10)

ΔH=max(H(q))-min(H(q))

(11)

MF-DCCA的計算流程如圖1所示。

圖1 MF-DCCA和MF-A-DCCA計算流程Fig.1 Calculation steps of MF-DCCA and MF-A-DCCA

1.2 多重分形非對稱去趨勢相關性分析法

MF-DCCA只能辨識序列間相關性的多重分形特征，但無法進一步刻畫當位移序列處于上升或下降趨勢時，其與環(huán)境量間相關性的多重分形特征，故引入MF-A-DCCA[15-16]來實現(xiàn)這一目的，算法步驟如下。

步驟1同MF-DCCA的步驟1，得到兩序列的累計離差和y(i)(t)，i=1，2。

步驟2選取位移序列x(1)(t)、位移累計離差和y(1)(t)、某一環(huán)境量實測序列累計離差和y(2)(t)，采取MF-DCCA步驟2的方法，得到2h個子區(qū)間。

步驟3利用式(12)(13)對步驟2中的子區(qū)間進行趨勢擬合，L為線性擬合函數(shù)，a為常數(shù)部分，b為擬合斜率，參數(shù)取值范圍：v=1，2,…,2h；i=1，2；k=1，2,…,s。

Lxv(1)(k)=bxv(1)k+axv(1)

(12)

Lyv(i)(k)=bxv(i)k+axv(i)

(13)

步驟4對子區(qū)間進行去趨勢處理，并計算協(xié)方差F2(s，v)：

(14)

步驟5考察式(12)的斜率值bxv(1)，當bxv(1)為正時，表明位移序列在該子區(qū)間為上升趨勢，記為正趨勢，反之，則記為負趨勢。按式(15)(16)計算當位移處于不同趨勢時，位移與環(huán)境量的q階相關系數(shù)：

(15)

(16)

式中：M+為正趨勢序列的數(shù)目，M-為負趨勢序列的數(shù)目。當bxv(1)≠0時，二者之和即為子區(qū)間數(shù)：

(17)

(18)

步驟6若不同趨勢下位移與環(huán)境量的多重分形相關關系仍存在，則有式(5)所示的冪律關系。同理兩邊取對數(shù)：

(19)

(20)

|ΔHA|越大，表明不同趨勢下位移與環(huán)境量相關關系的多重分形非對稱程度越高，趨勢變化對多重分形強度的影響越顯著。當|ΔHA|為正時，說明上升趨勢下的位移與環(huán)境量相關性的持續(xù)性更強；當|ΔHA|為負時，說明下降趨勢下的位移與環(huán)境量相關性的持續(xù)性更強。

1.3 多尺度滑動時間窗

傳統(tǒng)多重分形分析易在區(qū)間劃分中產生偽波動和冗余數(shù)據(jù)，若將該部分忽略，會造成數(shù)據(jù)信息的丟失；若采用逆序處理法，則會打亂原數(shù)據(jù)的順序，影響信息獲取。鑒于此，結合大壩位移監(jiān)測數(shù)據(jù)的特殊性，采用滑動時間窗[17]對區(qū)間劃分方式進行優(yōu)化。

圖2 滑動時間窗Fig.2 Sliding time window

滑動時間窗通常指長度一定、起始時間和結束時間變化的窗口。給定窗口初始長度為s，兩端均以相同步長r推進，即一端加入r個數(shù)據(jù)，另一端依次分離出r個數(shù)據(jù)，運行原理如圖2所示?；瑒哟翱谠诒３謹?shù)據(jù)總量不變的同時，能夠實時更新窗口數(shù)據(jù)，運行一次得到子區(qū)間的數(shù)目為(N-s)/r+1，通過改變s的取值，使滑動窗的長度隨s的改變而改變，建立多尺度滑動時間窗模型。取r=1，并用改進后的式(21)～(23)代替式(3)(4)、式(15)(16):

(21)

(22)

(23)

2 實例應用

圖3 測點布置Fig.3 Layout of measuring points

某水電站主要由混凝土重力壩、壩后溢流式廠房、埋設于壩內的輸水系統(tǒng)、泄水底孔及過木筏等組成。大壩壩頂高程280 m，最大壩高78 m，壩頂全長253 m，水平位移引張線測點布設于壩頂處，測點布置見圖3。選取河床壩段7～10號對應的引張線測點EX5～EX8自動化監(jiān)測序列為樣本，樣本時間跨度自2008年1月1日至2014年12月31日，規(guī)定向下游位移為正、向上游位移為負，實測位移過程線和環(huán)境量如圖4所示。

圖4 實測位移過程線和環(huán)境量Fig.4 Hydrographs of measured displacement and environmental quantities

表1 主成分貢獻率

由圖4(a)可知，各測點監(jiān)測數(shù)據(jù)同升同降，同方向、同高程、不同位置各測點位移間存在關聯(lián)性。傳統(tǒng)方法基于權重思想[18-19]建立各測點的權重關系，但計算量大且權重確定方法有待驗證。而運用PCA[20-21]進行數(shù)據(jù)處理時，可保證在有效信息丟失最少的情況下，提取最能表征大壩多測點變形整體關系的主成分綜合效應量，該量是各測點效應量的線性組合，可將其視為單測點效應量建模分析，降低了分析的復雜度。表1為PCA的主成分貢獻率。

根據(jù)文獻[20-21]，PCA主成分的閾值為85%，由表1可知，第一主成分貢獻率1已達到PCA的閾值，能夠在反映有效信息的同時將數(shù)據(jù)最簡化，故選擇第一主成分表征多測點位移信息。根據(jù)SPSS分析給出的評分系數(shù)，帶入式(24)求解多測點位移主成分。

F=m1Y1+m2Y2+m3Y3+m4Y4

(24)

式中：F為主成分；Y1，Y2,Y3,Y4分別為位移；m1，m2,m3,m4分別為評分系數(shù)。

2.1 水平位移與環(huán)境量的多重分形相關性分析

將MF-DCCA應用到水平位移與環(huán)境量的多重分形相關性分析中，在matlab環(huán)境下運行MF-DCCA程序，分別計算該混凝土重力壩水平位移與環(huán)境量的波動函數(shù)Fq(s)、廣義相關Hurst指數(shù)H12(q)、標度函數(shù)τ(q)和多重分形譜f(α)，如圖5、圖6所示。

圖5 位移與水位的多重分形相關性分析Fig.5 Multifractal correlation analysis of displacement and water level

圖6 位移與溫度的多重分形相關性分析Fig.6 Multifractal correlation analysis of displacement and temperature

a.結合圖5(a)和圖6(a)，在q一定時，波動函數(shù)Fq(s)與s近似呈正相關關系。整體來看，q越小，F(xiàn)q(s)受波動影響程度越大，位移與環(huán)境量相關性的波動隨著s的增大而逐漸有序，并趨于均化。

b.結合圖5(b)和圖6(b)，位移與環(huán)境量的廣義相關Hurst指數(shù)H12(q) 并非恒定不變，均隨著q的增加而減小，說明位移的多重分形復雜性與水位、溫度的影響均有關。位移與環(huán)境量的H12(q)均大于0.5，說明位移與環(huán)境量間的相關性具有明顯的正持續(xù)性特征，總體表現(xiàn)出良好的長程相關性，即位移序列數(shù)值的上升或下降伴隨著影響因子的上升或下降。當q較大或較小時，位移與溫度間的H12(q)變化趨勢相較于與水位間的H12(q)變化趨勢更劇烈，說明位移與溫度間的多重分形相關性對波動的敏感性更強。但從定量的角度來看，位移與水位間H12(q)明顯大于位移與溫度間的H12(q)，說明位移與水位間的正持續(xù)性強于位移與溫度間的正持續(xù)性，水位變化對位移的影響更顯著。

c.進一步分析圖5(c)和圖6(c)，τ(q)呈非線性遞增趨勢，且圖5(c)的上凸程度和非線性程度均強于圖6(c)，說明位移與環(huán)境量間的相關性具有多重分形特征，且位移與水位的這一特征要強于位移與溫度的組合，進一步印證了前述分析。

d.結合圖5(d)和圖6(d)，位移與環(huán)境量的相關性多重分形譜圖像呈單峰凸分布，刻畫了不同時刻局部變化的多樣性。α主要集中于圖像的兩側，α分布的不均勻也印證了位移與環(huán)境量間的多重分形相關性。相關性多重分形譜基本對稱，整體協(xié)同性較好，發(fā)展狀態(tài)穩(wěn)定，其中位移與溫度的相關性多重分形譜呈現(xiàn)不明顯的右鉤狀，說明在該情況下小波動的影響略占優(yōu)勢。進一步計算位移與環(huán)境量間的分形特征統(tǒng)計量，即位移-水位：ΔH=1.255 2，Δα=1.596 4；位移-溫度：ΔH=0.476 8，Δα=0.713 8，位移與環(huán)境量間的ΔH和Δα均不為0。由于Δα越大，相關性的分形奇異性越強，ΔH越大，相關性的分形強度越強，故可知此時水平位移與水位間的多重分形相關性顯著強于位移與溫度間的多重分形相關性，水位變化對水平位移的發(fā)展影響更顯著。加之水位的漲落變化在周期性上不如溫度的變化穩(wěn)定，屬于高頻非線性影響量，故其與位移間的多重分形相關性往往更復雜，波動性更強。

2.2 水平位移與環(huán)境量的非對稱多重分形相關性分析

基于MF-A-DCCA進一步探究當位移序列處于正趨勢或負趨勢時，其與環(huán)境量間相關性的多重分形特征，分析成果見圖7、圖8，非對稱分形特征統(tǒng)計量見表2。

圖7 位移與水位的非對稱多重分形相關性分析Fig.7 Correlation analysis of asymmetric multifractal of displacement and water level

圖8 位移與溫度的非對稱多重分形相關性分析Fig.8 Correlation analysis of asymmetric multifractal of displacement and temperature

表2 非對稱分形特征統(tǒng)計量

由圖7、圖8和表2可知，隨著q的變化，當位移分別處于正趨勢和負趨勢時，其與環(huán)境量的廣義相關Hurst指數(shù)曲線互不重疊，說明此時位移與環(huán)境量間相關性的多重分形特征具有明顯的非對稱性；無論是位移原序列、正趨勢還是負趨勢，其與環(huán)境量間的H12(q)并非定值，而是隨著q的變化而變化，且均大于0.5，說明位移與環(huán)境量的相關性無論在位移序列上升還是下降時都具備多重分形特性，同時兼具良好的長程相關性。從表2的結果可知，位移與環(huán)境量的正趨勢分形強度均大于負趨勢分形強度，說明當位移序列處于上升趨勢時，其與環(huán)境量間的正持續(xù)性比處于下降趨勢時的正持續(xù)性強；進一步對比表3同一序列類型下位移-水位、位移-溫度的計算結果可知，無論原序列、正趨勢還是負趨勢，前者的分形奇異性和分形強度都明顯高于后者，說明即使在不同的變化趨勢下，水位變化對位移多重分形特征的影響仍占優(yōu)勢，與前文結論對應。

圖9 位移與環(huán)境量相關性的非對稱度量Fig.9 Asymmetric metric of displacement and environmental quantities

3 結論

a.運用多重分形理論刻畫了混凝土重力壩水平位移與環(huán)境量相關性的多重分形特征，并研究了其非對稱性，利用廣義相關Hurst指數(shù)、相關性多重分形譜等分析了多重分形相關性的強度，剖析了在環(huán)境量的影響下，大壩變形性態(tài)的動力學特征和演變規(guī)律。

b.引入多尺度滑動時間窗模型消除由子區(qū)間劃分而產生的偽波動，相較于文獻[9]中未優(yōu)化的情形，優(yōu)化后所繪圖像平滑性較好，趨勢性明顯，減小了偽波動的影響，提高了多重分形分析的信息利用率。

c.實例分析表明，該混凝土重力壩水平位移與環(huán)境量的相關性具有較強的多重分形特征，且這種特征在位移處于不同趨勢時具有非對稱性，位移變化從整體波動到局部趨勢均表現(xiàn)出良好的記憶性和長程相關性。水平位移的多重分形特征受環(huán)境量的影響，其中水位的頻繁漲落起主導作用，在環(huán)境量的影響下，該混凝土重力壩頂水平位移性態(tài)良好，分析結果符合實際工程情況，表明本文的方法在混凝土重力壩水平位移的分析中具有較好的可行性和適應性。