基于粒子濾波算法的風力發(fā)電塔地震動力響應預測

徐亞洲， 任倩倩， 于明陽， 時文浩

(西安建筑科技大學 土木工程學院，西安 710055)

隨著電子計算機技術的發(fā)展，有限元數(shù)值分析模型和結構動力試驗成為結構動力學中愈加成熟的領域，因此動力學模型更新和對動態(tài)系統(tǒng)結構性能的合理預測仍是非常重要的研究方向[1-2]。利用測量系統(tǒng)響應進行動力模型更新，在結構振動控制、結構健康監(jiān)測以及可靠性評估等方面都有廣泛的應用，但在建立結構計算模型及其激勵的數(shù)學模型的過程中，可能存在由于外部原因導致的建模誤差及參數(shù)不確定性影響，因此能夠準確和適當?shù)亓炕祟惒淮_定性對于結構動力模型的更新至關重要[3-5]。土木工程結構中的有限元模型往往與許多來源的不確定性和建模誤差有關，尤其是對于復雜的土木工程結構，在建模過程中有限元模型相比實際結構進行了大量的假設和簡化，當此類有限元模型用于預測一些感興趣的量時，這些不確定性因素會影響土木工程結構的設計和性能評估[6-8]。這些不確定性因素促使研究人員通過概率有限元模型更新方法來考慮潛在的結構不確定性。

基于經(jīng)典概率的理論分析模型是目前工程中常用的不確定性分析模型[9]，但在統(tǒng)計信息不充分的條件下，經(jīng)典概率分析往往具有一定的局限性。近年來，貝葉斯估計理論由于克服了經(jīng)典概率理論中難以處理的小樣本問題成為概率不確定性分析研究的新熱點[9]，貝葉斯的基本思路是根據(jù)狀態(tài)的先驗信息和測量信息來不斷遞推更新狀態(tài)的后驗概率密度函數(shù)[10]。毛文貴等[11]運用基于遺傳智能采樣技術的改進貝葉斯理論，可以很好地處理軸承轉子系統(tǒng)不平衡識別中的模型不確定性參數(shù)的影響。Pepi等[12]建立了位于特爾尼的彎曲斜拉橋有限元模型，該有限元計算模型考慮了模型不確定性以及對連接和構件的有限認識，在貝葉斯框架下進行環(huán)境振動試驗以確定用于模型更新的主要動態(tài)參數(shù)，并用于評價基于有限元模型計算結果的準確性。Sedehi等[13]在貝葉斯層次模型的基礎上，提出了漸近近似和最大后驗估計方法，該方法可以有效解決由外界因素引起的參數(shù)不確定性的影響。

在貝葉斯理論估計計算框架中，由于粒子濾波(particle fliter,PF)方法具有較好的處理復雜系統(tǒng)的能力而被廣泛研究及應用[14]。萬春風等[15-17]基于粒子濾波算法對單自由度和多自由度結構進行參數(shù)損傷識別研究，指出了粒子濾波算法在參數(shù)識別方面相較于其他濾波方法的優(yōu)越性。一些學者[18-21]采用粒子濾波方法在線預測疲勞裂紋擴展情況并更新剩余壽命，并通過試驗結果驗證了粒子濾波對疲勞裂紋擴展和剩余壽命在線預測的可行性和合理性。樊學平等[22-23]采用粒子濾波算法建立了隨時間更新的動態(tài)模型，并結合極值監(jiān)測數(shù)據(jù)實現(xiàn)橋梁結構可靠度指標的動態(tài)預測。

風力發(fā)電塔對動力荷載的高敏感性引起了研究者的關注[24-26]，且目前關于風力發(fā)電塔的抗震性能研究多采用有限元模擬的手段，但在構造結構數(shù)值分析模型的過程中，總是存在簡化模型所帶來的建模誤差和參數(shù)不確定性。振動臺試驗作為目前獲得結構動力響應分析的主要手段，被廣泛應用于結構性能評價與分析中，試驗結果對于校正能夠準確描述結構響應的有限元數(shù)值分析模型至關重要。因此本文提出了一個用于結構動力響應預測的概率貝葉斯估計計算框架，即首先在貝葉斯估計原理基礎上，建立地震激勵下結構有限元分析模型。然后采用粒子濾波方法，充分利用結構先驗信息并結合振動臺試驗觀測的風力發(fā)電塔動力響應數(shù)據(jù)信息選擇重要性函數(shù)，實現(xiàn)風力發(fā)電塔動力響應預測，并以有限元數(shù)值解與試驗值之間的差值來量化有限元模型的不確定性誤差。最后通過風力發(fā)電塔振動臺試驗結果，對計算方法的合理性及有效性進行驗證。

