一種基于不連續伽遼金方法求解多區域目標散射問題的優化預處理器

張慧雯 黃曉偉 吳比翼 盛新慶

（北京理工大學信息與電子學院應用電磁研究所，北京 100081）

引 言

在電磁計算應用中，由介質和金屬結合的多區域結構有著廣泛的應用價值.典型的介質金屬結合的多區域目標，例如天線-天線罩系統、微波天線陣列等，無法通過建立單一等效模型得到準確的電磁計算結果. 而一些曾經使用單一區域模型來計算電磁特性的目標，也可以通過細化區域結構的方式得到更精確的結果，比如飛行器、導彈模型等. 因此，研究多區域目標的高效、高精度電磁計算方法是十分必要的.

為了提高多區域目標的電磁計算效率，我們使用面積分(surface integral equation, SIE)[1]矩量法(method of moment, MoM)[2]以更高效地求解分塊均勻目標的散射問題，因為其相較于體積分方程具有更少的未知數. SIE 方程一般使用Rao-Wilton-Glisson(RWG)基函數進行離散[3]，但RWG 基函數的連續性使得其很難處理多區域目標的區域結合邊界[4]. 通常的解決方法是在邊界處增加區域連接模型(connectregion modeling, CRM)以保證邊界的連續性條件，但這需要共型區域剖分[5]. 然而當目標具有精細幾何結構時，共型剖分是十分繁瑣且耗費時間的[6-7].

為解決上述問題，區域分解技術(domain decomposition method, DDM)被前人提出，以在不同區域中使用非共型且相互獨立的幾何剖分[8]SIE 方程的區域分解可以分為兩大類：封閉式區域分解(體分解) (volume-based SIE-DDM)和非封閉式區域分解(面分解)(surface-based SIE-DDM)[9]. 封閉式區域分解以體積為基本區域單元，在子區域的連接處需要施加Robin傳輸條件(TCs)以保證子區域電流的連續性[10-11]，而這將帶來更多的未知數. 非封閉式區域分解方法中，不連續伽遼金 (discontinuous Galerkin, DG)方法的研究最為深入[12-13]. 在不同面區域的邊界處，DG 方法使用半RWG 基函數進行剖分，同時將電流的連續性條件直接加在半RWG 的阻抗矩陣元素上. DG 方法無需增加未知數，因此該方法的非共型剖分可以更容易地運用在復雜目標中. 近年來，學者已將DG 方法成功運用于完美電導體[14-15]、均勻介質[16]和阻抗邊界表面[17-18]等目標的計算中，而對于多區域結構的目標，該方法還未有深入的研究.

對于電大目標的仿真，預處理技術是一種加速迭代收斂的有效方法[19-20]. 近場預處理器(near-field preconditioner, NFP)是多層快速多極子方法(multilevel fast multipole algorithm, MLFMA)[21]常用的預處理方式. 對于復雜的多區域目標，其近場矩陣元素非常多，增加了預處理矩陣求逆的運算時間. 距離稀疏預處理器(distance sparse preconditioner,DSP)使用距離限制減少了預處理器的矩陣元素，可以明顯減少計算時間和內存需求[22]. 進一步，我們提出了優化距離稀疏預處理器(optimized DSP, ODSP)，不同積分算子的矩陣塊，其數值大小隨距離變化的趨勢不盡相同. 算子的預處理矩陣塊可以保留更少的矩陣元素而并不會大程度削弱其預處理的效果. 對于多區域目標，O-DSP 對不同的子區域矩陣塊使用不同的構建方法，使之合成為完整的預處理矩陣. 與DSP 相比，O-DSP 保持了相近的加速收斂效果，而且構建的預處理矩陣元素數目大幅降低.

1 理論與公式

1.1 多區域目標的電磁散射計算

考慮圖1 所示的交界面模型，這里 Ωi表示不同的體積區域，可以是均勻介質或金屬. 金屬 δi=0，均勻介質 δi=1. 任意兩個空間區域的交界面用S ij表示，默認法相方向n?ij以區分面的外表面和內表面. 若交界面的兩個空間區域至少一個是金屬區域，則只需要在表面上定義電流，用 Δij=0表示；當交界面的兩個空間區域均為介質時，還需要定義磁流，用Δij=1表示.

圖1 任意交界表面與其兩側體積區域示意圖Fig. 1 An arbitrary surface with its internal and external regions

1.2 DG 方法

在式(1)、式(2)所表示的原方程中所有邊界上建立方程(10)，并乘以 ?jαijαii0/ki加到原方程中，用同樣的基函數進行離散，以消去式(7)L算子元素中的第五項基函數是半RWG 的計算.

對于復雜目標，使用MLFMA 來加速計算[23]. 對于DG 建立的離散方程同樣適用. 唯一的區別在于邊界的半RWG 需要進行一個三角形的積分計算.

