CN107657137B - 一种有理函数逼近的分数阶电磁反常扩散三维模拟方法 - Google Patents

Abstract

本发明涉及一种有理函数逼近的分数阶电磁反常扩散三维模拟方法，目的在于计算分数阶科尔‑科尔模型的三维时域感应‑极化双场响应。主要包括基于频域有理函数逼近法，构建科尔‑科尔模型分数阶传递函数和n阶有理逼近函数，将误差函数实、虚部绝对值之和作为目标函数；通过辅助变量法，实现目标函数的线性化，采用线性规划方法获得最佳逼近有理函数；采用部分分式展开法和拉普拉斯逆变换获得电导率的时域形式；将其代入Maxwell方程，基于有限差分方法推导电磁场的迭代方程，实现分数阶科尔‑科尔模型三维电磁响应数值计算。本发明有益效果在于，快速准确地模拟了分数阶柯尔‑柯尔模型的三维时域电磁响应，为研究极化介质中电磁反常扩散提供了理论依据。

Description

一种有理函数逼近的分数阶电磁反常扩散三维模拟方法

技术领域

本发明涉及一种地球物理正演模拟领域中的含极化效应的时域电磁场反常扩散数值模拟方法，尤其适用时间域快速计算分数阶科尔-科尔模型三维空间感应-极化双场电磁响应。

背景技术

瞬变电磁法(Transient-electromagnetic method简称TEM)是一种以电磁感应原理为理论基础，采用人工场源为激励源的时域电磁探测方法。瞬变电磁法具有对低阻体反应灵敏，探测精度高等优势，已逐渐成为一种极具发展前景的方法，被广泛应用于水文和矿产探测领域中。然而在应用中心回线或重叠回线进行实际测量时，电磁波在传播过程存在反常扩散现象，具体表现为电磁响应符号的反转，而符号反转现象给数据的解释与反演带来了很大的困难。经过国内外专家学者大量研究发现符号反转现象是由激发极化效应引起的。激发极化效应(Induced Polarization effect，简称IP effect)是由地下导体在外加电场时阴阳离子定向移动产生的电化学反应，具体表现为，外加电场时，阴阳离子定向移动，产生极化电场，外加电场撤掉后，阴阳离子重新恢复杂乱无章的初始状态，极化电场随之衰减。国内外专家学者针对极化效应对电磁响应的影响进行了大量的研究。

中国专利CN105676295A公开了一种基于SQUID的磁源激发极化-感应的联合探测系统与方法，系统利用低温SQUID传感器和接收电路组成接收装置对磁源激发极化-感应电磁响应进行联合探测，首先铺设发射回线，将无感电阻串联在发射回线中，接着实现发射单元与接收单元同步传输，最后基于迪拜(整数阶)模型极化特性，将一次场剔除后，进行电阻率-深度成像。

中国专利CN105893678A公开了一种时域有限差分的三维感应-极化双场数值模拟方法，首先采用拉普拉斯逆变换获得迪拜(整数阶)模型电导率的时域表达式，构建电导率参数的e指数辅助方程，通过梯形积分法获得欧姆定律时域离散递推表达式，再将其代入无源Maxwell旋度方程中，基于三维时域有限差分方法推导电场和磁场的迭代方程，进而完成三维模型的感应-极化双场电磁响应数值计算。该方法可有效实现存在解析解的迪拜模型正演计算，而无法直接应用于分数阶科尔-科尔模型的计算中来。

中国专利CN106776478A公开了针对流体扩散浓度问题的一种反常扩散中的基于分步计算的离散分数阶差分方法，首先获取相关参数及初边条件；其次，采用离散分数阶差分方法离散控制方程，得到时间空间反常扩散方程的离散格式；然后将两个gamma函数的比值看作一个参数，利用递推关系，结合分布计算的思路和整体考虑的思想，将两个gamma函数比值转化为多个小数乘积的形式，代入离散控制方程，得到扩散浓度的数值结果。然而计算模型为一维模型，未实现对三维空间模型的数值模拟。

欧洲专利EP2102688公开了一种基于伪随机二进制序列(PRBS)电流源来诱导极化效应的步骤和装置，首先在理论上采用科尔-科尔模型来模拟极化介质的频散特性，接着于频域内进行数值计算与特性分析，最后提出了一种诱导极化效应发生的装置与探测方法。

