地球物理学报  2020, Vol. 63 Issue (4): 1340-1356   PDF    
综合InSAR和应变张量估计2016年MW7.0熊本地震同震三维形变场
袁霜1, 何平1, 温扬茂2, 许才军2     
1. 中国地质大学(武汉)地球物理与空间信息学院, 湖北省近地表多尺度成像重点实验室, 武汉 430074;
2. 武汉大学测绘学院, 武汉 430079
摘要:利用合成孔径雷达(Synthetic Aperture Radar,SAR)影像提取与地质活动相关的三维地表形变场,对深入理解地质灾害的形成机制及其潜在灾害风险评估非常重要.目前,利用SAR影像的同震三维形变场提取主要利用单个像素点的多次观测构建观测方程,然后基于加权最小二乘(Weighted Least Squares,WLS)方法分解从而获得同震三维形变场,因此该方法缺乏对相邻像素点空间相关性的约束.考虑相邻同震位移点的应力连续性,研究学者提出了顾及大地测量应变张量和卫星形变观测的SAR同震三维形变场方法(Extended Simultaneous and Integrated Strain Tensor Estimation from geodetic and satellite deformation Measurements,ESISTEM).本文以2016年MW7.0熊本地震为例,收集了覆盖此次地震的ALOS-2卫星升降轨影像,利用传统差分InSAR(DInSAR)方法和子孔径雷达干涉测量(Multiple Aperture InSAR,MAI)方法分别对升降轨SAR影像对进行处理,得到视线向(LOS)形变和方位向形变,最后利用ESISTEM方法获取此次地震的三维同震形变场.此外,利用GPS和野外考察观测对本文的三维形变场结果进行结果精度分析.研究结果表明,与传统WLS方法相比,ESISTEM方法不仅能有效抑制奇异像素点对形变结果的干扰,同时对近断层的失相干信号能进行较好的恢复,更有助于解释地表破裂区的地震形变特征和掌握地震发生机制.本文确定的三维同震形变场结果显示主形变区发生在Futagawa断层中部和Hinagu断层最北端,最大水平位移为2 m,抬升为0.55 m.断层破裂以NE-SW走向的右旋走滑为主兼有部分正断成分.应变张量分析表明发震断层处受到了明显的收缩力和剪切力的作用.
关键词: InSAR      同震三维形变场      应变张量      熊本地震     
Integrated InSAR and strain tensor to estimate three-dimensional coseismic displacements associated with the 2016 MW7.0 Kumamoto earthquake
YUAN Shuang1, HE Ping1, WEN YangMao2, XU CaiJun2     
1. Hubei Subsurface Multi-scale Imaging Key Laboratory, Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China;
2. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China
Abstract: The three-dimensional (3D) surface displacements associated with different geohazard phenomenon derived from Synthetic Aperture Radar (SAR) images significantly improve our ability to understand their formation mechanisms and assess the potential disaster risk. At present, the most common approach for 3D coseismic deformation field derived from SAR images is using the Weighted Least Squares (WLS) method to resolve the 3D components based on more than three observations at each pixel point. However, this approach will be no constraints on the spatial correlation between the adjacent pixel points. Therefore, a new approach, i.e., the Extended Simultaneous and Integrated Strain Tensor Estimation from geodetic and satellite deformation Measurements (ESISTEM), has been proposed to reveal the 3D coseismic displacement by adding a constraint equation of stress continuity between the adjacent pixel points. In this study, we collect both the ascending and descending ALOS-2 images for the 2016 Kumamoto earthquake sequence, and then generate the displacement both in line of sight (LOS) and azimuth direction by using the traditional differential InSAR (DInSAR) method and the Multiple Aperture InSAR (MAI) method, respectively. After that, we utilize the ESISTEM method to recover 3D coseismic deformation field for this event. In addition, we verify the uncertainty of the 3D displacements with GPS and field work observations, respectively. In compare to the traditional WLS approach, our results show the ESISTEM approach can not only effectively reduce the disturbance of singular displacement, and recover the displacement in incoherence region, which is helpful to the coseismic 3D displacement characteristics analyzed. Our final 3D displacement map demonstrates that the main deformation zone occurs in the central section of the Futagawa fault and the northernmost part of the Hinagu fault, and the maximum horizontal and uplift displacement are 2 m and 0.55 m, respectively. The subsurface rupture is dominated by NE-SW-trending dextral strike-slip with minor normal-slip component. The strain tensor indicates that the seismic area sustains an obvious contraction and shear forces.
Keywords: InSAR    Three-dimensional coseismic displacements    Strain Tensor    Kumamoto earthquake    
0 引言

对于地震形变监测来说,真实的同震三维形变场能更准确解译地震形变特征和发震断层的几何学特征.目前基于合成孔径雷达(Synthetic Aperture Radar,SAR)数据观测,常用的解算三维形变场的策略有三种:(1)SAR卫星升降轨视线向(Line of Sight, LOS)+方位向观测量,其中LOS观测主要指差分干涉测量(Differential Interferometric Synthetic Aperture Radar, DInSAR),而方位向观测方法常用的有子孔径干涉测量(Multiple Aperture Interferometry,MAI)和像素偏移估计(Pixel Offset-Tracking,POT). Fialko等(2001)融合DInSAR+POT方法成功获取了1999年MW7.1 Hector Mine地震的同震三维形变场. Jung等(2011)融合DInSAR+MAI方法成功获取了2007年夏威夷Kilauea火山喷发引起的三维形变场.(2)仅升降轨多视向DInSAR融合分解,该策略由Wright等在2004年提出,王永哲等(2012)利用该方法解算了2009年MW6.3拉奎拉地震的同震三维形变场.由于SAR卫星一般为极轨右视模式,因此该方法在南北向分量的分辨率仍然不足,导致并不常用.目前日本空间局(JAXA)的ALOS-2卫星具有同时左右视观测能力,为该方法应用提供了可能(Morishita et al., 2016).(3)DInSAR+GPS,Gudmundsson等于2002年提出将DInSAR和GPS数据进行融合,利用模拟退火算法得到了冰岛Reykjanes Peninsula地区的连续三维形变速率场.Samsonov和Tiampo(2006)对Gudmundsson的融合模型进行了优化,提出了能量函数最小化的解析解法,并获取了美国南加州地区的地表三维形变速率场.但这类方法都需要将GPS观测值进行内插至DInSAR观测格网内,并且强烈依赖于GPS台站的数量和分布(Hu et al., 2014).上述三类算法主要基于单个像素点的多次观测构建观测方程,通过最小二乘算法分解获得同震三维形变场,所以缺乏对相邻像素点空间相关性的约束.Guglielmino等(2011)提出一种顾及大地测量应变张量和卫星形变观测的三维形变场分解(Simultaneous and Integrated Strain Tensor Estimation from geodetic and satellite deformation Measurements,SISTEM)方法,该方法避免对GPS进行内插处理,在分解过程,利用周围点的观测值和应变张量梯度关系恢复待求点的三维形变量.Luo和Chen(2016)考虑到SISTEM方法中的待求点在实际观测中同时存在着InSAR值,因此能够对待求点三维形变量提供附加约束,从而提出了扩展的SISTEM方法(Extended SISTEM, ESISTEM),并成功应用于尼泊尔2015年MW7.9 Gorkha地震的三维形变场提取.

