基于GPU并行的重力、重力梯度三維正演快速計(jì)算及反演策略

陳召曦，孟小紅，郭良輝，劉國(guó)峰

1 地質(zhì)過程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，中國(guó)地質(zhì)大學(xué)，北京 100083

2 地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室，中國(guó)地質(zhì)大學(xué)，北京 100083

3 地球物理與信息技術(shù)學(xué)院，中國(guó)地質(zhì)大學(xué)，北京 100083

1 引 言

隨著重力測(cè)量設(shè)備的不斷發(fā)展，重力勘探方法已經(jīng)廣泛應(yīng)用于油氣勘探、金屬礦勘查、地球深部構(gòu)造及地殼研究、劃分大地構(gòu)造單元及工程與環(huán)境地球物理勘探等領(lǐng)域［1-2］.特別是近幾年來(lái)，具有更高分辨率的航空、海洋重力梯度測(cè)量?jī)x器與配套技術(shù)的快速發(fā)展，使得重力勘探方法也愈來(lái)愈引起地球物理學(xué)界的重視.相應(yīng)地，重力（重力梯度）數(shù)據(jù)量也愈來(lái)愈成為海量級(jí)別.針對(duì)海量重力（重力梯度數(shù)據(jù)）實(shí)際資料快速估算解釋的要求，國(guó)內(nèi)外學(xué)者已經(jīng)開展了很多工作，特別是帶有不同約束條件的物性反演研究［3-21］.物性反演是將觀測(cè)區(qū)域?qū)?yīng)的地下半空間離散化成規(guī)則的長(zhǎng)方體單元，通過反演方法確定各離散單元的物性值，由物性的分布確定場(chǎng)源的實(shí)際分布情況.近幾年來(lái)隨著計(jì)算水平的提高，這種方法已逐漸成為國(guó)內(nèi)外重力數(shù)據(jù)反演的主要方向［22］.

在基于物性模型的重力（重力梯度）數(shù)據(jù)約束反演或者聯(lián)合反演中，模型體的正演計(jì)算是反演過程的重要組成部分.不同規(guī)模模型體的正演計(jì)算效率是主要制約反演充分發(fā)揮作用的瓶頸問題.將網(wǎng)格觀測(cè)單元與地半空間模型單元中心一一對(duì)應(yīng)，采用等效格架技術(shù)，只需計(jì)算一遍模型正演，并且只計(jì)算每層的一個(gè)網(wǎng)格單元，使模型正演計(jì)算幾乎消失［22-23］.上述方法只適用于模型的剖分與網(wǎng)格采取某種對(duì)應(yīng)關(guān)系時(shí).基于小波變換，采用系數(shù)矩陣壓縮與重構(gòu)的方式，減少存儲(chǔ)空間，提高了反演效率.但這種方法損失信號(hào)［8，24］.

另外一種思路是采用并行計(jì)算技術(shù).以往大規(guī)模數(shù)據(jù)并行運(yùn)算只局限于專業(yè)以及昂貴的并行計(jì)算機(jī)上.如今隨著GPU通用計(jì)算的實(shí)現(xiàn)，這一現(xiàn)狀得到改變.GPU在處理能力和存儲(chǔ)器帶寬相對(duì)CPU有明顯優(yōu)勢(shì)，在成本和功耗上也不需要付出太大代價(jià)，從而為相同的大規(guī)模數(shù)據(jù)運(yùn)算問題提供了新的解決方案［25-29］.鑒于正演計(jì)算程序具有很好的并行改進(jìn)能力，本文利用NIVIDIA CUDA編程平臺(tái)，詳細(xì)分析了并行計(jì)算中的幾個(gè)關(guān)鍵問題（代碼設(shè)計(jì)、存儲(chǔ)器優(yōu)化、精度問題），從而實(shí)現(xiàn)了基于GPU并行的重力、重力梯度三維快速正演計(jì)算方法.試驗(yàn)結(jié)果表明，該并行技術(shù)計(jì)算結(jié)果和單線程計(jì)算結(jié)果一致，可大大提高效率.最后討論了基于GPU并行計(jì)算的重力、重力梯度數(shù)據(jù)三維反演策略.

2 正演理論

常規(guī)重力測(cè)量只觀測(cè)到重力位的鉛垂一次導(dǎo)數(shù)，即Δg或Uz，而重力梯度測(cè)量可以得到重力位的一次導(dǎo)數(shù)Gx，Gy，Gz在x，y，z方向上的變化率，即重力位的二次導(dǎo)數(shù).重力梯度張量由9個(gè)重力梯度分量構(gòu)成：

