翟 暢，林中朝，趙勛旺，張 玉

(西安電子科技大學(xué) 陜西省超大規(guī)模電磁計算重點實驗室，陜西 西安 710071)

隨著雷達(dá)技術(shù)的不斷發(fā)展，人們對于雷達(dá)散射截面(Radar Cross Section，RCS)的分析要求越來越高，而半空間電大尺寸電磁問題的仿真分析就是其中極具挑戰(zhàn)的一類問題。因受到環(huán)境、技術(shù)、資金等問題的限制，目前一些科研機(jī)構(gòu)較難通過實驗測量的方式獲得電大尺寸目標(biāo)的雷達(dá)散射截面數(shù)據(jù)，因此大多數(shù)研究都選擇仿真計算的方式。在過去幾十年間，為了分析此類半空間電大尺寸問題，各種電磁方法層出不窮。在上世紀(jì)90年代，大多數(shù)研究均采用高頻電磁算法進(jìn)行仿真分析[1-4]，如物理光學(xué)法(Physical Optics，PO)或彈跳射線法(Shooting and Bouncing Ray，SBR)等。此類方法簡單易行，計算快速，但是無法保證計算精度。因此，人們開始研究如矩量法(Method of Moments，MoM)、有限元法(Finite Element Method，F(xiàn)EM)等高精度方法[5-9]，但此類高精度方法僅適用于電小尺寸問題。當(dāng)目標(biāo)電尺寸較大時，計算內(nèi)存會非線性增長，即使利用多層快速多極子(Multi Level Fast Multipole Algorithm，MLFMA)這類快速算法，也很難在有限的計算資源下解決如海上艦船這類千波長以上的電大尺寸問題。因此，需要一種合理的方法來在有限資源下解決此類電大尺寸問題。

筆者采用一種“分而治之”的基本思想作為指導(dǎo)，將復(fù)雜問題細(xì)化分割成若干簡單問題，再對其一一解決，這就形成了區(qū)域分解方法[10-12]。近幾年，人們對利用區(qū)域分解方法來解決電大尺寸電磁問題做了大量研究工作。新墨西哥大學(xué)的PENG等[13]提出了一種基于積分方法的內(nèi)罰型區(qū)域分解算法(IP-IE-DDM)，通過在交界面處強加罰函數(shù)的方式，保證區(qū)域電流連續(xù)，無須引入非重疊型區(qū)域分解算法(IE-NDDM)中的人工交界面，但其中沒有討論到半空間問題，并且強加罰函數(shù)可能無法有效地消除交界面反射波。通過對基于非重疊型區(qū)域分解方法的研究發(fā)現(xiàn)[14]，其需要較為復(fù)雜的模型前處理，并且人工虛擬交接面的建立會引入額外的未知量，導(dǎo)致內(nèi)存消耗增加。通過總結(jié)探索，根據(jù)多層快速多極子方法所特有的八叉樹結(jié)構(gòu)，令其自適應(yīng)地劃分區(qū)域，并在交界面處采取嚴(yán)格阻抗求解的方式來保證區(qū)域間電流連續(xù)，可達(dá)到與整體解基本一致的計算精度。將這種方法稱為基于八叉樹的區(qū)域分解算法(OT-DDM)。該方法既有傳統(tǒng)區(qū)域分解算法的優(yōu)勢，避免了復(fù)雜的區(qū)域劃分過程，又具有較高的精度，因而可以滿足半空間電大尺寸雷達(dá)散射截面計算的需求。

半空間問題相對自由空間問題，關(guān)鍵技術(shù)難點是如何處理分層介質(zhì)格林函數(shù)。分層介質(zhì)格林函數(shù)每個分量都是一個復(fù)雜的索末菲積分，傳統(tǒng)處理方式采用離散復(fù)鏡像方法，但復(fù)鏡像的引入會導(dǎo)致加法定理的收斂性和穩(wěn)定性變差。因此，采用一種分層級處理方法，將半空間場作用劃分為3個層級進(jìn)行處理：第1層級對近相互作用，利用嚴(yán)格的半空間并矢格林函數(shù)計算；第2層級對遠(yuǎn)相互作用中的直射場作用，利用多層快速多極子“聚合、轉(zhuǎn)移、配置”操作計算；第3層級對遠(yuǎn)相互作用中的反射場作用，利用實鏡像方法計算。相比復(fù)鏡像方法引入多個復(fù)鏡像，實鏡像方法不但能夠保證算法的穩(wěn)定性和收斂性，同時僅引入一個實鏡像，有效地降低了計算量和存儲量。

