蔣甫玉，謝磊磊，常文凱，黃 巖，張作宏

1.河海大學(xué)地球科學(xué)與工程學(xué)院，南京 210098 2.江蘇省地質(zhì)勘查技術(shù)院，南京 210048

?

三度體重力矢量的有限單元法正演計算

蔣甫玉1，謝磊磊1，常文凱1，黃 巖2，張作宏2

1.河海大學(xué)地球科學(xué)與工程學(xué)院，南京 210098 2.江蘇省地質(zhì)勘查技術(shù)院，南京 210048

以重力位在場源內(nèi)部滿足泊松方程為依據(jù)，以重力矢量滿足第三類邊界條件為切入點，推導(dǎo)了與三度體重力矢量滿足的邊值問題相對應(yīng)的變分問題，進(jìn)而利用有限單元法實現(xiàn)了對變分問題的求解。立方體模型試驗結(jié)果表明：文中提出的新的系數(shù)矩陣存儲方式較之傳統(tǒng)方式能夠更有效地節(jié)約存儲空間，且為利用預(yù)條件共軛梯度技術(shù)更加快速地求解線性方程組提供了保障；重力矢量的計算精度與邊界長度及單元網(wǎng)格的邊長息息相關(guān)，其計算效率則主要取決于所要計算的節(jié)點總數(shù)和大型稀疏線性方程組求解算法的優(yōu)劣；一般情況下，當(dāng)單元的邊長小于場源體邊長的1/10、邊界長度大于場源體長度的7.5倍時，能夠獲得理想的結(jié)果。

變分問題；重力矢量；有限單元法；數(shù)據(jù)存儲；計算精度；三度體

0 引言

在地球物理勘探中，模型體的正演在位場異常的解釋中有著重要的意義，對其產(chǎn)生異常特征的認(rèn)識是掌握場與場源對應(yīng)關(guān)系的切入點，是進(jìn)行位場反演、地質(zhì)解釋的基礎(chǔ)，歷來深受地球物理學(xué)家的關(guān)注。重力矢量異常不僅反映重力在垂直方向上的變化，而且提供了更加豐富更加細(xì)致乃至于能反映地球時變的中高頻重力場信息[1-2]。

針對三度體重力異常的正演計算，目前主要是通過已有的解析公式直接計算、逼近目標(biāo)體的方法來實現(xiàn)。這類方法以牛頓萬有引力定律為出發(fā)點，對具有簡單形態(tài)和均勻密度分布的地質(zhì)體直接利用解析公式計算獲得；對具有復(fù)雜形態(tài)和非均勻密度分布的地質(zhì)體采用不同的方式對復(fù)雜形體進(jìn)行分割，使之轉(zhuǎn)化為一系列簡單形體的組合，計算簡單形體的重力異常，通過求和得到整個復(fù)雜形體的異常，如基于直立矩形柱體、直立多邊形棱柱體、多面體模型的三維正演技術(shù)[3-13]。一般而言，利用直立矩形柱體模擬復(fù)雜的三維非均勻密度模型是最簡單的，如果矩形體足夠小，則每個矩形體均可看成具有單一的密度分布。該方法雖然理論上簡單，但是在實際應(yīng)用中卻稍顯笨重。一方面是因為實際的地質(zhì)體通常難以用矩形柱體進(jìn)行建模；另一方面，如果與某個矩形柱體相鄰的其他柱體與該矩形柱體具有相同的密度分布，則完全沒必要將其他柱體切割成矩形柱體進(jìn)行重復(fù)計算。為改進(jìn)直立矩形柱體的上述缺點，Talwani[4]提出了直立多邊形棱柱體模型，通過對一系列無限薄的水平薄板的疊加模擬實際地質(zhì)體，每個薄板的形狀用多邊形近似。但該方法對具有任意形狀的地質(zhì)體模擬卻稍顯不足。為此，Okabe[7]提出了更加實用的針對具有任意形狀地質(zhì)體的多面體法，用一系列不同的多邊形圍成的多面體來逼近任意形狀的地質(zhì)體，其逼近的程度取決于圍成多面體多邊形的數(shù)量與頂點的選取。由此可見，基于均勻直立矩形柱體、直立多邊形棱柱體、多面體模型模擬復(fù)雜地質(zhì)形體的三維正演技術(shù)包括了兩個方面的近似：一是分割方式與實際形體的擬合程度；二是數(shù)值積分代替解析積分的近似程度。這種基于簡單規(guī)則幾何體剖分的重力異常計算方法，無論是計算精度還是計算效率均不高，易導(dǎo)致對后續(xù)的反演處理與解釋產(chǎn)生嚴(yán)重影響。

