有效抑制由靜態或慢動組織引起的雜波是診斷超聲彩色血流成像技術研究的關鍵問題之一,不充分的雜波抑制將導致多普勒平均頻率估計趨同于雜波成分特性,從而導致潛在的血流信息描述誤差。本文在研究廣義壁濾波器設計原理基礎之上,比較了三種不同類型壁濾波器——投影初始化有限沖擊響應濾波器(Prj-ⅡR)、多項式回歸濾波器(Pol-Reg)與基于特征值分解(Eigen-based)濾波器的設計方法差異,并對其雜波抑制性能進行了仿真對比分析。血流平均頻率一維自相關估計結果表明,投影初始化ⅡR壁濾波器與多項式回歸壁濾波器雜波抑制性能相近,在平穩雜波環境下可實現血流平均速度的精確估計,當多普勒采樣向量包含樣本數據增加時,改變設計參數可進一步改善雜波抑制能力,而基于特征值分解的壁濾波器則可對非平穩雜波進行抑制,進一步改善低速血流的估計精度,而且在多普勒矢量包含樣本數小于10的情況下,計算復雜度并無顯著增加。
引用本文： 王錄濤, 肖軍, 柴華. 超聲血流成像中的壁濾波器設計方法分析與性能比較. 生物醫學工程學雜志, 2015, 32(4): 773-778. doi: 10.7507/1001-5515.20150141 復制
引言
在彩色超聲血流成像技術中，由靜止或慢動組織反射引起的雜波信號強度通常比血流信號強度高40~100 dB，在多普勒回波信號中處于主導地位。即使當采樣容積位于血管中心位置處，由于采樣容積的三維特性、合成聲束旁瓣泄露、聲波多重反射以及距離模糊等因素的影響，血流回波信號中也包含著大量的雜波信號[1-2]。雜波信號的存在導致血流參數估計產生偏差，從而嚴重降低了系統的血流鑒別能力，因此充分抑制雜波成分對于血流信號特征的準確提取至關重要。
根據雜波與血流信號頻域特性的差別，一般采用壁濾波器對雜波成分進行抑制。常用的壁濾波器是一種高通濾波器，包括有限沖擊響應(finite impulse response,FIR)濾波器、無限沖擊響應(infinite impulse response,IIR)濾波器等。在彩色超聲血流成像系統中，受制于組織與血流的非平穩性，一般僅有8~16個回波采樣數據可用于壁濾波處理，有限的采樣數據導致FIR濾波器的幅頻特性無法滿足高強度雜波抑制需求。與FIR濾波器相比，相同階數的IIR濾波器具有更加陡峭的過渡帶滾降特性與阻帶抑制特性，但當采樣數據有限時，IIR濾波器的幅頻響應由暫態響應主導。為增強阻帶抑制能力與減小過渡帶帶寬，一般采用階躍初始化、投影初始化等方法可消除暫態響應的影響[3]。多項式回歸濾波器(polynomials regression,Pol-Reg)采用正交多項式對雜波分量進行最小二乘擬合來消除雜波成分的影響，為壁濾波器的設計提供了新的途徑[4]。針對雜波信號的非平穩特性，Kruse等[5]提出了基于特征值分解(eigen-based)壁濾波器設計方法，根據采樣數據的二階統計特性動態改變壁濾波器的截止頻率，實現動態雜波干擾下低速血流信號的提取[6]。
本文從廣義壁濾波設計的基本原理出發，深入分析三種不同類型的壁濾波設計原理與設計參數對頻域響應的影響，最后采用一維自相關估計法對兩種不同雜波情況下的三種濾波方法的血流平均頻率估計偏差與標準差進行量化評估，并對計算復雜度進行分析。
1 廣義壁濾波器
彩色超聲血流成像利用超聲聲束對待成像區域進行掃描，在一系列離散掃描方向形成有限個發射聲束，通過接收血紅細胞散射的多普勒信號估計血流動力學參數[3]。在深度為k處獲得的解調后的采樣復多普勒信號向量可表示為
| $\boldsymbol{x} = {\left[ {x\left( 0 \right),x\left( 1 \right), \cdots ,x\left( {N - 1} \right)} \right]^T}$ |
其中N為包含的數據樣本數。
為避免雜波多普勒信號與血流多普勒信號的相互調制，壁濾波器一般是線性的。因此，廣義壁濾波器可看作N維復向量空間CN上的線性變換，可用簡單的矩陣乘法表示[4]：
| $\boldsymbol{y} = \boldsymbol{A}\boldsymbol{x}$ |
其中A為M×N維系數矩陣，濾波器輸出向量y=[y(0),y(1),…,y(M-1)]T為M×1維矩陣。設A中第n行第k列元素為a(n,k)，則濾波器輸出向量元素y(n)為
| $y\left( n \right) = \sum\limits_{k = 0}^{N - 1} {a\left( {n,k} \right)x\left( k \right),n = 0, \cdots ,M - 1} $ |
廣義濾波器是線性的但不是時不變的，因此不能利用濾波器的單位沖擊響應的傅里葉變換描述廣義濾波器的頻域特性，但可采用輸入為單位復諧波信號時的輸出功率對廣義濾波器的頻域特性進行描述。當輸入為單位復諧波信號時，廣義濾波器輸出功率為[4]
| ${y_\omega }\left( n \right) = \sum\limits_{k = 0}^{N - 1} {a\left( {n,k} \right){e^{jk\omega }} = {A_n}\left( { - \omega } \right),n = 0, \cdots ,M - 1} $ |
其中An(ω)為A的第n行元素的傅里葉變換。廣義濾波器的頻率響應可表示為
| ${H_0}\left( \omega \right) = \frac{1}{M}\sum\limits_{n = 0}^{M - 1} {{{\left| {{y_\omega }\left( n \right)} \right|}^2} = \frac{1}{M}\sum\limits_{n = 0}^{M - 1} {{{\left| {{A_n}\left( { - \omega } \right)} \right|}^2}} } $ |
2 壁濾波器設計方法分析
2.1 投影初始化有限沖擊響應(infinitely impulse response with projection initialization,Prj-IIR)濾波器
K階IIR濾波器可用下述差分方程進行描述[4]：
| $y\left( n \right) = - \sum\limits_{k = 1}^K {{a_k}y\left( {n - k} \right) + \sum\limits_{k = 0}^K {{b_k}x\left( {n - k} \right)} } $ |
在僅有有限個采樣數據情況下，IIR濾波器的性能由暫態響應主導。為分析不同初始化方法對IIR濾波器暫態響應的影響，可將IIR濾波器用狀態方程進行表示。選取狀態向量v(n)=[v1(n),v2(n),…,vk(n)]T，則IIR濾波的狀態空間描述為：
| $\begin{array}{l} \boldsymbol{v}\left( {n + 1} \right) = \boldsymbol{F}\boldsymbol{v}\left( n \right) + \boldsymbol{q}\left( n \right)\\ y\left( n \right) = {\boldsymbol{g}^T}\boldsymbol{v}\left( d \right) + dx\left( n \right) \end{array}$ |
其中
| $\begin{array}{l} \boldsymbol{F} = \left[ {\begin{array}{*{20}{c}} 0&1&0& \cdots &0\\ 0&0&1& \cdots &0\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0&0&0& \cdots &1\\ { - {a_K}}&{ - {a_{K - 1}}}&{ - {a_{K - 2}}}& \cdots &{ - {a_1}} \end{array}} \right],\boldsymbol{q} = \left[ {\begin{array}{*{20}{c}} 0\\ 0\\ \vdots \\ 0\\ 1 \end{array}} \right]\\ \boldsymbol{g} = \left[ {\begin{array}{*{20}{c}} {{b_K} - {b_0}{a_K}}\\ {{b_{K - 1}} - {b_0}{a_{K - 1}}}\\ {{b_1} - {b_0}{a_1}} \end{array}} \right],d = {b_0} \end{array}$ |
當輸入與輸出均為N×1維向量時，IIR濾波輸出與輸入關系可用矩陣向量方程表示：
| $\boldsymbol{y} = \boldsymbol{B}\boldsymbol{v}\left( 0 \right) + \boldsymbol{C}\boldsymbol{x}$ |
其中
| $\boldsymbol{B} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{g}^T}}\\ {{\boldsymbol{g}^T}\boldsymbol{F}}\\ {{\boldsymbol{g}^T}{\boldsymbol{F}^{N - 1}}} \end{array}} \right]\boldsymbol{C} = \left[ {\begin{array}{*{20}{c}} d&0& \cdots &0\\ {{\boldsymbol{g}^T}\boldsymbol{q}}&d& \cdots &0\\ \vdots & \vdots & \vdots & \vdots \\ {{\boldsymbol{g}^T}{\boldsymbol{F}^{N - 2}}\boldsymbol{q}}&{{\boldsymbol{g}^T}{\boldsymbol{F}^{N - 3}}\boldsymbol{q}}& \cdots &d \end{array}} \right]$ |
在式(9)中，通過選擇合適的初始化向量v(n)，可以有效消除暫態效應，使濾波器的幅頻響應更趨近于穩態響應。常用的初始化方法有零初始、階躍初始化與投影初始化三種。其中投影初始化方法認為暫態響應位于由矩陣B的列向量張成的子空間內，通過約束輸出信號在暫態子空間的投影為零，即PBy=0，PB=B(BTB)-1BT為投影矩陣，可得初始化向量：
| $v\left( 0 \right) = - \boldsymbol{B}{\left( {{\boldsymbol{B}^T}\boldsymbol{B}} \right)^{ - 1}}{\boldsymbol{B}^T}\boldsymbol{C}\boldsymbol{x}$ |
將式(11)代入式(9)可得Prj-IIR濾波器系統傳遞函數為：
| $\boldsymbol{y} = \left( {\boldsymbol{I} - \boldsymbol{B}{{\left( {{\boldsymbol{B}^T}\boldsymbol{B}} \right)}^{ - 1}}{\boldsymbol{B}^T}} \right)\boldsymbol{C}\boldsymbol{x}$ |
令A=(I-B(BTB)-1BT)C，根據廣義濾波器頻率響應定義可求得Prj-IIR濾波器頻域響應。與零初始化方法與階躍初始化方法相比，Prj-IIR通過去除在由B的列向量張成的暫態子空間內的輸出信號成分，縮小過渡帶寬，同時增加阻帶衰減。
2.2 Pol-Reg濾波器
Pol-Reg濾波器利用一組給定階數的正交多項式基函數在時域對雜波信號進行最小二乘擬合，然后通過從超聲多普勒回波信號中減去擬合的雜波成分以獲得期望的血流信號[4]。由多項式基函數張成的雜波空間κ是輸入復信號空間CN的子空間，從CN到κ的投影變換Pκ則給出了雜波成分的最小二乘估計，多項式回歸濾波器系數矩陣可表示為[4]
| $\boldsymbol{A} = \boldsymbol{I} - {P_\kappa } = \boldsymbol{I} - \sum\limits_{k = 0}^{K - 1} {{\boldsymbol{b}_k}\boldsymbol{b}_k^H} $ |
其中bk為張成K維雜波空間的正交基向量，I為單位矩陣。根據式(5)，多項式回歸濾波器頻域響應可表示為
| ${H_0}\left( \omega \right) = 1 - \frac{1}{N}\sum\limits_{k = 0}^{K - 1} {{{\left| {{B_k}\left( \omega \right)} \right|}^2}} $ |
其中
| ${B_k}\left( \omega \right) = \sum\limits_{n = 1}^N {{b_k}\left( n \right){e^{ - in\omega }}} $ |
2.3 Eigen-based濾波器
Eigen-based濾波器采用Karhunen-Loeve變換將解調后的采樣復多普勒信號向量x展開為一組正交基函數的加權和[5]：
| $\boldsymbol{x} = \sum\limits_{k = 1}^N {{\gamma _k}{\boldsymbol{e}_k}} $ |
其中，ek是第k個特征向量，γk是第k個展開權值，滿足
| $E\left\{ {{\gamma _k}{\gamma _l}} \right\} = \left\{ \begin{array}{l} {\lambda _k},k = 1\\ 0,k \ne l \end{array} \right.$ |
其中λk是與ek相對應的第k個特征值。將λk按降序排列，根據雜波向量的時域與頻域特性自適應選取與前Kc個特征值相對應的特征向量張成雜波子空間фc：
| $\widehat {\bf{c}} = \sum\limits_{k = 1}^{{K_c}} {{\gamma _k}{{\bf{e}}_k} \Leftrightarrow {\bf{C}} = \left[ {\begin{array}{*{20}{c}} |&|&{}&|\\ {{{\bf{e}}_1}}&{{{\bf{e}}_2}}& \cdots &{{{\bf{e}}_{{K_c}}}}\\ |&|&{}&| \end{array}} \right]} $ |
則Eigen-based濾波器輸出可由фc的正交補空間與采樣復多普勒信號向量的矩陣乘積給出：
| $\boldsymbol{y} = \left( {\boldsymbol{I} - \sum\limits_{k = 1}^{{K_c}} {{\boldsymbol{e}_k}\boldsymbol{e}_k^H} } \right)\boldsymbol{x}$ |
由式(5)、(19)可知，Eigen-based濾波器的響應由采樣復多普勒信號向量分解得到的雜波特征向量決定。當組織運動狀態發生變化引起雜波成分發生改變時，雜波特征向量的維數Kc會相應改變，導致Eigen-based壁濾波器雜波截至特性產生相應變化，因而具有自適應性。
對采樣復多普勒信號向量的協方差矩陣進行特征值分解，可獲得式(16)中的特征向量與特征值：
| ${\boldsymbol{R}_x} = E\left\{ {x{x^H}} \right\} = \sum\limits_{k = 1}^N {{\lambda _k}{\boldsymbol{e}_k}\boldsymbol{e}_k^H} $ |
在實際應用中，Rx由沿掃描深度方向相鄰L個采樣復多普勒信號向量估計得出：
| ${\boldsymbol{R}_x} \approx \frac{1}{L}\sum\limits_{l = 1}^L {{\boldsymbol{x}_l}\boldsymbol{x}_l^H} $ |
其中，L需滿足平穩性假設。
3 仿真分析
3.1 雜波抑制性能比較
為了對三種壁濾波器的雜波抑制性能進行量化評估，采用Yu[6]提出的回波合成模型產生多普勒回波信號。仿真共產生兩組多普勒回波信號，每組包含的雜波分量的帶寬不同，以對不同雜波環境下的各種壁濾波器的性能進行比較。壁濾波處理后的采樣數據采用標準一維自相關估計法[7]估計血流平均多普勒頻率。為分析壁濾波技術對不同血流流速估計的影響，設計32組歸一化平均多普勒頻率在0~0.5內均勻分布的血流信號。仿真參數如表 1所示。
表1
仿真參數設置
Table1.
Parameters used for the simulation
圖 1給出了雜波1、采樣多普勒向量信號長度N為8時，經1 000次獨立仿真實驗得到的血流平均多普勒頻率估計偏差與標準差。其中Prj-IIR為投影初始化Chebyshev濾波器，通帶歸一化截止頻率為0.1，階數為3；Pol-Reg的正交基采用Legendre多項式基函數，階數為3，以確保阻帶截止頻率與Prj-IIR相近。對于Eigen-based壁濾波器，按照Yu[8-9]提出的基于頻率比較方法自適應選取雜波空間維數。從圖 1可以看出，三種壁濾波器均有效抑制了位于阻帶內的多普勒頻率成分，導致在該頻率范圍內的頻率估計值與理論值相比產生顯著的差異，而在通帶范圍內，雜波頻率成分得以充分抑制，血流頻率成分得以保留，隨著血流多普勒頻率的增加，頻率估計產生的偏差與標準差則迅速降低并趨近于零。對比經三種壁濾波器處理后平均頻率估計結果可知，Prj-IIR與Pol-Reg性能相近，而Eigen-based壁濾波器在抑制雜波的同時可盡量保證血流成分不被衰減，因而可在較低血流多普勒頻率(250 Hz)時達到穩態估計結果。
圖1
血流平均頻率估計(雜波1，N=8)
Figure1.
Blood mean frequency estimation with clutter signal 1 (package size N=8)
為比較多普勒向量信號采樣長度對壁濾波器性能的影響，圖 2給出了將N從8增加到16的估計結果。為保證近似阻帶帶寬，Prj-IIR濾波器的階數增加為4，通帶歸一化截止頻率為0.05；Pol-Reg階數為4階。比較圖 1與圖 2可以發現，三種不同類型壁濾波器均可對1型強度、中心頻率與帶寬的雜波信號進行有效抑制，通過增加多普勒向量信號采樣長度可以顯著降低過渡帶寬，增強阻帶抑制能力。
圖2
血流平均頻率估計(雜波1，N=16)
Figure2.
Blood mean frequency estimation with clutter signal 1 (package size N=16)
圖 3給出了雜波2、采樣多普勒向量信號長度N為8時的估計結果。濾波器參數設計同圖 1。從圖 3中可以看出2型雜波嚴重影響了平均頻率估計結果。由于雜波帶寬為1型雜波的2.5倍，在不改變設計參數的情況下，Prj-IIR與Pol-Reg濾波器無法對雜波進行有效抑制，導致部分低頻雜波成分混入血流多普勒信號，造成血流平均頻率估計低于理論值，因而在高頻時出現顯著的負估計偏差，估計標準差也明顯高于圖 1中所示的針對雜波1的估計結果。而Eigen-based壁濾波器根據輸入信號的二階統計特性自適應選擇雜波空間維數，阻帶帶寬與圖 1相比顯著增加以確保雜波成分的有效抑制，因而可實現通帶內平均多普勒頻率的準確估計。
圖3
血流平均頻率估計(雜波2，N=8)
Figure3.
Blood mean frequency estimation with clutter signal 2 (package size N=8)
圖 4為將N從8增加到16時，針對雜波2的平均頻率估計結果。Prj-IIR與Pol-Reg設計參數同圖 2。對比圖 3與圖 4可知，增加多普勒向量信號采樣長度可以改善Prj-IIR與Pol-Reg的雜波抑制性能，減小平均頻率估計誤差，但與圖 1和圖 2的處理結果相比，通帶標準差仍存在較大差異，因而無法有效抑制雜波對平均頻率估計的影響。對于Eigen-based壁濾波器，N的增加使得過渡帶變窄，產生較大平均頻率，估計偏差的頻率范圍進一步減小。
圖4
血流平均頻率估計(雜波2，N=16)
Figure4.
Blood mean frequency estimation with clutter signal 2 (package size N=16)
3.2 計算復雜度比較
Prj-IIR、Pol-Reg為非自適應濾波器，濾波操作為矩陣乘法運算，因此，系數矩陣A的維數決定了算法的復雜度。Eigen-based壁濾波器的計算包括數據樣本的協方差矩陣估計、特征值分解與濾波運算三部分，完成濾波操作所需的浮點運算次數(float-point operations per second,FLOPS)由多普勒向量維數N、軸向采樣體積數、用以估計數據樣本的協方差矩陣的多普勒向量數L與雜波空間維數Kc共同決定。圖 5示出了三種濾波算法計算量比較，軸向采樣體積數為150。從圖 5可知，Prj-IIR與Pol-Reg計算量相同，當N≤10時，Eigen-based的計算量與Prj-IIR、Pol-Reg相比有所增加，但隨著N的增加，計算所需FLOPS則顯著增加。Kc對Eigen-based型壁濾波器的計算量影響不明顯，當Kc從2增加到5時，計算量略有增加，當N≤10時，增加量基本可以忽略。
圖5
浮點計算量比較
Figure5.
Flop counts comparison
4 結論
雜波抑制是彩色超聲血流信號處理的一個重要環節，壁濾波的雜波抑制性能直接影響血流參數估計精度。Prj-IIR、Pol-Reg與Eigen-based壁濾波器均通過計算輸入多普勒向量在雜波信號基向量上的投影，并從輸入信號中減去該投影值，實現雜波抑制。在阻帶帶寬相同時，Prj-IIR與Pol-Reg雜波抑制性能相近。對于Prj-IIR濾波器，增加階數的同時保持低的通帶截止頻率可增加阻帶帶寬，而Pol-Reg壁濾波器的阻帶帶寬僅取決于設計階數，Eigen-based壁濾波器根據雜波統計特性動態改變雜波截止特性。針對靜態或均勻運動組織，三種壁濾波技術均可有效抑制雜波，而增加采樣多普勒向量樣本有助于進一步增強阻帶抑制能力與降低過渡帶寬。對于組織空時非平穩運動，Eigen-based壁濾波技術通過動態改變雜波截止特性可實現雜波成分的充分抑制，而對于非自適應濾波器，由于阻帶截止頻率由組織運動特性的先驗知識確定，當組織運動特性發生變化導致雜波帶寬改變時，雜波成分的泄漏或血流成分的過渡抑制都將導致血流參數估計相對實際值產生較大偏差。Prj-IIR與Pol-Reg的計算復雜度為O(N2)，Eigen-based壁濾波器的計算復雜度為O(N3)，但當多普勒向量維數N≤10時，兩者計算所需FLOPS相差不顯著。雖然Eigen-based可有效提高血流信號的檢測能力，但在臨床診斷應用中，采樣多普勒向量協方差矩陣的構建、雜波子空間維數的準確選取等問題仍亟待解決。
引言
在彩色超聲血流成像技術中，由靜止或慢動組織反射引起的雜波信號強度通常比血流信號強度高40~100 dB，在多普勒回波信號中處于主導地位。即使當采樣容積位于血管中心位置處，由于采樣容積的三維特性、合成聲束旁瓣泄露、聲波多重反射以及距離模糊等因素的影響，血流回波信號中也包含著大量的雜波信號[1-2]。雜波信號的存在導致血流參數估計產生偏差，從而嚴重降低了系統的血流鑒別能力，因此充分抑制雜波成分對于血流信號特征的準確提取至關重要。
根據雜波與血流信號頻域特性的差別，一般采用壁濾波器對雜波成分進行抑制。常用的壁濾波器是一種高通濾波器，包括有限沖擊響應(finite impulse response,FIR)濾波器、無限沖擊響應(infinite impulse response,IIR)濾波器等。在彩色超聲血流成像系統中，受制于組織與血流的非平穩性，一般僅有8~16個回波采樣數據可用于壁濾波處理，有限的采樣數據導致FIR濾波器的幅頻特性無法滿足高強度雜波抑制需求。與FIR濾波器相比，相同階數的IIR濾波器具有更加陡峭的過渡帶滾降特性與阻帶抑制特性，但當采樣數據有限時，IIR濾波器的幅頻響應由暫態響應主導。為增強阻帶抑制能力與減小過渡帶帶寬，一般采用階躍初始化、投影初始化等方法可消除暫態響應的影響[3]。多項式回歸濾波器(polynomials regression,Pol-Reg)采用正交多項式對雜波分量進行最小二乘擬合來消除雜波成分的影響，為壁濾波器的設計提供了新的途徑[4]。針對雜波信號的非平穩特性，Kruse等[5]提出了基于特征值分解(eigen-based)壁濾波器設計方法，根據采樣數據的二階統計特性動態改變壁濾波器的截止頻率，實現動態雜波干擾下低速血流信號的提取[6]。
本文從廣義壁濾波設計的基本原理出發，深入分析三種不同類型的壁濾波設計原理與設計參數對頻域響應的影響，最后采用一維自相關估計法對兩種不同雜波情況下的三種濾波方法的血流平均頻率估計偏差與標準差進行量化評估，并對計算復雜度進行分析。
1 廣義壁濾波器
彩色超聲血流成像利用超聲聲束對待成像區域進行掃描，在一系列離散掃描方向形成有限個發射聲束，通過接收血紅細胞散射的多普勒信號估計血流動力學參數[3]。在深度為k處獲得的解調后的采樣復多普勒信號向量可表示為
| $\boldsymbol{x} = {\left[ {x\left( 0 \right),x\left( 1 \right), \cdots ,x\left( {N - 1} \right)} \right]^T}$ |
其中N為包含的數據樣本數。
為避免雜波多普勒信號與血流多普勒信號的相互調制，壁濾波器一般是線性的。因此，廣義壁濾波器可看作N維復向量空間CN上的線性變換，可用簡單的矩陣乘法表示[4]：
| $\boldsymbol{y} = \boldsymbol{A}\boldsymbol{x}$ |
其中A為M×N維系數矩陣，濾波器輸出向量y=[y(0),y(1),…,y(M-1)]T為M×1維矩陣。設A中第n行第k列元素為a(n,k)，則濾波器輸出向量元素y(n)為
| $y\left( n \right) = \sum\limits_{k = 0}^{N - 1} {a\left( {n,k} \right)x\left( k \right),n = 0, \cdots ,M - 1} $ |
廣義濾波器是線性的但不是時不變的，因此不能利用濾波器的單位沖擊響應的傅里葉變換描述廣義濾波器的頻域特性，但可采用輸入為單位復諧波信號時的輸出功率對廣義濾波器的頻域特性進行描述。當輸入為單位復諧波信號時，廣義濾波器輸出功率為[4]
| ${y_\omega }\left( n \right) = \sum\limits_{k = 0}^{N - 1} {a\left( {n,k} \right){e^{jk\omega }} = {A_n}\left( { - \omega } \right),n = 0, \cdots ,M - 1} $ |
其中An(ω)為A的第n行元素的傅里葉變換。廣義濾波器的頻率響應可表示為
| ${H_0}\left( \omega \right) = \frac{1}{M}\sum\limits_{n = 0}^{M - 1} {{{\left| {{y_\omega }\left( n \right)} \right|}^2} = \frac{1}{M}\sum\limits_{n = 0}^{M - 1} {{{\left| {{A_n}\left( { - \omega } \right)} \right|}^2}} } $ |
2 壁濾波器設計方法分析
2.1 投影初始化有限沖擊響應(infinitely impulse response with projection initialization,Prj-IIR)濾波器
K階IIR濾波器可用下述差分方程進行描述[4]：
| $y\left( n \right) = - \sum\limits_{k = 1}^K {{a_k}y\left( {n - k} \right) + \sum\limits_{k = 0}^K {{b_k}x\left( {n - k} \right)} } $ |
在僅有有限個采樣數據情況下，IIR濾波器的性能由暫態響應主導。為分析不同初始化方法對IIR濾波器暫態響應的影響，可將IIR濾波器用狀態方程進行表示。選取狀態向量v(n)=[v1(n),v2(n),…,vk(n)]T，則IIR濾波的狀態空間描述為：
| $\begin{array}{l} \boldsymbol{v}\left( {n + 1} \right) = \boldsymbol{F}\boldsymbol{v}\left( n \right) + \boldsymbol{q}\left( n \right)\\ y\left( n \right) = {\boldsymbol{g}^T}\boldsymbol{v}\left( d \right) + dx\left( n \right) \end{array}$ |
其中
| $\begin{array}{l} \boldsymbol{F} = \left[ {\begin{array}{*{20}{c}} 0&1&0& \cdots &0\\ 0&0&1& \cdots &0\\ \vdots & \vdots & \vdots & \vdots & \vdots \\ 0&0&0& \cdots &1\\ { - {a_K}}&{ - {a_{K - 1}}}&{ - {a_{K - 2}}}& \cdots &{ - {a_1}} \end{array}} \right],\boldsymbol{q} = \left[ {\begin{array}{*{20}{c}} 0\\ 0\\ \vdots \\ 0\\ 1 \end{array}} \right]\\ \boldsymbol{g} = \left[ {\begin{array}{*{20}{c}} {{b_K} - {b_0}{a_K}}\\ {{b_{K - 1}} - {b_0}{a_{K - 1}}}\\ {{b_1} - {b_0}{a_1}} \end{array}} \right],d = {b_0} \end{array}$ |
當輸入與輸出均為N×1維向量時，IIR濾波輸出與輸入關系可用矩陣向量方程表示：
| $\boldsymbol{y} = \boldsymbol{B}\boldsymbol{v}\left( 0 \right) + \boldsymbol{C}\boldsymbol{x}$ |
其中
| $\boldsymbol{B} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{g}^T}}\\ {{\boldsymbol{g}^T}\boldsymbol{F}}\\ {{\boldsymbol{g}^T}{\boldsymbol{F}^{N - 1}}} \end{array}} \right]\boldsymbol{C} = \left[ {\begin{array}{*{20}{c}} d&0& \cdots &0\\ {{\boldsymbol{g}^T}\boldsymbol{q}}&d& \cdots &0\\ \vdots & \vdots & \vdots & \vdots \\ {{\boldsymbol{g}^T}{\boldsymbol{F}^{N - 2}}\boldsymbol{q}}&{{\boldsymbol{g}^T}{\boldsymbol{F}^{N - 3}}\boldsymbol{q}}& \cdots &d \end{array}} \right]$ |
在式(9)中，通過選擇合適的初始化向量v(n)，可以有效消除暫態效應，使濾波器的幅頻響應更趨近于穩態響應。常用的初始化方法有零初始、階躍初始化與投影初始化三種。其中投影初始化方法認為暫態響應位于由矩陣B的列向量張成的子空間內，通過約束輸出信號在暫態子空間的投影為零，即PBy=0，PB=B(BTB)-1BT為投影矩陣，可得初始化向量：
| $v\left( 0 \right) = - \boldsymbol{B}{\left( {{\boldsymbol{B}^T}\boldsymbol{B}} \right)^{ - 1}}{\boldsymbol{B}^T}\boldsymbol{C}\boldsymbol{x}$ |
將式(11)代入式(9)可得Prj-IIR濾波器系統傳遞函數為：
| $\boldsymbol{y} = \left( {\boldsymbol{I} - \boldsymbol{B}{{\left( {{\boldsymbol{B}^T}\boldsymbol{B}} \right)}^{ - 1}}{\boldsymbol{B}^T}} \right)\boldsymbol{C}\boldsymbol{x}$ |
令A=(I-B(BTB)-1BT)C，根據廣義濾波器頻率響應定義可求得Prj-IIR濾波器頻域響應。與零初始化方法與階躍初始化方法相比，Prj-IIR通過去除在由B的列向量張成的暫態子空間內的輸出信號成分，縮小過渡帶寬，同時增加阻帶衰減。
2.2 Pol-Reg濾波器
Pol-Reg濾波器利用一組給定階數的正交多項式基函數在時域對雜波信號進行最小二乘擬合，然后通過從超聲多普勒回波信號中減去擬合的雜波成分以獲得期望的血流信號[4]。由多項式基函數張成的雜波空間κ是輸入復信號空間CN的子空間，從CN到κ的投影變換Pκ則給出了雜波成分的最小二乘估計，多項式回歸濾波器系數矩陣可表示為[4]
| $\boldsymbol{A} = \boldsymbol{I} - {P_\kappa } = \boldsymbol{I} - \sum\limits_{k = 0}^{K - 1} {{\boldsymbol{b}_k}\boldsymbol{b}_k^H} $ |
其中bk為張成K維雜波空間的正交基向量，I為單位矩陣。根據式(5)，多項式回歸濾波器頻域響應可表示為
| ${H_0}\left( \omega \right) = 1 - \frac{1}{N}\sum\limits_{k = 0}^{K - 1} {{{\left| {{B_k}\left( \omega \right)} \right|}^2}} $ |
其中
| ${B_k}\left( \omega \right) = \sum\limits_{n = 1}^N {{b_k}\left( n \right){e^{ - in\omega }}} $ |
2.3 Eigen-based濾波器
Eigen-based濾波器采用Karhunen-Loeve變換將解調后的采樣復多普勒信號向量x展開為一組正交基函數的加權和[5]：
| $\boldsymbol{x} = \sum\limits_{k = 1}^N {{\gamma _k}{\boldsymbol{e}_k}} $ |
其中，ek是第k個特征向量，γk是第k個展開權值，滿足
| $E\left\{ {{\gamma _k}{\gamma _l}} \right\} = \left\{ \begin{array}{l} {\lambda _k},k = 1\\ 0,k \ne l \end{array} \right.$ |
其中λk是與ek相對應的第k個特征值。將λk按降序排列，根據雜波向量的時域與頻域特性自適應選取與前Kc個特征值相對應的特征向量張成雜波子空間фc：
| $\widehat {\bf{c}} = \sum\limits_{k = 1}^{{K_c}} {{\gamma _k}{{\bf{e}}_k} \Leftrightarrow {\bf{C}} = \left[ {\begin{array}{*{20}{c}} |&|&{}&|\\ {{{\bf{e}}_1}}&{{{\bf{e}}_2}}& \cdots &{{{\bf{e}}_{{K_c}}}}\\ |&|&{}&| \end{array}} \right]} $ |
則Eigen-based濾波器輸出可由фc的正交補空間與采樣復多普勒信號向量的矩陣乘積給出：
| $\boldsymbol{y} = \left( {\boldsymbol{I} - \sum\limits_{k = 1}^{{K_c}} {{\boldsymbol{e}_k}\boldsymbol{e}_k^H} } \right)\boldsymbol{x}$ |
由式(5)、(19)可知，Eigen-based濾波器的響應由采樣復多普勒信號向量分解得到的雜波特征向量決定。當組織運動狀態發生變化引起雜波成分發生改變時，雜波特征向量的維數Kc會相應改變，導致Eigen-based壁濾波器雜波截至特性產生相應變化，因而具有自適應性。
對采樣復多普勒信號向量的協方差矩陣進行特征值分解，可獲得式(16)中的特征向量與特征值：
| ${\boldsymbol{R}_x} = E\left\{ {x{x^H}} \right\} = \sum\limits_{k = 1}^N {{\lambda _k}{\boldsymbol{e}_k}\boldsymbol{e}_k^H} $ |
在實際應用中，Rx由沿掃描深度方向相鄰L個采樣復多普勒信號向量估計得出：
| ${\boldsymbol{R}_x} \approx \frac{1}{L}\sum\limits_{l = 1}^L {{\boldsymbol{x}_l}\boldsymbol{x}_l^H} $ |
其中，L需滿足平穩性假設。
3 仿真分析
3.1 雜波抑制性能比較
為了對三種壁濾波器的雜波抑制性能進行量化評估，采用Yu[6]提出的回波合成模型產生多普勒回波信號。仿真共產生兩組多普勒回波信號，每組包含的雜波分量的帶寬不同，以對不同雜波環境下的各種壁濾波器的性能進行比較。壁濾波處理后的采樣數據采用標準一維自相關估計法[7]估計血流平均多普勒頻率。為分析壁濾波技術對不同血流流速估計的影響，設計32組歸一化平均多普勒頻率在0~0.5內均勻分布的血流信號。仿真參數如表 1所示。
表1
仿真參數設置
Table1.
Parameters used for the simulation
圖 1給出了雜波1、采樣多普勒向量信號長度N為8時，經1 000次獨立仿真實驗得到的血流平均多普勒頻率估計偏差與標準差。其中Prj-IIR為投影初始化Chebyshev濾波器，通帶歸一化截止頻率為0.1，階數為3；Pol-Reg的正交基采用Legendre多項式基函數，階數為3，以確保阻帶截止頻率與Prj-IIR相近。對于Eigen-based壁濾波器，按照Yu[8-9]提出的基于頻率比較方法自適應選取雜波空間維數。從圖 1可以看出，三種壁濾波器均有效抑制了位于阻帶內的多普勒頻率成分，導致在該頻率范圍內的頻率估計值與理論值相比產生顯著的差異，而在通帶范圍內，雜波頻率成分得以充分抑制，血流頻率成分得以保留，隨著血流多普勒頻率的增加，頻率估計產生的偏差與標準差則迅速降低并趨近于零。對比經三種壁濾波器處理后平均頻率估計結果可知，Prj-IIR與Pol-Reg性能相近，而Eigen-based壁濾波器在抑制雜波的同時可盡量保證血流成分不被衰減，因而可在較低血流多普勒頻率(250 Hz)時達到穩態估計結果。
圖1
血流平均頻率估計(雜波1，N=8)
Figure1.
Blood mean frequency estimation with clutter signal 1 (package size N=8)
為比較多普勒向量信號采樣長度對壁濾波器性能的影響，圖 2給出了將N從8增加到16的估計結果。為保證近似阻帶帶寬，Prj-IIR濾波器的階數增加為4，通帶歸一化截止頻率為0.05；Pol-Reg階數為4階。比較圖 1與圖 2可以發現，三種不同類型壁濾波器均可對1型強度、中心頻率與帶寬的雜波信號進行有效抑制，通過增加多普勒向量信號采樣長度可以顯著降低過渡帶寬，增強阻帶抑制能力。
圖2
血流平均頻率估計(雜波1，N=16)
Figure2.
Blood mean frequency estimation with clutter signal 1 (package size N=16)
圖 3給出了雜波2、采樣多普勒向量信號長度N為8時的估計結果。濾波器參數設計同圖 1。從圖 3中可以看出2型雜波嚴重影響了平均頻率估計結果。由于雜波帶寬為1型雜波的2.5倍，在不改變設計參數的情況下，Prj-IIR與Pol-Reg濾波器無法對雜波進行有效抑制，導致部分低頻雜波成分混入血流多普勒信號，造成血流平均頻率估計低于理論值，因而在高頻時出現顯著的負估計偏差，估計標準差也明顯高于圖 1中所示的針對雜波1的估計結果。而Eigen-based壁濾波器根據輸入信號的二階統計特性自適應選擇雜波空間維數，阻帶帶寬與圖 1相比顯著增加以確保雜波成分的有效抑制，因而可實現通帶內平均多普勒頻率的準確估計。
圖3
血流平均頻率估計(雜波2，N=8)
Figure3.
Blood mean frequency estimation with clutter signal 2 (package size N=8)
圖 4為將N從8增加到16時，針對雜波2的平均頻率估計結果。Prj-IIR與Pol-Reg設計參數同圖 2。對比圖 3與圖 4可知，增加多普勒向量信號采樣長度可以改善Prj-IIR與Pol-Reg的雜波抑制性能，減小平均頻率估計誤差，但與圖 1和圖 2的處理結果相比，通帶標準差仍存在較大差異，因而無法有效抑制雜波對平均頻率估計的影響。對于Eigen-based壁濾波器，N的增加使得過渡帶變窄，產生較大平均頻率，估計偏差的頻率范圍進一步減小。
圖4
血流平均頻率估計(雜波2，N=16)
Figure4.
Blood mean frequency estimation with clutter signal 2 (package size N=16)
3.2 計算復雜度比較
Prj-IIR、Pol-Reg為非自適應濾波器，濾波操作為矩陣乘法運算，因此，系數矩陣A的維數決定了算法的復雜度。Eigen-based壁濾波器的計算包括數據樣本的協方差矩陣估計、特征值分解與濾波運算三部分，完成濾波操作所需的浮點運算次數(float-point operations per second,FLOPS)由多普勒向量維數N、軸向采樣體積數、用以估計數據樣本的協方差矩陣的多普勒向量數L與雜波空間維數Kc共同決定。圖 5示出了三種濾波算法計算量比較，軸向采樣體積數為150。從圖 5可知，Prj-IIR與Pol-Reg計算量相同，當N≤10時，Eigen-based的計算量與Prj-IIR、Pol-Reg相比有所增加，但隨著N的增加，計算所需FLOPS則顯著增加。Kc對Eigen-based型壁濾波器的計算量影響不明顯，當Kc從2增加到5時，計算量略有增加，當N≤10時，增加量基本可以忽略。
圖5
浮點計算量比較
Figure5.
Flop counts comparison
4 結論
雜波抑制是彩色超聲血流信號處理的一個重要環節，壁濾波的雜波抑制性能直接影響血流參數估計精度。Prj-IIR、Pol-Reg與Eigen-based壁濾波器均通過計算輸入多普勒向量在雜波信號基向量上的投影，并從輸入信號中減去該投影值，實現雜波抑制。在阻帶帶寬相同時，Prj-IIR與Pol-Reg雜波抑制性能相近。對于Prj-IIR濾波器，增加階數的同時保持低的通帶截止頻率可增加阻帶帶寬，而Pol-Reg壁濾波器的阻帶帶寬僅取決于設計階數，Eigen-based壁濾波器根據雜波統計特性動態改變雜波截止特性。針對靜態或均勻運動組織，三種壁濾波技術均可有效抑制雜波，而增加采樣多普勒向量樣本有助于進一步增強阻帶抑制能力與降低過渡帶寬。對于組織空時非平穩運動，Eigen-based壁濾波技術通過動態改變雜波截止特性可實現雜波成分的充分抑制，而對于非自適應濾波器，由于阻帶截止頻率由組織運動特性的先驗知識確定，當組織運動特性發生變化導致雜波帶寬改變時，雜波成分的泄漏或血流成分的過渡抑制都將導致血流參數估計相對實際值產生較大偏差。Prj-IIR與Pol-Reg的計算復雜度為O(N2)，Eigen-based壁濾波器的計算復雜度為O(N3)，但當多普勒向量維數N≤10時，兩者計算所需FLOPS相差不顯著。雖然Eigen-based可有效提高血流信號的檢測能力，但在臨床診斷應用中，采樣多普勒向量協方差矩陣的構建、雜波子空間維數的準確選取等問題仍亟待解決。
