image
Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 
Online first  /  Current issue  /  All issues  /  Top access  /  RSS

Progress in Geophysics >

2025 , Vol. 40 >Issue 1: 304 - 317

DOI: https://doi.org/10.6038/pg2025HH0284

Research on quantitative characterization method of near-well anomalies in transient electromagnetic logging

  • YuXuan LIU , 1, 2 ,
  • YiRen FAN , 1, 2, * ,
  • XuFei HU 3 ,
  • ShaoGui DENG 1, 2 ,
  • DongYue ZHAO 1, 2
Expand
  • 1 School of Geosciences, China University of Petroleum, Qingdao 266580, China
  • 2 CNPC Key Laboratory for Well Logging, China University of Petroleum, Qingdao 266580, China
  • 3 School of Civil Engineering and Architecture, Wuyi University, Jiangmen 529020, China

Received date: 2023-10-25

  Online published: 2025-03-13

Copyright

Copyright ©2025 Progress in Geophysics. All rights reserved.
Fold

Abstract

Transient electromagnetic wave logging is one of the important methods for subsurface media detection, especially suitable for detecting anomalies near wells. However, the electromagnetic response of wellbore anomalies is diverse and complex, demanding a precise positioning and quantitative characterization method. This method needs to accurately identify key attributes of the anomalies, such as their position, shape, size, and electrical properties, and provide quantitative descriptions. In this study, the three-dimensional Finite Difference Time Domain(FDTD)method in the time domain is employed to investigate the transient electromagnetic field responses of homogeneous formations, horizontally layered anomalies, and three-dimensional block anomalies. Factors such as anomaly resistivity and position are examined to understand their influence on the measurement responses. By analyzing the responses generated by scattering bodies, important information regarding the extension of layered media and the position, size, and response magnitude of three-dimensional isolated anomalies is revealed. A layered wellbore anomaly electrical profile inversion technique is proposed based on the full-area apparent resistivity inversion method. To address the issue of late signals in the full-area apparent resistivity inversion method that cannot accurately reflect the true variation of formation resistivity, a direct vertical partition inversion combining gradient optimization algorithm is proposed. Additionally, a smoke ring inversion algorithm is employed for three-dimensional block anomaly electrical profile inversion. The improved detection method achieves quantitative characterization of wellbore anomalies, significantly improving the accuracy of positioning and electrical property extraction of isolated anomalies near wells, with an overall characterization error of less than 5%. This improved detection method is of great significance for petroleum and solid mineral exploration, aiding in determining the distribution and reserves of oil, gas, and solid mineral resources. It can guide decision-making in exploration activities and optimize resource development plans, thereby enhancing the accuracy and reliability of wellbore anomaly detection. Simultaneously, accurate detection of underground anomalies contributes to understanding the characteristics of subsurface structures, guiding engineering design and construction, and improving the efficiency and safety of engineering projects.

Key words： Transient electromagnetic wave logging; Finite Difference Time Domain(FDTD); Wellbore anomalies; Apparent resistivity inversion; Quantitative anomaly characterization

Cite this article

YuXuan LIU , YiRen FAN , XuFei HU , ShaoGui DENG , DongYue ZHAO . Research on quantitative characterization method of near-well anomalies in transient electromagnetic logging[J]. Progress in Geophysics, 2025 , 40(1) : 304 -317 . DOI: 10.6038/pg2025HH0284

0 引言

瞬变电磁波测井是一种时间域的地球物理勘探方法(朴化荣, 1990; 李貅, 2002; 沈建国等, 2016)，利用脉冲电流激励发射线圈，并通过测量在地层中产生的感应涡流的扩散过程来提取地层的电阻率信息.这种技术已经在近地表矿产资源勘探、地下水勘察、隧道超前地质预报以及矿井水害探测等领域得到了广泛应用(严良俊等, 2003; 孙怀凤, 2013; 武强, 2014).在理解地下介质的物理性质、水文地质特征以及评价和开发矿产资源方面，井旁异常体的探测和识别至关重要(Irvin, 1987; 张红权等, 2019; Yi et al., 2023).异常体可能是地下矿体、地下水体或其他地质构造的表现，其准确定位和定量表征对于资源勘探、环境监测和工程建设具有重要意义.在这方面，瞬变电磁波测井技术展现出了巨大的潜力(赵越和许枫, 2019; 郭嵩巍和闫强, 2020).然而，瞬变电磁波井旁异常体探测技术仍然面临挑战和问题(黄海昆等, 2022; 柳顺彬等, 2023).井旁异常体的电磁响应多样且复杂，与仪器距离、异常体尺寸和介质性质等因素密切相关.因此，进一步研究异常体与仪器之间的关系，探索电磁信号与异常体特征的定量关联，以提高异常体的准确定位和精确识别能力显得尤为重要.
当前，资源勘探、环境监测和工程建设等领域对地下信息的需求迫切，急需一种能够精确定位和定量表征异常体的方法.这种方法应能够准确识别异常体的位置、形状、尺寸和电性质等关键属性，并能够提供定量化的描述.这种方法应能够准确识别异常体的位置、形状、尺寸和电性质等关键属性，并能够提供定量化的描述.本文通过瞬变电磁波测井技术，深入研究了井旁异常体的响应特征，并结合适当的反演算法和模型，提出了一种瞬变电磁波测井井旁异常体定量表征方法.该方法克服了异常体识别中的困难，实现了异常体的精确定位和准确表征.这种改进的探测方法对石油和固体矿产勘探具有重要意义，同时为地下管道、隧道和其他地下设施的规划和建设提供关键支持.准确探测地下异常体有助于了解地下结构特性，指导工程设计和施工，提高工程建设的效率和安全性.此外，该研究成果对地下防御工程的构建、国家安全和城市防护具有重要意义.

