微波成像技術(shù)及其算法綜述

許 會，陳艷玲

（沈陽工業(yè)大學(xué) 先進(jìn)在線檢測技術(shù)省重點(diǎn)實(shí)驗(yàn)室，沈陽 110870）

近些年來，微波成像理論在科學(xué)、技術(shù)方面的重要作用已經(jīng)受到越來越多的關(guān)注。微波成像技術(shù)以微波在各種復(fù)雜媒質(zhì)中的傳播和散射的研究為基礎(chǔ)，從待測目標(biāo)的散射場反演目標(biāo)的幾何形狀、空間位置、電磁特性等參數(shù)。微波成像技術(shù)作為電磁逆散射技術(shù)的重要內(nèi)容，可實(shí)現(xiàn)對目標(biāo)的無損檢測，能同時對目標(biāo)進(jìn)行幾何成像和物理成像，且具有較高的分辨率，因而微波成像技術(shù)在地球物理探測、生物醫(yī)學(xué)成像、環(huán)境監(jiān)測等領(lǐng)域具有重要用途［1－2］。但是由于背景介質(zhì)復(fù)雜、未知量大以及反演過程病態(tài)等問題，該項(xiàng)研究一直是電磁學(xué)領(lǐng)域的重要課題。筆者閱讀了近5年來的80余篇國內(nèi)外的文獻(xiàn)，發(fā)現(xiàn)其中針對微波成像細(xì)節(jié)方面，如參數(shù)估計方法、測量方案、優(yōu)化準(zhǔn)則等研究［3］的文獻(xiàn)占了總數(shù)的10%；算法研究［4－14］的文獻(xiàn)將近占據(jù)了40%；微波成像的應(yīng)用方面的文獻(xiàn)則占據(jù)50%左右，尤其是在生物醫(yī)學(xué)成像方面如對乳腺癌的檢測等應(yīng)用方面的文獻(xiàn)有著很大的比重［15－16］。盡管應(yīng)用于實(shí)踐是理論研究的最終目的，但是逆算法作為微波成像的重要內(nèi)容，起著關(guān)鍵的作用。文章針對微波成像算法問題進(jìn)行討論，對現(xiàn)存的微波成像算法做了簡要的概述，并對幾種具有代表性的成像算法的原理、特點(diǎn)及使用范圍進(jìn)行介紹。

微波成像技術(shù)的主要任務(wù)是求解微波照射被測媒質(zhì)時的逆散射問題，通過測量被測媒質(zhì)外部的散射場數(shù)據(jù)，重建被測媒質(zhì)內(nèi)部的復(fù)介電常數(shù)圖像。被測的散射場攜帶大量有關(guān)散射體的信息，利用關(guān)于散射目標(biāo)的先驗(yàn)知識，經(jīng)過適當(dāng)?shù)臄?shù)學(xué)處理之后可以提取出散射體本身所具有的某些特性，如散射體的形狀、介電常數(shù)的分布等。微波成像技術(shù)一般包含兩部分：一部分是正演數(shù)值模擬，即計算給定模型的電磁場分布；另一部分是逆過程，即根據(jù)給定的測量場重構(gòu)電參數(shù)的分布［1－4］。此技術(shù)可以檢測非金屬或金屬材料內(nèi)部缺陷的大小、形狀、位置以及物理檢測化學(xué)變化過程，在物體的無損傷性檢測其內(nèi)部缺陷的應(yīng)用中能起很大的作用。過去的幾十年里，在微波成像和電磁逆散射技術(shù)的相關(guān)領(lǐng)域內(nèi)應(yīng)用逆算法在效率、魯棒性、有效性方面取得了很大的進(jìn)展。

對于微波成像中正問題的求解可使用解析法或數(shù)值法進(jìn)行，其中解析法僅適用于少數(shù)典型物體的散射場的計算，而數(shù)值法可計算復(fù)雜形狀物體的散射。而對于逆問題，許多的算法被用來反演測得的散射數(shù)據(jù)。根據(jù)不同算法的特點(diǎn)，有三種分類方式［5－8］，如表1所示。下面對電磁逆散射的圖像重建問題進(jìn)行討論，介紹幾種具有代表性的成像算法。

1 ω－k算法

ω－k算法［10］是一種典型的譜域重建算法。與傳統(tǒng)的合成孔徑成像算法相比，該算法具有更高的精度和計算速度，適用于均勻散射背景下的成像，如飛機(jī)降落時對不明物體的偵查、對寬測繪帶星載SAR數(shù)據(jù)進(jìn)行精確成像處理等。ω－k算法公式為：

式中E（xn，ωk）表示在測試點(diǎn)xn處，頻率為ωk時的電場值，其中n＝1，2，3…，N，k＝1，2，3…，kmax；O（x，y）為目標(biāo)物體函數(shù)，（x，y）為成像點(diǎn)的坐標(biāo)；v為波在介質(zhì)中傳播的速度。

