基于子阵策略和区域分解矩量法的电磁散射分析方法

马 骥，王浩放，马蓓丽，张 帅

(1.中国电子科技集团公司第二十七研究所，郑州450047；2.西安电子科技大学，西安710071)

周期结构(periodic structure，PS)最早是由美国科学家Rittenhouse通过观察发丝制成的等间距栅格对日光的衍射现象提出的[1]。经过半个多世纪的研究，周期结构的频率/空间滤波特性在军事和民用的多个方面得到了广泛应用，对该结构的研究也遍及微波、红外及可见光等多个波段[2]。

如今，周期结构广泛应用于微波和光学工程中，如栅格(grating)、频率选择表面(FSS)及光子晶体(photonic crystals，PC)等。以FSS为代表的周期结构主要应用于各种机/弹/舰载雷达舱的隐身设计。雷达舱是各类武器系统的主要强散射源，对隐身性能要求较高的飞机和其他武器雷达系统来说，采用FSS雷达罩是实现天线舱雷达隐身的有效手段，因此研究其电磁特性非常重要。目前，在计算电磁学领域，处理周期结构最常用的方法是把整个结构看成无限周期，根据Floquet定理[3]或周期Green函数[4]，将计算区域限制在一个单元上，可极大地简化问题。然而，在工程实际中，并不存在无限周期结构，所有周期结构的尺寸都是有限的，所以，采用Floquet定理无法准确反映单元之间的互耦及边缘效应。为更准确地研究有限周期结构，需同时考虑所有单元，这往往面临电大尺寸有限周期问题。

矩量法(method of moments，MoM)[5]作为一种严格的积分方程方法，在电磁辐射和散射等电磁特性问题研究中的应用较广泛。由于MoM方法需要求解满秩的线性方程组，当未知量个数为N时，传统MoM的计算复杂度为O(N2) ，直接求解和迭代求解的计算复杂度分别为O(N3)和O(N2)。对于现代电磁工程中处理的复杂周期结构，由于阻抗矩阵条件数变差，基于迭代算法的MoM方法收敛速度太慢。为解决这一问题，近年来发展了基于区域分解的矩量法(subdomain method of moments, SMoM )[6-7]。SMoM方法将整个目标分成若干个区域，每个区域可独立求解。然而，受计算机资源的限制，MoM及SMoM在分析电大尺寸目标体散射问题时计算效率很低，对于超大型电磁问题的求解更是无从谈起。

为降低计算复杂度并减少内存需求量，近年来发展了许多高效算法分析有限大周期结构的电磁散射特性。第1类是基于MoM方法的快速算法，如多层快速多极子(multilevel fast multipole method, MLFMM)[8]、预修正快速傅里叶变换法(precorrected FFT, P-FFT)[9]、积分方程快速傅里叶变换法(IE-FFT)[10]及自适应积分方法(adaptive integral method，AIM)[11]等，此类算法将内存需求和计算复杂度降低到O(NlogN)。第2类是采用高阶基函数，如宏基函数(macro basis function，MBF)[12]、特征基函数法(characteristic basis function method, CBFM)[13]和子全域基函数法(sub-entire-domain，SED)[14]等，利用这些基函数可极大减少未知量的数目。第3类方法是子阵方法[15]，利用提取小型周期结构阵列中各单元的表面电流，来等效更大规模周期结构阵列中相似阵列环境下的单元表面电流，进而计算该大型周期阵列的散射场，在有效减小计算量的情况下获得较高的精度。上述算法中，子阵方法具有计算量低和计算效率高的优势，被广泛应用于阵列天线的辐射和散射特性分析中。

本文通过基于子阵策略的区域分解矩量法(SMoM-Subarray)来分析有限大周期结构的电磁散射问题。首先，选取一个合适的子阵，通过引入区域分解法(domain decomposition method, DDM)[16-17]的基本概念，将整个子阵分成若干个子区域；其次，利用三角形离散子区域表面并输出剖分后的子区域，考虑子区域之间的耦合，利用MoM方法获得该部分子区域表面的感应电流，来等效整个周期结构的表面电流；最后，利用整个周期结构的表面感应电流计算远区散射场，进而可获得目标的雷达散射截面(radar cross section, RCS)。

1理论分析

1.1电场积分方程(electric field integral equation, EFIE)

图1为整个目标分解为若干个子区域的示意图。图中，Ω为任意形状的入射平面波Ei所照射的导体。

图1 整个目标Ω分解为Q个子区域示意图Fig.1 Diagram of decomposing target domain Ω into Q sub domain

将整个目标分为Q个子区域，每个子区域上的EFIE可表示为

(1)

q=1,2,…,Q

(2)

其中，k为频率f的波数；η0为自由空间本征阻抗；Jq(r)为子区域ΩQ上的等效电流；r为球坐标矩阵；G(r,r′)是自由空间格林函数，定义为

(3)

1.2SMoM

对于EFIE的数值解，所有的子区域被离散为许多小三角形块，每个子区域上的等效电流Jq(r)可用Rao-Wilton-Glisson (RWG)基函数[18]fk(r)来展开为

