改進(jìn)的EMD方法及大壩強(qiáng)震反應(yīng)信號應(yīng)用分析

許亮華，郭永剛，杜修力

（1.中國水利水電科學(xué)研究院 工程抗震研究中心，北京 100048；2.北京工業(yè)大學(xué) 建筑工程學(xué)院，北京 100124）

1 研究背景

HHT（Hilbert-Huang transform）信號分析方法［1-2］是Norden E.Huang在1998年提出的一種全新的數(shù)據(jù)處理方法，方法核心由兩部分構(gòu)成，即經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒ǎ‥MD）和Hilbert變換。Hilbert變換得到的瞬時頻率、振幅以及時頻分布圖等可以反應(yīng)出數(shù)據(jù)的非平穩(wěn)變化規(guī)律。雖然HHT方法在非平穩(wěn)數(shù)據(jù)分析中得到廣泛應(yīng)用，但其方法本身仍然存在許多問題，如包絡(luò)線擬合以及邊界問題，沒有確切公式理論支持，模態(tài)混淆等。這些問題需要得到解決才能推動HHT方法的進(jìn)一步發(fā)展和應(yīng)用。大壩強(qiáng)震反應(yīng)信號是研究大壩結(jié)構(gòu)特性的珍貴數(shù)據(jù)［3～5］，可以從中獲得大壩的強(qiáng)震下非平穩(wěn)反應(yīng)規(guī)律、特性變化規(guī)律和損傷信息等有價值的信息，從而檢驗(yàn)和改進(jìn)大壩抗震設(shè)計，驗(yàn)證大壩抗震性能，保證大壩安全。大壩強(qiáng)震反應(yīng)數(shù)據(jù)是一種非平穩(wěn)數(shù)據(jù)，作者在應(yīng)用HHT方法分析處理中發(fā)現(xiàn)，常規(guī)的EMD方法對頻譜密集的數(shù)據(jù)分解不可避免產(chǎn)生模態(tài)混淆現(xiàn)象，導(dǎo)致HHT分析結(jié)果與實(shí)際偏差很大，分析結(jié)果無法作為結(jié)論進(jìn)行應(yīng)用。因此，應(yīng)用HHT方法首先必須抑制模態(tài)混淆現(xiàn)象發(fā)生。

2 常規(guī)EMD方法

信號的EMD過程［1-2］是一個將非平穩(wěn)信號進(jìn)行平穩(wěn)化處理的過程，分解目的是得到有近似窄帶特征且盡量滿足單分量函數(shù)［6］特征的固有模式函數(shù)（IMF），為下一步Hibert變換提供準(zhǔn)備。

EMD過程是一個篩選獲得IMF的過程，通過時間序列上下包絡(luò)的平均值確定瞬時平衡位置，進(jìn)而提取IMF。目前學(xué)者采用的EMD方法（為了與本文改進(jìn)方法名稱上區(qū)分，本文稱之為常規(guī)EMD方法）步驟為：首先搜索查找出信號x（t）的局部極大值和極小值，利用三次樣條曲線插值連接局部極大值和極小值點(diǎn)，分別得到極大值包絡(luò)線和極小值包絡(luò)線；其次對極大值包絡(luò)線線和極小值包絡(luò)線取平均，得到瞬時平均值m（t）；最后用原始數(shù)列x（t）減去瞬時平均值m（t），得到一個去掉低頻的新數(shù)列 h（t）。

檢查h（t）是否滿足IMF的2個條件：（1）極值點(diǎn)數(shù)目和過零點(diǎn)數(shù)目相等或最多相差1個；（2）在任意點(diǎn)，上下包絡(luò)線的平均值為0。

若h（t）滿足上述兩個條件，則將h（t）作為1個IMF；若不滿足，將h（t）序列按照上述步驟進(jìn)行重新“篩選”，直到滿足下式

滿足式（2）的兩個條件就得到第1個IMF，記為c1（t）。

