CN105893678A - 一种时域有限差分的三维感应-极化双场数值模拟方法 - Google Patents

Abstract

本发明涉及一种时域有限差分的三维感应‑极化双场数值模拟方法，目的在于快速计算三维模型的感应‑极化双场的电磁响应。主要包括采用逆拉普拉斯变换获得迪拜模型电导率的时域表达式，构建电导率参数的e指数辅助方程，通过梯形积分法获得欧姆定律时域离散递推表达式，由四维数值运算降低为三维运算，再将其代入无源Maxwell旋度方程中，基于三维时域有限差分方法推导电场和磁场的迭代方程，进而完成三维模型的感应‑极化双场电磁响应数值计算。本发明目的在于解决欧姆定律的时域卷积运算耗时长、内存占用大等问题，最终实现三维模型的感应‑极化双场电磁响应快速数值计算。

Description

一种时域有限差分的三维感应-极化双场数值模拟方法

技术领域

本发明涉及一种地球物理勘探领域中的时域电磁场数值模拟方法，尤其适用快速计算三维模型的感应-极化双场电磁响应。

背景技术

瞬变电磁法(Time-domain Electromagnetic，简称TEM)是基于电磁感应涡流效应为主的地下近地表探测方法。当低频时忽略位移电流和极化电荷产生的电磁场，仅利用地下介质的涡流效应观测二次瞬变场。时域激发极化法(Time-domain InducedPolarization，简称IP)是基于地下介质极化效应为主的地下浅层探测方法，主要采用电性源施加直流或交流电磁场激励下，忽略传导电流产生的电磁场以及高频介电常数引起的极化响应，仅利用时域电流断开后低频激励下地下介质产生的极化场，观测极化电荷产生的二次电位变化过程。在时域激发极化测量中总是去除感应场的耦合。

但是，无论是瞬变电磁还是时域激发极化方法，都是以电磁场麦克斯韦方程为理论依据并进行波场简化近似，实际探测中感应场和极化场同时存在，互相伴生。无论瞬变电磁还是时域激发极化方法，将两种电流场割裂开分别进行研究，是不符合实际电磁场的传播规律的，可能导致实测数据与理论计算结果不符，甚至导致解释结果出现错误。为此，只有准确计算地下结构的感应和极化双场电磁响应，理论模拟结果与实际测量数据互相吻合。

中国专利CN104408021A公开了一种电偶源三维时域有限差分正演成像方法，对海洋空气、海水和海底大地三分空间均建立麦克斯韦方程组和本构方程，采用时域有限差分法得到海水和海底大地中任意时刻电磁场的分布。

中国专利CN105277980A公开了一种高精度空间和时间任意倍数可变网格 有限差分正演方法，通过建立地下介质的正演速度模型，对正演速度模型中的声波波场进行二维网格离散化，同时对最佳匹配层边界条件进行网格离散化，进而通过声波波动方程进行时域有限差分正演模拟。

加拿大专利CA2388271公开了一种基于时域有限差分方法的电磁场计算方法，并提出了一种基于二维导体的网格划分方法和数据储存方式，但只涉及感应场的数值计算。

加拿大，D Marchant(2015，University of British Columbia)采用有限元方法，在时间域直接求解并模拟计算了三维空间下存在激发极化效应的电磁响应。

以上所述方法公布了基于时域有限差分(FDTD)方法的电磁响应计算方法，国内外专利还未涉及时域有限差分方法计算三维模型的感应和极化双场的电磁响应方面的，为此，本发明基于FDTD方法直接计算三维模型的感应和极化双场电磁响应，并解决欧姆定律的时域卷积计算耗时长、内存占用大等问题，最终实现三维模型的感应-极化双场电磁响应快速数值计算。

发明内容

本发明所要解决的关键问题是提供一种时域有限差分的三维感应-极化双场数值模拟方法，先获得迪拜模型的电导率时域形式，通过构建电导率参数e指数辅助方程，并采用梯形积分法获得欧姆定律的时域离散递推形式，最后基于FDTD方法，实现了时域内快速计算三维模型的感应-极化双场电磁响应。

