射线跟踪技术的一种加速算法

周晓平 谭凤杰 柳朝阳 邹 澎

(1.郑州大学信息工程学院，河南 郑州 450052； 2.郑州大学数学系，河南 郑州 450052)

引 言

电波传播的准确预测对微蜂窝无线网络规划、设计和优化、基站选址和架设高度的选定、无线定位等都起着非常重要的作用，因而成为现代移动通信技术研究的热点之一.预测模型主要分为两类：统计模型和确定性模型.统计模型是经过对大量的测量数据进行统计分析得出的经验公式.一方面测量需要耗费大量的人力、物力和财力，另一方面由于微小区、微微小区等的出现，使不同小区之间的统计相似性消失，致使统计模型失效.确定性模型主要有射线跟踪、时域有限差分、MEI等方法.其中射线跟踪法根据几何光学理论和一致性绕射理论，跟踪直射线、反射线、绕射线等，寻找出收发之间所有的射线路径.然后，根据电波传播理论可以计算出每条射线的幅度、相位、延迟和极化等参数，并将通过各路径传来的电波进行矢量叠加就可以确定接收到的场强.

相对于FDTD算法，射线跟踪技术计算量较小，相对于MEI算法射线跟踪技术能进行三维计算，并且能考虑地形、地貌、地物等实际因素的影响，模拟多径传播路径，进行矢量叠加，故本文采用射线跟踪技术进行微小区电波传播预测.对于复杂的建筑物环境，直接遍历求解射线与面、劈、顶点的交点时，计算量极大，需要采用合适的射线跟踪加速算法.近年来，国内外学者提出的主要的加速算法有二元空间分区[1]、空间体积分区[1,3-4]和角度的z缓存区[5](angular z-buffer，AZB)、镜像法[6-7,9]、并行计算模型[8-11]等加速算法.

角度z缓存区算法的基本思路是把源或等效源所看见的空间分成角域，把角域所包含的多面体面储存在该角域中，并且按照它们与源或等效源的距离远近进行排序，对于进入该角域的射线，按照距离远近顺序进行与多面体面相交测试，若存在合理的交点，就不需再与排在其后的多面体面做相交测试，这样就可以减少必须做的相交测试数量.但是，这种算法是在空间内进行分区，并且当遇到反射等效源及绕射等效源时必须重新分区，同时在处理绕射问题时需要运用射线基坐标系，导致编程实现时比较复杂并且计算量也比较大，限制了跟踪速度的提高.

本文对AZB进行一些改进，对多面体在投影平面(地面上)内进行二维分区，在三维空间内进行跟踪，采用比背面采集算法计算量小的算法来排除被多面体自身遮挡的面、棱，及采用一种新算法对多面体进行排序.最后运用此改进算法，给出了微小区电波传播预测实例，并将预测结果与实测结果进行比较，一致性良好，证明了该算法的可行性.

1 射线跟踪角度缓存区改进算法

下面从辐射空间分区、抛弃自身遮挡的面及棱、前后排序、垂直平面内遮挡测试、水平棱绕射线跟踪等算法入手，对传统的角度缓存区算法进行一些改进，使其更易于程序实现，计算量更小，运算速度更快.城市中的建筑物大多数为平顶的多面体，故认为所讨论的绕射棱都垂直于或平行于地面.

下文中多面体指建筑物，辐射源指广义辐射源，即包括真实辐射源和等效辐射源，接收点指广义接收点，即包括真实接收点和等效接收点.射线跟踪过程中某次反射点或绕射点是上一级辐射源(包括等效辐射源)的等效接收点，同时也是下一次反射、绕射的等效辐射源.

1.1 基本方法

一辐射源位于空间中某点，以辐射源在地平面上的投影O为坐标原点，x、y坐标轴位于地平面上，以多面体的高度方向为z轴，建立直角坐标系O(x，y，z).则接收点的投影及空间中的多面体的底面均在xoy平面内.

1.1.1 基于平面内斜率角的射线搜索算法

