王言英，徐继文，赫亮

(大连理工大学船舶工程学院，辽宁大连116024)

浮式海洋结构物的许多力学性能都对海洋环境运动因素产生响应.在诸多海洋环境因素中起主导作用的是海浪，所指的是风生浪，为重力波.波浪运动属于随机过程，通常被定义为平稳的各态历经的随机过程.工程计算与设计不仅关注那些周期性运动性能的幅值，还关注其运动能量在频域的分布.海浪能量谱密度函数是浮式海洋结构物的许多力学性能预报的重要输入信息(输入谱)，计入结构物固有的频率响应函数，可以得到相关性能的能量谱密度函数(输出谱)，该谱的各阶矩对应各种力学性能［1］.

数学上根据Wiener-Khintchine定理，谱函数可以由随机过程时历的自相关函数的Fourier变换得到，这就是目前广泛流行的快速傅里叶变换(FFT)算法的原理.在船舶与海洋工程领域，ITTC筛选了一批广泛采用的海浪能量谱密度函数的半理论半经验的所谓谱展式［2］.随机过程能量谱密度函数计算方法的进一步研究是工程设计发展的需要.就海浪谱而言，上述的谱展式大多为可用于各海域的通式，并未考虑特定海域环境的特殊性.实际上鉴于海域的位置不同，周围遮蔽陆地环境的不同以及水深的不同，海浪运动的能量大小及其频率特性的差异是显著的.如所周知，在工程设计中使用的谱展式的能量直接关系到结构物的可靠性和经济性.为此，根据结构物营运海域的海洋环境资料，建立或修正现用的海浪谱展式，以最大限度地符合当地海域的环境因素是必要的.另外，在必须自行确定环境或结构物的某些力学性能时，所拥有的过程子样的容量常常是受到限制或有限的.所以，寻求更多的计算方法，有效地利用小子样信息或者尽量挖掘有限容量子样的有益信息，以获得有足够置信度的谱函数是工程上需要的.

本文以欧洲北海North Alywn平台对于一次风暴的观测记录为样本，以FFT算法，对JONSWAP谱展式的谱峰升高因子和峰频控制因子进行了修正［3-4］;分别建立了谱分析的最大熵变换(MET)算法，小波变换(WT)算法［5］和自举变换(BT)算法［6］，并同FFT算法给出的结果进行了比较.

1 快速傅里叶变换(FFT)算法

为验证JONSWAP谱展式在欧洲北海海域的有效性，根据1977年11月北海North Alywn平台提供的一次风暴的记录7，应用FFT算法进行了能量谱密度函数分析计算，并同ITTC推荐的JONSWAP谱展式［2］进行了比较分析.该记录共包括409个子样，每个子样长20 min，采样间隔为0.2 s，有6 000个海面瞬时升高值的记录.

409个子样的能量谱密度函数如图1所示，各个谱的零阶、二阶与四阶矩载于图 2.比对JONSWAP谱展式，按最佳拟合得到的谱展式中的参变量的数值见表1.

根据实测谱和JONSWAP谱展式分别可以计算得到峰频海况的能量谱密度函数，进一步计算得到其各阶矩以及相应的统计特征值，诸如有义波高HS，谱峰圆频率 ωp，平均周期等.表2给出的两组分别由FFT和JONSWAP得到的数据表明两者符合较好，应用JONSWAP谱展式可以获得具有工程精度的统计特征预报值.鉴于JONSWAP谱展式本来就源于欧洲北海开发联合计划，同一海域的后报海浪谱分析结果，同谱展式给出的相应结果的符合是令人信服的.

图1 实测409个子样的能量谱密度函数Fig.1 The energy spectral density functions for real data with 409 samples

图2 实测409个子样的谱矩m0，m2和m4Fig.2 The spectral moments m0，m2，and m4for real data with 409 samples

表1JONSAWP谱展式的参变量Table 1 Parameters for JONSAWP formula

