乳腺癌檢測的三維時域微波斷層成像方法

劉廣東 張業(yè)榮

(南京郵電大學電子科學與工程學院,江蘇南京210003)

1.引 言

乳腺癌是國內(nèi)外危害女性健康的頑疾,其危害性大,患者的死亡率高。研究表明,如果乳腺腫瘤能在早期時被及時發(fā)現(xiàn)并加以治療,那么患者的五年成活率將明顯提高。早期檢測對降低乳腺癌患者的死亡率具有重要意義[1-2]。

比起乳腺癌常規(guī)檢測方法,如X線影像技術(shù)、計算機斷層攝影術(shù)、超聲成像、核磁共振成像和光學成像等方法,微波成像檢測技術(shù)具有輻射小、費用低、分辨率高和易普查等優(yōu)點,有望成為安全實用的常規(guī)或者輔助檢測手段[1-2]。

微波頻段下,正常乳房組織和惡性腫瘤組織的電特性參數(shù)差異明顯,它們的介電常數(shù)和電導率差異均在5倍以上,這為微波成像檢測乳腺癌提供了物理基礎(chǔ)[1-2]。

目前,微波成像方法主要有有源微波成像、無源微波成像和微波熱致超聲成像三種,且主要工作集中在有源微波成像。其中有源微波成像主要有微波斷層成像(microwave tomography imaging,MTI)和共焦微波成像(confocal microwave imaging,CMI)兩種方式[1-2]。

共焦微波成像方法借鑒了軍事領(lǐng)域中脈沖探地雷達技術(shù),由發(fā)射天線發(fā)射超寬帶脈沖,多個接收天線接收散射信號,基于目標的電磁參數(shù)與周圍環(huán)境的電磁參數(shù)的明顯不同,區(qū)分出強散射區(qū)域,從而判斷目標的位置[1-2]。

微波斷層成像方法是一種電磁逆散射方法,通過在散射體外部觀測到的電磁場來反演成像區(qū)域的電磁特征參數(shù)分布,從而判斷散射體目標的位置、形狀和尺寸分布等信息。所以它比共焦微波成像方法有更清晰的物理解釋[1-2]。

然而,電磁逆散射屬于不適定問題,是一個富有挑戰(zhàn)性的工作[3],非線性和病態(tài)性是其中兩個主要困難[3-11]。從處理非線性的方式來看,早期的線性化近似方法僅適用于弱散射體。對于高對比度反演問題,近年來主流方法是將重建問題轉(zhuǎn)化為最優(yōu)化問題迭代求解[3-4,10],其中,解決最優(yōu)化問題的方法主要有局部尋優(yōu)和全局尋優(yōu)兩大類[12]。無論是哪一類方法,一般都需要采用正則化方案來處理病態(tài)性[3-4,10]。

早期的電磁場數(shù)值計算多采用頻域方法,它對窄帶信號而言可行、經(jīng)濟。為了檢測早期的小尺寸腫瘤病灶,需要提高重建的分辨率,因此要求利用寬帶高頻脈沖信號作為激勵源。此時,時域方法顯示了更大的優(yōu)勢。時域電磁逆散射是利用目標的瞬態(tài)響應波形重建目標的物理特征。相比頻域方法,時域方法由于包含了更多的信息,可能獲得更準確的重建結(jié)果[4,10-11]。

得益于時域有限差分(finite-difference timedomain,FDTD)法[13]和計算機水平的快速發(fā)展,近十年來已經(jīng)成功建立了好幾種時域電磁逆散射算法,這其中主要有Takenaka等提出的正反演時間步進(forward-backward time-stepping,FBTS)算法[4]和Rekanos等提出的拉格朗日乘子算法[7]以及Winters等人提出的時域逆散射 (time-domain inverse scattering,TDIS)算法[10-11]。其中FBTS算法在二維乳腺癌檢測[4,14]、TDIS算法在二維和三維乳腺癌檢測[10-11]中均已獲得了重要進展[10-11],而拉格朗日乘子算法在這方面的應用還未見報道。

