2021年云南省漾濞Ms6.4地表二维形变分析

潘建平，赵瑞淇，袁雨馨，蔡卓言，付占宝

(重庆交通大学 智慧城市学院，重庆 400074)

据中国地震台网消息，北京时间2021年5月21日T21:48云南省大理州漾濞县(99.87°N，25.67°E)发生Ms6.4地震，震源深度8 km，地震致3人遇难，37人受伤，并造成不同程度损毁。漾濞地震发生在横断山滇西高山峡谷区，地质构造复杂，地震活动频繁，震中附近断裂构造复杂，断裂众多[1]。因此，确定本次地震的同震形变场，对认识该地区的地震孕震环境、同震滑动模型反演以及地震风险防范具有重要意义，过去对地表形变监测大都基于水准仪、全站仪、GNSS等设备的外业手段，不能在较短时间内获得大范围的、监测点分布密集的观测数据[2],因而不适用于地震这种大范围高烈度的自然灾害形变监测。

合成孔径雷达差分干涉测量(differential synthetic aperture radar interferometry，D-InSAR)具有大范围，高精度以及低成本的优势，可以监测到cm级或mm级的形变[3]。1993年，测量美国加州Landers地震形变场首次应用D-InSAR技术[4]，从此D-InSAR技术被广泛地应用到地震相关研究中，特别是在震间以及同震形变领域的应用[5]：张庆云[6]等采用D-InSAR技术，与ALOS-2和Sentinel-1A数据联合使用得到了2017年伊朗Mw7.3地震的同震形变场，并在此基础上反演出震发断层的详细几何参数。Elias[7]等采用D-InSAR技术提取希腊Samos岛近海M7.0级地震的同震形变，并联合GNSS偏移量进行了断层参数反演。Zhao L[8]等利用ALOS-2以及Sentinel-1A/B数据得到了菲律宾棉兰老岛4次Mw>6.0地震的InSAR同震形变场，并在InSAR监测结果的基础上对断层几何形状和同震滑动分布进行了几何学约束。Qingyun Zhang[9]等利用InSAR监测结果与强震数据和GPS数据结合得到了2016日本熊本地震的三维形变场。本次地震发生后，我国学者立即采用InSAR技术开展研究，王绍俊[10]等基于D-InSAR技术获得了此次地震的升降轨同震形变场，并以InSAR形变场和GNSS数据为约束，反演了断层滑动分布。杨九元[11]等利用Sentinel-1A/B数据的InSAR同震形变反演了此次地震详细的断层几何参数与滑动分布，并对该区域的地震危险性进行了评估。黎朕灵[12]等利用InSAR升降轨同震形变反演了该地震的震源滑动模型，并与GPS同震位移数据进行了联合反演。此外，另有其他学者利用地震学相关理论对本次地震开展了震源机制反演以及深部孕震背景探讨[13-14]。针对此次地震，文中利用欧空局Sentinel-1A升降轨数据，分别采用D-InSAR技术解算得到了视线向(line of sight,LOS)的同震形变场，并分别进行了形变方向和形变量级的分析。利用美国地质调查局USGS给出的震源机制解和魏云杰[15]等人基于梯度最速下降法得到的同震滑动分布结果模拟计算了本次地震并探讨了本次地震引发的周边断裂库仑应力变化并分析了发震危险性。针对本次地震的现有研究在InSAR数据源上仅利用LOS向的同震形变场，但是LOS向的形变场不能完全反应地表真实形变，因而不能为进一步开展震源参数和同震滑动分布反演提供较好的约束条件，因此文中联合升降轨同震形变场，使用直接解算法得到了垂直向和东-西向的同震形变场,为进一步开展此次地震震源参数和同震滑动模型精细化反演提供更可靠的约束条件。

1 研究方法

1.1 D-InSAR技术原理及数据处理

D-InSAR技术是在InSAR技术基础上发展而来，利用卫星重复轨道获取的两个时相的影像进行差分干涉，得到的干涉相位组成如式(1)所示。

φ=φref+φtop+φdef+φatm+φnoi.

(1)

式中:φref为参考椭球面相位;φtop为地形相位;φdef为形变相位;φatm为大气延迟相位;φnoi为噪声相位。文中利用瑞士SARMAP公司的Sarscape5.2.1软件采用二轨法来提取形变相位，利用卫星姿态参数以及美国宇航局发布的30 m分辨率的SRTM(shuttle radar topography mission)DEM分别去除参考椭球面相位和地形相位，采用10∶2的多视比抑制斑点噪声并提高信噪比，采用GOLDSTEIN滤波方法消除部分噪声相位的影响，这种滤波方法滤波器可变，提高了干涉条纹的清晰度。设置解缠阈值为0.2，使用最小费用流法(minimum cost flow, MCF)对干涉相位进行相位解缠[16]，最后通过相位转形变和地理编码得到了地理坐标系下的形变信息。

1.2 二维形变提取方法

