在过去的几十年里, 北京、天津和唐山地区(以下简称京津唐地区), 由于其地理位置的重要性以及地质构造的复杂性, 一直是许多地震研究关注的地区.该地区的地壳和上地幔的三维速度结构迄今已有许多具有参考价值和科学意义的研究结果^[1~10].近年来, 国际上地震定位(相对定位和绝对定位)方法有显著的进展, 例如双差地震定位方法的提出及其广泛应用与成功验证^{[11, 12]}.在我国, 运用一维速度模型, 双差地震定位方法在地震定位方面亦得到了一些有参考价值和科学意义的研究结果^[13~16].

事实上, 双差地震定位方法可以与常规地震层析成像方法相结合, 既联合各自的优点, 又克服各自的缺点和局限性.基于上述思想, Zhang和Thurber^[17]提出了基于双差地震定位的地震层析成像方法, 简称双差地震层析成像方法.与常规地震层析成像方法相比, 用双差地震层析成像方法得到的三维速度结构模型更精确可信.在常规地震层析成像方法中, 由于绝对到时数据的观测误差的影响, 使得联合反演的震源位置会有某种程度上的散开现象, 而在双差地震层析成像方法中由于使用了相对到时数据则去除了这些误差, 从而避免了地震定位中出现的弥散现象, 继而在反演速度模型时避免了某些失真部分.

本文利用双差地震层析成像方法研究京津唐地区更加精细的三维速度结构和地震的精确定位, 借以探索该地区的三维速度结构和地震活动性之间的关系.

根据射线理论, 对于微扰问题, 震源i到地震观测台站k的地震波观测到时和理论到时之间的残差r^ik和地震波走时对震源和速度结构的微扰有如下关系:

(1)

式中, T_kⁱ是震源i到地震观测台站k的体波观测到时, τⁱ是震源i的初始发震时刻, u是波的慢度(速度的倒数), ds表示沿射线路径的线元, 震源坐标为(x₁, x₂, x³), 积分路径从震源i到观测台站k.在式(1)中, 初始发震时刻、射线路径和波的慢度(模型参数)均为未知.

同样, 对于震源j到地震观测台站k地震波观测到时和理论到时之间的残差r_k^j和方程(1)类似:

(2)

因此, 对于同一个地震观测台站k, 震源i和震源j之间有:

(3)

如果震源i和震源j彼此相距足够接近, 则震源到同一个观测地震台站k的射线路径可以认为是近似相同的, 此时若射线经过区域的速度结构为已知, 则方程(3)可进一步写为

(4)

式中, dr_k^ij称为双差(double-difference)^[11], 即为震源i和震源j的观测到时之差和理论到时之差的残差, 即

(5)

式(5)即为双差地震定位方法(简称DD法)^[11].式中, 观测到时差(T_kⁱ-T_k^j)^obs可以通过地震目录中绝对到时计算得到.方程(3)能简化成方程(4)是有其隐含的路径假定限制的, 即由地震对所在的震源区附近的速度不均匀性所引起的射线路径异常与震源对(i, j)的相对位置近似无关.这一假定在震源彼此之间相距很近时才成立, 而对于彼此相距较远的震源对则不成立, 在此情况下会造成相距较远的地震对的相对定位有很大的系统性偏差.

Zhang和Thurber^[17]提出了一种基于双差地震定位方法的地震层析成像方法, 简称双差地震层析成像方法.利用到时数据直接求解方程(3), 考虑地震对之间的路径异常, 即该方程右端第三项和第六项.若两个震源相距很近, 它们到同一个台站的射线路径基本上是重叠的, 此时方程(3)中震源区以外的模型偏导数部分可以互相抵消, 震源区以外的速度结构在反演中通过地震波的绝对到时数据来确定.这样, 既使用了相对到时数据, 又使用了绝对到时数据, 经过联合反演同时得到了研究区的三维速度结构模型、地震的绝对定位和地震的相对定位, 这是以往的常规地震层析成像方法所不能做到的.

