周彬彬，蔡建國(guó)，馮 健

（東南大學(xué)土木工程學(xué)院，南京 210096）

基于ARMA模型和時(shí)空K riging插值聯(lián)合模擬大跨結(jié)構(gòu)的脈動(dòng)風(fēng)速時(shí)程

周彬彬，蔡建國(guó)，馮 健

（東南大學(xué)土木工程學(xué)院，南京 210096）

采用ARMA模型模擬具有時(shí)間和空間相關(guān)性的大跨結(jié)構(gòu)的隨機(jī)風(fēng)速時(shí)程。并對(duì)Whittle遞推算法進(jìn)行了改進(jìn)，求解出自回歸系數(shù)和滑動(dòng)回歸系數(shù)矩陣。利用MATLAB軟件編制了模擬大跨橋梁的風(fēng)速時(shí)程模擬程序，運(yùn)行結(jié)果表明：改進(jìn)的方法計(jì)算效率高，模擬精度符合要求。然后又利用時(shí)空Kriging插值法對(duì)控制點(diǎn)區(qū)域內(nèi)節(jié)點(diǎn)的脈動(dòng)風(fēng)速時(shí)程進(jìn)行模擬，其模擬精度較好且求解過(guò)程快，可用于大跨空間結(jié)構(gòu)、超高層建筑等含有大量節(jié)點(diǎn)的結(jié)構(gòu)風(fēng)場(chǎng)模擬。

數(shù)值模擬；脈動(dòng)風(fēng)速時(shí)程；ARMA模型；Whittle遞推算法；時(shí)空Kriging插值

隨著計(jì)算機(jī)科學(xué)技術(shù)的發(fā)展，隨機(jī)風(fēng)場(chǎng)的數(shù)值模擬技術(shù)得到了很大的提高。模擬脈動(dòng)風(fēng)速時(shí)程主要有諧波合成法和基于線(xiàn)性濾波技術(shù)的線(xiàn)性濾波器法（自回歸線(xiàn)性濾波器AR法和自回歸滑動(dòng)平均線(xiàn)性濾波器ARMA法）兩大類(lèi)。其中ARMA法能夠在保證模擬精度的前提下，高效地模擬出具有時(shí)間和空間相關(guān)性的隨機(jī)風(fēng)速時(shí)程［1］。而時(shí)空Kriging插值法是一種近年來(lái)才發(fā)展起來(lái)的基于隨機(jī)變量在空間和時(shí)間上的相關(guān)結(jié)構(gòu)而建立起來(lái)的統(tǒng)計(jì)方法［2］，它以建立有效的時(shí)空協(xié)方差（變異）函數(shù)為前提，根據(jù)變量的時(shí)空自相關(guān)性對(duì)待估點(diǎn)進(jìn)行插值［3］。

本文在ARMA模型用于模擬多變量穩(wěn)態(tài)隨機(jī)過(guò)程的基礎(chǔ)上進(jìn)一步簡(jiǎn)化求解算法并和時(shí)空Kriging插值法聯(lián)合用于對(duì)大跨結(jié)構(gòu)的脈動(dòng)風(fēng)速時(shí)程進(jìn)行模擬。隨后利用MATLAB軟件編制風(fēng)速模擬程序，運(yùn)行效果表明該聯(lián)合方法在保證精度的前提下，能顯著地提高大跨結(jié)構(gòu)風(fēng)速時(shí)程的模擬效率。

1 風(fēng)的基本特性

近地風(fēng)一般包含平均風(fēng)和脈動(dòng)風(fēng)兩部分。平均風(fēng)在給定的時(shí)間間隔內(nèi)風(fēng)力大小、方向等不隨時(shí)間而改變，而脈動(dòng)風(fēng)則隨時(shí)間和空間隨機(jī)地變化，在數(shù)學(xué)上屬于隨機(jī)過(guò)程范疇，工程中常常作為具有零均值的各態(tài)歷經(jīng)的高斯隨機(jī)過(guò)程來(lái)處理。風(fēng)場(chǎng)模擬主要是針對(duì)脈動(dòng)風(fēng)而言［4］。作用于結(jié)構(gòu)上任一點(diǎn)的風(fēng)速V（x，y，z，t），可以表達(dá)為平均風(fēng)速和脈動(dòng)風(fēng)速v（x，y，z，t）之和，即：

2 ARMA模型解法