1.3 預處理方法

多區域目標SIE 方程迭代計算的收斂速度較慢，因此需要進行合適的預處理. 預處理矩陣的定義如下：

圖2 多區域交界處示意圖Fig. 2 The junction of multiple regions

圖3 介質球體不同算子矩陣塊內元素模值隨距離約束的變化Fig. 3 The variation of modulus value with distance constraint of a dielectric sphere

數值實驗發現，在K 算子矩陣塊中使用式(13)的預處理方式和原DSP 的預處理方式在迭代速率上相近，但是預處理矩陣元素進一步減少. 結合之前L算子矩陣塊的預處理方法，簡稱這種方法為 ODSP.

矩陣塊中的未知數為N，則近場矩陣元素個數為O(N2). 式(12) 中距離判斷與介質波長有關，因此其矩陣元素數量為O(N)，由于有很多元素仍被約束在距離判斷內，O(N)的實際系數很大. 式(13) O-DSP 中矩陣每一行只會保留4 個或2 個元素，因此整體元素數將小于4N，最終整體的預處理矩陣元素數將會明顯減少，且介質表面未知數越多，減少效果越明顯.

2 數值結果

利用廣義最小殘差 (generalized minimum residual,GMRES)求解器來進行仿真，Krylov 子空間的大小為100，迭代殘差為1 0?3. 預處理求逆使用大規模并行求解器(multifrontal massively parallel solver, MUMPS)，計算平臺為2.2 GHz CPUs 和 1 TB 內存的 Intel Xeon E7-4850.

如圖4 所示，半徑為0.5 m 的介質球，相對介電常數為5，被分為4 個子區域，剖分尺寸為0.005 m，為介質波長的11/100，頻率為300 MHz 的平面波沿-Z方向入射. 利用DSP 和O-DSP 對 VV 極化方向的雙棧雷達散射截面(radar cross section, RCS)(記為VV-RCS)進行對比計算，計算結果如圖5 所示. 可以看出，DSP 和O-DSP 數值計算結果均很好地吻合了理論計算的Mie 求解，不會對求解結果的精度產生影響.

圖4 介質球的區域分解形式Fig. 4 Domain partitioning for the dielectric sphere

圖5 介質球的雙棧VV-RCSFig. 5 The bistatic VV-RCS of the dielectric sphere

L算子的距離約束不變，為介質波長的1 /2，令K算子矩陣塊內元素受不同的距離約束，對比其預處理內存和迭代次數，結果如圖6 所示. 可以看到，當距離閾值在介質波長的1/10～1/2 時，迭代次數幾乎不變，但元素個數隨距離的減小大大減少，說明去掉含奇異點殘留項的元素不會對迭代產生很大的影響.當距離閾值為介質波長的3/50 時，處于是否包含奇異點殘留元素的臨界，迭代次數發生了較為明顯的改變，因此有保留這些元素的必要性.

圖6 K 算子矩陣迭代次數和元素數量隨距離約束的變化Fig. 6 The iteration numbers and element numbers of K operator matrix block, as a function of distance constraint

表1 對幾種不同預處理器的數值效果進行了對比. 對比的稀疏策略(sparse strategies, ST)有如下四種：

表1 不同預處理器的DG 方法計算數值效果對比Tab. 1 Comparison of the n umerical performance of the DG solution with different preconditioner

由表1 可知，和NP 相比，所有的預處理器都增加了收斂速度，其中NFP，DSP，O-DSP 效果最好. 這三種收斂速度較好的預處理器中，O-DSP 矩陣元素數量最少. 和DS-2BDP 相比O-DSP 收斂速度極快，說明K算子矩陣塊中有奇異點殘留項的元素對于矩陣的收斂十分重要，即使數量并不多. 同樣使用ST3 作用在L算子矩陣塊，迭代速率明顯下降，說明L算子矩陣塊需要保留較多的矩陣元素.

圖7 為圖4 中的介質球子區域數量分別為4、8、16 和32 的分解情況. 不同分解區域數對DSP 和O-DSP 的迭代次數和元素比率的影響如表2 所示.區域數越多意味著邊界半RWG 越多，K算子矩陣塊會保留更少的元素. 從表2 中可以看出，兩個預處理器的迭代速率幾乎是一致的，說明子區域的數量并不會影響O-DSP 的效率.

圖7 介質球的不同區域分解Fig. 7 Domain partition schemes for a dielectric sphere

表2 不同分解區域數對DSP 和O-DSP 的迭代次數和元素比率的影響Tab. 2 Comparison of the iteration numbers and element ratios of the DG solution with DSP and O-DSP for different decomposed spheres

