陈 辉，张龙才

(桂林电子科技大学信息与通信工程学院，广西桂林541004)

0 引言

辐射电磁波在空间中自由传播，与媒质和媒质交界面发生相互作用，产生反射、折射、散射、绕射等现象称为电波传播。电波传播与实际地理模型构成的电磁环境，不仅在民用移动通信中起基础性作用，而且在当代电子、信息战中突显其巨大威力［1］。国际上最早是基于数理统计的方法对电波传播进行研究，其中典型的数学模型就是Okumura-Hata模型［2］以及后来修正过的COST 231模型。这些都属于经验模型，有一定的参考价值，但是不具有通用性。

后面出现的投射法(shooting and bouncing ray，SBR)是利用几何光学原理，将计算机图形学里的光线跟踪方法引入到无线电学来计算无线电信号的场强［3］。但科技工作者们发现，SBR与一致性绕射方法［4］(uniform theory of diffraction，UTD)在实际操作中，特别是出于计算机硬件的考虑，有一定的难度，射线跟踪不仅计算量大，而且对计算机的RAM，CPU主频要求都非常高，是一个非常耗时的算法。于是二叉树分区(binary space partitioning，BSP)、八叉树分区(Octal tree)、k 维树分区(K-D tree)［5］等一系列加速算法应运而生了，这些算法在一定程度上提高的计算速度。与此同时，基于开放消息传递(open message passing interface，openMPI)协议［6］的cluster集群计算也得到了一定的发展。

随着计算机硬件技术的发展，从硬件底层上就能实现并行计算的处理器将成为未来的发展方向。全球顶尖的GPU制造商NVIDIA于2007年6月成功研制出统一设备架构［7］(compute unified device architecture，CUDA)。新技术的推出，从底层硬件到上层开发平台都彻底改变了未来的硬件制造与软件开发方向，异构体并行计算已然成为未来高性能计算的主流。随着GPU可编程性的不断提高，利用GPU完成通用计算(general-purpose computing on graphics processing unit，GPGPU)的研究渐渐活跃起来，本文基于NVIDIA GeForce 9800GT(具有CUDA功能)硬件架构，在 Visual Studio 2010 IDE+CUDA Tool Kit 4.0开发平台下，采用串行射线跟踪与一致性绕射相结合的并行算法，实现城市小区电波传播的模拟预测。

1 CUDA加速运算物理机制

在WINTEL时代我们编写程序基本上没有考虑过硬件底层，也没有过多的考虑其耗时代价及RAM、寄存器等的使用情况。但现代应用软件程序量巨大，计算效率成了人们所关注的问题之一。

本次所用CUDA架构模块如图1所示，它是一种异构体架构，包括Host端(CPU端，本次实验是AMD Athlon 64 ×2 主频2.1 GHz，RAM为 3 GByte)与Device/Slave端(GPU端，本次实验是 NVIDIA GeForce 9800GT，GRAM为1 GByte)。由于 Device只进行并行通用处理，提高计算速度，还不能和外围设备交互，所以处理的结果要交给Host端，Host与Device之间通过PCI-E总线进行数据交流。其中GPU上含有多个 Block，每个 Block中，含有多个Thread(最小计算单元，为并行计算做了处理器的准备)，它们可以共用某一Block中的共享RAM。正因为这种底层构架才导致了编程模式从传统的CPU编程到异构体编程的转变。

基于CUDA技术就必须考虑硬件底层的资源问题，编写并行程序的主要的任务之一就是正确的分配NVIDIA GeForce 9800GT当中的Block，Thread执行相应的计算任务，确保各个Thread有条不紊的进行自己份内的任务，同时还需要各个Thread之间进行协调，确保步调一致。

图1 CUDA架构图Fig.1 Architecture of CUDA

2 软件算法的设计

2.1 并行算法

利用CUDA的目的是实现程序的并行化，在同一时刻启动多条射线跟踪以达到加速计算的目的。

