CN103279600B - 基于积分边界的介质粗糙面有限元电磁仿真方法 - Google Patents

Abstract

本发明公开一种基于积分边界的介质粗糙面有限元电磁仿真方法，主要解决积分方程类方法难以求解非均匀介质和近似吸收边界的近似性问题。其实现步骤是：采用蒙特卡洛方法获得介质粗糙面；通过人工边界对介质粗糙面进行截断处理并进行离散；在仿真区域内部，根据泛函方法获得仿真区域的有限元方程；在仿真区域外部，采用积分方程方法得到积分边界条件，并通过连续性条件形成电磁耦合方程；通过求解电磁耦合矩阵方程，最终获得介质粗糙面的雷达散射参数。本发明与现有技术相比，不仅比积分方程类方法能更容易地用于仿真非均匀介质电磁问题，而且与近似吸收边界相比能提高电磁仿真精度，可用于获取介质粗糙面无界区域问题的电磁散射参数。

Description

基于积分边界的介质粗糙面有限元电磁仿真方法

技术领域

本发明属于雷达电磁仿真领域，具体涉及电磁散射数值仿真，用于获取介质粗糙面等无界区域的电磁散射参数。

背景技术

随着雷达技术的快速发展，粗糙面的电磁散射研究在理论分析与实际应用中具有重要的意义。当机载或者星载雷达对粗糙面等背景中的目标进行电磁探测时，雷达的回波信号中包含了粗糙面背景的电磁散射信号，通过分析这些回波信号，进而可以得出粗糙面以及探测目标的某些特性。因而，粗糙面雷达散射特性的研究在国防领域和民用领域具有显著的学术价值和广泛的应用前景。

在过去的几十年中，许多电磁仿真技术被学者提出用以处理地海面的电磁散射问题，大致分为近似方法和数值方法。近似方法的优点是消耗内存低、分析速度快，然而近似方法都是基于特定的物理近似，其精度往往较低并且不具有通用性。虽然数值方法消耗内存较大，分析速度较慢，然而与近似方法相比保持了较高的仿真精度。基于积分方程方法的矩量法、快速多级子方法等等虽然仿真精度较高，但对于非均匀介质问题，传统的基于积分方程的数值方法难以处理。利用有限元方法讨论无界区域问题，往往需要引入截断边界条件，将无限大区域截断为有限大小的仿真区域。在以往的地海面问题中，完全匹配层、吸收边界条件往往用来作为人工边界，然而采用这些近似边界条件分析精度较低，而且往往需要设置在离模型较远的距离处，这样就大大增加了仿真问题所需的未知量。而作为最精确的边界条件，积分边界条件仅仅用于有限大小目标的电磁仿真中，未见基于积分边界的有限元方法应用于介质粗糙面雷达散射参数仿真中。

发明内容

本发明的目的在于针对上述已有技术的不足，提供一种基于积分边界的介质粗糙面有限元电磁仿真方法，以将有限元方法应用到介质粗糙面雷达散射参数的电磁仿真中并保持较好的仿真精度。

实现本发明目的的技术方案，包括如下步骤：

（1）根据介质粗糙面的功率谱密度和粗糙面参数，通过蒙特卡洛方法产生仿真的介质粗糙面；

（2）通过人工边界将介质粗糙面仿真区域进行截断处理，然后对介质粗糙面仿真区域进行三角面元离散，并对离散三角面元进行编号，将该编号作为仿真区域的面元编号；

（3）根据面元编号的仿真区域，获得介质粗糙面离散仿真区域的有限元方程和人工边界处的积分边界条件，并将获得的有限元方程与积分边界条件结合成一个电磁耦合方程组；

（4）将电磁耦合方程组转化为电磁耦合矩阵，并求解电磁耦合矩阵获得介质粗糙面仿真区域的雷达散射参数。

本方法由于采用积分边界条件作为人工边界，避免了近似吸收边界的吸收不完全性，提高了电磁数值仿真的准确性，保留了有限元方法的优点，并可以将积分边界设置在离介质粗糙面很近的距离上甚至在粗糙面表面；同时由于本发明需要对粗糙面仿真区域进行离散，所以特别适合用于介质粗糙面雷达散射参数的电磁仿真。

