任意形状水平电性源瞬变电磁全区视电阻率计算

商天新，张继锋*,冯 兵，刘长生

(1.长安大学 地质工程与测绘学院，陕西 西安 710054；2.湖南省飞机维修工程技术研究中心，湖南 长沙 410124)

0 引 言

目前，地面瞬变电磁法已经发展出多种发射装置，常见的有中心回线装置、重叠回线装置、大定源回线装置、电性源装置等。然而，回线源瞬变电磁采用晚期视电阻率进行近似计算时有一定的局限性[1]。因此，有必要采用全区视电阻率计算和表征地下电性参数分布。为完成地面上任意发射源全区视电阻率的计算，本文在冯兵等提出的采用均匀半空间电偶极子叠加定义全区视电阻率[2]的基础上开展研究。

电性源瞬变电磁法具有很多优势，近年来越来越受到重视。薛国强等通过对短偏移距瞬变电磁法深入研究，发现其具有勘探深度大、纵向分辨能力高、体积效应小、工作效率高的优势[3-6]；底青云等认为电性源短偏移瞬变电磁法是中国实现深部目标体精细探测的重要手段[7]；陈卫营等认为磁感应强度垂直分量(Bz)及其时间导数(∂Bz/∂t)仅对低阻体敏感，电场强度水平x分量(Ex)对低阻体和高阻体都比较敏感[8]；李展辉等提出了一种适应于任意形状水平接地电性源瞬变电磁法的一维正演方法，正演精度较高[9]。

视电阻率定义的发展主要经历了早期、晚期及全程视电阻率3个阶段。殷长春等较早地给出了电磁法中定义视电阻率的基本原则，并研究了全区视电阻率定义时的若干问题[10]；李貅系统地推导了在满足某种极限条件下的基于磁感应强度垂直分量Bz(∂Bz/∂t)的早期、晚期视电阻率计算公式[11]；白登海等通过对核函数特征的分析，提出一种用计算机迭代求取中心回线装置全程视电阻率的方法[12]；李建平等基于水平电偶极子垂直磁场叠加和数值迭代法计算任意形状回线源的全区视电阻率[13]；王华军通过推导均匀半空间下瞬变电磁场的平移特性，提出了计算回线源中心点全区视电阻率的平移算法[14]；陈清礼运用二分搜索算法计算了基于∂Bz/∂t定义的中心回线装置全区视电阻率[15]；崔江伟等计算电性源瞬变电磁的全区视电阻率时，依据磁感应强度垂直分量与电阻率的单调特性,先将∂Bz/∂t转化为磁感应强度垂直分量，再用二分搜索算法计算了瞬变电磁B场全程视电阻率[16]；赵越等研究了大回线源非中心点处基于瞬变电磁B场定义的全区视电阻率[17]；闫国翔等研究了接地长导线源基于磁场强度垂直分量(Hz)及磁场强度水平y分量(Hy)定义的全区视电阻率，对比分析了两种分量定义的视电阻率曲线特征及分辨高、低阻层的能力[1]。

综上可知：①目前，对回线源及直导线源视电阻率定义研究较多，对复杂发射源研究较少；②对电性源全区视电阻率的计算主要基于反函数的方法，通过研究均匀半空间场与电阻率单调性的关系，采用迭代的方法实现；③对回线源中心点数据处理研究较多，对非中心点数据处理研究较少。基于此，本文研究了任意发射源下磁场强度垂直分量及其时间导数的平移特性，并分别给出了基于上面两种分量定义的全区视电阻率计算方法。同时，通过模型计算与实测数据处理，进一步验证了该方法的有效性。

1 电性源瞬变电磁一维正演

地面任意电性源产生的电磁响应都可以通过对发射源进行电偶极子分割，然后对所有电偶极子积分得到，本节给出了不同角度水平电性源瞬变电磁场的正演方法。

1.1 水平电偶极子的垂直磁场

