基于域積分方程的對比源反演算法在三維微波成像中的應(yīng)用

謝玉芯 繆竟鴻 王學(xué)靜

（天津工業(yè)大學(xué)電子與信息工程學(xué)院 天津 300387）

1 引言

微波具有極強的穿透性，因此三維微波成像在無損探傷、醫(yī)學(xué)成像、地表探測、地球物理勘探等實際工程領(lǐng)域中的應(yīng)用具有巨大的潛力。眾所周知，微波成像是逆問題，具有非線性和不確定性，難以求解，因此迫切需要精確的成像算法。算法可以分為線性和非線性算法，目前實際應(yīng)用領(lǐng)域中采用的成熟算法多是線性算法，如合成孔徑聚焦（Synthetic Aperture Focussing Technique，SAFT）算法[1]和合成孔徑雷達（Synthetic Aperture Radar，SAR）算法[2,3]。由于采用了近似，線性算法重建圖像的分辨率和精確度是非常有限的，而且從圖像中無法獲得興趣目標的物理（材料）參數(shù)，如介電常數(shù)、電導(dǎo)率等。非線性算法，如對比源反演（Contrast Source Inverison，CSI）算法，則可以在迭代過程中逐步改善圖像，如興趣目標的位置、形狀及材料參數(shù)，最后取得與興趣目標最近似的圖像。正在進行的國家 863計劃——“繞月探測工程科學(xué)數(shù)據(jù)應(yīng)用與研究”就明確提出了采用非線性（優(yōu)化）算法對月球表面探測數(shù)據(jù)進行處理的研究，所以針對非線性算法的研究是非常必要的。對比源反演算法無須正演計算，在迭代過程中采用快速傅里葉（Fast Fourier）變換計算并矢格林（Dyadic Green）函數(shù)算子及其共軛算子，確保了反演過程的高效率及穩(wěn)定性[4,5]。驗證算法是否有效，通常采用仿真（理想）數(shù)據(jù)來檢驗。如果要進一步驗證算法在微波成像中應(yīng)用的可行性則需要對實測（非理想）數(shù)據(jù)成像，這種檢驗手段是非常必要的。文獻[6]已經(jīng)驗證了對比源反演在二維微波成像中應(yīng)用的可行性，本文通過該算法對三維微波實測數(shù)據(jù)[7]成像來檢驗算法的有效性。

2 實測數(shù)據(jù)的試驗設(shè)置及測量方式

微波實測數(shù)據(jù)集收集在表1中，共有2個電介質(zhì)目標：“twocubes”和“twospheres”。

表1 實測數(shù)據(jù)集Tab.1 The measured datasets

如圖1所示，在以均質(zhì)、各向同性的媒質(zhì)為背景的立方體“興趣區(qū)域”D中存在一個目標，發(fā)射天線和接收天線被放置在距D中心1.796m的球面上，分別形成圓環(huán)形區(qū)域T和M。對于發(fā)射天線，方位角φt在的0o～320o范圍內(nèi)以每次移動40o的方式變化，極角θt在 30o～150o的范圍內(nèi)以每次移動15o的方式變化，共有81個放射天線的位置。對于接收天線，極角θr=90o保持不變，即接收天線位于z=0的平面上，方位角φr和放射天線的方位角φt的距離不能小于 5 0o，因此接收天線的方位角φr=φt+40°+n×10°在 0°～350°的范圍內(nèi)以 360°為周期變化，因此對于每個發(fā)射天線，都有 27個接收天線與之對應(yīng)。接收天線都沿極角θr方向極化，即接收天線所測量的微波實測數(shù)據(jù)為Eθ(R, 90°,φ)，是沿極角θr方向極化的電場，在z=0（θr=90°）的平面上，極角θr方向又與-z方向相同，所以測量的微波（電場）實測數(shù)據(jù)的方向也可以看做是沿-z方向的。這些興趣目標都是由電介質(zhì)材料構(gòu)成的，對于這一類型的目標，本文中的算法用來重建這些目標的位置、形狀及介電常數(shù)的對比度值。

