基于形函數方法快速識別結構動態荷載的試驗驗證

張青霞，段忠東，Lukasz Jankowski，王 豐

(1.大連民族學院 土木建筑工程學院，大連 116600;2.哈爾濱工業大學 深圳研究生院，深圳 518055;3.哈爾濱工業大學 土木工程學院，哈爾濱 150090;4.Smart－Tech Centre，Institute of Fundamental Technological Research，Polish Academy of Sciences，02 －106，Warsaw，Poland)

準確的結構動態荷載信息是確定結構動態行為的關鍵因素之一，對于結構損傷識別、安全評定等具有重要的作用。此外，它為災難性事故的分析評估提供了科學依據，有益于預防和避免類似事故的重發生。荷載識別屬于結構動力學中的反問題分析，由結構的動力響應和其動態特性確定。

目前的荷載識別方法［1－4]中，大多數屬于頻域［5，6]、時域［7，8]或小波域［9]內分析，基本思想是通過計算實測的結構響應和提前獲取的系統脈沖響應函數的反卷積計算荷載，也就是所謂的反卷積方法［2]。此外，還有時間有限元法［10]，基于魯棒觀測器(robust observer)法［11]，逆結構濾波法(ISF)［12]等。每種方法具有不同的優缺點，但均面臨反問題分析中常存在的病態性問題，如解的唯一性、穩定性等。反卷積法以其形式簡單、可靈活采用時域和頻域內分析的特點應用較廣泛;但若測量時間長或采樣頻率高時，則導致脈沖響應函數(傳遞函數)組成的系數矩陣維數很大，從而使該方法數值求解困難。

文獻［8] 提出荷載形函數方法來改善荷載識別反卷積方法中的不足，利用有限元中的形函數概念逼近動態荷載曲線，從而把識別未知荷載的時間歷程轉變為識別數目有限的荷載形函數的權重或系數，大大降低了系數矩陣的維數，提高計算效率。另外，荷載形函數是連續的，在一定程度上能夠過濾噪聲的影響，提高逆問題對噪聲的魯棒性。本文對荷載形函數的相關參數設置進行了理論補充，并通過試驗驗證了該方法的有效性。

1 反卷積方法

假定所考慮的結構為線彈性而且在荷載作用期內為小變形即忽略幾何非線性，那么結構的響應線性依賴于外部荷載作用。設結構初始狀態為零，根據杜哈姆積分，結構的響應可由外部荷載通過卷積表示為:

式中:α為測點編號;fi為作用在結構上的第i個外部荷載，hαi(t)為測點α和荷載fi之間的脈沖響應函數。

實際應用中，結構的響應一般通過測量或數值模擬得到，從而是離散化的，因此考慮測量時間內所有的離散的結構測點響應，并按一定順序排列，式(1)可以表示為如下的矩陣形式:

式中yM匯集了所有測點在測量時間內的離散響應;f匯集了所有外部荷載的離散時間歷程;矩陣h是分塊矩陣，子矩陣是由測點和外部荷載之間的相關脈沖響應函數構成的Toeplize矩陣。

由式(1)、(2)知，通過解實測結構響應和系統脈沖響應函數計算反卷積可以得到未知荷載，即求解方程(2)。這是一個線性方程組，為保障所求解的唯一性，需要獨立測點(傳感器)的數目至少等于可能的未知荷載數目。在實際應用中，需要預先知道或假定荷載可能存在的位置信息以確定傳感器的數目和構造對應的脈沖響應函數。

解式(2)就是解線性方程組Ax=b或求解min Ax=b2，往往面臨病態問題。實測響應中一個小的擾動，比如無法避免的測量誤差，將會引起很大的荷載識別誤差，即方程的解對擾動非常敏感［3]。因此，有必要對解進行數值正則化，常用的正則化技術有截斷奇異值分解方法(TSVD)或 Tikhonov方法［3]。

