杜宜綱 張夢(mèng)一 陳思平 李 勇

1(深圳大學(xué)醫(yī)學(xué)部生物醫(yī)學(xué)工程學(xué)院，醫(yī)學(xué)超聲關(guān)鍵技術(shù)國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，廣東 深圳 518060)2(深圳邁瑞生物醫(yī)療電子股份有限公司，廣東 深圳 518057)3(香港中文大學(xué)計(jì)算機(jī)科學(xué)與工程系，香港)

基于低秩理論的超聲血流壁濾波器

杜宜綱1，2張夢(mèng)一3陳思平1*李 勇2

1(深圳大學(xué)醫(yī)學(xué)部生物醫(yī)學(xué)工程學(xué)院，醫(yī)學(xué)超聲關(guān)鍵技術(shù)國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，廣東 深圳 518060)2(深圳邁瑞生物醫(yī)療電子股份有限公司，廣東 深圳 518057)3(香港中文大學(xué)計(jì)算機(jī)科學(xué)與工程系，香港)

傳統(tǒng)的超聲血流壁濾波器采用單一的截止頻率，對(duì)于心跳、呼吸等因素所產(chǎn)生的不同組織運(yùn)動(dòng)濾波效果不佳。提出一種基于低秩理論的超聲血流壁濾波方法。通過(guò)分析超聲血流的特點(diǎn)，將其公式化并建立低秩模型，模型包括一個(gè)矩陣秩的最小化和一個(gè)稀疏矩陣問(wèn)題。松弛后轉(zhuǎn)化為最小化核范數(shù)與1范數(shù)的線性組合，成為凸優(yōu)化可解問(wèn)題，從而實(shí)現(xiàn)低秩濾波。這種濾波的創(chuàng)新性在于它是自適應(yīng)的，同時(shí)考慮了信號(hào)的低秩性與稀疏性?xún)蓚€(gè)特點(diǎn)，通過(guò)線性組合使得濾波效果達(dá)到最優(yōu)。通過(guò)低秩濾波，超聲血流仿真數(shù)據(jù)得到更加精確的血流信號(hào)。對(duì)比傳統(tǒng)的FIR濾波，分別采用12階、56階和92階的FIR濾波器，得到血流信號(hào)對(duì)比實(shí)際仿真數(shù)據(jù)的RMS誤差分別約為34%、16%和12%，而低秩濾波的RMS誤差則低于0.001%。低秩濾波相比傳統(tǒng)FIR濾波，不僅提高了精度，而且濾波后的信號(hào)長(zhǎng)度不會(huì)損失。但是低秩模型計(jì)算過(guò)于復(fù)雜，因此目前還很難應(yīng)用在實(shí)時(shí)的超聲系統(tǒng)中。

超聲血流成像；壁濾波；低秩模型；FIR濾波器

引言

超聲血流成像中最關(guān)鍵的一步是提取血流信號(hào)，因?yàn)檠旱姆瓷湎禂?shù)大約只有組織的1/500～1/10[1-4]。在超聲回波中，組織和血流信號(hào)混在一起，只有提取并得到信噪比高的血流信號(hào)，才能計(jì)算出高靈敏度和高精度的血流速度。這兩種信號(hào)在理想狀態(tài)下的區(qū)別在于組織信號(hào)不隨時(shí)間變化，而血流信號(hào)相反。傳統(tǒng)方法常采用高通濾波器[5]，例如FIR和IIR濾波器，通過(guò)設(shè)置一個(gè)截止頻率將變化速率慢的組織信號(hào)濾掉，剩下的是變化速率快的血流信號(hào)。

