格尔木地震台初至波到时残差分析*

2014-12-25杨广华杨欢庆

地震研究 2014年1期

杨广华，严浩，张博，杨欢庆

(1.青海省地震局格尔木地震台，青海格尔木816000;2.青海省地震局，青海西宁810001)

0 引言

地震波是震源辐射的弹性波。地震波的走时特别是体波的走时主要取决于介质的密度 (ρ)，杨氏模量 (E)和泊松比 (σ)，当三者稳定时波速也是稳定的。随着震中距的增大，地震波经过的介质是非常复杂的，尤其是地球内部结构的不对称性，介质的不均匀性，地核、地幔的起伏，莫霍面以上地壳厚度差异，地震孕育区的活动等诸多因素使ρ、E、σ处于不稳定状态，地震波速发生变化，地震通过这些区域的实际走时随之改变，造成了地震实际到时与理论到时有一定时差，即到时残差或称为波速异常。1956年宫本贞夫通过研究福井7.2级地震发现，震前在震中距400～600 km的首波异常达到± (2～3)s，这是首次观测到大震前波速异常现象 (冯德益，1981)。随着地震学的发展，利用P波走时残差研究更加广泛，不仅更多的用来研究地球内部结构，还用地震前的波速异常来预测地震 (冯德益等，1980;顾瑾平，盛国英，1983;冯锐等，1976;姜秀娥，陈非比，1981;戴维乐，1991)。此后全球地震工作者对波速异常进行了大量研究 (刁桂苓等，2005;张小涛等，2006;邹振轩等，2006;王林瑛等，2008;李永莉等，2009;龙海英等，2011;姚家骏等，2012)。

目前地震参数计算时大多使用的理论走时表为全球范围内的IASPEI91标准地震震相表、J-B走时表和中国地区走时表。由于受区域地震波速的影响和使用大范围的理论走时计算，使所记录到的地震波观测到时与中国地震台网 (CSN)地震目录提供的发震时刻通过理论计算的到时残差差异较大，有的甚至相差6～7 s。本文使用SSDP地震处理软件，并采用J-B走时表分析计算了格尔木地震台记录到的1 608个地震的初至波到时残差，并发现残差随方位、震中距、纬度都有不同程度的变化。

1 资料选取和计算方法

本文选取2001年1月至2013年12月格尔木地震台记录到的国内外1 608个初至波清晰，国内及边境为M≥5.0，境外为M≥6.0，震中距0.5°～171.1°，震源深度1～651 km，发生在各个方位的中强地震。用λ表示经度，Ψ表示纬度，θ表示方位角，Δ表示震中距，H表示深度，其中λ、Ψ、θ、Δ均使用角度制 (°)。

到时残差为观测到时减去理论到时。初至波到时残差δ定义为

式中，Pi为观测到时，P'i为理论到时，i为地震序号。Pi是台站测得的初至震相到时，直接从地震图上量取，可以认为是观测到时。P'i是以观测到时(《中国地震台网地震目录》提供的发震时刻)减去理论到时 (以《中国地震台网地震目录》提供的震中深度、位置至台站的震中距为依据，依J-B走时表计算的走时)。用式 (1)计算每一个地震的残差，计算时为了消除明显的错误，舍去偏离长时间±6 s以上的数据。

平均残差及残差的标准差σ的公式分别为

标准化残差δ*计算公式为

Pearson相关分析是一种计算两个变量间线性相关的方法，相关系数是用来分析判断直线相关的方向和密切程度的一种统计分析指标，其计算方法是两个变量的协方差与两变量的标准差的乘积之比，计算公式为

多元线性回归分析是研究一个因变量和多个自变量间定量线性相关关系的一种统计分析方法，多元线性回归模型一般表达式为

其中，yi为因变量，即被解释变量，a为常数项，xi为自变量，即解释变量，ai为自变量系数，即回归系数，ε为随机误差项。

图1 残差 (a)和标准化残差 (b)的频次分布Fig.1 Frequency distribution of residual(a)and standardized residual(b)

