基于CEEMD的船模非線性波浪力測量信號分析

王啟興，許 勇，歐勇鵬

(1. 海軍駐柳州434廠軍事代表室, 廣西 柳州 545000；2. 海軍工程大學(xué) 艦船學(xué)院, 湖北 武漢 430033)

0 引 言

船模試驗在船舶流體力學(xué)研究領(lǐng)域扮演著重要的角色，是考核理論計算方法和闡述船舶流體力學(xué)機理的關(guān)鍵。在船模的耐波性試驗中常需要測量船模在規(guī)則波中所受的波浪作用力。這種波浪力可分為兩大部分：一部分為具有零平均性質(zhì)的線性振蕩波浪力，其振蕩頻率為船模的遭遇頻率其中，及分別為入射波頻率、波數(shù)及浪向角，為船模航速），稱為一階波浪力；另一部分為非線性影響產(chǎn)生的非線性波浪力，稱為高階波浪力[1]。通常二階波浪力在高階波浪力中占主導(dǎo)，但當(dāng)船模遭遇頻率和其運動模態(tài)的共振頻率相近時，船模所受的波浪力非線性加劇，更高階的波浪力組分也會在高階波浪力中占有一定的比重。完整的二階波浪力中包括定常二階波浪力、低頻二階波浪力和倍頻二階波浪力[2]。綜上可見，船模在波浪中所受的波浪力組分非常豐富，且是非線性的，因此模型試驗中測得的波浪力信號也是非線性的。實測船模波浪力信號中往往還包含著高頻噪聲，其主要來源于測量設(shè)備中的交流電噪聲和拖車及其旋轉(zhuǎn)機械的振動噪聲[3]。測量儀器過熱時的零點漂移、儀器調(diào)零不精確等因素還可能使測量信號中含有低頻趨勢項的影響，使得測量波形發(fā)生偏移。另外，測量信號中往往還存在著隨機噪聲，這種噪聲由時間上隨機產(chǎn)生的大量起伏騷擾積累而造成的，比如試驗過程中等波時間不夠，水面上的隨機無序的破碎波對船模的累積作用就會產(chǎn)生隨機噪聲。高頻噪聲、低頻趨勢項、隨機噪聲和非線性的波浪力疊加在一起，使得實測信號不再是平穩(wěn)信號。圖1給出了某船模在規(guī)則波中頂浪航行時首搖力矩測量信號時歷曲線及對應(yīng)的幅值譜，可見該信號是非平穩(wěn)的。在頻率＜50 rad/s的頻段內(nèi)，可以清晰地觀察到信號中含有一階波浪力（與對應(yīng)）、定常2階波浪力（與=0 rad/s對應(yīng)）、低頻2階波浪力（介于0～之間）、2階倍頻波浪力（與2對應(yīng)）及更高倍頻波浪力組分，而在=315 rad/s（約為50 Hz）附近也有一定的信號組分存在，對應(yīng)于測量信號中的高頻噪聲。如何對這種含有高頻噪聲及隨機噪聲的非平穩(wěn)信號進(jìn)行分析，準(zhǔn)確提取波浪力的各階組分是本文的研究目的。