鉴于现有资料和地震台站的分布情况, 本文的研究区域确定为114°E~120°E, 37.5°N~41.5°N (简称京津唐地区).该地区主要位于华北遥测地震台网和首都圈数字地震台网所覆盖的区域内.研究中共收集了1993年1月至2001年12月间华北遥测地震台网分布在京津唐地区的46个台站记录的、发生在研究区域内的近震直达P波初动到时数据, 和2002年1月至2004年8月间首都圈数字地震台网分布在京津唐地区105个地震台站(部分台站与华北遥测地震台网的台站重叠)记录到的发生在研究区域内的近震直达P波初动到时数据.将华北遥测地震台网和首都圈数字地震台网相互重叠的台站合并后, 两个地震台网分布在研究区域内的共有112个地震台站(图 1).华北遥测地震台网和首都圈数字地震台网提供的这些高质量、大批量的地震观测数据为京津唐地区的双差地震层析成像研究创造了很好的条件.

图 1

Fig. 1

图 1 研究区域内地震震源与台站分布图 (a)研究区域内地震震中与台站分布, 图中三角形表示地震台站, 正方形表示主要城市, 粗线表示活动断层, 圆圈表示本研究所用的地震的初始震中; b)研究区域内地震震源沿纬度剖面的深度分布; (c)研究区域内地震震源沿经度剖面的深度分布. Fig. 1 Distribution of seismcc stations and earthquake hypocenters used in this study (a) Distribution of earthquake epicenters (circles) and seismcc stations (triangles). Squares show the major cities. Thick lines denote active faults; (b) Cross sectional viewof focal depth along latitude profile; (c) Cross sectional viewof focal depth along longitude profile.

经过筛选, 共筛选出M1.0~6.2范围内3983次地震(每次地震的观测记录数至少6个, 图 1), 43408条直达P波初动绝对到时数据.在挑选出来的数据中, 每个地震对之间最大的空间距离为10 km, 每次地震最多可以和20次地震组成地震对, 每个震相对的到时残差小于0.5 s, 地震对的平均空间距离为5.3 km, 若地震对有6个以上震相对则认为是"强匹配地震对"其平均空间距离为5.0 km在筛选出的地震中, 震源深度的分布在地表到地表以下约30 km深, 从这些地震对的P波绝对到时数据中共挑选出200665条P波相对到时数据, 与P波绝对到时数据联合用于本文中的双差地震层析成像计算.P波初动到时数据的读数误差平均为±0.2~±0.3 s.由于唐山地区地震分布密集, 为了保持整个研究区域内的地震分布的均衡, 对该地区的地震做了更严格的筛选, 即每次地震至少有8个观测记录, 使得筛选后的地震尽可能均匀地分布在整个研究区域内, 又保证在地震活动性较强的区域有足够多且记录足够好的地震.图 2是3983次地震在研究区内的二维射线分布图, 在5 km到15 km深处, 穿过北京地区的射线数约为4000~6000条, 而在天津和唐山地区约为300~5500条; 在20~30 km深处, 穿过北京地区的射线数最低约为2000条, 最高约为3500条; 天津和唐山地区约为3000条.可以看出, 京津唐地区是射线分布最为密集的地区.

图 2

Fig. 2

图 2 研究区域内射线二维分布情况 Fig. 2 Distribution of P-wave ray paths used in this study

我们根据所选用的地震资料的数量和质量、地震震中和台站分布情况确定节点间距, 以便使尽可能多的节点处有较好的分辨率.通过反复实验, 比较反演结果, 将研究区域中节点间距水平向按0.5°× 0.5°划分, 垂直向节点分别放置在0~40 km (步长5 km)处, 这样92%以上的节点周围有几百条以上射线穿过.初始一维速度模型(图 3)采用于湘伟等^[6]所得到的京津唐地区的一维初始速度模型, 35 km以下部分速度模型参考人工爆破方面的结果以及已有的地震层析成像方面的研究结果^[18].

图 3

