一種基于Casorati-奇異值分解的超快平面波超聲多普勒自適應時空域雜波抑制算法

徐依雯 楊 晨 徐 杰 焦 陽 崔崤峣

①(蘇州國科昂卓醫療科技有限公司 蘇州 215163)

②(中國科學院蘇州生物醫學工程技術研究所 蘇州 215163)

③(中國科學技術大學電信學院 合肥 230031)

④(復旦大學工程與應用技術研究院 上海 200433)

1 引言

超聲多普勒是經典的血流分析工具，因無害性、實時性、高性價比等優勢得到廣泛運用，對血液動力學研究及臨床醫學均有非常重要的意義[1–3]。雜波抑制是超聲多普勒成像的重要一環，由靜止或慢動組織產生的雜波是圖像偽影的主要成因，該領域許多工作都圍繞更高效更精準的雜波濾波器展開[4–7]。

最經典的雜波抑制技術是高通濾波，如無限脈沖響應(Infinite Impulse Response,IIR)濾波器和有限脈沖響應(Finite Impulse Response,FIR)濾波器[4,6,7]。在早期的硬件水平限制下，2維超聲圖像經聚焦聲束逐列掃描拼接獲得，導致傳統超聲血流圖像視野狹小且時間分辨率低，基于組織與血流信號頻譜不重疊的假設，抑制雜波時通常只針對低頻組織信號(如圖1所示)[4]。然而研究表明，組織和血流信號的頻譜有重疊部分，尤其當兩者運動速度相近時(腦部與心臟處)，經典高通濾波器難以提取血流信號[8–11]。1996年Allam等人首提特征估計濾波器(Eigen-based filters)[12]，此后在特征估計基礎上發展出許多雜波抑制方法，如Hankel-奇異值分解法

圖1 理想高通濾波器示意圖

(Hankel Singular Value Decomposition,Hankel-SVD)，利用雜波強度通常比血流信號強度高40～100 dB的原理，通過舍棄信號主要成分來實現雜波抑制[12–15]，對不同輸入有特異性輸出，具有空間適應能力。此類方法需逐點構建慢時信息矩陣，對數據量和信號平穩性要求較高，不僅運算量大，而且在低速血流、高速組織運動或噪聲干擾較強的條件下濾波效果不佳。以上兩類經典濾波器，只利用超聲回波中單點的慢時維度信息，又稱時間域濾波器[11]。

復合平面波成像技術的發展，設備運算能力的升級，使高幀頻寬視野2維超聲信號同步采集成為現實[16]，平面波提供了豐富的2維空間信息[17–20]，是對慢時信息的重要補充，與多普勒相結合，對低速微血流非常敏感，因此平面波多普勒非常適用于觀察血流速度低、組織和血管分布錯綜復雜的腦部[21]。2015年Demené等人[11]提出了經典Casorati奇異值分解(Casorati Singular Value Decomposition,Casorati–SVD)濾波算法，結合了超快復合平面波提供的空間信息和特征提取算法的空間適應性優勢，極大地提高了超聲多普勒對微血流的靈敏度，為許多新型應用研究提供了技術支持[22–25]。該算法的缺點在于濾波截止參數為固定經驗值，改變成像目標后調整參數優化圖像的工作量較大。2017年Song等人[26]提出一種基于Casorati–SVD的分區域自適應局部雜波濾波器(Block-wise Adaptive Local Clutter Filtering,BALCF)，分析奇異值變化趨勢設置閾值，該算法對不同圖像有較好的自適應性，缺點在于閾值選取標準單一，像素值易發生“跳變”導致圖像失真。2017年Arnal等人[27,28]也提出一些利用空間相關性的Casorati-SVD改進算法，通過遍歷所有濾波系數選出空間相關指數最高的組合進行成像，運算量較大，犧牲時間分辨率來換取圖像質量的提升。

本文在Casorati–SVD基礎上提出一種改進的自適應雜波抑制算法Casorati奇異值分解-頻率與能量分析法(Casorati-Singular Value Decomposition-Frequency&Energy,Casorati-SVD-F&E)，基于奇異向量頻譜和能量特征兩個標準自適應提取血流信號，并結合分區域濾波算法進一步優化抑噪效果。本文利用仿體、人體、動物上采得的不同回波數據驗證該方法的濾波效果，在不同的實驗中，本算法均對圖像信噪比有明顯的提升效果。

2 Casorati-SVD-F&E的原理與實現

2.1 經典Casorati-SVD壁濾波

SVD算法將其所作用的分解數據矩陣分解成有限個分量，并將分量從主到次排序，選擇合適分量重構則可提取目標分量所含信息[11]。經典的Hankle-SVD 方法，利用單點慢時信號構建Hankle數據矩陣，因此Hankle-SVD只提取了慢時信號所包含的時間域信息[14]，這也是傳統掃描模式的局限—2維空間數據不同步，難以利用空間信息[4]。

