地球物理学进展  2014, Vol. 29 Issue (1): 61-72   PDF    
均衡模型在莫霍面反演中的应用综述
邢健1,2, 郝天珧13, 徐亚1, 秦静欣    
1. 中国科学院地质与地球物理研究所, 北京 100029;
2. 中国科学院大学, 北京 100049;
3. 中国科学院遥感应用研究所, 北京 100094
摘要:莫霍面是重要的地质与地球物理界面, 对其深度分布的研究是地球物理领域的重要内容.莫霍面深度分布与地球均衡重力场之间有重要关系, 但长期以来部分研究人员在使用重测-均衡方法探究莫霍面深度时或多或少地忽视了参考均衡模型的重要性, 造成了均衡补偿模式理论发展与重测-均衡方法研究莫霍面深度之间的不同步.本文介绍了世界上常用的几种均衡补偿模式当前的发展状况, 对均衡解释中的一些争议问题和应用性问题进行了讨论, 并针对莫霍面深度反演问题, 指出了传统重测-均衡方法尚需进行的改进, 并从五个方面阐释了均衡理论对传统反演方法的贡献, 并探讨了下一步可继续探究的方向.
关键词莫霍面     均衡     重力    
Review of the application of isostatic models in the inversion of Moho depth
XING Jian1,2, HAO Tian-yao1 , XU Ya1, QIN Jing-xin3    
1. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences, Beijing 100094, China
Abstract: Moho discontinuity is an important boundary in the realms of geology and geophysics, the distribution of whose depth is of great significance in geophysics, which is also what geologists and other researchers from a variety of subjects are concerned about. In spite of the fact that numerous methods for dealing with the inversion of Moho depth are widely used, however, some researchers have for a long time more or less neglected referencing isostatic models during their researches in the distribution of Moho depth by means of gravimetric-isostatic methods, which results in the failure of synchronous development in the researches between isostatic compensation theory and gravimetric-isostatic methods as well as of taking advantage of geological and geophysical data presented as a way to ensure the correctness of the inversion results. In this paper, we brief readers on the status quo of the development of several isostatic compensation models including Pratt-Hayford model, Airy-Heiskanen model, Vening Meinesz model, Vening Meinesz-Moritz model, experimental isostatic models, dynamic compensation models and other isostatic models that are commonly used in the world and discuss some controversies over the isostatic compensation mechanism, the explanation of isostasy and its applicability in the research of the crust. In addition, apropos of the inversion of the distribution of Moho depth, we specifically focus on gravimetric-isostatic methods, list the weaknesses of traditional inversion methods and their corresponding improvement needed, expatiate in five parts on the contribution of isostatic theory on traditional inversion methods. The five parts in detail are the perfection of traditional isostatic models, application of new isostatic models, taking into consideration the inhomogeneity of density in the crust and the mantle, the research on dynamic compensation and methods in dealing with non-isostatic areas. The discussion of the contribution above shows that, globally, gravimetric-isostatic methods have been developed smoothly with more and more complicated models of variable densities, much closer than they did before in reflecting the reality under the ground, while at the same time a lot of models have not been applied yet in practice for the inversion of the Moho depth. In the end, for the sake of both improving the quality of inversion results and making full use of data that are not only globally covered but also favored on the part of their precision, possible auspicious directions, with dynamic factors and satellite data specifically mentioned for researchers interested in this area to embark on further explorations, are presented.
Key words: Moho     isostasy     gravity    

0 引 言

莫霍面是地壳与地幔的分界面.对莫霍面深度分布的研究可以获知区域范围内的地壳厚度变化和分布情况,对研究地球的深部构造特征和形成机制等问题有重要意义.

长期以来,研究人员通常采用地震与重测-均衡这两类方法研究莫霍面的深度分布,其中重测-均衡方法因为价格低廉、数据达到全球覆盖等原因而受到众多研究人员的青睐.但是在重测-均衡方法上,目前国内的许多研究对均衡模型的参考还远远不够,大多采用了非常简单的均衡模型,出现了均衡补偿模式研究与莫霍面反演研究在理论发展和实际工作中的不同步.考虑到当前对反演精度要求越来越高的实际情况,使用更为完善的均衡重力方法研究莫霍面深度分布是未来的趋势所在.

本文回顾了均衡补偿模式研究的发展情况及世界范围内莫霍面深度反演的现状,讨论了更先进的均衡理论引入莫霍面反演的必要性和可行性,并讨论了可能的发展方向.

1 均衡补偿模式发展情况

由流体静力平衡理论可知,在一定的深度下,单位面积的地壳块体应承载着相同的压力,称之为地壳均衡.众多地球物理资料已经证实了地壳均衡现象的存在,促使着越来越多的研究人员重视并研究均衡理论,发展均衡补偿模式.目前常用的均衡补偿模式包括局部补偿模式、区域补偿模式、试验补偿模式和动态补偿模式等,每种模式又包含各种不同的模型(Watts, 2001; 曾华霖, 2005; Bagherbandi, 2011a; Bagherbandi and Sj berg, 2012).

本文主要会对局部补偿模式中最常见的两种模型(Pratt-Hayford模型与Airy-Heiskanen模型)和区域补偿模式与全球补偿模式进行较详细的讨论,对其他均衡补偿模式的发展情况则进行简要说明,并厘清了均衡理论中一些令人困惑的问题.

1.1 局部补偿模式(Pratt-Hayford模型与Airy-Heiskanen模型)

Pratt-Hayford模型与Airy-Heiskanen模型是局部补偿模式中最典型的两种模型(曾华霖, 2005; Bagherbandi, 2011a).Pratt-Hayford模型选取的补偿深度D大约是113.7 km(李润江等, 1986),相当于把地壳划分为若干柱体,每个柱体具有不同的高度和密度,高度越高则密度越小,从而在补偿深度处达到静力平衡;Airy-Heiskanen模型的地壳正常厚度大约在30 km左右,认为宏观上地壳密度均一,高山地区存在“山根”,海洋地区有“反山根”,完成地壳均衡(曾华霖, 2005).由于补偿机制存在差异,二者也被分别称为密度补偿模式和深度补偿模式.由于使用简便,这两种补偿模式应用相当广泛(图 1图 2).

图 1 Pratt-Hayford模式示意图 Fig. 1 A schematic diagram of Pratt-Hayford model

图 2 Airy-Heiskanen模式示意图 Fig. 2 A schematic diagram of Airy-Heiskanen model

相比较之下,Pratt-Hayford模型补偿深度更深,更适合大区域深部构造研究(李润江等, 1986);Airy-Heiskanen模式从形式上符合地球物理原型.Airy-Heiskanen均衡模型的计算结果受地壳平均深度、壳幔密度差和模型剖分截面积的影响(赵玉合, 2007),其中壳幔密度差较地壳平均深度而言误差会对均衡补偿值造成的影响更为明显.另外,网格间距越小则保存细节信息越多,网格间距越大则高频信息越少,因此描述浅层地壳深度变化和均衡情况应该用小网格,描述深部地壳深度变化和均衡情况应该用大网格.

