地球物理学报  2017, Vol. 60 Issue (4): 1398-1410   PDF    
三维复杂地壳结构非线性走时反演
俞贵平1,2, 徐涛1,3 , 张明辉1,2, 白志明1, 刘有山1, 武澄泷1,2, 滕吉文1     
1. 中国科学院地质与地球物理研究所, 岩石圈演化国家重点实验室, 北京 100029;
2. 中国科学院大学, 北京 100049;
3. 中国科学院青藏高原地球科学卓越创新中心, 北京 100101
摘要: 中国大陆中西部乃至全球造山带普遍具有复杂地壳结构.随着矿产资源勘探和深部探测研究的深入,探测造山带及盆山耦合区下方地壳精细结构正逐渐成为当前面临的巨大挑战.人工源深地震测深方法正越来越清晰地揭示出不同构造域地壳速度结构的基本特征,然而传统的层状结构模型参数化方法难以准确描述复杂地质模型,通常情况下多忽略速度结构的精细间断面且采用层边界平滑处理,难以满足地壳精细结构成像的发展要求.针对上述困难,本文采用最近发展的块状结构建模方案构建三维复杂地壳模型,基于逐段迭代射线追踪正演走时计算方法,推导了走时对三角形界面深度以及网格速度的偏导数,开展了非线性共轭梯度走时反演方法研究.发展了利用直达波和反射波等多震相走时数据对界面深度和网格速度的多参数联合反演方法,并引入不同种类震相数据的权系数和不同类型参数偏导数归一化的方法.数值算例表明,基于块状结构的非线性共轭梯度走时反演方法适用于复杂地壳结构模型,在利用人工源走时数据反演复杂地壳精细结构领域具有良好的应用前景.
关键词: 地壳结构      走时反演      块状建模      射线追踪      非线性共轭梯度     
Nonlinear travel-time inversion for 3-D complex crustal velocity structure
YU Gui-Ping1,2, XU Tao1,3, ZHANG Ming-Hui1,2, BAI Zhi-Ming1, LIU You-Shan1, WU Cheng-Long1,2, TENG Ji-Wen1     
1. State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. CAS Center for Excellence in Tibetan Plateau Earth Science, Beijing 100101, China
Abstract: Complex crustal structures generally characterize global orogenic belts and Midwest China. With the deepening of mineral resources exploration and increasing detection of the deep earth, it is becoming a great challenge to use new methods to probe fine crustal structure beneath orogenic belts and basin-mountain coupling regions. Basic characteristics of crustal structure of different tectonic domains are becoming clearer and clearer in deep seismic soundings. However, it is difficult for the traditional layered structure modeling method, in which layer boundaries are smoothing and fine velocity discontinuities are often ignored, to describe complex geological models, making it difficult to meet the development requirements of fine structure imaging of crust. In view of the above difficulties, based on a newly developed block modeling scheme to describe three-dimensional complex crustal models and a corresponding segmentally iterative ray-tracing (SIRT) method, we derive the travel-time partial derivatives of triangular interface depth and grid velocity, and develop a 3-D nonlinear conjugate gradient travel-time inversion method. Block modeling scheme is able to construct any complex geological models in theory and can be used to efficiently build initial models by combining various priori velocity and interface information. In the inversion process, PRP (Polak-Ribière-Polyak) type of the conjugate gradient method is used to solve the constrained damping least squares problem. We develop the joint inversion of grid velocity and interface depth based on multi seismic phases like direct waves and reflected waves, and make a great improvement of the inversion resolution compared with the traditional method which is based on single phase. To improve the convergence accuracy of inversion results, strategies like different weighting factors for different seismic phases and normalization of travel-time partial derivatives of different parameters are introduced in the joint inversion process. Numerical examples are given and a special layered cross-cutting model is used to simply simulate the application form of our method for fine structure imaging of crust, showing that nonlinear conjugate gradient travel-time inversion method based on block modeling is suitable for complex crustal models, which has good application prospects in the field of fine crustal structure imaging based on travel-time data of artificial seismic sources.
Key words: Crustal structure      Travel-time inversion      Block modeling      Ray tracing      Nonlinear conjugate gradient     
1 引言

全球造山带及中国大陆西部盆山结合带普遍具有丰富的资源勘探前景,同时这些构造带域由于板块拼贴和碰撞作用往往具有复杂地壳结构.自20世纪50年代以来,人们将地震探测技术广泛应用于研究洋脊、造山带、克拉通和盆山结合带等不同构造区域的壳幔几何结构和物理属性,越来越清晰地揭示出大型速度间断面和地质体的展布形态和基本特征 (Thybo et al., 2013; Prodehl et al., 2013),建立起更加精细的区域地震学参考模型 (Teng et al., 2013, 2014; Zhang et al., 2014),为认识全球壳幔分界面 (Moho) 的展布形态和典型构造域地壳结构的基本特征提供了重要约束.

