利用Chebyshev谱方法模拟一维大地电磁响应

童孝忠，何 婷，赵理芳，李爱勇

(1.中南大学 地球科学与信息物理学院，湖南 长沙 410083；2.中南大学 有色资源与地质灾害探查湖南省重点实验室，湖南 长沙 410083；3.江苏省有色金属华东地质勘查局 八一四队，江苏 镇江 212005)

1 引 言

大地电磁测深(Magnetotelluric，简称MT)是以天然电磁场为场源来研究地球内部电性结构的一种重要的地球物理手段，在深部地球动力学研究、矿产与油气资源勘探、地热资源勘探等领域得到广泛应用[1]。但大地电磁勘探过程中，需要将电磁信号转换为勘探目标的电性变化，进而达到查明断裂、构造层以及岩体的空间分布特征。这一过程离不开大地电磁勘探数据的处理与定性定量解释，其中包括数据预处理、正演模拟与反演成像等。大地电磁正反演解释需要以勘探对象的实际电性结构变化为基础，可根据实际地质资料针对性地开展一维、二维和三维正反演解释。

大地电磁测深的正演问题求解与其他地球物理勘探方法的正演问题类似，不外乎有物理模拟和数值模拟两种方法，但物理模拟对制作者要求较高，并且很难实现复杂模型的模拟，而伴随着计算机的快速发展，大数据量的存储和大型线性方程组的求解已成为可能，所以数值模拟对于复杂模型更容易实现且模拟精度较高。大地电磁测深正演模拟的数值方法主要有3种：有限差分法[2-5]、有限单元法[6-8]和积分方程法[9，10]。前两种方法常常被应用于二维地电模型的数值模拟，而第3种方法主要被应用于三维地电模型的数值计算。

随着计算机技术和计算能力的迅速发展，以谱方法和谱元素法为代表的高精度算法在电磁学、计算流体力学、地球物理学等学科的偏微分方程数值解法中展现了强大活力，得到了高度关注[11-15]。为了消除Runge现象，引入Chebyshev点来代替等距点，进而产生了Chebyshev谱方法[16]。将Chebyshev谱方法用于偏微分方程求解取得了一些研究成果，但局限于定解问题的求解，被应用于地球物理的数值模拟仍处于起步阶段[17]，有待进一步深入研究。因此，本文选择高精度的Chebyshev谱方法模拟一维大地电磁响应，详细推导正演算法，并编写Matlab计算程序。

2 Chebyshev谱方法

2.1 Chebyshev点

在等距点插值过程中，会出现Runge现象，即插值函数在区间的边界处出现震荡。为了消除Runge现象，可以引入Chebyshev点来代替等距点。Chebyshev点等于第一类Chebyshev多项式的根，在区间[-1,1]内的Chebyshev点可以定义为[18]

(1)

可以把这些Chebyshev点理解为半个单位圆上等距的点在横轴上的投影，N=8的情况如图1所示。

图1 一维网格的Chebyshev结点示意图Fig.1 A sketch of Chebyshev nodes on a 1D grid

2.2 Chebyshev求导矩阵

选取计算区间为[-1,1]，坐标被离散化为xj=-cos(jπ/N)，且j=0,1,…,N。若用x表示N+1维向量(x0,x1,…,xN)T，用u代表这些位置上的函数值组成的N+1维向量(u0,u1,…,uN)T，则可以定义DN为(N+1)×(N+1)的Chebyshev求导矩阵，并使得下式成立：

u′(x)=DNu

(2)

根据多项式插值，可以得到任意N的Chebyshev求导矩阵DN中的每个元素的表达式[16]：

(3)

(4)

(5)

(6)

这里：

(7)

由式(3)～式(6)，可以直观地给出Chebyshev求导矩阵DN的结构，如图2所示。这样，便得到了用于求一阶导数的Chebyshev求导矩阵DN。

图2 Chebyshev求导矩阵结构示意图Fig.2 A sketch of Chebyshev differentiation matrix

DN是用于求一阶导数的Chebyshev求导矩阵，要得到用于求二阶导数的Chebyshev求导矩阵，可以直接对DN进行平方处理，即

(8)

3 大地电磁正演的Chebyshev谱算法

3.1 正演方程组导出

在一维大地介质中，电场Ex所满足的微分方程为[1]

