首都圈地区位于华北北部, 主要包括华北盆地、燕山和太行山隆起区三个构造单元(图 1), 区内地质构造复杂, 新生代构造活动强烈, 历史上曾发生过多次破坏性地震, 如1679年三河-平谷8.0级地震、1976年唐山7.8级和1998年张北6.2级地震, 给国家经济和人民生命财产造成了巨大损失.30多年来, 我国一直将该区列为地震重点防御区, 并陆续建立了较密集的数字地震台网.

图 1

Fig. 1

图 1 首都圈地区构造背景和107个地震台站分布图 图中灰度背景表示地形; 黑色实心正方形代表甚宽带台、黑色实心三角形代表宽频带台、白色黑边三角形代表短周期台; 白色曲线代表海岸线和省界, 灰色曲线为主要活动断层: F1.太行山山前断裂, F2.蔚县-延庆断裂, F3.紫荆关断裂, F4.夏垫断裂, F5.沧东断裂, F6.宁河断裂, F7.唐山断裂, F8.滦县西断裂; 黑色直线示出了图 7和图 8中垂直剖面的位置. Fig. 1 Geological background of Chinese capital region and distribution of 107 seismic stations The gray background shows topography. Black solid squares denote very-broadband stations, black solid triangles denote broadband stations, and white triangles with black sidelines denote short-period stations. White curved lines denote coastlines and boundaries between provinces, gray curved lines show major active faults: F1, Taihangshan Front Fault; F2, Yuxian-Yanqing Fault; F3, Zijinguan Fault; F4, Xiadian Fault; F5, Cangdong Fault; F6, Ninghe Fault; F7, Tangshan Fault; F8, West Luanxian Fautt. Black lines show the locations of the vertical cross sections in Fig. 7 and Fig. 8.

鉴于首都圈地区重要的地理位置、复杂的地质构造及较强的地震活动, 我国学者围绕该区深部结构开展了大量的地球物理观测和研究, 大多数研究主要利用人工地震测深和区域天然地震资料反演地壳结构^[1~5].在上地幔结构的研究上, 早期有些学者利用当时的华北电信传输台网或模拟台网记录的区域震及远震资料, 对该区的地壳、上地幔速度结构开展了一些研究工作^[6~9], 并取得了一些有意义的结果.近期, 随着地学界对华北克拉通破坏这一重大科学问题的关注, 使得该区三维速度结构的研究再次成为热点, Huang等^[10]利用新建首都圈数字地震台网记录的区域震和远震, 采用联合反演的方法得到该区地壳、上地幔三维速度结构, 结果显示首都圈地区三个构造单元的速度结构存在明显的差异, 在华北盆地、渤海湾下50~150 km的范围内存在大尺度强低速异常区, 而该低速异常区的下面存在明显的高速块体, 这种高、低速异常在空间位置上的倒置关系说明华北克拉通东部岩石圈可能发生了拆沉.最近, Dng等^[11]利用华北密集台阵及首都圈数字地震台网资料反演的该区地壳、上地幔结构也在华北盆地深部发现了高速块体下沉的残骸.这些地震层析成像的新成果使人们对首都圈地区的地壳、上地幔结构及深部动力学过程有了更深人的认识, 但到目前为止应用地震层析成像对首都圈深部结构的研究基本都是基于传统的射线理论.

有限频率层析成像是近年由Dahlen等提出的用于反演地球深部结构的新方法^{[12, 13]}, 该方法考虑了地震波在非均匀介质中的绕射作用所产生的波前恢复效应, 以及不同频率的散射波相互干涉对走时的影响, 建立了走时变化与速度异常之间更为真实的对应关系; 而此法在实际应用中的主要特点是对宽频带波形资料的分频带处理, 从而可以充分地利用宽频带资料, 提高数据的利用率^{[14, 15]}. Montell等^{[16, 17]}利用该方法反演长周期的P、PP波资料得到了有限频率全球模型, 发现了一些新的地幔柱, 并给出了存在深地幔柱的证据.Hung等^[14]应用该方法反演ICEMELT和HOTSPOT流动台阵记录的宽频带资料获得了冰岛地幔柱的清晰图像, 并通过与基于射线理论的传统走时层析成像方法的结果对比指出, 有限频率层析成像方法反演得到的速度扰动的幅值是射线方法的2~3倍, 很好地改善了速度模型的分辨率.最近, Ren等^[18]利用Namche Barwa和Bhutan流动台阵以及LSA地震台记录的695个震级大于5.5级的远震事件的波形数据, 采用有限频率层析成像反演得到的研究结果显示, 西藏东南部之下的地幔岩石圈可能发生了拆沉.

