巴姆地震三维形变场获取

张庆云, 张景发, 李永生, 李兵权, 罗三明, 马庆尊

(1.中国地震局第一监测中心, 天津 300180； 2.应急管理部国家自然灾害防治研究院, 北京 100085)

随着合成孔径雷达干涉测量技术(interferometric synthetic aperture radar, InSAR)的不断发展,在地震发生后,基于合成孔径雷达(synthetic aperture radar,SAR)数据快速获取地震同震形变场,可以对地震造成的破坏程度和范围进行直观描述。虽然使用差分干涉测量(differential interferometric synthetic aperture radar,D-InSAR)进行同震形变场获取的技术已经相对成熟,但对于地形复杂地区,仍然会存在地形误差影响；而且传统D-InSAR技术获取的同震形变场是卫星视线向(line of sight,LOS)结果,无法反映同震形变的全貌,因此三维形变场的解算对深入了解发震构造以及地震造成的形变具有重要意义。

目前结合InSAR技术获取三维形变场的方法主要有:①不同传感器升降轨InSAR数据分析,该方法的缺点是对数据要求较高,而且由于SAR独特成像原理,导致其对南北向不敏感,所以结合升降轨的三维解算在南北向存在较大误差,同时目前不同传感器的入射角和方位角相差不大,也会导致解算存在误差[1-2]；②InSAR技术联合外部数据全球定位系统(global positioning system,GPS)数据、或者根据模型模拟得到的三维结果进行三维解算,该方法缺点是有些地区没有GPS覆盖,模型模拟会引入一定程度的模型误差[3-6]；③不同InSAR技术联合进行三维解算,D-InSAR技术获取LOS向形变,多孔径干涉雷达测量技术(multiple aperture interferometry,MAI)可以获取方位向形变场。偏移跟踪技术(offset-tracking)可以获取地震的方位向和距离向形变结果,综合使用升降轨数据的LOS向形变以及方位向或者距离向形变可以获取三维形变结果[7-10]。

经过中外学者多年努力,目前InSAR三维形变场获取技术已日益成熟趋于完善,在同震、地面沉降和滑坡等各类地质灾害形变监测中得到广泛应用。He等[11]、张庆云等[12-13]利用InSAR技术获取了2017年伊朗地震的三维形变场；Chen等[14]、姚佳明等[15]和董龙凯[16]利用InSAR技术进行了矿区三维形变场解算研究；Shi等[17]和王晨兴等[18]基于InSAR技术研究了滑坡三维形变特征,并取得了可靠的分析结果；杨丽叶等[19]利用InSAR技术分析了错郎玛冰川的三维形变特征；Pepe[20]和马艳鸽[21]分析了InSAR技术获取三维形变场的方法原理。虽然已有的研究成果已经相对成熟,但是仍存在不足之处,如MAI技术和D-InSAR技术对SAR影像相干性要求相对较高,导致地震、矿区等形变量较大区域,存在严重失相干；以及在地形复杂区域,地形误差如何准确去除等。失相干以及各类误差的去除对三维形变场解算结果的准确性都有重要影响,这些都是当前仍需要解决的主要问题。

针对InSAR技术获取的形变场中存在失相干及地形等相关误差去除问题,以常用的第3种三维形变场解算方法为基础,以2003年12月26日巴姆地震为例,研究地形误差去除方法、完整形变场获取技术,基于D-InSAR和MAI技术进行三维形变场解算,分析地震三维形变场结果。

1 同震形变场获取

2003年12月26日,伊朗东南部巴姆(Bam)古城发生Mw6.6地震,震中位置在29.01°N,58.26°E,地震造成四万多人死亡,五万多人受伤,成千上万人无家可归,巴姆古城内几乎所有的建筑物都被损毁,地震造成了严重的人员伤亡和财产损失。

地震发生后,Zare等[22]通过分析伊朗国家强震台网获取的强震记录,认为Bam断层长约10 km,且经过巴姆市；Tatar等[23]分析了巴姆地震余震分布特征；Wang等[24]使用先进合成孔径雷达(advanced synthetic aperture radar,ASAR)数据分析了巴姆地震发震特征；季灵运等[25]、王洪友[26]使用ASAR数据获取了巴姆地震三维形变场,但是上述研究存在一定的局限性,形变场中存在一定程度的失相干现象和地形误差。在已有研究基础上,对地形误差进行进一步去除,同时对失相干区域进行恢复,然后基于多种技术联合的方法进行三维形变场解算。图1为巴姆地震区域构造图。

