探地雷达数值模拟及全波形反演

作者简介

内容简介

第1章 绪论
　　1.1 探地雷达正反演研究的目的及意义
　　探地雷达（ground penetrating radar，GPR）是一种对地下异常体、地下结构特性或物体内部不可见目标体进行定位的电磁无损探测技术（李大心，1994）。GPR作为一种新兴的地球物理手段，可以在不破坏地表结构的情况下获取地下未知区域的目标体信息，有探测速度快、探测精度高、抗干扰能力强、结果直观、安全、成本低等特点，已发展成为工程、水文环境、地质灾害、考古、军事等应用领域的一种重要近地表地球物理方法（粟毅等，2006；曾昭发等，2006）。
　　近年来，随着GPR勘探领域的不断延伸、勘探环境更加复杂、勘探目标也愈加精细，探测要求不再局限于目标体位置、埋深的粗略估计，而需要能够直接给出工程勘探中更为关心的目标材质、大小、尺寸等定量指标。需要从大量的雷达数据中，得到符合工程勘探要求的精确物性参数指标，从而挖掘更多的有效信息。GPR数据处理手段由传统的定性分析处理逐步变为更精确的定量分析，进一步从观测数据中推断出物性参数。而反演本质上是选取合适的物性参数，能最优地拟合地球物理观测数据。反演过程中需要通过正演建模得到观测数据，因此正演是反演的基础。研究能够处理任意复杂地质结构的高效、高精度的正演求解器也成为GPR反演的难点和热点问题。
　　Maxwell方程是描述与刻画电磁波波动现象、传播规律与物理过程的有力工具，也是GPR波动方程正反演探求复杂探测对象奥秘的有效手段。目前计算电磁学中计算方法众多，从求解方程的形式上可以分为积分方程（integral equation，IE）类方法和微分方程（partial differential equation，PDE）类方法，IE类方法具有很高的精度，但是难以处理复杂介质的电磁问题，因此应用范围受限，考虑到GPR的探测对象特点，本书中主要介绍PDE方法。从计算域进行分类，探地雷达正演可以分为时间域方法和频率域方法，两者的计算结果可由傅里叶变换进行转换，比较而言，由于时间域方法结果更为直观，故使用更为广泛。然而，时域和频域的Maxwell方程是理解电磁现象的互补工具：时域方程对于研究瞬态和动力学是必不可少的，频域方程对于研究稳态和精确处理分散材料也至关重要。因此，频率域方法适用于单频问题，时间域方法在处理宽频带、时变介质时具有更大优势，它们各具优缺点。GPR正演模拟的常用时域数值方法主要包括时域有限差分法（finitedifference timedomain，FDTD）和时域有限单元法（finite element time domain，FETD）。FDTD方法是计算电磁学中最流行的方法，其原理简单、思路清晰、空间和时间离散处理方案容易，易于并行实现，然而FDTD方法建模时，由于不能与非结构化网格结合，对不规则目标体、复杂的物性分界面拟合效果不好，存在台阶化误差，且时间步长受稳定性条件的限制，为保证结构复杂的细微异常体的模拟精度与解的稳定性，时间步长必须足够小。FETD能够与非结构化网格结合，空间离散更加灵活，可以较好地拟合目标体的形状，但该方法需要在每一个时间步上求解一次组装后的大型方程组，计算复杂度高，并行计算实现困难（冯德山等，2017）。为了处理电磁问题时既有FETD方法空间离散的灵活性，同时又能够像FDTD方法一样采用显式迭代的计算方案节省计算资源，通过数值通量概念的引入，本书将时域间断伽辽金法（discontinuous Galerkin time domain，DGTD）引入到GPR正演中。该方法在近十年来得到了迅速的发展和应用，在计算流体力学、电磁学、地球物理学等学科的PDE求解中展现了强有力的活力，得到了高度关注（Hesthaven and Warburton，2011；葛德彪和魏兵，2019）。相对于其他几种成熟的时域计算方法，DGTD算法在GPR正演领域具有较好的优势，但目前可查阅的文献资料较少，仍然不成熟。与时间域类似，频率域正演方法主要有有限差分法（finitedifference method，FDM）与有限单元法（finite element method，FEM）；由于GPR正演是宽频带电磁问题，频率域方法在GPR中应用较少，然而考虑到频率域GPR全波形反演算法具备的效率优势，研究合适的频率域GPR正演算法势在必行。
　　探地雷达的反演是根据地面或者井中接收天线的观测数据推断地下介质的物性参数，如介电常数、电导率、电磁波波速等。从数学的角度来讲，反演是指在选定的正演数学物理模型的情况下，建立观测数据与正演数据的误差泛函，然后利用寻优算法求解误差泛函的极小值，得到介质物性参数的分布，因此，正演数学物理模型的选取对反演至关重要。GPR的全波形反演（full waveform inversion，FWI）以电磁波PDE方程预测数据，充分利用电磁波形中的走时、相位与振幅等信息，具有揭示复杂地质背景下构造与储层物性的潜力，其理论分辨率可以达到1/3～1/2波长，是重构复杂地质结构、异常体的有力工具并已成为当前研究的热点。GPR全波形反演可以视为一种非线性优化问题，从实现途径上分类主要包括时间域、频率域以及拉普拉斯域。时间域全波形反演便于对数据进行各种预处理，能更好地适应多节点并行计算，由于采用时间域的雷达记录作为观测数据，数据量巨大，需要巨大的存储空间；频率域全波形反演仅使用几个离散频率的采样点信息，节省了计算的存储资源，多激励源计算效率高，在预处理上受限；拉普拉斯域方法反演精度不高，其结果适合作为全波形反演的初始模型。从最优化方法上一般分为两类：全局最优化算法和局部最优化算法。全局最优化算法也称为完全非线性方法，包括模拟退火法、遗传算法、粒子群算法以及神经网络算法等。全局算法理论上可以解决GPR全波形反演问题，并且在实现上不依赖于初始模型，然而存在着早熟收敛、计算量巨大、计算结果不够稳定等缺点，在二、三维反演中，参数空间的大小限制了该类方法在GPR全波形反演中的应用。局部优化算法也称为拟线性反演方法，包括牛顿法、高斯牛顿法、非线性共轭梯度法、拟牛顿法、最速下降法等。局部优化算法主要通过局部搜索算法解决全波形反演问题，其优点是简单、灵活且易于实现，但它们都是局部收敛方法，迭代结果明显依赖于初始猜测的选取，容易陷入局部极值，很难捕获到满意的全局最优解。
　　目前地面探地雷达理论模型及小数据量雷达剖面的全波形反演研究已经初有成效，但是由于该技术计算效率上的局限性以及实际资料的复杂性，仍局限于实验室阶段。受计算成本大、目标函数非线性强、对初始模型的依赖性强以及实测数据子波未知等因素的影响，全波形反演精度和效率达不到实际应用的要求。因此，寻求高效、简洁、实用的GPR正反演方法，就成了当前计算地球物理学领域研究的迫切任务。本书将重点介绍探地雷达正演模拟及全波形反演新算法，可为今后开展高效率、高精度的三维GPR反演技术以及资料解释奠定扎实基础。
　　1.2 探地雷达正演模拟的研究现状及趋势
　　GPR正演模拟是通过已知地下介质或物体内部参数的分布，运用数学手段模拟电磁波在其中的传播规律或特征。从电磁学角度而言，GPR正演模拟属于计算电磁学的一种应用方向。它作为一种深入认知地球内部或人工构筑物的手段，得到了迅速发展。下面从时间域和频率域分别介绍GPR正演研究现状。
　　1.2.1 时间域正演模拟方法
　　1.FDTD研究现状
　　在时间域模拟方法中，FDTD具有直接时域计算、思路清晰、编程简单等特点，得到了广泛的应用。1966年，Yee（1966）首次提出电场、磁场分量在时间空间上交错采样计算方法，并成功地模拟了电磁脉冲与理想导体作用的电磁问题；Taflove（1980）于1980年采用该技术实现了三维金属腔的电磁入射模拟，并首次给出了FDTD的缩写，并于1995年出版专著详细地介绍了FDTD算法（Taflove，1995），国内王长清和祝西里（1994）、高本庆（1995）及葛德彪和闫玉波（2005）分别出版了电磁波FDTD算法原理介绍与应用研究的专著。FDTD也是GPR模拟中应用最广泛的一种方法，国内外有不少研究成果。研究方法主要集中于算法实现、复杂介质建模、吸收边界、天线特性和实际应用模型分析等方面。其中算法实现上：冯德山（2003）采用C语言编制了二维FDTD的探地雷达程序，并给出了所有源代码；Giannopoulos（2005）编制了GPR二、三维正演模拟软件“gprMax”，可进行规则雷达模型正演；Irving和Knight（2006）采用Matlab语言开发了基于卷积完全匹配层的二维GPR正演程序，使用TM/TE模式分别模拟了地面GPR数据和钻孔雷达数据；Warren等（2016，2019）采用Python语言开发了开源的gprMax3.0版本，并采用CUDA对FDTD算法进行加速；冯德山等（2017）采用Matlab语言开发了二、三维FDTD算法GPR正演程序。针对复杂介质：Xu和McMechan（1997）、方广有等（1998）运用2.5维的FDTD算法对Debye型色散介质的地电模型进行了脉冲GPR的模拟计算；Bergmann等（1998）采用时间2阶空间4阶FDTD算法实现了二维色散、衰减介质的GPR正演模拟；Teixeira等（1998）开展了三维洛伦兹色散介质GPR的FDTD正演模拟；Cassidy和Millington（2009）应用FDTD算法实现了三维磁损耗介质的GPR正演模拟；刘四新和曾昭发（2007）用含传导电流项的Maxwell方程开展三维Debye型频散介质的GPR正演；Giannakis和Giannopoulos（2014）提出了一种基于FDTD的色散介质分段线性递归卷积方法；Jiang等（2013）、郭士礼（2013）、李静（2014）采用FDTD算法分析了随机介质模型中GPR信号的响应特征。吸收边界条件中：岳建华和何兵寿（1999）实现了基于Mur吸收边界的FDTD GPR正演模拟；刘四新等（2005）将Berenger PML边界条件（Berenger，1994）应用到了FDTD算法的GPR模拟中；Giannopoulos（2008）提出了一种改进的FDTD复频移PML（complex frequency shifted PML，CFS-PML）实现方法；冯德山等（2010）、李静等（2010）分别将各向异性PML（uniaxial anisotropy PML，UPML）边界条件（Gedney，1996）应用到GPR数值模拟中；Li等（2012）采用基于CFS-PML边界条件的高阶FDTD算法进行了三维GPR数值模拟；Giannopoulos（2018）提出了一种基于多极完美匹配层有限差分时域电磁建模算法；Zhang B等（2019）提出了一种新的非分裂域卷积完全匹配层方案，该方法采用优化项和调整因子来寻求最优本构系数。天线性能分析中：Jiao等（2000）、Kwan等（2004）、刘立业等（2006）、Shaari等（2010）采用FDTD法开展了GPR天线特性分析研究；冯晅等（2011）应用FDTD算法研究了全极化GPR模型对不同目标的响应特征；Warren和Giannopoulos（2011）应用Taguchi优化方法建立了商用探地雷达天线FDTD模型；王芳芳和张业荣（2010）应用FDTD算法开展了超宽带（ultra wide band，UWB）雷达的数值模拟；田钢等（2011）使用三维FDTD算法分析了地面以上物体对GPR信号的反射干扰特征；Diamanti等（2008）用FDTD数值模拟和试验验证了GPR在砖混拱桥中环形分离的可行性；杨峰等（2008）、张鸿飞等（2009）、刘新荣等（2010）、肖建平等（2017）基于FDTD算法建立了衬砌病害模型，分析了衬砌病害的GPR响应特征；Scheers等（2000）、Giannakis等（2015）通过FDTD算法