張昆， 何鈺嫻， 嚴加永

(1. 中國地質(zhì)科學院礦產(chǎn)資源研究所,國土資源部成礦作用和資源評價重點實驗室，北京　100037；

2. 國家地質(zhì)實驗測試中心，北京　100037)

大地電磁測深三維反演系統(tǒng)及應用

張昆1， 何鈺嫻2， 嚴加永1

(1. 中國地質(zhì)科學院礦產(chǎn)資源研究所,國土資源部成礦作用和資源評價重點實驗室，北京100037；

2. 國家地質(zhì)實驗測試中心，北京100037)

摘要：大地電磁測深是基于電磁感應原理，在地表采集正交的電、磁場來研究地殼乃至上地幔電性結(jié)構的地球物理方法，后期衍生出采集高頻天然場信號的音頻大地電磁測深和采集人工源信號的可控源大地電磁法。此類方法在深部構造研究、礦產(chǎn)資源勘查、油氣資源探測等多領域得到了廣泛應用。近幾年來，雖然其野外工作趨近三維勘探水平，但其采集數(shù)據(jù)的處理和三維反演研究仍然存在不足，限制了地質(zhì)解釋的準確度和效率。針對這些問題，我們開發(fā)了大地電磁測深的三維反演軟件系統(tǒng)，系統(tǒng)主要分為靜位移校正、三維反演的初始模型構建和三維反演3個模塊，有效地解決了數(shù)據(jù)處理中的靜位移問題和數(shù)據(jù)解釋中的三維反演問題，并且其中涉及的很多算法改進了現(xiàn)有技術，在計算精度、計算效率和成本節(jié)約方面有很大優(yōu)勢，具有良好的應用前景。

關鍵詞：大地電磁測深;三維反演系統(tǒng);靜位移校正;模型構建

1大地電磁測深原理

大地電磁測深法(MT)是基于電磁感應原理以天然電磁場為場源，在地表測量相互正交的電場和磁場來研究地殼和上地幔電性分布的地球物理方法。它以天然交變電磁場為場源，當交變電磁場以博得形式在地下介質(zhì)中傳播時，由于電磁感應作用，地面電磁場的觀測值將包含地下介質(zhì)電阻率分布的信息。而且，由于電磁場的集膚效應，不同周期的電磁場信號具有不同的穿透深度，因此，研究大地對天然電磁場的頻率響應，可以獲得地下不同深度介質(zhì)電阻率分布的信息。研究大地電磁場對了解地球內(nèi)部信息有著重要的作用。

在20世紀50年代初期，卡尼亞(Cagniard) 提出了在平面電磁波場中水平均勻?qū)訝罱橘|(zhì)的大地電磁法初始模型，研究了標量阻抗的計算方法，開創(chuàng)了電磁測深研究的先河。然而地球的層狀介質(zhì)模型假設與實際的地球模型相差很遠，不能很好地解決實際地質(zhì)問題。60年代中后期大地電磁法在理論上取得了一些重大的進展，考慮到了地球各向異性以及二維、三維等地下介質(zhì)的不均勻性，建立了大地電磁場張量阻抗的概念及相應的數(shù)據(jù)處理方法、技術[1]。例如，坎特韋爾(Cantwell)提出考慮多維介質(zhì)的張量阻抗估算方法，斯維福特(Swift)、 莫爾森(Morrison)、西姆斯(Sims)、沃扎夫(Vozzof) 等相繼研究了張量阻抗的性質(zhì)、分析方法和計算過程[1]。70年代以后隨著計算技術的發(fā)展，國內(nèi)外對大地電磁理論方法進行了深入研究，提出了各種一維及二維正反演方法、三維正演計算、場的影響及改正方法。并且，針對大地電磁法場源的隨機性和天然電磁場信號的不足，提出了音頻大地電磁法(AMT)和可控源大地電磁法(CSAMT)?？偟膩碚f，大地電磁測深及衍生方法基本分為兩部分，野外數(shù)據(jù)采集和室內(nèi)資料處理解釋，而數(shù)據(jù)處理解釋主要包含阻抗估算、去噪、靜位移校正、正反演計算等[2]。

