基于譜表示與AR模型的風場模擬對比分析

劉瑞莉， 徐公勇， 宋淳宸， 邱順冬

(1. 西南交通大學，四川成都 610031； 2. 成都興城人居地產(chǎn)投資集團有限公司， 四川成都 610031； 3.建研科技股份有限公司，北京 100000)

隨著計算機技術(shù)的快速發(fā)展，蒙特卡洛模擬在工程實際上得到越來越廣泛的應(yīng)用。與頻域方法相比，蒙特卡洛模擬方法在考慮結(jié)構(gòu)非線性、系統(tǒng)隨機性以及其它相關(guān)問題等方面具有明顯的優(yōu)勢。例如，其被廣泛應(yīng)用于考慮結(jié)構(gòu)和氣動力非線性的大跨橋梁抖振響應(yīng)分析[1]。此外，蒙特卡洛模擬方法也常被用作檢驗頻域方法計算精度的標準。

目前，蒙特卡洛模擬方法主要有譜表示方法和線性濾波法。譜表示方法由于具有嚴密的理論推導以及容易使用等特點而得到廣泛的應(yīng)用，而快速傅里葉變換(FFT)算法的引入使得該類方法的模擬效率得到顯著提高[2-4]。最初，多點平穩(wěn)隨機過程的譜表示方法為單指標方法，模擬出的樣本不具備各態(tài)歷經(jīng)性[5]。隨后，Deodatis[6]將其拓展到了多指標的譜表示方法，以實現(xiàn)各態(tài)歷經(jīng)性的要求。為了進一步提高模擬效率，研究者從減少功率譜矩陣分解時間方面提出了多種改進的譜表示方法[7-9]。近期，Peng等[10]提出了一種基于隨機場的風場模擬方法，該方法不僅可避免功率譜矩陣分解，同時可利用二維FFT極大提高模擬效率。

線性濾波法又稱白噪聲濾波法，即將均值為零的高斯白噪聲通過設(shè)計的濾波器，從而獲得具有目標統(tǒng)計特性的隨機過程。該類方法包含自回歸(Autoregressive, AR)模型、移動平均 (Moving averages, MA)模型以及自回歸和移動平均模型(Autoregressive and moving-average, ARMA)模型[11-13]。與ARMA模型相比，AR模型的參數(shù)相對較少，因此更加便于使用，其已被廣泛應(yīng)用于多點風速時程的模擬。Novak等[14]對基于AR模型的風速模擬誤差進行了詳細分析。舒新玲等[15]采用AR模型模擬了大跨屋蓋的隨機脈動風速時程，并提出VC與Matlab混合編程模擬的快速實現(xiàn)方法。李春祥等[16]基于AR模型模擬了一幢高度為200 m超高層建筑的風速時程，并驗證了其可行性。李明等[17]基于西堠門大橋風場模擬的算例，對AR模型的定階進行了討論，并認為階數(shù)過高與過低均將影響模擬精度，選擇合適的模型階數(shù)即可得到滿意的精度又可節(jié)省計算時間。

基于上述研究基礎(chǔ)，本文對譜表示方法與基于AR模型的風場模擬方法進行對比研究。首先，分別對譜表示方法和基于AR模型的模擬方法進行了介紹。隨后，以某橋梁風場模擬為算例，對兩種方法在模擬精度和效率方面進行了對比。最后，給出了相應(yīng)的結(jié)論。

1　經(jīng)典譜表示方法

考慮一個均值為零的n變量平穩(wěn)隨機過程X(t)=[X1(t),X2(t),...,Xn(t)]T，設(shè)其功率譜密度矩陣可表示如下：

(1)

式中：ω為頻率；Sjk(ω),j,k=1, 2, …,n;j≠k為互功率譜密度函數(shù)，一般為復(fù)數(shù)；當j=k時，Sjj(ω)為自功率譜密度函數(shù)。功率譜矩陣元素具有如下性質(zhì)：

Sjj(ω)=Sjj(-ω),j=1, 2,…,n

(2)

j,k=1, 2, …,n;j≠k

(3)

其中*表示復(fù)數(shù)共軛。由于功率譜矩陣S(ω)為非負定的Hermitian矩陣，因此其可進一步分解如下：

S(ω)=H(ω)HT*(ω)

(4)

式中：H(ω)為下三角矩陣，其元素Hjk(ω)可表示為如下的極坐標形式：

Hjk(ω)=|Hjk(ω)|eiθjk(ω)

(5)

(6)

式中：Im和Re分別表示復(fù)數(shù)的虛部和實部；上述多點隨機過程被分解為n個完全不相干的子多點隨機過程的疊加，且每個子過程的元素之間完全相干。因此，多點隨機過程X(t)可采用下式模擬[6]：

j=1, 2,…,n

(7)

式中：Δω=ωu/N；ωu為截止頻率；N為頻率離散點數(shù)；ωml=(l-m/n)Δω,m=1, 2, ...,n；φml為在[0,2π]上服從均勻分布的隨機相位角。值得說明，當ωml退化為ωl=lΔω時，上述經(jīng)典的雙指標模擬公式即退化為單指標模擬公式。

2　基于AR模型的模擬方法

設(shè)上節(jié)的n變量平穩(wěn)隨機過程X(t)可表示為如下的多點AR模型形式[9]：

(8)

式中p為模型階數(shù)；Φr為n×n的自回歸系數(shù)矩陣；L為n×n的系數(shù)矩陣；Y(t)=[Y1(t),Y2(t),...,Yn(t)]T是均值為0的高斯白噪聲隨機過程，且滿足下式：

