CN103279589B - 基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法 - Google Patents

Abstract

本发明公开了一种基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法，步骤如下：建立旋转对称体的母线，剖分离散得剖分线段；对剖分线段建立二叉树分组关系，并建立近场组、远场组之间的索引关系；建立旋转对称体电场积分方程，确定激励向量和总模式数；使用矩阵抽取法将每个最细层远场组模式阻抗矩阵分解成两个带条状矩阵分别存储；通过层组的索引将远场组阻抗矩阵分为父层的近场和父层的远场，并逐层嵌套处理，将远场阻抗矩阵分解成小矩阵形式；依次确定各个模式的模式阻抗矩阵，逐个求解模式阻抗矩阵的方程，由互易定理计算得到远场信息，完成仿真。该方法内存消耗低、方程迭代速度快，可高效仿真电大尺寸旋转对称体的电磁散射特性。

Description

基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法

技术领域

本发明涉及电磁仿真技术，特别是一种基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法。

背景技术

旋转对称体是指能由母线绕固定轴旋转一周得到的物体。作为一类重要的雷达目标，如导弹体，天线罩等，人们已经对其电磁散射特性进行了相当多的研究。特别是在跟踪和识别空间目标上，旋转体的电磁散射分析占据重要地位。利用旋转对称矩量法分析和研究各类旋转体的电磁散射问题一直是电磁散射研究中重要的课题。

Andreasen,M.G在1965年首先提出了旋转对称矩量法。文中将入射平面波利用傅立叶级数展开为相互正交的柱面波形式，利用各模式间的正交性，分别求解单一模式下的感应电流，然后进行线性叠加，从而求得散射场的分布(M.Andreasen,"Scatteringfrombodiesofrevolution,"AntennasandPropagation,IEEETransactionson,vol.13,pp.303-310,1965.)。但对建立的积分方程的求解采用的是点配法，计算简单但精度不高。J.R.MautzandR.F.Harrington利用等效电流的概念和电场边界条件建立了金属旋转对称体散射的电场积分方程(J.RMautzandR.F.Harrington,Generalizednetworkparametersforbodiesofrevolution,TechnicalReport,TR-687.SyracuseUniversity.Syracuse.NY13244.June1968)，将等效电流展开为关于φ的傅立叶（Fourier）级数和关于t的分段三角函数，然后利用数值方法中的矩量法（MoM）求解，解决了旋转金属体的散射问题。但当金属体形成谐振腔后，这种方法就失效了。J.R.Mautz.andR.F.Harrington又利用等效电流和等效磁流的概念和边界条件建立了金属散射体的混合场积分方程(J.R.MautzandR.F.Harrington,H-field,E-fieldandcombined-fieldsolutionsforconductingbodiesofrevolution,AEU321978,157-164)，解决了金属形成谐振腔后不能准确求解的问题。不论金属体是否形成了谐振腔，都能得到准确解。

在基于表面积分方程离散后形成的整个阻抗矩阵是稠密矩阵，用传统的矩量法直接求解时需要o(N²)存储量和计算复杂度（N为未知量的个数）。这严重阻碍了求解电大尺目标的散射问题。旋转对称体矩量法虽然将原来的问题进行了降维，但随着电尺寸的增加内存也会急剧上升。SDGedney,和RMittra使用快速傅里叶算法来加速旋转对称体矩量法计算（S.D.GedneyandR.Mittra,"TheuseoftheFFTfortheefficientsolutionoftheproblemofelectromagneticscatteringbyabodyofrevolution,"AntennasandPropagation,IEEETransactionson,vol.38,pp.313-322,1990.）。FFT算法加速了模式阻抗矩阵的形成，但却要一次存储所有模式阻抗矩阵，内存消耗很大。特别是斜入射时，内存的消耗为o(αN²)，其中α为总的模式数，N为母线剖分线段数。本发明将其利用旋转对称体模式阻抗矩阵远场低秩性，构造逐层嵌套压缩方法，用小的简易矩阵来等效原先的矩阵，内存和计算复杂度降低为O(NlogN)。

