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

回顾上一章的内容，首先是探地雷达的使用范围和其主要面向的地球物理勘探对象，其次是Maxwell方程及FDTD基础知识，本章的内容包括：1、基于C++的TE波波动方程实现

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

一、基本波动方程实现：

使用C++实现，当然，可以选择的语言很多，比如Fortran、C/C++、JAVA、Python和Matlab等，但是后两者在网格剖分较细（元胞数量较多）时，这里仅限于在串行计算方式下讨论，并行计算方式下，Python可以使用PyCuda，Matlab也有相应的并行转换方式。下一章再做讨论。

大道至简，任何一个程序都是由多个函数（模块）构成的，FDTD也不例外，但是，要想使用FDTD实现GPR正演，则还需要根据实际情况进行分析，这里列举主要模块如下：

1、模型创建（网格剖分）函数；

2、FDTD计算函数；

3、激励源函数；

4、输出函数，其输出内容包括：波场快照、接收记录、模型等

5、边界处理函数（本部分在本章第二节讨论，主要涉及简单吸收边界和cpml吸收边界）

将上述函数组合，就可以实现一个基本的FDTD程序。

注意：本章C++程序均由Visual Studio编写,之前的博客中有讨论过简单数组、指针数组和vector的区别：二维数组创建方式比较-CSDN博客https://blog.csdn.net/faltas/article/details/132620446

但这里为了简便和可理解性，使用Vector容器进行编写。

首先是要导入的库和要定义好的参数：

//editor：WangBoHua
//2024.6.10
#include<iostream>
#include<fstream>
#include<vector>
#include<cmath>
using namespace std;
//光速
const double c = 3e8;
//圆周率
const double pi = 3.1415926;
//介电常数
const double eplison0 = 8.854e-12;
//磁导率
const double mu0 = 4 * pi * 1e-7;
//模型X、Y方向大小
const double UnderGroundX = 10;
const double UnderGroundY = 10;
//网格数量 
const int cellNumForX = 1000;
const int cellNumForY = 1000;
//边界厚度
const int cpmlBound=10
//Mur边界参数
const double Sc=0.3

1、模型创建函数

//editor：WangBoHua
//2024.6.10
//模型相对介电常数
vector<vector<double>> eplison(cellNumForX + 1, vector<double>(cellNumForY + 1, 6.0));
//模型电导率
vector<vector<double>> sigma(cellNumForX + 1, vector<double>(cellNumForY + 1, 1e-5));
void Model_Build(vector<vector<double>> eplison,vector<vector<double>> sigma){printMedia(eplison);
}

这里设计函数Model_Build，传参是eplison和sigma，当然，最开始也创建了一个均质模型用作参考。

2、FDTD计算函数：

该部分是本程序的核心内容，根据第一章的差分格式，可以按如下方式编写：

