基于射線模型的超聲聲速層析算法研究進展?

朱 翔 田 超

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

0 引言

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

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

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

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

1 射線前向模型

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

1.1 直線模型

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

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

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

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

1.2 彎曲射線模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 聲速重建算法

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

2.1 第一到達時間提取

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

AIC 算法通過檢測第一到達時間前后信號的差異性，可以估計第一到達時間[51]。數(shù)學(xué)模型可以表示為

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

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

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

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

圖5為利用FBP 算法對乳腺數(shù)值仿體重建的結(jié)果?？梢钥闯?，F(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 算法可以重建出組織的初始聲速分布。但由于實際生物組織的非均勻性，聲波在不同組織界面處會發(fā)生折射，彎曲射線模型相較于直線模型可以更好地描述該種情況下聲波的傳播軌跡。如圖6所示，基于彎曲射線模型重建結(jié)果相對于直線模型重建結(jié)果(圖5)更加準確。下面將主要介紹基于彎曲射線模型的迭代重建算法。

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

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

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

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

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

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

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

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

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

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

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

2.3.2 代數(shù)迭代算法

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

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

2.3.3 正則化迭代算法

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

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

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

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

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

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

3 結(jié)論與展望

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

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

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

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

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