黃巍浩,何思遠,張云華,李婉聰,楊澤望,劉 建
(1.武漢大學電子信息學院,湖北武漢 430072;2.中國人民解放軍63921部隊,北京 100094)
靶標設計和優化對提升武器系統性能有著重要的支撐作用,靶標雷達特性預估和分析則是靶標設計優化的重要基礎。多次散射結構廣泛存在于各類靶標外形結構和靶場應用之中,如靶場角反射器及其陣列的使用、飛機進氣道結構、艦船復雜上層建筑形成的多次散射結構等[1-2]。如何準確、高效地實現對含多次散射復雜靶標雷達特性的快速預估和分析,仍然是電磁散射建模領域的重要研究課題[3-4]。高頻電磁散射建模方法具有計算效率高、物理機理清晰、工程應用性強的優點,非常適合于大型靶標寬帶雷達特性的快速計算和分析[5]。
然而,在靶場尤其是大型靶標設計的實際工程應用中,為了獲取合理的設計結果,需要大量重復計算和優化,現有高頻計算方法依然難以滿足靶標設計高效性的需求,靶標雷達特性快速預估和分析面臨諸多困難和挑戰。這些挑戰和困難主要來自以下幾個方面。首先,對電大尺寸靶標的高分辨一維距離像(HRRP),二維雷達圖像(SAR)等目標寬帶雷達特性的需求加劇了電磁計算的困難[6]。靶標寬帶電磁散射特性仿真必須克服雷達寬帶信號多個頻點和目標多姿態采樣所帶來計算量大的挑戰[7]。其次,為了滿足靶標特性計算精度的需求,含多次散射結構(靶場角反射器及其陣列、腔體等)在電磁建模過程中,需要計入三次及三次以上的電磁高階耦合作用,開展三次以上的高頻射線追蹤。這類多次的射線追蹤過程極其復雜,計算非常耗時,導致靶標雷達特性預估面臨計算效率低下的問題[8]。第三,含多次散射結構靶標涉及復雜多次電磁散射機理,對其寬帶雷達特性形成機制的分析同樣面臨巨大挑戰。多次散射形成的等效視在中心(距離像位置或者散射中心位置)往往偏離目標區域[9],對多次散射機理和等效視在中心位置的分析,目前還缺乏有效的分析手段。
為了提升電大尺寸復雜目標高頻方法計算效率,研究者們采用加速算法和并行計算技術改進高頻仿真過程中最耗時的射線追蹤模塊。Suk 等人提出利用多分辨率網格算法減少射線管的數量從而減少計算量[10];Glassner 等人提出空間劃分算法,利用面元之間的空間關系,通過減少射線追蹤過程中求交判斷的次數來減少計算量[11];隨后KDTree 和八叉樹等算法在此基礎上不斷發展[12];國內浙江大學團隊利用GPU 平臺同時計算眾多射線管,通過并行計算來提高計算效率[13-14]。
除了上述加速算法和GPU 并行技術,CPU 多核技術近年來迅猛發展,為并行電磁計算的發展創造了發展空間和條件[15-16]。本文針對含多次散射復雜靶標雷達特性的快速計算難題,提出了一種CPU 多核并行技術、GPU 硬件加速技術和KDTree 遍歷技術相結合的靶標電磁散射加速計算方法。在此基礎上,建立了帶腔艦船、角反射器陣列與艦船等含多次散射復雜靶標的電磁散射(高頻)模型,可以滿足大型靶標進行快速評估的需求。
針對多次散射特性分析困難,本文進一步提出了多次散射距離像位置快速預測方法。通過多次散射射線分集結合射線路徑相位理論預估,預測了多次散射等效視在中心的位置。將多次散射等效視在中心理論預估位置與高頻仿真雷達特性進行比較,揭示了腔體等含多次散射結構復雜的多路徑散射機理與作用過程。通過具體數值算例,建立了靶標多次散射雷達特性與其多次散射幾何結構之間的映射關系。
本文結構安排如下:首先針對靶標雷達特性仿真效率低下的問題,提出了基于三種計算機加速技術的靶標電磁散射高頻建模加速方法。然后,針對多次散射特性分析困難,基于電磁場射線理論相位預估,發展了一種多次散射距離像位置快速預測方法。最后,通過數值仿真和分析,表明了本文方法的有效性。
圖1 展示了傳統的靶標雷達特性高頻建模流程。