Fig. 3

图 3 用于双差地震层析成像反演的一维初始速度模型 Fig. 3 The starting 1-D velocity models used for double difference tomographic inversion in this study

我们在双差地震层析成像中使用了两种类型的到时数据:也震波绝对到时数据和相对到时数据.为了把这两种数据联合在一个反演体系中, 在双差地震层析成像方法中使用了分级加权方法.在联合反演的前2次迭代过程中, 赋予绝对到时数据较大的加权, 以便建立一个"大"尺度的计算结果.具体的权重比为:绝对到时数据的权重为1, 相对到时数据的权重为0.1.在反演的后2次迭代过程中, 赋予相对到时数据更大的权重, 进一步改进地震定位和震源区附近的速度结构, 所选用的权重比为:相对到时数据的权重为1, 绝对到时数据的权重为0.1.

我们用3983次地震的43408条直达P波初动绝对到时数据做了检测板实验^[19], 对反演结果的分辨率进行了检测, 结果如图 4所示.图中蓝色代表正的速度异常, 红色代表负的速度异常, 扰动值取为正常值的±3%.由图可见, 研究区内各深度层上均能正确恢复正负速度的相对变化, 特别是研究区内的重点地区北京、天津和唐山, 由于有较多台站和事件的分布, 使得每个深度层上这些地区的分辨率更高.密集数字地震台网记录的大量高精度数据为反演结果的可靠性提供了良好的保证.

图 4

Fig. 4

图 4 研究区域内不同深度处的检测板结果分布 Fig. 4 Map view of results of the checkerboard resolution test for the P-wave structure at different depths in the study area

除检测板的检验外, 每个节点附近通过的射线数也可以作为衡量解的可靠性的一个估计.穿过全部节点的射线数分布如图 5所示, 有92%的节点有足够的射线数(射线数≥10)参与反演.穿过每个网格节点的射线数介于0~7000之间, 其中, 穿过北京、天津、唐山地区的节点的射线数大于3000.

图 5

Fig. 5

图 5 研究区域内不同深度的射线分布 Fig. 5 Distribution of the number of the P-wave rays passing through each grid node (hit counts) The depth of the layer is shown at the lower right corner of each map. The hit count scale is shown at the bottom.

经过四次迭代, 用双差地震层析成像方法得到了以百分比表示的京津唐地区各个深度的P波速度扰动(图 6).图 6显示了在不同深度处速度结构的横向变化, 可以看出在研究区域内的中上地壳中, P波速度结构存在明显的横向不均匀性.

图 6

Fig. 6

图 6 研究区域内各个深度层上的P波速度扰动 Fig. 6 P-wave velocity image at each depth slice (in percentage from the average velocity)

将上述结果与京津唐地区的地质构造图(图 1)进行对比后可以看出, 在5 km到10 km深度处, P波相对速度扰动异常带的走向与各构造单元内地表地质构造和活动断裂的方向具有很好的一致性, 总的来说浅层速度图像很好地反映了地表地质、地形的变化, 即山脉和隆起区的下部表现为高速异常, 凹陷和山间盆地等处的下部则表现为低速异常.这与以往的地震层析成像结果^[1~10]基本是一致的, 但是在P波相对速度的高速异常和低速异常的交界带地区, 本文的高、低速速度异常对比更为剧烈.唐山地区虽位于华北平原, 但却是张(家口)渤(海)地震带上一个隆起的块体^[20], 在上地壳层位上(图 6中5 km和10 km深度层)都显示为高速异常.在深部, 在35 km深度层上, 太行山区则表现为低速异常, 而在燕山和华北平原的大部地区呈现高速.这些结果与Pn波层析成像结果^[7~10]是类似的.