1 瞬变电磁波井旁异常体探测响应

首先，均匀地层的zz分量感应电动势表达式(吴柏志等, 2022)为：
$V(t)=\frac{N S \mu M}{\pi^{3 / 2} t} \theta^3 \mathrm{e}^{-\theta^2 L^2}，$
其中，N表示线圈匝数，S表示线圈面积，磁矩M=IS，I表示通电电流，L表示发射接收线圈源距，$\theta=$ $\left(\frac{\mu \sigma}{4 t}\right)^{1 / 2}$.由于：
$\mathrm{e}^{-\theta^2 L^2}=1-\theta^2 L^2+\frac{\left(\theta^2 L^2\right)}{2}-\frac{\left(\theta^2 L^2\right)}{3!}+\cdots，$
当时间t取较大值，即测量晚期时，公式(1)可以简化为：
$V(t) \cong \frac{N S \mu M}{\pi^{3 / 2} t} \theta^3 .$
基于时域有限差分方法(3D-FDTD)，建立三维井旁异常体探测模型进行正演计算，总结其响应规律.首先研究均匀无限厚地层瞬变电磁场探测响应特征，在实验中，设定仪器源距为1.016 m，并分别考虑地层电阻率为1、5、10、20、50、100、200、500 Ω · m的情况，图 1展示了相应地层电阻率下的zz分量感应电动势曲线.可以发现，地层电阻率越高，电磁波在地层中的传播速度越快，因此感应电动势的峰值出现的时间越早，并且峰值的幅度随着电阻率的增加而增大.
图1 不同电阻率的均匀地层感应电动势曲线

Fig 1 Induced electromotive force curve of homogeneous layers with different resistivities

为了研究层状井旁异常体在瞬变电磁场中的探测响应，采用了三层地层模型.根据地层情况的不同，将其分为高阻夹层异常体和低阻夹层异常体两类.在模拟过程中，假设两侧介质具有相同的电阻率，并采用垂向牵拉仪器模拟不同深度点的瞬变电磁场测量响应.假设仪器源距为1.016 m，表 1给出了地层模型中低阻夹层与高阻夹层情况的电阻率参数.当仪器位于高阻夹层中间时，图 2展示了层状异常体模型的瞬变电磁场总场及散射场响应.由图 2a、b可以看出，总场测量响应的早期信号在低阻井旁异常体处出现异常，可以指示目标体的存在，但晚期信号受两侧围岩影响大.需要注意的是，早期信号纵向分辨率高，径向上难以直接通过响应判断目标体延伸.相比而言，当剔除了井两侧的背景介质的场贡献后，井旁的异常体散射响应更能够明显地表现出层状异常体的存在，但是其在早期的分辨率仍弱于总场信号.然而，对比图 2b、c，相比于低阻夹层，高阻夹层时总场响应的早期信号具有较高的分辨率和分辨能力，但散射体贡献更易受两侧低阻围岩影响，导致其延伸难以定性判断.若夹层厚度为5 m，图 2e、f给出了高阻井旁异常体瞬变电磁场总响应与散射体响应贡献谱图.我们可以得到与之前类似的结论，即总场响应分辨率更高，散射体贡献谱图则有助于识别目标体的延伸.
表1 地层参数

Table 1 Stratigraphic parameters

地层模型 R1/(Ω·m) R2/(Ω·m) R3/(Ω·m)
高阻夹层 1 10 1
低阻夹层 10 1 10
图2 层状模型瞬变电磁场总场及井旁异常体散射体响应贡献

(a) 低阻夹层厚度为2 m时总场响应; (b) 低阻夹层厚度为2 m时散射体贡献; (c) 高阻夹层厚度为2 m时总场响应; (d) 高阻夹层厚度为2 m时散射体贡献; (e) 高阻夹层厚度为5 m时总场响应; (f) 高阻夹层厚度为5 m时散射体贡献.

Fig 2 Total field and scattered response contribution from the anomaly in the layered model near the wellbore

(a) Total field response with a 2 m thickness of low resistivity interlayer; (b) Scattered response contribution with a 2 m thickness of low resistivity interlayer; (c) Total field response with a 2 m thickness of high resistivity interlayer; (d) Scattered response contribution with a 2 m thickness of high resistivity interlayer; (e) Total field response with a 5 m thickness of high resistivity interlayer; (f) Scattered response contribution with a 5 m thickness of high resistivity interlayer.

当井/孔周附近存在孤立的异常体时，同样采用3D-FDTD方法模拟孤立块状异常体瞬变电磁波响应，从总测量响应中剔除背景介质的贡献，获得三维孤立块状异常体引起的散射部分的贡献.通过分析散射响应，可以得到以下结论：
(1) 图 3a展示了散射响应的二维谱图.异常体的位置和大小在谱图中表现为特定的异常信号区域.这个区域的纵向位置与异常体的位置一致，中心点对应于异常体的中心，半幅点对应于异常体的上下边界.
图3 井旁三维块状异常体响应

