重力异常是地下各岩性不均匀变化与各物性界面起伏的综合反映，重力密度界面反演对于研究区域构造基底起伏、莫霍面起伏以及岩石圈厚度具有非常重要的意义，一直以来都是地球物理反演中的重要问题(Telford et al., 1990；Kaufman et al., 1992；曾华霖等，2005).反演过程中不同的基准面深度与界面密度差会反演出不同的结果，因此确定界面基准面以及界面密度差具有十分重要的意义.

近些年来，众多地球物理学家对如何选取界面密度差和基准面深度展开了研究，并且取得了许多研究成果.方剑(1999)使用统计模型，对洋中脊、海洋和大陆的地壳及上地幔密度分别赋值进行反演；冯锐(1985)和方剑等(1999)在计算沉积层与地壳平均密度时，除了参照各地区在做重力反演时曾选用过的数据外，还根据地震纵波的实测值使用层速度-密度转换公式得到界面密度差；柯小平等(2006)对青藏高原地区的莫霍面结构采用了变密度反演，浅层的密度由岩芯及野外标本实测得到，中深层的密度可由地震波速-密度的经验公式转换得到.然而不同的区域地质情况相差很大，对于不同的区域使用统计模型中统一的密度差是不合适的；此外在野外工作较少的地区，浅层岩芯及标本和中深层的地震资料较少，速度-密度转换公式也只是经验公式，在实际应用中会带来误差.求取重力基准面的传统方法是根据重力异常的功率谱估计界面的平均深度值，然而在功率谱的推导过程中有很多的假设和近似，在最后根据斜率求取也是根据个人经验，因此功率谱只能够大致估计界面的埋深，导致反演结果与已知深度点无法拟合的很好(穆石敏等，1982；侯重初等，1988；李成立等，1998；徐科军等，2003；闫建平等，2008；戴伟铭，2010)；或者是在反演的区域选取不同参考深度进行运算，对每一种反演结果与已知深度点进行比较，选取误差较小的作为最终结果，这种算法计算量大以及所选取的不同参考深度数量有限并且在很大程度上依赖于个人经验，这无疑造成了重力反演结果的不准确性.

上述方法没有或者很少使用控制点深度约束反演过程，使最终反演结果可能与控制点相差较大.因此，本文采用一种密度界面上已知若干控制点深度的迭代反演法，通过控制点来计算最合适的界面密度差和基准面深度，逐次迭代使最终得到的界面反演深度与已知控制点相差最小，实现了三维单一密度界面起伏的模型计算，并对南海地区的实测剖面进行了反演，取得了较为理想的结果.

1 三维单一密度界面的正演计算

正演计算是反演过程的重要组成部分，概括起来主要分为两大类，即频率域和空间域.频率域主要是基于快速傅里叶变换的Parker反演法，计算速度很快.在空间域中对于形状不规则的地下介质，可将其分割为不同的几何单元分别求取异常，再相加得到总异常.根据分割单元方式，正演方法可以分为线元法、面元法、体元法、表面积法及三角棱柱组合法等，较为常用的为长方体组合模型.

本文采用长方体组合模型，如图 1所示，在坐标原点处产生的重力异常为

图 1

Fig. 1

图 1 小立方体模型 Fig. 1 Single prism model

现在的界面反演方法大都是在频率域中利用Parker-Oldenburg公式迭代反演，其优点是速度快，但存在对重力的“向下延拓”过程而产生高频噪声，为了消除高频噪声，需要采取低通滤波，在滤波的同时也把数据中有用的高频数据滤掉了(Oldenburg et al., 1974；柯小平等，2005).因此本方法结合压缩质面法与等效源法(刘元龙等，1977；冯娟等，2014)，在空间域中进行反演运算，使运算过程中的误差尽量减小.

如图 2所示，h为地表，m为重力基准面，起伏的界面可以划分成许多小立方体.在反演过程中将小立方体压缩到压缩面上，根据重力异常值反演其高度，使用控制点信息调整界面密度差和基准面深度，计算出界面起伏深度，再对其正演值与实际重力异常值的偏差再次反演，直到到达一定次数或者偏差值达到指定范围为止，具体流程如下.

图 2

Fig. 2

图 2 地下小立方体的划分 Fig. 2 Subsurface prism

如图 3所示，假定已知反演区域内m行n列的布格重力异常值G(i，j)，(i=1，2，…m;j=1，2，…，n)，以及反演深度内已知深度的e个控制点深度h(i)，(i=1，2，…，e)，其对应的重力异常值为Gn(i)，(i=1，2，…，e)，将地下网格也划分成m行n列，坐标与重力异常值坐标一一对应.

图 3

Fig. 3

图 3 重力异常值 Fig. 3 The gravity anomaly

(1)初始基准面及界面密度差的建立

设初始基准面的深度为m₀，初始密度差为p₀，根据无限水平板公式：

2πFp₀(h(i)－m₀)=Gn(i)(i=1，2，3，…，e)

F为万有引力常数，根据这e个方程，使用最小二乘法求取m₀，p₀.

(2)求取初始面密度矩阵 σ ₁

第一次反演将m₀作为压缩面，将代表起伏界面的小立方体压缩到m₀上，使其成为面密度不均匀的物质面.每一块的面密度σ₁(i，j)可看为(depth₁(i，j)－m₀)·p，depth₁(i，j)为点(i，j)处起伏界面深度，则depth₁(i，j)－m₀为每个小立方体的高，也就是相对于压缩面的起伏，p为界面密度差.

被压缩到m₀面上的小长方形(i，j)在重力观测点(k，l)所产生的重力异常值Δg为