大地電磁測深主要用于解決基于導電性差異的結(jié)構和構造研究，主要是[3]：研究沉積表層的層理和各項異性，確定結(jié)晶基底和變質(zhì)基底的深度，估計在研究地殼和上地幔時沉積表層的屏蔽影響； 建立固結(jié)地殼上部(高電阻)的地電模型，估計孔隙度和滲透性，分出異常導電的(特別是石墨化的)地帶； 估計在比較深的(導電良好的)地殼層位上流質(zhì)的含量； 研究海洋地殼層理； 研究海洋和大陸的軟流圈的導電性，估計局部溶化的程度； 建立全球的地幔導電率模型； 圈定礦體礦脈的范圍和走向，大體確定礦體深度，為有電磁差異的礦藏開發(fā)作前期準備；尋找地下水走向及深度。

由于MT具有工作效率高、勘探深度范圍大、垂向分辨能力、水平方向分辨能力高、地形影響小、高阻層的屏蔽作用小等特點，很快在深部構造、油氣探查、地震預報、地熱、煤田探測、固體礦產(chǎn)深部找礦、水文、工程、環(huán)境地質(zhì)、海洋地質(zhì)調(diào)查等得到廣泛應用，對它的研究也越來越深入。

2三維反演研究現(xiàn)狀

目前，大地電磁測深(MT)資料的三維正、反演問題，已成為國際地球內(nèi)部電磁感應領域研究的前沿課題。國外從20 世紀70 年代中期，就有關于三維電磁正演模擬的研究。隨著有限差分法、有限元法、積分方程法、邊界元法等應用，MT 二維、三維模擬和反演都取得了長足的發(fā)展。近年來，隨著計算機內(nèi)存和速度的倍增，在三維正演方面的研究已趨于成熟，交錯網(wǎng)格有限差分法成為主導的計算方法;隨著三維正演的發(fā)展，MT 的三維反演研究也日趨升溫，反演方法眾多，主要有共軛梯度法極大似然反演、非線性共軛梯度反演、擬線性近似反演、快速松弛反演、貝葉斯統(tǒng)計反演和人工神經(jīng)網(wǎng)絡反演等。

而隨著MT三維正演理論的發(fā)展，首先出現(xiàn)了積分方程法(Integral Equation, IE)，由于其快速的計算效率，曾經(jīng)引起許多學者的關注，使三維大地電磁模擬算法取得巨大的進步，在 90 年代初期，IE成為比較其他方法的標準。直到20世紀末，史密斯(Smith)于1996年提出了一種帶有超低頻阻抗校正的三維交錯網(wǎng)格有限差分迭代算法，用以模擬三維大地電磁場分量[4,5]，成為之后大地電磁場三維正、反演研究的重要基礎。此外，譚捍東等系統(tǒng)論述了大地電磁三維交錯網(wǎng)格有限差分數(shù)值模擬算法的實現(xiàn)過程，改進了正演算法，提出了簡潔的邊界條件[6]。徐凱軍等用積分方程法實現(xiàn)了均勻半空間三維大地電磁響應的數(shù)值模擬。

