彭江英，吳　興，韓文文，柳康偉

(成都理工大學(xué) a.地球物理學(xué)院,b.信息科學(xué)與技術(shù)學(xué)院,成都　610059)

?

基于有限單元法的重力梯度張量正演計算與分析

彭江英a，吳興a，韓文文b，柳康偉a

(成都理工大學(xué) a.地球物理學(xué)院,b.信息科學(xué)與技術(shù)學(xué)院,成都610059)

摘要：由于重力梯度張量測量相對于重力異常測量有許多優(yōu)點(diǎn)，因此對重力梯度的研究十分必要。在重力梯度各張量的正演計算中，復(fù)雜地質(zhì)體的計算式難以直接推導(dǎo)，而利用有限元技術(shù)在復(fù)雜形體體積積分的優(yōu)勢，可以較為快速和簡便地進(jìn)行復(fù)雜模型的重力梯度全張量正演計算。通過模型的建立和正演計算，分析了球體模型梯度張量異常的平面和剖面特征，以及重力梯度張量與球體模型位置的規(guī)律；最后進(jìn)一步利用復(fù)雜模型的重力梯度正演模擬，說明了重力梯度識別地下地質(zhì)體位置的優(yōu)勢和不足。

關(guān)鍵詞：重力梯度全張量；有限單元法；正演計算；異常識別

0前言

重力資料是區(qū)域地質(zhì)研究的基礎(chǔ)，不僅對于解析區(qū)域大地構(gòu)造格架、板塊運(yùn)動、形成與演化等大尺度的基礎(chǔ)性研究有參考價值，而且對于斷裂分布、盆地劃分、沉積厚度等重要油氣地質(zhì)問題也有相當(dāng)程度的指導(dǎo)意義。地球重力場可用重力位(重力標(biāo)量)、重力位梯度(重力矢量)和位的高階導(dǎo)數(shù)來表示，重力矢量各分量的梯度是重力位的二階導(dǎo)數(shù)[1]。重力梯度研究和測量開始于19世紀(jì)七十年代，重力梯度張量包含了對于大地測量和地球物理學(xué)中許多問題都十分重要的局部重力場信息[2]，它被廣泛應(yīng)用于大地測量學(xué)[3]、地球物理學(xué)、地質(zhì)學(xué)、地震學(xué)[4]、海洋學(xué)[5]、航空和空間科學(xué)[3,5-7]等領(lǐng)域。與傳統(tǒng)的重力異常測量相比，重力梯度張量測量有更高的靈敏度、短波信息豐富，包含多個信息不同的重力梯度張量分量，對慣性加速度不敏感等。

由于重力梯度測量相對于重力異常測量有上述各種優(yōu)越性，對重力梯度的研究顯得十分重要和必要的。重力梯度張量測量在國外已經(jīng)十分成熟并進(jìn)入了商業(yè)應(yīng)用階段，而在國內(nèi)還處于起步階段，但重力梯度張量處理和解釋已經(jīng)得到了開展。除了用重力梯度儀直接觀測重力梯度張量，還可由重力異常數(shù)據(jù)經(jīng)過一定的算法求得重力梯度張量，更直接的辦法是直接用模型正演來計算簡單模型研究復(fù)雜地質(zhì)體。國內(nèi)、外學(xué)者已經(jīng)研究和提出了眾多重力梯度的求取算法，例如樣條函數(shù)插值法[8]、差商方法[9]、頻率域中的傅里葉變換法[10]和余弦變換法[11]以及有限單元方法[12]等。作者介紹了重力梯度張量概念，利用有限元技術(shù)推導(dǎo)出梯度張量的計算式，計算并分析了球體模型張度張量異常的平面和剖面特征，總結(jié)了重力張量與球體模型位置關(guān)系的規(guī)律，最后進(jìn)一步建立復(fù)雜模型，通過計算與分析認(rèn)識到重力梯度在識別地下地質(zhì)體位置的優(yōu)勢。

1重力梯度張量計算

1.1重力梯度張量概念

(1)