对于局部补偿模式,有两种常见的误解:一种是将Pratt-Hayford模型和Airy-Heiskanen模型对立起来,另一种是不能辩证看待局部补偿模式与实际情况的关系.对于前者,应注意到这两种模型都是基于平面地球近似的,而且由前人实验可知(殷秀华和刘铁胜等, 1993),这两类均衡补偿模型的结果相容,异常不存在空间和数值上的巨大变化;对于后者,应该指出,Pratt-Hayford模型和Airy-Heiskanen模型从各自的角度反映了地球内部真实的质量分布情况(曾融生, 1984).一方面,大量重力和地震等地球物理资料表明,地形越高的地区壳幔密度差越大,这一点与Pratt-Hayford模型假定相似;另一方面,地震测深资料说明一般情况下地形越高的地方地壳越厚,这和Airy-Heiskanen模型假定相符.实际的地壳均衡模式应为两者按一定比例结合的成果.Heiskanen(童永春, 1985)曾指出,从全球范围来看,地壳均衡可认为是由Airy-Heiskanen模型和Pratt-Hayford模型共同补偿完成的,前者大致占63%,后者约为37%.但Pratt-Hayford模型和Airy-Heiskanen模型各自并不完全符合实际情况.譬如,Pratt-Hayford模型假设整个补偿物质有唯一的质量中心,故推断出补偿水平面深度是一个常数.但实际上补偿物质质量中心有多个,不同地区的补偿水平面深度会差异很大.1970年时有研究表明最好的均衡模型需要在60 km、120 km、250 km和400 km深度处分别设定四个质量中心.另外,Airy-Heiskanen模型的地壳假定都是相同密度时在不同地区与实际情况会有差异.为此,20世纪60年代时伍拉德将可变地壳密度和上地幔密度模型用于Airy-Heiskanen均衡模型的修正,使之更符合地球物理的实际.一种可行的做法是把Pratt-Hayford模型和Airy-Heiskanen模型结合起来.前人(Bagherbandi, 2011a)曾经尝试在陆区使用Airy-Heiskanen模型,在海区选用Pratt-Hayford模型,取得了较好的效果.

局部补偿模式在构造时遵循了这样一个假设:地下任意点补偿密度直接决定其上地形高度,地面负载不影响临近块体.换言之,局部补偿模式在研究时只考虑了地幔的浮力作用,没有考虑地壳之间的相互作用,更没有考虑很老的岩浆岩侵入体对地壳均衡的影响.实际情况下,地壳存在内聚力,可保证柱体位置的稳定性.一般地,若地壳质量存在未补偿者,则会引起地壳内强烈的张性作用,促使地壳缓慢形变,形变到一定时候发生地壳断裂,也即地震.局部模式的另一个问题是,没有考虑到参与地形势场补偿的波长实际上是存在一定范围的,只是假定所有波长都可被补偿,这与实际情况并不相符(Bagherbandi, 2011a).此外,局部补偿模式无法解释部分地区“山越高地壳未必越厚”和“山越高岩石密度未必越低”的现象(张赤军, 2003).以上问题使得局部补偿模式在应用上受到很大的局限.

1.2 区域补偿模式

Vening Meinesz模型和全球补偿模式中的典型模型——Vening Meinesz-Moritz模型都属区域补偿模式.后者是在前者基础上的进一步发展.

1931年,Vening Meinesz提出了一种补偿模型(Meinesz, 1931).他在Airy-Heiskanen模型的基础上进行了修改,引入了“弹性板”模型的概念,认为地壳作为一个均质弹性板浮于粘性地幔之上,补偿质量在纵向和横向上都起到作用.在这个模型中,莫霍面上下地层的密度差恒定,莫霍 面深度可变(Bagherbandi, 2011a),而且这种“弹性板”模型较Hertz(1884)的传统模型(图 3)相比,考虑了间接效应(使用地壳密度填平了地表上外力作用点造成的作用点周围的地表凹陷部分).这种补偿模式考虑到了相邻地壳柱体之间的相互作用,指出的弹性弯曲效应涉及到周边,较局部补偿模式更符合实际情况,人们也将其称之为区域补偿模式(或弯曲模式).

图 3 Hertz模型(左)与考虑了间接效应的Vening Meinesz模型(右)的比较(Abd-Elmotaal,1993) Fig. 3 Comparison between Hertz Model (left) and Vening Meinesz Model (right) which considers the indirect effect (Abd-Elmotaal,1993)

图 4展示了Airy-Heiskanen模式与Vening Meinesz模式的对比,其中ρc为地壳密度,ρm为地幔密度,H为地形高度,t为Vening Meinesz补偿模式下的莫霍面深度,T0是正常地壳深度.图 3展示了Hertz模型和考虑了间接效应的Vening Meinesz模型,其中P为地表的力作用,f为最大弯曲,r为离负载起始位置的水平距离.显然,最大弯曲点在水平面上与受力点是重合的,水平方向上离受力点越远,板的弯曲度就更小,而描述这种板的弯曲度与受力点的水平距离关系的解是许多学者的研究重点.

图 4 Airy-Heiskanen模式与Vening Meinesz模式的比较(Bagherbandi et al., 2011) Fig. 4 Comparison between Airy-HeiskanenModel and Vening Meinesz Model (Bagherbandi et al., 2011)

对于上述问题,虽然Hertz(1884)曾给出过一个解(现在一般称为“传统解”),但该解没有考虑到间接效应,不适合地球物理直接使用.Meinesz(1931)在Hertz工作的基础上考虑了间接效应,后来众多学者参与公式简化的工作,指出板的弯曲方程解可由Kelvin函数表示(Abd-Elmotaal, 1993; 张赤军等, 2007),称为“精确解”.但Kelvin函数存在补偿半径无穷远处发散问题,而且计算起来过于复杂,为此Vening Meinesz考虑简化问题,将补偿半径限制为区域度的2.914倍,由此给出了“近似解”,也叫Vening Meinesz解.这种“近似解”避开了函数发散问题,且节省了机时,实验结果证明这种近似效果很好(Abd-Elmotaal, 1993;Suenkel,1985)基于调和分析,在线性估计条件下对比了Airy-Heiskanen模型和Vening Meinesz模型的补偿深度,发现如果Airy-Heiskanen模型使用更大的补偿深度,对Vening Meinesz模型进行泊松平滑从而补偿深度变浅,则此时二者等价; 方剑等(2006)基于区域补偿模式建立地形-均衡补偿重力的球谐关系式中,也得出了类似的结论,表明Airy-Heiskanen模型即Vening Meinesz模型在弹性刚度或弹性厚度为0时的特殊情况.以上工作将区域补偿模式与局部补偿模式联系了起来.

Vening Meinesz提出的区域补偿模式与局部补偿模式的相同点是都采用了平面近似.后来Moritz(1990)对Vening Meinesz的均衡模型进行了概括,引入了全球球面近似,将模型应用于球形海平面,使得区域补偿模式被推广为全球补偿模式,对应均衡模型被称为Vening Meinesz-Moritz补偿模型.有研究表明,Vening Meinesz-Moritz模型与地震方法得到的模型CRUST2.0吻合度较好(Bagherbandi, 2011a).

总体来看,区域补偿模式考虑到了地壳力学强度,较传统均衡模型更符合实际情况(黄方等, 2012).但应当注意,即使是当前区域补偿模式中比较先进的Vening Meinesz-Moritz补偿模型也存在局限性(Bagherbandi, 2011a).一方面,目前的Vening Meinesz-Moritz模型计算时还是使用了地球平均半径,也即将地球视为一个标准球体进行研究,如何将模型应用于椭球体上是下一步要解决的问题;另一方面,该模型假定地形势仅由地下质量来补偿,但事实上诸如构造运动和冰后期反弹等动力学因素也会起到补偿作用.为此,有的研究者考虑为此在计算均衡异常时除布格异常项和均衡改正项外,再增加一个校正项,包含板块构造运动、冰后期反弹、地幔对流与热补偿吸引以及其他影响.具体的做法还有待研究.

1.3 其他均衡模型

Parker(1973)曾撰文讨论如何利用傅立叶变换计算物质层在观测点上引发的重力改正值.文中Parker利用地形改正项进行了讨论,保证观测面为平面,高度在测区内地形最高点之上,且地形层可以扩展到顶底界面不平坦、层密度不均匀的情形.显然这种做法不局限于地形改正项的处理,可以用于构建均衡模型,称为Parker模型.Parker模型从本质上讲属于局部补偿模式(Bagherbandi, 2011a).与后来的Vening Meinesz-Moritz模式相比较,Parker的模型依然与Vening Meinesz一样使用平面近似,没有采用全球球面近似.Parker的模型后来由Oldenburg(1974)发展,形成了Parker-Oldenburg方法.这种方法在傅氏域进行计算,速度很快,常用于莫霍面深度估测,使用非常广泛.

