刘广东，葛新同

(1.阜阳师范学院物理与电子工程学院，安徽阜阳236037; 2.阜阳师范学院数学与统计学院，安徽阜阳236037)



时域反演德鲁色散媒质的电磁逆散射技术

刘广东1，葛新同2

(1.阜阳师范学院物理与电子工程学院，安徽阜阳236037; 2.阜阳师范学院数学与统计学院，安徽阜阳236037)

摘要:德鲁(Drude)经验模型常用于描述等离子体、金属等媒质的电色散特性．利用宽带的时域测量数据直接反演电参数，相比单频(频域)技术而言，具有信息量大、成像分辨率高的优势．时域直接反演色散媒质电参数的主要困难在于它们是频率相关的．为了克服该困难，本文提出了一种时域电磁(EM)逆散射新技术:转而同时反演德鲁模型的4类频率无关的模型参数．该技术的主要环节为: (1)描述为含正则化项的约束最小化问题; (2)转化为无约束最小化问题; (3)解析导出梯度; (4)分别利用时域有限差分(FDTD)法、共轭梯度(CG)法迭代求解正演、反演子问题．在一维(1-D)、二维(2-D)两个数值算例中，所需的测量数据也由FDTD仿真值代替，并加入了加性高斯白噪声(AWGN)．反演结果初步证实了该技术的性能．

关键词:电磁逆散射;德鲁色散媒质;正则化;时域有限差分法;共轭梯度法

1　引言

实验表明电磁场作用下自然界中大部分媒质，如生物组织、土壤、等离子体、金属等，其电磁特性均与工作频率f有关［1］．这些媒质称作色散媒质，电磁参数的频率相关性称作色散特性．拟合获得模型参数后，便可用于描述特定媒质的色散特性，如德拜(Debye)模型、洛伦兹(Lorentz)模型、德鲁(Drude)模型就是最为常见的三大类，其中德鲁(Drude)模型适用于等离子体、金属等媒质，这一大类媒质常称作德鲁色散媒质［2］．为处理这三大类色散媒质的电磁辐射、散射等问题，魏兵等人提出通用的时域有限差分(finite-difference time-domain，FDTD)方案［2］和通用的吸收边界［3］，王飞等人近来提出通用的Newmark-FDTD新方法［4］，都为研究这些色散媒质的时域逆散射(亦常称作反演、重建等)问题奠定了基础．

近年来，面向色散媒质的时域反演技术已经取得一些进展，主要有:对于Debye色散媒质，Winters等人提出同时反演光频相对介电常数ε∞、静态相对介电常数εs、静态电导率σs共3类模型参数的时域逆散射技术［5］，刘广东等人引入正则化技术对抗逆问题的病态特性［6］，Papadopoulos等人补充弛豫时间τ，将反演参数增加到4类［7］;对于Lorentz色散媒质，Papadopoulos等人提出同时反演光频相对介电常数ε∞、静态相对介电常数εs、固有频率ω0、碰撞频率ζ共4类模型参数的时域逆散射技术［8］．

然而，对于Drude色散媒质，未见到相关反演方法的报道．为此，本文探索在时域同时反演Drude色散媒质的4类模型参数，巧妙克服了时域直接反演电参数所面临的主要困难:频率相关性．该项工作有望为冶金、能源、材料、地球物理等领域提供参考．

2　问题描述

预设条件: (1)所有媒质均为线性、各向同性的无磁媒质; (2)已知问题空间V的边界，假设其间存在电色散媒质，其色散特性满足Drude经验模型［2］，但其模型参数未知(因此其电参数亦未知) ; (3)已知背景媒质的类型和电参数分布; (4)假设测量系统采用阵列天线，其发射阵共I个阵元(i =1，2，…，I)，接收阵共K个阵元(k =1，2，…，K)，已知各阵元的位置和激励源．

依次激活各发射阵元(每次一个)，产生的入射波辐照问题空间的散射目标，K个接收阵元同时测量，获得的时域总电场记为．本文的目的就是利用这些测量电场反演问题空间的电参数分布．由于前文假定目标媒质为无磁的Drude电色散媒质，其复值相对介电常数由一般的多极Drude经验模型确定［2］

