据中国地震台网测定，北京时间2013年4月20日8时2分在四川省雅安市芦山县(30.3°N，103.0°E)发生 M_S7.0 级地震，震源深度13 km.这是继2008年汶川地震后龙门山地区再一次发生严重破坏性大地震，给当地人员和经济造成重大损失.地震发生后不久，各科研机构动用多种手段对此次地震进行了快速分析，初步研究了震源机制解、破裂过程等震源参数(王卫民等，2013；谢祖军等，2013;刘成利等，2013)，结果表明芦山地震是发生在龙门山断裂带南段上的一次逆冲地震，但是不同机构所给出的地震断层参数存在一定差异.地震发生后，中国地震局地质研究所与来自中国地震局地壳应力研究所和四川省地震局的科研人员一起奔赴震区进行地震应急科考，以确定地震地表破裂带的分布和同震位移分布，分析芦山7.0级地震的发震构造.但野外的科考工作并没有发现明显的地震地表破裂带，此次地震属典型的盲逆断层型地震(徐锡伟等，2013)，其发震断层尚未出露地表，隐伏在地下.对于没有地表破裂的地震，其震源和发震断层的参数不易精确确定，通过分析D-InSAR形变场和采用数值反演，则有可能定位震源，揭示隐伏发震断层(Belabbès et al., 2008).在这方面最为成功的一个震例是巴姆地震，意大利学者通过D-InSAR技术探测出该地震是由一条未知隐伏断裂诱发的，而非震源机制解给出的邻近Bam断裂(Stramondo et al., 2005).汶川地震后，许多学者利用ALOS卫星提供的L波段的观测数据获取了全面的形变场(单新建等，2009；屈春燕等，2008)，为深入理解形变过程和发震机理提供了重要的观测资料.然而该卫星于2011年4月停止运行，另外在很多震例中都成功获取过地震形变场的ENVISAT也已经停止运行，目前仅有RadarSAT-2卫星在该地区有存档数据.加拿大MDA公司调整卫星姿态后于5月10日过境拍摄了能够和存档数据组成干涉像对的影像，虽然C波段雷达系统对于像龙门山地区这样地形起伏较大、植被茂密的地区来说，保持相干性具有很大挑战，但这也是唯一可以利用InSAR手段获取芦山地震形变场的机会.本文利用该卫星数据获取了芦山地震InSAR形变场，并利用观测数据进行了同震滑动分布的反演，给出了发震断层参数，结果可为广大科研人员做进一步的深入研究提供参考. 2 芦山地震构造背景

芦山地震发生在龙门山断裂带的南段，龙门山断裂带处于青藏高原东缘与四川盆地的交汇处，是南北构造带中部一个重要的构造单元，具有十分复杂的地质构造和演化历史(张培震，2008).本次震中位于双石—大川断裂东侧，大川—双石断裂为龙门山推覆构造带南段一条区域性大断裂，野外地质调查结果表明，其走向N43°E、倾向NW、倾角45°~65°不等，其东南侧为开阔的中新生代陆相盆地，西北侧为古生代地层组成的中高山区(杨晓平等，1999).图 1为利用中国地震台网中心公布的芦山地震余震绘制的平面分布图，从图中可以看到，地震余震在龙门山前缘双石—大川出露线两侧大致呈椭圆状分布.而余震在断裂带的西北与东南两侧都有，说明双石—大川断裂是发震断裂的可能性不大，发震断层的地表行迹应该还在其东边，但在本次地震中并未出露地表.徐冲(2013)的地震滑坡分布研究成果表明了双石—大川断裂与龙门山山前褶皱之间的区域在本次地震中经历了较强烈变形，也初步印证了发震断裂不是双石—大川断裂.从中国地震局发布的芦山地震烈度分布图来看，等烈度线围绕震中大致呈椭圆状分布，这说明此次地震破裂为点状破裂，受双石—大川断裂控制不明显(中国地震局，2013).根据芦山地震主震震中位于双石—大川断裂下盘，震源深度大于10 km，盆地内震中附近不存在相应的地表大断裂等特点，推断本次芦山地震发震断裂应为盆地内的隐伏断裂.

