平稳随机脉动风场模拟的高维概率空间选点法

徐军 李哲帆 张洋 李扬

(湖南大学 土木工程学院，湖南 长沙 410082)

现有结构设计规范通常忽视风荷载的动力特征，而将其等效为静力荷载。对于高层建筑、大跨网壳结构、大跨度桥梁结构等风敏感性较强的结构，风荷载的动力效应往往不可忽视。风荷载下结构动力响应最主要来源于脉动风速场的作用[1]，而脉动风速场具有显著的时空分布特征和随机性，通常被视为在不同空间位置相关的平稳随机向量过程，或称为时-空随机场[2]。合理地描述和模拟随机脉动风速场的时域动力特征，对于高层建筑和大跨结构等复杂工程结构的抗风分析和保障其风灾安全性具有重要意义。

目前，模拟平稳随机向量过程的主要方法有：谱表达法[3]、线性滤波法、本征正交分解法[4]以及小波变换方法[5]等。其中，由于谱表达法仅涉及一系列谐波叠加，其算法相对简单且具有严谨的理论基础，在工程实际中被广泛采用。谱表达法的起源可追溯到20世纪40年代对一维单变量白噪声的描述[6]。此后，Shinozuka[7]将其引入模拟一维多变量平稳随机向量过程；Yang[8]通过引入FFT算法进一步提高了平稳随机向量过程的模拟效率；Shinozuka等[9]从高斯性、无偏性和遍历性的角度进一步详细阐述谱表达法的原理[9]；为生成遍历性更好的样本函数，Deodatis等[10]引入了双索引频率，使得谱表达法成为模拟随机向量过程的经典方法，被广泛地应用到风场、地震、波浪荷载的模拟中。

平稳随机脉动风速场模拟的谱表达法中往往涉及数目庞大的随机变量，通常采用蒙特卡洛模拟的方式产生基本随机变量样本点来得到多个空间点上的脉动风速过程。一方面，为精确地捕捉平稳随机脉动风速场的统计特性，所需的蒙特卡洛模拟样本数目往往较大，造成后续结构随机风振响应和可靠度分析的计算量难以令人接受；另一方面，蒙特卡洛模拟结果本质上并非完备的概率集合，因此不能在概率密度层次精细地描述风场的概率特性[11]。为克服此问题，Chen等[12]建议用一类随机谐波函数的谱表达法对随机风场进行模拟，仅需要少量的随机变量便可较好地模拟平稳随机脉动风速场特性；李杰等[13]从物理随机系统角度出发，提出了随机脉动风速场的物理模型，将若干关键物理参数描述为随机变量；刘章军等[14]提出基于随机函数的降维模拟方法，结合谱表达方法仅用1-2个基本随机变量即可在二阶统计意义上对平稳随机脉动风速场进行精确描述。在上述研究基础上，本文拟提出一类高维概率空间选点法，实现基于谱表达法的平稳脉动风速场精确与高效模拟，为工程结构随机风振响应与可靠度分析奠定基础。

1 平稳随机脉动风速场的谱表达

假设X(t)=[X1(t),X2(t),…,Xn(t)]T是nV-1D零均值平稳随机向量过程，其双边功率谱密度函数矩阵SX(ω)可表示为

(1)

式中，n为空间点数目，ω为圆频率，Sii(ω)为自谱密度函数，Sij(ω)为空间点i处的风速分量和空间点j处的风速分量互谱密度函数。一般地，SX(ω)为正定Hermitian矩阵，满足以下关系：

Sii(ω)=Sii(-ω)，i=1,2,…,n

(2)

(3)

(4)

为避免直接对谱密度函数矩阵进行乔列斯基分解而造成存储困难,文献[15]中引入如下空间相干函数矩阵：

(5)

β(ω)=B(ω)B*T(ω)

(6)

那么，互谱密度函数Sij(ω)可进一步表示为

(7)

根据式(2)-(4)及式(7),功率谱密度函数矩SX(ω)可进一步分解为[16]：

SX(ω)=K(ω)B(ω)B*T(ω)K*T(ω)

(8)

式中,K(ω)为如下对角矩阵,

(9)

由Fourier-Stieltjes积分可知：

(10)

式中，Z(ω)为增量正交零均值复向量过程，且满足

E[dZ(ω)]=0n×1,dZ(-ω)=dZ*(ω)

(11)

E[dZ(ω)dZ*T(ω′)]=SX(ω)dωdω′δ(ω-ω′)

(12)

式中，E表示期望算子，“*”是共轭符号，δ为狄拉克函数。

那么,第i个复过程增量dZi(ω),i=1,2,…,n,表示为

(13)

式中，Δω=ωu/N为频率间隔，ωu为频率上限，N为频率的截断项数；Pgh=Ugh-iVgh表示零均值标准化的复随机变量，Ugh和Vgh为实正交随机变量，且满足

E[Ugh]=E[Vgh]=0,E[UghVrs]=0;

E[UghUrs]=E[VghVrs]=(δgrδhs)/2;

g,r=1,2,…,n;h,s=1,2,…,N

(14)

为各态历经性，选用了如下双索引频率：

