一種超聲彈性成像試驗數(shù)據(jù)的快速仿真算法

崔少國 LIU Pauls 劉東權(quán)(四川大學(xué)計算機學(xué)院，四川成都 610065)

2(鄖陽醫(yī)學(xué)院計算機中心，湖北十堰 442000)

3(聲泰特成都科技有限公司，四川成都 610041)

引言

超聲彈性成像是一種生物組織成像的新模式，是當(dāng)前超聲醫(yī)學(xué)工程學(xué)領(lǐng)域的前沿和熱點。該成像技術(shù)可以探測出人體組織與病理變化相關(guān)的硬度改變，對腫瘤(癌癥)的早期檢測和診斷具有重要意義。

彈性成像的研究自1991 年被Ophir 等人提出以來，19 年中取得了快速發(fā)展，在臨床中已逐漸顯示出重要診斷價值，比如在乳腺癌、前列腺癌、甲狀腺結(jié)節(jié)、肝纖維化、動脈粥樣硬化等的檢測中取得了重要的應(yīng)用［1-4］。其成像的基本原理是:對組織施加一個動態(tài)或者準(zhǔn)靜態(tài)的壓力，組織將產(chǎn)生一個響應(yīng)。在相同的應(yīng)力下硬的組織形變較小，而軟的組織形變較大。使用信號處理和圖像處理相關(guān)技術(shù)計算出組織的應(yīng)變大小，進而重構(gòu)出組織的彈性模量。對組織的應(yīng)變或彈性模量進行成像就是彈性成像［2］。

1 彈性成像試驗數(shù)據(jù)仿真

在彈性成像算法設(shè)計、性能評估和成像參數(shù)選擇的研究中，由于從實際超聲系統(tǒng)中獲取的數(shù)據(jù)具有事先不可預(yù)知的噪聲或其它干擾因素，所以一種新算法的驗證通常都需要在仿真的數(shù)據(jù)上進行預(yù)實驗，然后再進行體模實驗和在體實驗。彈性成像試驗數(shù)據(jù)的仿真在超聲彈性成像研究中具有重要的意義。

當(dāng)前，超聲彈性成像試驗數(shù)據(jù)的仿真一般采用有限元分析(FEA)和Field Ⅱ聲場模擬程序來實現(xiàn)。這種方法把組織仿真成在某種基質(zhì)中分布著與基質(zhì)聲學(xué)特性有顯著差異的散射子，各散射子的位置服從均勻隨機分布，而其反射強度服從高斯正態(tài)分布，從而建立計算機人體模型［5］(以下簡稱計算機體模)。對于壓縮前回波信號，依據(jù)線性聲學(xué)理論，使用Field Ⅱ程序在仿真的計算機體模上生成。對于壓縮后回波信號，先使用有限元分析出各散射子在一定應(yīng)力作用下運動到達的新位置，使用它對計算機體模中散射子位置進行重新調(diào)整，然后再使用Field Ⅱ程序生成壓縮后回波信號。該仿真方法與實際成像過程很相似，結(jié)果也較準(zhǔn)確，但仿真數(shù)據(jù)生成速度較慢。一般來說，仿真一幀射頻回波信號需要6 小時左右。在彈性成像研究中，通常需要仿真不同組織模型、不同應(yīng)變量和不同信噪比的序列回波數(shù)據(jù)幀，數(shù)據(jù)量較大，使用該方法來仿真需要較長時間。

2 一種快速仿真算法

筆者在彈性成像研究過程中探討出了一種快速仿真彈性成像試驗數(shù)據(jù)的算法，該算法主要包括以下三步:參考幀數(shù)據(jù)生成，壓縮后回波數(shù)據(jù)生成和數(shù)據(jù)加噪處理。

2.1 參考幀數(shù)據(jù)生成

參考幀就是壓縮過程中采集的序列幀中的第一幀。對參考幀回波數(shù)據(jù)的仿真可以使用兩種方法:基于計算機體模和Field Ⅱ的方法和使用超聲線性系統(tǒng)理論的方法。實際應(yīng)用中可以選擇其中任一種方法。

2.1.1 使用計算機體模和Field Ⅱ聲場模擬程序

