航空瞬變電磁中MPI并行計算分析

陳 輝， 杜興忠， 楊正剛， 江曉濤

(中國電建集團 貴陽勘測設(shè)計研究院有限公司，貴陽 550081)

航空瞬變電磁中MPI并行計算分析

陳 輝， 杜興忠， 楊正剛， 江曉濤

(中國電建集團 貴陽勘測設(shè)計研究院有限公司，貴陽 550081)

由于航空瞬變電磁的數(shù)據(jù)量很大，外加瞬變電磁反演方法耗時較大地電磁更多，為了在后期資料處理時減少不必要的耗時，對正演計算中的頻時轉(zhuǎn)換的方法做了一定研究，采用基于拉普拉斯尺度變換性質(zhì)的Gaver-Stehfest變換做頻時轉(zhuǎn)換，適當減少了編寫程序代碼的行數(shù)及程序計算的時間消耗。針對現(xiàn)在的計算機很多都是多核的特點，研究了航空瞬變電磁一維反演的MPI(Message Passing Interface)并行計算，大大地節(jié)約了運算時間, 分析了并行計算時不同的多核CPU和不同的進程數(shù)時的加速比及并行效率。

航空瞬變電磁； 頻時轉(zhuǎn)換； G-S變換； 并行； MPI

0 引言

航空瞬變電磁法(Airborne Transient Electromagnetic Method, ATEM)是依托飛機運載，發(fā)射脈沖電流，在發(fā)射機斷開電流時，觀測二次場隨時間的變化特性，進而推測地層中不同物質(zhì)的分布特點，從而解釋各種地質(zhì)問題的方法[1-2]。該方法能夠克服復(fù)雜地形帶來的難以施工的困難，具有低成本、勘探速度快、可大面積勘測、通用性好等優(yōu)點，有利于探測較深目標體和相對低阻覆蓋區(qū)的探測[1,3]。

很多學(xué)者對航空瞬變電磁一維、2.5維和三維正反演方面做了很多的研究，同時2.5維和3維航空瞬變電磁反演也達到較高的計算精度，但其計算效率太低，目前仍沒有得到較好地解決，還在理論研究階段。

在國內(nèi)、外，并行處理技術(shù)的各種算法研究取得了一定的成果。Newman等[4]在電磁反演中成功運用并行算法，完成了多節(jié)點并行運算；陳金窗等[5]在微機網(wǎng)絡(luò)上成功實現(xiàn)了 2.5 維 CSAMT 正演的并行計算；匡斌等[6]實現(xiàn)了三維有限差分深度偏移的并行算法； Siripunvaraporn等[7]完成了對三維大地電磁正演和靈敏度矩陣的并行計算；殷文[8]運用MPI并行技術(shù)對時頻分布做了改進和應(yīng)用研究，使運算效率大大提高；譚捍東等[9]通過頻點并行的方式成功實現(xiàn)了大地電磁三維正演的并行計算，得到了很高的加速比，驗證了并行算法的穩(wěn)定性以及高效性；劉羽[10-12]基于PC機群實現(xiàn)了二維大地電磁Occam反演的并行計算研究；管健和[13]對電磁場有限單元法解釋分布式做了并行計算研究；謝麗等[14]在基于MPI的分布式網(wǎng)絡(luò)并行編程環(huán)境下，成功實現(xiàn)了瞬變電磁響應(yīng)計算的并行算法研究；Siripunvaraporn等[15]研究了大地電磁反演中垂直磁場傳遞函數(shù)反演和并行實現(xiàn)；李焱和胡祥云[16-19]分析了大地電磁的MPI并行計算；張帆[20]利用基于GPU并行計算和基于MPI的并行計算相結(jié)合的計算模式對直流電法和大地電磁法已有的三維正演串行算法進行了并行化處理；張清等[21]實現(xiàn)了基于GPU的疊前時間偏移走時計算的并行算法；李小康[22]研究了頻率域航空電磁法有限單元二維正演的并行計算；柳建新[23]探究了基于MPI瞬變電磁測深一維反演的并行算法；史維等[24]對長時窗長偏移距瞬變電磁測深正反演方法做了并行研究，減少了瞬變電磁法正反演的計算時間；劉鵬茂[25]研究了基于MPI的大地電磁二維正則化反演并行算法；許洋[26]研究了基于MPI的半航空電磁反演的并行算法。

