MF-DFA在大壩安全監(jiān)測序列分析和整體性態(tài)識別中的應(yīng)用

胡 江，蘇懷智，馬福恒，李子陽

(1.南京水利科學(xué)研究院水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室， 江蘇 南京 210029； 2.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098)

變形和應(yīng)力是水壓力、溫度、滲流、材料和地質(zhì)特性等多因素的綜合反映，其中變形觀測直觀可靠，國內(nèi)外普遍將其作為大壩監(jiān)測最主要效應(yīng)量[1-3]。大壩與基巖有機(jī)組合成為一個復(fù)雜的動力系統(tǒng)，多因素的共同作用使大壩變形觀測時間序列(位移和撓度傾斜等)具有顯著的時空非線性特征。依據(jù)實(shí)測變形時間序列，借助先進(jìn)的數(shù)學(xué)模型，分析其中蘊(yùn)含的復(fù)雜物理機(jī)理，是大壩安全監(jiān)控的重要手段。

一般采用統(tǒng)計(jì)模型、灰色系統(tǒng)理論和模糊數(shù)學(xué)模型等分析大壩原型觀測時間序列，但由于時間序列的復(fù)雜非線性動力特性，很難用參數(shù)形式描述其中的數(shù)字特征。蘇懷智等[4-5]應(yīng)用分形學(xué)分析法挖掘觀測時間序列的內(nèi)在規(guī)律性，采用Hurst指數(shù)H和分形維數(shù)D分析大壩結(jié)構(gòu)的工作性態(tài)和發(fā)現(xiàn)病害隱患，然而采用單一的Hurst指數(shù)從整體上識別整個時間序列的結(jié)構(gòu)特征過于籠統(tǒng)，缺乏對時間序列中的小波動的細(xì)節(jié)分析。

多重分形消除趨勢波動分析 (multifractal detrended fluctuation analysis，MF-DFA) 方法是在單分形的基礎(chǔ)上提出的，是定量描述復(fù)雜序列非線性演化過程及其形成自相似結(jié)構(gòu)特征的有效手段，最先被廣泛應(yīng)用于自然現(xiàn)象、生命科學(xué)和醫(yī)學(xué)等多個領(lǐng)域[6-7]。通過MF-DFA，不僅可以根據(jù)廣義Hurst指數(shù)檢測序列的長程相關(guān)性和確定標(biāo)度不變性，還可以根據(jù)多維分形譜區(qū)別不同時間序列演變規(guī)律的差異。近年來，MF-DFA方法的潛力被逐步挖掘，Telesca等[8]將MF-DFA方法應(yīng)用到大地電磁場的研究中，由不同場地的大地電磁數(shù)據(jù)序列的多重分形分析，辨識了環(huán)境因素對序列波動的影響，由此解釋序列的動態(tài)變化。本文運(yùn)用MF-DFA方法研究大壩監(jiān)測時間序列的長程相關(guān)性和多重分形特征，通過多測點(diǎn)的異同分析，探索大壩系統(tǒng)整體監(jiān)控模型，以此對大壩健康進(jìn)行診斷。

1 MF-DFA方法

消除趨勢波動分析(detrended fluctuation analysis，DFA)方法最初是由Peng等[9]在研究DNA序列機(jī)理時提出的一種分析冪律函數(shù)關(guān)系的方法，該方法可有效濾去各階趨勢成分，適合分析非平穩(wěn)時間序列的長程相關(guān)性。在DFA方法的基礎(chǔ)上，Kantelhardt等[10]提出了MF-DFA方法，MF-DFA方法提高了非平穩(wěn)序列波動奇異性的多重分形結(jié)構(gòu)度量的準(zhǔn)確性，很快被被運(yùn)用于時間序列蘊(yùn)含的多重分形性質(zhì)的研究。

假定長度為N的原位監(jiān)測序列為{x1,x2,…,xN}，MF-DFA方法的計(jì)算過程如下：

a. 計(jì)算序列的累積離差：

