电性源短偏移距瞬变电磁地形影响特征分析

2022-03-10齐彦福孙乃泉戚志鹏周建美张明晶

吉林大学学报（地球科学版） 2022年1期

齐彦福，李貅，孙乃泉，戚志鹏，周建美，张明晶，邢涛

1.长安大学地质工程与测绘学院，西安 710054 2.长安大学地球物理场多参数综合模拟实验室，西安 710054 3.山东省交通规划设计院，济南 250031 4.北京探创资源科技有限公司，北京 100071

0 引言

瞬变电磁法是一种建立在电磁感应原理基础上的时间域电磁探测方法[1]。该方法首先发射大功率一次电磁信号，然后在断电间歇观测导电大地中感应涡流产生的二次电磁场，通过分析二次场的空间和时间分布特征获得地下电性结构信息[2-3]。根据发射源装置形式的不同，瞬变电磁系统可以分为磁性源瞬变电磁系统和电性源瞬变电磁系统。磁性源瞬变电磁系统采用不接地回线圈作为发射源，然后利用接收线圈观测感应电压，因此磁性源瞬变电磁系统也被称为回线源瞬变电磁系统。该方法不受接地条件的限制，施工方便，对高阻覆盖层具有较强的穿透能力[4]。电性源瞬变电磁系统采用接地长导线作为发射源，然后采用接地电极和接收线圈观测电场和感应电压，根据收发距的不同可分为长偏移距瞬变电磁[5]和短偏移距瞬变电磁[6]两种工作形式。其中：长偏移距瞬变电磁系统利用数千米长的接地导线向地下发射电流，在3～6倍探测深度的偏移距范围观测电磁响应，其探测深度可达几千米至十几千米[5]；而短偏移距瞬变电磁系统在小于2倍探测深度的偏移距范围内观测纯二次场[4]。相比于前者，短偏移距瞬变电磁系统的收发距更小，收发布置具有更强的灵活性，可根据测区的地形条件灵活布置；由于收发距较小，观测信号具有更高的信噪比，同时对发射机的功率和性能要求极大降低，使得野外施工更加方便[4]。

为此，前人针对地形影响开展了大量工作[14-16]。Li等[13]模拟了起伏地表条件下磁性源地面瞬变电磁响应，结果表明地形对观测响应影响极大，其导致观测响应发生复杂畸变。薛国强等[12]通过分析地形影响的基本规律，针对磁性源中心回线装置开展了地形影响快速校正方法研究，并成功应用于理论和实测数据中。Tang等[17]利用积分方程法模拟了层状大地起伏地表模型水平电偶极子源长偏移距瞬变电磁剖面响应，并分析了地形影响特征。然而，与长偏移距瞬变电磁系统相比，短偏移距瞬变电磁系统由于收发距小，场源影响和地形相互耦合将使得观测响应更加复杂。此外，当前电磁探测技术正在向三维观测模式转变，全面了解三维地形影响的空间分布规律成为提高瞬变电磁数据解释精度的关键，因此开展系统的短偏移距瞬变电磁地形影响特征研究十分必要。

作为分析地形影响的重要工具，三维数值模拟技术近年来发展迅速。目前，三维瞬变电磁正演方法主要分为三大类：1)频时转换算法。该方法首先采用有限差分、有限元、有限体积或积分方程等数值算法模拟较宽频带内的频率域电磁响应，然后利用频时转换技术[18]将频率域响应转换到时间域。2)基于显式差分格式的时间域直接求解方法。时域有限差分算法是此类方法的代表，该算法采用向前差分的时间离散格式，结合交错网格技术建立显式时间迭代方程，通过交替求解电磁场的方式实现瞬变电磁正演[19-20]。然而该方法需要满足Courant稳定性条件，对网格尺寸和时间步长的限制比较严格。3)基于隐式差分格式的时间域直接求解方法。该方法采用后退欧拉离散格式，可以获得无条件稳定的隐式差分方程，极大地放宽了对时间步长和网格尺寸的限制，具有更高的稳定性[21]。根据空间离散方式的不同可分为时间域有限体积法[22]和时间域有限元法[23-24]。与前者相比，时间域有限元法可以采用非规则四面体网格进行模型剖分，因此其在处理起伏地形和复杂地电结构时更具优势。

