基于稀疏脈沖采樣的低復(fù)雜度血流速度估計(jì)算法

馬碧云 吳港 劉嬌蛟 范逸華

基于稀疏脈沖采樣的低復(fù)雜度血流速度估計(jì)算法

馬碧云 吳港 劉嬌蛟 范逸華

（華南理工大學(xué) 電子與信息學(xué)院，廣東 廣州 510640）

雙模超聲廣泛用于醫(yī)學(xué)臨床診斷，其中B模式脈沖用于成像，多普勒脈沖則用于血流速度估計(jì)。數(shù)據(jù)采集時(shí)間在兩種模式之間共享。為了提高B模式圖像的更新頻率，需要減少多普勒脈沖數(shù)量，即發(fā)射稀疏多普勒脈沖進(jìn)行血流速度估計(jì)。然而現(xiàn)有的適應(yīng)稀疏脈沖采樣算法，如迭代自適應(yīng)算法、稀疏貝葉斯法以及基于陣列虛擬拓展的子空間類方法，計(jì)算開(kāi)銷巨大，難以滿足實(shí)時(shí)成像的要求，且在稀疏度大的情況下會(huì)產(chǎn)生明顯的偽影。為此，文中提出了一種基于稀疏脈沖采樣的低復(fù)雜度血流速度估計(jì)算法。根據(jù)超聲多普勒回波信號(hào)是由血紅細(xì)胞的散射產(chǎn)生，具有強(qiáng)相干、信源個(gè)數(shù)時(shí)變的特點(diǎn)，文中首先從子空間角度解析了偽影的成因，并驗(yàn)證了包含均勻脈沖的稀疏發(fā)射脈沖排布方式可以有效地抑制偽影；然后以均勻脈沖回波構(gòu)建協(xié)方差矩陣，并進(jìn)行空間平滑獲取特征值，以較大特征值的個(gè)數(shù)和相互的比值作為標(biāo)準(zhǔn)，判斷血流不同時(shí)刻的頻率分布特征；最后以此頻率分布特征為標(biāo)準(zhǔn)，自適應(yīng)采用B-MUSIC算法或TBVAM算法進(jìn)行血流速度估計(jì)，以降低算法的復(fù)雜度。Matlab仿真和人體實(shí)測(cè)數(shù)據(jù)的實(shí)驗(yàn)結(jié)果表明，該算法在極大地減小計(jì)算復(fù)雜度的同時(shí)，可以獲得較為連續(xù)、清晰且偽影抑制效果較佳的血流速度估計(jì)結(jié)果。

醫(yī)學(xué)超聲成像；血流速度估計(jì)；多普勒超聲信號(hào)；稀疏采樣

超聲多普勒技術(shù)可以實(shí)現(xiàn)非介入式的血流速度估計(jì)［1］，廣泛用于臨床診斷。醫(yī)學(xué)超聲系統(tǒng)采用雙模發(fā)射，其中B模式脈沖用于成像，多普勒脈沖用于血流速度估計(jì)。數(shù)據(jù)采集時(shí)間在兩種模式之間共享，需要合理地安排兩者發(fā)射間隔和發(fā)射頻次，即發(fā)射排布。為了提高B模式圖像的更新頻率，需要減少多普勒脈沖數(shù)量，即發(fā)射稀疏脈沖完成血流速度估計(jì)。因此，基于稀疏多普勒脈沖的血流速度高精度估計(jì)算法成為研究的熱點(diǎn)［2］。

