谱分析-相关函数法

1、海浪谱分析相关函数法一、 基本概念已经提出的海浪频谱很多，其中大部分是由观测到的波要素连同某些假定推导出来的，大部分则利用定点波面记录通过特殊的谱分析方法得到。后一方法是目前得到海浪谱的主要手段。在固定点连续记录到波面，通常认为它是弱平稳的过程，其相关函数为： （1.1）由已有理论可知此过程的单侧谱为 （1.2）假定海浪为具有各态历经性的平稳随机过程，可利用过程中的现实（一次波面记录）的离散值计算相关函数 （1.3）式中，N为样本容量；为乘积的个数。由此相关函数并参照式（1.2）可得谱的估计值为 （1.4）另一方面，我们定义谱密度函数 （1.5）对于离散值， （1.6）代入式（1.5），可得

2、（1.7）当时，上式变为 （1.8）而 （1.9）式（1.8）右侧称为周期图，它可通过对样本实行离散傅里叶变化得到。因此估计谱通常有两种途径，其一通过相关函数，其二通过周期图。在每一途径中又可采用不同的方法。不管用何法，都要对实测记录取离散值，并进行中心化处理。采样间隔的选取，非常重要。在图（1.1）中，细线代表谱中圆频率为的组成波，今按时间间隔读取波面值，连接这些离散值得粗线所示的圆频率为的波动。容易推想，许多高频率的波动可表现为同一低频的波动。设定义圆频率 （1.10）则可证明频率的波动，由于离散化的结果均变现为频率的波动。设都是整数，则 （1.11）图1.1 折叠频率说明图这意味着利用波

3、面记录离散值进行谱估计时，将使谱于频率间隔内的能量都全部叠加到间隔内，这不仅使谱值的分布范围缩窄至频率范围内，而且得到的谱值不是真实的。称为Nyquist圆频率或折叠圆频率，而 （1.12）叫做Nyquist频率或折叠频率（Folding frequency）。因此，估计谱的分布范围取决于的值。事实上，海浪谱通常集中在较窄的频域内，通常可从记录中选择最短波的频率作为谱估计的频率上界，把大于的高频部分切去不计，故叫做切割频率（Cut-off frequency）。选取时应使 （1.13）海浪可视为随机过程，但可供使用的定点波形记录具有下列局限性：记录次数是有限的；记录长度是有限的；计算时使用以一

4、定间隔读取的波面数值。理论上可证明，即使计算本身无误差，由此记录得到的结果并非谱的真值，而是对真值的某种估计，故为了说明估计值接近真值的程度，尚需利用一些统计上的特征量（偏度、方差、置信度等）加以描述。二、 由相关函数估计频谱2.1 计算相关函数设采样时距为，则式（1.3）的相关函数可改写成 （2.1）这样便得到的m+1个值，它们等间隔地分布着，并分别位于。2.2 估算谱粗值将式（2.1）代入前文公式以数值积分计算谱值。由于折叠的影响，谱值系在范围内进行计算。等间隔地取m+1个频率。我们的目的就是计算谱于这些频率所具有的值，令代表频率对应的粗谱值，得 （2.2）如数值积分中采用梯形公式，谱值为

5、 （2.3）此处采用的频率间隔为 （2.4）代入式（2.3），得 （2.5）2.3 谱的平滑以上估计出的谱值是不精确的，由它们给出的谱曲线参差不齐。因为样本容量N是有限的，故式（2.1）计算相关函数时，对于小的，乘积的个数越多，从而的值较可靠；而对于大的，的可靠性较差。为了改进精确度，可令不同的具有不同的权。这种权函数的形式很多，其中一种常用的为 （2.6）此权函数乘以式（2.2）中的，最后可得谱值 （2.7）对于两个端点频率，可取 （2.8）另一种常用的权函数为 （2.9）同上，由此权函数可得谱值： （2.10）对于两个端点频率，系数均为0.5.由上可见，此处所谓改进谱的质量，实际上是采用特

6、定的系数，对谱的粗值进行平滑，而权函数称为延时窗，前者叫做哈明（Hamming）窗，后者叫做哈宁（Hanning）窗。此权函数的傅里叶变换 （2.11）称为谱窗。图2.1为由数值模拟方法得到的相关函数和谱的示例，估计谱时采用了哈明窗。图2.1 由数值模拟得到的相关函数和谱2.4 确定置信度设都是符合标准正态分布的随机变量，且，则的分布成分布（具有k个自由度）： （2.12）又设平滑后的谱值为，谱的真值为，已证明随机量遵从自由度为k的分布。根据估计谱值的概率分布，可利用置信界限来表示估计值的可靠性。设给定置信水平为（以%表示），则可由分布确定上界和下界，使估计值落入此界限内的概率为。设以代表，为

