劉志慧， 徐興平， 牛懷磊， 王 騰

(1.中國石油大學(xué)(華東) 石油工程學(xué)院，山東 青島 266580； 2.中國石油大學(xué)(華東) 機(jī)電工程學(xué)院，山東 青島 266580；3.中國石油集團(tuán)海洋工程(青島)有限公司，山東 青島 266555)

隨著深海資源的開發(fā)，海洋立管作為海洋平臺中生產(chǎn)運(yùn)輸?shù)年P(guān)鍵設(shè)備，受到眾多研究者的關(guān)注。準(zhǔn)確分析海洋立管渦激振動的特性，是設(shè)計(jì)海洋立管的依據(jù)，也是損傷分析的基礎(chǔ)。實(shí)際監(jiān)測中，真實(shí)信號受到環(huán)境的影響，產(chǎn)生較大的干擾易導(dǎo)致信號失真，因此需對測量數(shù)據(jù)降噪處理，進(jìn)行信號識別。對立管渦激振動非線性信號的辨識，主要集中在利用傅里葉變換識別渦激振動頻率，多數(shù)的研究手段是基于線性和平穩(wěn)的假設(shè)進(jìn)行的，因此有必要引入新的方法進(jìn)行研究。Huang等[1]提出一種非參數(shù)的數(shù)據(jù)驅(qū)動算法，即經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，EMD)，用來解決現(xiàn)實(shí)世界的非線性非平穩(wěn)信號的分解和時(shí)頻分析。EMD自提出后，在很多領(lǐng)域都得到了廣泛的應(yīng)用[2-5]。廣大學(xué)者在研究中發(fā)現(xiàn)EMD方法存在一個(gè)最主要問題模態(tài)混淆現(xiàn)象：在同一個(gè)本征模態(tài)函數(shù)(intrinsic mode function，IMF)分量中存在不同幅值或頻率的信號，或者在不同的IMF分量中存在非常相似的信號。Wu等[6]在研究白噪聲信號的統(tǒng)計(jì)特征基礎(chǔ)上，提出總體平均經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition,EEMD)方法，依據(jù)數(shù)據(jù)本身的時(shí)間尺度特征進(jìn)行時(shí)頻分析，分解的過程保留了數(shù)據(jù)本身的特性，在處理信號的過程中能夠避免傅里葉變換帶來的缺陷，如出現(xiàn)虛假頻率和多余信號分量等[7-9]。EEMD成功解決了模態(tài)混合問題，Yi等[12-13]提出置信指數(shù)算法和最大能量法選擇敏感的IMF，用于研究軸承故障診斷。Prosvirin等[14]在研究滾動軸承故障診斷時(shí)將碰摩的存在度比率和基于Kullback-Leibler散度的相似性度量結(jié)合到IMF質(zhì)量度量中來選擇有意義的信號主模式。

EEMD雖然有效解決了模態(tài)混合問題，但在完成分解過程后生成有限的IMF集合，它的數(shù)量遠(yuǎn)遠(yuǎn)大于有實(shí)際物理意義的IMF數(shù)量。如何自適應(yīng)地確定海洋立管渦激振動包含特定振動信息的具有實(shí)際物理意義的IMF成為面臨的困難。因此，本文基于EEMD對海洋立管渦激振動特征進(jìn)行識別，基于響應(yīng)數(shù)據(jù)分解，考慮光滑度、相似度、存在度建立最優(yōu)降噪光滑數(shù)學(xué)模型，以獲得具有實(shí)際物理意義的有效IMF分量。

1 基于EEMD的最優(yōu)降噪光滑數(shù)學(xué)模型

為獲得具有實(shí)際物理意義的有效IMF分量，需要首先對立管渦激振動響應(yīng)進(jìn)行降噪處理，對經(jīng)過EEMD分解的IMF進(jìn)行信號重構(gòu)，在此基礎(chǔ)上考慮信號的相似度、光滑度和存在度，建立最優(yōu)降噪光滑算法的數(shù)學(xué)模型，基于此進(jìn)行有效IMF的選擇。該模型由總體平均經(jīng)驗(yàn)?zāi)B(tài)分解、重構(gòu)信號、目標(biāo)尋優(yōu)三部分組成。

1.1 總體平均經(jīng)驗(yàn)?zāi)B(tài)分解

