黄埃散漫风萧索，云栈萦纡登剑阁

……

风吹仙袂飘飘举，犹似霓裳羽衣舞

——《长恨歌》

当你沉浸在“此恨绵绵无绝期”的爱恨情仇中时可曾注意到白居易描写的风的变幻？

什么？你才没工夫注意这些小事呢？各位结构大师，怎么可以连风荷载的变化都注意不到呢！

知己知彼百战不殆，抗风需先知风，今天小编带你领（模）略（拟）那让人捉摸不住的风。

通常认为空间一点的风速为平均风与脉动风之和，平均风可由经验公式的指数律分布函数和对数律分布函数两种方式来描述，小编不再一一赘述。而关于脉动风，小编可就要与大家一起捉摸一下了。

一、脉动风速的模拟方法有哪些？

经过社会主义的洗礼，我们都晓得实践是检验真理的唯一标准，真理的获得一定是从实践出发并回到实践中去。事实上我们对于风速的认知确实起源于“实践”：对实际风速的监测结果。与现场监测相似的是风洞试验，二者都是比较可靠的获取样本的方法，但都需要投入大量的时间或资金，随着计算机技术的飞速发展和人们对数值分析方法的深入研究，人工计算机模拟荷载的随机输入被广泛应用。人工计算机模拟风荷载可以考虑场地、风谱特征、建筑物的特点等条件的任意性，使模拟得到的荷载尽量接近结构的实际风力。目前主要的两种模拟风速的方法是：线性滤波法和谐波叠加法。前者基于线性滤波技术，也称为时间系列法，如状态空间法、自回归(Auto Regressive ，AR) 法、滑动平均(Moving Average，MA) 法、自回归滑动平均(Auto Regressive Moving Average ，ARMA)法等; 后者基于三角级数求和，也称为频谱表示法，如Constant Amplitude Wave Superposition(CAWS)法、Weighted Amplitude Wave Superposition(WAWS)法等。

今天小编就从线型滤波法说起，如果没有接触过频域时域相关的知识，你会像小编刚学习时一样面对佶屈聱牙的专业术语不知所措，这里小编先多嘴一句：线型滤波法的关键在于滤波功能，在风速模拟中可以将均值为零的白噪声随机序列通过滤波器，使其输出为具有指定谱特征的随机过程。而线性滤波法中的自回归 ( Auto-Regressive, AR)模型因其计算量小、速度快,广泛用于随机振动和时间系列分析中。

二、自回归（AR）模型是怎样的呢？

将t时刻的脉动风速时程 v (t) 表示为由之前若干个时刻时程的线性组合，再加上一个t时刻的独立随机过程。其模拟的脉动风速时程具有时间相关性，用如下公式表示：

式中：v(t) 为 t 时刻的脉动风速；φk 为自回归系数，p为AR模型的阶数；v(t-k∆t)为t时刻之前k个时刻的脉动风速；∆t为脉动风速的时间步长； N(t)是均值为 0、方差为（σN）^2的独立随机过程。

将上式两端同时右乘v(t-j∆t)得：

式中：j=1,2…,p。对上式左右两端同时取数学期望，结合自相关函数的性质：随机过程x(t)的自相关函数定义为x(t)x(t+τ)的平均值。如果该过程是平稳随机过程，则其平均值E[x(t)x(t+τ)]就与时间的绝对值无关，而只取决于时间差，在这里即表示为：

因N(t)的均值为 0 且与脉动风速v(t)独立，于是可得相关函数R(j·∆t)与自回归系数φk之间的如下关系：

式中的自回归系数φk由下面的方程组确定：

其中，R(j·∆t)为脉动风速在t=j∆t的相关函数，其值可由目标功率谱密度函数经傅里叶变换求得，可按下式计算：

将式(1)两端同时右乘v(t)，可得：

将式(1)代入式(7)，可得

对上式左右两端取数学期望，可得：

结合自相关函数的性质可得：

在线性滤波法的计算中，AR 模型的阶数p是影响整个计算过程准确度和效率的决定性因素。当模型阶数选择较少时，运行速度快，但可能忽略掉对v(t)有较大影响的项数，导致计算精度不足；当模型阶数选择较多时，计算速度慢，并且距离 t 时刻越远，自回归系数越小，这样虽然增加了阶数，但计算的精度却没有得到明显提升，因此就需要选择合适的模型阶数。通常采用残差平方极小化法确定模型阶数，定义AR模型的残差平方和为：

从 p=1 阶模型开始，逐渐增加模型阶数并计算其相应的残差平方和，比较两相邻模型的残差平方和大小，当两模型之间的残差平方和变化不显著时，就可以确定计算时需要的模型阶数p的值。

三、道理我都听明（糊）白（涂）了，那要怎么实现自回归模型呢？

小编以matlab为例来实现AR模型的数值模拟：

1、首先选取一个风速的功率密度谱，小编以风工程中最常用的Davenport功率密度谱为例：

式中：S(n)—脉动风速功率密度；

ω—脉动风频率；

x=1200n/V10；

k—地面粗糙度系数；

是Davenport于1961年在不同地点、不同条件下测得的90多次强风记录基础上，获得风速的相关曲线并建立数学表达式，然后将相关函数通过傅里叶变换得到脉动风速功率谱密度函数的表达式。大多数国家（包括我国）建筑荷载规范都采用此水平风速谱公式，功率谱密度函数本质上表示的是某一频率上脉动风的能量大小的分布。

2、将目标功率谱代入式(6)确定相关函数R(τ)，将计算得到的相关函数代入式(5)求出自回归系数φk。

Matlab步骤：1）将功率谱密度函数设为一个m函数，便于谱的替换；

2）建立m函数对式（6）进行数值积分得到相关函数的数值解；

3）循环语句建立自相关函数矩阵，求解式（5）得到自回归系数矩阵；

3、再将求得的自回归系数代入式(10)，便可求得σN。

4、将求得的自回归系数φk和σN代入式(1)，并假定初始时刻之前的脉动风速为0，即可求出一点的脉动风速时程。

求解思路就是这样了，具体在计算机上采用何种软件、语言去实现就要看个人喜好了，小编取标准平均风速为12m/s，5阶AR模型的模拟结果如下图：

文中只展示了部分程序，回复“AR模型”可以获取matlab中实现模拟单点脉动风速AR模型的全部程序哦~

线性滤波法中的自回归模拟方法的介绍就到这里了，另一种应用范围很广的谐波叠加法且听小编下回分解。

参考文献：

[1]Ghosh R , Joshi Y . Proper Orthogonal Decomposition-Based Modeling Framework for Improving Spatial Resolution of Measured Temperature Data[J]. IEEE Transactions on Components, Packaging and Manufacturing Technology, 2014, 4(5):848-858.

[2]陈小波. 海上风力发电塔脉动风速时程数值模拟[J]. 中国电机工程学报, 28(32).

[3]舒新玲, 周岱. 风速时程AR模型及其快速实现[J]. 空间结构, 2003, 9(4):27-32.

[4]莫继华, 何炎平, 李勇刚, et al. 近海风电机组单桩式支撑结构疲劳分析[J]. 上海交通大学学报（自然版）, 2011, 45(04): 565-569.