阅读说明：本技术 一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法 (Three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation ) 是由 嵇艳鞠 王世鹏 赵雪娇 吴琼 黎东升 关珊珊 于 2019-11-12 设计创作，主要内容包括：本发明涉及一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法,通过将复电导率模型引入频域Maxwell方程组后,电磁场扩散方程中含有复频变量的负分数次幂项,先进行频-时变换得到含有Caputo分数阶微分项的时间域控制方程；再采用Alikhanov超收敛插值逼近方法,对电场控制方程中Caputo分数阶导数进行超收敛逼近,获得分数阶微分项的非均匀步长离散近似表达式,从而完成时间域分数阶微分项的稳定、高精度直接求解；最后基于有限差分算法对控制方程进行离散,推导出电场和磁场迭代方程,最终实现了三维时域电磁慢扩散的高精度数值模拟。本发明目的在于克服分数阶微分求解的弱奇异不稳定性及误差较大问题,实现三维时域电磁慢扩散的高精度数值模拟。(The invention relates to a three-dimensional electromagnetic slow diffusion numerical simulation method based on hyperconvergence interpolation approximation.A complex conductivity model is introduced into a frequency domain Maxwell equation set, an electromagnetic field diffusion equation contains a negative fraction power term of a complex frequency variable, and frequency-time transformation is carried out to obtain a time domain control equation containing a Caputo fractional order differential term; then, performing hyperconvergence approximation on the Caputo fractional derivative in the electric field control equation by adopting an Alikhanov hyperconvergence interpolation approximation method to obtain a non-uniform step discrete approximate expression of a fractional derivative term, thereby completing the stable and high-precision direct solution of the time domain fractional derivative term; and finally, dispersing the control equation based on a finite difference algorithm, deducing an electric field and magnetic field iterative equation, and finally realizing high-precision numerical simulation of three-dimensional time domain electromagnetic slow diffusion. The invention aims to solve the problems of weak singular instability and large error of fractional order differential solution and realize high-precision numerical simulation of three-dimensional time domain electromagnetic slow diffusion.)

一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法

技术领域

本发明涉及一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法，适用于时域电磁慢扩散数值模拟，尤其是对磁源电磁慢扩散的数值模拟。

背景技术

磁源瞬变电磁探测方法，选用线圈作为发射源，加载阶跃电流激发地下良导体产生涡流，通过接收机采集二次感应场获得地下丰富的电性信息，由于其具有探测深度大、经济、简便的优点，已经被广泛应用于地质资源探测中。超导量子干涉仪(SQUID，Superconducting Quantum Interference Device)和原子磁力仪的逐渐成熟与实际应用，显著提高了勘探的精度，获得了实际地质中更加准确的电磁衰减信号，不符合理论-5/2幂律的瞬变电磁慢扩散现象逐渐被观测，随着对地质结构特征的不断研究，传统的单一尺度、单一参数的地质电导率模型已经不能满足高精度地质探测的要求。建立更符合地质构造的物理模型是提高电磁探测分辨率的关键技术之一。

三维时域有限差分(Three-dimensional Finite-difference Time-domain)方法是电磁场数值计算的重要方法之一，它的基本思想是场量对时间和空间的一阶偏导数用中心差分近似,通过在时域的递推模拟波的传播过程,从而得出场分布，与一维和二维差分方法相比，三维差分更适用于观察域广、结构复杂的地质探测分析中，目前已被广泛应用于电磁场数值模拟计算中。

Chester J.Weiss等在美国布拉索斯郡利用TEM47观测到时域电磁慢扩散响应，并排除了慢扩散是由分层结构或横向异性导致的可能性，得出慢扩散响应是由粗糙介质引起的一种分数阶扩散。Mark E.Everett等人提出介质具有非均匀性，即粗糙度，并将粗糙参数引入到频率域电导率的表达式中。Detwiler等研究了裂隙中溶质的迁移问题，发现裂隙介质中溶质迁移是裂隙内部的泰勒扩散和裂隙表面的常次扩散的双尺度耦合扩散；在地质问题中，大尺度粗糙介质裂隙介质中电磁波的传播也具有相同的双耦合性质，即常次扩散和次扩散的双尺度耦合扩散。国内主要对电导率随空间随机变化的结构进行研究，很少对电导率随时间变化的慢扩散开展研究。