1 粒子濾波算法

在存在結構建模不確定性的情況下，需要確定并建立合理的結構模型，并基于實測結果進行加權修正。粒子濾波算法具有顯著處理參數(shù)不確定性估計的優(yōu)越性，目前被廣泛用于結構損傷識別和結構不確定性分析中。該算法主要是應用蒙特卡洛方法產(chǎn)生大量的隨機樣本(粒子)，然后結合觀測值利用樣本權重預測下一時刻的系統(tǒng)狀態(tài)或參數(shù)概率分布，最后根據(jù)算法不斷遞推更新得到結構參數(shù)的最優(yōu)貝葉斯估計值。

與其他濾波方法一樣，粒子濾波算法需要基于估計問題構建狀態(tài)空間模型，主要由狀態(tài)方程、觀測方程和初始狀態(tài)信息組成。狀態(tài)方程描述了系統(tǒng)狀態(tài)量隨時間的演變過程，觀測方程描述了系統(tǒng)適時輸出與系統(tǒng)當時狀態(tài)的模型，函數(shù)表達式如下

狀態(tài)方程為

Xk=f(Xk-1)+ωk

(1)

觀測方程為

Zk=h(Xk)+νk

(2)

初始狀態(tài)信息為

p(X0Z0)=p(X0)

(3)

式中:Xk和Xk-1分別為k時刻和k-1時刻的狀態(tài)值；Zk為k時刻的觀測值；f(·)為表達系統(tǒng)當前時刻與前一時刻的變化關系的已知非線性函數(shù)；h(·)為狀態(tài)值和觀測值之間的關系非線性函數(shù)；wk為過程噪聲，描述了系統(tǒng)的各種不確定性影響因素；vk為觀測噪聲，代表了觀測過程中的不確定性觀測誤差。

1.1 貝葉斯估計理論

貝葉斯的基本步驟包括預測和更新[27]，預測為利用模型前一時刻的先驗已知信息對當前時刻進行預測，更新是結合當前時刻的實際觀測進行遞推，從而獲得待估值在當前時刻的后驗概率密度，在貝葉斯更新過程中，狀態(tài)變量的更新是為了更好地表示真實的結構，其中更新過程是由先驗已知信息和所研究的結構測量中包含的測量信息完成的。狀態(tài)量中包含了所有的不確定性，因此，假設前一時刻的后驗概率p(Xk-1Z1:k-1)已知，則基于貝葉斯估計算法，可得到后驗概率密度p(XkZ1:k)的估計式

狀態(tài)預測方程

(4)

狀態(tài)更新方程

(5)

式(4)和式(5)組成了貝葉斯估計的計算框架。首先若已知k-1時刻的狀態(tài)值Xk-1，則可通過式(1)預測得到k時刻Xk的先驗概率密度p(Xk-1Z1:k-1)，基于觀測方程(式(2))可得到似然概率p(ZkX1:k)。在貝葉斯計算中，直接計算后驗概率密度函數(shù)牽涉到復雜的積分計算，對于非線性動態(tài)系統(tǒng)的狀態(tài)估計問題求解困難。因此，引入了基于隨機抽樣計算的蒙特卡洛方法來計算后驗概率。

1.2 蒙特卡洛原理

基于大數(shù)定理，蒙特卡洛方法[28]將所求解的隨機問題發(fā)生的概率用大量試驗中發(fā)生的頻率來估計，當統(tǒng)計數(shù)據(jù)比較大時，可以認為該隨機事件發(fā)生的頻率能無限接近期望值，從而可使貝葉斯后驗概率中的積分運算用大量樣本點的平均回報來逼近。

如對本研究來說，可以在狀態(tài)的后驗概率密度中大量的隨機采樣，則狀態(tài)的后驗概率密度可表達如式(6)。作為一種復雜問題有效的處理求解方式，該方法回避了工程應用中積分計算的困難，用狀態(tài)發(fā)生的頻率來表示概率，這在積分運算中有著非常重要的意義。

(6)

式中，δ(·)為狄拉克函數(shù)。然后待估狀態(tài)的期望值計算如式(7)所示

