蘇 秦，顧宗靜，林中朝，趙勛旺，張 玉

(西安電子科技大學(xué) 天線與微波技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710071)

雷達(dá)散射截面(Radar Cross Section，RCS)的計(jì)算是伴隨著雷達(dá)目標(biāo)的檢測和識別技術(shù)發(fā)展起來的，它不僅可以應(yīng)用于現(xiàn)有的各種目標(biāo)，又可以用來預(yù)估和優(yōu)化未來的信息化系統(tǒng)，在工程應(yīng)用領(lǐng)域具有重要作用．尤其是近幾十年來，隱身技術(shù)與精確目標(biāo)識別技術(shù)日新月異，這對電大尺寸復(fù)雜目標(biāo)的RCS計(jì)算提出了更高的要求．對于電尺寸很大的散射體，通常采用物理光學(xué)法(Physical Optics method，PO)等高頻算法求解，然而，高頻算法難以保證足夠的精度．傳統(tǒng)的高精度電磁場數(shù)值算法如矩量法(Method of Moment，MoM)和有限元等算法往往只能解決幾十個(gè)波長的電磁問題．多層快速多極子算法(MultiLevel Fast Multipole Algorithm，MLFMA)的求解規(guī)?？梢赃_(dá)到上百波長[1]，但對于上千波長規(guī)模的問題在資源有限的情況下仍然很難求解．這迫使人們必須找到一種高效可實(shí)現(xiàn)的數(shù)值仿真方法．

區(qū)域分解方法(Domain Decomposition Methods，DDM)是使用“分而治之”的思想解決問題的方法．它先對整體目標(biāo)進(jìn)行分區(qū)，依據(jù)其物理或電磁散射特性劃分為多個(gè)子區(qū)域獨(dú)立求解，然后通過考慮子區(qū)域間的耦合關(guān)系，最終得到原問題的解．由于每個(gè)子區(qū)域處于獨(dú)立求解空間，并且具有相對較少的未知量，這就使得大型電磁仿真分析的快速實(shí)現(xiàn)成為可能．在有限元與有限差分領(lǐng)域，DDM已經(jīng)成功實(shí)現(xiàn)[2-3]．而在積分方程領(lǐng)域，如MoM，近年來對于DDM的研究也越來越受到學(xué)者們的關(guān)注．一種廣泛采用的重疊型DDM中[4]，相鄰兩個(gè)子區(qū)域間存在重疊部分，這種方法可以使得電磁流更加逼近真實(shí)電磁流．但是重疊區(qū)的引入給建模帶來了不便，同時(shí)區(qū)域延伸也導(dǎo)致了計(jì)算量增加．而另一種非重疊型DDM則是在子區(qū)域之間引入人工虛擬面，將原問題分解為很多閉合的子區(qū)，每次迭代通過修正虛擬面上電磁流來保證相鄰子區(qū)域場的連續(xù)性[5]．這種方法節(jié)省了幾何建模的工作量，具有較高的靈活性．然而對于用區(qū)域分解方法解決千波長規(guī)模電磁問題，鮮有有效工作內(nèi)容發(fā)表．這主要是因?yàn)镈DM的并行方案復(fù)雜，并且在計(jì)算電大尺寸問題時(shí)，為了保證足夠的精度，通常需要用MLFMA來降低內(nèi)存消耗與計(jì)算復(fù)雜度，而MLFMA的八叉樹結(jié)構(gòu)，使得其在并行化過程中的負(fù)載均衡難以實(shí)現(xiàn)．

針對這一問題，文中給出一種高效的基于積分方程的并行非重疊區(qū)域分解方法(Integral Equation-based Nonoverlapping Domain Decomposition Method，IE-NDDM)．與文獻(xiàn)[6]中所有子區(qū)域同時(shí)計(jì)算的并行方案相比，并行IE-NDDM將所有計(jì)算機(jī)資源集中對單個(gè)子區(qū)域求解，避免了子區(qū)域大小劃分不均等所導(dǎo)致的負(fù)載不均衡，同時(shí)待求解問題的規(guī)模也大大提升．每個(gè)子區(qū)域內(nèi)部以及子區(qū)域間耦合量的求解均用并行MLFMA進(jìn)行加速．改進(jìn)的平面波自適應(yīng)劃分策略使得并行MLFMA的數(shù)據(jù)劃分更加均衡，進(jìn)一步提高了并行IE-NDDM的計(jì)算能力和效率．子區(qū)域間的耦合量計(jì)算采用場迭代的方式，降低了計(jì)算過程中的內(nèi)存消耗．文中最后給出一個(gè)實(shí)際應(yīng)用中的千波長艦船的散射特性分析，數(shù)值結(jié)果表明，該方法能夠有效應(yīng)用于電大尺寸散射體的快速精確仿真．

1 IE-NDDM原理

