CN107766666B - 一种基于分数阶差分法的三维时域电磁反常扩散模拟方法 - Google Patents

Abstract

本发明涉及一种基于分数阶差分法的三维时域电磁反常扩散模拟方法，根据复杂岩石结构的非均质性、不连续性、孔隙性等特征，将岩石的粗糙度与极化率参数引入，建立分数阶电导率模型，推导频域电磁场的分数阶扩散方程；通过频时转换获得电磁场的时域分数阶微分‑积分表达式；采用Riemann‑Liouville分数阶积分和有限差分方法，对扩散方程的积分和微分项进行时域离散，构建电、磁场的时域迭代公式；加载初始条件和边界条件，实现了三维时域电磁反常扩散的数值模拟。本发明有益效果在于，提出了一种复杂岩石结构的分数阶电导率模型，可以准确描述地下随机介质的慢扩散和多孔极化介质的极化反常扩散现象，为研究复杂地质构造的电磁波传播机理提供理论依据。

Description

一种基于分数阶差分法的三维时域电磁反常扩散模拟方法

技术领域

本发明涉及一种基于分数阶差分法的三维时域电磁反常扩散模拟方法，适用于时域电磁反常扩散的三维仿真，尤其是对电磁慢扩散与极化效应的高精度三维数值模拟。

背景技术

时间域瞬变电磁法(Time domain Transient electromagnetic methods)利用长导线源或回线源向地下输出时变电流，激励大地介质产生感应电磁场，通过测量电场或磁场信号，探测地下介质的电性差异和结构。该方法具有分辨低阻能力强、抗电磁干扰性好等优势，在工程地质勘查、资源探测、构造带勘探等领域发挥着重要作用。随着高精度传感测量技术的快速进步，近些年来，学者们在复杂地质条件下进行时域电磁测量时，观察到了一些与经典电磁场扩散不相符合的现象，并把这些不遵守经典Fick扩散定律或布朗运动规律的扩散行为统称为电磁反常扩散(Electromagnetic Anomalous Diffusion)。时域电磁勘探中常见的反常扩散有两种，一种是极化(InducedPolarization)效应，另一种则是在测量数据晚期出现的区别于经典幂律衰减的慢扩散现象(Electromagnetic Subdiffusion)。反常扩散出现的原因是通常我们假设电导率是各向同性均匀分布的，而实际地质构造是非均匀的，含有孔隙和裂缝的复杂介质。因此只有通过重新定义电导率模型，才能有效的模拟这些电磁反常扩散现象，为后期的数据解释提供有效的依据。

目前国内外关于电磁反常扩散的研究主要有两种：一种是极化效应的研究，以三维Cole-Cole模型的复杂极化体电磁响应数值计算为主，Marchant等人(2015)在频率域进行欧姆定律离散，再采用数字滤波方法进行频时变换，最后基于整数阶传统有限差分方法计算了多孔极化介质的感应场和极化场响应；另一种则是随机介质慢扩散研究，MarkE.Everett团队(2009，2012，2015)，基于可控源电感应磁法，推导了随机介质模型的频率域分数阶扩散方程，并基于G-S变换算法，计算了时域随机介质模型的电磁响应。而国内关于随机介质的研究只适用于电导率随空间随机变化的结构，不能计算随时间或频率变化的随机地质构造。

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

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

以上所述方法公布了国内外关于反常扩散的研究方法。但在电磁反常扩散方面，目前还没有关于在时间域直接进行分数阶三维有限差分运算的研究。且国内外的研究均是将极化效应与慢扩散效应的研究分开进行的，实际上这两种现象均属于电磁反常扩散现象，可以用同一个电导率模型来有效解释。因此，提出一种能够同时模拟极化效应与慢扩散现象的电导率模型，并且能在时间域直接实现分数阶有限差分方法的三维数值模拟，对后期测量数据的有效解释具有重大的意义。

发明内容

