地球物理学报  2016, Vol. 59 Issue (2): 476-483   PDF    
垂直重力梯度反演Moho面的频谱域公式及其应用
叶周润1,2, 柳林涛2, 梁星辉2, 郎骏健2    
1. 合肥工业大学土木与水利工程学院, 测绘工程系, 合肥 230009;
2. 中国科学院测量与地球物理研究所大地测量与地球动力学国家重点实验室, 武汉 430077
摘要: 通过求解引力相等原则下的Fredholm积分方程,可以得到不规则单一密度界面(Moho面)的起伏.本文充分参考了前人的理论研究,推导出扰动垂直重力梯度确定Moho面深度的频谱域表达式,该式具有二次项迭代精度.运用此公式进行了全球Moho面的恢复计算,并将该结果与CRUST1.0模型和GEMMA Moho模型进行了对比和验证.
关键词: 垂直重力梯度     Moho面反演     频谱域公式    
Formula of Moho inversion in the spectral domain using vertical gravity gradient and its application
YE Zhou-Run1,2, LIU Lin-Tao2, LIANG Xing-Hui2, LANG Jun-Jian2     
1. School of Civil Engineering, Institute of Geomatics Engineering, Hefei University of Technology, Hefei 230009, China;
2. State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China
Abstract: In classical gravity research, the determination of Moho is mainly based on the Vening Meinesz-Moritz (VMM) method. This approach assumes that terrain and seawater can be considered as the load on the elastic crust. The Moho depth or crustal thickness can be obtained by solving the Fredholm integral equation which requires a form of expression for gravity gradient observations. However, there is no such accurate expression available now.
In this work, the Fredholm integral equation of Moho depth is converted into the one in the spectrum domain. By considering the second-item influence, we provide a simple and general expression that is applicable for the vertical gravity gradients. In the investigation, we utilize the new methodology to calculate the global Moho depth. The method of the spherical harmonic synthesis is applied to raw gravity gradient data and step-by-step stripping of Moho signal from the GOCO03s model and the prior crust information, respectively. The global crustal thickness is from 6 km to 79 km. The land crust is usually thicker than the normal value. The minimum thickness of the crust is located in the ocean. Especially, because of the collision between in the India and the Eurasian plates, the Tibetan plateau has a very thick crust. With respect to the African region, the crust in the north continent is thinner than the one in the south. Finally, by comparison with CRUST1.0 and GEMMA model, it is demonstrated that our Moho solutions have a good consistency on the global scale with standard deviations of 4.6 km and 3.4 km, respectively.
Key words: Vertical gravity gradient     Moho inversion     Spectral domain formula    
1 引言

重力和重力梯度扰动反映了物质的密度变化和界面的起伏.在经典重力学中,Moho面(地壳厚度)的确定主要依据Vening Meinesz-Moritz(VMM)方法(Moritz,1990).在该方法中,地形或海水等质量被认为是一种加在不断裂而有弹性的地壳上的负载,通过求解引力相等原则下的Fredholm积分方程即可反推出地壳的厚度.Moritz(1990)首先进行了理论推导,但在该表达式中重力异常求导的方向并没有严格指向球心,从而不太符合大尺度的界面反演(Moritz,1990汪汉胜等,1993).Sjöberg(2009)汪汉胜等(1993)修正了导数指向的不足,并分别推导出不同形式的顾及二次项影响的迭代算法.另外,在FFT正演位场界面的基础上(Parker,1973),Oldenburg(1974)提出确定起伏密度界面深度的迭代法,该方法被广泛应用于局部层面反演中,如海洋深度计算和青藏地区的地壳厚度计算(Oldenburg,1974Braitenberg et al.,2000柯小平,2006Shin et al.,2007).此外,重力反演地壳界面的主要进展还包括:(1)伴随GPS 定位技术的发展,重力扰动逐渐取代重力异常并被应用于地球物理学研究中(Sjöberg,2013).Tenzer和Bagherbandi(2012)使用重力扰动数据得到了结果更优的全球Moho面.Sjöberg(2013)推导出联合重力异常和重力扰动反演Moho面的公式.(2)随着GOCE卫星的成功发射,空间分辨率更高的卫星重力梯度观测值被用于全球地壳厚度的恢复研究(Ebbing et al.,2011Hirt et al.,2012).Reguzzoni等(2013)首先给出一阶近似下的球面迭代算法.Bagherbandi和Eshagh(2012)采用离散最小二乘方法模拟反演了伊朗地区的地壳厚度.(3)地壳先验模型的联合应用.Tenzer等(2009)首先提出基于地壳先验信息的逐步剥离策略,该方法经过地形、海水、冰层、沉积层和地壳残余密度差改正后,可以有效地减小密度差对地壳界面恢复的影响.它被认为是重震数据的一种融合解算方法,用来改善地震先验模型在资料欠缺地区的不足.

