钢闸门振动信号频谱分析

10.4　钢闸门振动信号频谱分析

水工钢闸门在地震荷载，水流脉动压力、启闭过程的摩擦和水中漂浮物冲击等作用下产生振动，振动响应信号较为复杂(一般为随机振动)，振动频率成分丰富，振动强度激烈，对闸门造成较大危害。为此，必须对闸门振动信号进行频谱分析，以了解闸门的动态特性(固有频率，主振及阻尼等)、闸门的动力响应及振动信号的频率成分等，进而采取一些防震、隔震及消震措施，以保证闸门的安全运行。

在第6章闸门动态检测中曾介绍了闸门动力检测的原理与方法，其中包括对闸门的振动信号的谱分析方法。关于频谱分析问题，本节主要结合闸门的动态检测介绍频谱分析的基本概念，各类信号的频谱分析的原理、数学上的傅立叶变换，数字信号分析仪及系统等内容。

10.4.1　频谱分析基础

频谱分析基础知识包含谱、频谱及频谱分析的基本概念。

10.4.1.1　谱

谱在自然界中普遍存在，如音乐中有乐谱;自然当中有光谱(按其波长大小排列);将复杂的谱音分解为单纯的音调，并按其频率的高低排列，便得到音响谱，等等。通俗说法是，把复杂成分的物质分解成简单的物质成分，并把其具有的特征量按其一定的顺序排列起来，就构成所谓的谱。

10.4.1.2　频谱

频谱这一词汇，在振动分析中用得较多，亦即把振动的时域信号波形分解成许多不同频率的谐波，按各谐波频率的高低排列起来，所形成的波形称为频谱图或简称频谱。

例如，在闸门检测中，所测的振动波形是一个宽带随机波形，若将之分解，使其由许多谐波组成，而这些谐波分量又都具有不同的振幅与相位。以振幅或相位作纵坐标，以频率作横坐标，这就得到了各种振幅和相位以频率为变量的频谱函数，作出的曲线或谱线，即为该振动波形的频谱或频谱图。这样，原测的时域信号便变成了频率域信号了。如图10-8第4项第4栏所示。

10.4.1.3　频谱分析

所谓频谱分析就是把时间域的各种动态信号变换到频率域的分析。分析结果得到以频率为横坐标的各种动态参量的谱线与曲线。如振动幅值谱、相位谱及功率谱等。

(1)幅值谱——动态信号中所包含的各次谐波幅度(振幅)值的全体，它表征振动幅值随其频率的分布情况。

(2)相位谱——表示各次谐波相位值的全体。它表征振动相位随其频率的分布情况。

(3)功率谱——表示各次谐波能量(或功率)的全体，也称功率谱密度。它表征各谐波能量(或功率)随频率分布情况。

随机振动信号主要是用功率谱密度作频谱分析;对确定性信号，一般用幅值、相位谱进行分析。

通过频谱分析可以解决振动信号的频率结构、闸门的动态特性(f、A、C)及闸门的动力响应等情况，同时还可用来对实测波形进行修正，从而得到真实的振动波形。

10.4.2　确定性信号的频谱分析

确定性信号需要做频谱分析的有周期波形信号及非周期波形信号，下面分别予以介绍。

10.4.2.1　周期性信号的傅氏谱

求解任意形式的周期信号分析问题，依高等数学知识，借助于傅立叶(Fourier)分析法，把周期波展开成傅氏级数形式，即可将这个复杂的波形分解成许多不同频率的正弦和余弦谐波之和。

设任意周期函数X(t)的周期为T=(ω为角频率)，则其傅氏级数展开式为:

式中: a_n=X(t)cos nωt d t(n= 1，2，…，∞)

　　　b_n=X(t)sin nωt d t(n= 0，1，2，…，∞)

　　　a₀=X(t)dt

项表示周期T内的振幅平均值;而a_ncos nωt+ b_nsin nωt表示n次谐波，n= 1时正弦波叫做基波。

