朱 翔 田 超

(中國(guó)科學(xué)技術(shù)大學(xué)工程科學(xué)學(xué)院 合肥 230027)

0 引言

超聲層析成像是一種根據(jù)穿透物體聲波信號(hào)反演物體內(nèi)部聲學(xué)特性的成像技術(shù)，其可以分為反射型層析成像和透射型層析成像[1?2]。反射型層析成像通過(guò)探測(cè)被組織反射的聲波信號(hào)，可以重建組織的阻抗特性。透射型層析成像通過(guò)探測(cè)穿透組織的聲波信號(hào)，可以反演出組織的聲學(xué)特性如聲速度和聲衰減。1974年，Greenleaf等[3]利用透射型超聲層析成像實(shí)現(xiàn)了組織聲學(xué)吸收的定量表征。超聲層析成像提供的生物組織聲速分布信息可以幫助醫(yī)生區(qū)分正常組織和病灶部位[4?6]，在乳腺癌早期診斷等方面有重要作用[7?8]。

生物組織聲速圖像的重建通常利用透射超聲波信號(hào)，如圖1所示。環(huán)形超聲換能器包圍待成像生物組織，換能器陣元依次發(fā)射超聲波信號(hào)，對(duì)面換能器陣元接收穿過(guò)組織的超聲信號(hào)，并進(jìn)行數(shù)據(jù)分析和反演重建[2,9?10]。除利用透射信號(hào)外，也有同時(shí)結(jié)合反射信號(hào)和透射信號(hào)進(jìn)行組織聲速分布重建的相關(guān)研究[11?14]。目前超聲層析成像圖像重建方法主要分為兩類，分別為基于射線模型[15?18]和基于全波模型[19?23]的重建算法?；谏渚€模型的重建算法利用超聲波的高頻近似[24]，忽略了超聲傳播過(guò)程中可能的衍射和散射過(guò)程，認(rèn)為聲波在組織內(nèi)近似沿射線傳播?；谌Ｐ偷闹亟ㄋ惴ㄍㄟ^(guò)求解波動(dòng)方程來(lái)代替射線模型的高頻近似，前向過(guò)程考慮到了聲波的反射、衍射以及散射等物理現(xiàn)象，可以得到較射線類重建算法更優(yōu)的圖像分辨率。雖然基于射線模型的聲速圖像分辨率低于全波模型，但由于其數(shù)學(xué)原理簡(jiǎn)單，計(jì)算速度相對(duì)較快，在實(shí)際臨床應(yīng)用方面具有重要價(jià)值[25?28]。此外，基于射線模型的重建結(jié)果可以作為全波重建的初始值[21,29]，可以加快全波重建方法的收斂速度。

圖1 透射超聲層析成像原理圖Fig.1 Principle of ultrasound transmission tomography

本文將首先介紹基于射線模型的超聲層析成像前向過(guò)程，然后是基于射線模型的反問(wèn)題求解方法，主要包括聲波第一到達(dá)時(shí)間提取和圖像重建反問(wèn)題數(shù)學(xué)模型的建立和求解。

1 射線前向模型

聲波傳播的射線模型可以分為直線模型和彎曲射線模型兩種。直線模型忽略了聲波傳播過(guò)程中的折射行為，認(rèn)為聲波完全沿直線傳播。彎曲射線模型考慮了聲波在不同界面處的折射過(guò)程，模型相對(duì)更加精確。本節(jié)將分別從直線模型和彎曲射線模型兩個(gè)方面進(jìn)行展開(kāi)討論。

1.1 直線模型

在直線模型的假設(shè)條件下，超聲層析成像中的前向過(guò)程與計(jì)算機(jī)斷層成像(Computed Tomography,CT)中的前向模型類似，如圖2所示。其中，由發(fā)射器出射聲波在接收器上的投影可以寫為雷登變換形式：

圖2 超聲層析成像中聲波直線投影模型Fig.2 Straight ray projection model in ultrasound transmission tomography

