王泰馬國慶，李 野，林 松

1.吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026 2.吉林省電力勘測設(shè)計(jì)院，長春 130021

0 引言

多年來，重力勘探被廣泛地應(yīng)用于大地構(gòu)造、勘探礦產(chǎn)和油氣源，以及工程與環(huán)境調(diào)查等研究領(lǐng)域[1-3]。在重力勘探研究中，三維重力反演已經(jīng)成為一種重要和有效的處理和解釋手段。通常情況下，重力數(shù)據(jù)的三維反演往往是病態(tài)的，即需要一些先驗(yàn)信息或約束條件來保證解的唯一性和準(zhǔn)確性?；谡齽t化方法，地球物理學(xué)家們做了很多工作來解決反演的不穩(wěn)定性和非唯一性問題。Li等[4-5]采用目標(biāo)函數(shù)的深度權(quán)函數(shù)進(jìn)行磁化率反演和密度反演, 通過賦予隨深度增加的棱柱體單元更多的權(quán)重來抵消核函數(shù)的固有衰減(即“趨膚效應(yīng)”)。隨著重力數(shù)據(jù)采集量和精度的提高[6-8]，快速處理大規(guī)模高精度數(shù)據(jù)在處理與解釋方面越來越常見。解決大規(guī)模三維重力數(shù)據(jù)反演計(jì)算問題面臨兩個(gè)困難：需要足夠的計(jì)算機(jī)內(nèi)存來儲存反演過程中生成的矩陣；減少反演求解過程中矩陣向量乘法的計(jì)算用時(shí)。對于這樣的問題，特別是在航空測量成千上萬個(gè)數(shù)據(jù)點(diǎn)的情況下，尋求一種快速、有效的算法對于快速解決三維反演問題十分重要。

共軛梯度算法(conjugate gradient algorithm, CG)是一種解決大型線性重力反演的通用迭代算法。對其使用有效的預(yù)處理因子，預(yù)處理共軛梯度算法(PCG)已經(jīng)在很多領(lǐng)域研究[9-12]中證明了其良好的性能和應(yīng)用效果。然而，不同的預(yù)處理因子在提高反演速率和減少迭代次數(shù)方面的能力有所不同。在迭代過程中，準(zhǔn)備(分解)預(yù)處理因子也需要額外的計(jì)算時(shí)間。因此，反演計(jì)算效率由迭代次數(shù)和總計(jì)算時(shí)間綜合決定。隨著圖形處理單元(GPU)的計(jì)算能力和可編譯能力的不斷增強(qiáng)，利用GPU輔助進(jìn)行科學(xué)計(jì)算逐漸成為研究的熱點(diǎn)。CUDA[13](統(tǒng)一計(jì)算設(shè)備架構(gòu))是顯卡廠商N(yùn)VIDIA推出的硬件和軟件架構(gòu)平臺，使GPU成為一種并行計(jì)算工具，解決復(fù)雜的計(jì)算問題，而不只是完成當(dāng)初的圖像渲染。如今，CUDA已經(jīng)成為最適合并行計(jì)算的手段之一：Moorkamp[14]等使用CUDA進(jìn)行了標(biāo)量和張量數(shù)據(jù)的正演并行計(jì)算；Chen等[15]利用CUDA對重力及梯度數(shù)據(jù)進(jìn)行了三維相關(guān)成像的并行計(jì)算；Liu等[16]提出了CUDA下基于概率成像的三維磁化率反演。

本文為提高三維重力數(shù)據(jù)反演效率，基于CPU和GPU給出一種并行的預(yù)處理共軛梯度算法。首先簡單地介紹三維重力數(shù)據(jù)正反演的原理，給出本文的并行算法的流程；進(jìn)而通過不同網(wǎng)格大小的單一模型試驗(yàn)比較出CPU/GPU算法相對于串行CPU算法的加速能力，以及用組合模型試驗(yàn)比較基于兩種不同預(yù)處理算子的并行預(yù)處理共軛梯度算法與傳統(tǒng)共軛梯度算法在重力反演中的效果和加速能力；最后，將該算法應(yīng)用于美國路易斯安那州Vinton鹽丘實(shí)測重力數(shù)據(jù)三維反演中，進(jìn)一步證明本文并行預(yù)處理快速方法的高效性和適用性。

1 方法原理

