回顾上一章的内容,首先是探地雷达的使用范围和其主要面向的地球物理勘探对象,其次是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()
![](https://img-blog.csdnimg.cn/direct/734e1dc000c94615ae2ef308ad2157a1.png)
从图中可以看到,不同时间的波都是在以同心圆形式传播的,那就证明我们这里写的基础版本是没错的。
二、边界效应处理
通俗的讲,由于我们计算机的内存大小是固定的,但是工区是很大的,这就导致无法使用计算机完全模拟雷达波在地下的传播过程,同时在模型边界处会出现非物理性质的反射,对模拟波场和接收记录有着非常大的影响,继续第一节的运算,将时间放大些看看波场会怎样:
未处理边界的情况下,雷达波波前在截断边界处发生全反射并向波场内部传播。
边界条件处理方法一般有两种,一种是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后吸收效果较差,可以看看波长快照:
![](https://img-blog.csdnimg.cn/direct/c15089435e3841fc9afc8195ff0cba37.png)
从图中可以清晰看到,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]:
![](https://img-blog.csdnimg.cn/direct/4c47fe03b41c4e498b5bffaa6cec4844.png)
由图可知,边界、角点均未出现明显的反射,相较于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