沈 磊， 董春暉,2

(1.山東交通職業(yè)學(xué)院 公路與建筑學(xué)院，濰坊261206； 2.大連海事大學(xué) 道路與橋梁工程研究所，大連 116026)

基于GPU的航空瞬變電磁法響應(yīng)模擬

沈 磊1， 董春暉1,2

(1.山東交通職業(yè)學(xué)院 公路與建筑學(xué)院，濰坊261206； 2.大連海事大學(xué) 道路與橋梁工程研究所，大連 116026)

針對(duì)航空瞬變電磁法正演算法結(jié)構(gòu)與GPU的并行性運(yùn)算相結(jié)合，推導(dǎo)出了航空瞬變電磁法一維正演公式，在2007年NVIDIA公司推出的CUDA平臺(tái)上實(shí)現(xiàn)了一維正演算法在GPU上的調(diào)用。為驗(yàn)證CUDA C語(yǔ)言程序的正確性，設(shè)計(jì)了幾個(gè)模型，并與C代碼模擬的響應(yīng)做了對(duì)比，結(jié)果表明，基于GPU的CUDA C代碼所編的程序在航空瞬變電磁法的模擬中擁有明顯優(yōu)勢(shì)，速度提高非常明顯，在復(fù)雜的地球物理學(xué)模型計(jì)算和反演中，擁有巨大的發(fā)展?jié)摿Α?/p>

GPU； CUDA； 航空瞬變電磁法； 一維正演

0 引言

隨著地球物理勘探工作的深度和廣度不斷提高，數(shù)據(jù)處理結(jié)果的精度要求也越來(lái)越高。巨大的數(shù)據(jù)量和較高的數(shù)據(jù)處理精度對(duì)計(jì)算機(jī)處理技術(shù)提出了越來(lái)越高的要求。提高計(jì)算精度和效率，吸引了更多研究者地注意。近年來(lái)，由于計(jì)算機(jī)圖形處理器(Graphic Processing Unit，GPU)成本和功耗比較低且并行處理能力強(qiáng)大的特點(diǎn)，使它在通用計(jì)算方面得到了極大地發(fā)展[1]。

2007年，NVIDIA公司發(fā)布了一種用于GPU通用運(yùn)算的C擴(kuò)展語(yǔ)言，被稱(chēng)為CUDA(Compute Unified Device Architecture)平臺(tái)。CUDA是一種不需要借助圖形學(xué)API(Application Programming Interface，應(yīng)用程序接口)，就可以使用類(lèi)C語(yǔ)言進(jìn)行GPU通用計(jì)算的開(kāi)發(fā)環(huán)境和軟件體系。

航空瞬變電磁法是一種使用直升機(jī)或固定翼飛機(jī)作為機(jī)載工具，利用物質(zhì)電導(dǎo)率(或電阻率)的差異來(lái)探測(cè)地球電磁響應(yīng)的非常實(shí)用的地球物理探測(cè)方法，由于它擁有快速、高效的掃描廣大區(qū)域的能力，且空間密度大大高于典型的地面電磁法，所以很受歡迎。

筆者編寫(xiě)了基于GPU的航空瞬變電磁法一維正演模擬程序，并且證明了它是高效、可行的。

1 CUDA C程序

CUDA C與C語(yǔ)言非常類(lèi)似[2]，僅僅是增加了少量說(shuō)明符。在用CUDA C編寫(xiě)的程序中， CPU為“主機(jī)(host)”,GPU為“裝置(device)”。在使用GPU時(shí)，我們編寫(xiě)了能在GPU上運(yùn)行的數(shù)據(jù)轉(zhuǎn)換代碼和函數(shù)定義。首先分配GPU上的內(nèi)存，數(shù)據(jù)要在GPU上成功運(yùn)算要傳送到GPU。GPU上的計(jì)算通過(guò)一個(gè)定義了特殊說(shuō)明“__global__”的函數(shù)傳送。這個(gè)函數(shù)叫做“kernel”函數(shù)，然后當(dāng)“主機(jī)”調(diào)用它的時(shí)候，GPU執(zhí)行“kernel”函數(shù)。

