谱分析-相关函数法

一、 基本概念

已经提出的海浪频谱很多，其中大部分是由观测到的波要素连同某些假定推导出来的，大部分则利用定点波面记录通过特殊的谱分析方法得到。后一方法是目前得到海浪谱的主要手段。

在固定点连续记录到波面t，通常认为它是弱平稳的过程，其相关函数为：

REtt （1.1） 由已有理论可知此过程的单侧谱为 SRe02it dt （1.2）

假定海浪为具有各态历经性的平稳随机过程，可利用过程中的现实（一次波面记录）的离散值x1,x2,...,xn计算相关函数

N1ˆtRxnxn,0,1,2,...,mN （1.3） n1ˆtRˆtR式中，N为样本容量；N为乘积xnxn的个数。由此相关函数并参照式（1.2）可得谱的估计值为

2mˆˆ SRteitt,t

0（1.4）

另一方面，我们定义谱密度函数

ˆlim1X2 （1.5） S

T2T对于离散值，

Xxneintt

n1N（1.6）

代入式（1.5），TNt，可得

ˆ1 S2Nttintxetn2Nn1N2xenn1N2int

（1.7）

当t1时，上式变为

ˆ1 S12Nxenn1N2in ， （1.8）

而

ˆtSˆt S1（1.9）

式（1.8）右侧称为周期图，它可通过对样本实行离散傅里叶变化得到。

因此估计谱通常有两种途径，其一通过相关函数，其二通过周期图。在每一途径中又可采用不同的方法。不管用何法，都要对实测记录取离散值，并进行中心化处理。采样间隔的选取，非常重要。在图（1.1）中，细线代表谱中圆频率为1的组成波，今按时间间隔t读取波面值，连接这些离散值得粗线所示的圆频率为21的波动。容易推想，许多高频率的波动可表现为同一低频的波动。 设定义圆频率

Nt （1.10） 则可证明频率2N,4N,...的波动，由于离散化的结果均变现为频率N的波动。

设r,k都是整数，tkt，则

aei2raeit （1.11）

N

图1.1 折叠频率说明图

这意味着利用波面记录离散值进行谱估计时，将使谱于频率间隔

...,5N,3N,3N,N,N,3N,3N,5N,...内的能量都全部叠加到

间隔N,N内，这不仅使谱值的分布范围缩窄至频率范围N,N内，而且得到的谱值不是真实的。

N称为Nyquist圆频率或折叠圆频率，而

fNN212t （1.12） 叫做Nyquist频率或折叠频率（Folding frequency）。

因此，估计谱的分布范围取决于t的值。事实上，海浪谱通常集中在较窄的频域内，通常可从记录中选择最短波的频率fc作为谱估计的频率上界，把大于fc的高频部分切去不计，故fc叫做切割频率（Cut-off frequency）。选取t时应使

t12fc （1.13）

海浪可视为随机过程，但可供使用的定点波形记录具有下列局限性：<1>记录次数是有限的；<2>记录长度是有限的；<3>计算时使用以一定间隔读取的波面数值。理论上可证明，即使计算本身无误差，由此记录得到的结果并非谱的真值，而是对真值的某种估计，故为了

说明估计值接近真值的程度，尚需利用一些统计上的特征量（偏度、方差、置信度等）加以描述。

二、 由相关函数估计频谱 2.1 计算相关函数

设采样时距为t，则式（1.3）的相关函数可改写成

1NRtxtntxtn （2.1） Nn1t,0,1,2,...,m这样便得到R的m+1个值，它们等间隔地分布着，并分别位于

0,t,2t,...,mt。

2.2 估算谱粗值

将式（2.1）代入前文公式以数值积分计算谱值。由于折叠的影响，谱值系在f0~fN范围内进行计算。等间隔地取m+1个频率

f0,f1,f2,...fmfN。我们的目的就是计算谱于这些频率所具有的值，

令Ln代表频率fn对应的粗谱值，得 Ln2RcosdRtcos2ftt （2.2） nn00mt2m如数值积分中采用梯形公式，谱值为

m1211LnR0Rtcos2fntRmtcos2fnmtt 221 （2.3）

n0,1,2,...,m此处采用的频率间隔为

12mt （2.4）

n1fNnfm2tffNm代入式（2.3），得

m12t1n1LnR0RtcosRmtcosn （2.5） m212n0,1,2,...,m2.3 谱的平滑

以上估计出的谱值Ln是不精确的，由它们给出的谱曲线参差不齐。因为样本容量N是有限的，故式（2.1）计算相关函数时，对于小的