(7)

1.3 粒子濾波算法

(8)

作為本文研究問題所要求解的對象，故不能在待求解的后驗概率中采樣。一般采取從已知重要性密度函數(shù)q(XkZ1:k)中選取的N個樣本近似表示

(9)

其中,ω(Xk)為樣本權重，即

(10)

依據(jù)式(10)，式(9)可進一步化簡為

(11)

其中，歸一化后的權重值為

(12)

基于馬爾可夫假設與簡化，歸一化后粒子的權值(重要性權值)在每個時間步長上有如下遞推關系

(13)

通過權重值的不斷計算，這些權重被附加到每個相應的粒子上。但是這些粒子并沒有被修改，只是它們的相對權重被改變了，也就是說，當所有的粒子偏離觀測值時它們不會被拉回觀測值，重要性采樣只改變了它們的相對權重。這個缺點可能會導致單個粒子的權重占比非常大，而其他粒子的權值很小。因此，需要進行重采樣步驟，定義有效粒子數(shù)量為

(14)

當有效粒子的數(shù)量小于采樣數(shù)量N時，則表示粒子樣本退化嚴重，需要進行重采樣，即高重要性權值的粒子會被多次復制，低權值的粒子少量復制或被去除以保持粒子數(shù)量不變。重采樣后每個后驗粒子的權值相同，均為1/N。

重要性密度函數(shù)一般選取為狀態(tài)轉移密度函數(shù)，即

(15)

將式(15)代入式(13)，可重新得到權值更新方程的表達式

(16)

2 動力響應預測的粒子濾波實現(xiàn)

建立有限元模型進行風力發(fā)電塔動力響應分析的過程中，由于理想假設或簡化容易導致參數(shù)不確定性或模型不確定性，繼而導致有限元模擬值與實際測量值之間存在較大區(qū)別。因此本章基于第1章關于粒子濾波算法的計算原理的介紹，利用粒子濾波算法中的狀態(tài)空間模型將這些不確定信息考慮進來，作為先驗已知信息實現(xiàn)結構動力響應的預測過程。

2.1 力學模型

對于建筑結構，常用的力學模型分為層模型、桿模型及桿系-層間剪切模型。根據(jù)風力發(fā)電塔有限元分析方法的需要，目前有4種比較常見的建模方式。第一種為IEC規(guī)范[29]中所建議的單自由度模型簡化方法，即將塔體頂部葉輪、機艙的質量及塔筒50%的質量集中到塔筒頂部；第二種為我國GB 50011—2010《建筑抗震設計規(guī)范》[30]中常用的建議多自由度模型簡化方法，即將塔體頂部葉輪、機艙的質量集中到塔體頂部質量點處，各個筒段的質量分別集中到對應質量點處；第三種為簡化的有限元模型[31]，即葉輪、機艙采用剛體單元，將質點位置按照剛體運動耦合到塔筒頂部中心；第四種為精細化有限元模型[32]，即將葉片按照實際結構原尺寸建模從而考慮葉片與塔體相互作用對地震作用下結構性能的影響，每個葉片單元的質量集中在一個參考點，輪轂和機艙用集中質量點表示，通過剛性連接連接到塔頂部分。本節(jié)采用的第二種建模方法，即將風力發(fā)電塔簡化為多自由度層剪切結構進行分析。

以剪切型動力學模型為例，地震激勵下其動力運動方程為

(17)

結構的阻尼采用瑞麗阻尼，即

C=αM+βK

(18)

式中，α，β為阻尼系數(shù)。

2.2 動力方程求解

在結構動力方程的求解過程中，Newmark-β法應用較為廣泛，其主要思想是假定ti～ti+1時間段內的加速度變化呈線性規(guī)律，其假定為

(19)

式中，a，b為可調整參數(shù)。當a=1/2，b=1/6時，即為線性加速度法；當a=1/2，b=1/4時，即為平均常加速度法；當a=1/2，b=0時，即為中心差分法。本節(jié)采用線性加速度法計算。由式(19)可求得t+Δt時刻的速度和加速度

(20)

其增量表示形式為

(21)

將式(21)代入式(17)，可得到

(22)

通過求解式(22)可得到t+Δt時刻的位移，將其代入式(20)即可求得結構加速度及速度響應。

2.3 結構動力響應預測實現(xiàn)

基于試驗觀測數(shù)據(jù)，粒子濾波算法的結構動力響應預測過程的主要步驟如下。

