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工具的基础上修改了代码，可以同时显示灰度图与波形图，首先看一下效果：
废话不多说，直接上代码，代码不懂的地方，多运行几次就很清楚了。
"""
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_signal
x = outputdata[:,i] + i
y = np.linspace(0,tw,len(outputdata))
axes[1].plot(x, y, 'k',linewidth = 0.2)
# c = i*space_signal
c = i
axes[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():
# 生成字符串类型的标识符，IP
return uuid.uuid4()
def splitext(p):
# 分离文件名与扩展名
return os.path.splitext(p)[1].lower()
class TextEdit(QTextEdit):
# 定义一个文本编辑类
def canInsertFromMimeData(self, source):
# 文本编辑器中插入图片
if source.hasImage():
return True
else:
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 = None
layout.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 = path
self.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:
return
text = 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 = path
self.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