中国专利CN107766666A公开了一种基于分数阶差分法的三维时域电磁反常扩散模拟方法,通过频时转换获得电磁场的时域分数阶微分-积分表达式；采用Riemann-Liouville分数阶积分和有限差分方法,对扩散方程的积分和微分项进行时域离散,构建电、磁场的时域迭代公式；加载初始条件和边界条件,实现了三维时域电磁反常扩散的数值模拟。实现了三维时域电磁反常扩散的数值模拟。

中国专利CN 106776478A公开了一种反常扩散中的基于分步计算的离散分数阶差分方法，将时间空间反常扩散方程的离散格式中的两个gamma函数的比值看做一个参数,利用gamma函数的递推关系,结合数值模拟中的分布计算的思路和整体考虑的思想减小了gamma函数计算的限制，扩展模拟的点数，提高模拟的效率。

中国专利CN108897052A公开了一种基于分数阶线性近似的三维时域电磁慢扩散模拟方法，将广义电导率引入Maxwell方程组,将分数阶微分项进行线性近似，并进行频时转换得到近似的时间域的整数阶微分，完成从分数阶到整数阶的转换，达到模拟三维时域电磁慢扩散现象的目的。

以上所述方法公布了国内外关于反常扩散的研究方法。但对电磁反常扩散中慢扩散现象，目前几乎没有关于在时间域进行分数阶三维有限差分运算的研究，且对于分数阶算子在时间域直接运算具有弱奇异性导致结果不稳定及误差大等问题，如何在电磁慢扩散领域中准确、稳定的进行分数阶差分运算，是本领域技术人员迫切解决的一个技术问题。

发明内容

本发明所要解决的技术问题在于提供一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法，通过将复电导率表达式引入频域Maxwell方程组产生含有复频变量的负分数次幂项，进行频-时变换获得含有Caputo分数阶微分项的时间域控制方程。采用Alikhanov超收敛插值逼近方法，对电场控制方程中Caputo分数阶导数进行超收敛逼近，得出分数阶微分项的非均匀步长离散近似公式，完成时间域分数阶微分项的稳定、高精度直接求解；最后基于有限差分算法对控制方程空间和时间偏导项进行离散，并推导出电场和磁场各分量迭代方程。最终实现了三维时域电磁慢扩散的稳定、高精度数值模拟。本发明目的在于克服分数阶微分求解的弱奇异不稳定性及误差较大问题，实现三维时域电磁慢扩散的高精度数值模拟。

本发明是这样实现的，一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法：

1)、将复电导率表达式引入Maxwell方程组作为电磁场迭代控制方程，引入复电导率后的电磁扩散过程早期与经典扩散一致，晚期衰减缓慢，符合实际观测中的电磁拖尾现象；

2)、将控制方程中的复频变量负分数次幂转化为正分数次幂，并进行频-时转换得到含有Caputo分数阶微分项的时间域电磁场控制方程；

3)、对电场控制方程中的电场分量进行二次插值多项式近似，再对等式两端进行求导，代入Caputo分数阶导数超收敛插值逼近公式，得出分数阶微分项的非均匀步长离散近似公式，完成分数阶微分方程的直接求解；

4)、采用有限差分方法对控制方程各偏导项进行离散，递归得出电磁场各分量的控制方程；

5)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置计算域电导率、磁导率和人工介电常数，计算初始场；

6)、加载C-PML边界条件并使用GPU加速，并在观测时间内开展电、磁场各分量迭代运算；

7)、迭代计算结束后，提取电磁场各分量响应并进行成图，对结果进行分析处理；

其中步骤1)中，复电导率表达式为：

σ(ω)＝σ₀+ησ₀(iω)^-β (1)

式(1)中σ₀为直流电导率，η为权重系数，i为虚数单位，ω为角频率，β为粗糙度，β和η的取值均大于0小于1；这种电导率定义可以模拟电磁场慢扩散现象，当权重系数等于零时，此电磁扩散过程退化为经典的电磁扩散过程；

引入复电导率表达式后电磁场的控制方程可以表示为：

其中ε为真空中介电常数，μ为磁导率，E为电场，H为磁场，

为哈密顿算子；

进一步地，步骤2)中将控制方程中的复频变量负分数次幂转化为正分数次幂采用的方法是将公式(2)两端同时进行微分处理，即同时乘iω，可以得到：

其中步骤2)中对频率域控制方程进行频-时变换得到的时间域控制方程为:

其中H_m为磁场，