图 1

Fig. 1

图 1 芦山地震地质构造及余震分布图 Fig. 1 Aftershock distribution of Lushan earthquake

从表 1各机构给出的震源机制解来看，本次地震震源深度在13~22 km之间，显现出浅震的特点，破裂面参数均比较接近，只是哈佛大学给出的震中位置更往东一些.截止到4月29 日零时，共记录 余震5700 余次，其中M_L2.0~2.9级667次，M_L3.0~3.9级237次，M_L4.0~4.9级42次，M_L5.0~5.9级7次(吕坚等，2013).地震余震区长轴呈东北至西南展布约40 km，破裂尺度小(杜方等，2013).

表 1(Table 1)

表 1 不同机构给出的芦山地震震源即破裂面参数 Table 1 Source parameter of the Lushan Earthquake from different agency 来源 震源位置 破裂面参数 震级Mw 最佳深度/km 经度/(°E) 纬度/(°N) 走向/(°) 倾角/(°) 滑动角/(°) Harvard 103.18 30.21 210 38 96 6.6 21.8 USGS 102.888 30.308 216 47 93 6.5 19 CENC 103.0 30.3 220 35 95 7.0 13

表 1 不同机构给出的芦山地震震源即破裂面参数 Table 1 Source parameter of the Lushan Earthquake from different agency

地表发生形变是三维的，但根据InSAR成像特点，其只能观测卫星视线向(Line of Sight)上的形变值.InSAR监测到的每个地表分辨单元视线向形 变值与地表三维形变之间的关系为(Hanson，2001)

本次地震发生时，正值ENVISAT、ALOS停止服务而欧空局、日本新一代雷达卫星尚未升空之际，这也使得数据选择的余地很小.而TerraSAR的存档数据与震后不久拍摄数据因时间相隔太长且为x波段在山区完全失相干，导致干涉处理失败.因此RadarSAT-2卫星数据也就成了唯一选择.表 1为本次干涉处理中用到的数据基本情况，两景数据为震前存档数据，一景为震后3周左右获得，可以组成两组同震干涉对.采用GAMMA软件为数据处理 平台(Werner and Wegmüller，2000)，整个D-InSAR处理采用两轨法，运用二通模式来测量地面形变量还需要有精确DEM数据.通常的处理方法是将精确的DEM数据用模拟算法生成干涉相位，再将二通模式得到的干涉相位与DEM生成的干涉相位作差分，就可以得到真实的地面形变产生的干涉相位(王超等，2002).因RadarSAT-2不再提供raw数据，因此只能从SLC数据开始处理，精密轨道数据已经包含其中.在处理过程中，地形消除采用美国宇航局发布的SRTM 3 sec数字高程模型(DEM)(Gorokhovich and Voustianiouk， 2006).处理中SAR干涉图和DEM模拟干涉图都没有去除平地效 应，在进行差分时，将地形影响和平地效应同时去除.

表 2(Table 2)

表 2 RADARSAT-2数据参数 Table 2 SAR data used in the InSAR analysis pair Inclination Master(yymmdd) Slave(yymmdd) ΔT/(days) B⊥/m h2π/cm coseismic-1 Ascending 20120702 20130510 312 -4 2504 coseismic-2 Ascending 20130110 20130510 120 196 51 注：表中参数为卫星飞行方向(升降轨)及接收日期，时间间隔ΔT，垂直基线B⊥和模糊高.

表 2 RADARSAT-2数据参数 Table 2 SAR data used in the InSAR analysis

图 2a给出了由2013年1月和5月组成的干涉对得到的缠绕的干涉图，也就是缠绕在一个相位周期(-π~π)之间的InSAR观测值，能够反映雷达数据覆盖范围内的同震形变场.

图 2

Fig. 2

图 2 芦山地震同震干涉形变场 Fig. 2 The interferometric fringes map of coseismic deformation Lushan earthquake

