基于MPI瞬變電磁測深一維反演并行算法探究

柳建新，劉鵬茂，劉 穎，童孝忠

(中南大學地球科學與信息物理學院，長沙 410083)

基于MPI瞬變電磁測深一維反演并行算法探究

柳建新，劉鵬茂，劉 穎，童孝忠

(中南大學地球科學與信息物理學院，長沙 410083)

在瞬變電磁測深反演中運用并行技術可以減少計算時間，提高反演的運算效率。MPI(MessagePassingInterface)是目前最重要的并行編程工具，它具有移植性好、功能強大、效率高等多種優點。這里基于在Windows系統下使用FORTRAN和MPICH2相結合的開發工具，編寫瞬變電磁并行算法程序，對瞬變電磁一維采用直接反演法，通過理論模型對該算法進行試算，計算結果證明了該算法的正確性、高效性和穩定性。

瞬變電磁;直接反演;并行計算;MPI

0 前言

二十世紀以來，計算機技術的廣泛應用，極大地促進了科學的發展、社會的進步，推動了計算機的不斷更新和完善。隨著計算機在硬件方面的快速發展，已經大大提高了計算機微處理器的處理速度。在礦產地球物理勘探中，信息技術是勘探的支撐，而對高性能計算的需求和依賴，是勘探信息技術的前進方向。因此，促使了地球物理領域并行計算技術的飛速發展。

并行計算可以加快速度，在更短的時間內解決相同的問題，或者在相同的時間內解決更多更復雜的問題;也可以節省投入，以較低的成本完成同量的任務。因此，并行計算技術將會成為地球物理數據處理的主流，在國內，并行處理各種算法研究取得了一定的成果，陳金窗、戴光明等［17］在微機網絡上實現了2.5維CSAMT正演的并行計算;譚捍東等［18］實現了三維模型的大地電磁正演并行計算，隨后又實現了三維快速松弛反演的并行計算，證明了此方法的有效性和穩定性;劉羽［12］實現了二維大地電磁Occam反演的并行計算;陳露軍、王緒本等［19］以電磁場理論應用為背景，把并行思想和電磁勘探中的電磁計算問題結合起來，研究了在曙光TC－4000微機集群上，基于的分布式并行編程環境下的求解三維電磁場正演模型的問題;王月英［19］實現了三維波動方程的有限元并行正演模擬，采用有限元并行算法，克服了單機對數據容量和計算速度的限制，可以較串行算法模擬更廣泛區域的三維空間;這些都大大地推進了國內并行計算在電磁領域的發展。在國外，Newman和 Alumbaugh［20］在1995年運用并行算法來進行電磁成像;1997年美國洛倫茨貝克利國家實驗室的謝干權等［18］人把并行計算技術運用到全局積分局部微分并行域分解的新算法中進行電磁場的反演，并取得了一定的效果;2000年，Zyserman和 Santos［19］實現了大地電磁三維有限元的并行計算。

A.G.Nekut［7］在 1987 年發表了一篇摘要，主要介紹了瞬變電磁一種簡單的一維直接反演方法。這種方法是把瞬變電磁法測深數據轉換為電導率，測深剖面圖的方法是基于對衰變斡旋電流的近似圖像表示方法，它可以為任何時域電磁測深數據提供一維反演。作者是在Windows系統上使用FORTRAN和MPICH2相結合的開發工具，編寫了瞬變電磁測深的一維直接反演并行程序，通過理論模型對該算法進行試算，計算結果證明了該算法的正確性、可行性和穩定性。可以預知，當計算量比較大時，特別是對于瞬變電磁二維、三維的計算時，這種并行技術優越性可以得到更好的體現。

1 MPI簡介

MPI(MessagePassingInterface)是為了統一不同廠商的消息傳遞API(ApplicationInterface)，有來自高性能計算領域的專家和MPP廠商的代表組成的委員會制定的工業標準，是消息傳遞并行程序設計的標準之一。它能為高性能并行計算提供一個方便的環境，具有移植性好、功能強大、效率高等多種優點。MPI是一個庫，而不是一門語言。按照并行語言的分類，可以把FORTRAN+MPI或C+MPI，看作是一種在原來串行語言基礎之上擴展后得到的并行語言，MPI庫可以被FORTRAN77/C/FORTRAN90、C++調用。

MPI并行程序設計平臺由標準消息傳遞函數及相關輔助函數構成，多個進程通過調用這些函數，進行通信。MPI程序采用的是SPMD執行模式，即一個程序同時啟動多份，形成多個獨立進程，在不同的處理機上運行。每個進程開始執行時。將獲得一個唯一的序號。

