那雪璐 李海洋? 張文棟

(1 中北大學信息與探測處理山西省重點實驗室太原 030051）

(2 中北大學動態(tài)測試省部共建實驗室太原 030051）

0 引言

乳腺癌是我國發(fā)病率最高的女性惡性腫瘤。早期乳腺癌的篩查主要依賴于定期體檢，檢查手段有乳腺鉬靶、超聲檢查、X線檢查等。其中，超聲檢查以其價格低廉、無輻射等優(yōu)點，已經成為乳腺癌診斷的重要工具。

乳腺是一個從外部到內部由皮膚、脂肪層、腺體等組織層層分布或相互交錯組成的復雜的生物結構體。不同乳腺組織的平均聲速不同[1]，而病變組織一般具有更大的聲速，因此本文采用聲速作為反演參量可以精確判斷病變組織的位置及大小。相較于二維成像，三維成像的優(yōu)勢在于能夠從各個角度全方面地呈現(xiàn)乳腺組織內部構造，但三維成像的算法復雜度較高，故目前乳腺超聲成像的研究成果更多集中在二維成像，對于乳腺結構的三維成像研究開展較少。

圍繞乳腺組織的二維超聲成像也可稱為超聲層析成像，目前超聲成像算法主要有基于射線理論的反演和全波形反演兩大類。理論上，基于波動方程的全波形反演[2-4]考慮了運動學和動力學特征，能夠充分利用信號所攜帶的振幅、相位等信息，成像分辨率明顯高于基于射線理論的成像方法。然而，在實際應用中，全波形反演成像仍存在一些困難，例如：對初始模型依賴程度較高、反演時目標函數(shù)與速度擾動之間強烈的非線性關系、計算量較大等，以上因素使得全波形反演未能得到廣泛應用。射線反演成像[5-8]基于高頻近似理論，將聲波傳播路徑近似為射線，使得反演數(shù)學模型簡單、計算量小，在滿足聲波波長遠小于待反演目標尺寸的前提下具有較好的超聲成像質量。實際聲波傳播過程中，影響聲波傳播的不僅僅局限于射線上，而在射線鄰近區(qū)域內的生物組織結構也會對聲波傳播造成影響，這個區(qū)域被稱為菲涅爾區(qū)域[9-12]。與傳統(tǒng)射線相比，采用菲涅爾區(qū)域描述聲波傳播更加符合實際聲波傳播過程。Fang等[13]將菲涅爾區(qū)域應用于乳腺成像，并分析了成像誤差，證明了菲涅爾區(qū)域成像算法在乳腺成像方面的可行性。

本文主要內容包括分析了環(huán)形傳感陣列菲涅爾區(qū)域成像法基本原理，設計了基于柱面?zhèn)鞲嘘嚵械娜S成像算法；通過計算三維超聲反演誤差，實現(xiàn)了對乳腺內部聲速反演精確度的評價；最后，本文通過對菲涅爾區(qū)域法成像分辨率進行分析，研究了不同入射頻率對乳腺內部聲速反演的影響。

1 理論與方法

1.1 環(huán)形傳感陣列與菲涅爾區(qū)域法

相比于射線理論，菲涅爾區(qū)域法是將聲波傳播看作是一個體積范圍內所有點共同作用的結果，而這個體積是傳播的直達波與介質中散射波相互干涉的區(qū)域。如圖1(a)所示，s點是發(fā)射傳感器位置，g點是接收傳感器位置，r點是空間介質中任意一點。

圖1 菲涅爾區(qū)域Fig.1 Fresnel region

圖1中黑色實線為射線理論傳播路徑，紅色實線是聲波在不均勻介質中的實際傳播路徑。定義Δt為擾動時間，表示聲波在實際路徑(s點→r點→g點)的傳播時間(tsr+trg)與在射線理論路徑(s點→g點)的傳播時間tsg之間的差值，即Δt=tsr+trg-tsg。利用最大擾動時間Δtmax定義菲涅爾區(qū)域分布范圍，計算得到[9]

其中，T是入射聲波的周期。當散射點滿足擾動時間小于(3×T)/8時，該點就在菲涅爾區(qū)域范圍之內，表示該點是導致接收點g點接收到的聲波時間發(fā)生擾動的其中一個因素。

以圖1中s點與g點為例，兩點之間透射時間最小的路徑是直線路徑，此時Δt等于0，而菲涅爾區(qū)域內其他r點產生的Δt是小于(3×T)/8的一個數(shù)值。因此，可采用單調線性遞減加權函數(shù)來表征菲涅爾區(qū)域內的點對第一到達時間的影響，加權函數(shù)定義[13]如下：

其中，f是入射聲波頻率。在均勻介質中，由于均勻介質內部聲速分布恒定，具有相同擾動時間的散射點呈對稱分布，所以均勻介質中發(fā)射點與接收點之間的菲涅爾區(qū)域呈橢圓形狀。當介質中出現(xiàn)不均勻聲速分布，菲涅爾區(qū)域相應地產生彎折現(xiàn)象，如圖2所示。

圖2 非均勻介質造成的菲涅爾區(qū)域彎折Fig.2 Curved Fresnel region caused by the inhomogeneous medium

如圖3，環(huán)形傳感陣列以逐一發(fā)射、多點接收的掃查模式，采用透射法實現(xiàn)對生物組織內部慢度的反演成像。檢測過程中，為覆蓋檢測區(qū)域，選擇與發(fā)射傳感器位置相對的傳感器作為接收傳感器，對每一對發(fā)射與接收傳感器分別采集透射聲信號并提取透射時間，再由式(2)根據(jù)介質擾動時間計算菲涅爾區(qū)域內任一點的加權值，進一步利用式(3)迭代更新反演介質內擾動慢度分布，其基本原理是利用通過任一網格點的所有菲涅爾區(qū)域加權值和計算的透射時間相對誤差計算每次迭代的聲速慢度變化值，以此不斷修正聲速慢度。

圖3 環(huán)形傳感陣列Fig.3 Ring sensor array

其中，i、j和k分別表示收發(fā)探頭對、網格點和迭代的后綴，是由如式(4)表示的時間殘差：

其中，Tobs和Tcal表示相對于當前模型觀測的透射時間和由超聲數(shù)據(jù)提取出的透射時間,本文分別由有限差分法[14]和赤池信息準則法[15]計算得到。

隨著迭代次數(shù)不斷增加，當ΔT開始增大時停止更新慢度。由此，采用菲涅爾區(qū)域法，建立基于環(huán)形傳感陣列的乳腺二維層析成像方法。

1.2 柱面?zhèn)鞲嘘嚵信c乳腺三維成像

