基于InSAR约束的2022年泸定MW6.6地震同震滑动分布及库仑应力变化

王 欣 李水平,2,3 宋舜跃

1 合肥工业大学土木与水利工程学院，合肥市屯溪路193号，230009 2 武汉引力与固体潮国家野外科学观测研究站，武汉市洪山侧路40号，430071 3 中国地质大学(武汉)地球物理与空间信息学院，武汉市鲁磨路388号，430074

据中国地震台网测定，2022-09-05四川泸定县发生MS6.8地震，震中位于29.59°N、102.08°E(图1)，震源深度16 km。泸定地震是该地区自1955年康定7.5级地震以来的又一次强震，也是近年来造成人员伤亡最严重的一次地震，受到社会的广泛关注。GCMT、USGS、IPGP、GFZ、CENC等机构发布了此次地震的震源机制解(表1)，张龙等[1]、韩炳权等[2]也开展了相关研究。虽然所用数据和方法不同，但研究结果总体认为泸定地震断层的走向为NNW-SSE。

表1 泸定地震破裂参数Tab.1 Rupture parameters of the Luding earthquake

此次泸定地震发生在鲜水河断裂磨西段附近，距龙门山断裂带、大凉山断裂带、安宁河断裂带和玉龙希断裂带距离较近，其中鲜水河断裂带、龙门山断裂带和安宁河断裂带交会构成了“Y”字形断裂带。鲜水河断裂带地处川西高山区，是巴颜喀拉块体的西南边界和川滇块体的北边界，地质构造复杂，始于甘孜东谷一带，向南经过炉霍、道孚、康定、磨西一线，至石棉安顺场一带逐渐消失，全长约350 km，总体走向为N320°～330°W，呈略向NE凸出的弧形，是青藏高原东南缘的一条全新世活动强烈的左旋走滑型断裂带[3-4]。根据构造地貌学，晚第四纪以来，鲜水河断裂磨西段的滑动速率为9.6～13.4 mm/a[5]；根据时序InSAR方法，鲜水河断裂带整体滑动速率为8.8～17.9 mm/a，闭锁深度为7.6～18.5 km[6]；根据GPS或水准测量方法估计出的鲜水河断裂北西段滑动速率为8～15 mm/a，南段滑动速率为5～14 mm/a[7]。18世纪以来，整条鲜水河断裂带几乎都发生过地表破裂，且复发频繁，有历史记录的7级以上地震共8次[4](图1)。已有研究表明，2008年汶川地震和2013年芦山地震后，鲜水河断裂磨西段的库仑应力有一定程度的加载[8]，说明该地区未来地震危险性较高。

红色五角星为CENC震中位置，黑色震源球为GCMT、USGS、IPGP和GFZ等机构给出的震源机制解，黄色震源球为历史上鲜水河断裂带7级以上地震，深蓝色震源球为龙门山断裂带上的芦山和汶川地震，圆圈为余震，蓝色箭头为GPS速度场图1 泸定地震区域构造背景Fig.1 Tectonic background of Luding seismic area

泸定地区地形复杂，海拔高差大，容易发生滑坡等地质灾害，很难进行常规大地测量以及野外现场考察工作。本文采用D-InSAR获取2022-09-05泸定地震的同震形变场，首先利用贝叶斯方法[9]中的均匀滑动模型搜索断层的先验几何参数，然后采用非负最小二乘方法[10]获得断层面的同震滑动分布，最后通过计算同震库仑应力变化对区域地震危险性进行评估，计算震间应变场对发震断层构造进行讨论。

1 泸定地震InSAR同震形变场

本文使用泸定地震前后Sentinel-1A升降轨影像数据，参数见表2。使用InSAR处理软件GMTSAR，采用二轨法进行SAR数据差分干涉处理，具体流程见图2。

图2 D-InSAR处理流程Fig.2 D-InSAR processing flow chart

表2 Sentinel-1A SAR影像干涉对参数Tab.2 Parameters of Sentinel-1A SAR image interference

图3为升降轨干涉图和同震形变场。由图3(d)可见，2022-08-26～09-19升轨T26同震形变场最大视线向沉降量为15 cm，最大视线向抬升量为13.8 cm；由图3(e)可见，降轨T135同震形变场最大视线向沉降量为7.4 cm，最大视线向抬升量为12.7 cm；由图3(f)可见，2022-08-26～09-07升轨T26同震形变最大视线向沉降量为7.7 cm，最大视线向抬升量为14.9 cm。升降轨沉降量和抬升量存在差异，主要是升降轨观测角度不同所致。本文结果与张龙等[1]和韩炳权等[2]的结果大致相同。需要说明的是，相比于2022-08-26～09-19升轨T26图像，2022-08-26～09-07升轨T26图像存在更大的观测噪声，且存在将一些未发生形变区域错误解缠成有相位区域的情况，因此其可靠性相对较差。为提高后续断层滑动分布反演效率，本文对升降轨同震形变场进行四叉树降采样处理[11]，最终得到1 355 个2022-08-26～09-19升轨T26视线向形变观测数据、1 567个2022-08-26～09-07升轨T26视线向形变观测数据和1 382个降轨T135视线向数据。

