超聲微血管成像中壁濾波器的設計直接影響著血流成像的分辨率。本文通過仿真數據實驗和人體腎臟實體數據成像實驗對比了傳統多項式回歸壁濾波器算法以及基于奇異值分解（SVD）的Full-SVD算法和RS-RSVD算法（一種基于隨機奇異值分解的隨機下采樣算法）。實驗結果表明，傳統的多項式回歸壁濾波器算法濾波效果有限，而Full-SVD算法和RS-RSVD算法可以更好地將微血流信號從組織或噪聲信號中提取出來。RS-RSVD算法在隨機分塊16次時與Full-SVD算法信噪比相同，對比噪聲比降低了2.05 dB，執行時間降低了90.41%。RS-RSVD算法能夠提升運行效率，更有利于高幀頻超聲微血管的實時成像。
引用本文： 王寶宇, 張淼, 劉瑞麟, 張石. 應用于超聲微血管成像的壁濾波算法對比研究. 生物醫學工程學雜志, 2022, 39(4): 740-748, 758. doi: 10.7507/1001-5515.202203032 復制
引言
在醫學影像學領域，對血流狀態的定量和定性分析是非常重要的工作。多普勒技術是當前超聲設備獲取人體血流信號的主要手段，可以實現對血流信號的識別以及血流運動速度的估計。通常高速、高能量的血流識別較為容易，而低速微細血流（血管直徑小于1 mm）由于血流速度與組織運動接近，造成血流信號的提取困難。因為超聲單次發射覆蓋范圍較窄，獲得一幀圖像所需的時間較長，所以為了保證圖像的實時性，多普勒采樣點有限（多為8～16之間），使得低速血流信號很難被精確分離，造成血流估計誤差大，成像分辨率低。隨著超聲快速成像技術的發展，高幀頻超聲成像可以在短時間內采集到數量巨大的時空樣本，這使微血管成像得到改善，極大地提高了多普勒對微細血流信號的敏感性[1-2]。在超聲成像系統中實現信號分離的模塊被稱作雜波濾波器（clutter filter），即壁濾波器。圖1展示了超聲彩色血流成像系統框架，壁濾波器的設計和實現是超聲血流信號處理的關鍵性技術問題。
圖1
超聲彩色血流成像系統框架
Figure1.
Framework of ultrasonic color blood flow imaging system
為了提高血流信號的檢測能力，學者們提出了很多辦法，主要目標是抑制血流中的雜波及組織信號，基于頻域分析的方法是最早采用的運動識別分離手段。其中有限脈沖響應（finite impulse response，FIR）和無限脈沖響應（infinite impulse response，IIR）濾波器[3-4]是最簡單的濾波方法，它們均是沿時間維度進行高通濾波，然而IIR濾波器需要較長的穩定時間，FIR濾波器則需要高階才能從血流中識別雜波信號。由于血流信號的運動速度與頻率信息并非具有分布一致性[5]，這使基于時間維度濾波的高通濾波器提取血流信號更加困難。
奇異值分解（singular value decomposition，SVD）已經應用于醫學超聲成像中，該方法被應用于組織和雜波的去除以及運動估計[6-7]。與僅基于時間的高通濾波器相比，基于奇異值的雜波濾波器可以利用組織和血液的散斑強度和時空特征的差異來更好地抑制雜波。在執行SVD后，組織信息通常保留在高奇異值中，而血液信號由于低的散射強度和時空相干性，通常存在于較低的奇異值中。然后，通過奇異值的映射關系進行低秩矩陣分離，可以區分出組織信號和血流信號。圖2給出了組織信號和血流信號的頻率分布及信號強度對比。在超聲回波數據中，通常組織信號比血流信號強40～100 dB且具有更高的時間和空間相干性[8-10]。
圖2
組織信號和血流信號對比
Figure2.
Comparison of wall signal and blood flow signal
Yu等[11-12]學者提出一種基于矩陣束（matrix pencil，MP）的方法，首先利用單個采樣數據模擬多個采樣數據集生成Hankel矩陣，并對Hankel進行了SVD，即Hankel-SVD。Forsberg等[13]將SVD和ARMA模型相結合得到SVD-ARMA方法。Demené等[14]提出將基于特征的雜波濾波與平面波成像相結合，用于高靈敏度功能超聲成像。Mauldin等[15]提出了基于主成分分析（principal component analysis，PCA）的濾波器設計的框架。Song等[16]為提升小血管血流分辨率提出了塊自適應局部雜波濾波算法，但該方法執行時間是超聲系統無法承擔的。為加速奇異值分解雜波濾波器算法，Song等[17]同時提出了隨機奇異值分解算法，即RSVD算法。
由于SVD可以處理時間或空間數據，為了對射頻數據進行SVD處理，深度和橫向空間通常以重疊的形式組合為2D數據空間，對新的2D數據進行SVD處理，即Full-SVD算法。該方法同時包含數據時空特性，因此可以獲得更好的濾波性能。RSVD算法是使用QR分解得到近似Full-SVD算法的空間基向量以進行濾波處理。本文根據RSVD算法的高效性，結合隨機下采樣方案，降低系統處理時間，提出了基于隨機奇異值分解的隨機下采樣算法，即RS-RSVD算法，在盡可能高的信噪比下加速濾波處理。
本文首先介紹傳統的多項式回歸壁濾波器算法，進一步對Full-SVD算法和RS-RSVD算法進行詳細闡述，并通過實驗從信噪比、對比噪聲比和執行運算時間等方面對三種算法進行綜合對比，分析并總結出更適用于微血管成像的壁濾波算法。人體的微血管分布廣泛，但結構復雜，使用仿體很難進行模擬，因此需要實際采集人體微細血流數據進行驗證，本次實驗實體數據由東軟醫療系統股份有限公司提供及授權。
1 方法
1.1 傳統多項式回歸壁濾波器算法
高通濾波器是一種傳統的壁濾波器，主要包括三種類型，FIR、IIR和回歸型。多項式回歸壁濾波器[18-19]是目前商用的主流壁濾波器，其根據多項式的正交分解和矩陣運算進行濾波處理。多項式回歸壁濾波器采用正交勒讓德多項式，從原始I/Q或RF數據中刪除擬合多項式，以保留高頻分量及血流信號部分。
 |
