CN113960674A - 一种广域电磁法二维反演方法 - Google Patents

Abstract

本发明提供一种广域电磁法二维反演方法，包括以下步骤：S1、对实测的广域电磁法电阻率曲线进行一维有源反演，得到各测点的一维有源反演模型；S2、根据一维有源反演模型，计算其一维MT响应，得到MT电阻率曲线，即，相当于将广域电磁法电阻率曲线校正到MT电阻率曲线；S3、对各个校正后的MT电阻率曲线，合成MT剖面，进行MT二维反演，得到二维反演的地电断面图；S4、利用步骤S3得到的二维反演的地电断面图进行地质构造的解译。本发明通过一维模型的过渡，将复杂的广域电磁法(WFEM)电阻率曲线校正成MT(大地电磁测深法)电阻率曲线，进而减少广域电磁法二维反演计算工作量。

Description

一种广域电磁法二维反演方法

技术领域

本发明涉及地球物理技术领域，具体是一种广域电磁法二维反演方法。

背景技术

广域电磁法(WFEM)是近几年快速发展起来的一种人工源电磁勘探方法，主要采用E-E的工作方式，即采用电性源，接收电场的方式进行工作。近年来，E-E广域电磁法在油气、矿产、地热、水文、工程等方面的应用越来越多。E-E广域电磁法计算的视电阻率曲线和CSAMT卡尼亚视电阻率曲线、MT电阻率曲线等形态差异较大，也就是说，广域电磁法电阻率，不是CSAMT(可控源音频大地电磁法)电阻率，也不是MT(大地电磁法)的电阻率。但现阶段广域电磁法二维反演方法尚不成熟，很多情况下是将广域电阻率曲线直接用MT的方法进行二维反演。理论分析表明，在低阻背景区，上述反演方法基本能得到较为吻合的地电断面图，但在高阻区，则基本上得不到好的结果，甚至得到的是错误结果。同时从理论上来讲，有源的二维反演，计算量较大，在计算方法没有突破的情况下，反演时间将很长，如何快速反演，也是当前地球物理工作者应该考虑的问题。

发明内容

针对上述现有技术的不足，本发明提供了一种广域电磁法二维反演方法，旨在较快、较准确地对广域电磁法数据进行二维反演，得到较为接近实际情况的地电断面图。

为了解决上述技术问题，本发明采用了如下的技术方案：

一种广域电磁法二维反演方法，包括以下步骤：

S1、对实测的广域电磁法电阻率曲线进行一维有源反演，得到各测点的一维有源反演模型；

S2、根据一维有源反演模型，计算其一维MT响应，得到MT电阻率曲线，即，相当于将广域电磁法电阻率曲线校正到MT电阻率曲线；

S3、对各个校正后的MT电阻率曲线，合成MT剖面，进行MT二维反演，得到二维反演的地电断面图；

S4、利用步骤S3得到的二维反演的地电断面图，进行地质构造的解译。

进一步的，步骤S1具体包括：

S1.1、实测电场迭代计算全域电阻率，所述步骤S1.1具体包括：

S1.1.1、观测点的大地电阻率ρ计算公式如下：

其中r为收发距，即电偶源至观测点的距离，E_mn1为水平电偶源中任意方向的电场分量，I为发射端的供电电流；L为偶极矩；

为观测方位角，是水平电偶源X轴与观测点径向的夹角，即发射端与观测点的夹角，ɑ为发射电极与测量电极的夹角；D＝(1+ikr)e^-ikr，其中k为波数，且k²＝ω²με-iωμσ，ω为发射电流圆频率，μ为大地导磁率，ε为大地介电常数，σ为大地电导率，σ＝1/ρ，ρ为大地电阻率；i为纯虚数；

公式右侧也含有大地电阻率ρ这一未知数，不能一次计算出准确的大地电阻率ρ，需要进行迭代计算才能最终收敛在准确的大地电阻率ρ，保证公式两侧电阻率小于给定的误差限制，迭代过程如下：

S1.1.2、令D＝0，根据公式(1)计算出电阻率初值ρ⁽⁰⁾：

S1.1.3、以ρ⁽⁰⁾为大地电阻率的初值，带入公式(1)的右侧，求得下一次迭代视电阻率ρ⁽¹⁾；

S1.1.4、判断|ρ⁽¹⁾-ρ⁽⁰⁾|/ρ(0)，是否不大于给定误差限制，若是，则停止迭代，得到全域视电阻率；若否，则返回S1.1.3；

最终得到广域电阻率曲线；

S1.2、一维广域电阻率曲线反演拟合

通过括号+下标0表示贝塞尔函数零阶积分，括号+下标1表示贝塞尔函数一阶积分，水平层状地层中任意方向的电场分量E_mn′为：