1 ATEM正演計算

對于航空瞬變電磁正演計算，其裝置一般是中心回線裝置，航空瞬變電磁響應(yīng)計算公式為式(1)[27]。

(1)

式中：J1代表一階第一類貝塞爾函數(shù)；a為發(fā)射線圈的半徑長度(m)；λ表示積分的變量；RTE代表反射系數(shù)：

(2)

并且

(3)

(4)

(5)

式中：ω表示角頻率；kn表示第n層的波數(shù)；σn表示第n層的電導(dǎo)率；一般情況下，通常取大地的導(dǎo)磁率μn等于自由空間的導(dǎo)磁率μ0；在準靜態(tài)情況下取μ0=λ。

聯(lián)立以上所有方程式可以完成航空瞬變電磁法一維正演模擬，其計算流程如圖1所示。

圖1 航空瞬變電磁一維正演流程圖Fig.1 ATEM flow chart of one dimensional forward

由圖1可以看出，在上述計算步驟中，涉及到了漢克爾變換和頻時轉(zhuǎn)換。漢克爾變換數(shù)值濾波計算方法比較成熟，經(jīng)過長時間的推敲，擬采用Guptasarma等[28]研究的數(shù)字濾波方法。

查閱文獻[28～30]可知，漢克爾積分變換的一般形式為：

(6)

式中：K(λ)為待變換函數(shù)；Ji為i階第一類貝塞爾函數(shù)；r為一個已給定的值；K(λ)為被積函數(shù)。根據(jù)Guptasarma和Singh的研究得到漢克爾變換的計算公式為式(7):

(7)

其中：λi=(1/r)×10[a+(i-1)s]；i=1、2、…、n;Wi為濾波系數(shù)。

我們利用Guptasarma和Singh研究給出的漢克爾變換濾波系數(shù)，即一階有47或120個系數(shù)，零階有61或120個系數(shù)。

2 頻時轉(zhuǎn)換

2.1 G-S變換

(8)

式中：n是決定于計算機位數(shù)的正偶整數(shù)；Km是G-S變換系數(shù)，由式(9)算出。

Km= (-1)m+n/2·

(9)

式中：m=1、2、…、n，求和下限i1是(m+1)/2的整數(shù)部分。

[30]，我們?nèi)=12,根據(jù)式(9)求得的G-S變換系數(shù)如表1所示。

根據(jù)航空瞬變電磁中心回線裝置，假設(shè)一個5層地電模型HKH，其模型參數(shù)設(shè)定為：電阻率從上到下分別為：100 Ω·m、10 Ω·m、50 Ω·m、5 Ω·m、100 Ω·m；層厚從上到下分別為：200 m、100 m、50 m、500 m、 ∞ m；發(fā)射線圈的半徑為8 m；發(fā)射線圈中供電電流強度為1 A；發(fā)射線圈高度為50 m；

表1 G-S變換系數(shù)

根據(jù)航空瞬變電磁中心回線裝置電磁響應(yīng)計算式(1)和圖1的一維正演計算流程，得到模型HKH的正演瞬變階躍響應(yīng)如圖2所示。

2.2 基于Laplace尺度變換性質(zhì)的G-S變換

在地球物理瞬變電磁法時頻轉(zhuǎn)換中，我們輸入的時間序列通常是按某個公比輸入，傳統(tǒng)的G-S變換要求已知這一串時間數(shù)值，并用讀文件的形式讀入程序參與計算，這樣對于寫代碼和程序計算時會產(chǎn)生額外的時間消耗。這里我們基于拉氏尺度變換性質(zhì)提出一種改進的G-S變換，它只需要輸入初始時刻及時間序列的公比，即可直接在程序中算出所需的N個時刻的值，這樣減少了程序代碼的行數(shù)及讀文件的時間消耗，具有一定的實用價值。其具體做法是：