1.1 重力數(shù)據(jù)正反演理論

根據(jù)牛頓萬有引力定律，三維重力數(shù)據(jù)的正演是計(jì)算由地下密度分布為ρ引起的在地表觀測的重力異常g。為簡單起見，將三維地下空間劃分為有限多的密度一定的長方體單元來計(jì)算重力數(shù)據(jù)。在笛卡爾坐標(biāo)系中，我們定義x軸和y軸為水平方向，z軸為垂直向下，其分別代表東向、北向和深度方向。根據(jù)Haz[17]提出的積分解，由密度為ρ的單一棱柱體單元C在觀測點(diǎn)O=(x0,y0,z0)處產(chǎn)生的重力異常為

(1)

(2)

考慮到有n個(gè)觀測點(diǎn)和m個(gè)長方體模型單元，其線性關(guān)系也可以用矩陣形式表示為

gn×1=Gn×mρm×1。

(3)

式中：g是地表觀測重力異常值；ρ是地下空間密度分布；G是重力正演核函數(shù)。

在三維重力反演中，地下網(wǎng)格個(gè)數(shù)遠(yuǎn)大于地面觀測點(diǎn)的個(gè)數(shù)，也就是說，反演問題是欠定的，需要引入正則化約束來解決反演過程的非唯一性和不確定性[18]，恢復(fù)出大多數(shù)合理的解。我們用正則化約束Tikhnov參數(shù)方程P(ρ)為

Pα(ρ)=φ(ρ)+αs(ρ)=

(4)

式中：Wd是數(shù)據(jù)誤差加權(quán)矩陣；Wz是深度加權(quán)矩陣形式，計(jì)算公式為Wz=(z+z0)-1/2，用于抵消核函數(shù)矩陣G隨深度的固有衰減，克服趨膚效應(yīng)，z是長方體單元平均埋深，z0取決于模型離散化的塊體大小和重力數(shù)據(jù)觀測高度；s(ρ)是穩(wěn)定函數(shù)；φ(ρ)是觀測數(shù)據(jù)與預(yù)測數(shù)據(jù)的擬合函數(shù)，為正則化參數(shù)。φ(ρ)值對反演結(jié)果有很大影響，本文選取一種自適應(yīng)方法選擇正則化參數(shù)，即隨著反演結(jié)果的收斂逐漸減小穩(wěn)定函數(shù)的貢獻(xiàn)值。初次迭代計(jì)算得到的正則化參數(shù)表達(dá)式為

(5)

隨后的迭代計(jì)算中，通過

(6)

計(jì)算正則化參數(shù)的值。

根據(jù)Pilkington[10]研究，式(4)又可以被轉(zhuǎn)換成矩陣形式：

?AM=b，

(7)

ATAM=ATb。

(8)

1.2 共軛梯度算法和預(yù)處理策略

共軛梯度算法是一種求解線性方程組的快速迭代算法，其收斂速度比最速下降法和牛頓法下降法要快。CG算法只進(jìn)行矩陣向量乘積和向量內(nèi)積運(yùn)算，需要低存儲和低計(jì)算成本，被廣泛地應(yīng)用于地球物理及其他領(lǐng)域反演中。重力反演的基本CG算法為

fork=0;m0=0 andr0=AT(b-Am0)

whiler0≠0

k=k+1,

ifk=1,

p1=r0

pk=rk-1+βkpk-1,

end

qk=Apk,

mk=mk-1+akpk,

rk=rk-1-akATqk,

end

(9)

式中；k為迭代次數(shù)，rk為目標(biāo)函數(shù)的梯度；pk和αk分別表示迭代的搜索方向和步長；βk和qk是迭代計(jì)算的中間變量。在算法中可以發(fā)現(xiàn)，系數(shù)矩陣A的逆矩陣不需要被計(jì)算。這對于大規(guī)模數(shù)據(jù)可以節(jié)省很長時(shí)間和很大存儲空間。隨著觀測數(shù)據(jù)量增多或者離散網(wǎng)格數(shù)加大的三維實(shí)際反演問題，矩陣A也變得相當(dāng)大。因此，反演的迭代次數(shù)和計(jì)算時(shí)間應(yīng)該被考慮到。眾所周知，CG算法收斂速度的快慢受系數(shù)矩陣特征值分布控制。為了減少迭代次數(shù)進(jìn)而加快收斂速度，國內(nèi)外許多研究人員將式(8)變換為式(10)，使特征值分布盡量集中在對角線上。這種改善特征值分布的方法叫做預(yù)處理方法，應(yīng)用于CG算法可得到預(yù)處理共軛梯度法(PCG),算法如式(11)。