(a) 异常体距离接收线圈2 m; (b) 异常体距离接收线圈7 m; (c) 两块孤立异常体散射体.

Fig 3 3D block anomaly response near the wellbore

(a) Anomaly body at a distance of 2 m from the receiver coil; (b) Anomaly body at a distance of 7 m from the receiver coil; (c) Scattered response of two isolated anomaly bodies.

(2) 利用不同深度位置的异常体进行测量响应的二维谱分析可以准确探测异常体的位置、大小和规模.异常区域的时间轴分布范围主要受异常体径向位置的影响.
(3) 移动异常体到距离发射轴右侧7 m的位置后，进行相同测量方法的二维谱分析，如图 3b所示，仍然可以清晰地识别异常体的纵向中心位置，表明散射响应仍然有效.
(4) 注意到异常体的轴向位置不变，但在时间轴上向后推移.同时，与异常体距离接收线圈2 m时相比，散射响应区域明显减小，但仍然能够清晰地辨认异常体的中心位置.
综上所述，散射响应的二维谱分析能够准确探测异常体的位置和影响，即使异常体位置发生变化.
此外，当距离仪器右侧2 m处存在两个孤立异常体时，它们具有相同的电阻率、大小和径向距离.一个位于接收线圈上方，另一个位于下方，它们距离接收线圈的垂直距离均为3 m.通过排除背景介质引起的场贡献，得到了图 3c中的二维谱.从图中清晰可见，该方法能直观检测异常体的位置、大小和响应规模.井旁的两个异常体在图中的异常区域明显分离，可以清晰识别它们的存在.

2 瞬变电磁波视电阻率提取方法

瞬变电磁波测井的测量信号为时域感应电动势信号.在本节中，将重点研究如何从感应电动势曲线中提取地层电阻率.这一过程主要包括两种方法：晚期视电阻率量化公式方法和全区视电阻率反演方法，下面将对这两种方法进行详细的介绍.

2.1 晚期视电阻率量化公式方法

根据公式(3)，定义瞬变电磁波测井的晚期视电阻率ρlate time，并给出无限厚各向同性地层的晚期视电阻率ρlate-time与感应电动势信号V(t)之间存在的关系(袁习勇等, 2020)：
$\rho_{\text {late-time }}=\left(V(t) * \frac{8 \pi^{3 / 2} t^{5 / 2}}{N S \mu^{5 / 2} m}\right)^{-2 / 3}，$
在测量时域感应电动势的同时，可以利用上述关系得到一条时域视电阻率曲线.由于该曲线是在满足晚期感应电动势简化条件的情况下获得的，因此将其称为晚期视电阻率曲线.实际上，可以使用公式(4)处理整个测量周期的感应电动势数据，只是由此得到的测量早期的视电阻率可能与地层真实电阻率存在较大差异，只有在测量晚期才能等于地层的真实电阻率.这是因为晚期感应电动势的衰减满足以下关系：
$V_{\text {late-time }}(t) \propto \frac{\sigma^{3 / 2}}{t^{5 / 2}}，$
晚期的感应电动势仅是时间t和地层电导率σ的函数.
尽管通过晚期视电阻率计算公式(4)可以得到一条时域视电阻率曲线，但测量早期不满足该算法的适用条件，导致早期时刻点的视电阻率与地层真实情况差异很大.针对这个问题，不同领域如近地表和矿井等发展了各自的全区视电阻率算法.本节提出了针对瞬变电磁测井的全区视电阻率反演方法.
首先，将公式(1)改写为：
$V(t)=\frac{N S \mu M}{\pi^{3 / 2} t L^3} * Y(\varphi)，$
其中，$Y(\varphi)=\varphi^3 \mathrm{e}^{-\varphi^2}，\varphi=\theta L$.
然后，构建全区视电阻率反演的代价函数：
$F(\varphi)=Y(\varphi)-V_{\mathrm{a}}(t) \cdot \frac{\pi^{3 / 2} t L^3}{N S \mu M}，$
其中，Va(t)表示实测感应电动势数据.
根据Halley反演方法(Meng et al., 2018)，确定迭代步长：
$\varphi_{n+1}=\varphi_n-\frac{F\left(\varphi_n\right)}{F^{\prime}\left(\varphi_n\right)-\frac{F^{\prime \prime}\left(\varphi_n\right) F\left(\varphi_n\right)}{2 F^{\prime}\left(\varphi_n\right)}}，$
当目标函数达到所需精度时(如0.0001)，停止参数的更新，此时，全区视电阻率ρalltime由式(9)给出：
$\rho_{\text {all-time }}=\frac{\mu L^2}{4 \varphi^2 t}，$
在下一节中，将以几种典型的地层模型阐述晚期视电阻率公式和全区视电阻率反演算法的处理效果.

2.2 均匀介质视电阻率提取效果分析

首先，以均匀地层为例，说明晚期视电阻率公式(4)的计算效果.利用该式处理图 1的感应电动势数据，生成时域视电阻率曲线.图 4a中显示了不同电阻率均匀地层的晚期视电阻率曲线，表明只有在时间足够长以满足晚期视电阻率公式条件时才能准确获得地层电阻率值.根据涡流场在高阻层中扩散速度较快的特性，高阻层能在较早的时刻满足晚期条件.因此，随着地层电阻率的增大，视电阻率也能在更早的时刻趋近于地层的真实电阻率.
图4 均匀地层的视电阻率曲线

