一种时域有限差分法计算电磁散射的瞬态场远场外推方法

摘要

一种时域有限差分法计算电磁散射的瞬态场远场外推方法，步骤有：1、通过连接边界将计算区域划分成总场区和散射场区两个区域；2、在散射场区设置外推边界，将外推边界上的电场和磁场外推得到远区场；3、设置瞬态入射波，在每个时间步依次更新电场和磁场，直到电磁场分布趋于稳态；4、将远场写成频率因子和积分因子的乘积，对积分因子进行逆傅里叶变换，得到时域表达式；5、设置积分因子远场时域响应数组；6、积分因子远场时域响应经傅里叶变换得到其频域响应；7、利用远场和入射波的频域响应，得到RCS随频率的变化。该方法将远场写成频率因子和积分因子之积的形式，积分因子的频时变换和时频变换形式较为简单，简化了计算，提高了效率。

权利要求

1.一种时域有限差分法计算电磁散射的瞬态场远场外推方法，其特征在于：该方法包含以下步骤：步骤1：通过连接边界将计算区域划分成总场区和散射场区两个区域，总场区的电磁场包含了入射场和散射场，而散射场区只包含散射场；

步骤2：在散射场区设置外推边界，将外推边界上的电场和磁场外推得到远区场；

步骤3：设置瞬态入射波，在每个时间步依次更新电场和磁场，直到电磁场分布趋于稳态；

步骤4：将远场写成频率因子和积分因子的乘积，对积分因子进行逆傅里叶变换，得到时域表达式；

步骤5：设置积分因子远场时域响应数组，将每个时间步外推边界电磁场对远场的响应叠加到积分因子远场时域响应数组中；

步骤6：积分因子远场时域响应经傅里叶变换得到其频域响应；

步骤7：利用远场和入射波的频域响应，得到RCS随频率的变化；

其中，步骤4所述的将远场写成频率因子和积分因子的乘积，其三维情形表达式为：其中 是雷达波接收天线极化方向，jke-jkr/(4πr)是频率因子，I是积分因子，Es是散射电场强度，积分因子I的表达式为其中 是远场接收天线磁场极化方向，r是源点到场点的距离，S是外推边界面，是真空波阻抗；r′为外推边界上的位置矢量；

其二维情形为横磁波即TM波，表达式为：

其中 是z方向的坐标矢量， 是频率因子，I′是积分因子，表达式为其中l为外推边界，二维情形下外推边界是一个矩形，因此具有线积分的形式。

2.根据权利要求1所述的一种时域有限差分法计算电磁散射的瞬态场远场外推方法，其特征在于：在步骤1中所述的“计算区域”，是指一个三维情形的长方体，或是二维情形的矩形区域，电磁场的采样或迭代就在此区域进行；所述的“划分成总场区和散射场区两个区域”，是指用一个三维情形的长方体，或是二维情形矩形的连接边界，将计算区域一分为二，连接边界以内是总场区，连接边界以外是散射场区，散射体位于总场区。

3.根据权利要求1所述的一种时域有限差分法计算电磁散射的瞬态场远场外推方法，其特征在于：在步骤2中所述的“外推边界”，是指位于散射场区的一个三维情形的长方体，或是二维情形矩形的边界，将总场区包含在内，由于外推边界位于散射场区，因此外推边界上的电磁场只有散射场，不包含入射场，由其外推得到的远场就只有散射场的贡献，直接计算雷达散射截面即RCS。

4.根据权利要求1所述的一种时域有限差分法计算电磁散射的瞬态场远场外推方法，其特征在于：在步骤3中所述的“电磁场分布趋于稳态”，是指总场区和散射场区的电磁场上不发生变化。

5.根据权利要求1所述的一种时域有限差分法计算电磁散射的瞬态场远场外推方法，其特征在于：在步骤5中所述的“将每个时间步外推边界电磁场对远场的响应叠加到积分因子远场时域响应数组中”，是指每个时间步外推边界上不同位置对积分因子的远场贡献体现在不同的时刻，因此需要预先准备好待填充的积分因子远场时域响应数组，在每个时间步将对应的远场贡献叠加到积分因子远场时域响应数组中。

