地球物理学报  2015, Vol. 58 Issue (5): 1634-1644   PDF    
基于Landsat8光学影像的巴基斯坦Awaran Mw7.7地震形变监测及参数反演研究
冯光财1, 许兵1, 单新建2, 李志伟1, 张国宏2    
1. 中南大学地球科学与信息物理学院雷达遥感研究室, 长沙 410083;
2. 地震动力学国家重点实验室, 中国地震局地质研究所, 北京 100029
摘要:2013年9月24日巴基斯坦俾路支省(Balochistan)境内的阿瓦兰县(Awaran)发生了Mw7.7级地震.本文利用覆盖该地区的Landsat 8数据,基于影像配准的方法获取了该次地震的同震形变场,并运用地统计的方法对形变结果进行精度评定.针对传统四叉树算法中近场和远场中采样密度的不均匀性,以及噪音区域对数据降采样和反演结果收敛性的影响,本文提出了改进的四叉树算法对点的密度和形变梯度进行合理兼顾.最后利用光学影像获取的形变结果和数据的精度水平,基于Okada弹性半空间形变模型反演了该地震的震源参数和断层滑动分布.结果表明,地震断层北倾47°,滑动以左旋走滑为主,断层的西南部兼具少量的倾滑运动分量,断层滑动主要集中分布在断层面0~15 km深度范围,最大滑动量达10 m.反演获得的地震标量矩为4.68×1020N·m,震级约为Mw7.75级.本文的研究结果可以为该地区的地壳应力变化研究和地震灾害评估提供依据,同时为Landsat 8光学影像应用于地震的形变研究提供参考.
关键词巴基斯坦Awaran地震     光学影像     地震形变     断层滑移量    
Coseismic deformation and source parameters of the 24 September 2013 Awaran, Pakistan Mw7.7 Earthquake derived from optical Landsat 8 satellite images
FENG Guang-Cai1, XU Bing1, SHAN Xin-Jian2, LI Zhi-Wei1, ZHANG Guo-Hong2    
1. Laboratory of Radar Remote Sensing, School of Geoscience and Info-Physics, Central South University, Changsha 410083, China;
2. State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China
Abstract: On September 24th 2013, an earthquake of Mw7.7 occurred North-Northeast (NNE) of Awaran, Balochistan Province, Pakistan, causing extensive damages to cities and counties along the Hoshab fault. At least 825 people were killed and hundreds more were injured. However, there were quite limited SAR and geodetic data available, because the ENVISAT/ASAR satellite was out of service and ALOS/PALSAR stopped working at that time. Fortunately, this event was completely mapped by the Apr. 2013 launched Landsat8 optical satellites, which is usually for the change detection of the earth. So this earthquake provides a unique chance to study how to map the coseismic deformation of great earthquake with the Landsat 8 optical images.
In this study, we utilize the Landsat 8 images to map the coseismic deformation based on the cross-spectrum correlation method, which is usually used in SAR pixel offset tracking measurement. A linear model is used to remove the long wavelength phase ramp caused by the orbital error of satellite. Our results indicate there is a curved oblique strike-slip faulting along the pervious Hoshab fault and the maximum horizontal slip is up to 12 m in the center of fault trace.In order to weigh the different datasets in modeling, we propose a geostatistic approach to assess accuracy of the North-South (NS) and East-West (EW) offsets derived from the Landsat 8 optical images. The results show their accuracy is about 0.6 m and 0.5 m respectively. In order to reduce the effects from topography change, atmosphere and low coherence area, we mask the measurements on the water and cloud regions. We propose a modified quadtree resampling method that not only takes account of density of samples but also gradient of deformation to overcome the limitation of conventional quadtree resampling method near and far earthquake rupture area. We keep both the largest and smallest windows as 128×128 and 8×8 pixels, respectively, to ensure sufficient point density in the non-deforming area and proper samples near the surface rupture. We locate the curved fault traces and construct the fault geometry based on the discontinuity in surface offsets. We use a two-step approach, combining the simulated annealing (SA) optimization and subsequent derivative-based Levenberg-Marquardt inversion, to solve for the fault parameters in nonlinear modeling. It is a simplification to assume uniform slip on fault planes since heterogeneity in fault displacement and rake is observed on all scales of faults. Therefore, we construct a complex curved fault plane with dip angle estimated above to match the fault trace we mapped and divide it into 34×5 subfaults (7 km×7 km) to estimate variations in slip on the fault.
Based on Okada half-space elastic model,the slip distribution of this earthquake is inferred from the Landsat 8 measurements. The inversion results show that the left-lateral strike-slip dominates the earthquake fault slip concentrating at the depth less than 15 km, with a slight dip-slip component and the maximum slip is up to 10 m. The inverted seismic moment is 4.68×1020 Nm, corresponding to Mw7.75. We also find that some thrust faults in the Kirthar range and strike-slip faults within the Chaman fault system have been brought closer to failure by this earthquake.This research can help the earthquake disaster evaluation in this area, and also shed some lights on the earthquake study with the Landsat 8 optical images in future.
Key words: Pakistan earthquake     Optical satellite images     Earthquake deformation     Fault slip    
 1 引言

