戴云峰，周志芳，強(qiáng)建科，劉 冰

（1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098；2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083）

基于GPU實(shí)現(xiàn)漢克爾變換并行計(jì)算

戴云峰1，周志芳1，強(qiáng)建科2，劉 冰2

（1.河海大學(xué) 地球科學(xué)與工程學(xué)院，江蘇 南京 210098；2.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙 410083）

地球物理勘探技術(shù)日新月異，地球物理勘探數(shù)據(jù)的處理和解釋對(duì)高性能計(jì)算機(jī)的要求越來越高。相比于地震勘探，重力、磁法、電法勘探中的并行計(jì)算研究還都處于起步階段。基于GPU的并行計(jì)算能夠提供強(qiáng)大的計(jì)算能力和存儲(chǔ)器帶寬，同時(shí)具有良好的可編程性、較低的成本和較短的開發(fā)周期。這里實(shí)現(xiàn)了瞬變電磁法一維正演計(jì)算中漢克爾變換基于GPU的并行計(jì)算，比較了漢克爾變換串行算法和并行算法的計(jì)算耗時(shí)，基于GPU技術(shù)的并行計(jì)算相比串行計(jì)算，獲得了很高的加速比。

圖形處理器；漢克爾變換；并行計(jì)算

0 前言

計(jì)算機(jī)科學(xué)技術(shù)的高速發(fā)展和大型工程計(jì)算的需求，使得并行計(jì)算（Parallel Computing）技術(shù)已經(jīng)成為目前乃至今后工程計(jì)算領(lǐng)域中一項(xiàng)重要技術(shù)。并行計(jì)算是相對(duì)于串行技術(shù)而言的，并行計(jì)算可以分為時(shí)間和空間上的計(jì)算，時(shí)間上的并行就是指流水線技術(shù)，空間上的并行則是指用多個(gè)處理器并發(fā)地執(zhí)行計(jì)算。

隨著地球物理勘探工作的深度和廣度不斷提高，數(shù)據(jù)處理結(jié)果的精度要求也越來越高。高性能計(jì)算技術(shù)對(duì)于地球物理勘察數(shù)據(jù)處理有著重要的推動(dòng)作用，這是因?yàn)椋?］：①地球物理勘察的數(shù)據(jù)量巨大；②地球物理勘察數(shù)據(jù)中蘊(yùn)涵著大量可并行計(jì)算成份。

計(jì)算機(jī)圖形處理器（Graphics Processing U－nit，GPU，俗稱顯卡）的技術(shù)定義是“具有集成轉(zhuǎn)換、照明、三角形設(shè)置／剪切、渲染引擎并且每秒至少可處理一千萬個(gè)多邊形的單芯片處理器”。GPU通用計(jì)算技術(shù)雖然在石油勘探領(lǐng)域中略有使用，但是還不能滿足石油工業(yè)的發(fā)展需要。地震勘探為了獲得更加精確的成像效果，由二維勘探向三維勘探發(fā)展，以及由疊后處理向疊前處理發(fā)展，隨之增加的數(shù)據(jù)處理量向高性能計(jì)算技術(shù)提出了更高的要求。GPU計(jì)算技術(shù)已經(jīng)用于加速三維有限元地震波數(shù)值模擬以及地震屬性的提取，鐘勇［2］應(yīng)用GPU加速計(jì)算提取三維地震圖像的結(jié)構(gòu)張量，由結(jié)構(gòu)張量導(dǎo)向進(jìn)行地震圖像擴(kuò)散濾波增強(qiáng)預(yù)處理；劉紅偉等［3］針對(duì)疊前逆時(shí)偏移計(jì)算量大的問題，使用GPU實(shí)現(xiàn)算法加速，比傳統(tǒng)的CPU計(jì)算速度提高了一個(gè)數(shù)量級(jí)。GPU計(jì)算技術(shù)在電法勘探中的使用較少，電磁法數(shù)據(jù)處理中包含了很多大規(guī)模矩陣的計(jì)算，涉及到矩陣乘法和線性方程求解，計(jì)算量較大。大地電磁正反演研究的熱點(diǎn)開始向三維問題轉(zhuǎn)移，交錯(cuò)網(wǎng)格法等已成為大地電磁三維正演計(jì)算的主要方法。直流電法和大地電磁法的三維正演，都存在需要求解大型稀疏對(duì)稱系數(shù)矩陣線性方程組的問題，張帆［4］使用基于GPU并行計(jì)算和基于MPI的并行計(jì)算相結(jié)合的計(jì)算模式，對(duì)這兩種電磁法已有的三維正演串行算法進(jìn)行并行化處理，提高了運(yùn)算效率。

