李元瑞，張勇德

中國(guó)醫(yī)科大學(xué) 生物醫(yī)學(xué)工程系，遼寧沈陽(yáng) 110000

基于零相位與光流法的混合超聲彈性成像算法的研究

李元瑞，張勇德

中國(guó)醫(yī)科大學(xué) 生物醫(yī)學(xué)工程系，遼寧沈陽(yáng) 110000

鑒于超聲彈性成像傳統(tǒng)的時(shí)域互相關(guān)算法在估計(jì)組織位移時(shí)效率低下這一問(wèn)題，本文提出了一種基于零相位與光流法的混合超聲彈性成像方法。首先利用零相位法來(lái)獲得組織壓縮前后的位移場(chǎng)，將結(jié)果作為光流法的初始位移；然后采用經(jīng)典光流法對(duì)為位移場(chǎng)進(jìn)行矯正，抑制由于組織橫向形變所帶來(lái)的去相關(guān)；最后采用最小二乘法對(duì)位移場(chǎng)進(jìn)行計(jì)算獲得了組織的軸向應(yīng)變。本文利用仿真數(shù)據(jù)對(duì)算法進(jìn)行了驗(yàn)證，并且定性分析了不同壓縮量下的算法的性能。實(shí)驗(yàn)結(jié)果表明，本論文中提出的零相位和光流法的混合超聲彈性成像算法所得縱向位移圖和縱向應(yīng)變圖與有限元結(jié)果基本一致，證明了算法的正確性。

超聲彈性成像；光流法；零相位迭代；去相關(guān)；最小二乘法

引言

病理變化一般都伴隨著組織硬度的變化，組織硬度變化可以從另一個(gè)角度反映腫瘤的活性。乳腺導(dǎo)管癌約比正常皮下脂肪組織硬23 dB，浸潤(rùn)性癌癥有28 dB；在前列腺癌方面，癌癥組織硬度比正常組織高10 dB以上。目前臨床上大部分醫(yī)學(xué)影像設(shè)備主要是對(duì)體內(nèi)組織結(jié)構(gòu)進(jìn)行成像，有時(shí)候醫(yī)生只依靠反映組織結(jié)構(gòu)信息的圖像，是無(wú)法做出最終的診斷的。檢查組織軟硬程度的傳統(tǒng)方式為手工觸診，然而手工觸診一般僅限于身體淺表處病灶的檢測(cè)，而位于身體深部的腫塊幾乎不能通過(guò)手工觸診來(lái)檢測(cè)，而且檢測(cè)結(jié)果在很大程度上依賴于醫(yī)師的主觀判斷。超聲彈性成像這一概念最早是由Ophir等[1-2]于1991提出的，通過(guò)外部或內(nèi)部激勵(lì)使組織產(chǎn)生微小形變，計(jì)算壓縮前后回波信號(hào)互相關(guān)系數(shù)，得到位移分布，通過(guò)計(jì)算位移梯度得到局部組織應(yīng)變。目前，根據(jù)所施加壓縮形變的激勵(lì)方式不同，可以將超聲彈性成像分為靜態(tài)彈性成像和動(dòng)態(tài)彈性成像[3-4]，準(zhǔn)靜態(tài)超聲彈性成像技術(shù)由于實(shí)現(xiàn)簡(jiǎn)單，不需要額外的輔助設(shè)備，而且易于實(shí)時(shí)成像，因此一直以來(lái)是彈性成像領(lǐng)域的研究熱點(diǎn)。目前，超聲彈性成像技術(shù)在國(guó)內(nèi)還仍處于理論研究階段，研究水平和國(guó)外還有很大的差距，其原因主要在于國(guó)內(nèi)的許多研發(fā)單位無(wú)法獲得進(jìn)行彈性評(píng)估所需的射頻回波數(shù)據(jù)。因此提高我國(guó)超聲彈性成像技術(shù)的研發(fā)水平，在臨床上普及具備該功能的超聲診斷設(shè)備，探索和擴(kuò)大其應(yīng)用范圍，對(duì)于提高我國(guó)醫(yī)療水平具有十分重要的意義。

1 材料與方法

1.1 一般資料

