基于CEEMD的船模非线性波浪力测量信号分析

2018-10-15王启兴欧勇鹏

舰船科学技术 2018年9期

王启兴，许勇，欧勇鹏

(1. 海军驻柳州434厂军事代表室, 广西柳州 545000；2. 海军工程大学舰船学院, 湖北武汉 430033)

0 引言

船模试验在船舶流体力学研究领域扮演着重要的角色，是考核理论计算方法和阐述船舶流体力学机理的关键。在船模的耐波性试验中常需要测量船模在规则波中所受的波浪作用力。这种波浪力可分为两大部分：一部分为具有零平均性质的线性振荡波浪力，其振荡频率为船模的遭遇频率其中，及分别为入射波频率、波数及浪向角，为船模航速），称为一阶波浪力；另一部分为非线性影响产生的非线性波浪力，称为高阶波浪力[1]。通常二阶波浪力在高阶波浪力中占主导，但当船模遭遇频率和其运动模态的共振频率相近时，船模所受的波浪力非线性加剧，更高阶的波浪力组分也会在高阶波浪力中占有一定的比重。完整的二阶波浪力中包括定常二阶波浪力、低频二阶波浪力和倍频二阶波浪力[2]。综上可见，船模在波浪中所受的波浪力组分非常丰富，且是非线性的，因此模型试验中测得的波浪力信号也是非线性的。实测船模波浪力信号中往往还包含着高频噪声，其主要来源于测量设备中的交流电噪声和拖车及其旋转机械的振动噪声[3]。测量仪器过热时的零点漂移、仪器调零不精确等因素还可能使测量信号中含有低频趋势项的影响，使得测量波形发生偏移。另外，测量信号中往往还存在着随机噪声，这种噪声由时间上随机产生的大量起伏骚扰积累而造成的，比如试验过程中等波时间不够，水面上的随机无序的破碎波对船模的累积作用就会产生随机噪声。高频噪声、低频趋势项、随机噪声和非线性的波浪力叠加在一起，使得实测信号不再是平稳信号。图1给出了某船模在规则波中顶浪航行时首摇力矩测量信号时历曲线及对应的幅值谱，可见该信号是非平稳的。在频率＜50 rad/s的频段内，可以清晰地观察到信号中含有一阶波浪力（与对应）、定常2阶波浪力（与=0 rad/s对应）、低频2阶波浪力（介于0～之间）、2阶倍频波浪力（与2对应）及更高倍频波浪力组分，而在=315 rad/s（约为50 Hz）附近也有一定的信号组分存在，对应于测量信号中的高频噪声。如何对这种含有高频噪声及随机噪声的非平稳信号进行分析，准确提取波浪力的各阶组分是本文的研究目的。

经验模态分解算法（Empirical Mode Decomposition，EMD）[4]集平稳化、层次处理能力于一体，为提取船模非线性波浪力各组分提供了新思路。该方法基于信号的局部特征时间尺度，将非平稳信号分解，得到一系列不同频率的分量——本征模态函数（IMF）。相对于原非平稳、非线性信号，不同时间尺度的IMF平稳性得到很大改善。分解出最初IMF序列代表原始序列的高频组分，随着分解的进行，相应IMF的频率变小、周期增长。这些IMF可作为原序列的一组完全或几乎正交的展开基，并与最终的分解残差能够很好的重构原时间序列[5]。船模在规则波中所受波浪力的各阶组分对应着不同频率的分量，具备良好的局部尺度特征，从理论上而言，应用EMD方法进行处理，各阶波浪力组分应该能够很好的与分解得到的IMF相对应。但由于实测数据因随机因素或测量仪器等影响会带有随机误差，使得很多情况下EMD直接对信号分解提取的规律并不准确，容易出现模态混叠现象，即相同频率的信号分量被分解到不同的模态分量当中，或者不同频率的信号分量被分解到同一个模态分量当中[5]。集合经验模态分解（Ensemble Empirical Mode Decomposition，EEMD）[6]对EMD方法进行改进，通过加入不同随机白噪声来平衡测量数据中随机误差的影响，进而从数据中提取更准确的规律，但也因为如此，增大了EEMD分解后信号的重构误差。Yeh等对EEMD算法进行了改进，提出互补集合经验模态分解方法[7]，与EEMD算法不同的是所加入的白噪声正负成对，能够极大地提高重构精度，同时减少迭代次数并提高运算效率。本文基于CEEMD算法，并辅以合适的其他信号分析手段，研究船模非线性波浪力测量信号各阶组分的提取方法。

1 EMD和CEEMD算法原理

1.1 EMD算法原理