附图说明

图1是本发明的实现流程图；

图2是本发明中对介质粗糙面仿真区域进行截断后的示意图；

图3是用本发明与现有矩量法仿真的介质粗糙面表面总场对比图；

图4是用本发明与现有矩量法仿真的介质粗糙面雷达散射系数对比图。

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤1：根据介质粗糙面的功率谱密度和粗糙面参数，通过蒙特卡洛方法产生仿真的介质粗糙面：

（1.1）通过实验获得所要仿真介质粗糙面的均方根高度、相关长度等粗糙面参数和功率谱密度，并选取一系列振幅独立的高斯随机变量的高斯谐波；

（1.2）根据获得的介质粗糙面的均方根高度、相关长度等粗糙面参数和功率谱密度，选取N个离散采样点x_-N/2+1，…，x_-1，x₀，x₁，…，x_N/2，在离散点处通过蒙特卡洛方法将大量的振幅独立的高斯谐波进行叠加，然后作逆傅里叶变换生成所要仿真的介质粗糙面。

所述的蒙特卡洛方法是通过功率谱密度产生随机粗糙面的数值仿真方法，具体内容可见TsangL、JAKong和KHDing所写著作《ScatteringofElectromagneticWaves:NumericalSimulations》。

步骤2：通过人工边界将介质粗糙面仿真区域进行截断处理。

参照图2，本步骤的具体实现如下：

（2.1）将介质粗糙面上方的人工边界设置在介质粗糙面轮廓上以减小分析区域，并将下方人工边界设置在离介质粗糙面水平面约0.5λ处；

（2.2）在介质粗糙面仿真区域上方的人工边界处引入锥形入射波，其入射场Φ^inc(r)表示如下

Φ^inc(r)＝exp[-jk·r(1+w(r))]·exp[-(x-yctgθ_inc)²/g²]，

其中w(r)＝exp[2(x-yctgθ_inc)²/g²-1]/(kgsinθ_inc)²，k是空间波数，是空间位置矢量，g是锥形入射波因子，θ_inc为入射角。

步骤3：对介质粗糙面仿真区域进行三角面元离散，并对离散三角面元进行编号，将该编号作为仿真区域的面元编号：

（3.1）将产生的介质粗糙面导入模型处理有限元商用软件ANSYS，用边长为0.1λ左右的三角面元剖分介质粗糙面仿真区域，并导出介质粗糙面仿真区域数据文件，获得介质粗糙面三角面元离散仿真区域；

（3.2）引入整数数组n(e,i)对介质粗糙面仿真区域的离散三角面元编号进行记录，在联系数组n(e,i)中，e是单元编码，i是结点的局部编码，n(e,i)的值是结点的全局编码。

步骤4：根据仿真区域单元编号，获得介质粗糙面离散仿真区域的有限元方程和人工边界处的积分边界条件，并将获得的有限元方程与积分边界条件结合成为一个电磁耦合方程。

（4.1）通过泛函方法，在仿真区域内部获得介质粗糙面的等效有限元方程：

δF(Φ)＝0，

其中,F(Φ)是仿真区域电磁场的泛函表示式，其具体表达式如下

$\left(\mathrm{\Φ}\right)=\frac{1}{2}{\∫\∫}_{{\mathrm{\Ω}}_{\mathrm{in}}}[\frac{1}{\mathrm{\ρ}}{\left(\frac{\∂\mathrm{\Φ}}{\∂x}\right)}^2\frac{1}{\mathrm{\ρ}}{\left(\frac{\∂\mathrm{\Φ}}{\∂x}\right)}^2{k_0}^2\mathrm{\υ}{\mathrm{\Φ}}^2]\mathrm{d\Ω}-{\∫\∫}_{{\mathrm{\Ω}}_{\mathrm{in}}}\mathrm{f\Φd\Ω},$

$+{\∫}_{{\mathrm{\Γ}}_{\mathrm{UBIE}}^{-}}\mathrm{\Φ}{\mathrm{\ψ}}_U\mathrm{d\Γ}+{\∫}_{{\mathrm{\Γ}}_{\mathrm{DBIE}}^{-}}\mathrm{\Φ}{\mathrm{\ψ}}_D\mathrm{d\Γ}$

