基于旋轉(zhuǎn)橢球體的等軸狀場(chǎng)源體重力異常場(chǎng)模型構(gòu)建方法

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

（1.哈爾濱工程大學(xué) 自動(dòng)化學(xué)院，哈爾濱 150001；2.青島農(nóng)業(yè)大學(xué) 機(jī)電工程學(xué)院，青島 266109）

基于旋轉(zhuǎn)橢球體的等軸狀場(chǎng)源體重力異常場(chǎng)模型構(gòu)建方法

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

（1.哈爾濱工程大學(xué) 自動(dòng)化學(xué)院，哈爾濱 150001；2.青島農(nóng)業(yè)大學(xué) 機(jī)電工程學(xué)院，青島 266109）

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

旋轉(zhuǎn)橢球體；重力異常場(chǎng)；正演模型；球體；等軸狀場(chǎng)源

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

另一方面，隨著重力梯度測(cè)量技術(shù)的不斷提升，針對(duì)重力梯度測(cè)量不受外界信號(hào)干擾、不受測(cè)區(qū)環(huán)境影響、不向外界發(fā)射信號(hào)的無(wú)源特征[7-10]，孫嵐等提出了基于重力梯度探測(cè)的潛艇探測(cè)方法，將潛艇假設(shè)為一種繞長(zhǎng)軸旋轉(zhuǎn)的旋轉(zhuǎn)橢球體并進(jìn)行重力梯度正演計(jì)算[7]，然而由于只針對(duì)特定尺寸驗(yàn)證了該方法的可行性，并沒(méi)有給出具體的計(jì)算表達(dá)式，無(wú)法普遍應(yīng)用于其它尺寸的潛艇。推導(dǎo)繞長(zhǎng)軸旋轉(zhuǎn)型旋轉(zhuǎn)橢球體重力梯度異常場(chǎng)正演計(jì)算表達(dá)式可以為該方法的實(shí)現(xiàn)提供理論依據(jù)。

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

1 兩種類(lèi)型的旋轉(zhuǎn)橢球體

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

規(guī)定繞長(zhǎng)軸旋轉(zhuǎn)的旋轉(zhuǎn)橢球體為旋轉(zhuǎn)橢球體II，如圖1(b)所示，其球面方程為

當(dāng)a=c時(shí)，為球體。

圖1 兩種類(lèi)型的旋轉(zhuǎn)橢球體Fig.1 Two types of the rotational ellipsoid

2 旋轉(zhuǎn)橢球體重力異?；颈磉_(dá)式

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

3 旋轉(zhuǎn)橢球體的重力異常場(chǎng)分布特征

由于旋轉(zhuǎn)橢球體與球體同屬于球狀體類(lèi)型，因此無(wú)論是基本計(jì)算表達(dá)式還是重力異常場(chǎng)的分布特征都應(yīng)該具有相似之處。因此將上述旋轉(zhuǎn)橢球體的推導(dǎo)結(jié)果與球體進(jìn)行比較可以驗(yàn)證推導(dǎo)結(jié)果的正確性。

3.1 基本表達(dá)式及比較

將兩種旋轉(zhuǎn)橢球體與球體的重力異常場(chǎng)基本表達(dá)式比較可知，基本表達(dá)式在結(jié)構(gòu)形式上是相似的，都包含3個(gè)部分：即常數(shù)部分、質(zhì)量部分以及與外部被吸引點(diǎn)和球心的坐標(biāo)位置有關(guān)的部分。但各個(gè)結(jié)構(gòu)組成部分的系數(shù)互不相同，由于旋轉(zhuǎn)橢球體自身結(jié)構(gòu)還與其自身的長(zhǎng)、短半軸尺寸密切相關(guān)，因此旋轉(zhuǎn)橢球體的基本表達(dá)式比球體更加復(fù)雜。

3.2 重力異常場(chǎng)分布特征比較

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

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

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

圖2 三種球狀體的重力異常Vz、重力梯度異常Vzz和Vxz、重力垂直三階導(dǎo)數(shù)異常Vzzz的平面等值線(xiàn)圖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 重力異常場(chǎng)各分量極值比較表Tab.1 Extremum values of the gravity anomaly field’s components

4 結(jié) 論

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

[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] 孫嵐, 李厚樸, 邊少鋒, 等. 基于重力梯度的潛艇探測(cè)方法研究[J]. 海洋測(cè)繪, 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] 錢(qián)東, 劉繁明, 李艷, 等. 導(dǎo)航用重力梯度基準(zhǔn)圖構(gòu)建方法的比較研究[J]. 測(cè)繪學(xué)報(bào), 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] 吳太旗, 歐陽(yáng)永忠, 陸秀平, 等. 重力匹配導(dǎo)航的影響模式分析[J]. 中國(guó)慣性技學(xué)報(bào), 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

國(guó)家自然科學(xué)基金（60834005）

劉繁明（1963—），男，教授，從事水下潛器定位技術(shù)、微弱信號(hào)測(cè)量與處理技術(shù)、無(wú)源導(dǎo)航定位技術(shù)。E-mail: Hrblfm407@hrbeu.edu.cn