基于逐步最小二乘算法的風荷載識別*

付康勝，呂中榮，汪利

中山大學(xué)航空航天學(xué)院，廣東廣州510006

在我國經(jīng)濟發(fā)展較好的沿海城市，為了提高土地利用效率，超高層建筑往往被當作首選目標。而，在沿海地區(qū)臺風頻發(fā)，風荷載對建筑的影響甚是嚴重。超高層建筑柔性較高，低階固有頻率往往很低，同時風荷載也屬于低頻荷載，建筑在風荷載作用下很容易發(fā)生共振。玻璃幕墻在強風作用下容易發(fā)生損毀，高層建筑在受到繞流風作用時，建筑邊緣處易發(fā)生旋渦的脫落和氣流的再附，從而引起橫風向的振動，這對建筑的舒適性影響尤其嚴重［1］。因此，需要對作用于建筑上的風荷載進行測量，在設(shè)計時使其固有頻率盡可能地避開風荷載頻率?？紤]到建筑表面的復(fù)雜性，基于現(xiàn)有技術(shù)很難把荷載直接精確地測量出來。近年來有學(xué)者提出基于動力響應(yīng)反演風荷載，動力響應(yīng)數(shù)據(jù)的測量精度也遠高于荷載的測量精度。所以，利用實測動力響應(yīng)數(shù)據(jù)間接識別風荷載能獲得更加接近真實風荷載時程的數(shù)據(jù)。

國內(nèi)外的學(xué)者對結(jié)構(gòu)荷載反演進行了一些研究。譚棟等［2］利用新興群智能算法對作用于梁上的沖擊荷載進行識別，發(fā)現(xiàn)該方法抗噪性能較好。陳雋等［3］在結(jié)構(gòu)物理參數(shù)未知前提下，根據(jù)風荷載作用特點，利用荷載歸一化平均法，獲得作用在結(jié)構(gòu)上的風荷載時程。方明新［4］基于Kalman 濾波理論和超高層建筑的風振特性，研究了利用有限的風響應(yīng)數(shù)據(jù)反演超高層脈動風荷載的時程。Hwang 等［5］對高聳混凝土煙囪進行了風洞試驗，利用卡爾曼濾波方法和有限的動力響應(yīng)測量數(shù)據(jù)反演風荷載。Azam 等［6］提出一種對偶卡爾曼濾波方法識別作用在米蘭摩天大樓的沖擊荷載。Lourens 等［7］提出了一種增廣卡爾曼濾波方法識別作用在梁上的沖擊荷載，得到較好的識別結(jié)果。郅倫海等［8］基于遺忘因子最小二乘法和離散卡爾曼濾波方法識別高層建筑的脈動風荷載，并用風洞實驗驗證該方法的可靠性。葉靜嫻等［9］基于Newmark 方法和最小二乘法構(gòu)建目標函數(shù)，反演出作用在海洋平臺的波浪荷載。李正農(nóng)等［10］提出一種基于現(xiàn)場實測高層建筑加速度響應(yīng)反演脈動風荷載的方法，并將脈動風荷載與平均風荷載疊加，利用Newmark 方法進行時程分析，得到的加速度與實測加速度擬合得很好。

以上的工作主要基于卡爾曼濾波方法對荷載進行反演，需要假設(shè)荷載為狀態(tài)變量，給定其狀態(tài)演化誤差的協(xié)方差。不合理的協(xié)方差將得不到滿意的識別結(jié)果。本文提出一種逐步最小二乘算法反演風荷載，它基于精細積分法和最小二乘法構(gòu)建目標函數(shù)，逐個時間步識別風荷載。算例表明，該方法較增廣卡爾曼濾波方法精度高。

1 基于Davenport譜的風荷載

任意高度的風速通?？梢员硎鰹槠骄L速v 和脈動風速vf的和［11］，即

同樣地，對應(yīng)的風壓w 可以表示成平均風壓wˉ和脈動風壓wf的疊加，即

其中平均風壓表示為

式中w0為建筑物當?shù)鼗撅L壓；μs為結(jié)構(gòu)體型系數(shù)；μγ為重現(xiàn)期調(diào)整系數(shù)，模擬作用于高層建筑的風荷載時取1.1；βz為風壓高度變化系數(shù)。

脈動風為隨機過程，通常利用統(tǒng)計的方法進行描述。用于描述風荷載的風速譜有多種，我國現(xiàn)在規(guī)范使用的是Davenport 譜。Davenport 風速譜描述的風速功率譜Sv( f)不隨高度變化，其表達式為

其中脈動風速的方差為

脈動風壓的方差為

因此，脈動風壓譜與脈動風速譜關(guān)系為

