倪 東，陳思平，汪天富

醫(yī)學(xué)超聲關(guān)鍵技術(shù)國家地方聯(lián)合工程實驗室，廣東省生物醫(yī)學(xué)信息檢測與超聲成像重點實驗室，深圳大學(xué)醫(yī)學(xué)院，深圳518060

醫(yī)學(xué)超聲因具有無損、無輻射、實時、便捷等特點，被廣泛應(yīng)用于臨床診斷、治療和介入手術(shù)引導(dǎo)中. 在醫(yī)學(xué)超聲成像研究中，為優(yōu)化超聲探頭設(shè)計，Jensen［1-3］提出一種用于任意形、不等指長換能器的脈沖壓力場的模擬方法，在工業(yè)界和學(xué)術(shù)界均取得巨大成功. 該方法產(chǎn)生的超聲模擬圖像雖然較為逼真，但是往往需要幾個小時才能產(chǎn)生一幅圖像.

近年來，越來越多的研究者開始關(guān)注實時超聲成像模擬，Maggie 等［4-5］提出從真實超聲體數(shù)據(jù)中提取超聲圖像紋理來構(gòu)建超聲模擬圖像. Henry［6］和倪東［7］等提出利用虛擬的超聲探頭，從已獲取的寬視野三維超聲數(shù)據(jù)中獲取二維超聲圖像. 但是，構(gòu)建寬視野三維超聲數(shù)據(jù)需要魯棒的圖像配準(zhǔn)和拼接算法，或者需要利用定位器來確定不同三維超聲數(shù)據(jù)間的位置關(guān)系. Shams 等［8］提出通過合成離線獲取的散射圖和在線獲取的反射圖來獲取超聲圖像，其中散射圖通過Field II 軟件獲取. 這種方法的一個局限性是即使在一個有20 個CPU 的計算機服務(wù)器上，產(chǎn)生散射圖也需約30 h. Song［9］提出用光線跟蹤算法模擬聲束，并用于心臟超聲圖像的模擬. 這種方法用理想化的光線模擬聲束，光線不具有面積，降低了模擬圖像的真實性. 文獻［4-9］的另一個局限在于不能實時調(diào)整超聲成像的重要參數(shù)，如聲波頻率和焦距. 但在實際超聲掃查過程中，選擇最優(yōu)成像參數(shù)是一個重要的技巧.

斑點是由于聲波的隨機散射造成的，是超聲成像的固有特性，在病理診斷中有一定作用. 但由于聲波在人體內(nèi)傳輸?shù)膹?fù)雜性，完全模擬超聲斑點非常耗時. 近年已有研究利用不同的概率分布函數(shù)來描述超聲斑點模式［10-14］. Johan［10］指出斑點可以點擴散函數(shù)來顯示.

本研究提出一種基于曲光線跟蹤算法的超聲成像模擬新方法. 利用曲光線跟蹤算法和基本聲學(xué)模型，產(chǎn)生能量反射圖和衰減圖;用多條曲光線模擬聚焦聲束，通過調(diào)整光線參數(shù)模擬聲束頻率和焦距的變化;運用瑞利概率分布函數(shù)產(chǎn)生超聲斑點模式，并利用高斯型點擴散函數(shù)進行顯示. 實驗結(jié)果表明，該算法能滿足實時模擬超聲成像的要求.

1 聲束聚焦模型

在現(xiàn)代超聲成像系統(tǒng)中，探頭孔徑由一組陣元組成. 未聚焦的換能器形成的波束較寬，導(dǎo)致橫向分辨率較低. 因此，聲束聚焦技術(shù)被廣泛應(yīng)用于現(xiàn)代超聲成像系統(tǒng)中. Christensen［15］指出聲束在近場具有不規(guī)則形狀，難以定義近場聲束模型. 但在遠(yuǎn)場，聲束輪廓隨著距離的變化進行線性傳播，并偏移固定的角度，如圖1(a). 偏移角φd反比于探頭直徑D，正比于聲波波長λ，定義為

文獻［15］還指出聚焦區(qū)域橫向分辨率與聲束寬度相關(guān)，可近似為1.22λ/D. 提高聲波頻率可增加橫向分辨率，超聲圖像質(zhì)量與聲波頻率緊密相關(guān)，而在遠(yuǎn)場，超聲圖像因橫向分辨率的降低而顯得模糊. 基于以上聲波傳播的基本模型，本研究提出曲光線跟蹤模型.

圖1 聲束聚焦及曲光線跟蹤模型Fig.1 Illustration of the ray tracing model

