2017, Vol. 37
2. 中国地震局地壳应力研究所武汉科技创新基地，武汉市洪山侧路40号，430071
野外考察、Iknos卫星影像解译、InSAR和GPS观测数据反演等结果表明, 2001年昆仑山地震地表破裂带由走向N45°~50°E张破裂、走向N60°~75°E张剪切走滑破裂、走向约N100°剪切走滑破裂、右阶斜列阶区鼓包和左阶斜列阶区张性裂陷等基本地表破裂单元组合而成, 是一条整体走向N100°±10°E的纯剪切左旋走滑破裂带地震地表破裂带。西起库水浣湖东岸90.257°E, 东止于青藏公路东94.795°E, 整体长426 km[1]。在长期的GPS观测结果里，昆仑断裂以大约1 cm/a的滑移速率吸收了南侧羌塘块体和北侧柴达木块体的相对运动，对区域构造动力学过程具有重要的影响[2]。另外，印度板块和欧亚板块的持续挤压和青藏高原的东向挤出被认为是昆仑地震发生的主要构造成因[3]。Lesserre等[4]基于InSAR图像反演了昆仑地震同震位错模型，显示昆仑地震同震最大位错滑移量为8 m，并且大部分位错滑动都分布在20 km以上的地壳浅层。
地震发生后会产生各种不同的震后效应, 如震后余滑、震后孔隙弹性回弹、下地壳和上地幔震后粘弹性松弛等。这些不同的震后活动机制表现出不同的特征, 反映了地球深部物质的不同属性[5]。地震震后粘弹性松弛现象为认识岩石圈的流变结构提供了一种可行的方法，断层破裂相当于是对地球介质加载了一种应力。地震发生后，具有粘滞性的下地壳和上地幔的应力状态发生变化，随后以粘弹性松弛的方式释放并作用到上地壳而形成明显的地表形变。反之，基于观测的震后地表形变，可以构建合理的模型反演下地壳和上地幔的流变结构，为研究其他动力学过程提供依据[6]。已有的研究表明，昆仑山地震震后形变可以用震后余滑和粘弹性松弛两种形变机制来解释，而不同地区的粘滞系数一般在1016~1021 Pas之间。Shen等[6]结合两种震后形变机制，以同震破裂区下倾边的余滑和下地壳的粘弹性松弛来模拟昆仑山地震震后的GPS观测形变。Wen等[7]也考虑两种机制共同产生影响，并将地表下划分不同区域，假定震后余滑发生在厚度为15 km的弹性上地壳，粘弹性松弛发生在中-下地壳以及地表 70 km以下的上地幔中。本文在前人研究的基础上，采用震后2~6 a的InSAR数据，基于不同的深度分层结构，细致反演昆仑山地区下地壳和上地幔的粘滞系数，得到该地区岩石圈的流变结构，对青藏高原北部区域的变形方式以及区域构造活动有了更深刻的认识，对研究昆仑断裂的动力学演化有重要意义。
1 InSAR形变场与采样数据本文使用的T133D轨道北部的数据来源于Wen等[7]处理后的结果，源数据包含了5个轨道40多万个点，所有的ASAR数据均为level 0原始格式数据，干涉处理软件为JPL/Caltech的ROI_PAC开源软件，并使用3″分辨率的SRTM DEM来移除地形的影响，之后采用基于能量谱的局部自适应滤波对图像进行干涉滤波来降低干涉相位噪声，最后采用SNAPHU软件解缠得到差分干涉相位图。InSAR形变图中的量级和误差大小采用一维方差-协方差函数来描述，即
从Wen等[7]震后粘弹性松弛模拟的结果来看，地震地表破裂带南侧的点较为稀薄，远场的模拟效果比近场的更加理想，并且采用震后余滑模型对近场的模拟效果更好。刁法启[3]提出昆仑山地震震后余滑主要影响的是近场，震后下地壳和上地幔的粘弹性松弛效应则是远场震后形变的主要驱动机制。考虑到反演的可行性和计算效率，选择数据较为密集的轨道T133D，范围控制约为93°~94°E、36.5°~37.5°N之间, 选择经度分辨率为0.02°，纬度分辨率控制为0.01°，采样点数量为16 064。
图 1是5个相邻的降轨道ENVISAT图像, 448、176、405、133、362对应的黑色矩形框代表着5个相邻的SAR降轨影像的数据范围，其中红色矩形代表本文所采样的数据范围。
|
图 1 昆仑山地震地表破裂带位置背景图(震源机制解取自美国地质调查局(USGS)发布的CMT解) Fig. 1 Location map of surface rupture zone of Kunlun mountains earthquake(The focal mechanism is from the CMT catalog by USGS) |
震后粘弹性松弛模型中通常使用的两种粘弹性结构为Maxwell体和Burgers体。对于昆仑山可可西里无人区的地震来说，Maxwell体和Burgers体所引起的空间分布非常接近，所以用Maxwell体来模拟本次震后形变是可行的[8]。
2.2 反演方法为了模拟由弹性上地壳以下粘弹性半空间松弛效应引起的震后形变，采用汪荣江[9]的计算程序包PSGRN/PSCMP对Maxwell体介质本构关系进行计算。该软件已被成功应用到过去的一些研究中，虽然该软件基于地球平面分层模型，忽略了地球的曲率影响，但即使是最大的地震，空间域影响范围最大也只扩张到了3 000 km。而这种水平长度下弯曲的高度差只占了不到3%，极小的几何偏差对于形变场的影响是非常有限的，而且相对于Pollitz球面分层模型来说，最后计算得到的形变差别也不会超过几个百分点[9]。
该软件基于半解析平面分层模型，考虑重力影响模拟分层介质中地震的同震及震后形变。软件由两部分组成：第一个程序包PSGRN用来计算给定分层粘弹性半空间中基本位错源时变格林函数；第二个程序包PSCMP以这些格林函数作为基础数据库，自动将地震破裂面离散成许多点位错，再通过线性叠加方法计算同震及震后形变。对于时域格林函数的获取，在快速傅里叶变换过程中引入反混叠技术，以保证计算震后形变的数值稳定性。
根据加权均方差的算法，模拟值与观测值加权中误差为：
|
(1) |
最优模型即使得该中误差最小：
|
(2) |
式中，n为观测值个数，y为观测值，p为观测值的权重，ηc为下地壳粘滞系数，ηm为上地幔粘滞系数，t为粘滞系数的个数，f(ηc, ηm)为震后模拟的形变值，i为观测值编号。权值的大小根据观测值点到震中的距离来定夺，由于远场受震后粘弹性松弛的影响较大，所以距离越远权重越大。下地壳和上地幔的粘滞系数的常用对数以步长0.2在17~21之间变化，因此对于每一种固定了上地壳和总地壳厚度的模型来说，会得到441个计算结果。
根据地质学和地震学的研究结果[10]，分层模型定为弹性的上地壳、粘弹性的下地壳和粘弹性的上地幔，备选的上地壳厚度为10 km、15 km和20 km，地壳总厚度为65 km、70 km、75 km。粘度系数是待求量，基于温度结果和GPS观测得到的应变率数据，青藏高原地区下地壳粘滞系数数量级在1019~1021 Pas之间，在模拟计算的时候，需要把范围扩大一定量级，以得出更可靠的结果，所以对下地壳和上地幔设定的格网搜索范围为1017~1021 Pas。地震破裂断层参数以及地震波波速和介质密度均取自于USGS的结果(表 1)，由GPH于2015-03-18制作发布，断层总长定为432 km，宽度36.64 km，总共分为36×8个子断层。参数设定好后，最后计算得出指定时间和指定位置的震后粘弹性松弛形变。
|
表 1 昆仑山地区地表下波速和密度结构 Tab. 1 Layered crustal structure of the Kunlunshan region |
本文对每一种固定地壳厚度分层模型的441个结果进行格网搜索，最终得到加权中误差最小的分层模型为：上地壳厚度20 km，地壳总厚度为70 km。对应的下地壳粘滞系数为2.5×1019 Pas，上地幔粘滞系数在(1~6.3)×1019 Pas之间，而上地壳厚度为15 km的模型中误差与最小中误差之差仅为0.01，因此考虑最优模型的上地壳厚度在15~20 km之间。图 2给出中误差最小的分层模型，其中左图是地震波波速，右图是介质密度，图中的两条紫色的线表示上地壳与下地壳以及下地壳和上地幔的分界线。图 3给出了该模型下不同粘滞系数下加权中误差分布图，图中黑色点代表的是格网搜索中的最优下地壳和上地幔的粘滞系数组合的其中一个，对应的值均为2.5×1019 Pas，对应的加权中误差的值为0.90 cm，这个结果被应用于图 4的视线向形变模拟中。从图 3可以看出，模型值与观测值的拟合程度对上地幔粘度不是特别敏感，拟合的加权中误差先是快速增大后来快速减小，当下地壳的粘滞系数达到1×1019 Pas的时候，中误差的变化变得逐渐平缓；到了1×1020 Pas左右，加权中误差趋于平稳而不再变化。
|
图 2 最优分层介质模型参数 Fig. 2 Plot of parameters of optimal layered medium model |
|
图 3 在图 2分层模型下不同的下地壳和上地幔粘度对应的加权中误差的分布 Fig. 3 Plot of weighted RMS for different viscosity of lower crust and upper mantle under layeled model of Fig. 2 |
|
图 4 观测值、模拟值的视线向形变以及两者的残差 Fig. 4 The observed, simulated line-of-sight displacement and residuals between them |
在地震破裂断层参数一定的条件下，最佳的粘滞系数受不同的地壳分层厚度影响，模拟计算最终得到了不同的结果(表 2)。
|
|
表 2 下地壳和上地幔粘滞系数及加权中误差 Tab. 2 Viscosities and weighted RMS of lower crust and upper mantle |
图 4给出采用最优模型投影到InSAR的LOS方向的模拟结果，该模型上地壳厚度20 km，地壳厚度为70 km，下地壳与上地幔的粘滞系数均取自图 3中的最优值, 加权中误差为0.90 cm。从残差图中可以看出，模拟的视线向形变量基本小于实际观测的InSAR形变场。可能的原因之一是，即便过了数年，震后的余滑效应对远场靠近震中位置的地表位移仍有一定影响。同时，震后粘弹性松弛的影响范围也是有限度的，远离震区的地表受到其他地区的地震以及本身板块运动和挤压的影响更大。另一个可能原因是，形变场在柴达木盆地内，该地区常年风力强盛，风力侵蚀和堆积现象严重，气候的不稳定也对该地区的地表结构产生了一定的影响。
3 讨论昆仑山地震发生在可可西里无人区，该地区海拔高、环境比较恶劣，这些外界因素给野外GPS观测造成很大的困难，所以该地区的GPS震后观测资料相对较少，而InSAR卫星观测弥补了这个缺陷。本文使用震后2~6 a的InSAR观测数据，获得了丰富的地表水平和垂直形变。虽然数据源相同，但由于重采样的点不同以及同震滑动分布模型的差异，导致粘滞系数的反演结果与文献[7]有差异。
Ryder等[11]根据1997年位于昆仑山断层最西边的西南方向矩震级7.6的玛尼地震，利用震后25 d的InSAR数据，采用标准线性流变体得出下地壳瞬时粘滞系数为4×1018 Pas，而同时依据震后同一InSAR影像的震后1 145 d数据，反演得到下地壳粘滞系数数量级为1019 Pas左右。刁法启[3]以2001年昆仑山地震震后4个月的GPS地表观测位移为约束，在粘弹性松弛的模拟中，得到下地壳的粘滞系数为约1017 Pas。Ryder等[8]采用昆仑山震后2~5 a的InSAR视线向观测以及震后第1年的GPS时间序列得出青藏高原东北部下地壳瞬时粘度为9×1017 Pas，稳态粘度为1×1019 Pas。所以说，震后初期下地壳的粘滞系数偏小，震后变形场的地球深部物质的粘度系数的反演与所选区域和震后形变场的观测时间段相关，岩石圈流变性是随时间变化的, 构造区域和震后形变场的观测时段的不同必然导致所反演的粘度系数具有差异性。随着时间的推移，震区的物质构造逐渐稳定，地壳地幔的粘性系数越来越大并趋向于稳定，反映到震后形变场表现为应变率缩减。
震后形变反演出的粘滞系数需要与外部数据相互验证，如冰后回弹，而根据冰后回弹得到的岩石圈粘滞系数较大，一般为1020~1021 Pas[12]。孙玉军等[13]以滑动摩擦、脆性破裂和蠕变等3种强度机制为约束，基于地震波波速得到的上地幔温度、气象台记录的地表温度结果和GPS观测得到的应变率数据，计算得到了中国大陆和邻区岩石圈三维流变结构，得到青藏高原中地壳厚度约为20 km，下地壳厚度为20~30 km，即韧性层厚度约在40~50 km之间，等效粘滞系数为1019~1021 Pas，与本文结果对比，验证了不同的观测弛豫时间的计算结果具有差异性。
4 结语强震释放能量巨大，引起的应力应变随着震后时间的延长可以传播很远的距离，这是根据震后形变研究震区流变结构的基础。本文根据昆仑山地震震后2~6 a的InSAR形变场，对2001年昆仑山地震的震后形变采用粘弹性松弛模型进行模拟，使用了多个分层结构模型。结果表明，昆仑山地区弹性上地壳的厚度为15~20 km，莫霍面的深度为70 km，下地壳的粘滞系数为(1~2.5)×1019 Pas，上地幔的粘滞系数在(1~6.3)×1019 Pas之间，模拟的形变较好地拟合了该区域的地表形变。该结果与青藏高原玛尼震区的震后粘弹性松弛反演结果基本一致。
| [1] |
徐锡伟, 于贵华, 马文涛, 等. 昆仑山地震(MW7.8)破裂行为、变形局部化特征及其构造内涵讨论[J]. 中国科学:地球科学, 2008, 38(7): 785-796 (Xu Xiwei, Yu Guihua, Ma Wentao, et al. Discussion on Rupture Behavior and Localization of Deformation of Kunlunshan Earthquake Mw7.8 and Its Tectonic Implication[J]. Science in China:Earth Sciences, 2008, 38(7): 785-796)
(  0) |
| [2] |
Wang Q, Zhang P Z, Freymueller J T, et al. Present-Day Crustal Deformation in China Constrained by Global Positioning System Measurements[J]. Science, 2001, 294(5542): 574-577 DOI:10.1126/science.1063647
( 0) |
| [3] |
刁法启.基于GPS观测的同震, 震后形变研究[D].北京: 中国科学院大学, 2011 (Diao Faqi.Co-Seismic and Post-Seismic Deformation Studies Based on GPS Observations[D].Beijing: University of Chinese Academy of Sciences, 2011) http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y2035439
( 0) |
| [4] |
Lasserre C, Peltzer G, Crampé F, et al. Coseismic Deformation of the 2001 MW7.8 Kokoxili Earthquake in Tibet, Measured by Synthetic Aperture Radar Interferometry[J]. Journal of Geophysical Research: Solid Earth, 2005, 110(B12): 1-18
( 0) |
| [5] |
谭凯, 王琪, 王晓强, 等. 震后形变的解析模型和时空分布特征[J]. 大地测量与地球动力学, 2005, 25(4): 23-26 (Tan Kai, Wang Qi, Wang Xiaoqiang, et al. Analytic Models and Space Time Distribution of Postseismic Deformation[J]. Journal of Geodesy and Geodynamics, 2005, 25(4): 23-26)
( 0) |
| [6] |
Shen Z K, Zeng Y, Wang M, et al. Postseismic Deformation Modeling of the 2001 Kokoxili Earthquake, Western China[C].EGS-AGU-EUG Joint Assembly, Nice, 2003 https://www.researchgate.net/publication/253764179_Postseismic_deformation_modeling_of_the_2001_Kokoxili_earthquake_western_China
( 0) |
| [7] |
Wen Y M, Li Z H, Xu C J, et al. Postseismic Motion after the 2001 MW7.8 Kokoxili Earthquake in Tibet Observed by InSAR Time Series[J]. Journal of Geophysical Research: Solid Earth, 2012, 117(B8)
( 0) |
| [8] |
Ryder I, Bürgmann R, Pollitz F. Lower Crustal Relaxation Beneath the Tibetan Plateau and Qaidam Basin Following the 2001 Kokoxili Earthquake[J]. Geophysical Journal International, 2011, 187(2): 613-630 DOI:10.1111/gji.2011.187.issue-2
( 0) |
| [9] |
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) |
| [10] |
许志琴, 杨经绥, 姜枚, 等. 青藏高原北部东昆仑-羌塘地区的岩石圈结构及岩石圈剪切断层[J]. 中国科学:地球科学, 2001, 31(B12): 1-7 (Xu Zhiqin, Yang Jingsui, Jiang Mei, et al. Lithospheric Structure and Shear Fault from Eastern Kunlun to Qiangtang in Tibet Plateau[J]. Sci China Earth Sci, 2001, 31(B12): 1-7)
( 0) |
| [11] |
Ryder I, Parsons B, Wright T J, et al. Post-Seismic Motion Following the 1997 Manyi (Tibet) Earthquake: InSAR Observations and Modelling[J]. Geophysical Journal International, 2007, 169(3): 1009-1027 DOI:10.1111/gji.2007.169.issue-3
( 0) |
| [12] |
Nakiboglu S M, Lambeck K. A Reevaluation of the Isostatic Rebound of Lake Bonneville[J]. Journal of Geophysical Research: Solid Earth, 1983, 88(B12): 10439-10447 DOI:10.1029/JB088iB12p10439
( 0) |
| [13] |
孙玉军, 董树文, 范桃园, 等. 中国大陆及邻区岩石圈三维流变结构[J]. 地球物理学报, 2013, 56(9): 2936-2946 (Sun Yujun, Dong Shuwen, Fan Taoyuan, et al. 3D Rheological Structure of the Continental Lithosphere Beneath China and Adjacent Regions[J]. Chinese Journal of Geophysics, 2013, 56(9): 2936-2946)
( 0) |
2. Wuhan Base of Institute of Crustal Dynamics, CEA, 40 Hongshance Road, Wuhan 430071, China