人工源宽角反射/折射地震探测 (或深地震测深) 是研究地壳速度结构的重要途径.自20世纪80年代以来,国内外专家已发展了一系列宽角反射/折射资料处理方法,并使之适应于复杂非均匀层状介质,旨在解决运动学 (走时) 和动力学射线 (振幅) 计算问题,以获得地壳速度结构等.当前,人工源宽角地震观测以二维剖面为主.二维地震剖面资料的处理方法主要包括以Seis8X系列软件为代表的正演走时拟合 (Červený et al., 1988; Červený, 2001;徐涛等, 2014, Xu et al., 2015) 和以Rayinvr为代表的走时反演 (Zelt and Smith, 1992; Zhang et al., 2013).三维宽角地震资料处理方法,主要包括利用初至波走时数据,基于离散介质程函方程数值解并反演上地壳结构的有限差分成像方法 (Vidale, 1990; Hole, 1992; Zelt and Barton, 1998; 赵烽帆等, 2014a),以及基于连续介质射线理论的走时反演方法 (Gajewski and Pšenčík, 1988; Rawlinson et al., 2001; Červený, 2001; Pšenčík and Červený, 2002).后者以Anray95软件为代表,能适用于三维各向同性或各向异性介质 (Gajewski and Pšenčík, 1988) 以及复杂层状结构和块状结构的地震射线方法研究 (Červený and De Castro, 1993; Pšenčík and Červený, 2002).例如,CELEBRATION 2000实验 (Guterch et al., 2003),在中欧地区采集了多条二维观测剖面的人工源地震资料 (środa et al., 2006; Hrubcova et al., 2008),同时开展了部分三维壳幔结构成像研究 (Malinowski et al., 2008, 2009),并在反演过程中通过逐步增加偏移距范围来约束模型的不同深度和逐步减小模型划分单元来获得更优解,以及开展直达波、反射波震相走时联合反演等.

人工源宽角反射/折射地震探测资料的处理方法主要基于射线理论,模型参数化方法多采用层状结构模型 (Červený et al., 1988; Červený, 2001; Gajewski and Pšenčík, 1988; Zelt and Smith, 1992; Rawlinson et al., 2001; 周龙泉等, 2006).例如,Rawlinson等 (2001)提出的三维宽角地震资料走时反演方法,虽然可以实现速度和界面的联合反演,但却仍然局限于相对简单的层状结构模型.层状结构模型简单,正演射线追踪计算以及走时反演方便,在人工源地震资料处理过程中占据重要的位置.但是层状结构模型参数化方法往往忽略壳内速度结构的精细间断面且采用层边界平滑处理,难以准确地描述复杂地质模型,无法满足地壳精细结构成像的发展要求.因此,发展复杂地壳结构的模型参数化方法,以及该基础上的射线正演及走时反演方法对于实现精细地壳结构成像有重要意义.

针对上述问题,我们提出了块状结构建模方案,并采用三角形网格化描述地质界面,理论上可以描述任意复杂地质模型 (徐涛等, 2004; Xu et al., 2006, 2010, 2014);发展了逐段迭代射线追踪方法 (Xu et al., 2006, 2010, 2014; 李飞等, 2013),从而高效地计算直达波、反射波等震相走时.该方法首先应用在石油地震勘探领域,在复杂地壳结构成像领域有良好的应用前景.

本文在块状模型基础上,发展了利用直达波和反射波等多震相走时对界面深度和网格速度的多参数联合反演方法,通过非线性共轭梯度反演,实现复杂地壳速度结构走时反演.

2 模型参数化 2.1 块状结构模型

针对层状结构很难描述复杂地质体的问题,本文采用块状结构建模方案,即三维地质体被看作是由大小不等、形状各异的地质块体组成的集合体,结合曲面三角形网格来描述地质界面,理论上可以构建任意复杂地质模型 (徐涛等, 2004; Xu et al., 2006, 2008, 2010, 2014; 李飞等, 2013),图 1为块状结构构建伸展盆地中的堑垒构造,模型由18个地质块构成,6676个曲面三角形构建模型界面,三角形由2700个离散点来控制.为了克服三角形面片法向量方向的不连续导致地震射线路径的畸变,我们重新定义三角形顶点的法向量,并由此重定义界面任意点的法向量,使得法向量的方向在整个界面上是光滑连续的 (徐涛等, 2004; Xu et al., 2006, 2010, 2014).

图 1 块状结构构建三维复杂堑垒构造地质模型 (a) 模型由不同的地质块组成;(b) 界面由三角形面片拼接而成. Fig. 1 3-D Complex horst and graben structure model constructed by block modeling (a) Model composed of different blocks; (b) Triangulated interfaces.
2.2 非均匀速度结构分布

块状结构模型实质为分块连续介质,不同块体之间的分界面为速度间断面.根据实际地质情况以及研究的需要,可以在不同块体内分别定义均匀、非均匀、各向异性等不同介质.对于常见的非均匀介质,给定一组矩形离散网格节点,我们通过三线性插值方法来获得任意点的速度 (Zhao et al., 1992, 1994; 李飞等, 2013; Xu et al., 2014).如图 2所示,对于某块体内部的任意点P,首先找到该点所在矩形网格,并由周围8个节点位置处的速度三线性插值获得P点的速度:

图 2 基于矩形网格的三线性速度插值 Fig. 2 Trilinear velocity interpolation based on rectangular grids

(1)

式中,v(i+l, j+m, k+n) 分别代表 8个节点的速度值;xiyjzk分别为矩形速度网格在xyz三个方向的位置坐标.

3 射线追踪

Um和Thurber (1987)提出了伪弯曲法射线追踪,用于解决连续介质中的两点射线追踪问题,但却无法适应存在强速度间断面的情况.我们提出了逐段迭代射线追踪方法,用来追踪存在强速度间断面的模型 (Xu et al., 2006, 2010, 2014).该方法属于弯曲法,在扰动修正路径点的过程中,对落在块体内部或界面上的路径点,采用一阶显式扰动修正的方案,修正公式分别如3.1节和3.2节所示.

3.1 块体内部路径点的修正

我们采用伪弯曲法 (Um and Thurber, 1987) 对位于块体内部的射线路径点进行修正,如图 3所示.对于连续的三个路径点,固定起始点Pk-1和终止点Pk+1,中间点Pk的修正包括方向向量n和距离R两个参数的修正:

图 3 伪弯曲法修正落在介质内部的中间点示意图 Fig. 3 Illustration of mid-point modification of pseudo-bending method when the mid-point isinside a medium

