基于EEMD有效成分優(yōu)選的超聲多普勒血流測速法

楊本賢，何冰冰，張榆鋒**，聶建云，姚瑞晗，劉亞杰

(1.云南大學(xué) 信息學(xué)院，云南 昆明 650500；2.昆明醫(yī)科大學(xué)第三附屬醫(yī)院（云南省腫瘤醫(yī)院），云南 昆明 650118)

大量臨床研究結(jié)果顯示[1-6]，頸動(dòng)脈粥樣硬化是導(dǎo)致腦血管疾病的主要原因，多普勒彩色血流成像技術(shù)可探究動(dòng)脈粥樣硬化引起的徑向血流速度剖面的變化[1].血流速度分布剖面估計(jì)關(guān)鍵在于對血流多普勒信號(hào)進(jìn)行時(shí)頻分析的準(zhǔn)確性[7].傳統(tǒng)傅里葉變換無法同時(shí)展現(xiàn)信號(hào)在時(shí)域和頻域的完整信息[8]，所以短時(shí)傅里葉變換(Short-time Fourier Transform,STFT)是估計(jì)多普勒血流信號(hào)頻率的常用方法[9].但STFT法使用固定窗函數(shù)，無法兼顧時(shí)間和頻率分辨率[10].小波的多分辨率分析具有良好的空間域和頻域局部化特性，適用于非平穩(wěn)信號(hào)[11-12].文獻(xiàn)[13-16]分別使用小波變換對多普勒信號(hào)進(jìn)行分析并取得了良好的效果.小波變換需預(yù)先確定小波基函數(shù)和分解尺度，但是缺乏系統(tǒng)規(guī)范的最佳小波基選取方法.與小波變換相比，Huang等[17]提出經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)是一種完全基于數(shù)據(jù)驅(qū)動(dòng)的自適應(yīng)分解方法，更適合處理非平穩(wěn)信號(hào).但EMD本身存在如模式混疊、端點(diǎn)效應(yīng)、停止條件等不足[18-19].2009年Wu等[20]對EMD進(jìn)行改進(jìn)，提出集合經(jīng)驗(yàn) 模 態(tài) 分 解(Ensemble Empirical Mode Decomposition,EEMD).EEMD法能有效處理非線性、非平穩(wěn)的信號(hào)，具有直接、后驗(yàn)和自適應(yīng)等優(yōu)點(diǎn)[21-22].因此，為提高速度剖面測量精度，本文提出基于集合經(jīng)驗(yàn)?zāi)B(tài)分解，并根據(jù)歸一化波動(dòng)指數(shù)(Normalized Fluctuation Index, NFI)選擇有效血流成分的新方法(EEMD_N).在仿真實(shí)驗(yàn)中搭建血流模型并獲取血流多普勒信號(hào)，對血流多普勒信號(hào)進(jìn)行EEMD分解，得到本征模態(tài)函數(shù)(Intrinsic Mode Function,IMF)組；然后計(jì)算IMF的NFI，使用傅里葉函數(shù)擬合法找到NFI最優(yōu)閾值，根據(jù)NFI最優(yōu)閾值選擇有效血流成分計(jì)算血流速度，與使用原信號(hào)直接計(jì)算血流速度的STFT法進(jìn)行對比，EEMD_N法能夠有效改善靠近血管壁低速血流的測量精度；最后，將EEMD_N法應(yīng)用于人體頸動(dòng)脈臨床試驗(yàn)，進(jìn)一步驗(yàn)證了EEMD_N法的有效性.

1 方法原理

