CN102129523B - 基于mda和mlssm的分析复杂目标电磁散射的方法 - Google Patents

Abstract

本发明公开了一种基于MDA和MLSSM的分析复杂目标电磁散射的方法。该方法首先对待求物体建立八叉树结构，通过MDA方法来获得离散电场积分方程形成的稠密阻抗矩阵中对应于远场作用部分的稀疏低秩表示，然后将阻抗矩阵中低秩表示的部分用改进的MLSSM表示出来。最后基于改进的MLSSM表示，构造一种适用于迭代求解的快速计算矩阵矢量乘的算法。本发明能够大大减少用电场积分方程方法分析复杂目标电磁散射的计算时间和内存需求。

Description

基于MDA和MLSSM的分析复杂目标电磁散射的方法

技术领域

本发明涉及一种积分方程方法求解电磁问题里加速矩阵矢量乘计算的快速技术，特别是一种基于MDA和MLSSM的分析复杂目标电磁散射的方法。

背景技术

近年来，分析任意形状三维目标的电磁散射特性，已经受到越来越多国内外研究学者的关注。其应用主要表现在：雷达系统的设计与目标识别、军用武器的隐身与反隐身、复杂环境中的电磁兼容问题等等。由于雷达工作在微波频段，常见军用目标，如战车、导弹、飞机等等不但具有复杂的几何结构，而且具有超大的电尺寸，增加了散射特性分析的复杂度。因此，如何高效计算这类大型复杂目标的电磁散射特性，如目标的雷达截面积RCS 参数的提取，一直是从事雷达总体设计以及隐身与反隐身研究的学者、工程师们所共同关心的问题。

对于目标电磁散射特性的分析，传统的解析方法仅局限于解决几何形状为球体或无限长圆柱体等简单目标。随着计算机技术与计算电磁学理论的迅猛发展，出现了许多可用于分析任意形状复杂目标电磁散射特性的数值方法。这些方法大致可划分为基于高频近似的高频RCS预估方法和基于微分方法或积分方程的严格数值方法两大类。

基于高频近似的高频方法，如几何光学法(GO)、物理光学法(PO)、物理绕射理论(PTD)、弹跳射线方法(SBR)等等，由于具有计算速度快，消耗机器内存少的优点被广泛应用于各类复杂目标电磁散射特性的分析。然而，由于高频近似条件的引入，使得高频方法的计算结果精度低，这是制约高频方法发展及其应用的一个重要的因素。用于散射特性分析的微分方程法，包括时域有限差分法(FDTD)和有限元法(FEM)，虽然得到稀疏矩阵，但是难以精确拟合目标的曲面边界。并且需要在传播空间内进行网格划分，存在网格数值色散。对于开域问题，微分方程的求解必须施加吸收边界条件，常用于分析目标的电磁散射特性。积分方程法又可分为体积分方程法和表面积分方程法两种。体积分方程法适合于非均匀介质材料的电磁散射特性分析，而表面积分方程法适合于理想导体或均匀介质材料的电磁散射特性分析。

尽管基于积分方程的矩量法具有严格的理论模型，但是其生成的矩阵为满阵。假定N为未知量个数，则存储该稠密矩阵将要耗费O(N²)的内存量。同时，如果利用直接法来求解矩量法的阻抗矩阵方程，其计算复杂度为O(N²)。对于电大尺寸的目标, 问题变得越发严重。正是由于受到计算机存储量的限制，长期以来矩量法仅仅局限于分析低频区到谐振区范围内目标的电磁散射特性。近几十年来，由于计算机技术与计算电磁学理论的迅猛发展，矩量法(MoM)作为一种严格的数值方法，重新引起了人们极大的关注。特别是各种高效的快速积分方法的提出，大幅度地降低了矩量法分析所需的计算量以及内存量。其中自由空间里多层快速多极子(MLFMA)是其中最有代表性的快速算法。MLFMA是依赖于Green函数，因此应用到多层媒质的情形，MLFMA就需要做修改。而基于矩阵低秩分解的代数类快速算法是不依赖于Green函数，因此很容易直接应用到多层媒质的情形。低秩类方法主要有ACA，UV，MDA。其中MDA主要利用等效源的思想，最初是应用在二维的情况(E.Michielssen and A. Boag, “A multilevel matrix decomposition algorithm for analyzing scattering from largestructures”, IEEE Trans. Antennas Propag., vol. 44,no. 8, pp. 1086-1093, August 1996.)，后来由J.M. Rius推广到三维的情况(J.M. Rius, J. Parr′on, E. beda, J.R.Mosig, “Multilevel Matrix Decomposition Algorithm forAnalysis of Electrically Large Electromagnetic Problemsin 3-D”, Microwave and Optical Technology Letters,Vol. 22, No.3, pp. 177-182, 5th August 1999.)。相对于其它方法而言，MDA在计算平面结构的目标在计算时间和内存需求上有优势。