式中，G為萬有引力常數(shù)。若對引力位U求一階導(dǎo)數(shù)，則可得到重力三分量，如式(2)所示。

(2)

圖1　笛卡兒坐標(biāo)系下地質(zhì)體空間位置示意圖Fig.1　Spatial location of geological body in Cartesian coordinates

由于重力方向就是Z軸方向，所以所謂的重力異常就是引力位在Z方向的一階導(dǎo)數(shù)，如式(3)所示。

(3)

而對引力位U求二階導(dǎo)數(shù)，則重力梯度九分量(又稱重力梯度全張量)可表示為式(4)。

由于計算點(diǎn)在地質(zhì)體外部，g滿足自由空間的微分方程，見式(5)。

▽·g=0,▽×g=0

(5)

將式(2)代入式(5)，可以得到重力梯度各分量之間的關(guān)系如式(6)所示。

(6)

由式(6)可知，重力梯度全張量中只有5個獨(dú)立的分量(由此可推導(dǎo)出其余分量)，見式(7)。

(Uyx,Uxz)=(Uzx,Uyz)=Uzy

(7)

將式(1)代入式(4)，可推導(dǎo)出重力梯度全張量的表達(dá)式[14]如式(8)所示。

(8)

1.2重力梯度張量的有限元計算

在正演計算中，一些簡單規(guī)則的模型體(如球體、垂直棱柱體等)引起的重力梯度計算式可以由式(8)推算出來并直接計算；而對于復(fù)雜的組合模型或?qū)嶋H地質(zhì)體來說，不容易推算出其計算式，可以用簡單模型去疊加模擬逼近或?qū)⑵浜唵位笤儆嬎?。采用有限單元法就可以將地下?fù)雜的地質(zhì)體劃分為由有限多個有限大小的單元體組成的離散結(jié)構(gòu)(或稱有限網(wǎng)格)，有限單元之間相交的點(diǎn)稱為結(jié)點(diǎn)。然后利用每一個單元求解在計算點(diǎn)產(chǎn)生的重力梯度異常值，最后進(jìn)行疊加計算就可求得地下復(fù)雜地質(zhì)體在觀測平面的重力梯度異常。

為了采用有限單元離散復(fù)雜地質(zhì)體，對于空間問題，常用的單元有四面體單元、長方體單元、任意六面體單元等，其中，四面體單元是最早被提出的最簡單的空間單元，如圖10所示，S(x,y,z)為單元內(nèi)任意一點(diǎn)，S點(diǎn)與4個結(jié)點(diǎn)的連線把原四面體分割為4個小四面體，用V、V1、V2、V3、V4分別表示四面體1234、S234、S341、S412、S123的體積，令

(9)

(10)

圖2　四面體有限單元示意圖Fig.2　Sketch map of tetrahedron elemnt

由于V=V1+V〗2+V3+V4，有

L1+L2+L3+L4=1

(11)

在四面體單元中，采用的是體積坐標(biāo)(L1,L2,L3)，其Hammer積分形式[15]如式(12)所示。

F(b,b,a,b)+F(b,b,b,a)]

(12)

式中，權(quán)系數(shù)和積分點(diǎn)坐標(biāo)關(guān)系[15]如表1所示。

表1　Hammer積分的權(quán)系數(shù)和積分點(diǎn)坐標(biāo)關(guān)系Tab.1　Weights and integration point coordinates of Hammer integral

將式(8)從直角坐標(biāo)系轉(zhuǎn)換到體積坐標(biāo)系，見式(13)。

(13)

式中，|J|為坐標(biāo)映射的雅可比行列式[15]，見式(14)。

(14)

由上述計算式可推算出線性四面體網(wǎng)格的重力梯度全張量表達(dá)式[12]，見式(15)。

2模型正演及分析

對較為常見的典型三度體(球體、水平圓柱體和棱柱體等)進(jìn)行正演模擬計算，得出重力梯度張量異常曲線，并分析模型的重力梯度張量各個分量的曲線特征，同時分析所建立模型的不同體積元參數(shù)(如板狀體的上頂邊埋深及板狀體厚度、傾角、球體的半徑及球心的埋深、水平圓柱體的截面半徑、中軸線埋深、水平延伸長度、棱柱體的各邊長及頂界面埋深等)，對重力梯度張量異常曲線的響應(yīng)特征及變化規(guī)律。

