引言

青藏高原块体的地壳及地幔结构已成为目前国际地学界的研究热点(McNamara et al，1997; 曾融生等，1998),如中美合作的INDEPTH计划,其宗旨在于除了探测青藏高原块体的深部构造外,还在于寻求大陆地球动力学的力源及形成机制(Zhao,1993). 进而使与之相邻的青藏高原北缘(雷建设等,2002)、 东北缘(陈九辉等，2005)及东部地区(Wang et al,2007)也成为不少学者的研究课题,学者们试图解释构造应力场的继承性和传导性,得到了较好的结果. 华北(郭飚等,2004； 嘉世旭,张先康,2005； 王志铄等,2008； 吕作勇,吴建平,2010)和华南块体(江为为，2007)的地壳乃至上地幔结构已有较多且更为详细的研究报道. 作为上述地质块体过渡区的黄土高原及邻区的地壳乃至上地幔结构目前未见详细的研究报道,仅有关于渭河断陷盆地的部分结果(丁韫玉等，2000). 深入研究这一地区的深部构造将为我们系统理解大陆地球动力学的形成机理，汾渭盆地地裂缝形成的区域构造背景，以及陕北油气田的可持续开发将有着重要的科学意义.

如上所述,该区已有的地震层析成像研究主要是渭河断陷盆地 (33.5°—35.8°N,107°—110°E)的地壳速度结构,其结果并未穿透莫霍面的深度 (该区莫氏面的深度一般为40 km). 为了研究黄土高原大区域的地壳结构,本文选取了相对较大的区域,即31.0°—39.0°N,103°—113°E的矩形区域,基本上覆盖了黄土高原和汾渭盆地的大部分地区.

采用地震走时联合反演方法来重建黄土高原及邻区的三维地壳速度结构,其成像深度可达到50 km. 在地震走时联合反演中弯曲地震射线的追踪采用了改进型最短路径方法(Bai et al,2007),而非线性反演算法则采用共轭梯度求解带约束的阻尼最小二乘法问题(Bai,Greenhalgh,2005b),同时对速度模型和震中位置进行了更新.

1.   反演方法简述

地震走时联合反演进行速度模型和震源位置的同时更新是目前常用的反演方法,流行的方法较多,例如双差地震成像法(double-difference seismic tomography method,简写为DDT) (Zhang,Thurber,2003). 其本质依然是寻求一种实际意义上的非线性反演算法. 要解决这一问题,需要作以下两方面的深入研究: ① 正演算法应能够较好地处理弯曲地震射线的追踪或是地震射线在速度不连续界面的反射、 折射及透射等; ② 反演算法应具有非线性反演的特点,即能够反映弯曲射线路径上地震走时关于速度和震源位置的一阶或二阶导数变化量的信息. 鉴于上述考虑,我们提出了一种非线性的地震走时联合反演算法,并进行了数值模拟实验,最终取得了较好的成像效果. 非线性的地震走时联合反演算法主要包括3方面的内容: 弯曲地震射线追踪,多地震联合精确定位和非线性反演.

1.1   弯曲地震射线的追踪

弯曲地震射线的追踪与传统的最短路径方法不同,我们采用基于单元模型的方法，即改进型最短路径算法(modified shortest path method，简写为MSPM) (Bai et al,2007). 其优点在于: ① 模型单元化分时引入了主节点和次级节点的概念; ② 定义了三次线性拉格朗日插值函数; ③ 计算速度与传统的最短路径算法相比至少可提高一个数量级； ④ 易于跟踪计算后续震相(Bai et al,2009, 2010; 唐小平,白超英,2009a, b). 这样就可以充分利用MSPM算法计算速度快、 稳定性好、 全局解的优势，同时又避免了传统算法(shortest path method，简写为SPM)不便于跟踪计算后续震相的不足.

1.2   拟非线性的反演算法