2016年4月15日16时25分(UTC)日本熊本县发生了MW7.0地震,震中位于(32.791°N,130.754°E),震源深度10 km,而主震前几天内发生了两次MW>6前震(USGS)(也称为熊本地震序列).此次地震序列造成了100多人死亡,还引发了大面积山体滑坡和严重的桥梁建筑破坏(Hashimoto et al., 2017; Scott et al., 2018).熊本地震发震区地质构造背景复杂,引发多段地表破裂,吸引了研究学者的广泛关注.地震发生后,Shirahama等(2016)进行了野外调查观测,结果表明熊本地震序列发生在长期以右旋走滑为主的Futagawa和Hinagu断层上,如图 1所示,并且Futagawa断层兼有正断成分(Toda et al., 2016).另外此次地震造成了约34 km的地表破裂,破裂终止于日本九州岛的阿苏火山口(Aso caldera)西部(Lin and Chiba, 2017).Himematsu和Furuya(2016)以及Fukahata和Hashimoto(2016)利用InSAR技术反演了此次地震的同震滑动分布特征.Moya等(2017)以及Scott等(2018)利用激光雷达数据估算了该地震沿Futagawa断层发生的同震位移空间分布.郑傲等(2017)以及金明培等(2017)利用近场强震记录反演了此次地震的震源破裂过程.蒋生淼等(2018)利用远场体波和InSAR测量对此次地震进行有限断层联合反演.然而,对此次地震的同震三维形变场及近断层破裂研究仍显不足.

图 1 2016年熊本地震序列的地形图 五角星为主震震中位置;圆点为前震震中位置;大的沙滩球为主震;小的沙滩球为前震;三角形为阿苏火山(Aso caldera);方块为熊本的位置;实线表示断层线;实线框表示升降轨ALOS-2数据覆盖范围. Fig. 1 Topographic map of the 2016 Kumamoto earthquake sequence The star shows the epicenter of the main shock; the dots show the epicenter of the foreshocks; the large beachball shows the main shock; the small beachballs shows the foreshocks; the triangle show the Aso caldera; the square shows the location of Kumamoto; the solid lines show the faults; the solid boxes show the ascending and descending ALOS-2 data coverage, respectively.

为揭示此次地震精确的同震三维地表破裂,本文收集了ALOS-2卫星的升降轨SAR数据,分别利用DInSAR和MAI技术处理获取LOS向和方位向观测值,再利用ESISTEM方法来获取此次地震的同震三维形变场,并与基于单个像素点的最小二乘分解结果进行比较.同时,本文收集了6个同震GPS点观测(Yue et al., 2017)和野外考察观测(Shirahama et al., 2016)对本文结果进行了不确定度分析.最后,比较分析震间应变率场与同震三维形变场的联系.

1 研究区域和数据预处理 1.1 研究区构造背景

日本岛链位于欧亚板块、太平洋板块和菲律宾板块之间的交汇处,活跃的板块相对运动导致日本是世界上最容易发生地震的国家之一.现代大地测量观测表明太平洋板块和菲律宾板块分别以8~9 cm·a-1和5~6 cm·a-1的速度向欧亚板块俯冲(Imanishi et al., 2012Kato et al., 2016),图 2a是日本九州岛及周边区域在欧亚参考框架下的震前GPS速度场,数据来自Kreemer等(2014),可发现太平洋板块和菲律宾向欧亚板块俯冲时,在欧亚板块的阻挡下日本发生逆时针旋转的整体运动.此次地震发生在日本九州岛中部的熊本县(图 1),这是自2011年3月11日在日本境内发生9.0级大地震后,又一次7级以上的地震,也是九州岛内首次发生如此大规模的地震.图 2b是扣除九州岛整体刚性运动的GPS无旋转速度场,以Futagawa-Hinagu断层带走向为界,西北块体内的站点一致向东北方向运动,越靠近断层处运动速率越小,表现为块体应变积累过程.东南块体的东侧内的站点以北西向运动为主,靠近断层附近的站点运动速率减小,且运动方向转变为南西方向,表明块体内部存在挤压应变积累;东南块体南部内的站点向南东向运动,东南块体整体表现为逆时针旋转运动.西北块体和东南块体明显的差异运动,表明断层两侧附近的地块存在右旋走滑运动.熊本地震序列的两次前震发生在Hinagu断层的最北端,主震发生在Futagawa断层上,发震断层的东北侧是活跃的阿苏火山群.Futagawa-Hinagu断层系统是九州岛主要的活动断层系统之一,是一种右旋走滑断层系统.Hinagu断层和Futagawa断层的走向分别为205°和235°(Shirahama et al., 2016).此次地震在Futagawa断层破裂长度约为28 km,最大右旋错动约220 cm;发生在Hinagu断层最北端的破裂长度约6 km,最大右旋错动达到70~80 cm(Shirahama et al., 2016).

图 2 (a) 九州岛地区GPS速度场和(b)九州岛地区无旋转GPS速度场 Fig. 2 GPS velocity field (a) and after removing the whole rotating component (b) in Kyushu island area
1.2 实验数据预处理

本文研究利用日本空间局(JAXA)L波段的ALOS-2 PALSAR影像数据,成像模式为StripMap(分辨率3~10 m,幅宽70 km),包括一对升轨数据和一对降轨数据(表 1).与C波段或X波段SAR数据相比,波长更长的L波段数据有利于在植被茂密和地形起伏的山区保持较高的相干性(Kobayashi et al., 2011).另外表 1的升降轨数据显示了较短的时间和空间基线,也有利于保持高的相干性.表 1中的升降轨道像对都跨越了整个熊本地震序列,即同时包含两次MW>6前震和MW7.0主震的形变信息.

表 1 本次研究所采用的ALOS-2卫星数据 Table 1 The ALOS-2 data used in this study
1.2.1 LOS向同震位移

本节采用传统的二轨DInSAR方法和GAMMA干涉处理软件对PALSAR原始数据进行处理.将表 1中的两对SLC影像分别进行配准和重采样.然后将主辅影像进行复共轭相乘得到干涉图,并且利用1弧度分辨率的SRTM数据(Farr et al., 2007)去除地形相位.为降低干涉图噪声的影响,在数据处理过程中分别对数据进行了多视处理和Goldstein滤波(Goldstein and Werner, 1998)处理.最后,利用SNAPHU相位解缠方法(Chen and Zebker, 2000)进行解缠处理,并将干涉图地理编码至WGS-84坐标系.