EMD方法的实质是提取本征模态函数IMF的过程。IMF分量必须满足以下2个条件：1）其极值个数和过零点数相同或最多相差一个；2）其上下包络线关于时间轴局部对称。分解方法如下[8]：

2）对信号中计算得到的极大（小）极值点利用3次样条法进行插值，并拟合成信号的上下包络线，计算2条曲线的均值序列，求出。

1.2 CEEMD算法原理

CEEMD算法实质上是EMD算法的一种辅助算法，其核心仍是EMD，但在进行EMD前在原始信号中添加了正负对形式的白噪声，其过程如下[9]：

2）对各组信号采用EMD方法进行分解，集合中的信号都可以分解为IMF分量，其中第 i个信号的第个IMF分量为。

3）通过多组分量组合得到分解结果如式（3）所示，残余分量记作。

2 船模非线性波浪力测量信号的处理方法

根据船模非线性波浪力测量信号的特点，本文提出如下信号处理流程：

1）低通滤波

相对于高频噪声的频率而言，船模在规则波中的遭遇频率是低频的[10]，其频率往往低于5 Hz。大量的试验数据分析表明，5阶以上的高阶波浪力量值相对于1阶、甚至2阶倍频波浪力量值而言，可以忽略，因此有实际意义的高阶波浪力的频率相对于高频噪声的频率而言也是低频的。本文采用低通滤波方法去除实测波浪力信号中的高频噪声。低通滤波的具体实施方法见文献[3]，高频截断频率取6。采用低通滤波的作用过滤高频噪声，减少后续CEEMD分解时的迭代次数，提高运算效率。

2）CEEMD求解

对过滤高频噪声的船模波浪力测量信号进行CEEMD分解，获得一系列由高频到低频排列的本征模态函数IMF。分解所得IMF中的前几个模态往往对应于随机噪声占主导的模态。另外，当信号中多倍频的高阶波浪力为小量，且量值与这些随机信号相当时，也会与随机噪声耦合在一起，存在于这几个模态当中。因此，必须要对这些模态进行甄别，为提取与船模非线性波浪力各组分相对应的有效IMF提供方便。

3）有效IMF识别

随机信号的自相关函数反映了信号与其自身在不同时间点的相似程度。对于随机噪声而言，其自相关函数在零点处最大，其余点处立即衰减为0，而不含噪信号自相关函数值在零点处最大，其余点随着时间差的增大而缓慢衰减至0[11]。本文利用随机噪声和不含噪信号自相关函数的性质差异来对噪声模态和有效模态进行识别。步骤如下：分别求取各IMF的归一化自相关函数及其方差，设定自相关函数方差的阈值，当模态分量的自相关函数方差小于阈值要求时则将相关IMF视为噪声，当自相关函数方差大于方差阈值时，则判定为有效IMF。

4）各阶波浪力提取

采用FFT-FS频谱细化技术[10]准确计算各IMF的特征频率，根据各特征频率与遭遇频率的倍数关系，依次确定高阶倍频波浪力模态及1阶波浪力模态。将1阶波浪力模态之后的各IMF及残余分量求和，所得结果为2阶低频波浪力与定常项的混合模态，采用FFT-FS方法求混合模态在频率为0处的幅值，该幅值即为定常力，将混合模态减去定常力幅值则可得到2阶低频波浪力组分。除纵荡波浪力外，其余5个方向的2阶定常波浪力就是这个定常力，对于纵荡方向而言，提取2阶定常波浪力时要扣除静水阻力。

3 船模非线性波浪力测量信号处理实例

3.1 试验概况

在武汉理工大学长条形拖曳水池开展了规则波中一约束模（记作船模a）及一半约束模（记作船模b）近距离并行顶浪航行时约束模所受波浪力的测量试验。试验主要目的是研究波浪中近距离并行航行两船模之间的水动力干扰规律。船模a重173 kg，水线长3.651 m，水线宽0.462 m，设计吃水0.156 m；船b重226 kg，水线长4.352 m，水线宽0.524 m，设计吃水0.148 m，模型其他主要参数及三维轮廓参见文献[12]。两船模中对中布置，横向间距（两船模重心在横向的距离）为0.35（为船模a的水线长），船模a位于船模b的√左舷，两船模航速对应的傅汝德数（，为重力加速度）为0.15。规则入射波波长范围为2～10 m，波高为70 mm。试验中船模b的纵荡、横荡及首摇3个自由度进行约束，而垂荡、横摇及纵摇3个自由度均自由。船a各自由度均约束，其受到的波浪力一部分为入射波及入射波遭遇该船时反射的绕射波对其的激励作用，另一部分为船b在摇荡过程中产生的辐射波及入射波经该船体反射的绕射波传播至船a并对其产生的激励作用，这种力反映了船模b对船模a的水动力干扰。采用两六分力天平对约束模的波浪力进行测量，测力装置原理如图2所示。两六分力天平通过法兰分别固连于船底板的水平基座上并使受力面平行于该基座的上表平面，天平的中轴线均和船体中纵剖面重合，天平的上表平面通过法兰连接在两竖直的导航杆上。通过调节导航杆与拖车之间的高度使船模吃水为设计吃水，拖车启动前，对六分力天平进行清零处理，使得静水浮力、船体自重及安装过程中船模施加于天平的作用力平衡。当拖车在规则波中匀速运动时，由于船模a各自由度运动均约束，因此船模a处于平衡状态。将两天平测的波浪力及力矩简化到船模的重心处并进行合成，合成结果则为船模a所受波浪作用力。

