复杂电磁环境三维仿真及可视化系统

张 霖，沈月伟，吴迎年，2，穆 兰，王海波，吕彦东

(1.北京航空航天大学电气信息与自动化学院，北京 100191;2.北京信息科技大学自动化学院，北京 100192)

0 引言

近几十年来，随着无线通信和现代战争技术的发展，预测无线信号的传播范围、雷达的探测覆盖区域等变得越来越重要。准确的电磁特性数据是复杂电磁环境(EME，electromagnetic environment)预测过程中不可或缺的组成部分，而EME仿真作为一种简单、灵活、高效的技术，已成为获取电磁特性数据的最常用的方法之一。EME仿真过程中需要考虑的因素包括电磁干扰，电磁脉冲，电磁辐射对人员、军械和挥发性材料危害，以及雷电和沉积静电等自然现象的综合［1］。

复杂EME(CEME)是一种动态多变的复杂巨系统，尤其是在几十、几百公里这样大区域内，电磁传播将受到复杂地形地貌、复杂气象条件、各种形状材质建筑物、动态移动目标等各种因素的影响，仿真的快速性和准确性受到了严峻的挑战。产生CEME的主要因素包含电波源辐射传播，以及辐射传播因素中的地形因素等。对于较为复杂的电磁计算问题，目前人们已提出了多种确定性方法，如时域有限差分算法(FDTD，finite difference time domain)［2］、矩量法(MoM，method of moments)［3］、射线跟踪(RT，ray tracing)［4］等，这些方法对于计算大范围的CEME仿真问题都不能获得较好的效果。

抛物方程(PE，parabolic equation)是近年来得到广泛关注的电磁计算方法之一［5，6］，相比其他诸多电磁计算方法，PE方法在远距离大范围的对流层电波传播特性预测和电大尺寸目标的雷达散射截面(RCS)计算上显示了其他方法所不具备的独特性能［7］:PE本身体现了电波折射与绕射效应，不用再计算电波传播过程中各种传播机制的区域与繁琐的判据;对复杂地形与不均匀大气折射特性的处理简洁而精确，计算大区域场强衰减特性时，与其他方法相比计算精度高并且速度快。从目前国内外相关研究来看，基于标量的二维抛物方程(2DPE)方法的相关研究与应用已日趋成熟［6］，三维抛物方程(3DPE)方法的理论研究和应用在电磁场计算精度要求较高的情况下需求迫切。作为大范围CEME三维仿真的有效方法之一，3DPE法的计算时空复杂度仍然很大，很难获得满意的效果，P3DTEMPER［8］和虚拟 3D 方法［9］能大大提高计算效率，但是其计算范围还是比较小，也没有采取有效的加速方法。在2DPE的基础上，提出了采用并行度高的准三维抛物方程(quasi-3DPE)方法来解决3D空间电磁分布的问题，使其计算速度显著加速，计算范围得到扩展。

采用PE方法实现了CEME仿真，第1部分介绍了2DPE法的推导、初始场、边界条件及不规则地形的处理方法;第2部分介绍了准三维方法的构建方法，包括场源建模、地理环境构建，对其精度进行了分析，并介绍了并行仿真方法;第3部分介绍了电磁环境三维可视化系统的设计与实现;第4部分进行了总结。

1 二维抛物方程方法研究

2DPE是进行EME仿真的一种基本方法，quasi-3DPE法就是基于这种方法进行的。下面将对于2DPE的基本推导和离散过程进行分析。

1.1 基本推导过程

二维标量波动方程表示为［6］

式中，k代表真空中传播的波数;ψ代表电场或磁场的一个任意分量;n代表媒质的折射系数。假设该选择的距离范围内n保持不变，并设定恰当的边界条件来处理该范围内的截断边界。一般的方法是把空气层和地表面的分界面作为下部边界条件，而将上部边界扩展到无限远处。电波沿着一定的角度向外传播，选取x轴正方向为传播方向，竖直方向用z来表示，如图1所示。

图1 二维面内的电波传播方向示意图

经过一系列的近似，只考虑电波的前向传播问题，可得标准抛物方程(SPE)［6］为

1.2 离散过程

采用有限差分(FD)方法对SPE进行离散，这里采用的是隐式的C－N格式［6］，可以得到

将其转换成矩阵形式［6］

其中系数矩阵Am是一个三对角的形式［6］

1.3 二维电磁传播衰减值变化图

对于二维EME传播问题，采用2DPE方法进行仿真，计算了在不同的方位角上的传播衰减值，所得的效果如图2所示。

图2 不同方位角上的传播衰减值变化图

