劉廣東 張業(yè)榮

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

1.引 言

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

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

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

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

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

微波斷層成像方法是一種電磁逆散射方法,通過在散射體外部觀測到的電磁場來反演成像區(qū)域的電磁特征參數(shù)分布,從而判斷散射體目標(biāo)的位置、形狀和尺寸分布等信息。所以它比共焦微波成像方法有更清晰的物理解釋[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ì)算多采用頻域方法,它對窄帶信號而言可行、經(jīng)濟(jì)。為了檢測早期的小尺寸腫瘤病灶,需要提高重建的分辨率,因此要求利用寬帶高頻脈沖信號作為激勵源。此時,時域方法顯示了更大的優(yōu)勢。時域電磁逆散射是利用目標(biāo)的瞬態(tài)響應(yīng)波形重建目標(biāo)的物理特征。相比頻域方法,時域方法由于包含了更多的信息,可能獲得更準(zhǔn)確的重建結(jié)果[4,10-11]。

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

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

2.成像方法

建立如圖1所示的三維半球乳房模型和直角坐標(biāo)系。乳房浸入無耗的匹配液體介質(zhì)中。本文假設(shè)所有介質(zhì)的電磁參數(shù)為各向同性,本構(gòu)關(guān)系為線性,相對磁導(dǎo)率μ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,上標(biāo)T表示轉(zhuǎn)置(下文同)。依次激勵其中的一根發(fā)射天線元,產(chǎn)生的電磁總場由N根接收天線元接收。

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

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

這里

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

和初始條件

其中

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

這里

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

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

式中,剩余差為

歸一化因子為

在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)是相應(yīng)的 Lagrange乘子,由 Kuhn-Tucker定理[18]知,p*必是以下 Lagrange函數(shù)

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

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

和初始條件

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

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

式中:

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

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

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

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

②初始化:取k=0,給定p0,仿真計(jì)算接收天線上的場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,計(jì)算激勵源Jpm(r n,tp),仿真計(jì)算成像區(qū)域函數(shù)wpm(p0,r,tp),作時間變量替換tp=T-t,得到成像區(qū)域的乘子罰函數(shù) wm(p0,r,t),r∈V1。計(jì)算Fréchet導(dǎo)數(shù) g0,計(jì)算搜索方向 d0,計(jì)算步長 α0。

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

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

⑥計(jì)算標(biāo)量因子βkPRP和搜索方向dk.

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

3.仿真算例

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

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

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

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

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

另外,在真實(shí)的乳房組織中還存在乳腺導(dǎo)管和小葉等組織,它們降低了腫瘤的檢測效果,為了考察它們對反演算法的影響,取乳房組織的電特性參數(shù)在正常值的1±10%倍范圍內(nèi)均勻隨機(jī)變化[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)分別表示相應(yīng)電參數(shù)的真實(shí)分布以及在迭代次數(shù)分別為k=0、25、50時的重建結(jié)果。為了檢驗(yàn)成像方法檢測腫瘤的效果,各子圖(d)中附加的黑色圓周表示真實(shí)腫瘤的輪廓。

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

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

4.結(jié) 論

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

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

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

[1] 劉廣東,張業(yè)榮.An overview of active microwave imaging for early breast cancer detection[J].南京郵電大學(xué)學(xué)報(bào)(自然科學(xué)版),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].微波學(xué)報(bào),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è)榮,聶在平,阮穎錚.用變形玻恩迭代法反演電導(dǎo)率的二維非均勻分布[J].電子學(xué)報(bào),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] 徐偉強(qiáng),王敏錫.用求解二維導(dǎo)體柱電磁逆散射[J].電波科學(xué)學(xué)報(bào),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] 王文珂,黃卡瑪,劉智明,等.遺傳算法在生物組織電導(dǎo)率重構(gòu)中的應(yīng)用研究[J].電波科學(xué)學(xué)報(bào),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] 鄧小波,聶在平,趙延文,等.利用相位感應(yīng)測井?dāng)?shù)據(jù)重建電導(dǎo)率和介電常數(shù)[J].電波科學(xué)學(xué)報(bào),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] 傅 林,黃卡瑪,向勝昭.生物組織磁聚焦電導(dǎo)率成像原理及反演算法[J].電波科學(xué)學(xué)報(bào),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].物理學(xué)報(bào),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)