基于旋轉橢球體的等軸狀場源體重力異常場模型構建方法

劉繁明，劉惠敏,2，荊 心

（1.哈爾濱工程大學 自動化學院，哈爾濱 150001；2.青島農業大學 機電工程學院，青島 266109）

基于旋轉橢球體的等軸狀場源體重力異常場模型構建方法

劉繁明1，劉惠敏1,2，荊 心1

（1.哈爾濱工程大學 自動化學院，哈爾濱 150001；2.青島農業大學 機電工程學院，青島 266109）

在中心埋藏淺、測量精度高的場合，采用球體作為有限大小、近于等軸狀的場源體的模型進行重力異常正演會造成較大誤差。提出一種基于旋轉橢球體的等軸狀場源體重力異常模型。通過改變長、短軸長度以匹配場源體形態，降低正演誤差。在已有繞短軸旋轉型旋轉橢球體的引力位的基礎上，推導旋轉橢球體完整的重力異常基本表達式。將結果與球體重力異常比較，分析其間的特征差異，驗證了結果的正確性。

旋轉橢球體；重力異常場；正演模型；球體；等軸狀場源

基于微重力測量技術探測自然形成的地下空穴已經成為城市、交通、廢棄煤礦、大壩和核電站等重要設施地面塌陷防治和穩定性檢測的一種重要的非破壞性手段[1-6]。通常，地下空穴可視為等軸狀分布，測量往往集中在面積、深度以及網格間距都較小的區域內。按照現有的重力場正演理論，簡單地將地下空洞等效成球體、棱柱體、圓柱體等常用規則形體，在高精度和小范圍的測量條件下，勢必造成正演結果與實際情況不符，誤差所占比例高，不利于做出正確的解釋。旋轉橢球體也是一種等軸狀規則形體，可分成繞短軸旋轉型和繞長軸旋轉型兩種基本類型。由于通過改變長、短軸的長度可以得到無限個不同的旋轉橢球體，利于更準確地匹配實際地下空穴的狀態，獲得更接近實際的正演結果，因此可以將旋轉橢球體用于地下空穴正演模型進行重力異常場分析。

另一方面，隨著重力梯度測量技術的不斷提升，針對重力梯度測量不受外界信號干擾、不受測區環境影響、不向外界發射信號的無源特征[7-10]，孫嵐等提出了基于重力梯度探測的潛艇探測方法，將潛艇假設為一種繞長軸旋轉的旋轉橢球體并進行重力梯度正演計算[7]，然而由于只針對特定尺寸驗證了該方法的可行性，并沒有給出具體的計算表達式，無法普遍應用于其它尺寸的潛艇。推導繞長軸旋轉型旋轉橢球體重力梯度異常場正演計算表達式可以為該方法的實現提供理論依據。

目前，關于旋轉橢球體重力場理論只局限于“地球橢球體”這類典型的繞短軸旋轉型旋轉橢球體，主要用于地球形狀學和大地測量學研究領域，并且在已有文獻中，只推導出其引力位及一、二階導數，無法滿足上述兩個方面對旋轉橢球體的重力場的分析要求。因此，本文在現有公式的基礎上，進一步推導出繞長軸選旋轉型旋轉橢球體的重力場正演計算公式。另外，雖然目前尚沒有用于直接測定重力位三階導數的探測儀器，但是由于重力位三階導數可以由重力異常換算得到，并且它對淺而小的異常源較重力位二階導數反映更為突出，因此本文將兩種旋轉橢球體的引力位進一步推導至三階導數。

1 兩種類型的旋轉橢球體

在右手直角坐標系中，XOY平面為觀測平面，Z軸垂直向下。均質且密度為ρ的旋轉橢球體的球心位于(0, 0, h)處，坐標原點為其在地面上的投影點。規定繞短軸旋轉橢球體為旋轉橢球體I，如圖1(a)所示，其球面方程為

規定繞長軸旋轉的旋轉橢球體為旋轉橢球體II，如圖1(b)所示，其球面方程為

當a=c時，為球體。

圖1 兩種類型的旋轉橢球體Fig.1 Two types of the rotational ellipsoid

2 旋轉橢球體重力異常基本表達式

盡管旋轉橢球體屬于簡單形體，但其引力位的計算與球體、圓柱體、矩形棱柱體等規則形體不同。首先，利用直接積分的方法求出旋轉橢球體對內點的引力位，；然后，根據橢球體引力的特征以及麥克勞林定理，將內點引力位問題轉化成外點引力位問題[11]；最后將引力位在X、Y和Z方向上求一、二、三階導數，得到旋轉橢球體I在外部空間任意點P(x, y, z)所引起的重力異常、重力梯度異常以及重力垂直三階導數異常的基本表達式，如式(3)至(10)所示：

