CFI中基于動態(tài)區(qū)域劃分的非平穩(wěn)雜波抑制方法*

熊秀娟，肖 磊，陳 波，彭 勇，王耀彬

(西南科技大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院，四川 綿陽 621010)

0 引言

在臨床醫(yī)學(xué)領(lǐng)域中，CFI能夠?qū)崟r準(zhǔn)確無創(chuàng)地探測血流速度，是診斷心血管疾病的重要技術(shù)[1]。在CFI中，探頭接收到的回波信號包括血流信號、雜波信號（由血管搏動和組織慢速移動引起），通常雜波信號強(qiáng)度比血流信號強(qiáng)度高出 40～80 dB[2]，這給準(zhǔn)確地估計血流速度帶來了極大困難。因此，為了得到真實可靠的血流速度，必須對回波信號中的雜波進(jìn)行充分地抑制。

目前常用的雜波抑制器有：低階FIR濾波器[3]、基于投影初始化的IIR濾波器[4]和回歸濾波器[5]。其中低階FIR濾波器、投影初始化IIR濾波器和回歸濾波器均屬于靜態(tài)濾波器，它們的通帶截止頻率、阻帶衰減等特性較為固定。當(dāng)雜波信號是平穩(wěn)信號時，這些濾波器能夠獲得較為理想的效果。但在實際臨床診斷中，由于人體呼吸、脈搏等因素造成組織加速運動，導(dǎo)致雜波信號屬于非平穩(wěn)隨機(jī)信號，靜態(tài)濾波器對此類信號的處理效果不理想。后來有學(xué)者提出非平穩(wěn)雜波抑制法[1]，該方法對雜波抑制效果較好，但是在雜波信號較弱的區(qū)域會造成誤消除，導(dǎo)致血管內(nèi)血流信息的損失，對血流速度的估計造成極大誤差。

本文在非平穩(wěn)雜波抑制法的基礎(chǔ)上提出一種基于動態(tài)區(qū)域劃分的非平穩(wěn)雜波抑制方法。該方法首先利用能量特性對回波信號進(jìn)行動態(tài)區(qū)域劃分，然后結(jié)合非平穩(wěn)雜波抑制法和多項式回歸法對不同區(qū)域的信號分別進(jìn)行處理。實驗結(jié)果表明，該方法的算法復(fù)雜度低，滿足超聲診斷設(shè)備的實時性要求；同時，該算法能較好地抑制雜波并保證流速度剖面的完整性，是一種快速有效的雜波抑制方法。

1 原理與方法

1.1 動態(tài)區(qū)域劃分和局部平滑

多普勒回波信號是在同一探測位置發(fā)射K次脈沖波得到的一維向量矩陣。將軸向采樣容積為M的K個復(fù)解調(diào)信號排成二維矩陣，則對回波信號處理的過程即是對二維矩陣進(jìn)行運算。其中，信號構(gòu)成如圖1所示。通常，信號矩陣的行向量被稱為慢時信號，列向量被稱為快時信號[2]。

圖1 信號構(gòu)成圖

復(fù)解調(diào)回波信號可表示為矩陣X：

其中K是慢時信號方向的采樣容積(脈沖重復(fù)數(shù))，M是快時信號方向的采樣容積。

回波信號慢時信號方向的平均能量E_X表示為：

同理，快時信號方向上不同深度的慢時信號能量E_Xs(m)可由式(3)表示：

其中j是慢時信號方向上的第j個采樣點，m是快時信號方向上的第m個采樣點。由于不同區(qū)域的組織產(chǎn)生的回波信號強(qiáng)度不同，以式(4)的方式對信號進(jìn)行劃分，可以區(qū)分不同組織區(qū)的信號。

通常在動態(tài)組織區(qū)，雜波信號的強(qiáng)度比血流信號強(qiáng)度高出 40～80 dB，因此可以由式(5)再次對動態(tài)組織區(qū)的信號進(jìn)行劃分：

其中E_M是動態(tài)組織區(qū)快時信號上慢時方向的平均能量。

靜態(tài)組織區(qū)遠(yuǎn)離血管，整體趨于靜止，所以回波信號趨于平穩(wěn)，且不攜帶任何速度信息；血流區(qū)分布在血管內(nèi)部，加之血管中的血流速度呈拋物線分布且相對較穩(wěn)定，因此血流區(qū)的信號相對平穩(wěn)；雜波區(qū)受呼吸、脈搏等影響，造成血管及血管附近組織具有加速度，因此雜波區(qū)信號屬于非平穩(wěn)信號。