將c1（t）從原始數(shù)列中分離出來。將r1（t）作為新的信號重復(fù)上述步驟處理。經(jīng)過處理后，原始數(shù)列x（t）依次被分解出n個IMF和一個殘余項rn（t）其中殘余項rn（t）代表了原始數(shù)列的趨勢。

通過EMD分解得到的IMF在特點(diǎn)上非常適合作Hilbert變換，從而構(gòu)造解析信號得到瞬時頻率。以上是常規(guī)EMD分解過程，具體流程見圖1。

3 基于濾波的EMD分解

在非平穩(wěn)信號中，間歇信號、噪聲信號和奇異信號存在，系統(tǒng)模態(tài)密集等因素會使EMD處理的IMF產(chǎn)生模態(tài)混淆，原因在于篩選得到的IMF并不能有效滿足單分量函數(shù)特征。在提高信噪比，剔除奇異點(diǎn)，以及采取手段抑制端點(diǎn)效應(yīng)的同時，如果信號保持窄帶特性，那么EMD篩選得到的IMF必然保持窄帶特性，這樣可以盡可能避免模態(tài)混淆或降低模態(tài)混淆程度。為了解決常規(guī)EMD分解會產(chǎn)生模態(tài)混淆問題，研究人員將濾波處理和EMD方法結(jié)合起來，根據(jù)信號的傅立葉譜頻帶特征進(jìn)行濾波處理，使得EMD分解時信號具有窄帶特性。

3.1 頻帶濾波與EMD方法結(jié)合應(yīng)用其他學(xué)者的思路 許多方法普遍都是采取預(yù)濾波處理的思路（圖3），對信號先行帶通濾波［7-8］。首先將信號根據(jù)不同頻帶參數(shù)帶通濾波獲得幾個帶通子信號，對每個子信號分別進(jìn)行EMD分解，篩選出各自的IMF；最后從所有帶通信號的IMF中提取認(rèn)為有用的IMF。

這種處理過程比較繁雜，生成了許多冗余IMF，失去了EMD方法的自適應(yīng)性優(yōu)勢。這種處理思路對于批量數(shù)據(jù)的信號處理是不可行的。因?yàn)槊總€通道信號都要選擇許多帶通濾波的參數(shù)，再加上人為篩選判斷IMF的過程，使求解效率大大降低，無法滿足工程應(yīng)用的要求。同時這樣處理結(jié)果可能無法保證分解信號的完備性和正交性。這種濾波方法特點(diǎn)都是EMD前進(jìn)行信號濾波，未觸及EMD過程。

3.2 本文改進(jìn)的EMD過程 基于前面濾波思路的不足，本文對于濾波思路進(jìn)行改進(jìn)，將濾波處理過程和EMD過程結(jié)合起來，在EMD分解的篩選過程中對信號進(jìn)行濾波處理。通過改進(jìn)，這種含濾波的EMD分解方法在保證了信號窄帶特性同時發(fā)揮了EMD方法自適應(yīng)分解的優(yōu)點(diǎn)。

3.2.1 篩選流程的改進(jìn) 首先對信號r（t）進(jìn)行高通濾波獲得R（t），將R（t）為EMD處理的信號，獲取到一個IMF后，獲取殘余項時采用r（t）-c1（t）=r1（t）處理。將r1（t）項作為新的r（t）重復(fù)以上過程，直到滿足篩選停止條件。改進(jìn)的EMD流程見圖2。

對于高頻噪聲較強(qiáng)且噪聲頻帶分布較寬的信號，可以先進(jìn)行低通濾波將信號噪聲濾除提高信噪比后，再采用以上方法進(jìn)行處理。

3.2.2 濾波參數(shù)設(shè)置過程信號r（t）高通濾波獲得R（t）過程中濾波參數(shù)的設(shè)定可以采用以下兩種方法設(shè)置。