利用DInSAR方法处理得到熊本地震主要破裂区域的LOS向同震形变场如图 3.图 3显示了形变场沿着近NE-SW走向的Futagawa-Hinagu断层迹线两侧展布.在靠近断层位置存在条带状失相干区,原因在于大的地表破裂造成形变梯度过大从而难以解缠.在升轨干涉图(图 3a)中,沿断层两侧形变性质相反,其中西北盘形变值为正(沿着视线方向抬升),最大形变值约2.4 m,东南盘形变为负(沿视线向沉降),呈现出典型走滑型地震在弹性形变场下产生的反对称特征(Segal, 2010胡俊等,2013).在降轨干涉图(图 3b)中,断层两侧并没有明显的形变性质相反特征,西北盘主要表现有两个沿视线向的沉降中心,在火山口附近表现为沿视线向抬升.虽然此次地震为走滑型地震,但是在降轨形变图中并没有表现出典型走滑型地震的反对称形变特征,主要是降轨卫星飞行方向几何与此次断层走向近似平行,导致形变信息在干涉图中的投影很小,反映了一维LOS向形变的局限性,不能准确解译地表真实的形变状况.一维LOS向形变uLOS与地表三维形变矢量的投影关系可简单描述如下:

图 3 熊本地震序列的同震LOS向(左列)和方位向(右列)形变场 (a, c)升轨形变场; (b, d)降轨形变场. Fig. 3 The coseismic displacements of the Kumamoto earthquake sequence in the LOS direction (left) and azimuth direction (right) (a, c) The ascending displacements; (b, d) The descending displacements.

(1)

其中,unueuv分别是南北、东西和垂向的地表形变分量;φ是卫星轨道的方位角;θ是卫星入射角;δuLOS是LOS向观测值的误差.

1.2.2 方位向同震位移

在本节中,利用MAI方法对表 1中升降轨数据进行处理以获得熊本地震序列的方位向形变场.MAI方法由Bechor和Zebker(2006)首次提出,作为DInSAR方法一维视线向局限性的补充.到目前为止,MAI已被应用于冰川运动、火山活动和地震等领域的方位向形变测量中(王晓文,2017).首先将表 1中的两对SLC影像分别完成配准和重采样,再进行带通滤波分割出前、后视子孔径SLC影像(Barbot et al., 2008);然后将前、后视子孔径SLC影进行干涉,得到前、后视干涉图,再对前、后视干涉图进行干涉得到初始MAI干涉图;其中基线误差导致MAI干涉图会出现平地和地形相位残余,我们利用多项式模型来消除这两种相位残余(Hu et al., 2014),最后将方位向形变进行地理编码得到熊本地震的方位向形变场.

本文获取的方位向形变场如图 3(c—d)所示.和LOS向形变场相似,升降轨方位向形变都沿着发震断层分布,并且在断层附近出现失相干现象,最大正负形变值分别约为2 m和-2.2 m.升降轨方位向形变在发震断层两侧均呈现出形变性质相反的特征.在升轨方位向形变图(图 3c)中,断层西北盘地表朝向卫星飞行的方向运动,断层东南盘远离卫星飞行的方向运动;除了断层附近,东南盘出现了失相干区域,这是由于MAI方法将全孔径SLC分割为子孔径SLC后,雷达回波信号方位向带宽受到了损失,所以相对全孔径SLC干涉得到的DInSAR测量,MAI测量中更容易受相位失相干的影响.在降轨方位向形变图(图 3d)中,断层西北盘地表远离卫星飞行的方向运动,断层东南盘朝向卫星飞行的方向运动.方位向形变uAZO与地表三维形变矢量的投影关系可简单描述如下:

(2)

其中,δuAZO是方位向观测值的误差.

2 ESISITEM原理

假设一次地震活动使得地球表面发生了变形,我们定义一个任意点P和其周围N个实验点EPs,P点的位置为(xeP, xnP, xvP),实验点EPs的位置和位移分别是(xe(i)EPs, xn(i)EPs, xv(i)EPs)和(uEPs, ei, uEPs, ni, uEPs, vi),其中i=1, 2, …, N.基于弹性理论,用实验点EPs去估计任意点P的位移(ue, un, uv)的方程(Guglielmino et al., 2011),即应变张量估计方程表示为

(3)

其中

(4)

(5)

(6)

其中,l是包含了12个参数的待求向量,其中前3个参数(ue, un, uv)为P点的三维形变分量,后9个参数是应变矢量,其中εω分别代表应变张量和刚体旋转张量;uEPs是由N个实验点EPs的位移组成的向量;AEPs是一个3N×12的线性矩阵,AEPs=[AEPs, e1, AEPs, n1, AEPs, v1, …, AEPs, eN, AEPs, nN, AEPs, vN]T,其中Δxj(i)=xj(i)EPs-xjP(j=e, n, v; i=1, 2, …, N)表示第i个EPs点与任意点P在j方向上的距离.

假设在DInSAR观测下,存在实验点Q在LOS向位移uQ, LOS上的单位矢量为(SeQ, SnQ, SvQ),位置和位移分别是(xeQ, xnQ, xvQ)和(uQ, e, uQ, n, uQ, v),其中i=1, 2, …, N.

实验点Q的位移(uQ, e, uQ, n, uQ, v)与P点三维位移之间的关系为

(7)

实验点Q的LOS向位移uQ, LOS与三维位移之间的关系为

(8)

把式(7)代入式(8)有

(9)

上式表示实验点Q的LOS向位移对P点的约束,当有N个实验点EPs的LOS向位移作为约束时,约束方程为

(10)

其中

(11)

(12)

由于DInSAR技术仅能测量沿卫星LOS向的位移(Wright et al., 2004),导致进行三维形变场分解时对南北向的约束不足,因此需要补充方位向观测.本文利用DInSAR和MAI观测的位移矢量来构建三维位移场.同上述DInSAR观测类似,在MAI观测下,N个实验点EPs的方位向位移uEPs, AZO对P点的约束方程为

(13)

如果在P点上存在DInSAR和MAI观测,则可以将P点的DInSAR和MAI观测位移(uLOSPuAZOP)引入到应变张量估计中建立新的观测方程(Guglielmino et al., 2011),则有

(14)

其中uLOSP和的uAZOP单位向量分别为(SeLOS, SnLOS, SvLOS)和(SeAZO, SnAZO, SvAZO).上式右侧系数阵为

(15)

结合方程(10)、(13)和(14),可以得到

(16)

将上式简写为

(17)

因此,任意点P的三维位移场是式(17)的最小二乘解,即

(18)

其中,W是观测量权阵,本文定权方法基于实验点EPs到任意点P的距离和观测数据误差.待求的未知向量l除包括三维形变位移外,还包含了6应变张量(ε11, ε12, ε13, ε22, ε23, ε33)和3个刚体旋转张量(ω1, ω2, ω3)参数.因此,ESISTEM方法可以获得三个应变场分量,即体积胀缩应变分量σ、差分旋转幅度Ω2和最大剪切应变M,分别表示为

(19)

(20)

(21)

其中,λmaxλmin是应变张量的最大和最小特征值(Guglielmino et al., 2011).

Guglielmino等(2011)在实验中得到均匀分布的实验点EPs对点P位移的估计性能最佳.由于InSAR测量与三维形变分量的线性叠加关系,在本文中,优先选择均匀分布的实验点EPs的观测位移进行约束.周围约束点EPs的数量太少则会使模型病态,而数量过多时,性能没有明显的改善.Guglielmino等(2011)分析了水平位移和垂直位移的均方根误差与EPs数量的关系,得出超过60个EPs点后均方根误差开始趋向平缓.

