CN104778286B - 掠海飞行器电磁散射特性快速仿真方法 - Google Patents

Abstract

本发明公开了一种掠海飞行器电磁散射特性快速仿真方法。首先建立物体的离散模型，确定抛物线的轴向为x轴，采用网格对物体沿抛物线的轴向方向进行离散，形成垂直于x轴的若干个切面，确定每个切面所切物体的边界点并判断所有节点的位置；在x轴方向使用CN差分格式获取相邻两个切面间的关系，在y轴、z轴方向采用RPIM构造形函数及空间导数，并且在散射体表面根据入射波和经过海平面反射波的关系方程以及散度方程，联立构造出矩阵方程；依次对各个切面上的节点电场值进行递推，求解最后一个切面和前一个切面的矩阵方程，根据这两个切面的电场值对x轴方向的差分确定磁场，根据互易定理确定雷达散射截面积。本发明能够快速分析半空间环境下的电磁散射特性。

Description

掠海飞行器电磁散射特性快速仿真方法

技术领域

本发明属于目标电磁散射特性数值计算技术领域，特别是一种掠海飞行器电磁散射特性快速仿真方法。

背景技术

各种计算方法包括从高频近似方法到低频数值方法，以及关于计算各种媒质散射特性的方法，都正在被广泛研究。在最近几十年内，随着计算电磁学的发展和应用，目前已出现多种分析散射问题的方法，对电磁散射的分析方法大致主要从两种途径进行，即高频近似方法、低频数值方法，二者均有其特点和局限。计算RCS的低频方法都是数值方法，按其原理分为两类，一类是微分方程法，如有限元法(FEM)和时域有限差分(FDTD)；另一类是积分方程法，如:矩量法(MOM)。有限元法主要是求解微分方程问题，此方法适应性强，能适应各种形状的边界，但它的缺点是计算量大，不适合无界区域，花费时间长。时域有限差分法对于时域的分析是有效的，但由于计算量大的原因，目前也只能用于电小尺寸物体。作为计算方法中最具代表性的矩量法，理论上，它可以适用于求解任意目标在任意入射场下的散射问题。然而在求解过程中，该算法的可行性受到实际情况的限制。一方面，在计算大型目标时，矩量法的矩阵求解复杂度太大，可能会碰到矩阵求逆不收敛的问题；另一方面，在高频区，高频局部散射特性使得目标各部分之间的相互影响较小，从而降低了应用此类过于复杂算法的必要性。因此，这类方法也常用于求解电小尺寸的物体。近年来，国内外学者开始将抛物线方程方法应用于处理电磁散射问题.该算法把波动方程简化为抛物线方程，将散射目标等效为一系列的面元或线元，然后通过散射体上的边界条件和场的空间递推方式求解抛物线方程，把三维问题转化为一系列的二维问题来计算，通过近场——远场转换得到远区散射场，进而计算目标的双站RCS。PE方法在数值方法和解析方法之间架起了一座桥梁。数值方法如矩量法(MOM)，FDTD给出了Mxawell方程的精确解；解析方法则基于射线理论或物理光学理论。

电磁计算的数值方法如有限元法(FEM)，时域有限差分方法(FDTD)等微分方法去解决半空间电大尺寸目标的散射问题时，存在对传播空间的离散导致巨大的未知量等缺陷。另一方面积分方程中将半空间格林函数与矩量法相结合，并引入快速多极子对较大尺寸目标进行计算，但是由于受到现有计算机内存和计算时间的限制，还是无法较好处理电大尺寸复杂目标的散射问题。

由上可知，精确的数值方法解决电大尺寸物体在半空间环境下的散射时存在着困难，而通过引入半空间的反射系数到PE方法中，可以快速计算掠海飞行器的散射问题，同时将无网格方法来进行对复杂目标的建模。

发明内容

本发明的目的在于提供一种掠海飞行器电磁散射特性快速仿真方法，该方法不依赖于目标的规则化网格剖分，同时加入半空间的影响不会提高未知量的增减，能够快速得到电磁散射特性参数。

实现本发明目的的技术解决方案为：一种掠海飞行器电磁散射特性快速仿真方法，步骤如下：

步骤1、建立物体的离散模型，确定抛物线的轴向方向作为x轴，采用网格对物体沿抛物线的轴向方向进行离散处理，形成垂直于x轴的若干个切面，通过求解剖分的三角形网格与切面交点确定每个切面所切物体的边界点，再通过四面体网格来判断所有节点的位置；