表2 实测谱与谱展式给出子样特征值的比较Table 2 Comparison of eigenvalue between measured andspactral expansion

目前在海洋工程领域，由于海浪观测资料的贫乏，工程设计中大都使用JONSWAP谱展式.为了保证结构物营运的可靠性和经济性，充分利用目标海域的风和浪的观测资料，以FFT算法对JONSWAP谱展式中参变量进行校验确认是必要的.

2 最大熵变换(MEM)算法

熵在信息论中是衡量随机事件不定性程度的量，Burg首先提出了最大熵法计算频谱，其主要思路是把相关函数外推至无穷后再进行频域变换.

有N个事件，若每个事件的出现概率为pi(i=1，2，…，n)，且n个事件相互独立，则单位时间间隔内的平均信息量，熵H为

其中，S(ω)为概率密度函数，为正态分布的连续型随机变量的能量谱密度函数.由信息论知，当随机事件是以等概率可能性出现时其熵值达到最大.在连续型随机变量呈标准正态分布的概率密度时，其熵达到极大值.根据变分原理和自回归分析理论，可得x的谱密度函数的表达式为

其中，m=2，3，…，M);

在估计最大熵谱时，模型阶数的选择是一个关键问题.常用的判阶准则有信息论准则(AIC)，自回归准则(CTA)和最终预测误差准则(FPE)，数值试验结果表明，当信噪比较高时上述3种方法确定的阶数M基本一致.根据FPE准则，对于中心化的序列有

取(PEM)M达到最小时的M值作为最大熵谱的最佳阶数;也有采用容量比定阶数，即M=N/20.

对欧洲北海North Alwyn平台给出的风浪实测资料，分别以FFT和MEM算法进行了谱分析计算.为避免频率折叠的影响，通常从波面记录中选取最短波的频率，大于此频率的波的能量小到可以忽略，取大于此频率的fc为截止频率.此外，根据取样原理当 Δt≤1/2fc时，可用取样序列 x(jΔt)(j=1，2，…，N)，恢复原来的连续函数 X(t)，使 fc≤1/2Δt.图3为以不同采样容量N和阶数M计算得到的MEM估计谱的谱型;图4为以不同采样容量N计算得到的FFT估计谱的谱型.表3给出的是考虑不同采样容量N和阶数M得到的谱特征参数均方根波高Hrms，有义波高 HS，平均过零周期 TZ，平均周期T1，峰周期 TP，峰频 ωP，谱宽系数 ξ，以及谱零阶、二阶与四阶矩 m0，m2，m4.

图3N和M对MEM谱形的影响Fig.3 Effect of N and M on MEM spectral profile

图4 N对FFT谱形的影响Fig.4 Effect of N on FFT spectral profiles

采用某一时段波面记录的样本序列直接计算谱值，这等于使用了矩形窗.因为矩形窗对应的谱窗其旁瓣效应较大，计算结果会出现虚假的谱峰或使谱线上下波动，这就是所谓的谱泄漏.为了减少这种泄漏，在使用FFT法计算时先采用余弦半钟式数据窗.由表3可以看出，用FFT法计算的结果与MEM的计算结果相比基本是一致的，由谱的特征值与谱矩而得到的特征波要素均较吻合.就此算例而言，当N=60 000和M=36时，FFT和MEM算法给出的谱特征参数基本一致，而且当子样容量减少到6 000的1/2～1/4时，选择适当的阶数仍然可以保持足够的谱特征参数精度.

表3 快速傅里叶与最大熵谱估计特征参数计算结果的比较Table 3 Comparison of statistical characteristic values given by FFT and MEM algorithms

3 小波变换(WT)算法

随机信号x(t)的小波变换为

ψab为变换与扩展小波，由小波母函数ψ(t)得到:

式中:fb为带宽参数，fc为小波中心频率.其逆变换为

其中，系数C由式(8)确定.