3 三维同震形变场

利用上述的试验数据和ESISTEM方法,即可得到此次地震的三维同震形变位移.此外为了比较分析,本文同时利用传统基于像素点分解方法进行三维形变位移分解.在本文的分解之前,所有的像素点格网化为100 m×100 m的单元大小.

3.1 ESISTEM的三维同震形变场

ESISTEM方法解算过程中,经过测试选择N=100作为适当约束点数量阈值(Luo and Chen, 2016).图 4(a—c)是ESISTEM方法解算的同震三维形变场,可看出在地表形变梯度过大的失相干区域,仍然能够得到连续、有效的形变信息.图 4a东西形变分量主要表现为西北盘向东运动,最大形变区位于断层中部,最大形变量约2 m.图 4b南北形变分量表现为南北方向的扩张性拉伸,断层西北盘主要呈现出两个北向运动的形变中心,最大形变量为0.9 m.断层东南盘为南向运动,最大位移约1.2 m.图 4c垂直形变分量显示断层西北盘最大沉降约2 m,主要形变中心在断层中部,东南盘最大隆升约0.55 m.阿苏火山的最大升降约50 cm,沉降可能是此次右旋地震向西位移,阿苏火山岩浆室运移的结果(Himematsu and Furuya, 2016).另外在阿苏火山口北部,存在一个局部异常变形信号(黑框),在东西向形变中表现为向西运动,约0.4 m,在南北向形变中表现为北运动,最大位移2 m,垂向上为轻微的正位移.Himematsu和Furuya(2016)研究分析表明不能用同震山体滑坡解释,可能是小断层导致的局部变形,这些小断层被阿苏火山厚厚的火山灰掩盖.

图 4 熊本地震序列的三维形变场和残差 E-W、N-S和U-D分别表示东西,南北和垂直方向的位移方向;(a)—(c) ESISTEM方法的三维形变场;(d)—(f)WLS方法的三维形变场;(g)—(i)ESISTEM与WLS方法的三维形变场残差;箭头为水平位移矢量;黑色虚线表示4条剖面;虚线黑框为局部异常信号;(a)中黑色方块为GPS台站;(g)—(i)中实线黑框(S1—S6)用于后面3.4节做不确定度分析. Fig. 4 The three-dimensional displacements and residuals of the Kumamoto earthquake sequence The E-W, N-S and U-D indicate displacement directions in the east-west, north-south, and vertical directions, respectively; (a)—(c) The three-dimensional displacements of the ESISTEM method; (d)—(f) The three-dimensional displacements of the WLS method; (g)—(i) The residuals of displacements recovered by the ESISTEM and WLS methods; the arrows indicate horizontal displacement vectors; the black dashed lines indicate four displacement profiles; the black dashed boxes indicate the patchy localized signals; In (a), the black squares are GPS stations; In (g)—(i), the black solid boxes(S1—S6) used for later uncertainty analysis in the section 3.4.

综上所述,ESISTEM方法确定的三维同震形变场中(图 4c),水平位移分量在发震断层两侧存在明显的差异运动,西北盘一致向东北方向运动,东南盘向南东向运动,表明了2016年熊本地震序列以NE-SW向右旋走滑为主,这与震前GPS无旋转速度场结果一致;南北盘横向错动量和走滑量几乎一致,表明断层倾角近乎垂直;Futagawa断层西北盘受南部地块俯冲导致地表沉降,东南盘轻微隆升,表明Futagawa断层运动以右旋走滑为主兼有正断成分.

3.2 基于单像素点分解的三维同震形变场

利用1.2节中的升降轨LOS向形变和升降轨方位向形变位移组成观测方法,即联合式(1)和(2)式,然后利用加权最小二乘可解算出每个像素点的三维位移(Fialko et al., 2001He et al., 2018),即传统WLS方法(本文中WLS方法即指基于单像点的三维形变分解),表达式可简写为

(22)

其中和观测量Σ是观测量的协方差矩阵,σdLOSσdAZ是LOS向和方位向形变的误差.

本节利用WLS方法获取了2016年熊本地震序列的三维形变场(图 4(d—f)).从结果图可看出WLS法和ESISTEM方法获取的三维形变场变化趋势基本一致,表现为沿断层两侧相反方向的错动,并且主要形变区在震中Futagawa-Hinagu断层带附近.在图 4(d—f)中,可看出三维形变结果中存在大量的失相干区域,特别是在震中附近和破裂断层两侧,限制了对地震形变特征的解释.除了失相干区域,图 4d东西向形变分量主要表现为向东运动,断层北部向东最大位移约2 m,断层南部向西最大位移约0.9 m;图 4e南北向形变分量显示东南盘为南向运动,最大位移约1.3 m,断层西北盘主要呈现出两个北向运动的形变中心,最大形变量为1 m;图 4f中黑色箭头表示水平位移,颜色表示垂直形变场,垂向上断层北部最大沉降约2 m,南部最大隆升约0.6 m,这与Scott等(2018)的结果一致.

图 4(g—i)是ESISTEM与WLS的三维形变残差结果,可得出ESISTEM方法恢复了断层附近以及东南盘的失相干信息,断层附近地表破裂带形变信息的恢复,有利于分析发震断层机制和解释地震形变特征.除了失相干区域,其他区域的形变近似一致(残差都小于10-3 m).

3.3 三维同震形变场特征比较

此次地震事件以Futagawa断层为地表主要形变中心,另外还有Hinagu断层最北端和阿苏火山西部两个形变区.为了研究跨断层三维同震形变场的特征,我们沿垂直断层地表破裂走向做了四条剖面(图 4a—c)来进行分析,剖面结果如图 5.从同震形变场剖面图(图 5)可看出此次地震的最大形变都发生在断层两侧,其中最大形变区在剖面C-C′附近,即Futagawa断层中部,并且对应地表发生了明显错动.另外,ESISITEM方法(灰圆点)对应的形变场剖面都呈现了良好连续性,WLS方法(黑圆点)对应的形变场剖面出现了失相干现象.总体而言,二种方法的形变场剖面趋势一致.

图 5 ESISTEM和WLS的三维形变场分别沿着A-A′、B-B′、C-C′和D-D′剖面的比较结果 Fig. 5 Comparison of the three-dimensional displacements in the three different components (E-W, N-S, U-D) along the profiles A-A′, B-B′, C-C′, D-D′ derived from the ESISTEM (gray points) and WLS (black points) methods

A-A′剖面位于阿苏火山西部,可看出断层两侧的水平位移明显大于垂直位移,与此次地震以走滑为主的运动性质一致;WLS方法东西向形变在跨断层两侧连续,表明地表无明显破裂.B-B′剖面的东西向和南北向形变以断层为中心对称分布,形变主要集中在断层附近10 km的范围内,且断层两侧形变衰减极快;垂向形变在断层10 km的范围内也发生了较快的衰减,西北盘地表沉降量远大于东南盘的隆升量.C-C′剖面位于主断裂的中部,垂向形变与B-B′剖面中表现一致,沉降远大于隆升,表明了断层地表以沉降运动为主.C-C′和D-D′剖面的东西向和南北向形变与B-B′剖面相同,且断层两侧地表都发生明显错动和破裂.WLS方法和本文方法在四条形变剖面上都很吻合,但本文方法比WLS方法平滑和连续.