(1)

式中：xave為序列的均值。

b. 把序列Yi分割成長度為s的Ns=int(N/s)個互不重疊的等長區(qū)間(圖1)，由于長度N通常不是尺度s的整數(shù)倍，序列尾部可能剩余部分。為了不丟棄這一序列部分，從序列尾部重復(fù)這一分割過程。因此，共分割得到2Ns個區(qū)間。

圖1 消除趨勢波動分析中去趨勢程序示意圖

c. 采用最小二乘擬合法計(jì)算序列2Ns個區(qū)間中每個區(qū)間的局部趨勢。定義時間尺度s對應(yīng)的消除趨勢時間序列Ys,i，與原始時間序列和擬合序列不同，消除趨勢時間序列Ys,i表示為相對局部平均的波動:

Ys,i=Yi-yv,i

(2)

式中：yv,i是第v個區(qū)間的擬合多項(xiàng)式。圖1表示了s分別為100和200時采用二次多項(xiàng)式擬合(即DFA2)的情形。擬合過程還可采用線性、三次式或更高階多項(xiàng)式，相應(yīng)地，被記為DFA1、DFA3和高階DFA。確定2Ns個區(qū)間中每個區(qū)間v的均方差，當(dāng)v=1,2,…,Ns時,有

(3)

當(dāng)v=Ns+1,Ns+2,…,2Ns時

(4)

d. 計(jì)算所有區(qū)間的平均值得到整個序列的q階波動函數(shù):

(5)

式中：下標(biāo)q可以取0之外的任何的實(shí)值，當(dāng)q=0時，由洛必達(dá)法則有

(6)

e. 在smin≈6和smax≈N/4的范圍內(nèi)變換尺度s的值，確定消除趨勢波動函數(shù)Fq(s)和尺度s大小間的標(biāo)度關(guān)系，即

Fq(s)-sh(q)

(7)

一般地，指數(shù)h(q)隨q不同而變化，對于非平穩(wěn)序列，h(2)等于Hurst指數(shù)，因此，將h(q)定義為廣義Hurst指數(shù)。對于單分形序列，h(q)不隨q不同而變化。這是因?yàn)閷τ谒蟹指顓^(qū)間v，均方差F2(v,s)相等，且由式(4)表示的平均過程計(jì)算得到的波動函數(shù)也相等。只有當(dāng)大、小波動明顯不同時，h(q)才會隨q顯著變化。當(dāng)q為正值時，由大的均方差F2(v,s)(即與相應(yīng)的擬合值偏差大)的區(qū)間v決定均值Fq(s)；相反地，當(dāng)q為負(fù)值時，由小的均方差F2(v,s)的區(qū)間v決定均值Fq(s)。換句話說，q為正值，h(q)刻畫了大波動區(qū)間的分形行為；q為負(fù)值，h(q)則刻畫了小波動區(qū)間的分形行為[10]。

圖2 大壩垂線的布置方案(單位：m)

f. 對于多重分形序列，MF-DFA的最后結(jié)果是一系列的分形指數(shù)h(q)，是q的減函數(shù)。根據(jù)標(biāo)準(zhǔn)多重分形分析，由MF-DFA方法計(jì)算得到的h(q)與標(biāo)度指數(shù)τ(q)的解析式[10]為

τ(q)=qh(q)-1

(8)

式(8)表明，MF-DFA中由式(7)定義的廣義Hurst指數(shù)h(q)與傳統(tǒng)的多重分形標(biāo)度指數(shù)τ(q)直接相關(guān)。另一種描述多重分形序列的方法是奇異譜f(α)，通過Legendre轉(zhuǎn)換得到f(α)和τ(q)間的關(guān)系

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

(9)

式中：α為奇異強(qiáng)度，或稱為H?lder指數(shù)；而奇異譜f(α)則表示了由α描述的序列子集的維數(shù)。由式(8)可得

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

(10)

