王書明，李德山，胡 浩

中國地質(zhì)大學(xué)（武漢）地球物理與空間信息學(xué)院，武漢 430074

1 引 言

近年來，隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，大地電磁測(cè)深法（MT）這種20世紀(jì)50年代初期提出的地球物理勘探方法在基礎(chǔ)理論、資料處理解釋以及儀器設(shè)備等方面都取得了長(zhǎng)足的進(jìn)步.然而，MT處理解釋方法還存在許多難題沒有解決，如MT三維解釋還不完善，對(duì)于復(fù)雜三維地質(zhì)構(gòu)造的精確解釋仍處于發(fā)展探索階段.對(duì)于復(fù)雜三維地質(zhì)構(gòu)造而言，一種值得研究的解決方法是把它看作由地表局部三維構(gòu)造和地下深部三維構(gòu)造兩部分組成，地球物理勘探的目標(biāo)體一般是深部三維構(gòu)造.這種情況下，研究表層三維構(gòu)造MT響應(yīng)特征，把它看作一種地質(zhì)噪聲加以除掉，無疑可以突顯深部勘探體的響應(yīng)信號(hào)，提高三維MT勘探精度.

大地電磁勘探中，近地表局部非均勻電導(dǎo)率分布使下置區(qū)域電導(dǎo)率構(gòu)造引起的電磁響應(yīng)發(fā)生畸變，隨著大地電磁信號(hào)周期增加，局部近地表構(gòu)造的影響逐漸衰減，最終與區(qū)域電導(dǎo)率構(gòu)造產(chǎn)生的感應(yīng)響應(yīng)相比變得可以忽略［1］.區(qū)域電流經(jīng)過近地表非均勻構(gòu)造時(shí)會(huì)發(fā)生畸變，如果這種畸變不太嚴(yán)重，所觀測(cè)的磁場(chǎng)水平分量接近區(qū)域磁場(chǎng)水平分量，電磁場(chǎng)畸變主要?dú)w于電場(chǎng)畸變，這種畸變主要由區(qū)域電流經(jīng)過近地表非均勻構(gòu)造時(shí)在其邊界上形成的累積電荷引起.這類畸變不但可以在大地電磁勘探中觀測(cè)到，還可以發(fā)生在直流電阻率法（DC）和可控源電磁法（CSEM）中.

對(duì)于電場(chǎng)畸變的進(jìn)一步分析會(huì)發(fā)現(xiàn)，雖然觀測(cè)電場(chǎng)的幅度被近地表非均勻體嚴(yán)重畸變，但是電場(chǎng)和磁場(chǎng)向量之間的相位關(guān)系實(shí)質(zhì)上并未受這種畸變的太大影響，相位張量與近地表非均勻體不存在時(shí)相比差別不大，大地電磁畸變分析就是要從發(fā)生畸變的觀測(cè)數(shù)據(jù)中恢復(fù)區(qū)域信息.以往，對(duì)于近地表非均勻構(gòu)造影響的研究多集中于近地表為三維，區(qū)域是一維或二維這種構(gòu)造［2－14］.這時(shí)，區(qū)域電磁場(chǎng)水平分量之間振幅和相位關(guān)系可以用具有兩個(gè)分量成分的阻抗張量表示，而對(duì)于更為普遍、一般的三維／三維構(gòu)造（近地表非均勻構(gòu)造和區(qū)域構(gòu)造均為三維），這種畸變影響研究甚少.本文將論證，在三維／三維構(gòu)造下，可以根據(jù)觀測(cè)阻抗張量計(jì)算相位張量，消除近地表非均勻構(gòu)造畸變影響，恢復(fù)區(qū)域構(gòu)造信息.

2 大地電磁相位張量基本理論

根據(jù)大地電磁阻抗張量定義［15－18］：

利用（2）（3）式確定阻抗張量分量，表示電磁場(chǎng)分量之間具有線性關(guān)系：

阻抗張量Zxx和Zxy表達(dá)式為：

式中，1和2表示大地電磁兩種極化方式.相類似，利用式（3），阻抗張量分量Zyx和Zyy可以表示為：