3.4 三维形变场不确定度估计

本节对ESISTEM方法和WLS方法获取的三维形变场结果进行不确定度分析,将从两个方面进行分析,一是计算这两种方法结果残差的均值和标准差来进行精度分析,二是利用GPS数据和野外考察资料进行结果的准确度评估.

3.4.1 精度分析

本小节根据离发震断层的距离远近,分别在断层两侧选取三个区域来计算同震三维形变场残差(图 3gi)的均值和标准差,结果如表 2所示,每个区域的大小为20×20像素.西北盘于垂直断层方向由近到远选取S1、S2和S3三个区域,东南盘垂直断层方向由近到远选取S4、S5和S6三个区域(图 3gi中的实线黑框),断层两盘的三个区域关于断层大致对称,且都选在影像相干区域.从表 2可看出西北盘区域(S1—S3)分别在三个形变场分量上的均值相同,标准差几乎没有差异且都小于1 cm;东南盘区域(S4—S6)分别在三个形变场分量上的均值和标准差都几乎没有差异.这些统计结果表明了三维形变场残差的精度跟断层远近没有关系,ESISTEM方法和WLS方法在影像相干区域的效果具有高的吻合度,也证明了ESISTEM方法的合理性.对于盲破裂区域的研究这两种方法都可以成功获取地表三维形变场,区别在于当地表发生在高形变失相干区域时ESISTEM方法能有效地恢复失相干区域信息.

表 2 同震三维形变场残差比较 Table 2 Residual comparison of three-dimensional coseismic displacements
3.4.2 准确度分析

除比较WLS方法和ESISTEM方法获取的同震三维形变场远场内符合精度外,本文同时利用了此次地震的GPS同震位移和野外考察资料进行结果的外符合精度(即准确度)评估.

本文首先收集了同震三维形变场覆盖区域的6个GPS站点数据(Yue et al., 2017)(如图 4d),与本文两种方法构建的三维形变位移比较分析.图 6(a—b)分别显示了GPS测量、WLS方法和ESISTEM方法的水平和垂向位移矢量,它们对应的位移残差如图 6(c—d).位移残差结果显示水平位移在西北方向上表现出系统性偏移,垂向位移的残差很小,总体来说GPS测量和本文两种测量方法的位移有良好的一致性.表 3给出了GPS测量与本文两种方法在三维形变场上的RMSEs(均方根误差),其中GPS与WLS结果之间的RMSEs在东西向、南北向和垂向上分别为14.26、20.92和8.84 cm,与ESISTEM方法结果之间对应的RMSEs在东西向、南北向和垂向上分别为14.23、19.81和8.54 cm.从比较结果看出,ESISTEM方法在三个方向上的结果精度都比WLS方法的要好,表明了ESISTEM方法更能解算出近场的形变特征.但是总体来说与GPS验证站的结果差异比较大,可能是因为不够精确的InSAR测量和GPS观测值先验方差所引起的.

图 6 InSAR测量和GPS测量在GPS站点上的同震三维位移比较 水平位移(a)和垂直位移(b)来自GPS测量(黑色),ESISTEM方法(红色)和WLS(蓝色);InSAR测量与GPS测量在水平(c)和垂直(d)向上的位移残差. Fig. 6 Comparison of three-dimensional coseismic displacements between InSAR measurements and GPS measurements at the corresponding GPS sites Horizontal (a) and vertical (b) displacements from GPS(black), ESISTEM(red) and WLS (blue) methods; residual displacements between InSAR measurements and GPS measurements in horizontal (c) and vertical (d) directions.
表 3 与GPS和野外调查测量对比的三维位移精度分析(单位:cm) Table 3 Accuracy analyzed for three-dimensional displacement component compared with GPS and field work measurement (unit: cm)

另外我们收集了来自震区的野外调查结果(Shirahama et al., 2016),如图 7所示.图 7a是Futagawa和Hinagu断层上表面破裂的水平位错,在Futagawa断层中部测量到大部分水平位错量约为2.0 m,最大水平位错为2.2 m,并沿着西南和东北两端水平位错量逐渐减小.图 7f是Futagawa断层上表面破裂的垂直位错,在Futagawa断层西南部和中部南侧以抬升为主,东北部以沉降为主,Futagawa断层中部测量到最大抬升1.2 m.图 7b7c分别是本文WLS和ESISTEM对应的水平位错,整体可看出ESISTEM结果(c)与野外调查结果(a)更相似,如(c)中黑框区域,且ESISTEM恢复了破裂断层上的失相干水平位错;(d)、(e)分别是(b)、(c)与野外调查结果(a)的残差,可看(c)的残差明显小于(b)的残差.图 7(g)(h)分别是WLS和ESISTEM对应的垂直位错,经过对比得出ESISTEM能恢复破裂断层上的失相干垂直位错;(i)、(j)分别是(g)、(h)与野外调查结果(f)的残差,整体比较可看出(h)的残差小于(g)的残差.表 3给出了本文WLS和ESISTEM与野外调查结果在水平和垂直位错上的RMSEs.WLS方法与野外调查在水平和垂直位错上的RMSEs分别是46.33和47.11 cm,ESISTEM与野外调查在水平和垂直位错上的RMSEs分别是18.94和25.65 cm.本文两种方法与野外调查结果对比得出ESISTEM能够对近断层的失相干信号进行一定程度的恢复,并且结果比WLS方法结果更接近真实的地表三维形变场.

图 7 Futagawa-Hinagu断层带(灰线)表面破裂的水平位错和垂直位错 (a)、(b)和(c)是分别来自野外调查,WLS方法和ESISTEM方法的水平位错;(d)和(e)是野外调查分别与WLS和ESISTEM的水平位错残差;(f)、(g)和(h)是分别来自野外调查,WLS方法和ESISTEM方法的垂直位错;(i)和(j)是野外调查分别与WLS和ESISTEM垂直位错的残差;实线黑框是对比点位. Fig. 7 Horizontal and vertical offsets of surface rupture along the Futagawa-Hinagu fault systems (gray lines) (a), (b), (c) are the horizontal offsets from field work, WLS method, and ESISTEM method, respectively; (d) is the residual of (a) and (b); (e) is the residual of (a) and (c); (f), (g), (h) are the vertical offsets from field work, WLS method, and ESISTEM method, respectively; (i) is the residual of (f) and (g); (j) is the residual of (g) and (h); the black solid boxes are the contrast point.
4 讨论

