基于室內動靜結合分割的射線跟蹤加速方法

黃一航，江 虹，韓 賓

（西南科技大學信息工程學院，四川綿陽 621010）

（*通信作者電子郵箱491938847@qq.com）

0 引言

信息通信技術的快速發展，特別是5G時代帶來了通信頻率的提高、室內基站以及智能天線等技術的應用，使得對預測電磁波傳播的精確性有了更高的要求。

射線跟蹤（Ray Tracing，RT）算法是基于幾何光學理論、幾何繞射理論和電磁理論的確定性預測模型［1］，廣泛應用于電磁傳播仿真中。傳統的射線跟蹤算法在復雜的密閉三維環境模型下應用［2］能夠提高算法精確性，但是成倍增加的無用求交次數會導致計算效率大幅度降低［3］。因此，業界一直在尋求如何大幅度降低無用交點的求解次數，以在保證射線跟蹤算法精確性的同時能提高運算的效率。

1 射線跟蹤算法工作機制

射線跟蹤算法是通過模擬電磁波的傳播路徑，來確定多徑信道在發射與接收機之間傳播過程中所有可能存在的路徑［4］。基本思想是將從源點輻射出的電磁波模擬為一條條射線，電磁在各自獨立的射線管［5］內傳播。依次沿著各條射線路徑進行跟蹤，判斷射線是否碰到物體表面或者被接收點所接收。若射線遇到障礙物會發生反射、透射、繞射和散射現象，根據幾何光學理論和一致性繞射原理［6］計算出各空間傳播機制發生后的場強。跟蹤射線直到場強衰減到預設閾值［7］，停止對本條射線跟蹤。求得到達接收點的所有射線后，使用矢量疊加的方法計算出輻射源傳遞到目標位置的信號值。

在目標接收點處的總場強（單位：V/m）可表示為：

其中：接收點處矢量合并后的總輻射場強用Etotal表示；l為射線路徑總條數；Ei為第i條射線路徑末場的矢量場強。矢量場強Ei由式（2）［8］計算：

其中：Einc是射線在第一個反射、繞射或透射節點處的矢量輻射場；n為出現反射的總次數；Rh是在第h次反射時的并矢反射系數；m為出現繞射的總次數；Df是在第f次繞射時的并矢繞射系數；u為出現透射的總次數；Tt是在第t次透射時的并矢透射系數；As是經過反射透射或繞射后的擴散因子；rq是跟蹤的射線的第q個傳播節點到第q+1 個節點的距離，q為跟蹤的信道總的空間傳播機制節點數。

計算出終點場強Etotal后，接收點處的接收功率Pr可由式（3）［9］計算：

其中：Pt為發射天線輻射功率；Gt和Gr分別為發射天線和接收天線增益；λ為由電磁波頻率計算出的發射天線工作波長；E0為發射點處初始電場強度。

由以上理論分析可知，影響射線跟蹤算法預測精度和計算效率的因素主要有以下方面：

1）射線源發射的射線數量和射線之間的角度，決定了算法的預測精度。理論上發射間隔越小、跟蹤的射線數量越多以及射線之間分布得越均勻，模擬電磁波傳播的仿真算法就會越精確。本文使用相同大小正方形無縫覆蓋單位球面的方式建立發射源［10］。

其中：（Vx，Vy，Vz）是需要跟蹤的射線方向向量vt；Δθ是相鄰射線間的角度，遍歷n和m值即可得到所有射線的方向。因此間隔角度直接決定了射線源發射的射線數量。

因為跟蹤一條完整的射線需要多次遍歷所有的三維環境數據，所以當跟蹤的射線數量越多，功率預測值雖然會更精確，但也將導致算法運行時長大幅增加。

2）空間內物體的數量與環境復雜程度決定了求交點運算次數的多寡。跟蹤某一射線的基本步驟，是將該射線與檢索范圍內所有的物體表面進行求交點計算，找到交點坐標后再計算出相應的反射向量。因此跟蹤的射線為了找到一個正確的碰撞點需要進行多次求交計算以預先找到多個交點坐標，而在這些交點中只有距起點Pstart最近的碰撞點是有效的碰撞點Pc和反射向量。即信號在傳播過程中第一個碰撞到的障礙物交點是有效的，其余的求交點計算都是無用的運算。有效交點坐標Pc可以由式（5）計算得到：