計算機體模使用文獻［4］中介紹的方法生成，采用組織分立散射模型，各散射子的位置服從均勻隨機分布，而其反射強度服從高斯正態(tài)分布。聲場計算使用Tupholme 和Stepanishen 提出的算法。

探頭在聲場中位置為→r 處產(chǎn)生的聲壓被描述為空間沖擊響應(yīng)［6］，表示如下:

這里，δ 是沖擊函數(shù)，c 是超聲波在組織中傳播速度，s 是探頭的表面積。那么探頭發(fā)射波的聲壓場為［7］

這里Vn(t)是探頭的表面速度，ρ 是組織密度，“* ”代表卷積。

因此，從探頭接收到的超聲回波信號可模擬為［8］

式中，Vpe是回波沖擊響應(yīng)，fm項是由于組織的密度變化和聲速擾動引起的散射信號，hpe是脈沖回波的沖擊響應(yīng)，分別為

2.1.2 使用超聲線性系統(tǒng)理論

使用超聲簡單線性系統(tǒng)理論，參考幀回波信號可以用如下的公式產(chǎn)生:

這里f(x，y)是射頻回波信號，p(x，y)是系統(tǒng)點擴散函數(shù)(PSF)，h(x，y)是散射函數(shù)，散射子的位置服從均勻隨機分布，散射子的散射強度服從高斯分布，“* ”代表卷積。二維點擴散函數(shù)p(x，y)在軸向(x方向)上可以模擬成高斯調(diào)制的余弦函數(shù)，在側(cè)向(y 方向)上可以模似成一個高斯包絡(luò)，表示如下

其中，ω0是探頭的中心角頻率，σx和σy分別是軸向和側(cè)向的標(biāo)準(zhǔn)方差。散射函數(shù)h(x，y)可以用服從均勻分布的隨機函數(shù)生成。

2.2 壓縮后回波數(shù)據(jù)生成

如果組織內(nèi)各區(qū)域彈性模量一致，當(dāng)在相同應(yīng)力作用下，組織各部分發(fā)生相同的應(yīng)變，這種稱為均勻應(yīng)變模型。此模型主要用來仿真均勻的無腫瘤人體組織，用來研究彈性成像算法的信噪比SNRe。如果組織中某些區(qū)域彈性模量與周圍區(qū)域不同，當(dāng)在均勻應(yīng)力作用下，不同彈性模量的組織將產(chǎn)生不同的應(yīng)變，這種模型稱為非均勻應(yīng)變模型。此模型主要用來仿真均勻組織中存在一硬癌塊或軟囊腫，用來研究彈性成像算法的對比度(contrast)或?qū)Ρ榷仍肼暠菴NRc。

算法1 是仿真組織壓縮后回波信號的快速生成算法。設(shè)均勻組織各點的彈性模量m 是常數(shù)C，非均勻組織中各點的彈性模量m 是e(i，j)(i，j 是組織區(qū)域二維空間坐標(biāo))，模量大小隨空間位置坐標(biāo)(i，j)而變化。設(shè)組織為平面應(yīng)變狀態(tài)，組織的各部分均勻受力，各點的應(yīng)力為F。根據(jù)彈性力學(xué)及胡克定律，各點的應(yīng)變

根據(jù)線彈性體中位移連續(xù)性原則，則各點的位移d(i，j)為

設(shè)RF_pre 和RF_post 是組織壓縮前和壓縮后的回波RF 信號，IQ_pre 和IQ_post 是其對應(yīng)的基帶I/Q信號，則算法可以簡單描述如下:

算法1:組織壓縮后回波信號快速生成算法

1)設(shè)計預(yù)仿真組織的彈性模量分布圖。對于均勻組織各區(qū)域設(shè)成相同的模量，即m(i，j)=C(C是常數(shù));對于非均勻組織根據(jù)預(yù)仿真的模型，不同的空間位置設(shè)成不同的彈性模量，即

2)根據(jù)式(7)由彈性模量圖計算組織應(yīng)變分布圖。

3)根據(jù)式(8)由應(yīng)變圖計算各點在應(yīng)力作用下的位移。

4)根據(jù)位移圖求出插值位置圖，使用式p(i，j)= d(i，j)+i 來完成。

