徐 強(qiáng)，劉 博，陳健云，李 靜，徐舒桐

（1. 大連理工大學(xué) 海岸與近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；2. 大連理工大學(xué) 工程抗震研究所，遼寧 大連 116024）

1 研究背景

強(qiáng)風(fēng)和地震等災(zāi)害會(huì)對(duì)建筑物造成不同程度的損傷，導(dǎo)致結(jié)構(gòu)力學(xué)性能發(fā)生改變，進(jìn)而引起結(jié)構(gòu)頻率、阻尼等參數(shù)的改變［1-2］。精確地識(shí)別出結(jié)構(gòu)模態(tài)參數(shù)對(duì)結(jié)構(gòu)健康監(jiān)測(cè)與損傷識(shí)別十分重要［3-5］。

近些年，通過結(jié)構(gòu)振動(dòng)數(shù)據(jù)進(jìn)行系統(tǒng)辨識(shí)及損傷識(shí)別的研究較為廣泛，研究方法主要包括最小二乘識(shí)別［6-8］、卡爾曼濾波［9-10］、小波分析［11-12］及基于人工智能的遺傳算法［13］等。如羅鈞等［7］提出了基于約束最小二乘的識(shí)別方法對(duì)剪切型框架結(jié)構(gòu)的損傷進(jìn)行定位和定量，但測(cè)試過程中存在的噪聲對(duì)識(shí)別結(jié)果造成了一定影響；何浩祥等［9］提出基于靜動(dòng)力凝聚方法的擴(kuò)展卡爾曼濾波對(duì)梁橋結(jié)構(gòu)的模態(tài)進(jìn)行識(shí)別，進(jìn)而診斷結(jié)構(gòu)的損傷程度；Pahlo 等［12］將單元模態(tài)應(yīng)變能的小波系數(shù)差應(yīng)用于簡(jiǎn)支梁的損傷識(shí)別中，根據(jù)小波變換系數(shù)的變化對(duì)結(jié)構(gòu)進(jìn)行損傷定位，但該方法需要建立精準(zhǔn)的有限元模型。最小二乘辨識(shí)是實(shí)際工程中最為常見的一種參數(shù)識(shí)別方法，但系統(tǒng)的輸入輸出信息難免會(huì)存在噪聲，且多為白噪聲，其在整個(gè)頻率上的功率譜相同，而輸入輸出信號(hào)在不同頻段有不同的功率譜。即不同頻段輸入輸出時(shí)程信息的信噪比不同，識(shí)別結(jié)果的可信度也不同［14］。而最小二乘方法對(duì)于各頻段識(shí)別結(jié)果的權(quán)值是均等的，分頻加權(quán)最小二乘在最小二乘方法的基礎(chǔ)上進(jìn)行改進(jìn)。選擇可信度高的頻率段給予較大的權(quán)重，可信度低的給予較小的權(quán)重，提高了識(shí)別精度［15］。

本文采用切比雪夫?yàn)V波器將系統(tǒng)的輸入和輸出信息平均分解到不同頻率段上，分別基于各頻段能量和普通克里金插值法作為權(quán)重，提出分頻段加權(quán)最小二乘識(shí)別法，利用一個(gè)五自由度數(shù)值模型驗(yàn)證參數(shù)識(shí)別的準(zhǔn)確性。將上述識(shí)別方法應(yīng)用于Koyna混凝土重力壩在地震作用下的線彈性和彈塑性損傷模型中，通過壩體觀測(cè)點(diǎn)在不同時(shí)刻識(shí)別的模態(tài)參數(shù)，識(shí)別壩體的損傷的區(qū)域，并采用有限元數(shù)值模擬對(duì)識(shí)別結(jié)果進(jìn)行驗(yàn)證。

2 分頻段加權(quán)最小二乘法識(shí)別

2.1 最小二乘法識(shí)別地震荷載下單自由度系統(tǒng)動(dòng)力方程為：

即

式中：m、c、k分別為結(jié)構(gòu)的質(zhì)量、阻尼和剛度；ag（t）為地震加速度時(shí)程；ξ、ωn分別為結(jié)構(gòu)的阻尼比和固有角速度。