1) 角扇区内多面体的搜索：设辐射源辐射空间内的第i个多面体底面的顶点与坐标原点的连线的最大、最小斜率分别为αi，max、αi，min，根据所有多面体的斜率将xoy平面分成若干个角扇区，尽量保证角扇区内的多面体的完整性，设第j个角扇区边线的斜率角为βj，max、βj，min.通过比较αi，max、αi，min和βj，max、βj，min的关系可以判断出某角扇区内包含的多面体.为了将不含多面体的角形区域排除掉，角扇区不一定相等.

2) 自身遮挡面、棱的判断：某反射面的法线向量用km表示，反射面辐射空间中的某多面体面的法线向量用kn表示，如果km·kn≥0，则该多面体面被多面体自身遮挡，否则，该面可以被源照射到，要搜寻的反射点和绕射点都位于该面及其棱上.通过这种方法可以抛弃多面体自身遮挡的面及棱，至少可以减小三分之一的计算量.

入射波能照射到的面和其上平行于地面及垂直于地面的棱为有效面、棱.

3) 排序算法：对于角扇区内多面体m与多面体n，设多面体m底面顶点与坐标原点之间连线斜率最大的连线OAm，max对应的角度为αm，max，斜率最小的连线OBm，min对应的角度为αm，min.对于多面体n有OAn，max、αn，max，OBn，min，αn，min.对于接收点，αm，max=αm，min，点Am，max与Bm，min重合.

如果αn，max≥αm，max，αm，min≥αn，min，则求出连线An，maxBn，min与连线OAm，max(或OBm，min)的交点C到源点O的距离dOC及dOBm，max(或dOBm，min).如果dOC≥dOAm，max(dOBm，min)，则多面体m比多面体n离源近，否则多面体n离源近.

如果αn，max≥αm，max，αn，min≥αm，min，则求出点Am，max和Bn，min到源点O的距离dOAm，max和dOBn，min，如果dOBn，min≥dOAm，max，则多面体m离源近，反之，多面体n离源近.

否则，两多面体之间无关，不需排序.

1.1.2 辐射源与接收点之间射线的遮挡测试

在接收点所在的角扇区中，对排列在辐射源与接收点之间的多面体，按照从前往后的顺序，首先在xoy平面内找到某多面体上可能遮挡射线的多面体面，然后仅对此面在垂直平面内做遮挡测试.具体算法是：对于某多面体，如果辐射源投影和接收点投影之间连线的斜率不在某多面体面投影的两个端点与辐射源投影之间连线的斜率之间，则该多面体面不遮挡辐射源到接收点之间的射线，否则该多面体面可能遮挡射线，只需对这种多面体面做遮挡测试，这样至少可以减小一半的计算量.

求可能遮挡射线的多面体面与射线的交点，如果交点不在该多面体面内，则该多面体面没有挡住射线，继续对下一个多面体面进行判断；如果交点在该多面体面内，则该多面体挡住了射线，抛弃该射线.只要有一个多面体面挡住了射线，则排在该多面体后面的多面体无需再做遮挡测试.这样可以更进一步减少计算量.如果辐射源与接收点之间的建筑物都没有挡住射线，则该射线有效.

1.2 直射线、反射线、绕射线的跟踪

1.2.1 直射线、一次入射线、地面反射线的跟踪

直射线、一次入射线的辐射源为发射天线，辐射区域为全空间，直射线的接收点为接收天线，入射线的接收点在多面体面或棱上.应用1.1.1的算法搜索有效反射面、绕射棱，并将其存入数据库中.

1) 对于直射线:找到接收天线所在的角扇区，根据1.1.2的算法对排列在发射天线与接收天线之间的多面体按顺序进行遮挡测试，如果没有多面体挡住直射线，则直射线有效.

2) 对于入射线:该数据库供一次反射线、绕射线跟踪使用.

3) 对于地面反射线:找到接收天线所在的角扇区，根据镜像原理，求出地面反射点，首先判断地面反射点是否在该角扇区内的某多面体底面内，如果在，则地面反射点无效，否则根据1.1.2的算法对排列在发射天线与地面反射点间的多面体依次进行入射线、反射线遮挡测试，如果没有多面体挡住射线，则地面反射有效.

