空間角度復合應變估計中的運動偽影產生因素研究

潘宗輝，劉志，李富峰，何瓊

清華大學 醫學院 生物醫學工程系，北京 100084

引言

中風是全世界死亡率第二的疾病，而在中國則是頭號殺手[1]。其中，頸動脈粥樣硬化斑塊的破裂是導致缺血性中風發生的主要原因之一[2]，因此，對斑塊的易損性進行早期評估有助于預防缺血事件的發生。在組織學上，易損斑塊通常具有大的脂質核、薄的纖維帽、斑塊內出血、炎癥以及新生血管等特征，而穩定斑塊則通常有鈣化、厚的纖維帽以及沒有脂質核[3]。因此，根據斑塊的不同成分特征，我們可以對斑塊的易損性進行有效評估[4]。

目前已經有多種影像技術可以用于頸動脈斑塊易損性的評估，如X射線電子計算機斷層掃描（Computed Tomography，CT）、磁共振成像（Magnetic Resonance Imaging，MRI）和超聲成像（Ultrasound Imaging）等[5]。CT會對人體產生電離輻射危害，MRI掃描價格較為昂貴且耗時長[6-8]，超聲成像具有掃描價格相對便宜、實時成像、無電離輻射危害等優勢，近年來，超聲彈性成像技術被提出用來識別易損斑塊[9]。彈性成像方法首先被應用于血管內超聲（Intravascular Ultrasound，IVUS）來區分穩定和易損的冠狀動脈斑塊[10-12]。雖然IVUS彈性成像具有較高的靈敏度和特異度，但是它具有侵入性且費用昂貴。鑒于IVUS的局限性，有研究者開始發展非侵入性血管彈性成像技術，例如頸動脈彈性成像[13-16]。頸動脈彈性成像通過運動估計方法算出斑塊的位移和應變（率）分布[17-18]，然而低精度的側向估計降低了短軸掃描下的成像性能。為了提高側向估計的精度，近年來有研究者已經提出了側向插值算法和在側向點擴散函數上加振蕩的方法[19-20]。

空間角度復合（Spatial Angular Compounding，SAC）成像是另一種提高側向估計精度的方法[21-24]。在多角度平面波SAC中，軸向（豎直）位移和側向（水平）位移都是通過多個角度的軸向位移復合得到的。如圖1理想情況所示，我們分別在物體形變前和形變后采集多個角度的射頻（Radiofrequency，RF）數據，并通過二維互相關算法分別得到沿各個角度方向的軸向位移。隨后，最終的軸向位移uver和側向位移uhor可以通過最小二乘法從這些軸向位移中重建得到，其具體公式如下：

其中AT是矩陣A的轉置，θi（i=0,1，...，n）是轉向角，n是NSA，uax，θi則代表沿偏轉角θi的軸向位移。

上述情況考慮的是準靜態彈性成像過程，即對物體施加一個靜態壓縮，分別在壓縮前后對物體成像，進而得到位移估計。然而在實際情況中，物體在多角度發射過程中也有運動和形變，這可能會導致運動偽影的產生。圖1分別展示了基于理想情況和實際情況下的SAC成像過程，我們設定空間復合的角度數量為3（即NSA=3），PRF保持不變，物體的高度為h。在理想情況下，我們在對物體施加2%的壓縮前后各采3個角度的數據，獲得各個角度下的位移估計，進而通過角度復合得到更準確的結果。在實際情況下，在多個角度發射過程中，應變率固定為σ的軸向應變被施加到物體上表面，用來模擬物體在SAC過程中的運動和形變。當σ較小或PRF足夠高時，可以認為不同偏轉角下的數據是同時獲得的，因此相鄰角度的發射過程之間沒有運動偽影。然而在實際的頸動脈彈性成像中，當PRF不夠高或NSA過多時，斑塊的運動和形變將不能忽略，可能會產生運動偽影。

圖1 基于理想情況和實際情況下的SAC成像

為了研究PRF和NSA這兩個因素對SAC中運動偽影的影響，我們在仿真中分別模擬了理想情況和實際（有偽影）情況。如圖1所示，在理想情況下，施加2%形變前，物體在三個角度的發射過程中位移始終為0，在施加2%形變后，物體的位移也保持不變（2%×h）。在實際情況下，施加2%形變前，物體在每兩個相鄰角度的發射之間發生了一段位移（σ/PRF×h），同理可知，施加2%形變后，物體在每兩個相鄰角度之間也會發生位移（σ/PRF×h），疊加2%的形變后就是物體形變后在每個角度下的位移值。

1 方法

1.1 仿真設置

