CN106294920A - 分析介质目标电磁散射特性的时域准显式方法 - Google Patents

Abstract

本发明公开了一种分析介质目标电磁散射特性的时域准显式方法。步骤如下：建立时域体积分方程；将待求未知量进行时间上用时间基函数进行离散，空间上用非共形的四面体网格进行离散；形成待求解的矩阵方程，未知电流为介质瞬态体电流；利用准显式方法高效求解矩阵方程，得到介质的瞬态体电流系数，再根据体电流系数确定瞬态电磁散射参量。本发明中利用准显式方法求解的非共形时域体积分方程方法可以灵活的处理待求目标的网格离散，形成的相对稠密的迭代矩阵可以利用准显式方法的优势加速求解。

Description

分析介质目标电磁散射特性的时域准显式方法

一技术领域

本发明属于电磁仿真技术领域，特别是一种分析介质目标电磁散射特性的时域准显式方法。

二背景技术

介质目标的电磁散射特性分析是电磁问题中的非常重要研究领域，雷达目标的形状和体积等物理量都是通过对雷达散射截面等参数进行计算得出的，且雷达散射截面积是雷达***对目标“可观测性”的一个重要指标。因此，对于各种目标散射特性的研究在这些应用领域具有特别重要的现实意义。

近年，随着宽频带电磁散射***的快速发展，瞬态电磁散射特性的分析越来越引起科研学者和工程人员的关注。相比于其它方法，时域体积分方程方法非常适合于介质目标瞬态电磁散射特性的分析,尤其适合非均匀介质目标瞬态电磁散射特性的分析(NoelT.Gres,Arif A.Ergin and Eric Michielssen,“Volume-integral-equation-based analysis oftransient electromagnetic scattering from three-dimensional inhomogeneous dielectricobjects,”Radio Science,vol.36,no.3,pp.379–386,2001.)。但是当分析的介质目标存在高度不均匀介电常数，或者是存在多尺度的时候，网格的处理成为了普通时域体积分方程面临的难题，虽然可以利用非共形网格的时域积分方程方法进行分析，但是由于迭代矩阵相对稠密，求解效率较低。

三发明内容

本发明的目的在于提供一种更加高效和准确地分析非均匀或存在多尺度的介质目标电磁散射特性的时域准显式方法。

实现本发明目的的技术解决方案为：一种分析介质目标电磁散射特性的时域准显式方法，步骤如下：

步骤1，建立介质时域体积分方程；

步骤2，将待求未知量采用四阶拉格朗日基函数进行时间上的离散，并采用非共形的四面体单元进行空间上的离散；

步骤3，形成准显式方法待求解的矩阵方程，未知电流为介质瞬态体电流；

步骤4，利用预测-校正方法求解矩阵方程，得到介质的瞬态体电流系数，再根据体电流系数确定瞬态电磁散射参量。

本发明与现有技术相比，其显著优点是：(1)可以更加灵活和准确地分析非均匀或存在多尺度的介质目标瞬态电磁散射特性，并且对离散网格具有鲁棒性；(2)可以高效地求解矩阵方程。

四附图说明

图1是本发明中算例模型示意图。

图2是本发明实施例中介质目标在不同频率点处的双站雷达散射截面结果图。

图3是本发明实施例中介质目标上某一点处的电流随时间的变化曲线图。

五具体实施方式

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

本发明分析介质目标电磁散射特性的时域准显式方法，步骤如下：

步骤1，建立介质时域体积分方程；

令电磁波照射到介质目标上，在介质体内产生感应体电流J，根据介质的电场边界条件，即总电场等于入射电场与散射电场之和，得到介质时域体积分方程TD-VIE，如下：

E^inc(r,t)+E^sca(r,t)＝E^tot(r,t) (1)

其中，E^inc表示照射在介质目标上的电磁波的入射电场，E^tot表示总电场，E^sca表示介质目标在电磁波照射后产生的散射电场，t为当前时刻，瞬态散射电场的表达形式为：

