基于PE模型的三维广域电波传播预测方法

肖金光，周新力，宋斌斌，张 烨

（1.海军航空工程学院 研究生管理大队，山东 烟台264001；2.海军航空工程学院 电子信息工程系，山东 烟台264001）

0 引 言

军民应用日渐关注开放空间的电磁环境，其中的难点是地形和海上大气波导条件下电波传播问题。电波传播计算方法中的抛物方程 （parabolic equation，PE）法可以采用分步傅里叶变换 （split step fourier transform，SSFT）法步进求解，适于广域电磁环境求解［1－3］。美国基于PE 开发了 “高级折射效应预测系统” （advanced refractive effects prediction system，AREPS）和 “对流层电磁抛物方程程序” （tropospheric electromagnetic parabolic equation routine，TEMPER）［2］，但 均 为 二 维。而 与 数 字 高 程 模 型（digital elevation model，DEM）数据相结合的可视化更加便于非电磁专业人员对电磁环境的理解和应用。直接的3DPE技术在电磁环境中的应用尚有大量问题有待进一步研究［4］。文献 ［5］基于OpenMP线程并行计算多辐射源的二维传播。而PE的原理决定了在进行广域电波传播计算时需要大量的大矩阵分配，基于线程共享内存机制的OpenMP极易发生内存竞争，需要复杂的线程间内存共享分配控制以保证程序稳定。文献 ［6，7］基于高级传播模型 （advanced propagation model，APM）研究了雷达三维探测范围计算，并分别基于OpenGL和VTK 软件包显示，但在二维参数准备和三维数据生成方面可以进一步完善，且第三方软件显示便于可视化却不利于数据分析研究和应用。

针对广域三维电波传播预测在地理信息数据获取、并行计算、三维数据生成以及仿真研究的便利性方面存在的问题，本文提出了基于二维PE模型进行三维电磁环境电波传播预测和应用方法，从DEM 数据中抓取地形剖面并配伍相关参数，批处理并行运行PE，由多角度二维剖面生成三维电磁数据，多背景下参数等价表征法应用，设计了软件实现流程，并通过数值实验验证了方法的可行性。

1 基于抛物方程模型的电波传播计算

1.1 抛物方程模型［2，3］

笛卡尔坐标系下，时谐因子为exp（－iwt），二维标量波方程忽略后向传播并作近轴近似，可得标准抛物方程（standard parabolic equation，SPE）为

显然式 （2）可以步进迭代解算，两指数项分别表征折射和绕射效应。F 、F－1为正、反傅里叶变换，变换点数N由奈奎斯特准则确定

式中：zmax和pmax——需输出高度和变换域最大值。若β是以辐射源点为顶点的最大计算域锥角，相邻两点高度差Δz和距离上步进Δx 需满足

初始场由天线方向图的傅里叶逆变换求得，上边界采用虚部增量法，避免边界上的强反射 “污染”计算域［8］。下边界为光滑表面时，可由镜像原理将FFT 进一步简化为正弦变换或余弦变换，但对于粗糙表面或地形等阻抗边界，需要引入满足SPE解的辅助函数

式中：α——边界阻抗。根据Dirichlet边界条件，z ＝0 时w ＝0，场分量为u 的电磁波在阻抗边界上的传播就等效于为w 在光滑阻抗边界上的传播。引入式 （6）的差分形式求解PE，称为离散混合傅里叶变换 （discrete mixed Fourier transform，DMFT）法，边界阻抗与电波入射余角及地表类型有关。

1.2 剖面上任意位置数据的获得

由式 （4）、式 （5）可知，计算结果为剖面上的均匀网格，如果要求输出点（rout，zout）不在网格点上，直接线性插值不符合电波传播特性。

对于水平方向，将最后一个步进取为输出水平距离rout与Δx 的余数即可，这对于步进运算的SSFT 和DMFT 来说都是很容易实现的［9］

对于水平方向，输出点（rout，zout）在垂直方向的偏移量Δzshift为

则由傅里叶变换的性质可知

1.3 表征参数的等价换算

辐射源在自由空间距离x／m 处的场强E0／dBμ 为

式中：Pt／dBW 、Gt／dB 为辐射功率、发射天线增益。传播因子Fa定义为接收点场强E 和E0之比［2］

则传输损耗L／dB 为

根据辐射源的发射功率和收发天线增益，由式 （11）、式 （12）可 知 点 （x，z）上 信 号 的 场 强E／dBμ 和 功 率Pr／dBW 为

式中：Gr／dB ——接收天线增益。

2 三维广域电波传播算法