Reguzzoni等(2013)虽然给出了基于垂直重力梯度反演全球地壳厚度的迭代法,但在其公式推导中的截断误差只考虑了一阶近似.同时,关于Moho面信号的提取采用了精度较差的点质量正演模型. 为改善上述提及工作的不足,本文在Moritz,Wang,Sjberg和Reguzzoni等研究思路的基础上,推导出顾及二次项影响的频谱域迭代公式,并在正演计算过程中采用了精度更高的球谐谱方法.

2 计算方法 2.1 理论推导

图 1所示,在大地水准面以下存在一个深度为D的不规则界面.假设该层面的正常厚度为D0,同时设定该界面与外界的密度差为Δρc/m,则其产生的重力位扰动可以表达为(通常忽略离心力的影响)(Sjöberg,2009):

图 1 不规则界面示意图Fig. 1 Schematic of irregular interface

式中,G为万有引力常数;(r,θ,λ)和(r′,θ′,λ′)分别表示计算点和流动点的球心距离、余纬和经度,且两点间的距离;R是地球平 均半径;dσ′=sinθ′dθ′dλ′; 并且有cosψ=cosθcosθ′+sinθsinθ′cos(λ-λ′).这里指出,后续的公式推导皆基于球形近似.因为椭球近似下的各种函数表达复杂,而且在通常情况下球近似的误差也很小,可以满足大部分的大尺度的物理大地测量研究(Tsoulis,2001Novák and Grafarend,2006Wild and Heck,2008Tenzer et al.,2009Grombein et al.,2010Liang et al.,2014),例如,利用Tesseroid单元体进行GOCE卫星的地形校正.同时考虑到垂直重力梯度较其他分量在地球科学中的应用更广,并且它在Moho面反演过程中的数学处理较其他分量简易很多,故本文只针对该分量进行了相关的理论推导工作.

根据Heiskanen和Moritz(1967)可知,当r≥R时,式(1)中关于l-1(r,θ,λ)的表达式可以用球谐函数展开为:

式中,Pn(cosψ)表示面谐函数.将式(2)代入式(1),则有:

根据定义,垂直重力梯度为重力位在球心方向的二阶导数,则可以得到:

考虑到定积分的可加性质,(4)式中的积分项可以改写为:

为了后续表示的方便,我们定义:由部分计算得到的扰动垂直重力梯度记为Γrr0;类似,Γrr1和Γrr2分别表示从部分计算得到的扰动垂直重力梯度值.Γrr0包括正常界面和扰动界面产生的重力梯度值,可以通过信号提取的方法从测量值中得到(Tenzer et al.,2009Reguzzoni et al.,2013);对于Γrr2,当正常界面厚度确定时,该部分实际是一个规则体产生的重力梯度.在进行全球计算时可以用Tesseroid单元体的空间域或频谱域方法得到相应数值(Novák and Grafarend,2006Wild and Heck,2008Grombein et al.,2010);对于Γrr1,这里包含待求的扰动面深度信息项D.经过上面的数学技巧处理之后,我们的问题便转换为如何从部分的表达式找到一个关于Γrr1与未知Moho面深度D的线性数学关系.

定义D/R=τ,并对r′进行积分,对应的部分的积分表达式可以改写为:

