Gprmax浅水域三维地质雷达数值模拟研究

前言

浅水域地下不良地质体的探测一直是工程勘察的难点，地质雷达具有仪器轻便、操作简洁、分辨率高的优势，在浅水域勘察中具有很大的应用前景。目前，二维地质雷达已经有不少应用，三维地质雷达理论上探测效果更好，本文利用开源的Gprmax软件进行三维地质雷达的数值模拟研究，优化了Python语言批量运行Gprmax的程序，可以根据工程需要自由的修改代码，同时基于Python语言开发了编写Gprmax输入in文件的小软件，用Matlab语言开发了三维地质雷达可视化的程序，也可根据工程需要自由的修改Matlab代码，代码仅供参考。

文章目录

Gprmax浅水域三维地质雷达数值模拟研究
- - 前言
- 1、优化的Python程序运行Gprmax
- 2、Pyqt开发的Gprmax输入in文件编辑小软件
- 3、Paraview与matlab三维可视化
- - 未完待续......三维模拟速度实在是太慢了，跑了一天还没出结果，后面再加上。

1、优化的Python程序运行Gprmax

利用Gprmax自带的工具函数绘图，调节图像的大小、字体、颜色、显示方式不是很方便，而且Gprmax中提供的B—scan没有波形变面积图显示直观。在tools.B-scan工具的基础上修改了代码，可以同时显示灰度图与波形图，首先看一下效果：
图1 管线模型
图2 空洞模型
废话不多说，直接上代码，代码不懂的地方，多运行几次就很清楚了。

"""
python运行gprmax
读取.in文件
运行api函数模拟
shangxiang
2023.6.2
"""import os
import tkinter as TK
import numpy as np
import matplotlib.pyplot as plt
from gprMax.gprMax import api
from tools.outputfiles_merge import get_output_data, merge_files
from tkinter import filedialog# 弹出窗口读取in文件
root = TK.Tk()
root.withdraw()
# 读取in文件，获取文件路径
filename = filedialog.askopenfilename()
# 把Gprmax需要用到的各种文件名设置好
merge_filename = filename[:-3]
merged_filename = filename[:-3]+'_merged.out'
ez_filename = filename[:-3]+'.txt'# 开始调用gprmax正演
# n：仿真次数（A扫描次数）->B扫描，gpu，使用gpu加速
# api(filename, n=61, geometry_only=False, gpu=[0])
# 将多个out文件合成merge_files
# merge_files(merge_filename, removefiles=True)# 获取回波数据
rxnumber = 1
rxcomponent = 'Ez'
outputdata, dt = get_output_data(merged_filename, rxnumber, rxcomponent)
# 保存回波数据
np.savetxt(ez_filename, outputdata, delimiter=' ')# 画图，不用gprmax自带的plot_Bscan工具，自编一个函数用来绘图
# 画两张图，灰度图和波形变面积图
def mpl_plot(filename, outputdata, dt, rxnumber, rxcomponent,scale):(path, filename) = os.path.split(filename)fig,axes = plt.subplots(1,2,squeeze=[1200,800],figsize=(10,6),sharex='all',sharey='all',constrained_layout=True)fig.suptitle( 'Trace number',  fontsize=20)im = axes[0].imshow(outputdata, extent=[0, outputdata.shape[1], outputdata.shape[0] * dt *1e9, 0], interpolation='nearest', aspect='auto', cmap='seismic', vmin=-np.amax(np.abs(outputdata)), vmax=np.amax(np.abs(outputdata)))axes[0].tick_params(labelsize=20)axes[0].xaxis.tick_top()axes[0].set_ylabel( 'Time[ns]',fontsize=20)axes[0].grid(which='both', axis='both', linestyle='-.')cb = fig.colorbar(im,ax=axes[0])cb.set_label('Field strength [V/m]')tw = outputdata.shape[1]             # 时间窗（与in文件一致）trace_number = len(outputdata[0])outputdata = scale*outputdata/np.max(outputdata)for i in range(trace_number):# x = outputdata[:,i]+i*space_signalx = outputdata[:,i] + iy = np.linspace(0,tw,len(outputdata))axes[1].plot(x, y, 'k',linewidth = 0.2)# c = i*space_signalc = iaxes[1].fill_betweenx(y, c, x,where=x>c, facecolor = 'black')axes[1].set_xlim(-2, trace_number+2)axes[1].set_ylim(tw,0)axes[1].tick_params(labelsize=20)    axes[1].xaxis.tick_top() # 保存图片为PDF或者PNG# savefile = os.path.splitext(filename)[0]# fig.savefig(path + os.sep + savefile + '.pdf', dpi=None, format='pdf', #             bbox_inches='tight', pad_inches=0.1)# fig.savefig(path + os.sep + savefile + '.png', dpi=150, format='png', #             bbox_inches='tight', pad_inches=0.1)return plt
scale = float(3)            # 波形压缩量# space_signal = float(1)   # 波形间隔(按实际情况变更)
B_plt = mpl_plot(merged_filename, outputdata, dt, rxnumber, rxcomponent,scale)
B_plt.show()