6.根据权利要求1所述的一种时域有限差分法计算电磁散射的瞬态场远场外推方法，其特征在于：在步骤7中所述的“利用远场和入射波的频域响应，得到RCS随频率的变化”，是指根据RCS的定义，将RCS写成如下形式：对应的三维情形：公式(5)

对应的三维情形：公式(6)

其中f是频率，I和I′是三维和二维情形的积分因子远场频域响应，c是光速，Ei是入射波。

说明书

一种时域有限差分法计算电磁散射的瞬态场远场外推方法

一、技术领域

[0001] 本发明提供一种时域有限差分法计算电磁散射的瞬态场远场外推方法，它具体涉及时域有限差分法(finite-difference time-domain，FDTD)的电磁仿真方法，属于电磁场仿真技术领域。二、背景技术

[0002] 时域有限差分法即FDTD是电磁仿真的一种重要方法，在电磁领域具有广泛应用。该方法的基本原理是将麦克斯韦方程组离散成差分格式，利用Yee氏网格在空间和时间上的交错特性，实现电场和磁场的交替更新，求解得到电磁场随空间和时间分布。在用FDTD计算目标电磁散射时，首先将计算区域划分成总场区和散射场区，通过连接边界加入入射电磁波，根据等效原理将散射场区的电磁场外推至远场，可求得目标的雷达散射截面(radar cross section，RCS)。

[0003] 如果选取时谐形式的入射场，可以在电磁场分布趋于稳态时，求解得出电磁场分布的幅度和相位信息，但用这种方法只能得出一个频率的解，不能完全发挥FDTD方法时域计算的优势。为得到目标电磁散射的频域响应，往往采用瞬态场作为入射场，求解得到远场瞬态响应，再通过傅里叶变换可以得到频域响应。传统方法三维瞬态场远场外推需要进行远场时域响应的差分计算，计算较为复杂。二维瞬态场远场外推还需要对比二维和三维的频域远区场表达式，以三维瞬态场方法为前提导出结果。

[0004] 本发明提出的瞬态场外推方法将远场表示成频率因子与积分因子的乘积，只对积分因子进行变换，频率因子保持不变，避免了差分计算，简化了时域到频域的变换，同时也将二维与三维的外推方法统一起来，二维瞬态场外推也不需要以三维瞬态场外推结果为前提。频率因子不参与频时与时频变换，而积分因子进行远场外推得到时域响应，进一步经过傅里叶变换得到远场频域响应，再与频率因子相乘，最终得到各频率上RCS。三、发明内容

[0005] 本发明的目的是针对FDTD计算中传统瞬态场远场外推方法的形式复杂、计算量大、三维和二维外推方法存在关联的特点，提出一种新的瞬态场远场外推方法，即一种时域有限差分法计算电磁散射的瞬态场远场外推方法，该方法具体包含以下步骤：

[0006] 步骤1：通过连接边界将计算区域划分成总场区和散射场区两个区域，总场区的电磁场包含了入射场和散射场，而散射场区只包含散射场；

[0007] 步骤2：在散射场区设置外推边界，将外推边界上的电场和磁场外推得到远区场；

[0008] 步骤3：设置瞬态入射波，在每个时间步依次更新电场和磁场，直到电磁场分布趋于稳态；

[0009] 步骤4：将远场写成频率因子和积分因子的乘积，对积分因子进行逆傅里叶变换，得到时域表达式；

[0010] 步骤5：设置积分因子远场时域响应数组，将每个时间步外推边界电磁场对远场的响应叠加到积分因子远场时域响应数组中；

[0011] 步骤6：积分因子远场时域响应经傅里叶变换得到其频域响应；

[0012] 步骤7：利用远场和入射波的频域响应，得到RCS随频率的变化。