根據(jù)場(chǎng)論可知，Gxy＝Gyx，Gxz＝Gzx，Gyz＝Gzy，Gxx＋Gyy＋Gzz＝0.因此，重力梯度張量的9個(gè)分量中僅有5個(gè)是獨(dú)立的.

重力梯度異常反映了地下密度突變引起的重力異常的變化，突出場(chǎng)源體的細(xì)節(jié)，具有比重力異常本身高的分辨率，綜合利用各梯度張量分量信息能夠提高地質(zhì)解釋的準(zhǔn)確性.

單個(gè)直立長(zhǎng)方體模型的重力、重力梯度正演計(jì)算解析公式如下［30-32］：

將所有直立長(zhǎng)方體模型異常在地面某點(diǎn)的異常響應(yīng)累加求和（圖1），就可以得出地下半空間在該點(diǎn)的異常響應(yīng)值，可記為：

其中，Δgi為第i個(gè)觀測(cè)點(diǎn)處的重力（重力梯度）異常值，Gij為第i個(gè)觀測(cè)點(diǎn)處第j個(gè)模型單元所對(duì)應(yīng)的核函數(shù)，ρj為第j個(gè)模型單元的密度值.

圖1 地下剖分模型示意圖（a）地下剖分單元示意圖；（b）單個(gè)直立長(zhǎng)方體模型.Fig.1 Interpretation models of the subsurface

寫成矩陣形式：

其中，d為數(shù)據(jù)向量，G為核矩陣，m為模型向量.正演問題就是已知上述G和m，求d，而反演問題則是已知G和d，求m.

后文的基于物性模型的重力（重力梯度）數(shù)據(jù)三維正演，觀測(cè)數(shù)據(jù)是二維的，而地下剖分網(wǎng)格單元是三維的.因此根據(jù)公式（5）計(jì)算觀測(cè)面上的重力（重力梯度）異常需要五重循環(huán).當(dāng)測(cè)點(diǎn)數(shù)據(jù)量大，剖分網(wǎng)格單元個(gè)數(shù)多時(shí)，每次循環(huán)都要計(jì)算單個(gè)G以及G和m的乘積，計(jì)算量相當(dāng)大.而在三維物性反演中，需要反復(fù)迭代進(jìn)行正演計(jì)算，因此需要更多的計(jì)算時(shí)間，以致無(wú)法實(shí)現(xiàn)海量重磁數(shù)據(jù)的快速反演.為縮短計(jì)算時(shí)間，需尋求快速算法.因此本文首先采用GPU并行算法進(jìn)行正演計(jì)算加速試驗(yàn)研究.

3 正演計(jì)算GPU并行策略

3.1 GPU簡(jiǎn)介

新近發(fā)展的NVIDIA GPU（Graphics Processing Unit）通用計(jì)算技術(shù)，在處理能力和存儲(chǔ)器帶寬方面相對(duì)CPU計(jì)算技術(shù)有明顯優(yōu)勢(shì)，在成本和功耗上也不需要付出太大代價(jià)，已在諸多領(lǐng)域得到了廣泛應(yīng) 用.GPU 以 CUDA（Compute Unified Device Architecture）作為新的并行計(jì)算設(shè)備的軟硬件體系，其編程模型將CPU作為主機(jī)（host），負(fù)責(zé)邏輯性強(qiáng)的事務(wù)和串行計(jì)算；GPU作為設(shè)備（device），負(fù)責(zé)執(zhí)行高度并行的線程化任務(wù).運(yùn)行在GPU上的CUDA并 行 計(jì) 算 函 數(shù) 定 義 為 kernelFunction〈〈〈dimGrid，dimBlock〉〉〉（Parameters）.其中，dimGrid為并行計(jì)算所用的線程塊（block）數(shù)目，dimBlock為每個(gè)block中所包含的線程（thread）數(shù)目.同一線程塊內(nèi)在線程執(zhí)行過程中可以通過共享存儲(chǔ)器（Shared memory）共享信息，為了便于并行控制，CUDA提供了四個(gè)內(nèi)建變量（blockIdx，gridDim，threadIdx，blockDim）來(lái)描述線程塊和線程網(wǎng)格的維數(shù)及索引.這些內(nèi)建變量可以有多維形式，我們使用一維線程和一維線程塊，完成正演計(jì)算.

3.2 正演計(jì)算GPU并行策略

