海浪能量譜密度函數(shù)計算方法

王言英，徐繼文，赫亮

(大連理工大學船舶工程學院，遼寧大連116024)

浮式海洋結構物的許多力學性能都對海洋環(huán)境運動因素產生響應.在諸多海洋環(huán)境因素中起主導作用的是海浪，所指的是風生浪，為重力波.波浪運動屬于隨機過程，通常被定義為平穩(wěn)的各態(tài)歷經的隨機過程.工程計算與設計不僅關注那些周期性運動性能的幅值，還關注其運動能量在頻域的分布.海浪能量譜密度函數(shù)是浮式海洋結構物的許多力學性能預報的重要輸入信息(輸入譜)，計入結構物固有的頻率響應函數(shù)，可以得到相關性能的能量譜密度函數(shù)(輸出譜)，該譜的各階矩對應各種力學性能［1］.

數(shù)學上根據(jù)Wiener-Khintchine定理，譜函數(shù)可以由隨機過程時歷的自相關函數(shù)的Fourier變換得到，這就是目前廣泛流行的快速傅里葉變換(FFT)算法的原理.在船舶與海洋工程領域，ITTC篩選了一批廣泛采用的海浪能量譜密度函數(shù)的半理論半經驗的所謂譜展式［2］.隨機過程能量譜密度函數(shù)計算方法的進一步研究是工程設計發(fā)展的需要.就海浪譜而言，上述的譜展式大多為可用于各海域的通式，并未考慮特定海域環(huán)境的特殊性.實際上鑒于海域的位置不同，周圍遮蔽陸地環(huán)境的不同以及水深的不同，海浪運動的能量大小及其頻率特性的差異是顯著的.如所周知，在工程設計中使用的譜展式的能量直接關系到結構物的可靠性和經濟性.為此，根據(jù)結構物營運海域的海洋環(huán)境資料，建立或修正現(xiàn)用的海浪譜展式，以最大限度地符合當?shù)睾Ｓ虻沫h(huán)境因素是必要的.另外，在必須自行確定環(huán)境或結構物的某些力學性能時，所擁有的過程子樣的容量常常是受到限制或有限的.所以，尋求更多的計算方法，有效地利用小子樣信息或者盡量挖掘有限容量子樣的有益信息，以獲得有足夠置信度的譜函數(shù)是工程上需要的.

本文以歐洲北海North Alywn平臺對于一次風暴的觀測記錄為樣本，以FFT算法，對JONSWAP譜展式的譜峰升高因子和峰頻控制因子進行了修正［3-4］;分別建立了譜分析的最大熵變換(MET)算法，小波變換(WT)算法［5］和自舉變換(BT)算法［6］，并同F(xiàn)FT算法給出的結果進行了比較.

1 快速傅里葉變換(FFT)算法

為驗證JONSWAP譜展式在歐洲北海海域的有效性，根據(jù)1977年11月北海North Alywn平臺提供的一次風暴的記錄7，應用FFT算法進行了能量譜密度函數(shù)分析計算，并同ITTC推薦的JONSWAP譜展式［2］進行了比較分析.該記錄共包括409個子樣，每個子樣長20 min，采樣間隔為0.2 s，有6 000個海面瞬時升高值的記錄.

409個子樣的能量譜密度函數(shù)如圖1所示，各個譜的零階、二階與四階矩載于圖 2.比對JONSWAP譜展式，按最佳擬合得到的譜展式中的參變量的數(shù)值見表1.

根據(jù)實測譜和JONSWAP譜展式分別可以計算得到峰頻海況的能量譜密度函數(shù)，進一步計算得到其各階矩以及相應的統(tǒng)計特征值，諸如有義波高HS，譜峰圓頻率 ωp，平均周期等.表2給出的兩組分別由FFT和JONSWAP得到的數(shù)據(jù)表明兩者符合較好，應用JONSWAP譜展式可以獲得具有工程精度的統(tǒng)計特征預報值.鑒于JONSWAP譜展式本來就源于歐洲北海開發(fā)聯(lián)合計劃，同一海域的后報海浪譜分析結果，同譜展式給出的相應結果的符合是令人信服的.

圖1 實測409個子樣的能量譜密度函數(shù)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首先提出了最大熵法計算頻譜，其主要思路是把相關函數(shù)外推至無窮后再進行頻域變換.

有N個事件，若每個事件的出現(xiàn)概率為pi(i=1，2，…，n)，且n個事件相互獨立，則單位時間間隔內的平均信息量，熵H為

