作为地球物理学研究中的一个重要领域，重力学是研究地球本体，特别是地壳和上地幔结构与物性特征的主要方法之一.它主要根据在地表观测的重力异常，通过正演、反演方法来研究地下物质的结构和构造状态(王谦身等，2003；滕吉文等，2004).重力异常是不同深度物质产生重力场的叠加效果，如何从重力异常中提取必要的信息，并反演出地球内部的密度分布，是重力学的重要研究目标.重力学研究内容主要包括两个方面：

(1)确定剩余重力异常，即从布格重力异常中划分出地质体密度变化或岩性不同而产生的异常.

(2)确定引起剩余重力异常的地质体的几何参数和密度(侯遵泽和杨文采，1997).20世纪80年代以来，我国学者在重力学方面做出了众多成果.王懋基等(1997)展望了重力勘探的发展.管志宁等(2002)讨论了21世纪重磁位场的发展方向.刘元龙等(1987)给出了利用三维重力方法计算地壳构造的方法，计算了我国一些区域的地壳厚度及莫霍界面分布.方剑(1999)利用卫星重力资料反演了全球的地壳厚度和岩石圈厚度，得出了全球地壳厚度图.

野外观测得到的重力异常，不仅包含高频的噪声及观测误差影响，更是从浅到深多种地质体密度差异引起的叠加异常(王谦身等，2003).因此，在应用重力数据推断地下物性结构时，往往要对重力异常进行分离，提取不同尺度研究对象的异常.分离重力异常的主要方法是滤波，包括空间域滤波(王宜昌和杨辉，1986;Pawlowski，1995)和频率域滤波(Pawlowski and Hansen，1990;杨文采等，2001).在得到研究对象所引起的重力异常后，通常用正演拟合和反演方法来获得地下的密度结构.

正演拟合方法是通过实测资料建立初始模型，并计算其异常，根据理论异常与实测异常的差异不断调整初始模型，使理论异常和实测异常能够很好的拟合，从而获得最终的密度模型.随着计算机技术的不断发展，即使对于很复杂的密度模型也能很快地计算出理论异常.更重要的是，计算机不仅能进行理论计算，而且能够实时的显示出物性模型和位场异常的计算结果，及时给解释者提供大量直观信息.通过图像化的人机交互界面，解释者可以方便的修改物性模型.交互模拟解释系统给解释者提供了很大的便利，提高了解释者的工作效率(楼海，2001).

密度反演方法是通过建立重力异常与密度模型的对应关系，并构建偏导数矩阵，通过对比模型的正演异常和实测异常，使用一定的算法寻找新的模型参数估计，缩小拟合误差，对密度模型修正并多次迭代，从而获得最终的密度模型.密度反演方法计算速度比较快，但是，在某些情况下，模型与异常之间是非线性关系，最优解不是唯一的；另外，数学意义上的最优解也不一定符合实际地质情况(楼海，2001).

本文重点介绍重力异常的分离方法，以及利用重力异常通过正演拟合和反演手段构建密度模型的方法.

实测重力异常数据是叠加异常，包含了从地表到深部所有密度不均匀引起的重力效应.引起重力异常的主要地质因素有：地球深部莫霍面的起伏;地壳、岩石圈地幔物质密度的横向变化等；结晶基岩内部的密度变化；结晶基底顶面的起伏；沉积岩的构造和成分变化和其他密度不均匀因素如高密度金属矿床等(王谦身等，2003;曾华霖，2005).重力解释通常把实测重力异常看作由区域异常和局部异常组成，区域异常是由分布范围广、相对深的地质因素引起的，局部异常是由比区域范围小的研究对象引起的，二者的概念是相对的，与研究对象的尺度范围有关(王谦身等，2003;曾华霖，2005).