地面上某一測(cè)點(diǎn)的重力（重力梯度）正演計(jì)算，需要對(duì)地下半空間所有的長(zhǎng)方體單元異常響應(yīng)值一一計(jì)算并進(jìn)行疊加，即在地下半空間X方向、Y方向和Z方向完成三重模型域循環(huán).對(duì)于面積性觀測(cè)數(shù)據(jù)，還要包括在X、Y方向測(cè)點(diǎn)方向的兩重?cái)?shù)據(jù)域循環(huán).鑒于該過程中計(jì)算每個(gè)點(diǎn)核矩陣（G）的獨(dú)立性，可以將五重循環(huán)同時(shí)進(jìn)行并行化.但當(dāng)模型體數(shù)目變大時(shí)，如100×100×100，模型域循環(huán)數(shù)目就已達(dá)到百萬(wàn)，并且還會(huì)受到顯存大小的限制.因此，這里只將三重模型域循環(huán)在GPU設(shè)備上進(jìn)行并行計(jì)算.一種的最直接方法是按照循環(huán)方式采用線程塊和線程為多維的形式進(jìn)行并行索引，將模型域三重循環(huán)（X、Y、Z 方向循環(huán)）分別利用blockIdx.x和blockIdx.y和threadIdx.x有效索引表示，如圖2所示.然而當(dāng)模型數(shù)據(jù)體很大時(shí)，這種方法會(huì)在分配線程的過程中會(huì)產(chǎn)生很多的線程空閑，從而造成線程資源的巨大浪費(fèi).為了更有效地利用線程塊內(nèi)的每一個(gè)線程，我們將三維模型域循環(huán)利用一維線程塊和線程進(jìn)行有效索引.圖3是三維數(shù)據(jù)體映射為一維數(shù)據(jù)體示意圖，可以根據(jù)所對(duì)應(yīng)的映射關(guān)系，根據(jù)內(nèi)建變量blockDim.x，blockIdx.x和threadIdx.x完成X、Y和Z方向循環(huán)變量的有效索引.

具體有效索引表示如下：

線程索引：index＝blockIdx.x＊blockDim.x＋threadIdx.x；

Z方向循環(huán)索引表示：zindex＝index%nz；

臨時(shí)變量：tmp＝（index-zindex）／nz；

Y方向循環(huán)索引表示：yindex＝tmp%ny；

X方向循環(huán)索引表示：xindex＝（tmp-yindex）／ny.

圖2 二維線程塊與一維線程有效索引表示Fig.2 The index representation for 2Dblocks and 1Dthreads

圖3 三維數(shù)據(jù)映射一維數(shù)據(jù)體示意圖Fig.3 The mapping relations between 3Ddata and 1Ddata

圖4 優(yōu)化的并行歸約流程圖Fig.4 The flow chart of optimized parallel reduction

這樣，每個(gè)線程完成以后，需要對(duì)求得的異常響應(yīng)值進(jìn)行求和運(yùn)算.為了充分利用共享存儲(chǔ)器（shared memory）訪問速度快等優(yōu)勢(shì)，采用經(jīng)過優(yōu)化的并行歸約算法（reduction）.優(yōu)化的并行歸約算法在編程技巧和存儲(chǔ)器使用等方面做了優(yōu)化改進(jìn)（包括避免取模、分支要求，循環(huán)完全展開等），有關(guān)其算法進(jìn)行優(yōu)化的詳細(xì)內(nèi)容請(qǐng)參考 CUDA programming guide［25，28］.如圖4所示，在第一層次（淺綠色），首先將每個(gè)block內(nèi)的線程進(jìn)行優(yōu)化并行歸約到第一個(gè)線程上；在第二層次（橘紅色），將所有的block得到的值傳遞到一個(gè)或幾個(gè)block內(nèi)的線程進(jìn)行優(yōu)化并行歸約，一直循環(huán)計(jì)算得到一個(gè)疊加值，從而完成并行計(jì)算.

4 GPU正演計(jì)算結(jié)果分析

這里采用的數(shù)值試驗(yàn)硬件環(huán)境為：CPU配置為Intel Q9550.主頻2.83GHz，4G內(nèi)存.GPU顯卡型號(hào)為NVIDIA Quadro 2000.192個(gè)處理核心，993MB顯存.軟件環(huán)境為：CUDA驅(qū)動(dòng)版本號(hào)為CUDA 4.0 Windows 764位操作系統(tǒng)，Visual Studio 2008集成開發(fā)編譯環(huán)境.以下計(jì)算時(shí)間為計(jì)算循環(huán)十次取平均所得的統(tǒng)計(jì)時(shí)間.