2.1單個球體模型

對球體模型重力梯度張量進(jìn)行討論時，設(shè)球體半徑為r，球心埋深為h，剩余密度值均取為1 g/cm3。對不同球體半徑及不同球心埋深所確定的模型進(jìn)行正演計算，并對其重力梯度張量異常的特征進(jìn)行對比分析。對于不同埋深、不同大小的球體模型，其重力梯度張量(Uxx、Uyy、Uzz、Uxy、Uxz和Uyz)等值線如圖3—圖5所示，圖3—圖5中黑色虛線為零值線，紅色虛線表示球體在水平面上的投影。

從重力梯度張量異常圖形可以歸納其各自的特點(diǎn)，由于目標(biāo)體的模型具有對稱性，其重力梯度張量的各分量等值線也具有對稱性，既分別關(guān)于X、Y軸對稱，也關(guān)于中心原點(diǎn)對稱。梯度張量Uxx和Uyy分量等值線上有三個極值點(diǎn)，其中球心在水平面投影點(diǎn)(記為“O點(diǎn)”)處為極小值，等值線呈類橢圓形狀向外擴(kuò)散，越靠近O點(diǎn)等值線變得越密，呈現(xiàn)出明顯的梯級帶；軸線方向上遠(yuǎn)離O點(diǎn)異常值逐漸變大，并出現(xiàn)兩個形態(tài)大小對稱的異常極大值；零等值線呈靠近O點(diǎn)彎曲的形態(tài)。梯度張量Uxy分量平面等值線形態(tài)為“四葉草”形，等值線零值在X、Y軸上，曲線在靠近X、Y軸時出現(xiàn)明顯的梯級帶，梯度張量等值線有四個極值點(diǎn)，在1、3象限取得異常極大值，在2、4象限取得異常極小值，在極值點(diǎn)沿半徑延長線像遠(yuǎn)離原點(diǎn)方向上，呈發(fā)散的形態(tài)，且等值線變稀疏。梯度張量Uxz和Uyz分量等值線形態(tài)成“扇貝”形，等值線零值在X軸或Y軸上，平面等值線有一正一負(fù)兩個極值，從兩個極值點(diǎn)靠近軸線的等值線變密，呈現(xiàn)明顯的異常梯級帶，遠(yuǎn)離O點(diǎn)方向，等值線變稀疏，呈波紋擴(kuò)散形態(tài)。梯度張量Uzz為以O(shè)點(diǎn)為中心的一系列同心圓，在O點(diǎn)處取得極大值，零等值線范圍大于球體在水平面投影的邊界，從外側(cè)到O點(diǎn)張量等值線先變密，在球體投影邊界附近呈現(xiàn)明顯的梯級帶，繼續(xù)向O心靠近等值線開始變稀疏。

圖3　球體重力梯度張量等值線圖(r=150 m，h=300 m)Fig.3　Contour map of gravity gradient tensors base on sphere(r=150 m,h=300 m)(a)Uxx;(b)Uyy;(c)Uzz;(d)Uxy;(e)Uxz;(f)Uyz

圖4　球體重力梯度張量等值線圖(r=200 m，h=300 m)Fig.4　Contour map of gravity gradient tensors base on sphere(r=20 m,h=300 m)(a)Uxx;(b)Uyy;(c)Uzz;(d)Uxy;(e)Uxz;(f)Uyz

圖5　球體重力梯度張量等值線圖(r=150 m，h=400 m)Fig.5　Contour map of gravity gradient tensors base on sphere(r=150 m,h=400 m)(a)Uxx;(b)Uyy;(c)Uzz;(d)Uxy;(e)Uxz;(f)Uyz