(4)

其中，N1，N2，…，NQ-1，NQ为基函数的个数；I1,k，I2,k，…，IQ-1,k，IQ,k为Q个子区域的当前系数。将式(4)代入式(1)，应用Galerkin方法[19]，得到每个子区域的结果为

(5)

其中，I1，I2，…，IQ-1，IQ为列向量N1×1,N2×1,…，NQ-1×1,NQ×1在Q个子区域上基函数的未知幅度；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

对于大规模的周期性结构，如200×200阵列，需填充40 000个矩阵元素并求解上表面电流，而SMoM程序难以在单个工作站运行。为解决这个问题，本文引入子阵技术，基于子阵策略，提出SMoM-Subarray方法，不仅可大幅降低内存，还可减少总计算时间。

在不失一般性的情况下，我们使用3×3子阵列的表面电流来推导出大的5×5阵列的表面电流。图2为由3×3子阵列推导5×5阵列等效表面电流示意图。

图2 由3×3子阵列推导5×5阵列等效表面电流示意图Fig.2 Equivalent surface current array 5×5 obtained from subarray 3×3

图3为本文提出的SMoM-Subarray算法流程图，程序详细执行步骤为：

(9)

(10)

(11)

4)检查Q个子区域的最大电流偏差是否小于规定的偏差

(12)

将该过程从式(10)到式(12)进行T次迭代，直至达到期望的精度。

沿x轴排列定义为行，用p表示；沿y轴排列定义为列，用q表示；3×3元子阵的各子区域(单元)电流用In,sub来表示；5×5元平面阵各子区域(单元)的电流用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轴方向单元的间距；θ和φ为原始阵列单元对应的球坐标分量

利用整个有限大周期结构的表面电流，可计算得到远区散射场，并进而求得有限大周期结构的RCS。

2数值结果与分析

为验证SMoM-Subarray的准确性和有效性，考虑圆金属与金属十字阵列2个算例，所有算例都在主频为2.4 GHz和内存为16 G的个人计算机上完成。

算例1 分别采用MoM，SMoM，SMoM-Subarray方法计算如图4所示的10×10圆金属阵列在平面波照射下的RCS。

图4 10×10圆金属阵列 Fig.410×10 round metal array

其中，圆单元的直径为0.1 m；相邻2个圆单元的间距为d=0.5λ0；λ0为自由空间中的波长。阵列采用三角形面片进行剖分后，三角形数目为3 200；总未知量数目为4 100。选取7×7圆金属阵列为子阵，将该子阵分为49个子区域，子阵采用三角形面片进行剖分后，三角形数目为1 568，未知量数目为2 009，入射波频率为1 GHz，入射角θ变化范围为-90°～90°，极化方式分别为θ极化和φ极化。图5为10×10圆金属阵列单站RCS计算结果。由图5可见，3种方法针对2种极化入射波的RCS计算结果吻合得很好，证明了本文提出的SMoM-Subarray方法的正确性。

(a)θ polarized incident wave

(b)φ polarized incident wave

算例2分别采用MoM，SMoM，SMoM-Subarray方法计算如图6所示的20×20金属十字阵列在平面波照射下的RCS。其中，十字单元的尺寸为1.6 cm×0.4 cm；相邻2个十字单元的间距为d= 2/3λ0。阵列采用三角形面片进行剖分后，三角形数目为20 000，总未知量数目为23 600。选取7×7金属十字阵列为子阵，将该子阵划分为49个子区域，子阵采用三角形面片进行剖分后，三角形总数为2 450，未知量数目为2 891。平面入射波频率为10 GHz；入射角θ变化范围为-90°～90°；极化方式分别为θ极化和φ极化。图7为20×20金属十字阵列单站RCS计算结果。由图7可见，3种方法针对2种极化入射波的RCS计算结果吻合得很好，同样证明了本文所提出算法的正确性。

图6 20×20金属十字阵列 Fig.620×20 cross-shaped metal array

(a)θ polarized incident wave

(b)φ polarized incident wave

表1列出了2个阵列模型通过MoM，SMoM，SMoM-Subarray方法的未知量数目、计算时间及所消耗的内存。

表1 不同计算方法的比较Tab.1 Comparison of different calculation

由表1可知，本文所提出的SMoM-Subarray方法不仅大大降低了内存消耗，同时在不失精度的情况下分别使圆金属阵列和十字金属阵列的计算时间减少了88.2%和90.6%，证明了该方法的高效性。

3结论

本文提出了一种快速有效分析周期结构的SMoM-Subarray方法，该方法将区域分解法和子阵方法相结合，把求解区域分为若干个非重叠的子阵，并利用子阵结构的几何重复性，只需分析其中的某个子阵，大大提高了计算效率。通过圆金属与金属十字阵列2个算例验证了SMoM-Subarray方法的正确性。同时计算给出，与MoM、SMoM方法相比，SMoM-Subarray方法在不失精度的情况下计算时间减少了88.2%和90.6%，证明了该方法的高效性。