ABSTRACT

Because of the existence of shallow inhomogeneous bodies, static shift occurs in magnetotelluric sounding curves, which makes the interpretation of magnetotelluric data very difficult, or even impossible. It is very important to identify static shift in the interpretation of magnetotelluric data. In this paper, in order to compute the two-dimensional magnetotelluric responses of static shift model, non-uniform grids finite difference method was adopted for numerical results. From the electric field and magnetic field based on the differential equation, the linear equation for magnetotelluric numerical modeling was derived by non-uniform meshes finite difference method. By the numerical simulation for the geo-electric model with shallow inhomogeneous body, the character of static shift can be obtained and the magnetotelluric responses computed by non-uniform girds finite difference method are effective.

Keywords:Magnetotelluric, Two-Dimensional Model, Static Shift, Finite Difference Method, Non-Uniform Meshes

基于非均匀网格有限差分法的大地电磁静位移模拟

邓小康^1,2，童孝忠^3*

¹中南勘测设计研究院水能资源利用关键技术湖南省重点实验室，湖南 长沙

²湖南工业大学土木工程学院，湖南 株洲

³中南大学地球科学与信息物理学院，湖南 长沙

收稿日期：2018年12月3日；录用日期：2018年12月17日；发布日期：2018年12月24日

摘 要

在大地电磁勘探中，视电阻率曲线总会存在一些静态现象，而静态问题又一直是电磁测深领域一个棘手的问题，它直接影响到电磁测量结果的解释，能否正确识别静态是实际勘探工作能否顺利进行的一个关键，同时也关系到测深工作的成败。为了计算二维静位移地电模型的大地电磁响应，本文采用非均匀网格有限差分法进行了数值模拟。从电磁场满足的微分方程出发，利用非均匀网格有限差分法导出了二维大地电磁正演计算的线性方程组。通过对浅部不均匀体模型的模拟计算，得到了大地电磁的静位移响应规律，同时也验证了非均匀网格有限差分正演算法的准确性。

关键词 :大地电磁，二维模型，静位移，有限差分法，非均匀网格

Copyright © 2018 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

大地电磁测深(Magnetotelluric，简称MT)是以天然电磁场为场源来研究地球内部电性结构的一种重要的地球物理手段。当交变电磁场在地下介质中传播时，由于趋肤深度的作用，不同频率的信号具有不同的穿透深度，在地面上观测大地电磁场，它的频率域响应将反映着地下介质电性的垂向分布情况。因此，研究大地电磁的频率域响应，可以获得地下不同深度介质的电阻率信息 [1] 。然而，当近地表存在局部导电不均匀体时，在外电场的作用下，电流流过不均匀体时，不均匀体表面会产生积累电荷，由此产生一个与外电场成正比的附加电场，导致实测电场产生畸变，根据实测电磁场计算的视电阻率也随之产生畸变，该畸变为静态效应或静位移 [2] [3] [4] 。

大地电磁正演模拟的数值方法主要有3种：有限单元法 [5] [6] [7] 、有限差分法 [8] [9] [10] [11] 和积分方程法 [12] [13] ，前两者经常用于二维数值模拟，后者主要用在三维数值模拟。而有限差分法作为一种经典的数值模拟方法，依靠其相对简单的理论基础及简便的计算过程，在地球物理正演计算中应用非常广泛。但传统有限差分法的网格采用均匀剖分方式，必然存在着计算精度与计算速度的矛盾。若将传统的网格剖分方式应用到大地电磁正演过程中，对于衰减极快的高频电磁波，很容易造成采样的不足；而对于衰减极慢的低频电磁波，却会出现过采样的情况 [14] 。因此，本文选择非均匀网格有限差分法模拟大地电磁的静位移响应，详细推导了非均匀网格有限差分正演算法，并对浅地表不均匀体模型进行了试算分析。

2. 非均匀网格差分正演算法

2.1. 边值问题

在二维地下介质中，大地电磁场 $E_x$ 、 $H_x$ 满足的偏微分方程为 [15]