有限元法作為地球物理數(shù)據(jù)正演的方法之一，最大優(yōu)勢在于可計算任意形狀、任意密度的復(fù)雜地質(zhì)體。隨著計算機(jī)技術(shù)的發(fā)展，CAD(computer aided design)/CAE(computer aided engineering)的集成及ANSYS等大型有限元軟件提供了強(qiáng)大的有限元建模及網(wǎng)格剖分功能。對于三維變密度的復(fù)雜地質(zhì)體模型，只要將連續(xù)的求解區(qū)域離散為一組有限個、且按一定方式互相聯(lián)結(jié)在一起的單元組合體，再將不同單元賦以不同密度值，就可以模型化形狀復(fù)雜的求解域，通過疊加原理，很快就能計算出不同模型在整個空間的重力矢量。目前，基于有限單元法對三維地質(zhì)體重力矢量的正演研究還未見相關(guān)文獻(xiàn)報道。但是僅僅針對重力異常及二維常密度地質(zhì)體梯度張量的正演計算卻取得了可喜的進(jìn)展[14-18]，如Cai和Wang[15]以引力位在場源外滿足拉普拉斯方程、在場源內(nèi)滿足泊松方程以及相關(guān)邊界條件為出發(fā)點，利用有限元法能真實地再現(xiàn)復(fù)雜情況下地球物理場的分布，有效地處理不均勻介質(zhì)和求解邊界形狀復(fù)雜問題的獨特性能，實現(xiàn)了利用有限元法計算非均勻復(fù)雜形體重力異常的快速算法。其模型試驗表明，有限元法受密度分布的非均勻性影響遠(yuǎn)遠(yuǎn)小于傳統(tǒng)的數(shù)值積分法，而且在網(wǎng)格單元相同時，極大地提高了重力異常正演計算的精度；更為可貴的是相對傳統(tǒng)的數(shù)值積分法，計算效率得到了極大的提高。

本文依據(jù)重力位的泊松方程，得到與計算三度體重力矢量邊值問題對應(yīng)的變分問題，進(jìn)而利用有限元法實現(xiàn)了對三度體重力矢量的正演計算。

1 三度體重力矢量的變分問題

1.1 邊值問題

在位場勘探中，引力源一般位于地下，在地表上部的無源空間內(nèi)滿足拉普拉斯方程，在場源內(nèi)部滿足泊松方程，即：

(1)

(2)

式中：U為重力位；G為引力常數(shù)；ρ是場源體與圍巖的密度差。對式(2)分別在x，y和z方向求導(dǎo)，則有

(3)

如果給定區(qū)域Ω上的ρ分布以及邊界S上重力矢量的某類邊界條件，例如第三類邊界條件，則求解區(qū)域內(nèi)的Ui與以下的邊值問題相對應(yīng)：

(4)

式中：n為邊界S的外法向方向；α、β是常量參數(shù)。

1.2 變分問題

利用有限單元法求解邊值問題式(3)中的Ui，與之相適應(yīng)的變分問題為對泛函F(Ui)求極值(推導(dǎo)過程見附錄A)，即

(5)

需要指出的是,在對式(4)求解的過程中，參數(shù)α,β是未知的。在本文的研究中做以下近似處理：將邊界區(qū)域S取得足夠遠(yuǎn)，這時場源在邊界產(chǎn)生的重力位可看作是質(zhì)量集中于場源體質(zhì)心的質(zhì)點在邊界上的重力位，可表示為

(6)

式中：c為常數(shù)；r為質(zhì)心到邊界單元面中心的距離。則對式(6)在x,y,z方向上求偏導(dǎo)數(shù)可得

(7)

進(jìn)一步地，Ui對邊界面法線方向n的偏導(dǎo)數(shù)為

(8)

由此可知

(9)