2013年巴基斯坦Awaran地震发生于当地时间9月24日16时29分48秒,震中位于巴基斯坦俾路支省(Balochistan)境内的阿瓦兰县(Awaran)北北东方向66 km处,矩震级Mw7.7,截止到震后3个月已造成大约825人丧生,700多人受伤(USGS,2013a; 科学网, 2014Wikipedia,2013a).主震4天后(2013年9月28日),另一个震级Mw6.8的地震,又至少造成45人丧生(USGS,2013b).这次地震发生后,世界上很多著名的地学科研机构都进行了关注.USGS在地震发生后迅速用地震波资料给出了该次地震的震源机制,显示震中位于(26.951°N,65.501°E),震源深度约为15 km,地震的主破裂面走向216.5°,北倾44.7°,是一个左旋走滑为主兼少许逆冲的地震(USGS,2013a).Global CMT 也给出了类似的机制( Global CMT,2013),但是与USGS给出的地震中心位置有少许差别(26.64°N,65.04°E)(如表 1所示).此外,中国和欧洲的地震机构对这次地震也有相关报道.

表 1 不同地震机构估计的2013年巴基斯坦Awaran地震震源参数表 Table 1 Earthquake source parameters of the 2013 Pakistan Awaran earthquake from different research institutions

2013年Mw7.7巴基斯坦Awaran地震区域的板块运动非常活跃.一方面,坚硬的印度板块(India plate)以约40 mm/a的速度向北移动碰撞亚欧大陆(Eurasia plate),通过地壳物质的挤压和紧缩在印度板块的北缘和西边缘生成了一系列的山脉,如喜马拉雅(Himalayan),昆仑山(Karakoram),帕米尔(Pamir)和兴都库什山(Hindu Kush)山脉,也引发了许多的Mw7级以上强震(Wikipedia,2013a)(如图 1所示).另一方面,阿拉伯板块(Arabia plate)也以约14 mm/a的速度往北-东北方向滑向亚欧板块,在莫克兰俯冲带(Makran trench)发生冲撞,阿拉伯板块在该处俯冲至亚欧板块之下,两个板块边界之间碰撞会引发地震和海啸.历史上有记录的如1945年11月27日的8级逆冲型地震,引起巨大的海啸,袭击了阿曼海湾和阿拉伯海域,导致4000多人死亡(USGS,2013a). 2013年9月24日的地震就发生在印度板块向北挤压和阿拉伯板块向下俯冲欧亚板块的转换断裂Chaman fault的分支断裂Hoshab fault上.Chaman fault南起阿拉伯板块、亚欧板块和印度板块在巴基斯坦莫克兰的交汇处,经巴基斯坦俾路支省(Balochistan),阿富汗的东北部喀布尔(Kabul)交汇于帕米尔断层系统(Pamir),全长约900 km,是世界上几大著名的转换断裂之一.初步的地震波资料显示该地震机制是一个左旋走滑为主兼有少量倾滑的地震,发震机制与Chaman fault也非常吻合(Wikipedia,2013b; Shuhab et al., 2013).

图 1 2013巴基斯坦Awaran地震地质构造背景图.红色为板块运动碰撞的边界线,印度板块向北东方向以43 mm/a的速度挤压亚欧板块,阿拉伯板块以14 mm/a的速度俯冲于亚欧板块之下,因此沿着板块碰撞边界发生了数以万计的大地震.右图中的地震机制图是自1976年以来Global CMT记录的里氏Mw5级以上的地震,紫色的五角星显示的 是自1900年以来USGS记录的大于Mw7级的大地震Fig. 1 Tectonic setting of 2013 Pakistan earthquake. Red lines are the boundary of plate collision in left subplot. Due to Indian plate underthrust beneath the Eurasian plate at a velocity of 43 mm/a towards northeast, and Arabian plate dives to Eurasian plate with a speed of 14 mm/a,tens of thous and s of earthquakes have occurred at the plate collision boundary. Right subplot shows a figure of earthquake mechanisms recording earthquakes of magnitude great than Mw5.0 by Global CMT since 1976. Purple stars are earthquakes of magnitude great than Mw5.0 recorded by USGS since 1900

