大地電磁三維交錯網(wǎng)格有限差分?jǐn)?shù)值模擬的并行計算研究

李 焱，胡祥云，楊文采，魏文博，方 慧，韓 波，彭榮華

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

2 中國國土資源航空物探遙感中心，北京 100083

3 中國地質(zhì)科學(xué)院，北京 100037

4 中國地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院，北京 100083

1 引 言

大地電磁三維數(shù)值模擬一直是國際地球內(nèi)部電磁感應(yīng)領(lǐng)域研究的前沿和熱點(diǎn)課題.目前，大地電磁二維反演解譯已用于生產(chǎn)實踐［1－3］，三維反演仍處于研究和試驗階段［4－6］，實際資料的三維反演尚未大范圍推廣應(yīng)用.主要原因是三維反演計算規(guī)模很大，耗費(fèi)時間很長，普通的微機(jī)難以承受，而反演中絕大部分計算時間花費(fèi)在正演或者與正演有關(guān)的運(yùn)算上，包括模型修正量、雅可比矩陣或其與向量的乘積的計算等，因此采用并行計算來加快正演計算速度無疑具有現(xiàn)實意義.MPI（Message Passing Interface）是目前國內(nèi)外在高性能計算機(jī)系統(tǒng)中最廣泛使用的并行編程環(huán)境，它具有移植性好、功能強(qiáng)大、效率高、有多種不同的免費(fèi)、高效、實用的實現(xiàn)版本、幾乎所有的并行計算機(jī)廠商都提供對它的支持等多種優(yōu)點(diǎn)，是目前最重要的并行編程工具［7－8］.

自20世紀(jì)70年代中期開始，已經(jīng)有學(xué)者致力于三維大地電磁的正演研究［9］，目前成熟的三維正演方法主要有：有限元法［10－12］、有限差分法［13－15］、積分方程法［16－18］、有限體積法、邊界元法等，其中Mackie等發(fā)展的交錯網(wǎng)格有限差分法因迭代收斂穩(wěn)定、計算精度高、能反映較復(fù)雜地電模型已成為主導(dǎo)的正演計算方法，眾多三維反演方法均以該方法為正演核心［4，6，19－24］.但三維正演數(shù)值模擬需要求解大型稀疏對稱系數(shù)矩陣線性方程組，計算量非常大，消耗的時間非常長，尤其是網(wǎng)格剖分較密、頻點(diǎn)數(shù)較多時.基于此，本文開展了大地電磁三維交錯網(wǎng)格有限差分?jǐn)?shù)值模擬并行計算研究.大地電磁三維正演數(shù)值模擬是按照不同頻率來計算的，各頻率之間求取電磁場值的過程是相互獨(dú)立的.每個頻率都需要求解大型線性方程組一次，當(dāng)參加計算的頻率個數(shù)較多時，可以將多個頻率的計算任務(wù)平均分配到不同計算節(jié)點(diǎn)去并行執(zhí)行.基于該并行思想，本文在曙光TC－5000A高性能并行平臺上，以MPI為并行編程環(huán)境，實現(xiàn)了大地電磁三維交錯網(wǎng)格有限差分?jǐn)?shù)值模擬的并行算法.通過對兩個理論模型進(jìn)行試算和分析，證明了該并行算法的正確性和高效性，為進(jìn)一步大地電磁三維反演并行算法研究奠定了重要基礎(chǔ).

2 MPI和高性能并行平臺簡介

MPI是由全世界工業(yè)、科研和政府部門聯(lián)合建立的一個消息傳遞編程標(biāo)準(zhǔn)，其目的是為基于消息傳遞的并行程序設(shè)計提供一個高效、可擴(kuò)展、統(tǒng)一的編程環(huán)境.它是目前最為通用的并行編程方式，也是分布式并行系統(tǒng)的主要編程環(huán)境.MPI標(biāo)準(zhǔn)中定義了一組函數(shù)接口用于進(jìn)程間的消息傳遞.編程人員通過調(diào)用MPI的庫程序來編寫并行程序.MPI是一種標(biāo)準(zhǔn)或規(guī)范的代表，而不特指某一個對它的具體實現(xiàn)，一個正確的MPI程序，可以不加修改地在所有的并行機(jī)上運(yùn)行［7－8］.

