基于演化相位谱的脉动风速模拟

阎 启，李 杰，2

(1.同济大学 建筑工程系，上海 200092;2.同济大学 土木工程防灾国家重点实验室，上海 200092)

结构在风荷载作用下的计算一般有频域和时域两种方法。频域方法速度快，但是只能用于线性阶段的计算;而对于输电线路、大跨桥梁、高层建筑等具有明显非线性特征的结构，时域分析能够更为准确的把握结构响应特征。

时域分析采用的风速时程通常由数字模拟的方法生成［1]，其中谱表现法(Spectral representation method)是一种较为成熟的数字模拟方法。经历了从1954年［2]至今的发展，谱表现法从最初只能模拟一维、单变量、平稳的随机过程，已发展为可以广泛的模拟多维、多变量、非平稳随机过程的方法［3－5]，并在脉动风场模拟方面得到了广泛的应用［1，6]。谱表现法的实质，是利用具有随机相位的谐波叠加、对具有目标功率谱的随机过程进行模拟。但仔细分析发现:谱表现法在基础上至少存在以下三点问题。首先，谱表现法本质上是基于二阶矩的模拟。在进行随机过程模拟时，各谐波的幅值是利用功率谱密度进行计算，对目标随机过程的模拟精度也是通过功率谱密度的符合程度来体现。由于功率谱密度在本质上属于随机过程的二阶数值特征，因此，所模拟随机过程的概率信息不会高于二阶。其次，谱表现法存在无法重现样本的困难。由于功率谱本质上属于集合特征［7]，因此基于功率谱模拟的随机过程无法重现现实中的随机过程样本。再次，谱表现方法中需要用到初始相位，通常是在不同频率处取区间［0，2π)内均匀分布并且相互独立的随机变量作为初始相位［3－5]。这种做法使得在实际风场模拟中随机变量的数目巨大(400～600个)，而不同频率相位之间、相位与能量之间的可能的内在联系并未被考虑。上述问题带来的后果是:应用谱表现法得到的风速时程进行结构随机响应分析，不仅计算工作量大，也很难得到响应的概率密度信息，进而，造成了结构可靠度分析的困难。

随机Fourier函数模型的提出为动力激励的建模与模拟提供了新的途径［8]。对地震动、脉动风速以及海浪等动力激励的时程，应用Fourier变换，可以将其分解为幅值部分和相位部分，称为Fourier幅值谱和Fourier相位谱。幅值谱是对动力激励能量的一种谱分解，基于大气湍流的物理机制，可以得到脉动风速随机Fourier幅值谱的物理模型［9]。相位谱则主要控制时程的波形、反映过程的非高斯性［10]。作为随机Fourier谱中十分重要的一部分，建立合理的相位谱模型，对于脉动风速的模拟以及结构抗风可靠度的计算，具有至关重要的意义。

本文首先从概念上明确了相位谱和相位谱主值的区别，基于对风场湍流物理背景的分析，提出了相位演化速度的概念，并由此建立了相位的零点演化时间Te和相位谱一一对应的关系，依据大量实际风速测量数据，建立了Te的概率模型。结合脉动风速的Fourier幅值谱、应用Fourier逆变换进行了脉动风速的模拟。与实测结果的对比，表明了本文建议方法的可行性。

1 相位谱和相位谱主值

设脉动风速时程为u(t)，其离散Fourier变换为:

其中n为自然频率，T为脉动风速的持时，d t为采样间隔，与采样频率Fs互为倒数，F(n)称为脉动风速时程u(t)的Fourier谱。乘以是为了取单边谱。F(n)是一组复数，可以写为幅值和相位的形式:

图1 实测脉动风速时程及其Fourier幅值谱和相位谱Fig.1 Measured fluctuating wind speed and corresponding Fourier amplitude spectrum and phase spectrum

2 演化相位谱

2.1 相位的演化速度

由于测量条件的限制，很难得到某一时刻空间中各点的湍流速度。Taylor冻结假定将以空间变量为背景的湍流假设为空间一点以时间为变量的湍流［12]。根据这一假定，同一空间点、不同时刻的脉动风速时程(自变量为时间)和同一时刻、不同空间点的脉动风场分布(自变量为位置)是相同的。这里，不同空间点指沿主风方向长度为L的范围，L=UT，U为平均风速，T是测量脉动风速的时距。自Taylor冻结假定提出后，几乎所有的湍流测量都以此为基础［13]。

众所周知，风速仪记录到的是空间一点处不同气体质点的流动速度，风速记录减去平均风速即得到脉动风速时程。由Taylor假定可以推论:如果风速仪以和平均风速同样的速度沿主风方向前进，那么将记录到同一气体质点在空间不同位置的速度，该速度与前述脉动风速时程相同。此时，脉动风速时程可视为是同一气体质点在一以平均风速值大小前进的参考系内的纵向“振动”速度。