MPI的兩種最基本的并行程序設計模式，即對等模式和主從模式。

(1)對等模式即參與運算的各進程地位相同，計算程序一致，只是處理的數據不同。

(2)主從模式則是在程序中設置主進程和從進程，主進程向備從進程發送數據，從進程進行計算。然后將計算結果返回主進程，主進程收到結果后發送結束標志，使從進程退出執行。

2 直接反演方法

2.1 方法原理

應用中心回線測深裝置，在地表上放置半徑為R的一個環形線圈中通有的電流為I，在t=0時刻突然斷電，一個磁力計在環形線圈中心測量中心點垂直方向上的磁場H(t)，它的像源深度D(t)可以用H(t)的變化趨勢來反映［17］，其關系如式(1)。

D=R{［H(t)/H0］－2/3－ 1} (1)式中 H0=I/2R。

依照此方法，通過與放置在距水平良導面正上方δ的電磁源進行類比，可得象源的深度D與探測深度(穿透深度)δ的關系式(2)。

對于電導率為δ的均勻大地介質，瞬變場的有效勘探深度可寫為式(3)。

式(3)等價于頻率域中趨膚深度的公式。

因此，只要給出了地下均勻介質的電導率σ，通過公式(1)、公式(2)、公式(3)，就可以得到H(t)的近似值。

為了校正像源場合瞬變場在均勻半空間條件下的誤差，可以在穿透深度上加一個校正系數F。對于均勻半空間模型，校正系數可以用如下函數來表示:

或

Fˉ(δ/R)與δ/R的關系，我們可以通過均勻半空間中計算 δ/δ得到一組Fˉ(δ/R)和 δ/R，然后通過三次樣條插值法，得到層狀介質情況下的Fˉ。

因此，給定任何一個H(t)，通過公式(1)、公式(2)、公式(3)聯立算出場的有效勘探深度，即

電導率的階數n對應數據的n階連續微分。

一階視電阻率為:

二階視電阻率為:

也可以間接表示為:

圖1 (δ/R)與δ/R的關系圖Fig.1 Thediagramof(δ/R)and δ/R

3 并行算法實現

該瞬變電磁一維并行算法程序中，正演方法采用數值濾波方法，系數為801。反演方法采用的是直接反演法，并把反演結果存儲在 S0、S1、S2.dat文件中。通過對反演解釋中微分算字符的處理，把線性坐標下的微分，轉化為雙曲坐標下的微分，取得的數值更加精確，反演精度S0－S2精度逐漸升高。實現的并行算法程序是在雙核單機下模擬多進程進行，操作系統使用的是Windows系統，PC機配置是 GenuineIntel(R)CPUT1600，主頻1.66GHz，1G 內存。

圖2 (δ/R)與δ/R的關系圖(σ=0.01，10－8＜t＜1s)Fig.2 Thediagramof)F(δ/R)and δ/R(σ =0.01，10－8＜t＜1s)

算例一:對電阻率為三層的理論模型進行并行計算，層數(NC)=3;每層的電導率:σ1=0.01、σ2=0.001、σ3=0.005;每層的深度:H(1)=100、H(2)=100;電流(AMP)=25;天線長度(RAD)=100。并行計算結束需要的時間如表1所示。

表1 理論模型1并行計算不同進程運行時間Tab． 1 Different processes running time of theoreticalmodels 1 parallel computing

該并行計算程序對直接反演方法沒有做任何修改，反演結果如圖3所示。

算例二:層數(NC)=3;每層的電導率:σ1=0.001、σ2=0.01、σ3=0.001;每層的深度:H(1)=100、H(2)=100;電流(AMP)=1;天線長度(RAD)=200。并行計算結束需要的時間如表2所示。

表2 理論模型2并行計算不同進程運行時間Tab． 2 Different processes running time of theoreticalmodels 2 parallel computing

從計算時間結果可以看出，應用并行計算可以加快速度，在更短的時間內解決相同的問題。反演結果如圖4所示。

圖3 模型1電導率—深度反演結果圖Fig． 3 The inversion results map of electrical conductivity-depth of model1

圖4 模型2電導率—深度反演結果圖Fig． 4 The inversion results map of electrical conductivity-depth of model2

4 結論

(1)瞬變電磁一維正演計算，常采用數值濾波方法，即運用漢克爾變換計算頻率域一維響應，然后通過余弦變換法從頻率域轉換到時間域中。數值濾波法具有精度高，時間帶寬長，不同濾波系數之間無相關，因此特別適合并行計算。

(2)通過對均勻半空間模型的數值解與解析解對比，計算其相對誤差，驗證正演算法的正確性。