(15)

令Big(ωgh)=ξig(ωgh)-iΨig(ωgh)，由式(15)可知-ωgh=ωg(-h)，式(10)可进一步表示为

{ξig(ωgh)(Ughcosωght+Vghsinωght)+

Ψig(ωgh)(Ughsinωght-Vghcosωght)}

(16)

(i=1,2,…,n)

空间相干矩阵β(ω)一般是实对称矩阵，其下三角矩阵β(ω)亦为实值矩阵，式(16)可进一步简化为

[Ughcosωght+Vghsinωght]

(17)

由此可见：平稳脉动风速随机向量过程的模拟即可转化为一系列正交随机变量调制的谐波函数叠加，进而可获得多个空间点处的脉动风速时程序列。一般来说，式(16)和(17)中涉及的正交随机变量(即Ugh和Vgh)一般取为标准正态分布随机变量。由此可见：平稳脉动风速随机向量过程的谱表达中涉及成千上万个标准正态分布随机变量(2×n×N个)，为典型高维随机问题。传统的蒙特卡洛模拟所需样本数量往往过大，模拟效率低下，亦会对后续结构随机动力响应分析效率造成严重负担。为结合概率密度演化方法进行结构随机风振响应与可靠度分析，需产生脉动风速场典型代表样本及其赋得概率，有必要对高维正态空间离散代表点的选取问题展开专门研究，以期产生少量的代表点数目，构成一个完备的概率集，提高后续结构随机动力响应的计算效率。

2 高维正态空间离散代表点选取法

借鉴部分分层抽样的思想[17]，本文提出一类新的高维正态空间离散代表点选取方法来解决这一问题，由如下3步构成：

步骤1 首先，将高维标准正态空间剖分为一系列2维正交子空间；

步骤2 其次，在每个正交子空间中应用现有低维概率空间选点方法进行子空间离散代表点选取和赋得概率的计算；

步骤3 最后，将子空间中的离散代表点和赋得概率通过随机映射的方式获得原始高维标准正态空间的离散代表点及其赋得概率。将上述高维标准正态空间离散代表点及其赋得概率代入至式(17)中，则可生成代表性风速过程样本进行有限元计算。

2.1 高维概率空间分解

由式(17)可见：平稳随机脉动风速场谱表达模拟中共包括2Q(Q=n×N)个独立标准正态随机变量(即Ugh和Vgh)，涉及高维概率空间问题。鉴于直接处理高维随机问题的困难，本文首先将原始高维正态空间分解为一系列正交的低维子空间。根据文献[18]中对正交子空间的定义，可以将高维正态空间Ωr划分为Q=n×N个二维正交子空间，即：

(18)

式中，Ωi,2表示第i个2维标准正态子空间，例如，1个10维标准空间可以剖分为5个2维正交子空间。

2.2 正交子空间中离散代表点的选取

在获得一系列2维正交子空间后，对正交子空间采用数论选点方法确定其离散代表点及赋得概率。首先，采用华罗庚-王元的数论方法获得二维子空间上的均匀分布点集uk,k=1,2,…，nsel[19]:

(19)

(k=1,…，Fm)

式中，uk=[u1,k,u2,k],Fm=Fm-1+Fm-2,m=2,3,…;F0=F1=1为斐波那契数列，int(·)表示取数值的整数。本文取m=13,点集总数目nsel=F13=377，二维均匀分布点集如图1所示。

图1 二维均匀分布点集Fig.1 Uniform points in the 2-dimensional unit square

(20)

(j=1,2;k=1,2,…,nsel)

式中，Φ-1指的是标准正态分布的逆累积分布函数。

图2所示为图1中点集转化至标准正态空间后所得到的子空间初始离散代表点结果。

图2 二维标准正态空间点集Fig.2 Points in the 2-dimensional standard normal space

(21)

(22)

图3所示为子空间中离散代表点所对应的Voronoi域剖分，计算所得的子空间离散代表点赋得概率如图4所示。

图3 子空间中Voronoi域剖分Fig.3 Voronoi dissection in subspaces

图4 子空间离散代表点的赋得概率Fig.4 Assigned probabilities for points set in subspaces

(23)

(j=1,2)

若A成立，I{A}=1；否则I{A}=0。图5所示为初始离散代表点和重整后离散代表点的对比图，可见：重整后点集的位置与初始点集位置略有不同。

图5 初始点集与重整后点集对比Fig.5 Comparison between the initial points set and the one after rearrangement

2.3 高维离散代表点及赋得概率

(24)

图6 点集随机配对Fig.6 Random permutation of the point set

由于低维子空间均正交，可认为高维离散代表点的赋得概率为上述随机配对的低维正交子空间离散代表点赋得概率的乘积：

(25)

进一步,由归一化条件得到高维离散代表点赋得概率为:

(26)

3 框架-剪力墙结构平稳随机脉动风场模拟