另外，式(2)中方程組和未知數的數目與時間步數成正比，當采樣頻率高或測量時間較長時，比如離線荷載識別中，線性方程組可能會很大，導致計算操作難以實現。利用提出的荷載形函數方法可以改善荷載識別方法中存在的病態問題，減小計算工作，提高識別效率。

2 荷載形函數方法

2.1 基本原理

荷載形函數的基本思想源于有限元中的形函數。不同的荷載形式決定了不同類型的形函數類型，如荷載只作用在一點時，其時間歷程是一個曲線，可以采用梁單元的形函數，本質上是采用曲線插值的方法對荷載曲線進行逼近。如果結構上作用一個線荷載，其時間歷程就是一個曲面，需要選用板的形函數進行曲面插值。對于面荷載，則需要借用實體單元的形函數進行插值。因此選擇形函數類型時，需要根據經驗判斷荷載的類型。

本文假定荷載的時間歷程是一個曲線，如圖 1，可以比擬為一個“時間梁”的變形，其中荷載的時間軸比擬為

圖1 荷載的時間歷程曲線Fig.1 The load time history

“梁”長度方向的位置坐標。只考慮梁的彎曲變形，即每個結點有豎向和轉動兩個位移。假如把荷載歷程按時間劃分為四段，那么該“時間梁”單元有五個結點，十個形函數，定義為荷載形函數，如圖2，其意義為“時間梁”某個結點在豎向或轉動方向發生單位位移引起的荷載歷程。記荷載形函數矩陣N=［N1N2…Nnl] ，nl是形函數的個數。

圖2 荷載形函數Fig.2 Load shape function

根據形函數的意義可知，動態荷載時間歷程f(t)可由若干個形函數近似表示:

式中:α=［α1α2… αnl]T，αi是對應第 i個形函數的“時間梁”結點的位移，稱為該形函數的系數。

2.2 形函數的選取

圖2給出了形函數的形式，具體表達式參考有限元中梁單元形函數。在實際應用中，荷載歷程一般被均勻劃分為若干段。設相鄰每段兩個結點的時間長度定義為形函數的半個周期Tf/2，則形函數的頻率為ff=1/Tf，兩個結點之間為lnt=Tf/(2Δ)=fs/(2ff)個時間步，其中fs為采樣頻率，采樣時間間隔為Δ=1/fs。若分析的離散荷載時間為t，共含nt個時間步，則形函數的個數nl:

由上可知，形函數內時間步數取決于采樣頻率fs和形函數的頻率ff，而形函數個數取決于形函數的頻率ff和測量時間t。測量時間一定時，形函數個數由形函數頻率決定，形函數的頻率選取會影響識別精度。結構動態荷載往往含有多個頻率成分，為能夠準確逼近動態荷載曲線，要求形函數的頻率大于荷載所分析的最高截斷頻率。形函數的頻率小于荷載的截斷頻率時，它不能準確地逼近荷載;而形函數的頻率較高時，取極限類似單位脈沖，識別結果對噪聲敏感。另外，較高的形函數頻率意味著形函數個數多，這時系數矩陣的維數大，計算效率低。因此，選取合適的形函數頻率是基于形函數法識別荷載的關鍵之一。對未知荷載而言，鑒于線性結構的響應線性依賴于外部荷載，確定其所需的形函數頻率可以通過對結構實測響應的頻譜分析決定。

設F(ω)代表響應的傅里葉變換，定義下式:

對于一般荷載，可取r(ωc)=95%對應的頻率ωc為截斷頻率，視ωc為荷載可能的最高頻率，即把ωc作為逼近未知荷載所需的形函數的參考頻率。實際操作中根據荷載形式，r(ωc)的取值可適當增減。那么兩個相鄰結點含的時間步lnt=fs/2ωc，代入式(4)可確定所需的形函數個數nl。

2.3 荷載識別

實際應用中，根據式(3)，式(2)中離散的荷載歷程可以表示為:

式中:Nf=Inf×nf?N 是 nfnt× nfnl的矩陣，N 是 nt× nl維的離散的荷載形函數矩陣。是nfnl維的列向量，由各個荷載的形函數的系數組成，αi(j)是荷載i的第j個形函數的系數。