图 7为研究区域内沿115.5°E, 116.5°E, 117.5°E, 118.5°E (图 7a)与沿40.5°N, 39.8°N, 39.0°N (图 7b)垂直剖面上的P波相对速度扰动分布(以百分比表示).可以看出, 在浅部, P波相对速度高速异常区和低速异常区的交界地带很好地反映了隆起区与平原、盆地的交界, 如在115.5°E剖面上39.5°N处, 117.5°E剖面上40.0°N处(图 7a) 40.5°N剖面上118.0°E, 39.8°N剖面上118.0°E处(图 7b), 均为平原或凹陷区(低速异常区)与隆起区(高速异常区)的交界地区, 两者非常一致.此外, 从116.5°E, 117.5°E和118.5°E垂直剖面上(图 7a)可以看出, 北京、天津和唐山地区的下部在15~20 km深度有一个尺度约50~90 km的高速异常带, 与周围速度对比剧烈.在研究区域的西部, 在115.5°E剖面上(图 7a), 燕山隆起和平原地区在20~30 km深度都呈现出大范围的低速异常, 而唐山地区(118.5°E垂直剖面, 图 7a)在20~30 km深度范围内也呈现出低速异常, 并且与上层的速度结构对比剧烈.39.8°N剖面中(图 7b)唐山以东地区作为张渤地震带上一个隆起的块体表现为高速异常, 太行山隆起区下部的高速异常区从地表向东延伸直到北京下部的约20 km处, 与周围的速度结构形成强烈对比.

图 7

Fig. 7

图 7 研究区域沿不同经度(a)和不同纬度(b)的垂直剖面上的P波相对速度扰动分布 图中BJ表示北京, TJ表示天津, TS表示唐山, ZJK表示张家口.(a)中每个剖面的上边缘经度位置即为沿经度方向的P波相对速度扰动的垂直剖面, 从上到下依次为沿115. 5°E, 116. 5°E, 117. 5°E和118. 5°E的垂直剖面上的P波相对速度扰动分布(以百分比表示); (b)中每个剖面的上边缘纬度位置即为沿纬度方向的P波相对速度扰动的垂直剖面.从上到下依次为沿40. 5°N、39. 8°N和39. 0°N的垂直剖面上的P波相对速度扰动分布(以百分比表示). Fig. 7 Vertical cross-sections of P-wave velocity perturbations at each longitude (115. 5°E, 116. 5°E, 117. 5°Eand 118. 5°E) (a) and each latitude (40. 5°N, 39. 8°N and 39.0°N)(b) BJ, Beijing; TJ, Tianjin; TS, Tangshan; ZJK, Zhangjiakou

本文用双差地震层析成像方法联合反演了京津唐地区的三维P波速度和震源参数(绝对位置和相对位置), 在得到京津唐地区P波三维速度结构的同时, 也得到了京津唐地区所用的3983次地震中2809次地震的重新定位结果^[21].

用所得到的P波三维速度结构对地震做重新定位后, 地震的均方根残差由初始的1.2 s降为0.3 s, 有了显著的改善.在整个研究区域内, 地震重新定位后得到的2809次地震的震源位置测定误差在E-W方向为1.0 km, 在N-S方向为1.1km, 深度方向为15 km.

为进一步探讨地震活动性与P波速度结构之间的关系, 我们将用双差地震层析成像方法联合反演得到的震级M_L≥3.0, 震源深度在5~15 km范围内的地震的震源投影在10 km深度处的P波相对速度扰动分布图上(图 8a)将震级M_L≥3.0, 震源深度在10~20 km范围内的地震投影在15 km深度处的P波相对速度扰动图上(图 8b).根据臧绍先等^[2]对华北地区地震活动的统计分析, 该区大多数较大地震都发生在研究区的中上地壳, 主要在10~15 km深度范围内.我们将发生在研究区域内, 公元前780~2002年间震级M≥6.0的历史地震的震源位置分别投影在10 km深度和15 km深度处的P波相对速度扰动分布图上(图 8a, 8b).

图 8

Fig. 8