2 曲光線模型

在傳統(tǒng)光線跟蹤模型中，理想光線用單個直線表示. 由于聲束具有一定寬度及聚焦特性，我們提出用多條曲線模擬聲束. 對于探頭發(fā)出的每根掃描線如圖1(b)均用多條曲線來模擬. 定義標(biāo)量n 為曲線的標(biāo)號，T 為探頭掃描線的方向矢量，單位矢量N 為方向矢量的正交矢量，標(biāo)量s 為聲波傳播距離，曲光線的方程定義為

其中，n = 0，±1，±2，…，±L;δn(s)是用于調(diào)整曲線形狀的函數(shù)，可由式(3)計得，

這里，A 是孔徑;F 是焦距;L 是曲線在矢量T 兩側(cè)的最大個數(shù);f 是聲波頻率;β 是用于控制曲線形狀的常量.

以下分析式(2)和式(3)的物理特性. 設(shè)n= L 對應(yīng)于圖1(c)中最外面的一條曲線，當(dāng)s ＜F時，隨著傳播距離s 的增加，通過控制β 值，δn(s)趨于最小值F，曲光線的曲率可由參數(shù)β 進行調(diào)整，與β 值成反比;當(dāng)s ≥F 時，聲束處于遠(yuǎn)場，將式(2)對s 求導(dǎo)數(shù)，求得曲線的斜率，即

因此，本研究提出的曲光線模型在近場是曲線，其寬度隨著傳播距離的增加而減小，而在遠(yuǎn)場是直線，直線的斜率與聲束在遠(yuǎn)場的偏移角的正弦值相同. 在實際應(yīng)用中，聲束偏移角φd通常較小(＜5°)［16］，則

可見曲光線模型與聲束聚焦模型一致. 在式(3)中，當(dāng)s = F 時，δn(s)值最小，此時它就是聲束聚焦處的寬度.

圖2 曲光線跟蹤模型示意圖Fig.2 Illustration of different effects in ray tracing process

3 反射圖模擬

圖3 為為超聲圖像模擬示意圖. 其中，圖3(a)是合成的CT 圖像，圖3(b)是對目標(biāo)光線跟蹤得到的灰度圖M，在此基礎(chǔ)上，引入基本聲學(xué)模型，計算能量反射圖R 如圖3(c)和衰減圖S 如圖3(d).

聲波在傳播過程中，由一種介質(zhì)到達另一種介質(zhì)，在兩種介質(zhì)的分界面上，由于聲阻抗的不同，聲波會發(fā)生方向和能量的變化:一部分聲波被反射回原介質(zhì)中，稱為反射波;另一部分聲波透過界面在另一種介質(zhì)中繼續(xù)傳播，稱為折射波. 這種現(xiàn)象將導(dǎo)致聲波能量的衰減. 由于超聲成像主要利用反射波檢測信號，并且折射波的模擬非常困難和耗時［4］，因此，本研究僅考慮聲波的反射和由此引起的能量衰減.

聲波反射與聲阻抗密切相關(guān). 聲阻抗Z = ρc.其中，ρ 為介質(zhì)密度;c 是聲速. 聲速在人體軟組織中的傳播速度約為1 540 m/s. 設(shè)聲速為常數(shù)，聲阻抗正比于組織密度. 利用文獻［17］提出的CT圖像亨氏單位(Hounsfield unit)與組織密度之間存在映射關(guān)系，可由圖像灰度求取組織聲阻抗，并根據(jù)聲波能量反射系數(shù)公式(Z2-Z1)2/(Z2+Z1)2求得反射圖.

聲波衰減與聲波所經(jīng)過的組織和聲波頻率相關(guān)，文獻［15］指出，聲波能量I 隨傳輸距離的變化呈指數(shù)衰減，即

其中，γ 是與組織性質(zhì)相關(guān)的衰減系數(shù);f 是聲波頻率;τ 是經(jīng)驗常數(shù)，本研究設(shè)τ = 1.

圖3 超聲圖像模擬中間結(jié)果示意圖Fig.3 Intermediate results

4 斑點模式模擬

斑點是超聲圖像由于聲波隨機散射造成的固有特性，在病理診斷中有一定作用. 但要完全從物理機制上模擬斑點模式是非常困難和耗時的，不適用于超聲成像的實時模擬. 瑞利分布是描述超聲斑點模式最常用的概率分布之一，我們采用瑞利概率分布函數(shù)產(chǎn)生超聲斑點模式，并利用高斯型點擴散函數(shù)進行顯示，起到平滑斑點噪聲的效果. 為達到實時模擬的效果，斑點模擬和顯示均用圖形處理器(graphics processing unit，GPU)進行加速.