把式(7)代入式(2)得到:

式(8)中的脈沖響應函數矩陣、荷載形函數矩陣均可以提前構造。形函數的個數nl一般遠小于測量時間步數nt，因此把式(2)直接計算荷載轉化為先求解式(8)中荷載形函數的系數，然后計算荷載，使未知數個數大大減少，顯著提高計算效率。另外，荷載形函數是連續光滑的，根據式(8)計算荷載在一定程度上能削減噪聲的影響。

利用式(8)計算荷載，其系數矩陣的維數得以減小，采用截斷奇異值分解方法(TSVD)對解進行數值正則化更簡單快捷，即:

式中:ui是酉矩陣 U=［u1，…，unfnt] 的第 i個列向量，vi是酉矩陣V=［v1，…，vnfnl] 的第i個列向量，酉矩陣U，V由式(8)中系數矩陣B的奇異值分解得到。σi是矩陣B的第i個正奇異值。q指選用前q個正奇異值求解。奇異值截斷數目可以根據L－曲線計算。通過L－曲線的拐點來權衡響應殘差和荷載光滑度。

根據公式(9)，系數矩陣很小的奇異值會放大結構測量中的響應誤差的影響，從而使識別的結果嚴重偏離真實值。利用反卷積方法識別荷載時，公式(2)中系數矩陣由脈沖響應函數組成，而脈沖響應的振動一般比較劇烈，頻率成分含量比較豐富，特別是高頻成分，系數矩陣中小的奇異值可認為是由脈沖響應中的高頻成分導致的，因此在計算反卷積時需要通過截斷奇異值分解等正則化方法去掉較小的奇異值，以改善解的識別精度。當采用荷載形函數識別荷載時，系數矩陣由脈沖響應函數組成的矩陣和形函數矩陣相乘得到(見式(8))，該系數矩陣中的每列元素相當于對應每個荷載形函數引起的結構響應。因為荷載形函數是連續光滑的，它引起的結構響應也相對光滑，相對于單位脈沖引起的結構響應而言，它含有高頻成分很少或幾乎沒有，使系數矩陣中只存在少數幾個或不存在特別小的奇異值，因此采用荷載形函數識別荷載時，系數矩陣只截斷少數幾個奇異值或不需截斷奇異值。

另外，已知結構的有限元模型，荷載形函數方法結合“移動時間窗”［13]可以推廣到在線識別，也就是式(8)中的實測響應yM被替代為，其中，yM(n)是第n個時間窗的采樣數據;假定前面時間窗內的荷載已被識別，該時間窗內的結構初始狀態和非零初始狀態引起的橋體的自由振動根據結構模型易求。

3 數值算例

首先借助一個兩跨連續梁模型(圖3)來驗證所提出的荷載識別方法的有效性。

梁每跨5 m，工字型組合截面55 mm×100 mm×7.2 mm ×4.5 mm。彈性模量 E=210 GPa，密度 ρ=7 800 kg/m3，截面面積為1 200 mm2，受彎慣性距為1 943 800 mm4。有限元梁模型被均勻劃分為20個單元，在距左端2.26 m、5.26 m的截面底部分別布一個應變傳感器，記為 S1，S2。取前兩階瑞利阻尼比為0.01。

假定結構承受兩個豎向荷載f1，f2(圖4)，分別作用在距梁左端4.5 m、6.5 m截面上。結構動力響應計算采用Newmark逐步積分方法，Newmark參數α=0.25，β=0.5。積分步長0.001 s，即采樣頻率1 000 Hz。測量時間為0.9 s。考慮響應含有5%的高斯隨機噪聲。

假定外部荷載未知，采用荷載形函數法進行識別，并與荷載識別的基本方法—反卷積法(解式(2))進行比較。鑒于荷載識別等反問題往往是病態的，采用TSVD正則化方法求解。奇異值截斷的個數k即正則化程度通過L－曲線確定。

