0 引言

密度分界面与区域地质构造、流体储层结构、不整合界面空间展布等有密切关系，利用重力异常估算密度分界面的起伏变化在油气资源普查、矿产资源勘探、深部构造研究、全球环境与气候变化研究中具有重要的意义(曾华霖，2005；屈念念等，2018；孙少波等，2018；冯旭亮等，2019).

经过国内外学者几十年的研究，逐步形成了较为成熟的物性界面深度反演方法体系，概括起来主要分为空间域和频率域两大类.其中空间域方法中比较典型的有线性回归法、压缩质面法和等效源法等(刘元龙和王谦身，1977；刘元龙等，1987；曾华霖，2005；冯娟等，2014a).线性回归法认为界面起伏平缓时，重力异常与界面的起伏近似线性关系，通过线性回归计算实现界面深度反演，是一种快速估计密度分界面深度的近似方法.压缩质面法和等效源法都是空间域迭代反演方法，通过不断的正演和修正界面模型，使计算的重力异常值与实际测量值之残差小到满足要求，从而得到较准确的密度分界面.近年来，为了提高反演密度界面的准确度，众多学者在反演算法中引入模型约束信息，包括模型几何形态约束、已知深度点约束等.在几何约束方面，可以采用光滑约束来提高反演结果的空间连续性，也可以采用全变差约束方式提高算法应对断裂等引起的界面形态突变的能力(冯旭亮等, 2014, 2018).在已知深度约束信息利用方面，通常可以利用约束信息建立初始模型，也可以在迭代过程中引入已知深度信息来控制反演迭代过程(林振民和阳明，1985；张盛和孟小红，2013；Chen et al., 2015).这两类约束中，已知深度约束较几何特征约束更符合实际地质情况，是最为重要的约束信息.频率域中的Parker-Oldenburg迭代反演方法具有计算快速的优势(Parker, 1973, 1974；Oldenburg，1974；关小平，1991)，因此被广泛采用.经典的Parker-Oldenburg方法层间物性只能为常数或横向二维变化，无法在深度方向变化.而实际的地层密度在水平方向，特别是深度方向是变化的，为了提高反演的效果，许多重力界面反演文献采用了变密度模型，例如线性变化模型、多项式模型、指数变化模型等(Barbosa et al., 1999；Garcia-Adbeslem，2000；冯锐，1986；Chai and Hinze, 1988；柴玉璞和贾继军，1990；柯小平等，2006；张盛和孟小红，2013；冯娟等，2014b；陈光希等，2019).

在约束密度界面深度反演中，已知深度信息多来自地震测深剖面解释结果.由于地震数据采集时间、数据精度、处理参数等存在差异，地震解释结果之间往往存在一定的系统性测线误差，层位之间很难形成闭合.在利用这些约束点信息进行约束反演时，如果系统性误差不消除，必定对反演结果产生影响，进而影响反演结果的准确评价.此外，在大面积重力数据解释时，所涉及的研究区往往跨越不同地质构造环境，不同地质构造环境下产生的地层密度结构差异较大，在利用重力资料进行某一密度界面深度反演过程中，不可能用同一密度模型参数来表示密度空间变化.以往为进行不同构造单元下物性界面反演，采用先分区单独反演再拼接的方式，这种方法不仅耗时耗力，而且最终结果的拼接很难进行，效果不够理想.为此，本文通过梳理变密度界面正反演算法，引入已知深度约束信息校验和水平密度分区，形成了一套约束变密度界面反演策略，并在南海研究区莫霍面深度反演应用中得到了验证.

1 反演算法 1.1 密度界面正演

密度界面正演算法是反演的基础，本研究采用速度快、精度高、且能够兼容多种纵向与横向变密度模型的频率域正演算法，以指数变密度模型为例(图 1)，假设指数变密度模型Δρ(z)=Δρ₀e^-μz，其重力异常频率域计算公式为(柯小平等，2006；张盛和孟小红，2013)：

(1)

式中Δρ₀(x, y)为密度差分布，μ为衰减系数，z为深度，F[Δg]为重力异常的频谱，G为万有引力常数，h(x, y)为界面起伏，z₀为界面平均深度，k为频率，其计算式为

(2)

其中δ_x和δ_y分别为x和y方向的网格间距，m_x和m_y分别为x和y方向的网格点距数，M和N分别为x和y方向的总点数.

1.2 迭代反演算法

在众多密度界面反演算法中，空间域迭代法具有比频率域迭代法更加稳定可靠等优势，为此本研究采用带约束信息的空间域迭代法进行密度界面反演，流程如下：

(1) 输入观测点和约束点数据：给出网格节点处界面引起的重力异常值Δg_obs，已知约束点平面坐标和深度值(x_j, y_j, z_j)，界面平均深度z₀，界面密度差水平分布Δρ₀，以及密度差随深度变化系数，根据约束点与重力异常的线性关系建立初始模型h⁰.

(2) 模型正演：根据公式(1)以及界面密度差Δρ₀、衰减系数μ和初始模型h⁰正演观测面重力异常Δg_cal，并计算其与观测值的差值.