对于观测面每一个重力异常值G(k，l)，可得

对于地面上所有重力异常值来说，可写出下面的线性方程组：

(3)求取初始模型 depth ₁

由于初始的基准面和密度差是由控制点使用无限水平板公式求得的，当界面起伏较大时会有较大误差，因此，需要对基准面和界面密度差进行改进.

设改进后的基准面为m₁，密度差为p₁，使用e个控制点所对应的的面密度Ω₁(1)，Ω₁(2)，…，Ω₁(e)和控制点深度建立方程组：

由最小二乘法求出m₁和p₁，则起伏界面深度depth₁(i，j)可由各个小立方体的高度和基准面深度求得：depth₁(i，j)=m₁+σ₁(i，j)/p₁，求得初始模型 depth₁.

(4)迭代

迭代的目的就是改变基准面深度和界面密度差，使结果更加逼近真实情况.使用小立方体正演depth ₁得到重力异常值 g ₁，将实际观测异常 G减去g₁作为重力异常值G₂，由于此时的反演是在 depth₁上的改进，因此以depth ₁为压缩界面，m₁为基准面，p₁为密度差，重复上述计算得到面密度 σ ₂.

e个控制点对应的的面密度Ω₂(1)，Ω₂(2)，…，Ω₂(e)，此时计算新的基准面深度时所用的面密度 Ω是第一次的面密度Ω₁加上改进的面密度Ω ₂，即Ω(i)=Ω₁(i)+Ω₂(i)，(i=1，2，…，e)，设改进后的基准面深度为m₂，界面密度差为p₂，则有：

通过最小二乘法，得到m₂与p₂，以m₂为新的基准面深度，p₂为新的界面密度差，σ=σ₁+σ₂为面密度，求得第二次的界面深度模型 depth ₂.

(5)继续迭代，直到到达一定的迭代次数或者反演结果符合以下两个标准

标准1： 反演结果的正演重力异常值与原始重力异常值的偏差达到指定范围.

标准2： 控制点与反演结果的差值达到指定范围.

3 模型试验

下面以三维密度界面上存在若干控制点为例进行反演实验.图 4为合成的起伏界面深度模型，基准面为30 km，密度差为0.5 g/cm³，黑点为控制点位置，图 5为按照上述方法迭代五次之后反演得到的界面起伏，反演结果与模型的平均偏差为2×10^-6 km，标准差为0.005 km.迭代过程中每次得到的基准面深度和密度差如图 6、图 7所示，可看出第一次使用无限板公式获取基准面深度和密度差时，所得基准面和密度差与真实值相差较大，但是在此后的迭代中会逐步逼近真实值，反演结果与模型吻合很好.

图 4

Fig. 4

图 4 深度模型(黑点为控制点位置，单位：km) Fig. 4 Moho model(black dots are control points，unit:km)

图 5

Fig. 5

图 5 反演结果(单位：km) Fig. 5 Inversion result(unit:km)

图 6

Fig. 6

图 6 每次迭代得到的基准面深度 Fig. 6 Reference plane in each iteration

图 7

Fig. 7

图 7 每次迭代得到的密度差 Fig. 7 Density contrast in each iteration

南海范围广阔，构造演化复杂，其莫霍面可从20多公里变化到中央海盆的10 km左右，起伏较大，前人已做了大量工作(闫贫等，2002；Wang et al., 2006；阮爱国等，2011；Qiu et al., 2011；李家彪等；2011；卫小东等，2012；丘学林等，2012).本次选以南海(图 8)为例进行带控制点的三维密度界面反演，将地震结果作为控制点，控制点点距20 km，共594个，位置如图 9所示.反演结果如图 10所示，得出的莫霍基准面深度为23.625 km，密度差为0.5402 g/cm³，控制点反演结果与地震结果的平均偏差为-2.1×10^-14 km，标准差为1.56 km，对反演结果使用小立方体正演得到的重力值与原始重力异常值的平均偏差为-0.0123 mGal，标准差为0.0284 mGal.同时使用Parker法对不同平均深度和密度差进行多次实验，选出最好的反演结果(图 11)，其平均偏差为2.23 km，标准偏差为2.56 km.此外，如图 12所示，抽取不同位置的两条地震剖面比较两种反演结果，可看出带控制点的三维界面反演方法反演结果优于传统Parker法.由偏差分析和剖面比较可知，带控制点的三维界面反演方法相比于传统的Parker法与控制点拟合的更好，结果更接近于真实值，而且相比于Parker法的多次反演选最好结果，带控制点的三维界面反演方法也更加简便.

图 8

Fig. 8

图 8 研究区范围 Fig. 8 The research aera

图 9

Fig. 9

图 9 南海重力异常 Fig. 9 The gravity anomaly of south China sea

图 10

Fig. 10

图 10 带控制点的三维界面反演结果 Fig. 10 The result of three-dimensional interface inversion with control points

图 11

Fig. 11

图 11 Parker反演结果 Fig. 11 The Parker inversion result

图 12

Fig. 12

图 12 两种反演结果的比较 Fig. 12 The coparision of two inversion results

本文提出的带控制点的三维密度界面反演方法使用小立方体模型，避免了传统的Parker法向下延拓时产生的误差，在每次迭代过程中使用控制点约束基准面深度和密度差，使反演结果与控制点拟合较好.从本文的模型计算和实例反演来看，已知控制点数据的控制约束作用比较明显，从而避免了用不同参数分别计算，再将不同反演结果与控制点比较来选取最优结果，减少了多解性，提高了反演效率.