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

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

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

?

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

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

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

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

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

0前言

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

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

1重力梯度張量計算

1.1重力梯度張量概念

(1)

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

(2)

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

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

(3)

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

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

▽·g=0,▽×g=0

(5)

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

(6)

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

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

(7)

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

(8)

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

在正演計算中，一些簡單規則的模型體(如球體、垂直棱柱體等)引起的重力梯度計算式可以由式(8)推算出來并直接計算；而對于復雜的組合模型或實際地質體來說，不容易推算出其計算式，可以用簡單模型去疊加模擬逼近或將其簡單化后再計算。采用有限單元法就可以將地下復雜的地質體劃分為由有限多個有限大小的單元體組成的離散結構(或稱有限網格)，有限單元之間相交的點稱為結點。然后利用每一個單元求解在計算點產生的重力梯度異常值，最后進行疊加計算就可求得地下復雜地質體在觀測平面的重力梯度異常。

為了采用有限單元離散復雜地質體，對于空間問題，常用的單元有四面體單元、長方體單元、任意六面體單元等，其中，四面體單元是最早被提出的最簡單的空間單元，如圖10所示，S(x,y,z)為單元內任意一點，S點與4個結點的連線把原四面體分割為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)

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

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

(12)

式中，權系數和積分點坐標關系[15]如表1所示。

表1　Hammer積分的權系數和積分點坐標關系Tab.1　Weights and integration point coordinates of Hammer integral

將式(8)從直角坐標系轉換到體積坐標系，見式(13)。

(13)

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

(14)

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

2模型正演及分析

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

2.1單個球體模型

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

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

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

2.2復雜組合模型

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

圖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　組合模型中各棱柱體的幾何參數Tab.2　Geometric parameter of each prism in composite pattern

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

圖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結論

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

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

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

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

參考文獻：

[1]王謙身.重力學[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]寧津生,羅志才,晁定波.衛星重力梯度測量的研究現狀及其在物理大地測量中的應用前景[J].武漢測繪科技大學學報,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]張永志,夏朝龍,王衛東,等.日本9.0級地震區重力梯度的時空分布[J].大地測量與地球動力學,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]高春春,陸洋,史紅嶺.衛星重力梯度測量在海洋科學中的應用分析[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]舒晴,周堅鑫,尹航.航空重力梯度儀研究現狀及發展趨勢[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]孟嘉春,蔡喜楣.衛星重力梯度測量及其應用前景探討[J].地球物理學報,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]汪炳柱.用樣條函數法求重力異常二階垂向導數和向上延拓計算[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]楊輝 王宜昌.復雜形體重力異常高階導數的正演計算[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]張鳳旭,孟令順,張鳳琴,等.重力位譜分析及重力異常導數換算新方法——余弦變換[J].地球物理學報,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]曹書錦,朱自強,魯光銀,等.基于單元細分H-自適應有限元全張量重力梯度正演[J].地球物理學進展,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].北京:地質出版社,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]王勖成,邵敏.有限單元法基本原理和數值方法[M].北京:清華大學出版社,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)；地質調查項目(12120113095100)

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

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

中圖分類號：P 631.1

文獻標志碼：A

DOI：10.3969/j.issn.1001-1749.2016.02.02