2、Pyqt开发的Gprmax输入in文件编辑小软件

Gprmax中in文件是记事本直接打开的文件，在记事本中进行编辑。为了方便，基于Pyqt开发了编辑in文件的小软件，与记事本差不多，但是使用更加方便。首先看一下软件界面。
在这里插入图片描述
软件添加了很多菜单栏，功能比记事本更加强大。

代码如下：

"""
一个简易的文本编辑器用pyqt编辑
此程序功能是用于Gprmax的in文件编辑
作者 shangxiang
时间 2023.6.7
"""# 导入所需的库
from PyQt5.QtGui import *
from PyQt5.QtWidgets import *
from PyQt5.QtCore import *
from PyQt5.QtPrintSupport import *import os
import sys
import uuid# 初始化参数
# 定义界面所需参数和函数
FONT_SIZES = [7, 8, 9, 10, 11, 12, 13, 14, 18, 24, 36, 48, 64, 72, 96, 144, 288]
IMAGE_EXTENSIONS = ['.jpg','.png','.bmp']
HTML_EXTENSIONS = ['.htm', '.html']def hexuuid():# 生成字符串类型的标识符，IPreturn uuid.uuid4()def splitext(p):# 分离文件名与扩展名return os.path.splitext(p)[1].lower()class TextEdit(QTextEdit):# 定义一个文本编辑类def canInsertFromMimeData(self, source):# 文本编辑器中插入图片if source.hasImage():return Trueelse:return super(TextEdit, self).canInsertFromMimeData(source)# def insertFromMimeData(self, source):#     # 获取文本区间的文本#     cursor = self.textCursor()#     document = self.document()#     if source.hasUrls():#         for u in source.urls():#             file_ext = splitext(str(u.toLocalFile()))#             if u.isLocalFile() and file_ext in IMAGE_EXTENSIONS:#                 image = QImage(u.toLocalFile())#                 document.addResource(QTextDocument.ImageResource, u, image)#                 cursor.insertImage(u.toLocalFile())#             else:#                 break#         else:#             return#     elif source.hasImage():#         image = source.imageData()#         uuid = hexuuid()#         document.addResource(QTextDocument.ImageResource, uuid, image)#         cursor.insertImage(uuid)#         return#     super(TextEdit, self).insertFromMimeData(source)class MainWindow(QMainWindow):def __init__(self, *args, **kwargs):super(MainWindow, self).__init__(*args, **kwargs)layout = QVBoxLayout()# 垂直布局self.editor = TextEdit()self.editor.setAutoFormatting(QTextEdit.AutoAll) # 启用自动创建格式化功能self.editor.selectionChanged.connect(self.update_format) # 当选择文本发生变化时，调用更新函数font = QFont('Times', 12)# 设置文本框的字体及字体大小self.editor.setFont(font)self.editor.setFontPointSize(12)self.path = Nonelayout.addWidget(self.editor)container = QWidget()container.setLayout(layout)self.setCentralWidget(container)self.status = QStatusBar()self.setStatusBar(self.status)file_toolbar = QToolBar("文件")file_toolbar.setIconSize(QSize(14, 14))self.addToolBar(file_toolbar)file_menu = self.menuBar().addMenu("&文件")open_file_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'cil-folder-open.svg')), "打开文件", self)open_file_action.setStatusTip("从本地磁盘中读取文件..")open_file_action.triggered.connect(self.file_open)file_menu.addAction(open_file_action)file_toolbar.addAction(open_file_action)save_file_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'cil-save.svg')), "保存", self)save_file_action.setStatusTip("保存到本地磁盘..")save_file_action.triggered.connect(self.file_save)file_menu.addAction(save_file_action)file_toolbar.addAction(save_file_action)saveas_file_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'folders.svg')), "另存为", self)saveas_file_action.setStatusTip("另存为文件..")saveas_file_action.triggered.connect(self.file_saveas)file_menu.addAction(saveas_file_action)file_toolbar.addAction(saveas_file_action)print_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'printer.svg')), "打印", self)print_action.setStatusTip("打印本页..")print_action.triggered.connect(self.file_print)file_menu.addAction(print_action)file_toolbar.addAction(print_action)edit_toolbar = QToolBar("编辑")edit_toolbar.setIconSize(QSize(16, 16))self.addToolBar(edit_toolbar)edit_menu = self.menuBar().addMenu("&编辑")undo_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'arrow-forward-up.svg')), "撤回", self)undo_action.setStatusTip("撤回上一个操作..")undo_action.triggered.connect(self.editor.undo)edit_menu.addAction(undo_action)redo_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'arrow-forward-up.svg')), "重做", self)redo_action.setStatusTip("重做撤回的操作..")redo_action.triggered.connect(self.editor.redo)edit_toolbar.addAction(redo_action)edit_menu.addAction(redo_action)edit_menu.addSeparator()cut_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'scissors.svg')), "剪切", self)cut_action.setStatusTip("剪切选定内容..")cut_action.setShortcut(QKeySequence.Cut)cut_action.triggered.connect(self.editor.cut)edit_toolbar.addAction(cut_action)edit_menu.addAction(cut_action)copy_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'cil-copy.svg')), "复制", self)copy_action.setStatusTip("复制选定内容..")cut_action.setShortcut(QKeySequence.Copy)copy_action.triggered.connect(self.editor.copy)edit_toolbar.addAction(copy_action)edit_menu.addAction(copy_action)paste_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'sticky.svg')), "粘帖", self)paste_action.setStatusTip("从剪贴板粘帖..")cut_action.setShortcut(QKeySequence.Paste)paste_action.triggered.connect(self.editor.paste)edit_toolbar.addAction(paste_action)edit_menu.addAction(paste_action)select_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'checkbox.svg')), "全选", self)select_action.setStatusTip("全选所有文字..")cut_action.setShortcut(QKeySequence.SelectAll)select_action.triggered.connect(self.editor.selectAll)edit_toolbar.addAction(select_action)edit_menu.addAction(select_action)edit_menu.addSeparator()wrap_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'arrow-down-circle.svg')), "自动换行", self)wrap_action.setStatusTip("当文字长度超过边框大小时自动换行..")wrap_action.setCheckable(True)wrap_action.setChecked(True)wrap_action.triggered.connect(self.edit_toggle_wrap)edit_toolbar.addAction(wrap_action)edit_menu.addAction(wrap_action)format_toolbar = QToolBar("格式")format_toolbar.setIconSize(QSize(16, 16))self.addToolBar(format_toolbar)format_menu = self.menuBar().addMenu("&格式")self.fonts = QFontComboBox()self.fonts.currentFontChanged.connect(self.editor.setCurrentFont)format_toolbar.addWidget(self.fonts)self.fontsize = QComboBox()self.fontsize.addItems([str(s) for s in FONT_SIZES])self.fontsize.currentIndexChanged[str].connect(lambda s: self.editor.setFontPointSize(float(s)) )format_toolbar.addWidget(self.fontsize)self.bold_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'cil-bold.svg')), "加粗", self)self.bold_action.setStatusTip("加粗选定内容..")self.bold_action.setShortcut(QKeySequence.Bold)self.bold_action.setCheckable(True)self.bold_action.toggled.connect(lambda x: self.editor.setFontWeight(QFont.Bold if x else QFont.Normal))format_toolbar.addAction(self.bold_action)format_menu.addAction(self.bold_action)self.italic_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'cil-italic.svg')), "斜体", self)self.italic_action.setStatusTip("将选定内容设为斜体..")self.italic_action.setShortcut(QKeySequence.Italic)self.italic_action.setCheckable(True)self.italic_action.toggled.connect(self.editor.setFontItalic)format_toolbar.addAction(self.italic_action)format_menu.addAction(self.italic_action)self.underline_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'cil-underline.svg')), "下划线", self)self.underline_action.setStatusTip("将选定内容加下划线..")self.underline_action.setShortcut(QKeySequence.Underline)self.underline_action.setCheckable(True)self.underline_action.toggled.connect(self.editor.setFontUnderline)format_toolbar.addAction(self.underline_action)format_menu.addAction(self.underline_action)format_menu.addSeparator()self.alignl_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'justify-left.svg')), "靠左对齐", self)self.alignl_action.setStatusTip("将文本靠左对齐..")self.alignl_action.setCheckable(True)self.alignl_action.triggered.connect(lambda: self.editor.setAlignment(Qt.AlignLeft))format_toolbar.addAction(self.alignl_action)format_menu.addAction(self.alignl_action)self.alignc_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'text-center.svg')), "居中对齐", self)self.alignc_action.setStatusTip("将文本居中对齐..")self.alignc_action.setCheckable(True)self.alignc_action.triggered.connect(lambda: self.editor.setAlignment(Qt.AlignCenter))format_toolbar.addAction(self.alignc_action)format_menu.addAction(self.alignc_action)self.alignr_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'justify-right.svg')), "靠右对齐", self)self.alignr_action.setStatusTip("将文本靠右对齐..")self.alignr_action.setCheckable(True)self.alignr_action.triggered.connect(lambda: self.editor.setAlignment(Qt.AlignRight))format_toolbar.addAction(self.alignr_action)format_menu.addAction(self.alignr_action)self.alignj_action = QAction(QIcon(os.path.join('shang_gprmax_qicon', 'justify.svg')), "对齐", self)self.alignj_action.setStatusTip("分散对齐文本..")self.alignj_action.setCheckable(True)self.alignj_action.triggered.connect(lambda: self.editor.setAlignment(Qt.AlignJustify))format_toolbar.addAction(self.alignj_action)format_menu.addAction(self.alignj_action)format_group = QActionGroup(self)format_group.setExclusive(True)format_group.addAction(self.alignl_action)format_group.addAction(self.alignc_action)format_group.addAction(self.alignr_action)format_group.addAction(self.alignj_action)format_menu.addSeparator()self._format_actions = [self.fonts,self.fontsize,self.bold_action,self.italic_action,self.underline_action]self.update_format()self.update_title()self.show()def block_signals(self, objects, b):for o in objects:o.blockSignals(b)def update_format(self):self.block_signals(self._format_actions, True)self.fonts.setCurrentFont(self.editor.currentFont())self.fontsize.setCurrentText(str(int(self.editor.fontPointSize())))self.italic_action.setChecked(self.editor.fontItalic())self.underline_action.setChecked(self.editor.fontUnderline())self.bold_action.setChecked(self.editor.fontWeight() == QFont.Bold)self.alignl_action.setChecked(self.editor.alignment() == Qt.AlignLeft)self.alignc_action.setChecked(self.editor.alignment() == Qt.AlignCenter)self.alignr_action.setChecked(self.editor.alignment() == Qt.AlignRight)self.alignj_action.setChecked(self.editor.alignment() == Qt.AlignJustify)self.block_signals(self._format_actions, False)def dialog_critical(self, s):dlg = QMessageBox(self)dlg.setText(s)dlg.setIcon(QMessageBox.Critical)dlg.show()def file_open(self):path, _ = QFileDialog.getOpenFileName(self, "Open file", "", "HTML documents (*.html);Text documents (*.txt);All files (*.*)")try:with open(path, 'rU') as f:text = f.read()except Exception as e:self.dialog_critical(str(e))else:self.path = pathself.editor.setText(text)self.update_title()def file_save(self):if self.path is None:return self.file_saveas()text = self.editor.toHtml() if splitext(self.path) in HTML_EXTENSIONS else self.editor.toPlainText()try:with open(self.path, 'w') as f:f.write(text)except Exception as e:self.dialog_critical(str(e))def file_saveas(self):path, _ = QFileDialog.getSaveFileName(self, "Save file", "", "HTML documents (*.html);Text documents (*.txt);All files (*.*)")if not path:returntext = self.editor.toHtml() if splitext(path) in HTML_EXTENSIONS else self.editor.toPlainText()try:with open(path, 'w') as f:f.write(text)except Exception as e:self.dialog_critical(str(e))else:self.path = pathself.update_title()def file_print(self):dlg = QPrintDialog()if dlg.exec_():self.editor.print_(dlg.printer())def update_title(self):self.setWindowTitle("%s - shang_Gprmax的in文件编辑" % (os.path.basename(self.path) if self.path else "Untitled"))def edit_toggle_wrap(self):self.editor.setLineWrapMode( 1 if self.editor.lineWrapMode() == 0 else 0 )if __name__ == '__main__':app = QApplication(sys.argv)app.setApplicationName("shang_Gprmax的in文件编辑")window = MainWindow()window.resize(1300,750)app.exec_()

