基于改進經驗小波變換的海洋平臺結構模態參數自動識別方法

冷建成, 刁凱欣, 龐 哲, 馮慧玉

(1. 東北石油大學 機械科學與工程學院, 黑龍江 大慶 163318;2. 東北石油大學 三亞海洋油氣研究院, 海南 三亞 572025)

海洋平臺作為油氣資源開發的主要設施,長期受到環境、人為、老化因素的影響,會產生不同程度的結構損傷,損傷的不斷累積會導致部分結構失效,對平臺的安全運營產生極大隱患。為此在平臺上布置結構健康監測系統,可以輔助管理人員把握平臺的整體結構參數及損傷演化規律,是保障平臺健康服役的有效手段之一。由于振動信號易于測量和采集,而且監測過程無需中斷結構的正常運行或需要特定的激勵設備,可實現實時監測并在整體上評估結構損傷,因此振動監測技術具有很高的應用價值,而模態參數識別是該技術的基礎[1]。

模態參數識別是指通過分析結構動力響應數據,獲取固有頻率、阻尼比、模態振型等參數的過程[2]。海洋平臺在環境激勵、作業激勵、溫度變化的作用下,其振動響應信號屬于幅值、頻率波動的非線性、非平穩信號,以傅里葉變換為基礎的頻域分析方法處理此類信號具有一定的局限性[3]。

相關學者提出希爾伯特-黃變換(Hilbert-Huang transform,HHT)[4]、小波變換(wavelet transform,WT)[5]等時頻域分析方法,在非線性、非平穩信號處理上具有很好的優越性。HHT中經驗模態分解(empirical mode decomposition,EMD)算法缺乏理論基礎,容易出現端點效應和模態混疊;WT擁有嚴格的理論基礎,但需要設定小波基函數等參數,且容易受到噪聲影響,自適應性較差[6]。

Gilles[7]在HHT算法與WT算法的基礎上,提出了經驗小波變換(empirical wavelet transform,EWT),該方法基于傅里葉頻譜劃分頻帶分割邊界,并構建小波濾波器將信號分解為若干子頻帶。由于EWT算法以傅里葉頻譜為基礎確定分割邊界,當信號信噪比較低時,頻譜的“毛刺”會影響邊界的劃分,導致方法失效。萬熹等[8]基于Burg算法構造自回歸(autoregression, AR)功率譜,用于替換傅里葉頻譜進行邊界劃分,該方法能夠抑制噪聲對頻譜模態峰值的干擾,提高了頻譜分割的準確度,成功應用于香港汀九斜拉橋的模態參數識別。Fang等[9]基于Yule-Walker算法構造AR功率譜進行頻譜劃分,解決了橋梁GNSS數據噪聲過大的問題,實現了某橋梁的模態參數識別。Amezquita-Sanchez等[10]使用多信號分類(MUSIC)算法進行頻譜分割,并將該方法應用于樂天大廈的模態參數識別。黃書亭等[11]提出了一種結合數學形態學的改進經驗小波變換方法,并將該方法應用到南京長江大橋的模態參數識別中。以上算法主要針對頻譜進行降噪來提升頻譜分割的準確性,但是頻譜降噪與邊界自動劃分的結合方面需要進一步深入研究。

為此,本文基于EMD算法的自適應性及EWT算法的模態分解正交性,引入奇異值譜和尺度空間算法,提出一種改進經驗小波變換方法,并結合隨機減量技術和希爾伯特變換,以實現環境激勵下結構的模態參數自動識別。

1 IEWT方法

1.1 經驗小波變換

經驗小波變換以傅里葉頻譜為基礎,確定分割邊界,并以此構造小波濾波器組將信號分解為一系列的調幅調頻單分量信號,其基本原理如下所述。

首先,預設單分量信號的數目為N,對信號x(t)進行傅里葉變換,并對其進行正則化,獲得傅里葉頻譜x(ω),其中ω∈[0,π];其次,遍歷頻譜x(ω)的局部極大值,并對其進行降序,提取前N個極大值,以相鄰極大值之間的最小值作為頻譜的分割邊界x=ωn;再次,根據分割邊界,以ωn為中心構造寬度為2τi過渡段,傅里葉頻譜被劃分為N個連續區間Λn,如式(1)所示,并以此構造經驗尺度函數φn(ω)與經驗小波函數ψn(ω),如式(2)～式(3)所示;最終,信號x(t)分解為調頻調幅組分,如式(4)～式(5)所示。