美国专利US7529627公开了一种对海底深处的频散介质进行电磁测量的方法。利用科尔-科尔复电阻模型对频散介质进行等效，计算了层状模型的电磁响应，并提出了一套磁偶极子测量设备与探测方法。通过对主、二次场的信号进行修正，对背景电阻率和地下极化特性进行瞬态场分析，从而提高了极化介质的预测精度。然而计算模型为简单层状模型，未实现对含极化效应的三维异常体模型的数值计算。

以上所述方法公布了基于科尔-科尔模型在频域内对层状极化介质进行数值计算的方法以及基于整数阶迪拜模型的时域电磁响应计算方法，国内外专利还未涉及基于分数阶科尔-科尔复电阻率模型在时域内直接模拟三维电磁反常扩散过程的方法，为此，本发明基于时域有限差分法方法直接计算三维异常体模型的分数阶科尔-科尔模型的感应-极化双场电磁响应，采用有理函数逼近方法和线性规划的方法避免了对Maxwell方程组分数阶微分方程以及函数极值复杂的求解过程，并且提高了三维感应-极化双场电磁响应的计算效率。

发明内容

本发明所要解决的关键问题是提供一种有理函数逼近的分数阶电磁反常扩散三维模拟方法，通过频域有理函数逼近法和线性规划方法求解分数阶科尔-科尔模型传递函数的最佳逼近有理传递函数；采用部分分式展开法和拉普拉斯逆变换获得电导率时域形式，发明可以克服有限差分方法的分数阶微分方程和拟合过程函数极值求解困难的问题，适用于任意分数阶科尔-科尔模型的感应-极化双场快速计算。

本发明是这样实现的，一种有理函数逼近的分数阶电磁反常扩散三维模拟方法包括：

1)、基于频域有理函数逼近方法，构建科尔-科尔模型的分数阶传递函数和n阶有理传递函数，根据计算精度需求，合理设置有理传递函数的分子、分母阶数m、n，拟合频带以及频点值，将误差函数的实部与虚部绝对值之和作为目标函数；

2)、利用辅助变量法，实现目标函数的线性化，并且基于绝对值特征关系获得线性约束条件，应用线性规划求解频域误差极小化问题，提取误差函数计算结果；

3)、判断误差函数值是否达到精度要求，若未达到要求，则反复调整有理传递函数的分子、分母阶数或扩大目标函数的参数倍数，重复步骤2)，直到满足精度要求为止；

4)、提取最佳有理传递函数系数，获得有理传递函数表达式，采用部分分式展开法和拉普拉斯逆变换获得有理传递函数时域形式，得到科尔-科尔模型电导率的时域表达式；

5)、将科尔-科尔模型电导率的时域表达式代入欧姆定律卷积方程中，基于同底数幂相乘的原理，对时间变量t与积分变量τ进行分离，获得欧姆定律的线性积分表达式；

6)、对计算时间进行离散，采用梯形积分法对每一个时间间隔进行数值积分，将所有时间间隔的积分值进行叠加，再通过提取最后一个积分项构建欧姆定律的时域离散递推表达式；

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

8)、采用非均匀三维Yee氏网格对求解域进行空间剖分，设置电性参数，时间步长，空间间隔和迭代步数；

9)、基于三维异常体模型，进行电场E(t)、磁场H(t)和电流密度J(t)的交替迭代计算，加载CPML吸收边界条件，完成三维时域感应-极化双场电磁响应快速数值计算；

进一步地，步骤1)中，针对复杂多孔极化介质频散特性，可以通过0<c<1的任意分数阶科尔-科尔模型进行等效，提取电导率表达式中的分数阶项作为待拟合分数阶传递函数G(jω)：

其中ω是角频率，σ_∞是频率无穷大时极化介质的电导率，τ是时间常数，η和c分别是极化率与频散系数，0<η<1，0<c<1；

接着构建n阶有理传递函数表达式，根据计算精度需求，合理设置有理传递函数中分子最高阶数m和分母最高阶数n，通常情况下，为了取得良好的拟合效果，取m＝n或m＝n-1；再次，选定拟合频段为[ω_L,ω_H]，从中选取k个频点(ω₁,ω₂,ω₃,…ω_k)并分别计算G(jω)各个频点上实部值与虚部值；最后，基于频域误差极小化方法，将待拟合分数阶传递函数与有理传递函数误差函数的实、虚部绝对值之和作为目标函数。