非线性走时成像问题可由局部线性化后的循环解近似给出,在反演过程中带先验信息的约束或物理条件的约束可使反演解更逼进真实解. (Zhou等1992)系统地对比了几种流行的反演算法(例如，singular value decomposition； algebraic reconstruction technique； least squares QR-factorization； conjugate gradient),提出了一种用共轭梯度法求解带约束阻尼最小二乘问题的算法(damped minimum norm,least squares and constrained problem with CG approach,简写为DMNCLS-CG). 通过适当简化,DMNCLS-CG可退化成Scales算法(Scales,1987). 随后我们将带有一阶导数变化信息的Jacob矩阵用于反演算法中去求解非线性地震走时反演问题 (Bai,2004; Bai,Greenhalgh,2005a， b). 结果表明，由于受先验信息约束或物理约束解的非唯一性得到较大地改善，该反演算法对走时数据中的随机误差并不敏感. 实践证明它十分适用于稀疏矩阵的反演问题,且对计算机内存要求不高.

1.3   多地震同时定位

利用MSPM算法的优势,我们研发了一种多地震同时定位算法(multiple event location，简写为MEL),它十分适用于地震台网数据的处理(Bai,Greenhalgh,2006). 其优点在于: ① 地震定位中射线追踪(正演过程)只需进行一次; ② 可同时进行多次地震的联合定位而不需要增加额外的计算用时; ③ 由于引入了辅助震源,具有偏导意义的Jacob矩阵在正演过程中即可求出; ④ 解的稳定性好,且对初始值的依赖程度不高; ⑤ 速度快,其定位时间主要取决于一次正演所用的时间.

1.4   地震走时的联合反演成像

带约束的阻尼最小二乘法用于求解地震走时联合反演问题时可写成

上述模型空间的未知参数向量 m 既包括未知速度参数又包括未知震源位置参数. Jacob导数矩阵也相应包含两部分的贡献: 即走时关于速度和震源位置的变化率. 因而从震源到某一地震台站的地震走时误差可由式(2)表示:

式中，t_ij表示第j个台站接收到的第i个震源的地震走时,x_ik是第i个震源在3个方向上的位置坐标(k=1,2,3),而v_k则是长度为N的模型空间中第k个未知速度参数. 因此,Jacob导数矩阵中走时关于速度的变化率可由(Bai和 2005b)文章中公式(10)求出,而走时关于震源位置的变化率可仿照(Bai和Greenhalgh 2006)文章中公式(1)解得. 这样就得到了参加反演所用的Jacob导数矩阵. 结合MSPM，DMNCLS-CG，MEL，以及带有一阶偏导信息的Jacob分析矩阵,就可以形成一套用于地震非线性走时联合反演成像的算法,同时更新速度模型和震源位置.

2.   研究区地质构造背景

黄土高原主要覆盖在鄂尔多斯块体及其周缘断陷盆地上,并与4个地质块体相连． 其东北向为华北地块,东南向为华南地块，西南向为青藏高原板块，西北向为阿拉善块体(图1). 青藏高原受到来自印度洋板块北偏东向的俯冲,华北及华南地块则主要受到太平洋板块向西的推挤. 因此,黄土高原主要间接受到上述两个力源的联合作用,故可称之为应力调整区. 由于其力源的多样性和局部调整性,使之地壳结构更具有局部区域性(横向和纵向变化差异显著) 的特点. 黄土高原的地震构造特点是存在类似于U形的地震带,其右翼为汾渭断裂带,左翼为南北地震带的北端,并与祁连山断裂带东端交汇. 而U型断裂带内为汾渭断陷盆地,向北与鄂尔多斯块体及燕山构造带相遇. 其南部则为秦岭造山带,再往南则与川北地块相连.

图  1  研究区基本构造示意图

实线表示活动断层; 虚线表示不活动断层

Figure  1.  Basic tectonic structure of the study region

Solid lines show active faults,dash lines denote inactive faults

下载: 全尺寸图片 幻灯片

3.   地震资料预处理