其中

为观测方位角，是水平电偶源X轴与观测点径向的夹角，即发射端与观测点的夹角；m为空间频率；r为收发距；ɑ为发射电极与测量电极的夹角；h为各层的分层厚度；N为地层的层数；σ各层的为大地电导率；

其中i为纯虚数；ω为发射电流圆频率；μ为大地导磁率；I为发射端的供电电流；L为电偶源极距的长度；k为波数；

R^*和R为空间频率特性函数；

利用公式(3)计算出水平层状地层水平电偶源的电场分量E_mn′，并进一步使用步骤S1.1迭代计算出对应的全域视电阻率，得到层状分层模型的广域电阻率-频率曲线，再通过采用正则化反演方法，从实测的广域电阻率曲线反演得到得到实测广域电阻率-频率曲线所对应的光滑水平电性分层模型，即得到所述各测点的一维有源反演模型。

进一步的，所述步骤S4具体为：

根据步骤S3得到的二维反演地电断面图，对比步骤S1得到的一维有源反演地点断面图，进行地质构造的解译，得到更接近实际情况的地质解译成果。

本发明的有益效果是：

1、本发明能够在广域电阻率曲线一维反演的基础上，得到校正后的MT电阻率曲线，进而利用成熟的MT二维反演方法，得到二维反演结果，极大地减小计算工作量，为广域电磁法的二维反演提出了新的方法；

2、本发明的二维反演结果也可以作为完全有源的广域电磁法二维反演的初始模型，即使今后在有源的广域电磁法二维反演算法上有所突破，采用本发明也可以快速使模型逼近真实的模型，减少耗时较长的有源广域电磁法二维反演次数。

附图说明

图1为本发明一种广域电磁法二维反演方法的流程图；

图2为本发明广域电磁法的坐标系示意图；

图3为有限收发距下，同一水平层状模型的广域电阻率和MT电阻率对比曲线图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

在进行具体实施例的说明之前，先对本发明的背景解释一下：

物探的数据处理的方法，一般是基于均匀半空间的前提下进行的，再拓展到水平层状地层型、二维、三维等复杂模型。

电偶源即发射端，实际工作中，发射端的尺寸相对于测量端与发射端之间的距离很小，可以认为发射端就是电偶源，观测点为测量电极MN中点(如图2所示)。而偶极矩则是大小等于发射电极距AB、方向与发射端电流I方向一致的矢量。本发明中的所有公式符号均通用。

如图1所示，本发明实施例提供一种广域电磁法二维反演方法，包括以下步骤：

S1、对实测的广域电磁法电阻率曲线进行一维有源反演，得到各测点的一维有源反演模型。步骤S1具体包括：

S1.1、实测电场迭代计算全域电阻率

CSAMT是在MT基础上发展起来的，采用人工源提高了信噪比，很快成为一种主要的勘探手段。但由于其电阻率计算公式成立的条件较为苛刻，很多情况下低频段的测量、计算结果，不是地下结构的正确反映，导致CSAMT一直以来存在勘探深度较浅的问题。随着计算技术的发展，目前可以精确计算出单分量的电磁响应，既可以计算出远区的电磁响应，也可以计算出过渡区和近区的电磁响应。在此基础上，何继善院士提出了广域电磁法，即可以在“广大区域”进行观测的电磁法，并获得2018年度国家技术发明一等奖。

广域电磁法野外施工时，除了测量供电电流和MN之间的电压之外，还测量两个发射电极A、B以及两个测量电极M、N的坐标。这样AB与MN夹角α可以通过测量两个发射电极A、B以及两个测量电极M、N的坐标得到。

如图2所示，测量出水平电偶源中任意方向的电场分量E_mn1。

对于均匀半无限空间，水平电偶源在地表处圆柱坐标系r、

z的电场分量E_r和

分别为：

由公式

可知E_mn1应满足下式：

即

其中，D＝(1+ikr)e^-ikr；I为发射端的供电电流；L为偶极矩；ρ为观测点的电阻率；k为波数；r为收发距，即电偶源至观测点的距离；

为观测方位角，是水平电偶源X轴与观测点径向的夹角，即发射端与观测点的夹角；i为纯虚数；z为高度。

电场分量E_mn1可由现有的电场强度测试仪测出。

由现有公式k²＝ω²με-iωμσ，σ＝1/ρ，μ为大地导磁率；ω为发射端的发射电流圆频率可知，D参数中含有电阻率ρ和频率参数ω，实质就是电磁响应的特征。根据公式(3)可知E_mn1与D参数有关，也就是说，观测E_mn1也可以进行频率测深。