這里需要注意，阻抗張量分量是復(fù)數(shù)函數(shù)，分解為實(shí)部和虛部，有以下公式：

其中，det（X）為X的矩陣行列式.

3 積分方程多網(wǎng)格法

3.1 基本原理

積分方程法（IE）是電磁數(shù)值模擬與反演的有力工具［20－21］，這種方法基于麥克斯韋方程組簡(jiǎn)化為關(guān)于異常體內(nèi)過載電流je的積分方程組.積分方程法可以有效應(yīng)用于大地電磁數(shù)值模擬，背景模型為層狀電導(dǎo)率分布σb，異常體電導(dǎo)率為σa（r）＝σb＋Δσ（r）.

這種模型中電磁場(chǎng)可以表示為背景場(chǎng)和異常場(chǎng)之和：

其中，背景場(chǎng)是由給定的場(chǎng)源在背景電導(dǎo)率分布σb中產(chǎn)生，而異常場(chǎng)由異常電導(dǎo)率分布感應(yīng)電流產(chǎn)生.由此，電磁場(chǎng)可由下列積分表達(dá)式獲得：

為了加快三維電磁計(jì)算速度、提高資料處理效率，同時(shí)保持一定解釋精度，十幾年來，人們發(fā)展了多種不同程度Born近似方法［22－23］，如擴(kuò)展Born近似法［24］、似線性近似（QL）、似解析近似（QA）［25－26］等.本文在QL基礎(chǔ)上，采用了多網(wǎng)格方法，處理方法是：利用粗網(wǎng)格和細(xì)網(wǎng)格分別剖分模型.正演數(shù)值模擬分兩步：第一步在粗網(wǎng)格基礎(chǔ)上應(yīng)用嚴(yán)格積分方程法確定電磁場(chǎng)，繼而利用該積分方程法模擬結(jié)果計(jì)算電反射張量；第二步對(duì)由粗網(wǎng)格計(jì)算的電反射張量進(jìn)行插值，在細(xì)網(wǎng)格基礎(chǔ)上利用QL方法計(jì)算電磁場(chǎng).試驗(yàn)證明，這種技術(shù)可以大大加快計(jì)算過程，同時(shí)保持?jǐn)?shù)值模擬的計(jì)算精度.

在QL方法中，我們可以把異常體電導(dǎo)率看作某種已知背景電導(dǎo)率分布上的一種擾動(dòng).這種情況下，電磁問題的解包含兩部分：（1）線性部分，非均勻體對(duì)源場(chǎng)的直接散射，不考慮散射電流之間的耦合；（2）非線性部分，由異常電導(dǎo)率和非均勻構(gòu)造中未知散射場(chǎng)聯(lián)合效應(yīng)組成.本文所用的QL方法基于這樣的假設(shè)，異常場(chǎng)與背景場(chǎng)的絕對(duì)值之間為線性比例關(guān)系［27］：

得到λ（rc）后，再引入精細(xì)網(wǎng)格Σf描述同一模型中的電導(dǎo)率分布，通過線性插值可以確定新網(wǎng)格上的λ（rc）值.利用（20）式可以計(jì)算具有精細(xì)剖分的新網(wǎng)格Σf單元中心處異常電場(chǎng)Ea（rf），其中rf表示精細(xì)剖分單元位置：

最后，對(duì)于精細(xì)網(wǎng)格可以利用麥克斯韋方程組積分方程的離散形式計(jì)算觀測(cè)數(shù)據(jù).

3.2 精度驗(yàn)證

構(gòu)建模型：背景構(gòu)造為5層水平電阻率分布，各層電導(dǎo)率和厚度分別為10、1000、5、10000、20Ωm，以及1、2、2、95km，底層為均勻半空間.在第三層背景構(gòu)造中，存在一個(gè)三維異常體，其長(zhǎng)寬均為3km，高度與第三層層高相同，異常體電阻率為2Ωm.分別利用正常網(wǎng)格積分方程法和本文實(shí)現(xiàn)的多網(wǎng)格積分方程計(jì)算，異常體正常網(wǎng)格剖分單元為x、y方向均為300m，垂直方向?yàn)?00m；粗網(wǎng)格剖分單元為x、y方向均為1000m，垂直方向?yàn)?00m.接收測(cè)線位于異常體的正上方，結(jié)果如圖1、圖2所示，二者相對(duì)誤差保持在3%之內(nèi).