ω－k算法流程如圖1所示。算法步驟如下：

圖1 ω－k算法流程框圖

（2）作Y方向一維傅里葉變換。

表1 微波成像算法分類列表

（3）進(jìn)行空間移位、插值等處理，其中（x0，y0）為目標(biāo)的中心坐標(biāo)。

（4）將處理過的信號作二維傅里葉反變換，得到的幅度為空間分布圖像。

（5）對反變換后的矩陣中的復(fù)數(shù)數(shù)據(jù)元素逐個取模。

2 局部形狀函數(shù)算法

局部形狀函數(shù)算法［11］（LSF）是一種典型的空間域非線性迭代算法。LSF是一種充分的、一般的方法，可用于包含多個強(qiáng)散射體和任意形狀、大小的散射體的情況。由于該方法原理簡單，適用于強(qiáng)散射體的成像，故多被用來對金屬散射目標(biāo)進(jìn)行重建。但散射矩陣T的求取也成為該算法的一個難點(diǎn)。

如圖2所示，LSF方法是通過將散射問題轉(zhuǎn)化成求解二進(jìn)制函數(shù)問題來將非線性問題線性化。散射區(qū)域V被離散化為N 個單元，每一單元用Vi表示，i＝1，2，3…，N。假設(shè)子區(qū)域Vi劃分的足夠小且每一單元Vi里具有相同的介電常數(shù)和電導(dǎo)率，S為金屬散射體所在的區(qū)域。

圖2 LSF算法的二維散射模型

則二進(jìn)制形狀函數(shù)γi如下所示：

這里假定γi為二進(jìn)制變量，取值為0或1，含有金屬散射體的區(qū)域的形狀函數(shù)設(shè)為1，不含金屬散射體的區(qū)域的形狀函數(shù)設(shè)為0。但是在實(shí)際迭代和求解最優(yōu)值時，需要將γi設(shè)定為一個在0到1之間連續(xù)變化的變量，然后通過逆算法得出關(guān)于γi的圖像。則LSF算法的重建步驟如下：

（1）對參數(shù)設(shè)置初值。求取散射域內(nèi)的已知的總場分布，得到目標(biāo)物體函數(shù)γi的分布作為初始分布，或者假設(shè)初始值由其它先驗(yàn)知識提供。

（2）求解重建的目標(biāo)函數(shù)，可以用矩量法或其它算法如時域有限差分法（FDTD）進(jìn)行正向問題的計算，得到檢測域和成像域的總場分布。

（3）比較由重建的目標(biāo)函數(shù)得到的檢測場分布和實(shí)際的檢測場分布，若兩者相當(dāng)接近，低于某種準(zhǔn)則所規(guī)定的誤差可停止迭代過程，否則，用新的計算結(jié)果更新成像域分布，重新迭代至滿足誤差限。

（4）重建二進(jìn)制形狀函數(shù)γi來構(gòu)建散射體圖像。

3 波恩近似迭代算法

波恩近似迭代算法［12］（BIM）屬于一種典型的傳統(tǒng)型的空間域非線性迭代算法。此算法在每一次迭代中均需進(jìn)行正演數(shù)值模擬計算和線性化過程，適用于對電介質(zhì)和導(dǎo)電媒質(zhì)目標(biāo)的重建。當(dāng)被測媒質(zhì)介電常數(shù)與背景介質(zhì)的介電常數(shù)相比差別不大于10∶1的時候，可以利用波恩近似來簡化逆散射問題，不適用于對強(qiáng)散射體成像的計算且反演速度較慢，精度低。變形波恩迭代法（DBIM）是對BIM的改進(jìn)。與BIM相比，DBIM收斂快，但抗隨機(jī)噪聲的能力低，數(shù)值實(shí)現(xiàn)不方便，只適用于小尺寸、低對比度目標(biāo)的反演，且精度不高。

設(shè)入射場為單位TM波，則電磁散射方程為：

式中ω為角頻率；μ0為自由空間的磁導(dǎo)率；r＝（x，y）為場點(diǎn)；r′為源點(diǎn)；O（r′）為目標(biāo)函數(shù)；為0階第二類Hankel函數(shù)。Es為散射場，總場為入射場和散射場之和，如式（4）所示：

式中Ei為入射場，E為總場。將式（4）代入式（3），得：

從上面過程中得到波恩近似迭代法的反演方程。

4 對比源算法