旋转对称体矩量法中的格林函数是以电流呈正弦分布的圆环为单位源产生的场，称为模式格林函数。斜入射时需要的模式数随着入射波的倾斜角和计算对象在柱坐标系下最大截面半径的电尺寸的增大而增大，模式格林函数也随之剧烈振荡，高斯积分需要的高斯积分点也随之增加，高模式下的阻抗矩阵生成速度变慢。文献1（Y.WenMing,etal.,"ClosedFormModalGreen'sFunctionsforAcceleratedComputationofBodiesofRevolution,"AntennasandPropagation,IEEETransactionson,vol.56,pp.3452-3461,2008.）使用闭式模式格林函数加速旋转对称体的计算是一种常用的克服积分振荡的方法，但该方法在分析电大尺寸的旋转对称体时会出现发散。

发明内容

本发明的目的在于提供一种快速稳定的基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法，该方法内存消耗低、矩阵填充和方程求解快。

实现本发明的技术解决方案是：一种基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法，步骤如下：

第1步，建立旋转对称体的母线，并将该母线按十分之一介质波长剖分离散得到剖分线段；

第2步，对母线的剖分线段建立二叉树分组关系，建立近场组和远场组之间的索引关系；

第3步，建立旋转对称体电场积分方程，使用旋转对称体矩量法离散电场积分方程，确定激励向量和总模式数；

第4步，使用矩阵抽取法将每个最细层远场组模式阻抗矩阵分解成两个带条状矩阵分别存储，近场组模式阻抗矩阵仍采用旋转对称体矩量法直接确定；

第5步，根据第4步得到的最细层两个带条状矩阵，通过层组的索引将远场组阻抗矩阵分为父层的近场和父层的远场，然后再将父层的远场分为父层的父层的近场和父层的父层的远场，这样逐层嵌套处理，接着再使用奇异值分解将远场阻抗矩阵分解成小矩阵的形式；

第6步，重复第4步～5步，依次确定各个模式的模式阻抗矩阵；

第7步，通过第第4～6步得到的嵌套型的远场阻抗矩阵，使用迭代方法逐个求解模式阻抗矩阵的方程，由互易定理计算得到远场信息，完成整个仿真过程。

本发明与现有技术相比其显著效果是：（1）相比于传统的旋转对称体矩量法，不用计算填充所有的基函数相互之间的作用元素，仅需要计算近场作用元素和抽取出来的远场作用元素，计算时间和计算速度提升；（2）相比于直接使用矩阵抽取方法计算，本发明方法构造的矩阵具有嵌套形式，内存消耗进一步降低，矩阵矢量乘的速度加快；（3）相比于其他类的快速算法，本发明方法构造方式简单，数学基础牢靠，计算精度可控。

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

附图说明

图1是本发明旋转对称体目标示意图。

图2是本发明二叉树近场组和远场组示意图。

图3是本发明二叉树分层示意图。

图4是本发明的嵌套处理的逐层处理流程图。

图5是本发明实施例1的计算结果图。

具体实施方式

下面结合附图和具体实施例对本发明的具体步骤进一步详细描述。

本发明是针对电大尺寸旋转对称体电磁散射特性分析的仿真平台，它基于旋转对称体矩量法和矩阵嵌套压缩方法，以达到对待求解快速求解目的，具体步骤如下：

第1步：建立旋转对称体的母线，并将该母线按十分之一介质波长剖分离散得到剖分线段。

结合图1，首先使用Ansys等商用软件建立旋转对称体的母线，并按照十分之一介质波长剖分离散，其中ε_r表示介质的相对介电常数，λ表示入射电磁波在空气中传播时的波长。旋转对称体表面的等效电磁流将展开成轴向方向φ的傅立叶（Fourier）级数和母线方向t的分段三角函数。通过这一步我们得到立体模型的一维剖分基本参数：剖分线段数量，剖分线段长度，剖分线段编号，剖分线段端点的编号和坐标。

第2步：对母线的剖分线段建立二叉树分组关系，建立近场组和远场组之间的索引关系；

