付康勝，呂中榮，汪利

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

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

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

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

1 基于Davenport譜的風(fēng)荷載

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

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

其中平均風(fēng)壓表示為

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

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

其中脈動(dòng)風(fēng)速的方差為

脈動(dòng)風(fēng)壓的方差為

因此，脈動(dòng)風(fēng)壓譜與脈動(dòng)風(fēng)速譜關(guān)系為

由脈動(dòng)風(fēng)壓譜通過諧波疊加法模擬脈動(dòng)風(fēng)壓時(shí)程，即

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

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

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

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

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

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

基于此，振動(dòng)方程（10）可以簡化為

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

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

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

其中矩陣AC和BC為

本文僅觀測(cè)模態(tài)坐標(biāo)下的加速度時(shí)程，觀測(cè)方程為

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

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

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

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

同樣，離散的觀測(cè)方程為

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

3 風(fēng)荷載反演

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

3.1 增廣卡爾曼濾波方法

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

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

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

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

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

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

其中增廣觀測(cè)矩陣為

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

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

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

增廣卡爾曼濾波由5個(gè)步驟組成：

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

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

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

3）計(jì)算卡爾曼增益矩陣

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

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

3.2 逐步最小二乘算法

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

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

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

其中

式中w 為權(quán)重矩陣，本文取為xi+1協(xié)方差矩陣的逆矩陣。當(dāng)系統(tǒng)自由度較大時(shí)，Am矩陣往往會(huì)具有較高的條件數(shù)，可以對(duì)式（41）進(jìn)行1-范數(shù)均衡預(yù)處理來降低待求逆矩陣的條件數(shù)［12］，利用矩陣指數(shù)的無窮積分代替其求逆過程，并且可通過精細(xì)積分法快速求解無窮積分過程［13］，該方法的精確性遠(yuǎn)高于利用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 層的風(fēng)荷載Davenport功率譜取為

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

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

圖2 風(fēng)荷載時(shí)程Fig.2 Wind load time history

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

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

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

圖4 增廣卡爾曼濾波反演風(fēng)荷載結(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

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

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

5 結(jié) 論

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

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