然而，現(xiàn)實(shí)并非如此理想。在超聲成像中，血流成像和常規(guī)二維灰階成像的發(fā)射波形是不同的。血流成像注重靈敏度，發(fā)射波形要長(zhǎng)一些，從而增加信噪比。灰階成像更重視分辨率，因此發(fā)射波形要短一些，從而可以直接提升縱向深度的分辨率。由于這兩種成像發(fā)射波形不同，對(duì)于傳統(tǒng)彩超顯示方式(血流與二維灰階圖同時(shí)顯示)，血流的發(fā)射不能是連續(xù)的，中間需要留出時(shí)間發(fā)射接收用來(lái)生成灰階圖像的信號(hào)。此外，傳統(tǒng)的彩超采用的是聚焦波發(fā)射逐線掃描方式。每根血流的掃描線發(fā)射次數(shù)不能太多，否則會(huì)降低圖像的幀率。一般來(lái)說(shuō)，血流的一條成像線大概需要8～16次發(fā)射[5-6]。在同一個(gè)空間成像位置上，將多次發(fā)射的回波信號(hào)沿時(shí)間方向做濾波后，得到血流信號(hào)，然后可以采用傳統(tǒng)的自相關(guān)法[7]計(jì)算出血流的速度。這里每組離散信號(hào)只有8～16個(gè)，而FIR濾波器需要很高的階數(shù)才能滿足設(shè)計(jì)需要，因此很難用于血流信號(hào)的提取。假設(shè)血流信號(hào)長(zhǎng)度為10，F(xiàn)IR階數(shù)為6，這樣FIR濾波器系數(shù)個(gè)數(shù)為7，原始信號(hào)通過(guò)高通濾波器之后的長(zhǎng)度只有4。濾波后的信號(hào)長(zhǎng)度越短，采用自相關(guān)法計(jì)算出來(lái)的平均速度值穩(wěn)定性(精度)越差。只有6階的FIR濾波器的頻率響應(yīng)如圖1所示，其效果顯然難以達(dá)到相應(yīng)的精度要求。IIR濾波器則不需要太高的階數(shù)，通常采用2～3階的濾波器。Bjrum等對(duì)比了幾種不同的IIR濾波器[5]。從綜合性能上看，投影初始化的IIR高通濾波器表現(xiàn)最佳。這也是當(dāng)前超聲血流成像中最常采用的壁濾波器。采用這種濾波器，達(dá)到類(lèi)似的效果，不但其階數(shù)幾乎可以比FIR低10倍，而且濾波后的信號(hào)長(zhǎng)度也不會(huì)損失，因此還可以提高速度計(jì)算的穩(wěn)定性。

圖1 不同階數(shù)FIR濾波器的頻率響應(yīng)Fig.1 Frequency response of different orders’ FIR filters

圖2 聚焦波與平面波發(fā)射聲場(chǎng)的比較結(jié)果。聚焦波的發(fā)射焦點(diǎn)為60 mm。采用Field II[9-10]仿真發(fā)射聲場(chǎng)，64個(gè)陣元，探頭中心頻為率5 MHz，陣元間距為0.22 mmFig.2 Compared results of the focused wave and plane wave emitting fields. The focal depth is 60 mm. Simulations are made by Field II[9-10]. The aperture has 64 elements, the center frequency is 5 MHz, and the pitch is 0.22 mm

超聲血流成像的幀率主要受限于上述這種彩超和灰階成像聚焦波交替掃描的方式，而每條血流成像線的連續(xù)掃描次數(shù)也受到幀率限制不能太高，因此采用投影初始化IIR濾波器是這種掃描方式的一個(gè)合適的選擇。近年來(lái)，由于平面波超快成像的發(fā)展，成像的掃描方式從逐線發(fā)展到部分區(qū)域乃至整體一次掃描成像，這使得血流的顯示幀率有了很大的提高[8]，因此血流成像時(shí)同一個(gè)位置掃描次數(shù)也不再受到幀率的限制，這樣高階的FIR濾波器也可以用在血流信號(hào)的提取上。但是，由于采用的是平面波掃描，沒(méi)有發(fā)射聚焦，因此發(fā)射聲場(chǎng)的旁瓣會(huì)很高。如圖2所示，采用了Field II[9-10]仿真發(fā)射聲場(chǎng)，即使不在焦點(diǎn)深度的聚焦波聲場(chǎng)旁瓣也比平面波低很多，加之血流信號(hào)的幅值本身就比組織信號(hào)低10～27dB[1]，所以血管內(nèi)的回波很容易受到血管壁以及組織的回波干擾產(chǎn)生雜波信號(hào)。如圖3所示，采用了Field II函數(shù)calc_scat_multi[11]生成通道數(shù)據(jù)，平面波的雜波信號(hào)要比聚焦波的大很多。通過(guò)仿真計(jì)算，得到聚焦波發(fā)射時(shí)血管里的雜波平均要比平面波的小16.2 dB，相當(dāng)于平面波的雜波信號(hào)強(qiáng)度是聚焦波的6～7倍。如果無(wú)法抑制這些雜波信號(hào)，那么平面波成像雖然可以提高幀率，但是血流的靈敏度和精度比起傳統(tǒng)方法都會(huì)有較大損失。