4 三維模型正演計(jì)算

4.1 模型1

為了檢驗(yàn)相位張量算法效果，本文選擇一個(gè)具有代表性的區(qū)域大地電磁勘探模型，該模型背景構(gòu)造為四層電阻率分布，以及兩個(gè)立方狀異常體.圖3是異常體三維分布圖，圖中上下兩個(gè)異常體的水平尺度分別為20km×20km和30km×30km.兩個(gè)異常體水平位置分別為：上異常體分別沿x、y向自－20km開始至0km止；下異常體分別沿x、y向自－10km開始至20km止.層狀背景模型的電阻率自頂部起依次是10，1000，50，10000Ωm，其下是電阻率為20Ωm的均勻半空間，代表導(dǎo)電性的軟流圈.背景層厚度分別為1000，2000，2000，95000m.上邊的異常體由第二層隆起形成，從而在電阻率為10Ωm的頂層形成電阻率為1000Ωm、600m厚的異常.第二個(gè)異常體電阻率為2Ωm、厚度為2000m，位于第三層之中.數(shù)值模擬分別沿x、y向，自－30km開始至30km止，以間隔5km布設(shè)接收器，共169個(gè)，共接收四個(gè)電磁場(chǎng)水平分量Ex，Ey，Hx，Hy.所記錄的頻率等對(duì)數(shù)間隔分布于0.001和10Hz之間.

根據(jù)IE原則，剖分單元區(qū)域分為兩個(gè)，分別對(duì)應(yīng)上下異常體.上異常體的初始網(wǎng)格剖分單元水平尺度沿x、y方向均為2.5km，垂直方向?yàn)?0m；下異常體初始網(wǎng)格剖分單元水平尺度沿x、y方向均為2.5km，垂直方向?yàn)?50m.則上異常體被剖分為768個(gè)單元（8×8×12），下異常體被剖分為1152個(gè)單元（12×12×8）.上異常體的粗網(wǎng)格剖分單元水平尺度沿x、y方向均為5km，垂直方向?yàn)?00m；下異常體粗網(wǎng)格剖分單元水平尺度沿x、y方向均為5km，垂直方向?yàn)?00m.則上異常體被剖分為96個(gè)單元（4×4×6），下異常體被剖分為144個(gè)單元（6×6×4）.在CPU主頻為1.5GHz的個(gè)人電腦進(jìn)行數(shù)值模擬測(cè)試，IE法計(jì)算時(shí)間約為12min，多網(wǎng)格QL法的計(jì)算時(shí)間約為23s，顯著減小了數(shù)值模擬時(shí)間.對(duì)不同頻率和不同測(cè)點(diǎn)，對(duì)比兩種數(shù)值模擬結(jié)果，最大相對(duì)誤差小于3%，多網(wǎng)格法的計(jì)算精度可以接受.

圖4是阻抗張量各分量振幅.顯然，對(duì)于這種一般三維電阻率分布模型，阻抗張量分量｜Zxy｜和｜Zyx｜能夠正確反映地下兩個(gè)異常體的信息.｜Zxy｜和｜Zyx｜圖中右上區(qū)塊反映了下異常體的感應(yīng)信號(hào)，左下區(qū)塊則反映上異常體二次感應(yīng)信號(hào).圖5是相位張量各分量.可以看到，相位張量分量Pxx和Pyy也能夠正確反映地下異常體的信息.Pxx和Pyy圖中近似橢圓狀反映下異常體感應(yīng)信息，左下拉伸扭曲處則對(duì)應(yīng)上異常體.相對(duì)淺部異常體，相位張量能更好反映深部區(qū)域構(gòu)造信息.在0.001和10Hz頻率范圍內(nèi)測(cè)試了十幾個(gè)頻點(diǎn)，結(jié)果穩(wěn)定，本文僅給出0.2Hz模擬結(jié)果.