对于如图1所示的A(x1,y1)→B(x2,y2)水平电偶极子，频率域下地表产生的磁场强度垂直分量表达式[18]为

图1 水平电偶极子示意图

(1)

(2)

(3)

(4)

(5)

(6)

(7)

1.2 地面任意水平电性源在地表的垂直磁场

对如图2所示的电性源，可先将其剖分为Q1、Q2、…、QM等M段长直导线，每一段长直导线又可进一步剖分为N个有限长度的电偶极子，此时产生的磁场强度垂直分量可以通过对所有电偶极子求和近似得到，其表达式为

图2 电性源剖分示意图

(8)

(9)

(10)

式中：Δli为第i个电偶极子的长度；y′i为接收点到第i个电偶极子的垂直距离；Ri为接收点到第i个电偶极子的水平距离。

要计算层状介质下时间域的磁场强度垂直分量(Hz(t))及其时间导数(∂Hz(t)/∂t)，可利用汉克尔积分变换及频时变换[11]处理。本文运用140点的J1汉克尔系数做积分运算，运用160点的正弦系数做磁场强度时间导数运算，运用300点的余弦系数做磁场强度的计算。

1.3 算法的正确性验证

为验证编制的层状地层正演程序的正确性，可以与均匀半空间下水平电偶极子时间域叠加得到的磁场强度垂直分量Hz(t)及其时间导数∂Hz(t)/∂t的解析解进行对比。电偶极子解析解表达式[18]为

(11)

(12)

(13)

(14)

式中：erf(u)为误差函数；u为瞬变场参数；σ为电导率；t为时间。

本次正演程序设计模型如下：发射导线沿y方向布设，长为1 km，发射电流为1 A，发收距为2 km，均匀半空间下电导率取0.01 S·m-1。发射源共剖分为50段，分别采用本文中已讨论过的数值方法和解析方法叠加计算垂直分量，并做误差对比，结果如图3所示。垂直分量的数值解与解析解的相对误差均在1%以内。这表明数值解具有很高的精度，有利于后面全区视电阻率的计算。

图3 垂直分量解析解与数值解对比

2 平移特性分析

2.1 Hz(t)平移特性

从均匀半空间下水平电偶极子时间域下的磁场强度垂直分量[式(11)]可以看出，假设对电导率和时间同时增加K倍，记为σ*、t*，则增加K倍瞬变场参数(u*)的表达式为

(15)

同一个测点都存在以下关系，即

=Hz(σ,R,t)

(16)

(17)

(18)

式中：ka为磁场强度垂直分量平移轨迹斜率。

式(16)～(18)表明，对于每一个电偶极子，在保持装置参数不变的条件下，它们具有同样的性质。当均匀半空间的电导率增加为原来的K倍，同时将观测时长延长为原来的K倍，它们的磁场强度保持不变，即平移特性。在单对数坐标系下，其平移轨迹斜率为0。

(19)

因此，上述对于电偶极子讨论过的关于磁场强度的平移特性对于任意形状的电性源应同样成立。

2.2 ∂Hz(t)/∂t平移特性

对于磁场强度垂直分量时间导数，对电导率和时间同时增加K倍，根据式(12)可以推导出

(20)

(21)

(22)

式中：kb为磁场强度垂直分量时间导数平移轨迹斜率。

式(20)～(22)表明，对每一个电偶极子，将电导率及时间同时增加K倍，观测到的磁场强度时间导数是原来的1/K，此即关于∂Hz(t)/∂t的平移特性。在双对数坐标系下，其平移轨迹斜率为-1。

(23)

因此，上述对于电偶极子讨论过的关于磁场强度时间导数的平移特性对于任意形状的电性源应同样成立。

3 全区视电阻率计算

视电阻率实质是地下不均匀体和地形的综合反映，当地下为非均匀介质时，仍然用均匀大地的公式进行计算，即把测量的感应电动势或磁场看做是某均匀大地模型的结果，该均匀大地的电阻率可称为视电阻率，这是平移算法的理论基础。