$\frac{\partial }{\partial y}\left(\frac{1}{\text{i}\omega \mu }\frac{\partial E_x}{\partial y}\right)+\frac{\partial }{\partial z}\left(\frac{1}{\text{i}\omega \mu }\frac{\partial E_x}{\partial z}\right)+\sigma E_x=0$ (1)

和

$\frac{\partial }{\partial y}\left(\frac{1}{\sigma }\frac{\partial H_x}{\partial y}\right)+\frac{\partial }{\partial z}\left(\frac{1}{\sigma }\frac{\partial H_x}{\partial z}\right)+\text{i}\omega \mu H_x=0$ (2)

式(1)和式(2)可统一表示成

$\nabla \cdot \left(\tau \nabla u\right)+\lambda u=0$ (3)

对于TE极化模式，有

$u=E_x$ ， $\tau =\frac{1}{\text{i}\omega \mu }$ ， $\lambda =\sigma$ (4)

而对于TM极化模式，有

$u=H_x$ ， $\tau =\frac{1}{\sigma }$ ， $\lambda =\text{i}\omega \mu$ (5)

为了求解大地电磁正演问题的亥姆霍兹方程，即式(3)，我们还必须给出相应的边界条件，如图1所示。

图1. 二维地电模型的边界示意图

1) TE极化模式的外边界条件

a) 上边界 $z=z_{\min }$ 离地面足够远，使异常场在边界上为零，即以该处的 $u$ 为1单位

${u|}_{z_{\min }}=1$ (6)

b) 下边界 $z=z_{\max }$ 为第三类边界条件，即

${\left(\frac{\partial u}{\partial z}+\sqrt{-\text{i}\omega \mu \sigma }u\right)|}_{z=z_{\max }}=0$ (7)

c) 取左右边界 $y=y_{\min }$ 、 $y=y_{\max }$ 离局部不均匀体足够远，电磁场在 $y_{\min }$ 、 $y_{\max }$ 上左右对称，其上的边界条件是

${\frac{\partial u}{\partial y}|}_{y=y_{\min }}={\frac{\partial u}{\partial y}|}_{y=y_{\max }}=0$ (8)

2) TM极化模式外边界条件

a) 上边界 $z=z_{\min }$ 直接取在地面上，并以该处的 $u$ 为1单位，则有

${u|}_{z=z_{\min }}=1$

2) 下边界 $z=z_{\max }$ 的边界条件，同TE极化模式，见式(7)。

3) 左右边界 $y=y_{\min }$ 和 $y=y_{\max }$ 的边界条件，同TE极化模式，见式(8)。

2.2. 差分线性方程组导出

将二维地电模型离散化，如图2所示。下面以求解TM极化模式的电磁场为例，详细推导差分正演算法。

图2. 二维地电模型离散化

对于研究区域的内部节点，在单元 $\left(i,j\right)$ 中心处，采用有限差分近似计算偏导数：

${\left[\frac{\partial }{\partial y}\left(\tau \frac{\partial u}{\partial y}\right)\right]}_{i,j}\approx \frac{1}{\Delta y}\left[{\left(\tau \frac{\partial u}{\partial y}\right)}_{i+\frac{1}{2},j}-{\left(\tau \frac{\partial u}{\partial y}\right)}_{i-\frac{1}{2},j}\right]\approx \frac{{\tau }_{i,j}+{\tau }_{i+1,j}}{2}\frac{u_{i+1,j}-u_{i,j}}{\Delta y^2}-\frac{{\tau }_{i,j}+{\tau }_{i-1,j}}{2}\frac{u_{i,j}-u_{i-1,j}}{\Delta y^2}$ (9)

和

${\left[\frac{\partial }{\partial z}\left(\tau \frac{\partial u}{\partial z}\right)\right]}_{i,j}\approx \frac{1}{\Delta z}\left[{\left(\tau \frac{\partial u}{\partial y}\right)}_{i,j+\frac{1}{2}}-{\left(\tau \frac{\partial u}{\partial y}\right)}_{i,j-\frac{1}{2}}\right]\approx \frac{{\tau }_{i,j}+{\tau }_{i,j+1}}{2}\frac{u_{i,j+1}-u_{i,j}}{\Delta z^2}-\frac{{\tau }_{i,j}+{\tau }_{i,j-1}}{2}\frac{u_{i,j}-u_{i,j-1}}{\Delta z^2}$ (10)