基于二维PE建立三维电磁环境，原理是在以辐射源为中心的多个角度二维剖面上，进行PE电波传播计算，生成三维数据，进行数据应用。

2.1 二维剖面参数准备

二维剖面参数包括地形／海洋数据、大气折射率剖面、辐射源参数、输出要求及PE 参数等，完成基于DEM 抓取地形剖面数据并进行地面类型和大气折射率剖面的配伍，将各剖面的输入参数存为约定格式文件，便于并行计算时独立运行。

2.1.1 输出参数设置

设定输出的最大垂直高度、最大水平距离和水平、垂直分辨率。除了PE模型本身的计算误差限制外，电波传播预测的精度与可获得的输入参数分辨率正相关，与计算中采用的水平垂直分辨率负相关，最终计算精度是与计算速度的妥协产物。

2.1.2 地形剖面的生成

地形剖面数越多，计算的空间精度越高，但是至少应该包含源点到地图四角所在的剖面。DMFT 算法对下边界阻抗边界的处理中的介电常数等参数与地面／水面 （海水、淡水）类型有关，因此剖面中地面高度和地表类型是两个重要参数。风作用下的水面相对于电波波长呈光滑或粗糙表面特性，主要参数为风速和风向。

首先是要在关注区域按照地形、水面 （湖泊、河流、海洋）等不同地表类型进行划分，这就相当于在DEM 数据上叠加与地形或水面区域相应的矢量数据。MATLAB软件自带的地图工具箱，具备高程数据和矢量数据两种类型数据的叠加功能。第二步是通过视在距离的确定，可以选取视距范围内的地形急剧变化点，提取从源点经过这些点的地形剖面，就转化为某角度地形剖面的提取。第三步是搜索与前述地形或水面矢量数据交叠区域，从而可以确定距离、高程、地表类型三者的对应关系。此时剖面的分辨率为DEM 的分辨率，在进行PE 计算时，PE 网格 （垂直FFT 点间距和水平步长决定）与DEM 的网格不一定匹配，需要采用双线性插值法产生PE网格上对应的高程数据。

以SRTM90 地 图 （网 址http：／／srtm.csi.cgiar.org）为例，取经纬度范围为：东经119°～125°，北纬35°～40°。区分陆地和海洋的矢量数据表征的就是陆海分界线，即为存储这一分界线的地理坐标。如果辐射源在经纬度坐标（121，37）处，以方位角0°方向为例，利用MATLAB 地图工具箱抓取地形剖面，地形切面如图1所示，其它任意角度的剖面图同理可得。

图1 分类地形剖面

在PE单剖面计算中，并不必知道经纬度，只需要与距离 （相对于源点）相关联的高程、地形类型、风速 （海洋）。图1只区分了陆地和海洋，进一步的地表类型细分可以提高计算准确度。

2.1.3 辐射源参数

辐射源参数包含功率、天线方向图、天线仰角、极化方式、安装位置等。如果实际计算中只有水平和垂直方向上的天线方向图，可以根据电波方向图一般具有对称性或准对称性，以地形剖面与垂直方向图的夹角大小进行量化。

2.1.4 大气折射指数剖面

式（2）中第一个指数项表征的是大气折射对电波传播的作用。大气折射指数变化的绝对值过小，为了便于反映大气折射作用，使用大气折射率M 表征，可以由气象测量获得。主要有标准大气和大气波导两类典型情况，后者对电波的异常传播影响较明显。如果大气折射率剖面在传播路径上变化，则必须与距离关联。

2.2 基于批处理的多剖面计算

根据设定的参数，基于PE 计算二维切面上的场强分布，涉及大量的数学公式，FORTRAN 语言可以进行高效的科学计算。计算结果存储为FORTRAN 直接访问文件或.mat文件。由于各个切面间的参数输入和计算是相对独立的，可以并行运算以缩短计算时间。

通常认为并行计算有基于进程的和基于线程的两类。批处理 （Batch）是一种简化的脚本语言，由系统内嵌的命令解释器解释运行，能够直接执行命令行，且windows程序都可以放在批处理文件中启动运行，通过Call命令顺序执行 （执行完毕才继续调用下一个程序）或Start命令并发执行 （依次调用但不等待前一程序执行完毕）。显然，通过Start命令执行多个单剖面PE 程序的方法，既简单又高效的实现了并行运算。

本文基于批处理的方式，利用操作系统自身的调度机制并发计算各个二维剖面的电波传播，是一种基于进程的并行计算方式。相比于多线程共享内存的OpenMP，优势在于防止了线程间内存的直接竞争，在内存不足时，利用系统自身的调度机制，简单有效解决了内存竞争问题，大幅降低了编程工作量，而且随着CPU 多核化和操作系统对多核的利用率越来越高，这种基于批处理的并行计算在简便性上将更加突出。