图3 InSAR升降轨干涉图和同震形变场Fig.3 InSAR interferogram of ascending and descending and coseismic deformation field

2 断层滑动分布反演

2.1 均匀滑动分布

由于泸定地震缺乏发震断层的先验几何参数，因此在进行分布式滑动分布反演前，需要使用均匀滑动模型获得断层倾角、走向、震源深度等参数。本文采用Bagnardi等[9]开发的GBIS MATLAB软件包，基于升降轨InSAR数据，使用贝叶斯算法进行断层参数估计。在参数估计过程中，对初始模型参数空间不加以严格约束，如将断层走向设为0°～360°、倾角设为-90°～90°，可以保证搜索结果的可靠性。迭代106次后，最终得到断层模型参数的后验概率分布见图4。由图可见，断层最优长度为19.5±0.06 km、宽度为15.29±0.15 km、深度为15.61±0.13 km、走向为166.86°±0.09°、倾角为74.23°±0.26°、走滑分量为0.905±0.005 m、倾滑分量为-0.095±0.001 m，表明本次地震以左旋走滑为主。这与韩炳权等[2]GCMT和IPGP等给出的结果基本一致(表1)。

图4 断层模型参数的后验概率分布Fig.4 The posterior probability distribution of fault model parameters

2.2 分布式滑动分布

由于均匀滑动分布无法得到精细的分布式滑动断层模型，因此采用非负最小二乘法求解断层滑动与地表形变的关系[12]，反演的目标函数为：

‖W(Gs-d)‖2+k2‖Ls‖2=min

(1)

式中，W为观测值权阵，G为格林函数，s为每个断层单元上的滑移量，d为地表观测值，k2为平滑因子，L为拉普拉斯二阶有限差分算子。

由于平滑因子对分布式滑动分布反演结果影响较大，因此固定其他参数以测试平滑因子的大小，最终确定平滑因子为35。为构建断层面的精细同震破裂模型，对断层长度和深度分别进行适度扩展，以消除反演结果的边缘效应，并将断层面划分为1 km×1 km的矩形单元。

建立分布式断层滑动模型是同震反演的关键步骤。由于2022-09-07的升轨数据质量存在问题，为测试其对反演结果的影响，对本次地震升降轨数据进行以下2种实验：1)使用2022-08-26～09-07升轨T26和降轨视线向数据进行断层滑动分布反演；2)使用2022-08-26～09-19升轨T26和降轨视线向数据进行断层滑动分布反演。基于2种实验方式分别构建分布式滑动分布模型，然后根据数据模拟效果评估反演结果。

以2022-08-26～09-07升轨T26和降轨T135同震形变场为约束进行滑动分布反演，将断层走向设置为NNW-SSE，对应的数据模拟结果见图5。由图可见，升轨的观测值和模拟值显著不同，存在明显的残差，残差主要集中在震中东侧；降轨的观测值和模拟值吻合较好。升轨数据在震中东侧产生明显残差的原因可能是该地区植被茂盛，且地形梯度较大，使得升轨图像中存在较为明显的大气误差。

图5 升降轨同震形变场观测值、模拟值和残差分布Fig.5 Observed and simulated values and residual distribution of the coseismic deformation of ascending and descending

利用2022-08-26～09-19升轨T26和降轨T135同震形变场作为约束，反演同震滑动分布，数据模拟结果见图6。由图可见，升轨平均残差为1.2 cm，降轨平均残差为1.8 cm，残差图中不存在明显的系统误差，表明同震滑动模型能较好地模拟地表InSAR观测数据。

图6 升降轨同震形变场观测值、模拟值和残差分布Fig.6 Observed and simulated values and residual distribution of the coseismic deformation of ascending and descending

图6中对应的二维和三维同震滑动分布模型见图7。反演结果显示，断层滑动主要以走滑为主，兼少量逆冲分量，此结果与GCMT、USGS等给出的结果大致相同。地震破裂深度主要集中在0～17 km，破裂沿走向延伸约55 km，大于0.8 m的较大量级滑动主要发生在0～7 km深度处，断层最大滑动量为1.12 m，对应深度为1 km，释放的总地震矩为1.02×1019Nm，相当于MW6.64，此结果与GCMT计算的结果大致相同。由表1可见，不同类型观测资料约束的地震震级和发震断层参数基本一致。余震活动主要发生在同震破裂区的中部、深部和边缘扩展区，总体上看，余震分布与同震破裂具有较好的空间互补关系，说明泸定地区余震活动是由地震同震破裂触发的应力加载所致。