$E^{sca}(r,t)={\∫}_{V^{\′}}\▿\▿\·\frac{{\∂}_t^{-1}J(r^{\′},t-|r-r^{\′}|/c)}{4\mathrm{\π}\mathrm{\ϵ}|r-r^{\′}|}{\mathrm{dV}}^{\′}-{\∫}_{V^{\′}}\frac{\mathrm{\μ}{\∂}_{t}J(r^{\′},t-|r-r^{\′}|/c)}{4\mathrm{\π}|r-r^{\′}|}{\mathrm{dV}}^{\′}---\left(2\right)$

将式(2)代入式(1)，则式(1)重新改写为：

$E^{inc}(r,t)=\frac{{\∂}_t^{-1}J(r^{\′},t)}{{\mathrm{\ϵ}}_0({\mathrm{\ϵ}}_r-1)}+\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫}_{V^{\′}}\frac{{\∂}_{\mathrm{\τ}}J(r^{\′},t-R/c)}{R}{\mathrm{dV}}^{\′}-\frac{1}{4{\mathrm{\π\ϵ}}_0}{\∫}_{V^{\′}}\▿\▿\·\frac{{\∂}_t^{-1}J(r^{\′},t-R/c)}{R}{\mathrm{dV}}^{\′}---\left(3\right)$

其中，V表示四面体单元，μ₀表示自由空间的磁导率，ε表示介电常数，ε₀表示自由空间的介电常数，ε_r为介质体的相对介电常数，r为场的位置坐标，r′为源的位置坐标，R为场源点之间的距离，c表示真空中的光速，表示对时间函数的积分，表示对时间函数的求导，▽为梯度算子。

步骤2，将待求未知量采用四阶拉格朗日基函数进行时间上的离散，并采用非共形的四面体单元进行空间上的离散，具体如下：

介质目标的瞬态感应体电流可离散表示如下：

$J(r,t)=\underset{n=1}{\overset{N_V}{\Σ}}\underset{l=1}{\overset{N_t}{\Σ}}{I}_n^l{T}_l\left(t\right)f_n\left(r\right)---\left(4\right)$

其中：

式中，f_n(r)为半个SWG基函数，T_l(t)为四阶拉格朗日时间基函数，为第n个未知量在第l时刻的待求瞬态电流系数，n为空间未知量编号、l为时间步数编号，N_V为空间未知量个数、N_t为时间步数。

步骤3，形成准显式方法待求解的矩阵方程，未知电流为介质瞬态体电流，具体如下：

为了得到一阶常微分方程的形式，将式(3)在时间上进行两次求导：

${\∂}_t^2{E}^{inc}(r,t)=\frac{{\∂}_{t}J(r^{\′},t)}{{\mathrm{\ϵ}}_0({\mathrm{\ϵ}}_r-1)}+\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫}_{V^{\′}}\frac{{\∂}_t^{3}J(r^{\′},t-R/c)}{R}{\mathrm{dV}}^{\′}-\frac{1}{4{\mathrm{\π\ϵ}}_0}{\∫}_{V^{\′}}\▿\▿\·\frac{{\∂}_{t}J(r^{\′},t-R/c)}{R}{\mathrm{dV}}^{\′}---\left(6\right)$

将式(6)空间上采用伽辽金测试，时间上采用点测试，得准显式的时域体积分的矩阵方程形式：

$G{\∂}_t{I}_i=V_i^{\exp }-\underset{j=0}{\overset{i}{\Σ}}{Z}_{i-j}^{\exp }{I}_j---\left(7\right)$

其中

$V_i^{\exp }={\∫}_V{f}^v\left(r\right)\·{\∂}_t^2{E}^{inc}(r,i\mathrm{\Δ}t)dV---\left(8\right)$