(a)晚期视电阻率曲线; (b) 全区视电阻率曲线.

Fig 4 Resistivity curve of homogeneous layers

(a) Late-time resistivity curve; (b) Regional resistivity curve.

图 4b展示了上述地层的全区视电阻率曲线.与图 4a进行对比可以发现，经过全区视电阻率反演后，在任意时刻都能获得准确的地层电阻率.这表明全区视电阻率反演方法能够克服晚期视电阻率公式在早期测量数据处理时存在的误差问题，并提供更为准确的地层电阻率结果.
在上述视电阻率提取的例子中，并未考虑到噪声的影响，而实际测量中通常存在噪声.假设噪声阈值为1e-8V，图 5a展示了加入阈值噪声前后的视电阻率提取结果.模拟中，源距为1.016 m，地层电阻率为10 Ω · m.可以观察到，加入噪声后，两种方法的晚期信号提取结果与地层真实值之间存在较大的差异，受到噪声的严重影响.特别是对于晚期视电阻率公式法，其早期和晚期效果均严重偏离地层真实情况，因此该方法基本不能使用.相比之下，全区视电阻率反演提取方法在测量早期和中间时间段的提取效果都非常稳定.图 5b进一步展示了考虑实际噪声时的反演结果，其中实际噪声类似于高斯白噪声.可以观察到，全区视电阻率曲线的早期信号能够准确确定地层电阻率.整体而言，该方法在早期阶段的提取结果可达到95%以上，尽管受到晚期影响，整体电阻率的定量表征误差在10%以内.
图5 考虑噪声时的视电阻率提取结果

(a)噪声阈值为1e-8V;(b) 实际噪声情况下.

Fig 5 Resistivity extraction results considering noise

(a) Noise threshold of 1e-8V;(b) Results with actual noise conditions.

2.3 层状介质视电阻率提取效果分析

对于给定的三层层状模型，其中中间层和两侧围岩的电阻率分别为1 Ω · m和10 Ω · m，采用全区视电阻率提取算法循环处理每个深度点的测量电动势谱，从而得到全区视电阻率的二维成像图.当两侧围岩为半空间无限厚地层，中间层层厚分别为2 m、5 m和10 m时，经过处理后获得的二维视电阻率谱如图 6所示.
图6 不同层厚低阻夹层所提取的全区视电阻率二维图

(a)中间层厚为2 m; (b) 中间层厚为5 m; (c)中间层厚为10 m.

Fig 6 2D map of regional resistivity extracted with different thicknesses of low resistivity interlayers

(a) Interlayer thickness of 2 m; (b) Interlayer thickness of 5 m; (c) Interlayer thickness of 10 m.

经过全区视电阻率反演后，观察到提取的视电阻率图像在早期阶段能够直观地反映中间层和围岩的真实电阻率以及层厚.早期信号能够准确地表示地层界面位置和电阻率纵向剖面分布.随着测量时间的增加，视电阻率谱开始发散，表明两侧围岩的影响逐渐增大，此时无法准确确定纵向电阻率剖面分布.较厚的中间层会在视电阻率响应中受到围岩的影响较晚.在测量信号的晚期阶段，电磁波扩散范围增大，夹层的影响较小，此时视电阻率响应基本呈直线.
同样地，假设中间层为高阻夹层，夹层厚度分别为2.0 m和10 m时.经过提取后的层状介质视电阻率二维谱如图 7所示.与低阻中间层类似，通过早期视电阻率信号可以准确获取地层纵向电性剖面分布和地层界面位置.当涡流遇到两侧低阻界面时，视电阻率开始下降；在测量晚期，涡流已经进入两侧围岩中，此时视电阻率主要反映围岩的电阻率.与低阻中间层相比，高阻中间层的视电阻率谱异常区域更短，即狭长的中间层在横向时间上的延伸更短.
图7 不同层厚高阻夹层所提取的全区视电阻率二维图

(a) 中间层厚为2 m; (b) 中间层厚为10 m.

Fig 7 2D map of regional resistivity extracted with different thicknesses of high resistivity interlayers

(a) Interlayer thickness of 2 m; (b) Interlayer thickness of 10 m.

3 层状井旁异常体定量表征方法

根据上述研究结果可知，视电阻率中-晚期信号受到地层层厚的显著影响，并随时间推移，实测的视电阻率已无法准确反映地层电阻率的真实变化.为准确提取层状井旁异常体分布，可以利用视电阻率二维图所提供的已知信息作为反演初值，并结合梯度寻优算法进行直接的地层纵向剖分反演.
这种反演方法的目标是通过利用视电阻率图像中获取的已知信息作为初始值，并运用梯度寻优算法进行优化，以获得地层的纵向电阻率剖面分布.该方法能够最大限度地利用视电阻率图像中的信息，提高对地下目标体分布的准确性.通过多次迭代优化过程，可以逐步改善反演结果，实现对层状井旁异常体的准确提取.

3.1 梯度反演算法

