劉興霞，張利軍，趙玉祥

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

0 引言

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

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

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

2 算法的設計與實現(xiàn)

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

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

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

2.2 射線管任務在GPU 上的劃分

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

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

2.3 包圍盒的線索添加

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

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

2.4 BVH 樹在GPU 上的無堆棧遍歷

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

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

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

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

3 實驗驗證及結(jié)果討論

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

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

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

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

為了驗證本文方法對復雜目標散射場計算的有效性，計算了一個簡化的船模型，如圖5所示，模型的長為1 m，寬0.2 m，高0.2 m，入射波頻率為5 GHz，HH 極化。

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

表1 簡化船模型不同的方法的計算時間（秒）比較Tab.1 Comparison of computational time（s）for different methods of ship model

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

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

圖5 簡化的船模型Fig.5 Simplified ship model

圖6 簡化船模型的RCSFig.6 The RCS of Simplify the ship model

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

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

由于在計算該問題時所消耗的時間很長，只選取在入射角θ=45°，方位角φ=0°時這一個角度的計算時間作比較。

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

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

表2 復雜船模型不同的方法的計算時間比較Tab.2 Comparison of computational time for different methods of complex ship model

4 結(jié)論

本文改進了電大目標電磁散射彈跳射線算法，在采用BVH 樹的二叉樹加速數(shù)據(jù)結(jié)構(gòu)的基礎上，提出了在GPU 實現(xiàn)射線管的分裂和追蹤更新的CUDA 并行計算SBR 算法。通過實驗算例驗證表明，該方法可有效提高SBR 方法處理電大尺寸目標電磁散射問題的的計算效率，并且隨著GPU 性能的提升，計算效率還有進一步提升的空間；同時，該方法也為傳統(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］.北京：機械工業(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］張舒，褚艷利，趙開勇，張鈺勃.GPU 高效能運算之CUDA［M］.北京：中國水利水電出版社，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］李軍，馬東立，武哲.用射線法研究復雜目標的電磁散射特 征 ［J］.北 京 航 空 航 天 大 學 學 報，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.