奇異譜f(α)的形狀和延伸長度包含了分析序列的分布特征的顯著信息。通常奇異譜f(α)有一個下凹的曲率，分形譜寬度Δα(Δα=αmax-αmin)增加，意味著分布的奇異性增大。

2 應(yīng)用實(shí)例

2.1 工程背景和數(shù)據(jù)來源

某工程屬 I 等樞紐工程，主要由碾壓混凝土重力壩、壩頂開敞式溢洪道、泄水底孔、左岸輸水建筑物等組成。碾壓混凝土重力壩最大壩高 113.0 m，壩頂長308.5 m，壩頂高程為179.0 m，水庫正常蓄水位173.0 m，屬不完全調(diào)節(jié)水庫。主體工程于2000 年12 月18 日下閘蓄水。

為確保大壩及地下廠房的安全運(yùn)行，在主要建筑物表面(或內(nèi)部)布置了變形、滲流、溫度和應(yīng)力應(yīng)變等較為全面的監(jiān)測項(xiàng)目。大壩變形監(jiān)測項(xiàng)目主要有正垂線(PL)、倒垂線(IP)、引張線和視準(zhǔn)線，垂直位移監(jiān)測采用水準(zhǔn)測量。大壩正、倒垂線用于監(jiān)測壩頂及壩體內(nèi)部的水平順河向、垂直水流向位移量，布置方案如圖2所示。為提高可比性，本文選取懸掛點(diǎn)均位于壩頂179.0 m高程的4號壩段(溢流壩段)、2號壩段(擋水壩段)和5號壩段(擋水壩段)的正垂線PL7、PL5和PL6的自動化監(jiān)測時間序列進(jìn)行多重分形分析。為方便分析，順河向水平位移以向下游為正，向上游為負(fù)?？捎玫臅r間序列自2002年10月開始，選取2003年1月1日至2008年12月31日間的3個測點(diǎn)的順河向位移序列進(jìn)行分析，實(shí)測的環(huán)境量和位移時間序列如圖3所示。

圖3 環(huán)境量和正垂線順河向位移實(shí)測時間序列

2.2 統(tǒng)計(jì)描述與正態(tài)檢驗(yàn)

從圖3可以看出，變形呈一定程度的增加趨勢，但并非單調(diào)增加，而是有漲有落，且變形波動過程呈一定的周期性。計(jì)算PL5、 PL6和PL7測點(diǎn)徑向位移序列的標(biāo)準(zhǔn)差、峰度和偏度等基本統(tǒng)計(jì)量，結(jié)果見表1，可看出，各序列的偏度都不為 0，且峰度都小于0，因此均屬于非正態(tài)分布。PL7和PL5測點(diǎn)位移序列分布的偏度均大于0，這表示其數(shù)據(jù)分布形態(tài)與正態(tài)分布相比為右偏；PL6測點(diǎn)位移序列分布的偏度小于0，其數(shù)據(jù)分布形態(tài)與正態(tài)分布相比為左偏。同時，所有測點(diǎn)序列的峰度均小于0，這表示測點(diǎn)序列的整體分布與正態(tài)分布相比較為平緩，為寬頂峰。由此可知，測點(diǎn)序列為呈右偏或左偏的非正態(tài)分布。

表1 正垂線徑向位移時間序列的基本統(tǒng)計(jì)量

2.3 多重分形分析結(jié)果

根據(jù)式(1)～(4)計(jì)算波動函數(shù)Fq(s)，參照文獻(xiàn)[10]，這里波動函數(shù)的階數(shù)取值范圍為-6≤q≤6；對于時間尺度s，取值需適中，取值太小沒有實(shí)際意義，過大，則分割的區(qū)間數(shù)太少，結(jié)果將與單分形結(jié)果相似，一般地，時間尺度s的范圍為5≤s≤N/4，根據(jù)大壩變形監(jiān)測資料的周期性波動趨勢，這里選取s={6,12,24,48,96,192,384}。圖4為采用MF-DFA1計(jì)算得到的q=±6、±4、±2、0時溫度和庫水位及PL5、PL6和PL7測點(diǎn)徑向位移的Fq(s)-s對數(shù)圖。整體來說，小的時間尺度對應(yīng)的q階波動函數(shù)的差異較大的時間尺度顯著，小尺度能夠識別局部時間內(nèi)的大小波動；相反，大尺度在一定程度上均化了大小波動的幅度。