1.2.2 反射线的跟踪

有效反射面将空间分成了两个区域，一个是辐射源所在的区域，另一个是辐射源镜像所在的区域.连接辐射源对于一次有效反射面的镜像的投影与该反射面投影的两端点，则两条连线与反射面的投影线之间的区域(辐射源所在的那一部分)就是反射面的角扇区，在此角扇区内的多面体能被反射射线照射到.

1) 一次有效反射：如果接收天线位于反射面的角扇区内，则根据镜像原理，求出一次反射点，根据1.1.2的算法依次对入射线、反射线进行遮挡测试，如果都未被遮挡，则一次反射有效.

2) 一次反射线成为二次反射或绕射的入射线：应用1.1.1的算法确定一次反射能照射到的有效多面体面及棱，并存入数据库中，成为二次反射、绕射的有效辐射源.供二次反射线、绕射线跟踪使用.

1.2.3 绕射射线的跟踪

1) 绕射棱的角扇区

对于垂直棱，将1.2.1中的有效垂直棱当作新的辐射源，他所在的多面体拐角的外角区域为其辐射空间，当外角超过90°时，将其划分成几个角扇区，否则角扇区就取为外角.

对于水平棱，如图1所示，求出源到水平棱的垂足G1，根据垂足的位置，可以将绕射线所能照射到的多面体所在的区域，分成2种情况：

图1 绕射线照射空间

② 如图1(b)所示，垂足G1在水平棱两端点之外，过离垂足G1近的端点A作垂直于水平棱AB的平面AA′CC′，该平面将全空间分为两个区域.取端点A为等效辐射源，由Keller 锥体2可知，辐射区域为水平棱AB所在的区域，对该区域再分角扇区.

位于角扇区中的某多面体，如果xoy平面内的顶点都满足km·kn≤0且高度低于该水平棱，则不能被绕射线照射到，否则多面体能被水平棱绕射线照射到，如图1所示.把不能被照射到的多面体排除掉，可以进一步减小计算量.

2) 跟踪射线

① 一次有效绕射：对于接收天线所在的角扇区，根据几何绕射理论，求出一次绕射点，根据1.1.2的算法依次对入射线、绕射线进行遮挡测试，如果都未被遮挡，则一次绕射有效.

② 一次绕射线成为二次反射或绕射的入射线：应用1.1.1的算法确定一次绕射能照射到的有效多面体面及棱，并存入数据库中，成为二次反射和绕射的有效辐射源.供二次反射线、绕射线跟踪使用.

1.3 多次反射、绕射射线的跟踪

与上述算法类似，可以找到两次以上的反射面、绕射棱，追踪到多次反射、多次绕射、及多次反射绕射的混合射线.

2 场强精确算法

寻找到所有的射线路径之后，根据辐射理论、反射定理及几何绕射理论，可以求得直射、反射及绕射场强.

移动通信基站天线的架设方式比较复杂，不但有三副天线而且带有下倾角及±45°极化角.其他用途天线的架设可以看成是它的特例，所以以移动通信基站为例研究场强的精确算法.

为了精确计算场强，必须考虑发射天线的下倾角、极化、波在传播过程中的极化面的改变等实际因素，由此需要建立以建筑物为参考的建筑物直角坐标系，以发射天线为参考的三个发射天线球面坐标系、直角坐标系，以射线为参考的若干个射线基坐标系[3].

算法思路是：对于远区的接收点，一般认为发射天线辐射球面波，所以首先将入射场强由发射天线球坐标系转换到发射天线直角坐标系，然后再变换到追踪到的入射线的射线基坐标系，根据反射场强公式[5]、绕射场强公式[4]分别计算出反射射线基坐标系、绕射射线基坐标系下的反射场、绕射场，由于不同的射线路径的射线基坐标方向不同，故需将场强再变换到发射天线直角坐标系，完成一副发射天线的场强合成，再将合场强转换到建筑物直角坐标系，最终将三副天线的场强统一到建筑物直角坐标系进行合成.

