許 會(huì)，陳艷玲

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

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

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

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

1 ω－k算法

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

式中E（xn，ωk）表示在測(cè)試點(diǎn)xn處，頻率為ωk時(shí)的電場(chǎng)值，其中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）將處理過(guò)的信號(hào)作二維傅里葉反變換，得到的幅度為空間分布圖像。

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

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

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

如圖2所示，LSF方法是通過(guò)將散射問題轉(zhuǎn)化成求解二進(jìn)制函數(shù)問題來(lái)將非線性問題線性化。散射區(qū)域V被離散化為N 個(gè)單元，每一單元用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)值時(shí)，需要將γi設(shè)定為一個(gè)在0到1之間連續(xù)變化的變量，然后通過(guò)逆算法得出關(guān)于γi的圖像。則LSF算法的重建步驟如下：

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

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

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

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

3 波恩近似迭代算法

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

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

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

式中Ei為入射場(chǎng)，E為總場(chǎng)。將式（4）代入式（3），得：

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

4 對(duì)比源算法

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

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

圖3 CSIM法的二維散射模型

設(shè)p，q為X－Y 平面內(nèi)的位置向量（x，y），n代表發(fā)射源的序號(hào)（因?yàn)樵谀嫔⑸涞臏y(cè)量中需要多個(gè)發(fā)射源），背景的復(fù)介電常數(shù)為εb，則電參數(shù)對(duì)比度χ可以寫為式（6），對(duì)比源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時(shí)，這兩項(xiàng)均為1，且因?yàn)槟繕?biāo)方程的約束，CSIM反演問題的解的非惟一性大大降低。分別對(duì)對(duì)比源wn和電參數(shù)對(duì)比度χ進(jìn)行更新，輪流依次迭代，可以實(shí)現(xiàn)對(duì)式（8）的極小化。

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

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

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

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

6 非精確牛頓算法

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

圖4 PSO算法的流程圖

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

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

（3）運(yùn)用截?cái)嗟腖andweber算法，進(jìn)行正則化操作，計(jì)算有限的序列hn。

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

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

7 結(jié)論

除了上述介紹的方法外，微波成像還有很多算法，并且不斷有新算法提出。隨著人們對(duì)實(shí)際問題的重視和研究的深入，智能算法開始越來(lái)越多地應(yīng)用于微波成像問題中，一系列的混合算法也被提出，如牛頓算法與對(duì)比源算法相結(jié)合、線性采樣方法與蟻群優(yōu)化算法相結(jié)合等。同時(shí)，研究人員也開始把目光更多地關(guān)注于實(shí)際數(shù)據(jù)的重構(gòu)。在未來(lái)時(shí)期內(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.