作者在本文中介紹了圖形處理器GPU并行計(jì)算的原理，以及CUDA編程模型，實(shí)現(xiàn)了漢克爾變換基于GPU的并行計(jì)算，比較了漢克爾變換串行計(jì)算和并行計(jì)算的計(jì)算耗時(shí)，探討GPU在瞬變電磁法一維正演計(jì)算中的應(yīng)用前景。

1 圖形處理器GPU并行計(jì)算的原理

1.1 GPU與CPU結(jié)構(gòu)比較

計(jì)算機(jī)中的處理器包括CPU（Central Processing Unit，中央處理器）和GPU圖形處理器（Graphic Processing Unit，圖形處理器）。CPU提高單個(gè)芯片性能的主要傳統(tǒng)手段是提高處理器的工作頻率，以及增加指令級(jí)并行?；贕PU的通用計(jì)算技術(shù)，它主要應(yīng)用圖形處理芯片（GPU）的高性能并行處理能力實(shí)現(xiàn)科學(xué)計(jì)算。與CPU相比，GPU的設(shè)計(jì)中能使更多晶體管用于數(shù)據(jù)處理（ALU，Arithmetic Logic Unit即算術(shù)邏輯單元），而非用于數(shù)據(jù)緩存（DRAM）和流控制（Control），見圖1。

圖1 CPU與GPU晶體管使用方式的對(duì)比Fig.1 Comparison of transistors used in CPU and GPU

GPU專門用于解決可表示為數(shù)據(jù)并行計(jì)算的問題，在許多數(shù)據(jù)元素上并行執(zhí)行的程序，具有極高的計(jì)算密度（數(shù)學(xué)運(yùn)算與存儲(chǔ)器運(yùn)算的比率）。由于所有數(shù)據(jù)元素都執(zhí)行相同的程序，因此對(duì)精密流控制的要求不高。由于在許多數(shù)據(jù)元素上運(yùn)行，且具有較高的計(jì)算密度，因而可通過計(jì)算隱藏存儲(chǔ)器訪問延遲，而不必使用較大的數(shù)據(jù)緩存。數(shù)據(jù)并行處理會(huì)將數(shù)據(jù)元素映射到并行處理線程，許多處理大型數(shù)據(jù)集的應(yīng)用程序都可使用數(shù)據(jù)并行編程模型來加速計(jì)算。

1.2 基于GPU的CUDA編程模型

CUDA（Compute Unified Device Architecture，統(tǒng)一計(jì)算設(shè)備架構(gòu)）是一種新型的硬件和軟件架構(gòu)，用于將GPU作為數(shù)據(jù)并行計(jì)算設(shè)備并在GPU上進(jìn)行計(jì)算的發(fā)放和管理，而無需將其映射到圖像 API（Application Programming Interface，應(yīng)用程序編程接口）。隨著以CUDA為代表的GPU通用計(jì)算API的普及，GPU的含義有可能從圖形處理器（Graphic Processing Unit）擴(kuò)展為通用處理器（General Purpose Unit）。GPU擅長的是圖形類的或者是非圖形類的高度并行數(shù)值計(jì)算，GPU可以容納成百上千條沒有邏輯關(guān)系的數(shù)值計(jì)算線程，它的優(yōu)勢(shì)是無邏輯關(guān)系數(shù)據(jù)的并行計(jì)算。GPU數(shù)值計(jì)算的優(yōu)勢(shì)主要是浮點(diǎn)運(yùn)算，它執(zhí)行浮點(diǎn)運(yùn)算快是靠大量并行，但是這種數(shù)值運(yùn)算的并行性在面對(duì)邏輯判斷執(zhí)行時(shí)卻發(fā)揮不了作用。