3、Paraview与matlab三维可视化

在paraview中可以直接绘制三维可视化的模型，操作比较简单，下图是几个模型。
在这里插入图片描述

在paraview中也可以绘制三维的波场快照，如下：

在matlab中编程，将三维测线上雷达数据合成一个三维图，效果如下：

直接上matlab代码：

close all
clear
clc% tic
% 此程序是读取matlab文件夹下面的内容
% 绘制探地雷达C-scan图
% 作者：shangxiang
% 时间：2021年10月12日
% 地点：物探楼314% 读取文件夹中的文件名
path = 'E:\你的文件路径\';
file = dir(fullfile(path));
file_names = {file.name};
file_names(1:2) = [];% 创建三维数组
nz = size(dlmread([path file_names{1} '\'  'EData.txt']),1);
nx = size(dlmread([path file_names{1} '\'  'EData.txt']),2);
ny = length(file_names);data_3D = ones(ny, nx, nz);% 构造时间增益函数
xa = 1:nz;
ya = 1.3.^(xa*1e-2) - 0.5;
% 设置最大增益倍数
ymax = 25;
ya(ya>ymax) = ymax;for i = 1:length(file_names)file_name = file_names(i);file_name_root = [path file_name{1} '\' 'EData.txt'];
%     disp(file_name_root);data = (dlmread(file_name_root).*ya')';data_average = mean(data,1);% 去直达波
%    data(:,1:700) = data(:,1:700) - data_average(:,1:700);
%     disp(size(data));data_reshape = reshape(data,[1,nx,nz]);data_3D(i,:,:) = data_reshape;
end
tic
data_3D = gpuArray(single(data_3D));
%  dt = 1.9258332015464704e-11% data_3D(:,:,:) = 1
% data = dlmread(file_name_root);
imagesc(data');
colormap(gray);figure(2);
xslice = [];
yslice = 1:1:5;
zslice = [];
[x,y,z] = meshgrid(1:nx,1:ny,1:nz);
S = slice(x,y,z,data_3D,xslice,yslice,zslice);
set(gca,'zdir','reverse');
nd = linspace(0,nz,7);
set(gca,'ZMinorTick','off');
% set(gca,'ztick',[0 10 20 30 40 50 60]);
set(gca,'ztick',nd); 
set(gca,'zticklabel',{'0','10','20','30','40','50','60'}); xlabel('x_direction(m)','Interpreter','none');
ylabel('y_direction(m)','Interpreter','none');
zlabel('Time (ns)');
colormap(jet);
% colorbar
% lim = caxis;
caxis([-0.5 0.5]);
alpha(0.6);
% alpha('color');
alphamap('vdown');
alphamap('decrease',.2);
shading flat
toc