E-mail Alert Rss
 

物探与化探, 2020, 44(5): 1161-1171 doi: 10.11720/wtyht.2020.0067

方法研究·信息处理·仪器研制

球坐标系密度界面反演方法及在华南大陆的应用

王祥,1,2, 郭良辉1

1.中国地质大学(北京) 地球物理与信息技术学院,北京 100083

2.中国地质调查局 昆明自然资源综合调查中心,云南 昆明 650000

Density interface inversion method in spherical coordinates and its application in the South China mainland

WANG Xiang,1,2, GUO Liang-Hui1

1.School of Geophysics and Information Technology, China University of Geosciences (Beijing), Beijing 100083, China

2.Kunming Natural Resources Comprehensive Survey Center of China Geological Survey, Kunming 650000, China

责任编辑: 王萌

收稿日期: 2020-02-17   修回日期: 2020-04-21   网络出版日期: 2020-10-20

基金资助: 国家自然科学基金面上项目.  41774098

Received: 2020-02-17   Revised: 2020-04-21   Online: 2020-10-20

作者简介 About authors

王祥(1991-),男,中国地质调查局昆明自然资源综合调查中心助理工程师,主要研究方向为重磁数据处理和反演算法。Email: real_shuan_wang@163.com

摘要

密度界面反演方法在油气勘探、推断区域构造、获取地壳结晶基底面、莫霍面起伏形态等方面具有重要意义。现有的密度界面反演方法大多基于笛卡尔坐标系统,当涉及到大区域乃至全球尺度的密度界面反演时,地球曲率的影响将不可忽视,需考虑基于球坐标系Tesseroid模型的密度界面反演方法。然而,受计算精度和效率制约,已有的基于Tesseroid模型的密度界面反演方法并不能很好地适用于地表重力观测数据的反演计算。笔者基于前人研究,给出了一种适用于地表观测数据的球坐标系密度界面反演方法。该方法首先将常规的球坐标系高斯—勒让德重力积分公式进行简化,提高了重力正演计算效率。随后,引入并改进了前人的自适应剖分方案,提高了重力正演计算精度。在此基础上,采用Cordell迭代优化算法,得到了适用于地表观测数据的球坐标系密度界面反演方法。通过模型数据试验,对本文球坐标系密度界面反演方法进行了验证。结果表明,笔者对高斯—勒让德积分公式加以改进和引入改进的自适应剖分方案后,很好地克服了计算精度和效率对地表观测重力计算的掣肘,并且,基于球坐标系的密度界面反演结果优于基于笛卡尔坐标系的密度界面反演结果。在华南大陆实际数据试验中,利用文中方法得到的华南大陆的莫霍面深度与前人得到的莫霍面深度具有较高的吻合度,莫霍面自西向东逐渐抬升,西部抬升剧烈,东部抬升平缓,沿武陵山—贵桂交界一带呈现明显的NNE向莫霍面梯级带,验证了文中方法的科学有效性。

关键词: 重力 ; 球坐标系 ; 密度界面反演 ; Tesseroid ; 华南大陆

Abstract

The density interface inversion method has been playing an important role in the oil and gas exploration, regional structure inference studies as well as crustal crystal basement surface and Moho undulations researches. Most of the density interface inversion methods are generally based on the Cartesian coordinate system. When large regional or even global scale data are dealt with, the influence of earth curvature cannot be ignored, and the density interface inversion method based on Tesseroid model of spherical coordinate system needs to be considered. However, due to the limitations of calculation accuracy and efficiency, the existing density interface inversion method based on Tesseroid cannot be well applicable to the surface gravity observation data. In this paper, on the basis of previous studies, a density interface inversion method of spherical coordinate system suitable for surface observation data is proposed. Firstly, the gravity Gauss-Legendre integral formula in the spherical coordinate system is simplified to improve the forward calculation efficiency. Then, an optimized adaptive subdivision algorithm is introduced to enhance the calculation accuracy. According to the previous forward calculation and by using Cordell iterative optimization algorithm, the authors propose a density interface inversion method for the surface observation data in the spherical coordinate system. The proposed inversion method in this paper can be verified through the synthetic data test. The inversion results show that the proposed method can overcome the limitation of calculating precision and efficiency of the surface observation data. In addition, the inversion results based on spherical coordinate system are better than those based on cartesian coordinate system. Finally, tests on real data from South China mainland verify the feasibility of the presented methods. The results show that Moho depth rises gradually from the west to the east, with the western part uplifting dramatically and the eastern part uplifting gently. Between Wuling Mountain and Guizhou-Guangxi border, there is an obvious NNE-Moho step.

Keywords: gravity ; spherical coordinate system ; density interface inversion method ; Tesseroid ; South China mainland

