周蘭庭,柳志坤

(河海大學水利水電學院，江蘇 南京 210098)

大壩變形監(jiān)測數據是大壩安全運行的直觀體現，其主要受溫度、水壓、時效三大因素影響[1-2]。由于大壩整體處于一個復雜的動態(tài)時空中，由此得到的位移序列也具有較強的非線性和不確定性，增加了資料分析的困難度，本文將引入多重分形理論[3-5]來細致刻畫位移序列的內在規(guī)律，從多重分形的視角描述大壩位移的變化性態(tài)。

分形是自然界的一種常見現象，具有整體與局部的自相似性，可以分為數學分形與統(tǒng)計分形兩類，數學分形多由迭代產生，是理想化的嚴格自相似，并且具有無窮的層次結構，而統(tǒng)計分形則具有弱化的且層次結構有限的自相似性，相比于數學分形而言，其實際研究意義更突出，如海岸線輪廓、污染物濃度等。在大壩變形分形特征領域，賴道平等[6]應用重標度極差分析法(rescaled range analysis，R/S)對大壩變形資料進行了單分形分析，用Hurst指數、分形維數刻畫了大壩的運行狀況；YIN等[7]在單分形基礎上，通過對比雙對數散點圖提出確定時序統(tǒng)計循環(huán)長度的方法，后在大壩監(jiān)測領域得到推廣；田振華等[8]利用R/S方法分析監(jiān)測資料，通過分形維數將監(jiān)測數據與大壩損傷聯系起來，從分形維數的角度定量刻畫大壩的結構損傷程度；趙卿等[9]利用去趨勢波動分析法(detrended fluctuation analysis，DFA)對大壩多測點進行分析，從單分形角度刻畫了位移序列的變化特征，取得了較好的效果。上述單分形理論往往只在宏觀上刻畫位移時序的分形特性，但實際大壩位移由于受多種因素的影響導致其實測數據往往具有復雜的多重分形特性，僅靠單分形理論無法完整刻畫其內部規(guī)律，于是在單分形理論的基礎上進一步引入了多重分形理論，該方法在大壩變形性態(tài)分析方面得到一定應用。胡江等[10]利用多重分形去趨勢波動分析法(multifractal detrended fluctuation analysis，MF-DFA)對大壩位移序列和環(huán)境量進行了分析，對大壩變形多重分形特征和整體性態(tài)進行了辨識；韋武昌等[11]采用滑動平均優(yōu)化傳統(tǒng)的MF-DFA方法，并在某雙曲拱壩的實測資料上進行了應用，對該壩的變形性態(tài)進行了分析；楊景文等[12]基于多重分形理論對大壩多個裂縫序列進行了多重分形分析，并探討了裂縫序列多重分形特征的成因，均取得較好的效果。

目前已有研究主要基于實測數據進行單分形和多重分形分析，從而刻畫大壩變形的演化規(guī)律和發(fā)展趨勢，而在位移序列的多元多重分形分析、多重分形非對稱性以及算法優(yōu)化等方面探索較少。故本文以大壩實測位移數據為基礎，引入滑動時間窗口對傳統(tǒng)多重分形子區(qū)間劃分方法進行優(yōu)化，基于MF-DFA和多元多重分形去趨勢波動分析法(multivariate multifractal detrended fluctuation analysis，MV-MFDFA)進行單測點和多元測點的聯合多重分形分析，并用非對稱多重分形去趨勢波動分析法(asymmetric multifractal detrended fluctuation analysis，A-MFDFA)進一步分析位移測點不同變化趨勢下的分形特征，實現大壩位移從整體波動到局部趨勢的多重分形描述和解讀，為多重分形與大壩位移性態(tài)分析的結合提供更詳盡的信息。

1 多重分形計算方法

1.1 多重分形去趨勢波動分析法