本次研究应用有限频率层析成像方法反演首都圈数字地震台网记录的远震波形资料, 得到该区一个新的上地幔三维P波速度结构模型, 并基于此模型对上地幔速度结构特征及其构造意义进行探讨.

传统的走时层析成像方法是基于射线理论, 即将地震波看作无限高频的光波, 用一条狭窄的射线来代表地震波的传播路径, 地震波走时表示为慢度(速度的倒数)沿几何射线路径的线积分, 而射线路径以外的速度扰动对走时没有任何影响.

实际中, 由于一定体积的地震震源、地震波沿几何射线路径的衰减以及记录仪器的限制, 观测到的地震波都是有限频带宽度的, 而有限频带宽度的地震波对偏离几何射线路径的三维速度结构敏感^{[12, 13]}. Dahlen等^{[12, 13]}针对有限频带宽地震波的这一性质, 结合体波传播理论和线性Born单次散射理论来模拟通过地球内部任意一点的散射波到达台站的波形扰动, 利用不同频率的散射波与未经散射的直达波之间相互干涉的结果, 建立了有限频率走时变化与中心几何射线附近第一菲涅耳带内三维速度扰动之间的关系, 并由此提出了有限频率层析成像:

(1)

其中, δt是散射波与直达波的走时差; δc(x)/c(x)是速度扰动量; K(x, w)是三维Fréchet灵敏度核, 代表了散射点x处的速度异常对走时变化的贡献.从上式可以看出, K不仅是位置x的函数, 同时还是频率w的函数, 即不同频率的地震波对不同空间范围内的速度结构敏感, 上式将地震波的走时差更精确地表示为灵敏度核函数在射线周围一定空间内的体积分, 这正是有限频率层析成像与建立在射线理论下的传统走时层析成像最本质的差别.因此, 确定在一定速度模型下的有限频率灵敏度核函数是有限频率层析成像中的关键问题.利用"旁轴近似"可得到一个效率较高的计算灵敏度核的公式^[12]:

(2)

其中, ΔT=T'+T"-T, T'、T"与T分别为震源到散射点、散射点到台站, 以及震源到台站的走时; c与c_r分别为初始模型中散射点x与接收点r处的波速; R、R‘和R"分别是"震源到台站"、"震源到散射点"和"台站到散射点"射线的几何扩散因子; |S_syn(w)|²则是震源时间函数的功率谱.从(2)式可以发现, 当散射点x位于中心射线路径上时灵敏度核为零, 也就是说有限频率走时完全不受中心射线路径上速度扰动的影响, 反而对射线周围的速度扰动较为敏感, 这同射线理论下只与中心射线路径上的速度结构有关完全相反.从不同频段的灵敏度核(图 2)可以看出, 波的频率越高其灵敏度核越窄, 且强度越大, 最接近于传统的射线(图 2a)而频带范围越低, 灵敏度核越宽, 受离射线越远处速度不均匀性的影响越大(图 2(b, c)), 即采样范围越大, 但这时的灵敏度却会降低.由此可见, 不同频段的三维灵敏度核在三维模型空间内的相互交叉, 可加大数据采样的空间范围, 降低数据不均匀覆盖的影响, 对异常的估计可能更加充分.

图 2

Fig. 2

