基于发射电流全波形的地面瞬变电磁法反演方法及装置

技术领域

本发明涉及地球物理探测领域，具体涉及一种基于发射电流全波形的地面瞬变电磁法反演方法及装置。

背景技术

地面瞬变电磁法多采用“双极性”发射电流波形，即一个发射周期内既正向供电，也负向供电。如图1所示，正向和负向供电半周期均包含电流开启、稳定和关断阶段，以及测量阶段。“双极性”波形中，发射电流每一次开启和关断都会在地下介质中感应产生涡流，并且发射电流前一次开启或关断产生的涡流必然会部分抵消后一次关断或开启产生的涡流。受上述现象影响，正向和负向供电两个半周期内分别测量的电磁响应幅值存在差异。为了确保连续两个供电半周期内测量的电磁响应绝对值差异在极小范围之内，瞬变电磁法通常在多个发射周期之后才开始记录数据。此外，为了降低随机干扰影响，瞬变电磁法数据测量通常采用多次叠加技术，即需要连续观测叠加次数规定发射周期数目的电磁信号。因此，瞬变电磁法观测数据必然受发射电流波形，及其周期数目的影响。

目前，地面瞬变电磁法数据处理软件采用斜阶跃电流波形，即只考虑关断时间影响，而忽略了发射电流波形中稳定时间和开启时间，以及周期数目影响。然而，在高电导率勘探背景下，伴随电流开启产生的涡流信号强、衰减慢，即使稳定时间长达数十毫秒也未必充分衰减，依然会部分抵消关断后瞬变电磁响应，且在晚期尤为明显。因此，只考虑关断时间影响的地面瞬变电磁法数据处理方法不适用于晚期数据，并制约了进一步提高数据处理解释精度。

发明内容

本发明解决的主要技术问题在于，在考虑关断时间影响的基础上，还考虑发射电流波形中稳定时间和开启时间，以及周期数目影响，以提高地面瞬变电磁法数据一维反演解释精度。

为了实现上述目的，根据本发明的一个方面，提供了一种基于发射电流全波形的地面瞬变电磁法反演方法，包括以下步骤：

步骤一：基于Tikhonov正则化方式构建一维反演目标函数，所述反演目标函数

式中，||W

步骤二：采用高斯牛顿法对反演目标函数进行最优化，输出最终反演模型；

所述采用高斯牛顿法对反演目标函数进行最优化的过程，包括：

步骤S1：输入瞬变电磁法的观测数据和工作参数；

步骤S2：输入初始反演模型、参考模型和控制参数；

步骤S3：根据所述工作参数，计算所述初始反演模型的瞬变电磁法正演响应；

步骤S3包括：

S3.1：计算下阶跃电流波形激发的磁感应强度脉冲响应或阶跃响应；

S3.2：根据所述下阶跃电流波形激发的磁感应强度脉冲响应或阶跃响应，采用卷积算法计算发射电流基本波形激发的磁感应强度脉冲响应或阶跃响应；

所述发射电流基本波形分为两类：第一类只包含关断时间，即斜阶跃波形；第二类包括多个周期数目不同且由开启时间、稳定时间、关断时间和占空比确定的全波形，全波形参数取决于实际电流波形，由步骤S1输入；

步骤S4：采用扰动法计算当前反演模型(即第k次反演迭代获得的反演模型)的雅克比矩阵；

步骤S5：根据所述观测数据、所述参考模型、所述控制参数和所述雅克比矩阵生成高斯牛顿法正规方程，并求解所述高斯牛顿法正规方程，获得模型更新量；

步骤S6：根据所述模型更新量更新所述当前反演模型，并计算更新后反演模型的正演响应和均方根拟合差；

步骤S7：判断当前迭代次数或所述均方根拟合差是否满足所述控制参数中的迭代终止条件，若是，则停止迭代并输出最终反演模型，否则，返回步骤S4开启下一次迭代。