以下計(jì)算中所用模型皆可視為三維復(fù)雜模型.因?yàn)?，這里所用的觀測(cè)地面下剖分而成的長(zhǎng)方體單元模型物性值是隨機(jī)的.因此，不同長(zhǎng)方體單元隨機(jī)物性值的組合可以代表不同形態(tài)的復(fù)雜模型.對(duì)于特殊的簡(jiǎn)單模型，由于剖分的長(zhǎng)方體單元大部分值為零，在計(jì)算過程中可以忽略.這里只針對(duì)一般情況進(jìn)行分析.

4.1 精度分析

圖5 利用不同數(shù)據(jù)類型在GPU上計(jì)算與CPU上計(jì)算的得到的相對(duì)誤差（a）單精度數(shù)據(jù)類型；（b）雙精度數(shù)據(jù)類型.Fig.5 The relative difference between GPU and CPU results（a）Single－precision data type.（b）double－precision data type.

目前GPU的單精度計(jì)算性能要遠(yuǎn)遠(yuǎn)超過雙精度計(jì)算性能［25］.具體對(duì)于重力異常正演計(jì)算來(lái)說(shuō)，利用單精度數(shù)據(jù)類型進(jìn)行計(jì)算要比雙精度數(shù)據(jù)類型快四倍還要多，但是利用單精度數(shù)據(jù)計(jì)算并不能達(dá)到要求.這是由于重力、重力梯度正演公式中有反正切函數(shù)和對(duì)數(shù)函數(shù)的存在，在利用單精度數(shù)據(jù)計(jì)算過程中舍入誤差會(huì)變得更加明顯，通過誤差傳遞從而影響最終結(jié)果.分別在GPU設(shè)備上用單精度數(shù)據(jù)類型和雙精度數(shù)據(jù)類型進(jìn)行了計(jì)算，并與該程序在CPU上的雙精度數(shù)據(jù)類型計(jì)算結(jié)果進(jìn)行對(duì)比（以重力正演計(jì)算32×32數(shù)據(jù)個(gè)數(shù)為例）.從圖5中可以看出，當(dāng)數(shù)據(jù)使用單精度數(shù)據(jù)類型時(shí)候，最大相對(duì)誤差為1×10-4，并不能滿足計(jì)算精度要求.正如M.Moorkamp所指出，當(dāng)數(shù)據(jù)類型為雙精度時(shí)候最大相對(duì)誤差為1.5×10-12，基本滿足計(jì)算要求［33］.但是，隨著計(jì)算設(shè)備的改進(jìn)（如GTX 500Series）和計(jì)算能力的提高（CUDA 4.0），利用GPU得到的結(jié)果與利用CPU計(jì)算得到的結(jié)果是一致的.因此，在重力正演計(jì)算過程中必須采用雙精度數(shù)據(jù)類型或者更高數(shù)據(jù)類型進(jìn)行計(jì)算，并且對(duì)精度要求已經(jīng)不成問題.

4.2 共享存儲(chǔ)器優(yōu)化

有效地采用共享存儲(chǔ)器隱藏延遲，可以在一定程度上提高加速比.本文在異常值疊加過程中采用經(jīng)過優(yōu)化的并行歸約算法.將地下剖分成不同規(guī)模的長(zhǎng)方體單元，以計(jì)算地面上某一點(diǎn)的重力異常為例，對(duì)采用并行歸約和未采用并行歸約的算法進(jìn)行了時(shí)間對(duì)比（去掉數(shù)據(jù)在CPU與GPU之間傳輸所消耗的時(shí)間影響）.從圖6中可以看出，當(dāng)數(shù)據(jù)規(guī)模小于64×32×32時(shí)，采用并行歸約算法比未采用并行歸約算法的速度要慢，這時(shí)GPU處理性能由于數(shù)據(jù)量小，性能并沒有體現(xiàn).當(dāng)數(shù)據(jù)規(guī)模大于（64×64×32）時(shí)，GPU處理明顯加快，數(shù)據(jù)規(guī)模再進(jìn)一步變大時(shí)，加速比基本上呈線性增長(zhǎng)，其最高可達(dá)15倍左右.

圖6 異常值疊加過程中的優(yōu)化效率對(duì)比Fig.6 Comparsion of reduction and non－reduction for anomaly summation

4.3 不同數(shù)據(jù)規(guī)模加速比