其中：D是根據相似三角原理計算得到的射線起點Pstart到交點之間的距離；v是射線起點與三維面內任意一點的向量；vt是跟蹤射線的方向向量；n是三維面的法向量。

反射向量vRe由式（6）計算：

在整個射線跟蹤算法運算中，只有能被接收點接收，且在射線傳播路徑中碰撞到正確三維表面的求交點計算才是有效計算，其余的交點計算都是無用運算。所以大量無用求交點運算導致了較低的計算效率。

在上述影響因素中，一般在達到精度要求后射線的間隔角大小是固定不變的，即發射源發射的射線數量通常不變。因此環境復雜程度（障礙物分布情況）主要決定了求交運算的次數。而隨著環境復雜度的提高，無用求交運算次數所占的比例會大幅度增加。所以，減少算法中大量的無用求交運算是提高計算效率的主要途徑。

2 基于三維空間分割的加速方法

在傳統的射線跟蹤算法中，正在跟蹤的射線向量需要與環境模型內所有物體的三維表面進行求交點的運算［11］。隨著環境復雜度的提高，繞射發生次數隨之增加，使得算法需要跟蹤更多的射線。并且環境復雜度和總射線數量的提高都會使得算法求交運算的次數大幅提升。

為此，文獻［12］中提出了一種基于動態分區的射線跟蹤加速方法。該方法根據環境模型內的障礙物分布情況動態劃分區域，可減少射線與建筑物求交點次數；但是該方法只考慮了二維模型，不適合三維環境模型。因此，文獻［13］中提出了使用正向發射與反向鏡像混合的方法對原始射線跟蹤算法進行改進；但該加速方法只是單純改進了算法的計算方式，沒有考慮到環境復雜度提高對算法計算效率的負增益。

因此，本文提出了基于靜態空間分割與動態空間分割相結合的方法對三維環境模型進行空間分割，旨在大幅減少求交運算次數，有效提高算法的計算效率。

2.1 靜態三維空間分割的加速方法

對三維環境模型進行靜態空間分割，是將一個大區域平均劃分為多個小分區以分攤區域中物體的數量，形成n級區域分割方式：一級分割將大區域平均劃分為四部分；二級分割是將已經分割后的區域再次平均劃分；后續的高級分割以此類推。射線跟蹤每次只在射線當前所在分區Nold中檢索碰撞點，若在該分區沒有碰撞點，則射線離開當前分區Nold進入到下一個分區Nnew，并在Nnew內繼續檢索碰撞點，直到跟蹤的射線到達接收位置或信號衰落權重超過預設閾值為止。如圖1所示為空間分割方法流程。

對大區域進行空間分割后，檢索碰撞點的算法將不需要和大區域內所有的物體表面進行求交點運算，只需在當前射線所在的小分區中進行碰撞檢索，這樣可以大幅度減少跟蹤一條射線所需檢索的碰撞點的次數。

因此，正確判斷射線當前所在的分區，是靜態空間分割方法的關鍵。如式（7）所示：

其中：Nold為射線當前所在的分區序號；Nnew為射線將要進入的分區序號；n為空間分割等級；Pc為射線碰撞點坐標；x和y是射線三維方向向量中的x值與y值。

圖1 空間分割方法流程Fig.1 Method flow of space segmentation algorithm

當射線的碰撞點Pc在分區的邊界上時，可根據式（7）計算得出射線下一次范圍檢索的正確分區序號。

理論上對空間分割得越細，求交點運算計算量減少得越多，算法的效率也會越高。但是每細分一次區域不可避免會增加8 條分區的邊界面，且一個三維物體通常為6 個表面，因此，當一個小分區中包含物體的表面數量少于等于6 條時就不能再細分區域。

