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

         回顾上一章的内容，首先是探地雷达的使用范围和其主要面向的地球物理勘探对象，其次是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不同是为了便于差分式计算
    //Hz
    vector<vector<double>>Hz(cellNumForX, vector<double>(cellNumForY, 0.0));
	//Ex
    vector<vector<double>>Ex(cellNumForX, vector<double>(cellNumForY + 1, 0.0));
	//Ey
    vector<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边界条件需要存储前一时刻场量

，我做过测试，在0~0.5之间吸收效果较好，大于0.5后吸收效果较差，可以看看波长快照：

Mur边界不同时刻波长快照

从图中可以清晰看到，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);
	//上方Cpml
	for (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]);
	}
	//下方Cpml
	for (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]);
	}
	//左侧Cpml
	for (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]);
	}
	//右侧Cpml
	for (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]：

CPML边界吸收效果

        由图可知，边界、角点均未出现明显的反射，相较于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