在有限元分析軟件FEMLAB 3.5a（Comsol Inc.,Burlington, MA, USA）中，我們構建了一個三維數值仿體。該仿體的尺寸為40 mm×40 mm×6 mm（深度×寬度×厚度），其中心嵌入一個直徑為4.05 mm的圓柱形異物。圓柱形異物和周圍背景的楊氏模量分別為80 kPa和25 kPa。

許多學者都采用該模型對應變成像的基本參數進行研究，這種結構簡單的模型也比較適合對基本成像方法或參數進行早期研究[18]。為了模擬頸動脈彈性成像過程中斑塊的運動和形變，我們對物體施加了一個軸向壓縮。頸動脈斑塊的軸向應變率一般不超過2 s-1，因此仿真中該軸向壓縮的應變率設為 2 s-1[14]。

我們用Field II軟件模擬了一個具有192陣元的線陣探頭，其有效陣元數為128[25-26]，并將仿體放置在探頭下方10 mm處進行成像。探頭的陣元間距設為0.20 mm，中心頻率設為6.25 MHz，采樣率設為100 MHz。在控制探頭發射過程中，我們設置了多組不同的PRF和NSA。PRF的變化范圍是200 Hz～10 kHz，基本涵蓋了超聲機器常用的PRF范圍，具體包括 200 Hz，500 Hz，800 Hz，1 kHz，1.5 kHz，2 kHz，3 kHz，4 kHz，5 kHz，6 kHz，8 kHz，10 kHz。我們在設置發射角度時都包含0°以提高軸向估計的精度，且正負角度對稱，因此NSA都是奇數。常用的SAC角度有3和7，且PRF較低時設置9個或者更多的角度數量會帶來更大的偽影問題。因此，NSA 取值為 1（0°），3（0°，±15°），5（0°，±7.5°，±15°）和 7（0°，±5°，±10°，±15°）。

1.2 超聲仿真數據收集和處理

首先，我們以Field II軟件對數值仿體進行仿真計算，獲得仿體形變前的超聲通道數據。接著，我們按照圖1介紹的方法對仿體施加軸向壓縮產生形變。仿體初始的散射子密度設為12.5個/mm3，其位置和幅度均為隨機分布。形變后仿體中散射子的位置是通過FEMLAB仿真得到的理論位移加上形變前散射子的位置獲得，而其幅度保持不變。然后，我們再次用Field II軟件對形變后仿體進行仿真計算，獲得相應的超聲通道數據。最后，我們采用一個基于漢寧窗和f-number為1.5的延時疊加（Delay and Sum，DAS）算法分別對形變前后的超聲通道數據進行波束合成計算，獲得波束合成后的RF數據。

1.3 應變計算

我們利用基于光流法的SAC（窗長：2.0mm×2.0 mm）來分別對形變前后的RF數據進行處理，以得到軸向和側向位移估計分布。隨后，我們再用窗長為0.5 mm×0.5 mm的二維中值濾波器來去除位移場中的離群值。最后，我們用Savitzky-Golay（SG）差分器對濾波后的軸向和側向位移進行差分，來得到相應的軸向和側向應變圖像[27]。這里，SG差分器在軸向和側向的長度分別為1.0 mm和2.0 mm。

1.4 成像性能評估

我們用SNR和CNR作為應變估計結果的評估標準。這里，SNR的定義如下所示：

其中，es是感興趣區域（Region of Interest，ROI）中軸向或側向應變的平均值，而σs是應變估計的標準差(Standard Deviation，SD)。

CNR的定義為：

2 結果

數值仿體在2%壓縮量的形變下，設定不同PRF和NSA所得到的軸向應變圖像如圖2所示，而側向應變圖像如圖3所示。軸向應變方面，我們可觀察到，沒有使用SAC的軸向應變圖像（圖2b）的周圍背景，不如在PRF=5 kHz和NSA = 3時使用SAC（圖2c）得到的背景均勻。此外，同樣是NSA=3的SAC，PRF=5 kHz時的軸向應變圖像（圖2c）背景，比PRF = 0.5 kHz時的（圖2d）更接近理論結果（圖2a）。當PRF=0.5 kHz時，與NSA=3的軸向應變圖像（圖2d）相比，NSA=7的圖像（圖2e）中心不再是圓形，且周圍背景更不均勻。側向應變方面，我們可觀察到，沒有使用SAC的側向應變圖像（圖3b）的異物變小，不如在PRF=5 kHz和NSA=3時使用SAC得到的側向應變圖像（圖3c）。圖3c的結果更接近理論結果（圖3a）。此外，同樣是NSA=3的SAC，PRF=0.5 kHz的側向應變圖像（圖3d）出現了多個低應變的藍色區域，與理論結果（圖3a）的差異較大，且周圍背景不如PRF=5 kHz（圖3c）時均勻。當PRF=0.5 kHz時，與NSA = 3的側向應變圖像（圖3d）相比，NSA=7時圖像（圖3e）的低應變藍色和高應變紅色區域出現了更大程度的位置偏移。

