2022, Vol. 42
2. 中国地震局地震大地测量重点实验室，武汉市洪山侧路40号，430071;
3. 湖北省地震局，武汉市洪山侧路48号，430071;
4. 河北省地震动力学重点实验室，河北省三河市学院街465号，065201
据中国地震台网中心(CENC)测定，2021-03-30西藏那曲市双湖县(87.68°E, 34.38°N)发生MS5.8地震，震源深度10 km。美国地质调查局(USGS)提供的体波矩张量解显示，该地震震级为MW5.56，震源深度4 km。该地震震中位于北羌塘盆地，邻近琵琶湖-吐坡错断裂与琵琶湖-映天湖断裂交汇处[1]。该地区1976年以来共发生38次MW≥5地震(https://www.globalcmt.org/CMTsearch.html)，其中大部分地震表现为正断性质，近年发生的最大地震为2020-07-23西藏尼玛MS6.6地震(距离本次地震约155 km)(图 1)。该地震震中附近5 km范围内平均海拔约5 023 m，交通不便，缺少GNSS等观测数据及地质资料。本文对Sentinel-1A数据进行DInSAR处理，获取同震形变场，通过基于序列蒙特卡罗采样的贝叶斯方法反演发震断层的几何参数[2]和同震滑动分布，可为快速、准确认识双湖MS5.8地震的震中位置及同震特征提供参考。
|
图 1 2021年双湖MS5.8地震构造背景 Fig. 1 Tectonic setting of the 2021 Shuanghu MS5.8 earthquake |
本文采用欧空局官方网站(https://sentinel.esa.int/)提供的Sentinel-1A干涉宽幅数据，具体信息见表 1。
|
表 1 Sentinel-1A数据相关参数 Tab. 1 Related parameters of sentinel-1A data |
基于ISCE(InSAR scientific computing environment)软件对Sentinel-1A数据进行DInSAR处理[3]。在具体处理流程中，距离向和方位向按5 ∶1进行多视处理，生成干涉图；利用精密轨道文件(https://scihub.copernicus.eu/gnss/#/home)和SRTM(shuttle radar topography mission)数字高程模型[4]，分别消除轨道和地形影响；为提高信噪比，采用Goldstein方法进行滤波[5]；相位解缠采用软件SNAPHU(statistical-cost, network-flow phase-unwrapping algorithm)[6]；通过地理编码，获取WGS84坐标系下双湖地震的LOS同震形变场，其中对相关性低值区进行掩膜处理。从双湖地震的同震干涉图(图 2)可以看出，形变长轴呈近NS向展布，影响范围约10 km×10 km。
|
图 2 双湖地震同震干涉图 Fig. 2 Coseismic interferograms of Shuanghu earthquake |
为获取发震断层相关的几何形状参数，利用基于序列蒙特卡罗采样的贝叶斯方法进行反演[2, 7-8]。本文采用四叉树算法对LOS同震形变场进行降采样处理[9]，最终得到426个数据点，其中T114轨道228个，T121轨道198个，计算过程共迭代1.9×106次。图 3为断层几何参数1D和2D的后验概率密度函数(posterior probability density function, PDF)，可表征几何参数的离散程度。选取最大后验解(maximum a posterior solutions, MAP)作为最佳断层参数：东偏移量-0.452 km、北偏移量7.025 km，表示的是参考GCMT震中(87.7° E, 34.32° N)的偏移量，代表断层上边界中心点所处的位置，经坐标转换后为87.695° E、34.383° N，相应深度为3.662 km，断层长度、宽度分别为4.271 km和6.461 km，断层走向、倾角、滑动角分别为47.810°±2.218°、36.664°±2.499、-66.598°±3.258°。具体参数见表 2。
|
图 3 双湖地震矩形震源的1D和2D后验概率密度分布 Fig. 3 1D and 2D PDF of a rectangular source for Shuanghu earthquake |
|
|
表 2 反演的断层几何形状参数 Tab. 2 Inversion parameters of fault geometry |
基于上述断层模型得到的LOS残差分布如图 4所示，图中浅灰色散点表示采样点，黑色直线与矩形分别表示断层上边界和断层面的位置。从图中可以看出，升降轨均显示以沉降为主，最大沉降量分别约5.8 cm和4.8 cm。对于T114升轨，残差范围为-0.9~0.9 cm，残差均方根(RMS)为0.4 cm；对于T121降轨，残差范围为-0.6~0.7 cm，残差RMS为0.3 cm，表明反演获取的断层几何参数较为合理。
|
图 4 断层几何参数反演获得的LOS形变残差 Fig. 4 LOS deformation residuals of fault geometry inversion |
在已知断层几何参数的基础上，将断层离散为若干子断层，断层滑动分布反演可简化为求解每个子断层对应的滑动。利用基于Python语言的Pyrocko-GF软件生成格林函数库，然后通过基于序列蒙特卡罗采样的贝叶斯方法反演滑动分布[2, 10]，子断层滑动分布与相应格林函数库的乘积即为地表位移。每个子断层的滑动量搜索范围设置为0~1.0 m。
在贝叶斯方法中，正则化可通过高斯先验p(s, α)实现，其包含一个非对角项的协方差矩阵，Cholesky分解等价于一个大小为M×M的拉普拉斯有限差分算子L(M为子断层个数)。p(s, α)可表示为[2]：
| $ p(\mathit{\boldsymbol{s}}, \alpha ) = \frac{{{{\left| {{\mathit{\boldsymbol{L}}^{\rm{T}}}} \right|}^{1/2}}}}{{{{\left( {2\pi {\alpha ^2}} \right)}^{M/2}}}} \times \exp \left[ { - \frac{1}{{2{\alpha ^2}}}{{(\mathit{\boldsymbol{Ls}})}^{\rm{T}}}(\mathit{\boldsymbol{Ls}})} \right] $ | (1) |
正则化后的后验概率密度函数PDF可表示为[11]：
| $ \begin{array}{*{20}{c}} {p\left( {\mathit{\boldsymbol{m}}, {\sigma _1}, {\sigma _2}, \cdots , {\sigma _K}\mid {\mathit{\boldsymbol{d}}_{{\rm{obs}}, 1}}, {\mathit{\boldsymbol{d}}_{{\rm{obs}}, 2}}, \cdots , {\mathit{\boldsymbol{d}}_{{\rm{obs}}, K}}} \right)}\\ { \propto p(\mathit{\boldsymbol{m}}) \times p(s\mid \alpha )\prod\limits_{k = 1}^K p \left( {{\mathit{\boldsymbol{d}}_{{\rm{obs}}, k}}\mid \mathit{\boldsymbol{m}}, {\sigma _k}} \right)} \end{array} $ | (2) |
式中，m为已知断层参数，σK为第K个超参数，dobs, K为第K个观测数据，平滑因子α可作为待求参数，s为子断层待求的滑动矢量。
在反演断层滑动分布时，参考表 2结果对断层参数进行初始化设置，其中走向取47.810°，倾角取36.664°，延伸断层长、宽分别为15 km和19 km，各子断层大小为1 km×1 km，平均滑动角固定为-66.598°。
反演采用的数据共426个，与前文反演断层几何参数所使用的数据个数一致。从同震滑动分布(图 5)可以看出，主要滑动区集中在深度3.9~7.5 km，最大滑动量约0.42 m，平均滑动角结果表明，双湖地震同震破裂主要表现为正断拉张兼少量左旋走滑性质。反演得到的震中位置为87.715° E、34.365° N，震源深度为5.763 km。假设地壳的剪切模量为30 GPa，计算得到地震释放的地震矩约为2.189×1017 Nm，矩震级为MW5.5，略小于USGS提供的震级。
|
图 5 基于InSAR数据反演的同震滑动分布 Fig. 5 Coseismic slip distribution constrained by InSAR |
图 6为同震滑动分布模型的残差分布，图中小矩形表示断层位置，大矩形表示拓展长宽后得到的断层面。从图中可以看出，T114升轨的残差范围为-1.2~1.5 cm，T121降轨的残差范围为-0.5~0.9 cm，升降轨残差RMS分别为0.7 cm和0.3 cm，表明滑动分布结果具有合理性与可靠性。
|
图 6 断层滑动分布模型的LOS形变残差 Fig. 6 The LOS deformation residuals of the fault slip distribution model |
表 3为不同机构提供的震源参数，由于采用的数据和方法不同，获得的参数存在差异。反演结果表明，发震断层属于具有正断性质的隐伏断层，震中位于琵琶湖-吐坡错断裂与琵琶湖-映天湖断裂的交汇处，地震可能受这2个断裂共同作用。Sentinel-1卫星对南北方向不敏感，当南北向分量近似忽略时，通过升降轨数据可以获取近东西向和垂直向的形变分量，称为2.5D形变场[12]，结果如图 7所示，图中灰色矩形框表示滑动分布的断层面，黑色和红色五角星分别表示CENC和InSAR计算的震中位置，后者更靠近形变中心。从图 7可以看出，东向运动最大值约2.8 cm，西向运动最大值约2.2 cm；最大沉降量约6.4 cm，位于断层上盘附近。2.5D形变场结果表明，双湖地震以正断特性为主，发震断层属于东倾断层。此次地震反演得到的断层倾角较小，表明地震的成因不仅受自身重力影响，还可能受到区域拉张作用的影响，具体情况需要作进一步研究。双湖地区是羌塘盆地最具油气资源潜力的地区之一[13]，本次地震确定的地下隐伏断层相关参数有利于区域内部结构的探测，可为相关资源的开采提供参考。
|
|
表 3 不同机构给出的震源参数 Tab. 3 Focal parameters retrieved by different agencies |
|
图 7 2.5D形变场 Fig. 7 2.5D deformation filed |
本文基于Sentinel-1A干涉宽幅数据，利用DInSAR技术提取双湖地震的同震形变场，并通过SMC方法反演发震断层参数及同震滑动分布，得到以下结论：
1) T114升轨、T121降轨的LOS形变场均显示双湖地震以沉降为主，升降轨LOS最大沉降量分别约5.8 cm和4.8 cm。
2) 反演得到发震断层走向约47.810°±2.218°，倾角约36.664°±2.499。
3) 断层滑动分布结果表明，双湖地震主要破裂区集中在3.9~7.5 km深度，最大滑动量约0.42 m，反演确定的震中位置为87.715° E、34.365° N，震源深度为5.673 km，平均滑动角约-66.598°±3.258°，表明本次地震以正断拉张为主，兼少量左旋走滑性质。假设地壳的剪切模量为30 GPa，计算得出的矩震级为MW5.5，略小于USGS的结果。
致谢: 感谢欧洲空间局(ESA)提供Sentinel-1卫星升降轨SAR数据。
| [1] |
张培震, 邓启东, 张国民, 等. 中国大陆的强震活动与活动地块[J]. 中国科学: 地球科学, 2003, 33(增1): 12-20 (Zhang Peizhen, Deng Qidong, Zhang Guomin, et al. Active Tectonic Blocks and Strong Earthquakes in the Continental of China[J]. Science China: Earth Sciences, 2003, 33(S1): 12-20)
(  0) |
| [2] |
Vasyura-Bathke H, Dettmer J, Steinberg A, et al. The Bayesian Earthquake Analysis Tool[J]. Seismological Research Letters, 2020, 91(2A): 1 003-1 018 DOI:10.1785/0220190075
( 0) |
| [3] |
Rosen P A, Gurrola E, Sacco G F, et al. The InSAR Scientific Computing Environment[C]. The 9th European Conference on Synthetic Aperture Radar, Nuremberg, 2012
( 0) |
| [4] |
Farr T G, Rosen P A, Caro E, et al. The Shuttle Radar Topography Mission[J]. Reviews of Geophysics, 2007, 45(2)
( 0) |
| [5] |
Goldstein R M, Werner C L. Radar Interferogram Filtering for Geophysical Applications[J]. Geophysical Research Letters, 1998, 25(21): 4 035-4 038 DOI:10.1029/1998GL900033
( 0) |
| [6] |
Chen C W, Zebker H A. Phase Unwrapping for Large SAR Interferograms: Statistical Segmentation and Generalized Network Models[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(8): 1 709-1 719 DOI:10.1109/TGRS.2002.802453
( 0) |
| [7] |
Moral P D, Doucet A, Jasra A. Sequential Monte Carlo Samplers[J]. Journal of the Royal Statistical Society, 2006, 68(3): 411-436 DOI:10.1111/j.1467-9868.2006.00553.x
( 0) |
| [8] |
Zheng A, Yu X W, Xu W B, et al. A Hybrid Source Mechanism of the 2017 MW6.5 Jiuzhaigou Earthquake Revealed by the Joint Inversion of Strong-Motion, Teleseismic and InSAR Data[J]. Tectonophysics, 2020, 789: 228 538 DOI:10.1016/j.tecto.2020.228538
( 0) |
| [9] |
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) |
| [10] |
Heimann S, Vasyura-Bathke H, Sudhaus H, et al. A Python Framework for Efficient Use of Pre-Computed Green's Functions in Seismological and other Physical Forward and Inverse Source Problems[J]. Solid Earth, 2019, 10(6): 1 921-1 935 DOI:10.5194/se-10-1921-2019
( 0) |
| [11] |
Fukuda J, Johnson K M. A Fully Bayesian Inversion for Spatial Distribution of Fault Slip with Objective Smoothing[J]. Bulletin of the Seismological Society of America, 2008, 98(3): 1 128-1 146 DOI:10.1785/0120070194
( 0) |
| [12] |
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) |
| [13] |
尹青, 伊海生, 夏国清. 羌塘盆地双湖地区侏罗系岩石物性特征[J]. 地球物理学进展, 2015, 30(1): 285-292 (Yin Qing, Yi Haisheng, Xia Guoqing. Characteristics of the Petrophysical Parameters of the Rocks in Shuanghu Region, Qiangtang Basin[J]. Progress in Geophysics, 2015, 30(1): 285-292)
( 0) |
2. Key Laboratory of Earthquake Geodesy, CEA, 40 Hongshance Road, Wuhan 430071, China;
3. Hubei Earthquake Agency, 48 Hongshance Road, Wuhan 430071, China;
4. Hebei Key Laboratory of Earthquake Dynamics, 465 Xueyuan Road, Sanhe 065201, China