(2)

其中,V为中间路径点Pk处的速度梯度向量,Vmid为起始点和终止点的中间点Pmid处的速度值,V为中点Pmid处的速度梯度向量.

3.2 界面上路径点的修正

采用逐段迭代法对界面上的路径点进行修正,如图 4所示.界面两侧为非均匀速度分布,我们采用一阶显式近似修正公式.对于连续的三个路径点Pk-1PkPk+1,其中Pk是界面上的路径点,Pk-1Pk+1分别为界面两侧的路径点,固定起始点和终止点,则射线路径走时T是中间点的函数,可以近似表示为

图 4 逐段迭代法修正落在界面上的中间点示意图 Fig. 4 Illustration of mid-point modification of segmentally-iterative method when the mid-point is located on an interface

(3)

其中v0v3分别为起始点和终止点处的速度,v1v2为界面中间点两侧的速度.

根据稳定走时原理,修正后的中间点位置处,满足地震波走时对空间位置的偏导数为零,即

(4)

对 (4) 式进行泰勒展开并保留一阶小量,最终可以得到中间路径点的一阶修正量

(5)

其中,Uij的具体表达式详见相关射线追踪方法研究 (李飞等, 2013; Xu et al., 2014).

4 非线性共轭梯度走时反演 4.1 非线性共轭梯度方法

反演思路为采用PRP (Polak-Ribière-Polyak) 型非线性共轭梯度法求解带约束的阻尼最小二乘问题,并建立如下目标函数:

(6)

其中,g(m) 和dobs分别为理论计算的走时数据和实际观测的走时数据;λ为正则化参数,一般为正值并随着迭代次数的增加而减小;mp为先验模型参数,CdCm分别为估计的数据和模型协方差矩阵.相应的Jacobian矩阵和近似Hessian矩阵形式如下:

(7)

其中,矩阵G的维数为 (M1+M2)×(N1+N2),M1M2分别为直达波和反射波震相的射线数,N1N2分别为速度网格和界面顶点的参数个数.

非线性共轭梯度反演流程可以归纳如下:

(1) 令i=0,计算目标函数的梯度γ(i)和下降方向u(i),如果梯度向量‖γ(i)‖ < ε,则迭代结束,其中

(8)

(9)

(2) 计算步长α(i),使得目标函数Ψ(m(i)+α(i)u(i)) 下降,定义v(i)=G(i)u(i),有

(10)

(3) 更新模型参数

(11)

(4) 计算新的梯度向量γ(i+1),并判断是否满足‖γ(i+1)‖ < ε,满足则退出循环,其中

(12)

(5) 更新PRP型共轭梯度参数

(13)

(6) 更新下降方向

(14)

(7) 判断u(i+1) Tγ(i+1)≥0或i=M是否成立 (M为待反演参数的总个数),若其中一项或两项均成立则转步骤 (8),否则转步骤 (9);

(8) 更新下降方向 (也可根据实际情况提前重置搜索方向为负梯度方向)

(15)

(9) i=i+1,转步骤 (2).

4.2 偏导数矩阵

反演过程中,构建Jacobian矩阵需计算射线走时对界面三角形节点位置坐标以及矩形网格节点速度的偏导数 (Zelt and Smith, 1992; Rawlinson et al., 2001; 周龙泉等, 2006; 白超英等, 2011; 黄国娇和白超英, 2013; Li et al., 2014).

(1) 走时对三角形节点深度的偏导数

在界面深度反演时,本文固定界面三角形节点的xy坐标,只反演z坐标一个参数.走时对界面节点深度偏导数一般化为三个部分 (Zelt and Smith, 1992; Rawlinson et al., 2001)

(16)

其中,zj为界面节点深度,hint为路径点处垂直于界面的位移,zint为路径点的深度,如图 5所示.

图 5 走时对三角界面节点深度的偏导数示意图 Fig. 5 Illustration of the calculation of travel-time derivative with respect to reflector depth change

考虑到本研究模型界面采用的是三角网格,具体的偏导数公式 (Li et al., 2014) 可以化为

(17)

其中,Pk-1PkPk+1分别为射线路径上连续的三个路径点,Pk为落在三角形Aj-1AjAj+1上的中间路径点;ulPk点在三角形Aj-1AjAj+1中的面积坐标 (Xu et al., 2006);WzZ轴方向向量,Wt为垂直平面三角形的法向量,WnPk点处重新定义的法向量 (Xu et al., 2006).一般情况下,WnWt方向不一定相同,WkWk+1分别为射线入射方向和出射方向的单位向量.

(2) 走时对网格节点速度的偏导数

沿着射线路径L,走时t为慢度 (速度倒数) 沿射线轨迹的积分,其离散的近似表达式为

(18)

其中,v(k)v(k+1)为离散射线路径段 (直线) 的起始点和终止点处的速度.则走时对网格节点速度的偏导数为

(19)

或者

(20)

假定某路径点P(k)(x, y, z) 的速度值为v(k),周围是8个网格速度点,如图 2所示.根据线性插值速度分布定义 (式 (1)),则式 (21) 中,除射线路径点P(k)周围的8个节点以外,路径点速度v(k)对网格节点速度的偏导数βi(k)为零,即

(21)

其中,l=0,1;m=0,1;n=0,1.

4.3 联合反演中不同类型偏导数间的平衡

开展多震相走时联合反演,可以有效提高反演结果的分辨率 (田有等, 2009; Bai et al., 2010; 白超英等, 2011; Huang et al., 2012; 黄国娇和白超英, 2013).如果各震相走时偏导数权重一致,由于不同震相的走时拾取误差和射线覆盖区域不同,会导致模型空间收敛不均衡,影响整体的收敛精度和速度.因此,需要对不同震相的走时偏导数给以不同的权重 (白超英等, 2011),从而保证迭代收敛的均衡性.本文在直达波和反射波联合反演速度结构时,保持反射波走时偏导数不变,通过调整直达波走时偏导数的权系数 (0~1) 来获得最佳方案.