在三維反演方面，史密斯和布克(Booker)提出了快速松弛算法(Rapid Relaxation Inversion,簡稱 RRI)，通過解一個近似一維的反問題來計算在每個測點下面的電阻率擾動，把大地電磁三維反演問題轉(zhuǎn)化為一維反演問題以提高計算效率[7]，卻降低了計算精度。麥克(Mackie)和馬登(Madden)將共軛梯度算法(Conjugate Gradient Inversion，簡稱 CGI)引入大地電磁場的三維反演中,通過極大似然反演公式使模型響應和觀測數(shù)據(jù)之間不同的方差權重之和達到極小化[8]，是一種計算效率較高的算法，但是收斂速度較慢。此外，在OCCAM法基礎上，將反演問題從模型空間轉(zhuǎn)換到數(shù)據(jù)空間的反演算法[9]雖然在一定程度上提高了反演速度，但是對計算機CPU和內(nèi)存的要求很高。紐曼(Newman)首先實現(xiàn)了阻抗資料的三維非線性共軛梯度并行反演，并首先發(fā)表論文，通過線性搜索和預處理代替求解目標函數(shù)的二次導數(shù)，大大地減少了計算時間[10]，但基于網(wǎng)格節(jié)點的并行結(jié)構可能對計算精度有較大影響。譚捍東等實現(xiàn)了求最小構造的三維快速松弛反演算法[11]，并成功對日本Kayabe地區(qū)實測資料和新疆土屋銅礦床MT資料進行了反演計算。林昌洪等在麥克和馬登所做工作的基礎上，實現(xiàn)了阻抗資料的三維共軛梯度反演算法，并對日本Kayabe地區(qū)實測資料進行了反演計算[2]。

3三維反演系統(tǒng)開發(fā)

針對數(shù)據(jù)處理解釋中的靜位移校正和反演兩個部分，大地電磁測深三維反演系統(tǒng)主要包括三維初始模型構建、靜位移校正和三維反演等功能模塊。

3.1　三維靜位移校正

小的、淺部不均勻體可能會產(chǎn)生不同于MT記錄頻段內(nèi)所有頻率相對應的隨機響應[12]，這些小地質(zhì)體使得電阻率測深曲線在對數(shù)坐標上沿垂直頻率軸方向產(chǎn)生明顯的平行移動[13]，相位曲線保持不變的隨機響應被稱為 MT 的靜態(tài)偏移效應[14]。這一現(xiàn)象的產(chǎn)生主要是由于地表淺部二維和三維不均勻體邊界的積累電荷所產(chǎn)生的二次電場[15]的作用[16]。

前人對靜位移效應的校正方法的研究工作主要分為5類：空間濾波法[13]、表層不均勻體計算法、綜合資料解釋法[17]、獨立測量校正法[18]、相位改正方法[17]。但現(xiàn)有方法的不通用性以及改正結(jié)果的不確定性，并且，這些方法都是基于大地電磁場二維理論發(fā)展而來的，沒有將更為準確的三維大地電磁場理論應用于靜位移的研究中，而是僅在三維測量等實際應用上提到了三維空間概念。

因此，大地電磁測深三維反演系統(tǒng)以靜位移的三維理論和實測數(shù)據(jù)為基礎，通過定量分析大地電磁場的實測高頻數(shù)據(jù)(視電阻率和阻抗相位)，判斷是否存在靜位移，并依托任意維度的反演計算校正靜位移。該靜位移校正方法為計算機自動處理技術[19]，避免了額外的野外工作和室內(nèi)處理工作，將靜位移校正成本降為0，并且具有較好的校正效果，如圖1所示。

圖1　靜位移校正方法的理論驗證示例

由理論合成模型(圖1-a、b)模擬得到的正、反演結(jié)果(圖1-e、f、g)可見，當淺表存在不均勻體時，地表測點存在靜位移現(xiàn)象(圖1-c、d)，通過三維反演系統(tǒng)采用的靜位移校正方法處理得到的數(shù)據(jù)更為接近真實的電磁場響應，其反演結(jié)果(圖1-g)與正常數(shù)據(jù)的反演結(jié)果(圖1-e)十分相似，而未校正靜位移的數(shù)據(jù)反演結(jié)果與正常結(jié)果相差較大，對地下異常體的反映不準確。

3.2　三維初始模型構建