2.3 三维数据生成

以辐射源点所在高度轴为中心的多个二维剖面在空间上构成了分辨率有限的三维数据。已计算剖面上任意距离和高度点的场强可以根据式 （7）～式 （9）获得，不在已计算剖面上数据可以插值获得。方法为：①根据该点相对于与源点的方位角，确定处于那两个剖面之间；②进行坐标系转换，由于PE模型忽略了后向传播，此时决定电波强度的主要机制是电波的球面扩展，即主要因素是电波的传播距离，因此进行数据插值时最适宜的坐标系是球坐标系；③搜索需要查询点空间上紧邻的点，空间一点都会落在两个剖面各紧邻4个点围成的六面体中；最后采用三次样条插值得对应点的场强值。

2.4 数据应用

由式 （11）～式 （14）可知，场强、传播因子、传播损耗和功率四者一一对应。空间上某一点的传播因子、场强、接收功率 （此时还与接收天线参数有关）3 个参数的门限值，均可以换算为电波传播损耗门限。也就是说雷达探测范围、安全区域 （对人员、装备等）、通信距离覆盖（相对于特定灵敏度的接收机）、被探测或干扰的界定等应用，均可以等价的用对应的单路径传播损耗门限进行表征。在水平／垂直切面上该门限表现为等值线，遍历所有的输出高度，可得该门限在三维空间中的等值面，覆盖区域成为一个立体覆盖空间，所含区域就是探测范围、不安全区域、通信覆盖区域和干扰区域等。针对辐射源按照相应标准，建立等价门限库，便于开展上述应用。一个典型的应用是文献 ［10］提取外围轮廓线并基于OPEGL 研究了虚拟战场雷达探测范围和波束可视化方法。

2.5 软件流程设计

FORTRAN 在科学计算、MATLAB 在数据分析和可视化上各有优点，借助MATLAB 引擎和.mat文件，两者间的数据交互和程序调用非常方便。在软件实现中，使用FORTRAN 语言进行电磁计算，通过MATLAB 引擎调用地图工具箱抓取地形剖面，用Intel Visual FORTRAN 实现交互界面，PE计算结果存储为非格式化直接访问文件，提高内部访问效率［11］，同时通过MATLAB引擎存为.mat文件，便于利用MATLAB 丰富的函数库进行数据分析和图形绘制［12］。软件流程如图2所示。

图2 软件流程

3 数值实验

参数设置：选 （121，37.53）的局部高点为3GHz辐射源安放点，对应高程为221 m，全向天线的架设高度为20m，水平极化，标准大气，海面风速5m／s。以该辐射源点为中心选取地形剖面，选取16个典型方位角上的剖面，分别为0°～360°上均匀间隔30°的12个剖面，以及源点与地 图 四 角 所 在 的 剖 面 方 位 角 为：51.62°、128.12°、197.68°、342.48°。方位角0°方向上二维剖面上电波的传播损耗 （dB）如图3所示。

图3 方位角0°方向上电波传播损耗

不妨设某应用中的单路径传播损耗门限值为180dB，则在高度250m 和50m 上覆盖区域即为曲线包含区域 （白色五角星为辐射源所处位置）如图4、图5所示。

图4 高度250m 上损耗180dB对应的覆盖区域

图5 高度50m 上损耗180dB对应的覆盖区域

图4、图5白色线条包含范围清晰的表征了特定高度上电波180dB传播损耗所对应的覆盖区域，与地形条件紧密相关的电波直射、绕射等机制导致了该覆盖区域不规则分布。由两图可知最大作用距离在海洋方向上要大于陆地方向，这与辐射源点所处的地形有关系，例如图4中250m 高度略高于天线安装高度，辐射源180°方向上存在山峰阻挡，因此在该方向上作用距离较短。而图5中50m 高度低于天线安装高度，该高度上天线安装位置所在地局部高点附近为陆地，无有效的计算数据，因此看起来辐射源在覆盖区域外部。

4 结束语

针对广域三维电波传播计算和应用中存在的问题，本文提出了基于二维PE电波传播计算的三维电磁环境电波传播预测和应用方法，给出了广域三维电磁计算中从DEM 数据中提取地形剖面并与地表类型、大气折射率剖面、距离的配伍方法，提出了使用批处理技术并行运行各剖面上的PE计算，研究了由二维数据生成三维数据的原则与方法，对不同背景下的应用等价为单路径损耗表征，软件流程设计中充分利用了FORTRAN 和MATLAB 在高效科学计算和数据分析显示上的优势，实现的广域三维电波传播预测软件计算速度快、精度设置灵活、数据交互和分析应用方便、可视化友好，数值实验验证了方法的可行性。对于多个辐射源形成的综合场强、频谱分布和干扰问题将在下一步进行研究。

