CN103425816A - 快速获取金属旋转对称体电磁散射特性的矩阵抽取方法 - Google Patents

Abstract

本发明公开了一种快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，建立金属旋转对称体的母线，按十分之一介质波长离散；采用二叉树结构对剖分线段进行分组，确定近场、远场组；使用具有傅里叶级数形式的旋转对称基函数展开未知散射电流并确实总模式数；对金属旋转对称体边界建立混合场积分方程，使用矩量法离散混合场积分方程，得到模式阻抗矩阵元素计算公式；按照近、远场的索引关系将模式阻抗矩阵分成近场模式阻抗矩阵和远场模式阻抗矩阵；将得到的各个模式电流系数累加得到金属旋转对称体表面的电流，由互易定理计算得到雷达散射截面积。本发明基于矩阵抽取方法，内存消耗低，矩阵填充和方程求解快。

Description

快速获取金属旋转对称体电磁散射特性的矩阵抽取方法

技术领域

本发明涉及电磁应用技术领域，特别涉及一种快速获取金属旋转对称体电磁散射特性的矩阵抽取方法。

背景技术

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

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

SD Gedney,和R Mittra使用快速傅里叶方法来加速金属旋转对称体矩量法计算（S.D.Gedney and R.Mittra,"The use of the FFT for the efficient solution of the problemof electromagnetic scattering by a body of revolution,"Antennas and Propagation,IEEE Transactions on,vol.38,pp.313-322,1990.）。FFT方法加速了模式矩阵的形成，但却要一次存储所有模式矩阵，内存消耗很大。

金属旋转对称体矩量法中的格林函数是以电流呈正弦分布的圆环为单位源产生的场，称为模式林函数。斜入射时需要的模式数随着入射波的倾斜角和计算对象在柱坐标系下最大截面半径的电尺寸的增大而增大，模式格林函数也随之剧烈振荡。使用闭式模式格林函数加速金属旋转对称体的计算是一种常用的克服积分振荡的方法。但该方法在仿真电大尺寸的金属旋转对称体时会出现发散（Y.Wen Ming,et al.,"Closed Form Modal Green's Functionsfor Accelerated Computation of Bodies of Revolution,"Antennas and Propagation,IEEETransactions on,vol.56,pp.3452-3461,2008.）。

现有仿真金属旋转对称体的方法存在以下问题：

（1）FTT方法加速金属旋转对称体计算但对内存的消耗大，仿真的金属旋转对称体的规模有限，只对部分结构十分有效；

（2）推导闭式模式格林函数来加速金属旋转对称体计算，操作繁琐，限制条件多；

（3）对于快速计算电大尺寸旋转对称目标斜入射时的散射场，高次模式矩阵填充效率低。

发明内容

针对现有技术中存在的以上问题，本发明的目的在于提供一种快速稳定的仿真电大尺寸金属旋转对称体目标电磁散射特性的方法，该方法基于矩阵抽取方法，内存消耗低，矩阵填充和方程求解快。

实现本发明目的的技术方案为：

一种快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，其包括步骤如下：

第一步，建立金属旋转对称体的母线，按十分之一介质波长离散，形成了母线的剖分线段；

第二步，采用二叉树结构对剖分线段进行分组，确定近场、远场组；

第三步，使用具有傅里叶级数形式的旋转对称基函数展开未知散射电流；

第四步，对金属旋转对称体边界建立混合场积分方程，使用矩量法离散混合场积分方程，得到模式阻抗矩阵方程；按照近、远场组将模式阻抗矩阵分成近场模式阻抗矩阵和远场模式阻抗矩阵；

第五步，采用矩量法直接计算近场组第α个模式矩阵；使用矩阵抽取方法填充远场组第α个模式阻抗矩阵；计算第α个模式对应激励向量；使用迭代方法求解第α个模式矩阵方程，得到第α个模式散射电流系数；上述求解过程从第0个模式数开始，直到α＝Mod；

第六步，将第五步得到的各个模式电流系数累加得到金属旋转对称体表面的电流，由互易定理计算得到雷达散射截面积。

本发明与现有技术相比，其显著优点：（1）避免了一次存储所有模式阻抗矩阵，降低了内存消耗低；（2）由于仅需要计算一小部分远场互作用的矩阵元素，因而即使不使用闭式格林函数也可以高效的生成阻抗矩阵，而高次模的低秩性更明显，十分适合使用本发明方法计算；（3）由于远场矩阵的低秩特性，迭代求解是矩阵矢量乘的速度加快，整个求解时间大大缩减；（4）方法有坚实的数学基础，鲁棒性好，易于编程实现，可以在个人计算机上仿真近千个波长的电大金属旋转对称体的电磁散射特性。

附图说明