圖4 波動函數(shù)和時間尺度的對數(shù)關(guān)系

圖5為廣義Hurst指數(shù)h(q)隨q的變化情況，h(q)-q曲線能反映序列是否存在多重分形行為。當(dāng)h(q)不隨q變化時，表明序列是單分形；反之，序列為多重分形。由圖5可以看出，由環(huán)境量序列和各測點(diǎn)位移序列計(jì)算得到的廣義指數(shù)h(q)均隨q增加而單調(diào)減少，這意味著環(huán)境量序列和各測點(diǎn)位移序列都存在較顯著的多重分形特征。當(dāng)q=2時，計(jì)算得到的溫度和庫水位序列h(2)的值分別為1.31和1.46，各測點(diǎn)位移序列h(2)的值分別為1.38、1.43和1.5。

圖5說明位移序列存在明顯的分形特征，因此有必要進(jìn)行多重分形譜分析。由式(8)和式(9)分別計(jì)算得到多重分形標(biāo)度指數(shù)τ(q)和分形譜寬度Δα、奇異譜高度Δf(α)，得到的τ(q)-q和f(α)-α曲線分別如圖6和圖7所示。其中，溫度和庫水位序列的分形譜寬度Δα分別為0.612和1.227，奇異譜高度Δf(α)分別為0.715和1.122；PL5、PL6和PL7測點(diǎn)序列的分形譜寬度Δα分別為0.614、0.728和0.679，奇異譜高度Δf(α)分別為0.72、0.79和0.69。

圖5 廣義Hurst指數(shù)和q的關(guān)系

圖6 位移日測值序列質(zhì)量指數(shù)譜曲線

圖7 位移日測值序列多重分形奇異譜曲線

2.4 結(jié)果討論和大壩的整體性態(tài)識別

從圖4(a)(b)可以看出，溫度序列的Fq(s)-s曲線在s=96處扭曲，對應(yīng)時間尺度大約為一個季度；水位序列在s=192處扭曲，對應(yīng)時間尺度大約為半年周期。從圖4(c)(d)(e)可以看出，PL5測點(diǎn)序列的Fq(s)-s曲線都在s=96處扭曲；PL6測點(diǎn)序列有所不同，在s=24處接近月周期附近扭曲；PL7測點(diǎn)曲線扭曲發(fā)生在s=384接近年周期處。同時，由圖3(b)可以看出，各測點(diǎn)位移序列大小波動存在一定的差異，F(xiàn)q(s)-s曲線較好地反映了測點(diǎn)序列局部時間內(nèi)和較大尺度范圍內(nèi)的日測值的波動情況。

對測點(diǎn)序列進(jìn)行統(tǒng)計(jì)分析，其中水壓因子取3次方，溫度取半年周期，時效取對數(shù)形式[1]，統(tǒng)計(jì)模型分析得到2007年典型測點(diǎn)年變幅的結(jié)果(表2)，計(jì)算得到PL5、PL6和PL7測點(diǎn)序列的溫度分量約占總位移年變幅的61.7%、63.5%和63.0%，溫度變化對壩體上下游方向水平位移變化影響最大，其次為庫水位變化；時效分量相對較小，表明壩體變形趨于穩(wěn)定。

表2 統(tǒng)計(jì)模型分析結(jié)果 mm