PATAM=PATb。

(10)

fork=0;m0=0 andr0=AT(b-Am0)

whiler0≠0

zk=Pr

k=k+1,

ifk=1,

p1=z0

pk=zk-1+βkpk-1,

end

qk=Apk,

mk=mk-1+akpk,

rk=rk-1-akATqk,

end

(11)

其中，P是為提高收斂速度的預(yù)處理算子。

可以看出，相比于CG算法(9)，PCG算法(11)只是多了一步zk=Pr，若P為單位矩陣，則PCG算法退化為CG。理想情況下，P是(ATA)-1的近似，這樣PATA近似等于單位矩陣。為了加速CG算法的收斂速度，PATA特征值分布必須比原對稱正定矩陣ATA更加集中，改善了方程的條件數(shù)，會提高迭代速度。因此，如何選取一個(gè)有效的預(yù)處理因子對于減少反演迭代次數(shù)是非常重要的。PCG方法在很多研究中已經(jīng)證明其效率和良好的性能。本文選取兩種不同的預(yù)處理算子應(yīng)用于三維重力反演。

一種最常見和簡單的方法是對角預(yù)處理算子(DP)。這種方法是使算子P等于重力反演中第一次迭代的最小二乘解，在對角預(yù)處理共軛梯度算法中，選取P等于ATA的對角元素，即P=diag(diag(ATA))。

另一種預(yù)處理算子可以從系數(shù)矩陣直接近似推導(dǎo)，為對稱逐次超松弛(SSOR)預(yù)處理方法。SSOR預(yù)處理方法相比于對角預(yù)處理方法包含了更多的系數(shù)矩陣(地球物理信息)的信息。其分解方式是將矩陣ATA分解為

ATA=D-LAA-LAAT，

(12)

