南海是西太平洋的边缘海盆之一，地处欧亚板块、菲律宾海板块和印度-澳大利亚板块的交汇地带，其形成和打开受特提斯构造域和太平洋构造域作用影响^[1-2]。南海的形成和演化过程一直是地球科学的研究热点，不同地质历史时期的构造活动均会在岩石圈中留下痕迹，研究南海及邻区岩石圈结构特征对认识南海的形成和演化规律具有重要意义。近年来，随着地球物理数据的积累，已有多种反演南海莫霍面深度的成果陆续发表，前人在南海及邻区进行了声呐浮标探测^[3]、双船扩展剖面探测^[4]和海底地震仪探测剖面^[5-10]，并计算测量剖面上的莫霍面深度。Braitenberg等^[11]考虑海洋岩石圈弹性形变，利用卫星重力数据计算南海海盆的莫霍面深度; 秦静欣等^[12]利用空间重力异常数据，联合地形和地震剖面数据，通过三维带控制点界面反演算法，计算得到南海及邻区的莫霍面深度。由于采用的数据来源及计算方法不同，南海及邻区的莫霍面显示有不同的变化范围和深度特征。本文在前人研究基础上，采用EGM2008重力场模型联合GOCE卫星重力场模型，计算得到研究区内的高精度重力异常数据，进一步研究南海地区的重力场特征。同时联合海底地形、沉积层厚度及地震探测剖面等资料，计算南海及邻区的莫霍面深度，为综合解释南海及邻区地壳的构造属性和海盆的扩张机制提供依据。

1 数据来源及处理方法 1.1 数据来源

本文研究范围为106°~120°E、8°~22°N。受研究区内地面/船测重力数据质量的限制，根据EGM2008重力场模型计算的自由空气重力异常的短波长精度有待提高^[13-14]。研究发现，由GRACE和GOCE重力卫星观测数据解算得到的重力场模型，在波长小于100 km时的精度较高^[14-15]。为提高研究区的重力场精度，本文将EGM2008重力场模型中的前200阶替换为GO_CONS_GCF_2_SPW_R2^[15](基于GOCE重力卫星观测数据解算得到)的前200阶，替换后的位系数模型计算得到的重力异常为研究区域的自由空气重力异常。在计算过程中对数据进行UTM投影重采样为10 km×10 km的格网数据。

1.2 布格改正

自由空气重力异常是地球内部物质密度分布不均匀性在地表的综合反映。在研究区内，自由空气重力异常可认为是海水、(新、中生代)沉积层、地壳、地幔岩石圈和岩石圈以下物质密度异常的总体响应。利用重力异常反演莫霍面深度，需要从观测重力异常中分离出莫霍面起伏所引起的重力异常。首先对自由空气重力异常数据采用完全布格校正的方法消除海底地形的影响，解算得到研究区的布格重力异常，地形数据采用GEBCO数据^[16]，校正密度为2 670 kg/m^3[11-12]。研究区经过完全布格改正之后的布格重力异常见图 1。

1.3 沉积层重力异常

相关研究表明^[11]，南海海盆及邻区基底覆盖有0.5~12 km厚度不等的新、中生代沉积层。由于基底上、下界面密度差在有些区域可达300 kg/m³，沉积层内部同样存在密度分布不均匀的异常体，因此重力场分离时有必要在观测重力异常中扣除由沉积层内密度异常引起的重力异常。沉积层厚度采用美国NOAA公布的分辨率为5′×5′的格网数据^[17]，密度分布(ρ_z)采用指数密度^[18]分布模型拟合，拟合公式为ρ_z = Φ₀e^－czρ_w+(1－Φ₀e^－cz)ρ_sg，其中, Φ₀为初始沉积物孔隙度，z为沉积层厚度，c为转换参数，ρ_w为海水密度，取值1 030 kg/m³，ρ_sg为填充颗粒物密度。Φ₀、c、ρ_sg的取值根据研究区内海底地震仪探测剖面的P波波速，通过波速-密度转换关系^[19]计算得到：Φ₀ = 0.54, c = 4.75×10^－4m^－1, ρ_sg = 2 650 kg/m³，自由空气重力异常通过地形改正和沉积层密度异常改正后的剩余重力异常见图 2。

1.4 地幔岩石圈重力异常

