引言

据中国地震台网中心（2021）正式测定，北京时间2021年5月21日，云南省西部大理白族自治州漾濞彝族自治县发生M_S6.4地震，震中位置为（99.87°E，25.67°N），震源深度为8.0 km，此次地震给人民的生命财产带来极大的损失。漾濞地震是一次典型的前-主-余震型序列地震，该序列从5月18日持续到5月31日。漾濞地区具有频繁发生中强地震的记录，在此次地震震中100 km范围内曾发生上百次M5.0—7.0地震（国家地震局震害防御司，1995）。

漾濞地震发生后，国内外科研机构和学者对该地震的同震形变场、震源机制及整个地震序列等开展了一系列研究。张克亮等（2021）通过连续的全球卫星导航系统（Global Navigation Satellite System，缩写为GNSS）数据处理获得震区50 km范围内的同震形变场并反演获得地震的破裂滑动分布，但由于台站分布稀疏，近场观测有效数据较少，获得的断层面滑动分布结果较为简单；朱俊文等（2021）利用双轨差分雷达干涉法获取了地震同震形变场，他们通过实验获取的合成孔径雷达干涉测量（interferometric synthetic aperture radar，缩写为InSAR）数据虽然在空间上覆盖了震中区域，但在时间上涵盖了部分前震和许多余震，得到的同震形变场是前-主-余震共同作用的结果；杨九元等（2021）的反演结果表明，发震断层以右旋走滑运动为主，正断层分量较少，认为地震破裂发生在维西—乔后—巍山断裂的隐伏次生断裂或独立的未知隐伏主断层上，同震库仑应力的变化也表明部分区域具有较高的破裂风险。国内在获取同震地震形变场时主要采用的方法有传统大地测量技术、GNSS （宁津生，王正涛，2006）、电子激光测距（electronic distance measurement，缩写为EDM）（王若柏等，2000）和InSAR （Peltzer，Rosen，1995）等技术。传统的大地测量技术获取震源位置主要依赖于观测台站与震源的相对位置，难以实现震源的有效定位。我国已经在诸多地方布设了高密集的观测台站网络，但在一些人迹罕至、环境恶劣地区，仍存在着大面积的空白区域。近几十年，随着InSAR数据的免费使用及其在地震、滑坡等地表形变监测方面的良好表现，卫星大地测量技术极大地弥补了传统大地测量技术的不足。

目前研究多集中于通过震源机制解给出的断层几何参数来反演断层面的滑移分布或单一数据源构建模型，例如InSAR数据仅能观测到一维形变量即雷达视线向（line-of-sight，缩写为LOS）的形变，即地表三维形变在视线向的投影，并不能获取真实的地表三维形变，因此，并不能很好地分析整个发震断层的真实时空解析度。鉴于此，为弥补InSAR数据单一视线向的不足，本文基于欧空局的Sentinel-1A卫星震前、震后升降轨数据，采用合成孔径雷达差分干涉测量技术（differential interferometry synthetic aperture radar，缩写为 DInSAR）技术重建漾濞M_S6.4地震高精度、高密度同震形变场，从而获取形变区域真实的同震形变位移。结合Okada弹性半空间位错模型（Okada，1985）进一步拟合漾濞地震的断层几何参数。联合InSAR同震形变场结果以及震后大地测量观测数据进一步反演发震断层面的滑动分布。最后，结合余震的分布区域和模拟的断层迹线走向以及同震库仑应力减小的展布方向，确定发震断层的确切位置，为后续该区域的地学研究提供参考。

1.   研究区与数据来源

1.1   研究区概况

