馬 驥，王浩放，馬蓓麗，張 帥

(1.中國電子科技集團(tuán)公司第二十七研究所，鄭州450047；2.西安電子科技大學(xué)，西安710071)

周期結(jié)構(gòu)(periodic structure，PS)最早是由美國科學(xué)家Rittenhouse通過觀察發(fā)絲制成的等間距柵格對日光的衍射現(xiàn)象提出的[1]。經(jīng)過半個多世紀(jì)的研究，周期結(jié)構(gòu)的頻率/空間濾波特性在軍事和民用的多個方面得到了廣泛應(yīng)用，對該結(jié)構(gòu)的研究也遍及微波、紅外及可見光等多個波段[2]。

如今，周期結(jié)構(gòu)廣泛應(yīng)用于微波和光學(xué)工程中，如柵格(grating)、頻率選擇表面(FSS)及光子晶體(photonic crystals，PC)等。以FSS為代表的周期結(jié)構(gòu)主要應(yīng)用于各種機(jī)/彈/艦載雷達(dá)艙的隱身設(shè)計(jì)。雷達(dá)艙是各類武器系統(tǒng)的主要強(qiáng)散射源，對隱身性能要求較高的飛機(jī)和其他武器雷達(dá)系統(tǒng)來說，采用FSS雷達(dá)罩是實(shí)現(xiàn)天線艙雷達(dá)隱身的有效手段，因此研究其電磁特性非常重要。目前，在計(jì)算電磁學(xué)領(lǐng)域，處理周期結(jié)構(gòu)最常用的方法是把整個結(jié)構(gòu)看成無限周期，根據(jù)Floquet定理[3]或周期Green函數(shù)[4]，將計(jì)算區(qū)域限制在一個單元上，可極大地簡化問題。然而，在工程實(shí)際中，并不存在無限周期結(jié)構(gòu)，所有周期結(jié)構(gòu)的尺寸都是有限的，所以，采用Floquet定理無法準(zhǔn)確反映單元之間的互耦及邊緣效應(yīng)。為更準(zhǔn)確地研究有限周期結(jié)構(gòu)，需同時考慮所有單元，這往往面臨電大尺寸有限周期問題。

矩量法(method of moments，MoM)[5]作為一種嚴(yán)格的積分方程方法，在電磁輻射和散射等電磁特性問題研究中的應(yīng)用較廣泛。由于MoM方法需要求解滿秩的線性方程組，當(dāng)未知量個數(shù)為N時，傳統(tǒng)MoM的計(jì)算復(fù)雜度為O(N2) ，直接求解和迭代求解的計(jì)算復(fù)雜度分別為O(N3)和O(N2)。對于現(xiàn)代電磁工程中處理的復(fù)雜周期結(jié)構(gòu)，由于阻抗矩陣條件數(shù)變差，基于迭代算法的MoM方法收斂速度太慢。為解決這一問題，近年來發(fā)展了基于區(qū)域分解的矩量法(subdomain method of moments, SMoM )[6-7]。SMoM方法將整個目標(biāo)分成若干個區(qū)域，每個區(qū)域可獨(dú)立求解。然而，受計(jì)算機(jī)資源的限制，MoM及SMoM在分析電大尺寸目標(biāo)體散射問題時計(jì)算效率很低，對于超大型電磁問題的求解更是無從談起。

為降低計(jì)算復(fù)雜度并減少內(nèi)存需求量，近年來發(fā)展了許多高效算法分析有限大周期結(jié)構(gòu)的電磁散射特性。第1類是基于MoM方法的快速算法，如多層快速多極子(multilevel fast multipole method, MLFMM)[8]、預(yù)修正快速傅里葉變換法(precorrected FFT, P-FFT)[9]、積分方程快速傅里葉變換法(IE-FFT)[10]及自適應(yīng)積分方法(adaptive integral method，AIM)[11]等，此類算法將內(nèi)存需求和計(jì)算復(fù)雜度降低到O(NlogN)。第2類是采用高階基函數(shù)，如宏基函數(shù)(macro basis function，MBF)[12]、特征基函數(shù)法(characteristic basis function method, CBFM)[13]和子全域基函數(shù)法(sub-entire-domain，SED)[14]等，利用這些基函數(shù)可極大減少未知量的數(shù)目。第3類方法是子陣方法[15]，利用提取小型周期結(jié)構(gòu)陣列中各單元的表面電流，來等效更大規(guī)模周期結(jié)構(gòu)陣列中相似陣列環(huán)境下的單元表面電流，進(jìn)而計(jì)算該大型周期陣列的散射場，在有效減小計(jì)算量的情況下獲得較高的精度。上述算法中，子陣方法具有計(jì)算量低和計(jì)算效率高的優(yōu)勢，被廣泛應(yīng)用于陣列天線的輻射和散射特性分析中。

