CN103217675A - 多个不共轴旋转对称体电磁散射特性的仿真方法 - Google Patents

Abstract

本发明公开了一种多个不共轴旋转对称体电磁散射特性的仿真方法，步骤如下：建立每个旋转对称体模型及相应的等效球面，分别在每个旋转对称体上建立第一、二局部直角坐标系；以平面波为激励源，确定激励向量；分别建立各个旋转对称体的第一局部坐标系下的散射矩阵；为每两个等效球面建立第三局部直角坐标系，确定该两个等效球面之间的传输矩阵；建立第一、二、三局部直角坐标系之间的旋转矩阵；建立方程组，使用迭代法得到各个等效球面上的等效电磁流；确定远区的散射场，得到雷达散射截面积。本发明对多个旋转对称体的电磁散射仿真提供了快速高效的解决方案，节省了空间资源，可以准确、快速的对多个不共轴旋转对称体的电磁散射特性进行仿真。

Description

多个不共轴旋转对称体电磁散射特性的仿真方法

技术领域

本发明涉及电磁仿真技术领域，特别是一种多个不共轴旋转对称体的电磁散射特性的仿真方法。

背景技术

旋转对称体是指绕一轴线旋转对称的三维目标，是一种常见的雷达目标，如各种多弹头问题，导弹群问题等，被广泛的关注和研究。旋转对称体矩量法可以高效的分析各类旋转对称体的电磁散射，相比于传统的基于RWG基函数剖分建模的分析方法内存消耗和计算时间都要少很多。旋转对称矩量法由Andreasen,M.G在1965年首先提出(M.Andreasen,"Scattering from bodies of revolution,"Antennas and Propagation,IEEETransactions on,vol.13,pp.303-310,1965.)。文中将入射平面波利用傅立叶级数展开为相互正交的柱面波形式，利用各模式间的正交性，分别求解单一模式下的感应电流，然后进行线性叠加求得散射场的分布。旋转对称体矩量法实际上是将原三维问题降成了两维半维问题，从而大大降低的计算复杂度。然而，对于多个旋转对称体群的散射问题，我们无法直接在一个坐标系中建立旋转对称体群中所有旋转对称体的母线且保持它们共轴，除非它们本来就公用一条旋转对称轴。因此直接使用传统的旋转对称体矩量法不能对多个不共轴的旋转对称体实现降维，因而无法发挥高效性能。

发明内容

本发明的目的在于提供一种仿真速度快、内存消耗低、精度高的多个旋转对称体电磁散射特性的仿真方法。

实现本发明的技术解决方案是：一种多个不共轴旋转对称体电磁散射特性的仿真方法，所述多个不共轴旋转对称体包括至少两个不共轴旋转对称体，步骤如下：

第1步，建立模型和局部直角坐标系：分别建立各个旋转对称体的第一局部直角坐标系，并建立各个旋转对称体的等效球面，然后建立各个等效球面的第二局部直角坐标系；选择其中一个旋转对称体的第二局部直角坐标系为全局坐标系；确定各个第二局部直角坐标系原点在全局坐标系下的坐标，以还原各个旋转对称体的实际空间位置；

第2步，以平面波为激励源，确定激励向量：在每个旋转对称体的第二局部直角坐标系下，得到等效球面上等效入射电磁流即激励向量，每个模式的激励向量乘以相位因子，还原各个等效球面的实际入射电磁流；

第3步，分别建立各个旋转对称体的第一局部坐标系下的散射矩阵：散射矩阵描述等效球面上的外界等效入射电磁流经过内部旋转对称体作用后在等效球面上产生的等效散射电磁流之间的关系；

第4步，为每两个场等效球面、源等效球面建立第三局部直角坐标系，在相应的第三局部直角坐标系下确定该场等效球面、源等效球面之间的传输矩阵；

第5步，建立第一局部直角坐标系、第二局部直角坐标系、第三局部直角坐标系之间的关系，确定各个局部直角坐标系之间的旋转矩阵；

第6步，根据第2步～第5步的信息建立求解方程组，使用迭代法得到各个等效球面上的等效电磁流；

第7步，由第6步中等效球面上的等效电磁流确定远区的散射场，得到雷达散射截面积。

本发明与现有技术相比其显著效果是：（1）仿真速度快和内存消耗低，适合电大尺寸的旋转对称体群的电磁散射问题的仿真分析；（2）建模方便，存储的信息少，分析不受各个目标相互之间距离和目标本身姿态的限制；（3）理论可靠，易于实现，对不共轴旋转对称体的电磁散射特性的仿真与设计有重要的意义。

附图说明

图1是本发明多个旋转对称体的模型结构示意图。

图2是本发明旋转对称体在第一局部直角坐标系下的母线和等效球面的母线。

图3是本发明第一局部直角坐标系和第二局部直角坐标系之间的关系。

图4本发明各个局部坐标系之间的位置关系及空间平移变量的设置。

图5本发明三角基函数示意图。

图6是本发明传输矩阵作用示意图。

图7本发明实施例1中多个旋转对称体弹头模型示意图。

图8本发明实施例1中双站雷达散射截面积。

具体实施方式