圖2（a）是沒有進行空間分割加速的無線傳播三維建模俯視圖，使用的是原始射線跟蹤算法對傳播過程進行跟蹤仿真。內部矩形部分是三維環境中各物體的三維模型俯視圖，內部線條是跟蹤算法跟蹤到的能從信號發射端到接收端的無線三維信道模型，該模型可以作為后續對改進加速方法模型精確度的對比驗證模型。圖2（b）是靜態一級空間分割的三維俯視效果圖，使用一級射線跟蹤加速方法對無線信道進行跟蹤。在同一個三維環境模型下和圖2（a）的原始跟蹤算法作對比驗證，可以看出兩張俯視效果圖中無線信道三維模型的數量和布局都一致，所以各級的靜態加速方法沒有降低算法的精確性。

圖2 環境空間分割俯視圖Fig.2 Vertical views of environmental space division

2.2 動態結合靜態的三維空間分割的加速方法

在三維環境下，雖然單獨使用靜態空間分割可以大幅提高計算效率。但在某些特殊情況下，靜態分割方法提升算法效率的作用會失效。如圖3 所示，當跟蹤射線的相鄰兩個碰撞點間橫跨了多個小分區時，由于算法會遍歷橫跨的各小分區中所有的物體面，所以無用求交運算的比例將逐漸接近傳統的射線跟蹤算法，使得在該情況下算法的計算效率并沒有獲得提高。

圖3 射線一次性穿越多個小分區情況Fig.3 Ray traversing multiple small divisions at once

因為射線跟蹤是在三維環境中進行，所以上述特殊情況發生的原因大概率是跟蹤射線的該段傳播高度過高，超越了其所在小分區物體的最高高度，不能在小分區內形成碰撞點。所以兩個相鄰碰撞點之間可能將跨越多個分區。

為了改善以上問題，本文在靜態空間分割的基礎上設計了一種基于高度劃分的動靜結合分割方法。實際密閉環境下的障礙物都會有一個最高高度，所以可以根據如下方式設計動靜結合的空間分割，以進一步提高算法效率：

1）原始三維環境區域內所有物體高度的集合為h，則根據式（8）動態空間分割的閾值hd應是該區域內最高物體的高度值。

2）假設密閉空間高度為H，將整個三維環境模型根據分割閾值hd動態分割為上下兩部分：上半部分空間合并劃分為一整塊空白區域；下半部分空間依舊使用靜態空間分割的方式劃分區域。

3）與靜態空間分割類似，實現該方法需要正確識別跟蹤的射線當前所在的分區序號。如式（9）所示：

其中：Nnew為射線將要進入的分區序號；為射線碰撞點所在的高度；z為射線三維方向向量中高度值；u為上半部分分區序號；d為下半分區中某靜態小分區序號，根據后續遍歷法判斷射線具體進入的小分區。

動靜結合的分割方法可以最大限度地提高算法的計算效率。當只有靜態空間分割時，跟蹤某些角度的射線時，檢索一個碰撞點可能一次會穿越多個小分區，如此造成的無用求交檢索也會較多，導致計算效率大幅降低。因此本文使用對高度的動態分割方法結合靜態的空間分割改進了算法的不足。如圖4 所示，上升的射線進入動態分割的上半分區，只需檢索一個上半分區便可替代原本靜態方法需要檢索的多個小分區，找到下一個碰撞點位置，大大減少了無用面與射線的求交運算次數。在圖4 的三維環境模型下，此加速方法的算法效率在靜態法已經提升的基礎上還會有9.8%左右的提高。

圖4 靜態與動態空間分割Fig.4 Static and dynamic space divisions

2.3 分區實現

本文的空間分割方法使用了四叉樹的原理［14］。四叉樹是一棵重力平衡樹，通過對目標根節點空間進行平均四分，直至滿足預測條件停止，最終形成一棵層次分明的樹。空間分割步驟如下：

1）對真實環境建立模型后，獲取環境模型的二維投影數據，通過坐標平移，使投影數據中最小位置點置于坐標原點，形成初始大區域。

