2024, Vol. 44
2. 自然资源部生态地质与灾害防控重点实验室，西安市雁塔路126号，710054;
3. 黄土科学全国重点实验室，西安市雁塔路126号，710054
据中国地震台网中心测定，2022-01-08青海省门源县发生6.9级地震，震源深度10 km，震中37.77°N、101.26°E，是一次左旋走滑型地震[1]。此次地震位于青藏高原东北部，地形复杂，地震监测能力较弱，传统地震形变监测手段，如GNSS和水准测量难以获得全面的形变信息。合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)因具有空间分辨率高、覆盖范围广、全天时、全天候、精度高等优势，被广泛应用于同震形变场特征观测和构造机理研究[2]。已有研究基于Sentinel-1数据获取门源地震视线向(LOS)一维形变并对其进行反演，而在地震形变监测中，描述地表形变和发震断层几何学特征最准确的方式是获取真实的同震三维形变场[3]。Liu等[4]研究获取2022年门源地震三维同震形变，但未对该事件的断层模型进行解算。而在地震研究中，需利用InSAR同震形变场数据，结合弹性半空间均匀滑移模型进行初步非线性优化[5]，以定义断层破裂机制和几何结构。再以此为先验信息，使用分布式滑动模型开展受三维同震形变场约束的线性反演[6]，获得发震断层的精细滑动分布。
本文运用多孔径InSAR(MAI)技术[7]和像素偏移量追踪(POT)技术[8]获取方位向(azimuth，AZI)形变，再联合升降轨LOS向形变位移信息，采用SM-VCE方法[3]获取地表三维形变结果，并反演断层倾角、走向等几何参数以及同震滑动分布特征。该项研究对了解地块断裂带运动特征与整体地震危险性评估具有重要意义。
1 研究区地质构造背景与数据源 1.1 研究区概况2022年门源MW6.6地震发生在青藏高原东北缘的祁连山地块内部，是我国地震活动与构造变形非常显著的地区之一。该地区构造格局十分复杂，发育有长约1 000 km的祁连-海原断裂带。在该地震动力学背景下，需要进一步关注祁连山地块及周边地区的强震风险。
冷龙岭断裂整体走向约110°，位于门源盆地北部，具有明显的线性特征，全长约120 km。断裂呈左旋走滑特征，为全新世活动断裂[9]。托莱山断裂走向约115°，全长约300 km，在该断裂最东端15 km范围内，走向呈东西向，约89.1°[10]。表 1为USGS和GCMT提供的详细震源参数，图 1为该区域构造背景。
|
表 1 2022年门源地震震源机制解 Tab. 1 The focal mechanism solution of the 2022 Menyuan earthquake |
|
图 1 研究区域构造背景 Fig. 1 Tectonic setting of the study area |
收集覆盖研究区的C波段Sentinel-1升降轨遥感影像以及L波段ALOS-2降轨影像。L波段相较于C波段波长更长，因此具有更强的穿透能力，能更有效地抑制去相关噪声。具体影像参数见表 2。
|
|
表 2 卫星影像参数 Tab. 2 Satellite image parameters |
本文利用InSAR技术获取LOS向和方位向同震形变，并采用SM-VCE算法构建门源MW6.6地震三维同震位移。再以上述数据为约束，利用两步法反演得到发震断层最优几何参数与精细断层滑动分布，具体流程如图 2所示。
|
图 2 技术流程图 Fig. 2 Technology flowchart |
使用GAMMA软件平台对Sentinel-1卫星Track128(升轨)和Track33(降轨)数据进行D-InSAR处理[2]。采用Goldstein自适应滤波方法去除噪声相位[11]，滤波窗口设置为64×64；对影像进行距离向: 方位向为10 ∶2的多视处理以提高信噪比；采用最小费用流法(MCF)进行解缠[12]，设定相干性阈值为0.55；应用通用型卫星雷达大气改正系统(GACOS)减小大气延迟的影响[13]；将相位值转变为形变值并进行地理编码。使用D-InSAR技术处理ALOS-2数据时，设置多视比为3 ∶10，其余处理方式与Sentinel-1数据一致，得到WGS-84地理坐标系下门源MW6.6地震升降轨干涉条纹图(图 3)。
|
图 3 2022年门源地震同震干涉条纹图 Fig. 3 Coseismic interference fringe of the 2022 Menyuan earthquake |
整体来看，升降轨形变区域表现出较好的相干性。升降轨同震形变场沿冷龙岭断裂呈对称分布，地震震中位于呈蝴蝶状分布的条纹图中心，上下两侧断层破裂迹线明显，周围干涉条纹分布非常密集。结合定量分析可知，Sentinel-1升轨数据LOS向最大抬升形变量约为0.53 m，最大沉降形变量约为0.58 m；Sentinel-1降轨数据LOS向最大抬升形变量约为0.74 m，最大沉降形变量约为0.6 m。ALOS-2降轨数据LOS向最大抬升形变量约为1.04 m，最大沉降形变量约为0.71 m。根据雷达卫星成像特点，可以确定发震断层活动特征与左旋走滑断层运动特征一致。
为了更好地了解门源地震同震形变场的空间分布特征，从Sentinel-1卫星升降轨同震形变场和ALOS-2卫星降轨形变场中提取未经掩膜处理的跨断层剖线i-i′和ii-ii′(剖线位置见图 3)，获取剖线对应的地形信息和形变信息。由图 4可知，在升降轨形变剖线中存在明显的上升和沉降中心，且近场形变梯度较大；降轨与升轨形变剖线具有相反的上升和沉降趋势。
|
图 4 LOS向同震形变场剖面 Fig. 4 Coseismic deformation field profile in LOS direction |
D-InSAR方法仅能获得卫星LOS向一维形变信息，难以有效反映真实的三维地表形变信息，尤其是南北向分量[3]。相比之下，MAI和POT方法都可以获取卫星方位向同震形变场，从而弥补仅有卫星LOS向观测值的局限。使用POT技术时，利用D-InSAR预处理中配准完成的Sentinel-1影像，在偏移量追踪程序中设置搜索窗口大小为128像素×64像素，设置互相干阈值为0.1，过采样因子为2，采用最小二乘拟合去除轨道偏移量。最后对地表形变偏移量结果进行中值滤波，通过地理编码后获得门源地震距离向和方位向的二维同震形变场。使用MAI方法时，首先将1幅影像重新分解为前视和后视2幅影像；然后对主、从影像的前视和后视影像分别进行影像配准、去平地相位和Goldstein滤波等一系列处理，得到其初始干涉图；最后对前视和后视干涉图进行差分处理及地理编码，得到MAI干涉图。
图 5(a)、5(d)为基于D-InSAR和MAI技术得到的ALOS-2降轨数据LOS向和方位向形变场，图 5(b)、5(c)为基于POT技术得到的Sentinel-1升降轨LOS向同震形变场。
|
图 5 LOS和方位向同震形变 Fig. 5 Coseismic deformation in LOS and azimuth directions |
基于生成的一维和二维SAR形变观测数据，采用SM-VCE方法获得此次地震的三维地表同震形变。相对于单点解算的传统加权最小二乘方法，基于窗口的SM-VCE方法可以获得更高精度的三维形变场[3-4]。由于SM方法会受到邻近点的约束，因此确定窗口大小为60像素×60像素，在该空间尺度下，POT观测中的空间高频噪声也可以被很好地抑制，从而使三维同震形变结果更加可靠。
由图 6可知，SM-VCE方法解算的三维同震形变场进一步证明，2022年门源地震发震断层以左旋走滑运动为主，地震破裂走向为NWW-SEE，最大垂直位移发生在断层附近。从整体上看，东西向(图 6(a))形变幅度比南北向(图 6(b))和垂直向(图 6(c))大得多，断层两侧水平位移在空间格局上总体呈对称分布，断层迹线北侧最大形变量接近1.26 m，断层南侧最大形变量约为1.58 m。南北向和垂直向位移场空间分布较为复杂，垂直向位移场在空间上分布极不均匀，最大沉降量仅为0.24 m左右，远小于东西向形变。
|
图 6 基于SM-VCE解算的三维同震形变场 Fig. 6 3D coseismic deformation field based on SM-VCE method |
由于InSAR形变场具有高空间分辨率，门源地震2组SAR数据形变场覆盖量大、相干性好且复杂，为了降低反演过程中数据冗余的影响，需要对其进行降采样处理。本文采用四叉树降采样方法[14]，采样后升轨和降轨数据分别保留344和484个数据点。
利用Sentinel-1卫星升降轨同震形变信息作为反演约束条件，使用大地测量贝叶斯反演软件包(GBIS)[5]，基于弹性半空间均匀滑动模型[15]，参考同震形变场和震源机制解给定的断层几何参数和相关运动学参数设置合理的搜索范围，其中初步设定2个断层的走向角分别为105°~115°和80°~90°，倾角均限制为70°~90°，并采用马尔科夫链蒙特卡洛搜索法进行搜索，当观测值拟合残差最小时，得到断层模型的最佳拟合解。其中，一组断层模型的最优走向角为109.23°、倾角为88.28°；另一组断层模型的最优走向角为86.6°、倾角为84.16°。这与三维形变场显示的断层迹线方向和野外地质调查结果[1, 10]基本一致。非线性反演的总地震矩约为9.6×1018 Nm，矩震级为MW6.62。图 7为InSAR观测模拟残差模型，由图可知，升降轨模拟模型与InSAR观测结果拟合度较好。
|
图 7 均匀式滑动模型拟合结果 Fig. 7 Fitting results of uniform slip model |
将上述非线性反演获取的震源几何参数和位置作为线性反演的先验信息，以三维地表形变为约束，求解运动参数的线性关系，可以更细致地反演发震断层内部几何形态。本文利用德国地学中心基于地壳分层介质模型的最速下降法(SDM)，通过最小二乘优化算法反演断层面滑动分布[6]。
将三维形变场数据的3个形变分量分别进行四叉树降采样处理，在东西、南北和垂直方向上分别保留1 015、430和739个数据点。反演中采用Crust1.0模型[16]，设置介质泊松比为0.25。本文综合分析多种数据，认为此次地震的发震主断层为冷龙岭断裂，次级断层为托莱山断裂，因此建立双断层模型进行反演。基于InSAR技术获取的形变场已经揭示地表破裂情况，将断层上边界设定为0，在此基础上将发震主断层长度(沿走向方向) 扩展至25 km，宽度(沿倾向方向)扩展至15 km，次级断层长度扩展至11 km，宽度扩展至15 km，然后细分为1 km×1 km的375和165个离散小块，根据输入的参数计算每个子单元上的滑动分布。实验中采用穷举法不断调整三维同震数据的权重，最终以东西向：南北向：垂直向为6 ∶2 ∶2的权重比进行计算。为提高反演精度，利用数据拟合度和模型粗糙度之间的折中曲线进行平滑约束，最终确定平滑因子a为0.08(图 8)。
|
图 8 最佳平滑因子曲线 Fig. 8 Optimal smoothing factor curve |
在InSAR三维同震形变场约束下，线性反演获取的2022年门源地震断层面最优滑动分布结果如图 9所示，图中箭头方向表示断层平面上划分的子断层的相对运动矢量方向，表明该断层呈现明显的左旋走滑特征。分布式滑动模型反演结果显示，震中位置为37.79°N、101.28°E。同震滑动主要集中在断层浅部区域，主断层深度范围为0~8 km，最大滑动量位于深度4 km左右，达到3.54 m，次级断层最大滑动量位于深度约3 km处，约1.59 m。反演获得的地震矩震级为MW6.67，对应地震矩为1.15×1019 Nm，与地震学结果较为一致。为检验反演结果的可信度，对三维形变场的观测数据和模拟数据进行对比分析，其拟合程度最高为98.71%。这表明线性反演的滑动分布结果可以准确反映发震断层的几何形态，也验证了三维形变场用于门源地震线性反演的合理性与可行性。东西、南北和垂直向同震形变场的拟合结果在形变量级和分布形态上与观测结果基本一致，残差分布主要集中在断层附近(图 10)，主要原因包括：1)断层近场区域形变梯度大且复杂，导致相位无法完全解缠，存在不连续现象；2)断层附近被积雪覆盖，存在部分失相干现象。
|
图 9 断层滑动分布结果 Fig. 9 Fault slip distribution results |
|
图 10 分布式滑动模型拟合结果 Fig. 10 Fitting results of the distributed slip model |
1) 2022年门源地震地表破裂沿NWW-SEE方向延伸，运动方式以左旋走滑为主。形变主要由南北向和东西向水平位移分量所贡献，最大地表水平形变为1.58 m。本次地震的主要滑动发生在冷龙岭断裂西段且延伸至托莱山断裂东段，破裂达到地表，与野外现场地质考察结果一致。
2) 2022年门源地震反演获取的地震矩为1.15×1019 Nm，相应矩震级为MW6.67，震中位置为37.79°N、101.28°E，主要滑动分布在地下0~8 km范围内，最大滑动量在深度4 km左右，达到3.54 m。门源盆地曾在1986年、2016年和2022年发生6级以上地震，导致冷龙岭断裂与托莱山断裂次生断裂繁多，地质构造复杂，在此地震学背景下，应关注该地区未来发生中强地震的可能。
| [1] |
李智敏, 盖海龙, 李鑫, 等. 2022年青海门源MS6.9地震发震构造和地表破裂初步调查[J]. 地质学报, 2022, 96(1): 330-335 (Li Zhimin, Gai Hailong, Li Xin, et al. Seismogenic Fault and Coseismic Surface Deformation of the Menyuan MS6.9 Earthquake in Qinghai, China[J]. Acta Geologica Sinica, 2022, 96(1): 330-335)
(  0) |
| [2] |
朱建军, 李志伟, 胡俊. InSAR变形监测方法与研究进展[J]. 测绘学报, 2017, 46(10): 1 717-1 733 (Zhu Jianjun, Li Zhiwei, Hu Jun. Research Progress and Methods of InSAR for Deformation Monitoring[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1 717-1 733)
( 0) |
| [3] |
刘计洪, 胡俊, 李志伟, 等. InSAR三维同震地表形变监测——窗口优化的SM-VCE算法[J]. 测绘学报, 2021, 50(9): 1 222-1 239 (Liu Jihong, Hu Jun, Li Zhiwei, et al. Estimation of 3D Coseismic Deformation with InSAR: An Improved SM-VCE Method by Window Optimization[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(9): 1 222-1 239)
( 0) |
| [4] |
Liu J H, Hu J, Li Z W, et al. Three-Dimensional Surface Displacements of the 8 January 2022 MW6.7 Menyuan Earthquake, China from Sentinel-1 and ALOS-2 SAR Observations[J]. Remote Sensing, 2022, 14(6)
( 0) |
| [5] |
Bagnardi M, Hooper A. Inversion of Surface Deformation Data for Rapid Estimates of Source Parameters and Uncertainties: A Bayesian Approach[J]. Geochemistry, Geophysics, Geosystems, 2018, 19(7): 2 194-2 211 DOI:10.1029/2018GC007585
( 0) |
| [6] |
Wang R J, Diao F Q, Hoechner A. SDM——A Geodetic Inversion Code Incorporating with Layered Crust Structure and Curved Fault Geometry[C]. EGU General Assembly, Vienna, 2013
( 0) |
| [7] |
Bechor N B D, Zebker H A. Measuring Two-Dimensional Movements Using a Single InSAR Pair[J]. Geophysical Research Letters, 2006, 33(16)
( 0) |
| [8] |
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 DOI:10.1029/1999GL900138
( 0) |
| [9] |
Guo P, Han Z J, An Y F, et al. Activity of the Lenglongling Fault System and Seismotectonics of the 2016 MS6.4 Menyuan Earthquake[J]. Science China: Earth Sciences, 2017, 60(5): 929-942 DOI:10.1007/s11430-016-9007-2
( 0) |
| [10] |
梁宽, 何仲太, 姜文亮, 等. 2022年1月8日青海门源MS6.9地震的同震地表破裂特征[J]. 地震地质, 2022, 44(1): 256-278 (Liang Kuan, He Zhongtai, Jiang Wenliang, et al. Surface Rupture Characteristics of the Menyuan MS6.9 Earthquake on January 8, 2022, Qinghai Province[J]. Seismology and Geology, 2022, 44(1): 256-278 DOI:10.3969/j.issn.0253-4967.2022.01.016)
( 0) |
| [11] |
Goldstein R M, Werner C L. Radar Interferogram Filtering for Geophysical Applications[J]. Geophysical Research Letters, 1998, 25(21): 4 035-4 038 DOI:10.1029/1998GL900033
( 0) |
| [12] |
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): 2 374-2 383 DOI:10.1109/TGRS.2006.873207
( 0) |
| [13] |
Xiao R Y, Yu C, Li Z H, et al. Statistical Assessment Metrics for InSAR Atmospheric Correction: Applications to Generic Atmospheric Correction Online Service for InSAR(GACOS) in Eastern China[J]. International Journal of Applied Earth Observation and Geoinformation, 2021, 96
( 0) |
| [14] |
Simons M, Fialko Y, Rivera L. Coseismic Deformation from the 1999 MW7.1 Hector Mine, California, Earthquake as Inferred from InSAR and GPS Observations[J]. Bulletin of the Seismological Society of America, 2002, 92(4): 1 390-1 402 DOI:10.1785/0120000933
( 0) |
| [15] |
Okada Y. Surface Deformation Due to Shear and Tensile Faults in a Half-Space[J]. Bulletin of the Seismological Society of America, 1985, 75(4): 1 135-1 154 DOI:10.1785/BSSA0750041135
( 0) |
| [16] |
L aske G, Masters G, Ma Z T, et al. Update on CRUST1.0——A 1-Degree Global Model of Earth's Crust[C]. EGU General Assembly, Vienna, 2013
( 0) |
2. Key Laboratory of Ecological Geology and Disaster Prevention, MNR, 126 Yanta Road, Xi'an 710054, China;
3. State Key Laboratory of Loess Science, 126 Yanta Road, Xi'an 710054, China