，乘积的个数越多，从而Rt的值较可靠；而对于大的，Rt的

可靠性较差。为了改进精确度，可令不同的Rt具有不同的权。这种权函数的形式很多，其中一种常用的为

0.540.46cos,m D （2.6）m0,,mtmm此权函数乘以式（2.2）中的Rt，最后可得谱值

Sn0.23Ln10.54Ln0.23Ln1 （2.7） 对于两个端点频率，可取

S00.54L00.46L1Sn0.46Lm10.54Lm （2.8）

另一种常用的权函数为

11cos,m2 D （2.9） m0,,mtmm同上，由此权函数可得谱值：

Sn0.25Ln10.50Ln0.25Ln1 （2.10） 对于两个端点频率，系数均为0.5.

由上可见，此处所谓改进谱的质量，实际上是采用特定的系数，对谱的粗值进行平滑，而权函数D称为延时窗，前者叫做哈明

（Hamming）窗，后者叫做哈宁（Hanning）窗。此权函数的傅里叶变换

QfDcos2fd （2.11） 称为谱窗。

图2.1为由数值模拟方法得到的相关函数和谱的示例，估计谱时采用了哈明窗。



图2.1 由数值模拟得到的相关函数和谱

2.4 确定置信度

设X1,X2,...都是符合标准正态分布的随机变量，且mx0，则

i2YX12X2...Xk2的分布成2分布（具有k个自由度）：

f21k2k222k2122e （2.12）

ˆf，谱的真值为Sf，已证明随机量又设平滑后的谱值为SˆfSf遵从自由度为k的2分布。根据估计谱值的概率分布，可kS利用置信界限来表示估计值的可靠性。设给定置信水平为（以%表示），则可由2分布确定上界和下界，使估计值落入此界限内的概率为。

ˆfSf，为了确定置信界限，我们依2分布求设以2k代表kS出两个正数a和b，使

12 （2.13） 1P2kb2P2ka从而

Pa2kb

式（2.13）中的a和b可利用已编制的2分布概率表查得。上式可写为

ˆfSPakb Sf或

kˆˆffSfkS P S （2.14）

ba上式表明，对于给定的置信水平，置信上上下界分别为

ˆfkSˆfS1b （2.15）

ˆfkSˆfS2a以上为频率f对应的置信界限，各频率的置信界限构成置信带。

式（2.14）中的自由度k，因使用的延时窗而异。在海浪谱估计中，常使用Tukey导出的结果，即 k2采用与上式最接近的整数。 2.5 参量的选取

谱估计涉及到一系列参量的选取，如样本长度trNt，最大推移乘积个数m等，它们的选取直接影响到谱的质量。

<1>理论上，样本越长，统计特征值越稳定，但其计算工作量较大。另外，海浪并不是严格平稳的，记录时间太长可能会使平稳性受到影响。一般对小的波浪，样本可短些；周期大的波浪，样本宜长些，通常可取10~20min，波数不宜于100个。Arhan在北海北部连续测波22h45min，把记录分成80段，每段长17min，分析得各段的有效波高随时间的变化曲线，发现在波浪迅速成长时，如时段大于17min，波浪不平稳。Haver在挪威沿海收集了384组波浪资料，每段长约17min，发现有15~30%的时段不平稳。故认为时段长于17min是不可取的。

<2>如上文所述，t的选取必须充分小，以避免折叠影响，同时

t过大时失去信息过多，会使估计得谱变形。但如t过小，增加样本

N1 （2.16） m4容量，且所得序列数据相关性增大。通常t值应满足式（1.13）即 t12fc （2.12）

fc为切割频率。

合田良实建议取t110~120TH，即在一个有效波周期内采用

1310~20个样。也有人建议取t0.12T。如由此估计所得谱在fN附近的谱值明显地大于零，应缩小t。反之，如fN远小于处的谱值已接近于零，可适当地加大t。

<3>最大推移乘积个数m对谱估计结果有相当影响。估计谱的质

ˆf相对于真值Sf的均方误差来度量，容易证明： 量要用估计值S2ˆˆf （2.17） ESfSf b2fDS此处

ˆfSf （2.18）偏差 bfES

ˆfESˆfESˆf （2.19）方差 DS 2ˆf表示估计值的离散程度，两bf表示估计值偏离真值得程度；DSbf正比于1m2，者的值越小，表明谱估计结果的质量越高。一般来说，ˆf正比于mN。故m增大时，bf即偏度减小，但DSˆf即方差DS增大。另一方面，m增大时，由式（2.16）知自由度k减小，使估计谱的置信带加宽。同时由式（2.4）知f减小，即谱的分辨率提高。能否通过式（2.17）对m值进行优选，国外学者曾对此做过研究，但尚无可供使用的结果，只能根据经验或通过比较选取。图2.2为谱估计结果随m的变化情况。图中Hs是由谱矩求得的平均波高。m太小时，因点据稀少（谱估计的频率间距mt），不能确切地显示谱形，

