对于海面舰船目标电磁散射特性而言，除了一次散射占较大比重外，舰船的舰面布置复杂，桅杆、开口结构等间相互作用，是改变舰船雷达波多次散射特征复杂的重要原因。因此，需要通过采用射线追踪法计算多次散射贡献。在射线追踪过程中，需要记录射线的轨迹，计算射线在分界面处产生的反射场信息。基于物理光学法(physical optics, PO)理论，分界面处等效电磁流可以由射线的入射场、反射场共同确定，进而再通过Stratton-Chu公式求得散射场。因此，该高频方法被标记为“PO-SBR”法^[9]。

由于海面舰船的超电大规模，其电磁射线管数量巨大，需要分批次处理。基于区域射线束的任务分配策略首先应确定计算区域。由于目标处于远场区，射线追踪范围即为目标在垂直于入射方向上的射线管平面投影。建立以${\boldsymbol{o}} = ( r\sin \theta \cos \varphi ,r\sin \theta \sin \varphi , r\cos \theta )$为原点，以方向矢量${\boldsymbol{u}} = ( - \cos \theta \cos \varphi , - \cos \theta \sin \varphi , \sin \theta )$，${\boldsymbol{v}} = \left( { - \sin \varphi ,\cos \varphi ,0} \right)$，${\boldsymbol{w}} = ( - \sin \theta \cos \varphi , - \sin \theta$ $\sin \varphi , \cos \theta )$为坐标轴的局部坐标系，$\theta $，$\varphi $分别为入射俯仰角与方位角。则空间中任意点${\boldsymbol{r}}$在射线束平面uov的投影可确定为：