步驟1首先需要構造狀態(tài)空間模型，根據(jù)Newmark-β方法求得有限元簡化模型結構的動力響應數(shù)值解作為趨勢項，用來構造本文的狀態(tài)方程；使用監(jiān)測儀器測得結構的實際動力響應，然后利用狀態(tài)量及振動臺試驗動力響應觀測量的信息構造觀測方程；初始狀態(tài)先驗信息基于有限元數(shù)值計算結果進行概率統(tǒng)計給出。

步驟2基于狀態(tài)空間模型，根據(jù)初始先驗信息的概率分布函數(shù)，利用蒙特卡洛抽樣方法可得到狀態(tài)變量的初始參數(shù)樣本，即初始先驗粒子{x0,i,i=1,2,3,…,N}。

步驟3在每一個時間步長上進行粒子更新，利用狀態(tài)方程可獲得下一個時刻的后驗粒子

(i=1,2,3,…,N)

步驟4基于k時刻的實際測量值Zk，可通過似然概率密度函數(shù)更新每個后驗粒子的權值，此權值衡量了該粒子產(chǎn)生的模擬數(shù)據(jù)接近實際觀測數(shù)據(jù)的近似程度。似然概率密度函數(shù)可近似為正態(tài)分布[33]，即

(23)

式中:Zmn,i為第i個粒子的模擬觀測數(shù)據(jù)；R為量測噪聲方差，可通過狀態(tài)量及觀測值的數(shù)據(jù)信息近似估計。

步驟5根據(jù)式(12)進行權重歸一化。

步驟6根據(jù)式(14)計算有效粒子的數(shù)量，當有效粒子的數(shù)量小于初始采樣數(shù)量N時，則進行重采樣，重采樣后每個粒子的權值相同，均為1/N。

以此不斷遞推更新，求得結構響應的最優(yōu)預測值。

具體流程圖如圖1所示。

圖1 粒子濾波算法預測過程流程圖

3 試驗驗證

本文開展了2 MW陸上風力發(fā)電塔振動臺試驗，基于狀態(tài)空間模型和試驗觀測值對有限元計算結果的不確定性進行修正，說明本文所提用于結構動力響應預測的粒子濾波框架的合理性。

3.1 試驗模型信息

根據(jù)量綱協(xié)調原理和結構動力方程[34-35]確定相似關系，按照相似原則及抗彎剛度等效設計模型。試件與原型結構的幾何相似比為1/20，模型材料為Q345鋼。塔呈錐形，直徑從底部的200 mm線性減小到頂部的130 mm。塔筒為4段，段與段之間用10 mm厚法蘭板連接，每段高度從下到上分別為1 m,0.85 m,1 m,1 m。由于薄板在軋制和焊接過程中制作困難，根據(jù)抗彎剛度等效原則，4個管段對應的管壁厚度分別為4 mm,3 mm,3 mm和3 mm。機艙和轉子的總質量為157.5 kg，4段塔筒的質量分別為85 kg，58 kg，51 kg，35 kg。試驗模型的示意圖如圖2所示。

按照抗震設計規(guī)范中輸入地震波的選取原則，從太平洋地震工程研究中心[36]數(shù)據(jù)庫選取了Chi-Chi波、Westmorland波、El-Centro波和Taft波作為模型結構振動臺輸入地震激勵，將垂直于葉片的方向定義為地震激勵輸入的方向，圖3給出了4條地震動記錄水平方向加速度時程。將地震動加速度峰值(peak ground acceleration，PGA)統(tǒng)一調幅至0.07g，其加速度反應譜如圖4所示。

(a) 試驗模型

(a)

(c)

地震波的輸入工況分別設置為0.07g，0.2g和0.4g。在不同幅值地震動輸入前后采用PGA為0.03g的白噪聲激勵來識別試驗模型的動力特性。本試驗加速度相似比為2，頻率相似比為6.324 6，故實際臺面輸入地震激勵需按照相似關系進行調整。

在風力發(fā)電塔結構抗震設計中，塔頂?shù)膭恿憫强拐鹦枨蟮闹匾笜?，因此通過合理布置測點測量了塔頂加速度響應及位移響應。不同地震動強度幅值輸入前后利用白噪聲掃頻，并對結果進行處理得到模型結構的動力特性。利用傳遞函數(shù)對試驗模型的固有頻率和振型進行識別，得到結構前兩階自振頻率分別為3.22 Hz,22.46 Hz，結構一階阻尼比為1.5%。