CPU與GPU協(xié)同工作，各司其職。CPU負(fù)責(zé)邏輯性強(qiáng)的事務(wù)處理和串行計(jì)算(如數(shù)據(jù)準(zhǔn)備和設(shè)備初始化)；GPU則專(zhuān)注于執(zhí)行高線(xiàn)程化的并行處理。CPU與GPU擁有相互獨(dú)立的存儲(chǔ)器地址空間。一旦確定了程序中的并行部分，我們就可以考慮把這部分計(jì)算工作交給GPU來(lái)執(zhí)行。本程序模型[3]如圖1所示。

圖1 CUDA程序模型Fig.1 CUDA programming model

2 航空瞬變電磁一維正演模型

2.1 一維正演公式推導(dǎo)

正演模型的收發(fā)系統(tǒng)為中心回線(xiàn)裝置。一維地電模型[4]如圖2所示。

圖2 一維地電模型Fig.2 One dimensional geoelectric model

根據(jù)參考文獻(xiàn)[5]，推導(dǎo)出的回線(xiàn)中心處磁場(chǎng)的垂直分量表達(dá)式為：

(1)

式中：I為源電流；a是發(fā)射線(xiàn)圈的半徑；ρ是接收線(xiàn)圈半徑；rTE是水平層狀介質(zhì)的反射系數(shù)；λ是波長(zhǎng)；u0是自由空間的磁導(dǎo)率；ui為大地的磁導(dǎo)率，通常我們?nèi)〈蟮氐拇艑?dǎo)率為自由空間磁導(dǎo)率，ui=u0=4π×10-7H/m；h是發(fā)射線(xiàn)圈的高度；z為接收線(xiàn)圈距地面的距離,在航空電磁法中，z=-h,模型共有n層;hi是每一層的厚度;σi是地層的電導(dǎo)率。

在準(zhǔn)靜態(tài)條件下：

(2)

(3)

(4)

(5)

μ0=λ

(6)

(7)

由頻譜分析理論可以知道，時(shí)間域電磁場(chǎng)與頻率域電磁場(chǎng)通過(guò)傅立葉變換相互關(guān)聯(lián)[6]。

(8)

或

(9)

為了減少橫向的不均勻性，實(shí)際工作中我們一般測(cè)量的是測(cè)量磁場(chǎng)的垂直分量的時(shí)間導(dǎo)數(shù)Hz。

(10)

或

(11)

現(xiàn)在我們得到公式：

(12)

又因?yàn)?/p>

(13)

(14)

顯然，時(shí)間域響應(yīng)是含有貝塞爾函數(shù)的雙重?zé)o窮積分。含有貝塞爾函數(shù)的積分為式(15)。

(15)

其中:Jν是ν階第一類(lèi)貝塞爾函數(shù);實(shí)數(shù)ν>-1。則可以進(jìn)行快速漢克爾變換，成為如下形式[4]。

(16)

式中：ω是快速漢克爾變換系數(shù)。本設(shè)計(jì)中所使用的漢克爾變換系數(shù)是根據(jù)參考文獻(xiàn)[5]求取的，Δ=ln10/10。

筆者發(fā)現(xiàn)FHT(快速漢克爾變換)后的雙重?zé)o窮積分非常適用于并行計(jì)算。為此，編寫(xiě)了一個(gè)CUDA C語(yǔ)言程序，為了提高數(shù)值計(jì)算的穩(wěn)定性，tanh(uihi)使e-uihi用分解運(yùn)算。使用式(17)變換直接分離實(shí)部與虛部。

(17)

2.2 CUDA C程序驗(yàn)證

為驗(yàn)證CUDA C語(yǔ)言程序的正確性和精度，筆者計(jì)算了幾種典型地電模型的瞬變電磁響應(yīng)。

由式(15)可知，發(fā)射電流強(qiáng)度為一個(gè)常數(shù)倍數(shù)，所以本設(shè)計(jì)不在考慮發(fā)射電流的影響，取發(fā)射電流強(qiáng)度I=1A。

2.2.1 不同發(fā)射線(xiàn)圈半徑的瞬變電磁響應(yīng)

在均勻半空間模型的條件下，地層電導(dǎo)率設(shè)為0.1 s/m，飛行高度為30 m，發(fā)射線(xiàn)圈半徑不同的條件下垂直方向航空瞬變電磁Hz響應(yīng)的比較(圖3)。取雙對(duì)數(shù)坐標(biāo)，線(xiàn)圈半徑分別為10 m、50 m、100 m，采樣時(shí)間為0.100 000 0 E-07 s到0.158 500 0 E-01 s取對(duì)數(shù)平均。