對比源反演算法［6－7］（CSIM）是一種典型的改進(jìn)型算法，將反演問題轉(zhuǎn)化為求解成本泛函的極小值問題，從而形成重構(gòu)對比源和對比度的迭代序列。該算法由于無須正演計算，亦無須人為地選擇正則化參數(shù)，所以反演過程穩(wěn)定，適用于均勻背景介質(zhì)和層狀背景介質(zhì)中。乘法正則化－對比源算法［7］（multiplicative regularized contrast source inversion method，MR－CSIM）是在CSIM上的改進(jìn)，能夠得到更好的分辨率和精確度。

在圖3的散射模型中，S為散射體所在的區(qū)域，T為成像區(qū)域，V為散射區(qū)域。若用E表示總場，Ei表示入射場，Es表示散射場，則總場可以如式（4）所示。

圖3 CSIM法的二維散射模型

設(shè)p，q為X－Y 平面內(nèi)的位置向量（x，y），n代表發(fā)射源的序號（因?yàn)樵谀嫔⑸涞臏y量中需要多個發(fā)射源），背景的復(fù)介電常數(shù)為εb，則電參數(shù)對比度χ可以寫為式（6），對比源w則如式（7）所示。

則該算法成本泛函的數(shù)據(jù)方程和目標(biāo)方程的線性組合為：

式中ω為角頻率；μ0為自由空間的磁導(dǎo)率；G（p，q）為關(guān)于背景介質(zhì)的格林函數(shù)。

式（8）中右側(cè)第一項(xiàng)為度量數(shù)據(jù)方程的歸一化誤差，第二項(xiàng)為度量目標(biāo)方程的歸一化誤差。當(dāng)wn＝0時，這兩項(xiàng)均為1，且因?yàn)槟繕?biāo)方程的約束，CSIM反演問題的解的非惟一性大大降低。分別對對比源wn和電參數(shù)對比度χ進(jìn)行更新，輪流依次迭代，可以實(shí)現(xiàn)對式（8）的極小化。

5 粒子群優(yōu)化算法

粒子群優(yōu)化算法［13］（PSO）是一種啟發(fā)式全局優(yōu)化技術(shù)，一種基于群智能的進(jìn)化計算技術(shù)，由Eberhart博士和Kennedy博士發(fā)明，源于對鳥群捕食的行為研究。PSO同遺傳算法類似，是一種基于迭代的優(yōu)化工具。系統(tǒng)初始化為一組隨機(jī)解，通過迭代搜尋最優(yōu)值。但是并沒有遺傳算法用的交叉以及變異。通過利用個體間的協(xié)作和競爭來實(shí)現(xiàn)對問題最優(yōu)解的搜索。PSO的優(yōu)勢在于簡單容易實(shí)現(xiàn)，采用實(shí)數(shù)求解并且需要調(diào)整的參數(shù)較少，是一種通用的全局搜索算法。

PSO初始化為一群隨機(jī)粒子（隨機(jī)解），然后通過迭代找到最優(yōu)解。假設(shè)在一個D維的搜索空間中，隨機(jī)初始化m個粒子組成一個群落，其中第j個粒子的空間位置為xj＝（xj1，xj2，…，xjD），j＝1，2，…，m。每個粒子在搜索空間中以一定的速度飛行，速度定義為vj＝（vj1，vj2，…，vjD），可根據(jù)經(jīng)驗(yàn)進(jìn)行動態(tài)調(diào)整。在每一次迭代中，粒子通過跟蹤兩個“極值”來更新自己。第一個就是粒子本身所找到的最優(yōu)解，記為gj＝（gj1，gj2，…，gjD），這個解叫做個體極值，另一個極值是整個種群目前找到的最優(yōu)解，記為zj＝（zj1，zj2，…，zjD），這個極值是全局極值。在找到兩個極值時，粒子可根據(jù)如下的迭代公式來更新自己的速度和位置。

式中1≤j≤m；1≤d≤D；u是介于0～1之間的隨機(jī)數(shù)；c1，c2是學(xué)習(xí)因子，通常c1＝c2＝2。每一維粒子的速度都會被限制在一個最大速度vmax內(nèi)，當(dāng)vmax較大時，粒子飛行速度大，有利于全局搜索，但有可能錯過最優(yōu)解；當(dāng)vmax較小時，粒子在特定區(qū)域內(nèi)精細(xì)搜索，但容易陷入局部最優(yōu)。PSO的算法流程如圖4所示。

6 非精確牛頓算法

非精確牛頓算法［14］（Inexact Netwon Method，INM）是一種典型的確定性算法，由兩個嵌套的循環(huán)組成。外面的循環(huán)是使非線性散射方程線性化，里面的循環(huán)則是運(yùn)用截斷的Landweber算法求解已經(jīng)被線性化了的方程。該算法適用于包含強(qiáng)散射體的電磁逆散射問題。該算法的迭代步驟為：

圖4 PSO算法的流程圖

（1）給定初始值x0，置n＝0。