3 计算实例及结果分析

根据前面所论述的算法，对两个运营中的移动通信基站(GSM1800)周围小区中的射线进行跟踪及进行场强预测.由于所研究场景的占地面积较大，无法拍摄到全景图，故文中只列出了场景示意图.假设传播环境为平坦地、建筑物的表面是光滑平面.计算中，建筑物εr=4.5，σ=0.025 S/m，地面εr=5，σ=0.002 S/m.将位于所讨论环境中的建筑物及接收天线称作目标，软件采用Matlab实现，运算流程如图2所示.

图2 主程序流程图

为了验证算法的有效性，测试点应选择在没有汽车停放及基底场强小的地方.测量设备为SRM3000.

3.1 实例1-郑州大学柳园22号楼顶的通信基站

该基站周边建筑物均为四面体，如图4所示，基站安装在建筑物1的楼顶上，三副天线参数：发射功率20 W，增益18 dB，下倾角6°，极化±45°，架高35.1 m，第三副北偏东58°，距桅杆140 cm，第一副北偏西58°，距桅杆145 cm，第二副南偏东71°，距桅杆148 cm，接收点高1.03 m.测试频率为1 839 MHz.

该基站的西北角300 m左右处有另一个基站，且两个基站之间都是低矮建筑，致使本基站西北角附近基底场比较大，不能在此选点.在建筑物1的西、南、东、北方向各选10个点进行测量与计算，图3示出40个点的预测场强与实测结果对比曲线，图4示出其中四个测试点(第10、14、24、39点)的射线跟踪路径结果.

图3 实例1-预测场强与实测结果的对比

(a) 10点坐标(-59，0，1.03)

(b) 14点坐标(-21，-21，5，1.03)

(c) 24点坐标(21，6，1.03)

(d) 39点坐标(-17，33，1.03)图4 实例1-四个测试点的射线跟踪结果

3.2 实例2-郑州大学南校区图书馆周边小区

该基站周围建筑物1、4为12面体，2、3、5为6面体，如图6所示，基站位于建筑物1楼顶，三副天线参数为：发射功率20 W，增益18 dB，下倾角8°，±45°极化，架高28 m，第一副正西，第二副东偏北60°，第三副东偏南60°，三副天线距桅杆均为35 cm，接收点高1.1 m.测试频率为1 835 MHz.

该基站周边环境比较复杂，距该基站300 m左右的北方、西南方各有一基站，导致这些方向基底场比较大，不能选点.且本基站周围停放的大、小汽车较多，致使可选的测试点不多.在该基站东面、南面各选12个点进行测试与计算，图5示出24个点的预测场强与实测结果对比曲线，图6示出其中两个测试点的射线跟踪结果.

3.3 结果分析

每个观察点都有很多射线到达，如图4、6所示，有很多点没有直射线到达，只有绕射线，如图4(a)、(c)所示，计算场强与测试场强差别很小，说明多径反射、多径绕射对场强的贡献很大，不可忽略.

很多观察点没有三次以上的反射线和绕射线到达，如图4、图6所示，说明三次以上的反射线和绕射线出现的可能性较小，并且经过多次反射、绕射衰减使其对场强的贡献较小，可以忽略.

从图3、图5可以看出，实测值和计算值之间的误差在0.06～4.9 dB之间，实例1标准误差为2.834 6 dB，实例2标准误差为3.588 dB，误差较小，符合工程上的精度要求，说明运用该角度缓存区改进算法可以准确有效地预测小区环境中场强分布.

图5 实例2-预测场强与实测结果的对比

(a) 1点坐标(-43.43，12.93，1.1)

(b) 14点坐标(11.97,-31.07,1.1)图6 实例2-两个测试点的射线跟踪结果

为了验证该算法的运算速度，运用传统的角度缓存区算法对这两个基站也进行了计算，结果表明改进算法的运算速度提高了1倍左右.

