基于重力的莫霍面反演方法比较

柴永杰

〔深圳市勘察测绘院(集团)有限公司,广东 深圳 518028〕

0 引 言

地球外部重力场属于地球的空间物理场范畴,其表征了地球内部物质的密度、体积等物理和几何量的整体空间分布情况。因此,重力方法理论上来说可以作为一种手段来研究地球的内部物质分布。随着重力卫星GRACE[1]和GOCE[2]的成功发射,人们获取了史无前例的高分辨率、高精度和(几乎)全球覆盖的静态和时变重力场信息,从而为采用重力方法研究地球内部物质分布注入了新的活力和手段。Moho面为地壳与地幔之间的分界面。该面一个显著的特点是密度跳跃(也称不连续),因而采用重力方法研究Moho面的基本原理即是根据该界面密度跳跃特性所产生的重力异常建立起重力异常与Moho面参数之间的函数关系。

确定Moho面深度具有很广泛的地学应用。例如可以更好地理解和认识地震触发的机制、地热的流向、热平衡及其演化和板块构造等地学过程。在定义与扩展大陆架方面,Moho面的深度也是重要的依据[3]。在大地测量方面,Moho深度也有很多重要的应用,例如Moho面深度可以用来计算均衡重力异常。均衡重力异常相对于自由空气重力异常和布格重力异常更加平滑,因此为了减小重力异常内插和外推的误差通常采用均衡重力异常用于内插和外推[4]。另一方面,为了提高GOCE重力梯度观测值向下延拓反演平均海平面上重力异常值的稳定性和精度,可以利用地形-均衡信息对GOCE沿轨观测的重力梯度值进行归算,而该归算需要用到Moho面的深度信息。地形-均衡信息还可以用来构建高阶重力场模型例如EGM2008[5]的高阶部分的球谐系数主要来源于地形-均衡的位系数。在地球物理领域,Moho面结合地形信息可以用来反演岩石圈的有效弹性厚度[6];皮娇龙等[7]研究了青藏高原区域Moho面起伏对通道流模型的动力学响应。在地质方面Moho面深度与成矿之间也存在关联,如陈安国等的研究所述[8]。

1 研究进展

现阶段主要有地震和重力两种方法来确定Moho面的深度。当然也可采用重力和地震联合的方法。地震方法的基本原理是基于地震波在Moho面处的反射和折射方程来建立数学模型,而重力方法则是基于Moho面的密度不连续性所引起的重力异常大小来建立观测方程。地震和重力方法是基于不同的原理和假设,因此这两种方法得出的结果必然存在差异[9]。

随着新一代重力卫星GRACE[1]和GOCE[2]的成功发射,人类跨入了卫星重力遥感时代,从而可以快速、高效地获取高分辨率、高精度和(几乎)全球覆盖的均匀重力场信息。在此卫星重力遥感时代背景下,采用重力方法研究地壳结构能够有效的克服地震数据的不均匀性和分辨率低等缺点,尤其在地震数据稀疏甚至完全缺失的情况下。

采用重力模型反演Moho面深度,国内外学者做出了大量的工作。本文认为这些方法可以概括为4种:① 以地壳均衡假设为基础导出的算法,例如Vening Meinesz-Moritz(VMM)方法[9-10];② 以FFT算法为基础的Parker-Oldenburg 区域算法[11-12];③ 以球谐分析和球谐综合为基础的Parker-Oldenburg扩展算法[13],且该算法也能利用FFT技术;④ Moho引力信号提取的直接方法其包括一阶泰勒线展开性化法、凝聚法。需要指出的是,以上各种算法之间也是相互联系和一脉相承的,并非完全独立。本文对以上各种算法进行了回顾和分析,并指出其优缺点。

1.1 地壳均衡假设

地壳均衡假设是指地壳处于静态平衡状态,即地壳的密度和厚度在不同区域内达到平衡。Moho面是地球的地壳和地幔之间的界面,是地壳和地幔密度和速度差异的表现。地壳均衡假设和Moho面之间存在密切的关系。根据地壳均衡假设,地壳的密度和厚度是在不同的地质时间尺度上达到平衡的,而Moho面是地壳和地幔之间的物理界面,随着地质时间的演化,其深度和形态也在不断变化。因此,研究地壳均衡假设可以帮助我们更好地理解Moho面的形成和演化过程。此外,地壳均衡假设也可以用于解释地球内部的一些现象,例如地震波传播的速度和路径,地壳形变等。通过对地壳均衡假设和Moho面的研究,我们可以更深入地了解地球内部的结构和演化过程。

