并行FDTD方法在海面及其上方漂浮目标复合电磁散射中的应用

贾春刚　郭立新, 2　刘伟　尤立志

(1.西安电子科技大学物理与光电工程学院,西安 710071;2.西安电子科技大学 综合业务网理论及

关键技术国家重点实验室,西安 710071;3.中航工业雷达与电子设备研究院

航空电子系统射频综合仿真航空科技重点实验室,无锡 214063)



并行FDTD方法在海面及其上方漂浮目标复合电磁散射中的应用

贾春刚1郭立新1, 2刘伟1尤立志3

(1.西安电子科技大学物理与光电工程学院,西安 710071;2.西安电子科技大学 综合业务网理论及

关键技术国家重点实验室,西安 710071;3.中航工业雷达与电子设备研究院

航空电子系统射频综合仿真航空科技重点实验室,无锡 214063)

摘要利用基于图形处理器(Graphics Processing Unit， GPU)的并行时域有限差分(Finite Difference Time Domain, FDTD)法计算一维粗糙海面及其上方二维漂浮目标的复合电磁散射.采用各向异性完全匹配层(Uniaxial Perfectly Matched Layer， UPML)吸收边界作为截断边界,为了便于并行程序的设计,在整个计算区域使用UPML吸收边界差分公式进行迭代.利用异步通信技术来隐藏主机和设备之间的通信时间,同时使用片上的共享存储器提高读取速度,进一步对程序进行优化,得到很好的加速比,证明了该方法的计算高效性.通过与串行FDTD法以及串行矩量法获得的数值结果进行比较,验证了该并行方法的正确性,进而研究了海面上方类似舰船漂浮目标的电磁散射特性,讨论了入射角、海面风速以及目标吃水深度对双站散射系数的影响.

关键词复合散射;海面;漂浮目标;FDTD;GPU

DOI10.13443/j.cjors.2015010501

Application of parallel FDTD to EM scattering from a target floating on sea surface

JIA Chungang1GUO Lixin1,2LIU Wei1YOU Lizhi3

(1.SchoolofPhysicsandOptoelectronicEngineering,XidianUniversity,Xi’an710071,China;2.StateKeyLaboratoryofIntegratedServicesNetworks,XidianUniversity,Xi’an710071,China;3.AviationKeyLaboratoryofScienceandTechnologyonAISSSRadarandAvionicsInstituteofAVIC,Wuxi214063,China)

Abstract In this paper, GPU-based FDTD algorithm is applied to study the electromagnetic(EM) scattering from two-dimensional(2-D) target floating on one-dimensional(1-D) rough sea surface with Pierson-Moskowitz(PM) spectrum. The FDTD lattices are truncated by uniaxial perfectly matched layer (UPML), and the finite-difference equations are employed in the whole computation domain for the parallelization convenient to carry out. Also, the parallelism design is limited to the iteration of the near field that is extremely time consuming. To improve the performance, asynchronous transfers is implemented to mask the memory transfers time and the shared memory is used to achieve high memory bandwidth. Using compute unified device architecture(CUDA) technology, significant speedup ratios are achieved, which demonstrates the efficiency of GPU accelerated the FDTD method. The validation of our method is verified by comparing the numerical results with these obtained by sequential FDTD executing on CPU as well as method of moments (MOM), which shows favorable agreements. Furthermore, our parallel implementation is employed to study the influences of the incident angle, the wind speed, the depth of the target on the EM scattering from the target and a sea surface composite model.

Keywords composite scattering; sea surface; floating target; FDTD; GPU

引言