一般而言，对电磁波测井资料的反演可以转化为求实测资料与模拟资料拟合差最小值，即最小二乘反演.为了抑制噪声对反演结果的影响，在反演过程中引入了正则化参数.反演所用代价函数C(m)如下：
$C(\boldsymbol{m})=\frac{1}{2}\left\{\left\|S(\boldsymbol{m})-\boldsymbol{d}^{\mathrm{obs}}\right\|^2+\lambda\left\|\boldsymbol{m}-\boldsymbol{m}_{\text {ref }}\right\|^2\right\}, $
其中，第一项是模拟仪器响应向量S(m)与实测电磁波测井资料向量 d obs之差的L2范数，m表示数据拟合程度.$S(\boldsymbol{m})=\left\{S_1(\boldsymbol{m})，S_2(\boldsymbol{m})，\cdots，S_{M-1}(\boldsymbol{m})\right.$ ， $\left.S_M(\boldsymbol{m})^T\right\}，\boldsymbol{d}^{\text {obs }}=\left\{d_1，d_2，\cdots，d_{M-1}，d_M\right\}^T$ . M表示测量数据向量的大小，$\boldsymbol{m}=\left\{m_1，m_2，\cdots，m_{N-1}，m_N\right\}$，N表示待反演参数的个数，T表示向量或矩阵的转置.代价函数的第二项是正则化项，主要用于抑制测量噪声和减小非线性反演中的病态问题.正则化项由正则化参数(和模型向量与参考向量 m ref之差组成.
为了解决非线性反演问题，可以对代价函数进行二阶Taylor展开，取二次模型的前三项，将代价函数转化为如下形式：
$\begin{aligned}C(\boldsymbol{m}+\Delta \boldsymbol{m})= & C(\boldsymbol{m})+\nabla C(\boldsymbol{m}) \cdot \Delta \boldsymbol{m}+ \\& \frac{1}{2}(\Delta \boldsymbol{m})^{\mathrm{T}} \cdot \nabla \nabla C(\boldsymbol{m}) \cdot \Delta \boldsymbol{m}，\end{aligned}$
其中，Δm表示模型向量 m附近的扰动.$\nabla C(\boldsymbol{m})$表示代价函数C(m)的梯度：
$\begin{aligned}\nabla C(\boldsymbol{m})= & {\left[\begin{array}{c}\frac{\partial C(\boldsymbol{m})}{\partial m_1} \\\vdots \\\frac{\partial C(\boldsymbol{m})}{\partial m_N}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial S_1}{\partial m_1} & \cdots \frac{\partial S_M}{\partial m_1} \\\vdots & \vdots \\\frac{\partial S_1}{\partial m_N} & \cdots \frac{\partial S_M}{\partial m_N}\end{array}\right]\left[\begin{array}{c}S_1-d_1 \\\vdots \\S_M-d_M\end{array}\right]+} \\& \lambda\left(\boldsymbol{m}-\boldsymbol{m}_{\mathrm{ref}}\right)，\end{aligned}$
J(m)为雅克比矩阵.$\nabla \nabla C(\boldsymbol{m})$表示维数为N*N阶的Hessian矩阵，具体表示为：
$\nabla \nabla C(\boldsymbol{m})=\left[\begin{array}{cccc}\frac{\partial^2 C(\boldsymbol{m})}{\partial m_1^2} & \frac{\partial^2 C(\boldsymbol{m})}{\partial m_1 \partial m_2} & \cdots & \frac{\partial^2 C(\boldsymbol{m})}{\partial m_1 \partial m_N} \\\frac{\partial^2 C(\boldsymbol{m})}{\partial m_2 \partial m_1} & \frac{\partial^2 C(\boldsymbol{m})}{\partial m_2^2} & \cdots & \frac{\partial^2 C(\boldsymbol{m})}{\partial m_2 \partial m_N} \\\vdots & & & \vdots \\\frac{\partial^2 C(\boldsymbol{m})}{\partial m_N \partial m_1} & \frac{\partial^2 C(\boldsymbol{m})}{\partial m_N \partial m_2} & \cdots & \frac{\partial^2 C(\boldsymbol{m})}{\partial m_N^2}\end{array}\right]+\lambda \boldsymbol{I}，$
其中，I为单位矩阵.对上述式子进行整理，可以得到如下形式：
$\begin{aligned}\nabla \nabla C(\boldsymbol{m})= & \boldsymbol{J}^{\mathrm{T}}(\boldsymbol{m}) \cdot \boldsymbol{J}(\boldsymbol{m})+\lambda \boldsymbol{I}+\nabla \boldsymbol{J}^{\mathrm{T}}(\boldsymbol{m}) \cdot \\& \left(S(\boldsymbol{m})-\boldsymbol{d}^{\text {obs }}\right)，\end{aligned}$
其中$\nabla \boldsymbol{J}^{\mathrm{T}}(\boldsymbol{m})$为含仪器各个子响应对待反演参数的二阶导数.在实际应用中，由于计算二次导数的复杂性，通常忽略Hessian矩阵的最后一项贡献，简化代价函数为：
$\nabla \nabla C(\boldsymbol{m})=\boldsymbol{J}^{\mathrm{T}}(\boldsymbol{m}) \cdot \boldsymbol{J}(\boldsymbol{m})+\lambda \boldsymbol{I} .$
为了最大化为了最大化降低反演多解性，通常在迭代过程中对反演参数施加地质和岩石物理约束.在反演迭代过程中，通过不断优化使数据拟合差逐渐减小.迭代终止条件通常分为三种情况：首先，在数据拟合程度达到预设阈值时，停止迭代；其次，当迭代次数达到设定上限时，也终止迭代；最后，如果连续两次迭代的模型向量差变化很小，小于设定阈值，表示反演结果已趋于稳定，此时再次迭代将不具有实际意义，也会终止迭代.