在實際計算中，式(9)中α第二項的分母i代表x，y，z，可表示為x=rcosθcosD,y=rcosθsinD,z=rsinθ，如圖1所示。

將式(9)代入式(5)中即可得到與邊值問題式(4)相對應(yīng)的變分問題：

(10)

進(jìn)而可利用有限單元法來求解式(10)的變分問題。

圖1 計算重力矢量的區(qū)域Ω和邊界SFig. 1 Region Ω and boundary S used in computing gravity vector

2 數(shù)值實現(xiàn)

2.1 單元插值函數(shù)及積分實現(xiàn)

取如圖2a所示的六面體母單元，取8個角點為節(jié)點，其編號及坐標(biāo)如圖中所示。則可構(gòu)造如下形函數(shù)：

(11)

其中：(ξj,ηj,ζj)是j點(j=1, 2, …, 8)的坐標(biāo)；ξ,η,ζ∈[-1,1]。單元中Ui的插值函數(shù)是

(12)

式中，Ui,j是六面體單元頂點j的待定Ui值。

a. 母單元及其坐標(biāo)；b. 子單元。圖2 六面體母單元及子單元Fig.2 Mother units and subunits of hexahedron

(13)

式中：a，b，c為六面體子單元的棱邊長；x0，y0，z0是子單元的中點(圖2b)。據(jù)此，對式(10)第一個積分中的第一項計算可表示為

(14)

式(10)中第一個積分中的第二項可表示為

(15)

式中：

對式(10)中的第二個積分進(jìn)行計算，則可得到以下結(jié)果：

(16)

若單元e的底面位于邊界面上，則K2e具有如下的形式：

由此，將各個單元的節(jié)點組成的矩陣相加，得

三是社交性。 不再像音樂電臺那樣，用戶是被動且封閉的信息受傳者，在音樂的網(wǎng)絡(luò)傳播過程中創(chuàng)作者和受眾、創(chuàng)作者之間、受眾之間可以及時進(jìn)行交流。 Beats1作為蘋果音樂下的分支，間接地享受到網(wǎng)絡(luò)傳播模式帶來的良好社交性。

(17)

式中：K和P是由所有節(jié)點構(gòu)成的擴(kuò)展矩陣。故而，由式(17)對式(10)取極值，且考慮總體合成矩陣K的對稱性，有

式中，δ表示變分。由于δUi的任意性，得線性方程組：

(18)

解線性方程組(18)可得各個節(jié)點的重力矢量Ui。

2.2 數(shù)據(jù)存儲方式

對任何地球物理問題的有限元分析最終都?xì)w結(jié)為求解大型線性方程組問題，即對具有式(18)形式的方程組的求解。對于三維問題，一般情況下系數(shù)矩陣K的維數(shù)非常大。例如，將求解區(qū)域分為60×60×60的網(wǎng)格單元，則有226 981個網(wǎng)格節(jié)點，形成的矩陣K大小為226 981×226 981。如此巨大的矩陣，常規(guī)存儲方式很難實現(xiàn)其數(shù)值計算。因而系數(shù)矩陣的存儲結(jié)構(gòu)及方式成為有限元求解效率的關(guān)鍵之一。考慮到在地球物理問題中系數(shù)矩陣的稀疏性及對稱性，很多學(xué)者提出了有效的節(jié)約存儲空間的存儲方法，如半帶寬、變帶寬、coordinate storage、CRS(compressed row storage)等[19-20]。為進(jìn)一步節(jié)約存儲空間，筆者根據(jù)求解重力矢量三維有限元的系數(shù)矩陣特點，提出一種適合迭代算法的改進(jìn)一維變帶寬壓縮存儲方法。因系數(shù)矩陣的對稱性，在存儲過程中，僅存儲矩陣的下三角部分。假設(shè)矩陣A有如下的形式：

則筆者提出的改進(jìn)存儲方法可按下列步驟實現(xiàn)：

1)用一維矩陣B按行依次存儲各非零元素，非零元素的總1個數(shù)為Nd。在B的Nd+1位置，存放矩陣A的行數(shù)。

由此，A的存儲可表示為表1。

表1 系數(shù)矩陣非零元素存儲