进一步地，步骤2)中，对于任意的绝对值|x|都有u>0、v>0满足条件|x|＝u+v，x＝u-v，因此构造辅助变量U＝[u₁,u₂,…,u_k],V＝[v₁,v₂,…,v_k]，S＝[s₁,s₂,…,s_k],Z＝[z₁,z₂,…,z_k]，并令它们满足关系式：

其中α(ω)，β(ω)为有理传递函数的分子的实数函数，

和ψ(ω)、为分母的实数函数；α_i,β_i,

ψ_i分别是α(ω_i)、β(ω_i)、

和ψ(ω_i)的缩写形式；根据以上条件可以将非线性目标函数转化为线性目标函数，并且基于绝对值特征关系获得线性约束条件，线性规划目标函数与约束条件方程为：

应用线性规划求解频域误差极小化问题，提取误差函数结果；并且按照步骤3)对逼近效果进行评估，最终获得最佳逼近有理传递函数的分子、分母系数。

进一步地，步骤4)中，分别提取最佳有理传递函数分子、分母系数，通过部分分式展开法获得最佳有理传递函数的频域分式叠加形式，再通过拉普拉斯逆变换，得到时域表达式

其中r₁,r₂,…,r_n为各分式前系数，p₁,p₂,…,p_n为弛豫时间；

结合科尔-科尔模型电导率完整表达式，最终获得电导率时域表达式：

其中r是三维空间坐标；

进一步地，步骤5)中，已知Maxwell方程中欧姆定律表达式，将科尔-科尔模型电导率的时域表达式代入欧姆定律中获得欧姆定律时域卷积形式：

当激励源为磁性源激励且发射波形为阶跃波关断时，有E(t)＝0，t≤0，欧姆定律时域卷积表达式可写为：

为消除卷积的时间交错相乘项，基于同底数幂相乘的性质，将含时间变量t的

指数项从积分函数

中分离出来，成为两个指数函数，并且分别定义为κ_j(τ),

最后对积分变量τ单独进行积分，构建欧姆定律线性积分公式为：

式中J(r,t)表示电流密度；E(r,τ)表示电场。

本发明与现有技术相比，有益效果在于，基于分数阶科尔-科尔模型对时域电磁探测中的感应-极化双场同时进行理论计算，计算结果更符合实际极化介质的频散特性；通过频域有理函数逼近法和线性规划方法求解科尔-科尔模型分数阶传递函数的最佳逼近有理传递函数；发明可以克服有限差分方法对分数阶微分方程和拟合过程中函数极值求解困难的问题，适用于任意分数阶科尔-科尔模型的感应-极化双场快速计算。

附图说明

图1是有理函数逼近的分数阶电磁反常扩散三维数值模拟算法整体示意图；

图2是频域误差极小化方法进行分数阶系统逼近的原理图；

图3是最佳有理传递函数与科尔-科尔模型分数阶传递函数伯德图；

图4是不同拟合阶数、不同频散系数的拟合效果对比图；

图5是三维异常体模型的示意图；

图6是三维感应-极化双场电磁响应计算迭代步骤流程图；

图7是本发明分数阶科尔-科尔模型三维异常体的磁场响应计算结果与辅助微分方程方法计算结果对比图。

具体实施方式

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

参见图2结合图1所示，一种有理函数逼近的分数阶电磁反常扩散三维模拟方法，包括：

1、基于频域有理函数逼近方法，构建科尔-科尔模型的分数阶传递函数和n阶有理传递函数，根据计算精度需求，合理设置有理传递函数的分子、分母阶数m、n，拟合频带以及拟合频点值，将误差函数的实部与虚部绝对值之和作为目标函数；

针对复杂多孔极化介质频散特性，可以通过0<c<1的任意分数阶科尔-科尔模型进行等效，获得电导率频域表达式为：

提取电导率表达式中的分数阶项作为待拟合的分数阶传递函数G(jω)：

其中ω是角频率，σ_∞是频率无穷大时极化介质的电导率，τ是时间常数，η和c分别是极化率与频散系数，0<η<1，0<c<1。

接着构建n阶有理传递函数为：