2 准三维抛物方程方法

对于大范围的EME仿真，3DPE方法虽然已取得了一定的研究成果［10，11］，但目前仍存在着处理边界条件复杂，计算量大，仿真时间过长等问题。而2DPE法却能高效地预测仿真在对流层大气2D空间面内的电磁传播问题［12］，由此将3D电磁传播空间以场源为中心，剖分成一个个的2D面，在每个面内再用2DPE法来并行计算电磁传播，得到组合的三维仿真可视效果，称此方法为准三维方法，如图3所示。这种方法的复杂度降低，效率高，且易于并行。

图3 准三维方法仿真示意图

2.1 场源建模

PE的一般求解过程是步进迭代过程，对于SPE，每个剖面内都需要解决式(2)的初始值问题。无论标量3DPE还是2DPE法，获得初始场的方法都是依据天线的远场方向图通过近场/远场变换原理得到。首先考虑如何从3D的天线方向图等效出各个辐射角度上的2D辐射初始场。天线方向图是角度θ和Φ的2D函数，从标量3DPE出发，根据天线的近/远场变换原理，自由空间中，电磁场量的一个场分量与天线方向图的关系为

式中，B(θ，φ)为天线的辐射方向图;k为波数;u(0，y，z)为 3D 初始场。

因为上述剖分依据的是自由空间中场分量的远场方向图，其等效性仅限于标量3DPE和2DPE法之间，即在实际传播环境中应用时，quasi-3DPE法忽略了极化效应和各分量的耦合。实现上，对傅里叶变换式可以先按采样定理进行离散化，再采用快速傅里叶变换进行数值求解，最终获得离散的初始场。

2.2 地理环境建模

地理环境建模主要应该考虑到不规则的地形、多样的地貌因素对电磁传播的影响，尤其是不规则的地形起伏对于电波的反射、绕射等传播机制有重要影响。对于quasi-3DPE法，按照剖分方位角从2D数字地图描述的地形表面中采用合适的算法抽取出1D的地形起伏线，如图4(b)所示，每条地形起伏线即是2D仿真剖面中的下边界。可以按需求选择阻抗边界条件或良导体边界条件，对于上边界处理可以采用PML吸收边界条件或加窗函数或增大计算域高度减少对所关注仿真域的数值反射等三种方式。关于地形起伏数据，采用数据存储格式相对简单的数字高程模型(DEM，digital elevation model)。为了在PE法仿真中能随意抽取地图中的地形剖面，在对USGS DEM数据存储格式的分析基础上，采用内插方法编制了数据抽取程序，基于实际地理环境中地形起伏变化的连续性，如图4所示。

图4 双线性插值获取地形剖面示意图

接着采用分块内插法中简单而又有效的双线性插值算法来计算数字地图中非节点处的高程值。

2.3 分析与比较

2.3.1 误差分析

在quasi-3DPE法中，每一个剖面上的计算方法都是2DPE法，关于2DPE法的计算精度已有诸多文献论述［6］，2DPE 法的应用已较成熟可靠，2DPE 法的广泛应用也正是quasi-3DPE法的一定精度保证。quasi-3DPE可以看作是对于3DPE法的2D组合近似，其中忽略了剖面间的横向电波绕射和散射，且没有考虑去极化效应，故其精度不如3DPE法高。P3DTEMPER［8］已经与3DPE进行了对比，在此基础上分析地形的横向变化对准三维方法仿真的影响，这里在理论上先定性地给出影响准三维方法精度的因素。

首先分析不规则地形下的3DPE法。假设传播环境为无源介质，时谐因子为exp(－iwt)，引入了位函数来表示任意电场或磁场，这是因为电磁学中任意电场或磁场都可以表示成一个横电场和横磁场的叠加。假设ψ=ψe+ψm为电位或磁位的任意标量场值，那么它自动满足标量波动方程

用 eik0xψ(x，y，z)代换上式中的 ψ，然后按照由波动方程推到PE的过程，可以得到用电磁位函数表示的3DPE形式

式中，

注意这里的抛物型方程本质上不同于标量形式的3DPE，因为ψ是电磁位函数，所以式(9)本质上表达的是三维矢量抛物方程形式。然后经过一系列的近似和转换，可得到

从式(10)与式(2)可知若不考虑大气因素，影响电磁在横向上传播的因素主要是地形的高度变化率，即可以得知准三维方法的误差来源除大气之外的主要因素是不规则地形的高度变化率，即在地形高度变化剧烈处误差当增大和积累。

2.3.2 与其他方法的比较