本文采用柱面?zhèn)鞲嘘嚵凶鳛槌暀z測信號激勵與接收裝置。如圖4所示，該裝置可看作由N個x-y平面內環(huán)形傳感陣列沿著z軸等間隔分布構成，且各環(huán)形傳感器陣元在z軸呈直線分布(如圖3紅色虛線所示)。檢測過程中，將柱面?zhèn)鞲嘘嚵锌傮w浸在水中并放置在檢測平臺上，患者俯趴在檢測平臺使乳房下垂并完全放置在柱面?zhèn)鞲嘘嚵兄小Ｔ趚-y-z坐標空間中，乳腺是一個沿著z軸分布的類圓錐體，沿水平方向可將其離散成數(shù)片在x-y平面內的冠狀面，結合柱面?zhèn)鞲嘘嚵械奶攸c，即可將乳腺三維成像問題轉為多個對環(huán)形傳感陣列的二維層析成像問題?？梢?，柱面?zhèn)鞲嘘嚵锌赏ㄟ^控制N個環(huán)形傳感陣列依次發(fā)射與接收得到不同位置的乳腺冠狀面圖像，然后將所有的二維乳腺冠狀面進行數(shù)據(jù)整合，實現(xiàn)乳腺三維超聲成像。

圖4 柱面?zhèn)鞲嘘嚵惺疽鈭DFig.4 Schematic diagram of cylindrical sensor array

根據(jù)上述對于聲速慢度成像的二維及三維成像理論分析，利用柱面?zhèn)鞲嘘嚵械膬?yōu)勢，將三維問題降維，分解為多個二維成像問題，能夠從一定程度上減弱三維成像相較二維成像的算法計算量大的問題。實現(xiàn)乳腺三維聲速成像的具體流程如圖5所示。

圖5中，首先根據(jù)柱面陣中環(huán)形傳感陣列的分布將乳腺劃分為數(shù)個冠狀面；然后基于1.1節(jié)菲涅爾區(qū)域法實現(xiàn)乳腺冠狀面的二維慢度成像；最后，結合柱面?zhèn)鞲嘘嚵械姆植继卣髡蠑?shù)據(jù)呈現(xiàn)乳腺組織三維成像。

圖5 三維成像流程圖Fig.5 The flow chart of three-dimensional imaging