式(1)中，p(ρ,θ)為射線在角度θ情況下的聲波第一到達(dá)時(shí)間，其中ρ表示坐標(biāo)原點(diǎn)到射線的距離，s(x,y)為待測(cè)物體在空間坐標(biāo)(x,y)下的聲慢度(即聲速度的倒數(shù))分布，δ(·)為狄拉克沖激函數(shù)。實(shí)際超聲換能器陣元近似以扇形束發(fā)射，可以先將射線的扇形束坐標(biāo)(β,γ)轉(zhuǎn)換成平行束坐標(biāo)(ρ,θ)，即ρ=Rsin(γ)和θ=β ?γ，然后再利用上述雷登變換進(jìn)行計(jì)算。

直線模型可以近似描述聲波的前向過(guò)程，但最終重建出來(lái)的圖像在形狀和數(shù)值上會(huì)存在較大誤差，主要是因?yàn)橹本€模型忽略了聲波在界面處的折射過(guò)程，最終導(dǎo)致射線積分路徑錯(cuò)誤。為此，可以引入彎曲射線模型校正直線模型存在的誤差。

1.2 彎曲射線模型

近年來(lái)，基于彎曲射線模型的迭代重建算法不斷發(fā)展，通過(guò)追跡聲波的折射路徑可以較好地校正直線模型引入的誤差。相關(guān)研究表明彎曲射線追跡算法能夠補(bǔ)償介質(zhì)中高達(dá)20% 的聲速變化[30?31]。然而在典型的乳腺癌診斷中，乳腺中不同組織的聲速變化大概為8%[32]。因此，可以利用彎曲射線模型較準(zhǔn)確地描述組織聲速差異引起的聲波折射行為。

準(zhǔn)確高效地追跡聲波在組織中的傳播路徑，是基于彎曲射線模型圖像重建方法要解決的關(guān)鍵問(wèn)題。目前基于彎曲射線模型的射線追跡方法可以分為兩大類，分別為全局法[33]和局部法[34]。全局法是基于程函方程的波前求解方法，而局部法則主要是基于射線方程的兩點(diǎn)射線追跡方法，屬于邊值問(wèn)題(Boundary value problem)。

1.2.1 基于程函方程的全局射線追跡

全局求解方法的關(guān)鍵在于程函方程的求解，主要為了確定聲波的波前位置。其中程函方程可以表示為[35]

式(2)中，T表示聲波到達(dá)時(shí)間的空間分布，v(x,y)表示聲速的空間分布。程函方程的數(shù)值解法有很多，Vidale[36?37]將有限差分方法應(yīng)用至程函方程的求解之中。Klime?[38]隨后在此基礎(chǔ)上做了改進(jìn)，將該方法推廣至各網(wǎng)格點(diǎn)處到達(dá)時(shí)間以及慢度向量的求解，其可以同時(shí)計(jì)算出網(wǎng)格中心位置的到達(dá)時(shí)間和慢度向量，網(wǎng)格內(nèi)其他位置的到達(dá)時(shí)間和慢度向量值可以利用四階拉格朗日插值方法求得。該方法后來(lái)被美國(guó)卡爾馬諾斯癌癥研究所(Karmanos Cancer Institute)運(yùn)用至超聲層析聲速重建的射線前向追跡過(guò)程中，用于重建乳腺組織的聲速分布[32,39]。

近年來(lái)，程函方程的數(shù)值求解方法不斷發(fā)展，快速行進(jìn)方法(Fast marching method FMM)以其較高的計(jì)算精度和計(jì)算速度成為求解程函方程的重要手段[40?42]。圖3顯示了點(diǎn)源發(fā)射超聲波信號(hào)經(jīng)過(guò)圓形高聲速區(qū)域時(shí)的聲波傳播軌跡。其中黑色虛線為利用FMM 方法計(jì)算出的等飛行時(shí)間分布曲線，即聲波傳播波前；黃色曲線為考慮非均勻聲速分布時(shí)聲波的實(shí)際傳播路徑；藍(lán)色直線為不考慮非均勻聲速分布時(shí)的直線傳播路徑。由圖3可以發(fā)現(xiàn)聲波經(jīng)過(guò)高聲速區(qū)域時(shí)，其波前會(huì)超前于背景聲波的波前，產(chǎn)生畸變。