P3DTEMPER方法［8］是基于二维的TEMPER模型而提出的一种伪三维方法，它大大提高了计算效率，然而它忽略了散射和绕射的作用，通过与全三维方法的对比，它的计算结果是有效的;虚拟3D策略［9］是基于二维APM来构造三维的雷达作用范围的一种方法，它能及时准确形象地展示复杂环境影响下雷达的三维作用范围，允许交互调整雷达和目标参数，从而为用户决策和规划提供支持。以上两种方法与quasi-3DPE方法类似，都是用二维模型来模拟全三维的方法，但是相比较而言还存在着以下几点区别:①这两种方法采用的二维模型与本文不同;②这两种方法都未能实现并行加速计算;③应用的领域也是不同的;④计算的传播范围上，这两种方法都不如本文的范围广。

David Mendez等［13］采用“随机微分方程”的方法研究传播特性，它把导波因子建模为随机变量，在计算椭圆极化波传播时能获得更好的结果。与quasi-3DPE方法相比，这种方法在计算时需要经过几次复杂的数学变换，导致其计算复杂度非常大，相应的会使计算效率大幅下降，尤其在计算复杂的大范围电磁传播问题时，很难达到满意的效果。与上述方法相比，提出的quasi-3DPE方法是实用且有效的，能对于复杂问题获得较好的结果。

2.3.3 与商业软件之比较

Wireless InSite(WI)是REMCOM公司开发的一款对CEME进行仿真预测分析的软件，它基于UTD/GTD理论，采用射线跟踪(RT，ray tracing)方法建立传播模型，可用于对城市、郊区、室内等规则区域，山脉、植被区等非规则区域，机场、大型船舰等复杂平台的EME预测。关于准三维方法的计算精度和速度，这里与WI软件进行初步比较，两方法中场源都为高斯型，辐射功率都设为0 dbm，最大传播仰角10°，地形数据都采用同一数据文件，发射天线高度都为150.0 m，频率都为1.0 GHz。地形表面都设置为PEC，大气都设为真空。其他参数和计算时间见表1。

表1 二维剖面内抛物方程法和射线跟踪法计算时间比较

用基于准三维方法开发的软件和WI软件得到路径损耗的可视化效果，如图5所示。

图5 两种方法的路径损耗值宏观显示示意图

其中暖色调代表电磁能量高，冷色调处电磁能量低，即前者路径损耗小，后者路径损耗值大。从准三维方法结果中抽取场源辐射最强角度处的2D剖面数据，将它和WI对应同一2D面上的结果进行比较。下面分别抽取了距离场源10 km和15 km处0～600 m高度的路径损耗值［14］，对比曲线如图6所示。这里所求路径损耗值的计算公式依据文献［6］

式中，Lp为考虑了自由空间损耗的路径损值;λ为波长。从上述路径损耗值的比较来看，quasi-3DPE与RT法的计算结果吻合。但从可视化效果上来看，准三维法中组合出的三维显示效果不如WI，要解决此问题需要增多计算剖面，组合出更细致的显示效果，增多剖面势必会增大计算量，在单机上会遭遇计算困难。但由于准三维法的并行度高，只要高性能机的计算能力可扩展，增多的计算量将不会称为瓶颈。

图6 距离发射机不同距离的损耗值比较

2.4 并行仿真方法

准三维方法中将三维空间剖分成一个个单独的二维剖面，并利用2DPE法进行剖面内的电磁传播仿真。假设应用准三维法时共需要计算n个剖面，每个剖面的的计算时间为Ts(1)，串行计算n个剖面，理论上需要时间nTs(1)。因为从准三维法仿真模型的结构上来看，各计算过程中各剖面之间没有耦合，假设并行机的计算能力无限，每个剖面都有一个处理器和相应独占的内存来计算，那么理论上准三维剖面中的每个剖面都可以同时在Ts(1)时间内完成求解计算，即并行计算时间为Ts(1)，所以其理想情况下的并行加速比为n，相应地，其理想情况下的并行效率为100%。

实际上，在并行机的计算资源有限的情况下，不可能每个剖面都能够获得足够计算资源，所以计算中就产生了串行部分，加速比下降。另外在每个剖面进行计算时，其所需的地形数据可能存在同一个存储文件或数据库中，需要串行的读取，或者在计算过程中需要向相同的文件或数据库存储计算结果，需要串行的写入。所有这些产生串行计算部分的因素都会造成计算时间的延长和并行效率的降低，实际情况下不可能出现百分之百的并行效率。

因为准三维方法计算时各剖面间无耦合，与三维抛物方程并行求解过程相较，并行度高。在并行机上计算无需通信，也避免了复杂并行算法的设计和实现。