圖3 采用Field II仿真平面波和聚焦波成像。黑色部分為模擬血管無(wú)反射物，血管壁為強(qiáng)反射點(diǎn)，采用Field II函數(shù)calc_scat_multi[11]生成通道數(shù)據(jù)，然后分別做平面波波束合成和動(dòng)態(tài)接收聚焦波束合成。平面波成像時(shí)血管內(nèi)平均雜波信號(hào)為-106.9 dB，動(dòng)態(tài)接收聚焦成像后血管內(nèi)平均雜波信號(hào)為-123.1dB(這里只模擬了二維血管，因此雜波信號(hào)很低。血管形狀在Field II仿真中為長(zhǎng)方體)。(a)平面波；(b)聚焦波Fig.3 Plane wave and focused wave imaging generated by Field II, where the blood vessel is simulated without scatterers and the vessel wall has the stronger scattering points. The channel data are generated by Field II function “calc_scat_multi”[11]. The mean clutter signals are around -106.9 dB for plane wave imaging and around -123.1 dB for dynamic receive focused imaging (Only 2D vessel is simulated so that the clutter signals are very low. The vessel shape is cubic in Field II simulation). (a)Plane wave；(b) Focused wave

平面波血流成像通常還是延用了傳統(tǒng)的FIR或者IIR壁濾波器。為增加信噪比，有些研究采用了編碼發(fā)射平面波技術(shù)[8,12]，還有的研究采用了造影劑進(jìn)一步提高平面波的成像效果[13]。本研究希望設(shè)計(jì)出一種新型的超聲血流壁濾波器，從而直接克服上述由于血流信號(hào)弱、信噪比低、雜波干擾大造成的血流精度和靈敏度的損失。這種濾波器基于一個(gè)低秩模型，是針對(duì)平面波發(fā)射而設(shè)計(jì)的。這種低秩模型可以用來(lái)區(qū)分?jǐn)?shù)據(jù)中數(shù)值較小的變化量與數(shù)值較大的常量，它們分別對(duì)應(yīng)了能量較弱的血流信號(hào)和相對(duì)變化很小且能量較強(qiáng)的雜波信號(hào)。因?yàn)檫@些雜波信號(hào)大多來(lái)自周?chē)\(yùn)動(dòng)較小的組織或者血管壁的強(qiáng)反射，所以雜波信號(hào)隨時(shí)間變化較小。雖然血流的信號(hào)很弱，但它隨時(shí)間變化較大。可見(jiàn)，血流與雜波的這些特點(diǎn)恰好符合了低秩模型理論。

從數(shù)學(xué)本質(zhì)來(lái)說(shuō)，低秩濾波實(shí)際上是將一個(gè)矩陣分解為一個(gè)低秩矩陣和一個(gè)稀疏矩陣之和，因此也被稱(chēng)為矩陣低秩稀疏分解或者魯棒主成分分析(RPCA)[14]。從理論和算法的層面上看，低秩濾波問(wèn)題和壓縮感知、矩陣秩最小化一脈相承。壓縮感知是對(duì)傳統(tǒng)奈奎斯特采樣定理的重大突破，是信號(hào)處理和計(jì)算機(jī)視覺(jué)領(lǐng)域的研究熱點(diǎn)[15-16]。當(dāng)采樣信號(hào)具有稀疏性時(shí)，壓縮感知能以低于奈奎斯特采樣率來(lái)進(jìn)行采樣，在采樣矩陣具有約束等距性(restricted isometry property)時(shí)，可以準(zhǔn)確重構(gòu)原信號(hào)。從數(shù)學(xué)本質(zhì)上來(lái)說(shuō)，信號(hào)的稀疏性由其對(duì)應(yīng)向量的0范數(shù)表示。在常見(jiàn)的壓縮感知算法里，0范數(shù)被其凸包絡(luò)(即1范數(shù))逼近，從而將信號(hào)恢復(fù)轉(zhuǎn)化成一個(gè)低復(fù)雜度的凸優(yōu)化問(wèn)題。矩陣秩最小化相當(dāng)于將壓縮感知從向量擴(kuò)展到矩陣。矩陣的秩相當(dāng)于奇異值向量的0范數(shù)，因此可以用其凸包絡(luò)1范數(shù)逼近，即矩陣的核范數(shù)(nuclear norm)[17-18]。矩陣秩最小化問(wèn)題的應(yīng)用包括協(xié)同濾波、推薦系統(tǒng)，即著名的Netflix問(wèn)題[19]。低秩濾波可以看成以上兩類(lèi)問(wèn)題的衍生和結(jié)合：在矩陣分解的過(guò)程中，同時(shí)考慮稀疏特征和低秩特征。只要低秩矩陣奇異值向量和稀疏矩陣的非零元素滿足一定條件，最小化核范數(shù)和1范數(shù)的線性組合能以很高的概率恢復(fù)原矩陣[20]。當(dāng)前，應(yīng)用在超聲血流成像中的低秩濾波還不多見(jiàn)。1997年Ledoux等首先提出了基于奇異值分解(SVD)的血流信號(hào)雜波濾波方法[21]，Tao等提出了基于PCA的血流壁濾波器[22]，Mauldin等提出了一個(gè)SVD濾波器，從而可以過(guò)濾掉圖像中的靜態(tài)偽影[23]。另外，Yu等提出了基于一維漢克(Hankel)矩陣SVD以及特征值的壁濾波器[24-26]。這些濾波器目前還都未能實(shí)時(shí)地應(yīng)用在商用機(jī)上，并且上述低秩濾波的研究并不是針對(duì)平面波成像而進(jìn)行的。本研究主要是針對(duì)平面波信號(hào)本身比較弱的特點(diǎn)，設(shè)計(jì)出更加魯棒的壁濾波器，從而實(shí)現(xiàn)平面波成像時(shí)更精確的血流信號(hào)提??；將低秩濾波信號(hào)與傳統(tǒng)的FIR濾波結(jié)果進(jìn)行比較，進(jìn)一步驗(yàn)證低秩模型的精確性。

