探地雷达频率域2.5维正演

戴世坤 欧阳振崇* 周印明 张钱江 李 昆赵东东 陈轻蕊 凌嘉宣

①(中南大学地球科学与信息物理学院 长沙 410083)

②(中南大学有色金属成矿预测与地质环境监测教育部重点实验室 长沙 410083)

③(东方地球物理公司综合物化探处 涿州 072751)

④(桂林理工大学地球科学学院 桂林 541006)

1 引言

探地雷达(Ground Penetrating Radar, GPR)是一种利用高频(106～109Hz)电磁波来确定介质内部物质分布规律的地球物理方法，因其具有抗干扰能力强、探测分辨率高、场地适应能力强、操作简单且为无损探测等优势而被广泛地应用于工程勘察及地质调查。GPR正演是雷达理论研究的重点之一，随着勘探要求更加精细化和高效，大范围高效、高精度的数值模拟和反演成像成为现在GPR技术的目标[1,2]。

波动方程正演模拟方法，因其能包含雷达波的运动学特征和动力学特征被广泛应用于GPR正演模拟中，主要包括时域有限差分法(Finite-Difference Time Domain method, FDTD)和有限单元法(Finite Element Method, FEM)两种。在FDTD法的应用方面，自1996年Yee[3]提出著名的Yee氏网格后，FDTD法被广泛应用于GPR数值模拟中。Bergmann等人[4]应用FDTD法开展了不均匀非线性和衰减介质的GPR正演模拟；Irving等人[5]采用PML吸收边界研究了TE和TM模式的2维GPR数值模拟；刘四新等人[6]对比了GPR2维频散介质与非频散介质中GPR信号的区别；冯德山等人[7]研究了FDTD数值模拟中不分裂卷积完全匹配层对倏逝波的吸收效果。

FDTD法原理简单，易于编程实现，但要求模型规则剖分，对复杂问题适应性差。FEM因具有网格剖分灵活，适用于物性参数分布复杂和几何特征不规则模型的优势，被引入到GPR数值模拟领域。底青云等人[8]推导GPR有限元方程，开展了一系列典型模型的正演模拟；杜华坤等人[9]基于优化的Delaunay三角剖分，采用线性插值进行了FEM的GPR2维正演，提高了FEM对复杂模型模拟的适应性和数值模拟精度；石明等人[10]采用Delaunay非结构化网格有限元法进行了各向异性介质GPR有限元正演；王洪华等人[11]实现了PML边界条件在2阶电磁波动方程GPR时域有限元模拟中的应用，验证了PML边界条件在复杂地电结构电磁波传播模拟中具有良好的吸收效果。

近年来GPR数值模拟大都在时间域内研究，时间域GPR波的传播满足的波动方程模拟结果符合雷达波传播的运动学特征，但在表现波的动力学特征方面存在不足，文献[12]指出频率域波形反演在地震中的重要地位，首次研究了频率域波形反演，对频域波形反演存在的问题进行了初次探讨，本文从频率域出发，研究了2.5维GPR数值模拟方法，这样可以很好地保留雷达波传播的运动学特征和动力学特征，准确地研究雷达波在频率域的传播特性，为GPR全波形反演提供重要基础。本文采用有限单元法，推导基于行波分解的吸收边界条件的GPR有限元方程，实现频率域2.5D探地雷达正演模拟。重点分析和总结在雷达频率不同相对介电常数和不同收发距波数选取的规律；设计均匀全空间和半空间模型，将数值解与解析解对比验证了算法的正确性；另外，雷达2.5维频率域正演在不同波数之间的计算具有高度并行性，通过设计Open MP并行，探究不同线程下算法并行的效率，验证了算法的高效性。

2 方法理论

2.1 控制方程

2.2 有限单元法

3 数值算例

本文算法的程序代码采用Fortran语言编写，测试平台配置为CPU-Inter(R) Core(TM) i9-7980XE，主频2.60 GHz(18核，36线程)，内存64 GB,64位操作系统。

3.1 波场参数分析

图1 dx=0.1 m处不同相对介电常数的电磁场谱随波数变化特征

图2 dx=0.5 m处不同相对介电常数电磁场谱随波数变化特征

图3 dx=1 m处不同相对介电常数电磁场谱随波数变化特征

图4 dx=5 m处不同相对介电常数电磁场谱随波数变化特征

图5 dx=10 m处不同相对介电常数电磁场谱随波数变化特征

结论：频率为100 MHz，计算距离范围为10 m时，波数最大值选取103即可将谱的能量全包含其中，而且谱的能量主要分布在波数0～10的范围内，可对该范围内波数适当加密，使傅里叶逆变换精度更高。

3.2 正确性验证

3.2.1 全空间模型