其中，S(ω)為概率密度函數(shù)，為正態(tài)分布的連續(xù)型隨機變量的能量譜密度函數(shù).由信息論知，當隨機事件是以等概率可能性出現(xiàn)時其熵值達到最大.在連續(xù)型隨機變量呈標準正態(tài)分布的概率密度時，其熵達到極大值.根據(jù)變分原理和自回歸分析理論，可得x的譜密度函數(shù)的表達式為

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

在估計最大熵譜時，模型階數(shù)的選擇是一個關鍵問題.常用的判階準則有信息論準則(AIC)，自回歸準則(CTA)和最終預測誤差準則(FPE)，數(shù)值試驗結果表明，當信噪比較高時上述3種方法確定的階數(shù)M基本一致.根據(jù)FPE準則，對于中心化的序列有

取(PEM)M達到最小時的M值作為最大熵譜的最佳階數(shù);也有采用容量比定階數(shù)，即M=N/20.

對歐洲北海North Alwyn平臺給出的風浪實測資料，分別以FFT和MEM算法進行了譜分析計算.為避免頻率折疊的影響，通常從波面記錄中選取最短波的頻率，大于此頻率的波的能量小到可以忽略，取大于此頻率的fc為截止頻率.此外，根據(jù)取樣原理當 Δt≤1/2fc時，可用取樣序列 x(jΔt)(j=1，2，…，N)，恢復原來的連續(xù)函數(shù) X(t)，使 fc≤1/2Δt.圖3為以不同采樣容量N和階數(shù)M計算得到的MEM估計譜的譜型;圖4為以不同采樣容量N計算得到的FFT估計譜的譜型.表3給出的是考慮不同采樣容量N和階數(shù)M得到的譜特征參數(shù)均方根波高Hrms，有義波高 HS，平均過零周期 TZ，平均周期T1，峰周期 TP，峰頻 ωP，譜寬系數(shù) ξ，以及譜零階、二階與四階矩 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

采用某一時段波面記錄的樣本序列直接計算譜值，這等于使用了矩形窗.因為矩形窗對應的譜窗其旁瓣效應較大，計算結果會出現(xiàn)虛假的譜峰或使譜線上下波動，這就是所謂的譜泄漏.為了減少這種泄漏，在使用FFT法計算時先采用余弦半鐘式數(shù)據(jù)窗.由表3可以看出，用FFT法計算的結果與MEM的計算結果相比基本是一致的，由譜的特征值與譜矩而得到的特征波要素均較吻合.就此算例而言，當N=60 000和M=36時，F(xiàn)FT和MEM算法給出的譜特征參數(shù)基本一致，而且當子樣容量減少到6 000的1/2～1/4時，選擇適當?shù)碾A數(shù)仍然可以保持足夠的譜特征參數(shù)精度.

表3 快速傅里葉與最大熵譜估計特征參數(shù)計算結果的比較Table 3 Comparison of statistical characteristic values given by FFT and MEM algorithms

3 小波變換(WT)算法

隨機信號x(t)的小波變換為

ψab為變換與擴展小波，由小波母函數(shù)ψ(t)得到:

式中:fb為帶寬參數(shù)，fc為小波中心頻率.其逆變換為

其中，系數(shù)C由式(8)確定.

小波譜分析的基本思想是將一時間序列的方差分解成許多分量，其中每一個分量都是在一特定時間的一個特定尺度.在處理有限長度時間序列時，考慮到傅里葉變換內在的有限性，可以應用小波譜分析方法.小波滿足的能量守恒方程為

那么，頻率譜密度函數(shù)S(f)可以從下式推導得出

當采用波數(shù)k=2的Morlet小波時，中心頻率fc=1.0.

圖5 子樣2 000～6 000的JONSWAP譜與常規(guī)譜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，F(xiàn)FT 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譜與常規(guī)譜的比較，圖6為同上子樣容量的情況下JONSWAP譜與小波譜的比較，圖7為子樣容量為1 000情況下JONSWAP譜，快速傅里葉譜與小波譜的比較，圖8的2種算法在不同子樣容量情況下譜給出的有義波高結果的相對偏差.所謂常規(guī)的譜分析方法是基于Wiener-Kintchine定理的通過時間序列自相關函數(shù)的傅里葉變換得到的譜函數(shù)［8］.關于有義波高相對偏差計算分別采用實測子樣和根據(jù)統(tǒng)計子樣特征模擬的時間序列得到的.

同F(xiàn)FT算法相比小波算法更適合于小子樣的譜分析，而且對于時間-頻率域的信號處理具有良好特性.在船舶與海洋工程領域，基于有限容量或小子樣的波浪和結構物力學性能響應的信號處理，海浪的時間-頻率譜的計算，以及探討頻率隨時間的變化趨勢的計算，小波分析算法是一個值得推崇的方法.