7、了确定置信界限，我们依分布求出两个正数a和b，使 （2.13）从而式（2.13）中的a和b可利用已编制的分布概率表查得。上式可写为或 （2.14）上式表明，对于给定的置信水平，置信上上下界分别为 （2.15）以上为频率对应的置信界限，各频率的置信界限构成置信带。式（2.14）中的自由度k，因使用的延时窗而异。在海浪谱估计中，常使用Tukey导出的结果，即 （2.16）采用与上式最接近的整数。2.5 参量的选取谱估计涉及到一系列参量的选取，如样本长度，最大推移乘积个数m等，它们的选取直接影响到谱的质量。理论上，样本越长，统计特征值越稳定，但其计算工作量较大。另外，海浪并不是严格平稳的，记录时间太

8、长可能会使平稳性受到影响。一般对小的波浪，样本可短些；周期大的波浪，样本宜长些，通常可取1020min，波数不宜于100个。Arhan在北海北部连续测波22h45min，把记录分成80段，每段长17min，分析得各段的有效波高随时间的变化曲线，发现在波浪迅速成长时，如时段大于17min，波浪不平稳。Haver在挪威沿海收集了384组波浪资料，每段长约17min，发现有1530%的时段不平稳。故认为时段长于17min是不可取的。如上文所述，的选取必须充分小，以避免折叠影响，同时过大时失去信息过多，会使估计得谱变形。但如过小，增加样本容量，且所得序列数据相关性增大。通常值应满足式（1.13）即 （

9、2.12）为切割频率。合田良实建议取，即在一个有效波周期内采用1020个样。也有人建议取。如由此估计所得谱在附近的谱值明显地大于零，应缩小。反之，如远小于处的谱值已接近于零，可适当地加大。最大推移乘积个数m对谱估计结果有相当影响。估计谱的质量要用估计值相对于真值的均方误差来度量，容易证明： （2.17）此处偏差 （2.18）方差 （2.19）表示估计值偏离真值得程度；表示估计值的离散程度，两者的值越小，表明谱估计结果的质量越高。一般来说，正比于，正比于。故m增大时，即偏度减小，但即方差增大。另一方面，m增大时，由式（2.16）知自由度k减小，使估计谱的置信带加宽。同时由式（2.4）知减小，即谱

10、的分辨率提高。能否通过式（2.17）对m值进行优选，国外学者曾对此做过研究，但尚无可供使用的结果，只能根据经验或通过比较选取。图2.2为谱估计结果随m的变化情况。图中是由谱矩求得的平均波高。m太小时，因点据稀少（谱估计的频率间距），不能确切地显示谱形，峰值偏小。随着m增大谱形显露。m再增大，谱形振动，谱峰尖突，故可由此选定适合的m值。另一种考虑是由各态历经性，当增大到一定值是，趋近于零，可选相关函数近于零时相应的m值。第三，不同m值对应的谱的特征值也不同，可选择谱特征值（谱宽参数）趋于稳定时所对应的m值。图2.2 谱估计结果随m值得变化 上述通过估计结果来选取m值比较麻烦。国家海洋局海洋调查规

11、范建议选取。笔者由数模结果建议选取；。2.6 由互相关函数估计互谱研究任意两个随机过程之间的相互关系时，需计算其互相关函数和互谱。如前文所述，经中心化后任意两个随机过程和的互相关函数为 （2.20）令 （2.21） （2.22）则可得各粗谱值：粗同相谱 （2.23）粗转相谱 （2.24）粗互谱 （2.25）式中，。采用Hamming窗或Hanning窗进行平滑，得 （2.26）式中，Hamming和Hanning窗的系数分别为0.23和0.25。3、 程序实现3.1 程序一：随机信号的自相关函数、协方差函数的实现% 随机信号的建立a=randn(2000,1); % 随机信号wc=0.45,0

12、.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,

13、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 程序二：自相关函数法、傅立叶变换法及最

14、大熵估计法% 傅立叶变换、自相关函数法以及最大熵估计法对离散数据进行谱分析% 离散数据的建立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(

15、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);% en

16、d图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(

17、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(横

18、荡);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

19、(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