【作 者】辛莉莉，張勇德

中國醫(yī)科大學(xué)公共基礎(chǔ)學(xué)院，沈陽市，110122

0 引言

通常，正常組織和病灶組織之間存在彈性差異，而這種差異在彈性成像模式中清晰可見，超聲彈性成像這種能反應(yīng)組織力學(xué)屬性的性能彌補(bǔ)了醫(yī)學(xué)成像模態(tài)的不足[1-2]，在近些年得到了迅速發(fā)展。

超聲彈性成像按照組織激勵(lì)方式的不同可分為準(zhǔn)靜態(tài)彈性成像和聲輻射力彈性成像等[3]多種方式，其中準(zhǔn)靜態(tài)彈性成像是最基本的方法，很多其他方式的超聲彈性成像也是用同樣或類似的方法進(jìn)行位移估計(jì)或者應(yīng)變估計(jì)[4]。在準(zhǔn)靜態(tài)超聲彈性成像中，當(dāng)用探頭對(duì)組織施加壓力時(shí)，組織作為一種線彈性體[5]，不僅有平移運(yùn)動(dòng)，而且還有形變?；谖灰铺荻鹊臏?zhǔn)靜態(tài)超聲彈性成像算法可分為兩步：假設(shè)局部組織運(yùn)動(dòng)只存在平移，首先計(jì)算壓縮前、后射頻回波信號(hào)在時(shí)域上的相對(duì)位移；根據(jù)相對(duì)位移利用差分求出組織應(yīng)變[6]。該類算法一般都假設(shè)壓縮后的信號(hào)是壓縮前信號(hào)的時(shí)移，忽略了形變。當(dāng)信號(hào)形變較大時(shí)，會(huì)出現(xiàn)嚴(yán)重的去相關(guān)，進(jìn)而導(dǎo)致位移估計(jì)過程中會(huì)引入去相關(guān)噪聲，而位移到應(yīng)變的差分過程會(huì)進(jìn)一步放大噪聲[7]。因此，該類算法的應(yīng)變圖像質(zhì)量較差，信噪比低。

仿射變換是一種六參數(shù)組織運(yùn)動(dòng)模型，它綜合了平移、旋轉(zhuǎn)、縮放、切變等線性變換[8-9]，既包含了位移，又包含了形變，更符合局部組織運(yùn)動(dòng)，而且基于該模型的算法可以同時(shí)得到位移和應(yīng)變信息，避免了位移到應(yīng)變的差分過程對(duì)噪聲的放大。

光流法是運(yùn)動(dòng)估計(jì)的重要算法之一，光流計(jì)算是利用相鄰兩幀圖像中像素強(qiáng)度信息的時(shí)域變化和相關(guān)性來確定各個(gè)像素的運(yùn)動(dòng)[10-11]。光流法本身有良好的并行性，并且可得到亞像素單位的位移[12]。該研究綜合了仿射變換和光流法，在先驗(yàn)估計(jì)的思想下[13]，在小范圍內(nèi)進(jìn)行運(yùn)動(dòng)追蹤，提高了算法成像質(zhì)量和運(yùn)行速度。

1 方法

1.1 光流法

光流法實(shí)際上是一種根據(jù)要檢測圖像像素灰度隨時(shí)間的變化來了解物體的移動(dòng)速度和方向的手段[10]。根據(jù)圖像灰度一致性假設(shè)，可推導(dǎo)出光流的基本約束方程。在光流基本方程的基礎(chǔ)上，HORN和SCHUNCK[11]提出了全局平滑性約束；LUCAS和KANADE[10]引入了局部平滑約束，即假設(shè)在一個(gè)很小的空間鄰域內(nèi)光流保持恒定。在一個(gè)小的空間鄰域上，假設(shè)鄰域內(nèi)數(shù)據(jù)點(diǎn)個(gè)數(shù)為n，對(duì)不同的點(diǎn)賦予不同的權(quán)值，光流估計(jì)誤差為：