本文基于非结构时间域有限元算法开展电性源短偏移距瞬变电磁地形影响研究。利用非结构四面体网格的灵活性精细刻画起伏地形，基于后退欧拉离散格式的稳定性合理设定时间步长，提高正演效率。首先简要介绍非结构时间域有限元正演理论，并对我们所编写的三维正演程序进行算法精度验证，然后着重分析地形影响特征，通过模拟大量典型起伏地表模型的三维电磁响应分析起伏地形对短偏移距瞬变电磁观测响应的影响规律。

1 非结构时间域有限元正演理论

我们从时间域麦克斯韦方程组出发，通过消去磁场并忽略位移电流项建立电场扩散方程：

(1)

式中：μ为磁导率；e(r,t)为电场；r为空间位置(x,y,z)矢量；t为时间；σ为电导率；js(r,t)为源电流密度。首先采用非结构四面体网格对整个计算区域进行剖分，利用四面体网格的灵活性精细刻画起伏地形和复杂地电结构。然后利用矢量插值基函数[25]对电场进行空间离散，则任意四面体单元内的电场可表示为

(2)

(3)

式中：A为质量矩阵；C为刚度矩阵；S为电流源项。对于第k个单元：

式中：V(k)为第k个单元的体积；dV为体积分。接下来，我们采用二阶后退欧拉离散格式对式(3)进行时间离散，得到无条件稳定的隐式差分方程：

围绕三聚氰胺奶粉事件这一案例，设计了多层级渐进性的探究性课题供学生小组探究（如图1所示）。多层级渐进性探究型教学模式具有5个特征，分别是“知识由浅入深”“能力培养由低到高”“主动性由弱到强”“探究点的力度逐渐变大”“探究点的难度逐渐增加”，这5个特征以一种循序渐进的方式把学生引入到探究式学习中[5]。

(3A+2ΔtC)e(m+2)=A(4e(m+1)-

e(m))-2ΔtS(m+2)。

(5)

式中：e(m)为m时刻的电场；S(m)为m时刻的源项；Δt为时间步长。为了模拟随地表起伏的复杂场源，我们将接地导线分解成若干段首尾相连的短导线，每段导线近似为一个电偶极子。每个电偶极子可以表示为

js(r,t)=δ(r-rs)I(t)fdl。

(6)

式中：δ为脉冲函数；rs为电偶极子的位置；I(t)为电流强度；f为电流方向；dl为电偶极子长度。

本文采用Dirichlet边界条件[24]，假设计算区域外边界的切向电场分量为0，即

(7)

2 数值精度验证

为了检验本文算法的精度，首先设计了如图1所示的100 Ω·m均匀半空间模型，发射源长度为1 000 m，中心点坐标为(0 m,-500 m, 0 m)，接收点坐标为(0 m, 500 m, 0 m)，发射电流20 A。利用非规则四面体网格将整个计算区域剖分成249 486个单元。为了精细刻画场源附近电磁场的剧烈变化，同时提高观测点处的插值精度，对发射源和接收点处的网格进行了局部加密。核心计算区域为3 000 m×3 000 m×500 m，为了满足Dirichlet边界条件，在每个方向分别设置20 km的扩边。采用我们编写的非结构时间域有限元正演模拟程序模拟了阶跃电流波形dBz/dt和ex响应。计算时间被分成了7段，每段100个时间步长，各段时间步长分别为0.1，0.5，2.5，12.5，62.5，312.5和1 560.0 μs。图2展示了本文三维时间域有限元正演结果和一维解析解的对比情况，可以看到三维正演结果与一维解析解吻合得非常好，最大误差均在5%以内，该结果有效地验证了本文算法的准确性。

