劉興霞，張利軍，趙玉祥

（天水師范學(xué)院電子信息與電氣工程學(xué)院，甘肅 天水 741001）

0 引言

彈跳射線（SBR）方法［1］是一種通用且高效的計(jì)算電大復(fù)雜目標(biāo)雷達(dá)散射截面的估算方法，其混合了幾何光學(xué)法和物理光學(xué)法；由于SBR 方法考慮到了目標(biāo)幾何結(jié)構(gòu)間的多次反射情況，因此特別適合計(jì)算此類(lèi)多次反射場(chǎng)的問(wèn)題。但是在使用SBR 方法計(jì)算電大尺寸目標(biāo)時(shí)，由于目標(biāo)幾何模型復(fù)雜，面元數(shù)量大，且考慮多次反射等因素，射線管數(shù)量十分巨大，射線管的追蹤和更新將十分耗時(shí)，運(yùn)算效率不高。本文采用一種層次包圍盒［2］（BVH）的二叉樹(shù)加速數(shù)據(jù)結(jié)構(gòu)［3］，提出了在圖形處理器（GPU）端實(shí)現(xiàn)射線管的分裂和追蹤更新的計(jì)算統(tǒng)一設(shè)備構(gòu)架（CUDA）并行計(jì)算SBR 算法.

1 BVH 的光線追蹤和基于GPU 的并行計(jì)算技術(shù)

SBR 方法與圖形學(xué)中的光線跟蹤算法的處理思路十分相似，主要包括射線管分裂，追蹤更新和電磁計(jì)算這三個(gè)步驟［4］。BVH 樹(shù)是光線跟蹤算法中最基本的加速技術(shù)之一。根據(jù)目標(biāo)幾何模型的面元的空間分布的情況，使用一個(gè)包圍盒將臨近的一個(gè)或多個(gè)面元包裹住，進(jìn)而通過(guò)構(gòu)造樹(shù)狀的層次結(jié)構(gòu)來(lái)將目標(biāo)幾何模型的面元用包圍盒完全包裹住。當(dāng)射線管與當(dāng)前的包圍盒不相交，那么射線管與該包圍盒里的三角面元或者其內(nèi)部的小包圍盒也不相交；當(dāng)線管與當(dāng)前的包圍盒相交，那么才進(jìn)行三角面元的求交計(jì)算或者往下遍歷其內(nèi)部的包圍盒。因此BVH 樹(shù)能夠顯著的減少每根射線管與三角面元求交的計(jì)算量，提高計(jì)算效率。而基于GPU 的CUDA 并行計(jì)算在處理大規(guī)模的相同的問(wèn)題時(shí)能夠顯著的提高計(jì)算效率［5］，對(duì)于這種射線管的數(shù)量巨大，并且單根射線管之間時(shí)相互獨(dú)立的，射線管的追蹤都是依照相同的步驟執(zhí)行的問(wèn)題模型，非常適合采用CUDA 處理，來(lái)提高其運(yùn)算速度［6］。

2 算法的設(shè)計(jì)與實(shí)現(xiàn)

2.1 基于GPU 的SBR方法的實(shí)現(xiàn)過(guò)程

在本文中，加速樹(shù)的構(gòu)建在CPU 端完成，將其轉(zhuǎn)化為線性數(shù)組后拷貝至GPU 上，射線管的生成，追蹤，RCS求解都是在GPU 端上完成的。

基于GPU 的SBR 方法與傳統(tǒng)的SBR 方法的流程是一樣的，主要區(qū)別在于在本文中是在GPU上實(shí)現(xiàn)了射線管的分裂，追蹤更新和電磁計(jì)算。在CPU 端構(gòu)建所需的BVH 樹(shù)結(jié)構(gòu)，主要是通過(guò)SAH算法將面元進(jìn)行劃分，具體的建樹(shù)方法在文獻(xiàn)［3］和文獻(xiàn)［7］中有詳細(xì)的描述，然后將樹(shù)結(jié)構(gòu)和其他需要的參數(shù)拷貝至GPU 端。SBR 方法的主要流程如下：在GPU 上生成并初始化射線管，如圖1 上所示；然后追蹤射線管，對(duì)于與目標(biāo)面有交點(diǎn)的射線管求交并更新其交點(diǎn)和電磁信息，對(duì)于與目標(biāo)面沒(méi)有交點(diǎn)的射線管則不計(jì)算，直到所有的射線管都離開(kāi)目標(biāo)表面，如圖1中間所示；最后對(duì)于發(fā)生更新的射線管計(jì)算其雷達(dá)散射截面（RCS），如圖1下所示，具體的電磁計(jì)算公式在文獻(xiàn)［8］中有詳細(xì)的討論。這里，主要討論帶有加速樹(shù)的SBR 方法在GPU 上的執(zhí)行步驟。