式中，Φ表示场值，ψ表示场法向导数，Ω_in表示仿真区域内部，代表粗糙面轮廓上的积分边界，表示介质粗糙面下方的积分边界，上标-表示人工边界内侧，式中符号ρ和υ分别表示仿真区域的介电常数ε_r和磁常数μ_r，当入射波为横磁波时，ρ＝μ_r,υ＝ε_r；当入射波为横电波时，ρ＝ε_r,υ＝μ_r；

（4.2）通过自由空间积分方程方法，得到粗糙面仿真区域上方人工边界Γ_UBIE的积分边界条件，表示如下

$\mathrm{\Φ}\left(r^{\′}\right)={\mathrm{\Φ}}^{\mathrm{inc}}\left(r^{\′}\right)+{\∫}_{{\mathrm{\Γ}}_{\mathrm{UBIE}}^{+}}[\mathrm{\Φ}\left(r\right)\frac{\∂G_U(r,r^{\′})}{\∂n}-G_U(r,r^{\′})\frac{\∂}{\∂n}\mathrm{\Φ}\left(r\right)]\mathrm{d\Γ},$

其中，G_U(r,r′)为上方自由空间中的格林函数，上标+表示人工边界外侧。

（4.3）通过介质空间积分方程方法，得到粗糙面仿真区域下方人工边界Γ_DBIE的积分边界条件，表示如下

$\mathrm{\Φ}\left(r^{\′}\right)={\∫}_{{\mathrm{\Γ}}_{\mathrm{DBIE}}^{+}}[\mathrm{\Φ}\left(r\right)\frac{{\∂G}_D(r,r^{\′})}{\∂n}-G_D(r,r^{\′})\frac{\∂}{\∂n}\mathrm{\Φ}\left(r\right)]\mathrm{d\Γ},$

其中，G_D(r,r′)为下方介质空间中的格林函数；

（4.4）根据连续性边界条件，将仿真区域内部有限元方程和上下人工边界的积分边界条件耦合成电磁耦合方程组：

$\mathrm{\δ}\left\{\begin{array}{c}\frac{1}{2}{\∫}_{{\mathrm{\Ω}}_{\mathrm{in}}}[\frac{1}{\mathrm{\ρ}}{\left(\frac{\∂\mathrm{\Φ}}{\∂x}\right)}^2+\frac{1}{\mathrm{\ρ}}{\left(\frac{\∂\mathrm{\Φ}}{\∂y}\right)}^2-{k_0}^2\mathrm{\υ}{\mathrm{\Φ}}^2]\mathrm{d\Ω}-{\∫}_{{\mathrm{\Ω}}_{\mathrm{in}}}\mathrm{f\Φd\Ω}\\ +{\∫}_{{\text{\Γ}}_{\mathrm{UBIE}}^{-}}\mathrm{\Φ}{\mathrm{\ψ}}_U\mathrm{d\Γ}+{\∫}_{{\mathrm{\Γ}}_{\mathrm{DBIE}}^{-}}\mathrm{\Φ}{\mathrm{\ψ}}_D\mathrm{d\Γ}\end{array}\right\}=0,$

$\mathrm{\Φ}\left(r^{\′}\right)-{\∫}_{{\mathrm{\Γ}}_{\mathrm{UBIE}}^{-}}[\mathrm{\Φ}\left(r\right)\frac{\∂G_U(r,r^{\′})}{\∂n}-G_U(r,r^{\′})\frac{{\mathrm{\ρ}}^{+}}{{\mathrm{\ρ}}^{-}}\frac{\∂}{\∂n}\mathrm{\Φ}\left(r\right)]\mathrm{d\Γ}={\mathrm{\Φ}}^{\mathrm{inc}}\left(r^{\′}\right),$