CUDA編程模型將CPU作為主機(jī)（Host），GPU作為協(xié)處理器（co－processor）或者設(shè)備（Device）。在一個(gè)系統(tǒng)中可以存在一個(gè)主機(jī)和若干個(gè)設(shè)備。GPU計(jì)算的模式就是，在異構(gòu)協(xié)同處理計(jì)算模型中，將CPU與GPU結(jié)合起來加以利用，CPU負(fù)責(zé)進(jìn)行邏輯性強(qiáng)的事務(wù)處理和串行計(jì)算，GPU則擅長于執(zhí)行高度并行化的線程任務(wù)。CUDA假設(shè)CUDA線程可在物理獨(dú)立的設(shè)備上執(zhí)行，此類設(shè)備作為運(yùn)行C語言程序的主機(jī)協(xié)處理器操作，見圖2。例如，當(dāng)內(nèi)核在GPU上執(zhí)行，而C語言程序的其它部份在CPU上執(zhí)行。

圖2 異構(gòu)編程（串行代碼在主機(jī)上執(zhí)行，并行代碼在設(shè)備上執(zhí)行）Fig.2 Heterogeneous programming（serial code running on the host，parallel code executes on the device）

在CUDA編程模型中，將threadIdx設(shè)置為一個(gè)包含三個(gè)元素的向量，因而可使用一維、二維或三維索引標(biāo)識(shí)線程，構(gòu)成一維、二維或三維線程塊。在NVIDIA Tesla架構(gòu)中，一個(gè)線程塊最多可以擁有512條線程。但一個(gè)內(nèi)核可能由多個(gè)大小相同的線程塊執(zhí)行，因而線程總數(shù)應(yīng)等于每個(gè)塊的線程數(shù)乘以塊的數(shù)量。這些塊將組織為一個(gè)一維或二維線程塊網(wǎng)格，見圖3。

圖3 線程塊網(wǎng)絡(luò)Fig.3 Thread blocks network

2 基于GPU的漢克爾變換并行計(jì)算

目前，基于GPU計(jì)算的并行算法主要來源于兩個(gè)方面：①將串行算法改寫成并行算法；②從數(shù)學(xué)模型出發(fā)直接構(gòu)造并行算法。作者在本文根據(jù)瞬變電磁法一維正演計(jì)算的算法，用C語言編制了正演的串行計(jì)算程序，并研究使用CUDA的并行算法對(duì)串行算法進(jìn)行改寫，以提高數(shù)據(jù)處理的速度。在windows操作系統(tǒng)下，本文作者所有程序都是在Visual Studio C＋＋2005環(huán)境下運(yùn)行。實(shí)驗(yàn)使用的是NVIDIA Geforce 8500GT顯卡，它的計(jì)算能力為1.1。

2.1 漢克爾變換表達(dá)式

水平層狀介質(zhì)上瞬變電磁一維正演公式可表示為一個(gè)雙重積分：內(nèi)層含有零階和一階貝塞爾函數(shù)的漢克爾型積分，外層為正弦或余弦積分。在積分公式中，由于貝塞爾函數(shù)隨自變量緩慢地震蕩衰減，所以不采用常規(guī)的數(shù)值積分方法計(jì)算，采用數(shù)字濾波法具有良好的效果［5］。

漢克爾變換的一般表達(dá)式［6］：

其中 Ji是i階第一類貝塞爾函數(shù)；K（m）是積分變換函數(shù)；m是被積變量。

數(shù)字濾波器有一組濾波系數(shù)，兩個(gè)常數(shù)：位移a和抽樣間隔s。位移a決定了輸入函數(shù)抽樣起始點(diǎn)，間隔s決定抽樣時(shí)的時(shí)間間隔。由

可得出

其中 mi＝ （1／r）＊10［a＋（i－1）s］（i＝1、2、…、n）；Wi為濾波器系數(shù)。

2.2 漢克爾變換串行算法性能分析及并行算法實(shí)現(xiàn)

零階漢克爾變換的串行程序，使用clock＿t和clock（）函數(shù)來獲得程序各部份的運(yùn)行時(shí)間。計(jì)算所耗時(shí)間與計(jì)算過程中所取時(shí)間點(diǎn)個(gè)數(shù)的關(guān)系見圖4，橫坐標(biāo)為所取的時(shí)間點(diǎn)數(shù)，縱坐標(biāo)為對(duì)應(yīng)的計(jì)算耗時(shí)（單位：S），分別測(cè)得200、300、500、700、900、1 000個(gè)時(shí)間點(diǎn)的程序計(jì)算耗時(shí)。在測(cè)試運(yùn)行單元時(shí)間時(shí)，先清除上一次運(yùn)行的解決方案，然后再生成新的解決方案。