其中m为有理传递函数中分子最高阶数，n为分母最高阶数，c_i,d_j(i＝1,2,…n,j＝1,2,…m)为有理传递函数分子、分母系数；D(ω)，C(ω)为传递函数的分子、分母表达式，α(ω)，β(ω)为有理传递函数的分子的实数函数，

和ψ(ω)为分母的实数函数。

根据计算精度需求，合理设置m、n的值，通常情况下，为了取得良好的拟合效果，取m＝n或m＝n-1；接着，选定拟合频段为[ω_L,ω_H]，从中选取k个频点(ω₁,ω₂,ω₃,…ω_k)并分别计算G(jω)各个频点上实部值与虚部值；最后采用如图2所示的频域误差极小化方法求解分数阶传递函数与有理传递函数的误差表达式为：

式(4)等号左右两边同时乘C(jω_i)，并将G(jω_i)展开获得误差函数：

其中α_i,β_i,

ψ_i分别是α(ω_i)、β(ω_i)、

和ψ(ω_i)的缩写形式。

为了克服最小二乘拟合方法求解函数极值时复杂困难的问题，取为误差函数各频点实部绝对值与虚部绝对值之和作为目标函数：

2、利用辅助变量法，实现目标函数的线性化，并且基于绝对值特征关系获得线性约束条件，应用线性规划求解频域误差极小化问题，提取误差函数计算结果；

数学上，对于任意的绝对值|x|都有u>0、v>0满足条件|x|＝u+v，x＝u-v，因此构造辅助变量U＝[u₁,u₂,…,u_k],V＝[v₁,v₂,…,v_k]，S＝[s₁,s₂,…,s_k],Z＝[z₁,z₂,…,z_k]，并令它们满足关系式：

根据以上条件可以将非线性目标函数转化为线性目标函数，并且基于绝对值特征关系获得线性约束条件，线性规划目标函数与约束条件方程为：

应用线性规划求解频域误差极小化问题，提取误差函数结果。

3、判断误差函数是否达到精度要求，若未达到要求，则反复调整有理传递函数的分子、分母阶数或扩大目标函数的参数倍数，重复步骤2)，直到满足精度要求为止；

通过调整参数对分数阶传递函数进行多次逼近后获得最佳逼近有理传递函数，最佳逼近有理传递函数与科尔-科尔模型分数阶传递函数伯德图如图3所示，通过对比分析可以得出在幅值和相位的拟合效果均很好。另外，还给出不同拟合阶数、不同频散系数情况下的拟合效果对比图，如图4所示，通过分析可以得出，相同拟合阶数下，随着频散系数变小，拟合效果逐渐变差；而相同频散系数下，有理传递函数的阶数越高，拟合效果越好。因此对于频散系数较大时，我们可以选择较低的拟合阶数，当频散系数较小时，可以选择较高的拟合阶数。

4、提取最佳有理传递函数系数，获得分数阶科尔-科尔模型的频域表达式，采用部分分式展开法和拉普拉斯逆变换获得有理传递函数时域形式，得到科尔-科尔模型电导率的时域表达式；

首先，提取最佳逼近有理传递函数系数c_i,d_j(i＝1,2,…n,j＝1,2,…m)，采用部分分式展开法将传递函数表达为：

其中r₁,r₂,…,r_n为各分式前系数，p₁,p₂,…,p_n为弛豫时间。

利用逆拉普拉斯变换获得时域表达式：

结合科尔-科尔模型电导率完整表达式，最终获得电导率时域表达式σ(r,t)：

其中r是三维空间坐标；

5、将科尔-科尔模型电导率的时域表达式代入欧姆定律卷积方程中，基于同底数幂相乘的原理，对时间变量t与积分变量τ进行分离，获得欧姆定律的线性积分表达式；

已知磁性源激励发射电流为阶跃波关断时，Maxwell方程中欧姆定律表达式，将科尔-科尔模型电导率的时域表达式σ(r,t)代入时域欧姆定律中，获得欧姆定律时域卷积积分表达式：

为消除卷积的时间交错相乘项，基于同底数幂相乘的原理，将含时间变量t的

指数项从积分函数

中分离出来，成为两个指数函数，并且分别定义为κ_j(τ),

为：

其中r是三维空间坐标；

最后对积分变量τ单独进行积分，构建欧姆定律线性积分公式如下：

式中J(r,t)表示电流密度；E(r,τ)表示电场。