為了分析梯度張量的各個分量等值線與球體埋深和半徑的關(guān)系，選取不同的模型進(jìn)行實(shí)驗，并以零值線和極值點(diǎn)為參照進(jìn)行分析。圖3和圖4顯示了球心埋深相同、球體半徑不同時重力梯度張量的等值線特征，對比各個分量曲線可以發(fā)現(xiàn)：當(dāng)球體半徑改變時，等值線的零值線和極值點(diǎn)的位置不變，只有異常幅度(極值點(diǎn)的大小)隨著半徑增大而增大。圖3和圖5顯示了球體半徑相同、球心埋深不同時重力梯度張量的等值線特征，對比各個分量曲線可以發(fā)現(xiàn)：當(dāng)球心埋深增大時，除了異常幅度減小外，各個分量(除Uzz外)的極值點(diǎn)間距隨之增大，而Uxx、Uyy和Uzz分量的零值線在主剖面(過O點(diǎn)的剖面)上間距也在增大。通過上述分析，說明了球體模型的重力梯度等值線對反映球體的埋深十分靈敏，而對球體的規(guī)模大小(半徑大小)的反映卻不太明顯。

2.2復(fù)雜組合模型

為了模擬實(shí)際地質(zhì)情況，反映重力梯度全張量在識別地下復(fù)雜地質(zhì)體位置的可靠性，設(shè)計了一個組合復(fù)雜模型，該模型由多個不同埋深、不同大小和不同密度的垂直棱柱體組成，模型具體參數(shù)見表2，圖8(a)和圖8(b)分別顯示了組合模型位置在XOY平面和YOZ平面的投影。圖8(c)為組合模型的重力異常圖，從圖8中可以看到，模型A1和A2形成了獨(dú)立的封閉異常，且異常范圍和幅度較大，很容易辨別；疊加的模型B1和B2受到模型A1和A2的影響只是引起了等值線小幅度的向內(nèi)彎曲，不能形成封閉曲線，很難分辨出這兩個模型的存在；模型B3由于距離模型A1和A2較遠(yuǎn)，能夠形成封閉曲線，較易識別；而模型C1、C2和C3由于規(guī)模較小，在模型A1和A2的影響下，在等值線上呈現(xiàn)小范圍的異常凸點(diǎn)，不易識別。

圖6　不同埋深的梯度張量主剖面曲線Fig.6　Main sections of gravity gradient tensors base on different buried depth(a)Uxx(Uyy)；(b)Uzz；(c)Uxz(Uyz)；(d)Uxy

圖7　不同半徑的梯度張量主剖面曲線Fig.7　Main sections of gravity gradient tensors base on different radii(a)Uxx(Uyy)；(b)Uzz；(c)Uxz(Uyz)；(d)Uxy

圖8　組合模型在XOY平面和YOZ平面的投影以及重力異常等值線Fig.8　Projections of composite pattern in XOY plane and YOZ plane,and gravity anomaly contour of composite pattern(a)XOY;(b)YOZ;(c)重力異常等值線

表2　組合模型中各棱柱體的幾何參數(shù)Tab.2　Geometric parameter of each prism in composite pattern

應(yīng)用重力梯度全張量正演計算，可得到組合模型的各個梯度分量異常，如圖9所示。由于重力梯度張量對高頻信息最為靈敏[1]，所以其對重力異常中較難識別的規(guī)模小、異常幅度大的高頻信號識別效果十分明顯。從圖9中重力梯度各個分量中可以看到，規(guī)模最小、不易識別的模型C1、C2和C3在圖9中以模型梯度張量的特征獨(dú)立的呈現(xiàn)了出來，能夠很地識別其位置；除了重力異常中獨(dú)立的B3異常外，模型B1和B2也能在重力梯度張量中呈現(xiàn)獨(dú)立的特征形態(tài)；而在重力異常中最容易識別的A1和A2異常，在重力梯度張量中雖能基本呈現(xiàn)應(yīng)有的梯度張量特征，異常曲線卻受到了其他異常的干擾而變得畸形、影響識別。上述分析說明，重力梯度全張量在識別復(fù)雜組合模型異常體的位置方面效果明顯，尤其是對淺部的、規(guī)模較小的異常識別效果更加突出，體現(xiàn)了重力梯度張量相對于重力異常的優(yōu)越性；但這也反映了重力梯度全張量在實(shí)際資料應(yīng)用中的一個弊端，那就是由于重力梯度全張量對高頻或甚高頻信息十分靈敏，在實(shí)際重力資料處理中噪聲的存在對張量分析效果的影響也會很大。