圖1 雷達目標特性獲取流程圖
針對靶標建模對象,首先需要建立對應的高精度CAD 模型,并將目標表面剖分為若干面元,后續根據目標的點面信息進行電磁計算。在對模型預處理之后,將目標散射貢獻按照高頻散射機理進行分解,可將復雜目標電磁散射貢獻分解為邊緣繞射場、面元散射及多次射線場等貢獻之和。在給定計算姿態下,通過射線追蹤確定電磁場的傳播路徑,并結合物理光學法計算該姿態下的總射線場。進一步可結合時頻變換獲取寬帶雷達特性數據,包括高分辨一維距離像數據以及合成孔徑雷達二維圖像數據。上述建模流程中,復雜多次射線追蹤和大量的目標寬姿態寬頻帶采樣是電磁計算耗時的兩個主要原因。
在射線追蹤過程中,需要所有射線和所有面元兩兩進行照亮判斷,此時的時間復雜度為O(N2),N為面元數量。針對電大尺寸復雜靶標高頻區散射,N往往數以萬計,此時要獲取雷達寬帶特性數據的計算代價是巨大的。尤其是當目標存在多次散射結構時,每一次反射點記錄需要重復射線與面元兩兩判斷的過程,成倍地增加了計算時間。考慮到射線追蹤的過程中,各個射線之間的追蹤互不相關,非常適合大規模并行計算。本文通過GPU 平臺來對多次射線追蹤進行加速。GPU 上有大量的計算單元,在每個計算單元上創建線程進行射線的追蹤運算,實現所有射線并行追蹤的功能,此時射線追蹤的時間復雜度降為O(N)。同時,利用目標各面元的空間位置關系,進一步引入KD-Tree 數據結構,僅對射線與發射方向附近的面元進行求交測試,此時射線追蹤的時間復雜度進一步下降到O(log2N)與O(N)之間。
傳統仿真串行地依次計算每一個姿態和頻點下的回波數據,使得目標寬頻帶和寬姿態采樣的計算非常耗時。考慮到各個姿態下的計算流程相同,且計算過程中的中間數據互不相關,本文采用功能更為強大的CPU 來創建指定數量的多線程實現并行計算。由CPU 創建的每一個子線程在射線追蹤環節都將訪問GPU,之后繼續在CPU 端完成后續的射線場計算。而GPU 作為計算過程中的共享資源,在被其中一個子線程訪問時,需要阻止其他線程的訪問。本文利用OpenMP 編程模型來協調線程和管理共享資源,使得雷達目標多姿態并行計算順利進行。如圖2 所示,通過OpenMP 中的critical指導語句設立臨界區,臨界區內的內容同一時刻只能有一個線程訪問,將訪問GPU 的代碼塊放置在臨界區內,可以保證同一時刻只有一個線程在訪問GPU,從而實現了共享資源的保護。

圖2 多姿態并行計算模型
在上述過程中,子線程數量并非越多越好,需要確定合理的子線程數量。由于不同線程之間對共享資源GPU 存在競爭關系,當創建一定數量的線程之后,GPU 將達到滿負載工作狀態,整體效率不會再因為線程數量增多而增加。本文根據線程在GPU 和CPU 上運行的時間長短來創建合理數量的線程。如果把GPU 端進行射線追蹤的時間記為A,把CPU端后續射線場計算部分的時間記為B,那么最合理的線程數N可以由式(1)得到:
圖3 展示了兩種不同情況的用時占比。第一種情況下,GPU 上用時與CPU 上用時相等,此時創建的最佳線程數為2。第二種情況下CPU 上用時情況為GPU 上的2 倍,最佳線程數為3。兩種情況下,GPU 上線程運行狀況如圖4 所示。可以看出,此時GPU 已經達到滿負載工作,無法通過創建更多的線程來提高系統效率。

圖3 GPU用時不同占比的兩種情況