其中，Ix，Iy，It是圖像灰度在x，y，t方向上的偏導(dǎo)數(shù)，u，v是像素在x,y方向上的位移，由此可知，式(1)只有兩個(gè)未知運(yùn)動(dòng)參數(shù)，即光流法是假設(shè)組織在施壓時(shí)的局部運(yùn)動(dòng)是平移運(yùn)動(dòng)。

在準(zhǔn)靜態(tài)超聲彈性成像中，當(dāng)用探頭對(duì)組織施加壓力時(shí)，以探頭為參考點(diǎn)，與探頭直接接觸的區(qū)域相對(duì)位移較小，滿足光流法的基本約束條件，可直接用光流法進(jìn)行計(jì)算。隨著深度的增加，組織與探頭的相對(duì)位移也增加，尤其是遠(yuǎn)離探頭的邊緣區(qū)域，形變較大，此時(shí)如果直接用光流法進(jìn)行估計(jì)，成像質(zhì)量差。

1.2 改進(jìn)措施

為了補(bǔ)償壓縮后信號(hào)的形變影響，CESPEDES等[14]提出對(duì)壓縮后的信號(hào)進(jìn)行拉伸，從而增加壓縮前后信號(hào)的相關(guān)性。但是該補(bǔ)償措施是在組織局部的運(yùn)動(dòng)為平移的假設(shè)條件下進(jìn)行的，忽略了組織的形變，即使對(duì)信號(hào)進(jìn)行補(bǔ)償，最后估計(jì)的運(yùn)動(dòng)也存在誤差。

1.2.1 仿射變換模型

仿射變換是一種組織運(yùn)動(dòng)模型，可用來估計(jì)組織運(yùn)動(dòng)[15-17]。假設(shè)點(diǎn)A是某一窗口的中心點(diǎn)，B是窗口內(nèi)的除中心點(diǎn)的任意一點(diǎn)，A’是窗口形變后的中心點(diǎn)，B形變后是B’。B到B’的形變可以用仿射變換來表示[17]：

其中，xA和yA，xB和yB分別表示點(diǎn)A和點(diǎn)B在x軸和y軸的坐標(biāo)，Δx=xB-xA，Δy=yB-yA表示A和B在x方向和y方向的距離，εxx， εxy， εyx， εyy分別表示橫向應(yīng)變，橫向剪切應(yīng)變，軸向剪切應(yīng)變和軸向應(yīng)變。該算法可同時(shí)得到位移和應(yīng)變信息。

把仿射變換等式代入到式 (1) 中即可得到新算法的基本方程：

此時(shí)，加權(quán)系數(shù)全部為1，利用最小二乘法對(duì)6個(gè)未知參數(shù)進(jìn)行求解。

1.2.2 先驗(yàn)估計(jì)

由于基本約束方程中的假設(shè)，光流法在大位移情況時(shí)會(huì)估計(jì)失敗[12]。為了解決該問題，有學(xué)者在Lucas-Kanade光流法的基礎(chǔ)上利用金字塔分層技術(shù)[18-19]來求解位移。首先把數(shù)據(jù)進(jìn)行金字塔分層，從最頂層開始求解光流，將其結(jié)果作為下一層的初始值，重復(fù)上述過程直至最底層，這種由粗到細(xì)的方法很好地將大位移分解成小位移。

有些學(xué)者在仿射變換的基礎(chǔ)上進(jìn)行組織的運(yùn)動(dòng)估計(jì)，如LIU等[15]用Powell方法迭代求出位移和應(yīng)變；MERCURE等[16]基于拉格朗日斑點(diǎn)模型，利用Levenberg-Marquardt方法迭代求解，實(shí)現(xiàn)無創(chuàng)血管超聲彈性成像；PAN等[17]提出了兩步光流法可同時(shí)得到位移和應(yīng)變信息。