随着雷达遥感和目标检测技术的飞速发展,因其在海洋遥感、目标识别和雷达探测等军事以及民用领域具有重要的应用价值[1-2].海面与目标复合电磁散射的研究越来越得到国内外学者的广泛关注.目前为止,针对这种随机问题的解决方法总结起来主要分为两大类:高频近似方法和数值算法.其中,以高频方法为例:基尔霍夫近似[3]、小斜率近似[4]、微扰法[5]以及双尺度[6]等,由于上述方法只限于单次散射,无法考虑多次散射与多径传播、边缘衍射、相位干涉产生的Bragg谐振散射等复杂过程,用来解决单纯粗糙面散射问题.对于粗糙面与目标复合散射,目标与粗糙面复杂的相互作用,因此数值方法广泛地用于计算上述复合电磁散射模型,如矩量法[7]、快速多极子[8]、时域有限差分法(Finite Diffe-rence Time Domain, FDTD)方法[9]等.数值方法由于直接对Maxwell微分或积分方程进行数值求解,因而适用范围广,并且计算精确度高,但受计算机内存和时间的限制,计算效率不高.因此,针对数值方法的加速算法相继被提出,如基于信息传递接口(Message Passing Interface，MPI) 的并行矩量法[10]、并行FDTD 法[11],有效地减少了计算时间.

然而,上述并行方法的加速比受限于计算机CPU核数,对计算机的硬件要求很高.近些年,基于GPU的统一计算设备架构(Compute Unified Device Architecture, CUDA)技术广泛地用于电磁场的数值模拟中[12-13].与并行的CPU技术相比,GPU由于其强大的计算能力能够轻松地实现很高的加速比.FDTD方法在求解粗糙面与目标复合散射问题时,有其自身的优点,由于采用统一剖分建模,电磁参数已经反映在每一个网格的电磁场计算中,因此它很容易处理复杂介质电磁散射模型.并且电场(磁场)分量的迭代公式只与它本身上一时刻的场值和周围网格磁场(电场)上半个时间步的值有关,而与计算区域内其它场量没有直接的关系,因此它具有天然的并行能力.基于GPU 的并行FDTD方法虽然已经应用于很多的电磁场数值模拟中,但在粗糙面电磁散射特性的研究中尚未涉及到.因此,本文首次提出利用FDTD方法结合GPU并行技术解决粗糙面与目标复合模型电磁散射问题.与传统的FDTD方法相比,该方法在处理电大尺寸粗糙面与目标复合散射问题时具有明显优势.

本文采用FDTD方法结合基于GPU的并行技术讨论了一维海面与上方漂浮类舰船目标的复合电磁散射特性.第二节介绍了海面的生成以及复合散射的FDTD方法计算模型,并且详尽地描述了基于GPU的FDTD方法加速算法在解决粗糙面与目标复合散射时的并行以及优化方案;第三节从计算准确性和计算效率两个方面验证了本文基于GPU的FDTD方法解决粗糙面与目标复合散射问题的有效性,在此基础上,进一步讨论了入射角、海面风速以及目标吃水深度对海面与上方漂浮类舰船目标复合电磁散射特性的影响;最后,给出了相应的结论以及今后工作的研究重点.

1理论分析

1.1复合散射模型

图1以TM入射波为例,给出了一维海面及其与上方二维漂浮类似舰船目标复合的电磁散射模型,入射波入射矢量ki,与y轴正方向夹角为入射角θi,ks为散射波矢量,与y轴正方向夹角为散射角θs,h为舰船目标的吃水深度,L为海面的模拟长度y=f(x)为海面轮廓,本文采用Pierson-Moskowitz (PM)谱来描述海面轮廓的起伏,功率密度函数为[14]

(1)

式中： 空间波数K=2π/L;U19.5为海面高度19.5 m处的风速；α=8.10×10-3以及β=0.74分别为经验常数; 重力加速度常数gc=9.81 m/s2.图1复合散射模型中,当入射波为平面波时,总场边界需设置均匀的等效面电磁流,而总场边界的左右两端被吸收边界截断,导致等效面电磁流在左右两端处有突变,为了克服这一缺点,本文引入了高斯窗函数[15]

(2)

式中:xcen和ycen为连接边界处的中心点坐标;T为常数,决定了高斯窗函数的宽度,当边缘值为中心点处的10-3时,满足关系cos(θi)/T=2.6/ρm,其中ρm为中心点到边缘处的最小距离,对于一维粗糙面而言,与y方向上的的变化没有关系,即y=ycen.

图1　一维海面与舰船目标复合电磁散射模型

