童倬慧，賀佳

(湖南大學(xué) 土木工程學(xué)院 建筑安全與節(jié)能教育部重點(diǎn)實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410082)

土木工程結(jié)構(gòu)在服役期間受各種因素的影響，可能出現(xiàn)損傷，而結(jié)構(gòu)損傷的發(fā)生、發(fā)展可以通過結(jié)構(gòu)參數(shù)的變化(如剛度、阻尼)得到體現(xiàn)。因此，利用參數(shù)識(shí)別方法獲取結(jié)構(gòu)參數(shù)，有利于判斷結(jié)構(gòu)狀態(tài)和評(píng)估工作性能。根據(jù)數(shù)字信號(hào)處理方式的不同，系統(tǒng)識(shí)別方法可分為3類：頻域、時(shí)域和時(shí)頻域方法[1-4]。頻域方法主要基于結(jié)構(gòu)模態(tài)參數(shù)，如頻率、振型等展開識(shí)別，此類方法能較好的反映結(jié)構(gòu)的整體狀態(tài)，但是往往對(duì)于局部損傷不敏感。時(shí)域方法直接利用時(shí)域數(shù)字信號(hào)以及各種優(yōu)化算法開展參數(shù)識(shí)別，包括最小二乘估計(jì)、EKF、深度學(xué)習(xí)方法等等。時(shí)頻域方法則通過各種變換，如小波變換、希爾伯特變換等，利用時(shí)域和頻域的聯(lián)合信息開展參數(shù)識(shí)別。由于本文提出的是基于EKF 的參數(shù)識(shí)別方法，因此，以下重點(diǎn)關(guān)注基于EKF 的時(shí)域方法。EKF 方法從系統(tǒng)狀態(tài)方程出發(fā)，將待識(shí)別的結(jié)構(gòu)參數(shù)引入狀態(tài)向量，利用泰勒級(jí)數(shù)進(jìn)行線性化處理，并通過計(jì)算增益矩陣，對(duì)擴(kuò)展的狀態(tài)向量進(jìn)行最小方差估計(jì)。目前，已有許多學(xué)者開展了基于EKF 的參數(shù)識(shí)別研究。例如，利用自適應(yīng)因子，YANG 等[5-6]提出了基于EKF 的時(shí)變參數(shù)識(shí)別方法。LⅠU 等[7]提出了結(jié)合EKF和無跡卡爾曼濾波(UKF)適用于強(qiáng)非線性系統(tǒng)的兩步識(shí)別方法。ZHANG 等[8]將l1范數(shù)正則化引入EKF 方法，提出了復(fù)雜結(jié)構(gòu)的局部損傷識(shí)別方法。HUANG等[9]在ZHANG基礎(chǔ)上將l1范數(shù)正則化改為lp范數(shù)正則化，減少了識(shí)別所需觀測(cè)量。齊夢(mèng)晨等[10]利用EKF 和全局迭代實(shí)現(xiàn)了免模型識(shí)別結(jié)構(gòu)非線性恢復(fù)力。XU 等[11]利用EKF 和多項(xiàng)式擬合同時(shí)識(shí)別了結(jié)構(gòu)非線性恢復(fù)力和結(jié)構(gòu)質(zhì)量。XⅠAO 等[12]提出了基于自適應(yīng)EKF 的軌道不平順及高鐵橋梁頻率識(shí)別方法。雖然以上方法均能實(shí)現(xiàn)結(jié)構(gòu)參數(shù)的有效識(shí)別，但是均需已知作用于結(jié)構(gòu)的外部激勵(lì)，這從某種程度上限制了方法的普適性。因此，一些學(xué)者進(jìn)一步開展了基于EKF 的未知激勵(lì)下的結(jié)構(gòu)參數(shù)識(shí)別研究。例如，基于全局迭代EKF 算法，XU 等[13]提出了子結(jié)構(gòu)參數(shù)和荷載識(shí)別方法。EFTEKHAR 等[14]構(gòu)造了未知外激勵(lì)的狀態(tài)方程，采用兩階段EKF 方法識(shí)別了未知外激勵(lì)和結(jié)構(gòu)狀態(tài)。LEⅠ等[15]運(yùn)用整體卡爾曼濾波法(GEKF)識(shí)別未知地震作用下的結(jié)構(gòu)參數(shù)。LⅠU等[16]提出了一種模態(tài)EKF 算法，識(shí)別了結(jié)構(gòu)參數(shù)和外荷載。LEⅠ等[17]提出了一種基于EKF 的非線性恢復(fù)力免模型識(shí)別方法。HE 等[18]通過引入投影矩陣，提出了基于EKF 的結(jié)構(gòu)參數(shù)和荷載識(shí)別方法。利用子結(jié)構(gòu)理論，HE 等[19]進(jìn)一步將以上方法應(yīng)用于子結(jié)構(gòu)參數(shù)和邊界力識(shí)別。以上提出的未知激勵(lì)下的結(jié)構(gòu)參數(shù)識(shí)別方法，計(jì)算過程相對(duì)較為復(fù)雜，對(duì)未知外激勵(lì)進(jìn)行最優(yōu)估計(jì)時(shí)，都基于最小二乘原理。不同于以上方法，本文在理論上提出了一種相對(duì)清晰直觀的方法，即通過引入初等變換矩陣將未知外激勵(lì)所在的方程進(jìn)行消元，得到不含未知外激勵(lì)的系統(tǒng)觀測(cè)方程組，進(jìn)而通過EKF 識(shí)別結(jié)構(gòu)參數(shù)，并進(jìn)行系統(tǒng)狀態(tài)估計(jì)，最后通過狀態(tài)方程獲取未知外激勵(lì)。該方法可以有效縮減在計(jì)算后驗(yàn)估計(jì)的狀態(tài)向量時(shí)，對(duì)應(yīng)方程組的維度。方法的推導(dǎo)過程具體介紹如下，并通過線性和非線性數(shù)值算例驗(yàn)證了該方法的有效性。