在超聲彈性成像算法研究中，可供研究的數(shù)據(jù)一般有3種來(lái)源：① 通過(guò)仿真軟件得到射頻仿真數(shù)據(jù)；② 仿體數(shù)據(jù)，通過(guò)制作得到的仿體，利用超聲射頻回波信號(hào)采集系統(tǒng)獲??；③ 在體數(shù)據(jù)，從實(shí)際的病患某些器官組織獲得數(shù)據(jù)。仿真數(shù)據(jù)由于組織結(jié)構(gòu)簡(jiǎn)單，而且包裹物和嵌入物彈性模量可以根據(jù)需求設(shè)置，適合用于評(píng)價(jià)算法性能。本次實(shí)驗(yàn)為了得到仿真數(shù)據(jù)，設(shè)計(jì)開(kāi)發(fā)了用于超聲彈性成像的射頻回波信號(hào)仿真系統(tǒng)。利用Matlab建立了組織仿體模型，產(chǎn)生散射點(diǎn)，并利用Comsol軟件對(duì)相應(yīng)的組織模型進(jìn)行有限元分析，得到相應(yīng)壓縮量下組織的應(yīng)變和位移分布，從而計(jì)算出仿體模型在被壓縮后各散射點(diǎn)的位置分布[5-9]，流程圖，見(jiàn)圖1。

圖1 仿真流程圖

有限元分析方法應(yīng)用在生物力學(xué)中，可以模擬出復(fù)雜的生物組織并對(duì)各部分結(jié)構(gòu)有更加深入的研究，本文就是利用有限元方法來(lái)模擬人體中有腫瘤或者囊腫的病變部位。通過(guò)有限元分析和求解，將經(jīng)過(guò)網(wǎng)格劃分后的仿體模型進(jìn)行有限元分析求解，得到相應(yīng)的仿真結(jié)果。圖2給出了1%壓縮量下，一個(gè)嵌入物仿體模型的仿真位移圖和軸向應(yīng)變圖，并以偽彩色顯示。圖2a、2b分別是三維位移圖和應(yīng)變圖。

1.2 實(shí)驗(yàn)方法

本文利用計(jì)算機(jī)仿真的射頻回波信號(hào)對(duì)混合算法的性能進(jìn)行了定性分析，算法實(shí)現(xiàn)的硬件平臺(tái)為Dell T5600圖形工作站，6核CPU，主頻2.0 GHz，內(nèi)存4 GB，編程語(yǔ)言為Matlab R2014a。算法主要參數(shù)設(shè)置如下：零相位法窗長(zhǎng)129，窗間間隔16，迭代次數(shù)5次，光流法平滑項(xiàng)因子α=640。

圖2 仿體模型

初始位移場(chǎng)的獲得，利用Field-II獲得的壓縮前后的射頻回波數(shù)據(jù)，將每一條掃描線壓縮前的射頻數(shù)據(jù)按每隔16個(gè)數(shù)據(jù)點(diǎn)設(shè)置為一個(gè)窗口，壓縮后的射頻數(shù)據(jù)的軸向窗口長(zhǎng)度要大于壓縮前的窗口，并且加上上一次迭代得到的位移偏移量，每隔一個(gè)窗口計(jì)算一次壓縮前射頻數(shù)據(jù)與壓縮后的射頻數(shù)據(jù)的復(fù)數(shù)互相關(guān)，在零相位時(shí)得到壓縮前后射頻數(shù)據(jù)的位移。

零相位迭代算法其中一條掃描線的某一窗口的位移求解示意圖（圖3）。迭代次數(shù)少得到的零相位對(duì)應(yīng)下的位移不精確，迭代次數(shù)太多影響計(jì)算效率，在本文介紹的零相位迭代算法中設(shè)置迭代次數(shù)為50[10-11]。從圖中可以看出，雖然相位逐漸逼近零相位但是并沒(méi)有達(dá)到零相位，但位移與相位已呈現(xiàn)出線性關(guān)系，本文基于這一點(diǎn)提出了對(duì)相位與位移進(jìn)行線性擬合，減少迭代次數(shù)，進(jìn)行5次迭代并保存每次迭代所得到的壓縮前后射頻數(shù)據(jù)的位移與相位。設(shè)相位為y，位移為x，對(duì)5次迭代的結(jié)果進(jìn)行直線擬合y=kx+b求得斜率與截距，在零相位時(shí)可以求得壓縮前后射頻數(shù)據(jù)的位移x=-b/k，較之于之前的50次迭代次數(shù)可以得到較好的位移結(jié)果并且節(jié)省了幾倍的時(shí)間。

零相位迭代算法即沒(méi)有做任何噪聲處理并且沒(méi)有考慮組織位移的連續(xù)性，因而得到的位移場(chǎng)噪聲比較明顯，光流法本身有平滑機(jī)制，所以本文采用直接利用零相位迭代算法得到的位移作為光流法的初始值[12]。