從式（2）可以看出，結(jié)構(gòu)的振動(dòng)情況是由阻尼比ξ和固有角速度的平方ω2n共同決定的。定義結(jié)構(gòu)的單位質(zhì)量剛度＜k＞=k/m，結(jié)構(gòu)的單位質(zhì)量阻尼＜c＞=2ξωn。

對(duì)于單自由度體系，應(yīng)用最小二乘方法識(shí)別系統(tǒng)參數(shù)，相當(dāng)求解以下方程組：

式（3）可以寫成：

其中：

式中N為采樣點(diǎn)數(shù)。

參數(shù)θ的最小二乘估計(jì)解為：

對(duì)于多自由度系統(tǒng)（如圖1）：

式中：q為多自由度系統(tǒng)節(jié)點(diǎn)個(gè)數(shù)；下標(biāo)i代表第i個(gè)節(jié)點(diǎn)；ag(tN)為地震加速度時(shí)程；N為采樣點(diǎn)數(shù)。

采用式（6）可識(shí)別多自由度系統(tǒng)的單位質(zhì)量剛度＜ki＞和單位質(zhì)量阻尼＜ci＞。

2.2 分頻段方法本文采用切比雪夫帶通濾波器將結(jié)構(gòu)的輸入輸出信號(hào)根據(jù)頻率的高低平均分成p個(gè)頻段區(qū)間進(jìn)行濾波，所得到的濾波后的每一個(gè)輸入輸出信號(hào)時(shí)程應(yīng)用式（6）—式（9）計(jì)算多自由度系統(tǒng)的單位質(zhì)量剛度＜ki＞j和單位質(zhì)量阻尼＜ci＞j，j=1，…，p，下標(biāo)j代表第j個(gè)頻段區(qū)間。然后加權(quán)求和。

切比雪夫?yàn)V波器的幅度與頻率的關(guān)系可用下式表示：

圖1 多自由度系統(tǒng)

式中：?為帶通波動(dòng)系數(shù)，為濾波器在截止頻率ω0的放大率；為n階切比雪夫多項(xiàng)式。

將以上低通濾波器進(jìn)行如下變換可得到帶通濾波器：

式中：f0為帶通濾波器中心頻率，分別為帶通濾波器的上、下邊界頻率；BW為濾波器的帶寬，BW=fh-fl。

2.3 加權(quán)方法（1）基于能量加權(quán)方法。取輸入輸出信號(hào)每個(gè)頻段區(qū)間j濾波后的時(shí)程求方差和MDj（代表信號(hào)的能量），通過方差和求各個(gè)頻率區(qū)間的權(quán)重λj：

式中：Xm,j為j頻段第m個(gè)樣本信號(hào)幅值；為j頻段的平均幅值。

式中：θ?為多自由度系統(tǒng)單位質(zhì)量剛度＜ki＞和單位質(zhì)量阻尼＜ci＞的估計(jì)值。

（2）基于普通克里金模型計(jì)算權(quán)重方法。普通克里金插值綜合考慮了空間異質(zhì)性及依賴性，在模型的模擬和預(yù)測(cè)方面取得了良好的體現(xiàn)［16-19］。定義無(wú)偏估計(jì)條件：

式中：E（·）為期望；θ0為估計(jì)值。

定義目標(biāo)J：

式中：Var（·）為方差；φ為拉格朗日乘數(shù)。

通過計(jì)算下式可求出各個(gè)頻率區(qū)間的權(quán)重λj：

2.4 模態(tài)識(shí)別方法通過以上方法識(shí)別多自由度系統(tǒng)單位質(zhì)量剛度＜ki＞和單位質(zhì)量阻尼＜ci＞的估計(jì)值，從而得到系統(tǒng)單位質(zhì)量陣、系統(tǒng)單位剛度陣和系統(tǒng)單位阻尼陣，進(jìn)行模態(tài)計(jì)算識(shí)別。當(dāng)阻尼比較小時(shí)，可直接對(duì)系統(tǒng)單位剛度陣＜K＞做特征值求解得到結(jié)構(gòu)固有頻率和振型。