考慮5%的高斯噪聲，反卷積法識別荷載的L－曲線(圖5)的拐點對應k=470，也就是小于最大奇異值0.14%的奇異值都被截斷，識別的相應荷載(見圖6)與真實值相比有明顯振蕩，說明直接利用反卷積方法識別荷載對噪聲敏感。

利用形函數進行荷載識別時，荷載形函數通過測點響應進行傅里葉變換確定。兩個測點響應的各頻率含量相似，以測點S1的傅里葉變換為例進行分析，如圖7。選取ωmax=100 Hz，根據式(6)得到的形函數的頻率ωc=14.43 Hz，見圖8，根據此確定每個荷載需要62個形函數。考慮5%的高斯噪聲，得到的L－曲線如圖9。該曲線沒有典型的L－曲線拐點，說明解式(8)條件數好，k=1即可得到不錯的識別結果，如圖10，與真實值非常接近。

從上面的分析可知荷載識別的主要計算工作來自系數矩陣的奇異值分解，對一個m×n維的矩陣(m≥n)的計算量為O(mn2)階［14]。基本反卷積法的系數維數矩陣為1 802×1 802，利用形函數方法后變為1 802×124，因此這里形函數方法能使計算量降低幾乎兩個數量級。基本反卷積方法識別荷載對噪聲敏感，需要通過合適的奇異值截斷來進行，而且仍存在明顯的識別誤差。形函數方法不但降低計算量，而且能夠在一定程度上過濾噪聲的影響，識別結果抗噪能力強。

4 試驗驗證

圖11 試驗裝置Fig.11 Experimental setup

試驗裝置如圖11所示，試件為一個懸臂鋁梁，長136.15 cm，矩形截面 2.7 cm ×0.31 cm。固定端夾在一個穩定框架上。梁的彈性模量為70 GPa，密度是2 700 kg/m3。采用壓電激勵器APA(Amplified Piezo Actuator)對試件施加一個純彎矩激勵，結構動態響應利用貼在梁上的三個壓電應變片(PVDF)測量，記為傳感器S1～S3。試驗中設計了光滑連續荷載、三角形荷載，見圖12。為保證測得的樣本數據含有必需的結構響應信息，采樣頻率為2 500 Hz。兩種工況對應的測點傳感器的響應見圖12。

圖12 兩種工況下的實測結構響應和激勵Fig.12 Measured responses and excitations

利用移動時間窗進行荷載在線識別，每個窗內采用形函數方法識別。首先以測點S1為代表，對各工況的前2 000個點的實測響應進行傅里葉變化，根據式(6)計算形函數的參考頻率，從而確定每個窗內形函數的個數。兩種工況取lnt分別為30、40，定每個窗含的時間步數分別為420、400。移動時間窗前后重疊一半。兩種工況的總識別時間步分別為8 821、8 801。

這里激勵器施加的荷載是一個純彎矩，因此，至少需要一個傳感器才能準確識別荷載。分別利用單個測點S1，兩個測點S1－S2，三個測點S1－S3進行識別。為避免贅述，本文僅針對工況1給出了由不同測點數目識別的結果，見圖13。其中當沒有截斷奇異值時，由S1識別的荷載尾部發散，使荷載計算收斂的奇異值截斷數目為4;利用多個測點識別時不需奇異值截斷。

由圖13可以看出，利用單個傳感器可以識別荷載，但識別精度稍差。利用兩個傳感器和三個傳感器識別的結果基本上是一致的，識別精度更高。采用不同傳感器識別的荷載誤差曲線見圖14，考慮所有測量時間內的荷載識別相對誤差 F計算－F實測2/F實測2分別為51.83%，32.46%，26.6%。針對工況二只給出了利用兩個測點S1－S2識別的結果，見圖15，未采用奇異值截斷，荷載識別相對誤差為22.93%。