大部分异常分离方法的前提是不同异常在空间频率有差异，一般来说，由深部场源引起的宽缓异常具有低频的性质；而由浅部场源引起的局部异常具有高频特征；观测误差和近地表密度不均匀引起的重力效应则表现为高频干扰(王谦身等，2003;曾华霖，2005).因此，重力异常主要采用滤波的方法分离，包括平滑、高次导数和解析延拓等空间域滤波(王宜昌和杨辉，1986;Pawlowski，1995)以及傅里叶变换、小波分解等频率域滤波(Pawlowski and Hansen，1990;杨文采等，2001).将空间域位场数据变换到频率域进行数据处理，可以充分发挥位场数据的横向分辨率优势，突出位场异常横向变化的整体规律和性质.利用频率域的特殊性质，可以将位场的计算公式在频率域中大幅度的简化，有效地提高计算速度和效率(陈石和张健，2006).

局部异常与区域异常具有相对性，因此，随着研究范围的不同，研究对象引起的局部异常在布格异常中的比例也会发生变化，重力异常的分离在不同的研究范围中有不同的应用.

在重力勘探中，研究对象主要是矿体、油气藏等异常体，引起重力异常的主要地质因素如莫霍面的起伏、结晶基岩内部的密度变化、结晶基底顶面的起伏等不在研究范围内；对盆地进行研究时，深度可到结晶基底，但仍有深部因素未被考虑.在上述情况下，研究对象引起的局部异常在总异常中的比例较小，因此滤波用于压制其他地质因素产生的异常，突出浅层研究对象的异常(刘清泉，1986;李盛汉，1997;王勤等，2004;刘彦等，2012).

在壳幔密度结构研究中，研究对象是地壳和上地幔的密度结构.绝大部分引起重力异常的主要地质因素在此研究范围内.在此范围以下，地幔和地核形状规则，密度变化平缓(曾华霖，2005)，加之引力与距离的平方成反比，因此其产生的异常小且变化缓慢.王谦身和杨新社(1997)通过中国大陆三维地震P波速度分层模型建立了7个密度分层模型，并计算了各层的重力异常，他们指出：由于“浅层烙印效应”的存在，综合重力异常主要反映浅层部分特征.可见，在壳幔结构研究中，局部异常在总异常中所占比例很大，其分布和数值都与总异常类似(王谦身和杨新社，1997)，滤波主要用于消除测量误差和近地表密度不均匀等影响(贺日政等，2001).

采用经各项改正后的重力异常资料，研究地壳、上地幔的密度结构，是重力学的重要研究对象(王谦身等，2003).利用重力异常构建密度结构分为正演拟合和反演两类方法.由密度结构计算重力异常，是正演拟合和反演方法的基础工作.因此，我们首先介绍计算重力异常的基本方法，进而对正演拟合和反演构建密度结构的方法进行说明.

根据模型的不同，正演计算分为三维模型和2.5维模型两种情况.

三维模型中各单元的密度和形状在三个维度上均可发生变化，可以很好的反应真实的密度分布，直角坐标系内重力异常正演计算的基本公式为

其中V为重力位，G为万有引力常数，p=(x，y，z)，q=(ξ，η，ζ)分别为观测点和异常体内点的坐标；r_pq为p与q点的距离，表达式为r_pq=[(ξ-x)²+(η-y)²+(ζ-z)^2]］^1/2，ρ为异常体内密度函数.

三维模型可以较好的反映实际密度分布，可以用于模拟面状的重力异常数据，但是计算量较大.在实际工作中，测量面状重力异常的工作量较大，一般会采用布格异常图、卫星重力异常等数据。实测重力数据多为线状，一般采用2.5维模型进行研究，即各单元的密度和形状在两个维度上可发生变化，而在另一维度上是恒定的，可能为一维尺度远大于另二维尺度的柱状体或无限延伸的二度体，本文中我们只讨论二度体的情况.

在具有均一密度的情况下，设二度体的走向与y轴平行，则直角坐标系内二度体重力异常正演计算的基本公式为

其中V为重力位，G为万有引力常数，p=(x，z)为观测点的坐标；r_pq为p与q点的距离，表达式为r_pq=[(ξ-x)²+(ζ-z)²]^1/2，ρ为异常体密度.

对于具有任意形状的异常体q而言，其体积分难以计算，一般选择将其近似为规则形体进行计算，根据近似方式的不同，可以分为有限单元法与边界单元法两类.

有限单元法的基本思想是将复杂形体进行分割，使之转化为一系列简单形体的组合，计算简单形体的重力异常再求和即可得整个复杂形体的异常，常用的有点元法，线元法和面元法等.