2013年Mw7.7巴基斯坦Awaran地震发生在人烟稀少的山区,当地的GPS连续观测站以及传统大地测量观测都非常稀疏,因此该地震的研究主要依靠卫星遥感的手段进行.目前正处在新老SAR卫星交接的时间窗口,一直以来用于地震形变研究的ENVISAT和ALOS-1卫星都已经相继不能工作,新的卫星要到2014年才能发 射.而在轨的只有高分辨率的SAR卫星,如TerraSAR和COSMO-SkyMed卫星,但是它们都主要关注城市区域,并且X波段对于地震的大形变非常敏感,失相干也非常严重,数据覆盖范围也非常有限.该地震区域处于干燥的山区,很少有云遮挡,气象条件非常 好,适合光学影像成像,因此2013年发射的L and sat 8 光学卫星影像就成为了该地震重要的观测数据(USGS,2013c).关于如何用光学影像监测地震形变,美国加州理工学院在这方面已经做了不少工作,他们基于ENVI软件利用IDL语言编写了COSI- Corr程序包,并已经成功用于Aster,Spot,Quickbird 等高分辨率光学影像的处理,获取水平方向的同震形变(Leprince et al., 2007).除此之外,他们还广泛地把这种方法应用于冰川和滑坡的形变研究当中(Caltech,2007).然而ENVI和IDL语言包都是商业软件并且价格昂贵,因此该软件的应用和推广受到一定的限制和约束.

为了研究该地震的震源参数,本文首先利用L and sat 8光学卫星影像,通过ROI_PAC软件的偏移量配准算法获取2013年9月24号巴基斯坦Awaran地震形变.然后通过形变信息的信噪特征,基于地统计的变异函数方法对形变结果进行精度评定.接着根据形变信息提取该地震的断层几何特征,包括破裂长度、地表破裂轨迹等参数.最后结合Okada弹性形变模型,反演该地震的震源参数.通过本文的研究将为该地区的地壳应力变化和地震灾害评估提供依据,同时为最新发射的L and sat 8光学影像广泛应用于其他大地震的研究提供参考. 2 数据处理 2.1 光学影像处理

L and sat 8被称为地球资源卫星数据连续性任务(L and sat Data Continuity Mission,LDCM),是美国NASA于2013年2月11日发射的第8颗L and sat系列卫星.与之前的卫星相比,L and sat 8 的光线、热量感应器精准度更高,同时其传感器的空间分辨率也有很大的提升,如B and 8的空间分辨率为15 m,比之前L and sat 7 ETM+的30 m更好地监测地表的变化(USGS,2013c).因此,我们利用L and sat 8第8波段影像监测2013年Mw7.7 巴基斯坦Awaran地震.地震发生之后,我们收集了地震区域震前和震后的L and sat 8数据各两景,轨道(Path):154,分幅(Row):41和42(见表 1),两景数据拼接之后完全覆盖该地震破裂带的全部区域(如图 2所示).数据的气象条件非常好,基本没有云覆盖(影像的整体云量均小于1%),时间相隔26天,保证了震前震后数据的相干性.由于USGS提供的L and sat 8数据经过了几何校正和基于SRTM DEM的粗略正射校正,并投影到UTM坐标系下,但是未提供卫星轨道参数,在本文的数据处理过程中无法进行严密的影像正射校正操作.具体的数据包括三个步骤:

图 2 L and sat 8影像覆盖图.图中蓝色虚线覆盖的范围为L and sat 8影像覆盖,浅紫色五角星分别显示的是2013年9月24号发生的Mw7.7级的主震和2013年9月28日发生的Mw6.8级的余震(地震机制球来自Global CMT). 图中紫色的圆点代表地震发生后3个月内发生大于3级的余震.红线表示通过L and sat8影像配准偏移量算法得到的2013年Mw7.7巴基斯坦Awaran地震的地表破裂线.黑色的线为该区域的断裂分布,该地震就发生在Hoshab fault上Fig. 2 Coverage of L and sat 8 images for the 2013 Pakistan Earthquake. Blue dash line shows the L and sat 8 images data coverage. Light purple stars denote earthquakes occurred on 24 September and 28 September 2013,respectively(focal mechanism is from Global CMT). Purple circles denote aftershocks of magnitude great than Mw3.0 within 3 months after the 24 September 2013 earthquake. Red line represents the fault of 24 September 2013 earthquake retrieved by offset-tracking method using L and sat 8 image. Dark lines are faults in this area and the main shock occurred on the Hoshab fault

(1)影像预处理 :将同一地理范围的影像提取出来,以保证震前震后的影像大小一致.

(2)计算影像偏移量 :利用ROI_PAC软件,采用交叉频谱相关(cross-spectrum correlation)算法计算影像偏移量.配准运算窗口大小设置为64×64像素,偏移量估计步长设置为5×5像素(最终输出的偏移量结果的分辨率约为90 m).

(3)误差后处理 :由于卫星在成像过程中姿态的变化,导致偏移量结果中存在条带误差,我们利用文献(Leprince et al., 2007)提出的Destripe Model去除该误差.此外,由于本文数据处理流程没有进行严密的整体配准,从而导致偏移量存在一个整体性的轨道趋势误差,我们屏蔽掉该地震的形变区域,利 用多项式拟合的方法去掉跟轨道相关的趋势,我们就可以获得该地震的同震形变(Feng et al., 2012).同时,为了进一步去除噪音,我们利用中值滤波进一步降噪.