4 自舉變換(BT)算法

自舉(Bootstrap)的意思是毋需借助外界的幫助而依靠自身的主動性來推進和發(fā)展.近年來，自舉算法已獲得廣泛用于估算標準差、置信區(qū)間、偏差與預報誤差.在一些統(tǒng)計應用中，興趣集中在根據(jù)一隨機子樣的概率分布估算統(tǒng)計量，該概率分布并不確切知道，而且用以估算那個量的統(tǒng)計采樣分布也并不確切知道.自舉方法允許人們通過在計算機上多次模擬原始試驗，實現(xiàn)近似概率計算.原始的自舉算法要求依據(jù)獨立的和相等的數(shù)據(jù)重構子樣.在這種情況下，人們可以從數(shù)據(jù)中通過隨機抽樣更換其位置，構建出人為的重復子樣.對于時間序列分析，原始方法不能捕捉到附近觀測的相關結構，新的自舉方法可以做到這一點.對于有模型逼近，相關結構被模擬成少數(shù)已知參數(shù)和具有獨立的誤差.鑒于波浪產生的機理不是很清楚與明確的，在無模型自舉方法中多采用區(qū)間自舉逼近［9］.

應用區(qū)間自舉逼近時，觀測子樣被分成一些區(qū)間，旨在從原始數(shù)據(jù)序列中捕捉相關因素.本文采用的是移動區(qū)間自舉逼近方法.應用移動區(qū)間自舉逼近選擇一最佳區(qū)間長度L和移動尺度M是必要的.現(xiàn)假定一觀測序列 X1，X2，…，XN，系嚴格取自固定的觀測序列{Xn，n∈Z}，典型的移動區(qū)間自舉逼近算法是，選擇L和M，區(qū)間數(shù)為B=fix((N-L)/M)+1，對于 m=1，2，…，B，有

有些情況下只能以容量較小的子樣為依據(jù)從事信號的譜分析，Welch方法是一種調制的周期圖法，將時間序列劃分為重疊或非重疊的區(qū)段，用該方法必須在頻率分辨率和譜方差之間做出權衡.較長的區(qū)段長度或較少的數(shù)據(jù)點的重合都會提高譜的分辨率，而較少的區(qū)段則會降低其方差.但是對于固定的波浪數(shù)據(jù)長度，區(qū)段的數(shù)量同區(qū)段長度和重合點數(shù)總是呈逆反關系的［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算法的使用，區(qū)間的子樣容量總是取為2的冪，這個數(shù)總是主宰著譜的分辨率，這是指能夠分辨2個相鄰譜峰和其他頻率的譜峰.以下的計算其區(qū)間長度均取為L=512和不同的尺度M.

為了解決短數(shù)據(jù)長度的譜分辨率同方差的矛盾，必須增加波浪時間數(shù)據(jù)長度，引進一個區(qū)段自舉方法.這個方法是對Welch方法的修正，可以稱之為Welch方法或者區(qū)段自舉方法.在各個區(qū)段所用的自舉方法是相同的，新的自舉子樣是從B區(qū)段分布的數(shù)據(jù)點隨機地移至Q區(qū)段的數(shù)據(jù)組成的，自舉的Welch譜是Q區(qū)段譜的平均.為了比較的方便，區(qū)段的尺度同其移動的尺度是一樣的，在該區(qū)段應用Welch方法，再建區(qū)段數(shù)為80.計算結果見圖10，其中(a)、(b)為子樣容量為3 000的結果，(c)、(d)為子樣容量為6 000的結果.從這4張圖很難看出由于移動尺度的不同所導致的譜的差別，說明當重新取樣的數(shù)量足夠大時，譜并不隨移動尺度而變化.考慮到自舉Welch譜對移動尺度不敏感，進行了M=0的數(shù)值試驗，結果再次表明譜型同移動尺度幾乎無關，所以，對于自舉Welch方法并不需要區(qū)段數(shù)據(jù)的重疊.另外，數(shù)值試驗還表明，當數(shù)據(jù)長度小到1 000時，Welch譜已經嚴重扭曲，而自舉Welch譜仍然保持光滑.

5 結束語

對于有限時間序列信號，為了保證譜估計有足夠精度的分辨率和方差，最大熵算法(MEM)和小波變換算法(WT)是有效的方法.對于小容量時間序列信號，為了得到譜峰具有高分辨率和小方差，自舉(BT)算法是實現(xiàn)時間序列自身擴展行之有效的方法.對于新購建的時間序列從事譜分析，采用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.