文中采用某三十层框架-剪力墙结构的平稳随机脉动风场来验证本文提出的方法的正确性。该框架-剪力墙结构立面图和平面图如图7所示(黑色区域代表剪力墙和框架柱)。模型底部柱高为6 m，标准层高度为3 m，迎风宽度15 m，坐标系xoz与风向平行。将该结构上的随机脉动风场近似为45个空间点处的随机脉动风速矢量过程。其空间分布如图7所示。表1所示为该模型中用于ANSYS建模的单元、尺寸参数和材料参数信息。

表1 模型相关参数Table 1 Related parameters of the model

图7 结构模型立面图与平面图Fig.7 Elevation and plan layout of structural model

本例中，采用Davenport谱作为双边自功率谱密度函数，定义如下[21]：

(27)

(i=1,2,…,n)

采用本文所建议的高维标准正态空间选点方法，可以模拟结构上45个空间点(n=45)的脉动风速矢量代表过程。表2所示为模拟随机脉动风速矢量过程所需的参数。

表2 随机脉动风速场模拟所需相关参数Table 2 Related parameters required to simulate a random fluctuating wind speed vector

选择高度为6、42、78 m的空间点1、19和37的仿真结果来说明本文方法的正确性和效率。图8所示为3个空间点的平稳风速时程代表性样本过程，可以清楚地看到：3个风速时程过程样本的值域范围大致相同，且样本在该值域范围内波动，呈现平稳特性。

图8 风速时程曲线代表性样本Fig.8 Representative samples of fluctuating wind speeds

由国内外相关文献可知:一般选择均值误差、标准差误差、功率谱密度函数平均值误差3项指标来考察平稳风速过程的统计特性,从而进一步验证本文方法的准确性。那么，第q点模拟样本过程的均值、标准差和功率谱密度函数平均值可表示为

(28)

(29)

(30)

(31)

(32)

图9 均值模拟值与目标值的比较Fig.9 Comparison between simulated mean values and target values

图10 标准差模拟值与目标值比较Fig.10 Comparison between simulated standard deviation values and target values

图11 功率谱密度函数平均值与目标值比较Fig.11 Comparison between the mean of simulated PSD values and target values

从图中不难发现：采用本文方法模拟的平稳随机脉动风场的均值、标准差及功率谱密度函数平均值与目标值均吻合良好，表明本文方法具有较高的精度。此外，由于本文方法仅采用带有赋得概率的377个高维标准正态空间离散代表点来模拟平稳随机脉动风场，亦表明本文方法具有较高的计算效率。进一步地，采用平均相对误差进行定量分析来表明本文方法的准确性。那么，第q个空间点脉动风速随机过程均值、标准差和功率谱密度平均值的平均相对误差可定义为[22]：

(33)

(34)

(35)

式中，Nt表示时间步数，Nω表示频率步数。同时，采用相同样本数目的蒙特卡洛模拟和文献[3]中的降维谱表达法，与本文方法进行对比，结果见表3。

表3 3种方法平均相对误差比较Table 3 Comparison of average relative errors among three methods %

从表中可以发现：相较于这两种方法，本文方法计算得到的均值、标准差及功率谱密度平均值的精度均较高，相对误差更小。特别地，均值的精度可以提高1-2个量级，这意味着：可利用较少数目样本模拟随机脉动风场，满足精度要求，从而进一步降低后续结构随机风振响应与可靠度分析的计算量。

将上述平稳随机脉动风场导入至通用有限元软件(ANSYS)，对结构模型进行确定性风振响应分析，可以获得每个风速矢量样本相对应的结构动力响应。考察结构的层间位移角

(36)

式中，xl(t)表示在第l层在x方向上的水平位移，hl表示第l层的高度。那么，可方便地计算得到层间位移角的均值与标准差。图12所示为顶层层间位移角的均值与标准差。进一步地，将这些带有概率信息的响应过程代入至概率密度演化方法中，即可方便地获得结构随机风振响应概率密度演化全过程及动力可靠度信息。

图12 顶层层间位移角响应的均值和标准差Fig.12 Mean and standard deviation of the top inter-storey drift

考虑到结构形式的多样化，笔者还选取了K8单层球面网壳进一步验证了本文方法的有效性，限于篇幅，本文不再赘述。

4 结语

采用谱表达方法模拟平稳随机脉动风场随机向量过程往往涉及高维标准正态随机变量，对数值模拟的精度和效率提出挑战。鉴于此，本文借鉴部分分层抽样思想，提出了一种新的高维正态空间选点方法,这一方法采用了Voronoi域剖分、点集重整和随机配对等技术手段获得高维正态空间离散代表点及其赋得概率，将其代入至平稳随机脉动风场谱表达模型中，即可得到带有完备概率信息的风场代表性样本过程。

采用在30层框架-剪力墙结构的平稳随机风速随机向量过程模拟验证了本文方法的准确性与高效性，研究结果表明：本文方法能够准确地获得平稳脉动风速随机向量过程的统计信息。值得注意的是：相较于蒙特卡洛模拟方法，本文方法能用较少的样本数构成完备的概率集，可大幅减少后续的结构随机风振响应与可靠度分析所需的有限元分析计算量，展现了这一方法在工程实际运用中的优越性。

若考虑随机脉动风场的非高斯性与非平稳性，拟进一步拓展这一方法的适用性。