由a_n cos nωt+ b_n sin nωt=

式中:A_n=+表示n次谐波的振幅;

φ_n=arctan[]表示n次谐波的相位角。

根据A_n、φ_n即可作出波形x(t)的振幅谱及相位谱。周期波形的频谱图是离散的线谱，如图10-11所示。

此外，还可根据欧拉(Euler)公式，将傅立叶级数表示成复数形式:

式中: F(nω)=x(t)e^－jnωt dt　　　(10-63a)

显然F(nω)是个复数，可表示为:

式中: F(nω)=φ_n= arctan()

根据F(nω)、φ_n与频率ω的关系，可画出复数的振幅谱和相位谱，如图10-11所示。

图10-11　周期波形频谱图

如果x(t)不是周期函数，那么在时间t的有限区间［t₀，t₀+ T］内，傅立叶级数及傅立叶系数的公式仍然成立:

上述表明，在有限时间区间上一个复杂波可分解成无限多个简谐波。

10.4.2.2非周期波形的傅立叶谱

非周期波形包括冲击、脉冲和瞬态振动等波形，它们的频谱不能用线谱来表示，因此，非周期波形不能直接展开成傅立叶级数。但可将其看成一个周期为无穷大的周期波形来研究。

周期波的频率是离散的线谱(见图10-11)，而相邻谱线的间隔是波形的重复频率，Δf=。显然，周期T越长，相邻间隔就越小，谱线的密度则越大。当T→∞时，Δf→0。亦即非周期波(含冲击波及瞬态振动等)的频谱是连续频谱，如图10-12所示。

图10-12　非周期波频谱图

研究连续频谱的主要目的是了解其各次谐波的振幅随频率增加而变化的状况，即各谱线长度之间的关系，其次才是谱线长度本身。其振幅随频率增加而变化的状态可用谱线的包络线来表示。如上所述，令T→∞，则上述的无穷级数就变成了无穷积分，这就是傅立叶级数变换。

由图10-12可见，形状相同但周期不同的波形，它们的频谱的包络线具有相似的形状，只是高度不同而已。且周期增加一倍，包络线的高度则按比例降低一半。如果将各谱线长度都乘以该波形的周期T，就得到一条与周期无关的包络线G(jω)，即频谱特性曲线。频谱特性曲线的函数表达式称为频谱函数，它反映振幅随频率的变化规律。

若将周期波的傅立叶级数展开式的复数形式的振幅F(nω)乘以其周期T，则得到频谱特性曲线的积分表达式，即

当T→∞时，上式中x(t)变成非周期波形，此时，Δω→dω，即对nω求和变成对ω积分，并将nω记为jω，于是

由上式可得:

由欧拉公式可推得:

这样，就可以根据非周期函数计算其频谱。傅立叶积分可看做^jωt傅立叶级数的推广，它与x(t)= G(jω)e·dω一起是复数形式的傅立叶变换对。要求x(t)满足绝对可积条件。

10.4.3　随机信号的频谱分析

随机信号的频谱分析通常又称之为谱密度分析。它的频谱是连续谱。用连续来说明随机波的幅值、相位及功率谱等随频率变化的连续分布情况。

频谱分析的数学基础是傅立叶变换，为此，我们首先给读者介绍傅立叶变换。

10.4.3.1　傅立叶变换

傅立叶变换可分为有限离散傅立叶变换(DFT)和快速傅立叶变换(FFT)。下面分别叙述。

1.有限离散傅立叶变换(DFT)

前面已介绍，频谱分析的数学基础是傅立叶变换，即

求解任意时间函数x(t)的傅立叶变换的前提是该时间函数有解析解。但实际测试得到的振动记录曲线是一条连续的波形线，很难用数学表达式描述，而且也不适用数字计算机计算。为此，必须对所记录的连续波形线进行离散采样和模数转换，得到该波形信号的无穷离散数值序列x₀，x₁，…，x_∞。显然，由于计算机的容量有限，不可能接受无限个离散采样数值，只能对有限个(n)的离散数值序列{x_n}进行傅立叶变换，称之为有限离散傅立叶变换(DFT)。