川滇菱形地块位于青藏高原东南缘，因受到印度板块、欧亚板块以及刚性的四川盆地连续碰撞阻断而造就了强烈的区域性大型走滑断层活动（徐锡伟等，2003），漾濞地震的震中则位于川滇菱形地块的西南边界，该边界也是维西—乔后—巍山断裂与红河断裂带的交界区域。本文以漾濞震区为研究区，该区域附近断层复杂繁多，是导致区域地震频发的重要原因。红河断裂带是一条典型的右旋走滑断裂，由北至南分为北、中、南三段，总长度为1 000 km （邓起东，1996；张建国，2009）。维西—乔后—巍山断裂因其与红河断裂带具有非常相似的运动特征，被认为是红河断裂带的北延部分，南接红河断裂带，北连金沙江断裂（常祖峰等，2016），距离此次云南漾濞地震的震中仅约10 km。图1为云南漾濞区域构造背景，蓝色虚框为Sentinel-1A升轨数据覆盖区域，黄色虚框为Sentinel-1A降轨数据覆盖范围。

图  1  云南漾濞区域构造背景

Figure  1.  Tectonic background of Yangbi region in Yunnan Province

下载: 全尺寸图片 幻灯片

1.2   研究数据

Sentinel系列卫星具有高空间分辨率、高光谱、多功能等优点。本文选取空间分辨率为5 m×20 m的Sentinel-1A干涉宽幅（Interferometric Wide swath，缩写为IW）模式的斜距单视复数产品，极化方式为同极化；从中选用主震前后4景升、降轨合成孔径雷达（synthetic aperture radar，缩写为SAR）影像作为数据源，其参数详见表1；时间基线均为12 d，这是由于较短的时间间隔有助于增强主、从影像间的相干性，获得高时空分辨率的同震形变场。选用覆盖研究区的AW3D30 （ALOS World 3D 30 m） DEM数据作为参考DEM，该数据由日本宇宙航空研究开发机构（Japan Aerospace Exploration Agency，2021）提供。

表  1  合成孔径雷达影像数据参数

Table  1.  Data parameters of synthetic aperture radar image

轨道方向 （轨道号）	获取时间		雷达入射角/°	轨道方位角/°	空间基线/m	时间基线/d
	震前	震后
升轨（T99）	2021-05-20	2021-06-01	39.25	−12.45	19.20	12
降轨（T135）	2021-05-10	2021-05-22	39.35	−167.58	−50.35	12

下载: 导出CSV 

| 显示表格

2.   同震形变场建立

合成孔径雷达差分干涉测量技术（DInSAR）是通过对同一地区前后两次获取的不同SAR图像干涉结果进行差分处理的一种方法。本文基于双轨DInSAR技术对地震前后两景影像进行处理，具体流程如图2所示。首先，为减少轨道误差的影响，引入精密轨道文件进行去噪等一系列数据预处理之后，对震前震后的SAR图像进行配准、计算干涉相位、生成干涉图、去除平地效应等处理得到去平干涉图。其次，采用Goldstein滤波方法（李佳等，2011）对上一步生成的差分干涉图进行滤波，这种滤波方法的滤波器是可调的，有助于提高干涉条纹的清晰度、降低空间或者时间基线引起的失相干噪声，并且在一定程度上抑制干涉图噪声。最后，基于最小费用流（minimum cost flow）相位解缠方法（Orlin，1993）实现绝对相位恢复，结合地面控制点（ground control point，缩写为GCP）的相位，进行多项式拟合，并对干涉图和解缠后的相位进行轨道精炼和重去平，将优化后的相位结果进行地理编码。

图  2  DInSAR技术流程

Figure  2.  Technical process of DInSAR

下载: 全尺寸图片 幻灯片

