本发明专利技术公开了一种基于间断有限元的探地雷达三维正演方法，包括构建用于地质勘探的正演模型并对正演模型进行网格剖分；基于网格剖分结果对电磁波在正演模型中传播的麦克斯韦方程组进行求解，求解时引入基于间断有限元算法，采用数值通量实现相邻单元之间的信息交换，同时定义插值基函数，并通过显示时间离散进行时间递推，得到正演模型中的电场强度和磁场强度。该方法具有提高正演精度、计算内存需求少、易于并行计算，以及更适用于起伏界面的正演等技术效果。
【技术实现步骤摘要】
本专利技术属于地理勘探，具体涉及一种基于间断有限元的探地雷达三维正演方法。
技术介绍
1、探地雷达作为一种能对地质目标进行快速高效探测的地球物理方法，是目前工程和环境地球物理领域的研究热点之一，在工程勘查、煤矿开采、地下管线探测和隧道质量检测领域获得广泛应用。目前，探地雷达成像技术已较为成熟。然而，现有的一维和二维、三维算法尚不能满足复杂地下构造的模拟需求，限制了针对复杂模型的探地雷达数据解释。
2、目前探地雷达中最常用的数值模拟算法包括：时域有限差分法(fdtd)(李静等，2010；崔凡等，2022)、有限体积法、有限元法(fem)等(冯德山等，2018)。时域有限差分法是一种计算较为简单的方法。该方法采用差分逼近空间导数，得到显式的半离散形式可以灵活地选取不同的时间离散方法，其优势是可以很好地处理介质内部电磁场的不连续性、易于编程。有限体积法是一种将计算区域划分为一系列不重复的控制体积，通过在每个控制体积上对待解微分方程进行积分得到离散方程的方法，计算区域内场的守恒可以由每个控制体积内场变量的守恒来保证。有限元法可以灵活地选取六面体或非规则四面体网格，因而能很好地拟合复杂边界。该方法通过增加单元自由度可相对容易地构造高阶格式，并且可以在不同的单元上采取不同的阶数，目前各种自适应有限元算法已经非常成熟，可以实现网格和阶数的自适应。然而，这种传统的有限元算法需要组装和求解大型线性方程组，因而会占用大量计算资源、限制计算效率。采用一种高精度高效率并节省内存的正演算法对提高探地雷达正演、反演、数据解释都至关重要。
3、探地雷达中主流的正演算法都有各自的优势，但每种方法本身都存在无法避免的缺陷。时域有限差分法采用局部一维多项式进行逼近，导致其只能在一维基础上求解高维问题。进而，网格结构的限制使其不能很好地处理含复杂模型边界。对于有限体积法来说，想要通过构造高阶格式以提高算法精度，容易导致方程变得异常复杂且不稳定，这会破坏算法的几何灵活性。有限元算法在求解过程中，需要组装和求解大型线性方程组，会占用大量计算资源、限制计算效率。
技术实现思路
1、鉴于上述，本专利技术的目的是提供一种基于间断有限元的探地雷达三维正演方法，具有提高正演精度、计算内存需求少、易于并行计算，以及更适用于起伏界面的正演等技术效果。
2、为实现上述专利技术目的，本专利技术实施例提供的一种基于间断有限元的探地雷达三维正演方法，包括以下步骤：
3、构建用于地质勘探的正演模型并对正演模型进行网格剖分；
4、基于网格剖分结果对电磁波在正演模型中传播的麦克斯韦方程组进行求解，求解时引入基于间断有限元算法，采用数值通量实现相邻单元之间的信息交换，同时定义插值基函数，并通过显示时间离散进行时间递推，得到正演模型中的电场强度和磁场强度。
5、优选地，对正演模型进行网格剖分时，采用非结构四面体进行空间离散，即将正演模型对应的计算域ω剖分为k个互不重叠的四面体单元dk，k为四面体单元索引，取值为k＝1…k，并假设在每个四面体单元中电导率σ、磁导率μ、介电常数ε均为常数。
6、优选地，对麦克斯韦方程组进行求解时，假设第k个四面体单元dk的残差加权积分为0，即
7、
8、
9、其中，w为加权函数，e、h分别为电场和磁场强度，js为外加源电流密度，t为时间，为梯度算子，dv为对四面体单元中的体积v进行积分。
10、优选地，引入基于间断有限元算法求解时，结合矢量公式对(1)和(2)式中的旋度项进行变换，并进一步整理得到间断有限元算法伽辽金弱形式为：
11、
12、
13、其中，n为单元界面外法向向量，ds为对单元界面的面积s进行积分；
14、根据间断有限元法基本原理，需要采用数值通量(e*，h*)实现相邻单元之间的信息交换，利用数值通量(e*，h*)代替(3)和(4)式中的面积分项，得到：
15、
16、
17、优选地，采用mercerat和glinsky给出的数值通量(e*，h*)的格式，即：
18、
19、
20、其中，z、y分别为当前四面体单元dk的阻抗与导纳z+、y+分别为在单元交界面上与dk相邻的四面体单元的阻抗与导纳，e+、h+分别为与dk相邻的四面体单元上的电场值与磁场值，为单元dk上的电场值。
21、优选地，引入基于间断有限元算法时，对(5)和(6)式中的电磁场进行空间离散，则单个四面体单元内任意位置的电磁场在x、y、z三个方向的分量表示为：
22、
23、
24、其中，np为单元上插值节点个数，其取值与插值基函数的阶数相关，j为单元上插值节点索引，i＝1、2、3分别对应为x、y、z三个方向，r为任一点的位置矢量，ei(r,t)为t时刻第i个方向的电场分量，ei,j(t)为第j个节点在第i个方向的电场值，lj(r)为任一点处的插值基函数，ei＝[ei1,ei2,…,ei,np]t和hi＝[hi1,hi2,…,hi,np]t分别是单元插值节点上的电场值和磁场值，l＝[l1,l2,…,lnp]t为单元插值基函数，hi(r,t)为t时刻第i个方向的磁场分量，hi,j(t)ei,j(t)为第j个节点在第i个方向的磁场值，上标t表示转置；
25、将(9)和(10)式中带入(5)和(6)式，并将插值基函数设为加权函数w，得到：
26、
27、
28、
29、
30、
31、
32、其中，m为质量矩阵，sx、sy、sz为刚度矩阵，jx、jy、jz为外加源项，fxe、fye、fze、fxh、fyh、fzh为离散数值通量积分。
33、优选地，所述插值基函数采用定义在标准四面体上的层级型正交基函数，首先在标准四面体ωs＝{(r,s,t)|r,s,t≥-1,r+s+t≤-1}上定义n阶proriol-koornwinder-dubiner多项式，表示为：
34、
35、其中，m＝1,2,…,np代表np个节点,pi(α,β)为i阶jacobi多项式，当α＝β＝0时pi(α,β)为legendre多项式，c＝t，r,s,t为标准四面体上节点的坐标；
36、然后根据多项式ψm(r,s,t)构造范德蒙德矩阵v：
37、
38、其中，r1，r2，...rnp为插值节点坐标，ψnp(rnp)为插值节点rnp处的多项式；
39、最后根据范德蒙德矩阵v构造插值基函数l＝[l1，…，lnp]。
40、优选地，引入基于间断有限元算法求解时，显示时间离散算法采用显式低存储五级四阶runge-kutta，即采用显式低存储五级四阶runge-kutta进行时间递推，包括：
41、首先，将(11)-(16)式改写成(19)式标本文档来自技高网
...
【技术保护点】
1.一种基于间断有限元的探地雷达三维正演方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的基于间断有限元的探地雷达三维正演方法，其特征在于，对正演模型进行网格剖分时，采用非结构四面体进行空间离散，即将正演模型对应的计算域Ω剖分为K个互不重叠的四面体单元Dk，k为四面体单元索引，取值为k＝1…K，并假设在每个四面体单元中电导率σ、磁导率μ、介电常数ε均为常数。
3.根据权利要求2所述的基于间断有限元的探地雷达三维正演方法，其特征在于，对麦克斯韦方程组进行求解时，假设第k个四面体单元Dk的残差加权积分为0，即
4.根据权利要求3所述的基于间断有限元的探地雷达三维正演方法，其特征在于，引入基于间断有限元算法求解时，结合矢量公式对(1)和(2)式中的旋度项进行变换，并进一步整理得到间断有限元算法伽辽金弱形式为：
5.根据权利要求4所述的基于间断有限元的探地雷达三维正演方法，其特征在于，采用Mercerat和Glinsky给出的数值通量(E*，H*)的格式，即：
6.根据权利要求4所述的基于间断有限元的探地雷达三维正演方法，其特征在于，引入基于间断有限元算法时，对(5)和(6)式中的电磁场进行空间离散，则单个四面体单元内任意位置的电磁场在x、y、z三个方向的分量表示为：
7.根据权利要求6所述的基于间断有限元的探地雷达三维正演方法，其特征在于，所述插值基函数采用定义在标准四面体上的层级型正交基函数，首先在标准四面体Ωs＝{(r,s,t)|r,s,t≥-1,r+s+t≤-1}上定义N阶Proriol-Koornwinder-Dubiner多项式，表示为：
8.根据权利要求6所述的基于间断有限元的探地雷达三维正演方法，其特征在于，引入基于间断有限元算法求解时，显示时间离散算法采用显式低存储五级四阶Runge-Kutta，即采用显式低存储五级四阶Runge-Kutta进行时间递推，包括：
...
【技术特征摘要】
1.一种基于间断有限元的探地雷达三维正演方法，其特征在于，包括以下步骤：
2.根据权利要求1所述的基于间断有限元的探地雷达三维正演方法，其特征在于，对正演模型进行网格剖分时，采用非结构四面体进行空间离散，即将正演模型对应的计算域ω剖分为k个互不重叠的四面体单元dk，k为四面体单元索引，取值为k＝1…k，并假设在每个四面体单元中电导率σ、磁导率μ、介电常数ε均为常数。
3.根据权利要求2所述的基于间断有限元的探地雷达三维正演方法，其特征在于，对麦克斯韦方程组进行求解时，假设第k个四面体单元dk的残差加权积分为0，即
4.根据权利要求3所述的基于间断有限元的探地雷达三维正演方法，其特征在于，引入基于间断有限元算法求解时，结合矢量公式对(1)和(2)式中的旋度项进行变换，并进一步整理得到间断有限元算法伽辽金弱形式为：
5.根据权利要求4所述的基于间断有限元的探地雷达三维正演方法，其特征在于，采...
【专利技术属性】
技术研发人员：薛舒杨，殷长春，刘午扬，李俊，黄贤阳，
申请(专利权)人：吉林大学，
类型：发明
国别省市：
还没有人留言评论。发表了对其他浏览者有用的留言会获得科技券。