一种基于保散度ADI-FDTD的探地雷达数值求解方法及系统

本发明涉及探地雷达中电磁波的高效数值模拟领域，具体涉及一种基于保散度adi-fdtd的探地雷达数值求解方法及系统。

背景技术：

1、探地雷达（ground penetrating radar，gpr）是一种广泛使用的地球物理技术，它利用天线发射和接收电磁（em）波，能够探测地下介质中的物质属性和分布模式。当em波穿过地下时，其路径、场强、和波形受介质的电特性和几何特性的影响，通过分析传播时间、幅度、波形，和接收信号的其他参数（通常显示在示波器上）并将它们与预先计算的模型进行比较，可以推断地下界面或地质构造的空间位置、结构和分布。通过em模拟生成与地下条件相对应的反射波信号数据库，便于正演建模。麦克斯韦方程控制em场通过这些介质的传播，而求解这些方程是精确地进行探地雷达模拟的关键。

2、fdtd方法在时间和空间域中离散麦克斯韦方程的广泛采用的方法。通过这样做，它使得能够沿时间轴沿着计算空间em场，有效地捕获瞬态em场信息。该方法利用空间网格将目标介质的电磁特性纳入时域有限差分法中，通过在不同区域分配适当的电磁参数，可以模拟包括非均匀性、各向异性、色散特性和非线性的各种复杂结构。

3、尽管其实用性，传统的显式fdtd方法有一个关键的限制：时间步长受到courant-friedrichs-lewy（cfl）稳定性条件的限制。在gpr模拟中，土壤和目标的色散特性，结合某些区域的不规则结构特征，需要使用精细，均匀的yee网格以获得数值精度。 然而，这导致了两个主要挑战：

4、1）时间步长必须足够小以满足cfl稳定性条件，导致计算时间较长；

5、2）均匀细网格离散化增加了未知数的数量，显著增加了内存使用量。

6、为了解决这些挑战，几种无条件稳定的fdtd方法，如adi-fdtd和lod-fdtd方法，dp-adi-fdtd方法是一种保发散、无条件稳定的方法，它在求解麦克斯韦方程组时，会失去保发散特性，导致无源区电荷积累，被引入以解决em粒子在细胞（empic）问题。 该方法后来被推广用于模拟多极德拜色散介质，尽管其优点，但基于ade方法的数值实现仍然是复杂和繁琐的。

7、此外，还提出了子网格技术，通过消除整个计算域中对统一yee网格的需求，来提高gpr和其他em模拟的计算效率。这些方法允许在需要高分辨率的区域进行更精细的离散化，而在不太关键的区域使用更粗糙的网格，从而防止不必要的过采样。以前的子网格方法可以分为两种类型：

8、1）一种方法采用满足整个域上的细网格的稳定性条件的时间步长，其中空间插值用于在粗网格区域和细网格区域之间交换em场信息；

9、2）另一种方法对粗网格区域和细网格区域利用不同的时间步长，每个时间步长满足其各自的cfl稳定性条件，采用时空插值在网格之间交换em信息。

10、最近，无条件稳定的fdtd方法已应用于细网格区域，使得整个域能够使用满足粗网格cfl稳定性条件的时间步长来求解。

11、这些子网格技术允许更大的时间步长，与传统的混合亚网格方法相比，该方法大大提高了计算效率，而且采用统一的时间步长，消除了时间插值的需要，简化了数值实现，然而，现有的方法仍然面临着一些挑战，包括：

12、1）主要关注二维模型，对复杂色散介质的适用性有限；

13、2）缺乏统一的，任意比例的电磁交换公式在粗-细网格界；

14、3）缺乏的无条件稳定的fdtd方法中使用的混合子网格方法的发散保持属性。

技术实现思路

1、针对现有技术的不足，本发明提出了一种基于保散度adi-fdtd的探地雷达数值求解方法及系统。

2、本发明的目的可以通过以下技术方案实现：

3、本发明的第一方面，涉及一种基于保散度adi-fdtd的探地雷达数值求解方法，包括以下步骤：

4、列出探地雷达的麦克斯韦方程组的微分时域；

5、将所述麦克斯韦方程组的微分时域离散成矩阵形式；

