童孝忠，何 婷，趙理芳，李愛勇

(1.中南大學 地球科學與信息物理學院，湖南 長沙 410083；2.中南大學 有色資源與地質災害探查湖南省重點實驗室，湖南 長沙 410083；3.江蘇省有色金屬華東地質勘查局 八一四隊，江蘇 鎮(zhèn)江 212005)

1 引 言

大地電磁測深(Magnetotelluric，簡稱MT)是以天然電磁場為場源來研究地球內部電性結構的一種重要的地球物理手段，在深部地球動力學研究、礦產與油氣資源勘探、地熱資源勘探等領域得到廣泛應用[1]。但大地電磁勘探過程中，需要將電磁信號轉換為勘探目標的電性變化，進而達到查明斷裂、構造層以及巖體的空間分布特征。這一過程離不開大地電磁勘探數據的處理與定性定量解釋，其中包括數據預處理、正演模擬與反演成像等。大地電磁正反演解釋需要以勘探對象的實際電性結構變化為基礎，可根據實際地質資料針對性地開展一維、二維和三維正反演解釋。

大地電磁測深的正演問題求解與其他地球物理勘探方法的正演問題類似，不外乎有物理模擬和數值模擬兩種方法，但物理模擬對制作者要求較高，并且很難實現復雜模型的模擬，而伴隨著計算機的快速發(fā)展，大數據量的存儲和大型線性方程組的求解已成為可能，所以數值模擬對于復雜模型更容易實現且模擬精度較高。大地電磁測深正演模擬的數值方法主要有3種：有限差分法[2-5]、有限單元法[6-8]和積分方程法[9，10]。前兩種方法常常被應用于二維地電模型的數值模擬，而第3種方法主要被應用于三維地電模型的數值計算。

隨著計算機技術和計算能力的迅速發(fā)展，以譜方法和譜元素法為代表的高精度算法在電磁學、計算流體力學、地球物理學等學科的偏微分方程數值解法中展現了強大活力，得到了高度關注[11-15]。為了消除Runge現象，引入Chebyshev點來代替等距點，進而產生了Chebyshev譜方法[16]。將Chebyshev譜方法用于偏微分方程求解取得了一些研究成果，但局限于定解問題的求解，被應用于地球物理的數值模擬仍處于起步階段[17]，有待進一步深入研究。因此，本文選擇高精度的Chebyshev譜方法模擬一維大地電磁響應，詳細推導正演算法，并編寫Matlab計算程序。

2 Chebyshev譜方法

2.1 Chebyshev點

在等距點插值過程中，會出現Runge現象，即插值函數在區(qū)間的邊界處出現震蕩。為了消除Runge現象，可以引入Chebyshev點來代替等距點。Chebyshev點等于第一類Chebyshev多項式的根，在區(qū)間[-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求導矩陣

選取計算區(qū)間為[-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)便能得到地下介質剖分節(jié)點處的電場值，再近似計算大地電磁輔助場值，即可得到一維地電模型的視電阻率和相位。

圖3 Chebyshev譜方法計算形成的系數矩陣Fig.3 The coefficient matrix formed by Chebyshev spectral method

3.2 大地電磁響應計算

(17)

采用差分法近似計算偏導數，則有

(18)

其中，L是近地表節(jié)點1與節(jié)點2間的距離。

4 程序實現與算法測試

4.1 Matlab程序實現

根據上面推導的正演算法，下面給出Chebyshev譜方法計算一維大地電磁響應的Matlab程序代碼，主函數程序如下：

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

%輸入參數

% Length：計算區(qū)域的深度

% 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譜方法進行正演近似計算，取計算區(qū)域的長度為2 km，且計算網格按Chebyshev點分別剖分為N=200、N=100、N=50和N=20。

6 結 語

1) 從大地電磁場滿足的邊值問題出發(fā)，推導了大地電磁一維正演計算的Chebyshev譜算法，并利用Matlab工具編寫了計算一維大地電磁響應的Chebyshev譜方法近似計算程序。

2) 通過對均勻半空間模型的模擬計算，并與解析結果作對比，驗證了Chebyshev譜方法計算大地電磁響應的準確性和精確性。

3) 通過計算一維D型地電模型的大地電磁響應，說明了Chebyshev譜方法的有效性，同時分析了剖分的Chebyshev點數目對視電阻率值和相位值的影響，給出了相應的經驗公式。