[0013] 其中，在步骤1中所述的“计算区域”，是指一个长方体(三维情形)或矩形(二维情形)区域，电磁场的采样、迭代就在此区域进行。所属的“划分成总场区和散射场区两个区域”，是指如图2所示，用一个长方体(三维情形)或矩形(二维情形)的连接边界，将计算区域一分为二，连接边界以内是总场区，连接边界以外是散射场区，散射体位于总场区。

[0014] 其中，在步骤2中所述的“外推边界”，是指位于散射场区的一个长方体(三维情形)或矩形(二维情形)边界，将总场区包含在内，由于外推边界位于散射场区，因此外推边界上的电磁场只有散射场，不包含入射场，由其外推得到的远场就只有散射场的贡献，可直接计算雷达散射截面即RCS。

[0015] 其中，在步骤3中所述的“电磁场分布趋于稳态”，是指总场区和散射场区的电磁场基本上不发生变化。

[0016] 其中，在步骤4所述的“将远场写成频率因子和积分因子的乘积”，其三维情形是指如下表达式：

[0017]

[0018] 其中 是雷达波接收天线极化方向，jke-jkr/(4πr)是频率因子，I是积分因子，Ei是入射电场强度，积分因子I的表达式为

[0019]

[0020] 其中 是远场接收天线磁场极化方向，r是源点到场点的距离，S是外推边界面， 是真空波阻抗；

[0021] 其二维情形以横磁(Transverse Magnetic，TM)波为例，表达式如下：

[0022]

[0023] 其中 是z方向的坐标矢量， 是频率因子，I′是积分因子，表达式为

[0024]

[0025] 其中l为外推边界，二维情形下外推边界是一个矩形，因此具有线积分的形式。

[0026] 其中，在步骤5中所述的“将每个时间步外推边界电磁场对远场的响应叠加到积分因子远场时域响应数组中”，是指每个时间步外推边界上不同位置对积分因子的远场贡献体现在不同的时刻，因此需要预先准备好待填充的积分因子远场时域响应数组，在每个时间步将对应的远场贡献叠加到积分因子远场时域响应数组中。

[0027] 其中，在步骤7中所述的“利用远场和入射波的频域响应，得到RCS随频率的变化”，是指根据RCS的定义，将RCS写成如下形式：

[0028] (三维情形)                    (5)

[0029] (二维情形)                     (6)

[0030] 其中f是频率，I和I′是三维和二维情形的积分因子远场频域响应，c是光速，Ei是入射波。

[0031] 本发明所述方法的优点是：该方法将远场写成频率因子和积分因子之积的形式，频率因子不参与变换，而积分因子的频时变换和时频变换形式又较为简单，因此简化了计算，提高了效率。四、附图说明

[0032] 图1是本发明所述方法流程图；

[0033] 图2是本发明计算区域的划分示意图。五、具体实施方式

[0034] 下面将结合附图和实施例对本发明做进一步的详细说明。

[0035] 如图1所示，本发明一种时域有限差分法计算电磁散射的瞬态场远场外推方法，它包括以下步骤：

[0036] 步骤1：通过连接边界将计算区域分成总场区和散射场区，总场区的电磁场包含了入射场和散射场，而散射场区的电磁场只包含散射场；

[0037] 步骤2：在散射场区设置封闭的外推边界，用于远场外推；

[0038] 步骤3：设置瞬态入射波，在每个时间步依次更新电场和磁场，直到电磁场分布趋于稳态；

[0039] 步骤4：将远场写成频率因子和积分因子的乘积，将积分因子进行逆傅里叶变换，得到时域表达式；

[0040] 步骤5：设置远场积分因子时域响应数组，将每个时间步外推边界电磁场对远场的响应叠加到远场积分因子时域响应数组中；

[0041] 步骤6：远场积分因子时域响应经傅里叶变换得到远场积分因子频域响应；

[0042] 步骤7：将频率因子和远场积分因子频域响应乘积，再结合入射场的频域响应，得到RCS随频率的变化。

[0043] 针对上述步骤的进一步说明如下：

[0044] 步骤1：计算场区划分