地壳的密度异常引起的重力异常一般被认为是高频异常，布格异常扣除中、新生代沉积层密度异常引起的重力异常之后，得到的剩余异常主要包含莫霍面起伏引起的重力异常及地幔岩石圈和岩石圈以下物质密度异常引起的重力异常，通过解算该剩余异常的对数振幅谱与波数的对应关系(图 3)，分离出莫霍面起伏引起的重力异常。根据前人研究成果^[11-12]和海底地震仪探测剖面^[3-10]的结果，研究区域莫霍面深度约在7~40 km，图 3显示的25 km对应研究区的平均莫霍面深度。图 3中另外可见平均深度为173 km和55 km的异常波段，平均深度173 km的波段对应地幔岩石圈以下密度异常引起的重力异常，主要是长波长异常。由于南海海盆是新生代大陆破裂、海底扩张形成的边缘海盆，在距今约15.5 Ma时停止扩张^[1-2]，而大陆破裂、海底扩张形成的非稳态的温度场环境引起的岩石圈内的密度异常，是产生平均深度为55 km重力异常波段的主要原因。对于平均深度为173 km和55 km对应的重力异常项，主要采用滤波器(式(3))来扣除。

为控制反演过程中的参数，同时评价反演质量，本文搜集在研究区内的声呐浮标探测^[3]、双船扩展剖面探测^[4]和海底地震仪探测剖面^[5-10]，计算得到71个莫霍面深度资料作为反演质量的检核点，检核点的位置分布见图 1。

2 Parker-Oldenburg反演算法及参数设置

Parker ^[20]提出频率域内计算起伏界面重力异常的公式，Oldenburg^[21]在Parker公式的基础上，将正演公式变为迭代反演形式，提出单界面模式的位场迭代反演公式，即