表 2 文中采用的L and sat 8影像信息表 Table 2 Information of L and sat 8 images used in this study
2.2 地震形变和断层几何特征分析

通过2.1节中的数据处理,可以获得覆盖整个同震形变场的东西向(E/W)和南北向(N/S)的地表形变(如图 3所示).形变沿着地震造成的NEE走向的地表破裂迹线从北向南呈弧形分布.距断层越近,形变越大,上盘形变明显大于下盘,这是走滑型兼逆冲 型地震的主要形变特征.如图 3所示,在东西(E/W)方向,形变主要集中在地表破裂迹的中部和西南部,北部形变较小,最大形变约为7 m;在南北(S/N)方向,形变主要集中在中部和东北部,最大形变达到 10 m.由此可见地表破裂形变从北到南经历了一个 形变逐渐增大,到达最大值后又渐渐变小的过程,因此地表破裂中间部分形变最大.破裂的上盘形变明显大于下盘,说明这次地震的断层向北倾且倾角并非直立,否则上下盘形变会呈对称分布.该特征与表 1中USGS和Global CMT给出的约40°倾角的震源参数结果吻合.两幅影像中远离形变中心的位置仍存在少许形变,这主要是由于影像中仍残余部分与地形相关的误差以及条带和拼接误差的缘故.

图 3 2013年巴基斯坦Awaran地震水平方向的同震地表形变.(a)东(+)西(-)向形变;(b)南(-)北(+)方向形 变.黑色虚线范围用来评定形变结果的精度,为了使结果更有代表性,我们分别在破裂断层的上下部分各选择一个区域,用统计结果的平均值来代表整幅影像结果的精度Fig. 3 Horizontal co-seismic deformation of 2013 Pakistan Awaran earthquake.(a)Deformation in east(+) and west(-)direction;(b)Deformation in south(-) and north(+)direction. The region circled by dark dash line on the right is used to assess the accuracy of deformation calculation results. To make the results be more representative,two parts in upper and lower regions are selected,which are used to calculate the mean value,representing the accuracy of the deformation results

形变结果除了为地震参数反演提供形变观测数据外,还可以在缺少野外调查的情况下,为本次地震提供地震的断层几何参数,如地表破裂轨迹,断层的长度.从该地震的形变信息来看,从北到南沿破裂方向长度约为200 km,地表破裂呈弧形分布,走向角约从240°减少至200°.因此该地震的破裂面用单一或者简单分段的破裂面描述是不够精确的,需要用变走向角的弧形破裂面描述. 2.3 误差分析和精度评定

光学影像偏移量配准技术获取的形变误差来源有很多种,如数据处理误差、地形相关误差、大气误差(云朵覆盖)等.数据处理误差受到配准算法精度限制.我们一般认为光学影像配准最高可以达到 1/50个像素(Leprince et al., 2007),对应到L and sat 8第8波段影像而言,15 m空间分辨率对应30 cm的精度水平,这是误差的主要来源.大气误差一般指影像受云朵的影响,从而导致影像不能真实反映地表特征.这些区域容易识别并且由于配准相干性较低往往会被屏蔽掉,因此对结果的影响可以忽略不计.跟地形相关的误差是由于影像的正射纠正误差引起的.我们可以获取的L and sat 8影像是经过各种纠正过后的1级产品,无法借助外部DEM数据进行地形纠正,因此该误差也是一个重要源头,误差水平跟DEM的高度有一定的相关性.另外,影像的拼接误差,沿方位向的条带误差对于个别影像而言可能也比较明显.

对于每一对影像而言,各种误差来源对其总误差的贡献量也是不一样的.同时由于数据处理误差和地形相关误差在空间范围内是相关的,因此需要建立形变误差的方差-协方差结构函数来对不同影像对进行精度评定,同时获取同一幅形变图中不同像素点处误差的相关关系(周辉,2012).误差结构函数作为衡量形变结果精度的指标,可以用来对不同的数据源(如东西向形变和南北向形变,不同卫星传感器数据等)进行适当定权来优化模型参数反演.

本文采用自协方差函数来描述形变结果的误差分布.在东西和南北方向形变图中选取非形变区域,为了使选取的区域更有代表性,我们分别在断层的上下盘各选取一个区域(图 3中黑色虚线所围部分).由于图 3形变结果的长波段形变已经扣除,因此认为该区域均值接近于0,假定该区域误差分布满足二阶平稳特性(即该区域的误差分布特征可以代表整个形变图的误差分布特征),此时任意两点之间的协方差只与距离有关,与两点之间的方向、位置无关.本文根据地统计学中的变异函数理论(周辉,2012),分别采用变差函数和协变差函数对形变场取样,进而估算两点之间的空间相干性.对任一距离hi来说,变差函数通常可以表示为