5)使用式(9)對壓縮前信號在新的位置上進行插值，生成壓縮后RF 回波信號RF_post。

6)使用希爾伯特變換生成復(fù)信號，Analytic_post=Hilbert(RF_post)。

7)生成基帶I/Q 信號，IQ_post = Analytic_post×exp(j ×w0 ×p)。

8)如果還需生成壓縮過程中的下一幀數(shù)據(jù)，則RF_pre =RF_post，轉(zhuǎn)4);否則結(jié)束。

這里插值方法可以選用計算代價較小的線性插值，也可以選用計算代價較大但較精確的三次樣條插值。這里采用線性插值技術(shù)，插值函數(shù)如下:

2.3 數(shù)據(jù)加噪處理

為了仿真真實超聲系統(tǒng)中獲取的不同信噪比的回波數(shù)據(jù)，需要對以上產(chǎn)生的信號進行加噪處理。

根據(jù)信噪比的定義

式中，Ps是信號的功率，Pu是噪聲的功率，則

因此，所需要加的高斯白噪聲為

這里，Rand(m，n)是產(chǎn)生m 行n 列均值為0 方差為1 的高斯隨機信號值函數(shù)。加噪算法可簡單描述如下:

算法2:數(shù)據(jù)加噪算法

1)設(shè)定數(shù)據(jù)的信噪比SNR，計算原信號的功率Ps。

2)根據(jù)式(10)計算需加的噪聲功率Pu。

3)根據(jù)式(11)計算所需的高斯白噪聲u(m，n)。

4)對仿真的回波數(shù)據(jù)f(i，j)加噪處理，fnoise(i，j)= f(i，j)+ u(i，j);

5)結(jié)束。

3 算法驗證

為了驗證以上提出的快速仿真算法的正確性和高效性，筆者編程實現(xiàn)了上述仿真算法，并生成了一系列數(shù)據(jù)集。然后選用一標(biāo)準(zhǔn)超聲彈性成像程序，在仿真的數(shù)據(jù)集上運行，以驗證能否產(chǎn)生預(yù)期的應(yīng)變圖及成像效果。

3.1 參考幀生成

實驗對兩種生成參考幀的方法都進行了驗證。計算機體模是仿真40 mm ×40 mm ×10 mm 并均勻分布有200 000 個散射子的人體組織，探頭中心頻率5 MHz，系統(tǒng)采樣頻率是20 MHz，使用Field Ⅱ和式(1)～式(4)生成128 條A-lines，幀大小128 像素×1039 像素。使用超聲線性系統(tǒng)理論產(chǎn)生參考幀時，使用式(5)和式(6)。

3.2 應(yīng)變模型設(shè)計及壓縮后回波數(shù)據(jù)生成

本實驗中預(yù)仿真兩種組織應(yīng)變模型:均勻組織應(yīng)變模型和軟背景中包含一球形硬包容物非均勻組織應(yīng)變模型，模型的彈性模量分布如圖1 所示(彈性模量用256 級灰度值表示，平均模量映射成128)。圖1(a)是均勻組織彈性模量分布圖，組織內(nèi)各點的彈性模均是30 kPa。圖1(b)是非均勻組織應(yīng)變模型的彈性模量分布圖，圖中硬物的直徑為20 mm，其彈性模是300 kPa，背景的彈性模量是30 kPa，硬物是背景硬度的10 倍。使用算法1 對每種模型產(chǎn)生20 幀應(yīng)變序列。

實驗針對均勻組織應(yīng)變模型仿真不同整體應(yīng)變的回波序列，包括0.01%，0.1%，0.5%，1%，2%，3%，4%，5%，用來驗證算法產(chǎn)生的數(shù)據(jù)是否滿足彈性成像應(yīng)變?yōu)V波器。對非均勻模型仿真了背景應(yīng)變?yōu)?%、硬物應(yīng)變?yōu)?.1%的數(shù)據(jù)。

圖1 預(yù)仿真的組織模型彈性模量分布圖。(a)均勻組織;(b)含一硬包容物的非均勻組織Fig.1 Simulated tissue modulus distributing images.(a)uniform model;(b)one-lesion model