图 8 (a)研究区内地震重新定位后的震源(ML≥3.0, 5km≤深度≤15km, 灰色圆圈)和历史地震(M≥6.0, 白色五角星)在10 km深度处P波相对速度扰动(以百分比表示)图上的分布; (b)研究区内地震重新定位后的震源(ML≥3.0, 10km≤深度≤20km, 灰色圆圈)和历史地震(M≥6.0, 白色五角星)在15 km深度处P波相对速度扰动(以百分比表示)图上的分布 Fig. 8 (a) Map view of P-wave velocity perturbations (in percentage) and distribution of both relocated earthquake hypocenters (ML≥3.0, 5 km≤focal depth≤15 km, gray circles) and historic earthquake hypocenters (M≥6.0, white stars) at depth 10 km; (b) Map view of P-wave velocity perturbations (in percentage) and distribution of both relocated earthquake hypocenters (ML≥3.0, 10 km≤focal depth≤20 km, gray circles) and historic earthquake hypocenters (M≥6.0, white stars) at depth 15 km

由震源在10 km和15 km深度P波相对速度扰动图上的投影(图 8a, 8b)可以看出, 重新定位后, 大部分地震的震源投影位于P波低速异常区与高速异常区的对比强烈的交界地带, 较少地震的震源投影位于偏P波低速异常区的过渡带(黄色区域)上, 较多地震的震源投影位于偏P波高速异常区的过渡带(浅蓝色区域)上.历史地震的震源投影(白色五角星)也有类似的分布, 其中, 唐山地震、滦县地震、三河-平谷地震的震源投影均位于高速异常区的过渡带内.值得注意的是, 在10 km深度上(图 8a), 北京地区、唐山地区的M_L≥3.0的地震的震源位置投影全部沿着P波高速异常区的边缘分布, 震源的展布方向与高速异常带的走向非常一致.

根据研究区域内地震活动性的分布, 我们在研究区内做了三个剖面来探讨地震活动性和P波速度异常之间的关系(图 9).在各个剖面的P波相对速度扰动分布图中, 我们选取了震中位置距离该剖面在纬度上相差±0.5°范围内M_L≥2.5地震, 将其震源投影在该剖面上(图中灰色圆圈).

图 9

Fig. 9

图 9 研究区域内沿三条剖面的P波相对速度扰动和地震活动性分布, 剖面位置如(d)所示 Fig. 9 Distribution of P-wave velocity perturbations (in percentage) and seismicity along three vertical cross-sections, (d) the locations of cross-sections

在剖面AA'的P波速度扰动分布和地震震源分布图(图 9a)中, 0~10 km深度范围内P波速度异常和浅层的地表、地形特征很吻合.在10~20 km深度处沿唐山到天津一线存在一个沿水平方向尺度约100 km的高速异常区, 而在20~30 km范围内唐山及其以北地区则存在大范围的低速异常区, 这与Huang等^{[8, 9]}的研究结果相一致.在Huang等^{[8, 9]}的工作中, 使用了非天然地震和定点爆破资料, 在本文的工作中使用的全部是天然地震, 但是取得了相同的地震层析成像结果, 这说明:吏用双差地震层析成像方法提高了地震的定位精度, 更加精确的地震定位使得我们仅用天然地震数据也可以得到更好的地震层析成像结果.从剖面AA'的P波相对速度变化中可以看出, 唐山-河间-磁县地震活动断裂带在上地壳均处于不连续的低速异常区内, 而在10~20 km深度范围内则表现为一条断续的南浅北深的高速异常带, 在浅层的低速异常体中和深部的高速异常体边缘分布有地震.唐山地区地震的最大深度与该地区中下地壳的低速异常的上缘相一致.此外石家庄以北地区的震源深度在10~20 km的地震均分布在该深度范围内的低速异常的边缘地带