在物理上，这样的空气质点振动速度可视为是一系列不同尺度、不同频率涡旋振动的叠加。涡旋的特征振动速度可以表示为［14]:

对所有频率范围内特征振动速度的平方求和，便得到脉动风速时程的能量总和即方差σ2:

特征速度为v(n)的涡旋从时间t0至t1行进的距离与其波长之比，表征了该时间间隔内涡旋变化的周期数，每个周期对应2π的相位变化，如图2所示。不同尺度、即不同频率涡旋在这一时间间隔τ内的相位改变Δφ(n)可写为:

其中:τ=t1－t0;l(n)为涡旋的波长，表征了涡旋的尺度大小;k(n)为波数，与l(n)互为倒数;波数与自然频率的关系为:

这里U为平均风速。

对式(5)关于时间t求导，便可得到不同频率涡旋的相位演化速度:

图2 涡旋相位变化示意图Fig.2 Schematic diagram of the change of eddy’s phase angle

2.2 相位的零点演化时间

不同频率涡旋振动速度不同，因而相位演化速度也不同。一般来说，低频、大尺度的涡旋相位变化慢，高频、小尺度的涡旋相位变化快。可以设想，真实脉动风速可视为是由一簇具有相同初始相位的涡旋(谐波)经过时间Te的演化后叠加而成的。最简单的初始相位取值是零相位，如图3所示。记录到的某一段风速时程，可以看作是一批具有相同初始零相位的涡旋经过时间Te演化而来。我们称Te为零点演化时间，单位为s。

图3 相位演化时间Te的示意图Fig.3 Schematic diagram of phase evolution time Te

应用实测的风速时程，可以识别出Te的具体数值。本文采用10 Hz采样、10分钟持时的脉动风速样本进行识别，识别原理如下。

首先根据式(1)和(2)求取样本风速时程的Fourier幅值谱和相位谱φ(n)，应用式(7)，可得到不同频率涡旋的相位演化速度Δ·φ(n);以时程开始时刻为起点沿时间轴t向反方向推进(为计算方便，t仍取正值)，用下式求各频率涡旋的相位:

求得相位φ(n，t)后需要换算到［0，2π)的主值区间中，得到相位主值φm(n，t)。若某一时刻所有涡旋相位主值均为零，则对应此时的t值即为所要识别的零点演化时间Te。

考虑理论抽象与真实问题背景之间的差距，本文推荐采取下述零点演化时间识别方法。

设最低频涡旋的初始相位为φ(n1，0)，特征速度为v(n1)，以时程开始时刻为起点向反方向推进，用下式求得一系列时间点T0(T0仍为正值)，使得此时相位主值 φm(n0，T0)=0

这里p为零点个数搜索范围，是自然数的子集，从1开始直到一个计算量允许的数，每个p值对应一个T0。在各个时间点T0，用下式计算前10个频率点(对应频率为1/600 Hz到1/60 Hz)的相位值:

受计算量的限制，搜索范围设置为p≤106。在所搜索范围内，如果T0满足最大偏差条件:

以及加权约束条件

那么就认为该T0值为所要识别的Te。其中φmc(ni，T0)是为了进行识别而转化到(－π，π] 区间的相位谱主值，转化方式如下:

得到Te后，令所有谐波的初相位为零，从该时刻出发，利用公式(7)的相位演化速度，可以重建经过时间Te后的真实时程起点演化相位谱Φ(n，Te):

进而，结合实测的Fourier幅值谱，利用逆Fourier变换即可得到重建后的时程。

理论上，重建时程应当和原时程一致，但是由于采用了松弛识别准则，重建时程和原时程在细节上会有一定出入，但基本信息保持一致。图4为一段重建时程和原时程的比较。可以看到，重建时程与原始时程的波形基本一致，但在高频部分有一定差别。两段时程的相关系数为0.81(完全相同时为1)，可以认为基本达到了重建的目标。

图4 脉动风速重建时程和原始时程的比较Fig.4 Comparison of rebuilt time history and measured time history of fluctuating wind speed

2.3 零点演化时间的统计建模

由上节可知，任意一个风速时程样本都可以看作是由一簇谐波自零相位出发经过时间Te演化而来。对于脉动风速随机过程，Te显然是一个随机变量。若能够统计得到Te的概率分布，结合Fourier幅值谱，便可以得到模型脉动风速时程的样本集合，从而，可以用随机Fourier函数反映脉动风速过程。