2 仿真計算

本文采用聲學工具箱k-Wave實現(xiàn)乳腺中聲波時域仿真與信號采集。按照1.2節(jié)要求，設計如圖6所示的三維模型。定義乳腺病變標號如圖6(b)所示，與表1參數(shù)相對應，且全文病變標號一致。

圖6 乳腺三維仿真模型Fig.6 The three-dimensional breast model

表1 三維乳腺模型仿真參數(shù)Table 1 Simulation parameters of three-dimensional breast model

設置如表2所示的仿真參數(shù)：入射信號中心頻率為1 MHz，網格步長為0.4 mm，網格步長小于1/3波長，滿足了空間采樣定理。時間步長設為0.075 μs，滿足了時間采樣定理：dt=(CFL×dx)/cmax，其中，CFL(Courant-Friedrichs-Lewy)數(shù)設為0.3，dt表示時間步長，dx表示網格步長，cmax表示乳腺組織中最大的聲速。本文依次將每個陣元作為發(fā)射端，位于發(fā)射陣元對立面的覆蓋約270°的陣元作為接收信號，實現(xiàn)一發(fā)多收，滿足環(huán)形傳感陣列的成像需求。組成柱面環(huán)形傳感陣列的一個環(huán)形陣列如圖7(a)所示，聲波傳播如圖7(b)所示，發(fā)射聲信號波形如圖7(c)所示。

圖7 基于k-Wave軟件的環(huán)形傳感陣列及時域波形Fig.7 Ring sensor array and time domain waveform based on k-Wave toolbox

表2 仿真參數(shù)Table 2 Simulation parameters

綜上，本文利用k-Wave工具箱建立了乳腺三維模型，并且可獲得冠狀面、橫切面和矢狀面的乳腺切片圖像。采用一發(fā)多收的方式，基于柱面?zhèn)鞲嘘嚵袑崿F(xiàn)了聲波時域仿真與透射信號采集，為乳腺三維成像提供了原始數(shù)據(jù)。

3 結果與討論

3.1 三維反演質量分析

距離乳腺上表面19.6 mm、24 mm處的冠狀面反演圖像如圖8(a)、圖8(b)所示。基于柱面三維傳感陣列，每隔0.4 mm獲得二維層析成像整合成三維成像，獲取病變1中心位置橫切面及矢狀面反演圖像，如圖8(c)所示。

由圖8中的乳腺反演圖像可清晰辨認出病變，證明菲涅爾區(qū)域法可實現(xiàn)乳腺組織三維反演，且基于柱面?zhèn)鞲嘘嚵械娜橄偃S成像可以獲取任意位置的空間信息，方便了乳腺早期病變診斷。為了進一步說明反演質量，獲取如圖8(a)、圖8(b)中的反演一維聲速剖面與實際聲速對比，如圖9所示。

圖8 乳腺三維反演圖像Fig.8 Three-dimensional inversion image of breast

由圖9看出，本文采用菲涅爾區(qū)域法可以實現(xiàn)乳腺組織聲速的精確測量。為更好地分析菲涅爾區(qū)域法實現(xiàn)乳腺組織反演成像質量，需要對反演聲速圖像進行定量評估。本文分別對病變成像聲速、尺寸誤差進行計算，定義病變聲速測量絕對誤差為

圖9 乳腺反演成像中的一維聲速剖面Fig.9 One-dimensional sound velocity profile in breast inversion imaging

其中，ssext表示病變反演聲速最值，ssmod表示模型中病變的聲速。

為實現(xiàn)病變尺寸的定量測量，結合病變聲速成像特點，提出公式(6)作為病變尺寸測量線：

其中，ssgland表示背景腺體聲速。

利用式(5)、式(6)分別計算各病變聲速絕對誤差和尺寸反演誤差以評估反演精度，誤差如表3、表4所示。

表3 聲速誤差Table 3 Sound velocity error

表4 尺寸誤差Table 4 Dimension error

觀察病變1、病變2、病變3可見：相對于背景腺體的聲速差越大其聲速反演誤差越大；觀察病變4、病變5、病變6可見：病變尺寸越小，聲速反演誤差越大。這是由于穿過小尺寸病變的菲涅爾體積更少，迭代時的加權樣本更少，增加了小尺寸病變的聲速反演誤差。