图1 均匀半空间模型

图2 均匀半空间模型正演响应和误差分析

3 典型山峰和山谷地形对电性源瞬变电磁响应影响特征分析

3.1 接收点处地形对观测响应的影响

本文首先设计了如图3所示的3个半空间模型来分析接收点附近地形对观测响应的影响特征。山峰和山谷地形的最大起伏均为100 m，地下半空间电阻率为100 Ω·m。发射源沿平行x轴方向布设，长度为1 000 m，发射源的中心点坐标为(0 m,-500 m, 0 m)，发射20 A的阶跃电流。采用面积性观测方式进行测量，点距和线距均为20 m，共71×71=5 041个测点，同时观测dBz/dt和ex响应。山峰地形的峰顶和山谷地形的谷底分别位于(0 m, 500 m,-100 m)和(0 m, 500 m, 100 m)处。

a.山峰模型；b.水平地表模型；c.山谷模型。

图4和图5展示了断电后整个观测区域内不同时刻的dBz/dt和ex响应，图6和图7展示了不同偏移距的dBz/dt和ex剖面响应。根据电性源瞬变电磁扩散理论[27]，在水平地表条件下，当导线源中电流突然关断后，在发射源附近的地下介质中会瞬间产生一个与导线源同向的感应主电流，并在四周很大范围内形成回流电流，以维持原有磁场；随着时间的推移，感应主电流开始向下扩散，且主电流的范围逐渐扩大；根据毕奥萨伐尔定律，在发射源两侧的赤道方向会产生两个dBz/dt响应中心，而且中心位置随着主电流下移逐渐向远离发射源方向移动(图 4f—j)。相比之下，ex响应的分布形态比较简单，最大响应始终位于发射源附近(图5f—j)。

当存在起伏地形时，观测响应变得十分复杂，尤其是早期观测信号出现严重畸变。对于早期dBz/dt响应，当响应中心未到达山峰地形时，山峰地形对dBz/dt响应表现为“吸引作用”(图4a)，在峰顶处出现明显正异常(图6b)，同时在峰顶远离发射源一侧出现了反号现象(图4a和图6c)。随着时间的推移，反号现象迅速消失，但山峰地形对dBz/dt响应的”吸引作用”依然存在(图4b)。当dBz/dt响应中心到达并经过地形时，山峰地形对dBz/dt响应表现为“排斥作用”(图4c—e)，产生负异常(图6b)。然而山谷地形恰恰相反，早期dBz/dt响应中心未到达山谷地形时，山谷地形对dBz/dt响应表现为“排斥作用”(图4k)，在谷底处出现明显负异常(图6j)，同时在谷底靠近发射源一侧出现了反号现象(图4k和图6i)。随着时间的推移，反号现象同样迅速消失。当dBz/dt响应中心到达并经过谷底时，山谷地形对dBz/dt响应表现为“吸引作用”(图4m—o)，产生正异常(图6j)。对于ex响应，不论山峰地形还是山谷地形均未出现变号现象。山峰地形对早期ex响应表现为“吸引作用”(图5a—b)，产生正异常(图7b)，而对晚期ex响应起到“排斥作用”(图5c—e)，产生负异常(图7b)。山谷地形完全相反，山谷地形对早期ex响应表现为“排斥作用”(图5k—l)，产生负异常(图7j)，而对晚期ex响应表现为“吸引作用”(图5m—o)，产生正异常(图7j)。图4e、j、o和图5e、j、o显示晚期dBz/dt和ex响应受到地形的影响仍然十分严重，而且相比于dBz/dt响应，地形对电场ex的影响更明显。此外，我们还可以发现不论是dBz/dt响应还是ex响应，地形对早期和晚期信号的影响规律具有相反特征，dBz/dt响应在早期甚至出现变号现象。可以看出，接收点处地形对早期和晚期观测响应均有严重影响，如果在数据解释时不考虑地形的影响，将导致解释结果存在巨大偏差。