其中,xi表示影像中随机采样点位置,f(xi)表示在影像随机采样点xi上的形变值,hi为距离间隔,N表示给定距离hi时采样点对的个数.同样,协变差函数可以表示为

在选取的非形变范围1和2内,针对每个距离h,随机选择1000个点对个数,根据公式(1)、(2)分别计算其变差函数和协变差函数.当γ(h)趋于稳定时,该值称为基台值,表示该形变场的方差大小.文中使用一个非负协方差函数来拟合协变差函数C(h)=a·e-h/b,设函数模型分别代表形变图的误差分布(图 4),最终估计得到的东西和南北方向的协方差函数如表 3所示.
图 4 形变结果的精度分析图.图中实线分别代表东西和南北方向形变图中范围1和范围2的协方差函数.另外通过该范围的形变计算了各个区域的中误差φ.图 4表 3明显反映,南北方向形变的误差水平大于东西方向的 Fig. 4 Uncertainty of deformation results. Solid lines represent fitting covariance functions of E/W(east/west) and N/S(north/south)directions in area 1 and area 2,respectively. The mean square error(MSE)φ is calculated for each area. This figure and Table 1 show that the error level in N/S direction is obviously greater than that in E/W direction

图 4表 3可以看到,东西向形变的协方差 在距离大于10 km时均逐渐趋于0,表示影像中距离超过10 km的两个像素点之间是不相关的.南北向形变的协方差在距离大于20 km时才逐渐趋于 0,说明在南北方向有非常明显跟空间相关的误差.这个跟我们选取的区域是吻合的,我们在南北向形变场中选取的区域有非常明显的跟地形相关的误差,因此反映在协方差函数上空间相关性就明显大于东西方向形变.由表 3可知,南北方向的中误差大于东西方向的,前者中误差约为59.8 cm,略大于后者的中误差(约为49 cm).

表 3 方差和协方差估计结果参数列表 Table 3 Parameters of variance and covariance estimation
3 地震参数反演 3.1 数据降采样

由于偏移量配准技术获取的形变信息是空间连续的,对于这次地震而言,仅仅东西向形变数据点个数已达到107级,再加上南北向形变数据,总的数据点个数达到千万数量级.如果全部数据点都参与反演,会大大增加反演的时间,降低反演效率.另外,形变场中大量噪音点也会造成计算结果很难收敛,因此有必要对明显噪音区域预处理.本文提出改进的四叉树算法有以下特点:(1)数据的预处理.对于东西和南北向形变场中明显的噪音区域,比如大海区域,山脊区域(地形相关误差非常明显),余震形变区域等进行掩膜处理.该预处理将可以起到数据质量控制的作用,使噪音区域不参与四叉树的降采样处理,这样也避免了四叉树结果在局部区域和边界区域出现异常“聚集”现象;(2)控制窗口大小的降采样.在传统四叉树处理的基础上添加最大窗口和最小窗口控制条件,使得在超过窗口上下阈值的情况下进行均匀采样,这样处理可能会在近场区域引入一定的误差,但它使采样结果的分布更为合理(如图 5所示),更重要的是通过该处理可以避免传统四叉树方法在远场区域采样点数过少和断层附近区域采样过密的缺陷,使采样点的分布和点的数目有一个合理的折中.在本次地震中,选取中误差阈值为55 cm,控制最大窗口为128×128像素,最小窗口为8×8像素,在南北和东西形变场中都提取出2268个数据点参与地震参数的反演,数据个数仅占原来数据总数约0.1%(如图 5所示).

图 5 基于改进四叉树算法的数据降采样.(a)南北方向形变改进的四叉树降采样图;(b)形变的矢量图.四叉树算法选用最小窗口为8×8,最大窗口为128×128,左图中的空白区域因为噪音较大,在四叉树降采样之前已经被屏蔽,未参与降采样处理.经过本文改进的四叉树降采样后,形变观测的数据点降到了2286,减小到原来数据点总数的0.1%,大大提高 了数据的反演效率Fig. 5 Data down-sampling based on modified quadtree partition.(a)Map of modified down-sampling quadtree in N/S direction;(b)Vector graph of deformation. The minimum and maximum sizes of windows are 8×8 and 128×128,respectively. The blank area in left subplot is masked before quadtree resampling procedure due to high noise level. After down-sampling process,the number of data points is reduced to 2286,which is only 0.1% of original dataset, and this will improve the efficiency of inversion
3.2 非线性震源参数反演

在地震的震源参数反演中,一般用9个参数来描述地震参数模型,即断层长度、宽度、深度、倾角、走向角、中心经纬度以及沿走向、倾向的滑动量.这些震源参数和同震形变之间是非常复杂的非线性关系.为了获得更好的参数估计效率和精度,我们往往会借助地震调查获得的地表破裂信息,初步约束断 层起始点的经纬度、走向、断层的长度等显性参数(可以直接通过野外调查或者高分辨率光学影像获取的参数),然后用非线性反演方法估计断层倾角、断层宽度等隐性参数.由于非线性反演过程运算量较大,反演过程中往往选择简单的断层模型来限制反演过程中需要估计的断层参数数目.

