2021年双湖MS5.8地震InSAR同震形变场及断层滑动分布反演

李承涛 李 琦 谭 凯 鲁小飞

1 中国地震局地震研究所，武汉市洪山侧路40号，4300712 中国地震局地震大地测量重点实验室，武汉市洪山侧路40号，4300713 湖北省地震局，武汉市洪山侧路48号，4300714 河北省地震动力学重点实验室，河北省三河市学院街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

1 InSAR数据处理

本文采用欧空局官方网站(https:∥sentinel.esa.int/)提供的Sentinel-1A干涉宽幅数据，具体信息见表1。

表1 Sentinel-1A数据相关参数

基于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 Shuanghuearthquake

2 断层参数和滑动分布反演

2.1 断层参数反演

为获取发震断层相关的几何形状参数，利用基于序列蒙特卡罗采样的贝叶斯方法进行反演[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。

表2 反演的断层几何形状参数

基于上述断层模型得到的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

2.2 断层滑动分布反演

在已知断层几何参数的基础上，将断层离散为若干子断层，断层滑动分布反演可简化为求解每个子断层对应的滑动。利用基于Python语言的Pyrocko-GF软件生成格林函数库，然后通过基于序列蒙特卡罗采样的贝叶斯方法反演滑动分布[2,10]，子断层滑动分布与相应格林函数库的乘积即为地表位移。每个子断层的滑动量搜索范围设置为0～1.0 m。

在贝叶斯方法中，正则化可通过高斯先验p(s,α)实现，其包含一个非对角项的协方差矩阵，Cholesky分解等价于一个大小为M×M的拉普拉斯有限差分算子L(M为子断层个数)。p(s,α)可表示为[2]：

(1)

正则化后的后验概率密度函数PDF可表示为[11]：

p(m,σ1,σ2,…,σK|dobs,1,dobs,2,…,dobs,K)

(2)

式中，m为已知断层参数，σK为第K个超参数，dobs,K为第K个观测数据，平滑因子α可作为待求参数，s为子断层待求的滑动矢量。

直方图中红色直线和2D相关图中的红点表示最大后验解，蓝色表示高概率区域图3 双湖地震矩形震源的1D和2D后验概率密度分布Fig.3 1D and 2D PDF of a rectangular source for Shuanghu earthquake

在反演断层滑动分布时，参考表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×1017Nm，矩震级为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 thefault slip distribution model

3 讨 论

表3为不同机构提供的震源参数，由于采用的数据和方法不同，获得的参数存在差异。反演结果表明，发震断层属于具有正断性质的隐伏断层，震中位于琵琶湖-吐坡错断裂与琵琶湖-映天湖断裂的交汇处，地震可能受这2个断裂共同作用。Sentinel-1卫星对南北方向不敏感，当南北向分量近似忽略时，通过升降轨数据可以获取近东西向和垂直向的形变分量，称为2.5D形变场[12]，结果如图7所示，图中灰色矩形框表示滑动分布的断层面，黑色和红色五角星分别表示CENC和InSAR计算的震中位置，后者更靠近形变中心。从图7可以看出，东向运动最大值约2.8 cm，西向运动最大值约2.2 cm；最大沉降量约6.4 cm，位于断层上盘附近。2.5D形变场结果表明，双湖地震以正断特性为主，发震断层属于东倾断层。此次地震反演得到的断层倾角较小，表明地震的成因不仅受自身重力影响，还可能受到区域拉张作用的影响，具体情况需要作进一步研究。双湖地区是羌塘盆地最具油气资源潜力的地区之一[13]，本次地震确定的地下隐伏断层相关参数有利于区域内部结构的探测，可为相关资源的开采提供参考。

表3 不同机构给出的震源参数

图7 2.5D形变场Fig.7 2.5D deformation filed

4 结 语

本文基于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数据。