經(jīng)驗?zāi)B(tài)分解算法（Empirical Mode Decomposition，EMD）[4]集平穩(wěn)化、層次處理能力于一體，為提取船模非線性波浪力各組分提供了新思路。該方法基于信號的局部特征時間尺度，將非平穩(wěn)信號分解，得到一系列不同頻率的分量——本征模態(tài)函數(shù)（IMF）。相對于原非平穩(wěn)、非線性信號，不同時間尺度的IMF平穩(wěn)性得到很大改善。分解出最初IMF序列代表原始序列的高頻組分，隨著分解的進(jìn)行，相應(yīng)IMF的頻率變小、周期增長。這些IMF可作為原序列的一組完全或幾乎正交的展開基，并與最終的分解殘差能夠很好的重構(gòu)原時間序列[5]。船模在規(guī)則波中所受波浪力的各階組分對應(yīng)著不同頻率的分量，具備良好的局部尺度特征，從理論上而言，應(yīng)用EMD方法進(jìn)行處理，各階波浪力組分應(yīng)該能夠很好的與分解得到的IMF相對應(yīng)。但由于實測數(shù)據(jù)因隨機因素或測量儀器等影響會帶有隨機誤差，使得很多情況下EMD直接對信號分解提取的規(guī)律并不準(zhǔn)確，容易出現(xiàn)模態(tài)混疊現(xiàn)象，即相同頻率的信號分量被分解到不同的模態(tài)分量當(dāng)中，或者不同頻率的信號分量被分解到同一個模態(tài)分量當(dāng)中[5]。集合經(jīng)驗?zāi)B(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）[6]對EMD方法進(jìn)行改進(jìn)，通過加入不同隨機白噪聲來平衡測量數(shù)據(jù)中隨機誤差的影響，進(jìn)而從數(shù)據(jù)中提取更準(zhǔn)確的規(guī)律，但也因為如此，增大了EEMD分解后信號的重構(gòu)誤差。Yeh等對EEMD算法進(jìn)行了改進(jìn)，提出互補集合經(jīng)驗?zāi)B(tài)分解方法[7]，與EEMD算法不同的是所加入的白噪聲正負(fù)成對，能夠極大地提高重構(gòu)精度，同時減少迭代次數(shù)并提高運算效率。本文基于CEEMD算法，并輔以合適的其他信號分析手段，研究船模非線性波浪力測量信號各階組分的提取方法。

1 EMD和CEEMD算法原理

1.1 EMD算法原理

EMD方法的實質(zhì)是提取本征模態(tài)函數(shù)IMF的過程。IMF分量必須滿足以下2個條件：1）其極值個數(shù)和過零點數(shù)相同或最多相差一個；2）其上下包絡(luò)線關(guān)于時間軸局部對稱。分解方法如下[8]：

2）對信號中計算得到的極大（?。O值點利用3次樣條法進(jìn)行插值，并擬合成信號的上下包絡(luò)線，計算2條曲線的均值序列，求出。

1.2 CEEMD算法原理

CEEMD算法實質(zhì)上是EMD算法的一種輔助算法，其核心仍是EMD，但在進(jìn)行EMD前在原始信號中添加了正負(fù)對形式的白噪聲，其過程如下[9]：

2）對各組信號采用EMD方法進(jìn)行分解，集合中的信號都可以分解為IMF分量，其中第 i個信號的第個IMF分量為。

3）通過多組分量組合得到分解結(jié)果如式（3）所示，殘余分量記作。

2 船模非線性波浪力測量信號的處理方法

根據(jù)船模非線性波浪力測量信號的特點，本文提出如下信號處理流程：

1）低通濾波

相對于高頻噪聲的頻率而言，船模在規(guī)則波中的遭遇頻率是低頻的[10]，其頻率往往低于5 Hz。大量的試驗數(shù)據(jù)分析表明，5階以上的高階波浪力量值相對于1階、甚至2階倍頻波浪力量值而言，可以忽略，因此有實際意義的高階波浪力的頻率相對于高頻噪聲的頻率而言也是低頻的。本文采用低通濾波方法去除實測波浪力信號中的高頻噪聲。低通濾波的具體實施方法見文獻(xiàn)[3]，高頻截斷頻率取6。采用低通濾波的作用過濾高頻噪聲，減少后續(xù)CEEMD分解時的迭代次數(shù)，提高運算效率。

2）CEEMD求解

對過濾高頻噪聲的船模波浪力測量信號進(jìn)行CEEMD分解，獲得一系列由高頻到低頻排列的本征模態(tài)函數(shù)IMF。分解所得IMF中的前幾個模態(tài)往往對應(yīng)于隨機噪聲占主導(dǎo)的模態(tài)。另外，當(dāng)信號中多倍頻的高階波浪力為小量，且量值與這些隨機信號相當(dāng)時，也會與隨機噪聲耦合在一起，存在于這幾個模態(tài)當(dāng)中。因此，必須要對這些模態(tài)進(jìn)行甄別，為提取與船模非線性波浪力各組分相對應(yīng)的有效IMF提供方便。