（2）通過計算Fre′chet導(dǎo)數(shù)，進(jìn)行線性化，從而得到線性方程。

（3）運(yùn)用截斷的Landweber算法，進(jìn)行正則化操作，計算有限的序列hn。

（4）通過設(shè)置xn＋1＝xn＋hn來更新數(shù)據(jù)，且n＝n＋1。

（5）若║xn＋1－xn║≤δ，δ為終止條件，則迭代終止。否則，置n＝n＋1，轉(zhuǎn)到步驟（2）繼續(xù)進(jìn)行迭代。

7 結(jié)論

除了上述介紹的方法外，微波成像還有很多算法，并且不斷有新算法提出。隨著人們對實(shí)際問題的重視和研究的深入，智能算法開始越來越多地應(yīng)用于微波成像問題中，一系列的混合算法也被提出，如牛頓算法與對比源算法相結(jié)合、線性采樣方法與蟻群優(yōu)化算法相結(jié)合等。同時，研究人員也開始把目光更多地關(guān)注于實(shí)際數(shù)據(jù)的重構(gòu)。在未來時期內(nèi)，智能算法和混合算法將成為解決微波成像問題的主流算法，微波成像技術(shù)的研究將不斷得到完善，其應(yīng)用也將更加廣泛。

［1］Hu G H，Qu F L，Zhang B.Direct and inverse problems for electromagnetic scattering by a doubly periodic structure with a partially coated dielectric［J］.Mathematical Methods in the Applied Sciences，2010，33（2）：147－156.

［2］Cakoni F.Recent developments in the qualitative approach to inverse electromagnetic scattering theory［J］.Journal of Computational and Applied Mathematics，2007，204（2）：242－255.

［3］Winters D W，Van Veen B D，Hagness S C.A sparsity regularization approach to the electromagnetic inverse scattering problem［J］.IEEE Transactions on Antennas and Propagation，2010，58（1）：145－154.

［4］YU Chun，YUAN Meng－qing，ZHANG Yang－jun，et al.Microwave imaging in layered Media：3－D image reconstruction from experimental data［J］.IEEE Transactions on Antennas and Propagation，2010，58（2）：440－448.

［5］Dubois A，Belkebir K，Catapano I，et al.Iterative solution of the electromagnetic inverse scattering problem from the transient scattered field［J］.Radio Science，2009，44（1）：199－201.

［6］Gilmore C，Mojabi P，Lovetri J.Comparison of an enhanced distorted born iterative method and the multiplicative－regularized contrast source inversion method［J］.IEEE Transactions on Antennas and Propagation，2009，57（8）：2341－2351.

［7］Bozza G，Pastorino M.An inexact newton－based approach to microwave imaging with in the contrast source formulation［J］.IEEE Transactions on Antennas and Propagation，2009，57（4）：1122－1132.

［8］Pastorino M.Stochastic optimization methods applied to microwave imaging：areview［J］.IEEE Transactions on Antennas and Propagation，2007，55（3）：538－548.

［9］Brignone M，Bozza G，Randazzo A，et al.A hybrid approach to 3Dmicrowave imaging by using linear sampling and ACO［J］.IEEE Transactions on Antennas and Propagation，2008，56（10）：3224－3232.

［10］Cadalli N，Munson D C.A simulation study of theωk SAR algorithm squinted case with application to runway imaging［C］.Proceedings of the IEEE Interational Conference on Application Specific Standard Product.Istanbul，Turkey：2000.

［11］Weedon W H，Chew W C.Time domain inverse scattering using the local shape function（LSF）method［J］.Inverse Problems，1993，（9）：551－564.

［12］Moghaddam M，Chew W C.Comparison of the born iterative method and tarantola's method for an electromagnetic time－domain inverse problem［J］.Interational Journal of Imaging Systems and Technology，1994（3）：318－333.

［13］Donelli M.，F(xiàn)ranceschini G.，Martini A.，et al.An integrated multiscaling strategy based on a particle swarm algorithm for inverse scattering problems［J］.IEEE Transactions on Geoscience and Remote Sensing，2006，44（2）：298－311.

［14］Bozza G，Estatico C，Pastorino M，et al.An inexact newton method for microwave reconstruction of strong scatterers［J］.IEEE Antennas and Wireless Propagation Letters，2006，5（1）：61－64.

［15］Golnabi A H，Meaney P M，Geimer S，et al.Micro－wave Imaging for Breast Cancer Detection and Therapy Monitoring［C］.IEEE Topical Conference on Biomedical Wireless Technologies，Networks，and Sensing Systems.Phoenix，Arizona，USA，2011.

［16］Ashtari A.Using apriori information for regularization in breast microwave image reconstruction［J］.IEEE Transactions on Biomedical Engineering，2010，57（9）：2197－2208.