${\[G\]}_{mn}=\frac{1}{{\mathrm{\ϵ}}_0\left({\mathrm{\ϵ}}_r\right(r)-1)}{\∫}_V{f}_m^v\left(r\right)\·f_n^v\left(r^{\′}\right)\mathrm{dV}---\left(9\right)$

式中，

$g_j(i\mathrm{\Δ}t,R)=\frac{T_j(i\mathrm{\Δ}t-R/c)}{R}---\left(11\right)$

其中，为建立场源之间联系的时域阻抗矩阵，G为稀疏的Gram矩阵，表示积分主值项上场源基函数作用形成的矩阵，表示第i个时间步的激励，Δt表示每个时间步长，I_i和I_j分别是第i和第j个时间步的待求未知量的系数，i、j为时间步的编号，f^v(r)为半个SWG基函数，v表示介质，和分别表示第m个和第n个空间基函数，S为积分区域，和分别表示第m个和第n个介质三角形的外法向量，g_j(iΔt,R)为时域格林函数、T_j为时间基函数。

步骤4，利用预测-校正方法求解矩阵方程，得到介质的瞬态体电流系数，再根据体电流系数确定瞬态电磁散射参量，采用预测-修正型线性k步法，具体如下：

在假设式(7)矩阵方程中的I_i(i＝0:j)已知的情况下可以利用式(7)来计算 的值。为了得到I_j的值，应该对式(7)进行时间积分。为此，要用到预测-修正型线性k步法，并且用逐次超松弛(SOR)来提高预测-修正的精确性和稳定性。生成的时间步进如下：

假定在第j时间步，j＝k,…,N_t-1，设置初始值：

$I_i={\stackrel{\·}{I}}_i=0,i=0,\mathrm{...},k-1;$

(4.1)计算式(7)中等号右边在一个时间步长内不会发生变化的固定部分

$V_j^{fixed}=V_j^{\exp }+\underset{i=0}{\overset{j-1}{\Σ}}{Z}_{j-i}^{\exp }{I}_i$

(4.2)通过预测系数p、I_j与的过去值预测I_j：

$I_j^{\left(0\right)}=\underset{l=1}{\overset{k}{\Σ}}\[{\left\{p\right\}}_l{I}_{j-1+l-k}+{\left\{p\right\}}_{k+l}{\tilde{I}}_{i-1+l-k}\]$

上式中，l、k表示时刻的变量；

(4.3)将预测的I_j值及式(7)右边的固定部分代入式(7)计算

$G{\stackrel{\·}{I}}_j^{\left(0\right)}=V_j^{fixed}+Z_0^{\exp }{I}_j^{\left(0\right)}$

其中，表示第j时刻电流导数的系数在预测后的初始值，表示场源基函数距离在一个时间步内的基函数对形成的矩阵；

(4.4)重复操作以下步骤，直到结果I_j收敛，具体如下：

(a)用修正系数c修正I_j：

$I_j^{\left(q\right)}=\underset{l=1}{\overset{k}{\Σ}}\[{\left\{c\right\}}_l{I}_{j-1+l-k}+{\left\{c\right\}}_{k+l}{\tilde{I}}_{j-1+l-k}\]+{\left\{c\right\}}_{2k+1}{\stackrel{\·}{I}}_j^{(q-1)}$

(b)将修正的值与式(7)右边的固定部分代入式(7)计算

$G{\stackrel{\·}{I}}_j^{\left(q\right)}=Z_0^{\exp }{I}_j^{\left(q\right)}+V_j^{fixed}$

其中，表示第j时刻电流导数的系数在第q次修正后的值；

(c)对应用逐次超松弛SOR：

$I_j^{\left(q\right)}={\mathrm{\αI}}_j^{\left(q\right)}+(1-\mathrm{\α})I_j^{(q-1)}$