1 理論推導(dǎo)

對(duì)于有m個(gè)自由度，n個(gè)待識(shí)別結(jié)構(gòu)參數(shù)和r個(gè)未知外激勵(lì)的結(jié)構(gòu)而言，其運(yùn)動(dòng)平衡微分方程為：

其中，M表示結(jié)構(gòu)的質(zhì)量矩陣，這里假設(shè)為已知；

引入2m+n維的狀態(tài)向量Z來表示結(jié)構(gòu)狀態(tài)：

式(1)可以進(jìn)一步寫成：

其中：

w(t)表示系統(tǒng)噪聲向量，設(shè)為均值為0 的白噪聲，其協(xié)方差陣E(w(t)w(t)T)=Q(t)。表示影響矩陣μ中與未知外激勵(lì)直接相關(guān)的r×r維子矩陣。

基于EKF 原理，系統(tǒng)狀態(tài)的先驗(yàn)估計(jì)可以表示為：

一般而言，加速度信號(hào)較容易獲取，且相對(duì)可靠，因此，本文假設(shè)實(shí)測(cè)部分加速度，觀測(cè)點(diǎn)個(gè)數(shù)設(shè)為l，則離散化后的觀測(cè)方程可表示為：

其中：yk表示t=k×Δt時(shí)刻實(shí)測(cè)的l個(gè)結(jié)構(gòu)加速度響應(yīng)。L表示與加速度傳感器位置有關(guān)的(l×m)維觀測(cè)矩陣，vk表示測(cè)量噪聲，這里假設(shè)均值為0 的高斯白噪聲，其協(xié)方差矩陣。

可以看出，在觀測(cè)方程(6)中，將未知外激勵(lì)系數(shù)矩陣LΦ左乘初等行變換矩陣T，可使得未知外激勵(lì)向量的系數(shù)矩陣轉(zhuǎn)化成行最簡(jiǎn)型矩陣：

Matlab 中可以通過rref 函數(shù)根據(jù)Gauss-Jordan消元法找到矩陣行最簡(jiǎn)型，從而求出用于消元的初等行變換矩陣T。

不難看出，式(8b)中不含未知外激勵(lì)，因此可以采取EKF，利用分離后的觀測(cè)方程式(8b)和式(5)對(duì)Zk進(jìn)行識(shí)別，假設(shè)式(8a)中觀測(cè)噪聲Rk,2。

根據(jù)HE 等[18]推導(dǎo)出的先驗(yàn)誤差協(xié)方差矩陣Pk+1|k，可得

其中：

根據(jù)EKF 原理，系統(tǒng)狀態(tài)的后驗(yàn)估計(jì)值需以其先驗(yàn)估計(jì)值為基礎(chǔ)，同時(shí)考慮觀測(cè)值與估計(jì)值之間的偏差，這里采用的第k+1 步測(cè)量值為不含未知外激勵(lì)的yk+1,2，則

其中：Kk+1表示EKF增益矩陣，可由下式求得：

此時(shí)，系統(tǒng)的后驗(yàn)估計(jì)誤差協(xié)方差矩陣可由下式計(jì)算得到：

基于以上獲取的系統(tǒng)狀態(tài)的后驗(yàn)估計(jì)值，作用于結(jié)構(gòu)的未知外激勵(lì)未知外激勵(lì)可由下式求得：

圖1 識(shí)別算法流程圖Fig.1 Flow chart of proposed approach for identification

2 數(shù)值算例

2.1 兩層剪切型模型