本发明是这样实现的，一种时域有限差分的三维感应-极化双场数值模拟方法包括：

1)、基于迪拜模型(频率相关系数c＝1时的Cole-Cole复电阻率模型)，先获得迪拜模型电导率的频域形式，再通过逆拉普拉斯变换，得到电导率的时域表达式σ(t)；

2)、将电导率的时域表达式σ(t)代入欧姆定律中，获得欧姆定律时域卷积积分表达式，结合电导率σ(t)近似于αe^-βt形式的特点，构建电导率参数的e指数辅助方程r(t)，代入欧姆定律时域卷积积分表达式，将其转换成线性积分形式；

3)、将计算时间进行密集剖分，利用梯形积分法，获得欧姆定律时域离散递推表达式；

4)、将欧姆定律时域离散递推表达式代入Maxwell方程中，基于三维时间域有限差分方法，在时间和空间上进行中心差分离散，推导电场E(t)、磁场H(t)的迭代方程；

5)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置电性参数，时间步长，并对各个网格进行电导率参数赋值；

6)、计算初始场，对三维模型进行电场E(t)、磁场H(t)的迭代，加载狄利克雷边界条件，完成电磁场数值计算；

7)、计算结束后，提取磁场或电场的各分量响应进行成图。

进一步地，步骤2中，根据欧姆定律的频域微分形式，先通过频时变换得到欧姆定律时域卷积积分形式；再将电导率时域表达式代入，获得欧姆定律时域卷积积分表达式，如式(1)所示；根据电导率σ(t)近似于αe^-βt的形式的特点，构建电导率参数的e指数辅助方程r(t)，如式(2)所示，将r(t)的表达式代入欧姆定律的卷积积分表达式中，得到欧姆定律的线性积分形式，如式(3)所示：

式中J(t)是电流密度，E(t)为电场值，t为时间；电导率绝对值表达式为τ₁为时间常数、η为极化率，σ_∞为频率趋近无穷大时的电导率值。

进一步地，步骤3中，将计算时间进行密集剖分，利用梯形积分法，得到积分求和ξ(tⁿ)，提取ξ(tⁿ)最后一项构建递推公式，如式(4)所示，进而实现降维计算；最后将式(4)代入式(3)，并对时间进行离散，最终得到n时刻欧姆定律时域离散递推表达式，如式(5)所示：

式中n表示第n个时刻，r(tⁿ)、E(tⁿ)、ξ(tⁿ)为第n个时刻所对应的数值，Δt⁽ⁿ⁾为第n个时间步长；

接着在时刻对式(4)进行离散，再代入(5)式中，进行整理，得到时刻欧姆定律时域离散递推表达式，如式所示：

进一步地，步骤4中，利用步骤3中的推导结果，将式(6)中的项，采用平均值近似法进行计算，有代入无源Maxwell方程，采用时域有限差分方法，构建出电场E(tⁿ⁺¹)与E(tⁿ)的迭代方程；对于中的ξ(tⁿ)项，先按照步骤3中式(4)的递推公式进行计算，获得n时刻的ξ(tⁿ)值，再将其代入电场迭代公式，最终获得n+1时刻电场的迭代方程，如式(8)所示：

式中有

(i,j,k)为Yee式网格的坐标点，为n+1时刻x轴方向电场分量，E_x ⁽ⁿ⁾(m)为n时刻x轴方向电场分量，为时刻对应坐标点上的y轴方向的磁场分量，为时刻对应坐标点上的z轴方向的磁场分量，Δy_j、Δz_k为Yee氏网格的y、z方向的步长。

进一步地，步骤5中，根据空间、时间的先后顺序进行离散，先对三维Yee网格的σ_∞(x,y,z)赋值，将每个网格σ_∞(x,y,z)代入到电导率时域表达式中，得到σ(x,y,z,t)；再结合电流密度的计算时刻特点，设定时间步长Δt⁽ⁿ⁾，将时间密集剖分成N+1个时刻，计算时间半步长再对时间进行时刻剖分，对应为t⁽⁰⁾、t⁽ⁿ⁾，共2N+1个时刻，在2N+1个时刻分别对r(t)、进行赋值后，存储r(t)的2N+1个值。对于的数据，先舍弃n时刻的值，再保留时刻的值，这样只需存储N个数据。