［1］WANG Hongguang，ZHANG Rui，KANG Shifeng，et al.Overview on parabolic equation model research for atmospheric ducting propagation ［J］.Equipment Environmental Engineering，2008，5 （1）：11－15 （in Chinese）. ［王红光，张蕊，康士峰，等.大气波导传播的抛物方程模型研究综述 ［J］.装备环境工程，2008，5 （1）：11－15.］

［2］Sirkova I.Brief review on PE method application to propagation channel modeling in sea environment ［J］.Central European Journal of Engineering，2012，2 （1）：19－38.

［3］XIAO Jinguang，ZHOU Xinli，LIU Xiaodi.Application of cured wave spectral estimation in wave propagation with parabolic equation model［J］.Chinese Journal of Radio Science，2014，29 （3）：559－566 （in Chinese）.［肖金光，周新力，刘晓娣.弯折谱估计在PE模型电波传播中的应用研究 ［J］.电波科学学报.2014，29 （3）：559－566.］

［4］Nunes da，Silva M，Costa E，et al.Analysis of the effects from lateral variations of irregular terrain based on a three－dimensional parabolic equation ［C］／／Proceedings of the 4th European Conference on Antennas ＆Propagation，2010：2－1.

［5］SHENG Nan，LIAO Chen，ZHANG Qinghong，et al.Approach for modeling 2－D radio wave propagation of multi－radiation sources based on OpenMP ［J］.Chinese Journal of Radio Science，2013，28 （4）：611－615 （in Chinese）. ［盛楠，廖成，张青洪，等.基于OpenMP 的多辐射源二维电波传播预测方法 ［J］.电波科学学报，2013，28 （4）：611－615.］

［6］ZHANG Jingzhuo，YUAN Xiujiu，ZHAO Xuejun，et al.3D Visualization of radar detection range based on advanced propagation model ［J］.Computer Engineering，2012，38 （4）：281－283 （in Chinese）.［张敬卓，袁修久，赵学军，等.基于APM 的雷达探测范围三维可视化 ［J］.计算机工程，2012，38 （4）：281－283.］

［7］ZHANG Jingzhuo，YUAN Xiujiu，ZHAO Xuejun，et al.3D detection range of radar in complex environment［J］.Journal of Computer Applications，2011，31 （10）：2738－2741 （in Chinese）.［张敬卓，袁修久，赵学军，等.复杂环境下雷达三 维 探 测 范 围 ［J］. 计 算 机 应 用，2011，31 （10）：2738－2741.］

［8］XIAO Jinguang，ZHOU Xinli，LIU Xiaodi.Image increment method on solving wave propagation problems with SSFT of PE［J］.Computer Engineering and Design，2014，35 （9）：3219－3223 （in Chinese）.［肖金光，周新力，刘晓娣.PE 模型SSFT 解电波传播问题的虚部增量法 ［J］.计算机工程与设计，2014，35 （9）：3219－3223.］

［9］ZHANG Qinghong，LIAO Cheng，SHENG Nan，et al.Non－uniform mesh techniques for parabolic equation ［J］.Chinese Journal of Radio Science，2013，28 （4）：636－640 （in Chinese）. ［张青洪，廖成，盛楠，等.抛物方程方法的非均匀网格技术研究［J］.电波科学学报，2013，28 （4）：636－640.］

［10］MA Bole，ZHANG Xiushan，ZHAO Binghua，et al.Study of visualization of radar range and beam in virtual battlefield［J］.Computer Engineering and Applications，2013，49 （9）：263－266 （in Chinese）.［马伯乐，张秀山，赵冰化，等.虚拟战场雷达作用范围与波束可视化方法研究 ［J］.计算机工程与应用，2013，49 （9）：263－266.］

［11］BAI Haibo.Authoritative guide on FORTRAN program designing ［M］.Beijing：China Machine Press，2013：279－281（in Chinese）. ［白海波.FORTRAN 程序设计权威指南［M］.北京：机械工业出版社，2013：279－281.］

［12］CHEN Jie.Bible of MATLAB ［M］.4th ed.Beijing：Electronic Industry Press，2013：844－848 （in Chinese）. ［陈杰.MATLAB宝典 ［M］.4版.北京：电子工业出版社，2013：844－848.］