基于多维AR模型的桥梁随机风场模拟
张田，夏禾，郭薇薇
(北京交通大学 土木建筑工程学院，北京，100044)
摘要：基于自然风特性，通过考虑结构节点间风速时程的空间相关性，采用多维AR模型模拟主梁和桥墩节点随机脉动风速时程，利用FPE准则和AIC准则确定模型阶数，并对模拟过程中的自回归顺序、功率谱密度等问题进行研究。对兰新二线铁路白杨河大桥采用多维AR模型模拟各节点的脉动风速时程，结果表明：当AR模型阶数为4时，模拟功率谱与目标功率谱吻合较好；当自回归顺序颠倒时，模拟功率谱明显偏离目标功率谱。
关键词：桥梁；脉动风；随机风场模拟；AR模型；定阶
中图分类号：U24 文献标志码：A 文章编号：1672-7207(2012)03-1114-08
Simulation of bridge stochastic wind field using multi-variate Auto-Regressive model
ZHANG Tian, XIA He, GUO Wei-wei
(School of Civil Engineering, Beijing Jiaotong University, Beijing 100044, China)
Abstract: Based on the natural wind properties and the correlativity of nodal wind speed time history, a multi-variate auto-regressive (MVAR) model was used to simulate the time history of fluctuating wind for the main beam and bridge pier, and FPE and AIC rules were employed to determine the order of MVAR model, then related problems such as auto-regressive sequence and power spectrum density were discussed. Taking the Baiyang River Bridge on the Second Lanzhou-Urumqi railway as an example, the wind speed time histories were simulated with the MVAR model. The results show that the simulated power spectrum density is consistent with the target one when the model order is equal to 4, however, it shows obvious deviation when the auto-regressive sequence is reversed.
Key words: bridge; fluctuating wind; stochastic wind field simulation; Auto-Regressive (AR) model; order determination
要在时域内对桥梁进行风致振动分析，进而分析车辆-桥梁系统的动力响应，首先必须得到桥梁上脉动风速的时间历程，即首先要解决桥梁脉动风速场的随机模拟[1]。根据前人的试验研究，风速场的非平稳性常表现初始阶段，因此，自然风的脉动近似处理为各态历经的平稳高斯随机过程[2]。目前模拟多变量随机过程主要使用的是基于线性滤波技术的回归方法和基于三角级数的谐波叠加法，另有小波模拟[3]、本征正交分解法(POD)[4]等方法。其中，谐波叠加法算法具有简单明确、精度高和稳定性好等特征，在工程实际中已经得到了广泛的应用，常用的有Shinozuka等[5]提出的谐波合成法、Deodatis等[6]基于频率双索引的谐波合成法以及Cao等[7]提出的特殊情况下应用的快速谱分析方法，但对模拟风速点较多时计算量巨大。而线性滤波技术是广泛应用于随机振动和时间序列分析中的一种分析手段，近年来，线性滤波法中的自回归(AR)模型因其计算量小和速度快，被广泛应用于脉动风场模拟中，目前多用于平面大尺度的空间大跨房屋建筑结构的脉动风场模拟[8]，而对于立面大尺度线形结构的桥梁风场模拟应用还比较少。AR模型将均值为0的白噪声随机系列通过线性滤波器，使其输出为具有指定谱特征的平稳随机过程。一般常用的方法是先生成一系列互不相关的脉动风速时程，然后再考虑各点的之间的空间相关性，因此，互相关函数未考虑时滞的影响[9-12]；本文作者基于自然风特性，认为脉动风场是联合平稳的，考虑桥梁结构不同节点之间的空间相关性和时滞的影响，采用多维AR模型模拟具有随机性、时空相关性的脉动风速时程，利用最终预测误差(FPE)定阶准则或信息论(AIC)准则识别适宜的模型阶数。
1 随机脉动风速时程的多维AR模型
风速观测记录表明瞬时风速包含2种成分：周期在10 min以上的平均风和周期在几秒的脉动风[8]。在满足工程计算精度要求的前提下，可对风速时程模拟作以下假定：(1) 任意一点处平均风速不随时间改变；(2) 脉动风速时程是各态历经、零均值平稳随机过程；(3) 各点风速时程间具有空间相关性；(4) 不同高度处的脉动风速同相，即脉动风速的互相关函数偶对称。
风速时程本质上是随机时间系列，桥梁作为大尺度的线形空间结构，其风速时程可用多维AR模型模拟。m个相关脉动风速时程v(X, Y, Z, t)=[v1(x1, y1, z1, t) v2(x2, y2, z2, t) …vm (xm , ym, zm, t)]可由下式生成：
(1)
式中：
，
，
，
为空间第i点坐标(i=1, 2, 3, …, m)；p为AR模型阶数；?t为模拟风速时程的时间步长；φk为AR模型自回归系数矩阵，为m×m阶方阵，k=1, 2, …, p；N(t)为零均值、具有给定方差的独立随机过程。
为了方便起见，v(X, Y, Z, t)简写为v(t)。根据风速时程模拟的假定，将式(1)两边同时右乘
，得：
(2)
式中：j=0, 1, 2, …, p。对上式两边同时取数学期望，且考虑到N(t)的均值为0，与随机过程v(t)独立，并结合自相关函数的如下性质：
(3)
(4)
其中“E[ ]”为数学期望算子。于是可以得到相关函数
与自回归系数
之间的关系如下：
( j=1, 2, …, p) (5)
将式(1)两边同右乘v T(t)，取数学期望，则有
(6)
写成矩阵形式，有：
(7)
式中：
，I为m阶单位矩阵；RN为m×m阶协方差矩阵；Op为pm×m阶矩阵，其元素全部为0；R为(p+1)m×(p+1)m阶自相关托普利兹(Toeplitz)矩阵，写成分块矩阵的形式，则有：
(8)
式中：
为m×m阶方阵(s=1, 2, …, p+1；t=1, 2, …, p+1；j=0, 1, 2, …, p)；其值可由维纳-辛钦(Wiener-Khintchine)公式求得：
(9)
其中：f为脉动风速频率，Hz；
在i＝h时为脉动风速自谱密度函数；在i≠h时为脉动风速互谱密度函数，可由脉动风速自谱密度函数
和相关性函数
确定(i=1, 2, 3, …, m；h=1, 2, 3, …, m)。
风场中，不同点风速的互功率谱可以写为：
(10)
式中：
为相位角对f的函数，根据风速时程模拟的假设(3)和(4)，
。
脉动风的功率谱由强风观测记录得到，方法有2种：一种是把强风观测记录经过相关分析，先获得风速的相关曲线，建立相关曲线的数学表达式，再通过傅里叶变换求得；另一种是把强风记录通过超低频滤波器，直接得到风速的功率谱曲线。目前国内外学者提出了许多脉动风功率谱，其中Kaimal随高度变化的水平风谱以及Panofsky竖向风谱列入了我国《公路桥梁抗风设计规范》(JTG D60-01—2004)[13]。本文选用这2种谱作为目标谱，其具体形式为：
(11)
(12)
式中：
；
；
；z0为地面粗糙高度，m；
为周围建筑物平均高度，m；下标u和w分别表示水平向和竖向；
和
分别为脉动风的水平顺风向及竖直方向的功率谱密度函数；
为气流摩阻速度，也称剪切速度，m/s；K为无量纲常数，K≈0.4；z为地面或水面以上的高度，m；
为高度z处的平均风速，m/s，可按对数律来计算，即
(13)
其中：v1为标准高度z1(一般为10 m)处的平均风速。
空间相关性系数反映了不同点位上风谱的空间相关程度，对于迎风面尺度较大的结构需要考虑脉动风的空间相关性。对于侧向和竖向尺寸较大的结构，例如高桥墩桥梁结构，Davenport[14]建议的空间相关系数r(f)为：
(14)
式中：
和
分别为第i点和第h点的平均风速；(xi, zi)和(xh, zh)分别为空间点i和h的坐标(i=1, 2, 3, …, m；h=1, 2, 3, …, m)；Cx和Cz分别为空间任意两点左右、前后的衰减系数，通过实测或实验确定，Simiu等[15]建议：Cx=16，Cz=10。
独立随机过程向量N(t)可写为：
(15)
其中：
，
为均值是0、方差为1且相互独立的正态随机过程(i=1, 2, 3, …, m)；L为m阶下三角矩阵，通过对m×m阶协方差矩阵RN的乔里斯基(Cholesky)分解确定：
(16)
2 AR模型的定阶
研究表明[16]：低阶的AR模型即可较好的模拟随机脉动风速时程。因此，应从AR(1)模型开始，逐步增加模型阶数。模型定阶方法有如下2种。
(1) 最终预测误差定阶准则
最终预测误差函数FPE定义为：
(17)
其中：
为残差，根据下述方法计算。
计算样本自协方差函数：
，p=0, 1, 2, …, P (18)
其中：t为样本数；P为与t有关的较大正整数。
按如下Levioson递推公式计算
：
(19)
使用时可能遇到以下情况：① 若FPE(k)(k=1, 2, …)一直上升，则可判定p=1作为模型的阶数；② 若FPE(k)一直下降，则可认为过程不是有限阶AR过程，只能考虑p阶自回归逼近问题；③ 若FPE(k)一直下降很快，到某个整数p时下降缓慢，也可以把该p作为模型阶数；④ 若FPE(k)剧烈地上下摆动，则可以通过增大样本容量再进行定阶。
(2) 信息论准则
该准则是极大似然原理应用于时间序列的统计假设试验而得出的。该准则指出，最优阶次p的选择应使式(2)为最小值。
(20)
实际上，也可以表述为：当欲从一组可供选择的模型中选择一个最佳模型时，选取使AIC准则函数值为最小的模型是适宜的。AR(p)模型的AIC准则函 数为：
(21)
除了用上述方法求得外，还可以用最小二乘估计或极大似然估计法在自有变化参数个数p的前提下对数据拟合一个AR模型的残差方差估计。即从AR(1)开始，AR(1), AR(2), …, AR(p)，一共建立p个模型，计算各个模型的AIC准则函数值，比较2个相邻模型间AIC准则函数值的大小，若AR(p-1)与AR(p)间的AIC准则函数值变化不显著时，即可确定模型的阶数为p。
3 工程实例
以兰新铁路第二双线的白杨河大桥为例，模拟该桥桥墩与桥面主梁上的脉动风速时程。白杨河大桥位于线路的三十里风区，为(40+3×64+40)预应力混凝土连续梁桥，最高墩高49 m，其模拟风速点的划分如图1所示。其中桥墩顶点与主梁处节点重合，即桥墩顶点处的风速时程与桥墩相同位置处的风速时程相同。风速时程模拟时的主要参数如表1所示。
根据文中给出的定阶准则，计算出AR模型取各阶数时FPE函数值和AIC准则函数值列于表2。从表2可知：该工程实例脉动风速模拟的最优阶数为p=4。然后采用上文的理论模型编程模拟该桥梁各节点脉动风速时程，其中主梁的节点7，8和17，以及桥墩节点21和41的瞬时脉动风速时程曲线如图2所示，其中
为该点的平均风速。从图2可以看出：模拟脉动风速的随机性较好；并计算比较了节点8，17和41的脉动风速模拟功率谱与目标功率谱，结果如图3所示。由图3可以看出：3个节点位置处模拟功率谱与目标功率谱趋势及值的大小较为一致，该方法模拟脉动风场是行之有效的。同时计算了主梁节点8的自相关函数，主梁节点8与主梁节点7的互相关函数以及主梁节点8与主梁节点17的互相关函数；桥墩节点21的自相关函数，桥墩节点21与桥墩节点41的互相关函数，如图4所示。由图4可以看出：相邻点的相干性较强，相距越远，相干性越弱；并且桥墩上风速点的相干性减弱比主梁上风速点快。
表1 模拟脉动风速时程曲线时的计算参数
Table 1 Calculation parameters for simulation of time history curves of fluctuating wind
图1 白杨河大桥桥型及风速点划分示意图(单位：cm)
Fig.1 Wind velocity points distribution and bridge type of Baiyang River Bridge
表2 AR模型的定阶
Table 2 Order determination of AR model
图2 模拟点的瞬时风速时程曲线
Fig.2 Instantaneous wind velocity curves of simulated points
图3 模拟节点与目标功率谱比较
Fig.3 Comparisons of power spectra between simulated points and target
图4 相关函数比较
Fig.4 Comparison of correlation functions
4 关于AR模型模拟脉动风速时程的相关讨论
4.1 协方差矩阵RN
由于AR模型中N(t)为零均值、具有给定方差的随机过程，即N(t)中包含2个方面信息：一方面表示方差；另一方面表示独立的随机过程。由于空间点数不止一个，表示方差的部分由协方差矩阵RN 组成。然而求解RN时，各种文献和著作中表达不一致，这会使人
感到困惑。如文献[17]中，
，而文献[18]中却为
，实际上
上述2种表达都是正确的，这是鉴于对AR模型的表达不同而引起，只要前后保持统一即可。
4.2 自回归的顺序问题
从AR模型的表达式可以看出：是第1个自回归参数对应最后一个输出，第2个自回归参数对应倒数第2个输出，以此类推，最后一个自回归参数对应倒数第p个输出。该顺序问题似乎无需强调，但之所以在此单独提出，是因为作者在用该方法编制程序时偶然发现，若将以上自回归顺序颠倒，模拟出的脉动风速时程样本、自相关估计及其自功率谱估计和正确顺序时的表现出明显不同的特点。自回归顺序颠倒时的风速时程样本及其自相关估计明显地表现出“黑”，即所得风速及其自相关估计变化频率非常快，如图5所示。自回归顺序颠倒时风速时程样本的功率谱密度较自回归顺序正确时明显偏离目标功率谱，而且在高频时会发生突变，如图6所示。
图5 自回归顺序颠倒后模拟节点8的风速时程曲线
Fig.5 Instantaneous wind velocity curve of simulated point with wrong auto-regressive sequence
图6 自回归顺序颠倒后的功率谱密度
Fig.6 Power spectrum density with wrong auto-regressive sequence
4.3 功率谱密度及归一化问题
功率谱密度函数(PSD)是通过均方值的谱密度来描述数据的频率结构。一个实的过程X(t)，
自功率谱密度函数
是它的自相关函数
的Fourier变换：
(22)
其Fourier逆变换即自相关函数：
(23)
根据其定义，引用Euler公式及式
，得：
(24)
(25)
由样本函数x(t)的Fourier变换计算遍历性过程X(t)的功率谱的公式为：
≥0 (26)
如果过程X(t)是不具遍历性的平稳过程，则其功率谱需用集合平均方法得到，即：
≥0 (27)
因此，自功率谱密度函数为实的非负偶函数。
上述的功率谱密度函数的定义范围是
，故称为双边谱，而工程上负频率无实际的物理意义，为此又定义了单边谱
，文献[19]给出的单边谱表达式如下：
(28)
即在
时，
；对于该问题，本文认为当时，
。因为前面已经证明自功率谱密度函数为实的非负偶函数，而偶函数只有在
时才有正负频率上对应的2个点；时谱密度值只有一个，即在处，谱密度值PSD应保持不变，如图7所示，“l”表示取该点的值。
图7 双边自谱与单边自谱示意图
Fig.7 Diagram of double-side and single-side auto-power spectrum
由于功率谱的纵坐标量纲为对象变量量纲的平方，但是实际中规范和一些经典的风速谱函数常选用折算频率
(各符号含义同式(11))作为横坐标，而不是通常的f作为横坐标，纵坐标也做了相应的折算，这样相当于把目标值归一化了，以便于比较分析，因此，求模拟的风速时程序列的功率谱函数时要考虑是否进行归一化处理，这与模拟时所采用的目标风谱有关。
5 结论
(1) 采用较低阶的多维AR模型，可以考虑桥梁各节点风速时程间的空间相关性，较好地模拟出各节点的风速时程曲线，且计算速度较快。
(2) 利用AR模型定阶准则，可以准确地确定模型的最优阶数。对本文的工程实例，选择风速时程模拟的最优模型阶数为p=4。
(3) 模拟出的节点风速时程与目标功率谱吻合较好，也通过了相关性检验，证明了所提出模拟方法的可靠性。
(4) 对AR模型的相关问题进行讨论，尤其是自回归顺序问题，当自回归顺序颠倒时，模拟的风速时程样本会表现出明显的“黑”，即风速变化频率非常快。此外，单边谱与双边谱之间的关系问题，即当时，应引起注意。
参考文献：
[1] 夏禾, 张楠. 车辆与结构动力相互作用[M]. 第2版. 北京: 科学出版社, 2005: 199-210.
XIA He, ZHANG Nan. Dynamic interaction of vehicles and structure[M]. 2nd ed. Beijing: Science Press, 2005: 199-210.
[2] 陈英俊, 于希哲. 风荷载计算[M]. 北京: 中国铁道出版社, 1998: 1-18.
CHEN Ying-jun, YU Xi-zhe. Wind load calculation[M]. Beijing: China Railway Publishing House, 1998: 1-18.
[3] 陈艾荣, 王毅. 基于小波方法的随机脉动风模拟[J]. 同济大学学报: 自然科学版, 2005, 33(4): 427-431.
CHEN Ai-rong, WANG Yi. Simulation of random fluctuating wind speed based on wavelet method[J]. Journal of Tongji University: Natural Science, 2005, 33(4): 427-431.
[4] 胡亮, 李黎, 樊剑. 基于特征正交分解的桥梁风场模拟[J]. 武汉理工大学学报: 交通科学与工程版, 2008, 32(1): 16-19.
HU Liang, LI Li, FAN Jian. Digital simulation of bridge wind fields based on spectral representation method with proper orthogonal decomposition[J]. Journal of Wuhan University of Technology: Transportation Science & Engineering, 2008, 32(1): 16-19.
[5] Shinozuka M, Jan C M. Digital simulation of random processes and its applications[J]. Journal of Sound and Vibration, 1972, 25(1): 111-128.
[6] Deodatis G. Simulation of ergodic multivariate stochastic processes[J].Journal of Engineering Mechanics, ASCE, 1996, 122(8): 778-787.
[7] CAO Ying-hong, XIANG Hai-fan, ZHOU Ying. Simulation of stochastic wind velocity field on long-span bridges[J]. Journal of Engineering Mechanics, ASCE, 2000, 126(1): 1-6.
[8] 赵臣, 张小刚, 吕伟平. 具有空间相关性风场的计算机模拟[J]. 空间结构, 1996, 2(2): 21-25.
ZHAO Chen, ZHANG Xiao-gang, L? Wei-ping. Computer modeling on wind histories with spatial correlativity[J]. Spatial Structures, 1996, 2(2): 21-25.
[9] 黄国辉, 徐国彬. 基于数值模拟法的张力膜结构风振响应分析[J]. 空间结构, 2005, 11(3): 56-59.
HUANG Guo-hui, XU Guo-bin. Analysis of wind induced response of tensile membrane structure based on digital simulation method[J]. Spatial Structures, 2005, 11(3): 56-59.
[10] 杨帆, 杜文风, 杨国忠, 等. 输电铁塔结构脉动风模拟AR法[J]. 科技资讯, 2008, 6(7): 42-43.
YANG Fan, DU Wen-feng, YANG Guo-zhong, et al. AR model method of fluctuating wind simulation for transmission tower structures[J]. Science & Technology Information, 2008, 6(7): 42-43.
[11] 刘学利, 王肇民. 高耸结构空间相关风场的模拟研究[J]. 四川建筑, 2000, 20(4): 45-47.
LIU Xue-li, WANG Zhao-min. Spatial correlation wind field simulation for high-rise structures[J]. Sichuan Architecture, 2000, 20(4): 45-47.
[12] 董军, 邓洪洲, 刘学利. 高层建筑脉动风荷载时程模拟的AR模型方法[J]. 南京建筑工程学院学报, 2000(2): 20-25.
DONG Jun, DENG Hong-zhou, LIU Xue-li. AR model method for time-history imitation of pulse wind load acting on high-rise buildings[J]. Journal of Nanjing Architectural and Civil Engineering Institute, 2000(2): 20-25.
[13] JTG/T D60-01—2004, 公路桥梁抗风设计规范[S].
JTG/T D60-01—2004, Wind resistant design specification for highway bridges [S].
[14] Davenport A G. The application of statistical concepts to the wind loading of structures[J]. Proc Inst Civ Engrs, 1961, 19(2): 449-472.
[15] Simiu E, Scanlan R H. Wind effects on structures[M]. New York: Wiley, 1978: 42.
[16] Lannuzzi A, Spinelli P. Artificial wind generation and structural response[J]. Journal of Structural Engineering, ASCE, 1987, 113(12): 43-57.
[17] 舒新玲, 周岱. 风速时程AR模型及其快速实现[J]. 空间结构, 2003, 9(4): 26-32.
SHU Xin-ling, ZHOU Dai. AR model of wind speed time series and its rapid implementation[J]. Spatial Structures, 2003, 9(4): 26-32.
[18] 李元齐, 董石麟. 大跨度空间结构风荷载模拟技术研究及程序编制[J]. 空间结构, 2001, 7(3): 3-11.
LI Yuan-qi, DONG Shi-lin. Random wind load simulation and computer program for large-span spatial structures[J]. Spatial Structures, 2001, 7(3): 3-11.
[19] 陈英俊, 甘幼琛, 于希哲. 结构随机振动[M]. 北京: 人民交通出版社, 1993: 77-78.
CHEN Ying-jun, GAN You-chen, YU Xi-zhe. Random vibration of structures[M]. Beijing: China Communications Press, 1993: 77-78.
(编辑 陈爱华)
收稿日期：2011-03-10；修回日期：2011-07-23
基金项目：国家自然科学基金资助项目(51078029)；中央高校基本科研业务费专项资金资助(C11JB00540)
通信作者：张田(1986-)，男，湖北荆州人，博士研究生，从事桥梁抗风研究；电话：13401052635；E-mail: saghb@126.com