进一步地，步骤6中，在三维模型的基础上，建立地面中心回线发射与接收线圈的系统，计算出初始场后，代入到控制方程中，计算x、y、z方向的电场与x、y方向的磁场，为了确保计算结果的唯一性，采用散度方程计算z方向的磁场，应用狄利克雷边界条件，完成三维模型的感应-极化双场电磁响应 数值计算。

本发明与现有技术相比,有益效果在于，将时域电磁探测中的感应和极化双场同时进行理论计算，通过构建辅助方程，推导欧姆定律的时域离散递推表达式，对计算场量进行降维处理，避免了卷积和的求取，大大减少了内存的使用，可以更高效、准确地模拟三维模型的感应-极化双场电磁响应。

附图说明

图1是三维有限差分数值算法的整体示意图；

图2是欧姆定律时域卷积形式的离散算法示意图；

图3是三维异常体模型的示意图。

图4是迪拜模型均匀半空间的磁场响应图；

图5是迪拜模型层状大地的磁场响应图；

图6是迪拜模型三维异常体的磁场响应图；

图7是欧姆定律时域卷积运算降维前后的计算时间对比图。

图8是欧姆定律时域卷积运算降维前后的占用内存对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。以均匀半空间模型、层状的大地模型、三维异常体模型为例，进行时域有限差分的三维感应-极化双场电磁响应计算。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。参见图1结合图2所示，一种时域有限差分的三维感应-极化双场数值模拟方法，包括：

1、基于迪拜模型(频率相关系数c＝1时的Cole-Cole复电阻率模型)，先获得迪拜模型电导率的频域形式，再通过逆拉普拉斯变换，得到时域电导率表达式σ(t)；迪拜模型的频域表达式如式(1)所示

式(1)中极化率为η，时间常数为τ₁；经拉普拉斯逆变换后得到迪拜模型时域电导率表达式为

式(2)又可表示为

式(3)中δ(t)、u(t)分别是单位脉冲函数和单位阶跃函数。

2、将电导率的时域表达式σ(t)代入欧姆定律中，获得欧姆定律时域卷积积分表达式，结合电导率σ(t)近似于αe^-βt形式的特点，构建电导率参数的e指数辅助方程r(t)，代入欧姆定律时域卷积积分表达式，将其转换成线性积分形式；

首先，将欧姆定律的频域微分形式进行频时变换后得到时域卷积形式，表达式为：

将σ(t)代入到式(4)中得欧姆定律时域卷积积分表达式为：

根据电导率σ(t)近似于αe^-βt的形式的特点，构造辅助方代入欧姆定律时域卷积积分表达式将其转换成线性积分形式，如式(6)所示：

3、对计算时间进行密集剖分，利用梯形积分法，获得欧姆定律时域离散递推表达式：

1)、将时间轴划分成N+1个时刻，则有t⁽⁰⁾、t⁽¹⁾、t⁽²⁾～t⁽ⁿ⁾这样的时间点，Δt⁽ⁿ⁾＝t⁽ⁿ⁾-t^(n-1)，Δt⁽ⁿ⁾是可变步长，随时间的增加逐渐加长；

采用梯形积分法

将各个时间步长的积分用梯形积分法近似后的表达式为

对式(9)进行拆分，提取最后一项，构建可递推公式，实现数据的降维计算，表达式为

2)、特殊时间点递推公式的求解方法，当t＝t¹时，ξ(t¹)表达式为

由于电磁场计算中E(t⁰)＝0，所以忽略掉r(t⁰)·E(t⁰)的乘积项，之后仍按照梯形积分法来近似各个时间步长上的线性积分。

3)、时域有限差分方法中以时刻为电场的观察点，为了与时域有限差分方法更好的结合，继续推导时刻方程为

