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

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

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

?

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

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

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

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

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

0 引言

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

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

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

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

1 三度體重力矢量的變分問(wèn)題

1.1 邊值問(wèn)題

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

(1)

(2)

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

(3)

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

(4)

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

1.2 變分問(wèn)題

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

(5)

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

(6)

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

(7)

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

(8)

由此可知

(9)

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

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

(10)

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

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

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

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

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

(11)

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

(12)

式中，Ui,j是六面體單元頂點(diǎn)j的待定Ui值。

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

(13)

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

(14)

式(10)中第一個(gè)積分中的第二項(xiàng)可表示為

(15)

式中：

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

(16)

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

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

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

(17)

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

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

(18)

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

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

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

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

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

由此，A的存儲(chǔ)可表示為表1。

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

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

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

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

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

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

3 模型試驗(yàn)

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

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

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

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

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

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

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

4 計(jì)算精度分析

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

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

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

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

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

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

5 結(jié)論

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

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

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

附錄 A

針對(duì)式(4)的邊值問(wèn)題，構(gòu)造泛函

(A-1)

則其變分為

(A-2)

根據(jù)格林公式

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

(A-3)

又(A-2)式的第二項(xiàng)可寫(xiě)為

(A-4)

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

(A-5)

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

(A-6)

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

(A-7)

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

[2] 李?yuàn)檴? 吳曉平, 吳星. 重力矢量及張量信息在地球重力場(chǎng)中的應(yīng)用[J]. 測(cè)繪學(xué)院學(xué)報(bào), 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é)報(bào): 地球科學(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] 駱遙, 姚長(zhǎng)利. 復(fù)雜形體重力場(chǎng)、梯度及磁場(chǎng)計(jì)算方法[J]. 地球科學(xué): 中國(guó)地質(zhì)大學(xué)學(xué)報(bào), 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ì)算及反演策略[J]. 地球物理學(xué)報(bào), 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] 張力. 坑道直流聚焦超前探測(cè)有限元數(shù)值模擬研究[D]. 長(zhǎng)沙：中南大學(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共軛梯度法及其在地電場(chǎng)計(jì)算中的應(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] 吳小平, 徐果明, 李時(shí)燦. 利用不完全Cholesky共軛梯度法求解點(diǎn)源三維地電場(chǎng)[J]. 地球物理學(xué)報(bào), 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é)基金項(xiàng)目 (BK20140844)；江蘇省地質(zhì)礦產(chǎn)勘查局科研技改項(xiàng)目(2014-KY-15)

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

10.13278/j.cnki.jjuese.201504301

P631.1

A

蔣甫玉，謝磊磊，常文凱，等. 三度體重力矢量的有限單元法正演計(jì)算.吉林大學(xué)學(xué)報(bào):地球科學(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.