优选地，步骤S1中，所述观测数据包括：观测时间序列对应的按接收探头等效面积和发射电流强度归一化的磁感应强度脉冲响应或阶跃响应和相应的测量标准差；

所述工作参数包括：工作装置、发射回线尺度和位置、发射电流波形参数、测点位置和观测时间序列；

其中，所述发射电流波形参数包括：基频、占空比、开启时间、稳定时间和关断时间。

优选地，步骤S2中，所述初始反演模型为一维层状介质，所述输入初始反演模型包括：输入层状介质层数、各层介质厚度和电导率；

通过输入参考模型的形式输入先验信息；

所述控制参数包括：目标均方根拟合差、最大迭代次数、正则化因子初值及更新参数。

优选地，步骤S31包括：

铺设于水平层状介质表面的矩形发射回线源通以一定频率f的电流，在测点(x,y,0)激发的频率域磁感应强度沿z轴方向的分量B

其中，

式中，E

采用12阶高斯—勒让德积分和Hankel变换分别求解公式三至公式六中的外层和内层积分，并获得频率域磁感应强度；

对所述频率域磁感应强度的虚部分别开展正弦或余弦变换，获得下阶跃电流波形激发的磁感应强度脉冲响应

优选地，步骤S32包括：

将计算得到的所述下阶跃电流波形激发的磁感应强度脉冲响应或阶跃响应与任一基本波形中电流对时间的微分做卷积运算，这一过程表示为公式九：

式中，f

采用梯形积分方法求解公式九，获得该基本波形激发的磁感应强度脉冲响应或阶跃响应。

优选地，步骤S4中，扰动量根据当前反演模型m

优选地，步骤S5中，对于反演目标函数，高斯牛顿法第k+1次迭代的正规方程为公式十：

式中，δm

采用预条件共轭梯度法求解公式十，获得δm

优选地，步骤S6包括：

根据模型更新量采用公式十一对当前反演模型进行更新：

m

式中，a为模型更新步长，δm

更新模型后，再采用步骤S3中的方法计算该更新后的反演模型的瞬变电磁法正演响应，并根据公式十二计算均方根拟合差x

式中，M为数据数目，s

优选地，步骤S7包括：

判断所述均方根拟合差是否满足小于或等于目标拟合差，或达到设定迭代次数，若满足，则停止迭代并输出最终反演模型；否则，开启下一次迭代。

根据本发明的另一方面，提供了一种基于发射电流全波形的地面瞬变电磁法反演装置，包括：

反演目标函数构建模块，用于基于Tikhonov正则化方式构建一维反演目标函数，所述反演目标函数

式中，||W

反演目标函数优化模块，用于采用高斯牛顿法对反演目标函数进行最优化，输出最终反演模型；

所述反演目标函数优化模块包括：

第一输入单元，用于输入瞬变电磁法的观测数据和工作参数；

第二输入单元，用于输入初始反演模型、参考模型和控制参数；

正演响应计算单元，用于根据所述工作参数，计算所述初始反演模型的瞬变电磁法正演响应；

其中，所述正演响应计算单元包括：

第一计算单元，用于计算下阶跃电流波形激发的磁感应强度脉冲响应或阶跃响应；

第二计算单元，用于根据所述下阶跃电流波形激发的磁感应强度脉冲响应或阶跃响应，采用卷积算法计算发射电流基本波形激发的磁感应强度脉冲响应或阶跃响应；

所述发射电流基本波形分为两类：第一类只包含关断时间，即斜阶跃波形；第二类包括多个周期数目不同且由开启时间、稳定时间、关断时间和占空比确定的全波形，全波形参数取决于实际电流波形，由所述第一输入单元输入；

雅克比矩阵计算单元，用于采用扰动法计算当前反演模型的雅克比矩阵；