我們首先建立三維半球乳房模型;接著以Rekanos等人的拉格朗日乘子算法[7]為基礎(chǔ),應用泛函分析和變分法,導出閉式的拉格朗日(Lagrange)函數(shù)關(guān)于乳房內(nèi)各介質(zhì)電參數(shù)的Fréchet導數(shù);其次,考察了噪聲的影響,并加進了一階的吉洪諾夫(Tikhonov)正則化方案[15],改進為三維時域微波斷層成像算法;進而運用FDTD法和Polak-Ribière-Polyak(PRP)非線性共軛梯度(conjugate gradient,CG)法[12],對一個數(shù)值算例進行迭代計算,最后給出定量的相對介電常數(shù)和電導率分布圖像。

2.成像方法

建立如圖1所示的三維半球乳房模型和直角坐標系。乳房浸入無耗的匹配液體介質(zhì)中。本文假設(shè)所有介質(zhì)的電磁參數(shù)為各向同性,本構(gòu)關(guān)系為線性,相對磁導率μr為1。

圖1 三維半球乳房模型

圖1 中,T和R分別表示收發(fā)分置(即雙站)的發(fā)射和接收陣列天線[16],M根發(fā)射天線元和N根接收天線元分別位于r=r m和r=r n位置,其中m=1,2,…,M,n=1,2,…,N。這里,位置矢量r=(x,y,z)T,上標T表示轉(zhuǎn)置(下文同)。依次激勵其中的一根發(fā)射天線元,產(chǎn)生的電磁總場由N根接收天線元接收。

V 1表示待重建的乳房區(qū)域,V 2表示已知電磁參數(shù)的匹配液體和胸壁區(qū)域。整個區(qū)域的相對介電常數(shù)記為εr(r),電導率記為σ(r)。電參數(shù)分布進一步簡記為

加在第m根發(fā)射天線上的激勵源為

這里

其中:δ是Dirac函數(shù);Ι(t)為時間因子,自由空間的特征阻抗η=,其中 μ0和ε0分別為自由空間的磁導率和介電常數(shù)。產(chǎn)生相應的電磁場矢量um(p,r,t)滿足算子方程

和初始條件

其中

類似于文獻[17],包含邊界條件的偏微分算子￡定義為

這里

同時滿足約束條件式(4)和(5)。

其中與剩余差有關(guān)的項為

式中,剩余差為

歸一化因子為

在Jm(r,t)激勵下的時域測量場為

這里采用了一階 Tikhonov正則化[15],γ為待定的正則化參數(shù),▽為Hamilton算子。

引入Lagrange乘子罰函數(shù)

這里,r∈V1,t∈[T,0],這樣,利用罰函數(shù)方法將原問題轉(zhuǎn)化為無約束的優(yōu)化問題。設(shè) p*是原問題的解,且w m(p,r,t)是相應的 Lagrange乘子,由 Kuhn-Tucker定理[18]知,p*必是以下 Lagrange函數(shù)

的穩(wěn)定點,其必要條件是 Lagrange函數(shù)L的一階變分δL|p=p*=0.

由于電磁場矢量u m(p,r,t)的時間t∈[0,T]是正演步進,而乘子罰函數(shù)wm(p,r,t)的時間t∈[T,0]是反演步進,這樣不便于數(shù)值計算的統(tǒng)一實現(xiàn)。若對乘子罰函數(shù)w m(p,r,t)作時間變量替換tp=T-t,產(chǎn)生相應的函數(shù) wpm(p,r,tp)。并經(jīng)過一些計算[7],一方面,可以得到,wpm(p,r,tp)滿足以下算子方程

和初始條件

其中,tp∈[0,T],激勵源為

