1 引 言

GOCE(Gravity Field and Steady-state Ocean Circulation Explorer)是一颗用于恢复地球重力场和探测静态洋流的重力梯度卫星，其优于100km的空间分辨率观测值能改善地球重力场的中高频信息，因此在精化全球大地水准面^[1]的同时也可以为地球物理学研究带来更多信息。GOCE被认为可以捕捉到壳幔边界信号^[2]，对岩石圈的密度敏感^[3]，因此在盆地密度、地质构造、Moho面^{[4 ,5, 6]}等具体研究中具有巨大潜力。

壳幔起伏和岩石圈密度是重要的地球内部结构研究，而Moho面的确定可以为地壳均衡和岩石圈有效弹性厚度研究提供相应依据^[7]，并且重点地区(如青藏高原等)的Moho面结构可以为此区域的地球动力学研究提供参考^{[8 ,9]}。相对于地震方法，利用卫星重力确定地壳厚度可以弥补地面测量重力和地震资料的缺陷^{[10, 11, 12]}，尤其是重力梯度可以改善特定频段信号在地球物理应用中的精度^[13]，因此在GOCE卫星发射后非洲等地区成为研究热点^[3]。另外重力测量值可以和地震资料进行联合反演，使解算结果相互约束从而提高反演精度。同时重力资料的分析方法可以用于其他星球，如月球等无地震数据的星球地壳厚度和地质构造研究^[14]。本文的研究重点主要在于如何从GOCE观测值中提取出用于Moho面反演的初值。在理论上借鉴和改进了前人的研究思路：①伴随GPS定位技术发展，利用重力扰动(纯重力异常)分析扰动源更加方便^[15]，相对于重力异常，重力扰动用于Moho面反演也显现出更高精度^{[15, 16, 17, 18]}，因此本文研究中所提取的信号都是扰动重力梯度值；②比点质量模型正演精度更高的空间域Tesseroid单元体和频谱域球谐分析与综合方法被应用于重力梯度正演^{[19, 20, 21]}；③重力观测值结合地壳先验模型的逐步剥离方案被应用于重力梯度信号的提取中^{[22, 23]}，该处理策略减少了Moho面反演中由于地壳密度不均带来的误差影响。在具体计算中采用了最新CRUST1.0先验模型^[24]。本文最后运用上述方法提取的扰动垂直重力梯度值，进行了全球Moho面的恢复工作，并与现有模型进行了对比验证。

2 方法介绍 2.1 Moho面信息提取策略

如图 1所示，Moho面以上的地球结构大体分层主要有陆地岩石(1)、海水(2)、冰层(3)、沉积层(4)和除上述之外的密度不均匀地壳(5)。

如果要用Vening Meinesz-Moritz^{[15, 16, 17 ,18]}方法解算出地壳厚度，首先需要得到如图 2所示的一个单密度界面。

为此，本文采取结合地壳先验模型的逐步剥离方案，这种处理策略得到的最终结果被认为和由地震数据解算的Moho面具有很强的相关系数^{[22 ,23]}。对于GOCE观测数据而言，则可采用式(1)进行有效信息逐步提取

式中，i,j∈1,2,3；Γ_ij表示梯度分量；δΓ_ij表示上地幔以上的扰动重力梯度；Γ_ij^T、Γ_ij^B、Γ_ij^I、Γ_ij^S、Γ_ij^C分别表示由陆地地形产生的重力梯度及海水、冰层、沉积层和不均匀地壳相对于岩石密度产生的补偿重力梯度。

在壳幔信号分离中，依托于GOCE重力梯度模型，根据文献^[25]结论在此扣除17阶以内主要存在于上地幔以下的长波信息。对于重力梯度正演部分，采用球坐标系下的空间域Tesseroid和频域球谐方法^{[19, 20, 21]}，以克服意大利GEMMA Moho团队的点质量模型计算误差^[26]。在地壳信息校正部分，采用2013年7月最新发布的CRUST 1.0模型。该模型的陆地地形，海洋和冰层资料来自ETOPO1^[27]，同时沉积层和地壳密度的空间域分辨率也从2°提高到1°^[24]。

2.2 GOCE模型计算扰动重力梯度值

从GOCE球谐位系数模型恢复扰动重力位信息可以采用式(2)的形式

式中，r、θ和λ分别表示计算点的球心距离、地心纬度和经度；GM代表万有引力常数和地球总质量的乘积(3.986005e14m³/s²)；R是地球平均半径(6371km)；C_nm和S_nm是GOCE模型位系数和正常位系数相减后得到的扰动位系数；P_nm是完全正则化的连带勒让德函数；N_min和N_max分别为起始和最大阶数。