下面结合附图，以一个转对称体群的电磁散射问题的仿真分析为例（如图1所示），对实现本发明的具体步骤作进一步阐述：

第1步，建立模型和局部直角坐标系：分别建立各个旋转对称体的第一局部直角坐标系，并建立各个旋转对称体的等效球面，然后建立各个等效球面的第二局部直角坐标系；在各个第一局部直角坐标系中均建立相应旋转对称体的母线以及相应等效球面的第一母线，在各个第二局部直角坐标系中均建立相应等效球面的第二母线；选择其中一个旋转对称体的第二局部直角坐标系为全局坐标系；设置第二局部直角坐标系原点在全局坐标系下的坐标，以还原各个旋转对称体的实际空间位置；

我们按照以下方案建立旋转对称体群中的各个不同的旋转对称体的母线：

（1.1）以每个旋转对称体的轴线线段的中点为原点、轴线方向为坐标系的z轴方向，建立与各个旋转对称体一一对应的第一局部直角坐标系x'y'z'，如图1中的各个旋转对称体上建立的x'y'z'坐标系；

（1.2）以各个第一局部直角坐标系x'y'z'的原点为圆心，建立相应旋转对称体的等效球面，各个等效球面均完全包围对应的旋转对称体；以实际空间中的垂直方向为z轴方向，分别建立与各个第一局部直角坐标系共原点的第二局部直角坐标系xyz，且所有的第二局部直角坐标系的x轴、y轴、z轴分别相互平行；若旋转对称体在实际空间中轴线方向为垂直向上，那么第一局部直角坐标系和对应的第二局部直角坐标系重合，第一局部直角坐标系和第二局部直角坐标系之间的关系如图3所示；选择其中一个旋转对称体的第二局部直角坐标系为全局坐标系；

（1.3）在各x'oz'平面内建立相应旋转对称体的母线及完全包围该旋转对称体的等效球面的第一母线，如图2所示；在各xoz平面内建立相应等效球面的第二母线；

（1.4）设置第二局部直角坐标系原点在全局坐标系下的坐标，以还原各个旋转对称体的实际空间位置。如图4所示。

第2步，以平面波为激励源，确定激励向量。根据等效原理，空间中源产生的场可以由一个包围源的闭合面上的等效电磁流表示。通过建立等效球面，无需得到每个旋转对称体表面的入射场，而只需要确定等效球面上等效入射电磁流作为激励向量。平面波的表达式为

其中E为电场，H为磁场，

k为波矢量，r为观测点的位置向量，η表示自由空间中的波阻抗；

在等效球面表面建立PMCHW方程，激励向量为等效球面表面的等效入射电磁流，具体公式如下：

$\left[\begin{array}{c}{J}_H^{\mathrm{inc}}\\ M_H^{\mathrm{inc}}\end{array}\right]=\left[\begin{array}{c}-\hat{n}\×H_H^{\mathrm{inc}}\\ \hat{n}\×E_H^{\mathrm{inc}}\end{array}\right]---\left(1\right)$

下标H表示等效球面，上标inc表示入射场，

表示等效球面表面的入射电场，

表示等效球面表面的入射磁场，

表示等效球面表面的等效入射电流，

表示等效球面表面的等效入射磁流，表示等效球面的法向方向；

激励向量在每个旋转对称体的第二局部直角坐标系下求得，等效球面上的等效入射电磁流由旋转对称体基函数展开如下：

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

$\left(2\right)$

$M\left(r\right)=\underset{\mathrm{\α}=-\∞}{\overset{\∞}{\mathrm{\Σ}}}\underset{n=1}{\overset{N-1}{\mathrm{\Σ}}}[b_{\mathrm{\αn}}^t{f}_{\mathrm{\αn}}^t\left(r\right)+b_{\mathrm{\αn}}^{\mathrm{\φ}}{f}_{\mathrm{\αn}}^{\mathrm{\φ}}\left(r\right)]$

其中J(r)表示等效球面表面任意一点r点的电流，M(r)表示r点的磁流，表示电流J(r)使用旋转对称体基函数

展开对应于第α个模式的第n个基函数切向方向

的展开系数；表示电流J(r)使用旋转对称体基函数

展开对应于第α个模式的第n个基函数周向方向

的展开系数；

表示磁流M(r)使用旋转对称体基函数

展开对应于第α个模式的第n个基函数切向方向

的展开系数；

表示磁流M(r)使用旋转对称体基函数展开对应于第α个模式的第n个基函数周向方向的展开系数；α表示旋转对称体基函数的模式数编号，n表示基函数编号；N表示旋转对称体基函数对应的未知量的个数；旋转对称体基函数

表达式为：

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

$f_{\mathrm{\αn}}^{\mathrm{\φ}}\left(r\right)=\widehat{\mathrm{\φ}}\frac{T_n\left(t\right)}{\mathrm{\ρ}}{e}^{\mathrm{j\α\φ}},$ $n=1,\·\·\·,N---\left(3\right)$

其中T_n(t)表示三角基函数，是一维局部基函数，T_n(t)定义在两条相连接的剖分线段上，我们分别称这两条线段为前段和后段，其表达式为