2.2 射線管任務(wù)在GPU 上的劃分

由于射線管是二維分布的，在GPU 上采用二維的線程塊和網(wǎng)格來(lái)對(duì)應(yīng)處理射線管，每個(gè)線程處理一根射線管，若干個(gè)線程組成一個(gè)線程塊BLOCK，然后若干個(gè)線程塊組成一個(gè)網(wǎng)格GRID。因此可以將射線管的追蹤計(jì)算完整的映射到GPU的計(jì)算模型上。計(jì)算任務(wù)的尺寸是有SBR 方法分裂的射線管數(shù)量決定的；BLOCK 的大小取決于GPU 硬件的限制，由程序員來(lái)優(yōu)化，一般來(lái)說(shuō)為了更好的利用GPU 資源，應(yīng)該盡量讓BLOCK 的大小是32的整數(shù)倍［9］。一旦確定了計(jì)算問(wèn)題尺寸和BLOCK 的尺寸，就可以計(jì)算出GRID 的尺寸。GRID 在某個(gè)維度上的尺寸大?。ㄒ詘 軸為例）：DRID 在x 軸上的尺寸= ｛問(wèn) 題在x 軸上的尺寸+每個(gè)BLOCK 在x 軸上的尺寸- ｝1 ÷每個(gè)BLOCK 在x 軸上的尺寸

采用這種確定GRID 的方法可以保證待求問(wèn)題都會(huì)處理，不會(huì)因?yàn)檎麛?shù)求除運(yùn)算只取整數(shù)部分而出現(xiàn)待求問(wèn)題的邊界得不到處理的異常。但是這樣會(huì)導(dǎo)致實(shí)際的線程數(shù)大于所需的線程數(shù)，在內(nèi)核的處理中，加入邊界的判斷，對(duì)于超出邊界的線程不參與計(jì)算即可。對(duì)于無(wú)法一次加載所有計(jì)算數(shù)據(jù)的任務(wù)，可將總的射線管追蹤任務(wù)分解成更細(xì)小的塊來(lái)處理，每次GPU 只會(huì)加載其中的一塊任務(wù)進(jìn)行處理，處理完畢后，求出該塊上的總的RCS并將其寫(xiě)入顯存中，然后釋放顯卡資源；依次處理完所有的塊并將所有塊的總的RCS 求和，即可得到最終所需結(jié)果。

2.3 包圍盒的線索添加

二維的BVH 樹(shù)如圖2 左所示。BVH 樹(shù)的節(jié)點(diǎn)用線索來(lái)增強(qiáng)，每個(gè)樹(shù)的節(jié)點(diǎn)的線索指向該節(jié)點(diǎn)的相鄰的樹(shù)節(jié)點(diǎn)。樹(shù)的線索是在樹(shù)的構(gòu)建中遞歸生成的［9］，其流程如下：對(duì)一個(gè)樹(shù)的葉節(jié)點(diǎn)分割后，該節(jié)點(diǎn)的兩個(gè)孩子節(jié)點(diǎn)繼承其父節(jié)點(diǎn)的線索的值，然后子節(jié)點(diǎn)在分割面上的線索相互指向另一子節(jié)點(diǎn)。射線管在BVH 樹(shù)遍歷時(shí)，可以在需要時(shí)通過(guò)當(dāng)前遍歷的節(jié)點(diǎn)的線索直接訪問(wèn)與其空間相鄰的節(jié)點(diǎn)，因此消除了待遍歷節(jié)點(diǎn)的堆棧，不僅減少了堆棧數(shù)據(jù)的壓入彈出，而且在計(jì)算多次反射時(shí)有效的減少了不必要的節(jié)點(diǎn)遍歷。如圖2右所示，對(duì)于二維的BVH 樹(shù)結(jié)構(gòu)，節(jié)點(diǎn)6有4個(gè)線索分別指向節(jié)點(diǎn)2和7，還有2個(gè)線索為空。對(duì)于三維的BVH 樹(shù)，共有6個(gè)線索，分別指向包圍盒的6個(gè)側(cè)面。