考虑到:(1)τ的最大值是100 km/6371 km≈0.0157(100 km表示最大Moho面深度,6371km表示地球平均半径),故进行有限级数展开即可满足计算精度.(2)据Sjöberg(2009)指出,τ级数展开的三次 项最大近似误差是1003/63712 km≈25 m,略去此项可满足大多数情况的实际应用.因此对1-(1-τ)n+3展开到二次项即可,具体展开如下:

引入变量hrrΓ=(Γrr0+Γrr2)/(GΔρc/m),综合式(6)和(7),则式(4)可改写为:

根据球谐谱分析理论,式(8)可以等价转换为:

式中,hnmΓrr和τnm分别表示hΓrr和τ经过球谐分析得到的谱系数;而(τ2)nm是关于τ的二次项的球谐谱系数.继续引入Dnm/R=τnm,则由式(9)可以得到关于Moho面深度的频谱系数Dnm的表达式:

考虑到有如下数学关系式成立(Heiskanen and Moritz,1967):

式中,Yn(θ,λ)是球谐分析与综合理论中的球面函数;V(θ,λ)V(θ′,λ′)是对应积分点和流动点在球面上的普通变量函数.综合式(10)和(11),进行球谐系数恢复计算,则可以得到扰动垂直重力梯度恢复Moho面的频谱域公式:

式中,D(θ,λ)D(θ′,λ′)分别为计算点和流动点的界面深度.虽然(12)式是个迭代公式,但由于它有类似于Molodensky公式的快速收敛性质,Sjöberg(2009)认为该类型公式在具体实用计算时可以把它当作直接表达式.

2.2 Moho面计算流程

图 2所示,本文进行Moho面反演的基本流程如下:首先从原始观测数据中扣除球谐函数低阶项的影响,由此得到的结果主要包含Moho面以上的扰动信息;再结合地壳先验密度模型和重力梯度正演方法,通过扣除(补偿)计算得到单一密度界面的扰动值;之后加上正常地壳厚度的重力梯度值即反演初值;最后利用前面介绍的频谱域反演公式便可解算出Moho面深度.

图 2 Moho面反演流程Fig. 2 Process of Moho inversion

图 2中所提及的Moho面信号提取策略介绍如下:如图 3a所示,Moho面以上的地球结构大体分层(依次)有陆地岩石(1)、海水(2)、冰层(3)、沉积层(4)和密度不均匀的地壳(5).如果我们要用Vening Meinesz-Moritz方法解算出地壳厚度,首先需要得到如图 3b所示的一个单密度界面.为此,这里采取结合地壳先验模型的逐步剥离方案,这种处理策略得到的最终结果被认为和地震数据解算的Moho面具有很强的相关性(Tenzer et al.,2009).对于GOCE观测数据而言,则可采用(13)式进行有效信息逐步提取:

图 3 (a) Moho面上部地壳结构示意图; (b) Moho面示意图Fig. 3 (a) Crustal structure above Moho; (b) Schematic of Moho

式中,ΓrrM表示上地幔以上的扰动重力梯度;ΓrrT、ΓrrB、ΓrrI、ΓrrS和ΓrrC分别表示由陆地地形产生的重力梯度及海水、冰层、沉积层和不均匀地壳相对于岩石密度产生的补偿重力梯度.

2.3 重力梯度球谐合成和正演公式

从位系数模型恢复垂直扰动重力梯度的信息可 以采用(14)式(罗志才,1996Tsoulis,2001Novák and Grafarend,2006Zhu,2007):

式中,GM代表万有引力常数和地球总质量的乘积(3.986005×1014 m3·s-2);,且CnmSnm是GOCE模型位系数和正常位系数相减后得到的扰动位系数;Pnm(cosθ)是完全正则化的连带勒让德函数;NminNmax分别为起始和最大阶数.

Tsoulis(2001)Novák和Grafarend(2006)可知,频谱域的重力梯度正演公式和式(14)在形式上一致,但Vnmβ部分需要从扰动物质(如地形和海水等)计算得到.为方便表达,定义:

式中,i∈{1,2,3};ρ表示扰动物质密度; ρ表示地球平均密度(5500 kg·m-3);huphlw分别表示积分体高程上下限,具体数值以球面(大地水准面)为基准,即球面以上为正,反之为负;Hnmi表示球谐谱系数,可由式(17)计算得到,将Hnmi代入式(16)即可得到最终重力梯度正演计算所需的谱系数.