總體平均經(jīng)驗(yàn)?zāi)B(tài)分解利用白噪聲頻譜均勻分布的特點(diǎn)，對原始信號多次加入不同的高斯白噪聲進(jìn)行EMD分解。增加的白噪聲與不同尺度的組成成分均勻地填充在整個(gè)時(shí)頻空間。當(dāng)信號被添加在這個(gè)均勻分布的背景中時(shí)，不同尺度的信號將自動投影到由背景中白噪聲建立的參考尺度上。增加的白噪聲用足夠多的試驗(yàn)進(jìn)行平均抵消，只有信號分量在總體均值中保留下來。EEMD具體步驟如下。

步驟1原始信S(t)中多次添加自適應(yīng)性白噪聲信號Bi(t)，得到加噪信號xi(t)

xi(t)=S(t)+σiBi(t)

(1)

式中:i為添加白噪聲的次數(shù)；σi為第i次添加的白噪聲的標(biāo)準(zhǔn)差。

步驟2EMD分解，將多次分解的結(jié)果取平均得到IMFs

(2)

S(t)經(jīng)過EEMD分解，可表示為

(3)

式中：N為添加白噪聲的總次數(shù)，本研究中取100；fim，j為第i次添加白噪聲后EMD分解得到的第j個(gè)固有模態(tài)函數(shù)分量；rn(t)為原始信號的趨勢分量。

1.2 重構(gòu)信號(濾波算法)

信號EEMD分解后，固有模態(tài)函數(shù)按頻率高低排列。高頻信號一般為噪聲引起，頻率最低的余項(xiàng)為趨勢項(xiàng)，有效信號頻率往往處于高頻和低頻之間，因此利用帶通濾波的方法，同時(shí)去掉若干個(gè)高頻和低頻IMF分量，再以其余分量重構(gòu)信號[15]，得到降噪后的信號s(t)，可表示為

(4)

1.3 目標(biāo)尋優(yōu)

考慮重構(gòu)信號與原始信號的相似度、重構(gòu)信號算法的光滑程度，構(gòu)建最優(yōu)化區(qū)間，根據(jù)存在度進(jìn)行目標(biāo)尋優(yōu)。

1.3.1 相似度

相似度指的是濾波算法得到的重構(gòu)信號和原始信號的相似程度，重構(gòu)信號曲線越接近原始信號曲線，說明相似程度越高，用均方誤差來表示

(5)

式中，N1為信號長度。

均方誤差越小，表示重構(gòu)信號曲線越接近原始信號，兩個(gè)曲線的相似程度越高,較小的均方誤差可以使重構(gòu)信號保存原始監(jiān)測信號的真實(shí)性。

1.3.2 光滑度

光滑度指的是曲線的曲折程度。曲線f(x)光滑是指其在區(qū)間[a，b]內(nèi)連續(xù)，具有一階連續(xù)導(dǎo)數(shù)，且切線隨著切點(diǎn)的移動連續(xù)轉(zhuǎn)動即曲率相等。

若曲線連續(xù)，則有

f(x0)=g(x0)

(6)

切線連續(xù)，則一階導(dǎo)數(shù)相等

f′(x0)=g′(x0)

(7)

曲率相等，則有

(8)

將式(7)代入式(8)可得

f″(x0)=g″(x0)

(9)

將二階導(dǎo)數(shù)按離散公式展開，左右兩側(cè)分別為

(10)

(11)

將式(10)和式(11)代入式(9)可得

f(x0+2h)-f(x0-2h)-2[f(x0+h)-f(x0-h)]=0(12)

因此，對于任意一條曲線，任意連接點(diǎn)x0處光滑度的指標(biāo)可表示為

MS,f|x=x0=f(x0+2h)-f(x0-2h)-

2[f(x0+h)-f(x0-h)]

(13)

式中:h為步長;MS,f越接近于零，則曲線f(x)在除左右兩個(gè)端點(diǎn)之外的x0鄰域內(nèi)越光滑。

1.3.3 存在度

基于EEMD得到的IMF建立濾波算法，考慮算法的相似度指標(biāo)和光滑度指標(biāo)，構(gòu)造目標(biāo)函數(shù)[16-17]為

{f}={λ(EMS,f)+(1-λ)(MS,f)}

(14)

式中:λ為相似度的影響因子；1-λ為光滑度影響因子。

考慮20%浮動，最優(yōu)化區(qū)間下限設(shè)置為目標(biāo)函數(shù)的最小值，上限設(shè)置為目標(biāo)函數(shù)最小值的1.2倍,即[min{f},1.2min{f}] 。FIM,i的存在度用第i個(gè)IMF在最優(yōu)化區(qū)間中所有算法中出現(xiàn)的頻率表現(xiàn)，可以表示為

(15)

式中，p為IMF出現(xiàn)的頻率。