在地震形变监测中,真实的同震三维形变场能更准确掌握地震发生机制和地震形变特征.常用的三维形变测量手段有地表重复点观测技术(如全球导航定位技术和水准测量等)以及机载激光雷达(LiDAR)技术和InSAR技术等.GPS技术可以获取连续的地表形变信息,但全球连续观测站分布仍然有限,其密度受限于接收机设备和维护的高昂成本,导致空间分辨率低(>10 km);目前美国南加州的SCIGN网和日本的GEONET网是世界上最密集的GPS监测网(Samsonov and Tiampo, 2006).LiDAR技术能够快速获取高精度的地表三维形变,但是它的费用昂贵,并且相对InSAR技术来说覆盖面积小,例如Moya等(2017)利用LiDAR数据也获取了此次熊本地震的三维形变场,覆盖面积为80 km2.InSAR技术具有高空间分辨率(几十到几米)、高精度(厘米到毫米)、低成本、全天候和观测面积广(几十到几百公里)等巨大优势,所以InSAR技术是获取同震三维形变最重要的手段之一.

目前利用InSAR监测同震三维形变场最常用的方法是融合DInSAR技术和方位向观测技术,然后通过加权最小二乘算法重建三维形变场,即传统WLS方法.但是WLS方法缺乏对相邻像素点空间相关性的约束,所以对InSAR监测中的失相干区域无能为力.Guglielmino等(2011)的SISTEM方法,融合DInSAR和GPS技术,并且加入了应变张量估计模型,不仅避免了对GPS进行内插,还有效地恢复了三维形变场中的失相干信息;但是仅考虑了GPS点与InSAR观测点之间的约束关系,忽视了InSAR观测点之间的相互约束.后面Luo和Chen(2016)在SISTEM方法的基础上考虑了InSAR观测点之间的相互约束,提出了ESISTEM方法.本文仅基于ALOS-2卫星的SAR数据,利用ESISTEM方法来获取此次熊本地震序列的三维形变场.日本JAXA的ALOS-2卫星SAR数据能够完整地涵盖熊本地震的形变区域,并且ALOS-2卫星的L波段更有利于监测此次地震严重的、长达34 km的地表破裂情况.虽然此次地震已经被多种观测手段进行了研究,例如Shirahama等(2016)进行了野外调查观测,Moya等(2017)以及Scott等(2018)利用了LiDAR技术,Himematsu和Furuya(2016)利用POT技术获取了同震三维形变场,但是本文基于SAR数据并顾及了应变张量估计,对近场失相干区域进行了恢复和识别;另外本文的研究结果方便与已有的野外调查数据以及GPS数据进行比较验证.

为了直观地比较分析ESISTEM方法对熊本地震三维形变场的解算能力,将其与传统WLS方法同时进行解算,比较了这两种方法结果残差的标准差,并利用野外调查数据和GPS数据验证了该方法的可靠性和精度.残差标准差对比得出ESISTEM方法与WLS方法在影像相干区域的精度相当,对于盲破裂研究区域,这两种方法都可以有效地获取地表三维形变场,例如Luo和Chen(2016)利用ESISTEM方法成功获取了2015年MW7.9尼泊尔Gorkha地震高精度的三维形变场,据报道此次地震中断层没有破裂到地表(Feng et al., 2015);He等(2018)利用WLS方法获取了2017年伊朗Ezgeleh地震盲断层的高精度三维形变场.但是对于地表发生破裂的失相干区域,ESISTEM方法能有效地恢复失相干区域的形变信息,更利于对同震三维形变场的特征进行分析.本文收集的GPS站点数据都位于影像相干区域内,从GPS测量与这两种方法在三维形变场上的RMSEs结果(表 3)看出,在相干区域内ESISTEM方法能更有效地解算出近场的形变特征.收集的野外调查结果位于近断层处,由于发生了地表破裂,所以包含了影像失相干区域,从野外调查与这两种方法获取的水平和垂直位错对比(图 7)和RMSEs结果(表 3),表明了ESISTEM方法在相干区域的精度比WLS方法的高,并且能恢复出可靠的近断层失相干区域信息.

本文利用ESISTEM方法除了得到同震三维形变场外,还可以得到三个无量纲的应变场分量(图 8),即体积胀缩应变、差分旋转幅度和最大剪切应变.首先这三个应变场直观地反映了发震断层的位置,和断层为NE-SW走向.图 8a显示了此次地震的体积胀缩应变分量,揭示了体内的膨胀和收缩;(图 8a)中膨胀系数在断层处呈现为负值,最大数值约2.8×10-1,表明断层处发生收缩,由于断层附近西北盘向北东向运动,东南盘向南西向运动,所以断层近场地区发生膨胀现象.差分旋转幅度如图 8b,在Futagawa断层东部差分旋转幅度最大值约9.4×10-2.图 8c是最大剪切应变分量,在断层处及其附近剪切作用明显,最大值约9.8×10-1.

图 8 熊本地震序列的应变场 (a)体积胀缩应变分量;(b)差分旋转幅度应变分量;(c)最大剪切应变分量. Fig. 8 The strain field of the Kumamoto earthquake sequence (a) Dilatation strain component; (b) Differential rotation magnitude; (c) Maximum shear strain component.

地壳构造背景应变率场与浅源地震活动的发生被认为有着重要的联系(如Riguzzi et al., 2012; Bird and Kreemer, 2015).本文收集了日本九州岛的震前GPS观测数据(图 2),基于多尺度球面小波方法(Tape et al., 2009)解算区域的构造背景应变率场.图 9给出了主应变率、第一剪应变率和面应变率场结果.结果显示九州岛地区主压应变率方向从北到南由北西向逐渐转为南西向,发震断层附近以近东西向的主压应变率为主,这与ESISTEM方法得到的体积胀缩应变分量结果一致,主压应变率高值区分布在发震断层的北东区域,造成主压应变率分布的主要构造动力来源于菲律宾板块的俯冲挤压;九州岛地区南部则以北西向的主张应变率为主,与北部地区的主应变率场呈现差异.图 9a为第一剪应变率场,由于此次地震发震断层的走向近北东向,选择第一剪应变率能更有效地反映该地块的剪切变形.图 9a中整体呈现为负值,只有最西南部表现为弱的正值,而第一剪应变率负值代表北东向右旋剪切变形,正值代表北东向左旋剪切变形(江在森等,2003),所以九州岛地区以右旋剪切变形为主,最显著的右旋剪切变形发生在发震断层的北东侧.图 9b为面应变率场,可看出九州岛地区整体表现为面压缩特征,只有最西南部表现为面膨胀特征,面压缩率最低值位于发震断层的北东侧.本文ESISTEM方法获取的同震位移场和多尺度球面小波方法解算的震前GPS速度场都显示了熊本地震发震区域存在显著的运动差异,且ESISTEM方法获取的应变场和多尺度球面小波方法解算的应变率场都显示了熊本地震发震区域位于剪切应变和面应变作用强的区域附近,这与江在森等(2006)的研究得出强震分布区通常存在于显著的地壳差异运动带,尤其是与区域主干断裂的构造活动背景相一致的剪应变率高值区或其边缘的结论相符.

图 9 九州岛地区GPS应变率分布 (a)第一剪应变率;(b)面应变率;十字线为主应变率. Fig. 9 The strain rate field from GPS in Kyushu island area (a) The first shear strain rate; (b) Dilatation strain rate; the cross indicate principal strain rate.
5 结论