PDF (4704KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

王祥, 郭良辉. 球坐标系密度界面反演方法及在华南大陆的应用. 物探与化探[J], 2020, 44(5): 1161-1171 doi:10.11720/wtyht.2020.0067

WANG Xiang, GUO Liang-Hui. Density interface inversion method in spherical coordinates and its application in the South China mainland. Geophysical and Geochemical Exploration[J], 2020, 44(5): 1161-1171 doi:10.11720/wtyht.2020.0067

0 引言

地下存在密度差异的岩石或地层接触带上往往形成明显的密度分界面,其形态可以通过密度分界面两边介质的密度差异和形态不规则所引起的重力异常反演得到。密度界面反演方法在推断区域构造、地壳结晶基底面分布情况、获取莫霍面深度及其起伏形态等方面具有重要意义。目前,密度界面的反演方法众多,总体上可分为空间域密度界面反演方法和频率域密度界面反演方法。空间域密度界面反演方法大体包括直接迭代法[1,2,3,4,5]、压缩质面法[6,7,8]和非线性反演方法[9,10,11]。频率域的密度界面反演方法则是通过引入快速傅里叶变换(FFT)进而提高计算效率[12,13,14,15,16]。为了得到更加精确的反演效果,一些学者在常密度界面反演方法的基础上,引入变密度界面模型进行反演[17,18]。尽管已有的密度界面反演方法已经比较成熟,然而,这些方法的绝大多数都是针对基于笛卡尔坐标系统的小范围研究区域。然而,当涉及到区域性乃至全球尺度问题时,地球曲率问题将不容忽视[19,20,21,22,23],此时,基于笛卡尔坐标系统的界面反演方法已不再适用,需要引入基于球坐标系统的Tesseroid模型[24]进行反演计算。

当前,基于Tesseroid模型的重力反演研究主要集中在三维密度成像和莫霍面反演。Liang[25]将Tesseroid模型的重力三维反演方法应用到月球GRAIL重力数据反演中,得到了月球的三维密度结构。Zhang[26]利用高斯—勒让德积分法构建核函数并使用重力梯度全张量对月球GRAIL重力异常进行反演,得到了月球的三维密度分布。同样,Uieda[23]也利用Tesseroid模型对卫星重力数据进行球坐标系快速非线性反演,在反演计算过程中,作者利用参考莫霍面深度、正则化参数和莫霍面起伏重力异常对反演结果加以约束,得到了南美洲大陆的莫霍面深度。前人[23,25-27]所使用的重力反演方法在针对卫星和探测器等高空观测数据时能够取得很好的反演效果,然而,由于邻近效应的存在(当近距离观测时,误差会随着观测距离的减小而急剧增大),这些方法并不能很好地适用于近地表重力观测数据的反演解释。

实现球坐标系地表观测数据的密度界面高精度反演的前提是实现其高精度的正演计算。目前,Tesseroid模型的重力正演计算方法主要有泰勒(Taylor)级数积分法[20,27-30]和高斯—勒让德积分法(Gauss Legendre Quadrature,缩写GLQ)[19,21-23,26],结合前人的研究,两种方法在应用于近地表观测时,其计算精度都会大打折扣,其中,泰勒级数积分法尤为严重,但GLQ方法的计算效率又要略微低于泰勒级数积分法。受限于球坐标系重力正演精度和效率,前人在反演时多数采用的是卫星和探测器等高空数据[23,25,27],而非地表观测数据。由于卫星数据精度较实际地表观测数据低,会给反演结果带来误差。然而,目前鲜有适用于球坐标系地表观测数据的密度界面反演方法。

本文基于Asgharzadeh[19]的球坐标系GLQ重力正演计算方法,给出了一种适用于地表观测数据的球坐标系密度界面反演方法。首先,对重力正演公式进行了改进,将三重GLQ积分计算(3D-GLQ)简化为两重GLQ积分计算(2D-GLQ),提高了计算效率。同时,为了提高计算精度,在正演过程中引入了Tesseroid网格自适应剖分方案[22]。最后,利用Cordell迭代优化算法[2,31]进行迭代反演。通过理论模型和华南大陆的实测数据试验,对本文方法进行了验证。

1 方法与原理

1.1 球坐标系Tesseroid单元体重力正演方法

1.1.1 重力正演公式

Tesseroid单元体(图1)概念最初由Anderson[24]引入,与笛卡尔坐标系统中的直立棱柱体模型类似,它是由6个边界面((Sr1,Sr2)、(Sφ1,Sφ2)、(Sλ1,Sλ2))包裹而成的扇形球壳块体。考虑到实际观测中的经纬度情况,文中采用的是局部北向的球坐标系结构[26]

图1

图1   Tesseroid单元体空间几何示意

Fig.1   The sketch of Tesseroid model


图1Q(r',φ',λ')为场源内任意一点,P(r,φ,λ)为观测点。根据位场理论,Tesseroid单元体在P点产生的引力位V由下式表示[19,20]:

V(r,φ,λ)=λ1λ2φ1φ2r1r2r'2cosφ'ldr'dφ'dλ',

式中:G为万有引力常数(6.67×10-11 m3/kg/s2);ρ为单元体剩余密度(kg/m3);lPQ两点间的欧氏距离。如果ψ为球半径r'r的径向夹角,则重力异常g(r,φ,λ)r可由下式表示[19,20]:

g(r,φ,λ)r=-V(r,φ,λ)r=λ1λ2φ1φ2r1r2r'2cosφ'(r-r'cosψ)l3dr'dφ'dλ',

其中:

l=r2+r'2-2rr'cosψ,cosψ=sinφsinφ'+cosφcosφ'cos(λ'-λ),

将其在r方向上进行积分计算[20],三重体积分就简化为两重面积分,得到如下表达式:

g(r,φ,λ)r= rλ1λ2φ1φ2r1r2r'2cosφ'(r-r'cosψ)l3dr'dφ'dλ'= λ1λ2φ1φ2cosφ'[r'3l-l(r'+3rcosψ)-  r2(3cos2ψ-1)ln(l+r'-rcosψ)]|r'=r1r'=r2dφ'dλ',

当Tesseroid单元体区间范围足够小时,可以利用GLQ积分公式对式(3)的积分计算结果进行近似,进一步可得到:

g(r,φ,λ)r(φ2-φ1)(λ2-λ1)4·  i=1j=1ωiωjcosφ'r×  [r'3l-l(r'+3rcosψ)-   r2(3cos2ψ-1)ln(l+r'-rcosψ)]|r'=r1r'=r2,

其中:

λsi=(λ2+λ1)2+λ'si(λ2-λ1)2,ωi=2nλP-1(λ'si)P'(λ'si)φsj=(φ2+φ1)2+φ'sj(φ2-φ1)2,ωj=2nφP-1(φ'sj)P'(φ'sj)

式中:λsiφsj分别为积分区间内经度、纬度方向的高斯节点;λ'siφ'sj为各坐标方向上(-1,1)区间的高斯节点;P-1(λ'si)为高斯节点λ'si-1阶GLQ多项式的值;P'(λ'si)为高斯节点λ'si阶GLQ多项式的一阶导数值。

式(4)中r方向上的计算精度主要受Tesseroid单元体大小、GLQ多项式阶次和式(2)中l的大小影响。需要特别强调的是,由于引力公式的局限性,当近距离观测时,观测精度对l的大小比较敏感。本文选择GLQ积分法的优势在于该方法可以将Tesseroid表面的积分计算点移至Tesseroid单元体内部,避免了计算时奇异点的出现,该方法相较于泰勒级数积分法精度已经有了明显的提升,但也略微增大了计算量。