3.3 數(shù)據(jù)加噪

使用算法2 對產(chǎn)生的各幀仿真數(shù)據(jù)添加均值為零的不同能量高斯白噪聲，使其信噪比分別為5、10、15、20、25、30、35、40、45、50 dB。

3.4 成像算法的選擇

實驗中選用英國劍橋大學(xué)Joel. E Lindop 于2008 年提出的基于相位的WPS(weighted phase separation)實時彈性成像算法［9］來驗證仿真的數(shù)據(jù)，該算法已被驗證是正確可行的，能在實際超聲系統(tǒng)上進行彈性成像。位移估計使用式(12)迭代進行

其中，WA(t)= Ax(t)Ay(t + dk)，

應(yīng)變計算使用低通數(shù)字差分濾波器SG-I DD［10］，方法為

式中，s 為應(yīng)變，Δk= d(n + k)- d(n - k)

3.5 時間效率和數(shù)據(jù)仿真效果檢測方法

為了對新算法產(chǎn)生彈性成像試驗數(shù)據(jù)的時間效率和仿真效果進行檢測，本研究中使用了三種方法產(chǎn)生彈性成像試驗數(shù)據(jù)并進行比較。方法1 是:參考幀使用計算機體模和Field Ⅱ生成，壓縮后的回波數(shù)據(jù)先使用有限元分析，進而對體模中的散射子的位置進行調(diào)整，然后使用Field Ⅱ產(chǎn)生回波，數(shù)據(jù)加噪使用算法2;方法2 是:參考幀使用計算機體模和Field Ⅱ生成，壓縮后的回波數(shù)據(jù)使用提出的算法1、算法2 生成;方法3 是:參考幀使用簡單線性系統(tǒng)理論卷積生成，壓縮后的回波數(shù)據(jù)使用提出的算法1、算法2 生成。運行算法程序所使用的計算 機 配 置 是:Pentium (R) 4，CPU 2.4 GHz，RAM 1GB。

彈性成像試驗數(shù)據(jù)效果的比較，主要以在該數(shù)據(jù)集上產(chǎn)生的對應(yīng)彈性圖像的信噪比來衡量。信噪比越高說明仿真與實際越接近，仿真越準(zhǔn)確。該項實驗使用均勻應(yīng)變模型的數(shù)據(jù)，數(shù)據(jù)沒有進行加噪處理。

4 結(jié)果

4.1 應(yīng)變圖像及性能

圖2 是在仿真的均勻組織應(yīng)變模型的數(shù)據(jù)集上使用WPS 算法產(chǎn)生的應(yīng)變圖。成像中窗長為60 個采樣點，窗口重疊率為0.75，迭代3 次。運動追蹤在1D 軸向進行，軸向應(yīng)變計算采用M =4 的低通數(shù)字差分濾波器SG-I DD 生成。圖2(a)是由仿真的1%應(yīng)變的無噪聲數(shù)據(jù)集產(chǎn)生，圖2(b)是由仿真的1%應(yīng)變且信噪比為10 dB 的數(shù)據(jù)集產(chǎn)生。應(yīng)變映射成256 級灰度值，平均應(yīng)變映射成128(下同)。

圖2 仿真的均勻組織應(yīng)變模型數(shù)據(jù)集上產(chǎn)生的應(yīng)變圖像。(a)無噪聲;(b)SNR =10 dBFig. 2 Elastography obtained from the simulated echoes of uniform strain profile. (a)no noise;(b)SNR =10 dB

圖3 是在仿真的非均勻組織應(yīng)變模型(含一硬包容物)的數(shù)據(jù)集上使用WPS 算法產(chǎn)生的應(yīng)變圖。成像參數(shù)與方法同圖2。圖3(a)是由仿真的1%整體應(yīng)變的無噪聲數(shù)據(jù)集產(chǎn)生，圖3(b)是由仿真的1%整體應(yīng)變且信噪比為10 dB 的數(shù)據(jù)集產(chǎn)生。

圖3 仿真的非均勻應(yīng)變模型數(shù)據(jù)集上產(chǎn)生的應(yīng)變圖像。(a)無噪聲;(b)SNR =10 dBFig. 3 Elastography obtained from the simulated datum of one-lesion strain profile. (a)no noise;(b)SNR =10 dB