由脈動風壓譜通過諧波疊加法模擬脈動風壓時程，即

φj為[0，2π]之間的一個隨機數(shù)，A 為風壓作用面積，Δf表示頻率步長。

2 精細積分法求解結(jié)構(gòu)動力響應(yīng)

對于一般的多自由度系統(tǒng)，其動力學(xué)方程為

在進行荷載反演之前，需要建立觀測值與荷載值的關(guān)系。本文采用精細積分法進行數(shù)值求解并建立觀測加速度與荷載的關(guān)系。實際中，在求解結(jié)構(gòu)動力響應(yīng)時，自由度數(shù)目往往會特別大，而真實響應(yīng)往往只與少數(shù)低階模態(tài)相關(guān)，因此可以采用振型疊加法對式（10）進行解耦和降維。工程應(yīng)用中，各節(jié)點的質(zhì)量和剛度難以準確確定，但系統(tǒng)前幾階固有頻率、阻尼比、振型可以較為容易地識別。比如，結(jié)構(gòu)的固有頻率可以通過對實際測量點的時域數(shù)據(jù)進行傅里葉變換，獲得對應(yīng)的自功率譜，再通過拾取峰值點來獲得；結(jié)構(gòu)振型則可以通過互功率譜與自功率譜的比值確定［10］。取前幾階振型，形成矩陣Φ，對應(yīng)的頻率矩陣為Ω2，根據(jù)振型的定義，有

位移yi+1與模態(tài)坐標下的位移ui+1滿足

根據(jù)質(zhì)量歸一化準則，假設(shè)阻尼為比例阻尼，有

基于此，振動方程（10）可以簡化為

其中對角矩陣Γ 的第i 個對角分量為2ξiωi，ξi表示第i階模態(tài)阻尼比，ωi表示第i階固有頻率，u(t)為模態(tài)響應(yīng)。

以模態(tài)坐標下的位移、速度構(gòu)建狀態(tài)變量

由振動方程（14）可得連續(xù)狀態(tài)方程

其中矩陣AC和BC為

本文僅觀測模態(tài)坐標下的加速度時程，觀測方程為

其中對應(yīng)的矩陣G和矩陣J為

由于狀態(tài)方程（16）為線性定常方程，根據(jù)精細積分法，其解可表達為

根據(jù)式（20）的結(jié)果，構(gòu)建離散狀態(tài)方程

其中狀態(tài)轉(zhuǎn)移矩陣A和B矩陣為

同樣，離散的觀測方程為

由式（22）和（23）可知，通過第i 時刻的位移、速度和已知的荷載時程可以求解出第i + 1 時刻的位移、速度和加速度。循環(huán)這個過程，可以獲得所有離散時間點位移、速度和加速度。

3 風荷載反演

荷載反演是典型的反問題，它通過觀測動力響應(yīng)反求荷載時程。本文以模態(tài)坐標下的加速度為測量數(shù)據(jù)，主要考慮兩種荷載反演方法-增廣卡爾曼濾波方法和逐步最小二乘算法。

3.1 增廣卡爾曼濾波方法

基于已推導(dǎo)出的離散狀態(tài)方程（21），引入滿足零均值高斯白噪聲假定的模型噪聲wi，則可以得到狀態(tài)方程

令相鄰時刻的荷載滿足以下關(guān)聯(lián)方程

ηi是一個隨機過程，滿足零均值高斯白噪聲假定。合并式（24）和式（25）的狀態(tài)向量，得到增廣的狀態(tài)向量

以及增廣狀態(tài)方程為

其中增廣系統(tǒng)矩陣Aa和增廣噪聲向量ξi為

在觀測方程中，取vi為第i 時刻的觀測噪聲，滿足零均值高斯白噪聲假定，則有增廣觀測方程

其中增廣觀測矩陣為

卡爾曼濾波在最小方差無偏估計下定義為最優(yōu)的線性狀態(tài)估計器，上文給出的離散狀態(tài)方程（30）和離散觀測方程（29），涉及到的噪聲之間互不相關(guān)。用矩陣Q、R 和S 分別表示模型噪聲方差強度矩陣、觀測噪聲方差強度矩陣和荷載方差強度矩陣，則應(yīng)滿足以下關(guān)系

其中δi-j是克羅內(nèi)克函數(shù)。

增廣狀態(tài)變量的估計誤差協(xié)方差矩陣Pi|j定義為

增廣卡爾曼濾波由5個步驟組成：

1）用第i 時刻的狀態(tài)向量預(yù)測第i + 1 時刻的狀態(tài)向量

2）用第i 時刻的估計誤差協(xié)方差矩陣Pi|i預(yù)測第i + 1時刻的估計誤差協(xié)方差矩陣