最終的超聲圖像模擬定義為

其中，S、M 和R 分別為能量衰減圖、灰度圖和反射圖;P 為高斯卷積函數(shù);A 為圖像灰度值的瑞利概率分布值;δx和δy是高斯卷積函數(shù)的參數(shù)，分別對應(yīng)橫向分辨率和縱向分辨率. 超聲圖像模擬示意圖見圖3(e).

5 實驗結(jié)果分析

為驗證算法的可行性，我們對合成圖像、解剖結(jié)構(gòu)模型以及CT 圖像進行了超聲成像模擬，并比較了用不同的聲波頻率和聚焦位置獲得的模擬結(jié)果. 圖4 中的白色箭頭表示焦距的位置.

首先，將自行合成的體模作為超聲成像模擬的輸入，圖4(a)中(i)是分辨率體模圖像，圖4(b)中(i)是另一種體模圖像. 圖4(a)其他圖像采用不同的聲波頻率對輸入體模進行成像:(ii)為0.5 MHz，(iii)為1 MHz，(iv)為3 MHz. 結(jié)果顯示，隨著頻率的增加，圖像的橫向分辨率得到提高，但由于高頻率導(dǎo)致聲波傳輸能量的衰減增多，灰度逐步變小. 圖4(b)采用不同的聲波頻率和焦距進行成像:(ii)為0.5 MHz，1.5 cm;(iii)為0.5 MHz，3 cm;(iv)為1 MHz，3 cm. 從結(jié)果可見，隨著頻率的升高圖像的橫向分辨率得到提高，且在焦點位置的圖像分辨率高于其他位置的圖像.

圖4 合成體模超聲成像模擬圖像Fig.4 Effect of acoustic frequency and beam focus on phantom data

為進一步說明本算法的可行性，將本研究結(jié)果與超聲成像模擬的權(quán)威軟件——Field II 生成的結(jié)果進行對比. 在文獻［1］ 中介紹了兩種超聲體模:腎臟體模和胎兒體模. 由于文獻［1］將這兩種體模描述為散射子的分布，我們首先基于解剖結(jié)構(gòu)信息將散射子分布圖轉(zhuǎn)為灰度圖，并將其作為超聲成像模擬的輸入. 圖5(a)和圖6(a)是Field II 軟件的生成結(jié)果，在圖5(b)中，由于聲束聚焦位置在上端，下面的圖像略顯模糊，而圖5(c)中，由于聲束聚焦位置在下端的圖像區(qū)域，則顯示了較好的細(xì)節(jié)，圖像也較清晰. 圖(6)顯示的是胎兒體模的超聲模擬圖，由于焦距的不同，圖6(b)和圖6(c)中的胎兒手指清晰度有較大差異.

圖5 腎臟體模超聲成像模擬圖Fig.5 Comparison of simulated kidney image

圖6 胎兒體模超聲成像模擬圖Fig.6 Comparison of simulated fetus image

最后，我們將模擬結(jié)果與真實的超聲圖像進行了對比. 把一個實際的肝臟體模(CIRS Model 057)分別進行CT 成像和超聲3 維成像，并將超聲三維圖像與CT 圖像進行配準(zhǔn). 選取CT 圖像的1 個切面作為超聲成像模擬的輸入，并將模擬結(jié)果與配準(zhǔn)后相同位置的真實超聲圖像進行對比. 圖7(a)是1幅CT 圖像，圖7(b)是用本算法產(chǎn)生的模擬結(jié)果，圖7(c)是用超聲掃描儀獲取的真實超聲圖像，兩者之間的相似性說明本算法可行.

另外，我們還對模擬成像的速度進行了評估.模擬軟件在臺式計算機上進行測試(CUP 為Penium 4 dual core 2，2.8 GHz;內(nèi)存為2 Gbit;GPU 為Nvidia 8800). 通過采用GPU 對斑點的生成和顯示進行加速，對于128 ×300 (128 根掃描線，每線300 個采樣點)的超聲圖像，模擬成像速度可達約20 幀/s. 因此，本研究在超聲引導(dǎo)下的虛擬手術(shù)和計劃，超聲診斷培訓(xùn)等方面有較好應(yīng)用前景.

圖7 肝臟體模超聲成像模擬圖Fig.7 Simulated ultrasound image from CT volume

結(jié) 語

