【作 者】辛莉莉，張勇德

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

0 引言

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

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

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

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

1 方法

1.1 光流法

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

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

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

1.2 改進(jìn)措施

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

1.2.1 仿射變換模型

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

其中，xA和yA，xB和yB分別表示點A和點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)變。該算法可同時得到位移和應(yīng)變信息。

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

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

1.2.2 先驗估計

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

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

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

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

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

1.3 評價指標(biāo)

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

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 仿真實驗

利用Comsol Multiphysics5.2軟件構(gòu)建一個大小為100 mm×10 mm×50 mm的三維模型，中心處有一個半徑為5 mm的球，球形區(qū)域的彈性模量為基質(zhì)區(qū)域的10倍，泊松比為0. 495。由于壓縮量過大沒有實際應(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，散射點個數(shù)為50 000，最終的仿真數(shù)據(jù)大小為6 501×449。信號保存之后在MatLab離線處理。本研究只對比兩種算法的軸向運動，即沿著超聲探頭發(fā)射方向。

2 結(jié)果

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

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

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

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

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

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

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

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

3 總結(jié)

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