譜分析-相關(guān)函數(shù)法

1、海浪譜分析相關(guān)函數(shù)法一、 基本概念已經(jīng)提出的海浪頻譜很多，其中大部分是由觀測到的波要素連同某些假定推導(dǎo)出來的，大部分則利用定點(diǎn)波面記錄通過特殊的譜分析方法得到。后一方法是目前得到海浪譜的主要手段。在固定點(diǎn)連續(xù)記錄到波面，通常認(rèn)為它是弱平穩(wěn)的過程，其相關(guān)函數(shù)為： （1.1）由已有理論可知此過程的單側(cè)譜為 （1.2）假定海浪為具有各態(tài)歷經(jīng)性的平穩(wěn)隨機(jī)過程，可利用過程中的現(xiàn)實(shí)（一次波面記錄）的離散值計(jì)算相關(guān)函數(shù) （1.3）式中，N為樣本容量；為乘積的個(gè)數(shù)。由此相關(guān)函數(shù)并參照式（1.2）可得譜的估計(jì)值為 （1.4）另一方面，我們定義譜密度函數(shù) （1.5）對(duì)于離散值， （1.6）代入式（1.5），可得

2、（1.7）當(dāng)時(shí)，上式變?yōu)?（1.8）而 （1.9）式（1.8）右側(cè)稱為周期圖，它可通過對(duì)樣本實(shí)行離散傅里葉變化得到。因此估計(jì)譜通常有兩種途徑，其一通過相關(guān)函數(shù)，其二通過周期圖。在每一途徑中又可采用不同的方法。不管用何法，都要對(duì)實(shí)測記錄取離散值，并進(jìn)行中心化處理。采樣間隔的選取，非常重要。在圖（1.1）中，細(xì)線代表譜中圓頻率為的組成波，今按時(shí)間間隔讀取波面值，連接這些離散值得粗線所示的圓頻率為的波動(dòng)。容易推想，許多高頻率的波動(dòng)可表現(xiàn)為同一低頻的波動(dòng)。設(shè)定義圓頻率 （1.10）則可證明頻率的波動(dòng)，由于離散化的結(jié)果均變現(xiàn)為頻率的波動(dòng)。設(shè)都是整數(shù)，則 （1.11）圖1.1 折疊頻率說明圖這意味著利用波

3、面記錄離散值進(jìn)行譜估計(jì)時(shí)，將使譜于頻率間隔內(nèi)的能量都全部疊加到間隔內(nèi)，這不僅使譜值的分布范圍縮窄至頻率范圍內(nèi)，而且得到的譜值不是真實(shí)的。稱為Nyquist圓頻率或折疊圓頻率，而 （1.12）叫做Nyquist頻率或折疊頻率（Folding frequency）。因此，估計(jì)譜的分布范圍取決于的值。事實(shí)上，海浪譜通常集中在較窄的頻域內(nèi)，通?？蓮挠涗浿羞x擇最短波的頻率作為譜估計(jì)的頻率上界，把大于的高頻部分切去不計(jì)，故叫做切割頻率（Cut-off frequency）。選取時(shí)應(yīng)使 （1.13）海浪可視為隨機(jī)過程，但可供使用的定點(diǎn)波形記錄具有下列局限性：記錄次數(shù)是有限的；記錄長度是有限的；計(jì)算時(shí)使用以一

4、定間隔讀取的波面數(shù)值。理論上可證明，即使計(jì)算本身無誤差，由此記錄得到的結(jié)果并非譜的真值，而是對(duì)真值的某種估計(jì)，故為了說明估計(jì)值接近真值的程度，尚需利用一些統(tǒng)計(jì)上的特征量（偏度、方差、置信度等）加以描述。二、 由相關(guān)函數(shù)估計(jì)頻譜2.1 計(jì)算相關(guān)函數(shù)設(shè)采樣時(shí)距為，則式（1.3）的相關(guān)函數(shù)可改寫成 （2.1）這樣便得到的m+1個(gè)值，它們等間隔地分布著，并分別位于。2.2 估算譜粗值將式（2.1）代入前文公式以數(shù)值積分計(jì)算譜值。由于折疊的影響，譜值系在范圍內(nèi)進(jìn)行計(jì)算。等間隔地取m+1個(gè)頻率。我們的目的就是計(jì)算譜于這些頻率所具有的值，令代表頻率對(duì)應(yīng)的粗譜值，得 （2.2）如數(shù)值積分中采用梯形公式，譜值為