6、对计算时间进行离散，采用梯形积分法对每一个时间间隔进行数值积分，将所有时间间隔的积分值进行叠加，再通过提取最后一个积分项构建欧姆定律的时域离散递推表达式；

将时间轴剖分成l+1个时刻，则有t⁽⁰⁾,t⁽¹⁾,t⁽²⁾～t^(l)这样的时间点，Δt^(l)＝t^(l)-t^(l-1)，Δt^(l)是可变步长，随时间的增加逐渐加长；在每个时间间隔应用梯形积分公式：

则所有时间间隔的积分值进行叠加获得欧姆定律离散表达式：

再通过提取最后一个积分项构建欧姆定律的时域离散递推表达式：

综合式(16)和式(17)获得欧姆定律时域离散递推表达式为：

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

无源Maxwell旋度方程组表达式及其本构关系表达式为：

D(r,t)＝ε(r)E(r,t) (21a)

B(r,t)＝μ(r)H(r,t) (21b)其中，B(r,t)、H(r,t)和D(r,t)分别是磁感应强度、磁场和电感应强度；ε(r)和μ(r)分别是介电常数和磁导率。

采用三维时域有限差分方法，将Maxwell方程组在空间和时间进行差分离散，将微商转换为差商。对函数的时间和空间的一阶偏导数取中心差分近似，推导迭代公式。

例如：选取观察点(x,y,z)为E_x的节点，即(i+1/2,j,k)点在时间t＝l+1/2的时间节点观察，电场迭代方程为：

其中，

同样，选取观察点(x,y,z)为H_x的节点，即(i,j+1/2,k+1/2)点在l时刻进行观察，磁场迭代方程为

8、采用非均匀三维Yee氏网格对求解域进行空间剖分，设置电性参数，时间步长，空间间隔和迭代步数；

1)、采用非均匀Yee式网格划分计算区域，设置网格数目为101×101×50，其中x、y方向上的网格数目均为101个，z方向上网格数50个，根据时域电磁响应在近源处幅值大且变化快，离源较远时幅值小且变化慢的特点，满足网格步长近源处密集，离源远处稀疏的条件，三个方向上的最小和最大网格步长均为10米和120米；设置三维异常体模型(如图5所示)的电性参数，磁导率均设置为真空磁导率，背景场电导率设置为0.01西门子/米，异常体放置在中心位置，埋深为100米，尺寸为610米×610米×600米，电导率设置为0.1西门子/米；频散系数c＝0.5，极化率η＝0.1，时间常数τ＝0.01秒。

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

其中，α取值范围为0.1～0.2。

3)、设置发射-接收装置为中心回线方式，发射电流大小为30安培，电流波形为阶跃信号，发射线圈半径为7.5米。

9、基于三维异常体模型，进行电场E(t)、磁场H(t)和电流密度J(t)的交替迭代计算，加载CPML吸收边界条件，完成三维时域感应-极化双场电磁响应快速数值计算；

在三维模型的基础上，建立地面中心回线发射-接收系统，计算出初始场后，代入到控制方程中，每次迭代先求解x、y、z方向的电场后，再迭代计算当前时刻的电流密度，最后计算x、y方向的磁场，为了确保计算结果的唯一性，采用散度方程

计算z方向的磁场，最后加载CPML吸收边界条件，完成三维感应-极化双场电磁响应计算，具体的迭代步骤如图6所示。

计算成图结果详见图7。

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

Claims (5)

1.一种有理函数逼近的分数阶电磁反常扩散三维模拟方法，其特征在于，包括如下步骤：

1)、基于频域有理函数逼近方法，构建科尔-科尔模型的分数阶传递函数和n阶有理传递函数，根据计算精度需求，合理设置有理传递函数的分子、分母阶数m、n，拟合频带以及频点值，将误差函数的实部与虚部绝对值之和作为目标函数；

2)、利用辅助变量法，实现目标函数的线性化，并且基于绝对值特征关系获得线性约束条件，应用线性规划求解频域误差极小化问题，提取误差函数计算结果；

3)、判断误差函数值是否达到精度要求，若未达到要求，则反复调整有理传递函数的分子、分母阶数或扩大目标函数的参数倍数，重复步骤2)，直到满足精度要求为止；