然而，由于原始射頻回波信號(hào)的數(shù)據(jù)量大，上述方法的實(shí)時(shí)性都不高。該研究采用了先驗(yàn)估計(jì)的思想，利用相鄰點(diǎn)已知位移作為當(dāng)前點(diǎn)位移估計(jì)的先驗(yàn)值[13]，縮小搜索范圍，既滿足了光流法的小位移條件下的假設(shè)，又大大提高成像速度。但由于算法中存在窗移，因此需要對(duì)位移和應(yīng)變數(shù)據(jù)進(jìn)行二維線性插值。算法流程如圖1所示。

圖1 算法流程圖Fig.1 The flow diagram of algorithm

為了方便定量分析，新提出的算法與一種局部圖像形變算法[20]的成像質(zhì)量進(jìn)行對(duì)比。對(duì)比算法首先用塊匹配法獲得初始整數(shù)位移，然后假設(shè)組織局部形變?yōu)榉律渥儞Q，利用光流法對(duì)初始位移進(jìn)行修正。

1.3 評(píng)價(jià)指標(biāo)

本文主要從應(yīng)變成像質(zhì)量和算法速度兩方面對(duì)算法性能進(jìn)行評(píng)價(jià)。通常利用信噪比（SNRe）和對(duì)比度噪聲比（CNRe）來評(píng)價(jià)應(yīng)變圖像質(zhì)量和病灶的檢測能力[15,17,21]，用時(shí)長來評(píng)價(jià)算法速度。

SNRe定義為：

其中，es和σs是在均勻彈性模量區(qū)域的平均值和標(biāo)準(zhǔn)差，在該文中均勻彈性模量區(qū)域有嵌入物目標(biāo)區(qū)和背景區(qū)。除了SNRe外，CNRe也常用于在彈性成像中， CNRe的值越高表示檢測病變的能力越好。

CNRe定義為：

其中，eB和σB別是背景區(qū)域的均值和標(biāo)準(zhǔn)差， eT和σT是目標(biāo)區(qū)域的均值和標(biāo)準(zhǔn)差。

1.4 仿真實(shí)驗(yàn)

利用Comsol Multiphysics5.2軟件構(gòu)建一個(gè)大小為100 mm×10 mm×50 mm的三維模型，中心處有一個(gè)半徑為5 mm的球，球形區(qū)域的彈性模量為基質(zhì)區(qū)域的10倍，泊松比為0. 495。由于壓縮量過大沒有實(shí)際應(yīng)用意義，所以本文仿真了2%，5%，8%三種壓縮量的數(shù)據(jù)。設(shè)置邊界條件，上表面分別指定位移為-1 mm、-2.5 mm和-4 mm，即代表壓縮量為2%、5%和8%，下表面固定，其他表面為自由邊界條件。FieldⅡ模擬超聲場，獲取壓縮前后的射頻回波數(shù)據(jù)。仿真數(shù)據(jù)的中心頻率為7.5 MHz，采樣頻率為100 MHz，散射點(diǎn)個(gè)數(shù)為50 000，最終的仿真數(shù)據(jù)大小為6 501×449。信號(hào)保存之后在MatLab離線處理。本研究只對(duì)比兩種算法的軸向運(yùn)動(dòng)，即沿著超聲探頭發(fā)射方向。

2 結(jié)果

在MatLab平臺(tái)上用仿真射頻回波數(shù)據(jù)對(duì)兩種算法進(jìn)行實(shí)驗(yàn)，其中本研究算法中的窗口大小為3.08 mm×2.87 mm，窗口之間的重疊率為63%。從圖2和圖3中可以看出，在2%和5%壓縮量下兩種算法估計(jì)的位移和應(yīng)變差別不太明顯。在8%壓縮量情況下，如圖4(c)所示，對(duì)比算法得到的位移圖中間部分出現(xiàn)大的噪聲斑塊；應(yīng)變圖中間硬物部分整體邊界不規(guī)則，背景區(qū)域噪聲過多，整體位移和應(yīng)變成像質(zhì)量差。從圖2，圖3和圖4無法直接得到應(yīng)變圖成像質(zhì)量的好壞，因此，需要對(duì)兩種算法估計(jì)得到的應(yīng)變圖成像質(zhì)量進(jìn)行定量分析，結(jié)果如表1所示。