如图2所示，用一个圆柱将旋转对称体包围住，该圆柱体定义为第零层组节点，把该圆柱体二等分形成第一层组节点，以此类推，直到最细层的圆柱高为1～2个波长之间。按照1个波长离散成10条母线剖分线段度量，则每个最细层组包含10～20个剖分线段。相邻的组为近场组，不相邻的组为远场组，再在此分组的基础上建立二叉树，如图3所示。近场组由旋转对称矩量法直接求解，远场组按照逐层嵌套的方法计算填充。

第3步：建立旋转对称体电场积分方程，使用旋转对称体矩量法离散电场积分方程，确定激励向量和总模式数；

3.1建立积分方程：

（1）电场积分方程EFIE

$\hat{t}\left(r\right)\·\underset{S}{\∫\∫}(1+\frac{1}{k^2}\▿\▿\·)J\left(r^{\′}\right)\frac{e^{-\mathrm{jkr}}}{4\mathrm{\πr}}{\mathrm{dr}}^{\′}=-\frac{j}{\mathrm{\ω\μ}}\hat{t}\left(r\right)\·E^{\mathrm{inc}}\left(r\right)---\left(1\right)$

（2）磁场积分方程MFIE

$\frac{J\left(r\right)}{2}+\frac{1}{4\mathrm{\π}}\underset{S-\mathrm{\δS}}{\∫\∫}\hat{n}[(r-r^{\′})\×J\left(r^{\′}\right)][1+\mathrm{jkr}]\frac{e^{-\mathrm{jk}|r-r^{\′}|}}{{|r-r^{\′}|}^3}{\mathrm{dr}}^{\′}=\hat{n}\×H^{\mathrm{inc}}\left(r\right)---\left(2\right)$