需要指出的是，筆者提出的存儲方法一方面避免了CRS等存儲方法中需要單獨存入對角線元素位置及列號的缺點，特別當(dāng)A數(shù)組維數(shù)較大時，存儲的列號數(shù)據(jù)要遠(yuǎn)遠(yuǎn)大于本文中IA的數(shù)據(jù)。此外，就存儲的數(shù)據(jù)量而言，CRS為2Nd+N，本文為2Nd+1，較之CSR存儲方法，筆者提出的存儲方法IA的數(shù)據(jù)小，進(jìn)一步地節(jié)約了存儲空間。另一方面，在計算A的轉(zhuǎn)置矩陣時是非常方便的，只要將B中的數(shù)據(jù)按列，以IA所示的間隔數(shù)依次存儲即可得到。由此可見，利用筆者提出的存儲方法更易于得到上三角陣，這為利用預(yù)條件處理技術(shù)求解線性方程組提供了方便。

2.3 預(yù)條件處理技術(shù)

對式(18)所示的大型稀疏線性方程組，利用不完全Cholesky共軛梯度(ICCG)迭代方法能夠取得較好的效果，將其應(yīng)用于地球物理正演具有較大潛力。對于ICCG法，文獻(xiàn)[21-22]中均有詳細(xì)的描述，文中不再贅言，僅給出其迭代過程。

令r0=-P-KUi0,q0=(CCT)-1r0, 則

其中：j=0,1,2,…為迭代次數(shù)，例如Uij即表示第j次迭代的Ui值；C為下三角矩陣，是不完全Cholesky分解因子，即K≈CCT。需要指出的是，在ICCG方法中，C矩陣與K矩陣具有完全一樣的稀疏性，僅對角線元素值不同，從而能夠簡單快速地計算出C矩陣；此外，在存儲過程中，不需要另外的存儲空間來存儲C矩陣。

3 模型試驗

前文從理論上說明了利用有限元法實現(xiàn)對重力矢量計算的可行性。下面以頂面埋深為20 m、邊長為40 m的理想立方體模型為例，研究利用有限元法對重力矢量進(jìn)行正演計算的可靠性，探討網(wǎng)格剖分間距大小及邊界距離立方體中心遠(yuǎn)近對計算結(jié)果的影響。計算過程中立方體的剩余密度取為1.0 g/cm3。網(wǎng)格的節(jié)點編號如圖3所示。

Nz、Ny分別為z方向和y方向一條直線內(nèi)的節(jié)點總數(shù)。圖3 單元立方體編號示意圖Fig.3 Schematic of cell cube node number

圖4給出了利用有限單元法計算的重力矢量。由圖4可以看出異常具有明顯的對稱性，這與立方體模型的重力矢量異常分布相符，證明了計算結(jié)果的正確性。圖4中：Ux異常等值線呈橢圓形分布，南高北低，有兩個符號相反的極值中心，極值大小分別為0.986 4和-0.986 4 g.u.；而Uy異常則沿東西向呈對稱狀分布，其余特征與Ux類似，不再贅述；Uz異常等值線呈近圓形狀分布，有一個極大值，位置在立方體中心的地表投影處，大小為2.465 9 g.u.，向四周等值線逾稀疏。

誤差限為0.001時，不同立方體單元邊長及單元總數(shù)、總體系數(shù)矩陣非零元素數(shù)、占用空間、計算用時等參數(shù)見表2。表2中的計算用時均為程序在主頻為3.10 GHz、內(nèi)存為3.25 GB的微機(jī)上運行時間，邊界長度指的是2倍邊界單元的外表面到場源體中心的垂直距離。從表2中可以看出：隨著單元數(shù)及節(jié)點數(shù)的增加，總體系數(shù)矩陣中非零元素數(shù)倍增，相應(yīng)的存儲空間亦呈逐漸增大的趨勢；因采用前文中所述及的存儲方法，整體系數(shù)矩陣及其聯(lián)系數(shù)組所占用的空間并不大，當(dāng)節(jié)點總數(shù)為226 981時，所占用的存儲空間僅為23.49 Mb，一般的微機(jī)均能滿足存儲要求。此外，從ICCG迭代次數(shù)看，計算Ux的迭代次數(shù)要大于Uz，但二者差別不大，在3次左右；從計算效率看，隨計算節(jié)點總數(shù)的增加，計算用時愈來愈多，且計算Ux的用時要略大于Uz。