大地電磁場三維正、反演的計算量巨大，以往技術的反演計算消耗時間基本以“天”為計量單位，效率較低。而且，間接反演是目前使用的主流方法，此類算法需要首先建立初始模型，才能夠進行反演計算，而目前的現(xiàn)有算法對初始模型具有一定的依賴性，初始模型的構建對反演結(jié)果的準確性有一定的影響。因此，需要設計一種初始模型構建方法在提高反演計算精度的條件下，同時提高計算效率。此時低維度反演成為解決這一問題的有效手段。

大地電磁測深數(shù)據(jù)的一維或二維反演具有很高的計算效率和一定的計算精度(準確性)，因此，使用一維或二維反演結(jié)果作為三維反演的初始模型，能夠在三維反演計算之初便較好地擬合實測數(shù)據(jù)，能夠在很大程度上減少反演的迭代次數(shù)和計算時間，此外，一維或二維反演結(jié)果具備一定的準確性和合理性，能夠輔助三維反演達到更高的計算精度和準確性。

因此，大地電磁測深三維反演系統(tǒng)以大地電磁場一維和二維正反演、三維插值和“印?！崩碚摓榛A，使用一維或二維反演結(jié)果構建三維反演初始模型的框架；基于低維度反演結(jié)果，通過三維距離差值重構三維反演初始模型；通過“印?！狈ㄐ拚Ｐ?，從而構建最優(yōu)的三維反演初始模型[20]。

由圖2可見，一維反演模型的三維插值結(jié)果(圖2-a)較為粗糙，與層狀理論合成模型有一定的出入，但以其結(jié)果作為三維反演的初始模型得到的三維電阻率模型基本上呈現(xiàn)出傾斜層狀結(jié)構，電性分布與理論合成模型基本一致。此外，該大規(guī)模網(wǎng)格算例中均勻背景作為初始模型的三維反演迭代次數(shù)約50次，計算時間約8小時，但是用一維反演結(jié)果作為初始模型的三維反演迭代次數(shù)為15次，計算時間約2小時，大幅提高了計算效率。

圖2　初始模型構建及反演結(jié)果對比

3.3　三維反演

近年來，大地電磁測深資料的三維反演問題一直是國際地球內(nèi)部電性結(jié)構研究的前沿課題。大地電磁數(shù)據(jù)三維反演，由于數(shù)據(jù)量和計算量大，對計算機硬件資源的要求高，并且需要花費大量的計算時間。但國內(nèi)對大地電磁場正、反演的研究工作滯后，對三維正、反演問題的研究工作開展相對較少，國內(nèi)目前還缺乏一種高效的、能夠適用于普通計算機的且便于推廣的、準確反映地下三維電性結(jié)構的反演方法。而我國大地電磁測深的工作量較大，很多工作已經(jīng)達到三維勘探標準，如此龐大的數(shù)據(jù)量需要在較短時間內(nèi)進行反演解釋，在時間和經(jīng)濟成本上存在著很大的問題。

因此，大地電磁測深三維反演系統(tǒng)用預處理非線性共軛梯度方法構建三維MT反演問題[21, 22]，在總結(jié)前人工作成果[4, 5, 10]的基礎上，改進了現(xiàn)有的三維反演技術，提出了新的預處理因子。此外，系統(tǒng)優(yōu)化了反演算法，添加了并行計算構架，大幅減少了對計算機硬件資源的占用率和計算時間，使得大地電磁測深數(shù)據(jù)的三維反演能夠順利的從大型服務器推廣到普通計算機，乃至筆記本電腦。而計算時間由以往的“天”計量單位縮減到“分鐘”，計算效果如圖所示[23]。

圖3為安徽泥河鐵礦三維反演電阻率模型(圖3-b、c、d)，反演使用的數(shù)據(jù)為EH4采集的AMT阻抗數(shù)據(jù)，共使用測點224個，每個測點的使用頻點2~10個，網(wǎng)格設置為30×36×40，迭代次數(shù)為79次，反演計算時間約90 min，最終擬合差為8.24。可見，四層電性結(jié)構與地質(zhì)鉆孔信息完全符合，并且圈定了侵入巖隆起部位，即為賦礦部位。