超快平面波技術提供的2維空間同步射頻(Radio Frequency,RF)數據包含了可用的空間信息[11]，從取一組RF數據中25個相鄰點的慢時信號，做低通和高通處理得到低頻組織信號和高頻血流信號，分別計算自相關性，如圖2所示。高頻血流信號相關系數普遍偏低，說明相鄰位置的血流信號幾乎沒有相關性；而低頻組織信號的相關系數普遍較高，說明相鄰位置的組織信號具有顯著的相似性，圖2表明了空間信息特性。

圖2 相鄰25個像素點慢時信號相關矩陣

Casorati數據矩陣是一種同時利用空間和時間信息的數據矩陣，適用于平面波技術，其構建方法如圖3所示，每幀同步空間數據變形成矩陣S的一列，因此矩陣S同時包含空間信息和時間信息。

圖3 Casorati數據矩陣重構

將原始數據變形成Casorati矩陣S后，對S進行奇異值分解[11]，將S分解成3個矩陣乘積的表示形式，分解過程如式(1)所示

其中，Δ為奇異值矩陣，U為空間奇異向量矩陣，V為時間奇異向量矩陣。根據SVD特性及矩陣變換理論，S又可表示為加權信號分量之和，如式(2)所示，每個信號分量包含權重即奇異值λi，空間奇異向量Ui，以及調制Ui的時間奇異向量V i

奇異值標志信號分量能量大小，階數i越小奇異值越大，因此高強度組織分量集中于低階奇異向量，低強度噪聲分量集中于高階奇異向量，如圖4所示，利用這一分布特性，設計帶通濾波矩陣If可實現對血流信號的提取，如式(3)所示。

圖4 SVD奇異值分布示意圖

2.2 Casorati-SVD-F&E自適應壁濾波

良好的成像效果需要合適的濾波系數配合，為了更好地抑制雜波，須準確分離屬于雜波的信號分量，本文在Casorati-SVD的基礎上，提出一種自適應分離組織雜波和高頻噪聲的算法，分別從平均多普勒頻率和信號分量能量兩個角度出發，計算帶通濾波矩陣的截止波濾參數。

2.2.1頻譜分析法設計低階閾值抑制組織雜波

組織信號具有強度高而時變頻率低的特點，利用多普勒頻移的不同區分組織信號與血流信號十分有效。對一組數據做Casorati-SVD分解后，得到奇異值降序排列的時間奇異向量，主要包含原始數據中的慢時信息，對時間奇異向量做譜分析如圖5所示，隨著階數i增大，多普勒頻率向高頻移動，因此分析多普勒頻率計算閾值，能有效抑制低頻組織雜波且保證該閾值具有很強的普適性，平均多普勒頻率計算方法為

圖5 時間奇異向量頻譜圖

2.2.2能量分析法設計高階閾值分離高頻噪聲

高頻噪聲常與血流信號混合，區別在于不同分量中兩者占比不同。在有效探測區域內，血流信號的能量應強于噪聲信號，因此，以奇異值均值為閾值，奇異值高于閾值的分量視為血流分量，奇異值低于閾值的視為噪聲分量。以此方法選出的高閾值，在大部分圖像中都能較好地抑制高頻噪聲，但是在深度深、噪聲大的區域，噪聲占比明顯提高，此時需要依據噪聲強度對閾值進行校正：當高閾值過高時，從當前高閾值向低階檢索，直至兩閾值之差能滿足預設的要求，閾值差可根據成像目標調整，這種高閾值校正方法在實際應用中更加合理。

2.2.3分區域滑動濾波算法

線陣平面波所成圖像視野較大，淺表和深層區域噪聲水平不同，如全圖統一濾波，則對不同噪聲無法取得最佳抑制效果，而分區域濾波方法在處理局部噪聲方面優勢明顯。因此采用分區域算法，對不同深度不同噪聲水平的區域自適應濾波，取得更精準的抑噪效果。分區域濾波方法如圖6所示。

圖6 分區域濾波方法示意圖

權衡濾波效果和整體運算量，選擇合適的濾波區域參數(分區大小和重疊率)，噪聲抑制效果會有明顯提升。

2.2.4濾波效果評估方法

3 實驗系統及測量結果

為驗證SVD-F&E自適應壁濾波算法的雜波抑制效果，利用超快速平面波發射接收平臺美國Verasonics公司的Le64數據采集系統，分別在仿體、手臂動脈和家兔腦部采集回波數據，其他數據處理過程一致，選用不同的雜波抑制算法：經典Casorati-SVD固定經驗閾值算法(簡稱固定閾值法，fixed thresholds)、BALCF算法和本文的SVD-F&E算法進行濾波處理，將3種算法的結果進行對比。

3.1 仿體實驗