圖4 地震動記錄的加速度反應譜

3.2 模型動力響應預測

基于本文方法，對模型結構動力響應進行預測。

有限元數(shù)值分析：基于2.1節(jié)的介紹，本研究將風力發(fā)電塔動力系統(tǒng)簡化為5個自由度的剪切動力學模型來描述，如圖2(b)所示，根據(jù)質量點處的截面尺寸計算出截面抗彎剛度?；贛ATLAB有限元軟件，計算5個自由度簡化模型的質量和抗彎剛度矩陣。根據(jù)結構實測的阻尼比及自振頻率，采用瑞麗阻尼計算阻尼矩陣，最后利用2.2節(jié)中介紹的Newmark-β法求得有限元模型的塔頂加速度響應和位移響應。

狀態(tài)空間方程建立：根據(jù)12個地震工況作用下求得的5個自由度結構的動力響應數(shù)值計算結果，選用塔頂加速度響應和位移響應為狀態(tài)待估量，可估計出每個地震工況作用下的狀態(tài)噪聲方差；利用數(shù)值計算結果及對應工況下的振動臺試驗動力響應觀測量信息估計出測量噪聲方差；基于數(shù)值計算結果進行概率統(tǒng)計給出狀態(tài)量的初始信息，并基于初始狀態(tài)信息選取200個粒子。利用第2章中介紹的粒子濾波算法框架進行風力發(fā)電塔動力響應預測。

以振動臺試驗測得的觀測值為參考，基于粒子濾波算法對有限元模型結果進行修正預測。定義考慮和不考慮觀測值修正后的計算值與試驗實測值之間差值的絕對值作為動力響應的計算偏差。則考慮觀測值修正后(PF)，和不考慮觀測值修正(finite element method,FEM)的加速度和位移響應的偏差最大值計算結果，如圖5所示。結果表明有限元計算的動力響應偏差隨地震動PGA幅值的增大而增大。以Chi-Chi地震動激勵下塔頂加速度響應為例，當PGA值為140 gal時，加速度偏差為2.21 m/s2，而PGA值為800 gal時，加速度偏差為14.89 m/s2。考慮觀測值修正之后有限元計算結果偏差顯著減小，以塔頂位移響應為例，考慮觀測值修正之后位移偏差減小近一倍。

(a) 塔頂加速度響應偏差

(b) 塔頂位移響應偏差

為進一步說明本文方法的有效性及可行性，選取了El-Centro波和Chi-Chi波激勵作用下剪切動力學模型結構的加速度響應時程進行分析，圖6給出了粒子濾波算法修正前后加速度時程曲線對比。同一地震激勵下，隨輸入地震動幅值的增大，有限元計算結果相對于試驗值的偏差逐漸增大?？紤]觀測值修正后粒子濾波預測的結果精度較好，預測值與試驗實測值近似相等，具有很好的一致性，進一步驗證了本文所提校正不確定有限元分析模型的概率貝葉斯框架和粒子濾波預測算法的合理性，表明粒子濾波算法可以利用有限元分析結果不確定性的所有先驗信息，從而實現(xiàn)對包含變異性結構動力響應的穩(wěn)步預測。

4 結 論

本研究基于粒子濾波算法，提出了一個用于校正不確定有限元數(shù)值模型的概率貝葉斯計算框架對有限元計算的動力響應結果進行預測。利用某2 MW風力發(fā)電塔縮尺模型的振動臺試驗驗證了此框架的有效性和可行性。從分析結果可以得到以下結論：

(1) 有限元計算結果具有顯著的不確定性，且其計算誤差隨地震輸入幅值的增大而增大，考慮觀測值修正之后預測結果不確定性顯著減小。

(a)

(d)

(2) 粒子濾波模擬預測的動力響應值精度較高，模擬預測值與試驗實測值的時程曲線近似相等，變化規(guī)律具有很好的一致性。

(3) 粒子濾波算法與試驗觀測手段相結合，可以完成高效的結構動力響應預測，繼而為結構不確定性及可靠度分析提供參考。

結構不確定性分析應成為工程結構分析的一個重要部分，粒子濾波算法可以利用有限元分析結果不確定性的所有先驗信息，從而實現(xiàn)對結構動力響應的穩(wěn)步預測。