當(dāng)時(shí)間點(diǎn)的個(gè)數(shù)小于500點(diǎn)時(shí)，即使測(cè)試的是整個(gè)程序運(yùn)行的時(shí)間，由于運(yùn)行速度很快，加上計(jì)時(shí)函數(shù)的精度限制，所得的計(jì)算耗時(shí)幾乎為“0”值，見圖4。在這種情況下，使用并行的算法來實(shí)現(xiàn)漢克爾變換反而得不償失。根據(jù)CUDA實(shí)現(xiàn)效果原則，計(jì)算量太小使用CUDA不劃算。作者在對(duì)漢克爾變換的串行算法與并行算法進(jìn)行計(jì)算性能對(duì)比時(shí)，所取時(shí)間點(diǎn)個(gè)數(shù)為512（單個(gè)線程塊最多擁有的線程數(shù)）。

作者在本文中進(jìn)行了漢克爾變換的并行計(jì)算，計(jì)算出N點(diǎn)的零階漢克爾變換（一階變換原理相同，不再重復(fù)）。使用取對(duì)數(shù)間隔函數(shù)，獲得的對(duì)數(shù)間隔時(shí)間點(diǎn)，即計(jì)算式中r的值，以數(shù)組的形式作為漢克爾變換計(jì)算的輸入。對(duì)于單個(gè)的時(shí)間點(diǎn)計(jì)算涉及累加的操作，但各個(gè)時(shí)間點(diǎn)之間是沒有影響的，故具有較好的并行性。因?yàn)橹皇沁M(jìn)行漢克爾變換并行計(jì)算，數(shù)據(jù)量不大，所以只使用了一個(gè)線程塊，每個(gè)線程塊中定義與時(shí)間點(diǎn)相同個(gè)數(shù)的線程（每個(gè)線程塊Block最多可以定義512條線程Thread），每個(gè)線程負(fù)責(zé)進(jìn)行一個(gè)時(shí)間點(diǎn)的漢克爾變換。漢克爾變換的系數(shù)為單精度，計(jì)算結(jié)果和串行計(jì)算漢克爾變換的結(jié)果相同。數(shù)值計(jì)算的精度取決于積分區(qū)間的長度n，抽樣點(diǎn)的位置λi和加權(quán)系數(shù)Wi。

所取時(shí)間點(diǎn)數(shù)為512的漢克爾變換串行計(jì)算與并行計(jì)算性能的對(duì)比，時(shí)間單位：ms，見圖5。串行計(jì)算用時(shí)為15ms，而并行內(nèi)核函數(shù)計(jì)算的用時(shí)為0.027 937ms，即使包含計(jì)算結(jié)果從顯存輸出到內(nèi)存中計(jì)算耗時(shí)（圖中并行計(jì)算1代表）也僅僅只有0.525 206ms。

圖5 漢克爾變換串行與并行運(yùn)行耗時(shí)對(duì)比（單位：ms）Fig.5 Time－consuming comparison of Hankel transform serial and parallel operation（unit：ms）

表1 漢克爾變換并行計(jì)算相對(duì)串行計(jì)算加速比Tab.1 Hankel transform parallel computing speedup relative serial

由此可見，對(duì)于數(shù)據(jù)量比較大的漢克爾變換，使用并行計(jì)算可以獲得很大的性能提高，見表1。但是并行計(jì)算的過程相對(duì)于串行算法的計(jì)算過程略顯繁瑣，而且CUDA的編程模型很大程度上依賴于硬件的模型，特別是CUDA存儲(chǔ)器使用的限制，這就給編程帶來了一些麻煩。鑒于設(shè)備內(nèi)寄存器和共享存儲(chǔ)器的容量有限，在進(jìn)行大型科學(xué)計(jì)算的過程中一定要考慮內(nèi)存的溢出。GPU單精度性能遠(yuǎn)遠(yuǎn)超過雙精度，作者在本文中使用過的雙精度需要截?cái)喑蓡尉?，在程序運(yùn)行時(shí)出現(xiàn)了一些警告，但這并不影響程序的運(yùn)行，可能計(jì)算精度會(huì)稍微差一點(diǎn)。使用共享存儲(chǔ)器使得漢克爾變換計(jì)算過程中線程間的通信延遲最小，性能大幅度提高。

3 結(jié)論

