經(jīng)驗(yàn)?zāi)B(tài)分解算法（EMD）原理分析綜述目錄TOC\o"1-3"\h\u18855經(jīng)驗(yàn)?zāi)B(tài)分解算法（EMD）原理分析綜述 165631.1瞬時(shí)頻率的概念 186101.2固有模式函數(shù) 2324471.3經(jīng)驗(yàn)?zāi)B(tài)分解算法（EMD） 31.1瞬時(shí)頻率的概念一直以來(lái)，對(duì)信號(hào)頻率的定義都是建立在傅里葉變換的基礎(chǔ)上的，信號(hào)的頻率定義為跨越整個(gè)數(shù)長(zhǎng)度且恒定振幅的正弦或余弦函數(shù)的頻率。而這個(gè)定義就意味必然要求一段完整的信號(hào)才有頻率，無(wú)法定義瞬時(shí)頻率。自從有了希爾伯特變換之后，頻率有了新的定義，對(duì)于一個(gè)隨機(jī)序列x(t)，定義希爾伯特變換如下：y其中，P表示柯西（Cauchy）主值。X(t)和Y(t)可以構(gòu)成復(fù)共軛對(duì)，可以得到一個(gè)解析信號(hào)，Z(t)，如下：z(t)=x(t)+iy(t)=a(t)其中，a(t)=[x2a(t)為隨時(shí)間變化的幅度信息，θtω=但是，hilbert變換得到瞬時(shí)頻率并不一定是有意義的，從上面的公式可以看出，對(duì)于一個(gè)自然的非平穩(wěn)非線性信號(hào)，同一時(shí)刻的hilbert變換得到的瞬時(shí)頻率只有一個(gè)。如果一個(gè)自然信號(hào)本身是多種物理信號(hào)的疊加，那么理論上這個(gè)信號(hào)在同一時(shí)刻應(yīng)該具有多個(gè)瞬時(shí)頻率，通過(guò)hilbert變換得到的單一瞬時(shí)頻率在這時(shí)顯然是沒(méi)有物理意義的。為了得到有意義的瞬時(shí)頻率，要求希爾伯特變換處理的信號(hào)是單一模式。因此，希爾伯特變換的瞬時(shí)頻率定義一直沒(méi)能得到比較合理應(yīng)用，Huang等創(chuàng)新性設(shè)計(jì)出的EMD算法，將每個(gè)非線性非平穩(wěn)信號(hào)分解成固有模式分量（IMF），使得信號(hào)分解出的成分進(jìn)行希爾伯特變換得到的瞬時(shí)頻率是有意義的。1.2固有模式函數(shù)對(duì)希爾伯特變換，在物理上得到有意義的瞬時(shí)頻率的必要條件是：函數(shù)對(duì)稱于局部零均值，且有相同的極值點(diǎn)和過(guò)零點(diǎn)?；谙柌刈儞Q的特性，Huang等人在設(shè)計(jì)EMD分解出的IMF分量的時(shí)候，要求IMF分量滿足兩個(gè)條件：（1）在整個(gè)數(shù)據(jù)序列中，極值點(diǎn)的數(shù)量與過(guò)零點(diǎn)的數(shù)量必須相等，或最多相差一個(gè)；（2）在任何時(shí)間點(diǎn)上，被它的局部最大值與局部最小值定義的包絡(luò)均值必須是零。第一個(gè)條件可以滿足IMF分量在每個(gè)時(shí)間尺度內(nèi)只有單一的震動(dòng)模式，沒(méi)有其他諧波分量。第二個(gè)條件保證IMF分量的均值為零，能量分布在時(shí)間軸上下對(duì)稱，這樣，IMF分量就是一個(gè)調(diào)幅調(diào)頻的非線性非平穩(wěn)信號(hào)。下圖就是一個(gè)典型的IMF波形，來(lái)源于EMD算法分解MIT-BIH心率失常數(shù)據(jù)庫(kù)（MIT-BIHArrhythmiaDatabase）mitdb100數(shù)據(jù)產(chǎn)生的IMF5分量：圖2-1EMD分解mitdb100產(chǎn)生的IMF5Fig2-1EMDdecompositionofIMF5generatedbymitdb1001.3經(jīng)驗(yàn)?zāi)B(tài)分解算法（EMD）自然界中的非線性非平穩(wěn)信號(hào)往往是多個(gè)物理信源產(chǎn)生信號(hào)的疊加，包含不止一個(gè)震蕩模式，經(jīng)驗(yàn)?zāi)B(tài)分解就是通過(guò)后驗(yàn)的，自適應(yīng)的方式將不同信源產(chǎn)生的單一震蕩模式分開(kāi)。EMD分解的本質(zhì)是通過(guò)信號(hào)或者數(shù)據(jù)的時(shí)間尺度經(jīng)驗(yàn)的鑒別出數(shù)據(jù)中存在的固有震蕩模式，并以此為根據(jù)進(jìn)行數(shù)據(jù)分解。EMD算法的原理，首先第一步觀察數(shù)據(jù)，可以通過(guò)觀察數(shù)據(jù)定義局部時(shí)間尺度：同連續(xù)局部最大最小值確定的時(shí)間間隔或者是通過(guò)局部過(guò)零點(diǎn)確定的時(shí)間間隔。但是由于數(shù)據(jù)疊加的原因，高頻信號(hào)疊加在低頻信號(hào)上成為騎波，騎波的擾動(dòng)不能在過(guò)零點(diǎn)上檢測(cè)，通過(guò)局部過(guò)零點(diǎn)定義時(shí)間尺度局限性較強(qiáng)。所以必須采用極值之間的時(shí)間尺度來(lái)定義固有模式分量的時(shí)間尺度；這樣，不僅提供了精細(xì)的震蕩模式劃分，而且可以用于非零均值的數(shù)據(jù)。以相鄰局部極大值與極小值作為時(shí)間尺度，為將相同或相似時(shí)間尺度的信號(hào)從原始信號(hào)剝離，就需要找到這些振蕩波作為騎波所騎的趨勢(shì)分量。本論文以研究領(lǐng)域中經(jīng)典的非線性非穩(wěn)定信號(hào)，全球年地表氣溫異常數(shù)據(jù)（GlobalSurfaceTemperatureanomaly，GSTA），為例如下圖所示：圖2-2GSTA數(shù)據(jù)Fig2-2GSTAdata圖2-3中，紅色曲線為GSTA數(shù)據(jù)，綠色曲線為局部均值數(shù)據(jù)，也就是原始數(shù)據(jù)在每一個(gè)時(shí)刻的局部均值；在原始信號(hào)中，每一個(gè)相鄰的極大值和極小值之間的數(shù)據(jù)為一個(gè)時(shí)間尺度，此時(shí)間尺度中的信號(hào)是疊加在局部均值信號(hào)上的半個(gè)基本震蕩模式的周期，EMD方法的思路就是將這些時(shí)間尺度內(nèi)的基本震蕩模式與它所騎的局部均值分離，得到固有模式分量（IMF）。在Huang等提出的EMD方法實(shí)現(xiàn)中，原始數(shù)據(jù)極大值點(diǎn)和極小值點(diǎn)通過(guò)三次樣條插值得到原始信號(hào)的上下包絡(luò)，上下包絡(luò)相加求平均，得到局部均值。經(jīng)過(guò)多年的實(shí)踐檢驗(yàn)，這是一種比較合適的求局部均值的方式。得到局部均值之后，原始信號(hào)x(t)減去局部均值m1,理論上就可以得到最高頻率的IMF分量h1了：hh1如下圖所示：圖2-3理論上的IMF1分量Fig2-3ThetheoreticalIMF1component理論上，上圖就是一個(gè)標(biāo)準(zhǔn)的IMF分量，但實(shí)際上，這個(gè)h1一般不能滿足IMF分量的第二個(gè)條件，即這個(gè)h1的均值一般不為零，這是因?yàn)樵谌螛訔l進(jìn)行數(shù)據(jù)上下包絡(luò)擬合的時(shí)候，存在著過(guò)沖和欠沖現(xiàn)象，擬合的曲線的局部最大值和最小值相對(duì)于原始信號(hào)會(huì)產(chǎn)生偏移，使得局部峰值高于原始信號(hào)在該處的局部極大值，或者局部極小值低于原信號(hào)在該處的局部極小值。過(guò)沖和欠沖現(xiàn)象會(huì)使得分解出的IMF分量的極值點(diǎn)產(chǎn)生偏移和產(chǎn)生新的極值點(diǎn)，從而對(duì)IMF分量的特性和物理意義造成一定的破壞。為了減小過(guò)沖和欠沖帶來(lái)的影響，一方面可以通過(guò)改變求取局部均值的方式，得到更加精確的局部均值，目前效果較好的是建立在EMD方法基礎(chǔ)上的局域波概念；另一方面，Huang通過(guò)一個(gè)篩選過(guò)程，減少由于過(guò)沖和欠沖帶來(lái)的IMF波形不符合定義的問(wèn)題，減小雜波分量并且使得波形更加對(duì)稱。通過(guò)一步步的篩選，得到的最高頻單一模式分量的均值逐漸趨近于零。篩選過(guò)程如下：對(duì)上述得到的h1求上下極值點(diǎn)，采用三次樣條插值得到h1的上下包絡(luò)，上下包絡(luò)求平均得到h1的局部均值m11，h1減去m11得到新的篩選一次的IMF分量h11：一般情況下，h11也不滿足IMF定義要求，需要繼續(xù)進(jìn)行篩選：Huang等在定義篩選次數(shù)的時(shí)候經(jīng)驗(yàn)的判斷當(dāng)經(jīng)過(guò)十次這樣的篩選后，得到的h110就是第一個(gè)IMF分量。人們發(fā)現(xiàn)迭代次數(shù)不能太多，當(dāng)?shù)螖?shù)太多，波形在越來(lái)越規(guī)則的同時(shí)也會(huì)逐漸丟失調(diào)幅信息，成為純調(diào)頻的信號(hào)，IMF分量損失了它本來(lái)含有的物理意義。有學(xué)者研究幾種經(jīng)驗(yàn)?zāi)B(tài)分解的迭代停止準(zhǔn)則，用于自適應(yīng)的確定迭代次數(shù)。最簡(jiǎn)單的就是設(shè)置一個(gè)閾值A(chǔ)，當(dāng)h1k的均值小于閾值A(chǔ)時(shí)，就認(rèn)為h1k為篩選出的第一個(gè)IMF分量：meanHuang等人后來(lái)提出了用相鄰兩個(gè)篩選結(jié)果之間的標(biāo)準(zhǔn)差的大小作為篩選停止準(zhǔn)則，標(biāo)準(zhǔn)差（Standarddeviation，SD）的計(jì)算公式如下：SD=如果SD小于一定閾值時(shí)，篩選停止，得到最終的IMF分量。這里，閾值一般設(shè)置在0.2~0.3之間。下圖為篩選四次，十次和一百次的第一個(gè)IMF分量?？梢悦黠@看出，十次篩選的結(jié)果比四次篩選的結(jié)果更加關(guān)于時(shí)間軸上下對(duì)稱。尤其是在1900年附近，四次篩選的結(jié)果明顯存在趨勢(shì)波的影響，極小值點(diǎn)幾乎在零值附近，而十次篩選之后，相較而言極小值點(diǎn)與該位置的極大值點(diǎn)上下對(duì)稱。而當(dāng)篩選一百次后，IMF分量過(guò)分平滑，并且失去了幅度調(diào)制，信號(hào)幅值變化逐漸失去。圖2-4GSTA經(jīng)過(guò)三次篩選的IMF1分量Fig2-4IMF1componentsofGSTAafterthreescreenings圖2-5GSTA經(jīng)過(guò)十次篩選的IMF分量Fig2-5GSTA'sIMFcomponentsaftertenscreenings圖2-6GSTA經(jīng)過(guò)100次篩選得到的IMF1分量Fig2-6IMF1componentobtainedbyGSTAafter100screenings確定第一個(gè)IMF分量后，將它從原始信號(hào)中剝離，得到剩余分量進(jìn)行下一輪的篩選得到下一個(gè)低頻的IMF分量，整個(gè)經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）過(guò)程如下：（1）對(duì)一個(gè)信號(hào)x(t)，先找出他的所有極大值點(diǎn)和極小值點(diǎn)，然后利用三次樣條插值函數(shù)求得信號(hào)極大值和極小值的包絡(luò)，并求兩個(gè)包絡(luò)的均值，得到一條均值曲線m1(t)，令h1(t)=x(t)-m1(t)，得到h1(t)。（2）對(duì)得到的h1(t)進(jìn)行（1）操作，得到h2(t)，并重復(fù)10次，得到第一個(gè)IMF分量e1(t)，并令r1=x(t)-e1(t)，r1為對(duì)應(yīng)于第一個(gè)IMF分量的剩余分量。（3）對(duì)r1(t)進(jìn)行（1），（2），得到e2(t)，e3(t)……em(t)和r(t)，r(t)為剩余分量。結(jié)果，EMD將信號(hào)x(t)分解成如下形式：x其中k為IMF分量的階數(shù)，而m為總的IMF分量個(gè)數(shù)，r(t)為剩余分量。這里值得一提的是，最后的IMF分量如何確定，一般由于三次樣條插值的特性，當(dāng)極小值或者極大值點(diǎn)的個(gè)數(shù)小于3個(gè)點(diǎn)時(shí)，無(wú)法形成極大值或者極小值包絡(luò)，因此，當(dāng)最后的剩余分量不能產(chǎn)生上下包絡(luò)時(shí)，就意味著它是最后的剩余分量了；但是經(jīng)過(guò)Huang等人的大量經(jīng)驗(yàn)性的實(shí)驗(yàn)，發(fā)現(xiàn)IMF的個(gè)數(shù)N與信號(hào)的長(zhǎng)度L存在如下關(guān)系：N=在Huang等的EMD算法實(shí)現(xiàn)中，使用的是第二種方式，通過(guò)數(shù)據(jù)長(zhǎng)度計(jì)算出IMF分量的個(gè)數(shù)，剩余的能量分量就是剩余分量r。如下為Huang等設(shè)計(jì)的整個(gè)EMD算法流程：圖2-7EMD算法流程Fig2-7EMDalgorithmflow下圖為GSTA數(shù)據(jù)經(jīng)過(guò)EMD分解后，對(duì)得到的每個(gè)IMF分量做hilbert變換并計(jì)算相應(yīng)的頻率，得到GSTA數(shù)據(jù)HHT變換后的hilbert時(shí)頻分析圖，從圖中可以清晰的看到，在每個(gè)年份點(diǎn)上，都有對(duì)應(yīng)的頻率分量，而且不止一個(gè)頻率，每個(gè)時(shí)間點(diǎn)上的頻譜都是獨(dú)立的頻率，不存在頻譜的泄露，這是HHT算法與小波以及傅里葉變換最大的區(qū)別。圖2-8GSTAHHT時(shí)頻圖Fig2-8GSTAHHTtime-frequencydiagram圖2-9中由于GSTA時(shí)間軸以年為單位，頻率極低，所以時(shí)頻圖中縱軸使用hilbert變換得到θt綜上，EMD算法采用hilbert時(shí)頻譜，克服了頻率分辨率和時(shí)間分辨率之前權(quán)衡的難題；同時(shí)EMD算法分解信號(hào)是完全后驗(yàn)的，不需要先驗(yàn)的參數(shù)或者基，對(duì)于ECG信號(hào)這類存在個(gè)體較大差異的信號(hào)具有很好的普適性，近些年來(lái)在心電信號(hào)處理領(lǐng)域得到的廣泛應(yīng)用。張淼等將EMD分解與心電信號(hào)準(zhǔn)周期相結(jié)合，提出EMD-QPCE方法，在分解得到的IMF中消除母體心電并通過(guò)準(zhǔn)周期成分提取得到胎兒心電，與傳統(tǒng)方法相比，提升了靈敏度和準(zhǔn)確度[71]。ZhangMia