由于龙门山断裂带位于成都平原到青藏高原之间的交界部分，其地形起伏明显，上下两盘的地形相差较大，在上盘山区出现雷达影像的叠掩、前缩、阴影等现象，另外，该地区的植被覆盖密集，雷达回波中存在较大的噪声，使得C波段雷达卫星在这些区域失相干严重.因此在上盘大部分地区没有得到干涉图像，只是在芦山县境内有6个不完整的干涉条纹(见图 2a红色方框内)，最外侧一个条纹距离CENC所确定震中距离大约30 km，可以大致推断此次地震形变较大区域的范围.而在山前平原地带相干性较好，得到了干涉条纹，但条纹的形状显得有些复杂，且影响范围一直到100 km外的眉山市附近.为了厘清下盘条纹是否是由地震形变引发，我们用同样的方法处理了由2012年7月和2013年5月 两景图像组成的干涉像对(参数见表 2中的coseismic-2)，结果见图 2b，这对影像由于时间间隔更长，失相干现象更加严重，即使是在形变最大的芦山县境内也没有得到干涉条纹，但是在下盘平原地区却表现出和干涉像对1类似的情形.根据逻辑对比法(Massonnet and Feigl， 1995)，我们有理由排除DEM误差的影响，因为尽管这两个干涉像对的模糊高相差近50倍，但它们却表现出相似的条纹形状.只是在双流县附近(104.3°E，30.3°N)条纹方向转向平原方向，而没有向山区断层方向形成包络，有可能是受到大气等其他因素的影响，有待考证.

为了去除差分干涉图中的轨道残余相位，本文采用最小二乘方法拟合轨道残差双线性多项式模型(Wang et al., 2007)，从而削弱轨道残差的影响.多项式模型的表达式为

图 3

Fig. 3

图 3(a)经过改正后的解缠图像；(b)形变剖面图；(c)形变图和(d)四叉树重采样点 Fig. 3(a)Unwrapped map after correction; (b)Profile of displacement；(c)Displacement; (d)Distribution of sample points using quad-tree method

在成都平原内地形相对平坦，图像相干性好，相位稳定，因此在远离断层、形变可能很小的中部附近选一高相干点做为相位解缠的起始点，在确定了景的解缠起始点及缠起值后，使用枝切法进行解缠(Goldstein et al., 1988).通过枝状切线连接残差点，获取积分路径.当积分路径封闭时，积分沿着包含单个残差点的路径进行.当封闭积分路径内的正负残差相等时，封闭路径内相位也是连续变化的.高山地区由于几乎完全失相干，就算解缠后其值也是不可靠的，反而会给后面的反演工作添加干扰因素.利用该算法首先找出低于设置相干性阈值的区域和有残点的区域，对这些区域不做解缠，就只对平原地区相干性高的区域进行解缠.要提及的是，对缠绕差分相位进行多次窗口逐渐减小的自适应滤波，可大大改善相位解缠质量.在得到解缠相位后，就可以计算出视线向位移.

图 3c为解缠后得到的形变位移场，由于地面运动是三维的，而InSAR技术仅能获取一维的视线向形变，因此仅仅依靠一景干涉图无法确定地面运动的确切方向.以本次干涉所用的升轨为例，在处理中约定目标靠近卫星为正(红色)，远离卫星为负(蓝色).图中红色(正值)代表视线向缩短，朝向卫星运动；蓝色代表视线向拉长，远离卫星运动.芦山县附近在缠绕的干涉图中有6个条纹的地区，在解缠的过程中设定相干性阈值为0.9，有4个条纹的相干性低于这一阈值而未被解缠.在能够解缠的天全县附近拉一剖面AA′，(剖面位置见图 3a中AA′，剖面图见图 3b)，反映了从断层附近到远处的形变值变化趋势.近断层相对初始点的LOS位移为7 cm.假设视线向的形变都是由垂直形变造成，那么这个最大垂直形变量应该在7.7 cm左右.汶川地震后江在森研究团队在该地区架设了GPS连续观测台站，震后不久地震局地质所甘卫军研究团队就发布了基于 GPS连续观测站历史数据和地震当天全天观测数据的同震位移结果，其中芦山站(LS05)靠近断层，GPS 得到的同震垂直位移为7.6 cm，两者还是较为吻合的. 4 同震滑动分布反演 4.1 四叉树重采样