图 2 远震(震中距Δ=70°) P波走时在三个不同频段内的有限频率Fréchet核 (a)高频1.0~2.0Hz; (b)中频0.1~1.0 Hz; (c)低频0.05~0.1Hz; R为到地心的距离, 弧度表示距离.在射线平面内Féchet核表现为弯曲的椭圆形, 几何射线路径(黑色曲线)上的灵敏度为零, 最大灵敏度出现在几何射线两边的旁瓣上, 即图中的红色和蓝色区域, 它们展开的宽度正比于(其中, λ为波长, L为传播距离), 图的右侧给出了色标. Fig. 2 Finite-frequency Frechet kernels for travel-times of teleseismic (Epicentral distance Δ=70°) P phases in three different frequency bands (a) High-frequency, 1.0~2.0 Hz; (b) Intermediate-frequency, 0.1~1.0 Hz; (c) Low-frequency, 0.05~0.1 Hz; R denote geocentric distances, radians denote distances. In the ray plane cross sections, Frechet kernels exhibit bent elliptical shape with zero sensitivity right on the geometrical ray paths (black curved lines). The maximum sensitivity is in the two side-lobed, red and blue regions, and their widespread cross-path scales are proportional to (where λ is the wavelength and L is the propagation distance).

在地震走时层析成像中, 观测走时数据是有限的, 而所要求解的速度扰动量在模型空间内的变化却是连续的.通常, 可将此反演问题写成如下形式:

(3)

其中, d_i和e_i分别代表第i个走时数据及其观测误差(i=1, 2, …N); x是三维模型空间D中的位置矢量; m(x)为模型函数; 数据核g_i(x)表示第i个走时数据对模型函数的敏感程度.在有限频率层析成像中, g_i(x)是Frfechet灵敏度核在无限小立方体中的三维体积分Kd³x^[14]; 而在射线理论中, g_i(x)则与通过无限小立方体的射线长度有关.应用三维网格节点对模型空间离散化后, 可将(3)式写成如下的矩阵方程^{[14, 15, 18, 19]}:

(4)