3 數(shù)值算例

為了對(duì)本文的方法進(jìn)行驗(yàn)證，設(shè)計(jì)一個(gè)五自由度剪切型框架結(jié)構(gòu)，由于計(jì)算單位質(zhì)量剛度，取節(jié)點(diǎn)質(zhì)量為1 kg，阻尼為0.1 N/（s/m），剛度為100 N/m，在水平峰值加速度0.474g、豎直峰值加速度0.312g的Koyna 地震動(dòng)作用下（如圖2），通過有限元數(shù)值分析得到結(jié)構(gòu)各節(jié)點(diǎn)的水平位移、速度及加速度響應(yīng)。分別采用三種方法對(duì)結(jié)構(gòu)的模態(tài)參數(shù)進(jìn)行識(shí)別：（1）方法一。對(duì)系統(tǒng)輸入輸出信號(hào)進(jìn)行最小二乘識(shí)別；（2）方法二。采用切比雪夫帶通濾波器將結(jié)構(gòu)的輸入輸出信號(hào)在1 ～100 Hz 的頻率范圍內(nèi)平均分成5 個(gè)頻率段，進(jìn)行最小二乘參數(shù)識(shí)別，然后以各頻率段的能量作為權(quán)值，進(jìn)行加權(quán)計(jì)算；（3）方法三?；谄胀死锝鸩逯捣椒ㄇ蟪鰴?quán)值，對(duì)信號(hào)進(jìn)行分頻率段加權(quán)最小二乘辨識(shí)。將識(shí)別的結(jié)構(gòu)一階固有頻率與有限元數(shù)值計(jì)算結(jié)果進(jìn)行對(duì)比，如表1所示。

圖2 Koyna地震波的加速時(shí)程

表1 固有頻率識(shí)別結(jié)果

從表1可以看出，上述方法均可以準(zhǔn)確地對(duì)結(jié)構(gòu)的固有頻率進(jìn)行識(shí)別，在對(duì)輸入和輸出信號(hào)進(jìn)行加權(quán)分頻后，固有頻率識(shí)別值的精確度有明顯提高。其中基于普通克里金插值的分頻段加權(quán)最小二乘的識(shí)別誤差均在1%以內(nèi)，最小誤差為0.28%，即基于普通克里金加權(quán)的分頻段最小二乘法對(duì)結(jié)構(gòu)固有頻率的識(shí)別效果最佳。

圖3 一階固有頻率時(shí)程

為了能夠?qū)⒋朔椒ㄓ糜谧R(shí)別時(shí)變系統(tǒng)，對(duì)不同時(shí)刻的模態(tài)參數(shù)識(shí)別情況進(jìn)行描述，將時(shí)間段0.5 s內(nèi)的50個(gè)數(shù)據(jù)點(diǎn)作為識(shí)別點(diǎn)，即采用本文提出的方法識(shí)別結(jié)構(gòu)前9.5 s的一階固有頻率時(shí)程。通過最小二乘識(shí)別和基于普通克里金加權(quán)最小二乘識(shí)別的結(jié)構(gòu)前5 s 固有頻率時(shí)程曲線如圖3 所示。從圖3可以看出，對(duì)信號(hào)進(jìn)行分頻加權(quán)后，固有頻率時(shí)程曲線更加平穩(wěn)，連續(xù)性更好。在0 ～3 s 時(shí)由于地震波處于上升段，幅值非平穩(wěn)較大，固有頻率的識(shí)別值波動(dòng)較大，與理論值有較大的差距，3 s后識(shí)別值逐漸趨于穩(wěn)定，且最后識(shí)別結(jié)果與有限元計(jì)算結(jié)果較為相近。

將結(jié)構(gòu)的輸入、輸出信號(hào)平均分成8個(gè)頻率區(qū)間，分別采用上文的三種方法對(duì)結(jié)構(gòu)的固有頻率進(jìn)行識(shí)別，得到的結(jié)果如表2 所示。表2 與表1 對(duì)比可知，在增加信號(hào)的頻率區(qū)間后，分頻加權(quán)最小二乘方法的識(shí)別精度有明顯的提高。

