胡 江,蘇懷智,馬福恒,李子陽
(1.南京水利科學(xué)研究院水文水資源與水利工程科學(xué)國家重點實驗室, 江蘇 南京 210029; 2.河海大學(xué)水利水電學(xué)院,江蘇 南京 210098)
變形和應(yīng)力是水壓力、溫度、滲流、材料和地質(zhì)特性等多因素的綜合反映,其中變形觀測直觀可靠,國內(nèi)外普遍將其作為大壩監(jiān)測最主要效應(yīng)量[1-3]。大壩與基巖有機(jī)組合成為一個復(fù)雜的動力系統(tǒng),多因素的共同作用使大壩變形觀測時間序列(位移和撓度傾斜等)具有顯著的時空非線性特征。依據(jù)實測變形時間序列,借助先進(jìn)的數(shù)學(xué)模型,分析其中蘊(yùn)含的復(fù)雜物理機(jī)理,是大壩安全監(jiān)控的重要手段。
一般采用統(tǒng)計模型、灰色系統(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)性和多重分形特征,通過多測點的異同分析,探索大壩系統(tǒng)整體監(jiān)控模型,以此對大壩健康進(jìn)行診斷。
消除趨勢波動分析(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方法的計算過程如下:
a. 計算序列的累積離差:
(1)
式中:xave為序列的均值。
b. 把序列Yi分割成長度為s的Ns=int(N/s)個互不重疊的等長區(qū)間(圖1),由于長度N通常不是尺度s的整數(shù)倍,序列尾部可能剩余部分。為了不丟棄這一序列部分,從序列尾部重復(fù)這一分割過程。因此,共分割得到2Ns個區(qū)間。
圖1 消除趨勢波動分析中去趨勢程序示意圖
c. 采用最小二乘擬合法計算序列2Ns個區(qū)間中每個區(qū)間的局部趨勢。定義時間尺度s對應(yīng)的消除趨勢時間序列Ys,i,與原始時間序列和擬合序列不同,消除趨勢時間序列Ys,i表示為相對局部平均的波動:
Ys,i=Yi-yv,i
(2)
式中:yv,i是第v個區(qū)間的擬合多項式。圖1表示了s分別為100和200時采用二次多項式擬合(即DFA2)的情形。擬合過程還可采用線性、三次式或更高階多項式,相應(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. 計算所有區(qū)間的平均值得到整個序列的q階波動函數(shù):
(5)
式中:下標(biāo)q可以取0之外的任何的實值,當(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不同而變化。這是因為對于所有分割區(qū)間v,均方差F2(v,s)相等,且由式(4)表示的平均過程計算得到的波動函數(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方法計算得到的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)增加,意味著分布的奇異性增大。
某工程屬 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)測項目。大壩變形監(jiān)測項目主要有正垂線(PL)、倒垂線(IP)、引張線和視準(zhǔn)線,垂直位移監(jiān)測采用水準(zhǔn)測量。大壩正、倒垂線用于監(jiān)測壩頂及壩體內(nèi)部的水平順河向、垂直水流向位移量,布置方案如圖2所示。為提高可比性,本文選取懸掛點均位于壩頂179.0 m高程的4號壩段(溢流壩段)、2號壩段(擋水壩段)和5號壩段(擋水壩段)的正垂線PL7、PL5和PL6的自動化監(jiān)測時間序列進(jìn)行多重分形分析。為方便分析,順河向水平位移以向下游為正,向上游為負(fù)。可用的時間序列自2002年10月開始,選取2003年1月1日至2008年12月31日間的3個測點的順河向位移序列進(jìn)行分析,實測的環(huán)境量和位移時間序列如圖3所示。
圖3 環(huán)境量和正垂線順河向位移實測時間序列
從圖3可以看出,變形呈一定程度的增加趨勢,但并非單調(diào)增加,而是有漲有落,且變形波動過程呈一定的周期性。計算PL5、 PL6和PL7測點徑向位移序列的標(biāo)準(zhǔn)差、峰度和偏度等基本統(tǒng)計量,結(jié)果見表1,可看出,各序列的偏度都不為 0,且峰度都小于0,因此均屬于非正態(tài)分布。PL7和PL5測點位移序列分布的偏度均大于0,這表示其數(shù)據(jù)分布形態(tài)與正態(tài)分布相比為右偏;PL6測點位移序列分布的偏度小于0,其數(shù)據(jù)分布形態(tài)與正態(tài)分布相比為左偏。同時,所有測點序列的峰度均小于0,這表示測點序列的整體分布與正態(tài)分布相比較為平緩,為寬頂峰。由此可知,測點序列為呈右偏或左偏的非正態(tài)分布。
表1 正垂線徑向位移時間序列的基本統(tǒng)計量
根據(jù)式(1)~(4)計算波動函數(shù)Fq(s),參照文獻(xiàn)[10],這里波動函數(shù)的階數(shù)取值范圍為-6≤q≤6;對于時間尺度s,取值需適中,取值太小沒有實際意義,過大,則分割的區(qū)間數(shù)太少,結(jié)果將與單分形結(jié)果相似,一般地,時間尺度s的范圍為5≤s≤N/4,根據(jù)大壩變形監(jiān)測資料的周期性波動趨勢,這里選取s={6,12,24,48,96,192,384}。圖4為采用MF-DFA1計算得到的q=±6、±4、±2、0時溫度和庫水位及PL5、PL6和PL7測點徑向位移的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)境量序列和各測點位移序列計算得到的廣義指數(shù)h(q)均隨q增加而單調(diào)減少,這意味著環(huán)境量序列和各測點位移序列都存在較顯著的多重分形特征。當(dāng)q=2時,計算得到的溫度和庫水位序列h(2)的值分別為1.31和1.46,各測點位移序列h(2)的值分別為1.38、1.43和1.5。
圖5說明位移序列存在明顯的分形特征,因此有必要進(jìn)行多重分形譜分析。由式(8)和式(9)分別計算得到多重分形標(biāo)度指數(shù)τ(q)和分形譜寬度Δα、奇異譜高度Δf(α),得到的τ(q)-q和f(α)-α曲線分別如圖6和圖7所示。其中,溫度和庫水位序列的分形譜寬度Δα分別為0.612和1.227,奇異譜高度Δf(α)分別為0.715和1.122;PL5、PL6和PL7測點序列的分形譜寬度Δα分別為0.614、0.728和0.679,奇異譜高度Δf(α)分別為0.72、0.79和0.69。
圖5 廣義Hurst指數(shù)和q的關(guān)系
圖6 位移日測值序列質(zhì)量指數(shù)譜曲線
圖7 位移日測值序列多重分形奇異譜曲線
從圖4(a)(b)可以看出,溫度序列的Fq(s)-s曲線在s=96處扭曲,對應(yīng)時間尺度大約為一個季度;水位序列在s=192處扭曲,對應(yīng)時間尺度大約為半年周期。從圖4(c)(d)(e)可以看出,PL5測點序列的Fq(s)-s曲線都在s=96處扭曲;PL6測點序列有所不同,在s=24處接近月周期附近扭曲;PL7測點曲線扭曲發(fā)生在s=384接近年周期處。同時,由圖3(b)可以看出,各測點位移序列大小波動存在一定的差異,F(xiàn)q(s)-s曲線較好地反映了測點序列局部時間內(nèi)和較大尺度范圍內(nèi)的日測值的波動情況。
對測點序列進(jìn)行統(tǒng)計分析,其中水壓因子取3次方,溫度取半年周期,時效取對數(shù)形式[1],統(tǒng)計模型分析得到2007年典型測點年變幅的結(jié)果(表2),計算得到PL5、PL6和PL7測點序列的溫度分量約占總位移年變幅的61.7%、63.5%和63.0%,溫度變化對壩體上下游方向水平位移變化影響最大,其次為庫水位變化;時效分量相對較小,表明壩體變形趨于穩(wěn)定。
表2 統(tǒng)計模型分析結(jié)果 mm
測點序列的h(2)值大于1且接近于1.5,這說明測點序列的波動相關(guān)性和布朗運(yùn)動有一定相似性,換句話說,監(jiān)測序列存在較大的不確定性。同時,各序列h(2)值接近,由文獻(xiàn)[4]可知,整體上大壩結(jié)構(gòu)處于彈性狀態(tài)。由h(q)-q曲線q>0和q<0兩段可以區(qū)分各測點的大小波動的差異。圖5(a)表明溫度序列的大波動標(biāo)度無顯著變化,庫水位序列的大波動標(biāo)度相對較大;溫度和庫水位序列尤其是庫水位序列的小波動標(biāo)度單調(diào)遞減較顯著。圖5(b)表明PL7測點序列波動較顯著,其次為PL6測點系列,PL5測點序列波動較??;各測點序列的大波動標(biāo)度單調(diào)遞減較明顯,這是由于蓄水和運(yùn)行初期,環(huán)境量和壩體自身引起的,各測點遞減的趨勢基本相同,表明整體變化趨勢相同;PL6和PL7測點序列的小波動標(biāo)度單調(diào)遞減趨勢也基本相同;PL5測點序列h(q)序列小波動的標(biāo)度較其他兩測點小,表明PL5測點序列的波動較小。
f(α)-α曲線可以反映位移序列在分形結(jié)構(gòu)上不均勻分布的性質(zhì)。由圖7(a)可見溫度和庫水位序列尤其是庫水位序列的奇異譜呈明顯的右鉤狀,小波動的影響在庫水位序列中占顯著優(yōu)勢。由圖7(b)可見各測點序列尤其是PL5和PL7的奇異譜基本對稱;而PL6測點序列奇異譜則呈不明顯的右鉤狀,即對于PL6測點序列,小波動的影響在日測值序列中稍占優(yōu)勢。如前所述,分形譜寬度Δα的值可以反映大壩性態(tài)演變復(fù)雜度的空間差異,各測點分形譜寬度Δα相差不大,相對而言,PL6測點序列的奇異譜譜寬Δα最大,PL5次之,PL7最小。從圖2可以看出,PL6測點高程較低,水位波動相對壩段壩高在3個測點中最大;PL7測點位于河床,水位波動相對壩段壩高在3個測點中最小。初步推斷,各壩段不同的運(yùn)行機(jī)理是導(dǎo)致壩頂變形測點序列多重分形特征差異的根本原因。
由以上分析和討論可以看出,整體上,大壩結(jié)構(gòu)變形處于正常狀態(tài),壩體處于彈性狀態(tài);對比來看,岸坡5號壩段的PL6測點序列表現(xiàn)出較為復(fù)雜的性態(tài),河床4號壩段的PL7測點序列較為穩(wěn)定。綜上,統(tǒng)計模型方法能夠得到監(jiān)測效應(yīng)量的主要影響因子;MF-DFA方法可以識別監(jiān)測效應(yīng)量序列不同時間尺度上的分形特征和大壩整體性態(tài)時空差異和演變規(guī)律。
a. 各垂線測值序列分別在月、季度和年周期時間尺度上,波動函數(shù)存在一個轉(zhuǎn)折點,反映了各測點測值時間序列在不同時間尺度上的波動特征,說明各壩段運(yùn)行機(jī)理具有一定的差異。環(huán)境量的波動函數(shù)表明溫度和水位存在較為顯著的季度和年周期波動,同時廣義Hurst指數(shù)說明水位大小波動較為顯著;通過統(tǒng)計模型分析表明,溫度分量和水位分量占變形量年變幅的大部分。綜合統(tǒng)計模型分析和MF-DFA方法分析結(jié)果,整體上,溫度和庫水位的大波動決定了壩體變形幅度,水位頻繁漲落控制各壩段尤其是岸坡壩段變形的小波動。
b. 廣義Hurst指數(shù)h(q)均隨q的增加而減小,大壩變形測值序列呈顯著的多重分形特征。由Hurst指數(shù)大于1且接近于1.5說明運(yùn)行初期大壩整體上處于正常線彈性狀態(tài)。測點序列波動具有較大的不確定性,這與大壩運(yùn)行初期受到較為復(fù)雜的環(huán)境、水文和運(yùn)行管理等多重因素影響有關(guān)。
c. 綜合環(huán)境量和測點序列的奇異譜形狀和寬度,可以刻畫大壩不同壩段和局部運(yùn)行性態(tài)的動力學(xué)特征,其物理意義明確,為此通過大壩多個測點的奇異譜寬度分析可反映大壩長期性態(tài)演變的時空差異。進(jì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é)報,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.