基于插值尺度函数法的探地雷达三维正演方法技术

本发明专利技术公开了一种基于插值尺度函数法的探地雷达三维正演方法，包括设置三维地质模型并给出模型参数分布；设定三维正演方法的参数；设置CFS‑PML边界条件；从新测点开始进行正演计算；对设置的边界条件进行空间离散和时间离散；计算正演参数并对场值进行更新；重复上述步骤直至所有的测点均完成正演计算，保存所有数据并完成探地雷达的三维正演。本发明专利技术方法较常规FDTD算法具有更高的精度和效率；而且较FDTD算法具有更好的色散特性，能够使RK4‑ISFM实际计算中具有较好的数值色散特性，节省计算机运行时间和内存。

3-D Forward Modeling of Ground Penetrating Radar Based on Interpolated Scale Function Method

【技术实现步骤摘要】
基于插值尺度函数法的探地雷达三维正演方法
本专利技术具体涉及一种基于插值尺度函数法的探地雷达三维正演方法。
技术介绍
探地雷达已成为工程检测和环境地球物理勘探领域中的一种重要浅层方法。鉴于日益精细化的勘探要求，需要给出工程中最关心的异常体大小、形态和位置等准确量化指标。目前二维雷达探测因为勘探成本低、探测效率高而被广泛采用，但二维探测能探测地质体的存在与否，但难于区分旁侧异常体的干扰。三维雷达探测却能更加全面而细致地了解探测对象的准确信息，被越来越广泛地应用。通过对已知地质模型的三维正演，可以丰富雷达模型数据库，熟悉并了解典型地质体的雷达图像回波特征，指导雷达探测图谱解释的可靠性和准确性的提高。目前，诸多学者已经成功地使用不同算法开展了GPR数值模拟研究：Bergmann等、Zhan应用时域有限差分法(FiniteDifferenceTimeDomain,FDTD)进行了GPR数值模拟，Diamanti、冯德山对算法改进、边界条件等内容进行了深入的研究；底青云&王妙月、冯德山等应用有限单元法开展了二维GPR正演；Liu等将时域伪谱法应用于频散介质的GPR模拟，与FDTD算法相比，该方法需要更少的未知数，适于求解多维大尺度问题；Tiao等应用时域间断有限元(DiscontinuousGalerkinMethod,DGM)进行Debye型色散介质中的GPR数值模拟；方宏远等利用Hamilton系统的辛分块龙格库塔方法开展了复杂地下结构模型的雷达剖面特征及波传播规律研究；冯德山等将样条小波有限元引入到探雷达数值模拟中。Liu等对典型的不利地质条件，如溶洞、断层、裂隙岩体、裂隙网络、含水体等进行了探地雷达正演模拟。Solla等通过FDTD算法对GPR模型进行重构，提出了对历史桥梁的评价，利用得到的数据来解释桥梁内部不同材料的组成和分布。在所有的GPR正演算法中，时域有限差分法具有直接时域计算、编程容易、节约计算时间等优点而备受青睐。Irving&Knight基于MATLAB开发了GPR二维(2D)FDTD正演程序。Giannopoulos编制了GPR正演软件“GprMax”，得到许多用户青睐，进一步推广了FDTD在GPR数值模拟中的应用。但上述文献大多局限于2D探地雷达数值模拟，三维(3D)探地雷达数值模拟由于受计算量大的限制，仍不成熟，处于发展阶段。Wang&Tripp应用FDTD算法求解了3D介质中的电磁波传问题，提出了一种新的优化二阶差分方案来近似空间导数，有效地节约计算时间和内存需求。Champagne等应用Fortran90编写了三维频率域有限差分法电磁感应成像程序。李静等将高阶FDTD法用于3D探地雷达数值模拟，截断边界处采用了不分裂复频移完全匹配层(ComplexFrequencyShiftedPerfectlyMatchedLayer,CFS-PML)，有效地提高了计算精度与掠角波、低频波、隐失波的吸收效果。但是目前对于三维(3D)探地雷达数值模拟，依然存在精度不高和三维正演效率较低的问题。
技术实现思路
本专利技术的目的在于提供一种可靠性高、计算精度高且正演效率较高的基于插值尺度函数法的探地雷达三维正演方法。本专利技术提供的这种基于插值尺度函数法的探地雷达三维正演方法，包括如下步骤：S1.设置三维地质模型，并给出模型的参数分布；S2.设定三维正演方法的参数；S3.设置CFS-PML边界条件；S4.从新测点开始，进行正演计算；S5.对步骤S3设置的边界条件进行空间离散；S6.对步骤S3设置的边界条件进行时间离散；S7.计算正演参数，并对场值进行更新；S8.重复步骤S4～S7直至所有的测点均完成正演计算后，保存所有数据，从而完成探地雷达的三维正演。步骤S2所述的三维正演方法的参数，具体包括设定天线位置、网格空间步长和视窗长度，并加载激励源。步骤S3所述的设置CFS-PML边界条件，具体为采用如下步骤设置边界条件：采用辅助微分方程法来实现该边界条件；伸缩坐标中三维Maxwell方程公式：式中si为坐标伸缩因子，si＝κi+σi/(αi+jωε0)αi和σi为正实数，κi为大于1的实数，i＝x,y,z；对si进行变换得到：引入辅助变量：其中为与空间偏导数相关的辅助变量，为与空间偏导数相关的辅助变量；从而得到：由频率域变回时间域后得到：对辅助变量定义公式进行分解可得：并由频率域变回时间域后得：步骤S5所述的对步骤S3设置的边界条件进行空间离散，具体为采用ISFM的导数对Maxwell方程进行空间离散：规定函数u在空间离散，记为(-i-1/2)Δx且i∈N，N为阶数；则u的空间导数为式中ui+1/2＝u((i+1/2)Δx)；则空间导数为：式中l0∈N，a(l)为ISFM的导数系数。步骤S6所述的对步骤S3设置的边界条件进行时间离散，具体为采用RK4对Maxwell方程进行时间离散：递推方程中的时间偏导数采用RK4格式求解；RK4格式如下：式中L为算子且w＝Hx，n表示时间步，h(i)(i＝1,2,3,4)表示RK4的每一步函数值。本专利技术提供的这种基于插值尺度函数法的探地雷达三维正演方法，提出了一种新的探地雷达三维数值模拟新算法—RK4-ISFM，该算法将插值小波尺度函数的导数代替中心差分格式来离散Maxwell方程组的空间导数，并采用四阶龙格库塔算法计算时间导数，由于插值小波尺度函数及四阶龙格库塔算法有效提高了GPR波动方程的空间和时间离散精度，因此，插值小波尺度算法较常规FDTD算法具有更高的精度和效率；而且尽管插值小波尺度算法的Courant稳定性条件强化了对时间步长的限制，但在选取同样时间步长条件下，它可在一个波长内取8个网格而不产生数值频散，较FDTD算法具有更好的色散特性；该特性使RK4-ISFM实际计算中，即使在网格采样点数较少，空间步长较大的情况下，仍具有较好的数值色散特性，节省计算机运行时间和内存。附图说明图1为本专利技术方法的方法流程示意图。图2为本专利技术方法的2～5阶次插值尺度函数示意图。图3为本专利技术方法的θ＝π/2、时，绘出的FDTD及ISFM两种算法的相对相速与网格数变化关系示意图。图4为本专利技术方法的包含三种异常体下的三维雷达模型示意图。图5为本专利技术方法的ISFM算法得到的三维正向切片视图。图6为本专利技术方法的x＝0.25m，0.45m，0.65m，0.85m，1.05m五个位置处的3D剖面示意图。图7为本专利技术方法的x＝0.65m，y＝0.3m处的3D雷达切片示意图。图8为本专利技术方法的去除直达波场与背景波场后提取异常成果示意图。图9为本专利技术方法的在真空介质中，不同网格数下，FDTD及ISFM两种算法的单道雷达波形与误差对比示意图。具体实施方式如图1所示为本专利技术方法的方法流程示意图：本专利技术提供的这种基于插值尺度函数法的探地雷达三维正演方法，包括如下步骤：S1.设置三维地质模型，并给出模型的参数分布；S2.设定三维正演方法的参数；具体包括设定天线位置、网格空间步长和视窗长度，并加载激励源；S3.设置CFS-PML边界条件；具体为采用如下步骤设置边界条件：采用辅助微分方程法来实现该边界条件；伸缩坐标中三维Maxwell方程公式：式中si为坐标伸缩因子，si＝κi+σi/本文档来自技高网...