應(yīng)變是對(duì)得到的位移場(chǎng)進(jìn)行數(shù)值差分進(jìn)行的，這種梯度運(yùn)算對(duì)位移估計(jì)噪聲有放大的作用，最小二乘應(yīng)變估計(jì)實(shí)現(xiàn)的目標(biāo)是控制噪聲進(jìn)一步放大，在實(shí)現(xiàn)上采用分段線性曲線擬合所估計(jì)得到的位移場(chǎng)[13-14]。

圖3 位移-相移圖

2 結(jié)果

2.1 仿真數(shù)據(jù)結(jié)果分析

組織受壓后，壓縮后的射頻數(shù)據(jù)中每個(gè)點(diǎn)的位移都會(huì)發(fā)生一定程度的變化，通過(guò)Matlab畫(huà)出129條掃描線中的第60條掃描線的位移圖（圖4）。圖4a中縱坐標(biāo)是位移量，橫坐標(biāo)為射頻數(shù)據(jù)的點(diǎn)數(shù)。截取射頻數(shù)據(jù)6500個(gè)點(diǎn)，在1/100壓縮量的仿真數(shù)據(jù)中，可能存在的最大位移是-60，組織受壓后，各個(gè)點(diǎn)都有相應(yīng)的位移變化，由于仿真數(shù)據(jù)中中間區(qū)域是球形囊腫，彈性模量一致所以在掃描線中的中間區(qū)域表現(xiàn)為位移基本一致；而圖4b中，采用零相位與光流法的混合算法抗噪能力強(qiáng)，得到的單條掃描線的位移圖沒(méi)有毛刺，位移隨組織深度增加而加大，同時(shí)外基質(zhì)在外力作用下位移變化較快，而嵌入物部分位移變化緩慢，由此可以驗(yàn)證實(shí)驗(yàn)結(jié)果與理論分析一致，證明了算法的有效性。由圖4a可以看出，經(jīng)典的零相位法由于沒(méi)有位移連續(xù)性的約束，結(jié)果容易受到信號(hào)去相關(guān)因素的影響，缺乏抗噪能力，因此位移曲線存在著大量的噪聲；而本文提出的算法獲得的位移曲線非常平滑，如圖4b所示，具有很好的抗噪聲能力。圖4c～d為1%壓縮量下，零相位迭代算法與本文算法得到的縱向位移圖，圖中藍(lán)色代表位移較大，紅色代表位移小[15-17]。圖4c中雖然噪聲明顯但依然可以清晰的看到在圖形中間與仿體中間嵌入物的區(qū)域，零相位迭代算法依賴于算法中上一次計(jì)算的結(jié)果，因?yàn)樽铋_(kāi)始的初始位移并沒(méi)有引入先驗(yàn)值，而迭代算法每經(jīng)過(guò)一次迭代又會(huì)放大一次錯(cuò)誤，所以位移圖中出現(xiàn)一些跳變。在零相位與光流法的混合算法中，首先通過(guò)零相位迭代算法得到位移場(chǎng)，代入光流法中進(jìn)行計(jì)算得到最后的精確位移。

圖4 129條掃描線中的第60條掃描線的位移圖

2.2 最小二乘法應(yīng)變圖

本文是對(duì)離散位移場(chǎng)進(jìn)行局部最小二乘擬合來(lái)計(jì)算應(yīng)變場(chǎng)，針對(duì)1%、2%及4%等不同壓縮量下獲得的應(yīng)變場(chǎng)與經(jīng)典的零相位法進(jìn)行了分析和比較。1%壓縮量不同算法縱向應(yīng)變圖有限元仿真的應(yīng)變圖對(duì)比，見(jiàn)圖5；2%壓縮量的不同算法的縱向應(yīng)變圖對(duì)比，見(jiàn)圖6；本文所采用的零相位與光流法結(jié)合的方法在2%壓縮量的條件下，運(yùn)算速度與去除應(yīng)變圖噪聲的優(yōu)勢(shì)明顯好于零相位迭代算法；4%壓縮量不同算法的縱向應(yīng)變圖對(duì)比，見(jiàn)圖7。從圖7中可以看出4%壓縮量的應(yīng)變圖，當(dāng)應(yīng)變較大時(shí)尤其在邊緣處，檢測(cè)結(jié)果較差，因?yàn)樵诖髩嚎s量的情況下，邊緣形變大，窗口的相位匹配不準(zhǔn)確而且零相位迭代法的初值依賴于上次迭代的結(jié)果，導(dǎo)致位移估計(jì)錯(cuò)誤得不到理想的應(yīng)變圖，但與零相位迭代算法相比，本文的混合算法有更高的精確度，這也是使用光流法對(duì)零相位迭代算法進(jìn)行橫向位移修正的優(yōu)勢(shì)所在。

