一种基于无网格法的线源时域电磁响应数值计算方法技术

本发明专利技术涉及一种基于无网格法的线源时域电磁响应数值计算方法，尤其是可以克服传统数值计算方法中对于网格的依赖，适用于复杂地形下时域电磁探测的数值模拟。本发明专利技术基于瞬变电磁法满足的控制方程和定解条件，建立了二维线源边值问题的泛函，利用罚因子法加载本质边界条件，提出旁轴近似方程消除截断边界处的反射波，采用Crack-Nicolson格式进行时间离散，得到递推方程。利用等参单元思想将局部坐标中形状规则的单元离散为节点任意分布的不规则求解对象。采用LU分解方法求解递推方程，最终得到求解区域内各个节点的场值。计算结果表明，该方法形函数光滑性好，模拟精度高，最大误差不超过1×10-3，实现了电磁法高精度的数值计算。

【技术实现步骤摘要】

本专利技术涉及一种地球物理勘探领域的电磁场数值计算方法，尤其是基于无网格法的线源时域电磁响应数值计算方法。
技术介绍
时间域电磁法(Timedomainelectromagneticmethods)或称瞬变电磁法(Transientelectromagneticmethods)，是一种建立在电磁感应原理基础上的人工源电磁探测方法。目前该方法被广泛应用于金属矿勘探、煤矿水文地质调查和工程勘查等领域。电磁场数值模拟是数据处理和反演解释的基础，在地球物理数值计算中具有基础性、全局性作用。现有瞬变电磁数值计算方法主要有积分方程法(IEM)，有限差分法(FDM)，有限单元法(FEM)等。CN201410125866.5公开了一种频率域正演方法及装置，该方法通过建立17点格式的差分公式，构造稀疏矩阵，读入子波参数和速度模型，频率循环得到单频波场，对所有频率波场求反傅里叶变换，得到正演结果。该方法增加了差分公式的适用性，通过加载PML(最佳匹配层)边界条件，防止边界反射的波场回折到地表。CN201110459223.0公开了一种电磁场仿真分析方法，该方法采用棱边单元对三维模型进行网格剖分，引入满足inf-sub条件的标量乘子空间，加入有约束条件的集成后总体矩阵方程求解，得到电磁场分析结果。美国专利US11756384公开了一种基于三维时域正演建模基础上的反演方法。该方法基于梯度的波形反演理论，通过匹配建模的数据来估计模型参数，在时域建模之后通过离散傅里叶变换(FFT)转换到频率域，完成了三维频域波形反演。以上所述方法均是基于网格实现的，在计算的过程中都不可避免地涉及到网格的剖分问题，对于单元的形状也有一定的要求，特别是对于物性参数分布复杂和几何特征分布不规则的地电模型适应性差。
技术实现思路
本专利技术所要解决的技术问题在于提供一种基于无网格法的线源时域电磁响应数值计算方法，可以克服传统正演方法中对于网格的依赖，适用于复杂地形下时域电磁探测的数值模拟。本专利技术是这样实现的，一种基于无网格法的线源时域电磁响应数值计算方法包括：1)在计算区域内进行节点、电性参数、支持域、背景网格、初始场值设置；2)在计算区域内对背景网格进行等参数变换，在每个背景单元中分别计算高斯积分点及其权系数；3)对所有背景网格进行外循环，对所有高斯积分点内循环，搜索背景网格中局部定义域内的有效节点并计算局部定义域内节点处形函数；4)加载强加边界条件得到线性方程组，加载吸收边界条件得到边界阻尼矩阵，将边界阻尼矩阵加入到线性方程组中；5)采用LU分解法求解线性方程组，得到各个节点各个时刻的场值。进一步地，步骤3中局部定义域内利用滑动最小二乘法建立形函数。进一步地，步骤4中利用罚因子法处理计算模型的强加边界条件，建立和瞬变电磁边值问题的等价泛函数，求取泛函数的变分，得到刚度矩阵K，阻尼矩阵K'和右端项S，形成线性方程组：K′E·+KE=S]]>进一步地，步骤4中根据旁轴近似方程加载吸收边界条件，将左行波作用于右行入射的波，右行波作用于左行入射的波，下行波作用于上行入射的波，相互抵消后消除截断边界处的反射波，得到边界阻尼矩阵：F=1v×(∫ΓlN1N2(nx+ny)dΓ+∫ΓdN1N2(nx+ny)dΓ+∫ΓrN1N2(nx+ny)dΓ)]]>N1和N2均为无网格法边界处两个节点的形函数，nx和ny为边界外法线的方向余弦，v为地下介质中的电磁波传播速度。进一步地，步骤5中线性方程组的时间变量的离散采用Crack-Nicolson格式，令K′+F＝K”，最终形成的递推方程为：(2K”+ΔtKt+Δt)Et+Δt＝ΔtSt+Δt+StΔt+2K”Et-ΔtKtEt其中Et为当前时刻的场值，Et+Δt为待求的下一时刻的场值，S为线性方程组的右端项，Δt为离散的时间步长，K为刚度矩阵。本专利技术与现有技术相比，有益效果在于：本专利技术针对传统的网格化数值计算方式的不足，提出了无网格法模拟的方案。以二维线源瞬变电磁法为例，利用变分原理详细推导了等价线性方程组的形成过程，采用Crack-Nicolson方案离散时间变量，形成最终的递推方程，求取了线源各个节点各个时刻的电磁响应，完整地展现了无网格法电磁数值模拟的全部过程。由于无网格法摆脱了传统方法中对于网格的依赖，理论上节点可以任意布置，为含有起伏地形等复杂情况下电磁法的数值模拟提供新的参考。在实践中同时发现无网格法形函数光滑性好，模拟精度高，为实现电磁法高精度反演打下了坚实的基础。附图说明图1是无网格法高斯点、场节点、影响域，局部定义域与背景单元示意图；图2是二维线源时域电磁探测响应计算模拟示意图；图3是截断边界处的旁轴近似方程示意图；图4是任意四边形下的等参数变换示意图；图5是二维线源时域电磁响应计算方法流程图；图6是本专利技术一个实施例节点处的无网格法数值解与理论解对比图；图7是本专利技术一个实施例一条剖面中无网格法数值解与理论解对比图。具体实施方式为了使本专利技术的目的、技术方案及优点更加清楚明白，以下结合实施例，对本专利技术进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本专利技术，并不用于限定本专利技术。实施例参见图5结合图1所示，一种基于无网格法的线源时域电磁响应数值计算方法，包括：1)设置模型参数，包括节点设置、电性参数设置、背景网格设置、初始场值设置；如图1所示的无网格法高斯点、场节点、影响域ΩI，局部定义域Ωx与背景单元示意图，包括背景积分网络1中的求解域Ω，在一个背景积分网络1的求解域Ω内包括用黑点表示的节点3，以及用空心圆表示的高斯积分点或计算点2，影响域ΩI是指一个节点的影响区域，影响域ΩI的中心为场节点，而局部定义域Ωx则为计算点进行最小二乘拟合所选择的区域，中心常常为一个计算点。在实际计算中常用节点的影响域取代计算点的局部定义域，即在构造一个计算点的无网格形函数时，若计算点位于场节点的影响域中，则该场节点属于构造该形函数的场节点，也即如果场节点的影响域包含该计算点，则该场节点参与构造该计算点的形函数。2)在计算区域内对背景网格进行等参数变换，在每个背景单元中分别计算高斯积分点及其权系数；3)对所有背景网格进行循环；4)对本文档来自技高网...

