Fast solution for electromagnetic scattering of dynamic group targets under sea background based on a hybrid method of CM and PO
-
摘要:
群目标广泛存在于电磁工程中，对其开发高效、准确的仿真方法具有重要意义。针对动态群目标在各个时刻下单元排布可能不同的特点，首先提出了按照单元分组的快速多极子算法，使用特征模方法降低群目标未知量，采用预存储的策略来实现动态群目标的高效、准确分析，该方法称为动态群目标的特征模快速分析方法；然后，考虑海面背景对群目标电磁散射的影响，将以上提出的动态群目标的特征模快速分析方法与高频的物理光学法进行结合，提出了一种高低频混合方法。数值算例表明，本文提出的高低频混合方法能够在保证较高计算精度的同时实现动态群目标和海面背景整体的高效电磁特性分析，同时还实现了计算效率和精度之间的合理平衡。
-
关键词:
- 动态群目标 /
- 海面背景 /
- 高低频混合方法 /
- 特征模理论 (TCM) /
- 快速多极子算法 (FMA) /
- 物理光学 (PO) 法
Abstract:Group targets are widely present in electromagnetic engineering, and the development of efficient and accurate simulation methods for them is of great significance. Considering the characteristic that the unit arrangement of dynamic group targets may vary at different times, this paper first proposes a fast multipole algorithm (FMA) based on unit grouping. The method uses the theory of characteristic mode (TCM) to reduce the unknowns of group targets and uses a pre-storage strategy to achieve efficient and accurate analysis of dynamic group targets; Then, in order to consider the influence of sea background on electromagnetic scattering of group targets, a hybrid method of high and low frequency is proposed by combining the fast analysis method based on TCM of dynamic group targets proposed above with the high-frequency physical optics (PO) method. Numerical examples demonstrate that the proposed hybrid high-low frequency method can achieve efficient electromagnetic characteristic analysis of dynamic group targets and the sea background as a whole while ensuring high computational accuracy. This approach thereby achieves a reasonable balance between computational efficiency and precision.
-
0 引 言
在电磁工程中，很多常见目标，如无人机、舰船等，通常不是以单一目标的形式出现，而是以群目标的形式出现，并且它们的位置和姿态是随着时间动态变化的。现有的针对群目标的快速分析方法大多只能分析静止不动的群目标，很难兼顾外部环境的影响。
矩量法(method of moments, MoM)的计算精度高，然而它建立的系数矩阵为满阵，存在高内存需求和高复杂度的弊端。针对此问题，研究人员提出了诸如快速多极子算法(fast multipole algorithm, FMA)[1]或多层FMA(multilevel FMA, MLFMA)[2]、自适应积分方法(adaptive integration method, AIM)[3]和快速傅里叶变换(fast Fourier transform, FFT)方法[4]等一系列加速算法。
基于区域分解方法(domain decomposition method, DDM)[5]思想，研究人员提出了一类新的求解方法——全域基函数法[6]，其中比较典型的有特征基函数方法[7-8]、复合基函数方法[9]、子全域基函数方法[10]等。特征模式是由MoM阻抗矩阵建立的广义特征方程进行求解而得到的，因此其只与目标的结构、材料等固有属性相关，与外加的激励源无关[11]。因此，可作为全域基函数用于准确、高效地分析重复性结构的电磁散射或辐射问题[12-13]。
在动态群目标的快速分析方面，Zhang等提出了对每个单元采用局域八叉树分组的方式[14]，在此基础上，又进一步提出了使用全局八叉树分组和局域八叉树分组相结合的双八叉树结构对局域八叉树分组方法进行了改进[15]。Li等[16]提出了基于等效球面的区域分解方法用于高效分析运动的旋转对称体群目标的电磁散射和雷达成像。
针对海面这类电尺寸较大目标，通常可选择精度较低但非常高效的高频方法对其进行分析。高频类方法主要可分为两类：一类是以物理光学(physical optics, PO)法[17]、边缘等效电磁流法 (equivalent edge current method, EEC)[18]、物理绕射理论(physical theory of diffraction, PTD)[19]为代表的电流类方法；另一类是以几何光学(geometric optics, GO)法[17]、几何绕射理论(geometrical theory of diffraction, GTD)[20]、一致绕射理论(uniform theory of diffraction, UTD)[21]等为代表的射线类方法。
海面环境下复杂目标的电磁散射必然会受粗糙背景的影响，要解决这类问题通常可选择高低频混合方法。现有的高低频混合方法主要有两类，一类是基于射线混合的方法。该类方法是将高频的弹跳射线(shooting and bouncing ray, SBR)方法与低频的MoM结合形成的SBR-MoM混合方法[22]，这类方法计算效率高但是精度较低。另一类是基于电流混合的方法。该类方法是将高频的PO方法与低频的MoM结合形成的PO-MoM混合方法[23]，耗时较长但计算精度相对较高。本文在第二类高低频混合方法的基础上提出了一种新的混合方法，使用精确且高效的特征模方法来分析动态群目标，使用PO法来分析海面背景，利用迭代的方式来考虑两个区域之间的耦合作用，以求在不损失较多精度的前提下实现对整个场景的高效电磁分析。
1 动态群目标的特征模快速分析方法
1.1 群目标的特征模方法
假设待分析的群目标中共有P个重复性单元，应用MoM可建立如下矩阵方程[24]：
{\boldsymbol{ZI}} = {\boldsymbol{V}} (1) 式中： {\boldsymbol{Z}} 为MoM的系数矩阵； {\boldsymbol{I}} 为待求的解向量； {\boldsymbol{V}} 为与激励有关的右边向量。当基函数按照单元顺序排列时，系数矩阵和右边向量又可以表示为：
{\boldsymbol{Z}} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{Z}}_{11}}}&{{{\boldsymbol{Z}}_{12}}}& \cdots &{{{\boldsymbol{Z}}_{1P}}} \\ {{{\boldsymbol{Z}}_{21}}}&{{{\boldsymbol{Z}}_{22}}}& \cdots &{{{\boldsymbol{Z}}_{2P}}} \\ \vdots & \vdots &{}& \vdots \\ {{{\boldsymbol{Z}}_{P1}}}&{{{\boldsymbol{Z}}_{P2}}}& \cdots &{{{\boldsymbol{Z}}_{PP}}} \end{array}} \right]_{NP \times NP}} (2) {\boldsymbol{V}} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{V}}_1}} \\ {{{\boldsymbol{V}}_2}} \\ \vdots \\ {{{\boldsymbol{V}}_P}} \end{array}} \right]_{NP \times 1}} (3) 式中， N 为单个单元中的未知量个数。由于单元都是相同的，它们的特征模式也是相同的，所以这些特征模式可以通过求解广义特征值方程得到。当单元是纯金属结构时，广义特征值方程为[11]
{{\boldsymbol{X}}_{11}}{{\boldsymbol{J}}_n}{\text{ = }}{\lambda _n}{{\boldsymbol{R}}_{11}}{{\boldsymbol{J}}_n} (4) 式中： {{\boldsymbol{R}}_{11}} 和 {{\boldsymbol{X}}_{11}} 分别为矩阵 {{\boldsymbol{Z}}_{11}} 的实部和虚部； {\lambda _n} 和 {{\boldsymbol{J}}_n} 分别为特征值和特征模式。这些特征模式可以被用作全域基函数来实现矩阵方程的降阶，具体过程为[13]
{\boldsymbol{\tilde Z\tilde I}} = {\boldsymbol{\tilde V}} (5) 式中：
{\boldsymbol{\tilde Z}} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{Z}}_{11}}{\boldsymbol{B}}}&{{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{Z}}_{12}}{\boldsymbol{B}}}& \cdots &{{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{Z}}_{1P}}{\boldsymbol{B}}} \\ {{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{Z}}_{21}}{\boldsymbol{B}}}&{{{\boldsymbol{B}}^T}{{\boldsymbol{Z}}_{22}}{\boldsymbol{B}}}& \cdots &{{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{Z}}_{2P}}{\boldsymbol{B}}} \\ \vdots & \vdots &{}& \vdots \\ {{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{Z}}_{P1}}{\boldsymbol{B}}}&{{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{Z}}_{P2}}{\boldsymbol{B}}}& \cdots &{{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{Z}}_{PP}}{\boldsymbol{B}}} \end{array}} \right]_{KP \times KP}} (6) {\boldsymbol{\tilde V}} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{V}}_1}} \\ {{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{V}}_2}} \\ \vdots \\ {{{\boldsymbol{B}}^{\text{T}}}{{\boldsymbol{V}}_P}} \end{array}} \right]_{KP \times 1}} (7) {\boldsymbol{B}} 为由 K 个特征模式组成的矩阵，并且 K \ll N 。因此，原始矩阵方程的维度得到了大幅度的缩减，从而可以快速进行求解。当得到降阶的解向量 {\boldsymbol{\tilde I}} 后，就可以利用特征模式将其还原为原始的解向量，有
{\boldsymbol{I}} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{B}}{{{\boldsymbol{\tilde I}}}_1}} \\ {{\boldsymbol{B}}{{{\boldsymbol{\tilde I}}}_2}} \\ \vdots \\ {{\boldsymbol{B}}{{{\boldsymbol{\tilde I}}}_P}} \end{array}} \right] (8) 模式数K的选取原则为：针对单个单元 i 按照MoM计算出其感应电流 {I_{i,{\text{MoM}}}} ，按照式(5)和(8)计算其感应电流{I_{i,{\text{CM}}}}，以 {I_{i,{\text{MoM}}}} 为标准，二者相对误差小于5%对应的模式数就是最终计算需要的模式数K。在这个过程中采用ARPACK库中隐式重启阿诺尔德方法(implicitly restarted Arnoldi method, IRAM)计算式(4)提取特征模式。IRAM需要计算多次矩阵向量积，采用MLFMA加速这一过程可以将计算复杂度从 o({N^2}) 降到o(N\lg N)[25]。
从式(6)可以看出，阻抗矩阵 {{\boldsymbol{Z}}_{ij}} 与矩阵 {\boldsymbol{B}} 的乘积可以分解为若干矩阵矢量乘，从而可以被MLFMA加速计算。首先，假设单元i和单元j相互作用产生的降阶矩阵块 {{\boldsymbol{\tilde Z}}_{ij}} 的矩阵元素为
{{\boldsymbol{\tilde Z}}_{ij}} = \left[ {\begin{array}{*{20}{c}} {{c_{11}}}&{{c_{12}}}& \cdots &{{c_{1K}}} \\ {{c_{21}}}&{{c_{22}}}& \cdots &{{c_{2K}}} \\ \vdots & \vdots &{}& \vdots \\ {{c_{K1}}}&{{c_{K2}}}& \cdots &{{c_{KK}}} \end{array}} \right] (9) 在矩阵块降阶过程中 {{\boldsymbol{Z}}_{ij}} 右乘矩阵 {\boldsymbol{B}} 中的第 t 个向量的过程可以使用MLFMA加速，即[13]
{{\boldsymbol{Z}}_{ij}}{{\boldsymbol{J}}_t} = \left[ {\begin{array}{*{20}{c}} {{a_{11}}}&{{a_{12}}}& \cdots &{{a_{1N}}} \\ {{a_{21}}}&{{a_{22}}}& \cdots &{{a_{2N}}} \\ \vdots & \vdots &{}& \vdots \\ {{a_{N1}}}&{{a_{N2}}}& \cdots &{{a_{NN}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {b_1^{(t)}} \\ {b_2^{(t)}} \\ \vdots \\ {b_N^{(t)}} \end{array}} \right] (10) \begin{split} \sum\limits_{n = 1}^N {{a_{mn}}b_n^{(t)}} \approx & \sum\limits_{q \in {G_p}} {\sum\limits_{n \in {G_q}} {{a_{mn}}b_n^{(t)}} } \\ & + \mathop{{\int \int}\mkern-18mu \bigcirc} {{{\boldsymbol{R}}_{mp}}({\boldsymbol{\hat k}}) \cdot \sum\limits_{q \notin {G_p}} {{\varGamma _{pq}}({{{\boldsymbol{\hat r}}}_{pq}} \cdot {\boldsymbol{\hat k}})} } \sum\limits_{n \in {G_q}} {{{\boldsymbol{F}}_{qn}}({\boldsymbol{\hat k}})} b_n^{(t)}{{\text{d}}^2}{\boldsymbol{\hat k}} \\ =& {{\textit{z}}_{mt}} + {d_{mt}}\\[-1pt] \end{split} (11) 式中： {G_p} 为场基函数的所有近场组； {G_q} 为源基函数的所有近场组；\displaystyle \mathop{{\int \int}\mkern-18mu \bigcirc} {\; \cdot \;} {{\mathrm{d}}^2}{\boldsymbol{\hat k}} 为单位球面上的积分； {{\boldsymbol{R}}_{mp}} 为配置因子； {\varGamma _{pq}} 为转移因子； {{\boldsymbol{F}}_{qn}} 为聚合因子； {{\textit{z}}_{mt}} 为近场作用得到的部分； {d_{mt}} 为通过聚合、转移和配置作用得到的部分。关于MLFMA的更多细节可参考文献[13, 25-26]，这里不再赘述。然后，降阶矩阵元素 {c_{st}} 可以通过以下方式计算：
\begin{split} {c_{st}} =& \left[ {\begin{array}{*{20}{c}} {b_1^{(s)}}&{b_2^{(s)}}& \cdots &{b_N^{(s)}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\textit{z}}_{1t}} + {d_{1t}}} \\ {{{\textit{z}}_{2t}} + {d_{2t}}} \\ \vdots \\ {{{\textit{z}}_{Nt}} + {d_{Nt}}} \end{array}} \right] \\ =& (b_1^{(s)}{{\textit{z}}_{1t}} + b_2^{(s)}{{\textit{z}}_{2t}} + \cdots + b_N^{(s)}{{\textit{z}}_{Nt}}) \\ & + (b_1^{(s)}{d_{1t}} + b_2^{(s)}{d_{2t}} + \cdots + b_N^{(s)}{d_{Nt}}) \\ =& \sum\limits_{m = 1}^N {b_m^{(s)}\sum\limits_{q \in {G_p}} {\sum\limits_{n \in {G_q}} {{a_{mn}}b_n^{(t)}} } } \\ &+ \sum\limits_{m = 1}^N {b_m^{(s)}} \mathop{{\int \int}\mkern-18mu \bigcirc} {{{\boldsymbol{R}}_{mp}}({\boldsymbol{\hat k}}) \cdot \sum\limits_{q \notin {G_p}} {{\varGamma _{pq}}({{{\boldsymbol{\hat r}}}_{pq}} \cdot {\boldsymbol{\hat k}})} } \sum\limits_{n \in {G_q}} {{{\boldsymbol{F}}_{qn}}({\boldsymbol{\hat k}})} b_n^{(t)}{{\text{d}}^2}{\boldsymbol{\hat k}} \end{split} (12) 1.2 按照单元分组的FMA
在不同时刻下各单元的位置和姿态可能会发生变化，常规的全局八叉树分组无法保证每个时刻下的分组信息不变[14]，如图1所示。为解决这个问题，介绍一种新的按照单元分组的FMA(unit-grouping-based FMA, UGFMA)。如图2所示，在UGFMA中不再需要建立父层与子层之间繁琐的索引关系。分组以后第i个单元与第j个单元上的基函数相互作用就可以按照图3所示的路线进行。
![]()
图 1 动态群目标的全局八叉树分组示意图Fig. 1 Schematic of global octree grouping for dynamic group targets![]()
图 3 UGFMA中不同单元间基函数的作用路线示意图Fig. 3 Schematic of the action path of basis functions between different units in UGFMA各单元间通常都会保持一定的安全距离以避免相撞。因此，图3中的3个向量之间一般都会满足以下条件：
\left| {{{\boldsymbol{r}}_{pq}}} \right| > \left| {{{\boldsymbol{r}}_{mp}} + {{\boldsymbol{r}}_{qn}}} \right| (13) 因此，格林函数可以应用加法定理展开为
G({\boldsymbol{r}},{\boldsymbol{r'}}) = \frac{{{{\text{e}}^{ - {\text{j}}{k_0}\left| {{{\boldsymbol{r}}_m} - {{\boldsymbol{r}}_n}} \right|}}}}{{4{\text{π}}\left| {{{\boldsymbol{r}}_m} - {{\boldsymbol{r}}_n}} \right|}} = - \frac{{{\text{j}}{k_0}}}{{16{\text{π}}}}{\mathop{{\int \int}\mkern-18mu \bigcirc} {\text{e}} ^{ - {\text{j}}{k_0}({{\boldsymbol{r}}_{mp}} + {{\boldsymbol{r}}_{qn}})}}{\alpha _{pq}}({{\boldsymbol{\hat r}}_{pq}} \cdot {\boldsymbol{\hat k}}){{\text{d}}^2}{\boldsymbol{\hat k}} (14) 相比于常规八叉树分组的情况，UGFMA转移因子的计算和存储会简单很多。
1.3 动态群目标特征模方法中的加速技术
事实上，多个重复性结构在运动过程中的自耦矩阵块是不变的，因此，影响动态群目标的特征模方法效率的关键因素在于如何快速填充数量众多的互耦降阶矩阵块。多个单元间的互耦降阶矩阵块的快速填充问题可以抽象为一个场单元(设单元号为1)和一个源单元(设单元号为2)在不同排布条件下互耦降阶矩阵块 {{\boldsymbol{Z}}_{12}} 的快速计算，如图4所示，计算并存储相应排布下的互耦降阶矩阵。然而，这种策略是不现实的。假设在①、②和③过程的采样次数分别为A、B和C，则三者组合后总的采样次数为A×B×C，这将是一个异常庞大的数字。
![]()
图 4 影响互耦矩阵块的3种不同排布示意图Fig. 4 Schematic of three factors affecting the mutual coupling matrix block如前所述，组成动态群目标的单元间相距都比较大，因此它们之间的作用都可认为是远场作用。即在式(11)中， {G_p} = 0 ， {G_q} = 0 且 {{\textit{z}}_{mt}} = 0 。因此，在应用UGFMA之后，特征模方法中任意一个互耦降阶矩阵块 {{\boldsymbol{Z}}_{ij}},\;i \ne j 中的元素都可以按照以下方式计算：
\begin{split} \sum\limits_{n = 1}^N {{a_{mn}}b_n^{(t)}} =& {a_{m1}}b_1^{(t)} + {a_{m2}}b_2^{(t)} + \cdots + {a_{mN}}b_N^{(t)} \\ \approx& \mathop{{\int \int}\mkern-18mu \bigcirc} {{\boldsymbol{R}}_{mp}^{\text{S}}({\boldsymbol{\hat k}}) \cdot \sum {\varGamma _{pq}^{\text{S}}({{{\boldsymbol{\hat r}}}_{pq}} \cdot {\boldsymbol{\hat k}})} } {\boldsymbol{F}}_{qn}^{\text{S}}({\boldsymbol{\hat k}})b_n^{(t)}{{\text{d}}^2}{\boldsymbol{\hat k}} \\ =& {{d'_{mt}}} \end{split} (15) \begin{split} {c_{st}} =& \left[ {\begin{array}{*{20}{c}} {b_1^{(s)}}&{b_2^{(s)}}& \cdots &{b_N^{(s)}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{d'_{1t}}}} \\ {{{d'_{2t}}}} \\ \vdots \\ {{{d'_{Nt}}}} \end{array}} \right] \\ =& b_1^{(s)}{d_{1t}} + b_2^{(s)}{d_{2t}} + \cdots + b_N^{(s)}{d_{Nt}} \\ = &\sum\limits_{m = 1}^N {b_m^{(s)}} \mathop{{\int \int}\mkern-18mu \bigcirc} {{\boldsymbol{R}}_{mp}^{\text{S}}({\boldsymbol{\hat k}}) \cdot \varGamma _{pq}^{\text{S}}({{{\boldsymbol{\hat r}}}_{pq}} \cdot {\boldsymbol{\hat k}})} {\boldsymbol{F}}_{qn}^{\text{S}}({\boldsymbol{\hat k}})b_n^{(t)}{{\text{d}}^2}{\boldsymbol{\hat k}} \end{split} (16) 为了与常规的MLFMA进行区分，在UGFMA的聚合、转移和配置因子中增加了上标S，它们可以分别表示为：
{\boldsymbol{R}}_{mp}^{\text{S}} = \int_{{S _m}} {({\boldsymbol{\bar I}} - {\boldsymbol{\hat k\hat k}}){{\boldsymbol{f}}_m}({\boldsymbol{r}}){{\text{e}}^{ - {\text{j}}{\boldsymbol{k}} \cdot {{\boldsymbol{r}}_{mp}}}}} {\text{d}}S (17) \varGamma _{pq}^{\text{S}} = \frac{{\omega {\mu _0}{k_0}}}{{{{(4{\text{π}})}^2}}}{T_L}({{\boldsymbol{\hat r}}_{pq}} \cdot {\boldsymbol{\hat k}}) (18) {\boldsymbol{F}}_{qn}^{\text{S}} = \int_{{S _n}} {({\boldsymbol{\bar I}} - {\boldsymbol{\hat k\hat k}}){{\boldsymbol{f}}_n}({\boldsymbol{r'}}){{\text{e}}^{ - {\text{j}}{\boldsymbol{k}} \cdot {{\boldsymbol{r}}_{qn}}}}} {\text{d}}S' (19) 式中， {{\boldsymbol{f}}_m}({\boldsymbol{r}}) 和 {{\boldsymbol{f}}_n}({\boldsymbol{r'}}) 分别为场基函数和源基函数，它们分别定义在 {S _m} 和 {S _n} 上。综合式(15)～(19)，可以得出如图5所示的结论。这样就为下一步动态相关量的预存储创造了条件。
![]()
图 5 UGFMA中各变量与群目标运动之间的关系Fig. 5 The relationship between variables in UGFMA and the motion of group targets多个单元在运动过程中，所有时刻下、所有单元的自耦降阶矩阵块都是完全相同的。因此，只需要提前计算好一个单元的自耦降阶矩阵块，然后输出至数据库1中。此外，还需要提前计算单个单元在各种可能姿态下的模式聚合因子和模式配置因子，并存储至数据库2中。
2 海面背景下动态群目标的高低频混合方法——CM-PO算法
如图6所示，在本文中设置群目标为特征模快速分析方法的求解区域，简记为CM区域；而海面背景为PO法的求解区域，简记为PO区域。
![]()
图 6 群目标和海面组成的整体目标示意图Fig. 6 Schematic of the overall target composed of group targets and sea1) 首先不考虑相互耦合，用平面波分别独立照射PO区和CM区，则PO区的感应电流可利用PO法求解得到：
{{\boldsymbol{J}}_{{\text{PO}}}} = 2{\delta _0}{\boldsymbol{\hat n}} \times {{\boldsymbol{H}}^{{\text{inc}}}} (20) 式中， {\delta _0} 为遮挡因子。CM区的感应电流可根据上文中介绍的特征模快速分析方法进行求解，这里不再赘述。
2) 接着用PO的表面电流作为激励源在CM区产生散射电场 {{\boldsymbol{E}}_{{\text{CM}}}} 。散射电场公式由Huygens原理表示如下：
{\boldsymbol{E}}_{{\text{CM}}}^1({\boldsymbol{r}}) = \iint_S {\Big( { - {\text{j}}k{\eta _0}\overline{\overline {\boldsymbol{G}}} ({\boldsymbol{r}},{\boldsymbol{r'}}) \cdot {{\boldsymbol{J}}_{{\text{PO}}}}} \Big)}{\text{d}}S ' (21) 式中：上标中的1表示第一次耦合作用； \overline{\overline {\boldsymbol{G}}} 为并矢格林函数。
同时，用CM区的电流分布作为激励源在PO区产生散射磁场。散射磁场公式由Huygens原理表示如下：
{\boldsymbol{H}}_{{\text{PO}}}^{\text{1}}({\boldsymbol{r}}) = \iint_S {\nabla \times \overline{\overline {\boldsymbol{G}}} }({\boldsymbol{r}},{\boldsymbol{r'}}) \cdot {{\boldsymbol{J}}_{{\text{CM}}}}{\text{d}}S ' (22) 3) 然后将步骤2中算得的耦合场与最初入射场一起作为两个区域新的入射场，重新计算出更新后的电流 {\boldsymbol{J}}_{{\text{CM}}}^{\text{1}} 和 {\boldsymbol{J}}_{{\text{PO}}}^{\text{1}} ，如图7所示。至此，第一次相互作用结束。
![]()
图 7 CM区和PO区之间的耦合迭代过程Fig. 7 The coupling iteration process between CM area and PO area4) 重复以上过程直到电流达到稳定，即满足
\frac{{\left| {\left[ {{{\boldsymbol{I}}_{i + 1}}} \right] - \left[ {{{\boldsymbol{I}}_i}} \right]} \right|}}{{\left| {\left[ {{{\boldsymbol{I}}_i}} \right]} \right|}} < \varepsilon (23) 式中， \varepsilon 为预先设置的精度阈值。
当CM区域的目标为动态群目标时，可选择文献[13]中的CM方法，或选择采用本文第1节提出的方法。其中，前者的高低频混合方法被简记为“CM-PO”，后者的高低频混合方法则被简记为“CM-PO-Dynamic”，纯CM区的计算方法记为“CM-Dynamic”。
3 数值算例
本文算例所使用的计算机为Inter(R) Core(TM) i5-10400F CPU @ 2.9 GHz，32 GB RAM，64位Windows操作系统。以下算例采用高低频混合方法计算时，收敛标准式(23)中的精度阈值设置为 \varepsilon = 0.05 。
3.1 不同海面环境下的导弹群目标
本算例考虑在不同海况环境下由8枚导弹组成的导弹群目标，如图8所示。单枚导弹的尺寸为2.15 m×0.25 m×0.32 m。入射波频率为600 MHz，导弹单元的低阶未知量个数为2 298，导弹群总的低阶未知量个数为18 384。对单个导弹模型提取了130个模式，导弹群总的高阶未知量个数为1 040。
![]()
图 8 导弹群和海面组成的混合场景(上)以及对应的RCS(下)Fig. 8 A mixed scenario consisting of missile groups and sea (top) and corresponding RCS (bottom)为了评估不同海况对其上方飞行的群目标散射的影响，在本算例中固定群目标不动，海面(用粗糙面表示)分别考虑3种海况，即1级海况、3级海况和5级海况，同时令导弹群目标的最低点位置距离海面的高度为1 m。所有情况下海面的尺寸均为10 m×10 m，且粗糙面的未知量都是10 325。分别采用CM-PO方法、FEKO中的PO方法以及FEKO中的MoM计算以上3种情况下导弹群与海面组成的混合场景的VV极化双站雷达散射截面(radar cross section, RCS)。平面波入射角度为 \theta {\text{ = 4}}{{\text{5}}{\text{°}} }、 \varphi = {0{\text{°}} } ，接收角度为 \theta {\text{ = 4}}{{\text{5}}{\text{°}} }、\varphi = {0{\text{°}} }\sim {360{\text{°}} } 。在图8中分别给出了1级、3级和5级海况情况下不同方法计算的RCS的对比。
以FEKO(MoM)的结果作为参考值来评估其他两种方法的精度，从图8可以看出，不同海况下，CM-PO方法和FEKO(MoM)计算的结果具有较好的一致性，而FEKO(PO)方法和FEKO(MoM)方法计算的结果拟合较差。在计算效率方面，表1统计了3种方法的计算资源消耗情况。
表 1 不同方法计算资源消耗对比(3种海况的平均值)Tab. 1 Comparison of computational resource consumption between different methods (average of 3 cases)方法 峰值内存 耗时 FEKO(PO) 6.62 MB 3.50 s CM-PO 490.54 MB 8.39 min FEKO(MoM) 6.16 GB 29.61 min 从表1可以发现：尽管FEKO(MoM)计算结果最准确，但是需要消耗非常多的计算资源；FEKO(PO)尽管计算资源消耗最少，但计算的精度又太差；CM-PO方法平衡了计算效率与计算精度，实现了在不损失较大精度的前提下准确分析混合场景目标的效果。
3.2 海面环境下的动态角反射器阵列
本算例考虑在海面环境下由5个角反射器组成的动态阵列。在此算例中，保持海面不动，而角反射器阵列运动。单个角反射器的尺寸为3.1 m×3.1 m×3.1 m。入射波频率为300 MHz，角反射器单元的低阶未知量为25 700，角反射器阵列的总低阶未知量个数为128 500。对单个角反射器模型提取250个特征模式，角反射器阵列的总的高阶未知量个数为1 250。在动态角反射器阵列中，不同时刻下两两角反射器单元质心之间的距离范围为3.5～7 m。每个角反射器在实际海面中，可以认为其受海风及海面起伏的影响会以轴线为基准产生 \pm {45{\text{°}}} 的倾斜，因此数据库建立的方式为：俯仰角 \varphi {\text{ = }} - {45{\text{°}}}\sim {\text{4}}{{\text{5}}{\text{°}}} 按0.3°等间隔建库，共315个采样，时长共为1.25 h，所占存储内存为2.38 GB。
分别采用不同的方法计算图9中3个时刻下角反射器阵列场景VV极化双站RCS、角反射器阵列与海面(用粗糙面表示，参与计算海域用黑色虚线框标记出，是尺寸为17 m×12 m风速为4 m/s的PM谱海面，未知量个数为95 765)组成的混合场景VV极化双站RCS，以及角反射器阵列与海面组成的混合场景VH极化双站RCS，结果分别展示在图10、图11和图12中。
![]()
图 9 3个时刻下动态角反射器阵列和海面组成的混合场景Fig. 9 A mixed scene consisting of a dynamic corner reflector array and the sea at 3 time points![]()
图 10 3个时刻下动态角反射器阵列场景的VV极化RCSFig. 10 VV polarization RCS of dynamic corner reflector array scene at 3 time points![]()
图 11 3个时刻下动态角反射器阵列和海面组成的混合场景的VV极化RCSFig. 11 VV polarization RCS of a mixed scene composed of a dynamic corner reflector array and the sea at 3 time points![]()
图 12 3个时刻下动态角反射器阵列和海面组成的混合场景的VH极化RCSFig. 12 VH polarization RCS of a mixed scene composed of a dynamic corner reflector array and the sea at 3 time points平面波入射角度为 \theta {\text{ = 4}}{{\text{5}}{\text{°}} }、\varphi = {0{\text{°}} } ，接收角度为 \theta {{ = 0 \sim 9}}{{\text{0}}{\text{°}} }、\varphi = {180{\text{°}} } 。为了测试算法的性能，每个时刻下各角反射器单元的姿态和位置都是以随机的方式生成的。
表2给出了3个时刻下各方法计算VV极化RCS的结果与FEKO(MoM)给出的VV极化RCS结果之间的相对误差。FEKO的MoM方法是加MLFMA的MoM方法，下同。
表3给出了3个时刻下各方法计算VH极化RCS结果与FEKO (MoM)给出的VH极化RCS结果之间的相对误差。
表 2 VV极化下不同方法与FEKO (MoM) RCS结果间的相对误差Tab. 2 Relative error between different methods under VV polarization and RCS results provided by FEKO (MoM)dB 方法 时刻1 时刻2 时刻3 FEKO (PO)(角反阵+海面) 4.04 5.59 4.72 CM-PO-Dynamic(角反阵+海面有耦合) 2.63 2.03 2.58 CM-PO-Dynamic(角反阵+海面无耦合) 3.89 4.86 5.46 CM-Dynamic(角反阵) 1.54 0.89 1.02 表 3 VH极化下不同方法与FEKO (MoM) RCS结果间的相对误差Tab. 3 Relative error between different methods under VH polarization and RCS results provided by FEKO (MoM)dB 方法 时刻1 时刻2 时刻3 FEKO (RL-GO)(角反阵+海面) 11.62 11.45 11.45 CM-PO-Dynamic(角反阵+海面有耦合) 7.42 7.70 6.51 表4中给出了3个时刻下各个方法的平均计算效率。
表 4 不同方法3个时刻下的平均计算效率Tab. 4 Average computational efficiency at 3 time points of different methods方法 计算总时间 峰值内存 FEKO (MoM)(角反阵) 31 min 4.0 GB CM-Dynamic(角反阵) 17 s 1.1 GB FEKO (MoM)(角反阵+海面有耦合) 2.5 h 6.0 GB FEKO (PO)(角反阵+海面) 11 s 49.3 MB FEKO (RL-GO)(角反阵+海面) 33 s 74.8 MB CM-PO-Dynamic(角反阵+海面有耦合) 27.6 min 1.9 GB CM-PO-Dynamic(角反阵+海面无耦合) 14.3 min 1.5 GB 综合图10~12以及表2~3可以发现，针对由5个角反射器组成的动态阵列采用MoM分析，本文通过预建库、读库计算的CM-Dynamic方法实现了动态群目标的多时刻快速分析，虽然RCS结果上牺牲了一定的精度但仍然在可接受范围内。针对角反阵和海面的混合场景VV极化RCS的计算，如果整个场景都使用PO方法进行计算，RCS结果的误差会比较大，采用高低频混合的CM-PO-Dynamic方法在不考虑海面和角反阵耦合的情况下RCS结果的误差仍然很大，而考虑耦合作用的CM-PO-Dynamic方法能在PO方法的基础上明显提高计算的精度，且在分析这些时刻下都保持了较小且稳定的计算误差。此外本文的方法在计算交叉极化时精度相比于考虑多次弹跳的纯高频方法有所提高。通过表4可见，在计算效率方面，尽管FEKO (MoM)计算的结果最准确，但是它单次仿真时间在2.5 h左右；本文提出的CM-PO-Dynamic方法尽管损失了少量精度，但其仿真时间缩短到了27.6 min，更加合理地实现了计算效率和计算精度的平衡。表5给出了采用所提出的方法计算角反射器阵列在不同海情和入射角下满足式(23)所需要的收敛步数，可以看出随着海浪变高或者入射角变小，所需要的收敛步数变多。
表 5 不同海情和入射角下式(23)所需要的收敛步数Tab. 5 The required convergence steps for equation (23) under different sea conditions and incident angles海情 \theta = {0{\text{°}}} \theta = {30{\text{°}}} \theta = {60{\text{°}}} \theta = {90{\text{°}}} 3 m/s PM谱海面 2 1 1 1 9 m/s PM谱海面 3 3 2 1 15 m/s PM谱海面 4 3 2 2 4 结 论
本文针对海面环境下动态群目标的电磁散射问题，提出了一种新的高低频混合方法。在动态群目标的处理上采用特征模方法以及按照单元分组的FMA来分析，采用PO法来分析海面背景，能够实现电大粗糙海面的快速分析。将以上两种算法进行有机结合，以多次迭代的方式考虑海面和动态群目标之间的耦合，使整体的分析效率和精度实现了平衡。由于结合各算法的特点合理地分配了各算法的计算任务，该混合方法可以在高效计算的同时可实现对整个场景更加全面、合理的动态仿真。
-
![]()
图 1 动态群目标的全局八叉树分组示意图
Fig. 1 Schematic of global octree grouping for dynamic group targets
![]()
图 3 UGFMA中不同单元间基函数的作用路线示意图
Fig. 3 Schematic of the action path of basis functions between different units in UGFMA
![]()
图 4 影响互耦矩阵块的3种不同排布示意图
Fig. 4 Schematic of three factors affecting the mutual coupling matrix block
![]()
图 5 UGFMA中各变量与群目标运动之间的关系
Fig. 5 The relationship between variables in UGFMA and the motion of group targets
![]()
图 6 群目标和海面组成的整体目标示意图
Fig. 6 Schematic of the overall target composed of group targets and sea
![]()
图 7 CM区和PO区之间的耦合迭代过程
Fig. 7 The coupling iteration process between CM area and PO area
![]()
图 8 导弹群和海面组成的混合场景(上)以及对应的RCS(下)
Fig. 8 A mixed scenario consisting of missile groups and sea (top) and corresponding RCS (bottom)
![]()
图 9 3个时刻下动态角反射器阵列和海面组成的混合场景
Fig. 9 A mixed scene consisting of a dynamic corner reflector array and the sea at 3 time points
![]()
图 10 3个时刻下动态角反射器阵列场景的VV极化RCS
Fig. 10 VV polarization RCS of dynamic corner reflector array scene at 3 time points
![]()
图 11 3个时刻下动态角反射器阵列和海面组成的混合场景的VV极化RCS
Fig. 11 VV polarization RCS of a mixed scene composed of a dynamic corner reflector array and the sea at 3 time points
![]()
图 12 3个时刻下动态角反射器阵列和海面组成的混合场景的VH极化RCS
Fig. 12 VH polarization RCS of a mixed scene composed of a dynamic corner reflector array and the sea at 3 time points
表 1 不同方法计算资源消耗对比(3种海况的平均值)
Tab. 1 Comparison of computational resource consumption between different methods (average of 3 cases)
方法 峰值内存 耗时 FEKO(PO) 6.62 MB 3.50 s CM-PO 490.54 MB 8.39 min FEKO(MoM) 6.16 GB 29.61 min 
下载: 导出CSV
表 2 VV极化下不同方法与FEKO (MoM) RCS结果间的相对误差
Tab. 2 Relative error between different methods under VV polarization and RCS results provided by FEKO (MoM)
dB 方法 时刻1 时刻2 时刻3 FEKO (PO)(角反阵+海面) 4.04 5.59 4.72 CM-PO-Dynamic(角反阵+海面有耦合) 2.63 2.03 2.58 CM-PO-Dynamic(角反阵+海面无耦合) 3.89 4.86 5.46 CM-Dynamic(角反阵) 1.54 0.89 1.02
下载: 导出CSV
表 3 VH极化下不同方法与FEKO (MoM) RCS结果间的相对误差
Tab. 3 Relative error between different methods under VH polarization and RCS results provided by FEKO (MoM)
dB 方法 时刻1 时刻2 时刻3 FEKO (RL-GO)(角反阵+海面) 11.62 11.45 11.45 CM-PO-Dynamic(角反阵+海面有耦合) 7.42 7.70 6.51
下载: 导出CSV
表 4 不同方法3个时刻下的平均计算效率
Tab. 4 Average computational efficiency at 3 time points of different methods
方法 计算总时间 峰值内存 FEKO (MoM)(角反阵) 31 min 4.0 GB CM-Dynamic(角反阵) 17 s 1.1 GB FEKO (MoM)(角反阵+海面有耦合) 2.5 h 6.0 GB FEKO (PO)(角反阵+海面) 11 s 49.3 MB FEKO (RL-GO)(角反阵+海面) 33 s 74.8 MB CM-PO-Dynamic(角反阵+海面有耦合) 27.6 min 1.9 GB CM-PO-Dynamic(角反阵+海面无耦合) 14.3 min 1.5 GB
下载: 导出CSV
表 5 不同海情和入射角下式(23)所需要的收敛步数
Tab. 5 The required convergence steps for equation (23) under different sea conditions and incident angles
海情 \theta = {0{\text{°}}} \theta = {30{\text{°}}} \theta = {60{\text{°}}} \theta = {90{\text{°}}} 3 m/s PM谱海面 2 1 1 1 9 m/s PM谱海面 3 3 2 1 15 m/s PM谱海面 4 3 2 2
下载: 导出CSV
-
[1] COIFMAN R,ROKHLIN V,WANDZURA S. The fast multipole method for the wave equation:a pedestrian prescription[J]. IEEE antennas and propagation magazine,1993,35(3):7-12. doi: 10.1109/74.250128
[2] SONG J,LU C C,CHEW W C. Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects[J]. IEEE transactions on antennas and propagation,1997,45(10):1488-1493. doi: 10.1109/8.633855
[3] EWE W B,LI E P,CHU H S,et al. AIM analysis of electromagnetic scattering by arbitrarily shaped magnetodielectric object[J]. IEEE transactions on antennas and propagation,2007,55(7):2073-2079. doi: 10.1109/TAP.2007.900263
[4] ZWAMBORN P,VAN DEN BERG P M. The three dimensional weak form of the conjugate gradient FFT method for solving scattering problems[J]. IEEE transactions on microwave theory and techniques,1992,40(9):1757-1766. doi: 10.1109/22.156602
[5] PENG Z,WANG X,LEE J F. Integral equation based domain decomposition method for solving electromagnetic wave scattering from non-penetrable objects[J]. IEEE transactions on antennas and propagation,2011,59(9):3328-3338. doi: 10.1109/TAP.2011.2161542
[6] 盛新庆. 计算电磁学要论[M]. 2版. 安徽:中国科学技术大学出版社,2008. [7] LAVIADA J,LAS-HERAS F,PINO M R,et al. Solution of electrically large problems with multilevel characteristic basis functions[J]. IEEE transactions on antennas and propagation,2009,57(10):3189-3198. doi: 10.1109/TAP.2009.2028603
[8] LUCENTE E,MONORCHIO A,MITTRA R. An iteration-free MoM approach based on excitation independent characteristic basis functions for solving large multiscale electromagnetic scattering problems[J]. IEEE transactions on antennas and propagation,2008,56(4):999-1007. doi: 10.1109/TAP.2008.919166
[9] MATEKOVITS L,VECCHI G,DASSANO G,et al. Synthetic function analysis of large printed structures:The solution space sampling approach[C]//IEEE Antennas and Propagation Society International Symposium. 2001 Digest. Held in conjunction with:USNC/URSI National Radio Science Meeting (Cat. No. 01CH37229). IEEE,2001,2:568-571.
[10] LU W B,CUI T J,QIAN Z G,et al. Accurate analysis of large-scale periodic structures using an efficient sub-entire-domain basis function method[J]. IEEE transactions on antennas and propagation,2004,52(11):3078-3085. doi: 10.1109/TAP.2004.835143
[11] HARRINGTON R,MAUTZ J. Theory of characteristic modes for conducting bodies[J]. IEEE transactions on antennas and propagation,1971,19(5):622-628. doi: 10.1109/TAP.1971.1139999
[12] GUAN L,HE Z,DING D,et al. Efficient characteristic mode analysis for radiation problems of antenna arrays[J]. IEEE transactions on antennas and propagation,2018,67(1):199-206.
[13] WU C,GUAN L,GU P,et al. Application of parallel CM-MLFMA method to the analysis of array structures[J]. IEEE transactions on antennas and propagation,2021,69(9):6116-6121. doi: 10.1109/TAP.2021.3069427
[14] ZHANG H L,SHA Y X,GUO X Y,et al. Efficient analysis of scattering by multiple moving objects using a tailored MLFMA[J]. IEEE Transactions on antennas and propagation,2019,67(3):2023-2027. doi: 10.1109/TAP.2019.2891226
[15] ZHANG H L,SHA Y X,HE X Y,et al. Efficient algorithm for scattering by a large cluster of moving objects[J]. IEEE access,2019,7:124948-124955. doi: 10.1109/ACCESS.2019.2937598
[16] LI M,HU Y,CHEN R,et al. Electromagnetic modeling of moving mixed conductive and dielectric BoRs with an effective domain decomposition method[J]. IEEE transactions on antennas and propagation,2020,68(12):7978-7985. doi: 10.1109/TAP.2020.2998914
[17] SPENCER R,HYDE G. Studies of the focal region of a spherical reflector:geometric optics[J]. IEEE transactions on antennas and propagation,1968,16(3):317-324. doi: 10.1109/TAP.1968.1139187
[18] MICHAELI A. Equivalent edge currents for arbitrary aspects of observation[J]. IEEE transactions on antennas and propagation,1984,32(3):252-258. doi: 10.1109/TAP.1984.1143303
[19] UFIMTSEV P Y. Elementary edge waves and the physical theory of diffraction[J]. Electromagnetics,1991,11(2):125-160. doi: 10.1080/02726349108908270
[20] KELLER J B. A geometrical theory of diffraction in calculus of variations and its applications[M]. New York:Mc Graw-Hill,1958.
[21] PATHAK P,BURNSIDE W,MARHEFKA R. A uniform GTD analysis of the diffraction of electromagnetic waves by a smooth convex surface[J]. IEEE transactions on antennas and propagation,1980,28(5):631-642. doi: 10.1109/TAP.1980.1142396
[22] FAN T Q,GUO L X,LIU W. A novel OpenGL-based MoM/SBR hybrid method for radiation pattern analysis of an antenna above an electrically large complicated platform[J]. IEEE transactions on antennas and propagation,2015,64(1):201-209.
[23] DENG J Y,GUO L X. An efficient octree-based Mom-PO method for analysis of antennas on large platform[J]. IEEE antennas and wireless propagation letters,2015,14:819-822. doi: 10.1109/LAWP.2015.2399372
[24] HARRINGTON R F. Field computation by moment methods [M]. New York:Macmillan,1968.
[25] DAI Q I,WU J,GAN H,et al. Large-scale characteristic mode analysis with fast multipole algorithms[J]. IEEE Transactions on antennas & propagation,2016,64(7):2608-2616.
[26] 樊振宏. 电磁散射分析中的快速方法[D]. 南京:南京理工大学,2007. FAN Z H. Fast methods in electromagnetic scattering analysis[D]. Nanjing:Nanjing University of Science and Technology,2007. (in Chinese)
计量
- 文章访问数: 127
- HTML全文浏览量: 20
- PDF下载量: 59