1 算法理論

1.1 半空間多層快速多極子算法

半空間多層快速多極子可以模擬無限大平坦地面、海面等環(huán)境中目標(biāo)的電磁特性。與自由空間多層快速多極子相比，其難點在于分層介質(zhì)格林函數(shù)是并矢形式，其每個分量都是一個復(fù)雜的索末菲積分[15]。傳統(tǒng)方法中會引入復(fù)鏡像來求解索末菲積分，但復(fù)鏡像會給算法帶來穩(wěn)定性及收斂性降低的風(fēng)險[16]。因此在保證精度的前提下，將基函數(shù)和權(quán)函數(shù)之間的作用劃分為近相互作用和遠(yuǎn)相互作用，并將遠(yuǎn)相互作用再次細(xì)分為直射場作用及反射場作用。遠(yuǎn)相互作用是指源點和場點之間不經(jīng)過反射面的直接相互作用與源點和場點之間通過半空間分界面處的反射作用。圖1給出半空間中源點和場點之間的遠(yuǎn)相互作用示意。

圖1 半空間中源點和場點遠(yuǎn)相互作用示意圖

對半空間電大尺寸目標(biāo)進(jìn)行分析，阻抗形式如下：

Z=Znear+Zfar，i+Zfar，

(1)

其中，Znear是一個稀疏矩陣，代表近相互作用，通過矩量法直接填充，采用離散復(fù)鏡像技術(shù)(Discrete Complex Image Technology，DCIT)計算格林函數(shù)；Zfar，i代表直接遠(yuǎn)相互作用，其中i代表第i層的媒質(zhì)參數(shù)；Zfar代表具有分界面影響的遠(yuǎn)相互作用，下劃線“_”表示分界面的影響作用。分界面影響的遠(yuǎn)相互作用矩陣元素為

(2)

(3)

(4)

因此，半空間多層快速多極子算法僅僅只需要在自由空間多層快速多極子算法上額外計算式(4)中的實鏡像轉(zhuǎn)移因子與實鏡像聚合因子，不需要引入額外的計算量。

1.2 基于八叉樹的區(qū)域分解算法

本節(jié)給出一種新穎的基于八叉樹結(jié)構(gòu)的區(qū)域分解方法。與傳統(tǒng)區(qū)域分解算法相比，該方法最大的特點是無須創(chuàng)建人工交界面，并且不需要在子區(qū)域交界處強加電流連續(xù)性條件，在節(jié)省更多計算資源的同時具有較高的計算精度。

通過分析多層快速多極子算法特點可以看出，多層快速多極子是一種通過將基函數(shù)在空間內(nèi)分層分組，從而將基函數(shù)之間耦合作用轉(zhuǎn)化為非空組之間耦合作用的一種算法。其核心思路與區(qū)域分解算法吻合。基于此理念，筆者提出的OT-DDM方法并不需要對模型進(jìn)行額外的區(qū)域劃分處理，而是利用多層快速多極子算法在空間中形成的八叉樹結(jié)構(gòu)，自適應(yīng)地對未知量進(jìn)行分組來實現(xiàn)區(qū)域的劃分，從而大大降低了模型處理的工作量。

以一個PEC物體為例，假設(shè)目標(biāo)的整個表面S利用多層快速多極子八叉樹分成了兩個開放的子區(qū)域S1和S2，如圖2所示。

圖2 區(qū)域分解交界面示意圖

(5)

各個子區(qū)域上基函數(shù)個數(shù)與原始模型上的基函數(shù)個數(shù)滿足

N1+N2=N。

(6)

矩陣方程ZI=V可以寫成如下形式：

(7)