MPI并行程序可分為兩種基本模式，即對等模式和主從模式.對等模式是指MPI程序的各個進(jìn)程的功能、地位相同或相近，MPI程序的代碼也應(yīng)該是相近的，所不同的只是處理的對象和操作的數(shù)據(jù).主從模式是指MPI程序的各個進(jìn)程所起的作用和地位并不相同，一個或者一些進(jìn)程完成一類任務(wù)，而另外的進(jìn)程完成其它的任務(wù)，這些功能或者地位不同的進(jìn)程所對應(yīng)的代碼也有較大的差別.

本文開發(fā)的大地電磁三維交錯網(wǎng)格有限差分?jǐn)?shù)值模擬并行程序均在中國地質(zhì)大學(xué)（武漢）曙光TC5000A高性能計算平臺上運(yùn)行.該平臺采用X86刀片集群服務(wù)器架構(gòu)，整套系統(tǒng)包括92個計算節(jié)點(diǎn)、2臺SMP 8路計算節(jié)點(diǎn)、5個I／O、雙作業(yè)調(diào)度系統(tǒng)和一臺集群管理維護(hù)節(jié)點(diǎn).節(jié)點(diǎn)間的通信連接采用20G的Infiniband連接，管理網(wǎng)絡(luò)采用1000M以太網(wǎng)交換機(jī)連接，MPI層消息傳遞延遲小于1.5μs.SMP服務(wù)器配置128G的海量內(nèi)存，每一個計算節(jié)點(diǎn)配置內(nèi)存為16G，管理節(jié)點(diǎn)和作業(yè)調(diào)度節(jié)點(diǎn)配置16G內(nèi)存、每一個I／O節(jié)點(diǎn)配置32G內(nèi)存；整個系統(tǒng)內(nèi)存總?cè)萘渴?936G.采用業(yè)界主流的x4DDR Infiniband作為通信網(wǎng)絡(luò)互聯(lián)全部節(jié)點(diǎn)，點(diǎn)對點(diǎn)單向網(wǎng)絡(luò)帶寬達(dá)到線速20Gb／s；提供適用于AMD多核平臺的全套編譯、調(diào)試軟件以及數(shù)學(xué)函數(shù)庫，支持標(biāo)準(zhǔn)的Fortran／C／C＋＋編程，支持OpenMP、MPI以及OpenMP和MPI的混合并行編程.曙光TC5000A高性能計算平臺功能強(qiáng)大、運(yùn)算快捷、存儲量大，完全滿足本文的計算要求.

3 大地電磁三維交錯網(wǎng)格有限差分?jǐn)?shù)值模擬

在大地電磁研究的頻率范圍內(nèi)，位移電流的作用可以忽略［25］.通常取電磁場隨時間變化的因子為e－iωt，麥克斯韋方程組的積分形式如下：

其中，σ為電導(dǎo)率，μ為空氣中的磁導(dǎo)率，ω為角頻率，J為傳導(dǎo)電流密度，要解上述方程組需要將連續(xù)形式的積分方程組轉(zhuǎn)化成離散形式.對研究區(qū)域離散化，即把研究區(qū)域剖分成若干個小的長方體單元.假設(shè)沿x、y和z軸方向分別剖分成Nx、Ny和Nz段，網(wǎng)格間距分別為Δx（i）（i＝1，…，Nx）、Δy（j）（j＝1，…，Ny）和Δz（k）（k＝1，…，Nz）.長方體網(wǎng)格單元的6個電磁場分量的采樣點(diǎn)位置取法為磁場H取在長方體單元棱上的中點(diǎn)處，電場E取在長方體單元面上的中心處，這樣可以自動保證電磁場分布滿足能量守恒定律.以編號為（i，j，k）的長方體網(wǎng)格單元為例，如圖1所示.

將積分形式的麥克斯韋方程組（1）、（2）離散化后即可獲得大型稀疏線性方程組：

圖1 交錯采樣網(wǎng)格示意圖Fig.1 Sketch of Staggered-grid

其中，b是包含源場和與已知邊界條件值有關(guān)的項，x是求解域內(nèi)部的未知磁場分量，A是對稱的大型稀疏系數(shù)矩陣.為了保證迭代收斂穩(wěn)定，對A進(jìn)行不完全Cholesky分解，采用雙共軛梯度法對方程（3）求解，就可以得到各網(wǎng)格單元采樣點(diǎn)處所有的磁場值H，進(jìn)而求得電場E.

