基于旋转椭球体的等轴状场源体重力异常场模型构建方法

刘繁明，刘惠敏,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