1 方法

1.1 理論模型

b0,i=f0,i+c0,i

(1)

式中：f0,i為第i幅圖的血流信號(hào)(flow)，c0,i為第i幅圖的組織信號(hào)或者叫做雜波信號(hào)(clutter)；b0,i是第i次超聲成像后得到的圖像，顯示血管和其他組織結(jié)構(gòu)。

因此，在式(1)中，f0,i在非血管處的值是0，而c0,i在血管處的值要比在非血管處小一些，也就是混在血流中的雜波信號(hào)要比組織信號(hào)小一些。計(jì)算血流速度，理論上至少需要兩個(gè)不同時(shí)刻的圖像信息。設(shè)每幅圖像含有m個(gè)離散像素點(diǎn)，那么

(2)

式(1)可以被寫(xiě)成矩陣的形式

B0=F0+C0

(3)

式中：B0為一個(gè)m行n列的矩陣，每列數(shù)據(jù)表示一幅圖像，每幅圖像含有m個(gè)元素；F0為血流信號(hào)矩陣；C0為組織(雜波)信號(hào)矩陣。

上述模型的特點(diǎn)包括：組織(雜波)信號(hào)在時(shí)間上具有較強(qiáng)的相關(guān)性，血流信號(hào)在時(shí)間上是變化的，血流信號(hào)在空間上是稀疏的，血流信號(hào)在空間上具有連續(xù)性?；谶@些特點(diǎn)，可進(jìn)行進(jìn)一步分析，即利用組織信號(hào)在時(shí)間上的相關(guān)性，得到C0是一個(gè)低秩(low rank)矩陣；傳統(tǒng)的FIR或者IIR壁濾波器沿時(shí)間方向進(jìn)行濾波(逐點(diǎn))，而全新的低秩模型不僅利用信號(hào)在時(shí)間上的特點(diǎn)，還利用信號(hào)在空間上的特點(diǎn)；血流信號(hào)是稀疏的，因此‖F(xiàn)0‖0較小(‖F(xiàn)0‖0是F0的0范數(shù)，表示F0的非零個(gè)數(shù)之和)。將這些分析公式化，有

subject toC+F=B0

(4)

式中，C和F分別代表了計(jì)算出來(lái)的組織和血流信號(hào)(式(3)中的C0和F0則是實(shí)際的組織和血流信號(hào))。

由于矩陣秩的最小化和稀疏矩陣問(wèn)題是指數(shù)復(fù)雜度的(NP hard)[14]，需要尋找凸包絡(luò)(convex envelop)，從而將其轉(zhuǎn)化為一個(gè)凸優(yōu)化可解的問(wèn)題，如圖4所示。

圖4 函數(shù)g(x)是f(x)的凸包絡(luò)Fig.4 The function g(x) is the convex envelope of f(x)

根據(jù)相關(guān)參考文獻(xiàn)[15]，1范數(shù)可視為0范數(shù)的凸包絡(luò)，而核范數(shù)(nuclear norm)則是矩陣秩(rank)的凸包絡(luò)[17]。1范數(shù)是所有元素絕對(duì)值的和，核范數(shù)是矩陣所有奇異值的和。這樣，優(yōu)化問(wèn)題即式(4)的目標(biāo)函數(shù)松弛(relax)后，可得到

subject toC+F=B0

(5)

式中，‖C‖*表示C的核范數(shù)，‖F(xiàn)‖1表示F的1范數(shù)。