$\mathrm{\Φ}\left(r^{\′}\right)-{\∫}_{{\mathrm{\Γ}}_{\mathrm{DBIE}}^{-}}[\mathrm{\Φ}\left(r\right)\frac{{\∂G}_D(r,r^{\′})}{\∂n}-G_D(r,r^{\′})\frac{{\mathrm{\ρ}}^{+}}{{\mathrm{\ρ}}^{-}}\frac{\∂}{\∂n}\mathrm{\Φ}\left(r\right)]\mathrm{d\Γ}=0,$

其中，连续性边界条件表示如下

$\mathrm{\Φ}{|}_{{\mathrm{\Γ}}^{+}}=\mathrm{\Φ}{|}_{{\mathrm{\Γ}}^{-}},$

$\frac{1}{{\mathrm{\ρ}}^{+}}\frac{\∂\mathrm{\Φ}}{\∂n}{|}_{{\mathrm{\Γ}}^{+}}=\frac{1}{{\mathrm{\ρ}}^{-}}\frac{\∂\mathrm{\Φ}}{\∂n}{|}_{{\mathrm{\Γ}}^{-}},$

式中，Γ⁺表示从仿真区域外部靠近截断边界，Γ^-是从仿真区域内部靠近截断边界。

步骤5：将步骤4中获得的电磁耦合方程组转化为电磁耦合矩阵，并求解电磁耦合矩阵，得到介质粗糙面仿真区域的雷达散射参数。

（5.1）根据线性插值，得到线型插值基函数N(x,y)，从而获得仿真区域面单元和边界单元处的场Φ(x,y)和场法向导数的插值形式：

${\mathrm{\Φ}}^e(x,y)=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{N}_i^e(x,y){\mathrm{\Φ}}_i^e,$

${\mathrm{\Φ}}^s(x,y)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{N}_i^s(x,y){\mathrm{\Φ}}_i^s,$

$\frac{{\∂\mathrm{\Φ}}^s}{\∂n}(x,y)=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}{N}_i^s(x,y)\frac{{\∂\mathrm{\Φ}}_i^s}{\∂n},$

其中，上标e表示面元编号，s是边界单元编号；

（5.2）将上述仿真区域各单元的场和场法向导数代入电磁耦合方程组中，获得仿真区域的电磁耦合矩阵：

$\left[\begin{array}{ccc}{K}_{\mathrm{II}}& K_{\mathrm{IU}}& K_{\mathrm{ID}}\\ K_{\mathrm{UI}}& K_{\mathrm{UU}}& 0\\ K_{\mathrm{DI}}& K_{\mathrm{DD}}& 0\end{array}\right]\left\{\begin{array}{c}{\mathrm{\Φ}}_I\\ {\mathrm{\ψ}}_U\\ {\mathrm{\ψ}}_D\end{array}\right\}=\left\{\begin{array}{c}0\\ b_U\\ 0\end{array}\right\},$

式中，ρ⁺表示仿真区域外部的电磁常数，ρ^-表示仿真区域内部的电磁常数，K是子矩阵元素，b是入射波的列矩阵，下标I表示仿真区域内部，U表示上方人工边界，D表示下方人工边界；

（5.3）通过LU矩阵分解方法求解电磁耦合矩阵方程，得到仿真区域各结点的场值和人工边界上的场法向导数值，并将人工边界处的场值和场法向导数值代入双站雷达散射参数方程，得出介质粗糙面的双站雷达散射参数σ(θ)：

$\mathrm{\σ}\left(\mathrm{\θ}\right)=10\lg \left(2\mathrm{\πr}{\left|\frac{{\mathrm{\Φ}}^{\mathrm{scat}}}{{\mathrm{\Φ}}^{\mathrm{inc}}}\right|}^2\right),$

其中，θ是雷达的散射角度，r是仿真位置离原点的距离，Φ^scat是散射场，Φ^inc表示入射场。

本发明的试验仿真效果可以通过以下试验进一步说明

（1）试验仿真条件

介质粗糙面电磁模型的参数选择如下：观测频率f＝300MHz，入射波波长λ＝1.0m，离散的三角单元边长约为0.1m，粗糙面的长度L为25.6m，其他参量例如入射角θ_inc,锥形波因子g,粗糙面的相对介电常数ε_r,均方根高度rms,以及相关长度lc均在图中给出。由于粗糙面的轮廓的随机性，采用30个样本进行平均从而得到稳定的双站雷达散射参数。