(3) 计算模型修正量并更新反演模型：根据Δg_obs、Δg_cal和约束点信息，按照公式(3)修改模型参数：

(3)

其中，为无限平板修正公式，为深度调整系数，β取2，Δh^cons是根据约束点深度与重力异常按一定关系拟合得到的修正项(张盛和孟小红，2013)，λ为约束项权系数，在反演初期设置为1，随着迭代次数的增加逐步降为0.

(4) 模型正演：根据密度模型参数和界面起伏hⁿ⁺¹，正演计算观测面上的重力异常Δg_cal，并计算与观测值的差值.计算重力异常均方差rms=.

(5) 判断迭代终止条件：采用如下三个迭代终止条件判断迭代是否终止：(a)重力异常均方差rms小于给定误差限eps；(b)均方差rms不再减小；(c)迭代次数达到最大值n_max.满足其中任何一个条件迭代终止，输出反演结果.若上述三个条件均不满足，则转步骤2至5.

2 优化策略

为获得更加真实可靠的反演结果，根据实际需求，本研究在约束变密度反演算法之上增加了已知深度约束信息校验和水平密度分区两方面的应用策略.

2.1 约束信息校验

众所周知，约束信息的准确性对于约束反演结果及其评价影响巨大.在约束变密度界面反演中，约束信息通常来自不同时期获得的地震剖面解释结果，这些结果由于解释人员地质认识、处理参数和数据资料的差异，剖面之间往往存在不同程度的系统性误差，具体表现为在不同测线相交位置解释结果存在较大差异.如果没有在约束反演之前将这些约束信息的系统误差消除，必将影响最终反演结果及其评价的可信度.为此本研究提出在界面反演之前对约束信息的精度进行校验.

本研究采用经典的交叉线法进行约束信息内符合精度评估，并通过对约束信息按照测线进行调差，提高约束信息的准确度.与典型的航磁剖面调平不同的是这些约束测线在空间上没有明确的测量线和联络线，此时可以按照一定规则构建虚拟联络线(或趋势面)来进行调差(刘增群，1992)，具体计算过程如下：

(1) 利用资料，沿垂直测线方向构建虚拟联络线(或趋势面).

(2) 计算某条约束测线与联络线(或趋势面)的重合点并记录坐标.

(3) 利用插值计算测线和联络线(或趋势面)重合点处的值，根据差值均方误差作为精度评判结果，公式为

(4)

(4) 依测线差值平均值的半值对相交测线值进行调整.

(5) 如此反复进行，直至系统差小于给定值，或达到给定的最大调差次数，输出测线调差结果和最后估计精度.

2.2 密度分区策略

为应对横跨不同地质构造背景区域重力数据资料解释问题，避免单独分区反演产生的拼接问题，本研究引入横向密度分区策略.对于莫霍面深度反演，根据OBS等地震测深结果可以看出，陆壳一般存在较厚的沉积层和基底，壳幔密度差相对较大；而洋壳区域，沉积层和基岩通常较薄，地幔上涌，致使界面密度差相对较小，可见不同地壳类型的地层密度结构差异较大，有必要根据需求在反演时进行横向密度分区处理.

密度分区主要包括三个步骤：(1)根据研究区大地构造单元分区成果和地壳密度结构特征，按需求对研究区进行密度分区规划，形成分区图.(2)根据分区结果，对每一分区分别构建密度模型，之后拼接整合得到研究区分区密度差模型.(3)为消除分区密度边界的台阶状突变，有必要对分区密度模型进行一定程度的光滑滤波处理，从而避免由此引起的反演界面畸变.

3 模型实验 3.1 约束点校验与调差

以南海西北次海盆研究区OBH1996(Qiu et al., 2001)、OBS2006-1(吴振利等, 2011)、OBS2006-2(敖威等, 2012)三条莫霍面深度解释结果为例说明约束点校验(图 2).从三条解释剖面可以看出各测线交叉点处莫霍面深度解释结果均有较大差异，最大相差5 km，相对误差超过了30%.用交叉线方法对约束测线进行调差并估算其内符合精度.测线调差之后，测线交点处的均方差由2.62 km降到了1.56 km，最大绝对误差降到了3 km，约束信息内符合精度得到有效提高.

3.2 密度分区反演实验

为了说明密度分区在界面反演中的作用，我们设计了如下横向变密度界面模型，该模型包括两个起伏为5 km的界面凸起(图 3a)，分属于不同的密度分区内(图 3b)，取密度界面平均深度为10 km进行正演，得到该模型的重力异常(图 5c)，为了消除分区引起的边界畸变，我们对密度分区进行了高斯滤波得到光滑后的密度分区(图 3d).分别采用密度差为0.5 g·cm^-3、0.3 g·cm^-3和分区密度进行反演，得到图 3e、f、g、h的反演结果.