這樣，可將一個(gè)非凸優(yōu)化問(wèn)題(即式(4))轉(zhuǎn)化成一個(gè)凸優(yōu)化問(wèn)題(即式(5))。具體來(lái)說(shuō)，式(5)屬于凸優(yōu)化問(wèn)題中的半正定優(yōu)化問(wèn)題(semi-definite programming)[27]，這類(lèi)問(wèn)題可以通過(guò)一些通用軟件包來(lái)計(jì)算。例如，cvx、sedumi、mosek等屬于多項(xiàng)式復(fù)雜度(polynomial complexity)。然而，在矩陣維度較大的情況下，這些通用軟件包會(huì)遇到求解不穩(wěn)定的情況。在這種情況下，Lin等提出了專(zhuān)門(mén)針對(duì)這類(lèi)凸優(yōu)化問(wèn)題(即式(5))的軟件包[28-29]。相比通用軟件包，這類(lèi)專(zhuān)用軟件包更穩(wěn)定，收斂更快，并且可以將誤差控制在可接受范圍內(nèi)。

1.2 仿真方法

首先，模擬一組隨時(shí)間變化的超聲血流圖像。其中，組織為強(qiáng)反射，圖像值為高；血管中的雜波為較弱信號(hào)，其圖像值為中；血流信號(hào)的強(qiáng)度最低，其圖像值最小。在血管中，圖像值為雜波信號(hào)疊加血流信號(hào)F0。這組圖像中只有血流信號(hào)F0隨時(shí)間變化，其余保持常數(shù)，具體參數(shù)可參考表1。

表1 血流仿真參數(shù)Tab.1 Parameters for blood flow simulation

信號(hào)由隨機(jī)函數(shù)乘以相應(yīng)的反射系數(shù)產(chǎn)生，血管外的組織信號(hào)與血管內(nèi)的雜波信號(hào)統(tǒng)稱(chēng)為C0，將其與血流信號(hào)F0疊加得到用來(lái)做濾波實(shí)驗(yàn)的超聲血流圖像B0。本研究建立了兩種不同的血流仿真數(shù)據(jù)：第1種采用了固定的血流速度，在血管內(nèi)各個(gè)位置上的速度恒定不變；第2種模擬了隨位置變化的血流速度，即越靠近血管中心軸的速度越大，越靠近血管壁的速度越小。這兩種不同的血流仿真數(shù)據(jù)分別采用了兩種不同的掃描幀率。分別采用3種不同階數(shù)的FIR濾波器(具體階數(shù)請(qǐng)參考表1)和低秩濾波器進(jìn)行濾波，得到血流信號(hào)。由于隨時(shí)間變化的超聲圖像維度較大，采用Lin等提出的專(zhuān)門(mén)針對(duì)解這類(lèi)維度較大問(wèn)題的軟件包(Inexact ALM[29])進(jìn)行計(jì)算，得到低秩濾波后的血流信號(hào)[28-29]。濾波結(jié)果對(duì)比預(yù)先合成的血流信號(hào)F0，采用均方根誤差RMSE(root mean square error)，衡量幾種不同階數(shù)FIR濾波和低秩濾波的效果。RMS誤差的計(jì)算公式如下：

(6)

式中，F(xiàn)0是預(yù)先合成的血流信號(hào)，F(xiàn)是采用傳統(tǒng)的和本研究提出的方法(FIR和低秩)濾波后得到的血流信號(hào)。

2 結(jié)果

對(duì)前述兩種不同的血流仿真數(shù)據(jù)(血流速度1和2見(jiàn)表1)，分別進(jìn)行了3種不同階數(shù)的FIR濾波和低秩濾波。圖5為血流速度恒定(血流速度1)時(shí)的濾波結(jié)果，圖6為血流速度隨位置變化(血流速度2)時(shí)其中一個(gè)空間位置上的血流濾波結(jié)果，表2、3分別表示濾波結(jié)果與實(shí)際血流信號(hào)(預(yù)先合成的信號(hào))之間的RMS誤差。為了驗(yàn)證其穩(wěn)定性，每種方法分別計(jì)算了6次，每次采用隨機(jī)函數(shù)產(chǎn)生略有不同的超聲圖像以及血流信號(hào)，其反射系數(shù)始終保持不變。3種不同階數(shù)的FIR濾波在通過(guò)多次計(jì)算后，在不同仿真血流速度和掃描幀率下，RMS誤差保持在34%、16%和12%左右，浮動(dòng)較小。采用低秩濾波，經(jīng)過(guò)多次計(jì)算后，RMS誤差保持在0.000 7%左右。

圖5 血流速度恒定為0.25 m/s時(shí)FIR與低秩濾波的結(jié)果比較。F0為實(shí)際血流信號(hào)，其他為濾波后的血流信號(hào)(血流信號(hào)強(qiáng)度基于表1所述的血流反射系數(shù)，無(wú)物理單位)Fig.5 Fixed blood flow velocity is 0.25 m/s, comparisons of FIR and low rank filtered results. F0 is the true blood signal and the rest are filtered signals. (The blood signal is based on the blood reflection coefficient and no physical unit)

