摘要：地-井瞬变电磁(TEM)方法是在地表发射，探头沿钻孔(井)测量瞬变响应的一种井中物探方法.由于接收探头沿钻孔深入地下，因此能够获得更加可靠的地下目标体信息，尤其当存在低阻覆盖层、浅部矿化等地质干扰，或者勘查深部规模不大之良导矿体时，地-井TEM方法的优势更加明显.相对于磁性源系统，电性源系统发射功率大，辐射面积广，更适合地形复杂地区之深部找矿.本文的目的是研究电性源地-井TEM的近似解释方法，首先给出了层状介质电偶极子在地下的TEM响应，进而通过电偶极子叠加的方式获得了电性源条件下的TEM响应.利用感应磁场与均匀半空间电阻率的单调关系通过反函数定理进行了全域视电阻率定义.理论模型的视电阻率计算结果显示，不同深度测点的视电阻率曲线首支不同，而尾支基本一致.这说明了全域视电阻率是测点周围有限范围内介质的综合反映，所以不同深度测点的视电阻率曲线首支所反映的范围是不同的.绘制了不同时间道视电阻率剖面曲线，以期显示地下电性分布规律.计算结果表明全域视电阻率定义能够基本可靠地反映地下信息，从而为该技术在矿区深部找矿的应用提供有力支持.
关键词：电性源瞬变电磁法； 地-井瞬变电磁法； 三分量响应； 全域视电阻率； 深部探测
1 引言
我国矿区找矿的目标已指向“第二深度空间”，对当前的矿区物探方法深部找矿提出了更高的要求(滕吉文，2010).地-井装置瞬变电磁法是将发射源布置在井孔上方或附近地面上，而沿钻孔置入测量装置，使得接收传感器更接近于产生异常响应的深部导体，从而获得比地面方法响应更强的异常(牛之琏，1992；蒋邦远，1997).由于井中受近地表的干扰小，且利用钻孔进行探测，接收装置更加接近深部矿体，因而能获得钻孔周围数百米范围内的有用地质信息，从而提高了见矿率和找矿效果.在地面电磁法工作因矿体深度大，或者受电性干扰因素(如导电覆盖、浅部硫化物、地表矿化地层等)影响大的地区，地-井TEM法的优势就更加突出(Tom et al.,1987).国外对于井中瞬变电磁法的研究较早，主要集中在20世纪80年代(Barnett,1984；Duncan,1987；Tom and Guido,1987；Lindsay et al., 1987；Fullagar et al., 1987；Zhang and Xiao，2001；James et al．，2005；Nigel et al., 2014).我国在此后也出现了相关的研究(牛之琏，1992；蒋邦远等，1998；张兆京和陈卫，1988；蒋慎君和陈卫，1987).近年来张杰等总结出了一套实用的井中TEM数据采集技术流程，以及从定性分析到半定量解释技术，采用矢量交汇技术实现了地-井TEM井旁异常的快速定位(张杰等，2007，2013，2014)；孟庆鑫等应用时域有限差分，计算了均匀半空间和低阻覆盖下均匀半空间中的低阻板状导体的地-井TEM响应并对响应特征和规律进行了研究分析(孟庆鑫等，2012)；杨毅等实现了基于等效涡流的地-井瞬变电磁纯异常反演(杨毅等，2014).
但是，以上研究均是基于磁性源瞬变电磁法，即发射源为不接地回线.在地形较好地区磁性源装置效果明显，然而作者所在课题组实际开展磁性源地-井TEM工作的钻孔多在地形十分复杂之山区.在此类地区布设五方位回线源发射框(200 m×200 m或400 m×400 m)难度很大，退而求其次所布设的回线源难以达到预期效果.
因此，本文中提出利用电性源发射系统的思路.目前电性源瞬变电磁方法研究中以地面方法较多，其中地面电性源瞬变电磁法中主要以长偏移瞬变电磁法(LOTEM)为主，该方法研究较早，成果较多，已经形成较成熟的方法技术(Strack，1992).薛国强等在前人研究的基础上提出了短偏移瞬变电磁法(SOTEM)，并证实其深部探测能力(薛国强等，2013，2014；陈卫营等，2013，2015).此外，电性源地-空瞬变电磁法已成功应用在地热调查、火山结构、地下水盐渍化及地下水监测等领域(Mogi et al.，1998，2009；Allah et al.，2011；Ito et al.，2013;李肃义等，2013；嵇艳鞠等，2013；Wang et al.，2013).因此，电性源瞬变电磁方法在深部勘探中具备较大的潜力(周平等，2007；李貅等，2015).目前国外商品化的多功能电法仪如GDP、V8系统以及我国自主研制的大功率多功能电法系统均能满足电性源瞬变电磁大功率发射的勘探需求(于生宝等，2014； 底青云等，2015).
本文提出的电性源地-井瞬变电磁法是在地表布置接地发射源向地下供电，将探头置于井(钻孔)中接收三分量瞬变响应的装置方式.在这种发射源条件下，发射源的布设位置可选范围较大(周围数公里范围内)，而且一次布源，可同时兼顾发射源辐射区的多个钻孔，可在很大程度上降低施工难度，提高工作效率.但是受限于理论方法的复杂性，对于电性源地-井瞬变电磁法的解释技术研究较少.目前实际应用多采用回线源发射，还未见到电性源地-井装置在矿区深部找矿的实际应用案例.本文旨在通过研究电性源地-井瞬变电磁的视电阻率解释方法，以期获得较为直观、可靠的初步解释方法，为下一步反演解释提供参考信息，解决地形复杂地形条件下地-井瞬变电磁法实际应用提供技术支持.
本文中给出了电偶极子在地下介质中的TEM响应，因此电性源地-井装置条件下井中TEM响应，便可通过电偶极子叠加的方式获取.电性发射源在地下产生的感应磁场与电阻率的关系在一定范围内是单调的，因此通过反函数关系利用感应磁场进行全域视电阻率定义.文中对系列模型进行视电阻率计算，计算结果表明利用本文提出的全域视电阻率定义方法能够可靠反应地下信息，为电性源地-井装置瞬变电磁法提供了一种有效的解释方法.
2 理论方法
2.1 一维层状介质电偶极子地下TEM响应
图1 N层一维理论模型(Key, 2009)
Fig.1 The N-layered 1D model (Key, 2009)
如图1建立层状地电模型，每层的参数由其电阻率值σi和顶面z坐标值zi确定.采用右手直角坐标系，并假定谐变因子为e-iω t，则电磁场所满足的控制方程为
(1)
(2)
其中，μ为磁导率，ω为角频率，E、B分别为电场和磁场，Js为源电流密度.
在库伦规范条件下，可引入磁矢量势
(3)
(4)
对于图1所示的一维地电模型，电偶极子产生的磁矢量势(Key, 2009)可以通过汉克尔变换获得
(5)
其中，J0为0阶第一类贝塞尔函数，r为水平收发距.求得(5)式后，由(3)、(4)式即可获得电磁响应.
电磁场所满足的控制方程组中将(3)和(4)式代入(2)式，并在x和y方向进行傅里叶变换得
(6)
其中
和ky为波数，仅需求解(6)式即可求得(5)式的解.
对于水平电偶极子，取其矢量势为
并令
(7)
(8)
在第i层中，矢量势的形式为
(9)
Λz,i= cieγi(z-zi+1)+die-γi(z-zi)
(10)
其中zi为第i层顶面深度，γ2=λ2-iω μ σ.地下介质中有递归函数Ri=ai/bi和Si=ci/di.R和S代表横电场和横磁场的反射系数(Ward and Hohmann, 1988).将切向连续条件施加到(9)式和(10)式可得到递归函数的表达式
(11)
其中
(12)
以及
(13)
其中，
(14)
hi=zi+1-zi为第i层厚度.在外边界上有RN=SN=0.于是得到反射系数递推公式
(15)
(16)
(17)
(18)
在均匀半空间条件下，地下介质只需计算第二层系数，有
为验证一维正演公式(数字滤波法)计算的准确性，利用均匀半空间上电偶极子响应的解析解(Kaufman，1987)进行了计算，以Hz分量为例计算频域响应分析其计算精度和误差分布特征.计算参数为：均匀半空间电阻率100 Ωm，发射源AB=100 m，发射电流1 A，AB中心O点坐标(0 m,0 m,0 m)，A(0 m,-50 m,0 m)，B(0 m,50 m,0 m)，测点坐标(500 m,50 m,0 m).其计算结果如图2所示.从图2可以看出，在瞬变电磁常用的[10-5, 1011]rad·s-1频段上，数值解和解析解吻合均较好.在低频时数值计算结果较为稳定，而高频时虚部出现了震荡现象.但由于此时频域磁场响应和正余弦滤波系数均已接近于0，对时-频转换的精度影响不大.图3给出了通过时-频转换得到的磁场值及其解析解的对比，并绘制了误差曲线.结果表明，数值解得到的相对误差不超过1%，能够满足计算需求.
图2 均匀半空间频域响应计算结果对比 (a) 实部; (b) 虚部.
Fig.2 Numerical and analytic solution of half space in frequency domain (a) Real part; (b) Imaginary part.
图3 时域计算结果和精度图 (a) 数值解和解析解对比; (b) 相对计算误差.
Fig.3 TEM response and the accuracy of half space in time domain (a) Numerical and analytic solution of half space in time domain; (b) Relative calculation error curve.
2.2 全域视电阻率定义
电性源条件下井中的瞬变响应可以通过电偶极子在地下介质中的瞬变响应叠加获得.公式推导可以看出，电性源地-井装置条件下，井中的瞬变响应与模型参数关系复杂，不能直接以显式关系式表达.为了解钻孔中瞬变响应与均匀半空间电阻率的关系，选择距离发射源不同位置的两测点(P1、P2，如图4)作为研究对象，通过计算不同电阻率不同时间道的瞬变响应曲线分析其变化规律.发射源中心为坐标原点(0 m,0 m,0 m)，P1测点坐标为(300 m,200 m,-200 m)，P2测点坐标为(500 m,300 m,-600 m).
图4 发射源与计算测点位置关系
Fig.4 Tx and survey stations map
图5、6分别为P1、P2点三分量感应电动势随电阻率变化的曲线对比.图中可以看出感应电动势值随电阻率变化的曲线不是单调函数，利用感应电动势值来定义视电阻率需要采取分段的方式来处理，较为复杂.
图7、8分别为P1、P2两测点三分量感应磁场随电阻率变化的曲线对比.两图中可以看出，三分量感应磁场曲线均是电阻率的单调函数，均可利用以定义视电阻率.本文中选择Z分量响应来定义.图中可以看出，测点距离发射源位置不同，其不同时间道对于电阻率的变化率是不一样的.因此，基于感应磁场定义的视电阻率必须同时考虑位置、时间道.本文采取了全域视电阻率的定义方法(李貅，2015).该定义方法是一种空间上不分远近，时间上不分早晚的定义方法，适合于地-井装置的视电阻率计算.
电性源地-井装置下，井中的瞬变电磁响应可以统一地表达为
(19)
其中，x、y、z为测点坐标，t为采样时间点，ρ为介质的电阻率，I0为阶跃电流的大小，a为发射源几何参数，f代表观测场值关于装置参数(发射源长度、发射电流、观测位置等)和模型参数的函数.此公式为基于均匀半空间的计算公式.将装置参数和观测参数看作函数中的参量，当装置固定不变时，略去这些参量不写，即有
(20)
如果能求取(20)式的反函数，即可获得地-井瞬变电磁视电阻率定义式.考虑到(20)式的反函数表达式求取困难，我们采一种等价方法求取(20)式的反函数值.取
求取(20)式的反函数值等价于求解方程
(21)
假定ρ0为g(ρ)零点ρs的一个近似点，将g(ρ)在ρ0处作Taylor展开，有
(22)
忽略二阶以上高阶项，并进行移项整理，可构造迭代格式有
(23)
式(23)即为迭代求取视电阻率值的基本迭代格式.由于这种视电阻率定义方法，在计算过程中考虑了测点空间坐标、时间等参数，因此属于全域的视电阻率定义方法.这种定义方法实际上是将视电阻率定义问题转变成为了一个数值迭代求解非线性方程的过程.为研究迭代格式(23)的收敛性和收敛速度，我们首先不加证明地列出局部收敛性定理和收敛阶定义：
局部收敛性定理 设x*为φ(x)的不动点，φ′(x)在x*的某个邻域内连续，且
则迭代法
(24)
局部收敛；
收敛阶定义 设迭代过程xk+1=φ(xk)收敛于方程x=φ(x)的根x*，如果当k→∞时迭代误差ek=xk-x*满足渐进关系式
则称该迭代过程是p阶收敛的.特别地，p=1(|C|<1)时称为线性收敛.
迭代格式(23)实际上是迭代法(11)式的一种特殊形式，此时
于是有φ(ρ)在ρs处的导数
图5 P1点感应电动势三分量响应随电阻率变化曲线
Fig.5 Three component ∂B/∂t curves with different resistivity at P1 station
图6 P2点感应电动势三分量随电阻率变化曲线
Fig.6 Three component ∂B/∂t curves with different resistivity at P2 station
图7 P1点感应磁场三分量随电阻率变化曲线
Fig.7 Three component B curves with different resistivity
图8 P2点感应磁场三分量随电阻率变化曲线
Fig.8 Three component B curves with different resistivity
(25)
值的注意的是，由于g(ρ)的定义为场值差的二次方，当方程(5)和(8)严格满足时，则必有g′(ρ)=0，即ρs至少为方程(21)的2重根，此时无法直接确定式(25)的值.假定ρs为方程(21)的m重根,并令g(ρ)=(ρ-ρs)mq(ρ)，且q(ρs)≠0，则有
(26)
于是
迭代法(10)式在ρs附近局部收敛.
下面讨论迭代法(23)式的收敛速度，据(23)式容易写出
(27)
于是，有
(28)
根据收敛阶定义可以判断，迭代法(23)式线性收敛.
3 数值模拟结果分析
为检验全域视电阻率计算方法的正确性，选择系列一维模型进行数值计算.一维模型主要选择D、G、H、K型断面，模型参数如表1所示.发射源AB=100 m，发射电流1 A，AB中心O点坐标(0 m,0 m)，A(0 m,-50 m)，B(0 m,50 m),钻孔坐标(500 m,50 m).钻孔均为直孔，测点深度0～-500 m，点距20 m.
表1 一维层状模型参数表
Table 1 Parameters of 1D layers
由表1中可知D型模型第一层电阻率为200 Ωm，厚度为300 m，第二层电阻率为10 Ωm.图9为D型模型计算结果，其中图9a为不同深度测点正演Z分量衰减曲线，图9a中随着深度增加，衰减曲线早期道逐渐变得平缓，这是由于感应磁场在深部测点上早期变化较小产生的，距离越近磁场变化越剧烈，距离越远变化越不明显.图9b为不同时间道随深度变化剖面曲线图，图中早期道0.1 ms和0.125 ms曲线无明显异常显示，而自0.199～1.99 ms曲线在300 m深度界面处出现明显异常，响应幅值随深度逐渐变大并在达到极值后变小.曲线在5.01 ms之后再无明显异常显示.图9c为不同深度测点计算所得视电阻率曲线，不同颜色曲线代表不同深度处测点计算结果.图9c中可以看出，地表测点视电阻率曲线首支接近第一层电阻率值200 Ωm，之后曲线下降直到接近第二层电阻率值10 Ωm，为典型的D型曲线.而随着深度增加，视电阻率曲线首支越来越不明显，曲线一直随时间减小并最终达到第二层电阻率值，尾支明显.深度360～500 m测点首支早期道视电阻率值仅为几十欧姆米，但是曲线保持了自高到低的变化趋势.分析这种现象，其原因主要是视电阻率值所反映的地下信息不同造成的.地面上测点的视电阻率曲线完全反映了自地面至地下的断面信息，因此曲线完整，与断面完全相符.而地下测点的视电阻率曲线则主要反应测点周围一定范围内的电阻率变化情况，尤其是当测点处于低阻层内时，该测点视电阻率曲线所反映信息则以低阻层为主，计算所得视电阻率值会偏低.但是这种差异仅在早期明显，随时间增加，其差异逐渐变小.
因此基于以上分析，本文将不同深度同一时间道对应的视电阻率值绘制曲线，则该曲线能够反映地下介质沿钻孔的变化情况.图9d中给出了0.199～3.16 ms的视电阻率随深度变化曲线以及理论模型曲线.图中可以看出0.199 ms曲线在100 m以浅幅值接近200 Ωm，之后随深度增加逐渐变小，最终保持在约60 Ωm.该曲线为典型D型曲线，与实际地点断面相符.图中可以看出随着时间增加D型曲线越来越不明显，3.16 ms的曲线已看不出电阻率值随深度的变化.这种现象与之前的分析是一致的.
图10为G型模型计算结果，图a、b分别为Z分量衰减曲线和剖面曲线结果.图c为不同深度测点视电阻率曲线，图中可以看出地表测点0 m以及深度较小测点(如-60 m、-100 m)视电阻率曲线与实际模型较为接近，首支接近第一层电阻率值，位置趋向第二层电阻率值.而随深度增大视电阻率曲线虽呈“G型”，但首支逐渐抬高偏离真实值，其原因亦是由于视电阻率值所反映的范围的缘故.当测点位于地表附近，其视电阻率值主要反映地表以下介质，因此与实际吻合.然而当测点位于地下一定深度，则视电阻率值反映的范围改变，尤其是当测点靠近高阻层时，视电阻率值由于受高阻影响而偏大.图10c中可以看出在早期时间道中视电阻率曲线首支变化明显，而随着时间道增加，尾支趋于一致.图10d中为不同深度测点不同时间道计算的视电阻率剖面曲线，可以看出0.1～2.51 ms之间的曲线均可以反映地下介质电阻率变化情况，而之后电阻率无明显变化趋势，尤其是到100 ms之后视电阻率剖面曲线呈直线并趋向第二层介质的电阻率值.
图9 D型模型计算结果 (a) TEM Z分量响应曲线; (b) Z分量响应随深度剖面曲线; (c) 不同深度测点视电阻率曲线; (d) 视电阻率曲线随深度变化曲线.
Fig.9 Calculation results of D type model (a) Z component TEM decay curves of different depth; (b) Z component TEM plot curves of different times; (c) Apparent resistivity curves of different depth; (d) Apparent resistivity plot curves of different times.
图10 G型模型计算结果 (a) TEM Z分量响应曲线; (b) Z分量响应随深度剖面曲线; (c) 不同深度测点视电阻率曲线; (d)视电阻率曲线随深度变化曲线.
Fig.10 Calculation results of G type model (a) Z component TEM decay curves of different depth; (b) Z component TEM plot curves of different times; (c) Apparent resistivity curves of different depth; (d) Apparent resistivity plot curves of different times.
图11为H型模型计算结果.图11a为正演Z分量衰减曲线，与D型曲线类似，曲线早期道随测点深度增加而逐渐变缓，这是由于深部测点距离场源远，磁场早期变化缓慢造成的.图11b可以看出早期道(0.01～0.05 ms)曲线无异常显示，随着时间道增加(0.1～1.25 ms)异常较为明显，之后(1.99～3.16 ms)异常幅值很小.图11c可以看出不同深度测点计算结果均能显示明显“H”型特征，但是视电阻率曲线首支不同.前文中已分析其原因，这是由于测点深度不同视电阻率值反应的范围不同，当测点深度增加并逐渐靠近低阻层时，早期响应的视电阻率值也会偏低，晚期受其影响较小.图11d中可以明显看出测点靠近低阻层时早期响应对应的视电阻率值会呈现低值.值得注意的是，当测点位于低阻层下方时，低阻层的影响仍然存在，使得早期视电阻率曲线偏低.
图12为K型模型计算结果.图12a为正演Z分量衰减曲线.图12b中在200～300 m深度范围内出现微弱异常显示，这主要是二次场对于高阻层的反应不灵敏.图12c可以看出不同深度测点计算结果中，仅浅部测点能显示明显“K”型特征，深部测点视电阻率曲线仅表现为“D”型特征.这是由于随着测点深度增大并逐渐靠近高阻层时，由于高阻层上下均为低阻层，当测点位于高阻层附近时，该测点处反映的综合信息中低阻地层仍然占有较大的比例，使得计算所得视电阻率值整体偏低.当深度继续增加进入第三层低阻层(-300～-500 m)后，视电阻率曲线主要反映低阻信息，其首支曲线又逐渐变低.因此对于“K”型地层，图12d中给出了0.1～0.797 ms视电阻率剖面曲线，图中曲线略呈“K”型，幅值较小，异常不明显.这是由于图中同时绘制了理论模型曲线，在一定程度上压缩了视电阻率曲线，使得视电阻率对比不太明显.图13是单独绘制的曲线图，图中可以看出大致的曲线特征，早期时间道的视电阻率曲线K型特征明显，随时间增大曲线特征逐渐模糊.其原因主要是高阻层对于电磁场的聚流作用弱于低阻层，电磁场在高阻环境衰减较快，而在低阻环境较缓.
图11 H型模型计算结果 (a) TEM Z分量响应曲线; (b) Z分量响应随深度剖面曲线; (c) 不同深度测点视电阻率曲线; (d) 视电阻率曲线随深度变化曲线.
Fig.11 Calculation results of D type model (a) Z component TEM decay curves of different depth; (b) Z component TEM plot curves of different times; (c) Apparent resistivity curves of different depth; (d) Apparent resistivity plot curves of different times.
图12 K型模型计算结果 (a) TEM Z分量响应曲线; (b) Z分量响应随深度剖面曲线; (c) 不同深度测点视电阻率曲线; (d) 视电阻率曲线随深度变化曲线.
Fig.10 Calculation results of K type model (a) Z component TEM decay curves of different depth; (b) Z component TEM plot curves of different times; (c) Apparent resistivity curves of different depth; (d) Apparent resistivity plot curves of different times.
图13 K型模型视电阻率曲线随深度变化曲线图
Fig.13 Apparent resistivity plot curves of different times of K type model
4 结论
通过前文理论公式、计算方法介绍以及理论模型数值模拟及分析，初步获得了以下基本结论：
(1) 本文针对地-井TEM方法的应用现状和实际需求，研究了电性源地-井TEM装置，该装置是在地面通过接地发射源发射，在井(钻孔)中接收三分量瞬变响应的装置.通过接地发射源可以使用大功率发射机，从而能够覆盖较大的工作范围.能够实现同一发射源辐射矿区多个钻孔，适合于我国目前矿区深部找矿的现状，可以使得野外施工更加方便，更加高效.
(2) 本文基于反函数思想提出了电性源地-井瞬变电磁法的全域视电阻率定义方法，理论模型计算结果表明，文中提出的全域视电阻率曲线能够反映地下电性沿钻孔变化情况.
(3) 计算结果表明，井中视电阻率曲线与地面视电阻率曲线有所不同，地面视电阻率曲线反映自地面至地下一定深度范围内的电性变化情况，而钻孔中的视电阻率值主要反映测点周围的电性，因此钻孔中的视电阻率曲线反应的不是整体地电断面，而是测点周围一定范围内的介质.这个差异在早期道对应的视电阻率曲线中非常明显，随着时间道增加视电阻率曲线与地面计算结果趋于一致.本文利用此特点，将不同深度早期道视电阻率绘制成剖面曲线，能够反映地电断面的变化情况.