最優(yōu)化區(qū)間中，若FIM，a在所有算法中存在度最高，則判斷FIM，a為有效IMF，對應(yīng)僅包含F(xiàn)IM，a的算法為最優(yōu)光滑降噪模型。

2 仿真信號驗(yàn)證

為證明構(gòu)建的數(shù)學(xué)模型的有效性，采用仿真信號模擬的方法進(jìn)行驗(yàn)證。

由于渦激振動信號頻率交錯，很難用單一的模型表示，考慮到渦激振動信號的特點(diǎn)(順流向渦激振動的頻率是橫流向頻率的二倍)，采用兩個(gè)周期信號與間歇信號，噪聲信號相統(tǒng)一的復(fù)合信號模型來仿真模擬。

x(t)=s1(t)+s2(t)+c(t)+n(t)

(16)

周期性信號分別取為

s1(t)=10cos(20πt)

(17)

s2(t)=8cos(10πt+π/2)

(18)

間歇信號取為間斷隨機(jī)信號；噪聲信號取強(qiáng)度為0.5高斯白噪聲，各個(gè)信號及其合成信號如圖1所示。

圖1 渦激振動合成信號Fig.1 Composite signal of VIV

2.1 數(shù)值模擬渦激振動信號有效IMF選取

對數(shù)值模擬渦激振動信號x(t)進(jìn)行總體經(jīng)驗(yàn)?zāi)B(tài)分解，得到10個(gè)模態(tài)分量FIM,1-FIM,10和趨勢項(xiàng)R10，如圖2所示。

圖2 仿真模擬信號的EEMD分解Fig.2 EEMD decomposition of simulation signals

針對10個(gè)IMF，構(gòu)建28個(gè)濾波算法BP1-28，進(jìn)行最優(yōu)降噪光滑計(jì)算，結(jié)果如表1和圖3所示。

表1 濾波算法性能數(shù)據(jù)Tab.1 Performance data of filtering algorithm

圖3 濾波算法性能指標(biāo)Fig.3 Performance index of filtering algorithm

考慮到信號處理工作主要是進(jìn)行渦激振動正余弦信號的提取，因此相似度影響因子和光滑度影響因子都取0.5，即算法相似度和光滑度的影響相當(dāng)。

由表1數(shù)據(jù)可以看出：BP14，BP15，BP19，BP20，BP23，BP24這6個(gè)濾波算法目標(biāo)函數(shù)值小于0.02，其中BP14，BP15比其余濾波算法的值大一個(gè)數(shù)量級；濾波算法BP19，BP20比BP23，BP24的計(jì)算結(jié)果略大；濾波算法BP23計(jì)算結(jié)果最小，算法包含IMF6～I(xiàn)MF8。從圖3可以看出，濾波算法BP14，BP15，BP19，BP20，BP23，BP24目標(biāo)函數(shù)比其他算法要小，在0～0.05。從圖4可見：在最優(yōu)化區(qū)間中，IMF6和IMF7存在度最高。

圖4 最優(yōu)化區(qū)間存在度Fig.4 The ratio of existence for IMFs

因此從目標(biāo)函數(shù)值最小的角度考慮，建立的最優(yōu)降噪光滑數(shù)學(xué)模型為算法BP23，包含IMF6～I(xiàn)MF8?？紤]存在度，建立的最優(yōu)降噪光滑數(shù)學(xué)模型為算法BP24，包含有效本征模態(tài)函數(shù)IMF6～I(xiàn)MF7?？梢娍紤]存在度建立的最優(yōu)數(shù)學(xué)模型不一定對應(yīng)目標(biāo)函數(shù)計(jì)算結(jié)果最小的算法，這是區(qū)別于一般最優(yōu)降噪光滑數(shù)學(xué)模型的地方。

2.2 結(jié)果驗(yàn)證

為考察判斷結(jié)果的準(zhǔn)確性，從本征模態(tài)函數(shù)與原信號中周期性信號的相關(guān)度、頻率比較，以及能量貢獻(xiàn)度3個(gè)角度進(jìn)行驗(yàn)證。

由數(shù)理統(tǒng)計(jì)可知，數(shù)據(jù)相關(guān)性可用相關(guān)系數(shù)表示

(19)

式中:cov(FMI,i,x)為IMF與周期性信號x的協(xié)方差；D(FMI,i)為IMF的方差；D(x)為周期性信號x的方差。

相關(guān)系數(shù)|ρf|越趨近于1，說明該IMF和信號線性相關(guān)程度越高，二者的圖形越相似。IMF和信號s1,s2相關(guān)度如圖5所示，相關(guān)系數(shù)見表2。

