基于CUDA 技術(shù)城市小區(qū)電波傳播并行計(jì)算的研究

陳 輝，張龍才

(桂林電子科技大學(xué)信息與通信工程學(xué)院，廣西桂林541004)

0 引言

輻射電磁波在空間中自由傳播，與媒質(zhì)和媒質(zhì)交界面發(fā)生相互作用，產(chǎn)生反射、折射、散射、繞射等現(xiàn)象稱(chēng)為電波傳播。電波傳播與實(shí)際地理模型構(gòu)成的電磁環(huán)境，不僅在民用移動(dòng)通信中起基礎(chǔ)性作用，而且在當(dāng)代電子、信息戰(zhàn)中突顯其巨大威力［1］。國(guó)際上最早是基于數(shù)理統(tǒng)計(jì)的方法對(duì)電波傳播進(jìn)行研究，其中典型的數(shù)學(xué)模型就是Okumura-Hata模型［2］以及后來(lái)修正過(guò)的COST 231模型。這些都屬于經(jīng)驗(yàn)?zāi)Ｐ?，有一定的參考價(jià)值，但是不具有通用性。

后面出現(xiàn)的投射法(shooting and bouncing ray，SBR)是利用幾何光學(xué)原理，將計(jì)算機(jī)圖形學(xué)里的光線(xiàn)跟蹤方法引入到無(wú)線(xiàn)電學(xué)來(lái)計(jì)算無(wú)線(xiàn)電信號(hào)的場(chǎng)強(qiáng)［3］。但科技工作者們發(fā)現(xiàn)，SBR與一致性繞射方法［4］(uniform theory of diffraction，UTD)在實(shí)際操作中，特別是出于計(jì)算機(jī)硬件的考慮，有一定的難度，射線(xiàn)跟蹤不僅計(jì)算量大，而且對(duì)計(jì)算機(jī)的RAM，CPU主頻要求都非常高，是一個(gè)非常耗時(shí)的算法。于是二叉樹(shù)分區(qū)(binary space partitioning，BSP)、八叉樹(shù)分區(qū)(Octal tree)、k 維樹(shù)分區(qū)(K-D tree)［5］等一系列加速算法應(yīng)運(yùn)而生了，這些算法在一定程度上提高的計(jì)算速度。與此同時(shí)，基于開(kāi)放消息傳遞(open message passing interface，openMPI)協(xié)議［6］的cluster集群計(jì)算也得到了一定的發(fā)展。

隨著計(jì)算機(jī)硬件技術(shù)的發(fā)展，從硬件底層上就能實(shí)現(xiàn)并行計(jì)算的處理器將成為未來(lái)的發(fā)展方向。全球頂尖的GPU制造商N(yùn)VIDIA于2007年6月成功研制出統(tǒng)一設(shè)備架構(gòu)［7］(compute unified device architecture，CUDA)。新技術(shù)的推出，從底層硬件到上層開(kāi)發(fā)平臺(tái)都徹底改變了未來(lái)的硬件制造與軟件開(kāi)發(fā)方向，異構(gòu)體并行計(jì)算已然成為未來(lái)高性能計(jì)算的主流。隨著GPU可編程性的不斷提高，利用GPU完成通用計(jì)算(general-purpose computing on graphics processing unit，GPGPU)的研究漸漸活躍起來(lái)，本文基于NVIDIA GeForce 9800GT(具有CUDA功能)硬件架構(gòu)，在 Visual Studio 2010 IDE+CUDA Tool Kit 4.0開(kāi)發(fā)平臺(tái)下，采用串行射線(xiàn)跟蹤與一致性繞射相結(jié)合的并行算法，實(shí)現(xiàn)城市小區(qū)電波傳播的模擬預(yù)測(cè)。

1 CUDA加速運(yùn)算物理機(jī)制

在WINTEL時(shí)代我們編寫(xiě)程序基本上沒(méi)有考慮過(guò)硬件底層，也沒(méi)有過(guò)多的考慮其耗時(shí)代價(jià)及RAM、寄存器等的使用情況。但現(xiàn)代應(yīng)用軟件程序量巨大，計(jì)算效率成了人們所關(guān)注的問(wèn)題之一。