3.2 效果分析

为验证上述反演算法的准确性及处理效果，考虑三层介质情况，其中中间夹层为高阻井旁异常体.假设中间夹层厚度分别为2 m、5 m和10 m.图 8展示了不同时间点处提取的全区视电阻率曲线.
图8 高阻井旁异常体不同时刻全区视电阻率曲线

(a) 中间层厚为2 m; (b) 中间层厚为5 m; (c)中间层厚为10 m.

Fig 8 Regional resistivity curves at different time intervals for high resistivity anomalies near the wellbore

(a) Interlayer thickness of 2 m; (b) Interlayer thickness of 5 m; (c) Interlayer thickness of 10 m.

在模型中，高阻夹层的电阻率设为10 Ω · m，两侧围岩的电阻率设为1 Ω · m.根据图中的曲线可以观察到以下情况：初始时刻，涡流受低阻围岩的影响较小，因此视电阻率与地层真实电阻率基本相等.然而，界面附近的视电阻率曲线受围岩影响较大.随着时间的增加，中间层的视电阻率甚至可能大于地层真实值.在测量信号晚期，围岩的影响变得更加显著，导致中间层的视电阻率曲线远小于地层真实值.随后，这种影响逐渐减弱，在晚期，视电阻率远小于地层真实电阻率.
对于不同的中间夹层厚度，地层厚度越大，受两侧围岩影响越小.在早期和中间时间段，视电阻率值能够较准确地反映中间夹层的地层真实信息.总体而言，基于均匀介质模型提取的视电阻率信息为后续处理提供了视电阻率及界面的初始信息，但无法提供准确的地层纵向电性剖面.
经过采用本节提供的一维反演方法后，对不同时刻测量点进行处理，得到的地层电阻率纵向剖面如图 9所示.从图中可以观察到，经过一维反演后的地层电阻率与地层的真实情况一致.在地层相对较薄的情况下，晚期的反演结果略微偏离地层的真实值.同时，地层厚度越大，反演结果的准确性越高.
图9 一维反演后的地层电阻率纵向剖面图

(a) 中间层厚为2 m; (b) 中间层厚为5 m; (c)中间层厚为10 m.

Fig 9 Vertical cross-sections of lithospheric resistivity after 1 D inversion

(a) Interlayer thickness of 2 m; (b) Interlayer thickness of 5 m; (c) Interlayer thickness of 10 m.

图 10展示了井旁异常体为低阻夹层时各个时刻的视电阻率曲线.与高阻夹层相比，低阻夹层的电阻率更接近地层的真实值，因此更容易确定中间层的信息.
图10 低阻井旁异常体不同时刻全区视电阻率曲线

(a) 中间层厚为2 m; (b) 中间层厚为5 m; (c)中间层厚为10 m.

Fig 10 Regional resistivity curves at different time intervals for low resistivity anomalies near the wellbore

(a) Interlayer thickness of 2 m; (b) Interlayer thickness of 5 m; (c) Interlayer thickness of 10 m.

经过反演处理后，图 11展示了不同厚度的薄层到厚层的反演结果，并与地层的真实模型相比较.可以看出，反演结果与地层真实模型吻合良好，反演精度达到95%以上.
图11 一维反演后的地层电阻率纵向剖面图

(a) 中间层厚为2 m; (b) 中间层厚为5 m; (c)中间层厚为10 m.

Fig 11 Vertical cross-sections of lithospheric resistivity after 1D inversion

(a) Interlayer thickness of 2 m; (b) Interlayer thickness of 5 m; (c) Interlayer thickness of 10 m.

4 三维井旁异常体层状介质电性剖面反演方法

4.1 “烟圈”反演方法