在本文中,我们选取单一的断层模型(Okada,1985),沿着断层走滑方向给出断层参数的初始值,然后再约束这9个参数的大致解空间.首先利用模拟退火算法,获取估计参数的初始最优值,然后基于该初始值用条件约束规划算法在约束空间里面获得待估参数的最终最优结果.通过模拟退火算法,我们可以获得接近结果的最优值,然后再用约束规划算法最终获得整个解空间的最优结果,这样可以很好地克服局部最优解,做到“双重保障”(Feng et al., 2010; 周辉等,2013).通过该非线性反演组合,我们发现该地震的最优倾角为47°,走向角为223°,这与USGS和Global CMT的地震机制解吻合. 3.3 断层滑移量反演

基于3.2节非线性反演确定的断层参数,沿着地表破裂轨迹(断层的走向)建立复杂的弧形断层模型.新建立模型的断层单元的走向是沿地表破裂方向变化的,从北到南走向角约从240°减少至200°,断层的长度为200 km,宽度为35 km,深度约为25 km,延伸到地表,并把断层按7 km×7 km划分为34×5个子单元,分别求解每个子单元上的滑动量,最终获取整个断层面的滑动分布.这样就大大增加了待求解参数的个数,为了避免解空间的“震荡”,需要在空间上对滑移量进行一定的平滑约束.同时,为了提高反演结果的收敛率,固定断层倾角为47°,介质泊松比设为0.25.考虑到本文使用的光学影像偏移量形变在东西向和南北向的形变精度水平不一样,在反演过程中,我们利用2.3节获取的中误差水平进行定权.最后,根据模型粗糙度与光滑度之间的 折中曲线,我们选取的粗糙度为4.1 cm/km(Jonsson et al., 2002).参数设定完成后,我们反演出2013年巴基斯坦Awaran地震的断层滑动分布结果如图 6所示. 4 结果与讨论

图 6所示,断层面的滑动分布比较集中,主要位于沿走向从北到南方向70~180 km、倾向距地表 15 km以内的区域,最大断层滑动分布接近地表.在该区域内,断层滑动量存在两个峰值,第一个峰值位于断层中部,靠近地表,最大滑动量达10 m,第二个峰值出现在断层面靠近南端约40 km靠近地表处,最大滑动量达5 m.整个断层面滑动量以左旋滑动为主,在滑动量比较集中的区域也兼有少许倾滑量,且数值较小,最大倾滑量约为2.5 m,远远小于最大走滑量9.0 m.本文在反演过程中没有固定滑动角,结果显示平均滑动角约为10°,与GlobalCMT和USGS给出的滑动角较为接近(如表 1所示).

图 6 基于L and sat 8形变观测反演的走滑和倾滑同震滑移量分布结果.图中每个单元代表 km×7 km,结果显示该地震是左旋走滑为主兼有少量倾滑分量的地震,滑移量主要集中在靠近地表 0~15 km区域,最大滑移量达到10 m,靠近 地表,断层靠近南边也有一个小的峰值,走滑形变分量达到5 mFig. 6 Strike-slip and dip-slip deformation of the preferred model. Each rectangular element denotes an area of 7 km×7 km. The results show that strike-slip dominates the earthquake fault slippage concentrating in the depth above 15 km,with a little dip-slip component. The maximum slip is up to 10 m. Near ground surface,a small slip peak is located south near the fault line and its strike-slip deformation is up to 5 m

通过本文反演得到的断层滑动分布模型可知,这次地震的滑动量主要是以左旋走滑为主,兼具少许倾滑分量.该地震相对其他大地震,如1999年集集地震和2008年汶川地震(Feng et al., 2010)来说,断层破裂方式较为特别(呈弧形破裂轨迹),本文采用真实断层走向的曲面断层模型,因此可以最大限度地模拟真实地震的地表形变.本文的数据拟合残差约35 cm,该误差水平跟光学影像的形变监测能力是吻合的.根据反演获得的模型分别正演了在东西和南北向的形变值,并计算了模型的残差值(如图 7所示).在残差图中(如图 7c和d7),我们可以看到的误差源有条带误差、跟地形相关的误差(部分数值达到0.5 m),另外也可以看到非常明显的拼接误差.这里需要说明的是,在3.1节数据降采样处理过程中,我们已经把这些明显误差区域给屏蔽了,因此这些误差源对模型的影响是非常有限的,因此这些误差可以非常直观地反映到模型的残差图中.比较东西和南北向的残差图(图 7c和7d),我们可以发现,跟地形相关的误差主要反映在南北向上,因此在将来利用L and sat 8进行地表形变监测时,研究如何对地形引起的误差进行改正也是非常有必要 的.本研究中参与反演的数据拟合度达98%,证明本文得到的断层滑动分布结果具有很好的可靠性.本文 选取介质的剪切模量μ=30 GPa,根据震级计算公式,计算得到该地震的标量地震矩为 4.68×1020N·m,震级为Mw7.75级,介于Global CMT和USGS给出的矩震级之间.