比較式(4)、(5)和(17)、(18),容易看出wpm(p,r,tp)和um(p,r,t)滿足相似形式的偏微分算子方程,所不同的僅有激勵源項,這樣可以方便地利用相似的數(shù)值計算方法求解。另一方面,可以得到Lagrange函數(shù)L 關(guān)于p的Fréchet導數(shù)為

式中:

這樣,有了閉式的Fréchet導數(shù),本文的反演問題可以利用基于梯度的優(yōu)化方法進行迭代處理。本文采用PRP非線性共軛梯度法[12]。設(shè)k表示迭代次數(shù),這里,(k=0,1,…,kmax),kmax表示總共迭代計算次數(shù),第k次迭代的估計特征參數(shù)分布為pk,則第k+1次迭代的更新公式為

這里 ,gk、dk、βkPRP和αk分別為第k 次迭代的 Fréchet導數(shù)、搜索方向、標量因子和步長,其中步長αk可以通過求解線搜索問題式(26)得到[7]。

測量的電磁場u^m(r n,t)可以通過實驗得到或者利用數(shù)值方法仿真代替,并且設(shè)定重建參數(shù)的迭代初值pk(k=0)以后,可以利用上述算法迭代求解重建參數(shù)pk,直到達到所需要的迭代次數(shù)或者計算精度為止。總結(jié)起來,三維時域微波斷層成像算法的計算步驟如下:

①給定參數(shù)T、M、N、總共迭代次數(shù)k max或者閾值eps＞0以及入射脈沖 I(t),t∈[0,T],實驗測量或者仿真計算場u^m(r n,t),其中t∈[0,T],m=1…M,n=1…N。

②初始化:取k=0,給定p0,仿真計算接收天線上的場u m(p0,r n,t)和成像區(qū)域的場u m(p0,r,t),其中,t∈[0,T],r∈V 1,m=1,…,M,n=1,…,N 。對(rn,t)和um(p0,rn,t)作時間變量替換t=T-tp,計算激勵源Jpm(r n,tp),仿真計算成像區(qū)域函數(shù)wpm(p0,r,tp),作時間變量替換tp=T-t,得到成像區(qū)域的乘子罰函數(shù) wm(p0,r,t),r∈V1。計算Fréchet導數(shù) g0,計算搜索方向 d0,計算步長 α0。

③令k=k+1,更新估計參數(shù)pk,計算接收天線上的場um(pk,rn,t)和成像區(qū)域的場um(pk,r,t),其中,t∈[0,T],r∈V 1,m=1,…,M,n=1,…,N。

⑤計算Fréchet導數(shù) gk,如果‖gk‖ ≤eps或者k=則停;否則進入下一步。

⑥計算標量因子βkPRP和搜索方向dk.

⑦計算步長αk,轉(zhuǎn)③步。

3.仿真算例

為了驗證本文成像算法的性能,數(shù)值算例中,取乳房半球的半徑為50 mm,外表2 mm厚為皮膚層,其余為脂類乳房組織。胸壁厚為30 mm,未被乳房覆蓋的表層也接2 mm厚的皮膚層。設(shè)乳房組織中含有兩個半徑均為3 mm的球狀腫瘤1和2。腫瘤1的中心坐標分別為x1=45 mm、y 1=70 mm和z1=45mm;腫瘤2的中心坐標分別為 x2=70 mm、y2=70 mm和z2=45 mm。

發(fā)射陣列天線取最簡單的情形,只有1根天線元,即M=1,貼在乳房頂部放置。接收天線陣列按照乳房高度均分5層,由上向下分別在乳房表面圓周上均勻放置 4、5、6、7和 8根天線元,即 N=30。這里,我們忽略對天線元具體形狀的建模,即分別視為發(fā)射和接收點。

激勵源式(3)中的時間因子I(t)取為[19]

式中:τ=100 ps;t0=4·τ;中心頻率 f=3.2 GHz。

在頻率 f下,用于仿真測量數(shù)據(jù)的各種介質(zhì)的電特性參數(shù)如表1所示[16]。