（2）试验仿真结果分析

仿真实验1，利用本发明和现有矩量法MoM对横磁波TM极化和横电波TE极化情况下介质粗糙面的表面场进行了电磁仿真，结果如图3，其中图3(a)给出了横磁波TM极化情况下单个介质粗糙面样本表面场与矩量法的对比结果，图3(b)给出了横电波TE极化情况下单个介质粗糙面样本的表面场与矩量法的对比结果。

从图3(a)和图3(b)的对比可以看出，本发明仿真的介质粗糙面表面场分布与现有矩量法的结果吻合较好，可以看出本发明在介质粗糙面表面场的仿真中具有良好的准确性。

仿真实验2，利用本发明和现有矩量法MoM对横磁波TM极化和横电波TE极化情况下介质粗糙面的双站雷达散射参数进行电磁仿真，结果如图4，其中图4(a)为横磁波TM极化情况下介质粗糙面的双站雷达散射参数与矩量法的对比结果，图4(b)为横电波TE极化情况下介质粗糙面的双站雷达散射参数与矩量法的对比结果。

从图4(a)和图4(b)的对比可以看出，本发明仿真的介质粗糙面的双站雷达散射参数同样与矩量法获得的结果基本一致，可以看出本发明可以有效地用以获得介质粗糙面的雷达散射参数。

从图3和图4的实例对比可以看出，本发明在介质粗糙面的电磁散射数值仿真中结果保持较好的仿真精度，无论是表面场还是雷达双站散射参数，本方法的仿真结果与现有矩量法相比数据曲线均基本吻合，能有效地获取介质粗糙面的雷达散射参数。

以上描述仅是本发明的一个具体实例，显然对于本领域内的专业人员来说，在了解了本发明内容和原理后，可以进行形式上和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。

Claims (6)

1.一种基于边界积分的介质粗糙面有限元电磁仿真方法，包括如下步骤：

(1)根据介质粗糙面的功率谱密度和粗糙面参数，通过蒙特卡洛方法产生仿真的介质粗糙面；

(2)通过人工边界将介质粗糙面仿真区域进行截断处理，然后对介质粗糙面仿真区域进行三角面元离散，并对离散三角面元进行编号，将该编号作为仿真区域的面元编号；

(3)根据面元编号的仿真区域，获得介质粗糙面离散仿真区域的有限元方程和人工边界处的积分边界条件，并将获得的有限元方程与积分边界条件结合成一个电磁耦合方程组:

(3a)通过在仿真区域内部利用面元编号结合泛函方法获得的有限元方程：

$\mathrm{\δ}\left\{\begin{array}{c}\frac{1}{2}\∫{\∫}_{{\mathrm{\Ω}}_{in}}\[\frac{1}{\mathrm{\ρ}}{\left(\frac{\∂\mathrm{\Φ}}{\∂x}\right)}^2+\frac{1}{\mathrm{\ρ}}{\left(\frac{\∂\mathrm{\Φ}}{\∂y}\right)}^2-{k_0}^2{\mathrm{\υ\Φ}}^2\]d\mathrm{\Ω}-\∫{\∫}_{{\mathrm{\Ω}}_{in}}f\mathrm{\Φ}d\mathrm{\Ω}\\ +{\∫}_{{\mathrm{\Γ}}_{UBIE}^{-}}{\mathrm{\Φ\ψ}}_{U}d\mathrm{\Γ}+{\∫}_{{\mathrm{\Γ}}_{DBIE}^{-}}{\mathrm{\Φ\ψ}}_{D}d\mathrm{\Γ}\end{array}\right\}=0,$

式中，Φ表示场值，ψ表示场法向导数，Ω_in表示仿真区域内部，代表粗糙面轮廓上的积分边界，表示介质粗糙面下方的积分边界，上标-表示人工边界内侧，式中符号ρ和υ分别表示仿真区域的介电常数ε_r和磁常数μ_r，当入射波为横磁波时，ρ＝μ_r,υ＝ε_r；当入射波为横电波时，ρ＝ε_r,υ＝μ_r；