由圖3可見(jiàn)，在地電模型相同的情況下，發(fā)射線(xiàn)圈半徑越小，瞬變電磁響應(yīng)的衰減越快，在實(shí)驗(yàn)所用的10 m～50 m內(nèi)，由于迅速衰減，在觀測(cè)時(shí)間內(nèi)，線(xiàn)圈半徑越大響應(yīng)值越大，但在線(xiàn)圈半徑為100 m時(shí)，與50 m的響應(yīng)有交叉，可以推斷，線(xiàn)圈半徑越大，其早期響應(yīng)值越小。在晚期響應(yīng)中，曲線(xiàn)顯然是趨近于平行的。

2.2.2 不同電導(dǎo)率的瞬變電磁響應(yīng)特征

依舊采用均勻半空間模型，飛行高度為30 m，發(fā)射半徑為50 m，模型電導(dǎo)率分別取10 S/m、 0.1 S/m、0.01 S/m。其垂直方向航空瞬變電磁響應(yīng)如圖4所示。采樣時(shí)間為0.100 000 0E-05 s到0.158 500 0E-01 s取對(duì)數(shù)平均。

由圖4可見(jiàn)，地下電導(dǎo)率越小，進(jìn)入晚期的雙對(duì)數(shù)坐標(biāo)下的負(fù)斜率直線(xiàn)越早，在早期觀測(cè)中的響應(yīng)值越大，反之亦然。晚期曲線(xiàn)趨近于平行。

圖3 不同發(fā)射線(xiàn)圈半徑瞬變電磁響應(yīng)特征曲線(xiàn)Fig.3 Airborne transient electromagnetic response characteristic curves for different transmitting coil radius

圖4 不同電導(dǎo)率瞬變電磁響應(yīng)特征曲線(xiàn)Fig.4 Airborne transient electromagnetic response characteristic curves for different conductivity

2.2.3 均勻半空間下不同飛行高度的瞬變電磁響應(yīng)特征

在均勻半空間條件下，發(fā)射線(xiàn)圈半徑為50 m，飛行高度分別取10 m、30 m、100 m，為了使趨勢(shì)更加明顯，我們?cè)黾恿瞬蓸訒r(shí)間，由0.100 000 0E-04 s到0.158 500 0E+01 s取對(duì)數(shù)進(jìn)行采樣，垂直方向的航空瞬變電磁響應(yīng)曲線(xiàn)如圖5。

圖5 不同飛行高度的瞬變電磁響應(yīng)特征曲線(xiàn)圖Fig.5 Airborne transient electromagnetic response characteristic curves for different flight altitude

由圖5可見(jiàn)，不同飛行高度的早期為平行的斜線(xiàn)，但高度越大，早期響應(yīng)值越小，經(jīng)過(guò)中間的漸變過(guò)程后，晚期趨近于同一條負(fù)斜率直線(xiàn)，這與使用相同的地電模型是相對(duì)應(yīng)的。

2.2.4 不同參數(shù)地電模型的瞬變電磁響應(yīng)特征。

圖6是D型地電模型利用CUDA C語(yǔ)言程序的計(jì)算結(jié)果。地層介質(zhì)的導(dǎo)電率分別是0.01 S/m，10 S/m，第一層的厚度是100 m，發(fā)射線(xiàn)圈半徑為50 m，飛行高度30 m。

圖6 D型地電模型航空瞬變電磁響應(yīng)特征曲線(xiàn)圖Fig.6 Airborne transient electromagnetic response characteristic curves for Type D geoelectric model

圖7是三層地電模型利用CUDA C語(yǔ)言程序的計(jì)算結(jié)果。K型地電模型的地層電導(dǎo)率分別為0.1 S/m、0.01 S/m、10 S/m。Q模型的地層電導(dǎo)率分別為0.01 S/m、0.1 S/m、10 S/m。H模型的地層電導(dǎo)率為0.1 S/m、10 S/m、0.01 S/m。地層厚度為100 m,發(fā)射線(xiàn)圈半徑為50 m,飛行高度30 m。