其中，位置矢量r∈V，角频率ω= 2πf，j为虚数单位，W为极的总数，ε∞、ε0分别表示光频相对介电常数、真空介电常数，σs表示静态电导率，分别表示第w极的等离子体角频率、碰撞频率．该模型较为通用，适用的几个特例: (1)当W = 1时简化为单极情形; (2)当= 0时退化为常规的非色散媒质［9］; (3)当σs= 0时简化为文献［2］的(3)式．由式(1)知，Drude电色散媒质的相对介电常数εr(r，ω)、等效电导率σeff(r，ω)两类电参数均与频率f有关，在时域难以直接反演，这也是和常规(非色散)媒质的主要差别［6］．然而，同时也不难发现，Drude模型的4类模型参数ε∞(r)、σs(r)、、ζw(r)却是频率无关的，因此，本文转而反演这4类模型参数．为简化后文表述，省略自变量(r)，并定义反演参数向量，其中上角标T表示转置(后文类同)，对问题空间的每一个像素(点)，需要同时反演p1，p2，…，p2 +2W共2 +2W个未知数．

选用辅助微分方程(auxiliary differential equation，ADE)法引入Drude媒质的色散特性［10］，可以得到:当激活第i个发射阵元时，无源空间的时域电场Ei、磁场Hi、第w极色散电流满足麦克斯韦(Maxwell)方程组

和一组(w =1，2，…，W)辅助微分方程(ADEs)

其中，∂t表示对时间自变量t(文中t已被省略)的偏微分算子，为哈米尔顿(Hamilton)算子［9］，μ0为真空磁导率．可见，待反演的各电磁参数都显式呈现在前述方程中，有利于后文逆问题的求解．

3　求解逆问题

电磁逆散射问题属于一类不适定问题，求解面临的两个主要困难分别是问题的非线性和病态特性［5～9］．

3.1描述为约束最小化问题

首先，为克服电磁逆散射问题的非线性困难，本文将该问题描述为关于反演参数p的约束最小化问题，其目标泛函F满足方程:

同时满足约束条件:方程组(2)～(4)．

方程(5)的右手边:第一项描述计算电场Ei，k和测量电场之间的差异，T表示测量时间，符号中右下、上角标“2”分别表示欧氏(Euclidean)范数、平方运算;第二项为加性的吉洪诺夫(Tikhonov)正则化项［11］，以抑制逆问题的病态特性，这也是本文与文献［5，7，8］的不同之处，其中γl表示与p的第l个分量pl对应的正则化参数．

3.2转化为无约束最小化问题

其次，借助拉格朗日(Lagrange)乘子罚函数法［6］，将约束最小化问题转化为无约束最小化问题，其增广目标泛函Fa如下:

3.3变分法推导梯度

利用变分法［12］，求解约束最小化问题，即等价于求解变分方程

其中，δ为一阶变分算子．结合方程(6)，为方便表述，令

其中

3.4选取梯度算法

至此，剩下的步骤只需从已有的多种非线性梯度算法中做出选择［13］．研究显示: PRP(Polak-Ribière-Polyak)共轭梯度(conjugate gradient，CG)算法的总体性能最佳［14］，因此，本文也选取这一算法．若得到第m步迭代的反演参数pm，则下一步迭代的更新公式为［13］

其中，迭代步数指标m = 1，2，…，M，M为迭代步总数，步长αm通过求解一维的线搜索问题获得［6］，方向向量dm为

4　仿真结果与讨论

为了仿真检验本文技术的性能，后文设计2个数值算例，分别针对一维(1-D)、二维(2-D)问题，探究可能影响反演效果的下列因素: (1)问题维数; (2)天线工作模式; (3)测量视角; (4)极总数; (5)参数分布类型; (6)背景媒质类型; (7)散射强度; (8)散射体位置; (9)散射体尺寸．

补充说明: (1)激励源采用调制的高斯脉冲源，详见文献［7］的式(23) ; (2) FDTD解算器采用均匀网格，空间步长设为Δ，选取Δ兼顾数值精度和数值色散要求Δ≤c0/fmax/10，其中c0为真空光速，fmax为激励源的上限频率［10］，一维、二维问题的计算区周围分别采用5层、6层卷积完全匹配层(convolution perfectly matched layer，CPML)吸收边界［15］，选取时间步长Δt应满足(Courant Friedrichs Lewy，CFL)稳定条件，其中nD表示问题的物理维度［10］; (3)条件所限，反演所需的测量值也采用FDTD仿真替代，但离散化采用双倍精细网格，考虑噪声影响时，噪声模型选用加性高斯白噪声(additive white Gaussian noise，AWGN)，其信噪比(signal to noise ratio，SNR)设为20dB; (4)取正则化参数γl= 0.001［6］，取正则化参数γl=0时则视为无正则化项; (5)以迭代总步数为迭代终止判据，取值为M =60，并定义第m步迭代的均方根误差(mean square error，MSE) e为［9］