步骤2、构造矩阵方程，在x轴方向使用CN差分格式获取相邻两个切面间的关系，在y轴、z轴方向采用RPIM构造形函数及空间导数，并且在散射体表面根据入射波和经过海平面反射波的关系方程以及散度方程，联立构造出矩阵方程；

步骤3、令x轴方向为待求的散射方向，依次对各个切面上的节点电场值进行递推求解，通过不断更新边界点的信息以及方程的右边向量来求解下一个切面上各个离散节点处的电场值；

步骤4、对最后一个切面和前一个切面的电场值进行后处理：求解最后一个切面和前一个切面的矩阵方程，得到两个切面相同离散节点处的电场值，根据这两个切面的电场值对x轴方向的差分确定磁场，并通过电场磁场求得电磁流，根据互易定理确定雷达散射截面积。

本发明与现有技术相比，其显著优点为：（1）建立模型简单：在垂直于抛物线轴向的切面上，不需要再建立类似于FDTD的等间距规则网格，只要确定一些离散点的信息即可。（2）方程形成简单：将一个三维问题转化为一系列的二维问题进行求解，通过形函数构造矩阵方程，矩阵形成快捷简便。（3）半空间环境因素的引入并没有使未知量增加，可以快速得到电磁散射特性参数。

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

附图说明

图1是本发明某一切面上未知量分布的示意图。

图2是本发明能量沿抛物线轴向传播示意图。

图3是本发明离散节点支撑域的示意图。

图4是本发明前后两个切面边界点有交差情况处理的示意图。

图5是本发明入射场方向与矢量抛物线轴向方向示意图。

图6是本发明某掠海飞行器示意图。

图7是本发明实施例中某掠海飞行器双站RCS曲线图。

具体实施方式

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

结合附图1～6，本发明掠海飞行器电磁散射特性快速仿真方法，步骤如下：

步骤1、建立物体的离散模型，确定抛物线的轴向方向作为x轴，采用网格对物体沿抛物线的轴向方向进行离散处理，形成垂直于x轴的若干个切面，通过求解剖分的三角形网格与切面交点确定每个切面所切物体的边界点，再通过四面体网格来判断所有节点的位置，具体步骤如下：

步骤1.1、在每个切面上任意选取规则分布的参考点，所选参考点用于在x轴方向上进行差分以及在y轴、z轴方向进行无网格形函数的构造；

步骤1.2、对散射体进行三角面元的面剖分，确定轴方向每个切面的方程，通过几何关系求解三角面元与切面的交点，并将交点标记为每个切面上散射体的边界点；

步骤1.3、对散射体进行四面体的体剖分，通过判别参考点是否处于四面体内部区分参考点处于散射体内部或者散射体外部，并对这些不同位置的参考点进行标记。

首先，在每个(y-z)切面上选取一些分布均匀的参考点，这些参考点用作于两个切面间的插值以及构造形函数时使用。

用三角形面元对物体进行面剖分，获取物体表面的一些离散的节点信息。垂直于x轴即为抛物线轴向，形成很多切面，这些切面与三角形相交，通过节点的几何信息求解出与切面的交点，将这些交点作为散射体在当前切面的边界点。同时对散射体进行四面体的体剖分，对每个切面上的参考点进行循环判断，看该点是否处于某个四面体的内部，如果该点处于四面体的内部则认为该点为散射体的内部点，否则认为该点处在空气层。认为离空气盒边界一定距离的点为PML层内的参考点。

通过上面的方法可得到各个切面上物体边界的节点，结合每个面上散射体外的参考点，构成了一个切面上总的未知量,各个切面的未知量分别由每个面上散射体外部固有的离散参考点和边界点相加得到。某个切面上未知量的分布示意图如图1所示，根据各个点的几何位置关系以及坐标关系确定出点所在的位置的属性，具体判断准则如下所示：

第1：离切面的上下左右边缘1个波长的节点都设置为PML的性质；

第2：由上述方法找到的交点即为物体的边界点，由边界点连成的轮廓线内的所有节点为物体的内部节点，这些参考点不作为当前面的未知量；

第3：其余的节点即为空气层的离散节点。

以上即可完成目标的建模，为下面的矩阵构造以及求解奠定了基础。

步骤2、构造矩阵方程，在x轴方向使用CN差分格式获取相邻两个切面间的关系，在y轴、z轴方向采用RPIM构造形函数及空间导数，并且在散射体表面根据入射波和经过海平面反射波的关系方程以及散度方程，联立构造出矩阵方程，具体步骤如下：

首先，我们给出三维标量波动方程：

其中，E代表电场分量，k为波数。取x轴方向为抛物线的轴方向，定义沿x方向传播的波函数，如图2所示：

u(x，y，z)＝e^-ikxE(x，y，z) (2)