圖1 實測數(shù)據(jù)的測量方式Fig.1 The measurement setup

3 對比源反演算法

在三維電磁逆散射問題中，電場和磁場均為矢量，即在直角坐標系下的x、y和z方向上都有分量。介電常數(shù)、電導(dǎo)率及磁導(dǎo)率等電磁參數(shù)會影響電磁場的分布，在本文中忽略磁導(dǎo)率非均質(zhì)的影響，只考慮電場分布的問題。文中首先介紹在三維空間下的域積分方程，這些域積分方程被稱為電場的“數(shù)據(jù)”方程和“目標”方程[4]，表達式為

式中，R′和R分別代表了源點和場點，E表示總的電場強度，上標“sc”和“in”分別代表了散射場和入射場，G表示在真空中電磁波的三維格林函數(shù)，k為波數(shù)，x為對比度值函數(shù)（contrast function），即興趣目標與背景媒質(zhì)的電磁參數(shù)的比值，可被定義為

式中，0ε為真空的介電常數(shù)；εr(R′)和rε分別為目標和背景媒質(zhì)的相對介電常數(shù)；σe(R′)為目標的電導(dǎo)率，所以可以看出實數(shù)部分Re()χ代表了介電常數(shù)的對比度值，虛數(shù)部分Im()χ與角頻率ω有關(guān)，代表了電導(dǎo)率與背景媒質(zhì)介電常數(shù)的對比度值。本文中的興趣目標是由電介質(zhì)材料構(gòu)成的，因此重建對比度值的實數(shù)部分應(yīng)該遠遠大于其虛數(shù)部分。

為了能夠在算法中緊湊、簡潔地表示“數(shù)據(jù)”和“目標”方程，就需要用積分算子來表示積分方程。文中首先定義“對比源”，也就是對比度值和總場的積，即

4 重建結(jié)果

根據(jù)上述試驗設(shè)置的實際情況，式（1）中的“數(shù)據(jù)”方程要進行相應(yīng)的修正，即微波（電場）數(shù)據(jù)只沿-z方向極化。如表1所示，實測數(shù)據(jù)包含多個頻率（從3～8GHz，共21個頻率），文中用對比源反演算法對多個頻率成像，以便結(jié)合不同頻率數(shù)據(jù)包含的信息，取得更好的重建結(jié)果。常用的一種方法是依次用算法對這些頻率下的數(shù)據(jù)進行重建，且用前一個頻率的結(jié)果（對比源w對比度值x）作為算法對下一個頻率數(shù)據(jù)重建的初始值，稱這種方法為頻率跳變（Frequency Hopping，F(xiàn)H）方法。使用表2中收集的反演參數(shù)，可以得到下文的重建結(jié)果。對于所有的重建結(jié)果，興趣區(qū)域D均被離散為32×32，即32個子域。

表2 實測數(shù)據(jù)的反演參數(shù)Tab.2 The inversion parameters of the measured dataset

“twocubes”如圖 2所示，目標是由兩個介電常數(shù)為2.3（相對于背景媒質(zhì)，其介電常數(shù)的對比度值x=1.3），邊長均為25mm的電介質(zhì)立方體構(gòu)成，兩個立方體的中心位置分別為（12.5，-12.5，37.5）mm和（-12.5，12.5，62.5）mm。在這興趣區(qū)域D的體積為 1 12× 1 12× 1 12mm3，中心位置為（0，0，50）mm。首先用對比源反演算法及頻率跳變方法對目標“twocubes”進行重建，選用 3、4、5、6、7和8GHz共6個頻率的數(shù)據(jù)，算法對每個頻率數(shù)據(jù)的迭代次數(shù)均為100。

圖2 目標“twocubess”的實際結(jié)構(gòu)Fig.2 The configuration of “twocubes”