表2 八階頻率區(qū)間的固有頻率識(shí)別結(jié)果

4 實(shí)例分析

4.1 Koyna 混凝土重力壩結(jié)構(gòu)模態(tài)參數(shù)識(shí)別Koyna 混凝土重力壩是強(qiáng)震作用下發(fā)生破壞的實(shí)例之一，被很多學(xué)者用于結(jié)構(gòu)的動(dòng)力響應(yīng)分析以及結(jié)構(gòu)抗震性能評(píng)價(jià)［20-22］。壩體高度為103.0 m，壩頂和壩底的寬度分別為14.8 m 和70 m。本構(gòu)模型采用混凝土線彈性模型，混凝土彈性模量為31 GPa，泊松比為0.2，密度為2643 kg/m3，有限元模型如圖4所示。計(jì)算荷載為重力和地震荷載。在壩體上游面取5 個(gè)監(jiān)測(cè)點(diǎn)，高程分別為：9.267、31.06、45.6、76.5 和103 m，通過有限元數(shù)值模擬計(jì)算出其在Koyna地震作用下的水平和豎直位移、速度及加速度響應(yīng)，并將響應(yīng)信號(hào)分別代入到模態(tài)識(shí)別的三種方法中，得到識(shí)別的結(jié)構(gòu)模態(tài)參數(shù)。將結(jié)構(gòu)水平向響應(yīng)作為輸出信號(hào)所識(shí)別的壩體五階固有頻率，與有限元模態(tài)分析計(jì)算的壩體十階固有頻率進(jìn)行對(duì)比，發(fā)現(xiàn)所識(shí)別的五階固有頻率與有限元模態(tài)分析的第一階、第二階、第三階、第五階和第八階頻率存在著對(duì)應(yīng)關(guān)系，結(jié)果如表3所示。從表3可以看出，采用分頻段加權(quán)最小二乘識(shí)別方法，對(duì)結(jié)構(gòu)固有頻率的識(shí)別誤差明顯小于最小二乘識(shí)別。其中基于普通克里金加權(quán)的分頻段最小二乘的識(shí)別精確度最高，這一結(jié)果也與第3節(jié)的數(shù)值驗(yàn)證的結(jié)果相一致。分頻段加權(quán)最小二乘識(shí)別對(duì)結(jié)構(gòu)第二階和第五階固有頻率的識(shí)別誤差相對(duì)較大，其他階固有頻率的識(shí)別誤差均在3%以內(nèi)。以上分析說(shuō)明，基于普通克里金法加權(quán)的分頻段最小二乘識(shí)別可以更為準(zhǔn)確地識(shí)別結(jié)構(gòu)的固有頻率。

圖4 Koyna重力壩有限元模型

4.2 Koyna 重力壩在地震響應(yīng)下的損傷區(qū)域辨識(shí)結(jié)構(gòu)發(fā)生損傷會(huì)導(dǎo)致剛度和模態(tài)發(fā)生變化，已有學(xué)者將其應(yīng)用于結(jié)構(gòu)的損傷定位［3］。本文采用分頻加權(quán)最小二乘方法分別對(duì)結(jié)構(gòu)的剛度和模態(tài)進(jìn)行識(shí)別，得出其對(duì)于模態(tài)的識(shí)別效果要優(yōu)于剛度，因此基于其對(duì)固有頻率和振型的識(shí)別結(jié)果對(duì)結(jié)構(gòu)的損傷進(jìn)行描述。

Koyna 重力壩混凝土本構(gòu)模型采用混凝土塑性損傷模型，其密度、彈性模量及泊松比與4.1 節(jié)中線彈性模型相同，膨脹角為36.31°，阻尼比為0.05，初始抗壓強(qiáng)度和抗拉強(qiáng)度分別為24.1 MPa和2.9 MPa，有限元模型如圖4 所示。計(jì)算荷載為重力和地震荷載。在壩體的上游面取與4.1 中線彈性模型相同的監(jiān)測(cè)點(diǎn)，將壩體在垂直方向分成5 個(gè)區(qū)域，計(jì)算各觀測(cè)點(diǎn)在Koyna 地震作用下的位移、速度及加速度響應(yīng)信號(hào)，通過基于普通克里金加權(quán)的分頻最小二乘方法識(shí)別結(jié)構(gòu)的模態(tài)參數(shù)。