[0045] 如图2所示，连接边界将计算区域划分成了总场区和散射场区。总场区的电磁场包含了入射场和散射场，而散射场区只包含散射场。总场区与散射场区在连接边界是不连续的，需要在时域差分过程中考虑入射场的影响。而在总场区的内部和散射场区的内部，都各自满足麦克斯韦方程组，不需要对差分格式进行改动。在计算区域的边界，通过吸收边界加以截断，以模拟无限大空间。

[0046] 步骤2：设置外推边界

[0047] 如图2所示，在散射场区设置外推边界，用于远场外推。外推边界位于散射场区，通过等效原理可以由外推边界上的场得到外推边界外的场，如外推边界外的电场Es表达形式为

[0048]

[0049] 其中 是虚单位，ω＝2πf是角频率，ε是介电常数，A是电矢量位函数，F是磁矢量位函数，是电标量位函数，三个位函数表达式分别是

[0050]

[0051] 其中μ是磁导率，是外推边界外法向，E和H分别是外推边界上的电场和磁场，G＝e-jk|r-r′|/(4π|r-r′|)是格林函数，r为场点位置矢量，r′为外推边界上的位置矢量，S是外推边界。设r＝|r|，对于远场有r＞＞1，方向的远场格林函数形式为 进一步通过对应关系 可以将远区电场写成如下形式：

[0052]

[0053] 积分中的第二项只有沿散射方向的分量，没有横向分量，而RCS计算过程中只考虑横向分量，因此第二项在后续计算中不予可虑。

[0054] 步骤3：设置瞬态入射场,进行差分计算

[0055] 瞬态入射波具有一定的带宽，可经一次仿真得出频域响应。本发明采用微分高斯脉冲作为瞬态入射波，表达形式如下式所示：

[0056]

[0057] 其中τ是脉冲宽度，t0是时间偏移量。具体仿真过程中，参数设为t0＝0.8τ，τ＝60Δt，其中Δt是时间步长。该瞬态入射波的频域响应是

[0058]

[0059] 将上式写成频率f的表达式，如下所示：

[0060]

[0061] 在每个时间步依次更新电场和磁场的空间分布，直到电磁场分布趋于稳态。更新方法是将麦克斯韦方程组的微分形式离散成差分形式，再将电场与磁场表示成显式形式，如下所示：

[0062]

[0063] 式中，Hn-1/2和Hn+1/2是时刻(n-1/2)Δt和(n+1/2)Δt的磁场强度，En和En+1是时刻nΔt和(n+1)Δt的电场强度，σ是电导率，σm是导磁率。通过式(13)的差分方程，即可求出电场和磁场在计算区域中的分布情况。

[0064] 步骤4：将远场写成频率因子和积分因子的乘积，并将积分因子进行逆傅里叶变换得到时域形式

[0065] 设 是雷达波接收天线极化方向，则考虑到式(9)，该方向上的远场电场可表示为[0066]

[0067] 其中jke-jkr/(4πr)是频率因子，I是积分因子，Ei是入射电场强度，积分因子I的表达式为

[0068]

[0069] 其中 是远场接收天线磁场极化方向，r是源点到场点的距离，S是外推边界面， 是真空波阻抗。根据k＝ω/c，其中c是光速，积分因子I可写为：

[0070]

[0071] 将上式进行逆傅里叶变换，写成由磁场和电场两部分贡献的和的形式：

[0072] I(t)＝IH(t)+IE(t)                    (17)

[0073] 其中IH(t)表示由切向磁场引起的远区场，IE(t)表示由切向电场引起的远区场，表达式为。

[0074]

[0075]

[0076] 二维情形的远场有类似的形式，以横磁(Transverse Magnetic，TM)波为例，表达式如下所示：

[0077]

[0078] 其中 是z方向的坐标矢量， 是频率因子，I′是二维情形下的积分因子，表达式为

[0079]

[0080] 其中l为外推边界。二维情形下外推边界是一个矩形，因此具有线积分的形式。

[0081] 步骤5：计算积分因子的远场时域响应

