探地雷达正演模拟，基于时域有限差分方法，一

     声明：本博客中的公式均是在Word中使用AxMath写好后截图使用的，欢迎引用，但请标注来源。

本系列会有四篇博客：

第一篇内容：

1、基础知识掌握

2、Maxwell方法差分求解原理

第二篇内容：

1、基于C++的TE波波动方程实现

2、边界问题的产生及处理

第三篇内容：

1、基于实际操作流程的GPR正演模拟

2、简单并行化加速GPR正演模拟

第四篇内容：

1、GPR剖面获取及成像

2、简单数据处理及分析

   探地雷达（Ground-Penetrating-Radar，GPR）是一种使用高频电磁波探测近地表地下结构的地球物理勘探方法，该方法对地表无破坏，操作简便，已被广泛应用于考古、城市道路建设、岩土工程调查等领域，有着非常良好的使用效果。但是，当前公开的GPR实际资料较少且难以形成庞大的知识库作为先验性息以备研究人员对实测剖面进行识别，所以，使用GPR正演模拟获得一套资料就显得尤为重要了。

        当前，基于Maxwell方程的GPR正演模拟主要使用的方法有：

        1、有限元方法（论文和开源软件都有）

        2、时域有限差分方法(GPRMax软件已经集成，但该软件没有界面，使用需要通过命令行程序执行)

        其中，时域有限差分方法（Finite-Difference-Time-Domain，FDTD）由K.S.Yee于1966年提出，是一种非常成熟且概念简单的正演模拟方法，因此，本系列博客决定基于FDTD实现GPR正演模拟。

一、时域有限差分方法

        首先从Maxwell旋度方程及其本构方程讲起：

      Maxwell旋度方程：

本构关系：

在上述公式中，E是电场强度（V/m），H是磁场强度（H/m），D是电通量密度（C/m^2）,B是磁通量密度（Wb/m^2），J是电流密度（A/m^2），介电常数,磁导常数 。

        将旋度方程转换为差分式是比较容易理解的，这里为了方便表示，使用了矩阵形式表示二维微分方程形式的TE波波动方程：

其中，

参照参考文献1中的Yee元胞结构：

        通过Yee元胞离散后，就可以对上述波动方程公式进行空间取样。

        这时要注意：虽然文献1中的差分取样是中心差分，该方式的好处是可以得到二阶计算精度（具体可以通过Taylor公式计算），但是计算机中无法通过上述方式实现，只能是在规定点取样，所以，理论上是二阶精度，但是，计算中仍然是一阶精度。

        首先是系数矩阵问题，公式详见参考文献1P18-P19：

由n状态更新到n+1状态就可以：

上述内容为使用FDTD实现GPR正演模拟所需要的基础知识，重新梳理一遍：通过Yee元胞对Maxwell方程进行离散，得到电磁场的更新方程。

参考文献：葛德彪,闫玉波. 电磁波时域有限差分法第2版．西安:西安电子科技大学出版社,2005：133-137,37,118．

Shen, H, Y., Li, X, S., Duan, R, F., et al., (2023), Quality evaluation of ground improvement by deep cement mixing piles via ground-penetrating radar. Nature Communications,14,34-48. http://doi.org/10.1038/s41467-023-39236-4

冯德山,戴前伟,何继善等. 探地雷达GPR正演的时域有限差分实现(英文)[J].地球物理学进展,2006,(02):630-636.

李静,刘津杰,曾昭发等. 基于变换光学有限差分探地雷达数值模拟研究[J].地球物理学报,2016,59(06):2280-2289.