一种电性源瞬变电磁电场响应成像方法

摘要

本发明公开了一种电性源瞬变电磁电场响应成像方法，属于瞬变电磁地球物理勘探方法技术领域。目的是解决传统的瞬变电磁快速成像方法和基于最小二乘的反演方法，在进一步提高瞬变电磁响应的解释精度方面遇到瓶颈的问题本发明提出一种对电性源瞬变电磁法电场响应进行成像的新方法，方法具体包括：(1)电性源瞬变电磁电场响应的虚拟子波解析解推导；(2)虚拟子波数值求解方法；(3)电性源瞬变电磁电场响应的速度分析成像。本发明的成像方法基于扩散场与波动场之间的数学积分变换，计算量小、计算速度快，为瞬变电磁场响应的拟波场精细解释奠定了基础。

说明书

技术领域

本发明具体涉及一种电性源瞬变电磁电场响应成像方法，属于瞬变电磁地球物理勘探方法技术领域。

背景技术

瞬变电磁法(Transient electromagnetic method，TEM)被广泛应用于煤田水患、矿产资源和地下资源探测中，是一种重要的地球物理勘探方法。电性源瞬变电磁法属于瞬变电磁法的一个分支，其采用电性源向地下注入电磁场信号，在源关断期间采集来自地下地质体的感应信号以获得其分布信息。

瞬变电磁场信号的处理和解释是瞬变电磁法的关键技术之一。通过地面采集的电磁场信号，分析地下介质的电性分布是一个典型的地球物理反问题。通常，瞬变电磁场响应通过视电阻率、“烟圈”等效反演和S反演等快速成像方法或者基于最小二乘的OCCAM、Marquadt-Levenberg反演获得地下介质的电阻率分布。快速成像方法的成像精度有限，而基于最小二乘的反演算法受到瞬变电磁场响应计算的复杂性和计算量大的限制，发展精度较高的高维反演难度较大。为此，探索新的瞬变电磁解释技术具有重要的现实意义。

发明内容

因此，本发明目的是目的在于提出一种对电性源瞬变电磁法电场响应进行成像的新方法。

具体的，本发明的方法包括以下步骤：

步骤A电性源瞬变电磁响应虚拟子波解析解的推导；

步骤A具体包括：

步骤A1将场值进行拉普拉斯变换；

步骤A2进行非线性变换s＝p²，得到波动场与扩散场之间的等价表达式；

步骤A3反拉普拉斯变换；

步骤B电性源瞬变电磁响应虚拟波场数值求解；

首先对公式一所示的积分方程进行离散；

其中t为时间，(x,y,z)为空间坐标，q被称为虚拟时间，其量纲为s^1/2；

得到线性方程组：d＝Gm， 公式二

采用正则化的方法求解公式二的线性方程组，引入的模型向量的二范数作为正则化项，并采用L曲线法求取最佳正则化因子；

步骤C基于虚拟子波峰值速度的电性源瞬变电磁数据成像；

在步骤A和步骤B的基础上，对电性源瞬变电磁场电场响应进行成像，分析虚拟子波峰值随虚拟时间的传播速度，估算地下介质的电阻率，将基于所提取的虚拟子波的峰值虚拟时刻的传播速度进行成像。

进一步的，所述步骤A的推导过程具体为：

在均匀半空间下，单位偶极源激发所产生的轴向电场响应为：

其中，c²＝ρ/μ，μ＝4π10^-7H/m，erf为误差函数。令

则公式三可以表示为：

对F₁(t)关于时间t进行拉普拉斯变换，得到

与F₁(t)所对应的虚拟波场的拉普拉斯变换为：

根据拉普拉斯变换恒等式

其中H表示Heavi side函数，令a＝r/c，对公式九进行拉普拉斯反变换，得到

对公式六中的F₂(t)进行拉普拉斯变换，我们得到

与F₂(t)所对应的虚拟波场的拉普拉斯变换为：

对式公式十三进行拉普拉斯反变换，可以得到：

将公式十一和公式十二进行结合，得到与公式三所表示的均匀半空间下轴向电场解析式所对应的虚拟波场的解析表达式：

化简，得到：

进一步的，所述步骤B中对公式一所示的积分方程进行离散具体为：

积分区间表达为[a,b]，采用中点公式，将其离散为n个积分子区间，每个子区间的中点为q₁,q₂,...,q_n，则

其中，

假设

公式一可以离散为：

其中，

G_i,j＝G(t_i,q_j)Δq>

m_j＝m(q_j)>

则我们可以得到线性方程组：

d＝Gm 公式二。

进一步的，所述步骤B中，采用奇异值分解法、Tikhonov零阶正则化、一阶正则化、截断SVD或阻尼SVD进行正则化的求解。