本次所用CUDA架構(gòu)模塊如圖1所示，它是一種異構(gòu)體架構(gòu)，包括Host端(CPU端，本次實(shí)驗(yàn)是AMD Athlon 64 ×2 主頻2.1 GHz，RAM為 3 GByte)與Device/Slave端(GPU端，本次實(shí)驗(yàn)是 NVIDIA GeForce 9800GT，GRAM為1 GByte)。由于 Device只進(jìn)行并行通用處理，提高計(jì)算速度，還不能和外圍設(shè)備交互，所以處理的結(jié)果要交給Host端，Host與Device之間通過(guò)PCI-E總線(xiàn)進(jìn)行數(shù)據(jù)交流。其中GPU上含有多個(gè) Block，每個(gè) Block中，含有多個(gè)Thread(最小計(jì)算單元，為并行計(jì)算做了處理器的準(zhǔn)備)，它們可以共用某一Block中的共享RAM。正因?yàn)檫@種底層構(gòu)架才導(dǎo)致了編程模式從傳統(tǒng)的CPU編程到異構(gòu)體編程的轉(zhuǎn)變。

基于CUDA技術(shù)就必須考慮硬件底層的資源問(wèn)題，編寫(xiě)并行程序的主要的任務(wù)之一就是正確的分配N(xiāo)VIDIA GeForce 9800GT當(dāng)中的Block，Thread執(zhí)行相應(yīng)的計(jì)算任務(wù)，確保各個(gè)Thread有條不紊的進(jìn)行自己份內(nèi)的任務(wù)，同時(shí)還需要各個(gè)Thread之間進(jìn)行協(xié)調(diào)，確保步調(diào)一致。

圖1 CUDA架構(gòu)圖Fig.1 Architecture of CUDA

2 軟件算法的設(shè)計(jì)

2.1 并行算法

利用CUDA的目的是實(shí)現(xiàn)程序的并行化，在同一時(shí)刻啟動(dòng)多條射線(xiàn)跟蹤以達(dá)到加速計(jì)算的目的。

如圖2所示，圖2為本次基于CUDA異構(gòu)體模式下并行編程的總體框架圖。圖2中包括初始化模塊和Kernel Core模塊。初始化模塊分為CPU，GPU空間分配，CPU數(shù)據(jù)生成，kernel函數(shù)的執(zhí)行參數(shù)配置。Kernel Core模塊分為輻射射線(xiàn)生成與跟蹤、繞射射線(xiàn)跟蹤、矢量疊加，最后做一些動(dòng)態(tài)顯示、釋放內(nèi)存等工作。

圖2 并行計(jì)算總體框圖Fig.2 Block diagram of parallel computing

2.2 并行處理器資源分配

NVIDIA GeForce 9800GT(支持CUDA計(jì)算)擁有112個(gè)計(jì)算core［8］，開(kāi)發(fā)者需合理分配流處理器，從而達(dá)到最佳計(jì)算效率。并行射線(xiàn)跟蹤任務(wù)主要由3個(gè)kernel函數(shù)來(lái)完成(注:雖然目前Fermi構(gòu)架的GPU支持多個(gè)kernel函數(shù)并發(fā)執(zhí)行，但是9800GT還沒(méi)有這一功能，多個(gè)Kernel還是串行的)。

第1個(gè)kernel函數(shù)用于跟蹤天線(xiàn)(本次實(shí)驗(yàn)假設(shè)天線(xiàn)為理想的全向單偶極子天線(xiàn))輻射出的射線(xiàn)，存儲(chǔ)好接收射線(xiàn)、繞射點(diǎn)信息。NVIDIA GeForce 9800GT硬件本身可支持最多512個(gè)threads/block，65 535個(gè)blocks/grid。結(jié)合硬件本身的限制，分配方法如下:在球坐標(biāo)系下，不同的θ角(天頂角，0°≤θ＜＜180°)的射線(xiàn)分給不同的 block，在同一 θ角下的不同 φ 角(方位角，0°≤φ≤360°)分配給 block中不同的thread。其示意圖如圖3所示。

圖3 并行處理器資源分配框圖Fig.3 Block diagram of parallel processors resource allocation

第2個(gè)kernel函數(shù)將繞射點(diǎn)作為新的輻射源，然后跟蹤每一條新的射線(xiàn)。分配方法是:一個(gè)block負(fù)責(zé)一個(gè)繞射點(diǎn)的所有繞射射線(xiàn)跟蹤，此block中的所有thread負(fù)責(zé)繞射點(diǎn)繞射出的各條射線(xiàn)，m個(gè)block完成m個(gè)繞射點(diǎn)的跟蹤。