x代表原始信號，y代表血流信號，c代表雜波信號；K代表雜波維度；bk為正交向量基。本文使用勒讓德多項式是通過斯密特正交化得到{1, x, x2, 
, xn}，通過創建一個變換矩陣A，A矩陣的每一列定義為序列x的冪次方，序列x為[? n : n]，n為時間序列長度，用于壁濾波器的構建。
多項式回歸壁濾波器頻域響應可表示為：
 |
其中：
 |
1.2 Full-SVD算法
Full-SVD雜波濾波器是利用接收到的超聲數據豐富的時間和空間信息，根據投影到奇異值域時信號不同分量的不同特征，將血流信號、組織雜波和噪聲進行魯棒分離[20]。一般認為，組織呈高回聲，不能快速移動，而血液呈低回聲并快速振蕩。奇異值分解能夠抑制組織雜波和血液信號重疊的頻譜，這是傳統的高通濾波器無法實現的。
在彩色血流成像中，超聲回波數據為圖3所示的三維數據，假設超聲回波數據S是一個x*y*t的三維矩陣，x為橫向尺寸，y為軸向尺寸，t對應于時間維度。將它重新定義為一個維度為xy*t的重塑矩陣M，利用M進行SVD分解。
圖3
超聲回波數據
depth代表超聲檢測深度，slow-time代表連續時間采集序列，channel代表超聲成像區域的橫向距離
Figure3. Ultrasonic echo datadepth represents ultrasonic detection depth, slow-time represents continuous time acquisition sequence and channel represents transverse distance of ultrasonic imaging area
 |
其中，U為m*n左奇異向量，Σ = diag(σ1…σn)是一個對角矩陣，對角線上的元素為奇異值，V為n*n右奇異向量。具體的展開形式：
 |
σn稱為A的奇異值，它們按照遞減的順序排列：
 |
可以得到一個秩k，該秩接近最佳奇異值ε。通過選取了一個合適的秩σk[21-22]，
 |
 |
 |
 |
即AT表示組織雜波。血流信號表示為：
 |
在Full-SVD算法中，利用奇異值的映射關系找到血流信號、噪聲信號、組織信號的分布情況，假設原始信號經過奇異值分解后有N個奇異值，根據組織信號與血流信號的分布特點可以找到信號轉折點k，即前k階被認為是回波信號主要成分，即組織信號。血流信息主要存在于σk + 1
σn個奇異值中，通過截取該部分的奇異值，即可得到濾除組織和噪聲的血流信號，以達到濾波效果。
1.3 RS-RSVD算法
RS-RSVD算法是根據文獻[17]中所提出的RSVD算法結合圖4所示的隨機采樣方案進行濾波處理。假設在超聲血流系統中，接收到的矩陣信號是272*288*41的三維數據集（在實際系統中采集數據遠大于此）用于一幀成像，經過降維重排，三維數據集降維處理重塑為一個二維78 336*41的數據集，當采集到的數據是16位定點數時，78 336*41的數據量大約需要6 423 552 B（即6.12 MB）的內存資源，假設隨機采樣分4組，則每組數據量降低為19 584*41（即1.53 MB）。本次實驗中，我們通過將聚焦點拉長（數據顯示深度2～3倍距離），實現較大覆蓋范圍的發射聲場，在一定寬度均勻波束面范圍內采集并合成回波數據。成像幀頻會受超聲的探測深度和發射角度的影響，表1給出了超聲探測深度和成像角度對應的高幀頻超聲系統的采集幀頻。
圖4
隨機下采樣
Figure4.
Random down sampling
表1
深度和角度對幀頻的限制
Table1.
Limitation of depth and angle on frame rate
將隨機采樣的回波信號形成重塑矩陣M1（xy*t），假設組織雜波逼近前k階奇異值，Ω（t*k）為滿足N（0，1）的標準正態分布的隨機矩陣。
 '/> |
經過矩陣乘法，
為xy*k二維矩陣。
 '/> |
qr表示QR分解，可以找到列空間近似標準正交基。
 |
*表示共軛轉置，本文通過一次迭代循環得出Q [23-24]。
 |
組織雜波信號T表示為：
 |
血流信號F表示為：
 |