為了驗(yàn)證本文提出方法的有效性，首先考慮2層剪切型結(jié)構(gòu)。結(jié)構(gòu)質(zhì)量m1=m2=200 kg，各層剛度k1=k2=150 kN/m；采用瑞利阻尼構(gòu)造阻尼矩陣，即C=α1M+α2K，這里取α1=0.73，α2=9.8×10-4。在結(jié)構(gòu)頂層作用隨機(jī)外激勵(lì)，相應(yīng)的結(jié)構(gòu)響應(yīng)根據(jù)狀態(tài)空間法計(jì)算求出，計(jì)算的時(shí)間步長(zhǎng)為0.001 s。假設(shè)觀測(cè)加速度響應(yīng)，考慮5%的噪聲影響。系統(tǒng)噪聲協(xié)方差矩陣R=I，觀測(cè)噪聲協(xié)方差矩陣Q=10-4I，待識(shí)別的結(jié)構(gòu)參數(shù)為結(jié)構(gòu)各層剛度ki(i=1,2)，瑞利阻尼系數(shù)α1和α2，其初始值均假設(shè)為對(duì)應(yīng)真實(shí)值的50%。

利用本文提出的方法，結(jié)構(gòu)參數(shù)識(shí)別結(jié)果如表1所示。可以看出，本文提出的方法可以有效識(shí)別結(jié)構(gòu)參數(shù)。為進(jìn)一步考察參數(shù)的收斂情況，圖2給出了部分參數(shù)的識(shí)別結(jié)果，為便于清晰顯示其收斂情況，圖中采用了不同比例的時(shí)間坐標(biāo)。不難看出，各參數(shù)可以快速、穩(wěn)定地收斂于真實(shí)值附近。由于篇幅限制，這里僅給出了部分參數(shù)的收斂結(jié)果。除了能有效識(shí)別結(jié)構(gòu)參數(shù)以外，本文提出的方法還可以識(shí)別作用于結(jié)構(gòu)的未知外激勵(lì)，其識(shí)別結(jié)果如圖2(c)所示，為便于清晰比較，圖中僅給出了9～10 s 的識(shí)別結(jié)果。從圖中可以看出，荷載識(shí)別值與真實(shí)值吻合得較好。

圖2 2層剪切型結(jié)構(gòu)參數(shù)和未知外激勵(lì)識(shí)別結(jié)果Fig.2 Ⅰdentification of unknown excitation and parameters of the 2-floor shear structure

表1 2層剪切型結(jié)構(gòu)參數(shù)識(shí)別結(jié)果Table 1 Ⅰdentification of the 2-floor shear structure

2.2 平面桁架模型

為了進(jìn)一步驗(yàn)證算法的可行性，本節(jié)以一個(gè)平面桁架結(jié)構(gòu)為例，如圖3 所示。桁架參數(shù)取值為：水平桿長(zhǎng)2 m，斜桿長(zhǎng) 2m，各桿橫截面面積S=2×10-4m2，材料彈性模量E=2×108Pa。采用一致質(zhì)量矩陣來模擬結(jié)構(gòu)的質(zhì)量分布，結(jié)構(gòu)阻尼采用瑞利阻尼假定，其阻尼系數(shù)為α1=6.15，α2=1.46×10-4。為模擬損傷，這里假設(shè)6 號(hào)和7 號(hào)桿件的剛度分別退化了15%和5%。在結(jié)構(gòu)4 號(hào)節(jié)點(diǎn)的豎直方向施加隨機(jī)激勵(lì)，對(duì)應(yīng)的結(jié)構(gòu)響應(yīng)根據(jù)狀態(tài)空間法求出，計(jì)算的時(shí)間步長(zhǎng)為0.001 s。假設(shè)觀測(cè)第1，4，5，7，8 和9 個(gè)自由度上的加速度響應(yīng)，并考慮5%的觀測(cè)噪聲影響。待識(shí)別的結(jié)構(gòu)參數(shù)為各桿件的剛度，類似的，系統(tǒng)噪聲協(xié)方差矩陣R=I，觀測(cè)噪聲協(xié)方差矩陣Q=10-4I，初始值均取對(duì)應(yīng)真實(shí)值的50%。

圖3 平面桁架結(jié)構(gòu)Fig.3 Planar truss structure

運(yùn)用本文提出的方法，結(jié)構(gòu)參數(shù)識(shí)別結(jié)果見表2，由表可見，各參數(shù)識(shí)別結(jié)果誤差較小。圖4給出了部分參數(shù)的收斂情況，可以看出，待識(shí)別的參數(shù)能夠快速、穩(wěn)定地收斂于真實(shí)值附近。同樣地，本文提出的方法還可以有效識(shí)別作用于結(jié)構(gòu)的未知外激勵(lì)，如圖4(c)所示。圖中僅給出了5.0～5.5 s 的識(shí)別結(jié)果，很明顯，識(shí)別值與真實(shí)值吻合得較好。