如图2所示，图2为本次基于CUDA异构体模式下并行编程的总体框架图。图2中包括初始化模块和Kernel Core模块。初始化模块分为CPU，GPU空间分配，CPU数据生成，kernel函数的执行参数配置。Kernel Core模块分为辐射射线生成与跟踪、绕射射线跟踪、矢量叠加，最后做一些动态显示、释放内存等工作。

图2 并行计算总体框图Fig.2 Block diagram of parallel computing

2.2 并行处理器资源分配

NVIDIA GeForce 9800GT(支持CUDA计算)拥有112个计算core［8］，开发者需合理分配流处理器，从而达到最佳计算效率。并行射线跟踪任务主要由3个kernel函数来完成(注:虽然目前Fermi构架的GPU支持多个kernel函数并发执行，但是9800GT还没有这一功能，多个Kernel还是串行的)。

第1个kernel函数用于跟踪天线(本次实验假设天线为理想的全向单偶极子天线)辐射出的射线，存储好接收射线、绕射点信息。NVIDIA GeForce 9800GT硬件本身可支持最多512个threads/block，65 535个blocks/grid。结合硬件本身的限制，分配方法如下:在球坐标系下，不同的θ角(天顶角，0°≤θ＜＜180°)的射线分给不同的 block，在同一 θ角下的不同 φ 角(方位角，0°≤φ≤360°)分配给 block中不同的thread。其示意图如图3所示。

图3 并行处理器资源分配框图Fig.3 Block diagram of parallel processors resource allocation

第2个kernel函数将绕射点作为新的辐射源，然后跟踪每一条新的射线。分配方法是:一个block负责一个绕射点的所有绕射射线跟踪，此block中的所有thread负责绕射点绕射出的各条射线，m个block完成m个绕射点的跟踪。

第3个kernel函数将接收点处所有的射线的场强进行矢量叠加。

2.3 相关并行算法的伪码

现实中的电波传播是在三维空间中进行的，与常用的C/C++体系不同，CUDA的build-in内建变量[8]中给开发者提供了三维变量类型，不再需要开发者自己定义。比如三维空间中的点坐标、方向矢量等变量的数据类型在CUDA程序中可直接定义为float3类型。另外CUDA include头文件中也提供了一些常见三维计算函数供开发者直接调用，不必自写常用函数。如cutil_math.h［8］中提供了许多操作float3类型变量的函数。本次所需要调用的函数如下。

float dot(float3 a，float3 b)//点积

float3 cross(float3 a，float3 b)//叉积

float length(float3 v)//长度

float normalize(float3 v)//归一化

float3 reflect(float3 i，float3 n)//

反射向量计算，长度与入射向量相同。本次任务将射线及交点坐标等常规函数定义如下。

float3 in_origin//入射线起点

float3 in_direction//入射线方向

float3 reflect_origin//反射点

float3 reflect_direction//反射方向

float3 diffract_origin//绕射点

float3 building_norm//建筑物法线，用的是外法线

前面描述了多个流处理器 (stream processor，SP)如何有效的分配，但是当涉及到单个流处理器在进行计算的时候，还会有传统串行计算。下面基于单个SP分析其算法，其流程如图4所示。

图4 传统串行射线跟踪算法流程图Fig.4 Flow chart of traditional ray-tracing algorithm

从图4中可以看出射线跟踪算法耗时最多，计算最繁琐的要属射线与三维建筑物判交、求交了。利用遮挡测试［9］来提高计算效率，采用dot(float3 in_dirction，float3 building_norm)函数的结果作为判别依据，如果其值小于0，可能与建筑物相交继续往下判断。然后根据直线与平面是否相交以及交点是否在平面内的数学原理来计算相交点坐标并判断相交点是否在建筑物表面上［10］。如果交点落在建筑物平面上，则在建筑物表面内发生反射，相应的交点为反射点，反射线的矢量通过reflect(float3，in_direction，float3，building_norm)计算求得。如果交点落在建筑物棱上则发生绕射，相交点为绕射点。存储绕射点，子绕射射线的方向矢量由绕射前射线与建筑物的夹角以及划分射线的间隔所决定。

主程序的部分伪码如下。