将式(9)和式(10)代入到式(3)，可得

${\left[\frac{\partial }{\partial y}\left(\tau \frac{\partial u}{\partial y}\right)\right]}_{i,j}+{\left[\frac{\partial }{\partial z}\left(\tau \frac{\partial u}{\partial z}\right)\right]}_{i,j}+\lambda u_{i,j}=0$

即

$\left(\frac{{\tau }_{i,j}+{\tau }_{i+1,j}}{2}\frac{u_{i+1,j}-u_{i,j}}{\Delta y^2}-\frac{{\tau }_{i,j}+{\tau }_{i-1,j}}{2}\frac{u_{i,j}-u_{i-1,j}}{\Delta y^2}\right)+\left(\frac{{\tau }_{i,j}+{\tau }_{i,j+1}}{2}\frac{u_{i,j+1}-u_{i,j}}{\Delta z^2}-\frac{{\tau }_{i,j}+{\tau }_{i,j-1}}{2}\frac{u_{i,j}-u_{i,j-1}}{\Delta z^2}\right)+\lambda u_{i,j}=0$ (11)

或写成

$\begin{array}{l}\frac{{\tau }_{i,j}+{\tau }_{i+1,j}}{2\Delta y^2}{u}_{i+1,j}+\frac{{\tau }_{i,j}+{\tau }_{i-1,j}}{2\Delta y^2}{u}_{i-1,j}+\frac{{\tau }_{i,j}+{\tau }_{i,j+1}}{2\Delta z^2}{u}_{i,j+1}+\frac{{\tau }_{i,j}+{\tau }_{i,j-1}}{2\Delta z^2}{u}_{i,j-1}\\ \left(-\frac{{\tau }_{i-1,j}+2{\tau }_{i,j}+{\tau }_{i+1,j}}{2\Delta y^2}-\frac{{\tau }_{i,j-1}+2{\tau }_{i,j}+{\tau }_{i,j+1}}{2\Delta z^2}+\lambda \right)u_{i,j}=0\end{array}$ (12)

令

$u=\left(\begin{array}{c}{u}_1\\ u_2\\ ⋮\\ u_{\left(i-1\right)\times Nz+j}\\ u_{\left(i-1\right)\times Nz+j+1}\\ ⋮\\ u_{Ny\times N_z-1}\\ u_{Ny\times N_z}\end{array}\right)=\left(\begin{array}{c}{u}_{1,1}\\ u_{1,2}\\ ⋮\\ u_{i,j}\\ u_{i,j+1}\\ ⋮\\ u_{Ny,N_z-1}\\ u_{Ny,N_z}\end{array}\right)$ (13)

推广到所有节点，将磁场所满足的4个边界条件待入，可得

$Ku=p$ (14)

这时，方程组含有 $Nz\times Ny$ 个方程以及 $Nz\times Ny$ 个未知数，该线性方程组的系数矩阵同样具有稀疏形式，如图3所示。求解上述线性方程组(14)即可得到各单元中心处的磁场值，从而可以进一步计算模型响应的视电阻率和相位。

2.3. 电磁响应计算

当计算出各单元中心的 $u$ 值后，再利用数值方法求出场值沿垂向的偏导数 $\frac{\partial u}{\partial z}$ ，它相当于 $\frac{\partial E_x}{\partial z}$ 或 $\frac{\partial H_x}{\partial z}$ ，代入到以下两式便可计算视电阻率和阻抗相位。

对于TE极化模式，有

$\begin{array}{l}{Z}_{TE}=E_x/\frac{1}{i\omega \mu }\frac{\partial E_x}{\partial z}\\ {\rho }_a^{TE}=\frac{1}{\omega \mu }{\left|Z_{TE}\right|}^2\\ {\phi }^{TE}=\arctan \frac{\mathrm{Im}\left[Z_{TE}\right]}{\mathrm{Re}\left[Z_{TE}\right]}\end{array}\}$ (15)