圖5 相位張量各分量（頻率：0.2Hz，單位：°）Fig.5 Each component of phase tensor（frequency：0.2Hz，unit：°）

4.2 模型2

與上述模型1相似，模型2具有相同的背景電導(dǎo)率分布，以及同樣的接收配置和頻率，圖6為模型2異常體電阻率三維分布圖.與模型1不同之處在于，模型2頂層中沒有電阻性異常體，取而代之的是一個(gè)非均勻?qū)樱呻娮杪蕿?，100和1000Ωm的單元隨機(jī)分布組成，單元尺寸為2500m×2500m×20m，該非均勻?qū)幽M經(jīng)常出現(xiàn)在勘探區(qū)域的近地表非均勻體，正是這些非均勻體引起了MT觀測(cè)數(shù)據(jù)中的靜位移.

剖分單元區(qū)域分為兩個(gè)，分別對(duì)應(yīng)上下異常體.上異常體為頂層非均勻分布，剖分單元水平尺度沿x、y方向均為500m，垂直方向?yàn)?m.與上異常體組成單元相比，其剖分單元已經(jīng)較大，不再考慮粗網(wǎng)格近似，則上異常體被剖分為57600個(gè)單元（120×120×4）.下異常體初始網(wǎng)格剖分單元水平尺度沿x、y方向均為2.5km，垂直方向?yàn)?50m，下異常體被剖分為1152個(gè)單元（12×12×8）.下異常體粗網(wǎng)格剖分單元水平尺度沿x、y方向均為5km，垂直方向?yàn)?00m，下異常體被剖分為144個(gè)單元（6×6×4）.類似于模型1，多網(wǎng)格法加快計(jì)算速度的同時(shí)，保持了可以接受的數(shù)值模擬精度.

圖7是頻率為0.2Hz時(shí)，模型中不存在地表非均勻?qū)忧闆r下阻抗張量各分量振幅.對(duì)于這種一般三維電阻率分布模型，與模型1相似，阻抗張量分量｜Zxy｜和｜Zyx｜能夠準(zhǔn)確反映地下異常體的信息.圖8是相位張量各分量.可以看到，相位張量分量Pxx和Pyy也能夠正確反映地下異常信息.圖9、圖10對(duì)應(yīng)地表存在非均勻分布模型2，這種模型下，阻抗張量各分量被地表異常層強(qiáng)烈干擾，與不存在非均勻?qū)幽Ｐ拖啾?，已?jīng)不能夠直接從｜Zxy｜和｜Zyx｜分量看到下異常信號(hào).然而，相位張量分量Pxx和Pyy仍然能夠清楚顯示下異常體感應(yīng)信號(hào)，雖然與非均勻?qū)硬淮嬖跁r(shí)的信號(hào)比較，這兩個(gè)分量信號(hào)受到表層干擾有一定的畸變.在0.001和10Hz頻率范圍內(nèi)測(cè)試了十幾個(gè)頻點(diǎn)，結(jié)果顯示，隨著頻率的增加，表層三維構(gòu)造對(duì)下置勘探體的影響越來越大，在0.5Hz以上尤為明顯，但是在大地電磁主要勘探頻段內(nèi)，相位張量分量能夠清晰地反映出深部勘探體的感應(yīng)信號(hào).限于篇幅，本文僅給出0.2Hz模擬結(jié)果.

圖6 真實(shí)異常電阻率分布三維示意圖Fig.6 3Dfigure of true anomalous resistivity distribution

5 結(jié) 論