a—d.山峰模型响应；e—h.水平地表模型响应；i—l.山谷模型响应。

上述复杂的地形影响规律是由观测点高程变化和地形与周围介质电性差异共同引起的[28]。不同高度接收位置的电磁响应数值模拟结果(图4和图5)表明，接收点高度的增加将导致早期响应被增强，晚期响应被削弱；而接收点高度降低将产生完全相反的影响。与水平地表模型相比，山峰地形接收点高程明显高于水平地面，而山谷地形接收点高程低于水平地面，这成为影响地形顶点处电磁响应产生畸变的主要因素。除了影响观测点高程之外，地形本身还是一个处在地表附近的异常体，其中山峰是空气中的低阻体，而山谷则是埋藏在地下的高阻体。当电源供电时，在地下空间中会产生与供电电流方向相反的回流电流，由于地形与周围介质间存在电性差，因此在地形两侧会产生积累电荷。当突然断电时，地下会产生与供电电流方向相同的主电流，此时地形两侧的积累电荷也会放电，产生散射电流场。由于山峰地形相当于一个低阻体，地形内部散射电流方向与半空间产生的主电流方向相反，因此在山峰地形靠近源一侧产生与半空间同向的感应电压，使得响应增强，而在远离源一侧产生与半空间反向的感应电压，由于早期散射电流场很强导致感应电压dBz/dt出现变号。同时，根据电流的连续性，山顶地表会产生散射回流，且与主电流方向相同，因此导致早期电场ex响应增强。同理，对于山谷地形，由于是高阻体，产生的散射电流场与半空间产生的主电流方向相同，因此在山谷地形靠近源一侧产生与半空间反向的感应电压，产生变号现象，而在远离源一侧产生与半空间同向的感应电压，使得响应增强。同时，谷底地表会产生与主电流方向相反的散射回流，导致早期电场ex响应被削弱。到了晚期，积累电荷放电结束，由其产生的局部畸变也随之消失，此时接收点高程的影响逐渐占主导地位。

a—e.山峰模型响应；f—j.水平地表模型响应；k—o.山谷模型响应。

3.2 发射源和接收点之间的地形对观测响应的影响

在实际工作过程中，经常会遇到在发射源和测区之间存在地形起伏的情况，为了分析该种情况下地形对观测响应的影响，我们在图6和图7中展示了图3模型中不同偏移距的dBz/dt和ex剖面响应，其中y=500 m剖面恰好经过地形中心，以y=850 m和y=1200 m两条测线来分析发射源和接收点之间的地形对观测响应的影响。可以看出不论是dBz/dt响应还是ex响应，当观测点距离地形较近时(图6c、g、k和图7c、g、k)，观测响应会受到很大影响，但是随着偏移距的增加这种影响迅速减弱(图6d、h、l和图7d、h、l)，因此当观测位置距离地形较远时，发射源和接收点之间地形的影响可以忽略。

3.3 发射源处地形对观测响应的影响