3.2 实测数据分析

3.2.1 CEEMD和EMD处理结果的比较

图3及图4分别给出了当入射波波长λ为4.5 m时船a横摇力矩测量信号的EMD及CEEMD分解结果。在分解前，对测试信号进行低通滤波处理，滤波后的横摇力矩时历曲线如图3第1个子图所示。从图3可知，EMD分解得到6个IMF和1个残差分量。前3个模态存在明显的模态混叠现象：比如在采样时间t位于0～2.5 s，10～11 s，14～16 s的区间内时，IMF2模态中含有IMF1的分量，又比如在区间0～5 s及14～16 s内，IMF3中含有IMF2的分量。这种模态混叠现象使得原模态在局部发生突变，对应的信号也变得不再平稳，违背了EMD算法能够改善不同时间尺度IMF平稳性的初衷。采用FFT-FS频谱细化算法求得IMF1～IMF3的幅值谱如图5所示，该图能够更直观地观察到信号在频率尺度的混叠现象：比如IMF1峰值频率为2=9.92 rad/s，而在该频率附近5～15 rad/s的区间范围内，也有一定比重的信号分量存在；又比如IMF2和IMF3的峰值频率均为=4.96 rad/s，且在该频率附近频率区间内，也存在一定的信号分量。综上可见，这3个IMF不是期望获得的船模非线性波浪力重要组分的提取结果。究其原因主要是由于原测量信号中存在着随机噪声（这一点将在下文对CEEMD结果分析时进一步说明），导致被分解信号的极值点移位，无法正确拟合出信号中的上下包络线，求出的上下包络线的均值存在严重失真，从而使频率尺度产生混叠。因此对原测量信号直接采用EMD方法是无法正确提取船模非线性波浪力中各阶组分的。

表 1 各归一化自相关函数的方差Tab. 1 Normalized autocorrelation function variance of the IMFs and the residual term

从图4可知采用CEEMD分解得到9个IMF分量及1个残余项，对这10个分量求取归一化自相关函数并计算方差，结果如图6及表1所示。从图6可知，前4个模态的自相关函数在零点处取得最大值，而在其余点则迅速衰减，并接近于0，非常接近随机噪声自相关函数的特性。对于其余6个分量，其自相关函数在零点处也取得最大值，但其余点随着时间差的增大缓慢变化下降至0，符合不含噪信号的特性。再从各分量归一化自相关函数的方差来看，按照文献[5]推荐的方差阈值0.03，前4个分量的方差值均小于0.03，判定这4个模态为随机噪声占主导的模态；其余6个分量方差值均大于0.03，可以判定为船模波浪力的有效模态。值得一提的是，在CEEMD分解过程中虽然人为地加入了高斯白噪声，但这种高斯白噪声和这4个噪声模态有本质区别，这4个模态中的随机噪声真实地存在于测量信号中，而人为加入的高斯白噪声由于正负成对，在合成过程中会相互抵消，不会残留在分解结果中。为进一步说明这个问题，将CEEMD所得的各分量进行重构，并与原分解信号取差值，可以得到重构误差如图7所示。从图可见重构误差可达到10–15的量级，远小于随机噪声占主导的4个模态的量级，因此分解得到的前4个IMF确实为信号中的真实分量。正是由于这些随机噪声的存在导致EMD运算所得的IMF存在模态混叠现象。而CEEMD算法由于辅助加入频率平均的高斯白噪声，使得含随机噪声的船模非线性波浪力信号变得连续，在求取IMF的过程中可以抑制随机噪声引起的极值点移位及包络线拟合偏差等问题，从而抑制模态混叠现象。图8给出了图3中IMF5，IMF6幅值谱，从图可见这2个模态没有出现模态混叠现象，根据其峰值频率可以判定IMF5对应于横摇力矩测量数据中的2阶倍频波浪力，IMF6对应于1阶波浪力。按第2节所提方法对后续模态及残差进行处理可以得到2阶定常波浪力及2阶低频波浪力组分。