3）有效IMF識別

隨機信號的自相關(guān)函數(shù)反映了信號與其自身在不同時間點的相似程度。對于隨機噪聲而言，其自相關(guān)函數(shù)在零點處最大，其余點處立即衰減為0，而不含噪信號自相關(guān)函數(shù)值在零點處最大，其余點隨著時間差的增大而緩慢衰減至0[11]。本文利用隨機噪聲和不含噪信號自相關(guān)函數(shù)的性質(zhì)差異來對噪聲模態(tài)和有效模態(tài)進(jìn)行識別。步驟如下：分別求取各IMF的歸一化自相關(guān)函數(shù)及其方差，設(shè)定自相關(guān)函數(shù)方差的閾值，當(dāng)模態(tài)分量的自相關(guān)函數(shù)方差小于閾值要求時則將相關(guān)IMF視為噪聲，當(dāng)自相關(guān)函數(shù)方差大于方差閾值時，則判定為有效IMF。

4）各階波浪力提取

采用FFT-FS頻譜細(xì)化技術(shù)[10]準(zhǔn)確計算各IMF的特征頻率，根據(jù)各特征頻率與遭遇頻率的倍數(shù)關(guān)系，依次確定高階倍頻波浪力模態(tài)及1階波浪力模態(tài)。將1階波浪力模態(tài)之后的各IMF及殘余分量求和，所得結(jié)果為2階低頻波浪力與定常項的混合模態(tài)，采用FFT-FS方法求混合模態(tài)在頻率為0處的幅值，該幅值即為定常力，將混合模態(tài)減去定常力幅值則可得到2階低頻波浪力組分。除縱蕩波浪力外，其余5個方向的2階定常波浪力就是這個定常力，對于縱蕩方向而言，提取2階定常波浪力時要扣除靜水阻力。

3 船模非線性波浪力測量信號處理實例

3.1 試驗概況

在武漢理工大學(xué)長條形拖曳水池開展了規(guī)則波中一約束模（記作船模a）及一半約束模（記作船模b）近距離并行頂浪航行時約束模所受波浪力的測量試驗。試驗主要目的是研究波浪中近距離并行航行兩船模之間的水動力干擾規(guī)律。船模a重173 kg，水線長3.651 m，水線寬0.462 m，設(shè)計吃水0.156 m；船b重226 kg，水線長4.352 m，水線寬0.524 m，設(shè)計吃水0.148 m，模型其他主要參數(shù)及三維輪廓參見文獻(xiàn)[12]。兩船模中對中布置，橫向間距（兩船模重心在橫向的距離）為0.35（為船模a的水線長），船模a位于船模b的√左舷，兩船模航速對應(yīng)的傅汝德數(shù)（，為重力加速度）為0.15。規(guī)則入射波波長范圍為2～10 m，波高為70 mm。試驗中船模b的縱蕩、橫蕩及首搖3個自由度進(jìn)行約束，而垂蕩、橫搖及縱搖3個自由度均自由。船a各自由度均約束，其受到的波浪力一部分為入射波及入射波遭遇該船時反射的繞射波對其的激勵作用，另一部分為船b在搖蕩過程中產(chǎn)生的輻射波及入射波經(jīng)該船體反射的繞射波傳播至船a并對其產(chǎn)生的激勵作用，這種力反映了船模b對船模a的水動力干擾。采用兩六分力天平對約束模的波浪力進(jìn)行測量，測力裝置原理如圖2所示。兩六分力天平通過法蘭分別固連于船底板的水平基座上并使受力面平行于該基座的上表平面，天平的中軸線均和船體中縱剖面重合，天平的上表平面通過法蘭連接在兩豎直的導(dǎo)航桿上。通過調(diào)節(jié)導(dǎo)航桿與拖車之間的高度使船模吃水為設(shè)計吃水，拖車啟動前，對六分力天平進(jìn)行清零處理，使得靜水浮力、船體自重及安裝過程中船模施加于天平的作用力平衡。當(dāng)拖車在規(guī)則波中勻速運動時，由于船模a各自由度運動均約束，因此船模a處于平衡狀態(tài)。將兩天平測的波浪力及力矩簡化到船模的重心處并進(jìn)行合成，合成結(jié)果則為船模a所受波浪作用力。