仿體實驗結果如圖7所示。圖7(d)是仿體實驗設置示意圖，水槽中傾斜放置一根外徑3 mm，內徑2 mm的管子，以注射泵推入抗凝牛血，采用128陣元線陣探頭L11-5v，中心頻率7.8 MHz，圖像幀率250 Hz采集一組2200幀平面波回波數據，經過波束合成、11角度復合、IQ解調和壁濾波處理，得到功率多普勒圖像。

圖7 仿體實驗與功率多普勒圖像

3.2 人體手臂動脈實驗

以手臂尺動脈為實驗對象，采用128陣元線陣探頭(vermon)如圖8(e)所示，陣列平行于尺動脈延伸方向放置，以中心頻率15 MHz，圖像幀率500 Hz采集一組1750幀多角度平面波回波數據，經過波束合成、7角度復合、IQ解調和壁濾波處理，得到功率多普勒圖像。

圖8 手臂動脈實驗與功率多普勒圖像

3.3 家兔大腦實驗

以成年家兔大腦為實驗對象，垂直兔腦中動脈放置探頭如圖9(d)所示，對兔腦冠狀面進行成像。采用線陣128陣元線陣探頭(vermon)，中心頻率15 MHz，圖像幀率800 Hz采集一組2100幀多角度平面波回波信號，經過波束合成、7角度復合、IQ解調和壁濾波處理，得到功率多普勒圖像。

圖9 家兔大腦功率多普勒圖像

4 討論

將BALCF和固定閾值法，與Casorati-SVDF&E的兩種用法(全局濾波、分區域濾波)進行對比，做信噪比曲線圖如圖10所示，圖10(a)—圖10(c)分別描述仿體實驗、人體手臂實驗、兔腦實驗中4種方法所處理圖像的SNR和CNR。從CNR角度來看，SVD-F&E的全局濾波和分區域濾波相比BALCF和固定閾值兩種方法均有不同程度的提升；從SNR角度來看，在仿體和手臂實驗中，全局SVDF&E比前兩種方法略有不如，原因是全局濾波時，為抑制血管周邊的局部高亮組織噪聲，犧牲了部分血流信號，只考慮有效信號和噪聲的信噪比SNR會略有下降，但因為抑制了組織噪聲，對比信噪比CNR有所提升。通過分析第3節的3組實驗，基于Casorati-SVD的分區域自適應-F&E雜波抑制算法，通過針對局部噪聲區域采用最優化自適應參數濾波，減小了局部噪聲，大幅提高了SNR和CNR。

圖10 信噪比曲線

分區域滑動濾波能大幅提高圖像信噪比，也因計算量增大影響計算速度，因此不同應用場景需要選擇合適的濾波區域參數，本文算法還可以通過調整單次濾波區域大小和滑動重疊比來滿足不同的成像需求(成像速度和圖像質量)：單次濾波區域越大，計算量越小，成像越快，當單次濾波區域為圖像全局時，成像速度很快，能滿足實時成像的需求；單次濾波區域越小，重疊比越高，則圖像質量提高，但相應的成像速度會下降。該算法有兩個注意點：一是單次濾波區域不能太小，nx·n z需大于慢時方向幀數nt，若區域太小，SVD分解時數據量過小，導致SVD分解運算結果極不準確；二是濾波區域滑動重疊比一般不得小于50%，否則相鄰區域濾波參數相差過大，導致鄰近像素值階躍，最終圖像中出現“斷層”現象。

本文提出的Casorati-SVD-F&E算法，通過多普勒頻率和信號能量兩個標準確定截止濾波參數，既具有對不同目標自適應調節的優點，又能根據經驗修正濾波參數分析算法，提高了雜波抑制算法的普適性。算法不足之處有二，一是為抑制局部噪聲所采用的分區域濾波極大地增加了運算代價，降低時間分辨率；二是分區域濾波要優化圖像，濾波區域大小、重疊比等參數目前依賴經驗設置，自適應濾波參數優化的工作量較大。下一步研究方向為濾波區域參數的優化，針對不同圖像能更快地確定最優濾波區域參數，取得更佳濾波效果。

5 結論

超快平面波技術發展迅速且應用前景廣泛，基于平面波的多普勒血流成像也有很強的臨床應用需求，面對成像速度和成像質量兩方面的要求，利用Casorati數據矩陣實現的雜波抑制算法具有速度快、濾波效果好等優勢。本文提出一種改進的自適應壁濾波算法Casorati-SVD-F&E，更充分地利用了平面波所提供的2維空間信息，并在不同目標(仿體、手臂動脈、家兔大腦)上驗證了該算法抑制噪聲提高信噪比的能力，與傳統壁濾波算法相比，具有空間自適應性，能根據信號特征自動濾波，輔以經驗參數校正法，可更高效地進行雜波抑制，提高多普勒血流圖像的質量。