對(duì)比分析了32×32個(gè)測(cè)點(diǎn)、不同模型規(guī)模重力、重力梯度數(shù)據(jù)正演計(jì)算的加速比.從圖7可以看出，和上述共享存儲(chǔ)器優(yōu)化中的加速比形態(tài)類似，當(dāng)數(shù)據(jù)規(guī)模小于8×8×1時(shí)，采用GPU并行算法比單線程CPU算法的速度要慢，這時(shí)GPU處理性能由于數(shù)據(jù)量小，性能并沒有體現(xiàn)出來(lái).由于重力和重力梯度正演在計(jì)算過程中只有很少的系統(tǒng)開銷，當(dāng)數(shù)據(jù)規(guī)模大于8×8×1時(shí)，GPU計(jì)算速度明顯加快，一直到32×32×4.由于計(jì)算能力達(dá)到飽和，加速比大體上呈線性增長(zhǎng)，最高可達(dá)60倍.對(duì)比重力、重力梯度數(shù)據(jù)，由于并行計(jì)算中的大部分時(shí)間用在計(jì)算公式（4）中反正切函數(shù)和自然對(duì)數(shù)函數(shù)上，因此正演公式的復(fù)雜度決定著計(jì)算效率的快慢.由于重力正演公式比重力梯度正演公式復(fù)雜度高，因此前者的加速效果更明顯.如圖8，對(duì)于Uxx和Uxy正演計(jì)算加速比可達(dá)50倍左右.

圖7 不同規(guī)模模型重力正演GPU與單線程CPU計(jì)算時(shí)間對(duì)比（測(cè)點(diǎn)數(shù)為32×32）Fig.7 Comparison of GPU and CPU computing time for gravity modeling（32×32observational points）

圖8 不同規(guī)模模型重力梯度（Uxx和Uxy）正演GPU與單線程CPU計(jì)算時(shí)間對(duì)比（測(cè)點(diǎn)數(shù)為32×2）Fig.8 Comparison of GPU and CPU computing time for gravimetry modeling（32×32observational points）

5 重力（重力梯度）數(shù)據(jù)GPU并行三維反演策略

基于上述GPU正演快速算法的有效實(shí)現(xiàn)，本節(jié)首先討論GPU相關(guān)成像三維物性反演思路，然后探討基于GPU的重力（重力梯度）數(shù)據(jù)三維約束反演策略.

5.1 GPU相關(guān)成像三維物性反演策略

根據(jù)測(cè)區(qū)范圍將地下半空間剖分為大小相等、規(guī)則排列的長(zhǎng)方體網(wǎng)格單元.首先，根據(jù)三維相關(guān)成像原理［5］，利用實(shí)測(cè)重力（重力梯度）異常計(jì)算地下網(wǎng)格單元節(jié)點(diǎn)q的相關(guān)系數(shù)Cq，該相關(guān)系數(shù)在-1與1之間變化.其次，根據(jù)實(shí)測(cè)區(qū)域的密度范圍，將相關(guān)系數(shù)線性加權(quán)轉(zhuǎn)化后，作為該網(wǎng)格節(jié)點(diǎn)相鄰的網(wǎng)格單元的初始密度值，利用直立長(zhǎng)方體無(wú)奇點(diǎn)正演公式進(jìn)行重力（重力梯度）正演計(jì)算.最后，利用實(shí)測(cè)重力（重力梯度）異常與正演計(jì)算結(jié)果的殘差進(jìn)行相關(guān)成像，將殘差的相關(guān)成像的線性加權(quán)結(jié)果作為模型的修改量，再對(duì)模型進(jìn)行正演計(jì)算.如此迭代反復(fù)，直到殘差達(dá)到合理要求為止.通過分析不難發(fā)現(xiàn)，對(duì)于三維相關(guān)成像求取某一網(wǎng)格節(jié)點(diǎn)的相關(guān)系數(shù)值，不受其他網(wǎng)格節(jié)點(diǎn)的影響.對(duì)于正演計(jì)算求取某一觀測(cè)點(diǎn)的異常值，同樣不受其它觀測(cè)點(diǎn)值的影響.可以采用GPU并行計(jì)算技術(shù)提高計(jì)算效率.將三維相關(guān)成像和正演計(jì)算算法改造成GPU并行程序，并將相關(guān)數(shù)據(jù)及參數(shù)傳遞到GPU上，使得三維相關(guān)成像和正演計(jì)算在GPU進(jìn)行.最后將計(jì)算結(jié)果傳回主機(jī)端內(nèi)存中［34］.

5.2 重力（重力梯度）數(shù)據(jù)三維約束物性反演GPU加速策略

重力（重力梯度）異常三維約束物性反演是基于反演理論、在最小二乘意義下使目標(biāo)函數(shù)達(dá)到極小的線性或非線性反演.在引入足夠的、適當(dāng)?shù)募s束條件下，使反演目標(biāo)函數(shù)最小化，能夠得到接近于實(shí)際地質(zhì)情況的物性分布和幾何形態(tài).即：