圖3安徽泥河鐵礦三維反演電阻率模型[23]

Fig.3The 3D resistivity inversion model of Nihe iron deposit in Anhu[23]

1. 全新統(tǒng)； 2. 更新統(tǒng)； 3. 楊灣組； 4. 浮山組； 5. 雙廟組； 6. 正長斑巖； 7. 粗安玢巖； 8. 斷裂； 9. 研究區(qū)

4結(jié)語

(1) 大地電磁測深的三維反演系統(tǒng)雖然僅涉及大地電磁測深、音頻大地電磁測深和可控源音頻大地電磁測深勘探工作中的數(shù)據(jù)處理解釋部分，但具有較高的技術水平，在計算效率和計算精度上較以往技術有一定提高。此外，系統(tǒng)中涉及的靜位移校正方法為獨有的原創(chuàng)技術，在校正效果和節(jié)約成本方面有較大優(yōu)勢。

(2) 應用實例表明，本系統(tǒng)反演的地下電性結(jié)果與已知情況吻合較好，說明其具有較高的反演精度。當前，地調(diào)項目中設置了大量的大地電磁測深工作，但缺乏有效的反演軟件，嚴重制約了該方法的應用效果。本系統(tǒng)的成功開發(fā)，可為地調(diào)項目中的大地電磁數(shù)據(jù)的反演、解釋提供支撐，具有良好的應用前景和推廣價值。

(3)但該系統(tǒng)目前還處于試用階段，各個模塊還沒有有機的整合起來，目前仍然缺少一個優(yōu)化的人機交互界面控制系統(tǒng)的操作，需要進一步完善圖形界面，以方便地質(zhì)工作者使用。

參考文獻：

[1]嚴家斌.大地電磁信號處理理論及方法研究[D].湖南長沙：中南大學，2003.

[2]林昌洪.大地電磁張量阻抗三維共軛梯度反演研究[D].北京：中國地質(zhì)大學(北京), 2009.

[3]石應駿. 大地電磁測深法教程[M].北京：地震出版社,1985.

[4]Smith J T.Conservative modeling of 3-D electromagnetic fields, Part I: Properties and error analysis[J]. Geophysics, 1996, 61: 1308-1318.

[5]Smith J T.Conservative modeling of 3-D electromagnetic fields, Part II: Biconjugate gradient solution and an accelerator[J]. Geophysics, 1996,61: 1319-1324.

[6]譚捍東,余欽范,Booker J,等.大地電磁法三維交錯采樣有限差分數(shù)值模擬[J]. 地球物理學報, 2003, 46: 5.

[7]Smith J T,Booker J R. Rapid inversion of two- and three- dimensional magnetotelluric data[J]. Journal of Geophys Research, 1991, 96: 3905-3922.

[8]Mackie R L,Madden T R. Three-dimensional magnetotelluric inversion using conjugate gradients[J].Geophysics Jonrnal International, 1993, 115: 215-229.

[9]Siripunvaraporna W,Egbert G. Three-dimensional magnetotelluric inversion: data-space method[J]. Physics of the Earth and Planetary Interiors, 2005, 150: 3-14.

[10]Newman G A,Alumbaugh D L.Three-dimensional magnetotellutric inversion using non-linear conjugate gradients[J]. Geophys J Int, 2000,410:424.

[11]譚捍東,余欽范,魏文博. 大地電磁法三維快速松弛反演[J]. 地球物理學報, 2003, 46(6)：850-855.

[12]Swift C M. A magnetotelluric investigation of an electrical conductivity anomaly in the southwestern United States[D].cambridge:M I T, 1967:48-67.