常規(guī)算法中，使用如Gauss消元法、Gauss-Seidel迭代法對(duì)線(xiàn)性方程組（18）進(jìn)行求解，得到Φi，Ψj后，代入式（5）即可進(jìn)行風(fēng)速模擬。但對(duì)于大跨橋梁、大跨空間結(jié)構(gòu)、超高層建筑這樣含有大量節(jié)點(diǎn)的結(jié)構(gòu)求解過(guò)程是繁瑣的，精度也很難控制。本文使用Whittle遞推算法用于求解自回歸系數(shù)Φi［7］進(jìn)而推導(dǎo)出了滑動(dòng)回歸系數(shù)Ψj。算法如下：

式中：Φi，j為p階AR模型在階次為i時(shí)的第j個(gè)系數(shù)，為m×m階矩陣，i，j＝1，2，…p利用得到的自回歸系數(shù)Φi，推出滑動(dòng)回歸系數(shù)Ψj：

至此，ARMA模型的所有系數(shù)矩陣均已求出。

3 時(shí)空K riging插值

3.1 時(shí)空協(xié)方差函數(shù)

設(shè)U（s，t）∈Rk×T是定義在時(shí)空域的隨機(jī)過(guò)程，其中Rk代表k維的歐氏空間，T代表時(shí)間，同時(shí)定義（si，ti）為時(shí)空?qǐng)鲋腥我鈽颖军c(diǎn)的位置，Δs為樣本點(diǎn)間的空間距離，Δt為時(shí)間距離。當(dāng)U（s，t）滿(mǎn)足二階平穩(wěn)

或固有假設(shè)時(shí)，其協(xié)方差函數(shù)記為［8］：

對(duì)應(yīng)的變異函數(shù)記為：

式中σ2為U（s，t）的方差。

本文采用積和模型來(lái)擬合風(fēng)速的時(shí)空變異結(jié)構(gòu)如下：

式中：Cst為時(shí)空協(xié)方差；Cs為空間協(xié)方差；Ct為時(shí)間協(xié)方差；γst，γs，γt分別為對(duì)應(yīng)變異函數(shù)；Cst（0，0），Cs（0），Ct（0）分別是對(duì)應(yīng)的基臺(tái)值［2］。

由變異函數(shù)性質(zhì)得：

實(shí)驗(yàn)中往往采用式（26）計(jì)算樣本變異函數(shù)：

式中：U（si，ti）為樣本觀測(cè)值；N（Δs，Δt）為相距（Δs，Δt）的樣本對(duì)數(shù)［9］。

當(dāng)Δt＝0時(shí)，式（27）為純空間域的樣本變異函數(shù)估計(jì)值：

當(dāng)Δs＝0時(shí)，式（28）為純空間域的樣本變異函數(shù)估計(jì)值：

3.2 K riging插值

Kriging插值是一種線(xiàn)性的無(wú)偏估計(jì)方法，要求估計(jì)誤差的方差最?。?0］。本文采用普通Kriging進(jìn)行時(shí)空擴(kuò)展，即：

U*（s0，t0）為時(shí)空點(diǎn)處（s0，t0）的估計(jì)值，λi為鄰近觀測(cè)值的加權(quán)系數(shù)［2］。引入拉格朗日系數(shù)μ可得參數(shù)λi的計(jì)算式為：

將時(shí)空Kriging插值與ARMA模型聯(lián)合運(yùn)用于大跨橋梁、大跨空間結(jié)構(gòu)、超高層建筑的風(fēng)速時(shí)程的模擬，即先使用ARMA模型對(duì)少量的結(jié)構(gòu)控制節(jié)點(diǎn)處進(jìn)行風(fēng)速時(shí)程模擬，再利用時(shí)空Kriging插值對(duì)控制節(jié)點(diǎn)處臨近的節(jié)點(diǎn)進(jìn)行風(fēng)速時(shí)程的模擬。

4 工程算例