由于积分方程方法得到的阻抗矩阵是稠密的，F. X. Canning提出Multilevel Simple Sparse Method (MLSSM) 矩阵稀疏化表示的方法。最初是F. X. Canning 在利用直接方法来求解电场积分方程时存储阻抗矩阵 (1. F. X. Canning and K. Rogovin, "Auniversal matrix solver for integral-equation based problems," IEEE Antennas and Propagation Magazine, vol.45, February, 2003; 2. A. Zhu, R. J. Adams, and F.X. Canning, "Modified simply sparse method for electromagnetic scattering by PEC," IEEE Antennas and Propagation Society International Symposium, Washington, DC, 2005, pp. 427-430.) 。随着目标电尺寸增大，待求的未知量增多，迭代求解过程中矩阵矢量乘的计算复杂度是求解的关键技术。基于MLSSM稀疏表示的迭代解法还没有报道。

发明内容

本发明所解决的技术问题在于提供一种基于MDA和MLSSM的分析复杂目标电磁散射的方法。

实现本发明目的的技术解决方案为：一种基于MDA和MLSSM的分析复杂目标电磁散射的方法，包括以下步骤：

步骤1、对入射平面波参数进行设置，根据实际物理背景，对入射平面波的频率和入射角度进行设置，入射平面波的频率为微波频段，入射角度在球坐标下为(θ,φ)，其中θ∈[0°,180°], φ∈[0°,360°]；

步骤2、利用Ansys软件建立待求目标的几何模型，之后利用平面三角形进行离散剖分；

步骤3、采用八叉树结构对剖分后的目标模型进行分组，用一个立方体将目标体包围住，该立方体定义为第零层的第一个且是最后一个组结点，把该立方体等分为八个子立方体结点形成第一层组结点，然后再对每个子立方体进行与上一步相同的细分，并以此类推直到最底层立方体的电尺寸达到所需合适的大小为止；

步骤4、根据步骤3得到的八叉树结构，首先进行Morton编码，然后将每个立方体相邻的立方体组设定为近场区组，之后设定远场区，所述远场区为包含该立方体的父层立方体区域的近场区域中除掉本层区域的近场区组；

步骤5、根据立方体电尺寸的大小，在其表面设置等效源，利用等效源对阻抗矩阵的远作用子块进行低秩表示，并直接计算阻抗矩阵的近作用子块；

步骤6、利用QR和SVD对步骤5得到的低秩子块表示进一步压缩处理，得到阻抗矩阵改进的MLSSM表示；

步骤7、根据步骤6得到改进的MLSSM表示，构造迭代求解的快速矩阵矢量乘算法；

步骤8、求解步骤7中的阻抗矩阵方程，得到模型表面电流系数，并且根据电流系数得到目标表面的电流分布，从而得到目标的电磁散射特性参数RCS。

本发明与现有技术相比，其显著优点为：1）利用本方法最终形成的MLSSM表示降低了计算机内存的需求，从而在有限的计算机资源下处理更大点尺寸问题；2）在迭代法求解的时候，能实现快速的矩阵矢量乘，从而减少迭代求解时间。因此，本发明可以在有限计算资源下更快处理更大的问题。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明基于MDA和MLSSM的快速矩阵矢量乘技术的流程图。

图2是sphere-cone在Ansys下的剖分示意图。

图3是三维的八叉树结构的建立过程示意图，从左至右分别是第零层，第一层，第二层。

图4 阻抗矩阵的示意图。

图5 根据等效原理建立远场矩阵子块的低秩表示的过程示意图。

图6 是MLSSM表示的示意图，其中(a)是第四层改进的MLSSM表示中的U₄，，V₄示意图； (b) 第三层改进的MLSSM表示中的U₄，

，V₄，U₃，

，V₃示意图。

图7是sphere-cone模型随电尺寸变化的矩阵矢量乘时间曲线示意图。

图8是sphere-cone模型随电尺寸变化的内存需求曲线示意图。

图9是入射波频率为8 GHz ,入射角是θ=0°,φ=0°时, sphere-cone模型的双站RCS曲线图。

具体实施方式

结合图1，本发明的一种基于MDA和MLSSM的分析复杂目标电磁散射的方法，包括以下步骤：