选取上述研究区内1984—1993年发生的242次M_S≥3.0地震的54个地震单台报告,共收集到可用的5 000多条P波到时资料. 这54个地震台站主要分布在陕西、 山西、 甘肃及宁夏等地区,同时也采用了河南和湖北个别地震台站的资料(图2). 这些地震主要是浅源地震,震源深度Z≤30 km. 在走时资料的选取中采用了相关台站的互相关检验,剔除了误差较大的走时资料. 射线的覆盖率在研究区的中北端和中南端较差,特别检测板(checkerboard)分辨率实验(Inoue et al,1990)表明,在陕西延安与榆林之间有一个分辨率似椭圆形的空白区 (图3). 几种不同单元尺度的检测板分辨率实验表明,采用50 km×50 km×50 km单元进行模型参数化其分辨率最佳,这一单元尺度也较为接近地震台网的平均间距.

图  2  地震(实心圆)、 台站(十字)及随后展示剖面(虚线)分布图

Figure  2.  Distribution of earthquakes (solid circles) and stations (crosses).Dash lines show position of vertical profiles analyzed in the following

下载: 全尺寸图片 幻灯片

图  3  检测板试验. 十字表示台站; 虚线矩形表示最佳分辨区; 实线椭圆表示分辨空区

Figure  3.  Checkerboard test result. Crosses stand for stations; dash square represents best resolution area; ellipse stands for a resolution gap

下载: 全尺寸图片 幻灯片

4.   黄土高原三维P波速度的重建

由P波走时时距曲线采用Wiechert-Herglotz一维反演方法(Bullen，1963),可得到一维近似层状P波速度模型. 该模型可作为反演的初始速度模型(图4). 在最终的反演中,经过20次迭代走时均方根残差由开始的1.637 s收敛于0.402 s,相对下降超过75%. 此走时均方根残差终止值也与地震单台到时拾取误差相应. 反演前后震源位置的平均改正量为2.74 km,其结果如下.

图  4  原始数据时距曲线(a)和用Wiechert-Herglotz方法得到的一维速度模型(b)

Figure  4.  Raw time-distance curve for P-wave arrivals (a) and the deduced 1-D velocity model from relocated dataset,which was used as the initial model in the following 3-D inversion (b)

下载: 全尺寸图片 幻灯片

4.1   随深度变化的横向P波速度场

研究区包括不同的地质构造单元,因此讨论其局部速度乃至速度结构随深度的变化是十分必要的,它为我们了解不同构造单元的属性以及划分其构造边界提供了科学依据.

4.1.1   浅部地壳P波速度特征

此深度范围内的显著特征是存在着与浅部地震带相一致的U字型高速异常区,图5在Z= 5,10 km的平面上大体相似. 其中在山西永济和陕西彬县分别有一低速异常体插入. 秦岭及南部山区与北部地区则主要表现为P波低速异常. 限于篇幅,图5仅给出了5 km平面上的P波波速相对于背景值的百分比变化. 一般而言,研究区的东北为高速异常区,而西南则主要表现为低速异常,这种高、 低速异常幅值的绝对变化可达15%. 上述浅部P波速度分布特征可延续至Z=15 km,但其异常形态和幅度均有所变化.

图  5  P波速度在Z=5 km水平面上相对于背景值的百分比变化空心圆表示台站所在地; 实线表示速度等值线; Δv/v₀表示速度扰动

Figure  5.  The percentage P-wave variation against the background velocity at the depth of 5 kmOpen circles represent stations; dash lines stand for equal velocity contours; Δv/v₀ represents velocity perturbation

下载: 全尺寸图片 幻灯片

4.1.2   中部地壳P波速度特征

此深度范围内其醒目的速度异常结构是上述U字型高速异常区发展成Y字型高速异常带,其下部已延伸至陕西安康(图6). 这一高速异常特征在Z=25 km平面上尤为显著. 在此地壳深度范围内秦岭及南部山区和北部地区依然表现为P波低速异常. 地震大多发生在高速异常区的边缘或高、 低速异常区的交汇处. 在研究区西南部出现了走向北西后转至近东西的高、 低速相间的异常区带，而在研究区东南部则出现了走向北东的高速异常带. 另外，宁夏盐池的低速异常区有所扩展.

图  6  P波速度在Z=20 km平面上相对于背景值的百分比变化空心圆表示台站所在地; 实线表示速度等值线; Δv/v₀表示速度扰动

Figure  6.  The percentage P-wave variation against the background velocity at depth of 20 kmOpen circles represent stations; dash lines stand for equal velocity contousr; Δv/v₀ represents velocity perturbation