(9)

式中,σ为地下介质的电导率(即电阻率的倒数)，其单位为S/m；μ为地下介质的磁导率，其值取为4π×10-7H/m。

利用Chebyshev谱方法求解，首先需要将空间变量离散化为向量z=(z0,z1,…,zN)T，相应地，Ex(z)被离散化为向量E=(E0,E1,…,EN)T、σ(z)被离散化为向量σ=(σ0,σ1,…,σN)T。于是，根据Chebyshev求导矩阵可得：

(10)

(11)

因此，电场所满足的微分方程便能转换为代数方程形式：

(12)

令

这时，式(12)可以写成：

(13)

根据上边界条件：电场在空气中不衰减，取E0=1，则有

(14)

(N+1,:)

(15)

其中I为单位矩阵。于是有

(16)

加入边界条件后的正演方程组(16)含有N+1个方程和N+1个未知数。这有别于有限差分法或有限单元法正演所形成的系数矩阵[19-25]，具有非稀疏形式，如图3所示。求解线性方程组(16)便能得到地下介质剖分节点处的电场值，再近似计算大地电磁辅助场值，即可得到一维地电模型的视电阻率和相位。

图3 Chebyshev谱方法计算形成的系数矩阵Fig.3 The coefficient matrix formed by Chebyshev spectral method

3.2 大地电磁响应计算

(17)

采用差分法近似计算偏导数，则有

(18)

其中，L是近地表节点1与节点2间的距离。

4 程序实现与算法测试

4.1 Matlab程序实现

根据上面推导的正演算法，下面给出Chebyshev谱方法计算一维大地电磁响应的Matlab程序代码，主函数程序如下：

function [Ex,rho_a,phase]=MT1D_spectral(Length,S)

%输入参数

% Length：计算区域的深度

% S： 电导率

%输出参数

% Ex：电场

% rho_a：视电阻率

% phase：相位

mu=4e-7*pi;

a=0;

b=Length;

N = length(S)-1;

fre=logspace(-3,3,40);

for i=1:length(fre)

[D,xi] = cheb(N);

D=D/((b-a)/2);

D2 = D^2;

Omega=2*pi*fre(i);

D2=D2+sqrt(-1)*Omega*mu*diag(S);

% 上边界条件

D2(1,:)=0;

D2(1,1)=1;

% 下边界条件

k=sqrt(-sqrt(-1)*Omega*mu*S(N+1));

I=eye(N+1);

D2(N+1,:)=D(N+1,:)+k*I(N+1,:);

f = zeros(N-1,1);

f=[1;f;0];

Ex(:,i)=D2f;

x_new=(a+b)/2+xi*(b-a)/2;

Ex_g=Ex(1,i);

Hy_g=(Ex(2,i)-Ex(1,i))/(x_new(2)-x_new(1));

Hy_g=Hy_g/(sqrt(-1)*mu*Omega);

rho_a(i)=abs(Ex_g/Hy_g)^2/mu/Omega;

phase(i)=-atan(imag(Ex_g/Hy_g)/real(Ex_g/Hy_g))*180/pi;

end

4.2 算法测试

图4 均匀半空间模型中电场的Chebyshev谱方法数值解Fig.4 Chebyshev spectral solution of electric field in homogeneous half-space model

图5 二层D型地电模型Fig.5 Two-layer geo-electric model of D-type

5 一维模型试算分析

选取二层D型地电模型，其模型参数为σ1=0.01 S/m，σ2=0.1 S/m和h1=1 000 m，如图5所示。采用Chebyshev谱方法进行正演近似计算，取计算区域的长度为2 km，且计算网格按Chebyshev点分别剖分为N=200、N=100、N=50和N=20。

6 结 语

1) 从大地电磁场满足的边值问题出发，推导了大地电磁一维正演计算的Chebyshev谱算法，并利用Matlab工具编写了计算一维大地电磁响应的Chebyshev谱方法近似计算程序。

2) 通过对均匀半空间模型的模拟计算，并与解析结果作对比，验证了Chebyshev谱方法计算大地电磁响应的准确性和精确性。

3) 通过计算一维D型地电模型的大地电磁响应，说明了Chebyshev谱方法的有效性，同时分析了剖分的Chebyshev点数目对视电阻率值和相位值的影响，给出了相应的经验公式。