解釋MT數(shù)據(jù)困難之一是來自深部的勘探目標(biāo)信息被近地表非均勻電導(dǎo)率分布畸變，為了克服這個(gè)困難，本文研究并實(shí)現(xiàn)了相位張量正演數(shù)值算法.模擬結(jié)果表明，與傳統(tǒng)的阻抗張量各分量振幅相比，本文實(shí)現(xiàn)的相位張量實(shí)際應(yīng)用價(jià)值在于，這種方法可以有效壓制地表異常電阻率分布的干擾，突出反映深部勘探體信息.尤其是，這種方法可以應(yīng)用于更為一般和普遍的三維／三維地質(zhì)構(gòu)造（地表非均勻構(gòu)造是三維，深部區(qū)域構(gòu)造亦為三維），而無須像傳統(tǒng)方法所假設(shè)的那樣，深部區(qū)域?yàn)橐痪S或二維構(gòu)造，從而增加了方法的實(shí)用性.為了加快計(jì)算，減小算法對(duì)計(jì)算機(jī)內(nèi)存的需求，正演算法采取了積分方程多網(wǎng)格技術(shù)，數(shù)值試驗(yàn)表明這是一種能夠保持一定精度的快速有效的數(shù)值模擬技術(shù).下一步的工作重點(diǎn)是，基于這種正演算法，實(shí)現(xiàn)三維相位張量反演技術(shù).根據(jù)相位張量正演結(jié)果初步判斷，對(duì)于地表構(gòu)造復(fù)雜或帶有地形起伏的實(shí)際勘探區(qū)域，相應(yīng)的反演技術(shù)應(yīng)該可以得到更可靠的區(qū)域構(gòu)造.

致 謝 本項(xiàng)研究的部分工作是在美國尤他大學(xué)CEMI課題組完成，在此期間得到了Michael S.Zhdanov教授的指導(dǎo)和幫助，作者表示誠摯的感謝！

（

）

［1］ Jiracek G R.Near－surface and topographic distortions in electromagnetic induction.Surv.Geophys.，1990，11（2－3）：163－203.

［2］ Bruton P.Analysis of broadband magnetotelluric data and an application to the Irish Variscides［Ph.D.thesis］.Galway：National University of Ireland，1994.

［3］ Smith J T.Understanding telluric distortion matrices.Geophys.J.Int.，1995，122（1）：219－226.

［4］ 魏勝，王家映，羅志瓊.全MT張量阻抗分解及其應(yīng)用.／／應(yīng)用地球物理學(xué)進(jìn)展.武漢：中國地質(zhì)大學(xué)出版社，1998：38－45.

Wei S，Wang J Y，Luo Z Q.Full MT impedance tensor separation and its application.／／Applied Geophysics Advance（in Chinese）.Wuhan：China University of Geosciences Press，1998：38－45.

［5］ 王書明.表面局部三維大地電磁曲線畸變校正：MT畸變校正阻抗張量分解技術(shù).西北地震學(xué)報(bào)，1998，20（4）：1－11.

Wang S M.The correction of magnetotelluric curve distortion caused by surficial local three－dimension inhomogeneities：the impedance tensor decomposition technique for the correction of MT curves distortion.NorthwesternSeismologicalJournal（in Chinese），1998，20（4）：1－11.

［6］ 晉光文，孫潔，王繼軍.大地電磁（MT）阻抗張量的正則分解及其初步應(yīng)用.地震地質(zhì)，1998，20（3）：243－249.

Jin G W，Sun J，Wang J J.Canonical decomposition of the magnetotelluric（MT）impedance tensor and its preliminary application.SeismologyandGeology（in Chinese），1998，20（3）：243－249.

［7］ 趙國澤，湯吉，劉鐵勝等.山西陽高—河北容城剖面大地電磁資料的初步解釋——阻抗張量分解技術(shù)及其應(yīng)用.地震地質(zhì)，1996，18（1）：66－74.

Zhao G Z，Tang J，Liu T S，et al.Preliminary interpretation of MT data along profile from Yanggao to Rongcheng：application of decomposition of MT impedance tensor.SeismologyandGeology（in Chinese），1996，18（1）：66－74.

［8］ Ritter P.Separation of local and regional information in geomagnetic response functions using hypothetical event analysis［Ph.D.thesis］.Edinburgh：University of Edinburgh，1996.

［9］ McNeice G W，Jones A G.Multisite，multifrequency tensor decomposition of magnetotelluric data.Geophysics，2001，66（1）：158－173.

［10］ 晉光文，孔祥儒.大地電磁阻抗張量的畸變與分解.北京：地震出版社，2006.

Jin G W，Kong X R.MT Impedance Distortion and Separation（in Chinese）.Beijing：Seismological Press，2006.

