2022, Vol. 42
2. 防灾科技学院, 河北省三河市学院街465号, 065201
2021-03-30西藏双湖县黄水湖西岸发生MW5.7地震，中国地震台网中心(CENC)给出震中位置为34.38°N、87.68°E，震源深度为10 km。双湖县位于青藏高原腹地的羌塘块体中部，自1970年有地震记录以来，羌塘块体中部共发生过14次5级以上地震，本次地震是周边100 km范围内的首个5级以上地震(图 1(a))。由于自然环境恶劣、交通限制等原因，羌塘块体中部测震台和GNSS连续站非常稀疏，地质学资料有限，距离震中200 km范围内除震中西侧约1 km处有一条近SN向、倾角不明的无名正断层外(本文称黄水湖断层)，尚未有其他已知断层[1-2]。受限于震中区域地质构造、余震和同震GPS等资料较少，发震断层的确定较难，CENC未给出此次地震的震源机制解，而GCMT和USGS给出的震源机制解差异较大，发震断层面选择不统一。
|
断裂数据来源于文献[1-2];震源机制解源于哈佛大学GC-MT;历史地震源于CENC;GPS速度场源于文献[3] 图 1 羌塘块体中部历史地震及断裂分布 Fig. 1 The distribution of the faults and historical earthquakes of the central Qiangtang block |
InSAR可大范围、高精度地获取地震同震形变场，是确定中强地震发震构造的重要手段之一。本文利用Sentinel-1A SAR-C地震前后升轨和降轨4景数据获取2021年双湖地震同震形变场，基于半空间位错模型确定发震构造几何学参数，计算不同深度上的同震静态库仑应力变化，并结合羌塘块体构造背景讨论本次地震的发震机制及其对区域未来地震危险性的影响，为深入理解羌塘块体中部构造活动性和地壳变形机制提供参考依据。
1 区域构造背景GPS速度场(图 1(b))显示，青藏高原内部各次级块体呈现不同的地壳形变机制[4-5]。班公湖-怒江缝合带以北的羌塘块体从西向东滑移速率逐渐减小，羌塘块体中部至东部滑移速率相差可达10.3 mm/a[6]，表现出明显的空间差异性；班公湖-怒江缝合带以西在上地壳物质持续向东挤出的构造作用下，沿羌塘内部小型正断层和共轭走滑断裂发生弥散型变形特征，构造活动以EW向延伸扩展为主；班公湖-怒江缝合带以东则直接向东流出，使得块体东边界不断向东扩张，发生刚性块体变形，呈现明显的物质快速转移运动特征[4]。
目前，对羌塘块体的断裂活动性仍缺乏定量研究，特别是地处青藏高原腹地的块体中部，仅Wang等[7]利用时序InSAR技术发现羌塘块体中部存在持续向东的滑移运动，速率约为4 mm/a。2021年双湖发震发生于羌塘块体中部的黄水湖西岸，其发震构造和变形机制不明，变形特征是呈弥散型还是刚性整体向东挤出均有待进一步研究。
2 InSAR数据处理本文利用Sentinel-1A影像(表 1)，基于两轨差分干涉获取本次双湖地震的干涉影像，提取相应的同震形变场。基于开源雷达数据处理平台GMTSAR[8]进行以下数据处理：1)采用20 ∶4的多视比例；2)通过GACOS移除影像中的大气误差[9]；3)根据精密轨道和SRTM提供的30 m DEM去除轨道误差和平地相位；4)采用Goldstein滤波方法和最小费用流解缠方法分别对干涉图进行滤波和相位解缠处理；5)地理编码。
|
表 1 干涉对信息 Tab. 1 Interferometry pairs information |
主震升降轨干涉图和同震形变场如图 2所示。可以看出，形变场以近SN向分界呈三瓣式(L1、L2、L3)分布。L1部分在升降轨形变场中最大视线向(line of sight，LOS)形变量分别为0.53 cm、-0.46 cm，二者符号相反、量值接近；升轨时远离、降轨时靠近L2部分，指示L1存在NWW向水平运动。L2部分在升降轨形变场中皆呈现负值，LOS向形变量分别为-6.78 cm、-4.73 cm，表明该区域整体为垂直沉降运动。L3部分在升降轨形变场中LOS形变量分别为-1.65 cm、1.45 cm，符号相反、量值接近；升轨时接近、降轨时远离L2部分，指示L3存在近东向水平运动。图 2(g)显示，沿AB测线得到的剖面T114A在从L1向L2过渡时(L1与L2之间的虚线)，其梯度明显大于剖面T121D，可能是由于T121D的成像视角与干涉图长轴方向近乎平行，导致真实地表形变在形变场中的投影分量较少，限制了LOS形变的观测精度。
|
暖色代表接近卫星方向或抬升运动,冷色代表远离卫星方向或沉降运动 图 2 双湖地震同震形变场 Fig. 2 Coseismic deformation field of the Shuanghu earthquake |
本次双湖事件形变场长轴为近SN走向，主要以垂直向形变为主。因InSAR对SN方向形变分量不敏感，因此在忽略SN向形变分量的情况下，利用不同视角的升降轨同震形变场提取震区EW水平位移场和垂直向位移场[10](图 2(c)、2(f))更适用于本次地震。计算得出，EW向和垂直向形变场的残差(RMS)分别为0.4 cm和1.3 cm，在InSAR精度范围内，说明2.5D形变场具有一定的可靠性。从图 2(c)水平向形变场来看，L1表现为西向水平位移，最大位移值为2.5 cm；L2、L3则为东向水平位移，最大位移值为5.8 cm，且出现在L2区间；从图 2(f)垂直向形变场来看，最大沉降量达11.8 cm，出现在L2区间；L1、L3部分呈微弱抬升，最大抬升量为2 cm。此外，对比3组形变场来看，在L1与L2的过渡区域，与T121D形变场相比，T114A表现出更为明显的近SN向长轴走向。由此推断，造成T114A和T121D沿测线形变量变化梯度差距较大(图 2(g))的可能原因是T121D成像角度与形变场近SN向长轴几乎垂直，对形变场水平向的位移并不敏感，且L1部分本身水平位移量较小。
3 断层滑动分布反演为精细地确定发震断层位错滑移特征，采用以下步骤进行计算：1)以同震形变场为输入数据，利用蒙特卡洛算法搜索最佳断层参数，获得发震断层的均匀滑动模型；2)以最佳均匀滑动模型参数为先验条件约束，利用Okada位错模型反演同震滑动分布。
3.1 均匀滑动反演为提高最佳断层参数的搜索效率，确保有效收敛，使用四叉树采样算法[11]对升降轨同震形变场进行降采样处理。处理后得到494个升轨、422个降轨形变点位，以此为蒙特卡洛搜索约束[12]，获得收敛性较好的马尔科夫链(图 3)，并统计特征最佳的均匀滑动模型断层参数(表 2)。
|
红线表示最佳断层模型参数 图 3 蒙特卡洛-马尔科夫链联合概率密度分布 Fig. 3 Distribution of joint probability density of Monte Carlo-Markov chain |
|
|
表 2 均匀滑动模型 Tab. 2 Uniform slip model |
表 2中最佳均匀滑动模型的倾角为51°，走向为36.6°，走向与主震干涉图NNE向长轴走向一致；断层走滑量为0.39 m，倾滑量为-0.61 m，指示双湖地震可能是一次正断为主、兼右旋走滑的地震，但已有研究并未表明该区域存在右旋走滑断层[13]。因此，有必要构建滑动分布模型，进一步正演精细的三维形变特征以确定发震断层几何学参数。
3.2 断层滑动分布反演以均匀滑动模型得到的最佳断层模型参数为先验条件，基于Okada位错模型[14]进行同震断层滑动分布反演[15]，使用Crust1.0模型，反演约束方程为：
| $ \left( {\begin{array}{*{20}{c}} \mathit{\boldsymbol{G}}\\ {{\gamma ^2}{\nabla ^2}} \end{array}} \right)m = \left( {\begin{array}{*{20}{l}} \mathit{\boldsymbol{d}}\\ 0 \end{array}} \right) $ | (1) |
式中，G为联系观测值与弹性位错模型的格林函数矩阵，d 为观测值，▽2为拉普拉斯算子，γ2为平滑因子。
根据均匀滑动模型，设定以下初始参数：断层长18 km、宽20 km、倾角50°、走向30°、滑动角范围-130°~-50°。将断层按1 km×1 km划分为360个子断层，并且设定平滑因子γ2为0.07。最终得到的断层滑动模型为：震中34.37°N、87.71°E，震源深度6.51 km，走向33.0°、倾角50°、平均滑动角-74°、最大滑动量0.26 m(图 4)。从图 4(a)~4(c)可见，破裂集中在沿走向8~16 km(滑动量大于10 cm)、地下深度5~15 km处，断层滑动以倾滑为主，兼有微弱左旋分量。假定剪切模量为30 GPa，计算得到地震矩张量为4.5×1017Nm，相当于震级MW5.74。
|
图(c)中,视线方位角为150°,高度角为20°,加粗实线为断层上盘;图(d)中,箭头表示水平运动方向 图 4 断层滑动分布模型及模拟三维形变场 Fig. 4 Fault slip distribution model and the simulated 3D deformation field |
图 4(d)是基于滑动分布模型重建的矢量形变场。水平向位移显示，断层西侧向西移动，东侧则微弱东移，呈左旋走滑特征；垂直向形变场存在明显的NNE向分界线，以沉降运动为主。显然，模拟的矢量形变场整体形变特征与2.5D形变场基本一致，且更加精细地展示了断层在不同维度上的滑动量值。
图 5是基于均匀滑动模型、滑动分布模型得到的正演模型及其拟合残差。可以看到，T114A的2个模型与观测值拟合度更好，RMS均为0.3 cm，残差主要出现在形变场沉降区，主要是因黄水湖干扰而缺乏部分近场观测信号导致的，滑动分布模型以清晰的NNE向分界线指示了发震构造行迹。T121D的均匀滑动模型和滑动分布模型拟合的RMS分别为0.3 cm和0.6 cm，但均匀滑动模型的指示意义相对不佳。联合T114A得到的断层走向参数约束后，可以改善T121D滑动分布模型在断层走向的拟合效果(图 5(g)、5(i))，进一步表明T121D近乎平行于断层走向的飞行方向导致该方向上的形变获取受限，进而使得在该方向上的形变量低于T114A。
|
黑色矩形框表示滑动分布模型断层地表投影,加粗实线表示断层上盘,(a)~(e)和(f)~(j)分别对应T114A和T121D 图 5 同震形变场观测值、模型值与残差 Fig. 5 The observation values, model values and residual of coseismic deformation field |
均匀滑动模型和滑动分布模型获得的震源机制参数见表 3。可以看出，InSAR反演的模型反映的是地震破裂面的总体特征，与USGS、GCMT给出的远震体波模型结果基本一致，仅震源深度相差较大，这是由于远震体波模型地震矩张量反演中波形拟合误差随深度变化，深度结果误差相对较大[16-17]。
|
|
表 3 2021年双湖MW5.7地震震源参数 Tab. 3 Focal parameters of the 2021 Shuanghu MW5.7 earthquake |
测震学和大地测量学结果都显示，2021年双湖MW5.7地震为一次正断为主、兼具走滑分量的事件。但受特定方位的远震地震资料信噪比低等影响，USGS和GCMT给出的节面I的滑动角以及InSAR均匀滑动模型和滑动分布模型的滑动角分别表现出右旋或左旋走滑，干扰了走滑方向的确定。已有地质学资料显示，羌塘块体中部已知断层未见右旋走滑性质[13]。由于震中区域缺少GPS连续观测站，测震台也非常稀疏(图 1)，主震后2个月内仅发生8次余震(表 4，数据来源于CENC地震目录)，无法提供更多资料约束发震节面。
|
|
表 4 2021年双湖MW5.7主震及余震信息 Tab. 4 2021 Shuanghu MW5.7 main shock and aftershocks information |
为了解决上述问题，采用静态库仑应力转移与同震形变场、余震的对应关系来讨论不同性质发震节面的差异，进而判定发震构造性质，评估本次地震对区域未来地震危险性的影响。考虑到InSAR滑动分布模型给定的左旋节面解与GCMT节面I参数更接近，为对比分析，采用Crust1.0模型、摩擦系数0.4、杨氏模量30 GPa、泊松比0.25，分别以InSAR左旋节面解(走向33°、倾角50°，滑动角-74°)和GCMT右旋节面解(走向179°、倾角54°，滑动角-112°)为接收面计算不同深度上的静态库仑应力，结果分别见图 6(a)~6(g)、6(h)~6(n)。
|
(a)~(g)为InSAR结果;(h)~(n)为GCMT结果;绿色、黄色小球分别表示不同震级的余震分布,序号与表4对应 图 6 不同接收面对应的不同深度上同震静态库仑应力变化与主余震分布 Fig. 6 Coseismic static Coulomb stress changes at different depths corresponding to different receivers and distribution of main shock and aftershocks |
从图 6(a)~6(g)看出，5 km深度处，库仑应力呈现明显的对称特征，应力卸载区呈EW向，应力加载区呈SN向，东侧卸载区与南侧加载区分界明显，呈NNE向，与左旋节面走向一致，左旋节面给定的主震即落于该分界处；最大应力增加为0.6 bar，超过应力触发阈值0.1 bar，存在较大地震风险。7 km深度处，应力卸载区主要沿左旋节面走向进一步扩张，与左旋节面给定的主震位置和SAR观测到的主要形变区范围(图 5)基本一致，同时南侧应力加载明显。考虑到InSAR左旋节面给定的主震深度为6.51 km，且在7 km处释放了28.9%的力矩，推断该深度为主震应力释放的主要深度。10 km深度处，SN向应力加载区范围增大，EW向应力卸载区范围减小，且呈东西两瓣沿左旋节面方向反向对称；最大应力增加达0.9 bar，远大于应力触发阈值0.1 bar，发生于2021-04-02、震源深度为9 km的4级余震即落于东侧卸载区与北侧加载区的分界处；南侧应力加载也较明显，仍存在较大地震风险。15 km深度处，SN向应力加载区进一步扩张并连通，主震所在区域应力持续卸载。在20~30 km深度处仍以应力加载为主，最大应力在20 km深度时增加到0.4 bar，而其他7次余震都集中发生于该深度区间，且均落于库仑应力大于0.1 bar的应力加载区或边界。
从图 6(h)~6(n)看出，应力变化分布与左旋节面解结果存在较大差异，主要表现在4个方面: 1) 5~10 km深度，应力卸载区呈NW向扩张，与右旋节面走向不符；2) 最大应力卸载区与SAR观测到的主要形变区范围(图 5)相差较大，且出现在7 km深度处，与表 3中GCMT模型给出的主震深度均不一致；3) GCMT给出的主震深度为16.4 km，但在15 km深度处未出现明显的应力卸载区，相反在主震东侧应力加载明显，而所有余震均未出现在该区域范围；4)除2021-04-02的4级余震震源深度为10 km左右，其余7次余震震源深度都集中在25~32 km，但在25 km、30 km深度处余震所在位置绝大部分对应应力卸载区。
综上分析，由左旋节面作为接收面获得的静态库仑应力变化与SAR观测得到的形变场、主余震分布更加契合，也符合羌塘块体区域已有应力场及地质学研究结果[6, 13]。可以认为，2021年双湖MW5.7地震的发震断层应为正断兼左旋走滑性质。从图 5、图 6(a)~6(g)看出，主震和2021-04-02的4级余震发生于10 km深度以上的浅部地壳，且分布于黄水湖断层南段；其他余震均发生于25 km以下深度，且大都散落于黄水湖断层两侧，集中在断层面地表投影范围内。分析图 6中不同深度上库仑应力变化及余震分布，认为黄水湖断层南段短时间内地震风险性较小，北段地震风险性需结合更多资料进一步分析。
本次地震发生于羌塘块体中部的NNE向地堑区，InSAR形变场特征反映这是一次典型的地堑扩张活动。从区域状态来看[6]，羌塘块体整体张、压、剪应变都较强烈，为扩张型地块。从西至东，块体SN向扩张率逐渐增加，最大达24.40×10－9/a；而EW向扩张率逐渐减小，低至3.21×10－9/a；块体中部的主压、主张、最大剪切和面应变率分别为14.29×10－9/a、25.13×10－9/a、35.41×10－9/a和6.84×10－9/a，表明羌塘块体持续向东延展。Yu等[18]指出，羌塘块体东边界可能已向东延伸至95°E，而不是Taylor等[2]认为的93°E。因此，本次地震进一步指示，羌塘块体仍持续向东扩张，因EW向拉伸作用导致近SN向小规模正断层呈弥散型变形和地堑发育。
5 结语1) 本次双湖地震震中位置为34.37°N、87.71°E，震源深度6.51 km，发震断层走向33°、倾角50°、倾向东、滑动角-74°，以倾滑为主兼少量左旋走滑分量，最大滑动量0.26 m；
2) 短时间内，黄水湖断层南段地震风险性较小，北段则需结合更多资料进一步分析；
3) 本次地震是在羌塘块体持续向东扩张的背景下，受EW向拉伸作用使得黄水湖正断层发生的一次弥散型变形活动，SN向地堑得到进一步扩张。
| [1] |
Tapponnier P, Xu Z Q, Roger F, et al. Oblique Stepwise Rise and Growth of the Tibet Plateau[J]. Science, 2001, 294(5 547): 1 671-1 677
(  0) |
| [2] |
Taylor M, Yin A. Active Structures of the Himalayan-Tibetan Orogen and Their Relationships to Earthquake Distribution, Contemporary Strain Field, and Cenozoic Volcanism[J]. Geosphere, 2009, 5(3): 199-214 DOI:10.1130/GES00217.1
( 0) |
| [3] |
Wang M, Shen Z K. Present-Day Crustal Deformation of Continental China Derived from GPS and Its Tectonic Implications[J]. Journal of Geophysical Research: Solid Earth, 2020, 125(2)
( 0) |
| [4] |
李海兵, 潘家伟, 孙知明, 等. 大陆构造变形与地震活动——以青藏高原为例[J]. 地质学报, 2021, 95(1): 194-213 (Li Haibing, Pan Jiawei, Sun Zhiming, et al. Continental Tectonic Deformation and Seismic Activity: A Case Study from the Tibetan Plateau[J]. Acta Geologica Sinica, 2021, 95(1): 194-213)
( 0) |
| [5] |
张培震, 邓起东, 张竹琪, 等. 中国大陆的活动断裂、地震灾害及其动力过程[J]. 中国科学: 地球科学, 2013, 43(10): 1 607-1 620 (Zhang Peizhen, Deng Qidong, Zhang Zhuqi, et al. Active Faults, Earthquake Disasters and Their Dynamic Processes in Mainland China[J]. Science China: Earth Sciences, 2013, 43(10): 1 607-1 620)
( 0) |
| [6] |
李延兴, 杨国华, 李智, 等. 中国大陆活动地块的运动与应变状态[J]. 中国科学: 地球科学, 2003, 33(增1): 65-81 (Li Yanxing, Yang Guohua, Li Zhi, et al. Movement and Strain State of Active Blocks in Mainland China[J]. Science China: Earth Sciences, 2003, 33(S1): 65-81)
( 0) |
| [7] |
Wang H, Wright T J, Jing L Z, et al. Strain Rate Distribution in South-Central Tibet from Two Decades of InSAR and GPS[J]. Geophysical Research Letters, 2019, 46(10): 5 170-5 179 DOI:10.1029/2019GL081916
( 0) |
| [8] |
Sandwell D, Mellors R, Tong X P, et al. Open Radar Interferometry Software for Mapping Surface Deformation[J]. Eos, Transactions American Geophysical Union, 2011, 92(28): 234-235 DOI:10.1029/2011EO280002
( 0) |
| [9] |
Yu C, Li Z H, Penna N T, et al. Generic Atmospheric Correction Model for Interferometric Synthetic Aperture Radar Observations[J]. Journal of Geophysical Research: Solid Earth, 2018, 123(10): 9 202-9 222 DOI:10.1029/2017JB015305
( 0) |
| [10] |
He P, Wen Y M, Ding K H, et al. Normal Faulting in the 2020 MW6.2 Yutian Event: Implications for Ongoing E-W Thinning in Northern Tibet[J]. Remote Sensing, 2020, 12(18)
( 0) |
| [11] |
Jónsson S, Zebker H, Segall P, et al. Fault Slip Distribution of the 1999 MW7.1 Hector Mine, California, Earthquake, Estimated from Satellite Radar and GPS Measurements[J]. Bulletin of the Seismological Society of America, 2002, 92(4): 1 377-1 389 DOI:10.1785/0120000922
( 0) |
| [12] |
Bagnardi M, Hooper A. Inversion of Surface Deformation Data for Rapid Estimates of Source Parameters and Uncertainties: A Bayesian Approach[J]. Geochemistry, Geophysics, Geosystems, 2018, 19(7): 2 194-2 211 DOI:10.1029/2018GC007585
( 0) |
| [13] |
Taylor M, Yin A, Ryerson F J, et al. Conjugate Strike-Slip Faulting along the Bangong-Nujiang Suture Zone Accommodates Coeval East-West Extension and North-South Shortening in the Interior of the Tibetan Plateau[J]. Tectonics, 2003, 22(4)
( 0) |
| [14] |
Okada Y. Surface Deformation Due to Shear and Tensile Faults in a Half-Space[J]. Bulletin of the Seismological Society of America, 1985, 75(4): 1 135-1 154 DOI:10.1785/BSSA0750041135
( 0) |
| [15] |
Wang R, Schurr B, Milkereit C, et al. An Improved Automatic Scheme for Empirical Baseline Correction of Digital Strong-Motion Records[J]. Bulletin of the Seismological Society of America, 2011, 101(5): 2 029-2 044 DOI:10.1785/0120110039
( 0) |
| [16] |
Dziewonski A M, Woodhouse J H. An Experiment in Systematic Study of Global Seismicity: Centroid-Moment Tensor Solutions for 201 Moderate and Large Earthquakes of 1981[J]. Journal of Geophysical Research: Solid Earth, 1983, 88(B4): 3 247-3 271 DOI:10.1029/JB088iB04p03247
( 0) |
| [17] |
Wright T J, Parsons B E, Jackson J A, et al. Source Parameters of the 1 October 1995 Dinar(Turkey) Earthquake from SAR Interferometry and Seismic Bodywave Modelling[J]. Earth and Planetary Science Letters, 1999, 172(1-2): 23-37 DOI:10.1016/S0012-821X(99)00186-7
( 0) |
| [18] |
Yu J S, Zhao B, Xu W B, et al. Oblique Fault Movement during the 2016 MW5.9 Zaduo Earthquake: Insights into Regional Tectonics of the Qiangtang Block, Tibetan Plateau[J]. Journal of Seismology, 2020, 24(3): 693-708 DOI:10.1007/s10950-020-09930-7
( 0) |
2. Institute of Disaster Prevention, 465 Xueyuan Street, Sanhe 065201, China