圖7 三層地電模型航空瞬變電磁響應(yīng)特征曲線(xiàn)圖Fig.7 Airborne transient electromagnetic response characteristic curves for three layers geoelectric model

2.3 CUDA C程序與C程序效果對(duì)比

為了比較運(yùn)算性能，我們又用相同的運(yùn)算方法，使用C語(yǔ)言程序運(yùn)算了一次。使用當(dāng)前主流計(jì)算機(jī)平臺(tái)，CPU型號(hào)AMD Athlon II X2 240；GPU型號(hào)Nvidia GeForce 9500 GT，我們?cè)O(shè)計(jì)了以下三個(gè)實(shí)驗(yàn)。

2.3.1 運(yùn)算精度對(duì)比

在確保GPU上的高性能運(yùn)算是精確的同時(shí)，我們構(gòu)建一個(gè)均勻半空間模型。它的導(dǎo)電率為0.1 S/m.發(fā)射線(xiàn)圈半徑為50 m，飛行高度為30 m。圖8為由C語(yǔ)言源代碼雙精度浮點(diǎn)數(shù)計(jì)算值與CUDA C的單精度浮點(diǎn)數(shù)計(jì)算值比較。

從圖8中，我們可以看到兩程序數(shù)據(jù)的結(jié)果是一致的，CUDA C程序的精確度是可以信賴(lài)的。

圖8 均勻半空間航空瞬變電磁響應(yīng)數(shù)值模擬曲線(xiàn)GPU與CPU對(duì)比圖Fig.8 Comparison chart the airborne transient electromagnetic response characteristic on GPU and CPU

首先定義相對(duì)誤差，

GPUS為均勻半空間模型在使用GPU運(yùn)算的響應(yīng)，CPUS為均勻半空間模型使用CPU運(yùn)算的響應(yīng)，則相對(duì)誤差曲線(xiàn)如圖9所示。

圖9 均勻半空間航空瞬變電磁響應(yīng)數(shù)值模擬GPU與CPU運(yùn)算相對(duì)誤差曲線(xiàn)Fig.9 The relative error curves of GPU and CPU in the numerical simulation of airborne transient electromagnetic response in homogeneous half-space model

可見(jiàn)，在均勻半空間比較計(jì)算比較穩(wěn)定的情況下，計(jì)算誤差在10-6左右，且相對(duì)比較穩(wěn)定，因此計(jì)算精度滿(mǎn)足要求。

2.3.2 運(yùn)算速度對(duì)比

(1)不同地電模型的運(yùn)行時(shí)間。在正演運(yùn)算中，地電模型層數(shù)與運(yùn)算時(shí)間直接相關(guān)。僅通過(guò)改變地電模型層數(shù)，我們得到了兩程序執(zhí)行時(shí)間如表1所示。Z表示地電模型層數(shù)，GT表示GPU執(zhí)行時(shí)間，CT表示CPU執(zhí)行時(shí)間，圖10所示為GPU和CPU代碼執(zhí)行時(shí)間的對(duì)比結(jié)果。事實(shí)上，結(jié)論非常明顯，隨著地電模型層數(shù)的增加，GPU執(zhí)行時(shí)間相對(duì)CPU進(jìn)步一縮短了。

圖10 GPU與CPU代碼執(zhí)行時(shí)間對(duì)比圖Fig.10 Execution time comparison chart about CPU program and GPU program

地層數(shù)ZGPU耗時(shí)GT/msCPU耗時(shí)CT/msCT/GT278132817.036125418733.5011171782845.78162191151552.58212651507856.90