4、将离散后的欧姆定律代入Maxwell方程中，利用三维时间域有限差分方法，在时间和空间上进行中心差分离散，推导得到电场E(t)、磁场H(t)的迭代方程：

1)、将步骤3的中拆分成两项，两者形式相同，表达式为

对(14)式单独进行递推计算，对(15)式采用平均值近似法把替换成项，从而构建出E(tⁿ⁺¹)与E(tⁿ)的迭代方程。

2)、将离散欧姆定律的卷积递推公式代入无源Maxwell方程，推导E、H迭代公式。无源Maxwell旋度方程组表达式及其本构关系表达式为

将旋度方程在x、y、z方向展开表达式如下：

进一步地，采用三维时域有限差分方法，将Maxwell方程组在空间和时间进行差分离散，将微商转换为差商。令f(x,y,z,t)代表直角坐标系下的E、H的某一分量，在时间和空间上离散有式(19)

f(x,y,z,t)＝f(iΔx,jΔy,kΔz,nΔt)＝f⁽ⁿ⁾(i,j,k) (19)

对函数的时间和空间的一阶偏导数取中心差分近似，推导迭代公式。

例如：选取观察点(x,y,z)为E_x的节点，即点在时间的时间节点观察，将式(12)、(14)、(15)代入后推导表达式为

式中有

同理可以推导出其他方向的分量。

同样，选取观察点(x,y,z)为H_x的节点，即点在n时刻进行观察，获得表达式为

其中同理，可以继续推导磁场其他方向的分量。

5、采用非均匀三维Yee氏网格对计算区域进行剖分，设置电性参数，时间步长，并对各个网格进行电导率参数赋值：

1)、划分非均匀Yee式网格，根据时域电磁响应在近源处幅值大且变化快，离源较远时幅值小且变化慢的特点，满足网格步长近源处密集，离源远处稀疏的条件，设置网格数目为101×101×50，其中x、y方向上的网格数目均为101个， z方向上网格数50个，三个方向上的最小和最大网格步长均为10m和120m；在计算区域内磁导率均设置为真空磁导率，均匀半空间模型σ_∞设置为0.01S/m；层状大地模型上层50m，上层设置σ_∞为0.01S/m，下层2300m，σ_∞设置为0.1S/m；三维异常体模型，如图3所示，背景场尺寸2860×2860×2350(m³)，σ_∞设置为0.01S/m，异常体尺寸1600×1600×1700(m³)，σ_∞设置为0.1S/m。

2)、初始时刻t⁽⁰⁾按照公式进行取值，其中，μ₁为最上层的磁导率，这里取真空磁导率，σ₁为最顶层的电导率σ_∞，Δ₁为Yee氏网格中最小的空间步长；在时间域有限差分方法中，一般采用

式(22)中，α的取值范围是0.1～0.2。

迭代前就对Δt⁽ⁿ⁾进行赋值并完成迭代所需的所有中间参数的计算，将时间t预先分成几个较大的时间段，近似地依照式(22)，为每一段的Δt赋值。这样的处理方式加快了迭代速度。算例中，初始时刻为1微秒，总时长为10毫秒，一共分成了15段，前14段每段有100个时间步长，最后一段有500个时间步长，其中，时间步长Δt最小为0.05微秒，最大为10微秒。

3)、根据空间、时间的先后顺序进行离散，先对三维Yee网格的σ_∞(x,y,z)赋值，将每个网格σ_∞(x,y,z)代入到电导率时域表达式中，得到σ(x,y,z,t)；再结合电流密度的计算时刻特点，设定时间步长Δt⁽ⁿ⁾，将时间密集剖分成N+1个时刻，计算时间半步长再对时间进行时刻剖分，对应为t⁽⁰⁾、t⁽ⁿ⁾，共2N+1个时刻，在2N+1个时刻分别对r(t)、进行赋值后，存储r(t)的2N+1个值。对于的数据，先舍弃n时刻的值，再保留时刻的值，这样只需存储N个数据。

6、计算初始场，对三维模型进行电场E(t)、磁场H(t)的迭代，加载狄利克雷边 界条件，完成电磁场数值计算：