測點(diǎn)序列的h(2)值大于1且接近于1.5，這說明測點(diǎn)序列的波動相關(guān)性和布朗運(yùn)動有一定相似性，換句話說，監(jiān)測序列存在較大的不確定性。同時，各序列h(2)值接近，由文獻(xiàn)[4]可知，整體上大壩結(jié)構(gòu)處于彈性狀態(tài)。由h(q)-q曲線q>0和q<0兩段可以區(qū)分各測點(diǎn)的大小波動的差異。圖5(a)表明溫度序列的大波動標(biāo)度無顯著變化，庫水位序列的大波動標(biāo)度相對較大；溫度和庫水位序列尤其是庫水位序列的小波動標(biāo)度單調(diào)遞減較顯著。圖5(b)表明PL7測點(diǎn)序列波動較顯著，其次為PL6測點(diǎn)系列，PL5測點(diǎn)序列波動較?。桓鳒y點(diǎn)序列的大波動標(biāo)度單調(diào)遞減較明顯，這是由于蓄水和運(yùn)行初期，環(huán)境量和壩體自身引起的，各測點(diǎn)遞減的趨勢基本相同，表明整體變化趨勢相同；PL6和PL7測點(diǎn)序列的小波動標(biāo)度單調(diào)遞減趨勢也基本相同；PL5測點(diǎn)序列h(q)序列小波動的標(biāo)度較其他兩測點(diǎn)小，表明PL5測點(diǎn)序列的波動較小。

f(α)-α曲線可以反映位移序列在分形結(jié)構(gòu)上不均勻分布的性質(zhì)。由圖7(a)可見溫度和庫水位序列尤其是庫水位序列的奇異譜呈明顯的右鉤狀，小波動的影響在庫水位序列中占顯著優(yōu)勢。由圖7(b)可見各測點(diǎn)序列尤其是PL5和PL7的奇異譜基本對稱；而PL6測點(diǎn)序列奇異譜則呈不明顯的右鉤狀，即對于PL6測點(diǎn)序列，小波動的影響在日測值序列中稍占優(yōu)勢。如前所述，分形譜寬度Δα的值可以反映大壩性態(tài)演變復(fù)雜度的空間差異，各測點(diǎn)分形譜寬度Δα相差不大，相對而言，PL6測點(diǎn)序列的奇異譜譜寬Δα最大，PL5次之，PL7最小。從圖2可以看出，PL6測點(diǎn)高程較低，水位波動相對壩段壩高在3個測點(diǎn)中最大；PL7測點(diǎn)位于河床，水位波動相對壩段壩高在3個測點(diǎn)中最小。初步推斷，各壩段不同的運(yùn)行機(jī)理是導(dǎo)致壩頂變形測點(diǎn)序列多重分形特征差異的根本原因。

由以上分析和討論可以看出，整體上，大壩結(jié)構(gòu)變形處于正常狀態(tài)，壩體處于彈性狀態(tài)；對比來看，岸坡5號壩段的PL6測點(diǎn)序列表現(xiàn)出較為復(fù)雜的性態(tài)，河床4號壩段的PL7測點(diǎn)序列較為穩(wěn)定。綜上，統(tǒng)計(jì)模型方法能夠得到監(jiān)測效應(yīng)量的主要影響因子；MF-DFA方法可以識別監(jiān)測效應(yīng)量序列不同時間尺度上的分形特征和大壩整體性態(tài)時空差異和演變規(guī)律。

3 結(jié) 論

a. 各垂線測值序列分別在月、季度和年周期時間尺度上，波動函數(shù)存在一個轉(zhuǎn)折點(diǎn)，反映了各測點(diǎn)測值時間序列在不同時間尺度上的波動特征，說明各壩段運(yùn)行機(jī)理具有一定的差異。環(huán)境量的波動函數(shù)表明溫度和水位存在較為顯著的季度和年周期波動，同時廣義Hurst指數(shù)說明水位大小波動較為顯著；通過統(tǒng)計(jì)模型分析表明，溫度分量和水位分量占變形量年變幅的大部分。綜合統(tǒng)計(jì)模型分析和MF-DFA方法分析結(jié)果，整體上，溫度和庫水位的大波動決定了壩體變形幅度，水位頻繁漲落控制各壩段尤其是岸坡壩段變形的小波動。