2）檢索環境模型中所有物體的高度確定最大值，該值作為動態分割的動態閾值。

3）對動態分割后的下半部分區域做靜態空間分割。初始平均劃分為2×2排列的4個子分區并排序（如圖5所示），然后統計各分區中物體數量。

4）若子分區內的物體總數大于預設閾值，則該子分區返回步驟3）繼續進行靜態分割，直到所有分區內物體數量小于閾值。

5）動態與靜態空間分割完成后，判斷環境中每個物體分別所屬的小分區序號，采用遍歷法判斷。

上述空間分割方法實現的重要一步，是確定空間內各物體所屬小分區，一旦出現位置誤判，最終計算所得的仿真結果將會和真實數據相去甚遠。本文使用遍歷法來精確判斷物體所在的小分區：物體棱邊的平面線段與分割邊界做交點計算判斷交點所在分區，若無交點則判斷線段在分割邊界的哪一側，進而確定物體所在小分區。

以圖5 所示的一級靜態分割為例，使用式（10）判斷物體的棱邊線段所屬的小分區，進而判斷物體所在的分區序號。N為線段所在的分區序號，X為指向坐標x軸正方向的分割邊界單位向量，Y為指向坐標y軸正方向的分割邊界單位向量，L為兩分割線交點指向棱邊線段上任意一點的方向向量。

更高級分區的相關判斷以及式（9）中d序號的確立沿用上述公式的衍生公式。

圖5 靜態空間分割排序示意圖Fig.5 Schematic diagram of sorting for static space division

2.4 方法分析

2.4.1 時間復雜度

根據上述加速方法實現流程，分析得到加速方法跟蹤一條射線的計算時長T由以下幾個部分構成：1）時間T1，根據大環境中物體的分布情況動態分割區域，且將各物體的表面分配到各自所在小分區中的耗時；2）時間T2，信道射線跟蹤過程中與物體表面進行的求交點運算以及進一步計算射線入射、出射向量性質的耗時；3）時間T3，跟蹤的射線穿越邊界時判斷進入的下個分區的耗時；4）時間T4，判斷跟蹤的射線是否能被接收區域接收的耗時；5）時間T5，計算接收區域信道傳播特性的耗時。則可知：

然而，由于需要跟蹤的射線數量龐大，加速方法是一個多層循環運算，且T1和T5只需一次計算即可得到結果，所以T1和T5在時間復雜度中所占的比例極小，可以忽略不計。

假設大環境中總的物體表面個數為n、動態分割出了g個小分區，即可知各個小分區中平均擁有的物體表面個數為n/g。因為T2中主要的運算是射線遍歷所在小分區中所有物體表面進行求交點等相關計算，若跟蹤一條射線與一個物體表面做相關計算耗時為k1，那么分區中一次求交點計算的耗時為(n/g) ×k1。假設檢索到一個射線碰撞點平均需要遍歷v個小分區，且一條射線從開始到停止跟蹤平均需要檢索m個碰撞點，則時間T2為：

T3中的主要運算是當射線與分區中物體都沒有交點時，根據射線方向向量一次判斷進入的下一個分區。假設相關的操作耗時為k2，則時間T3為：

T4中的主要運算是每次檢索到碰撞點后，判斷射線是否能被接收區域所接收。假設相關的操作耗時為k3，則時間T4為：

相對于T3和T4，T2的耗時占了絕大部分時間，這首先是因為求交點等相關運算的計算過程復雜，而其他耗時部分的主要運算過程相對比較簡單，因此使得k1遠大于k2和k3；其次，則因為T2同時與大環境中總的物體表面個數n、發生一次空間傳播機制平均需要遍歷的分區個數v以及重復檢索空間傳播機制的次數m都相關，且物體表面個數n不與T3和T4相關，且其數值遠大于v和m的值。因此，可以忽略T3和T4，則加速方法的時間復雜度T可表示為：

2.4.2 空間復雜度

從加速方法的實現流程可知，加速方法的空間復雜度可以類似于時間復雜度，運算時的內存占用絕大部分是跟蹤的射線與分區中的物體表面進行的求交點計算，且循環多個分區。因此，整個加速方法的空間復雜度S可表示為：