第3個(gè)kernel函數(shù)將接收點(diǎn)處所有的射線(xiàn)的場(chǎng)強(qiáng)進(jìn)行矢量疊加。

2.3 相關(guān)并行算法的偽碼

現(xiàn)實(shí)中的電波傳播是在三維空間中進(jìn)行的，與常用的C/C++體系不同，CUDA的build-in內(nèi)建變量[8]中給開(kāi)發(fā)者提供了三維變量類(lèi)型，不再需要開(kāi)發(fā)者自己定義。比如三維空間中的點(diǎn)坐標(biāo)、方向矢量等變量的數(shù)據(jù)類(lèi)型在CUDA程序中可直接定義為float3類(lèi)型。另外CUDA include頭文件中也提供了一些常見(jiàn)三維計(jì)算函數(shù)供開(kāi)發(fā)者直接調(diào)用，不必自寫(xiě)常用函數(shù)。如cutil_math.h［8］中提供了許多操作float3類(lèi)型變量的函數(shù)。本次所需要調(diào)用的函數(shù)如下。

float dot(float3 a，float3 b)//點(diǎn)積

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

float length(float3 v)//長(zhǎng)度

float normalize(float3 v)//歸一化

float3 reflect(float3 i，float3 n)//

反射向量計(jì)算，長(zhǎng)度與入射向量相同。本次任務(wù)將射線(xiàn)及交點(diǎn)坐標(biāo)等常規(guī)函數(shù)定義如下。

float3 in_origin//入射線(xiàn)起點(diǎn)

float3 in_direction//入射線(xiàn)方向

float3 reflect_origin//反射點(diǎn)

float3 reflect_direction//反射方向

float3 diffract_origin//繞射點(diǎn)

float3 building_norm//建筑物法線(xiàn)，用的是外法線(xiàn)

前面描述了多個(gè)流處理器 (stream processor，SP)如何有效的分配，但是當(dāng)涉及到單個(gè)流處理器在進(jìn)行計(jì)算的時(shí)候，還會(huì)有傳統(tǒng)串行計(jì)算。下面基于單個(gè)SP分析其算法，其流程如圖4所示。

圖4 傳統(tǒng)串行射線(xiàn)跟蹤算法流程圖Fig.4 Flow chart of traditional ray-tracing algorithm

從圖4中可以看出射線(xiàn)跟蹤算法耗時(shí)最多，計(jì)算最繁瑣的要屬射線(xiàn)與三維建筑物判交、求交了。利用遮擋測(cè)試［9］來(lái)提高計(jì)算效率，采用dot(float3 in_dirction，float3 building_norm)函數(shù)的結(jié)果作為判別依據(jù)，如果其值小于0，可能與建筑物相交繼續(xù)往下判斷。然后根據(jù)直線(xiàn)與平面是否相交以及交點(diǎn)是否在平面內(nèi)的數(shù)學(xué)原理來(lái)計(jì)算相交點(diǎn)坐標(biāo)并判斷相交點(diǎn)是否在建筑物表面上［10］。如果交點(diǎn)落在建筑物平面上，則在建筑物表面內(nèi)發(fā)生反射，相應(yīng)的交點(diǎn)為反射點(diǎn)，反射線(xiàn)的矢量通過(guò)reflect(float3，in_direction，float3，building_norm)計(jì)算求得。如果交點(diǎn)落在建筑物棱上則發(fā)生繞射，相交點(diǎn)為繞射點(diǎn)。存儲(chǔ)繞射點(diǎn)，子繞射射線(xiàn)的方向矢量由繞射前射線(xiàn)與建筑物的夾角以及劃分射線(xiàn)的間隔所決定。

主程序的部分偽碼如下。

1-2行:在內(nèi)存、顯存上開(kāi)辟存儲(chǔ)空間;3行:在host與device之間復(fù)制數(shù)據(jù);4-5行:配置第一個(gè)kernel函數(shù)的維度(由于硬件所限，本文grid采用二維，block采用三維)

6行:調(diào)用并執(zhí)行單偶極子天線(xiàn)發(fā)出的輻射射線(xiàn)跟蹤函數(shù)kernel_1。

7行:釋放內(nèi)存、顯存上所有占用存儲(chǔ)單元，本次實(shí)驗(yàn)主要程序是前兩個(gè)kernel函數(shù)，所以下面主要介紹前兩個(gè)kernel的偽碼。

輻射射線(xiàn)跟蹤內(nèi)核函數(shù)偽碼部分如下。

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