式中，${\boldsymbol{r'}}$为${\boldsymbol{r}}$投影于uov平面的位置矢量。为简化投影计算，构建以目标在各坐标轴上最大坐标及最小坐标为基准构建的轴对齐包围盒(AABB)，由包围盒8个顶点在uov的投影可得到射线追踪区域。

基于区域射线束分割策略，对剖分后的射线进行编号，则第$ i\left( {i = 0,1,2,\cdots} \right) $个射线中心的局部坐标为：

式中，${u_{\max }}$、${u_{\min }}$、${v_{\max }}$、${v_{\min }}$分别表示投影在局部坐标轴${\boldsymbol{u}}$、${\boldsymbol{v}}$处的边界值；$\varDelta$表示网格剖分密度；$\left\lceil \cdot \right\rceil $表示向上取整；$\left\lfloor \cdot \right\rfloor $表示向下取整；$\bmod \left( \cdot \right)$表示取模运算。

由于弹跳射线法以单一射线为最小计算单位，通过对照射目标射线束的追踪(ray tracing, RT)实现散射计算，以追踪区域中射线集合为单位进行划分并封装为子任务。需要保证散射体在射线管平面的投影全部落在射线管平面区域中，需要根据设定的射线管密度将射线管平面剖分为大小均匀的矩阵网络，通常设定网格剖分密度$\varDelta \leqslant {\lambda \mathord{\left/ {\vphantom {\lambda {10}}} \right. } {10}}$，在射线管平面上构成$M \times N$个射线管，且每个矩形网格的中心为射线管发射起点，射线追踪示意图如图1所示。

图  1  射线追踪示意图

由于目标处于远场区，射线追踪范围即为目标在垂直于入射方向上的射线管平面投影。由于GPU中的全局内存空间有限，全局内存空间无法完整地存储所有射线管的每一级射线所包含的信息。在射线追踪的过程中，每一批次的核函数最多调用$ T $个射线追踪任务。

在$M \times N$个射线管构成的矩形网格中，对矩形网格进行按照如下分割规则进行划分：不失一般性，假设$M \leqslant N$，交替选取未被划分的矩形网格中最长边${M_i}$、最短边${N_i}$，其中下标$ i $表示其为第$ i $次划分的边长；特别地，取决于第一次划分从长边或短边开始，有${M_0} = M$或${N_0} = N$。将选取的边长作为矩形的一条边长，在网格中截取一个极大的矩形网格分区，其面积等于T。每次选取的矩形网格分区至少能够覆盖${M_i} \times 1$或${N_i} \times 1$个网格单元，将尺寸为$M \times N$的矩形网格分割成若干个矩形分区，每一矩形分区中的所有射线管将作为一个批次的射线管任务，如图2所示。第i次GPU全局内存利用率为：

式中，$\left\lfloor \cdot \right\rfloor $表示向下取整。数学理论证明，上述方法能够使得${\varsigma _i}$最大，从而使得整体内存利用率$\varsigma = \prod {{\varsigma _i}}$最优，提高算法的整体计算效率。在程序调用核函数时，将矩形分区的其中一个包含射线管数目$ T $的分区作为一个批次。每批次调用核函数之前，将当前分区的射线管位置、长度、宽度等数据存储到全局内存中，并在核函数调用时将此分区射线管的位置数据指针传入核函数中。在核函数中使用线程编号就可以定位每一个线程所对应的射线管的位置，这样就完成了射线管映射到每一个线程的过程。

图  2  大规模射线管的任务分割方法

考虑到中央处理器(CPU)与GPU并行计算能力存在极大的差距，实际程序中繁重的射线追踪及电磁计算任务将全部交由GPU处理，CPU端仅对程序进行整体控制，无需执行高耗时操作，因此采用MPS并发控制GPU计算^[10]。

基于MPS的Multi-GPU框架由任务调度线程及GPU控制线程组成。通过读取主机中挂载GPU设备信息进而为创建thread线程对象实现Multi-GPU控制，各线程利用cudaSetDevice( )指定对应的设备加以控制，MPS将为线程分配唯一的上下文(Context)，该上下文始终指向绑定的GPU且包含此设备运行所需的内存、模块(Modules)、数据流(Streams)、本地内核(Kernel)缓存等信息，CUDA根据上下文信息处理并控制对应GPU的运行。主程序完成线程分配后执行系统调度，负责将划分后的任务集合中的追踪任务分配至不同设备进行处理，Multi-GPU并行框架如图3所示。

图  3  Multi-GPU并行框架

客观上不同GPU间性能存在一定差异，如若单元间无法实现协同调配，可能导致部分GPU空转闲置的现象发生，程序性能急剧退化。现有的负载均衡算法大体可以分为动态及静态两类。静态负载均衡算法依照固化完成的任务分配策略对系统进行调度，其运行时不依赖运算单元状态，算法结构简单，但无法适应环境变化所带来的负载失衡。

动态负载均衡算法随运算单元状态调整任务分配，为保证代码健硕性，避免线程冲突，利用原子操作的不可分割性构建自旋锁对临界区代码进行保护。C++标准库中所带的原子量(Atomic)为实现提供了便利，其中比较交换操作(CAS)操作函数中的compare_exchange_strong()与store()函数为多线程提供的原子操作，保证线程安全。其中前者用于比较原子量数值是否符合预期，后者用于为原子量赋值，同时每个线程都将拥有atomic<bool> isComplete任务状态标志位。调度过程主要包括：1)调度器对各线程isComplete标志位进行轮询，通过调用compare_exchange_strong()函数获取标志位数值，当检测到GPU完成任务后，调度器检查任务集合并重新推送新任务至线程，并通过store()函数清除isComplete位。2)当线程接到标志位清除信号，读入待处理的追踪任务并开始追踪计算，计算完成后置位isComplete。3)调度器检测到该标志位后重新推送新任务至线程。循环执行上述调度过程。若任务池中已无待处理数据，调度器等待所有isComplete置位后获取散射场计算结果并销毁线程。

为了展示基于Multi-GPU并行技术的高频方法的计算能力，首先对本方法的精确度进行验证，随后比较不同目标尺寸电磁散射计算单GPU与多GPU 的加速效果，最后利用基于Multi-GPU的高频方法对2 km×2 km大小的航母舰队电磁散射进行分析。其中，运行平台基本硬件参数为：Intel Xeon E5-2620 V4 CPU，双路Nvidia Titan XP GPU。舰船模型主要包括：1)简化船模型，尺寸27 m(L)×3.2 m(W)×5 m(H)；2)货船，尺寸140 m(L)×18 m(W)×17 m(W)；3)监测船，尺寸132.8 m(L)×20 m(W)×23.3 m(H)；4)航母，尺寸347.5 m(L)×87.8 m(W)×58.4 m(H)，如图4所示。

图  4  舰船模型示意图

