EMD方法原理及其在ECG信號處理中存在的問題分析綜述目錄TOC\o"1-3"\h\u21799EMD方法原理及其在ECG信號處理中存在的問題分析綜述 1234951.1引言 1156781.2經(jīng)驗模態(tài)分解算法（EMD）原理 256391.2.1瞬時頻率的概念 2155021.2.2固有模式函數(shù) 3320041.2.3經(jīng)驗模態(tài)分解算法（EMD） 429311.3EMD算法的缺陷 12278751.4EEMD分解及存在的問題 16281081.5變分模式分解（VMD）及其缺陷 181.1引言1998年，Huang等人提出一種截然不同的信號分析方法，HHT（希爾伯特-黃變換）。它包含分解信號的經(jīng)驗模態(tài)分解（EMD，Empiricalmodedecomposition）和進行時頻分析的Hilbert變換[58]兩個部分。不同以往的信號分析手段，EMD算法將一段非線性非穩(wěn)定信號分解為一系列固有模態(tài)函數(shù)（IMF，Intrinsicmodefunctioncomponents），IMF作為原始信號的成分，不需要一組先驗的基對信號擬合。EMD分解出來的信號成分之所以稱之為固有模式分量，是因為這些IMF經(jīng)過Hilbert變換得到的頻率是有意義的，通過希爾伯特變換就可以得到信號的能量-頻率-時間分布，直觀看到信號的瞬時頻率值和對應頻率分量的能量大小。由此，HHT方法對頻率的定義與傅里葉變換和小波變換有了質(zhì)的不同，擺脫了基于先驗基測不準原理的束縛，在各行各業(yè)中得到廣泛應用，對很多工程應用中涉及的非線性問題，提供了不錯的解決方案，表現(xiàn)出了良好的處理非線性非平穩(wěn)信號的性能。王文波等基于模態(tài)單元比例萎縮的EMD方法設(shè)計出一套針對脈沖星信號的消噪算法，更加有效去除脈沖星信號中的噪聲的同時更好的保留了原始信號中有用的細節(jié)[67]；李珊等利用EMD算法對燃燒器火焰靜電信號進行去噪，獲得信號能量熵表征火焰帶電顆粒運動的復雜度，取得了不錯的分析效果[68]；周宣赤等利用HHT方法對2007年Notohanto地震數(shù)據(jù)進行分析，探究地震內(nèi)在的特點和規(guī)律，用于對地震的預防檢測中[69]；王翰等采用EMD方法和指標能量相結(jié)合的方法，對水輪機尾水管動態(tài)特征信息進行提取并進行驗證，結(jié)果表明該算法提取水輪機特征向量能力明顯[70]。但是經(jīng)過多年的算法實踐，經(jīng)驗模式分解（EMD）算法存在的一系列問題也從理論和應用兩個層面暴露了出來。本章簡單介紹EMD，EEMD和VMD算法原理，及其在心電信號處理方面的應用以及還存在的問題，并分析產(chǎn)生缺陷的原因。1.2經(jīng)驗模態(tài)分解算法（EMD）原理1.2.1瞬時頻率的概念一直以來，對信號頻率的定義都是建立在傅里葉變換的基礎(chǔ)上的，信號的頻率定義為跨越整個數(shù)長度且恒定振幅的正弦或余弦函數(shù)的頻率。而這個定義就意味必然要求一段完整的信號才有頻率，無法定義瞬時頻率。自從有了希爾伯特變換之后，頻率有了新的定義，對于一個隨機序列x(t)，定義希爾伯特變換如下：y其中，P表示柯西（Cauchy）主值。X(t)和Y(t)可以構(gòu)成復共軛對，可以得到一個解析信號，Z(t)，如下：z(t)=x(t)+iy(t)=a(t)其中，a(t)=[x2a(t)為隨時間變化的幅度信息，θtω=但是，hilbert變換得到瞬時頻率并不一定是有意義的，從上面的公式可以看出，對于一個自然的非平穩(wěn)非線性信號，同一時刻的hilbert變換得到的瞬時頻率只有一個。如果一個自然信號本身是多種物理信號的疊加，那么理論上這個信號在同一時刻應該具有多個瞬時頻率，通過hilbert變換得到的單一瞬時頻率在這時顯然是沒有物理意義的。為了得到有意義的瞬時頻率，要求希爾伯特變換處理的信號是單一模式。因此，希爾伯特變換的瞬時頻率定義一直沒能得到比較合理應用，Huang等創(chuàng)新性設(shè)計出的EMD算法，將每個非線性非平穩(wěn)信號分解成固有模式分量（IMF），使得信號分解出的成分進行希爾伯特變換得到的瞬時頻率是有意義的。1.2.2固有模式函數(shù)對希爾伯特變換，在物理上得到有意義的瞬時頻率的必要條件是：函數(shù)對稱于局部零均值，且有相同的極值點和過零點?；谙柌刈儞Q的特性，Huang等人在設(shè)計EMD分解出的IMF分量的時候，要求IMF分量滿足兩個條件：（1）在整個數(shù)據(jù)序列中，極值點的數(shù)量與過零點的數(shù)量必須相等，或最多相差一個；（2）在任何時間點上，被它的局部最大值與局部最小值定義的包絡均值必須是零。第一個條件可以滿足IMF分量在每個時間尺度內(nèi)只有單一的震動模式，沒有其他諧波分量。第二個條件保證IMF分量的均值為零，能量分布在時間軸上下對稱，這樣，IMF分量就是一個調(diào)幅調(diào)頻的非線性非平穩(wěn)信號。下圖就是一個典型的IMF波形，來源于EMD算法分解MIT-BIH心率失常數(shù)據(jù)庫（MIT-BIHArrhythmiaDatabase）mitdb100數(shù)據(jù)產(chǎn)生的IMF5分量：圖2-1EMD分解mitdb100產(chǎn)生的IMF5Fig2-1EMDdecompositionofIMF5generatedbymitdb1001.2.3經(jīng)驗模態(tài)分解算法（EMD）自然界中的非線性非平穩(wěn)信號往往是多個物理信源產(chǎn)生信號的疊加，包含不止一個震蕩模式，經(jīng)驗模態(tài)分解就是通過后驗的，自適應的方式將不同信源產(chǎn)生的單一震蕩模式分開。EMD分解的本質(zhì)是通過信號或者數(shù)據(jù)的時間尺度經(jīng)驗的鑒別出數(shù)據(jù)中存在的固有震蕩模式，并以此為根據(jù)進行數(shù)據(jù)分解。EMD算法的原理，首先第一步觀察數(shù)據(jù)，可以通過觀察數(shù)據(jù)定義局部時間尺度：同連續(xù)局部最大最小值確定的時間間隔或者是通過局部過零點確定的時間間隔。但是由于數(shù)據(jù)疊加的原因，高頻信號疊加在低頻信號上成為騎波，騎波的擾動不能在過零點上檢測，通過局部過零點定義時間尺度局限性較強。所以必須采用極值之間的時間尺度來定義固有模式分量的時間尺度；這樣，不僅提供了精細的震蕩模式劃分，而且可以用于非零均值的數(shù)據(jù)。以相鄰局部極大值與極小值作為時間尺度，為將相同或相似時間尺度的信號從原始信號剝離，就需要找到這些振蕩波作為騎波所騎的趨勢分量。本論文以研究領(lǐng)域中經(jīng)典的非線性非穩(wěn)定信號，全球年地表氣溫異常數(shù)據(jù)（GlobalSurfaceTemperatureanomaly，GSTA），為例如下圖所示：圖2-2GSTA數(shù)據(jù)Fig2-2GSTAdata圖2-3中，紅色曲線為GSTA數(shù)據(jù)，綠色曲線為局部均值數(shù)據(jù)，也就是原始數(shù)據(jù)在每一個時刻的局部均值；在原始信號中，每一個相鄰的極大值和極小值之間的數(shù)據(jù)為一個時間尺度，此時間尺度中的信號是疊加在局部均值信號上的半個基本震蕩模式的周期，EMD方法的思路就是將這些時間尺度內(nèi)的基本震蕩模式與它所騎的局部均值分離，得到固有模式分量（IMF）。在Huang等提出的EMD方法實現(xiàn)中，原始數(shù)據(jù)極大值點和極小值點通過三次樣條插值得到原始信號的上下包絡，上下包絡相加求平均，得到局部均值。經(jīng)過多年的實踐檢驗，這是一種比較合適的求局部均值的方式。得到局部均值之后，原始信號x(t)減去局部均值m1,理論上就可以得到最高頻率的IMF分量h1了：hh1如下圖所示：圖2-3理論上的IMF1分量Fig2-3ThetheoreticalIMF1component理論上，上圖就是一個標準的IMF分量，但實際上，這個h1一般不能滿足IMF分量的第二個條件，即這個h1的均值一般不為零，這是因為在三次樣條進行數(shù)據(jù)上下包絡擬合的時候，存在著過沖和欠沖現(xiàn)象，擬合的曲線的局部最大值和最小值相對于原始信號會產(chǎn)生偏移，使得局部峰值高于原始信號在該處的局部極大值，或者局部極小值低于原信號在該處的局部極小值。過沖和欠沖現(xiàn)象會使得分解出的IMF分量的極值點產(chǎn)生偏移和產(chǎn)生新的極值點，從而對IMF分量的特性和物理意義造成一定的破壞。為了減小過沖和欠沖帶來的影響，一方面可以通過改變求取局部均值的方式，得到更加精確的局部均值，目前效果較好的是建立在EMD方法基礎(chǔ)上的局域波概念；另一方面，Huang通過一個篩選過程，減少由于過沖和欠沖帶來的IMF波形不符合定義的問題，減小雜波分量并且使得波形更加對稱。通過一步步的篩選，得到的最高頻單一模式分量的均值逐漸趨近于零。篩選過程如下：對上述得到的h1求上下極值點，采用三次樣條插值得到h1的上下包絡，上下包絡求平均得到h1的局部均值m11，h1減去m11得到新的篩選一次的IMF分量h11：一般情況下，h11也不滿足IMF定義要求，需要繼續(xù)進行篩選：Huang等在定義篩選次數(shù)的時候經(jīng)驗的判斷當經(jīng)過十次這樣的篩選后，得到的h110就是第一個IMF分量。人們發(fā)現(xiàn)迭代次數(shù)不能太多，當?shù)螖?shù)太多，波形在越來越規(guī)則的同時也會逐漸丟失調(diào)幅信息，成為純調(diào)頻的信號，IMF分量損失了它本來含有的物理意義。有學者研究幾種經(jīng)驗模態(tài)分解的迭代停止準則，用于自適應的確定迭代次數(shù)。最簡單的就是設(shè)置一個閾值A(chǔ)，當h1k的均值小于閾值A(chǔ)時，就認為h1k為篩選出的第一個IMF分量：meanHuang等人后來提出了用相鄰兩個篩選結(jié)果之間的標準差的大小作為篩選停止準則，標準差（Standarddeviation，SD）的計算公式如下：SD=如果SD小于一定閾值時，篩選停止，得到最終的IMF分量。這里，閾值一般設(shè)置在0.2~0.3之間。下圖為篩選四次，十次和一百次的第一個IMF分量?？梢悦黠@看出，十次篩選的結(jié)果比四次篩選的結(jié)果更加關(guān)于時間軸上下對稱。尤其是在1900年附近，四次篩選的結(jié)果明顯存在趨勢波的影響，極小值點幾乎在零值附近，而十次篩選之后，相較而言極小值點與該位置的極大值點上下對稱。而當篩選一百次后，IMF分量過分平滑，并且失去了幅度調(diào)制，信號幅值變化逐漸失去。圖2-4GSTA經(jīng)過三次篩選的IMF1分量Fig2-4IMF1componentsofGSTAafterthreescreenings圖2-5GSTA經(jīng)過十次篩選的IMF分量Fig2-5GSTA'sIMFcomponentsaftertenscreenings圖2-6GSTA經(jīng)過100次篩選得到的IMF1分量Fig2-6IMF1componentobtainedbyGSTAafter100screenings確定第一個IMF分量后，將它從原始信號中剝離，得到剩余分量進行下一輪的篩選得到下一個低頻的IMF分量，整個經(jīng)驗模態(tài)分解（EMD）過程如下：（1）對一個信號x(t)，先找出他的所有極大值點和極小值點，然后利用三次樣條插值函數(shù)求得信號極大值和極小值的包絡，并求兩個包絡的均值，得到一條均值曲線m1(t)，令h1(t)=x(t)-m1(t)，得到h1(t)。（2）對得到的h1(t)進行（1）操作，得到h2(t)，并重復10次，得到第一個IMF分量e1(t)，并令r1=x(t)-e1(t)，r1為對應于第一個IMF分量的剩余分量。（3）對r1(t)進行（1），（2），得到e2(t)，e3(t)……em(t)和r(t)，r(t)為剩余分量。結(jié)果，EMD將信號x(t)分解成如下形式：x其中k為IMF分量的階數(shù)，而m為總的IMF分量個數(shù)，r(t)為剩余分量。這里值得一提的是，最后的IMF分量如何確定，一般由于三次樣條插值的特性，當極小值或者極大值點的個數(shù)小于3個點時，無法形成極大值或者極小值包絡，因此，當最后的剩余分量不能產(chǎn)生上下包絡時，就意味著它是最后的剩余分量了；但是經(jīng)過Huang等人的大量經(jīng)驗性的實驗，發(fā)現(xiàn)IMF的個數(shù)N與信號的長度L存在如下關(guān)系：N=在Huang等的EMD算法實現(xiàn)中，使用的是第二種方式，通過數(shù)據(jù)長度計算出IMF分量的個數(shù)，剩余的能量分量就是剩余分量r。如下為Huang等設(shè)計的整個EMD算法流程：圖2-7EMD算法流程Fig2-7EMDalgorithmflow下圖為GSTA數(shù)據(jù)經(jīng)過EMD分解后，對得到的每個IMF分量做hilbert變換并計算相應的頻率，得到GSTA數(shù)據(jù)HHT變換后的hilbert時頻分析圖，從圖中可以清晰的看到，在每個年份點上，都有對應的頻率分量，而且不止一個頻率，每個時間點上的頻譜都是獨立的頻率，不存在頻譜的泄露，這是HHT算法與小波以及傅里葉變換最大的區(qū)別。圖2-8GSTAHHT時頻圖Fig2-8GSTAHHTtime-frequencydiagram圖2-9中由于GSTA時間軸以年為單位，頻率極低，所以時頻圖中縱軸使用hilbert變換得到θt綜上，EMD算法采用hilbert時頻譜，克服了頻率分辨率和時間分辨率之前權(quán)衡的難題；同時EMD算法分解信號是完全后驗的，不需要先驗的參數(shù)或者基，對于ECG信號這類存在個體較大差異的信號具有很好的普適性，近些年來在心電信號處理領(lǐng)域得到的廣泛應用。張淼等將EMD分解與心電信號準周期相結(jié)合，提出EMD-QPCE方法，在分解得到的IMF中消除母體心電并通過準周期成分提取得到胎兒心電，與傳統(tǒng)方法相比，提升了靈敏度和準確度[71]。ZhangMiao等提出一種集總經(jīng)驗模式分解自適應閾值去噪方法（IEMD-ATD），基于IMF的能量和本征周期，開發(fā)了一種新的固有模式函數(shù)（IMF）分組方法。自適應劃定高頻噪聲，ECG信號和低頻噪聲能量為主的IMF分量，并推導了一種自適應濾波方法對高頻IMF分量進行降噪；結(jié)果表明，IEMD-ATD對ECG去噪在改善信噪比和相關(guān)系數(shù)方面具有更好的性能[72]。HuiXiong等使用隨機置換的方式處理EMD分解ECG信號得到的IMF分量，重構(gòu)后得到干凈的ECG信號，經(jīng)過MIT-BIH數(shù)據(jù)庫數(shù)據(jù)測試，有效抑制了噪聲，提高信噪比，并降低了均方誤差[73]。1.3EMD算法的缺陷雖然EMD算法擁有良好的分解信號的性能和時頻分析的性能，但是還是存在一系列影響算法性能的問題需要解決。EMD算法存在的主要缺陷如下：（1）首先由于EMD算法是經(jīng)驗性的模式篩選，沒有嚴格的數(shù)學公式推導，這使得EMD算法的理論化和進一步的研究困難重重；傳統(tǒng)方法基于的正交性等數(shù)學特性在EMD算法中推導和驗證存在不小困難。（2）局部均值的確定和迭代次數(shù)的確定。EMD算法要提取單一模式分量首先很重要的一點就是局部均值的概念和如何確定局部均值的問題，Huang等采用的三次樣條求上下包絡的方法目前是比較成熟和有效的方法，但存在過沖和欠沖的問題不得不進行迭代篩選，篩選次數(shù)的確定也是問題，目前篩選次數(shù)的確定都是根據(jù)IMF的兩個標準給出的；也有學者希望改進局部均值的求解方法，目前比較好的手段是局域波的方法，但這些方法也沒有嚴格的數(shù)學證明，或精準的物理意義。（3）端點效應。在EMD方法分解過程中，查找極大值點和極小值點時，端點也會被視為極大值極小值或者通過端點計算出虛擬的極大值極小值點。由于不是真實存在的極大值或極小值，三次樣條插值產(chǎn)生的曲線在端點處必然存在著誤差，這個誤差會在每一次篩選過程中向內(nèi)延伸，對整個數(shù)據(jù)的分解產(chǎn)生污染。對于這一現(xiàn)象，通過一個模擬信號進行檢驗：(a)原始信號(b)EMD分解圖2-9端點效應對比Fig2-9Endeffectcomparison從上圖可以看出，原始信號在左右端點處值為零，IMF2和IMF3在左端點處的值都是負值，IMF2在右端點的值是正值，而IMF3在左端點的值是負值，這是由于EMD分解過程中，在信號端點處理極值包絡時存在較大的誤差造成。（4）模態(tài)混疊問題。由于EMD算法得到單一模式分量的原理，在原始信號中存在不同模式的信號在同一層的時間尺度上的情況。這種情況與原始信號的特性有關(guān)，主要有以下三種情況：①某一模式分量信號存在間斷性；②信號的某一部分具有沖擊性，心電信號就屬于典型的帶有沖擊性大振幅的信號；③信號本身存在較強的噪聲干擾。模態(tài)混疊有可能是本來屬于不同模態(tài)的頻率成分出現(xiàn)在了同一模態(tài)中，也有可能是本屬于同一模態(tài)中的頻率成分出現(xiàn)在了不同模態(tài)中。不論是哪一種情況，都會造成IMF分量分解不準確，從而導致分解產(chǎn)生的IMF分量失去它本來的物理意義。下面通過構(gòu)建一個典型的間斷波模擬信號說明模式混疊現(xiàn)象對EMD方法造成的影響。構(gòu)建一個信號y，它有兩部分組成：(a)原始信號(b)EMD分解的模式混疊(c)y及上下包絡圖2-10EMD模式混疊及成因Fig2-10EMDmodealiasinganditscauses在上面的例子中，由于y2這種間斷波的存在，曲線上包絡不能反映正常趨勢量，導致局部均值函數(shù)產(chǎn)生很大改變。IMF1如上圖所示，除了兩個間斷高頻脈沖，顯然攜帶了低頻載波的能量成分，明顯發(fā)生了模態(tài)混疊現(xiàn)象，IMF2在5s和6s中間存在著高頻抖動，也存在著模態(tài)混疊，IMF分量已經(jīng)不能正確反映信號本身的物理意義了。下圖為EMD方法分解MIT-BIH心律失常數(shù)據(jù)庫100得到的IMF1，IMF2，IMF3分量。（a）IMF1（b）IMF2（c）IMF3圖2-11EMD分解mitdb100Fig2-11EMDdecompositionofmitdb100從上面可以得到，EMD分解mitdb100數(shù)據(jù)的三個IMF分量中，都包含了高頻噪聲分量和心電信號QRS波分量。無論是在時頻分析還是去噪方面，這樣的分解結(jié)果并不能達到更高精度的性能標準。1.4EEMD分解及存在的問題為消除經(jīng)驗模態(tài)分解存在的模式混疊問題，Huang和Wu等借助噪聲輔助分析的手段，對EMD算法進行改進，通過在原始信號中添加高斯白噪聲，使用EMD加入白噪聲的原始信號并進行集總平均，得到IMF分量，這種方法被稱為集總經(jīng)驗模式分解[74]（ensembleempiricalmodedecomposition，EEMD）。Patrick等發(fā)現(xiàn)，使用EMD方法分解高斯白噪聲時，EMD方法類似于一個二進制濾波器[75]，分解產(chǎn)生的IMF分量中心頻率以二倍的速度下降，頻率范圍以二倍速度上升，利用原始信號中加入高斯白噪聲，可以按二進制特性分解原始信號為一系列不同尺度的信號，抑制模態(tài)混疊的發(fā)生。EEMD分解步驟如下：（1）在原始信號x(t)中加入同長度的白噪聲N(t)，得到信號X(t)：X(t)=x(t)+N(t)(2（2）對得到的新信號X(t)進行EMD分解，得到一組IMF分量：X(t)=（3）重復（1）、（2）得到m組IMF分量，對這些IMF分量求集總平均：IMFIMFj(t)就是EEMD分解信號得到的單一模式分量。通過分解上面模式混疊模擬信號，可以驗證EEMD在消除模式混疊方面的能力：圖2-12EEMD分解混疊模擬信息Figure2-12EEMDdecompositionofaliasingsimulationinformation從EEMD分解混疊模擬信號來看，確實將間斷高頻分量與趨勢分量分開，IMF3，IMF4就是間斷高頻成分y2的成分；IMF7，IMF8是趨勢成分y1的成分。但是，從EEMD分解的高頻成分可以看出，EEMD雖然理論上可以通過集總平均將加入的噪聲消除，但在實際應用中并不能完全消除，導致在高頻IMF成分IMF3，IMF4中混入了噪聲分量，同時，在針對高脈沖的心電信號時，模式混疊的現(xiàn)象依然存在。圖2-13EEMD分解MIT-BIHmitdb100數(shù)據(jù)庫的IMF1分量Figure2-13EEMDbreaksdowntheIMF1componentoftheMIT-BIHmitdb100database上圖就是EEMD分解的MIT-BIH心率失常數(shù)據(jù)庫中的心電信號得到的IMF1分量，可以看出IMF1分量中不僅包含了原始心電信號中的高頻噪聲，同時包含了心電信號R波的能量分量。這就說明EEMD方法在面對心電信號這樣的陡脈沖信號時，集總的白噪聲能量不足以形成將噪聲與高能量脈沖分開的上下包絡信息。1.5變分模式分解（VMD）及其缺陷對于經(jīng)驗模式分解算法缺乏理論支撐的問題，KonstantinDragomiretskiy和DominiqueZosso在2014年提出變分模式分解[76]，將基于極值點三次樣條擬合進行單一模式分量篩選的模式通過數(shù)學表達代替，采用維納濾波等自適應濾波器的原理，將篩選過程變?yōu)樽顑?yōu)化問題，采用求