（1）在水平層狀介質(zhì)上瞬變電磁測(cè)深一維正演計(jì)算中，對(duì)于漢克爾變換的計(jì)算，基于GPU技術(shù)的并行計(jì)算相比串行計(jì)算獲得很高的加速比。并行計(jì)算的過程相對(duì)于串行算法的計(jì)算過程略顯繁瑣，而且CUDA的編程模型很大程度上依賴于硬件的模型，但是GPU在并行計(jì)算方面的優(yōu)勢(shì)以及超強(qiáng)的浮點(diǎn)運(yùn)算性能可以大大縮短計(jì)算時(shí)間。作者在本文中僅僅實(shí)現(xiàn)了漢克爾變換不同時(shí)間點(diǎn)的并行計(jì)算，僅僅是整個(gè)正演程序的一部份，對(duì)于整個(gè)串行正演程序更深入的并行算法有待進(jìn)一步研究。

（2）相比于地震勘探，重力、磁法、電法勘探中的并行計(jì)算研究還都處于起步階段，我們需要基于GPU計(jì)算平臺(tái)開發(fā)更多應(yīng)用于地球物理方面的并行計(jì)算軟件。研究并行算法和編程模型，在地球物理勘探數(shù)據(jù)處理中具有廣闊的前景。

［1］ 劉羽，王家映.地球物理數(shù)據(jù)處理與并行計(jì)算［J］.桂林工學(xué)院學(xué)報(bào)，2004，24（4）：412.

［2］ 鐘勇，陳磊.基于GPU計(jì)算的三維地震斷層解釋［J］.石油物探，2011，50（2）：160.

［3］ 劉紅偉，李博，劉洪，等.地震疊前逆時(shí)偏移高階有限差分算法及 GPU實(shí)現(xiàn)［J］.地球物理學(xué)報(bào)，2010，53（7）：1725.

［4］ 張帆.基于MPI和GPU直流電法和大地電磁法三維正演的并行算法研究［D］.北京：中國地質(zhì)大學(xué)，2011.

［5］ 唐寶山.瞬變電磁法中心回線裝置一維正反演研究［D］.北京：中國地質(zhì)大學(xué)，2008.

［6］ 劉桂芬.回線源層狀大地航空瞬變電磁場(chǎng)的理論計(jì)算［D］.吉林：吉林大學(xué)，2008.

［7］ 趙改善.地球物理高性能計(jì)算的新選擇：GPU計(jì)算技術(shù)［J］.勘探地球物理進(jìn)展，2007，30（5）：399.

［8］ 趙改善.可重構(gòu)計(jì)算技術(shù)及其在地球物理中的應(yīng)用前景［J］.勘探地球物理進(jìn)展，2007，30（4）：309.

［9］ HE CHUAN，SUN CHUANWEN，LU MI，et al.Prestack Kirchhoff time migration on hight performance reconfigurable computing platform［J］.Expanded Abstract of 75th Annual international SEG Meeting ，2005：902.

［10］HE CHUAN，QIN GUAN，ZHAO WEI.Highorder finite difference modeling on reconfigurable computing platform［J］.Expanded Abstract of 75th Annual international SEG Meeting 2005：755.

［11］張舒.模式識(shí)別并行算法與GPU高速實(shí)現(xiàn)研究［D］.成都：電子科技大學(xué)，2009.

［12］趙國澤，陳小斌，湯吉.中國地球電磁法新進(jìn)展和發(fā)展趨勢(shì)［J］.地球物理學(xué)進(jìn)展，2007，22（4）：1171.

［13］NVIDIA CUDA計(jì)算統(tǒng)一設(shè)備架構(gòu)編程指南2.0［S］.

［14］羅延鐘，張勝業(yè)，王衛(wèi)平.時(shí)間域航空電磁法一維正演研究［J］.地球物理學(xué)報(bào)，2003，46（5）：719.

［15］KNIGHT JH，RAICH AP.Transient electro－magnetic calculations using the Gaver－Stehfest inverse Laplace transform method［J］.Geophysics，1982，47（1）：47.

TP 317.4

A

10.3969／j.issn.1001－1749.2012.05.20

1001—1749（2012）05—0614—05

水利部公益性行業(yè)科研專項(xiàng)經(jīng)費(fèi)項(xiàng)目資助（201001020）

2011－11－24 改回日期：2012－04－02

戴云峰（1988－），男，碩士，從事水文地質(zhì)方面研究。