圖2 2%位移和應(yīng)變計(jì)算結(jié)果Fig.2 Displacement and strain calculated results at 2%

圖3 5%位移和應(yīng)變計(jì)算結(jié)果Fig.3 Displacement and strain calculated results at 5%

圖4 8%位移和應(yīng)變計(jì)算結(jié)果Fig.4 Displacement and strain calculated results at 8%

表1 信噪比和對(duì)比度噪聲比結(jié)果Tab.1 The result of SNRe and CNRe

計(jì)算信噪比和對(duì)比度噪聲比時(shí)，在應(yīng)變圖的目標(biāo)區(qū)域和背景區(qū)域選取兩個(gè)等大矩形，如圖2(e)所示，大小為4.61 mm×8.00 mm，定量指標(biāo)結(jié)果如表1所示。在2%和5%壓縮量時(shí)，新提出的算法估計(jì)得到的應(yīng)變圖在目標(biāo)區(qū)域的信噪比和對(duì)比度信噪比都高于引用算法，但兩種算法的背景區(qū)域的信噪比相差不大，說明在小應(yīng)變情況下，雖然兩種算法對(duì)病灶的檢測能力不同，但應(yīng)變圖的背景區(qū)域成像質(zhì)量相差不大。在8%壓縮量時(shí)，新提出的算法應(yīng)變圖定量指標(biāo)明顯優(yōu)于對(duì)比算法，其中目標(biāo)區(qū)域的信噪比高出3.98 dB，背景區(qū)域的信噪比高出3.94 dB，表明該研究算法在8%大壓縮量下應(yīng)變成像質(zhì)量較好；對(duì)比度噪聲比高8.29 dB，說明本研究算法對(duì)病灶的檢測能力較強(qiáng)。對(duì)比上述結(jié)果，說明本文算法可以在大壓縮量情況下相對(duì)準(zhǔn)確地估計(jì)出軸向位移和應(yīng)變結(jié)果。

在相同的計(jì)算環(huán)境下（Dell，內(nèi)存8 GB，處理器Inter(R) Xeon(R)CPU E5-2665 0 @2. 40 GHz 2.40 型號(hào)PrecisionT7600）對(duì)兩種算法進(jìn)行實(shí)驗(yàn)，隨著壓縮量的增加，對(duì)比算法的計(jì)算時(shí)間明顯增加，如表2所示，主要由于對(duì)比算法中第一步利用塊匹配來估計(jì)整數(shù)位移，當(dāng)壓縮量增加時(shí)，為了相對(duì)準(zhǔn)確地估計(jì)出初始位移，塊匹配搜索參數(shù)也要增大，所以對(duì)比算法在未知壓縮量情況下使用時(shí)存在劣勢。本文算法計(jì)算時(shí)間基本穩(wěn)定在8.3 s左右，其遠(yuǎn)遠(yuǎn)小于引用算法，因此，本文算法在實(shí)時(shí)顯示上更有優(yōu)勢。

表2 兩種算法在不同壓縮量下的時(shí)長Tab.2 The time of two algorithm under different compression

3 總結(jié)

本文闡述了一種基于仿射變換的光流超聲彈性成像算法，引用了先驗(yàn)估計(jì)的思想，不僅克服了光流法只能進(jìn)行小位移估計(jì)的缺點(diǎn)，而且削弱了基于位移梯度的算法中由于模型不精確造成位移估計(jì)誤差，同時(shí)估計(jì)出位移和應(yīng)變信息。仿真實(shí)驗(yàn)證明，本研究算法在大壓縮量情況時(shí)，應(yīng)變圖信噪比、對(duì)比度噪聲比和算法時(shí)長與對(duì)比算法相比都有明顯提升，對(duì)實(shí)時(shí)超聲彈性成像的應(yīng)用有良好的推進(jìn)作用。