红色线表示研究区范围内断层分布；红色五角星表示主震位置；绿色圆圈表示余震分布；蓝色虚线框为获取的ASAR数据分布情况；Gowk和Mayband为断层名称图1 2003年12月26日巴姆Mw6.6地震区域构造Fig.1 Regional structure of Bam Mw6.6 earthquake on December 26, 2003

1.1 LOS向形变场获取

地震发生后,欧洲航天局ENVISAT卫星获取了地震震前和震后SAR数据,由于该区域植被稀疏且以平原为主,因此适合用SAR数据分析其形变场。选择合适的SAR数据,获取地震的同震形变场,表1为主要数据参数。

使用两轨差分方法,以GAMMA软件为InSAR处理平台,使用航天飞机雷达地形测绘使命(shuttle radar topography mission,SRTM)获取的分辨率为90 m的数字地形高程模型(digital elevation model,DEM)对干涉图进行去地形处理,再通过滤波、相位解缠以及轨道误差校正,最后得到地理编码后的同震形变结果(图2)。从图2(a)可以看出,降轨获取的形变场结果左下角存在一定的地形影响，而且由于地震震级较大,极震区破坏严重,导致升轨和降轨获取的同震形变场中极震区都存在一定的失相干现象,失相干区域的存在会影响后续三维形变场解算的准确度。

图2 基于ASAR数据获取的巴姆地震同震原始形变场Fig.2 Coseismic original deformation field of Bam earthquake based on ASAR data

为保证获取可靠的形变场,针对降轨中存在地形影响的区域进行误差去除。首先勾选出存在地形影响的区域,然后对存在地形误差区域的参考DEM与形变场经纬度之间建立函数模型,进行地形误差去除。最终使用的模型为

Z=[Ddisp-0.021 237-0.002 076LlatDdem-(-0.002 617 8)LlonDdem-(-9.989×10-6)Llat2Ddem-(-1.608 2×10-5)Llon2Ddem-(-2.580 5×10-5)Llon2Ddem-(-0.106 15)Ddem]

(1)

式(1)中:Z为去除地形误差之后点的形变量；Ddisp为未去除地形误差时点的形变量；Llat为点的纬度；Llat为点的经度；Ddem为点的高程。

针对极震区失相干现象,采用基于非线性支持向量机(support vector regression,SVR)方法进行失相干恢复,首先对SAR数据获取的形变场进行分类,对形变场中非0区域进行数据预处理,建立训练集和测试集；对形变场中0区域进行一定的数据预处理；然后结合训练集、测试集以及0结果进行核函数的选择,进行SVR模型训练和模型预测,最终获取完整的形变场。

图3、图4为经过失相干恢复和地形误差去除后的升降轨形变场。图3为升轨数据获取的失相干恢复后的巴姆地震同震形变场,可以发现失相干恢复后的结果保证了形变场的完整性。图4(a)为失相干恢复后的降轨形变场,可以看出,失相干恢复后,保持了形变场的完整性,但是形变场中存在明显的地形相关误差；图4(b)为存在地形误差的降轨数据缠绕相位；图4(c)为使用函数模型去除地形误差后降轨数据形变场；图4(d)为去除地形误差后缠绕相位。从图4(b)、图4(d)对比可以看出,去除地形误差前,存在地形误差的区域与周边区域存在很大的相位差。经过去地形误差处理后,整体形变场除中心主形变区外,周边区域的形变量都趋于0。从图4(d)可以看出，原本存在相位差的地形区域进行处理后相位差的变化。最终使用进行失相干恢复并去除地形误差后的降轨形变场[图4(c)]及升轨数据形变场(图3)进行后续三维分析处理。

图3 2003年12月26日巴姆Mw6.6地震升轨数据LOS向形变场Fig.3 Los direction deformation field of the orbit ascending data of Bam Mw6.6 earthquake on December 26, 2003