其中，Xi為第i個子區(qū)域Si的電磁流系數(shù)矩陣，Vi為第i個區(qū)域Si的電壓矩陣。當(dāng)i=j時，Zij為第i個子區(qū)域Si的自阻抗矩陣；當(dāng)i≠j時，Zij為第i個區(qū)域Si與第j個區(qū)域Sj的互阻抗矩陣。

OT-DDM采用基于 Krylov 子空間的迭代法，對該方程進(jìn)行求解，選用GMRES迭代求解器。與整體求解不同的是，在迭代過程中，每輪迭代時，GMRES求解器所給出的解都會根據(jù)區(qū)域數(shù)進(jìn)行拆分，只保留處理區(qū)域所需的那一部分電流，從而將求解整個矩陣向量乘問題轉(zhuǎn)化為處理各個子區(qū)域矩陣向量乘問題。

OT-DDM方法對整體模型進(jìn)行區(qū)域劃分之后，相鄰區(qū)域之間會存在一層鋸齒狀重合區(qū)域，如圖3所示。

圖3 特殊阻抗元素計算示意圖

為了保證區(qū)域間電流連續(xù)，需要對該部分區(qū)域所產(chǎn)生的奇異積分進(jìn)行特殊處理。采用部分阻抗元素來修正通過對近場檢驗而得到的ΔBi(Eij，Hij)，則有

(8)

(9)

(10)

由此在計算互阻抗元素時，即可將格林函數(shù)的梯度運算轉(zhuǎn)化為檢驗函數(shù)上的散度計算，從而將難以處理的高階奇異性轉(zhuǎn)化為能夠解析處理的低階奇異性，降低了計算復(fù)雜度，保證了計算的精度。

2 計算平臺

筆者使用了兩種計算平臺。

平臺I：西電曙光高性能集群，由16個計算節(jié)點組成，每個節(jié)點配置4顆18核的Intel(R)Xeon(R)Gold 6140 CPU @2.30 GHz處理器，以及512 GB RAM和3 TB磁盤。

平臺II：DELL工作站，該工作站配置2顆6核64位Intel Xeon E5-2620 2.0 GHz CPU、512 GB RAM和12 TB磁盤。

3 數(shù)值算例

首先通過一個簡單算例驗證OT-DDM的準(zhǔn)確性、有效性及可靠性。在驗證算例中，使用NASA杏仁體作為目標(biāo)，與商業(yè)軟件FEKO的多層快速多極子整體解對比，證明所提算法的可靠性；對比分析了OT-DDM與IE-DDM的計算結(jié)果，證明了筆者提出的方法對低散射體計算效果更好。再通過一個500波長艦船案例分析算法的資源消耗，表明本算法相對整體多層快速多極子內(nèi)存需求更低。最后，以一個海面上電大尺寸艦船計算問題來證明本算法的可擴(kuò)展性。

3.1 半空間杏仁體散射特性

以一個0.96 m×0.34 m×0.04 m的低雷達(dá)散射截面目標(biāo)的雙站雷達(dá)散射截面為例，說明OT-DDM具有較高的精度。此處，垂直極化的平面波沿θinc=60°，φinc=0°方向入射，頻率為600 MHz，杏仁體距離地面(相對介電常數(shù)和磁導(dǎo)率分別為εr=6.0-j1.0，μr=1.0)1個波長(0.5 m)，如圖4所示。

圖4 半空間杏仁體仿真模型

為了表明OT-DDM具有較高的數(shù)值精度，采用IE-NDDM、OT-DDM以及商業(yè)軟件FEKO進(jìn)行計算，收斂精度設(shè)置均為0.005，CFIE取0.6。同時，為了說明半空間環(huán)境對目標(biāo)散射特性的影響，在此將自由空間仿真結(jié)果也一同給出，其對比結(jié)果如圖5所示。由圖5可見，OT-DDM與多層快速多極子的計算結(jié)果在-30 dB以下仍然吻合良好，驗證了OT-DDM算法的準(zhǔn)確性與可靠性。

(a)θscat=89°，0°≤ φscat ≤360°平面雙站雷達(dá)散射截面計算結(jié)果