小波谱分析的基本思想是将一时间序列的方差分解成许多分量，其中每一个分量都是在一特定时间的一个特定尺度.在处理有限长度时间序列时，考虑到傅里叶变换内在的有限性，可以应用小波谱分析方法.小波满足的能量守恒方程为

那么，频率谱密度函数S(f)可以从下式推导得出

当采用波数k=2的Morlet小波时，中心频率fc=1.0.

图5 子样2 000～6 000的JONSWAP谱与常规谱Fig.5 JONSWAP spectrum and traditional spectra with sample length of 2 000 to 6 000

图6 子样6 000～2 000的JONSWAP谱与小波谱Fig.6 JONSWAP spectrum and wavelet spectra with sample length of 6000-2000

子样1 000的JONSWAP谱，快速傅里叶谱与小波谱Fig.7 JONSWAP spectrum，FFT and WT spectrum with 1 000 sample points

图8 快速傅里叶与小波算法给出的有义波高的相对误差Fig.8 Relative errors of with FFT and WT Algorithm

对欧洲北海North Alwyn平台给出的风浪实测资料，分别以 FFT和 WT算法进行了谱分析计算［5］.图5为子样容量分别为6 000，4 000和2 000情况下JONSWAP谱与常规谱的比较，图6为同上子样容量的情况下JONSWAP谱与小波谱的比较，图7为子样容量为1 000情况下JONSWAP谱，快速傅里叶谱与小波谱的比较，图8的2种算法在不同子样容量情况下谱给出的有义波高结果的相对偏差.所谓常规的谱分析方法是基于Wiener-Kintchine定理的通过时间序列自相关函数的傅里叶变换得到的谱函数［8］.关于有义波高相对偏差计算分别采用实测子样和根据统计子样特征模拟的时间序列得到的.

同FFT算法相比小波算法更适合于小子样的谱分析，而且对于时间-频率域的信号处理具有良好特性.在船舶与海洋工程领域，基于有限容量或小子样的波浪和结构物力学性能响应的信号处理，海浪的时间-频率谱的计算，以及探讨频率随时间的变化趋势的计算，小波分析算法是一个值得推崇的方法.

4 自举变换(BT)算法

自举(Bootstrap)的意思是毋需借助外界的帮助而依靠自身的主动性来推进和发展.近年来，自举算法已获得广泛用于估算标准差、置信区间、偏差与预报误差.在一些统计应用中，兴趣集中在根据一随机子样的概率分布估算统计量，该概率分布并不确切知道，而且用以估算那个量的统计采样分布也并不确切知道.自举方法允许人们通过在计算机上多次模拟原始试验，实现近似概率计算.原始的自举算法要求依据独立的和相等的数据重构子样.在这种情况下，人们可以从数据中通过随机抽样更换其位置，构建出人为的重复子样.对于时间序列分析，原始方法不能捕捉到附近观测的相关结构，新的自举方法可以做到这一点.对于有模型逼近，相关结构被模拟成少数已知参数和具有独立的误差.鉴于波浪产生的机理不是很清楚与明确的，在无模型自举方法中多采用区间自举逼近［9］.

应用区间自举逼近时，观测子样被分成一些区间，旨在从原始数据序列中捕捉相关因素.本文采用的是移动区间自举逼近方法.应用移动区间自举逼近选择一最佳区间长度L和移动尺度M是必要的.现假定一观测序列 X1，X2，…，XN，系严格取自固定的观测序列{Xn，n∈Z}，典型的移动区间自举逼近算法是，选择L和M，区间数为B=fix((N-L)/M)+1，对于 m=1，2，…，B，有

有些情况下只能以容量较小的子样为依据从事信号的谱分析，Welch方法是一种调制的周期图法，将时间序列划分为重叠或非重叠的区段，用该方法必须在频率分辨率和谱方差之间做出权衡.较长的区段长度或较少的数据点的重合都会提高谱的分辨率，而较少的区段则会降低其方差.但是对于固定的波浪数据长度，区段的数量同区段长度和重合点数总是呈逆反关系的［6］.为了检验分析的鲁棒性，从总的记录中随机地选择了2个时间序列，见图9.

