王言英，徐繼文，赫亮

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

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

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

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

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

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

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

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

圖1 實(shí)測409個子樣的能量譜密度函數(shù)Fig.1 The energy spectral density functions for real data with 409 samples

圖2 實(shí)測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 實(shí)測譜與譜展式給出子樣特征值的比較Table 2 Comparison of eigenvalue between measured andspactral expansion

目前在海洋工程領(lǐng)域，由于海浪觀測資料的貧乏，工程設(shè)計中大都使用JONSWAP譜展式.為了保證結(jié)構(gòu)物營運(yùn)的可靠性和經(jīng)濟(jì)性，充分利用目標(biāo)海域的風(fēng)和浪的觀測資料，以FFT算法對JONSWAP譜展式中參變量進(jìn)行校驗(yàn)確認(rèn)是必要的.

2 最大熵變換(MEM)算法

熵在信息論中是衡量隨機(jī)事件不定性程度的量，Burg首先提出了最大熵法計算頻譜，其主要思路是把相關(guān)函數(shù)外推至無窮后再進(jìn)行頻域變換.

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

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

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

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

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

對歐洲北海North Alwyn平臺給出的風(fēng)浪實(shí)測資料，分別以FFT和MEM算法進(jìn)行了譜分析計算.為避免頻率折疊的影響，通常從波面記錄中選取最短波的頻率，大于此頻率的波的能量小到可以忽略，取大于此頻率的fc為截止頻率.此外，根據(jù)取樣原理當(dāng) Δt≤1/2fc時，可用取樣序列 x(jΔt)(j=1，2，…，N)，恢復(fù)原來的連續(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

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

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

3 小波變換(WT)算法

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

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

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

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

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

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

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

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

4 自舉變換(BT)算法

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

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

有些情況下只能以容量較小的子樣為依據(jù)從事信號的譜分析，Welch方法是一種調(diào)制的周期圖法，將時間序列劃分為重疊或非重疊的區(qū)段，用該方法必須在頻率分辨率和譜方差之間做出權(quán)衡.較長的區(qū)段長度或較少的數(shù)據(jù)點(diǎn)的重合都會提高譜的分辨率，而較少的區(qū)段則會降低其方差.但是對于固定的波浪數(shù)據(jù)長度，區(qū)段的數(shù)量同區(qū)段長度和重合點(diǎn)數(shù)總是呈逆反關(guān)系的［6］.為了檢驗(yàn)分析的魯棒性，從總的記錄中隨機(jī)地選擇了2個時間序列，見圖9.

圖9 兩組波面升高序列Fig.9 Wave elevation values of two sets

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

5 結(jié)束語

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