3.2 500波長半空間艦船雷達(dá)散射截面計算

以一個150 m×16.5 m×25.9 m艦船作為案例分析。當(dāng)平面波頻率為1 GHz時，模型對應(yīng)的電尺寸為500λ×55λ×86.25λ。其總網(wǎng)格量為11 188 432，總未知量為16 782 648。利用八叉樹算法劃分區(qū)域，選第3層，得到9個子區(qū)域。其模型與區(qū)域劃分如圖6(a)和圖6(b)所示。入射波沿φ=0°，θ=60°方向入射，極化方式為垂直極化，觀察角為-89°≤θscat≤ 89°，φscat=0°。上半空間為空氣，下半空間媒質(zhì)參數(shù)εr=80.0，σ=1.0 S/m。使用OT-DDM和MLFMA進(jìn)行對比計算，收斂精度均設(shè)為0.01，CFIE取0.6。MLFMA采用平臺I進(jìn)行計算，使用16個節(jié)點，每個節(jié)點使用12核，共計192核。OT-DDM采用平臺II進(jìn)行計算，使用24核。計算結(jié)果如圖6(c)所示?？梢钥闯?，兩種算法計算結(jié)果吻合良好。其資源消耗如表1所示，可以看出OT-DDM相對于MLFMA只消耗了約17.5%的內(nèi)存即可完成計算。

(a)整體模型

表1 各區(qū)域計算資源消耗表

由表1可以看到，區(qū)域分解最大內(nèi)存消耗為10 190 MB，而整體解則需要58 283 MB。此例中對艦船的模型劃分為9個區(qū)域，當(dāng)區(qū)域劃分更多時，則所需內(nèi)存消耗會進(jìn)一步降低。

3.3 千波長半空間艦船雷達(dá)散射截面計算

本例使用與上一算例中相同的艦船，將頻率提升為2 GHz，其余參數(shù)均不改變，其整體電尺寸大小為千波長，此時整體解在平臺I中已經(jīng)無法進(jìn)行計算。因此使用OT-DDM算法在平臺I中用9個節(jié)點，每個節(jié)點使用16核進(jìn)行計算，計算結(jié)果如圖7所示。

圖7 -89°≤θscat≤89°，φscat=0°平面的雙站雷達(dá)散射截面計算結(jié)果

4 總 結(jié)

為了在有限資源下解決半空間電大尺寸目標(biāo)雷達(dá)散射截面的散射計算問題，筆者提出一種OT-DDM算法。該算法利用多層快速多極子形成的八叉樹結(jié)構(gòu)進(jìn)行自適應(yīng)的區(qū)域劃分，避免了人工交界面的構(gòu)造，利用區(qū)域邊界上1/4阻抗來保證區(qū)域間電流連續(xù)，使得結(jié)果更加精確。筆者引入實鏡像源來計算半空間遠(yuǎn)相互作用，保證算法具有良好的收斂性和穩(wěn)定性。

筆者提出的算法具有3個優(yōu)點：(1)模型自動劃分，無須創(chuàng)建人工交界面；(2)有效地消除了區(qū)域邊界間的反射波，保證了區(qū)域間電流連續(xù)，使得計算精度與整體解基本一致；(3)在處理半空間問題時，沒有引入其他未知量，保證了計算資源消耗最少。

通過與商業(yè)軟件FEKO整體解進(jìn)行對比，證明了所提出算法的正確性；與傳統(tǒng)IE-NDDM進(jìn)行對比，證明了本算法具有較高的精度。同時，對比了筆者提出的算法與多層快速多極子整體解所需計算資源，證明本算法適宜在有限計算資源的條件下，求解半空間電大尺寸目標(biāo)雷達(dá)散射截面的散射問題。最后，以一個千波長半空間艦船雷達(dá)散射截面計算為例，驗證了所提出算法的計算能力。OT-DDM算法打破了對半空間電磁散射問題計算的電尺寸限制，能夠在有限資源的情況下計算更大規(guī)模的問題，為高頻段電大尺寸半空間雷達(dá)散射截面分析提供了更多可能與思路。