步骤1、对入射平面波参数进行设置，根据实际物理背景，对入射平面波的频率和入射角度进行设置；入射平面波的频率为微波频段，入射角度在球坐标下为(θ,φ)，其中θ∈[0°,180°], φ∈[0°,360°]。

步骤2、利用Ansys软件建立待求目标的几何模型，之后利用平面三角形进行离散剖分；利用平面三角形进行离散剖分时，采用基于RWG基函数的三角形网格对复杂目标表面进行剖分，每平方电波长内剖分的三角形数目大于120。

步骤3、采用八叉树结构对剖分后的目标模型进行分组，用一个立方体将目标体包围住，该立方体定义为第零层的第一个且是最后一个组结点，把该立方体等分为八个子立方体结点形成第一层组结点，然后再对每个子立方体进行与上一步相同的细分，并以此类推直到最底层立方体的电尺寸达到所需合适的大小为止；最底层立方体的电尺寸为0.2～0.4个入射波波长。

步骤4、根据步骤3得到的八叉树结构，首先进行Morton编码，然后将每个立方体相邻的立方体组设定为近场区组，之后设定远场区，所述远场区为包含该立方体的父层立方体区域的近场区域中除掉本层区域的近场区组；

步骤5、根据立方体电尺寸的大小，在其表面设置等效源，利用等效源对阻抗矩阵的远作用子块进行低秩表示，并直接计算阻抗矩阵的近作用子块；所述等效源的个数Q按照以下标准确定：

$Q\≥10\sqrt{k\frac{S_s{S}_o}{d}}$

其中S_s，S_o，d分别是远作用的时候场所在的立方体和源所在的立方体的电尺寸，以及它们之间的距离，k是入射波的波数，每一个等效源就是一个RWG基函数，而且这些等效源均匀分布在立方体的表面。利用等效源对阻抗矩阵的远作用子块进行低秩表示具体为：

$\left[Z_{\mathrm{mn}}\right]=\left[U\right]{\left[V\right]}^T$

其中[Z_mn]∈C^m×n是阻抗矩阵中对应于远场作用的子块，[U]∈C^m×r，[V]∈C^n×r，m、n、r为正整数，其中r小于min(m，n)。

步骤6、利用QR和SVD对步骤5得到的低秩子块表示进一步压缩处理，得到阻抗矩阵改进的MLSSM表示；阻抗矩阵改进的MLSSM表示为：

$Z_l={\hat{Z}}_l+U_l{Z}_{(l-1)}{V}_l^H$

其中Z_L=Z是阻抗矩阵，

表示压缩后第二层表示，V_l=(U_l)^T是块对角的酉矩阵， l表示的是层数，0≤l≤L。

步骤7、根据步骤6得到改进的MLSSM表示，构造迭代求解的快速矩阵矢量乘算法；构造迭代求解的快速矩阵矢量乘算法按照以下函数执行：

其中x,y分别为输入向量和输出向量，y₁，y₂，y₃为临时变量，U_l和

是MLSSM表示中的矩阵，l表示的是层数。

步骤8、求解步骤7中的阻抗矩阵方程，得到模型表面电流系数，并且根据电流系数得到目标表面的电流分布，从而得到目标的电磁散射特性参数RCS。

下面以图2所示sphere-cone结构为例，对本发明的具体步骤作进一步详细描述。

根据本发明分析图2所示sphere-cone 结构的基于Matrix Decomposition Algorithm和Multilevel Simple Sparse Method的分析复杂目标电磁散射的快速方法，其具体操作步骤如下：

第一步，入射波参数的设置。根据实际物理背景，设置入射波为平面波。这里为了验证我们快速方法的有效性, 对于相同的sphere-cone结构, 入射波的频率设置为2 GHz, 4 GHz, 8 GHz, 得到电尺寸逐渐增大的结构。入射角度是θ=0°,φ=0°

第二步，目标的几何建模。根据目标模型的几何参数，用Ansys建立目标的几何模型。这里建立的模型的时候按十分之一波长来剖分sphere-cone。

第三步，采用八叉树结构对剖分后的sphere-cone结构进行分组，先用一个立方体将第一层sphere-cone包围住，该立方体定义为第零层的第一个且是最后一个组结点，把该立方体等分为八个子立方体结点形成第一层组结点，然后再对每个子立方体进行与上一步相同的细分，并以此类推直到最底层立方体的尺寸达到所需要求为止，一般为0.2～0.4个电波长。

