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

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

1 引 言

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

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

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

2 Chebyshev譜方法

2.1 Chebyshev點

在等距點插值過程中，會出現(xiàn)Runge現(xiàn)象，即插值函數(shù)在區(qū)間的邊界處出現(xiàn)震蕩。為了消除Runge現(xiàn)象，可以引入Chebyshev點來代替等距點。Chebyshev點等于第一類Chebyshev多項式的根，在區(qū)間[-1,1]內(nèi)的Chebyshev點可以定義為[18]

(1)

可以把這些Chebyshev點理解為半個單位圓上等距的點在橫軸上的投影，N=8的情況如圖1所示。

圖1 一維網(wǎng)格的Chebyshev結(jié)點示意圖Fig.1 A sketch of Chebyshev nodes on a 1D grid

2.2 Chebyshev求導(dǎo)矩陣

選取計算區(qū)間為[-1,1]，坐標被離散化為xj=-cos(jπ/N)，且j=0,1,…,N。若用x表示N+1維向量(x0,x1,…,xN)T，用u代表這些位置上的函數(shù)值組成的N+1維向量(u0,u1,…,uN)T，則可以定義DN為(N+1)×(N+1)的Chebyshev求導(dǎo)矩陣，并使得下式成立：

u′(x)=DNu

(2)

根據(jù)多項式插值，可以得到任意N的Chebyshev求導(dǎo)矩陣DN中的每個元素的表達式[16]：

(3)

(4)

(5)

(6)

這里：

(7)

由式(3)～式(6)，可以直觀地給出Chebyshev求導(dǎo)矩陣DN的結(jié)構(gòu)，如圖2所示。這樣，便得到了用于求一階導(dǎo)數(shù)的Chebyshev求導(dǎo)矩陣DN。

圖2 Chebyshev求導(dǎo)矩陣結(jié)構(gòu)示意圖Fig.2 A sketch of Chebyshev differentiation matrix

DN是用于求一階導(dǎo)數(shù)的Chebyshev求導(dǎo)矩陣，要得到用于求二階導(dǎo)數(shù)的Chebyshev求導(dǎo)矩陣，可以直接對DN進行平方處理，即

(8)

3 大地電磁正演的Chebyshev譜算法

3.1 正演方程組導(dǎo)出

在一維大地介質(zhì)中，電場Ex所滿足的微分方程為[1]

(9)

式中,σ為地下介質(zhì)的電導(dǎo)率(即電阻率的倒數(shù))，其單位為S/m；μ為地下介質(zhì)的磁導(dǎo)率，其值取為4π×10-7H/m。

利用Chebyshev譜方法求解，首先需要將空間變量離散化為向量z=(z0,z1,…,zN)T，相應(yīng)地，Ex(z)被離散化為向量E=(E0,E1,…,EN)T、σ(z)被離散化為向量σ=(σ0,σ1,…,σN)T。于是，根據(jù)Chebyshev求導(dǎo)矩陣可得：

(10)

(11)

因此，電場所滿足的微分方程便能轉(zhuǎn)換為代數(shù)方程形式：

(12)

令

這時，式(12)可以寫成：

(13)

根據(jù)上邊界條件：電場在空氣中不衰減，取E0=1，則有

(14)

(N+1,:)

(15)

其中I為單位矩陣。于是有

(16)

加入邊界條件后的正演方程組(16)含有N+1個方程和N+1個未知數(shù)。這有別于有限差分法或有限單元法正演所形成的系數(shù)矩陣[19-25]，具有非稀疏形式，如圖3所示。求解線性方程組(16)便能得到地下介質(zhì)剖分節(jié)點處的電場值，再近似計算大地電磁輔助場值，即可得到一維地電模型的視電阻率和相位。

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

3.2 大地電磁響應(yīng)計算

(17)

采用差分法近似計算偏導(dǎo)數(shù)，則有

(18)

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

4 程序?qū)崿F(xiàn)與算法測試

4.1 Matlab程序?qū)崿F(xiàn)

根據(jù)上面推導(dǎo)的正演算法，下面給出Chebyshev譜方法計算一維大地電磁響應(yīng)的Matlab程序代碼，主函數(shù)程序如下：

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

%輸入?yún)?shù)

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

% S： 電導(dǎo)率

%輸出參數(shù)

% 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譜方法數(shù)值解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型地電模型，其模型參數(shù)為σ1=0.01 S/m，σ2=0.1 S/m和h1=1 000 m，如圖5所示。采用Chebyshev譜方法進行正演近似計算，取計算區(qū)域的長度為2 km，且計算網(wǎng)格按Chebyshev點分別剖分為N=200、N=100、N=50和N=20。

6 結(jié) 語

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

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

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