与多震相走时联合反演相比,速度和界面等多参数联合反演的实现通常是比较困难的.难点主要有两个方面.一是联合反演导致参数增加,反演难度整体上升,且界面反演依赖的反射波射线路径对界面形态非常敏感,易出现射线覆盖的阴影区,多解性强.二是如何实现偏导数矩阵中走时对界面深度偏导数和走时对速度偏导数的平衡,从而确保收敛方向不会一致地收敛到偏导数矩阵中系数大的一方.本文对界面深度偏导数乘以因子ω,从而保证偏导数矩阵中走时关于速度的偏导数和走时关于界面深度的偏导数在总值上相当 (McCaughey and Singh, 1997),其中ω的表达式为:

(22)

式中,Aij表示第i条射线对第j个网格速度的走时偏导数,Bij表示第i条射线对第j个界面顶点深度的走时偏导数,S表示总射线数,N1N2分别表示待反演的速度参数和界面深度参数的个数.

5 数值试验

为了验证本文反演方法的有效性,我们开展了三维体波走时反演数值模拟试验,并采用检测板来检验反演结果的精度.反演过程中增加了模型空间约束:mpmmaxmnewmpmmax,其中,mp为初始模型参数,mnew为迭代更新后的模型参数,Δmmax即为约束值的大小.对射线覆盖较差的区域,增加约束可以保证反演结果的合理性.速度结构反演时,检测板的速度扰动值均为±200 m·s-1,故速度约束值取200 m·s-1;界面深度反演时,理论模型的起伏界面与初始水平界面的深度差范围为±3 km,故深度约束值取3 km.以下各模型的地表面均为水平面,模型尺度、速度网格和界面顶点划分等基本参数详见表 1.

表 1 三维体波走时反演模型参数 Table 1 Model parameters in 3-D body-wave travel-time inversion

在速度检测板恢复实验中,理论速度模型由常梯度速度加正负相间的速度扰动构成,初始速度模型设置为常梯度模型.其中常梯度模型横向速度均匀,纵向速度随深度增加线性递增,各常梯度模型的速度区间详见表 1.值得指出的是,由于地层互切模型Model 3只反演第二层速度结构,故表 1给出的是第二层介质的速度和界面参数,其他模型均为单层模型,各单层的速度或界面参数详见表 1;在Model 1b中,为了提高对异常边缘的精细刻画能力,检测板是密集网格剖分下以多个节点 (2×2×1) 为单位的正负扰动,相应图件做了插值处理,其他模型的检测板均是以单个网格节点 (1×1×1) 为单位的正负扰动,为直观展示反演后检测板的恢复程度,作图时均未做插值处理.

5.1 速度结构反演

(1) 直达波反演

直达波震相通常表现为初至波震相,在地壳结构成像中具有重要地位.为合理检验反演算法的可行性,采用理想固定观测系统 (炮点位于地下z=-45 km,接收器位于地表) 以减少对反演结果的影响,并以Model 1a (图 6b) 和Model 1b (图 6b) 为例,开展了直达波反演速度分辨率测试试验.Model 1a和Model 1b的模型尺度、观测系统和初始模型均相同,但网格划分和检测板形态不同,有关参数详见表 1.

图 6 直达波反演速度结构分辨率测试 (a) 均匀观测系统平面分布图 (炮点数9×9,接收器数8×8),其中蓝色五角星表示炮点,红色三角形表示接收器,炮点位于z=-45 km,接收器位于地表,灰色虚线表示速度网格线; (b) 两炮激发时,初始模型直达波射线路径. Model 1a检测板在z=-22.5 km处的原始切片 (c) 和反演切片 (e),以及在y=25 km处的原始剖面 (d) 和反演剖面 (f). Model 1b检测板在z=-22.5 km处的原始切片 (g) 和反演切片 (i),以及在y=23 km处的原始剖面 (h) 和反演剖面 (j). Fig. 6 Resolution tests for velocity inversion based on direct waves (a) Plane distribution of uniform observation system (shot points 9×9, receivers 8×8). Blue stars denote shot points and red triangles denote receivers. Shot points are located at z=-45 km and receivers are located on the surface. Gray dashed lines denote velocity grid lines; (b) Direct wave ray paths of initial model under two blasting; (c)(e) Initial slice and inversion slice of Model 1a checkerboard at z=-22.5 km; (d)(f) Initial profile and inversion profile of Model 1a checkerboard at y=25 km; (g)(i) Initial slice and inversion slice of Model 1b checkerboard at z=-22.5 km; (h)(j) Initial profile and inversion profile of Model 1b checkerboard at y=23 km.

Model 1a的反演结果如图 6e6f所示.检测板中间恢复良好,由中心向外围逐渐模糊,恢复程度与射线覆盖次数和射线的方位角宽度呈正相关,说明所采用的正、反演方法是有效的.在以单个网格节点为单位的正负扰动检测板模型中,Model 1a的速度网格分辨率已达到最高 (9×9×7),由于射线覆盖有限,继续增加网格点密度,检测板将无法恢复.在实际工作中,为提高对异常边缘的精细刻画能力,通常会进一步加密网格点,并开展以多个网格节点为单位的正负扰动检测板试验.为此,我们设计了以2×2×1个网格节点为单位的正负扰动检测板模型Model 1b,反演结果如图 6i6j所示,检测板整体恢复良好,但收敛速度和精度均有下降.该方案适用于速度变化平稳、异常区块较大的情况,但对于速度变化剧烈的复杂结构,以单个网格节点为单位的正负扰动检测板试验仍然是最可靠的.