圖6 血流速度隨位置變化(變化范圍0.05～1 m/s)時(shí)某一特定位置上的FIR與低秩濾波的結(jié)果。F0為實(shí)際血流信號(hào)，其他為濾波后的血流信號(hào)(血流信號(hào)強(qiáng)度基于表1所述的血流反射系數(shù)，無(wú)物理單位)Fig.6 Blood flow velocity varied spatially in the range from 0.05 m/s to 1 m/s. The figure shows the results of FIR and low rank filtered results at a specific location. F0 is the true blood signal and the rest are filtered signals (The blood signal is based on the blood reflection coefficient and no physical unit)

表2 血流速度恒定為0.25 m/s時(shí)不同階數(shù)的FIR和低秩濾波結(jié)果的RMS誤差(%)

Tab.2 RMS errors (%) of the FIR and low rank filtered results for a fixed blood flow velocity of 0.25 m/s

測(cè)量組別濾波方法FIR-12FIR-56FIR-92低秩135.016.212.50.00073235.416.211.80.00073334.316.112.00.00072433.215.211.60.00068534.015.511.20.00076635.616.912.60.00074

表3 血流速度隨位置變化(變化范圍0.05～1 m/s)時(shí)不同階的FIR和低秩濾波結(jié)果的RMS誤差(%)

Tab.3 RMS errors (%) of the FIR and low rank filtered results when the blood flow velocity varied spatially from 0.05～1 m/s

測(cè)量組別濾波方法FIR-12FIR-56FIR-92低秩134.515.911.80.00071234.515.911.90.00068334.516.011.90.00070434.115.711.70.00072534.616.011.80.00069634.715.911.80.00072

3 討論

低秩濾波的主要?jiǎng)?chuàng)新性在于它是自適應(yīng)的，而傳統(tǒng)FIR不同的階數(shù)頻率響應(yīng)區(qū)別很大，如圖1所示。由于階數(shù)不夠，可能會(huì)導(dǎo)致一些低速的血流被濾掉。然而，即使FIR的階數(shù)足夠高，它也是建立在雜波頻率為零的基礎(chǔ)上的。在現(xiàn)實(shí)中，由于心跳、呼吸或者探頭顫抖等影響，導(dǎo)致雜波的頻率可能隨圖像中的位置發(fā)生改變。這種變化是非常隨機(jī)的，所以FIR采用單一的濾波截止頻率，無(wú)論多少階的濾波器都會(huì)造成顧此失彼的結(jié)果。而本研究采用低秩模型，考慮到了圖像的低秩性以及血流信號(hào)的稀疏性?xún)蓚€(gè)特點(diǎn)，通過(guò)矩陣秩的最小化和稀疏矩陣的線性組合，使濾波效果達(dá)到最優(yōu)。

通過(guò)比較FIR和低秩濾波結(jié)果，可以看出：對(duì)于不同的血流速度以及掃描頻率，低秩濾波的效果都明顯好于FIR的效果。無(wú)論FIR還是低秩濾波，它們相對(duì)于實(shí)際信號(hào)的RMS誤差在6次仿真中得到的結(jié)果都相對(duì)比較穩(wěn)定。階數(shù)越高，F(xiàn)IR的濾波效果越好，但是濾波后由于初始化損失的信號(hào)也就越多。例如，92階的FIR濾波后，將會(huì)損失92個(gè)時(shí)間上的信號(hào)。因此，雖然FIR可以通過(guò)增加階數(shù)來(lái)提高濾波效果，但是也會(huì)造成信號(hào)長(zhǎng)度的損失。對(duì)于采用求平均的方法計(jì)算血流速度，信號(hào)長(zhǎng)度的縮短會(huì)降低計(jì)算結(jié)果的穩(wěn)定性。

此外，仿真數(shù)據(jù)與現(xiàn)實(shí)中的超聲圖像還存在一定差距。首先，血管中除了存在血流信號(hào)與雜波信號(hào)外，還有一些系統(tǒng)噪聲，這些噪聲有時(shí)可能要大過(guò)雜波信號(hào)，并且這些噪聲不像上述的雜波信號(hào)那么穩(wěn)定，它們的隨機(jī)性更高，混在血流信號(hào)中難以濾去。正如前面所述，人體中的組織信號(hào)不是恒定不變的，血管壁和組織都會(huì)隨心臟跳動(dòng)產(chǎn)生位置上的改變，這樣一來(lái)，血管中的雜波信號(hào)也會(huì)隨之發(fā)生變化。未來(lái)還需要建立更加復(fù)雜的仿真數(shù)據(jù)，用來(lái)模擬系統(tǒng)噪聲以及波動(dòng)的雜波信號(hào)，解這種問(wèn)題的低秩模型也會(huì)更加復(fù)雜。