圖9　組合模型的重力梯度張量等值線圖Fig.9　Contour map of gravity gradient tensors base on composite pattern(a)Uxx；(b)Uyy；(c)Uzz；(d)Uxy；(e)Uxz；(f)Uyz

3結(jié)論

1)針對復(fù)雜模型正演計算的復(fù)雜性，利用有限單元法能很好地運(yùn)用離散網(wǎng)格單元的疊加計算進(jìn)行重力梯度全張量正演，并且有較高的精度。

2)在球體模型參數(shù)的識別方面，重力梯度全張量能夠很好地識別模型水平面投影中心位置及中心埋深，對球體規(guī)模(半徑大小)的識別不太明顯。

3)重力梯度全張量能夠減小模型體之間的相互干擾和影響，能很好地識別掩蓋在埋深和規(guī)模較大的模型異常中的埋深和規(guī)模較小的模型異常，壓制背景場而突出局部異常。

4)由于重力梯度全張量對高頻信息很敏感，因此噪聲的存在對實(shí)際資料處理和解釋帶來很大的影響，不利于噪聲較大的實(shí)際資料處理，這也是重力梯度需要進(jìn)一步研究和改善的一個方面。

參考文獻(xiàn)：

[1]王謙身.重力學(xué)[M].北京:地震出版社,2003.

WANG Q S.Gravitology [M].Beijing:Seismological press,2003.(In Chinese)

[2]JAMES WHILE,ANDREW JACKSON,DIRK SM-IT,et al.Spectral analysis of gravity gradiometry profiles [J].Geophysics,2006,71(1):11-22.

[3]寧津生,羅志才,晁定波.衛(wèi)星重力梯度測量的研究現(xiàn)狀及其在物理大地測量中的應(yīng)用前景[J].武漢測繪科技大學(xué)學(xué)報,1996,21(4):309-314.

NING J S,LUO Z C,CHAO D B.The present situation on satellite gravity gradiometry and its vistas in the application of physical geodesy [J].Journal of Wuhan Technical University of Surveying and Mapping,1996,21(4):309-314.(In Chinese)

[4]張永志,夏朝龍,王衛(wèi)東,等.日本9.0級地震區(qū)重力梯度的時空分布[J].大地測量與地球動力學(xué),2013,33(6):2-4.

ZHANG Y Z,XIA C L,WANG W D,et al.Temporal and space contribution of gravity gradients in Japan Mw9.0 earthquake area [J].Journal of geodesy and geodynamics,2013,33(6):2-4.(In Chinese)

[5]高春春,陸洋,史紅嶺.衛(wèi)星重力梯度測量在海洋科學(xué)中的應(yīng)用分析[J].海洋測繪,2012,32(5):77-81.

GAO C C,LU Y,SHI H L.Application analysis of satellite gravity gradiometry in ocean science [J].Hydrographic surveying and charting,2012,32(5):77-81.(In Chinese)

[6]舒晴,周堅鑫,尹航.航空重力梯度儀研究現(xiàn)狀及發(fā)展趨勢[J].物探與化探,2007,31(6):485-488.

SHU Q,ZHOU J X,YIN H.Present research situation and development trend of airborne gravity gradiometer [J].Geophysical &geochemical exploration,2007,31(6):485-488.(In Chinese)

[7]孟嘉春,蔡喜楣.衛(wèi)星重力梯度測量及其應(yīng)用前景探討[J].地球物理學(xué)報,1991,34(3):369-376.

MENG J C,CAI X M .Approach on satellite gravity gradiometry and its vistas of applications [J].Chinese journal of geophysics,1991,34(3):369-376.(In Chinese)

[8]汪炳柱.用樣條函數(shù)法求重力異常二階垂向?qū)?shù)和向上延拓計算[J].石油地球物理勘探,1996,31(3):415-422.