Gran等［3-4］提出了使用血流迭代自適應(yīng)算法（BIAA）進(jìn)行血流速度估計(jì)，這種方式可以處理任意排布的多普勒發(fā)射脈沖，但需要進(jìn)行迭代運(yùn)算，計(jì)算量大，且需使用回歸濾波器去除雜波，會(huì)產(chǎn)生雜散頻率分量而降低頻譜估計(jì)的質(zhì)量。為了進(jìn)一步降低復(fù)雜度和提高計(jì)算精度，Lorintiu等［5］提出使用塊稀疏貝葉斯學(xué)習(xí)（BSBL）方法進(jìn)行血流速度估計(jì)，可采用隨機(jī)排布多普勒發(fā)射脈沖，然而其計(jì)算開(kāi)銷巨大，且在稀疏度大的情況下會(huì)產(chǎn)生明顯的偽影。為了抑制稀疏脈沖產(chǎn)生的偽影，Cohen等［6］利用多輸入多輸出雷達(dá)系統(tǒng)中嵌套陣列的研究成果，使用嵌套排布（NEST）的多普勒發(fā)射脈沖進(jìn)行血流速度估計(jì)，提出了基于傅里葉變換的NEST和基于NEST的旋轉(zhuǎn)不變性子空間算法（NESPRIT），用于血流速度估計(jì)。NEST具有較低的復(fù)雜性，并且在同等脈沖個(gè)數(shù)的情況下，由于使用了嵌套排布擴(kuò)展了孔徑，相較于以往的方法有更高的分辨率。NESPRIT由于采用了虛擬陣列技術(shù)，實(shí)現(xiàn)了較高分辨率的血流速度估計(jì)，但運(yùn)算復(fù)雜度較高，僅在發(fā)射嵌套排布的稀疏多普勒脈沖時(shí)才可獲得較好的估計(jì)效果，而其他的稀疏排布方式（如超嵌套排布、互質(zhì)排布、四階嵌套排布等）會(huì)產(chǎn)生嚴(yán)重的偽影，無(wú)法取得較好的估計(jì)結(jié)果，受限于其發(fā)射多普勒脈沖間隔，在相同稀疏發(fā)射情況下比經(jīng)典Welch和Capon等方法［7-8］的計(jì)算復(fù)雜度高。

針對(duì)上述問(wèn)題，為了在相同條件下獲得比Welch等方法較高的分辨率，文中首先從子空間的角度解釋了偽影的成因，并驗(yàn)證了一種抑制偽影的發(fā)射脈沖排布方式；然后根據(jù)血流回波不同時(shí)刻的頻率分布特征，提出了一種可用于稀疏脈沖的自適應(yīng)估計(jì)算法，并通過(guò)Matlab仿真和人體實(shí)測(cè)數(shù)據(jù)驗(yàn)證了該算法的有效性。

1　血流速度估計(jì)模型

對(duì)回波信號(hào)進(jìn)行協(xié)方差運(yùn)算，得到協(xié)方差矩陣

2　子空間類算法與偽影成因

基于譜估計(jì)的MUSIC算法，對(duì)協(xié)方差矩陣進(jìn)行特征值分解，可以得到噪聲子空間和信號(hào)子空間，即

于是可以通過(guò)頻率范圍內(nèi)的搜索得到頻譜，也就是血流速度譜，即

上述方法稱為B-MUSIC算法。

于是使用與B-MUSIC算法相似的步驟進(jìn)行特征值分解，求得噪聲子空間為

接著對(duì)其進(jìn)行譜搜索，可得到血流速度譜

該方法稱為T(mén)BVAM算法。

然而，除了嵌套排布外，稀疏發(fā)射脈沖排布使用各類算法均會(huì)出現(xiàn)偽影［6］。對(duì)于B-MUSIC算法，由于采樣容積［10］的存在，回波還有多種頻率分量，即多個(gè)信源。使用虛擬方式出現(xiàn)偽影，是由于頻率分量過(guò)多，導(dǎo)致其估計(jì)是有偏的［11］。

若僅采用物理發(fā)射脈沖，從子空間角度出發(fā)，要抑制偽影最簡(jiǎn)單的方式就是讓發(fā)射排布在前段有足夠數(shù)量的單位間隔均勻脈沖。如果將流型矩陣簡(jiǎn)單劃分為兩部分，即

圖1　第二峰值處的隨均勻脈沖長(zhǎng)度的變化

3　自適應(yīng)血流速度估計(jì)

血流速度隨著心跳呈現(xiàn)周期變化，當(dāng)血流速度較小時(shí)，其回波頻率分量會(huì)較為集中［12］，此時(shí)偽影與實(shí)際血流速度估計(jì)譜很容易區(qū)分，在源數(shù)少的情況下，用B-MUSIC算法的估計(jì)性能比用虛擬陣列方法更優(yōu)［13-14］。文中結(jié)合B-MUSIC和TBVAM算法的各自優(yōu)勢(shì)，提出了一種可以自適應(yīng)血流流速變化的低復(fù)雜度稀疏脈沖血流速度估計(jì)算法BSEA。

首先，為了區(qū)分不同時(shí)刻的流速擴(kuò)散情況，對(duì)均勻線陣處使用空間平滑，通過(guò)較大特征值的個(gè)數(shù)和相互的比值進(jìn)行血流流速的判斷。