（1）自適應(yīng)設(shè)置。這種方法的濾波參數(shù)是在EMD分解過程中自適應(yīng)設(shè)置并濾波，不需要人為確定濾波參數(shù)。根據(jù)每次篩選進(jìn)程中剩余成分r（t）的傅立葉譜，程序自行調(diào)整濾波參數(shù)。為便于描述，以一個仿真信號的傅立葉譜圖（見圖4）為例進(jìn)行說明。過程為：①r（t）傅立葉變換，求出傅立葉譜；②傅立葉譜平滑處理，得到平滑后的傅立葉譜；③查找出傅立葉譜中一系列的局部極大值和極小值頻率點(diǎn)，如例圖4中所示的A、B、C和D共4個局部極大值和a、b、c、d和e共5個局部極小值；④找出局部極大值中頻率最大的兩個點(diǎn)C和D點(diǎn)，通過C、D兩點(diǎn)確定出處兩點(diǎn)間的極小值點(diǎn)d點(diǎn)。將d點(diǎn)的頻率值作為本次高通濾波的最小頻率參數(shù)；⑤對r（t）采用Fourier濾波器高通濾波得到R（t）。

對傅立葉譜平滑處理是為了便于查找極大值極小值，平滑程度要根據(jù)具體信號處理調(diào)節(jié)。如果信號是密集頻率信號，譜平滑的次數(shù)要減少，否則，原本相鄰的兩個譜峰可能會合并成一個譜峰，濾波時便不能區(qū)分兩個頻帶。

（2）預(yù)先設(shè)定濾波參數(shù)。信號頻譜十分密集并且信噪比不高時，傅立葉譜平滑后劃分的頻帶可能也不滿足要求，EMD結(jié)果仍然有模態(tài)混淆現(xiàn)象，這時可以采用預(yù)先設(shè)定參數(shù)的方法。過程為：①對原始信號進(jìn)行傅立葉分析，獲得信號的傅立葉譜圖。根據(jù)譜圖特點(diǎn)先劃分頻帶，獲得每個頻帶的頻率參數(shù)，按頻率大小順序依次設(shè)置濾波參數(shù)并保存；②在每次篩選IMF前，程序依次從設(shè)置的濾波參數(shù)中提取濾波參數(shù)對r（t）進(jìn)行高通濾波。

4 實(shí)例分析

4.1 仿真信號分析 為討論改進(jìn)前后的分析效果，本文用最簡單的簡諧波構(gòu)建仿真信號（見圖5中No.1通道）。構(gòu)成的各頻率成分是15Hz的間歇簡諧波及2、5和6Hz簡諧波（分別見圖5中No.2-5通道），構(gòu)成的仿真信號含間歇信號具有密集模態(tài)的特點(diǎn)。

分別用常規(guī)EMD方法和本文改進(jìn)的EMD方法對仿真信號進(jìn)行分解。比較看出，常規(guī)EMD方法分解的imf1-4波形（見圖6）與仿真信號成分No.2-5波形（見圖5）明顯不同，說明已經(jīng)出現(xiàn)模態(tài)混淆。而本文改進(jìn)的EMD方法分解得到的imf1-4（見圖7），只是在端點(diǎn)處有輕微端點(diǎn)效應(yīng)，間歇信號兩端產(chǎn)生過渡信號。

對得到的EMD分解結(jié)果進(jìn)一步進(jìn)行Hilbert分析后，得到瞬時頻率曲線。常規(guī)EMD方法得到的瞬時頻率（見圖8）波動大，特別是構(gòu)成信號中5和6Hz簡諧波的頻率特征因?yàn)槟B(tài)混淆，出現(xiàn)了錯誤的結(jié)果。而改進(jìn)的EMD方法得到的瞬時頻率（見圖9），反應(yīng)了正確的頻率分布規(guī)律。這表明改進(jìn)的EMD方法有效解決了模態(tài)混淆問題。

本文瞬時頻率分布圖中每個IMF的瞬時頻率曲線的顏色深淺與相應(yīng)的IMF瞬時幅值相關(guān)，瞬時幅值最大時瞬時頻率用純黑色表示，瞬時幅值等于0時瞬時頻率用純白色表現(xiàn)。