图1是本发明具体实施例的金属旋转对称体目标示意图。

图2是本发明具体实施例的与待求金属旋转对称体共轴并能将其完全包围的圆柱结构示意图。

图3是本发明具体实施例的二叉树近场组和远场组示意图。

图4是本发明具体实施例的快速仿真金属旋转对称体电磁散射特性的方法流程图。

图5是采用本发明具体实施例的快速获取金属旋转对称体电磁散射特性的矩阵抽取方法对金属圆柱进行仿真示意图。

图6是本发明具体实施例的计算时间随未知量增加的复杂度曲线图。

图7是本发明具体实施例的内存消耗随未知量增加的复杂度曲线图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明，使本发明的上述及其它目的、特征和优势将更加清晰。

本发明一种快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，包括步骤如下：

第一步，建立金属旋转对称体的母线，按十分之一介质波长离散，形成了母线的剖分线段；

第二步，采用二叉树结构对剖分线段进行分组，确定近场、远场组；

第三步，使用具有傅里叶级数形式的旋转对称基函数展开未知散射电流；

第四步，对金属旋转对称体边界建立混合场积分方程，使用矩量法离散混合场积分方程，得到模式阻抗矩阵方程；按照近、远场组将模式阻抗矩阵分成近场模式阻抗矩阵和远场模式阻抗矩阵；

第五步，采用矩量法直接计算近场组第α个模式矩阵；使用矩阵抽取方法填充远场组第α个模式阻抗矩阵；计算第α个模式对应激励向量（说明书中相应位置加以说明）；使用迭代方法求解第α个模式矩阵方程，得到第α个模式散射电流系数；上述求解过程从第0个模式数开始，直到α＝Mod；

第六步，将第五步得到的各个模式电流系数累加得到金属旋转对称体表面的电流，由互易定理计算得到雷达散射截面积。

图1是本发明具体实施例的金属旋转对称体目标示意图。图2是本发明具体实施例的与待求金属旋转对称体共轴并能将其完全包围的圆柱结构示意图。图3是本发明具体实施例的二叉树近场组和远场组示意图。图4是本发明具体实施例的快速仿真金属旋转对称体电磁散射特性的方法流程图。

结合图1至图4，本发明是针对电大金属旋转对称体的仿真平台，它基于旋转对称矩量法和矩阵抽取方法，以达到对待求解散射体快速仿真的目的。具体步骤如下：

第一步：建立金属旋转对称体的母线，按十分之一介质波长离散，形成了母线的剖分线段；

首先使用Ansys等商用软件建立金属旋转对称体的母线，并按照λ/10对母线进行离散，其中λ表示波长，ε_r表示相对介电常数。金属旋转对称体表面的等效电磁流将展开成周向

方向向的傅立叶（Fourier）级数和母线方向的分段三角函数。通过这一步我们得到立体模型的一维剖分基本参数：剖分的线段数，线段长度，线段编号，线段端点的编号和坐标。

第二步：采用二叉树结构对剖分线段进行分组，确定近场、远场组；

用一个虚拟圆柱将金属旋转对称体包围住，如图2所示，该圆柱体为第零层组节点，对该圆柱体二等分，形成的每个子圆柱体为第一层组节点，然后再对每个子圆柱体二等分，直到圆柱高为1到2个波长之间，最后形成的子圆柱体为最细层组节点，每个组节点即为一组；

同一层相邻的组为近场组，不相邻的组为远场组，再在此分组的基础上建立二叉树，子层的远场组为父层的近场组。

第三步：使用具有傅里叶级数形式的旋转对称基函数展开未知散射电流；

在金属旋转对称体表面建立局部坐标系，

是表面法向单位矢量，

是沿方位角方向的单位矢量，

是沿母线方向的单位矢量，满足

n表示基函数编号，m表示测试函数编号；利用旋转对称体结构特性，将表面r'点散射电磁流J(r')表示为：

