2012, Vol. 55
2009年4 月6 日，意大利当地时间3 时32 分(北京时间9 时32 分)，意大利中部L′Aquila地区(图 1)发生了Mw6.3 级强烈地震(文中称为L′Aquila地震)，震中位于(42.334°N，13.334°E)，震源深度为8.8km(http://earthquake.usgs.gov/earthquakes/eqinthenews/2009/us2009fcaf[2010-9]).此次地震造成了300 多人死亡，1000 多人受伤，66000多人无家可归，成千上万的房屋倒塌，给当地造成了巨大的经济损失和人员伤亡[1].
L′Aquila地区位于两条大的山脉之间(图 1)，在山区由于雷达影像的叠影、前缩、阴影以及茂密的地表植被的影响，使得在得到的雷达回波中存在较大的噪声，C 波段(波长5.6cm，如Envisat卫星)和X 波段(波长3.1cm，如COSMO-SkyMed)雷达影像在该区域的干涉效果受到一定的影响，使得部分区域无法得到有意义的信号[1]；而对于L 波段的ALOSPALSAR 雷达系统而言，由于其较长的波长(23.6cm)以及较大的雷达波入射角(~34°)，可以有效地克服上述问题，从而获取到覆盖L′Aquila地震震中区域的高质量同震干涉地表形变场.本研究中所采用的Envisat卫星数据同时包含有升降轨数据(表 1)，再结合不同入射倾角的ALOS卫星得到的同震地表形变场，相当于提供了L′Aquila地震的三维地表形变场，从而有效地降低LOS(Line-Of-Sight，视线向)的模糊问题[2]，可以对地震震源机制及破裂过程反演等提供更为可靠的数据约束[3].
|
图 1 L'AqUla地区地质和地震分布图.五角星为震中位置，蓝色圆点为地震后一年内的M4〜5级(蓝色小点）和M5〜6级(蓝色大点）的余震，黑色线条为Paganica断层，其他粉色线条为区域内的活动断层，黑色三角形为行政区域分布.左上图中方框所在区域为SAR影像覆盖区，其中粉色为EnVst卫星、黑色为ALOS卫星覆盖范围 Fig. 1 Geologic setting and seismicity of the L'Aquila region. The red star is the location of epicenter and blue circles represent aftershocks (small circles: M4〜5，big circles: M5〜6). The black line is Paganica faultand pink lines are other active faults in this region. The box in insetmap means area covered by SAR images，where the pink boxes indicate the area covered by ASAR data, and the black box indicates area covered by PALSAR data |
|
表 1 研究中所选用的SAR卫星数据参数 Table 1 The SAR satellite data used in this work |
在利用InSAR 等地表形变观测资料来确定地震震源机制过程中，为了获取更为精确的断层滑动模型，需要将断层面上的非均匀滑动在整个破裂面上进行分段或区域平均，即将发震断层面离散成一系列的小矩形元或者三角元等.在这个离散过程中，通常采用的是均匀剖分[1, 3]和人为设定剖分方式(例如上小下大的形式)[3]等方法.但是按照这两种离散剖分方式反演得到的滑动分布在断层面的某些位置上会呈现出太多的细节以及给出一些基于数据本身所不能完全确定的滑动(伪滑动)特征；或者，得到的滑动分布太过于平滑以至于无法得到更为精细的分布特征.
在本研究中，首先采用二通差分干涉法对不同波长、不同入射倾角的升降轨Envisat和ALOS 卫星的ASAR 和PALSAR 数据进行干涉处理，得到了高质量、全覆盖的同震地表干涉形变场.在进行震源机制反演过程中，为了提高反演计算效率并提取精细的数据特征，采用了结合四叉树和均匀采样的方法对原始数据进行降采样，具体为对近场数据进行四叉树采样，远场数据利用均匀采样，采样数据较好地保留了形变场特征.最后结合GPS 同震地表形变观测资料，利用弹性半空间位错模型对观测数据进行震源机制联合反演，即先采用矩形位错均匀滑动模型反演确定了断层的震中位置、长度、深度、走向等几何参数和滑动量、滑动角等运动参数，然后采用断层自动剖分技术[4]对发震断层进行最优剖分，得到L′Aquila地震的最优断层面剖分模型，在此基础上采用三角位错分布式滑动反演来确定此次地震断层面上的精细的空间破裂分布，从而进一步加深对L′Aquila地震的发震机理和震源破裂过程的认识.
2 地震构造背景L′Aquila地区在上几个世纪经历了多次地震:1461年发生了一次M6.4级地震，该地震与这次的L′Aquila地震有极大的相似性；1703年又发生了一次M6.7级强震，震中位于L′Aquila地区NW 向30km，靠近图 1北部的余震区；1762、1916和1958年又分别发生震级为5.2~5.9 级地震[5].该地区地质结构复杂，活动构造运动主要由欧亚板块和非洲板块的相互作用所主导.近年来，多位学者分别利用三角测量、DInSAR 和GPS 数据对Apennines山脉的地壳扩张速率进行了估计，这几类不同的数据都得到了相似的结果，表明该区域地壳扩张速率为3~5mm/a[6-7].
L′Aquila地震后许多研究机构对此次地震进行了较为细致的研究，主要包括有野外地质调查[8]和利用各类观测数据进行的震源机制解释.地震震源机制解结果有:Anzidei 等[9]、Cheloni 等[10] 和Cirella等[11]分别利用GPS数据反演得到此次地震的主断层Paganica断层是一条SW 走向的正断层；Pino & DiLuccio[12]利用地震波数据得到了此次地震的震源机制解，表明SW 倾向断层是此次地震活动的主要原因；而Bindi等[13]、Cirella等[11]分别利用强地面位移(strongmotion)数据也给出了该地震的震源机制解；此外，Walters等[1]利用Envisat的ASAR 数据，以及Atzori 等[14] 利用Envisat 的ASAR 数据及COSMO-SkyMed卫星数据反演得到发震断层的几何参数及滑动分布，表明发震断层为SW 倾向断层.
在L′Aquila地震主震发生后的一周内，该区域又发生了7次Mw>5级和11次Mw4.0~ 4.9级的余震，其中最大的一次余震(Mw5.6)发生在4月17号，震源深度为5 km，位于主震东南15 km 处(http://www.emsc-csem.org/Earthquake/index.php[2010-10]). 主震的地震矩张量表明L′Aquila地震的发震断层是一个正断层.余震分布显示了该发震断层是一个具有NW-SE 结构和SW倾向的活动断层.
3 InSAR 同震形变场 3.1 数据处理意大利L′Aquila 地震后，欧空局(EuropeanSpaceAgency，ESA)迅速在其官网上公布了此次地震的一系列SAR 影像以供免费下载研究，提供的SAR 影像有:升降轨ASAR 影像各6 幅，其中震后影像各一幅，其余为震前影像；升轨PALSAR 震前影像3幅.本文为了更精细地研究本次地震的几何结构和运动学特征，还从日本宇宙航空研究开发机构(Japan AerospaceeXploration Agency，JAXA)申请了3幅震后ALOS卫星的PALSAR 影像.图 1中显示了本文选取的SAR 影像覆盖区域与震中位置，从中可以看出这三幅SAR 影像同时覆盖此次震中所在区域.本文所选取的SAR 数据具体参数见表 1.
在获取到卫星SAR 数据后，我们首先采用二通法[15]对ASAR 和PALSAR 影像进行差分干涉处理来获取地震的地表同震形变场.本研究中所使用的InSAR 处理平台为瑞士的Gamma 软件[16]，采用ESA 提供的DOR 精密轨道以及JAXA 提供的精密轨道来修正轨道误差以及利用90 m 分辨率的SRTM DEM[17]来去除地形影响；为进一步降低干涉相位的噪声，提高经过地形改正后的干涉图信号质量，本文采用了基于能量谱的局部自适应滤波[18]对干涉图进行滤波；最后采用枝切法(branch-cutmethod)来解缠得到差分干涉相位.完成整个干涉处理后，本文得到了三幅经过地理编码的同震形变图(图 2).
|
图 2 DInSAR干涉图(a)Envisat降轨Track 079干涉图(每一个条纹的变化代表LOS向2. 8 cm的形变）；(b)Envisat升轨Track129干涉图(每一个条纹的变 化代表LOS向2. 8 cm的形变）；(c) ALOS升轨Tack 638干涉图(每一个条纹的变化代表LOS向11. 8 cm的形变）.白线为Paganica断层，虚线为用于图 3的垂直断层剖面. Fig. 2 DInSAR interferograms (a) Interferogram of Envisat descending Track 079 (Each fringe represents 2. 8 cm LOS displacement) ； (b) Interferogram of Envisat ascending Track 129 (Each fringe represents 2. 8 cm LOS displacement) ； (c) Interferogram of ALOS ascending Track 638 (Each fringe represents 11. 8 cm LOS displacement). The white line shows the Paganica fault and the dash line which cro^s the fault trace will be used in Fig. 3. |
由于Envisat和ALOS 卫星的SAR 传感器的波长分别为5.6cm(C 波段)和23.8cm(L 波段)，因此图 2中ASAR数据的干涉条纹较PALSAR条纹更为密集，同时可以看出ALOS 的干涉图像更清晰，覆盖区域更全面，这是因为L 波段在地形复杂区域成像效果较C 波段要好，图 2中的失相关区域主要由植被覆盖茂密造成的.此外，从图 2 中可以看出L′Aquila地震地面形变的变形趋势，差分干涉条纹以Paganica断层断裂带为中心环绕分布，形变主要集中在一个长约35km，宽约30km的NW-SE向区域内，并且距离断层位置越近条纹越密集，距离断层位置越远条纹越稀疏.其中，Track079的干涉图(图 2a)中下盘有9个清晰的干涉条纹，上盘有2个清晰的干涉条纹；Track129 的干涉图(图 2b)中下盘存在8个清晰的干涉条纹；Track638的干涉图(图 2c)中下盘存在2个清晰的干涉条纹.比较断层上盘、下盘的椭圆变形区域，发现其具有明显的不对称性，表现为下盘条纹密集，上盘条纹较稀疏，这正是正断层地震造成的地表形变的一个显著特点，即下盘的地表形变要比上盘的大得多.干涉图中的不对称性还表明在东北向存在较大的形变梯度，并且升轨的不对称椭圆形变区域也表明了形变区是由一个走向约为NW-SE 向、倾向为SW 向的正断层造成的.同时，图 2 中的右边两幅升轨干涉图显示仅一个清晰的椭圆形变区，另一幅降轨干涉图显示了不对称的两个椭圆形变区，这是因为升降轨的卫星方位角的不同造成的.提取图 2中垂直断层迹线上的观测数据进行分析(图 3)，结果显示了LOS向的最大位移量为25cm，而远离断层的两侧位置，随距离增加，LOS向形变量逐渐减少，远场形变量接近于0cm，这与地震活动的物理机制相一致.
|
图 3 断层剖面(图 2中虚线位置）上的LOS形变及其拟合值(a)Envisat降轨Track 079断层剖面；(b) Envisat升轨Track129断层剖面；(c) ALOS升轨Track638断层剖面. Fig. 3 The observed and fitted LOS displacement of selected profile (dashed lines in Fig. 2)(a) Profile for Envisat descending Track 079； (b) Profile for Envisat ascending Track 129 ；c) Profile for ALOS ascending Track 638. |
在InSAR 差分干涉图中，同时包含有大气、轨道、DEM 误差、热噪声等多种误差，这些误差对干涉相位的影响方式和量级各不相同，常规的误差估计方法难以研究干涉图的误差大小.本文假定差分干涉图中的误差(主要为残余大气效应和轨道误差)的统计特征在整个影像中具有相同的空间结构[19-20]，则可以使用1D 协方差函数来描述每个Track中误差的特征(包括量级和空间尺度).研究表明，一个简单的指数函数可以很好地表征观测值的协方差函数[19-20]:
|
(1) |
其中Ci，j为相距r的像素i和j之间的协方差函数，σ2 是整个影像的方差，α 是误差衰减距离.从表 1中可以看到，Envisat卫星Track079 和Track129轨道的同震位移场的中误差均为0.68cm，方差-协方差衰减距离分别为7.3km 和3.9km.结果表明，Track079的干涉图中误差和衰减距离与Walters等[1]的类似；而Track129 的干涉图结果本文中误差较Walters等[1]的0.48cm 大，但其衰减距离更短.ALOS卫星的Track638的同震位移场的中误差为2.16cm，协方差衰减距离为10.2km.比较本文中ALOS卫星和Envisat卫星数据获取的同震位移场的中误差结果可知，前者近似为后者的3倍，这可能是由于两者波长之间的差别所造成的.由式(1)获取到的三幅不同轨道的差分干涉图的中误差将用于后续联合反演中的定权.
4 模型反演首先通过对干涉图进行降采样处理来获取适当大小的InSAR 数据集，同时为了提高反演结果的可靠性，在反演过程中我们还从文献[10]中引入了56个GPS水平形变观测值和42 个GPS 垂直形变观测数据(图 4)，这些GPS 观测数据的东西向(EW)和南北向(SN)水平形变的中误差为3.9 mm，垂直方向(V)形变的中误差为7.8mm.在此基础上，我们采用弹性半空间位错模型来进行发震断层的几何参数和滑动分布的联合反演，联合反演过程中的各类观测值的相对权比Pi根据其相应的中误差σi由式(2)确定[21]:
|
(2) |
|
图 4 GPS同震地表形变场 (a)GPS水平形变场；(b)GPS垂直形变场.黑色矢量为观测得到的同震地表形变， 红色矢量为单一断层模型均勻滑动反演结果以及蓝色矢量为分布式滑动反演结果. Fig. 4 GPS coseismic displacements of the earthquake (a) GPS horizontal deformation field；b) GPS vertical deformation field. Black，red and blue vectors are observed， simulated displacements .^rom single-fault uniform model and distributed slip model, respectively. |
文中取σ0=2.16cm，由此计算得到的相对权比为P079D∶P129A∶P638A∶PGPS_EW∶PGPS_SN∶PGPS_V=10.2∶10.2∶1∶30.7∶30.7∶7.7.在具体的反演过程中，包括以下两个步骤:一是由非线性反演方法和矩形位错模型[22]来反演均匀滑动分布模型的断层几何参数；二是在获取到断层几何参数的基础上，利用线性反演方法和三角位错模型[23]提取断层面上的精细滑动分布结构，从而确定发震断层的精细运动学特征.
4.1 混合采样图 2 中的同震形变场包含有几十万个观测数据，考虑到计算效率和反演的可行性，首先需要对三幅InSAR 干涉图进行降采样生成一个点数适当的观测数据集.在InSAR 观测数据采样中，常用的有均匀采样和四叉树采样[24].其中，均匀采样能有效降低部分误差较大观测区域结果对整体结果的影响，但是不能最大地保留图像的空间特征.而四叉树采样在高形变梯度区域采样量大，在低形变梯度区域采样数据小，能最大限度地保留图像的空间特征和采样密度.
理论上远场区域形变可假定为零或忽略不计，观测数据结果梯度非常小，而近场形变量大，观测数据结果梯度较大.已有研究表明[24]，由于受到轨道、DEM 误差、大气等误差的影响，远场观测的部分区域的噪声会大于信号，产生较大的形变梯度区域，对这些区域采样过密，会影响反演的结果，对远场数据利用均匀采样方法可避免引入过多的噪声数据；而近场区域是主要的形变区域，信号远远大于噪声的量级，利用四叉树算法能最大限度地保留观测数据特征和采样密度.
考虑到这两种采样方法的特点，为了处理上的方便，本文对远场数据进行均匀采样，对近场数据利用四叉树进行采样，采样结果如图 5.最后Track079、Track129和Track368三个轨道的原始干涉图经采样后的数据点数分别为1254、1282 和1667个，且近场点数密度大，远场点数密度小，采样结果更为合理，并且能有效地代表整个形变场的特征和量值.
|
图 5 结合四叉树与均匀采样得到的点位分布 Fig. 5 The distribution of sampled points by combination method of quad-tree and uniform sample method |
利用4.1节中采样得到的三幅干涉图的数据以及GPS同震形变资料[10]，本节首先采用均匀滑动模型来获取L′Aquila地震的断层几何参数，寻求理论形变量和实际观测值的最佳拟合.断层几何参数和运动参数是理解地壳同震形变的主要依据之一，具体包括有断层位置(经度、纬度)、长度、深度(埋深和底界深度)、走向和倾角等；断层的运动参数包括滑移量和滑动角.由于此次地震规模较小，这里采用单断层模型来进行反演.虽然InSAR 数据处理过程中使用了精密轨道信息，但是干涉图中仍然存在残余的轨道误差，因此本文在反演模型中加入了线性函数来估计轨道误差.模型反演中总共需要估计的参数有9个断层参数和9个轨道参数.反演模型所用到的断层参数初始值及范围由图 2 来确定.由于地表形变与断层几何参数之间是非线性的关系，本文中采用单纯型算法[25]来搜索断层参数的最优解.在单纯型算法中，还需要对断层的深度、长度和宽度进行约束(给出上下限)以保证这些参数值具有物理意义，此外未对其他参数进行约束.
反演过程中，需要解如下方程:
|
(3) |
其中d为InSAR、GPS 同震形变观测值，G为设计矩阵，m为包括断层位置、深度、长度、滑动量和轨道参数等的参数集，ε 为观测误差.反演的最终目的是为了使观测数据和模型最佳拟合，从而达到目标函数最小.
本文采用okinv程序[25, 26]来进行反演，均匀滑动模型反演得到的参数值见表 2，表明发震断层是一个以正倾滑为主兼有少量右旋走滑的断层，其走向近似为NW-SE方向，倾角为51.5±0.4°，由表 2可知，本文震中位置与GPS 确定的结果[9]相一致.断层走向同Walters等的结果[1]相一致，但大于其它研究机构的结果，这是由于相比GPS 观测数据，InSAR 数据在走向角的约束上具有更强的能力[1].断层的埋深为2.9km，意味着造成该破裂的发震断层并没有出露于地表.此外，本文确定的断层平均滑动量为0.607±0.005 m，滑动量同Cheloni等[9]利用GPS 数据得到的0.62 m 相一致，略小于Walters等[1]利用InSAR 数据确定的0.66m；反演得到的滑动角为-106.3±0.5°，与Walters等[1]和Atzori等[14]单独利用InSAR 反演得到的-105°和-103°相一致，而比USGS研究机构根据地震波数据发布的-112°略小.反演给出的地震矩为2.90×1018 N·m(Mw6.24)，与基于InSAR 数据[1, 14]以及体波资料得到的研究结果相当.
|
|
表 2 不同研究给出的反震断层几何和运动学参数 Table 2 The various source parameters for the L Aquila earthquake |
一般来说，在非线性反演中还没有很好的方法可以将观测值的误差传递给模型参数，本文使用蒙特卡洛方法，通过利用100 组合成的带误差观测数据来估计断层参数的不确定性.图 6给出了由非线性反演获得的断层参数的不确定性，表明本文确定的断层参数具有较高的可信度.图 6中的散点还显示了断层参数之间的最优权衡(trade-off)，从图中可以看到，最大深度(Maxdepth)和最小深度(Mindepth)与滑动量(Slip)间有极强的相关性；断层的经度(Lon)与断层的纬度(Lat)、断层走向(Strike)、倾角(Dip)和滑动角(Rake)间亦存在着相关性；此外断层的纬度(Lat)与断层走向(Strike)和滑动角(Rake)间有很强的相关性.
|
图 6 非线性反演的不确定性分析，基于蒙特卡洛方法的断层参数精度估计 Fig. 6 Uncertainty analysis for the nonlinear inversion and parameter precision estimation based on the Monte Carlo method |
利用表 2反演得到的断层几何参数进行正演，得到的GPS同震地表形变场见图 4、拟合得到的干涉图和残差图见图 7.从图 4、7可以看出，单断层的均匀滑动模型反演获取的结果总体上与观测值拟合良好，其中GPS东西向、南北向和垂直向形变的均方中误差分别为0.7cm、1.1cm 和0.6cm；三幅残差图(图 7)的均方中误差分别为1.1cm、1.2cm 和2.0cm，并且残差图中靠近断层线两侧部分出现了约1.5个条纹的误差，这主要是由于形变的中心位置形变梯度较大，破裂情况复杂，表明单一断层均匀模型在同震形变场的细节描述上的能力不足，需要进行进一步的分布式滑动反演.
|
图 7 均匀滑动模型拟合的干涉图和残差图 (a) Track 079拟合干涉图；(b) Track 129拟合干涉图；(c) Track 638拟合干涉图；(d) Track 079残差干涉图；(e) Track 129残差干涉图；(D Track 638残差干涉图.拟合干涉图和残差图中Track 079和Track 129的每一个条纹的变化代表LOS向2. 8 cm的形变，Track 638的每一个条纹的变化代表LOS向11. 8 cm的形变. Fig. 7 The fitted and residual interferograms based on the uniform slip model (a) The fitted interferogram Dor Envisat Track 079? (b) The fitted interferogram for Track 129 ；c) The fitted interferogram for ALOS Track 638? (d) The residual interferogram for Envisat Track 079； (e) The residual interferogram for Track 129? (() The residual interferogram for ALOS Track 638. Each fringe represents 2. 8 cm and 11. 8 cm for the Envisat and ALOS data respectively. |
在获取到断层的几何参数后，就可以利用非均匀滑动模型来获取断层面上的精细滑动分布.在对L′Aquila地震的研究上，已有研究是将断层面均匀地离散成相同大小(1km×1km)的断层片[1, 8, 12].在这种情况下，反演得到的滑动分布通常都不是观测数据所对应的最优结果.这表现在反演得到的滑动分布在断层面的某些位置上呈现出太多的细节以及给出了一些数据所不能确定的滑动(伪滑动)特征；或者，得到的滑动分布太过于平滑以至于无法得到较为细致的分布特征.在大多数情况下，由于发震断层的位置以及观测数据的分布，反演给出的结果同时具有以上特征.因此为了给出L′Aquila地震的最优断层剖分以及分布式滑动分布模型，即每个断层片上滑动都是可以通过足够的观测数据所确定的，这里采用了断层自动剖分方法[4]来建立滑动分布模型.
在分布式滑动模型反演中，根据设计矩阵(格林函数)G可以计算得到模型的分辨率矩阵R:
|
(4) |
对于同震滑动分布而言，分辨率矩阵R中的每一行表示净滑移在相关断层片中的平滑程度[27].即如果断层模型合理，R中每行的滑动将是一个中心位于对应断层片的δ 函数.反之，该滑动量将分布到其他相邻的断层片上.基于模型分辨率矩阵R，可以确定反演得到的断层模型是否适定.在此基础上，利用高斯曲线拟合确定模型分辨率矩阵与断层元最佳尺度之间的关系，即确定需要继续进行剖分的断层元进行更进一步的剖分，从而生成一个新的断层模型来更好地拟合观测数据.
在本研究中，首先将断层的长度和宽度分别拓展为30km 和20km，初始时采用三角元将断层剖分为16个断层片，然后根据式(4)计算初始断层模型的分辨率矩阵，以确定断层元的最佳尺度大小，对尺寸大于相应尺度的断层片采用四叉树采样法[21]进行进一步的剖分后重新计算其矩阵分辨率矩阵.重复上述过程，直到各断层片的尺寸都满足各自最佳断层尺度的要求，最终得到的滑动分布模型见图 8.
图 8中的基于断层自动剖分技术给出的最优断层模型由322个不同尺寸的三角元组成，并且可以看到，这些三角元中，尺寸较小的三角元主要分布在断层的近地表位置，而尺寸较大的三角元分布在深部位置；同震滑动分布主要发生在5~14km 深度范围内，平均滑移量为0.34m，平均滑动角为-96.0°，最大滑移量达1.07 m，滑动角为-102.8°，深度为6.1km.分布式滑动分布结果显示了断层为正倾滑和少量的右旋走滑，与野外调查结果一致[13].为了估计分布式滑动反演获取的断层面滑动分布的精度，与均匀滑动分布模型的精度确定方法类似，这里也对原始观测数据生成100组带随机扰动误差的数据集，通过这些数据计算出相应的滑动分布结果，从而估计出模型的精度，得到的同震滑动分布的误差如图 8b.从图 8b中的误差分布图中可以看出，断层滑移的误差分布比较均匀，平均误差为0.01 m，说明分布式滑动模型反演得到的滑动分布是可靠的，其最大误差分布在断层面的右下方位置，达0.03m.
|
图 8 同震滑动分布(a)及其误差图(b) Fig. 8 Coseismal slip distribution (a) and its error (b) |
为了研究分布式滑动反演得到的断层精细结构与观测值之间的拟合效果，本文利用反演得到的滑动分布结果计算了拟合的GPS同震地表形变场(图 4)和干涉图及残差图(图 9).从图 4中可以看出，相比单断层的均匀滑动模型反演，分布式滑动反演的拟合结果更为理想，其GPS 东西向、南北向和垂直向形变的均方中误差分别为0.6cm、0.7cm 和0.6cm，其GPS东西方向和南北向形变的均方中误差分别减小了0.1cm 和0.4cm.从图 9中可以看到，分布式滑动反演结果与观测值之间拟合得非常好，模型生成的干涉条纹非常清晰，同时残差结果也很小，Track079、Track129 和Track638 的数据拟合模型残差中误差分别为0.8cm、1.0cm 和2.0cm，其中Track079和Track129的残差均方中误差分别减小了0.3cm 和0.2cm.分布式滑动分布模型的数据拟合率达到了99.3%.从图 9的残差干涉图中可以看到，残差分布显示出很强的随机性，这主要与大气延迟、电离层异常、DEM 误差等因素相关.
|
图 9 分布式滑动模型拟合的干涉图和残差图 (a) Track 079拟合干涉图；(b) Track129拟合干涉图；(c) Track 638拟合干涉图；(d) Track 079残差干涉图；(e) Track 129残差干涉 图；(D Track 638残差干涉图.拟合干涉图和残差图中Track 079和Track129的每一个条纹的变化代表LOS向2. 8 cm的形变，Tack 638的每一个条纹的变化代表LOS向11. 8 cm的形变. Fig. 9 The fitted and residual interferograms base on the distributed slip model (a) The fitted interferogram for Envisat Track 079； (b) The fitted interferogram for Track 129 ；c) The fitted interferogram for ALOS Track 638； (d) The residual interferogram for Envisat Track 079； (e) The residual interferogram for Track 129； (f) The residual mterferogram for ALOS Track 638. Each fringe represents 2. 8 cm and 11. 8 cm for the Envisat and ALOS data respectively. |
此外，本文提取了干涉图中垂直于断层线的一条剖面线上的观测结果和两种不同模型的拟合结果来进行比较(图 3)，结果显示出两个反演模型都能较好得解释同震InSAR 形变场，且拟合残差都较小，在0值附近，由于分布式滑动模型顾及了断层面上滑动的不均匀性，对观测资料的拟合更好.图 4中三个轨道的拟合结果中ASAR 数据的拟合效果较PALSAR 数据的拟合效果更好，符合上面数据处理部分获取的PALSAR 干涉结果误差大于ASAR干涉结果的结论.在这三个轨道的拟合结果中，远场结果拟合情况均要好于近场结果，这主要是由于近场破裂情况更为复杂所造成的.
5 讨 论干涉图和地震学结果均显示2009 年L′Aquila地震破裂是沿着SW 倾向正断层，该断层位于L′Aquila的NE 向.Cirella等[11]利用强位移和GPS数据联合反演此次地震的破裂过程，结果显示破裂传播是沿着SE 向区域、破裂速度在断层走向上为2km/s，破裂持续时间~6.5s，这个结果与本文反演得到的断层长度13.36km 相一致.
利用断层自动剖分技术获取了地震的最优滑动几何模型由322个三角断层片组成，断层片数目小于按照文献[1, 10, 14]中断层面离散方式(1km×1km)给出的600个，并且在这些三角元中，可以分辨的最短边长为0.4km，优于1km×1km，最长边长为6.3km.此外，由降采样原理可知，相对于地表形变变化平缓区域，高形变梯度区具有更为密集的采样点，即与远场区域相比，震中位置具有更多的采样点，这意味着经过四叉树采样后的观测数据能够对浅层滑动提供更多的约束，这在最优滑动模型(图 8)中得到了证实.从图 8 中可以看到，与均匀剖分技术不同，基于断层自动剖分技术给出的滑动分布模型的断层上部区域具有更多更小尺度的三角元，这表明这些位置上的滑动分布具有更好地分辨率，即这些区域得到了更多的数据约束.
Cheloni等[10]利用布设在L′Aquila震中附近的36个永久GPS站观测数据研究此次地震的同震滑动分布和初始的余震滑动分布，结果显示最大滑动量为1.1m，同震的地震矩为3.9×1018 N·m，最大滑动量结果与本文的1.07 m 非常接近，但是地震矩较本文的分布式滑动结果3.43×1018 N·m 大，主要原因可能是两者数据源的差别所造成的.
Walters等[1]联合两轨ASAR 数据反演得到断层的滑动量为0.6m，地震矩为2.80×1018 N·m，而Atzori等[14]联合两轨ASAR数据和一轨COSMO-SkyMed升轨X 波段SAR 数据研究得到断层滑动值为40~60cm，主要位于在断层从西北边缘到南部6km 的深度，在7km 深度，滑动量达到最大值0.9m，地震矩为2.7×1018 N·m，Cirella等[11]利用强位移数据反演得到的最大滑移量为1m，这三者与本文反演得到的最大滑动深度6.1km、最大滑移量1.07m 相一致，本文结果更好地支持了强震动位移反演结果，但是比较地震矩，本文地震矩为3.43×1018 N·m，大于Walters等[1]和Atzori等[14]研究的地震矩结果，与基于GPS 数据反演给出的结果[8, 9]更为接近.这可能是由于本文利用的PALSAR数据跨度时间更长，因此反演得到的地震矩结果中包含有更多的震后分量；此外，采用的GPS 数据也对干涉数据产生了约束作用.
地质野外考查结果发现该地震在地表部分区域发生了7~10cm 的破裂[8]，而根据已有的GPS[9-11]、InSAR[1, 14]反演得到的滑动分布结果均表明，断层主位错没有达到地表.本研究得到的最佳拟合分布式滑动结果与野外考察结果[8]更为一致，表明Paganica断层的最浅深度靠近地表，仅在部分区域出现了一些较大的地表滑动(图 8 中的断层上端的蓝色位置)，这个地表破裂可能是由于地下的滑动能量传递到地表的脆性区域造成的.
本文确定的震中位置与其他几家机构确定的震中位置之间存在一些差别.根据本文和已有SAR研究结果[1, 14]显示，此次震中位置没有位于高速滑动区域.地震活动是由于应力的集中释放，但由于地质构造的复杂性以及地球内部脆性和刚性的非均匀分布，造成的破裂及滑动程度的不一致，即震中与高速滑动区域间并没有必然联系，详细的解释地震活动应力释放与高速滑动区域间的关系，依赖于对地下断层的地质构造、脆弱程度以及应力释放的传播方式的了解.
许多学者认为Paganica断层与相邻断层比较在地形活动上更弱[28]，但来自当地的和遥感观测结果揭露了此次地震发生在Paganica断层，该断层恰好位于已认识的两个主要活动断层之间[10, 28-30]，表明此次构造加载非单一的主地震带.这次观测使得在L′Aquila以及相似地质区域确定潜在危险地震断层这一问题更加突出，应该提高多种测量手段对震间应变积累进行研究.
6 结 论通过联合Envisat数据和ALOS 卫星SAR 数据来反演确定L′Aquila地震的震源机制解，有以下几点结论:
(1) 本文在顾及形变场特征的情况下采用了结合四叉树与均匀采样的混合采样方法，联合ASAR数据、ALOS 数据以及GPS 数据来约束L′Aquila地震的震源机制，丰富了观测数据类型，并在此基础上利用基于弹性半空间矩形位错模型来反演L′Aquila地震的断层参数和滑动量，较好地拟合了观测到的GPS地表形变场和同震干涉形变场.
(2) 作为此次L′Aquila 地震的发震断层Paganica断层，是一条沿着SW 倾向的正断层，断层面上的最大滑动量为1.07 m，整个断层位于地下.本文均匀滑动反演和分布式滑动反演给出的地震矩结果分别为2.90×1018 N·m 和3.43×1018N·m，研究结果与地震学等的研究成果相当一致.
(3) 断层面剖分方式对最终的滑动分布特征有较大的影响，本文基于模型分辨率矩阵，采用断层自动剖分技术给出了L′Aquila地震发震断层的最优断层模型，该模型由322个三角元组成，其中的最短边长为0.4km，最长边长为6.3km，该模型在部分地表区域有10~20cm 的破裂.
(4) 弱活动性断层在地震危险性评估中同样值得重视，在L′Aquila及相似区域存在潜在的地质活动灾害的可能性有赖于增加各种测量手段对断层间的应变积累进行更为深入的研究.
致谢感谢两位审稿专家对本文提出的宝贵意见.感谢欧空局(ESA)提供的Envisat卫星数据(Cat-1No.7195)和日本宇宙航空研究开发机构(JAXA)提供的ALOS卫星数据(RA3No.608).
| [1] | Walters R J, Elliott J R, D'Agostino N, et al. The 2009 L'Aquila Earthquake (Central Italy): a source mechanism and implications for seismic hazard. Geophys. Res. Lett. , 2009, 36: L17312. DOI:10.1029/2009GL039337. |
| [2] | Fialko Y, Sandwell D, Simons M, et al. Three-dimensional deformation caused by the Bam, Iran, earthquake and the origin of shallow slip deficit. Nature , 2005, 435(7040): 295-299. DOI:10.1038/nature03425 |
| [3] | Zhang G H, Qu C Y, Shan X J, et al. Slip distribution of the 2008 Wenchuan Ms 7. 9 earthquake by joint inversion from GPS and InSAR measurements: a resolution test study. Geophys. J. Int. , 2011, 186(1): 207-220. |
| [4] | Barnhart W D, Lohman R B. Automated fault model discretization for inversions for coseismic slip distributions. J. Geophys. Res. , 2010, 115: B10419. DOI:10.1029/2010JB007545. |
| [5] | Stucchi M, Camassi R, Rovida A, et al. DBMI04, il database delle osservazioni macrosismiche dei terremoti italiani utilizzate per la compilazione del catalogo parametrico CPTI04 (in Italian), Quaderni Geofis. Ist. Naz. di Geofis. e Vulcanol., Rome, 2007. |
| [6] | Hunstad I, Selvaggi G, D'Agostino N, et al. Geodetic strain in peninsular Italy between 1875 and 2001. Geophys. Res. Lett. , 2003, 30(4). DOI:10.1029/2002GL016447. |
| [7] | D'Agostino N, Avallone A, Cheloni D, et al. Active tectonics of the Adriatic region from GPS and earthquake slip vectors. J. Geophys. Res. , 2008, 113(B12413). DOI:10.1029/2008JB005860. |
| [8] | Emergeo Working Group. Rilievi geologici di terreno effettuati nell'area epicentrale della sequenza sismica dell'Aquilano del 6 aprile 2009. Ist. Naz. di Geofis. e Vulcanol., Rome. 2009 (Available at http: //www.earth-prints.org/handle/2122/5036). |
| [9] | Anzidei M, Boschi E, Cannelli V, et al. Coseismic deformation of the destructive April 6, 2009 L'Aquila earthquake (central Italy) from GPS data. Geophys. Res. Lett. , 2009, 36: L17307. DOI:10.1029/2009GL039145. |
| [10] | Cheloni D, D'Agostino N, D'Anastasio E, et al. Coseismic and initial post-seismic slip of the 2009 Mw6. 3 L'Aquila earthquake, Italy, from GPS measurements. Geophys. J. Int. , 2010, 181(3): 1539-1546. |
| [11] | Cirella A, Piatanesi A, Cocco M, et al. Rupture history of the 2009 L'Aquila (Italy) earthquake from non-linear joint inversion of strong motion and GPS data. Geophys. Res. Lett. , 2009, 36: L19304. DOI:10.1029/2009GL039795. |
| [12] | Pino N A, Di Luccio F. Source complexity of the 6 April 2009 L'Aquila (central Italy) earthquake and its strongest aftershock revealed by elementary seismological analysis. Geophys. Res. Lett. , 2009, 36: L23305. DOI:10.1029/2009GL041331. |
| [13] | Bindi D, Pacor F, Luzi L, et al. . 3, 2009 L'Aquila earthquake: source, path and site effects from spectral analysis of strong motion data. Geophys. J. Int. , 2009, 179(3): 1573-1579. DOI:10.1111/gji.2009.179.issue-3 |
| [14] | Atzori S, Hunstad I, Chini M, et al. Finite fault inversion of DInSAR coseismic displacement of the 2009 L'Aquila earthquake (central Italy). Geophys. Res. Lett. , 2009, 36: L15305. DOI:10.1029/2009GL039293. |
| [15] | Massonnet D, Rossi M, Carmona C, et al. The displacement field of the Landers earthquake mapped by radar interferometry. Nature , 1993, 364(6433): 138-142. DOI:10.1038/364138a0 |
| [16] | Werner C L, Wegmüller U, Strozzi T, et al. GAMMA SAR and interferometry processing software. ERS ENVISAT Symp. Gothenbury, Sweden , 2000: 16-20. |
| [17] | Farr T G, Rosen P A, Caro E. The Shuttle Radar Topography Mission. Rev. Geophys. , 2007, 45. DOI:10.1029/2005RG000183. |
| [18] | Goldstein R M, Werner C L. Radar interferogram filtering for geophysical applications. Geophys. Res. Lett. , 1998, 25(21): 4035-4038. DOI:10.1029/1998GL900033 |
| [19] | Hanssen R F. Radar Interferometric: Data Interpreatation and Error Analysis. New York: Kluwer Academic Publishers, 2002 . |
| [20] | Parsons B, Wright T, Rowe P, et al. The 1994 Sefidabeh (eastern Iran) earthquakes revisited: new evidence from satellite radar interferometry and carbonate dating about the growth of an active fold above a blind thrust fault. Geophys. J. Int. , 2006, 164(1): 202-217. DOI:10.1111/gji.2006.164.issue-1 |
| [21] | Jónsson S, Zebker H, Segall P, et al. Fault slip distribution of the 1999 Mw7. 1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements. Bull. Seism. Soc. Am. , 2002, 92(4): 1377-1389. |
| [22] | Okada Y. Surface deformation due to shear and tensile faults in a half-space. Bull. Seism. Soc. Am. , 1985, 75(4): 1135-1154. |
| [23] | Meade B J. Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space. Comput. Geosci. , 2007, 33(8): 1064-1075. DOI:10.1016/j.cageo.2006.12.003 |
| [24] | Lohman R B, Simons M. Some thoughts on the use of InSAR data to constrain models of surface deformation: noise structure and data downsampling. Geochem. Geophys. Geosyst. , 2005, 6: Q01007. DOI:10.1029/2004GC000841. |
| [25] | Clarke P J, Paradissis D, Briole P, et al. Geodetic investigation of the 13 May 1995 Kozani-Grevena (Greece) earthquake. Geophys. Res. Lett. , 1997, 24(6): 707-710. DOI:10.1029/97GL00430 |
| [26] | Wright T J, Parsons B E, Jackson J A, et al. Source parameters of the 1 October 1995 Dinar (Turkey) earthquake from SAR interferometry and seismic bodywave modelling. Earth Planet. Sci. Lett. , 1999, 172(1-2): 23-37. DOI:10.1016/S0012-821X(99)00186-7 |
| [27] | Du Y, Aydin A, Segall P. Comparison of various inversion techniques as applied to the determination of a geophysical deformation model for the 1983 Borah Peak earthquake. Bull. Seism. Soc. Am. , 1992, 82(4): 1840-1866. |
| [28] | Boncio P, Lavecchia G, Pace B. Defining a model of 3D seismogenic sources for seismic hazard assessment applications: the case of central Apennines (Italy). J. Seismol. , 2004, 8(3): 407-425. DOI:10.1023/B:JOSE.0000038449.78801.05 |
| [29] | D'Agostino N, Giuliani R, Mattone M, et al. Active crustal extension in the Central Apennines (Italy) inferred from GPS measurements in the interval 1994-1999. Geophys. Res. Lett. , 2001, 28(10): 2121-2124. DOI:10.1029/2000GL012462 |
| [30] | Galli P, Galadini F, Pantosti D. Twenty years of paleoseismology in Italy. Earth Sci. Rev. , 2008, 88(1-2): 89-117. DOI:10.1016/j.earscirev.2008.01.001 |