某橋?yàn)橹骺缈鐝竭_(dá)1 428 m的三跨連續(xù)鋼箱梁懸索橋。本文對(duì)橋主梁處（圖1）的脈動(dòng)風(fēng)速時(shí)程進(jìn)行模擬。數(shù)值模擬中采用Kaimal風(fēng)速譜，10 m處平均風(fēng)速為31.2m／s，地面粗糙度α＝0.142；模擬計(jì)算步長(zhǎng)Nc＝10 000，步長(zhǎng)Δt＝0.5 s，模擬時(shí)程為T(mén)0＝5 000 s，頻率取值范圍為0～2 Hz。橋面縱向相鄰位置等間距50m，共模擬29個(gè)點(diǎn)。選取部分節(jié)點(diǎn)的風(fēng)速模擬曲線(xiàn)見(jiàn)圖2及模擬譜與目標(biāo)譜的對(duì)比見(jiàn)圖3?？梢钥闯瞿M譜與目標(biāo)譜吻合度高，模擬效果較好，計(jì)算耗時(shí)為0.2 s。本文又采用傳統(tǒng)的Gauss消元法進(jìn)行計(jì)算對(duì)比，其耗時(shí)約為0.47 s。使用改進(jìn)的算法要比傳統(tǒng)算法效率提高57%。

圖1 主梁風(fēng)速模擬點(diǎn)位示意（單位：cm）Fig.1 Positions for simulated points of wind speed on main girders（unit：cm）

圖2 14、15、16節(jié)點(diǎn)風(fēng)速時(shí)程Fig.2 Simulated wind speed time-history series at points 14，15，16

圖3 14、15、16節(jié)點(diǎn)模擬風(fēng)速譜與目標(biāo)譜對(duì)比Fig.3 Comparison between target power spectrum and simulated power spectrum at points 14，15，16

利用已求解出控制點(diǎn)的風(fēng)速時(shí)程數(shù)據(jù)，運(yùn)用式（27）擬合出的空間變異函數(shù)表達(dá)式為式（31），運(yùn)用式（28）擬合出的時(shí)間變異函數(shù)，表達(dá)式為（32）：

γs（Δs）＝30.85＋

求出k1，k2，k3值后將式（31）和式（32）代入式（24）得時(shí)空變異模型如圖5。由求出的時(shí)空變異函數(shù)，根據(jù)測(cè)點(diǎn)14、16處與測(cè)點(diǎn)15處的時(shí)空關(guān)系，利用測(cè)點(diǎn)14、16處的風(fēng)速時(shí)程，對(duì)測(cè)點(diǎn)15處的風(fēng)速進(jìn)行時(shí)空Kriging插值，將所得的結(jié)果與已模擬的風(fēng)速時(shí)程及功率譜進(jìn)行對(duì)比如圖6和圖7。從結(jié)果來(lái)看，時(shí)空Kriging插值的結(jié)果與已模擬的結(jié)果重合度較高，且整個(gè)插值過(guò)程耗時(shí)少。

圖4 空間和時(shí)間的變異函數(shù)Fig.4 Spatial variogram and temporal variogram

圖5 時(shí)空變異函數(shù)模型Fig.5 Spatial-temporal variogram model

圖6 15節(jié)點(diǎn)插值風(fēng)速與模擬風(fēng)速對(duì)比Fig.6 Comparison between interpolated wind speed time-history and simulated wind speed time-history at point 15

圖7 15節(jié)點(diǎn)插值風(fēng)速譜與模擬風(fēng)速譜對(duì)比Fig.7 Comparison between interpolated power spectrum and simulated power spectrum at point15

5 結(jié) 論

本文推導(dǎo)出自回歸（AR）階數(shù)p大于滑動(dòng)回歸（MA）階數(shù)q時(shí)，ARMA模型用于模擬大跨結(jié)構(gòu)風(fēng)速時(shí)程的公式。并針對(duì)求解矩陣的特點(diǎn)，對(duì)其解法進(jìn)行了簡(jiǎn)化。采用改進(jìn)后的Whittle遞推算法使得ARMA模型模擬大跨結(jié)構(gòu)風(fēng)速時(shí)程的效率得到進(jìn)一步提高，模擬精度也符合工程技術(shù)的要求。然后又將基于隨機(jī)變量在空間和時(shí)間上的相關(guān)結(jié)構(gòu)而建立起來(lái)的時(shí)空Kriging插值法用于結(jié)構(gòu)控制點(diǎn)臨近區(qū)域的風(fēng)速時(shí)程的模擬。模擬過(guò)程快速，精度滿(mǎn)足工程要求，可將此聯(lián)合方法推廣至復(fù)雜的大跨空間結(jié)構(gòu)、超高層建筑等含有大量節(jié)點(diǎn)結(jié)構(gòu)的風(fēng)速時(shí)程的模擬。

［1］曾憲武，韓大建.大跨度橋梁風(fēng)場(chǎng)模擬方法對(duì)比研究［J］.地震工程與工程振動(dòng)，2004，24（1）：135－140.