从反演结果可以看出，在未采用密度分区情况下进行反演，均无法同时反演出两个界面的真实起伏情况，当采用合理水平密度分区之后，反演结果与给定模型基本一致.采用光滑处理之后的密度分区进行反演，在不影响界面反演效果的同时，能够消除密度突变边界引起的界面阶梯状变形.模型实验表明，在研究区存在明显的横向密度差异时，有必要根据不同地质构造单元密度特征进行横向密度分区反演，才能够获得更为真实的反演结果.

4 南海莫霍面反演

南海是位于中国大陆以南，太平洋以西的半封闭海域，地形地貌复杂，海底平原，海底火山、海底高原、峡谷、海槽和海沟等均有发育，南海油气资源丰富，是世界上主要的产油区之一，对该区进行构造单元研究具有重要意义.我们将上述反演策略应用于南海莫霍面深度反演中，以期获得新的认识.

4.1 数据来源

本研究采用的重力数据来源于Getech公司提供的大陆边缘卫星测高重力数据，其精度约1 mGal，空间分辨率约2 km网度.使用美国国家海洋和大气管理局(NOAA)国家地球物理数据中心(NGDC)发布的ETOPO1数字高程模型(数据网度为1′×1′)进行了布格改正，得到研究区布格重力异常.为了提高重力异常分离可信度，本研究采用剥层法结合优化滤波方法进行异常分离，其中沉积层厚度数据来源于美国国家海洋和大气管理局(NOAA)发布的世界海洋和边缘海域的总沉积厚度数据，数据网度为5′×5′，最终获得了研究区莫霍面重力异常，如图 4所示.

4.2 莫霍面深度反演

(1) 约束点校验

本研究收集整理研究区范围内以往实施的OBS、OBH和声纳浮标等测深解释成果剖面12条(吴招才等, 2017)，共计168个已知深度点作为已知深度约束信息(图 5).为了提高约束信息的准确度，我们以CRUST1.0(Laske et al., 2013)提供的莫霍面深度构建趋势面为基准，对这些约束信息进行校验.调差之前，研究区约束测线与趋势面重合点处最大绝对误差达到12km，最大相对误差在50%以上，利用均方误差评价约束点内符合精度为3.7 km.进行调平处理之后，约束点内符合精度降为2.2 km，最大误差小于5 km.从深度与重力异常的关系(图 6)也可以看出，经过约束信息校验，数据点更加紧凑.

(2)密度模型分析

南海发育面积广泛，地质构造背景复杂，根据研究区已有密度资料，研究区主要存在海底面、基底面和莫霍面三个密度界面，其中海底面为海水层与海底沉积层的密度界面，密度差通常能达到1.4 g·cm^-3，通过完全布格校正予以消除.基底与沉积层的密度界面之间的密度差在0.2~0.4 g·cm^-3之间，与地层年龄与基底性质有关，采用剥层法去掉该部分异常.莫霍面是地壳与上地幔的分界面，重力异常的区域场随该界面的起伏而发生变化.横向上，研究区大致可以分为洋壳、过渡壳和陆壳三种地壳类型，不同地区的同时代地层密度差异较大，存在明显的横向密度分区特征(图 7).

参照广州海洋地质调查局1:200万地质地球物理图集中的南海地壳速度结构图和大地构造格架图(陈洁和温宁，2010)，将研究区分为洋壳、陆壳、过渡壳三种地壳类型.根据收集的地震剖面解释结果，采用波速密度转换经验公式，估计了研究区沉积层、基底层和上地壳密度，对分别建立密度差模型，经平滑处理得到分区密度差模型(图 8)，衰减系数0.002.

4.3 反演结果及评价

采用校验之后的约束信息以及分区密度差，设置平均深度18 km，经过12次迭代反演得到莫霍面深度，最终均方误差为1.7 mGal.图 9为反演结果与约束点剖面对比，可以看出，尽管约束剖面局部细节与反演结果局部存在一定差异，但是整体上界面反演结果与地震测深结果吻合较好，地震测深与反演结果相关系数为0.95，平均深度差值为0.3 km，反演结果准确可靠.

从反演结果可知(图 10)，南海研究区莫霍面界面深度为9 km至32 km，基本反映了陆域到海域地壳从陆壳、过渡壳、洋壳的过渡特征，整体上与已有结果保持一致(秦静欣等，2011).在莺歌海盆地、曾母盆地等油气资源产区莫霍面深度较浅，反映了盆地下部莫霍面隆起，上地幔上涌，为盆地提供了热源，具有很好的生烃条件，与该盆地地区广泛发育的油气资源相关.

5 结论

本文通过梳理变密度界面正反演算法，提出了在已知深度信息约束下的变密度界面反演策略，该策略主要包括对深度约束点的校验、水平变密度分区模型和约束变密度反演算法三个方面.其中已知约束点校验能够有效消除约束信息之间的系统误差，提高约束信息的可信度；在反演过程中引入水平密度分区以应对不同构造背景密度界面模型，无需后期反演结果拼接.最后将本策略应用于南海莫霍面深度反演研究中，结果显示该策略能够在已知信息约束下，获得较为可信的深度反演结果，验证了方法的可行性.

致谢  感谢审稿专家提出的建设性意见，感谢德克萨斯大学达拉斯分校任丽博士对本文提出的改进建议.