RSVD算法通過一個k階的標準正態分布矩陣結合重塑矩陣M1生成 
，進行QR分解生成一個近似標準正交基，此標準正交基即為組織信號，通過一次迭代計算恢復組織信號，用原始信號減去組織信號得到濾除雜波信號的血流信號。圖5展示了Full-SVD、RSVD不進行迭代以及RSVD迭代一次的實驗誤差對比。根據實驗結果可以得到，通過一次QR分解迭代計算的擬合曲線和Full-SVD濾波結果基本重合。文獻[17]表明分塊下采樣與交叉下采樣會使生成圖像產生偽影。所以本文利用RSVD算法結合圖4所示隨機下采樣信號抽取方案，可以減小QR分解所計算矩陣大小以此縮短矩陣分解的執行時間，從而降低微血流成像執行的處理時間。
圖5
Full-SVD、RSVD不迭代、RSVD迭代一次對比
Figure5.
Comparison of Full SVD, RSVD-no iteration and RSVD- one iteration
2 實驗結果
為了驗證不同濾波器的濾波效果，本文使用Intel Xeon E5-2620 CPU、主頻2.1 GHz、內存32 GB的工作站進行實驗。在算法仿真方面使用丹麥理工大學Jenson團隊開發的Field II超聲平臺[25-26]模擬血流與雜波信號進行仿真實驗。為充分驗證各濾波器算法的濾波效果，使用東軟醫療系統股份有限公司提供的人體部分腎臟血流實體數據進行實驗。根據實驗結果，對比了傳統多項式回歸壁濾波器算法、Full-SVD壁濾波器算法、RS-RSVD壁濾波器算法對血流信號的提取效果。
2.1 仿真數據成像實驗
由于呼吸、心跳、組織運動會對血管造成影響，仿真實驗將上述因素考慮其中，使用Field II生成符合人體的微細血流散射體。設置血液密度1.06 × 103 kg/m3，血液黏度0.004 kg/（m·s），血管半徑1 mm，假設血流平均速度4 cm/s，心率62次/分，其中探頭選用凸陣探頭，聲速1 540 m/s，中心頻率為7.5 MHz，采樣頻率設置為105 MHz，陣元數為60，脈沖重復頻率（pulse repetition frequency，PRF）為3 500 Hz。聚焦中心設置在x = 0，y = 0，z = 4 cm，通過200次發射采集到的微細血流數據大小為1 158*200。
圖6a展示了模擬微細血管，1 158個不同深度連續采樣點和200個連續時間序列，圖6b顯示了模擬微細血流信號和組織信號回波數據，結果顯示在400～800采樣區間為組織信號與血流信號的分布區域。
圖6
單點數據集微細血流仿真實驗回波
a. 功率多普勒成像；b. 時域
Figure6. Echo of micro-flow simulation experiment of single point data seta. power Doppler imaging; b. time domain
圖7分別展示了多項式回歸壁濾波器和Full-SVD、RS-RSVD壁濾波器算法針對微細血流的提取效果。當多項式回歸壁濾波器算法在連續采集時間序列大于155時，由于實驗平臺無法創建矩陣規格大于156*156的勒讓德多項式，實驗將降低連續時間采集序列，將實驗數據降低為1 158*155進行實驗。
圖7
仿真數據實驗血流提取成像
a. 多項式回歸壁濾波器；b. Full-SVD；c. RS-RSVD
Figure7. Blood flow extraction imaging of simulation data experimentala. polynomial regression wall filter; b. Full-SVD; c RS-RSVD
圖8a展示了模擬微細血流回波信號的奇異值分布情況，在Full-SVD壁濾波器中選取（4，?368）為轉折點進行濾波。圖8b展示了多項式回歸壁濾波器、Full-SVD、RS-RSVD對微細血流信號提取效果，根據實驗結果，在仿真實驗中三種濾波器算法均能有效地將血流信號提取出來。
圖8
Full-SVD算法分解仿真數據的奇異值分布及三種濾波算法提取血流信號強度對比
a. 仿真實驗奇異值排布；b. 多項式回歸壁濾波器、Full-SVD、RS-RSVD濾波效果對比
Figure8. The singular value distribution of simulation data decomposed by Full-SVD algorithm and comparison of the intensity of blood flow signals extracted by the three filtering algorithmsa. singular value arrangements of simulation experiment; b. comparison of filtering effects of polynomial regression wall filter, Full-SVD and RS-RSVD
2.2 實體數據成像實驗
實體數據成像實驗采用人體腎臟部分實體回波數據來驗證不同濾波器算法的濾波效果。實驗探頭選用SL14-3H探頭類型，中心頻率8.17 MHz，陣元數量192，聲速1 540 m/s，采樣頻率30 MHz，線密度0.5，感興趣區（region of interest，ROI）為13.96 mm*12.80 mm，連續時間序列為41，數據大小為272*128*41。
實體數據的連續時間序列是41，實驗設定勒讓德多項式大小為41*41，多項式回歸壁濾波器得到的實際濾波效果如圖9所示。結果表明，多項式回歸濾波器選取的階數適用于仿體實驗，但并不適用于實體數據濾波處理，并且由于多項式回歸濾波器的頻率響應關于零頻率軸對稱，需求的多項式階數較大，所以在濾除雜波信號時，錯誤地將血流信號濾除，結果無法有效提取微細血流信號。
圖9
多項式回歸壁濾波器處理的功率多普勒血流圖像
Figure9.
Power Doppler blood flow image processed by polyno mial regression wall filter
圖10a為Full-SVD壁濾波器算法根據實體數據得到的奇異值分布情況，根據奇異值排布曲線的轉折點來確定組織雜波截止點A。假設數據中存在高斯噪聲分布，噪聲對應的高階奇異值在對數下呈線性分布。對超過15階奇異值曲線進行線性擬合，可以得到奇異值曲線偏離擬合曲線的轉折點B。確定A、B兩點即為組織和噪聲在奇異值映射范圍內的截止階數。根據截止點A、B進行雜波抑制，得到如圖10b所示的部分腎臟微細血流功率多普勒圖像。實驗結果表明，相比多項式回歸型濾波器算法，Full-SVD算法可以有效提高微細血流信號的提取。
圖10
Full-SVD算法分解實體數據的奇異值分布和Full-SVD處理完成的功率多普勒血流圖像
a. 腎臟回波數據奇異值分布；b. Full-SVD處理完成的功率多普勒血流圖像
Figure10. The singular value distribution of entity data decomposed by Full-SVD algorithm and the power Doppler flow image completed by Full-SVD processinga. singular value distribution of kidney echo data; b. power Doppler blood flow image processed by Full-SVD
為驗證RS-RSVD算法的成像效果，我們調用MATLAB工具箱中的rand函數，對原始矩陣272*128*41進行隨機抽取，得到四個大小為136*64*41的隨機數組，對四個數組分別進行RSVD濾波實驗，得到如圖11a所示的濾波效果。RS1-RSVD、RS2-RSVD、RS3-RSVD、RS4-RSVD分別為四次抽取進行RSVD的濾波效果。將四個136*64*41矩陣數進行組合并恢復至原始矩陣位置，得到原始圖像272*128*41的超聲回波數據濾波結果，此即為濾除組織雜波的血流信號。圖11b為恢復原始圖像的功率多普勒血流圖像。假設算法在并行超聲系統中，當ROI為固定區域時，抽取分塊的數量將由并行系統可同時執行的流水線決定，當系統可以同時進行16通道算法處理，則原始矩陣可被分為16個68*32*41大小的矩陣。
圖11
RS-RSVD算法分塊濾波和恢復的功率多普勒血流圖像
a. RS-RSVD算法分塊濾波的功率多普勒圖像；b. 恢復的功率多普勒血流圖像
Figure11. RS-RSVD algorithm block filtering and restored power Doppler blood flow imagea. power Doppler images processed by RS-RSVD algorithm block filtering; b. restored power Doppler blood flow imaging
除了對三種算法形成的功率多普勒圖像進行視覺觀察外，本實驗將三種算法的信噪比（signal to noise ratio，SNR）和噪聲比（contrast to noise ratio，CNR）進行比較，結果見圖12。圖12藍色區域為分別選取的部分ROI，x = 1.5～5.0 mm，y =2.0～5.5 mm。SNR和CNR的具體表達式如下：
圖12
Full-SVD和RS-RSVD區域選取
Figure12.
Full-SVD and RS-RSVD area selection
 |