//editor：WangBoHua
//2024.6.10
void FDTD(int source,int receive,double frequency){//注意，这里Ex、Ey的大小和Hz不同是为了便于差分式计算//Hzvector<vector<double>>Hz(cellNumForX, vector<double>(cellNumForY, 0.0));//Exvector<vector<double>>Ex(cellNumForX, vector<double>(cellNumForY + 1, 0.0));//Eyvector<vector<double>>Ey(cellNumForX + 1, vector<double>(cellNumForY, 0.0));//系数参数vector<vector<double>>CA(cellNumForX + 1, vector<double>(cellNumForY + 1, 0.0));vector<vector<double>>CB(cellNumForX + 1, vector<double>(cellNumForY + 1, 0.0));vector<vector<double>>CP(cellNumForX + 1, vector<double>(cellNumForY + 1, 0.0));vector<vector<double>>CQ(cellNumForX + 1, vector<double>(cellNumForY + 1, 0.0));    //系数计算for (int i = 0; i < cellNumForX + 1; i++) {for (int j = 0; j < cellNumForY + 1; j++) {CA[i][j] = double((2 * eplison[i][j] * eplison0 - sigma[i][j] * Courant) / (2 * eplison[i][j] * eplison0 + sigma[i][j] * Courant));CB[i][j] = double((2 * Courant) / (2 * eplison[i][j] * eplison0 + sigma[i][j] * Courant));CP[i][j] = double((2 * mu0 - sigma[i][j] * Courant) / (2 * mu0 + sigma[i][j] * Courant));CQ[i][j] = double((2 * Courant) / (2 * mu0 + sigma[i][j] * Courant));}}//FDTD更新for(int n=0;n<simulation_steps,n++){//Hz更新for (int i = 0; i < cellNumForX; i++) {for (int j = 0; j < cellNumForY; j++) {Hz[i][j] = CP[i][j] * Hz[i][j] - CQ[i][j] * ((Ey[i + 1][j] - Ey[i][j]) / dx - (Ex[i][j + 1] - Ex[i][j]) / dy);}}//Ex更新for (int i = 0; i < cellNumForX; i++) {for (int j = 1; j < cellNumForY; j++) {Ex[i][j] = CA[i][j] * Ex[i][j] + CB[i][j] * (Hz[i][j] - Hz[i][j - 1]) / dy;         }}//Ey更新for (int i = 1; i < cellNumForX; i++) {for (int j = 0; j < cellNumForY; j++) {Ey[i][j] = CA[i][j] * Ey[i][j] - CB[i][j] * (Hz[i][j] - Hz[i - 1][j]) / dx;}}//输出波场快照打印if(n%100==0){printSlice(Hz,n);}    
}//清理内存，防止程序崩溃Ex.shrink_to_fit();Ey.shrink_to_fit();Hz.shrink_to_fit();}

3、激励源函数:

GPR一般使用类似地震勘探中的雷克子波，所以这里激励源也使用雷克子波：

//editor：WangBoHua
//2024.6.10
double Ricker(double time, double frequency) {double ricker = (1 - 2 * pow(pi * (time - 1 / frequency) * frequency, 2)) * exp(-pow(pi * (time - 1 / frequency) * frequency, 2));//这是地震里常用的雷克子波return ricker;
}

4、输出函数

其实之前的代码中已经有了，如第一节中的printMedia输出模型文件，第二节中的printSlice输出不同时刻波场快照，以及之后我们要在FDTD函数内集成的输出接受记录的语句；这里为了简便，直接写了：

//editor：WangBoHua
//2024.6.10
//1、模型打印
void printMedia(vector <vector<double>> media) {fstream ofs;string Filename = "C:/Users/L/Desktop/Simulation/out/Media.data";ofs.open(Filename, ios::out);ofs.clear();if (ofs.is_open()) {for (int i = 0; i < cellNumForX ; i++) {for (int j = 0; j < cellNumForY ; j++) {ofs << media[i][j] << "\t";}ofs << endl;}cout << "相应模型创建完毕" << endl;}else {cout << "wrong!!!check the path!!!" << endl;}ofs.close();
}
//2、波长快照输出
void printSlice(vector<vector<double>> photo, int n) {fstream ofs;char filename[200];sprintf_s(filename, "%s%d%s", "C:/Users/12981/Desktop//", int(n * Courant * 1e9), "ns.data");ofs.open(filename, ios::out);ofs.clear();if (ofs.is_open()) {for (int i = 0; i < cellNumForX; i++) {for (int j = 0; j < cellNumForY; j++) {ofs << photo[i][j] << '\t';}ofs << endl;}}else {cout << "Check the path you have,I can't open the f**king floder!!! ";}ofs.close()
}

接收记录是这样的：

//editor：WangBoHua
//2024.6.10	
string filename = //输出文件地址fstream ofs;ofs.open(filename, ios::out);ofs.clear();if(ofs.is_open()){FDTD更新if(n%Sample==0){for(int k=0;k<cellNumForX;k++){ofs<<Ey[0][k]<<'\t';}ofs<<endl;}    }ofs.close();

好了，一个没有处理边界条件的最简单FDTD就可以通过组装上述函数获得了，根据参考文献一，在特定位置添加线激励源：