從圖2 和圖3 可以看出，在用提出的快速算法所仿真的數(shù)據(jù)集上可以生成預(yù)期的應(yīng)變圖像。圖2產(chǎn)生的應(yīng)變圖像與預(yù)仿真的均勻組織應(yīng)變模型相符，圖3 也與預(yù)仿真的的非均勻組織應(yīng)變模型相符(較硬的圓形區(qū)域里應(yīng)變約為背景的1/10)。在無或有噪聲的仿真數(shù)據(jù)上采用相同的成像算法，和成像參數(shù)所生成的應(yīng)變圖有明顯差別，在有噪聲的仿真數(shù)據(jù)上產(chǎn)生的應(yīng)變圖噪聲較大，這與彈性成像的實際情況相符。

圖4 顯示了應(yīng)變圖像中最中間的一條線(軸向)上的應(yīng)變數(shù)據(jù)。從圖4 中也可以看出，在仿真的無噪聲數(shù)據(jù)集和加噪聲數(shù)據(jù)集上產(chǎn)生的應(yīng)變數(shù)據(jù)是有效的，與彈性成像理論相符。

圖5 顯示了在仿真的各種應(yīng)變數(shù)據(jù)集上生成的應(yīng)變圖像的信噪比。從圖中可以看出，信噪比曲線滿足應(yīng)變?yōu)V波器特性，類似于一帶通濾波器。而且在低信噪比(10 dB)的數(shù)據(jù)集上生成的應(yīng)變圖像信噪比明顯低于無噪聲數(shù)據(jù)集。

4.2 時間效率與仿真效果比較

各種產(chǎn)生彈性成像試驗數(shù)據(jù)方法的時間性能如表1。從表1 中可以看出，采用基于計算機體模和Field II 方法，產(chǎn)生1 幀形變數(shù)據(jù)大約需要6.3 h，而采用新方法(算法1)，產(chǎn)生1 幀相同的數(shù)據(jù)只需要約35 s，速度大大提高。

表1 3 種方法產(chǎn)生彈性成像試驗數(shù)據(jù)所需時間(h/幀)Tab.1 Comparison time of three methods for producing strain echoes(h/frame)

各種方法產(chǎn)生的彈性成像試驗數(shù)據(jù)的效果見表2。表2 顯示，在新方法(方法3)產(chǎn)生的試驗數(shù)據(jù)集上所獲得的應(yīng)變圖像的信噪比略低于傳統(tǒng)的基于有限元和Field II 方法(方法1)，尤其在大應(yīng)變時較為明顯。方法2 和方法3 的信噪比幾乎沒有差別。

表2 3 種方法數(shù)據(jù)集上產(chǎn)生的彈性圖像信噪比比較Tab.2 Comparison of SNR obtained from data by using three methods

圖4 一條應(yīng)變數(shù)據(jù)線。(a)均勻應(yīng)變;(b)10 dB 信噪比數(shù)據(jù)的均勻應(yīng)變;(c)含一硬包容物的應(yīng)變;(d)10 dB 信噪比數(shù)據(jù)且含一硬包容物的應(yīng)變Fig.4 A strain line. (a)uniform strain;(b)uniform strain with a 10 dB SNR;(c)one-lesion strain;(d)onelesion strain with a 10 dB SNR

5 討論和結(jié)論

圖5 各種不同應(yīng)變的仿真數(shù)據(jù)集上產(chǎn)生的應(yīng)變圖像信噪比Fig.5 SNRe of strain images obtained from various strain simulation data

本研究提出的一種超聲彈性成像試驗數(shù)據(jù)的仿真算法是可行的、高效的，能靈活快速地生成各種應(yīng)變、各種信噪比、各種組織模型的仿真超聲回波數(shù)據(jù)。從成像的可行性來看，使用新算法產(chǎn)生的試驗數(shù)據(jù)能夠生成預(yù)期的均勻組織的彈性圖像和含有硬包容物組織的彈性圖像;從時間效率來看，新算法產(chǎn)生彈性成像試驗數(shù)據(jù)比Field Ⅱ方法快得多;從產(chǎn)生的試驗數(shù)據(jù)效果來看，與使用Field Ⅱ方法產(chǎn)生的數(shù)據(jù)相比除了在較大應(yīng)變時略低外，其它幾乎無差別。