ZENG Xian-wu，HAN Da-jian.A comparative study on simulation methods for wind field of long-span bridges［J］.Earthquake Engineering And Engineering Vibration，2004，24（1）：135－140.

［2］李 莎，舒 紅，徐正全.東北三省月降水量的時(shí)空克里金插值研究［J］.水文，2011，31（3）：31－35.

LI Sha，SHU Hong，XU Zheng-quan.Study on spatialtemporal Kriging interpolation ofmonthly precipitation in three provincesofnortheast China［J］.Journal of China Hydrology，2011，31（3）：31－35.

［3］李 莎，舒 紅，徐正全.利用時(shí)空Kriging進(jìn)行氣溫插值研究［J］.武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2012，37（2）：237－240.

LI Sha，SHU Hong，XU Zheng-quan.Interpolation of temperature based on spatial-temporal Kriging［J］.Geomatics and Information Science of Wuhan University，2012，37（2）：237－240.

［4］羅俊杰，韓大建.大跨度結(jié)構(gòu)三維隨機(jī)脈動(dòng)風(fēng)場(chǎng)的模擬方法［J］.沖擊與振動(dòng)，2008，27（3）：87－91.

LUO Jun-jie，HAN Da-jian.Simulation method for 3－d stochastic wind field around long-span structures［J］.Journal of Vibration And Shock，2008，27（3）：87－91.

［5］《公路橋梁抗風(fēng)設(shè)計(jì)指南》編寫(xiě)組.公路橋梁抗風(fēng)設(shè)計(jì)指南［M］.北京：人民交通出版社，1996.

［6］李春祥，談雅雅，李錦華.基于ARMA模型模擬高架橋的脈動(dòng)風(fēng)速時(shí)程［J］.振動(dòng)與沖擊，2009，28（6）：46－51.

LI Chun-xiang，TAN Ya-ya，LI Jin-hua.Simulation offluctuating wind speed time series applied on overpass bridges with resorting to ARMA model［J］.Journal of Vibration and Shock，2009，28（6）：46－51.

［7］賀擁軍，周緒紅，王海順.超大跨屋面結(jié)構(gòu)風(fēng)速時(shí)程的數(shù)值模擬研究［J］.湖南大學(xué)學(xué)報(bào)（自然科學(xué)版），2008，35（9）：1－5.

HE Yong-jun，ZHOU Xu-hong，WANG Hai-shun.Numerical simulation of wind speed on super-long-span roof structures［J］.Journal of Hunan University：Natural Science，2008，35（9）：1－5.

［8］De Cesare L，Myers D，Posa D.Estimating and modeling space time correlation structures［J］.Stat.Probab.Lett.，2001，51（1）：9－14.

［9］De Cesare L，Myers D，Posa D.Product-sum covariance for space-time modeling，an environmental application［J］.Environmetrics，2001，12：11－23.

［10］徐愛(ài)萍，胡 力，舒 紅.空間克里金插值的時(shí)空擴(kuò)展與實(shí)現(xiàn)［J］.計(jì)算機(jī)應(yīng)用，2011，31（1）：273－276.

XU Ai-pin，HU Li，SHU Hong.Extension and implementation from spatial-only to spatio temporal Kriging interpolation［J］.Journal of Computer Applications，2011，31（1）：273－276.

Simulation of fluctuating w ind speed time-history series for large-span structures based on ARMA model and spatial-tem poral K riging interpolation

ZHOU Bin-bin，CAIJian-guo，F(xiàn)ENG Jian
（College of Civil Engineering，Southeast University，Nanjing 210096，China）

ARMA model was employed to simulate fluctuating wind speed time-history series for large-span structures with time and space correlativity.Whittle iteration method was improved to derive matrices of auto-regressive and moving-regressive coefficients.The corresponding computational program was developed with MATLAB software.Numerical examples indicated that these two methods are efficient and speedy.In addition，the spatial-temporal Kriging interpolation method was also proposed to simulate the fluctuating wind speed time-history series of unknown pointes near the simulated points based on ARMA model.Simulation results showed that the approach is swift and accurate and available to simulate fluctuating wind speed time-history series for complex long-span space structures and high-rise buildings containing plenty of joints.

numerical simulation；fluctuating wind speed time-history series；ARMA model；Whittle iteration；spatial-temporal Kriging interpolation

TU311

A

2013－01－10 修改稿收到日期：2013－03－13

周彬彬男，博士生，1986年6月生

馮 健男，教授，1963年10月生