Ey[source][receive] += double(Courant / eplison0) * Ricker(Courant * n, frequency) / (dx * dy);

为了验证我们计算结果的正确性（主要验证波长快照中的雷达波是否以球面波形式传播），本章提供了简单的Python绘图程序：

# editor：WangBoHua
# 2024.6.10
import numpy as np
import matplotlib.pyplot as plt
# 绘制波场快照
Snapshot=np.genfromtxt("C:/Users/12981/Desktop/10ns.data")
Snapshot2=np.genfromtxt("C:/Users/12981/Desktop/12ns.data")
Snapshot3=np.genfromtxt("C:/Users/12981/Desktop/14ns.data")
Snapshot4=np.genfromtxt("C:/Users/12981/Desktop/16ns.data")
Snapshot5=np.genfromtxt("C:/Users/12981/Desktop/18ns.data")
Snapshot6=np.genfromtxt("C:/Users/12981/Desktop/20ns.data")
fig=plt.figure(figsize=(16,9),dpi=1000,layout="constrained")
# x方向长度
# y方向长度
xlength = np.linspace(0, 10, 1000, dtype='float')
ylength = np.linspace(0, 10, 1000, dtype='float')
ax1 = fig.add_subplot(2, 3, 1)
ax1.pcolormesh(xlength, ylength, Snapshot, vmin=-100,vmax=100, cmap='jet', rasterized=True)
ax1.tick_params(direction='in', width=0.5)
ax1.invert_yaxis()
ax1.set_xlabel("x/m")
ax1.set_ylabel("z/m")
ax1.xaxis.set_ticks_position('top') 
ax2 = fig.add_subplot(2, 3,  2)
ax2.pcolormesh(xlength, ylength, Snapshot2, vmin=-100,vmax=100, cmap='jet', rasterized=True)
ax2.tick_params(direction='in', width=0.5)
ax2.invert_yaxis()
ax2.xaxis.set_ticks_position('top') ax3 = fig.add_subplot(2, 3,  3)
ax3.pcolormesh(xlength, ylength, Snapshot3, vmin=-100,vmax=100, cmap='jet', rasterized=True)
ax3.tick_params(direction='in', width=0.5)
ax3.invert_yaxis()
ax3.xaxis.set_ticks_position('top') ax4 = fig.add_subplot(2, 3,4)
ax4.pcolormesh(xlength, ylength, Snapshot4, vmin=-100,vmax=100, cmap='jet', rasterized=True)
ax4.tick_params(direction='in', width=0.5)
ax4.invert_yaxis()
ax4.xaxis.set_ticks_position('top') ax5 = fig.add_subplot(2, 3,  5)
ax5.pcolormesh(xlength, ylength, Snapshot5, vmin=-100,vmax=100, cmap='jet', rasterized=True)
ax5.tick_params(direction='in', width=0.5)
ax5.invert_yaxis()
ax5.xaxis.set_ticks_position('top') ax6 = fig.add_subplot(2, 3,  6)
ax6.pcolormesh(xlength, ylength, Snapshot6, vmin=-100,vmax=100, cmap='jet', rasterized=True)
ax6.tick_params(direction='in', width=0.5)
ax6.invert_yaxis()
ax6.xaxis.set_ticks_position('top') plt.show()

从图中可以看到，不同时间的波都是在以同心圆形式传播的，那就证明我们这里写的基础版本是没错的。

二、边界效应处理

通俗的讲，由于我们计算机的内存大小是固定的，但是工区是很大的，这就导致无法使用计算机完全模拟雷达波在地下的传播过程，同时在模型边界处会出现非物理性质的反射，对模拟波场和接收记录有着非常大的影响，继续第一节的运算，将时间放大些看看波场会怎样：

未处理边界的情况下，雷达波波前在截断边界处发生全反射并向波场内部传播。