DFA[13]源于核糖核酸行為機理的研究，其優(yōu)勢在于能夠消除非平穩(wěn)時序的趨勢部分，但其對于非平穩(wěn)時序的分形特征刻畫不夠完整，后人在DFA的基礎上引入多重分形的概念，將DFA改進為MF-DFA[14-15]，經典的MF-DFA計算步驟如下。

(1)

步驟 2給定一時間尺度s，將序列y(t)以s為標準進行等分，總共可劃分為m個等長連續(xù)且不重疊的子區(qū)間，m=int(N/s)。從實際數據來看，N未必能夠整除s，即可能出現尾部數據冗余，故此處采用逆序處理法來充分利用數據信息，即在正序劃分的基礎上，從序列的末端開始重復前述操作，由此得到2m個子區(qū)間。

步驟 3對各個子區(qū)間進行趨勢擬合，并在原基礎上減去趨勢部分，得到相應殘差序列，記為zv(t)。

(2)

步驟 4計算殘差序列zv(t)的均方差F2(s,v)

(3)

步驟 5對F2(s,v)數據集取均值，按式(4)計算序列的q階波動函數Fq(s)

(4)

q的取值與Fq(s)受波動影響的程度有關，原則上可取任意非零實數。特別地，當q=0時，有式(5)的極限形式：

(5)

傳統(tǒng)的MF-DFA易在區(qū)間的劃分中產生偽波動，對后續(xù)分析造成干擾，此外，時序長度未必能夠整除子區(qū)間長度，造成數據的冗余，若將冗余部分忽略，會造成數據信息的丟失，若采用逆序處理法，則會打亂原來數據的順序，影響信息的獲取。鑒于此，考慮大壩變形數據的特殊性，此處采用滑動窗口[16-17]對傳統(tǒng)MF-DFA劃分子區(qū)間的方式進行優(yōu)化，通過一定長度的窗口在序列上按照一定的步長進行滑動取值，以減少數據的偽波動并充分利用數據信息。設窗口長度為s，序列長度為N，取滑動步長為1，則運行一次得到的子區(qū)間數目為N-s+1，并用式(6)代替式(4)(5)。

(6)

通過上述步驟可以得到某一尺度s對應的q階波動函數Fq(s)，變換s的取值，重復上述步驟可得一系列的s-Fq(s)點值，若該時間序列存在長程相關性，則Fq(s)與s間有式(7)所示的冪律關系：

Fq(s)∝sh(q)

(7)

對上式兩邊取對數，得到式(8)的形式：

lgFq(s)=h(q)lgs+lgb

(8)

式中：Fq(s)為序列的q階波動函數；h(q)為相應的廣義Hurst指數；b為常系數。作lgFq(s)-lgs雙對數散點圖，并進行擬合，其斜率即為廣義Hurst指數h(q)，若h(q)為一固定常數，表明序列為單分形，無多重分形特征；若h(q)為q的非線性減函數，則位移序列具有多重分形特征，當h(q)<0.5時，位移序列表現為具有反持續(xù)性的記憶過程，當h(q)=0.5時，位移序列表現為不相關的隨機過程，當h(q)>0.5時，位移序列表現為具有正持續(xù)性的記憶過程，另外，當h(q)>1時，位移序列表現為具有強非平穩(wěn)性的長程正相關過程。

時間序列的分形強度和分形奇異性通常可由多重分形譜f(α)來表征，按式(9)～(11)求解：

τ(q)=qh(q)-1

(9)

α=τ′(q)

(10)

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

(11)

式中：τ(q)為Renyi指數，也叫標度函數，若其為q的非線性上凸函數，則位移序列具有多重分形特征，若其為q的線性函數，則位移序列表現為單分形特征，故該特征也常作為判斷序列是否具有多重分形的判據。α為奇異強度，f(α)為多重分形譜，當α-f(α)呈單峰上凸狀，形似二次函數時，說明位移序列具有多重分形特征；當α-f(α)聚集到一個點附近時，說明序列時單分形的。此外，常用分形譜寬度Δα和分形強度Δh來定量描述位移序列的多重分形特性，Δα主要反映大壩變形的奇異性和性態(tài)演變的空間差異性，Δh主要反映序列的分形強度，并與參數值成正比。

