大跨度桥梁脉动风场模拟的插值算法

祝志文， 黄 炎

(1.湖南大学 土木工程学院，长沙 410082；2. 湖南大学 风工程与桥梁工程湖南省重点实验室，长沙 410082)

大跨度桥梁脉动风场模拟的插值算法

祝志文1，2， 黄 炎1

(1.湖南大学 土木工程学院，长沙 410082；2. 湖南大学 风工程与桥梁工程湖南省重点实验室，长沙 410082)

为减少谐波合成法中功率谱矩阵分解的计算量，基于双索引频率和谱矩阵特点，提出了谱解矩阵双轴插值算法和递归插值算法。两种方法均使风场模拟计算效率大幅度提高，且得到的模拟风场均与传统谐波合成法在统计意义上相符；研究了插值节点数量及分布形式对风场模拟的影响，认为在谱解矩阵的频率轴向采取前密后疏分段插值，双索引频率轴向采用均匀分布插值节点的形式，更适合于风场模拟。通过对一实际大跨度桥梁风场的模拟，验证了方法的有效性。

风场模拟；大跨度桥梁；谐波合成法；插值算法

在时域内进行大跨度桥梁抖振响应分析，首先需要在脉动风场离散空间点上，生成所需的具有特定要求的脉动风速时程。目前，脉动风场数值模拟方法主要有基于三角级数的谐波合成法和基于线性滤波技术的回归方法等[1]。谐波合成法以功率谱作为权系数，与一系列带随机相位的三角级数的加权和来逐渐逼近随机过程，适用于指定谱特征的平稳高斯随机过程。该方法有恒幅谐波叠加法和加权振幅谐波叠加法[2]。加权振幅谐波叠加利用FFT(Fast Fourier Transform)算法可以明显减小机时，但由于频率均匀分布导致模拟曲线出现周期性。为此，SHINOZUKA等[3]引入双索引频率，将频率微扰量均匀分布在频率增量内，这样既可采用FFT算法，又可保证模拟曲线的各态历经性。CAO等针对大跨度桥梁风场离散点沿桥向可均匀分布的特点，导出了谱分解的显示表达式，减少了部分余弦项，提高了计算效率，但该方法不能考虑三维空间风场三个方向之间的相关性。DING等提出谱分解的三次拉格朗日多项式插值近似建议，但未给出插值节点的具体布置形式。BAO等[4]引入BP神经网络方法拟合谱分解函数曲线，减少了谱分解的计算量。XU等[5]根据实际大跨度桥梁风速监测点少的特点，采用克里金插值，模拟了桥梁抖振分析所需其他未测风场离散点的非平稳脉动风。

本文基于双索引频率和谱矩阵特点，提出了谱解矩阵双轴插值算法和递归插值算法。通过对一实际大跨度桥梁风场的模拟，验证了方法的合理性。

1 脉动风场谐波合成理论

1.1 传统谐波合成方法

根据随机过程理论[6-9]，模拟所生成的样本必须能准确描述随机过程的概率特性。对于平稳随机过程，脉动风速的特性可用功率谱和相关函数予以描述，功率谱密度函数揭示了随机振动过程的能量按频率分布的规律，相关函数则反映各点脉动风速之间在时间或空间的相互关系。三维空间n个风场离散点的功率谱矩阵见式(1)。

(1)

在工程实际中，三维空间点之间的互谱密度通常用实函数表达，因而互谱矩阵S(ω)为实对称矩阵。Sjk(ω)是互谱矩阵S(ω)的元素，其中ω为风的脉动圆频率，j,k=1,2,…,n分别表示所模拟的第j个和第k个风场离散空间点。当j=k时为脉动风速的自谱密度，当j≠k时为脉动风速的互谱密度。互谱密度可由自谱密度和三维空间相干函数根据式(2)确定

(2)

其中，脉动风速的三维空间相干函数可表示为

Cohjk(ω)=

(3)

对互谱矩阵S(ω)进行Cholesky分解，即S(ω)=H(ω)H(ω)T，H(ω)为下三角矩阵，也是实数矩阵，H(ω)T是其转置矩阵。

(4)

(5)

(6)

互谱矩阵S(ω)和谱分解矩阵H(ω)都是双索引频率的矩阵，若用四维矩阵表示则为一个n×n×N×n的四维矩阵，其中，第三维对应频率轴，第四维对应双索引频率轴。显然，随着模拟点n的增多，这个四维矩阵呈幂级数增大，不仅存储量会非常巨大，而且要进行n×N次n×n阶矩阵的Cholesky分解，计算量相当大，甚至很难实现。

1.2 谱解法的插值方法

1.3 谱解法的递归算法