1.1.2 自适应剖分方案

为了进一步提高计算精度,本文还引入了自适应剖分方案,前人[22,32]的剖分方案揭示了在满足同一精度要求的情况下,可以根据Tesseroid单元体中的各个区域与观测点P(r,φ,λ)距离不同而采取不一样的网格大小进行剖分,即近距离采用小网格,远距离采用大网格(图2)。

图2

图2   Tesseroid单元体网格剖分示意

a—常规剖分的网格边长大小平均,LφLλ为Tesseroid单元体上表面纬度和经度方向上边长;b—自适应剖分方案的剖分采用的是在观测点近距离处小网格剖分,远距离处采用大网格剖分

Fig.2   Sketch of discretization of Tesseroid

a—routine discretization shown that the tesseroid had a uniform division. Lφ, Lλ are the dimensions of the tesseroid’s top surface;b—adaptive discretization shown that a fine division of the tesseroid close the computation point and a coarser division further away


笔者针对已经改进的GLQ积分公式,对前人的自适应剖分方案也进行了简化改进。对于一个径向、纬向和经向的计算范围分别为[r1,r2]、[φ1,φ2]、[λ1,λ2]的Tesseroid单元体,其几何中心点坐标为(rc,φc,λc),其与观测点P(r,φ,λ)之间的距离为d。对一个Tesseroid单元体是否进行进一步的剖分,可以由式(5)来进行判定。

dLW,

其中:

d=r2+rc2-2rrccosψccosψc=sinφsinφc+cosφcosφccos(λ-λc)Lφ=πr2(φ2-φ1)180,Lλ=πr2cosφ1(λ2-λ1)180,L=max(Lφ,Lλ),

式中:Lφ,Lλ根据Tesseroid单元体的弧长与周长的比例得来;W为根据不同的计算精度要求而给定的距离判定值。

基于此,适合本文2D-GLQ积分公式的网格加密方案可以分为以下几步:

1) 计算d,Li的值,判定关系式(5)是否成立,如果成立,直接计算重力值。

2) 如果式(5)不成立,直接将当前单元体剖分成网格大小为Δφ'×Δλ'N×N个单元体(式(6)),则新生成的所有小单元体都符合式(5)的判定条件。将新生成的所有小的Tesseroid单元体在观测点处的重力值进行累加就得到了加密剖分前的Tesseroid单元体重力异常值。

3) 进行下一个Tesseroid单元体的计算,重复步骤1),2)。

Δφ'=(φ2-φ1)N,Δλ'=(λ2-λ1)N,

其中:

N=CWLd,

式中:C为向上取整符号。对前人自适应剖分方案改进后,公式更加简化,对计算精度的控制也变得更加细腻,而且还大大减少了W值的循环判定次数,在一定程度上节省了计算时间。笔者将2D-GLQ积分公式与改进的自适应剖分方案相结合的重力正演方法称为GLQ-Plus。

1.2 球坐标系密度界面正演方法

相较于传统的笛卡尔直角坐标系统,球坐标系统中所选择的参考面为球面,如图3所示,参考球面可以位于Tesseroid顶界面,也可以位于底界面。这样,区域内由m×n个Tesseroid单元体所组成的计算区域,其在观测点P处的重力异常gP为各个单元网格块体在P点处所引起重力异常的叠加。

图3

图3   正演计算模型剖面示意

Fig.3   The forward model of tesseroid


gP=m=1Mn=1N[g(r,φ,λ)r]mn

1.3 球坐标系密度界面反演方法

基于Cordell迭代算法流程,本文球坐标系密度界面迭代反演流程如下图4所示:

图4

图4   Cordell迭代算法流程

Fig.4   The flow chart for Cordell iterative approach


1)利用Bott[33]公式对Tesseroid单元体的初始厚度进行估算:

H1,mn=Kgmnobs,

式中:K= 12π; gmnobsPmn点处的重力异常观测值。

2)由式(7)中可以得到重力异常计算初值 g1,mncalc。重力异常观测值与计算值的均方误差RMS表示如式(9)所示,如果RMS符合要求,则输出区域密度界面深度。

RMS=m=1Mn=1N(gmnobs-gmncalc)2M×N,

3)如果RMS不符合要求,则计算下一次重力正演计算前的深度修改值(式(10)),然后继续执行第2步运算,如此循环迭代,直到RMS满足要求为止。

Hs+1,mn=Hs,mngmnobsg1,mncalc

2 理论模型数据试验

2.1 Tesseroid单元体模型试验与评价

为了验证文中方法的精确性,根据前人的研究[34,35,36],我们建立球壳体模型(图5),P为位于极轴上的观测点,O为球心,R2为球壳外表面半径,R1为球壳内表面半径,D为观测点P与球心之间的距离。

图5

图5   球壳体模型示意

a—球壳内视;b—球壳侧视

Fig.5   The images of spherical shell models

a—inside view of a shell;b—side view of a shell


利用球壳体模型求取重力精确解析值,将本文方法和常规方法球壳体的重力计算值与球壳体的理论重力值进行对比,以此达到验证计算精度和效率的目的。球壳体的理论重力值gtrue的计算关系式:

gtrue=43πR23-R13D2,

此次模型试验中,ρ取0.3 g/cm3,R2取6.371×106 m,R1取6.361×106 m,球壳厚度为104 m,D-R2为观测点P与球壳体外球面之间的距离(观测高度H),笔者通过不同的观测高度来验证各方法对于观测高度的敏感度(图6)。从图中可看出,对于不同的网格剖分大小,同阶次的Taylor和2D-GLQ随着观测高度的提升,两种方法的计算结果都越来越接近异常理论值(True Value),网格越小,精度也越高。当观测高度为零,也就是地表观测时,2D-GLQ方法的计算精度始终大大高于Taylor方法,并且,随着观测高度的增大,在一定的观测高度区间内,2D-GLQ方法的计算结果会略微高于异常理论值,这些都是2D-GLQ计算特点所决定的[19]。然而,当网格大小仅为3'×3'时,仍然难以达到地表观测的高精度计算要求,所以,需要进一步提高计算精度。