上述反问题,通常由于缺少对误差或模型的认识,可能同时存在超定或欠定方面的问题,即表现为混定问题.对于射线覆盖较差的网格边缘,可以考虑单独设置较大的网格来提高反演结果的可靠性和稳定性.网格大小的设计需要综合考虑射线覆盖次数和模型的复杂程度,其中射线覆盖情况决定了反演结果的分辨率.

(2) 多震相联合反演

我们以Model 2a (图 7a) 为例,在直达波或反射波独立反演速度结构的基础上,开展了多震相走时联合反演研究.Model 2a在z=-45 km处有一个水平界面.模型尺度和观测系统见表 1图 7b.类似人工源地震观测,炮点和接收器均位于地表.由于模型和观测系统均是对称的,这里仅以y=75 km处垂直剖面为例对反演结果进行评价.

图 7 直达波和反射波反演速度结构 (c)—(g) 均为y=75 km处的检测板剖面.(a) 两炮激发时,初始模型直达波和反射波射线路径; (b) 均匀观测系统平面分布图 (炮点数9×9,接收器数8×8),其中蓝色五角星表示炮点,红色三角形表示接收器,炮检点均位于地表,灰色虚线表示速度网格线; (c) 原始检测板剖面; (d) 直达波独立反演; (e) 反射波独立反演; (f) 直达波和反射波权重系数相同时的联合反演; (g) 直达波走时对速度的偏导数乘以0.4的权重系数后的联合反演. Fig. 7 Velocity inversion based on direct waves and reflected waves (c)—(g) All the checkerboard profiles are located at y=75 km. (a) Direct wave and reflected wave ray paths of initial model under two blasting; (b) Plane distribution of uniform observation system (shot points 9×9, receivers 8×8). Blue stars denote shot points and red triangles denote receivers. Shot points and receivers are located on the surface. Gray dashed lines denote velocity grid lines; (c) Original checkerboard profile; (d) Velocity inversion based on direct wave; (e) Velocity inversion based on reflected wave; (f) Velocity inversion based on direct wave and reflected wave with the same weight; (g) Velocity inversion based on direct wave and reflected wave, with a weighting factor of 0.4 for direct wave travel-time derivative of velocity.

图 7a为两炮激发时,直达波 (回折波) 和反射波射线路径图,图中直达波集中于近地表,射线密集且各方向交错充分,直达波独立反演剖面如图 7d所示,浅层速度结构恢复良好,与实际射线覆盖情况相对应.图 7e为反射波独立反演剖面,检测板的整体形态得到了恢复,但浅层结构恢复程度较低,这是由于浅层射线的方位角较窄造成的.直达波或反射波独立反演的结果均不够理想,为了提高反演结果的分辨率,我们开展了直达波和反射波联合走时反演试验,结果如图 7f所示.检测板的浅层和深部都得到了恢复,相比于单一震相反演有明显优势.但是,深部周缘的恢复程度并未达到反射波独立反演时的效果,这是由于浅层直达波射线覆盖过密造成的浅层和深部收敛不均衡.

为了平衡直达波 (回折波) 和反射波对模型的约束,经测试决定:对直达波走时偏导数乘以0.4的权重系数,且保持反射波走时偏导数不变.考虑不同震相走时的权重后,联合反演迭代后期收敛依然明显,图 7g为40次迭代后的反演结果,而图 7d—7f对应的反演过程在迭代后期已陷入局部极小,因此这里展示20次迭代的结果.可以看到,调整不同震相的权重系数后,检测板得到了很好的恢复,充分发挥了两种震相联合反演的优势.事实上采用二次反演的组合策略同样可以达到类似的效果,比如以反射波独立反演的结果作为联合反演的初始速度模型,或者基于联合反演结果利用反射震相进行二次反演等,其目的均是为了充分发挥反射波的深部探测优势和直达波的浅层探测优势.

5.2 界面深度反演

我们以起伏界面模型Model 2b (图 8b) 为例,开展了反射波反演界面深度研究.试验采用常梯度速度模型,模型尺度和界面参数见表 1,观测系统和理论界面形态如图 8a8b所示.反演时初始水平界面深度为z=-42 km,Model 2b界面最高点和最低点相差6 km.反演界面与初始水平界面的深度差如图 8d所示,图中界面的大部分区域得到了很好的恢复,只有射线覆盖差的边缘恢复程度较低.由于初始深度设置合理,左右两侧边界未出现明显假象.

图 8 反射波反演界面深度 (a) 均匀观测系统平面分布图 (炮点数3×3,接收器数7×7),其中绿色五角星表示炮点,红色三角形表示接收器,炮检点均位于地表,灰色虚线为界面顶点网格线,蓝色实线为粗略表示的射线路径; (b) 两炮激发时理论模型反射波射线路径; (c) 理论模型界面与初始水平界面的深度差; (d) 反演界面与初始水平界面的深度差; (e) 反演界面与理论模型界面的深度差. Fig. 8 Interface depth inversion based on reflected waves (a) Plane distribution of uniform observation system (shot points 3×3, receivers 7×7). Green stars denote shot points and red triangles denote receivers. Shot points and receivers are located on the surface. Gray dashed lines denote grid lines of interface vertexes; (b) Reflected wave ray paths of real model under two blasting; (c) Depth difference between theoretical model interface and initial level interface; (d) Depth difference between inversion interface and initial level interface; (e) Depth difference between inversion interface and theoretical model interface.
5.3 界面深度和速度结构联合反演