圖2 二維BVH 樹(shù)結(jié)構(gòu)Fig.2 Two dimension BVH tree structure

2.4 BVH 樹(shù)在GPU 上的無(wú)堆棧遍歷

由于BVH 樹(shù)的遍歷要在GPU 上運(yùn)行，因此必須將BVH 樹(shù)的數(shù)據(jù)傳遞到GPU 端。使用按層次遍歷樹(shù)的方法來(lái)將其先在CPU 端轉(zhuǎn)化為線性數(shù)組，然后在拷貝到GPU 上。將樹(shù)的節(jié)點(diǎn)數(shù)據(jù)存儲(chǔ)為線性數(shù)組不僅便于將樹(shù)的數(shù)據(jù)從CPU 端拷貝到GPU 端，而且由于GPU 在并行處理時(shí)是按線程束來(lái)讀寫(xiě)和處理數(shù)據(jù)的，當(dāng)樹(shù)的數(shù)據(jù)在內(nèi)存中是連續(xù)排列時(shí)，對(duì)于連續(xù)的待讀取數(shù)據(jù)在滿(mǎn)足一定條件下可以將幾次的數(shù)據(jù)讀取操作合并為一次數(shù)據(jù)讀取操作，從而提高了數(shù)據(jù)的讀取速度。轉(zhuǎn)化后的樹(shù)結(jié)構(gòu)在GPU 端的存儲(chǔ)方式如圖3 所示，數(shù)據(jù)都被存儲(chǔ)到一塊連續(xù)的顯存空間中。此時(shí)使用節(jié)點(diǎn)編號(hào)來(lái)替代地址來(lái)尋找節(jié)點(diǎn)信息，并且將線索中指向其他節(jié)點(diǎn)地址的指針替換為指向其他節(jié)點(diǎn)的編號(hào)。

BVH 樹(shù)在GPU 上無(wú)堆棧遍歷過(guò)程如圖3 所示，CPU 與GPU 上樹(shù)的無(wú)堆棧遍歷的主要區(qū)別在于CPU 上樹(shù)的節(jié)點(diǎn)的數(shù)據(jù)在內(nèi)存中是隨機(jī)分布的，在GPU 上樹(shù)的節(jié)點(diǎn)的數(shù)據(jù)在顯存中是連續(xù)的。

圖3 二維BVH 在GPU 上無(wú)堆棧遍歷Fig.3 Two dimension BVH on GPU without a stack traversal

以二維BVH 樹(shù)為例，對(duì)于如圖2 左所示的樹(shù)結(jié)構(gòu)，射線從上方入射到該樹(shù)結(jié)構(gòu)，最后與節(jié)點(diǎn)9中的三角面元相交，簡(jiǎn)要說(shuō)明其遍歷過(guò)程如下：射線從根節(jié)點(diǎn)1開(kāi)始遍歷，接著遍歷中間節(jié)點(diǎn)3和葉子節(jié)點(diǎn)6；與葉子節(jié)點(diǎn)6中的三角面元求交，無(wú)交點(diǎn)，此時(shí)判斷射線的出射面，通過(guò)該節(jié)點(diǎn)的線索來(lái)到并遍歷中間節(jié)點(diǎn)2，接著遍歷中間節(jié)點(diǎn)4和葉子節(jié)點(diǎn)9，在葉子節(jié)點(diǎn)9處射線與其內(nèi)部的三角面元求交，存在交點(diǎn)，并更新該射線管相關(guān)信息。該射線追蹤的一次追蹤完畢。