1-2行:在内存、显存上开辟存储空间;3行:在host与device之间复制数据;4-5行:配置第一个kernel函数的维度(由于硬件所限，本文grid采用二维，block采用三维)

6行:调用并执行单偶极子天线发出的辐射射线跟踪函数kernel_1。

7行:释放内存、显存上所有占用存储单元，本次实验主要程序是前两个kernel函数，所以下面主要介绍前两个kernel的伪码。

辐射射线跟踪内核函数伪码部分如下。

1-3行:block、thread的索引值;

4-6行:设置index辐射射线的方位角与天顶角，并由θ和φ惟一确定了射线;

7-14:跟踪射线，若发生反射则求入射角、反射系数、反射向量，经过的路径长度及损耗，将反射点作为新的次级辐射源并继续跟踪，若发生绕射则存储绕射点及射线与障碍物边缘的夹角等参数，若射线被接收则存储该射线并结束，这部分与传统的射线跟踪一致。

绕射射线跟踪内核函数伪码部分如下。

判定该射线的传播路径，若反射则继续跟踪，若被接收则计算场强，存储该射线信息并结束。

3 实验结果与分析

与传统的经验模型不同，利用射线跟踪技术对小区电波场强预测时需要小区地理、地形数据的支持。地理数据包括建筑物、道路、树木、路灯、汽车、草坪等，这些施都会影响射线跟踪的传播路径。如果完全贴近现实世界，那样计算量无法估计，所以对三维世界进行必要的简化是合理的。我们主要关注建筑物与道路对电波传播的影响。其实电波传播的实质就是电磁波与三维空间中各种介质相互作用的过程。这种相互作用从物理学的角度讲，包括直射、反射、透射、绕射与散射等。本文采用长方体包围盒来简化建筑物的各种形态，用平行四边形来表示地面。

为了给人以直观的印象，利用标准图形OpenGL库进行可视化显示［11］。考虑计算量与计算机硬件配置，我们使用标准的多边形与多面体来描绘三维空间地理数据。先在视口中放置一个长方形(其规格为100 m×100 m)作为地面，然后在地面上放置4个建筑物，在正视图当中按逆时针算起，分别为标号为1-4。建筑物的几何结构分别为40 m×30 m×90 m，20 m ×20 m ×60 m，30 m ×20 m ×50 m，60 m×40 m×40 m。简单三维地理数据的正视图如图5所示。另外，设计了一个天线放在建筑物1的顶上，其具体的坐标位置为(30，90，0)。

图5 简单三维地理数据的正视图Fig.5 Front view of a simple 3 dimensional geography

以上介绍了三维建筑物的几何结构参数，下面将介绍与电波传播有关的电磁参数，本次实验采用的天线为单偶极子天线，其辐射功率为43 dBm，使用频段为1 800 MHz，极化方式为垂直极化，另外建筑物的相对介电常数εr=5.5，建筑物电导率σ=0.001，地面的相对介电常数 εr=3.4，地面的电导率为σ=0.002。为了方便观察，通过编写程序实现自由控制整个视口的旋转与缩放，从而可以从不同的视角来观察三维空间地理模型。经过旋转处理，得到的俯视图如图6所示。

利用前面的算法编写并行程序进行射线跟踪与CUDA加速之后，信号强度的分布图如图7所示。

从7图中可以看出，与天线处于视距，没有障碍物阻挡时，信号能够直达接收点，信号强度好，当电波被建筑物阻挡时，相对于天线不可见部分的信号强度明显较弱。这也类似于光学里面的阴影效应。实际当中，绕射现象是存在的，但是本次实验结果显示绕射现象不明显，这个跟绕射算法有关系。

图6 简单的三维地理数据的俯视图Fig.6 Vertical view of a simple 3 dimensional geography

图7 射线跟踪之后分布Fig.7 Electronic strength distribution after ray tracing

串行与并行算法耗时对比图如图8所示。

图8 串行与并行算法耗时对比图Fig.8 Comparison chart of computing time between serial and parallel algorithm

