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

謝玉芯 繆竟鴻 王學靜

（天津工業大學電子與信息工程學院 天津 300387）

1 引言

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

2 實測數據的試驗設置及測量方式

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

表1 實測數據集Tab.1 The measured datasets

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

圖1 實測數據的測量方式Fig.1 The measurement setup

3 對比源反演算法

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

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

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

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

4 重建結果

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

表2 實測數據的反演參數Tab.2 The inversion parameters of the measured dataset

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

圖2 目標“twocubess”的實際結構Fig.2 The configuration of “twocubes”

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

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

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

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

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

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

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

圖7 目標“twospheres”的實際結構Fig.7 The configuration of “twospheres”

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

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

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

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

5 結論

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

致謝：本文使用的實測數據由法國馬賽菲涅耳（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]. 電子與信息學報, 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]. 電子與信息學報, 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] 何勁, 羅迎, 張群, 等. 隨機線性調頻步進雷達波形設計及成像算法研究[J]. 電子與信息學報, 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.