如圖8 所示，區域數量為4 的介質球內嵌不同個數半徑為0.15 m、剖分尺寸均為0.005 m 的金屬球.頻率為300 MHz 的平面波沿-Z方向入射，介質球內嵌不同個數的導體塊模型的雙棧VV-RCS 結果如圖9所示. 從圖9 可以看出，使用DG 和O-DSP 的預處理器計算的數值結果很好地吻合了商業軟件FEKO 使用共型剖分計算的結果. 導體塊數量分別為2，4 和8 時三種預處理器迭代次數和元素個數對比如表3所示. 對比DSP，O-DSP 的迭代數有輕微的增加，而矩陣元素數量明顯減少. 由此說明，O-DSP 仍然適用于多區域的目標.

圖8 介質球內嵌不同個數的導體塊模型Fig. 8 Different numbers of conducting spheres embedded in one dielectric sphere

圖9 介質球內嵌不同個數導體塊模型的雙棧VV-RCSFig. 9 The bistatic VV-RCS for the model of multiple conducting spheres embedded in one dielectric sphere

表3 導體塊數量不同情況下三種預處理器迭代次數和元素個數對比Tab. 3 Comparison of the iteration numbers and element numbers by 3 different preconditioners with different numbers of conducting spheres

圖10 所示為一個簡化的射頻芯片模型，包含了所有射頻芯片內部的基本結構. 介質長方體模型尺寸為20 mm×20 mm×5 mm ，介電常數為 2.5，內嵌間距為 1 mm 的 5 片金屬片，其中第一層和第五層用 4個金屬圓柱相接作為整體的支撐，第一層和第二層之間由過孔相連，第三層和第四層之間是介電常數為4.2 的介質體. 剖分尺寸為細金屬柱0.3 mm，其余部分0.5 mm. 一共分為 17 個 DG 區域，總未知數為 85 970 個.

頻率為 40 GHz 的平面波沿-Z方向入射，使用CTF 方程計算，使用O-DSP 作為預處理器， VV極化方向的雙棧 RCS 如圖11 所示，驗證了方程的計算正確性. 迭代收斂速度如圖12 所示， DSP 和O-DSP 的迭代次數分別是 149 和153 ，略有差異，兩者明顯快于NP 的收斂速度. DSP 的矩陣元素數量是 78 061 756，O-DSP 為 46 922 730，數量明顯減少.

圖11 射頻芯片模型的雙棧VV-RCSFig. 11 The bistatic VV-RCS of the radio frequency chip model

圖12 使用不同預處理器求解射頻芯片模型的迭代收斂情況Fig. 12 Convergence histories for the radio frequency chip model under different preconditioners

圖13 所示為四旋翼無人機模型的區域分解示意圖，四個旋翼和支架是金屬的，機身是相對介電常數為4 的介質體. 飛行器長寬約為2 000 mm 和830 mm.12 個子區域的總未知數為415 223 . 介質處的剖分為6 mm，一般導體的剖分為8 mm，旋翼處由于太細剖分尺寸為6 mm.

圖13 四旋翼無人機模型的區域分解示意圖Fig. 13 Illustration of subdomain partition scheme of a fourrotor aircraft model

頻率為 4 GHz 的平面波沿-Z方向入射， VV極化方向的雙棧 RCS 如圖14 所示. 可以看出，O-DSP的 DG 方法和FEKO 符合良好.

圖14 四旋翼無人機模型的雙棧VV-RCSFig. 14 The bistatic VV-RCS of the four-rotor aircraft model

圖15 為利用NP，DSP 和O-DSP 求解四旋翼無人機模型的迭代收斂情況. 可以看出，DSP 和O-DSP均快速收斂，二者速度相近. 表4 為 DSP 和O-DSP預處理器計算四旋翼無人機迭代次數和矩陣元素數量對比，其中O-DSP 矩陣元素數量明顯減少.

圖15 使用不同預處理器求解四旋翼無人機模型的迭代收斂情況Fig. 15 Convergence histories for the four-rotor aircraft model under different preconditioners

表4 DSP 和O-DSP 預處理器計算四旋翼無人機迭代次數和矩陣元素數量對比Tab. 4 Comparison of the iteration numbers and element numbers of the DG solution with DSP and O-DSP for the fourrotor aircraft model

3 結 論

本文使用DG 方法來求解多區域目標問題，并使用O-DSP 來提高仿真效率. 優化的預處理矩陣根據積分微分算子的不同數值特性，分塊進行更恰當的預處理構建. 數值算例證明了預處理器的靈活性和準確性. O-DSP 相比DSP 擁有相近的迭代次數和更稀疏的結構特性，由于K算子的預處理矩陣塊使用了一種更精準的判斷方式保留了恰當的元素，因此ODSP 是一種更加有效的計算多區域目標散射的預處理方式.