2 残差相关性分析检验

2.1 正态分布检验

由式 (1)、(4)分别计算初至波到时残差δ和标准化残差δ*，用出现的统计次数绘制残差和标准化残差的频次分布，如图1所示，从图中可看出残差δ和标准化残差δ*均服从正态分布。图1a中，残差δ多数集中在δ∈［0，5.0］区间内，且在3.3 s附近频次最高，3.1 s附近次之，说明格尔木地震台记录的国内外中强地震初至波到时比理论到时平均慢2.3 s，存在一定的偏差;图1b中，对残差标准化后，即标准化残差δ多数集中在δ∈［-3，3］区间内，在0 s附近频次最高，标准化后地震初至波实际到时更接近理论到时。

2.2 残差影响因素相关性分析

由式 (5)计算残差δ与λ、Ψ、θ、Δ、H间Pearson相关系数，具体结果见表1。由r(δ，Ψ)=0.429可知δ与Ψ是正相关关系，同时显著性 (双侧)结果为0，由于在0.01水平 (双侧)上显著相关，且显著性＜0.01，得出δ和Ψ间呈正线性相关性;r(δ，θ)=-0.262 和 r(δ，Δ)=-0.117 可知，δ与θ、Δ呈负相关关系，同时显著性 (双侧)结果都为0，得出δ与θ、Δ间有负线性相关性;r(δ，λ)=0.003和 r(δ，H)=0.010，r ＜ 0.3，且显著性 (双侧)结果都大于0.01，表示δ与λ、H间基本无线性相关关系。

表1 残差与影响因素间的Pearson相关系数Tab.1 Pearson correlation coefficient between residuals and influence factors

2.3 回归模型和回归系数的显著性分析

在95%置信水平下，选择进入线性回归方法，以δ为因变量，以Ψ、θ、Δ为自变量进入回归模型，得到多元线性回归模型的方差分析和回归系数，具体结果见表2。分析表2可得到，回归平方和为909.848，残差平方和为2 708.640，总平方和为3 618.488，F统计量的值为179.598，Sig.=0＜0.05(Sig.是显著性检验结果的检验值，越接近0越显著)，说明在给定置信水平下，回归方程能通过F检验，可以认为建立的回归方程有效。同时对其进行显著性检验发现，因变量δ对自变量Ψ、θ、Δ的回归非标准化回归系数分别为0.032、-0.003和0.011;对应的显著性检验的t值分别为19.584、-7.844和9.624，回归系数B的显著性水平Sig.=0，均小于0.05，可以认为自变量Ψ、θ、Δ对因变量δ均有显著影响。通过上述分析得到格尔木地震台记录的国内外中强地震初至波到时残差δ的回归方程为

抽取所选地震中的任何一个地震的参数代入式(7)，计算的残差结果都更加接近0，通过实际检验，表明建立的回归方程有效、可信。

表2 回归模型方差和系数Tab.2 Variance and coefficient of the regression model

3 计算结果及分析

3.1 残差与震中方位角的关系

表3以90°为间隔给出格尔木地震台不同方位上的平均残差值。由表3可见，残差受方位角和地区的制约较为明显。其原因除了与震源环境和介质有关外，主要与地震波的传播路径和介质的不均匀性有关。袁健和高龙生 (1990)通过P波到时残差的研究，得出中国东部地区的上地幔存在严重的横向不均匀性，似乎有一条北北东走向的深层构造带将中国东部地区分成两部分，西部速度高，东部速度偏低;薛峰和赵永 (1998)通过对中国速报台网走时残差分析，认为中国东部存在的欧亚板块与太平洋板块的交结带，是一个大的地震波阻隔带，它的存在会改变波的频谱，降低波速，减小波的幅度;贾淑娟 (1983)用P波走时残差研究深部结构，认为上地幔顶部的橄榄石晶粒在受力状态下会发生定向排列，这种定向排列的结果造成了地震波速度的各向异性，沿压力轴的方向速度增大，沿张应力轴的方向速度减少;张志军和任康 (1994)通过对地表和外空重力场的研究得出核幔存在较大的起伏，起伏幅度在±(1.7～±15)km之间，核幔形状的起伏，说明了地球深处并非处于流体静力平衡状态，相反，它存在着大尺度的动力过程，即在液核之内质量发生了迁移或流动，从而引起地震波走时异常。