以長(zhǎng)度為128的單位間隔均勻線性排列為例，對(duì)1 s內(nèi)血流信號(hào)進(jìn)行空間平滑后以軟閾值算子算法［15］進(jìn)行信源個(gè)數(shù)估計(jì)，結(jié)果如圖2所示，其信源個(gè)數(shù)隨時(shí)間變化的規(guī)律與血流速度變化規(guī)律一致。

圖2　軟閾值算子算法的信源個(gè)數(shù)估計(jì)結(jié)果

從圖2可以看出：對(duì)于血流速度較快區(qū)域，信源個(gè)數(shù)較多，即血流頻率分量較多，進(jìn)行特征分解會(huì)有較多較大的特征值；對(duì)于血流速度較慢的區(qū)域，頻率分類較少，只有少量的大特征值。因此，可以對(duì)最大特征值與第二大特征值進(jìn)行比較，若

由式（6）進(jìn)行特征值分解求特征值；

else

end if

4　仿真結(jié)果分析

圖3　最小冗余排布使用B-MUSIC算法的估計(jì)結(jié)果

對(duì)27個(gè)發(fā)射脈沖的最小冗余排布發(fā)射排序加上一個(gè)長(zhǎng)度為16的發(fā)射間隔均勻線陣，其總的物理發(fā)射脈沖數(shù)為37。同時(shí)，為了防止由于物理發(fā)射脈沖數(shù)目的增加而提升估計(jì)性能，文中隨機(jī)在27個(gè)發(fā)射脈沖的最小冗余排序中加入發(fā)射序列直至其達(dá)到37個(gè)物理脈沖。對(duì)這兩種發(fā)射排布分別使用B-MUSIC算法與TBVAM算法進(jìn)行血流速度估計(jì)，結(jié)果如圖4所示。

從圖4（a）和4（c）可以看出，采用B-MUSIC算法估計(jì)時(shí)，前端添加均勻線陣的偽影抑制效果明顯優(yōu)于隨機(jī)添加發(fā)射脈沖。從圖4（b）和4（d）可以看出，TBVAM算法相較于B-MUSIC算法，具有更好的偽影抑制效果，且估計(jì)精度更高，同時(shí)前端添加均勻脈沖會(huì)比隨機(jī)添加發(fā)射脈沖有更好的偽影抑制性能。

采用前端添加均勻線陣的最小冗余排布，分別使用BSBL算法與BSEA算法進(jìn)行血流速度估計(jì)，結(jié)果如圖5所示。從圖中可知，BSEA算法相對(duì)于同樣可以用于稀疏陣列的BSBL算法，偽影抑制效果更優(yōu)，且算法的估計(jì)精度更高，可與TBVAM算法的性能進(jìn)行比較。此外，同樣的仿真環(huán)境和稀疏發(fā)射脈沖排布，使用BSBL算法耗時(shí)為1 211.71 s，TBVAM算法耗時(shí)為28.41 s，而使用BSEA算法耗時(shí)為15.03 s，由此可見(jiàn)，使用BSEA算法可以有效地降低計(jì)算復(fù)雜度。

圖5　最小冗余排布添加均勻脈沖后的估計(jì)結(jié)果

最后采用體超聲多普勒回波的開(kāi)源數(shù)據(jù)［6］驗(yàn)證文中BSEA算法的實(shí)用性。該數(shù)據(jù)來(lái)自于一位33歲的健康男性，探測(cè)部位分別為頸部主動(dòng)脈和肝臟靜脈?；夭〝?shù)據(jù)采集方式為：用B-K8556 3.2 MHz線性陣列探頭進(jìn)行掃描，脈沖重復(fù)發(fā)射頻率為3.5 kHz，總脈沖個(gè)數(shù)為8 000，每個(gè)脈沖快拍數(shù)為128，采樣頻率15 MHz。首先對(duì)回波中由血管壁、組織反射所產(chǎn)生的低頻雜波部分進(jìn)行壁濾波［6，16］，此處壁濾波選取巴特沃茲高通濾波器，其截止頻率為80 Hz；然后進(jìn)行Hilbert變換，將實(shí)信號(hào)變?yōu)閺?fù)信號(hào)以便進(jìn)行血流速度估計(jì)；接著采用27個(gè)發(fā)射脈沖的最小冗余排布發(fā)射排序加上一個(gè)長(zhǎng)度為16的發(fā)射間隔均勻線陣作為發(fā)射排布，采用BSEA算法進(jìn)行血流速度估計(jì)，結(jié)果如圖6所示。實(shí)驗(yàn)結(jié)果表明，文中所提出的BSEA算法在實(shí)際應(yīng)用中，尤其對(duì)于血流頻率分量較多的動(dòng)脈回波數(shù)據(jù)，可以獲得較為連續(xù)、清晰且偽影抑制效果較佳的血流速度估計(jì)結(jié)果，具有實(shí)用性。