圖2 數值仿體的軸向應變圖像

圖3 數值仿體的側向應變圖像

不同PRF和NSA下軸向應變和側向應變的SNR和CNR如圖4所示。

圖4 不同PRF和NSA下軸向應變和側向應變的SNR和CNR

其中，NSA=1指的是只有一個波束發射和接收，沒有使用SAC的情況。軸向和側向應變的SNR和CNR隨著PRF的增大而增大，當PRF達到4 kHz時保持穩定。當PRF＞1 kHz時，使用SAC（NSA=3，5或7）得到的應變的SNR和CNR比沒有使用SAC（NSA = 1）的高。當PRF＜1 kHz時，NSA越大，SNR和CNR越低，即NSA=1，3，5，7的SNR和CNR依次變低。而當PRF達到4 kHz后，不同NSA的SNR和CNR則較為接近。

3 討論

在基于SAC的應變估計中，物體的運動和形變可能會導致運動偽影，而采用合適的成像參數可以減小或避免運動偽影，提高應變估計性能。之前雖然有學者對基于SAC位移和應變估計的最優參數進行了研究[22-22,24,28]，但是據我們所知，目前還缺乏關于其偽影問題的系統研究。本研究通過仿真實驗研究了PRF和NSA對SAC運動偽影的影響，并有望在進一步的在體實驗驗證后提出具有臨床指導意義的成像參數，以此來減小SAC的運動偽影，提升頸動脈彈性成像評估易損斑塊的性能。

在圖4中，NSA=1表現為一條水平直線。這是因為NSA = 1代表只有一個波束發射和接收，沒有SAC，即不存在運動偽影，因此SNR和CNR不會隨PRF變化。此外，當NSA不變時，使用SAC（即NSA=3，5或7）的SNR和CNR隨著PRF增大而增大，在當PRF達到4 kHz后則趨于定值。這是因為物體在多角度發射和接收過程中會產生形變，從而產生偽影。隨著PRF增大，每個發射事件之間的間隔時間會縮短，物體在每個發射事件下的形變量也會變小，運動偽影會隨著減小，從而應變圖像質量會變好，因此SNR和CNR會提高。

當PRF小于1 kHz時，使用SAC的SNR和CNR低于沒有使用SAC（即NSA=1）的結果；當PRF達到1 kHz時，趨勢開始發生反轉。一個合理的推測是，當PRF較小時，SAC的每個發射事件的時間比較長，物體在每個發射事件下的形變量較大，因此使用SAC時的運動偽影較大，而隨著PRF逐漸提高，運動偽影的問題才得到補償。

此外，當PRF小于1 kHz時，角度數量越小時（NSA越小），SNR和CNR越大，隨著PRF不斷增加，最終不同角度數量的結果接近一致。這是因為當PRF較低時，NSA越大，就會產生越多角度之間的形變，因此使用SAC時的運動偽影越大；當PRF較高時，各角度之間的形變逐漸縮小，最終導致NSA對SNR和CNR的影響逐漸消失。

當PRF＞4 kHz時，NSA為3，5和7的應變結果相似，但在相同PRF下，NSA為3的有效幀頻高于NSA為5或7的有效幀頻。根據此初步發現，本研究建議將NSA設置為3，可以獲得較好的應變估計。考慮到未來應用到商用超聲成像系統上時，其存儲空間有限，在保證成像質量的前提下，采用較低的PRF可以有效節省存儲空間，也可以節省數據傳輸和處理的時間，本研究建議可獲得較好應變估計的最低PRF是4 kHz。

在實際情況中，頸動脈斑塊的形變和運動可能更為復雜，因此本研究所提出的方法有其局限性。在仿真實驗中，我們設定的壓縮方向和大小都是恒定的，然而頸動脈斑塊的運動方向和形變可能會隨時間變化。因此，在實際的臨床應用中，超聲數據采集時間不能過長，盡可能避免斑塊在采集期間有太多復雜運動。未來我們也會用在體實驗對仿真的結論做進一步的驗證。

4 結論

本研究通過仿真實驗，分析了PRF和NSA這兩項參數對SAC應變估計中產生運動偽影的影響。對于頸動脈彈性成像（斑塊的應變率一般小于2 s-1）[14-15]，基于PRF＞4 kHz和NSA=3的SAC可以獲得運動偽影較小的應變圖像。