现场监测脉动风速的APES法幅值谱相位谱估计*

王向杰， 淡丹辉， 闫兴非， 程 纬

(1.同济大学土木工程学院 上海, 200092； (2. 上海市建设设计研究总院 上海, 200125)

引 言

对于超高、大跨结构而言，风荷载引起的结构振动可对结构的安全造成严重威胁，风荷载成为此类结构的设计控制因素[1]。对风荷载的监测、分析、建模，以及对结构风致效应的分析，不仅成为结构前期规划设计的重点难题之一，也成为建成结构服役期安全性能评估关心的一项重要内容[2]。脉动风速谱是描述其能量频域分布特征的一种重要参数。

脉动风速谱包括功率谱和幅值谱两种。已有规范中给出的均为前者，但只能反应脉动风速幅值的二阶数据特征，不能表征其他概率特性，也丢失了包括相位在内的更丰富的信息，难以完全还原脉动风的随机性原貌[3]。因此使用已知脉动风功率谱产生模拟风速时程，常常需要用到谐波合成法[4-5]或线性回归滤波器法[6]，这些方法计算量大，而且质量不佳。采用幅值谱和相位谱的方式来表征脉动风速，不仅可以保留脉动风速的全面统计特性，从而得到更加接近真实情况的脉动风速模拟时程，而且可以大大降低模拟计算的工作量。李杰等[7]利用实测数据拟合得到了基于Fourier变换的随机脉动风幅值谱，证明了幅值谱在反映脉动风统计特性方面的优势，但由于Fourier变换的频率分辨率较差，对幅值和相位的估计精度也较差，因此得到的幅值谱质量不佳。

近年来，得益于结构监测技术的发展，很多工程结构上安装有风速风向监测系统，笔者采用上海闵浦二桥监测系统的风速监测数据，利用APES方法在频率分辨率和幅值相位估计方面的优势[8-11]，对其进行风速的幅值和相位谱估计，并对估计结果进行统计分析。

1 基于APES算法的脉动风速特征表示

1.1 脉动风速的信号模型

在平稳风条件下，可以将实测得到的紊流风速看成平稳的高斯过程，由多个频率成分的谐振分量与平稳噪声成分叠加而成，即

(1)

其中：u(t)为水平方向的脉动风速；ωr为第r个谐振成分的频率;U(ωr)和φ(ωr)分别为该成分的幅值和相位；w(ωr)为该频率处的噪声。

在任意感兴趣的频率ω处，风速可以被写成一个经幅值函数a(ω)调制的正弦项与一个噪声残余项wω(t)之和

u(t)=U(ω)sin(ωt+φ(ω))+wω(t)

(2)

其中：U(ω)，φ(ω)分别为定义在频率ω处的幅值谱函数和相位谱函数；wω(t)包含了除频率ω外所有的其他频率点的衰减正弦项及一个随机噪声成分ω(n)。

1.2 基于APES的脉动风幅值谱相位谱估计

APES法是一种高精度的非参数幅值和相位谱估计方法，这类非参数谱估计方法起初多用于对复信号、非阻尼信号的处理[8-9]。Stoica等[12]给出了APES和Capon方法的实信号实现形式。由于脉动风速信号为实信号，故笔者采用了后者改进的APES谱估计方法来对脉动风速信号进行谱估计，其过程如下。

将式(2)进一步写为

u(t)=UI(ω)cos(ωt)+UQ(ω)sin(ωt)+wω(t)

(3)

其中：UI(ω)=U(ω)cosφ(ω)；UQ(ω)=U(ω)sinφ(ω)。 设通过连续的风速监测，获得脉动风的时间序列u(t)的N个观测数据u(n)，n=0,1,…,N-1。则每一个观测数据都可以写成如下形式

u(n)=UI(ω)cos(ωn)+UQ(ω)sin(ωn)+wω(n)

(4)

构造观测数据的M×L阶Hankel阵

Y=[yl]

(5)

设

则有如下信号模型

yl=AM(ω)Bω(l)θ(ω)+nl

(6(a))

(6(b))

(7)

其中：nl={wω(l),wω(l+1),…,wω(l+M-1)}T。

(8(a))

约束条件为

hT(ω)AM(ω)=cT

(8(b))

式(8(a))等价于使滤波后的信号的信噪比(SNR)最大化，即

(9)

其中：噪声协方差定义为Q(ω)=E(nl·nl)。

(10)

式(9)的优化问题进一步可转化为

(11)

Stoica等[10]给出的最优滤波器估计器如下

(12)

由式(8(b))，有

即

注意到G的对称性，即

(i= 1,2,…,M)

(13)

根据式(3)，可以给出最终的脉动风速的幅值谱和相位谱估计为

(14)

2 APES算法与Fourier变换在信号处理中的比较

已有研究表明，APES谱估计方法是一种特殊形式的离散Fourier变换(discrete fourier transform，简称DFT)，前者对频率和幅值的识别精度要优于后者。为了比较二者在谱识别方面的性能表

现，给定一段测试信号进行谱识别。该信号包5个谐波分量，其频率分别为：0.053，0.056，0.077，0.081，0.093 Hz，幅值真实值分别为12.3，10.0，5.3，10.0，10.0 m。信号长度为600 s，取样频率为1 Hz。在信号时程中加入一个均值为0，标准差为0.5的随机噪声。分别用APES法和DFT进行处理，谱图在图1中给出，识别得到的结果见表1。

图1 APES算法与DFT识别结果对比Fig.1 Amplitude spectrum identified by APES and DFT

表1 测试信号的幅值谱识别结果

由表1和图1可见，通过计算结果可知，APES法以0误差的方式识别得到了全部5个频率成分，而DFT的识别结果相对精度较低，其误差最大可达到1.2%。对于幅值识别结果，APES法同样具有明显识别精度优势，其最大相对误差仅为1.1%，而DFT的幅值识别精度很差，且识别结果均小于真实值，其相对误差最大达到23.2%，最小为4.9%。从图1可见，对于前两个频率成分非常接近的信号，DFT识别得到的两个峰值在顶端几乎重合，而APES法得到的两个峰值在底部仍然能够有较大的距离，这反映了后者比前者具有更好的频率分辨率。

由此可见，采用APES方法，不仅可以得到高质量的风谱幅值估计，而且因为更高的频率分辨率，可以得到比Fourier法幅值谱更密的频率线，使风谱在全频率域内的谱线分布更趋合理。另外，由于风谱估计的常用时长为600 s，其包含各频率信号的完整周期数最小只有31.8个周期，在这样的短的时长内，APES法的幅值识别结果比Fourier法更可靠。

3 闵浦二桥风速监测数据处理

上海闵浦二桥是一座公轨两用一体化双层特大桥，位于上海市黄浦江上游闵行-奉贤段。主桥为独塔双索面双层斜拉桥，主跨为251.4 m，锚跨为147 m+38.25 m。上层为二级公路双向4车道，桥面宽度为18 m；下层为双线轻轨(上海轨道交通5号线闵奉段)，最小功能宽度为10 m。闵浦二桥结构健康监测系统监测，综合考虑桥梁结构特点，选择了多种监测目标，共有测点172个，其中风速风向测点有2个，分别位于塔顶与主跨跨中。塔顶为150 m，跨中为29 m(吴淞零点高程)，均小于当地梯度风高度，采样频率为4 Hz，采用不间断采样方式监测作用在结构上的风速。

图2 闵浦二桥鸟瞰图Fig.2 Aerial view of minpu second bridge

图3 闵浦二桥风速仪布置图Fig.3 The layout of anemorumbometers in Minpu Second Bridge

采用2016年1～3月之间的3个月(2 184 h)的风速监测数据为研究对象，根据中国交通部规范《公路桥梁抗风设计规范》(JTG/T D60-01-2004)[13]规定的10 min为平均风时距，将其分解为13 104个样本，每个样本包括塔顶与跨中两组数据。

(15)

对于取样频率4 Hz，时距10 min的子样本来说，其计算公式为

(16)

在进行正式风谱估计之前，对得到的监测数据样本进行筛选处理，筛选原则为[16]:a.两处风速原始数据中风速超过64 m/s的数据点小于50个;b.跨中与塔顶两处的主风速均大于6 m/s；c.跨中主风速小于塔顶主风速；d.跨中平面风速小于塔顶平面风速。

通过以上筛选，共得到1 014个符合平稳风要求的子样本。按下式计算各子样本的水平面内脉动风速时程[17]

(17)

图4(a)和(b)给出由计算得到的其中两个时段内的脉动风风速时程。

图4 典型脉动风速时程曲线Fig.4 Typical fluctuating wind velocity time history

4 平稳脉动风幅值谱与相位谱估计结果及统计分析

4.1 平稳脉动风幅值谱和相位谱估计

针对得到的1 041个平稳风速数据样本，分别用APES法与DFT估计各个样本的幅值谱和相位谱。图5给出其中一个样本的幅值谱和相位谱估计曲线。由图可见，在频率的对数坐标轴上，APES法的谱曲线比Fourier谱更加平滑，在低频段，前者比后者的数据点多；在高频段，后者数据则过于密集，反而掩盖了总体趋势信息。总体来看，APES幅值谱曲线在全频段内的质量要好于Fourier幅值谱。这主要是因为前者可以实现任意指定频率线上的谱估计，因此可以在低频区加密频率线，使频率点在对数轴上合理分布；而DFT则需要将时域信息平均地投射到事先划定好的等间隔频率线上，其频率点只能是一个等差数列，低频区缺乏数据，而在高频区则显得数据过密。

图5 某一时段脉动风速样本的幅值谱相位谱估计Fig.5 Amplitude spectrum and phase spectrum estimation of a certain wind velocity sample

4.2 幅值谱相位谱的估计结果统计

将上述两种方法估计得到的所有样本的幅值、相位谱作进一步统计分析。用上述两种方法估计的塔顶、跨中两处脉动风速的幅值均值谱与标准差谱分别在图6和图7给出。定义两种方法估计得到的幅值均值谱差异为

(18)

其中：UAPES(ω)，UFourier(ω)分别为APES法与Fourier变换得到的幅值谱。

Fig.6 Mean of the measured fluctuating wind amplitude spectrum obtained using the two method

图7 两种算法得到的脉动风幅值标准差谱

Fig.7 Standard deviation of measured fluctuating wind amplitude spectrum obtained using the two methods

差异如图8所示。相位均值谱在图9给出。

通过对图6～9的分析可知：两种方法得到的幅值谱均值规律一致：低频区域幅值大，高频区域幅值小，在双对数坐标系下，谱线呈近似线性下降趋势，这与已有文献的结论相同，即脉动风的能量集中在低频区域；跨中的幅值谱均值在低频区和高频区均明显小于塔顶的幅值，而在中间频率区间两者非常靠近；两种方法得到的风速幅值谱标准差在低频区大，高频区小，而且塔顶大于跨中，这从另一个侧面反映紊流能量集中于低频区。

DFT得到的幅值谱均值在全频率区段均明显小于APES法风幅值谱均值估计结果，最大差异达到23.37%,最小也有0.63%；这与前文给出的二种方法性能对比结论一致，有理由相信，APES法得到的脉动风幅值均值谱比Fourier幅值均值谱更准确。

两种方法估计得到的相位谱塔顶和跨中的均值均于0。

图8 两种方法得到的脉动风幅值均值谱差异

图9 两种方法得到的脉动风速幅值均值谱

Fig.9 Mean of fluctuating wind phase spectrum obtained using the two method

4.3 幅值谱相位谱的分布估计

为了进一步了解脉动风幅值谱与相位谱在各频率点处的分布类型，图10、图11分别给出塔顶频率为0.005，0.1，1 Hz处实测风幅值与相位数据的分布频数图，3个频率分别代表低频区、中频区和高频区。由于跨中与塔顶的分布类型相同，APES算法较DFT结果更为准确，故此处只给出塔顶脉动风谱幅值与相位的APES法计算结果。图10中横轴为无量纲化的幅值，定义为各频率处幅值与该频率处均值的比值。

由图10可知，在各频率处幅值分布近似一致，均为先急剧上升再急剧下降的趋势，与对数正态分布较为吻合。

由图11可知，各相位区间的频数大致相等，表明脉动风速的相位基本符合零均值且均匀分布在[-π, π]区间的规律。

根据图10,11可知，将塔顶与跨中的脉动风幅值谱均值以及在特定频率点处的幅值分布情况表示在图12中。

Fig.10 Histogram of amplitude at the top of the tower using APES algorithm

图11 塔顶处相位频数分布图

Fig.11 Histogram of phase at the top of the tower using APES algorithm

图12 塔顶与跨中的幅值谱均值与其分布规律

Fig.12 The mean of fluctuating wind amplitude spectrum and its distribution at the top of the tower and middle of the span

5 结束语

平稳脉动风的幅值谱和相位谱可以很好地描述风场概率特性。为了提高谱估计的质量，笔者采用APES法对闵浦二桥现场实测的风速数据进行了幅值谱与相位谱估计，并与DFT估计结果进行对比。结果表明，APES具有更高的频率分辨率和幅值识别精度，是进行现场风谱估计的有效工具。通过对估计得到的幅值谱与相位谱进行统计分析，发现幅值谱均值在双对数坐标轴中近似呈线性关系，各频率处幅值均呈对数正态分布；相位在[-π, π]上符合均值为0的均匀分布，其规律与已有研究一致。本研究结果可以为桥梁所在地区风荷载估计提供方法参考，也为此地区的结构抗风设计提供依据。