谱分析-相关函数法

海浪谱分析—相关函数法

一、 基本概念

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

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

()()()[]τηητ+=t t E R （1.1）

由已有理论可知此过程的单侧谱为

()()dt e R S t

i ωτπω-∞

?

=

2

（1.2）

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

()()()t R t R

m x x N t R N n n n ?-=?=-=?∑-=+ννννννν??,...,2,1,0,1?1

（1.3） 式中，N 为样本容量；ν-N 为乘积n n x x ν+的个数。由此相关函数并参照式（1.2）可得谱的估计值为

()()t t e t R S t i m ?

（1.4）

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

()()221lim ?ωπωX T

S

T ∞→=（1.5） 对于离散值，

()t e x X t n i N

n n ?=?-=∑ωω1（1.6）

代入式（1.5），t N T ?=，可得

()21

2

1

221?∑∑=?=??=??=N

n t

n i n N

n t

n i n

e x N t t e x t

N S

ω

ωππω（1.7）

当1=?t 时，上式变为

()πωπωω

<=∑=，21

1

21?N

n n

i n e x N

S （1.8）

而

()()t S t S ??=ωω1

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

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

t N ?=πω（1.10）

则可证明频率,...4,2N N ωωωω±±的波动，由于离散化的结果均变现为频率()N ωωω<的波动。 设k r ,都是整数，t k t ?=，则

()t i r i ae ae N

ωωωη==+2（1.11）

图1.1折叠频率说明图

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

()()()(),...5,3,3,,,3,3,5...,N N N N N N N N ωωωωωωωω----内的能量都全部叠加到

间隔()N N ωω,-内，这不仅使谱值的分布范围缩窄至频率范围()N N ωω,-内，而且得到的谱值不是真实的。

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

()t f N N ?==212πω（1.12）

叫做Nyquist 频率或折叠频率（Folding frequency ）。

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

)c f t 21≤?（1.13）

海浪可视为随机过程，但可供使用的定点波形记录具有下列局限性：<1>记录次数是有限的；<2>记录长度是有限的；<3>计算时使用以一定间隔读取的波面数值。理论上可证明，即使计算本身无误差，由此记录得到的结果并非谱的真值，而是对真值的某种估计，故为了说明估计值接近真值的程度，尚需利用一些统计上的特征量（偏度、方差、置信度等）加以描述。

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

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

()()()m

t t x t t x N t R N n n n ,...,2,1,0,11

=?=?+-=?∑-=νντνννν

（2.1） 这样便得到()τR 的m +1个值，它们等间隔地分布着，并分别位于

t m t t ???=,...,2,,0τ。

2.2 估算谱粗值

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

N m f f f f f ==,...,,021。我们的目的就是计算谱于这些频率所具有的值，

令n L 代表频率n f 对应的粗谱值，得

()()()t t f t R d R L n

m

n

t

m n ???==

∑∑=?=νπνπττωτπ

ντ

2cos 2

cos 2

（2.2） 如数值积分中采用梯形公式，谱值为

()()()()()()

m n t

t m f t m R t f t R R L n m n n ,...,2,1,02cos 212cos 02121

1=?????????+??+=∑-=πνπνπν（2.3）

此处采用的频率间隔为

t

m n f n f t

m m f f N N ?=

?=?=

=?21

21

（2.4） 代入式（2.3），得

()()()()

m n n t m R m n t R R t L m n ,...,2,1,0cos 21cos 02121

1=??

?

????+?+?=∑-=ππννπν（2.5） 2.3 谱的平滑

以上估计出的谱值n L 是不精确的，由它们给出的谱曲线参差不齐。因为样本容量N 是有限的，故式（2.1）计算相关函数时，对于小的ν，乘积的个数越多，从而()t R ?ν的值较可靠；而对于大的ν，()t R ?ν的可靠性较差。为了改进精确度，可令不同的()t R ?ν具有不同的权。这种权函数的形式很多，其中一种常用的为

()()??

????=>≤+=t m D m m m m

ττττττπττ,,0,cos 46.054.0（2.6）此权函数乘以式（2.2）中的()t R ?ν，最后可得谱值

()1123.054.023.0+-++=n n n n L L L S ω（2.7）

对于两个端点频率，可取

()()m

m n L L S L L S 54.046.046.054.01100+=+=-ωω（2.8）

另一种常用的权函数为

()()

??

????=>≤???

? ??+=t m D m m m m ττττττπττ,,0,cos 121（2.9） 同上，由此权函数可得谱值：

()1125.050.025.0+-++=n n n n L L L S ω（2.10）

对于两个端点频率，系数均为0.5.

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

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

()()()ττπτd f D f Q 2cos ?∞

∞-=（2.11）

称为谱窗。

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

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

2.4 确定置信度