这与基于射线理论走时层析成像得到的矩阵方程具有同样的形式, 只是其中某些元素的含义有所不同.在有限频率层析成像中d为观测走时向量, 由N个相对走时构成, 即[δt₁, δt₂, …, δt_N]^T, 它可以通过对地震波采用分频带多道互相关处理得到; m=[m₁, m₂, …, m_L为要求解的相对于初始模型的速度扰动向量, 维数L是节点数.G是系数矩阵, 其元素G_ij是第i个相对走时的灵敏度核函数在节点j处一定空间范围内的体积分, 而此相对走时的灵敏度核就是与各自走时相应的单个灵敏度核的简单差值^[14].所以应用台网中一系列地震台站记录同一震相的相对走时数据和相应的灵敏度核可以建立矩阵方程(4), 由此可对台网之下的速度结构进行约束.

本文在反演计算时, 采用阻尼最小二乘法对模型进行约束, 得到模型的阻尼最小二乘解:

(5)

其中, I为L×L的单位矩阵; θ²为阻尼因子, 一般采用模型L₂范数与走时方差的折衷分析方法来确定合适的阻尼因子^[20], 用LSQR方法^[21]进行反演.

研究中我们选取了首都圈数字地震台网(图 1) 2003年9月~2005年12月记录的300多个远震事件的波形资料, 为了保证波形清晰且有较高的信噪比, 所用事件的震级基本大于5.5级, 远震的震中及方位分布见图 3a, 可以看出具有比较完备的方位角覆盖.首都圈数字地震台网107个台站中有将近一半是宽频带或甚宽带数字地震台, 全部为三分量记录, 由于本次只研究P波速度结构, 所以仅对垂向分量的波形记录进行处理.在上述方法中我们看到, 建立有限频率层析成像矩阵方程(4)主要是通过数据处理得到相对走时向量及灵敏度核的系数矩阵.

图 3

Fig. 3

图 3 本次研究中应用的远震事件震中分布和模型空间离散化示意图 (a)图的中心在本文研究区的中心(39.5°N, 117.0°E), 同心圆上标示的数字30°、60°、90°为震中距, 黑色圆点为震中位置; (b)图中等间距的网格示出了离散化的模型空间, 黑色粗边方框示出了本文的研究区(7.5°~41.5°N, 113.5°~120.5°E), 曲线代表海岸线. Fig. 3 Distribution of teleseismic events used in this study and model space parameterization (a) Themap center is the center of the present study region (39. 5°N, 117. 0°E), concentric circles denote epicentral distances of 30°, 60° and 90°, black dots denote earthquakes.(b) Equally spaced grids show the parameterized model space, the thick black frame shows the present study region (37.5°~41.5°N, 113. 5°~120. 5°E), curved lines denote coastlines.

本次研究采用分频带多道互相关确定相对走时向量.根据首都圈数字地震台网三种不同地震仪的频带范围, 在分频带处理中分别取高频段为1.0~2.0 Hz、中频段为0.1~1.0 Hz及低频段为0.05~0.1Hz.首先, 利用巴特沃斯带通滤波器将地震波分为上述高、中、低三个频段, 对于甚宽带台和宽频带台记录的远震波形在三个频段内都需要进行处理, 而对于短周期台只在高频段进行分析, 图 4是经过分频处理后不同频段内波形的一个例子.然后, 对不同频段的波形资料进行互相关处理^[22].在进行互相关处理时需要确定时窗, 该时窗由起始拾取时间T₀、窗长和插人时长(在T₀之前的时间)三个参数共同确定.在拾取T₀时需先确定一个参考道, 其他道的T₀则按照与参考道一致的方式(即都以波峰或都以波谷为标准)由程序自动拾取, 并将拾取的结果标在分频带波形上, 最后再人工进行检查, 若发现不合理的拾取结果则通过手动调整; 在确定窗长和插人时长时需根据T₀的整体位置和具体的频段范围适当设置, 窗长的设定一般要比频段内的最长周期大些(例如对低频段可取30 s左右, 而对高频段可取5 s左右), 插入时长的选择要保证时窗能够包含被计算的地震波(本文多取T₀前半个至一个周期).时窗确定后, 用分频带多道互相关计算各地震道相对于参考道的远震P波震相的相对走时(图 4).本文最终得到18499个高质量的P波相对走时数据, 其中高频段数据6212个, 中频段数据7229个, 低频段数据5058个.

图 4

Fig. 4

图 4 不同频段的波形及多道互相关(MCCC)示意图 (a)高频1.0~2.0Hz; (b)中频0.1~1.0Hz; (c)低频0.05~0.1Hz; 短竖线标示了多道互相关的拾取时间. Fig. 4 Waveforms in different frequency bands and multi-channel cross-correlation (MCCC) (a) High-frequency, 1.0~2.0 Hz; (b) Intermediate-frequency, 0.1~1.0 Hz; (c) Low-frequency, 0.05~0.1 Hz. Short vertical lines showthe timepicks of themulti-channel cross-correlation.

在计算远震P波灵敏度核前, 先要对模型空间进行离散化处理.本次研究以(9.5°N, 117.0°E)为中心, 沿经度、纬度和深度方向划分37×35×31的三维网格, 节点间距均取约35 km.在远震有限频率层析成像中, 所取模型空间离散化的范围比研究区的范围要大很多^[14], 原因在于使用的是相距不远(与震中距相比)的两个台站间的相对走时, 并认为此相对走时只是由这两个台站下的速度结构异常引起, 而在更远的区域与此相对走时对应的两条射线的路径几乎重合, 它们的灵敏度核之差几乎为零.但是, 如果模型空间太小, 模型外的这两个灵敏度核包括的区域会有很大不同, 其差值不为零, 会导致模型外的速度异常也包括在研究区域内, 所以模型空间离散化的范围应大些(如图 3b所示), 使模型外的灵敏度核尽可能一致.图 2是与本文频带划分相应的不同频段灵敏度核的一个例子.这样, 对于每个频段由不同台站间的相对走时数据和相应的灵敏度核的差值建立矩阵方程(4), 并对此方程进行反演求取模型.

此外, 因为本文是用远震数据反演上地幔结构, 为了尽量减少地壳速度结构对反演的上地幔速度结构的影响, 我们还对相对走时数据进行了地壳校正和高程校正, 在校正中本文采用一维地壳模型CRUST2.0, 并在反演中对每个台站加入一个自由项(台站校正项), 来吸收浅层结构中无法反演出的异常所引起的走时变化.

本次研究通过射线密度分布、检测板分辨率分析以及恢复测试等多种分辨率分析方法检测了反演结果的可靠性.图 5给出了以3×3×3网格的输入模型进行检测得到的不同深度的检测板分辨率测试结果.检测板输入模型是在IASP91模型^[23]的基础上加上正负相间的扰动构成的, 输入的速度扰动在水平和垂直方向上均逐渐从-6%变化到6%(图 5a).可以看到各深度层上大部分区域基本能正确恢复正负速度的相对变化, 但异常的幅度很难恢复到输入模型给定的值, 且随着深度的增加能恢复的幅值不断降低, 到385 km深度层上一些节点的速度扰动幅值只能恢复到原来的30%以内(图 5b)特别是渤海湾及周边地区, 由于基本没有台站分布, 浅层上地幔的分辩率较差, 我们认为异常不太可靠, 因此后面对这些异常不做太多的解释.

图 5

Fig. 5

图 5 检测板分辨率测试结果 (a)输入模型, 输入的速度扰动在水平和垂直方向上均逐渐从-6%变化到6%; (b)输出模型; 层深度标在各图的右下角; 工色和蓝色分别代表低速和高速, 图的右侧给出了速度扰动的色标; 图中曲线代表海岸线. Fig. 5 Results of checkerboard resolution test (a) Input models, the velocity perturbations change from -6% to 6% in horizontal and vertical directions; (b) Output models. The depth of each layer is shown at the lower right corner of each map. Red and blue colors denote low-and high-velocities, respectively.The velocity perturbation scale is shown on the right. Curved lines denote coastlines.

本次研究采用阻尼最小二乘法对模型进行约束, 在反演中我们进行了大量的试算, 采用模型L₂范数与走时方差的折衷分析方法确定了合适的阻尼因子为10, 此时相应的走时方差减小65%.利用LSQR算法^[21]完成最终的模型求解, 得到各节点处相对于IASP91^[23]初始模型的速度扰动量.图 6是反演得到的不同深度层上的速度结构平面图, 图 7 (a~d)给出了穿过一些重要构造部位的速度结构剖面(剖面位置示于图 1中).反演揭示的P波速度变化达到6%.

图 6

Fig. 6

图 6 不同深度的P波速度扰动 层深度标在各图的右下角; 工色和蓝色分别代表低速和高速, 白色圆点代表大于6.0级的地震, 图的下方给出了速度扰动的色标和震级图例; 黑色粗线代表主要活动断层, 黑色细线代表海岸线; 图中BJ表示北京、TJ表示天津、TS表示唐山、ZJK表示张家口. Fig. 6 P-wave velocity perturbation images at each depth slice The depth of each layer is shown at the lower right corner of each map. Red and blue colors denote low-and high-velocities, respectively. White dots denote earthquakes with M≥6.0. The velocity perturbation and earthquake magnitude scales are shown at the bottom. Thick black curved lines show major active faults, thin black curved lines show coastlines. BJ, Beijing; TJ, Tianjin; TS, Tangshan; ZJK, Zhangjiakou.

图 7

Fig. 7

图 7 P波速度扰动剖面图 (a)~(d)有限频率层析成像模型, 弧度表示距离; (e)~(h)射线层析成像模型[10]; 红色和蓝色分别代表低速和高速; 地形图标注在剖面上方, H表示高程; 图的右侧给出了速度扰动的色标; 黑线代表410 km间断面; 图中NCB表示华北盆地, YS表示燕山隆起区, THS表示太行山隆起区, BHB表示渤海湾, F1表示张家口一蓬莱断裂带, F2表示太行山山前断裂带; 剖面位置示于图 1中. Fig. 7 Vertical cross sections of P-wave velocity perturbations (a)~(d) Finite-frequency tomographic model, radians denote distances.(e)~(h) Ray tomographic model [10]. Red and blue colors denote low-and high-velocities, respectively. The surface topography along each profile is shown above each cross section. H denote altitudes. The velocity perturbation scale is shown on the right. Black lines denote the 410 km discontinuities. NCB, North China basin; YS, Yanshan uplift; THS, Taihangshan uplift; BHB, Bohai bay; F1, Zhangjiakou-Penglai fault zone; F2, Taihangshan front fauit zone. The locations of the vertical cross sections are shown in Fig. 1.

结果显示, 在70 km深度的速度图上(图 6a), 渤海湾周边、华北盆地、燕山隆起区北部及太行山隆起区的大部表现为低速; 而燕山隆起区的大部、张家口一带和胶州湾等地呈现高速.在此深度层上, 张家口一蓬莱断裂带和太行山山前断裂带两边的速度结构差异较为明显, 反映出这两条区内最主要的断裂带可能已切入岩石圈.105 km和140 km深度层上的速度结构特征大体相似(图 6(b, c)), 以张家口一蓬莱断裂带为界, 东北部的燕山隆起区整体为高速, 西南部的太行山隆起区、华北盆地大部及渤海湾内则主要表现为低速异常, 但低速异常的内部还存有小范围的不均匀高速岩体.张家口-蓬莱断裂带和太行山山前断裂带作为速度结构边界在这两个深度上仍清楚可见.175 km到280 km深度层上(图 6(d~f)), 华北盆地及渤海湾高速区的范围逐渐扩大, 而燕山隆起区内也出现了小范围的低速区.在175 km以下的深度层上(图 6(d~i))张家口-蓬莱断裂带和太行山山前断裂带作为深部构造边界已变得模糊不清, 整个区域的速度不均匀性也逐渐减弱.

从速度结构剖面图可见(图 7(a~d)), 以张家口一蓬莱断裂带(图中F1)和太行山山前断裂带(图中F2)为构造边界, 区内三个构造单元在150 km以上具有明显不同的速度结构特征:燕山隆起区下表现为主体的高速, 且高速异常可从地表一直延伸到200多公里深(图 7(b~d); 太行山隆起区内的速度不均匀性较为明显, 在主体的低速结构内分布有小范围的高速块体, 总的说来该区北部速度偏低(图 7d)南部速度略高(图 7a)在华北盆地、渤海湾下存在较大范围的强低速异常, 华北盆地下的低速异常从50 km延伸到约150 km深(图 7(a, b)), 而渤海湾内的低速异常则可追踪到200 km以下(图 7c).此外, 在华北盆地的低速异常区下还隐约可见高速异常的分布(图 7(a, b).

上地幔软流圈具有较低的地震波速度、高电导率、高温等特征, 通过这些特征可以推测软流圈顶部的埋深, 进而解释岩石圈厚度的变化.在华北东部依据大地电磁测深并结合大地热流值数据推测的上地幔高导层的顶部埋深为60~80 km, 呈NNE向展布^[24], 而本文揭示的华北盆地和渤海湾下大范围的强低速异常的顶面深度为50~70 km (图 7(a~c)), 与上述推测高导层的埋深较为接近, 且展布方向相似.由此可见, 华北盆地和渤海湾下软流圈的顶面埋深可能小于70 km, 而在许多大陆地区软流圈顶面的深度也就是岩石圈的平均厚度应为110 km左右, 这说明华北盆地和渤海湾下的岩石圈厚度比平均厚度减薄40 km以上.

魏文博等^[25]大地电磁测深的研究结果也显示华北地区东部岩石圈厚度约为50~80 km, 其减薄的区域和厚度与本文的结果较为接近.最近, 走时地震层析成像的结果也揭示华北盆地、渤海湾岩石圈的减薄, 并结合中国大陆地幔层析成像结果^[26]推测岩石圈减薄可能与太平洋板块在东亚大陆下俯冲滞留所引起的一系列深部过程相关^[10].地壳泊松比的分布特征也反映了华北盆地下部上地幔软流圈物质的侵人^[27], 而区内明显的Pn波低速异常也说明下部软流圈可能与地壳部分接触^{[28, 29]}.

张家口一蓬莱断裂带在地表构造上是燕山隆起区与太行山隆起区、华北盆地区的分界线, 在上地幔浅层速度结构上是高、低速异常的变化带, 也是岩石圈减薄的边界带.可以看到区内大部分的历史强震都发生在这条构造边界带上(图 6(a, b)), 因此我们推测该带上强震的发生还可能与上地幔速度结构不均匀性有关.已有研究表明, 在部分熔融的高温软流圈的隆起部位其介质与周围介质存在温度差异, 这种温度差异可以产生附加热应力, 在隆起部位侧面顶部的张应力和剪应力都最大, 是最容易破坏和发生地震的地区^[30].按此观点, 张家口一蓬莱断裂带作为区内岩石圈速度结构的变化带, 其两侧不均匀的软流圈差异活动所引起的热应力有可能是导致大地震频发的原因之一.

本次研究是首次应用远震有限频率层析成像反演首都圈地区上地幔结构, 我们自然想到与所用资料相近的射线层析成像方法的结果进行比较^[10].总体说来, 两种方法揭示的速度结构主体是类似的, 但也存在可观的差别(图 7).两种方法都揭示了华北盆地下浅层上地幔强低速体的存在, 且异常体的位置和展布情况大体相似, 但异常大小和幅值有些差别(图 7(a, b, e, f)); 在渤海湾下, 有限频率方法显示的低速异常体延伸更深(图 7(c, g)); 在燕山隆起区下, 射线方法展示的高速异常体更完整, 有限频率方法展示的高速异常范围更大, 追踪更深(图 7(b~d), (f~h); 太行山隆起区北部都显示主体的低速异常(图 7(d, h)), 而在南部, 有限频率方法显示的高速异常范围更大(图 7(a, e).射线层析成像的结果清楚地显示了华北盆地、渤海湾强低速岩体下高速异常的存在(图 7(e~g)), Huang等^[10]将华北盆地、渤海湾下的上地幔这种强低速和强高速异常在空间上的倒置关系解释成具有重要动力学意义的华北岩石圈拆沉, 为华北岩石圈减薄的拆沉机制提供了地震学深部证据, 但在本文有限频率层析成像方法的结果中, 相应位置上的高速异常只是隐约可见(图 7(a~c)), 与基于射线层析成像的结果相比其空间范围和幅值都小一些.

为了进一步检验有限频率层析成像结果中华北盆地和渤海湾下低、高速构造的可靠性, 本文对其进行恢复测试(图 8).在不改变模型空间离散化方式的条件下, 我们依据由实际资料反演获得的三维P波速度模型, 在各节点上输人正负速度扰动构建了与这两个低、高速构造形状相似的构造体作为输人模型, 输人的速度扰动分别为正负6%(图 8(a~c)), 然后采用有限频率方法对合成数据进行反演.结果显示:低、高速体的几何形状恢复情况良好, 但速度扰动的幅值没有得到充分恢复, 特别是高速体的幅值在很多区域只恢复到输人值的40%左右(图 8(d~f)), 渤海湾下的高速体更是几乎没有得到恢复(图 8f), 而华北盆地南部的低速区恢复也欠佳(图 8e), 这可能是华北盆地南部台站较少、数据采样密度较低所致.总体上讲, 恢复测试的结果表明我们所得成像结果中的主要构造特征是可靠的.

图 8

Fig. 8

图 8 华北盆地和渤海湾下低、高速异常恢复测试结果 (a)~(c)输入模型; (d)~(f)输出结果; 弧度表示距离; 红色和蓝色分别代表低速和高速; 图的右侧给出了速度扰动的色标图中NCB表示华北盆地, YS表示燕山隆起区, THS表示太行山隆起区, BHB表示渤海湾; 剖面位置示于图 1中. Fig. 8 Results of restoring tests for low and high velocity anomalies beneath the North China basin and Bohai bay (a)~(c) Inputmodels; (d)~(f) Output results. Radians denote distances. Red and blue colors denote low-and high-velocities, respectively. The velocity perturbation scale is shown on the right. NCB, North China basin; YS, Yanshan uplift; THS, Taihangshan uplift; BHB, Bohai bay. The locations of the vertical cross sections are shown in Fig. 1.

通过结果对比可以看到, 有限频率层析成像方法的反演结果受数据覆盖的影响较大, 在数据采样和分辨率都较好的区域, 有限频率方法与射线方法的结果具有大体一致的趋势.它的主要特点是通过对宽频带资料采用分频带处理, 利用了不同频带的信息可以提高反演结果的分辨率, 所以只有在应用宽频带资料的情况下才能更充分体现其理论上的优势.在本次研究的首都圈地区, 渤海湾内没有布设地震台站, 而华北盆地大多数又是短周期台(图 1).在对波形资料进行分频带互相关处理计算相对走时的过程中, 短周期台只有高频段才有资料, 但高频段波形资料的复杂性又使得在互相关处理中很多资料不符合要求被舍弃, 这样可用的高频资料相对较少.经统计本文仅有2016个高频段相对走时数据取自于华北盆地内的短周期台, 其数目只占全部数据量(18499个)的11%, 在华北盆地区所用资料比基于射线理论的走时层析成像所用资料还要少, 这可能是在有重要构造意义的华北盆地、渤海湾内采样密度较低、分辨率较差的主要原因.

本文采用有限频率层析成像方法反演首都圈数字地震台网记录的远震资料获得该区上地幔三维P波速度结构模型, 并与射线层析成像结果进行了比较, 取得了以下认识:

(1)首都圈地区三个构造单元具有明显不同的速度结构特征, 其差异到150 km深仍然存在:燕山隆起区表现为相对的高速, 太行山隆起区整体以低速为主并存在小范围的高速块体, 华北盆地、渤海湾下的浅层上地幔存在大范围的强低速异常, 该异常的延伸深度在华北盆地下为50~150 km, 在渤海湾内则可达到200 km以下, 其顶面埋深为50~70 km, 与大地电磁测深并结合大地热流值数据推测的上地幔高导层的顶部埋深较为接近, 我们认为这是华北盆地、渤海湾下岩石圈明显减薄的显示.另外, 此低速异常下还隐约可见高速异常的分布, 随着深度的增加各构造单元的速度结构差异逐渐减弱.

(2)张家口-蓬莱断裂带作为上地幔浅层速度结构的变化带以及岩石圈减薄的边界带, 区内大部分历史强震都发生在这条构造边界带上, 由此看来该带上强震的发生不仅与地壳结构的不均匀性有关, 还可能与上地幔速度结构相关, 有着较深的构造背景.

(3)与射线层析成像方法相比, 有限频率层析成像方法的反演结果受数据采样的影响较大, 在数据覆盖和分辨率都较好的区域, 两种方法得到的结果具有大体一致的趋势, 但在有重要构造意义的华北盆地、渤海湾下速度结构特征的差异依然可见, 在这些地区缺少地震台站或宽频带地震台站可能是最为主要的原因.所以在现有资料基础上进一步补充应用宽频带资料, 可望获得更高分辨率的有限频率层析成像三维上地幔速度模型.

美国罗德岛大学YangShen教授和国立台湾大学洪淑蕙教授为本研究提供了有限频率层析成像方法的程序, 中国地震台网中心提供了远震波形数据, 匿名审稿专家对本文提出了有益的建议, 在此一并致谢!