下载: 全尺寸图片 幻灯片

4.1.3   下部地壳P波速度特征

该区速度不连续界面(莫氏面)的平均深度为Z=30—40 km,所以有必要讨论接近莫霍界面深度上的P波速度分布. 在此深度上速度分布特征则与其上截然 不同,主要表现为东部大面积高速异常区和西部块状连接的高速异常区带,其间由S型的低速异常带自南向北弯延穿过(图7),同时山西永济下的楔形低速异常插入体随即消失. 而Z=35 km平面又与Z=30 km平面上P波速度分布存在较大的差异,说明该区莫霍速度不连续界面的局部特征较为显著.

图  7  P波速度在Z=30 km平面上相对于背景值的百分比变化空心圆表示台站所在地; 实线表示速度等值线; Δv/v₀表示速度扰动

Figure  7.  The percentage P-wave variation against the background velocity at depth of 30 kmOpen circles represent stations; dash lines stand for equal velocity contours; Δv/v₀ represents velocity perturbation

下载: 全尺寸图片 幻灯片

4.1.4   上地幔顶部P波速度特征

此深度范围内P波速度的分布特征除了横向局部变化外,纵向变化也较为明显. 基本上各层(Z=40,45,50 km) 内的速度变化特征继承性差,即纵向变化较大. 这一特点与地壳内速度场分布截面不同. 这可能与上地幔物质的流变性强有关,因此更具有局域性的特点. 图8仅给出了Z=45 km平面上的速度相对于背景场的变化,从中可明显看出整体上存在着近北西—南东走向的高、 低速相间的异常带分布的基本特征. 此特性可能与上地幔塑性介质在围压下的剪切应力条带分布有关(白超英,秦保燕,1990).

图  8  P波速度在Z=45 km平面相对于背景值的百分比变化空心圆表示台站所在地; 实线表示速度等值线; Δv/v₀: 速度扰动

Figure  8.  The percentage P-wave variation against the background velocity at depth of 45 kmOpen circles represent stations; dash lines stand for equal velocity contours; Δv/v₀ represents velocity perturbation

下载: 全尺寸图片 幻灯片

4.2   垂直剖面P波速度分布特征

为了详细讨论黄土高原及汾渭盆地的特征,并兼顾资料的覆盖率,我们选取了如下4条垂直剖面: 一条是沿汾渭断裂带走向,一条垂直于汾渭断裂带走向,另外两条分别沿经向和纬向的垂直剖面(图2).

4.2.1   沿汾渭断裂带走向垂直剖面的速度特征

该垂直剖面过陕西汉中、 周至、 西安、 山西蒲县、 韩城、 临汾,全长近850 km(图9). 其速度分布特征主要是浅部地壳陕西汉中—周至的V字型低速异常区,并且此异常区可延伸至地壳底部. 而在陕西西安和山西临汾下各存在一个反U字型高速异常隆起区,这些高速异常隆起区可追溯至深部 (上地幔顶部). 这种类型的速度异常结构有可能是水平向挤压和局部地幔上涌力的联合作用. 另外,从整体上来说,地壳结构存在由西南向东北逐步变浅的趋势,但不十分明显. 这一特征也与该区大范围内的地壳深度相吻合.

图  9  沿汾渭断裂带走向垂直剖面P波速度分布(图中附有速度等值线)

Figure  9.  P-wave velocity vertical profile along the Fenwei belt. Lines represent equal velocity contours

下载: 全尺寸图片 幻灯片

4.2.2   垂直于汾渭断裂带剖面的速度特征

该剖面过甘肃景泰、 静宁、 靖远、 峡门、 陕西陇县、 宁陕、 安康、 胡北竹溪，全长近920 km (图10). 此垂直剖面上的速度分布相对而言较为复杂,其地壳浅部主要分布特征是位于陕西陇县下的U字型低速凹陷异常(与陕西汉中—周至的V字型低速异常区相应,见图9); 陕西宁陕下的上高速、 下低速的异常体(与西安下的情形类似,见图9). 此外,中深地壳及上地幔顶部的速度结构主要表现为成对的局部高速隆起与低速凹陷异常.