2006年，本研究小组于江苏某地建立了国内第一座风场观测台阵，经近5年观测，得到了大量的风速数据［15，16]。采用这一台阵实测的 10 m、20 m、28 m 和 43 m高度处各200条10分钟风速时程，根据上节的识别方法，识别了800条样本脉动风速时程的零点演化时间Te。对800个Te值做统计并用不同概率统计模型进行拟合，发现Gamma分布可以很好的对Te的分布进行描述，如图5(a)所示。

图5 相位演化时间Te的概率分布拟合Fig.5 Probability distribution fitting of evolution time Te

Gamma分布是一个两参数连续分布函数，一般用来模拟随机事件的等候时间，其概率密度函数为［17]:

其中，θ为尺度参数、k为形状参数，两者都必须为正;Γ为Gamma函数。

应用实测数据做Gamma分布的参数估计时，形状参数 k可用下式估算［17]:

其中:

即s为实测变量均值的对数值与实测变量对数值的均值之差。Gamma分布的均值为kθ，因此容易得到尺度参数θ的计算公式:

根据前述800个样本统计，得到的模型参数值为:k≈1.1，θ≈8.2 ×108。进而，应用其他1 200 条 10 分钟风速记录对该分布参数进行了检验，所得效果仍然非常理想，见图5(b)。

3 脉动风速模拟

3.1 随机Fourier函数

前已指出:谱表现法是基于功率谱密度的脉动风速模拟方法。虽然许多学者提出了不同的脉动风速功率谱密度模型［18]，但从本质上来说，由于功率谱密度是二阶矩，谱表现方法仅适用于平稳随机过程［7]。针对谱表现方法的固有局限性，我们提出应采用随机函数模型进行风速时程模拟的观点［9]。随机Fourier函数的基本表达式如下式:

式中，η和γ为影响脉动风速随机过程的可测物理因素。γ为确定性变量，η为随机变量。当需要计算具体样本时，可取η的具体实现值。

前已述及，Fourier谱可以分解为幅值谱和相位谱，其中Fourier幅值谱是对脉动风速能量的一种谱分解，无论对于平稳随机过程或是非平稳随机过程都适用［7]。根据大气湍流的物理图景和经典湍流理论，我们已经建立了随机Fourier幅值谱的双线性模型［9]，其中基本的随机变量为地面粗糙度z0和10分钟平均风速。本文的脉动风速时程模拟中的幅值谱即基于该双线性模型，模型的具体细节请参阅文献［9] 。

考虑到本文重点在于Fourier相位谱，因此进行脉动风速模拟时可给定幅值谱，如图6所示。其中，在剪切含能区满足“－1”幂次规律，惯性子区满足“－5/3”幂次规律，两个子区交点在频率0.12 Hz。平均风速取10 m/s、湍流强度取 0.15。

图6 Fourier幅值谱双线性模型Fig.6 Bilinear model of Fourier amplitude spectrum

3.2 演化相位谱检验

由于零点演化时间Te的概率模型是由800个时程样本统计得到，为进行集合样本检验，也等概率选取800个Te样本值，分别生成演化相位谱。这里等概率取样的含义是:等间距分割概率分布函数，并取分割位置概率分布函数所对应的样本值。图7为Gamma分布的概率分布函数与等概率取值示意。

图7 Gamma分布概率分布函数与等概率取值示意Fig.7 CDF of Gamma distribution and schematic diagram of equally probability sampling

对所选取的800个Te的样本值，应用式(14)生成演化相位谱并换算到［0)的主值区间。图8比较了不同频率位置实测的800条样本相位谱和所生成的800条演化相位谱主值的概率分布直方图，图中直线为均匀分布的概率密度函数。可以看到，在各个频率点处，实测和演化所得相位谱主值均基本符合［0，2π)内均匀分布。这与已有相关研究结果［10]是相符的。

3.3 脉动风速模拟

图8 不同频率处相位谱主值概率分布直方图比较Fig.8 Histogram of the distribution of phase-spectrum main value at different frequency

其中，Re表示取实部，Fs为采样频率，d n是频率间隔，与采样持时T互为倒数。除以是将能量调整至与原单边谱相同。

图9为任意选取的四个Te值计算出的演化相位谱主值及其分布，图10为用其生成的脉动风速时程。两图从上到下，Te均依次取值为 1.131e8 s、4.036e8 s、7.129e8 s、12.997e8 s。可见，不仅相位谱主值具有［0，2π)内均匀分布的特征，所生成的脉动风速时程也与常见的风速观测结果有很好的相似性。

图9 不同Te所生成的相位谱主值及其分布Fig.9 Main value of phase spectrum generated by different Te and its distribution

3.4 模拟风速的检验

一般认为，良态风主风向脉动速度的概率分布可用正态分布近似描述［19]。比较两个概率密度函数相近程度的一个有效指标是相对熵［20]。实际风速样本的概率密度p(u)和其对应正态分布概率密度函数f(u)的相对熵为:

图10 不同Te所生成脉动风速时程Fig.10 Fluctuating wind speed history generated by different Te

实际风速样本的概率密度p(u)用直方图代替。由于每个时间截口脉动速度u的均值为零，因此正态分布的概率密度函数为:

其中σ为截口样本的标准差。当两个概率密度函数完全相同时，相对熵为零。当两个概率密度函数相差越大时，相对熵绝对值越大。

计算得到的相对熵如图11所示。可以看到，相对熵的值较小但不为零，说明各时间截口的脉动风速分布近似符合正态分布但又略有区别，这也与已有实测结果相吻合［19]。由此可以判断，本文提出的脉动风速模拟方法是合理的。

图11 脉动风速样本概率密度与正态分布的相对熵Fig.11 Relative entropy of sampled fluctuating wind speed and normal distribution

4 结论

依据分析，本文提出任何一段脉动风速时程都可以认为是由一簇具有零相位的谐波经过时间Te演化而来，根据800条实测风速记录统计得到了Te的概率模型。这一概率模型经1 200条实测风速记录检验正确。结合Fourier幅值谱，由Te的样本值可以得到相位谱的样本，进而进行脉动风速模拟。

与传统方法相比，相位零点演化时间Te的概念将原本无规律的各谐波相位联系在一起，且具有物理意义。演化相位谱的变量仅为Te，不仅大幅度减少了变量数目，也使得采用精细化的概率密度分析方法进行结构抗风可靠度的计算成为可能。

［1] Kareem A. Numerical simulation of wind effects:A probabilistic perspective［J] .Journal of Wind Engineering and Industrial Aerodynamics，2008，96:1472－1497.

［2] Rice S O.Mathematical analysis of random noise［C] .In Selected Papers on Noise and Stochastic Processes，edited by N.Wax.Dover.1954:133－294.

［3] Shinozuka M.Simulation of multivariate and multidimensional random process［J] .Journal of the Acoustical Society of America，1971，49(1):357－367.

［4] Shinozuka M，Jan C M.Digital simulation of random process and its application ［J] .Journal of Sound and Vibration，1972，25(1):111－128.

［5] Shinozuka M，Deodatis G.Simulation of stochastic processes and fields［J] .Probabilistic Engineering Mechanics，1997，12(4):203－207.

［6] Shinozuka M，Yun C B，Seya H.Stochastic methods in wind engineering［J] .Journal of Wind Engineering and Industrial Aerodynamics，1990，36:829 －843.

［7] Li J，Chen JB.Stochastic dynamics of structures［M] .John Wiley＆ Sons(Asia)Pte Ltd.2009.

［8] 李 杰.随机动力系统的物理逼近［J] .中国科技论文在线，2006，1(2):95 －104.

［9] 李 杰，阎 启.结构随机动力激励的物理模型——以脉动风速为例［J] .工程力学，2009，26 Sup II，175－183.

［10] Seong SH，Peterka J A.Experiments on Fourier phases for synthesis of non-Gaussian spikes in turbulence time series［J] . Journal of Wind Engineering and Industrial Aerodynamics，2001，89:421 －443.

［11] 金 星，廖振鹏.地震动相位特性的研究［J] .地震工程与工程振动，1993，13(1):7－13.

［12] Taylor G I.The spectrum of turbulence［J] .Proceedings of the Royal Society of London，Series A，1938，164:476.

［13] Bahraminasab A，Niry M D，Davoudi J，et al.Taylor's frozen-flow hypothesis in Burgers turbulence［J] .Physical Review，2008，E77.065302(R).

［14] Hinze J Q. Turbulence［M] . New York:McGraw-Hill，1975.

［15] 阎 启，谢 强，李 杰.风场长期观测与数据分析［J] .建筑科学与工程学报，2009，26(1):37－42.

［16] 李 杰，阎 启，谢 强，等.台风“韦帕”风场实测及风致输电塔振动响应［J] .建筑科学与工程学报，2009，26(2):1－8.

［17] Choi S C，Wette R.Maximum likelihood estimation of the parameters of the Gamma distribution and their bias［J] .Technometrics，1969，11(4):683 －690.

［18] Dyrbye C，Hansen SO.Wind loads on structures［M] .New York:John Wiley＆ Sons，1997.

［19] Yim J Z，Chou C R，Huang W P.A study on the distributions of the measured fluctuating wind velocity components［J] . Atmospheric Environment，2000，34:1583－1590.

［20] Sobezyk K， Trebicki J. Maximum entropy principle in stochastic dynamics［J] . Probabilistic Engineering Mechanics，1990，5(3):102 －110.