圖5 IMFs和原始信號的相關(guān)度Fig.5 Correlation between IMFs & periodic signals

表2 IMFs與信號s1與s2、仿真信號相關(guān)度Tab.2 Correlation between IMFs and signals

由圖5可以看出IMF6和周期性信號s1相關(guān)度最高，IMF7和周期性信號s2最相關(guān)。由表2可知，IMF6與s1信號相關(guān)度為0.968 5，IMF7與s2信號相關(guān)度為0.976 9，二者數(shù)值在均接近于1，因此IMF6與信號s1形狀最接近；IMF7與信號s2形狀最接近，且IMF6和IMF7與仿真信號的相關(guān)度均在0.6以上。

將IMF6與IMF7做頻譜分析，可得IMF6的頻率為10 Hz，與s1的頻率一致；IMF7的頻率為5 Hz，與s2的頻率一致，如圖6所示。這說明基于構(gòu)建的最優(yōu)降噪光滑模型得到的有效IMF頻率與原始信號頻率保持一致。

圖6 IMF6與IMF7的頻譜分析圖Fig.6 Spectrum analysis of IMF6 and IMF7

在時(shí)域內(nèi)，將本征模態(tài)函數(shù)的幅值平方對時(shí)間進(jìn)行積分，可得能量譜密度，表示每個(gè)頻率在整個(gè)時(shí)間長度內(nèi)的能量積累[18]，表達(dá)式為

(20)

根據(jù)式(20)計(jì)算IMFs能量，可以得到IMFs能量分布圖譜，如圖7所示。從圖中可以看出，IMF6和IMF7在整個(gè)信號中占主要地位，經(jīng)計(jì)算可知其能量占總能量的97.98%。

圖7 IMFs能量分布Fig.7 Energy distribution of IMFs

上述從相關(guān)性分析、頻譜分析、能量譜分析3個(gè)角度證明：利用本文構(gòu)建的算法優(yōu)選的IMF6，IMF7分別和模擬渦激振動順流向特點(diǎn)的周期信號s1(t)和橫流向特點(diǎn)的周期性信號s2(t)相關(guān)性最大，且頻率對應(yīng)相等，在整個(gè)信號中能量占比高達(dá)97.98%。這說明重構(gòu)信號去除了仿真信號中的干擾信號，能夠準(zhǔn)確識別原始信號中的周期信號，保存了原始信號的真實(shí)性，驗(yàn)證了由構(gòu)建的最優(yōu)降噪光滑算法得到的有效本征模態(tài)函數(shù)的準(zhǔn)確性。

3 試驗(yàn)數(shù)據(jù)分析

為了驗(yàn)證最優(yōu)光滑數(shù)學(xué)模型的適用性，本文對光滑立管模型進(jìn)行試驗(yàn)與分析。本次試驗(yàn)在中國石油大學(xué)(華東)波流循環(huán)水槽中進(jìn)行(見圖8)，水槽尺度為16.0 m×0.8 m×1.4 m，最大水深1.0 m，流速0.05～0.70 m/s。

圖8 試驗(yàn)水槽Fig.8 Laboratory flume

試驗(yàn)選用外徑12 mm，長度為1 500 mm的有機(jī)玻璃管作為光滑立管模型，結(jié)構(gòu)特征數(shù)據(jù)見表3。試驗(yàn)過程中水深保持 700 mm，在立管模型中部750 mm處沿順流向和橫流向分別粘貼應(yīng)變片，以測量順流向和橫流向的立管響應(yīng)。試驗(yàn)采用東華測試DH8302動態(tài)信號測試分析系統(tǒng)采集應(yīng)變數(shù)據(jù)，采樣頻率設(shè)置為200 Hz。試驗(yàn)流速的選取參考王海青等[19]的研究結(jié)果，依據(jù)鎖振區(qū)立管的約化速度為4.5～7.0，設(shè)定水槽流速為0.110～0.233 m/s，流速間隔為0.03 m/s。

表3 立管結(jié)構(gòu)特征數(shù)據(jù)Tab.3 Riser structure characteristic data

以0.17 m/s流速下的試驗(yàn)為例，橫流向和順流向渦激振動應(yīng)變響應(yīng)如圖9所示。

圖9 渦激振動應(yīng)變響應(yīng) Fig.9 Strain response of VIV

將試驗(yàn)測量得到的應(yīng)變信號數(shù)據(jù)進(jìn)行EEMD分解，得到順流向和橫流向渦激振動信號的分解結(jié)果如圖10、圖11所示。

圖10 順流向信號EEMD分解Fig.10 EEMD decomposition of signals for IF