AB為極值中心主剖面。圖4 立方體模型重力矢量(單位：g.u.)Fig.4 Gravity vector of cube model (unit: g.u.)

單元總數(shù)節(jié)點總數(shù)總體系數(shù)矩陣非零元素數(shù)占用空間/Mb迭代次數(shù)UxUz計算用時/sUx本文算法棱柱體算法Uz本文算法棱柱體算法邊界長度201031331155610.121191.22.21.12.1200m,不1020392611181210.90181536.072.832.672.2同立方體5403689219202417.0228252781.44528.52761.64518.6邊長/m4503132651178780113.64322810942.225432.810821.625056.3立方體邊10020392611181210.90191648.598.447.197.5長5m,不150303297913916812.992320441.81025.3423.21011.8同邊界200403689219202417.0228252781.49875.62761.69859.3長度/m250503132651178780113.64322810847.630895.810799.230786.4300603226981307836123.49363242290.0109865.042100.2109400.0

較之傳統(tǒng)算法中的棱柱體法而言[15]，若將本文模型中的立方體按照文中提出的有限元單元法計算過程中所劃分的單元數(shù)量，計算同等數(shù)量的節(jié)點重力矢量，本文的計算用時要遠(yuǎn)小于棱柱體法，詳見表2。究其原因是因為在本文的方法中，假如求解區(qū)域為M個單元，場源體為N個單元，二者的比值M/N=Q為10～100，則可以同時計算出所有節(jié)點處的重力矢量，所需的積分次數(shù)則為M。而棱柱體算法中，計算任意一點的重力矢量均需要進(jìn)行N次積分運算，若求解區(qū)域的節(jié)點數(shù)為P，則所要進(jìn)行的積分總數(shù)為N·P。這就表明比值N·P/M=P/Q可以用來作為衡量棱柱體算法與本文提出的方法之間計算效率的標(biāo)準(zhǔn)。由此可見，文中提出的方法就計算效率而言，要遠(yuǎn)遠(yuǎn)高于傳統(tǒng)的棱柱體算法。

4 計算精度分析

考慮到前文中立方體模型Ux與Uy的呈轉(zhuǎn)置關(guān)系的特點，文中僅分析Ux和Uz的計算精度。為方便精度分析，選取過Ux和Uz極值中心的主剖面AB(圖4)作為分析對象。圖5給出了邊界長度為200 m，立方體單元邊長分別為20、10、5和4 m的有限元法計算結(jié)果及其與理論值的相對誤差和相對百分比誤差，其中模型的重力矢量理論值按文獻(xiàn)[23]推導(dǎo)的理論計算公式計算。由圖5a可以看出：有限元法的計算的Uz值小于理論值，而Ux總體上呈正值端小于理論值、負(fù)值端略大于理論值、0值處與理論值相等的特征；此外，隨著單元立方體邊長的減小，其與理論值的擬合程度逾來逾高。由圖5b可以看出：Uz的相對誤差絕對值在邊界處較小，而在極值處則與理論值有較大偏差，總體呈兩端平直，中部略微下凹的特征；而Ux則在總體上呈線性特征。從百分比誤差看，Ux和Uz的百分比誤差絕對值均在端部較大，而在觀測平面的0值附近，二者的百分比絕對值誤差均小于5%(圖5c)。

圖6給出了單元立方體邊長為5 m，邊界長度分別為100、150、200、250和300 m時的有限元計算結(jié)果及其與理論值的相對誤差和相對百分比誤差。由圖6a可以看出，隨著邊界長度的增加，Ux和Uz與理論值的擬合程度逐漸增高。這在相對誤差曲線(圖6b)上表現(xiàn)得尤為明顯：Uz誤差曲線在邊界長度較小時呈“V”型特征，隨著邊界長度的增加，這種特征逐漸消失，曲線呈逐漸被拉直的特征；Ux誤差曲線則呈過0點的斜率大于0的近直線特征，邊界長度越大，直線的斜率越小。從百分比誤差曲線(圖6c)看：Ux和Uz均呈現(xiàn)出倒“U”型的特征，其中Ux的“U”型曲線陡度大于Uz，即其在靠近端部的百分比誤差絕對值要大于Uz；隨著邊界長度的增加，二者的倒“U”型曲線底部越來越平直，百分比誤差絕對值越來越小。需要特別指出的是，當(dāng)邊界長度增加到200 m時，再增加邊界長度時，Ux和Uz端部的百分比誤差減小得并不多，表明當(dāng)邊界長度達(dá)到5倍場源體邊長時，邊界長度對Ux和Uz端部值的影響越來越小。為進(jìn)一步說明邊界長度對有限元法計算結(jié)果的影響，文中統(tǒng)計了在0～50 m范圍內(nèi)的百分比誤差數(shù)值，如表3所示。