[13]Berdichevsky M N,Vanyan L L,Kuznetsov V A,et al.Geoelectric model of the Baikal region[J]. Physics of the Earth and Planetary Interiors, 1980, 22: 1-11.

[14]Pellerin L,Hohmann G W. Transient electromagnetic inversion :A remedy for magnetotelluric static shifts[J]. Geophysics, 1990, 55: 1242-1250.

[15]Torres-Verdin C,Bostick F X.Principles of spatial surface electric field filtering in magnetotellurics[J]. Geophysics, 1992, 57: 603-622.

[16]Park S K,Orange A S,Madden T R.Effects of three-dimensional structure on magnetotelluric sounding curves[J]. Geophysics, 1983, 48: 1402-1405.

[17]Berdichevsky M N, Dmitriev V I. Distortions of magnetic and electrical fields by near-surface lateral inhomogeneities[J]. Acta Geodaet Geophys, 1976, 11: 447-483.

[18]Andrieux P,Wightman W E. The so-called static correction in magnetotelluric measurements[M]. Georgia,Atlanta,United States：54th Annual International Meeting SEG, 1984: 43-44.

[19]張昆.用于大地電磁場靜位移校正的模型設計軟件[P].中國人民共和國國家版權局, 2014:2014SR009601.

[20]張昆.基于一維反演的大地電磁場三維反演初始模型設計軟件[P].中國人民共和國國家版權局, 2013:2013SR160125.

[21]張昆.大地電磁場非線性共軛梯度三維反演軟件 [P].中國人民共和國國家版權局, 2012: 2012SR057553.

[22]張昆，董浩，嚴加永，等.一種并行的大地電磁場非線性共軛梯度三維反演方法[J].地球物理學報, 2013, 56(11)：3902-3911.

[23]張昆，嚴加永，呂慶田，等. 安徽泥河玢巖鐵礦電磁法探測試驗[J]. 地質(zhì)學報, 2014, 88(4)：496-506.

3D Inversion System of Magnetotelluric Sounding and its Application

ZHANG Kun1, HE Yu-xian2, YAN Jia-yong1

(1.InstituteofMineralResources,ChineseAcademyofGeologicalSciences,MLRKeyLaboratoryofMetallogenyandMineral

Assessment,Beijing100037,China; 2.NationalResearchCenterforGeoanalysis,ChineseAcademyofGeological

Sciences,Beijing100037,China)

Abstract:Magnetotelluric is a geophysics exploration method for studying the crust and upper mantle electrical structure, which measures the surface orthogonal electric and magnetic fields based on the electromagnetic induction and natural electromagnetic fields. Thereafter, the audio-frequency magnetotelluric method which collects the high frequency domain of natural field and controlled source audio-freguency magnetotelluric method which collects artificial source are derived. These methods are widely used in deep structure study, mineral resources exploration, oil and gas exploration, and so on. Though the field work achieves 3D level in recent years, 3D data processing and inversion are insufficient and the efficiency and accuracy of interpretation is limited. To solve these problems, we develop the 3D inversion system of magnetotelluric sounding which contains static shift correction, initial model design and 3D inversion. The system effectively solves the static shift of data process and 3D inversion problem. A lot of algorithms in the system are improved which have advantages on calculation accuracy, efficiency and cost savings. And the system has a good application prospect.

Key words:magnetotelluric sounding; 3D inversion system; static shift correction; model design

作者簡介：張昆(1983—)，博士，助理研究員，主要從事電磁測深方法研究。Email:zhangkun1010@163.com。

基金項目：中央級公益性行業(yè)專項“大地電磁場靜位移效應機理分析與校正方法探索(編號：K1318)”和中國地質(zhì)調(diào)查“破碎蝕變巖型金礦深部找礦勘查技術集成(編號：12120114053301)”項目聯(lián)合資助。

收稿日期：2014-10-08； 改回日期： 2014-12-05。

中圖分類號：P631

文獻標志碼：A

文章編號：2095-8706(2015)01-0025-06