混凝土壩強震觀測資料的ARX模型分析

鄭東健,羅德河,仇建春

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室,江蘇 南京 210098;2.河海大學(xué)水利水電學(xué)院,江蘇 南京 210098; 3. 珠江水利科學(xué)研究院,廣東 廣州 510611)

我國是世界上遭受地震災(zāi)害比較嚴(yán)重的國家之一。為了減少地震帶來的損害,需要科學(xué)地認(rèn)識強震作用及其對水工結(jié)構(gòu)工程的影響,對水工結(jié)構(gòu)進行強震觀測就是解決上述問題的重要途徑。迄今為止,世界范圍內(nèi)的混凝土壩已有不少因遭遇強烈地震而破壞的實例。1962年3月19日我國廣東新豐江水庫發(fā)生了MS6.1級強烈地震,使其混凝土大頭壩頭部發(fā)生水平裂縫[1];1967年12月11日印度的Koyna重力壩遭遇了MS6.5級水庫誘發(fā)地震,使其下游折坡點附近的上、下游壩面出現(xiàn)大量的水平裂縫[2];1971年2月9日和1994年1月17日,美國Pacoima拱壩分別遭遇了MS6.6級和MS6.8級強烈地震,造成左岸壩體與重力墩之間橫縫張開[3-4];1990年6月21日伊朗的Sefid Rud大頭壩遭遇了MS7.6級大震,壩體受到了嚴(yán)重?fù)p傷,上、下游壩面出現(xiàn)大量的水平裂縫[5];1999年9月21日臺灣石崗閘壩遭遇MS7.3級強烈地震,北段三跨泄洪道斷塌,斷裂處錯動達8 m[6];2008年5月12日我國沙牌碾壓混凝土壩遭遇MS8.0級汶川大地震,拱壩僅壩頂電梯井房發(fā)生破壞,其他部位表現(xiàn)良好,沒有明顯的損傷開裂等現(xiàn)象發(fā)生。以上工程震害相對明顯,可以通過目視檢查發(fā)現(xiàn)。但結(jié)構(gòu)內(nèi)部、水下等隱蔽部位損傷情況常難以目視檢查,而一般損傷情況下滲流和變形觀測資料也變化不明顯,如何及時把握地震過程結(jié)構(gòu)狀態(tài)變化是目前面臨的巨大挑戰(zhàn)。隨著我國西南地震高烈度區(qū)一批高壩大庫的完建運行,及時分析強震觀測資料,在線快速評估地震過程大壩結(jié)構(gòu)狀態(tài)變化,對確保大壩安全有重要意義。目前DL/T 5416—2009《水工建筑物強震動安全監(jiān)測技術(shù)規(guī)范》第9條加速度記錄的處理分析,規(guī)定根據(jù)強震觀測系統(tǒng)“各個測點的記錄和預(yù)存的抗震設(shè)計計算求得的各個測點能抗御的最大加速度值,對大壩進行安全評估”(以下簡稱“規(guī)范法”),并規(guī)定“應(yīng)對加速度記錄進行常規(guī)處理分析”,其內(nèi)容包括校正加速度記錄、速度和位移時程、5個阻尼值(0,0.02,0.05,0.1,0.2)反應(yīng)譜和傅里葉譜等。SL 486—2011《水工建筑物強震動安全監(jiān)測技術(shù)規(guī)范》第7條加速度記錄的處理分析,規(guī)定與DL/T 5416—2009類似,這些條款對水工建筑物的震后評價起重要作用。但由于工程的復(fù)雜性以及運行環(huán)境的影響,設(shè)計模型和參數(shù)與工程實際情況會存在差異,根據(jù)設(shè)計計算的測點能抗御最大加速度評估大壩震后安全有其不足。為此,本文直接應(yīng)用混凝土壩強震觀測資料,在強震加速度記錄常規(guī)處理分析基礎(chǔ)上,研究了強震資料時間序列ARX(autoregressive exogenous)模型,通過跟蹤ARX模型的時變參數(shù)變化,提出混凝土壩結(jié)構(gòu)狀態(tài)震后評價方法,并采用遺忘因子提高模型效率,實現(xiàn)混凝土壩強震觀測資料的在線快速評估。