从图8中可以看出，当跟踪128×64条射线时，串行射线跟踪与并行射线跟踪所耗时相差不多，基本持平，甚至串行还比并行稍微用时少一点，这是因为CUDA中host与device之间数据的传输消耗了很多时间，对于跟踪较少射线时，并行算法的优势得不到体现，此时的效率比=并行耗时/串得耗时=1.07。当跟踪256×128条射线以及512×256条射线时，并行射线跟踪明显优于串行射线跟踪，这时大量的射线跟踪可以将host与device之间数据传输所消耗时间淹没掉，这时计算时间起主导作用，此时的效率比分别为:0.356，0.458，这说明基于 GPU 的通用计算需要良好的分配，才能达到较好的提速效果。

4 结论

CUDA技术是一个新生的并行计算技术，相对于传统的多核、Cluster集群并行计算更加优越。由实验结果可证明，相对于传统的全CPU编程，基于CUDA技术的无线电波传播模拟预测可以大幅度地缩减计算耗时，提高计算效率。从而也改变了人们的编程思维方式，有利于并行计算的发展。

［1］杨超，徐江斌，吴玲达.硬件加速的虚拟电磁环境体可视化［J］.北京邮电大学学报，2011，34(1):55-59.YANG Chao，XU Jiangbin，WU Lingda.Hardware Accelerated Volume Visualization in Virtual Electromagnetic Environment［J］.Journal of Beijing University of Posts and Telecommunications，2011，34(1):55-59.

［2］OKUMURA Y，OHMORI E，KAWANO T.Field Strength and its Variability in VHF and UHF Land Mobile Radio Service［J］.Rev Elec Commun，1969，16(1):825-873.

［3］ZHONG Ji，LI Binhong，WANG Haoxing.Efficient Raytracing Methods for Propagation Prediction for Indoor Wireless Communications［J］.IEEE Trans on AP，2001，43(2):41-49.

［4］KANATAS A，KOUNTOURIS I，KOSTARAS G.A UTD Propagation Model in Urban Microcellular Environments［J］.IEEE Transactions on Vehicular Technology，1997，46(1):185-193.

［5］SHEVTSOV Maxim，SOUPOKOV Alexei，KAPUSTIN Alexander.Highly Parallel Fast KD-Tree Construction for Interactive Ray Tracing of Dynamic Scenes［J］.Computer Graphics Forum，2007，26(3):395-404.

［6］KOMATITSCH Dimitri，ERLEBACHER Gordon.High-Order Finite-element Seismic Wave Propagation Modeling with MPI on a large GPU Cluster［J］.Journal of Computational Physics，2010，229(20):7692-7714.

［7］JASON Sanders，EDWARD Kandrot.GPU高性能编程CUDA实践［M］.聂学军，译.北京:机械工业出版社，2011.JASON Sanders，EDWARD Kandrot.CUDA By Example:an Introduction to General Purpose GPU Programming［M］.NIE Xuejun，Trans.Beijing:China Machine Press，2010.

［8］张舒.GPU高性能运算之CUDA［M］.北京:中国水利水电出版社，2009.ZHANG Shu.CUDA-high performance computing by GPU［M］.Beijing:China Water-Power Press，2009.

［9］郭建光.移动通信系统中微小区射线跟踪算法研究［M］.北京:北京邮电大学出版社，2010.GUO Jianguang.Microcellular Ray-tracing Algorithm in Mobile Communication System［M］.Beijing:Beijing U-niversity of Posts and Telecommunications Press，2010.

［10］袁正午.移动通信系统终端射线跟踪定位理论与方法［M］.北京:电子工业出版社，2007.YUAN Zhengwu.Theory and Methods of Ray-tracing and Location in Mobile Communication System［M］.Beijing:Publishing House of Electronics Industry，2007.

［11］OpenGL Architecture Review Board.OpenGL 编程指南［M］.第6版.徐波，译.北京:机械工业出版社，2008.OpenGL Architecture Review Board.OpenGL Programming Guide:The Official Guide to Learning OpenGL，Version 2.1［M］.Sixth Edition.XU bo，Trans.Beijing:China Machine Press，2008.