由于海面与目标的复合电磁散射属于半空间问题,图 1中,只在计算区域的上半部分加入了虚拟的外推边界以及连接边界.为了在有限的区域内模拟无限的空间电磁散射问题,在计算区域最外层的截断边界处加入吸收边界,本文使用了各向异性完全匹配层吸收边界[16],为了便于并行化设计,在整个计算区域使用吸收边界的差分公式进行迭代;连接边界[16]将计算区域划分成总场区和入射场区,并应用惠更斯原理,在此边界引入入射波;当近场迭代达到稳定之后,在外推边界处的近场做近远场变换实现远场的外推[16],得到远区散射场,进而获得对粗糙面长度归一化的双站散射系数:

(3)

1.2FDTD方法计算海面与目标复合散射的GPU并行方案

该部分介绍了并行平台以及CUDA编程模型.详尽地阐述了基于GPU的FDTD方法计算粗糙面电磁散射的并行策略,在此基础上,利用异步传输技术以及共享存储器对其进一步优化,提高并行效率.

1.2.1CUDA编程模型

英伟达推出了CUDA编程模型的GPU显卡之后,由于其相对简单的并行化设计,凭借其强大的计算能力在各个领域到了广泛的推广和应用.串行程序在主机(CPU)上执行,具有高度并行化的数据则在设备上(GPU)进行处理.并行程序首先在主机上串行执行,包括CPU内存和GPU显存的分配、初始化以及内存和显存的释放;在设备上由成千上万线程并行执行的函数称为核函数,不同核函数在主机上串行执行.在CUDA编程模型中,主机和设备上的存储器在物理上是分离的.

1.2.2并行化设计

图2所示为并行FDTD方法计算粗糙面电磁散射的并行化设计的流程图.首先如前文所述,利用模特卡洛方法模拟PM海面的轮廓建立复合散射模型,然后CPU根据计算问题的未知量个数分配内存和显存大小.FDTD方法计算时,在近场的迭代很耗时并且具有很好的并行特性,所以并行化处理在此处进行.由于整个网格线程间需要进行数据交换,所以需要重新调用一个核函数进行同步,但全局存储器的读写速度很慢,所以应尽量使用少的核函数,本算法使用5个核函数来实现近场迭代功能,包括IncidentHKernel()(入射波磁场分量的更新)、IncidentEKernel()(入射波电场分量的更新)、ConnectionKernel()(连接边界处引入入射波)、eKernel()(电场分量的更新)、hKernel()(磁场分量的更新).当GPU完成近场迭代时,远区散射场很容易由CPU平台得到.

图2　并行FDTD方法计算海面与目标复合散射流程图

1.2.3并行优化方案

1) 共享存储器

以TM入射波为例,在Ez(i,j)分量进行迭代时,需要前一时刻的Hy(i,j)和Hy(i-1,j)以及Hx(i,j)和Hx(i,j-1),而其相邻的Ez(i-1,j)分量进行迭代时,需要前一时刻的Hy(i-1,j)和Hy(i-2,j)以及Hx(i-1,j)和Hx(i-1,j-1),所以可以看出在当前Ez分量进行迭代时,需要从全局存储器中读取两次的前一时刻的Hx和Hy分量,由于在对全局存储器读写时的访问延迟很高,所以极大地影响了并行效率.因此,利用片上的共享存储器提高读取速度,只需要从全局存储器中读一次数据到共享存储器中,然后数据就可以在共享存储器中进行重复的读取,大大降低了线程间的通信时间.当磁场分量Hx和Hy进行迭代时,不仅当前线程块的电场的值Ez要从全局存储器拷贝到共享存储器里,而且其相邻的右边线程块最左边一列的值和相邻的下边线程块最上面一行的值也要传到共享存储器;当电场分量进行迭代时,当前线程块的磁场分量Hx的值,以及相邻下边线程块的最上一行的值都要传到共享存储器.磁场分量Hy当前线程块以及相邻右边线程块最左边一列的值都要传到共享存储器中.

2) 异步传输