图7 泸定地震滑动分布模型Fig.7 Slip distribution model of Luding earthquake

3 同震库仑应力变化

地震的破裂滑动会使区域应力场发生改变，进而影响余震活动，因此余震活动分布可以验证破裂滑动分布结果的正确性。分析震中CFS的变化是预测震中附近发生地震可能性的重要方法之一。利用反演的断层滑动分布参数，运用Coulomb3.3软件包计算泸定地震附近地区的库仑应力变化。

CFS=Δτ+μ′Δσn

(2)

式中，Δτ为剪切应力变化，μ′为有效摩擦系数，Δσn为法向应力变化。

将泊松比设为0.25，有效摩擦系数设为0.4，分别计算5 km、10 km、15 km和20 km深度处的同震库仑应力变化，结果见图8。由图可见，随着深度增加，库仑应力在震中附近及靠南地区发生显著变化，5～15 km深度范围内大部分余震分布在库仑应力加载区域，符合地震应力传输和同震库仑应力触发机制；小部分余震活动分布在应力卸载区域，表明此次地震未能释放断层面积累的全部应力。余震分布趋势为NNW-SSE向，与断层滑动分布走向大致相同。泸定地震发生后，库仑应力的增加主要集中在破裂断层的西北段、东南段以及西南段，位于鲜水河断裂带、安宁河断裂带北段、玉龙希断裂中北段。上述地区未来发生地震的危险性很大，尤其需要注意安宁河断裂带北段，该断裂带在1480年和1536年发生7.5级地震、在1952年发生6.75级地震，是需要重点关注的地震空区[13]。

红色五角星为震中位置，白色实线为断层位置，白色虚线为断层俯视图，黑色实线为断层，黄色圆圈为不同深度的余震分布图8 泸定地震引起的CFS变化Fig.8 Variation of CFS caused by Luding earthquake

4 震源区孕震特征

泸定地震发生在鲜水河断裂带东南端，毗邻鲜水河断裂带、龙门山断裂带、安宁河断裂带的交会区。分析震源区震前长期应变特征，可为识别强震孕育阶段以及预测强震提供重要参考。本文利用球面多尺度小波变换方法[14]，以GPS速度场为约束，建立泸定地震震源区地壳应变场模型。GPS速度场数据来源于Wang等[15]给出的鲜水河地区近20 a的GPS观测结果。球面多尺度小波分解时，球面网格的最大分辨率(qmax)对结果的影响较大，根据数据的分布密度，经过多次实验，取qmax=7。此时，最大网格尺度为87.2 km，分辨率低于该阈值的区域设定为不确定区域，在计算应变参数时，忽略不确定区域的结果。图9(a)为鲜水河断裂及邻区主应变率和最大剪应变率分布。可以看出，泸定地震震源区以NE-SW向的挤压应变为主，主压应变率量值约为(20～30)×10-9/a，明显小于鲜水河断裂中段和北段，表明鲜水河断裂挤压变形从西北段到东南段逐渐减弱，此结果与藏东高原整体的东向挤出背景基本一致。最大剪切应变率基本沿鲜水河-安宁河断裂带展布，并且在鲜水河断裂带、龙门山断裂带和丽江-小金河断裂带的交会区附近存在一个应变率高值过渡区，该过渡区的位置与泸定地震的位置基本一致。图9(b)为震源区面膨胀率分布。由图可见，此次泸定地震震源北侧和东侧以挤压应变为主，西南侧以拉张应变为主，震中所在的位置正好位于面应变率从挤压应变到拉张应变的转换区域，该应变转换区的存在可能与多个不同活动块体在该地区交会并产生的相互作用有关。

图9 泸定地震震源区最大剪切应变率和面膨胀率分布Fig.9 Distribution of maximum shear strain rate and surface dilatation rate in the focal area of Luding earthquake

5 结 语

1)泸定地震升降轨InSAR同震形变场中，升轨最大视线向形变量为-15.0 cm，降轨最大视线向形变量为12.7 cm。

2)基于视线向同震形变场反演的滑动分布模型显示，发震断层以左旋走滑为主，兼少量逆冲，断层走向沿NNW-SSE，约为167°，沿走向破裂约55 km，倾角约为74°，断裂主要集中在0～17 km深度处，最大滑动量约为1.12 m，对应深度为1 km。泸定地震释放的总地震矩为1.02×1019Nm，对应矩震级为MW6.64。

3)同震库仑应力结果显示，泸定地震使鲜水河断裂带南东段、安宁河断裂带北段、玉龙希断裂中北段处于应力加载状态，其中安宁河断裂带北段库仑应力加载状态明显，未来发生地震的风险较大，需要开展持续监测和危险性评估。

4)震间应变场表明，此次泸定地震震源区位于拉张应变和挤压应变的转换区域，该应变转换区的存在可能与不同活动块体在该地区的交会有关。