峰值偏小。随着m增大谱形显露。m再增大，谱形振动，谱峰尖突，故可由此选定适合的m值。另一种考虑是由各态历经性，当增大到

R趋近于零，一定值是，可选相关函数近于零时相应的m值。第三，

不同m值对应的谱的特征值也不同，可选择谱特征值（谱宽参数）趋于稳定时所对应的m值。

图2.2 谱估计结果随m值得变化

上述通过估计结果来选取m值比较麻烦。国家海洋局《海洋调查规范》建议选取mN25~30t0.5~1.0s。笔者由数模结果建议选取mN15~20t0.5s；mN20~40t1.0s。 2.6 由互相关函数估计互谱

研究任意两个随机过程之间的相互关系时，需计算其互相关函数和互谱。如前文所述，经中心化后任意两个随机过程Xt和Yt的互相

关函数为

1RxytNNn1xtytnnt1Nxtntytn （2.20） Ryxt Nn1t,0,1,2,...,m令

A1xyt2RxytRyxt B1xyt2RxytRyxt 则可得各粗谱值： 粗同相谱

Cˆ2tm1nnxynAxy02Axytcos1Axymt 1m粗转相谱

Qˆ4tm1xynBxytsinnm 1粗互谱

SˆxynCˆxynIQˆxyn 式中，nnnmt。

采用Hamming窗或Hanning窗进行平滑，得

SxynSˆxyn112SˆxynSxyn1

Sxy012Sˆxy02Sxy1 Sxym2Sˆxym112Sˆxymn1,2,3,...,n1式中，Hamming和Hanning窗的系数分别为0.23和0.25。2.21） 2.22） 2.23） 2.24） 2.25） 2.26）

（ （ （ （ （ （

三、程序实现

3.1 程序一：随机信号的自相关函数、协方差函数的实现

%% 随机信号的建立

a=randn(2000,1); % 随机信号 wc=[0.45,0.65];N=79;

%% 基于窗函数的fir滤波器

window=blackman(N+1); % 调用格式：w=blackman(n),根据长度n产生一个布拉克曼窗w；

h=fir1(N,wc,window); % 调用格式：fir1(n,Wn,'ftype',Window),n 为阶数,Wn是截止频率；

x=filter(h,1,a); % [Y,Zf] = filter(B,A,X,Zi)，输入X为滤波前序列，Y为滤波结果序列，B/A 提供滤波器系数，B为分子， A为分母 %% 时域信号

subplot(2,2,1),plot(x),title('\\bf时域信号') grid on

%% 自相关函数

[c,n]=xcorr(x,10,'coeff');

subplot(2,2,2),stem(n,c,'filled'),title('\\bf自相关函数') grid on

%% 协方差函数

[b,m]=xcov(x,10,'coeff');

subplot(2,2,3),stem(m,b,'filled'),title('\\bf协方差函数'); grid on

%% 概率密度函数

subplot(2,2,4),pwelch(x,33,32,[],500),title('\\bf概率密度函数'); % [Pxx,f]=pwelch(x,window,noverlap,nfft,fs)； grid on %% end

图3.1 随机信号的自相关函数的实现

3.2 程序二：自相关函数法、傅立叶变换法及最大熵估计法

%% 傅立叶变换、自相关函数法以及最大熵估计法对离散数据进行谱分析 %% 离散数据的建立

t=0:0.001:1;e=randn(size(t)); % 噪声

data=4*sin(25*2*pi*t)+4*sin(50*2*pi*t)+e; %% 离散数据的截取

subplot(2,2,1);plot(1000*t(1:50),data(1:50));xlabel('\\ittime(mm)'); title('\\bf一元时间序列直观图'); %% 自相关函数法

Fs=1000;NFFT=1024;Cx=xcorr(data,'unbiased');Cxk=fft(Cx,NFFT);Pxx=abs(Cxk); t=0:round(NFFT/2-1);k=t*Fs/NFFT;P=10*log10(Pxx(t+1));subplot(2,2,2); plot(k,P);

title('\\bf自相关函数法'); xlabel('\\it频率（Hz）'); %% 傅立叶变换法 Y=fft(data,512)