实验均衡模式出现于1970年(Dorman and Lewis, 1970).这种模式认为补偿密度通过地形高和一个各向同性核函数的褶积与地形线性相关,故补偿质量的引力与地形和均衡响应函数的褶积相关,其中这种均衡响应函数可由反褶积获知,并且需要地形和补偿质量引力的调和分析(Abd-Elmotaal, 2005).因此他们转而利用实验数据直接算得均衡响应函数,以求得研究区内均衡改正值和密度展布.这种模式不需要对补偿机制、密度和深度进行假设.该模式提出后, Lewis and Dorman(1970)利用美国的重力和地形数据建立了均衡模型.后来,Banks采用了区域补偿思想对该理论进行了完善,所求均衡响应函数和地球内部密度梯度相关联(张赤军等, 2007).1990年,Moritz将此方法扩展到了球面域,并指出Pratt-Hayford模式可归结为实验均衡模式的一种特殊情况(张赤军等, 2007).但实验均衡模式的合理性没有得到公认.一方面,诸如美国等地实验均衡模式补偿深度过深的问题没有得到合理解释(张赤军等, 2007);另一方面,将补偿密度与地形高直接相关的做法有待商榷,因为目前已经知道地形高并不是影响补偿作用的唯一因素.

动态补偿模式是一种综合了传统补偿模式和实验均衡模式的方法,于20世纪80年代后期逐渐受人重视.这类模式具体包括弹性板模式和粘弹性板模式等,建模时考虑了区域补偿与局部补偿的异同以及地壳和上地幔的热结构变化及其对其他地球物理参量的影响.该模式计算精度较高,效果比较好.

1.4 均衡重力异常的正负与补偿程度的关系

虽然均衡作为一种理论已被多数地学工作者接受,但在对均衡的补偿机制和解释上不同的研究人员理解不一,导致长期以来很多文献中“正均衡异常”“负均衡异常”与“补偿不足”“补偿过剩”等术语之间关系的混乱,甚至导致同一地区不同的地球动力学图件完全相反(马杏垣, 1989; 袁学诚, 1996).曾华霖和万天丰(1999)研究了产生不同动力学解释结果的原因,澄清了“补偿”一词的概念,正确地将均衡重力异常的正负与补偿程度进行了关联.

首先是补偿、补偿不足与补偿过剩的概念.以Airy-Heiskanen模式为例,均衡重力研究中的“补偿”一词应指山下的质量亏损对地表山体质量的补偿.因此,莫霍面凹陷过深表明这种补偿过剩,莫霍面凹陷过浅表明这种补偿不足.

其次是正均衡异常与负均衡异常的概念.均衡重力异常的计算是在布格重力异常的基础上增加一个均衡校正项,该均衡校正项与上载地形和均衡面深度相关,可使研究区域均衡时均衡重力异常逼近0.因此,当实际莫霍面埋深较深时,均衡校正项不能“填满”质量亏损区域,导致均衡重力异常为负;反之,实际莫霍面埋深较浅时,均衡重力异常会变成正值(图 5).

图 5 均衡异常正负示意(曾华霖等,1999) Fig. 5 Demonstration of positive & negative isostatic anomalies (Zeng et al., 1999)

因此,正均衡异常应对应补偿不足,负均衡异常应对应补偿过剩.

1.5 全球范围内的地壳均衡

资料显示,地球上90%以上的地区处于地壳均衡状态(胡敏章和李建成, 2010).因此从宏观上看地球保持静力平衡.但是,因为非均衡在一些地方是事实存在的,故不能简单地认为地表地形和莫霍深度存在必然关系(Aitken, 2010).地球物理资料表明,全球范围内地壳不均衡的地区主要是板块运动剧烈或是火山活动频繁的地区(李姗姗等, 2004; 黄建平等, 2006),这一般是年轻的高山和深海沟地区.地壳不均衡的实质是重力等位面与等压面的不一致,产生的原因有多种,包括地表物质堆积和莫霍面抬升等(徐德琼等, 1983).均衡调整力的大小受地表负载水平宽度和岩石圈厚度的影响,一般而言构造运动先于均衡调整,故区域均衡重力异常主要反映了现代构造运动导致的均衡失调(殷秀华等,1982).