圖3 基于程函方程的聲波傳播路徑追跡Fig.3 Ray tracing based on the Eikonal equation

1.2.2 基于射線方程的局部射線追跡

基于程函方程的射線追跡是一種基于全局空間到達(dá)時(shí)間計(jì)算的方法。在三維成像時(shí)，該方法的計(jì)算量會(huì)顯著增加并最終影響重建速度。利用基于射線方程的兩點(diǎn)位置彎曲射線追跡方法，可以有效減少計(jì)算量，提高計(jì)算速度。射線方程可以表示為

式(3)中，射線的起始點(diǎn)和終點(diǎn)分別為r(0)=a,r(l)=b。式(3)表明射線路徑r(s)可以由離聲源距離s、射線的起始位置a、終止位置b和介質(zhì)的折射率n確定。求解二階常微分方程式(3)的方法主要包括射線彎曲(Ray-bending)和射線鏈接(Ray-linking)兩種方法。

在射線彎曲方法中，上述求解問(wèn)題可以簡(jiǎn)化成已知兩點(diǎn)位置的邊值問(wèn)題。通過(guò)固定射線的兩個(gè)端點(diǎn)，迭代計(jì)算射線的傳播路徑，直到路徑的變化小于一定范圍[43?44]。據(jù)文獻(xiàn)[45]報(bào)道，對(duì)于聲速分布較為簡(jiǎn)單的組織，射線彎曲比射線鏈接更有效；但對(duì)于復(fù)雜的組織，射線彎曲的效率較低，有可能失敗。

射線鏈接方法常用于兩點(diǎn)之間的射線追跡問(wèn)題。當(dāng)射線的初始位置設(shè)置于發(fā)射點(diǎn)時(shí)，反復(fù)地調(diào)整射線的初始方向，直到射線的終點(diǎn)在一個(gè)很小范圍內(nèi)被接收點(diǎn)接收。將兩點(diǎn)位置的邊值問(wèn)題轉(zhuǎn)換成初值問(wèn)題進(jìn)行求解。Denis 等最早利用射線鏈接方法獲取二維聲波的彎曲射線路徑[46]，總結(jié)出射線方程的不同求解算法，并用標(biāo)準(zhǔn)數(shù)值仿體測(cè)試不同追跡算法的性能，最終得出混合步長(zhǎng)算法(Mixed-step algorithm)相比其他算法更有優(yōu)勢(shì)的結(jié)論[47]。后來(lái)，Javaherian 等[48]將射線鏈接方法應(yīng)用至半球換能器陣列和任意三維換能器陣列中聲波路徑的追跡，將發(fā)射和接收陣元之間的射線追跡問(wèn)題簡(jiǎn)化為求解最優(yōu)起始方向的問(wèn)題，并構(gòu)建了相應(yīng)的目標(biāo)方程，將問(wèn)題轉(zhuǎn)化為

式(4)中，γ(e,r)表示發(fā)射點(diǎn)e到終點(diǎn)r的方向向量；γ(e,p)[de]表示發(fā)射點(diǎn)e在初始方向de的情況下，通過(guò)射線追跡得出聲波路徑到達(dá)終點(diǎn)p的方向向量。利用優(yōu)化方法如高斯–牛頓法等求解目標(biāo)方程，最終可以獲得射線的最優(yōu)初始方向de。

雖然射線鏈接方法常用于實(shí)際生物組織中的射線追跡，但依然需要不斷地調(diào)整初始方向，使最終獲取的方向和終點(diǎn)方向重合。對(duì)此，Qu 等[49]提出了不需要反復(fù)優(yōu)化初始方向的射線鏈接方法。其基本思想為：首先定義射線的初始方向指向相應(yīng)的接收陣元，然后根據(jù)空間折射率分布以及射線追跡方程，得到下一個(gè)步長(zhǎng)?s的空間位置，直至追跡至接收陣元。該接收陣元的位置可能會(huì)偏離實(shí)際物理接收陣元，此時(shí)稱其為虛擬接收陣元。其中，射線追跡方程可以表示為