从表3可以看出，发生在格尔木东北方向的地震初至波到时残差平均值最大为3.401 8 s，残差的标准差为±1.039，残差离散程度小，初至波残差整体偏大。说明在此方向的地球内部存在低速层或低速带，造成了地震波速变慢。同样发生在格尔木西南方向的地震初至波到时残差平均值为2.412 1 s，略大于总平均残差2.310 6 s，残差离散程度较小，初至波残差整体偏大。说明在此方向的地球内部也可能存在低速区，造成了地震波速变慢。在东南和西北方向的平均残差都小于总平均残差，且离散程度较大，说明这两个方向的初至波残差受低速层的影响小，主要是受地下介质和传播路径影响。

表3 δ-θ对应关系Tab.3 Corresponding relation ofδ-θ

3.2 残差与震中距的关系

理论研究和大量实际资料证明，地震波在岩层中传播速度与岩石地质年代、岩性、埋藏深度、密度、孔隙度、压力、温度等因素有关，或与岩石的弹性性质有关。一般来讲，随着震中距增大，路径就越复杂，走时残差变化越明显 (薛峰，赵永，1998)。为了分析这个问题，以10°为间隔计算震中距为0.5°～180°的初至波到时平均残差δ与震中距Δ的关系，具体结果见表4。从表4可以看出，在格尔木周围200 km范围内的地震初至波到时残差平均值最小为0.1 s，接近于理论值，地震波主要在地壳内传播，波速受介质的影响较大。震中距在2°～10°时，初至波到时残差平均值为3.461 2 s，地震波主要在地壳、沿莫霍面和上地幔传播，波速主要是受地壳厚度、莫霍面深度起伏变化和地壳内低速区的影响，因此在格尔木周围震中距在2°～10°内，某些区域的地壳或上地幔存在低速区。

震中距在37°～38°的地震共109个，平均残差最大为3.632 1 s，初至波到时大于理论到时。此震中距范围内，地震波的透射深度在900～1 000 km(国家地震局震害防御司，1992)。布伦地球分层模型把地幔分为上、下地幔，自莫霍界面至1 000 km左右为上地幔，1 000～2 900 km深度为下地幔。根据地震波的理论透射深度，初至波通过上、下地幔的交汇层，从格尔木地震台的观测资料显示，此时的初至波残差最大，因此说明，在上、下地幔的交汇处存在低速区或低速带，造成格尔木台初至波残差增大。

震中距在10°～40°范围内，残差随着震中距的增大而增大，在此范围内的地震波在地幔内传播，波速主要是受地幔内低速层和传播路径的影响。震中距在40°～180°范围内，残差随着震中距的增大而减小，地震波穿过下地幔和地核，波速受地幔低速区或低速带的影响明显减弱，尤其是震中距大于105°以后，残差平均值较小。

表4 δ-Δ对应关系Tab.4 Corresponding relation ofδ-Δ

3.3 残差与不同方位角的不同震中距的关系

在不同方位角的不同震中距初至波残差的明显变化，更能反映某地区介质的不均匀性和地震波传播路径的差异性。以10°为间隔计算在不同方位角θ，震中距0.5°～180°内发生地震的初至波平均残差δ与震中距Δ的关系见表5。由表5可以看出，格尔木地震台在不同方位角的初至波残差，随着震中距的变化有明显偏差。在0°＜θ＜90°和270°＜θ＜360°范围内，且 30°≤Δ ＜105°，初至波残差分别为3.600 7 s和3.939 5 s，在两个方位都为最大值，说明地震波受地球内部构造的影响非常大，震中距在此范围内的地震波主要在地幔传播，在残差与震中距的关系中已经说明在上、下地幔的交汇处存在低速区或低速带，对比90°＜θ＜180°和 180°＜ θ＜270°范围内，且 30°≤Δ ＜105°，初至波平均残差都比较小。通过以上分析表明，在格尔木以北地幔内存在东西走向的低速层或一条低速带，受此影响，两个方位角的地震初至波速变慢。