t表示r点的切向分量；ρ表示r点到z轴的垂直距离；φ表示r点的周向角；e^jαφ表示表示傅里叶级数展开对应于第α个模式的指数项；N表示旋转对称体母线剖分段数，即旋转对称体基函数对应的未知量的个数；

表示第n个三角基函数对应的前段的起点切向分量，表示第n个三角基函数对应的前段的终点切向分量即后段的起点切向分量，

表示第n个三角基函数对应的后段的终点切向分量；Δ_n表示前段的长度，Δ_n+1表示后段的长度；三角基函数示意图如图5所示。

做伽辽金测试，生成激励向量，测试函数的表达形式为：

(5)

其中下标m表示测试函数编号，β表示测试函数模式数编号，

表示对应于切向方向

第β个模式的第m个测试函数，

表示对应于周向方向

的第β个模式第m个测试函数，N表示旋转对称体基函数对应的未知量的个数，Mod表示需要的总模式数；

对公式(1)两边使用公式(5)进行测试得到等效球面H的激励向量V_HS：

$\left[V_H\right]=\left[\begin{array}{c}{V}_H^J\\ V_H^M\end{array}\right]---\left(6\right)$

表示等效球面H表面上等效电流J对应的激励向量，

表示等效球面H表面上等效磁流M对应的激励向量，其每个元素为：

$\left[V_{\mathrm{m\β},H}^J\right]=\left[\begin{array}{c}\⟨W_{\mathrm{m\β},}^t{J}_H^{\mathrm{inc}}\left(r\right)\⟩\\ \⟨W_{\mathrm{m\β},}^{\mathrm{\φ}}{J}_H^{\mathrm{inc}}\left(r\right)\⟩\end{array}\right]$ $m=\mathrm{1,2},\·\·\·,N$

$\left(7\right)$

$\left[V_{\mathrm{m\β},H}^M\right]=\left[\begin{array}{c}\⟨W_{\mathrm{m\β},}^t{M}_H^{\mathrm{inc}}\left(r\right)\⟩\\ \⟨W_{\mathrm{m\β},}^{\mathrm{\φ}}{M}_H^{\mathrm{inc}}\left(r\right)\⟩\end{array}\right]$ $\mathrm{\β}=\mathrm{1,2},\·\·\·,\mathrm{Mod}$

其中<·>表示内积，表示等效球面H表面上等效电流J对应于第β个模式的第m个基函数的激励向量，

表示等效球面H表面上等效磁流M对应于第β个模式的第m个基函数的激励向量；

每个模式的激励向量V_H，β都乘以对应的相位因子

即：

$V_{H,\mathrm{\&beta;}}=e^{-j(k_x{x}_t+k_y{y}_t+k_z{z}_t)}{V}_{H,\mathrm{\&beta;}}$ $\mathrm{\&beta;}=1,\&CenterDot;\&CenterDot;\&CenterDot;,\mathrm{Mod}---\left(8\right)$

还原各个等效球面的实际入射电磁流，其中k_x、k_y、k_z分别为波矢量在第二局部直角坐标系的x轴、y轴、z轴坐标分量，(x_t,y_t,z_t)为第1步中所述该第二局部直角坐标系原点在全局坐标系下的坐标。

第3步，分别建立各个旋转对称体的第一局部坐标系下的散射矩阵。散射矩阵描述等效球面上的外界等效入射电磁流经过内部旋转对称体作用后在等效球面上产生的等效散射电磁流之间的关系，最终使得等效球面可以完全描述其内部包含的旋转对称体的散射作用。

每个散射矩阵S包括三部分：等效球面上的模式β对应的等效入射电磁流对旋转对旋转对称体的作用描述为Z_β ^PH；旋转对称体自身的散射作用描述为

该阻抗矩阵的建立和传统的旋转对称体矩量法求解单个旋转对称体方法相同；由旋转对称体表面自身产生的模式β对应的散射场在等效球面上得到的模式β对应的等效散射电磁流的作用描述为Z_β ^HP，则S矩阵可以写成：

$S_{\mathrm{\&beta;}}=U^{-1}{Z}_{\mathrm{\&beta;}}^{\mathrm{HP}}{\left(Z_{\mathrm{\&beta;}}^{\mathrm{pp}}\right)}^{-1}{Z}_{\mathrm{\&beta;}}^{\mathrm{PH}}---\left(9\right)$

${Z_{\mathrm{\&beta;}}}^{\mathrm{PH}}=\left[\begin{array}{cc}{L}_{\mathrm{PH}}\left(J_{\mathrm{\&beta;},H}\right)& {\mathrm{\&eta;K}}_{\mathrm{PH}}\left(M_{\mathrm{\&beta;},H}\right)\end{array}\right]$

${Z_{\mathrm{\&beta;}}}^{\mathrm{HP}}=\left[\begin{array}{c}\hat{n}\&times;K_{\mathrm{HP}}\left(J_{\mathrm{\&beta;},P}\right)\\ \hat{n}\&times;L_{\mathrm{HP}}\left(J_{\mathrm{\&beta;},P}\right)\end{array}\right]---\left(10\right)$

其中U^-1为电流系数求解矩阵：

