基于域积分方程的对比源反演算法在三维微波成像中的应用

谢玉芯 缪竟鸿 王学静

（天津工业大学电子与信息工程学院 天津 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.