表3 0～50 m范圍內(nèi)的相對百分比誤差統(tǒng)計

a.與理論值的對比；b.相對誤差；c.相對百分比誤差。圖5 不同單元立方體邊長的Uz和UxFig.5 Calculated Uz and Ux using different length of cell cube

a. 與理論值的對比；b. 相對誤差；c. 相對百分比誤差。圖6 不同邊界長度的Uz和UxFig.6 Calculated Uz and Ux using different boundary length

由以上可以看出，利用有限單元法計算重力矢量的計算精度主要受剖分網(wǎng)格單元的邊長大小及邊界長度的影響。一般情況下，邊界處的百分比誤差較大。究其原因，其一是因為在計算邊界單元的貢獻(xiàn)時，假定了場源在邊界產(chǎn)生的重力位是質(zhì)量集中于場源體質(zhì)心的質(zhì)點在邊界上的重力位；其二是因為邊界取的仍不夠遠(yuǎn)；其三則因為在本文的研究中單元立方體節(jié)點間插值函數(shù)取的是線性插值，相對精度不是太高。若要獲得理想的計算結(jié)果，一方面要考慮網(wǎng)格單元的大小、網(wǎng)格單元總數(shù)，另一方面需要權(quán)衡計算時間，這在實際應(yīng)用中尤為重要。結(jié)合本文的計算結(jié)果，筆者建議對場源體的剖分盡量細(xì)化，單元的邊長至少應(yīng)小于場源體邊長的1/10，而邊界長度則至少應(yīng)為場源體長度的7.5倍。

5 結(jié)論

1)文中提出的總體系數(shù)矩陣存儲方式相對于半帶寬、變帶寬、coordinate storage、CRS等傳統(tǒng)存儲方式能夠進(jìn)一步地節(jié)約存儲空間；且在利用ICCG法解方程組時，更易于得到系數(shù)矩陣的轉(zhuǎn)置矩陣，這為更快速地得到線性方程組的解提供了保障。

2)重力矢量正演在采用同一插值函數(shù)的前提下，其精度主要取決于單元總數(shù)和單元網(wǎng)格的尺寸大小。其計算效率則與所要計算的節(jié)點總數(shù)息息相關(guān)。此外，對大型稀疏線性方程組的求解算法的優(yōu)劣成為提高計算效率與否的關(guān)鍵。

3)利用有限單元法進(jìn)行重力矢量的正演計算，當(dāng)單元的邊長小于場源體邊長的1/10、邊界長度大于場源體長度的7.5倍時，能夠獲得較為理想的結(jié)果，除邊界附近百分比誤差較大外，其余位置的百分誤差一般小于1%。就計算效率而言，有限單元法要遠(yuǎn)遠(yuǎn)高于傳統(tǒng)的棱柱體算法。

附錄 A

針對式(4)的邊值問題，構(gòu)造泛函

(A-1)

則其變分為

(A-2)

根據(jù)格林公式

令U=δUi，V=Ui，則(A-2)式的第一項可寫為

(A-3)

又(A-2)式的第二項可寫為

(A-4)

將(A-3)和(A-4)式代入(A-2)式中，可得

(A-5)

考慮到(A-5)式中的第一項因重力矢量滿足的泊松方程，其值為0，第三項中由于在邊界處無密度體分布，故而其值亦為0，因此，有

(A-6)

將(A-1)式中的I(Ui)以及第三類邊界條件代入(A-5)式中，即可得到與(4)式表述的邊值問題相對應(yīng)的變分問題，即