其中，φd＝ ‖Wd（Gm －d）‖2為數(shù)據(jù)擬合函數(shù)，φm＝‖Wm（m－mref）‖為模型目標(biāo)函數(shù).Wd為數(shù)據(jù)空間加權(quán)矩陣，Wm為模型空間加權(quán)矩陣.μ為拉格朗日算子.

對(duì)（6）式，通過采用不同的加權(quán)因子Wd和Wm以及不同的計(jì)算策略，可以得到不同種意義上的解.對(duì)于式（6），求解方式有很多種，首先為模型空間迭代法，給定初始Wd和Wm以及m0采用公式：

進(jìn)行迭代求解.上述兩種直接求解法，需要構(gòu)造顯式矩陣，并且涉及到矩陣求逆.只是在最初的數(shù)據(jù)量很小的二維物性反演中適用.當(dāng)數(shù)據(jù)量很大時(shí)，需要采取相應(yīng)的計(jì)算技巧才能實(shí)現(xiàn).Li和Oldenburg結(jié)合小波壓縮算法存儲(chǔ)核矩陣，對(duì)目標(biāo)函數(shù)采用對(duì)數(shù)障礙法求解［11］，實(shí)現(xiàn)了大數(shù)據(jù)重磁數(shù)據(jù)的三維反演計(jì)算.Pilkington［18］提出了利用預(yù)條件共軛梯度法進(jìn)行三維磁化率成像.共軛梯度算法沿著由已知點(diǎn)處梯度構(gòu)造的共軛方向搜索目標(biāo)函數(shù)的極小值，理論上經(jīng)過有限步迭代，算法就可以收斂.另外，非線性的廣義隨機(jī)算法模擬退火算法等，可使反演求解過程穩(wěn)定，約束條件容易結(jié)合，但計(jì)算速度和維數(shù)困難同樣制約其發(fā)揮作用，采取針對(duì)性措施（等效存儲(chǔ)幾何格架）后，使三維反演進(jìn)入實(shí)用化階段.采用預(yù)條件共軛梯度法對(duì)目標(biāo)函數(shù)式（6）進(jìn)行求解的思路如下：

加入約束因子后，目標(biāo)函數(shù)（6）可寫為

圖9 基于相關(guān)成像算法的三維物性反演流程圖Fig.9 Flow chart of 3Dproperty inversion based on correlation imaging algorithm

圖10 GPU相關(guān)成像三維物性反演算法流程圖Fig.10 Flow chart of 3Dproperty inversion based on correlation imaging algorithm using GPU

原問題可變成求解：

將重力（重力梯度）異常三維約束算法任務(wù)劃分，確定程序中的并行部分.對(duì)于公式（7）和（8）涉及到矩陣與向量計(jì)算，矩陣求逆，矩陣與矩陣計(jì)算，在GPU加速計(jì)算中，已經(jīng)能夠很好地實(shí)現(xiàn)，并能取得很高的加速比，關(guān)鍵在于如何合理地利用好計(jì)算機(jī)內(nèi)存和顯卡內(nèi)存.當(dāng)反演數(shù)據(jù)以及反演的模型空間規(guī)模變得很大時(shí)，雖然預(yù)條件共軛梯度法中的矢量運(yùn)算在一定程度上減輕了內(nèi)存使用量，但仍然耗費(fèi)很長(zhǎng)的時(shí)間.NVIDIA的CUBLAS是在CUDA基礎(chǔ)上實(shí)現(xiàn)的一種并行的BLAS，并且實(shí)現(xiàn)了對(duì)CUDA過程的封裝.這樣很容易將預(yù)條件共軛梯度移植到GPU上進(jìn)行計(jì)算.相關(guān)的共軛梯度GPU加速的計(jì)算案例已證明了該方法的有效性.

6 結(jié)論與討論

利用GPU能夠很好地并行實(shí)現(xiàn)重力、重力梯度正演快速計(jì)算，并且可以獲得很高的加速比.數(shù)值試驗(yàn)顯示，利用雙精度數(shù)據(jù)類型或者更高精度數(shù)據(jù)類型進(jìn)行計(jì)算，可達(dá)到精度要求.文章還討論了GPU并行計(jì)算在兩種反演方法中的策略，為快速三維反演技術(shù)提供了借鑒.本文取得的認(rèn)識(shí)如下：

（1）由于GPU單元在傳輸上具有很高的帶寬，所以在數(shù)據(jù)傳輸方面不會(huì)占有很長(zhǎng)的時(shí)間.因此，真正能夠?qū)μ岣咚俣扔行У氖敲總€(gè)線程所處理的復(fù)雜度，復(fù)雜度越大，加速比就越高.由于重力正演公式比重力梯度正演公式復(fù)雜度高，因此前者的加速效果更明顯.