1.1 多普勒血流計(jì)算血紅細(xì)胞隨血液流動(dòng)而移動(dòng)，使得血紅細(xì)胞與超聲波之間存在一定的相對運(yùn)動(dòng)，換能器接收到的回波信號(hào)與發(fā)射信號(hào)的頻率會(huì)出現(xiàn)偏差，這個(gè)偏差稱為多普勒頻移[23].超聲探頭接收到的回波信號(hào)攜帶了大量的多普勒頻移信息，根據(jù)克里斯琴·約翰·多普勒提出的多普勒效應(yīng)，可通過回波中攜帶的多普勒頻移測得血流速度[24].若血流速度為v，超聲波速為c，超聲波的發(fā)射頻率為f0，超聲探頭與血液流動(dòng)方向的夾角為θ(30°≤θ≤60°），根據(jù)多普勒頻移表達(dá)式可以得到超聲探頭接收到的回波頻率為：[25]

多普勒頻移與血流速度的關(guān)系為：

根據(jù)多普勒頻移fs，可得血流速度計(jì)算公式為：

1.2 EEMD分解血流多普勒信號(hào)是非平穩(wěn)信號(hào)，具有明顯的非平穩(wěn)性和隨機(jī)性.Wu等[20]將高斯白噪聲加入EMD法的原信號(hào)中，利用白噪聲的零均值特性，進(jìn)行多次平均后來相互抵消噪聲的影響，以達(dá)到抑制甚至完全消除高斯白噪聲的目的[26].基于這種思想提出的EEMD法本質(zhì)是一種多次疊加高斯白噪聲的EMD分解.

EEMD法的分解步驟如下[27]：

步驟1在原始信號(hào)x(t)中分別加入Z次均值為0、幅值標(biāo)準(zhǔn)差為常數(shù)的高斯白噪聲nz(t)，即：

步驟2對xz(t)分別進(jìn)行EMD分解，得到K個(gè)IMF分量和一個(gè)余項(xiàng)rz(t):

其中，C(z,k)(t)為第z次加入高斯白噪聲后分解得到的第k個(gè)IMF.

步驟3利用不相關(guān)隨機(jī)序列的統(tǒng)計(jì)均值為0的原理，將步驟2對應(yīng)的IMF進(jìn)行總體平均運(yùn)算，消除多次加入高斯白噪聲對真實(shí)IMF的影響，得到EEMD分解后的IMF及余項(xiàng)r(t)為：

其中，ck(t)為對原始信號(hào)進(jìn)行EEMD分解后得到的第k個(gè)IMF.最終得到K個(gè)IMF分量和一個(gè)余項(xiàng)r(t)的關(guān)系為：

1.3 EEMD_N法彩色多普勒血流成像系統(tǒng)采用收發(fā)一體的超聲探頭以固定頻率發(fā)射脈沖掃描血流感興趣區(qū)獲得回波.取同一探測位置的T次回波，可得到含有頻移信息的血流多普勒信號(hào).在本文提出的EEMD_N方法中，對血流多普勒信號(hào)進(jìn)行EEMD分解后，分別計(jì)算每個(gè)IMF的NFI，再根據(jù)NFI最優(yōu)閾值選擇有效血流成分計(jì)算血流速度.其中，NFI表示信號(hào)的波動(dòng)強(qiáng)度，可用于衡量信號(hào)的變化[28].由于血流速度呈拋物線型分布，血流中心區(qū)域信號(hào)的波動(dòng)通常會(huì)比血流邊緣區(qū)域信號(hào)的波動(dòng)要強(qiáng)烈.

對血流多普勒信號(hào)xi(t)進(jìn)行EEMD分解得到3個(gè)IMF分量和一個(gè)余項(xiàng).計(jì)算所有IMF1的NFI，與原信號(hào)的NFI進(jìn)行對比，結(jié)果如圖1所示，二者的NFI變化基本一致，并且血流中心區(qū)域信號(hào)的NFI一般要大于血流邊緣區(qū)域信號(hào)的NFI.血流速度在管腔中心最大，隨著遠(yuǎn)離管腔中心而靠近管壁時(shí)不斷減小，直至管壁附近趨近于零，表現(xiàn)出層流特性[29].血流中心區(qū)域IMF1的血流速度與原信號(hào)的結(jié)果基本一致，所以血流中心區(qū)域選擇IMF1計(jì)算血流速度.在血流邊緣區(qū)域，原信號(hào)和IMF1計(jì)算的血流速度都已偏離標(biāo)準(zhǔn)血流速度曲線.根據(jù)血流的層流特性，在邊緣區(qū)域，血流速度逐漸減小并趨近于零，IMF2和IMF3的血流速度在血流邊緣區(qū)域與標(biāo)準(zhǔn)血流速度曲線接近，并且IMF2的血流速度一般要大于IMF3的血流速度.EEMD_N法通過設(shè)置NFI閾值?1、?2，從IMF1、IMF2和IMF3中選擇有效IMF計(jì)算血流速度.

圖1 原信號(hào)與IMF 1的NFI對比Fig.1 The comparison between original signals and the NFI of IMF1.

EEMD_N法的具體步驟如下：

步驟1自上管壁至下管壁，提取I個(gè)探測位置的血流多普勒信號(hào)，作為EEMD分解的原信號(hào)

步驟2使用公式（4）～（8）對原始血流多普勒信號(hào)xi(t)進(jìn)行EEMD分解，獲得K個(gè)IMF分量和一個(gè)余項(xiàng)ri(t).

步驟3計(jì)算IMF的波動(dòng)指數(shù)

其中，T為IMF分量的時(shí)間采樣點(diǎn)數(shù)目，且T≥t≥1，b(t)指第t個(gè)時(shí)間采樣點(diǎn)的信號(hào)強(qiáng)度.對波動(dòng)指數(shù)進(jìn)行歸一化，得到NFI值：

其中，γmax和 γmin分別表示波動(dòng)指數(shù)的最大與最小值.

步驟4將IMF1的NFI與?1、 ?2進(jìn)行對比:

其中， φ(i)為判定結(jié)果，H代表待判信號(hào)為高頻分量，L1、L2代表待判信號(hào)為低頻分量.

步驟5由于血流是低回聲區(qū)域，血流信號(hào)幅值較低，易受噪聲干擾，使用NFI閾值會(huì)存在判定異常的情況.如圖1中的P1、P2點(diǎn)，使用式（11）后會(huì)將P1、P2誤判.為了避免誤判，定義誤判區(qū)域S0，以待判定信號(hào)的空間位置為S0中心，計(jì)算S0內(nèi)所有血流多普勒信號(hào)的NFI的均值為:

其中，m是S0內(nèi)血流多普勒信號(hào)的數(shù)量.

步驟6將均值與?1、 ?2使用式（11）進(jìn)行對比，式(11)可進(jìn)一步寫為:

步驟7可得EEMD_N法血流多普勒信號(hào)fimfs(i)的表達(dá)式為：

步驟8將fimfs(i)帶入公式（3），得到EEMD_N法的血流速度表達(dá)式：

圖2為EEMD_N法的流程圖.

圖2 EEMD_N法的流程圖Fig.2 Flow chart of the EEMD_N method.

2 實(shí)驗(yàn)與結(jié)果

基于Field II超聲仿真軟件建立頸動(dòng)脈仿真模型，使用EEMD_N法測量該模型的血流速度剖面，并與傳統(tǒng)STFT法的測量結(jié)果進(jìn)行比較，以驗(yàn)證EEMD_N法的有效性.為進(jìn)一步驗(yàn)證該方法的臨床適用性，將EEMD_N法應(yīng)用于人體頸總動(dòng)脈的臨床試驗(yàn)，并分析臨床試驗(yàn)與仿真實(shí)驗(yàn)的一致性.

2.1 仿真實(shí)驗(yàn)用Field II超聲仿真軟件搭建頸動(dòng)脈血流模型，其幾何參數(shù)見表1.考慮到人體血流速度隨心臟搏動(dòng)而改變，該模型管腔中心線上的血流速度(即最大血流速度vmax)隨時(shí)間變化，變化過程如圖3所示.在徑向方向上，管腔中心線處血流速度最大，隨著靠近血管壁，速度不斷減小，直至管壁附近趨近于0.整體血流速度剖面自上管壁至下管壁呈拋物線型分布：

其中，r為血流徑向距離，R為血管幾何半徑.

掃描頸動(dòng)脈血流模型獲取仿真血流多普勒信號(hào)，多普勒聲學(xué)和換能器的參數(shù)見表1.在圖3所示的脈動(dòng)血流速度曲線中選取4個(gè)典型時(shí)刻測量血流速度剖面，其中F1、F3處的血流速度變化最快，非平穩(wěn)性最為顯著；F2是收縮期峰值，血流速度達(dá)到最大；F4為舒張期，血流速度已明顯下降.

圖3 仿真頸動(dòng)脈中心血流速度曲線Fig.3 The blood flow velocities at the centerline of the carotid artery in simulations

表1 仿真與臨床試驗(yàn)參數(shù)Tab.1 Parameters in computer simulations and the clinical experiments

對于上述4個(gè)典型時(shí)刻，在血管腔內(nèi)I=100個(gè)探測位置上提取20個(gè)時(shí)間采樣長度的血流多普勒信號(hào).其中，圖4給出了F1時(shí)刻5個(gè)探測位置的血流多普勒信號(hào).當(dāng)徑向距離為0時(shí)，探測位置為管腔中心，此處血流速度最快且信號(hào)頻率最高；隨著血流遠(yuǎn)離中心靠近血管壁，血流多普勒信號(hào)的頻率會(huì)隨著血流速度的減小而降低；當(dāng)徑向距離為-5 mm和5 mm時(shí),靠近血管壁，血流速度趨近于0，信號(hào)頻率最低.另外，由于F1時(shí)刻血流速度變化最快，血流多普勒信號(hào)頻率變化顯著，以0處的結(jié)果為例，信號(hào)前半段較為平緩，后半段波動(dòng)明顯.

圖4 F1時(shí)刻5個(gè)徑向位置的血流多普勒信號(hào)Fig.4 The Doppler signals at five different radial positions at F1 time

本文中，EEMD算法中輔助噪聲強(qiáng)度設(shè)置為0.2，添加高斯白噪聲次數(shù)設(shè)置為300[30].對圖4中徑向距離0位置的血流多普勒信號(hào)進(jìn)行EEMD分解，結(jié)果如圖5所示，IMF根據(jù)頻率從高到低依次排列，IMF1包含原信號(hào)主要成分，由高頻信號(hào)構(gòu)成，IMF2、IMF3由低頻信號(hào)構(gòu)成.

圖5 F1時(shí)刻徑向距離0的血流多普勒信號(hào)EEMD分解結(jié)果Fig.5 The EEMD-based decomposed results for the Doppler signals at the radial position of 0 at F1 time

作為EEMD_N法的關(guān)鍵參數(shù)，誤判區(qū)域和NFI閾值的設(shè)置與EEMD_N的性能密切相關(guān).為確定最佳誤判區(qū)域與NFI閾值，使用EEMD_N法基于不同的誤判區(qū)域和NFI閾值測量F1～F44個(gè)時(shí)刻的血流速度剖面，并計(jì)算測量結(jié)果與標(biāo)準(zhǔn)血流速度剖面的歸一化均方根誤差(Normalized Root Mean Squared Error, NRMSE).如圖6所示，當(dāng)誤判區(qū)域比為8.4%時(shí)為最優(yōu)誤判區(qū)域，此時(shí)NRMSE最小.將誤判區(qū)域比設(shè)置為8.4%，使用不同的?1、?2測量血流速度剖面，其與標(biāo)準(zhǔn)血流速度剖面的NRMSE如圖7所示.使用傅里葉函數(shù)分別擬合?1、?2與平均NRMSE值如圖8所示，當(dāng) ?1=0.175、?2=0.11時(shí)為最優(yōu)NFI閾值，平均NRMSE值最小.使用最佳的誤判區(qū)域比和NFI閾值，基于EEMD_N法測量的速度剖面如圖9所示，此時(shí)測量的血流速度剖面最為準(zhǔn)確，EEMD_N法性能最好；傳統(tǒng)STFT法的平均NRMSE值為0.926 8，EEMD_N法的平均NRMSE為0.843 0；與傳統(tǒng)STFT法相比，誤差減小9.04%.仿真實(shí)驗(yàn)結(jié)果表明EEMD_N法能有效改善血流邊緣區(qū)域的血流速度分布，提高靠近血管壁的低速血流檢測精度.

圖6 不同誤判區(qū)域下的平均NRMSEFig.6 The average NRMSE with the different rate of the judged regions

圖7 誤判區(qū)域比為8.4%時(shí)，不同閾值下的平均NRMSEFig.7 The average NRMSEs for the different thresholds with the judged region rate of 8.4%.

圖8 基于傅里葉函數(shù)的平均NRMSE擬合曲線Fig.8 The Fourier function-based fitting curves of the average NRMSE

圖9 EEMD_N法與傳統(tǒng)STFT在F1～F4時(shí)刻的血流速度分布結(jié)果對比Fig.9 The comparison of blood flow velocity profiles between the EEMD_N and the STFT methods from F1 to F4

2.2 臨床試驗(yàn)為了驗(yàn)證EEMD_N法在臨床上的適用性，采集人體頸總動(dòng)脈的血流多普勒信號(hào)（數(shù)據(jù)來源于10名受試者，其中5名男性、5名女性，年齡：20～40歲，體重：45～80 kg，試驗(yàn)前所有受試對象都已被告知試驗(yàn)相關(guān)信息并自愿簽署知情同意書）.為避免血管壁信號(hào)的干擾，采用八階巴特沃斯高通濾波器對血流多普勒信號(hào)濾波，歸一化截止頻率為0.18[31].對濾波后的血流多普勒信號(hào)進(jìn)行EEMD分解，計(jì)算IMF1的NFI如圖10所示，臨床血流多普勒信號(hào)的NFI在血流中心區(qū)域較大，血流邊緣區(qū)域NFI較小，與仿真信號(hào)的NFI分布一致.基于仿真實(shí)驗(yàn)測定的最優(yōu)誤判區(qū)域比和NFI閾值，使用多普勒信號(hào)進(jìn)行EEMD分解，計(jì)算IMF1的NFI血流速度剖面.擬合STFT和EEMD_N法測量的血流速度剖面，結(jié)果如圖11所示.在10位受試者的試驗(yàn)結(jié)果中，基于STFT和EEMD_N的擬合曲線的決定系數(shù)R2的均值分別為0.797 2±0.021 6、0.913 4±0.015 9，說明EEMD_N的速度剖面更接近拋物線，受噪聲干擾更小.臨床試驗(yàn)與仿真實(shí)驗(yàn)得到的結(jié)果一致，與傳統(tǒng)STFT法對比，EEMD_N法能夠有效改善血流速度剖面的測量精度，特別是靠近血管壁的低速血流.

圖10 人體頸動(dòng)脈超聲血流多普勒信號(hào)進(jìn)行EEMD分解后IMF1的NFIFig.10 The NFIs of IMF1 decomposed from the Doppler flow signals on the human carotid artery by the EEMD method

圖11 基于STFT與EEMD_N法的人體頸動(dòng)脈血流速度分布對比Fig.11 The comparison for the blood flow velocity profiles of the human carotid artery measured by the STFT and the EEMD_N methods

3 結(jié)論

EEMD_N法基于EEMD法進(jìn)行血流速度分布的計(jì)算，能夠克服傳統(tǒng)方法在應(yīng)對非平穩(wěn)血流信號(hào)的局限性，解決傳統(tǒng)時(shí)頻分析法在估計(jì)血流速度分布時(shí)血流邊緣區(qū)域存在精確度不高的問題.通過計(jì)算NFI，設(shè)置NFI閾值和誤判區(qū)域來劃分血流中心和邊緣區(qū)域，使用NFI閾值判斷結(jié)果選擇有效血流成分，在仿真實(shí)驗(yàn)中與使用原信號(hào)計(jì)算血流速度剖面的傳統(tǒng)STFT法相比，EEMD法測量的血流速度剖面NRMSE減小9.04%.人體頸動(dòng)脈臨床試驗(yàn)中，EEMD_N法同樣具有可行性.由仿真實(shí)驗(yàn)和臨床試驗(yàn)的結(jié)果可知，EEMD_N法能夠?qū)拷鼙诘牡退傺鞲泳珳?zhǔn)的測量，能達(dá)到提取精確血流信號(hào)的目的.但EEMD_N法仍待改進(jìn)，首先，為消除白噪聲的干擾，需要Z足夠大，這樣分解信號(hào)的時(shí)間變長，方法的實(shí)時(shí)性會(huì)大大受到限制.其次，EEMD_N法對于血流中心區(qū)域的血流速度較傳統(tǒng)方法沒有明顯改善.如何解決這兩個(gè)問題是下一步研究的重點(diǎn).