[0082] 设置积分因子远场时域响应数组，将每个时间步外推边界电磁场对远场的响应叠加到积分因子远场时域响应数组中。

[0083] 由IH(t)和IE(t)的表达式(18)和(19)可以看出，由于r′在积分的外推边界面上，要得到t时刻的积分因子表达式需要首先求得磁场H和电场E在一段时间内的值，而一般的FDTD计算过程只保留当前时间步的值，使用起来不变。本发明采用的方法是首先设置积分因子的远场响应数组，将每个时间步外推边界上各个网格对远场的贡献叠加到远场响应数组中，这样就不需要保存电场和磁场在一段时间内的值，极大减少了内存开销。将式(18)的积分离散化，如下所示：

[0084]

[0085] 其中IHm(t)表示成外推边界面上每个网格的贡献，将之表示成网格中点处的值与网格面积的乘积，如下所示：

[0086]

[0087] 其中m表示外推边界上网格序号， 表示第m个网格中点处时刻的磁场强度，δ是网格尺寸，则δ2是网格面积。式(23)还可以写成如下形式：

[0088]

[0089] 上式可由当前时间步的磁场得到不同时刻远场贡献。将空间、时间离散化，如下所示：

[0090]

[0091] 其中s＝cΔt/δ是时间因子，nx,ny,nz是r′m在x,y,z方向的坐标到原点的网格数，n是仿真当前时间步数。本发明中电场在n时间步采样，磁场在n+1/2时间步采样，因此式(25)表达式中含有n+1/2。用FDTD中的离散网格序号、时间步长表示式(25)，得到

[0092]

[0093] 其中nHf＝1/2-(sxnx+syny+sznz)/s。式(26)意味着将n+1/2时间步的磁场切向分量赋予n+nHf时间步的远场积分因子，但在程序实现过程中远场积分因子只在整数时间步取值，而n+nHf一般来说不是整数。本发明采用的方法是将该时刻的积分因子按一定比例分配给相邻的两个整数，越临近的整数所占比比重越大。设[nHf]是nHf的整数部分，{nHf}是nHf的小数部分，则可将IHm按照1-{nHf}和{nHf}的比例分配给n+[nHf]和n+[nHf]+1时间步的远场积分因子。考虑到数组的下标都取正整数值，因此再加上充分大的时间步偏移量nτ，满足数组下标为正的要求，得到：

[0094]

[0095]

[0096] 对于外推边界上电场对远场的贡献，需要对式(25)进行修改，两边公式中不出现1/2，因为电场是在整数时间步取值。由时间步n的切向电场求得远区场的积分因子时域响应公式如下所示

[0097]

[0098]

[0099] 其中IEm表示第m个网格上的切向电场对远区的贡献，nEf＝-(sxnx+syny+sznz)/s，[nEf]表示nEf的整数部分，{nEf}表示nEf的小数部分。

[0100] 步骤6：远场积分因子时域响应傅里叶变换

[0101] 远场积分因子时域响应经傅里叶变换得到远场积分因子频域响应，变换公式为[0102]

[0103] 由于I(t)只在整数时间步上取值，因此可以将积分写为

[0104]

[0105] 其中N是积分因子远场的时间步总数。将式(32)写成频率f的形式，可得到积分因子远场频域响应，即

[0106]

[0107] 步骤7：计算RCS随频率的变化

[0108] 要求的RCS，还需要求得入射波的频域响应，计算方法如下：

[0109]

[0110] 式中M是仿真的步数，Ei(n)是第n个时间步上的入射电场强度。再根据式(14)，即可得到RCS随频率的变化：

[0111]

[0112] 上式用到了k＝2πf/c，进一步求得RCS的值

[0113]

[0114] 对于二维情形，RCS有类似的表达式，为

[0115]

[0116] 应当指出，本实例仅列示性说明本发明的应用方法，而非用于限制本发明。任何熟悉此种使用技术的人员，均可在不违背本发明的精神及范围下，对上述实施例进行修改。因此，本发明的权利保护范围，应如权利要求书所列。