简单的并行求解实现方式可以采用在MPI并行程序实现时将每个剖面的计算任务分配到一个进程。并行编程模型简单，这也是准三维法在并行机上容易实现的优点。若为了实现更深层次的并行化，可以采用多级的并行结构，如图7所示。

图7 多级并行计算结构示意图

任务分解层次上，按角度进行剖分，作为单独的任务进行分发并行。对于每个剖面内比较耗时的部分包括不规则地形剖面的抽取，以及求解模型的线性方程组计算和结果的存储。对于这些可以采用多进程或多线程进行并行化。其中对于求解2DPE模型中三对角方程组的矩阵计算，曾在曙光TC4000L集群上采用并行分裂法对十万阶的三对角线性方程组进行过求解效率的测试，实验表明［5］采用MPI多进程求解的并行效率不高，这主要归因于MPI进程级的通信开销较大。对于准三维方法的多级化并行求解，将来可以采用粗粒度上多进程与细粒度上多线程相结合的方法来获得更高效率。

3 电磁环境三维可视化系统设计

前面介绍了CEME计算及其并行策略，为了更直观的显示出计算所得的数据，将设计并实现一个空间电磁环境可视化原型系统［15，16］，并且通过应用实例来验证此原型系统。

3.1 系统功能设计

空间电磁环境可视化系统通过设置电磁环境参数，计算求解电磁模型，将地形数据和电磁数据经过数据建模生成统一的数据模型，再通过解析映射后生成可视化结构，对其绘制渲染实现空间电磁环境的可视化与人机交互。

系统功能如图8所示，空间电磁环境可视化系统主要实现电磁环境参数设置、电磁模型求解、数据统一建模、解析映射、可视化显示和人机交互等功能。

图8 系统功能设计图

3.2 系统架构

电磁环境可视化系统采用三层架构设计，从下至上分别为:数据访问层、业务逻辑层、表示层。由图8可看到，空间电磁环境可视化系统首先设置电磁环境参数，从而计算出电磁数据。对地形数据与电磁数据进行读入映射转换，生成统一的空间电磁环境数据模型，然后对其经过解析与可视化映射生成空间电磁环境可视化结构，并对可视化结构根据设定的位置坐标信息、光源及材质、顶点颜色属性等，把三维电磁环境的所有顶点经过世界变换、观察变换、投影变换为计算机显示器窗口上的顶点，再根据顶点颜色填充整个平面，控制计算机显示器屏幕每点的颜色，完成三维电磁环境的可视化显示从而进行绘制与渲染，以电磁场强分布态势、电磁等值线、电磁波传播路径等来表现环境中的电磁分布情况，使用者根据已显示的图像与系统进行人机交互，对渲染图像可进行缩放、旋转、平移等操作，从多角度观察电磁场分布情况。

3.3 系统的类设计

根据系统的三层架构设计，系统解决方案下分三个命名空间 EM由数据访问(DA)、业务处理(BT)、使用者界面(UI)，分别对应系统的三层架构，每个类里又包括多个属性和方法。其中，处于底层的DA主要负责数据的读取和存储等相关操作，包括读取参数设置，读取DEM地形和空间剖面能量数据，读取WI的各文件存储路径及接收机位置、地形等;BT针对具体问题的业务逻辑处理，如数据映射转换操作、数据优化和等值线绘制等;UI实现空间可视化效果的显示及参数设置，并可实现显示或隐藏地形、电磁能量分布态势、电磁场等值线、电磁波传播路径，人机交互功能。主要包括系统的主界面、参数设置界面、二维平面俯视地形图用于设置发射机和接收范围的界面、帮助界面，如图9所示。

图9 电磁环境可视化系统主界面图

3.4 可视化效果

电磁场最主要的特征的便是能量的强弱，并且随着电磁波的传播，电磁能量也在变化，运用颜色、透明度等特征来表示电磁场能量的强弱这一属性，可以直观看出电磁场分布情况。以颜色表示能量值覆盖情况，根据覆盖点的颜色与能量标尺中的颜色对照即可得出该点的能量值，电磁能量分布态势图直观明了，方便全局的能量分布观测。电磁传播的各种可视化效果图如图10所示。图10(a)为空间剖面态势图，图10(b)与图10(c)分别为地面覆盖面电磁能量分布态势图和根据电磁能量分布绘制出的电磁场等值线图，白色点为发射机的位置;图10(d)为显示部分接收点的电磁波传播路径图［16］;图10(e)为电磁能量的点状分布，结合地形生成，宏观上展示了电磁的能量，达到电磁的三维模拟逼真显示效果，便于对电磁环境在整个空间的电磁能量强弱的判断;图10(f)为通过多个二维平面拟合出的三维电磁能量图。