在19世纪中期,Pratt在喜马拉雅山脉地区通过测量发现实测的垂线偏差小于可见的地形质量所计算的理论结果,因此他推断地球内部一定存在某种质量亏损来解释上述差异。另外,假定地形仅仅是简单的叠加在匀质的地壳之上,则布格重力归算将会扣除地形起伏所引起的大部分不规则重力场信号,因此布格重力异常理论上将会非常的平滑且其大小应该在零值附近摆动,然而事实却并不如此且布格重力异常具有系统性特征,即在山区表现为大的负值。对于以上理论与实际不符的情况,唯一的解释就是在山区地形以下有质量的亏损,换句话说,地形质量存在某种补偿。本文简单回顾一下3种经典的地壳均衡模型。

1.1.1 普拉特-海福特模型

Pratt首先提出了该模型的核心思想,后来Hayford推导了详细的数学计算公式并且把该模型广泛的应用于大地测量学中,因此通常把该模型称之为Pratt-Hayford模型。Pratt-Hayford模型(图1)的基本原理是位于补偿面以下的物质密度是均匀的,位于补偿面以上每个截面的密度随着地形高度的变化而变化。即位于补偿面以上每个截面的质量相等,因此高出海水面的地形越高,则该截面的密度越小;低于海水面的地形越深,则该截面的密度越大。

图1 Pratt-Hayford均衡模型

1.1.2 艾黎-海斯卡涅模型

该模型是Airy首先提出,为了便于大地测量学的应用,Heiskanen给出了精确的数学表达式,并对其进行了广泛的应用。因此通常把该模型称之为Airy-Heiskanen模型(图2)。该模型的基本思想是地形漂浮于密度较大的岩浆(地幔)上面,因此地形越高,则在重力的引导下其下沉得越深。在山区,其截面密度为常数ρc=2.67 g/cm3,下面密度较大的地幔密度为ρm=3.27 g/cm3,且定义地壳和地幔密度差为Δρc/m=ρm-ρc。根据Airy模型,在山区下面会形成山根,海洋下面会形成反山根。

图2 Airy-Heiskanen 均衡模型

地震研究表明,Airy-Heiskanen均衡模型比Pratt-Hayford均衡模型更加的接近实际地球状态,但是Pratt-Hayford均衡模型也存在一定的合理性,所以有些学者提出了联合均衡模型,即在陆地上采用Airy-Heiskanen均衡模型,在海洋上采用Pratt-Hayford均衡模型。尽管如此,实际情况也比联合均衡模型要更复杂,例如地球的曲率影响、地壳的弹性性质、冰川均衡调整(GIA)等因素都会对均衡状态有所影响。

1.1.3 维宁曼尼斯区域模型

通过考察Pratt-Hayford和Airy-Heiskanen模型,Vening Meinesz发现上述两种模型是局部补偿机制,也就说补偿仅仅沿着界面的垂直方向,显然这不符合实际情况,因此Vening Meinesz在Airy模型的基础上进行了改进,提出了区域补偿模型。局部补偿和区域补偿的差异如图3所示。根据Vening Meinesz假设,地形被视为一种负载施加于刚性的弹性地壳上面,也就说其考虑了地壳本身的弹性性质,因此部分地形负载不仅仅会通过垂直方向也会通过水平方向进行补偿。通过以上分析,可知Airy-Heiskanen模型是Vening Meinesz模型的一种特殊情况。

图3 区域补偿与局部补偿之间的差异

现实中地球均衡补偿机制非常复杂,涉及的因素很多,例如地质构造运动、冰川均衡调整等地球动力学过程。为了尽可能逼真的对实际地球均衡状态进行描述,在Moho面反演建模的过程中充分吸收Pratt-Hayford、Airy-Heiskanen、Vening Meinesz和Moritz的思想,即地形的均衡补偿是由于密度补偿、山根和反山根在全球尺度上的综合补偿所致,该思想可以抽象成是一种“广义均衡模型”。

1.2 Vening Meinesz-Moritz算法研究现状