一个连续变化的波形线x(t)，对其离散采样，其频谱分析采样图如图10-13所示。

图10-13　频谱分析采样图

设一次记录的波形样本，长为T，采样时间间隔Δt，采样点数为N(取偶数)，则T= NΔt。设N= 16个采样点，则其各点的数值构成离散值序列

{x_m= x(mΔt)}(m= 0，1，2，…，N－1)

m号采样点处的时刻为t_m= mΔt，该时刻的数据为x_m= x(mΔt)，于是就可以得到在［0，T］区间内x(t)的有限傅立叶级数。

根据式(10-61)，得

在t_m=mΔt(m= 0，1，2，…，N－1)处有:

式中:傅立叶系数a_k、b_k分别为:

进行有限的傅立叶变换计算时，当x(t)为周期函数时，T就是周期，当x(t)为非周期函数时，T就是样本长度，二者计算公式相同。

现用有限离散傅立叶变换方法对图10-13所示采样波形进行频谱分析。

首先对持续时间为T的波形离散采样，采样点数为N，采样时间间隔为Δt，Δt的选取按照采样定理，样本长度T= N·Δt，现将T视为样本的周期，则基频f₁==，各谐波频率f_k= k·f₁(k= 1，2，…)。任一时刻t采样点号m=，由式(10-71)可得

式中:

相应各频率f_k的谐波的振幅及相位为:

以f_k= k·f₁为横坐标，对应f_k的振幅A_k(或相位φ_k)为纵坐标就可得出该连续波形的功率谱和相位谱，这些频谱都是离散线谱。

2.快速傅立叶变换(FFT)

离散傅立叶变换已在计算机上实现了傅立叶变换，但若计算点数很多时，运算量很大，计算就出现麻烦。例如用复数形式表示一离散傅立叶变换为:

时，直接计算x(n)值的工作量太大，对于K个x(n)中的每一个必须作N次x(k)乘以，所以总共有N²次复数乘法运算，而且还要作N(N－1)次复数加法运算，当N较大时，即便使用大容量电子计算机也有困难。如N= 1024，则仅复数相乘的次数就为N²= 1048576，即一百多万次。可见即使离散傅立叶变换在计算机上能算也是不实用的。1965年柯立(J.W.Cooley)和杜开(Tukey)提出了DFT的快速方法，并编出了该方法的计算程序，称之为快速傅立叶变换(FFT)。也就是计算随机信号x(t)的有限傅立叶变换的快速方法。

快速傅立叶变换的具体算法很多，有时域分解法，频域分解法等。在此仅介绍一种时域分解法的基本原理。时域分解法是把离散时间序列{x_n}分割成若干较短序列(如先分成一半，再分成1/4……)，分别计算每个短序列的离散傅立叶变换，然后，想法“合并”，得到原始序列的离散傅立叶变换。具体如下:

首先，利用两个N/2点的离散傅立叶变换表示一个N点的离散傅立叶变换。

现将样本x(k)按顺序标号奇偶性分成两列，并考虑W2nk=Wnk，N N2故

式(10-78)表明:B(n)、C(n)实质上是两个周期为N/2的具有N/2个独立值的周期性频率函数。即