图6

图6   不同网格大小下Taylor和2D-GLQ在不同观测高度处的重力异常计算值与理论异常值的对比

a—1°×1°;b—30'×30';c—15'×15';d—3'×3'

Fig.6   The gravity anomaly calculated by the Taylor on the different cell size and 2D-GLQ methods based on the different height of observation surface


泰勒级数展开法的公式推导过程繁琐,目前它已有的重力正演计算公式最高阶次只有二阶,更高阶次(四阶、六阶等)的公式推导所带来的工作量一直让人望而却步,而GLQ方法可以很轻松地实现更高阶次的计算,配合自适应剖分方案,使得在实现同等精度情况下,计算耗时大大减少,实现了地表观测重力正演计算的更高精度、更高效率。我们用Er来表示本文方法重力计算值gGLQ-Plus与理论重力值gtrue的误差,则:

Er=gtrue-gGLQ-Plus,

由于不同的距离判定值W所对应的误差级数相差较大,所以,我们采取误差的10的对数形式来表示(图7a)。从分析结果可以看出,随着2D-GLQ公式阶次的提高,计算精度会有明显的提升。当W为64时,本文2D-GLQ公式经度和纬度计算阶次都为8或12时,重力异常误差分别为0.042 9、0.019 5 mGal,与理论异常的百分误差(=Er/gtrue)分别仅为0.017 2%、0.007 8%。但是,随着W的不断增大,精度提升速率也变缓。此外,在取得同样计算误差的前提下,我们计算了采取本文的二阶2D-GLQ公式时,常规剖分计算时间T1与自适应剖分方案计算时间T2的比值T1/T2(图7b),以此来验证本文引入自适应剖分方案后对计算效率的提升。计算机的中央处理器为Inter(R) Core(TM) i7-6700HQ,随机存取存储器(RAM)为8 GB,验证程序由MATLAB R2014a编译。从对比结果分析来看,当精度要求越高时,本文方法对时间的节省效果越明显,不仅如此,加上对GLQ公式的改进,计算效率在自适应剖分方案基础上有了更进一步的提升。

图7

图7   本文方法效率提升试验

a—,不同取值时,Er的十进对数值;b—不同的Er取值,时间耗费比例值T1/T2

Fig.7   The improvement of efficiency of computation by using the method of this paper

a—the Briggs logarithm of Er based on the different value of and ;b—the changing trend of T1/T2 when Er get reduced


此模型试验表明,相对于笛卡尔直角坐标系统,本文利用Tesseroid单元体模型,将球坐标系引入重力正演计算当中,从而避免了因忽略地球曲率对计算结果产生的影响[20]。并且,由于邻近效应的存在,现有的重力正演方法在计算精度和效率上,难以达到地表观测的要求,而本文通过对前人球坐标系重力计算公式和自适应剖分方案的改进,很好地克服了计算精度和效率对地表观测重力计算的掣肘。

2.2 密度界面模型试验

采用起伏密度界面模型(图8a)对本文重力反演方法的有效性进行验证。模型的起伏范围为31~39 km,网格大小为0.2°×0.2°,剩余密度值为-0.3 g/cm3,计算时将地表作为参考面。为了避免难以克服的边界效应,我们将模型区域四周剔除掉0.6°的宽度。

图8

图8   理论模型正演

a—密度界面模型深度分布图;b—GQL-Plus计算的重力异常

Fig.8   The depth of academic interface and its forward result

a—the density interface model graph; b— the gravity anomaly graph of density interface model forward result by the GQL-Plus method of this paper


图8b为通过笔者GQL-Plus算法得到的重力异常。我们直接把利用本文方法正演得到的重力异常值作为实际观测值进行接下来的反演计算,并与其他方法反演结果进行对比(图9)。由图9信息不难发现,基于球坐标系的反演方法得到的密度界面(图9b、c)起伏形态比基于笛卡尔坐标系反演得到的密度界面(图9a)形态更接近实际模型深度起伏形态,说明基于球坐标系的反演方法相较笛卡尔坐标系的反演方法更为稳定。各反演方法得到的界面深度与实际模型深度的差值(图9d~f)说明本文基于球坐标系的密度界面反演方法的反演结果与实际模型深度最为接近,本文反演方法的反演精度要远远高于其他两种方法的反演精度(图9g)。笔者反演方法在迭代5次过后,各参数的均方差(RMS)已经趋于稳定,当2D-GLQ公式为8阶,W取64时,笔者方法的深度RMS能低至8 m,异常RMS能低至0.003 542 mGal。

图9

图9   各方法的理论模型反演对比试验

a—笛卡尔坐标系反演结果;b—基于二阶泰勒级数法正演的反演结果;c—基于本文2D-GLQ正演方法的反演结果;d、e、f—为图(a)、(b)、(c)反演结果与理论模型深度的平面差值;g—图(a)、(b)、(c)对角剖面与理论深度值对比;h—深度RMS随迭代次数的走势图;i—布格异常RMS随迭代次数的走势图

Fig.9   Comparisons of the inversion results from the data of Figure 8(b)

a—Cartesian coordinates;b—two order Taylor;c—2D-GLQ inversion results; d,e,f—respectively are errors of (a), (b), and (c); g—shown the diagonal section comparisons of (a), (b) and (c);h—depth RMS;i—bouguer RMS


3 华南实际数据试验

华南大陆位于中国秦岭、淮河以南,横断山脉以东。本文研究中选取的华南大陆经度范围为104°E~122°E,纬度范围为21°N~32°N。前人[37,38]认为,华南大陆在早前寒武纪多块体复杂构造运动基础

上,又遭受全球大陆聚散和南北大陆离散拼合的联合构造应力作用,是欧亚板块、太平洋板块和印度洋板块的拼接挤压带,构造变形复杂,岩浆活动强烈。而到了中晚中生代(燕山期)时期,统一的华南陆块又经历了陆内构造改造和西太平洋板块西向俯冲,使得西北部陆块呈现自东向西推进的活化变形,而东南部陆块区域呈现出面状分布的岩浆侵入与构造变形。