图 7 形变数据拟合残差图.(a)和(b)分别为东西和南北方向模型拟合形变结果;(c)和(d)分别为东西和南北方向形变残差图
Fig. 7 Model prediction and residual for the L and sat observation.(a) and (b)are model prediction of deformation
in E/W and N/S direction,respectively.(c) and (d)are the residuals

Chaman fault是亚欧板块和印度板块之间碰撞的转换断裂,地质野外调查资料和高分辨率遥感影像显示该断裂平均滑移为24~35 mm/a(Ul-Hadi et al., 2012),通过GPS速度场观测也显示 平均滑移约为 18±1 mm/a(Ul-Hadi et al., 2012),因此不管是地质观测手段还是大地测量的手段都显示该断裂上累计了大量的应力.另一方面,发生在该断裂上 大于Mw6.0的地震非常有限(Ul-Hadi et al., 2008; Ul-Hadi et al., 2012),自1892年以来有记录的只 有三次,分别是1892年12月20发生在Chaman的 Mw6.8 Chaman地震,1935年5月31发生在Balochistan 的Mw8.0 Quetta地震(如图 1所示),近年野外调查研究发现该震级估计偏大,后被修正为Mw7.7,与本次地震的震级接近(Ul-Hadi et al., 2008),以及1978年3月16发生在Quetta地震北边(29.83°N,66.43°E)的Mw6.1地震.因此有学者认为在Chaman fault靠近巴基斯坦北部会有一个 “地震空区”(Seismic gap)(Ul-Hadi et al., 2012),预测该地震空区内将有大于Mw7级地震发生.同时历史地震记录显示发生在该断裂上5级左右的中小型地震较多,有学者通过InSAR监测了 2005年10月21发生在Chaman fault的一个Mw5.0 的震后形变,发现这次地震的震后形变以8 mm/a的速度持续1~2年,震后形变释放的能量是同震形变的8倍,因此解释了通过这种中小地震加上震后形变的模式释放了大量的能量,从而减小或者推迟了Mw7以上地震的发生(Furuya and Satyabala, 2008).主震发生后第4天(2013年9月28号),在主震的东 北方向靠近破裂的起始点位置又发生了一个Mw6.8 的余震,该余震跟主震的震源机制非常相似.另外,该次地震的余震主要聚集在破裂东北边缘(如图 2所示),衔接该分支断裂和Chaman fault主断裂的区域.同时,我们发现这处地震并没有沿着东北缘的Chaman fault主断裂破裂,而是向西南沿着弧形方向的Hoshab fault发生,足以说明该断裂在东北方向受到强烈阻挡,将继续累积更多应力,外加本次地震对Chaman fault主断裂的应力加强,因此使得该区域潜在地震灾害加大.从历史的地震记录来看,在 该次地震和1935年Mw8.0 Quetta地震之间约300~400 km 的区域,在过去100多年间都没发生过7级以上地震,因此进一步证明该次地震的东北缘沿着Chaman fault的主断裂存在重大的地震隐患,需要继续保持监测和警惕.

2013年Mw7.7巴基斯坦Awaran地震并非发生在前人(Ul-Hadi et al., 2012)研究指出的地震空区内,而是发生在目前了解甚少的Chaman fault以南的分支断裂Hoshab fault上,该断裂深部连接莫克兰俯冲带,地震发生前地质和地球物理学者一般认为该区域将会有较大倾滑分量,然而这次地震以左旋走滑为主,仅有少量倾滑分量.因此这次地震的发生将有助于对该区域的重新认识.另外一个值得大家思考的问题是传统的走滑型地震都是以高倾角型为主(>80°),如Mw7.6玛尼地震和Mw8.1可可西里地震,而本次地震倾角47°,是低倾角型的走滑型地震,这在世界范围内也是少见的,因此本次地震将对于我们研究Chaman fault转换带和莫克兰俯冲带之间的转换关系提供全新的视角. 5 结论

本文利用美国USGS提供的L and sat 8数据,基于光学影像偏移量配准算法,获取了2013年巴基斯坦Awaran地震的同震形变场和地震地表破裂几何特征.根据获得的东西和南北向形变的误差分布特征,使用地统计学中的变异函数模型对形变结果进行精度评定.最后结合Okada弹性形变模型和形变的精度水平反演了该次地震的断层滑动分布.结果显示该地震断层滑动以左旋走滑为主,兼具少量的倾滑分量.滑动主要集中在断层面距地表 15 km以上区域,靠近地表最大滑动量达10 m,另外在在弧形断层面南部也有一个小峰值,最大达5 m.反演获得的标量地震矩为4.68×1020N·m,震级为Mw7.75级.本文结果是采用L and sat 8光学影像对巴基斯坦地震震源参数的反演,结果揭示了该地区断层活动依然活跃,应采用多种观测手段对震后和震间形变进行观测,提高防震减灾的可行性.