3 旋轉橢球體的重力異常場分布特征

由于旋轉橢球體與球體同屬于球狀體類型，因此無論是基本計算表達式還是重力異常場的分布特征都應該具有相似之處。因此將上述旋轉橢球體的推導結果與球體進行比較可以驗證推導結果的正確性。

3.1 基本表達式及比較

將兩種旋轉橢球體與球體的重力異常場基本表達式比較可知，基本表達式在結構形式上是相似的，都包含3個部分：即常數部分、質量部分以及與外部被吸引點和球心的坐標位置有關的部分。但各個結構組成部分的系數互不相同，由于旋轉橢球體自身結構還與其自身的長、短半軸尺寸密切相關，因此旋轉橢球體的基本表達式比球體更加復雜。

3.2 重力異常場分布特征比較

為便于結果比較，將如下具體尺寸代入上述公式進行計算。旋轉橢球體 I：a=321/4m，c=21/2m；旋轉橢球體II：a=21/2m，c=4 m；球體：R=2 m。該尺寸使得兩種旋轉橢球體與球體的體積和剩余質量相同。令深度h=4 m，剩余密度ρ=2.0×103kg/m3，得到3種球狀體的重力異常場在 XOY觀測平面上各個分量的平面等值線圖，如圖2所示。

圖2中重力異常采用SI制中重力的分數單位g.u.（1g.u.=10-6m/s2），重力梯度異常采用單位 E（1E=10-9/s2），重力垂直三階導數異常采用單位nMKS(1 nMKS=10-9/(ms2))。圖2顯示兩種旋轉橢球體與球體產生的重力異常場的分布特征相似，其中，由于旋轉橢球體I在地面上的投影為圓形，它自身的形狀特點決定了其重力異常 Vz、重力垂直梯度異常 Vzz和重力垂直三階導數異常 Vzzz的平面等值線與球體一樣，是一系列同心圓，而旋轉橢球體II則是一系列同心橢圓。3種球狀體的重力梯度Vxz分量在Y軸上的值均為零，Y軸左右兩端值的絕對值相等，符號相反。這一點從表 1中Vxz分量的極大值和極小值也可以得到驗證。此外，雖然剩余質量相同，但由于球體的質量較為集中，因此產生的各項極值的絕對值也最大。如表1所示，在上述條件下，Vz、Vzz、Vxz和Vzzz這4個分量的極大值中，球體的極大值分別是橢球體I的1.13、1.28、1.24和1.51倍，是橢球體II的1.22、1.46、1.52和1.79倍。

通過比較，由于旋轉橢球體與球體同屬于球狀體，因而在正演計算表達式和具體分布特征上都存在相同的變化規律，證明前面的推導結果是正確的。由于所選尺寸較小，數值差異不大，當場源體尺寸較大時：旋轉橢球體I有a=321/2m，c=2 m；旋轉橢球體II有a=2 m，c=16 m；球體有R=4 m。同樣在h=4 m的情況下，球體的Vz、Vzz和Vzzz這3個分量的極大值分別是橢球體I的2、3.5、7.56倍，是旋轉橢球體II的3、5和7.72倍，說明將等軸狀場源體簡單地近似成球體會造成一定程度的誤差。這是由于球體質量集中，而旋轉橢球體則隨著長、短軸的變化得到了較之分散的質量分布，使得三者在數值及細節特征上具有明顯差異。這種差異會隨著體積、質量、長短軸與球體半徑比例等因素的不同而不同。因此，用可以任意改變長、短軸比的旋轉橢球體作為等軸狀場源體的重力異常場模型是合理的。

圖2 三種球狀體的重力異常Vz、重力梯度異常Vzz和Vxz、重力垂直三階導數異常Vzzz的平面等值線圖Fig.2 Contour maps for gravity anomaly Vz, gravity gradient anomaly Vzzand Vxz, and the third derivative of vertical gravity anomaly of three kinds of spheroids: (a) Rotational ellipsoid I; (b) Rotational ellipsoid II; (c) Sphere

表1 重力異常場各分量極值比較表Tab.1 Extremum values of the gravity anomaly field’s components

4 結 論

