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

(中國(guó)電建集團(tuán) 貴陽(yáng)勘測(cè)設(shè)計(jì)研究院有限公司，貴陽(yáng) 550081)

航空瞬變電磁中MPI并行計(jì)算分析

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

(中國(guó)電建集團(tuán) 貴陽(yáng)勘測(cè)設(shè)計(jì)研究院有限公司，貴陽(yáng) 550081)

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

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

0 引言

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

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

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

1 ATEM正演計(jì)算

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

(1)

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

(2)

并且

(3)

(4)

(5)

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

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

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

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

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

(6)

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

(7)

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

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

2 頻時(shí)轉(zhuǎn)換

2.1 G-S變換

(8)

式中：n是決定于計(jì)算機(jī)位數(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è)一個(gè)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ā)射線圈中供電電流強(qiáng)度為1 A；發(fā)射線圈高度為50 m；

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

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

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

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

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

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

(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變換做頻時(shí)轉(zhuǎn)換計(jì)算前面模型HKH，得到的瞬變階躍響應(yīng)如圖3(a)，它們的相對(duì)誤差見(jiàn)圖3(b)。

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

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

3 ATEM反演并行研究

3.1ATEM反演

在本次航空瞬變電磁一維反演的程序計(jì)算中，采用的是阻尼最小二乘法進(jìn)行反演解釋?zhuān)邱R奎特[33](Marquardt)通過(guò)改進(jìn)高斯-牛頓法的法方程組而得到的一種新的反演方法。對(duì)于目標(biāo)函數(shù)為

(11)

的極小問(wèn)題，阻尼最小二乘法的做法是：

(JTJ+λI)δ=ε

(12)

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

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

(13)

在用阻尼最小二乘法反演解釋時(shí)，主要阻尼因子的選取，當(dāng)阻尼因子λ=0時(shí)，阻尼最小二乘法就變?yōu)榱烁咚?牛頓法；當(dāng)阻尼因子λ=∞時(shí)，阻尼最小二乘法就變?yōu)榱俗钏傧陆捣?，所以阻尼最想二乘法修正向量的方向是在高?牛頓法和最速下降法的修正方向之間，是介于高斯-牛頓法和最速下降法之間的某種插值，而阻尼因子就是這種插值的權(quán)系數(shù)。

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

(14)

判別阻尼因子 的增大或減小的標(biāo)準(zhǔn)一般是通過(guò)設(shè)定r值的上限和下限來(lái)實(shí)現(xiàn)的，可以這樣設(shè)定：

1)當(dāng)rk<0.25，增大阻尼因子λ。

2)當(dāng)0.25

3)當(dāng)rk<0.75，減小阻尼因子λ。

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

(15)

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

(16)

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

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

(17)

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

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

根據(jù)反演計(jì)算流程，分別考慮對(duì)三層地電模型H和四層地電模型KH的一維反演解釋。對(duì)于 型地電模型，模型從上到下的電阻率分別為100Ω·m、5Ω·m和200Ω·m，各層之間的厚度均為70m，飛行高度為50m。利用H模型正演得到的響應(yīng)數(shù)據(jù)參與阻尼最小二乘反演，得到的反演結(jié)果如圖5所示。對(duì)于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可以看出，整個(gè)反演的效果是非常明顯的，但是在對(duì)于地層中存在高阻情況時(shí)，反演得到的結(jié)果并不是很理想，對(duì)高阻體的反演解釋不明顯。造成這樣的情況應(yīng)該是多方面的綜合影響結(jié)果，但最主要的可能有以下兩個(gè)方面的原因：①在航空瞬變電磁法勘探中，對(duì)于地下的低阻情況是比較敏感的，而對(duì)于地下存在的高阻情況卻不敏感，當(dāng)?shù)叵麓嬖诟咦梵w時(shí)，由于航空瞬變電磁對(duì)其不敏感的特點(diǎn)，理論上的感應(yīng)電動(dòng)勢(shì)響應(yīng)就不明顯，而其經(jīng)過(guò)反演后所得的反演結(jié)果模型就會(huì)和理論模型之間存在較大的差異；②在反演解釋方法的選擇上，本次實(shí)驗(yàn)選擇的反演解釋方法是阻尼最小二乘反演法，其實(shí)也可以選擇其他的反演方法，因?yàn)楦鞣N反演解釋方法有不同的特點(diǎn)，所以應(yīng)該根據(jù)實(shí)際情況來(lái)選擇合適的反演方法，當(dāng)然有時(shí)候?yàn)榱巳〉酶玫姆囱菪Ч?，還可以將多種反演方法配合使用。

3.2MPI簡(jiǎn)述及其性能評(píng)價(jià)

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

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

(18)

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

(19)

3.3 MPI并行計(jì)算分析

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

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

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

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

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

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

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

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

對(duì)于不同的核心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線程。分別用它們做并行計(jì)算，這里只考慮1進(jìn)程和8進(jìn)程的反演部分并行計(jì)算情況，得到的反演計(jì)算結(jié)果如表3所示。

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

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

4 結(jié)論

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

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

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

參考文獻(xiàn)：

[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)平. 頻率域航空電磁法特點(diǎn)及應(yīng)用[C].中國(guó)地球物理學(xué)會(huì)第二十六屆年會(huì)、中國(guó)地震學(xué)會(huì)第十三次學(xué)術(shù)大會(huì), 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] 陳金窗,戴光明.微機(jī)網(wǎng)絡(luò)并行計(jì)算及2.5維CSAMT正演的并行實(shí)現(xiàn)[J].物探化探計(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] 匡斌, 李心友, 王真理, 等. 三維有限差分深度偏移的并行算法和實(shí)現(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 的時(shí)頻分布的改進(jìn)及應(yīng)用[J]. 地球物理學(xué)進(jìn)展, 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ì)算[J]. 物探化探計(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機(jī)群的大地電磁Occam反演并行計(jì)算研究[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é)報(bào), 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]管建和. 電磁場(chǎng)有限元法解釋分布式并行計(jì)算的研究[D]. 北京: 中國(guó)地質(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ì)算的并行算法研究[J]. 長(zhǎng)江大學(xué)學(xué)報(bào)(自然科學(xué)報(bào)), 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]李焱, 胡祥云, 楊文采,等. 大地電磁三維交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬的并行計(jì)算研究[J]. 地球物理學(xué)報(bào), 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ì)算研究[J]. 地球物理學(xué)進(jìn)展, 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ì)算[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é)報(bào), 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]. 北京: 中國(guó)地質(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實(shí)現(xiàn)疊前時(shí)間偏移走時(shí)計(jì)算的并行算法[J]. 計(jì)算機(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的頻率域航空電磁法有限元二維正演并行計(jì)算研究[D]. 北京: 中國(guó)地質(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瞬變電磁測(cè)深一維反演并行算法探究[J]. 物探化探計(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]史維, 嚴(yán)良俊, 謝興兵. MPI環(huán)境的并行算法在瞬變電磁正反演中的研究[J]. 工程地球物理學(xué)報(bào), 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]. 長(zhǎng)沙: 中南大學(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é)報(bào), 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)啟珩, 楊為理. 信號(hào)與系統(tǒng)(上冊(cè))[M]. 北京：高等教育出版社, 2000. ZHENG L J, YING Q H, YANG W L. Signals and systems[M]. Beijing:CHEP, 2000. (In Chinese)

[32]張林波, 遲學(xué)斌, 莫?jiǎng)t堯,等. 并行計(jì)算導(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]陳樂(lè)壽,劉任,王天生. 大地電磁測(cè)深資料處理與解釋[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月物探化探計(jì)算技術(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