设计均匀全空间模型，σ =0.001 S/m , εr=1.0，节点301×301，电流1000 A，偶极距dl为1 m，源中心在原点，x方向源网格均0.02 m，源以外网格采用非均匀剖分最大0.05 m，边界取15个网格间距由0.05 m以1.5倍递增至1 m。x, z方向网格剖分相同，f=100 MHz，正波数范围(0.01, 1000)，波数55个。

图6中主剖面上电磁场数值解与解析解的形态、数值都相同，主剖面测线z=0.5 m处各节点的相对误差均小于1%，验证了本文算法的正确性。

图7为采用文献[15]中的算法计算主剖面电磁场的数值解与解析解对比图，由图可知，文献[15]算法在源附近的误差较大，表明本文算法在源的处理上要优于文献[15]中的算法，计算精度更高。

3.2.2 半空间模型

图6 主剖面电磁场数值解与解析解对比图

图7 主剖面电磁场其他算法数值解与解析解对比图

设计均匀半空间模型，空气层电导率σ =10−12S/m ，地下介质电导率σ =10−3S/m，相对介电常数 εr=1.0，频率f=100 MHz，网格节点601×601，x方向的源布设于地下0.5 m，源的其他参数和网格剖分与全空间相同，选取正波数范围(0.0001, 1000)，波数277个。

图8为主剖面(y = 0 m)数值解与解析解对比，第3列为剖面测线z=0.5 m的相对误差，由图可知，数值解与解析解的形态、数量级相同，测线各节点的相对误差均小于1.5%，同样验证了本文算法的准确性。

3.3 Open MP并行测试

本文算法耗时主要在波数循环计算上，每个波数均需求解1次方程组，但各波数计算相互独立，可采用并行方式提高计算效率。目前电磁法2.5D正反演应用较多的并行方法有MPI(Message Passing Interface)和OpenMP(Open Multi-Processing)。OpenMP使用线程间共享内存的方式协调并行计算，对原串行代码改动小，容易实现。本文采用OpenMP并行处理不同波数方程组的求解、傅里叶逆变换和辅助场计算。

2.5D GPR正演计算量大，将Pardiso求解器采用OpenMP并行求解，算法效率如表1。模型参数与3.2.1小节模型一致，网格节点301×301，波数277个。式中加速比 = 单机的单线程耗时/并行计算耗时。

表1中，随着并行线程个数增加，加速比增大，算法耗时减少，改造后的计算效率明显提高。2线程耗时比单线程耗时长，可能原因是2线程时发生了同步事件，耗时变长[16]；16线程耗时约80 s，比单线程计算快了3倍；增加到20线程时，耗时反倒增大，分析原因是当计算量一定时，线程数量增加，用于线程通信/线程调度的时间所占比例逐渐增大，计算效率降低。

结论：采用OpenMP将Pardiso求解器并行，随着并行线程个数增加，计算效率提高，但并行线程数并非越多越好，综合不同线程算法耗时和计算机资源的占用情况，本节线程数为16时，加速比最大，耗时最短，是当前条件下OpenMP并行选取的最佳线程数。

3.4 复杂模型算例

设计均匀半空间中存在两个异常体，异常体的模型参数如图9所示，源和网格等模型参数均与3.2.2小节中半空间模型相同。

图10为主剖面电磁场响应特征，图中可看出异常体的位置和形状，低阻异常体相对介电常数大，电磁波长小，电磁场波动明显；高阻体相对介电常数略小，异常反应不如低阻体灵敏，电磁场波动小，说明波场对介质不均匀性有一定程度的敏感度，说明本文提出的GPR正演算法对于复杂模型具有很好的适应性。

图8 主剖面半空间电磁场数值解与解析解的对比图

表1 整个程序OpenMP并行不同线程计算效率

图9 复杂模型示意图

4 结论

本文采用有限单元法，推导了基于行波分解的吸收边界条件的GPR有限元方程，保留雷达波传播的运动学特征和动力学特征，实现了频率域2.5D探地雷达正演模拟。主要有以下结论：

(2) 均匀全空间和半空间模型数值解与解析解除源处之外的剖面节点误差均小于1%，表明算法正确。

图10 主剖面电磁场分量实部和虚部

(3) 测试平台CPU-Inter(R) Core(TM) i9-7980XE，主频2.60 GHz(36核)，内存64 GB，64位操作系统，节点301×301、波数277时，并行最优线程数是16，此时算法耗时80 s，比一个线程计算时的速度快了3倍，验证了算法的高度并行性和高效性。

(4) 传统探地雷达的发射源在空气中，本文采用接地线源，这样做更方便能量的辐射，研究雷达波在频率域的传播特性，完善雷达的电磁理论，为GPR全波形反演提供重要基础。