圖3顯示了采用上述處理方法后所取得的重建結(jié)果的實數(shù)部分，圖4顯示了重結(jié)果的虛數(shù)部分。從上述圖中可以看到重建結(jié)果的實部（介電常數(shù)）的對比度值Re()χ遠遠大于虛部（電導(dǎo)率）的對比度值Im()χ。根據(jù)式（3），可以判斷出興趣目標是電介質(zhì)目標。

圖3 目標“twocubes”的重建結(jié)果的實數(shù)部分Fig.3 The real part of the inversion results of “twocubes”

圖4 目標“twocubes”的重建結(jié)果的虛數(shù)部分Fig.4 The imaginary part of the inversion results of“twocubes”

從圖3中還可以看出興趣目標在xy平面上的重建形狀是一個方形，重建的最大對比度值Re()χ約為 1.2，略低于實際的對比度值。如圖 3c和圖 3g所示，對于z方向上的不同位置，方形的位置發(fā)生了明顯的變化，方形的尺寸卻沒有大的變化，所以可以判斷出興趣目標是兩個在z方向上不同位置的立方體。圖 3a～圖 3d顯示了中心位置是（12.5，-12.5，37.5）mm的立方體，圖3e顯示了兩個立方體接觸的部分，圖 3f～圖 3i顯示了中心位置是（-12.5，12.5，62.5）mm 的立方體。圖 5和圖 6顯示了3～8GHz頻率數(shù)據(jù)在z=37.5mm和z=62.5mm位置的重建結(jié)果，隨著頻率的升高，可以看到重建的兩個立方體的位置，形狀變得越來越清楚，重建的介電常數(shù)的對比度值Re()χ也越來越準確。

圖5 目標“twocubes”在z=37.5mm位置的重建結(jié)果Fig.5 The inversion results of “twocubes” for z=37.5mm

圖6 目標“twocubes”在z=62.5mm位置的重建結(jié)果Fig.6 The inversion results of “twocubes” for z=62.5mm

“Twospheres”如圖 7所示，該目標由介電常數(shù)為 2.6（對比度值x=1.6），直徑為 50mm的兩個電介質(zhì)球體構(gòu)成。兩個球體沿著x軸排列，其中心位置分別為（-25，0，0）mm和（25，0，0）mm。因為“Twospheres”的尺寸較大，此時成像區(qū)域D的體積為 140×140×140mm3，中心位置為（0，0，0）mm。同樣采用對比源反演算法及頻率跳變方法對目標“Twospheres”進行重建，其中選用3，3.25，3.5，3.75和4GHz共5個頻率的數(shù)據(jù)，算法對每個頻率數(shù)據(jù)的迭代次數(shù)均為100。

圖7 目標“twospheres”的實際結(jié)構(gòu)Fig.7 The configuration of “twospheres”

圖8和圖9分別顯示了重建結(jié)果的實數(shù)和虛數(shù)部分。如圖所示，與重建結(jié)果的實部（介電常數(shù)）的對比度值Re()χ比較，虛部（電導(dǎo)率）的對比度值Im()χ接近于0，可以判斷出興趣目標也是電介質(zhì)目標。在圖8中，可以看到興趣目標在xy平面上的重建形狀是兩個圓，其位置和大小都重建得非常好，重建的最大對比度值Re()χ約為1.6，與實際的對比度值相符。對于z方向上不同位置，兩個圓的尺寸也相應(yīng)的變化，在圖8顯示的所有重建結(jié)果當(dāng)中，圖8e是在z=2.187 5mm的位置，重建的兩個圓的尺寸是最大的，其他位置上兩個圓的尺寸則相應(yīng)地減小，所以可以判斷出興趣目標是兩個球體。圖 10依次顯示了上述頻率數(shù)據(jù)在z=0位置的重建結(jié)果，圖10a中重建的最大的對比度值Re()χ約為1.1，圖10c中最大對比度值Re()χ約為1.4，圖10e中重建最大對比度值Re()χ約為1.6。隨著頻率的升高，可以看到重建的兩個球體的位置，形狀變得越來越清楚，重建的介電常數(shù)的對比度值Re()χ也越來越接近目標的實際對比度值。