圖4 兩種情況下GPU運行情況
本文提出了多次散射距離像位置快速預測方法。通過多次散射射線分集結合射線路徑相位理論預估,預測了多次散射等效視在中心的位置。
針對復雜目標多次散射,本文通過射線追蹤記錄電磁波與目標相互作用的反射點,獲取電磁波傳播過程中的射線路徑集合,計算每一條射線的光程。根據射線光程可以快速推導多次散射等效視在中心的位置,即距離像位置或者散射中心位置。一維距離像峰值的射線理論預測值X可以按式(2)計算:
式中L為光程,R為目標參考系原點到雷達的距離。為了說明光程的計算過程及其物理含義,本文選取多次散射典型結構,以矩形腔體為例,進行具體的分析。當雷達入射方向與腔體橫截面平行時,射線在一個平面內傳播,此時可以清晰簡潔地分析橫截面上射線的傳播路徑,并給出峰值預測結果。
以腔體口徑面為界,將射線分為腔體內部射線以及腔體外部射線,然后分別計算光程。圖5給出了腔體內部一條射線在口徑面某一點入射,經四次反射后返回雷達的傳播路徑。

圖5 腔體內部射線追蹤情況
按反射情況可以將腔體內部的射線分為5段,分別標記為L1,L2,L3,L4,L5。將L1,L3,L4段以腔內壁為軸翻轉,根據幾何反射原理,翻轉后的L′1,L2,L′3在同一條直線上,L′4,L5在同一條直線上,由此可以得到該射線長度為
由式(3)可以看出,射線光程在腔體內的長度僅由入射角度α和腔體深度D決定,與入射位置無關。所以對于在口徑面其他點入射的射線,光程均相等。而當反射次數增加時,同樣可以通過上述方式計算光程。
圖6 用于說明如何計算腔體內外總的射線長度。為方便分析,如圖建立二維坐標系,并將入射射線與口徑面的交點和出射射線與口徑面的交點記為y1,y2,與雷達的距離記為a,b,則腔體外射線長度可記為a+b。根據幾何關系,可以得到入射方向和反射方向平行,此時腔體外光程可以為y1y2中點到雷達距離的兩倍,所以需要先求得y1y2中點。根據幾何關系,有如下關系式:

圖6 自由空間中射線路徑圖
式中m,n為腔體底部被接觸點分割的兩部分,且二者之和為腔體底部長度H,可以求得
式(5)說明入射點與出射點的中點位置坐標僅受入射角度和腔體尺寸影響,從不同點入射的射線在腔體外的長度也相等。所以,針對矩形腔體,在同一入射姿態α下形成的所有4 次反射的射線路徑光程相同,一維距離像位置相同。
圖7 是預測腔體目標一維距離像的示意圖。A,B分別為射線入射時與口徑面的交點以及反射時與口徑面的交點在雷達徑向方向的投影。O點為目標坐標系原點在雷達徑向方向上的投影,記坐標原點與雷達之間的距離為R,OA,OB分別為O點和A,B點之間的距離,有R=a+OA=b+OB。一維距離像的峰值坐標可以由式(6)求出,其中OA+OB的值可以通過上文求出的y1y2中點和雷達間的距離求出。

圖7 自由空間中腔體目標成像示意圖
三面角和腔體是靶標設計中常見的具有多次散射的典型結構,本節首先通過本文加速方法,獲取三面角和矩形腔體兩種結構的雷達特性。并將矩形腔體距離像射線理論預估位置與高頻仿真距離像特性進行比較,兩者吻合一致。進一步,針對帶腔艦船、角反射器陣列與艦船復合目標開展了仿真與特性分析,仿真結果表明了本文方法的有效性。本文仿真環境如下:CPU 型號為Intel Core i7-9700,內存大小為16G,內含8 核。GPU 型號為NVIDIA GeForce GTX 1660Ti,顯存大小為8 192 MB。
三面角的幾何模型如圖8 所示,面元數為34 292,點數量為17 361。入射射線與z軸夾角記為俯仰角θ,入射射線在XOY平面上的投影與x軸夾角記為方位角φ。三面角雷達參數設置如下:θ角范圍為0°~90°,φ角為45°,頻率為3 GHz,極化方向為垂直極化。

圖8 三面角模型

圖9 三面角RCS結果對比
分別在無加速、利用GPU和KD-Tree數據結構加速射線追蹤過程以及CPU 和GPU 混合編程三種情況下記錄所需用時,加速情況如表1所示。此時線程調用GPU 的時間與在CPU 端運行的時間相近,所以通過OpenMP的協調,與不調用OpenMP相比提高了約兩倍的計算效率。

表1 不同方法計算三面角RCS曲線耗時對比
矩形腔體算例尺寸如圖10所示。