Δh=max[h(q)]-min[h(q)]

(12)

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

(13)

1.2 多元多重分形去趨勢波動分析法

1.1節(jié)所述MF-DFA主要是針對單個測點的一元時間序列的多重分形分析，對于多個測點組成的多元時間序列集，此處引入MV-MFDFA[18-19]進行分析，設一個多元時序集{xt,i}，其中t取值范圍為[1,N]，i表示測點個數，計算步驟如下。

步驟 4除去趨勢部分，得到測點多元時間序列的局部方差F2(s,v)：

(14)

式中最大范數按式(15)計算：

‖X-Y‖=‖(x1,x2,...,xn)-(y1,y2,...,yn)‖=

max{|x1-y1|,...,|xn-yn|}

(15)

步驟 5按1.1節(jié)中式(6)求解q階波動函數Fq(s)。

采用1.1節(jié)中的方法，可以得到多元序列整體的廣義Hurst指數h(q)，當h(q)>0.5時，系統(tǒng)呈現出正相關性，即系統(tǒng)今后的增減與過去的增減正相關；當h(q)<0.5時，系統(tǒng)呈現出負相關性，即系統(tǒng)今后的增減與過去的增減負相關；當h(q)=0.5時，系統(tǒng)具有較強的隨機性，相關性不明確。多元序列整體的多重分形譜計算方法同1.1節(jié)所述。

1.3 非對稱多重分形去趨勢波動分析法

傳統(tǒng)的多重分形理論只能辨識序列的多重分形特征，但無法捕捉序列在不同趨勢上的分形特征，此處引入A-MFDFA[20-21]進一步加強時序的非對稱多重分形特性研究。該方法前五步同MF-DFA，擬從步驟6開始介紹。

步驟 6通過線性擬合估計原序列x(t)各個子區(qū)間的斜率bv，當bv>0時，記為正趨勢，當bv<0時，記為負趨勢，并按式(16)～(17)計算正、負趨勢序列的q階波動函數。

其中

分別表示正、負趨勢序列的數目，當bv≠0時，二者之和為子區(qū)間數目。

若各類序列具有長程相關性，則有形如上文式(7)的冪律關系，同理兩邊取對數，可得式(18)的形式：

(18)

通過擬合得到相應的h(q)、h+(q)、h-(q)，若正、負趨勢序列的廣義Hurst指數圖重合，即h+(q)=h-(q)，說明位移序列在多重分形特征上具有對稱性，反之，則具有非對稱性，即在位移序列處于不同趨勢時，其表現出的多重分形特征也不盡相同。非對稱多重分形分析正、負趨勢序列的多重分形譜計算方法同1.1節(jié)中所述。

2 工程實例

某水電站樞紐布置如圖1所示，選取圖中河床壩段⑦號、⑧號、⑨號對應的EX5、EX6、EX7測點自2008年1月1日至2014年12月31日的壩頂水平位移監(jiān)測數據，對各測點部分缺失數據進行3次樣條插值，保證數據每天1個測值，去噪后的各測點實測過程線如圖2所示，并規(guī)定向下游位移為正、向上游位移為負。

圖1 某水電站樞紐布置圖

圖2 位移實測過程線

2.1 數據特征描述

在研究大壩位移的多重分形特征之前，應先了解樣本位移數據的基本統(tǒng)計特征，此處對EX5～EX7測點數據進行描述性統(tǒng)計如表1所示。

表1 EX5～EX8測點位移描述性統(tǒng)計