為保證雜波濾波器自適應(yīng)于組織運動且不造成血流信息的損失，本文算法采用空間平滑法處理動態(tài)組織區(qū)，如式(6)：

其中N是快時信號方向上參與空間平滑的區(qū)間半徑。

1.2 雜波頻率估計[1]

復(fù)解調(diào)回波信號x(m，k)用復(fù)指數(shù)形式表示為：

x(m，k)=A(m，k)ejφ(m，k)(7)其中 A(m，k)是信號的幅值，φ(m，k)是信號的相位。通常在多普勒回波信號中，雜波信號的相位可以用一個低階多項式表示：

由此可見，對雜波信號的瞬時頻率進(jìn)行估計，就是求出多項式系數(shù)ad的過程。

引入矩陣R和向量φ，對擬合多項式(9)求解得多項式系數(shù)向量a：

又由式(9)可得雜波瞬時頻率矩陣W=[w(m，1)w(m，2)… w(m，K)]T的表達(dá)式：

將式(10)得到的系數(shù)向量a帶入式(11)就可以進(jìn)一步求得雜波的瞬時頻率矩陣W。對于D的取值，參考文獻(xiàn)[5]通過分析低階多項式擬合的原理并進(jìn)行大量的實驗后得出：當(dāng)D≤4時，多項式能夠有效地擬合低頻信號。在超聲彩色血流成像中，當(dāng)D=3時，擬合多項式能夠精確地描述雜波的瞬時頻率。

1.3 局部非平穩(wěn)雜波抑制

多普勒回波信號經(jīng)過動態(tài)區(qū)域劃分和雜波瞬時頻率估計后就可進(jìn)行混頻處理。其中混頻是為了將非平穩(wěn)雜波信號的頻率搬移到零頻附近。為保證血流信號點的完整性，本文算法只對雜波區(qū)信號進(jìn)行混頻處理。利用前文得到的雜波瞬時頻率，混頻信號可表示為：

混頻過程可用式(14)表示：

通過混頻處理，雜波信號的非平穩(wěn)性大大降低，此時通過常用的靜態(tài)濾波器就可以對雜波信號進(jìn)行充分抑制。由于血流區(qū)信號可能含有較弱的雜波，因此本文算法選用能夠避免數(shù)據(jù)點損失的多項式回歸濾波器對雜波信號進(jìn)行抑制，這樣既保證了雜波信號的有效去除，又保證了血流信息的完整性。算法流程如圖2所示。

圖2 算法流程圖

2 仿真實驗

2.1 仿真模型與參數(shù)

為了分析和驗證所提出的雜波抑制方法的有效性，本文在計算機(jī)上用MATLAB R2012a軟件平臺進(jìn)行仿真，運行環(huán)境為 Centrino2（1.66 GHz），采用丹麥技術(shù)大學(xué)學(xué)者Jenson等研制的Field ii軟件包仿真超聲多普勒回波信號。仿真所用的超聲血流模型如圖3所示，仿真參數(shù)如表1所述。圖4是通過仿真模型得到的一條回波信號波形圖。

2.2 仿真實驗結(jié)果

為驗證本文算法對雜波的抑制效果，將其與投影初始化濾波器、多項式回歸濾波器、SVD濾波器和非平穩(wěn)濾波器進(jìn)行比較。圖5所示的是一組（16條）回波信號經(jīng)過各濾波器后自相關(guān)估計得到的相應(yīng)的血流速度剖面波形圖，圖6所示的是218組相鄰回波信號通過各濾波器后自相關(guān)估計、編碼映射得到的血流速度圖。兩圖直觀地展示了不同濾波器對雜波的抑制性能。

圖3 超聲血流模型圖

表1 Field II仿真參數(shù)設(shè)置

圖6 編碼映射血流圖

由圖5(e)可以看出，采用本文算法進(jìn)行雜波抑制后，自相關(guān)估計得出的血流速度信息幾乎僅存在于血管內(nèi)部（約在掃描深度24 mm～34 mm之間，血管理論直徑為8 mm），且速度分布相對比較完整，這與實際情況基本相符。由圖5(a)～5(d)可以看出，經(jīng)過其他濾波器處理后的回波信號經(jīng)自相關(guān)估計后得到的血流速度信息不僅存在于血管內(nèi)部，也存在于組織區(qū)，這顯然與實際情況不符。