在三维模型的基础上，建立地面中心回线发射与接收线圈的系统，计算出初始场后，代入到控制方程中，计算x、y、z方向的电场与x、y方向的磁场，为了确保计算结果的唯一性，采用散度方程计算z方向的磁场，应用狄利克雷边界条件，完成三维感应-极化双场电磁响应计算。最后提取z分量的磁场响应。

7、计算结束后，提取磁场或电场的各分量响应进行成图。

计算成图结果详见图4-图6。

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

Claims (5)

1.一种时域有限差分的三维感应-极化双场数值模拟方法其特征在于，包括如下步骤：

1)、基于迪拜模型(频率相关系数c＝1时的Cole-Cole复电阻率模型)，先获得迪拜模型电导率的频域形式，再通过逆拉普拉斯变换，得到电导率的时域表达式σ(t)；

2)、将电导率的时域表达式σ(t)代入欧姆定律中，获得欧姆定律时域卷积积分表达式，结合电导率σ(t)近似于αe^-βt形式的特点，构建电导率参数的e指数辅助方程r(t)，代入欧姆定律时域卷积积分表达式，将其转换成线性积分形式；

3)、将计算时间进行密集剖分，利用梯形积分法，获得欧姆定律时域离散递推表达式；

4)、将欧姆定律时域离散递推表达式代入Maxwell方程中，基于三维时间域有限差分方法，在时间和空间上进行中心差分离散，推导电场E(t)、磁场H(t)的迭代方程；

5)、采用非均匀三维Yee氏网格对计算区域进行剖分，设置电性参数，时间步长，并对各个网格进行电导率参数赋值；

6)、计算初始场，对三维模型进行电场E(t)、磁场H(t)的迭代，加载狄利克雷边界条件，完成电磁场数值计算；

7)、计算结束后，提取磁场或电场的各分量响应进行成图。

2.按照权利要求1所述的一种时域有限差分的三维感应-极化双场数值模拟方法，其特征在于：

步骤2中，根据欧姆定律的频域微分形式，先通过频时变换得到欧姆定律时域卷积积分形式；再将电导率时域表达式代入，获得欧姆定律时域卷积积分表达式，如式(1)所示；根据电导率σ(t)近似于αe^-βt的形式的特点，构建电导率参数的e指数辅助方程r(t)，如式(2)所示，将r(t)的表达式代入欧姆定律的卷积积分表达式中，得到欧姆定律的线性积分形式，如式(3)所示：

$\stackrel{\→}{J}\left(t\right)={\mathrm{\σ}}_{\mathrm{\∞}}\·\stackrel{\→}{E}\left(t\right)-{\∫}_0^t\widehat{\mathrm{\σ}}(t-\mathrm{\τ})\stackrel{\→}{E}\left(\mathrm{\τ}\right)d\mathrm{\τ}---\left(1\right)$

$r\left(t\right)=e^{\frac{t}{{\mathrm{\τ}}_1(1-\mathrm{\η})}}---\left(2\right)$

$\stackrel{\→}{J}\left(t\right)={\mathrm{\σ}}_{\mathrm{\∞}}\·\stackrel{\→}{E}\left(t\right)-\widehat{\mathrm{\σ}}\left(t\right)\·{\∫}_0^{t}r\left(\mathrm{\τ}\right)\stackrel{\→}{E}\left(\mathrm{\τ}\right)d\mathrm{\τ}---\left(3\right)$

式中J(t)是电流密度，E(t)为电场值，t为时间；电导率绝对值表达式为τ₁为时间常数、η为极化率，σ_∞为频率趋近无穷大时的电导率值。

3.按照权利要求1所述的一种时域有限差分的三维感应-极化双场数值模拟方法，其特征在于：

步骤3中，将计算时间进行密集剖分，利用梯形积分法，得到积分求和ξ(tⁿ)，提取ξ(tⁿ)最后一项构建递推公式，如式(4)所示，进而实现降维计算；最后将式(4)代入式(3)，并对时间进行离散，最终得到n时刻的欧姆定律时域离散递推表达式，如式(5)所示：