本研究提出一種實時模擬超聲成像的新方法，使用曲光線跟蹤模型模擬聲束，可實時調(diào)整聲波的頻率和焦點位置，同時利用概率分布函數(shù)產(chǎn)生超聲斑點，并用GPU 進行加速顯示. 將本研究算法產(chǎn)生的結(jié)果與Field II 軟件和實際超聲掃描儀獲得的結(jié)果進行了對比，結(jié)果的相似性說明了本研究算法的可行性，本算法的速度也達到了實時性的要求.下一步，我們將基于聲學(xué)理論，研究斑點噪聲的概率分布模式、聲波散射和折射現(xiàn)象的模擬，并研究臨床病例CT 數(shù)據(jù)的超聲成像模擬，最終將研究成果用于超聲診斷培訓(xùn)、多模態(tài)圖像配準(zhǔn)以及超聲引導(dǎo)下的虛擬手術(shù)等領(lǐng)域.

/References:

［1］Jensen J A. Field:a program for simulating ultrasound systems ［J］. Medical ＆ Biological Engineering ＆ Computing，1996，34(1):351-353.

［2］Jensen J A，Gandhi D，O'Brien W D Jr. Ultrasound fields in an attenuating medium ［C］// Proceedings of the IEEE Ultrasonics symposium. Vancouver (Canada):IEEE Press，1993，2:943-946.

［3］Jensen J A，Nikolov S. Fast simulation of ultrasound images ［C］// Proceedings of the IEEE Ultrasonics symposium.Vancouver (Canada):IEEE Press，2000:1721-1724.

［4］Magee D，Zhu Y，Ratnalingam R，et al. An augmented reality simulator for ultrasound guided needle placement training ［J］. Journal of Shenzhen University Science and Engineering，2007，45(10):957-967.

［5］Zhu Y，Magee D，Ratnalingam R，et al. A training system for ultrasound-guided needle insertion procedures［C］// Medical Image Computing and Computer-Assisted Intervention.Berlin:Springer，2007:566-574.

［6］Henry D，Troccaz J，Bosson J，et al. Ultrasound imaging simulation:application to the diagnosis of deep venous thromboses of lower limbs ［C］// Medical Image Computing and Computer-Assisted Intervention. Berlin:Springer，1998:1032-1040.

［7］NI Dong，CHAN Wing-ying，QIN Jing，et al. A virtual reality simulator for ultrasound guided organ biopsy training［J］. IEEE on Computer Graphics and Application，2011，31(2):36-48.

［8］Shams R，Hartley R，Navab N. Real-time simulation of medical ultrasound from CT images ［C］// Medical Image Computing and Computer-Assisted Intervention. Berlin:Springer，2008:734-741.

［9］Song M，Haralick R，Sheehan F. Ultrasound imaging simulation and echocardiographic image synthesis ［C］//Proceedings of International Conference on Image.Vancouver:IEEE Press，2000:420-423.

［10］Johan M T. Speckle formation，analysis and processing applied to ultrasound tissue characterization ［J］. Physics for Medical Imaging Applications，2007，240(3):151-176.

［11］Palhano Xavier de Fontes，Guillermo Andrade Barroso，Pierrick Coupé，et al. Real time ultrasound image denoising ［J］. Journal of Real-time Image Processing，2011，6(1):15-22.

［12］Goodman J W. Some fundamental properties of speckle［J］. Journal of the Optical Society of America B，1976，66(11):1145-1150.

［13］CHEN Si-ping，ZENG Si-ning，WANG Tian-fu，et al.Numerical simulation on acoustic radiation force induced by ultrasonic coded excitation ［J］. Journal of Shenzhen University Science and Engineering，2011，28(2):165-171.(in Chinese)陳思平，曾斯寧，汪天富，等. 超聲編碼激勵產(chǎn)生聲輻射力的數(shù)值模擬［J］. 深圳大學(xué)學(xué)報理工版，2011，28(2):165-171.

［14］Ogier A，Hellier P. A modified total variation denoisting method in the context of 3D ultrasound images ［C］//Medical Image Computing and Computer-Assisted Intervention.Berlin:Springer，200:70-77.

［15］Christensen D A. Ultrasound Bioinstrumentation ［M］.1st ed. New York:John Wiley ＆ Sons Inc，1988:89-112.

［16］Angelsen B A J，Waag R C. Ultrasound imaging:waves，signals，and signal processing ［J］. Journal of the Acoustical Society of America，2007，121(4):1820.

［17］Wein W，Brunke S，Khamene A，et al. Automatic CTultrasound registration for diagnostic imaging and imageguided intervention ［J］. Medical Image Analysis，2008，12(5):577-585.