根據(jù)拉氏變換的尺度變換性質(zhì)[31]：

那么可以將G-S變換式(8)改寫為：

(10)

式(10)即為基于拉氏變換尺度變換性質(zhì)的G-S變換。

圖2 HKH地電模型一維正演電磁響應(yīng)Fig.2 The airborne electromagnetic forward response of HKH model

分別用傳統(tǒng)的G-S變換和基于拉氏尺度變換性質(zhì)的G-S變換做頻時轉(zhuǎn)換計算前面模型HKH，得到的瞬變階躍響應(yīng)如圖3(a)，它們的相對誤差見圖3(b)。

圖3 兩種G-S變換得到的階躍響應(yīng)及它們的相對誤差Fig.3 The step response value of two Gaver-Stehfest transforms and their relative errors(a)階躍響應(yīng)；(b)相對誤差

由圖3中可以看出，基于拉氏尺度變換性質(zhì)的G-S變換得到的結(jié)果和傳統(tǒng)G-S變換的結(jié)果雖然存在一定的差異，但差異非常小，相對誤差最大為0.268%，最小為0%，完全能夠滿足電磁響應(yīng)精度要求，說明這種基于拉氏尺度變換性質(zhì)的G-S變換是正確可行的。

3 ATEM反演并行研究

3.1ATEM反演

在本次航空瞬變電磁一維反演的程序計算中，采用的是阻尼最小二乘法進行反演解釋，它是馬奎特[33](Marquardt)通過改進高斯-牛頓法的法方程組而得到的一種新的反演方法。對于目標函數(shù)為

(11)

的極小問題，阻尼最小二乘法的做法是：

(JTJ+λI)δ=ε

(12)

其中：J表示雅可比矩陣；λ表示阻尼因子(它是一個正數(shù)，用以控制校正的方向和步長)；δ是修正量即δ=pk+1-pk，從式(12)中可以明顯看出它是λ的函數(shù)即δ=δ(λ)，ε是殘差向量。所以阻尼最小二乘法的迭代公式為式(13)。

pk+1=pk+(JTJ+λI)-1ε

(13)

在用阻尼最小二乘法反演解釋時，主要阻尼因子的選取，當阻尼因子λ=0時，阻尼最小二乘法就變?yōu)榱烁咚?牛頓法；當阻尼因子λ=∞時，阻尼最小二乘法就變?yōu)榱俗钏傧陆捣ǎ宰枘嶙钕攵朔ㄐ拚蛄康姆较蚴窃诟咚?牛頓法和最速下降法的修正方向之間，是介于高斯-牛頓法和最速下降法之間的某種插值，而阻尼因子就是這種插值的權(quán)系數(shù)。

Fletcher[34]提出根據(jù)函數(shù)φ(p)在迭代過程中實際的減小量比上，用理想二次型函數(shù)表示φ(p)的變化量，即：

(14)

判別阻尼因子 的增大或減小的標準一般是通過設(shè)定r值的上限和下限來實現(xiàn)的，可以這樣設(shè)定：

1)當rk<0.25，增大阻尼因子λ。

2)當0.25

3)當rk<0.75，減小阻尼因子λ。

當需要增大阻尼因子λk時，可以令λ(k+1)=βλk，當需要減小阻尼因子λk時，可以令λk+1=λk/β。其中的β是一個大于“1”的正數(shù)，一般取為2～10之間的正數(shù)，它可以用式(15)計算出來。

(15)

在地球物理學(xué)反演算法中，雅克比矩陣J的計算精度與反演結(jié)果的準確性有直接關(guān)系，同樣的對于雅克比矩陣的計算速度和計算精度問題，一直都是反演問題中的重點和難點。設(shè)從模型空間映射到數(shù)據(jù)空間的函數(shù)為φ，模型參數(shù)用p表示，則雅克比矩陣J可以表示為：

(16)

式中：m表示觀測數(shù)據(jù)的個數(shù)有m個；n表示模型參數(shù)的個數(shù)有n個。

