混凝土壩強(qiáng)震觀測(cè)資料的ARX模型分析

鄭東健,羅德河,仇建春

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

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

1 ARX模型

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

本文所采用時(shí)間序列ARX模型如下:

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

為判斷混凝土壩結(jié)構(gòu)在不同時(shí)間段的結(jié)構(gòu)參數(shù),在式(1)基礎(chǔ)上,將結(jié)構(gòu)參數(shù)表示為時(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ù)的同時(shí),為實(shí)現(xiàn)結(jié)構(gòu)時(shí)變參數(shù)的實(shí)時(shí)在線識(shí)別,跟蹤混凝土壩結(jié)構(gòu)系統(tǒng)不確定性規(guī)律的參數(shù)變化且提高計(jì)算效率,可引入遺忘因子f(t),結(jié)合遞推最小二乘算法實(shí)現(xiàn)時(shí)變參數(shù)的識(shí)別,其主要公式如下[10]:

(5)

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

P(0)=αI

(6)

(7)

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

2 實(shí)例驗(yàn)證

為驗(yàn)證ARX模型時(shí)變參數(shù)對(duì)混凝土壩結(jié)構(gòu)狀態(tài)評(píng)價(jià)的有效性,采用印度Koyna地震波對(duì)西部某混凝土拱壩進(jìn)行結(jié)構(gòu)動(dòng)力有限元分析,并將拱壩典型位置動(dòng)力響應(yīng)作為仿真強(qiáng)震資料。該壩為混凝土雙曲拱壩,水庫(kù)正常蓄水位1 880 m,電站裝機(jī)容量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é)點(diǎn)六面體實(shí)體單元C3D8R進(jìn)行網(wǎng)格劃分,整個(gè)模型共有61 909個(gè)單元,68 448個(gè)節(jié)點(diǎn),其中壩體單元4 708個(gè)。

圖1 某拱壩有限元模型

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

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

(8)

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

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

圖2 水平向地震加速度時(shí)程曲線

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

圖3 拱壩強(qiáng)震測(cè)點(diǎn)仿真布置

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

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

圖4 測(cè)點(diǎn)響應(yīng)時(shí)間序列

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

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

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

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

3 結(jié) 語(yǔ)

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