【技术保护点】
一种基于无网格法的线源时域电磁响应数值计算方法，其特征在于：包括如下的步骤：1)在计算区域内进行节点、电性参数、支持域、背景网格、初始场值设置；2)在计算区域内对背景网格进行等参数变换，在每个背景单元中分别计算高斯积分点及其权系数；3)对所有背景网格进行外循环，对所有高斯积分点内循环，搜索背景网格中局部定义域内的有效节点并计算局部定义域内节点处形函数；4)加载强加边界条件得到线性方程组，加载吸收边界条件得到边界阻尼矩阵，将边界阻尼矩阵加入到线性方程组中；5)采用LU分解法求解线性方程组，得到各个节点各个时刻的场值。

【技术特征摘要】
1.一种基于无网格法的线源时域电磁响应数值计算方法，其特征在于：
包括如下的步骤：
1)在计算区域内进行节点、电性参数、支持域、背景网格、初始场值设置；
2)在计算区域内对背景网格进行等参数变换，在每个背景单元中分别计算
高斯积分点及其权系数；
3)对所有背景网格进行外循环，对所有高斯积分点内循环，搜索背景网格
中局部定义域内的有效节点并计算局部定义域内节点处形函数；
4)加载强加边界条件得到线性方程组，加载吸收边界条件得到边界阻尼矩
阵，将边界阻尼矩阵加入到线性方程组中；
5)采用LU分解法求解线性方程组，得到各个节点各个时刻的场值。
2.按照权利要求1所述的一种基于无网格法的线源时域电磁响应数值计算
方法，其特征在于：
步骤4中，从线源瞬变电磁满足的控制方程和边界条件出发，建立和时域
电磁边值问题的等价泛函数，利用罚因子法处理计算模型的强加边界条件，求
取泛函数的变分，得到刚度矩阵K，阻尼矩阵K'和右端项S，E为待求的场值，
形成线性方程组：
K′E·+KE=S]]>3.按照权利要求2所述的一种基于无网格法的线源时域电磁响应数值计算
方法，其特征在于：
步...

【专利技术属性】
技术研发人员：嵇艳鞠，黄廷哲，黄婉玉，关珊珊，吴琼，林君，
申请(专利权)人：吉林大学，
类型：发明
国别省市：吉林;22