Pratt模型是以地壳密度的亏损或者剩余来补偿地形的质量。而Airy模型是以地壳厚度的亏损或者剩余去补偿地形的质量。Pratt和Airy均衡补偿模型都是一种局部(Local)补偿模式。因此,从物理角度分析,这两种补偿模型都忽略了地壳张力等其他因素对地壳均衡的影响,为此Vening Meinesz(1931)引入岩石圈的挠曲模型对Airy均衡模型进行了进一步的拓展,提出了区域(Regional)补偿模型。后来,Moritz[10]进一步将Vening Meinesz区域均衡模型扩展成全球(Global)均衡模型。Sjöberg[9]对Moritz[10]全球均衡模型进行了数学上的重新表述,并将求解Moho面深度问题归纳为求解第一类非线性Fredholm积分方程,并进一步推导了顾及非线性二次项的直接反演公式。Moritz[10]和Sjöberg[9]方法都是以全球均衡思想为基础,但不同的是前者给出的是迭代解,而后者给出的是基于球近似导出的直接解。

Vening Meinesz-Moritz(VMM)算法虽然成功的反演了不同区域的Moho面深度,但是本文认为该算法存在3个缺点:① 当包括到非线性项的二次项时该反演方程包括一个二维空间积分项,即任意一点的积分理论上都需要进行全球积分导致算法的计算效率较低;② 当计算点和积分点之间的球面距离为零的时候上述空间积分会出现奇异值,因此还需要进一步近似处理才能消除奇异值的影响(例如,计算点附近采用平面近似);③ 当顾及非线性项的三次项及更高次项的时候导出的算法会更加的复杂,因此限制了该算法的反演效率和精度。

1.3 Parker-Oldenburg算法研究现状

Parker[12]提出了一种重力场正演的快速算法,该算法的基本思想是在平面直角坐标系下(平面近似)将任意界面所引起的重力异常表达为傅里叶级数。接着,Oldenburg[11]对Parker的正演算法进行了简单的数学变换并提出了采用重力异常迭代求解界面深度的反演算法。同时,为了提高迭代收敛的速度和稳定性,提出了一个滤波器,后来人们通常称该算法为Parker-Oldenburg(P-O)算法。本文认为该算法的一个显著优点是采用了FFT技术,从而大大提高了计算效率,同时也提高了对分辨率较高的格网数据的处理能力。同时,该算法的一个明显缺点是其基于平面近似,因此只能用来进行小范围的局部研究,即范围不能过大,否则地球的曲率对反演结果影响不容忽视。最初,Parker-Oldenburg算法只能用于二维重力数据的反演。后来,Gomez-Ortiz and Agarwal[14]对Parker-Oldenburg算法进行了扩展,将其由二维扩展到三维,因此大大提高了其应用范围。从均衡的角度看,本文认为Parker-Oldenburg算法也可以看成在Vening Meinesz区域均衡模型基础上导出的算法,只不过算法适用于平面近似的区域,但是其补偿范围比Airy均衡模型要广泛。

1.4 Parker-Oldenburg扩展算法研究现状

正如前面所论述,虽然Parker-Oldenburg计算效率高,但是由于其是基于平面近似导出的算法,因此该算法只能用于平面近似区域。为了能够将FFT技术也能应用于全球Moho面的快速反演,Chen and Tenzer[13]借鉴Parker-Oldenburg的思想,在球近似的情况下提出了一种新的算法,暂且称之为“Parker-Oldenburg扩展算法”,从而有效提高了反演全球Moho面深度的计算效率。和Vening Meinesz-Moritz方法比较,“Parker-Oldenburg扩展算法”表现出了明显优势,主要体现在4点:① 避免了任意一点的Moho面求解都需要费时的全球积分计算;② 不存在求解方程的奇异值问题;③ 该算法的二次以及更高次项都可以利用FFT技术进行计算且一次、二次及更高次项的公式具有一致性,即不存在阶次越高算法越复杂的情况;④ 该算法可以顾及Moho面的横向变化(二维变化)从而有效的提高反演结果的精度。但是本文认为“Parker-Oldenburg扩展算法”仍然存在以下缺点:① 平均Moho面深度和密度差选择不够精确会导致迭代不收敛;② 该算法不能够根据陆地和海洋的特征分别选择先验参数,因此导致反演的结果在海洋和陆地上均具有一定的系统性误差。