重力数据来源于世界地质图委员会(CGMW)公布的世界重力网格数据库 WGM 2012,将其网格化得到华南大陆布格重力异常数据(图10a),经度范围为104°E~122°E,纬度范围为21°N~32°N,数据网度为0.2°×0.2°。WGM 2012数据库是由国际重力测量局依据高分辨率的地球重力模型 EGM2008 和高程模型 ETOPO1采用球谐方法计算得到的,计算中考虑了实际的地球模型和大多数地表物质质量分布。WGM2012布格重力异常数据不仅包含了莫霍面起伏引起的重力异常,还包含了地壳中浅部地质体和地幔干扰所产生的重力异常,所以,为了得到主要由莫霍面起伏所引起的区域重力异常,我们进行参数测试,最终选择带宽为800 km的低通滤波处理,去除浅部高频异常的干扰,得到网格间距为0.2°×0.2°的区域重力异常数据(图10b),异常范围大致为-250~10 mGal。笔者在反演过程中对剩余密度值和参考深度进行了测试,最终选择剩余密度值为-0.28 g/cm3、参考深度值为30 km。本次莫霍面反演采用6阶2D-GLQ公式,W值为32的重力正演计算方法,利用Cordell迭代算法进行反演计算,迭代4次以后,结果趋于稳定,得到了华南大陆莫霍面深度(图10c)。反演结果的理论重力异常与区域实际重力异常RMS为0.582 mGal。

图10

图10   华南大陆实测数据试验

a—华南大陆重力异常滤波前异常;b—华南大陆重力异常低通滤波后区域异常;c—本文得到的华南大陆莫霍面;d—反演迭代收敛曲线;e—华南大陆莫霍面[39];f—c与e的莫霍面深度差值

Fig.10   Comparisons of the inversion results from measured data in South China

a— The Bouguer anomaly in South China;b—The Moho gravity anomaly in South China by using the low-pass filter;c—The Moho depth in South China by the inversion method of this paper;d—The trend chart of depth RMS;f—The Moho depth in South China[39];g—The difference of (c) and (e)


本文方法得到的华南大陆莫霍面深度在26~46 km之间,平均深度为35 km,这与前人[39]得到的该区域莫霍面深度范围(图10e)大致吻合,但局部地方也存在一些差异(图10f)。将反演得到的莫霍面深度(图10c)与重力异常作对比发现,莫霍面深度的分布规律与异常的分布规律基本一致,即异常高值区域,对应莫霍面凹陷,异常低值区,对应莫霍面隆起。莫霍面等值线呈EW向展布,总体自西向东逐渐递减。在东部临海地区莫霍面相对较浅,而四川、云南地区莫霍面深度较深,沿武陵山—贵桂交界一带呈现明显的NNE向莫霍面梯级带,莫霍面从西向东抬升5 km。

华南大陆莫霍面西深东浅的起伏特征,与该区域经历了西太平洋俯冲[40,41]与喜马拉雅造山及青藏高原隆升作用[38]有关。四川和云南地区莫霍面自西向东抬升最为剧烈,推测与该区域扬子板块受青藏高原隆升挤压,使得地壳厚度由西向东剧烈减薄有关。沿武陵山—贵桂交界一带以东地区,莫霍面平缓抬升,反应该地区地壳大面积减薄,前人[40,41]认为中国东南大陆边缘存在古太平洋俯冲带,该俯冲带可能位于东海西缘—东沙东缘深大断裂带,越靠近该俯冲带,地壳呈现逐渐减薄的趋势,并且厚度变化越趋于平缓,这与本文研究区域反演所得到的东部沿海地区莫霍面深度特征相吻合。

4 结论

笔者利用了Tesseroid模型构建方法,将球坐标系常规重力正演计算中的高斯—勒让德积分公式加以改进,并将引入的Tesseroid网格自适应剖分方案公式和流程进行了简化,提高了球坐标系重力正演计算效率和计算精度。其次,在实现了球坐标系密度界面重力高精度正演的基础上,我们应用Cordell迭代反演方法,得到了适用于地表观测数据的球坐标系密度界面高精度反演方法。在Tesseroid单元体模型试验中,对本文GLQ-Plus方法的地表观测计算精度和计算效率进行了验证,结果表明,受计算精度和计算效率的制约,常规的球坐标系重力正演方法并不能很好地适用于地表观测的重力正演计算,而本文对高斯—勒让德积分公式加以改进和引入改进的自适应剖分方案后,情况得到了明显的改善。在密度界面模型试验中,将笛卡尔坐标系的密度界面反演方法和球坐标系的密度界面反演方法进行了对比,该模型试验不仅表明了基于球坐标系的密度界面反演结果比笛卡尔坐标系的密度界面反演结果更加稳定可靠,还验证了本文给出的球坐标系密度界面反演方法可以实现地表重力观测数据的高精度反演(深度RMS<8 m)。在华南大陆实测数据试验中,我们得到了华南大陆的莫霍面深度,并与前人得到的莫霍面深度进行了对比,验证了本文方法的科学有效性。

参考文献

孙德梅, 闵志.

三维密度界面反演的一个近似方法

[J]. 物探与化探, 1984,8(2):89-98.

URL     [本文引用: 1]

在区域地球物理研究工作中,应用重力资料反演密度界面,尤其是计算莫氏面的深度,已被广大学者所重视。目前某些计算方法,如线性公式法、压缩质面法等已被应用。我们认为,在利用小比例尺区域重力资料反演密度界面时,现有的一些方法各有其一定的局限性。

Sun D M, Min Z.

A technique for approximate inversion of three-dimensional density interface

[J]. Geophysical and Geochemical Exploration, 1984,8(2):89-98.

URL     [本文引用: 1]