圖5 1%壓縮量不同算法縱向應(yīng)變圖有限元仿真的應(yīng)變圖

2.3 仿真數(shù)據(jù)結(jié)果分析

從仿真結(jié)果與算法運(yùn)行結(jié)果對(duì)比，嵌入物大小清晰可辨；周?chē)|(zhì)處顏色區(qū)域分布大致相似，基本上驗(yàn)證了算法的可行性，以及對(duì)噪聲的抑制功能?；诹阆辔慌c光流法的混合算法每次運(yùn)行都可以得到相同的實(shí)驗(yàn)結(jié)果說(shuō)明了具有可靠的實(shí)驗(yàn)重復(fù)性，在實(shí)際的臨床應(yīng)用中，實(shí)時(shí)顯示是一個(gè)重要需求，所以對(duì)算法計(jì)算效率的提高，節(jié)約時(shí)間成本顯得尤為重要。本實(shí)驗(yàn)中所采用的PC設(shè)備性能參數(shù)：① CPU，Inter（R）xeon（R）E5-26200；② 內(nèi)存，RAM 4.00 GB。從算法的計(jì)算時(shí)長(zhǎng)來(lái)看，所有程序代碼在MATLAB中運(yùn)行。1%壓縮量，零相位迭代算法需要469.29 s，而本文中利用位移-相位信息進(jìn)行擬合的算法只需要48.23 s，速度提升了近10倍。

圖6 2%壓縮量的不同算法的縱向應(yīng)變圖

圖7 4%壓縮量不同算法的縱向應(yīng)變圖

3 討論

本文是對(duì)超聲彈性成像算法的研究，改進(jìn)方案能有效提高計(jì)算效率，在大幅度減少迭代次數(shù)的前提下，通過(guò)直線擬合節(jié)省運(yùn)行時(shí)間。此外，零相位迭代算法沒(méi)有考慮橫向應(yīng)變的影響，用零相位迭代算法得到的位移場(chǎng)作為光流法的位移初始值，得到了經(jīng)過(guò)橫向位移修正后的縱向位移。

生物組織發(fā)生病變通常伴隨著組織硬度上的變化，即組織力學(xué)屬性會(huì)發(fā)生變化。而檢查組織硬度基本要靠觸診，超聲成像技術(shù)可對(duì)組織硬度屬性進(jìn)行成像為臨床疾病診斷提供更加豐富的功能信息，超聲彈性成像的實(shí)時(shí)性對(duì)于其臨床使用而言同樣具有非常重要的意義，因而探索計(jì)算速度快成像效果好的準(zhǔn)靜態(tài)超聲彈性成像有著廣泛的臨床應(yīng)用價(jià)值。零相位法由于采用牛頓迭代法僅通過(guò)少量迭代便可獲取組織位移，相比互相關(guān)算法具有較高的成像速度，然而該方法只考慮了組織的軸向運(yùn)動(dòng)，無(wú)法解決組織橫向形變所帶來(lái)的去相關(guān)問(wèn)題。光流法通過(guò)求解偏微分方程可以獲得組織的二維位移場(chǎng)，而且算法的效率很高。用光流法進(jìn)行組織運(yùn)動(dòng)估計(jì)時(shí)，隨著組織深度的增加，組織的位移逐漸加大，不準(zhǔn)確的初始位移將導(dǎo)致巨大的計(jì)算誤差甚至計(jì)算失敗。本文主要對(duì)準(zhǔn)靜態(tài)情況下的超聲彈性成像技術(shù)進(jìn)行了深入的研究，在了解了超聲彈性成像的基礎(chǔ)知識(shí)之后。使用超聲射頻數(shù)據(jù)，實(shí)現(xiàn)了改進(jìn)了計(jì)算效率的零相位迭代算法的位移估計(jì)，經(jīng)過(guò)高斯濾波的平滑濾波后，并且結(jié)合光流法的優(yōu)點(diǎn)，最后利用最小二乘法進(jìn)行了組織的應(yīng)變估計(jì)，得到與有限元仿真基本一致的結(jié)果。

本研究就現(xiàn)階段而言還需進(jìn)一步完成的工作有如下幾個(gè)方面：① 零相位迭代算法在大壓縮量下的應(yīng)變估計(jì)，由于壓縮量過(guò)大導(dǎo)致相位匹配不準(zhǔn)，檢測(cè)結(jié)果較差，需要找到解決的辦法；② 本文算法得出的應(yīng)變圖嵌入物邊界不光滑，對(duì)比度不好；③ 本文提出的零相位與光流法混合算法并沒(méi)有計(jì)算射頻數(shù)據(jù)的橫向位移，在以后的工作中要進(jìn)一步進(jìn)行完善。