可以看出，兩種工況下利用提出的荷載形函數識別方法均能較好地識別未知荷載。識別的荷載歷程曲線與實際曲線吻合較好，能夠反映實際荷載信息。與多個傳感器相比，利用單個傳感器識別荷載，結構的動態響應信息相對較少，由于噪聲和測量誤差的影響，識別精度相對較低。移動時間窗在線識別荷載時，利用反卷積方法需要每個窗內重新計算L－曲線，確定奇異值截斷數目，這很耗時，難于實現。而荷載形函數在一定程度上能夠過濾噪聲的影響，不需要或只需少數幾個奇異值截斷，避免計算L－曲線。這進一步驗證在線荷載識別中形函數方法的重要角色。

5 結論

本文介紹并理論補充了基于形函數的荷載識別方法，通過數值算例和一個懸臂梁的試驗驗證了該方法的有效性，并得到以下主要結論:

(1)針對反卷積方法識別荷載中存在的測量時間長或采樣頻率高時數值求解困難的問題，利用有限的形函數逼近未知荷載，將識別離散的荷載歷程轉化為識別有限的荷載形函數系數，提高計算效率。

(2)荷載形函數是連續光滑的，在一定程度上能夠過濾噪聲的影響，提高逆問題對噪聲的魯棒性。

(3)通過實測結構響應的頻譜分析可以確定逼近未知荷載所需的形函數個數。

(4)形函數與移動時間窗相結合，每個窗內不需或只需截斷少數幾個奇異值就可以有效地識別未知荷載，從而實現荷載的在線識別。

［1] 瞿偉廉，王錦文.振動結構動態荷載識別綜述［J] .華中科技大學學報(城市科學版)，2004，21(4):94－99.

［2] Inoue H，Harrigan JJ，Reid SR.Review of inverse analysis for indirect measurement of impact force［J] .Applied Mechanics Reviews，2001，54(6):503 －524.

［3] Jacquelin E，Bennani A，Hamelin P.Force reconstruction:analysis and regularization of a deconvolution problem［J] .Journal of Sound and Vibration，2003，265(1):81－107.

［4] Uhl T.The inverse identification problem and its technical application［J] .Archive of Applied Mechanics， 2007，77(5):325－337.

［5] 智 浩，文祥榮，繆龍秀，等.動態載荷的頻域識別方法［J] .北方交通大學學報，2004(4):5－10.

［6] Inoue H，Ishida H，Kishimoto K.et al.Measurement of impact load by using analysis technique［C] .JSME Int.J.Series 1，1991，34(4):453 －458.

［7] Jankowski?.Off－line identification of dynamic loads［J] .Structural and Multidisciplinary Optimization，2009，37(6):609－623.

［8] Zhang Q，Jankowski?，Duan Z.Identification of coexistent load and damage based on virtual distortion method［C] .In:Proceedings of the 4th European Workshop on Structural Health Monitoring，Cracow，Poland，DEStech，2008:1121－1128.

［9] Doyle J F.A wavelet deconvolution method for impact force identification［J] .Experimental mechanics ，1997，37(4):403－408.

［10] 張 方，朱德懋，張福祥.動荷載識別的時間有限元模型理論及其應用［J] .振動與沖擊，1998，17(2):1－4.

［11] Ha QP，Trinh H.State and input simultaneous estimation for a class of nonlinear systems［J] .Automatica，2004，40:1779－1785.

［12] Allen M S，Carne T G.Delayed，multi－step inverse structural filter for robust force identification［J] .Mechanical Systems and Signal Processing，2008，22:1036－1054.

［13] Hansen P C.Numerical tools for analysis and solution of Fredholm integral equations of the first kind［J] .Inverse Problems 8，1992:849 －872.

［14] Zhang Q，Jankowski?，Duan Z.Identification of coexistent load and damage［J] .Structural and Multidisciplinary Optimization，2010，41(2):243－253.

［15] Dahlquist G，Bj?rck ?.Numerical methods in scientific computing，volⅡ［M] .2008，to be published by SIAM，URL http://www.mai.liu.se/～ akbjo/NMbook.html.