代码分享：gprMax钻孔地质雷达波场模拟

前言

gprMax模拟地面地质雷达被广泛使用，但是在钻孔内进行地质雷达的模拟较少。本博文尝试利用gprMax进行钻孔地质雷达的模拟，代码仅供大家借鉴。

1、钻孔地质雷达模拟效果

为了显示电磁波在钻孔内传播的状态，利用paraview将模拟的结果动态展示。

1.1、发射源在孔口

1.2、发射源在孔底

1.3、发射源在孔中

2、模型参数设置

假设钻孔间距为20m，钻孔深度为30m，背景介质设置为花岗岩，介电常数为6，电导率为0.005；设置两个方形充水空洞，介电常数为81，电导率为0.003。模型图如下：

雷达天线频率设置为32MHz，gprMax模拟的具体参数设置如下：

#title: krast_carbonatite_Bscan_2D
#domain: 21.000 30.000 0.050
#dx_dy_dz: 0.050 0.050 0.050
#time_window: 50e-8

#material: 6 0.005 1 0 granite
#material: 81 0.003 1 0 krast

#waveform: ricker 5 3.2e7 my_ricker

#hertzian_dipole: z 0.55 29.5 0 my_ricker

#python:
for i in range(1,31):
    print("#rx: 20.450 {} 0".format(i-0.5))
#end_python:

#src_steps: 0 0 0
#rx_steps: 0 0.000 0

#box: 0 0 0 21.000 30.000 0.050 granite
#box: 3.000 5.000 0 8.000 10.000 0.050 krast
#box: 11.000 19.000 0 16.000 24.000 0.050 krast

cylinder: 15.500 15.000 0 15.500 15.000 2.000 0.500 krast

geometry_view: 0 0 0 21.000 30.000 0.010 0.010 0.010 0.010 krast_carbonatite n

#python:
for i in range(1, 51):
    print("#snapshot: 0 0 0 21.000 30.000 0.050 0.050 0.050 0.050 {} snapshot{}".format(i*1e-8, i))
#end_python:

本参数设置的不是很合理，但是不影响模拟出结果。

3、python代码

运行模拟的主程序代码由python编辑，代码如下：

"""
python运行gprmax
读取.in文件
运行api函数模拟
"""

import os
import numpy as np
import matplotlib.pyplot as plt
from gprMax.gprMax import api
from tools.outputfiles_merge import get_output_data, merge_files

# 文件路径+文件名
dmax = r"D:\Learnfile\gprmaxSTU\guanxian3D"  # 项目目录
# filename = os.path.join(dmax, 'guanxian_granite_Bscan_2D.txt')

# # 正演  n：仿真次数（A扫描次数）->B扫描
# api(filename, n=98, geometry_only=False )  # geometry_only：仅几何图形
merge_files(dmax+'/'+'PVC03', removefiles=True)
7
# # # 获取回波数据
# # # A B扫描时out文件名不一样
filename = os.path.join(dmax+'/'+'PVC03'+'_merged.out')
rxnumber = 1
rxcomponent = 'Ez'
outputdata, dt = get_output_data(filename, rxnumber, rxcomponent)

# # # 保存回波数据
np.savetxt('shangxiang.txt', outputdata, delimiter=' ')

# # 提取振幅
# def get_ez_value(data):
#     data = abs(data)
#     get_Evalue = data.sum(axis = 0)
#     return get_Evalue

# # print(get_ez_value(outputdata))
# plt.plot(np.arange(1, outputdata.shape[1]+1, 1), get_ez_value(outputdata))
# plt.show()

# # 提取走时
# def get_firsttime_value(data):
#     [a, b] = np.shape(data)
#     door=np.zeros(b)  
#     trvp = []
#     for i in range(b):
#         door[i] = np.max(abs(data))*0.02
#         for j in range(a):
#             if abs(data[j][i]) > door[i]:
#                 trvp = np.append(trvp, data[j][i])
#                 break
#     return trvp

# a = get_firsttime_value(outputdata)
# print(a.shape)
# plt.plot(np.arange(1, outputdata.shape[1]+1, 1), get_firsttime_value(outputdata))
# plt.show()

## B扫描绘图
# from tools.plot_Bscan import mpl_plot
# plt = mpl_plot(filename, outputdata, dt*1e9, rxnumber, rxcomponent)
# plt.ylabel('Time [ns]')
# plt.show()

# # # # A扫描绘图
# # # from tools.plot_Ascan import mpl_plot
# # # from gprMax.receivers import Rx
# # # outputs = Rx.defaultoutputs
# # # outputs = ['Ez']
# # # print(outputs)
# # # plt = mpl_plot(filename, outputs)
# # # plt.show()
# # # # #
# # # #
# # # # A扫描图
# # # outputdata[1:200,]=0    ## 通过置零消除天线耦合波
# # # output = outputdata[:,0]  # 第i道A扫信号：序号从0开始
# # # plt.plot(output)
# # # plt.show()
# # # #
# # # #
# 堆叠波形
# space_signal = float(0.01)   # 信号间隔(按实际情况变更)
# tw = 10              # 时间窗（与in文件一致）
# trace_number = len(outputdata[0])
# for i in range(trace_number):
#     plt.plot(outputdata[:,i]+(i+1)*space_signal,np.linspace(0,tw,len(outputdata)),color='m')
# # plt.xticks(range(space_signal,trace_number*space_signal+1,space_signal),range(1,trace_number+1))
# plt.xlim(0, space_signal*(trace_number+2))
# plt.ylim(0, tw)
# # plt.xlabel('trace_number')
# plt.ylabel('Time [ns]')
# ax = plt.gca()          # 获取句柄
# ax.invert_yaxis()       # y轴反向
# ax.xaxis.tick_top()     # x轴放在上方
# plt.show() 

注：程序运行出现bug，注意文件目录与路径的调整。