Λn=[ωn-1,ωn],n=1,2,…,N

(1)

(2)

(3)

式中,β(x)=x4(35-84x+70x2-20x3)。

(4)

(5)

1.2 改進經驗小波變換

由上述經驗小波變換原理可知,傅里葉頻譜的分割邊界主要依賴局部極大值的選擇,當所測信號信噪比較低,傅里葉頻譜中模態峰值附近存在若干緊密分布的“毛刺”,會影響分割邊界的確定,導致頻譜分割失效。

為減少噪聲對信號分解的影響,提出改進經驗小波變換(improved empirical wavelet transform,IEWT)方法,該方法結合了奇異值分解法及尺度空間法,在分割邊界確定方面對EWT方法進行改進,實現傅里葉頻譜降噪和分割邊界自適應劃分。

(1) 奇異值分解法

根據Brincker等[12]提出的頻域分解法,構造待測信號的功率譜密度矩陣,對其進行奇異值分解(singular value decomposition,SVD),如式(6)所示,能夠將多自由度系統轉換為單自由系統的疊加,在密集模態的模態參數識別中具有較好的效果。

(6)

由于奇異值對應單自由度系統的功率譜密度函數,且第一奇異值曲線包含了大量的模態信息,選擇該曲線構造奇異值譜,能夠有效降低噪聲對模態峰值干擾。

(2) 尺度空間法

針對頻譜中的模態峰值判定問題,Liutkus[13]提出了尺度空間峰值拾取(scale-space peak picking,SSPP)方法,通過使用不同尺寸的核(加權窗口)平滑數據,以迭代的方法確定峰值。

SSPP方法在迭代過程中,會建立一個C準則用于峰值估計,當發現潛在峰值時,C準則的值會變大,迭代結束后,可基于C準則直接進行峰值拾取,具體步驟如下:

步驟1初始參數輸入

設定待測信號v是維數為N×1的列向量,C準則的初始值為N×1的零向量,目標峰值數為K,峰值閾值為ρ,其中ρ∈[0,1]。

步驟2峰值估計準則初始化

檢出信號v中所有局部極大值,存儲于集合P中,如式(7)所示。若v(n)為局部極大值(n為位置索引),則v(n)應大于相鄰數據,如式(8)所示。

P←localmaxima(v)

(7)

v(n)∈P?v(n)=max[v(n-1),v(n),v(n+1)]

(8)

由于待測信號v未經平滑處理,包含較多噪聲,且集合P中的局部極大值點是基于待測信號v求得,需要進一步處理。在第一次迭代中,將所有的局部極大值存儲于C準則中,如式(9)所示。

?p∈P,C(p)←v(p)

(9)

另外,將所選峰值也存儲于集合O中,即O←P。

步驟3峰值迭代求解

選擇一個長度尺度s對待測信號v進行平滑,降低噪聲對信號的影響,獲得新的平滑信號v。在實際應用中,數據平滑可由信號v與加權窗口(如長度為s的標準化漢明窗)的卷積實現。

對信號進行數據平滑后,重新檢測局部極大值,得到新的集合P,如式(7)及式(8)所示。將新的局部極大值與集合O中已確定的極大值點相關聯,若某些局部極大值點被再次選中,將在后續處理中增加這些點的C準則。

由于數據平滑后,局部極大值點會發生漂移,需要確定集合O中哪些元素最接近集合P,從而建立近鄰集合I,如式(10)所示。

(10)

如式(10)所示,確定了集合O與集合P的近鄰集合I,該集合即為平滑后的被再次被選中的局部極大值點。依據近鄰集合I,增加對應點的C準則,如式(11)及式(12)所示。

?p∈P, ΔC(I(p))←v(I(p))s2

(11)

?n,C(n)←C(n)+ΔC(n)

(12)

不斷重復上述過程,并每次迭代都增加尺度s,并增加C準則的值,一般迭代30次即可。

步驟4峰值結果輸出

將C準則進行降序排序,依據預設的目標峰值數K或峰值閾值ρ進行峰值選擇。

若求解峰值數K已知,選取C準則中的前K個值對應的點作為峰值輸出結果;若峰值閾值ρ已知,將C準則進行歸一化處理,選取大于閾值ρ的點作為峰值輸出結果。

通過以上方法的改進,EWT方法中傅里葉頻譜被替換為信噪比較高的奇異值譜,尺度空間法自適應地確定模態峰值與分割邊界,結合小波濾波器組將信號分解為噪聲分量與若干固有模態函數(intrinsic mode function,IMF)分量。