基于反射波独立反演速度或界面的研究,我们以起伏界面模型Model 2c (图 9a) 为例,进一步开展了速度和界面联合反演试验.Model 2c的模型尺度与Model 2a和Model 2b相一致,观测系统同Model 2a (图 7b),仍以y=75 km处垂直剖面为例对反演结果进行分析.由于联合反演的参数更多、物理关系更复杂,相比Model 2a和Model 2b,我们适当降低了Model 2c的速度网格和界面顶点的划分密度,具体参数详见表 1.

图 9 反射波联合反演速度与界面 各检测板剖面均位于y=75 km.(a) 原始检测板剖面; (b) 理论模型界面与初始水平界面的深度差; (c) 速度与界面直接联合反演; (f) 先单独反演界面再单独反演速度; (i) 基于先后独立反演结果的速度与界面直接联合反演; (l) 基于先后独立反演结果的速度与界面深度偏导数归一化后联合反演; (d)(g)(j)(m) 分别为相应的反演界面与初始水平界面的深度差; (e)(h)(k)(n) 分别为相应的反演界面与理论模型界面的深度差.图中黄色实线表示理论模型界面,紫色虚线表示初始界面,红色点线表示联合反演得到的界面. Fig. 9 Joint inversion of velocity and interface based on reflected waves All the checkerboard profiles are located at y=75 km. (a) Original checkerboard profile; (b) Depth difference between theoretical model interface and initial level interface; (c) Direct joint inversion; (f) Independent velocity inversion after independent interface inversion; (i) Direct joint inversion based on independent inversion result; (l) Joint inversion based on independent inversion result after partial derivative normalization; (d)(g)(j)(m) are corresponding depth difference between inversion interface and initial level interface; (e)(h)(k)(n) are corresponding depth difference between inversion interface and theoretical model interface. In these profiles, yellow lines denote theoretical model interface. Purple dashed lines denote initial interface and red dot lines denote interface of joint inversion.

图 9c—9e分别为速度与界面直接联合反演的速度剖面、反演深度界面与初始深度和理论深度界面的差值.速度和界面的基本形态均得到了恢复,但恢复程度较低.为提高联合反演的精度,我们选择先固定速度,反演界面,再以新的界面作为固定界面,独立反演速度,结果如图 9f—9h所示;再以此作为初始模型进行速度与界面的联合反演,并试图提高浅层速度结构的成像精度,但图 9i—9k所示的反演结果并不理想.为此,我们又在界面与速度先后独立反演的基础上,开展了速度与界面深度偏导数归一化后的联合反演成像,结果如图 9l—9n所示.不难发现,近地表的速度结构特征得到了恢复,偏导数归一化效果明显.当然,对复杂模型进行速度和界面联合反演时,不同类型偏导数归一化的作用并不总是明显,迭代过程仍可能陷入局部极小值,此时可以进一步引入优化步长方法,如抛物搜索法 (Vigh et al., 2009).

5.4 地层互切模型速度反演

基于上述单层模型的速度和界面反演试验均取得了明显成效,通过开展类似的直达波 (回折波) 和反射波多震相走时对界面深度和网格速度的多参数联合反演,可以有效获得地壳的基本结构特征.本文以地层互切模型Model 3(图 10a) 为例,进一步开展了复杂块状结构模型走时反演研究,并将浅层速度结构和界面深度作为先验信息来约束深部地质结构成像.观测系统见图 10c图 10a给出了两炮激发时初始模型第二层界面的反射波射线路径.图 10b为初始速度模型剖面,各块体垂向速度均为常梯度,理论速度模型则是在常梯度模型的基础上为第二层加上检测板.模型尺度和第二层速度参数详见表 1.

图 10 地层互切模型速度结构反演 (b)(d)(e) 剖面均位于y=25 km. (a) 两炮激发时,初始模型第二层界面的反射波射线路径; (b) 初始速度模型剖面,其中各块体的速度在垂向上均为常梯度; (c) 均匀观测系统平面分布图 (炮点数9×9,接收器数8×8),其中蓝色五角星表示炮点,红色三角形表示接收器,炮检点均位于地表; (d) 第二层原始检测板剖面; (e) 第二层检测板反演剖面. Fig. 10 Velocity inversion of layered cross-cutting model (b)(d)(e) All the profiles are located at y=25 km. (a) Reflected wave ray paths from the second interface of initial model under two blasting; (b) Initial velocity model profile, in which different blocks′ velocities are all constant gradient in vertical; (c) Plane distribution of uniform observation system (shot points 9×9, receivers 8×8). Blue stars denote shot points and red triangles denote receivers. Shot points and receivers are located on the surface; (d) Original checkerboard profile of the second layer; (e) Inversion checkerboard profile of the second layer.

在测试过程中,我们固定第一层的速度结构和界面深度,以及第二层界面深度,只反演第二层速度结构.第二层界面为水平界面与大型起伏界面的互切组合,层位较深,反演第二层介质参数主要是利用反射震相,由于界面左右两侧倾角较大,在现有模型尺度下无法观测到反射震相,因此两侧速度结构理论上无法得到有效约束.图 10d为第二层速度原始检测板剖面,图 10e则为相应的反演剖面.经过反演,检测板中间部分恢复良好,两侧模糊,反演结果与实际观测相匹配.

6 讨论和结论

