2017, Vol. 37
Steketee[4]最早把位错理论引入地震学，假设均匀平面半空间地球介质内存在一个位错面并伴随一定的错动，并依此导出垂直走滑位错对应的位移格林函数。随后众多学者依据平面半空间地球模型研究了同震变形问题[5-8]。Okada[9]和Okubo[10-11]给出了完整、简洁、易于编程的地震位错引起的地表同震形变计算公式，适用于计算任意剪切与张裂位错引起的同震形变(含位移、应变、大地水准面变化、重力变化)。Sun等[12-15]和Pollitz[16-17]基于球对称地球模型分别提出球体位错理论，可计算任意类型地震在全球任意位置引起的同震形变，完成了从平面半空间位错理论向球体位错理论的转变。
在弹性位错理论的基础上考虑地球粘滞性构造的影响，即可研究地幔粘滞性松弛效应在地表引起的震后形变，即所谓的大地震弛豫效应。Wang等[18]基于自重可压缩层状半无限空间地球模型，提出一个研究震后形变的位错理论，并开发出一套可以计算地震引起的同震和震后位移、重力变化、应变以及应力变化的软件PSGRN/PSCMP，并应用于震后形变研究。Tanaka等[1-2]基于自重可压缩层状球体地球模型研究大地震后地幔粘滞性松弛效应导致的震后形变与重力变化，利用复平面路径积分的方法解决了地球只能有限分层的缺陷，并利用Poster-Widder公式保证了复平面数值积分的精度。该粘弹性球体位错理论同时考虑了地球的曲率、自重、可压缩性和层状结构，是迄今为止较为完备的震后形变理论之一。但Tanaka等编制的震后形变计算程序很难学，且近场结果中包含较大的计算误差，制约了该理论的推广。本文依据Tanaka的粘弹性球体位错理论[1-2]，编制出简便易学、便于推广的计算程序。新程序克服了原有程序近场计算精度不足的缺陷，可望得到广泛应用。
1 震后形变计算软件的设计思路弹性球体位错理论[12-15]和粘弹性球体位错理论[1-2]都是点源位错理论，其基本思想是把位错源视为尺寸可忽略的点源，并在此基础上推导出位错格林函数，建立位错模型与地表变形之间的关系。但实际地震产生的破裂面皆有一定尺寸，且破裂面上的位错量一般随空间变化，因此不可利用点源位错理论直接求解。为了解决点源位错理论和位错量随空间变化的实际震源之间的矛盾，一般做法是，先把整个破裂面分解成一定数量的子断层(图 1(a))，各子断层分别给定不同的位错量和位错方向，即可近似描述上述地震位错在破裂面上的空间变化(地震位错模型)；然后把每个子断层视为尺寸可忽略的点震源，即可通过点源位错理论计算该子断层位错在地表产生的形变；最后把整个破裂面上所有子断层位错在同一观测站产生的形变(含位移、重力变化、大地水准面变化)进行叠加，即可获得整个地震位错在该站产生的总体形变。上述编程思想即所谓的有限断层细分法[19]。
|
图 1 有限尺寸地震断层细分方案 Fig. 1 Model for a seismic fault |
Tanaka等根据上述有限断层细分法开发出一套计算软件，但该软件在具体应用时有两个方面的局限。其一，断层位错模型的子断层数目和大小一般事先给定，在具体形变计算时，没有考虑其尺寸是否可以忽略。研究表明，一个子断层是否可简化为点源，跟观测站与该子断层之间的相对位置有关。当观测站与子断层中心点之间的距离大于子断层边长的10倍时，将子断层视为点源进行计算可以保证精度，否则就会带来较大的计算误差，甚至错误[13, 20]。因此，对于一个具体的震例而言，当观测站距离断层面较近(子断层的震源距小于该子断层边长的10倍)时，利用Tanaka等的计算程序直接计算，就会出现较大的计算误差。其二，位错格林函数是震源深度的函数，对于不同的地震断层模型，各子断层的深度各不相同。因此，对于不同的震例，甚至同一地震的不同断层模型，位错格林函数都需要重新准备，才可计算该地震产生的地表形变。利用Tanaka等的计算软件，针对不同震例分别计算格林函数是一个相当繁杂且易出错的过程，导致其难以推广。
为克服Tanaka计算软件的上述局限，提高近场形变计算精度，并达到便于推广的目的，笔者开发出一套新的计算软件。该软件主要由3个部分组成：1)与32个震后时间点相匹配的32套离散格林函数数值框架，该框架基于PREM模型[21]建立；2)格林函数插值计算程序，可针对上述32套格林函数数值框架进行插值运算，输出任意震后时间点对应的格林函数数值结果；3)积分计算程序，调用上述格林函数数值结果，进行双二次样条插值运算，并对4类独立位错源对应的格林函数进行适当组合，计算出任意类型震源在地表产生的同震与震后形变(含位移、重力变化和大地水准面变化)。一般情况下，使用者只需按规定格式准备辅助文件，提供发震断层模型和观测点位置信息，即可运行积分计算程序进行同震与震后变形计算。
位错格林函数是地球介质对地震位错的响应，其计算是位错理论的核心。如上所述，利用Tanaka等设计的计算软件计算位错格林函数，是一个繁杂且易出错的过程。为解决该问题，本文设计了一个离散的格林函数数值框架(图 2)，任意深度的格林函数可根据该格林函数数值框架内插得到：基于PREM模型，针对震后某一具体时间点，利用Tanaka等的软件进行计算，获得16个深度的位错源对应的位错格林函数，组成一个离散的格林函数数值框架(图 2)。震源深度按照“浅层密集、深层稀疏”的原则构建，震源深度限定在0~200 km之间。同一深度的震源，按照“近场密集、远场稀疏”的原则采样。上述格林函数数值框架的设计原则，是根据位错格林函数近场(浅源)变化剧烈、远场(深源)变化平缓的特征确定。通过对上述16个深度的位错源对应的格林函数进行插值，即可获得0~200 km范围内任意深度的震源对应的位错格林函数。也就是说，震源深度低于200 km的震后位错格林函数(图 2中红色圆点)，可通过上述格林函数数值框架(图 2中的黑色圆点)内插得到。
|
图 2 二维离散位错格林函数数值框架示意图 Fig. 2 Sketch of a two-dimensional discrete frame of dislocation Green functions |
笔者开发的软件包包含32个震后时间点对应的格林函数(震后0、1、2、3、4、5、6、7、8、9、10、20、30、40、50、60、70、80、90、100、200、300、400、500、600、700、800、900、1 000、2 000、5 000和10 000 a)，每个震后时间点对应一套格林函数。在此基础上，依据震后形变随时间渐变的特征，设计一个插值计算程序，对上述32套格林函数进行插值运算，可获取任意震后时间点(比如1.5 a)对应的格林函数数值框架，然后运用积分计算程序进行相应的震后变形计算。软件可计算震后任意时间点对应的震后形变(含位移、重力变化、大地水准面变化)。本软件进一步完善了我们的前期研究成果[3]，省略了利用Tanaka等的程序自行计算位错格林函数的过程，达到简化计算过程、便于推广的目的。
此外，针对Tanaka等的计算程序近场计算精度不足的问题，在积分计算程序中设计一个自动判别功能。在计算某一个子断层位错在地表某一观测站引起的形变之前，先比较该观测站与子断层中心点之间的距离与子断层边长之间的关系。如果前者小于后者的10倍，程序就自动对该子断层进行细分(图 1(b))，然后分别计算细分后的微小位错在观测站产生的形变，最后利用子断层叠加法[19]进行处理，即可获得高精度的同震变形结果。
2 同震与震后位错格林函数的关系研究表明，对于球对称层状地球模型，如PREM [21]、1066A [22]，独立的位错源共有4个。为了方便表达，一般设定为水平走滑、水平引张、垂直引张和垂直倾滑位错(图 3)。应用上述独立点源位错对应的位错格林函数，可计算任意类型地震引起的同震与震后形变[1-2]。
|
图 3 4个独立震源，从左至右分别为水平走滑震源，垂直倾滑震源，水平引张震源和垂直引张震源 Fig. 3 Four independent seismic sources, from left to right are vertical strike-slip, vertical dip-slip, tensile on the vertical plane, and tensile on a horizontal plane, respectively |
当震后时间取为0时，对应的结果即为同震位错格林函数。为了探讨同震位错格林函数与震后位错格林函数之间的关系，以震源深度50 km的震源为例展开研究。假设地幔粘滞性因子为1× 1018 Pa·s，则震后0~10 a的重力变化格林函数如图 4所示, 图 4(a)、(b)、(c)和(d)分别为水平走滑、垂直倾滑、水平引张和垂直引张，震中距指其对数，黑色圆点为弹性球体位错理论格林函数[15]。由图 4可知，震后位错格林函数与同震位错格林函数在变化趋势上基本一致，其幅度随时间渐变。随着震后时间的推移，水平走滑震源和垂直引张震源引起的重力变化逐渐加大，但垂直倾滑震源和水平引张震源引起的变化逐渐减小。所有震源在某个时间点以后趋于稳定，不再变化。反之，当震后时刻靠近0时，震后格林函数逼近同震格林函数，反映出同震形变与震后形变在总体变化趋势上的一致性和渐变性。
|
图 4 4种独立震源的重力变化格林函数在震后0~10 a的变化趋势 Fig. 4 Dislocation Green functions of gravity changes for four independent seismic sources, the time that passes by after the event is set to be 0, 0.5, 1, 2, 4, 6, 8, and 10 years, respectively |
设震源周围地区地幔粘滞性因子为1×1019 Pa·s，震源深度为50 km，利用插值计算程序计算出震后4.5 a的格林函数数值结果，如图 5中红色线条所示。同时，利用Tanaka等的程序，计算震后4.5 a的格林函数数值结果，在图 5中以黑色三角形的形式离散给出。两种独立算法给出的格林函数完全一致，表明本文格林函数插值计算程序的正确性。
|
图 5 Tanaka等的程序直接计算得到的重力变化格林函数(黑色三角形)与插值算法得到的格林函数(红色实线)的比较 Fig. 5 Comparison of the interpolated (red curve) and calculated (black triangles) dislocation Green functions of gravity changes |
以2011年日本MW9.0地震引起的同震位移为例，与弹性位错理论配套计算程序[20]进行比较。选用Wei等[23]给出的地震断层模型，计算点为中国大陆的部分GPS观测站[24]。将震后时间设置为0，运行积分程序即可计算同震位移，配套计算程序同样可得到该地震的同震位移结果。图 6给出了2011年日本MW9.0地震在中国大陆部分GPS站引起的同震水平位移，图 6(a)中红色箭头为新程序的结果，蓝色箭头为配套计算程序[20]的结果，图 6(b)为2套理论结果之间的差异。经比较可知计算结果基本一致，二者的差异总体上只有同震信号的1%左右，应归因于2套理论基于的地球模型的差异以及计算误差所致。2套理论结果的一致性验证了新程序的正确性。
|
图 6 理论同震水平位移结果的比较 Fig. 6 Comparison of the results of theoretical co-seismic horizontal displacements |
此外，计算结果与观测结果的比较也是验证程序有效性的重要手段。图 7给出了中国大陆GPS观测到的2011年日本MW9.0地震的同震水平位移[24]，以及相应的理论计算结果, 红色箭头为理论结果，蓝色箭头为GPS观测结果。二者的大体一致性同样证实了新程序的正确性。总体上，理论计算结果在方向上与GPS观测结果一致，但位移量较观测结果稍小，这应该是地震断层模型[23]的地震矩偏小所致。
|
图 7 新程序计算结果与GPS观测结果的比较 Fig. 7 Comparison of the results calculated using our new code with GPS observations |
程序包中给出了32套位错格林函数数值框架，对应的地幔粘滞性系数为1×1020 Pa·s，震后时间点分别为0、1、2、3、4、5、6、7、8、9、10、20、30、40、50、60、70、80、90、100、200、300、400、500、600、700、800、900、1 000、2 000、5 000和10 000 a。如果震源附近的实际地幔粘滞性系数是1×1020 Pa·s，可以任选一套与上述震后时间点对应的格林函数，直接进行震后形变计算。
考虑到地幔粘滞性导致的震后形变与震源所在地区地幔粘滞性系数相关，开发了一个格林函数插值计算程序，用于计算目标震后时间点和特定地幔粘滞性系数对应的格林函数数值结果。格林函数插值计算程序的2个输入参数为震源周边地区实际地幔粘滞性因子与1×1019 Pa·s之间的比值，以及目标震后时间点(如3.5 a、10.4 a等)。具体地讲，如果实际地幔粘滞性因子为4.5×1019 Pa·s，计算震后3.5 a的震后格林函数，则第一个输入参数为4.5，第二个输入参数为3.5，即可获得对应的格林函数数值框架。然后运行积分计算程序，调用插值计算得到的格林函数数值框架，即可计算相应的震后形变。
此外，对于插值计算程序而言，内插的结果可靠，外插的结果则无法保证精度。因此，作为格林函数插值计算程序的2个输入参数有如下限制：2个参数皆为非负数，且其乘积必须小于或者等于1 000。
5.2 关于积分计算程序输入文件的格式及其代表的物理含义积分计算程序的2个输入文件分别为fault.dat(地震断层模型信息)与observation_point.dat(观测站的经纬度信息)，示例分别见表 1和表 2。使用者需根据不同的震例，修改2个文件的内容，即可运行积分计算程序，展开针对性的计算。
|
表 1 断层模型文件(fault.dat) Tab. 1 Fault model file (fault.dat) |
|
|
表 2 观测点位文件(observation_point.dat) Tab. 2 Observation point file (observation_point.dat) |
表 1给出了地震断层模型信息，一般依据大地测量、地震波等数据反演得到。为了表达断层位错在断层面上的空间变化，一般做法是把整个断层面按一定规律分解成一定数量的子断层，各子断层上的位错可近似认为一致。表 1中，第1行为说明信息，介绍各分量的物理含义；第2行为第1个子断层的几何模型及其上盘相对于下盘的滑动信息，第(N+1)行为第N个子断层的几何模型及其滑动信息，依次按此格式顺序排列。具体地说，Lat.和Longi.对应子断层中心点的纬度和经度，以(°)为单位，经度的取值范围为0~360，纬度为-90~90；Dislo.指子断层的相对位错量，以m为单位；Slip是位错面上盘相对于下盘的滑动角，即上盘的滑动角与水平右方向之间的夹角(从上盘方向看)，以(°)为单位；Dep.是子断层中心点的深度，以km为单位；Azim.是子断层的方位角，以(°)为单位；Dip是子断层的倾角，以(°)为单位；Length和Width分别为子矩形断层的长度和宽度，以km为单位。文件最后一行的99999是文件结束的标志。结束标志务必保留，否则积分计算程序无法正常运行。
表 2给出观测站(计算点)的经纬度信息，单位为(°)，经度的取值范围为0~360，纬度为-90~90，需要计算的各观测站的经纬度信息依照表 2的形式依次排列。同样，表 2的第1行为说明信息，第2行为第1个观测站的经纬度，第(N+1)行为第N个观测站的经纬度，最后一行为务必保留的99999为结束标志。
5.3 关于震源深度的限定地震断层模型各子断层中心点的深度务必小于200 km。由于位错格林函数数值框架的深度确定在0~200 km之间(图 2)，目标震源的深度必须限定在200 km以内。如果表 1中任意子断层的震源深度超过200 km，将会出现计算错误。
5.4 关于同震与震后计算结果的表达本程序包用方位角和位移量2个参数表达同震与震后水平位移计算结果，输出文件为DISPLACEMENT.dat(表 3)。方位角是指水平位移方向与正东方向的夹角，即从正东方向逆时针旋转到位移方向。大地水准面变化与重力变化为文件的第3列和第4列，前2列为观测点位置信息，与表 3类似。位移与大地水准面变化的单位为mm，重力变化的单位为μGal
|
|
表 3 同震与震后位移计算结果(DISPLACEMENT.dat) Tab. 3 Computation results of co-seismic and post-seismic (DISPLACEMENT.dat) |
基于Tanaka等[1-2]的粘弹性球体位错理论，编制了一套界面友好的计算程序，使用者只需按要求准备辅助文件，提供发震断层模型和计算点位信息，以及震源周边地区的地幔粘滞性因子，即可运行积分计算程序，计算任意地震在地表任意位置产生的同震与震后形变(含位移、重力变化、大地水准面变化)。
相对于Tanaka等[1-2]的计算程序，本文的新软件具有两大优势：可根据具体情况对断层模型进行自动细分，然后利用子断层叠加法[19]进行计算，明显提高了近场形变的计算精度；不用针对具体震例分别计算位错格林函数，极大地简化了计算过程的复杂度，因此具备推广的可能性。在我们前期研究成果[3]的基础上进行了进一步改进，可计算震后任意时间点对应的震后形变。
利用本文的新软件与弹性位错理论配套计算程序[20]，分别计算了日本MW9.0地震在中国大陆引起的远场同震位移，2套程序计算的位移结果之间的差异仅占信号的1%左右。此外，计算结果与GPS观测结果也大体一致，证明了新软件的正确性。
致谢: 感谢东京大学地震研究所田中爱幸博士提供位错格林函数计算程序，为本研究的位错格林函数数值框架的建立奠定了基础。
| [1] |
Tanaka Y, Okuno J, Okubo S. A New Method for the Computation of Global Viscoelastic Post-Seismic Deformation in a Realistic Earth Model (Ⅰ):Vertical Displacement and Gravity Variation[J]. Geophys J Int, 2006, 164(2): 273-289 DOI:10.1111/gji.2006.164.issue-2
(  0) |
| [2] |
Tanaka Y, Okuno J, Okubo S. A New Method for the Computation of Global Viscoelastic Post-Seismic Deformation in a Realistic Earth Model (Ⅱ):Horizontal Displacement[J]. Geophys J Int, 2007, 170(3): 1 031-1 052 DOI:10.1111/gji.2007.170.issue-3
( 0) |
| [3] |
Gao S H, Fu G Y, Liu T, et al. A New Code for Calculating Post-Seismic Displacements as Well as Geoid and Gravity Changes on a Layered Visco-Elastic Spherical Earth[J]. Pure and Applied Geophysics, 2017, 174(3): 1 167-1 180 DOI:10.1007/s00024-016-1453-2
( 0) |
| [4] |
Steketee J A. On Volterra's Dislocations in a Semi-Infinite Elastic Medium[J]. Can J Phys, 1958, 36(2): 192-205 DOI:10.1139/p58-024
( 0) |
| [5] |
Chinnery M A. The Deformation of Ground around Surface Faults[J]. Bull Seism Soc Am, 1961, 51(3): 355-372
( 0) |
| [6] |
Chinnery M A. The Stress Changes That Accompany Strike-Slip Faulting[J]. Bull Seism Soc Am, 1963, 53(5): 921-932
( 0) |
| [7] |
Maruyama T. Statical Elastic Dislocations in an Infinite and Semi-Infinite Medium[J]. Bull Earthq Res Inst, 1964, 42(2): 289-368
( 0) |
| [8] |
Matsu'ura M, Iwasaki T. Study on Coseismic and Postseismic Crustal Movements Associated with the 1923 Kanto Earthquake[J]. Tectonohysics, 1983, 97(1-4): 201-215 DOI:10.1016/0040-1951(83)90148-8
( 0) |
| [9] |
Okada Y. Surface Deformation Due to Shear and Tensile Faults in a Half-Space[J]. Bull Seism Soc Am, 1985, 75(4): 1 135-1 154
( 0) |
| [10] |
Okubo S. Potential and Gravity Changes Raised by Point Dislocations[J]. Geophys J Int, 1991, 105(3): 573-586 DOI:10.1111/gji.1991.105.issue-3
( 0) |
| [11] |
Okubo S. Potential and Gravity Changes Due to Shear and Tensile Faults in a Half-Space[J]. J Geophys Res, 1992, 97(B5): 7 137-7 144 DOI:10.1029/92JB00178
( 0) |
| [12] |
Sun W K, Okubo S. Surface Potential and Gravity Changes due to Internal Dislocations in a Spherical Earth-Ⅰ:Theory for a Point Dislocation[J]. Geophys J Int, 1993, 114(3): 569-592 DOI:10.1111/gji.1993.114.issue-3
( 0) |
| [13] |
Sun W K, Okubo S. Surface Potential and Gravity Changes Due to Internal Dislocations in a Spherical Earth-Ⅱ:Application to a Finite Fault[J]. Geophys J Int, 1998, 132(1): 79-88
( 0) |
| [14] |
Sun W K, Okubo S, Fu G Y. Green's Functions of Co-Seismic Strain Changes and Investigation of Effect of Earth's Spherical Curvature and Radial Heterogeneity[J]. Geophys J Int, 2006, 167(3): 1 273-1 291 DOI:10.1111/gji.2006.167.issue-3
( 0) |
| [15] |
Sun W K, Okubo S, Fu G Y, et al. General Formulations of Global Co-Seismic Deformations Caused by an Arbitrary Dislocation in a Spherically Symmetric Earth Model-Applicable to Deformed Earth Surface and Space-Fixed Point[J]. Geophys J Int, 2009, 177(3): 817-833 DOI:10.1111/gji.2009.177.issue-3
( 0) |
| [16] |
Pollitz F F. Coseismic Deformation from Earthquake Faulting in a Layered Spherical Earth[J]. Geophys J Int, 1996, 125(1): 1-14
( 0) |
| [17] |
Pollitz F F. Gravitational Viscoelastic Post-Seismic Relaxation on a Layered Spherical Earth[J]. J Geophys Res, 1997, 102(B8): 17 921-17 941 DOI:10.1029/97JB01277
( 0) |
| [18] |
Wang R J, Lorenzo-Martin F, 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[J]. Computers & Geosciences, 2006, 32(4): 527-541
( 0) |
| [19] |
Fu G Y, W K. Effects of Spatial Distribution of Fault Slip on Calculating Co-Seismic Displacement: Case Studies of the Chi-Chi Earthquake (MW7.6) and the Kunlun Earthquake (MW7.8)[J]. Geophys Res Lett, 2004, 31(21) http://onlinelibrary.wiley.com/doi/10.1029/2004GL020841/full
( 0) |
| [20] |
付广裕, 孙文科. 球体位错理论计算程序的总体设计与具体实现[J]. 地震, 2012, 32(2): 73-87 (Fu Guangyu, Sun Wenke. Overall Design and Specific Structures of the Computing Codes for Coseismic Deformations on a Layered Spherical Earth[J]. Earthquake, 2012, 32(2): 73-87)
( 0) |
| [21] |
Dziewonski A M, Anderson D L. Preliminary Reference Earth Model[J]. Phys Earth Planet Inter, 1981, 25(4): 297-356 DOI:10.1016/0031-9201(81)90046-7
( 0) |
| [22] |
Gilbert F, Dziewonski A M. An Application of Normal Mode Theory to the Retrieval of Structure Parameters and Source Mechanisms from Seismic Spectra[J]. Phil Trans R Soc A, 1975, 278(1 280): 187-269
( 0) |
| [23] |
Wei S J, Graves R, Helmberger D, et al. Sources of Shaking and Flooding during the Tohoku-Oki Earthquake: A Mixture of Rupture Styles[J]. Earth Planet Sci Lett, 2012, 333-334(6): 91-100
( 0) |
| [24] |
Wang M, Li Q, Wang F, et al. Far-Field Coseismic Displacements Associated with the 2011 Tohoku-Oki Earthquake in Japan Observed by Global Positioning System[J]. Chinese Science Bulletin, 2011, 56(23): 2 419-2 425 DOI:10.1007/s11434-011-4588-7
( 0) |