本发明所要解决的技术问题在于提供一种基于分数阶差分法的三维时域电磁反常扩散模拟方法，根据复杂岩石结构的非均质性、不连续性、孔隙性等特征，将岩石的粗糙度与极化率参数引入，建立分数阶电导率模型，推导频域电磁场的分数阶扩散方程；通过频时转换获得电磁场的时域分数阶微分-积分表达式；采用Riemann-Liouville分数阶积分和有限差分方法，对扩散方程的积分和微分项进行时域离散，构建电、磁场的时域迭代公式；加载初始条件和边界条件，实现了三维时域电磁反常扩散的数值模拟。本发明有益效果在于，提出了一种复杂岩石结构的分数阶电导率模型，可以准确描述地下随机介质的慢扩散和多孔极化介质的极化反常扩散现象，且通过对分数阶积分项的离散近似，实现了三维的时域电磁反常扩散高精度数值模拟，为研究复杂地质构造的电磁波传播机理提供理论依据。本发明是这样实现的，一种基于分数阶差分法的三维时域电磁反常扩散模拟方法：

1)、采用准静态条件下，无源Maxwell旋度方程作为电、磁场分量E_x、E_y、E_z、H_x、H_y的控制方程，散度方程

作为磁场垂直分量H_z的控制方程；

2)、在频率域电导率表达式中引入粗糙度α与极化率β(0<α<1,-1<β<1)；

3)、将频率域电导率表达式代入频域控制方程中，进行频-时转，产生分数阶积分项；

4)、采用Riemann-Liouville分数阶积分方程对分数阶积分项进行离散近似；

5)、基于有限差分算法对整个控制方程进行离散，推导出各电磁场分量的递归方程；

6)、代入初边条件，完成所有时刻的迭代计算，保存数据，绘制成图，进行结果分析；

进一步地，步骤2中，频率域电导率表达式中引入粗糙度α与极化率β(0<α<1,-1<β<1)后，其表达形式为：

σ(ω)＝σ₀-βσ₀(iω)^-a (1)

式(1)中σ(ω)为广义电导率表达式，σ₀为经典的电导率取值，α为粗糙度，β为极化率，粗糙度α表征介质的非均匀程度，极化率β表征介质孔隙内电导率与介质电导率的差异程度。当β大于零时表示孔隙内物质导电性优于介质本身，此时会产生极化现象；当β小于零时表示孔隙内物质导电性劣于介质本身，此时会产生电磁慢扩散现象；当β等于零时表示孔隙内物质导电性能与介质本身相同，那么此时介质即为均匀介质，与经典电磁扩散现象相符。

进一步地，步骤3中，将电导率表达式(1)代入频域电场的控制方程

将公式(2)进行时频转换，(2)中复频变量的负分数次幂(iω)^-α在时间域的表达式中为分数阶积分项：

进一步地，步骤4中，采用Riemann-Liouville分数阶积分方程对公式3中的分数阶积分项进行离散近似，时刻取t_n+1：

公式(4)中将(0～t_n+1)时间范围内的积分按时间步长的选取离散成(t₀～t₁)、(t₁～t₂)、……、(t_n～t_n+1)内积分的求和，E(i)表示t_i时刻的电场值，且在(t_i～t_i+1)范围内电场值用

近似，Г(α)表示gamma函数。

进一步地计算公式(4)中的积分可得

进一步地，步骤5中，基于有限差分方法对电场分量控制方程离散，推导递归方程时将公式(5)代入到公式(3)，在t_n+1/2时刻进行离散，以x方向电场为例：

进一步地，对公式(6)进行整理，得出E_x的递归关系式：

其中：

本发明与现有技术相比，有益效果在于：通过在电导率表达式中引入粗糙度α与极化率β,在时间域直接对分数阶积分进行离散近似，实现了在时间域基于分数阶有限差分算法的时域电磁反常扩散现象三维仿真。且两个参数的引用能够有效地解释电磁反常扩散中的慢扩散和极化现象。本方法为正确解释实测数据中的反常扩散现象提供了理论依据，对电磁法在矿产资源勘探等领域的应用和发展具有重大意义。

附图说明

图1是基于分数阶积分-差分法的三维时域电磁反常扩散模拟方法流程图；