4.1算例1:一维(1-D)问题

一维问题的几何模型类似于文献［8］，如图1所示:在厚度d =10mm的空气(视为真空)中，厚度为4d重建区为层状分布的2极(即W = 2)德鲁色散媒质，模型参数为光滑型的类正弦分布，其中ζ2随坐标z的分布分别如图2的子图(a)～(f)的黑色实线所示，其强散射的峰值分别为8.0、2.0×10－2S．m－1、1.8×1010rad．s－1、1.6×1010rad．s－1、2.0×1011Hz、1.8×1011Hz，弱散射的峰(同向峰值型)或谷值(异向峰值型)分别为2.0、0.5×10－2S．m－1、4.5×109rad．s－1、4.0×109rad．s－1、5.0×1010Hz、4.5×1010Hz;天线系统采用两边(一维情形下视为全视角)、双站(即收、发分离)的测量模式，发射阵列的两个阵元(I = 2)均距重建区d/2，接收阵列的两个阵元(K =2)均距重建区d/4．

FDTD离散化的空间、时间步长分别取为Δz = 0.5mm、Δt =0.5Δz/c0，T =2000Δt;迭代初值取自重建区的平缓区(假定先验知晓)，分别为4.0、1.0×10－2S．m－1、9.0×109rad．s－1、8.0×109rad．s－1、10.0×1010Hz、9.0×1010Hz，分别如图2的子图(a)～(f)的“+”形标记所示;考虑噪声且应用正则化条件下，分别经过1步和60步迭代后，反演结果分别如图2的子图(a)～(f)的“×”形标记、“．”形标记所示;归一化的目标泛函F、均方根误差e随迭代步数m的变化关系分别如图2的子图(g)、(h)所示(纵坐标采用了对数形式，后文同)，其中60步迭代的误差为0.0609．

另外，为了检验本文逆散射方法对抗逆问题的病态特性和噪声的性能，并考虑篇幅所限，图3给出了经过60步迭代后，光频相对介电常数ε∞的反演结果和真实分布的对比，子图(a)～(d)分别对应四种情形:无噪声且无正则化、有噪声且无正则化、无噪声且有正则化、有噪声且有正则化，反演误差依次为0.0658、0.0749、0.0347、0.0623．

敛散性分析:从图2的子图(g)、(h)可以看出，运用本文的迭代技术处理光滑型一维问题，算法收敛，这得益于正则化技术，抑制噪声影响和逆问题的病态特性，不过收敛速度呈现先快后慢的趋势．

误差分析: (1)同时反演4类(6种)参数共计6×80 = 480个未知数(其中80为离散的网格数)，加剧了逆问题的病态特性，是产生误差的原因之一; (2)从图2的子图(h)可以看出，当前的反演误差还不能满足一般的工程应用需求，可适当增加迭代步数或优化正则化参数，并在反演速度、反演精度之间折中选择; (3)ε∞、σs的反演效果优于ωp1、，其中σs的效果最好，这与激励源的频谱有关［9］; (4)异向峰值型的反演效果优于同向峰值型，产生原因是，对于异向峰值型，更容易获取不同目标的散射场; (5)反演目标的形状、位置信息基本准确，但从反演精度来说，左侧的强散射体相对高于右侧的弱散射体; (6)对比图3的子图(a)、(b)或(c)、(d)均表明，噪声影响也是产生误差的原因之一; (7)对比图3的子图(a)、(c)表明，正则化技术弱化了逆问题的病态特性; (8)对比图3的子图(c)、(d)表明，正则化技术抑制了噪声影响; (9)另外，正问题的求解精度也是影响逆问题病态特性及反演精度的原因之一．

4.2算例2:二维(2-D)问题

FDTD离散化的空间、时间步长分别取为Δx =Δy =Δ =1.0mm、Δt =0.5Δ/c0，T =1500Δt;选取的迭代初值和均匀的背景媒质相同(假设先验知晓)，分别如图5～8的子图(b)所示;考虑噪声且应用正则化条件下，经过60步迭代时，反演结果分别如图5～8的子图(c)所示，其中附加的黑色边框示意散射体S6的真实位置;此时y =6mm处反演参数、真实参数分布之间的对比分别如图5～8的子图(d)所示;反演误差e随迭代步数m的变化关系如图9所示，其中60步迭代的误差为0.1172．

敛散性分析:由图9得到，运用本文的迭代技术处理非光滑型二维问题，算法仍然是收敛的，不过收敛速度较一维光滑情形减缓．