Based upon an approximate formula of three-dimensional gravity anomalies and by using the successive approaching method, an inversion technique for the computation of density data in regional gravity survey has been obtained. The technique was applied to process gravity data from_the middlelower reaches of the Yantze River, and a map of Moho isobaths of the area was constructed.Acomprison of calculated results with seismic sounding data shows a fairly good agreement. Acomputer program for the technique is being developed for the interpretation of regional gravity data.

Cordell L, Henderson R G.

Iterative three-dimensional solution of gravity anomaly data using a digital computer

[J]. Geophysics, 1968,33(4):596-601.

DOI:10.1190/1.1439955      URL     [本文引用: 2]

Silva J B C.

Interactive gravity inversion

[J]. Geophysics, 2006: 1-9.

[本文引用: 1]

Zhou X.

Gravity inversion of 2D bedrock topography for heterogeneous sedimentary basins based on line integral and maximum difference reduction method

[J]. Geophysical Prospecting, 2013,61(1):220-234.

DOI:10.1111/j.1365-2478.2011.01046.x      URL     [本文引用: 1]

Based on the line integral (LI) and maximum difference reduction (MDR) methods, an automated iterative forward modelling scheme (LI-MDR algorithm) is developed for the inversion of 2D bedrock topography from a gravity anomaly profile for heterogeneous sedimentary basins. The unknown basin topography can be smooth as for intracratonic basins or discontinuous as for rift and strike-slip basins. In case studies using synthetic data, the new algorithm can invert the sedimentary basins bedrock depth within a mean accuracy better than 5% when the gravity anomaly data have an accuracy of better than 0.5 mGal. The main characteristics of the inversion algorithm include: (1) the density contrast of sedimentary basins can be constant or vary horizontally and/or vertically in a very broad but a priori known manner; (2) three inputs are required: the measured gravity anomaly, accuracy level and the density contrast function, (3) the simplification that each gravity station has only one bedrock depth leads to an approach to perform rapid inversions using the forward modelling calculated by LI. The inversion process stops when the residual anomalies (the observed minus the calculated) falls within an error envelope whose amplitude is the input accuracy level. The inversion algorithm offers in many cases the possibility of performing an agile 2D gravity inversion on basins with heterogeneous sediments. Both smooth and discontinuous bedrock topography with steep spatial gradients can be well recovered. Limitations include: (1) for each station position, there is only one corresponding point vertically down at the basement; and (2) the largest error in inverting bedrock topography occurs at the deepest points.

Silva J B C, Santos D F, Gomes K P.

Fast gravity inversion of basement relief

[J]. Geophysics, 2014,79(5):79-91.

[本文引用: 1]

刘元龙, 王谦身.

用压缩质面法反演重力资料以估算地壳构造

[J]. 地球物理学报, 1977,20(1):59-69.

[本文引用: 1]

Liu Y L, Wang Q S.

Inversion of gravity data by use of a method of “compressed mass plane” to estimate crustal structure

[J]. Chinese J. Geophys. , 1977,20(1):59-69.

[本文引用: 1]

江为为, 周立宏, 肖敦清, .

东北地区重磁场与地壳结构特征

[J]. 地球物理学进展, 2006,21(3):730-738.

URL     [本文引用: 1]

分析了东北地区的重、磁场特征,同时对研究区的布格重力异常和航磁异常据进行了小波分析计算.根据分析与计算可知,东北地区重力场以北东走向为主,表现出该地区重力场的主要趋势.根据磁场的分布特征,可将研究区分为六个区域:呼和浩特以北磁场相对平静区;赤峰一带正负磁场交互变化区;海拉尔以南磁场缓变区;齐齐哈尔—依春磁场剧烈变化区;长春—沈阳负磁场平缓区;牡丹江—丹东磁场区.利用重力资料,应用调和级数法对研究区的莫霍界面进行了反演计算,得到了该地区莫霍界面深度分布.根据磁力资料采用遗传算法反演计算了研究区居里界面的深度分布.同时对研究区的地壳结构特征进行了探讨.

Jiang W W, Zhou L H, Xiao D Q, et al.

The characteristics of crust structure and the gravity and magnetic fields in northeast region of China

[J]. Progress in Geophysics, 2006,21(3):730-738.

[本文引用: 1]

胡立天, 郝天珧.

带控制点的三维密度界面反演方法

[J]. 地球物理学进展, 2014,29(6):2498-2503.

DOI:10.6038/pg20140603      URL     [本文引用: 1]

利用重力资料反演密度界面对于研究盆地、区域构造和深部构造有着重要意义.但是密度基准面深度和界面密度差的不准确性降低了反演结果的可靠性.本文使用已知深度的控制点,在逐步迭代中计算出最合适的密度基准面深度和界面密度差,使反演结果和控制点拟合最好.理论模型和南海地区的莫霍反演结果表明本文方法可以计算出最合适的密度基准面深度和界面密度差,取得理想的反演结果.

Hu L T, Hao T Y.

The inversion of three-dimensional density interface with control points

[J]. Progress in Geophysics, 2014,29(6):2498-2503.

DOI:10.6038/pg20140603      URL     [本文引用: 1]

The inversion of density interface using gravity data has great values in basin research、regional tectonics and deep structure. However, the uncertainty of density datum plane and density contrast makes the inversion result unreliable. An inversion method of three-dimensional density interface with several known control points is presented in this paper. We calculate the datum plane and density contrast in each iteration to make the inversion result fit the control points best. The tests of synthetic model and south China sea gravity data prove that the method obtains suitable datum plane and density contrast and gets a good inversion result.

Marquardt D W.

An algorithm for least-squares estimation of nonlinear parameters

[J]. Journal of the Society for Industrial and Applied Mathematics, 1963,11(2):431-441.

DOI:10.1137/0111030      URL     [本文引用: 1]

Chakravarthi V, Sastry S R.

GUI based inversion code for automatic quantification of strike limited listric fault sources and regional gravity background from observed bouguer gravity anomalies

[J]. Journal Geological Society of India, 2014,83:625-638.

[本文引用: 1]

Mojica O F, Bassrei A.

Regularization parameter selection in the 3D gravity inversion of the basement relief using GCV: A parallel approach

[J]. Computers and Geosciences, 2015,82:205-213.