致谢 本研究所用L and sat 8数据来源于美国USGS;文中所有图件使用GMT5.0软件绘制而成.
参考文献
[1] Caltech, 2007. 《The Caltech COSI-Corr package》,http://www.tectonics.caltech.edu/slip_history/spot_coseis/download_software.html.
[2] Furuya M, Satyabala S P. 2008. Slow earthquake in Afghanistan detected by InSAR. Geophysical Research Letters,35: L06309, doi: 10.1029/2007GL033049.
[3] Feng G C, Ding X L, Li Z W, et al. 2012. Calibration of InSAR-derived coseimic deformation map associated with the 2011 Mw9.0 Tohoku-Oki earthquake. IEEE Geosci. Remote Sens. Lett., 9(2): 302-306.
[4] Feng G C, Hetland E A, Ding X L, et al. 2010. Coseismic fault slip of the 2008 Mw7.9 Wenchuan earthquake estimated from InSAR and GPS measurements. Geophys. Res. Lett., 37: L01302.
[5] Global CMT, 2013. 《Global CMT catalog earthquake》,http://www.globalcmt.org/CMTsearch.html[2014-01-06].
[6] Jonsson S, Zebker H, Segall P, et al. 2002. Fault slip distribution of the 1999 Mw7.1 Hector Mine, California Earthquake, estimated from satellite radar and GPS measurements. Bull. Seism. Soc. Am., 92(4): 1377-1389.
[7] Leprince S, Barbot S, Ayoub F, Avouac J P. 2007. Automatic and precise ortho-rectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. IEEE Transactions on Geoscience and Remote Sensing, 45(6): 1529-1558.
[8] Okada Y. 1985. Surface deformation due to shear and tensile fault in a half-space. Bull. Seism. Soc. Am., 75(4): 1135-1154.
[9] Sciencenet, 2013. 《Introduction to 2013 Pakistan Earthquake》,http://blog.sciencenet.cn/blog-51597-727793.html[2014-01-06].
[10] Shuhab Khan, 2013.《Earthquake Hazards and Studies of the Chaman Fault》,http://www.uh.edu/-sdkhan/research_chaman_fault.html [access 2013-11-26].
[11] USGS Earthquake Hazards Program. 2013a. M7.7-61 km NNE of Awaran, Pakistan, 2013-09-24 11:29:47 UTC [Online]. Available: http://comcat.cr.usgs.gov/earthquakes/eventpage/usb000jyiv summary[2013-10-25].
[12] USGS Earthquake Hazards Program. 2013b. M6.8-85 km NNE of Awaran, Pakistan, 2013-09-28 07:34:06 UTC [Online]. Available: http://comcat.cr.usgs.gov/earthquakes/eventpage/usb000k1gb summary USGS Landsat 8 dataset. 2013c. [Online]. Available: http://earthexplorer.usgs.gov[2013-10-25].
[13] Ul-Hadi S, Khan S D, Owen L A, et al. 2008. Slip-rates along the Chaman fault: Implication for transient strain accumulation and strain partitioning along the western Indian plate margin. Tectonophysics, 608, P389-400.
[14] Ul-Hadi S, Khan S D, Owen L A, Khan A S, 2012. Geomorphic response to an active transpressive regime: a case study along the Chaman strike-slip fault, Western Pakistan. Earth Surf. Process. Landforms, 38: 25-264.
[15] Wikipedia, 2013a. 2013 Pakistan earthquake Wiki introduction. [Online]. Available: http://en.wikipedia.org/wiki/2013_Pakistan_earthquake[2013-10-31].
[16] Wikipedia, 2013b.Chaman fault. 2013 [Online]. Available: http://en.wikipedia.org/wiki/Chaman_Fault[2013-10-31].
[17] Zhou H, 2012. Study of Burma Earthquake co-seismic deformation and its focal shock parameter by using InSAR data [Master's thesis] (in Chinese). Changsha: Central South University.
[18] Zhou H, Feng G C, Li Z W, et al. 2013. The fault slip distribution of the Myanmar Mw6.8 earthquake inferred from InSAR measurements. Chinese J. Geophys. (in Chinese), 56(9): 3011-3021."
[19] 科学网,2014. 2013年巴基斯坦地震的介绍[Online]. http:∥blog.sciencenet.cn/blog-51597-727793.html\[2014-01-06\].
[20] 周辉, 2012.利用InSAR资料研究缅甸地震同震形变和震源参数,长沙: 中南大学硕士学位论文.
[21] 周辉,冯光财,李志伟等.2013. 利用InSAR资料反演缅甸Mw6.8地震断层滑动分布.地球物理学报, 56(9): 3011-3021.