$ F[\Delta g] = - 2{\rm{ \mathsf{ π} }}G\rho {{\rm{e}}^{ - k{z_0}}}\sum\limits_{n = 1}^\infty {\frac{{{k^{n - 1}}}}{{n!}}} F[{h^n}\left( {\left( {\vec r} \right)} \right] $

(1)

$ F\left[{h\left( {\vec r} \right)} \right] = - \frac{{{\rm{ }}F[\Delta g]}}{{2{\rm{ \mathsf{ π} }}G\rho }}{{\rm{e}}^{k{z_0}}} - \sum\limits_{n = 2}^\infty {\frac{{{k^{n - 1}}}}{{n!}}F[{h^n}\left( {\left( {\vec r} \right)} \right]} {\rm{ }} $

(2)

式中，Δg为观测重力异常，h($\vec r$)为起伏界面深度，G为万有引力常数，ρ为密度参数，k为波数场的模值，z₀为观测平面，F[·]为二维Fourier变换。

根据式(1)和(2)，利用上文数据处理得到的剩余重力异常进行反演，计算南海及邻区的莫霍面深度。在计算过程中，为适应研究区地质构造环境，式(2)中的高频放大因子e^kz₀、观测平面z₀、密度差ρ和迭代的收敛条件需要相关的参数设置。

由于式(2)中的e^kz₀会使重力场的高频噪声放大，影响反演过程的收敛性，本文采用Oldenburg^[21]定义的低通滤波器HCF(k)来控制高频噪声。滤波器表达式为：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{HCK}}\left( k \right) = \\ \left\{ \begin{array}{l} 1, k/2{\rm{ \mathsf{ π} }} < {\rm{WH}}\\ \frac{1}{2}[1 + {\rm{cos}}(\frac{{k-2{\rm{ \mathsf{ π} WH}}}}{{2({\rm{SH-WH}})}})], {\rm{WH}} \le (k/2{\rm{ \mathsf{ π} }}) \le {\rm{SH}}\\ 0, k/2{\rm{ \mathsf{ π} }} > {\rm{SH}} \end{array} \right. \end{array} $

(3)

式中，k为波数，WH和SH为滤波参数，通过重力异常的对数振幅谱估计(图 3)。由图 3可以看出，研究区内的平均莫霍面深度为25 km，即式(2)中z₀ = 25 km，对应的波长范围为70~151 km，对应的高通和低通滤波参数分别为SH = 0.014 2和WH = 0.006 6。

通过对研究区OBS剖面提供的地震波速度计算发现，研究区莫霍面的密度差数值估计在300~550 kg/m³。本文采用迭代计算来获取研究区最佳的密度差ρ，具体计算过程则是通过设定密度差300~550 kg/m³以10 kg/m³的间隔变化，将不同密度差计算的莫霍面深度结果与地震剖面控制点的结果作差，均方根误差最小的一组对应的密度差ρ即为研究区内采用的最终密度差。通过26次反演计算，ρ = 420 kg/m³对应的均方根误差最小，为1.9 km。

迭代的收敛性控制则根据${A_n} = \left| {\frac{{{k^{n- 1}}}}{{n!}}F[{h^n}\left({\left({\vec r} \right)} \right]} \right|$，在每一次迭代过程中对所有的波数k计算A_n的最大值S_n，若S_n/S₁ < 0.01则迭代终止，输出最后的h(${\vec r}$)即为研究区域的莫霍面深度。在本文的计算过程中，Parker-Oldenburg迭代反演算法的最大迭代次数为21，即S₂₁ = 0.009 7。实际上，在迭代过程中S₁ = 5.6，S₂ = 0.51，在迭代算法中，高次项的影响量级在2~3级，计算流程见图 4。

3 计算结果及分析

本文搜集的地震剖面控制点共71个，地震观测数据结果的误差为±2 km^[3-5]，反演结果与控制点观测结果的平均偏差为1.59 km，标准偏差为1.9 km。经统计分析，71个控制点中有86%的点差值在±2 km以内，反演结果比较可靠。计算得到的莫霍面深度分布见图 5。

与前人^[12,22-23]得到的莫霍面深度结果相比，图 5显示的莫霍面在海区的深度有所不同。海盆区莫霍面深度约为8.0~13.0 km，与前人的计算结果约差0.5~1 km；在海盆北部洋陆过渡区，莫霍面深度与前人计算结果约差1~2 km。本文采用基于GOCE重力卫星观测数据解算得到的重力场模型GO_CONS_GCF_2_SPW_R2的前200阶替换EGM2008重力场模型的前200阶，提高了输入数据的重力异常精度。对比采用测高重力异常得到的莫霍面深度结果，测高卫星计算的莫霍面深度与地震剖面控制点计算结果的平均差值为1.67 km，替换后的重力异常使计算结果的可靠性略有提高。根据南海地质特征分离各层重力异常，尤其在扣除沉积层的重力异常时，采用指数形式的变密度公式拟合密度异常，提高剩余重力异常的精度。海盆区由于沉积层厚度较小，莫霍面深度与前人的计算结果差值小于1 km，而在沉积层厚度较大的区域，如莺歌海盆地，沉积层最大厚度大于13 km^[11,24]，莫霍面深度与前人结果相差1.5~3 km。但本文计算的莫霍面深度与海底地震仪探测OBH1996剖面^[24]计算得到的莫霍面深度更为相符，对比前人采用常数密度拟合沉积层密度异常的计算结果，变密度沉积层模型的拟合结果更为可靠。

图 5显示，研究区的莫霍面深度分布具有明显的区域化特征，莫霍面起伏较大，显示为经历复杂地质构造活动形成的边缘海盆地壳结构特征。若以莫霍面深度和地壳厚度作为划分陆壳、洋陆过渡性地壳和洋壳的依据，图 5反映出南海及邻区同时具有这3种类型的地壳，这与前人的研究结果一致。陆地区华南大陆和中印半岛莫霍面深度大于25 km，地壳厚度大于20 km，属于陆壳；西沙群岛、中沙群岛地区莫霍面深度表现为同心圆状隆升，地壳厚度大于20 km，属于陆壳；洋陆过渡壳主要分布在水深为0~3 000 m的区域，地壳厚度小于20 km，莫霍面深度约为13~25 km，莫霍面水平变化梯度大，等值线密集；海盆东侧边缘为马尼拉海沟俯冲带，莫霍面剧烈俯冲，莫霍面深度约为8.0~13.0 km，地壳厚度约为5.0~10.0 km，为洋壳。全区莫霍面深度最小处在西南次海盆，为8.0 km左右。图 5所示的莫霍面深度分布显示研究区内具有明显的海盆、海沟、岛弧等区域化构造特征。

4 结语

本文利用Parker-Oldenburg频率域反演算法，采用重力异常数据反演计算南海及邻区的莫霍面深度分布。采用EGM2008全球重力场模型联合GOCE重力场模型，计算研究区的自由空气重力异常，通过对自由空气重力异常进行布格改正，获取研究区的布格重力异常。在布格重力异常中扣除由沉积层密度异常、地壳内密度异常和地幔岩石圈及岩石圈以下密度异常引起的重力异常后，得到研究区的剩余重力异常，利用剩余重力异常反演南海及邻区的莫霍面深度。采用低通滤波器控制反演算法的收敛性，滤波器截止频率通过计算重力异常对数振幅谱估计。搜集研究区地震剖面计算的莫霍面深度数据来控制反演参数，同时评价反演质量，计算结果与研究区内地震波探测剖面的莫霍面深度平均偏差为1.59 km，标准偏差为1.9 km。结果显示，海盆区莫霍面深度分布约为8.0~13.0 km，陆地区莫霍面深度大于25 km，陆地区与海盆区之间存在洋陆过渡型地壳，莫霍面起伏较大，研究区内的莫霍面深度分布具有明显的区域化特征。