(A-7)

[1] 曾華霖. 重力場與重力勘探[M]. 北京: 地質(zhì)出版社, 2005. Zeng Hualin. Gravity Field and Gravity Exploration[M]. Beijing: Geological Publishing House, 2005.

[2] 李姍姍, 吳曉平, 吳星. 重力矢量及張量信息在地球重力場中的應(yīng)用[J]. 測繪學(xué)院學(xué)報, 2005, 22(2): 94-96. Li Shanshan, Wu Xiaoping, Wu Xing. Applications of Gravity Vector and Tensor in the Gravity Field[J]. Journal of Institute of Surveying and Mapping, 2005, 22(2): 94-96.

[3] Nagy D. The Gravitational Attraction of a Right Rectangular Prism[J]. Geophysics, 1966, 31(2): 362-371.

[4] Talwani M. Computer Usage in the Computation of Gravity Anomalies[J]. Methods in Computational Physics, 1973, 13(1): 343-389.

[5] Kwok Y K. Singularities in Gravity Computation for Vertical Cylinders and Prisms[J]. Geophysical Journal International, 1991, 104(1): 1-10.

[6] Kwok Y K. Gravity Gradient Tensors Due to a Polyhedron with Polygonal Facets[J]. Geophysical Prospecting, 1991, 39(3): 435-443.

[7] Okabe M. Analytical Expressions for Gravity Anomalies Due to Homogeneous Polyhedral Bodies and Translations into Magnetic Anomalies[J].Geophysics, 1979, 44(4): 730-741.

[8] Holstein H, Ketteridge B. Gravimetric Analysis of Uniform Polyhedral[J]. Geophysics, 1996, 61(2): 357-364.

[9] Commer M. Three-Dimensional Gravity Modelling and Focusing Inversion Using Rectangular Meshes[J]. Geophysical Prospecting, 2011, 59(5): 966-979.

[10] Li X, Chouteau M. Three-Dimensional Gravity Mo-deling in All Space[J]. Surveys in Geophysics, 1998, 19(4): 339-368.

[11] 蔣甫玉, 高麗坤, 黃麟云. 油氣模型的重力梯度張量研究[J]. 吉林大學(xué)學(xué)報: 地球科學(xué)版, 2011, 41(2): 545-551. Jiang Fuyu, Gao Likun, Huang Linyun. Study on Gravity Gradient Tensor of the Oil-Gas Model[J]. Journal of Jilin University: Earth Science Edition, 2011, 41(2): 545-551.

[12] 駱遙, 姚長利. 復(fù)雜形體重力場、梯度及磁場計算方法[J]. 地球科學(xué): 中國地質(zhì)大學(xué)學(xué)報, 2007, 32(4): 517-522. Luo Yao, Yao Changli. Forward Method for Gravity, Gravity Gradient and Magnetic Anomalies of Complex Body[J]. Earth Science: Journal of China University of Geosciences, 2007, 32(4): 517-522.

[13] 陳召曦, 孟小紅, 郭良輝, 等. 基于GPU并行的重力、重力梯度三維正演快速計算及反演策略[J]. 地球物理學(xué)報, 2012, 55(12): 4069-4077. Chen Zhaoxi, Meng Xiaohong, Guo Lianghui, et al. Three Dimensional Fast Forward Modeling and the Inversion Strategy for Large Scale Gravity and Gravimetry Data Based on GPU[J]. Chinese Journal of Geophysics, 2012, 55(12): 4069-4077.

[14] Dave A M, Matthew G K. Optimal, Scalable Forward Models for Computing Gravity Anomalies[J]. Geophysical Journal International, 2011, 187(1): 161-177.

[15] Cai Yongen, Wang Chiyun. Fast Finite-Element Cal-culation of Gravity Anomaly in Complex Geolo-gical Regions[J]. Geophysical Journal International, 2005, 162(3): 696-708.

[16] Cruz F A, Knepley M G, Barba L A. PetFMM: A Dynamically Load-Balancing Parallel Fast Multipole Library[J]. International Journal for Numerical Methods in Engineering, 2010, 85(4): 403-428.

[17] Farquharson C G,Mosher C R W.Three-Dimensional Modelling of Gravity Data Using Finite Differences[J]. Journal of Applied Geophysics, 2009, 68(3): 417-422.