由于InSAR得到的干涉形变场数据量很大，且形变数据在空间上是高度相关的，直接用于计算既浪费计算资源也是没有必要的.为了提高反演效率，需要对形变数据进行重采样.而四叉树采样方法在形变梯度大的区域采样多，在形变梯度低的区域采样数据少，能最大限度地保留图像的空间特征.因此采取四叉树法对解缠后的数据进行下采样(Jónsson et al., 2002)，在采样过程中屏蔽掉了东部边缘部分的数据，并对未能够解缠的4个条纹进行手工数字离散化，加入到采样后的数据集中(图 3d中黑色的点)，最终得到763个观测点.采样结果显示近场数据点密集，远场稀疏，能够反映整个形变场的特征. 4.2 均匀滑动反演

以上获取地震形变场仅仅是地震研究工作的基础，关键是以此为基础，获取震源参数，解释地震的几何学及运动学特征(Wang et al., 2012).下面首先采用均匀滑动模型来获取芦山地震的断层几何参数，寻求形变量的理论值与观测值之间的最佳拟合.断层几何参数和运动参数是理解地壳同震形变的主要依据之一，具体包括断层位置(经度、纬度)、长度、深度(埋深和底界深度)、走向和倾角等；断层的运动参数包括滑移量和滑动角.由于发震断层没有明显的地表出露，断层位置(Lon，Lat)、断层埋深(Depth)、断层宽(Width)、断层长(Length)、走向(Strike)、倾角(Dip)等7个非线性参数和滑动矢量的2个分量参数同时成为反演求解的对象.由于此次地震规模较小，这里采用单断层模型来进行反演.本文采用单纯型算法(Clarke and Paradissis， 1997； Wright et al., 1999)来搜索断层参数的最优解.在反演过程中对断层的走向、倾角、滑动角、深度和断层长宽进行约束(给出上下限)以保证这些参数值具有物理意义.根据Well等给出的震级与断层长度(Length)、宽度(Width)之间的经验公式(Wells and Coppersmith， 1994)，本研究在反演中将断层长度范围限定在20~40 km，根据徐锡伟等的研究成果(徐锡伟等，2013)，将断层的上埋深(MinDepth)限定在4~12 km，下埋深(MaxDepth)限定在15~30 km，倾角θ设定在30°~50°之间.滑动角(Rake)设定在50°~130°之间，但在反演的过程中滑动角的值始终打在边界上，也就是说该参数在均匀滑动模型反演中没有给出.反演得到的其他参数值及不同研究者的结果对比见表 3.

表 3(Table 3)

表 3 不同研究给出的发震断层几何和运动学参数 Table 3 The variable source parameters for the Lushan Earthquake 来源 走向/(°) 倾向/(°) 滑动/m 经度/(°E) 纬度/(°N) 长度/km 最小深度/km 最大深度/km 震级Mw 本文 210.8±1.8 40.2±1.6 0.94±0.06 103.068±0.002 30.2377±0.004 31.2±1.5 8.2±0.3 17.5±0.4 6.52 王为民等 (2013) 205 38.5 1.59 10.2 6.7 曾祥方等 (2013) 212 44、92 12 6.7 吕坚等 (2013) 209 46/94 14 6.6 Jiang et al. (2014) 208 43 0.7 102.938 30.295 22.5 7.7 13.5 6.6

表 3 不同研究给出的发震断层几何和运动学参数 Table 3 The variable source parameters for the Lushan Earthquake

本文确定的走向与吕坚等(2013)的结果接近，滑动量小于王为民等使用地震波反演得到的结果，但 大于Jiang等利用GPS反演得到的结果(Jiang et al.，2014).反演得到的最小埋深(min_depth)为8.2 km，表明造成破裂的发震断层未出露地表，反演给出的地 震矩为7.49×10¹⁸N·m(M_w6.52)，比王卫民等(2013)根据地震波研究结果略小.