與重疊型DDM相比，IE-NDDM僅在相鄰子區(qū)域切割處增加虛擬面，使各個(gè)子區(qū)域封閉．此處以金屬目標(biāo)劃分為兩個(gè)非重疊型子區(qū)域?yàn)槔M(jìn)行分析，對于多個(gè)子區(qū)域情況可以依此類推．

圖1 模型區(qū)域分解示意圖

以子區(qū)域Ω1為例，在平面波{Einc，Hinc}入射時(shí)，根據(jù)等效原理，其表面電場積分方程(Electric Field Integral Equation，EFIE)及磁場積分方程(Magnetic Field Integral Equation，MFIE)可以寫為

(1)

其中，{Es，Hs}代表金屬體表面等效源所產(chǎn)生的散射場．對于金屬體表面而言，只有電流J，因此，散射場用等效源的積分形式可以表示為

(2)

(3)

2 IE-NDDM的并行策略

對式(1)中的EFIE及MFIE線性組合起來得到混合場積分方程(Combined Field Integral Equation，CFIE)，并用伽遼金匹配來構(gòu)建出全局矩陣方程．對于n個(gè)子區(qū)域，方程ZI=V可以寫成如下形式:

(4)

其中，當(dāng)i=j時(shí)，Zij為第i個(gè)子區(qū)域Ωi的自阻抗矩陣; 當(dāng)i≠j時(shí)，Zij為子區(qū)域Ωj和子區(qū)域Ωi的互阻抗矩陣;Ii為子區(qū)域Ωi的電磁流系數(shù)矩陣，Vi為子區(qū)域Ωi的電壓矩陣．

(5)

(6)

圖2 并行IE-NDDM外迭代流程圖

MLFMA有兩種基本并行策略，分別是組劃分與平面波劃分策略．如果只使用組劃分策略，那么在高層會存在瓶頸，因?yàn)楦邔拥姆强战M較少，無法分配給較多的并行進(jìn)程．同理，如果只使用平面波劃分策略，那么在低層存在瓶頸．為了提高并行MLFMA的可擴(kuò)展性，需要混合使用這兩種并行策略．比較簡單的混合策略是在低層采用組劃分，在高層采用平面波劃分，通過設(shè)置中間過渡層實(shí)現(xiàn)兩種策略的結(jié)合．然而在過渡層，非空組數(shù)和每個(gè)非空組對應(yīng)的平面波數(shù)均為O(N0.5)．當(dāng)進(jìn)程數(shù)超過O(N0.5)時(shí)，該混合策略在過渡層出現(xiàn)瓶頸．為了解決過渡層的瓶頸，一種逐層漸變劃分策略(HiP)將組劃分與平面波劃分策略緊密結(jié)合，具有較高并行效率[8]．在此基礎(chǔ)上，文中對MLFMA采用一種改進(jìn)的平面波自適應(yīng)劃分(Adaptive Partition，AdP)并行策略[9]，并行編程模型為基于分布式內(nèi)存的消息傳遞接口(Message Passing Interface，MPI)．該策略同樣采用逐層漸變的方式對平面波方向進(jìn)行劃分．不同的是，在某一層對于不同的非空組，HiP對于平面波劃分的份數(shù)相同，而AdP根據(jù)非空組所在的進(jìn)程數(shù)進(jìn)行自適應(yīng)劃分，非空組所在的進(jìn)程數(shù)不同則平面波劃分份數(shù)也隨之不同．設(shè)并行進(jìn)程數(shù)為Np，將基函數(shù)(或權(quán)函數(shù))等分為Np個(gè)部分，每部分基函數(shù)被分配給一個(gè)進(jìn)程．然后，每個(gè)進(jìn)程根據(jù)所分配到的基函數(shù)建立相應(yīng)的子樹．如圖3(a)所示，第2層至第4層的八叉樹被分配到6個(gè)進(jìn)程中，其中進(jìn)程編號記為P0～P5，非空組中標(biāo)注的進(jìn)程編號就是該非空組所在的進(jìn)程．對于圖3(a)中的第3層，該層中的3個(gè)非空組分別存在于1個(gè)、2個(gè)、3個(gè)進(jìn)程中．每個(gè)非空組對應(yīng)的外向平面波和內(nèi)向平面波被均等地分配到其所在的進(jìn)程中．圖3(b)所示為圖3(a)中第2層非空組的平面波劃分示意圖，每個(gè)圓點(diǎn)代表一個(gè)平面波方向．

圖3 八叉樹的數(shù)據(jù)劃分