其中；D是系數(shù)矩陣；ATA對角元素組成的對角矩陣；LAA是ATA的下三角矩陣。w(0

(13)

雖然預(yù)處理共軛梯度法能夠從減少迭代次數(shù)的方面提高反演收斂速度，但是每一次迭代過程中構(gòu)建(分解)預(yù)處理算子需要花費(fèi)很多時(shí)間，尤其是處理大數(shù)據(jù)量的三維反演。筆者嘗試通過改進(jìn)并行算法來進(jìn)行三維快速反演。

1.3 并行算法的設(shè)計(jì)

近些年來，可編譯的NVIDIA GPU已經(jīng)被開發(fā)用于高性能計(jì)算，這吸引了研究人員和開發(fā)者更多的關(guān)注。當(dāng)處理大量地球物理數(shù)據(jù)時(shí)，擁有流多處理器集群的GPU具有極高的計(jì)算速度。在解決耗時(shí)問題的需求下，GPU已經(jīng)成功應(yīng)用于許多地球物理并行計(jì)算研究中。統(tǒng)一計(jì)算設(shè)備架構(gòu)(CUDA)是一種使用NVIDIA GPU作為數(shù)據(jù)并行計(jì)算設(shè)備的硬件和軟件架構(gòu)。CUDA技術(shù)最初只支持C++語言，需要高水平的編程。隨著GPU計(jì)算能力的提高，MATLAB在版本R2010b[19]后，提供并行計(jì)算工具箱(PCT)使用戶進(jìn)行快速并行計(jì)算更加方便。MATLAB中強(qiáng)大的GPU計(jì)算能力是在Tesla和Quadro GPU計(jì)算產(chǎn)品上開發(fā)的，需要使用最新的支持CUDA的NVIDIA顯卡，例如NVIDIA Tesla 10系列或20系列產(chǎn)品，支持計(jì)算能力在1.3或更高的計(jì)算能力。本文算法的設(shè)計(jì)平臺是MATLAB，用來討論CPU和GPU算法加速問題。

將上節(jié)提到的共軛梯度算法用串行程序分步計(jì)時(shí)后，發(fā)現(xiàn)每一步迭代過程中的更新矩陣r0、qk、rk是最耗時(shí)的部分。在CG算法的重力反演中主要耗時(shí)問題被歸納為矩陣向量乘積運(yùn)算和大系數(shù)矩陣A的內(nèi)積運(yùn)算。因此，我們將這部分用GPU計(jì)算，其他部分仍用CPU作常規(guī)的串行計(jì)算。為了體現(xiàn)GPU的加速能力和適用性，筆者先利用不同網(wǎng)格大小的單一模型進(jìn)行三維重力共軛梯度反演，通過記錄1 000次迭代的時(shí)間來比較CPU串行與CPU/GPU并行算法計(jì)算速度，比較結(jié)果見圖1。顯然CPU/GPU算法計(jì)算能提高計(jì)算效率，尤其是針對大數(shù)據(jù)量的重力數(shù)據(jù)反演。通過曲線可以看到，在沒有預(yù)處理算子的CG反演中，并行算法相比串行程序加快了大約4.5倍。本文程序運(yùn)行的硬件配置如下：

圖1 三維重力共軛梯度反演的CPU和CPU/GPU計(jì)算時(shí)間對比圖Fig.1 CPU and CPU/GPU computing time for the CG algorithm in 3D gravity inversion

CPU Intel Xeon E5-2620 6-core 2.0 GHz, 7.2 GT/s

GPU NVIDIA Tesla C2050

上節(jié)我們提到應(yīng)用不同預(yù)處理算子在預(yù)處理共軛梯度反演中可以減少迭代次數(shù)，在某種程度上使反演收斂快，減少耗時(shí)。然而，隨著數(shù)據(jù)量和網(wǎng)格數(shù)增大，構(gòu)建預(yù)處理算子需要更多額外的時(shí)間應(yīng)該被考慮。因此，本文給出GPU加速的PCG并行算法進(jìn)行三維重力數(shù)據(jù)反演，在減少迭代次數(shù)的同時(shí)，進(jìn)一步解決耗時(shí)問題。并行的PCG算法流程(圖2)為：1)初始化MATLAB環(huán)境；2)讀入觀測重力數(shù)據(jù)，開始計(jì)時(shí)；3)并行計(jì)算模型權(quán)Wm和深度權(quán)Wz，后用gather()函數(shù)將結(jié)果返回至CPU；4)初始化m0，計(jì)算A,b；5)用gpuArray()函數(shù)將以上數(shù)據(jù)傳輸?shù)紾PU，開始迭代，預(yù)處理因子P、搜索方向pk和搜索步長αk用GPU計(jì)算，其余都在CPU串行計(jì)算；6)達(dá)到迭代終止條件，停止計(jì)時(shí)，返回密度結(jié)果m到CPU，繪圖。其中g(shù)puArray()和gather()函數(shù)是由Parallel Computing Toolbox(PCT)工具箱提供的負(fù)責(zé)CPU和GPU之間數(shù)據(jù)收集和傳輸?shù)暮瘮?shù)。

圖2 重力反演的并行PCG算法流程圖Fig.2 Flowchart of the parallel PCG algorithm in gravity inversion

2 模型實(shí)驗(yàn)

為了比較并行預(yù)處理共軛梯度算法在三維重力數(shù)據(jù)反演的加速能力，我們首先模擬一個(gè)由5個(gè)長方體組成的復(fù)雜模型。各個(gè)模型參數(shù)如表1給出。三維地下空間被離散成20×20×10個(gè)小棱柱體單元，其單元大小為100 m×100 m×100 m(x×y×z)。圖3a展示了不同大小和密度的模型體的空間分布。對模型數(shù)據(jù)加入最大值5%的高斯噪聲，模擬真實(shí)數(shù)據(jù)，以驗(yàn)證算法的抗噪性，這樣地面觀測的含噪剩余重力異常如圖3b所示。

表1 組合模型參數(shù)表

我們將前面介紹的兩種不同的預(yù)處理方法應(yīng)用于模型數(shù)據(jù)的三維重力反演，比較其并行程序與CPU串行的CG算法的計(jì)算時(shí)間和結(jié)果，來驗(yàn)證并行預(yù)處理共軛梯度反演的高效性和有效性。3種方法重力反演結(jié)果如圖4。圖4a、b、c是反演結(jié)果在深度為400 m的水平切面圖，圖4d、e、f是在水平x方向?yàn)? 600 m的垂直切面圖。黑色矩形框?yàn)槊總€(gè)模型的真實(shí)位置?？梢园l(fā)現(xiàn)：對角預(yù)處理共軛梯度算法(DPCG)的反演結(jié)果與傳統(tǒng)的CG算法較為接近；而對稱逐次超松弛預(yù)處理共軛梯度算法(SSOR-PCG)反演出的密度值與真實(shí)密度值更為接近，圖像更為緊致，空間上更接近真實(shí)模型位置。