在航空瞬變電磁法反演程序?qū)崿F(xiàn)中，這里對雅克比矩陣行列中的每個元素的一階偏導(dǎo)數(shù)的計算，采用有限差分方法中的一階向前差分，即有：

(17)

綜上所述，用程序來實現(xiàn)阻尼最小二乘反演的具體流程見圖4。

圖4 阻尼最小二乘法反演解釋流程圖Fig.4 The flow chart of damping least squares inversion interpretation

根據(jù)反演計算流程，分別考慮對三層地電模型H和四層地電模型KH的一維反演解釋。對于 型地電模型，模型從上到下的電阻率分別為100Ω·m、5Ω·m和200Ω·m，各層之間的厚度均為70m，飛行高度為50m。利用H模型正演得到的響應(yīng)數(shù)據(jù)參與阻尼最小二乘反演，得到的反演結(jié)果如圖5所示。對于KH型地電模型，模型從上到下的電阻率分別為10Ω·m、100Ω·m、10Ω·m和100Ω·m，各層之間的厚度均為50m，飛行高度為50m。利用 模型正演得到的響應(yīng)數(shù)據(jù)參與阻尼最小二乘反演，得到的反演結(jié)果如圖6所示。

圖5 H型地電模型反演結(jié)果Fig.5 The inversion results of model H

圖6 KH型地電模型反演結(jié)果Fig.6 The inversion results of model KH

由圖5、圖6可以看出，整個反演的效果是非常明顯的，但是在對于地層中存在高阻情況時，反演得到的結(jié)果并不是很理想，對高阻體的反演解釋不明顯。造成這樣的情況應(yīng)該是多方面的綜合影響結(jié)果，但最主要的可能有以下兩個方面的原因：①在航空瞬變電磁法勘探中，對于地下的低阻情況是比較敏感的，而對于地下存在的高阻情況卻不敏感，當?shù)叵麓嬖诟咦梵w時，由于航空瞬變電磁對其不敏感的特點，理論上的感應(yīng)電動勢響應(yīng)就不明顯，而其經(jīng)過反演后所得的反演結(jié)果模型就會和理論模型之間存在較大的差異；②在反演解釋方法的選擇上，本次實驗選擇的反演解釋方法是阻尼最小二乘反演法，其實也可以選擇其他的反演方法，因為各種反演解釋方法有不同的特點，所以應(yīng)該根據(jù)實際情況來選擇合適的反演方法，當然有時候為了取得更好的反演效果，還可以將多種反演方法配合使用。

3.2MPI簡述及其性能評價

MPI(MessagePassingInterface)是消息傳遞方式并行程序設(shè)計的標準之一，它是一個函數(shù)庫規(guī)范，并不是一門并行語言，它的操作就像是庫函數(shù)的調(diào)用。我們完成用FORTRAN+MPI的并行語言編程，實現(xiàn)了航空瞬變電磁一維反演的并行計算研究。

對一個并行程序來說，效率評估的衡量參數(shù)有加速比與效率這倆最常用到的參數(shù)。先假設(shè)并行程序的每個進程獨享處理器資源，在這樣的前提下，設(shè)某個串行程序在某臺并行機單處理器上的執(zhí)行的時間為ts，對該程序并行化后，p個進程在p個處理器并行執(zhí)行所消耗的時間為tp，那么該并行程序在該并行機上的加速比sp可定義為[32]：

(18)

并行效率定義為式(19)。

(19)

3.3 MPI并行計算分析

由于反演中每次迭代都要計算雅克比矩陣J，都要重新修正初始模型，還要做正演計算，所以反演部分的耗時是最多的。因此，在設(shè)計其MPI并行程序主要還是對反演部分的雅可比矩陣計算做并行計算。采用主從模式，其并行計算的具體實現(xiàn)過程見圖7。

圖7 ATEM并行計算流程圖Fig.7 The flow chart of ATEM parallel computing