這種準(zhǔn)二維劃分的AdP對使用的進(jìn)程個(gè)數(shù)無特殊要求，組劃分與平面波劃分無縫結(jié)合．每個(gè)進(jìn)程都參與計(jì)算和通信，這使得各個(gè)進(jìn)程中的通信量和計(jì)算量更加均衡．在聚合-轉(zhuǎn)移-配置(MVP)過程中，向上聚合與向下配置時(shí)，進(jìn)程中外向平面波和內(nèi)向平面波的交換和疊加都在緩沖區(qū)中進(jìn)行，避免了插值過程中的通信．同時(shí)，利用MPI中的非阻塞通信技術(shù)，將進(jìn)程內(nèi)計(jì)算與進(jìn)程間通信相互重疊，進(jìn)一步提升了MLFMA的并行效率，這使得并行IE-NDDM在普通工作站上即可快速求解大規(guī)模電磁問題．關(guān)于AdP的大規(guī)模并行性能詳見文獻(xiàn)[9]．

3 數(shù)值算例

3.1 波音737飛機(jī)的散射特性

此處以波音737飛機(jī)的散射特性為例來證明并行IE-NDDM算法的正確性與高效性. 入射波面向機(jī)頭入射，極化方向?yàn)?+z．飛機(jī)模型劃分為5個(gè)子區(qū)域，如圖4(a)所示，其尺寸為 30.6 m× 29.0 m× 11.7 m．當(dāng)入射波頻率為 1 GHz 時(shí)，表1列出了兩種計(jì)算方式消耗的計(jì)算資源和計(jì)算時(shí)間．可以看出，并行IE-NDDM與FEKO整體解相比，內(nèi)存使用峰值降低86.7%，同樣使用16個(gè)CPU核數(shù)時(shí)，總計(jì)算時(shí)間也縮短了36.1%．圖4(c)～圖4(d)中實(shí)線為并行IE-NDDM數(shù)值計(jì)算結(jié)果，虛線為商業(yè)軟件FEKO用MLFMA整體求解的仿真結(jié)果，圖4(b)為并行IE-NDDM算法外迭代收斂曲線．可以看出，文中提出的并行IE-NDDM算法4步迭代收斂精度已達(dá)到 0.002 63，并且與整體解的結(jié)果曲線吻合良好，尤其在前后向幾乎完全重合．從而證明了并行IE-NDDM算法的準(zhǔn)確性，為計(jì)算更大規(guī)模的散射體提供了理論支持．圖4(c)和圖4(d)中同時(shí)給出了并行IE-NDDM計(jì)算得到的入射波頻率為 10 GHz 時(shí)飛機(jī)的RCS，此時(shí)飛機(jī)的電尺寸達(dá)到 1 000 波長．

表1 頻率為1 GHz時(shí)的計(jì)算資源列表

圖4 飛機(jī)模型及其RCS

3.2 艦船的散射特性分析

此處以艦船模型為例，用并行IE-NDDM算法對其散射特性進(jìn)行計(jì)算．該艦船模型如圖5(a)所示，其尺寸為 153.0 m× 16.5 m× 28.0 m，共劃分為7個(gè)子區(qū)域．入射方向?yàn)槊嫦虼^入射，極化方向?yàn)?+z．外迭代收斂殘差設(shè)為0.01．當(dāng)入射波的頻率為 300 MHz 時(shí)，艦船的RCS計(jì)算結(jié)果如圖5(b)～圖5(c)所示．可以看出，對于該模型并行IE-NDDM與整體解的結(jié)果曲線依然吻合良好．當(dāng)入射波的頻率為 2 GHz 時(shí)，相應(yīng)的電尺寸為 1 020λ× 110λ× 186.7λ．對整體模型剖分所產(chǎn)生的未知量個(gè)數(shù)為 66 171 474，此時(shí)用現(xiàn)有計(jì)算機(jī)資源已很難得到整體解．采用并行IE-NDDM算法求解，最大的子區(qū)域未知量個(gè)數(shù)為 14 813 784，用32個(gè)CPU核進(jìn)行計(jì)算，總耗時(shí)為 19.2 h．同時(shí)，占用的總內(nèi)存峰值僅為 105.6 GB．這表明在計(jì)算資源有限時(shí)，該并行IE-NDDM算法在計(jì)算時(shí)間和資源消耗上都擁有極大優(yōu)勢．此時(shí)艦船的RCS計(jì)算結(jié)果如圖6(a)～圖6(b)所示．

圖5 艦船模型以及在300 MHz時(shí)的RCS

圖6 艦船在2 GHz時(shí)的RCS

4 結(jié) 論

在非重疊型區(qū)域分解法的基礎(chǔ)上，采用并行多層快速多極子對子區(qū)域內(nèi)部以及子區(qū)域間耦合的計(jì)算進(jìn)行加速．改進(jìn)的自適應(yīng)劃分并行策略提高了多層快速多極子的并行效率，子區(qū)域間耦合計(jì)算用場迭代的方式降低了內(nèi)存．這使得電大尺寸電磁問題的精確高效求解成為可能，在資源有限的情況下，千波長目標(biāo)散射特性計(jì)算的工程難題得到了解決．下一步嘗試將天線加載到電大尺寸平臺上，解決更加復(fù)雜的電磁仿真問題．