E[Y(t)YT(s)]=Inδ(s-t)

(9)

式中δ(·)為狄拉克δ函數(shù)；In為n×n的單位矩陣。

將式兩邊同乘于X(t-qΔt)，然后取數(shù)學期望可得：

q=1,2,…,p

(10)

式中RXX(-qΔt)=E[X(t)XT(t-qΔt)]。上式可進一步表示為如下的矩陣形式：

(11)

(12)

由式(11)，可求得不同時刻的自回歸系數(shù)矩陣：

(13)

將式(8)兩邊同乘于YT(t)和XT(t)，并求期望后可分別得到：

RXY(0)=L

(14)

(15)

將式(14)代入式(15)，可得如下的表達式：

(16)

設(shè)L為下三角矩陣，其可以通過對等式右邊部分進行Cholesky分解得到。在求得Φr和L矩陣后，則可以通過該AR模型進行多點風速隨機過程的模擬。

3　數(shù)值算例

為了對比譜表示方法和基于AR模型的模擬方法在模擬精度和效率上的差異，對泰州長江大橋主梁上的順風向脈動風速進行了模擬。該橋為連接泰州市和揚中市的三塔懸索橋，其跨徑布置為390 m、1 080 m、1 080 m及390 m，如圖1所示。設(shè)模擬點的間距為7.2 m，則共需要模擬301個點。

圖1　泰州長江大橋示意

設(shè)主梁模擬點的自功率譜相等，且均用如下的雙邊Kaimal譜描述:

(17)

此外，采用如下的Davenport相干函數(shù)描述不同模擬點之間的相干性：

(18)

式中：λ=7為衰減系數(shù)；ξ表示兩點之間的距離。

為了對比，譜表示方法考慮雙指標和單指標兩種工況，基于多點AR模型的模擬方法選取較為常用的4階和6階模型兩種工況[16]。

3.1　模擬精度對比

在精度對比時，截止頻率均取為4π。譜表示方法的頻率離散點數(shù)為1 024，時間步長為0.25 s。在基于AR模型的模擬方法中，頻率離散點數(shù)為2 048，時間步長為0.2 s。限于篇幅，本文僅給出基于4階AR模型模擬方法得出的風速時程樣本，如圖2所示?？梢钥闯?，模擬點1與2之間的風速時程接近，而模擬點76則與前兩者差別較大，其原因是模擬點76與點1之間的距離已超500 m，相關(guān)性較小。

圖3給出了不同計算工況下模擬時程的估計功率譜與理論功率譜的對比?？梢钥闯觯p指標和單指標譜表示方法的模擬精度均高于AR模型的模擬精度。此外，雙指標譜表示方法的精度高于單指標譜表示方法，4階AR模型的模擬精度則高于6階AR模型。圖4給出了不同工況下的自/互相關(guān)函數(shù)對比圖?？梢钥闯?，雙指標譜表示方法的精度同樣明顯高于其它模擬方法。而單指標譜表示方法在相關(guān)函數(shù)的模擬精度上與基于AR模型的模擬方法差別不大。

3.2　模擬效率對比

在進行模擬效率對比時，所有工況的截止頻率和頻率離散點數(shù)均分別取為4π和1 024，模擬點數(shù)分別取32、64、128、256以及512五種工況。本算例的所有工況均在一臺處理器為Intel (R) Core (TM) i7-7500U，主頻為2.70GHz，內(nèi)存為8G的筆記本電腦上進行，計算結(jié)果如表1所示?？梢钥闯觯斈M點數(shù)為64時，基于4階AR模型的模擬方法效率為單指標譜表示方法的2倍，而為雙指標譜表示方法的80倍。因此，基于AR模型的模擬方法效率明顯高于譜表示方法，尤其是與雙指標譜表示方法相比。此外，4階和6階AR模型的模擬效率相差較小。在所需內(nèi)存上，當模擬點數(shù)分別為512和128時，單指標和雙指標譜表示方法即出現(xiàn)內(nèi)存不足的現(xiàn)象。由此可見，譜表示方法所需要的內(nèi)存明顯大于基于AR模型的模擬方法。

(a)模擬點1

(b)模擬點2

(c)模擬點76

(a) 單指標譜表示方法

(b) 雙指標譜表示方法

(c) 4階AR模型

(d) 6階AR模型

4　結(jié)論

本文首先對譜表示方法與基于AR模型的風場模擬方法分別進行了介紹，然后通過數(shù)值算例詳細對比了單指標譜表示方法、多指標譜表示方法、基于4階AR模型以及6階AR模型的風場模擬方法在模擬精度和效率上的差異，得出了以下結(jié)論：

(1)雙指標譜表示模擬方法的模擬精度明顯高于單指標譜表示方法以及基于AR模型的模擬方法。

(a)模擬點1

(b) 模擬點2

(c)模擬點1與2

(d) 模擬點1與76

模擬點數(shù)單指標譜表示雙指標譜表示4階AR模型6階AR模型 322.4924.020.881.07646.06237.812.953.2112822.41/11.9616.08256112.82/68.1080.49512//458.46491.73

備注：表中/表示電腦內(nèi)存不足。

(2)單指標譜表示方法在功率譜的模擬精度上高于基于AR模型的模擬方法，而在相關(guān)函數(shù)模擬精度上，兩者之間差別不明顯。

(3)基于AR模型的模擬方法在模擬效率以及消耗內(nèi)存上的性能均明顯優(yōu)于譜表示方法。因此，當風場模擬點數(shù)很大時，可考慮使用基于AR模型的模擬方法。