1.3 模態參數自動識別

在對結構進行實時監測分析時,需要重視數據分析的實時性及數據分析的準確性。對于實時監測的數據流,可采用滑動數據窗的形式對數據流進行分段截取,再進行后續的模態分析,其基本原理如圖1所示。其中,分析間隔Δt與數據的實時性有關,每次分析數據長度t與該段數據頻譜的分辨率有關。選擇合適的分析間隔Δt和每次分析數據長度t對于數據的實時性及后續分析的準確性尤為重要。

圖1 滑動數據窗原理Fig.1 Principle of sliding data window

實時監測數據流經過滑動數據窗分段后,采用IEWT方法獲得該段數據的IMF分量,進而采用隨機減量技術與希爾伯特變換相結合的方法,實現頻率及阻尼比等模態信息的識別。

(1) 隨機減量技術

隨機減量技術(random decrement technique,RDT)是一種從環境激勵下結構響應信號中提取自由衰減響應的處理方法[14]。

設y(t)為隨機激勵下某線性系統的單自由度振動響應信號,選取一個固定值A對y(t)的幅值進行截取(A的值通常為信號均方差的1.5倍),A與響應信號的交點為ti(i=1,2,…,N),獲取N段以ti為起點,長度為s的時間序列。

對以上N段時間序列進行疊加平均即可獲得y(t)的自由衰減響應信號x(t),如式(13)所示。

(13)

由于IEWT方法獲得的IMF分量只包含單一頻率,屬于單自由度振動響應信號,可采用RDT方法獲取對應的自由衰減響應信號。

(2) 希爾伯特變換

希爾伯特變換(Hilbert transform,HT)可對自由衰減信號進行處理,獲得幅值曲線和相位曲線,實現模態頻率及阻尼比的識別。

設xj(t)為RDT方法獲得的某自由衰減響應信號,如式(14)所示。

xj(t)=A0je-ξjωjtcos(ωdjt+φj)

(14)

對式(14)進行希爾伯特變換,得到相位函數θj(t)及幅值函數Aj(t),如式(15)～式(16)所示。

θj(t)=ωdjt+φj

(15)

Aj(t)=A0e-ξjωjt

(16)

對相位函數θj(t)求導可得第j階有阻尼圓頻率ωdj,對于實際的工程結構,阻尼比通常小于5%,故ωdj≈ωj,進而獲得第j階模態頻率fj,如式(17)所示。

(17)

由于幅值函數Aj(t)為指數函數,構造指數衰減曲線對幅值函數進行擬合,如式(18)所示。

y(t)=A0ebt,b=-2πfjξj

(18)

曲線擬合之后可獲得b的取值,進而獲得第j階阻尼比,如式(19)所示。

(19)

綜上,提出了基于IEWT的模態參數自動識別方法,基于實時監測工況下的數據流,首先采用滑動數據窗方法進行數據分段;其次采用IEWT將分段數據分解為若干IMF分量;再次采用RDT&HT方法對各IMF分量進行模態參數識別;最終獲得信號各階模態的頻率與阻尼比等信息,其基本流程如圖2所示。

圖2 模態參數自動識別流程圖Fig.2 Flow chart of automatic identification of modal parameters

2 仿真信號分析

2.1 自由振動響應信號

構造三自由度結構的自由振動響應信號s(t),該信號由三個單分量信號s1(t)、s2(t)、s3(t)及信噪比為10 dB的高斯噪聲n(t)組成,如式(20)所示。

(20)

式中:Ai為幅值;fi為頻率;ξi為阻尼比;θi為相角。各單分量信號的具體參數如表1所示。

表1 各分量信號參數Tab.1 Signal parameters of each component

以100 Hz的采樣頻率獲取30 s的數據進行波形展示,如圖3所示。

圖3 自由振動響應信號時程圖Fig.3 Time history chart of free vibration response signal

首先采用EWT算法對自由振動響應信號的傅里葉頻譜進行分析,如圖4所示。

圖4 傅里葉頻譜分割邊界Fig.4 Segmentation boundaries of Fourier spectrum

由圖4可見,由于噪聲干擾,傅里葉頻譜第二個峰值被過度分割,導致頻譜分割失效。下面采用IEWT方法進行分析:首先計算信號的互功率譜矩陣,對其進行奇異值分解,獲得該信號的奇異值譜;其次,使用尺度空間法對奇異值譜進行峰值識別,并確定頻譜的分割邊界;最后將分割邊界映射至傅里葉頻譜中,計算結果如圖5和圖6所示。