为Caputo分数阶导数，α为一非负实数，0为积分下限，t为积分上限，Caputo分数阶导数表达式为：

式(7)被称为左侧α阶Caputo分数阶导数，a为积分下限，x为积分上限，f(x)为一函数，n为一整数，且n-1≤α＜n，Γ为Gamma函数，其表达式为：

其中步骤3)中对电场控制方程中电场分量进行二次插值多项式近似及求导过程为：

在区间[t_k-1,t_k]上，利用三点(t_k-1,E(t_k-1)),(t_k,E(t_k)),(t_k+1,E(t_k+1))对E(t)做二次插值多项式近似可以得到：

式(9)的导数为：

其中E(t_k-1)、E(t_k)、E(t_k+1)为电场在不同时刻的电场值，Δt_k＝t_k-t_k-1，Δt_k+1＝t_k+1-t_k；

其中步骤3)中对控制方程中分数阶微分项进行超收敛近似直接求解的过程为：

Caputo分数阶超收敛格式为：

式(11)中

将式(10)代入式(11)得到任意步长下分数阶微分的超收敛离散近似格式：

其中步骤4)中采用有限差分方法对式(6)中偏导项进行离散可得：

其中Δt_n+1＝t_n+1-t_n，Δt_n＝t_n-t_n-1；

将电磁场控制方程进行展开，进行差分离散近似，得到电磁场各分量迭代方程，其中x方向电磁场迭代方程为：

本发明与现有技术相比，有益效果在于：在时域电磁三维数值模拟控制方程中引入复电导率模型，能够有效模拟电磁慢扩散现象；针对分数阶微分项在时间域直接求解具有弱奇异不稳定性、近似求解整体误差较大的问题，采用超收敛插值逼近方法对分数阶微分项进行直接求解，解决了求解过程中的震荡问题，达到了对三维时域电磁慢扩散现象的高精度数值模拟的目的。

附图说明

图1是一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法示意图；

图2是基于分数阶微分项超收敛直接求解的三维时域电磁慢扩散模拟得到的接收线圈感应电动势衰减曲线与数值积分解对比图及误差曲线；

图3是基于复电导率模型的三维时域电磁慢扩散模拟得到的接收线圈感应电动势衰减曲线与经典模型、广义电导率模型感应电动势衰减曲线对比图；

图4是基于复电导率模型的三维时域电磁慢扩散数值模拟得到地面二维电磁响应平面图；

图5是基于复电导率模型的三维时域电磁慢扩散数值模拟得到空中二维电磁响应平面图；

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例

参见图1，一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法，包括：

1)、推导电场和磁场控制方程，结合非均匀Yee式网格和DuFort-Frankel方法构建电磁场各分量显式差分格式，网格数目为117×117×58，其中x、y方向上的网格数均为117个，z方向上网格数58个，最小网格步长为10m，最大网格步长为120m，计算区域中除地面外其余5个面的最外层8个网格加载CFS-PML边界条件。

2)、在整个计算区域内设置电导率、磁导率、CFS-PML系数等参数，算例中模型为均匀半空间模型，电导率设置为σ₀＝0.05S/m，粗糙参数β＝0.6、权重参数η＝0.1，磁导率设置为真空磁导率。

3)、将计算区域的电性参数代入到发射部分在空中，接收部分在地面的电场响应表达式中，计算t₀、t₁时刻的x，y方向的电场响应和t₁时刻的x，y，z方向的磁场响应。

E_z＝0 (19)

其中I为发射电流，r为发射线圈半径，J₁为1阶贝塞尔函数，h为发射线圈距地面高度，iω为拉普拉斯变量，λ为汉克尔变换积分变量，式中

4)、将参与迭代的CPU序列转化为GPU序列，加速运算。

5)、计算当前时刻迭代时间步长。

Δ_min为最小网格步长，a取值为0.1。

6)、将电场和磁场初始场值代入电场迭代方程中在整个计算区域更新电场值E_x、E_y、E_z，并保存前一时刻与当前时刻电场值。

7)、将当前时刻电场值代入迭代方程中在整个计算区域更新磁场值H_x、H_y、H_z，并保存当前时刻磁场值。

8)、使用向上延拓方法将地面电磁场延拓到空中，求得空中接收线圈处电磁场值。

9)、判断是否完成全部迭代，若没有完成，重复5-8步骤，如果完成全部迭代，将计算结果从GPU释放出来，并对计算结果进行显示、成图，完成全部计算。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。