沈 磊， 董春暉,2

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

基于GPU的航空瞬變電磁法響應模擬

沈 磊1， 董春暉1,2

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

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

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

0 引言

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

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

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

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

1 CUDA C程序

CUDA C與C語言非常類似[2]，僅僅是增加了少量說明符。在用CUDA C編寫的程序中， CPU為“主機(host)”,GPU為“裝置(device)”。在使用GPU時，我們編寫了能在GPU上運行的數據轉換代碼和函數定義。首先分配GPU上的內存，數據要在GPU上成功運算要傳送到GPU。GPU上的計算通過一個定義了特殊說明“__global__”的函數傳送。這個函數叫做“kernel”函數，然后當“主機”調用它的時候，GPU執(zhí)行“kernel”函數。

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

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

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

2.1 一維正演公式推導

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

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

根據參考文獻[5]，推導出的回線中心處磁場的垂直分量表達式為：

(1)

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

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

(2)

(3)

(4)

(5)

μ0=λ

(6)

(7)

由頻譜分析理論可以知道，時間域電磁場與頻率域電磁場通過傅立葉變換相互關聯[6]。

(8)

或

(9)

為了減少橫向的不均勻性，實際工作中我們一般測量的是測量磁場的垂直分量的時間導數Hz。

(10)

或

(11)

現在我們得到公式：

(12)

又因為

(13)

(14)

顯然，時間域響應是含有貝塞爾函數的雙重無窮積分。含有貝塞爾函數的積分為式(15)。

(15)

其中:Jν是ν階第一類貝塞爾函數;實數ν>-1。則可以進行快速漢克爾變換，成為如下形式[4]。

(16)

式中：ω是快速漢克爾變換系數。本設計中所使用的漢克爾變換系數是根據參考文獻[5]求取的，Δ=ln10/10。

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

(17)

2.2 CUDA C程序驗證

為驗證CUDA C語言程序的正確性和精度，筆者計算了幾種典型地電模型的瞬變電磁響應。

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

2.2.1 不同發(fā)射線圈半徑的瞬變電磁響應

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

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

2.2.2 不同電導率的瞬變電磁響應特征

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

由圖4可見，地下電導率越小，進入晚期的雙對數坐標下的負斜率直線越早，在早期觀測中的響應值越大，反之亦然。晚期曲線趨近于平行。

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

圖4 不同電導率瞬變電磁響應特征曲線Fig.4 Airborne transient electromagnetic response characteristic curves for different conductivity

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

在均勻半空間條件下，發(fā)射線圈半徑為50 m，飛行高度分別取10 m、30 m、100 m，為了使趨勢更加明顯，我們增加了采樣時間，由0.100 000 0E-04 s到0.158 500 0E+01 s取對數進行采樣，垂直方向的航空瞬變電磁響應曲線如圖5。

圖5 不同飛行高度的瞬變電磁響應特征曲線圖Fig.5 Airborne transient electromagnetic response characteristic curves for different flight altitude

由圖5可見，不同飛行高度的早期為平行的斜線，但高度越大，早期響應值越小，經過中間的漸變過程后，晚期趨近于同一條負斜率直線，這與使用相同的地電模型是相對應的。

2.2.4 不同參數地電模型的瞬變電磁響應特征。

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

圖6 D型地電模型航空瞬變電磁響應特征曲線圖Fig.6 Airborne transient electromagnetic response characteristic curves for Type D geoelectric model

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

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

2.3 CUDA C程序與C程序效果對比

為了比較運算性能，我們又用相同的運算方法，使用C語言程序運算了一次。使用當前主流計算機平臺，CPU型號AMD Athlon II X2 240；GPU型號Nvidia GeForce 9500 GT，我們設計了以下三個實驗。

2.3.1 運算精度對比

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

從圖8中，我們可以看到兩程序數據的結果是一致的，CUDA C程序的精確度是可以信賴的。

圖8 均勻半空間航空瞬變電磁響應數值模擬曲線GPU與CPU對比圖Fig.8 Comparison chart the airborne transient electromagnetic response characteristic on GPU and CPU

首先定義相對誤差，

GPUS為均勻半空間模型在使用GPU運算的響應，CPUS為均勻半空間模型使用CPU運算的響應，則相對誤差曲線如圖9所示。

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

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

2.3.2 運算速度對比

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

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

地層數ZGPU耗時GT/msCPU耗時CT/msCT/GT278132817.036125418733.5011171782845.78162191151552.58212651507856.90

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

表2 實驗2執(zhí)行時間數據記錄

3 結論

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

[1] 丁科，譚營.GPU通用計算極其在計算智能領域的應用[J].智能系統(tǒng)學報，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高性能運算之CUDA[M].北京：中國水利水電出版社，2009. ZHANG S,CHU Y L. CUDA GPU high-performance computing[M].Beijing:China Water & Power Press,2009.(In Chinese)

[4] 朱凱光，林君，劉長勝，等.頻率域航空電磁法一維正演與探測深度[J].地球物理學進展，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].趙經祥，譯.北京：地質出版社，199. Misac N Nabighian.Explore geophysics-electromagnetic method[M].Translated by ZHAO jingxiang.Beijing:Geological Publishing House,1992.(In Chinese)

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

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

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

第39卷 第2期2017年3月物探化探計算技術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-)，男，講師，主要從事計算機軟件應用、工程巖土方面的教學與研究工作，E-mail：shenlei116@163.com。

1001-1749(2017)02-0170-06

P 631.2

A

10.3969/j.issn.1001-1749.2017.02.03