为了提高谱分解效率，将风场模拟式(5)中的|Hjm(ωml)|项写为矩阵形式

(7)

式(7)的第j行元素就是第j个脉动风场离散点模拟所需要的j个谱矩阵分解元素。该矩阵第m列元素仅取式(4)第m个双索引频率对应谱矩阵分解后的第m列。也就是说，对于每次功率谱矩阵分解，仅有一列为风场模拟所需，其余n-1列元素为多余元素。为了减少不必要的谱分解，根据递归算法，对第m个双索引频率对应的互谱矩阵[S(ωml)]进行式(8)形式分块，使其分块后的子矩阵[S11]为m×m阶方阵。

(8)

2 谱解法的应用

在风工程领域中，许多学者通过对自然风进行长期的观测，基于大量的实测数据，建立了不同的风谱模型，以反映自然风脉动的频率结构。本文采用我国规范《公路桥梁抗风设计指南》推荐的Simiu谱模拟顺风向脉动风。

(9)

本文针对一座实际大跨度桥梁开展了风场模拟。该桥中跨主缆跨度1 176 m，主缆中心距27 m；中跨加劲梁跨度1 005 m，加劲梁梁宽27.6 m、梁高8.65 m。该桥采用钢桁加劲梁，吊杆顺桥向中心距为14.5 m，相邻吊杆间设有两个桁架节段。两侧主塔分别高134 m、67 m。

传统的大跨度桥梁抖振响应分析，脉动风荷载作用往往仅考虑主要的风荷载作用构件加劲梁，忽略了主缆和桥塔上作用的风荷载，这样可能会影响大跨度桥梁抖振响应分析的精度，但由于仅在加劲梁上布设脉动风场模拟的空间离散点，因而脉动风场模拟的计算量明显减少，从而对计算资源的要求降低。本文从考虑大跨度桥梁全部主要风荷载作用构件的角度出发，风场离散空间点按照如下方式布置：沿顺桥向中心线布置在主桁上弦节点处，其编号为1～137；两侧主塔分别沿高度均匀布置15个离散点，其编号为138～152；上下游主缆风场离散点沿主缆布置，其编号为153～374。全桥共374个风场离散点，其空间分布见图1，其中，括号标注为横桥向对称分布于另一侧主缆上模拟点编号，这样便于获得风荷载作用的桥梁全部主要构件不同空间点的模拟风场，目的是为开展全桥全部主要构件在风荷载作用下响应分析。风场参数见表1。

表1 脉动风场模拟参数

图1 某悬索桥三维风场模拟点布置(单位：m)Fig.1 Spacial arrangement of wind field points on a suspension bridge (Unit: m)

为探讨功率谱矩阵在频率和双索引频率轴向插值节点的选取形式对谱分解矩阵插值后的影响，进而影响到风速样本模拟。本文在谱分解矩阵的频率和双索引频率轴向分别单独进行插值，各取5种工况(见表2)分别计算H(ω)与传统方法谱解矩阵对比，其中，在频率轴向为1 024个频率点，在双索引频率轴向为374个频率点，插值步距定义为插值节点间的频率点数量。由此为合理选择在频率轴向以及双索引频率轴上的插值节点布置形式提供依据。

表2 双轴插值10种工况

2.1 插值节点分布形式的讨论

2.1.1 频率轴向插值

2.1.2 双索引频率轴向插值

由图3可知，双索引频率轴的H(ω)值基本呈线性变化，这是由于双索引频率的本质是为了延长模拟样本周期性和保证各态历经性，而在频率增量内为均匀分布的微小扰动。5种插值工况与传统方法的谱分解矩阵均保持一致，所以在双索引频率轴，插值节点可以选用均匀分布形式，插值步距可以选用较大间隔。

2.2 合理选点的双轴插值法

图2 频率轴向插值结果Fig.2 Frequency axis interpolation

图3 双索引频率轴向插值结果Fig.3 Double-indexing frequency axis interpolation

图4 双轴插值法频率轴对比Fig.4 Frequency direction results

图5 双轴插值法双索引频率轴对比Fig.5 Double-indexing frequency direction results

如图6所示，双轴插值法与传统方法模拟得到的风速样本整体统计特性较为吻合，所存在的差异主要是由于谱分解矩阵的高阶项在高频频带部分趋近于0，从空间意义上说，风场空间相关性的高频频带部分较低频频带部分要弱，相距越远的离散点脉动风速相关性越小，而插值造成了在高频部分的震荡偏差。但通过后面的功率谱、相关函数对比，知其仍可以接受，也可进一步加密插值节点更趋近于传统方法。

图6 风场统计值Fig.6 Statistic of wind field