(3)對瞬變電磁測深一維采用直接反演法，進行并行計算，通過理論模型對實現的并行程序進行試算，計算結果證明了該算法的正確性、高效性和穩定性。

(4)由于反演算法中)F(δ/R)與δ/R的關系是確定的，為了確定)F(δ/R)與δ/R的關系，我們通過均勻半空間中計算)δ/δ得到一組)F(δ/R)與δ/R，采用三次樣條插值法，得到層狀介質情況下任意δ/R值對應的)F(δ/R)。通過對多個均勻模型的計算，我們發現，σ=0.01的均勻半空間模型計算的)F(δ/R)與δ/R關系，可以用來作為其它任何模型的插值基數。

隨著地球物理勘探向深度和廣度方向發展，對勘探精度的要求也越來越高，這就會使需要處理的數據量越來越大，大量的數據量使得并行技術將會在地球物理勘探領域中起到越來越重要的地位和作用。并行計算技術可以加快速度，節省投入，在瞬變電磁反演資料處理技術中，對于其二維、三維的計算，并行技術將會得到更好的應用和發展。

［1］ 羅省賢，何大可.基于MPI的網絡并行計算環境及應用［M］.成都:西南交通大學出版社，2001.

［2］ 李炎，胡祥云.基于MPI的一維大地電磁并行計算研究［J］.地球物理學進展，2010，25(5):1612.

［3］ 謝麗，胡文寶，陸輝.瞬變電磁響應計算的并行算法研究［J］.長江大學學報，2009(6):4.

［4］ 黃易，師學明.并行計算技術及其在勘探地球物理學中的現狀和展望［J］.地球物理學進展，2010，25(2):642.

［5］ 榮瑩，曹俊興.基于MPI的機群并行計算系統平臺構建［J］.物探化探計算技術，2005，27(4):89.

［6］ 牛之璉.時間域電磁法原理［M］.長沙:中南大學出版社，2007.

［7］ NUKET A G． Direct inversion of time． Domain in electromagneticdate． Geophysics， 1987，52( 10) : 1431．

［8］ 瞬變電磁測深法原理［M］.西安:西北工業大學出版社，1993.

［9］ JONES A G． On the equivalence of the“Bostick”and“Niblett”transformations in the magnetotelluric method［J］． Geophys， 1983， 53: 72．

［10］嚴良俊，胡文寶.中心回線瞬變電磁測深法快速電阻率成像方法及應用［J］.煤田地質與勘探，2002，30(6):58.

［11］郭文波，宋建平，韓俊明，等.中心回線瞬變電磁探測的一種快速解釋方法［J］.物探與化探，2006，30(2):616.

［12］ MACNAE J C，LAMONTAGNE Y． Data processingtechniques for EM sounding and prospecting［J］． Geophys，1984， 29( 1) : 464．

［13］都志輝.高性能計算并行編程計算－MPI并行程序設計［M］.北京:清華大學出版社，2001.

［14］孫家旭，張林波，遲學斌，等.網絡并行計算與分布式編程環境［M］.北京:科學出版社，1996.

［15］ MACNAE J C，LAMONTAGNE Y． Imaging quasi． layeredconductive structures by simple processing of transientelectromagnetic data． Geophysics， 1987，52: 45．

［16］石明娟，羅延鐘.西方地面瞬變電磁法理論的發展現狀［J］.國外地質勘探技術，1989(1):255.

［17］陳金窗，戴光明.微機網絡并行計算及2.5維CSAMT正演的并行實現［J］.物探化探計算技術，1997，19(2):103.

［18］ TAN HANDONG，TONG YOU，LIN CHANGHONG． Heparallel 3D magnetotelluric forward modeling algotithm［J］． Applied Geophysics， 2006，3( 4) : 197．

［19］FABIO I，ZYSERMAN JUAN E，SANTOS． Parallel finiteelement algorithm with domain decomposition forthree － dimensional magnetotelluric modeling［J］． J． Appl．Geophys， 2000， 44: 337．

［20］ NEW MAN G A，ALUMBAUGH D L． Three-dimensionalmassiverly parallel electromagnetic inversion［J］．Theory，Geophys，J ． Int， 1997， 128: 345．

P631.3+25

A

1001—1749(2011)05—0491—05

湖南省高校科技創新團隊(2008－244);湖南省重點實驗室資助項目(2010TC2007);湖南省重大專項(2008FJ1006)

2011－01－16 改回日期:2011－06－17

柳建新(1962－)，男，博士，教授，現主要從事大地電磁理論與研究工作。