3.2 實測數(shù)據(jù)分析

3.2.1 CEEMD和EMD處理結(jié)果的比較

圖3及圖4分別給出了當(dāng)入射波波長λ為4.5 m時船a橫搖力矩測量信號的EMD及CEEMD分解結(jié)果。在分解前，對測試信號進(jìn)行低通濾波處理，濾波后的橫搖力矩時歷曲線如圖3第1個子圖所示。從圖3可知，EMD分解得到6個IMF和1個殘差分量。前3個模態(tài)存在明顯的模態(tài)混疊現(xiàn)象：比如在采樣時間t位于0～2.5 s，10～11 s，14～16 s的區(qū)間內(nèi)時，IMF2模態(tài)中含有IMF1的分量，又比如在區(qū)間0～5 s及14～16 s內(nèi)，IMF3中含有IMF2的分量。這種模態(tài)混疊現(xiàn)象使得原模態(tài)在局部發(fā)生突變，對應(yīng)的信號也變得不再平穩(wěn)，違背了EMD算法能夠改善不同時間尺度IMF平穩(wěn)性的初衷。采用FFT-FS頻譜細(xì)化算法求得IMF1～I(xiàn)MF3的幅值譜如圖5所示，該圖能夠更直觀地觀察到信號在頻率尺度的混疊現(xiàn)象：比如IMF1峰值頻率為2=9.92 rad/s，而在該頻率附近5～15 rad/s的區(qū)間范圍內(nèi)，也有一定比重的信號分量存在；又比如IMF2和IMF3的峰值頻率均為=4.96 rad/s，且在該頻率附近頻率區(qū)間內(nèi)，也存在一定的信號分量。綜上可見，這3個IMF不是期望獲得的船模非線性波浪力重要組分的提取結(jié)果。究其原因主要是由于原測量信號中存在著隨機噪聲（這一點將在下文對CEEMD結(jié)果分析時進(jìn)一步說明），導(dǎo)致被分解信號的極值點移位，無法正確擬合出信號中的上下包絡(luò)線，求出的上下包絡(luò)線的均值存在嚴(yán)重失真，從而使頻率尺度產(chǎn)生混疊。因此對原測量信號直接采用EMD方法是無法正確提取船模非線性波浪力中各階組分的。

表 1 各歸一化自相關(guān)函數(shù)的方差Tab. 1 Normalized autocorrelation function variance of the IMFs and the residual term

