趙海霞　（南京航空航天大學金城學院，江蘇 南京 211156）

?

基于AR模型的空間脈動風速時程模擬方法研究

趙海霞（南京航空航天大學金城學院，江蘇 南京 211156）

空間脈動風速時程的模擬是橋梁抖振時域計算的前提。利用基于自回歸（AR）模型的線性回歸濾波器法，編制了模擬脈動風速的Matlab計算程序，模擬了一座大跨連續(xù)剛構(gòu)橋最大懸臂施工階段多點的水平和豎向脈動風速，并進行了功率譜和相關(guān)性檢驗，發(fā)現(xiàn)模擬效果較好。

脈動風速；AR模型；線性回歸濾波器法

1　研究背景

隨著我國橋梁工程的建設不斷向大跨度方向發(fā)展，橋梁的長細化使其剛度和阻尼不斷下降，導致結(jié)構(gòu)對風的敏感性不斷增加，對橋梁結(jié)構(gòu)進行風效應分析也逐漸成為橋梁結(jié)構(gòu)分析中的重要一環(huán)。由于時域分析方法可以直接計算頻域分析中難以處理的非線性響應，因而對大跨度橋梁的風載作用下的抖振行為進行非線性時域分析日益受到人們的關(guān)注。在抖振響應時域分析中，首先要根據(jù)目標功率譜函數(shù)人工模擬空間脈動風場，從而獲得離散的風速時程和離散的抖振力。

在人工模擬空間脈動風場方面，基于自回歸（Auto-Regressive，AR）模型的線性回歸濾波器法具有速度快、計算量小的特點，且精度較高，因此得到了廣泛應用[1][2]。該模型將均值為零的白噪聲隨機系列通過線性濾波器，使其輸出為具有指定譜特征的平穩(wěn)隨機過程[3]。本文首先詳細介紹了該方法的計算原理和過程，利用Matlab軟件編制了相應的計算程序；采用該方法模擬了一座大跨連續(xù)剛構(gòu)橋最大懸臂施工階段多點的水平和豎向脈動風速，并進行了功率譜和相關(guān)性檢驗。

2　自回歸（AR）模型

2.1求回歸系數(shù)

采用AR法推廣到模擬多維風速過程的技術(shù)，M個相關(guān)的隨機風過程[u（t）]=[u1（t），…，uM（t）]T可由下式生成：

式中：[u（t-kΔt）]=[u1（t-kΔt），…，uM（t-k Δt）]T；[N（t）]=[N1（t），…，NM（t）]T，Ni（t）為均值為0、具有給定協(xié)方差的正態(tài)分布隨機過程，i ＝1，…，M；[ψk]為M×M階矩陣，k=1，…，p；p 為AR模型的階數(shù)，一般取4或5。

對任一空間點i（i=1，…，M）具有時間差的隨機風過程ui（t）與ui（t-kΔt）的協(xié)方差可表示為：

由于ui（t）與ui（t-kΔt）為均值0的平穩(wěn)隨機風過程，其協(xié)方差的值僅為時間差的函數(shù)，式（2）可改寫為：

在式（1）同時右乘[u（t-jΔt）]=[u1（t-jΔt），…，uM（t-jΔt）]，j=1，…，p，并兩邊同時取數(shù)學期望（均值），考慮到[N（t）]的均值為0，且與隨機風過程ui（t）獨立，以及協(xié)方差Ru（jΔt）為偶函數(shù)，可得到協(xié)方差Ru（jΔt）與回歸系數(shù)ψk之間的關(guān)系，寫成矩陣形式，有：

其中

根據(jù)隨機振動理論[3]，功率譜密度與相關(guān)函數(shù)（協(xié)方差）之間符合維納-辛欽（Wiener-Khintchine）公式，即：

通過rjk（n）考慮風速時程的空間相關(guān)特性，其三維表達式為：

式中：Cx、Cy和Cz分別表示空間任意兩點左右、上下和前后的衰減系數(shù)，一般分別取為16、8和10[4]；（xi，yi，zi）、（xk，yk，zk）分別為空間i、k點的三維坐標，i、k=1，…，M；x、y、z分別為垂直于來流的水平方向、來流方向和豎向；分別表示第i點和第k點的平均風速。

求解式（4）給出的線性方程組，可以得到回歸系數(shù)矩陣[ψ]。

2.2求給定方差的隨機過程[N（t）]

對式（1）同時右乘[u（t）]=[u1（t），…，uM（t）]，并取期望，有

求出[RN]后，對其作喬利斯基（Cholesky）分解[RN]=[L][L]T，則

式中：[L]為下三角矩陣；[n（t）]=[n1（t），…，nM（t）]T，ni（t）是均值為0、方差為1且彼此相互獨立的正態(tài)隨機過程，i=1，…，M。

2.3求最終的M個隨機過程

求出回歸系數(shù)矩陣[ψ]及[RN]后，可按式（1）求解出M個空間相關(guān)的隨機風過程。