InSAR观测在地震同震形变场测量中占据重要地位,其中同震三维地表形变场对解释地震形变特征和掌握地震发生机制具有重要意义.由于地震同震形变场具有形变幅度大、形变梯度大和构造背景复杂等特点,导致InSAR观测在震中区域和主要破裂区域会出现严重的相位失相干现象,这时利用WLS方法来构建同震三维地表形变场仍然对失相干严重区域无能为力.本文基于SISTEM和ESISTEM方法的研究,成功地获取了2016年MW7.0熊本地震的同震三维地表形变场和三个应变场分量,并且将ESISTEM与WLS方法进行比较,表明ESISTEM能获取可靠的同震三维地表形变场,并且能够对近断层的失相干信号进行一定程度的恢复.我们的研究结果表明熊本地震主要发生在以右旋走滑为主的Futagawa断层中部和Hinagu断层最北端,其中Futagawa断层兼有少量正断成分,地震走向为NE-SW,断层倾角近乎垂直.东西向形变主要表现为西北盘的东向运动,最大形变量约2 m;南北向表现为南北方向的扩张性拉伸,两盘的横向错动量和走滑量几乎一致;垂向形变中西北盘沉降高达2 m,东南盘最大抬升0.55 m.而ESISTEM方法得到的应变场分量表明破裂断层处受到了明显的收缩力和剪切力的作用.同震形变场和应变场结果表明,Futagawa断层破裂是右旋剪切兼少量正断的应变释放结果.利用多尺度球面小波方法解算震前GPS速度场和应变率场,得到发震断层两侧存在明显的差异运动,显示为右旋走滑的运动特征,并且熊本地震发生在第一剪应变率负高值和面压缩率最低值区域的边缘.ESISTEM方法获取的同震位移场和应变场与多尺度球面小波方法获取的GPS速度场和应变率场都反映了熊本地震发生在有明显差异运动、剪切应变和面应变作用强的区域附近.