利用双轨法DInSAR技术对2021年云南漾濞地震前后两景SAR影像处理得到图3所示的干涉条纹图和同震形变场。图3a和3c为升降轨的干涉条纹图，整体连续光滑，干涉条纹以发震断层为中心呈同心圆形态分布，说明发震断层周围的位错量沿断层不均匀分布。由于升降轨观测方向不一致，发震断层两侧的形变符号相反，升轨影像时间包含震后10天的大量余震，使得处理后的形变场结果较差，两侧形变特征不明显只能识别出一个较大的形变中心。如图3b所示，该区域靠近卫星视线向，最大视线向形变量约为7.00 cm。从图3d可知，降轨影像完整地捕捉到漾濞地震的同震形变场，在震中东南部显示了两个主要的叶状形变场，西南方向远离卫星运动，最大视线向形变量约为−7.00 cm，东北方向靠近卫星运动，最大视线向形变量约为7.80 cm。此外，InSAR降轨数据得到的干涉条纹图以及同震形变场整体上呈椭圆形对称分布，形成了一条NW-SE走向的迹线，断层两侧形变差异明显，故推测该迹线很有可能为发震断层。

图  3  漾濞M_S6.4地震升（上）、降（下）轨干涉条纹图（a，c）及同震形变场（b，d）

Figure  3.  The interferograms （a，c） in ascending orbit （upper） and descending orbit （lower） and coseismic deformation field （b，d） of Yangbi M_S6.4 earthquake

下载: 全尺寸图片 幻灯片

为了进一步确定发震断层的形变特征，沿NE-SW方向对垂直于发震断层的升、降轨形变场绘制剖面（图4），可以观察到升降轨数据得到的剖面整体方向相反。InSAR同震形变场数据量较大，同时相邻观测值间具有很大的空间相关性。通过利用最近邻法对InSAR同震形变场进行降采样处理，InSAR观测值庞大的数据量得到了大幅筛减，同时保留了原始形变场的主要特征，从而提高了整体的反演效率。

图  4  升轨剖面AA′（a）和降轨剖面BB′（b）

Figure  4.  Ascending profile AA′ （a） and descending profile BB′ （b）

下载: 全尺寸图片 幻灯片

3.   断层滑动分布反演

3.1   断层几何参数

通过Sentinel-1A卫星升降轨同震形变场形变位置大致勾勒出一条NW-SE走向的迹线（图5）。为确定漾濞地震发震断层的位置，基于中国地震台网中心（2021）提供的2021年5月18日至2021年5月31日的震相资料，采用双差相对定位方法HypoDD （Zhang，Thurber，2003）对云南漾濞地震序列进行重定位。漾濞地震前后，共观测到4 451次地震，共有1万1 444条P波，8526条S波记录。由于漾濞地区台站间距较大且分布零散，可用于重定位的地震数量有限，参考胡鸿翔等（1986）和张中杰等（2005）的人工地震测深结果，采用层状速度结构模型对地震序列进行重定位。根据震中周围台站的接收函数h-k扫描结果（王长在等，2011），波速比设定为1.76。采用双差方法进行地震重定位前，将从中国地震台网中心（2021）下载的观测报告转为震相数据，根据区域台网数据选择合适的台站和震相数据进行走时处理，并将此作为重定位的初始位置。重定位后得到主震和4450次余震的震源位置参数，南北向平均误差约为0.20 km，东西向平均误差约为0.20 km，深度误差约为0.40 km。主震震中位置为（99.89° E，25.67° N），震源深度为13.29 km，除主震之外余震的震源深度主要集中在5—15 km （图5）。余震整体沿NW走向分布，与InSAR形变场结果相一致，其整体分布结果也表明，主震的发震断层走向是沿SE并且是由NW向SE的单向破裂。

图  5  2021年漾濞M_S6.4地震的余震精定位结果

（a） 小震精定位结果；（b） 沿BB′剖面的震源深度分布；（c） P波和S波的走时-震中距分布图

Figure  5.  Precise locating results of the aftershocks of M_S6.4 Yangbi earthquake in 2021

（a） Precise locating result based on small events；（b） Distribution of focal depths of the aftershocks along the seismic profile BB′；（c） Traveltime versus epicentral distance distribution of P and S waves

下载: 全尺寸图片 幻灯片