4.2 應(yīng)用分析 目前在我國在水工建筑物特別是大壩上獲得的強(qiáng)震動反應(yīng)信號還不多，每個強(qiáng)震反應(yīng)信號都是研究建筑物地震振動特征的珍貴資料。對于這些資料需要用各種方法，各個角度去分析應(yīng)用，充分體現(xiàn)出這些數(shù)據(jù)的價值。

新豐江水電站［9］位于廣東省河源縣東江支流新豐江上，距廣州市東北約160km。1962年3月19日庫區(qū)發(fā)生了震級Ms為6.1級、震中烈度為8度的強(qiáng)烈地震。這次地震導(dǎo)致大壩右岸壩段頭部108m高程出現(xiàn)長達(dá)82m的水平裂縫，地震后對大壩裂縫進(jìn)行了修復(fù)加固。為了研究大壩斷裂加固效果和裂縫上下兩部分在地震作用下反應(yīng)差異，設(shè)立了大壩裂縫研究臺陣，以14、16、17號墩作為觀測對象，在108m高程裂縫上、下各20cm布設(shè)拾振器。

該大壩裂縫研究臺陣，在1989年11月26日記錄到震中距為2km的4.6級地震。本文提取這次地震中14壩段108m高程縫上和下、順河向的強(qiáng)震反應(yīng)信號（圖10）進(jìn)行HHT分析舉例。

4.2.1 數(shù)據(jù)預(yù)處理 通過信號的傅立葉譜判斷該大壩在此次地震中的卓越頻率在10Hz以下，為了減少高頻影響，對數(shù)據(jù)進(jìn)行HHT處理之前先進(jìn)行了10Hz以下的低通濾波，增加低頻的信噪比，濾波后波形見圖11。

4.2.2 HHT分析 對預(yù)處理后的強(qiáng)震反應(yīng)數(shù)據(jù)采用改進(jìn)的EMD分解得到每個通道信號的IMF結(jié)果（見圖12和圖13），對分解結(jié)果進(jìn)行Hilbert分析，計算得到各個IMF的瞬時頻率和瞬時幅值，繪制得到瞬時頻率分布（為觀測強(qiáng)震階段，1～5s作了放大，分別見圖14和圖15）和邊際譜（見圖16）。

EMD分解方法分解間歇信號時，由于分解是利用信號的包絡(luò)，而包絡(luò)對突變處的處理會導(dǎo)致分解的各個IMF結(jié)果在間歇信號開始和結(jié)束兩個時間點(diǎn)處產(chǎn)生過渡信號。由IMF結(jié)果中的過渡信號獲得的瞬時頻率特征并非實(shí)際信息真實(shí)存在的瞬時頻率，得到的瞬時頻率在間歇信號兩端會有較大的波動變化。結(jié)構(gòu)強(qiáng)震反應(yīng)信號在其非平穩(wěn)性表現(xiàn)之一表現(xiàn)為各個頻率成分振動幅值都增大過程和衰減過程，在幅值開始增大和幅值衰減到脈動水平時的兩個時間點(diǎn)，相當(dāng)于間歇信號的首尾兩個時間點(diǎn)，得到的IMF結(jié)果同樣會有過渡信號產(chǎn)生。研究瞬時頻率變化必須針對每個頻率成分的強(qiáng)震動階段，在強(qiáng)震動階段信噪比強(qiáng)，得到的瞬時頻率才能夠比較準(zhǔn)確、接近實(shí)際值，從而反應(yīng)出大壩結(jié)構(gòu)強(qiáng)震動過程中的非平穩(wěn)變化規(guī)律。

4.2.3 HHT分析結(jié)論