图10 电磁传播的各种可视化效果图

4 结语

针对大范围的复杂电磁环境仿真问题，将三维电磁传播空间剖分成一个个的二维面，在每个面内再用2DPE法来并行计算电磁传播，这种准三维方法可大大简化计算过程，使计算的区域和复杂性都有很大提升，最后把得到的数据组合成三维仿真可视化效果显示，从而完成了EME三维可视化系统。这种方法将会对指导雷达基站选址、无线组网、优化配置无线通信网络资源、分析复杂环境中重要电子设备的EME效应等应用提供有效支撑。

尽管近几十年来EME并行仿真算法取得了较大进展，部分成果不仅成功应用在通信行业，在军事方面的应用优势也开始日益显现。但是面对复杂EME问题，尤其是不规则地形的大范围仿真时，还远远不能实现快速实时的计算和显示，无法满足实际工程的需求。并且对于陆地移动通信问题，尤其是在市区人口密集区，如何对环境建模变成了最大的考验，需要对这种情况提出一种合理的建模方法，相信随着仿真建模技术的不断发展，这种三维的实现方法将获得更多的应用。现代科学技术各学科之间的相互交叉、渗透，工程技术日趋集成化的特点，以及高性能计算机的飞速发展，正在进一步推动EME并行仿真的发展。

［1］Department of Defense Dictionary of Military and Associated Terms［Z］.2011.

［2］YU W，YANG X，LIU Y，et al.A Novel Hardware Acceleration Technique for High Performance Parallel FDTD Method［Z］.2011:3 154-3 157.

［3］TOPA T，KARWOWSKI A，NOGA A.Using GPU With CUDA to Accelerate MoM-Based Electromagnetic Simulation of Wire-Grid Models［J］.Antennas and Wireless Propagation Letters，IEEE，2011，10:342-345.

［4］ATHANAILEAS T E，ATHANASIADOU G E，TSOULOS G V，et al.Parallel Radio-Wave Propagation Modeling with Image-Based Ray Tracing Techniques［J］.Parallel Computing，2010，36(12):679-695.

［5］LIU D，ZHANG L，REN L，et al.Modeling and Simulation of Large-Domain Complex Electromagnetic Environment Based on Quasi-3D Method［C］//2011 International Conference in Electrics，Communication and Automatic Control Proceedings，New York，2012:1 357.

［6］LEVY M.Parabolic Equation Methods for Electromagnetic Wave Propagation［M］.London:Institution of Electrical Engineers，2000.

［7］胡绘斌.预测复杂环境下电波传播特性的算法研究［D］.长沙:国防科学技术大学，2006.

［8］AWADALLAH R S，GEHMAN J Z，KUTTLER J R，et al.Modeling Radar Propagation in Three-Dimensional Environments［J］.Johns Hopkins Apl Technical Digest，2004，25(2):101-111.

［9］CHEN P，WU L.3D Representation of Radar Coverage in Complicated Environment［J］.Simulation Modelling Practice and Theory，2008，16(9):1 190-1 199.

［10］SHEN Y，ZHANG L.Study on Parallel Methods of Three-Dimension Parabolic Equation for Complex Electromagnetic Environment Simulation［Z］.Genoa，Italy:2012.

［11］WU Y，ZHANG L，TAO F，et al.Three-Dimensional Electromagnetic Environment Modeling and High Performance Computing Analysis［Z］.Shanghai，China:2012.

［12］SHEN Y，ZHANG L，WU Y，et al.Methods to Improve Accuracy and Speed for the Quasi-3D Electromagnetic Environment Simulation［Z］.Shanghai，China，2012:53-59.

［13］MENDEZ D，REYES J A.Propagation of Electromagnetic Waves in Stochastic Helical Media［J］.2012:1-10.

［14］SHEN Y，ZHANG L，LIU D，et al.Comparisons of Ray-Tracing and Parabolic Equation Methods for the Large-Scale Complex Electromagnetic Environment Simulations［J］.International Journal of Modeling，Simulation，and Scientific Computing，2012:42-68.

［15］MU L，REN L，LIU D，et al.A Visualization Technology of Electromagnetic Environment Based on Parabolic Equation Method［C］.2011 International Conference in Electrics，Communication and Automatic Control Proceedings，Springer New York，2011.

［16］穆兰，任磊，吴迎年，等.空间电磁环境可视化系统的研究与应用［J］.系统仿真学报，2011(4):724-728.