在地震学研究中地震发生时产生的同震位移、应变以及应力变化有着非常重要的研究意义,同震位移可用于反演断层破裂过程,应力变化可用于地震成因以及地震触发等研究[1-3]. 地震发生后,周边某些断层上库仑应力增大,加大了断层破裂的几率,相反则减少了断层上的负载[4-7]. 在许多研究区域都找到了库仑应力与余震分布、地震触发的对应关系,例如1999年土耳其伊兹米特地震与杜赛地震等[8]. 因此,计算地震发生后的同震位移、应变、应力变化是一项重要的基础工作.
自从Steketee[9]首次将静力学位错理论引入地震学后,有关计算静态位移场的理论研究迅速发展起来. 根据介质模型的不同,主要有以下几类,首先是基于弹性半无限介质模型的计算方法[10-11],Okada[12-13]详细地评述了有关的研究工作,给出了半无限介质内位错点源和有限源静态位移场的解析解. 这种方法的优点是给出的结果为解析解,计算速度快,缺点是模型过于简单,特别是当地表存在低速层时,半无限介质模型给出的结果往往会存在较大误差[14]. 因此发展起来了基于均匀弹性层状介质的相关理论,Sato和Natsu′ura[14],Ben-Nenaham 和Singh[15],Singh[16],Wang等[17]给出了水平层介质中的计算方法. Israel等[18],Pollitz[19]给出了球体分层模型中的计算方法. 层状模型可以较好地模拟地球介质. 但模拟震后较长时间间隔内的位移、应变、应力变化时,黏弹性松弛效应的影响则不可忽略[20-22],弹性模型无法考虑这一点,而黏弹性层状介质模型可以更好地拟合这种震后形变. 随着计算能力的增强,利用有限元等方法,基于更复杂的横向不均匀介质模型计算应力变化,在近些年逐渐得到更广泛地应用[23-24]. 这些不同的模型有其自身的适用范围,可以根据不同的研究目的选择合适方法加以利用. 本文主要考虑均匀弹性水平层状介质模型下,静力学位移、应变、应力的计算.
在计算理论地震图时,Kennett[25-26]发展了反射和透射系数矩阵方法,Bouchon[27]发展了一种可以有效地计算波数积分的离散波数方法,Yao 和Harkrider[28]将两者相结合,利用广义反射透射系数矩阵和离散波数方法计算近场理论地震图. 谢小碧和姚振兴[29]把这种处理波动问题方法加以推广,用以处理静态问题,给出了均匀弹性分层介质中剪切位错点源在地面产生的静态位移场的解答,He 等[30]给出了拉张源所产生的静态位移场的解答. 但是他们均未给出相应的应变、应力的计算方法与结果. 在他们工作的基础上,我们进一步发展该方法,使之可以基于均匀分层弹性介质模型,根据广义反射透射系数矩阵和离散波数方法,计算同震应变、应力变化等.
谢小碧和姚振兴[29],He 等[30]计算同震位移时,均采用了简单的梯形积分方式. 积分间隔的选取将直接影响到计算精度与速度,较小的积分间隔才能得到较为精确的数值解,但是这样会增加计算时间,对于利用同震位移资料进行反演问题研究十分不利. 理论地震图计算中广泛采用的Filon 积分方法[31-33],可以使得在选取较大的积分间隔时,计算精度仍然很高. 在本文中,计算同震位移时,使用梯形积分与Filon积分相结合的方法,即提高了计算速度,又使得计算结果具有较高精度.
在本文主要考虑的均匀弹性水平层状介质模型中,前人研究主要存在以下一些问题:某些研究中震源类型并不全面,特别是拉张源的计算;某些研究中仅给出了位移的计算结果;某些研究中水平层的数量受到制约等[17]. Wang等[17]给出了一种全面的计算方法,使用了Thomson-Haskell矩阵进行计算,矩阵元素中存在导致溢出问题的指数项. Wang 等[17]通过一些方法,如正交归一化等避免这些问题的出现. 而广义反射透射系数矩阵摒弃了造成数值不稳定的指数项,解决了数值溢出等问题,因此我们利用广义反射透射系数矩阵方法计算同震位移、应变、应力变化.
2 方法 2.1 同震应变、应力以及库仑应力的计算根据先前的研究结果[29-30],在如图 1 所示的水平层状介质中,三分量的静力学位移解答可以表示成如下形式:
(1) |
其中wz,qr,vθ 分别表示垂向、径向和切向的静力学位移,Σ 为断层面元的面积,Δui,i=1,2,3分别对应走滑、倾滑和拉张断层的静力学位错的错距,Aim 为方向因子,qm,wm,vm分别为面谐矢量坐标系下的位移分量,Jm(kr),J′m(kr)分别为m阶贝塞尔函数及其导数.
根据位移解答求应变与应力的关键步骤为求位移分量的偏导数,即求解
对θ 的偏导数表示为
(2) |
根据方向性因子的表达式可得到相应的偏导数:
(3) |
(4) |
对r的偏导数表示为
(5) |
对z的偏导数表示为
(6) |
当接收点位于地面或介质内部时,qmwm和vm表达式各不相同. 当接收点位于地表时有
(7) |
其中RUFS 、TUFS 和RDSL 分别为在(z1,zs- )之间朝上入射时的广义反射、透射系数矩阵及在(zs+ ,zL+)之间向下入射时的广义反射系数矩阵. RU,LFS、TU,LFS 和RD,LSL 表示求vm时相应的量.
震源系数分别为
(8) |
此处Pm,SVm和SHm分别表示Pm+,Pm-,SVm+,SVm-和SHm+,SHm-,而
(9) |
Δ 为(λ+μ)/(λ+3μ).
当接收点位于介质内部时,分接收点的深度zR小于震源深度和大于震源深度两种情况. 当接收点位于震源以上时,即当zR <zs 时,有
(10) |
其中接收函数矩阵和接收函数分别为
(11) |
RUFR 为在(z1,zR-)之间向上入射的广义反射系数矩阵,而RDRS 和TURS 为在(zR+,zs+ )之间向下入射的广义反射系数矩阵和向上入射的广义透射系数矩阵.
当接收点位于震源以下时,即当zR >zs,有
(12) |
其中接收函数矩阵和接收函数分别为
(13) |
RUSR 和TDSR 为在(zs+ ,zR+)之间向上入射的广义反射系数矩阵和向下入射的广义透射系数矩阵,RRL D 为在(zR+,zL+)之间向下入射时的广义反射系数矩阵.
(14) |
其中err,eθθ,ezz,erθ,erz,eθz为应变分量,τrr,τθθ,τzz,τrθ,τrz,τθz为应力分量,urr,uθθ,uzz,urθ,urz,uθz为位移的偏导数,λ,μ 为弹性常数. 根据应力张量可以计算断层面上的库仑应力变化[4, 36-37],如下式所示:
(15) |
其中Δσ,Δτs,Δσn,μ′ 分别表示库仑应力,断层面上的剪应力,断层面上的正应力,以及等效摩擦系数.
2.2 同震位移计算中积分方式的改进同震位移计算中的关键为计算以下积分:
(16) |
为了提高计算速度,需要选取较大的积分间隔,如果此时仍使用梯形积分,将会带来较大的计算误差. 当kr较大时,Jm(kr)存在渐进展开式,可以采用Filon 积分方式,选取较大的积分间隔,仍会有较高的计算精度[31-33]. 因此在积分的计算过程中,采用分段积分的方式. 当kr较小时,选取较小的积分间隔,采用梯形积分. 当kr较大时,选取较大的积分间隔,采用纪晨和姚振兴改进的Filon积分方法[32],仅取贝塞尔函数Jm(kr)渐进展开式的零次项
(17) |
这样既保证了计算精度,又提高了计算速度.
3 数值检验为了检验公式的正确性及数值方法的精度,我们做了一系列的数值检验. 分别验证均匀弹性半无限介质模型中接收点位于不同深度时,以及分层介质模型中接收点位于地表时,点源产生的位移以及应力解的正确性. 在笛卡尔坐标系中,设定沿着断层走向为X轴正方向,逆时针旋转90°为Y轴正方向,垂直向上为Z轴正方向,断层左上角对应于地表的投影点为坐标原点,如图 2所示.
均匀弹性半无限介质参数设定为P波速度4.0km/s,S波速度2.7km/s,密度2.5g/cm3. 使用10层,每层厚度为3km,具有相同P 波速度,S 波速度以及密度的均匀弹性介质层,之下设置为半无限,用来模拟均匀弹性半空间. 震源参数设定如下:倾角为70°,断层面积0.01km×0.01km,震源深度5.0km. 接收点参数设定如下:接收点位于地表,在笛卡尔坐标系中,接收点横坐标3.0km,纵坐标从-10.0km 起始,到10.0km 截止,共取41个点. 分别计算了走滑断层、倾滑断层以及拉张断层所产生的位移以及应力结果,并将之与Okada给出的解析解[12-13]进行了对比. 图 3a,图 3b分别给出了位移,应力结果的对比图,由于接收点位于地表,相应的应力分量只有径向、切向正应力以及径向与切向之间的剪应力. 与解析解相比较,两者一致.
采用与3.1相同的介质设定,震源参数设定如下:倾角为70°,断层面积0.01km×0.01km,震源深度20.0km,接收点设定如下:接收点位于地下,深度16.0km,接收点横坐标3.0km,纵坐标从-10.0km 起始,到10.0km 截止,共取41 个点. 图 4a,图 4b 分别给出了位移,应力结果与Okada的解析解[12-13]的对比图,两者一致.
采用与3.1节相同的介质设定,震源参数设定如下:倾角为70°,断层面积0.01km×0.01km,震源深度10.0km,接收点设定如下:接收点位于地下,深度14.0km,接收点横坐标3.0km,纵坐标从-10.0km 起始,到10.0km 截止,共取41 个点. 图 5a,图 5b分别给出了位移,应力结果与Okada 的解析解[12-13]的对比图,两者一致. 通过以上数值检验可以看出,基于均匀弹性分层介质模型,根据广义反射透射矩阵和离散波数方法计算出的位移与应力结果与解析解一致,验证了公式的正确性以及计算的精度.
选取如表 1 所示的均匀弹性水平层状介质模型,将接收点设置于地表,计算点源产生的位移以及应力解,通过与Wang等(2003) 的计算结果比较,验证在分层介质中的数值计算精度. 震源参数设定如下:倾角为70°,断层面积0.01km×0.01km,震源深度5.0km. 接收点参数设定如下:接收点位于地表,在笛卡尔坐标系中,接收点横坐标3.0km,纵坐标从-10.0km 起始,到10.0km 截止,共取41个点. 计算了走滑断层与倾滑断层所产生的位移以及应力结果,图 6给出了位移,应力结果的对比图,由于接收点位于地表,相应的应力分量只有径向、切向正应力以及径向与切向之间的剪应力. 由图可见,两者一致.
根据计算出的同震应力张量,即可得出断层面上的静态库仑应力变化,用于地震触发研究[4]. 考虑一个断层,走向0°,倾角70°,滑动量1 m,滑动角45°,断层长度20km,断层宽度10km,计算库仑应力需要有断层面以及滑动方向,假设接收点处的断层面以及滑动方向与破裂断层一致,计算深度为地下5.0km,接收点横坐标从-20.0km 起始,到40.0km 截止,共取13个点,纵坐标从20.0km 起始,到40.0km 截止,共取5 个点,首先基于半无限介质模型,介质参数为P波速度6.2km/s,S波速度3.5km/s,密度2.9g/cm3,分别根据Okada[12-13]和本文的方法计算库仑应力变化. 图 7a分别给出了相应的计算结果,以及绝对误差,由图可见两者误差很小,再次验证了数值计算的精度. 考虑分层介质模型(表 2) ,计算得到库仑应力分布,图 7b给出其与半无限结果的比较图,两者在库仑应力大小上存在着显著差异,因此在计算应力触发时,考虑分层弹性半无限介质模型是有必要的.
为了提高同震位移的计算速度,利于反演问题研究,我们改进了积分方式. 通过不同的积分方法计算了半无限介质模型中点源产生的同震位移,考察了改进方式的优越性. 接收点设定于地表,坐标为(20km,250km). 震源参数设定为:震源深度7.0km,断层面积0.01km×0.01km,倾角为60°.
由计算结果(表 3) 可见,采用改进的积分方式(积分方式二)计算同震位移时,核函数计算次数减少了2/3左右,而误差均控在3%以内. 这样既保证了计算精度,又提高了计算速度,有利于反演问题研究.
本文进一步发展了基于均匀弹性水平层状介质,利用广义反射透射系数矩阵和离散波数计算同震位移的方法[29-30],使之可以计算相应的应变、应力以及同震库仑应力变化. 接收点可以位于地表以及地下,可以适应于剪切位错源和拉张源等. 根据应力张量计算同震库仑应力的变化,可用于地震触发研究,为判断余震分布提供参考. 谢小碧和姚振兴[29],He等[30]计算同震位移时,均采用了简单的梯形积分方式,本文使用了梯形积分与Filon 积分相结合的积分方式,既提高了同震位移计算的速度,又保证了计算精度,有利于反演问题研究.
[1] | Stein R S, King G C P, Lin J. Change in failure stress on the southern San Andreas fault system caused by the 1992 magnitude=7.4 Landers earthquake. Science , 1992, 258(5086): 1328-1332. DOI:10.1126/science.258.5086.1328 |
[2] | Bouchon M. The state of stress on some faults of the San Andreas system as inferred from near-field strong motion data. J. Geophys. Res. , 1997, 102(B6): 11731-11744. DOI:10.1029/97JB00623 |
[3] | 沈正康, 万永革, 甘卫军, 等. 东昆仑活动断裂带大地震之间的黏弹性应力触发研究. 地球物理学报 , 2003, 46(6): 786–795. Shen Z K, Wan Y G, Gan W J, et al. Viscoelastic triggering among large earthquakes along the east Kunlun fault system. Chinese J. Geophys. (in Chinese) (in Chinese) , 2003, 46(6): 786-795. |
[4] | King G C P, Stein R S, Lin J. Static stress changes and the triggering of earthquakes. Bull. Seism. Soc. Am. , 1994, 84(3): 935-953. |
[5] | 石耀霖, 曹建玲. 库仑应力计算及应用过程中若干问题的讨论—以汉川地震为例. 地球物理学报 , 2010, 53(1): 102–110. Shi Y L, Cao J L. Some aspects in static stress change calculation—case study on Wenchuan earthquake. Chinese J. Geophys. (in Chinese) (in Chinese) , 2010, 53(1): 102-110. |
[6] | Toda S, Stein R S, Reasenberg P A, et al. Stress transferred by the 1995 Mw=6.9 Kobe, Japan, shock: effect on aftershocks and future earthquake probabilities. J. Geophys. Res. , 1998, 103: 24543-24565. DOI:10.1029/98JB00765 |
[7] | 万永革. "地震静态应力触发"问题的研究. 北京: 中国地震局地球物理研究所, 2001. Wan Y G. The study about "Static stress triggering" (in Chinese). Beijing: Institute of Geophysics China Earthquake Administration, 2001. |
[8] | Barka A. The 17 August 1999 Izmit earthquake. Science , 1999, 258(5435): 1858-1859. |
[9] | Steketee J A. On volterra's dislocations in a semi-infinite elastic medium. Can. J. Phys. , 1958, 36(2): 192-205. DOI:10.1139/p58-024 |
[10] | Maruyama T. Static elastic dislocation in an infinite and semi-infinite medium. Bull. Earthquake Res. Inst. Tokyo Univ. , 1964, 42: 289-368. |
[11] | Press F. Displacements, strains, and tilts at teleseismic distances. J. Geophys. Res. , 1965, 70(10): 2395-2412. DOI:10.1029/JZ070i010p02395 |
[12] | Okada Y. Surface deformation due to shear and tensile faults in a half-space. Bull. Seism. Soc. Am. , 1985, 75(4): 1135-1154. |
[13] | Okada Y. Internal deformation due to shear and tensile faults in a half-space. Bull. Seism. Soc. Am. , 1992, 82: 1018-1040. |
[14] | Sato R, Matsu'ura M. Static deformations due to the fault spreading over several layers in a multilayered medium, I, displacement. J. Phys. Earth , 1973, 21(3): 227-249. DOI:10.4294/jpe1952.21.227 |
[15] | Ben-Menahem A, Singh S J. Multipolar elastic fields in a layered half space. Bull. Seism. Soc. Am. , 1968, 58(5): 1519-1572. |
[16] | Singh S J. Static deformation of a multilayered half-space by internal sources. J. Geophys. Res. , 1970, 75(17): 3257-3263. DOI:10.1029/JB075i017p03257 |
[17] | Wang R, Martin F L, Roth F. Computation of deformation induced by earthquakes in a multi-layered elastic crust: FORTRAN programs EDGRN/EDCMP. Comp. Geosci. , 2003, 29(2): 195-207. DOI:10.1016/S0098-3004(02)00111-5 |
[18] | Israel M, Ben-Menahem A, Singh S J. Residual deformation of real Earth models with application to the Chandler wobble. Geophysical Journal of the Royal Astronomical Society , 1973, 32(2): 219-247. DOI:10.1111/j.1365-246X.1973.tb06528.x |
[19] | Pollitz F F. Coseismic deformation from earthquake faulting on a layered spherical earth. Geophys. J. Int. , 1996, 125(1): 1-14. DOI:10.1111/gji.1996.125.issue-1 |
[20] | Nur A, Mavko G. Postseismic viscoelastic rebound. Science , 1974, 183(4121): 204-206. DOI:10.1126/science.183.4121.204 |
[21] | Wang R, Martin F L, Roth F. PSGRN/PSCMP—a new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory. Comp. Geosci. , 2006, 32(4): 527-541. DOI:10.1016/j.cageo.2005.08.006 |
[22] | Pollitz F F. Postseismic relaxation theory on the spherical earth. Bull. Seism. Soc. Am. , 1992, 82(1): 422-453. |
[23] | Masterlark T, Wang H F. Transient stress coupling between the 1992 Landers and 1999 Hector Mine, California, earthquakes. Bull. Seism. Soc. Am. , 2002, 92(2): 1470-1480. |
[24] | 张怀, 吴忠良, 张东宁, 等. 虚拟川滇—基于千万网格并行有限元计算的区域强震演化过程数值模型设计和构建. 中国科学D辑 , 2009, 39(3): 260–270. Zhang H, Wu Z L, Zhang D N, et al. Quasi Sichuan-Yunnan-Based on thousands of grid parallel finite element computations strong earthquake evolution process numerical model design and construction. Science Chian (D) (in Chinese) (in Chinese) , 2009, 39(3): 260-270. |
[25] | Kennett B L N. Reflections, rays, and reverberations. Bull. Seism. Soc. Am. , 1974, 64(6): 1685-1696. |
[26] | Kennett B L N. Seismic wave propagation in stratified media. New York: Cambridge U. Press, 1983 . |
[27] | Bouchon M. Discrete wave number representation of elastic wave fields in three-space dimensions. J. Geophys. Res. , 1979, 84(B7): 3609-3614. DOI:10.1029/JB084iB07p03609 |
[28] | Yao Z X, Harkrider D G. A generalized reflection transmission coefficient matrix and discrete wavenumber method for synthetic seismograms. Bull. Seism. Soc. Am. , 1983, 73(6A): 1685-1699. |
[29] | 谢小碧, 姚振兴. 计算分层介质中位错点源静态位移场的广义反射、透射系数矩阵和离散波数方法. 地球物理学报 , 1989, 32(3): 270–280. Xie X B, Yao Z X. A generalized reflection-transmition coefficient matrix method to calculate static displacement field of a stratified half-space by dislocation source. Chinese J. Geophys. (in Chinese) (in Chinese) , 1989, 32(3): 270-280. |
[30] | He Y M, Wang W M, Yao Z X. Static deformation due to shear and tensile faults in a layered half-space. Bull. Seism. Soc. Am. , 2003, 93(5): 2253-2263. DOI:10.1785/0120020136 |
[31] | Mallick S, Frazer L N. Practical aspects of reflectivity modeling. Geophysics , 1988, 55(10): 1355-1364. |
[32] | 纪晨, 姚振兴. 区域地震范围的宽频带理论地震图算法研究. 地球物理学报 , 1995, 38(4): 460–468. Ji C, Yao Z X. The study of the method for broadband regional synthetic seismogram. Chinese J. Geophys. (in Chinese) (in Chinese) , 1995, 38(4): 460-468. |
[33] | Chen X F, Zhang H M. An efficient method for computing Green's functions for a layered half-space at large epicentral distances. Bull. Seism. Soc. Am. , 2001, 91(4): 858-869. DOI:10.1785/0120000113 |
[34] | Cotton F, Coutant O. Dynamic stress variations due to shear faults in a plane-layered medium. Geophys. J. Int. , 1997, 128(3): 676-688. DOI:10.1111/gji.1997.128.issue-3 |
[35] | Aki K, Richards P G. Quantitative Seismology: Theory and Methods. San Francisco: W. H. Freeman and Co., 1980. |
[36] | 汪建军. 同震、震后和震间应力触发. 武汉: 武汉大学, 2010. Wang J J. Coseismic, Postseismic and Interseismic Stress Triggerings (in Chinese). Wuhan: Wuhan University, 2010. |
[37] | 周宇明, 单斌, 熊熊. 静态应力触发中影响库仑应力变化的参数敏感性分析. 大地测量与地球动力学 , 2008, 28(5): 21–26. Zhou Y M, Shan B, Xiong X. Parameters sensitivity analysis of Coulomb stress change in static stress triggering. Journal of Geodesy and Geodynamics (in Chinese) (in Chinese) , 2008, 28(5): 21-26. |