在計(jì)算速度上，我們選取Pilkington[10]在1997年提出的殘差記錄方式，通過下降曲線的比較來反映不同算法在三維重力反演中收斂的快慢。本文對所有算法選擇計(jì)算500次迭代，統(tǒng)計(jì)公式為:R= (g-Gρ)T(g-Gρ)/N，記錄的曲線如圖5。與傳統(tǒng)的CG算法相比，PCG算法殘差曲線明顯下降較快，表明其用較少的迭代次數(shù)即可達(dá)到相同的殘差。由圖5可以看到，迭代500次后，傳統(tǒng)的CG算法和并行的DCPG算法最小殘差能達(dá)到0.275，并行的SSOR-PCG算法能達(dá)到0.013。如同圖中黑色直線截取指示，當(dāng)反演達(dá)到同一殘差(0.04)的時(shí)候，SSOR-PCG算法只需要12次迭代計(jì)算，而DPCG算法和傳統(tǒng)CG算法分別需要46和100次迭代計(jì)算。由于系數(shù)矩陣A的變化，每次迭代計(jì)算中分解預(yù)處理矩陣也需要額外不可忽略的計(jì)算時(shí)間，因此計(jì)算效率需由迭代時(shí)間和總計(jì)算用時(shí)共同決定。表2給出了3個(gè)算法的計(jì)算用時(shí)及加速比。對于網(wǎng)格大小為20×20×10的三維重力反演計(jì)算，選取500次迭代時(shí)，并行的PCG算法相比串行的CG算法只有2.2倍。然而實(shí)際上迭代終止于同一迭代殘差(0.04)時(shí)候，相比于CG算法，并行DPCG可以達(dá)到4.7倍加速比，并行SSOR-PCG算法能獲得18.7倍的加速比。綜上，并行的SSOR-PCG預(yù)處理方法在三維重力數(shù)據(jù)反演中效果最好，反演出的密度更接近真實(shí)密度，模型更加緊致，在計(jì)算速度上可比傳統(tǒng)串行的CG算法獲得高達(dá)18.7倍的加速比，減少了計(jì)算耗時(shí)，提高了三維密度反演效率。

圖3 組合模型的空間分布圖(a)和含噪剩余重力異常(b)Fig.3 Three-dimensional distribution of synthetic model (a) and residual noise-contaminated gravity anomaly (b)

a,d為傳統(tǒng)共軛梯度算法；b,e為對角預(yù)處理共軛梯度算法；c,f為SSOR預(yù)處理共軛梯度算法。圖4 3種算法反演結(jié)果的深度切片和垂直切片F(xiàn)ig.4 Depth slices and vertical slices of results obtained by three inverse algorithms

黑色水平線表示同一迭代殘差下的迭代次數(shù)。圖5 組合模型數(shù)據(jù)反演迭代殘差Fig.5 Residual error curve for the synthetic gravity anomaly inversion

Table2ComparisonoftheruntimeandspeedupofCGandparallelPCG

方法500次迭代時(shí)間/s反演結(jié)果需要迭代次數(shù)運(yùn)行時(shí)間/s加速比CG?CPU543100108．61．0倍DPCG?GPU2504623．04．7倍SSOR?PCG?GPU242125．818．7倍

3 實(shí)際數(shù)據(jù)