表3 固有頻率識(shí)別對(duì)比結(jié)果

為了能夠?qū)⒋朔椒ㄓ糜跁r(shí)變系統(tǒng)，從而描述壩體的損傷情況，將監(jiān)測(cè)點(diǎn)在時(shí)間段0.5 s內(nèi)的50 個(gè)時(shí)程數(shù)據(jù)點(diǎn)作為識(shí)別數(shù)據(jù)，將識(shí)別的壩體前9.5 s 一階固有頻率時(shí)程曲線與有限元計(jì)算的損傷分布對(duì)比，如圖5 所示。從圖5 可以看出，固有頻率時(shí)程曲線在0.03 s 出現(xiàn)了下降趨勢(shì)，通過有限元計(jì)算結(jié)果可知在0.03 s時(shí)壩體后折坡處開始出現(xiàn)損傷區(qū)域；在3.3 s時(shí)，時(shí)程曲線出現(xiàn)了較大幅度的下降，對(duì)比有限元損傷分布可知在3.3 s 后，壩體在后折坡處的裂縫沿水平方向急速擴(kuò)展；在5.7 s 時(shí)程曲線不再下降，對(duì)比有限元損傷分布可知此時(shí)壩體幾乎形成了貫穿的主裂縫區(qū)；在7.2 s 后出現(xiàn)了較大幅度的震蕩，此時(shí)壩體已形成了貫穿的裂縫，裂縫上方的結(jié)構(gòu)為自由振動(dòng)狀態(tài)，因此固有頻率識(shí)別值波動(dòng)較大。通過上述分析可以得出：基于壩體水平方向響應(yīng)信號(hào)所識(shí)別的結(jié)構(gòu)一階固有頻率時(shí)程曲線，可以較好的對(duì)其損傷發(fā)展進(jìn)行描述。

圖5 一階固有頻率時(shí)程與損傷分布的對(duì)應(yīng)關(guān)系

圖6 水平向響應(yīng)識(shí)別的二階至五階固有頻率時(shí)程

重復(fù)上文步驟，得到結(jié)構(gòu)二至五階固有頻率時(shí)程如圖6 所示。對(duì)比圖6 與圖5 可以看出，二階至五階固有頻率時(shí)程和一階固有頻率時(shí)程具有相近的變化規(guī)律，并且固有頻率的階數(shù)越高，曲線在損傷區(qū)域的下降幅度越大，形成貫穿裂縫后曲線的波動(dòng)幅度也越大。

圖7 塑性損傷模型上游水平向歸一化振型

圖8 線彈性模型上游水平向歸一化振型

圖9 塑性損傷模型下游向歸一化振型

圖10 線彈性模型下游向歸一化振型

圖11 塑性損傷模型7觀測(cè)點(diǎn)上游向歸一化振型

圖12 線彈性模型7觀測(cè)點(diǎn)上游向歸一化振型

分別以監(jiān)測(cè)點(diǎn)在初始時(shí)刻與9.5 s 時(shí)的水平響應(yīng)作為輸出信號(hào)，采用基于普通克里金加權(quán)的分頻段最小二乘方法識(shí)別結(jié)構(gòu)的一階歸一化振型，并對(duì)其取絕對(duì)值，如圖7 所示。對(duì)4.1 節(jié)中的Koyna 重力壩線彈性模型進(jìn)行相同的振型識(shí)別，結(jié)果如圖8 所示。通過圖5 中的壩體損傷分布圖可知，壩體主要損傷區(qū)域?yàn)?5.6 ～76.5 m 的壩體后折坡處，從圖7 中可以看出，9.5 s 后識(shí)別的結(jié)構(gòu)一階振型在高程為45.6 和76.5 m 處出現(xiàn)了明顯的拐點(diǎn)，這是由此時(shí)壩體后折坡處出現(xiàn)了貫穿的裂縫，其上、下兩個(gè)監(jiān)測(cè)點(diǎn)的響應(yīng)信號(hào)發(fā)生變化，導(dǎo)致識(shí)別的結(jié)構(gòu)一階振型發(fā)生改變。而圖8 中未發(fā)生損傷的線彈性模型則沒有此變化，即根據(jù)壩體損傷前后歸一化振型識(shí)別值的變化，可以識(shí)別壩體的損傷區(qū)域。