主机与设备之间的数据交换非常耗时,尤其在计算电大尺寸问题时,主机与设备之间的数据交换量巨大,极大地影响了并行程序的计算效率.因此,本文利用CUDA模型中的异步功能来隐藏主机与设备之间的数据通信,即设备端的计算与主机与设备之间的数据交换可同时进行.异步传输是主机上的pinned memory 通过不同流执行cudaMemcpyAsync()函数实现主机与设备之间的数据通信.同时多个流并发执行,当前流的数据传输结束之后,核函数立刻启动,同时另一个流开始执行数据传输,通过数据传输与计算重叠的执行,能够很好地隐藏主机与设备之间的数据通信时间,进一步提高并行效率.由于在后面主机端进行近远场外推时需要电场分量Ez以及磁场分量Hx的值,需要在时间迭代稳定后将设备端的值传回主机,为了隐藏数据通信时间,在此用了6个异步流并发执行来实现异步传输功能.在拷贝时offset不仅包括了当前计算网格的数据长度,也包括了周围相邻网格(offset_boundary),这是因为当前网格的计算与周围相邻网格的值有关,如果相邻网格没有包括,这会导致当前一部分网格没有进行迭代,这与全局存储器和共享存储器之间的数据传输方式是一样的.

2数值仿真结果

本节将利用并行FDTD方法对海面与上方目标的复合电磁散射特性进行详细的讨论分析.为了保证FDTD的准确性和稳定性,空间步和时间步为Δx=Δy=Δ=λ/20与Δt=0.5×Δ/c,λ为入射波波长,c为真空中的光速.UPML层的厚度为10Δ.计算平台如表1所示.

表1　计算平台参数

为了验证本文方法的正确性,以无限长导体圆柱与PM海面复合散射模型为例,与串行FDTD方法以及矩量法(Method of Moment, MoM)进行对比.图3分别采用并行FDTD方法、串行FDTD方法和MoM计算TM 波和TE波入射下该模型的散射系数.其中,入射频率f=1.0 GHz,入射角θi=40°,海面的长度L=4 096Δ,与入射频率对应的海水介电常数εr=(72.0,84.0)(德拜模型)[17],海面风速U19.5=3.0 m/s,圆柱的高度H和半径r分别为4λ和λ.由图3可以看出,三种方法在大部分散射角范围内都有很好的一致性,其中并行FDTD方法与串行FDTD方法吻合得非常好,两者与MoM方法具有一定的差异,这是由FDTD方法和MoM算法本身导致的.

表2　不同计算网格的加速比

(a) TM

(b) TE图3　海面与其上方圆柱复合散射并行FDTD、串行FDTD以及MoM的比较

为了验证本文方法的效率,表2分别比较了不同网格数下的使用异步传输和没有使用异步传输的加速比.从表中可以看出在一定网格数的范围内,两种传输模式的加速比随着网格数的增加而逐渐增加,网格数增大到262 144Δ时,同步传输模式加速比减小,而异步传输模式下的加速比基本保持不变,并且在相同的计算网格的条件,异步传输模式的加速比比同步传输模式的加速比要大.随着计算网格数的增加,GPU的计算单元逐渐被充分利用,所以加速比随之增大,但当增大到262 144Δ时,由于数据量很大,同步传输模式下主机与设备之间数据交换时间的影响变得明显,所以加速比较小; 而在异步传输模式下,主机与设备的通信时间被很好地隐藏,所以加速比基本不变,同时相同网格下的异步传输模式的加速比比同步模式下的加速比大,也说明了这个问题.从表中也可以看出,本文算法较串行方法具有很高的计算效率,特别是在计算电大尺寸问题时优势明显.

在计算准确与高效的基础上,接下来,利用本文方法研究不同入射角下海面与上方漂浮类似舰船目标复合模型的电磁散射特性,如图4所示,粗糙面长度16 384Δ,海面风速U19.5=3.0 m/s,入射频率f=2 GHz,对应海水介电常数εr=(71.6,49.1),入射角θi分别为30°、45°、60°.吃水深度h=40Δ,从图4中可以看出,对于TM和TE两种入射情况,随着入射角的增大,镜向处散射减弱,并出现两个峰值,非镜向方向的散射增强,并且随着入射角的增大,非镜向方向有明显的双峰,特别当入射角θi=60°时,对于两种极化,θs=-60°以及θs=-10°附近峰值明显.