圖10 矩形腔體模型尺寸
腔體目標雷達參數設置如下:θ角為90°,φ角范圍為0°~90°,頻率為10 GHz,極化方向為垂直極化,商用軟件FEKO 所用方法為高頻方法,計算結果相比如圖11所示。

圖11 腔體目標RCS結果對比
選取RCS 曲線中的極大值點,計算當φ= 37°時的一維距離像。不模糊距離為3 m,成像分辨率為0.01 m。圖12 給出了腔體目標在該入射角度下一維距離像的數值結果。仿真結果表明,腔體目標在該入射角度下形成了一個強散射中心,距離像峰值坐標為0.43 m。用本文射線路徑相位理論預估的方法計算得到的結果為0.433 1 m,符合仿真結果,說明通過目標的幾何結構預測目標的一維距離像峰值坐標方法的可靠性。

圖12 腔體目標一維距離像(θ = 90°,φ = 37°)
選取RCS曲線中的其他極值點,分別計算φ角度為51.3°,60.2°,66.1°和73°情況下的一維距離像,不模糊距離為5 m,成像分辨率為0.01 m。結果如圖13所示。

圖13 不同入射角度下的腔體目標一維距離像
不同姿態下的射線長度不同,所以得到的距離像峰值位置也不同。由圖可見,隨著入射角增大,即入射射線與腔體口徑法向的夾角增大,射線經由腔體反射形成的強散射中心逐漸遠離參考原點。同時,口徑面在雷達方向的投影面積減小導致進入腔體的射線減少,一維距離像的峰值強度也逐漸減小。
在實際工程中,射線不一定在一個平面內傳播,多次反射的結構也可能是不規則的,此時同樣可以根據本文射線追蹤方法得到射線光程,進而預測目標的一維距離像。選取俯仰角θ= 82°,方位角φ= 37°的姿態進行射線追蹤和分析,目標姿態如圖14所示。

圖14 腔體目標姿態圖(θ = 82°,φ = 37°)
通過記錄所有射線路徑,可以計算得到所有射線長度均為1.824 0 m。坐標原點到雷達接收點之間的距離為0.471 7 m,根據射線路徑相位理論預估方法可計算目標一維距離像中峰值的理論坐標X= 0.440 3 m。該姿態下仿真得到的一維距離像結果如圖15 所示,圖中峰值坐標為0.44 m,與通過射線追蹤預測的坐標一致,說明了通過射線長度預測目標一維距離像峰值的普適性。

圖15 腔體目標一維距離像(θ = 82°,φ = 37°)
圖16所示為某一驅逐艦船模型,艦船長166.7 m、寬17.2 m、高33 m,船艙部件之間存在多種多次散射結構。在艦船的船頭、中部和船尾存在明顯的腔體部件,腔體部件長8 m、寬5 m、高3 m。后文分別以前腔、中腔和后腔指代。艦船除3個腔體部件以外關于x軸對稱。

圖16 帶腔艦船幾何模型
雷達參數設置如下:θ角為90°,φ角范圍為0°~360°,頻率為10 GHz,極化方向為垂直極化,計算結果如圖17 所示。φ角為0°時探測方向正對船頭,散射貢獻主要來自船體前部的平板結構。在90°和270°處RCS 達到峰值,此時主要散射貢獻來源于船側的平板結構。在俯仰角為90°的情況下,關于x軸對稱目標的全姿態RCS 曲線也具有對稱性,關于方位角為0°和180°對稱。后腔是唯一不關于x軸對稱的部件,在方位角為20°~110°的范圍內后腔散射貢獻主要來自于腔體內部的多次散射作用,在方位角為110°、200°和290°時,后腔散射貢獻主要來自于腔體外壁的一次作用,這些方位角對應的RCS 值略高于對稱姿態下的RCS 值。對含腔目標而言,如果射線追蹤的次數過少,則無法計入部分角度下腔體的散射貢獻,會影響對應角度下RCS的計算精度。
圖18為帶腔艦船目標在俯仰角為65°,方位角為37°姿態下的一維距離像,簡易艦船為不包含3個腔體部件的艦船模型。對帶腔艦船目標的幾個較強的峰值進行分析,形成峰值的等效散射源與目標幾何結構的對應關系如圖19 所示,射線理論預估和仿真結果對比如表2 所示。其中1,3,4 號峰值點分別由3個腔體部件形成,前腔和中腔開口方向與船頭一致,然而前腔截面偏離XOY面角度約3°,導致其距離像幅度較低。2 號峰值點為船體中部船艙各平板之間的三次作用形成。