本文通過基于子陣策略的區(qū)域分解矩量法(SMoM-Subarray)來分析有限大周期結(jié)構(gòu)的電磁散射問題。首先，選取一個合適的子陣，通過引入?yún)^(qū)域分解法(domain decomposition method, DDM)[16-17]的基本概念，將整個子陣分成若干個子區(qū)域；其次，利用三角形離散子區(qū)域表面并輸出剖分后的子區(qū)域，考慮子區(qū)域之間的耦合，利用MoM方法獲得該部分子區(qū)域表面的感應(yīng)電流，來等效整個周期結(jié)構(gòu)的表面電流；最后，利用整個周期結(jié)構(gòu)的表面感應(yīng)電流計(jì)算遠(yuǎn)區(qū)散射場，進(jìn)而可獲得目標(biāo)的雷達(dá)散射截面(radar cross section, RCS)。

1理論分析

1.1電場積分方程(electric field integral equation, EFIE)

圖1為整個目標(biāo)分解為若干個子區(qū)域的示意圖。圖中，Ω為任意形狀的入射平面波Ei所照射的導(dǎo)體。

圖1 整個目標(biāo)Ω分解為Q個子區(qū)域示意圖Fig.1 Diagram of decomposing target domain Ω into Q sub domain

將整個目標(biāo)分為Q個子區(qū)域，每個子區(qū)域上的EFIE可表示為

(1)

q=1,2,…,Q

(2)

其中，k為頻率f的波數(shù)；η0為自由空間本征阻抗；Jq(r)為子區(qū)域ΩQ上的等效電流；r為球坐標(biāo)矩陣；G(r,r′)是自由空間格林函數(shù)，定義為

(3)

1.2SMoM

對于EFIE的數(shù)值解，所有的子區(qū)域被離散為許多小三角形塊，每個子區(qū)域上的等效電流Jq(r)可用Rao-Wilton-Glisson (RWG)基函數(shù)[18]fk(r)來展開為

(4)

其中，N1，N2，…，NQ-1，NQ為基函數(shù)的個數(shù)；I1,k，I2,k，…，IQ-1,k，IQ,k為Q個子區(qū)域的當(dāng)前系數(shù)。將式(4)代入式(1)，應(yīng)用Galerkin方法[19]，得到每個子區(qū)域的結(jié)果為

(5)

其中，I1，I2，…，IQ-1，IQ為列向量N1×1,N2×1,…，NQ-1×1,NQ×1在Q個子區(qū)域上基函數(shù)的未知幅度；Z1,Z2,…，ZQ-1,ZQ為阻抗矩陣；V1,V2,…，VQ-1，VQ為激勵矢量；ΔV1，ΔV2，…，ΔVQ-1，ΔVQ為激勵矢量的變化量。它們可表示為

Zq,m,k=〈fq,m(r),L(fq,k(r))〉

(6)

(7)

(8)

1.3SMoM-Subarray

對于大規(guī)模的周期性結(jié)構(gòu)，如200×200陣列，需填充40 000個矩陣元素并求解上表面電流，而SMoM程序難以在單個工作站運(yùn)行。為解決這個問題，本文引入子陣技術(shù)，基于子陣策略，提出SMoM-Subarray方法，不僅可大幅降低內(nèi)存，還可減少總計(jì)算時間。

在不失一般性的情況下，我們使用3×3子陣列的表面電流來推導(dǎo)出大的5×5陣列的表面電流。圖2為由3×3子陣列推導(dǎo)5×5陣列等效表面電流示意圖。

圖2 由3×3子陣列推導(dǎo)5×5陣列等效表面電流示意圖Fig.2 Equivalent surface current array 5×5 obtained from subarray 3×3

圖3為本文提出的SMoM-Subarray算法流程圖，程序詳細(xì)執(zhí)行步驟為：

(9)

(10)

(11)

4)檢查Q個子區(qū)域的最大電流偏差是否小于規(guī)定的偏差

(12)

將該過程從式(10)到式(12)進(jìn)行T次迭代，直至達(dá)到期望的精度。

沿x軸排列定義為行，用p表示；沿y軸排列定義為列，用q表示；3×3元子陣的各子區(qū)域(單元)電流用In,sub來表示；5×5元平面陣各子區(qū)域(單元)的電流用Iβ表示，其中，β=1,2,3,…，24，25,β=(p-1)×5+q。各單元電流的等效建立過程為：