将式（2）带入式（1），可得：

可将其分解为：

其中微分算子Q为：

我们只取前向抛物线形式，并利用Q的泰勒展开式，可得小角度抛物线方程：

在三维情况下，标准矢量抛物线方程可表示为：

其中，分别为波函数在x轴、y轴、z轴方向的 分量，分别为电场在x轴、y轴、z轴方向的分量，k为波数，i为虚数。对x轴方向的 求导由CN差分获得，对y轴、z轴方向的求导采用RPIM构造形函数及其空间导数，电场u(x,y, z)通过形函数u(x,y,z)＝Φ(x,y,z)U_S(x,y,z)展开，U_S(x,y,z)为待求的电场系数,Φ(x,y, z)＝[Φ₁(x,y,z),Φ₂(x,y,z),...,Φ_N(x,y,z)]为形函数，如图3所示N为支撑域内离散节 点的个数，对u(x,y,z)关于y轴和z轴的求导可以通过对Φ(x,y,z)求导实现，上式可离散成 如下形式：

其中，Δx代表前后两个切面的间距，在PML媒质中，相应的矢量抛物线方程可表示为：

式中，σ()代表电损耗的函数，σ₀代表电损耗的系数，δ代表趋肤深度的系数。对x轴方向的求导由CN差分获得，对y轴、z轴方向的求导采用RPIM构造形函数及其空间导数，电场u(x,y,z)通过形函数u(x,y,z)＝Φ(x,y,z)U_S(x,y,z)展开，U_S(x,y,z)为待求的电场系数，Φ(x,y,z)＝[Φ₁(x,y,z),Φ₂(x,y,z),...,Φ_N(x,y,z)]为形函数，N为支撑域内离散节点的个数，对u(x,y,z)关于y和z的求导可以通过对Φ(x,y,z)求导实现，上式可离散成如下形式：

通过式（10）即可构造前后两个切面上电场值U_S,x(x,y,z),U_S，y(x,y,z),U_S,z(x,y,z)与U_S,x(x+Δx,y,z),U_S,y(x+Δx,y,z),U_S,z(x+Δx,y,z)的关系的矩阵方程。

矩阵方程金属边界条件的添加以及递推求解，具体步骤如下：

对于物体边界点，假设P为散射体表面上的点，n=(n_x,n_y,n_z)为P点的法向方向，在金属表面，切向电场为零，由即：

n_xE_y(P)-n_yE_x(P)＝0

n_xE_z(P)-n_zE_x(P)＝0（11）

n_yE_z(P)-n_zE_y(P)＝0

E_x(P)、E_y(P)、E_z(P)分别为P点电场在x轴、y轴、z轴方向上的分量；考虑海平面反射波的影响，式（11）中的电场分量表示为入射场、散射场和反射场之和，如下：

分别为P点入射场在x轴、y轴、z轴方向上的分量； 分别为P点散射射场在x轴、y轴、z轴方向上的分量； 分别为P点反射场在x轴、y轴、z轴方向上的分量；

定义场量为x轴方向传播波函数，则进行如下变换：

由式（11）、（12）、（13）得对应的三个方程：

在球坐标系下：

上式中的入射场和反射场分别为：

式（14）的三个方程并不是相互独立的，其系数矩阵的秩为2，没有定解，只有加上Maxwell的散度方程，才可构成系数矩阵秩为3的线性方程组，解具有唯一性。

将对应的抛物线方程代入，P点的三维坐标下的散度方程变为：

对电场u_x(x,y,z)、u_y(x,y,z)以及u_z(x,y,z)采用RPIM构造形函数及其空间导数；电场u(x,y,z)通过形函数u(x,y,z)＝Φ(x,y,z)US(x,y,z)展开，U_S(x,y,z)为待求的电场系数,Φ(x,y,z)＝[Φ₁(x,y,z),Φ₂(x,y,z),...,Φ_N(x,y,z)]为形函数，N为支撑域内离散节点的个数，对u(x,y,z)关于y和z的求导可以通过对Φ(x,y,z)求导实现，上式可离散成如下形式：

将式（14）与式（19）联立，构造系数矩阵秩为3的线性方程组，将耦合关系填入到矩阵方程中，即可完成金属边界条件的添加。

综上所述，构造最终矩阵方程：

步骤3、令x轴方向为待求的散射方向，依次对各个切面上的节点电场值进行递推求解，通过不断更新边界点的信息以及方程的右边向量来求解下一个切面上各个离散节点处的电场值；具体过程如下：

步骤3.1、将前一个切面各个离散的节点的电场值作为当前切面求解时的右边向量；