(d)检查收敛性，如果结果收敛，即χ^PECE为预设的修正收敛精度，那么 $I_j=I_j^{\left(q\right)},{\stackrel{\·}{I}}_j={\stackrel{\·}{I}}_j^{\left(q\right)};$

在步骤(c)中，α∈[0,1]是SOR参数，选择α＝1就是不使用SOR；选择α＜0.5将对更依赖于前面时刻的系数；选择α＞0.5将对更依赖于当前时刻的系数。

通过上述方法得到介质的瞬态体电流系数，再根据体电流系数确定瞬态电磁散射参量。

实施例1

为了验证本发明方法的准确性与高效性，下面给出了一个介质球壳的介质目标的瞬态电磁特性的分析，其中，球壳的内半径为0.25m，厚度为0.05m,相对介电常数为2，如图1。瞬态电磁散射的双站RCS的计算结果与解析值Mie计算的结果相比较吻合较好，图2所示，并验证了在(-0.239,-0.091,0.067)位置处的电流稳定性，如图3。

本算例中，入射电场采用调制高斯平面波，其表达式如下：

$E^{inc}(r,t)={\hat{P}}^{inc}\exp \[-{\left(\frac{\mathrm{\τ}-t_c}{\sqrt{2}\mathrm{\σ}}\right)}^2\]cos\left(2{\mathrm{\πf}}_c\mathrm{\τ}\right)---\left(12\right)$

其中，极化方向 ${\hat{P}}^{inc}=\hat{x},$ 传播方向 ${\hat{k}}^{inc}=\hat{z},\mathrm{\σ}=6/\left(2{\mathrm{\πf}}_{bw}\right),$ t_c＝7σ, $\mathrm{\τ}=t-r\·{\hat{k}}^{inc}/c,$ E^inc(r,t)的频谱的中心频率为f₀＝50MHz，最高频率为100MHz，f_bw为频带宽度。时间步长Δt＝0.3lm，总时间步N_t＝300，lm是光米(light meter)，即光在自由空间中传播1m距离所花的时间。表1给出了求解时间的数据对比，体现出本方法的高效性。

表1计算效率对比

	求解时间(分钟)
		隐式方法	146.87
准显式方法	28.88

综上所述，本发明与传统的时域体积分方程方法相比，基于非共形网格方法求解的时域体积分方程方法可以更加灵活的处理待求目标的网格离散，特别是对于不均匀的介质体或者是存在多尺度模型的情况。并针对稠密矩阵迭代较慢的情况，引入了准显式的求解方法，加速了***的求解速度。

Claims (5)

1.一种分析介质目标电磁散射特性的时域准显式方法，其特征在于，步骤如下：

步骤1，建立介质时域体积分方程；

步骤2，将待求未知量采用四阶拉格朗日基函数进行时间上的离散，并采用非共形的四面体单元进行空间上的离散；

步骤3，形成准显式方法待求解的矩阵方程，未知电流为介质瞬态体电流；

步骤4，利用预测-校正方法求解矩阵方程，得到介质的瞬态体电流系数，再根据体电流系数确定瞬态电磁散射参量。

2.根据权利要求1所述的分析介质目标电磁散射特性的时域准显式方法，其特征在于，步骤1中所述建立介质时域体积分方程，具体如下：

令电磁波照射到介质目标上，在介质体内产生感应体电流J，根据介质的电场边界条件，即总电场等于入射电场与散射电场之和，得到介质时域体积分方程TD-VIE，如下：

E^inc(r,t)+E^sca(r,t)＝E^tot(r,t) (1)

其中，E^inc表示照射在介质目标上的电磁波的入射电场，E^tot表示总电场，E^sca表示介质目标在电磁波照射后产生的散射电场，t为当前时刻，瞬态散射电场的表达形式为：