从表5可以看出，90°＜θ＜180°，0.9°≤Δ＜10°，残差值最大为 3.215 9;180°＜ θ＜270°，2°≤Δ＜20°，残差值最大为3.870 3。在这两个方位角，0.9°≤Δ＜20°内地震主要发生在青海、西藏、四川、甘肃4个省的一些地区，这些地区大多属于或靠近青藏高原，地壳厚度较大，地震初至波主要在地壳、沿莫霍面和上地幔传播。冯锐等(1981)利用地震面波研究了中国地壳结构，认为青藏和华北地区在地壳中的平均波速值较低，地壳内有低速层;吴建平等 (1998)用远震波形反演了高原整体的地幔深部P波速度模型，结果表明青藏高原地区的平均地壳厚度约68 km，上地幔Pn速度为8.15 km/s左右，P波的平均速度小于7.8 km/s，存在明显的上地幔P波低速异常区;丁志峰等 (2001)对青藏高原地震波三维速度结构的研究，认为青藏高原地壳内和上地幔存在明显的低速异常区;吴功建等 (1989)对青藏高原亚东—格尔木地学断面的研究认为青藏高原地壳的不均匀性显著，地壳厚度在70～80 km之间，壳内普遍存在低速层。因此，青藏高原地壳和地幔内存在低速异常区，从格尔木地震台的观测资料中也进一步反映出，发生在该地区的地震初至波到时残差较大，远远大于残差平均值。通过以上分析表明，发生在格尔木东南和西南方向地震的残差变化主要是受地壳厚度和低速异常区的影响，造成波速变慢残差增大。

表5 δ-θ-Δ对应关系Tab.5 Corresponding relation ofδ-θ-Δ

3.4 残差与纬度的关系

在残差影响因素相关性分析和回归系数的显著性分析中都表明了残差受纬度变化的相关性最为显著。为了反映初至波残差受纬度变化的影响，以15°为间隔计算了 (0°～90°)N 和 (0°～ -90°)N上发生地震的初至波平均残差δ与纬度φ的关系，具体见表6。由表6可以看出，北半球的初至波残差随着纬度的增大而增大，南半球的初至波残差随纬度变化不明显。

本文选取的1 608个地震，其中1 079个地震发生在北半球，占67.1%，震中距在0.5°～180°内地震的传播路径主要是沿大陆盖层路径传播，地震波在地壳、地幔和地核传播;529个地震发生在南半球，占32.9%，震中距36.4°～180°内地震的传播路径主要沿海洋盖层路径传播，地震波在地幔和地核内传播。在北半球格尔木地震台初至波残差受纬度φ变化而增大的主要原因是受传播路径的影响。由于青藏高原地壳巨厚、地壳和地幔存在低速区、格尔木以北地幔内存在东西走向的低速层或一条低速带，当地震波通过低速区或低速层、低速带时，波速减慢造成初至波残差增大。南半球发生的地震震中距都在36°以上，地震波主要在地幔、地核内传播，受地壳和区域低速层传播路径的影响较小，通过以上方位角和震中距的分析，方位角在90°～270°，震中距大于30°的地区，没有发现低速异常区，只是残差随震中距的增大而有所增大，残差在不同纬度的平均值都小于总的残差平均值，因此说明，发生在南半球的地震初至波残差，主要是随着震中距的增大，通过的介质呈现复杂，波速改变、残差增大。

表6 δ-Ψ对应关系Tab.6 Corresponding relation ofδ-Ψ

4 结论