3 实验计算与分析

本文原始重力梯度数据来源于GOCO03s模型,该模型具有低阶GRACE和中高频GOCE数据互补优势(Mayer-Gürr et al.,2012). 在球谐合成时,选取GRS80模型为正常重力场参数(Mortiz,1980),选定GOCE轨道平均值255 km为计算高度,并考虑到低阶信号主要包含地核和下地幔物质的密度不均匀信息,依据Bagherbandi和Sjöberg(2012)方剑(1999)结论,在此扣除16阶以内的球谐合成值.在地壳信息校正部分,采用2013年7月最新发布的CRUST1.0模型(Laske et al.,2013).岩石、海水、冰层、沉积层和地壳密度数据全部来自该模型,并取壳幔密度差为485 kg·m-3.综合GOCE观测数据的空间分辨率和地壳先验模型精度,文中计算结果分辨率设为1°×1°.考虑到空间域和球谐谱对应关系,式(12)和(14)中球谐函数阶次的最大值取180.

图 4展示的是全球扰动垂直重力梯度逐步剥离结果,对应的统计数值反映在表 1中.这里ΓrrM、ΓrrMT、ΓrrMTB、ΓrrMTBI、ΓrrMTBIS和ΓrrCS依次表示从原始数据扣除低阶项、陆地地形、海水补偿、冰层、沉积层和不均匀地壳的扰动重力梯度.从图 4a可知,虽然在球谐合成时进行了长波信息滤除,但该图中的扰动重力梯度图形轮廓变化依旧分明,尤其在青藏高原和安第斯山脉地区表现明显.这说明扰动信息的来源主要是陆地地形起伏、海水的密度补偿和地幔及其上部的密度异常等.如表 1所示,上面提及的扰动源造成的影响极值接近1.5E(1E=10-9 /s2).图 4b4c展示的是经过陆地地形和海水密度补偿校正后的结果.在陆地部分,经过地形高频信号扣除后,总体轮 廓趋于平缓,显示出大地水准面以下的主要地质构 造信号;在海洋部分,经过海水密度补偿产生的重力梯度平滑后,显示出明显的海沟构造信号.图 4d展示的是经过冰层校正的结果.由于南北极的冰层较厚,该项校正可以改善两极地区的地壳恢复精度.图 4e展示的是扣除沉积层影响后的结果.从该图可知,沉积层不仅在全球范围内分布广泛,而且极值区间也可达2E左右.这是由于虽然沉积层密度相对于岩石密度波动不大,但厚度可达数千米以上.所以此项改正不仅在地壳厚度研究中要加以重视,而且在地幔密度反演中也是一项重要研究.图 4f展示的是经过地壳密度差校正后的结果,也即是如图 3b所示的单一密度界面产生的扰动垂直重力梯度值.如表 1所示,由于地壳密度残差造成全球重力梯度影响最大值达到近4E,而由Moho面扰动产生的垂直重力梯度变化范围为[-8.65~4.88]E.从图 4f结果显示,在陆地地区该层面产生的扰动垂直重力梯度多为负值而海洋地区多为正值,并且陆地和海洋表现出明显的界限.

图 4 全球扰动垂直中立梯度逐步玻璃结果Fig. 4 Globel maps of the dtep-wise corrected vertical gradients of the gravity disturbaces

表 1 扰动垂直重力梯度逐步剥离结果统计 Table 1 Statistics of the step-wise corrected vertical gradients of the gravity disturbances