$J\left(r^{\′}\right)=\underset{\mathrm{\α}=-\mathrm{Mod}}{\overset{\mathrm{Mod}}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(I_{\mathrm{\αn}}^t{f}_{\mathrm{\αn}}^t\left(r^{\′}\right)+I_{\mathrm{\αn}}^{\mathrm{\φ}}{f}_{\mathrm{\αn}}^{\mathrm{\φ}}\left(r^{\′}\right))---\left(1\right)$

分别表示第α个模式数对应的方向和

方向第n个基函数，由以下公式给出，

分别是对应的基函数展开系数，N为总基函数个数；

$f_{\mathrm{\αn}}^t\left(r^{\′}\right)=\frac{T_n\left(t\right)}{\mathrm{\ρ}\left(t\right)}{e}^{\mathrm{j\α\φ}}\hat{t}\left(r^{\′}\right)$ $\left(2\right)$

$f_{\mathrm{\αn}}^{\mathrm{\φ}}\left(r^{\′}\right)=\frac{T_n\left(t\right)}{\mathrm{\ρ}\left(t\right)}{e}^{\mathrm{j\α\φ}}\widehat{\mathrm{\φ}}\left(r^{\′}\right)$

其中T_n(t)为第n个三角基函数，有公式给出；t为剖分线段上的一点在母线

方向的取值，ρ(t)是剖分线段上t处对应的半径值。

其中

和

表示第n条剖分线段的起点和终点，

和

表示第n+1条剖分线段的起点和终点，t为剖分线段上的一点在母线

方向的取值；

根据旋转对称体的最大半径R_max、入射频率f和入射俯仰角θ^inc确定总模式数Mod：

Mod_max＝k_ρR_max+1          (4)

其中

μ为空气磁导率，ε为空气介电常数。若入射角度垂直于旋转轴，则总模式数为1；

第四步：对金属旋转对称体边界建立混合场积分方程，使用矩量法离散混合场积分方程，得到模式阻抗矩阵；按照近、远场的索引关系将模式阻抗矩阵分成近场模式阻抗矩阵和远场模式阻抗矩阵；

建立积分方程：

（1）电场积分方程EFIE

$\hat{t}\left(r\right)\·\∫\underset{S}{\∫}(1+\frac{1}{k^2}\▿\▿\·)J\left(r^{\′}\right)\frac{e^{-\mathrm{jk}|r-r^{\′}|}}{4\mathrm{\π}|r-r^{\′}|}d{r}^{\′}=-\frac{j}{\mathrm{\ω\μ}}\hat{t}\left(r\right)\·E^{\mathrm{inc}}\left(r\right)---\left(5\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(6\right)$

其中k为入射波波数，J(r)表示场点r的散射电流，J(r')表示源点r'的散射电流，E^inc(r)表示场点r的入射电场，H^inc(r)表示场点r的入射磁场，

ω表示入射波角频率。

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

CFIE＝γEFIE+(1-γ)ηMFIE    (7)

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

使用公式（13）基函数和其共轭作为测试函数，应用矩量法基本原理对积分方程进行离散，得到一些列模式阻抗矩阵方程：

Z_αI_α=V_α        (8)

其中Z_α为模式α对应的模式阻抗矩阵，I_α为模式α对应的电流系数，V_α为模式α对应的激励向量。

模式阻抗矩阵包含四个子模块：

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

其中各个子模块的矩阵元素由下面公式给出：

$z_{\mathrm{mn}}^{\mathrm{pq}}=\mathrm{\γ}\left\{{\∫\∫}_{f_{\mathrm{\αm}}^{t,\mathrm{\φ}}}{\∫\∫}_{f_{\mathrm{\αn}}^{t,\mathrm{\φ}}}(f_{\mathrm{\αm}}^{t,\mathrm{\φ}}\left(r\right)\·f_{\mathrm{\αn}}^{t,\mathrm{\φ}}\left(r^{\′}\right)-\frac{1}{4\mathrm{\π}{k}^2}[\▿\·f_{\mathrm{\αm}}^{t,\mathrm{\φ}}\left(r\right)\left]\right[{\▿}^{\′}\·f_{\mathrm{\αn}}^{t,\mathrm{\φ}}\left(r^{\′}\right)\left]\right)\frac{e^{-\mathrm{jkr}}}{r}{\mathrm{dr}}^{\′}\mathrm{dr}\right\}$

$+(1-\mathrm{\γ})\mathrm{\η}\left\{\begin{array}{c}\frac{1}{2}{\∫\∫}_{f_{\mathrm{\αm}}^{t,\mathrm{\φ}}}{f}_{\mathrm{\αm}}^{t,\mathrm{\φ}}\left(r\right)\·f_{\mathrm{\αn}}^{t,\mathrm{\φ}}\left(r^{\′}\right)\mathrm{dr}\\ +\frac{1}{4\mathrm{\π}}{\∫\∫}_{f_{\mathrm{\αm}}^{t,\mathrm{\φ}}}{\∫\∫}_{f_{\mathrm{\αn}}^{t,\mathrm{\φ}}}{f}_{\mathrm{\αm}}^{t,\mathrm{\φ}}\left(r\right)\·\hat{n}\left(r\right)\×[(r-r^{\′})\×f_{\mathrm{\αn}}^{t,\mathrm{\φ}}\left(r^{\′}\right)]\·[1+\mathrm{jkr}]\frac{e^{-\mathrm{jkr}}}{r^3}{\mathrm{dr}}^{\′}\mathrm{dr}\end{array}\right\}---\left(10\right)$

其中pq表示t,φ，即代表选择基函数或测试函数的上标t,φ中任两种组合。对应于模式阻抗矩阵中的四个不同的子模块；

根据近、远场组关系可将（20）式分解为：

Z_α＝Z_α ^near+Z_α ^far        (11)

其中Z_α ^near表示近场模式阻抗矩阵，Z_α ^far表示远场模式阻抗矩阵；

第五步，采用矩量法直接计算近场组第α个模式矩阵；使用矩阵抽取方法填充远场组第α个模式阻抗矩阵；计算第α个模式对应激励向量；使用迭代方法求解第α个模式矩阵方程，得到第α个模式散射电流系数；上述求解过程从第0个模式数开始，直到α＝Mod；

5.1采用矩量法直接计算近场组第α个模式矩阵；

Z_α ^near近场模式阻抗矩阵使用传统的旋转对称体矩量法公式（22）逐个元素填充，并按照稀疏格式存储；

5.2使用矩阵抽取方法填充远场组第α个模式阻抗矩阵；

远场组之间的作用Z_α ^far采用矩阵抽取方法填充。而Z_α ^far又包含四个小矩阵Z_α ^ttfar,Z_α ^tφfar,Z_α ^φtfar和Z_α ^φφfar，先对Z_α ^ttfar使用矩阵抽取方法填充。用

矩阵

代表矩量法中两个远场组的相互作用阻抗矩阵Z_α ^ttfar。矩阵抽取方法通过两个满秩矩阵的乘积形式来构造近似矩阵

用

来近似估计

即

${\tilde{X}}^{\stackrel{\‾}{m}\×\stackrel{\‾}{n}}=P^{\stackrel{\‾}{m}\×r}{Q}^{r\×\stackrel{\‾}{n}}=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{p}_i^{\stackrel{\‾}{m}\×1}{q}_i^{1\×\stackrel{\‾}{n}}---\left(12\right)$

其中r为矩阵

的有效秩，

为场组中包含的未知量的个数，

为源组中包含的未知量的个数，

和

为两个满秩矩阵。因此，矩阵抽取近似方法只需要保存

和

这两个较小的矩阵即可，这样存储量就由

降到

计算复杂度由原来的O(N²)，经过新的矩阵矢量乘可以降低到2rN，其中N为未知量个数。

设C=[C₁,…,C_r]和L=[L₁,…,L_r]为包含从矩阵

中挑选出的行和列的索引的数组，p_k为矩阵P的第k列，q_k为矩阵Q的第k行，其中代表矩阵

的第C₁行，

为矩阵

在第k次循环得到的近似矩阵；

所述矩阵抽取方法如下：

初始化循环初始值：

初始化第一个行索引C₁=1，令X=0。

初始化近似误差矩阵的第一行： $\tilde{E}(C_1,:)=X(C_1,:).$

在第一行中寻找最大值从而确定第一个列索引L₁： $\left|\tilde{E}(C_1,L_1)\right|={\max }_l\left(\right|\tilde{E}(C_1,l)\left|\right).$

得到Q矩阵的第一行： $q_1=\tilde{E}(C_1,:)/\tilde{E}(C_1,L_1).$

初始化近似误差矩阵的第一列： $\tilde{E}(:,L_1)=X(:,L_1).$

得到P矩阵的第一列： $p_1=\tilde{E}(:,L_1).$

${\left|\right|{\tilde{X}}^{\left(1\right)}\left|\right|}^2={\left|\right|{\tilde{X}}^{\left(0\right)}\left|\right|}^2+{\left|\right|p_1\left|\right|}^2{\left|\right|q_1\left|\right|}^2.$

在第一列中寻找最大值从而确定第二个行索引C₂：

$\left|\tilde{E}(C_2,L_1)\right|={\max }_c\left(\right|\tilde{E}(c,L_1)\left|\right),c\≠C_1.$

第k次循环：

更新近似误差矩阵的第I_k行： $\tilde{E}(C_k,:)=X(C_k,:)-{\mathrm{\Σ}}_{i=1}^{k-1}{\left(p_i\right)}_{C_k}{q}_i.$

在第C_k行中找到最大值从而确定第k个列索引L_k：

$\left|\tilde{E}(C_k,L_k)\right|={\max }_l\left(\right|\tilde{E}(C_k,l)\left|\right),l\≠L_1,\·\·\·,L_{k-1}.$

得到Q矩阵的第ｋ行： $q_k=\tilde{E}(C_k,:)/\tilde{E}(C_k,L_k).$

更新近似误差矩阵的第L_k列： $\tilde{E}(:,L_k)=X(:,L_k)-{\mathrm{\Σ}}_{i=1}^{k-1}{\left(q_i\right)}_{L_k}{p}_i.$

得到P矩阵的第ｋ列： $p_k=\tilde{E}(:,L_k).$

${\left|\right|{\tilde{X}}^{\left(k\right)}\left|\right|}^2={\left|\right|{\tilde{X}}^{(k-1)}\left|\right|}^2+2{\mathrm{\Σ}}_{l=1}^{k-1}\left|p_l^T{p}_k\right|\·\left|q_l^T{q}_k\right|+{\left|\right|p_k\left|\right|}^2{\left|\right|q_k\left|\right|}^2.$

判断收敛误差：如果

迭代结束，否则找下一个行索引I_k+1： $\left|\tilde{E}(C_{k+1},L_k)\right|={\max }_i\left(\right|\tilde{E}(c,L_k)\left|\right),c\≠C_1,\·\·\·,C_k.$ 进入第k+1次循环，直至停止条件满足。

通过上述操作得到矩阵P和矩阵Q，从而得到了矩阵的近似矩阵，完成了一个子远场模式矩阵的计算。重复上述行列矩阵抽取方法的过程，完成所有远场子模式矩阵Z_α ^ttfar的填充。

其它三个矩阵Z_α ^tφfar,Z_α ^φtfar,Z_α ^φφfar按照Z_α ^ttfar得到的行列信息[C₁,…,C_r]和[L₁,…,L_r]进行抽取。重复上述矩阵抽取方法逐层计算各个远场组之间互作用远场模式阻抗矩阵，完成第α个模式的填充过程。

5.3计算第α个模式对应激励向量V_α；

设入射场为均匀平面波，设置入射角度为（θ^inc，φ^inc）及极化方式（垂直极化或水平极化）。入射电场矢量E^inc表达式为：E^inc＝E₀e^-jkr，E₀表示电场方向矢量，入射磁场矢量H^inc的表达式为：

其中η为空气的波阻抗，

表示传播方向。则激励向量V_α为：

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

其中激励向量中对应的元素为：

$V_{\mathrm{\αm}}^{t,\mathrm{\φ}}=\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}]$ $\left(14\right)$

$+(1-\mathrm{\γ})[{\∫}_{t_m}{\∫}_0^{2\mathrm{\π}}{f}_{\mathrm{\α},m}^{t,\mathrm{\φ}}\left(r\right)\·\hat{n}\×\hat{k}\×E^{\mathrm{inc}}\mathrm{\ρd\φdt}]$

其中，ρ表示r点距z轴的垂直距离；

5.4使用迭代方法求解第α个模式矩阵方程，得到第α个模式散射电流系数I_α；

GMRes为一种常用的迭代方法，这里被用来求解模式阻抗矩阵（20）；

重复上述5.1-5.4，逐个求解所有模式阻抗矩阵方程，得到各个模式电流系数；随着模式数的增高模式矩阵特性值整体变小，矩阵抽取行列时需要降低截断精度，否则计算效率下降，一般在模式数大于10以后可以取0.1以上。

第六步，将第五步得到的各个模式电流系数使用公式（12）累加得到金属旋转对称体表面的电流，由互易定理计算得到雷达散射截面积σ。

由互易定理可以推导得到每平方米的雷达散射截面积的θθ,θφ,φθ,φφ四种可能的情况下对应的公式：

$\frac{{\mathrm{\σ}}^{\mathrm{\θ\θ}}}{{\mathrm{\λ}}^2}=\frac{k^4}{16{\mathrm{\π}}^3{\mathrm{\η}}^2}{|V_0^{\mathrm{t\θ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_0^{\mathrm{t\θ}}+2\underset{\mathrm{\α}=1}{\overset{\mathrm{Mod}}{\mathrm{\Σ}}}(V_{\mathrm{\α}}^{\mathrm{t\θ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_{\mathrm{\α}}^{\mathrm{t\θ}}-V_{\mathrm{\α}}^{\mathrm{\φ\θ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_{\mathrm{\α}}^{\mathrm{\φ\θ}})\cos \left(\mathrm{\α}({\mathrm{\φ}}^{\mathrm{sca}}-{\mathrm{\φ}}^{\mathrm{inc}})\right)|}^2$

$\frac{{\mathrm{\σ}}^{\mathrm{\φ\θ}}}{{\mathrm{\λ}}^2}=\frac{k^4}{16{\mathrm{\π}}^3{\mathrm{\η}}^2}{\left|2\underset{\mathrm{\α}=1}{\overset{\mathrm{Mod}}{\mathrm{\Σ}}}(-V_{\mathrm{\α}}^{\mathrm{t\φ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_{\mathrm{\α}}^{\mathrm{t\θ}}+V_{\mathrm{\α}}^{\mathrm{\φ\φ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_{\mathrm{\α}}^{\mathrm{\φ\θ}})\sin \left(\mathrm{\α}({\mathrm{\φ}}^{\mathrm{sca}}-{\mathrm{\φ}}^{\mathrm{inc}})\right)\right|}^2$ $\left(15\right)$

$\frac{{\mathrm{\σ}}^{\mathrm{\θ\φ}}}{{\mathrm{\λ}}^2}=\frac{k^4}{16{\mathrm{\π}}^3{\mathrm{\η}}^2}{\left|2\underset{\mathrm{\α}=1}{\overset{\mathrm{Mod}}{\mathrm{\Σ}}}(V_{\mathrm{\α}}^{\mathrm{t\θ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_{\mathrm{\α}}^{\mathrm{t\φ}}-V_{\mathrm{\α}}^{\mathrm{\φ\θ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_{\mathrm{\α}}^{\mathrm{\φ\φ}})\sin \left(\mathrm{\α}({\mathrm{\φ}}^{\mathrm{sca}}-{\mathrm{\φ}}^{\mathrm{inc}})\right)\right|}^2$

$\frac{{\mathrm{\σ}}^{\mathrm{\φ\φ}}}{{\mathrm{\λ}}^2}=\frac{k^4}{16{\mathrm{\π}}^3{\mathrm{\η}}^2}{|V_0^{\mathrm{\φ\φ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_0^{\mathrm{\φ\φ}}+2\underset{\mathrm{\α}=1}{\overset{\mathrm{Mod}}{\mathrm{\Σ}}}(-V_{\mathrm{\α}}^{\mathrm{t\φ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_{\mathrm{\α}}^{\mathrm{t\φ}}+V_{\mathrm{\α}}^{\mathrm{\φ\φ}}\left({\mathrm{\θ}}^{\mathrm{sca}}\right)I_{\mathrm{\α}}^{\mathrm{\φ\φ}})\cos \left(\mathrm{\α}({\mathrm{\φ}}^{\mathrm{sca}}-{\mathrm{\φ}}^{\mathrm{inc}})\right)|}^2$

其中λ表示波长，(θ^sca,φ^sca)为远区场方向，上标

分别代表θ或φ，

表示散射场的极化方向，表示入射场的极化方向，

代表

或

表示

方向极波（θ或φ方向）俯仰角θ^sca时入射波α模式

方向分量激励向量，

表示α=0，

表示

方向极化波入射时对应的第α模式方向分量散射电流系数。

实施例：

图5是采用本发明具体实施例的快速获取金属旋转对称体电磁散射特性的矩阵抽取方法对金属圆柱进行仿真示意图，根据本发明所述方法对一个半径为1m高为100m的金属圆柱进行了仿真，入射场设置为均匀平面波，频率为0.3GHz，入射角度为（0°，0°），平行极化波，模式数为1。其结果与原始的旋转对称方法吻合很好，证明了该方法的正确性。图6是本发明具体实施例的计算时间随未知量增加的复杂度曲线图。图7是本发明具体实施例的内存消耗随未知量增加的复杂度曲线图，两者皆达到了O(N)，充分显示了其高效性。

Claims (6)

1.一种快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，其特征在于步骤如下：

第一步，建立金属旋转对称体的母线，按十分之一介质波长离散，形成了母线的剖分线段；

第二步，采用二叉树结构对剖分线段进行分组，确定近场、远场组；

第三步，使用具有傅里叶级数形式的旋转对称基函数展开未知散射电流；

第四步，对金属旋转对称体边界建立混合场积分方程，使用矩量法离散混合场积分方程，得到模式阻抗矩阵方程；按照近、远场组将模式阻抗矩阵分成近场模式阻抗矩阵和远场模式阻抗矩阵；

第五步，采用矩量法直接计算近场组第α个模式矩阵；使用矩阵抽取方法填充远场组第α个模式阻抗矩阵；计算第α个模式对应激励向量；使用迭代方法求解第α个模式矩阵方程，得到第α个模式散射电流系数；上述求解过程从第0个模式数开始，直到α＝Mod；

第六步，将第五步得到的各个模式电流系数累加得到金属旋转对称体表面的电流，由互易定理计算得到雷达散射截面积。

2.根据权利要求1所述的快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，其特征在于：所述第二步中采用二叉树结构对剖分线段进行分组为：用一个虚拟圆柱将金属旋转对称体包围住，该圆柱体为第零层组节点，对该圆柱体二等分，形成的每个子圆柱体为第一层组节点，然后再对每个子圆柱体二等分，直到圆柱高为1到2个波长之间，最后形成的子圆柱体为最细层组节点，每个组节点即为一组；

所述确定近场、远场组为：同一层相邻的组为近场组，不相邻的组为远场组，再在此分组的基础上建立二叉树，子层的远场组为父层的近场组。

3.根据权利要求1所述的快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，其特征在于，所述第三步的具体步骤如下：

根据旋转对称体的最大半径R_max、入射频率f和入射俯仰角θ^inc确定总模式数Mod；

在金属旋转对称体表面建立局部坐标系，

是表面法向单位矢量，

是沿方位角方向的单位矢量，

是沿母线方向的单位矢量，满足

利用旋转对称体结构特性，将表面散射电磁流J表示为：

$J=\underset{\mathrm{\α}=-\mathrm{Mod}}{\overset{\mathrm{Mod}}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(I_{\mathrm{\αn}}^t{f}_{\mathrm{\αn}}^t+I_{\mathrm{\αn}}^{\mathrm{\φ}}{f}_{\mathrm{\αn}}^{\mathrm{\φ}})---\left(1\right)$

分别表示第α个模式数对应的

方向和

方向第n个基函数，由以下公式给出，

分别是对应的基函数展开系数，N为总基函数个数；

$f_{\mathrm{\αn}}^t=\hat{t}\frac{T_n\left(t\right)}{\mathrm{\ρ}\left(t\right)}{e}^{\mathrm{j\α\φ}}$ $\left(2\right)$

$f_{\mathrm{\αn}}^{\mathrm{\φ}}=\widehat{\mathrm{\φ}}\frac{T_n\left(t\right)}{\mathrm{\ρ}\left(t\right)}{e}^{\mathrm{j\α\φ}}$

其中T_n(t)为第n个三角基函数，由以下公式给出；t为剖分线段上的一点在母线

方向的取值，ρ(t)是剖分线段上t处对应的半径值；

其中

和

表示第n条剖分线段的起点和终点，

和

表示第n+1条剖分线段的起点和终点；

根据旋转对称体的最大半径R_max、入射频率f和入射俯仰角θ^inc确定总模式数Mod：

Mod＝k_ρR_max+1          (4)

其中

μ为空气磁导率，ε为空气介电常数；若入射角度垂直于旋转轴，则总模式数为1。

4.根据权利要求1所述的快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，其特征在于，所述第四步的具体步骤如下：

建立积分方程：

（1）电场积分方程EFIE

$\hat{t}\left(r\right)\·\∫\underset{S}{\∫}(1+\frac{1}{k^2}\▿\▿\·)J\left(r^{\′}\right)\frac{e^{-\mathrm{jk}|r-r^{\′}|}}{4\mathrm{\π}|r-r^{\′}|}d{r}^{\′}=-\frac{j}{\mathrm{\ω\μ}}\hat{t}\left(r\right)\·E^{\mathrm{inc}}\left(r\right)---\left(5\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(6\right)$

其中k为入射波波数，J(r)表示场点r的散射电流，J(r')表示源点r'的散射电流，E^inc(r)表示场点r的入射电场，H^inc(r)表示场点r的入射磁场；

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

CFIE＝γEFIE+(1-γ)ηMFIE        (7)

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

使用公式(2)基函数和其共轭作为测试函数，应用矩量法基本原理对积分方程进行离散，得到一系列模式阻抗矩阵方程：

Z_αI_α=V_α         (8)

其中Z_α为模式α对应的模式阻抗矩阵，I_α为模式α对应的电流系数，V_α为模式α对应的激励向量；

模式阻抗矩阵包含四个子模块：

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

根据近、远场组将(9)式分解为：

Z_α＝Z_α ^near+Z_α ^far           (10)

其中Z_α ^near表示近场模式阻抗矩阵，Z_α ^far表示远场模式阻抗矩阵。

5.根据权利要求1所述的快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，其特征在于，所述第五步步骤如下：

5.1采用矩量法直接计算近场组第α个模式矩阵；

Z_α ^near近场模式阻抗矩阵使用旋转对称体矩量法公式逐个元素填充，使用稀疏存储；

5.2使用矩阵抽取方法填充远场组第α个模式阻抗矩阵；

远场模式阻抗矩阵Z_α ^far采用矩阵抽取方法填充，而Z_α ^far又包含四个小矩阵Z_α ^ttfar,Z_α ^tφfar,Z_α ^φtfar和Z_α ^φφfar，先对Z_α ^ttfar使用矩阵抽取方法填充，用

矩阵

代表矩量法中两个远场组的相互作用阻抗矩阵Z_α ^ttfar；矩阵抽取方法通过两个满秩矩阵的乘积形式来构造近似矩阵

用

来近似估计

即

${\tilde{X}}^{\stackrel{\‾}{m}\×\stackrel{\‾}{n}}=P^{\stackrel{\‾}{m}\×r}{Q}^{r\×\stackrel{\‾}{n}}=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{p}_i^{\stackrel{\‾}{m}\×1}{q}_i^{1\×\stackrel{\‾}{n}}---\left(11\right)$

其中r为矩阵

的有效秩，为场组中包含的未知量的个数，为源组中包含的未知量的个数，

和

为两个满秩矩阵；矩阵抽取方法只需要保存

和

这两个较小的矩阵；

设C=[C₁,…,C_r]和L=[L₁,…,L_r]为包含从矩阵

中挑选出的行和列的索引的数组，p_k为矩阵P的第k列，q_k为矩阵Q的第k行，其中

代表矩阵

的第C₁行，

为矩阵

在第k次循环得到的近似矩阵；

通过矩阵抽取方法得到矩阵P和矩阵Q，从而得到了矩阵的近似矩阵，完成一个子远场模式矩阵的计算；重复行列矩阵抽取方法的过程，完成所有远场子模式矩阵Z_α ^ttfar的填充；

其它三个矩阵Z_α ^tφfar,Z_α ^φtfar,Z_α ^φφfar按照Z_α ^ttfar得到的行列信息[C₁,…,C_r]和[L₁,…,L_r]进行抽取；重复矩阵抽取方法逐层计算各个远场组之间互作用远场模式阻抗矩阵，完成第α个模式的填充过程；

5.3计算第α个模式对应激励向量V_α；

5.4使用迭代方法求解第α个模式矩阵方程，得到第α个模式散射电流系数I_α；

重复上述5.1-5.4，逐个求解所有模式阻抗矩阵方程，得到各个模式电流系数；随着模式数的增高模式矩阵特性值整体变小，矩阵抽取行列时需要降低截断精度，否则计算效率下降，在模式数大于10以后可取0.1以上。

6.根据权利要求1或5所述的快速获取金属旋转对称体电磁散射特性的矩阵抽取方法，其特征在于，所述步骤5.2中的矩阵抽取方法如下：

初始化循环初始值：

初始化第一个行索引C₁=1，令Χ=0；

初始化近似误差矩阵的第一行： $\tilde{E}(C_1,:)=X(C_1,:);$

在第一行中寻找最大值从而确定第一个列索引L₁：

$\left|\tilde{E}(C_1,L_1)\right|={\max }_l\left(\right|\tilde{E}(C_1,l)\left|\right);$

得到Q矩阵的第一行： $q_1=\tilde{E}(C_1,:)/\tilde{E}(C_1,L_1);$

初始化近似误差矩阵的第一列： $\tilde{E}(:,L_1)=X(:,L_1).$

得到P矩阵的第一列： $p_1=\tilde{E}(:,L_1);$

${\left|\right|{\tilde{X}}^{\left(1\right)}\left|\right|}^2={\left|\right|{\tilde{X}}^{\left(0\right)}\left|\right|}^2+{\left|\right|p_1\left|\right|}^2{\left|\right|q_1\left|\right|}^2;$

在第一列中寻找最大值从而确定第二个行索引C₂：

$\left|\tilde{E}(C_2,L_1)\right|={\max }_c\left(\right|\tilde{E}(c,L_1)\left|\right),c\≠C_1;$

第k次循环：

更新近似误差矩阵的第I_k行： $\tilde{E}(C_k,:)=X(C_k,:)-{\mathrm{\Σ}}_{i=1}^{k-1}{\left(p_i\right)}_{C_k}{q}_i;$

在第C_k行中找到最大值从而确定第k个列索引L_k：

$\left|\tilde{E}(C_k,L_k)\right|={\max }_l\left(\right|\tilde{E}(C_k,l)\left|\right),l\≠L_1,\·\·\·,L_{k-1};$

得到Q矩阵的第ｋ行： $q_k=\tilde{E}(C_k,:)/\tilde{E}(C_k,L_k);$

更新近似误差矩阵的第L_k列： $\tilde{E}(:,L_k)=X(:,L_k)-{\mathrm{\Σ}}_{i=1}^{k-1}{\left(q_i\right)}_{L_k}{p}_i;$

得到P矩阵的第ｋ列： $p_k=\tilde{E}(:,L_k);$

${\left|\right|{\tilde{X}}^{\left(k\right)}\left|\right|}^2={\left|\right|{\tilde{X}}^{(k-1)}\left|\right|}^2+2{\mathrm{\Σ}}_{l=1}^{k-1}\left|p_l^T{p}_k\right|\·\left|q_l^T{q}_k\right|+{\left|\right|p_k\left|\right|}^2{\left|\right|q_k\left|\right|}^2;$

判断收敛误差：如果

迭代结束，否则找下一个行索引 $I_{k+1}:\left|\tilde{E}(C_{k+1},L_k)\right|={\max }_i\left(\right|\tilde{E}(c,L_k)\left|\right),c\≠C_1,\·\·\·,C_k.$ 进入第k+1次循环，直至停止条件满足。

Images