b. 廣義Hurst指數(shù)h(q)均隨q的增加而減小，大壩變形測值序列呈顯著的多重分形特征。由Hurst指數(shù)大于1且接近于1.5說明運(yùn)行初期大壩整體上處于正常線彈性狀態(tài)。測點(diǎn)序列波動具有較大的不確定性，這與大壩運(yùn)行初期受到較為復(fù)雜的環(huán)境、水文和運(yùn)行管理等多重因素影響有關(guān)。

c. 綜合環(huán)境量和測點(diǎn)序列的奇異譜形狀和寬度，可以刻畫大壩不同壩段和局部運(yùn)行性態(tài)的動力學(xué)特征，其物理意義明確，為此通過大壩多個測點(diǎn)的奇異譜寬度分析可反映大壩長期性態(tài)演變的時空差異。進(jìn)一步挖掘造成不同測點(diǎn)復(fù)雜程度差異的可能原因，如壩段的工作條件和材料的衰變等，可以指導(dǎo)大壩的運(yùn)行管理和維護(hù)。

參考文獻(xiàn)：

[1] 顧沖時，吳中如.大壩與壩基安全監(jiān)控理論和方法及其應(yīng)用[M].南京:河海大學(xué)出版社，2006.

[2] 遠(yuǎn)近，鄭東健，李波，等.模糊可拓模型在壩體變形評價中的應(yīng)用[J].水利水電科技進(jìn)展，2009，29(1)：23-25.(YUAN Jin,ZHENG Dongjian,LI Bo,et al.The application of fuzzy extension model in dam deformation evaluation[J].Advances in Science and Technology of Water Resources,2009,29(1):23-25.(in Chinese)

[3] 胡江，蘇懷智.再生縫對拱壩變形的影響分析與預(yù)測[J].水電能源科學(xué)，2008，26(5)：96-100.(HU Jiang,SU Huaizhi.Influence of regenerated crack to arch dam deformation and displacement forecast[J].Water Resources and Power,2008,26(5):96-100.(in Chinese)

[4] SU Huaizhi,HU Jiang,WU Zhongru.A study of safety evaluation and early-warning method for dam global behavior[J].Structural health Monitoring:An International Journal,2012,11(3):269-279.

[5] 賴道平，吳中如，周紅.分形學(xué)在大壩安全監(jiān)測資料分析中的應(yīng)用[J].水利學(xué)報(bào)，2004(1)：100-104.(LAI Daoping,WU Zhongru,ZHOU Hong.Application of fractal theory to analyze dam safety monitoring data[J].Journal of Hydraulic Engineering,2004(1):100-104.(in Chinese)

[6] TELESCA L,PIERINI J O,SCIAN B.Investigating the temporal variation of the scaling behavior in rainfall data measured in central Argentina by means of the detrended fluctuation analysis[J].Physica A,2012,391(4):1553-1562.

[7] 陳瑩，許有鵬，尹義星，等.長江干流日徑流序列的多重分形特征[J].地理研究，2008，27(4)：819-828.(CHEN Ying,XU Youpeng,YIN Yixing,et al.Multifractal characteristics of daily discharge series in the Yangtze River[J].Geographical Research,2008,27(4):819-828.(in Chinese)

[8] TELESCA L,LOVALLO M.,HSU H.L,et al.Analysis of site effects in magnetotelluric data by using the multifractaldetrended fluctuation analysis[J].Journal of Asian Earth Sciences,2012,54-55:72-77.

[9] PENG C K,BULDYREV S V,HAVLIN S S,et al.Mosaic organization of DNA nucleotides[J].Physical Review E,1994,49(2):1685-1689.

[10] KANTELHARDT J W,ZSCHIEGNER S A,KOSCIELNY-BUNDE E,et al.Multifractal detrended fluctuation analysis of nonstationary time series[J].Physica A,2002,316(1):87-114.