CN108897052B - 一种基于分数阶线性近似的三维时域电磁慢扩散模拟方法 - Google Patents

Abstract

本发明涉及一种基于分数阶线性近似的三维时域电磁慢扩散模拟方法，通过将广义电导率表达式引入Maxwell方程组产生分数阶积分项，将含有分数阶积分项的等式两端同时求导，将分数阶积分转为分数阶微分。将分数阶微分项进行线性近似，并进行频时转换得到近似的时间域的整数阶微分，完成从分数阶到整数阶的转换，避免了分数阶微积分运算对每一时刻电场值的存储。最后基于有限差分算法对控制方程各偏导项进行近似，并推导出电场和磁场各分量迭代方程。最终实现了三维时域电磁慢扩散的高精度、长时窗数值模拟。本发明目的在于可以克服分数阶微积分时间域直接运算具有弱奇异性及内存消耗巨大的问题，实现了分数阶电导率模型的快速模拟计算。

Description

一种基于分数阶线性近似的三维时域电磁慢扩散模拟方法

技术领域

本发明涉及一种地球物理勘探领域的三维时域电磁慢扩散模拟方法，尤其是对时域航空电磁慢扩散的数值模拟。

背景技术

随着地球物理勘探技术的快速发展，高精度的勘探仪器的研制，尤其是超导传感器(SQUID，Superconducting Quantum Interference Device)的逐渐成熟和应用，很大程度上提升了探测分辨率和深部勘探能力。随着勘探精度的提高，研究人员在进行时域电磁测量时，观察到了一些与电磁场经典扩散不相符的现象，其中一种被称为慢扩散现象(Sub-diffusion)。慢扩散主要是指感应电动势在地层结构均匀地区衰减速度慢于理论的5/2幂率衰减的现象，目前已经引起了地质工作者的高度重视。

时域航空瞬变电磁法(Time domain Airborne Transient electromagneticmethods)发射系统与接收系统均在空中，相比于地面瞬变电磁法(Ground Transientelectromagnetic methods)，航空电磁法能够更加适应复杂多变的地质环境，目前已被大量的用于金属矿、煤炭、油气的勘查，水文地质调查和环境监测等领域。

时域有限差分(Finite-difference Time-domain)方法是电磁场数值计算的重要方法之一，它的基本思想是用中心差商代替场量对时间和空间的一阶偏微商,通过在时域的递推模拟波的传播过程,从而得出场分布，目前已被广泛应用于电磁场数值模拟计算中。

在国外关于电磁场慢扩散研究方面，Mark E.Everett团队(2009，2012，2015)，基于可控源电磁感应法，推导了粗糙介质模型的频率域分数阶扩散方程，并基于G-S变换算法，计算了时域随机介质模型的电磁响应。而国内关于电磁慢扩散研究主要围绕电导率随空间随机变化的结构，并没有关于电导率算随时间变化的慢扩散研究。

中国专利CN107766666A公开了一种基于分数阶差分法的三维时域电磁反常扩散模拟方法,采用Riemann-Liouville分数阶积分和有限差分方法对扩散方程的积分和微分项进行时域离散，实现了三维时域电磁反常扩散的数值模拟。

中国专利CN 104392127A公开了一种基于离散分数阶差分的反常扩散模拟方法，采用描述扩散现象的一维经典扩散方程定义离散分数阶差分，利用离散分数阶差分将经典扩散方程离散化，根据初边界条件进行数值模拟。

中国专利CN 106776478A公开了一种反常扩散中的基于分步计算的离散分数阶差分方法，减小了gamma函数计算的限制，扩展模拟的点数，提高模拟的效率。

中国专利CN107657137A公开了一种有理函数逼近的分数阶电磁反常扩散三维模拟方法，基于频域有理函数逼近法，构建科尔-科尔模型分数阶传递函数和n阶有理逼近函数，采用部分分式展开法和拉普拉斯逆变换获得电导率的时域形式，实现分数阶科尔-科尔模型三维电磁响应数值计算。

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

发明内容

本发明所要解决的技术问题在于提供一种基于分数阶线性近似的三维时域电磁慢扩散模拟方法，通过将广义电导率表达式引入Maxwell方程组产生分数阶积分项，将含有分数阶积分项的等式两端同时求导，将分数阶积分转为分数阶微分。将分数阶微分项进行线性近似，并进行频时转换得到近似的时间域的整数阶微分，完成从分数阶到整数阶的转换，避免了分数阶微积分运算对每一时刻电场值的存储，大大降低了计算机内存的消耗。最后基于有限差分算法对控制方程各偏导项进行近似，并推导出电场和磁场各分量迭代方程。最终实现了三维时域电磁慢扩散的高精度、长时窗数值模拟。本发明目的在于可以克服分数阶微积分时间域直接运算具有弱奇异性及内存消耗巨大的问题，能够更高效准确地的模拟三维时域电磁慢扩散现象。

本发明是这样实现的，一种基于分数阶线性近似的三维时域电磁慢扩散模拟方法：

1)、将广义电导率表达式引入Maxwell方程组作为电磁场迭代控制方程，表征为含有复频变量负分数次幂项的电磁场迭代控制方程；

