Coseismic Deformation Analysis of Qinghai Maduo Ms 7.4 Earthquake Based on Optical Remote Sensing and SAR Images
-
摘要:
2021-05-22青海省果洛藏族自治州玛多县发生Ms 7.4地震,作为近年来少有的发生在巴颜喀拉块体内部的强震,研究其同震形变场特征是十分必要的。收集了玛多地震前后Sentinel-2和Landsat 8影像,利用光学像素追踪(pixel offset tracking,POT)技术获得了该地震东西向和南北向形变;基于地震前后Sentinel-1升、降轨影像,利用合成孔径雷达干涉测量技术获取了该地震升、降轨雷达视线向(line of sight,LOS)形变,利用合成孔径雷达(synthetic aperture radar,SAR)POT技术获得该地震距离向和方位向形变,联合解算得到其三维同震形变场,并对提取的同震形变场结果进行对比交叉验证。结果表明,此次玛多地震为左旋走滑型地震,同震形变以东西向水平运动为主,发震断裂为江错断裂。基于光学遥感影像,得到该地震东西向和南北向形变大约集中在±1.60 m和±0.60 m;基于SAR影像,得到升轨最大LOS向抬升和沉降量分别约为1.29 m和-1.12 m,降轨最大LOS向抬升和沉降量约分别为1.15 m和-1.26 m;解算的三维同震形变场中,东西向形变约为-2.00~1.70 m,南北向形变主要集中在-1.00~0.50 m,垂直向上沿断裂带两侧呈升降交替运动,形变约在±0.3 m。地震北侧形变量级相较于南侧更大,得到的地表破裂带长约176 km,在东南末端(34.48°N,99.04°E)、西北末端(34.76°N,97.61°E)和西段鄂陵湖南侧(34.74°N,97.75°E)存在分支破裂。基于光学遥感和SAR影像提取的玛多地震同震形变场具有一致特征,且多平台、多技术为获取该地震完整同震形变场补充了更多地表破裂带分支等特征,该研究为光学遥感和SAR影像在地震同震形变监测中的应用提供一些参考。
Abstract:ObjectivesAn Ms 7.4 earthquake occurred in Maduo, Qinghai Province,China on May 22, 2021, which is a rare powerful earthquake that occurred in Bayan Har block in recent years, it is essential to study the characteristics of its coseismic deformation field.
MethodsBased on Sentinel-2 and Landsat 8 images before and after Maduo earthquake, east-west and south-north two-dimensional coseismic deformation fields were obtained by optical pixel offset tracking (POT). Interferometric synthetic aperture radar (InSAR)was utilized to obtain line of sight (LOS) coseismic deformation based on Sentinel-1 ascending and descending images, and synthetic aperture radar (SAR) POT was utilized to obtain range and azimuth directions coseismic deformation of the earthquake, meanwhile, 3D coseismic deformation field of this earthquake is calculated, moreover, the results of coseismic deformation field extracted are compared and verified.
ResultsThe experimental results show that coseismic deformation field characteristics of Maduo earthquake based on optical and SAR images are great consistent. Maduo earthquake is a typical left-lateral strike-slip event, the coseismic deformation is dominated by east-west horizontal movements, and the seismogenic fault is Jiangcuo fault. Based on optical images, it is obtained that east-west and north-south deformations of the earthquake are about ±1.60 m and ±0.60 m,respectively. Based on SAR images, the maximum LOS uplift and subsidence of ascending are about 1.29 m and -1.12 m, and descending are about 1.15 m and -1.26 m. In the 3D coseismic deformation field, the east-west deformation is about -2.00 m to 1.70 m, north-south deformation is mainly concentrated in -1.00 m to 0.50 m, and it moves alternately up and down along both sides of the fault zone in vertical direction, the deformation is about ±0.3 m. The magnitude of deformation values on north side of the earthquake is larger than that on south side. Branch ruptures are shown at the end of the southeast (34.48°N, 99.04°E) and northwest (34.76°N, 97.61°E) of the surface rupture, and south side of Eling Lake (34.74°N, 97.75°E), and aftershock sequence is distributed near these branch ruptures.
ConclusionsInSAR and POT technologies complement each other, which provide an effective way to obtain the complete coseismic deformation field, and multi-platform, high-resolution optical and SAR images provide important datasets.
-
Keywords:
- Maduo earthquake /
- optical remote sensing image /
- SAR image /
- pixel offset tracking(POT) /
- InSAR
-
2021-05-22T02:04:00,中国青海省果洛藏族自治州玛多县发生Ms 7.4地震,震中(34.59°N,98.34°E)位于巴颜喀拉块体北部,是汶川地震之后中国发生震级最高的一次地震(http://www.ceic.ac.cn/)。地震常常造成严重的建筑物破坏及人员伤亡,据报道,此次地震造成青海省果洛藏族自治州7个县、40个乡镇35 521人受灾,同时部分房屋受损,交通干道变形塌陷,严重威胁人民群众的生命和财产安全。对地震灾情信息的快速获取是做好抗震救灾工作的关键,而地震同震形变场监测是其中必不可少的一环,地震同震形变场的获取是研究地震发生机制、地震参数反演、确定发震断层几何特征和运动性质的必要基础,对于揭示研究区域构造变形强度机制,防止次生地质灾害发生具有重要意义。
遥感技术的快速发展使其成为震害监测十分可靠、有效的手段。相较于传统水准测量和全球导航卫星系统(global navigation satellite system,GNSS)技术对离散点形变监测,合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)技术能够全天时、全天候、高精度、大范围地获取连续面的地表形变,对监测大面积受灾区域具有巨大优势[1],在地震同震形变监测中发挥着无可替代的作用。但是InSAR技术最大只能监测出相邻像元间低于半波长的形变[2],在地震近断裂带区域形变量超出最大形变梯度时,容易产生失相干现象,无法获取准确的形变量。针对该问题,许多学者利用基于影像幅度信息的像素追踪(pixel offset tracking,POT)技术来弥补InSAR技术这一缺陷[3],如文献[4]利用合成孔径雷达(synthetic aperture radar,SAR)POT技术得到了2002年Denali地震的同震形变;文献[5]基于EnviSat卫星升、降轨数据,利用InSAR技术和SAR POT技术获取了Bam地震的完整地表形变场;而基于光学遥感影像的POT技术也在Izmit地震[6]、Rudbar地震[7]、Kaikoura地震[8]取得了可靠的同震形变结果。事实证明,POT技术能够在InSAR技术严重失相干区域获取准确的形变信息,在形变量大且具有突变性的地震、滑坡、冰川等地质灾害监测中是一种非常有效的方法[9]。该技术获取形变的精度主要取决于影像的空间分辨率[10],如今高空间分辨率SAR卫星和光学遥感卫星的发展为其提供了重要数据保障。
对于2021年青海玛多地震,许多学者利用SAR影像和GNSS数据对该地震的同震形变及发震构造进行了研究[11-12],而基于光学遥感影像数据的研究较少。本文首先基于青海玛多地震前后Sentinel-2和Landsat 8影像,采用光学POT技术分别获得了该地震东西向和南北向二维形变场,并对形变结果中的误差源进行分析与改正;然后基于地震前后的Sentinel-1升、降轨影像,采用InSAR技术获取了该地震升、降轨视线向(line of sight,LOS)向形变,又采用SARPOT技术得到了该地震距离向和方位向形变;最后联合光学遥感影像和SAR影像解算了玛多地震三维同震形变场,并对形变结果进行交叉对比验证,分析该地震完整地表同震形变场特征,进一步证明光学遥感影像和SAR影像在地震同震形变监测中的可行性和重要性。
1 研究区域与数据
青藏高原位于喜马拉雅山脉以北,是印度板块与欧亚板块向北持续碰撞和挤压的结果,强烈的地质构造活动孕育了柴达木块体、巴颜喀拉块体、羌塘块体和拉萨块体等多个板内次级断块[13]。中国Ms 7.0浅源地震主要集中发生在青藏高原中部巴颜喀拉块体[14]。巴颜喀拉块体在地貌上呈西部狭长、东侧张开的倒三角形态[15],北部东昆仑断裂带和南部甘孜-玉树断裂带均是巨型左旋走滑型活动断裂,东部是右旋走滑逆冲型的龙门山断裂带和岷江断裂带,西部为阿尔金断裂带西南段尾端张性断层,是一条走滑拉张边界[14],块体边界主要为走滑型断裂带,走滑型地震是该区域主体地震,印度板块向北长期强烈俯冲推挤产生的应力应变积累是孕育这些边界断裂带发生强震的重要条件[16]。自1995年以来,巴颜喀拉块体边界发生多次Ms 7.0以上的中强地震,如2001年北部边界发生的昆仑山口地震[17],2011年南部边界发生的玉树地震[18],东部边界发生的2008年汶川地震[19]、2013年庐山地震[20]和2017年九寨沟地震[21],此次玛多Ms 7.4地震发生在巴颜喀拉块体北部,是近年来少有的发生在该块体内部的强震[22](图1)。巴颜喀拉块体北部东昆仑断裂带以南发育了一系列次级活动断裂,由北向南有玛多-甘德断裂、江错断裂、甘德南缘断裂和达日断裂等(图1),构成了该块体由西向东逐渐宽阔的“马尾状”断裂带[23]。众多研究成果表明,玛多地震的发震断裂为江错断裂,地表破裂特征主要体现为一系列雁行状张裂隙、剪切裂缝以及挤压隆起、裂陷等,表现为显著的左旋走滑运动性质[23-24]。玛多地震的发生说明巴颜喀拉块体构造活动仍在持续,作为青藏高原地震活动高发区具有重要研究价值。
本文收集了欧洲空间局(European Space Agency,ESA)Sentinel-2和美国陆地卫星计划Landsat 8两种光学遥感影像。Sentinel-2卫星任务包含Sentinel-2A和Sentinel-2B两颗卫星,可提供空间分辨率分别为10 m、20 m和60 m的光学影像,两颗卫星同时运行可将重访周期由10 d减少为5 d。玛多地震发生前后多为雨雪天气,光学影像质量受天气状况影响较大,因此本文从ESA官网获取了2020-06-24和2020-09-02两景震前Sentinel-2影像、2021-06-09和2021-07-19两景震后Sentinel-2影像,所选影像对时间跨度近1 a,且在主震形变区域云覆盖量较少。相较于Sentinel-2卫星,Landsat 8卫星每16 d实现一次全球覆盖,提供15 m和30 m空间分辨率的光学影像。本文从美国地质调查局(United States Geological Survey,USGS)网站下载了2020-09-17和2021-07-18两景地震前、后Landsat 8影像,影像整体云覆盖量较小,不到2%,所选影像对时间跨度和Sentinel-2影像一致(表1和图1)。另外,Sentinel-2影像10 m分辨率的波段包括可见光波段Band 2、Band 3和Band 4,以及近红外波段Band 8,研究表明Band 8影像在POT算法中结果最好[25]。因此,本文对于Sentinel-2影像选取了Band 8进行POT计算,而Landsat 8卫星选取了较高空间分辨率15 m的全色波段Band 8进行POT计算。
表 1 研究区光学遥感影像信息Table 1. Optical Remote Sensing Images Information Covering the Study Area影像日期 卫星 太阳高度角/(°) 太阳方位角/(°) 含云量/% 拼接域编号 波段(分辨率/m) 距主震天数/d 2020-09-17 Landsat 8 52.63 146.16 0.090 — Band 8 (15) 247 2021-07-02 Landsat 8 67.39 114.92 1.540 — Band 8 (15) 57 19.69 120.22 4.143 T47SLU Band 8 (10) 332 2020-06-24 Sentinel-2B 18.94 122.17 3.560 T47SMU Band 8 (10) 18.19 124.23 2.863 T47SNU Band 8 (10) 31.54 144.67 0.006 T47SLU Band 8 (10) 262 2020-09-02 Sentinel-2B 31.05 146.50 0.008 T47SMU Band 8 (10) 30.57 148.37 0.044 T47SNU Band 8 (10) 19.42 122.69 11.996 T47SLU Band 8 (10) 18 2021-06-09 Sentinel-2B 18.69 124.75 7.231 T47SMU Band 8 (10) 17.97 126.94 6.682 T47SNU Band 8 (10) 22.07 124.05 2.931 T47SLU Band 8 (10) 58 2021-07-19 Sentinel-2B 21.35 125.96 4.106 T47SMU Band 8 (10) 20.64 127.96 7.130 T47SNU Band 8 (10) 本文所用的SAR影像为ESA C波段Sentinel-1卫星影像,该卫星单星重访周期为12 d,Sentinel-1A与Sentinel-1B卫星协同观测重访周期为6 d,但遗憾的是2021-12-23 Sentinel-1B卫星因出现故障无法再提供雷达影像。本文共收集了玛多地震主震前后两对升、降轨干涉宽幅(interferometric wide swath,IW)模式Sentinel-1影像,其中震前为2021-05-20的Sentinel-1A影像,震后为2021-05-26的Sentinel-1B影像(表2和图1)。同时还收集了ESA提供的多普勒轨道学与无线电定位集成卫星(Doppler orbitography and radio-positioning integrated by satellite ,DORIS)精密轨道文件以及航天飞机雷达地形测绘使命(shuttle radar topography mission,SRTM)30 m分辨率的数字高程模型(digital elevation model,DEM)数据。
表 2 研究区SAR影像信息Table 2. SAR Images Information Covering the Study Area影像日期 卫星 升/降轨 轨道编号 垂直基线/m 入射角/(°) 方位角/(°) 时间间隔/d 2021-05-20 Sentinel-1A 升轨 99 53 33.8 -12.9 6 2021-05-26 Sentinel-1B 升轨 2021-05-20 Sentinel-1A 降轨 106 111 39.5 193.0 6 2021-05-26 Sentinel-1B 降轨 2 数据处理方法
2.1 光学遥感影像数据处理
光学遥感影像数据处理使用的是美国加州理工学院在ENVI软件开发的COSI-CORR(the coregistration of optical sensed images and correlation)模块,它可以实现两景光学遥感影像亚像素相关性计算,获取影像事件前后东西向和南北向形变,该模块主要包括统计相关算法和频率相关算法。统计相关算法最大化相关系数的绝对值,算法粗糙但稳健;频率相关算法是基于傅里叶偏移理论在频率域内采用相位相关技术进行匹配计算,相较于统计相关算法更精确[26]。因此,本文采用频率相关算法来提取玛多地震东西向和南北向同震形变场。首先对Sentinel-2影像完成拼接使其覆盖研究区域,然后在COSI-CORR软件频率相关器中[27],将初始滑动窗口和最终滑动窗口均设置为128×128像素,较大的滑动窗口尺寸能够获得更准确的位移,迭代次数设置为2,掩膜阈值设置为0.95,为使两种影像地面分辨率保持一致,最终能够获取相似分辨率的同震形变结果,Sentinel-2和Landsat 8影像移动步长分别设置为9×9像素和6×6像素[25]。
光学POT技术获取地表形变时的误差来源主要包括轨道趋势误差、条带误差、卫星姿态角误差、地形阴影误差和失相干噪声等[28]。针对本文中光学像素偏移量形变结果中的误差源,为保证形变结果精度,需要对各种误差进行后处理。首先对于轨道趋势误差,采用一次多项式曲面拟合模型进行去除[29];对于失相干噪声误差,在数据处理过程中通过设置信噪比阈值(0.95)来掩膜形变场中低于该阈值的区域,从而剔除由于失相干带来的高异常值;在Sentinel-2影像获取的东西向和南北向形变结果中出现明显沿着轨道运行方向均匀分布的线性信号,即条带误差,以及影像拼接所产生的误差,本文均采用传统的均值相减法来去除[30];而对于卫星姿态角误差,采用基于稳定区域的改进的均值相减法进行去除[25,31];然后通过中值滤波对形变结果进行进一步降噪处理,得到地震东西向和南北向二维同震形变场。误差去除前后对比结果如图2所示。由于Landsat 8影像处理过程中没有出现明显的条带误差,对于Landsat 8影像完成了轨道趋势误差、卫星姿态角误差和失相干噪声去除以及中值滤波处理。值得注意的是,光学POT结果存在地形阴影误差,在单一影像对且研究区域地势起伏不大时,该误差可以忽略,在有较多影像对时,可选取多个具有高相似地形阴影误差的影像对实现对该误差的掩膜处理[32]。
2.2 SAR影像数据处理
本文利用GAMMA软件对SAR影像进行数据处理[33],主要包括两部分:一是采用InSAR技术提取地震升、降轨LOS向同震形变场;二是采用POT技术获取地震距离向和方位向二维形变场。
InSAR数据处理中,采用二轨法对地震前后Sentinel-1影像进行差分干涉,首先对影像拼接裁剪使其覆盖研究区域,通过迭代强度匹配和频谱差异方法实现主(震前)从(震后)影像的粗配准和精配准[34],处理过程中对影像进行多视处理(10∶2)实现降噪,利用SRTM DEM和精密轨道数据消除干涉相位中的地形相位和平地效应;然后基于局部条纹频谱对干涉图进行自适应滤波处理[35],减少相位噪声,对低相干区域设置阈值(0.3)进行掩膜,采用最小费用流完成相位解缠[36],利用通用大气校正在线服务(generic atmospheric correction online service for InSAR, GACOS)进行大气校正[37-38],获得地表真实相位;最后将相位值转换成位移值并完成地理编码,得到玛多地震升、降轨LOS向同震形变场。
在进行SAR POT处理时,采用强度信息跟踪法来确定主从影像在距离向和方位向上的偏移量[39]。首先基于InSAR预处理中配准完成的Sentinel-1单视复数影像进行去斜处理[40];然后在GAMMA软件偏移量跟踪程序中设置搜索窗口大小为300×60像素,移动步长为20×4像素,互相干阈值设为0.1[41],考虑结果精度和计算效率,设置过采样因子为2[42],搜索主从影像之间最大互相关系数,计算对应像素之间的偏移量[43],并结合外部DEM数据,采用最小二乘拟合去除由于卫星轨道和成像姿态差异引起的轨道偏移量;最后将得到的地表形变偏移量结果进行中值滤波,地理编码后获得玛多地震距离向和方位向二维同震形变场。
本文在完成上述操作后,采用加权最小二乘解算玛多地震地表三维同震形变[8]。根据SAR影像成像几何关系,建立提取的形变矢量与地表三维形变矢量的关系,以形变场远场区域均方误差作为权重,基于最小二乘准则,解算得到玛多地震地表三维形变的最佳估计值。
3 结果与分析
3.1 光学影像同震形变分析
基于玛多地震前后Sentinel-2和Landsat 8影像,利用光学POT技术,分别得到该地震东西向和南北向形变场(图3)。不论是Sentinel-2影像还是Landsat 8影像,其东西向形变场显示出玛多地震震北向西、震南向东运动,表现为左旋走滑特征,南北向形变场显示出震区震北向北、震南向南运动,同时提取到该地震地表破裂带长约161 km,破裂带东南和西北两侧末端表现出分叉破裂。图3(a)和3(d)对距离地震发生后时间间隔最短的Sentinel-2影像(2020-06-24—2021-06-09)所提取的同震形变区域(黑色虚线)形变值进行统计,发现在地震北侧、距震源中心东西两侧存在两块较大形变区域,向西形变量集中在0.80~1.60 m,这两块形变区域外,该地震北侧向西水平位移量整体约为0.20~0.60 m;在地震南侧,存在西部、中部和东部3块明显形变区域,形变值集中在0.60~1.60 m,该地震南侧向东水平位移量整体约为0.40~1.20 m;而对于南北向形变,其形变量级较小,形变量主体范围在-0.60~0.60 m之间;玛多地震同震形变主要以东西向水平运动为主,且地震北侧相较于南侧形变量级更大。本文选取的两对Sentinel-2影像结果显示,距离震后58天(2020-07-19)相比距离震后18天(2020-06-09)的POT结果形变量明显减小(图3(a)和3(b)、图3(d)和3(e)),这可能和震后的上千次余震相关。
在利用光学POT技术获取玛多地震同震形变时,Sentinel-2影像对(2020-09-02—2021-07-19)和Landsat 8影像对(2020-09-17—2021-07-02)获取时间最为接近,因此基于这两组影像对Sentinel-2和Landsat 8影像的同震形变结果进行对比验证。图3(b)、3(c)、3(e)和3(f)分别是基于这两组影像的光学POT结果,图3(g)和3(h)分别为这两组影像形变场区域(黑色虚线)东西向和南北向形变差值。
对图3中形变值进行统计分析,可见在地震形变区域,Landsat 8影像的POT结果的高异常值多于Sentinel-2影像,主要原因为Sentinel-2影像比Landsat 8影像空间分辨率高;两种影像获取的形变值基本一致,东西向和南北向的形变差值平均值分别为-0.16 m和-0.14 m,标准差分别为0.94 m和0.81 m,在影像失相干区域边界处形变差值较大。另外,对这两种光学遥感影像获得的形变结果稳定区域(图3黑色实线)形变值进行了分析(表3),可以看出东西向形变结果优于南北向,平均值均较小,而标准差显示出Landsat 8影像POT结果相较于Sentinel-2影像更好,表明在本文所选稳定区域Sentinel-2影像的POT结果高异常值较多,经检查发现可能原因是所选区域在Sentinel-2震后影像(表1中2021-07-19影像)云覆盖量(4.106%)较Landsat 8震后影像(表1中2021-07-02影像)云覆盖量(1.540%)更大,导致基于地震前后光学遥感影像的POT结果受云覆盖影响误差较大,从而有更多的高异常值。
表 3 Sentinel-2与Landsat 8影像的POT结果对比验证/mTable 3. Comparison of Pixel Offset Tracking Results Between Sentinel-2 and Landsat 8 Images /m形变方向 Sentinel-2稳定区域 Landsat 8稳定区域 形变区域差值 平均值 标准差 平均值 标准差 平均值 标准差 东西向 -0.03 0.39 -0.05 0.34 -0.16 0.94 南北向 -0.09 0.58 0.06 0.46 -0.14 0.81 3.2 SAR影像同震形变分析
图4为基于Sentinel-1影像(2021-05-20—2021-05-26)的玛多地震同震形变场。图4(a)和4(d)分别为基于Sentinel-1升、降轨影像,利用InSAR技术获得玛多地震LOS向同震形变场,形变场呈长椭圆形,形似花生状,整体为SEE-NWW走向,提取到地表破裂带长约176 km,得到升轨LOS向最大抬升量约为1.29 m,最大沉降量约为-1.12 m,呈现出震区北升南降的形变特征;降轨LOS向最大抬升量约为1.15 m,最大沉降量约为-1.26 m,呈现出震区北降南升的形变特征。升、降轨形变场特征显示玛多地震为左旋走滑型地震,且相比升轨LOS向,降轨LOS向与玛多地震同震形变场走向大致平行,表现出更大的同震形变范围,相较于光学遥感影像提取的地表破裂带特征,降轨形变场另外在西段鄂陵湖南侧(34.74°N,97.75°E)显示出分支破裂。
图4(b)和4(e)、图4(c)和4(f)分别是基于SAR POT技术获得的玛多地震升、降轨距离向和方位向二维同震形变场。可以看到,与InSAR技术提取的LOS向形变场相比,SAR POT技术提取的距离向形变场在地震破裂带附近地表位移较大区域仍能够获取连续形变值,可以作为InSAR结果的有效补充,为获取地震完整地表同震形变场提供参考。然而由于SAR POT技术的精度主要与空间分辨率和地形高差有关[44],导致获取的方位向形变场并不理想,且降轨较之于升轨方位向形变结果显示出更大的误差,这可能是降轨影像垂直基线较大的缘故[45]。
在图4(a)中取4条剖面线(aa'、bb'、cc'、dd')对升、降轨LOS向和距离向形变结果进行分析(图5)。可以看出升轨和降轨形变趋势相反,震中区域(bb′)附近形变最小,距离向和LOS向具有相同形变趋势,由于InSAR和POT两种技术所得结果分辨率和精度的差异,SAR POT距离向结果出现抖动状特征。取4条剖面线上距离最近的LOS向形变点和距离向形变点作对比,可见升、降轨LOS向、距离向形变大小基本相同,地表位移量范围基本一致,得到升轨4条剖面线上两方向形变对比标准差分别为0.13 m、0.09 m、0.31 m和0.21 m,降轨4条剖面线上两方向形变对比标准差分别为0.15 m、0.17 m、0.44 m和0.16 m;并且升、降轨的LOS向和距离向结果显示出在距离a点约23.70 km、距离b点约24.50 km、距离c点约25.60 km、距离d点约22.00 km处地表产生断裂。
图 6 玛多地震三维同震形变场(a-c、e-g、i-k、m-o),以及玛多地震东西向和南北向同震形变对比分析(m、h、l、p)注:图(m-p)中正值表示向东、向北和垂直向抬升运动,负值表示向西、向南和垂直向沉降运动。Figure 6. Three Dimensional Coseismic Deformation Fields of Maduo Earthquake(a-c,e-g,i-k,m-o),and Comparison and Analysis of East-West and South-North Coseismic Deformations of Maduo Earthquake(m,h,l,p)InSAR升、降轨LOS向形变和SAR POT升轨方位向形变,分别采用未加权和加权的最小二乘解算得到的玛多地震地表三维同震形变场,在此基础上,图6(i)~6(k)和6(m)~6(o)分别为添加SAR POT升、降轨距离向形变参与解算得到的该地震地表三维同震形变场,其中Sentinel-1降轨方位向形变结果由于误差较大,未参与解算。图6显示提取到的地表三维形变结果符合该地震左旋走滑运动特征,东西向形变场显示该地震最大西向和东向形变分别位于震源中心东南侧约30 km和60 km处,形变量分别约为-2.00 m和1.70 m;在南北向形变场中,震区形变特征与光学POT南北向结果基本一致,呈现出震北向北、震南向南运动,形变值大约集中在-1.00~0.50 m;提取的地表垂直向形变主要发生在地震近断裂带区域,表现为沿断裂带两侧的升降交替运动,且形变量级较小,集中在-0.30~0.30 m之间。此外,在东西向和垂直向形变场西段鄂陵湖南侧(34.74°N,97.75°E)存在分支破裂。对比已有的野外现场考察的玛多地震地表破裂[23],本文三维形变场结果与震后地表产生裂隙、隆起、裂陷等引起的形变具有一致性特征。
在光学POT处理中,基于COSI-CORR软件获取地表形变的互相关匹配算法精度可以达到1/20个像素[46],理论上基于Sentinel-2和Landsat 8影像的POT结果精度分别可以达到0.50 m和0.75 m。相较于Landsat 8影像,使用精度更高且时隔主震最短的Sentinel-2影像(2020-06-24—2021-06-09)同震形变结果与基于SAR影像三维同震形变场中东西向和南北向形变进行验证(图6(d)、6(h)、6(l)、6(p))。基于升、降轨LOS向形变和升轨方位向形变解算的三维形变场与光学POT东西向和南北向形变差值平均值分别为0.001 m和0.071 m,标准差分别为0.674 m和0.605 m;增加升降轨距离向形变参与解算后,其东西向和南北向形变差值平均值分别为0.002 m和0.078 m,标准差分别为0.669 m和0.604 m,可以看出,增加距离向形变参与解算后,提高了三维形变场中东西向和南北向形变与光学POT形变结果的一致性,且南北向形变一致性程度高于东西向,而形变差值的平均值增大的原因很可能是由于距离向形变结果精度相对较低导致的。
本文基于光学遥感和SAR影像,利用POT技术和InSAR技术得到了玛多地震多个方向同震形变场,提取到该地震地表破裂带如图7所示。图7显示出玛多地震地表破裂主破裂带走向约在N101.17°E至N101.43°E。利用POT技术,基于Sentinel-2影像提取的地表破裂带长约161 km,破裂带东南端和西北端存在分叉,西北端分叉破裂走向约为N116.20°E和N87.28°E,东南端分叉破裂走向约为N85.44°E和N107.99°E;Landsat 8影像未完全覆盖研究区域,但提取到西北端分叉破裂走向(N110.08°E和N85.37°E)与Sentinel-2基本保持一致。此外,基于Sentinel-1影像得到的地表破裂带长约170 km,破裂带东南端和西北端分别约为N87.33°E、N108.84°E走向分叉和N112.10°E、N94.14°E走向分叉;而利用InSAR技术提取到地表破裂带长约176 km,破裂带东南端和西北端分叉破裂走向分别约为N81.90°E和N105.09°E、N109.39°E和N89.23°E,且SAR形变场另外在鄂陵湖南侧显示出分支破裂,走向约为N148.60°E(表4)。最后综合光学遥感和SAR影像提取的玛多地震地表破裂带长约176 km,可以看出该地震地表破裂带在东南末端(34.48°N,99.04°E)和西北末端(34.76°N,97.61°E)形成分叉,在主破裂带西段(鄂陵湖南侧(34.74°N,97.75°E))存在分支破裂。本文收集了玛多地震截至2021-10-01零时的Ms 2.0以上296次余震(图7),余震整体沿地表破裂带呈SEE-NWW向狭长条带状分布,定位地表破裂带和余震序列分布,确定该地震发震断裂为江错断裂,余震序列也可以看出在地表破裂带分支附近的展布。
表 4 玛多地震地表破裂信息表Table 4. Surface Rupture Information Table of Maduo Earthquake地表破裂带 POT InSAR Sentinel-2 Landsat 8 Sentinel-1 Sentinel-1 主破裂带 整体走向 N101.27°E N101.43°E N101.32°E N101.17°E 东南末端(34.48°N、99.04°E) 北侧分支走向 N85.44°E — N87.33°E N81.90°E 南侧分支走向 N107.99°E — N108.84°E N105.09°E 西北末端(34.76°N、97.61°E) 北侧分支走向 N116.20°E N110.08°E N112.10°E N109.39°E 南侧分支走向 N87.28°E N85.37°E N94.14°E N89.23°E 鄂陵湖南侧(34.74°N、97.75°E) 北侧分支走向 — — — N148.60°E 4 结语
本文以2021-05-22青海玛多Ms 7.4地震为研究对象,基于Sentinel-2和Landsat 8两种光学遥感影像,利用光学POT技术得到玛多地震东西向形变集中在±1.60 m,南北向形变集中在±0.60 m;基于Sentinel-1升、降轨SAR影像,利用InSAR技术得到该地震升轨最大LOS向抬升和沉降量分别约为1.29 m和-1.12 m,降轨最大LOS向抬升和沉降量分别约为1.15 m和-1.26 m;联合光学遥感和SAR影像获取的三维同震形变场中,东西向形变约为-2.00~1.70 m,南北向形变主要集中-1.00~0.50 m,垂直向形变约在-0.30~0.30 m之间。综合分析玛多地震同震形变场,确定此次玛多地震为左旋走滑型地震,定位发震断裂为江错断裂;同震形变以东西向水平运动为主,地震北侧形变值量级相较于南侧更大,南北向形变量较小,垂直方向上呈现出沿断裂带两侧升降交替运动,与该地震沿断
感谢ESA提供的Sentinel-2、Sentinel-1影像和DORIS精密轨道文件,USGS提供的Landsat 8影像和SRTM 30 m DEM数据,英国纽卡斯大学和长安大学提供的GACOS数据,以及中国地震台网中心提供的青藏高原历史地震目录、玛多地震余震序列和震源中心数据。http://ch.whu.edu.cn/cn/article/doi/10.13203/j.whugis20220615 -
图 6 玛多地震三维同震形变场(a-c、e-g、i-k、m-o),以及玛多地震东西向和南北向同震形变对比分析(m、h、l、p)
注:图(m-p)中正值表示向东、向北和垂直向抬升运动,负值表示向西、向南和垂直向沉降运动。
Figure 6. Three Dimensional Coseismic Deformation Fields of Maduo Earthquake(a-c,e-g,i-k,m-o),and Comparison and Analysis of East-West and South-North Coseismic Deformations of Maduo Earthquake(m,h,l,p)
表 1 研究区光学遥感影像信息
Table 1 Optical Remote Sensing Images Information Covering the Study Area
影像日期 卫星 太阳高度角/(°) 太阳方位角/(°) 含云量/% 拼接域编号 波段(分辨率/m) 距主震天数/d 2020-09-17 Landsat 8 52.63 146.16 0.090 — Band 8 (15) 247 2021-07-02 Landsat 8 67.39 114.92 1.540 — Band 8 (15) 57 19.69 120.22 4.143 T47SLU Band 8 (10) 332 2020-06-24 Sentinel-2B 18.94 122.17 3.560 T47SMU Band 8 (10) 18.19 124.23 2.863 T47SNU Band 8 (10) 31.54 144.67 0.006 T47SLU Band 8 (10) 262 2020-09-02 Sentinel-2B 31.05 146.50 0.008 T47SMU Band 8 (10) 30.57 148.37 0.044 T47SNU Band 8 (10) 19.42 122.69 11.996 T47SLU Band 8 (10) 18 2021-06-09 Sentinel-2B 18.69 124.75 7.231 T47SMU Band 8 (10) 17.97 126.94 6.682 T47SNU Band 8 (10) 22.07 124.05 2.931 T47SLU Band 8 (10) 58 2021-07-19 Sentinel-2B 21.35 125.96 4.106 T47SMU Band 8 (10) 20.64 127.96 7.130 T47SNU Band 8 (10) 表 2 研究区SAR影像信息
Table 2 SAR Images Information Covering the Study Area
影像日期 卫星 升/降轨 轨道编号 垂直基线/m 入射角/(°) 方位角/(°) 时间间隔/d 2021-05-20 Sentinel-1A 升轨 99 53 33.8 -12.9 6 2021-05-26 Sentinel-1B 升轨 2021-05-20 Sentinel-1A 降轨 106 111 39.5 193.0 6 2021-05-26 Sentinel-1B 降轨 表 3 Sentinel-2与Landsat 8影像的POT结果对比验证/m
Table 3 Comparison of Pixel Offset Tracking Results Between Sentinel-2 and Landsat 8 Images /m
形变方向 Sentinel-2稳定区域 Landsat 8稳定区域 形变区域差值 平均值 标准差 平均值 标准差 平均值 标准差 东西向 -0.03 0.39 -0.05 0.34 -0.16 0.94 南北向 -0.09 0.58 0.06 0.46 -0.14 0.81 表 4 玛多地震地表破裂信息表
Table 4 Surface Rupture Information Table of Maduo Earthquake
地表破裂带 POT InSAR Sentinel-2 Landsat 8 Sentinel-1 Sentinel-1 主破裂带 整体走向 N101.27°E N101.43°E N101.32°E N101.17°E 东南末端(34.48°N、99.04°E) 北侧分支走向 N85.44°E — N87.33°E N81.90°E 南侧分支走向 N107.99°E — N108.84°E N105.09°E 西北末端(34.76°N、97.61°E) 北侧分支走向 N116.20°E N110.08°E N112.10°E N109.39°E 南侧分支走向 N87.28°E N85.37°E N94.14°E N89.23°E 鄂陵湖南侧(34.74°N、97.75°E) 北侧分支走向 — — — N148.60°E -
[1] ZHANG J M, ZHU W, CHENG Y Q, et al. Landslide Detection in the Linzhi–Ya’an Section Along the Sichuan-Tibet Railway Based on InSAR and Hot Spot Analysis Methods[J]. Remote Sensing, 2021, 13(18): 3566.
[2] MASSONNET D, FEIGL K L. Radar Interferometry and Its Application to Changes in the Earth’s Surface[J]. Reviews of Geophysics, 1998, 36(4): 441-500.
[3] MICHEL R, AVOUAC J P, TABOURY J. Measuring Ground Displacements from SAR Amplitude Images: Application to the Landers Earthquake[J]. Geophysical Research Letters, 1999, 26(7): 875-878.
[4] ELLIOTT J L, FREYMUELLER J T, RABUS B. Coseismic Deformation of the 2002 Denali Fault Earthquake: Contributions from Synthetic Aperture Radar Range Offsets[J]. Journal of Geophysical Research: Solid Earth, 2007, 112(B6): B06421.
[5] FIALKO Y, SANDWELL D, SIMONS M, et al. Three-Dimensional Deformation Caused by the Bam, Iran, Earthquake and the Origin of Shallow Slip Deficit[J]. Nature, 2005, 435(7040): 295-299.
[6] MICHEL R, AVOUAC J P. Deformation Due to the 17 August 1999 Izmit, Turkey, Earthquake Measured from SPOT Images[J]. Journal of Geophysical Research: Solid Earth, 2002, 107(B4): ETG2-1-ETG2-6.
[7] AJORLOU N, HOLLINGSWORTH J, MOUSAVI Z, et al. Characterizing Near-field Surface Deformation in the 1990 Rudbar Earthquake (Iran) Using Optical Image Correlation[J]. Geochemistry, Geophysics, Geosystems, 2021, 22(6): e2021GC009704.
[8] 王乐洋, 邹阿健. 结合SAR和光学数据检索凯库拉Mw 7.8地震三维形变[J]. 武汉大学学报(信息科学版), 2024, 49(2): 303-312. WANG Leyang, ZOU Ajian. Retrieving 3D Coseismic Deformation of 2016 Mw 7.8 Kaikoura Earthquake Using SAR and Optical Data[J]. Geomatics and Information Science of Wuhan University, 2024, 49(2): 303-312.
[9] ZHANG C L, LI Z H, YU C, et al. An Integrated Framework for Wide-area Active Landslide Detection with InSAR Observations and SAR Pixel Offsets[J]. Landslides, 2022, 19(12): 2905-2923.
[10] YAO J M, YAO X, WU Z Q, et al. Research on Surface Deformation of Ordos Coal Mining Area by Integrating Multitemporal D-InSAR and Offset Tracking Technology[J]. Journal of Sensors, 2021, 2021(1): 6660922.
[11] 余鹏飞, 熊维, 陈威, 等. 基于GNSS和InSAR约束的2021年玛多MS7.4地震同震滑动分布及应用[J]. 地球物理学报, 2022, 65(2): 509-522. YU Pengfei, XIONG Wei, CHEN Wei, et al. Slip Model of the 2021 Ms 7.4 Madoi Earthquake Constrained by GNSS and InSAR Coseismic Deformation[J]. Chinese Journal of Geophysics, 2022, 65(2): 509-522.
[12] 刘计洪, 胡俊, 李志伟, 等. 融合哨兵1号和ALOS-2数据的2021年青海玛多地震高精度三维同震形变场研究[J]. 中国科学: 地球科学, 2022, 52(5): 882-892. LIU Jihong, HU Jun, LI Zhiwei, et al. Study on High-precision Three-dimensional Coseismic Deformation Field of Maduo Earthquake in Qinghai Province in 2021 by Fusing Sentinel 1 and ALOS-2 Data[J]. Scientia Sinica (Terrae), 2022, 52(5): 882-892.
[13] 李海兵, 潘家伟, 孙知明, 等. 大陆构造变形与地震活动: 以青藏高原为例[J]. 地质学报, 2021, 95(1): 194-213. Li Haibing, Pan Jiawei, Sun Zhiming, et al. Continental Tectonic Deformation and Seismic Activity: A Case Study from the Tibetan Plateau[J]. Acta Geologica Sinica, 2021, 95(1): 194-213.
[14] 邓起东, 程绍平, 马冀, 等. 青藏高原地震活动特征及当前地震活动形势[J]. 地球物理学报, 2014, 57(7): 2025-2042. DENG Qidong, CHENG Shaoping, MA Ji, et al. Seismic Activities and Earthquake Potential in the Tibetan Plateau[J]. Chinese Journal of Geophysics, 2014, 57(7): 2025-2042.
[15] 嘉世旭, 林吉焱, 郭文斌, 等. 巴颜喀拉块体地壳结构多样性探测[J]. 地球物理学报, 2017, 60(6): 2226-2238. JIA Shixu, LIN Jiyan, GUO Wenbin, et al. Investigation on Diversity of Crustal Structures Beneath the Bayan Har Block[J]. Chinese Journal of Geophysics, 2017, 60(6): 2226-2238.
[16] 高翔, 邓起东. 巴颜喀喇断块边界断裂强震活动分析[J]. 地质学报, 2013, 87(1): 9-19. GAO Xiang, DENG Qidong. Activity Analysis of Large Earthquakes in Boundary Faults Around the Bayankala Faulting Block[J]. Acta Geologica Sinica, 2013, 87(1): 9-19.
[17] ROBINSON D P, BROUGH C, DAS S. The Mw 7.8, 2001 Kunlunshan Earthquake: Extreme Rupture Speed Variability and Effect of Fault Geometry[J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B8): B08303.
[18] LI Z H, ELLIOTT J R, FENG W P, et al. The 2010 MW 6.8 Yushu (Qinghai, China) Earthquake: Constraints Provided by InSAR and Body Wave Seismology[J]. Journal of Geophysical Research: Solid Earth, 2011, 116(B10): B10302.
[19] FENG G C, HETLAND E A, DING X L, et al. Coseismic Fault Slip of the 2008 Mw 7.9 Wenchuan Earthquake Estimated from InSAR and GPS Measurements[J]. Geophysical Research Letters, 2010, 37(1): L01302.
[20] HUANG Y, QIAO X J, FREYMUELLER J T, et al. Fault Geometry and Slip Distribution of the 2013 Mw 6.6 Lushan Earthquake in China Constrained by GPS, InSAR, Leveling, and Strong Motion Data[J]. Journal of Geophysical Research: Solid Earth, 2019, 124(7): 7341-7353.
[21] ZHAO D Z, QU C Y, SHAN X J, et al. InSAR and GPS Derived Coseismic Deformation and Fault Model of the 2017 Ms7.0 Jiuzhaigou Earthquake in the Northeast Bayanhar Block[J]. Tectonophysics, 2018, 726: 86-99.
[22] 李志才, 丁开华, 张鹏, 等. GNSS观测的2021年青海玛多地震(Mw 7.4)同震形变及其滑动分布[J]. 武汉大学学报(信息科学版), 2021, 46(10): 1489-1497. LI Zhicai, DING Kaihua, ZHANG Peng, et al. Coseismic Deformation and Slip Distribution of 2021 Mw 7.4 Madoi Earthquake from GNSS Observation[J]. Geomatics and Information Science of Wuhan University, 2021, 46(10): 1489-1497.
[23] 潘家伟, 白明坤, 李超, 等. 2021年5月22日青海玛多Ms7.4地震地表破裂带及发震构造[J]. 地质学报, 2021, 95(6): 1655-1670. PAN Jiawei, BAI Mingkun, LI Chao, et al. Coseismic Surface Rupture and Seismogenic Structure of the 2021-05-22 Maduo(Qinghai) Ms7.4 Earthquake[J]. Acta Geologica Sinica, 2021, 95(6): 1655-1670.
[24] 王未来, 房立华, 吴建平, 等. 2021年青海玛多MS7.4地震序列精定位研究[J]. 中国科学: 地球科学, 2021, 51(7): 1193-1202. WANG Weilai, FANG Lihua, WU Jianping, et al. Study on Precise Location of Maduo MS7.4 Earthquake Sequence in Qinghai in 2021[J]. Scientia Sinica (Terrae), 2021, 51(7): 1193-1202.
[25] 贺礼家, 冯光财, 冯志雄, 等. 哨兵-2号光学影像地表形变监测: 以2016年Mw 7.8新西兰凯库拉地震为例[J]. 测绘学报, 2019, 48(3): 339-351. HE Lijia, FENG Guangcai, FENG Zhixiong, et al. Coseismic Displacements of 2016 Mw7.8 Kaikoura, New Zealand Earthquake, Using Sentinel-2 Optical Images[J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(3): 339-351.
[26] LEPRINCE S, BARBOT S, AYOUB F, et al. Automatic and Precise Orthorectification, Coregistration, and Subpixel Correlation of Satellite Images, Application to Ground Deformation Measurements[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(6): 1529-1558.
[27] YANG W T, WANG Y J, WANG Y Q, et al. Retrospective Deformation of the Baige Landslide Using Optical Remote Sensing Images[J]. Landslides, 2020, 17(3): 659-668.
[28] 柳林, 宋豪峰, 杜亚男, 等. 联合哨兵2号和Landsat 8估计白格滑坡时序偏移量[J]. 武汉大学学报(信息科学版), 2021, 46(10): 1461-1470. LIU Lin, SONG Haofeng, DU Yanan, et al. Time-series Offset Tracking of the Baige Landslide Based on Sentinel-2 and Landsat 8[J]. Geomatics and Information Science of Wuhan University, 2021, 46(10): 1461-1470.
[29] FENG G C, DING X L, LI Z W, et al. Calibration of an InSAR-derived Coseimic Deformation Map Associated with the 2011 Mw-9.0 Tohoku-Oki Earthquake[J]. IEEE Geoscience and Remote Sensing Letters, 2012, 9(2): 302-306.
[30] MICHEL R, AVOUAC J P. Coseismic Surface Deformation from Air Photos: The Kickapoo Step over in the 1992 Landers Rupture[J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B3): B03408.
[31] SCHERLER D, LEPRINCE S, STRECKER M R. Glacier-surface Velocities in Alpine Terrain from Optical Satellite Imagery—Accuracy Improvement and Quality Assessment[J]. Remote Sensing of Environment, 2008, 112(10): 3806-3819.
[32] DING C, FENG G C, LI Z W, et al. Spatio-temporal Error Sources Analysis and Accuracy Improvement in Landsat 8 Image Ground Displacement Measurements[J]. Remote Sensing, 2016, 8(11): 937.
[33] WEGMÜLLER U, WERNER C. Gamma SAR Processor and Interferometry Software[C]//Space at the service of Our Environment,Florence,1997.
[34] SCHEIBER R, MOREIRA A. Coregistration of Interferometric SAR Images Using Spectral Diversity[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(5): 2179-2191.
[35] GOLDSTEIN R M, WERNER C L. Radar Interferogram Filtering for Geophysical Applications[J]. Geophysical Research Letters, 1998, 25(21): 4035-4038.
[36] PEPE A, LANARI R. On the Extension of the Minimum Cost Flow Algorithm for Phase Unwrapping of Multitemporal Differential SAR Interferograms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(9): 2374-2383.
[37] YU C, PENNA N T, LI Z H. Generation of Real-time Mode High-resolution Water Vapor Fields from GPS Observations[J]. Journal of Geophysical Research: Atmospheres, 2017, 122(3): 2008-2025.
[38] YU C, LI Z H, PENNA N T. Interferometric Synthetic Aperture Radar Atmospheric Correction Using a GPS-based Iterative Tropospheric Decomposition Model[J]. Remote Sensing of Environment, 2018, 204: 109-121.
[39] STROZZI T, LUCKMAN A, MURRAY T, et al. Glacier Motion Estimation Using SAR Offset-tracking Procedures[J]. IEEE Transactions on Geoscience and Remote Sensing, 2002, 40(11): 2384-2391.
[40] WEGNÜLLER U, WERNER C, STROZZI T, et al. Sentinel-1 Support in the GAMMA Software[J]. Procedia Computer Science, 2016, 100: 1305-1312.
[41] 张成龙, 李振洪, 张双成, 等. 综合遥感解译2022年Mw 6.7青海门源地震地表破裂带[J]. 武汉大学学报(信息科学版), 2022, 47(8): 1257-1270. ZHANG Chenglong, LI Zhenhong, ZHANG Shuangcheng, et al. Surface Ruptures of the 2022 Mw 6.7 Menyuan Earthquake Revealed by Integrated Remote Sensing[J]. Geomatics and Information Science of Wuhan University, 2022, 47(8): 1257-1270.
[42] WERNER C, WEGMULLER U, STROZZI T, et al. Precision Estimation of Local Offsets Between Pairs of SAR SLCS and Detected SAR Images[C]//Proceedings of 2005 IEEE International Geoscience and Remote Sensing Symposium, 2005. IGARSS '05. Seoul, 2005: 4803-4805.
[43] HUANG L, LI Z. Comparison of SAR and Optical Data in Deriving Glacier Velocity with Feature Tracking[J]. International Journal of Remote Sensing, 2011, 32(10): 2681-2698.
[44] Yan Shiyong, GUO Huadong, LIU Guang, et al. Mountain Glacier Displacement Estimation Using a DEM-Assisted Offset Tracking Method with ALOS/PALSAR Data[J]. Remote Sensing Letters, 2013, 4(5): 494-503.
[45] ZHANG W T, ZHU W, TIAN X D, et al. Improved DEM Reconstruction Method Based on Multibaseline InSAR[J]. IEEE Geoscience and Remote Sensing Letters, 2021, 19: 4011505.
[46] KONCA A O, LEPRINCE S, AVOUAC J P, et al. Rupture Process of the 1999 Mw 7.1 Duzce Earthquake from Joint Analysis of SPOT, GPS, InSAR, Strong-Motion, and Teleseismic Data: A Supershear Rupture with Variable Rupture Velocity[J]. The Bulletin of the Seismological Society of America, 2010, 100(1): 267-288.