进一步的，所述步骤C中，电阻率和传播速度的获取方法具体为：

首先考虑均匀半空间情形，分别计算电阻率为200Ω·m,500Ω·m的均匀半空间下，单位偶极源激发所产生的下降沿阶跃响应经波场变换得到的虚拟子波；

根据虚拟子波所满足的波动方程，可知虚拟子波的传播速度为虚拟子波的峰值时刻等于；

通过计算峰值时刻随偏移距的变换关系，得到虚拟子波的传播速度，进而得到地下介质的电阻率，通过计算虚拟子波的峰值时刻对时间的偏导数，得到虚拟子波峰值时刻的传播速度。

本发明的有益效果在于：本专利提出一种对电性源瞬变电磁法电场响应进行成像的新方法，解决了传统的瞬变电磁快速成像方法和基于最小二乘的反演方法，在进一步提高瞬变电磁响应的解释精度方面遇到瓶颈的问题。该方法基于扩散场与波动场之间的数学积分变换，计算量小、计算速度快，为瞬变电磁场响应的拟波场精细解释奠定了基础。

附图说明

图1a至图1d为采用模型二范数作为正则化项的求解结果，其中图1a为零阶Tikhonov正则化解；图1b为TSVD解；图1c为DSVD解；图1d为不同截断点的TSVD解；

图2a至图2d为采用一阶模型粗糙度作为正则化项的求解结果,其中图2a为一阶Tikhonov正则化解；图2b为TSVD解；图2c为DSVD解；图2d为不同截断点的TSVD解；

图3a至图3d为采用二阶模型粗糙度作为正则化项的求解结果,其中图3a为二阶Tikhonov正则化解；图3b为TSVD解；图3c DSVD解；图3d为不同截断点的TSVD解；

图4a、图4b为不同电阻率均匀半空间的虚拟子波解析解，图4a中200 Ohm-m；图4b中500Ohm-m；

图5a、图5b为不同电阻率均匀半空间的虚拟子波数值解，图5a中200Ohm-m；图5b中500Ohm-m；

图6a、图6b为Q型模型虚拟子波数值解，图6a中ρ₁＝50Ohm-m,ρ₂＝200Ohm-m；图6b中ρ₁＝200Ohm-m,ρ₂＝100Ohm-m；

图7a、图7b为三层模型虚拟子波数值解，图7a中ρ₁＝200Ohm-m,ρ₂＝50Ohm-m,ρ₃＝200Ohm-m；图7b中ρ₁＝50Ohm-m,ρ₂＝500Ohm-m,ρ₃＝50Ohm-m。

具体实施方式

下面结合附图对本发明的具体实施方式进行说明：

瞬变电磁场数据的拟波场解释是瞬变电磁新解释技术的探索方向之一。瞬变电磁场在地下介质中的传播满足麦克斯韦方程组。受瞬变电磁信号的频带范围的限制，其以扩散方式传播，传播过程可以亥姆霍茨方程表征。Bragg(1968)和Filippi(1969)在数学上，从扩散方程和波动方程出发，推导了满足扩散方程的数学量与满足波动方程的数学量之间的数学积分关系式，若将该式应用于满足扩散场的电磁场E与满足波动方程的虚拟子波U的转换，可以得到

其中t为时间，(x,y,z)为空间坐标，q被称为虚拟时间，其量纲为s^1/2。Lee(1993)采用射线追踪的方法，利用式(1)对低频电磁场响应进行了成像，Lee(2006)利用频率域电磁场数据，获得虚拟子波U实现了频率域电磁场响应的成像处理。在国内，李貅等(2006)将回线源瞬变电磁场响应转换为虚拟子波，并采用克西霍夫偏移方法实现了回线源瞬变电磁响应的成像。

本发明提出一种新的电性源瞬变电磁响应拟波场成像方法。该方法基于式(1)所推导的扩散场与波场之间的积分转换，能够对地下介质的电阻率信息进行精确成像。具体包括：(1)电性源瞬变电磁响应虚拟子波解析解的推导；(2) 电性源瞬变电磁响应虚拟波场数值求解方法；(3)基于虚拟子波峰值速度的电性源瞬变电磁数据成像。

(模块1)：电性源瞬变电磁响应虚拟子波解析解的推导

式(1)所示积分变换为第一类Fredholm积分方程，无法直接对其进行反变换求取均匀半空间下与电磁场响应对应的虚拟子波的解析解。我们通过以下推导步骤对电性源瞬变电磁场响应的虚拟子波解析解进行推导：

1)将场值进行拉普拉斯变换；

2)进行非线性变换s＝p²，得到波动场与扩散场之间的等价表达式；

3)反拉普拉斯变换。

下面为具体推导过程：

