專利名稱：基于經驗模態分解的心肌超聲造影圖像生理參數測量方法
技術領域：
本發明涉及一種心肌超聲造影圖像的時間-強度曲線分析方法，具體涉及一種利用經驗模態分解進行時間-強度曲線的特征提取及利用希爾伯特邊際譜進行生理參數的測量方法。
背景技術：
心臟是人體的重要器官，隨著心臟每次收縮，將攜帶氧氣和營養物質的血流經主動脈輸送到全身，以供給各組織細胞代謝需要。由于脂質代謝不正常，血液中的脂質沉著在原本光滑的冠動脈內膜上，形成血栓。微血栓一旦形成，會使微血管發生嚴重狹窄甚至阻塞，產生微循環灌注異常，嚴重者發生心肌急性缺血、壞死、甚至危及生命。即便采用冠脈血運重建術進行治療，也大多存在心肌“低灌注”或“無灌注”異常，即心臟在缺血后再灌注時心肌存在微循環障礙。缺血心肌微循環灌注異常將嚴重影響受累心肌功能的恢復，導致病 人心臟功能減退。心肌超聲造影(MyocardialContrast Echocardiography, MCE)是一種影像新技術，通過將含有微氣泡的對比劑注入冠狀循環，來觀察二維超聲心動圖上的心肌顯影。由于微氣泡通過心肌時完全保持在血管內，且微泡的大小及變形性與紅細胞相當，因此可將其視作紅細胞流動的示蹤劑，故MCE可用來在跳動的心臟上估價心肌微循環狀況。目前臨床上常用定性法分析MCE，該方法雖簡便、直觀，但很大程度上依賴觀察者的經驗，缺乏客觀定量分析。超聲造影圖像感興趣區域的灌注強度曲線可以反映器官、組織的血流灌注情況。該技術可以在不破壞造影劑氣泡的前提下，實時動態地顯示組織的造影劑灌注，使精確的定量分析成為可能。通過對造影劑灌注的時間-強度曲線進行實時探測，獲得反映心肌血流灌注的生理參數，以便為心肌梗塞或心肌微循環診斷提供定量數據。經驗模態分解(Empirical Mode Decomposition,簡稱EMD)方法是由美國NASA的黃鍔博士提出的一種信號分析方法。它不會像傅立葉分解或小波分解那樣，預先設定好分解的基函數，而是會根據信號自身的時間尺度特征來進行信號分解。EMD方法在理論上可以應用于任何類型的信號的分解，尤其在處理非線性、非平穩信號上，具有非常明顯的優勢。利用心肌超聲造影圖像提取出的灌注強度曲線的內部時間尺度的變化做能量與頻率的解析，將該曲線展開成數個本征模態函數(Intrinsic Mode Function, IMF)的形式,再利用希爾伯特變換(Hilbert Transform, HT)獲得IMF的瞬時頻率和振幅,從而對提取出的時間-強度信號進行分析，從而獲得心肌灌注圖像的生理參數，為診斷提供可能。EMD是HHT (Hilbert-Huang Transform)算法的重要步驟,它認為信號是由許多個不同的模態組合而成的，可以是線性或非線性，分別對應不同過程，而且要求其相對于橫坐標軸是對稱的。要成為IMF內在模態函數有兩個條件
1)在整個函數中，極值點的數目與穿越零點的數目相等或者相差I;
2)在任何時刻，由局部極值包絡線所定義的包絡線局部均值為零。其中，第一個條件與傳統高斯平穩過程中窄頻寬的要求類似。第二個條件將整體性要求改變為局部性要求，使得瞬時頻率不會因為不對稱波形的存在而導致不必要的晃動。依托這兩個條件構建起來的EMD及HHT被認為是強有力地求解非線性、非平穩信號的自適應方法，是近年來對以傅立葉變換為基礎的線性及穩態譜分析的重大突破，并得到了廣泛的應用。臨床上經常使用擬合曲線對造影劑灌注的心肌參數進行定性測量，但是這樣的測量方法非常依賴于醫生的經驗，而且存在一定的測量誤差；同時由于擬合曲線只反映了時間-平均灰度曲線的變化趨勢信息，而曲線中蘊含的更多信息卻無法通過擬合曲線而表現，所以通常情況下健康的心肌區域和心肌梗塞區域的擬合曲線是基本類似的，測量的生理參數無法表現待診斷區域的性質，因此利用擬合曲線進行生理參數測量的臨床實用價值受到一定的限制。
發明內容
為了解決了當前心肌超聲造影圖像定性理解的主觀依賴性，本發明提出了一種基于經驗模態分解的心肌超聲造影圖像生理參數測量方法，利用經驗模態分解進行時間-強度曲線的特征提取，通過特征挖掘以及希爾伯特邊際譜，來進行心肌灌注情況下的生理參數測量方法，從而為基于心肌超聲造影圖像的心肌梗塞或心肌微循環障礙診斷提供了另一種分析途徑。本發明是通過以下技術方案實現的
一、圖像信號提取部分對心肌造影圖像的心肌區域按圓周順時針方向分為六個區域進行灌注信號提取，這六個區域分別代表前壁心肌區域、側壁心肌區域、后壁心肌區域、下壁心肌區域、后間隔心肌區域和前間隔心肌區域，分別提取每一個心肌分割區域的時間-強度曲線，曲線以橫時間為坐標，每個區域心肌的平均灰度值為縱坐標；
二、信號特征提取部分對提取出的時間-平均灰度值強度曲線進行EMD分解，得到擴展的信號本征模態函數，對得到的第一個本征模態函數進行Hilbert變換，得到瞬時譜參數，包括瞬時幅值、瞬時頻率和瞬時相位，對每個頻率點的幅值進行累加，得到Hilbert-Huang 邊際譜；
三、閾值判定部分利用Hilbert邊際譜的能量閾值對邊際譜進行閾值判定，如果邊際譜能量集中且大于閾值，可以得到造影劑灌注過程中的生理參數，包括呼吸頻率和心跳頻率參數。本發明具有如下優點
I、對心肌圖像進行六區域分割，可以兼顧臨床診斷及處理信息量的要求，利用區域內的像素平均值可以反映出該區域內的造影劑灌注情況，從而能夠利用時間-平均像素強度曲線來描述心肌超聲造影過程所攜帶的信息。經過該信息提取過程，可以很好地將造影過程中的生理參數信息反映出來。2、對提取出的時間-強度曲線進行經驗模態分解，得到能夠反映該曲線特征的本征模態函數，由于EMD會根據信號自身的時間尺度特征來進行信號分解，所以得到的IMF會含有豐富的頻率信息。利用HT的瞬時頻率概念得到希爾伯特譜，通過對其在時間域積分，利用邊際譜來研究在呼吸和心跳不同頻率內幅度值的變化情況，從而為心肌灌注的參數提取提供可能。
3、Hilbert邊際譜可以反映出信號能量在頻率軸上的分布情況，利用Hilbert邊際譜的能量閾值對邊際譜進行閾值判定，如果邊際譜能量集中且大于閾值，可以得到造影劑灌注過程中的呼吸頻率和心跳頻率參數，為心肌梗塞和心肌微循環提供定量依據。4、通過對SD大鼠的前壁心肌區域進行結扎后，通過心肌超聲造影進行觀察，并利用上述所給分析方法進行信號提取及特征提取。分析所得邊際譜可以看出，對應正常心肌區域的灌注生理參數，即呼吸頻率范圍為f 2Hz、心跳頻率為7 8Hz，與正常心肌情況下大鼠的基本生理參數一致。對于心肌梗塞區域，從灌注過程中得不到明確的心肌活動生理參數，與心肌壞死的前提一致，從而證明了本發明方法進行心肌灌注的生理參數測量的正確性。
圖I為基于經驗模態分解進行大鼠心肌超聲造影圖像時間-強度曲線的分析方法流程圖；
圖2為經驗模態分解流程 圖3為心肌的六區域分割 圖4為大鼠六個心肌區域的時間-強度曲線；
圖5為心肌區域I的時間-強度曲線的EMD分解結果；
圖6為心肌區域3的時間-強度曲線的EMD分解結果；
圖7為心肌區域I的IMFl中對應的各個瞬時參數；
圖8為心肌區域3的IMFl中對應的各個瞬時參數；
圖9為區域I的邊際譜；
圖10為區域3的邊際譜；
圖11為區域I的邊際譜閾值判斷 圖12為區域3的邊際譜閾值判斷圖。
具體實施例方式具體實施方式
一如圖I所示，本實施方式的基于經驗模態分解的心肌超聲造影圖像生理參數測量方法包括如下步驟首先對心肌造影圖像的心肌區域按順時針分為6個區域，分別提取每一個區域在心肌灌注過程中的時間-強度曲線；利用EMD方法對時間-強度曲線進行分解，得到對應不同頻率的本征模態函數；分別對本征模態函數做Hilbert變換，得到各種瞬時特征參數，并最終得到Hilbert邊際譜，根據能量閾值判斷，可以得到心肌造影圖像所對應的心肌灌注情況下的生理參數，包括呼吸頻率和心跳頻率范圍，從而為基于心肌超聲造影的心肌梗塞或心肌微循環診斷提供定量依據。具體步驟如下
步驟一圖像信號提取部分，為了滿足臨床上對于心肌梗塞診斷的要求，將心肌按圓周分為六個區域進行灌注分析，這六個區域分別代表心肌的不同部分。六個心肌區域按順時針分別標號為Region-I至Region-6,如圖3所示。Region-1,即順時針I點到3點的部分(I點為時鐘在I點時的位置，依次類推)對應實際心肌的前壁心肌區域；Region-2，即順時針3點到5點的部分對應心肌的側壁心肌區域；Region-3，即順時針5點到7點的部分對應心肌的后壁心肌區域；Region-4，即7點到9點的部分對應心肌的下壁心肌區域；Region-5，即9點到11點的部分對應心肌的后間隔心肌區域；Region-6，即11點到I點的部分對應心肌的前間隔心肌區域。心肌區域分割的目的是為了對每個區域進行分別處理，每個區域分別代表了心肌的不同部分，在臨床中通常僅對不同區域的心肌進行微循環灌注分析，通過計算每一個區域的平均灌注強度信息，來盡量反映該區域內的心肌灌注情況。實踐證明，將心肌進行分區處理能夠明顯提高灌注生理參數測量的準確性。時間-強度曲線提取是將分割出來的心肌區域進行時間-強度曲線繪制，由于心肌超聲造影圖像的幀順序反映了時間信息，而每一幀中的每個心肌區域中的平均像素強度反映了在該時刻上造影劑的灌注情況，所以以時間為橫坐標，每個區域心肌的平均灰度值為縱坐標可以做出時間-強度曲線，該曲線可以反映心肌內造影劑的灌注情況。直接使用平均灰度值繪制曲線的方法，能夠完全保留心腔區域灰度變化的信息，而且沒有引入心腔區域灰度變化的干擾、簡單直觀，后面的研究都是基于這種平均灰度值曲線的繪制方法進行的。
步驟二 信號特征提取部分，對提取出的時間-平均灰度值強度曲線進行EMD分解，得到擴展的信號本征模態函數。對經過提取的信號進行經驗模態分解(EMD)，得到擴展的特征曲線，用本征模態函數(IMF)表示。要成為IMF有兩個條件(I)極值及過零點相差不可以超過I個；(2)上下包絡線和為O。要完全嚴格的滿足這兩條中的第(2)條是不可能的，所以只能要求其盡量小。可以通過如下的分解方法得到IMF。經驗模態分解的流程圖如圖2所示，具體描述如下
為了 得到分解的本征模態函數IMF，定義. ^),其中1=1，…，L.代表了 IMF的序號，k=l，…K.代表了篩選次數。篩選過程從輸入信號開始，定義、(9 =忒4為了得到下一個IMF，我們將前一個IMF函數的余項作為輸入信號，對于給定信號χ( )的IMF篩選過程，包括以下步驟
(I)程序初始化，并搜索信號的振幅極大值及振幅極小值。(2)使用三次樣條插值得到上下包絡線，分別記為《+(O和.<(I)。并計算上下包絡的均值
、^(0+e~(t)
&^{ )=-；
2
(3)通過前一級信號與包絡均值計算出后一級信號
l^) =七觀(f);
檢查步驟(3)得到的信號是否滿足IMF的條件，即當包絡均值足夠接近零時迭代終止
卜.( )| < ε ；
(4)通過判定準則來判斷信號均值是否足夠趨近于零，如果不符合則重復步驟I至3。當滿足判定準則時，IMF函數Ci(Wj)定義為步驟3最后一次的結果，即
c 1(f) — ;
定義余項
= ^io ￠^) — cI (β)；
(5)計算下一個IMF需要重新初始化，將余項作為輸入信號，通過步驟I進行重復篩選，
權利要求
1.基于經驗模態分解的心肌超聲造影圖像生理參數測量方法，其特征在于所述方法包括如下步驟 一、圖像信號提取部分對心肌造影圖像的心肌區域按圓周順時針方向分為六個區域進行灌注信號提取，這六個區域分別代表前壁心肌區域、側壁心肌區域、后壁心肌區域、下壁心肌區域、后間隔心肌區域和前間隔心肌區域，分別提取每一個心肌分割區域的時間-強度曲線，曲線以橫時間為坐標，每個區域心肌的平均灰度值為縱坐標； 二、信號特征提取部分對提取出的時間-平均灰度值強度曲線進行EMD分解，得到擴展的信號本征模態函數，對得到的第一個本征模態函數進行Hilbert變換，得到瞬時譜參數，包括瞬時幅值、瞬時頻率和瞬時相位，對每個頻率點的幅值進行累加，得到Hilbert-Huang 邊際譜； 三、閾值判定部分利用Hilbert邊際譜的能量閾值對邊際譜進行閾值判定，如果邊際譜能量集中且大于閾值，可以得到造影劑灌注過程中的生理參數，包括呼吸頻率和心跳頻率參數。
2.根據權利I所述的基于經驗模態分解的心肌超聲造影圖像生理參數測量方法，其特征在于所述瞬時幅值、瞬時相位弓￠)和瞬時頻率A ￠)的表達式為
3.根據權利I所述的基于經驗模態分解的心肌超聲造影圖像生理參數測量方法，其特征在于所述邊際譜表達式為
4.根據權利I所述的基于經驗模態分解的心肌超聲造影圖像生理參數測量方法，其特征在于所述邊際譜有兩個大的能量幅值，且超過給定閾值之上，則該兩個能量閾值范圍對應的頻率點分別對應呼吸頻率點和心跳頻率點范圍。
5.根據權利I所述的基于經驗模態分解的心肌超聲造影圖像生理參數測量方法，其特征在于所述邊際譜只含有一個超過能量閾值的頻率段，則該能量閾值范圍對應的頻率點為呼吸頻率范圍。
全文摘要
基于經驗模態分解的心肌超聲造影圖像生理參數測量方法，涉及一種心肌超聲造影圖像的時間-強度曲線分析方法。本發明包括如下步驟一、對心肌區域按圓周分為六個區域進行灌注信號提取，分別提取每一個心肌分割區域的時間-強度曲線；二、對提取出的時間-平均灰度值強度曲線進行EMD分解，得到擴展的信號本征模態函數，對得到的第一個本征模態函數進行Hilbert變換，得到瞬時譜參數，對每個頻率點的幅值進行累加，得到邊際譜；三、利用邊際譜的能量閾值對邊際譜進行閾值判定。本發明為基于心肌超聲造影圖像的心肌梗塞或心肌微循環診斷提供了一種分析途徑，從而提高基于心肌超聲造影的生理參數測量的準確性，降低了醫生的主觀依賴性。
文檔編號G06T7/00GK102855623SQ20121025001
公開日2013年1月2日 申請日期2012年7月19日 優先權日2012年7月19日
發明者金晶, 杜國慶, 沈毅 申請人:哈爾濱工業大學