為了進(jìn)一步驗(yàn)證本文方法對實(shí)測數(shù)據(jù)的有效性，我們將表現(xiàn)效果最好的并行SSOR-PCG算法應(yīng)用于美國路易斯安那州西南部Vindon鹽丘實(shí)測地面重力數(shù)據(jù)，通過快速反演來恢復(fù)地下密度分布。許多學(xué)者對該地區(qū)進(jìn)行了研究[20-22]，其大部分結(jié)果表明該鹽丘的異常主要是由巖蓋引起，其平均剩余密度為0.55 g/cm3，深度范圍為200～600 m。對原始數(shù)據(jù)去除背景場后，在研究區(qū)域中心截取3 000 m×3 000 m 的區(qū)域，反演區(qū)域被離散成20×20×20=8 000個(gè)小長方體，每個(gè)長方體的大小為150 m×150 m×50 m。圖6a給出了處理后的Vindon鹽丘的剩余布格重力異常。根據(jù)Ennen[23]在2011年的研究表明，巖蓋的密度為2.75 g/cm3，周圍圍巖(主要是頁巖和砂巖)密度為2.20 g/cm3，據(jù)此我們給出密度約束上下限為0.00～0.55 g/cm3。圖6b給出三維快速反演結(jié)果在深度為250 m的切片，其良好地指示出異常的東西向范圍約為1 300 m，南北向范圍約為900 m。圖6c給出反演結(jié)果在東西和南北向的垂直交叉剖面，在空間上可以清晰地展示出異常的范圍，深度范圍為150～600 m。圖7為Marfurt等[24]在Vinton鹽丘地區(qū)通過垂直地震剖面(VSP)成像得到的二維垂直地震剖面以及經(jīng)時(shí)間深度轉(zhuǎn)換得到的解釋層位，將反演結(jié)果與地震層位資料對比，可以看到高密度體的深度范圍與地震層位解釋中巖蓋的埋深相吻合，進(jìn)一步驗(yàn)證了本文方法的可靠性。在計(jì)算速度上，應(yīng)用并行的SSOR-PCG算法進(jìn)行三維重力快速反演達(dá)到這一結(jié)果用22次迭代，只需75.6 s。綜上，實(shí)測數(shù)據(jù)的反演結(jié)果也證實(shí)了本文并行的預(yù)處理方法可以應(yīng)用于三維重力快速反演，提高了計(jì)算效率。

a.Vinton鹽丘的異常重力數(shù)據(jù)；b.快速反演結(jié)果在深度為250 m的切片；c.通過異常體中心的垂直交叉剖面，東向?yàn)? 650 m，北向?yàn)? 050 m。圖6 Vinton鹽丘重力異常及反演結(jié)果圖Fig.6 Gravity anomaly and inversion results of Vinton Salt dome

剖面(左)；層位解釋(右)。圖7 地震測量及層位解釋結(jié)果Fig.7 Result of seimic horizon interpretaion and measurement

4 結(jié)論

1) 本文將一種并行的預(yù)處理共軛梯度算法應(yīng)用于三維重力數(shù)據(jù)的光滑反演，來快速獲得地下密度分布。并行的PCG算法與傳統(tǒng)的CG算法相比，在減少迭代次數(shù)的同時(shí)縮短了計(jì)算總用時(shí)，提高了三維密度反演的效率。

2)含噪組合模型測試結(jié)果表明：對稱逐次超松弛預(yù)處理算子效果較好，反演結(jié)果更加緊致并貼近實(shí)際模型的三維密度分布。計(jì)算效率是由較少的迭代次數(shù)和總計(jì)算時(shí)間共同評價(jià)的。實(shí)驗(yàn)表明與2.0 GHz CPU上傳統(tǒng)CG算法的串行代碼相比，本文基于NVIDIA Tesla C2050 GPU并行的SSOR-PCG算法收斂更快，獲得約19倍的加速比。說明該方法既能加快反演收斂速度，又能提高反演精度。

3)通過美國Vinton鹽丘實(shí)測重力數(shù)據(jù)的三維重力快速反演計(jì)算，得到的反演結(jié)果深度與地震層位反演結(jié)果吻合，并有著較高的效率，驗(yàn)證了該方法在實(shí)際數(shù)據(jù)三維快速密度反演應(yīng)用方面的高效性和有效性。