Pyy2=Y.*conj(Y)/512; f2=1000*(0:256)/512; subplot(2,2,3);

plot(f2,Pyy2(1:257)); title('\\bf傅立叶频谱图'); xlabel('\\it频率（Hz）'); %% 最大熵估计法 subplot(2,2,4); Fs=500; NFFT=1024;

pyulear(data,20,NFFT,Fs); %% end

图3.2 离散数据的谱分析方法对比图

3.3 程序三：随机波的谱分析

读取试验数据：data=load('wave.txt');

图3.3 随机波幅值变化图

图3.4 随机波谱分析图

3.4 程序四：经典谱分析算法

基于中海油、能威设计有限公司和中国海洋大学工程学院合作的双船浮托整体弃置模型试验的船体六自由度数据，利用自相关系数法、傅立叶变换法及最大熵估计法进行谱分析。

%% 105工况分析

data1 = xlsread('工况105.csv','B6:J15005');data2 = xlsread('静水135.csv','B6:J6'); orig_a0 = zeros(4,3);orig_a0(1,:) = data2(1:3);orig_a0(2,:) = data2(4:6); orig_a0(3,:) = data2(7:9);disp_aa00 = data1;

a_loc = [360 -1527.5 427;360 -1642.5 427;360 -1642.5 327;360 -1527.5 327]; [c,~] = find(isnan(disp_aa00)==1);disp_aa00(c,:) = [];mainfunction % sway 横荡 subplot(3,2,1);

plot(disp_w_r(:,1));xlabel('\\itTime');ylabel('\\itAmplitude');title('\\bfsway(横荡)'); grid on

% surge 纵荡 subplot(3,2,3);

plot(disp_w_r(:,2));xlabel('\\itTime');ylabel('\\itAmplitude');title('\\bfsurge(纵荡)'); grid on

% heave 垂荡 subplot(3,2,5);

plot(disp_w_r(:,3));xlabel('\\itTime');ylabel('\\itAmplitude');title('\\bfheave(垂荡)'); grid on

% roll 横摇 subplot(3,2,2);

plot(rot_w_r(:,2));xlabel('\\itTime');ylabel('\\itAmplitude');title('\\bfroll(横摇)'); grid on

% pitch 纵摇 subplot(3,2,4);

plot(rot_w_r(:,1));xlabel('\\itTime');ylabel('\\itAmplitude');title('\\bfpitch(纵摇)'); grid on

% vaw 艏摇 subplot(3,2,6);

plot(rot_w_r(:,3));xlabel('\\itTime');ylabel('\\itAmplitude');title('\\bfvaw(艏摇)'); grid on

图3.5 船体六自由度图

%% 相关函数法、傅立叶变换法及最大熵法进行谱分析 % sway 横荡 subplot(2,2,1);

plot(disp_w_r(:,1));

xlabel('\\itTime');ylabel('\\itAmplitude'); title('\\bfsway(横荡)'); grid on

% 自相关函数法

Fs=1000;NFFT=1024;Cx=xcorr(disp_w_r(:,1),'unbiased'); Cxk=fft(Cx,NFFT);

Pxx=abs(Cxk);t=0:round(NFFT/2-1); k=t*Fs/NFFT;P=10*log10(Pxx(t+1)); subplot(2,2,2); plot(k,P);

title('\\bf自相关函数法');xlabel('\\it频率（Hz）'); grid on

% 傅立叶变换法

Y=fft(disp_w_r(:,1),512)

Pyy2=Y.*conj(Y)/512;f2=1000*(0:256)/512; subplot(2,2,3);

plot(f2,Pyy2(1:257));

title('\\bf傅立叶频谱图');xlabel('\\it频率（Hz）'); grid on

% 最大熵估计法 subplot(2,2,4); Fs=500; NFFT=1024;

pyulear(disp_w_r(:,1),20,NFFT,Fs); grid on

图3.6 船体横荡谱分析图

图3.7 船体纵荡谱分析图

图3.8 船体垂荡谱分析图

图3.9 船体横摇谱分析图

图3.10 船体纵摇谱分析图

图3.11 船体艏摇谱分析图

四、结论分析

通过对离散数据以及试验随机波的谱分析，我们可以得出，利用自相关函数可以进行波浪谱的估计，得到比较好的分析结果。现实中，无论是离散的正弦线性波的叠加还是现实中的随机波，都会存在一定的高频率波，我们一般只关心其低阶波，致使自相关函数法估计谱频率出现部分干扰，而傅立叶变换法能很好的弥补自相关函数法的不足。综合考虑相关系数法及傅立叶变换法进行谱估计，具有更高的可信度。