在块状模型基础上,本文推导了走时对三角网格界面顶点深度以及网格速度的偏导数,开展了非线性共轭梯度走时反演方法研究.发展了利用直达波和反射波等多震相走时数据对界面深度和网格速度的多参数联合反演方法,并引入了不同种类震相数据的权系数和不同类型参数偏导数归一化的方法.数值算例表明,联合反演时针对不同震相走时数据采用不同权系数及对不同类型参数偏导数进行归一化的方法,改善了反演结果的收敛精度;基于块状结构的非线性共轭梯度走时反演方法适用于复杂地壳结构模型,在利用人工源宽角地震走时数据反演地壳精细结构领域具有良好的应用前景.不仅如此,基于复杂块状模型的多震相、多参数联合走时反演方法同样可以应用在远震/近震体波走时层析成像领域 (Zhao et al., 1992, 1994; 赵烽帆等, 2014b),与Zhao等 (1992, 1994) 的方法相比,模型构建方法更加灵活,在刻画和约束典型构造域地质结构方面具有一定的优势.

谨以此文纪念中国科学院地质与地球物理研究所张忠杰研究员 (1964—2013).

参考文献
Bai C Y, Huang G J, Zhao R. 2010. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications. Geophys. J. Int., 183(3): 1596-1612. DOI:10.1111/j.1365-246X.2010.04817.x
Bai C Y, Huang G J, Li Z S. 2011. Simultaneous inversion combining multiple-phase travel times within 3D complex layered media. Chinese J. Geophys. , 54(1): 182-192. DOI:10.3969/j.issn.0001-5733.2011.01.019
Červený V, Klimeš L, Pšenčík I. 1988. Complete seismic-ray tracing in three-dimensional structures.//Doornbos D J. Seismological Algorithms. New York:Academic Press, 89-168.
Červený V, De Castro M A. 1993. Application of dynamic ray tracing in the 3-D inversion of seismic-reflection data. Geophys. J. Int., 113(3): 776-779. DOI:10.1111/gji.1993.113.issue-3
Červený V. Seismic Ray Theory. Cambridge: Cambridge University Press, 2001.
Gajewski D, Pšenčík I. 1988. Ray synthetic seismograms for a 3-D anisotropic lithospheric structure. Phys. Earth Planet. Interiors, 51(1-3): 1-23.
Guterch A, Grd M, Keller G R, et al. 2003. CELEBRATION 2000 seismic experiment. Stud. Geophys. Geod., 47(3): 659-669. DOI:10.1023/A:1024728005301
Huang G J, Bai C Y, Zhu D L, et al. 2012. 2D/3D seismic simultaneous inversion for the velocity and interface geometry using multiple classes of arrivals. Bull. Seismol. Soc. Am., 102(2): 790-801. DOI:10.1785/0120110155
Huang G J, Bai C Y. 2013. Simultaneous inversion of three model parameters with multiple classes of arrival times. Chinese J. Geophys. , 56(12): 4215-4225. DOI:10.6038/cjg20131224
Li F, Xu T, Wu Z B, et al. 2013. Segmentally iterative ray tracing in 3-D heterogeneous geological models. Chinese J. Geophys. , 56(10): 3514-3522. DOI:10.6038/cjg20131026
Li F, Xu T, Zhang M H, et al. 2014. Seismic traveltime inversion of 3D velocity model with triangulated interfaces. Earthq. Sci., 27(2): 127-136. DOI:10.1007/s11589-013-0025-0
Malinowski M, Grad M, Guterch A. 2008. Three-dimensional seismic modelling of the crustal structure between East European Craton and the Carpathians in SE Poland based on CELEBRATION 2000 data. Geophys. J. Int., 173(2): 546-565. DOI:10.1111/gji.2008.173.issue-2
Malinowski M, S'roda P, Grad M, et al. 2009. Testing robust inversion strategies for three-dimensional Moho topography based on CELEBRATION 2000 data. Geophys. J. Int., 179(2): 1093-1104. DOI:10.1111/gji.2009.179.issue-2
McCaughey M, Singh S C. 1997. Simultaneous velocity and interface tomography of normal-incidence and wide-aperture seismic traveltime data. Geophys. J. Int., 131(1): 87-99. DOI:10.1111/gji.1997.131.issue-1
Prodehl C, Kennett B, Artemieva I M, et al. 2013. 100 years of seismic research on the Moho. Tectonophysics, 609: 9-44. DOI:10.1016/j.tecto.2013.05.036
Pšenčík I, Červený V. Seismic Waves in Laterally Inhomogeneous Media. Basel: Birkhäuser, 2002.
Rawlinson N, Houseman G A, Collins C D N. 2001. Inversion of seismic refraction and wide-angle reflection traveltimes for three-dimensional layered crustal structure. Geophys. J. Int., 145(2): 381-400. DOI:10.1046/j.1365-246x.2001.01383.x
środa P, Czuba W, Grad M, et al. 2006. Crustal and upper mantle structure of the Western Carpathians from CELEBRATION 2000 profiles CEL01 and CEL04:seismic models and geological implications. Geophys. J. Int., 167(2): 737-760. DOI:10.1111/gji.2006.167.issue-2
Teng J W, Zhang Z J, Zhang X K, et al. 2013. Investigation of the Moho discontinuity beneath the Chinese mainland using deep seismic sounding profiles. Tectonophysics, 609: 202-216. DOI:10.1016/j.tecto.2012.11.024
Teng J W, Deng Y F, Badal J, et al. 2014. Moho depth, seismicity and seismogenic structure in China mainland. Tectonophysics, 627: 108-121. DOI:10.1016/j.tecto.2013.11.008
Thybo H, Artemieva I M, Kennett B. 2013. Moho:100 years after Andrija Mohorovicic. Tectonophysics, 609: 1-8. DOI:10.1016/j.tecto.2013.10.004
Tian Y, Zhao D P, Liu C, et al. 2009. A review of body-wave tomography and its applications to studying the crust and mantle structure in China. Earth Science Frontiers , 16(2): 347-360.
Um J, Thurber C. 1987. A fast algorithm for two-point seismic ray tracing. Bull. Seismol. Soc. Am., 77(3): 972-986.
Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions. Geophysics, 55(5): 521-526. DOI:10.1190/1.1442863
Vigh D, Starr E W, Kapoor J. 2009. Developing earth models with full waveform inversion. The Leading Edge, 28(4): 432-435. DOI:10.1190/1.3112760
Xu T, Xu G M, Gao E G, et al. 2004. Block modeling and shooting ray tracing in complex 3-D media. Chinese J. Geophys. , 47(6): 1118-1126.
Xu T, Xu G M, Gao E G, et al. 2006. Block modeling and segmentally iterative ray tracing in complex 3D media. Geophysics, 71(3): T41-T51. DOI:10.1190/1.2192948
Xu T, Zhang Z J, Zhao A H, et al. 2008. Sub-triangle shooting ray-tracing in complex 3D VTI media. J. Seism. Explor., 17(2): 131-146.
Xu T, Zhang Z J, Gao E G, et al. 2010. Segmentally iterative ray tracing in complex 2D and 3D heterogeneous block models. Bull. Seismol. Soc. Am., 100(2): 841-850. DOI:10.1785/0120090155
Xu T, Li F, Wu Z B, et al. 2014. A successive three-point perturbation method for fast ray tracing in complex 2D and 3D geological models. Tectonophysics, 627: 72-81. DOI:10.1016/j.tecto.2014.02.012
Xu T, Zhang Z J, Tian X B, et al. 2014. Crustal structure beneath the Middle-Lower Yangtze metallogenic belt and its surrounding areas:Constraints from active source seismic experiment along the Lixin to Yixing profile in East China. Acta Petrologica Sinica , 30(4): 918-930.
Xu T, Zhang Z J, Liu B F, et al. 2015. Crustal velocity structure in the Emeishan Large Igneous Province and evidence of the Permian mantle plume activity. Science China:Earth Science, 58(7): 1133-1147. DOI:10.1007/s11430-015-5094-6
Zelt C A, Smith R B. 1992. Seismic traveltime inversion for 2-D crustal velocity structure. Geophys. J. Int., 108(1): 16-34. DOI:10.1111/gji.1992.108.issue-1
Zelt C A, Barton P J. 1998. Three-dimensional seismic refraction tomography:A comparison of two methods applied to data from the Faeroe Basin. J. Geophys. Res., 103(B4): 7187-7210. DOI:10.1029/97JB03536
Zhang Z J, Bai Z M, Klemperer S L, et al. 2013. Crustal structure across northeastern Tibet from wide-angle seismic profiling:constraints on the Caledonian Qilian orogeny and its reactivation. Tectonophysics, 606: 140-159. DOI:10.1016/j.tecto.2013.02.040
Zhang Z J, Teng J W, Romanelli F, et al. 2014. Geophysical constraints on the link between cratonization and orogeny:Evidence from the Tibetan Plateau and the North China Craton. Earth-Science Reviews, 130: 1-48. DOI:10.1016/j.earscirev.2013.12.005
Zhao D P, Hasegawa A, Horiuchi S. 1992. Tomographic imaging of P and S wave velocity structure beneath northeastern Japan. J. Geophys. Res., 97(B13): 19909-19928. DOI:10.1029/92JB00603
Zhao D P, Hasegawa A, Kanamori H. 1994. Deep structure of Japan subduction zone as derived from local, regional, and teleseismic events. J. Geophys. Res., 99(B11): 22313-22329. DOI:10.1029/94JB01149
Zhao F F, Ma T, Xu T. 2014a. A review of the travel-time calculation methods of seismic first break. Progress in Geophysics , 29(3): 1102-1113. DOI:10.6038/pg20140313
Zhao F F, Zhang M H, Xu T. 2014b. A review of body wave traveltime tomography methods. Progress in Geophysics , 29(3): 1090-1101. DOI:10.6038/pg20140312
Zhou L Q, Liu F T, Chen X F. 2006. Simultaneous tomography of 3-D velocity structure and interface. Chinese J. Geophys. , 49(4): 1062-1067.
白超英, 黄国娇, 李忠生. 2011. 三维复杂层状介质中多震相走时联合反演成像. 地球物理学报, 54(1): 182–192. DOI:10.3969/j.issn.0001-5733.2011.01.019
黄国娇, 白超英. 2013. 多震相走时联合三参数同时反演成像. 地球物理学报, 56(12): 4215–4225. DOI:10.6038/cjg20131224
李飞, 徐涛, 武振波, 等. 2013. 三维非均匀地质模型中的逐段迭代射线追踪. 地球物理学报, 38(6): 823–832. DOI:10.6038/cjg20131026
田有, 赵大鹏, 刘财, 等. 2009. 体波走时层析成像方法及其在中国壳幔结构研究中的应用评述. 地学前缘, 16(2): 347–360.
徐涛, 徐果明, 高尔根, 等. 2004. 三维复杂介质的块状建模和试射射线追踪. 地球物理学报, 47(6): 1118–1126.
徐涛, 张忠杰, 田小波, 等. 2014. 长江中下游成矿带及邻区地壳速度结构:来自利辛-宜兴宽角地震资料的约束. 岩石学报, 30(4): 918–930.
赵烽帆, 马婷, 徐涛. 2014a. 地震波初至走时的计算方法综述. 地球物理学进展, 29(3): 1102–1113. DOI:10.6038/pg20140313
赵烽帆, 张明辉, 徐涛. 2014b. 地震体波走时层析成像方法研究综述. 地球物理学进展, 29(3): 1090–1101. DOI:10.6038/pg20140312
周龙泉, 刘福田, 陈晓非. 2006. 三维介质中速度结构和界面的联合成像. 地球物理学报, 49(4): 1062–1067.