新方法的提出對彈性成像仿真研究具有重要意義，為彈性成像算法驗證、成像各參數(shù)的選擇等研究提供了重要手段，為超聲彈性成像試驗數(shù)據(jù)仿真提供了一種新的高效快捷、方便可行的方法。該方法也能對任何復(fù)雜組織進行仿真，此時只需要設(shè)計好復(fù)雜組織的彈性量分布圖即可。該方法的局限性是對大應(yīng)變的試驗數(shù)據(jù)仿真較不準(zhǔn)確，這主要是因為當(dāng)應(yīng)變較大時，使用有限元方法對散射子位置進行估計要比使用插值法精確。但在實際彈性成像中一般應(yīng)用的應(yīng)變不會太大(一般小于0.01)，因為太大會引起嚴(yán)重的解相關(guān)，使當(dāng)前主流時延估計算法失效。在低應(yīng)變時新算法所獲得的數(shù)據(jù)效果略微偏低，但仍能實行有效的仿真任務(wù)，然而時間效率的大幅度提高是值得肯定的。另外一個局限性是此算法僅適用于1D 應(yīng)變模型。對于2D 應(yīng)變模型，必須考慮組織散射子在應(yīng)力作用下的三維運動。2D 應(yīng)變模型的快速仿真算法的研究是一個有價值的課題。

本研究提出了一種彈性成像試驗數(shù)據(jù)仿真的快速算法。通過成像實驗、應(yīng)變數(shù)據(jù)分析、信噪比實驗、時耗和數(shù)據(jù)效果比較等證實，用該算法生成的仿真數(shù)據(jù)與預(yù)期目標(biāo)相符，并且生成數(shù)據(jù)的速度大大提高。該算法的提出為彈性成像算法的驗證、性能的比較、參數(shù)優(yōu)化、成像效果的評價提供了有力工具，也為缺少真實超聲實驗環(huán)境的研究人員進行彈性成像前期研究提供了重要手段。

［1 ］ Ophir J，éspedes I，Ponnekanti H，et al． Elastography:a quantitative method for imaging the elasticity of biological tissue［J］. Ultrasonic Imaging，1991，13(3):111 – 134.

［2 ］ Céspedes I，Ophir J，Ponnekanti H，et al． Elasticity imaging using ultrasound with application to muscle and breast in vivo［J］. Ultrasonic Imaging，1993，15(2):73 –88.

［3 ］ Righetti R. Elastographic characterization of HIFU-induced lesions in canine livers［J］. Ultrasound in Medicine and Biology，1999，25(7):1099 –1113.

［4 ］ Konofagou E，Ophir J. Myocardial elastography—A feasibility study［J］. Ultrasound in Medicine and Biology，2002，28(4):1113 -1124.

［5 ］ Jensen JA. Computer phantoms for simulating ultrasound B-mode and CFM images［J］. Acoustical Imaging Symposium，1997，23:75 -80.

［6 ］ Jensen JA. Calculation of pressure fields from arbitrarily shaped，apodized，and excited ultrasound transducers ［J ］. IEEE Transactions on Ultrasonics， Ferroelectrics， and Frequency Control，1992，39(5):262 –267.

［7 ］ Jensen JA. Fast Simulation of Ultrasound Images［J］. IEEE Ultrasonics Symposium，2000，16:1721 -1724.

［8 ］ Jensen JA. Field:A program for simulating ultrasound systems［J］. Med. Biol. Eng. Comp.，1996，4:351 -353.

［9 ］ Lindop JE，Treece GE，Gee A，et al． Phase-based ultrasonic deformation estimation［J］. IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics，and Frequency Control，2008，55 (1):94–111.

［10］ Luo Jianwen，Bai Jing，He Ping，et al． Axial Strain Calculation Using a Low-Pass Digital Differentiator in Ultrasound Elastography ［J ］. IEEE Transactions on Ultrasonics，F(xiàn)erroelectrics，and Frequency Control，2004，51 (9):1119 –1127.