(a) TM

(b) TE图4　不同入射角的双站散射系数

图5给出了不同风速复合模型双站散射系数随散射角变化的结果.入射角θi=30°,风速U19.5分别为3 m/s、6 m/s、9 m/s,其他参数与图4中的相同.对两种入射而言,随着风速的增加,镜向附近的散射减弱.这是由于随着风速的增加,海面变得粗糙,均方根斜率增加,进而导致镜向方向散射减小.

图6研究了一维PM海面与上方漂浮导体舰船目标的复合散射特性随着舰船吃水深度变化的关系.其中,入射角θi=30°,舰船的吃水深度h分别为50Δ(0.375 m)、150Δ(1.125 m)、250Δ(1.875 m),其他参数与图4中相同.从图6中可以看出,对于两种入射情况来说,随着吃水深度的增加,非镜向方向散射减小明显,特别是在θs=-90°到θs=-30°之间,镜向附近变化不是很明显.这是因为随着舰船吃水深度的增加,目标在海面上方的部分减小,舰船与海面之间的耦合作用在减弱,进而导致非镜向方向的散射减弱.

(a) TM

(b) TE图5　不同风速的双站散射系数

(a) TM

(b) TE图6　不同吃水深度的双站散射系数

3结论

本文将基于GPU的并行时域有限差分法引入到一维PM海面与上方漂浮舰船目标的复合电磁散射问题的计算中.从计算准确性与效率两方面,验证了保证串行方法高计算精度的前提下,与传统串行方法以及MPI-based FDTD相比,在计算效率方面本文方法具有明显优势.在此基础上,利用该方法进一步分析了海面与上方漂浮目标复合电磁散射特性.发现随着入射角的增大,镜向附近的散射减弱,非镜向方向散射逐渐增强,甚至出现几个峰值,这是由于舰船目标的影响较大;随着风速的增加,镜向附近双站散射系数减小,非镜向方向的双站散射系数增大,这是由于风速的改变导致了海面粗糙面的变化;最后分析了舰船吃水深度对双站散射系数的影响,数值结果表明,随着吃水深度的增加,舰船与海面耦合作用减弱,导致非镜向方向散射减弱.由于本文方法具有很高的计算效率,所以在计算电大尺寸目标与大尺度粗糙面复合散射问题时具有明显优势,所以利用GPU-based FDTD方法解决更为复杂耗时的三维目标与二维粗糙面的复合电磁散射问题将是我们今后研究工作的重点.

参考文献

[1] VORONOVICH A G, ZAVOROTNY V U. Full-polarization modeling of monostatic and bistatic radar scattering from a rough sea surface[J]. IEEE transactions on antennas and propagation, 2013, 62(3): 1362-1371.

[2] LI Y, ZHANG Y Z, CHEN J, et al. Improved compact polarimetric SAR quad-pol reconstruction algorithm for oil spill detection[J]. IEEE geoscience and remote sensing letters, 2014, 11(6): 1139-1142.

[3] THORSOS E I. The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum[J]. The journal of the acoustical society of America, 1988, 83(1): 78-92.

[4] WANG Y Q, BROSCHAT S L. A systematic study of the lowest order small slope approximation for a Pierson-Moskowitz spectrum[J]. IEEE geoscience and remote sensing letters, 2011, 8(1): 158-162.

[5] THORSOS E I, JACKSON D R. The validity of the perturbation approximation for rough surface scattering using a Gaussian roughness spectrum[J]. The journal of the acoustical society of America, 1989, 86(1): 261-277.

[6] DURDEN S L, VESECKY J F. A numerical study of the separation wavenumber in the two-scale scattering approximation[J]. IEEE transactions on geoscience and remote sensing, 1990, 28(2): 271-272.