由式(6)可知，反演聲速最值與背景腺體聲速是影響尺寸估計的主要參數(shù)，而病變反演聲速最值是一個空間網格、病變實際聲速、病變實際尺寸、迭代截止條件等仿真參數(shù)綜合作用的結果，這加劇了本文基于仿真方法研究乳腺病變尺寸測量分析的難度。表4中病變1、病變2、病變3、病變4的成像尺寸誤差(1～1.4 mm之間)相近，而病變6的尺寸誤差最大，這與表3中的成像聲速誤差相對應：病變1、病變2、病變3、病變4的反演聲速誤差(0～4 m·s-1之間)相近，而病變6具有最大的聲速誤差，這表明病變尺寸成像誤差受到反演聲速質量的影響較大。而病變5的尺寸測量誤差最小則可能是由于病變5的尺寸設置處于中間值，本文的尺寸測量標準對于該尺寸的病變更有利。

3.2 成像分辨率

病變成像分辨率的定義有兩種：一是能夠成像的最小病變尺寸，二是兩病變間的最小分辨率。本文為研究菲涅爾區(qū)域法成像的最小分辨率，分別建立了單病變模型和雙病變模型，如圖10(a)、圖10(b)所示。其中，圖10(b)中兩個病變完全相同，參數(shù)設置如表5所示。

圖10 乳腺模型Fig.10 The breast model

表5 乳腺模型參數(shù)Table 5 The breast model parameters

對于單病變模型，不斷縮小病變尺寸(4 mm→2 mm→1 mm→0.6 mm→0.4 mm)進行成像，如圖11所示。

由圖11所示，該成像條件下菲涅爾區(qū)域法對單病變的最小分辨率為2 mm，2 mm以下病變成像無法明顯觀察到。

圖11 不同尺寸單病變模型成像Fig.11 Imaging of different lesion sizes in single lesion model

對于雙病變模型，不斷縮小兩病變間距

(6 mm→5 mm→4 mm→3 mm→2 mm→1 mm)進行成像。如圖12所示，當病變間距縮小至3 mm時，成像圖像中兩病變不能完全分離，證明該成像條件下菲涅爾區(qū)域法對雙病變的最小分辨率為4 mm。

圖12 不同間距雙病變模型成像Fig.12 Imaging with different lesion spacing in dual-lesion model

可見，相較單病變模型，雙病變模型對于病變間距成像分辨率更低。為進一步提高雙病變成像分辨率，分析各參數(shù)對成像分辨率的影響。顯而易見，網格步長的縮小會相應地提高成像分辨率，但隨著網格步長的縮小，求解聲速矩陣變大，運算復雜度增加，求解速度變慢；陣元個數(shù)增多也會相應地提高成像分辨率，但由于接收數(shù)據(jù)的增多，增加了成像數(shù)據(jù)量。由此，本文僅研究入射頻率對雙病變成像分辨率的影響，網格步長設為0.2 mm，環(huán)陣陣元個數(shù)設為128。在保證采樣定理滿足的情況下，逐步提高入射頻率(1 MHz→1.5 MHz→2 MHz→2.5 MHz→3 MHz→3.5 MHz)，對3 mm病變間距成像結果如圖13所示。

圖13 不同頻率下雙病變成像圖像Fig.13 Dual-lesion images at different frequencies

由于入射頻率的提高縮小了菲涅爾區(qū)域的大小，1 MHz、2 MHz、3 MHz時的菲涅爾成像區(qū)域分別為3.7500×10-7s、1.8750×10-7s、1.2500×10-7s，使得雙病變成像圖像更加精細，分辨率提高。但同時，隨著入射頻率的增加，依照式(1)計算的菲涅爾區(qū)域逐漸減小，則需要更多的陣元數(shù)以完全覆蓋成像區(qū)域，才能避免如圖13(d)、圖13(e)、圖13(f)中病變成像處出現(xiàn)區(qū)域分塊?？梢?，適當增加入射頻率對成像分辨率的提高有很大改善。

4 結論

本文依據(jù)柱面陣三維超聲成像原理，將乳腺三維超聲成像問題降為二維乳腺冠狀面成像問題，利用菲涅爾區(qū)域成像算法對乳腺二維冠狀面反演成像，進一步實現(xiàn)了由二維到三維的乳腺超聲成像。通過成像誤差計算，證明菲涅爾區(qū)域成像算法可以應用于乳腺三維超聲成像問題，能夠高質量地實現(xiàn)對乳腺三維結構的反演。同時，通過對成像分辨率分析，提出了實現(xiàn)高質量乳腺成像需要較高頻率的探頭，為實際中乳腺三維超聲成像問題的解決提供了一種可行的技術方案。