$E^{sca}(r,t)={\∫}_{V^{\′}}\▿\▿\·\frac{{\∂}_t^{-1}J(r^{\′},t-|r-r^{\′}|/c)}{4\mathrm{\π}\mathrm{\ϵ}|r-r^{\′}|}{\mathrm{dV}}^{\′}-{\∫}_{V^{\′}}\frac{\mathrm{\μ}{\∂}_{t}J(r^{\′},t-|r-r^{\′}|/c)}{4\mathrm{\π}|r-r^{\′}|}{\mathrm{dV}}^{\′}---\left(2\right)$

将式(2)代入式(1)，则式(1)重新改写为：

$E^{inc}(r,t)=\frac{{\∂}_t^{-1}J(r^{\′},t)}{{\mathrm{\ϵ}}_0({\mathrm{\ϵ}}_r-1)}+\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫}_{V^{\′}}\frac{{\∂}_{\mathrm{\τ}}J(r^{\′},t-R/c)}{R}{\mathrm{dV}}^{\′}-\frac{1}{4{\mathrm{\π\ϵ}}_0}{\∫}_{V^{\′}}\▿\▿\·\frac{{\∂}_t^{-1}J(r^{\′},t-R/c)}{R}{\mathrm{dV}}^{\′}---\left(3\right)$

其中，V表示四面体单元，μ₀表示自由空间的磁导率，ε表示介电常数，ε₀表示自由空间的介电常数，ε_r为介质体的相对介电常数，r为场的位置坐标，r′为源的位置坐标，R为场源点之间的距离，c表示真空中的光速，表示对时间函数的积分，表示对时间函数的求导，▽为梯度算子。

3.根据权利要求1所述的分析介质目标电磁散射特性的时域准显式方法，其特征在于，步骤2中所述将待求未知量采用四阶拉格朗日基函数进行时间上的离散，并采用非共形的四面体单元进行空间上的离散，具体如下：

介质目标的瞬态感应体电流离散表示如下：

$J(r,t)=\underset{n=1}{\overset{N_V}{\Σ}}\underset{l=1}{\overset{N_t}{\Σ}}{I}_n^l{T}_l\left(t\right)f_n\left(r\right)---\left(4\right)$

其中：

式中，f_n(r)为半个SWG基函数，T_l(t)为四阶拉格朗日时间基函数，为第n个未知量在第l时刻的待求瞬态电流系数，n为空间未知量编号、l为时间步数编号，N_V为空间未知量个数、N_t为时间步数。

4.根据权利要求1所述的分析介质目标电磁散射特性的时域准显式方法，其特征在于，步骤3中所述形成准显式方法待求解的矩阵方程，未知电流为介质瞬态体电流，具体如下：

为了得到一阶常微分方程的形式，将式(3)在时间上进行两次求导：

${\∂}_t^2{E}^{inc}(r,t)=\frac{{\∂}_{t}J(r^{\′},t)}{{\mathrm{\ϵ}}_0({\mathrm{\ϵ}}_r-1)}+\frac{{\mathrm{\μ}}_0}{4\mathrm{\π}}{\∫}_{V^{\′}}\frac{{\∂}_t^{3}J(r^{\′},t-R/c)}{R}{\mathrm{dV}}^{\′}-\frac{1}{4{\mathrm{\π\ϵ}}_0}{\∫}_{V^{\′}}\▿\▿\·\frac{{\∂}_{t}J(r^{\′},t-R/c)}{R}{\mathrm{dV}}^{\′}---\left(6\right)$

将式(6)空间上采用伽辽金测试，时间上采用点测试，得准显式的时域体积分的矩阵方程形式：

$G{\∂}_t{I}_i=V_i^{\exp }-\underset{j=0}{\overset{i}{\Σ}}{Z}_{i-j}^{\exp }{I}_j---\left(7\right)$

式中，

$V_i^{exp}={\∫}_V{f}^v\left(r\right)\·{\∂}_t^2{E}^{inc}(r,i\mathrm{\Δ}t)dV---\left(8\right)$