表1 介質(zhì)的電特性參數(shù)

定義第k次迭代時,重建過程的相對剩余誤差為[3]RREk

利用FDTD[13]法仿真計算代替測量電磁場u^m,每次迭代過程中,也用FDTD法計算u m和w m,取總共迭代次數(shù)kmax=50。計算空間離散化時,Yee元胞尺寸取為Δ=1 mm×1 mm×1 mm,時間采樣步長為 Δt=Δ/(2·c),整個計算空間在 x、y和z方向剖分的元胞數(shù)目分別為140、140和90,周圍用8層完全匹配層(perfectly matched layer,PML)吸收邊界截斷[13]。觀測時間取為 T=500·Δt。

成像的目標是重建乳房區(qū)域的電參數(shù)分布,k=0時,迭代計算的初值p0取等于均勻乳房組織的相應特征參數(shù)。

在實際的測量中一般存在噪聲污染,本文假定,在仿真的測量數(shù)據(jù)中加入均勻分布的隨機噪聲,其信噪比定義為[17]

另外,在真實的乳房組織中還存在乳腺導管和小葉等組織,它們降低了腫瘤的檢測效果,為了考察它們對反演算法的影響,取乳房組織的電特性參數(shù)在正常值的1±10%倍范圍內(nèi)均勻隨機變化[14]。

為了對抗噪聲污染和降低反演的病態(tài)性質(zhì),避免求解結(jié)果限于局部最優(yōu),這里加入一階的吉洪諾夫正則化,正則化參數(shù)取γ=0.001[7]。

通過腫瘤1和2中心的x-y橫截面、通過腫瘤1和2中心的 x-z橫截面和通過腫瘤2中心的y-z橫截面上得到的εr分布分別如圖2、3和 4所示;通過腫瘤 1和 2中心的x-y橫截面、通過腫瘤1和2中心的 xz橫截面和通過腫瘤2中心的 y-z橫截面上得到的σ分布分別如圖5、6(見 1234頁)和7所示 。其中各子圖(a)、(b)、(c)和(d)分別表示相應電參數(shù)的真實分布以及在迭代次數(shù)分別為k=0、25、50時的重建結(jié)果。為了檢驗成像方法檢測腫瘤的效果,各子圖(d)中附加的黑色圓周表示真實腫瘤的輪廓。

圖8給出相對剩余誤差隨迭代次數(shù)的變化關(guān)系。其中,第25和50次迭代的相對誤差分別為4.6×10-2%和4.2×10-3%。

圖8 不同迭代次數(shù)的相對剩余誤差

比較圖2～8,經(jīng)過分析發(fā)現(xiàn):1)從總體來看,算法是收斂的,收斂速度先快后慢。隨著迭代次數(shù)的增加,相對誤差逐漸減小,重建精度逐漸提高。當然,迭代次數(shù)越高,計算時間也越長,實際的工程應用中,可以選擇一個重建精度與實時重建的合理折中。2)皮膚、腫瘤和乳房組織的電特性參數(shù)差異較大,是主要的散射體。迭代次數(shù)k=50時,σ分布收到了較好的重建效果,較好地重現(xiàn)了它們的形狀、尺寸、位置。然而,εr分布和真實值尚有較大的差異,適當增加仿真時間步數(shù)或者適當加大迭代次數(shù)可以適當改善重建結(jié)果。3)對比皮膚和腫瘤,前者比后者的重建效果好;對比腫瘤1和腫瘤2,后者的重建效果不如前者。其原因是腫瘤比起皮膚位于較深位置,而且腫瘤2位于更深位置,距離天線較遠,信號在乳房組織中經(jīng)受衰減,來自腫瘤2的散射信號較弱所致,解決方法是適當降低激勵源的頻率。當然這也同時降低了重建的分辨率,在實際的工程應用中,應該選擇一個檢測深度和重建分辨率的合理折中,即選擇一個合理的激勵頻率。

