关于荣县地区三维同震形变场及荣县地震的反演分析作为近年来发展迅速的空间大地测量技术,合成孔径雷达干涉测量(InSAR)具有形变监测精度高、覆盖面积广、重访周期短且全天候全天时数据获取等突出优势,在地震周期形变监测研究中发挥着越来越重要的作用目前,国内外利用InSAR技术成功获取小型地震(震级小于MW5.0)同震形变的研究案例较少本文以2019-02-24、02-25发生在四川省荣县境内的地震为研究对象,其参考震中为29.47°N、104.48°E,震级均在5级以下荣县为我国页岩气开采示范区,其经济、社会意义重大,但该区地质结构不稳定,较易发生地震,且目前该区域InSAR观测研究仍为空白因此,本文利用覆盖荣县地区的1对升降轨数据和近5个月的SAR影像集,采用DInSAR、时间序列InSAR技术及均匀滑动模型,反演获取地震的三维同震形变、震源参数及近5个月的形变累计信息,为了解地震的发生机制及该区域的形变过程等提供重要参考1、研究区与实验数据1.1研究区概况本文研究区以荣县地震震中为中心,覆盖面积约22km×22km,该区域地处四川省东南部,覆盖自贡市荣县东部、内江市威远县西部和自贡市贡井区小部分辖区。
研究区属于长宁-荣县-威远国家级页岩气示范区,目前已有多处页岩气开采实验平台投入生产自2018年以来,该区域发生2.5级以上地震20余次1.2实验数据为研究地震对荣县地区地表造成的形变影响,选择欧洲空间局Sentinel-1AC波段SAR影像作为数据源,其极化模式为VV极化,距离向和方位向分辨率为5m×20m,选用的Sentinel-1A数据成像模式为TOPS模式,其中降轨影像11景,升轨影像2景,具体获取时间见表1原始影像幅宽约为250km,为提高数据处理效率,截取长宽均约为22km的区域作为研究区此外,为提高轨道精度,采用欧空局发布的精确轨道数据(PDGS),且在数据处理过程中引入外部数字高程模型(DEM)数据TanDEM-XDEM时间序列处理使用11景降轨影像,同震三维形变场解算及震源参数反演分别使用2景升轨、降轨影像,组成的干涉对为2019-02-14和2019-02-26(升轨)、2019-02-16和2019-02-28(降轨)时间序列研究选取2018-11-12的影像作为公共影像,其他影像为副影像分析11景SLC影像的时间基线距和垂直基线距,根据小基线对组合原则——较小的时间基线距阈值和垂直基线距阈值形成干涉图组合,将时间基线阀值设为90d,空间基线阀值设为250m。
表1Sentine-1A影像数据2、数据处理与方法2.1同震三维形变解算采用二轨DInSAR方法进行数据处理该方法的基本思路为:利用SAR影像对生成干涉图,引入外部DEM并将其转换为主影像雷达坐标系下的模拟SAR干涉纹图,即可将地形相位从干涉图中去除,获得所需形变相位二轨DInSAR法的基本原理可表示为[1]:φdef=φwrp−φflt−φtop−φatm−φnoise−2kπ (1)式中,各项自左向右依次为雷达视线(LOS)向地表形变相位、未解缠的差分干涉相位、平地相位、地形相位、大气相位和噪声相位,k为整周模糊度获得地表形变相位后,通过几何关系可以获得地表在各方向的形变分量单轨差分干涉测量难以获取三维地表(EW向、NS向和UD向)形变,即DInSAR技术存在LOS向模糊问题[2]首先采用Gamma软件[3]对研究区地震前后的2组数据(升轨和降轨)进行处理,即分别选取地震前、地震后的SAR影像作为主副影像,对SLC影像对进行配准、重采样等处理,获得升降轨干涉图,并对干涉图进行滤波处理;再利用30m分辨率的TanDEM-XDEM数据去除地形相位,并采用最小费用流解缠(MCF)算法对2组差分干涉图进行解缠处理;最后通过形变分量解算获得2个飞行方向的LOS向形变。
由升降轨数据分别解算的LOS向形变、三维地表形变、升降轨雷达传感器入射角和升降轨轨道方位角之间存在2个恒等式[4,5],但2个方程无法解算3个未知数考虑到本次地震断层走滑方向的形变量接近于0,且断层走向方位角γ可通过反演获取[6],因此加入约束方程:dNcosγ−dEsinγ=0 (2)将3个等式联立解算方程组,即可获得同震三维形变场2.2SBAS-InSAR方法处理时间序列分析为一种能有效克服DInSAR技术不足的InSAR策略,目前主要有永久散射体雷达干涉测量(PS-InSAR)和小基线子集雷达干涉测量(SBAS-InSAR)两种方法,综合考虑2种方法的特点及研究区的环境特性,本文选择SBAS-InSAR方法SBAS时间序列分析方法的基本原理为[7]:假设覆盖研究区域的SAR影像有N+1景,按时间先后顺序排列首先选取1景影像作为公共主影像,将其余N景影像配准到该主影像;随后依据时空基线组合原则确定干涉对组合,假设为L个;然后分别对每个组合进行差分干涉处理,得到M幅干涉纹图,M需满足关系式(N+1)/2≤M≤N(N+1)/2将第1景影像的获取时间作为初始时刻,将剩余N个时刻相对于初始时刻的差分相位作为未知参数,将数据处理得到的差分干涉相位作为观测量。
假设所得差分干涉图都已经过正确解缠,且差分干涉相位已经被校正至某个稳定的或形变已知的高相干点像元(参考像元)上,则可组成2个时间序列;对于第k幅差分干涉图(某2个时刻)中的任意像元(x,r),在不考虑大气延迟相位、地形残余相位和噪声的情况下,其相对于初始时刻在雷达LOS向的累计形变量与波长存在系数已知的等式关系将每对雷达干涉影像对经过配准、差分干涉处理后,其条纹相位值处于-π～π之间,采用最小费用流算法对其进行解缠处理然而解缠后的相位中依然含有干扰相位——大气延迟误差、轨道误差和DEM误差等,只有将其从解缠相位中去除才能得到由LOS向形变引起的相位部分[7]荣县、威远县降雨充沛、气候湿润,而大气中水汽含量较高会对时序InSAR分析产生较大影响,因此需对大气效应进行估计并加以去除本文采用的策略为:首先借助大气延迟相位变化在空间维低频和时间维高频的特性,通过时空维的带通滤波来消除部分大气延迟的影响;其次估算与地形相关的大气延迟相位,在保持地表非线性形变信号的情况下,最大限度地消除大气效应对最终沉降监测结果的影响轨道误差和DEM误差则通过时间和空间维滤波并进行最小二乘解算,即可得到对应相位值。
最后采用SVD反演算法获取荣县地区的地表形变参数——形变时间序列[7,8]3、结果与讨论3.1三维同震形变图1、2为利用二轨DInSAR技术获取的升降轨荣县地震LOS向同震形变场从图中可以看出,升降轨LOS向形变空间分布较为一致,表现为升降轨同震形变场在相同区域均有明显的抬升信号,其形状似椭圆形(长轴大体沿EW向);但在该区域两侧,升降轨结果表现有所不同,分析认为可能为雷达飞行方向不同造成视角不同所引起图中绿色圆圈标记为中国地震台网发布的3次地震的参考震中位置,可以发现,该位置大体沿椭圆形抬升区域分布,周围出现严重的地表形变现象,最大LOS向形变约17mm,最大沉降约12～15mm,表明本次荣县地震参考震中与InSAR技术监测到的形变中心位置基本吻合图1升轨雷达视线向形变图2降轨雷达视线向形变降轨同震形变场(图2)显示,形变影响范围可达10km×10km,被2条NS向分界线分割为左侧、中间和右侧3个部分,其中中间区域呈现抬升信号,存在明显的抬升中心,LOS向形变为12mm;而左、右两侧均为沉降信号,且右侧下降信号更为明显,有2个沉降中心,相比之下,左侧沉降趋势较缓和升轨同震形变场(图1)显示的抬升信号区域与降轨同震形变场相似,其LOS向形变为17mm;左侧沉降与降轨同震形变场相比,下沉现象更为明显,且存在1个沉降中心,LOS向形变为-12mm;升轨同震形变在右侧对应区域未获取到明显的形变信息。
为更加直观地了解荣县地震的震源机制及对该地区造成的影响,将前期处理得到的升轨、降轨DInSAR处理结果进行联合,使用§2.1中解算策略,获取该地区三维地表形变(图3、4)图3EW向同震形变图4UD向同震形变荣县地震三维同震形变场解算结果显示,椭圆形抬升区域及其右侧大部分沉降区域向E移动,最大偏移量为8mm;左侧沉降区域则向W运动,最大偏移量为9mm,该运动特征符合逆断层特点地震烈度是指地震时某一地区的地面和各类建筑物遭受到地震影响的强弱程度,烈度图可以反映地震区域内建筑物的破坏程度和地表变化状况[9]将四川省地震局发布的荣县地震烈度图[10]与InSAR监测结果进行对比可以看出,烈度覆盖范围与InSAR技术获取的地面形变区域空间分布一致,其中成佳场断裂总体沿NS向延伸三维同震形变中UD向分量结果表明,该区域存在左、右2条长度约8km的断裂,呈现中间抬升、两侧下沉的现象,走向也为NS向图5、6为跨地震断层沿AA′剖面的EW向和UD向位移,其中EW向形变量约为-6～6mm,UD向形变量为-10～15mm,因此地震引起的形变以UD向为主综合分析升降轨LOS向形变和三维同震形变可知,接近地震震中参考位置的椭圆形区域为受荣县地震影响最大的区域,结合该地区行政区划图可知,该区域分布在荣县高山镇,大致沿三新店-窑嘴上-营盘村一带分布,其UD向形变量最大约为18mm。
椭圆形抬升区域两侧为地震造成的下沉区域,主要在荣县正安乡、甘家湾和半边店村附近查阅相关信息可知,荣县地震对高山镇房屋、道路及水库等人工建筑造成不同程度的损坏,这与InSAR监测结果相吻合图5沿图3剖面AA′的EW向形变图6沿图4剖面AA′的UD向形变3.2震源参数反演由于InSAR数据量较大,为提升实验效率,需要对形变结果进行降采样处理,在保留形变信息的同时减少参与运算的数据点数本文借助GBIS软件进行反演工作,该软件可利用InSAR数据快速反演震源信息[11]结合前期升降轨形变结果,将迭代次数设置为100万次,由Okada弹性半空间模型求取模型最佳拟合解,并根据模型参数的不确定性对模型参数进行约束反演,相关参数见表2表2荣县地震反演结果由断层参数的后验概率密度函数分别获得2.5%和97.5%的最大后验概率解,结果表明,同震地表形变由断层引起,该断层长约9.5km,宽约2.1km,深度约1.8km,倾角-89.7°～-86.62°,走向351.12°～359.4°荣县地震的滑动方向与具有右旋分量的逆冲断层的滑动方向一致,沿走向方向为-0.2～-0.1m,沿倾向方向为0.05～0.08m。
断层深度远大于宽度表明断层未破裂至地表,走滑为负分量表明上盘运动方向与走滑方向相反,倾向滑动为正分量表明上盘相对隆起,这些现象表明荣县地震的发震断层符合逆冲断层的特性对比荣县地震烈度图可以发现,反演的断层走向与成佳场断裂走向一致3.3时间序列形变图7为由SBAS时间序列分析法解算的实验结果将2018-11-12视为未发生形变的基准,以各时间节点的累计形变量来展示此次荣县地震震前和震后的地表垂直动态变化分析结果表明,2018-11-12～2019-02-16时间段内,该区域基本处于微小形变状态,累计形变较小,此时震中附近最大累计量约为±5mm2019-02-16～2019-02-28跨越发震时刻(2019-02-24～25),由图7(i)可知,该时间段内震中周围累计形变量突增,最大值超过±20mm2019-02-28以后,该区域累计形变无较大变化,基本可以判定震后地表形变趋于稳定为更加直观地了解形变历史,选取图7中具有代表性的A、B两点分析其形变时间序列,A点位于抬升区域,B点位于沉降区域图8中红线为此次地震的发震时刻(2019-02-24～25),从图中可以看出,红线所处线段为形变速率最大的一段。