${\[G\]}_{mn}=\frac{1}{{\mathrm{\ϵ}}_0\left({\mathrm{\ϵ}}_r\right(r)-1)}{\∫}_V{f}_m^v\left(r\right)\·f_n^v\left(r^{\′}\right)dV---\left(9\right)$

式中，

$g_j(i\mathrm{\Δ}t,R)=\frac{T_j(i\mathrm{\Δ}t-R/c)}{R}---\left(11\right)$

其中，为建立场源之间联系的时域阻抗矩阵，G为稀疏的Gram矩阵，表示第i个时间步的激励，Δt表示每个时间步长，I_i和I_j分别是第i和第j个时间步的待求未知量的系数，i、j为时间步的编号，f^v(r)为半个SWG基函数_，v表示介质，和分别表示第m个和第n个空间基函数，S为积分区域，和分别表示第m个和第n个介质三角形的外法向量，g_j(iΔt,R)为时域格林函数、T_j为时间基函数。

5.根据权利要求1所述的分析介质目标电磁散射特性的时域准显式方法，其特征在于，步骤4中所述利用预测-校正方法求解矩阵方程，得到介质的瞬态体电流系数，再根据体电流系数确定瞬态电磁散射参量，采用预测-修正型线性k步法，具体如下：

假定在第j时间步，j＝k,…,N_t-1，设置初始值：

$I_i={\stackrel{\·}{I}}_i=0,i=0,...,k-1;$

(4.1)计算式(7)中等号右边在一个时间步长内不会发生变化的固定部分

$V_j^{fixed}=V_j^{exp}+\underset{i=0}{\overset{j-1}{\Σ}}{Z}_{j-i}^{exp}{I}_i$

(4.2)通过预测系数p、I_j与的过去值预测I_j：

$I_j^{\left(0\right)}=\underset{l=1}{\overset{k}{\Σ}}\[{\left\{p\right\}}_l{I}_{j-1+l-k}+{\left\{p\right\}}_{k+l}{\tilde{I}}_{j-1+l-k}\]$

上式中，l、k表示时刻的变量；

(4.3)将预测的I_j值及式(7)右边的固定部分代入式(7)计算

$G{\stackrel{\·}{I}}_j^{\left(0\right)}=V_j^{fixed}+Z_0^{\exp }{I}_j^{\left(0\right)}$

其中，表示第j时刻电流导数的系数在预测后的初始值，表示场源基函数距离在一个时间步内的基函数对形成的矩阵；

(4.4)重复操作以下步骤，直到结果I_j收敛，具体如下：

(a)用修正系数c修正I_j：

$I_j^{\left(q\right)}=\underset{l=1}{\overset{k}{\Σ}}\[{\left\{c\right\}}_l{I}_{j-1+l-k}+{\left\{c\right\}}_{k+l}{\tilde{I}}_{j-1+l-k}\]+{\left\{c\right\}}_{2k+1}{\stackrel{\·}{I}}_j^{(q-1)}$

(b)将修正的值与式(7)右边的固定部分代入式(7)计算

$G{\stackrel{\·}{I}}_j^{\left(q\right)}=Z_0^{\exp }{I}_j^{\left(q\right)}+V_j^{fixed}$

其中，表示第j时刻电流导数的系数在第q次修正后的值；

(c)对应用逐次超松弛SOR：

$I_j^{\left(q\right)}={\mathrm{\αI}}_j^{\left(q\right)}+(1-\mathrm{\α})I_j^{(q-1)}$

(d)检查收敛性，如果结果收敛，即χ^PECE为预设的修正收敛精度，那么 $I_j=I_j^{\left(q\right)},{\stackrel{\·}{I}}_j={\stackrel{\·}{I}}_j^{\left(q\right)};$

在步骤(c)中，α∈[0,1]是SOR参数，选择α＝1就是不使用SOR；选择α＜0.5将对更依赖于前面时刻的系数；选择α＞0.5将对更依赖于当前时刻的系数。