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

杜宜綱 張夢一 陳思平 李 勇

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

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

杜宜綱1，2張夢一3陳思平1*李 勇2

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

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

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

引言

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

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

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

圖2 聚焦波與平面波發(fā)射聲場的比較結(jié)果。聚焦波的發(fā)射焦點為60 mm。采用Field II[9-10]仿真發(fā)射聲場，64個陣元，探頭中心頻為率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濾波器是這種掃描方式的一個合適的選擇。近年來，由于平面波超快成像的發(fā)展，成像的掃描方式從逐線發(fā)展到部分區(qū)域乃至整體一次掃描成像，這使得血流的顯示幀率有了很大的提高[8]，因此血流成像時同一個位置掃描次數(shù)也不再受到幀率的限制，這樣高階的FIR濾波器也可以用在血流信號的提取上。但是，由于采用的是平面波掃描，沒有發(fā)射聚焦，因此發(fā)射聲場的旁瓣會很高。如圖2所示，采用了Field II[9-10]仿真發(fā)射聲場，即使不在焦點深度的聚焦波聲場旁瓣也比平面波低很多，加之血流信號的幅值本身就比組織信號低10～27dB[1]，所以血管內(nèi)的回波很容易受到血管壁以及組織的回波干擾產(chǎn)生雜波信號。如圖3所示，采用了Field II函數(shù)calc_scat_multi[11]生成通道數(shù)據(jù)，平面波的雜波信號要比聚焦波的大很多。通過仿真計算，得到聚焦波發(fā)射時血管里的雜波平均要比平面波的小16.2 dB，相當于平面波的雜波信號強度是聚焦波的6～7倍。如果無法抑制這些雜波信號，那么平面波成像雖然可以提高幀率，但是血流的靈敏度和精度比起傳統(tǒng)方法都會有較大損失。

圖3 采用Field II仿真平面波和聚焦波成像。黑色部分為模擬血管無反射物，血管壁為強反射點，采用Field II函數(shù)calc_scat_multi[11]生成通道數(shù)據(jù)，然后分別做平面波波束合成和動態(tài)接收聚焦波束合成。平面波成像時血管內(nèi)平均雜波信號為-106.9 dB，動態(tài)接收聚焦成像后血管內(nèi)平均雜波信號為-123.1dB(這里只模擬了二維血管，因此雜波信號很低。血管形狀在Field II仿真中為長方體)。(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]，還有的研究采用了造影劑進一步提高平面波的成像效果[13]。本研究希望設(shè)計出一種新型的超聲血流壁濾波器，從而直接克服上述由于血流信號弱、信噪比低、雜波干擾大造成的血流精度和靈敏度的損失。這種濾波器基于一個低秩模型，是針對平面波發(fā)射而設(shè)計的。這種低秩模型可以用來區(qū)分數(shù)據(jù)中數(shù)值較小的變化量與數(shù)值較大的常量，它們分別對應了能量較弱的血流信號和相對變化很小且能量較強的雜波信號。因為這些雜波信號大多來自周圍運動較小的組織或者血管壁的強反射，所以雜波信號隨時間變化較小。雖然血流的信號很弱，但它隨時間變化較大?？梢?，血流與雜波的這些特點恰好符合了低秩模型理論。

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

1 方法

1.1 理論模型

b0,i=f0,i+c0,i

(1)

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

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

(2)

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

B0=F0+C0

(3)

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

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

subject toC+F=B0

(4)

式中，C和F分別代表了計算出來的組織和血流信號(式(3)中的C0和F0則是實際的組織和血流信號)。

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

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

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

subject toC+F=B0

(5)

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

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

1.2 仿真方法

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

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

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

(6)

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

2 結(jié)果

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

圖5 血流速度恒定為0.25 m/s時FIR與低秩濾波的結(jié)果比較。F0為實際血流信號，其他為濾波后的血流信號(血流信號強度基于表1所述的血流反射系數(shù)，無物理單位)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)時某一特定位置上的FIR與低秩濾波的結(jié)果。F0為實際血流信號，其他為濾波后的血流信號(血流信號強度基于表1所述的血流反射系數(shù)，無物理單位)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ù)的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

測量組別濾波方法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)時不同階的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

測量組別濾波方法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 討論

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

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

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

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

4 結(jié)論

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

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

R318

A

0258-8021(2016) 04-0394-08

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