点元法是指将任意形体用三组平行于直角坐标面的截平面进行分割，使物体分成许多具有规则形状的长方体元.用解析方法计算出所有这些长方体在计算点产生的异常，并累加求和，就可以得到整个形体在计算点引起的异常值.点元法不适用于2.5维模型.

线元法是指用两组相互垂直的平行面，将复杂形体分割成一系列棱柱体，每个棱柱体以其中的线元来表示，用解析法求解线元方向的一重积分，然后在垂直线元的x，y 方向分别作数值积分，即可得出整个地质体的近似异常值.

面元法是以一组互相平行的铅垂面将复杂形体分成若干个直立薄片，每一片又用适当的多边形来逼近其形状，用解析法计算每一薄片在计算点的作用值，最后对所有薄片作用值进行数值积分，即可得到整个复杂形体引起的异常.

在上述三种方法中，点元法除适用于任意测网外，还可以模拟非均匀物性的情况，但是其计算量比较大.线元法中线元的划分有较多的限制，只能模拟密度在两个方向上的变化，精度不如点元法，但计算量比点元法小.面元法计算量最小，但是只能模拟密度在一个方向上的变化，因此精度最低(李焓等，2008).

边界单元法的基本思想是将复杂形体异常的体积分通过奥-高公式转化为面积分，再由格林公式转为线积分，而后累加求和得到整个形体的场值.常用的有多边形面多面体方法和三角形面多面体方法(李焓等，2008).

多边形面多面体方法是用一系列多边形来逼近复杂形体的形状，只要给出角点坐标，就可以计算出该形体的重力异常值.

三角形面多面体方法是任意形状和大小的三角形构成的多面体来逼近任意形体，只要知道多面体各角点的坐标，就能方便地用解析公式计算出多面体引起的重力异常.

有限元法和边界元法都是用规则形体的重力异常来表示任意形体的重力异常，前人给出了上述方法中涉及到的规则形体重力异常的计算公式。较简单的形体，如直立长方体、水平二度体等的计算公式，王谦身等(2003)和曾华霖(2005)给出了详细的计算过程；对于较复杂形体的计算公式，何昌礼和钟本善(1988)对三角形面多面体的计算公式，Okabe(1979)对多边形面多面体的计算公式分别进行了推导，本文不再一一列出.

正演拟合构建密度结构一般采用人机交互模拟来实现，其步骤可归纳为：

(1)给出地质体初始模型；

(2)输入给定的物性参数；

(3)初始模型理论异常的计算；

(4)对比理论与实测异常；

(5)根据对比结果修改初始模型，再重复上述步骤，直至拟合满足解释的要求为止(郑洪伟等，2006).

很多学者使用人机交互模拟对地壳密度结构进行了研究(施小斌等，2002;刘占坡等，2003;彭聪，2005;郑洪伟和贺日政，2006;唐新功等，2006，2008，2012;张季生等，2008;Sanchez et al.，2010;刘一峰等，2012;王谦身等，2013;Zhang et al.，2014;王芃等，2014).

人机交互模拟使用边界元法建立模型，单元的形状代表异常体的形状(Sanchez et al.，2011)，因此人机交互重力模拟需要较强的先验信息约束.

IGMAS(Interactive Gravity and Magnetic Application System)是常用的人机交互重力模拟软件之一，用于三维与2.5维重力、电磁的交互式正演模拟.该软件的优点是图形界面，可方便地进行人机对话，模型修改非常方便，实时显示模拟结果，并可加入钻井等已知的地质及地球物理信息(柯小平等，2009).

IGMAS的三维模型由一系列互相平行的剖面组成，每个剖面上有若干多边形，代表异常体在该剖面上的截面.同一异常体在不同剖面上的截面通过三角形面连接，组成异常体的三维形态(图 1).通过三角形面多面体的重力异常计算公式(Götze and Lahmeyer，1988)，计算模型在观测面产生的重力异常，通过人机交互重力模拟得到密度结构.

图 1

Fig. 1

图 1 IGMAS 3维模型示意图 Fig. 1 Sketch map of 3-D modeling in IGMAS