3 實(shí)驗(yàn)驗(yàn)證及結(jié)果討論

為了驗(yàn)證算法的運(yùn)算效率，本文做如下仿真驗(yàn)證。仿真所使用的計(jì)算機(jī)硬件為Intel i5 處理器，4 GB內(nèi)存，顯卡為GTX650（流處理器384個(gè)），1GB顯存。射線管分裂步長(zhǎng)為1／10個(gè)波長(zhǎng)。為了驗(yàn)證結(jié)果的正確性，使用FEKO 軟件的MLFMM 方法的計(jì)算結(jié)果與本文方法的計(jì)算結(jié)果來(lái)作對(duì)比。為了簡(jiǎn)化說(shuō)明，將在CPU 上運(yùn)行的原始SBR 方法簡(jiǎn)稱(chēng)為SBR，將在CPU 上運(yùn)行的帶BVH 加速樹(shù)的SBR方法簡(jiǎn)稱(chēng)為CPU-SBR，將在GPU 上運(yùn)行的帶有BVH 加速樹(shù)的SBR 方法簡(jiǎn)稱(chēng)為GPU-SBR。

具體參數(shù)為：三面角的參數(shù)長(zhǎng)寬高均為1m，入射電磁波的頻率為3GHz，HH 極化。

圖4左中的角度信息為入射角θ=45°，方位角φ 從0°到90°；圖4右中的角度信息為方位角φ=45°，入射角θ從0°到90°?？梢钥闯龆鄬涌焖俣鄻O子方法（MLFMM）的計(jì)算結(jié)果和CPU-SBR 與GPU-SBR 方法的計(jì)算結(jié)果符合的很好，該方法能夠適用于多次散射的電磁計(jì)算問(wèn)題。

圖4 三面角的RCSFig.4 RCS of the three corners

為了驗(yàn)證本文方法對(duì)復(fù)雜目標(biāo)散射場(chǎng)計(jì)算的有效性，計(jì)算了一個(gè)簡(jiǎn)化的船模型，如圖5所示，模型的長(zhǎng)為1 m，寬0.2 m，高0.2 m，入射波頻率為5 GHz，HH 極化。

圖6中的角度信息為方位角φ=0°，入射角θ從-90°到225°。從-90°到90°為簡(jiǎn)化船的頂部，MLFMM 的結(jié)果符合的很好，從90°到225°為船的底部，由于是平面結(jié)構(gòu)，使用SBR 法計(jì)算此類(lèi)情況時(shí)有一定的計(jì)算差異。

表1 簡(jiǎn)化船模型不同的方法的計(jì)算時(shí)間（秒）比較Tab.1 Comparison of computational time（s）for different methods of ship model

表1中，將不同方法計(jì)算消耗的總時(shí)間除以所有入射角度個(gè)數(shù)得到了單個(gè)角度下不同計(jì)算方法消耗的平均時(shí)間。從表1 可以看到，GPU-SBR 方法比CPU-SBR 方 法 快 了10 倍 左 右，GPU-SBR 方 法比SBR 方法快了198倍左右。本文方法在求解目標(biāo)的RCS時(shí)在保證準(zhǔn)確性的同時(shí)，還能夠獲得較好的加速效果。

為了進(jìn)一步驗(yàn)證本文方法的有效性，我們還計(jì)算了一個(gè)甲板上主船體部分模型的單站RCS。其幾何結(jié)構(gòu)如圖7所示，其尺寸為長(zhǎng)110m，寬26m，高44.4m。電磁波頻率為3GHz，HH 極化。

圖5 簡(jiǎn)化的船模型Fig.5 Simplified ship model

圖6 簡(jiǎn)化船模型的RCSFig.6 The RCS of Simplify the ship model

圖7 上半部分復(fù)雜艦船模型示意圖Fig.7 Schematic diagram of the complex ship model of the upper part