3 仿真與結果分析

3.1 驗證改進算法精確度

首先驗證改進算法相對于原始算法結果的精確性。改進算法可應用在任意三維空間環境模型，且復雜度越高改進算法提升效率比越高。對如圖6 所示的AB路徑進行功率預測值的對比仿真，各射線跟蹤算法運行時原始數據都相同。本次驗證的環境模型為12 m×12 m×4 m 的立方體結構，內部有18 個障礙物模型：9 個高度為2 m 的障礙物模型與9 個高度為2.5 m 的障礙物模型，共有105個需要遍歷的物體表面。射線發射源發射頻率60 MHz，發射源與接收點均為全向天線［15］，發射源間隔角為1.8°可向四周不同方向均勻發射2 萬條射線，物體墻面介電常數為4.5［16］。起點A坐標為（-5.5，3，1），止點B坐標為（3.5，-3，1），總長約為10.8 m。在該路徑上有55 個預測仿真點，間隔為0.2 m。仿真使用Intel i7 處理器，8 GB內存，Matlab仿真平臺。

圖6 用于驗證算法改進前后功率預測值的仿真路徑Fig.6 Simulation path for verifying power prediction values before and after algorithm improvement

圖6所示AB路徑上55個射線發射點，分別使用原始的射線跟蹤算法和改進算法之間做預測對比分析。從圖7 兩種算法的功率折線圖可以看出，原始跟蹤算法與改進跟蹤算法的折線圖幾乎重合。在同一個三維環境模型下，跟蹤路徑AB上各發射點發射出的模擬射線，使用改進算法的預測折線與使用原算法的預測折線差距非常小，僅在仿真位置39 位置處改進后的算法功率比原始算法高了0.03 dBm，偏差很小，誤差在合理范圍內，可視為精確仿真。

由圖7 的折線圖可以驗證，改進后的跟蹤算法幾乎沒有降低原始算法的預測精確度。

圖7 在路徑AB上使用原算法和改進算法的功率預測值對比Fig.7 Comparison of power prediction values using original algorithm and improved algorithm on path AB

3.2 驗證改進算法高效性

在同一個三維環境模型下驗證射線跟蹤加速方法的高效性。選取圖8 中的A、B、C三個位置點作為發射基站位置，依次使用無空間分割的原始算法、靜態多級分割方法以及動靜結合分割的改進方法運行仿真。

圖8 改進算法高效性驗證的基站天線分布Fig.8 Distribution diagram of base station antennas for improved algorithm efficiency verification

在驗證對比中除了各空間分割方法不同以外，其余的所有初始數據、方法都相同。多次運行記錄不同跟蹤算法的總運行時間與總求交次數的數據，結果如表1所示。

表1 高效性對比驗證結果Tab.1 Efficiency comparison and verification results

因為射線跟蹤算法以及本文的加速方法都是電磁仿真中的確定性仿真模型，所以多次重復運行同一個加速方法，總求交次數固定不變，而總運算時間會因為電腦CPU 占用等一些客觀因素的影響發生波動。因此需要對總運算時間做多次仿真，計算均值。圖9 為多次運行原始算法與本文的加速方法對圖8 環境模型中A基站做信道仿真的總運算時間波動折線圖。

根據圖9 所示的時間折線圖計算可知，使用原始算法多次重復仿真后其算法運算時間均值時間為23.172 2 s、標準差為0.056 4 s；使用本文的加速方法的運算時間均值為8.953 8 s、標準差為0.055 5 s。對于不同射線跟蹤加速方法的標準差相差很小，因此客觀因素對仿真時長的影響程度都相同，后續驗證可使用多次仿真的時間均值來表示各射線跟蹤方法運算的總時間。

圖9 仿真時間波動圖Fig.9 Fluctuation graph of simulation time

比較分析得出以下結論：

1）動態與靜態結合的空間分割加速效果最優，靜態空間分割次之，傳統的射線跟蹤算法計算效率最差。