设,...,21X X 都是符合标准正态分布的随机变量，且0=i

x m ，则

22

221...k X X X Y +++=的分布成2χ分布（具有k 个自由度）：

()

()()(

)2

1222

2221χχχ

--???

?????Γ??? ??=e

k f k k （2.12）

又设平滑后的谱值为()f S

?，谱的真值为()f S ，已证明随机量()()f S f S

k ?遵从自由度为k 的2χ分布。根据估计谱值的概率分布，可利用置信界限来表示估计值的可靠性。设给定置信水平为β（以%表示），则可由2χ分布确定上界和下界，使估计值落入此界限内的概率为β。

设以()k 2χ代表()()f S f S

k ?，为了确定置信界限，我们依2χ分布求出两个正数a 和b ，使

()()

()(

)

212

122β

χβ

χ-=

≥-=

≤b k P a k P （2.13） 从而

()()

βχ=≤≤b k a P 2

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

()()β=???

? ??≤≤b f S f S k a P ? 或

()()()β=??

? ??≤≤f S a k f S f S b k P ??（2.14）

上式表明，对于给定的置信水平β，置信上上下界分别为

()()()()f S a

k f S f S b

k f S ????2

1

==（2.15）

以上为频率f 对应的置信界限，各频率的置信界限构成置信带。

式（2.14）中的自由度k ，因使用的延时窗而异。在海浪谱估计中，常使用Tukey 导出的结果，即

??

?

??-=412m N k （2.16）

采用与上式最接近的整数。 2.5 参量的选取

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

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

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

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

容量，且所得序列数据相关性增大。通常t ?值应满足式（1.13）即

)c f t 21≤?（2.12） c f 为切割频率。

合田良实建议取()3

1201~101H T t =?，即在一个有效波周期内采用

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

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

量要用估计值()f S

?相对于真值()f S 的均方误差来度量，容易证明： ()()()()()[]

f S D f b f S f S E ??22+=??

????-（2.17） 此处

偏差()()[]

()f S f S

E f b -=?（2.18） 方差()[]()()[](

)

??

????-=2

???f S E f S E f S D （2.19） ()f b 表示估计值偏离真值得程度；()[]

f S

D ?表示估计值的离散程度，两者的值越小，表明谱估计结果的质量越高。一般来说，()f b 正比于21m ，

()[]

f S

D ?正比于N m 。故m 增大时，()f b 即偏度减小，但()[]

f S D ?即方差增大。另一方面，m 增大时，由式（2.16）知自由度k 减小，使估计谱的置信带加宽。同时由式（2.4）知f ?减小，即谱的分辨率提高。能否通过式（2.17）对m 值进行优选，国外学者曾对此做过研究，但尚无可供使用的结果，只能根据经验或通过比较选取。图2.2为谱估计结果随m 的变化情况。图中()s H 是由谱矩求得的平均波高。m 太小时，因点据稀少（谱估计的频率间距t

m ?=

?π

ω），不能确切地显示谱

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

可选相关函数近于零时相应的m 值。

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

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

上述通过估计结果来选取m值比较麻烦。国家海洋局《海洋调查规范》建议选取()()s

25=

?

=。笔者由数模结果建议选取

~

30

N

m0.1

t

~

5.0

()()s

40

N

~

m0.1

?

t

=。

20=

t

=；()()s

N

20

m5.0

?

15=

~

2.6 由互相关函数估计互谱

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

X和()t Y的互

相关函数为

()()()

()()()m

t t y t t x N t R t t

y t x N t R n N n n yx n

N n n

xy ,...,2,1,0,11

11

=?=?+-=??+-=

?∑∑-=-=νννννννν

νν

ν（2.20） 令

()()()[]

t R t R t A yx xy xy ?+?=?ννν2

1

（2.21） ()()()[]

t R t R t B yx xy xy ?-?=

?ννν2

1

（2.22） 则可得各粗谱值： 粗同相谱

()()()()()??

?????-+??? ???+?=∑-=t m A m n t A A t C xy n m xy xy n xy 1cos 202?1

1πννπων（2.23） 粗转相谱

()()??

? ????=∑-=m n t B t Q m xy

n xy πννπωνsin 4?11（2.24） 粗互谱

()()()n

xy n xy n xy Q I C S ωωω???-=（2.25） 式中，t

m n

n n ?=

?=πωω。

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

()()()()()()()()()()()()()()

1,...,3,2,1?21?22?21?21?110011-=-+=+-=+-+=-+-n n S S S S S S S S S S m xy m xy m xy xy xy xy n xy n xy n xy n xy ωαωαωωαωαωωαωαωαω（2.26）

式中，Hamming 和Hanning 窗的系数α分别为0.23和0.25。

三、程序实现

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船体艏摇谱分析图

四、结论分析

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