圖18 帶腔艦船與簡易艦船一維距離像對比圖

圖19 一維距離像峰值與目標幾何結構對應關系示意圖

表2 艦船目標一維距離像位置驗證
腔體結構不僅在某些角度下能形成強散射中心,腔體內部的多次散射還能影響其等效視在中心與部件幾何中心投影之間的位置關系。如圖20所示,前腔與中腔的開口方向一致,在該入射姿態下具有相同的光程,視在中心與幾何中心投影均相距5.2 m。后腔在該入射姿態下的散射貢獻來源于腔體內壁與底部的二次作用,視在中心與幾何中心投影距離為3.4 m。可以看出,多次散射使得目標的等效散射源偏離了原有目標空間范圍,偏離程度隨著射線光程的增加而增大。

圖20 腔體等效視在中心與腔體幾何結構關系示意圖
本文以圖21所示的角反陣列與艦船復合模型為例,快速評估復合模型的散射特性及角反陣列回波干擾信號的有效性。艦船長119 m、寬10.8 m、高12.88 m,角反射器陣列分別包含邊長為0.5 m、1 m、2 m 的角反射器各4 個,這些角反射器開口朝向各不相同,復合模型的網格數量為278 000。

圖21 角反射器陣列與艦船目標復合模型
當雷達參數為俯仰角60°,中心方位角45°,中心頻率10 GHz,極化方式為垂直極化,分辨率為1 m 時,獲取圖23 所示的復合模型SAR 圖像數據需對22 500個不同參數(150個頻點、150個方位角點)的目標RCS 進行采樣計算,不同方法的耗時如表3所示。本文提出的CPU多核并行技術、GPU硬件加速技術和KD-Tree 遍歷技術相結合的靶標高頻電磁散射加速方法與傳統高頻仿真方法相比,將普通性能計算機獲取單幅大型艦船與角反射器陣列復合散射SAR圖像的耗時由186 h降至0.35 h,可實現大型艦船靶標散射特性的單機快速評估。此時線程調用GPU 的時間遠小于在CPU 端運行的時間,所以通過OpenMP 的協調,與不調用OpenMP相比提高了約6倍的計算效率。

表3 不同方法計算單幅大型艦船與角反射器陣列復合散射SAR圖像的耗時對比
圖22為艦船與角反陣列復合散射的一維距離像特性,雷達參數設置如下:俯仰角為60°,方位角為45°,極化方向為垂直極化,不模糊距離為150 m,分辨率為1 m。圖中藍色點狀區域為目標在雷達視線方向的投影示意,依據距離像峰值與目標幾何結構的對應關系可知,6 個峰值點全部由開口朝向雷達的9個角反射結構三次散射形成。其中,1,4,5號峰值點分別由3對位于同一雷達分辨單元的角反射結構形成,而2,3,6號峰值點則分別由單個角反射器形成。圖23為艦船與角反陣列復合模型的SAR圖像,圖中的9個強散射中心同樣全部由這9個角反射結構的三次散射機理形成。

圖22 三面角陣列與艦船復合模型距離像

圖23 三面角陣列與艦船復合模型二維圖像
通過上述分析可以發現,在雷達觀測的一維距離像和SAR 圖像特性中,邊長0.5 m 的角反射器產生的強后向散射回波足以掩蓋艦船目標自身的回波信號,且由于雷達對角反射器三次散射機理的等效視在位置位于三面角頂點,因此,在艦船周圍放置合適尺寸及方向的角反陣列,可以完全改變艦船的雷達特性。
本文針對靶標特性建模的需求,提出了結合OpenMP 和GPU 混合編程的加速算法,提高了獲取雷達目標回波數據的效率。對于存在多次作用結構體的目標,本文計算了雷達目標包含高階散射場的總體回波信息,通過與FEKO 對比驗證了計算結果的準確性。通過對不同方法計算的用時說明了本文方法的高效性,并利用雷達目標特性數據獲得目標的一維距離像和SAR 圖像,根據射線理論預測一維距離像峰值坐標,解釋了距離像峰值的物理意義,雷達圖像所展示出的特征與目標幾何形狀有較好的關聯性,說明了計算結果的可靠性。