本次并行計算設(shè)計的算例為前面的四層地電模型 模型。分別考慮了在雙核單機下模擬多進程進行MPI并行的情況和不同處理器或者說不同線程相同進程進行MPI并行的情況，在這兩種情況下分別比較了它們的加速比和并行效率，并進行了一定的分析。

對于雙核單機下模擬多進程的情況，個人計算機的操作系統(tǒng)是Windows系統(tǒng)，配置是Pentium(R) Dual-Core CPU T4300 @ 2.10GHz處理器，內(nèi)存是3.00GB。反演的并行計算結(jié)束所需要的時間、加速比和并行效率表2所示。

并行計算得到的反演結(jié)果并沒有什么改變，如圖8所示。

表2 KH模型反演并行計算不同進程運行結(jié)果

圖8 并行計算得到反演結(jié)果模型Fig.8 Inverse model of parallel computing

從表2和圖8可以看出，基于MPI的并行計算加快了運行的速度，使得耗費的時間變短，并沒有改變反演的計算結(jié)果，說明這種并行是正確可行的。從表2中還可以看出，對同一種計算機或者說同一種處理器而言，此時增加并行的進程數(shù)，所耗費的時間就越短，加速比就越大，但是并行效率卻在減小。這是因為對于p個進程而言，在理想的情況下，可以得到的最大加速比為p即線性加速，得到的效率為100%，但是此時并沒有考慮到并行計算中的各種開銷。事實上，隨著進程數(shù)量的增加，由各種開銷引起的耗時也會越多，計算通信比將大幅度減小，因而造成效率隨著進程數(shù)增大而減小的現(xiàn)象。

對于不同的核心CPU而言，這里討論三種核心CPU情況，分別為Pentium(R) Dual-Core CPU T4300 @ 2.10GHz、Intel(R) CoreTMi5 CPU M480 @ 2.67GHz和Intel(R) CoreTMi7-2860QM CPU @ 2.50GHz，它們分別是雙核2線程、雙核4線程和4核8線程。分別用它們做并行計算，這里只考慮1進程和8進程的反演部分并行計算情況，得到的反演計算結(jié)果如表3所示。

表3 KH模型反演并行計算不同線程運行結(jié)果

從表3中可以看出，對于不同核心CPU的時候，處理器的線程數(shù)越大所消耗的時間相對越少，加速比相對越大，并行效率也相對越高,同樣它們經(jīng)過并行計算得到的反演結(jié)果并沒有被改變。

4 結(jié)論

1)通過對正演計算中的頻時轉(zhuǎn)換研究，得到基于拉普拉斯尺度變換性質(zhì)的G-S變換是正確可行的，它在保證滿足航空瞬變電磁響應(yīng)精度的前提下，減少了編寫程序代碼的行數(shù)及程序的計算時間消耗，具有一定的實用價值。

2)阻尼最小二乘法反演得到的反演解釋結(jié)果可以看出，對于三層和四層地電模型而言，總體反演效果比較明顯，但是在對于地層中存在高阻情況時，反演得到的結(jié)果并不是很理想，對高阻體的反演解釋不明顯。

3)對航空瞬變電磁法的反演部分做基于MPI的并行計算，和傳統(tǒng)反演計算相比，并行計算加快了程序運算的速度，減少了計算的耗費時間，并不會改變正演和反演計算的結(jié)果。不同的處理器并行計算的耗費時間不一樣，處理器的核心線程數(shù)越大，相對的耗時會越少，加速比相對越大，并行效率也相對越高。相同處理器不同進程數(shù)時，進程數(shù)量越多，相對耗費的總體時間越少，加速比相對越大，但是并行效率卻在相對變低。

參考文獻：

[1] NABIGHIAN M N. Electromagnetic methods in applied geophysics: volume 2, application, parts A and B [M]. Society of Exploration Geophysics, 1991.

[2] 王衛(wèi)平. 頻率域航空電磁法特點及應(yīng)用[C].中國地球物理學(xué)會第二十六屆年會、中國地震學(xué)會第十三次學(xué)術(shù)大會, 2010: 635-636. WANG W P. Frequency domain characteristics and its application of airborne electromagnetic method[C]. The Chinese geophysical society the 26th annual meeting and China earthquake to the 13th academic conference, 2010: 635-636.(In Chinese)