圖8 目標“twospheres”的重建結(jié)果的實數(shù)部分Fig.8 The real part of the inversion results of “twospheres”

圖9 目標“twospheres”的重建結(jié)果的虛數(shù)部分Fig.9 The imaginary part of the inversion results of“twospheres”

圖10 目標“twospheres”在z=0位置的重建結(jié)果Fig.10 The inversion results of “twspheres” for z=0

5 結(jié)論

本文主要提出了在對比源反演算法基礎(chǔ)上的擴展算法（頻率跳變方法）對多頻、多收發(fā)的三維實測微波數(shù)據(jù)成像的具體方法。從上述重建結(jié)果中可以看到，采用了頻率跳變方法的對比源反演算法對上述目標都取得了非常好的重建結(jié)果，目標的重建位置和形狀非常準確，且重建的介電常數(shù)的對比度值也與實際的對比度值接近，這些重建結(jié)果驗證了對比源反演算法的精確性及在三維微波成像中應(yīng)用的可行性。

致謝：本文使用的實測數(shù)據(jù)由法國馬賽菲涅耳（Fresnel）研究所微波研究小組所提供，在此作者向 Jean-Michel Geffrin和 Pierre Sabouroux表示感謝。

[1] J Miao, Y Xie, J Li, et al. Application of linear algorithms for reconstructing complicated targets from experimental data[C]. Proceedings of the International Symposium on Computer Science and Society, Kota Kinabalu, Malaysia, 2011: 26-28.

[2] 任笑珍, 楊汝良. 機載前視SAR三維成像算法研究[J]. 電子與信息學(xué)報, 2010, 32(6): 1361-1365.Ren Xiaozhen, Yang Ruliang. Study on threedimensional imaging algorithm for airborne forwardlooking SAR[J]. Journal of Electronics and Information Technology, 2010, 32(6): 1361-1365.

[3] 王金峰, 皮亦鳴, 曹宗杰. 一種機載 SAR層析三維成像算法[J]. 電子與信息學(xué)報, 2010, 32(5):1029-1033.Wang Jinfeng, Pi Yiming, Cao Zongjie. An algorithm for airborne SAR tomography 3D imaging[J]. Journal of Electronics and Information Technology, 2010,32(5): 1029-1033.

[4] 何勁, 羅迎, 張群, 等. 隨機線性調(diào)頻步進雷達波形設(shè)計及成像算法研究[J]. 電子與信息學(xué)報, 2011,33(9): 2068-2075.He Jin, Luo Ying, Zhang Qun, et al. Waveform design and imaging algorithm research of random frequency stepped chirp signal ISAR[J]. Journal of Electronics and Information Technology, 2011, 33(9): 2068-2075.

[5] A Abubarkar, P M van den Berg. Iterative forward and inverse algorithms based on domain integral equations for three dimensional electric and magnetic objects[J]. Journal of Computatinal Physics, 2004,195(1): 236-262.

[6] P Barrière, J Idier, J Laurin, et al. Contrast source inversion method applied to relatively high contrast objects[J]. Inverse Problems, 2011, 27(7): 075012.

[7] J M Geffrin, P Sabouroux, Continuing with the Fresnel database: experimental setup and improvements in 3D scattering measurement[J].Inverse Problems, 2009, 25: 024001.

[8] J Miao. Linear and nonlinear inverse scattering algorithms applied in 2-D electromagnetics and elastodynamics[M]. Kassel: Kassel University Press,2008.

[9] R E Kleinman, P M van den Berg. A contrast source inversion method[J]. Inverse Problems, 1997, 13(6):1607-1620.

[10] A Zakaria, C Gilmore, J Lovetri. Finite-element contrast source inversion method for microwave imaging[J]. Inverse Problems, 2010, 26(1): 115010.