步骤3.2、在当前切面所确定的边界点处，加入切向分量为0以及散度为0的边界条件，处于物体内部的节点电场值赋值为0，形成当前切面更新后的矩阵方程；

步骤3.3、求解步骤3.2中更新后的矩阵方程，方程的解即为当前切面各个离散的节点的电场值。

每个切面的未知量的个数是参考点的个数加上本切面边界点的个数，根据处于不同的位置，带入不同的离散方程，由前一个面的电场值求得下一个面的电场值，不断递推得到最后一个切面的电场值。

对于前后两个切面如果有重叠型区域的出现，如图4所示。对于同时处于两个切面的边界轮廓外的参考点直接将参考点上的场值赋值给下一个面的参考点；对于处于前一个切面边界轮廓内同时处于当前切面边界轮廓外的参考点视其为边界点使用阻抗边界条件进行处理；对于处于前一个切面边界轮廓外同时处于当前切面边界轮廓内的参考点则不视为未知量；对于当前面的边界点直接填入阻抗边界条件的方程。

步骤4、对最后一个切面和前一个切面的电场值进行后处理：求解最后一个切面和前一个切面的矩阵方程，得到两个切面相同离散节点处的电场值，根据这两个切面的电场值对x轴方向的差分确定磁场，并通过电场磁场求得电磁流，根据互易定理确定雷达散射截面积，具体步骤如下：

首先，已知一个切面上的各个点上的电场值E，由麦克斯韦方程求解出来H，其中j表示虚数，ω为角频率，μ为磁导率；

其次，由第一等效原理可得等效面上的电流以及磁流：

J＝e_n×H，M＝E×e_n (21)

其中，e_n为等效面的外法相分量，E为电场值，H为磁场值；

最后，远区散射场可以通过互易定理很方便地求得：

∫∫∫(E^sc·J₂-H^sc·M₂)dv＝∫∫_S(J·E₂-M·H₂)dS (22)

其中，E^sc为空间任意点处的散射电场，H^sc为空间任意点处的散射磁场，J为产生散射场的电流源，M为产生散射场的磁流源，E₂为入射波电场，H₂为入射波磁场，J₂为产生入射场的电流源，M₂为产生入射场的磁流源；

化简可得：

其中，为空间任意点处的散射电场的θ和方向的分量，μ₀为自由空间的磁导率，k₀为自由空间波数，π为圆周率，r为观察距离，为入射电场的θ和方向的分量，为入射磁场的θ和方向的分量；

上式都是在球坐标系下：

则散射场可表示成：

R^TM,R^TE分别为入射平面波为横电波和横磁波的反射系数，η为特性阻抗，ω为自由空间的相速度，μ₀为自由空间的磁导率，k₀为自由空间波数，k_p为入射波的传播方向，k_r为反射波的传播方向。

其中，

三维坐标系下，在(θ,φ)方向的双站RCS为：

其中E^s和Eⁱ分别表示散射场和入射场的电场分量，π为圆周率。

矢量抛物线方法充分考虑了极化的影响，将对波动方程的求解转换成对抛物线方程的求解，结合适当的边界条件，利用小角度矢量抛物线的形式，每个矢量抛物线方程计算出沿抛物线轴向方向大小不超过15°的锥形范围内的散射场。如图5所示，通过旋转抛物线的轴向方向来计算各个方向的散射场，然后通过近场远推获得远区的散射场，从而计算得到目标的双站RCS。

实施例1

本实施例进行了具有耗介质电磁散射的典型仿真，仿真在主频2.83GHz、内存3.5GB的个人计算机上实现，以半径为2.5m,高度为4m某圆柱飞行器为例，如图6所示，入射波频率为300MHz，入射波的方向θ=0°，，半空间介电常数为（2.3，-0.6），目标距离海平面0.2个波长，为了验证本发明方法的正确性，以商业软件FEKO仿真结果作为参照。图7为两种电磁散射特性仿真的RCS曲线图，从图中的曲线可以看出，本文方法与正确的数值结果吻合，并且时间上面具有明显的优势，此方法只需要10分钟左右的时间，而FEKO计算时间需要接近2个小时，说明本文方法能够快速仿真分析半空间环境下目标物体的电磁散射特性。

综上所述，本发明将复杂的三维问题分解为很多个二维的问题进行求解，并且不依赖于网格的规范性，半空间环境因素的加入对未知量没有直接的影响，便可对其进行快速的电磁散射仿真，其实现过程灵活自由，具有很强的实际工程应用价值。

Claims (4)

1.一种掠海飞行器电磁散射特性快速仿真方法，其特征在于，步骤如下：