(3b)根据仿真区域的编号得到仿真区域边界信息，通过积分方程方法得到在上方和下方人工边界处所满足的积分边界条件：

$\mathrm{\Φ}\left(r^{\′}\right)={\mathrm{\Φ}}^{inc}\left(r^{\′}\right)+{\∫}_{{\mathrm{\Γ}}_{UBIE}^{-}}\[\mathrm{\Φ}\left(t\right)\frac{\∂G_U(r,r^{\′})}{\∂n}-G_U(r,r^{\′})\frac{\∂}{\∂n}\mathrm{\Φ}\left(r\right)\]d\mathrm{\Γ},$

$\mathrm{\Φ}\left(r^{\′}\right)={\∫}_{{\mathrm{\Γ}}_{DBIE}^{+}}\[\mathrm{\Φ}\left(r\right)\frac{\∂G_D(r,r^{\′})}{\∂n}-G_D(r,r^{\′})\frac{\∂}{\∂n}\mathrm{\Φ}\left(r\right)\]d\mathrm{\Γ},$

其中，G_U(r,r′)为上方自由空间中的格林函数，G_D(r,r′)为下方介质空间中的格林函数；

(4)将电磁耦合方程组转化为电磁耦合矩阵，并求解电磁耦合矩阵获得介质粗糙面仿真区域的雷达散射参数。

2.根据权利要求1所述的基于边界积分的介质粗糙面有限元电磁仿真方法，其中步骤(1)所述的根据介质粗糙面的功率谱密度和粗糙面参数，通过蒙特卡洛方法产生仿真的介质粗糙面，按如下步骤进行：

(1a)通过实验获得所要仿真介质粗糙面的均方根高度、相关长度粗糙面参数和功率谱密度，并选取一系列振幅独立的高斯随机变量的高斯谐波；

(1b)根据获得的介质粗糙面的均方根高度、相关长度粗糙面参数和功率谱密度，通过蒙特卡洛方法对高斯谐波进行傅里叶变换产生仿真的介质粗糙面。

3.根据权利要求1所述的基于边界积分的介质粗糙面有限元电磁仿真方法，其中步骤(2)所述的通过人工边界将介质粗糙面仿真区域进行截断处理，是在介质粗糙面仿真区域上方和下方分别设置人工截断边界，将介质粗糙面仿真区域划分为内部区域和外部区域，并在仿真区域上方人工截断边界上引入锥形入射波。

4.根据权利要求1所述的基于边界积分的介质粗糙面有限元电磁仿真方法，其中步骤(2)所述的对介质粗糙面仿真区域进行三角面元离散，是将产生的介质粗糙面导入模型处理有限元商用软件ANSYS，剖分介质粗糙面仿真区域，并导出介质粗糙面仿真区域数据文件，获得介质粗糙面三角面元离散仿真区域。

5.根据权利要求1所述的基于边界积分的介质粗糙面有限元电磁仿真方法，其中步骤(2)所述的对离散三角面元进行编号，是将每个三角面元的顶点按右手原则进行编号，并引入整数数组n(e,i)对介质粗糙面仿真区域的离散三角面元编号进行记录，在联系数组n(e,i)中，e是单元编码，i是结点的局部编码，n(e,i)的值是结点的全局编码。

6.根据权利要求1所述的基于边界积分的介质粗糙面有限元电磁仿真方法，其中步骤(4)所述的将电磁耦合方程组转化为电磁耦合矩阵，是通过基函数离散电磁耦合方程组，得到电磁耦合矩阵，并根据LU矩阵求解方法求解电磁耦合矩阵，得出介质粗糙面的雷达散射参数σ(θ)：

$\mathrm{\σ}\left(\mathrm{\θ}\right)=10\lg \left(2\mathrm{\π}r\right|\frac{{\mathrm{\Φ}}^{scat}}{{\mathrm{\Φ}}^{inc}}{|}^2),$

其中，θ是雷达的散射角度，r是仿真位置离原点的距离，Φ^scat是散射场，Φ^inc表示入射场。