本研究提出的低秩模型是針對(duì)仿真數(shù)據(jù)的，而未來(lái)還需要針對(duì)實(shí)驗(yàn)數(shù)據(jù)做進(jìn)一步的改進(jìn)和優(yōu)化，最終的驗(yàn)證應(yīng)該基于超聲系統(tǒng)產(chǎn)生的實(shí)際數(shù)據(jù)。此外，低秩模型計(jì)算量很大，本研究提出的低秩模型需要的計(jì)算時(shí)間是FIR的3倍左右。如果增加系統(tǒng)噪聲并且考慮血管壁和組織運(yùn)動(dòng)產(chǎn)生的隨時(shí)間波動(dòng)的雜波信號(hào)，那么低秩模型的復(fù)雜度會(huì)增大很多。因此，需要并行計(jì)算才有可能實(shí)現(xiàn)實(shí)時(shí)濾波，從而將低秩模型應(yīng)用到實(shí)時(shí)的超聲系統(tǒng)中。

4 結(jié)論

在本文中，介紹了超聲血流成像中的發(fā)射波形、幀率、壁濾波等幾個(gè)關(guān)鍵參數(shù)的相互關(guān)系和影響，討論了血流信號(hào)的特點(diǎn)并對(duì)此提出了一種全新的低秩濾波方法，用來(lái)提取被雜波污染的血流信號(hào)。通過(guò)對(duì)比傳統(tǒng)的FIR濾波器，發(fā)現(xiàn)低秩濾波具有兩大主要優(yōu)勢(shì)：首先它可以提取更加精確的血流信號(hào)(對(duì)比FIR)，其次低秩濾波后的信號(hào)長(zhǎng)度保持不變(FIR濾波有信號(hào)長(zhǎng)度損失)。但是，低秩濾波比傳統(tǒng)的濾波器要復(fù)雜得多，尤其是如果針對(duì)更加復(fù)雜的仿真數(shù)據(jù)或者實(shí)際數(shù)據(jù)，低秩模型會(huì)變得更加復(fù)雜。由于低秩濾波的計(jì)算量太大，所以目前它還難以在實(shí)時(shí)的超聲成像系統(tǒng)上實(shí)現(xiàn)。

[1] Jensen JA. Estimation of blood velocities using ultrasound: A signal processing approach[M]//New York: Cambridge University Press, 1996:19-22.

[2] Nicholas D. Evaluation of backscattering coefficients for excited human tissue: results, interpretation and associated measurements[J]. Ultrasound in Medicine and Biology, 1982, 8(1):17-28.

[3] Fei DY, Shung KK. Ultrasonic backscatter from mammalian tissues[J]. Journal of the Acoustical Society of America, 1985, 78:871-876.

[4] Yuan YW, Shung KK. Ultrasonic backscatter from flowing whole blood. II: dependence on shear frequency and fibrinogen concentration[J]. Journal of the Acoustical Society of America, 1988, 84:1195-1200.

[6] Jensen JA, Medical ultrasound imaging[J]. Progress in Biophysics and Molecular Biology, 2007,93: 153-165.

[7] Kasai C, Namekawa K, Koyano A, et al. Real-time two-dimensional blood flow imaging using an autocorrelation technique[J]. IEEE Transactions on Sonics and Ultrasonics, 1985, 32(3): 458-464.

