2019, Vol. 39
2. 中国地震局地壳应力研究所武汉科技创新基地，武汉市洪山侧路40号，430071;
3. 中国地震局第二监测中心，西安市西影路316 号，710054
2008-05-12发生在青藏高原东缘龙门山推覆构造带中段的MW7.9汶川地震，是一次罕见的板内逆冲斜滑型特大地震[1]。该次事件自西南向北东方向破裂长达320 km，揭示出有2条近平行的破裂带：破裂出露地表长达240 km的北川断层(BCF)和地表破裂长约90 km的彭灌断层(PGF)。
汶川地震后，通过地震、大地测量或联合反演, 学者们提出不同的同震滑动模型[2-4]，所有这些模型都表明了汶川同震破裂的复杂性。同震造成的应力重新分配会诱发一系列导致可见地表形变的震后形变过程，如震后余滑、粘弹性松弛、孔隙回弹效应等，这3种震后形变机理作用的时空形变特征存在一定的差异，可以探测断层及深部岩石圈的流变性质[5]，对认识地震形变机理、评估地震风险具有重要的作用。不少学者应用不同时间跨度、不同空间展布的大地测量资料对汶川地震震后形变开展了初步研究[6-8]。这些研究采用震后短期的GPS、InSAR近场的观测资料初步揭示，震后粘弹性松弛效应在汶川地震震后形变中发挥了重要作用，而孔隙弹性回弹和震后余滑则占比有限[8]。尽管这些工作对川西高原下地壳上地幔的粘滞系数起到一定的约束作用，并初步揭示粘滞系数随时间变化的特征，但是青藏高原上地幔和下地壳流变性仍然没有得到很好的约束。
本文搜集中国大陆构造环境监测网络(简称陆态网络)震后7 a的GPS观测资料，以大地测量资料反演的同震滑动分布模型作为震后粘弹性松弛效应的应力驱动源，采用纵向分层的岩石圈结构，约束了川西高原的流变学参数。
1 震后GPS形变分析 1.1 GPS资料处理本文搜集了龙门山断裂带上盘59个GPS站点的观测数据，站点分布如图 1，其中红色三角形表示从1999年开始测量的震前布设的流动站，紫色三角形表示从2009年开始观测的震后布设的流动站，紫色五角星表示从2010年开始观测的震后布设的连续站；蓝色线为文献[3]采用的同震模型的地表迹线；虚线框为文献[4]反演的同震滑动模型区域，其中红色虚线框为超出文献[3]模型的破裂区域，红色方框为文献[4]反演的最大破裂集中区。采用GAMIT/GLOBK软件对上述GPS资料及中国境内和周边长期稳定的IGS站点联合解算，最终得到ITRF2008框架下的坐标时间序列，详细数据处理过程参见文献[9]。
|
图 1 龙门山断裂带构造背景、汶川地震分布及震后观测GPS站点分布 Fig. 1 The regional tectonic setting of the Longmenshan fault zone, the fault geometry of Wenchuan earthquake and distribution of GPS sites |
一般地，GPS坐标时间序列包含震间形变、同震及仪器变更引起的阶跃、震后形变及季节性信号等，其数学表达式为：
| $ \begin{gathered} {\text{y}}\left( {{t_i}} \right) = {a_1} + {a_2}t + {a_3}\sin \left( {2{\rm{π }}{t_i}} \right) + {a_4}\cos \left( {2{\rm{π }}{t_i}} \right) + \hfill \\ \;\;\;{a_5}\sin \left( {4{\rm{π }}{t_i}} \right) + {a_6}\cos \left( {4{\rm{π }}{t_i}} \right) + {a_7}H\left( {{t_i} - {t_{eq}}} \right) + \hfill \\ \;\;\;\;\;\;\;{a_8}H\left( {{t_i} - {t_{eq}}} \right){\text{ln}}\left( {1 + \left( {{t_i} - {t_{eq}}} \right)/{\tau _i}} \right) + \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\mathop \sum \limits_{j = 1}^N {d_j}\left( {H\left( {{t_i} - {t_j}} \right)} \right) \hfill \\ \end{gathered} $ | (1) |
式中，y(ti)是GPS站点在厉元ti时刻的位置，a1为常数偏移，a2是长期线性速度值，a3和a4描述1 a周期项振幅，a5和a6描述0.5 a周期项振幅，a7描述同震偏移项，a8为震后松弛振幅，H(ti-teq)为阶跃函数，dj为其他原因导致的偏移，τi为松弛时间常数。τi与a8存在较强的反相关性，经试验，τi对文中震后形变拟合影响可以忽略，故文中将τi固定为40 d。
本文针对GPS数据观测方式和起测时间的不同，分别采用插值和拟合两种方式获得各测站震间速度[9-10]，并按式(1)拟合得到各测站震后形变值。文中数据仅对连续观测站进行季节性信号估计，忽略区域站的季节性形变波动。且依据各测站起测时间的不同，最终采用2008~2015年、2009~2015年和2010~2015年共3个不同时间段的震后累积形变约束震后反演。
2 龙门山地区地壳结构反演本文采用文献[4]反演的同震滑动分布模型，该模型比其他滑动模型能更好地拟合GPS、InSAR等大地测量资料观测的同震形变。尽管文献[4]采用的流动GPS、三角点等观测的同震形变不可避免地包括早期震后形变，但即使考虑震后8个月的震后余滑，文献[4]反演的同震滑动依然能辨识滑脱层的凹凸体。为了验证结果，我们还用文献[3]的同震模型进行验证。
为了检验上地壳和中下地壳之间是否存在低速层，论文采用与文献[6]相似的地壳分层模型。地球分层模型如图 2，该模型认为低速层和中下地壳为粘弹性层，并且采用Maxwell流变体表示，其中黑色虚线表示横波波速，黑色实线表示纵波波速，红色实线表示固定为1.00×1020 Pa·s的地幔层的粘滞系数，红色虚线表示本文得到的最优粘弹性参数解。我们利用PSGRN/PSCMP程序[11]来计算川西高原岩石圈分层结构的低速层和中下地壳的粘弹性参数η1和η2。
|
图 2 松潘-甘孜地块岩石圈分层模型 Fig. 2 Schematic diagram of viscoelastic earth model representing the Songpan-Ganzi terrane |
利用本文所选上盘GPS站点观测得到的地表水平位移作为约束，采用格网搜索的方法搜索低速层和中下地壳粘滞系数η1和η2的最优值。搜索范围在1.0×1017~1.0×1021 Pa·s之间，最优参数对应GPS观测值与模拟值拟合残差平方和最小值处。图 3显示了粘弹性松弛模型的模型参数(低速层和中下地壳粘滞系数)误差，红色五角星表示搜索的最优值。本文所选模型约束得到的最优粘弹性参数显示，低速层的粘滞系数为2.51×1018 Pa·s，而中下地壳的粘滞系数为3.98×1018 Pa·s。
|
图 3 巴颜喀拉块体低速层和中下地壳最优拟合模型残差 Fig. 3 Data misfit for a range of viscoelastic forward models composed of viscoelastic thin low-velocity zone (η1), overlying the viscoelastic lower crust with viscosity (η2) |
图 4显示基于最优粘弹性参数模型正演计算的震后三维形变与GPS观测值，图中红色三角形、紫色五角星、紫色三角形分别表示3种不同的观测时段，分别为2008~2015年、2010~2015年、2009~2015年，彩色圆圈表示水准数据，彩色六边形表示GPS累积垂直形变。川西高原的粘弹性松弛效应可以较好地解释中远场GPS观测的水平震后形变，如龙日坝断裂以西及鲜水河断裂带沿线的GPS震后形变。但是在近场拟合的形变场和实测值存在一定的方向偏差，如靠近断层近场的SCMX(茂县)和J417两个测站。产生这种偏差主要是由于近场受到震后余滑的作用比较大，此外同震滑动分布模型可能对方向偏差具有一定影响。
|
图 4 震后形变观测值与粘弹性模型模拟值的比较 Fig. 4 Comparison between GPS observed postseismic displacements and modeled postseismic displacement from the best-fitting viscoelastic relaxation model |
由于GPS垂直形变误差较大，因而在确定最优粘滞系数时仅采用GPS水平观测值。但从另一方面看，垂直形变在震后形变机理区分中可以发挥重要的作用[5, 12]。为检验模型的合理性，本文基于文中最优粘弹性模型计算了部分GPS区域站的垂直形变，这些测站都具有较长的观测历史，并且垂直形变幅度较大，解算结果比较可靠。由图 4可见，GPS观测垂直形变主要以隆升为主，基本落在粘弹性松弛模拟值隆升区域的内部和正负交替区域。陆态网络茂县基准站自2010年以来，年均垂直速率约7.0 mm/a，考虑到震后形变幅度随时间而衰减的特性，该测站震后前2 a的震后形变至少14 mm。此外，我们还收集了同时段的水准数据[13]，结果显示，沿着龙门山断裂带中北段的水准点均观测到不同程度的隆升，最大值隆升达50 mm，该结果与顾及侧向非均匀的地球结构模拟的垂直形变相吻合[14]。稀疏的GPS垂直形变与近场的水准观测均显示，龙门山断裂带上盘存在显著的震后隆升，但由于观测噪声较大，本文未采用这些信息约束深部物质的流变结构。
3.2 同震模型差异性分析图 4所示的震后形变模拟值与观测值在中远场总体吻合，但可以发现，在断层东北端存在一定的偏差。青川(QICH)站点的震后计算值和观测值之间存在明显的方位角偏差。当采用沈正康等[5]的同震模型计算震后形变值时发现，在青川站点可以很好地拟合观测值(图 4红色虚线框中棕色向量)，说明同震模型的差异对靠近断层的部分测站模拟存在一定程度的影响。这可能与这两个模型在青川区域断层倾角差异有关。
另外，在断层东北端计算值明显偏大，如点H005、H006、H022等，其中尤以H022站点最为显著。由于文中所选用断层模型长达360 km，为避免因所用同震模型在东北端外延明显更长对结果造成影响，直接剔除同震模型东北端40 km(见图 1红色虚线框)，但仍然发现相同的现象，因此可以排除同震模型东北端相较于实际破裂线更长而直接造成影响。同时对比考虑龙门山断层两侧横向不均衡性的分层模型[14]计算得到的震后形变值发现，H022站点计算值仍然明显大于实测值，这也可能说明巴颜喀拉块体东北端相较于中部本身存在地层差异性。另外，鲜水河断裂带的GPS模拟值相较于观测值整体偏小，如点H345、JB35、SCLH，可以推测，该区域粘滞系数与北部地区不同，也可能说明巴颜喀拉块体一侧西南部与中部横向差异性的存在。
3.3 川西高原流变参数采用Maxwell流体模型估计的粘滞系数可以合理地解释震后累积形变，但不能拟合GPS观测随时间变化的震后形变过程。图 5显示位于川西高原的两个GPS区域站震后时间序列及粘弹性松弛效应模拟值。从模型计算得到的震后时间序列可以看到，震后早期计算值变化趋势明显小于实测值变化趋势，而2015年以后计算值变化趋势又趋向于大于观测值变化趋势。因而，下一步需要采用Burgers体同时考虑震后早期的快速瞬态形变及后期的线性变化。尽管如此，本文约束的粘滞系数对震后库仑应力变化仍有重要参考意义。
|
图 5 GPS站点时间序列拟合值和粘弹性松弛模拟值 Fig. 5 Time series fitted and times series modeled by viscoelastic relaxation |
基于汶川地震震后7 a的GPS观测值，用Maxwell体流变模型对川西高原地壳粘滞系数进行反演，获得川西高原低速层粘滞系数为2.51×1018 Pa·s，而中下地壳粘弹性参数为3.98×1018 Pa·s。虽然Maxwell流体模型不足以完全反映汶川震后长期形变，但仍然可以看到粘弹性松弛对于汶川震后的长期影响，而且在震后10 a仍在持续。后续可以结合震后余滑和Burgers体流变模型对汶川震后形变的细节进行进一步的约束。
| [1] |
Xu X, Wen X, Yu G, et al. Coseismic Reverse- and Oblique-Slip Surface Faulting Generated by the 2008 MW7.9 Wenchuan Earthquake, China[J]. Geology, 2009, 37(6): 515-518 DOI:10.1130/G25462A.1
(  0) |
| [2] |
张勇, 冯万鹏, 许力生, 等. 2008年汶川大地震的时空破裂过程[J]. 中国科学:地球科学, 2009, 52(10): 145-154 (Zhang Yong, Feng Wanpeng, Xu Lisheng, et al. Spatio-Temporal Rupture Process of the 2008 Great Wenchuan Earthquake[J]. Science China Earth Science, 2009, 52(2): 145-154)
( 0) |
| [3] |
Shen Z K, Sun J, Zhang P, et al. Slip Maxima at Fault Junctions and Rupturing of Barriers during the 2008 Wenchuan Earthquake[J]. Nature Geoscience, 2009, 2(10): 718-724 DOI:10.1038/ngeo636
( 0) |
| [4] |
Wang Q, Qiao X J, Lan Q G, et al. Rupture of Deep Faults in the 2008 Wenchuan Earthquake and Uplift of the Longmen Shan[J]. Nature Geoscience, 2011, 4(9): 634-640 DOI:10.1038/ngeo1210
( 0) |
| [5] |
Bürgmann R, Dresen G. Rheology of the Lower Crust and Upper Mantle: Evidence from Rock Mechanics, Geodesy, and Field Observations[J]. Annual Review of Earth and Planetary Sciences, 2008, 36(1): 531-567 DOI:10.1146/annurev.earth.36.031207.124326
( 0) |
| [6] |
Shao Z, Wang R, Wu Y, et al. Rapid Afterslip and Short-Term Viscoelastic Relaxation Following the 2008 MW7.9 Wenchuan Earthquake[J]. Earthquake Science, 2011, 24(2): 163-175 DOI:10.1007/s11589-010-0781-z
( 0) |
| [7] |
丁开华, 许才军, 温扬茂. 汶川地震震后形变的GPS反演[J]. 武汉大学学报:信息科学版, 2013, 38(2): 131-135 (Ding Kaihua, Xu Caijun, Wen Yangmao. Postseismic Deformation Associated with the 2008 Wenchuan Earthquake by GPS Data[J]. Geomatics and Information Science of Wuhua University, 2013, 38(2): 131-135)
( 0) |
| [8] |
Huang M H, Bürgmann R. Probing the Lithospheric Rheology across the Eastern Margin of the Tibetan Plateau[J]. Earth and Planetary Science Letters, 2014, 396: 88-96 DOI:10.1016/j.epsl.2014.04.003
( 0) |
| [9] |
Zhao B, Huang Y, Zhang C, et al. Crustal Deformation on the Chinese Mainland during 1998-2014 Based on GPS Data[J]. Geodesy and Geodynamics, 2015, 6(1): 7-15 DOI:10.1016/j.geog.2014.12.006
( 0) |
| [10] |
Wessel P. A General-Purpose Green's Function-Based Interpolator[J]. Computers & Geosciences, 2009, 35(6): 1247-1254
( 0) |
| [11] |
Wang R, Lorenzo-Martín 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) |
| [12] |
Pollitz F F, Wicks C, Thatcher W. Mantle Flow beneath a Continental Strike-Slip Fault: Postseismic Deformation after the 1999 Hector Mine Earthquake[J]. Science, 2001, 293(5536): 1814-1818 DOI:10.1126/science.1061361
( 0) |
| [13] |
董运洪, 罗三明, 韩月萍, 等. 汶川MS8.0地震北川附近同震垂直位移特征[J]. 地震学报, 2012, 34(5): 611-620 (Dong Yunhong, Luo Sanming, Han Yueping, et al. Co- and Post-Seismic Vertical Displacements of Wenchuan MS8.0 Earthquake near Beichuan[J]. Acta Seismologica Sinica, 2012, 34(5): 611-620 DOI:10.3969/j.issn.0253-3782.2012.05.003)
( 0) |
| [14] |
赵斌.利用震后GPS资料探测青藏高原东、南边界岩石圈流变结构[D].武汉: 武汉大学, 2017 (Zhao Bin. Probing the Rheological Structure across the Eastern and Southern Edges of the Tibetan Plateau Using GPS Observed Postseismic Displacements[D]. Wuhan: Wuhan University, 2017) http://cdmd.cnki.com.cn/Article/CDMD-10486-1017067304.htm
( 0) |
2. Wuhan Base of Institute of Crustal Dynamics, CEA, 40 Hongshance Road, Wuhan 430071, China;
3. The Second Monitoring and Application Center, CEA, 316 Xiying Road, Xi'an 710054, China