圖3 SMoM-Subarray算法流程圖Fig.3 Flow chart of SMoM-Subarray

1)第1行各單元的電流為

(13)

2)第2至第4行各單元的電流為

(14)

3)第5行各單元的電流為

(15)

其中，dx，dy分別為沿x軸和y軸方向單元的間距；θ和φ為原始陣列單元對應(yīng)的球坐標(biāo)分量

利用整個有限大周期結(jié)構(gòu)的表面電流，可計(jì)算得到遠(yuǎn)區(qū)散射場，并進(jìn)而求得有限大周期結(jié)構(gòu)的RCS。

2數(shù)值結(jié)果與分析

為驗(yàn)證SMoM-Subarray的準(zhǔn)確性和有效性，考慮圓金屬與金屬十字陣列2個算例，所有算例都在主頻為2.4 GHz和內(nèi)存為16 G的個人計(jì)算機(jī)上完成。

算例1 分別采用MoM，SMoM，SMoM-Subarray方法計(jì)算如圖4所示的10×10圓金屬陣列在平面波照射下的RCS。

圖4 10×10圓金屬陣列 Fig.410×10 round metal array

其中，圓單元的直徑為0.1 m；相鄰2個圓單元的間距為d=0.5λ0；λ0為自由空間中的波長。陣列采用三角形面片進(jìn)行剖分后，三角形數(shù)目為3 200；總未知量數(shù)目為4 100。選取7×7圓金屬陣列為子陣，將該子陣分為49個子區(qū)域，子陣采用三角形面片進(jìn)行剖分后，三角形數(shù)目為1 568，未知量數(shù)目為2 009，入射波頻率為1 GHz，入射角θ變化范圍為-90°～90°，極化方式分別為θ極化和φ極化。圖5為10×10圓金屬陣列單站RCS計(jì)算結(jié)果。由圖5可見，3種方法針對2種極化入射波的RCS計(jì)算結(jié)果吻合得很好，證明了本文提出的SMoM-Subarray方法的正確性。

(a)θ polarized incident wave

(b)φ polarized incident wave

算例2分別采用MoM，SMoM，SMoM-Subarray方法計(jì)算如圖6所示的20×20金屬十字陣列在平面波照射下的RCS。其中，十字單元的尺寸為1.6 cm×0.4 cm；相鄰2個十字單元的間距為d= 2/3λ0。陣列采用三角形面片進(jìn)行剖分后，三角形數(shù)目為20 000，總未知量數(shù)目為23 600。選取7×7金屬十字陣列為子陣，將該子陣劃分為49個子區(qū)域，子陣采用三角形面片進(jìn)行剖分后，三角形總數(shù)為2 450，未知量數(shù)目為2 891。平面入射波頻率為10 GHz；入射角θ變化范圍為-90°～90°；極化方式分別為θ極化和φ極化。圖7為20×20金屬十字陣列單站RCS計(jì)算結(jié)果。由圖7可見，3種方法針對2種極化入射波的RCS計(jì)算結(jié)果吻合得很好，同樣證明了本文所提出算法的正確性。

圖6 20×20金屬十字陣列 Fig.620×20 cross-shaped metal array

(a)θ polarized incident wave

(b)φ polarized incident wave

表1列出了2個陣列模型通過MoM，SMoM，SMoM-Subarray方法的未知量數(shù)目、計(jì)算時間及所消耗的內(nèi)存。

表1 不同計(jì)算方法的比較Tab.1 Comparison of different calculation

由表1可知，本文所提出的SMoM-Subarray方法不僅大大降低了內(nèi)存消耗，同時在不失精度的情況下分別使圓金屬陣列和十字金屬陣列的計(jì)算時間減少了88.2%和90.6%，證明了該方法的高效性。

3結(jié)論

本文提出了一種快速有效分析周期結(jié)構(gòu)的SMoM-Subarray方法，該方法將區(qū)域分解法和子陣方法相結(jié)合，把求解區(qū)域分為若干個非重疊的子陣，并利用子陣結(jié)構(gòu)的幾何重復(fù)性，只需分析其中的某個子陣，大大提高了計(jì)算效率。通過圓金屬與金屬十字陣列2個算例驗(yàn)證了SMoM-Subarray方法的正確性。同時計(jì)算給出，與MoM、SMoM方法相比，SMoM-Subarray方法在不失精度的情況下計(jì)算時間減少了88.2%和90.6%，證明了該方法的高效性。