对任意方向的电场分量E_mn1进行静态校正后，再通过迭代方法得到对应的广域电阻率，并最终得到广域电阻率-频率曲线。

S1.2、一维广域电阻率曲线反演拟合

实际勘探中，例如对水平层状地层进行勘探，假设水平层状地层在水平电偶源中，地表观测点处的电场分量E_r和

分别为：

其中，

通过用括号+下标0表示贝塞尔函数零阶积分，括号+下标1表示贝塞尔函数一阶积分可得到公式(5)和公式(6)；

则水平层状地层中任意方向的电场分量E_mn′为：

其中，ω为发射端的发射电流圆频率；μ为大地导磁率；m为空间频率；R^*和R为空间频率特性函数；h为各层的分层厚度；N为地层的层数；σ各层的为大地电导率；I为发射端的供电电流；L为电偶源极距的长度；i为纯虚数；k为波数。

已知层状电性分布的情况下，可以通过公式(7)计算出水平层状地层水平电偶源的电场分量E_mn′，再通过公式(4)进行迭代计算，可以计算出模型所对应的广域电阻率-频率曲线。

通过正则化反演方法，得到实测广域电阻率-频率曲线所对应的光滑电性分层模型。正则化反演方法较多，这里不做叙述。

一维有源反演的结果，通过网格化、插值、平滑等数学处理后，可以生成一维反演地电断面图。以上过程就是步骤S1。步骤S1也存在大量的数据预处理和计算工作，正演公式中的贝塞尔函数积分一般采用快速汉克尔变换方法进行计算，反演过程中要计算若干次正演。相对于MT的一维正演，广域电磁法的一维正演工作计算量高出几个数量级。

一维光滑反演结果，可通过网格化、插值、平滑等数学处理后，生成一维反演的断面图，但有时为压制局部畸变曲线的影响、提高深部的分辨率等原因，还是需要进行二维反演进行对比，以便更可靠地解译地下结构。

有源的广域电磁法二维反演方法复杂，详见中南大学何继善、戴世坤、李谛铨提出的有源的广域电磁法二维反演方法。有源场的复杂性决定了其计算工作量非常巨大，仅广域电磁法的一维正演工作计算量就高出MT的一维正演几个数量级。广域电磁法二维反演计算工作可以用海量计算来形容，预计比MT二维反演耗时高出2个数量级。

为类比计算工作量的情况，以一维正演时间为例，同样的地电分层模型，一维广域电磁法正演1000次的时间是11401.57毫秒，而一维MT正演1000次的时间是65.19毫秒，一维广域电磁法正演时间是一维MT正演的175倍。具体耗时长度随着模型的不同、计算机硬件的不同，时间大小上会有差别，但时间比值不会发生变化。反演过程中需要不断地进行正演，比较实测数据与模型响应之间的差别，所以正演时间比值可以认为与反演时间比值是一样的，即一维广域电磁法反演时间是一维MT反演时间的175倍。进一步推广到二维正反演上，可以认为有源的广域电磁法二维反演时间将比无源的MT二维反演时间，高出2个数量级。

为了提高反演效率，快速得到二维反演结果，本发明提出了将测量的广域电阻率校正为MT电阻率的方法，具体见步骤S2。

S2、根据一维有源反演模型，计算其一维MT响应，得到MT的电阻率和相位曲线，即，相当于将广域电磁法电阻率曲线校正到MT电阻率曲线。

理论分析可知，有限收发距条件下，广域电阻率曲线和MT电阻率曲线相差很大，用MT方法反演广域电阻率曲线将得到误差较大的结果。S1已经得到了实测广域电阻率曲线所对应的一维地电模型，采用该地电模型，可以计算出对应的MT电阻率曲线，这一过程相当于将广域电磁法电阻率曲线校正到MT电阻率曲线。

有限收发距下，广域电阻率曲线和MT电阻率曲线形态差异较大，见图3所示，在8000m的收发距下，某水平地层的广域电阻率曲线比MT电阻率曲线复杂得多。本步骤的实施，就是把广域电磁法电阻率曲线校正到MT电阻率曲线。

图3中计算条件如下：

收发距：8000m

给定一个5层电性的模型如下：

分层电阻率(欧姆米)：	150	350	50	500	3000
						分层厚度(m)：	80	160	320	600	∞

S3、对各个校正后的MT电阻率曲线，合成MT剖面，进行MT二维反演，得到二维反演的地电断面图。

得到各测点的MT电阻率曲线(相位曲线)后，可以合成MT测量剖面，并进行MT二维反演。由于MT正演计算相对简单，其计算量远远小于广域电磁法的正演，极大地简化了程序设计，缩短了运行时间。