式(5)中，r和s分別表示方向向量和離發(fā)射位置的距離，n表示折射率。dr/ds可以用后向差分的形式表示，這樣可以根據(jù)初始發(fā)射位置和初始發(fā)射方向，得到虛擬接收陣元的位置，如圖4所示。最終可以根據(jù)實(shí)際接收陣元的空間位置和聲波的到達(dá)時(shí)間信息，通過(guò)鄰近陣元插值的方法估計(jì)出虛擬接收陣元的到達(dá)時(shí)間，進(jìn)而反演空間聲速分布。該方法避免了射線鏈接方法需要迭代優(yōu)化初始方向的問(wèn)題，可以有效縮短重建時(shí)間。

圖4 虛擬接收陣元和實(shí)際接收陣元位置示意圖Fig.4 Schematic diagram showing the distribution of virtual receivers and actual receivers

2 聲速重建算法

第1節(jié)討論的基于射線模型的前向過(guò)程可以獲得理論上聲波透過(guò)組織后到達(dá)換能器的時(shí)間。通過(guò)比對(duì)該理論到達(dá)時(shí)間與實(shí)驗(yàn)到達(dá)時(shí)間，構(gòu)建合適的數(shù)學(xué)模型，可以反演聲速的空間分布。

2.1 第一到達(dá)時(shí)間提取

為了重建聲速分布圖，需要首先獲取聲波透過(guò)組織到達(dá)接收換能器陣元的第一到達(dá)時(shí)間分布數(shù)據(jù)，即投影數(shù)據(jù)。由于超聲層析成像中的數(shù)據(jù)量巨大，手動(dòng)提取聲波第一到達(dá)時(shí)間的方法是不切實(shí)際的，需要可以自動(dòng)提取聲波第一到達(dá)時(shí)間的算法。聲波第一到達(dá)時(shí)間自動(dòng)提取的方法最早出現(xiàn)于地球物理學(xué)科中的地震波成像[50]，隨后被移植到醫(yī)學(xué)超聲層析成像中。目前主要可以分為3 種方法，分別是赤池信息量準(zhǔn)則(Akaike information criterion AIC)算法[51]、互相關(guān)法[52]和能量法[53]。

AIC 算法通過(guò)檢測(cè)第一到達(dá)時(shí)間前后信號(hào)的差異性，可以估計(jì)第一到達(dá)時(shí)間[51]。數(shù)學(xué)模型可以表示為

式(6)中，σ21～k和σ2k+1～N分別為聲波信號(hào)采樣點(diǎn)[1,k] 和[k+1,N]兩段信號(hào)的方差值，N表示感興趣聲波信號(hào)區(qū)域內(nèi)的采樣點(diǎn)數(shù)。AIC取值最小的點(diǎn)可以作為聲波的第一到達(dá)時(shí)間點(diǎn)。為了解決超聲層析成像中的問(wèn)題，Li 等[54]在原有AIC 算法的基礎(chǔ)上，提出了改進(jìn)的到達(dá)時(shí)間自動(dòng)提取算法。改進(jìn)主要體現(xiàn)在兩個(gè)方面：(1)用權(quán)重平均模型代替了AIC 算法中的最優(yōu)模型，在求解精度上更接近于真實(shí)情況；(2)對(duì)投影數(shù)據(jù)進(jìn)行中值濾波去除異常值，有效抑制了投影噪聲引入的圖像偽影?；ハ嚓P(guān)法利用互相關(guān)運(yùn)算得到感興趣信號(hào)和參考信號(hào)到達(dá)時(shí)間的偏移量，可以計(jì)算出感興趣信號(hào)的絕對(duì)到達(dá)時(shí)間[52]。但是，當(dāng)感興趣信號(hào)和參考信號(hào)之間存在較大波形變化時(shí)，互相關(guān)方法的有效性會(huì)面臨嚴(yán)峻挑戰(zhàn)。對(duì)此，Qu 等[55]在AIC 提取到達(dá)時(shí)間的基礎(chǔ)上，結(jié)合相鄰換能器接收聲波信號(hào)的相似性，利用互相關(guān)的方法求出相鄰陣元接收聲波信號(hào)的相對(duì)時(shí)間差，并校正AIC提取數(shù)據(jù)，可以得到較好的結(jié)果。能量法通過(guò)計(jì)算聲波信號(hào)能量，將能量最先發(fā)生變化的地方，作為聲波的第一到達(dá)時(shí)間[53]。該方法需要設(shè)定合適的閾值探測(cè)能量最開(kāi)始波動(dòng)的位置，在實(shí)際應(yīng)用中易被噪聲干擾，較少采用。