圖6　BSEA算法的人體實(shí)測(cè)數(shù)據(jù)估計(jì)結(jié)果

5　結(jié)語(yǔ)

文中使用子空間方法分別基于直接稀疏發(fā)射脈沖排布以及虛擬化后的脈沖排布進(jìn)行血流速度估計(jì)，從子空間的角度說(shuō)明了當(dāng)稀疏發(fā)射排布中存在一個(gè)均勻線陣時(shí)有助于抑制偽影。采用虛擬脈沖進(jìn)行血流速度估計(jì)，既可以有效抑制偽影，也可以提高估計(jì)精度，但會(huì)大大增加運(yùn)算復(fù)雜度。因此，文中提出了一種可用于稀疏脈沖的自適應(yīng)估計(jì)算法BSEA，該算法可根據(jù)血流速度的時(shí)變性、流速較大時(shí)與流速較小時(shí)不同的頻率分布特征，自適應(yīng)采用B-MUSIC算法或TBVAM算法，以降低計(jì)算復(fù)雜度。Matlab仿真和人體數(shù)據(jù)血流速度估計(jì)結(jié)果表明，BSEA算法可在抑制偽影的同時(shí)，極大地減小計(jì)算復(fù)雜度。

［1］ 趙冬冬，劉雪松，周凡，等．便攜式三維成像聲納的非網(wǎng)格化稀疏陣列技術(shù)［J］．華南理工大學(xué)學(xué)報(bào)（自然科學(xué)版），2019，47（2）：24-30.

ZHAO Dongdong，LIU Xuesong，ZHOU Fan，et al．Nongrid sprse array technology for portable 3D imaging sonar［J］．Journal of South China University of Technology（Natural Science Edition），2019，47（2）：24-30.

［2］ JENSEN J A，ARENDT J．Spectral velocity estimation in ultrasound using sparse data sets［J］．The Journal of the Acoustical Society of America，2006，120（1）：211-220.

［3］ GRAN F，JAKOBSSON A，JENSEN J A．Adaptive spectral Doppler estimation［J］．IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics，and Frequency Control，2009，56（4）：700-714.

［4］ GUDMUNDSON E，JAKOBSSON A，JENSEN J A，et al．Blood velocity estimation using ultrasound and spectral iterative adaptive approaches［J］．Signal Processing，2011，91（5）：1275-1283.

［5］ LORINTIU O，LIEBGOTT H，F(xiàn)RIBOULET D．Compressed sensing Doppler ultrasound reconstruction using block sparse Bayesian learning［J］．IEEE Transactions on Medical Imaging，2016，35（4）：978-987.

［6］ COHEN R，ELDAR Y C．Sparse Doppler sensing based on nested arrays［J］．IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics，and Frequency Control，2018，65（12）：2349-2364.

［7］ WELCH P．The use of fast Fourier transform for the estimation of power spectra：a method based on time averaging over short，modified periodograms［J］．IEEE Transactions on Audio and Electroacoustics，1967，15（2）：70-73.

［8］ TABATABAEI MAJD S M M，MOHAMMADZADEH ASL B．Adaptive spectral Doppler estimation based on the modified amplitude spectrum Capon［J］．IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics，and Frequency Control，2021，68（5）：1664-1675.

［9］ LIU C L，VAIDYANATHAN P P．Remarks on the spatial smoothing step in coarray MUSIC［J］．IEEE Signal Processing Letters，2015，22（9）：1438-1442.

［10］ AMBROGIO S，JANE A，ELIZABETH G，et al．Pulsed wave Doppler measurements of maximum velocity：dependence on sample volume size［J］．Ultrasound in Medicine & Biology，2022，48（1）：68-77.