【技术保护点】
1.一种基于插值尺度函数法的探地雷达三维正演方法，包括如下步骤：S1.设置三维地质模型，并给出模型的参数分布；S2.设定三维正演方法的参数；S3.设置CFS‑PML边界条件；S4.从新测点开始，进行正演计算；S5.对步骤S3设置的边界条件进行空间离散；S6.对步骤S3设置的边界条件进行时间离散；S7.计算正演参数，并对场值进行更新；S8.重复步骤S4～S7直至所有的测点均完成正演计算后，保存所有数据，从而完成探地雷达的三维正演。

【技术特征摘要】
1.一种基于插值尺度函数法的探地雷达三维正演方法，包括如下步骤：S1.设置三维地质模型，并给出模型的参数分布；S2.设定三维正演方法的参数；S3.设置CFS-PML边界条件；S4.从新测点开始，进行正演计算；S5.对步骤S3设置的边界条件进行空间离散；S6.对步骤S3设置的边界条件进行时间离散；S7.计算正演参数，并对场值进行更新；S8.重复步骤S4～S7直至所有的测点均完成正演计算后，保存所有数据，从而完成探地雷达的三维正演。2.根据权利要求1所述的基于插值尺度函数法的探地雷达三维正演方法，其特征在于步骤S2所述的三维正演方法的参数，具体包括设定天线位置、网格空间步长和视窗长度，并加载激励源。3.根据权利要求1或2所述的基于插值尺度函数法的探地雷达三维正演方法，其特征在于步骤S3所述的设置CFS-PML边界条件，具体为采用如下步骤设置边界条件：采用辅助微分方程法来实现该边界条件；伸缩坐标中三维Maxwell方程公式：式中si为坐标伸缩因子，si＝κi+σi/(αi+jωε0)αi和σi为正实数，κi为大...

【专利技术属性】
技术研发人员：冯德山，王珣，杨道学，王向宇，曹岑，
申请(专利权)人：中南大学，
类型：发明
国别省市：湖南,43