式中n表示第n个时刻，r(tⁿ)、E(tⁿ)、ξ(tⁿ)为第n个时刻所对应的数值，Δt⁽ⁿ⁾为第n个时间步长；

接着在时刻对式(4)进行离散，再代入(5)式中，进行整理，得到时刻的欧姆定律时域离散递推表达式，如式所示：

4.按照权利要求1所述的一种时域有限差分的三维感应-极化双场数值模拟方法，其特征在于：

步骤4中，利用步骤3中的推导结果，将式(6)中的项，采用平均值近似法进行计算，有代入无源Maxwell方程，采用三维时域有限差分方法，构建出电场E(tⁿ⁺¹)与E(tⁿ)的迭代方程；对于式(6)中的ξ(tⁿ)项，先按照步骤3中式(4)的递推公式进行计算，获得n时刻的ξ(tⁿ)值，再将其代入电场迭代公式，最终获得n+1时刻电场的迭代方程，如式(8)所示：

$\begin{array}{c}{E}_x^{\left(n+1\right)}\left(m\right)=\frac{2}{{\mathrm{\α}}^{\left(n+1\right)}\left(m\right)+\frac{2\mathrm{\ϵ}}{{\mathrm{\Δt}}^{\left(n+1\right)}}}\·\[\frac{H_z^{\left(n+\frac{1}{2}\right)}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)-H_z^{\left(n+\frac{1}{2}\right)}\left(i+\frac{1}{2},j-\frac{1}{2},k\right)}{\stackrel{\‾}{{\mathrm{\Δy}}_J}}\\ -\frac{H_y^{\left(n+\frac{1}{2}\right)}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)-H_y^{\left(n+\frac{1}{2}\right)}\left(i+\frac{1}{2},j,k-\frac{1}{2}\right)}{\stackrel{\‾}{{\mathrm{\Δz}}_k}}\]\\ +\frac{2\widehat{\mathrm{\σ}}\left(t^{n+\frac{1}{2}}\right)}{{\mathrm{\α}}^{\left(n+1\right)}\left(m\right)+\frac{2\mathrm{\ϵ}}{{\mathrm{\Δt}}^{\left(n+1\right)}}}\mathrm{\ξ}\left(t^{\left(n\right)}\right)-\frac{{\mathrm{\α}}^{\left(n+1\right)}\left(m\right)+2\·{\mathrm{\β}}^{\left(n+1\right)}\left(m\right)-\frac{2\mathrm{\ϵ}}{{\mathrm{\Δt}}^{\left(n+1\right)}}}{{\mathrm{\α}}^{\left(n+1\right)}\left(m\right)+\frac{2\mathrm{\ϵ}}{{\mathrm{\Δt}}^{\left(n+1\right)}}}\·{E_x}^{\left(n\right)}\left(m\right)\end{array}---\left(8\right)$

式中

(i,j,k)为Yee氏网格的坐标点，为n+1时刻x轴方向电场分量，E_x ⁽ⁿ⁾(m)为n时刻x轴方向电场分量，为时刻对应坐标点上的y轴方向的磁场分量，为时刻对应坐标点上的z轴方向的磁场分量，Δy_j、Δz_k为Yee氏网格的y、z方向的步长。

5.按照权利要求1所述的一种时域有限差分的三维感应-极化双场数值模拟方法，其特征在于：

步骤5中，根据空间、时间的先后顺序进行离散，先对三维Yee网格的σ_∞(x,y,z)赋值，将每个网格σ_∞(x,y,z)代入到电导率时域表达式中，得到σ(x,y,z,t)；再结合电流密度的计算时刻特点，设定时间步长Δt⁽ⁿ⁾，将时间密集剖分成N+1个时刻，计算时间半步长再对时间进行时刻剖分，对应为t⁽⁰⁾、 t⁽ⁿ⁾，共2N+1个时刻，在2N+1个时刻分别对r(t)、进行赋值后，存储r(t)的2N+1个值，对于的数据，先舍弃n时刻的值，再保留时刻的值，这样只需存储N个的数据。