1 ARX模型

ARX模型可以看作系統(tǒng)的動力學(xué)方程,其參數(shù)反映了系統(tǒng)的動力特性[7]?；炷翂螐娬鹩^測資料是混凝土壩在受地震作用下的振動響應(yīng)。當(dāng)?shù)卣鹱饔玫膹姸容^小時,一般不會對結(jié)構(gòu)的狀況產(chǎn)生影響,此時結(jié)構(gòu)系統(tǒng)獲得的動力響應(yīng)較為平穩(wěn),強震觀測時間序列為平穩(wěn)時間序列;當(dāng)?shù)卣鹱饔脧姸容^大時,混凝土壩可能會進入非彈性階段并出現(xiàn)損傷,觀測資料為非平穩(wěn)時間序列。用強震觀測資料時間序列建立ARX模型,可以實現(xiàn)用較少參數(shù)反映序列所蘊含的結(jié)構(gòu)動力特性,識別混凝土壩結(jié)構(gòu)狀態(tài)的變化。

本文所采用時間序列ARX模型如下:

式中:y(t)為混凝土壩強震測點采集的信號;u(t)為混凝土壩基礎(chǔ)或周邊場地強震儀采集的信號;ξ(t)為與誤差和模型不確定性有關(guān)的白噪聲信號;t為時間；ai、bj為系數(shù);na、nb為系統(tǒng)模型的階次,可根據(jù)AIC準(zhǔn)則[8-9]確定。

為判斷混凝土壩結(jié)構(gòu)在不同時間段的結(jié)構(gòu)參數(shù),在式(1)基礎(chǔ)上,將結(jié)構(gòu)參數(shù)表示為時間的函數(shù),即:

y(t)=φT(t)θ(t)+ξ(t) (t=1,2,…,N)

(2)

其中φT(t)=(-y(t-1),-y(t-2),…,

-y(t-na),u(t),u(t-1),…,u(t-nb))

θ(t)=(a1(t),a2(t),…,ana(t),b0(t),b1(t),

b2(t),…,bnb(t))T

(3)

其中Y=(y(1),y(2),…,y(N))T

φ=(φT(1)φT(2) …φT(N))T

(4)

P(t)=[I-K(t)φT(t)]P(t-1)

式中：I為單位矩陣。

在采集到結(jié)構(gòu)系統(tǒng)輸入輸出數(shù)據(jù)的同時,為實現(xiàn)結(jié)構(gòu)時變參數(shù)的實時在線識別,跟蹤混凝土壩結(jié)構(gòu)系統(tǒng)不確定性規(guī)律的參數(shù)變化且提高計算效率,可引入遺忘因子f(t),結(jié)合遞推最小二乘算法實現(xiàn)時變參數(shù)的識別,其主要公式如下[10]:

(5)

式中：f(t)為遺忘因子,計算過程中需滿足0

P(0)=αI

(6)

(7)

式中:α為充分大的正實數(shù)(104～1010);ε為零向量或充分小的正實向量。

2 實例驗證

為驗證ARX模型時變參數(shù)對混凝土壩結(jié)構(gòu)狀態(tài)評價的有效性,采用印度Koyna地震波對西部某混凝土拱壩進行結(jié)構(gòu)動力有限元分析,并將拱壩典型位置動力響應(yīng)作為仿真強震資料。該壩為混凝土雙曲拱壩,水庫正常蓄水位1 880 m,電站裝機容量3 600 MW,壩頂高程1 885.00 m,壩基最低建基面高程1 580.00 m,最大壩高為305.0 m,壩頂寬度16.0 m,壩底厚度63.0 m。圖1為拱壩的有限元模型,模型上游、下游和地基取1倍的壩高。模型采用ABAQUS軟件8節(jié)點六面體實體單元C3D8R進行網(wǎng)格劃分,整個模型共有61 909個單元,68 448個節(jié)點,其中壩體單元4 708個。