DOI:10.1016/j.cageo.2015.06.013      URL     [本文引用: 1]

关小平.

利用Parker公式反演界面的一种有效方法

[J]. 物探化探计算技术, 1991,13(3):236-242.

[本文引用: 1]

Guan X P.

An effective and simple approach of subsurface inversion using Parker’s equation

[J]. Computing Techniques for Geophysical and Geochemical Exploration, 1991,13(3):236-242.

[本文引用: 1]

Parker R L.

The rapid calculation of potential anomalies

[J]. Geophysical Journal Royal Astronomical Society, 1972,37(3):447-455.

DOI:10.1111/j.1365-246X.1974.tb04096.x      URL     [本文引用: 1]

Oldenburg D W.

The inversion and interpretation of gravity anomalies

[J]. Geophysics, 1974,39:526-536.

DOI:10.1190/1.1440444      URL     [本文引用: 1]

Granser H.

Three-dimensional interpretation of gravity data from sedimentary basins using an exponential density-depth function

[J]. Geophysical Prospecting, 1987,35(9):1030-1041.

[本文引用: 1]

Guspi F.

Noniterative nonlinear gravity inversion

[J]. Geophysics, 1993,58(7):935-940.

DOI:10.1190/1.1443484      URL     [本文引用: 1]

Santos D F, Silva J B C, Martins C M, et al.

Efficient gravity inversion of discontinuous basement relief

[J]. Geophysics, 2015,80(4):95-106.

[本文引用: 1]

Shi L, Li Y H, Zhang E H.

A new approach for density contrast interface inversion based on the parabolic density function in the frequency domain

[J]. Journal of Applied Geophysics, 2015,116:1-9.

DOI:10.1016/j.jappgeo.2015.02.022      URL     [本文引用: 1]

Asgharzadeh M F, Von Frese R R B, Kim H R, et al.

Spherical prism gravity effects by Gauss-Legendre quadrature integration

[J]. Geophys. J. Int, 2007,169:1-11.

DOI:10.1111/gji.2007.169.issue-1      URL     [本文引用: 6]

Heck B, Seitz K.

A comparison of the Tesseroid prism and pointmass approaches for mass reductions in gravity field modelling

[J]. J. Geod., 2007,81(2):121-136.

DOI:10.1007/s00190-006-0094-0      URL     [本文引用: 6]

Roussel C, Verdun J, Cali J, et al.

Complete gravity field of an ellipsoidal prism by Gauss-Legendre quadrature

[J]. Geophys. J. Int., 2015,203:2220-2236.

DOI:10.1093/gji/ggv438      URL     [本文引用: 2]

Uieda L, Barbosa V C F, Braitenberg C.

Tesseroids: Forward-modeling gravitational fields in spherical coordinates

[J]. Geophysics, 2016: 41-48.

[本文引用: 3]

Uieda L, Barbosa V C F.

Fast nonlinear gravity inversion in spherical coordinates with application to the South American Moho

[J]. Geophys. J. Int., 2017,208:162-176.

DOI:10.1093/gji/ggw390      URL     [本文引用: 5]

Anderson E G.

The effect of topography on solutions of stokes' problem

[M]. School of Surveying, University of New South Wales, Kensington, NSW, Australia, 1976.

[本文引用: 2]

Liang Q, Chen C, Li Y.

3-Dinversion of gravity data in spherical coordinates with application to the GRAIL data

[J]. J. Geophys. Res. Planets, 2014,119:1359-1373.

DOI:10.1002/2014JE004626      URL     [本文引用: 3]

Zhang Y, Wu Y L, Yan J G, et al.

3D inversion of full gravity gradient tensor data in spherical coordinate system using local north oriented frame

[J]. Earth, Planets and Space, 2018,70(58):1-23.

DOI:10.1186/s40623-017-0766-4      URL     [本文引用: 3]

梁青.

月球重力异常特征与三维密度成像研究:

[D]. 武汉:中国地质大学(武汉), 2010.

[本文引用: 3]

Liang Q.

Gravity anomaly features and 3D density imaging of the moon:

[D]. Wuhan:China University of Geosciences (Wuhan), 2010.

[本文引用: 3]

Wild-Pfeiffer F.

A comparison of different mass elements for use in gravity gradiometry

[J]. J Geod., 2008,82:637-653.

DOI:10.1007/s00190-008-0219-8      URL    

Grombein T, Seitz K, Heck B.

Optimized formulas for the gravitational field of a Tesseroid

[J]. J. Geodesy, 2013,87:645-660.

DOI:10.1007/s00190-013-0636-1      URL    

Shen W B, Deng X L.

Evaluation of the fourth-order tesseroid formula and new combination approach to precisely determine gravitational potential

[J]. Stud. Geophys. Geod, 2016,60:1. (DOI: 10.1007/s11200-016-0402-y).

URL     [本文引用: 1]

Guo L H, Shi L, Meng X H, et al.

Apparent magnetization mapping in the presence of strong remanent magnetization: The space-domain inversion approach

[J]. Geophysics, 2015,81(2):J11-J24.

DOI:10.1190/geo2015-0082.1      URL     [本文引用: 1]

Li Z, Hao T, Xu Y, et al.

An efficient and adaptive approach for modeling gravity effects in spherical coordinates

[J]. Journal of Applied Geophysics, 2011,73(3):221-231.

DOI:10.1016/j.jappgeo.2011.01.004      URL     [本文引用: 1]