图2是粗糙度与极化率对接收感应电动势的影响；

图3是三维时域电磁反常扩散中，粗糙度为0.3极化率为-1时，在10毫秒处地面感应电动势平面图；

图4基于分数阶积分-差分法的三维时域电磁反常扩散接收线圈中心感应电动势对比；

具体实施方式

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

实施例

参见图1，一种基于分数阶差分法的三维时域电磁反常扩散模拟方法包括：

1)、选取准静态条件下、无源Maxwell方程作为电磁场分量的控制方程，散度方程

作为约束条件。采用非均匀Yee式网格进行刨分，网格数目为107×107×53，最小和最大网格步长分别为10m和120m。选取DuFort-Frankel方法计算各个时间步长。

2)在整个计算区域内设置电导率、磁导率、粗糙度、极化率，输入发射装置信息，算例中模型为均匀半空间模型，电导率设置为10西门子/米，磁导率设置为真空磁导率。

3)、计算初始时刻电磁响应值，设置边界条件，代入电磁场迭代方程中。

4)、计算分数阶积分项的值，代入电场迭代式中计算当前时刻电场值，计算下半个时刻磁场值。

5)、完成全部时刻的电磁场数值计算，保存计算结果，成像分析。

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

Claims (2)

1.一种基于分数阶差分法的三维时域电磁反常扩散模拟方法，特征在于，包括如下步骤：

1)、采用准静态条件下，无源Maxwell旋度方程作为电、磁场分量E_x、E_y、E_z、H_x、H_y的控制方程，散度方程

作为磁场垂直分量H_z的控制方程；

2)、在频率域电导率表达式中引入粗糙度α与极化率β，0＜α＜1，-1＜β＜1；

3)、将频率域电导率表达式代入频域控制方程中，进行频率域到时间域的转换，产生分数阶积分项；

4)、采用Riemann-Liouville分数阶积分方程对分数阶积分项进行离散近似；

5)、基于有限差分算法对整个控制方程进行离散，推导出各电磁场分量的递归方程；

6)、代入初始条件及边界条件，完成所有时刻的迭代计算，保存数据，绘制成图，进行结果分析。

2.按照权利要求1所述的一种基于分数阶差分法的三维时域电磁反常扩散模拟方法，其特征在于：

其中步骤2中，频率域电导率表达式中引入粗糙度α与极化率β，0＜α＜1，-1＜β＜1，后，其表达形式为：

σ(ω)＝σ₀-βσ₀(iω)^-α (1)

式(1)中σ(ω)为广义电导率表达式，σ₀为经典的电导率取值，β为极化率，粗糙度α表征介质的非均匀程度，极化率β表征介质孔隙内电导率与介质电导率的差异程度；当β大于零时表示孔隙内物质导电性优于介质本身，此时会产生极化现象；当β小于零时表示孔隙内物质导电性劣于介质本身，此时会产生电磁慢扩散现象；当β等于零时表示孔隙内物质导电性能与介质本身相同，那么此时介质即为均匀介质，与经典电磁扩散现象相符；

其中步骤3中，将电导率表达式(1)代入频域电场的控制方程

将公式(2)进行时频转换，(2)中复频变量的负分数次幂(iω)^-α在时间域的表达式中为分数阶积分项：

其中步骤4中，采用Riemann-Liouville分数阶积分方程对公式(3)中的分数阶积分项进行离散近似，时刻取t_n+1：

公式(4)中将(0～t_n+1)时间范围内的积分按时间步长的选取离散成(t₀～t₁)、(t₁～t₂)、……、(t_n～t_n+1)内积分的求和，E(t_i)表示t_i时刻的电场值，且在(t_i～t_i+1)范围内电场值用

近似，Γ(α)表示gamma函数；

进一步地计算公式(4)中的积分可得

其中步骤5中基于有限差分方法对电场分量控制方程离散，推导递归方程时将公式(5)代入到公式(3)，在t_n+1/2时刻进行离散，以x方向电场为例：

进一步地，对公式(6)进行整理，得出E_x的递归关系式：

其中：

Images