在沿着张家口-渤海地震带所做的北西西向剖面BB'(图 9b)中, 唐山地区仍然处于低、高速异常的交界地带, 其东南侧为高速异常区, 而西北侧为低速异常区, 但低速异常的幅度没有东南侧高速异常的幅度大.在唐山地区下部20~30 km范围内存在低速异常区, 低速异常的幅度达到6%, 该低速异常的顶部深度西北深东南浅.唐山地区的地震的震源位置处于低、高速异常带的交界部位, 震源深度的变化也表现为西北深东南浅, 两者相一致.北京地区下部在该剖面上表现为扰动约3%的高速异常.在该剖面上, 重新定位后张北地区震级M_L≥4.0的地震震源最大深度在15 km左右, 位于张北地区下部的高速异常体边缘.

为了更清楚地研究唐山地区地震活动性与P波速度异常之间的关系, 我们将唐山地区也做了剖面(TT'剖面, 图 9c).于湘伟等^[21]通过将地震台网、两种常规地震层析成像、双差法和双差地震层析成像不同方法确定的唐山地区地震震源参数做了详细对比后得出, 双差地震层析成像得到的震源参数(地震绝对位置和地震相对位置)的精度高于常规地震层析成像方法和双差法得到的震源参数的精度, 指出唐山地区的地震沿北东方向剖面在深度上呈现明显的三个小震群的特点.将双差地震层析成像联合反演得到的唐山地区地震的震源投影在P波速度剖面TT'(图 9c)上, 这三个小震群分布在P波低、高速异常区的交界地带和高速体内.靠近唐山地区下部低速异常体的地震的震源深度变化与展布与低速异常体上边缘的深度变化和走向很一致.唐山地区东北侧迁安地区的地震的震源深度分布与该地区10~20 km深度范围内的高速异常带的走向很一致, 唐山地区的地震最大深度为25.4 km, 位于唐山地区下部20~30 km范围内低速异常区的上边缘, 其中震级M_L≥4.0的地震的最大震源深度为19.4km (±1.1km), 该地震(M_L=4.0)震源位置位于唐山地区下部的低速异常区的上边缘(浅黄色区域)(该地震在地震目录中给出的震源深度为10 km).在以往的地震层析成像结果中^{[8, 9]}, 唐山地区下部的低速异常带上边缘的深度为近20 km, 而在本文得到的唐山地区下部的低速异常带的上边缘要深于20 km, 约在25 km, 与该地区地震震源的最大深度相一致.由于双差地震层析成像反演中采用了相对到时数据, 改进了震源区的精细速度结构, 因而更精细地刻画了震源区的速度结构特征.

本文联合使用华北遥测地震台网和首都圈数字地震台网所记录的大量高精度区域地震的P波绝对到时数据和相对到时数据, 用双差地震层析成像方法联合反演了京津唐地区精细的地壳三维P波速度结构模型和地震的绝对定位和相对定位.由反演结果得到以下主要结论:

(1)联合反演得到的京津唐地区的地震绝对定位结果显示了研究区域内地震震源分布的更加精细的图像, 沿主要断裂带的地震的震中呈现出更清晰的条带状分布, 同时给出了部分地震台网未确定震源深度的地震的震源参数.

(2)京津唐地区的三维速度结构图像在浅层上很好地反映了地表地质、地形和岩性的特征.在平原和凹陷的盆地处呈现P波低速速度异常, 而在隆起的山区或基岩出露区显示为P波高速速度异常.

(3)研究区域内震级M≥6.0历史地震的震源位置和重新定位后震级M_L≥3.0的地震在10 km深度和15 km深度处的P波相对速度扰动图上的投影都显示出相似的特点.在所统计的地震中, 绝大部分的地震的震源位置投影分布在低、高速异常的交界地带, 且偏高速体一侧, 极少数的地震分布在P波速度异常体内部.

(4)唐山地区三个小震群均分布在P波低、高速异常区的交界地带和高速体内, 地震最大深度为25.4 km, 位于唐山地区下部20~30 km范围内低速异常区的上边缘, 两者非常吻合.说明双差地震层析成像反演中采用了相对到时数据, 改进了震源区的精细速度结构, 更精细地刻画了震源区的速度结构特征.

感谢郑秀芬研究员在数据分析上给予的帮助和周元泽副教授及两位审稿专家对本文提出的宝贵意见.