[1] 劉財(cái), 楊寶俊, 馮晅,等. 論油氣資源的多元勘探[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 2016，46(4):1208-1220.

Liu Cai, Yang Baojun, Feng Xuan, et al. Multivariate Exploration Technology of Hydrocarbon Resources[J]. Journal of Jilin University(Earth Science Edition), 2016, 46(4):1208-1220.

[2] Nettleton L L. Gravity and Magnetics in Oil Prospe-cting[M]. New York：McGraw-Hill Companies, 1976.

[3] Hinze W J. The Role of Gravity and Magnetic Methods in Engineering and Environmental Studies[J]. Geotechnical and Environmental Geophysics, 1990, 1:75-126.

[4] Li Y, Oldenburg D W. 3-D Inversion of Magnetic Data[J]. Geophysics, 1996, 61(2):394-408.

[5] Li Y, Oldenburg D W. 3-D Inversion of Gravity Data[J]. Geophysics, 1998, 63(1):109-119.

[6]Bell R E, Anderson R, Pratson L. Gravity Gradio-metry Resurfaces[J]. Leading Edge, 1997, 16(1):55-59.

[7] Lee J B. FALCON Gravity Gradiometer Technology[J]. Exploration Geophysics, 2001, 32(4):247-250.

[8] 黃大年, 于平, 底青云,等. 地球深部探測關(guān)鍵技術(shù)裝備研發(fā)現(xiàn)狀及趨勢[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 2012, 42(5):1485-1496.

Huang Danian, Yu Ping, Di Qingyun, et al. Development of Key Instruments and Technologies of Deep Exploration Today and Tomorrow[J]. Journal of Jilin University(Earth Science Edition), 2012,42(5): 1485-1496.

[9] Canning F X, Scholl J F. Diagonal Preconditioners for the EFIE Using a Wavelet Basis[J]. IEEE Transactions on Antennas & Propagation, 1996, 44(9):1239-1246.

[10] Pilkington M. 3D Magnetic Imaging Using Conjugate Gradients[J]. Geophysics, 1997, 62(4):1132-1142.

[11] Chen R S, Yung E K N, Chan C H, et al. App-lication of Preconditioned CG-FFT Technique to Method of Lines for Analysis of the Infinite-Plane Metallic Grating[J]. Microwave & Optical Technology Letters, 2000, 24(3):170-175.

[12] Chen R S, Yung K N, Chan C H, et al. Application of the SSOR Preconditioned CG Algorithm to the Vector FEM for 3D Full-Wave Analysis of Electromagnetic-Field Boundary-Value Problems[J]. IEEE Transactions on Microwave Theory & Techniques, 2002, 50(4):1165-1172.

[13] NVIDIA Corporation. NVIDIA CUDA Programming Guide[R]. Santa Clara CA: NVIDIA Corporation, 2008.

[14] Moorkamp M, Jegen M, Roberts A, et al. Massively Parallel Forward Modeling of Scalar and Tensor Gravimetry Data[J]. Computers & Geosciences, 2010, 36(5):680-686.

[15] Chen Z, Meng X, Guo L, et al. GICUDA: A Parallel Program for 3D Correlation Imaging of Large Scale Gravity and Gravity Gradiometry Data on Graphics Processing Units with CUDA[J]. Computers & Geosciences, 2012, 46(3):119-128.

[16]Liu G, Meng X, Chen Z. 3D Magnetic Inversion Based on Probability Tomography and Its GPU Implement[J]. Computers & Geosciences, 2012, 48(9):86-92.

[18] Zhdanov M S. GeophysicalInverse Theory and Re-gularization Problems[M].Salt Lake: Elsevier Science, 2002:51-57.

[19]Szymczyk M, Szymczyk P. Matlab and Parallel Computing[J]. Image Processing & Communications, 2014, 17(4):207-216.

[20] Oliveira V C Jr, Barbosa V C F. 3-D Radial Gravity Gradient Inversion[J]. Geophysical Journal International, 2013, 195(2):883-902.

[21] Geng M, Huang D, Yang Q, et al. 3D Inversion of Airborne Gravity-Gradiometry Data Using Cokriging[J]. Geophysics, 2014, 79(4):G37-G47.

[22] Wang T H, Huang D N, Ma G Q, et al. Improved Preconditioned Conjugate Gradient Algorithm and Application in 3D Inversion of Gravity-Gradiometry Data[J]. Applied Geophysics,2017, 14(2): 301-313.

[23]Ennen C. Structural Mapping of the Vinton Salt Dome, Louisiana, Using Gravity Gradiometry Data[J]. Seg Technical Program Expanded Abstracts, 2011, 30(1):830-835.

[24] Marfurt K J, Zhou H W, Sullivan E C. Development and Calibration of New 3-D Vector Vsp Imaging Technology: Vinton Salt Dome, LA[R]. Houston:University of Houston, 2004.