步骤1、建立物体的离散模型，确定抛物线的轴向方向作为x轴，采用网格对物体沿抛物线的轴向方向进行离散处理，形成垂直于x轴的若干个切面，通过求解剖分的三角形网格与切面交点确定每个切面所切物体的边界点，再通过四面体网格来判断所有节点的位置；

步骤2、构造矩阵方程，在x轴方向使用CN差分格式获取相邻两个切面间的关系，在y轴、z轴方向采用RPIM构造形函数及空间导数，并且在散射体表面根据入射波和经过海平面反射波的关系方程以及散度方程，联立构造出矩阵方程；

步骤3、令x轴方向为待求的散射方向，依次对各个切面上的节点电场值进行递推求解，通过不断更新边界点的信息以及方程的右边向量来求解下一个切面上各个离散节点处的电场值；

步骤4、对最后一个切面和前一个切面的电场值进行后处理：求解最后一个切面和前一个切面的矩阵方程，得到两个切面相同离散节点处的电场值，根据这两个切面的电场值对x轴方向的差分确定磁场，并通过电场磁场求得电磁流，根据互易定理确定雷达散射截面积；

步骤1所述建立物体的离散模型，具体包括以下步骤：

步骤1.1、在每个切面上任意选取规则分布的参考点，所选参考点用于在x轴方向上进行差分以及在y轴、z轴方向进行无网格形函数的构造；

步骤1.2、对散射体进行三角面元的面剖分，确定轴方向每个切面的方程，通过几何关系求解三角面元与切面的交点，并将交点标记为每个切面上散射体的边界点；

步骤1.3、对散射体进行四面体的体剖分，通过判别参考点是否处于四面体内部区分参考点处于散射体内部或者散射体外部，并对这些不同位置的参考点进行标记。

2.根据权利要求1所述的掠海飞行器电磁散射特性快速仿真方法，其特征在于，步骤2中所述构造矩阵方程，具体包括以下步骤：

步骤2.1、在三维情况下，标准矢量抛物线方程表示为：

式中，分别为波函数在x轴、y轴、z轴方向的分量，分别为电场在x轴、y轴、z轴方向的分量，k为波数，i为虚数；其中，对y轴、z轴方向的求导采用RPIM构造形函数及空间导数，对x轴方向的求导采用差分；

步骤2.2、在PML媒质中，矢量抛物线方程表示为：

式中，σ()代表电损耗的函数，σ₀代表电损耗的系数，δ代表趋肤深度的系数；

对x轴方向的求导采用差分，对y轴、z轴方向的求导采用RPIM构造形函数及其空间导数；

步骤2.3、对于目标边界点，假设P为散射体表面上的点，n＝(n_x,n_y,n_z)为P点的法向方向，在金属表面切向电场为零，由将电场用各个分量来表示：

E_x(P)、E_y(P)、E_z(P)分别为P点电场在x轴、y轴、z轴方向上的分量；考虑海平面反射波的影响，式(3)中的电场分量表示为入射场、散射场和反射场之和，如下：

分别为P点入射场在x轴、y轴、z轴方向上的分量； 分别为P点散射场在x轴、y轴、z轴方向上的分量；分别为P点反射场在x轴、y轴、z轴方向上的分量；

定义场量为x轴方向传播波函数，则进行如下变换：

由式(3)、(4)、(5)得对应的三个方程：

式(6)是秩为2的方程组，不能唯一确定边界条件，引入散度方程使方程组具有唯一的解，P点的三维坐标下的散度方程变换为：

对电场u_x(x,y,z)、u_y(x,y,z)以及u_z(x,y,z)采用RPIM构造形函数及其空间导数；

综上所述，构造最终矩阵方程：

3.根据权利要求1所述的掠海飞行器电磁散射特性快速仿真方法，其特征在于，步骤3中所述对各个切面上的节点电场值进行递推求解，具体过程如下：

步骤3.1、将前一个切面各个离散的节点的电场值作为当前切面求解时的右边向量；

步骤3.2、在当前切面所确定的边界点处，加入切向分量为0以及散度为0的边界条件，处于物体内部的节点电场值赋值为0，形成当前切面更新后的矩阵方程；

步骤3.3、求解步骤3.2中更新后的矩阵方程，方程的解即为当前切面各个离散的节点的电场值。

4.根据权利要求1所述的掠海飞行器电磁散射特性快速仿真方法，其特征在于，

步骤4所述雷达散射截面积的表达式为：

三维坐标系下，在(θ,φ)方向的双站RCS为：

其中E^s和Eⁱ分别表示散射场和入射场的电场分量，π为圆周率。