本文设计了如图8所示的发射源附近地表起伏模型，其中山峰和山谷模型的最大地形起伏均为100 m。半空间电阻率为100 Ω·m。发射源平行于x轴方向，长度为1 000 m，沿起伏地表敷设，发射源中心点的水平坐标均为(0 m,-500 m)，发射20 A的阶跃电流。观测系统与图3相同。图9和图10展示了断电后整个观测区域内不同时刻的dBz/dt和ex响应。从图9可以看出观测源处地形对dBz/dt响应形态并没有明显影响，但是对早期响应幅值影响严重(图11a)。其中山峰地形使得早期dBz/dt响应增强，产生正异常(图11a)，而山谷地形的存在削弱了dBz/dt响应，产生负异常(图11a)。对于ex响应，山峰地形未改变观测响应的分布形态，只改变了响应幅值，产生了正异常(图11b),而山谷地形完全改变了早期ex响应的分布形态，并产生了明显的负异常(图10k、l、图11b)。这主要是由发射源和观测点相对高程变化以及电磁波传播路径改变引起的。对于山谷地形，发射源高程低于水平地表，感应电磁波到达接收点的路径增长，而且中间介质为导电大地，因此高频信号被吸收，导致早期响应幅值降低；对于山峰地形，虽然电磁波传播路径增长了，然而中间介质为空气，对电磁波的吸收更弱，因此早期信号增强。随着时间的推移，地形影响逐渐减弱。这是由于早期信号主要集中在地表，且高频信号占主导部分，具有较高的分辨率，此时对地形十分敏感；随着时间的推移，感应涡流逐渐向地下传播，高频成分迅速减弱，低频成分占主导地位，分辨率降低，地形的影响逐渐减弱。此外，我们还可以发现相对于dBz/dt响应，ex响应受发射源处地形的影响更加严重。与接收点处地形相比，发射源处地形的影响主要集中在早期，对晚期信号影响较小。

a.山峰模型；b.水平地表模型；c.山谷模型。

a—e.山峰模型响应；f—j.水平地表模型响应；k—o.山谷模型响应。

图11 发射源附近地表起伏模型水平位置(0 m, 500 m)处dBz/dt和ex响应

3.4 发射源和接收点附近同时存在地形对观测响应的影响

本文进一步考虑了在发射源和接收点附近同时存在地形的情况(图12)，地下电阻率为100 Ω·m的半空间。发射源和接收点处地形最大起伏均为100 m，发射和接收装置与图3和图8一致。从图13可以看出，dBz/dt响应的分布形态与只有接收点附近存在地形时(图4)基本一致。然而，受到发射源和接收点两处地形的影响，ex响应分布形态产生严重畸变，由于两处地形的相互耦合导致早期ex响应出现了变号现象；随着时间的推移，晚期ex响应分布形态与只有接收点附近存在地形时基本一致(图14和图5)。可见，接收点处地形对晚期观测信号的影响要强于发射源处地形的影响。

a.发射源位于山谷，接收点位于山峰；b.水平地表模型；c.发射源位于山峰，接收点位于山谷。

a—e.发射源位于山谷，接收点位于山峰；f—j.水平地表模型；k—o.发射源位于山峰，接收点位于山谷。

4 复杂地形影响特征分析

本文最后设计了如图15所示的复杂起伏地表模型，该模型为实际地形，地形数据来自于陕西省南部。地形最大高差约为450 m，地下半空间电阻率为100 Ω·m，电阻率为1 Ω·m、半径为200 m的良导球体埋藏在半空间中，球体中心水平位置位于(0 m, 600 m)处，顶面埋深约为110 m。发射源长度为1 000 m，沿平行x轴方向布设于y=-650 m处的山谷中，发射源的中心点水平坐标为(0 m,-650 m)，发射20 A的阶跃电流。采用面积性观测方式观测dBz/dt和ex响应。从图16和图17可以看出：相比于水平地表模型(图16k—o和图17k—o)，复杂地形使得观测响应产生严重畸变，早期的dBz/dt和ex响应均出现了变号现象(图16a，b，f，g；图17 a，b，f，g)；随着时间的推移，地形对dBz/dt响应的影响开始减弱，良导球体的影响逐渐凸显，即使存在地形的影响，但是仍然可以有效分辨地下良导目标(图16h—j)；然而，复杂地形对ex的早期和晚期响应均产生了严重影响，地形和异常体的相互耦合使得观测响应产生严重畸变，导致无法对地下目标体进行有效分辨，这给后续的数据解释造成了巨大困难(图17h—j)。这一对比结果表明地形对ex观测信号的影响更加严重。