地面回线源瞬变电磁法一维反演系统及其应用
A 1D inversion system of the ground-based loop-source transient electromagnetic method
通讯作者: 李建慧(1982-), 男, 教授, 博士生导师, 主要研究方向为电磁法数值计算与资料处理研究。 Email:ljhiiicumt@126.com
责任编辑: 沈效群
收稿日期: 2021-12-24 修回日期: 2022-04-28
基金资助: |
|
Received: 2021-12-24 Revised: 2022-04-28
作者简介 About authors
张文波(1975-), 男, 博士, 讲师, 长期从事电磁法勘探工作。 Email:
瞬变电磁法现阶段的资料处理解释仍以一维反演为主,因此开发一套功能齐全、高效稳定的一维反演系统对进一步提高国内瞬变电磁法的应用水平仍具有重要意义。本研究开发的地面回线源瞬变电磁法一维反演系统包括了基于高斯牛顿法的最小构造反演和Occam反演,也包括了基于阻尼最小二乘法的横向约束反演和空间约束反演。以内蒙古那仁宝力格煤田玄武岩岩体形态探测为例,将该反演系统的最小构造反演和Occam反演结果与商业软件IX1D进行了对比验证,发现不同反演方法获取的电阻率二维断面图中,玄武岩岩体形态相似、电阻率范围一致;结合钻孔资料,这些一维反演结果清晰反映了玄武岩岩体除岩浆上涌通道区域外的分布形态。将横向约束反演和空间约束反演应用于该实例,结果表明:相邻测点间玄武岩岩体电阻率差异缩小,玄武岩与沉积岩界面的连续性得到增强。
关键词:
The processing and interpretation of the data derived using the transient electromagnetic (TEM) method are still mainly conducted through one-dimensional (1D) inversion presently. Therefore, developing an efficient and s
Keywords:
本文引用格式
张文波, 张莹, 李建慧.
ZHANG Wen-Bo, ZHANG Ying, LI Jian-Hui.
0 引言
电磁法一维反演的常用方法有最小构造反演、Occam反演、横向约束反演(laterally constrained inversion,LCI)、空间约束反演(spatially constrained inversion,SCI),以及模拟退火法、深度学习等多种人工智能优化算法。一维最小构造反演是求取多层介质模型的最光滑解,即在一定拟合差条件下,使得模型粗糙度最小,因此,目标函数不仅包含数据拟合项,也包含模型正则化项[3-4]。Occam反演与最小构造反演类似,不同之处在于正则化因子的选取:Occam反演通过一元函数优化来选取正则化因子,而最小构造反演中正则化因子通常固定不变或者通过一些非常简单的算法选取[5],因此,Occam反演更加稳定,对正则化初始值依赖程度低,但其计算量偏大。LCI方法同时反演电阻率和层厚度,并在目标函数中加入相邻测点之间电阻率、层厚度和界面深度约束项[6],因此单条测线所有测点数据需整体反演。SCI与LCI类似,不仅考虑同一条测线相邻测点之间的约束,还需考虑其他测线相邻测点之间的约束[7],因此测区内所有测点数据需整体反演。
目前,国内外多个课题组均已开发了瞬变电磁法一维反演软件。如:加拿大英属哥伦比亚大学Oldenburg课题组在国际上率先实现了用于回线源瞬变电磁法数据处理的一维最小构造反演[8],开发了实用化的软件包EM1DTM,并成功应用于多个野外实例[9-10];德国科隆大学Tezkan课题组采用阻尼最小二乘法和Occam算法开展了回线源瞬变电磁法、长偏移距瞬变电磁法、射频大地电磁法等方法的单一反演或2种以上方法联合反演[11⇓⇓-14],并开发了EMUPLUS软件包;丹麦奥胡斯大学Auken课题组基于阻尼最小二乘法开发了瞬变电磁法一维LCI程序,用于解决反演电阻率和层界面横向连续性差等问题[15],还开发了SCI程序,用于解决反演电阻率和层界面水平方向连续性差的问题[7];Kirkegaard和Auken在并行计算、迭代法求解方程组、灵敏度矩阵近似等方面优化了算法,进而提升了反演效率[16]。上述反演算法和策略均集成于AarhusInv软件之中。此外,IX1D、Maxwell等商业软件也都包含瞬变电磁法一维反演模块。
国内多个课题组也开展了相关研究。吉林大学殷长春课题组先后开发了适用于航空瞬变电磁法资料处理解释的一维LCI和SCI程序[17-18],采用Occam算法对m序列发射波形多道瞬变电磁法数据开展了反演研究[19],还采用深度学习方法实现了航空瞬变电磁法的一维反演快速成像[20];北京大学黄清华课题组开展了地面回线源瞬变电磁法的一维最小构造反演研究,该算法适用于水平方向发生形变的复杂发射回线[21];中国矿业大学(北京)程久龙课题组开展了粒子群优化——阻尼最小二乘混合算法研究,并应用于矿井瞬变电磁法探测[22];山东大学孙怀凤等采用模拟退火法开展了回线源瞬变电磁法一维反演,并取得了良好效果[23]。
1 反演算法原理
图1
1.1 最小构造反演与Occam反演
基于Tikhonov正则化方式构建反演目标函数φ,有:
式中:右边第一项为实测数据与正演模型响应的拟合差函数,第二项为模型正则化函数,‖·‖表示l2范数;dobs为实测数据向量;dpre为模型参数向量m相应的正演模型响应向量,由瞬变电磁法一维正演计算;Wd为数据加权矩阵,此处为一个对角矩阵,其元素为实测数据标准差的倒数;Wm为模型粗糙度算子;λ为正则化因子;m*为参考模型,先验信息包含其中。这两种反演算法中,m为各层介质电导率。
对于目标函数式(1),采用高斯牛顿法[28]最优化时,第k+1次迭代的正规方程为:
图2
1.2 LCI和SCI
LCI和SCI中,构建如下目标函数φ:
式中:右边第一项为实测数据与正演模型响应的拟合差函数;第二项为模型电导率和层厚约束项,Rp和α分别为其约束矩阵和权重;第三项为模型层界面深度约束项,Rh和β分别为其约束矩阵和权重。其中,Rp和Rh具体形式可参考文献[6],其他参数与上文一致。
对于目标函数式(3),采用阻尼最小二乘法[28]最优化时,第k+1次迭代的正规方程为:
式中:I为单位矩阵,λ为阻尼因子,仍采用“冷却法”更新。图3为LCI和SCI流程图,合成正规方程后,其流程与最小构造反演一致。
图3
LCI和SCI的区别在于参与约束的测点不同。如图4所示,LCI中对于某一测点(红色),参与约束的为该测点所在测线中的2个相邻点(绿色),而SCI中对于某一测点(红色),参与约束的为该测点的所有邻近点(绿色)。本研究中,以某一测点为圆心,按特定搜索半径自动寻找各个方向与该测点距离最近的测点,并将这些测点作为约束测点。
图4
图4
LCI和SCI参与约束测点示意
Fig.4
The diagram showing the observation points used as constrained points in the LCI and SCI
2 应用实例
2.1 实例概述
内蒙古自治区那仁宝力格煤田勘探区内地势平坦,新近系火山喷出岩比较发育,火山岩浆喷溢时的通道穿越煤层时不但对煤层产生破坏,而且对煤的变质也会有不同程度的影响。为了保证今后煤田的顺利开采,瞬变电磁法被用于确定岩浆上涌通道,以及通道周边玄武岩与围岩的界面。
从火山口和火山锥的存在和分布形态看,岩浆以中心式喷发为主,裂隙式喷发次之;岩性组合主要有橄榄拉斑玄武岩和玄武岩。勘探区内,玄武岩体电阻率可达上千欧姆米,而砂岩和泥岩地层电阻率仅为几十欧姆米,二者电性差异明显[30]。这也是采用瞬变电磁法圈定玄武岩岩体在沉积岩中分布范围的前提。
如图5所示,野外工作采用了大回线源装置,发射回线边长为600 m×300 m,其中600 m边长沿测线方向布置。本次野外工作布设了11条测线,线距为100 m,点距为40 m,共计486个测点。160测线为控制测线,穿越了整个玄武岩露头,其他测线均位于玄武岩露头中心区域。对于160测线,铺设1次回线源只采集6个测点数据;对于其他测线,铺设1次回线源可采集2条相邻测线的12个测点数据。采用加拿大Phoenix公司生产的V8多功能电法仪,观测时间序列为分布于0.5~20 ms之间的30个时刻。
图5
图5
那仁宝力格煤田测点分布示意
Fig.5
The distribution diagram of measuring points in the Narenbaolige coalfield
2.2 单点反演
为了验证算法的正确性,首先将本研究开发的最小构造反演和Occam反演程序与商业软件IX1D对比。对于最小构造反演和Occam反演,初始模型为36层层状介质,每层电阻率均为101.5 Ω·m,目标拟合差为3。IX1D软件使用上述初始模型时,在多个测点无法拟合实测数据,因此使用了软件自动估算的初始模型。由于IX1D软件无法设置目标拟合差,其在多数测点的最终均方根拟合差(RMS)大于3(图6所示);而本研究开发的最小构造反演收敛性良好,所有测点均方根拟合差均在3以内。
图6
图6
IX1D软件和本研究最小构造反演结果的RMS对比
Fig.6
The comparison of RMS for IX1D software and the minimum-structure inversion method presented here
图7为160线由单点反演地电模型拼接而成的电阻率二维断面。图中,高阻异常体(玄武岩岩体)形态对称,呈“锅底状”分布,并且3种方法的反演电阻率分布范围基本一致。其中,IX1D反演结果在相邻测点之间跳跃性较强,玄武岩岩体与沉积岩界面连续性较差。这一现象与部分测点最终均方根拟合差偏大有关。此外,3种反演结果与钻孔揭露的玄武岩厚度吻合良好,仅在钻孔ZK3处差异明显。该钻孔终孔深度为716.6 m,未揭露下伏沉积岩,结合该钻孔位于高阻异常体中心,推断其极有可能位于火山通道分布范围之内。因此,一维反演能够获得玄武岩岩体除岩浆上涌通道区域外的分布形态,却无法直接反映岩浆上涌通道信息。这一现象与电磁场的传播规律有关,即随着探测深度的增大,电磁成像分辨率随之下降,也与玄武岩岩体“倒锥形”的真实形态有关。
图7
图7
160线单点反演电阻率断面
Fig.7
The section view stitched from the single-point 1D inversion models for survey line 160
图8
图8
测区最小构造反演电阻率平面(z=0 m)
Fig.8
The plane view of the 1D minimum-structure inversion results (z=0 m)
2.3 LCI反演
根据最小构造反演和Occam反演结果,研究区域电性结构相对简单,界面连续性良好,适合开展LCI反演。反演中,初始模型设置为三层层状介质,每层电阻率均为101.5 Ω·m,第一、第二层的厚度均为100 m;约束权重α和β依次设置为10、100和1 000;目标拟合差仍为3。在移动工作站(处理器为4核8线程、2.7 GHz主频)、Ubuntu操作系统条件下,单条测线一次反演耗时约为10 min,最终拟合差均小于3。
图9
图9
160线LCI反演电阻率断面
Fig.9
The section view stitched from LCI inversion models for survey line 160
图10为测区LCI反演结果俯视图(α=β=1 000),图中同一测线相邻测点之间的电阻率和层厚连续性良好,但测线之间相邻测点的电阻率差异依然较大。
图10
图10
测区LCI反演结果俯视图(z=0 m,α=β=1 000)
Fig.10
The plane view of the LCI results (z=0 m, α=β=1 000)
2.4 SCI反演
图11
图11
160线SCI反演电阻率断面
Fig.11
The section view stitched from SCI inversion models for survey line 160
图12
图12
测区SCI反演结果俯视图(z=0 m,α=β=10)
Fig.12
The plane view of the SCI results (z=0 m, α=β=10)
图13
图13
测区SCI反演结果俯视图(z=0 m,α=β=100)
Fig.13
The plane view of the SCI results (z=0 m,α=β=100)
图14
图14
测区SCI反演结果俯视图(z=0 m,α=β=1 000)
Fig.14
The plane view of the SCI results (z=0 m, α=β=1 000)
3 结论
本文详细叙述了笔者团队研发的地面回线源瞬变电磁法一维反演系统,并将其应用于内蒙古那仁宝力格煤田玄武岩岩体形态探测实例。该实例中,单点反演、LCI和SCI刻画了玄武岩岩体除岩浆上涌通道区域之外的分布形态。下一步,将以上述一维反演结果为初始模型,开展该实例的三维反演研究,力争获得更加精确的玄武岩岩体形态。
参考文献
An overview of a highly versatile forward and stable inverse algorithm for airborne, ground-based and borehole electromagnetic and electric data
[J].DOI:10.1071/EG13097 URL [本文引用: 1]
地球物理三维电磁反演方法研究动态
[J].
Status and prospect of 3D inversions in EM geophysics
[J].
Magnetotelluric inversion for minimum structure
[J].DOI:10.1190/1.1442438 URL [本文引用: 1]
基于遗传算法的CSAMT最小构造反演
[J].
The application of genetic algorithm to CSAMT inversion for minimum structure
[J].
Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data
[J].DOI:10.1190/1.1442303 URL [本文引用: 2]
Layered and laterally constrained 2D inversion of resistivity data
[J].DOI:10.1190/1.1759461 URL [本文引用: 2]
Quasi-3D modeling of airborne TEM data by spatially constrained inversion
[J].
Inversion of time-domain electromagnetic data for a horizontally layered earth
[J].DOI:10.1111/j.1365-246X.1993.tb06977.x URL [本文引用: 1]
An approximate inversion algorithm for time-domain electromagnetic surveys
[J].DOI:10.1016/S0926-9851(99)00023-3 URL [本文引用: 1]
Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit
[J].DOI:10.1190/geo2011-0194.1 URL [本文引用: 1]
1-D multimodel joint inversion of TEM-data over multidimensional structures
[J].DOI:10.1111/j.1365-246X.2008.03973.x URL [本文引用: 1]
Appraisal of a new 1D weighted joint inversion of ground based and helicopter-borne electromagnetic data
[J].DOI:10.1111/1365-2478.12091 URL [本文引用: 1]
Joint inversion of long-offset and central-loop transient electromagnetic data: Application to a mud volcano exploration in Perekishkul, Azerbaijan
[J].DOI:10.1111/1365-2478.12157 URL [本文引用: 1]
Innovative boat-towed transient electromagnetics—Investigation of the Furnas volcanic lake hydrothermal system, Azores
[J].DOI:10.1190/geo2019-0292.1 URL [本文引用: 1]
A resolution study of buried valleys using laterally constrained inversion of TEM data
[J].DOI:10.1016/j.jappgeo.2008.03.003 URL [本文引用: 1]
A parallel, scalable and memory efficient inversion code for very large-scale airborne electromagnetics surveys
[J].DOI:10.1111/1365-2478.12200 URL [本文引用: 1]
时间域航空电磁数据加权横向约束反演
[J].
Weighted laterally-constrained inversion of time-domain airborne electromagnetic data
[J].
航空电磁拟三维模型空间约束反演
[J].
Spatially constrained inversion for airborne EM data using quasi-3D models
[J].
多通道瞬变电磁m序列全时正演模拟与反演
[J].
Multi-transient EM full-time forward modeling and inversion of m-sequences
[J].
Fast imaging of time-domain airborne EM data using deep learning technology
[J].DOI:10.1190/geo2019-0015.1 URL [本文引用: 1]
A generic 1D forward modeling and inversion algorithm for TEM sounding with an arbitrary horizontal loop
[J].DOI:10.1007/s00024-016-1336-6 URL [本文引用: 1]
Transient electromagnetic inversion based on the PSO-DLS combination algorithm
[J].DOI:10.1080/08123985.2019.1627172 URL [本文引用: 1]
基于L1范数的瞬变电磁非线性反演
[J].
L1-norm based nonlinear inversion of transient electromagnetic data
[J].
1D inversion of multicomponent, multifrequency marine CSEM data: Methodology and synthetic studies for resolving thin resistive layers
[J].
Fourier cosine and sine transforms using lagged convolutions in double-precision (subprograms DLAGF0/DLAGF1)
[R].
基于Gaver-Stehfest算法的矩形发射回线激发的瞬变电磁场
[J].
Rectangular loop transient electromagnetic field expressed by Gaver-Stehfest algorithm
[J].
Three effective inverse Laplace transform algorithms for computing time-domain electromagnetic responses
[J].DOI:10.1190/geo2015-0174.1 URL [本文引用: 1]
/
〈 |
|
〉 |