2.2 基于直線模型的濾波反投影重建

超聲斷層成像中最簡(jiǎn)單的聲速重建方法是X射線CT成像中經(jīng)典的濾波反投影(Filter back projection,FBP)算法。FBP算法基于超聲的直線前向模型假設(shè)，數(shù)學(xué)模型可以表示為

式(7)中，p(ρ,θ)為直線角度為θ時(shí)換能器陣元接收聲波信號(hào)的第一到達(dá)時(shí)間，|ξ|為濾波函數(shù)，F(xiàn) 和F?為一維傅里葉變換和反變換，s(x,y)為最終重建的空間慢度分布。根據(jù)扇形束至平行束的轉(zhuǎn)換關(guān)系，利用FBP 算法可以重建出生物組織聲速度的空間分布。Greenleaf 等[56?58]最早開(kāi)展基于FBP 的超聲層析聲速重建，可以得到初步的成像結(jié)果。

圖5為利用FBP 算法對(duì)乳腺數(shù)值仿體重建的結(jié)果?？梢钥闯觯現(xiàn)BP重建結(jié)果存在較大誤差，這主要是由于FBP 算法中的聲波直線模型的假設(shè)所致。盡管如此，F(xiàn)BP重建結(jié)果可以作為彎曲射線迭代算法和全波反演算法的初始分布值，可以加速其收斂速度。

圖5 基于直線模型的FBP 圖像重建Fig.5 Image reconstruction using the FBP algorithm

2.3 基于彎曲射線模型的迭代重建

利用聲波的直線模型假設(shè)和FBP 算法可以重建出組織的初始聲速分布。但由于實(shí)際生物組織的非均勻性，聲波在不同組織界面處會(huì)發(fā)生折射，彎曲射線模型相較于直線模型可以更好地描述該種情況下聲波的傳播軌跡。如圖6所示，基于彎曲射線模型重建結(jié)果相對(duì)于直線模型重建結(jié)果(圖5)更加準(zhǔn)確。下面將主要介紹基于彎曲射線模型的迭代重建算法。

實(shí)踐證明，“因材施教、分層培養(yǎng)、工學(xué)交替、強(qiáng)化技能”的人才培養(yǎng)模式及課程體系是高職機(jī)制專業(yè)改革的一條可行之路，較好的解決了高職學(xué)生學(xué)習(xí)積極性不足、技能優(yōu)勢(shì)不明顯的問(wèn)題。至于工學(xué)交替實(shí)訓(xùn)及課程體系中存在的一些瑕疵，有待在實(shí)踐中予以完善。

圖6 基于彎曲射線模型的重建結(jié)果Fig.6 Image reconstruction based on a bent-ray model

2.3.1 數(shù)學(xué)模型

基于彎曲射線模型的迭代重建算法原理如圖7所示。待成像區(qū)域被離散成一定大小的網(wǎng)格，環(huán)形超聲換能器陣列中的一個(gè)陣元發(fā)射超聲信號(hào)，超聲波在經(jīng)過(guò)不同聲速區(qū)域時(shí)會(huì)發(fā)生折射，最終到達(dá)接收陣元。超聲波的折射過(guò)程可以利用彎曲射線表征。

圖7 基于彎曲射線模型的迭代重建原理Fig.7 Schematic diagram showing the principle of iterative reconstruction based on bent rays

圖7中線段aij表示i條聲波射線在第j個(gè)網(wǎng)格內(nèi)的截距，即權(quán)重因子。由此，每一個(gè)發(fā)射接收對(duì)之間的聲波射線可以用一個(gè)線性方程表示，即

式(8)中，p表示聲波射線的第一到達(dá)時(shí)間，s為待重建生物組織的慢度(聲速的倒數(shù))。式(8)也可以寫為矩陣形式，即