4-6行:設(shè)置index輻射射線(xiàn)的方位角與天頂角，并由θ和φ惟一確定了射線(xiàn);

7-14:跟蹤射線(xiàn)，若發(fā)生反射則求入射角、反射系數(shù)、反射向量，經(jīng)過(guò)的路徑長(zhǎng)度及損耗，將反射點(diǎn)作為新的次級(jí)輻射源并繼續(xù)跟蹤，若發(fā)生繞射則存儲(chǔ)繞射點(diǎn)及射線(xiàn)與障礙物邊緣的夾角等參數(shù)，若射線(xiàn)被接收則存儲(chǔ)該射線(xiàn)并結(jié)束，這部分與傳統(tǒng)的射線(xiàn)跟蹤一致。

繞射射線(xiàn)跟蹤內(nèi)核函數(shù)偽碼部分如下。

判定該射線(xiàn)的傳播路徑，若反射則繼續(xù)跟蹤，若被接收則計(jì)算場(chǎng)強(qiáng)，存儲(chǔ)該射線(xiàn)信息并結(jié)束。

3 實(shí)驗(yàn)結(jié)果與分析

與傳統(tǒng)的經(jīng)驗(yàn)?zāi)Ｐ筒煌?，利用射線(xiàn)跟蹤技術(shù)對(duì)小區(qū)電波場(chǎng)強(qiáng)預(yù)測(cè)時(shí)需要小區(qū)地理、地形數(shù)據(jù)的支持。地理數(shù)據(jù)包括建筑物、道路、樹(shù)木、路燈、汽車(chē)、草坪等，這些施都會(huì)影響射線(xiàn)跟蹤的傳播路徑。如果完全貼近現(xiàn)實(shí)世界，那樣計(jì)算量無(wú)法估計(jì)，所以對(duì)三維世界進(jìn)行必要的簡(jiǎn)化是合理的。我們主要關(guān)注建筑物與道路對(duì)電波傳播的影響。其實(shí)電波傳播的實(shí)質(zhì)就是電磁波與三維空間中各種介質(zhì)相互作用的過(guò)程。這種相互作用從物理學(xué)的角度講，包括直射、反射、透射、繞射與散射等。本文采用長(zhǎng)方體包圍盒來(lái)簡(jiǎn)化建筑物的各種形態(tài)，用平行四邊形來(lái)表示地面。

為了給人以直觀(guān)的印象，利用標(biāo)準(zhǔn)圖形OpenGL庫(kù)進(jìn)行可視化顯示［11］。考慮計(jì)算量與計(jì)算機(jī)硬件配置，我們使用標(biāo)準(zhǔn)的多邊形與多面體來(lái)描繪三維空間地理數(shù)據(jù)。先在視口中放置一個(gè)長(zhǎng)方形(其規(guī)格為100 m×100 m)作為地面，然后在地面上放置4個(gè)建筑物，在正視圖當(dāng)中按逆時(shí)針?biāo)闫?，分別為標(biāo)號(hào)為1-4。建筑物的幾何結(jié)構(gòu)分別為40 m×30 m×90 m，20 m ×20 m ×60 m，30 m ×20 m ×50 m，60 m×40 m×40 m。簡(jiǎn)單三維地理數(shù)據(jù)的正視圖如圖5所示。另外，設(shè)計(jì)了一個(gè)天線(xiàn)放在建筑物1的頂上，其具體的坐標(biāo)位置為(30，90，0)。

圖5 簡(jiǎn)單三維地理數(shù)據(jù)的正視圖Fig.5 Front view of a simple 3 dimensional geography

以上介紹了三維建筑物的幾何結(jié)構(gòu)參數(shù)，下面將介紹與電波傳播有關(guān)的電磁參數(shù)，本次實(shí)驗(yàn)采用的天線(xiàn)為單偶極子天線(xiàn)，其輻射功率為43 dBm，使用頻段為1 800 MHz，極化方式為垂直極化，另外建筑物的相對(duì)介電常數(shù)εr=5.5，建筑物電導(dǎo)率σ=0.001，地面的相對(duì)介電常數(shù) εr=3.4，地面的電導(dǎo)率為σ=0.002。為了方便觀(guān)察，通過(guò)編寫(xiě)程序?qū)崿F(xiàn)自由控制整個(gè)視口的旋轉(zhuǎn)與縮放，從而可以從不同的視角來(lái)觀(guān)察三維空間地理模型。經(jīng)過(guò)旋轉(zhuǎn)處理，得到的俯視圖如圖6所示。