图  10  垂直汾渭断裂带走向垂直剖面P波速度分布(图中附有速度等值线)

Figure  10.  P-wave velocity vertical profile perpendicular to the Fenwei belt. Lines represent equal velocity contours

下载: 全尺寸图片 幻灯片

4.2.3   沿纬向垂直剖面的速度特征

研究区内有许多近纬向的构造带,因此有必要讨论其沿纬向和经向速度随深度的分布特征. 考虑到地震射线的覆盖率,这里给出沿纬向 33.9°N的垂直剖面.

该纬向垂直剖面过甘肃成县、 陕西周至、 西安、 商县和河南卢氏(图11). 浅部主要特征是位于陕西周至、 商县和河南卢氏下分别有一个Y字型低速异常区,以及西安下的高速异常隆起. 而这些Y字型低速异常区可达到15 km深度,这是十分有趣的发现. 是否预示着活动断层的延伸深度，由于现有的射线分辨率无法识别这些Y字型低速异常区内的精细结构,所以还不能下确定性的结论,但至少说明我们所用的走时联合成像方法是有能力识别突变或高反差的速度结构,这一点在以前的数值模拟实验中已得到验证(Bai，Greenhalgh，2005b).

图  11  沿纬向(纬度为33.9°)垂直剖面P波速度分布(图中附有速度等值线)

Figure  11.  Vertical profile of P-wave velocity map along the latitude of 33.9°. Lines represent equal velocity contours

下载: 全尺寸图片 幻灯片

4.2.4   沿经向垂直剖面的速度特征

经向垂直剖面(111.0°E)过湖北郧县、 河南卢氏、 山西永济、 临汾和离石(图12). 浅部速度的主要特征是位于山西永济下的V字型低速凹陷异常,以及位于离石下面的反V字型高速异常隆起. 这种分布特征一直可延伸至深部. 这是典型的深浅速度结构贯通的实例,预示着较大区域内构造应力场的作用.

图  12  沿经向(111.0°E)垂直剖面P波速度分布(图中附有速度等值线)

Figure  12.  Vertical profile of P-wave velocity map along the longitude of 111.0°. Lines represent equal velocity contours

下载: 全尺寸图片 幻灯片

5.   讨论与结论

汾渭盆地及所处的黄土高原及邻区的地壳P波速度结构要比以往想象中的复杂. 作为中国大陆构造应力场调整区,其地壳乃至上地幔顶部速度结构理应具有局域性的特点,这一结论也在我们的研究中得到了验证. 本区的构造应力主要是来自青藏块体向北东的挤压,华北、 华南块体向西南、 向西北的推挤,以及局部地幔隆起(上涌力)的联合作用. 深入研究汾渭盆地及所处的黄土高原的深部构造特征,不仅为系统的理解和构建大陆地球动力学提供科学依据,而且为汾渭盆地地裂缝的构造成因提供科学证据,同时也为陕北油气田的持续开采提供大区域的构造背景.

三维地震走时联合成像结果表明，该研究区内浅部乃至中部的速度分布与浅部地震活动构造带有关,主要表现为U字型高速异常带; 而秦岭及南部山区和北部地区则表现为低速异常区(图5、 图6). 总体上表现为东北区高速,西南区低速的特点. 地震多发生在高速异常区的边缘或高、 低速异常区的交汇处. 莫霍面深度附近速度分布更具有局域性的特点,即差异性显著(图7). 多个垂直剖面速度分布结果说明，速度异常分布图像主要受制于水平挤压和地幔局部隆起的单一或联合作用,从而导致诸如: 高、 低速异常区相伴而生; 上高速、 下低速的异常体(图9、 图10).

由于受地震射线覆盖率(地震台网间距，地震实际分布，以及单一地震震相)的限制,我们无法进一步分辨速度异常体内的精细结构,所以也就无法确定诸如图11中的Y字型低速异常带的细微结构. 然而从以上的研究说明,采用实际意义上的非线性三维走时联合成像技术是有能力反演局部突变或者较大速度反差异常体的结构.