從圖4可知采用CEEMD分解得到9個IMF分量及1個殘余項，對這10個分量求取歸一化自相關(guān)函數(shù)并計算方差，結(jié)果如圖6及表1所示。從圖6可知，前4個模態(tài)的自相關(guān)函數(shù)在零點處取得最大值，而在其余點則迅速衰減，并接近于0，非常接近隨機噪聲自相關(guān)函數(shù)的特性。對于其余6個分量，其自相關(guān)函數(shù)在零點處也取得最大值，但其余點隨著時間差的增大緩慢變化下降至0，符合不含噪信號的特性。再從各分量歸一化自相關(guān)函數(shù)的方差來看，按照文獻(xiàn)[5]推薦的方差閾值0.03，前4個分量的方差值均小于0.03，判定這4個模態(tài)為隨機噪聲占主導(dǎo)的模態(tài)；其余6個分量方差值均大于0.03，可以判定為船模波浪力的有效模態(tài)。值得一提的是，在CEEMD分解過程中雖然人為地加入了高斯白噪聲，但這種高斯白噪聲和這4個噪聲模態(tài)有本質(zhì)區(qū)別，這4個模態(tài)中的隨機噪聲真實地存在于測量信號中，而人為加入的高斯白噪聲由于正負(fù)成對，在合成過程中會相互抵消，不會殘留在分解結(jié)果中。為進(jìn)一步說明這個問題，將CEEMD所得的各分量進(jìn)行重構(gòu)，并與原分解信號取差值，可以得到重構(gòu)誤差如圖7所示。從圖可見重構(gòu)誤差可達(dá)到10–15的量級，遠(yuǎn)小于隨機噪聲占主導(dǎo)的4個模態(tài)的量級，因此分解得到的前4個IMF確實為信號中的真實分量。正是由于這些隨機噪聲的存在導(dǎo)致EMD運算所得的IMF存在模態(tài)混疊現(xiàn)象。而CEEMD算法由于輔助加入頻率平均的高斯白噪聲，使得含隨機噪聲的船模非線性波浪力信號變得連續(xù)，在求取IMF的過程中可以抑制隨機噪聲引起的極值點移位及包絡(luò)線擬合偏差等問題，從而抑制模態(tài)混疊現(xiàn)象。圖8給出了圖3中IMF5，IMF6幅值譜，從圖可見這2個模態(tài)沒有出現(xiàn)模態(tài)混疊現(xiàn)象，根據(jù)其峰值頻率可以判定IMF5對應(yīng)于橫搖力矩測量數(shù)據(jù)中的2階倍頻波浪力，IMF6對應(yīng)于1階波浪力。按第2節(jié)所提方法對后續(xù)模態(tài)及殘差進(jìn)行處理可以得到2階定常波浪力及2階低頻波浪力組分。

3.2.2 共振工況下船模非線性波浪力測量信號分析

波浪中兩船近距離并行航行時，兩船內(nèi)側(cè)在船長方向形成一個前寬后窄的狹長型間隙。入射波、兩船搖蕩產(chǎn)生的輻射波，入射波遭遇船體以后產(chǎn)生的繞射波會在這一狹窄間隙內(nèi)相互干涉、衍射，形成復(fù)雜的波系干擾。由于兩船內(nèi)側(cè)的狹長型間隙入口寬，出口窄，對船間流體產(chǎn)生阻塞作用，使得一部分波浪難以從間歇出口溢出，另一部分波浪則通過船舷后發(fā)生反射，逐步從間隙出口傳播出去[13]。和單船在波浪中航行相比，這種阻塞作用和波浪的反射作用勢必引起船體波浪力的變化。特別是在一些特定頻率處，比如船體的共振頻率或流體自身的共振頻率，間隙內(nèi)的波面抬升顯著，作用于船體的波浪力非線性加劇，進(jìn)而加劇船體的搖動運動，給船舶安全航行帶來影響。因此，研究這些共振頻率下船體所受波浪力的非線性特征本身就具有重要的理論和應(yīng)用價值。

以本文試驗為例，當(dāng)入射波波長為4.5 m時，該波長與船模b船長基本相等，遭遇頻率=4.96 rad/s，船模b在此頻率處發(fā)生共振。采用本文信號處理方法對船模a在該工況下測得的縱蕩方向的波浪力進(jìn)行分析，為便于比較，對船模a單獨航行時的測量信號也進(jìn)行分析。當(dāng)有船模b的干擾存在時，從船模a所受的縱蕩方向波浪力測量信號中提取到了4個波浪力組分，包括3階倍頻縱蕩波浪力、1階縱蕩波浪力、2階低頻縱蕩波浪力及阻力（靜水阻力和縱向2階定常波浪力之和）。無干擾時，得到3個波浪力組分，包括1階縱蕩波浪力、2階低頻縱蕩波浪力及阻力。結(jié)果如圖9～圖12所示。從圖可知：1）有干擾時1階縱蕩波浪力的平均幅值約為無干擾時的1.5倍（見圖10），若船a縱蕩方向自由，勢必引起一階縱蕩運動響應(yīng)幅值的顯著增加；2）阻力增值顯著，有干擾及無干擾時阻力值分別為6.29 N及3.96 N（見圖12頻率為0時的幅值），實測的靜水阻力值為3.26 N，無干擾時阻力約為靜水阻力的1.2倍，有干擾時的阻力約為靜水阻力的1.9倍，這種顯著的阻力增值正是由于船間狹窄間隙對波浪產(chǎn)生阻塞，在船體縱向形成巨大的壓力差而形成的；3）有無干擾時2階低頻縱蕩波浪力差異不大，且均未出現(xiàn)2階倍頻縱蕩波浪力組分，但有干擾時，測量信號中存在一個3階倍頻波浪力組分，這說明在該工況下，2階倍頻波浪力量值要遠(yuǎn)小于3倍頻波浪力模態(tài)，這和常規(guī)情況下2階倍頻波浪力在高階波浪力中占主導(dǎo)的情況是不同的。