圖1 某拱壩有限元模型

為簡化計算,將壩基簡化為相同材料參數(shù)。結(jié)構(gòu)動力分析時采用動彈模,即在靜彈模的基礎(chǔ)上提高30%作為動彈模進行分析[14],采用無質(zhì)量地基進行模擬。為了模擬地震作用過程中壩體損傷的發(fā)生和發(fā)展,采用混凝土塑性損傷模型(CDP模型)進行分析,并通過損傷因子直觀的得到結(jié)構(gòu)的損傷位置與程度。壩體和壩基動彈性模量分別取32 GPa和20 GPa,密度分別取2.4 t/m3和2.7t/m3,泊松比分別取0.167和0.25;壩體混凝土材料塑性參數(shù)膨脹角取30°,偏心率取0.1,雙軸極限抗壓強度與單軸受壓極限強度之比fb0/fc0取1.16,拉伸子午面上與壓縮子午面上的第二應(yīng)力不變量之比K取0.666 7,黏性參數(shù)取0.000 5。

混凝土塑性損傷模型的損傷因子采用張勁公式法[15],該方法將規(guī)范提供的混凝土本構(gòu)模型與CDP模型統(tǒng)一起來,并提出損傷因子dk可表示為

(8)

式中:t、c分別為拉伸和壓縮情況;β為塑性應(yīng)變和非彈性應(yīng)變之比,在受壓情況時取值為0.35～0.7,在受拉情況時取值為0.5～0.95;εin為混凝土在拉壓情況下的非彈性階段應(yīng)變；E0為壩體初始(無損)彈性模量。dk=0表示無損傷,dk=1表示完全損傷。

進行混凝土拱壩損傷有限元模擬時,輸入的地震波采用有完整記錄的印度Koyna地震波,取其順河向地震波作為輸入波,地震加速度時程曲線如圖2所示。數(shù)值仿真模擬的采集時長取10 s,采樣頻率取50 Hz。為了分析不同地震強度下拱壩結(jié)構(gòu)狀態(tài)變化,將加速度峰值分別調(diào)整為0.1g和0.2g的兩條Koyna地震波時程曲線,作為結(jié)構(gòu)的地震激勵輸入,用ABAQUS有限元軟件模擬相應(yīng)地震強度下拱壩損傷的發(fā)展,并分析其與ARX模型參數(shù)變化的對應(yīng)關(guān)系。

圖2 水平向地震加速度時程曲線

為了提取各地震工況拱壩的動力響應(yīng)時間序列,在拱冠壩頂、1/2壩高、左右岸1/4拱壩處分別設(shè)置動力響應(yīng)采集測點,編號為1號～4號,具體位置如圖3所示。建立ARX時間序列模型時,采用下游河床壩基位置的一個測點響應(yīng)信息作為ARX時間序列模型的響應(yīng)輸入。

圖3 拱壩強震測點仿真布置

下面通過模擬兩種不同工況下的損傷情況,分析各測點加速度響應(yīng)時間序列建立的ARX模型參數(shù)和結(jié)構(gòu)狀態(tài)變化的對應(yīng)關(guān)系。

a. 峰值加速度為0.1g。當(dāng)結(jié)構(gòu)輸入加速度峰值為0.1g的地震激勵時,有限元模擬仿真的結(jié)果表明,壩體總體上保持完好,沒有發(fā)生損傷。壩基測點、1號、2號、3號測點加速度時程曲線見圖4。將壩基位置測得的加速度響應(yīng)信息作為輸入信息,各個測點的加速度響應(yīng)信息作為輸出建立ARX模型(na,nb分別取6和2)。圖5給出了1號測點ARX模型的時變參數(shù)辨識結(jié)果(2號～4號測點類似)。從圖5可以看出,測點的響應(yīng)時間序列計算的系統(tǒng)參數(shù)值在整個地震期間保持穩(wěn)定不變,說明結(jié)構(gòu)系統(tǒng)狀態(tài)沒有發(fā)生變化,即拱壩在小震情況下沒有發(fā)生損傷,與仿真分析結(jié)果一致。由于遞推迭代原因,迭代初始階段參數(shù)變化可能較大,可通過預(yù)先加1～2 s白噪聲激勵(幅值與地震波初始振幅接近)來消除。