考虑到便于后续空间域和频谱域方法比较，这里选用当地右手坐标系LNOF(local north-oriented frame)。如果要转换到此坐标系，则可运用式(3)表述^[28]

式(3)中对应的关于位的一阶和二阶导数详见文献^[28]。

2.3 空间域和频谱域正演扰动重力梯度 2.3.1 空间域Tesseroid单元体方法

在空间域计算中，笔者选取与地球表面地形非常吻合的Tesseroid单元体作为积分模型。根据文献^[20]模拟试验结论，相对于点质量模型，Tesseroid单元体正演精度可以提高3个量级；另外文献^[14]中试验表明0.25°的Tesseroid单元体对理论月球模型的拟合误差为0.001%，并且文献^[19]和^[20]中已经验证了该单元体用于扣除GOCE卫星扰动物质影响的可行性。

文中计算公式来自泰勒级数展开方法^[19]，若取Tesseroid单元体的中心点：r₀=(r₁+r₂)/2，φ₀=(φ₁+φ₂)/2，λ₀=(λ₁+λ₂)/2，并展开到4阶无穷小，则经过相应数学计算，可以得到

式中，计算基准面为半径R球面；G表示万有引力常数；ρ表示扰动物质密度； (r,φ,λ)、(r′,φ′,λ′)、(r₀,φ₀,λ₀)、(r₁,φ₁,λ₁)、(r₂,φ₂,λ₂)表示流动计算点、积分点、Tesseroid中心点、上限和下限点的球心距离、余纬和经度；Δr、Δφ、Δλ分别表示Tesseroid在球心、纬度和经度方向的间距；

分别表示在球心、纬度和经度方向的Γ_ij函数二阶导数位于Tesseroid中心点时的数值。o(Δ⁴)表示泰勒展开的截断误差。

由于从引力位表达式泰勒展开到二阶导数(引力梯度)比较繁琐，文献^[19]将核函数分开求导，则有

式中，φ表示纬度，且φ=90°－θ；i,j,k∈{1,2,3}；Δ_xi、Δ_xj、Δ_xk分别代表x、y、z方向的间距；并且当i=j时，δ_ij=1；i≠j时，δ_ij=0。
式(6)中其他表达式有 (10)式(7)—式(10)中，l和ψ分别表示积分点和计算点之间的球面距离和夹角；其他变量具体表达详见文献^[19]。

2.3.2 频谱域球谐分析与综合方法

设一物体的内部质点分布连续，则该物体外一点的引力位为

可将距离l函数进行球谐展开，然后对r′进行积分，进行级数3次项展开，并假定用球面代替大地水准面，则经过运算后有如文献^[29]形式的表达 式中 式中,i∈{1,2,3}；ρ表示地球平均密度(5500kg/m³)；h_up和h_lw分别表示积分体高程上下限，具体数值以球面(大地水准面)为基准，即球面以上为正，反之为负；H_nmⁱ表示球谐谱系数，可由式(15)计算得到，将H_nmⁱ代入式(14)即可得到最终计算需要谱系数。则上述就是由扰动物质通过球谐谱域计算扰动位的一般表达形式，之后的各重力梯度分量则可由式(3)得到。

2.4 Moho面恢复谱域计算公式

为了验证前述Moho面扰动重力梯度信息的提取结果，笔者利用Γ_zz分量进行了试验计算。从Γ_zz分量进行Moho面恢复的频谱域公式为

式中，dσ′=sinθ′dθ′dλ′；D(θ,λ)和D(θ′,λ′)分别为计算点和流动点的Moho面深度；h_nm^Γ为频谱域的Moho面解算初值。

3 试验计算及其结果分析

本文试验原始重力梯度数据来源于GOCO03S模型，该模型在低阶以GRACE为主，高阶采用GOCE从2009年11月开始历时18个月的中高频数据^[30]。在球谐合成时，正常重力场参考GRS80模型。选定GOCE轨道平均值255km为计算高度，并忽略空气微小值影响^[21] 。岩石、海水和冰层密度分别取值2670kg/m³、1025kg/m³920kg/m³。沉积层和地壳密度数据全部来自CRUST1.0模型。综合GOCE观测数据的空间分辨率和地壳先验模型精度，文中计算结果分辨率设为1°×1°。考虑到空间域和球谐谱对应关系，式(2)和式(12)中n最大值取180。