(由于）

所以，x(n)= B(n)+C(n)W^Nn只能给出x(n)的N/2个独立值，而

因为　　　

由此可见，具有N个样本点的时间序列傅立叶变换，可分成N/2个点的傅立叶变换。

同理，按上述分法，再将N/2点的傅立叶变换分成N/4，再将N/ 4点分成N/8点……以至分成m个两点的有限离散傅立叶变换(DFT)为止，求两点的傅立叶变换，只需进行加减换算。这就是快速傅立叶变换的原理。

用BASIC、FORTRAN……甚至C++等各种语言编写的FFT计算程序，已在很多频谱分析仪或作为计算机的动力分析软件中使用，读者可参考有关资料。

10.4.3.2　功率谱密度函数

功率谱函数也称自谱密度函数或简称功率谱。由于功率谱会突出该信号的主频率，故功率谱的主要用途是用来建立信号的频率结构，分析频率组成和相位量的大小。

1.功率谱密度函数计算表达式

由于功率谱密度函数一般可以从自相关函数推导出来。按照傅立叶变换理论可得到(参见前面公式):

特殊情况的τ= 0时，随机信号x(t)的自相关函数为:

由此可得:

式(10-79)表明随机信号的均方值等于G_xx(f) d f对全部频率求和，G_xx(f)d f就等于频率范围内的均方值，即G_xx(f)表示随机波信号x(t)在每单位频带的分量的均方值。由于G_xx(f)是频率f的函数，且是每单位频率上均方值的大小，因此它是均方值的“谱密度”。

在闸门振动中，功和能量与其振幅的平方成正比例，即x²(t)可看做闸门振动体系的能量或功的量度，所以功率谱G_xx(f)是随机振动在单位频带内的谐波分量的能量(功)按频率f分布的量度。

2.功率谱密度函数的计算方法

功率谱密度函数有两种数字计算方法。

(1)按上面所述，通过自相关函数R_xx(τ)作傅立叶变换得到相应的谱密度函数。称之为布克拉门-杜开(Blackman-Tukey)方法。

(2)对原始测试数据x(t)直接进行快速傅立叶变换得到相应的谱密度函数。称之为柯立-杜开(Cooley-Tukey)方法。

此外，还可以通过模拟滤波技术求得。此方法是在数字信号设备出现前常用的方法，现在仍有个别单位采用。但在数字信号设备(数模转换器及专门分析仪)出现及计算机广泛应用后，功率谱密度就主要通过上述两种计算方法求得。两种数字计算方法相比，从计算效率来看，方法(2)较好，当时域信号采样容量N越大时，越显出方法(2)比方法(1)计算效率高。所以，方法(2)的直接快速傅立叶变换在功率谱计算中得到广泛使用。

10.4.3.3功率谱密度函数的换算

现以自谱密度函数为例，用上述两种数字计算方法对功率谱密度函数换算分别给予说明。互谱密度函数也可采用相应方法。

1.由相关函数求得谱密度函数

根据谱密度函数与相关函数之间的关系式(10-46)，有

其中自相关函数可依式(10-42)得到

若相关函数的计算采用的是数字计算的方法，则G_xx(f)也可以由数字计算方法求得，都可采用数据序列直接计算得到。

2.用直接快速傅立叶变换计算谱密度函数(以自谱为例)

考虑一个平稳随机过程{x(t)}，对其第k个长度为T的时间记录进行有限傅立叶变换(DFT)，有

其谱密度定义为:

式中:期望值运算子E表示是对下标k的平均运算。这个式子定义的谱密度函数，等价于用相关函数的傅立叶变换所定义的相应的函数。(f，T)是x_k(f，T)的共轭复数。

去掉极限值和期望值运算，可得到G_xx(f)的原始估计

若考虑平稳随机过程{x(t)}的双边自谱密度，则有

同理，对于两个平稳的各态历经的随机振动过程{x(t)}和{y(t)}的双边互谱密度为

为求得(f)的离散值，一般采样柯立-杜开(Cooley-Tukey)方法，即对原始数据直接进行快速傅立叶变换(FFT)得到。

设对x(t)采样的点数为N，采样间距为h，则样本长度T= Nh，作FFT时离散频率取

那么，G_xx(f)的离散值G_k为

复数形式表示的离散傅立叶变换为

时，则

此时，是x_k的共轭复数。

故

此即柯立-杜开方法计算功率谱密度函数的公式。

x_k=即为x( t)的傅立叶变换x(if)=在离散频率f_k= kf₀=处的离散值。

显然可见:

(1){x_k}具有周期性，即

x_N+k= x_k

其周期为f_N=N·Δf= N·=。

(2)x_k是复数，其实部是偶函数，虚部是奇函数。因而x_k也具有这种奇偶函数的特性。又x_k以f_N=为周期，故N个复数只有N/2个是独立的。可见，当时域采样容量为N= 1024时谱线为512。

由N个x_n值可得到x_k(k= 0，1，2，…，N－1)。由x_k能否计算x_n(n= 0，1，2，…，N－1)呢?回答是可以的。这个计算过程称之为有限离散傅立叶变换，即

式中:x_k=Δf·x(f_k)。

当f_k= 0时，k= 0，则

正好是{x_n}的均值表达式。

10.4.3.4　加窗处理

用FFT方法计算功率谱密度函数时，由于各种条件的限制，我们实际测得并进行分析的样本函数记录仅仅是整个振动信号中的一小部分。这相当于用一高为1，长为T的矩形时间窗函数乘以原时间函数，即

式中:x(t)——原时间函数;

x(t)_r——实际分析用的样本记录;

b₀(t)——时间矩形窗函数。

时间矩形窗函数的截断必使窗外的信息全部失掉，时域中的这种信号损失将导致频域内附加一些频率分量，从而给傅立叶变换带来误差，这种现象称为泄漏误差。

在对钢闸门振动测试进行功率谱分析时，都要进行加窗处理。对原始数据x(t)进行加窗处理，即乘以一个窗函数，等于对原始数据进行不等的加权，结果会使计算出来的功率谱函数的值减小，因此，对计算结果要进行修正。即加窗后计算出来的功率谱密度函数应乘以修正系数K₀，见表10-5。

常用的窗函数除矩形窗以外，还有海宁窗(Hanning)、汉明窗(Hamming)、钟形窗、1/10余弦坡度窗等，如表10-5所示。具体选用什么窗要根据问题的需要及窗的特点而定。

表10-5　五种窗内容比较表

10.4.3.5　用FFT算法计算G_xx(f)的一般步骤

(1)根据采样定理，对信号x(t)进行采样得到数据序列{x_r}，采样点数最好为N= 2^p，p为任意整数;否则，应对原始数据截去一段或增加一些零点使其长度为N= 2^p，得到新的序列{x_r}，r= 0，1，…，N－1。目前，用得较多的是N= 1024，2048或4096。

(2)运用适当的窗函数ω()对{x_r}进行修正，得到

(3)用FFT计算{}的离散傅立叶变换得到序列{x_k}，在此

(4)利用式(10-89)

即柯立-杜开方法计算谱密度函数的估计值G'_k(f_k)。

(5)为修正窗函数的影响，用适当的比例因子K₀来修正这些谱密度函数的估计值，即根据表10-5中所列K₀值修正G'_k，最后便得到修正后的谱密度函数估计值G_k(f_k):

式中:N₀为采样容量，即

式中:B_e为最大有效分析带宽;e为随机误差，即e=。σ为标准差，μ为均值。

采用不同的时间窗，就有不同的K₀值，这就得到单边功率谱密度函数的估计值G_x(f)的N个离散值G_k(f_k)。

10.4.3.6　数字信号分析的仪器及其系统

在闸门动力测试中，对其振动信号的频谱分析的仪器及系统有:模拟式频率分析仪、数字信号分析仪、计算机信号分析系统和数模混合分析系统等。在电子计算机问世前，由于由测试系统所获得的振动信号多是模拟信号，早期的信号分析多采用人工的模拟量分析法，它比较直观、方便，仪器的功能单一，精度及灵活性方面均较差。1980年前，我们用过此法。随着测试设备的发展，电子计算机的出现为信号分析提供了更快、更方便、更精确的分析手段，出现了数字信号分析方法。此法是通过A/D转换器(模数转换器)将连续模拟信号转换成离散数据，再由电子计算机分析处理，最后由打印机输出数据，或经D/A转换器(数模转换器)将数据转换成模拟量，在示波器或X-Y绘图仪上绘出图像，如功率谱图、凝聚函数图等。目前我们对闸门动力测试分析全采用数字分析法。

自从1965年柯立和杜开提出了计算随机信号的有限傅立叶变换的快速方法(FFT)以来，随着大规模集成电路及电子计算机的飞速发展，数字信号分析仪及其系统便占据了主导地位，成了今后的发展方向。

数字信号分析仪及其分析系统按其控制方式大致可分为4类。

1.软件控制式专用信号分析仪

这类分析仪内装专用计算器，将各种信号分析功能事先编制好程序并录制在磁带上或装在计算器的存储器内，根据分析目的要求，调入相应的程序就可工作。这种分析仪的特点是功能较多，理论上可以无限扩展。但分析速度稍慢，操作不很方便。这类仪器的典型产品有TT08，TT08S，TT17S等。

2.硬件控制式专用分析仪

由于大规模集成电路的发展及快速乘法器等的出现，这类仪器中的FFT及控制逻辑均可用硬件实现，其主要特点如下:

(1)计算速度极快，1024点FFT运算时间在毫秒级。

(2)频率分辨率很高，带有频率细仪(ZOOM)功能的分析仪，其分辨率极高，在低频段很容易达到千分之几赫兹。

(3)操作、显示、存储、记录、再处理等非常方便。这类仪器一般把所有运算功能的程序都固化在ROM(输入、输出存储管理元件)中，操作时只要按相应的功能键，即可完成。CRT(阴极射线显示器)显示结果，一目了然。配备的X-Y记录仪或数字绘图仪记录方便。通用接口GP-IB还可以与其他计算机系统连接，便于进行数据的再处理或组成自动测试系统。

(4)分析仪内装有微型计算器(或微处理器)，可对整个分析仪进行控制管理和一般数据处理。

这类仪器的典型产品有HA I、TD4070(国产)、CF-900系列、SM-2100(日本产)、HP5423A、HP3582A、SD-375(美国产)、B＆K2081(丹麦产)等。

3.以微型计算机为主的信号分析系统

以通用微型计算机为主配以A/D和D/A转换板、必要的外围设备以及各种软件即可构成数字信号采集、分析与处理系统。这类系统的主要特点如下:

(1)软件丰富，通用性强。

(2)兼有信号数据采集、存储、分析、再处理等多种功能。

(3)便于同其他仪器或微机系统连接。系统中还配有GP-IB及RS-232C两种通用接口，它可以同带有GP-IB接口的微机、仪器直接相连，组成可编程的自动测试系统。

在钢闸门动力测试中，曾用过两种分析系统:早期(20世纪80年代)用过上海同济大学与浙江黄岩机械厂研制的，包括AF-1数据采样放大器、AD-1型模数转换器、配用APPLE-Ⅱ微机组成的TJAP数据处理系统和FFT分析软件系统等。绘图仪是日产的DXY-880型绘图仪;随着电子仪器及计算机的发展，我们现在常用的是XR-510多用记录仪，其分析系统由INV306型数据采集分析仪、便携式计算机和打印机组成。该系统由计算机直接控制测点扫描及数据采集，可进行自动调零，进行时域分析和频域分析等。该系统全部采用进口国际标准集成电路设计和制造，性能稳定、抗干扰性能好，使用方便可靠，尤其适用于大数据量的采集和分析处理。

4.小型机同FFT硬件相结合的综合性信号分析系统

这种分析系统兼有小型机，其软件丰富、编程灵活以及FFT硬件运算速度快等优点，相当于上述2、3种仪器系统相结合的产物。其主要特点是:分析功能齐全、处理速度快、能作各种谱分析等;既有专用程序，又可自编程序，也可作通用计算机使用。这是一种较高水平的综合性信号分析系统。典型产品有美国生产的HP-5451C傅立叶分析系统。

常用频谱分析仪(部分)及其基本性能参见表10-6。

表10-6　部分常用的频谱分析仪及其基本性能表

续表