利用前面的算法編寫(xiě)并行程序進(jìn)行射線(xiàn)跟蹤與CUDA加速之后，信號(hào)強(qiáng)度的分布圖如圖7所示。

從7圖中可以看出，與天線(xiàn)處于視距，沒(méi)有障礙物阻擋時(shí)，信號(hào)能夠直達(dá)接收點(diǎn)，信號(hào)強(qiáng)度好，當(dāng)電波被建筑物阻擋時(shí)，相對(duì)于天線(xiàn)不可見(jiàn)部分的信號(hào)強(qiáng)度明顯較弱。這也類(lèi)似于光學(xué)里面的陰影效應(yīng)。實(shí)際當(dāng)中，繞射現(xiàn)象是存在的，但是本次實(shí)驗(yàn)結(jié)果顯示繞射現(xiàn)象不明顯，這個(gè)跟繞射算法有關(guān)系。

圖6 簡(jiǎn)單的三維地理數(shù)據(jù)的俯視圖Fig.6 Vertical view of a simple 3 dimensional geography

圖7 射線(xiàn)跟蹤之后分布Fig.7 Electronic strength distribution after ray tracing

串行與并行算法耗時(shí)對(duì)比圖如圖8所示。

圖8 串行與并行算法耗時(shí)對(duì)比圖Fig.8 Comparison chart of computing time between serial and parallel algorithm

從圖8中可以看出，當(dāng)跟蹤128×64條射線(xiàn)時(shí)，串行射線(xiàn)跟蹤與并行射線(xiàn)跟蹤所耗時(shí)相差不多，基本持平，甚至串行還比并行稍微用時(shí)少一點(diǎn)，這是因?yàn)镃UDA中host與device之間數(shù)據(jù)的傳輸消耗了很多時(shí)間，對(duì)于跟蹤較少射線(xiàn)時(shí)，并行算法的優(yōu)勢(shì)得不到體現(xiàn)，此時(shí)的效率比=并行耗時(shí)/串得耗時(shí)=1.07。當(dāng)跟蹤256×128條射線(xiàn)以及512×256條射線(xiàn)時(shí)，并行射線(xiàn)跟蹤明顯優(yōu)于串行射線(xiàn)跟蹤，這時(shí)大量的射線(xiàn)跟蹤可以將host與device之間數(shù)據(jù)傳輸所消耗時(shí)間淹沒(méi)掉，這時(shí)計(jì)算時(shí)間起主導(dǎo)作用，此時(shí)的效率比分別為:0.356，0.458，這說(shuō)明基于 GPU 的通用計(jì)算需要良好的分配，才能達(dá)到較好的提速效果。

4 結(jié)論

CUDA技術(shù)是一個(gè)新生的并行計(jì)算技術(shù)，相對(duì)于傳統(tǒng)的多核、Cluster集群并行計(jì)算更加優(yōu)越。由實(shí)驗(yàn)結(jié)果可證明，相對(duì)于傳統(tǒng)的全CPU編程，基于CUDA技術(shù)的無(wú)線(xiàn)電波傳播模擬預(yù)測(cè)可以大幅度地縮減計(jì)算耗時(shí)，提高計(jì)算效率。從而也改變了人們的編程思維方式，有利于并行計(jì)算的發(fā)展。

［1］楊超，徐江斌，吳玲達(dá).硬件加速的虛擬電磁環(huán)境體可視化［J］.北京郵電大學(xué)學(xué)報(bào)，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實(shí)踐［M］.聶學(xué)軍，譯.北京:機(jī)械工業(yè)出版社，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高性能運(yùn)算之CUDA［M］.北京:中國(guó)水利水電出版社，2009.ZHANG Shu.CUDA-high performance computing by GPU［M］.Beijing:China Water-Power Press，2009.

［9］郭建光.移動(dòng)通信系統(tǒng)中微小區(qū)射線(xiàn)跟蹤算法研究［M］.北京:北京郵電大學(xué)出版社，2010.GUO Jianguang.Microcellular Ray-tracing Algorithm in Mobile Communication System［M］.Beijing:Beijing U-niversity of Posts and Telecommunications Press，2010.

［10］袁正午.移動(dòng)通信系統(tǒng)終端射線(xiàn)跟蹤定位理論與方法［M］.北京:電子工業(yè)出版社，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版.徐波，譯.北京:機(jī)械工業(yè)出版社，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.