其中增廣狀態(tài)向量的方差強度陣為

3）計算卡爾曼增益矩陣

4）利用預(yù)測的第i + 1 時刻的狀態(tài)向量和第i時刻的觀測加速度修正對狀態(tài)向量的預(yù)測，即

5）更新估計誤差協(xié)方差矩陣

3.2 逐步最小二乘算法

逐步最小二乘算法針對每一個時間步建立目標函數(shù)。首先，對于兩個相鄰時刻，已知第i 時刻的位移、速度和荷載和第i + 1時刻的觀測加速度，反演第i + 1 時刻的位移、速度和荷載。構(gòu)建的目標函數(shù)

式中λ 為權(quán)值，其作用是平衡模型誤差和觀測誤差。

對式（40）的xi+1和fi+1求變分，獲得線性方程組

其中

式中w 為權(quán)重矩陣，本文取為xi+1協(xié)方差矩陣的逆矩陣。當系統(tǒng)自由度較大時，Am矩陣往往會具有較高的條件數(shù)，可以對式（41）進行1-范數(shù)均衡預(yù)處理來降低待求逆矩陣的條件數(shù)［12］，利用矩陣指數(shù)的無窮積分代替其求逆過程，并且可通過精細積分法快速求解無窮積分過程［13］，該方法的精確性遠高于利用matlab直接求逆。

4 數(shù)值算例

考慮圖1 所示的16 自由度的剪切層模型，取各樓層具有相同的質(zhì)量m1= m2= m3= …= m16=3 500 kg 和相同的剛度k1= k2= k3= …= k16=1500 kN/m。

模型的各階模態(tài)阻尼比均取為2%。荷載的空間分布矩陣為

地面粗糙度指數(shù)α 取0.23。第10 層的風荷載Davenport功率譜取為

模擬風荷載的頻率范圍取為ω ∈[0，15]，參考利用諧波疊加法（9） 構(gòu)建風荷載時程，具體見圖2。

圖1 高層建筑簡化而來的剪切層模型Fig.1 Simplified shear layer model for high-rise buildings

圖2 風荷載時程Fig.2 Wind load time history

基于精細積分法求解加速度時程時，取時間長度為5 s，時間間隔為0.01 s，初始位移和速度均取為0?？紤]到測量加速度數(shù)據(jù)時存在誤差，給仿真的加速度添加一定水平的高斯白噪聲

圖4為基于增廣卡爾曼濾波方法反演的風荷載結(jié)果，其荷載協(xié)方差矩陣S 根據(jù)L 曲線方法已經(jīng)取到最優(yōu)（見圖5）；由于模型精度很高，其模型噪聲協(xié)方差矩陣Q 取為I × 10-20；一般來說，觀測噪聲協(xié)方差矩陣取決于傳感器精度。

圖3 模態(tài)加速度時程Fig.3 Time history of modal acceleration

圖4 增廣卡爾曼濾波反演風荷載結(jié)果（5%噪聲水平）Fig.4 The wind load is retrieved by the augmented Kalman filter(5%noise level)

圖5 L曲線法確定最優(yōu)荷載協(xié)方差矩陣SFig.5 L curve method is used to determine the optimal load covariance matrix S

簡化起見，本算例的觀測噪聲協(xié)方差矩陣直接由真實加速度和污染加速度計算出來。該方法的識別結(jié)果相對誤差為10.07%，本文相對誤差δ計算公式為

圖6為基于最小二乘法反演的風荷載結(jié)果，可以看出，僅考慮前四階的模態(tài)加速度數(shù)據(jù)，在5%噪聲干擾下，識別的風荷載很好地逼近真實荷載，識別結(jié)果相對誤差為5.66%。說明該方法抗噪性能良好。可以看出，本文所提的逐步最小二乘算法比增廣卡爾曼濾波方法具有更好的識別精度和抗噪能力。

5 結(jié) 論

本文提出了一種逐步最小二乘算法反演風荷載。該方法以模態(tài)加速度為觀測數(shù)據(jù)，針對每一個時間步建立目標函數(shù)，逐步識別荷載。與已有的增廣卡爾曼濾波方法相比，所提方法不需要荷載的協(xié)方差信息。算例表明，兩種方法在5%噪聲水平下，仍然獲得不錯的識別結(jié)果，說明兩種方法的抗噪性能都較好，但逐步最小二乘算法的識別精度優(yōu)于增廣卡爾曼濾波方法。

圖6 逐步最小二乘法反演風荷載時程結(jié)果（5%噪聲水平）Fig.6 Stepwise least square method is used to invert wind load time history(5%noise level)