大地電磁所用場源有兩種極化模式：Ex－Hy和Ey－Hx，Ex－Hy在水平面內(nèi)產(chǎn)生的電場和磁場分量值分別記為EX1、EY1、HX1和 HY1，Ey－Hx在水平面內(nèi)產(chǎn)生的電場和磁場分量值分別記為EX2、EY2、HX2和HY2，通過電磁場和阻抗張量關(guān)系式可以求得三維介質(zhì)的張量阻抗：

式（5）定義的響應(yīng)稱為XY模式響應(yīng)，（6）定義的響應(yīng)稱為YX模式響應(yīng).按照下面公式可以求出三維介質(zhì)的視電阻率和相位［25］：

4 三維正演數(shù)值模擬并行算法的總體設(shè)計

通過對三維交錯網(wǎng)格有限差分?jǐn)?shù)值模擬算法深入分析發(fā)現(xiàn)，三維正演計算中，頻率循環(huán)部分求解電磁場值為主要計算時間，占據(jù)了整個計算時間的90%以上，是需要并行化的部分.三維正演首先將積分形式的麥克斯韋方程組離散化，然后將方程組化簡去參加入邊界條件值，最后求解（3）式.對于給定的不同頻率值，方程組求出的磁場值是相互獨(dú)立的，不同頻率值相互之間對方程組沒有影響，根據(jù)這一特點(diǎn)我們可以將串行算法按頻率進(jìn)行粒度劃分，將每個頻率對應(yīng)的部分分配到不同節(jié)點(diǎn)上同時進(jìn)行計算，并行執(zhí)行.

程序采用主從并行模式，分為主進(jìn)程和子進(jìn)程，主進(jìn)程負(fù)責(zé)維護(hù)全局的數(shù)據(jù)結(jié)構(gòu)和任務(wù)的分配、參數(shù)的發(fā)送、計算結(jié)果的回收整合以及最后結(jié)果的輸出，子進(jìn)程從主進(jìn)程處接受參數(shù)執(zhí)行分配任務(wù)的計算并將結(jié)果發(fā)回主進(jìn)程.由于主進(jìn)程計算量不大，為了充分利用計算資源，不設(shè)置專用的控制節(jié)點(diǎn)，主進(jìn)程同時也作為子進(jìn)程參與計算.三維正演并行計算的基本思路是：啟動并行環(huán)境，主進(jìn)程讀入所有頻率值和模型網(wǎng)格剖分?jǐn)?shù)據(jù)，將其廣播給所有子進(jìn)程，主從進(jìn)程按照分配的頻率同時各自獨(dú)立的計算模型響應(yīng)，待全部計算完后，子進(jìn)程將計算結(jié)果發(fā)回主進(jìn)程，主進(jìn)程整合所有收到的結(jié)果，并將其輸出，結(jié)束并行環(huán)境.圖2為大地電磁三維正演并行計算簡化流程圖.

5 三維正演并行程序測試和計算效率分析

為了對開發(fā)的大地電磁三維正演并行算法的高效性進(jìn)行驗證，我們設(shè)計了兩個模型在曙光TC5000A高性能計算平臺上進(jìn)行試算，兩模型頻點(diǎn)數(shù)均為36，頻率從320Hz到0.005Hz.

5.1 模型一：3D／2D模型

如圖3所示：該模型為3D／2D模型，即區(qū)域構(gòu)造是2D，含有3D異常體.棱柱體大小為2km×2km×1km，電阻率為5Ωm，頂面埋深為1km，中心位于（0，－3000m，1500m），覆存于電阻率為100Ωm的二維構(gòu)造中，二維走向為y方向，基底電阻率為1000Ωm.三維交錯網(wǎng)格有限差分正演在x、y、z方向剖分網(wǎng)格單元數(shù)分別為Nx＝62，Ny＝56，Nz＝43.

圖2 大地電磁三維正演并行計算流程圖Fig.2 Flow chart of 3DMT forward modeling parallel computation

圖3 模型在y＝－1～1km間垂直方向斷面圖Fig.3 Vertical section of model at y＝－1～1km

