陳召曦，孟小紅，郭良輝，劉國峰

1 地質(zhì)過程與礦產(chǎn)資源國家重點實驗室，中國地質(zhì)大學，北京 100083

2 地下信息探測技術(shù)與儀器教育部重點實驗室，中國地質(zhì)大學，北京 100083

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

1 引 言

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

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

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

2 正演理論

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

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

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

單個直立長方體模型的重力、重力梯度正演計算解析公式如下［30-32］：

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

其中，Δgi為第i個觀測點處的重力（重力梯度）異常值，Gij為第i個觀測點處第j個模型單元所對應的核函數(shù)，ρj為第j個模型單元的密度值.

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

寫成矩陣形式：

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

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

3 正演計算GPU并行策略

3.1 GPU簡介

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

3.2 正演計算GPU并行策略

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

臨時變量：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

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

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

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

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

4.1 精度分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

加入約束因子后，目標函數(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

原問題可變成求解：

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

6 結(jié)論與討論

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

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

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

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

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

（5）本文GPU加速計算是針對重力、重力梯度數(shù)據(jù)開展試驗的.值得說明的是，此方法可以推廣到其它位場數(shù)據(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)成像.地球物理學報，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)成像.地球物理學報，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］劉天佑，楊宇山，李媛媛等.大型積分方程降階解法與重力資料曲面延拓.地球物理學報，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］劉天佑.位場勘探數(shù)據(jù)處理新方法.北京：科學出版社，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］姚長利，郝天珧，管志寧等.重磁遺傳算法三維反演中高速計算及有效存儲方法技術(shù).地球物理學報，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］姚長利，鄭元滿，張聿文.重磁異常三維物性反演隨機子域法方法技術(shù).地球物理學報，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］劉國峰，劉洪，李博等.山地地震資料疊前時間偏移方法及其 GPU實現(xiàn).地球物理學報，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］李博，劉國峰，劉洪.地震疊前時間偏移的一種圖形處理器提速實現(xiàn)方法.地球物理學報，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］劉國峰，劉欽，李博等.油氣勘探地震資料處理GPU／CPU協(xié)同并行計算.地球物理學進展，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］郭志宏，管志寧，熊盛青.長方體ΔΤ場及其梯度場無解析奇點理論表達式.地球物理學報，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.