式(9)中，p為聲波的第一到達(dá)時(shí)間向量；s為待重建組織的慢度向量；A為聯(lián)系兩者的系數(shù)矩陣，其包含了聲波射線傳播的路徑信息，和換能器的空間分布以及組織的慢度分布s直接相關(guān)。

聲速的重建過(guò)程需要引入射線軌跡校正步驟，在每次的慢度分布更新之后，需要根據(jù)更新后的聲速分布重新校正聲波的射線軌跡。由此，超聲層析成像中聲速重建的基本流程如圖8所示。

圖8 迭代重建算法流程圖Fig.8 Flow chart for iterative image reconstruction

在圖8中，為了得到穩(wěn)定的聲速分布，需要不斷比較基于彎曲射線模型的理論聲波到達(dá)時(shí)間和實(shí)驗(yàn)觀測(cè)到達(dá)時(shí)間。下面主要討論目前應(yīng)用于超聲聲速層析成像的迭代數(shù)學(xué)模型。

2.3.2 代數(shù)迭代算法

求解式(8)最簡(jiǎn)單的方法是利用X 射線CT 中代數(shù)迭代算法(Algebraic reconstruction technique,ART)對(duì)每條聲波射線經(jīng)過(guò)的網(wǎng)格慢度進(jìn)行校正[30]?？梢员硎緸?/p>

式(10)中，skj為第j個(gè)網(wǎng)格在第k次迭代時(shí)的慢度，pi為第i條射線的到達(dá)時(shí)間，aij為射線經(jīng)過(guò)網(wǎng)格的長(zhǎng)度信息，m為網(wǎng)格總數(shù)目。由于實(shí)際超聲換能器陣元數(shù)目較多，系統(tǒng)中的聲波射線數(shù)目龐大，式(10)所表示迭代過(guò)程相對(duì)較慢。Li 等[24]和Denis 等[46]采用聯(lián)合代數(shù)迭代算法(Simultaneous algebraic reconstruction technique,SART)加速傳統(tǒng)代數(shù)迭代算法，將每一個(gè)投影角度下的射線視為一組，可以提高重建速度。

2.3.3 正則化迭代算法

將基于彎曲射線模型的離散矩陣方程式(9)改寫為最小二乘形式，可以構(gòu)建出非線性目標(biāo)方程：

式(11)中，U(s)為關(guān)于慢度s的目標(biāo)方程，尋找合適的慢度分布使得最小二乘目標(biāo)方程取得最小。最常用的解法是高斯牛頓法[59]，可以不斷迭代更新慢度值s。但通常系數(shù)矩陣A為稀疏矩陣并且具有很大的條件數(shù)，使得原目標(biāo)方程具有病態(tài)性。當(dāng)噪聲增加時(shí)，將會(huì)引起重建結(jié)果的不穩(wěn)定。針對(duì)該問(wèn)題，可以通過(guò)添加正則化約束改善問(wèn)題的病態(tài)，穩(wěn)定解空間。由此，目標(biāo)方程可以改寫為

式(12)中，R(s)為正則函數(shù)；λ為正則參數(shù)，主要用于調(diào)節(jié)正則項(xiàng)的權(quán)重。Hormati 等[60]采用稀疏正則函數(shù)來(lái)表示R(s)，可以有效降低重建圖像的背景噪聲。Li 等[32]于2009年將TV 正則約束應(yīng)用于超聲聲速層析重建中，并與Tikhonov正則約束重建結(jié)果進(jìn)行了比較，得到TV 正則約束可以保持圖像邊界信息的結(jié)論。Intrator[61]和Huang等[62]針對(duì)TV正則和Tikhonov 正則各自的優(yōu)勢(shì)，將兩者進(jìn)行融合，可以兼顧Tikhonov 正則的噪聲抑制性能和TV正則的圖像邊界信息保持功能。

2.3.4 統(tǒng)計(jì)迭代算法