边界条件处理方法一般有两种，一种是Mur边界，该边界条件是一种类似插值计算的方法，另一种是PML边界条件及其衍生出来的UPML、CPML等，但是根据文献1的推导，UPML和CPML虽然在形式上不一样，但是二者是等价的，这里先从Mur边界开始，但公式推导这里不再展开，直接提供代码，如要对边界条件有更多的理解，文献1和文献2中的边界条件可以作为参考，尤其是文献1中的边界条件推导，非常扎实详细。

这里主要改变的是文献2中的Mur边界，该书中使用的是TM波，但无妨，TE波和TM波通过转换参数等就可以实现，但这里不做解读，直接使用。

		//editor：WangBoHua//2024.6.10//Mur-Boundary//左侧边界for (int i = 1; i < cellNumForY; i++){Hz[0][i] = Hz0[1][i] + float((Sc - 1) / (Sc + 1) * (Hz[1][i] - Hz0[0][i]));}//上边界for (int i = 1; i < cellNumForX; i++) {Hz[i][cellNumForY - 1] = Hz0[i][cellNumForY - 2] + float((Sc - 1) / (Sc + 1) * (Hz[i][cellNumForY - 2] - Hz0[i][cellNumForY - 1]));}//下边界for (int i = 1; i < cellNumForX; i++) {Hz[i][0] = Hz0[i][1] + float((Sc - 1) / (Sc + 1) * (Hz[i][1] - Hz0[i][0]));}//右侧边界for (int i = 1; i < cellNumForY; i++) {Hz[cellNumForX - 1][i] = Hz0[cellNumForX - 2][i] + float((Sc - 1) / (Sc + 1) * (Hz[cellNumForX - 2][i] - Hz0[cellNumForX - 1][i]));}//四周角点Hz[0][0] = 0.5 * (Hz[0][1] + Hz[1][0]);Hz[cellNumForX - 1][0] = 0.5 * (Hz[cellNumForX - 2][0] + Hz[cellNumForX - 1][1]);Hz[0][cellNumForY - 1] = 0.5 * (Hz[0][cellNumForY - 2] + Hz[1][cellNumForY - 1]);Hz[cellNumForX - 1][cellNumForY - 1] = 0.5 * (Hz[cellNumForX - 2][cellNumForY - 1] + Hz[cellNumForX - 1][cellNumForY - 2]);
//注意，这里的Mur边界条件需要存储前一时刻场量

$Sc \in [0,1]$ ，我做过测试，在0~0.5之间吸收效果较好，大于0.5后吸收效果较差，可以看看波长快照：

从图中可以清晰看到，Mur边界的确可以吸收绝大部分的能量，但是，计算区域中仍存在较为明显的反射波，对绘图工具中的能量定义是[-100,100]，而反射波能量清晰可见，这意味着反射波的能量还是很强的，说明Mur边界的吸收效果较差。

2、PML及CPML边界，参考文献1和文献2。

文献1中说明，PML的理论较为混乱，并且其分裂了计算区域，可以参考《并行编程原理与程序设计》后的声波方程计算代码，分裂计算区域的情况是在计算区域外镶边（添加一层CPML边界），计算较为困难，为此提出了卷积完美匹配层CPML，计算方式简单且高效。

