运用InSAR技术解算四川长宁MS6.0地震三维形变场及时间序列形变分析

于书媛, 方良好, 宴金旭, 倪红玉, 陈 靓, 丁 娟

(1. 安徽省地震局, 安徽 合肥 230031; 2. 安徽蒙城地球物理国家野外科学观测研究站, 安徽 蒙城 233527;3. 四川省地震局, 四川 成都 610041)

0 引言

北京时间2019年6月17日22时55分43秒,四川长宁发生6.0级地震。参考震中位于104.90°E,28.34°N,震源深度16 km。重庆、四川、贵州、云南等多省对此次地震有感。截至2019年6月26日8时,共记录到大于2.0级余震182次,其中大于MS4.0余震7次,最大余震震级5.4级。此次地震最大烈度达到Ⅷ度,截至2019年6月21日,地震共造成13人死亡,226人受伤。自长宁MS6.0地震发生后,相关地震行业专家从波形资料、地壳形变、重力场、地电场、地磁场、地质构造等方面对地震构造应力场特征、孕震机理、构造特征进行深入探索[1-3]。我国地震多发生在高原山区,震中附近地质构造复杂,植被覆盖,人不能至,传统手段无法获取连续、宏观的形变场。合成孔径雷达干涉测量技术(Interferometric Synthetic Aperture Radar,InSAR)自1993年首次应用于地震以来,已在全球地震研究中得到快速发展和成熟应用。本文利用覆盖长宁县地区的哨兵1A升降轨SAR数据(http://www.asf.allaska.edu/sar-data/palsar/),采用InSAR技术获取地震的三维同震形变场,并进一步研究震源位置,反演断层几何参数和近6个月的形变累积场,讨论其发震原因与双河背斜构造关系,弥补该地区活动构造研究的不足,为判定区域地震形势提供参考。

1 研究区及InSAR数据

1.1 研究区概况

本文研究区以长宁县地震震中为中心,覆盖面积约为30 km×30 km,长宁地震震中位于青藏高原东缘、云贵高原北缘地壳运动较为活动的区域,地貌位置表明西南侧是抬升的高原,北侧是相对沉陷的四川盆地。此次地震发生于川西南地区,其断裂活动程度及地震强度虽不能和横断山东部及其川滇地震带相比,但仍然是中强震活动较为频繁的区域,研究区构造情况如图1所示。

1.2 研究数据

欧洲航天局(ESA)2014年4月3日发射的Sentinel-1A卫星采用12天重访周期进行全球覆盖,其获取了长宁地震震前震后多景SAR影像。具体获取参数见表1。采用的外部参考DEM数据为NASA SRTM DEM,数据精度为30 m。

表1 升降轨差分干涉影像参数

时间序列研究选取降轨2019-07-26作为超级主影像,其他影像作为副影像,具体获取参数见表2。分析12景SLC影像的时间基线距和垂直基线距,在已有SAR影像数据集中形成基于不同主影像的时间序列干涉图子集,利用短基线合成孔径雷达干涉测量(Small Baseline Subset-InSAR,SBAS-InSAR)求解近6个月的震间形变累积、削弱信息。

表2 时间序列Sentine-1A影像参数

2 InSAR数据处理及三维同震形变场

常规单一的一对数据进行差分干涉测量(Differential InSAR,D-InSAR)处理,获得的是雷达视线向的形变,可能存在获取的差分干涉结果形变量差异较大的情况,因此单一方向的形变结果不能准确反映研究区的真实三维形变情况。本文结合升、降轨影像数据,采用D-InSAR方法获取沿LOS向形变信息并进行三维同震形变场联合解算,获取同震东西、南北、垂直三个方向上的形变量。

2.1 D-InSAR获取LOS向同震形变

本文运用震前震后各2景升轨、降轨Sentinel-1A幅宽约为250 km数据通过截取研究区经差分干涉处理获取震区三维同震形变场,其中升轨组成的干涉对为2019-06-09和2019-06-21、2019-06-14和2019-06-26,降轨组成的干涉对为2019-06-16和2019-06-28,时间基线均为12 d。利用商业软件ENVI的SARscape模块对Sentinel-1A数据进行SAR差分干涉技术(D-InSAR)处理[4],由于选用的三对干涉对具有较高的相干性,获取了长宁地震的形变场。从图2中可看出,升轨的雷达视线向形变场整体以沉降为主,降轨的视线向形变场整体以抬升为主,且升降轨卫星数据同区域形变场观测值符号相反,其形变场范围大致为18 km×15 km,形状似椭圆形(长轴大体呈北西向),区域地质构造与长宁背斜构造附近规模不大的次级断层对应。卫星LOS向最大沉降量(远离卫星传感器)分别是7.3 cm(升轨)和6.26 cm(降轨),最大抬升量(靠近卫星传感器)分别是5.5 cm(升轨)和8.12 cm(降轨)。综合本次地震的发震构造和升降轨LOS向沉降和抬升之间的关系,判定本次形变场是由断层的走滑和逆冲综合运动形成。根据形变场抬升量、沉降量分布及区域地质构造背景,确定一条北西—南东走向的断层迹线(图2)。结合发震地区1∶5万地质图和地质构造背景分析,长宁地震发生在四川盆地南缘盆山转换带,发震区内主要地质构造为褶皱及伴生断层,出露较少,走向和规模多变。本文形变场的整体呈NW-SE走向,其分布特征与长宁—双河复式大背斜主体构造一致,认为本次形变场推断断层水平位置与双河乡处存在的NW向伴生断层存在重合。图3为跨断层沿AA′剖面的升轨和降轨位移剖面图,其中升轨LOS向形变量约为-4～1 cm,降轨LOS向形变量约为-1～4 cm。

图2 D-InSAR雷达视线向形变Fig.2 D-InSAR LOS deformation

图3 沿图2剖面AA′的升降轨形变Fig.3 Ascending and descending deformation along the profile AA′ in figure 2

2.2 三维同震形变场解算

采用常规D-InSAR方法对升、降轨SAR数据对提取的雷达视线向形变量为一维方向的形变量,由于一维的形变量是由东西、南北、垂直三个维度方面投影得到的,因此通过D-InSAR技术处理得到的形变信息存在视线向模糊问题[5-8]。本文为解算上述3个方向的形变分量,设定dE、dN、dU为同震形变在东、北及垂直方向上的分量,排除卫星轨道误差、测量误差等因素,通过地表的空间几何与卫星视线和轨道的关系得到函数式,如式(1)。

dU·cosθ-dN·cosα·sinθ-dE·sinαsinθ=dLOS

(1)

式中:θ为雷达视线方向的入射角;α为方位角;dE、dN、dU分别为分解后同震形变在东南、西北、垂直方向上的分量;dLOS差分干涉测量的同震形变量。同时,为了解算同震形变三维方向的形变量,原则上需要3个视线方向的形变结果。因此,本文为解算更可靠的3D分解同震形变结果,结合2个升轨LOS向同震形变场和1个降轨LOS向同震形变场,采用加权最小二乘法反演(权重来自于LOS位移精度)。其中,LOS位移精度由相干性和波长导出,提供测量精度的估计值(即标准偏差值)。这个估计的标准偏差值越高,测量值精确度越低,精度公式表达为:

(2)

式中:γ表示干涉相干性;λ表示波长。

本文计算出的三维形变场EW向形变量(向E为正)、SN向形变量(向N为正)、垂直向形变量(向上为正)。长宁地震三维同震形变场解算结果显示,如图4(c)所示,地震形变以向北运动为主,向北最大偏移量达到17.8 cm。结合图4(a)、(b)所示,椭圆抬升区域及其西南侧大部分沉降区域向E移动,最大偏移量为8 cm;北东侧沉降区则向W运动,最大偏移量为7 cm,该同震形变场符合逆断层特征。由于地震烈度图可反映地震区域内建筑物的破坏范围、破坏程度以及地表变化情况[9]。本文将四川地震局发布的长宁6.0级地震烈度图与InSAR监测结果进行对比可以看出,烈度范围与InSAR技术获取的地表形变区域空间分布较为一致,发震断层总体沿NW向延伸。三维同震形变中UD向分量结果表明,该区域存在NW向共轭的2条长度约15 km的断裂,呈现中心抬升、两侧下沉的特点,形变场特征走向也为NW向。

图4 三维同震形变Fig.4 3D coseismic deformation

综合分析升降轨LOS向形变和三维同震形变特征可知,震中参考位置位于椭圆形变场的西南区域,结合研究区行政区划图和地震烈度图可知,椭圆形形变场区域为受到长宁地震影响最大的区域,该区域分布在长宁县双河镇、富兴乡,其UD向形变量最大约为8 mm。椭圆形抬升区域两侧为地震造成的下沉区域,主要在长宁双河镇附近,结合相关应急救援资料可知,长宁地震造成的双河镇伤亡人员数量、人工建(构)筑物的损坏,这与InSAR监测结果吻合。

3 同震滑动分布反演

本文运用两轨升降轨同震形变场等权为反演约束条件,基于Okada弹性半空间位错模型[10],反演地震断层滑动分布。同震滑动震源参数反演步骤包括:

(1)对升轨2019-06-09和2019-06-21、降轨2019-06-16和2019-06-28得到的D-InSAR形变场通过绘制矢量多边形图层,划定形变场和反演模型范围,形变场采样间隔是1 km,反演模型范围采样间隔为2 km,升降轨形变场分别采样得到1 320、976个数据点,两轨数据等权重值作为反演断层滑动的约束条件,震中形变场的降采样结果、模拟形变场和残差分布如图5所示。图5(a)、图5(b)、图5(c)和图5(d)、图5(e)、图5(f)分别为升轨和降轨的观测值、模拟值和残差,升降轨误差集中在-1～1 cm区间。

(2)将得到的升、降轨InSAR视线向形变场结果,采用非线性方法反演均匀滑动断层几何参数。为了准确进行反演,加入6个轨道参数进行线性估计轨道误差,Okada模型参数设置采用GCMT的震源机制解进行约束。其中,采用Levemberg-Marquardt最小二乘优化算法迭代,进行9个几何参数和6个轨道参数求解。

(3)最后基于上一步反演的断层参数采用线性反演方法并结合非负最小二乘算法求解活动断层精细滑动分布。断层参数反演结果如表3所列,发震断层二三维模型如图6所示。图6(a)、图6(b)分别为坐标系下二维同震滑动分布、二维断层滑动分布。通过断层几何参数反演结果表明本次地震的宏观震中为104.87°E,28.39°N,矩震级为MW5.8,破裂尺度模型断层长度20 km,宽度5 km,破裂深度主要集中在2～13 km范围内,震源深度10.2 km;断层倾角为56.42°,平均滑动角为44.37°,断层走向319.3°。此次地震时兼具逆冲和走滑的破裂模式,对比长宁地震烈度图可看出反演出断层的走向与烈度圈的方向一致。此外,如图6所示,基于DInSAR形变结果反演的地震形变宏观震中与多家研究结构基于地震波确定的震中位置具有比较大的差异,分析原因可能与台站的分布情况、观测手段及精度对研究机构的速报结果精度均有影响。同震形变场得到的宏观震中更能真实地反映中小震源发震位置,为地震监测预报中心提供更加准确的断层位置和几何信息。

图5 形变场采样及反演残差Fig.5 Sampling results and inversion residuals of deformation field

图6 长宁地震同震破裂模型的滑动分布Fig.6 Co-seismic slip distribution of rupture model in Changning earthquake

图7 震中位置Fig.7 Epicenter location

4 SBAS-InSAR震间形变分析

4.1 时间序列形变

时间序列InSAR观测结果可从形成场时空分布特征和断层活动状态进行活动断裂带地壳运动的精细结构研究,并可从断裂带及近场区域的宏观面状动态形变图像研究断层震前-震间-震后累积形变量、活动速率等[10-11]。本文结合研究区域的相干性和使用的数据量,选择SBAS-InSAR方法[12],计算平均形变结果。

4.2 连续空间形变场模拟

克里金插值是基于包含自相关(测量点之间的统计关系)的统计模型[13]。本文采用克里金插值方法对离散的形变场矢量点进行插值,图8为由SBAS时间序列方法解算和克里金插值方法得到的形变场结果。将2019-03-12作为未发生形变的时间基准,以11个时间节点的累积形变量来展示此次长宁地震震前、震后的地表LOS向形变动态演化过程。分析结果表明,2019-03-24—2019-06-16时间段内,该区域范围基本处于较为微小的形变量,形变场范围累积形变量较小,震中近场附近最大累积形变量约为±8 mm。2019-06-16—2019-06-28时间为发震前后时刻(2019-06-17)。由图8(h)可知,该时间震中附近累积形变量突增,最大值超过±20 mm。2019-06-28—2019-07-22期间,形变场累积形变量略有衰减。2019-08-03—2019-08-15期间,震中附近累积形变场又出现增长,最大抬升值超过30 mm,最小沉降值超过20 mm,分析可能由于本次地震的余震较为丰富是的震后地表形变处于一定时间的不稳定状态。

图8 研究区形变累积时间序列Fig.8 Accumulative time series of deformation in the study area

为较为直观地分析震间形变过程,选取图8(i)中具有典型形变特征的A、B两点分析形变时间序列,A点位于抬升区域,B点位于沉降区域。图8中红线为长宁地震的发震时刻(2019-06-17),从图9可看出,红线所处的时间段为形变速率最大的一段。地震发生前A点的形变波动表现为缓慢抬升,沉降量略大于抬升量,而B点沉降量始终大于抬升量,总体呈现沉降趋势,震中近场其他点呈现类似的形变趋势。

图9 图8中A点和B点形变时间序列Fig.9 Deformation time series at points A and B in figure 8

5 讨论

本次地震为典型的构造地震活动,区内仅有少量断层出露地表,地震发生于四川盆地边缘的长宁背斜构造规模不大的次级断层上,在已有的活动断层数据库中并无标识,给发震构造判识带来很大困难[14],是一次典型的盲断层破裂事件。发震区域位于断裂活动性较弱、分布式构造变形的川西南区域,西侧为长800 km的华蓥山褶皱带,形成北西—南东向的现今应力场[15],光学影像上没有显著的断层构造特征。DInSAR同震形变场表明,本次地震产生了18 m×15 km的形变场范围,卫星视线向最大沉降量是7.3 cm(升轨)、最大抬升量是8.12 cm(降轨),断层形变场走向为NW向,最大偏移量达到17.8 cm(向北)。反演结果表明,滑动量在断层深度2～6 km处达到最大,以逆断走滑为主,最大滑动量约0.06 m,位于沿断层面倾斜陷落深度3 km处,震级为MW5.8。

本研究得到的发震断层为逆断层,根据形变上升的主动盘位于西侧山脊区域,而被动盘位于东侧海拔较低的低洼区域,同震形变与地貌特征耦合。根据郭志[16]的研究显示长宁地震主震级几次震级较大余震的震源机制结果表明,发震断层以逆冲为主,结果与长宁—双河背斜的地质构造及应力背景较为一致,断层倾角为59°,略大于本研究中得到的发震断层的倾角(45°)。

截至2019年7月28日,长宁地震主震后发生余震6 888次,其中包含6次M>4.5的余震,相对于主震呈北西向条带分布,总长约20 km,正对应构造应力较为集中的背斜轴部区域,存在众多背斜伴生断层[2,17]。在此次地震震中南侧约15 km处曾发生2018年12月16日四川兴文M5.7和2019年1月3日珙县M5.3地震。易桂喜等[2]通过长宁6.0级地震序列16次MS≥3.6地震序列震源机制解统计认为发震断层高倾角,断层走向以NW-SE和NS走向略占优势,发震构造复杂,构造应力环境有利于断裂发生逆冲运动。尹欣欣等[17]利用四川长宁6.0级地震的主余震资料,利用CAP震源机制解方法对14次3级以上的地震序列做了震源机制解计算,认为应力场反演结果显示主震区域网格(28.3°N,104.8°E)断层为逆冲类型,破裂长度大致为35 km。孙凯等[18]采用SDM反演InSAR数据源,迭代计算发震断层面最优滑动分布,最大滑动量0.28 m,平均滑动角50°,与本文反演结果基本一致。总体来说,认为本次地震发生在上地壳浅部,发震断层滑动面上缘接近地表,推测由主震引起的次级断层活动触发短期内强余震频发,认为该区域未来需要进一步加强人口密集地区的地震风险防范工作。

6 结论

本文基于升降轨Sentinel-1A数据,采用DInSAR技术获取2019年6月17日四川长宁MS6.0地震升降轨同震形变场和加权最小二乘法反演三维同震形变场。分析结果表明,长宁县地震形变升轨LOS向形变场以沉降为主(最大沉降量7.3 cm),降轨LOS向形变场以抬升为主(最大抬升量8.12 cm),且升降轨卫星数据同区域形变场观测值符号相反;三维同震形变场以NS向为主,EW向形变量约为-7～8 cm,UD向形变量约为-4～7 cm;形变场空间分布特征与断裂走向为NW向,且与地震烈度图的空间分布特征较为一致。利用升降轨形变场采用均匀滑动模型反演结果表明断层几何参数反演结果表明本次地震的破裂尺度模型断层长度20 km,宽度5 km,破裂深度主要集中在2～13 km范围内,震源深度10.2 km;断层倾角为56.42°,平均滑动角为44.37°,断层走向319.3°,与CMT、USGS、中国地震局地质研究所、孙凯等结果基本一致[18]。SBAS-InSAR震间形变结果表明,该区域近6个月内区域累积形变场变化特征主要由长宁地震引起,形变场动态演化特征明显,震前时间段该区域范围基本处于较为微小的形变量,形变场范围累积形变量较小,临震时间段震中附近累积形变量突增,震后短期一段时间形变场累积形变量略有衰减,但最后一段时间震中附近累积形变场又出现增长,最大抬升值超过30 mm,最小沉降值超过20 mm,分析可能由于本次地震的余震较为丰富是的震后地表形变处于一定时间的不稳定状态。