关于MT二维反演的方法很多，一般采用正则化反演方法，这里不对MT二维反演方法进行叙述。

S4、得到二维反演的地电断面图，用于地质构造的解译。

通过二维反演得到二维反演地电断面图，和步骤S1得到的一维有源反演地电断面图，进行比对，共同分析、解译，得到更接近实际情况的地质解译成果。

实际工作中发现，一维反演地电断面图有时能反映更多的异常细节，而二维反演结果有时由于固有的二维模型光滑度控制等原因，压制了较多的异常细节，因此将一维和二维结果进行对比分析，再结合地质资料进行分析，可能发现更多的物探异常，有利于物探资料的地质解译。

本发明在一维反演部分用有源场进行反演，二维反演部分用无源场进行反演，两者之间的衔接就是一维模型的不同的正反演，通过一维模型的过渡，将复杂的广域电磁法(WFEM)电阻率曲线校正成MT(大地电磁测深法)电阻率曲线，进而减少广域电磁法二维反演计算工作量。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何属于本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (3)

1.一种广域电磁法二维反演方法，其特征在于，包括以下步骤：

S1、对实测的广域电磁法电阻率曲线进行一维有源反演，得到各测点的一维有源反演模型；

S2、根据一维有源反演模型，计算其一维MT响应，得到MT电阻率曲线，即，相当于将广域电磁法电阻率曲线校正到MT电阻率曲线；

S3、对各个校正后的MT电阻率曲线，合成MT剖面，进行MT二维反演，得到二维反演的地电断面图；

S4、利用步骤S3得到的二维反演的地电断面图，进行地质构造的解译。

2.根据权利要求1所述的一种广域电磁法二维反演方法，其特征在于，步骤S1具体包括：

S1.1、实测电场迭代计算全域电阻率，所述步骤S1.1具体包括：

S1.1.1、观测点的大地电阻率ρ计算公式如下：

其中r为收发距，即电偶源至观测点的距离，E_mn1为水平电偶源中任意方向的电场分量，I为发射端的供电电流；L为偶极矩；

为观测方位角，是水平电偶源X轴与观测点径向的夹角，即发射端与观测点的夹角，α为发射电极与测量电极的夹角；D＝(1+ikr)e^-ikr，其中k为波数，且k²＝ω²με-iωμσ，ω为发射电流圆频率，μ为大地导磁率，ε为大地介电常数，σ为大地电导率，σ＝1/ρ，ρ为大地电阻率；i为纯虚数；

公式右侧也含有大地电阻率ρ这一未知数，不能一次计算出准确的大地电阻率ρ，需要进行迭代计算才能最终收敛在准确的大地电阻率ρ，保证公式两侧电阻率小于给定的误差限制，迭代过程如下：

S1.1.2、令D＝0，根据公式(1)计算出电阻率初值ρ⁽⁰⁾：

S1.1.3、以ρ⁽⁰⁾为大地电阻率的初值，带入公式(1)的右侧，求得下一次迭代视电阻率ρ⁽¹⁾；

S1.1.4、判断|ρ⁽¹⁾-ρ⁽⁰⁾|/ρ⁽⁰⁾，是否不大于给定误差限制，若是，则停止迭代，得到全域视电阻率；若否，则返回S1.1.3；

最终得到广域电阻率曲线；

S1.2、一维广域电阻率曲线反演拟合

通过括号+下标0表示贝塞尔函数零阶积分，括号+下标1表示贝塞尔函数一阶积分，水平层状地层中任意方向的电场分量E_mn′为：

其中

为观测方位角，是水平电偶源X轴与观测点径向的夹角，即发射端与观测点的夹角；m为空间频率；r为收发距；α为发射电极与测量电极的夹角；h为各层的分层厚度；N为地层的层数；σ各层的为大地电导率；

其中i为纯虚数；ω为发射电流圆频率；μ为大地导磁率；I为发射端的供电电流；L为电偶源极距的长度；k为波数；

R^*和R为空间频率特性函数；

利用公式(3)计算出水平层状地层水平电偶源的电场分量E_mn′，并进一步使用步骤S1.1迭代计算出对应的全域视电阻率，得到层状分层模型的广域电阻率-频率曲线，再通过采用正则化反演方法，从实测的广域电阻率曲线反演得到得到实测广域电阻率-频率曲线所对应的光滑水平电性分层模型，即得到所述各测点的一维有源反演模型。

3.根据权利要求1所述的一种广域电磁法二维反演方法，其特征在于，所述步骤S4具体为：

根据步骤S3得到的二维反演地电断面图，对比步骤S1得到的一维有源反演地点断面图，进行地质构造的解译，得到更接近实际情况的地质解译成果。

Images