需要注意的是,岩石圈是具有弹性的,地壳均衡是在较大区域范围内实现的,因此相应的计算和研究也应该在区域范围内进行(杨宗仁, 1986;曾融生(1973)很早就指出,局部地形变化未必会得到补偿.王勇等(1996)曾对青藏高原进行计算,结果表明在较短波长范围内(小于50 km)的地质负荷几乎完全为地壳所支撑而非由局部补偿,这种情况下就不适宜采用均衡方法推估或内插邻近重力点.张赤军(2003)也指出,研究尺度不大的局部场时考虑补偿作用会导致场效应失真,地壳作为弹性体足以支撑小尺度重力场.

2 重测-均衡方法反演莫霍面深度
2.1 传统方法及尚需进行的改进

传统上对莫霍面深度的研究一般分为地震方法和重测-均衡方法两类(Bagherbandi, 2011b).地震方法非常常用,但存在着诸如分辨率太低、无法求解短波长特性和无法在三维域内直接提供地壳厚度分布等问题,而且价格昂贵,数据量太大,处理工作复杂.一般认为能提供更高分辨率和覆盖度的重力方法还是研究莫霍面深度的途径(许萍等, 2004; Aitken, 2010).

重测-均衡方法一般选择在地壳均衡较好的地区进行,使用布格重力异常数据,加以区域异常分离处理去除额外影响,只保留莫霍面对重力的贡献,并使用各种莫霍面深度计算方法反演深度(高玉文等, 2012; 郭良辉等, 2012; 姜文亮等, 2012).之所以使用布格重力异常,一方面是因为大尺度上布格重力异常值与莫霍面深度分布存在一定相关性,另一方面是因为布格重力数据包含了整个岩石圈的信息,可反映岩石圈横向和纵向密度不均匀性(Jiang and Zhang et al., 2011).这方面国内外的许多学者做了不少工作(许萍等, 2004).刘元龙等(1987)采用三维重力方法研究地壳构造,得到了我国部分区域的地壳厚度和莫霍面的分布情况; Molnar(1988)利用印度河-恒河平原与喜马拉雅山脉在布格异常梯度上的差异对印度板块的弯折进行了解释.马杏垣(1989)根据中国大陆的布格重力异常绘制了莫霍界面图; 方剑(1999)利用卫星重力数据,反演并绘制了全球范围内的地壳厚度分布图; 江为为等(2003)和宋学锋(2011)先后使用布格重力异常对冲绳海槽地区莫霍面分布进行了反演并给出了地质解释.

但单纯使用传统重力方法也存在一定问题,因为诸如重力和地形的频谱分析等常用的莫霍面建模方法一般不能明确地施以地震估测限制,限制了反演的准确性,而且使用的模型往往过于简化,诸如水平方向同质的设定不符合实际情况(Aitken, 2010).鉴于此,研究人员一般从以下两种途径入手考察解决方案.

第一种途径,也是目前的研究趋势,是地震与重力资料结合探究莫霍面深度(秦静欣等, 2009).这一般有两种情况:一种是用二维重力数据给地震结果增加限制,另一种是由地震数据得到莫霍面深度构建三维重力模型(Szafián and Horváth, 2006).已经有很多三维模型软件包可供研究人员使用(Gtze and Lahmeyer, 1988;Lefort and Agarwal(2002)研究对比了重力方法与地震方法得到的莫霍面深度信息,指出二者可能存在差异,不过适当修改质量分布计算的参考水平面和密度差异就可以使二者拟合.另外,当莫霍面深度变化较小时,用于计算莫霍面起伏的滤波技术结果与地震折射数据都非常好,但莫霍面地形变化大时,绝对莫霍面深度会和地震信息有差异.地壳变薄时,莫霍扰动会减小.胡敏章和李建成(2010)尝试在重力建模时引入CRUST2.0 提供的地球物理信息,使用真实莫霍面深度代替Airy-Heiskanen模式假设深度,并与传统Airy-Heiskanen模型进行了对比计算,证实了此举的可行性,并证明这种替换会在均衡重力异常中体现出更多的细节信息.Aitken(2010)在反演澳大利亚的莫霍面时摒弃了地形与莫霍面几何特征关系的假定,使用了基于空间重力异常限制反演的新方法制作大陆尺度的莫霍面模型,一方面可充分利用已有地震资料来限制反演,空间域的重力数据三维反演保证了重力数据的高精度和一般覆盖.Eshagh、Sjberg和Bagherbandi(2011a)对结合地震和重测数据的全球莫霍面模型也进行过研究.

第二种途径则是从采用更合适的地壳均衡模型的角度出发来考虑,这也是本文论述的重点.

2.2 均衡理论对传统反演方法的贡献

均衡理论能科学地反映大尺度情形下地壳的内部结构,较好地阐释地壳和地幔的相互关系,而莫霍面又是地壳与地幔的分界,因此在莫霍面深度反演中更好地引入均衡理论、参考利用或修正各种均衡模型是十分必要的.这一方面,具体的做法有以下几类.

2.2.1 完善传统均衡模型

这里的传统均衡模型主要指Airy-Heiskanen模型.因该模型形式上符合地球物理原型,使用比较方便,一些研究人员考虑直接在该模型基础上进行修正完善.

具体方法包括对一些参数的修正及考虑其他因素等.譬如,莫霍面反演中常用到压缩质面迭代法,该方法需要提供一个参考深度D,通过计算得到莫霍界面起伏值H,加在参考深度之上即得到莫霍面的深度.由于在同一区域内反演出的H基本不随所选的D的变化而变化,D的精确度对莫霍面反演深度的精确性影响很大.黄建平等(2006)在对中国及其邻区的莫霍面深度反演研究中采用了Airy-Heiskanen均衡模型,通过分析莫霍面起伏与地表地形数据的相关性,发现陆地的均衡系数接近于常用值4.45,但海洋上的计算结果与理论值2.73相差较大,由此给出了修正后的均衡系数值,并给出了全地区、大陆地区和海洋地区的参考深度与高程的关系式.但作者也指出该方法在地壳不均衡地区反演结果存在误差,且研究结论不一定适用于全球.胡敏章和李建成(2010)在研究Airy-Heiskanen模型时,考虑了沉积层和地壳内部载荷对均衡重力场的影响(图 6图 7),使得完善后的模型更贴近实际情况.图 6右图中水平实线为零异常线,实线a为Airy模式下莫霍面对沉积层响应导致的重力异常,虚线b为沉积层导致的重力异常(未考虑补偿作用),点划线c为同时考虑沉积层和均衡抵偿后的重力异常.图 7中密度异常体200×200 km,剩余密度为0.1 g/cm3.

图 6 沉积层展布(左)及其引发的重力效应(右)(胡敏章等,2010) Fig. 6 Distribution of the sediment layer (left) and its gravity effect (right) (Hu M Z et al., 2010)

图 7 地壳内部载荷及其引发重力效应(胡敏章等,2010) Fig. 7 Load in the inner crust and its gravity effect (Hu M Z et al., 2010)

另一种做法利用了Airy-Heiskanen模型与Vening Meinesz模型的联系,试图使Airy-Heiskanen模型解逼近Vening Meinesz模型解.Abd-Elmotaal(1993)在对中欧地区莫霍面分布进行研究时指出,给Airy-Heiskanen莫霍深度施加一个光滑因子,得到的结果与Vening Meinesz解很类似,二者之差几乎可忽略,且计算速度较Vening Meinesz解快很多.

2.2.2 应用更新的均衡模型

该方面的一个典型案例是Parker-Oldenburg方法的广泛应用.Parker模型出现之后, Oldenburg(1974)将其应用于重力异常的反演与解释,提出了一种以迭代方式反演地下单一密度界面的方法,常用于莫霍面深度分布的研究.Gómez-Ortiz and Agarwal(2005)Shin et al.(2006)基于Parker-Oldenburg方法,采用快速傅里叶变换技术估计了法国的Brittany和韩国的Ulleung盆地的莫霍面深度.Shin et al.(2006)基于Parker-Oldenburg方法,使用GRACE卫星的数据对西藏地区莫霍面起伏情况进行了研究.使用GRACE卫星数据,基于Parker-Oldenburg方法考察了南极洲地壳厚度(Bagherbandi, 2011b).

区域补偿模式和全球补偿模式在莫霍面深度分布或地壳厚度研究方面也得到了应用.Moritz(1990)建立了全球补偿模式后,给出了一种迭代求解莫霍面深度分布的方法.Abd-Elmotaal(1993)利用Vening Meinesz模型研讨了莫霍面深度分布问题,得到了传统解、精确解和估计解.Sjberg(2009)使用Vening Meinesz-Moritz模型,通过求解第一类非线性Fredholm方程将公式线性化,指出了基于重测数据的近似和实用的地壳厚度估计方法.

国外研究人员在使用更新的补偿模式的同时,也注重不同补偿模式对结果的影响对比.Bagherbandi(2011b)在研究青藏高原地壳厚度的时候比较了基于不同补偿模式的三种方法:基于Vening Meinesz-Moritz模型的Vening Meinesz-Moritz迭代法及Sjberg直接求解法,以及基于Parker模型的Parker-Oldenburg方法,并将反演结果与CRUST2.0进行了对比(图 8).结果指出,Sjberg直接求解法无需迭代和正则化,反演结果与CRUST2.0的数据偏差小,方差也小,反演效果要比Parker-Oldenburg方法和Vening Meinesz-Moritz迭代法效果好,另外,基于Vening Meinesz-Moritz的模型反演效果比Parker-Oldenburg的模型反演效果好,但这两种模型在迭代时都存在发散现象,这一方面是因为存在数据误差,另一方面是处理中要求把非线性反演问题强制线性化.但Bagherbandi的工作也有局限性,比如采用了地壳均匀假设、忽略了低密度沉积层、研究莫霍密度差时忽视了地幔密度的横向变化、未考虑热均衡等.

图 8 四种莫霍深度解的比较(Bagherbandi, 2011) (a) Sjberg直接求解法; (b)Parker-Oldenburg法; (c)Vening Meinesz-Moritz迭代法; (d)CRUST2.0解. Fig. 8 Comparison among four solutions in relation to Moho depth (Bagherbandi, 2011) (a) Sjberg’s direct solution; (b)Parker-Oldenburg solution; (c)iterative solution based on Vening Meinesz-Moritz method; (d)CRUST2.0 solution.
2.2.3 考虑模型地壳和地幔的密度不均匀性

应当注意,即使是比较新的均衡模型,如Vening Meinesz-Moritz模型等,都蕴含了壳幔密度差恒定的假定.国外方面,早在2002年, Lefort and Agarwal(2002)通过重力数据对法国的莫霍面起伏地形进行研究时就意识到,大范围内不应使用统一的壳幔密度差.Aitken(2010)在反演澳大利亚的莫霍面时采用了地壳与地幔分别满足水平方向异性的模型.Bagherbandi(2011a)指出,使用传统的0.6 g/cm3的壳幔密度差在一些地区(尤其是海区)是过大的.因此他在研究Vening Meinesz-Moritz模式时,尝试在海区和陆区使用不同的壳幔密度差,反演了全球的莫霍深度,并与CRUST2.0结果进行了比较,结果表明这种处理的反演效果更好(表 1)(Bagherbandi and Sjberg, 2011).国内方面, 许萍等(2004)考虑到了地壳和上地幔密度不均匀这一影响因素,对布格重力异常进行了修正,比较结果证实这一影响因素是不可忽略的.

一个案例是相对于标准地壳的微分密度模型的使用.传统的研究方法是建立参考地壳(或岩石圈),厚度是常数,每个单元使用平均密度,假定的要被模型化的结构密度从参考值中减去,重力异常由密度差得到.如果想达到足够好的效果,需要研究者获取研究区域里足够多的结构、地震和重力数据(Holliger and Kissling, 1992),否则只能通过诸如忽视密度纵向变化等方法来简化模型,这样得到的密度差和重力异常值很可能是不准确的.更好的方法是使用相对于标准地壳的微分密度模型.这种情形下,不同的边界仅由密度差异表征,允许出现密度随深度的变化.Lillie et al.(1994)对此曾进行过研究, Szafián and Horváth(2006)也基于该方法对Carpatho-Pannonian地区进行研究.

表 1 CRUST2.0与Vening Meinesz-Moritz模型(常密度/变密度)莫霍深度统计比较表(Bagherbandi and Sj oberg , 2011) Table 1 Statistics of Moho depths for models CRUST2.0 and Vening Meinesz-Moritz (Bagherbandi and Sj oberg , 2011)
2.2.4 对动力学补偿作用的研究

上文曾以Vening Meinesz-Moritz模型为例指出考虑动力学因素补偿作用的重要性.动力学补偿作用一般分为四项:板块构造运动Atm、冰后期反弹Apgr、地幔对流与热补偿吸引Amc以及其他影响Aε.板块构造运动指岩石圈的大尺度移动;冰后期反弹考虑了更新世冰盖融化的作用,这在北美和北欧地区是必须考虑的,目前已经有了末次冰期冰川均衡调整(GIA)模型(汪汉胜等, 2012); 地幔对流对地形势可以 起到一定的补偿作用;其他影响综合了其他动力学补偿作用,也包括原始数据误差造成的影响.这四项改正之和称为动力均衡改正项ΔgC,使用公式可表示为

ΔgC=Atm+Apgr+Amc+Aε .

有研究表明,在采用重力-均衡方法得到的均衡异常中增加这一附加项,得到的均衡异常会更有说服力(Bagherbandi, 2011a).但从总体上看,对这几个改正项的研究目前仍处于起步阶段,具体的研究方法还有待探索.相比较而言,这四项中目前比较热门的课题是第三项,即热均衡研究(Hasterok and Chapman, 2007a; Hasterok and Chapman, 2007b; Chen and Zhang et al., 2008).

热均衡理论认为,岩石圈拉张减薄使得热物质上涌,由此地壳会热胀冷缩,这种变化可由地表热流数据得到,研究人员可以借此计算岩石圈热结构及热均衡模式.大陆热均衡研究理论与方法由Hasterokl和Chapman(Hasterok and Chapman, 2007a; Hasterok and Chapman, 2007b)完善,并在北美进行了应用研究.按照该理论,热均衡问题可通过地表高程和热流的相关关系进行研究.同时,他们提出的均衡分析技术可很好地校正密度差异导致的均衡异常量,恢复在热均衡过程中地温梯度变化对区域高程变化产生的相关关系.具体做法可大致描述如下(陈石等, 2011):

首先,给出温度与深度的关系式,用于构造地壳分层理论参考模型.这里的关系式不是唯一的,考虑到地壳不同深度的热机制差异,岩石圈上部和下部应给出不同的关系式.譬如,认为岩石圈上部以热传导机制为主,下部热对流作用明显,则可分别给出公式

其中T表示温度(单位为K),q为热流,Δz为层位厚度,k为热导率,A为生热率,下标i为层位号(对于T和q,下标表示对应层号底部的相应物理量值),z为深度,Ta取1573K 此时 T z 为0.3 K/km .对于热导率k,与温度、压力等有关,可采用经验公式(McKenzieet al., 2005)

k(T)= 5.3 1+0.0015T +∑ 3 m=0 dm(T+273)m ,

这里d0=1.753×10-2,d1=-1.0365×10-4,d2=2.2451×10-7,d3=-3.4071×10-11.

其次,利用均衡分析技术校正观测高程ε,排除温度造成的地形缩放效应ΔεT和水平密度不均匀性造成的均衡调整效应Δερ,即

ε=εobsΔεTΔερ ,

其中

ΔεT=α∫zmax0[T(z)-Tref(z)]dz ,

εobs为观测高程,α为岩石热膨胀系数,zmax为岩石圈最大深度,T(z)是研究区地温,Tref(z)是标准参考地温结构.其中,

Δερ=hn 1- ρnc ρm -h 1- ρc ρm ,

这里h为厚度,ρ为密度,下标n表示参考模型,c为地壳,m为地幔.

经过上述修正后再进行均衡分析,得到的结果就考虑到了热均衡的影响.

陈石等(2011)使用该理论计算了组分均衡调节量,在华北地区的应用结果表明该方法可以使算得的莫霍面深度更为精确,也可看出整个均衡重力异常与上地幔密度分布较一致,梯度较大的区域莫霍面的变形也比较剧烈,恰好是华北地区各个活动块体单元边界.Jiménez-Munt et al.(2012)在研究欧亚碰撞带的地壳和岩石圈地幔几何一阶估计时引入了热分析的概念,在分离区域和局部异常时使用了拟合高程和融合了热分析的大地水准面异常数据的方法,并将分析结果与常见的三种研究方法(滤波法、地壳均衡、岩石圈均衡)进行了对比,指出前者最符合研究区地质特征.

2.2.5 对地壳不均衡地区的研究

如何得到地壳不均衡地区的莫霍面深度一直是研究难点.一些研究者考虑研究地壳不均衡地区重力异常值与其他影响因素的关系.许多研究人员在研究中意识到,地壳不均衡的地区空间重力异常与扰动地形产生的重力影响紧密相关,因此可通过计算地形重力效应和空间重力异常的相关频谱来估计空间重力异常(von Frese et al., 1997; Roman, 1999; von Frese et al., 1999; van der Veen et al., 2007).地形影响一旦减去,将异常上延到一定高度,异常值就可以反映莫霍面定义的密度分界影响,由此得到的重力异常值也可以用来与地震莫霍面估计值进行比较或加以限制(Dahl-Jensen et al., 2003).这种方法不受具体的均衡模式的限制,但假定了与地形相关的空间异常必然反映地形效应,同时要求对地形密度和壳幔密度差等进行假定,这可能会引发误差,不过在上延高度足够高的情况下误差可忽略.

一个案例是对格陵兰地区的地壳厚度和莫霍面深度进行的研究(Braun et al., 2007).格陵兰地区因露头和地震观测数据的稀缺,冰下地壳先验信息匮乏.研究人员对地形和空气重力异常进行了建模,由对均衡校正后的地形影响的反演得到了莫霍面和相应地壳厚度的估计,取得了比较好的效果(图 9,图中点位为主要冰核钻井点位).

图 9 格陵兰岛地壳厚度反演结果(据Braun,2007有修改) Fig. 9 The inversion result of crustal thickness in Greenland (Braun et al., 2007, edited)
3 讨 论

可以看出,近年来在国际上重测-均衡方法在理论上有了很大发展,并已被广泛应用于莫霍面深度分布的反演效果改善.这种应用既包括对传统均衡模型的进一步完善,也包括新均衡模型的投入使用.对于国内的研究人员而言,这也是在探索提高重测-均衡方法反演莫霍面深度分布效果方面可以选择的两种思路.

从模型的角度考虑,整体的趋势是模型从简单到复杂,从常密度模型到变密度模型,逐渐向实际情况靠拢.上文已述及,目前还有很多均衡模式没有应用于莫霍面深度反演中,已经投入使用的均衡模型也需要进一步改进.

对动力学补偿作用的研究将成为改善均衡模型的重要研究方向.如上文所述,为了提高反演的精度,动力均衡因素是应当考虑到的.实际工作中可参考具体研究区域来权衡各个动力因素对反演结果的影响,比如在北美和北欧等更新世为大冰盖所覆盖的地区,冰后期反弹因素很大,动力学作用补偿中必须考虑冰后期反弹项;地幔对流作用强烈的地区,必须考虑热均衡因素.

另外要重视卫星数据的重要性(Bagherbandi, 2011a; 林淼等, 2012).重测卫星提供了全球范围内的重力数据覆盖.研究重力场时常用GOCE卫星重力梯度测量(SGG)数据.GOCE数据的重力异常精度可达1mGal,提供地球重力场球谐系数达300阶,且SGG数据估计莫霍面深度不会存在截断或正则化过程导致的误差.近年来在国外有不少研究人员从事SGG数据下延与莫霍深度估计的结合工作,这可能将成为接下来的一个研究热点.

致 谢 感谢青岛海洋地质研究所的雷受旻教授,他为本文的完成提供了大量有益的建议与指导.中国科学院地质与地球物理研究所的张丽莉副研究员、黄松博士后、姬莉莉博士后、胡卫剑博士、刘丽华博士和胡立天博士对本文写作也提出了中肯的建议,在此一并致谢.同时感谢审稿专家对本文给出的意见和建议!

参考文献
[1] Abd-Elmotaal H. 1993. Vening Meinesz Moho depths: traditional, exact and approximated[J]. Manuscripta Geodaetica, 18(4): 171-181.
[2] Abd-Elmotaal H. 2005. Global isotropic isostatic response function in terms of zonal spherical harmonics[J]. A Window on the Future of Geodesy, 128: 471-476,   doi: 10.1007/3-540-27432-4_80.
[3] Aitken A R A. 2010. Moho geometry gravity inversion experiment (MoGGIE): A refined model of the Australian Moho, and its tectonic and isostatic implications[J]. Earth and Planetary Science Letters, 297(1-2): 71-83,   doi: 10.1016/j. epsl. 2010. 06.004.
[4] Bagherbandi M. 2011a. An Isostatic Earth Crustal Model: and Its Applications[D].   Stockholm, Sweden: Royal Institute of Technology.
[5] Bagherbandi M. 2011b. A comparison of three gravity inversion methods for crustal thickness modelling in Tibet plateau[J]. Journal of Asian Earth Sciences, 43(1): 89-97, doi: 10.1016/j. jseaes.2011.08.013.
[6] Bagherbandi M, Sjberg L E. 2011. Comparison of crustal thickness from two gravimetric-isostatic models and CRUST2.0[J]. Studia Geophysica et Geodaetica, 55(4): 641-666,   doi: 10.1007/s11200-010-9030-0.
[7] Bagherbandi M, Sjberg L E. 2012. Modelling the density contrast and depth of the Moho discontinuity by seismic and gravimetric-isostatic methods with an application to Africa[J]. Journal of African Earth Sciences, 68: 111-120,   doi: 10.1016/j. jafrearsci.2012.04.003.
[8] Braun A, Kim H R, et al. 2007. Gravity-inferred crustal thickness of Greenland[J]. Earth and Planetary Science Letters, 262(1-2): 138-158,   doi: 10.1016/j. epsl. 2007.07.050.
[9] Chen S, Zhang J, et al. 2008. Gravity inversion using the frequency characteristics of the density distribution[J].   Applied Geophysics, 5(2): 99-106.
[10] Chen S, Wang Q S, Xu W M, et al. 2011. Thermal isostasy of North China and its gravity isostasy and deep structure[J]. Chinese J. Geophys. (in Chinese), 54(11): 2864-2875,   doi: 10.3969/j. issn.0001-5733.2011.11.016.
[11] Dahl-Jensen T, Larsen T B, et al. 2003. Depth to Moho in Greenland: receiver-function analysis suggests two Proterozoic blocks in Greenland[J]. Earth and Planetary Science Letters, 205(3): 379-393,    doi: 10.1016/S0012-821X(02)01080-4.
[12] Dorman L M, Lewis B T R. 1970. Experimental isostasy 1. Theory of the determination of the earth’s isostatic response to a concentrated load[J]. Journal of Geophysical Research, 75(17): 3357-3365,    doi: 10.1029/JB075i017p03357.
[13] Fang J. 1999. Global crustal and lithospheric thickness inversed by using satellite gravity data[J]. Crustal Deformation and Earthquake, 19(1): 29-34,   doi: 10.3969/j. issn.1671-5942.1999.01.004.
[14] Fang J, Ma Z J, Xu Z H, et al. 2006. Frequency analyze of topography-isostasy compensative gravity and geoid anomaly[J]. Progress in Geophys.   (in Chinese), 21(1): 25-30.
[15] Gao Y W, Luo Y, Wen W. 2012. The compensation method for downward continuation of potential field from horizontal plane and its application[J]. Chinese J. Geophys. (in Chinese), 55(8): 2747,    doi: 10.6038/j. issn. 0001-5733. 2012. 08.026.
[16] Gómez-Ortiz D, Agarwal B N P. 2005. 3DINVER. M: a MATLAB program to invert the gravity anomaly over a 3D horizontal density interface by Parker-Oldenburg’s algorithm[J]. Computers & Geosciences, 31(4): 513-520,   doi: 10.1016/j. cageo. 2004.11.004.
[17] Gtze H J, Lahmeyer B. 1988. Application of three-dimensional interactive modeling in gravity and magnetics[J]. Geophysics, 53(8): 1096-1108, doi: 10. 1190/1.   1442546.
[18] Guo L H, Meng X H, Shi L, et al. 2012. Preferential filtering method and its application to Bouguer gravity anomaly of Chinese continent[J]. Chinese J. Geophys. (in Chinese), 55(12): 4078-4088,   doi: 10.6038/j. issn. 0001-5733.2012.12.020.
[19] Hasterok D, Chapman D S. 2007a. Continental thermal isostasy: 1. Methods and sensitivity[J]. Journal of Geophysical Research-Solid Earth, 112(B6): 6414, doi: 10.  1029/2006JB004663.
[20] Hasterok D, Chapman D S. 2007b. Continental thermal isostasy: 2. Application to North America[J]. Journal of Geophysical Research-Solid Earth, 112(B6): 6415, doi: 10.  1029/2006JB004664.
[21] Hertz H. 1884. über das Gleichgewicht schwimmender elastischer Platten[J].   Annalen der Physik, 258(7): 449-455.
[22] Holliger K, Kissling E. 1992. Gravity interpretation of a united 2-D acoustic image of the central Alpine collision zone[J]. Geophys J. Int., 111(2): 213-225,   doi: 10.1111/j. 1365-246X. 1992. tb00571.x.
[23] Hu M Z, Li J C. 2010. Computation of Airy-Heiskanen isostatic gravity anomaly with considering crust density model[J]. Journal of Geodesy and Geodynamics, 30(5): 48-52.
[24] Huang F, Liu Q Y, He L J. 2012. Tectono-thermal modeling of the Sichuan Basin since the Late Himalayan period[J]. Chinese J. Geophys. (in Chinese), 55(11): 3742-3753,   doi: 10.6038/j. issn. 0001-5733. 2012.11.021.
[25] Huang J P, Fu R S, Xu P, et al. 2006. Inversion of gravity and topography data for the crust thickness of China and its adjacency[J]. Acta Seismologica Sinica, 28(3): 250-258,   doi: 10.3321/j. issn: 0253-3782. 2006.03.004.
[26] Jiang W W, Liu S H, Hao T Y, et al. 2003. Character of geophysical field and crustal structure of Okinawa Trough and adjacent region[J]. Progress in Geophys. (in Chinese), 18(2): 287-292,    doi: 10. 3969/j. issn. 1004-2903. 2003. 02.017.
[27] Jiang W L, Zhang J F. 2012. Fine crustal structure beneath Capital area of China derived from gravity[J]. Chinese J. Geophys. (in Chinese), 55(5): 1646-1661,   doi: 10.6038/j. issn. 0001-5733. 2012.05.022.
[28] Jiang W L, Zhang J F, Lu X C, et al. 2011. Crustal structure beneath Beijing and its surrounding regions derived from gravity data[J]. Earthquake Science, 24(3): 299-310, doi: 10.  1007/s11589-011-0792-4.
[29] Jiménez-Munt I, Fernàndez M, Saura E, et al. 2012. 3-D lithospheric structure and regional/residual Bouguer anomalies in the Arabia-Eurasia collision (Iran)[J]. Geophysical Journal International, 190(3): 1311-1324,    doi: 10. 1111/j. 1365-246X. 2012. 05580.x.
[30] Li R J, Tong S J, Meng L S, et al. 1986. Preliminary study on the state of isostatic equilibrium in the northeast China[J]. Seismology and Geology, 8(3): 15-19.
[31] Li S S, Wu X P, Li J W. 2004. Inverse of the Mohorovicic’s surface based on Airy-Heiskanen’s isostatic model[J]. Science of Surveying and Mapping, 29(2): 25-27, doi: 10.3771/j. issn.1009-2307.2004.02.008.
[32] Lin M, Zhu J J, Tian Y M, et al. 2012. On the use of the geoid anomalies for geophysical interpretation over the area of Hunan[J]. Chinese J. Geophys. (in Chinese), 55(2): 472-483,   doi: 10.3969/j. issn. 0001-5733.2012.02.011.
[33] Liu Y L, Zheng J C, Wu C Z. 1987. Inversion of the three dimensional density discontinuity by use of a method of “compressed mass plane coefficient” based on the gravity data[J]. Acta Geophysica Sinica, 30(2): 186-196, doi: 10. 3321/j. issn: 0001-5733.1987.02.  009.
[34] Lefort J P, Agarwal B N P. 2002. Topography of the Moho undulations in France from gravity data: their age and origin[J]. Tectonophysics, 350(3): 193-213, doi: 10.  1016/S0040-1951(02)00114-2.
[35] Lewis B T R, Dorman L R M. 1970. Experimental isostasy 2. An isostatic model for the USA derived from gravity and topographic data[J]. Journal of Geophysical Research, 75(17): 3367-3386, doi: 10.  1029/JB075i017p03367.
[36] Lillie R J, Bielik M, Babuka V, et al. 1994. Gravity modelling of the lithosphere in the Eastern Alpine-Western Carpathian-Pannonian Basin region[J]. Tectonophysics, 231(4): 215-235, doi: 10.  1016/0040-1951(94)90036-1.
[37] Ma X Y. 1989. Lithospheric Dynamics Atlas of China[M]. Beijing: SinoMaps Press.
[38] McKenzie D, Jackson J, Priestley K. 2005. Thermal structure of oceanic and continental lithosphere[J]. Earth and Planetary Science Letters, 233(3): 337-349,    doi: 10.1016/j. epsl. 2005. 02.005.
[39] Meinesz F A V. 1931. Une nouvelle methode pour la reduction isostatique regionale de l’intensite de la pesanteur[J]. Bulletin Géodésique (1922-1941), 29(1): 33-51, doi: 10.1007/BF03030038.
[40] Molnar P. 1988. A review of geophysical constraints on the deep structure of the Tibetan Plateau, the Himalaya and the Karakoram, and their tectonic implications[J]. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 326(1589): 33-88,    doi: 10. 1098/rsta. 1988.0080.
[41] Moritz H. 1990. The Figure of the Earth: Theoretical Geodesy and the Earth’s Interior[M]. Karlsruhe, Germany: Wichmann.
[42] Oldenburg D W. 1974. The inversion and interpretation of gravity anomalies[J]. Geophysics, 39(4): 526-536, doi: 10.1190/1.   1440444.
[43] Parker R L. 1973. The rapid calculation of potential anomalies[J]. Geophysical Journal of the Royal Astronomical Society, 31(4): 447-455,    doi: 10.1111/j.1365-246X.1973.tb06513.x.
[44] Qin J X, Hao T Y, Xu Y, et al. 2009. The consideration of the Moho depth of China and adjacent areas[J]. Progress in Geophys. (in Chinese), 24(5): 1593-1601,   doi: 10.3969/j. issn.1004-2903.2009.05.007.
[45] Roman D R. 1999. An Integrated Geophysical Investigation of Greenland’s Tectonic History [Ph. D. thesis]. Ohio: The Ohio State University.
[46] Shin Y H, Choi K S, Sun K, et al. 2006. Three-dimensional forward and inverse models for gravity fields based on the Fast Fourier Transform[J]. Computers & Geosciences, 32(6): 727-738,   doi: 10.1016/j. cageo.2005.10.002.
[47] Sjberg L E. 2009. Solving Vening Meinesz‐Moritz inverse problem in isostasy[J]. Geophysical Journal International, 179(3): 1527-1536,    doi: 10.1111/j. 1365-246X. 2009. 04397.x.
[48] Song X F. 2011. The fault structure interpretation and Moho surface inversion in Okinawa trough[J]. Petroleum Geophysics, 9(2): 48-50.
[49] Suenkel H. 1985. An isostatic earth model[R]. Ohio, United States: DTIC Document.
[50] Szafián P, Horváth F. 2006. Crustal structure in the Carpatho-Pannonian region: insights from three-dimensional gravity modelling and their geodynamic significance[J]. International Journal of Earth Sciences, 95(1): 50-67, doi: 10.   1007/s00531-005-0488-x.
[51] Tong Y C. 1985. The initial exploration of Moho discontinuity and crustal isostatic characteristics in Jiangsu[J]. Geology and Prospecting, 21(9): 37-41.
[52] Van der Veen C J, Leftwich T, et al. 2007. Subglacial topography and geothermal heat flux: Potential interactions with drainage of the Greenland ice sheet[J]. Geophysical Research Letters, 34(12): L12501, doi: 10.   1029/2007GL030046.
[53] Von Frese R R B, Tan L, Potts L V, et al. 1997. Lunar crustal analysis of Mare Orientale from topographic and gravity correlations[J]. Journal of Geophysical Research, 102(E11): 25657-25675, doi: 10.  1029/97JE02216.
[54] Von Frese R R B, Tan L, Kim J W, et al. 1999. Antarctic crustal modeling from the spectral correlation of free-air gravity anomalies with the terrain[J]. Journal of Geophysical Research, 104(B11): 25275-25296, doi: 10.   1029/1999JB900232.
[55] Watts A B. 2001. Isostasy and Flexure of the Lithosphere[M].   Cambridge, United Kingdom: Cambridge Univ Pr.
[56] Wang H S, Jia L L, Wu P, et al. 2012. Effects of last-deglaciation on the historical relative sea levels of East Asia Seas and the implications[J]. Chinese J. Geophys. (in Chinese), 55(4): 1144-1153,    doi: 10.6038/j. issn.0001-5733.2012. 04.010.
[57] Wang Y, Xu H Z. 1996. The variations of lithospheric flexural strength and isostatic compensation mechanisms beneath the continent of China and vicinity[J]. Acta Geophysica Sinica, 39(S): 105-113.
[58] Xu D Q, Li Q X, Jiang J Z. 1983. The Mohorovicic discontinuity in the east China sea and its geologic significance[J]. Marine Science Bulletin, 2(5): 34-41.
[59] Yang Z R. 1986. Quantitative study of gravitational isostasy--taking Hu’nan province region for example[J]. Crustal Deformation and Earthquake, 6(4): 285-291.
[60] Yin X H, Liu T S. 1993. Isostatic anomalies of gravity and geologic structures of the surface and shallow crust[J].   Seismology and Geology, 15(2): 149-156.
[61] Yin X H, Shi Z H, Liu Z P, et al. 1982. Preliminary study on the isostatic gravitational anomaly in the north of northern China[J].   Seismology and Geology, 4(4): 27-34.
[62] Yuan X C. 1996. Atlas of Geophysics in China[M]. Beijing: Geological Publishing House.
[63] Zeng H L. 2005. Gravity Field and Gravity Exploration[M]. Beijing: Geological Publishing House.
[64] Zeng H L, Wan T F. 1999. Difference between two isostatic gravity anomaly maps of China[J]. Chinese J. Geophys. (in Chinese), 42(1): 127-134, doi: 10. 3321/j. issn: 0001-5733. 1999. 01. 015.
[65] Zeng R S. 1973. Gravity compensation of the Mohorovicic discontinuity and the basic model of crustal structure[J]. Chinese J. Geophys.   (in Chinese), 16(1): 1-5.
[66] Zeng R S. 1984. Introduction of Solid Geophysics[M]. Beijing: Science Press.
[67] Zhang C J. 2003. The gravity response in crust on different wavelength topography[J]. Progress in Geophysics (in Chinese), 18(2): 342-347,    doi: 10. 3969/j. issn. 1004-2903. 2003. 02.027.
[68] Zhang C J, Fang J, Bian S F, et al. 2007. Gravity (Crust) isostasy and topographic wavelength--on isostatic isn’t suitable for gravity estimation on Qumolangma Peak[J]. Geomatics and Information Science of Wuhan University, 32(3): 229-233, doi: 10.3969/j. issn.1671-8860.2007.03.011.
[69] Zhao Y H. 2007. Calculation and synthetical interpretation of isostatic gravity anomalies [M. S. Thesis]. Beijing: Institute of Geology and Geophysics Chinese Academy of Sciences.
[70] 陈石, 王谦身, 徐伟民,等. 2011. 华北地区热均衡, 重力均衡与深部构造[J]. 地球物理学报, 54(11): 2864-2875,    doi: 10. 3969/j. issn. 0001-5733. 2011. 11.016.
[71] 方剑. 1999. 利用卫星重力资料反演地壳及岩石圈厚度[J]. 地壳形变与地震, 19(1): 29-34,    doi: 10. 3969/j. issn. 1671-5942. 1999. 01.004.
[72] 方剑, 马宗晋, 许厚泽,等. 2006. 地形-均衡补偿重力、大地水准面异常频谱分析[J].   地球物理学进展, 21(1): 25-30.
[73] 高玉文, 骆遥, 文武. 2012. 补偿向下延拓方法研究及应用[J]. 地球物理学报, 55(8): 2747,    doi: 10. 6038/j. issn. 0001-5733. 2012. 08.026.
[74] 郭良辉, 孟小红, 石磊,等. 2012. 优化滤波方法及其在中国大陆布格重力异常数据处理中的应用[J]. 地球物理学报, 55(12): 4078-4088,   doi: 10. 6038/j. issn. 0001-5733. 2012. 12.020.
[75] 胡敏章, 李建成. 2010. 顾及地壳密度模型的 Airy-Heiskanen 均衡重力异常的计算[J].   大地测量与地球动力学, 30(5): 48-52.
[76] 黄方, 刘琼颖, 何丽娟. 2012. 晚喜山期以来四川盆地构造-热演化模拟[J]. 地球物理学报, 55(11): 3742-3753,    doi: 10. 6038/j. issn. 0001-5733. 2012. 11.021.
[77] 黄建平, 傅容珊, 许萍,等. 2006. 利用重力和地形观测反演中国及邻区地壳厚度[J]. 地震学报, 28(3): 250-258,    doi: 10. 3321/j. issn: 0253-3782. 2006. 03.004.
[78] 江为为, 刘少华, 郝天珧,等. 2003. 冲绳海槽及其邻域地球物理场与地壳结构特征[J]. 地球物理学进展, 18(2): 287-292,    doi: 10. 3969/j. issn. 1004-2903. 2003. 02.017.
[79] 姜文亮, 张景发. 2012. 首都圈地区精细地壳结构--基于重力场的反演[J]. 地球物理学报, 55(5): 1646-1661,    doi: 10. 6038/j. issn. 0001-5733. 2012. 05.022.
[80] 李润江, 佟淑娟, 孟令顺,等. 1986. 东北地区重力均衡异常特征的初步研究[J].   地震地质, 8(3): 15-19.
[81] 李姗姗, 吴小平, 李建伟. 2004. 基于 Airy-Hesikannen 均衡模型的莫霍面反演[J]. 测绘科学, 29(2): 25-27,    doi: 10. 3771/j. issn. 1009-2307. 2004. 02.008.
[82] 林淼, 朱建军, 田玉淼,等. 2012. 大地水准面异常在湖南地区的地球物理解释[J]. 地球物理学报, 55(2): 472-483,    doi: 10. 3969/j. issn. 0001-5733. 2012. 02.011.
[83] 刘元龙, 郑建昌, 武传珍. 1987. 利用重力资料反演三维密度界面的质面系数法[J]. 地球物理学报, 30(2): 186-196,    doi: 10. 3321/j. issn: 0001-5733. 1987. 02.009.
[84] 马杏垣. 1989. 中国岩石圈动力学地图集[M]. 北京: 中国地图出版社.
[85] 秦静欣, 郝天珧, 徐亚,等. 2009. 关于中国海陆莫霍面深度图编绘的思考[J]. 地球物理学进展, 24(5): 1593-1601,    doi: 10. 3969/j. issn. 1004-2903. 2009. 05.007.
[86] 宋学锋. 2011. 冲绳海槽地区断裂构造解释及莫霍面反演[J].   油气地球物理, 9(2): 48-50.
[87] 童永春. 1985. 江苏地区莫霍面与地壳均衡特征初探[J].   地质与勘探, 21(9): 37-41.
[88] 汪汉胜, 贾路路, Wu P,等. 2012. 末次冰期冰盖消融对东亚历史相对海平面的影响及意义[J]. 地球物理学报, 55(4): 1144-1153,    doi: 10. 6038/j. issn. 0001-5733. 2012. 04.010.
[89] 王勇, 许厚泽. 1996. 中国大陆及其邻区岩石层挠曲强度变化和均衡补偿机制[J].   地球物理学报, 39(S): 105-113.
[90] 徐德琼, 李全兴, 蒋家祯. 1983. 东海莫霍面及其地质意义[J].   海洋通报, 2(5): 34-41.
[91] 许萍, 傅容珊, 黄建平. 2004. 地震层析成像地壳上地幔密度异常与莫霍面深度[C]. // 李国营. 2004年重力学与固体潮学术研讨会暨祝贺许厚泽院士70寿辰研讨会.   武汉: 湖北科学技术出版社.
[92] 杨宗仁. 1986. 地壳重力均衡的定量研究--以湖南地区为例[J]. 地壳形变与地震, 6(4): 285-291.
[93] 殷秀华, 刘铁胜. 1993. 均衡重力异常和地壳表, 浅层地质结构[J].   地震地质, 15(2): 149-156.
[94] 殷秀华, 史志宏, 刘占坡,等. 1982. 华北北部均衡重力异常的初步研究[J].   地震地质, 4(4): 27-34.
[95] 袁学诚. 1996. 中国地球物理图集[M]. 北京: 地质出版社.
[96] 曾融生. 1973. 莫霍界面的重力补偿和地壳结构的基本模式[J].   地球物理学报, 16(1): 1-5.
[97] 曾融生. 1984. 固体地球物理学导论[M].   北京: 科学出版社.
[98] 曾华霖. 2005. 重力场与重力勘探[M]. 北京: 地质出版社.
[99] 曾华霖, 万天丰. 1999. 两幅全国均衡重力异常图的差异[J]. 地球物理学报, 42(1): 127-134,    doi: 10. 3321/j. issn: 0001-5733. 1999. 01.015.
[100] 张赤军. 2003. 地壳对不同波长地形在重力场中的响应[J]. 地球物理学进展, 18(2): 342-347,    doi: 10. 3969/j. issn. 1004-2903. 2003. 02.027.
[101] 张赤军, 方剑, 边少锋,等. 2007. 重力均衡与地形波长--兼论均衡理论不适宜于珠峰重力值的推估[J]. 武汉大学学报·信息科学版, 32(3): 229-233,    doi: 10. 3969/j. issn. 1671-8860. 2007. 03.011.
[102] 赵玉合. 2007. 均衡重力异常计算与综合解释--以黄海地区为例 [硕士论文]. 北京: 中国科学院地质与地球物理研究所.