$\begin{array}{cc}& \left[U_{\mathrm{mn}}\right]=[\&lang;W_{\mathrm{\&beta;m}},f_{\mathrm{\&alpha;n}}\&rang;]\begin{array}{c}n=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;,N\\ m=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,N\end{array}---\left(11\right)\end{array}$

其中f_αn表示旋转对称体基函数，n表示基函数编号，α表示模式数，W_βm表示旋转对称体测试函数，m表示测试函数的编号，β表示测试函数模式数，代入具体得到：

$\left[U_{\mathrm{mn}}\right]=\&lang;w_{\mathrm{\&beta;m}},f_{\mathrm{\&alpha;n}}\&rang;={\&Integral;}_{t_m}{\&Integral;}_0^{2\mathrm{\&pi;}}\mathrm{\&rho;}{f}_m\left(t\right)f_n\left(t\right)e^{\mathrm{j\&phi;}(\mathrm{\&alpha;}-\mathrm{\&beta;})}\mathrm{d\&phi;dt}$

$=\underset{p=1}{\overset{M_p}{\mathrm{\&Sigma;}}}\frac{\mathrm{\&pi;}}{4{\mathrm{\&rho;}}_p}{T}_{\mathrm{mp}}{T}_{\mathrm{np}}{\mathrm{\&Delta;}}_p---\left(12\right)$

其中T_mp、T_np表示对应于m、n第p条线段的三角基函数，M_p表示一个基函数包含的剖分线段数，Δ_p表示对应的剖分线段的长度；L、K表示积分算子，如下：

$L\left(\mathrm{\&chi;}\left(r^{\&prime;}\right)\right)=\mathrm{j\&omega;\&mu;}\&Integral;\underset{s}{\&Integral;}\mathrm{\&chi;}\left(r^{\&prime;}\right)G(r,r^{\&prime;})+\frac{1}{{\mathrm{\&omega;}}^2\mathrm{\&mu;\&epsiv;}}\&dtri;[{\&dtri;}^{\&prime;}\&CenterDot;\mathrm{\&chi;}\left(r^{\&prime;}\right)G(r,r^{\&prime;})]{\mathrm{ds}}^{\&prime;}$

$\left(13\right)$

$K\left(\mathrm{\&chi;}\left(r^{\&prime;}\right)\right)=\&Integral;\underset{s}{\&Integral;}\mathrm{\&chi;}(r\&prime;)\&times;\&dtri;G(r,r^{\&prime;}){\mathrm{ds}}^{\&prime;}$

χ(r')表示等效电流J(r')或等效磁流M(r')，它们将使用公式（2）离散，与各个模式电磁流分量相对应，ω表示角频率，ε表示电导率，μ表示磁导率；G(r,r')表示源点r'到场点r的格林函数其表达式为：

$G(r,r^{\&prime;})=\frac{e^{-\mathrm{jk}|r-r^{\&prime;}|}}{|r-r^{\&prime;}|}---\left(14\right)$

下标H表示等效球面，P表示旋转对称体，即积分算子L_PH(J_β,H)表示等效球面H外部等效面上的第β模式的电流J_β,H在旋转对称体P表面产生的β模式的激励电场，积分算子K_PH(M_β,H)表示等效球面H外部等效面上的第β模式的磁流M_β,H在旋转对称体P表面产生的β模式的激励磁场；积分算子K_HP(J_β,P)表示旋转对称体P上的β模式的电流J_β,P在等效球面H外部等效面上产生的β模式的激励磁场，积分算子L_HP(J_β,P)表示旋转对称体P上的β模式的电流J_β,P在等效球面H外部等效面上产生的β模式的激励电场；

对于旋转对称体基函数而言，每个模式数对应的系数求解矩阵都是相同的，由于三角基函数只和相邻的两个基函数有重叠，因此矩阵U是一个三带条型矩阵，可以通过追赶法求解得到电流系数求解矩阵U^-1。同时等效球面和旋转对称体的母线是共轴的，满足旋转对称体矩量法降维条件，因此散射矩阵是由Mod个独立的小矩阵组成的。

第4步，确定每两个场等效球面、源等效球面之间的传输矩阵。通过第3步散射矩阵的求解，已经将旋转对称体的散射作用转换到了等效球面上，问题转换成了空间中任意放置的均不相交的球面之间耦合作用的求解问题。使用传输矩阵描述一个源等效球面上的等效电磁流对一个场等效球面的作用。

(4.1)为了能够使用旋转对称体矩量法确定每两个场等效球面、源等效球面之间的传输矩阵，建立第三局部直角坐标系，第三局部直角坐标系的建立标准是以场等效球面的球心为坐标原点，源等效球面球心和场等效球面球心之间的连线为z''轴，z''轴方向指向源等效球面。x''oy''平面按照图3的旋转法则确定，该旋转法则在第五步中将给出。在x''oz''平面内建立场等效球面和源等效球面的母线，如图6所示。

(4.2)传输矩阵T^HH描述的是源等效球面上的等效电磁流在场等效球面上产生的二次散射电磁流之间的关系。由于对于每个源等效球面而言，等效球面之间传输矩阵的决定因素包括球的半径和两者之间的距离，若二者皆相同，那么传输矩阵也是相同的，仅仅需要构建一次。传输矩阵T^HH的公式如下：