3.2.2 共振工况下船模非线性波浪力测量信号分析

波浪中两船近距离并行航行时，两船内侧在船长方向形成一个前宽后窄的狭长型间隙。入射波、两船摇荡产生的辐射波，入射波遭遇船体以后产生的绕射波会在这一狭窄间隙内相互干涉、衍射，形成复杂的波系干扰。由于两船内侧的狭长型间隙入口宽，出口窄，对船间流体产生阻塞作用，使得一部分波浪难以从间歇出口溢出，另一部分波浪则通过船舷后发生反射，逐步从间隙出口传播出去[13]。和单船在波浪中航行相比，这种阻塞作用和波浪的反射作用势必引起船体波浪力的变化。特别是在一些特定频率处，比如船体的共振频率或流体自身的共振频率，间隙内的波面抬升显著，作用于船体的波浪力非线性加剧，进而加剧船体的摇动运动，给船舶安全航行带来影响。因此，研究这些共振频率下船体所受波浪力的非线性特征本身就具有重要的理论和应用价值。

以本文试验为例，当入射波波长为4.5 m时，该波长与船模b船长基本相等，遭遇频率=4.96 rad/s，船模b在此频率处发生共振。采用本文信号处理方法对船模a在该工况下测得的纵荡方向的波浪力进行分析，为便于比较，对船模a单独航行时的测量信号也进行分析。当有船模b的干扰存在时，从船模a所受的纵荡方向波浪力测量信号中提取到了4个波浪力组分，包括3阶倍频纵荡波浪力、1阶纵荡波浪力、2阶低频纵荡波浪力及阻力（静水阻力和纵向2阶定常波浪力之和）。无干扰时，得到3个波浪力组分，包括1阶纵荡波浪力、2阶低频纵荡波浪力及阻力。结果如图9～图12所示。从图可知：1）有干扰时1阶纵荡波浪力的平均幅值约为无干扰时的1.5倍（见图10），若船a纵荡方向自由，势必引起一阶纵荡运动响应幅值的显著增加；2）阻力增值显著，有干扰及无干扰时阻力值分别为6.29 N及3.96 N（见图12频率为0时的幅值），实测的静水阻力值为3.26 N，无干扰时阻力约为静水阻力的1.2倍，有干扰时的阻力约为静水阻力的1.9倍，这种显著的阻力增值正是由于船间狭窄间隙对波浪产生阻塞，在船体纵向形成巨大的压力差而形成的；3）有无干扰时2阶低频纵荡波浪力差异不大，且均未出现2阶倍频纵荡波浪力组分，但有干扰时，测量信号中存在一个3阶倍频波浪力组分，这说明在该工况下，2阶倍频波浪力量值要远小于3倍频波浪力模态，这和常规情况下2阶倍频波浪力在高阶波浪力中占主导的情况是不同的。

3.2.3 船模2阶波浪漂移力分析实例

使船体在波浪中产生漂移的波浪作用力称作波浪漂移力，是2阶波浪力在水平方向的分量，其成分包括定常波浪漂移力、低频波浪漂移力及2倍频波浪漂移力。波浪中近距离并行航行的两船，在定常波浪漂移力及低频波浪漂移力的作用下，会产生平均位置上的偏移，这种偏移运动可能引起两船在水平面内的相撞。对于波浪中两船横向补给作业而言，这种波浪漂移力是自动舵设计及两船航向稳定性设计的重要输入。因此两船近距离并行航行时所受波浪漂移力的研究具有重要的工程实际意义。仍以本文的试验为例，利用本文所提的信号分析方法分别对不同波长的入射波激励时船模a所受横荡波浪力及首摇力矩的测量信号进行分析，提取了各工况下横向定常漂移力及定常首摇力矩，结果如图13和图14所示。

图13给出了船模a受到的无因次横向定常波浪漂移力随波长船长比的变化曲线。从图可知，当＜1.5时，船模a受到的横向定常波浪漂移力为正值，即指向该船的左舷，表明此时船模b的干扰作用对船模a有排斥作用；当＞1.5时，船模a受到的横向定常波浪漂移力为负值，表明此时船模b对船模a的作用力为吸力。

图14给出了船模a受到的无因次定常摇首力矩随波长船长比的变化曲线，其变化规律和横向定常波浪漂移力变化规律基本一致，当＜1.5时，首摇力矩方向在水平面内为逆时针方向（绕竖直向上的铅垂轴），即船模a与船舶b首部相排斥，而尾部相吸引；当＜1.5时，首摇力矩方向在水平面内为顺时针方向，即两船首部相吸，尾部相斥，但这种作用比短波中要小得多。