第四步，对第三步建立的树形结构采用Morton编码。首先查找每一层非空的小正方体，其次通过树形结构得到每一个非空小正方体近场作用和远场作用小正方体的信息。这里我们采用Morton编码可以快速进行索引，从而提高计算效率。建立树形结构的目的是为了确定阻抗矩阵中需要满秩填充(近场作用)和低秩填充(远场作用)的子块。

第五步，利用矩阵分解算法（Matrix Decomposition Algorithm）对阻抗矩阵的远场作用子块进行填充。根据等效原理(如图5所示)，可以对远场作用子块表示为以下形式：

$\left[Z_{\mathrm{mn}}\right]=\left[Z_{\mathrm{mp}}\right]{\left[Z_{\mathrm{qp}}\right]}^{-1}\left[Z_{\mathrm{qn}}\right]$

其中Z_mn是待填充的子块，维数是m×n；Z_mp，Z_qp，Z_qn维数分别是m×p，p×q，q×n；p,q是等效源的数目，m,n>>p,q。也就是将待填充的子块Z_mn用Z_mp，Z_qp，Z_qn的乘积来低秩表示。

第六步，利用QR分解（QR Decomposition）和奇异值分解（Singular Value Decomposition —— SVD分解）对上面得到的Z_mn低秩表示进一步压缩。具体过程如下：

做SVD分解有

$\left[Z_{\mathrm{pq}}\right]=\left[U^{\′}\right]\left[w^{\′}\right]{\left[V^{\′}\right]}^T$

$\left[Z_{\mathrm{mn}}\right]=\left[Z_{\mathrm{mp}}\right]\left[V^{\′}\right]\left[\frac{1}{w^{\′}}\right]{\left[U^{\′}\right]}^T\left[Z_{\mathrm{qn}}\right]$

做QR分解有

$\left[Z_{\mathrm{mp}}\right]\left[V^{\′}\right]=Q_1{R}_1$  和  ${\left[Z_{\mathrm{qn}}\right]}^T{U}^{\′}=Q_2{R}_2$

从而有

$\left[R_1\right]\left[\frac{1}{w^{\′}}\right]{\left[{R_2}^{}\right]}^T=\left[U^{\′\′}\right]\left[w^{\′\′}\right]{\left[V^{\′\′}\right]}^T$

因此可以得到

$\left[Z_{\mathrm{mn}}\right]=\left[Q_1\right]\left[U^{\′\′}\right]\left[w^{\′\′}\right]{\left[V^{\′\′}\right]}^T{\left[Q_2\right]}^T={[U}_1\left]\right[w^{\′\′}]{\left[V_1\right]}^T=[U]{\left[V\right]}^T$

这样就得到压缩后的低秩表示矩阵。

第七步，利用上面得到远作用子块的低秩表示，得到阻抗矩阵Z改进的Multilevel Simple Sparse Method表示如下：

$Z_l={\hat{Z}}_l+U_l{Z}_{(l-1)}{V}_l^H$

其中Z_L=Z，

，V_l=(U_l)^T是块对角的酉矩阵。

构造改进的Multilevel Simple Sparse Method具体过程如下：

(1) 先对L层每一个非空场立方体L_i，根据八叉树结构找到l(3≤l≤L)层的远场作用源立方体l_j。而每个源立方体对应阻抗矩阵的一个低秩表示子块U(L_i,l_j)ω(L_i,l_j)V(L_i,l_j)，然后将U(L_i,l_j)ω(L_i,l_j)排成一行形成矩阵Z(L_i,l)。

(2) 将(1)得到的矩阵Z(L_i,l)排成一行形成矩阵Z(L_i)。对Z(L_i)做SVD操作得到Z(L_i)=U(L_i)S(L_i)V^H(L_i)。因此，对L层每一个非空场立方体L_i对应了一个矩阵U(L_i)。将这些子块按对角线排列就形成了块对角的酉矩阵U_L。令V_L=(U_L)^T就得到了V_L。

(3) 令Z_L-1(L_i,l_j)表示Z_L-1中与场立方体L_i和源立方体l_j相联系的子块，其中 l(3≤l≤L)。按照以下方式得到 $Z_{L-1}(L_i,l_j)=\hat{S}\left(L_i\right){\hat{V}}^H\left(L_i\right)V(L_i,l_j){\hat{U}}^C\left(L_i\right)$ 。

(4) 对于l(3≤l≤L-1), 按照上述(1)-(3)同样的步骤就可以得到

, U_l , V_l。在表示完成之后，只需要存储U_l，

 (2≤l≤L)。图5显示第四层和第三层改进的MLSSM表示中的U₄，

，V₄，U₃，