3.2.3 船模2階波浪漂移力分析實例

使船體在波浪中產(chǎn)生漂移的波浪作用力稱作波浪漂移力，是2階波浪力在水平方向的分量，其成分包括定常波浪漂移力、低頻波浪漂移力及2倍頻波浪漂移力。波浪中近距離并行航行的兩船，在定常波浪漂移力及低頻波浪漂移力的作用下，會產(chǎn)生平均位置上的偏移，這種偏移運動可能引起兩船在水平面內(nèi)的相撞。對于波浪中兩船橫向補給作業(yè)而言，這種波浪漂移力是自動舵設(shè)計及兩船航向穩(wěn)定性設(shè)計的重要輸入。因此兩船近距離并行航行時所受波浪漂移力的研究具有重要的工程實際意義。仍以本文的試驗為例，利用本文所提的信號分析方法分別對不同波長的入射波激勵時船模a所受橫蕩波浪力及首搖力矩的測量信號進(jìn)行分析，提取了各工況下橫向定常漂移力及定常首搖力矩，結(jié)果如圖13和圖14所示。

圖13給出了船模a受到的無因次橫向定常波浪漂移力隨波長船長比的變化曲線。從圖可知，當(dāng)＜1.5時，船模a受到的橫向定常波浪漂移力為正值，即指向該船的左舷，表明此時船模b的干擾作用對船模a有排斥作用；當(dāng)＞1.5時，船模a受到的橫向定常波浪漂移力為負(fù)值，表明此時船模b對船模a的作用力為吸力。

圖14給出了船模a受到的無因次定常搖首力矩隨波長船長比的變化曲線，其變化規(guī)律和橫向定常波浪漂移力變化規(guī)律基本一致，當(dāng)＜1.5時，首搖力矩方向在水平面內(nèi)為逆時針方向（繞豎直向上的鉛垂軸），即船模a與船舶b首部相排斥，而尾部相吸引；當(dāng)＜1.5時，首搖力矩方向在水平面內(nèi)為順時針方向，即兩船首部相吸，尾部相斥，但這種作用比短波中要小得多。

4 結(jié) 語

1）直接采用EMD算法對含有隨機噪聲的船模非線性波浪力測量信號進(jìn)行分析，會出現(xiàn)模態(tài)混疊現(xiàn)象，無法正確提取非線性波浪力中的各主要組分。

2）CEEMD算法能夠過濾船模非線性波浪力測量信號中的隨機噪聲，并在分解過程中抑制EMD算法出現(xiàn)的模態(tài)混疊現(xiàn)象。

3）當(dāng)船模非線性波浪力測量信號中高階倍頻波浪力的量值遠(yuǎn)小于信號中主成分量值時，CEEMD算法會將這些波浪力成分歸于隨機噪聲模態(tài)中，這有利于對船模非線性波浪力主要組分的提取，消除高階小量對主要組分提取精度的影響。

4）利用IMF歸一化自相關(guān)函數(shù)的方差值能夠有效識別噪聲模態(tài)和有效模態(tài)。

5）本文所建立信號處理方法能夠準(zhǔn)確提取船模非線波浪力測量信號中的重要波浪力組分，為船模在規(guī)則波中的阻力增值測量，船模1階、2階波浪力測量等試驗提供了一種可靠的數(shù)據(jù)分析方法。