$T^{\mathrm{HH}}=\left[\begin{array}{cc}-\hat{n}\&times;K_{\mathrm{HH}}\left(J\right)& -\frac{1}{\mathrm{\&eta;}}\hat{n}\&times;L_{\mathrm{HH}}\left(M\right)\\ -\frac{1}{\mathrm{\&eta;}}\widehat{\mathrm{\&eta;}}\&times;L_{\mathrm{HH}}\left(J\right)& -\hat{n}\&times;K_{\mathrm{HH}}\left(M\right)\end{array}\right]---\left(15\right)$

场等效球面的法线方向，η表示自由空间中的波阻抗，积分算子K_HH(J)表示源等效球面上的等效电流J在场等效球面上的感应磁场，积分算子L_HH(J)表示源等效球面上的等效电流J在场等效球面上的感应电场，积分算子L_HH(M)表示源等效球面上的等效磁流M在场等效球面上的感应电场，K_HH(M)表示源等效球面上的等效磁流M在场等效球面上的感应磁场；

每个传输矩阵描述的都是在各自对应的第三局部直角坐标系下源等效球面对场等效球面的作用，而每两个场等效球面、源等效球面所建立的第三局部直角坐标系都是不相同的，因此产生二次散射电磁流也是在不同的第三局部直角坐标系下表示出来的。

第5步，建立第一局部直角坐标系、第二局部直角坐标系、第三局部直角坐标系之间的关系。依次确定第一局部直角坐标系和第二局部直角坐标系之间的旋转矩阵，局部第二局部直角坐标系和第三局部直角坐标系之间的旋转矩阵。

以一个等效球面为源等效球面，建立以该等效球面球心为原点的3个局部直角坐标系。为了限定3个局部直角坐标系之间的关系，按照下面的旋转规则描述3个局部直角坐标系之间的关系：

第一局部直角坐标系和第三局部直角坐标系都是由第二局部直角坐标系旋转而来，如图3所示，x'y'z'为第一局部直角坐标系，xyz为第二局部直角坐标系。

（5.1）首先建立第二局部直角坐标系，确定由第二局部直角坐标系到第一局部直角坐标系的旋转角度，在第一局部直角坐标系的z'轴上任取一点G为旋转参考点，G点在第二局部直角坐标系中的坐标为(x₀,y₀,z₀)，则：

