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

清華大學(xué) 醫(yī)學(xué)院 生物醫(yī)學(xué)工程系，北京 100084

引言

中風(fēng)是全世界死亡率第二的疾病，而在中國(guó)則是頭號(hào)殺手[1]。其中，頸動(dòng)脈粥樣硬化斑塊的破裂是導(dǎo)致缺血性中風(fēng)發(fā)生的主要原因之一[2]，因此，對(duì)斑塊的易損性進(jìn)行早期評(píng)估有助于預(yù)防缺血事件的發(fā)生。在組織學(xué)上，易損斑塊通常具有大的脂質(zhì)核、薄的纖維帽、斑塊內(nèi)出血、炎癥以及新生血管等特征，而穩(wěn)定斑塊則通常有鈣化、厚的纖維帽以及沒有脂質(zhì)核[3]。因此，根據(jù)斑塊的不同成分特征，我們可以對(duì)斑塊的易損性進(jìn)行有效評(píng)估[4]。

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

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

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

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

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

為了研究PRF和NSA這兩個(gè)因素對(duì)SAC中運(yùn)動(dòng)偽影的影響，我們?cè)诜抡嬷蟹謩e模擬了理想情況和實(shí)際（有偽影）情況。如圖1所示，在理想情況下，施加2%形變前，物體在三個(gè)角度的發(fā)射過程中位移始終為0，在施加2%形變后，物體的位移也保持不變（2%×h）。在實(shí)際情況下，施加2%形變前，物體在每?jī)蓚€(gè)相鄰角度的發(fā)射之間發(fā)生了一段位移（σ/PRF×h），同理可知，施加2%形變后，物體在每?jī)蓚€(gè)相鄰角度之間也會(huì)發(fā)生位移（σ/PRF×h），疊加2%的形變后就是物體形變后在每個(gè)角度下的位移值。

1 方法

1.1 仿真設(shè)置

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

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

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

1.2 超聲仿真數(shù)據(jù)收集和處理

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

1.3 應(yīng)變計(jì)算

我們利用基于光流法的SAC（窗長(zhǎng)：2.0mm×2.0 mm）來分別對(duì)形變前后的RF數(shù)據(jù)進(jìn)行處理，以得到軸向和側(cè)向位移估計(jì)分布。隨后，我們?cè)儆么伴L(zhǎng)為0.5 mm×0.5 mm的二維中值濾波器來去除位移場(chǎng)中的離群值。最后，我們用Savitzky-Golay（SG）差分器對(duì)濾波后的軸向和側(cè)向位移進(jìn)行差分，來得到相應(yīng)的軸向和側(cè)向應(yīng)變圖像[27]。這里，SG差分器在軸向和側(cè)向的長(zhǎng)度分別為1.0 mm和2.0 mm。

1.4 成像性能評(píng)估

我們用SNR和CNR作為應(yīng)變估計(jì)結(jié)果的評(píng)估標(biāo)準(zhǔn)。這里，SNR的定義如下所示：

其中，es是感興趣區(qū)域（Region of Interest，ROI）中軸向或側(cè)向應(yīng)變的平均值，而σs是應(yīng)變估計(jì)的標(biāo)準(zhǔn)差(Standard Deviation，SD)。

CNR的定義為：

2 結(jié)果

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

圖2 數(shù)值仿體的軸向應(yīng)變圖像

圖3 數(shù)值仿體的側(cè)向應(yīng)變圖像

不同PRF和NSA下軸向應(yīng)變和側(cè)向應(yīng)變的SNR和CNR如圖4所示。

圖4 不同PRF和NSA下軸向應(yīng)變和側(cè)向應(yīng)變的SNR和CNR

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

3 討論

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

在圖4中，NSA=1表現(xiàn)為一條水平直線。這是因?yàn)镹SA = 1代表只有一個(gè)波束發(fā)射和接收，沒有SAC，即不存在運(yùn)動(dòng)偽影，因此SNR和CNR不會(huì)隨PRF變化。此外，當(dāng)NSA不變時(shí)，使用SAC（即NSA=3，5或7）的SNR和CNR隨著PRF增大而增大，在當(dāng)PRF達(dá)到4 kHz后則趨于定值。這是因?yàn)槲矬w在多角度發(fā)射和接收過程中會(huì)產(chǎn)生形變，從而產(chǎn)生偽影。隨著PRF增大，每個(gè)發(fā)射事件之間的間隔時(shí)間會(huì)縮短，物體在每個(gè)發(fā)射事件下的形變量也會(huì)變小，運(yùn)動(dòng)偽影會(huì)隨著減小，從而應(yīng)變圖像質(zhì)量會(huì)變好，因此SNR和CNR會(huì)提高。

當(dāng)PRF小于1 kHz時(shí)，使用SAC的SNR和CNR低于沒有使用SAC（即NSA=1）的結(jié)果；當(dāng)PRF達(dá)到1 kHz時(shí)，趨勢(shì)開始發(fā)生反轉(zhuǎn)。一個(gè)合理的推測(cè)是，當(dāng)PRF較小時(shí)，SAC的每個(gè)發(fā)射事件的時(shí)間比較長(zhǎng)，物體在每個(gè)發(fā)射事件下的形變量較大，因此使用SAC時(shí)的運(yùn)動(dòng)偽影較大，而隨著PRF逐漸提高，運(yùn)動(dòng)偽影的問題才得到補(bǔ)償。

此外，當(dāng)PRF小于1 kHz時(shí)，角度數(shù)量越小時(shí)（NSA越?。琒NR和CNR越大，隨著PRF不斷增加，最終不同角度數(shù)量的結(jié)果接近一致。這是因?yàn)楫?dāng)PRF較低時(shí)，NSA越大，就會(huì)產(chǎn)生越多角度之間的形變，因此使用SAC時(shí)的運(yùn)動(dòng)偽影越大；當(dāng)PRF較高時(shí)，各角度之間的形變逐漸縮小，最終導(dǎo)致NSA對(duì)SNR和CNR的影響逐漸消失。

當(dāng)PRF＞4 kHz時(shí)，NSA為3，5和7的應(yīng)變結(jié)果相似，但在相同PRF下，NSA為3的有效幀頻高于NSA為5或7的有效幀頻。根據(jù)此初步發(fā)現(xiàn)，本研究建議將NSA設(shè)置為3，可以獲得較好的應(yīng)變估計(jì)。考慮到未來應(yīng)用到商用超聲成像系統(tǒng)上時(shí)，其存儲(chǔ)空間有限，在保證成像質(zhì)量的前提下，采用較低的PRF可以有效節(jié)省存儲(chǔ)空間，也可以節(jié)省數(shù)據(jù)傳輸和處理的時(shí)間，本研究建議可獲得較好應(yīng)變估計(jì)的最低PRF是4 kHz。

在實(shí)際情況中，頸動(dòng)脈斑塊的形變和運(yùn)動(dòng)可能更為復(fù)雜，因此本研究所提出的方法有其局限性。在仿真實(shí)驗(yàn)中，我們?cè)O(shè)定的壓縮方向和大小都是恒定的，然而頸動(dòng)脈斑塊的運(yùn)動(dòng)方向和形變可能會(huì)隨時(shí)間變化。因此，在實(shí)際的臨床應(yīng)用中，超聲數(shù)據(jù)采集時(shí)間不能過長(zhǎng)，盡可能避免斑塊在采集期間有太多復(fù)雜運(yùn)動(dòng)。未來我們也會(huì)用在體實(shí)驗(yàn)對(duì)仿真的結(jié)論做進(jìn)一步的驗(yàn)證。

4 結(jié)論

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