//editor：WangBoHua
//2024.6.10
void E_CPML(vector<double>Eax,vector<double>Eay,vector<double>Ebx，vector<double>Eby) {double kmax = 6;double alphamax = 0.04;double m = 4;double sigmaMax = double((1 + m) / (dx 150*pi*sqrt(eplison0)));vector<double> kEx(cellNumForX + 1, 0.0), alphaEx(cellNumForX + 1, 0.0), sigmaEx(cellNumForX + 1, 0.0);vector<double>  kEy(cellNumForY + 1, 0.0), alphaEy(cellNumForY + 1, 0.0), sigmaEy(cellNumForY + 1, 0.0);//上方Cpmlfor (i = 0; i < CpmlBound; i++){sigmaEx[i] = sigmaMax * pow(double(CpmlBound - i - 1) / double(CpmlBound - 1), m);kEx[i] = 1 + (kmax - 1.e0) * pow(double(CpmlBound - i - 1) / double(CpmlBound - 1), m);alphaEx[i] = alphamax * pow(double(i) / double(CpmlBound - 1), ma);Ebx[i] = exp(-(sigmaEx[i] / kEx[i] + alphaEx[i]) * dt / eplison0);Eax[i] = sigmaEx[i] * (Ebx[i] - 1.e0) / (sigmaEx[i] * kEx[i] + kEx[i] * kEx[i] * alphaEx[i]);}//下方Cpmlfor (i = 0; i < CpmlBound; i++){k = cellNumForX - i;sigmaEx[k] = sigmaMax * pow(double(CpmlBound - i - 1) / double(CpmlBound - 1), m);kEx[k] = 1 + (kmax - 1.e0) * pow(double(CpmlBound - i - 1) / double(CpmlBound - 1), m);alphaEx[k] = alphamax * pow(double(i) / double(CpmlBound - 1), ma);Ebx[k] = exp(-(sigmaEx[k] / kEx[k] + alphaEx[k]) * dt / eplison0);Eax[k] = sigmaEx[k] * (Ebx[k] - 1.e0) / (sigmaEx[k] * kEx[k] + kEx[k] * kEx[k] * alphaEx[k]);}//左侧Cpmlfor (j = 0; j < CpmlBound; j++){sigmaEy[j] = sigmaMax * pow(double(CpmlBound - j - 1) / double(CpmlBound - 1), m);kEy[j] = 1 + (kmax - 1.e0) * pow(double(CpmlBound - j - 1) / double(CpmlBound - 1), m);alphaEy[j] = alphamax * pow(double(j) / double(CpmlBound - 1), ma);Eby[j] = exp(-(sigmaEy[j] / kEy[j] + alphaEy[j]) * dt / eplison0);Eay[j] = sigmaEy[j] * (Eby[j] - 1.e0) / (sigmaEy[j] * kEy[j] + kEy[j] * kEy[j] * alphaEy[j]);}//右侧Cpmlfor (j = 0; j < CpmlBound; j++){k = cellNumForY - j;sigmaEy[k] = sigmaMax * pow(double(CpmlBound - j - 1) / double(CpmlBound - 1), m);kEy[k] = 1 + (kmax - 1.e0) * pow(double(CpmlBound - j - 1) / double(CpmlBound - 1), m);alphaEy[k] = alphamax * pow(double(j) / double(CpmlBound - 1), ma);Eby[k] = exp(-(sigmaEy[k] / kEy[k] + alphaEy[k]) * dt / eplison0);Eay[k] = sigmaEy[k] * (Eby[k] - 1.e0) / (sigmaEy[k] * kEy[k] + kEy[k] * kEy[k] * alphaEy[k]);}}

//editor：WangBoHua
//2024.6.10
void H_CPML(vector<double>Hax,vector<double>Hay,vector<double>Hbx，vector<double>Hby){double kmax = 6;double alphamax = 0.04;double m = 4;double sigmaMax = double((1 + m) / (dx 150*pi*sqrt(eplison0)));for (i = 0; i < CpmlBound - 1; i++){sigmaHx[i] = sigmaMax * pow(double(CpmlBound - i - 1) / double(CpmlBound - 1), m);kHx[i] = 1 + (kmax - 1) * pow(double(CpmlBound - i - 1) / double(CpmlBound - 1), m);alphaHx[i] = alphamax * double(i + 1) / double(CpmlBound - 1);Hbx[i] = exp(-(sigmaHx[i] / kHx[i] + alphaHx[i]) * dt / eplison0);Hax[i] = sigmaHx[i] * (Hbx[i] - 1) / (sigmaHx[i] * kHx[i] + kHx[i] * kHx[i] * alphaHx[i]);}for (i = 0; i < CpmlBound - 1; i++){k = cellNumForX - i - 1;sigmaHx[k] = sigmaMax * pow(double(CpmlBound - i - 1) / double(CpmlBound - 1), m);kHx[k] = 1 + (kmax - 1) * pow(double(CpmlBound - i - 1) / double(CpmlBound - 1), m);alphaHx[k] = alphamax * double(i + 1) / double(CpmlBound - 1);Hbx[k] = exp(-(sigmaHx[k] / kHx[k] + alphaHx[k]) * dt / eplison0);Hax[k] = sigmaHx[k] * (Hbx[k] - 1) / (sigmaHx[k] * kHx[k] + kHx[k] * kHx[k] * alphaHx[k]);}for (j = 0; j < CpmlBound - 1; j++){sigmaHy[j] = sigmaMax * pow(double(CpmlBound - j - 1) / double(CpmlBound - 1), m);kHy[j] = 1 + (kmax - 1) * pow(double(CpmlBound - j - 1) / double(CpmlBound - 1), m);alphaHy[j] = alphamax * double(j +1) / double(CpmlBound - 1);Hby[j] = exp(-(sigmaHy[j] / kHy[j] + alphaHy[j]) * dt / eplison0);Hay[j] = sigmaHy[j] * (Hby[j] - 1) / (sigmaHy[j] * kHy[j] + kHy[j] * kHy[j] * alphaHy[j]);}for (j = 0; j < CpmlBound - 1; j++){k = cellNumForY - j - 1;sigmaHy[k] = sigmaMax * pow(double(CpmlBound - j - 1) / double(CpmlBound - 1), m);kHy[k] = 1 + (kmax - 1) * pow(double(CpmlBound - j - 1) / double(CpmlBound - 1), m);alphaHy[k] = alphamax * double(j + 1) / double(CpmlBound - 1);Hby[k] = exp(-(sigmaHy[k] / kHy[k] + alphaHy[k]) * dt / eplison0);Hay[k] = sigmaHy[k] * (Hby[k] - 1.e0) / (sigmaHy[k] * kHy[k] + kHy[k] * kHy[k] * alphaHy[k]);}}

//editor：WangBoHua
//2024.6.10
//FDTD函数中
void FDTD(int source,int receive,double frequency){/*if (i < CpmlBound - 1 || i >= cellNumForX - CpmlBound + 1) {PsiHzx[i][j] = Hbx[i] * PsiHzx[i][j] + Hax[i] * (Ey[i + 1][j] - Ey[i][j]) / dx;Hz[i][j] = Hz[i][j] - CQ[i][j] * PsiHzx[i][j];}if (j < CpmlBound - 1 || j >= cellNumForY - CpmlBound + 1) {PsiHzy[i][j] = Hby[j] * PsiHzy[i][j] + Hay[j] * (Ex[i][j + 1] - Ex[i][j]) / dy;Hz[i][j] = Hz[i][j] + CQ[i][j] * PsiHzy[i][j];}*//*					if (j < CpmlBound - 1 || j > cellNumForY - CpmlBound + 2) {PsiExy[i][j] = Ebx[j] * PsiExy[i][j] + Eax[j] * (Hz[i][j] - Hz[i][j - 1]) / dy;Ex[i][j] = Ex[i][j] + CB[i][j] * PsiExy[i][j];}*//*					if (i < CpmlBound - 1 || i > cellNumForY - CpmlBound + 2) {PsiEyx[i][j] = Eby[i] * PsiEyx[i][j] + Eay[i] * (Hz[i][j] - Hz[i - 1][j]) / dx;Ey[i][j] = Ey[i][j] - CB[i][j] * PsiEyx[i][j];}*///将上述语句分别添加至相应场量的更新循环中即可}

来看一看CPML的吸收效果，这里的能量范围仍是[-100，100]：

由图可知，边界、角点均未出现明显的反射，相较于Mur边界，CPML边界的吸收效果更好。

声明：欢迎上述代码和图片引用，但请标注公式和图片来源,创作不易，请多多支持，非常感谢！

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

梅中磊，李月娥，马阿宁等，MATLAB电磁场与微波技术仿真.清华大学出版社.

冯德山,探地雷达数值模拟及程序实现.

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. https://doi.org/10.1038/s41467-023-39236-4

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

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