在电性源瞬变电磁法测量中，若观源轴向方向上的电场分量。在均匀半空间下，单位偶极源激发所产生的轴向电场响应为：

其中，c²＝ρ/μ，μ＝4π10^-7H/m，erf为误差函数。令

则式(2)可以表示为：

对F₁(t)关于时间t进行拉普拉斯变换，得到

相应地，我们得到与F₁(t)所对应的虚拟波场的拉普拉斯变换为：

根据拉普拉斯变换恒等式

其中H表示Heavi side函数，令a＝r/c，对式(8)进行拉普拉斯反变换，可以得到

同样地，对式(5)中的F₂(t)进行拉普拉斯变换，我们得到

相应地，我们得到与F₂(t)所对应的虚拟波场的拉普拉斯变换为：

对式(12)进行拉普拉斯反变换，可以得到：

将式(10)和式(11)进行结合，即可得到与式(2)所表示的均匀半空间下轴向电场解析式所对应的虚拟波场的解析表达式：

化简，得到：

(模块2)：电性源瞬变电磁响应虚拟子波数值求解方法

本模块面向电性源瞬变电磁法虚拟子波提取。模块1推导了均匀半空间虚拟子波的解析解。然而在层状模型或者二、三维模型下，无法求得虚拟子波的解析解，需要采用数值方法进行求解。本模块给出了电性源瞬变电磁响应虚拟子波数值求解的最优方法。首先对式(1)所示的积分方程进行离散。在实际算例中，积分区间是有限的，将其表达为[a,b]。采用中点公式，将其离散为n个积分子区间，每个子区间的中点为q₁,q₂,...,q_n，则

其中，

若假设

此时，式(1)可以离散为：

其中，

G_i,j＝G(t_i,q_j)Δq>

m_j＝m(q_j)>

则我们可以得到线性方程组：

d＝Gm (22)

因此，虚拟子波的求解问题转换为了式(23)所示线性方程组进行求解的过程。积分表达式(1)的病态性使得线性方程组(23)的求解过程具有病态性，需要采用正则化的方法求解。为此，引入的模型向量的二范数作为正则化项，并采用L曲线法求取最佳正则化因子。正则化的求解过程采用奇异值分解法(Singular value decomposition,SVD)进行。考虑了Tikhonov零阶正则化、一阶正则化、截断SVD(truncated SVD,TSVD)和阻尼SVD(Damped SVD)等方法进行求解。下面为求解结果：

图1a至图1c分别给出了引入模型二范数作为正则化项，分别采用Tikhonov正则化，TSVD和DSVD的求解结果，图中纵坐标表示振幅，曲线表示真实虚拟子波，数据点表示数据求解结果。对于Tikhonov正则化和DSVD，正则化因子由L曲线获得；对于TSVD，其截断点的位置由L曲线获得，在图中均以λ标识为正则化因子。图1a中零阶Tikhonov正则化求解结果的晚期稳定，但是早期偏离真实值较大。图1c中DSVD求解结果在早期和晚期均不稳定。图1b中的TSVD求解结果的稳定性要优于另外两种方法，由于奇异值的截断，早期的解的不稳定性得到了压制，但是解的晚期出现了轻微的振荡。为此，计算了由L曲线获得的最佳正则化参数附近的TSVD解，将其绘制于图1b中。随着截断点后延，考虑的奇异值增多，所得到的TSVD解早期的稳定性降低。但是，所获得的虚拟子波的幅值更接近与真实值。相反地，随着截断点提前，考虑的奇异值减小，得到的TSVD解总体而言稳定性增加。总体而言，采用模型二范数作为正则化项时，所得到的正则化解早期的稳定性较差。

图2和图3分别一阶模型粗糙度和二阶模型粗糙度作为正则化项得到的虚拟子波提取结果。当采用高阶模型粗糙度时，需要采用广义SVD(Generalized SVD)方法进行求解。由图可知，当采用模型粗糙度作为正则化项时，Tikhonov方法和TGSVD方法均能得到的稳定的解。但是，采用阻尼奇异值分解得到的解的效果仍然不理想。通过对比Tikhonov方法和TGSVD方法所获得的解，认为TGSVD对虚拟子波的峰值和峰值时刻的提取效果更好。另外，通过对比图2d和图3d的TGSVD的提取结果可知，采用二阶模型粗糙度作为正则化项时，所提取的虚拟子波的峰值和峰值时刻更为准确。当采用一阶模型粗糙度作为正则化项时，所提取的虚拟子波的晚期振荡较弱。