本节介绍了一种用于处理复杂三维目标体的定量表征方法.首先，利用全区视电阻率反演理论和均匀模型假设，将实测总场响应转换成视电阻率二维图.接着，通过早期信号和一维反演技术来还原一维电阻率纵向剖面分布.在此基础上，模拟了在宏观背景条件下的瞬变电磁场电动势谱.然后，通过减去背景引起的测量响应贡献，得到了井旁异常体引起的响应.最终，结合烟圈理论，确定了目标体的大小和位置.
针对某一时刻t，假设发射源所处位置为零点，环形电流沿井轴方向的扩散深度可以通过式(16)表示：
$d=4 \sqrt{\frac{t \rho}{\pi \mu}}，$
其中，ρ表示地层的电阻率.
将扩散深度对时间求导得到环形电流的扩散速度，由式(17)给出：
$v=\frac{2}{\sqrt{\pi}} \sqrt{\frac{\rho}{\mu t}}，$
感应涡流在地层中的扩散速度与地层电阻率有关，电阻率越大，涡流的扩散速度越快.因此，地层电阻率可以表示为：
$\rho=\frac{v^2 \pi \mu t}{4}，$
其中扩散速度可通式(19)计算：
$v=\frac{\Delta d}{\Delta t}=\frac{d_j-d_i}{t_j-t_i}=\frac{4}{\sqrt{\pi \mu}}\left(\frac{\sqrt{t_j R_j}-\sqrt{t_i R_i}}{t_j-t_i}\right)，$
其中，ti和tj表示相邻的采样时间，tj>ti，ρi和ρj表示相邻两时刻的地层视电阻率，来源为上文提到的晚期视电阻率ρlate-time或全区视电阻率ρall-time.
一旦获得扩散速度，代入(19)即可计算地层的似电阻率：
$\rho_{\mathrm{r}}=4\left(\frac{\sqrt{t_j R_j}-\sqrt{t_i R_i}}{t_j-t_i}\right)^2 t_{i j}，$
ρr并不代表地层的真实电阻率，也不是经典意义上的地层视电阻率，而是一种中间结果，在此称之为“似电阻率”.其中，tij为ti和tj的算数平均值.
接下来，需要确定似电阻率对应的地层深度.将ρr对应的视深度表示为：
$H_{\mathrm{r}}=0.353 \frac{d_i+d_j}{2}，$
其中经验系数0.353是经过多次模拟后选择的值，可以根据具体地层情况进行调整.
烟圈反演的本质是基于时域视电阻率的一种变换方法，因此时域视电阻率曲线的质量直接影响着烟圈反演的结果.该方法不需要建立初始模型，反演结果可能存在与地层真实情况的误差，但适用于对原始数据进行初步解释和定性分析地层结构.同时，烟圈反演结果还可以作为其他定量反演方法的良好初始模型.

4.2 算法验证

上文中图 3给出了扣除宏观背景场后，井旁异常体散射贡献.对于图 3a中显示的孤立异常体，可以观察到目标体的纵向中心位置约为0 m，并使用半幅点法确定上下界面位置约为0.8 m，即纵向尺寸误差约为20%.然而，无法直接确定异常体与仪器轴之间的距离.因此，需要首先确定径向异常区域的左右边界.基于烟圈反演理论，在模拟背景介质情况下，可以模拟电磁波传播到左侧和右侧时刻所对应的径向距离，如图 11a所示.由图可以确定异常体左侧边界与右侧边界分别为1.2 m和3.5 m.通过这些信息，完成了井旁异常体位置及尺寸的确定.
同样，对于图 3c所示的两块异常体而言，需要进一步分析.在这种情况下，由于存在两个异常体，其位置和尺寸的确定变得更为复杂.针对这种情况，可以采取以下步骤进行处理：首先，针对每个异常体分别进行前述的纵向中心位置和纵向尺寸的确定，通过半幅点法获得上下界面位置，并计算纵向尺寸误差.接下来，需要确定每个异常体的径向边界.基于烟圈反演理论，模拟电磁波在背景介质中传播到左侧和右侧时的径向距离，如图 12b所示.由图可以直观提取两个异常体在左侧和右侧的径向距离.通过对比这些径向距离，可以确定每个异常体的左右边界位置.
图12 采用“烟圈”反演的异常体电性剖面示意图

(a) 单异常体电性剖面反演; (b) 多异常体电性剖面反演.

Fig 12 Schematic illustration of anomaly electrical profiling using the "smoke-ring" inversion method

(a) Single anomaly electrical profiling inversion; (b) Multiple anomaly electrical profiling inversion.

需要注意的是，对于存在多块异常体的情况，井旁异常体的提取精度可能会进一步降低.这是由于多个异常体之间的相互干扰和重叠导致的.此时，准确提取每个异常体的位置和尺寸将更具挑战性.因此，在处理这种复杂情况时，需要综合考虑多种信息和方法，如多步骤的反演算法、先验知识的引入以及对异常体之间相互作用的进一步研究，以提高井旁异常体的定量表征精度.尽管如此，本文建立的井旁异常体提取方法拟补了现有方法处理效果差和适用范围窄的局限，为后续的高维精细反演奠定了基础.

5 结论

(1) 数值模拟研究表明，基于散射场的瞬变电磁测井方法在无限厚均匀地层中有效识别层状异常体，特别是对于低阻层状异常体效果良好.然而，在处理三维块状井旁异常体时，总场响应对其不敏感，没有明显的变化.但在排除背景介质贡献后，可以成功获取三维井旁块状异常体引起的散射部分的贡献.
(2) 全区视电阻率反演方法提取了稳定和可靠的早期和中间时间段信息.通过该方法，提取的视电阻率图像在早期阶段展示了地层和岩层的电阻率以及层厚度.早期信号准确反映了地层界面和电阻率纵向分布.采用这一方法处理每个深度点的电动势谱，得到的视电阻率图像可用作反演的初值，结合梯度寻优算法进行直接反演，反演精度高达95%以上.基于背景场贡献扣除的三维目标体参数表征方法可直观描述井旁异常体的形状、大小和电性特征，整体误差小于10%.
(3) 综合研究结果表明，基于瞬变电磁波测井的井旁异常体定量表征方法能够准确定位和量化描述异常体，包括其位置、形状、尺寸和电性质等关键信息.这一方法提高了井旁异常体探测的准确性和可靠性，对于确定油气和固体矿产资源的分布和储量具有重要意义.

感谢审稿专家提出的修改意见和编辑部的大力支持！

References

Publishing order | Descend order by publishing year | Descend order by cited within
Guo S W , Yan Q . Application of transient electromagnetic methods in water exploration of the west of Ming'an Town in Wulateqianqi, Inner Mongolia. Science Technology and Engineering, 2020, 20 (7): 2564- 2572.

