大地電磁三維交錯網格有限差分數值模擬的并行計算研究

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

1 中國地質大學（武漢）地球物理與空間信息學院，武漢 430074

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

3 中國地質科學院，北京 100037

4 中國地質大學（北京）地球物理與信息技術學院，北京 100083

1 引 言

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

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

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

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

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

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

3 大地電磁三維交錯網格有限差分數值模擬

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

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

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

圖1 交錯采樣網格示意圖Fig.1 Sketch of Staggered-grid

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

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

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

4 三維正演數值模擬并行算法的總體設計

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

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

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

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

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

如圖3所示：該模型為3D／2D模型，即區域構造是2D，含有3D異常體.棱柱體大小為2km×2km×1km，電阻率為5Ωm，頂面埋深為1km，中心位于（0，－3000m，1500m），覆存于電阻率為100Ωm的二維構造中，二維走向為y方向，基底電阻率為1000Ωm.三維交錯網格有限差分正演在x、y、z方向剖分網格單元數分別為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軸方向的響應擬斷面圖.在XY模式和YX模式視電阻率擬斷面圖上，3D／2D模型的YX模式較好的反映出三維棱柱體的形態和范圍，特別是埋深和底界面，而XY模式反映的三維棱柱體明顯向下拉伸.XY模式和YX模式相位擬斷面圖都較好的反映出三維棱柱體的形態和范圍，其中XY模式橫向分辨率比YX模式橫向分辨率高，兩種模式對2D構造及基底形態反映不明顯.從XX模式和YY模式視電阻率擬斷面圖上可以看出兩種模式的視電阻率值都非常小，不能反映模型的形態.因為XX模式和YY模式視電阻率值都非常小，再加上正演模擬過程中解方程時所給定的收斂判別標準誤差，因此XX模式和YY模式的相位信息是不可靠的，也無法反映模型的形態.從這里還可以看出阻抗信息非對角元素的作用遠大于對角元素的作用.

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

設計的第二個模型如圖6所示：高阻三維棱柱體電阻率為1000Ωm，低阻三維棱柱體電阻率為10Ωm，它們大小相同，均為2km×2km×1km，頂面埋深為0.5km，中心分別為（3km，0，1km）和（－3km，0，1km）圍巖電阻率為100Ωm.三維交錯網格有限差分正演在x、y、z方向剖分網格單元數分別為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是高阻體和低阻體組合模型的視電阻率和相位兩種模式響應擬斷面圖.由于XX模式和YY模式的信息不能較好的反映出異常體，這里只給出XY模式和YX模式的響應擬斷面圖.在XY模式和YX模式視電阻率擬斷面圖上，較好的反映低阻體和高阻體，相比較而言右邊的高阻體響應不如左邊的低阻體明顯.兩種模式的視電阻率擬斷面圖都較好的反映了異常體的頂深和橫向范圍，其中XY模式比YX模式在橫向范圍上反映更準確一些，YX模式在橫向上有一定延伸.兩種模式的相位響應擬斷面圖均能較好的反映出低阻體和高阻體的整個范圍，尤其是能反映視電阻率無法反映的底界面.

5.3 并行程序正確性的驗證

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

5.4 并行效率分析

為了檢驗開發的三維正演并行程序的效率，采取不同的節點數來計算模型一和模型二，并通過并行加速比和并行效率來評價并行算法的效果.并行加速比＝單機串行算法的執行時間／N個進程并行算法的執行時間；并行效率＝并行加速比／參加并行計算的進程個數.表1是3D／2D模型正演測試時間統計表；表2是低阻體高阻體組合模型正演測試時間統計表.

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

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

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

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

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

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

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

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

6 結 語

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

（References）

［1］胡祥云，胡祖志，張榮等.油氣MT勘探COPROD－2S1模型數據的二維反演.天然氣工業，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等.西藏高原南部雅魯臧布江縫合帶地區地殼電性結構研究.地球物理學報，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］葉益信，胡祥云，金鋼燮等.大地電磁二維陡邊界反演應用效果分析.地球物理學進展，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］胡祖志，胡祥云，何展翔.大地電磁非線性共軛梯度擬三維反演.地球物理學報，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］楊迪琨，胡祥云.含噪聲數據反演的概率描述.地球物理學報，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等.大地電磁法三維快速松弛反演.地球物理學報，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］張林波，遲學斌，莫則堯等.并行計算導論.北京：清華大學出版社，2006.Zhang L B，Chi X B，Mo Z Y，et al.Introduction to Parallel Computing.Beijing：Tsinghua University Press，2006.

［8］都志輝 編著.高性能計算并行編程技術－MPI并行程序設計.北京：清華大學出版社，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］胡祖志，胡祥云.大地電磁三維反演方法綜述.地球物理學進展，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等.大地電磁法三維交錯采樣有限差分數值模擬.地球物理學報，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.