（1）通常由于放大效應(yīng)影響，同一次地震中高程高的監(jiān)測點(diǎn)監(jiān)測信號的幅值要比高程低的相同方向的監(jiān)測信號的幅值大。通過裂縫上下監(jiān)測點(diǎn)強(qiáng)震信號的波形振動幅度大小可以看出，裂縫上的振幅卻比裂縫下小（其它測點(diǎn)同樣是這種規(guī)律），這表明修復(fù)的大壩裂縫具有一定的緩沖吸能作用。放大效應(yīng)曲線在裂縫處會有所突變，并不是沿高程單純放大。

（2）大壩強(qiáng)震動的卓越頻率。從傅立葉譜分析裂縫上下的卓越頻率分別是5.811和5.859Hz。從邊際譜分析本次強(qiáng)震的大壩卓越頻率分別是5.762和5.859Hz。邊際譜和傅立葉譜分別見圖16、17。兩種譜圖的統(tǒng)計結(jié)果相近。頻率分析結(jié)果與過去觀測分析結(jié)果變化不大，這說明修復(fù)后大壩整體特性尚屬穩(wěn)定。

分析裂縫上下監(jiān)測點(diǎn)強(qiáng)震信號的邊際譜和傅立葉譜，發(fā)現(xiàn)各中心頻率略有不同，表明裂縫對大壩震動的頻率變化有一定影響。且各中心頻率的譜幅值比值不成相似的比例關(guān)系，說明裂縫對不同頻率的影響效果并不相同。具體裂縫對大壩振動的不同頻率成分影響是否存在一定的規(guī)律有待繼續(xù)深入研究。

（3）分別觀察裂縫上下強(qiáng)震信號的瞬時頻率分布圖，可以看出裂縫上下強(qiáng)震信號各頻率成分在強(qiáng)震階段頻率的瞬時頻率曲線沒有異常突變，基本處于水平波動變化，裂縫上下強(qiáng)震信號的各成分瞬時頻率分布規(guī)律相似。這些表明大壩在這次強(qiáng)震中處于彈性振動狀態(tài)，大壩修復(fù)后的裂縫在這次地震中處于安全狀態(tài)。

5 結(jié)語

本文對常規(guī)的EMD分解過程進(jìn)行改進(jìn)，在EMD分解進(jìn)程中采用濾波手段，這種改進(jìn)方式可以有效抑制模態(tài)混淆。通過仿真數(shù)據(jù)和實(shí)際大壩強(qiáng)震數(shù)據(jù)的分析表明改進(jìn)的EMD方法是有效可行的，可以用于大壩強(qiáng)震反應(yīng)數(shù)據(jù)的分析處理中。

［1］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proceeding of the Royal Society of Londom，1998，454：903-995.

［2］Huang N E，Shen Z，Long S R.A new view of nonlinear water waves：the Hilbert spectrum［J］.Annu.Rev.Flu?id Mech.，1999，31：417-557.

［3］DL/T5416-2009，水工建筑物強(qiáng)震動安全監(jiān)測技術(shù)規(guī)范［S］.

［4］蘇克忠，張力飛，等.大壩強(qiáng)震安全監(jiān)測［M］.北京：中國水利水電出版社，1995.

［5］郭永剛.水利水電工程強(qiáng)震監(jiān)測和強(qiáng)震監(jiān)測儀器［J］.地球物理學(xué)進(jìn)展，2005（2）：422-426.

［6］L.科恩.時—頻分析：理論與應(yīng)用［M］.白居憲譯，西安：西安交通大學(xué)出版社，1998.

［7］Yangj N，Lei Ying，Pan Shu-wen，et al.System identification of linear structures based on Hilbert-Huang spec?tral analysis：part 1 normal modes［J］.Earthquake Engng Struct Dyn.，2003，32（9）：1443-1467.

［8］郭淑卿.EMD的頻帶濾波篩分方法［J］.中國民航大學(xué)學(xué)報，2008，26（4）：30-33，39.

［9］陳厚群，蘇克忠，等.中國水工結(jié)構(gòu)重要強(qiáng)震數(shù)據(jù)及分析［M］.北京：地震出版社，2000.