從圖6(e)可以看出，本文算法得到的血管壁清晰、完整，血管截面速度分布均勻；從圖 6(b)和 6(d)可以看出，經(jīng)過回歸濾波器和非平穩(wěn)濾波器得到的血流速度分布不均勻，且組織區(qū)出現(xiàn)偽血流信息；從圖6(c)中可以看出，經(jīng)過SVD濾波器后得到的血流速度分布不完整，血管內(nèi)速度出現(xiàn)零點，且血管直徑明顯小于理論值，成像效果較差；從圖6(a)可以看出，投影初始化IIR濾波器由于采用固定頻率作為截止頻率而造成數(shù)據(jù)點損失較大，使得血管直徑明顯小于理論值（理論血管直徑為8 mm）。

此外，為了從客觀上對本文算法進(jìn)行評價，文章從最大速度估計、雜波血流比和算法運行時間三方面對算法進(jìn)行比較。比較結(jié)果如表2所示。

表2 算法比較表

從表2中可以看出，本文算法估計的最大速度值接近理論最大速度值0.972 m/s(本文算法的最大估計速度與多項式回歸濾波器相同，這是因為該算法采用多項式回歸法對血管中心區(qū)域信號進(jìn)行處理)，表明采用本文算法對雜波信號進(jìn)行抑制時，原始信號中的血流信號損失較?。唤?jīng)本文算法濾波后，信號的雜波血流比明顯低于其他幾種濾波算法，這表明該算法能夠有效地對原始信號中所含的雜波信號進(jìn)行抑制；在運行時間上，本文算法的運行時間遠(yuǎn)小于其他幾種算法，這是由于該算法對靜態(tài)組織區(qū)的信號不做處理導(dǎo)致的。綜上所述，本文算法在雜波抑制效果、血流信息的保留和算法時間復(fù)雜度上都能取得較為滿意的結(jié)果，由此證明此算法較其他雜波抑制算法更有效。

3 結(jié)論

本文提出的基于動態(tài)區(qū)域劃分的非平穩(wěn)雜波抑制方法一改以往采用單一濾波方法進(jìn)行雜波抑制的模式，通過分析回波信號的特點，首先采用動態(tài)區(qū)域劃分法將信號分割為不同部分，再根據(jù)各部分信號的特點做出相應(yīng)處理。本文算法采取對靜態(tài)組織區(qū)信號完全保留的方法，既保證濾波器不會在此區(qū)域引入噪聲，又保證了運行速度。同時，對雜波區(qū)和血流區(qū)采取不同的濾波方法，既保證了雜波的有效去除，又避免了血流信號的損失。此外，該算法采用矩陣進(jìn)行推理運算，易于工程實現(xiàn)。

另外，雜波瞬時頻率的估計是基于低階多項式來實現(xiàn)的，當(dāng)多項式階數(shù)取值適當(dāng)時，能較好地描述雜波頻率，但當(dāng)階數(shù)取值欠佳時，不能很好地對頻率進(jìn)行描述。因此，采取一種能精確描述雜波瞬時頻率的方法將是下一步研究的目標(biāo)。

[1]王沛東，沈毅，王艷.超聲彩色血流成像中非平穩(wěn)雜波的抑制[J].中國生物醫(yī)學(xué)工程學(xué)報，2008，27(2)：240-243.

[2]尤偉，汪源源.超聲彩色血流成像中基于動態(tài)區(qū)域劃分抑制雜波的方法[J].儀器儀表學(xué)報，2010，31(3)：594-599.

[3]JENSEN JA.Estimation of blood velocities using ultrasound[C].Cambridge，UK：Cambridge University Press，1999：195-226.

[4]CHORNOBOY E S.Initialization for improved IIR filter performance[J].IEEE Trans.Signal Processing，1992，40(5)：543-550.

[5]PARK G，KIM Y，SHIM H，et al.New adaptive clutter rejection based on spectral decomposition and tissue acceleration for ultrasound color Doppler imaging[C].IEEE Joint UFFC，EFTF and PFM Symposium，2013：1484-1487.