对于TM极化模式，有

图3. 差分法计算二维大地电磁响应形成的系数矩阵

$\begin{array}{l}{Z}_{TM}=-\frac{1}{\sigma }\frac{\partial H_x}{\partial z}/H_x\\ {\rho }_a^{TM}=\frac{1}{\omega \mu }{\left|Z_{TM}\right|}^2\\ {\phi }^{TM}=\arctan \frac{\mathrm{Im}[Z_{TM}]}{\mathrm{Re}[Z_{TM}]}\end{array}\}$ (16)

3. 静位移模型试算分析

选取的静位移模型为层状介质存在局部电性不均匀体，如图4所示。层状介质为三层A型地电模型，其模型参数为 ${\rho }_1=50\Omega \cdot \text{m}$ ， ${\rho }_2=200\Omega \cdot \text{m}$ ， ${\rho }_3=500\Omega \cdot \text{m}$ ， $h_1=450m$ 和 $h_2=1150m$ ；局部不均匀体的长和宽均为50 m，其电阻率值分别设置为 ${\rho }_0=5\Omega \cdot \text{m}$ 、 ${\rho }_0=50\Omega \cdot \text{m}$ 和 ${\rho }_0=500\Omega \cdot \text{m}$ 。

图4. 二维静位移地电模型

利用有限差分法模拟静位移模型，我们设置的二维正演网格剖分为 $50\times 40$ ，并选用40个记录频点(0.001~1000 Hz)，计算的TM视电阻率腻断面如图5所示。从模拟结果可以看出，近地表存在局部电性不均匀体时，视电阻率值受静态效应的影响严重。当近地表不均匀体为低阻时，其所在位置下方各频点视电阻率较正常值偏低，拟断面表现为下凹；当近地表不均匀体为高阻时，其所在位置下方各频点视电阻率较正常值偏高，拟断面表现为上凸。

图5. TM极化模式下高阻与低阻局部不均匀体的大地电磁响应拟断面图

为了分析局部电性不均匀体对同一测点各频点的影响情况，我们选择不均匀体上方的中心点为研究对象。图6给出了局部低阻不均匀体与局部高阻不均匀体的视电阻率响应对比，从图上可以看出：当存在局部高阻不均匀体时，视电阻率—频率曲线整体上移，电磁感应效应可以忽略不计；当存在局部低阻不均匀体时，视电阻率—频率曲线整体下移，但高频段视电阻率曲线的趋势发生了变化，这是由于低阻不均匀体的电磁感应引起的。

图6. TM极化模式下高阻与低阻局部不均匀的视电阻率–频率曲线

4. 结语

1) 针对各向同性介质中大地电磁场所满足的边值问题，详细推导了二维地电模型大地电磁场响应的有限差分正演算法；在保证计算精度的条件下，本文的非均匀网格剖分方式能减少正演计算时间，提高计算效率。

2) 通过对浅部不均匀体模型的模拟计算，得到了大地电磁的静位移响应规律，同时也验证了非均匀网格有限差分正演算法的准确性。

3) 高阻不均匀体使受影响测点的视电阻率整体向上移动而趋势不变，并且不均匀体电阻率越高，视电阻率向上平移的幅度越大。

4) 低阻不均匀体使受影响测点的视电阻率整体向下移动，并且由于电磁感应现象，单点视电阻率曲线在高频段趋势发生变化。

基金项目

水能资源利用关键技术湖南省省重点实验室开放研究基金(PKLHD201702)资助。

文章引用

邓小康,童孝忠. 基于非均匀网格有限差分法的大地电磁静位移模拟
Modeling of Magnetotelluric Static Shift Based on Non-Uniform Meshes Finite Difference Method[J]. 地球科学前沿, 2018, 08(08): 1353-1361. https://doi.org/10.12677/AG.2018.88148

参考文献