电大尺寸目标的瞬态电磁特性快速提取方法与流程

本发明属于目标电磁散射特性数值计算
技术领域：
，特别是一种电大尺寸目标的瞬态电磁特性快速提取方法。
背景技术：
：电磁计算的数值方法可分为频域数值方法和时域数值方法，频域数值方法包括矩量法(MOM)和有限元法(FEM)等，时域数值方法包括时域积分方程方法(TDIE)和时域有限差分方法(FDTD)等。频域方法处理的是稳态问题，一次计算只能分析单个频率的响应，而时域方法则可以处理瞬态问题，而且时域的激励信号常为宽带信号，可以通过一次计算获得宽频带的响应。因此，在分析目标的瞬态与宽频带的电磁散射特性时，采用时域数值方法有明显的优势。然而，时域积分方程方法(TDIE)和时域有限差分方法(FDTD)可以很好地解决电小尺寸物体的散射，但在计算电大物体的散射时，对计算机的配置要求过高。因此，这类方法不能够有效地提取瞬态电磁散射特性。迭代推进方法是用于求解目标散射问题的一种比较新型的方法，抛物线方程方法属于其中之一，它是波动方程的一种近似形式，假设电磁波能量在沿着抛物线轴向的锥形区域内传播。抛物线方程方法为求解电磁散射提供了一种准确、高效的计算方法，它的主要缺陷是只能对抛物线方向近轴区域内的电磁散射进行快速、准确地计算，不过这种限制可以通过旋转抛物线轴向来克服。Claerbout在文章“FundamentalsofGeophysicalDataProcessingwithApplicationtoPetroleumProspecting”中，使用了时域抛物线方程方法分析了地震波的传播问题；Popov在文章“ModelingEMtransientpropagationoverirregulardispersiveboundary”和“Electromagneticpulsepropagationovernonuniformearthsurface:Numericalsimulation”中又将时域抛物线方程方法推广到分析了二维电磁波的传播问题。但是，现有文献中没有三维时域抛物线方程方法提取瞬态电磁散射特性的相关报道。技术实现要素：本发明的目的在于提供一种电大尺寸目标的瞬态电磁特性快速提取方法，该方法大大简化了求解矩阵，能够快速得到电磁散射特性参数。实现本发明目的的技术解决方案为：一种电大尺寸目标的瞬态电磁特性快速提取方法，步骤如下：步骤1、建立散射体的离散模型，确定时域抛物线的轴向方向作为x轴，采用网格对散射体沿抛物线的轴向方向进行离散处理，形成垂直于x轴的若干个切面，再通过所有节点与四面体的几何关系判断节点是散射体的内部点、外部点或者边界点；步骤2、构造矩阵方程：将时域抛物线方程在时间上使用CN差分格式，在空间上写为交替方向隐式的差分形式，在迭代方向上的相邻两个面之间引入一个中间面，通过前一个面的散射场场值和前一个时间刻的散射场场值按行求解中间面的散射场场值，且各行之间求解相互独立；同样，通过中间面的散射场场值和前一个时间刻的散射场场值按列求解下一个面的散射场场值，且各列之间求解相互独立，最后在散射体表面根据切向总场分量为0的方程以及抛物线方程，联立构造出矩阵方程；步骤3、令x轴方向为待求的散射方向，依次对沿轴向方向的各个切面上的离散节点散射场场值进行递推求解，通过不断更新各个切面上的边界点的位置求解下一个切面上各个离散节点处的散射场场值；对各个切面上的离散节点散射场场值进行递推，求解最后一个切面的散射场场值；步骤4、在时间方向进行迭代，重复步骤3；根据所需的观察频点对最后一个切面的散射场场值进行傅立叶变换得到频域的散射场场值，由近远场转换求解目标散射体的双站雷达散射截面积。本发明与现有技术相比，其显著优点为：(1)方程形成简单：本发明在抛物线轴向方向使用了差分格式，可将一个三维问题转化为一系列的二维问题进行求解，矩阵形成快捷简便；(2)内存需求低：在每个切面上，使用了交替方向隐式的差分格式，降低了求解矩阵的维数，对计算内存的需求降低；(3)计算速度快：本发明交替方向隐式的差分格式使得求解矩阵为三对角矩阵，可通过追赶法求解，提高计算速度。下面结合附图对本发明作进一步详细描述。附图说明图1为本发明的电大尺寸目标的瞬态电磁特性快速提取方法在每个切面逐行/列求解示意图。图2为本发明实施例中观察频率为200MHz的散射体双站RCS曲线图。图3为本发明实施例中观察频率为300MHz的散射体双站RCS曲线图。图4为本发明实施例中观察频率为400MHz的散射体双站RCS曲线图。具体实施方式结合图1，本发明的电大尺寸目标的瞬态电磁特性快速提取方法，步骤如下：步骤1、建立散射体的离散模型，确定时域抛物线的轴向方向作为x轴，采用网格对散射体沿抛物线的轴向方向进行离散处理，形成垂直于x轴的若干个切面，再通过所有节点与四面体的几何关系判断节点是散射体的内部点、外部点或者边界点；具体为：步骤1.1、对散射体进行三角面元的面剖分，确定轴方向每个切面的方程，通过剖分网格的几何关系求解三角面元与切面的交点，与该交点距离最近的差分网格点标记为边界点并求出该点法向分量；步骤1.2、对散射体进行四面体的体剖分，通过判断某点是否处于四面体内部来区分该点处于散射体内部或者散射体外部，并对这些点进行标记；以上即可完成目标的建模，为下面的矩阵构造以及求解奠定了基础；步骤2、构造矩阵方程，将时域抛物线方程在时间上使用Crank-Nicholson差分格式，在空间上写为交替方向隐式的差分形式，在迭代方向(沿着抛物线轴向方向)上的相邻两个面之间引入一个中间面，通过前一个面的场值和前一个时间刻的场值按行求解中间面的散射场场值，且各行之间求解相互独立；同样，通过中间面的场值和前一个时间刻的场值按列求解下一个面的散射场场值，且各列之间求解相互独立；最后在散射体表面根据切向电场分量为0的方程以及抛物线方程，联立构造出矩阵方程，具体步骤如下：步骤2.1、将物体全部包围的闭合区域为计算区域，设该计算区域离散为N个切面，每个切面离散为M*M个离散节点，在三维情况下，时域抛物线方程为：∂2Π∂y2+∂2Π∂z2-2c∂2Π∂x∂t=0---(1)]]>其中，∏代表散射场场值；将式(1)写成交替方向隐式的差分形式：-ry1+2ry-ryΠp-1,qn+1/2,l+Πp-1,qn+1/2,l+1Πp,qn+1/2,l+Πp,qn+1/2,l+1Πp+1,qn+1/2,l+Πp+1,qn+1/2,l+1=rz1-2rzrzΠp,q-1n,l+Πp,q-1n,l+1Πp,qn,l+Πp,qn,l+1Πp,q+1n,l+Πp,q+1n,l+1-2Πp,qn,l+2Πp,qn+1/2,l---(2)]]>-rz1+2rz-rzΠp,q-1n+1,l+Πp,q-1n+1,l+1Πp,qn+1,l+Πp,qn+1,l+1Πp,q+1n+1,l+Πp,q+1n+1,l+1=ry1-2ryryΠp-1,qn+1/2,l+Πp-1,qn+1/2,l+1Πp,qn+1/2,l+Πp,qn+1/2,l+1Πp+1,qn+1/2,l+Πp+1,qn+1/2,l+1-2Πp,qn+1/2,l+2Πp,qn+1,l]]>(3)]]>其中，和分别表示l时间刻上第n个面上分布在第p列的连续的三个点的散射场场值，和分别表示l+1时间刻上第n个面上分布在第p列的连续的三个点的散射场场值，和分别表示l时间刻上n+1/2中间面上分布在第q行的连续的三个点的散射场场值，和分别表示l+1时间刻上n+1/2中间面上分布在第q行的连续的三个点的散射场场值，和分别表示l时间刻上第n+1面上分布在第p列的连续的三个点的散射场场值，和分别表示l+1时间刻上第n+1面上分布在第p列的连续的三个点的散射场场值；△t为时间迭代步长，△x、△y和△z分别为x、y、z方向上标准网格点的长度，c为光速，1≤n≤N，1≤p≤M，1≤q≤M，1≤l≤L，L为总时间刻个数；矢量抛物线方程由x、y、z三个方向上的标量抛物线方程构成：∂2Πx∂y2(x,y,z)+∂2Πx∂z2(x,y,z)-2c∂Πx∂x∂t(x,y,z)=0∂2Πy∂y2(x,y,z)+∂2Πy∂z2(x,y,z)-2c∂Πy∂x∂t(x,y,z)=0∂2Πz∂y2(x,y,z)+∂2Πz∂z2(x,y,z)-2c∂Πz∂x∂t(x,y,z)=0---(4)]]>按照上面的推导，将式(4)写成x、y、z三个方向的交替方向隐式格式；步骤2.2、对于目标边界点，假设P(xb,yb,zb)为散射体表面上的点，n＝(nx,ny,nz)为P点的法向方向，在金属表面，切向电场为零，由则对应的频域表达式为：nxEy(xb,yb,zb)-nyEx(xb,yb,zb)＝0nxEz(xb,yb,zb)-nzEx(xb,yb,zb)＝0(5)nyEz(xb,yb,zb)-nzEy(xb,yb,zb)＝0其中，Ex(xb,yb,zb)、Ey(xb,yb,zb)、Ez(xb,yb,zb)分别为P点电场在x轴、y轴、z 轴方向上的分量；则其对应的时域表达式如下所示：nxΠy(xb,yb,zb,t)-nyΠx(xb,yb,zb,t)=-[nxEyi(xb,yb,zb,t-xb/c)-nyExi(xb,yb,zb,t-xb/c)]nxΠz(xb,yb,zb,t)-nzΠx(xb,yb,zb,t)=-[nxEzi(xb,yb,zb,t-xb/c)-nzExi(xb,yb,zb,t-xb/c)]nyΠz(xb,yb,zb,t)-nzΠy(xb,yb,zb,t)=-[nyEzi(xb,yb,zb,t-xb/c)-nzEyi(xb,yb,zb,t-xb/c)]---(6)]]>其中，分别为入射电场在x轴、y轴、z轴方向上的分量；在球坐标系中，入射电场为：Ei(r,t)=E0cos[2πf0(t-r·k^/c)]e-(t-r·k^/c-τp)22σ2---(7)]]>其中，E0为入射电场的传播方向，f0为调制的中心频率，σ＝3/(πfbw)表示脉冲宽度，fbw为频带宽度，为入射方向，r·k^=xsinθcosφ+ysinθsinφ+cosθ,]]>θ,φ分别代表球坐标系下与z轴和xoy平面的夹角，τp表示入射场的延时，r表示点(x,y,z)到原点的距离；由于式(6)构成的矩阵秩为2，不足以唯一确定方程的解，所以由该方程结合抛物线方程即可构成秩为3的矩阵，此时便可求解边界上的场值；为了提高求解速度，将矩阵写成三对角矩阵，将式(6)分成3种情况：(1)对于x方向的法向量不为0，y、z方向的法向量任意的情况：nxΠy(xb,yb,zb,t)-nyΠx(xb,yb,zb,t)=-[nxEyi(xb,yb,zb,t-xb/c)-nyExi(xb,yb,zb,t-xb/c)]nxΠz(xb,yb,zb,t)-nzΠx(xb,yb,zb,t)=-[nxEzi(xb,yb,zb,t-xb/c)-nzExi(xb,yb,zb,t-xb/c)]---(8)]]>其中，Πx为散射场场值的x方向的分量，将式中的Πx近似为前一个面的Πx，此时联立x方向的抛物线方程，可求解出一组唯一的方程的解；(2)对于x方向的法向量为0，y的法向量不为0，z方向的法向量任意的情况：-nyΠx(xb,yb,zb,t)=-[-nyExi(xb,yb,zb,t-xb/c)]nyΠz(xb,yb,zb,t)-nzΠy(xb,yb,zb,t)=-[nyEzi(xb,yb,zb,t-xb/c)-nzEyi(xb,yb,zb,t-xb/c)]]]>(9)]]>其中，Πy为散射场场值的y方向的分量，将式中的Πy近似为前一个面的Πy，此时联立y方向的抛物线方程，可求解出一组唯一的方程的解；(3)对于x和y方向的法向量为0，z方向的法向量不为0的情况，-nzΠx(xb,yb,zb,t)=-[nxEzi(xb,yb,zb,t-xb/c)-nzExi(xb,yb,zb,t-xb/c)]-nzΠy(xb,yb,zb,t)=-[nyEzi(xb,yb,zb,t-xb/c)-nzEyi(xb,yb,zb,t-xb/c)]---(10)]]>此时联立x方向的抛物线方程，可求解出一组唯一的方程的解；步骤3、在某个固定的时间刻，令x轴方向为待求的散射方向依次对沿轴向方向的各个切面上的节点散射场场值进行递推求解，通过不断更新边界点的信息以及方程的右边向量来求解下一个切面上各个离散节点处的散射场场值；对各个切面上的节点散射场场值进行递推，求解最后一个切面散射场场值，具体过程如下：步骤3.1、将前一个时间刻各个离散的节点的散射场场值，以及前一个切面各个离散的节点的散射场电场值作为当前切面求解时的右边向量；步骤3.2、在当前切面所确定的边界点处，加入切向的电场分量为0的边界条件和抛物线方程，处于物体内部的节点电场值赋值为0，形成当前切面更新后的矩阵方程；步骤3.3、求解步骤3.2中更新后的矩阵方程，方程的解即为当前时间刻上当前切面各个离散的节点的散射场场值；每个切面上，根据处于不同的位置的离散节点，带入不同的离散方程，由前一个面的散射场场值求得下一个面的散射场场值，不断递推得到最后一个切面的散射场场值；步骤4、在时间方向进行迭代，重复步骤3；根据所需的观察频点对最后一个切面的散射场场值进行傅立叶变换得到频域场值，由近远场转换求解目标散射体双站雷达散射截面积RCS；对于指定的观察频点ω，时域抛物线方程与频域时域抛物线方程解的转换关系为：u(x,y,z,ω)=[∫-∞∞Π(x,y,z,t)e-iωtdt]/F~(ω)---(11)]]>其中，表示入射场的傅立叶变换，u(x,y,z,ω)表示某节点对应频点ω的频域场值，∏(x,y,z,t)表示节点(x,y,z)的时域散射场场值，i表示纯虚数；在频域中，利用近远场转换公式计算雷达散射截面积；三维坐标系下，在(θ,φ)方向的双站RCS为：σ(θ,φ)=limr→∞4πr2|us(x,y,z)|2|ui(x,y,z)|2---(12)]]>其中us和ui分别表示散射场和入射场的电场分量，π为圆周率；θ代表球坐标系下向量(x,y,z)与z轴的夹角，φ代表球坐标系下向量(x,y,z)与xoy面的夹角。使用旋转抛物线方程方法来克服抛物线方程的近轴限制，通过旋转抛物线的轴向方向来确定各个方向的散射场，然后通过近场远推获得远区的散射场，从而确定目标的双站RCS。下面结合具体实施例对本发明作进一步说明。实施例1本实施例进行了电磁散射的典型仿真，仿真在主频2.68GHz、内存8GB的计算机上实现，以半径5m的金属球为例，入射波的方向θ＝90°，φ＝0°，入射波中心频率为300MHz，带宽为600MHz，观察频点选择200MHz，300MHz，400MHz三个频点。为了验证本发明方法的正确性，与MIE结果作为参照；图2-图4为两种电磁散射特性仿真的RCS曲线图，图2的观察频率为200MHz，图3的观察频率为300MHz，图4的观察频率为400MHz，从图中的曲线可以看出，本发明的方法与正确的数值结果吻合。表1为本专利方法和时域有限差分方法(FDTD)的时间和内存对比：表1方法内存(GB)CPU时间(s)FDTD4.118000本专利方法1.5752通过表1中本方法与现有的FDTD方法对比，本专利方法不仅节省内存，还能够快速仿真分析目标物体的电磁散射特性。综上所述，本发明将复杂的三维问题分解为很多个二维的问题进行求解，追赶法解三对角矩阵更加快了计算速度，其实现过程灵活自由，具有很强的实际工程应用价值。当前第1页1 2 3