一种计算高速微动类球头柱体目标的电磁散射特性方法与流程

本发明属于电磁场数值计算领域，具体涉及一种计算高速微动的类球头柱体金属目标电磁散射特性的方法

背景技术：

随着航空器件的飞行速度越来越快，关于高速运动物体的电磁研究已经成为热点问题，研究高速目标的微动情况也逐渐成为必要的关注点。例如，微动情况包括旋转，章动和进动等。

以往的方法因坐标系限制、边界条件或者求解结果单一等缺点，不能适用于高速乃至超高速运动物体的电磁场值研究，也更不能适用于高速运动目标的微动研究。因此，非常有必要设计一种能够解决并分析高速微动目标电磁散射的技术。

对于运动目标电磁散射计算的数值方法，部分学者提出了等效面电流法，运动边界法等。这些方法因坐标系限制、边界条件或者求解结果单一，或者不能适应于三维目标复杂多样的形状和目标微动的散射场的求解等缺点，不能适用于高速乃至超高速运动物体的电磁场值研究；对于用矩量法和有限元法来计算高速微动目标的电磁散射，不能适应于电大尺寸目标，且耗时较长。因此，非常有必要设计一种能够解决并分析高速微动目标电磁散射的技术。

技术实现要素：

本发明解决的技术问题是：为解决现有技术中的高速目标微动状态下的电磁散射问题，本发明设计一种计算高速微动的类球头柱体金属目标电磁散射特性的方法。

本发明的技术方案是：一种计算高速微动类球头柱体目标的电磁散射特性方法，包括以下步骤：

步骤一：建立类球头柱体高速微动目标网格模型；

步骤二：设目标的速度为v，初始位置为原点；目标在运动时，设目标的滚转角为ψ，俯仰角为θ，偏航角为目标的章动频率为f1,进动频率为f2，旋转频率为f3；定义类球头柱体轴线方向为x轴方向，径向方向为y轴方向，与x轴、y轴均垂直的方向为z轴方向；俯仰角为θ目标坐标系x轴与水平面的夹角，偏航角为目标坐标系x轴在水平面上投影与地面坐标系xg轴之间的夹角，滚转角ψ为目标坐标系z轴与通过目标x轴的铅垂面间的夹角，目标向右滚为正，反之为负；

步骤三：将模型中的目标转换为相对静止的目标进行计算，依次建立四个坐标系，包括以下子步骤：

子步骤一：将目标水平放置，以目标的质心为原点，垂直方向为z轴，沿半球方向为-y轴，建立静止坐标系s；

子步骤二：将目标的速度v代入坐标系s中，得到新的坐标系p，且s坐标系到p坐标系的转换矩阵为

其中φv为运动方向相对于速度方向在xsys(s坐标系下的x轴y轴组成的平面)平面的投影与xs轴(s坐标系的x轴)之间的方位角，为运动方向相对于速度方向的仰角(以速度方向建立坐标系)；

子步骤三：将子步骤二中的坐标系p进行洛伦兹变换，得到相对静止的坐标系m，洛伦兹转换矩阵为设为lt；

上式(2)中vx,vy,vz分别为目标沿x轴,y轴,z轴的速度的分量，c为光速，β＝v/c；

子步骤四：将滚转角ψ，俯仰角θ，和偏航角代入子步骤三中的坐标系m中，得到与目标保持相对静止的运动坐标系s',由坐标系m到s'的转换矩阵为

由于章动,进动和旋转的影响，滚转角ψ，俯仰角θ，和偏航角是随时间变化的，变化规律为

其中ψ0为初始滚转角，θp为进动角，θm为最大章动角，为初始偏航角，ω1为章动角频率，ω2为进动角频率，ω3为旋转角频率，其中ω1＝2πf1，ω2＝2πf2，ω3＝2πf3。

子步骤五：矩阵rmic、矩阵lt和矩阵rini依次相乘得到矩阵b，b＝rini×lt×rmic，设b矩阵的逆矩阵为矩阵f；

步骤四，定义静止坐标系s下的入射波幅值e0，频率f以及入射波照射到目标上的入射角θis,经依次转换后的入射波的幅值为e0′，频率为fs′，入射角度为θi′s,

步骤五：运用fdtd算法,计算运动坐标系s'中电场与磁场