3.1 基于Hz(t)定义的全区视电阻率计算

图5 平移算法示意图

(24)

3.2 基于∂Hz(t)/∂t定义的全区视电阻率计算

磁场强度垂直分量时间导数的平移特性如图4(b)所示，其核函数为u的双值函数[19]。在双对数坐标系下，平移轨迹是一条斜率为-1的直线。一条平移轨迹上的所有平移点应该有相同的平移截距(b(t))，其表达式为

(25)

(26)

实测响应平移截距最大值[bobs(t)]max与理论平移截距最大值[b(t)]max存在：①[bobs(t)]max>[b(t)]max,存在一段无解区；②[bobs(t)]max=[b(t)]max，极大值点处存在唯一解；③[bobs(t)]max<[b(t)]max，此时早、晚期可按照前面的方法分别计算，但是极大值点处左、右两支会分别计算出两个不同的视电阻率，全区视电阻率不是连续变化的，会呈现出断接的现象。通常情况下，实测数据会出现①或③的情况。针对无解或断接现象，可以在平移截距极大值附近设置一段过渡区，采用无约束点的最小曲率原位迭代格式插值得到过渡区内平滑变化的视电阻率,从而得到全时段视电阻率曲线。无约束点的最小曲率原位迭代公式[20]为

(27)

其中

4 模型验证

为运用平移算法计算电性源在任意接收点的全区视电阻率，本次研究布置发射源如图6所示，线段ABCD为发射源，A、D为接地端，每段发射源均分别剖分为50个电偶极子。P1～P4为接收点。各点坐标如下：A(-500，1 000)、B(0，0)、C(1 000，0)、D(1 500，500)、P1(0，3 000)、P2(1 000，3 000)、P3(2 000，3 000)、P4(2 500，2 000)。选取多种地电模型计算基于Hz(t)及∂Hz(t)/∂t定义的全区视电阻率。

图6 发射源及接收点位置

将一、二层视电阻率分别记为ρ1、ρ2。两层地层模型全区视电阻率曲线分为两类：当ρ1<ρ2时，为高(G)型；当ρ1>ρ2时，为低(D)型。图7以P2接收点为例，计算了两层曲线的视电阻率变化，一层视电阻率ρ1=100 Ω·m，两层间厚度h1=150 m，ρ2/ρ1分别取10.00、6.00、2.00、1.00、0.50、0.10、0.05，发射电流为1.0 A。从图7可以看出，两种定义方式都大致能反映地层变化情况，但用磁场强度定义的全区视电阻率形态更平滑。用磁场强度时间导数定义的曲线过渡区内由最小曲率插值得到，故过渡区内有一定的误差。

图7 两层地电模型全区视电阻率曲线对比

图8以P2接收点为例，计算了三层H型和K型在不同中间层厚度中的全区视电阻率变化，发射电流为1.0 A。从图8可以看出，无论用哪种方式定义，均能反映出电性层的变化。首支趋于第一层的电阻率，尾支趋于最后一层的电阻率。从对中间电性层的反映来看，无论哪种定义方式，对中间低阻层的反映均要优于高阻层，这也符合电磁类方法的电性特征[21]。对比磁场强度时间导数和磁场强度，可以发现用磁场强度时间导数计算全区视电阻率得到的异常幅值要比磁场强度计算得到的更大。这说明用磁场强度时间导数定义全区视电阻率有更高的灵敏度且更趋近于真实地层变化。

H型模型中，一、二、三层视电阻率分别为100、10、100 Ω·m，一、二层间厚度为150 m，二、三层间厚度h2为变化值，分别取10、30、60、100 m；K型模型中，一、二、三层视电阻率分别为100、1 000、100 Ω·m，一、二层间厚度为150 m，二、三层间厚度h′2为变化值，分别取100、300、500、700 m

