基于AR模型的空间脉动风速时程模拟方法研究

赵海霞　（南京航空航天大学金城学院，江苏 南京 211156）



基于AR模型的空间脉动风速时程模拟方法研究

赵海霞（南京航空航天大学金城学院，江苏 南京 211156）

空间脉动风速时程的模拟是桥梁抖振时域计算的前提。利用基于自回归（AR）模型的线性回归滤波器法，编制了模拟脉动风速的Matlab计算程序，模拟了一座大跨连续刚构桥最大悬臂施工阶段多点的水平和竖向脉动风速，并进行了功率谱和相关性检验，发现模拟效果较好。

脉动风速；AR模型；线性回归滤波器法

1　研究背景

随着我国桥梁工程的建设不断向大跨度方向发展，桥梁的长细化使其刚度和阻尼不断下降，导致结构对风的敏感性不断增加，对桥梁结构进行风效应分析也逐渐成为桥梁结构分析中的重要一环。由于时域分析方法可以直接计算频域分析中难以处理的非线性响应，因而对大跨度桥梁的风载作用下的抖振行为进行非线性时域分析日益受到人们的关注。在抖振响应时域分析中，首先要根据目标功率谱函数人工模拟空间脉动风场，从而获得离散的风速时程和离散的抖振力。

在人工模拟空间脉动风场方面，基于自回归（Auto-Regressive，AR）模型的线性回归滤波器法具有速度快、计算量小的特点，且精度较高，因此得到了广泛应用[1][2]。该模型将均值为零的白噪声随机系列通过线性滤波器，使其输出为具有指定谱特征的平稳随机过程[3]。本文首先详细介绍了该方法的计算原理和过程，利用Matlab软件编制了相应的计算程序；采用该方法模拟了一座大跨连续刚构桥最大悬臂施工阶段多点的水平和竖向脉动风速，并进行了功率谱和相关性检验。

2　自回归（AR）模型

2.1求回归系数