进一步从整体统计特性上，对比模拟风场离散点的相关系数矩阵，其为对称矩阵，主对角线表示各点自相关系数为1，为方便对比，将上三角取为传统方法，下三角取为双轴插值法(见图7)。由图7(a)整体对比可见两者对称性很好，且与式(3)表示的任意频率点上各点相干函数矩阵波峰位置相契合。进一步取主梁(见图7(b))、主缆(见图7(c))以及一侧主缆与主梁(见图7(d))进行详细对比，在工程中一般认为相关系数在0.1～0.2以下则两者基本为不相关。从图7(b)可知主梁某离散点相邻约20点外的相关系数下降到了0.1以下，即可认为主梁相距超过约145 m则两点相关性很小；从图7(c)可知，主缆存在两个明显的波峰，中间波峰为自相关系数，另一个波峰为两侧主缆相邻点所对应的互相关系数；由图7(d)可知除开自相关系数的波峰，另一个小的波峰中间高，两侧低，其正好对应主梁与边跨、中跨主缆的互相关系数。综上可见，双轴插值法的谱分解矩阵及模拟风场样本在整体统计特性上均与传统方法基本保持一致。

图7 风场相关系数对比Fig.7 Correlation coefficient of wind field

为说明双轴插值谐波合成法风场模拟效果，以加劲梁所在平面为参考基准面，分别给出了加劲梁跨中第65点、第66点，右侧索塔103 m高处第150点，以及主缆跨中12 m高处第211点(布置见图1)的风速时程，虽然在时域内风速样本的统计特性无法看出，但显然风速时程在相邻的主梁第65点(见图8(a))、第66点(见图8(b))及主缆离散点第211点(见图8(d))的变化趋势较主塔第150点(见图8(c))更为相似。进一步在频域内对比模拟功率谱函数与目标功率谱函数，为包络住大跨度桥梁风场卓越频率，给出0.01～2 Hz频带的功率谱函数，并采用双对数坐标表达形式，其中模拟样本互功率谱为复数形式，图中幅值为其实部与虚部的模，由图9可知，模拟样本功率谱与目标功率谱在整个频带均较为一致。由于两点的相关性随距离增大而逐渐衰减，相关函数仅给出临近的互相关函数，同时，脉动风速经过一段时距后，也会与起始状态的相关性逐渐衰减，所以仅给出时滞在200 s以内的相关函数对比结果，由图10可知，模拟相关函数与目标函数变化趋势基本一致，均随着时间增加，相关性逐渐降低，约20 s后脉动风速样本与起始状态基本没有了相关性，即脉动风速流经约400 m后基本不相关，这与自然风积分尺度在150～450 m相吻合。由此说明了双轴插值方法模拟的风速样本的可靠性。

图8 风速时程曲线Fig.8 Wind speed samples

图9 自/互功率谱密度函数Fig.9 Auto-/cross-spectral density

2.3 风场模拟的递归插值法

所谓递归插值法，是根据“1.3”所述谱分解矩阵的特点，在双索引频率轴向采用递归算法。根据“2.2”所述，在频率轴向采用分段插值方法。也就是说，在所选的频率插值节点上，样本合成所需的谱分解矩阵与传统方法完全一致，而频率插值节点布置形式根据“2.1.1”讨论采用分段插值，插值后的谱分解矩阵与传统方法近似相同。由图6的均值和标准差对比可知，递归插值法模拟风场效果与双轴插值方法相同，这是因为在频率轴向，两种方法的插值节点布置相同，而在双索引频率轴向，谱分解矩阵值呈线性分布。也进一步说明，在高阶项高频部分出现的振荡偏差均由频率轴向插值导致，可减小频率轴向插值步距消除。由于递归插值法与双轴插值法的模拟效果相同，仅给出递归插值法的模拟功率谱与目标谱对比结果，由图11可知，递归插值法模拟风速样本功率谱与目标谱在整个频带上都吻合较好。

图10 自/互相关函数Fig.10 Auto-/cross-correlation function

图11 自/互功率谱密度函数Fig.11 Auto-/cross-spectral density

3 模拟时间对比

本文在CPU为3.40 GHz的PC机上，利用MATLAB编制了谐波合成传统方法、双轴插值法和递归插值法这三种方法，模拟了图1所示大跨度悬索桥主梁、主塔、缆索共374个风场离散点，其模拟时间对比见表3，双轴插值法的相对较短，且从存储量来看，谱矩阵在频率轴上若采用双精度浮点类型存储，传统方法为374×374×1 024的矩阵，需占用近800 MB内存，而双轴插值法为374×374×183的矩阵，只需150 MB内存。同时，由第“2”节讨论可知，双轴插值法以及递归插值法在风场模拟统计意义上，与传统方法基本保持一致，综上所述，本文建议采用双轴插值法。