利用一对升轨数据或者降轨数据提取到的形变信息只能反演单一的雷达视线向形变情况。单纯的升轨或降轨数据获取的沿LOS向的形变量相差可能较大，甚至完全相反，因此单一方向上的形变数据不能有效反映测区地表真实三维形变情况[17]。故本文利用升降轨数据提取到的形变信息进行联合解算。升降轨侧视成像几何如图1所示。利用常规D-InSAR技术提取到的形变是雷达视线方向上的一维形变，是地表东-西向、南-北向以及垂直向形变在视线向上的投影，由成像几何与投影关系可以得到式(2)。

dLOS=dUcosθ-dNcoaαsinθ-dEsinαsinθ.

(2)

式中:dLOS代表视线向形变；dU、dN、dE分别代表地表形变在垂直、南-北和东-西向上的分量；θ代表雷达平均入射角;α代表卫星平均飞行方位角。由式(2)可知要求解地表三维方向上的形变dU,dN和dE需要至少3个不同成像视角的LOS向形变信息，其联合求解式为：

(3)

目前所有的SAR卫星都是以近极地轨道的方式运行的，雷达视线方向对南-北向的形变不敏感，使获取的南-北向形变质量很差[17],且受制于该区域Sentinel-1A数据的成像几何，故文中研究内客忽略南-北向的形变，利用升轨和降轨两个不同成像视角的LOS向形变信息对东-西向和垂直向的形变进行求解，其求解式为：

(4)

(5)

图1 升降轨侧视成像几何图

2 研究区概况及数据选取

2.1 研究区概况

漾濞地震震中距离漾濞县城约7 km,距离大理市约40 km,其震中附近区域地质构造复杂，且发生过多次地质灾害。漾濞地震震中位于川滇菱形块体的西南边界，是维西-乔后-巍山断裂与红河断裂衔接处，震中区域由3条左阶右旋断裂组成，自西向东依次是维西-乔后-巍山断裂、红河断裂北段和红河断裂中段[10]，震中区域高程信息以及主要断裂带分布如图2所示，图2中黑色虚线框为文中同震形变解算的参考区域。

图2 震中主要断裂带分布及影像覆盖范围

2.2 数据选取

文中采用Sentinel-1A卫星TOPSAR(terrain observation with progressive scan,SAR)条带扫描模式获取数据，数据覆盖范围为250 km，分辨率为5 m×20 m。文中实验数据具体参数如表1所示。文中除了采用Sentinel-1A数据外，还采用欧空局提供的DORIS精密轨道数据AUX_POEORB进行轨道误差校正。升降轨影像覆盖范围在图2中用蓝色矩形框表示。

表1 影像信息

3 同震形变分析3.1 视线向同震形变分析

将D-InSAR相位信息转换为视线向的形变数据，得到了降轨和升轨两个视线向形变场，其中降轨表现为两个主要形变区域A和B，A和B大致呈椭圆形状对称分布，A表现为远离视线向的沉降，最大沉降量约为-7 cm，B区域表现为靠近雷达视线向的隆升，最大隆升量约为8 cm，这一结果与王绍俊[10]以及黎朕灵[12]等人的InSAR解算结果基本一致，降轨同震形变场如图3所示，在A与B之间形成了一个近北-西走向的没有明显形变的长轴，如图3中红色虚线所示，有可能为本次地震发震断层所在的空间位置，震中附近已知的维西-巍山断裂位于川滇块体西部边缘，北起雪龙山东麓，止于点苍山西南，走向北北西，但该北-西走向长轴位于维西-巍山断裂以西，不在其主断裂上，杨九元[11]等的研究认为此次地震可能破裂在维西-巍山断裂的隐伏分支断层或一个独立的未知的隐伏主断层上，但是准确的地震断层还需要进一步研究确认。

图3 降轨同震形变场

根据中国地震台网发布的相关数据显示，自2021年5月21日T21：48主震发生后到降轨辅影像成像时间5月22日T23：14止，震区陆续共发生31次余震，其中M3.0以上的余震29次，该时间段内的余震时间序列如图4(a)所示。将余震空间分布叠加到降轨形变场上，如图4(b)所示，不难发现余震发震位置与降轨形变场高度吻合，故文中得到的InSAR降轨形变场包含上述余震的贡献。

图4 余震时间分布及空间分布

为了更好地分析降轨同震形变场，绘制了aa*、bb*、cc*3条剖面线，其分别位于沉降区、隆升区、跨越沉降和隆升区。3条剖面线的位置已在图3中标明，其中剖面线aa*长度为12.7 km，剖面线bb*长度为10.3 km，剖面线cc*长度为12.5 km。沿3条剖面线的形变情况如图5所示。由剖面线形变曲线可知:剖面线aa*在2.5～4.5 km中形成了1个约5 cm的大型沉降漏斗，而且在4.3 km处形变量达到了-5.6 cm。在7.6～8.7 km中形成了一个-4 cm左右的小型沉降漏斗，在8.1 km处沉降量达到了-4.3 cm；剖面线bb*在1～6.8 km内隆升量均超过7 cm，隆升量较大，其余部分隆升量均小于7 cm；剖面线cc*在6 km处为北-西走向长轴所在位置，在6 km之前表现为沉降，沉降量在-0.3～-3.8 cm，6 km之后表现为隆升，隆升量在0～8 cm。