采用AR法推广到模拟多维风速过程的技术，M个相关的随机风过程[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、具有给定协方差的正态分布随机过程，i ＝1，…，M；[ψk]为M×M阶矩阵，k=1，…，p；p 为AR模型的阶数，一般取4或5。

对任一空间点i（i=1，…，M）具有时间差的随机风过程ui（t）与ui（t-kΔt）的协方差可表示为：

由于ui（t）与ui（t-kΔt）为均值0的平稳随机风过程，其协方差的值仅为时间差的函数，式（2）可改写为：

在式（1）同时右乘[u（t-jΔt）]=[u1（t-jΔt），…，uM（t-jΔt）]，j=1，…，p，并两边同时取数学期望（均值），考虑到[N（t）]的均值为0，且与随机风过程ui（t）独立，以及协方差Ru（jΔt）为偶函数，可得到协方差Ru（jΔt）与回归系数ψk之间的关系，写成矩阵形式，有：

其中

根据随机振动理论[3]，功率谱密度与相关函数（协方差）之间符合维纳-辛钦（Wiener-Khintchine）公式，即：

通过rjk（n）考虑风速时程的空间相关特性，其三维表达式为：

式中：Cx、Cy和Cz分别表示空间任意两点左右、上下和前后的衰减系数，一般分别取为16、8和10[4]；（xi，yi，zi）、（xk，yk，zk）分别为空间i、k点的三维坐标，i、k=1，…，M；x、y、z分别为垂直于来流的水平方向、来流方向和竖向；分别表示第i点和第k点的平均风速。

求解式（4）给出的线性方程组，可以得到回归系数矩阵[ψ]。

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且彼此相互独立的正态随机过程，i=1，…，M。

2.3求最终的M个随机过程

求出回归系数矩阵[ψ]及[RN]后，可按式（1）求解出M个空间相关的随机风过程。

将式（1）按时间间隔Δt离散化，分别考虑三种情况，即ui（t）为①偶函数；②奇函数；③当t<0时，ui（t）＝0。由此可得出不同的矩阵方程形式。其中，以第三种假设计算起来最为方便，其递推的矩阵表达式为：

从而得到M个具有时间、空间相关、时间间隔Δt的离散脉动风速时程向量。

为了避免模拟结果失真，Δt必须满足以下条件：

式中：ωup是截断频率（rad/s）。

3　风场模拟的Matlab程序

根据线性回归滤波器法的自回归（AR）模型的思路，所编制的Matlab程序流程如图1所示。

图1　自回归模型生成风速时程流程图

4　算例

以某连续刚构桥最大悬臂施工阶段为对象，选定的模拟点即风荷载加载节点共23点，其中主梁上17点，墩上6点。图2所示为部分模拟点位置示意图。主梁上的模拟点模拟水平脉动风速和竖向脉动风速，墩上的模拟点仅模拟水平脉动风速。

图2　风速模拟点示意图（单位：m）

顺风向风速模拟目标谱采用Kaimal水平脉动风速谱[5]：

竖直向风速模拟目标谱采用Panofsky谱：

式中：Su（n）、Sw（n）分别为水平向和竖向脉动风速功率谱；n为脉动风频率（Hz）；Z为有效高度；U（Z）为高度Z处的平均风速；u*为气流摩阻速度（亦称剪切速度）；K是无量纲常数，K≈0.4；Zd为零平面高度；为周围建筑物的平均高度；kd是地面阻力系数；Z0是地面粗糙长度。

考虑上述风速模拟的目标谱为非均方规一的单边功率谱，通过与实际脉动风速的均方差比较，可以换算得到圆频率表达的风谱形式：

模拟的其他参数有：场地类别为Ⅱ类；场地基本风速为24.1m/s[6]，其他高度处的风速按照抗风规范的指数分布，风速剖面无量纲幂指数α=0.16；地面粗糙长度Z0取为0.05；零平面高度Zd取为0；截至频率取为2Hz；频率等分数取为12000；采样时距取为0.1s，共取12000步，1200s。分别选取左悬臂端（点1）、距左悬臂端8m处（点2）、主梁0号块中点（点9）、右悬臂端处（点17）和墩左肢中部（点20）共五点作为展示对象。图3是用自回归（AR）模型模拟的这五点的水平风速时程，图4是其竖向风速时程（不含点20）。

图3　各点水平脉动风速

图4　各点竖向脉动风速

图5～图6是分别用改进的平均周期图法和最大熵法[5]对所模拟的点1风速时程做的谱估计；图7所示为点1的自相关函数检验，图8所示点1和点2的互相关函数的检验。从图中我们可以看出，模拟值与目标值之间相差很小，可见模拟的结果良好。

图5　点1功率谱密度函数（最大熵法）

图6　点1功率谱密度函数（多周期图法）

图7　点1自相关函数

图8　点1、2互相关函数

4　结语

本文详细介绍了基于自回归（Auto-Regressive，AR）模型的线性回归滤波器法的计算原理和过程，利用Matlab软件编制了相应的计算程序。采用该程序模拟了一座大跨连续刚构桥最大悬臂施工阶段多点的水平和竖向脉动风速，并进行了功率谱和相关性检验，发现模拟效果较好。因此本文的计算方法和程序具有一定的推广价值。

[1]舒新玲，周岱.风速时程AR模型及其快速实现[J].空间结构，2003，9 （4）：27-32.

[2]李元齐，董石麟.大跨度空间结构风荷载模拟技术研究及程序编制[J].空间结构，2001，7（3）：3-11.

[3]俞载道，曹国敖.随机振动理论及其应用[M].上海:同济大学出版社，1988.

[4]边建烽，魏德敏.大跨空间结构风速时程的数值模拟理论[J].暨南大学学报（自然科学版），2005，26（1）：87-90.

[5]JTG T D60-01-2004，公路桥梁抗风设计规范[S].

[6]张文明.大跨连续刚构桥最大悬臂施工阶段风致响应分析[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-），女，毕业于西安建筑科技大学，硕士；讲师，南京航空航天大学金城学院土木工程教研室。