（2）利用共享存儲(chǔ)器的存儲(chǔ)特點(diǎn)以及線程設(shè)計(jì)，對(duì)每個(gè)線程計(jì)算出的響應(yīng)值疊加采用經(jīng)過優(yōu)化過的并行歸約算法，在相同條件下能夠得到高達(dá)15倍的加速比.

（3）由于GPU在一定時(shí)間內(nèi)還會(huì)受到顯存的限制，因此在處理更大規(guī)模數(shù)據(jù)時(shí)，可以采用分塊的方式進(jìn)行.

（4）在以往的基于物性模型的三維重力、重力梯度反演過程中，由于核矩陣（G）計(jì)算時(shí)間較長(zhǎng)，須采用一次性計(jì)算完成并存儲(chǔ)核矩陣（G）.此加速算法策略可以無(wú)須一次性存儲(chǔ)G，這為基于物性模型的三維重力、重力梯度快速反演奠定了基礎(chǔ).

（5）本文GPU加速計(jì)算是針對(duì)重力、重力梯度數(shù)據(jù)開展試驗(yàn)的.值得說(shuō)明的是，此方法可以推廣到其它位場(chǎng)數(shù)據(jù)的類似計(jì)算.

（References）

［1］Nabighain M N，Grauch V J S，Hansen R O，et al.75th anniversary－The historical development of the magnetic method in exploration.Geophysics，2005，70（6）：33ND－61ND.

［2］Nabighain M N，Ander M E，Grauch V J S，et al.75th anniversary－Historical development of the gravity method in exploration.Geophysics，2005，70（6）：63ND－89ND.

［3］GuspíF，Introcaso B.A sparse spectrum technique for gridding and separating potential field anomalies.Geophysics，2000，65（4）：1154－1161.

［4］Meng X H，Guo L H，Chen Z X，et al.A method for gravity anomaly separation based on preferential continuation and its application.Applied Geophysics，2009，6（3）：217－225.

［5］郭良輝，孟小紅，石磊等.重力和重力梯度數(shù)據(jù)三維相關(guān)成像.地球物理學(xué)報(bào)，2009，52（4）：1098－1106.Guo L H，Meng X H，Shi L，et al.3Dcorrelation imaging for gravity and gravity gradiometry data.Chinese J.Geophys.（in Chinese），2009，52（4）：1098－1106.

［6］郭良輝，孟小紅，石磊.磁異常ΔT三維相關(guān)成像.地球物理學(xué)報(bào)，2010，53（2）：435－441.Guo L H，Meng X H，Shi L.3Dcorrelation imaging for magnetic anomalyΔTdata.Chinese J.Geophys.（in Chinese），2010，53（2）：435－441.

［7］劉天佑，楊宇山，李媛媛等.大型積分方程降階解法與重力資料曲面延拓.地球物理學(xué)報(bào)，2007，50（1）：290－296.Liu T Y，Yang Y S，Li Y Y，et al.The order－depression solution for large－scale integral equation and its application in the reduction of gravity data to a horizontal plane.Chinese J.Geophys.（in Chinese），2007，50（1）：290－296.

［8］劉天佑.位場(chǎng)勘探數(shù)據(jù)處理新方法.北京：科學(xué)出版社，2007.Liu T Y.New Data Processing Methods for Potential Field Exploration（in Chinese）.Beijing：Science Press，2007.

［9］Fedu M，Quarta T.Wavelet analysis for the regional－residual and local separation of potential field anomalies.Geophysical Prospecting，1998，46（5）：507－525.

［10］Last B J，Kubik K.Compact gravity inversion.Geophysics，1983，48（6）：713－721.

［11］Li Y G，Oldenburg D W.Fast inversion of large－scale magnetic data using wavelet transforms and a logarithmic barrier method.Geophys.J.Int.，2003，152（2）：251－265.

［12］Constable S C，Parker R L，Constable C G.Occam′s inversion：A practical algorithm for generation smooth models from electromagnetic sounding data.Geophysics，1987，52（3）：289－300.

［13］Li Y G，Oldenburg D W.3Dinversion of gravity data.Geophysics，1998，63（1）：109－119.

［14］Li Y G，Oldenburg D W.3Dinversion of magnetic data.Geophysics，1996，61（2）：394－408.

［15］Li Y G，Oldenburg D W.Joint inversion of surface and threecomponent borehole magnetic data.Geophysics，1996，65（2）：540－552.

［16］Portniaguine O， Zhdanov M S. Focusing geophysical inversion images.Geophysics，1999，64（3）：874－887.