图9 两组波面升高序列Fig.9 Wave elevation values of two sets

图10 N=3 000，6 000与图9对应图Fig.10 The corresponding diagram with Fig.9 then N=3 000，6 000

为方便FFT算法的使用，区间的子样容量总是取为2的幂，这个数总是主宰着谱的分辨率，这是指能够分辨2个相邻谱峰和其他频率的谱峰.以下的计算其区间长度均取为L=512和不同的尺度M.

为了解决短数据长度的谱分辨率同方差的矛盾，必须增加波浪时间数据长度，引进一个区段自举方法.这个方法是对Welch方法的修正，可以称之为Welch方法或者区段自举方法.在各个区段所用的自举方法是相同的，新的自举子样是从B区段分布的数据点随机地移至Q区段的数据组成的，自举的Welch谱是Q区段谱的平均.为了比较的方便，区段的尺度同其移动的尺度是一样的，在该区段应用Welch方法，再建区段数为80.计算结果见图10，其中(a)、(b)为子样容量为3 000的结果，(c)、(d)为子样容量为6 000的结果.从这4张图很难看出由于移动尺度的不同所导致的谱的差别，说明当重新取样的数量足够大时，谱并不随移动尺度而变化.考虑到自举Welch谱对移动尺度不敏感，进行了M=0的数值试验，结果再次表明谱型同移动尺度几乎无关，所以，对于自举Welch方法并不需要区段数据的重叠.另外，数值试验还表明，当数据长度小到1 000时，Welch谱已经严重扭曲，而自举Welch谱仍然保持光滑.

5 结束语

对于有限时间序列信号，为了保证谱估计有足够精度的分辨率和方差，最大熵算法(MEM)和小波变换算法(WT)是有效的方法.对于小容量时间序列信号，为了得到谱峰具有高分辨率和小方差，自举(BT)算法是实现时间序列自身扩展行之有效的方法.对于新购建的时间序列从事谱分析，采用Welch算法或者修正的Welch算法(自举Welch算法)是可行的.

［1］WANG Yanying.Waves and wave loads on offshore structures［M］.Dalian:Dalian Maritime University Press，2003:38-48.

［2］ITTC Report of the specialist committee of waves［C］//Proceedings of the 23 rd International Towing Tank Conference(ITTC).Venice，Italy，2002:497-544.

［3］XU Jiwen，HE Liang，WANG Yanying.Review of JONSWAP spectrum based on storm 149 from North Alwyn［J］.China Ocean Engineering，2003，17(2):283-288.

［4］WANG Yanying，XU Jiwen，HE Liang.A discussion on application of JONSWAP spectral function to wave data analysis for storm 149 from North Alwyn［C］//Proceedings of the 23rd International Towing Tank Conference(ITTC).Venice，Italy，2002.

［5］XU Jiwen，WANG Yanying.The application of wavelet algorithm to spectral analysis of oceanic waves and offshore structure responses［J］.China Ocean Engineering，2009，23(4):635-644.

［6］XU Jiwen，WANG Yanying.Bootstrap and its application to wave data analysis［C］//Proceedings of 3rd Asia-Pacific Workshop on Marine Hydrodynamics(APHydro).Shanghai，China，2006:280-282.

［7］LIU P C.Is the wind wave frequency spectrum outdated［J］.Ocean Engineering，2000，27:577-588.

［8］OPPENHEIM A V，SCHAFER R W，BUCK J R.Discretetime signal processing［M］.Beijing:Tsinghua University Press，2005:277-288.

［9］LIU R，SINGH K.Efficency and robustness in resampling［J］.Annals Statistics，1992，20:370-384.