误差分析: (1)从图9与图2的子图(h)对比中发现，60步迭代时，前者的反演误差明显大于后者，一是因为非光滑情形比光滑情形加剧逆问题的病态特性［9］，二是因为二维算例需要同时反演25×25×4 =2500个未知数，明显多于一维情形，显著加重了逆问题的不适性，改进措施是采用其它正则化，如自适应正则化［9］; (2)就参数类型而言，结论类似一维情形，σs的效果最好; (3)就散射强度而言，结论也类似一维情形，更容易检测强散射目标，如S6明显好于S2; (4)就散射体位置(深度)而言，理论上更容易识别浅层目标，而实际上S5不如S3，是因为受到强散射体S6的影响; (5) S1的反演精度好于S4，表明更容易发现大尺寸目标; (6) S2的反演精度不如S1，表明测量视角受限降低了反演精度．

综合一维、二维两个数值算例，不难看出:在噪声环境下、色散(或非色散)背景媒质中，时域全视(或视角受限)反演光滑(或非光滑)、单极(或多极)德鲁色散目标，反演算法收敛;重建的模型参数信息量丰富，重现的形状、位置、尺寸等目标信息基本准确，但反演精度有待提高．

5　结论

德鲁(Drude)模型可精确地描述等离子体、金属等媒质的电色散特性，因此广泛应用在冶金、能源、材料等领域．本文提出一种反演这类媒质模型参数的时域逆散射新技术，克服在时域难于直接反演媒质电参数的困难．在20 dB加性高斯白噪声(AWGN)环境下，通过一维(1-D)、二维(2-D)问题两个数值算例，并通过与无噪声、无正则化情形对比，仿真结果初步证实了该技术的良好性能．下一步研究拟利用电量、磁量的对偶关系，将该项反演技术推广到德鲁磁色散媒质．

参考文献

［1］陈西良，陈欣，朱智勇．不同极性聚合物材料的THz-TDS光谱测量研究［J］．红外与毫米波学报，2013，32(2) :150－159.Chen Xi-liang，Chen Xin，Zhu Zhi-yong．THz-TDS spectra study of polymer materials with different polarity［J］．Journal of Infrared and Millimeter Waves，2013，32(2) : 150－159.(in Chinese)

［2］魏兵，葛德彪，王飞．一种处理色散介质问题的通用时域有限差分方法［J］．物理学报，2008，57(10) : 6290－6297.Wei Bing，Ge De-biao，Wang Fei．A general method for finite difference time domain modeling of wave propagation in frequency-dispersive media［J］．Acta Physica Sinica，2008，57(10) : 6290－6297.(in Chinese)

［3］魏兵，李小勇，王飞，等．一种色散介质FDTD通用吸收边界［J］．物理学报，2009，58(9) : 6174－6178.Wei Bing，Li Xiao-yong，Wang Fei，et al．A finite difference time domain absorbing boundary condition for general frequency-dispersive media［J］．Acta Physica Sinica，2009，58(9) : 6174－6178.(in Chinese)

［4］王飞，魏兵，李林茜．色散介质电磁特性时域有限差分分析的Newmark方法［J］．物理学报，2014，63 (10) : 104101.Wang Fei，Wei Bing，Li Lin-qian．Newmark method for finite-difference time-domain modeling of wave propagation in frequency-dispersive medium［J］．Acta Physica Sinica，2014，63(10) : 104101.(in Chinese)

［5］Winters D W，Bond E J，Vanveen B D，et al．Estimation of the frequency-dependent average dielectric properties of breast tissue using a time-domain inverse scattering technique［J］．IEEE Transactions on Antennas and Propagation，2006，54(11) :3517 －3528.

［6］刘广东，张业荣．二维有耗色散介质的时域逆散射方法［J］．物理学报，2010，59(10) : 6969－6979.Liu Guang-dong，Zhang Ye-rong．Time-domain inverse scattering problem for two-dimensional frequency-dispersive lossy media［J］．Acta Physica Sinica，2010，59(10) :6969－6979.(in Chinese)

［7］Papadopoulos T G，Rekanos I T．Time-domain microwave imaging of inhomogeneous Debye dispersive scatterers［J］．IEEE Transactions on Antennas and Propagation，2012，60 (2) : 1197－1202.

［8］Papadopoulos T G，Rekanos I T．Estimation of the parameters of Lorentz dispersive media using a time-domain inverse scattering technique［J］．IEEE Transactions on Magnetics，2012，48(2) : 219－222.