Sblood代表血流信號的平均強度（圖12綠色方框中白線部分），Nnoise代表背景噪聲信號的平均強度（圖12中黃色方框）。
 |
Sbackground代表信號平均強度（圖12中綠色方框）。
表2給出了三種算法的運算時間、SNR和CNR，由于多項式回歸壁濾波器算法無法有效提取微細血流信號，本次實驗只對比了Full-SVD算法和RS-RSVD算法的SNR和CNR。運算時間包括數據讀取和矩陣運算。根據實驗結果，多項式回歸壁濾波器算法處理時間最短，Full-SVD所處理的時間最長達48.50 s左右，RS-RSVD在4次采樣時處理時間為19.03 s左右。RS-RSVD算法相比Full-SVD算法SNR基本不變，CNR有所下降。
表2
三種算法運行時間、SNR和CNR
Table2.
Running time, SNR and CNR of three algorithms
為評估RS-RSVD不同抽取次數對算法帶來的影響，分別進行了9次、16次、25次、36次抽樣處理實驗。不同采樣次數的執行時間、SNR和CNR如表3所示。根據實驗結果，在分塊次數達16次時，微血流信號的信噪比與Full-SVD算法接近，此時濾波的執行時間為（4.65 ± 0.50）s。
表3
RS-RSVD不同抽樣次數的運行時間、SNR、CNR
Table3.
Running time, SNR and CNR of RS-RSVD with different sampling times
3 討論
本文通過仿真數據和實際人體腎臟超聲回波數據進行算法的對比實驗。實驗結果表明，傳統的多項式回歸壁濾波器算法相比Full-SVD算法和RS-RSVD算法濾波效果較差，在提取微細血流信號時會出現錯誤情況。在高幀頻超聲成像系統中，回波信號的矩陣大小會隨著采樣頻率和ROI的增加而成倍增大，這使Full-SVD算法需要的資源空間更大，并且當矩陣數據量增加時奇異值分解所需的處理時間會呈倍數增加，使Full-SVD算法無法應用在實時的超聲檢測系統中，有一定的局限性。根據表3實驗結果，RS-RSVD算法相比多項式回歸濾波器算法具有可觀的視覺分辨率，與Full-SVD算法信噪比接近，在16次抽樣時雖然對比噪聲比降低了2.05 dB，但執行時間降低了90.52%。
4 結論
在超聲彩色血流成像中，壁濾波器的性能直接影響著血流信號的提取。臨床診斷中，微血流檢測具有重要的價值。近幾年，高性能的硬件被廣泛應用于超聲成像系統中，如圖形處理器（graphics processing unit，GPU）、現場可編程門陣列（field programmable gate array，FPGA）等。RS-RSVD算法為超聲并行處理提供了思路，并行處理減少成像時間將作為未來的研究重點。在超聲臨床醫學中，由于人體結構的差異和實際環境造成的影響，RS-RSVD算法還需在臨床上進行多次試驗，爭取早日應用于高幀頻超聲微血管成像系統中。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：王寶宇、張石提出研究思路，設計研究方案；張淼、劉瑞麟負責進行實驗；王寶宇、張淼負責采集和分析數據，同時負責論文起草及最終版本修訂。
倫理聲明：本研究實體數據由東軟醫療股份有限公司提供及授權。
引言
在醫學影像學領域，對血流狀態的定量和定性分析是非常重要的工作。多普勒技術是當前超聲設備獲取人體血流信號的主要手段，可以實現對血流信號的識別以及血流運動速度的估計。通常高速、高能量的血流識別較為容易，而低速微細血流（血管直徑小于1 mm）由于血流速度與組織運動接近，造成血流信號的提取困難。因為超聲單次發射覆蓋范圍較窄，獲得一幀圖像所需的時間較長，所以為了保證圖像的實時性，多普勒采樣點有限（多為8～16之間），使得低速血流信號很難被精確分離，造成血流估計誤差大，成像分辨率低。隨著超聲快速成像技術的發展，高幀頻超聲成像可以在短時間內采集到數量巨大的時空樣本，這使微血管成像得到改善，極大地提高了多普勒對微細血流信號的敏感性[1-2]。在超聲成像系統中實現信號分離的模塊被稱作雜波濾波器（clutter filter），即壁濾波器。圖1展示了超聲彩色血流成像系統框架，壁濾波器的設計和實現是超聲血流信號處理的關鍵性技術問題。
圖1
超聲彩色血流成像系統框架
Figure1.
Framework of ultrasonic color blood flow imaging system
為了提高血流信號的檢測能力，學者們提出了很多辦法，主要目標是抑制血流中的雜波及組織信號，基于頻域分析的方法是最早采用的運動識別分離手段。其中有限脈沖響應（finite impulse response，FIR）和無限脈沖響應（infinite impulse response，IIR）濾波器[3-4]是最簡單的濾波方法，它們均是沿時間維度進行高通濾波，然而IIR濾波器需要較長的穩定時間，FIR濾波器則需要高階才能從血流中識別雜波信號。由于血流信號的運動速度與頻率信息并非具有分布一致性[5]，這使基于時間維度濾波的高通濾波器提取血流信號更加困難。
奇異值分解（singular value decomposition，SVD）已經應用于醫學超聲成像中，該方法被應用于組織和雜波的去除以及運動估計[6-7]。與僅基于時間的高通濾波器相比，基于奇異值的雜波濾波器可以利用組織和血液的散斑強度和時空特征的差異來更好地抑制雜波。在執行SVD后，組織信息通常保留在高奇異值中，而血液信號由于低的散射強度和時空相干性，通常存在于較低的奇異值中。然后，通過奇異值的映射關系進行低秩矩陣分離，可以區分出組織信號和血流信號。圖2給出了組織信號和血流信號的頻率分布及信號強度對比。在超聲回波數據中，通常組織信號比血流信號強40～100 dB且具有更高的時間和空間相干性[8-10]。
圖2
組織信號和血流信號對比
Figure2.
Comparison of wall signal and blood flow signal
Yu等[11-12]學者提出一種基于矩陣束（matrix pencil，MP）的方法，首先利用單個采樣數據模擬多個采樣數據集生成Hankel矩陣，并對Hankel進行了SVD，即Hankel-SVD。Forsberg等[13]將SVD和ARMA模型相結合得到SVD-ARMA方法。Demené等[14]提出將基于特征的雜波濾波與平面波成像相結合，用于高靈敏度功能超聲成像。Mauldin等[15]提出了基于主成分分析（principal component analysis，PCA）的濾波器設計的框架。Song等[16]為提升小血管血流分辨率提出了塊自適應局部雜波濾波算法，但該方法執行時間是超聲系統無法承擔的。為加速奇異值分解雜波濾波器算法，Song等[17]同時提出了隨機奇異值分解算法，即RSVD算法。
由于SVD可以處理時間或空間數據，為了對射頻數據進行SVD處理，深度和橫向空間通常以重疊的形式組合為2D數據空間，對新的2D數據進行SVD處理，即Full-SVD算法。該方法同時包含數據時空特性，因此可以獲得更好的濾波性能。RSVD算法是使用QR分解得到近似Full-SVD算法的空間基向量以進行濾波處理。本文根據RSVD算法的高效性，結合隨機下采樣方案，降低系統處理時間，提出了基于隨機奇異值分解的隨機下采樣算法，即RS-RSVD算法，在盡可能高的信噪比下加速濾波處理。
本文首先介紹傳統的多項式回歸壁濾波器算法，進一步對Full-SVD算法和RS-RSVD算法進行詳細闡述，并通過實驗從信噪比、對比噪聲比和執行運算時間等方面對三種算法進行綜合對比，分析并總結出更適用于微血管成像的壁濾波算法。人體的微血管分布廣泛，但結構復雜，使用仿體很難進行模擬，因此需要實際采集人體微細血流數據進行驗證，本次實驗實體數據由東軟醫療系統股份有限公司提供及授權。
1 方法
1.1 傳統多項式回歸壁濾波器算法
高通濾波器是一種傳統的壁濾波器，主要包括三種類型，FIR、IIR和回歸型。多項式回歸壁濾波器[18-19]是目前商用的主流壁濾波器，其根據多項式的正交分解和矩陣運算進行濾波處理。多項式回歸壁濾波器采用正交勒讓德多項式，從原始I/Q或RF數據中刪除擬合多項式，以保留高頻分量及血流信號部分。
|
x代表原始信號，y代表血流信號，c代表雜波信號；K代表雜波維度；bk為正交向量基。本文使用勒讓德多項式是通過斯密特正交化得到{1, x, x2, , xn}，通過創建一個變換矩陣A，A矩陣的每一列定義為序列x的冪次方，序列x為[? n : n]，n為時間序列長度，用于壁濾波器的構建。
多項式回歸壁濾波器頻域響應可表示為：
|
其中：
|
1.2 Full-SVD算法
Full-SVD雜波濾波器是利用接收到的超聲數據豐富的時間和空間信息，根據投影到奇異值域時信號不同分量的不同特征，將血流信號、組織雜波和噪聲進行魯棒分離[20]。一般認為，組織呈高回聲，不能快速移動，而血液呈低回聲并快速振蕩。奇異值分解能夠抑制組織雜波和血液信號重疊的頻譜，這是傳統的高通濾波器無法實現的。
在彩色血流成像中，超聲回波數據為圖3所示的三維數據，假設超聲回波數據S是一個x*y*t的三維矩陣，x為橫向尺寸，y為軸向尺寸，t對應于時間維度。將它重新定義為一個維度為xy*t的重塑矩陣M，利用M進行SVD分解。
圖3
超聲回波數據
depth代表超聲檢測深度，slow-time代表連續時間采集序列，channel代表超聲成像區域的橫向距離
Figure3. Ultrasonic echo datadepth represents ultrasonic detection depth, slow-time represents continuous time acquisition sequence and channel represents transverse distance of ultrasonic imaging area
|
其中，U為m*n左奇異向量，Σ = diag(σ1…σn)是一個對角矩陣，對角線上的元素為奇異值，V為n*n右奇異向量。具體的展開形式：
|
σn稱為A的奇異值，它們按照遞減的順序排列：
|
可以得到一個秩k，該秩接近最佳奇異值ε。通過選取了一個合適的秩σk[21-22]，
|
|
|
|
即AT表示組織雜波。血流信號表示為：
|
在Full-SVD算法中，利用奇異值的映射關系找到血流信號、噪聲信號、組織信號的分布情況，假設原始信號經過奇異值分解后有N個奇異值，根據組織信號與血流信號的分布特點可以找到信號轉折點k，即前k階被認為是回波信號主要成分，即組織信號。血流信息主要存在于σk + 1σn個奇異值中，通過截取該部分的奇異值，即可得到濾除組織和噪聲的血流信號，以達到濾波效果。
1.3 RS-RSVD算法
RS-RSVD算法是根據文獻[17]中所提出的RSVD算法結合圖4所示的隨機采樣方案進行濾波處理。假設在超聲血流系統中，接收到的矩陣信號是272*288*41的三維數據集（在實際系統中采集數據遠大于此）用于一幀成像，經過降維重排，三維數據集降維處理重塑為一個二維78 336*41的數據集，當采集到的數據是16位定點數時，78 336*41的數據量大約需要6 423 552 B（即6.12 MB）的內存資源，假設隨機采樣分4組，則每組數據量降低為19 584*41（即1.53 MB）。本次實驗中，我們通過將聚焦點拉長（數據顯示深度2～3倍距離），實現較大覆蓋范圍的發射聲場，在一定寬度均勻波束面范圍內采集并合成回波數據。成像幀頻會受超聲的探測深度和發射角度的影響，表1給出了超聲探測深度和成像角度對應的高幀頻超聲系統的采集幀頻。
圖4
隨機下采樣
Figure4.
Random down sampling
表1
深度和角度對幀頻的限制
Table1.
Limitation of depth and angle on frame rate
將隨機采樣的回波信號形成重塑矩陣M1（xy*t），假設組織雜波逼近前k階奇異值，Ω（t*k）為滿足N（0，1）的標準正態分布的隨機矩陣。
| '/> |
經過矩陣乘法， 為xy*k二維矩陣。
| '/> |
qr表示QR分解，可以找到列空間近似標準正交基。
|
*表示共軛轉置，本文通過一次迭代循環得出Q [23-24]。
|
組織雜波信號T表示為：
|
血流信號F表示為：
|
RSVD算法通過一個k階的標準正態分布矩陣結合重塑矩陣M1生成 ，進行QR分解生成一個近似標準正交基，此標準正交基即為組織信號，通過一次迭代計算恢復組織信號，用原始信號減去組織信號得到濾除雜波信號的血流信號。圖5展示了Full-SVD、RSVD不進行迭代以及RSVD迭代一次的實驗誤差對比。根據實驗結果可以得到，通過一次QR分解迭代計算的擬合曲線和Full-SVD濾波結果基本重合。文獻[17]表明分塊下采樣與交叉下采樣會使生成圖像產生偽影。所以本文利用RSVD算法結合圖4所示隨機下采樣信號抽取方案，可以減小QR分解所計算矩陣大小以此縮短矩陣分解的執行時間，從而降低微血流成像執行的處理時間。
圖5
Full-SVD、RSVD不迭代、RSVD迭代一次對比
Figure5.
Comparison of Full SVD, RSVD-no iteration and RSVD- one iteration
2 實驗結果
為了驗證不同濾波器的濾波效果，本文使用Intel Xeon E5-2620 CPU、主頻2.1 GHz、內存32 GB的工作站進行實驗。在算法仿真方面使用丹麥理工大學Jenson團隊開發的Field II超聲平臺[25-26]模擬血流與雜波信號進行仿真實驗。為充分驗證各濾波器算法的濾波效果，使用東軟醫療系統股份有限公司提供的人體部分腎臟血流實體數據進行實驗。根據實驗結果，對比了傳統多項式回歸壁濾波器算法、Full-SVD壁濾波器算法、RS-RSVD壁濾波器算法對血流信號的提取效果。
2.1 仿真數據成像實驗
由于呼吸、心跳、組織運動會對血管造成影響，仿真實驗將上述因素考慮其中，使用Field II生成符合人體的微細血流散射體。設置血液密度1.06 × 103 kg/m3，血液黏度0.004 kg/（m·s），血管半徑1 mm，假設血流平均速度4 cm/s，心率62次/分，其中探頭選用凸陣探頭，聲速1 540 m/s，中心頻率為7.5 MHz，采樣頻率設置為105 MHz，陣元數為60，脈沖重復頻率（pulse repetition frequency，PRF）為3 500 Hz。聚焦中心設置在x = 0，y = 0，z = 4 cm，通過200次發射采集到的微細血流數據大小為1 158*200。
圖6a展示了模擬微細血管，1 158個不同深度連續采樣點和200個連續時間序列，圖6b顯示了模擬微細血流信號和組織信號回波數據，結果顯示在400～800采樣區間為組織信號與血流信號的分布區域。
圖6
單點數據集微細血流仿真實驗回波
a. 功率多普勒成像；b. 時域
Figure6. Echo of micro-flow simulation experiment of single point data seta. power Doppler imaging; b. time domain
圖7分別展示了多項式回歸壁濾波器和Full-SVD、RS-RSVD壁濾波器算法針對微細血流的提取效果。當多項式回歸壁濾波器算法在連續采集時間序列大于155時，由于實驗平臺無法創建矩陣規格大于156*156的勒讓德多項式，實驗將降低連續時間采集序列，將實驗數據降低為1 158*155進行實驗。
圖7
仿真數據實驗血流提取成像
a. 多項式回歸壁濾波器；b. Full-SVD；c. RS-RSVD
Figure7. Blood flow extraction imaging of simulation data experimentala. polynomial regression wall filter; b. Full-SVD; c RS-RSVD
圖8a展示了模擬微細血流回波信號的奇異值分布情況，在Full-SVD壁濾波器中選取（4，?368）為轉折點進行濾波。圖8b展示了多項式回歸壁濾波器、Full-SVD、RS-RSVD對微細血流信號提取效果，根據實驗結果，在仿真實驗中三種濾波器算法均能有效地將血流信號提取出來。
圖8
Full-SVD算法分解仿真數據的奇異值分布及三種濾波算法提取血流信號強度對比
a. 仿真實驗奇異值排布；b. 多項式回歸壁濾波器、Full-SVD、RS-RSVD濾波效果對比
Figure8. The singular value distribution of simulation data decomposed by Full-SVD algorithm and comparison of the intensity of blood flow signals extracted by the three filtering algorithmsa. singular value arrangements of simulation experiment; b. comparison of filtering effects of polynomial regression wall filter, Full-SVD and RS-RSVD
2.2 實體數據成像實驗
實體數據成像實驗采用人體腎臟部分實體回波數據來驗證不同濾波器算法的濾波效果。實驗探頭選用SL14-3H探頭類型，中心頻率8.17 MHz，陣元數量192，聲速1 540 m/s，采樣頻率30 MHz，線密度0.5，感興趣區（region of interest，ROI）為13.96 mm*12.80 mm，連續時間序列為41，數據大小為272*128*41。
實體數據的連續時間序列是41，實驗設定勒讓德多項式大小為41*41，多項式回歸壁濾波器得到的實際濾波效果如圖9所示。結果表明，多項式回歸濾波器選取的階數適用于仿體實驗，但并不適用于實體數據濾波處理，并且由于多項式回歸濾波器的頻率響應關于零頻率軸對稱，需求的多項式階數較大，所以在濾除雜波信號時，錯誤地將血流信號濾除，結果無法有效提取微細血流信號。
圖9
多項式回歸壁濾波器處理的功率多普勒血流圖像
Figure9.
Power Doppler blood flow image processed by polyno mial regression wall filter
圖10a為Full-SVD壁濾波器算法根據實體數據得到的奇異值分布情況，根據奇異值排布曲線的轉折點來確定組織雜波截止點A。假設數據中存在高斯噪聲分布，噪聲對應的高階奇異值在對數下呈線性分布。對超過15階奇異值曲線進行線性擬合，可以得到奇異值曲線偏離擬合曲線的轉折點B。確定A、B兩點即為組織和噪聲在奇異值映射范圍內的截止階數。根據截止點A、B進行雜波抑制，得到如圖10b所示的部分腎臟微細血流功率多普勒圖像。實驗結果表明，相比多項式回歸型濾波器算法，Full-SVD算法可以有效提高微細血流信號的提取。
圖10
Full-SVD算法分解實體數據的奇異值分布和Full-SVD處理完成的功率多普勒血流圖像
a. 腎臟回波數據奇異值分布；b. Full-SVD處理完成的功率多普勒血流圖像
Figure10. The singular value distribution of entity data decomposed by Full-SVD algorithm and the power Doppler flow image completed by Full-SVD processinga. singular value distribution of kidney echo data; b. power Doppler blood flow image processed by Full-SVD
為驗證RS-RSVD算法的成像效果，我們調用MATLAB工具箱中的rand函數，對原始矩陣272*128*41進行隨機抽取，得到四個大小為136*64*41的隨機數組，對四個數組分別進行RSVD濾波實驗，得到如圖11a所示的濾波效果。RS1-RSVD、RS2-RSVD、RS3-RSVD、RS4-RSVD分別為四次抽取進行RSVD的濾波效果。將四個136*64*41矩陣數進行組合并恢復至原始矩陣位置，得到原始圖像272*128*41的超聲回波數據濾波結果，此即為濾除組織雜波的血流信號。圖11b為恢復原始圖像的功率多普勒血流圖像。假設算法在并行超聲系統中，當ROI為固定區域時，抽取分塊的數量將由并行系統可同時執行的流水線決定，當系統可以同時進行16通道算法處理，則原始矩陣可被分為16個68*32*41大小的矩陣。
圖11
RS-RSVD算法分塊濾波和恢復的功率多普勒血流圖像
a. RS-RSVD算法分塊濾波的功率多普勒圖像；b. 恢復的功率多普勒血流圖像
Figure11. RS-RSVD algorithm block filtering and restored power Doppler blood flow imagea. power Doppler images processed by RS-RSVD algorithm block filtering; b. restored power Doppler blood flow imaging
除了對三種算法形成的功率多普勒圖像進行視覺觀察外，本實驗將三種算法的信噪比（signal to noise ratio，SNR）和噪聲比（contrast to noise ratio，CNR）進行比較，結果見圖12。圖12藍色區域為分別選取的部分ROI，x = 1.5～5.0 mm，y =2.0～5.5 mm。SNR和CNR的具體表達式如下：
圖12
Full-SVD和RS-RSVD區域選取
Figure12.
Full-SVD and RS-RSVD area selection
|
Sblood代表血流信號的平均強度（圖12綠色方框中白線部分），Nnoise代表背景噪聲信號的平均強度（圖12中黃色方框）。
|
Sbackground代表信號平均強度（圖12中綠色方框）。
表2給出了三種算法的運算時間、SNR和CNR，由于多項式回歸壁濾波器算法無法有效提取微細血流信號，本次實驗只對比了Full-SVD算法和RS-RSVD算法的SNR和CNR。運算時間包括數據讀取和矩陣運算。根據實驗結果，多項式回歸壁濾波器算法處理時間最短，Full-SVD所處理的時間最長達48.50 s左右，RS-RSVD在4次采樣時處理時間為19.03 s左右。RS-RSVD算法相比Full-SVD算法SNR基本不變，CNR有所下降。
表2
三種算法運行時間、SNR和CNR
Table2.
Running time, SNR and CNR of three algorithms
為評估RS-RSVD不同抽取次數對算法帶來的影響，分別進行了9次、16次、25次、36次抽樣處理實驗。不同采樣次數的執行時間、SNR和CNR如表3所示。根據實驗結果，在分塊次數達16次時，微血流信號的信噪比與Full-SVD算法接近，此時濾波的執行時間為（4.65 ± 0.50）s。
表3
RS-RSVD不同抽樣次數的運行時間、SNR、CNR
Table3.
Running time, SNR and CNR of RS-RSVD with different sampling times
3 討論
本文通過仿真數據和實際人體腎臟超聲回波數據進行算法的對比實驗。實驗結果表明，傳統的多項式回歸壁濾波器算法相比Full-SVD算法和RS-RSVD算法濾波效果較差，在提取微細血流信號時會出現錯誤情況。在高幀頻超聲成像系統中，回波信號的矩陣大小會隨著采樣頻率和ROI的增加而成倍增大，這使Full-SVD算法需要的資源空間更大，并且當矩陣數據量增加時奇異值分解所需的處理時間會呈倍數增加，使Full-SVD算法無法應用在實時的超聲檢測系統中，有一定的局限性。根據表3實驗結果，RS-RSVD算法相比多項式回歸濾波器算法具有可觀的視覺分辨率，與Full-SVD算法信噪比接近，在16次抽樣時雖然對比噪聲比降低了2.05 dB，但執行時間降低了90.52%。
4 結論
在超聲彩色血流成像中，壁濾波器的性能直接影響著血流信號的提取。臨床診斷中，微血流檢測具有重要的價值。近幾年，高性能的硬件被廣泛應用于超聲成像系統中，如圖形處理器（graphics processing unit，GPU）、現場可編程門陣列（field programmable gate array，FPGA）等。RS-RSVD算法為超聲并行處理提供了思路，并行處理減少成像時間將作為未來的研究重點。在超聲臨床醫學中，由于人體結構的差異和實際環境造成的影響，RS-RSVD算法還需在臨床上進行多次試驗，爭取早日應用于高幀頻超聲微血管成像系統中。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：王寶宇、張石提出研究思路，設計研究方案；張淼、劉瑞麟負責進行實驗；王寶宇、張淼負責采集和分析數據，同時負責論文起草及最終版本修訂。
倫理聲明：本研究實體數據由東軟醫療股份有限公司提供及授權。