基于余震精定位结果，结合万永革等（2008）提出的断层面与小震位置拟合方法，求解断层几何参数。该方法的基本原理是：首先利用成丛小震发生在主震断层面及其周围的原则，参考崔华伟等（2022）的区域构造应力参数，估计断层面的滑动角等；接着通过模拟退火算法和高斯-牛顿算法，利用小震密度求解主震断层面走向、倾角、位置并对其误差进行稳健估算。根据小震位置拟合的漾濞地震主断层面的几何参数为：走向316.69°，标准差为0.42°；断层倾角88.56°，标准差为1.03°；滑动角177.97°。综合考虑余震序列的整体分布（图1）、升降轨同震形变场（图3）以及美国哈佛大学（GCMT）和美国地质调查局（USGS）公布的震源机制解（表2），本文结果显示此次地震主震的震中位置为（99.89°E，25.67°N），震源深度为13.29 km，地震矩震级为M_W6.1，发震断层走向为NW-SE （316.69°），倾角为88.56°，滑动角为177.97°，为右旋走滑断裂。

表  2  不同机构给出的云南漾濞地震断层参数

Table  2.  Fault parameters of Yangbi，Yunnan earthquake given by different institutions

机构	震中位置			节面Ⅰ				节面Ⅱ			深度/ km	MW
	东经/°	北纬/°		走向 /°	倾角/°	滑动角/°		走向/°	倾角/°	滑动角/°
GCMT （2021）	100.01	25.62		45	80	7		314	83	170	15.10	6.10
USGS （2021）	100.01	25.73		43	75	−9		135	82	−165	9.00	6.09
本文	99.89	25.67						316.69	88.56	177.97	13.29	6.10

下载: 导出CSV 

| 显示表格

3.2   断层破裂滑动模型反演

地震的发生、地壳形变以及地表形变之间总是相互影响、相互制约的。基于弹性半空间位错模型（Okada，1985），以升降轨数据获得LOS向形变场为数据源，通过小震精定位拟合的断层几何参数，采用基于InSAR约束条件下的最小二乘原理和最速下降法（steepest descent method，缩写为SDM）进行发震断层的同震滑动分布反演。

断层倾角参考余震精定位结果，选定断层滑动模型长40 km，宽20 km，断层走向为316.69°，断层倾角为88.56°，断层顶部埋深为3.44 km，设置滑动角范围为−180°至180°，平滑因子升轨为0.25，降轨为0.50。为更好地获取断层面滑动分布细节，将断层沿走向和倾角离散成40×20个1 km×1 km的网络。地壳速度结构模型采用杨军等（2014）反演云南地区中小地震震源机制解所采用的莫霍面深度为55 km的模型（表3），基于线性反演算法获得漾濞地震断层滑动分布（图6）。

表  3  本文所用的地壳速度结构模型（杨军等，2014）

Table  3.  Crustal velocity structure model used in this study （Yang et al，2014）

地壳分层 /km	vP /（km·s−1）	vS /（km·s−1）	密度 /（kg·cm−3）
0—3	5.50	2.50	2500
3—18	6.10	3.45	2810
18—25	6.35	3.55	2880
25—39	5.70	3.35	2610
39—55	6.65	3.65	2970
≥55	7.80	4.30	3280

下载: 导出CSV 

| 显示表格

图  6  升（a）、降（b）轨断层滑动模型

Figure  6.  Fault slip models in ascending （a） and descending （b） orbits

下载: 全尺寸图片 幻灯片

由图6可见：断层滑动以走滑为主，破裂深度主要集中在5—15 km；升轨断层最大滑动量为0.80 m，对应的深度为8.85 km，平均滑动量为0.22 m，拟合度为0.86，矩震级为M_W6.41；降轨的断层最大滑动量为0.30 m，对应的深度为6.88 km，平均滑动量为0.05 m，拟合度为0.88，矩震级为M_W6.01。由InSAR获取的同震形变场和残差的分布（图7）可知，模拟得到的形变场与观测值具有较好的一致性，升降轨影像的残差较大，可能其中包含大量的余震余滑。

