2019, Vol. 62
2. 长安大学重磁方法技术研究所, 长安大学地质工程与测绘学院, 西安 710054
2. Institute of Gravity and Magnetic Technology, College of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China
当地质界面上下存在密度差时就构成了密度界面，如海底地形、沉积盆地基底和莫霍面.就密度界面的几何形态而言，一类重要的密度界面在剖面上呈Ⅴ字型，如海沟和海底峡谷地形、裂陷盆地中的半地堑以及俯冲带之下的莫霍面.此类密度界面形态的刻画对于区域构造研究和油气勘探等较为关键，对于物理海洋学、海洋生物学和海洋化学研究亦具有重要的意义(Smith and Sandwell, 1997).
密度界面上下的密度差较大时，可引起明显的重力异常，因此直接根据重力异常能反演得到密度界面起伏形态，通常可采用迭代的方法逐步消除剩余异常值来实现.这一反演计算方法可以在空间域进行(Bott, 1960; Cordell and Henderson, 1968; 林振民和阳明, 1985; 刘元龙等, 1987; Leäo et al., 1996; Chakravarthi and Sundararajan, 2004, 2007a, b; Zhou, 2013; Silva et al., 2014; Chakravarthi et al., 2016; Pallero et al., 2015, 2017)或者在频率域进行(Oldenburg, 1974; Pilkington and Crossley, 1986; 冯锐等, 1986; 王万银和潘作枢, 1993; Shin et al., 2006; 冯娟等, 2014; 张恩会等, 2015).重力场具有连续且光滑的特征，而以上方法仅根据重力异常进行计算，因此得到密度界面的形态与重力异常形态相关，亦表现为光滑的形态.这对于一些密度界面如克拉通盆地基底是合理的，但并不适用于本文提及的Ⅴ字型密度界面.对于Ⅴ字型密度界面，尽管以上方法能得到一个与实际界面接近的光滑模型，然而其并非一定是真实形态的模糊反映(Martins et al., 2011).
对于非光滑形态密度界面的反演，可利用正则化方法(Tikhonov and Arsenin, 1977)在反演时施加约束信息以提高反演的准确性，并能保证结果的稳定性.这一非线性反演方法通过建立一个稳定的泛函并利用最优化方法求解其极小化问题以达到反演的目的.此稳定泛函通常由两部分组成：(1)数据误差函数，用来保证反演结果吻合实测重力异常；(2)正则化函数，用以保证反演结果的稳定性并符合已知地质或地球物理特征.在密度界面反演中，正则化函数一般可由两项构成：已知深度约束函数和模型约束函数(Barbosa et al., 1997).前者利用待反演界面的已知深度作为约束，以保证反演结果接近真实深度，其几乎不影响界面的形态；后者利用不同范数建立约束函数，控制相邻剖分模型的变化规律以保证反演结果符合地质特征，其决定了密度界面的形态.其中，在保证非光滑密度界面形态方面，Barbosa等(1999)通过在Barbosa等(1997)建立的L2-范数形式的模型约束函数中加入加权矩阵而反演了非光滑形态的沉积盆地基底.Silva等(2010)将一阶熵测度最小化和零阶熵测度最大化结合起来，利用熵正则化反演了二维非光滑沉积盆地基底.Martins等(2011)将全变差函数作为正则化项，实现了非光滑沉积盆地基底三维反演.全变差函数的实质为L1-范数，与加权方法(Barbosa et al., 1999)和熵正则化方法(Silva et al., 2010)相比，全变差正则化方法更为简单，并且不需要基底最大深度等先验信息(Lima et al., 2011).在此基础上，冯旭亮等(2014)利用全变差函数进行了双界面模式的裂陷盆地基底二维重力反演，Xing等(2016)进行了多层界面几何结构的2.5维重力反演.
尽管L1-范数形式的模型约束函数可实现非光滑形态的密度界面反演，但其反演结果为矩形状的非光滑形态(冯旭亮，2016)，这与Ⅴ字型界面的特征并不相符.事实上，密度界面反演中模型约束函数常用的L2-范数和L1-范数可用Lp-范数统一表示.不同范数形式的模型约束函数反演得到的密度界面形态不同，说明当p取不同值时，反演结果会呈现不同的特征.基于此，本文利用Lp-范数建立模型约束函数，利用模型试验验证了p=5时能刻画Ⅴ字型密度界面的起伏形态，之后将该方法用于马里亚纳海沟南段挑战者深渊及其邻区莫霍面深度反演之中.
1 方法原理 1.1 目标函数密度界面起伏正则化反演的目标函数可写为以下形式(Barbosa et al., 1999; Martins et al., 2010; 冯旭亮等, 2014)：
|
(1) |
在密度界面反演时，通常将界面之上的密度层剖分为M个垂直并置的棱柱体，用棱柱体的底界面深度近似密度界面的深度，因此(1)式中m为M×1型的密度界面深度向量，其为待反演的参数. λh和λm为正则化参数；φg(m)是数据误差函数，可用实测重力异常与正演拟合重力异常的L2-范数的平方表示(Silva et al., 2010)：
|
(2) |
其中gobs和g(m)分别为实测重力异常和正演拟合重力异常.
(1) 式中，h(m)为已知深度约束函数，其形式为(Martins et al., 2011)：
|
(3) |
其中H为一个B×1型向量，表示了密度界面的已知深度，B为已知深度点的个数；W为B×M型矩阵，其每行只有一个非零元素，相当于单位矩阵(Martins et al., 2011).矩阵W的非零元素可取为1，第b行的非零元素与向量H的第b个元素有关，若该行中第i个元素为非零元素，则其所在位置处反演的密度界面深度与第b个已知深度点的水平距离最小，这样通过(3)式最小化可保证在某个已知深度点附近的反演结果接近真实深度.
(1) 式所示的目标函数中，τ(m)为模型约束函数，其决定了反演得到的密度界面的形态，用范数形式表示.为刻画Ⅴ字型密度界面的形态，本文采用Lp-范数建立τ(m)的表达式.Lp-范数的定义如下：
|
(4) |
其中p∈[1, ∞).为简单起见，通常采用Lp-范数的p次方形式：
|
(5) |
尽管Lp-范数有较好的通用性，但应用于正则化反演时，其数值计算很难实现，因为对于向量的任一元素mi，其导数为p|mi|p－1，当p=1时，在mi=0点不可微(Sun and Li, 2014).为解决这一问题，通常采用Lp-范数的近似形式.
Huber(1964)给出了表示向量长度的通用公式：
|
(6) |
其中ρ(m)为关于m的非负函数.
Ekblom(1973)给出了ρ(m)的一种表达式，可称为Ekblom范数：
|
(7) |
其中ε是阈值，当其值相对于m较小时，Ekblom范数与Lp-范数的特征近似.
Ekblom范数对任一元素mi的偏导数为
|
(8) |
可见Ekblom范数是二次可微的.在利用一些最优化方法(如牛顿法)求解目标函数极小化问题时，需要计算Hessian矩阵，要求目标函数是二次可微的，此时Ekblom范数满足目标函数的要求.显然Ekblom范数通用性较好.
利用Ekblom范数的定义，可将模型约束函数τ(m)写为：
|
(9) |
其中，mi和mj表示沿x方向和y方向相邻的两个棱柱体底面深度参数对(Martins et al., 2011; 冯旭亮等, 2014)，L为参数对的个数.
1.2 反演流程本文利用Polak-Ribière-Polyak(PRP)非线性共轭梯度法(Polak and Ribière, 1969)求解目标函数(1)的极小化问题，在Newman和Alumbaugh(2000)算法流程的基础上进行修改，由以下6个步骤组成：
第一步：输入实测重力异常数据gobs、重力观测点位数据(x, y, z)、密度界面最浅深度数据mcons、密度界面已知深度数据H及界面上下密度差Δρ.其中mcons用来约束反演结果，保证反演的密度界面深度大于该值，亦可作为上界面而实现双界面模式反演(冯旭亮等, 2014).设置迭代次数k=1.
第二步：利用Bott(1960)的方法得到密度界面反演的初值m(0).
第三步：计算r(k)=-∇φ(m(k))，u(k)=C(k)－1r(k)，其中C(k)=‖r(k)‖.
第四步：利用一维线搜索方法计算模型更新步长α(k)，使得φ(m(k)+α(k)u(k))最小.
第五步：计算m(k+1)=m(k)+α(k)u(k)，r(k+1)=- ∇φ(m(k+1))；如果|r(k+1)|小于设定的精度则停止并输出m(k+1)作为反演结果，否则进行下一步.
第六步：计算
从非线性共轭梯度法的算法流程来看，其关键在于目标函数梯度的计算.经推导，目标函数(1)的梯度为：
|
(10) |
其中A为雅克比矩阵，其为正演拟合重力异常g(m)对反演的密度界面深度m的偏导数，矩阵A的任一元素aij=∂gi/∂mj (Martins et al., 2011). J(m)为模型约束函数的τ(m)的梯度，其表达式为：
|
(11) |
其中R为一阶差分拉普拉斯算子，q为一个L×1型向量，其第l个元素的表达式为：
|
(12) |
为验证反演方法的准确性并为了清晰地显示反演效果，本文首先建立了两个二维Ⅴ字型密度界面进行试验，之后建立一个三维模型进一步试验方法的有效性.
2.1 简单二维模型简单形态二维模型如图 1中各图下半部分灰色实线所示.该密度界面由两个Ⅴ字型凹陷组成，左边凹陷为对称形式，右边凹陷非对称.给定界面上下的密度差为-0.35×103 kg·m-3，将界面之上的密度层剖分为61个相邻的棱柱体，柱体的水平尺寸为1 km，重力计算点位于每个柱体顶面中心点上，共61个计算点，则该界面引起的重力异常如图 1中各图上半部分灰色圆点所示.利用Lp-范数建立模型约束函数用于反演时，p是一个关键参数，其值影响了反演结果的形态.为确定适合Ⅴ字型密度界面反演的p值，本文令p为1~12分别进行反演.反演时界面上下的密度差为-0.35×103 kg·m-3，模型剖分参数与正演一致，给定阈值ε=10-4，正则化参数λh=0、λm=0.1(即没有采用已知深度约束)，反演得到的界面深度见图 1中各图下半部分中黑色实线.
|
图 1 p取不同值时简单二维模型反演结果 (a)-(l)分别为p取1至12时的反演结果. Fig. 1 Inversion results of 2D simple density interface with different values of p (a) to (l) show results when p is assigned 1 to 12, respectively. |
由图 1可以看出，p取1时反演的密度界面底部较为平坦，与理论模型差别较大；p取2~12时，反演结果均与理论模型较为接近，但其细节不同.p=2时，反演结果呈现光滑形态，尽管其表现出了理论模型的宏观特征，但其仅为理论模型的“模糊”反映，并没有反映出界面的Ⅴ字型特征.随着p取值的增大，反演结果的光滑性降低，与理论模型越来越接近.p=5时，反演结果为Ⅴ字型密度界面，其与理论模型形态一致，仅在模型最底部稍浅一些.随着p取值的进一步增大，反演结果又呈现光滑形态，Ⅴ字型特征逐渐变得不明显.需要提及的是，当p≥13时，反演不收敛，其原因有待进一步探讨.可见，对于Ⅴ字型密度界面的反演，p取5能取得较好的结果.
从图 1亦可以看出，利用Lp-范数建立的模型约束函数对于对称的Ⅴ字型密度界面的反演结果要稍优于非对称的Ⅴ字型密度界面，误差主要集中在非对称Ⅴ字型界面的最深处.推测其原因是重力异常本身对于较深的地质体反映较模糊，且非对称Ⅴ字型界面引起的重力异常并未呈现出明显的非对称性.
实际在反演密度界面时，界面上下的密度差不一定能准确地取为真实密度差，此时p取5能否刻画Ⅴ字型界面的特征？为此将密度差取为-0.385×103 kg·m-3(即比真实密度差大10%)对该简单二维模型进行反演，其他反演参数与图 1一致，结果如图 2所示.当p取4~7时，反演结果均为Ⅴ字型，但浅于理论模型；p≥8时，反演结果为光滑界面.此外，对密度差比真实值小10%(即密度差为-0.315×103 kg·m-3)的情形也进行了试算，结果表明p在1~10之间取值时，反演结果均为深于理论模型的光滑界面，且当p≥11时，反演不收敛.鉴于篇幅原因，这里只展示p取4~7的结果，如图 3所示.
|
图 2 密度差为-0.385×103 kg·m-3，p取不同值时二维简单模型反演结果 (a)-(l)分别为p取1至12时的反演结果.图中各曲线含义与图 1相同. Fig. 2 Inversion results of 2D simple density interface with different values of p when the density contrast is -0.385×103 kg·m-3 (a) to (l) show the results when p is assigned 1 to 12, respectively. Graphic symbols are the same as Fig. 1. |
|
图 3 密度差为-0.315×103 kg·m-3，p取不同值时二维简单模型反演结果 (a)-(d)分别为p取4至7时的反演结果.图中各曲线含义与图 1相同. Fig. 3 Inversion results of 2D simple density interface with different values of p when the density contrast is -0.315×103 kg·m-3 (a) to (d) show the results when p is assigned 4 to 7, respectively. Graphic symbols are the same as Fig. 1. |
需要说明的是，以上试验中给定阈值ε及正则化参数λh和λm进行反演，事实上从反演的方法原理来看，λh和λm的值决定了反演效果.其中λh是已知深度约束(可认为是绝对约束)的权重，其值不会影响反演结果的几何形态，因此反演时可取为1或更大的值.而λm是模型约束函数(可认为是相对约束)的权重，该值决定了反演界面的形态.为明确Ⅴ字型密度界面反演时λm的取值规律，本文亦对图 1所示的模型进行反演试验，计算时分别令λm为0.0001、0.01、0.1、0.5、1和2，其他参数与图 1一致，结果如图 4所示.
|
图 4 p=5，λm取值不同时简单二维模型反演结果 (a)-(f)分别为λm取0.0001，0.01，0.1，0.5，1和2时的反演结果.图中各曲线含义与图 1相同. Fig. 4 Inversion results of 2D simple density interface with different values of λm when p is 5 (a) to (f) show the results when λm is assigned 0.0001, 0.01, 0.1, 0.5, 1 and 2, respectively. Graphic symbols are the same as Fig. 1. |
当λm=0.0001时，模型约束函数τ(m)几乎不起作用，因此反演结果为光滑界面；随着λm增大，反演界面逐渐呈Ⅴ字型；当λm较大时(如λm=1或2)，尽管反演结果仍为Ⅴ字型，但界面最大深度处出现一定误差，且随着λm取值的增大而增大，推测其是由模型约束函数对反演结果的影响大于重力异常的影响而造成的.因此建议在反演时λm的值取0.1~0.5即可.
利用正则化方法反演Ⅴ字型密度界面时，根据(7)式的定义，阈值ε应尽量取较小的值以保证Ekblom范数与Lp-范数的特征近似，如上文反演时阈值ε取为10-4，但其值亦可能影响反演效果.为进一步说明这一问题，我们给ε取值为0.0001、0.01、0.1和1，p = 5，其他参数与图 1一致，对图 1所示的模型进行反演，结果如图 5所示.当ε取0.0001、0.01和0.1时，反演结果均与理论模型较为吻合，而当ε为1时，反演结果为光滑界面，此时ε已不满足(7)式的条件.可见ε在满足定义时，其值对反演结果影响不大.因此在使用时建议ε取值小于0.01.
|
图 5 p=5，ε取值不同时简单二维模型反演结果 (a)-(d)分别为ε取0.0001，0.01，0.1和1时的反演结果.图中各曲线含义与图 1相同. Fig. 5 Inversion results of 2D simple density interface with different values of ε when p is 5 (a) to (d) show the results when ε is assigned 0.0001, 0.01, 0.1 and 1, respectively. Graphic symbols are the same as Fig. 1. |
综上所述，在利用基于Lp-范数的正则化方法反演Ⅴ字型密度界面时，建议p取5，正则化参数λh为1，λm取0.1~0.5，阈值ε取值小于0.01.
2.2 复杂二维模型复杂形态二维密度界面模型如图 6下半部分灰色实线所示.该密度界面整体Ⅴ字型特征，界面最大深度为9.593 km，位于剖面长度x=144 km处，最大深度两侧的界面出现一系列的局部构造，使得该界面形态复杂化.令界面上下的密度差为-0.40×103 kg·m-3，将界面之上的密度层剖分为251个相邻的水平尺寸为1 km的棱柱体，计算其引起的重力异常.重力计算点位于每个柱体顶面中心点上，共251个计算点.为了消除或减弱边部效应，在正演及随后的反演计算时均将模型两个端点的深度值各向外延50个点，并在正演结果中加入均值为0、标准差为0.1×10-5 m·s-2的高斯白噪声以符合实际数据特征，得到该模型的理论重力异常如图 6上半部分灰色圆点所示.
|
图 6 采用3个已知深度点约束时复杂二维密度界面Lp-范数约束反演结果 Fig. 6 Inversion result of 2D complex shape density interface using Lp-norm constraint inversion with three known depth points |
采用3个已知深度点(图 6中黑色十字)作为约束，给定p=5，阈值ε=10-4，正则化参数λh=1、λm=0.1，密度差为-0.40×103 kg·m-3，模型剖面参数与正演参数一致进行反演，结果如图 6下半部分黑色虚线所示.反演结果与理论模型较为吻合，呈现了该复杂形态Ⅴ字型密度界面的起伏形态，仅在剖面长度x=150 km至156 km处，二者有一定的差别.理论模型中该段界面呈“N”字型，而反演结果接近水平界面，并没有反映出该段界面的局部构造.推测原因为该构造相对高差仅为0.5 km，而其埋深超过7 km，具有尺寸小、埋深大的特征，因此根据重力异常无法反演出其起伏特征.
在进行二维简单模型反演试验时，并没有采用已知深度点进行约束，反演结果与理论模型较为吻合.而对于复杂二维模型进行反演时，采用了3个已知深度点，亦得到了较准确的结果.综合以上结果来看，当界面上下的密度差取值准确时，已知深度点似乎没有起到约束的作用.为分析已知深度点在反演中的作用，我们对图 6所示的复杂二维模型进行了两方面的反演试验：一是自左向右将图 6中第1和第3个已知深度点的深度增大10%、将第2个已知深度点的深度减小10%，即已知深度点与真实深度存在误差时进行反演；二是不采用已知深度约束进行反演，结果分别如图 7和图 8所示.
|
图 7 已知深度约束误差10%时复杂二维密度界面Lp-范数约束反演结果 图中各曲线含义与图 6相同. Fig. 7 Inversion result of 2D complex shape density interface using Lp-norm constraint inversion with three known depth points whose errors are 10% Graphic symbols are the same as Fig. 6. |
|
图 8 不采用已知深度约束时复杂二维密度界面Lp-范数约束反演结果 图中各曲线含义与图 6相同. Fig. 8 Inversion result of 2D complex shape density interface using Lp-norm constraint inversion without known depth points Graphic symbols are the same as Fig. 6. |
当已知深度约束出现误差时，尽管反演结果呈现Ⅴ字型特征，但在已知深度点处反演结果和已知深度点吻合，而与理论模型差别较大.可见已知深度约束函数在反演中起到了约束的作用.而不采用已知深度约束的反演结果(图 8)与图 6所示的结果几乎一致，似乎反映了当密度差准确时已知深度约束没有起到作用.为进一步说明已知深度约束的作用，我们将已知深度点增加为4个，其中第4个已知深度点位于剖面长度x=150 km至156 km处局部“N”字型构造的最浅处.利用4个已知深度点作为约束，其他参数与图 6中一致进行反演，结果见图 9.与3个已知深度点约束结果(图 6)相比，采用4个已知深度点的反演结果将整个局部“N”字型构造准确地呈现出来.显然，当密度差为真实值时，已知深度约束函数将地球物理数据与已知地质信息结合起来，会有效地提高复杂密度界面反演的准确性.
|
图 9 采用4个已知深度点约束时复杂二维密度界面Lp-范数约束反演结果 图中各曲线含义与图 6相同. Fig. 9 Inversion result of 2D complex shape density interface using Lp-norm constraint inversion with four known depth points Graphic symbols are the same as Fig. 6. |
目前关于非光滑形态的密度界面反演方面，主要采用L1-范数形式的模型约束函数(相当于p=1的情形).为了对比Lp-范数的反演结果，本文也利用基于L1-范数的二维反演方法(冯旭亮等，2014)进行反演，结果如图 10所示.利用L1-范数作为约束的反演结果整体呈阶梯状的非光滑特征，其并没有呈现理论模型Ⅴ字型的特征，特别是界面深度最大的位置，反演结果与理论模型差别较大.另外在一些细节上，L1-范数反演结果也与理论模型差别较大，如剖面长度x=125 km附近，理论模型有一微小规模的凹陷，而该处反演结果几乎为水平界面.相比而言，Lp-范数的反演结果很好地刻画了这一局部特征.
|
图 10 复杂二维密度界面L1-范数约束反演结果 图中各曲线含义与图 6相同. Fig. 10 Inversion result of 2D complex shape density interface using L1-norm constraint inversion Graphic symbols are the same as Fig. 6. |
密度界面反演方法中，Bott(1960)提出的基于布格板重力正演公式建立直接迭代反演方法(以下简称直接迭代法)也是目前研究和使用较多的方法，本文亦利用该方法对图 6所示的模型进行反演，结果如图 11所示.可以看出，直接迭代法反演结果为光滑界面，其只能反映密度界面的近似形态.另外，该方法对于重力数据误差较为敏感，使得反演界面呈锯齿状跳跃(林振民和阳明, 1985).尽管这一问题可通过对重力数据进行光滑(林振民和阳明, 1985; Silva et al., 2014)来解决，但光滑参数不好选择，其可能会造成反演结果欠稳定或过度稳定(Silva et al., 2014).
|
图 11 复杂二维密度界面Bott(1960)方法反演结果 图中各曲线含义与图 6相同. Fig. 11 Inversion result of 2D complex shape density interface using Bott′s (1960) method Graphic symbols are the same as Fig. 6. |
三维密度界面模型如图 12b所示.界面在x= 8 km处深度最大，最大深度为7.887 km，两侧逐渐抬升，使得界面整体呈Ⅴ字型，但两侧存在局部的起伏变化.将界面之上的密度层剖分为17×19个垂直并置的棱柱体，棱柱体的尺寸为1 km×1 km，给定界面上下的密度差为-0.40×103 kg·m-3，计算其引起的重力异常，重力计算点位于每个柱体顶面中心点上，共17×19个计算点.与上文类似，我们将模型在x方向各向外延伸5个点距、y方向各向外延伸10个点距以消除或减弱边部效应.加入均值为0、标准差为0.1×10-5 m·s-2的高斯白噪声后，得到理论重力异常如图 12a所示.
|
图 12 三维Ⅴ字型密度界面模型及Lp-范数约束反演结果 (b)为密度界面模型，其引起的重力异常加入高斯白噪声后如(a)所示，单位为10-5 m·s-2，(c)为Lp-范数约束反演结果，红色圆点为已知深度约束点的位置. Fig. 12 Test using synthetic data produced by a 3D Ⅴ-shape density interface (a) Noise-corrupted gravity anomalies produced by the density interface in (b). The unit is 10-5 m·s-2. (b) Perspective view of the assumed interface. (c) Perspective view of the estimated interface using Lp-norm regularization. Red dots denote the locations of three known depths. |
采用3个已知深度点(图 12c中红色圆点)作为约束，给定p=5，阈值ε=10-4，正则化参数λh=1、λm=0.1，密度差为-0.40×103 kg·m-3，模型剖面参数与正演一致进行反演，结果如图 12c所示.反演结果整体呈Ⅴ字型，展现了理论模型的起伏形态.反演误差如图 13所示，大部分区域的误差小于0.5 km，最大误差为1.152 km，位于x=8 km、y=3 km处.经统计，平均误差为0.013 km，均方差为0.272 km.可见本文反演方法也适用于三维Ⅴ字型密度界面的形态的刻画.
|
图 13 三维Ⅴ字型密度界面模型Lp-范数约束反演结果误差 浅红色区域表示反演结果浅于理论模型，浅蓝色区域表示反演结果深于理论模型. Fig. 13 Perspective view of error between the estimated interface using Lp-norm regularization and the assumed interface The light red area represents the estimated interface is shallower than the assumed one, and the light blue area represents the estimated interface is deeper. |
三维密度界面反演方法中，直接迭代法(Bott, 1960)和Parker-Oldenburg法(Oldenburg, 1974)是研究和使用较多的方法，因此本文也用这两种方法对图 12a所示的重力数据进行反演，结果分别如图 14a和14b所示.两种方法的反演结果均为光滑界面，与理论模型(图 12b)的特征差别较大.此外就图 14a和14b的特征来看，这两种方法的反演结果也有一定的差别.直接迭代法反演结果在形态上稍窄一些，最大深度为6.641 km；而Parker-Oldenburg法反演结果更为宽缓和光滑，最大深度为5.199 km.图 14的试验结果表明，此类方法能较好的应用于光滑密度界面的反演并且具有计算速度快的优势，可应用于大区域、大数据量的反演，如Parker-Oldenburg法已广泛的应用于区域莫霍面深度反演(王万银和潘作枢, 1993; Prasanna et al., 2013; 张恩会等, 2015).
|
图 14 三维Ⅴ字型密度界面模型反演结果.其中(a)为直接迭代法反演结果，(b)为Parker-Oldenburg法反演结果 Fig. 14 Perspective view of the estimated interface using (a) Bott′s method and (b) Parker-Oldenburg method |
马里亚纳海沟是世界上最著名的、最典型的海沟之一，其南段号称“挑战者深渊”(Challenger Deep)(Fujioka et al., 2002; Gvirtzman and Stern, 2004)，最大水深超过10000 m，是地球最深的地方(Gvirtzman and Stern, 2004; Nakanishi and Hashimoto, 2011; 刘方兰和曲佳, 2013).马里亚纳海沟是由太平洋板块向西俯冲到菲律宾海板块之下形成的(Hall, 2002; 刘鑫等, 2017)，莫霍面形态较为复杂(邢健等, 2016; 胡立天等, 2016)，其为不连续的界面.若不考虑太平洋板块俯冲片的深部形态，则挑战者深渊之下的莫霍面连续且在剖面上应呈现Ⅴ字型.为检验本文反演方法的有效性及实用性，我们将其用于挑战者深渊及其邻区的莫霍面深度反演.
根据马里亚纳海沟及其周缘区域重力数据及地形数据质量评价的结果(邢健等, 2016; 胡立天等, 2016)，本文选择来自加利福尼亚圣迭戈分校和Scripps海洋研究所提供的覆盖全球最新的卫星重力数据S & S Global Anomaly V24.1(http://topex.ucsd.edu/cgi-bin/get_data.cgi)(图 15)为挑战者深渊及邻区重力基础数据，该数据网度为1′×1′.
|
图 15 挑战者深渊及其邻区自由空间重力异常图 Fig. 15 Free-air gravity anomalies in the Challenger Deep and adjacent regions |
海域卫星重力数据是海平面以下所有密度不均匀的综合反映.为反演莫霍面深度，需先从卫星重力异常中消除其他地质体的影响而得到“单纯”由莫霍面引起的重力异常.将胡立天等(2016)的成果进行简化，则马里亚纳海沟及其周缘区域主要存在4个密度层：海水(1.03×103 kg·m-3)、沉积层(2.30×103 kg·m-3)、地壳(2.70×103 kg·m-3)和上地幔(3.20×103 kg·m-3)，其形成了3个密度界面：海水与沉积层(-1.27 ×103 kg·m-3)、沉积层与地壳(-0.40×103 kg·m-3)及地壳与上地幔(-0.5×103 kg·m-3).本文采用逐层剥离的方法，首先利用加利福尼亚圣迭戈分校和Scripps海洋研究所提供的水深数据SRTM30_PLUS V11(http://topex.ucsd.edu/cgi-bin/get_srtm30.cgi)(数据网度为0.5′×0.5′，如图 16所示)正演消除海水的重力影响，然后利用来源于NOAA National Geophysical Data Center的沉积层厚度数据(http://www.ngdc.noaa.gov/mgg/sedthick，如图 17所示)正演消除沉积层重力影响，最终得到了挑战者深渊及其邻区的莫霍面引起的重力异常(图 18).
|
图 16 挑战者深渊及其邻区海底地形图 Fig. 16 Submarine topography in the Challenger Deep and adjacent regions |
|
图 17 挑战者深渊及其邻区沉积物厚度图 Fig. 17 Sedimentary thickness in the Challenger Deep and adjacent regions |
|
图 18 挑战者深渊及其邻区莫霍面重力异常图 Fig. 18 Gravity anomalies caused by the Moho in the Challenger Deep and adjacent regions |
以CRUST 1.0提供的莫霍面深度点(http://igppweb.ucsd.edu/~gabi/crust1.html)作为约束(图 19中白色五角星)，反演参数与模型试验中一致，得到挑战者深渊及其邻区莫霍面深度如图 19所示.在研究区中部存在一条NEE-近EW向的莫霍面凹陷带，深度为18~20 km，其清晰地反映了挑战者深渊之下的莫霍面特征.该凹陷带以南区域莫霍面深度为8~12 km，反映了西太平洋板块洋壳的性质.研究区西北角莫霍面深度较大，最大深度超过22 km，为马里亚纳岛弧的反映.
|
图 19 挑战者深渊及其邻区莫霍面深度图 白色五角星为CRUST 1.0提供的深度约束点，黑色实线为图 20断面位置. Fig. 19 Moho depth in the Challenger Deep and adjacent regions White stars denote the locations of known depths of the Moho derived from the CRUST 1.0. The black solid line marks the position of the cross section in Fig. 20. |
研究区莫霍面的断面形态如图 20所示.挑战者深渊海底地形(图 20中蓝色实线)在剖面上呈Ⅴ字型，显示了海沟的几何形态.用本文方法反演的莫霍面(图 20中红色实线)在挑战者深渊附近最深处略显光滑，整体近似地呈Ⅴ字型，而利用Parker-Oldenburg法反演的莫霍面(图 20中白色虚线)更为光滑，且整体起伏稍小.海沟深度最大位置与其下莫霍面深度最大位置(图 20中的黑色虚线)在空间上的连线可以近似的反映出板块俯冲片的形态(图 20中洋红色半透明曲面)，从而可近似的计算出俯冲板片倾角.根据本文反演结果计算出太平洋板块俯冲倾角大致为24°~28°之间，绝大部分区域倾角为25°，可见这一区域板块俯冲倾角相对较大.一般地，俯冲倾角与板块聚敛速率成反比(Luyendyk, 1970)，因此本文反演结果表明，太平洋板块与菲律宾板块汇聚速度相对较慢，基于其他资料的研究成果(邢健等, 2016)亦证实了这一特征.
|
图 20 挑战者深渊及其邻区海底地形及莫霍面起伏断面图(断面位置如图 19中黑色实线所示) 图中上部分界面为海底地形，下部分界面为本文方法反演得到的莫霍面.蓝色实线、红色实线和白色虚线分别为沿断面提取的海底地形、本文方法反演的莫霍面以及Parker-Oldenburg法反演的莫霍面.三维地形图和莫霍面图中黑色虚线分别表示了挑战者深渊地形及莫霍面最大深度连线，将其在空间上连接起来，可近似的表示太平洋板块俯冲片的形态(洋红色半透明区域). Fig. 20 Submarine topography and Moho relief along the cross section in Fig. 19 in the Challenger Deep and adjacent regions The upper interface is the submarine topography and the lower is the inverted Moho relief using the Lp-norm regularization method. The blue solid line, red solid line and white dashed line are the submarine topography, Moho relief using the Lp-norm regularization method and Moho relief using the Parker-Oldenburg method, respectively. The black dashed lines in the two interfaces mark the positions of the maximum depths of submarine topography and Moho, respectively. They, after linking in the space, can approximately delineate the morphology of Pacific plate subduction slab presented with the magenta translucent area. |
本文研究了Ⅴ字型密度界面的正则化反演方法.该方法建立了Lp-范数形式的模型约束函数，并与重力数据误差函数和已知深度约束函数结合形成了反演的目标函数，然后利用非线性共轭梯度法求解目标函数的极小化问题.二维简单模型测试结果表明，在利用基于Lp-范数的正则化方法反演Ⅴ字型密度界面时，p取5、正则化参数λh取1、λm取0.1~0.5、阈值ε取值小于0.01效果较好.二维复杂密度界面和三维密度界面测试以及挑战者深渊及邻区莫霍面反演结果验证了该方法实用性较好，能呈现板块俯冲作用引起的Ⅴ字型莫霍面特征.另外，反演结果表明研究区莫霍面起伏较大，反映了区内构造复杂，地壳性质多变.挑战者深渊及其下莫霍面反映了太平洋板块与菲律宾板块汇聚速度相对较慢.
二维和三维模型试算结果均能很好地刻画界面的Ⅴ字型特征，但挑战者深渊莫霍面反演结果却在最深处稍显光滑，并非严格的Ⅴ字型.推测其原因是重力数据比例尺较小、数据点较稀不足以精细地刻画莫霍面的形态，但与Parker-Oldenburg法相比，其结果更符合地质特征.当然这一区域由于构造复杂，海沟不同位置密度变化较大，其也是引起误差的原因.
需要提及的是，本文研究的非线性反演方法计算量大、速度较慢(效率远低于Parker-Oldenburg法)，并不适合大范围密度界面的反演，可考虑采用并行计算或其他措施提高计算效率.另外，从实际使用的角度考虑，可先采用Parker-Oldenburg法等其他速度较快的方法进行大区域密度界面的反演，之后将本文方法用于重点目标区密度界面起伏特征的精细研究，以提高地质-地球物理解释的精度.
致谢 感谢论文评审专家提出的宝贵意见，感谢本文编辑对文字的加工和修改.
Barbosa V C F, Silva J B C, Medeiros W E. 1997. Gravity inversion of basement relief using approximate equality constraints on depths. Geophysics, 62(6): 1745-1757. DOI:10.1190/1.1444275 |
Barbosa V C F, Silva J B C, Medeiros W E. 1999. Gravity inversion of a discontinuous relief stabilized by weighted smoothness constraints on depth. Geophysics, 64(5): 1429-1437. DOI:10.1190/1.1444647 |
Bott M H P. 1960. The use of rapid digital computing methods for direct gravity interpretation of sedimentary basins. Geophys. J. Roy. Astr. Soc., 3(1): 63-67. DOI:10.1111/j.1365-246X.1960.tb00065.x |
Chakravarthi V, Sundararajan N. 2004. Automatic 3-D gravity modeling of sedimentary basins with density contrast varying parabolically with depth. Comput. Geosci., 30(6): 601-607. DOI:10.1016/j.cageo.2004.03.014 |
Chakravarthi V, Sundararajan N. 2007a. INV2P5DSB-A code for gravity inversion of 2.5-D sedimentary basins using depth dependent density. Comput. Geosci., 33(4): 449-456. DOI:10.1016/j.cageo.2006.06.010 |
Chakravarthi V, Sundararajan N. 2007b. 3D gravity inversion of basement relief-A depth-dependent density approach. Geophysics, 72(2): 123-132. DOI:10.1190/1.2431634 |
Chakravarthi V, Kumar M P, Ramamma B, et al. 2016. Automatic gravity modeling of sedimentary basins by means of polygonal source geometry and exponential density contrast variation:Two space domain based algorithms. J. Appl. Geophys., 124: 54-61. DOI:10.1016/j.jappgeo.2015.11.007 |
Cordell L, Henderson R G. 1968. Iterative three-dimensional solution of gravity anomaly data using a digital computer. Geophysics, 33(4): 596-601. DOI:10.1190/1.1439955 |
Ekblom H. 1973. Calculation of linear best Lp-approximations. BIT Numerical Mathematics, 13(3): 292-300. DOI:10.1007/BF01951940 |
Feng J, Meng X H, Chen Z X, et al. 2014. The investigation and application of three-dimensional density interface. Chinese J. Geophys. (in Chinese), 57(1): 287-294. DOI:10.6038/cjg20140124 |
Feng R, Yan H F, Zhang R S. 1986. The rapid inversion of 3-D potential field and program design. Acta Geologica Sinica (in Chinese), 60(4): 390-403. |
Feng X L, Wang W Y, Liu F Q, et al. 2014. 2D gravity inversion of basement relief of rift basin based on a dual interface model. Chinese J. Geophys. (in Chinese), 57(6): 1934-1945. DOI:10.6038/cjg20140624 |
Feng X L. 2016. Study on the 3D gravity inversion of density interface relief in spatial domain [Ph. D. thesis](in Chinese). Xi'an: Chang'an University.
|
Fujioka K, Okino K, Kanamatsu T, et al. 2002. Morphology and origin of the Challenger Deep in the Southern Mariana Trench. Geophys. Res. Lett., 29(10): 1372. DOI:10.1029/2001GL013595 |
Gvirtzman Z, Stern R J. 2004. Bathymetry of Mariana trench-arc system and formation of the Challenger Deep as a consequence of weak plate coupling. Tectonics, 23(2): TC2011. DOI:10.1029/2003TC001581 |
Hall R. 2002. Cenozoic geological and plate tectonic evolution of SE Asia and the SW Pacific:Computer-based reconstructions, model and animations. Journal of Asian Earth Sciences, 20(4): 353-431. DOI:10.1016/S1367-9120(01)00069-4 |
Hu L T, Hao T Y, Xing J, et al. 2016. The Moho depth in the China Sea-West Pacific and its geological implications. Chinese J. Geophys. (in Chinese), 59(3): 871-883. DOI:10.6038/cjg20160310 |
Huber P J. 1964. Robust estimation of a location parameter. Ann. Math. Stat., 35(1): 73-101. |
Leão J W D, Menezes P T L, Beltro J F, et al. 1996. Gravity inversion of basement relief constrained by the knowledge of depth at isolated points. Geophysics, 61(6): 1702-1714. DOI:10.1190/1.1444088 |
Lima W A, Martins C M, Silva J B C, et al. 2011. Total variation regularization for depth-to-basement estimate:Part 2-Physicogeologic meaning and comparisons with previous inversion methods. Geophysics, 76(1): 113-120. DOI:10.1190/1.3524547 |
Lin Z M, Yang M. 1985. A computer method for gravity interpretation of two-dimensional density contrast interface with some known depths. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese), 28(3): 311-321. |
Liu F L, Qu J. 2013. Seafloor topography and bathymetric survey of the Challenger Deep of Mariana Trench. Marine Geology Frontiers (in Chinese), 29(4): 7-11. DOI:10.16028/j.1009-2722.2013.04.009 |
Liu X, Li S Z, Zhao S J, et al. 2017. Structure of the Mariana subduction system. Earth Science Frontiers (in Chinese), 24(4): 329-340. DOI:10.13745/j.esf.yx.2017-3-16 |
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. Chinese J. Geophys. (in Chinese), 30(2): 186-196. |
Luyendyk B P. 1970. Dips of downgoing lithospheric plates beneath island arcs. Geological Society of America Bulletin, 81(11): 3411-3416. DOI:10.1130/0016-7606(1970)812.0.CO;2 |
Martins C M, Barbosa V C F, Silva J B C. 2010. Simultaneous 3D depth-to-basement and density-contrast estimates using gravity data and depth control at few points. Geophysics, 75(3): 121-128. DOI:10.1190/1.3380225 |
Martins C M, Lima W A, Barbosa V C F, et al. 2011. Total variation regularization for depth-to-basement estimate:Part 1-Mathematical details and applications. Geophysics, 76(1): I1-I12. DOI:10.1190/1.3524286 |
Nakanishi M, Hashimoto J. 2011. A precise bathymetric map of the world's deepest seafloor, Challenger Deep in the Mariana Trench. Mar. Geophys. Res., 32(4): 455-463. DOI:10.1007/s11001-011-9134-0 |
Newman G A, Alumbaugh D L. 2000. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophys. J. Int., 140(2): 410-424. DOI:10.1046/j.1365-246x.2000.00007.x |
Oldenburg D W. 1974. The inversion and interpretation of gravity anomalies. Geophysics, 39(4): 526-536. DOI:10.1190/1.1440444 |
Pallero J L G, Fernández-Martínez J L, Bonvalot S, et al. 2015. Gravity inversion and uncertainty assessment of basement relief via particle swarm optimization. J. Appl. Geophys., 116: 180-191. DOI:10.1016/j.jappgeo.2015.03.008 |
Pallero J L G, Fernández-Martínez J L, Bonvalot S, et al. 2017. 3D gravity inversion and uncertainty assessment of basement relief via particle swarm optimization. J. Appl. Geophys., 139: 338-350. DOI:10.1016/j.jappgeo.2017.02.004 |
Pilkington M, Crossley D J. 1986. Determination of crustal interface topography from potential fields. Geophysics, 51(6): 1277-1284. DOI:10.1190/1.1442180 |
Polak E, Ribière G. 1969. Note sur la convergence de méthodes de directions conjuguées. Revue Fransçaise d'Informatique et de Recherche Opérationnelle, 16: 35-43. |
Prasanna H M I, Chen W, Íz H B. 2013. High resolution local Moho determination using gravity inversion:A case study in Sri Lanka. Journal of Asian Earth Sciences, 74: 62-70. DOI:10.1016/j.jseaes.2013.06.005 |
Shin Y H, Choi K S, Xu H Z. 2006. Three-dimensional forward and inverse models for gravity fields based on the Fast Fourier Transform. Comput. Geosci., 32(6): 727-738. DOI:10.1016/j.cageo.2005.10.002 |
Silva J B C, Oliveira A S, Barbosa V C F. 2010. Gravity inversion of 2D basement relief using entropic regularization. Geophysics, 75(3): I29-I35. DOI:10.1190/1.3374358 |
Silva J B C, Santos D F, Gomes K P. 2014. Fast gravity inversion of basement relief. Geophysics, 79(5): G79-G91. DOI:10.1190/geo2014-0024.1 |
Smith W H F, Sandwell D T. 1997. Global sea floor topography from satellite altimetry and ship depth soundings. Science, 277(5334): 1956-1962. DOI:10.1126/science.277.5334.1956 |
Sun J J, Li Y G. 2014. Adaptive Lp inversion for simultaneous recovery of both blocky and smooth features in a geophysical model. Geophys. J. Int., 197(2): 882-899. DOI:10.1093/gji/ggu067 |
Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. Washington DC: V. H. Winston & Sons.
|
Wang W Y, Pan Z S. 1993. Fast solution of forward and inverse problems for gravity field in a dual interface model. Geophysical Prospecting for Petroleum (in Chinese), 32(2): 81-87. |
Xing J, Hao T Y, Hu L T, et al. 2016. Characteristics of the Japan and IBM subduction zones:Evidence from gravity and distribution of earthquake sources. Chinese J. Geophys. (in Chinese), 59(1): 116-140. DOI:10.6038/cjg20160110 |
Xing J, Hao T Y, Xu Y, et al. 2016. Integration of geophysical constraints for multilayer geometry refinements in 2.5D gravity inversion. Geophysics, 81(5): G95-G106,. DOI:10.1190/GEO2015-0565.1 |
Zhang E H, Shi L, Li Y H, et al. 2015. 3D interface inversion of gravity data in the frequency domain using a parabolic density-depth function and the application in Sichuan-Yunnan region. Chinese J. Geophys. (in Chinese), 58(2): 556-565. DOI:10.6038/cjg20150218 |
Zhou X B. 2013. Gravity inversion of 2D bedrock topography for heterogeneous sedimentary basins based on line integral and maximum difference reduction methods. Geophysical Prospecting, 61(1): 220-234. DOI:10.1111/j.1365-2478.2011.01046.x |
冯娟, 孟小红, 陈召曦, 等. 2014. 三维密度界面的正反演研究和应用. 地球物理学报, 57(1): 287-294. DOI:10.6038/cjg20140124 |
冯锐, 严惠芬, 张若水. 1986. 三维位场的快速反演方法及程序设计. 地质学报, 60(4): 390-403. |
冯旭亮, 王万银, 刘富强, 等. 2014. 裂陷盆地基底双界面模式二维重力反演. 地球物理学报, 57(6): 1934-1945. DOI:10.6038/cjg20140624 |
冯旭亮. 2016.空间域密度界面三维反演方法研究[博士论文].西安: 长安大学地质工程与测绘学院.
|
胡立天, 郝天珧, 邢健, 等. 2016. 中国海-西太平洋莫霍面深度分布特征及其地质意义. 地球物理学报, 59(3): 871-883. DOI:10.6038/cjg20160310 |
林振民, 阳明. 1985. 具有已知深度点的条件下解二度单一密度界面反问题的方法. 地球物理学报, 28(3): 311-321. DOI:10.3321/j.issn:0001-5733.1985.03.009 |
刘方兰, 曲佳. 2013. 马里亚纳海沟水深探测及"挑战者深渊"海底地形特征. 海洋地质前沿, 29(4): 7-11. DOI:10.16028/j.1009-2722.2013.04.009 |
刘鑫, 李三忠, 赵淑娟, 等. 2017. 马里亚纳俯冲系统的构造特征. 地学前缘, 24(4): 329-340. DOI:10.13745/j.esf.yx.2017-3-16 |
刘元龙, 郑建昌, 武传珍. 1987. 利用重力资料反演三维密度界面的质面系数法. 地球物理学报, 30(2): 186-196. DOI:10.3321/j.issn:0001-5733.1987.02.009 |
王万银, 潘作枢. 1993. 双界面模型重力场快速正反演问题. 石油物探, 32(2): 81-87. |
邢健, 郝天珧, 胡立天, 等. 2016. 对日本俯冲带与IBM俯冲带俯冲特征的地球物理研究:来自重力与震源分布数据的启示. 地球物理学报, 59(1): 116-140. DOI:10.6038/cjg20160110 |
张恩会, 石磊, 李永华, 等. 2015. 基于抛物线密度模型的频率域三维界面反演及其在川滇地区的应用. 地球物理学报, 58(2): 556-565. DOI:10.6038/cjg20150218 |