IGMAS的2.5维模型仅有一条剖面，剖面由若干多边形组成，代表水平二度体在剖面上的截面.通过Won和Bevis的水平二度体的重力异常计算公式(1987)，计算模型在测线上产生的重力异常，通过人机交互重力模拟得到密度结构.

Sanchez等人利用IGMAS对加勒比板块和南美板块边界的岩石圈结构进行了研究(2011)(图 2).文中以宽角剖面、接收函数和地表地质等作为约束条件，建立了初始模型的几何结构；通过速度-密度经验关系式和岩石的化学组成得到初始密度值.最终的密度结构显示委内瑞拉东西部具有不同的结构特征：委内瑞拉西部有较长的加勒比板块出现，板块向南延伸至哥伦比亚北部(图 3a)，板块的倾角在100 km深度处由15°增加到32°；而在委内瑞拉东部，加勒比板块较短(图 3b)；从密度上看，南美北部的下地壳较轻，而加勒比板块的地壳虽然有异常厚度，但具有典型的洋壳密度.

图 2

Fig. 2

图 2 (a)南美北部边缘构造纲要图；(b)图 2a中的红色虚线区域的简化地质图，其中红色实线为四条剖面的位置(Sanchez et al.，2011) Fig. 2 (a)General tectonic of the northern margin of South America;(b)Simplified geologic map(red stippled line in Fig. 2a). Red lines show the position of main sections in the study area(Sanchez et al.，2011)

图 3

Fig. 3

图 3 委内瑞拉的岩石圈结构(Sanchez et al.，2011) (a)70°W剖面的密度结构与拟合情况；(b)64°W剖面的密度结构与拟合情况，剖面位置见图 2b. Fig. 3 Lithospheric structure in Venezuela(Sanchez et al.，2011) (a)Density profile in 70°W;(b)Density profile in 64°W，position of the profiles can be referred in Fig. 2b.

常用的重力反演主要包括确定物性界面深度及起伏和确定地质模型参数两种(王谦身等，2003;曾华霖，2005).前者一般用于计算基底界面或莫霍面起伏等情况，如Oldenburg根据Parker公式，提出了在已知密度差的条件下计算密度界面形状的方法(Oldenburg，1974)；后者则可以确定密度结构，是利用重力数据推断地下结构的常用方法.

为了适合于起伏地形、测网不规则、密度非均匀分布等复杂条件和避开异常源形体选择，研究者们发展了一种“密度线性反演”方法(王谦身等，2003).这种方法使用有限元法建立模型，如三维点元模型(王新胜等，2012)或2.5维线元模型(于鹏等，2007)，模型边界与异常体边界无关(Li et al.，1998).构建关于拟合误差的目标函数，使用一定算法求使目标函数取极小值的密度分布.方剑和许厚泽(1999)应用阻尼最小二乘反演得到了中国及邻区岩石层密度结构；杨辉(1998)应用奇异值分解(SVD)算法得到了Y盆地的基底密度分布；王新胜等(2012)应用代数重构技术得到了华北地区的密度结构；于鹏等(2007)应用快速模拟退火算法(VFSA)进行了模型实验并对徐闻地区的盖层结进行了反演.Li等(1998)提供了一种针对三维点元模型的反演方法，可以得到同时满足重力异常拟合误差和限制条件的平滑模型.这种方法因为可以加入限制条件与深度补偿而得到广泛的应用，张毅等(2012)应用这一方法研究了三峡地区中上地壳的密度结构.