圖4和圖5是3D／2D模型視電阻率和相位4種模式在剖面y＝0時沿x軸方向的響應(yīng)擬斷面圖.在XY模式和YX模式視電阻率擬斷面圖上，3D／2D模型的YX模式較好的反映出三維棱柱體的形態(tài)和范圍，特別是埋深和底界面，而XY模式反映的三維棱柱體明顯向下拉伸.XY模式和YX模式相位擬斷面圖都較好的反映出三維棱柱體的形態(tài)和范圍，其中XY模式橫向分辨率比YX模式橫向分辨率高，兩種模式對2D構(gòu)造及基底形態(tài)反映不明顯.從XX模式和YY模式視電阻率擬斷面圖上可以看出兩種模式的視電阻率值都非常小，不能反映模型的形態(tài).因為XX模式和YY模式視電阻率值都非常小，再加上正演模擬過程中解方程時所給定的收斂判別標(biāo)準(zhǔn)誤差，因此XX模式和YY模式的相位信息是不可靠的，也無法反映模型的形態(tài).從這里還可以看出阻抗信息非對角元素的作用遠(yuǎn)大于對角元素的作用.

5.2 模型二：低阻體和高阻體組合模型

設(shè)計的第二個模型如圖6所示：高阻三維棱柱體電阻率為1000Ωm，低阻三維棱柱體電阻率為10Ωm，它們大小相同，均為2km×2km×1km，頂面埋深為0.5km，中心分別為（3km，0，1km）和（－3km，0，1km）圍巖電阻率為100Ωm.三維交錯網(wǎng)格有限差分正演在x、y、z方向剖分網(wǎng)格單元數(shù)分別為Nx＝64，Ny＝44，Nz＝41.

圖6 低阻體和高阻體組合模型（a）模型在y＝－1km～1km處垂直方向斷面圖；（b）模型在z＝0.5km～1.5km間平面俯視圖.Fig.6 Low resistance and high resistance（a）Vertical section of model at y＝－1～1km；（b）Flat section of model at z＝0.5～1.5km.

圖7和圖8是高阻體和低阻體組合模型的視電阻率和相位兩種模式響應(yīng)擬斷面圖.由于XX模式和YY模式的信息不能較好的反映出異常體，這里只給出XY模式和YX模式的響應(yīng)擬斷面圖.在XY模式和YX模式視電阻率擬斷面圖上，較好的反映低阻體和高阻體，相比較而言右邊的高阻體響應(yīng)不如左邊的低阻體明顯.兩種模式的視電阻率擬斷面圖都較好的反映了異常體的頂深和橫向范圍，其中XY模式比YX模式在橫向范圍上反映更準(zhǔn)確一些，YX模式在橫向上有一定延伸.兩種模式的相位響應(yīng)擬斷面圖均能較好的反映出低阻體和高阻體的整個范圍，尤其是能反映視電阻率無法反映的底界面.

5.3 并行程序正確性的驗證

由于本文的并行計算完全是在原有串行算法上實現(xiàn)的，只是在相應(yīng)處增加了并行處理部分，即MPI處理語句和函數(shù)調(diào)用，對算法本身未做任何修改，因此，計算結(jié)果應(yīng)該和串行計算結(jié)果一致，這是檢驗并行程序設(shè)計正確性的標(biāo)準(zhǔn).圖9是3D／2D模型中串行計算和并行計算在異常體上方一測點(diǎn)上的XY模式視電阻率頻率圖，從圖中可以看出兩曲線完全重合，從而驗證了并行程序是正確的.

5.4 并行效率分析

為了檢驗開發(fā)的三維正演并行程序的效率，采取不同的節(jié)點(diǎn)數(shù)來計算模型一和模型二，并通過并行加速比和并行效率來評價并行算法的效果.并行加速比＝單機(jī)串行算法的執(zhí)行時間／N個進(jìn)程并行算法的執(zhí)行時間；并行效率＝并行加速比／參加并行計算的進(jìn)程個數(shù).表1是3D／2D模型正演測試時間統(tǒng)計表；表2是低阻體高阻體組合模型正演測試時間統(tǒng)計表.

表1 3D／2D模型正演測試時間統(tǒng)計表Table 1 Statistical runtime data for 3D／2Dforward modeling

表2 低阻高阻組合模型模型正演測試時間統(tǒng)計表Table 2 Statistical runtime data for combination model of low resistance and high resistance forward modeling

圖4 剖面y＝0時視電阻率響應(yīng)擬斷面圖（a）ρXX；（b）ρXY；（c）ρYX；（d）ρYY.Fig.4 Pseudo－section of apparent resistivity at y＝0

圖5 剖面y＝0時相位響應(yīng)擬斷面圖（a）φXX；（b）φXY；（c）φYX；（d）φYY.Fig.5 Pseudo－section of tensor phase at y＝0