將式（1）按時間間隔Δt離散化，分別考慮三種情況，即ui（t）為①偶函數(shù)；②奇函數(shù)；③當t<0時，ui（t）＝0。由此可得出不同的矩陣方程形式。其中，以第三種假設計算起來最為方便，其遞推的矩陣表達式為：

從而得到M個具有時間、空間相關(guān)、時間間隔Δt的離散脈動風速時程向量。

為了避免模擬結(jié)果失真，Δt必須滿足以下條件：

式中：ωup是截斷頻率（rad/s）。

3　風場模擬的Matlab程序

根據(jù)線性回歸濾波器法的自回歸（AR）模型的思路，所編制的Matlab程序流程如圖1所示。

圖1　自回歸模型生成風速時程流程圖

4　算例

以某連續(xù)剛構(gòu)橋最大懸臂施工階段為對象，選定的模擬點即風荷載加載節(jié)點共23點，其中主梁上17點，墩上6點。圖2所示為部分模擬點位置示意圖。主梁上的模擬點模擬水平脈動風速和豎向脈動風速，墩上的模擬點僅模擬水平脈動風速。

圖2　風速模擬點示意圖（單位：m）

順風向風速模擬目標譜采用Kaimal水平脈動風速譜[5]：

豎直向風速模擬目標譜采用Panofsky譜：

式中：Su（n）、Sw（n）分別為水平向和豎向脈動風速功率譜；n為脈動風頻率（Hz）；Z為有效高度；U（Z）為高度Z處的平均風速；u*為氣流摩阻速度（亦稱剪切速度）；K是無量綱常數(shù)，K≈0.4；Zd為零平面高度；為周圍建筑物的平均高度；kd是地面阻力系數(shù)；Z0是地面粗糙長度。

考慮上述風速模擬的目標譜為非均方規(guī)一的單邊功率譜，通過與實際脈動風速的均方差比較，可以換算得到圓頻率表達的風譜形式：

模擬的其他參數(shù)有：場地類別為Ⅱ類；場地基本風速為24.1m/s[6]，其他高度處的風速按照抗風規(guī)范的指數(shù)分布，風速剖面無量綱冪指數(shù)α=0.16；地面粗糙長度Z0取為0.05；零平面高度Zd取為0；截至頻率取為2Hz；頻率等分數(shù)取為12000；采樣時距取為0.1s，共取12000步，1200s。分別選取左懸臂端（點1）、距左懸臂端8m處（點2）、主梁0號塊中點（點9）、右懸臂端處（點17）和墩左肢中部（點20）共五點作為展示對象。圖3是用自回歸（AR）模型模擬的這五點的水平風速時程，圖4是其豎向風速時程（不含點20）。

圖3　各點水平脈動風速

圖4　各點豎向脈動風速

圖5～圖6是分別用改進的平均周期圖法和最大熵法[5]對所模擬的點1風速時程做的譜估計；圖7所示為點1的自相關(guān)函數(shù)檢驗，圖8所示點1和點2的互相關(guān)函數(shù)的檢驗。從圖中我們可以看出，模擬值與目標值之間相差很小，可見模擬的結(jié)果良好。

圖5　點1功率譜密度函數(shù)（最大熵法）

圖6　點1功率譜密度函數(shù)（多周期圖法）

圖7　點1自相關(guān)函數(shù)

圖8　點1、2互相關(guān)函數(shù)

4　結(jié)語

本文詳細介紹了基于自回歸（Auto-Regressive，AR）模型的線性回歸濾波器法的計算原理和過程，利用Matlab軟件編制了相應的計算程序。采用該程序模擬了一座大跨連續(xù)剛構(gòu)橋最大懸臂施工階段多點的水平和豎向脈動風速，并進行了功率譜和相關(guān)性檢驗，發(fā)現(xiàn)模擬效果較好。因此本文的計算方法和程序具有一定的推廣價值。

[1]舒新玲，周岱.風速時程AR模型及其快速實現(xiàn)[J].空間結(jié)構(gòu)，2003，9 （4）：27-32.

[2]李元齊，董石麟.大跨度空間結(jié)構(gòu)風荷載模擬技術(shù)研究及程序編制[J].空間結(jié)構(gòu)，2001，7（3）：3-11.

[3]俞載道，曹國敖.隨機振動理論及其應用[M].上海:同濟大學出版社，1988.

[4]邊建烽，魏德敏.大跨空間結(jié)構(gòu)風速時程的數(shù)值模擬理論[J].暨南大學學報（自然科學版），2005，26（1）：87-90.

[5]JTG T D60-01-2004，公路橋梁抗風設計規(guī)范[S].

[6]張文明.大跨連續(xù)剛構(gòu)橋最大懸臂施工階段風致響應分析[D].武漢：華中科技大學，2007.

[7]黃文梅.信號分析與處理：MATLAB語言與應用[M].國防電子科技大學出版社，2000.

TU442.5+9

A

1007-7359(2016)03-0058-04

10.16330/j.cnki.1007-7359.2016.03.020

趙海霞（1981-），女，畢業(yè)于西安建筑科技大學，碩士；講師，南京航空航天大學金城學院土木工程教研室。