綜合起來,盡管由于測量視角的受限和測量數(shù)據(jù)的有限增加了反演問題的病態(tài)性質(zhì),但是在噪聲環(huán)境下,在電特性參數(shù)取正常值的1±10%倍均勻隨機變化的乳房組織中,腫瘤目標仍被成功地探測。這可能得益于算法的魯棒性和正則化方案改善了逆問題的病態(tài)性質(zhì)。當然,正則化方案已有多種,正則化參數(shù)的選擇也是個困難的問題,最優(yōu)的正則化參數(shù)一般需要根據(jù)具體問題多次數(shù)值測試確定。

4.結(jié) 論

為了提高成像的分辨率,常利用高頻寬帶脈沖進行生物異常體檢測,在時域進行微波斷層成像是一種重要的生物醫(yī)學成像方法。用它進行早期的乳腺癌檢測具有良好的應用前景和重要的現(xiàn)實意義。

本文應用泛函分析和變分法,改進了拉格朗日乘子算法并得到閉式的Lagrange函數(shù)關(guān)于重建參數(shù)的Fréchet導數(shù),采用了一階的 Tikhonov正則化方案。在噪聲環(huán)境下,利用FDTD法和PRP非線性CG法對三維半球乳房模型進行了仿真測試,介電常數(shù)和電導率分布圖像均被定量地獲得,散射體的位置、尺寸和形狀等目標信息得以成功地再現(xiàn)。電導率和介電常數(shù)雙參數(shù)反演方法的研究,豐富了寬帶高頻電磁脈沖在有耗介質(zhì)里的傳播理論,使得電磁勘探方法的反演解釋方法得到一般化,具有一定的理論意義。初步的結(jié)果證實了改進后的時域微波斷層成像方法檢測乳腺癌的可行性和有效性,表明了正則化方案對反演問題的病態(tài)性質(zhì)起到了一定的抑制作用,為工程應用研究奠定了良好的理論基礎(chǔ)。

當然,本文的時域微波斷層成像算法并不局限于乳腺癌檢測,也可應用到目標識別、地球物理勘探和遙感等一般的電磁逆散射問題中。另外,生物肌體組織一般是頻率色散介質(zhì),特別是對超寬帶微波檢測,色散效應有重要影響[20],相關(guān)問題將另文討論。

[1] 劉廣東,張業(yè)榮.An overview of active microwave imaging for early breast cancer detection[J].南京郵電大學學報(自然科學版),2010,30(1):64-70.

LIU Guangdong,ZHANG Yerong.An overview of active microwaveimaging for early breast cancer detection[J].Journal of Nanjing University of Posts and Telecommunications(Natural Science),2010,30(1):64-70.(in Chinese)

[2] 田雨波,錢 鑒.微波近場成像檢測乳腺癌及其微波熱療[J].微波學報,2003,19(3):72-78.

TIAN Yubo,QIAN Jian.Breast cancer detection based on microwave near field imaging and its microwavehyperthermia[J].Journal of Microwaves,2003,19(3):72-78.(in Chinese)

[3] 張業(yè)榮,聶在平,阮穎錚.用變形玻恩迭代法反演電導率的二維非均勻分布[J].電子學報,1997,25(12):100-111.

ZHANG Yerong,NIE Zaiping,RUAN Yingzheng.Inversion of two-dimensional inhomogeneous conductivity distribution using the distorted born iterative method[J].Acta Electronica Sinica,1997,25(12):100-111.(in Chinese)

[4] TAKENAKA T,JIA H,TANAKA T.Microwave imaging of electrical property distributions by a forward-backward time-stepping method[J].J.Electromagn.Waves Appl.,2000,14(12):1609-1626.

[5] 徐偉強,王敏錫.用求解二維導體柱電磁逆散射[J].電波科學學報,2002,17(5):462-466.

XU Weiqiang,WANG Minxi.Electromagnetic inverse scattering of two-dimensional conducting objects by adaptive parameter real-coded genetic algorithm[J].Chinese Journal of Radio Science,2002,17(5):462-466.(in Chinese)