圖11 橫流向信號的EEMD分解Fig.11 EEMD decomposition of signals for CF

基于構(gòu)建的最優(yōu)降噪光滑算法，得到順流向和橫流向IMFs的存在度，如圖12所示。

圖12 渦激振動信號IMFs存在度Fig.12 The ratio of existence for IMFs in VIV

由圖12可知，順流向最優(yōu)化區(qū)間IMF3和IMF4存在度最高，橫流向最優(yōu)化區(qū)間IMF4存在度最高。因此判斷順流向有效本征模態(tài)函數(shù)為IMF3和IMF4，橫流向有效本征模態(tài)函數(shù)為IMF4。

為了判斷得到的有效本征模態(tài)函數(shù)是否具有物理意義，對其進(jìn)行傅里葉頻譜分析，結(jié)果見圖13。由頻譜分析結(jié)果可得：順流向?qū)?yīng)的本征函數(shù)模態(tài)IMF3頻率為6.25 Hz，IMF4頻率為3.125 Hz，橫流向?qū)?yīng)的本征函數(shù)模態(tài)IMF4頻率為3.125 Hz。

圖13 渦激振動頻譜圖Fig.13 Spectrum of VIV

從本文結(jié)果可知，順流向有兩個(gè)頻率參與到振動中去，分別是IMF3和IMF4；橫流向有一個(gè)頻率參與到振動中，對應(yīng)IMF4。從數(shù)據(jù)關(guān)系分析可得順流向IMF3的頻率是橫流向IMF4頻率的2倍，這和文獻(xiàn)[20]的研究結(jié)果(順流向主導(dǎo)頻率為橫流向主導(dǎo)頻率的二倍)是一致的。同時(shí)發(fā)現(xiàn)：順流向有效本征模態(tài)函數(shù)IMF4的頻率和橫流向IMF4頻率相等，這說明順流向振動受到橫流向振動的影響，這是由于漩渦脫落的過程中順流向受到漩渦脫落垂向的分力，橫流向受到漩渦脫落水平向分力的影響。

為進(jìn)一步驗(yàn)證結(jié)果的準(zhǔn)確性，利用IMFs能量分布進(jìn)行分析，結(jié)果如圖14所示。

圖14 IMFs能量譜Fig.14 Energy distribution of IMFs

由圖14可以看出，在0.17 m/s的流速下順流向和橫流向的能量對應(yīng)的頻段正好是相應(yīng)的有效本征模態(tài)函數(shù)對應(yīng)的頻段，說明本征模態(tài)函數(shù)選取是正確的。

為了驗(yàn)證該方法的通用性，對各工況進(jìn)行了同樣的分析，并將得到的有效IMFs及其能量占比匯總，如表4所示。

表4 各工況下有效IMFs能量占比Tab.4 Energy ratio of effective IMFs

從表4可以看出，在0.11～0.23 m/s的流速內(nèi)，基于構(gòu)造的最優(yōu)降噪光滑數(shù)學(xué)模型得到的順流向和橫流向有效IMFs占總能量的95%以上,證明該方法識別有效IMFs是準(zhǔn)確的。

4 結(jié) 論

文中采用EEMD算法對信號進(jìn)行分解，基于信號曲線的相似度、光滑度、存在度構(gòu)建了最優(yōu)降噪光滑數(shù)學(xué)模型，進(jìn)行立管渦激振動響應(yīng)參數(shù)識別。結(jié)果表明：

(1)構(gòu)建的最優(yōu)降噪光滑數(shù)學(xué)模型可以進(jìn)行有效IMF的優(yōu)選和信號參數(shù)識別。

(2)仿真信號識別表明，優(yōu)選的IMFs與原始信號頻率成分相關(guān)系數(shù)達(dá)到0.96以上，且頻率對應(yīng)相等，能量占原始信號總能量的97.98%。從相關(guān)度、頻譜分析、能量譜分析3個(gè)角度證明基于該算法選擇有效本征模態(tài)函數(shù)的方法是準(zhǔn)確的。

(3)渦激振動試驗(yàn)結(jié)果表明，經(jīng)過該算法優(yōu)選的IMFs能量占總能量的95%以上，順流向有效IMFs包含2個(gè)，橫流向有效IMFs包含1個(gè)，且順流向主要頻率是橫流向頻率的二倍關(guān)系?；跇?gòu)建的數(shù)學(xué)模型選取的IMFs含有明顯的渦激振動信息，這為后續(xù)渦激振動抑制、損傷分析的研究奠定了基礎(chǔ)。

Vol．41 No．12 2022