据中国地震台网中心测定，北京时间2016年11月25日22时24分在新疆维吾尔自治区阿克陶县发生6.7级中等强度地震，震中位置(74.04°E，39.27°N)，震源深度约11 km(http://www.csi.ac.cn/manage/eqDown/05LargeEQ/201611252224M6.7/zonghe.html).由于地震发生于高海拔地区，人烟稀少，且当地民众居所总体抗震性能较好，未造成重大人员伤亡，仅一人因房屋倒塌死亡，数人受伤.中国地震局启动Ⅲ级应急响应，除现场派出应急救援力量外，各研究机构也纷纷开展震源研究工作.中国地震局地质研究所由陈杰研究员组队赴震中区开展地震地质科学考察，科考组克服高寒积雪等不利因素，进行了野外调查和测量.初步分析认为2016年11月25日新疆阿克陶6.7级地震主震发生在NWW走向的木吉右旋走滑断裂上，并在距中国地震台网震中SEE方向约2.4 km处以及以东32.6 km两处分别发现了两处小规模右旋张扭性地震地表破裂带，但沿所有断层陡坎均未发现地表破裂及相关地质灾害(陈杰等，2016).该区域历史上强震活跃，方圆200 km范围内曾记录到多次中强地震，图 1a中给出了从全球质心矩张量项目(global CMT)(Ekström et al., 2012)中检索到的1976年至今16次6级以上地震的震源机制解(不包括本次地震).其中1974年8月11日乌恰M7.3地震(该次地震的震源机制解根据Burtman等人的研究成果)及2008年10月5日M_W6.7地震距离此次震中最近，均发生在其西北方向约20 km处东西走向的弧形构造上；而时间最近的一次中强震是2016年6月26日吉尔吉斯斯坦M_W6.4地震，在此次震中西北约60 km处(图 1b).地震具有由北向南迁移的趋势，是否对本次地震有触发关系？

图 1

Fig. 1

图 1 帕米尔东北缘及邻区构造图 KYTS—喀什叶城转换带，MPT—帕米尔主断裂带，PFT—帕米尔前缘逆冲推覆系，DF—达瓦孜断裂, KF—喀喇昆仑断裂，KAF—喀喇喀什(康西瓦)断裂，Muji fault—木吉断裂，Kongur fault—公格尔断裂系北段. (a)背景为SRTM，红色线条为区域内主要断层，蓝色箭头为GPS速度场(相对于欧亚大陆，数据来自Zhou等(2016))，黄色五星为本次震中，白色方框为图1b范围，黑色虚线框分别为Sentinel-1A升、降轨数据覆盖范围，黑色沙滩球为200 km范围内发生的大于6级地震震源机制解，红色圆圈为震级大于5级、震源深度小于60 km的浅源地震，紫色圆圈为震级大于5、震源深度小于300 km的中源地震; (b)Landsat 8遥感背景图，红色线条为近场断层，红色、黄色、蓝色沙滩球为本次阿克陶地震震源机制解，分别来自USGS、GCMT和CENC.黑色沙滩球为1970年以来本区域内发生的大于6级地震震源机制解.彩色圆点为余震分布，不同颜色代表不同深度; (c)余震分布剖面图，AA′和BB′剖面线的位置见图1b，此处不同颜色代表余震震级大小. Fig. 1 Tectonic map of the northeastern margin of Pamir and its adjacent area KYTS—Kashgar-Yecheng transfer system, MPT—Main Pamir Thrust, PFT—Pamir Frontal Thrust, DF—Darvaz fault, KF—Karakorum fault, KAF—Karakax fault. (a) Topography, faults in red line, GPS velocity in blue, yellow star denotes epicenter. The black dashed box is the coverage of Sentinel-1 data. The white rectangle presents the region in Fig. 1b; (b) Aftershock and focal mechanism, the red beach ball is from USGS, the yellow one is from GCMT and the blue one from CENC; (c) Profile of the aftershock, the location of the profile AA′, BB′see Fig. 1b.

据张勇等震后不久发布的震源破裂过程研究结果，破裂过程为单侧破裂，从震中开始向东传播，断层面上滑动分布主要分成两个部分，且断层东部破裂已经到达地表，同震滑动量最大可达30 cm(http://www.cea-igp.ac.cn/tpxw/275080.html).美国地质调查局(USGS)给出的地震初始破裂点为(74.021°E, 39.230°N)，震源深度17 km，矩震级为M_W6.6，震源机制解显示为1次具有拉张分量的右旋走滑错动事件(见表 1).中国地震台网中心(CENC)给出的震源机制解为右旋走滑兼具正断作用.显然，不同机构给出的地震破裂起始点位置、断层面解均有所不同，震中附近构造变形极为复杂，前面提到的发生在该地区附近的几次历史地震的震源情况均较为复杂，往往由多个子事件组成(Burtman and Molnar, 1993)，从已经发布的反演结果来看(Wang et al., 2017)，滑动发布中存在着一个破裂空区(slip gap)，那么本次地震是否也是由多次地震事件组成？厘清本次地震的发震断层及震源破裂细节重要而急迫，对深刻理解帕米尔高原现今构造变形状况以及动力学演化过程等问题也具有重要的科学意义.

表 1 Table 1

表 1 不同机构给出的阿克陶地震震源机制解 Table 1 Source parameter of the Aketao earthquake from different agency 来源 震中位置 节面Ⅰ 节面Ⅱ 震级 最佳深度 (km) 经纬度(°) 走向/倾角/滑动角(°) 走向/倾角/滑动角(°) Harvard 74.13 39.29 199/87/3 109/87/177 MW6.6 17.3 USGS 74.021 39.230 199/84/14 107/76/174 MW6.6 17.0 CENC 74.01 39.27 113/54/-168 16/80/-36 MW6.6 11

表 1 不同机构给出的阿克陶地震震源机制解 Table 1 Source parameter of the Aketao earthquake from different agency

本次地震发震区域交通不便、研究基础薄弱，且无明显地表破裂，发震构造信息尚不十分明确.InSAR可以获得大范围高空间分辨率的地表形变信息(Gabriel et al., 1989; Massonnet et al., 1993), 利用大地测量数据进行反演将有助于更好地认识地震滑动分布和断层构造特征(Dawson et al., 2007; Funning et al., 2005)，可以为浅源地震的深度提供很好的约束.研究滑动分布特征不仅可深入理解阿克陶地震的运动学特征，还可以为区域地震危险性评估提供依据.本研究同时获取了Sentinel-1A升降轨及Alos2升轨的同震形变场，采用一种结合先验知识的多视角D-lnSAR最小二乘迭代逼近解算法，求解出了三维形变场.又利用观测数据为约束反演得到了发震断层面上的精细同震滑动分布.在此基础上本文分析了2016年阿克陶地震的发震构造特征并给出其构造模型，讨论了发生于弧形构造带上地震的空间关系.

由印度板块向欧亚板块碰撞形成的喜马拉雅造山带的东、西两端分别有一个构造急剧转向的地区——构造结，是构造应力作用最强、隆升和剥露速率最快及新生代变质作用最强的地区(Burg and Podladchikov, 1999)，本次地震就发生在西构造结的前缘推覆构造上.位于青藏高原西北部的西构造结主要由向北突出的帕米尔构造弧组成.该构造弧北以阿莱山谷(Alai Valley)、南天山(South Tien Shan)为界, 在南北方向上由主喀喇昆仑逆冲断裂带(Main Karakoram thrust，MKT)、主帕米尔断层(Main Pamir Thrust，MPT)、帕米尔前缘断裂带(Pamir Frontal Thrust，PFT)等弧形活动褶皱逆断裂带及其间的推覆构造构成(图 1).频繁发生的中强地震沿活动褶皱逆断裂带形成了南北两条极为显著的中、浅源地震带(Burtman and Molnar, 1993; Fan et al., 1994).这2个地震带被认为是2个俯冲岩石圈板片的所在位置, 即欧亚板块以小角度向南俯冲于帕米尔之下，而印度板块以高角度俯冲于新都库什(Hindu Kush)之下(Fan et al., 1994; Burtman and Molnar, 1993).主帕米尔逆断裂带(MPT)弧型推覆构造在东西两侧由正向逆冲渐变为斜冲走滑大断裂，显示出东侧右旋、西侧左旋的断层活动特征.东缘以喀喇昆仑断裂(Karakorum fault, KF)、喀拉喀什断层(Karakax fault, KXF)和喀什—叶城走滑系统(Kashgar-Yecheng transfer system, KYTS)与西昆仑构造带和塔里木盆地为界(Cowgill，2010; 胡建中等，2008)，其中喀喇昆仑断层(KKF)总体走向N45°W，断层北段的右旋走滑滑动速率在6.9~10.8 mm·a^-1之间(Robinson，2009)；而喀拉喀什断层又称康西瓦断层，位于喀喇昆仑断层以东，它构成了青藏高原的北边界，Lin等(2008)通过卫星影像解译及在康西瓦附近的野外考察，认为该断层是一右旋走滑断层，走滑速率8~15 mm·a^-1，垂直滑动速率~2 mm·a^-1；最东缘的喀什—叶城转换带被认为是帕米尔和塔里木盆地的主要边界断层，吸收调节了两大块体的相对运动，晚新生代以来平均滑动速率11~15 mm·a^-1(Cowgill，2010).西则以达瓦孜断裂(Darvaz fault, DF)与塔吉克盆地(Tadjik depression)为界, 但对该断层的运动学特征至今未能很好确定，Burtman和Molnar(1993)认为断层左旋走滑速率为10~15 mm·a^-1.

在晚新生代，帕米尔的变形以前缘正向逆冲、两侧斜冲-走滑和内部拉张为特征.约50 Ma BP以来, 印度板块与欧亚板块的碰撞和持续的汇聚作用造成中亚地区强烈的陆内变形, 帕米尔构造带向北的逆冲推覆造成地壳的强烈缩短增厚，同时导致左右两侧的喀喇昆仑等断裂大规模右旋走滑、达瓦孜左旋走滑，而形成帕米尔弧形构造带(Burtman，2000; Strecker et al., 1995).现今的GPS观测资料表明，前缘地壳缩短、构造结内部拉张的运动趋势仍在继续.Zhou等(2016)利用2011年在帕米尔高原建立的一横一纵两条GPS测线，再结合Zubovich等(2010)、Ischuk等(2013)、Jouanne等(2014)前人的观测资料形成两条GPS剖面.南北向的剖面纵贯帕米尔高原不同的地质构造单元，旨在观测南北向的地壳缩短；东西向的剖面从塔吉克盆地到塔里木盆地，横穿整个帕米尔高原，以观测弧形构造内部的拉张.GPS速度场整体向北，但不同构造带间存在明显的挤压速率差异，帕米尔北缘断裂的断层挤压速率在~12 mm·a^-1, 而在帕米尔高原和塔里木盆地之间的东西向拉张在~10 mm·a^-1，而在帕米尔高原与塔吉克盆地之间，则主要在塔吉克盆地西部存在5~7 mm·a^-1的地壳缩短(Zhou et al., 2016).

帕米尔内部的构造变形以近东西向拉张伸展作用为主，具体表现形式为在东、西部分别发育有公格尔伸展系统和卡拉库尔地堑(Amidon and Hynek, 2010; Brunel et al., 1994; Strecker et al., 1995).其中东部的公格尔伸展系统长约250 km，走向NW—SE，是1条全新世活动断裂(Chevalier et al., 2011, 2015; 李文巧等，2011; 陈杰等，2011)，由北至南可分为木吉断层、公格尔正断层、塔合曼正断层和塔什库尔干正断层4段，东西向的伸展量在公格尔山附近最大，约为35 km(Robinson et al., 2007; 陈汉林等，2014).本次地震即发生在公格尔伸展系统北部的最北端的转换断层——木吉断层上(陈杰等，2016).关于该断层的运动特征，看法尚不一致.吴传勇等(2008)野外考察通过对断错地貌的测量及地貌面年龄的测定，认为木吉断层断面倾向NE，倾角50°~65°，为1条逆冲右旋走滑断层，全新世中期以来断层的平均垂直滑动速率为1.8~2.1 mm·a^-1.另外大部分学者认为木吉断层走向NWW，近直立，右旋走滑为主，兼具正断作用.从断层沿线发育壮观的断层三角面(陈杰等，2011; Robinson et al., 2004)可以支持这一论断，河流阶地测量位错量再通过宇宙核素测年后表明木吉断裂的右行走滑速率约为4.5~11 mm·a^-1 (Chevalier et al., 2011).断层走向在木吉盆地东南部突然发生改变，变为公格尔断裂公格尔段北部(又称昆盖山南麓断裂北段)起点.该段断层走向近SN，断面W倾，沿昆盖山山前发育，错断了一系列全新世冲沟、洪积扇或冰碛地貌(陈杰等，2011).

根据房立华博士(私人通讯)提供的余震精定位结果(时间截止到2016年11月26日18时5分，共有173条余震记录)，余震几乎全部分布在木吉断裂南侧，主要分布在地下3~15 km，余震带走向在主震位置处发生变化，从形态上可以分为东西两个部分.在主震以西近EW走向，余震剖面BB′(图 1c)显示余震较密集且以0轴大致呈对称分布，表现出走滑型地震余震分布特征；而在主震以东余震带折向SE，余震剖面AA′(图 1c)显示余震沿倾向SW的断层面分布.从余震分布剖面图来看，似乎本次地震在主震东西两侧的运动性质不同，是否是由两次事件引起？

Sentinel-1A雷达卫星在刚发射升空不久就捕捉到了美国纳帕6.1级地震(李永生等，2015)，近两年更是在尼泊尔M_W7.8(Elliott et al., 2016)、智利IllapelM_W8.3(Zhang et al., 2016)等大地震中发挥了重要作用.

本文选取了一对ALOS2升轨数据(T162)、Sentinel-1A数据升降轨各一对作为干涉像对，具体参数见表 2.采用GAMMA软件(Werner et al., 2000)来获取本次地震的同震形变场.ALOS2数据成像模式为精细模式，处理步骤同常规的差分干涉.而Sentinel-1A数据为Terrain Observation with Progressive Scans (TOPS)成像模式，覆盖范围广，数据量很大，而形变场只占影像的很小一部分.因此对于降轨T107的数据，只选取覆盖震中2个subswatch中的6个burst进行成像，而对于升轨T027的数据，震中上、下两景均只覆盖了一部分的形变场，只能从相邻两个subswatch中各选取3个burst，计算burst重叠区域的干涉图之差，并将差分相位分配到burst上，消除相位跳变，利用有效行与列拼合成一幅图像进行干涉处理.Sentinel-1A的TOPS成像模式造成SAR影像方位向多普勒频率变化较大, 干涉数据处理对影像配准要求很高(Yagüe-Martinez et al., 2016).为满足这一要求，除使用了Sentinel-1A精密轨道数据外，还采用了DEM辅助的影像配准过程(Nitti et al., 2008)，利用亚像素追踪(subpixel tracking)方法进行迭代精炼，最终配准精度优于0.001个像素.配准和地形相位去除中使用了Shuttle Radar Topography Mission (SRTM) V4.1(Jarvis et al., 2008).

表 2 Table 2

表 2 雷达干涉像对数据参数 Table 2 SAR data used in the InSAR analysis Satellite Track Inclination Master (YYYYMMDD) Slave (YYYYMMDD) ΔT (days) B⊥ (m) Incidence (°) Sentinel-1A T107 Descending 20161125 20161219 24 95 39.19 Sentinel-1A T027 Ascending 20161020 20161207 48 70.6 33.72 ALOS-2 T162 Ascending 20160217 20161207 294 68 36.28 注：每列分别为卫星飞行方向(升轨、降轨)、震前震后影像获取日期、时间间隔、垂直基线和入射角. Note: Columns show satellites directions (ascending or descending), pre- and post-seismic acquisition dates, time interval ΔT, and perpendicular baseline B⊥, incidence.

表 2 雷达干涉像对数据参数 Table 2 SAR data used in the InSAR analysis

最终得到的干涉形变场如图 2所示.干涉条纹沿木吉断裂展布，木吉断裂北侧山区地形陡峭导致失相干，但仍能分辨出形变的量值和形变覆盖范围.断裂南侧木吉盆地相干性较好，条纹清晰且较为密集，在CENC所定震中东、西两部分形成两个干涉条纹密集区，西部有四个条纹，分布范围较广，形变梯度小，而东部也有四个条纹，但条纹更加密集，ALOS2由于波长较长的缘故，只能体现出两个形变中心，且左边的成扫帚放射状，而东部的呈椭圆挤压状.在这三个轨道干涉像对条纹图的东西两部分各拉一条剖面线(位置见图 2a, 2b, 2c, 其中T027在两个形变椭圆中间地带多拉了一条剖面GG′)，西部形变椭圆的剖面线AA′和CC′(图 2d，2e)坡度减缓，最大形变10 cm，说明断层破裂较深，影响范围较大，而BB′和DD′更为陡峭一些，最大形变在12 cm，说明在这部分断层面深度较浅，影响范围较集中.升轨的T027在断层南侧区域视线向(Line of Sight，LoS)形变值为正值、朝向雷达传感器的方向运动，降轨的T107在断层南侧LoS向值为负值、远离传感器的方向运动，结合升、降轨雷达观测右视成像特点，据此可判断断层南盘相对北盘向西运动，表现为右旋走滑特征，这与木吉断裂的活动性质是一致的.另外断层南盘形变场条纹比北盘条纹密集，可以初步判断该断层南倾，但倾角大小则需要在滑动分布反演时确定.

图 2

Fig. 2

图 2 D-InSAR LoS向形变场及剖面图 Sentinel-1A数据每个条纹色周代表2.8 cm形变值，ALOS2每个条纹代表11.8 cm.(a) Sentinel-1A升轨T027；(b) Sentinel-1A降轨T107；(c) ALOS2升轨T162.红色、绿色散点分别为经过震中西部、东部的剖面线缓冲区1 km范围内的点，黑色实线和虚线分别为根据散点拟合出的剖面线；(d) AA′，BB′和GG′剖面值，剖面线位置见图2a，背景为经过剖面AA′的SRTM高程值；(e) CC′和DD′剖面值，背景为经过剖面BB′的SRTM高程值；(f) EE′和FF′剖面值，背景为经过剖面GG′的SRTM高程值. Fig. 2 The coseismic deformation map of LoS from D-InSAR Each fringe represents 2.8 cm for Sentinel-1 A and 11.8 cm for AL0S2. (a) Interferogram of track 27 for Sentinel-1 A; (b) Interferogram of track 107 for Sentinel-1 A; (c) Interferogram of trackl62 for ALOS2. Line-of-Sight (LoS) displacements along the profiles shown above. Dots are InSAR LoS changes extracted from the interferograms along the profiles with a buffer of 1 km; (d), (e), and (f) correspond to profiles AA′-BB′-GG′, CC′-DD′, and EE′-FF′respectively.

由于视线向模糊问题，InSAR所获取的形变是地表三维形变在视线向(LoS)上的投影，多数情况下这种一维形变不能完全反映地表的真实形变.近年来国内外学者在三维形变场提取方面开展了一些研究工作(Wright et al., 2004; Fialko et al., 2005; González et al., 2009; Hu et al., 2010, 2012; 洪顺英等，2013; 温少妍等，2016)，方法日趋成熟.而本次阿克陶地震获取了不同卫星平台的不同视角的三幅干涉像对，使得提取阿克陶三维形变场成为可能.本文采用王家庆等(2016)提出的一种多视角D-lnSAR间接解算法，即结合先验知识的最小二乘迭代逼近法，可在有效降低形变解算误差的同时，重建更为真实可靠的三维形变场.

根据SAR影像的成像几何关系，可以得到形变量图上任意一点的视线向形变关系式.用d_E、d_N、d_U分别表示地表位移在东西、南北、上下三个方向(E、N、U)上的分量.卫星视线向位移d_los是三个位移矢量在卫星视线向的投影(Hanssen，2001)，对右视的雷达观测数据分别有：

(1)

式中θ_inc、α_azi分别代表入射角和卫星飞行方位角.从(1)式可以看出，如果知道某一点的三个视线向形变观测值，即可求出该点的位移三分量.对于阿克陶地震，我们已经获取了Sentinel升、降轨及ALOS2升轨的三个视线向形变场，则可构建附有约束条件的三维分量矩阵方程为：

(2)

(3)

d_n×1为n组视线向观测值组成的矩阵，F_n×3为投影系数矩阵，由SAR成像几何条件确定；W_n×1为约束条件；B为条件系数阵.本文解算方法则是在式(3)的约束下，最小二乘迭代解算式(2).依据三维分量对LoS向形变的贡献程度以及解算精度的高低，从解算结果中选择可信度最高的分量作为约束条件，再用最小二乘解算出其他分量，依次得到三维形变分量.根据雷达对南北向形变不敏感的特征，在解算过程中，将南北向形变向极小逼近作为先验的初始约束条件，依据三维分量对视线向形变的贡献程度以及解算精度的高低，从解算结果中选择可信度最高的分量作为约束条件，再用最小二乘解算出其他分量；依次得到三维形变分量.因没有实际观测资料(如GPS、水准等)进行对比，因此很难直接对解算后的三维形变场可靠性进行定量评估.可利用误差估计公式计算直接解算模型与间接解算模型的误差放大倍数说明该方法在可靠性方面的改进.文中使用的三个不同模式下观测得到的形变场是相对独立的, 因此我们假设不同视线向形变场之间相互独立且具有相同的观测精度σ＝1，利用公式(4)和(5)计算得到各解算模型中地表三维分量的解算误差放大效应.

(4)

(5)

式中Q_d和Q_ls分别为直接解算和间接解算模型三维形变分量的协方差阵，其中F为数据的中心入射角和方位角得到投影系数矩阵，N_PP= F^TpF，N_BB= BN_PP^-1B^T，B为条件系数矩阵，p为观测数据的权重矩阵.将观测数据的方位角、入射角等参数代入(5)式, 即可量化估计解算误差的放大效应.

从表 3结果可以看出，两种解算方案中南北向分量的形变中误差都大于其他两个分量，尤其是在直接解算方案中，中误差被病态系数矩阵放大至48倍，而在间接解算方案中，则能有效抑制误差放大，提高可靠性.图 3即为解算结果.垂向分量可以看出(图 3a)，断层北盘呈隆升形态，形变量约10 cm，而断层南盘则存在两个沉降中心，尤其震中东部的沉降量值最大，可达20 cm，正断作用明显.东西向分量显示(图 3b)，北盘向东运动，量值约10 cm，南盘也分两个运动中心向西做拉张运动，最大量值约8 cm，表现出右旋走滑性质.南北向分量相对较小(图 3c)，且只有南盘有不到4 cm的北向运动，而北盘的南北向运动趋势不明显.从直观表现来看，分解后的结果能够更清晰地描述地震的形变特征，整体表现为右旋兼具正断作用，在震中东部分正断作用更加明显.

表 3 Table 3

表 3 不同解算方法下三维分量的中误差估计 Table 3 RMSE of full vector components resolved by different method σU σN σE 直接解算 4.9 48.3 2.1 间接解算 0.7 5.6 1.2

表 3 不同解算方法下三维分量的中误差估计 Table 3 RMSE of full vector components resolved by different method

图 3

Fig. 3

图 3 解算的三维形变场 (a)垂直向; (b)东西向; (c)南北向. Fig. 3 Three-dimension deformation maps decomposed from three pairs of interferograms Vertical (a), horizontal east (b) and horizontal north (c) components of the surface displacement field.

为了获取地震震源参数，通常是采用两步法策略，首先是基于矩形位错模型采用非线性反演方法确定断层的几何参数；然后通过线性反演方法来估计断层面上的精细滑动分布.本次地震从形变场的形态与地表已知断裂的分布情况，即可确定发震断裂的走向及倾向，结合震源机制解可以推断出破裂面.故本文先用“前向模拟”方式，根据震中位置、断层面尺寸、滑移矢量等参数来模拟干涉纹图，然后与实际观测形变场比较，通过对断层参数不断调整，最终找到最佳的参数组合，再将其作为线性反演时的参数输入估算断层面上的精细滑动分布.

假定断层模型可用弹性半空间中的一个有限断层面来表示(Okada，1985)，参考CENC的震源机制解，确定断层走向107°，倾角76°，滑动角-168°，滑动量0.8 m，上断层埋深3 km，下断层埋深17 km，断层出露地表中心参考点坐标(74.19°E，39.25°N)，断层长度25 km，计算有限矩形源在弹性半空间状态下位错引起的地表形变量.由该均匀滑动分布模型得到的地震矩为9.32×10¹⁸ N·m，换算为矩震级约M_W6.6.图 4为根据以上参数模拟得到的三分量形变场.东西向分量占了很大比重，南盘的西向运动约10 cm(图 4b)，在垂向上(图 4a)，虽然也有两个形变中心，但震中西部却呈隆升状态，东部的形变中心虽是沉降，但离震中位置太近.图 4c显示在断层的东部尾端有6 cm左右的北向挤压运动.我们期望的双沉降模式没有出现.

图 4

Fig. 4

图 4 单断层模型模拟的三维形变场 垂向、东西向、南北向运动分量分别如图(a), (b), (c)所示. Fig. 4 The 3-D displacement field for the Aketao earthquake inverting from single fault model Up, east, and north components are shown in (a), (b), (c) respectively.

图 5为单断层均匀滑动分布模型模拟得到的升、降轨缠绕LoS向形变场及残差.升轨模拟形变场在断层南盘震中东西两部分均有较大残差，可达两个条纹(图 5c)，而降轨模拟形变场在近断层附近残差更大，并且在真实形变场中观测到的两只牛眼状形变椭圆(图 2b)在这里并未出现，只有几乎对称的蝴蝶状条纹分布在断层两侧.而这些差别说明所观测到的变形模式并非仅仅是单一断层上的近乎纯走滑运动.这些差异需要额外的向下的运动量进行弥补，可通过改变滑动角的大小来增加正断的作用量.

图 5

Fig. 5

图 5 单断层模型(走向107°，倾角76°，滑动角-168°，滑动量0.8 m，断层长度25 km，上、下断点埋深为3 km和17 km)均匀滑动分布模拟的Sentinel-1A升、降轨LoS向形变场 (a)升轨; (b)降轨; (c)升轨残差; (d)降轨残差. Fig. 5 Single fault models of the Aketao earthquake. Calculated LoS deformation of an Okada elastic dislocation model with assumed fault parameters (strike 107°, dip 76°, rake -168°, slip 0.8 m, length 25 km, top depth 3 km, bottom depth 17 km) (a) For ascending, (b) for descending and (c), (d) are their residual respectively.

为求证这一假设，我们通过求解扩展断层面上的非均匀分布滑动来验证.采用汪荣江的Steepst Descent Method(SDM程序)进行断层面上精细滑动分布反演(Wang et al., 2013).该程序中的模型按照Aki传统的断层约定，确定了起点、终点位置后即可确定断层走向.按照木吉断裂的地表走向，断层模型的起始参考点设为(73.798°E，39.328°N)，断层面沿走向扩展到66 km长，总体走向107°，即可确定断层的地表迹线(图 6中白色迹线).滑动角允许在-180°~180°之间变化，断层面宽度为20 km，延伸到地表，并划分成1.0 km×1.0 km的离散网格，共计1320个，形成滑动分布断层模型.该单断层非均匀滑动分布模型模拟的形变场及残差如图 6所示，详细参数见表 4，断层面上的滑动分布结果如图 7所示.滑动分布分为两部分，震中东部的滑动量略大一些，最大滑动量可达0.8 m，深度10 km左右.总体上为右旋走滑兼少部分正断运动.在均匀滑动模型中期望的向下部分的正断作用在该模型中有所体现，总体模拟效果优于均匀滑动模型，但近断层部分仍有较大的残差存在，尤其对于降轨数据，最大仍有约4个条纹的残差存在，说明单断层模型并不能很好地解释观测数据.

图 6

Fig. 6

图 6 单断层模型非均匀滑动分布模拟的Sentinel-1A升、降轨LoS向形变场 (a)升轨; (b)降轨; (c)升轨残差; (d)降轨残差.白色线段代表断层在地表出露位置. Fig. 6 Single fault model with variable slip and rake for the Aketao earthquake (a) Model ascending interferogram; (b) Model descending interferogram; (c) Residual ascending interferogram; (d) Residual descending interferogram.White lines show location of the fault when projected up dip to the surface.

表 4 Table 4

表 4 基于弹性位错模型模拟结果 Table 4 Source parameters of the Aketao earthquake from elastic dislocation modeling 模型 分段 震中/(°) 走向/倾角/滑动角/(°) 滑动量/m 长度/km 宽度/km 矩震级 深度/km 均匀 1 74.19 39.23 107/76/-168 0.8 25 15 MW6.6 17.0 非均匀 1 74.37 39.20 107/76/-165 0.8 66 20 MW6.6 8.5 均匀 2 74.10 39.25 107/76/-168 0.8 13 15 MW6.4 17 74.36 39.20 107/76/174 0.8 12 15 MW6.4 17 非均匀 2 74.10 39.23 103/76/-167 0.6 38 20 MW6.5 10.1 74.37 39.16 107/54/-160 0.7 34 20 MW6.5 5.3

表 4 基于弹性位错模型模拟结果 Table 4 Source parameters of the Aketao earthquake from elastic dislocation modeling

图 7

Fig. 7

图 7 非均匀滑动分布模型断层面上的滑动分布 黑色箭头指示上盘断层面相对下盘运动方向. Fig. 7 Distribution of slip and rake on the fault for this model Black arrows indicate direction of motion of the hanging wall of the fault with respect to the foot wall.

如果地表观测到的形变模式是由两段断层破裂共同产生的，那么它们产生的形变是混叠在一起的，我们在模拟的时候也是无法将信号分离开来的.并且面对成倍增长的几何参数，如果采用非线性反演方法确定断层的几何参数，不仅计算量巨大而且结果可能也无法收敛.故仍采用跟单断层模型相同的步骤，假定两段断层的贡献是相同的，根据震中位置、断层面尺寸、滑移矢量等参数采用正演的方式来模拟干涉纹图，然后与实际观测形变场比较，找到最佳的参数组合，再通过线性反演方法来估计断层面上的精细滑动分布.

综合参考USGS和CENC的震源机制解，将断层1参数定为走向107°，倾角76°，滑动角-168°，断层出露地表中心参考点坐标(74.10°E，39.25°N)，断层长度13 km；断层2走向和倾角不变，滑动角设为174°，断层出露地表中心参考点坐标(74.36°E，39.20°N)，断层长度13 km；其余参数一致，即滑动量0.8 m，上断层埋深3 km，下断层埋深17 km，计算有限矩形源在弹性半空间状态下位错引起的地表形变量.两段断层模型计算得到的地震距分别为4.847×10¹⁸ N·m和4.474×10¹⁸ N·m，换算为矩震级约M_W6.4.图 8为根据以上参数模拟得到的三分量形变场.东西向分量占优势，南盘的西向运动约12 cm(图 8b)，在垂向上出现了两个沉降中心，但二者中间的震中附近出现了一个隆升形变区域(图 8a)，图 8c显示在断层的两端均有约6 cm左右的北向挤压运动.

图 8

Fig. 8

图 8 双断层均匀滑动模型模拟的三维形变场 垂向、东西向、南北向运动分量分别如图(a), (b), (c)所示. Fig. 8 The 3-D displacement field for the Aketao earthquake inverting from two segment uniform slip model Up, east, and north components are shown in (a), (b), (c) respectively.

图 9为双断层均匀滑动分布模型模拟得到的升、降轨缠绕LoS向形变场及残差.模拟得到的形变场在形态上更加接近真实观测，升、降轨均出现了两个椭圆形变中心，残差值也有所下降，降轨模拟形变场在断层东部分最大残差为三个条纹(图 5c).这可能是由于对于两段断层我们设置了相同的上、下断层埋深，而根据实际观测形变场的形态(图 2)推断，位于东侧的部分的埋深应该浅于西侧，因此模拟得到的条纹没有实际观测的那么密集.

图 9

Fig. 9

图 9 双断层模型均匀滑动分布模拟的Sentinel-1A升、降轨LoS向形变场 (a)升轨; (b)降轨; (c)升轨残差; (d)降轨残差.白色线段代表断层在地表出露位置. Fig. 9 Best fitting two fault uniform slip model of the Aketao earthquake (a) Model ascending interferogram; (b) Model descending interferogram; (c) Residual ascending interferogram; (d) Residual descending interferogram. White lines show locations of the faults when projected up dip to the surface.

为求证这一假设，我们再次通过求解扩展断层面上的非均匀分布滑动来验证.按照木吉断裂的地表走向，并参考余震分布在震中附近有一转折，模型两段断层的走向略有差异，断层1起始参考点(73.817°E，39.307°N)，断层面沿走向扩展到39 km长，总体走向103°；断层2起始参考点(74.273°E，39.227°N)，断层面沿走向扩展到34 km长，总体走向109°，即可确定断层的地表迹线(图 11中白色实线).滑动角允许在-180°~180°之间变化，断层面宽度为20 km，延伸到地表，并划分成1.0 km×1.0 km的离散网格，共计1460个，形成非均匀滑动分布断层模型.关于断层倾角的确定，采用了网格搜索的方法.首先认为两段断层顶部和底部倾角是相同的，在40°~90°范围内变动，步长为5°，得到一个均方差为最小时的倾角值；将西部的断层1倾角固定为该值，让东部的断层2顶部和底部倾角都是可变的，范围在50°~90°之间，步长为5°，得到一个均方差全局最小时的东部断层2的倾角取值，再将该值固定，重复同样的步骤即可确定整个断层的倾角值，最终西段断层倾角75°，东段倾角55°.

图 11

Fig. 11

图 11 双断层模型非均匀滑动分布模拟的Sentinel-1A升、降轨LoS向形变场 (a)升轨; (b)降轨; (c)升轨残差; (d)降轨残差.白色线段代表断层在地表出露位置. Fig. 11 Two fault variable slip model of the Aketao earthquake (a) Model ascending interferogram; (b) Model descending interferogram; (c) Residual ascending interferogram; (d) Residual descending interferogram.White lines show locations of the faults when projected up dip to the surface.

图 10为据此参数模拟的形变场三分量.东西向分量占据主导地位，木吉断裂南侧朝西运动，反映了断层的运动性质为右旋走滑.同时，在垂直向分量上这次仅出现了两个沉降中心，在均匀滑动模型(图 8a)震中部分出现的隆升已经消失，跟实际观测值分解得到的三分量形变场更加接近，说明该断层模型更加接近真实情况.图 11为双断层非均匀滑动分布模型模拟得到的升、降轨缠绕LoS向形变场及残差.依然是升轨数据拟合得较好，均方差在1.3 cm，而降轨数据拟合的均方差在2.2 cm，在近断层附近残差稍大，可能是由于采用的断层模型较为简单，跟实际的弯曲断层在位置上有所偏差.断层面上的滑动分布见图 12，模型参数见表 4.两个子断层面上各有一个滑动分布集中区，其中东部子断层面上滑动分布峰值为0.67 m，深度为5.3 km，而西部子断层面滑动分布峰值为0.61 m，深度为10.2 km.该滑动分布跟地表观测形变场的形态具有很好的对应关系.

图 10

Fig. 10

图 10 双断层非均匀滑动模型模拟的三维形变场 垂向、东西向、南北向运动分量分别如图(a), (b), (c)所示. Fig. 10 The 3-D displacement field for the Aketao earthquake inverting from two segment nonuniform slip model Up, east, and north components are shown in (a), (b), (c) respectively.

图 12

Fig. 12

图 12 非均匀滑动分布模型分别在东、西两个断层面上的滑动分布 黑色箭头指示上盘断层面相对下盘运动方向. Fig. 12 Distribution of slip on the west fault and distribution of slip on the east fault Black arrows indicate the direction of motion of the hanging wall of the fault with respect to the foot wall.

本次地震发生在帕米尔构造结东北部，该地区晚第四纪以来构造变形非常强烈.印度板块与欧亚大陆之间的汇聚作用是区域构造作用强烈的根本原因，西构造结向北的楔形挤入作用造成帕米尔及其前缘坳陷内的弧形构造，刚性的塔里木盆地基底也控制了变形的扩展, 将帕米尔的作用力传递到了盆地北缘与天山之间, 造成那里的强烈变形(沈军等，2001).在晚新生代，帕米尔的变形以前缘正向逆冲、两侧斜冲-走滑和内部拉张为特征.

而帕米尔内部的构造变形以西部卡拉库尔地堑和东部250 km长的公格尔拉张系近东西向拉张作用为主(Brunel et al., 1994)，晚新生代以来拉张系的东西向拉张量总体北大南小，在最北端的木吉盆地约为30 km，在慕士塔格山约为20 km，在南端的塔什库尔干盆地则小于3 km(Robinson et al., 2007).这一观察结果与帕米尔东西延伸由右旋走滑喀喇昆仑断层向北传播驱动的模型不一致，表明这种拉张作用是由于沿主帕米尔主断裂(MPT)的径向推力，或整个帕米尔地区的复合弯作用下的地形塌陷导致(Robinson et al., 2007).本次地震发震断层木吉断层就位于公格尔拉张系的北部，该断层是一条近NW走向以走滑运动为主的断裂，在卫星影像上显示为一条非常显著的线性构造及错断标志，同时在木吉盆地内直而陡峭的断裂面构造又显示出该断裂为正断性质(刘栋梁等，2011).木吉断层、公格尔山断层之间由小的阶区和转折相连，组成公格尔拉张系的北段，木吉断层右旋走滑兼具正断作用是公格尔拉张系形成应力场环境下的产物.

结合前人对帕米尔构造演化的研究及本次地震余震分布、InSAR反演结果等提出本次地震的发震构造模型如图 13所示.

图 13

Fig. 13

图 13 2016阿克陶地震发震构造示意图 (a)公格尔拉张系北段演化示意图，修改自Robinson等(2007), 图中虚线框为(b)图范围，且逆时针旋转90°; (b)发震构造三维示意图及沿倾向断层滑动分布，黄色线条为断层模型地表出露迹线，红色线条为断层分布. Fig. 13 Cartoon seismotectonics map for 2016 Aketao earthquake (a) Kinematic model for the development of the northern part of Kongur Shan extensional system modified from Robinson et al. (2007). The dashed frame is the range of Fig. 13b and rotated 90 degrees counterclockwise. (b) Landsat8 image overlayed on topography with fault trace at surface marked by the yellow line and distribution of slip for the two fault planes.

InSAR观测显示本次地震沿木吉地震形成约70 km左右长的形变场，且明显以震中为界分为东、西两部分，西部的形变条纹稀疏而东部密集，似乎是由两次地震事件叠加而成的.通过网格搜索断层倾角，在对比测试60个不同倾角组合的断层模型后，确定了两段式的非均匀滑动断层模型，断层在地表的出露迹线如图 13b所示.综合对比反演的均方差与模型的相关度，最后确定的两段式非均匀滑动模型反演结果显示，滑动分布有两个集中分布区，西段反演得到的地震矩7.2×10¹⁸ N·m，东段反演得到的地震矩为7.76×10¹⁸ N·m，前后有两次M_W6.5左右的地震叠加而成的形变场，东段的滑动值略大于西段.西段主要以走滑为主，而东段是走滑兼具正断作用.西段先右旋走滑运动，继而牵引带动东段走滑下坠运动.

在西段断层倾角较陡，并且以右旋走滑为主，而在该次右旋走滑事件的牵引作用下，东段也开始破裂，东段断层倾角较缓，在地形重力作用下断层上盘又开始塌陷，进而表现为正断作用.同时在帕米尔高原的总体北向南碰撞作用下，对断层的北盘又有挤压作用，断层北盘形成隆起.而发生于2008年3月12日的于田地震跟本次地震有着类似的构造几何类型.那次地震是在近南北走向、西倾的节面上发生的正断裂兼有左旋走滑分量的破裂模式.根据Furuya和Yasuda(2011)的研究结果，2008年M_S7.3于田地震可能存在着2个破裂面，即NW向的拉张型主破裂面和NE向的左旋走滑次破裂面，该模型较Elliott等(2010)的结果拟合效果更好.而发生在帕米尔前缘推覆构造带的1985年8月23日乌恰7.4级地震震源情况更为复杂，地表破裂表明发震构造可能为一弧形构造.主震可能由三次事件组成，一个以逆冲为主的事件发生在弧凸所在处；另两个以右旋走滑的事件分别发生在弧缘两侧，且走向、倾角等不一(Burtman and Molnar, 1993; 陈杰等，1997).这些震例表明，该地区的构造变形复杂，而本次地震极有可能是由两次地震事件组成的，一次以右旋走滑为主，而另一次走滑则带有很大的正断作用.

本文通过对比试验，测试不同断层模型反演的可靠性，并综合前人对帕米尔构造结及邻区的晚新生代构造与现今变形研究成果、本次地震震源机制解和已有研究，给出了2016年阿克陶地震的发震构造模型，得到了以下几点认识：

(1) 利用升、降轨Sentinel-1A，ALOS2数据得到了2016年阿克陶地震的LoS向同震形变场，形变场条纹主要分布在木吉断裂的南侧，位置与木吉断裂高度吻合，破裂长度约70 km，在震中东西两侧各形成一个形变中心，且西部条纹稀疏而东部密集，横跨断层的形变剖面显示LoS向的形变值最大约12 cm.采用一种结合先验知识的多视角D-lnSAR最小二乘迭代逼近解算法，求解出了三维形变场.结果显示垂直分量和东西向分量较大，断层运动状态为右旋走滑兼具正断作用.断层北盘隆升量约10 cm，而断层南盘则存在两个沉降中心，尤其震中东部的沉降量可达20 cm，正断作用明显.北盘东向运动量值约10 cm，南盘也分两个运动中心向西做拉张运动，最大量值约8 cm，表现出右旋走滑性质.南北向分量相对较小，南盘有不到4 cm的北向运动.

(2) 本次地震由InSAR形变场基本可以确定断层的地表走向等参数，故只有断层倾角、起点位置参数未知，先采用前向模拟的方式对比观测形变场确定发震断层参数，再将其作为线性反演时的参数输入估算断层面上的精细滑动分布.反演过程中最难确定的是断层倾角，为此设计了等间距的倾角序列，使用网格搜索法来寻找最佳的拟合参数，最终确定了两段式的非均匀滑动断层模型，基于弹性半空间形变模型来反演断层面上的精细滑动值.反演结果显示本次地震可能为两次地震事件叠加而成的，两次事件分别发生在东、西两段断层上.西段主要以走滑为主，倾角76°，东段的滑动值略大于西段，破裂峰值约0.8 m.而东段倾角较缓约55°，走滑兼具正断作用.

(3) 本次阿克陶地震是由印度板块与欧亚大陆之间的碰撞，西构造结向北的楔形挤入作用造成帕米尔及其前缘坳陷内的弧形构造，而发生在公格尔拉张系的一次构造地震，是该拉张系统自1895年塔什库尔干7级地震以来发生的最大一次地震，在其北部不远的东西走向的弧形构造上发生过多次地震，2008年10月5日乌恰M_W6.7地震距离此次震中最近，而时间最近的一次中强震是2016年6月26日吉尔吉斯斯坦M_W6.4地震，这可能表明该区的地震活动日趋活跃，未来仍有发生强震的可能性.