［11］ WANG M，NEHORAI A．Coarrays，MUSIC，and the Cramér-Rao bound［J］．IEEE Transactions on Signal Processing，2017，65（4）：933-946.

［12］ 高蓮，張榆鋒，章克信，等．病變血管內(nèi)血流分布模型的研究進(jìn)展［J］．系統(tǒng)仿真學(xué)報(bào)，2014，26（1）：146-151.

GAO Lian，ZHANG Yu-feng，ZHANG Ke-xin，et al．Research progress of blood flow distribution models in diseased vessels［J］．Journal of System Simulation，2014，26（1）：146-151.

［13］ VAIDYANATHAN P P，PAL P．Direct-MUSIC on sparse arrays［C］∥ Proceedings of 2012 International Conference on Signal Processing and Communications.Bangalore：IEEE，2012：1-5.

［14］ PAL P，VAIDYANATHAN P P．Gridless methods for underdetermined source estimation［C］∥ Proceedings of the 48th Asilomar Conference on Signals，Systems and Computers．Pacific Grove：IEEE，2014：111-115.

［15］ PAL P，VAIDYANATHAN P P．Soft-thresholding for spectrum sensing with coprime samplers［C］∥ Proceedings of 2014 IEEE the 8th Sensor Array Multichannel Signal Processing Workshop．A Coruna：IEEE，2014：517-520.

［16］ 王錄濤，王微，金鋼．基于Hankel-SVD的非平穩(wěn)超聲血流成像雜波抑制技術(shù)研究［J］．電子與信息學(xué)報(bào)，2015，37（4）：830-835.

WANG Lu-tao，WANG Wei，JIN Gang．Non-stationary clutter rejection based on Hankel-SVD for ultrasound color flow imaging［J］．Journal of Electronics and Information Technology，2015，37（4）：830-835.

Low Complexity Blood Flow Velocity Estimation Algorithm via Sparse Pulse Sampling

（School of Electronic and Information Engineering，South China University of Technology，Guangzhou 510640，Guangdong，China）

Dual-mode ultrasound is widely used in medical clinical diagnosis. The B-mode pulse is used for imaging and Doppler pulse is used for blood flow velocity estimation. The data collection time is shared between the two modes. To improve the update frequency of B-mode image, it is necessary to reduce the number of Doppler pulses, that is, to estimate the blood flow velocity by sparse Doppler emissions. However, the existing algorithms for sparse pulse sampling, such as iterative adaptive algorithm, sparse Bayesian algorithm and subspace method based on array virtual expansion, are huge in expense and can not meet the requirements of real-time imaging. What’s more, they will lead to obvious artifacts in the case of large sparsity. Therefore, this paper proposed a low complexity blood flow velocity estimation algorithm via sparse pulse sampling. Based on the fact that ultrasonic Doppler echo signal is generated by the scattering of red blood cells, so echoes are strong coherence signals with time-variation sources number, this paper firstly explained the cause of artifacts from the perspective of subspace, and verified that the sparse emission pulse arrangement with uniform pulse can effectively suppress artifacts. Then the covariance matrix was constructed with uniform pulse echo, and the eigenvalues were obtained after spatial smoothing. The frequency distribution characteristics of blood flow at different segments were derived by the number of larger eigenvalues and the ratio of each other. Finally, based on the frequency distribution characteristics, the B-MUSIC algorithm or TBVAM algorithm was adaptively used for blood flow velocity estimation to reduce the complexity of the algorithm. The experimental results with Matlab simulation and human body measurement data show that the algorithm can obtain continuous, clear blood flow velocity estimation results with well artifact suppression while reducing the computational complexity significantly.

medical ultrasound imaging；blood velocity estimation；Doppler ultrasound signal；sparse sampling

Supported by the General Program of the Natural Science Foundation of Guangdong Province （2021A1515011842，2022A1515011830，2023A1515011420）

10.12141/j.issn.1000-565X.220380

2022?06?17

廣東省自然科學(xué)基金面上項(xiàng)目（2021A1515011842，2022A1515011830，2023A1515011420）；廣州市科技計(jì)劃項(xiàng)目（202102080352）

馬碧云（1982-），女，博士，副教授，主要從事超聲檢測(cè)、超聲通信研究。E-mail：eebyma@scut.edu

TP391.41；R445.1

1000-565X（2023）05-0063-07