为估计断层参数的不确定性，首先基于协方差矩阵模拟100组大气噪声，并添加给观测资料，形成100组带误差的观测数据集，采用蒙特卡罗方法，利用反演程序来求取最佳折衷解(Parsons et al., 2006).图 4给出了由非线性反演获得的断层参数的不确定性分析结果，从图中可以看到，滑动量(slip)与最小埋深(min_depth)、最大埋深(max_depth)之间有极强的相关性，断层走向(strike)与最小埋深(min_depth)及滑动量(slip)之间也具有很强相关性.

图 4

Fig. 4

图 4 基于蒙特卡洛法的断层参数精度估计 Fig. 4 Uncertianty analysis of the fault parameter based on the Monte Carlo method

利用表 3反演得到的断层几何参数进行正演，拟合得到的升轨InSAR干涉图和残差图见图 5，不出意外，在影像的东部边缘，也就是彭山县以东地区的误差达到-4 cm，在分析解缠后的形变图中我们已经指出这部分地区为异常区，可以忽略.此外的大部分地区残差均不足1 cm，但是在天全县的近场地区残差也达到了3 cm，需要进一步进行分布式滑动反演.

图 5

Fig. 5

图 5 据均匀滑动模型模拟的干涉图(a)及残差图(b) Fig. 5 Model interferograms(a) and residuals(b)for calculations with slip on a single fault

当断层的位置和几何参数确定后，就可以利用非均匀滑动模型来获取断层面上的精细滑动分布.基于弹性位错理论(Okada，1985)，断层面的滑移量和观测值之间可以用如下线性关系表示:

基于均匀滑动反演结果，综合了余震精定位后的空间分布特征、震源机制解、地震地质等给出的发震构造模型，即走向取213°，倾角由39°线性渐变至43°，滑动角允许在60°~120°之间变化，将断层面扩 展到45 km×30 km，延伸到地表，并划分成1.5 km×1.5 km 的离散网格，共计600个，形成滑动分布断层模型.根据野外调查未发现明显地表破裂的结论，在断层面上接近地表深度处施加了位错为零的约束来计算芦山地震的滑动分布.最后得到残差较小的断层模型(表 4)，滑动分布结果如图 6所示.

表 4(Table 4)

表 4 模型参数 Table 4 Parameters of fault models 断层参数 震中位置 走向/(°) 倾角/(°) 长度/km 宽度/km 平均滑动角/(°) 最大滑动量/m 纬度/(°N) 经度/(°E) 深度/km 反演值 30.25 102.93 12.97 213 39~43 45 30 71 0.91

表 4 模型参数 Table 4 Parameters of fault models

图 6

Fig. 6

图 6 同震滑动分布 (a)滑动分布三维示意图;(b)断层面上的滑动分布. Fig. 6 Distributed slip for Lushan earthquake (a)3D view of the slip;(b)Slip distribution on the fault plain.

反演的位错模型相应的矩震级为M_w6.45，最 大同震位错为0.91 m，发生在震中附近区域约13 km 深度，断层滑移显示出一定的走滑分量，平均滑动角71°，整体上以逆冲滑动为主，但兼具左旋走滑特征，这与GPS的观测结果相吻合(Wu et al., 2013).矩震级与USGS给出的结果相近；矩震级M_w6.45的结果跟国内一些研究机构根据地震波反演的结果来比，中国科学院计算地球动力学重点实验室曾祥方等(2013)采用近远震波形联合反演方法得到的矩震级为M_w6.7级；中国科学院青藏高原研究所王卫民等(2013)计算得到 的地震矩为M_w6.7，释放的地震矩在1.54×10¹⁹N·m，本文结果稍小.