步骤六：将步骤五得到的电场与磁场通过矩阵rmic转换到坐标系m中，得到新的电场与磁场并将得到的电场与磁场外推对齐坐标得到新的电场与磁场

步骤七：将步骤六得到运动坐标系s'的电场与磁场经过两次转换(m→p→s)到s坐标系，即得到类球头柱体高速微动的电场与磁场

步骤八：将得到类球头柱体高速微动的电磁散射场进行能量计算和傅里叶变换、滤波的处理，得到其频率分布。

发明效果

本发明的技术效果在于：本发明通过visualstudio2015和fortran2016语言进行网格建模，计算与目标的平动相对静止的坐标系的入射波，加入射波，用fdtd计算s′坐标系的场值，并转换到表示目标自身运动的坐标系m下进行外推处理，最后转换到s坐标系下得到电磁散射场。本发明方法新颖，思路清晰，理论饱满，在任何速度下的微动目标的特征判断及目标识别方面具有理论指导作用和广阔的应用前景。本发明采用fdtd求微动目标散射特性，并继承了fdtd的计算简单，编程容易，可以计算电大尺寸的优点。和现有的方法相比，本发明计算速度快，可计算高速微动的电大尺寸目标，并适合任何速度下的微动目标的特征分析。

附图说明

图1：一种分析高速微动的类球头柱体金属目标电磁散射特性的计算方法的目标模型图；

图2：一种分析高速微动的类球头柱体金属目标电磁散射特性的计算方法的建模三视图；

图3：一种分析高速微动的类球头柱体金属目标电磁散射特性的计算方法的四个坐标系的关系图；

图4：一种分析高速微动的类球头柱体金属目标电磁散射特性的计算方法的边界示意图；

图5：一种分析高速微动的类球头柱体金属目标电磁散射特性的计算方法的流程图的示意图；

图6：一种分析高速微动的类球头柱体金属目标电磁散射特性的计算方法的外推前的网格对齐示意图；

图7：一种分析高速微动的类球头柱体金属目标电磁散射特性的计算方法的散射场的能量图；

图8：一种分析高速微动的类球头柱体金属目标电磁散射特性的计算方法的频谱图；

具体实施方式

参见图1—图8，一种计算类球头柱体高速微动目标电磁散射场的方法，包括以下步骤：

步骤一，建立类球头柱体高速微动目标网格模型；

步骤二，设目标的速度为v，初始位置为原点。目标在运动时，设目标的滚转角为ψ，俯仰角为θ，偏航角为目标的章动频率为f1,进动频率为f2，和旋转频率f3；

步骤三，建立四个坐标系。静止坐标系s，建立表示初始位置和方位的坐标系p，建立表示目标自身运动的坐标系m，建立与目标相对静止的运动坐标系s'。从坐标系s到坐标系p的转换矩阵设为目标方位矩阵rini。经过洛伦兹变换，将坐标系p转换到坐标系m，设转换矩阵设为矩阵lt。从坐标系m到坐标系s'的转换矩阵设为rmic。将矩阵rini，矩阵lt，矩阵rmic依次相乘得到从静止坐标系到与目标平动时相对静止的坐标系的转换矩阵，设为矩阵b；