The need to obtain more reliable Earth structures has been the impetus for conducting joint inversions of disparate geophysical datasets. For seismic arrival time tomography, joint inversion of arrival time and gravity data has become an important way to investigate velocity structure of the crust and upper mantle. However, the absence of an efficient approach for modeling gravity effects in spherical coordinates limits the joint tomographic analysis to only local scales. In order to extend the joint tomographic inversion into spherical coordinates, and enable it to be feasible for regional studies, we develop an efficient and adaptive approach for modeling gravity effects in spherical coordinates based on the longitudinal/latitudinal grid spacing. The complete gravity effects of spherical prisms, including gravitational potential, gravity vector and tensor gradients, are calculated by numerical integration of the Gauss-Legendre quadrature (GLQ). To ensure the efficiency of the gravity modeling, spherical prisms are recursively subdivided into smaller units according to their distances to the observation point. This approach is compatible with the parameterization of regional arrival time tomography for large areas, in which both the near- and far-field effects of the Earth's curvature cannot be ignored. Therefore, this approach can be implemented into the joint tomographic inversion of arrival time and gravity data conveniently. As practical applications, the complete gravity effects of a single anomalous density body have been calculated, and the gravity anomalies of two tomographic models in the Taiwan region have also been obtained using empirical relationships between P-wave velocity and density. (C) 2011 Elsevier B.V.

Bott M H P.

The use of rapid digital computing methods for direct gravity interpretation of sedimentary basin

[J]. Geophysical Journal Royal Astronomical Society, 1960,3:63-67.

DOI:10.1111/gji.1960.3.issue-1      URL     [本文引用: 1]

Heiskanen W A, Moritz H.

Physical geodesy

[M]. San Francisco:Freeman, 1967.

[本文引用: 1]

Vaníček P, Novák P, Martinec Z.

Geoid, topography, and the Bouguer plate or shell

[J]. J. Geodesy, 2001,75:210-215.

DOI:10.1007/s001900100165      URL     [本文引用: 1]

 Topography plays an important role in solving many geodetic and geophysical problems. In the evaluation of a topographical effect, a planar model, a spherical model or an even more sophisticated model can be used. In most applications, the planar model is considered appropriate: recall the evaluation of gravity reductions of the free-air, Poincaré–Prey or Bouguer kind. For some applications, such as the evaluation of topographical effects in gravimetric geoid computations, it is preferable or even necessary to use at least the spherical model of topography. In modelling the topographical effect, the bulk of the effect comes from the Bouguer plate, in the case of the planar model, or from the Bouguer shell, in the case of the spherical model. The difference between the effects of the Bouguer plate and the Bouguer shell is studied, while the effect of the rest of topography, the terrain, is discussed elsewhere. It is argued that the classical Bouguer plate gravity reduction should be considered as a mathematical construction with unclear physical meaning. It is shown that if the reduction is understood to be reducing observed gravity onto the geoid through the Bouguer plate/shell then both models give practically identical answers, as associated with Poincaré's and Prey's work. It is shown why only the spherical model should be used in the evaluation of topographical effects in the Stokes–Helmert solution of Stokes' boundary-value problem. The reason for this is that the Bouguer plate model does not allow for a physically acceptable condensation scheme for the topography.]]>

Karcol R.

Gravitational attraction and potential of spherical shell with radially dependent density

[J]. Stud. Geophys. Geod., 2011,55:21-34.

DOI:10.1007/s11200-011-0002-9      URL     [本文引用: 1]

Solutions to the direct problem in gravimetric interpretation are well-known for wide class of source bodies with constant density contrast. On the other hand, sources with non-uniform density can lead to relatively complicated formalisms. This is probably why analytical solutions for this type of sources are rather rare although utilization of these bodies can sometimes be very effective in gravity modeling. I demonstrate an analytical solution to that problem for a spherical shell with radial polynomial density distribution, and illustrate this result when applied to a special case of 5th degree polynomial. As a practical example, attraction of the normal atmosphere is calculated.

舒良树.

华南构造演化的基本特征

[J]. 地质通报, 2012,31(7):1035-1053.

[本文引用: 1]

Shu L S.

An analysis of principal features of tectonic evolution in South China Block

[J]. Geological Bulletin of China, 2012,31(7):1035-1053.

[本文引用: 1]

张国伟, 郭安林, 王岳军, .

中国华南大陆构造与问题

[J]. 中国科学:地球科学, 2013,43(10):1553-1582.

[本文引用: 2]

Zhang G W, Guo A L, Wang Y J, et al.

Tectonics of South China continent and its implications

[J]. Science China: Earth Sciences, 2013,43(10):1553-1582.

URL     [本文引用: 2]

郭良辉, 孟小红, 石磊, .

优化滤波方法及其在中国大陆布格重力异常数据处理中的应用

[J]. 地球物理学报, 2012,55(12):4078-4088.

DOI:10.6038/j.issn.0001-5733.2012.12.020      URL     [本文引用: 3]

在优选延拓法的理论基础上,研究提出基于格林等效层概念和维纳滤波器的优化滤波法,用于对重力异常数据进行去噪和分离.与传统向上延拓法和优选延拓法相比,优化滤波法分离异常与延拓高度无关,不需要已知延拓高度,具有一定的优势.理论重力模型数据的去噪和异常分离试验表明优化滤波法有效,异常分离效果优于传统向上延拓法和带通滤波法.利用优化滤波法对中国大陆重力异常数据去噪和异常分离,得到有效的布格重力异常和区域重力异常.以中国大陆深地震探测推断的莫霍面深度信息为约束,对区域重力异常数据进行密度界面约束反演,得到中国大陆莫霍面深度分布.本文方法为中国大陆深部探测和区域构造研究提供一定的技术支撑.

Guo L H, Meng X H, Shi L, et al.

Preferential filtering method and its application to Bouguer gravity anomaly of Chinese continent

[J]. Chinese J. Geophys. , 2012,55(12):4078-4088.

[本文引用: 3]

Hilde T W C, Uyeda S, Kroenke L.

Evolution of the western Pacific and its margin

[J]. Tectonophysics, 1977,38:145-165.

DOI:10.1016/0040-1951(77)90205-0      URL     [本文引用: 2]

Smith A D.

A plate model for Jurassic to recent intraplate volcanism in the Pacific Ocean basin

[G]// Foulger G R, Jurdy D M. Plates, plumes, and planetary processes. Geol Soc Am (Spec Papers), 2007,430:471-495.

[本文引用: 2]

/

京ICP备05055290号-3
版权所有 © 2021《物探与化探》编辑部
通讯地址:北京市学院路29号航遥中心 邮编:100083
电话:010-62060192;62060193 E-mail:whtbjb@sina.com