圖5 奇異值譜及分割邊界Fig.5 Singular value spectrum and boundaries

圖6 傅里葉頻譜分割邊界Fig.6 Segmentation boundaries of Fourier spectrum

相比圖4中的傅里葉頻譜,圖5中的奇異值譜噪聲較低,各階模態峰值突出,更便于進行峰值的識別與分割邊界的劃分。由圖6中映射的分割邊界可知,傅里葉頻譜將被分割邊界準確劃分,與圖4相比可見,IEWT方法具有很好的分割效果。

根據IEWT方法獲取的分割邊界,構造小波濾波器組對傅里葉頻譜進行分割,并對分割結果進行逆傅里葉變換獲得三個IMF分量,如圖7所示。

(a) 模態1

(b) 模態2

(c) 模態3圖7 模態分量Fig.7 Modal components

由以上自由振動響應信號的分析結果可知,IEWT方法能夠對低信噪比的密集模態信號進行頻譜分割邊界的自適應劃分,進而精準獲取信號的IMF分量。

2.2 ASCE Bechmark模型信號

為進一步評估IEWT方法在多模態信號處理中的效果,采用IASC-ASCE SHM Benchmark模型進行驗證[15],該模型為4層2×2跨的鋼結構模型。根據參數設置,選用12自由度無損結構進行分析,其中阻尼比為1%,以高斯白噪聲模擬外部激勵。選取x方向1#、2#、3#及4#傳感器數據進行分析,各傳感器分別布置于每層桁架的中點處,如圖8所示。

圖8 Benchmark模型結構示意圖(m)Fig.8 Structure diagram of Benchmark model (m)

基于Johnson等提供的MATLAB程序模擬該結構的動力響應,這里僅分析4#傳感器的時程曲線,如圖9所示。其中采樣頻率為1 000 Hz,采樣時長為30 s。

圖9 傳感器數據時程圖Fig.9 Time history chart of sensor data

分別采用EWT方法和IEWT方法對該信號進行處理,其分割邊界如圖10和圖11所示。

圖10 傅里葉頻譜分割邊界Fig.10 Segmentation boundaries of Fourier spectrum

圖11 奇異值譜及分割邊界Fig.11 Singular value spectrum and boundaries

由圖10可見,基于EWT的頻譜分割方法無法準確檢測到所有邊界,其中前兩階頻率被劃分至同一個區間,發生了模態混疊;在50 Hz附近的頻譜中,由于噪聲干擾,EWT方法密集在此處劃分了若干邊界,傅里葉頻譜被過度分割出若干無意義的頻帶。

相對應地,由圖11可知,IEWT方法將奇異值譜分割為若干頻帶區間,各階頻率被準確劃分至不同的頻帶中。第1、9頻帶用于隔離噪聲,不參與模態參數識別。對比圖10和圖11可知,IEWT方法具有很好的自適應性,能夠避免發生模態混疊及過度分割的問題。

將圖11中的分割邊界映射至傅里葉頻譜中,如圖12所示。構造小波濾波器組將傅里葉頻譜劃分至7個模態分量,其中前3階模態分量如圖13所示。

圖12 傅里葉譜分割邊界Fig.12 Segmentation boundaries of Fourier spectrum

(a) 模態1

(b) 模態2

(c) 模態3圖13 前三階模態分量Fig.13 The first three modal components

采用RDT方法獲取各階模態分量的自由衰減響應信號,其中前三階模態自由衰減響應如圖14所示。結合HT獲取幅值函數與相位函數,使用指數衰減擬合和直接求導方法得出各階模態頻率及阻尼比參數。

(a) 模態1

(c) 模態3圖14 前三階模態自由衰減響應Fig.14 Free attenuation response of the first three modals

為驗證基于IEWT的模態參數識別方法的準確性,采用WT方法[5]及AR-EWT方法[8]進行橫向對比,三種方法的模態參數識別結果匯總于表2中。為進一步評估各方法的模態參數識別精度,使用相對誤差對各參數進行評估,如式(21)所示。

(21)

式中:Er為相對誤差;x為識別值;X為理論值。

由表2可知,對比WT方法與AR-EWT方法,IEWT方法能夠全部準確識別出各階模態的頻率及阻尼比信息。對于固有頻率,IEWT方法的相對誤差小于1.40%,與理論值基本一致;對于阻尼比,IEWT方法對于前兩階的識別精度較高,且阻尼比識別與理論值較為接近。