${{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}=\frac{z_0}{\sqrt{{x_0}^2+{y_0}^2+{z_0}^2}}---\left(16\right)$

${{\tan \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}=\frac{y_0}{x_0}$

θ₁₂ ^rot和φ₁₂ ^rot为旋转参考角度，下标12表示由第二局部直角坐标系旋转得到第一局部直角坐标系，如图3所示，根据θ₁₂ ^rot和φ₁₂ ^rot，按照以下两步进行由第二局部直角坐标系旋转得到第一局部直角坐标系：

①以第二局部直角坐标系的z轴为旋转轴，顺时针旋转

使得第二局部直角坐标系的y轴经过G点在xoy平面上的投影点G'，得到第一局部直角坐标系x'轴的位置；

②以①中得到的x'轴为旋转轴顺时针旋转θ₁₂ ^rot，使得第二局部直角坐标系的z轴经过G点，得到第一局部直角坐标系的z'轴。

由以上两步得到了第一局部直角坐标系，设A点为第二局部直角坐标系中的一点，坐标为(x₁,y₁,z₁)，则A点在第一局部直角坐标系下对应的坐标(x₁',y₁',z₁')为：

$\left[\begin{array}{c}{x_1}^{\&prime;}\\ {y_1}^{\&prime;}\\ {z_1}^{\&prime;}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& {{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}& {{-\sin \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\\ 0& {{\sin \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\end{array}\right]\left[\begin{array}{ccc}{{\sin \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{-\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& 0\\ {{\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{\sin \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}\right]---\left(17\right)$

即第二局部直角坐标系到第一局部直角坐标系的旋转矩阵R¹²为：

$R^{12}=\left[\begin{array}{ccc}{{\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{-\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& 0\\ {{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}{{\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}{{\sin \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{-\sin \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\\ {{\sin \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}{{\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{\sin \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}{{\sin \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\end{array}\right]---\left(18\right)$

（5.2）还需要由第一局部直角坐标系旋转得到第二局部直角坐标系：

第二局部直角坐标系是一个由全局坐标系平移得到的局部直角坐标系，旋转第一局部直角坐标系满足第二局部直角坐标系的规定即可，旋转过程也可以描述为以下两步：

①以第一局部直角坐标系的x'轴为旋转轴，逆时针旋转θ₁₂ ^rot，使得第一局部直角坐标系的y'轴落在第二局部直角坐标系的xoy平面上，得到第二局部直角坐标系z轴的位置；

②以①中得到的z轴为旋转轴逆时针旋转使得第一局部直角坐标系的x'轴和第二局部直角坐标系规定的x轴重合，得到第二局部直角坐标系。

由第一局部直角坐标系到第二局部直角坐标系的矩阵R²¹为：

$R^{21}=\left[\begin{array}{ccc}{{\sin \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}{{\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{\sin \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}{{\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}\\ {{-\cos \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}{{\sin \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& {{\sin \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}{{\sin \mathrm{\&phi;}}_{12}}^{\mathrm{rot}}\\ 0& {{-\sin \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\end{array}\right]---\left(19\right)$

（5.3）由第二局部直角坐标系旋转到第三局部直角坐标系的旋转过程和由第二局部直角坐标系旋转到第一局部直角坐标系的过程一致，但旋转参考点变成了源等效球面的球心O_s，且要考虑场等效球面的第二局部直角坐标系在全局坐标系下的平移分量，设源等效球面球心O_s在全局坐标系下的坐标为(x_s,y_s,z_s)，场等效球面球心O_o在全局坐标系下的坐标为(x_o,y_o,z_o)，则旋转参考角度为：

${{\tan \mathrm{\&phi;}}_{23}}^{\mathrm{rot}}=\frac{y_o-y_s}{x_o-x_s}$

和

为旋转参考角度，下标23表示由第三局部直角坐标系旋转得到第二局部直角坐标系，下标32表示由第二局部直角坐标系旋转得到第三局部直角坐标系；

则第二局部直角坐标系到第三局部直角坐标系的旋转矩阵R³²为：

$R^{32}=\left[\begin{array}{ccc}{{\cos \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{-\cos \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& 0\\ {{\cos \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}{{\cos \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}{{\sin \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{-\sin \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\\ {{\sin \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}{{\cos \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{\sin \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}{{\sin \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\end{array}\right]---\left(21\right)$

第三局部直角坐标系到第二局部直角坐标系的旋转矩阵R²³为

$R^{23}=\left[\begin{array}{ccc}{{\sin \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}{{\cos \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{\sin \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}{{\cos \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}\\ {{-\cos \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}{{\sin \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& {{\sin \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}{{\sin \mathrm{\&phi;}}_{32}}^{\mathrm{rot}}\\ 0& {{-\sin \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}& {{\cos \mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\end{array}\right]---\left(22\right)$

第6步，根据第2步～第5步的信息建立求解方程组，使用迭代求解法求解各个等效球面上的等效电磁流，迭代求解精度设为1.0e-3。

由电磁场边界条件和零场等效原理在每个等效球面上建立方程：

其中下标H表示等效球面，1表示第一局部直角坐标系，2表示第二局部直角坐标系，3表示第三局部直角坐标系，ξ表示激励向量对应在第二局部直角坐标系下模式数编号，ζ表示由第二局部坐标系旋转得到的第一局部直角坐标系下对应的模式数编号，γ表示由第一局部坐标系旋转得到的第二局部坐标系下对应的模式数编号，ν表示由第二局部直角坐标系旋转得到的第三局部直角坐标系下对应的模式数编号，υ表示由第三局部直角坐标系旋转得到的第二局部直角坐标系下对应的模式数编号，ζ、γ、ν、υ取值范围均为-Mod,…,Mod，Mod为需要的总模式数，上标inc表示入射场且为已知量，sca表示散射场且为未知量，i表示场等效球面编号，j表示源等效球面编号，i＝1,2,…,N_obj，j＝1,2,…,N_obj，N_obj表示旋转对称体群中旋转对称体总数，S为散射矩阵，T为传输矩阵，R为旋转矩阵，J为电流，M为磁流。

使用GMRES迭代求解上述方程（23）。解该方程仅需要进行一系列小规模的独立模式矩阵矢量乘，仿真速度快，一般设置截断精度为1.0e-3。

第7步，由第6步中等效球面上的等效电磁流确定远区的散射场，最终求解雷达散射截面积。

每个等效球面上的等效散射电磁流系数已确定，和第2步类似，在第二局部直角坐标系下面的要给区域δ的远场E_δ ^sca需要乘以一个平移因子

${E_{\mathrm{\&delta;}}}^{\mathrm{sca}}={E_{\mathrm{\&delta;}}}^{\mathrm{sca}}\&CenterDot;e^{-j(k_x{x}_t+k_y{y}_t+k_z{z}_t)}---\left(24\right)$

由下列公式确定雷达散射截面积σ：

$\mathrm{\&sigma;}=\underset{r\&RightArrow;\&infin;}{\lim }{4\mathrm{\&pi;r}}^2\frac{{\left|E^s\right|}^2}{{\left|E^i\right|}^2}---\left(25\right)$

其中E^s为远区散射场，Eⁱ表示远区入射场。

实施例1

根据本发明所述的多个不共轴旋转对称体的电磁散射特性的仿真方法，我们对一个实际的多弹头的电磁散射进行仿真。

多个旋转对称体弹头如图7所示，以1号弹头的第二局部直角坐标系为全局坐标系，各个弹头以2m沿x，y两个方向为周期均匀分布。弹头的表面为抛物面，抛物面口径为0.500m，深度为1.5625m，焦距长度为0.01m。平面波入射角度为（0°，0°），工作频率为300MHz。我们使用本发明方法和多层快速多极子的仿真结果进行比较，如图8所示，结果相互吻合的很好，证明了本发明方法的正确性。多层快速多极子使用的是RWG基函数，未知量为113778，本发明使用线剖分，剖分线段数为37，两种方法皆按照0.05波长剖分（为了很少的逼近原模型剖分密度略大于0.1波长），本例子本发明方法中用到了21个模式数。多层快速多极子是现有分析仿真散射问题效率很高的一种方法，表1给出了本发明方法和多层快速多极子方法效率的对比，内存消耗约为多层快速多极子的1/30，矩阵填充时间约为多层快速多极子方法的1/50。迭代收敛步数约为多层快速多极子的1/3。

表1

方法	内存消耗（兆）	矩阵方程形成时间（秒）	迭代求解时间	迭代步数
					本发明方法	51.1	6.2	92.5	12
多层快速多极子	1224.2	338.5	378.5	45

综上所述，本发明提供一种多个不共轴旋转对称体的电磁散射特性的仿真方法，该方法使得传统的旋转对称体矩量法不仅仅可以高效的分析单一的旋转对称体目标，也可以高效的分析多个不共轴旋转对称体电磁散射问题，内存消耗和仿真时间都相比于现有基于三维剖分基函数矩量法方法低，仿真资源得到了很大地节省，可以准确、快速的仿真多个不共轴的旋转对称体的电磁散射。

Claims (5)

1.一种多个不共轴旋转对称体电磁散射特性的仿真方法，所述多个不共轴旋转对称体包括至少两个不共轴旋转对称体，其特征在于，步骤如下：

第1步，建立模型和局部直角坐标系：分别建立各个旋转对称体的第一局部直角坐标系，并建立各个旋转对称体的等效球面，然后建立各个等效球面的第二局部直角坐标系；选择其中一个旋转对称体的第二局部直角坐标系为全局坐标系；确定各个第二局部直角坐标系原点在全局坐标系下的坐标，以还原各个旋转对称体的实际空间位置；

第2步，以平面波为激励源，确定激励向量：在每个旋转对称体的第二局部直角坐标系下，得到等效球面上等效入射电磁流即激励向量，每个模式的激励向量乘以相位因子，还原各个等效球面的实际入射电磁流；

第3步，分别建立各个旋转对称体的第一局部坐标系下的散射矩阵：散射矩阵描述等效球面上的外界等效入射电磁流经过内部旋转对称体作用后在等效球面上产生的等效散射电磁流之间的关系；

第4步，为每两个场等效球面、源等效球面建立第三局部直角坐标系，在相应的第三局部直角坐标系下确定该场等效球面、源等效球面之间的传输矩阵；

第5步，建立第一局部直角坐标系、第二局部直角坐标系、第三局部直角坐标系之间的关系，确定各个局部直角坐标系之间的旋转矩阵；

第6步，根据第2步～第5步的信息建立求解方程组，使用迭代法得到各个等效球面上的等效电磁流；

第7步，由第6步中等效球面上的等效电磁流确定远区的散射场，得到雷达散射截面积。

2.根据权利要求1所述的多个不共轴旋转对称体电磁散射特性的仿真方法，其特征在于，第1步中所述建立模型和局部直角坐标系的具体过程如下：

（1.1）以每个旋转对称体的轴线线段的中点为原点、轴线方向为坐标系的z'轴方向，建立与各个旋转对称体一一对应的第一局部直角坐标系x'y'z'；

（1.2）以第一局部直角坐标系x'y'z'的原点为圆心，建立相应旋转对称体的等效球面，各个等效球面均完全包围对应的旋转对称体；以实际空间中的垂直方向为z轴方向，分别建立与各第一局部直角坐标系共原点的第二局部直角坐标系xyz，且所有的第二局部直角坐标系的x轴、y轴、z轴分别相互平行；选择其中一个旋转对称体的第二局部直角坐标系为全局坐标系；

（1.3）在每个第一局部直角坐标系的x'oz'平面内建立相应旋转对称体的母线及完全包围该旋转对称体的等效球面的第一母线，在每个第二坐标系的xoz平面均建立相应等效球面的第二母线；

（1.4）设置第二局部直角坐标系原点在全局坐标系下的坐标，以还原各个旋转对称体的实际空间位置。

3.根据权利要求1所述的多个不共轴旋转对称体电磁散射特性的仿真方法，其特征在于，第2步中所述的相位因子为其中k_x、k_y、k_z分别为平面波矢量在第二局部直角坐标系的x轴、y轴、z轴坐标分量，(x_t,y_t,z_t)为第二局部直角坐标系原点在全局坐标系下的坐标。

4.根据权利要求1所述的多个不共轴旋转对称体电磁散射特性的仿真方法，其特征在于，第4步中所述传输矩阵的具体构建过程如下：

（4.1）为每两个场、源等效球面建立第三局部直角坐标系：以场等效球面的球心为坐标原点，源等效球面球心和场等效球面球心之间的连线为z''轴，z''轴方向指向源等效球面；在x''oz''平面内建立场等效球面和源等效球面的母线；

（4.2）在每两个等效球面对应的第三局部直角坐标系下确定该两个等效球面之间的传输矩阵T^HH，公式如下：

$T^{\mathrm{HH}}=\left[\begin{array}{cc}-\hat{n}\&times;K_{\mathrm{HH}}\left(J\right)& -\frac{1}{\mathrm{\&eta;}}\hat{n}\&times;L_{\mathrm{HH}}\left(M\right)\\ -\frac{1}{\mathrm{\&eta;}}\hat{n}\&times;L_{\mathrm{HH}}\left(J\right)& -\hat{n}\&times;K_{\mathrm{HH}}\left(M\right)\end{array}\right]$

表示场等效球面的法线方向，η表示自由空间中的波阻抗，积分算子K_HH(J)表示源等效球面上的等效电流J在场等效球面上的感应磁场，积分算子L_HH(J)表示源等效球面上的等效电流J在场等效球面上的感应电场，积分算子L_HH(M)表示源等效球面上的等效磁流M在场等效球面上的感应电场，K_HH(M)表示源等效球面上的等效磁流M在场等效球面上的感应磁场。

5.根据权利要求1所述的多个不共轴旋转对称体电磁散射特性的仿真方法，其特征在于，第5步所述建立第一局部直角坐标系、第二局部直角坐标系、第三局部直角坐标系之间的关系，具体如下：

（5.1）确定由第二局部直角坐标系到第一局部直角坐标系的旋转矩阵

在第一局部直角坐标系的z'轴上任取一点G为旋转参考点，G点在第二局部直角坐标系中的坐标为(x₀,y₀,z₀)，则：

$\cos {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}=\frac{Z_0}{\sqrt{{x_0}^2+{y_0}^2+{z_0}^2}}$

$\tan {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}=\frac{y_0}{x_0}$

θ₁₂ ^rot和φ₁₂ ^rot为旋转参考角度，下标12表示由第二局部直角坐标系旋转得到第一局部直角坐标系；根据θ₁₂ ^rot和φ₁₂ ^rot，按照以下两步进行由第二局部直角坐标系旋转得到第一局部直角坐标系：

①以第二局部直角坐标系的z轴为旋转轴，顺时针旋转使得第二局部直角坐标系的y轴经过G点在xoy平面上的投影点G'，得到第一局部直角坐标系x'轴的位置；

②以①中得到的x'轴为旋转轴顺时针旋转θ₁₂ ^rot，使得第二局部直角坐标系的z轴经过G点，得到第一局部直角坐标系的z'轴；

第二局部直角坐标系到第一局部直角坐标系的旋转矩阵R¹²为：

$R^{12}=\left[\begin{array}{ccc}\cos {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& -\cos {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& 0\\ \cos {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\cos {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\sin {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& -\sin {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\\ \sin {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\cos {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& \sin {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\sin {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\end{array}\right]$

（5.2）确定由第一局部直角坐标系到第二局部直角坐标系的旋转矩阵

按照以下两步进行由第一局部直角坐标系旋转得到第二局部直角坐标系：

①以第一局部直角坐标系的x'轴为旋转轴，逆时针旋转θ₁₂ ^rot，使得第一局部直角坐标系的y'轴落在第二局部直角坐标系的xoy平面上，得到第二局部直角坐标系z轴的位置；

②以①中得到的z轴为旋转轴逆时针旋转

使得第一局部直角坐标系的x'轴和第二局部直角坐标系规定的x轴重合，得到第二局部直角坐标系；

第一局部直角坐标系到第二局部直角坐标系的旋转矩阵R²¹为：

$R^{21}=\left[\begin{array}{ccc}\sin {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\cos {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& \sin {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\cos {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}\\ -\cos {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\sin {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}& \sin {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\sin {{\mathrm{\&phi;}}_{12}}^{\mathrm{rot}}\\ 0& -\sin {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{12}}^{\mathrm{rot}}\end{array}\right]$

（5.3）确定第二局部直角坐标系和第三局部直角坐标系之间的旋转矩阵

设源等效球面球心O_s在全局坐标系下的坐标为(x_s,y_s,z_s)，场等效球面球心O_o在全局坐标系下的坐标为(x_o,y_o,z_o)，则旋转参考角度为：

${{\cos \mathrm{\&theta;}}_{23}}^{\mathrm{rot}}=\frac{z_o-z_s}{\sqrt{{(x_o-x_s)}^2+{(y_o-y_s)}^2+{(z_o-z_s)}^2}}$

$\tan {{\mathrm{\&phi;}}_{23}}^{\mathrm{rot}}=\frac{y_o-y_s}{x_o-x_s}$

则第二局部直角坐标系到第三局部直角坐标系的旋转矩阵R³²为：

$R^{32}=\left[\begin{array}{ccc}\cos {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& -\cos {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& 0\\ \cos {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\cos {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\sin {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& -\sin {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\\ \sin {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\cos {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& \sin {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\sin {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\end{array}\right]$

第三局部直角坐标系到第二局部直角坐标系的旋转矩阵R²³为：

$R^{23}\left[\begin{array}{ccc}\sin {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\cos {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& \sin {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\cos {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}\\ -\cos {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\sin {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}& \sin {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\sin {{\mathrm{\&phi;}}_{32}}^{\mathrm{rot}}\\ 0& -\sin {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}& \cos {{\mathrm{\&theta;}}_{32}}^{\mathrm{rot}}\end{array}\right]$

和

为旋转参考角度，下标23表示由第三局部直角坐标系旋转得到第二局部直角坐标系，下标32表示由第二局部直角坐标系旋转得到第三局部直角坐标系。

Images