，V₃示意图。根据上面得到改进的MLSSM表示，构造以下快速矩阵矢量乘算法。

其中x,y分别为输入向量和输出向量。如图7-8所示，相比MDA-SVD，现有的改进的MLSSSM表示在内存需求和矩阵矢量乘的时间有很大的优势。

第八步，利用广义最小余量法求解线性方程组的解得到电流系数。根据电流系数得到目标的散射参数RCS。

Claims (7)

1.一种基于MDA和MLSSM的分析复杂目标电磁散射的方法，其特征在于，包括以下步骤：

步骤1、对入射平面波参数进行设置，根据实际物理背景，对入射平面波的频率和入射角度进行设置，入射平面波的频率为微波频段，入射角度在球坐标下为(θ,φ)，其中θ∈[0°,180°], φ∈[0°,360°]；

步骤2、利用Ansys软件建立待求目标的几何模型，之后利用平面三角形进行离散剖分；

步骤3、采用八叉树结构对剖分后的目标模型进行分组，用一个立方体将目标体包围住，该立方体定义为第零层的第一个且是最后一个组结点，把该立方体等分为八个子立方体结点形成第一层组结点，然后再对每个子立方体进行与上一步相同的细分，并以此类推直到最底层立方体的电尺寸达到所需合适的大小为止；

步骤4、根据步骤3得到的八叉树结构，首先进行Morton编码，然后将每个立方体相邻的立方体组设定为近场区组，之后设定远场区，所述远场区为包含该立方体的父层立方体区域的近场区域中除掉本层区域的近场区组；

步骤5、根据立方体电尺寸的大小，在其表面设置等效源，利用等效源对阻抗矩阵的远作用子块进行低秩表示，并直接计算阻抗矩阵的近作用子块；

步骤6、利用QR和SVD对步骤5得到的低秩子块表示进一步压缩处理，得到阻抗矩阵改进的MLSSM表示；

步骤7、根据步骤6得到改进的MLSSM表示，构造迭代求解的快速矩阵矢量乘算法；

步骤8、求解步骤7中的阻抗矩阵方程，得到模型表面电流系数，并且根据电流系数得到目标表面的电流分布，从而得到目标的电磁散射特性参数RCS。

2.根据权利要求1所述的基于MDA和MLSSM的分析复杂目标电磁散射的方法，其特征在于，步骤2中利用平面三角形进行离散剖分时，采用基于RWG基函数的三角形网格对复杂目标表面进行剖分，每平方电波长内剖分的三角形数目大于120。

3.根据权利要求1所述的基于MDA和MLSSM的分析复杂目标电磁散射的方法，其特征在于，步骤3中最底层立方体的电尺寸为0.2～0.4个入射波波长。

4.根据权利要求1所述的基于MDA和MLSSM的分析复杂目标电磁散射的方法，其特征在于，步骤5中等效源的个数Q按照以下标准确定：

$Q\≥10\sqrt{k\frac{S_s{S}_o}{d}}$

其中S_s，S_o，d分别是远作用的时候场所在的立方体和源所在的立方体的电尺寸，以及它们之间的距离，k是入射波的波数，每一个等效源就是一个RWG基函数，而且这些等效源均匀分布在立方体的表面。

5.根据权利要求1所述的基于MDA和MLSSM的分析复杂目标电磁散射的方法，其特征在于，步骤5中利用等效源对阻抗矩阵的远作用子块进行低秩表示具体为：

$\left[Z_{\mathrm{mn}}\right]=\left[U\right]{\left[V\right]}^T$

其中[Z_mn] ∈C^m×n是阻抗矩阵中对应于远场作用的子块，[U] ∈C^m×r ，[V] ∈C^n×r，m、n、r为正整数，其中r小于min(m，n)。

6.根据权利要求1所述的基于MDA和MLSSM的分析复杂目标电磁散射的方法，其特征在于，步骤6中阻抗矩阵改进的MLSSM表示为：

$Z_l={\hat{Z}}_l+U_l{Z}_{(l-1)}{V}_l^H$

其中Z_L=Z是阻抗矩阵，

表示压缩后第二层表示，V_l=(U_l)^T是块对角的酉矩阵， l表示的是层数，0≤l≤L。

7.根据权利要求1所述的基于MDA和MLSSM的分析复杂目标电磁散射的方法，其特征在于，步骤7中构造迭代求解的快速矩阵矢量乘算法按照以下函数执行：

其中x,y分别为输入向量和输出向量，y₁，y₂，y₃为临时变量，U_l和是MLSSM表示中的矩阵，l表示的是层数。

Images