正规方程生成及求解单元，用于根据所述观测数据、所述参考模型、所述控制参数和所述雅克比矩阵生成高斯牛顿法正规方程，并求解所述高斯牛顿法正规方程，获得模型更新量；

迭代更新单元，用于根据所述模型更新量更新所述当前反演模型，并计算更新后反演模型的正演响应和均方根拟合差；

迭代终止单元，用于判断当前迭代次数或所述均方根拟合差是否满足所述控制参数中的迭代终止条件，若是，则停止迭代并输出最终反演模型，否则，返回所述雅克比矩阵计算单元开启下一次迭代。

本发明提供的技术方案具有以下有益效果：

提供了一种基于发射电流全波形的地面瞬变电磁法反演方法及装置，将用于地面瞬变电磁法数据处理的发射电流基本波形由斜阶跃波形，拓展为包含开启时间、稳定时间、关断时间和占空比的发射电流全波形，并在发射电流基本波形中纳入了周期数目参数。这一策略不仅有助于提高数据处理解释精度，也有助于完善仪器参数设计。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是单个周期“双极性”波形示意图；

图2是实施例1和实施例2中一种基于发射电流全波形的地面瞬变电磁法反演方法的流程图；

图3是实施例1和实施例2中发射电流基本波形示意图，其中图3(a)波形1(现有方法)，图3(b)波形2(本发明)，图3(c)波形3(本发明)，图3(d)波形4(本发明)；

图4是实施例1和实施例2中层状介质模型示意图；

图5是实施例1中基于图3的4种基本波形一维反演获取的地电模型，其中图5(a)波形1，图5(b)波形2，图5(c)波形3，图5(d)波形4；

图6是实施例2中野外观测数据(磁感应强度阶跃响应)曲线；

图7是对于图6的观测数据，基于图3的4种基本波形一维反演获取的地电模型；

图8是本发明提供的一种基于发射电流全波形的地面瞬变电磁法反演装置的结构图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

实施例1：理论模型

基于Tikhonov正则化方式构建一维反演目标函数，采用高斯牛顿法对其最优化。反演目标函数

式中，第一项为观测数据与正演模型响应的拟合差函数、第二项为模型正则化函数，||·||表示l

如图3所示，本实施例的反演中，发射电流基本波形先后采用了4种基本波形，并对每种波形都开展了10次反演计算。

如图2所示，每次反演计算的过程包括以下步骤：

步骤S1：输入瞬变电磁法的观测数据和工作参数；

采用理论模型时，观测数据包括两部分：即特定地电模型和工作参数的正演响应(由步骤S3中方法生成)与随机噪音。

地电模型和工作参数如下：地电模型为5Ωm均匀半空间；发射源为边长100m的方形回线，铺设于均匀半空间表面；装置类型为中心回线装置；发射电流为1A；发射电流基本波形为加拿大PROTEM仪器采用的、发射基频25Hz的双极性电流波形(占空比50％)，相应时间序列包含30个时刻，范围为6.8–6978μs。其中，25Hz发射基频相应的周期为0.04s，即0.25周期为10ms，并设置开启时间、稳定时间和关断时间依次为1ms、8.99ms和0.01ms。正演响应由8.25个周期的基本波形激发产生的感应电动势(磁感应强度脉冲响应)。

随机噪音按下述方法生成：首先生成高斯噪音，再乘以标准差作为随机噪音。标准差设置为某一时刻正演响应

将正演响应与随机噪音相加作为观测数据，用于反演。

表1理论模型正演响应标准差设置表

步骤S2：输入初始反演模型、参考模型和控制参数；

初始反演模型为39层层状介质，每层介质电阻率为20Ωm，厚度由5m阶梯式递增至40m(第1至10层为5m、第11至25层为10m、第26至35层为20m、第36至38层为40m、第39层介质为均匀半空间)，层厚度在反演中固定不变。无参考模型。目标均方根拟合差为1，最大迭代次数为30。正则化因子初值α

