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

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

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

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

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

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

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

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

0 前言

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

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

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

1 MPI簡介

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

MPI并行程序設(shè)計平臺由標(biāo)準(zhǔn)消息傳遞函數(shù)及相關(guān)輔助函數(shù)構(gòu)成，多個進(jìn)程通過調(diào)用這些函數(shù)，進(jìn)行通信。MPI程序采用的是SPMD執(zhí)行模式，即一個程序同時啟動多份，形成多個獨立進(jìn)程，在不同的處理機(jī)上運(yùn)行。每個進(jìn)程開始執(zhí)行時。將獲得一個唯一的序號。

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

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

(2)主從模式則是在程序中設(shè)置主進(jìn)程和從進(jìn)程，主進(jìn)程向備從進(jìn)程發(fā)送數(shù)據(jù)，從進(jìn)程進(jìn)行計算。然后將計算結(jié)果返回主進(jìn)程，主進(jìn)程收到結(jié)果后發(fā)送結(jié)束標(biāo)志，使從進(jìn)程退出執(zhí)行。

2 直接反演方法

2.1 方法原理

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

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

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

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

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

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

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

或

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

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

電導(dǎo)率的階數(shù)n對應(yīng)數(shù)據(jù)的n階連續(xù)微分。

一階視電阻率為:

二階視電阻率為:

也可以間接表示為:

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

3 并行算法實現(xiàn)

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

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

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

表1 理論模型1并行計算不同進(jìn)程運(yùn)行時間Tab． 1 Different processes running time of theoreticalmodels 1 parallel computing

該并行計算程序?qū)χ苯臃囱莘椒]有做任何修改，反演結(jié)果如圖3所示。

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

表2 理論模型2并行計算不同進(jìn)程運(yùn)行時間Tab． 2 Different processes running time of theoreticalmodels 2 parallel computing

從計算時間結(jié)果可以看出，應(yīng)用并行計算可以加快速度，在更短的時間內(nèi)解決相同的問題。反演結(jié)果如圖4所示。

圖3 模型1電導(dǎo)率—深度反演結(jié)果圖Fig． 3 The inversion results map of electrical conductivity-depth of model1

圖4 模型2電導(dǎo)率—深度反演結(jié)果圖Fig． 4 The inversion results map of electrical conductivity-depth of model2

4 結(jié)論

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

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

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

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

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

［1］ 羅省賢，何大可.基于MPI的網(wǎng)絡(luò)并行計算環(huán)境及應(yīng)用［M］.成都:西南交通大學(xué)出版社，2001.

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

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

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

［5］ 榮瑩，曹俊興.基于MPI的機(jī)群并行計算系統(tǒng)平臺構(gòu)建［J］.物探化探計算技術(shù)，2005，27(4):89.

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

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

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

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

［10］嚴(yán)良俊，胡文寶.中心回線瞬變電磁測深法快速電阻率成像方法及應(yīng)用［J］.煤田地質(zhì)與勘探，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并行程序設(shè)計［M］.北京:清華大學(xué)出版社，2001.

［14］孫家旭，張林波，遲學(xué)斌，等.網(wǎng)絡(luò)并行計算與分布式編程環(huán)境［M］.北京:科學(xué)出版社，1996.

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

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

［17］陳金窗，戴光明.微機(jī)網(wǎng)絡(luò)并行計算及2.5維CSAMT正演的并行實現(xiàn)［J］.物探化探計算技術(shù)，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

湖南省高?？萍紕?chuàng)新團(tuán)隊(2008－244);湖南省重點實驗室資助項目(2010TC2007);湖南省重大專項(2008FJ1006)

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

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