2）無分割的傳統射線跟蹤算法因為需要跟蹤的射線數量確定，且每次對碰撞點的檢索都是檢索全區域，所以在該情況下A、B、C三個基站點的求交次數相同。但由于三個基站所處的位置不同，所以經過的路徑不會相同，算法運行的總時間會有一定差異。

3）靜態空間分割若沒有依照環境密集度情況合理劃分，比如分割等級分的過高，算法的加速效果將變得不明顯，且反而有可能增加運算時間。而恰當的靜態分割加速方法對比傳統算法，在當前環境模型下的計算效率可以有61.4%以上的提高。

4）使用靜態與動態空間分割相結合的加速方法對比只使用靜態空間分割的算法，在當前環境模型下的計算效率在已經提升的基礎上還會有9.8%左右的提高。

當然，在減少算法計算時間的同時也需要保證計算的精度。對本文的三維射線跟蹤加速方法，因為射線與物體的相交點并不會發生改變，所以算法的精確度不會受到影響。

3.3 實驗分析

使用本文中動態結合靜態空間分割的射線跟蹤加速方法與文獻［12］中提出的改進算法作性能對比分析，都使用圖8所示的三維環境模型以及A、B、C三個坐標源點作為發射基站位置，且其余仿真參數都相同。對比結果如表2所示。

表2 本文的射線跟蹤加速方法與文獻［12］算法對比Tab.2 Comparison between the proposed ray tracing acceleration method and algorithm in literature［12］

根據上述結論比較分析后可知，本文使用的加速方法相對于文獻［12］算法在提升射線跟蹤算法運算效率方面有更好的加速指標，進一步減少了使用射線跟蹤算法對無線信道傳播仿真的時間。因此，本文的加速方法是一種更快速的仿真算法。

仿真驗證其他的如圖10 所示簡單三維環境模型，依次使用原始算法、靜態分割方法以及動靜結合分割的方法運行仿真，各完整運行50次后得到的平均運算時間數據如表3所示。根據動靜結合加速方法的空間分割規則，該模型只能做基本的一次分割，不能再繼續深入劃分。

圖10 簡單三維環境模型Fig.10 Simple 3D environment model

表3 最小加速效果驗證結果Tab.3 Verification results of minimum acceleration effect

從表3可得出以下結論：

1）動態與靜態結合的空間分割加速效果最優，靜態空間分割次之，傳統的射線跟蹤算法計算效率最差。

2）使用靜態空間分割的射線跟蹤加速方法對比原始算法計算效率提高了50.2%；而使用靜態與動態空間分割結合的加速方法對比只使用靜態空間分割的方法，計算效率在已經提高的基礎上還能提升8.9%。

3）和3.2 節中表1 的結論數據對比分析可以驗證，三維環境越復雜動靜結合的加速算法提升計算效率的比例也會越高。

上述結果表明本文所提出的改進算法不僅可以大幅度地提高原始算法的計算效率，而且幾乎沒有降低算法的預測精度。驗證了改進算法的高效性，也很好地改善了射線跟蹤算法模型中預測精度與計算效率之間的矛盾。

4 結語

本文使用射線跟蹤算法對空間環境中無線信號的傳播進行信道仿真與三維建模，并且分析了原始射線跟蹤算法在仿真計算過程中將產生大量的無用求交運算，導致計算效率過低的問題。因此，本文根據原始算法的不足之處，提出了一種提高計算效率的射線跟蹤加速方法。該方法結合了動態空間分割與靜態空間分割的方法，減少了跟蹤計算過程中的無用求交點次數，提高了射線跟蹤算法計算效率，降低了算法仿真運算時間。這能為電磁環境的實時仿真計算提供解決方法，是一種實用的射線跟蹤加速方法。

然而，本文設計的射線跟蹤改進加速方法有局限性，僅考慮了從減少求交次數的方面來減少仿真的運行時長。后續可以嘗試，在保證算法精確性的前提下動態減少射線發射端需要跟蹤的射線管數量，以達到加速射線跟蹤仿真計算的目的。