根據K-S準則，考慮樣本均值、標準差、偏度和峰度四個評價指標，各測點的位移均值均小于0，表現為向上游位移；從標準差來看，EX7標準差最大，EX5和EX6標準差接近，說明EX7位移序列偏離均值的離散程度比其他各測點大，而其他2測點的偏離程度和波動程度均較小，這一點也可從圖2中直觀印證；各測點的偏度均大于0，說明數據分布均為右偏態(tài)，即均值右側的離散度比左側強，表現為一定程度的右側長尾，同時各測點峰度均顯著大于0，說明數據高峰比正態(tài)分布的高峰更陡峭，處于尖峰狀態(tài)，綜合J-B統(tǒng)計量的結果，在峰度和偏度的共同作用下，各測點的位移序列均不服從正態(tài)分布；此外，對各測點進行長記憶性分析，所得長記憶性指標均大于0.5，說明各測點位移序列均具有較好的長程相關性，并非隨機游走序列。綜上所述，各測點位移序列均不符合正態(tài)分布特征，并呈現出Levy分形分布所特有的尖峰胖尾特點，表明位移序列具有明顯的分形特征，可用分形理論的相關方法進行研究。

2.2 多重分形分析

采用滑動時間窗口優(yōu)化MFDFA和MV-MFDFA，分別對測點和由測點構成的多元測點序列進行多重分形分析，此處波動階數q取值范圍為[-10,10]，尺度s的取值范圍為[10,100]，滑動窗口步長為1，分析成果如圖3～5所示。

圖3 廣義Hurst指數

圖3和圖4分別為測點整體及單個測點的廣義Hurst指數和標度函數τ(q)變化情況，由圖3可知，當q在[-10,10]之間變化時，無論是測點整體還是各測點分量序列，其廣義Hurst指數均非常數，而是隨著q的變化呈非線性遞減趨勢，說明從測點整體到內部各測點位移序列均具有明顯的多重分形特性，僅通過單分形理論進行描述是不充分的；在不同的波動階數q下，整體分析的廣義Hurst指數曲線集中在各單測點的下部波動，表現出較弱的多重分形性，但無論整體或分量的h(q)值仍明顯大于0.5，說明變形監(jiān)測序列從整體到局部分量均具有較好的記憶性和長程相關性，兼具非平穩(wěn)性和隨機性；當q<0時，h(q)主要受小波動方差的影響，當q>0時，h(q)主要受大波動方差的影響，由圖3可知，當q<0時，各測點h(q)值分布較離散，正持續(xù)性強弱差異明顯，而當q>0時，各測點h(q)值一致性較好，具有相近的正持續(xù)性，此外，當q>0時測點間h(q)值的差異比q<0時小，說明在水位等大波動因素的影響下，各測點遞減的趨勢基本相同，整體變化趨勢也相同，測點所在區(qū)域保持較穩(wěn)定的長程相關性，壩體狀態(tài)良好。在h(q)的基礎上，通過式(9)計算Renyi指數即標度函數τ(q)，如圖4所示，各測點標度函數一致性較好，中部呈上凸狀，滿足τ(0)=-1，且整體為非線性關系，這一點進一步印證了各位移測點具有多重分形特征。

圖4 標度函數

圖5為整體分析及各測點分量的多重分形譜，各多重分形譜圖像呈單峰凸分布，形似二次函數曲線，位移多重分形的局部標度不是常數，刻畫了不同時刻局部變化的多樣性，奇異強度α主要集中于圖像的兩側，反映了位移序列的分形結構分布不均勻，α分布的不均勻也印證了測點序列的多重分形特性；無論整體或測點分量的多重分形譜基本對稱，整體協同性較好，發(fā)展狀態(tài)穩(wěn)定，其中EX6的分形譜呈現不明顯的右鉤狀，說明對該點序列而言，小波動的影響略占優(yōu)勢；f(α)表征位移序列具有相同奇異指數α的子區(qū)間的分形維，與位移序列的分布特征和分形強度相關。由式(12)(13)計算多元序列及其分量的分形特征統(tǒng)計量如表2所示，從表2可以看到多元序列及其分量的分形強度Δh和多重分形譜寬度Δα均為最小，這是由于各測點位于同一高程的相鄰壩段，測點變形之間并非相互獨立的，而是存在一定的關聯性，這種關聯性削弱了多元測點整體的多重分形性，使得整體分析的分形特征統(tǒng)計量不會表現的比分量更強；整體來看，三個測點的Δα值接近，變形具有相似的奇異性和空間差異性，各測點水平位移的概率分布范圍相對接近，但分形強度Δh值略有差異，表現為EX5、EX6的分形強度均強于EX7，結合實際情況初步推斷，三個測點均位于河床中部同一高程的相鄰壩段，變形均受水位波動等因素的影響，從而表現出一定的奇異性和多重分形特征，但由于EX5、EX6所在壩段的運行機理與EX7所在壩段不同，使得測點間的多重分形特征存在一定的差異。