5、 （2.3）此處采用的頻率間隔為 （2.4）代入式（2.3），得 （2.5）2.3 譜的平滑以上估計(jì)出的譜值是不精確的，由它們給出的譜曲線參差不齊。因?yàn)闃颖救萘縉是有限的，故式（2.1）計(jì)算相關(guān)函數(shù)時(shí)，對(duì)于小的，乘積的個(gè)數(shù)越多，從而的值較可靠；而對(duì)于大的，的可靠性較差。為了改進(jìn)精確度，可令不同的具有不同的權(quán)。這種權(quán)函數(shù)的形式很多，其中一種常用的為 （2.6）此權(quán)函數(shù)乘以式（2.2）中的，最后可得譜值 （2.7）對(duì)于兩個(gè)端點(diǎn)頻率，可取 （2.8）另一種常用的權(quán)函數(shù)為 （2.9）同上，由此權(quán)函數(shù)可得譜值： （2.10）對(duì)于兩個(gè)端點(diǎn)頻率，系數(shù)均為0.5.由上可見，此處所謂改進(jìn)譜的質(zhì)量，實(shí)際上是采用特

6、定的系數(shù)，對(duì)譜的粗值進(jìn)行平滑，而權(quán)函數(shù)稱為延時(shí)窗，前者叫做哈明（Hamming）窗，后者叫做哈寧（Hanning）窗。此權(quán)函數(shù)的傅里葉變換 （2.11）稱為譜窗。圖2.1為由數(shù)值模擬方法得到的相關(guān)函數(shù)和譜的示例，估計(jì)譜時(shí)采用了哈明窗。圖2.1 由數(shù)值模擬得到的相關(guān)函數(shù)和譜2.4 確定置信度設(shè)都是符合標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量，且，則的分布成分布（具有k個(gè)自由度）： （2.12）又設(shè)平滑后的譜值為，譜的真值為，已證明隨機(jī)量遵從自由度為k的分布。根據(jù)估計(jì)譜值的概率分布，可利用置信界限來表示估計(jì)值的可靠性。設(shè)給定置信水平為（以%表示），則可由分布確定上界和下界，使估計(jì)值落入此界限內(nèi)的概率為。設(shè)以代表，為

7、了確定置信界限，我們依分布求出兩個(gè)正數(shù)a和b，使 （2.13）從而式（2.13）中的a和b可利用已編制的分布概率表查得。上式可寫為或 （2.14）上式表明，對(duì)于給定的置信水平，置信上上下界分別為 （2.15）以上為頻率對(duì)應(yīng)的置信界限，各頻率的置信界限構(gòu)成置信帶。式（2.14）中的自由度k，因使用的延時(shí)窗而異。在海浪譜估計(jì)中，常使用Tukey導(dǎo)出的結(jié)果，即 （2.16）采用與上式最接近的整數(shù)。2.5 參量的選取譜估計(jì)涉及到一系列參量的選取，如樣本長度，最大推移乘積個(gè)數(shù)m等，它們的選取直接影響到譜的質(zhì)量。理論上，樣本越長，統(tǒng)計(jì)特征值越穩(wěn)定，但其計(jì)算工作量較大。另外，海浪并不是嚴(yán)格平穩(wěn)的，記錄時(shí)間太