升轨同震形变场因受山区独特地理环境以及卫星成像几何的影响，噪声污染严重，形变特征不明显，但不难发现A和B两个形变区域明显扩大，且两个区域的形变方向与降轨结果相反，其中A区域表现为隆升，最大形变量为7 cm，B区域表现为沉降，最大形变量为-7 cm，如图6所示，相同区域表现出不同的形变方向，与升降轨不同的成像视角有关。自降轨辅影像成像时间2021年5月22日T23：14起至升轨辅影像成像时间6月1日止，震中附近区域发生多次余震，其中M3.0以上共7次(据中国地震台网)，故升轨解算结果应包含了上述余震贡献量，形变区域明显扩大与余震形变累积有关。

图5 剖面线形变曲线

图7 同震库仑应力变化

图6 升轨同震形变场

3.2 同震库仑应力变化计算

漾濞地震发生在维西-巍山断裂和红河断裂交会处，该区域的地震背景活动性较强，处于背景地震预测模型的低值区，地震风险较高[10]。本次地震的发震断层位于维西-巍山断裂以西，且周围区域断层分布密集，主要有红河断裂北段和红河断裂中段，且这些断层均为右旋走滑断层，因此有必要计算本次地震引起的同震库仑应力变化，对周围断层的地震危险性做出评估。

地震发生后，USGS给出了此次地震的震源机制解：走向135°、倾角82°、滑动角-165°；魏云杰等[15]给出了此次地震的平均同震滑动量0.19 m，发震断层以右旋走滑为主。选择上述参数作为同震库仑应力计算的输入参数，设置摩擦系数为0.4(走滑型断层的典型值)。基于弹性半空间位错模型，利用Coulomb3.3软件分别计算了深度为7.5 km和10 km的同震库仑应力变化，如图7所示。

同震库仑应力变化计算结果显示，此次地震对周边断层主要表现为库仑应力卸载作用，仅在断层东北部和东南部表现为加载作用，东北部对应为红河断裂北段，东南部没有已知断裂分布，因此此次地震引起的同震库仑应力变化对红河断裂北段具有加载作用，但是没有达到触发值(0.1 bar)。尽管该结果可以为地震危险性初步评估提供一定的参考，但是精确的震源参数以及同震滑动分布可以为同震库仑应力计算提供更可靠的数据源，进而可以实现对周围断裂发震危险性的精确评估，因此精确的GNSS数据和多维度的同震形变场可以为此提供更好的约束条件。

3.3 东-西向和垂直向同震形变分析

结合升降轨同震形变场解算结果，利用式(5)解算得到了东-西向和垂直向同震形变场，规定东-西向同震形变场形变量正值表示向东运动，负值表示向西运动，垂直向同震形变场形变量正值表示垂直方向的隆升，负值表示垂直方向的沉降。两个方向的形变场呈现出C和D两个主要形变区域，其中C区域面积大约为60 km2，D区域面积大约为40 km2，如图8所示。

图8 垂直向和东-西向同震形变场

图8(a)显示，C区域表现为向东运动，最大运动量达到了5.7 cm，D区域表现为向西运动，最大运动量达到了-7 cm。图8(b)显示，C区域表现为隆升，最大隆升量为3.9 cm，D区域表现为沉降，最大沉降量为-8 cm。

4 结束语

文中利用Sentinel-1A数据，采用二轨差分D-InSAR技术解算得到了2021年5月21日云南漾濞地震升轨和降轨的视线向同震形变场，并计算了此次地震引起的同震库仑应力变化。最后结合升降轨同震形变场直接解算得到了本次地震的东-西向和垂直向同震形变场，可以对进一步开展此次地震震源参数和同震滑动模型精细化反演提供更可靠的约束条件，主要形成了以下几点结论：①视线向同震形变场中，降轨的最大隆升量为8 cm，最大沉降量为-7 cm。升轨的最大隆升量为7 cm，最大沉降量为7 cm。②同震库仑应力变化计算结果显示此次地震仅对地震断层东北部的红河断裂北段有加载作用，但加载量较小，发震危险性较低。③东-西向和垂直向同震形变场呈现两个主要形变区域C和D，形变区域C整体表现为垂直向向上隆升，水平方向向东运动，水平方向最大形变量为5.7 cm，垂直向最大隆升量为3.9 cm。形变区域D整体表现为垂直向向下沉降，水平方向向西运动，水平方向最大形变量为-7 cm，垂直向最大沉降量为-8 cm。