[1] Ophir J,Cespedes I,Ponnekanti H,et al.Elastography: a quantitative method for imaging the elasticity of biological tissues[J].Ultrason imaging,1991,13(2):111-134.

[2] Ophir J, Alam SK,Garra B,et al.Elastography: ultrasonic estimation and imaging of the elastic properties of tissues[J].P I Mech Eng H-J Mec,1999,213(3):203-233.

[3] Krouskop TA,Wheeler TM,Kallel F,et al.Elastic moduli of breast and prostate tissues under compression[J].Ultrason Imaging,1998,20:260-274.

[4] Pan X,Gao J,Tao S,et al.A two-step optical flow method for strain estimation in elastography: Simulation and phantom study[J].Ultrasonics,2014,54(4):990-996.

[5] Hasan MK,Anas EMA,Alam SK,et al.Direct mean strain estimation for elastography using nearest-neighbor weighted least-squares approach in the frequency domain[J].Ultrasound Med Biol,2012,38(10):1759-1777.

[6] Alam SK,Feleppa EJ,Kalisz A,et al.Prostate elastography: preliminary in vivo results[J].Int Soc Opt Photonics,2005: 339-345.

[7] Dickinson RJ,Hill CR.Measurement of soft tissue motion using correlation between A-scans[J].Ultrason Med Bio,1982, 8(3):263-271.

[8] 劉鑫.基于超聲彈性成像技術(shù)的生物組織應(yīng)變估計(jì)算法仿真研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2012.

[9] 劉彥.人體組織的超聲彈性成像方法研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2011.

[10] 王勇萍.人體組織實(shí)時(shí)超聲彈性成像技術(shù)研究[D].哈爾濱:哈爾濱工業(yè)大學(xué),2012.

[11] 袁金偉.超聲彈性成像中的位移和應(yīng)變估計(jì)[D].浙江:浙江大學(xué),2012.

[12] 李偉方.超聲彈性成像全局位移場(chǎng)估計(jì)算法研究[D].武漢:武漢理工大學(xué),2013.

[13] 張耀南,李宏亮,康雁.基于二維壓縮擴(kuò)展和光流法的組織應(yīng)變估計(jì)方法[J].東北大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,539(9):22-26.

[14] Barron JL,Fleet DJ,Beauchemin SS,et al.Performance of optical flow techniques[A].Computer Vision and Pattern Recognition[C].Newyork:IEEE,1992:236-242.

[15] 祝彥平.基于互相關(guān)的超聲應(yīng)變成像研究[D].沈陽(yáng):中國(guó)醫(yī)科大學(xué),2012.

[16] 龍琪.基于圖形處理器的超聲仿真平臺(tái)[D].安徽:中國(guó)科技大學(xué),2010.

[17] 孫瑞超,唐亞男,李靈,等.實(shí)時(shí)超聲彈性成像原理與方法[J].中國(guó)醫(yī)療器械信息,2013,(6):1-9.

本文編輯 聶孝楠

Research on Hybrid Ultrasonic Elastography Algorithm Based on Zero Phase and Optical Flow

LI Yuan-rui, ZHANG Yong-de
Department of Biomedical Engineering, China Medical University, Shenyang Liaoning 110000, China

In view of the problem that the traditional time-domain cross-correlation algorithm of ultrasound elastography in estimating tissue displacement compute-intensive and inefficient, this paper proposed a hybrid method of ultrasound elastography based on phase zero estimation and optical flow method. Firstly, the zero phase was used to estimation method to obtain coarse displacement field, as a result of the initial displacement for optical flow; then the classical optical flow method Horn & Schunck (HS) optical flow method was used to correct for the displacement field; finally the least squares method was used to calculate the displacement field obtained axial strain organizations. In this paper, the simulation data algorithm was validated and the impact of simulation data under different compression on the experimental results was discussed. The experimental results showed that the longitudinal displacement diagram and the longitudinal strain map of the zero-phase and the optical flow method proposed in this paper were consistent with the finite element results, and the correctness of the algorithm was proved.

ultrasound elastography; optical flow; phase zero estimation; decorrelation; least squares method

R318.6

B

10.3969/j.issn.1674-1633.2017.06.016

1674-1633(2017)06-0061-05

2016-10-20

2016-12-20

省級(jí)科技項(xiàng)目“基于相關(guān)相位和光流法的二維實(shí)時(shí)超聲彈性成像方法研究”（L2012294）。

張勇德，副教授，主要研究方向?yàn)槌晱椥猿上瘛?/p>

通訊作者郵箱：zyd1977@sina.com