其中ω表示入射波角频率，μ表示空气磁导率，表示目标表面场点r的切向，表示目标表面法向，k为入射波波数，J(r)表示场点r的散射电流，J(r')表示源点r'的散射电流，E^inc(r)表示场点r的入射电场，H^inc(r)表示场点r的入射磁场。

若目标为开放结构，使用电场积分方程（1）计算；若目标为闭合结构使用混合场积分方程计算，所谓混合场积分方程CFIE为电场积分方程EFIE和磁场积分方程MFIE的线性叠加：

CFIE＝γEFIE+(1-γ)ηMFIE（3）

其中γ为混合参数，0≤γ≤1，η表示空气中的波阻抗。

3.2确定激励向量：

根据旋转对称特性，我们将目标表面的电流展开成一组基函数表示，该基函数沿着径向是分域基函数，而周向是完备正交的傅里叶级数：

$J\left(r\right)=\underset{\mathrm{\α}=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{a}_{\mathrm{\α},n}^t{f}_{\mathrm{\α},n}^t\left(r\right)+a_{\mathrm{\α},n}^{\mathrm{\φ}}{f}_{\mathrm{\α},n}^{\mathrm{\φ}}\left(r\right)---\left(4\right)$

α表示旋转对称体基函数的模式数，表示电流J(r)使用旋转对称体基函数展开对应于第α个模式的第n个基函数切向方向的展开系数，表示电流J(r)使用旋转对称体基函数展开对应于第α个模式的第n个基函数周向方向的展开系数。旋转对称体基函数表达式为：

$\begin{array}{c}{f}_{\mathrm{\α},n}^t\left(r\right)=f_n\left(t\right)e^{\mathrm{j\α\φ}}\hat{t}\left(r\right)\\ f_{\mathrm{\α},n}^{\mathrm{\φ}}\left(r\right)=f_n\left(t\right)e^{\mathrm{j\α\φ}}\widehat{\mathrm{\φ}}\left(r\right)\end{array},n=1,\·\·\·,N-1---\left(5\right)$

其中f_n(t)是径向展开函数，定义为：t表示r点的切向分量，T_n(t)表示三角基函数，ρ(t)表示r点到z轴的垂直距离，φ表示r点的周向角，e^jαφ表示傅里叶级数展开对应于第α个模式的指数项，N表示旋转对称体母线剖分段数，即旋转对称体基函数对应的未知量的个数。

根据矩量法原理，我们需要使用测试函数对积分方程两边进行测试形成阻抗矩阵，测试函数的表达式为：

$W_{\mathrm{\β},m}^t\left(r\right)=f_m{e}^{-\mathrm{j\β\φ}}\hat{t}\left(r\right),m=\mathrm{1,2},...,N-1$

$\left(6\right)$

$W_{\mathrm{\β},m}^{\mathrm{\φ}}\left(r\right)=f_m{e}^{-\mathrm{j\β\φ}}\widehat{\mathrm{\φ}}\left(r\right),\mathrm{\β}=1,\·\·\·,\mathrm{Mod}$

其中下标m表示基函数编号，β表示模式数编号，表示对应于第β个模式的第m个基函数切向方向的测试基函数，表示对应于第β个模式的第m个基函数周向方向的测试基函数，N表示旋转对称体母线剖分线段数，Mod表示旋转对称体需要的总模式数；

使用测试函数（6）对电场积分方程公式（1）两边测试并整理得到阻抗矩阵：

$Z=\left[\begin{array}{cc}{Z}^{\mathrm{tt}}& Z^{\mathrm{t\φ}}\\ Z^{\mathrm{\φt}}& Z^{\mathrm{\φ\φ}}\end{array}\right]---\left(7\right)$

每个子矩阵的元素（即表示子矩阵为Z^tt元素，表示子矩阵为Z^tφ元素，表示子矩阵为Z^φt元素，表示子矩阵为Z^φφ元素）的表达形式由如下公式给出：

$z_{\mathrm{mm}}^{\mathrm{tt}}=\underset{r=1}{\overset{M_r}{\mathrm{\Σ}}}\underset{s=1}{\overset{M_s}{\mathrm{\Σ}}}\left(T_r{T}_s\right[\sin {\mathrm{\γ}}_r\sin {\mathrm{\γ}}_s{G}_2+\cos {\mathrm{\γ}}_r\cos {\mathrm{\γ}}_s{G}_1]-\frac{1}{k^2}{\stackrel{\·}{T}}_r{\stackrel{\·}{T}}_s{G}_1)$

$z_{\mathrm{mn}}^{\mathrm{t\φ}}=j\underset{r=1}{\overset{M_r}{\mathrm{\Σ}}}\underset{s=1}{\overset{M_s}{\mathrm{\Σ}}}(\sin {\mathrm{\γ}}_r{T}_r{T}_s{G}_3+\frac{1}{k^2}\frac{\mathrm{\α}}{{\mathrm{\ρ}}_s}{\stackrel{\·}{T}}_r{T}_s{G}_1)$

$\left(8\right)$

$z_{\mathrm{mn}}^{\mathrm{\φt}}=j\underset{r=1}{\overset{M_r}{\mathrm{\Σ}}}\underset{s=1}{\overset{M_s}{\mathrm{\Σ}}}(\sin {\mathrm{\γ}}_s{T}_r{T}_s{G}_3+\frac{1}{k^2}\frac{\mathrm{\α}}{{\mathrm{\ρ}}_r}{T}_r{\stackrel{\·}{T}}_s{G}_1)$

$z_{\mathrm{mn}}^{\mathrm{\φ\φ}}=-\underset{r=1}{\overset{M_r}{\mathrm{\Σ}}}\underset{s=1}{\overset{M_s}{\mathrm{\Σ}}}(G_2-\frac{1}{k^2}\frac{{\mathrm{j\α}}^2}{{\mathrm{\ρ}}_r{\mathrm{\ρ}}_s}{G}_1)T_r{T}_s$

其中m表示测试函数编号，n表示基函数编号，α表示模式数，r表示测试函数对应的剖分线段编号、s表示基函数对应的剖分线段编号，T_r表示第r个剖分线段对应的三角基函数，T_s表示第s个剖分线段对应的三角基函数，M_r表示三角基函数T_r对应的母线剖分线段的个数，M_s表示三角基函数T_s对应的母线剖分线段的个数。表示三角基函数T_r在第r个剖分线段中心的散度，表示三角基函数T_s在第s个剖分线段中心的散度。其中表示第s个剖分线段的切向方向，表示第r个剖分线段的切向方向，表示第s个剖分线段的z轴方向，表示第r个剖分线段的z轴方向。ρ_r表示第r条剖分线段距z轴的垂直距离，ρ_s表示第s条剖分线段距z轴的垂直距离，G₁，G₂，G₃表示三种模式格林函数，其具体表达由下式给出：

$G_1={\mathrm{\Δ}}_r{\mathrm{\Δ}}_s{\∫}_0^{\mathrm{\π}}\cos \left({\mathrm{\α\φ}}^{\′}\right)\frac{e^{-\mathrm{jk}{R}_{\mathrm{rs}}}}{R_{\mathrm{rs}}}{\mathrm{d\φ}}^{\′}$

$G_2={\mathrm{\Δ}}_r{\mathrm{\Δ}}_s{\∫}_0^{\mathrm{\π}}\cos \left({\mathrm{\φ}}^{\′}\right)\cos \left({\mathrm{\α\φ}}^{\′}\right)\frac{e^{-\mathrm{jk}{R}_{\mathrm{rs}}}}{R_{\mathrm{rs}}}{\mathrm{d\φ}}^{\′}---\left(9\right)$

$G_3={\mathrm{\Δ}}_r{\mathrm{\Δ}}_s{\∫}_0^{\mathrm{\π}}\sin \left({\mathrm{\φ}}^{\′}\right)\sin \left({\mathrm{\α\φ}}^{\′}\right)\frac{e^{-{\mathrm{jkR}}_{\mathrm{rs}}}}{R_{\mathrm{rs}}}{\mathrm{d\φ}}^{\′}$

其中： $R_{\mathrm{rs}}=\sqrt{{({\mathrm{\ρ}}_r-{\mathrm{\ρ}}_s)}^2+{(z_r-z_s)}^2+2{\mathrm{\ρ}}_r{\mathrm{\ρ}}_s(1-\cos {\mathrm{\φ}}^{\′})}$ 表示场源点之间的距离，Δ_r表示测试函数对应的剖分线段的长度，Δ_s表示基函数对应的剖分线段的长度，φ表示周向角。公式（9）所示的三个格林函数我们称之为模式格林函数。

入射场为均匀平面波，设置入射角度为（θ，）及极化方式（垂直极化或水平极化）。入射电场矢量E^inc表达式为：E^inc＝E₀e^-jkr，E₀表示电场方向矢量，入射磁场矢量H^inc的表达式为：其中η为空气的波阻抗，表示传播方向。则激励向量V为：

$V=\left[\begin{array}{c}{V}^t\\ V^{\mathrm{\φ}}\end{array}\right]---\left(10\right)$

$V_{\mathrm{\β},m}^{t,\mathrm{\φ}}=-\frac{j}{\mathrm{\ω\μ}}{\∫}_{t_m}{\∫}_0^{2\mathrm{\π}}{f}_{\mathrm{\β},m}^{t,\mathrm{\φ}}\left(r\right)\·E^{\mathrm{inc}}\left(r\right)\mathrm{\ρd\φdt}$

其中V^t为激励向量母线切向方向的分量，V^φ为激励向量母线周向方向的分量，为激励向量V^t的元素，为激励向量V^φ的元素，ω表示入射波角频率，μ表示空气磁导率，表示母线切向方向和母线周向方向的β模式对应的第m个测试函数，E^inc(r)表示入射电磁场，ρ表示r点距z轴的垂直距离。

因为平面波可以用贝塞尔函数展开激励源，所以公式（10）可以由贝塞尔函数给出。

3.3确定最大模式数：

根据入射角度及旋转对称体的最大半径确定最大模式数：

Mod_max＝kR_max+1（11）

其中k为入射波波数，R_max为待求旋转对称体目标的最大半径。

建立相互独立的模式方程

Z_αI_α＝V_α（12）

其中α为模式数，I_α为第α个模式对应的电流系数解向量，-Mod_max≤α≤Mod_max垂直入射α取1，由于正负模式数之间存在对偶关系，一般仅需要求解非负模式方程即可。以上（7）～（10）给出的是电场积分方程的离散过程，磁场积分方程和混合场积分方程也可以完全类似的计算出来，这里不再赘述。

第4步，使用矩阵抽取法将每个最细层远场组模式阻抗矩阵分解成两个带条状矩阵分别存储，近场组模式阻抗矩阵仍采用旋转对称体矩量法直接确定；

若直接求解方程（12），其计算复杂度为o(N²)，当目标的电尺寸达到几百个电波长时，计算资源消耗大，使得难以在普通的个人计算机上完成。因此，我们使用第2步中的二叉树分组技术将目标分成各个小组，使用逐层嵌套压缩技术按组填充计算，假设目标被分成共L_max层，最细层即为第L_max层，如图3所示。

逐个模式数循环计算第L_max层组对应的模式阻抗矩阵，使用矩阵抽取算法将L_max层第α模式的阻抗矩阵分解成如下的表达形式：

$Z_{\mathrm{\α},L_{\max }}=Z_{\mathrm{\α},L_{\max }}^{\mathrm{near}}+Z_{\mathrm{\α},L_{\max }}^{\mathrm{far}}---\left(13\right)$

其中表示L_max层近场组第α个模式的阻抗矩阵，使用第三步中的公式直接计算得到。表示L_max层远场组第α个模式的阻抗矩阵，使用矩阵抽取算法计算得到的是一个稀疏矩阵的形式。

第5步，根据第4步得到最细层L_max层远场组第α个模式的阻抗稀疏矩阵通过层组的索引将远场组阻抗矩阵分为父层的近场和父层的远场，然后再将父层的远场分为父层的父层的近场和父层的父层的远场，这样逐层嵌套处理，接着再使用奇异值分解将远场阻抗矩阵分解成小矩阵的形式。

从最细层L_max开始循环递归到第2层，第l层远场第α个模式的阻抗矩阵Z_α,l可表示成如下嵌套形式：

$Z_{\mathrm{\α},l}={\tilde{Z}}_{\mathrm{\α},l}+P_{\mathrm{\α},l}{Z}_{\mathrm{\α},l-1}{Q}_{\mathrm{\α},l}^H---\left(14\right)$

其中Z_α,l是由第α个模式对应的第l+1层的远场组矩阵块降维后形成的第l层整个矩阵，当l＝L_max-1时， 表示第α个模式对应的第l层近场组之间的作用。方程右边第二项表示第α个模式对应的第l层的远场部分，P_α,l和分别是行变换矩阵和列变换矩阵，H表示共轭转置，并且在第l层中这两个矩阵都是通过对远场矩阵压缩以后得到的对角标准正交化矩阵。Z_α,l-1是计算得到的第α个模式对应的第l层的远场组矩阵块降维后形成的第l-1层整个矩阵，循环上述过程从而形成逐层嵌套压缩结构。结合附图4，给出逐层嵌套压缩结构形成的具体过程如下：

5.1）通过抽取法将最细层(第L_max层)的远场组第α个模式的阻抗矩阵分解为两个带条形矩阵，P、Q表示矩阵，上标表示矩阵的维数，下标表示层数。

5.2）将用奇异值分解得到三个矩阵相乘的形式：将整合得到矩阵P_l，将与作乘积形成一个的矩阵

5.3）将上述的进行奇异值分解得到新的三个矩阵相乘：将整合得到剩余的便是Z_α,l-1，这样就完成了第l层的处理。

5.4）对于Z_α,l-1可以分成两部分：表示l-1层近场的第l-1层远场的将进行奇异值分解得到将整理得到剩余的再一次进行奇异值分解，得到将整合得到P_l-1，剩下的部分便是Z_α,l-2，这样就完成了第l-1层的处理。即形成了第l层和第l-1层的嵌套压缩处理。

5.5）重复4）的过程直到第二层l＝2得到所有层的嵌套压缩关系（14）。

第6步，重复第4～5步，计算0到Mod个模式的模式阻抗矩阵。当模式阻抗矩阵的模式数α≥30时，则使用高震荡函数部分积分技术加速模式格林函数（9）的计算：

高震荡函数部分积分的原理如下，当α≥30时，具有以下形式的积分：

$I_1\left(\mathrm{\α}\right)={\∫}_a^{b}f\left(x\right)\cos \left(\mathrm{\αx}\right)\mathrm{dx}---\left(15\right)$

$I_2\left(\mathrm{\α}\right)={\∫}_a^{b}f\left(x\right)\sin \left(\mathrm{\αx}\right)\mathrm{dx}$

其中f(x)为积分核，b,a为积分上下限，积分上下限由模式格林函数的积分区间决定；α表示高震荡函数的模式即模式阻抗矩阵的模式数，这里对应为旋转对称体模式格林函数的模式数。I₁(α)，I₂(α)称为高振荡积分。这种形式的积分当α很大时可以由部分积分公式推出近似计算公式如下：

$I_1\left(\mathrm{\α}\right)={\∫}_a^{b}f\left(x\right)\cos \mathrm{\αxdx}$

$\≈\underset{k=0}{\overset{h-1}{\mathrm{\Σ}}}\frac{1}{{\mathrm{\α}}^{k+1}}[f^{\left(k\right)}\left(b\right)\sin (\frac{\mathrm{k\π}}{2}+\mathrm{\αb})-f^{\left(k\right)}\left(a\right)\sin (\frac{\mathrm{k\π}}{2}+\mathrm{\αa})]---\left(16\right)$

$I_2\left(\mathrm{\α}\right)={\∫}_a^{b}f\left(x\right)\sin \mathrm{\αxdx}$

$\≈\underset{k=0}{\overset{h-1}{\mathrm{\Σ}}}\frac{1}{{\mathrm{\α}}^{k+1}}[f^{\left(k\right)}\left(b\right)\cos (\frac{\mathrm{k\π}}{2}+\mathrm{\αb})-f^{\left(k\right)}\left(a\right)\cos (\frac{\mathrm{k\π}}{2}+\mathrm{\αa})]$

f^(k)(a)，f^(k)(b)表示f(x)的k次导数对应与x＝a,x＝b的值。h为截断项，一般取h＝5即可满足积分计算精度。

模式格林函数（9）具有高震荡积分的形式，因为当α＞30即模式数高于30时使用公式（16）计算，取截断h＝5。

第7步，通过第4～6步得到的嵌套型的远场阻抗矩阵，使用迭代方法逐个求解模式阻抗矩阵方程，由互易定理计算得到远场信息，完成整个仿真过程。

实施例1

根据本发明所述方法对一个半径为1m、高为20m的金属圆柱进行了仿真。入射波频率为3GHz，入射角度为（40°,0°），共需要43个模式数，母线剖分线段数为2200。其结果与原始的旋转对称算法吻合很好，证明了该方法的正确性，如图5所示。本发明还给出了本发明在时间消耗和内存消耗上的改变，表1为实施例1中计算效率与传统方法对比，从表中可以看出，内存消耗和计算时间都比传统的方法得到了节省。

表1

方法	内存消耗BOR	时间消耗BOR
			BoR-MoM	161.2M	69,512S
本发明方法	50.1M	18,214S

Claims (4)

1.一种基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法，其特征在于，步骤如下：

第1步，建立旋转对称体的母线，并将该母线按十分之一介质波长剖分离散得到剖分线段；

第2步，对母线的剖分线段建立二叉树分组关系，建立近场组、远场组之间的索引关系；

第3步，建立旋转对称体电场积分方程，使用旋转对称体矩量法离散电场积分方程，确定激励向量和总模式数；

第4步，使用矩阵抽取法将每个最细层远场组模式阻抗矩阵分解成两个带条状矩阵分别存储，近场组模式阻抗矩阵仍采用旋转对称体矩量法直接确定；

第5步，根据第4步得到的最细层两个带条状矩阵，通过层组的索引将远场组阻抗矩阵分为父层的近场和父层的远场，然后再将父层的远场分为父层的父层的近场和父层的父层的远场，这样逐层嵌套处理，接着再使用奇异值分解将远场阻抗矩阵分解成小矩阵的形式；

第6步，重复第4步～5步，依次确定各个模式的模式阻抗矩阵；

第7步，通过第4～6步得到的嵌套型的远场阻抗矩阵，使用迭代方法逐个求解模式阻抗矩阵的方程，由互易定理计算得到远场信息，完成整个仿真过程。

2.根据权利要求1所述的基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法，其特征在于，第2步中对母线的剖分线段建立二叉树分组关系，建立近场组和远场组之间的索引关系具体为：用一个圆柱将旋转对称体包围住，该圆柱体定义为第零层组节点，把该圆柱体二等分形成第一层组节点，以此类推，直到最细层的圆柱高为1～2个波长之间，所述波长为入射电磁波在空气中传播时的波长。

3.根据权利要求1所述的基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法，其特征在于，第5步中所述逐层嵌套处理具体过程为：

从最细层L_max开始循环递归到第2层，第l层远场第α个模式的阻抗矩阵Z_α,l可表示成如下嵌套形式：

$Z_{\mathrm{\α},l}={\tilde{Z}}_{\mathrm{\α},l}+P_{\mathrm{\α},l}{Z}_{\mathrm{\α},l-1}{Q}_{\mathrm{\α},l}^H$

其中α表示模式数的编号，l表示二叉树分层的编号，Z_α,l是由第α个模式对应的第l+1层的远场组矩阵块降维后形成的第l层整个矩阵，当l＝L_max-1时，Z_α,l＝Z_α,far，表示第α个模式对应的第l层近场组作用，方程右边第二项表示第α个模式对应的第l层的远场部分，P_α,l和分别是行变换矩阵和列变换矩阵，H表示共轭转置，并且在第l层中这两个矩阵都是通过对远场矩阵压缩以后得到的对角标准正交化矩阵，Z_α,l-1是计算得到的第α个模式对应的第l层的远场组矩阵块降维后形成的第l-1层整个矩阵，循环上述过程从而形成逐层嵌套压缩结构。

4.根据权利要求3所述的基于矩阵嵌套压缩的旋转对称体电磁散射特性仿真方法，其特征在于，第6步中所述确定各个模式的模式阻抗矩阵，当模式阻抗矩阵的模式数α≥30时，则使用高震荡函数部分积分技术加速模式格林函数的计算：

当α≥30时，具有以下形式的积分：

$I_1\left(\mathrm{\α}\right)={\∫}_a^{b}f\left(x\right)cos\left(\mathrm{\α}x\right)dx$

$I_2\left(\mathrm{\α}\right)={\∫}_a^{b}f\left(x\right)sin\left(\mathrm{\α}x\right)dx$

其中f(x)为积分核，b为积分上限，a为积分下限，I₁(α)，I₂(α)称为高振荡积分，由部分积分公式推出近似计算公式如下：

$\begin{array}{c}{I}_1\left(\mathrm{\α}\right)={\∫}_a^{b}f\left(x\right)cos\mathrm{\α}xdx\\ \≈\underset{k=0}{\overset{h-1}{\Σ}}\frac{1}{{\mathrm{\α}}^{k+1}}\[f^{\left(k\right)}\left(b\right)sin(\frac{k\mathrm{\π}}{2}+\mathrm{\α}b)-f^{\left(k\right)}\left(a\right)sin(\frac{k\mathrm{\π}}{2}+\mathrm{\α}a)\]\end{array}$

$\begin{array}{c}{I}_2\left(\mathrm{\α}\right)={\∫}_a^{b}f\left(x\right)\sin \mathrm{\α}xdx\\ \≈\underset{k=0}{\overset{h-1}{\Σ}}\frac{1}{{\mathrm{\α}}^{k+1}}\[f^{\left(k\right)}\left(b\right)\cos (\frac{k\mathrm{\π}}{2}+\mathrm{\α}b)-f^{\left(k\right)}\left(a\right)\cos (\frac{k\mathrm{\π}}{2}+\mathrm{\α}a)\]\end{array}$

f^(k)(a)表示f(x)的k次导数对应于x＝a的值，f^(k)(b)表示f(x)的k次导数对应于x＝b的值，h为截断项，取h＝5。