[3] MACNAE J C, NABIGHIAN M. Electrical and EM methods, 1980-2005 [J]. The Leading Edge, 2005, 24(s1): 42-45.

[4] NEWMAN G A, ALUMBAUGH D L. Three-dimensional massively parallel electromagnetic inversion[J].Geophysics, 1997, 128: 345-354.

[5] 陳金窗,戴光明.微機網(wǎng)絡(luò)并行計算及2.5維CSAMT正演的并行實現(xiàn)[J].物探化探計算技術(shù),1997,19(2):103-107. CHEN J C,DAI G M.Micro-computer networked computing and 2.5-d CSAMT forward parallel computing[J].Computing Techniques for Geophysical and Geochemical Exploration, 1997,19(2):103-107.(In Chinese)

[6] 匡斌, 李心友, 王真理, 等. 三維有限差分深度偏移的并行算法和實現(xiàn)[C]. 寸丹集-慶賀劉光鼎院士工作50周年學(xué)術(shù)論文集, 1998,12: 289-303. KUANG B, LI X Y, WANG Z L, et al. The design and implementation of parallel algorithm of 3-D Finite-Difference depth migration[C]. Cun Dan set - celebrate Liu Guangding academician work the 50th anniversary of the academic papers, 1998, 12: 289-303. (In Chinese)

[7] SIRIPUNVARAPORN W, UYESHIMA M,EGBERT G.Three-demensional inversion for Network-Magnetotelluric data[J].Earth Planets Space,2004,56:893-902.

[8] 殷文, 印興耀. 基于MPI 的時頻分布的改進及應(yīng)用[J]. 地球物理學(xué)進展, 2005, 20(1): 165-169. YIN W, YIN X Y. The amelioration and application of time frequency distributions on the basis of MPI[J]. Progress in Geophysics, 2005, 20(1):165-169. (In Chinese)

[9] TAN H D,TONG T,LIN C H.The parallel 3D magnetotelluric forward modeling algorithm[J]. Applied Geophysics,2006,3(4):197-202.