［9］刘广东，张开银．二维电磁逆散射问题的时域高斯-牛顿反演算法［J］．物理学报，2014，63(3) : 034102.Liu Guang-dong，Zhang Kai-yin．A time-domain Gauss-Newton inversion algorithm for solving two-dimensional electromagnetic inverse scattering problems［J］．Acta Physica Sinica，2014，63(3) : 034102.(in Chinese)

［10］葛德彪，闫玉波．电磁波时域有限差分法(第三版) ［M］．西安:西安电子科技大学出版社，2011.279－280.Ge De-biao，Yan Yu-bo．Finite-Difference Time-Domain Method for Electromagnetic Waves (3rded．)［M］．Xi’an: Xidian University Press，2011.279－280.(in Chinese)

［11］Tikhonov A N，Arsenin V Y．Solutions of ill-posed problems［J］．SIAM Review，1979，21(2) : 266－267.

［12］Sagan H．Introduction to the Calculus of Variations［M］．New York: McGraw-Hill Press，1969.10－110.

［13］Dai Y H．A family of hybrid conjugate gradient methods for unconstrained optimization［J］．Mathematics of Computation，2003，17(2) : 1317－1328.

［14］刘广东，张业荣．一种处理分层有耗色散介质的时域逆散射方法［J］．电子学报，2011，39(12) : 2856－2862.Liu Guang-dong，Zhang Ye-rong．An approach to the time-domain inverse scattering problem for the stratified frequency-dispersive lossy media［J］．Acta Electronica Sinica，2011，39(12) : 2856－2862.(in Chinese)

［15］Roden J A，Gedney S D．Convolutional PML (CPML) : An efficient FDTD implementation of the CFS-PML for arbitrary media［J］．Microwave and Optical Technology Letters，2000，27(5) : 334－339.

刘广东男，1972年生于江苏灌云．现为阜阳师范学院物电学院副教授、工学博士．研究方向为微波医学成像．

E-mail: liu－guang－dong@126.com

葛新同男，1966年生于安徽蒙城．现为阜阳师范学院数学与统计学院讲师、工学博士．研究方向为不确定最优控制及金融数学．

E-mail: gxtong01234@163.com

An Electromagnetic Inverse Scattering Technique in Time Domain for Drude Dispersive Media

LIU Guang-dong1，GE Xin-tong2
(1.School of Physics and Electronic Engineering，Fuyang Normal College，Fuyang，Anhui 236037，China; 2.School of Mathematics and Statistics，Fuyang Normal College，Fuyang，Anhui 236037，China)

Abstract:Drude empirical models are frequently used for description of dispersion characteristics of many media，such as plasmas，and metals．Reconstructed electrical properties by directly using wide-band measured data in time domain，are better than those by application of any single-frequency technique，in amount of information，and resolution of generated images．One of difficulties in time-domain reconstruction of dispersive characteristics is their frequency correlation．In order to overcome this difficulty，an electromagnetic (EM) inverse scattering technique in time domain is proposed，in which four kinds of frequency-independence parameters from a Drude model are estimated simultaneously．Main segments for the technology are: (1) formulating the inverse scattering problem as a constrained minimization problem with a term of regularization; (2) transforming resulting problem into an unconstrained minimization one; (3) deriving a set of closed gradients of its cost functional; (4) solving iteratively resulting forward and backward sub-problems by a finite-difference time-domain (FDTD) method and any conjugate gradient (CG) algorithm，respectively．In one-dimensional (1-D) and two-dimensional (2-D) numerical examples，necessary measurements are replaced by simulated fields based on the FDTD method，and it is assumed that they are corrupted by an additive white Gaussian noise (AWGN)．Numerical results preliminarily confirm performance of the inversion methodology．

Key words:electromagnetic (EM) inverse scattering; Drude dispersive media; regularization; finite-difference timedomain (FDTD) method; conjugate gradient (CG) algorithm

作者简介

基金项目:国家自然科学基金(No．51271059) ;安徽高校省级自然科学研究重点项目(No．KJ2014A193) ;安徽省科技计划项目(No．1501031114) ;阜阳师范学院自然科学研究项目(No．2014FSKJ14)

收稿日期:2015-05-22;修回日期: 2015-08-11;责任编辑:蓝红杰

DOI:电子学报URL: http: / /www．ejournal．org．cn10．3969/j．issn．0372-2112．2016．02．020

中图分类号:TN95; O441; O451

文献标识码:A

文章编号:0372-2112 (2016) 02-0385-07