圖5 多重分形譜

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

2.3 非對稱多重分形分析

在上述研究的基礎上，進一步引入多重分形的非對稱分析法，即A-MFDFA，并同樣采用滑動窗口進行優(yōu)化，該方法考慮到位移序列在上升或下降趨勢時的分形特征，故能夠更全面刻畫位移的分形特性和非對稱性，正、負趨勢序列的計算方法見前文所述，運行程序計算正、負趨勢序列及原序列的廣義Hurst指數、Renyi函數和多重分形譜，如圖6所示。

圖6 非對稱多重分形分析

圖6(a)(c)(e)為各測點不同趨勢下的廣義Hurst指數圖，圖6(b)(d)(f)為各測點不同趨勢下的多重分形譜。從圖6中可以看出，各測點位移原序列、正趨勢、負趨勢下的廣義Hurst指數和多重分形譜均存在不同程度的偏差，圖像并非重合的，說明位移序列在多重分形特征上具有明顯的非對稱性；位移序列的整體波動是正、負趨勢波動綜合作用的結果，正趨勢和負趨勢的h(q)均隨q的增加而減小，同樣表現出明顯的多重分形特征，且h(q)均大于0.5，說明位移序列無論是處于正趨勢還是負趨勢時其發(fā)展變化仍具有良好的記憶性和長程相關性；綜合各測點正趨勢和負趨勢的變化情況來看，q>0時h(q)的離散程度普遍高于q<0時的離散程度，說明水位等大波動因素對多重分形非對稱性的影響最顯著。

此處仍用多重分形普寬度Δα和分形強度Δh兩個參數來定量評價位移序列的非對稱多重分形特性，計算結果如表3所示。由表3可知，位移序列在各趨勢上的分形特征均不相同，有明顯的非對稱特點，其中，EX5和EX6各自正、負趨勢的分形特征接近，正負趨勢各有強弱但分形強度差異較小，二者位移發(fā)展走勢相近，相比之下由于所在壩段運行機理的差異，EX7測點正趨勢的分形強度要明顯大于負趨勢的分形強度，原序列分形強度與正趨勢分形強度接近，說明EX7位移序列在上升時的分形強度大于下降時的分形強度，即在上升時傾向于波動上升，而在下降時傾向于穩(wěn)步下降。

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

3 結 論

本文將多重分形理論應用到某混凝土重力壩的變形性態(tài)分析中，并利用滑動時間窗改進傳統(tǒng)多重分形子區(qū)間的劃分方式，聯合MF-DFA和MV-MFDFA對單測點和多測點序列進行多重分形特征辨識，并利用A-MFDFA進一步刻畫序列局部處于上升或下降狀態(tài)時的多重分形特征，分析結果表明，某大壩變形實測數據具有明顯的分形序列特征，在水位等大波動因素的影響下，單測點和多測點均保持穩(wěn)定的多重分形特征和長程相關性，壩體變形整體性態(tài)較好，處于穩(wěn)定的安全狀態(tài)。研究成果從變形的整體波動和局部趨勢兩個層面解析了大壩變形的多重分形特征及其非對稱性，從一個新的視角看待大壩變形的波動特征，實現對大壩變形性態(tài)的準確把握，有利于大壩運行的管理和維護。