圖7 剖面y＝0時視電阻率響應(yīng)擬斷面圖（a）ρXY；（b）ρYX.Fig.7 Pseudo－section of apparent resistivity at y＝0

圖8 剖面y＝0時相位響應(yīng)擬斷面圖（a）φXY；（b）φYX.Fig.8 Pseudo－section of tensor phase at y＝0

圖9 串行正演與并行正演結(jié)果對比圖Fig.9 The comparison of serial and parallel forward modeling

對表1和表2進(jìn)行分析，當(dāng)使用2個節(jié)點(diǎn)和4個節(jié)點(diǎn)時，并行效率很高.使用6個節(jié)點(diǎn)和9個節(jié)點(diǎn)時加速比雖然增大了，但增加的幅度不大，并行效率沒有2個和4個時高，這是因為隨著節(jié)點(diǎn)的增加，用于節(jié)點(diǎn)間的通信所占的比例逐漸增大，并行效率反而下降.因此如何減少節(jié)點(diǎn)間通信量是并行程序開發(fā)的關(guān)鍵.另外，一個節(jié)點(diǎn)的并行計算時間要比同一節(jié)點(diǎn)的串行計算時間要長，這是因為節(jié)點(diǎn)內(nèi)有少量的通信以及一些管理上的開銷.從表中總的可以看出并行計算的效果還是很明顯的，尤其是當(dāng)計算量很大時，并行優(yōu)越性更能體現(xiàn).

6 結(jié) 語

大地電磁三維反演推廣應(yīng)用的最大障礙是計算規(guī)模巨大，耗費(fèi)時間很長，普通微機(jī)往往難以承受.而反演中絕大部分計算時間花費(fèi)在正演或者與正演有關(guān)的運(yùn)算上，因此加快正演計算速度是提高反演計算速度的關(guān)鍵，并行計算為有效解決三維大地電磁計算問題提供了一個最優(yōu)途徑.本文在曙光TC5000A高性能計算平臺上，以MPI為并行編程環(huán)境，實現(xiàn)了大地電磁三維交錯網(wǎng)格有限差分?jǐn)?shù)值模擬并行算法.通過理論模型試算和分析，結(jié)果證明了該并行算法的正確性和高效性，很好的解決了三維正演計算速度問題，為進(jìn)一步的大地電磁三維反演并行算法研究奠定了重要基礎(chǔ).

（References）

［1］胡祥云，胡祖志，張榮等.油氣MT勘探COPROD－2S1模型數(shù)據(jù)的二維反演.天然氣工業(yè)，2004，25（9）：33－37.Hu X Y，Hu Z Z，Zhang R，et al.Two dimensional inversion of COPROD－2S1modeling dataset in oil and gas magnetotelluric exploration.Natural Gas Industry （in Chinese），2004，25（9）：33－37.

［2］譚捍東，魏文博，Unsworth M等.西藏高原南部雅魯臧布江縫合帶地區(qū)地殼電性結(jié)構(gòu)研究.地球物理學(xué)報，2004，47（4）：685－690.Tan H D，Wei W B，Unsworth M，et al.Crustal electrical conductivity structure beneath the Yarlung Zangbo Jiang structure in the southern Xizang plateau.Chinese J.Geophys.（in Chinese），2004，47（4）：685－690.

［3］葉益信，胡祥云，金鋼燮等.大地電磁二維陡邊界反演應(yīng)用效果分析.地球物理學(xué)進(jìn)展，2009，24（1）：668－674.Ye Y X，Hu X Y，Jing G X，et al.Application analysis of sharp boundary inversion of magnetotelluric data for 2D structure.Progress in Geophys （in Chinese），2009，24（1）：668－674.

［4］胡祖志，胡祥云，何展翔.大地電磁非線性共軛梯度擬三維反演.地球物理學(xué)報，2006，49（4）：1226－1234.Hu Z Z，Hu X Y，He Z Z.Pseudo－Three－Dimensional magnetotelluric inversion using nonlinear conjugate gradients.Chinese J.Geophys.（in Chinese），2006，49（4）：1226－1234.

［5］楊迪琨，胡祥云.含噪聲數(shù)據(jù)反演的概率描述.地球物理學(xué)報，2008，51（3）：901－907.Yang D K，Hu X Y.Inversion of noisy data by probabilist methodology.Chinese J.Geophys.（in Chinese），2008，51（3）：901－907.