GRAV3D是英属哥伦比亚大学地球物理反演中心(UBC-GIF)开发的三维重力正反演计算系统，反演中使用Li等(1998)的方法(http://www.eos.ubc.ca/research/ubcgif/iag/sftwrdocs/grav3d/background.htm).

Li等(1998)建立如下误差函数为

其中d ^obs=(g₁，g₂…g_n)^T为实测异常向量，d为理论异常向量，W_d为对角矩阵{1/σ₁，1/σ₂…，1/σ_n}，σ_i是第i个观测数据的标准差.密度分布首先要保证φ _d 足够小，但符合这一条件的密度分布很多，因此还需加入通过根据先验信息建立的参考模型，使密度分布尽量接近参考模型，并使其在三个方向平滑变化，因此建立如下目标函数为

w_s，w_x，w_y和w_z为权重函数，而a_s，a_x，a_y和a_z为影响不同分量在目标函数中相对重要性的系数，w(z)为深度权重函数，ρ₀为参考模型.权重函数可以加强或抑制某一分量的影响，参考模型可以代表先验信息，这些因素使得反演可以在限制条件下进行.

上述函数随深度衰减，会导致反演的密度分布集中于浅部，而深部的密度经常为0，为消除这一影响，Li等(1998)引入了深度权重函数为

其中β一般取2，而z₀与单元尺度有关.深度权重函数使深部和浅部的单元获得非0值的机会相同.

使用有限差分逼近，式(4)可以写为

其中ρ和ρ₀为长度m的向量，模型权重矩阵包含了式(4)中的权重函数与系数.

为解决反演问题，需要寻找一个密度模型ρ，使φ_m有最小值且数据误差在要求范围内，可以通过寻找函数φ(ρ)=φ_m+λ^-1(φ_d-φ^*_d)的最小值得到，式中φ^*_d是要求的精度范围，λ是拉格朗日乘子.最小值解通过广义子空间法迭代求出，避免了大型矩阵运算.如果想获得非负或非正密度模型，可以将正负性加入子空间算法中.

Welford等(2010)利用GRAV3D软件对爱尔兰大西洋型大陆边缘的岩石圈密度结构进行了研究(图 4).文中使用DNSC08重力异常数据，以海洋测量、沉积层厚度和地震资料作为限制建立层状参考模型，反演得到了区域密度结构.反演结果显示了与地震资料一致的莫霍面分布(图 5)；一些地区地壳中上地壳部分减薄主要由下地壳组成，这些地区的位置与蛇纹石化地幔上涌位置有关，因此缺少深源地震；通过沉积层厚度与均衡地壳理论厚度的对比，发现厚度有差异的地区可能继承了加里东基底的结构特征.反演结果提供了对爱尔兰大西洋型大陆边缘的独特理解，是前人地球物理研究的补充，可以为重建北大西洋裂谷演化史提供进一步的约束，并指导将来的地球物理研究工作.

图 4

Fig. 4

图 4 北大西洋地区等深图(Welford et al.，2010)，黑色线框范围是爱尔兰大西洋型大陆边缘，即为三维重力反演范围 Fig. 4 Bathymetric map of the North Atlantic region(Welford et al.，2010)，the black box outlines the Irish Atlantic continental margin，that is the area for 3-D gravity inversion

图 5

Fig. 5

图 5 爱尔兰大西洋型大陆边缘岩石圈密度结构剖面与莫霍面形态的栅状图(Welford et al.，2010) Fig. 5 Fence diagram of arbitrary slices in the density structure and the topography of Moho of the Irish Atlantic continental margin(Welford et al.，2010)

3.1     重力异常的分离是利用重力资料构建密度结构的前提，主要是通过空间域和频率域滤波，不同的研究范围有不同的重力异常分离要求，在构建壳幔密度结构时，主要用于消除测量误差和近地表密度不均匀等影响.

3.2    正演拟合和反演都包含重力异常的计算，任意形状异常体建模困难，且重力异常难以计算，通常采用规则形体的异常来近似，常用的建模方法包括有限元法和边界元法等.

3.3     以人机交互重力模拟方法为代表的正演拟合使用边界元法建立初始模型，通过不断调整初始模型计算重力异常，最终获得与观测异常拟合较好的模型.正演拟合可以得到精细的壳幔间断面等复杂结构，拟合过程中有大量的人为参与，受人为影响较大，为了得到合理的模型，需要很强的先验信息约束.

3.4     以线性反演方法为代表的重力反演使用有限元法建立初始模型，建立关于拟合误差和初始模型的目标函数，使用一定的算法获得使目标函数取极小值的密度分布.重力反演计算速度快，反演过程无人为参与.得到的密度结构具有平滑的特征，对精细结构反映较差.重力反演的非唯一性较强，反演结果强烈依赖初始模型.