图  5  射线追踪法精度验证

为了验证Multi-GPU加速的PO-SBR算法精度及性能，利用本方法与商业软件FEKO的数值方法——多层快速多极子(MLFMM)对如图4所示简化船模型及风驱海面进行电磁散射计算^[11]。其中，雷达入射波频率为400 MHz，入射方位角0°，俯仰角−89°～89°。限于FEKO难以应对电大尺寸目标电磁散射计算，海面大小为48 m (L)×16 m (W)，海面风速为5 m/s，相对介电常数${\varepsilon _r} = 4 - {\rm{j}}8$。

MLFMM 计算使用12线程，PO-SBR分别采用单GPU 与多GPU进行加速。其中，MLFMM平均每角度计算耗时395.2 s，单GPU下PO-SBR平均每角度计算耗时0.041 s，双GPU下PO-SBR平均每角度计算耗时0.028 s。可以看到GPU加速下的PO-SBR与MLFMM相比存在明显的速度优势，对于双GPU而言，由于目标场景尺寸较小，计算耗时在总耗时中占比低，所以相较于单GPU加速而言效果无法有效体现。同时发现，单GPU模式与双GPU模式的计算结果完全吻合，说明了计算精度的一致性。

为了进一步确认本文算法与方案的计算效率，采取如下两种场景进行分析说明。

场景1：货船编队，其中：舰船编队由6艘货船组成，舰船前后间隔200 m、左右间隔50 m。海面大小分别为：400 m (L)×400 m (W)，800 m (L)×800 m (W)，1000 m (L)×1000 m (W)，海面风速为5 m/s。船舰编队进行单站散射计算，其中入射波频率为3 GHz，入射方位角90°，俯仰角0°～89°。

场景2：大型航母与监测船编队，该战斗群由2艘航母及6艘监测船构成。电磁波工作频率为10 GHz，入射俯仰角为45°，方位角为0°～360°，海面尺寸为2000 m (L)×2000 m (W)，海面风速为5 m/s，场景如图6所示。

图  6  大型航母与监测船编队场景

针对场景1的仿真，不同海面对应不同的电尺寸规模，表1为不同海面尺寸下两种方法的计算耗时比较。其中，1000 m×1000 m的海面对应$ 10\;000\lambda \times 10\;000\lambda $规模大小，其计算总耗时为880 s，平均单角度耗时9.78 s，两种电磁极化结果如图7a所示。双GPU加速下平均加速比均达到甚至超过200%。双GPU模型下的加速比超过200%的原因为：上述Multi-GPU加速框架中独立的主线程使得任务剖分为异步处理过程，换而言之在Multi-GPU中任务剖分与射线追踪是同一时刻进行的，任务剖分耗时也由此节省。因此当计算量增大时，Multi-GPU的加速比也随之增加。

图  7  两种场景不同极化散射结果

针对场景2的仿真，其电尺寸规模远超过$10\;000\lambda $。采用单GPU与双GPU两种不同方法计算航母战斗群，最终RCS结果保持一致，如图7b所示。在计算耗时上，单GPU平均每角度计算耗时1 398.19 s，双GPU平均每个观测角度计算耗时653.59 s，加速比为213%。

在上述分析中可以发现：单GPU计算难以解决上述场景耗时高的问题。Multi-GPU由于其自生的可扩展性，通过扩充GPU数量可以保证对任意真实场景的散射计算时间在可容忍范围内。相比于常见的CPU计算集群能够用更少的设备达到更大的加速比。因此相较于CPU集群Multi-GPU方案在处理并行问题上更具性价比。

实际雷达工程中，受照射目标场景大小通常可达数万倍波长，对该场景的散射计算属于超电大尺寸问题，需要通过计算机并行技术加速算法实现。由于高频方法的局部性原则使其具有良好的并行性，在处理此类问题时，GPU相比于CPU展示出更高的性价比。现有的GPU加速通常为单设备加速，在处理超电大尺寸问题耗时仍难以满足要求，因此研究Multi-GPU加速的实现以动态提高算力是必要的。本文对海面船舰场景下的高频方法Multi-GPU并行加速技术的实现进行了研究。结果表明：双GPU相比于单GPU算法平均加速200%左右，可有效应用于超大场景下的电磁散射计算，具有明显的应用潜力。