图  7  升（上）、降（下）轨InSAR模型观测值（a）、模拟值（b）和二者残差图（c）

Figure  7.  Plot of observed （a），simulated （b） and residual values （c） for the ascending （upper） and descending （lower） orbits InSAR models

下载: 全尺寸图片 幻灯片

3.3   同震库仑应力计算

中强震同震滑移往往引起周围断层应力场的变化，导致未来的地震事件提前或延迟发生（Symithe et al，2013）。作为近45年来维西—乔后—巍山断裂附近发生的最大地震，计算和分析漾濞地震同震库仑应力变化是评估该地震对周边断层影响的重要手段。根据库仑破裂准则，库仑应力变化$\Delta \delta_{\mathrm{f}} $定义为

式中，Δτ代表断层面滑动方向上剪应力的变化，Δσ_n表示断层面上的法向正应力变化，μ′表示有效摩擦系数，取值范围为0—1 （Scholz，2002）。当库仑应力变化$\Delta \delta_{\mathrm{f}} $为正时，应力加载会促进断层的破裂，反之应力卸载则抑制断层的破裂。

基于反演的断层破裂滑动模型，利用Coulomb3.3软件，接收断层参数选用余震精定位的结果（走向316.69°，倾角88.56°，滑动角177.97°），采用有效摩擦系数0.40 （King et al，1994）计算已知断层上的同震库仑应力变化，结果如图8所示。结果表明：在7.50 km深的截面上，漾濞地震主要表现为库仑应力对周边断层区域的卸载作用，周围断层出现了少量的应力加载作用；在15 km 深的截面上，库仑应力对发震断层本身及其周围断层的应力变化影响不大。漾濞地震后，区域应力差水平降低，缓解了周围断裂的应力和应变积累。在距离主破裂带7.50 km处，余震展布方向与库仑应力减小方向一致，表明发震断层仍处于应力调整阶段，漾濞M_S6.4地震序列对周围的右旋走滑断裂具有显著的应力卸载作用，能够减轻周围断层地震发生的危险性。

图  8  升（上）、降（下）轨不同深度h处的断层库仑应力变化

Figure  8.  Coulomb stress changes of coseismic fault at different depth h in ascending （upper） and descending （lower） orbits

（a） h＝7.5 km；（b） h＝15 km

下载: 全尺寸图片 幻灯片

4.   讨论与结论

本次利用升、降轨同震观测数据的反演结果显示，在破裂断层中段区域缺乏断层近场数据，导致滑动反演时约束不足，但是对反演结果的整体准确性影响不大。由于InSAR原始影像覆盖时段内包含较多的余震，获得的InSAR同震形变场位移是观测期间所有形变量之和，因而反演得出的滑动分布最大滑动量和震级M_W6.01都相比USGS公布的结果略大。震后多学科综合科考调查初步成果（中国地震局地质研究所，2021）显示：野外地质调查未发现明显地表裂缝，同时维西—乔后断裂在此次地震过程中未见同震变形，这表明该断裂在漾濞地震中并未活动。对于主震震源位置，不同机构给出的结果略有不同，尤其是震源深度。反演深度不同可能与不同机构所采用的数据资料波长不同有关。USGS使用的是长周期波长，GCMT使用的是地幔波，中国地震局使用的是区域台网波形资料（王莹等，2021）。波长越短，对矩心的分辨率越高（高原等，1997）。使用区域台网波形资料反演的震源矩心深度分辨率要高于USGS和GCMT利用远震波形给出的矩心深度（易桂喜等，2017；Long et al，2019）。中国地震台网中心给出的速报结果为（99.87°E，25.67°N），震源深度为8.0 km，而本研究通过余震精定位给出的结果为（99.89°E，25.67°N），除主震外的震源深度范围为5—15 km，这与赵博等（2022）通过双差定位取得的结果一致。