为了研究分布式滑动反演得到的断层精细结构与观测值之间的拟合效果，本文利用反演得到的滑动分布结果计算了拟合的干涉图及残差图(图 7a，b).从图中可以看出，相比单断层的均匀滑动模型反演，分布式滑动反演的拟合结果更为理想.从模拟干涉图像中提取AA′和BB′两条位移剖面(位置如图 7a)，并给出相应的地形起伏剖面图.从跨断层剖面AA′显示的形变特征(图 7c)可以看出，该次地震是一次典型的褶皱地震(fold earthquake)，在震中附近可见背斜状隆起(Stein and King， 1984)，说 明在震中附近30 km范围内，NW-SE方向存在着地壳缩短运动.在宝兴县内有2cm左右的视线向沉降，在芦山县境内有20cm左右的视线向隆升，形变特征就如Stain提出的倾角为45°的盲断层地表位错模型(Stein and Yeats， 1989).遗憾的是本次InSAR观测未获得上盘的干涉图像.分析残差时仅考虑InSAR形变场有值的区域，从残差图中可以发现，模型对观测值在近场区域拟合较好.残差较大的观测点主要集中在图像边缘区域，这是由于解缠后的图像本身在该区域存在误差所致.

图 7

Fig. 7

图 7(a)分布式滑动模型模拟的干涉图;(b)残差图;(c)和(d)分别对应AA′及BB′的模拟位移剖面及地形分布 Fig. 7(a)Model interferograms;(b)Residual for the distributedslip calculation;(c) and (d)show the simulated deformation and topography profiles along the AA′ and BB′ lines respectively

本次芦山地震发生在龙门山断裂带南段，震后的现场考察没有发现明显的地表破裂带.从目前发表的有关余震精定位的论文结果来看，主震及大量余震与双石—大川断裂的关系并不明显(赵博等，2013；苏金蓉等，2013)，房立华等(2013)根据芦山地震序列震源深度分布特征采用最小二乘法拟合得到的断层在深部的展布形态，推测发震断层延伸至地表的位置在余震区东侧.我们将滑动分布模型在地表的投影叠加在L and sat遥感影像上，可以看出滑动分布范围与图 1中余震分布的范围相一致，断层的地表出露线位于双石—大川断裂东侧大约12 km处.反演过程中曾多次改变断层位置进行测试性计算，无论是将图 8中推测断层位置向东还是向西偏移，都会导致更大残差.因此我们认为该推测断层是发震断层的概率较大，走向大致与双石—大川断裂平行，在双石—大川断裂以东大约12 km，这些有待野外进一步验证.

图 8

Fig. 8

图 8 芦山地震滑动分布在地表投影白网格为空间滑动分布模型在地表的投影；红线为空间滑动模型与地表的交线，黑色五星为USGS所定震中，红色五星为CENC所定震中位置. Fig. 8 Projection map of the distributed slip on a ETM+ image The white net indicated the distributed slip projected on ground surface，the red line represents the intersection of the distributed slip model with the ground sruface. Black star represent the epicenter located by USGS，red star is CENC.

本文通过InSAR形变观测，受数据源的限制，虽然只获得了芦山地震的部分形变场，但对于认识本次地震的发震机制仍具有重要意义.芦山地震发生以来，一直对发震断层存在争议，根据InSAR处理结果，结合遥感解译，确定了初步的断层参数.以 InSAR观测值为约束，基于有限弹性均匀半空间模型，对断层的滑动分布进行了反演，结果表明，所采用的断层参数能够较好地拟合观测数据，发震断层走向213°，倾角39°~43°，最大滑动位于地表以下约13 km 深度位置，最大滑动量0.91 m，平均滑动角71°，整体上仍以逆冲滑动为主.推测在双石—大川断裂以东约12 km处展布一条隐伏断裂，为本次的发震断裂.

本文在模拟时采用了均匀弹性半空间位错模型，然而分层介质的垂向变化会对形变模拟产生较大影响，下一步的工作可考虑采用介质分层模型，并进一步研究同震和震后应力应变场演化等更为重要的科学问题.

致谢 地质所陈立春博士提供了雅安地区的活动断层数据，中国地震台网中心李卫东提供了余震目录，甘卫军研究小组在第一时间公布了芦山地震GPS同震结果，德国GMZ的汪荣江老师提供了反演程序，作者与利兹大学Richard博士进行了有益的讨论，审稿人在短时间内给出很好的建议，本文部分图件使用GMT软件绘制，在此一并表示感谢！