圖8（a）中的方位角φ=0°，入射角θ為從-90°到90°。圖8（b）的入射角θ=450，方位角φ 從-180°到180°。其中紅線表示的為本文計(jì)算方法得到艦船目標(biāo)的RCS結(jié)果，黑線表示的為只考慮一次散射的物理光學(xué)法（PO）法求解的RCS。通過(guò)計(jì)算可以看出，本文計(jì)算方法只需要786s，而傳統(tǒng)的SBR 方法需要一天以上的時(shí)間，這有效地說(shuō)明本文方法對(duì)于電大尺寸復(fù)雜目標(biāo)的電磁散射計(jì)算也有良好的擴(kuò)展性。另外，相對(duì)于PO 法計(jì)算效果和精度，采用本文SBR 方法做三次射線追蹤即可達(dá)到滿(mǎn)意的效果。

由于在計(jì)算該問(wèn)題時(shí)所消耗的時(shí)間很長(zhǎng)，只選取在入射角θ=45°，方位角φ=0°時(shí)這一個(gè)角度的計(jì)算時(shí)間作比較。

從表2可以看到，基于GPU 的CUDA 帶BVH加速樹(shù)的SBR 方法比在CPU 上使用BVH 樹(shù)的SBR 方法快了46倍左右。隨著目標(biāo)模型的尺寸和結(jié)構(gòu)復(fù)雜度的增加，其加速比也隨之增加。這有效地說(shuō)明本文方法對(duì)于電大尺寸復(fù)雜目標(biāo)的電磁散射計(jì)算也有良好的擴(kuò)展性。

圖8 不同入射角度下復(fù)雜艦船的RCSFig.8 The RCS of a complex ship with different incident angles

表2 復(fù)雜船模型不同的方法的計(jì)算時(shí)間比較Tab.2 Comparison of computational time for different methods of complex ship model

4 結(jié)論

本文改進(jìn)了電大目標(biāo)電磁散射彈跳射線算法，在采用BVH 樹(shù)的二叉樹(shù)加速數(shù)據(jù)結(jié)構(gòu)的基礎(chǔ)上，提出了在GPU 實(shí)現(xiàn)射線管的分裂和追蹤更新的CUDA 并行計(jì)算SBR 算法。通過(guò)實(shí)驗(yàn)算例驗(yàn)證表明，該方法可有效提高SBR 方法處理電大尺寸目標(biāo)電磁散射問(wèn)題的的計(jì)算效率，并且隨著GPU 性能的提升，計(jì)算效率還有進(jìn)一步提升的空間；同時(shí)，該方法也為傳統(tǒng)SBR 方法更新發(fā)展提供一種新的思路。

［1］Ding Jianjun，Chen Ru shan，Zhou Huan et al.An improverment for the acceleration technique based on monostatic bistatic equivalence for shooting and bouncing ray method［J］.Microwave and optical technology letters，2011，53（5）：259-266.

［2］張曉莉，王苗.數(shù)據(jù)結(jié)構(gòu)與算法［M］.北京：機(jī)械工業(yè)出版社，2008：97-141.

［3］Pharr M，Humphreys G.Physically based rendering：From theory to implementation ［M］.Morgan Kauf-mann，2004：208-245.

［4］Hermann Buddendick，Thomas F Eibert.Acceleration of Ray-Based Radar Cross Section Predictions Using Monostatic-Bistatic Equivalence［J］.IEEE Transactions on Antennas and Propagation 2010，58（2）：320-329.

［5］Tao Y，Lin H，Bao H.GPU-Based Shooting and Bouncing Ray Method for Fast RCS Prediction ［J］.IEEE Transactions on Antennas and Propagation，2010，58（2）：356-362.

［6］張舒，褚艷利，趙開(kāi)勇，張鈺勃.GPU 高效能運(yùn)算之CUDA［M］.北京：中國(guó)水利水電出版社，2009.

［7］Macdonald J D，Booth K S.Heuristics for ray tracing using space subdivision［J］.The Visual Computer，1990，6（3）：153-166.

［8］李軍，馬東立，武哲.用射線法研究復(fù)雜目標(biāo)的電磁散射特 征 ［J］.北 京 航 空 航 天 大 學(xué) 學(xué) 報(bào)，2010，26（5）：573-576.

［9］Popov S，Günther J，Seidel H P，et al.Stackless KD-tree traversal for high performance GPU ray tracing［J］.EURO Graphics，2007：1278-1283.