基于升、降轨InSAR同震形变场反演得到同震滑动破裂模型，但在以下两个方面还有很大的改进空间：首先，仅通过升、降轨InSAR观测数据作为约束，对于滑动分布结果的分析较为片面，缺少高精度的数据作为约束，下一步将把GNSS三维同震位移数据与实验结果进行融合，深化对本次地震的发震破裂、孕震机制的认识，为本地区地震风险评估提供更加可有力的数据支持；其次，本文反演所采用的断层模型具有局限性，与真实的发震断层有较大的差异，目前仅根据地震序列和相关机构给出的断层数据进行推测，简单地将断层设置为单一的走向或者倾角，模型相对简单。之后可以测试不同断层倾角和倾向对观测数据拟合效果的影响，也可以通过模型的正演进行结果比对。

漾濞地震为右旋走滑型地震，InSAR同震形变场也表现为NW倾向的断层走向。该地震的发震断层是距离维西—乔后—巍山断裂10 km处的一条未知断裂，联合反演得到的发震断层与维西—乔后—巍山断裂的总体走向和滑动形式基本一致。此外，根据余震序列分布勾勒的一条发震断层迹线几乎平行于维西—乔后—巍山断裂，断层的西侧有大量余震发生（图5）。对于此次地震的发震断层目前有两种观点：一种认为漾濞地震的发震断裂是维西—乔后—巍山断裂的隐伏分支断层，在川滇菱形地块向南挤出的构造背景下，地块西边界的红河断裂带和维西—巍山—乔后断裂发生右旋走滑，本次地震是维西—乔后—巍山断裂南段分支断裂右旋走滑的体现（朱俊文等，2021）；另一种观点认为该断层位于一条独立的、未知的隐伏走滑断层上，与维西—巍山—乔后断裂无关（杨九元等，2021）。本文结合余震的分布区域和模拟的断层迹线走向，以及同震库仑应力减小的展布方向，基本认定本次地震的发震断层是维西—乔后—巍山断裂的隐伏次生断裂。

本文结合升、降轨InSAR数据，基于DInSAR技术重建2021年5月21日云南漾濞升降轨同震形变场，通过双差定位方法获得了主震位置并利用最速下降法确定了断层几何参数。联合同震形变场和小震位置分布数据，构建断层破裂滑动模型，运用SDM方法反演本地地震的破裂分布，对断裂的时空解析度进行了分析，主要结论如下：

1） 基于多源InSAR数据获得2021年5月21日云南漾濞地震的精细同震形变场。升轨雷达视线向最大形变量约为5.00 cm，降轨雷达视线向向最大形变量约为7.80 cm。

2） 余震精定位的主震震中位置为（99.89°E，25.67°N），震源深度为13.29 km，除主震之外的震源深度主要集中在5—15 km。余震分布呈现NW走向，与InSAR同震形变场结果相一致。发震断层走向NW-SE （316.69°），断层倾角为88.56°，滑动角为177.97°。此次地震为右旋走滑，升轨断层最大滑动量为0.80 m，对应的深度为8.85 km，平均滑动量为0.22 m，矩震级为M_W6.41；降轨断层最大滑动量为0.30 m，对应的深度为6.88 km，平均滑动量为0.05 m，矩震级为M_W6.01。结合野外科考结果，地表没有明显的裂缝，地震未破裂至地表。

3） 经同震形变场特征分析以及余震展布拟合，认为此次地震的发展断层为维西—乔后—巍山断裂的次级断裂。以小震位置拟合的漾濞地震发震断层为接收断层，计算漾濞地震周边断裂的库仑应力变化，结果显示：在7.50 km深度漾濞地震对周边断层有应力卸载作用；在15 km深度时，断层本身的库仑应力和周围断层的应力差明显减小。这次地震的发生使该地区的应力得到释放，后续发生大地震的可能性不大。