提出了用旋轉橢球體作為等軸狀場源體的模型進行重力異常場正演計算的方法。通過對旋轉橢球體類型的劃分以及公式推導，完善了兩種類型的旋轉橢球體的重力異常場基本計算表達式。通過比較，在基本表達式的形式和分布特征上均具有與球體一致的規律，驗證了推導結果的正確性。同時，數值和局部細節特征體現的差異表明，在中心埋深淺、測量精度要求高的場合，將旋轉橢球體作為等軸狀場源體模型進行重力異常場分析是有必要的。同時，基本計算表達式的推導結果可為地下空穴和潛艇探測提供依據。

[1] Jekeli C, Abt T. Locating anomalies by gravity gradiometry using the matched filter with a probabilistic assessment[J]. Studia Geophysica et Geodaetica, 2011, 55(1): 1-19.

[2] Ardestani V E. Detecting, delineating and modeling the connected solution cavities in a dam site via microgravity data[J]. Acta Geod Geophys, 2013, 48: 123-138.

[3] Styles P, McGrath R, Thomas E, et al. The use of microgravity for cavity characterization in karstic terrains [J]. Journal of Engineering Geology and Hydrogeology, 2005, 38(3): 155-169.

[4] Jekeli C, Abt T. The statistical performance of the matched filter for anomaly detection using gravity gradients. OSU Report 494[R]. Columbus, Ohio, The Ohio State University, 2010.

[5] Dumrongchai P. Small anomalous mass detection from airborne gradiometry[D]. Columbus, Ohio: The Ohio State University, 2007.

[6] Leucci G, Giorgi L D. Microgravimetric and ground

penetrating radar geophysical methods to map the shallow karstic cavities network in a coastal area[J]. Exploration Geophysics, 2010, 41(2): 178-188.

[7] 孫嵐, 李厚樸, 邊少鋒, 等. 基于重力梯度的潛艇探測方法研究[J]. 海洋測繪, 2010, 30(2): 24-27. Sun Lan, Li Hou-pu, Bian Shao-feng, et al. Study on the submarine detection method based on the gravity gradient[J]. Hydrographic Surveying and Charting, 2010, 30(2): 24-27.

[8] Salah A M. Accurate and efficient regularized inversion approach for the interpretation of isolated gravity anomalies [J]. Pure and Applied Geophysics, 2014, 171(8): 1897-1937.

[9] 錢東, 劉繁明, 李艷, 等. 導航用重力梯度基準圖構建方法的比較研究[J]. 測繪學報, 2011, 40(6): 736-744. Qian Dong, Liu Fan-ming, Li Yan, et al. Comparison of gravity gradient reference map composition for navigation[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 736-744.

[10] 吳太旗, 歐陽永忠, 陸秀平, 等. 重力匹配導航的影響模式分析[J]. 中國慣性技學報, 2011, 19(5): 559-564. Wu Tai-qi1, Ouyang Yong-zhong, Lu Xiu-ping, et al. Analysis on effecting mode of several essential factors to gravity aided navigation[J]. Journal of Chinese Inertial Technology, 2011, 19(5): 559-564.

Modeling of isometric field-source in gravity anomaly field based on rotational ellipsoid

LIU Fan-ming1, LIU Hui-min1,2, JING Xin1
(1. College of Automation, Harbin Engineering University, Harbin 150001; 2. College of Mechanical and Electrical Engineering, Qingdao Agricultural University, Qingdao 266109, China)

The forward calculation of gravity anomalies would introduce larger errors when field-sources with finite sizes and isometric shapes are simplified as spheres. To solve this problem, a gravity anomaly model with isometric field-source body is proposed based on a rotating ellipsoid, which suggests that the spheres can be replaced by two kinds of rotational ellipsoids. The changing major and/or minor axis of the rotational ellipsoids disperse the mass’s compact distribution of the sphere, resulting in various different structures to match with the actual field-sources and reduce the errors produced in forward calculation. The complete calculation expressions about gravity anomaly are deduced and presented based on partial known formulas. Comparisons among the three spheroids are made to find their characteristic differences and to verify that the proposed method is correct and necessary.

rotational ellipsoid; gravity anomaly field; forward model; sphere; isometric field-source

P24

A

1005-6734(2016)01-0026-04

10.13695/j.cnki.12-1222/o3.2016.01.006

2015-10-14；

2016-01-20

國家自然科學基金（60834005）

劉繁明（1963—），男，教授，從事水下潛器定位技術、微弱信號測量與處理技術、無源導航定位技術。E-mail: Hrblfm407@hrbeu.edu.cn