图 5展示的是Moho面反演公式二次项的影响.总体而言,二次项的贡献最大可达14 km,尤其在重力梯度变化比较剧烈的青藏高原等地区.图 6展示的是全球Moho面分布图.本文解算得到的全球地壳厚度分布区间为[6~79] km,陆地地区的地壳厚度多大于正常地壳,最小的地壳厚度区域位于海洋,该区域结果多小于正常值.其中,青藏地区由于印度板块和欧亚板块的挤压,呈现明显的边缘浅和中部深特点;非洲地区则呈现出南部地壳厚于北部的特征.表 2展示的是本文解算结果分别与CRUST1.0和GEMMA Moho模型的对比统计,相应的模型对比偏差平均值是-3.9 km和-2.7 km,模型之间差异的标准差是4.7 km和3.5 km.与CRUST1.0模型的差异主要集中在非洲、南美和两极地区,这是由于该模型在这些区域缺少准确的地震数据.与GEMMA Moho模型的差异主要集中在青藏高原地区,该模型在此区域解算的地壳厚度值较大,最大值可达100 km左右,该区域值均大于本文和CRUST1.0模型结果.虽然存在局部差异,但 本文解算结果与CRUST1.0模型对比的平均值和标准偏差均优于Tenzer和Bagherbandi(2012)Reguzzoni等(2013).

图 5 Moho面反演公式二次项的影响Fig. 5 Effect of second item in the formula of Moho inversion

图 6 全球Moho面解算结果Fig. 6 Result of global Moho

表 2 模型结果对比统计 Table 2 Statistics of the comparison with other models
4 结论

本文的主要工作是推导了新的扰动垂直重力梯度确定界面深度的频谱域算法,并用该算法进行了全球地壳厚度的反演工作.在理论创新部分,将求解界面深度的Fredholm积分方程转为易于求解的频谱域等式进行解算,所推导的界面深度反演表达形式简洁而通用.由于考虑了级数展开二次项的影响,所得公式保证了理论上的精度.在实验部分,结合最新的地壳先验模型CRUST1.0,采用球谐谱正演重力梯度的方法对Moho面的信号进行了逐步剥离; 并将反演得到的全球Moho面深度分别与CRUST1.0 模型和GEMMA Moho模型进行了对比分析,从文中相关统计显示,本文解算结果和两种模型在全球范围内具有较良好的一致性.

Moho面的恢复计算在地球物理学研究中是比较重要和复杂的工作.本文的研究只是基于传统重力学的研究思路,在理论算法和正演模型方面提出了合理的改进措施,在全球反演过程的其他细节方面借鉴了前人的较成熟结论,比如低阶长波项的扣除,Moho面信号的提取和壳幔密度差的设定,等.并且在其他方面还需继续深入,比如,其他重力梯度分量的界面恢复公式的理论推导,解算结果的详细地球物理学解释和综合地震数据进行Moho面的联合结算,等.

致谢    感谢审稿专家的意见!