6、对于具有损耗项的多极德拜色散介质，构建电通量密度与电场强度的本构关系式；

7、根据z变换理论，将所述本构关系式从频域变换到z域，再通过z域与时域之间的变换，转换到时域，并代入到所述麦克斯韦方程组的微分时域离散成的矩阵形式中，得到基于时域的离散矩阵形式；

8、将所述基于时域的离散矩阵形式离散为两个子时间步长迭代，第一个子时间段，第二个子时间段，推导电磁场分量的数值迭代公式，n为迭代次数；分别计算第一个子时间段与第二个子时间段的电场和磁场强度表达式，获得有耗多极德拜模型zt-dp-adi-fdtd方法的基本数值离散化框架；

9、根据所述基本数值离散化框架计算电场与磁场各分量的数值离散化方程，并计算电场强度与磁场强度，推导出满足保散度的电场强度和磁场强度。

10、可选地，所述麦克斯韦方程组的微分时域表达式如下所示：

11、 (1)

12、 (2)

13、其中d为电通量密度，h为磁场强度，e为电场强度，μ0为自由空间磁导率，电通量密度d与电场e的本构关系表示为d = ε0εre， ε0为自由空间介电常数，εr为相对介电常数。

14、可选地，所述矩阵形式为：

15、 (3)

16、其中vdp= [edp_x, edp_y, edp_z, hdp_x,hdp_y, hdp_z]表示保散度的电场和磁场强度分量，vdp的上标n与n+1表示迭代步数；i是单位矩阵；edp_x, edp_y, edp_z分别为电场强度在x、y、z方向的强度，hdp_x,hdp_y, hdp_z分别为磁场强度在x、y、z方向的强度，矩阵a和b为空间偏导数矩阵，具体为：

17、,，其中,，，和分别表示场分量在x、y和z方向的偏导数。

18、可选地，所述构建电通量密度与电场强度的本构关系式：

19、 (4)

20、其中是电导率，是零频率和无限频率的相对介电系数之差，表示无穷大频率的相对介电系数；是极点弛豫时间，是灵敏度响应中的极数；表示虚部；表示角频率；表示频域下的电通量密度；表示频域下的电场强度。

21、可选地，所述的将所述本构关系式从频域变换到z域，再通过z域与时域之间的变换，转换到时域，包括以下步骤：

22、根据z变换理论，将(4)从频域变换到z域，得到：

23、 (5)

24、其中ψp和q为辅助变量，定义如下：

25、 (6)

26、 (7)

27、通过z域与时域之间的变换，将公式(5)、公式(6)与公式(7)转换到时域，代入(3)，得到：

28、 (8)

29、其中φp=ψp/ε0ε∞,= q/ε0ε∞，此外矩阵中a和b中的εr变为ε∞，是零频率和无限频率的相对介电系数之差；：z域下的ψp，：z域下的q。

30、可选地，所述有耗多极德拜模型zt-dp-adi-fdtd方法的基本数值离散化框架的计算过程包括以下步骤：

31、将公式(8)离散为两个子时间步长迭代，对于第一个子时间段：

32、 (9)

33、 (10)

34、对于第二个子时间段：

35、

36、其中 v = [ex, ey, ez, hx, hy, hz]， ex、ey、ez分别表示不服从散度特性的电场在x、y、z方向的分量，hx、hy、hz分别表示不服从散度特性的磁场在x、y、z方向的分量；、φp、vdp、的上标均表示迭代步数；

37、将（10）代入（11）得：

38、 (13)

39、定义辅助变量u = [ueuh] = [uex, uey, uez, uhx, uhy, uhz]，ue、uh分别表示电场和磁场的辅助变量，uex、uey、uez分别表示所述电场的辅助变量在x、y、z方向的分量，uhx、uhy、uhz分别表示所述磁场的辅助变量在x、y、z方向的分量；其中；则（13）重新表述为：

40、 (14)

41、将(9)前一个时间步并代入公式（12）得到：

42、 (15)

43、令un+1= vn+1+vn+1/2，则公式（15）重新表述为：

44、 (16)。

45、可选地，根据公式（14），获得ue和uh场分量的数值离散化方程：