8、長可能會(huì)使平穩(wěn)性受到影響。一般對(duì)小的波浪，樣本可短些；周期大的波浪，樣本宜長些，通常可取1020min，波數(shù)不宜于100個(gè)。Arhan在北海北部連續(xù)測波22h45min，把記錄分成80段，每段長17min，分析得各段的有效波高隨時(shí)間的變化曲線，發(fā)現(xiàn)在波浪迅速成長時(shí)，如時(shí)段大于17min，波浪不平穩(wěn)。Haver在挪威沿海收集了384組波浪資料，每段長約17min，發(fā)現(xiàn)有1530%的時(shí)段不平穩(wěn)。故認(rèn)為時(shí)段長于17min是不可取的。如上文所述，的選取必須充分小，以避免折疊影響，同時(shí)過大時(shí)失去信息過多，會(huì)使估計(jì)得譜變形。但如過小，增加樣本容量，且所得序列數(shù)據(jù)相關(guān)性增大。通常值應(yīng)滿足式（1.13）即 （

9、2.12）為切割頻率。合田良實(shí)建議取，即在一個(gè)有效波周期內(nèi)采用1020個(gè)樣。也有人建議取。如由此估計(jì)所得譜在附近的譜值明顯地大于零，應(yīng)縮小。反之，如遠(yuǎn)小于處的譜值已接近于零，可適當(dāng)?shù)丶哟?。最大推移乘積個(gè)數(shù)m對(duì)譜估計(jì)結(jié)果有相當(dāng)影響。估計(jì)譜的質(zhì)量要用估計(jì)值相對(duì)于真值的均方誤差來度量，容易證明： （2.17）此處偏差 （2.18）方差 （2.19）表示估計(jì)值偏離真值得程度；表示估計(jì)值的離散程度，兩者的值越小，表明譜估計(jì)結(jié)果的質(zhì)量越高。一般來說，正比于，正比于。故m增大時(shí)，即偏度減小，但即方差增大。另一方面，m增大時(shí)，由式（2.16）知自由度k減小，使估計(jì)譜的置信帶加寬。同時(shí)由式（2.4）知減小，即譜

10、的分辨率提高。能否通過式（2.17）對(duì)m值進(jìn)行優(yōu)選，國外學(xué)者曾對(duì)此做過研究，但尚無可供使用的結(jié)果，只能根據(jù)經(jīng)驗(yàn)或通過比較選取。圖2.2為譜估計(jì)結(jié)果隨m的變化情況。圖中是由譜矩求得的平均波高。m太小時(shí)，因點(diǎn)據(jù)稀少（譜估計(jì)的頻率間距），不能確切地顯示譜形，峰值偏小。隨著m增大譜形顯露。m再增大，譜形振動(dòng)，譜峰尖突，故可由此選定適合的m值。另一種考慮是由各態(tài)歷經(jīng)性，當(dāng)增大到一定值是，趨近于零，可選相關(guān)函數(shù)近于零時(shí)相應(yīng)的m值。第三，不同m值對(duì)應(yīng)的譜的特征值也不同，可選擇譜特征值（譜寬參數(shù)）趨于穩(wěn)定時(shí)所對(duì)應(yīng)的m值。圖2.2 譜估計(jì)結(jié)果隨m值得變化 上述通過估計(jì)結(jié)果來選取m值比較麻煩。國家海洋局海洋調(diào)查規(guī)

11、范建議選取。筆者由數(shù)模結(jié)果建議選??；。2.6 由互相關(guān)函數(shù)估計(jì)互譜研究任意兩個(gè)隨機(jī)過程之間的相互關(guān)系時(shí)，需計(jì)算其互相關(guān)函數(shù)和互譜。如前文所述，經(jīng)中心化后任意兩個(gè)隨機(jī)過程和的互相關(guān)函數(shù)為 （2.20）令 （2.21） （2.22）則可得各粗譜值：粗同相譜 （2.23）粗轉(zhuǎn)相譜 （2.24）粗互譜 （2.25）式中，。采用Hamming窗或Hanning窗進(jìn)行平滑，得 （2.26）式中，Hamming和Hanning窗的系數(shù)分別為0.23和0.25。3、 程序?qū)崿F(xiàn)3.1 程序一：隨機(jī)信號(hào)的自相關(guān)函數(shù)、協(xié)方差函數(shù)的實(shí)現(xiàn)% 隨機(jī)信號(hào)的建立a=randn(2000,1); % 隨機(jī)信號(hào)wc=0.45,0