［17］Portniaguine O，Zhdanov M S.3－D magnetic inversion with data compression and image focusing.Geophysics，2002，67（5）：1532－1541.

［18］Pilkington M. 3－D magnetic imaging using conjugate gradients.Geophysics，1997，62（4）：1132－1142.

［19］Nagihara S，Hall S A.Three－dimensional gravity inversion using simulated annealing：Constraints on the diapiric roots of allochthonous salt structures.Geophysics，2001，66（5）：1438－1449.

［20］Boulanger O，Chouteau M.Constraints in 3Dgravity inversion.Geophysical Prospecting，2001，49（2）：265－280.

［21］Fedi M，Rapolla A.3Dinversion of gravity and magnetic data with depth resolution.Geophysics，1999，64（2）：452－460.

［22］姚長(zhǎng)利，郝天珧，管志寧等.重磁遺傳算法三維反演中高速計(jì)算及有效存儲(chǔ)方法技術(shù).地球物理學(xué)報(bào)，2003，46（2）：252－258.Yao C L， Hao T Y，Guan Z N，et al. High－speed computation and efficient storage in 3－D gravity and magnetic inversion based on genetic algorithms.Chinese J.Geophys.（in Chinese），2003，46（2）：252－258.

［23］姚長(zhǎng)利，鄭元滿，張聿文.重磁異常三維物性反演隨機(jī)子域法方法技術(shù).地球物理學(xué)報(bào)，2007，50（5）：1576－1583.Yao C L，Zheng Y M，Zhang Y W.32Dgravity and magnetic inversion for physical properties using stochastic subspaces.Chinese J.Geophys.（in Chinese），2007，50（5）：1576－1583.

［24］Li Y G，Oldenburg D W.Fast inversion of large－scale magnetic data using wavelet transform and a logarithmic barrier method.Geophys.J.Int.，2003，152（2）：251－265.

［25］Zhang S，Chu Y L.GPU High－Performance Computing-CUDA.Beijing：China Water Power Press，2009.

［26］劉國(guó)峰，劉洪，李博等.山地地震資料疊前時(shí)間偏移方法及其 GPU實(shí)現(xiàn).地球物理學(xué)報(bào)，2009，52（12）：3101－3108.Liu G F，Liu H，Li B，et al.Method of prestack time migration of seismic data of mountainous regions and its GPU implementation.Chinese J.Geophys.（in Chinese），2009，52（12）：3101－3108.

［27］李博，劉國(guó)峰，劉洪.地震疊前時(shí)間偏移的一種圖形處理器提速實(shí)現(xiàn)方法.地球物理學(xué)報(bào)，2009，52（1）：245－252 Li B，Liu G F，Liu H.A method of using GPU to accelerate seismic pre－stack time migration.Chinese J.Geophys.（in Chinese），2009，52（1）：245－252.

［28］Sanders J，Kandrot E.CUDA by Example－An Introduction to General Purpose GPU Programming.Pearson Education Press，2010.

［29］劉國(guó)峰，劉欽，李博等.油氣勘探地震資料處理GPU／CPU協(xié)同并行計(jì)算.地球物理學(xué)進(jìn)展，2009，24（5）：1671－1678.Liu G F，Liu Q，Li B，et al.GPU／CPU co－processing parallel computation for seismic data processing in oil and gas exploration.Progress in Geophys （in Chinese），2009，24（5）：1671－1678.

［30］Nagy D.The gravitational attraction of a right rectangular prism.Geophysics，1966，31（2）：361－371.

［31］Li X，Chouteau M.Three－dimensional gravity modeling in all space.Surveys in Geophysics，1998，19（4）：339－368.

［32］郭志宏，管志寧，熊盛青.長(zhǎng)方體ΔΤ場(chǎng)及其梯度場(chǎng)無(wú)解析奇點(diǎn)理論表達(dá)式.地球物理學(xué)報(bào)，2004，47（6）：1131－1138.Guo Z H，Guan Z N，Xiong S Q.CuboidΔΤand its gradient forward theoretical expressions without analytic odd points.Chinese J.Geophys.（in Chinese），2004，47（6）：1131－1138.

［33］Moorkamp M，Jegen M，Roberts A，et al.Massively parallel forward modeling of scalar and tensor gravimetry data.Computers and Geosciences，2010，36（5）：680－686.

［34］Chen Z X，Meng X H，Guo L H，et al.GICUDA：A parallel program for 3Dcorrelation imaging of large scale gravity and gravity gradiometry data on graphics processing units with CUDA.Computers and Geosciences，2012，46：119－128.