[8] Udesen J, Gran F, Hansen KL, et al. High frame-rate blood vector velocity imaging using plane waves: simulations and preliminary experiments[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2008, 55(8):1729-1743.

[9] Jensen JA, Svendsen NB. Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 1992, 39(2):262-267.

[10] Jensen JA. Field: A program for simulating ultrasound systems[J]. Medical and Biological Engineering and Computing, 1996, 34(Sup1-Part1):351-353.

[11] Jensen JA. Users’ guide for the Field II program[R]. Technical University of Denmark, Release 3.20, 2011.

[12] Hansen KL, Udesen J, Gran F, et al. Fast blood vector velocity imaging using ultrasound: in-vivo examples of complex blood flow in the vascular system[C]//Proceedings of IEEE International Ultrasonics Symposium. New York: IEEE, 2008:1068-1071.

[13] Leow CH, Bazigou E, Eckersley RJ, et al. Flow velocity mapping using contrast enhanced high-frame-rate plane wave ultrasound and image tracking: methods and initial in vitro and in vivo evaluation[J]. Ultrasound in Medicine and Biology, 2015, 41(11):2913-2925.

[14] Wright J, Ganesh A, Rao S, et al. Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[C]//Advances in Neural Information Processing Systems. Red Hook: Curran Associates, 2009:2080-2088.

[15] Candes EJ, Wakin MB. An introduction to compressive sampling[J]. IEEE Signal Processing Magazine, 2008, 25(2):21-30.

[16] Wright J, Yang AY, Ganesh A, et al. Robust face recognition via sparse representation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2009, 31(2):210-227.

[17] Fazel M, Hindi H, Boyd SP. A rank minimization heuristic with application to minimum order system approximation[C]//Proceedings of IEEE American Control Conference. New York: IEEE, 2001, 6:4734-4739.

[18] Candes EJ, Recht B. Exact matrix completion via convex optimization[J]. Foundations of Computational mathematics,2009, 9(6):717-772.

[19] Candes EJ, Tao T. The power of convex relaxation: Near-optimal matrix completion[J]. IEEE Transactions on Information Theory, 2010, 56(5):2053-2080.

[20] Chandrasekaran V, Sanghavi S, Parrilo PA, et al. Rank sparsity incoherence for matrix decomposition[J]. SIAM Journal on Optimization, 2011, 21(2):572-596.

[21] Ledoux LA, Brands PJ, Hoeks AP. Reduction of the clutter component in Doppler ultrasound signals based on singular value decomposition: a simulation study[J]. Ultrasonic Imaging, 1997, 19:1-18.

[22] Tao Q, Wang Y, Fish P, et al. The wall signal removal in Doppler ultrasound systems based on recursive PCA[J].Ultrasound in Medicine and Biology, 2004, 30(3):369-379.

[23] Mauldin FW, Lin D, Hossack JA. A singular value filter for rejection of stationary artifact in medical ultrasound[C]//Proceedings of IEEE International Ultrasonics Symposium, New York: IEEE, 2010:359-362.

[24] Yu ACH, Cobbold R. Single-ensemble-based eigen-processing methods for color flow imaging-Part I. The Hankel-SVD filter[J].IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2008, 55(3):559-572.

[25] Yu ACH, Cobbold R. Single-ensemble-based eigen-processing methods for color flow imaging - Part II. The Matrix Pencil estimator[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2008, 55(3):573-587.

[26] Yu ACH, Lovstakken L. Eigen-based clutter filter design for ultrasound color flow imaging: a review[J]. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, 2010, 57(5):1096-1111.

[27] Boyd S, Vandenberghe L. Convex Optimization[M] //New York: Cambridge University Press, 2004:168-169.

[28] Lin Z, Liu R, Su Z. Linearized alternating direction method with adaptive penalty for low rank representation[C] //Advances in Neural Information Processing Systems. Red Hook: Curran Associates, 2011:612-620.

[29] Lin Z, Chen M, Ma Y. The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices[R]. University of Illinois, UILU-ENG-09-2214, 2010.

Clutter Removing Filter for Ultrasound Blood Flow Imaging Based on Low Rank Model

Du Yigang1,2Zhang Mengyi3Chen Siping1*Li Yong2

1(National-RegionalKeyTechnologyEngineeringLaboratoryforMedicalUltrasound,DepartmentofBiomedicalEngineering,SchoolofMedicine,ShenzhenUniversity,Shenzhen518060,Guangdong,China)2(ShenzhenMindrayBio-MedicalElectronicsCo.Ltd.,Shenzhen518057,Guangdong,China)3(DepartmentofComputerScienceandEngineering,TheChineseUniversityofHongKong,Shatin,HongKong,China)

The conventional ultrasound blood flow wall filter uses a fixed cut-off frequency, which is not effective when the tissue motion is different due to heart beat and breath. This paper presented an ultrasound clutter removing filter based on a low rank model. The characteristics of ultrasound flow signal was studied and formulated. The low rank model was comprised of a rank minimization and matrix sparsity problem. The convex optimization can be applied to solve it after relaxation. The novelty is that it is an adaptive filter due to the minimization of the combination of the nuclear norm and L1norm. Ultrasound blood flow data were simulated. The filtered signals were obtained by three different orders FIR filters and the low rank filter. The RMS errors for FIR filtering were around 34%, 16% and 12% respectively, and lower than 0.001% when using the low rank filter, which not only improved the accuracy a lot but also maintained the same length of the filtered signal as the original one′s, where the length of the FIR filtered signal was decreased compared to the original signal. However, the low rank model is much more complicated than the conventional method, and it is still difficult to be applied in a real-time ultrasound imaging system.

ultrasound blood flow imaging; clutter removing filter; low rank model; FIR filter

10.3969/j.issn.0258-8021. 2016. 04.002

2016-02-25， 錄用日期:2016-05-10

中國(guó)博士后科學(xué)基金(2014M562207)；國(guó)家自然科學(xué)基金(61372006，61427806)

R318

A

0258-8021(2016) 04-0394-08

*通信作者(Corresponding author)， E-mail: chensiping@szu.edu.cn