图9以P1、P2、P3、P4接收点为例，测试不同偏移距的接收点对同一地层的反应，一、二、三层视电阻率分别为100、10、100 Ω·m，一、二层间厚度为150 m，二、三层间厚度为100 m，发射电流为1.0 A。从图9可以看出，不同偏移距下用两种定义方式计算得到的全区视电阻率均差别很小，形态一致。本文所有计算均考虑了接收和发射导线的几何关系，故保证了任意电性源、任意接收点计算全区视电阻率的可靠性。

图9 H型模型不同偏移距下全区视电阻率曲线对比

5 应用实例

应用实例位于陕西省渭南市华州区华阳乡境内，区内连霍高速和韦罗高速交汇于此，郑西高铁、310国道穿境而过，交通便利。但在工作勘查区内，由于实行封山育林，植被发育，林木丛生，荆棘遍地，山大沟深，切割剧烈，为强烈侵蚀切割的中—陡坡地形，海拔一般为1 300～1 800 m，通行较为困难，工作难度较大。

勘查区地层岩性主要为太华群深变质片麻岩，包括长英质片麻岩、黑云斜长片麻岩、斜长角闪片麻岩。区内分布3个断裂构造带，分别为杏儿岭构造带、鹰嘴崖构造带与冰凌沟构造带，构造带主要为张扭性断裂。断裂带侵入的脉体主要有3类，分别为花岗质岩脉、花岗伟晶岩脉、石英方解石岩脉。其中，石英方解石岩脉含有金属硫化物，呈中低阻。本次研究的主要任务是通过短偏移距电性源瞬变电磁的数据采集与处理，了解勘查区内金属硫化物的赋存情况。

本次研究数据采集使用V8多功能电法仪，接收天线使用美国ZONGE公司生产的TEM-3接收线圈，接收面积为10 000 m2，发射源长度为1 810 m，发射电流为1.0 A，发射基频为8.33 Hz，偏移距为200 m，测线方向与发射源平行，观测点距为25 m。数据处理以Z2测线的一段为例，首先将测得的感应电动势转化为磁场强度时间导数，再按照平移算法计算全区视电阻率。

图10为实测数据处理结果，红色曲线为测区内地表所在的海拔，瞬变电磁探测存在一段盲区。横向上看，视电阻率基本随距离的增大而逐渐增大，在距离为2 650 m附近有比较明显的分界。结合金属硫化物呈中低阻的特点，可以推断距离2 400～2 650 m区段内为富金属硫化物的石英方解石岩脉；距离2 650～2 850 m区段内视电阻率较高，推断其为花岗片麻岩或混合岩类。纵向上看，视电阻率整体呈上低下高，推断石英方解石岩脉的赋存位置位于海拔700～1 200 m。这也与后期用SIP充电率法反演得到的结果基本保持一致。

图10 Z2测线视电阻率断面图

6 结 语

(1)用电偶极子叠加的方式可以实现地面任意电性源的一维正演，且误差均在1%以内，可对形态复杂的发射源进行正演计算。根据平移特性分析，提出了一种计算地面电性源全区视电阻率的方法。其中，基于磁场强度定义的视电阻率是单值且有解，基于磁场强度时间导数定义的视电阻率存在双解或无解情况，本文采用无约束点的最小曲率原位迭代格式插值解决了上述问题。

(2)模型计算表明，全区视电阻率曲线能在全时间范围内客观反映地电断面信息。两种定义方式都对多层介质的中间低阻层较为敏感，对高阻层不敏感。相比于磁场强度，用磁场强度时间导数定义的全区视电阻率有更大的异常幅值，更接近于地电断面的真实电阻率。全区视电阻率的计算与收发距无关，仅与层状介质的电性特征有关。

(3)实测数据的处理结果表明，本次研究采用平移算法计算出的电性源全区视电阻率有实际价值。在实际应用中，可以根据数据的具体格式选择本位的全区视电阻率计算方法。

(4)本文只针对磁场强度及其时间导数的垂直分量开展研究，对于水平分量的特性以及使用水平分量计算全区视电阻率需进一步讨论。