步骤S3：按照步骤S1输入的瞬变电磁法工作参数，采用步骤S3所述的方法计算基本波形激励下的磁感应强度脉冲响应

步骤S3这一步骤是区别于其他反演方法的关键，即在正演响应计算中纳入了发射电流全波形影响。步骤S3可细分两步，首先计算下阶跃电流波形激发的磁感应强度脉冲响应

铺设于水平层状介质表面的矩形发射回线源通以一定频率f的电流，在测点(x,y,0)激发的频率域磁感应强度z分量B

其中，

式中，I表示电流强度，W和L表示矩形发射回线的半长和半宽。z

r

采用12阶高斯—勒让德积分和Hankel变换分别求解公式(4)至公式(7)中的外层和内层积分，并获得频率域磁感应强度。

对上述频率域磁感应强度虚部分别开展正弦变换，获得下阶跃电流波形激发的磁感应强度脉冲响应

再采用卷积算法求解发射电流基本波形激发的电磁场。卷积算法中，首先根据式(3)计算下阶跃电流波形激发的磁感应强度脉冲响应(f

采用梯形积分方法求解式(8)。

步骤S4：采用扰动法计算当前反演模型的雅克比矩阵。本实施例中为39层介质，由于层厚度固定不变，需要对每层介质电阻率逐一扰动，扰动39次后即可获得完整的雅克比矩阵。扰动量为当前反演模型某一参数的10％。在一些实施例中，层状介质的层数可根据实际需要进行选择，对应的扰动次数及扰动量也可相应进行调整。

步骤S5：代入步骤S3和步骤S4分别生成正演响应和雅克比矩阵，再代入当前正则化因子，生成高斯牛顿法正规方程。采用预条件共轭梯度法求解正规方程，其终止条件为迭代次数达到200次或相对残差小于等于10

步骤S6：计算模型更新量和均方根拟合差，其中模型更新步长设置为1。

步骤S7：进行终止条件判断。这一理论模型研究中，均方根拟合差均在10次反演迭代之内小于目标值。

上述反演流程运行了10次。用于反演的观测数据含有随机噪音，因此每次反演获取的地电模型都略有差异。反演结果如图5所示，其中实线表示理论模型，虚线表示反演模型，采用波形1(现有方法)时，绝大多数反演地电模型与理论模型差异明显(见图5(a))；采用波形2、3和4(本发明)时，绝大多数反演地电模型与理论模型吻合(见图5(b)、5(c)和5(d))。因此，与现有基于斜阶跃电流波形的数据处理方法相比，本发明提出的基于发射电流全波形的地面瞬变电磁法反演方法具有明显优势，实现了预期目标。

实施例2：野外观测数据

步骤S1：输入瞬变电磁法的观测数据和工作参数；

观测数据采集地点为中国东海某滩涂，该地区处于海陆交替环境，浅部地层一般饱含海水，电阻率仅为0.1至0.8Ωm，深部地层电阻率略有增大。观测数据为一个测点的磁感应强度阶跃响应，由加拿大Crone仪器公司生产的超导量子干涉仪(SQUID)采集。发射源为边长400m×200m的矩形回线，回线中心坐标为(0,0,0)m，测点坐标为(-100,0,0)m。发射电流为25A；发射电流基本波形为加拿大Crone Pulse EM仪器采用的、发射基频0.833333Hz的双极性电流波形(占空比50％)，相应时间序列包含45个时刻，范围为-0.15ms–222.7ms。其中，0.833333Hz发射基频相应的周期为1.2s，即0.25周期为300ms，并设置开启时间、稳定时间和关断时间依次为0.5ms、299.5ms和1ms。反演中，只使用了第2至45时刻的观测数据，见图6。由于Crone Pulse EM仪器采集的观测数据不含标准差，本发明依据瞬变电磁法早期信号强、信噪比高，晚期信号弱、信噪比低的特点，为实测数据匹配了标准差，详见表2。

表2野外观测数据标准差设置表

步骤S2：输入初始反演模型、参考模型和控制参数；

初始反演模型为46层层状介质，每层介质电阻率为1Ωm，厚度由10m阶梯式递增至40m(第1至20层为10m、第21至30层为20m、第31至45层为40m、第46层介质为均匀半空间)，层厚度在反演中固定不变。无参考模型。目标均方根拟合差为1，最大迭代次数为30。正则化因子初值α

步骤S3：按照步骤S1输入的瞬变电磁法工作参数，采用步骤S3所述的方法计算基本波形激励下的磁感应强度阶跃响应b

步骤S4：采用扰动法计算当前反演模型的雅克比矩阵。本实施例中为46层介质，由于层厚度固定不变，需要对每层介质电阻率逐一扰动，扰动46次后即可获得完整的雅克比矩阵。

步骤S5：代入步骤S3和步骤S4分别生成的正演响应和雅克比矩阵，再代入当前正则化因子，生成高斯牛顿法正规方程。采用预条件共轭梯度法求解正规方程，其终止条件为迭代次数达到200次或相对残差小于等于10

步骤S6：计算模型更新量和均方根拟合差，其中模型更新步长设置为1。

步骤S7：进行终止条件判断。这一理论模型研究中，均方根拟合差均在10次反演迭代之内小于目标值。

反演结果如图7所示，采用波形1(现有方法)获取的反演结果与采用波形2、3和4(本发明)获取的反演结果在地下100至1000m存在明显差异，且采用波形2、3和4(本发明)获取的反演结果更加符合真实的地质情况。因此，本发明提出的基于发射电流全波形的地面瞬变电磁法反演方法实现了预期目标。

如图8所示，本实施例还提供了一种基于发射电流全波形的地面瞬变电磁法反演装置，包括：

反演目标函数构建模块01，用于基于Tikhonov正则化方式构建一维反演目标函数，反演目标函数

式中，||W

反演目标函数优化模块02，用于采用高斯牛顿法对反演目标函数进行最优化，输出最终反演模型；

反演目标函数优化模块02包括：

第一输入单元021，用于输入瞬变电磁法的观测数据和工作参数；

第二输入单元022，用于输入初始反演模型、参考模型和控制参数；

正演响应计算单元023，用于根据输入的工作参数，计算初始反演模型的瞬变电磁法正演响应；

正演响应计算单元023包括：

第一计算单元0231，用于计算下阶跃电流波形激发的磁感应强度脉冲响应或阶跃响应；

第二计算单元0232，用于根据下阶跃电流波形激发的磁感应强度脉冲响应或阶跃响应，采用卷积算法计算发射电流基本波形激发的磁感应强度脉冲响应或阶跃响应；

发射电流基本波形分为两类：第一类只包含关断时间，即斜阶跃波形；第二类包括多个周期数目不同且由开启时间、稳定时间、关断时间和占空比确定的全波形，全波形参数取决于实际电流波形，由第一输入单元021输入；

雅克比矩阵计算单元024，用于采用扰动法计算当前反演模型的雅克比矩阵；

正规方程生成及求解单元025，用于根据输入的观测数据、参考模型、控制参数和计算得到的雅克比矩阵生成高斯牛顿法正规方程，并求解所述高斯牛顿法正规方程，获得模型更新量；

迭代更新单元026，用于根据所述模型更新量更新所述当前反演模型，并计算更新后反演模型的正演响应和均方根拟合差；

迭代终止单元027，用于判断当前迭代次数或所述均方根拟合差是否满足所述控制参数中的迭代终止条件，若是，则停止迭代并输出最终反演模型，否则，返回雅克比矩阵计算单元024开启下一次迭代。

需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者系统不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者系统所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者系统中还存在另外的相同要素。

上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。词语第一、第二、以及第三等的使用不表示任何顺序，可将这些词语解释为标识。

以上仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是采用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。