本文通过对格尔木地震台记录到的国内外发生在不同方位、不同震中距、初至波清晰的1 608个M≥5地震，依据J-B走时表计算每个地震的理论到时，并与实际观测的初至波到时对比，发现残差差异较大;运用相关性分析和回归系数显著性分析对初至波残差的主要影响因素进行分析，建立回归方程，并对主要影响初至波残差变化较大的因素关系及造成的原因进行分析，得出如下结论:

(1)计算了格尔木地震台初至波残差平均值δ=2.310 6，说明格尔木地震台实际初至波到时比理论到时平均慢2.310 6 s。残差和标准化残差均服从正态分布。

(2)相关性分析结果得出，初至波残差与纬度、方位角、震中距之间有相关关系;与经度、深度之间基本无线性相关关系。说明格尔木地震台初至波残差变化受地震的纬度、方位角、震中距的影响。

(3)回归系数显著性分析结果得出，地震的纬度、方位角、震中距对初至波残差有显著影响，建立它们之间的回归方程为δ=1.734+0.011Δ-0.003θ+0.032Ψ，并通过实际检验，验证了所建立的回归方程有效、可信。

(4)对于发生在格尔木不同方位角的地震，统计结果得出，东北方向的地震初至波到时残差平均值最大为3.401 8 s，造成了地震波速变慢。

(5)对不同震中距的地震进行统计得出，震中距在2°～10°的地震，残差平均值为3.461 2 s，地震波速主要受地壳厚度、莫霍面深度起伏变化和地壳内低速层的影响，说明在格尔木周围某些区域的地壳或上地幔可能存在低速区;震中距在10°～40°的地震，残差随着震中距的增大而增大，地震波速主要是受地幔内低速层和传播路径的影响;震中距在37°～38°的地震，残差平均值最大为3.632 1 s，地震波速主要受地幔低速区或低速带的影响，说明上、下地幔的交汇处存在低速区或低速带;震中距在40°～180°的地震，残差随着震中距的增大而减小，地震波受地幔低速区或低速带的影响明显减弱。

(6)对不同方位角、不同震中距的地震进行统计得出，在不同方位的初至波残差，随着震中距的变化有明显偏差。东北和西北方位，30°≤Δ＜105°的地震，初至波残差最大，其主要原因可能是受地幔低速层的影响;西南和东南方位，0.9°≤Δ＜20°的地震，初至波残差最大，其主要原因可能是受地壳厚度、莫霍面起伏、低速区的影响。

(7)对于不同纬度的地震，统计结果得出初至波残差受纬度变化明显，在北半球残差随着纬度的增大而增大，其主要原因可能是受地壳厚度、地壳和地幔低速区的影响。

戴维乐.1991.菏泽5.9级地震前小震波速比时空变化特征［J］.西北地震学报，13(1):30-34.

刁桂苓，周仕勇，刘杰，等.2005.岫岩地震序列的平均波速比前兆异常［J］.华北地震科学，23(4):1-5.

丁志峰，何正勤，吴建平，等.2001.青藏高原地震波三维速度结构的研究［J］.中国地震，17(2):202-209.

冯德益，王周元，顾瑾平，等.1980.1976年四川省松潘—平武地震前波速比的异常变化［J］.地震学报，2(1):43-54.

冯德益.1981.地震波速异常［M］.北京:地震出版社.

冯锐，庞庆衍，傅征祥，等.1976.海城地震前后地震波速比的变化［J］.地球物理学报，19(4):295-305.

冯锐，朱介寿，丁韫玉，等.1981.利用地震面波研究中国地壳结构［J］.地震学报，3(4):335-340.

顾瑾平，盛国英.1983.多台求波速比方法的再探讨［J］.西北地震学报，5(4):95-98.

国家地震局震害防御司.1992.地震工作手册［M］.北京:地震出版社.

贾淑娟.1983.用P波走时残差研究深部结构［J］.山西地震，(1):17-22.

姜秀娥，陈非比.1981.用单台四震相法讨论唐山大震的波速异常［J］.地球物理学报，24(1):107-116.