46、

47、ψp、ue、 uh的上标均表示迭代步数，其中：

48、

49、求电场和磁场的数值迭代公式：

50、 (24a)

51、此外，和的迭代公式为：

52、

53、uex,、uey、uez、ex、ey、ez的上标均表示迭代步数，、、分别表示x、y、z方向上的坐标索引；

54、将条件代入公式（18）得：

55、 (25)

56、得到的计算公式为：

57、 (26)

58、将公式（24a）、公式（24b）、公式（24c）和公式（26）带入到公式（9）中，得到满足保散度的edp和hdp场的计算公式：

59、 (27)

60、h表示磁场，e表示电场，e与h的上标表示迭代步数。

61、本发明的第二方面，涉及一种基于保散度adi-fdtd的探地雷达数值求解系统，包括：

62、麦克斯韦方程组计算模块，用于列出探地雷达的麦克斯韦方程组的微分时域；将所述麦克斯韦方程组的微分时域离散成矩阵形式；

63、电通量密度与电场强度的本构关系式构建模块，对于具有损耗项的多极德拜色散介质，构建电通量密度与电场强度的本构关系式；

64、z域变换模块：用于根据z变换理论，将所述本构关系式从频域变换到z域，再通过z域与时域之间的变换，转换到时域，并代入到所述麦克斯韦方程组的微分时域离散成的矩阵形式中，得到基于时域的离散矩阵形式；

65、基本数值离散化框架构建模块：将所述基于时域的离散矩阵形式离散为两个子时间步长迭代，第一个子时间段，第二个子时间段，推导电磁场分量的数值迭代公式，n为迭代次数；分别计算第一个子时间段与第二个子时间段的电场和磁场强度表达式，获得有耗多极德拜模型zt-dp-adi-fdtd方法的基本数值离散化框架；

66、以及，保散度的电场强度和磁场强度求解模块：根据所述基本数值离散化框架计算电场与磁场各分量的数值离散化方程，并计算电场强度与磁场强度，推导出满足保散度的电场强度和磁场强度。

67、本发明的第三方面，涉及一种基于三维混合子网格的探地雷达数值求解方法，所述三维混合子网格包括若干粗网格与细网格，激励源放置在粗网格区域，当电磁波向粗网格和细网格分界面区域传播时，将分界面处粗网格中的电磁场值插值到相应的细网格分量中；

68、所述方法包括以下步骤：

69、采用时域有限差分法更新整个粗网格电场值和磁场值；

70、通过上述基于保散度adi-fdtd的探地雷达数值求解方法，计算细网格区域内的细网格电场值与磁场值；

71、借助空间线性插值将粗网格中的电场值与磁场值插值到处于同一界面的细网格中界面处；利用上述基于保散度adi-fdtd的探地雷达数值求解方法进行数值迭代，将细网格界面处的电场分量与磁场分量传递至细网格区域；

72、细网格区域内完成迭代后，将细网格电磁场在粗网格与细网格交界处以及相邻交界处的值插值回相应交界处的粗网格电磁场值。

73、本发明的第四方面，涉及一种探地雷达，包括存储有能够运行上述基于保散度adi-fdtd的探地雷达数值求解方法或上述的基于三维混合子网格的探地雷达数值求解方法的指令的存储介质，或者，上述的基于保散度adi-fdtd的探地雷达数值求解系统。

74、本发明的有益效果：

75、1）本发明基于z变换的dp-adi-fdtd方法（zt-dp-adi-fdtd）方法，由于其无条件稳定和保散特性，可以有效和精确地模拟有耗多级德拜色散介质中的电磁波传播，与ade方法相比，z变换技术将多级德拜模型从频域转换到z域，消除了存储来自多个先前时间步的场的需要，同时保留核心迭代方程。这导致更灵活和有效的数值实现。

76、2）本发明将显式fdtd方法与隐式zt-dp-adi-fdtd方法相结合，发展了一种混合亚网格技术，该方法在细网格区域采用zt-dp-adi-fdtd方法计算电磁场分量，在粗网格区域采用fdtd方法，这种混合子网格技术可以适应任意的奇偶网格比，进一步提高了建模过程的灵活性。