由于忽略了树木、线杆、停放的汽车、自行车等、建筑物表面的粗糙度等引起的绕射与散射、及电参数的近似等，使实测值与计算值有微小的误差.

4 结 论

本文基于平面内直线斜率对广义辐射源辐射空间进行分区、对广义辐射源与广义接收点之间的建筑物进行前后排序、对被建筑物自身遮挡和建筑物之间互相遮挡的面及棱做出判断等，从而实现对传统的角度缓存区算法的改进.在场强的计算中也考虑了天线架设的实际因素.运用本文算法对两个通信基站环境中的场强进行了计算，同时用传统的角度缓存区算法也进行了计算，结果显示，本文提出的改进算法运算速度提高了一倍.另外对环境中的场强也进行了测量，测量结果与计算结果误差小于5 dB，说明该算法可以搜索到天线到接收机之间所有的传播路径，而且运算效率较高.本文所开发的软件可以运用于电磁环境测评.

[1] 吴志忠.移动通信无线电波传播[M].北京:人民邮电出版社，2002．

[2] 汪茂光.几何绕射理论[M].西安:西安电子科技大学出版社，1985．

[3] 程 勇，吴剑锋，曹 伟.一种用于移动系统场强预测的准三维射线跟踪模型[J].电波科学学报，2002，17(2)：15l-159．

CHENG Yong，WU Jianfeng，CAO Wei.A quasi 3-D ray-tracingmodel for field prediction in mobile systems[J].Chinese Journal of Radio Science，2002，17(2)：151-159.(in Chinese)

[4] 杨晋生，胡自胜，陈为刚.射线跟踪中三维矢量数据库建立方法[J].电信科学，2012，28(2)：110-114．

YANG Jinsheng，HU Zisheng，CHEB Weigang.A method of building 3D vector database for ray tracing[J].Telecommunications Science，2012，28(2)：110-114.

[5] CATEDRA M F，PEREZ J，SAEZ DE ADANA F，et al.Efficient ray-tracing techniques for three-dimensional analyses of propagation in mobile communications：application to picocell and microell scenarios[J].IEEE Antennas and Propagation Magazine，1998，40(2)：15-27．

[6] 袁正午，杨青宏，沐 维.一种改进的镜像射线跟踪方法[J].计算机工程，2012，38(3)：60-62．

YUAN Zhengwu，YANG Qinghong，mu Wei.improved image ray tracking method[J].Computer Engineering，2012，38(3)：60-62.(in Chinese)

[7] 袁正午.移动通信系统终端射线跟踪定位理论与方法[M].北京:电子工业出版社，2007．

[8] 王利东，李朝奎，陶建军，等.基于.NET Remoting的射线跟踪并行计算模型[J].计算机应用，2011，31(10)：2603-2605．

WANG Lidong，LI Chaokui，TAO Jianjun，et al.Parallel ray tracing model based on NET remoting[J].Journal of Computer Application，2011，31(10):2603-2605.(in Chinese)

[9] ATHANAILEAS T E，ATHANASIADOU G E，TSOULOS G V，et a1.Parallel radio-wave propagation modeling with image-based ray tracing techniques[J].Parallel Computing，2010，36(12)：679-695．

[10] 刘海涛，黎滨洪，谢 勇.并行射线跟踪算法及其在城市电波预测的应用[J].电波科学学报，2004，19(5)：581-585．

LIU Haitao LI Binhong, XIE Yong，et a1.Parallel ray-tracing algorithm and its application for propagation prediction in urban microceilular environments[J].Chinese Journal of Radio Science，2004，19(5)：581-585.(in Chinese)

[11] 杨锦辉，张 文.一种基于三维射线弹跳法的并行电磁传播预测算法[J].装备指挥技术学院学报，2010，21(6)：91-96．

YANG Jinhui，ZHANG Wen.A new parallel radio propagation prediction algorithm based on 3D SBR technique[J].Journal of the Academy of Equipment Command & Technology，2010，21(6)：91-96.(in Chinese)