4)、提取最佳有理传递函数系数，获得有理传递函数表达式，采用部分分式展开法和拉普拉斯逆变换获得有理传递函数时域形式，得到科尔-科尔模型电导率的时域表达式；

5)、将科尔-科尔模型电导率的时域表达式代入欧姆定律卷积方程中，基于同底数幂相乘的原理，对时间变量t与积分变量τ进行分离，获得欧姆定律的线性积分表达式；

6)、对计算时间进行离散，采用梯形积分法对每一个时间间隔进行数值积分，将所有时间间隔的积分值进行叠加，再通过提取最后一个积分项构建欧姆定律的时域离散递推表达式；

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

8)、采用非均匀三维Yee氏网格对求解域进行空间剖分，设置电性参数，时间步长，空间间隔和迭代步数；

9)、基于三维异常体模型，进行电场E(t)、磁场H(t)和电流密度J(t)的交替迭代计算，加载CPML吸收边界条件，完成三维时域感应-极化双场电磁响应快速数值计算。

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

步骤1)中，针对复杂多孔极化介质频散特性，可以通过0<c<1的任意分数阶科尔-科尔模型进行等效，提取电导率表达式中的分数阶项作为待拟合分数阶传递函数G(jω)：

其中ω是角频率，σ_∞是频率无穷大时极化介质的电导率，τ是时间常数，η和c分别是极化率与频散系数，0<η<1、0<c<1；

接着构建n阶有理传递函数表达式，根据计算精度需求，合理设置有理传递函数中分子最高阶数m和分母最高阶数n，取m＝n或m＝n-1；再次，选定拟合频段为[ω_L,ω_H]，从中选取k个频点(ω₁,ω₂,ω₃,…ω_k)并分别计算G(jω)各个频点上实部值与虚部值；最后，基于频域误差极小化方法，将待拟合分数阶传递函数与有理传递函数误差函数的实、虚部绝对值之和作为目标函数。

3.按照权利要求1所述的一种有理函数逼近的分数阶电磁反常扩散三维模拟方法，其特征在于：

步骤2)中，对于任意的绝对值|x|都有u>0、v>0满足条件|x|＝u+v，x＝u-v，因此构造辅助变量U＝[u₁,u₂,…,u_k],V＝[v₁,v₂,…,v_k]，S＝[s₁,s₂,…,s_k],Z＝[z₁,z₂,…,z_k]，并令它们满足关系式：

其中α(ω)，β(ω)为有理传递函数的分子的实数函数，

和ψ(ω)为分母的实数函数；α_i,β_i,

ψ_i分别是α(ω_i)、β(ω_i)、

和ψ(ω_i)的缩写形式；根据以上条件可以将非线性目标函数转化为线性目标函数，并且基于绝对值特征关系获得线性约束条件，线性规划目标函数与约束条件方程为：

应用线性规划求解频域误差极小化问题，提取误差函数结果；并且按照步骤3)对逼近效果进行评估，最终获得最佳有理传递函数的分子、分母系数。

4.按照权利要求1所述的一种有理函数逼近的分数阶电磁反常扩散三维模拟方法，其特征在于：

步骤4)中，分别提取最佳有理传递函数分子、分母系数，通过部分分式展开法获得最佳有理传递函数的频域分式叠加形式，再通过拉普拉斯逆变换，得到时域表达式

其中r₁,r₂,…,r_n为各分式前系数，p₁,p₂,…,p_n为弛豫时间；

结合科尔-科尔模型电导率完整表达式，最终获得电导率时域表达式：

其中r是三维空间坐标。

5.按照权利要求1所述的一种有理函数逼近的分数阶电磁反常扩散三维模拟方法，其特征在于：

步骤5)中，已知Maxwell方程中欧姆定律表达式，将科尔-科尔模型电导率的时域表达式代入欧姆定律中获得欧姆定律时域卷积形式：

当激励源为磁性源激励且发射波形为阶跃波关断时，有E(t)＝0，t≤0，欧姆定律时域卷积表达式可写为：

为消除卷积的时间交错相乘项，基于同底数幂相乘的原理，将含时间变量t的

指数项从积分函数

中分离出来，成为两个指数函数，并且分别定义为κ_j(τ),

最后对积分变量τ单独进行积分，构建欧姆定律线性积分公式为：

式中J(r,t)表示电流密度，E(r,τ)表示电场。

Images