致谢  日本宇航探测研究机构JAXA提供PALSAR-2数据(ID:3048),GMT(Generic Mapping Tools)软件用于绘图.感谢罗海鹏博士提出的宝贵建议和帮助.
References
Barbot S, Hamiel Y, Fialko Y. 2008. Space geodetic investigation of the coseismic and postseismic deformation due to the 2003 MW7.2 Altai earthquake:implications for the local lithospheric rheology. Journal of Geophysical Research:Solid Earth, 113(B3): B03403. DOI:10.1029/2007JB005063
Bechor N B D, Zebker H A. 2006. Measuring two-dimensional movements using a single InSAR pair. Geophysical Research Letters, 33(16): L16311. DOI:10.1029/2006GL026883
Bird P, Kreemer C. 2015. Revised tectonic forecast of global shallow seismicity based on version 2.1 of the Global Strain Rate Map. Bulletin of the Seismological Society of America, 105(1): 152-166. DOI:10.1785/0120140129
Chen C W, Zebker H A. 2000. Network approaches to two-dimensional phase unwrapping:intractability and two new algorithms. Journal of the Optical Society of America A, 17(3): 401-414. DOI:10.1364/JOSAA.17.000401
Farr T G, Rosen P A, Caro E, et al. 2007. The shuttle radar topography mission. Reviews of Geophysics,, 45(2): RG2004. DOI:10.1029/2005RG000183
Feng G C, Li Z W, Shan X J, et al. 2015. Geodetic model of the 2015 April 25 MW7.8 Gorkha Nepal earthquake and MW7.3 aftershock estimated from InSAR and GPS data. Geophysical Journal International, 203(2): 896-900. DOI:10.1093/gji/ggv335
Fialko Y, Simons M, Agnew D. 2001. The complete (3-D) surface displacement field in the epicentral area of the 1999 MW7.1 Hector Mine earthquake, California, from space geodetic observations. Geophysical Research Letters, 28(16): 3063-3066. DOI:10.1029/2001GL013174
Fukahata Y, Hashimoto M. 2016. Simultaneous estimation of the dip angles and slip distribution on the faults of the 2016 Kumamoto earthquake through a weak nonlinear inversion of InSAR data. Earth, Planets and Space, 68: 204. DOI:10.1186/s40623-016-0580-4
Goldstein R M, Werner C L. 1998. Radar interferogram filtering for geophysical applications. Geophysical Research Letters, 25(21): 4035-4038. DOI:10.1029/1998GL900033
Guglielmino F, Nunnari G, Puglisi G, et al. 2011. Simultaneous and integrated strain tensor estimation from geodetic and satellite deformation measurements to obtain three-dimensional displacement maps. IEEE Transactions on Geoscience and Remote Sensing, 49(6): 1815-1826. DOI:10.1109/TGRS.2010.2103078
Hashimoto M, Savage M, Nishimura T, et al. 2017. Special issue "2016 Kumamoto earthquake sequence and its impact on earthquake science and hazard assessment". Earth, Planets and Space, 69: 98.
He P, Wen Y M, Xu C J, et al. 2018. High-quality three-dimensional displacement fields from new-generation SAR imagery:application to the 2017 Ezgeleh, Iran, earthquake. Journal of Geodesy, 99(4): 573-591.
Himematsu Y, Furuya M. 2016. Fault source model for the 2016 Kumamoto earthquake sequence based on ALOS-2/PALSAR-2 pixel-offset data:evidence for dynamic slip partitioning. Earth, Planets and Space, 68: 169. DOI:10.1186/s40623-016-0575-1
Hu J, Li Z W, Zhang L, et al. 2012. Correcting ionospheric effects and monitoring two-dimensional displacement fields with multiple-aperture InSAR technology with application to the Yushu earthquake. Science China Earth Sciences, 55(12): 1961-1971. DOI:10.1007/s11430-012-4509-x
Hu J, Li Z W, Ding X L, et al. 2013. Derivation of 3-D coseismic surface displacement fields for the 2011 MW9.0 Tohoku-Oki earthquake from InSAR and GPS measurements. Geophysical Journal International, 192(2): 573-585.
Hu J, Li Z W, Ding X L, et al. 2014. Resolving three-dimensional surface displacements from insar measurements:a review. Earth-Science Reviews, 133: 1-17. DOI:10.1016/j.earscirev.2014.02.005
Imanishi K, Ando R, Kuwahara Y. 2012. Unusual shallow normal-faulting earthquake sequence in compressional northeast Japan activated after the 2011 off the Pacific coast of Tohoku earthquake. Geophysical Research Letters, 39(9): L09306. DOI:10.1029/2012GL051491
Jiang S M, Yi L, Zhang X, et al. 2018. Joint inversion of teleseismic and co-seismic InSAR data for the rupture process of the 2016 Kumamoto earthquake in Japan. Acta Seismologica Sinica (in Chinese), 40(1): 13-23.
Jiang Z S, Ma Z J, Zhang X, et al. 2003. Horizontal strain field and tectonic deformation of China mainland revealed by preliminary GPS result. Chinese Journal of Geophysics (in Chinese), 46(3): 352-358.
Jiang Z S, Yang G H, Wang M, et al. 2006. On crustal movement in China continent and its relationship with strong earthquakes. Journal of Geodesy and Geodynamics (in Chinese), 26(3): 1-9.
Jin M P, Li Z L, Wang R J. 2017. Coseismic displacement field and slip model derived from near-source strong motion records of MW7.0 Kumamoto, Japan, earthquake. Acta Seismologica Sinica (in Chinese), 39(6): 819-830.
Jung H S, Lu Z, Won J S, et al. 2011. Mapping three-dimensional surface deformation by combining multiple-aperture interferometry and conventional interferometry:Application to the June 2007 eruption of Kilauea Volcano, Hawaii. IEEE Geoscience and Remote Sensing Letters, 8(1): 34-38.
Kato A, Nakamura K, Hiyama Y. 2016. The 2016 Kumamoto earthquake sequence. Proceedings of the Japan Academy, Series B, 92(8): 358-371. DOI:10.2183/pjab.92.359
Kobayashi T, Tobita M, Nishimura T, et al. 2011. Crustal deformation map for the 2011 off the Pacific coast of Tohoku Earthquake, detected by InSAR analysis combined with GEONET data. Earth, Planets and Space, 63(7): 621-625. DOI:10.5047/eps.2011.06.043
Kreemer C, Blewitt G, Klein E C. 2014. A geodetic plate motion and Global Strain Rate Model. Geochemistry, Geophysics, Geosystems, 15(10): 3849-3889. DOI:10.1002/2014GC005407
Lin A M, Chiba T. 2017. Coseismic conjugate faulting structures produced by the 2016 MW7.1 Kumamoto earthquake, Japan. Journal of Structural Geology, 99: 20-30. DOI:10.1016/j.jsg.2017.05.003
Luo H P, Chen T. 2016. Three-dimensional surface displacement field associated with the 25 April 2015 Gorkha, Nepal, earthquake:solution from integrated InSAR and GPS measurements with an extended SISTEM approach. Remote Sensing, 8(7): 559. DOI:10.3390/rs8070559
Morishita Y, Kobayashi T, Yarai H. 2016. Three-dimensional deformation mapping of a dike intrusion event in Sakurajima in 2015 by exploiting the right-and left-looking ALOS-2 InSAR. Geophysical Research Letters, 43(9): 4197-4204. DOI:10.1002/2016GL068293
Moya L, Yamazaki F, Liu W, et al. 2017. Calculation of coseismic displacement from lidar data in the 2016 Kumamoto, Japan, earthquake. Natural Hazards and Earth System Sciences, 17(1): 143-156. DOI:10.5194/nhess-17-143-2017
Riguzzi F, Crespi M, Devoti R, et al. 2012. Geodetic strain rate and earthquake size:New clues for seismic hazard studies. Physics of the Earth and Planetary Interiors, 206. DOI:10.1016/j.pepi.2012.07.005
Samsonov S, Tiampo K F. 2006. Analytical optimization of a DInSAR and GPS dataset for derivation of three-dimensional surface motion. IEEE Geoscience and Remote Sensing Letters, 3(1): 107-111.
Segall P. 2010. Earthquake and Volcano Deformation. New York: Princeton University Press.
Scott C P, Arrowsmith J R, Nissen E, et al. 2018. The M7 2016 Kumamoto, Japan, earthquake:3-D deformation along the fault and within the damage zone constrained from differential lidar topography. Journal of Geophysical Research:Solid Earth, 123(7): 6138-6155. DOI:10.1029/2018JB015581
Shirahama Y, Yoshimi M, Awata Y, et al. 2016. Characteristics of the surface ruptures associated with the 2016 Kumamoto earthquake sequence, central Kyushu, Japan. Earth, Planets and Space, 68: 191. DOI:10.1186/s40623-016-0559-1
Tape C, Musé P, Simons M, et al. 2009. Multiscale estimation of GPS velocity fields. Geophysical Journal International, 179(2): 945-971. DOI:10.1111/j.1365-246X.2009.04337.x
Toda S, Kaneda H, Okada S, et al. 2016. Slip-partitioned surface ruptures for the MW7.0 16 April 2016 Kumamoto, Japan, earthquake. Earth, Planets and Space, 68: 188. DOI:10.1186/s40623-016-0560-8
Wang X W. 2017. Ionospheric error correction and three-dimensional coseismic deformation estimation and fault slip inversion with InSAR and MAI[Ph. D. thesis] (in Chinese). Chengdu: Southwest Jiaotong University.
Wang Y Z, Li Z W, Zhu J J, et al. 2012. Coseismic three-dimensional deformation of L'Aquila earthquake derived form multi-platform DInSAR data. Geomatics and Information Science of Wuhan University (in Chinese), 37(7): 859-863.
Wright T J, Parsons B E, Lu Z. 2004. Toward mapping surface deformation in three dimensions using InSAR. Geophysical Research Letters, 31(1): L01607. DOI:10.1029/2003GL018827
Yue H, Ross Z E, Liang C R, et al. 2017. The 2016 Kumamoto MW=7.0 earthquake:A significant event in a fault-volcano system. Journal of Geophysical Research:Solid Earth, 122(11): 9166-9183. DOI:10.1002/2017JB014525
Zheng A, Wang M F, Zhang W B. 2017. Source rupture process of the 2016 Kumamoto prefecture, Japan, earthquake derived from near-source strong-motion records. Chinese Journal of Geophysics (in Chinese), 60(5): 1713-1724. DOI:10.6038/cjg20170509
胡俊, 李志伟, 张磊, 等. 2013. 多孔径InSAR技术电离层校正方法及二维形变场应用研究——以玉树地震为例. 中国科学:地球科学, 43(3): 457-468.
蒋生淼, 易磊, 张旭, 等. 2018. 2016年日本熊本地震破裂时空过程联合反演. 地震学报, 40(1): 13-23.
江在森, 马宗晋, 张希, 等. 2003. GPS初步结果揭示的中国大陆水平应变场与构造变形. 地球物理学报, 46(3): 352-358. DOI:10.3321/j.issn:0001-5733.2003.03.012
江在森, 杨国华, 王敏, 等. 2006. 中国大陆地壳运动与强震关系研究. 大地测量与地球动力学, 26(3): 1-9.
金明培, 黎朕灵, 汪荣江. 2017. 日本熊本MW7.0地震同震位移场和震源滑动模型反演. 地震学报, 39(6): 819-830.
王晓文. 2017.基于InSAR和MAI的电离层误差校正及同震三维形变场计算与断层滑动反演[博士论文].成都: 西南交通大学.
王永哲, 李志伟, 朱建军, 等. 2012. 融合多平台DInSAR数据解算拉奎拉地震三维同震形变场. 武汉大学学报:信息科学版, 37(7): 859-863.
郑傲, 王铭锋, 章文波. 2017. 利用近场强震记录反演2016年日本熊本县地震震源破裂过程. 地球物理学报, 60(5): 1713-1724. DOI:10.6038/cjg20170509