[10]劉羽. 大地電磁Occam反演中拉格朗日乘子搜索的并行計算[J]. 物探化探計算技術(shù), 2005, 28(4): 342-345. LIU Y. The parallel computation of l l SEARCHING in magnetotelluric occam inversion[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2005, 28(4): 342-345. (In Chinese)

[11]劉羽, 王家映, 孟永良. 基于PC機群的大地電磁Occam反演并行計算研究[J]. 石油物探, 2006, 45(3): 311-315. LIU Y, WANG J Y, MENG Y L. PC cluster based magnetotelluric 2-D Occam's inversion parallel calculation[J]. GPP, 2006, 45(3): 311-315. (In Chinese)

[12]劉羽. MT Occam并行反演方案及性能分析[J]. 武漢理工大學(xué)學(xué)報, 2007, 29(2): 136-140. LIU Y. Parallel MT Occam inversion scheme and its performance analysis[J]. Journal of Wuhan university of technology, 2007, 29(2): 136-140. (In Chinese)

[13]管建和. 電磁場有限元法解釋分布式并行計算的研究[D]. 北京: 中國地質(zhì)大學(xué), 2006. GUAN J H. Study of distributing and parallel computation of electromagnetic field finite method explaining[D]. Beijing:China University of Geosciences, 2006. (In Chinese)

[14]謝麗, 胡文寶, 陸輝. 瞬變電磁響應(yīng)計算的并行算法研究[J]. 長江大學(xué)學(xué)報(自然科學(xué)報), 2009, 6(4): 22-24. XIE L, HU W B, LU H. A study of parallel algorithm on the calculation of transient electromagnetic response[J]. Journal of Yangtze University (Nat Sci Edit) , 2009, 6(4): 22-24. (In Chinese)

[15]SIRIPUNVARAPORN W, EGBERT G. WSINV 3D MT:Vertical magnetic field transfer function inversion and parallel implementation[J]. Physics of the Earth and Planetary Interiors, 2009, 173: 317-329.

[16]李焱, 胡祥云, 楊文采,等. 大地電磁三維交錯網(wǎng)格有限差分數(shù)值模擬的并行計算研究[J]. 地球物理學(xué)報, 2012, 55(12): 4036-4043. LI Y, HU X Y, YANG W C, et al. A study on parallel computation for 3D magnetotelluric modeling using the staggered-grid finite difference method[J]. Chinese Journal of Geophysics, 2012, 55(12): 4036-4043. (In Chinese)

[17]李焱, 胡祥云, 金剛燮,等. 基于MPI的一維大地電磁并行計算研究[J]. 地球物理學(xué)進展, 2010, 25(5): 1612-1616. LI Y, HU X Y, KIMKANGSOP. Research of 1-D magnetotelluric parallel computation based on MPI[J]. Progress in Geophysics, 2010, 25(5): 1612-1616. (In Chinese)

[18]李焱, 胡祥云, 吳桂桔,等. 基于MPI的二維大地電磁正演的并行計算[J]. 地震地質(zhì), 2010, 32(3): 392-401. LI Y, HU X Y, WU G J,et al. Parallel computation of 2-d magnetotelluric forward modeling based on mpi[J]. Seismology and Geology, 2010, 32(3): 392-401.(In Chinese)

[19]胡祥云, 李焱, 楊文采,等. 大地電磁三維數(shù)據(jù)空間反演并行算法研究[J]. 地球物理學(xué)報, 2012, 55(12): 3969-3978. HU X Y, LI Y, YANG W C,et al. Three-dimensional magnetotelluic parallel inversion algorithm using data space method[J]. Chinese Journal of Geophysics, 2012, 55(12): 3969-3978. (In Chinese)

[20]張帆. 基于MPI和GPU直流電法和大地電磁法三維正演的并行算法研究[D]. 北京: 中國地質(zhì)大學(xué), 2011. ZHANG F. Research on parallel algorithm for three dimentional forward modeling of DC resistivity method and magnetotelluric method based on MPI and GPU[D]. Beijing:China University of Geoscience , 2011. (In Chinese)

[21]張清, 遲旭光, 謝海波,等. 基于GPU實現(xiàn)疊前時間偏移走時計算的并行算法[J]. 計算機系統(tǒng)應(yīng)用, 2011, 20(8): 42-46. ZHANG Q, CHI X G, XIE H B.,et al.Parallel algorithm based on the travel time computing of pre-stack time migration using GPU[J]. Computer Systems & Applications, 2011, 20(8):42-46. (In Chinese)

[22]李小康. 基于MPI的頻率域航空電磁法有限元二維正演并行計算研究[D]. 北京: 中國地質(zhì)大學(xué), 2011. LI X B. A MPI based parallel calculation investigation on two dimensional finite element modelling of AEM[D]. Beijing: China University of Geosciences, 2011. (In Chinese)

[23]柳建新, 劉鵬茂, 劉穎,等. 基于MPI瞬變電磁測深一維反演并行算法探究[J]. 物探化探計算技術(shù), 2011, 33(5): 491-495. LIU J X, LIU P M, LIU Y,et al. Research of parallel computing technology about the direct one-dimensional inversion of transient electromagnetic on MPI[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2011, 33(5): 491-495. (In Chinese)

[24]史維, 嚴良俊, 謝興兵. MPI環(huán)境的并行算法在瞬變電磁正反演中的研究[J]. 工程地球物理學(xué)報, 2012, 9(1): 75-80. SHI W, YAN L J, XIE X B. Parallel algorithm of MPI in TEM forward and inversion[J]. Chinese Journal of Engineering Geophysics, 2012, 9(1): 75-80. (In Chinese)

[25]劉鵬茂. 基于MPI大地電磁二維正則化反演并行算法研究[D]. 長沙: 中南大學(xué), 2012. LIU P M. Parallel regularization inversion of 2-D MT data based on MPI[D]. Changsha: Central South University, 2012. (In Chinese)

[26]許洋. 半航空電磁一維正反演研究[D]. 成都: 成都理工大學(xué), 2014. XU Y. Study about 1D forward and inversion of SAEM[D]. Chengdu: Chengdu University of Technology, 2014. (In Chinese)

[27]NABIGHIAN M N. Electromagnetic methods in applied geophysics: volume 1, theory [M]. Society of Exploration Geophysics，1988.

[28]GUPTASARMA D, SINGH B. New digital linear filter for Hankel J0 and J1 transforms [J]. GEOPHYSICAL PROSPECTING, 1997, 45(5): 745-762.

[29]KNIGHT J H and RAICHE A P. Transient electromagnetic calculations using the Gaver-Stehfest inverse Laplace transform method [J]. GEOPHYSICS, 1982, 47(1): 47-50.

[30]羅延鐘, 昌彥君. G-S變換的快速算法[J]. 地球物理學(xué)報, 2000, 43(5): 684-690. LUO Y ZH, CHANG Y J. A rapid algorithm for G-S transform[J]. Chinese Journal of Geophysics, 2000, 43 (5): 684-690. (In Chinese)

[31]鄭君里, 應(yīng)啟珩, 楊為理. 信號與系統(tǒng)(上冊)[M]. 北京：高等教育出版社, 2000. ZHENG L J, YING Q H, YANG W L. Signals and systems[M]. Beijing:CHEP, 2000. (In Chinese)

[32]張林波, 遲學(xué)斌, 莫則堯,等. 并行計算導(dǎo)論[M]. 北京: 清華大學(xué)出版社, 2006. ZHANG L B, CHI X B, MO Z Y. Introduction to parallel computing[M]. Beijing:Tsinghua University Press, 2006. (In Chinese)

[33]陳樂壽,劉任,王天生. 大地電磁測深資料處理與解釋[M]. 北京: 石油工業(yè)出版社, 1988. Chen L S, Liu R, Wang T S. The processing and interpretation of magnetotelluric sounding data [M]. Beijing: Petroleum Industry Press, 1988.(In Chinese)

[34]FLETCHER R.A Modified Marquardt Subroutine for Nonlinear Least Squares [M].Harwell:Rpt.AERE,1971.

第39卷 第2期2017年3月物探化探計算技術(shù)COMPUTING TECHNIQUES FOR GEOPHYSICAL AND GEOCHEMICAL EXPLORATIONVol.39 No.2Mar. 2017

MPI parallel computing in airborne transient electromagnetic method

CHEN Hui, DU Xingzhong, YANG Zhenggang, JIANG Xiaotao

(Guiyang Hydropower investigation design and Research Institute, CHECC, Guiyang 550081, China)

Becasue airborne transient electromagnetic data volume is huge and transient electromagnetic inversion method is more compared with the magnetotelluric time-consuming. It is necessary to apply parallel computing of the research on airborne transient electromagnetic to reduce the unnecessary time in data processing. With the advantages of most computers with multi-cores, parallel computing research can be easily used in airborne transient electromagnetic. We first analyzed the frequency-time conversion method in the forward calculation,using Gaver-Stehfest transform based on Laplace scale transform nature do frequency-time conversion. This appropriates to reduce the number of rows to write program code and program calculation of time consumption. And then studied airborne transient electromagnetic one dimensional inversion of MPI (Message Passing Interface) with parallel omputing, which greatly saves the operation time. The parallel computing is then analyzed with multicore cpus, different processes considering the speedup ratio and parallel efficiency.

airbone transient electromagnetic; frequency-time conversion; gaver-stehfest transform; parallel computing; message passing interface

2016-05-09 改回日期：2016-07-25

貴州省科學(xué)技術(shù)基金([2013]2297)

陳輝(1988-)，男，碩士，工程師，主要從事工程物探技術(shù)研究和應(yīng)用， E-mail：281636258@qq.com。

1001-1749(2017)02-0161-09

P 631.2

A

10.3969/j.issn.1001-1749.2017.02.02