通过上述分析可知，波场变换的数值解并不能准确的恢复出虚拟子波，而仅是虚拟子波的一个近似。求解病态反问题的难点在于反问题本身，而非所采用的正则化算法。因此，对病态反问题的理解和对解的认识往往比单纯求得反问题的稳定解更为重要。本模块通过采用均匀半空间下的瞬变电磁响应及其相对应的虚拟子波作为参照，给出了不同正则化方法的解的特征。在后续的成像 中，所关心的是所提取的虚拟子波的峰值时刻，因此在这些正则化方法中，采用二阶模型粗糙度并采用TGSVD所获得的解的特性更符合成像需要。

(模块3)：电性源瞬变电电磁法虚拟波场成像

本模块在模块1和模块2的基础上，对电性源瞬变电磁场电场响应进行成像。在瞬变电磁场数据解释中，成像的目的是获得地下介质的电阻率特性。均匀半空间下，虚拟子波随虚拟时间的传播速度与地下介质的电阻率相关。通过分析虚拟子波峰值随虚拟时间的传播速度可以估算地下介质的电阻率。为了借鉴地震数据处理中的成像技术，本模块将基于所提取的虚拟子波的峰值虚拟时刻的传播速度进行成像。

首先考虑均匀半空间情形，分别计算电阻率为200Ω·m,500Ω·m的均匀半空间下，单位偶极源激发所产生的下降沿阶跃响应经波场变换得到的虚拟子波。图4为偏移距1000m至7000m范围内虚拟子波的解析解，图5为相对应的数值解，图4中线段表示所拾取的虚拟子波的峰值时刻随偏移距的变化，图5中虚线对应图4中的线段。由于此处考虑的是虚拟子波的峰值时刻，因此对图4和图5的虚拟子波的幅值进行了归一化处理。根据虚拟子波所满足的波动方程，可知虚拟子波的传播速度为另外，均匀半空间下，虚拟子波的解析表达式表明，虚拟子波的峰值时刻等于

由此可知，通过计算峰值时刻随偏移距的变换关系，即可得到虚拟子波的传播速度，进而得到地下介质的电阻率。通过计算线段所示虚拟子波的峰值时刻对时间的偏导数，即可得到虚拟子波峰值时刻的传播速度。表1给出了不同电阻率的均匀半空间下虚拟子波的传播速度及相对应的电阻率。

表1

由表1可知，在不同的电阻率情形下，通过计算虚拟子波的峰值时刻的传 播速度，可以得到较为准确的地下介质的电阻率。在低阻介质中，虚拟子波的峰值时刻的传播速度较小，在高阻介质中虚拟子波峰值时刻的传播速度较大，这与电磁场在地下介质的扩散速度的基本规律是一致的。

数值算例

以数值算例来说明本专利所提出的电性源瞬变电磁电场响应成像方法及其应用效果。采用了两层模型和三层模型验证了本专利提出的成像方法的有效性。

两层模型算例。图6给出了Q型和D型模型虚拟子波的数值解。图6a所示虚拟子波由参数为ρ₁＝50Ω·m，h₁＝800m，ρ₂＝200Ω·m的Q型模型的瞬变电磁场响应生成，图6b所示虚拟子波由参数为ρ₁＝200Ω·m，h₁＝800m，ρ₂＝100Ω·m的D型模型的瞬变电磁场响应生成。图中线段表示所拾取虚拟子波的峰值时刻随偏移距的变化。与均匀半空间不同，Q型模型下所拾取的线段分为两段。由图6a计算两段虚拟子波的峰值时刻的传播速度分别为和计算得到相应的电阻率为53.89Ω·m和208.24Ω·m。由图6a计算两段虚拟子波的峰值时刻的传播速度为和相对应的两层模型的电阻率为200.52Ω·m和107.79Ω·m。计算得到的电阻率能够很好的与模型的真实电阻率吻合。

三层模型算例。图7a和图7b分别给出了H型模型和K型模型的虚拟子波的数值解。H型模型的参数为ρ₁＝200Ω·m，h₁＝800m，ρ₂＝50Ω·m，h₂＝200m,ρ₃＝200Ω·m，K型模型的参数为ρ₁＝50Ω·m，h₁＝800m，ρ₂＝500Ω·m，h₂＝200m,ρ₃＝50Ω·m。由图7可知，所拾取的虚拟子波的峰值时刻随偏移距的变换可以分为三条线段，H型模型相对应的虚拟子波的峰值时刻的传播速度分别为和计算得到的响应的电阻率为204.86Ω·m，53.72Ω·m和195.52Ω·m。K型模型相对应的虚拟子波的峰值时刻的传播速度分别为和计算得到的电阻率分别为53.74Ω·m，488.67Ω·m和190.52Ω·m。H型模型与K型模型各层的电阻率均得到了较好的恢复。由此可知，本发明所提出的成像方法能够对地下介质的电阻率进行较精确成像。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。