步骤四，将入射波通过转换矩阵经过三次坐标系转换(s→p→m→s')到s'坐标系；

步骤五，运用fdtd算法,计算得到运动坐标系s'中电场与磁场数据；

步骤六，将步骤五得到的电场与磁场通过矩阵rmic转换到坐标系m中，得到新的电场与磁场并将得到的电场与磁场外推对齐坐标得到新的电场与磁场数据；

步骤七，将步骤六得到运动坐标系s'的电场与磁场经过两次转换(m→p→s)到s坐标系，即得到类球头柱体高速微动的电场与磁场数据；

步骤八，将得到类球头柱体高速微动的电磁散射场进行能量计算，傅里叶变换和滤波的处理；

其中所述步骤一具体包括以下过程：

过程2.1，使用软件visualstudio2015和fortran2016语言建模，依据类球头柱体的实际长度，将模型离散化成一个个小的同等大小的长方体，长方体的长宽高设为δx',δy',δz',根据fdtd算法，δx',δy',δz'至少小于实际长度的十分之一，此次设δx'＝δy'＝δz'。确定整个fdtd计算区域和三维空间的总场区大小和散射场区大小，并确定散射边界以及吸收边界，总场边界；

过程2.2，由过程2.1确定的整个fdtd计算区域采用取整法，四舍五入，离散为一个个长宽高为δx',δy',δz'的长方体。

过程2.3，由过程2.2得到离散后fdtd计算区域以目标的质心所在的长方体为设为(0,0,0)，若空间第(i,j,k)个长方体在目标上，则设置为金属介质，否则设置为空气介质。依次类推，即可得到空间内的类球头柱体高速微动目标模型，并已将其按网格剖分为长方体元胞模型，得到.txt格式的平面模型数据；

所述步骤二具体包括以下过程：

过程3.1设目标的速度为v，初始位置为原点。目标在运动时，设目标的滚转角为ψ，俯仰角为θ，偏航角为目标的章动频率为f1,进动频率为f2，旋转。其中俯仰角为θ目标坐标系x轴与水平面的夹角，偏航角为目标坐标系x轴在水平面上投影与地面坐标系xg轴之间的夹角，滚转角ψ为目标坐标系z轴与通过目标x轴的铅垂面间的夹角，目标向右滚为正，反之为负，目标坐标系见图1。

所述步骤三具体包括以下过程：

过程4.1，将目标水平放置，以目标的质心为原点，垂直方向为z轴，沿半球方向为-y轴，建立坐标系s，可见图1。

过程4.2，将目标的速度v引入到由过程4.1得到坐标系s中，得到新的坐标系p，且s坐标系到p坐标系的转换矩阵为

式(1)中φv为运动方向相对于速度方向在xsys平面的投影与xs轴(s坐标系的x轴)之间的方位角，为运动方向相对于速度方向的仰角(以速度方向建立坐标系)；

过程4.3,将过程3.2得到坐标系p进行洛伦兹变换得到相对静止的坐标系m，洛伦兹转换矩阵为设为lt。

上式(2)中vx,vy,vz分别为目标沿x轴,y轴,z轴的速度的分量，c为光速，β＝v/c；

过程4.4，滚转角ψ，俯仰角θ，和偏航角引入到由过程4.4得到坐标系m中，得到与目标保持相对静止的运动坐标系s',由由坐标系m到s'的转换矩阵为

由于章动,进动和旋转的影响，滚转角ψ，俯仰角θ，和偏航角随时间变化为

其中ψ0为初始滚转角，θp为进动角，θm为最大章动角，为初始偏航角，ω1为章动角频率，ω2为进动角频率，ω3为旋转角频率，其中ω1＝2πf1，ω2＝2πf2，ω3＝2πf3。

过程4.5，矩阵rmic,矩阵lt,矩阵rini依次相乘得到含有运动初始状态和旋转，张动，进动的矩阵b，b＝rini×lt×rmic，设b矩阵的逆矩阵为矩阵f，通过三次矩阵转换将高速微动运动的目标化为相对静止目标来计算，这样fdtd计算将更为简单；

b和f在下面步骤中参数在s坐标系和s'坐标系中数据进行转换中参与计算，

所述步骤四具体包括以下过程：

过程5.1，确定静止坐标系s下的入射波，定义幅值e0，频率f以及入射角度θis,

过程5.2，通过转换矩阵经过三次坐标系转换(s→p→m→s')到s'坐标系，得到s'坐标系下入射波，入射波的幅值为e0′，频率为fs′，入射角度为θi′s,

6.根据权利要求1所述的一种计算类球头柱体高速微动目标电磁散射场的方法，其特征在于所述步骤五具体包括以下过程：

过程6.1，导入已经离散化的类球头柱体；

过程6.2，输入入射波的幅值e0′，频率fs′，入射角度θi′s,

过程6.3，运用三维fdtd算法的差分计算，电场和磁场相差半个时间步长，时间步长δt′＝δx′/2c，从而在时间上迭代求解，求得在s'坐标系各个时刻空间中的电场与磁场

所述步骤六具体包括以下过程：

过程7.1,将s'坐标系各个时刻空间中的电场与磁场用矩阵rmic转换到坐标系m中得到电场与磁场

过程7.2，在使用均值进行外推之前需要对齐网格，因此，插值出外推需要的网格点的电场值和磁场值；

过程7.3，将电场值和磁场值进行外推，得到新的电场与磁场

所述步骤七具体包括以下过程：

过程8.1，由于速度v的影响，导致s坐标系和s′坐标系的网格间隔和时间步长不一致。因此需将s'坐标系的时间步长δt'和空间步长δx',δy',δz'经过三次转换(s'→m→p→s)到s坐标系的时间步长为δt，空间步长为δx,δy,δz；

过程8.2，在m坐标系得到的电场与磁场经过两次转换(m→p→s)到s坐标系下的电场和磁场

过程8.3，得到的s坐标系下的电场和磁场进行能量计算，并进行离散傅里叶变换，滤波，得到其频率分布；

8.根据权利要求1所述的一种计算类球头柱体高速微动目标电磁散射场的方法，其特征在于通过visualstudio2015和fortran2016语言进行网格建模，计算与目标的平动相对静止的坐标系下的入射波，加入射波，用fdtd计算s'坐标系的场值，并转换到表示目标自身运动的坐标系m下进行外推处理，最后转换到s坐标系下得到电磁散射场。

如图1所示为目标模型图，有一个圆柱和半球组成，半球直径为1m，圆柱高1.5m，该目标的质心在直角坐标系的原点。目标运动方向为-y方向，速度v＝0.1c，目标的章动频率f1＝4hz,进动频率f2＝2.3hz，旋转频率f3＝0，初始滚转角ψ0＝0，进动角θp＝10°，最大章动角θm＝10°，初始偏航角入射波的初始角度为θi＝0°,

如题2所示为目标的数学模型的三个平面图，分别为x＝0面，y＝0面和z＝0面。选择空间步长δx'＝δy'＝δz'＝0.025m，将实际尺寸转换到网格数(i,j,k)，若空间某网格(i,j,k)的位置满足以下条件时

当-0.5833<y≤0.9167，

x²+z²≤1

当-1.0833≤y≤-0.5833,

x²+(y+0.5833)²+z²≤1

则该网格为金属介质，否则为空气介质；依次类推，即可得到空间内的剖分网格。

如图3所示为四个坐标系的之间的关系，其中s坐标系到p坐标系的转换矩阵为矩阵rini，p坐标系转换到m坐标系的转换矩阵为矩阵lt，m坐标系到s'坐标系的转换矩阵为rmic。通过这三个转换矩阵实现入射波，电场和磁场从静止坐标系到与运动目标相对静止的坐标系进行转换。

此处选择的入射波为调制高斯脉冲源，入射波入射角为θi,由相位不变性，和坐标系s和坐标系s'之间的时空转换关系，如下面的公式所示

便可算出运动坐标系s'下的入射波的入射角θi',和相位常数kis'。

如图4所示为，在图中，在fdtd计算中设置吸收边界，散射边界以及总场边界，各值如下表：

如图5为该方法的流程图，在总场边界处引入入射波，入射波为调制高斯脉冲，坐标系s中规定入射波幅值e0＝1，τ＝3.3299998×10^-9s，t0＝3.3299998×10^-9s,角度转换s'坐标系后的中入射波的参数参考图3得到幅值e0′＝1.005038，由于目标的微动，入射角度是一直在变化。

计算电场磁场的第n步计算x方向的电场公式如下：

其中

类似进行迭代，便可获得坐标系s'下的电场与磁场

如图6所示为在坐标系m下的对齐图，若图中p,p’点是所需外推的网格点，则p点的场值为p1,p2,p3,p4四点的平均值；同理p’点的场值为p5,p6,p7,p8点的平均值。

如图7所示为目标微动下的散射波的能量图。散射波的能量计算的公式

如图8所示为将图7进行傅里叶变化和滤波等处理后的的频谱图。由图可得到的目标的频谱特性和设置的微动的频率相符合。设f＝af1+bf2(a,b为整数,章动频率f1＝4hz,进动频率f2＝2.3hz)，图6中的频率的峰值点的频率代表f，峰值点边的括号里面的数字表示(a，b)。检测到的频率，为1.1hz，1.7hz，2.3hz，4hz，4.6hz，6.3hz和6.9hz。例如，1.1hz＝2f1+(-3)f2。