2)、将控制方程中的复频变量负分数次幂转化为正分数次幂，进行线性近似，得到含有整数次幂的控制方程；

3)、将控制方程转换由频率域转换到时间域，再采用有限差分方法对各项进行差分近似；

4)、推导电磁场各分量迭代格式；

5)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置计算域电导率、磁导率和人工介电常数，加载C-PML边界条件；

6)、在观测时间内开展电、磁场各分量迭代运算；

7)、迭代计算结束后，提取电磁场各分量响应，并对计算结果进行显示。

进一步地，步骤1)中，广义电导率表达式为：

σ(ω)＝σ'+(iω)^-βσ” (1)

其中β取值在0与1之间；

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

▽×H＝σ'E+σ”(iw)^-βE+ε(iw)E

(2)

▽×E＝-μ(iω)H (3)

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

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

▽×(iw)H＝σ'(iw)E+σ”(iw)^αE+ε(iw)²E (4)

其中α＝1-β；

采用线性近似方法将公式(4)中的复频变量分数阶次幂(iw)^α进行整数阶近似，由于α取值在0与1之间，因此(iw)^α可以近似表示为：

(iω)^α≈α(iω)¹+(1-α)(iω)⁰ (5)

将(5)代入(2)中整理可得：

▽×(iw)H＝(σ'+σ”α)(iw)E+(1-α)σ”E+ε(iw)²E (6)

进一步地，步骤3)中将控制方程转化到时间域可得：

由于瞬变电磁法主要检测感应电动势，因此研究人员正演时仅需提取分量，因此在迭代时可以将作为一个整体参与迭代，令则(7)、(8)可以表示为：

▽×E＝-μH_m (10)

采用有限差分方法对(9)中偏导项在n时刻进行差分近似可得：

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

进一步地，步骤4)中将(11)(12)代入公式(9)(10)并进行旋度展开，得到电磁场各分量迭代格式为，以x方向电场磁场为例：

本发明与现有技术相比，有益效果在于：在时域电磁三维数值模拟控制方程中引入广义电导率定义，能够有效模拟电磁慢扩散现象；针对分数阶微积分在时间域直接计算内存消耗巨大、弱奇异性、消耗时间长等问题，采用频域线性近似、时频转换的方法，将分数阶微积分近似为整数阶微分运算，达到了对三维时域电磁慢扩散现象的高效、精准数值模拟的目的。

附图说明

图1是一种基于分数阶微分的三维时域电磁慢扩散模拟方法示意图；

图2是基于分数阶微分的三维时域电磁慢扩散模拟得到的接收线圈感应电动势衰减曲线与解析解对比图；

图3是基于分数阶微分的三维时域电磁慢扩散模拟得到的接收线圈感应电动势衰减曲线与经典模型感应电动势衰减曲线对比图；

图4是是基于分数阶微分的三维时域电磁慢扩散模拟得到空中二维电磁响应平面图；

具体实施方式

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

实施例

参见图1，一种基于分数阶线性近似的三维时域电磁慢扩散模拟方法，包括：

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

2)在整个计算区域内设置电导率、磁导率、C-PML系数等参数，算例中模型为均匀半空间模型，电导率设置为σ’＝0.008西门子/米，σ”＝0.002西门子/米，α＝0.1，磁导率设置为真空磁导率。

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

E_z＝0 (17)

其中I为发射电流，r为发射线圈半径，J₀为零阶贝塞尔函数，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释放出来，并对计算结果进行显示，完成全部计算。

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

Claims (1)

1.一种基于分数阶线性近似的三维时域电磁慢扩散模拟方法，其特征在于，包括如下步骤：

1)、将广义电导率表达式引入Maxwell方程组作为电磁场迭代控制方程，表征为含有复频变量负分数次幂项的电磁场迭代控制方程；

2)、将控制方程中的复频变量负分数次幂转化为正分数次幂，进行线性近似，得到含有整数次幂的控制方程；

3)、将控制方程由频率域转换到时间域，再采用有限差分方法对各项进行差分近似；

4)、推导电磁场各分量迭代格式；

5)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置计算域电导率、磁导率和人工介电常数，加载C-PML边界条件；

6)、在观测时间内开展电、磁场各分量迭代运算；

7)、迭代计算结束后，提取电磁场各分量响应，并对计算结果进行显示；

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

σ(ω)＝σ'+(iω)^-βσ” (1)

其中β取值在0与1之间；

其中步骤2)中将控制方程中的复频变量负分数次幂转化为正分数次幂，再进行线性近似后得到的整数阶频域控制方程为:

其中ε为介电常数，μ为导磁率，E为电场，H为磁场，α＝1-β；

其中步骤3)中将控制方程转化到时间域可得：

由于瞬变电磁法主要检测感应电动势，因此研究人员正演时仅需提取分量，因此在迭代时可以将作为一个整体参与迭代，令则(4)、(5)可以表示为：

采用有限差分方法对(6)中偏导项在n时刻进行差分近似可得：

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

其中步骤4)中将电磁场控制方程进行展开，差分近似，得到的电磁场各分量迭代格式为，以x方向电磁场为例：