[6] 王文珂,黃卡瑪,劉智明,等.遺傳算法在生物組織電導率重構(gòu)中的應用研究[J].電波科學學報,2003,18(4):400-403.

WANG Wenke,HUANG Kama,LIU Zhiming,et al.Study on conductivity reconstruction of biological tissue by using genetic algorithm[J].Chinese Journal of Radio Science,2003,18(4):400-403.(in Chinese)

[7] REKANOS I T.Time-domain inverse scattering using Lagrange multipliers:an iterative FDTD-based optimization technique[J].J.Electromagn.Waves Applicat.,2003,17(2):271-289.

[8] 鄧小波,聶在平,趙延文,等.利用相位感應測井數(shù)據(jù)重建電導率和介電常數(shù)[J].電波科學學報,2005,20(2):236-240.

DENG Xiaobo,NIE Zaiping,ZHAO Yanwen,et al.Reconstruction of permittivity and conductivity using phase induction logging information[J].Chinese Journal of Radio Science,2005,20(2):236-240.(in Chinese)

[9] 傅 林,黃卡瑪,向勝昭.生物組織磁聚焦電導率成像原理及反演算法[J].電波科學學報,2006,21(2):249-254.

FU Lin,HUANG Kama,XIANG Shengzhao.Focused magnetic fields conductivity tomography of biological tissue and its inverse algorithm[J].Chinese Journal of Radio Science,2006,21(2):249-254.(in Chinese)

[10]WINTERS D W,BOND E J,BARRY D V,et al.Estimation of the frequency-dependent average dielectric properties of breast tissueusing a time-domain inverse scattering technique[J].IEEE Trans.Antennas Propag.,2006,54(11):3517-3528.

[11] WINTERS D W,SHEA J D,KOSMAS P,et al.Three-dimensional microwave breast imaging:dispersive dielectric properties estimation using patientzaspecific basis functions[J].IEEE Trans.Med.Imag.,2009,28(7):969-981.

[12] DAI Yuhong.A family of hybrid conjugate gradient methods for unconstrained optimization[J].Math.of Computation,2003,17(2):1317-1328.

[13] TAFLOVE A,HAGNESSS C.Computational Elec-trodynamics:The Finite-Difference Time-Domain Method[M].3rd ed.Norwood,MA:Artech House,2005.

[14] JOHNSON J E,TAKENAKA T,TANAKA T.Two-dimensional time-domain inverse scattering for quantitative analysis of breast composition[J].IEEE Trans.Biomed.Eng.,2008,55(8):1941-1945.

[15] TIKHONOV A N,ARSENIN V T.Solutions of Ill-Posed Problems[M].Washington,D.C.:Winston,1977.

[16] XIE Yao,GUO Bin,XU Luzhou,et al.Multistatic adaptive microwave imaging for early breastcancer detection[J].IEEE Trans.Biomed.Eng.,2006,53(8):1647-1657.

[17] TAKENAKA T,ZHOU H,TANAKA T.Inverse scattering for a three-dimensional object in the time domain[J].J.Opt.Soc.Am.A,2003,20(10):1867-1874.

[18] KUHN H W,TUCKER A W.Nonlinear Programming[M].Berkeley,California,University of California Press,1951:481-492.

[19] ZHOU Beibei,SHAO Wenyi,WANG Gang.The application of multi-look in UWB microwave imaging for early breast cancer detection using hemispherical breast model[C]∥IEEE EMBS,2005:1552-1555.

[20] 劉廣東,張業(yè)榮.二維有耗色散介質(zhì)的時域逆散射方法[J].物理學報,2010,59(10):6966-6976.

LIU Guangdong,ZHANG Yerong.Time-domain inverse scattering problem for two-dimensional frequency-dispersive lossy media[J].Acta Phys.Sin.,2010,59(10):6966-6976.(in Chinese)