2)相同地電模型的執(zhí)行時(shí)間。在航空瞬變電磁法一維正演模擬中，GPU程序和CPU程序執(zhí)行時(shí)間都非常短。執(zhí)行時(shí)間短會(huì)造成測(cè)量誤差，為了減少觀測(cè)誤差，我們做了以下改進(jìn)。在函數(shù)主體語(yǔ)句前加一個(gè)代碼：for(int dd=0;dd

表2 實(shí)驗(yàn)2執(zhí)行時(shí)間數(shù)據(jù)記錄

3 結(jié)論

筆者提出了基于GPU高性能運(yùn)算的CUDA C代碼編程的航空瞬變電磁法正演模擬法，通過(guò)模型實(shí)驗(yàn)，我們確認(rèn)基于GPU的CUDA C代碼所編的程序相比較C語(yǔ)言編程更有優(yōu)勢(shì)，同時(shí)，數(shù)據(jù)量越大，GPU耗時(shí)明顯越少，總而言之，GPU的并行計(jì)算能力使它在電磁法數(shù)值計(jì)算方面應(yīng)用前景廣闊。

[1] 丁科，譚營(yíng).GPU通用計(jì)算極其在計(jì)算智能領(lǐng)域的應(yīng)用[J].智能系統(tǒng)學(xué)報(bào)，2015，10(1)：1-11. DING K，TAN Y． A review on general purpose computing on GPU and its applications in computational intelligence[J]．CAAI Transactions on Intelligent Systems，2015，10( 1) : 1-11．(In Chinese)

[2] TAKEHIKO NAWATA，REIJI SUDA. APTCC:Auto Parallelizing Translator From C To CUDA[J].Procedia Computer Science,2011(4):352-361.

[3] 張舒，褚艷利.GPU高性能運(yùn)算之CUDA[M].北京：中國(guó)水利水電出版社，2009. ZHANG S,CHU Y L. CUDA GPU high-performance computing[M].Beijing:China Water & Power Press,2009.(In Chinese)

[4] 朱凱光，林君，劉長(zhǎng)勝，等.頻率域航空電磁法一維正演與探測(cè)深度[J].地球物理學(xué)進(jìn)展，2008，6(23)：1943-1946. ZHU K G，LIN J，LIU C S，et al. One-dimensional forwand and prospecting depth for airborne frequency domain electromagnetic method [J]. Progress in Geophysics，2008，6(23)：1943-1946. (In Chinese)

[5] 米薩克·N·納比吉安.勘察地球物理電磁法第一卷理論[M].趙經(jīng)祥，譯.北京：地質(zhì)出版社，199. Misac N Nabighian.Explore geophysics-electromagnetic method[M].Translated by ZHAO jingxiang.Beijing:Geological Publishing House,1992.(In Chinese)

[6] 劉繼東，林長(zhǎng)佑，楊長(zhǎng)福.瞬變電磁資料反演解釋方法[M].蘭州：蘭州大學(xué)出版社，2001 LIU J D，LIN Y, YANG C F. Transient electromagnetic information inversion method[M]. Lanzhou: Lanzhou University Press, 2001. (In Chinese)

[7] 牛之璉.時(shí)間域電磁法[M].長(zhǎng)沙：中南大學(xué)出版社，2007. NIU Z L. Time domain electromagnetic method[M]. Changsha: Central South University Press, 2007. (In Chinese)

[8] 樸化榮.電磁測(cè)深法原理[M].北京：地質(zhì)出版社，1990. PIAO H R. Principle of electromagnetic soundingmethod [M].Beijing:Geological publishing, 1990. (In Chinese)

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

Airborne TEM response simulation based on GPU

SHEN Lei1，DONG Chunhui1,2

(1. School of highway and architecture, Shandong Transport Vocational College, Weifang 261206,China； 2. Institute of Highway and Bridge Engineering,Dalian Maritime University，Dalian 116026,China)

An efficient algorithm for time-domain airborne EM forward simulation is presented. It is based on CUDA (compute unified device architecture) that was released by NVIDIA in 2007. We derivate the time domain airborne EM formula and then write the code with GPU language accordingly. The modeled results demonstrate that the parallel processing capabilities of graphics processors (GPUs). This significant improvement in performance is gained compared to the corresponding CPU-based solver, while maintaining the numerical accuracy. To the best of the authors' knowledge, this is the fastest time-domain solver for modeling the AEM simulation. It would have potential for complex geophysics model calculation and inversion.

GPU; CUDA; airbome transient electromagnetic method; one-dimensional forward

2016-03-18 改回日期：2016-03-30

沈磊(1978-)，男，講師，主要從事計(jì)算機(jī)軟件應(yīng)用、工程巖土方面的教學(xué)與研究工作，E-mail：shenlei116@163.com。

1001-1749(2017)02-0170-06

P 631.2

A

10.3969/j.issn.1001-1749.2017.02.03