参考文献
[1] Bagherbandi M, Sjöberg L E. 2012. Non-isostatic effects on crustal thickness:a study using CRUST2.0 in Fennoscandia. Physics of the Earth and Planetary Interiors, 200-201:37-44.
[2] Bagherbandi M, Eshagh M. 2012. Crustal thickness recovery using an isostatic model and GOCE data. Earth, Planets and Space, 64(11):1053-1057.
[3] Braitenberg C, Zadro M, Fang J, et al. 2000. The gravity and isostatic Moho undulations in Qinghai-Tibet plateau. J. Geodyn., 30(5):489-505.
[4] Ebbing J, Bouman J, Gradmann S, et al. 2011. Use of GOCE gravity gradient data for lithospheric modeling-a case study for the NE Atlantic margin.//AGU Fall Meeting Abstracts. AGU.
[5] Fang J. 1999. Global crustal and lithospheric thickness inversed by using satellite gravity data. Crustal Deformation and Earthquake (in Chinese), 19(1):26-31.
[6] Grombein T, Seitz K, Heck B, et al. 2010. Untersuchungen zur effizienten Berechnung topographischer Effekte auf den Gradiententensor am Fallbeispiel der Satellitengradiometriemission GOCE. KIT Scientific Publishing.
[7] Heiskanen W A, Moritz H. 1967. Physical geodesy. Bulletin Géodésique, 86(1):491-492.
[8] Hirt C, Kuhn M, Featherstone W E, et al. 2012. Topographic/isostatic evaluation of new-generation GOCE gravity field models. J. Geophys. Res., 117(B5):B05407.
[9] Ke X P. 2006. Crustal structures of Tibetan plateau from the research of gravity forward and inversion [Ph. D. thesis] (in Chinese). Wuhan:Institute of Geodesy and Geophysics, Chinese Academy of Sciences.
[10] Laske G, Masters G, Ma Z T, et al. 2013. Update on crust1. 0-a 1-degree global model of Earth's crust.//EGU General Assembly. Vienna, Austria.
[11] Liang Q, Chao C, Li Y G. 2014. 3-D inversion of gravity data in spherical coordinates with application to the GRAIL data. J. Geophys. Res., 119(6):1359-1373.
[12] Luo Z C. 1996. The theory and method for the determination of earth gravity field from satellite gravity gradient data [Ph. D. thesis] (in Chinese). Wuhan:Wuhan University of Surveying and Mapping Technology.
[13] Mayer-Gürr T, Rieser D, Höck E, et al. 2012. The new combined satellite only model GOCO03s.//Abstract, GGHS2012. Venice.
[14] Mortiz H. 1980. Geodetic reference system 1980 (GRS80). Bulletin Géodésique, 54(3):395-405.
[15] Moritz H. 1990. The Figure of the Earth:Theoretical Geodesy and the Earth's Interior. Karlsruhe:Wichmann.
[16] Novák P, Grafarend E W. 2006. The effect of topographical and atmospheric masses on spaceborne gravimetric and gradiometric data. Stud. Geophys. Geod., 50(4):549-582.
[17] Oldenburg D W. 1974. The inversion and interpretation of gravity anomalies. Geophysics, 39(4):526-536.
[18] Parker R L. 1973. The rapid calculation of potential anomalies. Geophys. J. Int., 31(4):447-455.
[19] Reguzzoni M, Sampietro D, Sansò F. 2013. Global Moho from the combination of the CRUST2.0 model and GOCE data. Geophys. J. Int., 195(1):222-237.
[20] Shin Y H, Xu H Z, Braitenberg C, et al. 2007. Moho undulations beneath Tibet from GRACE-integrated gravity data. Geophys. J. Int., 170(3):971-985.
[21] Sjöberg L E. 2009. Solving Vening Meinesz-Moritz inverse problem in isostasy. Geophys. J. Int., 179(3):1527-1536.
[22] Sjöberg L E. 2013. On the isostatic gravity anomaly and disturbance and their applications to Vening Meinesz-Moritz gravimetric inverse problem. Geophys. J. Int., 193(3):1277-1282.
[23] Tenzer R, Hamayun K, Vajda P. 2009. Global maps of the crust 2.0 crustal components stripped gravity disturbances. J. Geophys. Res., 114(B5):B05408.
[24] Tenzer R, Bagherbandi M. 2012. Reformulation of the Vening-Meinesz Moritz inverse problem of isostasy for isostatic gravity disturbances. Int. J. Geosci., 3(5):918-929.
[25] Tsoulis D. 2001. A comparison between the Airy/Heiskanen and the Pratt/Hayford isostatic models for the computation of potential harmonic coefficients. J. Geod., 74(9):637-643.
[26] Wang H S, Chen X, Yang H Z. 1993. An iterative method for inversion of deep large-scale single density interface by using gravity anomaly data. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese), 36(5):643-650.
[27] Wild F, Heck B. 2008. Topographic and isostatic reductions for use in satellite gravity gradiometry.//VI Hotine-Marussi Symposium on Theoretical and Computational Geodesy. Berlin Heidelberg:Springer, 49-55.
[28] Zhu L Z. 2007. Gradient modelling with gravity and DEM [Ph. D. thesis]. Columbus:The Ohio State University.
[29] 方剑. 1999. 利用卫星重力资料反演地壳及岩石圈厚度. 地壳形变与地震, 19(1):26-31.
[30] 柯小平. 2006. 青藏高原地壳结构的重力正反演研究[博士论文]. 武汉:中国科学院测量与地球物理研究所.
[31] 罗志才. 1996. 利用卫星重力梯度数据确定地球重力场的理论和方法[博士论文]. 武汉:武汉测绘科技大学.
[32] 汪汉胜, 陈雪, 杨洪之. 1993. 深部大尺度单一密度界面重力异常迭代反演. 地球物理学报, 36(5):643-650.