在與上游監(jiān)測(cè)點(diǎn)同一水平線的壩體下游處，取對(duì)應(yīng)的5個(gè)點(diǎn)作為下游的監(jiān)測(cè)點(diǎn)，分別以Koyna混凝土重力壩彈塑性損傷模型與線彈性模型豎直方向響應(yīng)作為輸出信號(hào)，重復(fù)上文對(duì)結(jié)構(gòu)在初始時(shí)刻和9.5 s 時(shí)的歸一化振型識(shí)別，結(jié)果如圖9 和圖10 所示。從圖9、圖10 可以發(fā)現(xiàn)，采用下游觀測(cè)點(diǎn)作為響應(yīng)信號(hào)，同樣可以根據(jù)壩體損傷前后歸一化振型的變化確定其損傷區(qū)域。對(duì)比圖7 和圖9 可看出，采用下游觀測(cè)點(diǎn)的豎直方向響應(yīng)信號(hào)作為系統(tǒng)的輸出信號(hào)，得到的歸一化振型在損傷區(qū)域的斜率變化更為明顯。

在壩體上游面選取7 個(gè)監(jiān)測(cè)點(diǎn)，其高程分別為9.267、31.06、55.3、71.5、76.5、92.8 和103 m。重復(fù)上文中基于普通克里金加權(quán)的分頻段最小二乘方法對(duì)結(jié)構(gòu)初始時(shí)刻和9.5 s的振型識(shí)別，圖11為基于壩體塑性損傷模型的識(shí)別結(jié)果，圖12 為基于壩體線彈性模型的識(shí)別結(jié)果。通過圖11 可以看出，在9.5 s 后結(jié)構(gòu)一階振型在高程為55.3 和71.5 m 處出現(xiàn)了明顯的拐點(diǎn)，即該區(qū)域?yàn)樽R(shí)別的結(jié)構(gòu)損傷區(qū)域。而上文通過5個(gè)監(jiān)測(cè)點(diǎn)所識(shí)別的損傷區(qū)域在45.6 ～76.5 m 之間。說(shuō)明適當(dāng)?shù)脑黾佑^測(cè)點(diǎn)，可以更為精準(zhǔn)的通過振型變化識(shí)別結(jié)構(gòu)的損傷區(qū)域。

5 結(jié)論

本文分別基于各頻段能量和普通克里金插值作為權(quán)值，提出了分頻段加權(quán)最小二乘對(duì)結(jié)構(gòu)的模態(tài)參數(shù)及損傷區(qū)域進(jìn)行識(shí)別。研究結(jié)果表明，在對(duì)系統(tǒng)的輸入輸信號(hào)進(jìn)行分頻加權(quán)后，可以明顯的提高其對(duì)結(jié)構(gòu)固有頻率識(shí)別的精確度，且基于普通克里金加權(quán)的分頻段最小二乘法的識(shí)別效果最佳。以Koyna混凝土重力壩彈塑性損傷模型在地震作用下的響應(yīng)信號(hào)為輸出信號(hào)，通過基于普通克里金加權(quán)的分頻段最小二乘法識(shí)別的結(jié)構(gòu)固有頻率時(shí)程曲線可以準(zhǔn)確地對(duì)結(jié)構(gòu)的損傷情況進(jìn)行描述。根據(jù)壩體結(jié)構(gòu)在損傷前后所識(shí)別的歸一化陣型的改變，可以對(duì)結(jié)構(gòu)的損傷區(qū)域進(jìn)行識(shí)別，且以下游監(jiān)測(cè)點(diǎn)豎直向響應(yīng)信號(hào)作為輸出信號(hào)的識(shí)別效果最佳。