从均衡的角度去考察,本文认为“Parker-Oldenburg扩展算法”也可以理解为以“广义均衡模型”为基础导出的算法。因为该算法在数学建模的时候,既考虑了Moho面深度的变化,同时也考虑了Moho面密度差的变化,换句话说,高出平均海水面的地形,既不是单一的Pratt模型提出的密度亏损来补偿,也不是单一的Airy模型提出的地壳厚度来补偿,而是密度和地壳联合补偿,同时该算法也包含了Moritz的全球均衡思想。即可以认为该算法综合了Pratt、Airy以及Moritz 3种均衡模型的思想,认为“广义均衡模型”是“Parker-Oldenburg扩展算法”的物理基础。众所周知,地壳的均衡状态是一个复杂的物理过程,其主要取决于岩石圈的负载、岩石圈有效的弹性厚度、硬度及其流变学特性、软流圈的黏度等因素。同时,冰川均衡调整、冰川融化、板块运动、地幔对流等地球动力学过程也对地壳均衡产生一定的影响。因此,本文认为基于地壳均衡模型的重力模型存在一定的物理缺陷,其主要体现在目前还难以模型化的地球动力学效应所引起的重力异常信号如何扣除,例如采用重力方法所反演的Moho面在俯冲带区域具有比较大的缺陷。

1.5 直接法研究现状

从信号恢复的角度看,不基于任何均衡假设,而是从重力观测值中利用地震、地球物理等辅助数据恢复Moho面所引起的重力异常信号,然后采用扣除各种重力异常之后的精化重力数据进行Moho面反演,本项目称之为“直接法”。直接法大体上可分为两步进行,首先将计算点至Moho面之间以及Moho面到地心之间的所有重力异常信号从观测值中进行扣除;然后把扣除各种重力异常信号之后的精化重力异常当作输入数据反演Moho面深度。

按照数学处理的不同,本文将其大致分为一阶泰勒线展开性化方法和凝聚法等。从数学形式上看,泰勒展开线性化方法就是对第一类非线性Fredholm积分方程利用泰勒展开定理进行线性化,但是其展开的位置选择平均Moho面深度处,从而可大大简化算法的复杂度和计算量。凝聚法的基本思想是将Moho面相对于平均Moho面的起伏压缩到平均Moho面深度所对应的球面上,然后导出相应的算法。需要指出的是直接法的一个核心问题是如何干净、有效的扣除非Moho面所引起的重力异常信号。

随着计算机技术的不断发展和重力数据获取的进一步完善,重力反演在Moho面深度研究中的应用越来越广泛,而且具有很大的潜力。首先,重力数据的精度和覆盖范围将会得到进一步提高。其次,多源数据融合将会成为一种趋势。重力反演方法可以与其他地球物理学方法结合使用,例如地震勘探、热流场分析、地表形貌等数据,这些数据的融合可以增加反演结果的可靠性和精度。第三,反演算法的改进和发展也将会推动重力反演Moho深度的研究,未来随着新的反演算法的发展,反演结果的准确性和可靠性将会进一步提高。

2 结 语

通过对重力反演Moho算法的研究现状剖析,未来可以对重力反演Moho面算法从如下几个方面进行改进。

(1)重力数据的优点是分辨率高、精度高和覆盖范围广;地震方法的优点是反演结果精度和可靠性高,缺点是数据采集成本高,导致地震数据分辨率和精度不均匀、某些区域分布稀疏甚至缺失导致这些区域的反演结果具有较大的不确定性。因此,第一个切入点是推导联合重力和地震优点的新算法。

(2)当前计算非Moho面所引起重力异常的时候采用的地壳密度通常为常密度(一维)或者横向变化的密度(二维),但是我们知道地壳真实密度是随着深度和位置三维变化,因此采用一维或者二维的密度值去替代三维的密度值必然会产生误差。因此,第二个切入点是联合多源数据精化地壳模型,然后推导顾及地形密度三维变化的重力场正演算法。

(3)基于地壳均衡模型的重力算法存在一定的物理缺陷,主要是地球动力学效应(例如地壳形变)难以模型化;重力算法根据研究范围的大小采用一定的近似:局部的采用平面近似,区域和全球的采用球近似。由于真实地球的形状更接近椭球,平面和球近似必然会产生近似误差。因此,第三个切入点是在数学上采用更加接近真实地球形状的椭球作为数学推导的基准面;在物理上联合GNSS、卫星重力和水文数据量化地壳形变与Moho面变化之间的关系。

重力反演Moho深度在地球科学中的应用将会更加广泛。随着对地球内部结构和演化的研究不断深入,重力反演Moho深度不仅可以用于地震、地质勘探等领域,还可以用于地球物理学、地球化学等领域的研究,有望在未来取得更多的研究成果。