DOI

Huang H K , Luo T T , He Y T , et al. Application of TEM to detect water accumulation area in a coal mine. Mineral Resources and Geology, 2022, 36 (5): 1004- 1010.

Irvin R J . Drillhole TEM surveys at Thalanga, Queensland. Exploration Geophysics, 1987, 18 (3): 285- 293.

DOI

Li X . Theory and Application of Transient Electromagnetic Sounding. Xi'an: Shaanxi Science and Technology Press, 2002

Liu S B , Li W X , Ran J L , et al. Application of S-ATEM for the survey of hidden disaster factors in coal mine. Xinjiang Geology, 2023, 41 (1): 113- 117.

Meng Q X , Hu X Y , Pan H P , et al. Apparent resistivity for transient electromagnetic induction logging and its correction in radial layer identification. Journal of Applied Geophysics, 2018, 151: 328- 342.

DOI

Piao H R . Principles of Electromagnetic Sounding Method. Beijing: Geological Publishing House, 1990

Shen J G , Meng C , Pi G Y . Transient electromagnetic logging theory—transient induction logging. Progress in Geophysics, 2016, 31 (2): 770- 774.

DOI

Sun H F. 2013. Three-dimensional transient electromagnetic response of water bearing structures in tunnels and prediction of water inrush sources[Ph. D. thesis] (in Chinese). Ji'nan: Shandong University.

Wu B Z , Yuan X Y , Guo T Z , et al. Method for extracting apparent resistivity from transient electromagnetic logging measurements and geosteering potential. Journal of China University of Petroleum (Edition of Natural Science), 2022, 46 (6): 99- 109.

Wu Q . Progress, problems and prospects of prevention and control technology of mine water and reutilization in China. Journal of China Coal Society, 2014, 39 (5): 795- 805.

Yan L J , Xu S Z , Hu W B , et al. All-time apparent vertical conductance interpretation method for central loop transient electromagnetic sounding. Journal of Zhejiang University (Science Edition), 2003, 30 (2): 236- 240.

Yi M J , Jeong S , Johmori A , et al. Three-dimensional inversion of airborne time-domain electromagnetic data for ground sources. Exploration Geophysics, 2023, 54 (4): 353- 361.

DOI

Yuan X Y , Deng S G , Hu X F , et al. Detection of the remote boundary using the transient electromagnetic logging method. Chinese Journal of Geophysics, 2020, 63 (7): 2751- 2761.

DOI

Zhang H Q , Chen G Y , Chen H , et al. Application of transient electromagnetic method in exploration of goaf of gypsum mine in Jiangxi Luotang. Progress in Geophysics, 2019, 34 (5): 2112- 2118.

DOI

Zhao Y , Xu F . 3D modeling of buried UXO detection in shallow sea using TEM. Progress in Geophysics, 2019, 34 (3): 1249- 1255.

DOI

郭 嵩巍 , 闫 强 . 瞬变电磁法在内蒙古乌拉特前旗明安镇西找水工程中的应用. 科学技术与工程, 2020, 20 (7): 2564- 2572.

DOI

黄 海昆 , 罗 腾腾 , 何 玉婷 , 等. 瞬变电磁法在某煤矿采空区积水探测中的应用. 矿产与地质, 2022, 36 (5): 1004- 1010.

李 貅 . 瞬变电磁测深的理论与应用. 西安: 陕西科学技术出版社, 2002

柳 顺彬 , 李 文祥 , 冉 军林 , 等. 半航空瞬变电磁在煤矿隐蔽致灾领域的应用. 新疆地质, 2023, 41 (1): 113- 117.

朴 化荣 . 电磁测深法原理. 北京: 地质出版社, 1990

沈 建国 , 孟 超 , 皮 光玉 . 瞬变电磁测井原理研究——瞬态感应测井. 地球物理学进展, 2016, 31 (2): 770- 774.

DOI

孙怀凤. 2013. 隧道含水构造三维瞬变电磁场响应特征及突水灾害源预报研究[博士论文]. 济南: 山东大学.

吴 柏志 , 袁 习勇 , 郭 同政 , 等. 瞬变电磁测井视电阻率提取方法及地质导向潜力. 中国石油大学学报(自然科学版), 2022, 46 (6): 99- 109.

武 强 . 我国矿井水防控与资源化利用的研究进展、问题和展望. 煤炭学报, 2014, 39 (5): 795- 805.

严 良俊 , 徐 世浙 , 胡 文宝 , 等. 中心回线瞬变电磁测深全区视纵向电导解释方法. 浙江大学学报(理学版), 2003, 30 (2): 236- 240.

袁 习勇 , 邓 少贵 , 胡 旭飞 , 等. 瞬变电磁波测井边界远探测方法研究. 地球物理学报, 2020, 63 (7): 2751- 2761.

DOI

张 红权 , 陈 国玉 , 陈 辉 , 等. 瞬变电磁法在江西罗塘石膏矿采空区勘查中的应用. 地球物理学进展, 2019, 34 (5): 2112- 2118.

DOI

赵 越 , 许 枫 . 瞬变电磁法探测浅海掩埋未爆物的三维仿真研究. 地球物理学进展, 2019, 34 (3): 1249- 1255.

DOI

Options
Outlines

/

〈 〉