图 3表示以Γ_zz为例的空间域与频谱域方法在随网格分辨率改变时的计算所需时间和两者结果差异变化图(648个10°×10°等间距采样点)。左侧纵轴表示计算时间，右侧纵轴表示两种方法绝对值之差的平均值。如图 3所示，随着网格分辨率的加密，空间域方法时间需求远高于频谱域方法，并且两者的计算时间差异呈非线性增长的方式。同时，计算方法的差异表现出随格网分辨率加密而减少的趋势。当格网点加密到1°×1°时，采样点绝对值之差的平均值仅为0.0038E(1E=10^-9/s²)。表 1所展示的是空间域Tesseroid和频谱域球谐分析与综合方法的结果比较。表 2所示为Moho面扰动重力梯度信息逐步提取的结果统计。表中T、B、I、SED(S)和CRD(C)分别代表陆地地形、海水、冰层、沉积层和不均匀地壳。GTBISC表示扣除长波影响后的扰动重力梯度(G)通过逐步剥离步骤得到的含Moho面信息的扰动重力梯度值。如表 1所示，两种方法的平均值差异在10^-3量级左右。此外，笔者统计了表 1中的差异值对于最终扰动值(表 2中GTBISC项)的相对平均幅度，其结果是：2.2%(Γ_xx)、3.4%(Γ_yy)和1.3%(Γ_zz)。

图 4表示在GOCE轨道平均高度扣除长波影响后的Γ_zz分量结果(1 E=10^-9/s²)。如图 4所示，虽然在球谐合成时进行了长波信息滤除，但图 4中的扰动重力梯度图形轮廓变化依旧分明，尤其在青藏高原和安第斯山脉地区表现明显。这说明扰动信息的来源主要是陆地地形起伏，海水补偿密度和地幔及其上部的密度异常等。如表 2 所示，上面提及的扰动源造成的Γ_zz影响接近1.5E。

图 5反映的是经过陆地地形、海水和冰层校正后的结果。图形轮廓相对比图 4已然发生较大改变，这些显著变化表明陆地和海洋地形是大成分扰动源。其实在地壳厚度粗略估算中，一般也只扣除陆地和海洋地形影响。这些在表 2中也有所体现，Γ_zz GTBI项极值较于G项数值已扩大3倍以上。

图 6和图 7表现的是沉积层和地壳补偿密度的正演结果。表 3是相应的结果统计。如图 6所示，沉积层不仅在全球范围内分布广泛，而且极值区间也可达2E左右。这是由于虽然沉积层密度相对于岩石密度波动不大，但厚度可达数千米以上。所以该项改正不仅在地壳厚度研究中要加以重视，而且在地幔密度反演中也是一项重要研究。考虑到Moho面实质是同一密度单层厚度反演，地壳中残余密度差影响要尽量消除，而图 7计算就是为解决此问题。如图 7所示，以青藏高原为例，此区平均密度值明显小于正常岩石密度。如表 3所示，由于地壳密度残差造成全球重力梯度影响最大值达到近4E。

图 8所示是经过全部校正后的结果，也即笔者采用重力梯度测量值进行Moho面反演的初始值。如表 2 所示，在GOCE平均轨道高度最终Moho面信息Γ_zz分量位于区间[-8.7，4.9]E。图 9表现的是采用文中计算的初始值解算的全球Moho面。结果显示陆地地区Moho面Γ_zz扰动多为负值，解算的Moho面多大于正常地壳厚度，海洋地区则反之。其中青藏地区由于印度板块和欧亚板块的挤压，呈现明显的边缘浅，中部深的特点。这些结果符合均衡理论，使得Moho面的深度和地形负载呈现强相关性。

4 结 论

依据球坐标系的空间域和频谱域正演方法，结合地壳先验模型，本文讨论了从GOCE重力梯度模型提取Moho面扰动信息的步骤和方法，并提供了适合该层面反演的最终结果。从文中的讨论和数据分析可以得到如下结论：

(1) 在文中1°×1°空间分辨率下，空间域和频谱域方法结果接近，同时也验证了程序和计算的正确性，这两种方法除了适合本文信息提取工作，同时也适用于其他扰动源改正。

(2) 在Moho面扰动信息剥离过程中，陆地地形和海水的影响最大，并且沉积层和地壳密度差的校正不可忽略。

(3) 在通过频域公式恢复的Moho面结果与地震模型和基于GOCE数据已解算的Moho面结果的对比后验证了文中的扰动信息提取计算的可靠性。

致谢：感谢斯图加特大学F.Wild博士、中国地质大学杜劲松博士、米兰理工大学D. Sampietro博士、武汉大学Robert Tenzer教授和中科院测地所方剑教授等在该论文上的帮助。