圖4 測點響應(yīng)時間序列

圖5 ARX模型參數(shù)辨識結(jié)果

圖6 峰值加速度為0.20g時壩體損傷狀態(tài)

b. 峰值加速度為0.20g。為了使拱壩出現(xiàn)損傷,將峰值加速度為0.2g的Koyna地震波作為地震動作用輸入,結(jié)構(gòu)的破壞情況如圖6(a)所示。從圖中可以看出拱壩的損傷區(qū)域集中在拱冠、兩岸壩肩及壩踵處,拱冠位置上下游面均發(fā)生破壞。其中拱壩的中上部位區(qū)域破壞的較為嚴(yán)重。在地震作用的后期,破壞區(qū)域逐漸向兩側(cè)和壩體下方擴展。本文模擬的大壩損傷結(jié)果與文獻[16]結(jié)果較為一致。三個損傷區(qū)域的拉伸損傷因子時程圖如圖6(b)所示,可以清晰地看出各區(qū)域的損傷情況,此時結(jié)構(gòu)發(fā)生破壞的時間不集中,并且損傷因子值呈階梯狀增長。如壩踵部位,4.26 s首先發(fā)生損傷,到4.9 s時損傷因子值為0.268,但損傷程度較小,到6.2 s時損傷因子值突變?yōu)?.775;左岸壩肩1/4拱位置在7.44 s時損傷因子突變?yōu)?.849;拱冠位置4.62 s時損傷初始,6.52 s時損傷因子值為0.191。將壩基位置測得的加速度響應(yīng)信息作為輸入信息,各個測點的加速度響應(yīng)信息作為輸出建立ARX模型。ARX模型時變參數(shù)的辨識結(jié)果如圖7所示，可以看出,1號、3號和4號測點ARX模型參數(shù)在4～5 s之間發(fā)生變化,尤其1號測點突變明顯;在6～7 s間各測點ARX模型參數(shù)再一次出現(xiàn)較明顯變化,2～4號測點突變最明顯;隨后各測點模型參數(shù)仍有不同程度變化,但幅度有所減小。對比圖6(b)和圖7可見模型參數(shù)突變時段損傷因子值也發(fā)生了突變,且不同位置測點識別出的結(jié)構(gòu)損傷發(fā)生時間與仿真模型損傷因子值突變時間基本一致,表明混凝土壩強震資料的ARX模型參數(shù)具有較好的結(jié)構(gòu)狀態(tài)識別能力,且壩頂?shù)葎恿憫?yīng)較大的位置測點對損傷初始判別更為敏感。

圖7 ARX模型參數(shù)辨識結(jié)果

3 結(jié) 語

我國許多高混凝土壩建在高地震烈度區(qū),在經(jīng)受各種程度地震動荷載作用下,這些高壩是否出現(xiàn)損傷直接關(guān)系到大壩安全和震后應(yīng)急措施。本文直接利用實測強震觀測資料,通過建立ARX模型，利用模型參數(shù)反映結(jié)構(gòu)動力特性的原理，結(jié)合遺忘因子遞推最小二乘算法,可以實現(xiàn)混凝土壩震后結(jié)構(gòu)狀態(tài)的快速評估?；炷凉皦蔚卣鹱饔梅抡娣治霰砻?ARX模型參數(shù)的時變過程與損傷因子的變化過程較為一致,驗證了時變參數(shù)識別方法的有效性。強震觀測資料的ARX模型分析法具有快速和便于在線評價的特點,可以作為規(guī)范法的補充。