从D-InSAR获取的LOS向形变场结果[图3、图4(c)]可以看出,地震造成的形变中心呈蝴蝶状,定义LOS向形变场中靠近卫星飞行方向为正值,远离卫星飞行方向为负值。从升轨同震形变场(图3)可以看出,断层左侧有两个明显的形变中心,左上形变中心表现为正值,左下表现为负值,相比之下断层右侧形变中心量级较小,但依然能看出两个特征不同的形变中心；从降轨同震形变场图4(c)可以看出,升降轨同震形变场特征类似,整体表现为蝴蝶状,但是受升降轨特征影响,降轨中断层右侧形变特征表现的更明显,最大形变量近0.3 m。

1.2 MAI技术获取方位向形变

MAI技术[27]是利用雷达波束多普勒频移的正负,分裂波束的InSAR处理算法,从而获得目标地物方位向形变的信息。将两幅单视复数影像根据多普勒频移的符号,分裂成前后视两对主从影像,并分别进行干涉。由于前后视干涉图的空间基线几乎相等,又受到了同等的大气影像,而且它们探测到的距离向的形变也是相同的,所以将前视和后视的干涉图进行相位差分,得到方位向形变的相位信息。

图5为MAI技术获取的巴姆地震方位向形变场,受数据影响,升轨数据对形变场的覆盖并不完整,降轨数据可以完整覆盖形变场。从图5可以看出,升降轨数据获取的方位向形变场在断层两侧都有明显的特征,升轨中[图5(a)],断层左侧表现为正值,断层右侧表现为负值,降轨数据与之相反[图5(b)],最大方位向形变量为1 m。

图5 MAI技术获取的2003年12月26日巴姆Mw6.6地震方位向形变场Fig.5 Azimuth deformation field of Bam Mw6.6 earthquake on December 26, 2003 acquired by MAI Technology

2 三维形变场解算

LOS向形变结果和方位向形变结果与地表真实东西向、南北向、垂直向形变结果的关系表达式为

(2)

(3)

通过获取升降轨LOS向和方位向形变场,根据升降轨情况,可以得到4个方程,转化为矩阵形式可表示为

(4)

根据最小二乘原理,三维形变场的最优解为

X=(BTPB)-1BTPL

(5)

式(5)中：P为观测值方差组成的权阵。

在SAR形变场结果中,选择合适的窗口,计算窗口中心点标准差,从而估计出数据的近似标准差。

根据式(4)、式(5)方程,使用多种InSAR技术联合的三维地表形变解算方法,获取巴姆地震的三维形变场,由于受升轨数据的影响,获取的三维形变场结果主要覆盖断层右侧信息,左侧信息不完全。从三维解算结果(图6)可以看出,在东西向上[图6(a)],定义向东为正,向西为负,主要表现为3个形变中心,断层左侧有一个向西形变中心,断层右侧有2个形变中心,其中右上为向西形变,右下为向东形变,最大向西和向西形变量为0.25 m；在南北向上[图6(b)],定义向北为正,向南为负,地震造成的形变沿断层方向断层左侧为向北运动,右侧向南运动,最大向南形变量为0.4 m；垂直向上[图6(c)],定义向上为正,向下为负,主要表现为三个形变中心,断层左侧为向上抬升,断层右侧两个形变中心,右上为向下沉降,右上为向上抬升,最大抬升形变量为0.24 m,最大下沉形变量为0.2 m。从三维形变场结果可以看出,巴姆地震的发震断层主要表现为近南北向的右旋走滑运动,同时存在一定的东西向水平旋转运动。获得的三维形变场结果与洪顺英[28]、季灵运等[25]和Fialko等[29]的研究成果一致。

3 结论

(1)引入非线性支持向量机回归(SVR)对形变场中失相干区域进行恢复,保证了形变场的完整性,为后续研究分析提供基础。

(2)构建DEM与形变场经纬度之间函数模型,对形变场中存在的地形误差进行去除,提高了形变场的可靠性。

(3)联合D-InSAR获取的LOS向同震形变场与MAI技术获取的方位向形变场,使用多种InSAR技术联合的方法进行三维形变场解算,对地震造成的地面形变进行更可靠的分析,通过巴姆地震的三维形变场解算结果可以看出,巴姆地震发震断层为近南北向的右旋走滑断层。完整三维形变场的获取,为后续发震断层研究、地质构造分析提供理论基础。