表3 三种方法的模拟时间对比

Tab.3 Simulation time of the three methods min

谐波合成法传统方法双轴插值递归插值模拟时间3391169

4 结 论

(1) 本文采用双轴插值法以及递归插值法对三维空间大跨度桥梁脉动风场进行了模拟，减少了谱分解次数，提高了谱分解效率。对比传统方法模拟的风速样本均值、标准差、相关系数矩阵、功率谱密度函数、相关函数，验证了这两种方法的可靠性。

(2) 讨论了插值节点分布形式对谱分解矩阵的影响，得出了在频率轴采用前密后疏的分段插值，在双参数频率轴向采用较大步距的等间距插值，得到的谱分解矩阵与传统方法较吻合。

(3) 谱分解矩阵值在双索引频率轴呈线性分布，若双轴插值法和递归插值法在频率轴向的插值节点布置相同，那么，两者模拟效果也相似。在谱分解矩阵高阶项高频部分出现的振荡误差均来自于频率轴向插值，可减小插值步距消除。

[ 1 ] CAO Y, XIANG H, ZHOU Y. Simulation of stochastic wind velocity field on long-span bridges[J]. Journal of Engineering Mechanics，2000,126(1):1-6.

[ 2 ] DING Q, ZHU L, XIANG H. An efficient ergodic simulation of multivariate stochastic processes with spectral representation[J]. Probabilistic Engineering Mechanics, 2011,26(2):350-356.

[ 3 ] LIANG J, CHAUDHURI S, SHINOZUKA M. Simulation of nonstationary stochastic processes by spectral representation[J]. Journal of Engineering Mechanics, 2007,133(6):616-627.

[ 4 ] 包龙生，刘克同，于玲,等. 大跨度桥梁空间脉动风场的数值模拟[J]. 沈阳建筑大学(自然科学版)，2010，26(2)：238-243. BAO Longsheng, LIU Ketong, YU Ling, et al , Numerical simulation of spatial fluctuating wind field on long-span bridges[J]. Journal of Shenyang Jianzhu University (Natural Sciences), 2010,26(2):238-243.

[ 5 ] XU Y, HU L, KAREEM A. Conditional simulation of nonstationary fluctuating wind speeds for long-span bridges[J]. Journal of Engineering Mechanics, 2014,140(1):61-73.

[ 6 ] DEODATIS G. Simulation of ergodic multivariate stochastic processes[J]. Journal of Engineering Mechanics, 1996,122(8):778-787.

[ 7 ] 李永乐，周述华，强士中. 大跨度斜拉桥三维脉动风场模拟[J]. 土木工程学报，2003，36(10):60-65. LI Yongle, ZHOU Shuhua, QIANG Shizhong. Simulation of three-dimensional fluctuating wind field for large span cable-stayed bridge [J].China Civil Engineering Journal,2003,36(10):60-65.

[ 8 ] 赵林.风场模式数值模拟与大跨桥梁抖振概率评价[D]. 上海：同济大学，2003.

[ 9 ] 丁泉顺.大跨度桥梁耦合为颤抖振响应的精细化分析[D]. 上海：同济大学,2001.

Interpolation algorithm for fluctuating wind field simulation of long-span bridges

ZHU Zhiwen1，2， HUANG Yan1

(1. College of Civil Engineering, Hunan University, Changsha 410082, China；2. Hunan Provincial Key Laboratory for Wind and Bridge Engineering, Hunan University, Changsha 410082, China)

In order to reduce computing time of power spectrum matrix decomposition in the harmonic superposition method, the biaxial-and recursive interpolation algorithms were respectively proposed based on the characteristics of double-indexing frequency and spectral matrix. The efficiency of the two methods was proved to be significantly improved, and from the statistical point of view, the simulated results were consistent with those using the traditional harmonic superposition method. The effects of the interpolation node number and arrangement form on the simulated results of wind field were also investigated. It was shown that a fine arrangement at front end and a coarse arrangement at rear end in frequency domain, and a uniform node distribution for double-indexing frequency are more suitable to wind field simulation. Through simulating a fluctuating wind field of a long-span suspension bridge, the effectiveness of the proposed methods was verified.

wind field simulation；long-span bridges；harmonic superposition method；interpolation algorithm

国家重点基础研究发展计划(2015CB057701；2015CB057702)；国家自然科学基金(51278191)；湖南省高校创新平台开放基金(13K016)

2015-10-23 修改稿收到日期： 2016-03-02

祝志文 男，博士，教授，博士生导师，1968年生

黄炎 男，博士生，1987年生

U441

A

10.13465/j.cnki.jvs.2017.07.024