[7] GUAN B, ZHANG J F, ZHOU X Y, et al. Electromagnetic scattering from objects above a rough surface using the method of moments with half-space Green’s function[J]. IEEE transactions on geoscience and remote sensing, 2009, 47(10): 3399-3405.

[8] GUO L X, WANG A Q, CHAI C. Parallel fast multiple method for electromagnetic scattering from one-dimensional large-scale two-layered rough surfaces for large angles of incidence[J]. IET microwaves antennas and propagation, 2011, 5(15): 1813-1821.

[9] LI J, GUO L X, JIAO Y C, et al. Investigation on wide-band scattering of a 2-D target above 1-D randomly rough surface by FDTD method[J]. Optics express, 2011, 19(2): 1091-1100.

[10]张连波, 郭立新, 苟雪银, 等. 三层粗糙面电磁散射的矩量法研究[J]. 西安电子科技大学学报, 2013, 40(6): 147-154.

ZHANG L B, GUO L X, GOU X Y, et al. Method of moment investigation on electromagnetic scattering from the three-layered rough interfaces[J]. Journal of Xidian University(natural science), 2013, 40(6): 147-154.(in Chinese)

[11]阎亚丽, 傅光, 龚书喜, 等.基于并行FDTD方法分析表面等离子波导的特性[J]. 电波科学学报, 2015,30(4): 668-672.

YAN Y L, FU G, GONG S X, et al. Analysis of a surface plasmonic waveguide using parallel finite difference time domain method[J]. Chinese journal of radio science, 2015, 30(4): 668-672. (in Chinese)

[12]SU X, WU J J, HUANG B, et al. GPU-accelerated computation for electromagnetic scattering of a double-layer vegetation model[J]. IEEE journal of selected topics in applied earth observations in remote sensing, 2013, 6(4): 1799-1806.

[13]JIANG W Q, ZHANG M, WEI P B, et al. CUDA-based SSA method in application to calculating EM scattering from large two-dimensional rough surface[J]. IEEE journal of selected topics in applied earth observations in remote sensing, 2014, 7(4): 1372-1382.

[14]THORSOS E I. Acoustic scattering from a “Pierson-Moskowitz” sea surface[J]. The journal of the acoustical society of America, 1990, 88(1): 335-349.

[15]FUNG A K, SHAH M R, TJUATJA S. Numerical solution of scattering from three-dimensional randomly rough surfaces[J]. IEEE transactions on geoscience and remote sensing, 1994, 32(5): 986-994.

[16]TAFLOVE A, HAGNESS S C. Computational electrodynamics: the finite-difference time-domain method[M]. Norwood: Academic, 2005: 297-365.

[17]KLEIN L A, SWIFT C T. An improved model for the dielectric constant of sea water at microwave frequencies[J]. IEEE transactions on antennas and propagation, 1977, 25(1): 104-111.

贾春刚(1986-),男,吉林人,博士研究生,研究方向为粗糙面及其与目标的并行加速计算.

郭立新(1968-),男,陕西人,西安电子科技大学物理与光电工程学院教授,研究方向为电磁/光波在复杂介质中的传播和散射.

刘伟(1979-),男,陕西人,西安电子科技大学副教授,主要研究方向为电磁散射、微波遥感与参数反演.

作者简介

中图分类号TN011

文献标志码A

文章编号1005-0388(2016)01-0116-07

收稿日期：2015-01-05

贾春刚, 郭立新, 刘伟, 等. 并行FDTD方法在海面及其上方漂浮目标复合电磁散射中的应用[J]. 电波科学学报，2016,31(1):116-122. DOI: 10.13443/j.cjors.2015010501

JIA C G, GUO L X, LIU W, et al. Application of parallel FDTD to EM scattering from a target floating on sea surface[J]. Chinese journal of radio science，2016,31(1):116-122. (in Chinese). DOI: 10.13443/j.cjors.2015010501

资助项目: 国家杰出青年科学基金(61225002); 航空科学基金与航空电子系统射频综合仿真航空科技重点实验室联合资助(20132081015)

联系人： 贾春刚 E-mail: cgjia@stu.xidian.edu.cn