WANG B Z .Computing the vertical second derivative and upward continuation of gravity anomaly by spline function method [J].Oil geophysical prospecting,996,31(3):415-422.(In Chinese)

[9]楊輝 王宜昌.復(fù)雜形體重力異常高階導(dǎo)數(shù)的正演計算[J].石油地球物理勘探,1998,33(2):278-283.

YANG H,WANG Y C .Forward computation of high order derivative of gravity anomaly relating to comlplex geologic body [J].Oil geophysical prospecting,1998,33(2):278-283.(In Chinese)

[10]KEVIN L.MIEKUS,JUAN HOMERO HINOJOSA.The complete gravity gradient tensor derived from the vertical component of gravity:a fourier transform technique [J].Journal of applied geophysics,2001,46:159-174.

[11]張鳳旭,孟令順,張鳳琴,等.重力位譜分析及重力異常導(dǎo)數(shù)換算新方法——余弦變換[J].地球物理學(xué)報,2006,49(1):244-248.

ZHANG F X,MENG L S,ZHANG F Q,et al.A new method for spectral analysis of the potential field and conversion of derivative of gravity-anomalies:cosine transform [J].Chinese journal of geophysics,2006,49(1):244-248.(In Chinese)

[12]曹書錦,朱自強(qiáng),魯光銀,等.基于單元細(xì)分H-自適應(yīng)有限元全張量重力梯度正演[J].地球物理學(xué)進(jìn)展,2010,25(3):1015-1023.

CAO S J,ZHU Z Q,LU G Y,et al.Forward modelling of full gravity gradient tensors based H-Adaptive mesh refinement [J].Progress in geophy,2010,25(3):1015-1023.(In Chinese)

[13]曾華霖.重力場與重力勘探[M].北京:地質(zhì)出版社,2005.

ZENG H L.Gravity field and gravity exploration [M].Beijing:Geological publishing house,2005.(In Chinese)

[14]MICHAEL S.ZHDANOV,ROBERT ELLISZ,SOUVIK MUKHERJEE.Three dimensional regularized focusing inversion of gravity gradient tensor component data [J].Geophysics,2004,69:925-937.

[15]王勖成,邵敏.有限單元法基本原理和數(shù)值方法[M].北京:清華大學(xué)出版社,1997.

WANG M C,SHAO M.Basic principle and numerical method of finite element method [M].Beijing:Tsinghua University press,1997.(In Chinese)

Forward calculation and analysis of gravity gradient tensors based on finite element method

PENG Jiang-yinga,WU Xinga,HAN Wen-wenb,LIU Kang-weia

(Chengdu University of Technology a.College of Geophysics,b.College of Geophysics Chengdu610059,China)

Abstract:Gravity gradient anomaly has many advantages relative to gravity anomaly,it is necessary to study gravity gradient tensors.In forward modeling of gravity gradient tensors,it is difficult to derive mathematical formulas of complex geological bodies.But finite element method has an advantage in volume integral,it can be used to quickly and easily do forward calculation for full tensors gravity gradient of complex geological bodies.By establishing models and forward calculating,we analyzed plane and section characteristics of gravity gradient anomaly of sphere model,and analyzed the relationship between gravity gradients and position of sphere model.Finally we utilized gravity gradients forward simulation of complex geological bodies,and illustrated the advantages and disadvantages of gravity gradients in finding out the position of geologic body.

Key words:full tensor gravity gradient;finite element method;forward calculation;anomaly recognition

收稿日期：2015-04-10改回日期：2015-05-25

基金項目:國家科技重大專項“十二五”(2011ZX05025-001-08)；地質(zhì)調(diào)查項目(12120113095100)

作者簡介：彭江英(1990-)，女，碩士，主要研究方向為深部地球物理與地球動力學(xué)，E-mail：wylxing@foxmail.com。

文章編號：1001-1749(2016)02-0151-09

中圖分類號：P 631.1

文獻(xiàn)標(biāo)志碼：A

DOI：10.3969/j.issn.1001-1749.2016.02.02