3 基于IEWT的海洋平臺結構模態參數識別

為驗證在實際應用中的可行性,選用實驗室自升式海洋平臺模型為研究對象,基于IEWT方法對其進行模態參數識別。

3.1 數據采集

該海洋平臺模型基于某自升式海洋平臺,以1∶40比例幾何縮放而成,整體結構由三個桁架樁腿和甲板組成,如圖15所示。甲板上布置6個加速度傳感器,其中1#、3#、5#傳感器對應x方向,2#、4#、6#傳感器對應y方向,傳感器布置位置如圖16所示。

圖15 海洋平臺模型Fig.15 Offshore platform model

圖16 傳感器布置圖Fig.16 Layout of acceleration sensors

傳感器選用891-2型低頻拾振器,數據采集設備選用cDAQ-9189機箱與NI 9231采集卡,采樣頻率為400 Hz。

3.2 模態參數識別

選取x方向傳感器數據進行分析,數據截取時長為60 s,其中1#傳感器數據如圖17所示。

圖17 1#傳感器數據時程圖Fig.17 Time history chart of sensor 1 data

基于IEWT方法獲取該信號的奇異值的分割邊界,如圖18所示。奇異值譜被劃分為5個頻譜區間,其中第1、5個頻譜區間為噪聲,對第2、3、4頻譜區間為模態區間。

圖18 奇異值譜及分割邊界Fig.18 Fourier spectrum and boundaries

將奇異值譜中的分割邊界映射至傅里葉頻譜中,如圖19所示。構造小波濾波器組進行頻譜分割,獲得三個模態分量,結合RDT及HT方法計算各階模態頻率及阻尼比,如表3所示。

圖19 傅里葉譜分割邊界Fig.19 Segmentation boundaries of Fourier spectrum

表3 模態參數識別結果Tab.3 Modal parameter identification results

由表3可知,IEWT方法成功識別出了前三階模態,其中模態頻率分別為4.18 Hz、4.85 Hz、6.35 Hz。為驗證該方法識別結果的可靠性,使用基于協方差的隨機子空間(stochastic subspace identification-covariance, SSI-COV)方法和自動頻域分解(automated frequency domain decomposition, AFDD)方法進行對比,不難發現,三種方法對于模態頻率的識別結果較為接近,也驗證了IEWT方法在實際應用中的可靠性。

海洋平臺等大型結構進行模態分析時,通常在環境激勵下進行,由于這種分析方法不會中斷結構的正常運營,又稱工作模態分析。在對結構進行工作模態分析時,結構阻尼比的識別結果會受到環境振動水平及分析時長的影響,具有一定的離散性[16]。目前大型結構的阻尼機理不明確,在低水平環境振動下增大采樣時長能夠降低阻尼比識別結果的離散性,而實際的海洋平臺環境激勵較為復雜,難以控制。相較于阻尼比,結構的模態頻率識別精度較高、離散性較低、可靠性強,可以作為結構的損傷指標。

此外,與SSI-COV等方法相比,IEWT具有較快的運算速度,完成60 s時長數據的模態參數識別用時僅為0.51 s,因此該方法在結構健康監測中具有較好的應用價值,可用于結構模態參數自動識別,為后續的結構損傷預警及損傷評估提供依據。

4 結 論

針對EWT算法在含噪信號傅里葉頻譜中分割邊界容易誤判的缺陷,本文利用奇異值譜在處理含噪信號中的優勢,結合尺度空間算法實現分割邊界的自適應劃分,對EWT算法進行改進。通過自由振動響應信號、Benchmark模型信號與自升式平臺模型實測信號對算法進行驗證,得到以下結論:

(1) 基于奇異值譜和尺度空間算法,對EWT算法進行改進,提出一種能夠實現模態參數自動識別的IEWT方法。

(2) 通過將IEWT方法應用于三自由度的自由振動響應信號和Benchmark模型信號,并分別與EWT、WT、AR-EWT方法進行對比,結果表明IEWT方法能夠實現頻譜的自適應劃分和模態參數的準確識別。

(3) 基于室內自升式海洋平臺模型,將IEWT方法應用到實際數據的模態參數識別中,對比SSI-COV及AFDD方法的識別結果,結果表明IEWT方法在實際應用中具有可靠性強、計算速度快的優點,在其他大型工程結構健康監測中具有一定的推廣應用價值。