12、.65;N=79;% 基于窗函數(shù)的fir濾波器window=blackman(N+1); % 調(diào)用格式：w=blackman(n),根據(jù)長度n產(chǎn)生一個(gè)布拉克曼窗w；h=fir1(N,wc,window); % 調(diào)用格式：fir1(n,Wn,ftype,Window),n 為階數(shù),Wn是截止頻率；x=filter(h,1,a); % Y,Zf = filter(B,A,X,Zi)，輸入X為濾波前序列，Y為濾波結(jié)果序列，B/A 提供濾波器系數(shù)，B為分子， A為分母% 時(shí)域信號(hào)subplot(2,2,1),plot(x),title(bf時(shí)域信號(hào))grid on% 自相關(guān)函數(shù)c,n=xcorr(x,

13、10,coeff);subplot(2,2,2),stem(n,c,filled),title(bf自相關(guān)函數(shù))grid on% 協(xié)方差函數(shù)b,m=xcov(x,10,coeff);subplot(2,2,3),stem(m,b,filled),title(bf協(xié)方差函數(shù));grid on% 概率密度函數(shù)subplot(2,2,4),pwelch(x,33,32,500),title(bf概率密度函數(shù)); % Pxx,f=pwelch(x,window,noverlap,nfft,fs)；grid on% end圖3.1 隨機(jī)信號(hào)的自相關(guān)函數(shù)的實(shí)現(xiàn)3.2 程序二：自相關(guān)函數(shù)法、傅立葉變換法及最

14、大熵估計(jì)法% 傅立葉變換、自相關(guān)函數(shù)法以及最大熵估計(jì)法對(duì)離散數(shù)據(jù)進(jìn)行譜分析% 離散數(shù)據(jù)的建立t=0:0.001:1;e=randn(size(t); % 噪聲data=4*sin(25*2*pi*t)+4*sin(50*2*pi*t)+e;% 離散數(shù)據(jù)的截取subplot(2,2,1);plot(1000*t(1:50),data(1:50);xlabel(ittime(mm);title(bf一元時(shí)間序列直觀圖);% 自相關(guān)函數(shù)法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自相關(guān)函數(shù)法);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）);% 最大熵估計(jì)法subplot(2,2,4);Fs=500;NFFT=1024;pyulear(data,20,NFFT,Fs);% en

16、d圖3.2 離散數(shù)據(jù)的譜分析方法對(duì)比圖3.3 程序三：隨機(jī)波的譜分析讀取試驗(yàn)數(shù)據(jù)：data=load(wave.txt);圖3.3 隨機(jī)波幅值變化圖圖3.4 隨機(jī)波譜分析圖3.4 程序四：經(jīng)典譜分析算法基于中海油、能威設(shè)計(jì)有限公司和中國海洋大學(xué)工程學(xué)院合作的雙船浮托整體棄置模型試驗(yàn)的船體六自由度數(shù)據(jù)，利用自相關(guān)系數(shù)法、傅立葉變換法及最大熵估計(jì)法進(jìn)行譜分析。% 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 船體六自由度圖% 相關(guān)函數(shù)法、傅立葉變換法及最大熵法進(jìn)行譜分析% sway 橫蕩subplot(2,2,1);plot(disp_w_r(:,1);xlabel(itTime);ylabel(itAmplitude);title(bfsway(橫蕩);grid on% 自相關(guān)函數(shù)法Fs=1000;NFFT=1024;Cx=xcorr(disp_w_r(:,1),unbiased