迭代方法中另一類重建算法是基于統(tǒng)計(jì)思想的圖像重建方法。將整個(gè)超聲換能器接收到的聲波數(shù)據(jù)看成滿足一定條件的統(tǒng)計(jì)分布，可以將重建問(wèn)題轉(zhuǎn)化為由隨機(jī)變量到達(dá)時(shí)間數(shù)據(jù)估計(jì)圖像慢度分布的問(wèn)題。其中極大似然估計(jì)(Maximum likelihood,ML)是一種運(yùn)用比較廣泛的算法，通過(guò)期望最大化方法(Expectation maximization,EM)求解極大似然模型，再通過(guò)迭代過(guò)程求解最終解。其中，迭代公式可以表示為[63]

式(13)中，為第j個(gè)網(wǎng)格第k+1次迭代的慢度值，pi為第i條聲波射線的第一到達(dá)時(shí)間，aij為射線經(jīng)過(guò)網(wǎng)格的長(zhǎng)度信息。Mer?ep 等[64]和Perez-Liva等[31]基于ML-EM 算法利用每一對(duì)發(fā)射接收陣元的投影數(shù)據(jù)校正上一次的慢度分布，可以得到更新的慢度空間分布。在噪聲比較大時(shí)，ML-EM方法直接迭代求解會(huì)出現(xiàn)不穩(wěn)定性。由此，Mer?ep 等[64]在ML-EM算法迭代過(guò)程中加入了基于中值的先驗(yàn)信息，最后形成后驗(yàn)分布即貝葉斯表達(dá)式。該思想類似于在最小二乘方程中添加正則約束。此外，Ali等[65]通過(guò)貝葉斯表達(dá)式加入了圖像的先驗(yàn)信息，最終可以得到慢度的迭代模型，進(jìn)而求解聲速的空間分布。先驗(yàn)信息的引入可以加速重建過(guò)程的收斂速度，穩(wěn)定解空間。

3 結(jié)論與展望

本文綜述了基于射線模型的超聲層析成像中聲速重建算法最新研究進(jìn)展，總結(jié)了直線和彎曲射線的前向模型，討論了聲波第一到達(dá)時(shí)間提取的不同算法。此外，綜述了圖像重建反問(wèn)題數(shù)學(xué)模型的建立和求解方法，主要包括基于直線模型的濾波反投影重建和基于彎曲射線模型的迭代重建?；趶澢渚€模型的迭代重建中主要討論了代數(shù)迭代算法、正則化迭代算法和統(tǒng)計(jì)迭代算法。為了解決數(shù)學(xué)模型中的病態(tài)問(wèn)題，可以在目標(biāo)方程中加入正則項(xiàng)或先驗(yàn)信息，以穩(wěn)定解空間，加速解的收斂速度。

通過(guò)超聲聲速重建算法的綜述，基于射線模型的聲速重建算法預(yù)計(jì)未來(lái)將著重向以下幾個(gè)方向發(fā)展：

(1)目前彎曲射線的追跡算法主要應(yīng)用于二維超聲換能器中，三維換能器中彎曲射線追跡問(wèn)題尚在研究中，開(kāi)發(fā)出能夠在三維超聲換能器中快速實(shí)現(xiàn)彎曲射線追跡是一個(gè)重要的研究方向。基于兩點(diǎn)的射線追跡算法可能可以成為三維射線追跡的重要突破口。

(2)從原始聲波信號(hào)中提取的第一到達(dá)時(shí)間是聲速重建中重要的輸入數(shù)據(jù)，它的精度直接影響最終的圖像重建質(zhì)量。盡管目前現(xiàn)存的第一到達(dá)時(shí)間提取算法具有較好的性能，但是在處理較大噪聲數(shù)據(jù)時(shí)，仍然存在提取到達(dá)時(shí)間錯(cuò)誤的情況。發(fā)展更魯棒的第一到達(dá)時(shí)間提取算法是未來(lái)發(fā)展的一個(gè)重要方向。

(3)最后是基于射線模型的聲速重建數(shù)學(xué)模型的建立問(wèn)題。由于實(shí)驗(yàn)中超聲換能器陣元數(shù)量有限，往往會(huì)出現(xiàn)待重建區(qū)域聲波射線密度不足的問(wèn)題，使得重建結(jié)果不穩(wěn)定，為此需要發(fā)展更加魯棒的數(shù)學(xué)模型。