［11］ 李洋，于鵬，張羅磊等.基于混合優(yōu)化算法的MT阻抗張量畸變分解方法.地球物理學(xué)報(bào)，2010，53（8）：1924－1930.

Li Y，Yu P，Zhang L L，et al.MT distortion decomposition of magnetotelluric impedance tensor based on hybrid optimization algorithm.ChineseJ.Geophys.（in Chinese），2010，53（8）：1924－1930.

［12］ 蔡軍濤，陳小斌，趙國澤.大地電磁資料精細(xì)處理和二維反演解釋技術(shù)研究（一）——阻抗張量分解與構(gòu)造維性分析.地球物理學(xué)報(bào)，2010，53（10）：2516－2526.

Cai J T，Chen X B，Zhao G Z.Refined techniques for data processing and two－dimensional inversion in magnetotelluricⅠ：tensor decomposition and dimensionality analysis.ChineseJ.Geophys.（in Chinese），2010，53（10）：2516－2526.

［13］ 楊長(zhǎng)福，林長(zhǎng)佑，徐世浙.大地電磁GB張量分解法的改進(jìn).地球物理學(xué)報(bào)，2002，45（增刊）：356－364.

Yang C F，Lin C Y，Xu S Z.The improvement of the Groom－Baily′s tensor decomposition in magnetotellurics.ChineseJ.Geophys.（in Chinese），2002，45（Suppl.）：356－364.

［14］ 王立鳳，晉光文，孫潔等.一種簡(jiǎn)單的大地電磁阻抗張量畸變分解方法.西北地震學(xué)報(bào)，2001，23（2）：172－180.

Wang L F，Jin G W，Sun J，et al.A simple decomposition method of distortion in magnetotelluric impedance tensor.NorthwesternSeismologicalJournal（in Chinese），2001，23（2）：172－180.

［15］ Berdichevsky M N，Dmitriev V L，Keller G V.Magnetotellurics in the Context of the Theory of Ill－Posed Problems.Frank：Society of Exploration Geophysicists，2002.

［16］ Berdichevsky M N，Dmitriev V I.Models and Methods of Magnetotellurics.Berlin：Springer－Verlag，2008.

［17］ 陳樂壽，劉任，王天生.大地電磁測(cè)深資料處理與解釋.北京：石油工業(yè)出版社，1989.

Chen L S，Liu R，Wang T S.Magnetotelluric Data Processing and Interpretation（in Chinese）.Beijing：Oil Industry Press，1989.

［18］ 陳樂壽，王光鍔.大地電磁測(cè)深法.北京：地質(zhì)出版社，1990.

Chen L S，Wang G E.Magnetotelluric Sounding（in Chinese）.Beijing：Geological Publishing House，1990.

［19］ Caldwell T G，Bibby H M，Brown C.The magnetotelluric phase tensor.Geophys.J.Int.，2004，158（2）：457－469.

［20］ Weidelt P.Inversion of two－dimensional conductivity structures.PhysicsoftheEarthandPlanetaryInteriors，1975，10（3）：282－291.

［21］ Homann G W.Three－dimensional induced polarization and electromagnetic modeling.Geophysics，1975，40（2）：301－324.

［22］ Born M.Optics.Berlin：Springer，1933.

［23］ Born M，Wolf E.Principles of Optics.6th ed.New York：Pergamon Press，1980.

［24］ Habasy T M，Groom R W，Spies B R.Beyond the Born and Rytov approximations：a nonlinear approach to electromagnetic scattering.J.Geophys.Res.，1993，98（B2）：1759－1775.

［25］ Zhdanov M S，F(xiàn)ang S.Quasi－linear approximation in 3－D electromagnetic modeling.Geophysics，1996，61（3）：646－665.

［26］ Zhdanov M S.Geophysical Inverse Theory and Regularization Problems.Amsterdam：Elsevier，2002.

［27］ Fang S，Atlas B，Gao G Z，et al.Fast 3Dmodeling of borehole induction measurements in dipping and anisotropic formations using a novel approximation technique.Petrophysics，2004，45（4）：335－349.