[18] 朱自強(qiáng), 曾思紅, 魯光銀, 等. 二度體的重力張量有限元正演模擬[J]. 物探與化探, 2010, 34(5): 668-672. Zhu Ziqiang, Zeng Sihong, Lu Guangyin, et al. Finite Element Forward and Simulation of the Two Dimensional Gravity Gradient Tensor[J]. Geophysical and Geochemical Exploration, 2010, 34(5): 668-672.

[19] 徐世浙. 地球物理中的有限單元法[M]. 北京: 科學(xué)出版社, 1994. Xu Shizhe. Finite Element Method in Geophysics[M]. Beijing: Science Press, 1994.

[20] 張力. 坑道直流聚焦超前探測有限元數(shù)值模擬研究[D]. 長沙：中南大學(xué), 2011. Zhang Li. Research on Numerical Simulation for DC Focusing Advanced Detection in Tunnel with Finite Element Method[D]. Changsha: Central South University, 2011.

[21] 吳小平, 徐果明. 不完全Cholesky共軛梯度法及其在地電場計算中的應(yīng)用[J]. 石油地球物理勘探, 1998, 33(1): 89-94. Wu Xiaoping, Xu Guoming. Incomplete Cholesky Conjugate Gradient Method and Its Application in Geoelectric Field Computation[J]. OGP, 1998, 33(1): 89-94.

[22] 吳小平, 徐果明, 李時燦. 利用不完全Cholesky共軛梯度法求解點源三維地電場[J]. 地球物理學(xué)報, 1998, 41(6): 848-855. Wu Xiaoping, Xu Guoming, Li Shican. The Calculation of Three-Dimensional Geoelectric Field of Point Source by Incomplete Cholesky Conjugate Gradient Method[J]. Chinese Joural of Geophysics, 1998, 41(6): 848-855.

[23] Nagy D, Papp G, Beneded J. The Gravitational Potential and Its Derivatives for the Prism[J]. Journal of Geodesy, 2000, 74: 552-560.

Forward Calculation of Three Dimensional Gravity Vector Using Finite Element Method

Jiang Fuyu1, Xie Leilei1, Chang Wenkai1, Huang Yan2, Zhang Zuohong2

1.SchoolofEarthSciencesandEngineering,HohaiUniversity,Nanjing210098,China2.GeologyExplorationTechnologyInstituteofJiangsuProvince,Nanjing210048,China

Variational problem of three dimensional gravity vector was deduced to meet the boundary value based on Poisson equation and the third boundary condition, and the solution of variational problem is further implemented by using the finite element method. The results of the cubic model test show that the proposed new coefficient matrix storage strategy is more effective to save storage space than a traditional approach; this, in turn, makes it possible to quickly solve liner equations by using the preconditioned conjugate gradient technology. The calculation precision of the gravity vector is closely related to the boundary length and unit grid; while the computational efficiency mainly depends on the total number of nodes and the algorithm used in solving a large sparse system of linear equation. In general, when the length of unit grid is less than 1/10 of the body length, and the boundary length is greater than 7.5 times of the length of the source, a desired result can be achieved.

variational problem; gravity vector; finite element method; data storage; calculation precision; three dimensional body

10.13278/j.cnki.jjuese.201504301.

2014-10-30

江蘇省自然科學(xué)基金項目 (BK20140844)；江蘇省地質(zhì)礦產(chǎn)勘查局科研技改項目(2014-KY-15)

蔣甫玉(1981--)，男，講師，博士，主要從事固體地球物理學(xué)研究，E-mail:jiangfy@hhu.edu.cn。

10.13278/j.cnki.jjuese.201504301

P631.1

A

蔣甫玉，謝磊磊，常文凱，等. 三度體重力矢量的有限單元法正演計算.吉林大學(xué)學(xué)報:地球科學(xué)版,2015,45(4):1217-1226.

Jiang Fuyu, Xie Leilei, Chang Wenkai, et al. Forward Calculation of Three Dimensional Gravity Vector Using Finite Element Method.Journal of Jilin University:Earth Science Edition,2015,45(4):1217-1226.doi:10.13278/j.cnki.jjuese.201504301.