［6］譚捍東，余欽范，Booker J等.大地電磁法三維快速松弛反演.地球物理學(xué)報，2003，46（6）：850－854.Tan H D，Yu Q F，Booker J，et al.Three－Dimensional rapid relaxation inversion for the magnetotelluric method.Chinese J.Geophys.（in Chinese），2003，46（6）：850－854.

［7］張林波，遲學(xué)斌，莫則堯等.并行計算導(dǎo)論.北京：清華大學(xué)出版社，2006.Zhang L B，Chi X B，Mo Z Y，et al.Introduction to Parallel Computing.Beijing：Tsinghua University Press，2006.

［8］都志輝 編著.高性能計算并行編程技術(shù)－MPI并行程序設(shè)計.北京：清華大學(xué)出版社，2001.Du Z H.High－Pwered Computing Parallel Programming Technique－MPI Parallel Program Design （in Chinese）.Beijing：Tsinghua University Press，2001.

［9］Hohmann G W.There－dimensional induced polarization and electromagnetic modeling.Geophysics，1975，40（2）：309－324.

［10］Rodi W L.A technique for improving the accuracy of finite element solutions for magnetotelluric data.Geophys.J.Roy.Astr.Soc.，1976，44（2）：483－506.

［11］Wannamaker P E，Stodt J A，Rijo L.Two－dimensional topographic responses in magnetotelluric modeled using finite elements.Geophysics，1986，51（11）：2131－2144.

［12］Mitsuhata Y，Uchida T.3Dmagnetotelluric modeling using the T－Ωfinite－element method.Geophysics，2004，69（1）：108－119.

［13］Mackie R L，Smith J T，Madden T R.There－dimensional electromagnetic modeling using finite difference equations：The magnetotelluric example.Radio Science，1994，29（4）：923－935.

［14］Smith J T.Conservative modeling of 3－D electromagnetic fields，Part I：Properties and error analysis.Geophysics，1996，61（5）：1308－1318.

［15］Smith J T.Conservative modeling of 3－D electromagnetic fields， Part II： Biconjugate gradient solution and an accelerator.Geophysics，1996，61（5）：1319－1324.

［16］Wannamaker P E， Hohmann G W，SanFilipo W A.Electromagnetic modeling of three－dimensional bodies in layered earths using integral equations.Geophysics，1984，49（1）：60－74.

［17］Wannamaker P E.Advances in three－dimensional magnetotelluric modeling using integral equations.Geophysics，1991，56（11）：1716－1728.

［18］Xiong Z H.Electromagnetic modeling of three－dimension structures by the method of system iteration using integral equations.Geophysics，1992，57（12）：1556－1561.

［19］Siripunvaraporn W，Egbert G.An efficient data－subspace inversion method for 2－D magnetotelluric data.Geophysics，2000，65（3）：791－803.

［20］Siripunvaraporn W， Uyeshima M， Egbert G. Threedimensional inversion for Network－Magnetotelluric data.Earth Planets Space，2004，56（9）：893－902.

［21］Siripunvaraporn W，Egbert G，Lenbury Y，et al.Threedimensional magnetotelluric inversion：data－space method.Physics of the Earth and Planetary Interiors，2005，150（1－3）：3－14.

［22］Siripunvaraporn W，Egbert G. WSINV3DMT： Vertical magnetic field transfer function inversion and parallel implementation.Physics of the Earth and Planetary Interiors，2009，173（3－4）：317－329.

［23］胡祖志，胡祥云.大地電磁三維反演方法綜述.地球物理學(xué)進(jìn)展，2005，20（1）：214－220.Hu Z Z， Hu X Y. Review of there dimensional magnetotelluric inversion methods.Progress in Geophys （in Chinese），2005，20（1）：214－220.

［24］Lin C H，Tan H D，Tong T.Three－dimensional conjugate gradient inversion of magnetotelluric full information data.Applied Geophysics，2011，8（1）：1－10.

［25］譚捍東，余欽范，Booker J等.大地電磁法三維交錯采樣有限差分?jǐn)?shù)值模擬.地球物理學(xué)報，2003，46（5）：705－711.Tan H D，Yu Q F，Booker J，et al.Magnetotelluric threedimensional modeling using the staggered－grid finite difference method.Chinese J.Geophys.（in Chinese），2003，46（5）：705－711.