圖4 平面桁架結(jié)構(gòu)參數(shù)和未知外激勵(lì)識(shí)別結(jié)果Fig.4 Ⅰdentification of unknown excitation and parameters of the planar truss structure

表2 平面桁架結(jié)構(gòu)參數(shù)識(shí)別結(jié)果Table 2 Ⅰdentification of parameters of the planar truss structure

2.3 裝有Dahl模型的非線性狀態(tài)結(jié)構(gòu)

為了討論本文提出的算法在非線性系統(tǒng)識(shí)別中的有效性，本節(jié)考慮一個(gè)底層安裝有MR阻尼器的5 層剪切型結(jié)構(gòu)，如圖5 所示。這里，采用Dahl模型來描述MR 阻尼器的非線性恢復(fù)力[20]，如下式所示：

圖5 裝有MR阻尼的5層模型Fig.5 Five-floor numerical model with MR damper

其中：kd為Dahl 模型剛度系數(shù)；cd為Dahl 模型黏滯阻尼系數(shù)；fd為可調(diào)庫倫摩擦力；f0為初始力；Δx(t)和Δ(t)分別表示結(jié)構(gòu)層間位移和層間速度，r(t)表示滯回位移，由以下微分方程確定：

其中：參數(shù)σd用于控制滯回曲線形狀。

系統(tǒng)參數(shù)取值如下：各層的質(zhì)量和剛度分別為200 kg 和350 kN/m，瑞利阻尼系數(shù)分別為α1=0.53，α2=1.29×10-3，Dahl 模型參數(shù)為kd=35 N/m，cd=250 (N·s)/m，fd=200 N,σd=2 000 s/m，f0=0 N。假設(shè)外激勵(lì)作用于結(jié)構(gòu)頂層，對(duì)應(yīng)的非線性響應(yīng)采用4 階龍格-庫塔法計(jì)算，時(shí)間步長(zhǎng)為0.001 s。取結(jié)構(gòu)第1，3，4 和5 層的加速度響應(yīng)作為觀測(cè)值，并考慮3%的噪聲影響。

待識(shí)別的參數(shù)包括Dahl非線性模型的4個(gè)參數(shù)(kd，cd，fd和σd)，以及剪切型結(jié)構(gòu)的剛度和阻尼系數(shù)，系統(tǒng)噪聲協(xié)方差矩陣R=I，觀測(cè)噪聲協(xié)方差矩陣Q=10-4I，初始值均取對(duì)應(yīng)真實(shí)值的50%。利用本文提出的方法，以上參數(shù)的識(shí)別結(jié)果如表3所示?？梢钥闯?，各參數(shù)的識(shí)別值與真實(shí)值均較為接近。圖6 給出了部分參數(shù)的收斂情況，很明顯，識(shí)別值能穩(wěn)定收斂于真實(shí)值附近。圖7(a)給出了未知外激勵(lì)的識(shí)別情況，可以看出，識(shí)別值與真實(shí)值吻合得較好。此外，基于以上識(shí)別的Dahl 模型參數(shù)，圖7(b)進(jìn)一步給出了該模型的非線性恢復(fù)力，不難發(fā)現(xiàn)，恢復(fù)力的估計(jì)值與真實(shí)值比較接近。

表3 未知激勵(lì)下Dahl模型參數(shù)識(shí)別結(jié)果Table 3 Ⅰdentification of the Dahl model under unknown excitation

圖6 Dahl模型部分參數(shù)識(shí)別結(jié)果Fig.6 Part of identified result of the Dahl model

圖7 未知外激勵(lì)和非線性恢復(fù)力的識(shí)別結(jié)果Fig.7 Ⅰdentification of unknown excitation and nonlinear restoring force

3 結(jié)論

1) 提出了一種基于EKF 的結(jié)構(gòu)參數(shù)和未知荷載的識(shí)別方法，利用初等行變換矩陣與觀測(cè)方程組構(gòu)成的矩陣方程進(jìn)行消元，從而使觀測(cè)方程中不含未知外激勵(lì)，并利用改進(jìn)后觀測(cè)方程和系統(tǒng)狀態(tài)方程對(duì)結(jié)構(gòu)狀態(tài)和外激勵(lì)進(jìn)行識(shí)別。

2) 以數(shù)值算例為主，通過剪切型結(jié)構(gòu)和桁架結(jié)構(gòu)驗(yàn)證了本方法對(duì)于線性系統(tǒng)識(shí)別的可行性，此外，通過裝有Dahl 模型的五層框架數(shù)值算例進(jìn)一步驗(yàn)證了本文方法對(duì)于非線性系統(tǒng)識(shí)別的有效性；相關(guān)的實(shí)驗(yàn)研究將在后續(xù)的研究工作中開展，以驗(yàn)證實(shí)際結(jié)構(gòu)的可行性。