印度板块与欧亚板块碰撞俯冲造成青藏高原隆升的同时，在其边界断裂带上形成一系列叠瓦状逆冲断层，由南向北依次展布主前缘逆冲断裂(MFT)、主边界逆冲断裂(MBT)以及主中央逆冲断裂(MCT)，并在深部汇聚于倾角平缓的主喜马拉雅逆冲断裂(MHT)(Lavé and Avouac，2000；Cattin and Avouac, 2000；Bettinelli et al., 2006；王卫民等，2015).这些断裂系构成了此地区强震的构造背景，历史上曾发生多次8级及以上强震.2015年4月25日06时11分26秒(UTC时间)，尼泊尔境内廓尔喀地区发生了M_W7.8地震，位于1505年M_W8.2地震与1934年M_W8.1地震之间的空区(图 1)，随后两小时内又发生了两次M_W6级以上的余震，最大余震发生在5月12日，矩震级达M_W7.3级.2015年尼泊尔地震共造成将近9000人丧生，20000以上人员受伤(Bilham，2015；Wei et al., 2018).震间形变研究(Avouac et al., 2015；Bollinger et al., 2016；Wei et al., 2018)结果表明此区域仍处于强闭锁状态，累积滑动亏损大于10m，证明震间积累的应变能在同震阶段并未完全释放，该地区仍有很强的地震危险性.

地震发生后，已有大量的文章基于大地测量数据、地震数据单独(Wang and Fialko, 2015；Feng et al., 2016；Wu et al., 2016)或联合(Galetzka et al., 2015；Avouac et al., 2015；Wei et al., 2018)研究了该地震的震源特征、破裂过程等.

综观这些文献，在同震滑动分布反演中，为了获得更为精细的断层滑动模型，需要将断层面上的非均匀滑动在整个破裂面上进行分段或区域平均，通常基于矩形位错理论将断层面剖分成相同大小的子断层(Wu et al., 2016；Feng et al., 2016；Yue et al., 2016；Zha and Dai, 2017)，或者子断层尺寸随深度而增加(Wang and Fialko, 2015)，但是按照这两种方法得到的滑动分布会表现出一些基于观测数据不能完全确定的滑动特征(伪滑动)，在断层面的某些位置呈现出太多的细节，或者滑动分布结果过于平滑以至于无法得到更为精细的特征(温扬茂等，2012).并且由于矩形位错元的几何性质，如果滑动区域具有弯曲的边界，那么将其划分为多个矩形单元来模拟时，断层面会出现空隙或重叠的情况.

其次在目前发布的该地震同震滑动分布模型各有差异，滑动模式从沿着走向扩展(Lindsey et al., 2015)，到沿着走向和倾向两个方向分布(Yagi and Okuwaki, 2015)，最大滑动量从4.9 m(Yadav et al., 2017)到8 m(Elliott et al., 2016；McNamara et al., 2017).这些差异可能是源于采用了不同的反演算法，观测数据源，断层面几何结构及剖分方式，但作者发现，在采用多源数据联合反演时，许多文献(Grandin et al., 2015；谭凯等，2016；Zhang et al., 2016, 2017；McNamara et al., 2017；Wei et al., 2018)仅是简单地根据经验或观测值的验前方差进行定权，实践证明，许多情况下这种定权方式是不够精确的(崔希璋等，2009)，它会影响反演结果的可靠性.

针对上述问题，本文基于三角位错的断层自动剖分方法(Barnhart and Lohman, 2010)来划分断层单元，同时通过联合GPS和InSAR数据，采用赫尔默特方差分量估计的方法(崔希璋等，2009；Xu et al., 2009, 2010；许才军等，2016)，确定GPS和InSAR数据的最佳权比，来建立2015年尼泊尔M_W7.8地震的同震滑动分布模型，并与已有研究结果进行对比分析.

2015年尼泊尔M_W7.8地震发生后，我们收集了来自JAXA(Japan Aerospace eXploration Agency)的降轨影像，即ALOS-2(Advanced Land Observing Satellite-2)卫星PALSAR-2雷达(Phased Array type L-band Synthetic Aperture Radar-2)的ScanSAR(Scanning Synthetic Aperture Radar)影像.通过使用双通道的差分InSAR方法和Caltech/JPL的ROI_PAC(Repeat Orbit Interferometry Package)软件(Rosen et al., 2004)，从PALSAR-2格式1.1级的SLC数据得到干涉影像.采用由JAXA提供的精密轨道产品对干涉图进行初始处理，之后通过使用ScanSAR的干涉处理器将五个子条带的干涉图对齐(Liang et al., 2013).并通过使用3弧秒(约90 m)精度的SRTM(Shuttle Radar Topography Mission)的DEM(Digital Elevation Model)模型，移除干涉图中的地形效应(Farr et al., 2007).然后使用自适应功率谱滤波器(Goldstein and Werner, 1998)对展平的干涉图进行滤波，并使用枝切法(Goldstein et al., 1988)对干涉图进行解缠.采用四叉树和线性函数方法移除干涉图中的残余轨道误差和地形相关的大气延迟.为了减小数据量以及它们之间的相关性，我们使用四叉树采样方法(Jónsson et al., 2002)对InSAR数据进行降采样.最后从T048A干涉图中提取了1180个LOS(Line-Of-Sight)位移点，用于本文反演的数据集.

为了尽可能多地包含GPS数据，我们在本文中采用来自Yadav等(2017)的50个GPS站的静态同震位移数据(GPS站点分布如图 1)，其中尼泊尔境内13个站点来自Galetzka等(2015)，中国境内33个站来自Wu等(2016)，用于本文的联合反演中.

图 1

Fig. 1

图 1 2015年尼泊尔MW7.8地震构造背景 红色五角星及沙滩球表示MW7.8级主震及MW7.3最大余震的震中及震源机制(来自USGS).浅红色填充的椭圆表示1505年MW8.2地震和1934年MW8.1地震破裂区域.红色实线代表研究区的断层迹线(Deng et al., 2003).红色锯齿线表示主前缘逆冲断裂(MFT).黑色沙滩球代表 1976—2015年5级以上地震分布(来自GCMT).绿色三角形表示本文选用的GPS站点分布.黑色虚线框表示ALOS-2 InSAR影像的空间覆盖范围. Fig. 1 Tectonic setting of the Nepal region The red stars and beach balls represent the epicenter and focal mechanism of the MW7.8 mainshock and the MW7.3 aftershock (from USGS). The two light red filled ellipses represent the 1505 MW8.2 and the 1934 MW8.1 earthquake rupture zone. The red solid lines represent the fault traces of the study area (from Deng et al., 2003). The red sawtooth line represents the Main Frontal thrust (MFT). The black beach balls represent the spatial distribution of earthquakes of magnitude ≥5 from 1976—2015 (from GCMT). The green triangles indicate the GPS sites used in this paper. The black dashed box indicates the spatial coverage of the ALOS-2 InSAR image.

图 2显示了50个GPS站点的同震水平位移，符合逆冲断层破裂的形变特征.在M_W7.8级主震以东、加德满都盆地以北的区域发生了较大的同震位移.KKN4站经历了最大水平位移1.83 m，显著地向南运动，并伴随着1.26 m的抬升.尼泊尔境内NAST、CHLM、KIRT站发生了1.3~1.5 m的水平位移.中国藏南地区的GPS站点表现为显著的向南运动，其中最大位移0.48 m发生在J041站.距离震中400 km以外的站点观测到微弱的位移(< 5 mm).

图 2

Fig. 2

图 2 基于方差分量估计定权和断层自动剖分的2015年尼泊尔MW7.8地震同震滑动分布以及GPS同震位移观测值与模型位移值，其中橙黄色箭头表示按验前方差定权的模拟位移值 Fig. 2 The coseismic slip distribution of the 2015 Gorkha (Nepal) earthquake based on HVCE-FR method, with the colored arrows indicating the GPS observations and synthetics, respectively. The orange arrows represent the GPS synthetics based on WPV-FR method

我们利用上述InSAR和GPS数据，采用Okada均匀弹性半空间位错模型，基于多峰值粒子群(Multiple peak Particle Swarm Optimization，MPSO)算法(冯万鹏和李振洪，2010)，运用PSOKINV程序包(Feng et al., 2013)，通过使观测值与模型值之间有最佳拟合，来估计断层的9个参数(断层位置，深度，长度、宽度、走向、倾角、滑动角和滑动量).考虑到InSAR和GPS数据之间的精度差异，本文采用赫尔默特方差分量估计(HVCE)的方法(崔希璋等，2009；Xu et al., 2009, 2010；许才军等，2016)来确定InSAR数据和GPS水平、垂直分量的权比.首先根据观测值的中误差进行初始权比估计，经过迭代计算，利用观测值的验后残差进行方差分量估计，求得InSAR、GPS水平向与垂直向观测值的单位权方差，并依据求得的单位权方差重新定权，如此反复进行，直到此三类观测值的单位权方差之比近似等于1，终止迭代计算(崔希璋等，2009；Xu et al., 2009, 2010；许才军等，2016).此时得到InSAR数据和GPS水平、垂直分量的相对权比为1:14.5:3.0，这既顾及GPS数据能够提供更高精度的形变信息的优点，又顾及InSAR数据能够提供更多观测量，更能反映大范围形变的整体特点(王乐洋等，2017；Wang et al., 2018a, b)，本文GPS观测数据精度更高，相应权比更大，对最终参数估计的贡献程度更大(王乐洋等，2012).

与Feng等(2016)类似，我们的结果表明，走向角与滑动角之间有很强的折中关系(图 3)，且GPS、InSAR拟合残差RMS值在相差1 cm范围以内波动，表明本文采用的数据不利于唯一确定走向角，因此我们将走向角固定为290°(Hayes et al., 2015; Feng et al., 2016; McNamara et al., 2017; Yang et al., 2018)进行反演，得到结果如表 1所示.为了评估断层参数的不确定性，本文采用蒙特卡洛方法，通过对观测数据产生100组带随机扰动误差的数据集，以此进行非线性反演，从而估计断层参数的精度(表 1).我们的反演结果显示此次地震发震断层以低倾角逆冲为主，兼少量右旋走滑(滑动角100.5°)，与已发表的结果有比较好的一致性.其中反演得到的断层位置与Feng等(2016)结果相近，断层尺寸与谭凯等(2016)、Yang等(2018)结果比较接近，倾角接近于GCMT、Feng等(2016)的结果，滑动角与USGS公布的结果吻合.

图 3

Fig. 3

图 3 均匀滑动反演的走向角与滑动角之间的折中关系 Fig. 3 The tradeoff between strike and rake angles for the 2015 Gorkha (Nepal) earthquake

表 1 Table 1

表 1 不同研究提供的发震断层参数 Table 1 The various source parameters for the 2015 Gorkha (Nepal) earthquake 位置 长度(km) 宽度(km) 震源机制 数据 来源 经度(°) 纬度(°) 深度(km) 走向(°) 倾角(°) 滑动角(°) 84.731(1) 28.231(1) 15(1) 290 7 101 远震体波 USGS 85.33(2) 27.91(2) 12(2) 287 6 96 远震体波 GCMT 85.9265(3) 27.5866(3) 14.9(3) 290 7 99 ALOS-2, RADARSAT-2 InSAR Hayes et al., 2015 85.392(4) 27.917(4) 8.5(4) 115 42 284.7 14 92.5 ALOS-2, Sentinel-1A, RADARSAT-2, GPS Feng et al., 2016 85.392(4) 27.917(4) 9.1(4) 180 92 290(fixed) 6 102.3 85.43(4) 27.93(4) 15(4) 140 60 291 7.6 95 ALOS-2/Sentinal-1A InSAR 李永生等, 2016 85.2280(5) 27.9688(5) 13.8(5) 80 44 285 9 - GPS, ALOS-2 InSAR 谭凯等, 2016 85.880±0.04(3) 27.583±0.02(3) 9.0±0.14(3) 114.9 40.1 288.8 6.6 114.3 GPS, ALOS-2 InSAR Yang et al., 2018 85.782±0.03(3) 27.596±0.03(3) 9.4±0.19(3) 94.5 41.6 290.0 9.3 116.8 85.402±0.004(4) 27.964±0.001(4) 9.2±0.2(4) 85.8±1.0 59.8±0.3 290.0fixed 5.9±0.3 100.5±0.3 GPS, ALOS-2 InSAR 本文a 85.401±0.004(4) 27.967±0.001(4) 9.3±0.2(4) 83.6±1.1 60.0±0.3 290.0fixed 5.0±0.3 100.0±0.4 本文b 注:(1)初始破裂点; (2)地震矩质心; (3)断层面左上角位置; (4)断层平面中心位置; (5)最大破裂点; a方差分量估计定权的非线性反演结果; b按数据中误差定权的非线性反演结果.

表 1 不同研究提供的发震断层参数 Table 1 The various source parameters for the 2015 Gorkha (Nepal) earthquake

表 1展示了方差分量估计与数据中误差定权下反演的断层参数差异，从本文结果来看，两者有些许差别，主要表现在后者反演得到的倾角更小，小于USGS发布的结果，而两者反演的其他参数则比较接近.从给出的InSAR、GPS水平向与垂直向拟合残差RMS来看，前者结果分别为0.095 m、0.044 m、0.072 m，优于后者结果0.103 m、0.052 m、0.077 m，证实方差分量估计的定权方法取得了相对更优的结果.与单纯按观测数据中误差定权相比，方差分量估计定权是更合理的，其参数估计具有无偏性(崔希璋等，2009)，从而获得的震源机制更加稳健.Ragon等(2018)分析了断层几何结构的不确定性在地震震源反演问题中的影响，表明断层几何结构的不确定性可能使得浅部滑动估计产生偏差，断层结构中假定倾角的不准确性可能是观测到的大陆地震浅部滑动亏损的另一个原因，并且大型逆冲事件中对破裂断层形态约束不足也可能对近地表滑动估计和海啸灾害评估产生影响.因此，精确地估计断层几何参数对研究震源破裂特征和震后灾害快速预警具有重要的意义.

由于基于均一断层模型的假设对于分布式滑动模型并不是最优的，更为精细的滑动模型需要采用分布式断层模型.在分布式滑动反演过程中，进一步对倾角δ与平滑因子α进行优化约束.采用不同的倾角值，通过调整不同的平滑因子来分析模型粗糙度和残差的变化.log函数模型f(δ, α)=log(ψ+ξ)可用于获得反演的最优粗糙度ψ和残差ξ(Feng et al., 2013；李永生等，2016).在归一化之后，粗糙度曲线表现为单调递增函数，残差曲线则为单调递减函数(图 4).图 4显示，在α²≈3.6，δ≈10.5°时(图 4b白色五角星)，获得log模型的全局最小值，表明最优倾角为10.5°，与Galetzka等(2015)、Wu等(2016)、Zha和Dai(2017)采用的倾角值相近，略大于USGS公布的结果.

图 4

Fig. 4

图 4 (a) 模型粗糙度与拟合残差的折中曲线，其中蓝色短横线表示归一化的模型粗糙度曲线，绿色点划线表示归一化的拟合残差，红色实线表示归一化的log模型; (b) log模型与倾角和平滑系数之间的折中解 Fig. 4 The tradeoff curve between roughness and residual. The blue dashed line indicates the normalized roughness curve, the green dotted line indicates the normalized residual curve, and the red solid line indicates the normalized log model; (b) A compromise map of log model with variations of dips and smoothing parameters

根据均匀滑动反演的发震断层参数结果，利用基于三角位错的断层自动剖分方法进行分布式滑动反演，设置断层走向290°，倾角10.5°，断层尺寸扩展为220 km×160 km，允许滑动方向在逆冲为主兼少量右旋走滑的约束下自由变化.反演过程中采用jRi方法(Barnhart and Lohman, 2010)选择适当的平滑因子，并施加拉普拉斯平滑约束，避免相邻断层片之间的滑动分布解出现振荡.同样地，考虑到InSAR和GPS水平向与垂直向之间的精度差异，这里仍采用赫尔默特方差分量估计方法(崔希璋等，2009；Xu et al., 2009, 2010；许才军等，2016)确定相对权比：

(1) 首先按照InSAR和GPS水平分量与垂直分量观测值验前方差确定初始权阵P₁，P₂，P₃；

(2) 根据断层自动剖分技术进行第一次平差，求得V_i^TP_iV_i；

(3) 进行第一次方差分量估计，求得InSAR和GPS水平、垂直分量观测值单位权方差的第一次估值, 依下式定权：

式中c为任意常数，这里取InSAR数据的单位权方差估值；

(4) 反复进行第(2)和(3)步，直至此三类观测值单位权方差之比近似等于1为止.

通过迭代计算，InSAR和GPS水平分量与垂直分量的相对权比最终收敛至1:8.4:4.0，与均匀滑动反演时确定的权比1:14.5:3.0有差别，主要是由于分布式滑动反演时对断层进行了更精细的剖分.

方差分量估计定权下的基于断层自动剖分的滑动分布结果如图 2.表 2展示了不同研究提供的该地震破裂参数.与已有的研究结果相似，我们的反演结果表明，同震破裂向东南单侧传播约130 km，滑动的优势方向与断层走向接近平行，主要滑动区域在加德满都北部、M_W7.8主震和M_W7.3最大余震之间.深部滑移主要集中在8~18 km，没有明显的地表破裂.最大滑动量5.5 m，位于深度约12 km处，释放地震矩为6.83×10²⁰ Nm，对应矩震级为M_W7.86.最大走滑位错1.2 m，最大倾滑位错5.5 m.同震滑动分布与余震有较好的相关性.余震主要发生在滑动量高值区以外(图 2黄色圆点)，也有少部分余震发生在滑动量高值区域，意味着震间积累的应变能在同震阶段没有完全释放.主震破裂的东北端，即M_W7.3最大余震区域形成了一个滑动空区(slip gap)(即“凹”型滑动缺损区域)(单新建等，2015；Wei et al., 2018)，发生了较多的余震.

表 2 Table 2

表 2 不同研究提供的地震破裂参数 Table 2 The various focal parameters of 2015 Gorkha (Nepal) earthquake 最大滑动量(m)/ 对应深度(km) 主要滑动深度 (km) 最大走滑量 (m) 最大倾滑量 (m) 地震矩 (1020Nm) 矩震级 (MW) 数据 来源 位错 模型 ~8/~16 12.5-22.5 1.2 GPS, ALOS-2, Radarsat-2 InSAR McNamara et al., 2017 三角 位错 5.5/~12 8-18 1.2 5.5 6.83 7.86 GPS, ALOS-2 InSAR 本文a 5.1/~12 8-18 1.0 5.1 6.24 7.83 本文b 6 7.2 7.84 远震波形, Sentinel-1 SAR Avouac et al., 2015 矩形位错 6.5/~16 10~20 7.4 7.8 5-Hz GPS速度波形, GPS, InSAR Galetzka et al., 2015 5.8/~8 6.08 7.79 ALOS-2 InSAR, GPS Wang and Fialko, 2015 5.7 7~15 1.2 5.7 7.55 7.89 ALOS-2 InSAR, GPS Feng et al., 2015 ~6/10.7 8~15 7.8 7.84 ALOS-2, Sentinel-1A, RADARSAT-2, GPS Feng et al., 2016 ~6/~15 12~18 5.3 6.5 7.8 ALOS-2/Sentinel-1A InSAR 李永生等, 2016 6.4/13.83 8~20 7.65 7.86 GPS, ALOS-2 InSAR 谭凯等, 2016 5.3/~15 8~22.5 1.5 5 7.85 GPS Wu et al., 2016 4.9 7.8 GPS Yadav et al., 2017 5.0/12~14 5.5 7.77 GPS, ALOS-2 InSAR Zha and Dai, 2017 ~7 InSAR, GPS, 强震仪和远震波形 Wei et al., 2018 6.73/~7 5~12 ~2 6.5 7.90 GPS, ALOS-2 InSAR Yang et al., 2018 6.0/11 8~20 0.7 6.0 7.43 7.88 本文c 5.3/10 8~18 1.2 5.2 6.36 7.84 GPS, ALOS-2 InSAR 本文d 4.9/10 8~18 1.2 4.9 5.94 7.82 本文e 注:a方差分量估计定权的断层自动剖分的结果; b验前方差定权的断层自动剖分的结果; c基于PSOKINV方法的反演结果; d均匀半空间的SDM方法反演结果; e顾及地壳分层结构的SDM方法反演结果.

表 2 不同研究提供的地震破裂参数 Table 2 The various focal parameters of 2015 Gorkha (Nepal) earthquake

Yue等(2016)对本次2015年尼泊尔M_W7.8地震的滑动反演结果揭示在滑动区域的东北角具有较大的滑动粗糙度和较长的上升时间，证明破裂在这部分受到阻止.Wei等(2018)认为这种滑动模式可能与断层结构有关.Zhang等(2017)通过利用GPS和InSAR数据反演本次地震的2-D倾角变化，表明滑动区域东北角的大倾角异常阻止东-东南破裂向深部传播.现有的地震学证据认为，大倾角异常的原因可能与印度板块的局部撕裂有关(Liang et al., 2011, 2012; Xu et al., 2011; Chen et al., 2015; Pei et al., 2016).

值得注意的是，大部分已有研究都基于单一的倾角对断层面建模，其结果也能很好地拟合观测数据(Galetzka et al., 2015; Feng et al., 2016; Yang et al., 2018).这是由于对于低倾角逆冲断层，抬升形变可以通过小倾角断层上的大滑动或大倾角断层上的小滑动，即滑动量与倾角之间的折中来解释(Zhang et al., 2017).而单一倾角的断层模型也会表现出破裂受阻挡的特性(Galetzka et al., 2015; Feng et al., 2016; Yang et al., 2018)，我们认为可能仍与印度板块的局部撕裂有关，其通过形成具有速度增强摩擦特性的韧性物质障碍体(ductile material barrier)来阻止动态破裂的传播(Yue et al., 2017; 刘静等，2015；姚华建和尹九洵，2015).

图 2和图 5分别显示GPS、InSAR LOS位移场的模型位移值与观测值之间的拟合有非常好的一致性，其拟合残差统计如图 6中红色直方图(具体数值见表 3)，残差分布比较集中，GPS E、N分量拟合残差标准差分别为0.003 m和0.004 m，InSAR LOS位移残差标准差为0.055 m.

图 5

Fig. 5

图 5 (a)(b)(c)和(d)(e)(f)分别代表基于方差分量估计与验前方差定权下的断层自动剖分的InSAR LOS位移场拟合情况，其中(a)(d)表示观测值，(b)(e)表示模型拟合值，(c)(f)表示拟合残差 Fig. 5 The top and bottom panel indicate the observed (a, d), synthetic (b, e) and residual (c, f) InSAR LOS displacement fields obtained by HVCE-FR and WPV-FR method, respectively

图 6

Fig. 6

图 6 不同方法的模型位移值与观测值间的拟合残差统计和正态分布拟合曲线 (a)(d)和(b)(e)分别表示GPS E、N分量的残差统计，(c)(f)为InSAR LOS位移的残差统计.红色直方图和绿色曲线为基于方差分量估计定权的断层自动剖分方法的结果.灰色直方图和灰色曲线为利用验前方差进行权估计的断层自动剖分方法的结果.蓝色直方图和蓝色曲线为基于PSOKINV方法的结果.白色直方图和黑色曲线为顾及地壳分层结构的SDM方法的结果.pdf.为概率密度函数(Probability Density Function) Fig. 6 The statistical histogram of the residuals between the observations and synthetics, with the colored curves indicating the normal distribution fitting curves The left and middle panel represent residuals of GPS E and N components, respectively. The right panel indicates the residuals of InSAR LOS displacement field. The red histogram and the green curve are the results given by HVCE-FR method. The gray histogram and the gray curve are the results derived by WPV-FR method. The blue histogram and the blue curve are the results based on the PSOKINV method. The white histogram and the black curve are the results of the SDM method that considering a layered crust structure. pdf., Probability Density Function.

表 3 Table 3

表 3 不同方法的滑动模型误差以及观测数据拟合残差统计 Table 3 Slip errors and fitting residuals statistics based on different methods 方法 滑动模型误差(m) GPS数据(E/N)拟合残差(m) InSAR数据拟合残差(m) 均值 最大值 标准差 均值 标准差 均值 标准差 HVCE-FRa 0.037 0.141 0.028 0.001/0.001 0.003/0.004 -0.007 0.055 WPV-FRb 0.039 0.116 0.019 -0.002/0.001 0.007/0.013 -0.012 0.080 PSOKINVc 0.040 0.299 0.055 0.007/0.004 0.010/0.014 0.007 0.075 Homogeneous SDMd 0.048 0.139 0.030 0.001/0.001 0.004/0.004 -0.007 0.069 Layered SDMe 0.053 0.128 0.029 0.001/-0.001 0.003/0.005 -0.008 0.067 注:a利用方差分量估计定权的三角位错断层自动剖分方法; b验前方差定权的三角位错断层自动剖分方法; cPSOKINV矩形位错方法; d均匀半空间的SDM矩形位错方法; e顾及地壳分层结构的SDM矩形位错方法.

表 3 不同方法的滑动模型误差以及观测数据拟合残差统计 Table 3 Slip errors and fitting residuals statistics based on different methods

为了评估滑动分布模型的精度，本文采用蒙特卡洛方法，对原始观测数据产生100组带随机扰动误差的数据集，通过这些数据反演相应的滑动分布模型，从而估计模型的精度，得到滑动分布误差如图 7d.可以看出滑动误差分布比较均匀，平均误差0.037 m，最大误差为0.141 m，标准差为0.028 m，证明滑动模型具有较好的无偏性和有效性，表明本文反演获得的滑动分布模型是可靠的.

图 7

Fig. 7

图 7 第一列(a，d，g，j)为利用方差分量估计方法定权下的基于断层自动剖分的同震滑动分布，滑动分布误差，走滑和倾滑分量.红色五角星表示MW7.8级主震及MW7.3级最大余震震中，绿色正方形表示加德满都，红色实线为主前缘逆冲断裂(MFT); 第二列(b，e，h，k)分别表示验前方差定权的基于断层自动剖分的同震滑动分布、滑动分布误差、走滑及倾滑分量; 第三列(c，f，i，l)为基于PSOKINV矩形位错的同震滑动分布、滑动分布误差、走滑及倾滑分量 Fig. 7 The left and middle panel represent the distributed slip model (a, b), slip error (d, e), strike slip (g, h) and dip slip (j, k) derived by HVCE-FR and WPV-FR method, respectively. The red stars indicate the epicenter of the MW7.8 mainshock and the MW7.3 aftershock, the green square indicates Kathmandu valley, and the red solid line is the Main Frontal Thrust (MFT). The right panel represent the distributed slip model (c), slip error (f), strike slip (i) and dip slip (l) based on the PSOKINV method

为了验证方差分量估计定权方法的优越性，本文计算了利用观测值验前方差定权的滑动分布结果，如图 7b.类似地，采用蒙特卡洛方法来估计滑动分布模型的精度，如图 7e.并与方差分量估计定权的滑动模型误差比较如表 3.从两者的滑动分布结果来看，大体滑动特征非常相似，两者滑动分布的空间互相关系数达98.4%.方差分量估计定权时的结果显示同震破裂沿东南方向传播至接近M_W7.3级最大余震区域，分别继续向北和东南传播约30 km，并形成一个滑动空区(slip gap).而单纯按照观测数据验前方差定权的滑动分布结果没有表现出这种显著的“滑动空区”特征，表明利用方差分量估计确定的权比能合理地体现各类观测数据对滑动分布结果的贡献程度与分辨能力，使得保留必要的滑动细节.验前方差定权下的最大滑动量为5.1 m，地震矩为6.24×10²⁰ Nm，最大走滑位错1.0 m，最大倾滑位错5.1 m，与Wu等(2016)结果5 m接近.方差分量估计定权情况下得到的最大滑动量为5.5 m，地震矩为6.83×10²⁰Nm，最大走滑位错为1.2 m，最大倾滑位错5.5 m，更接近于Feng等(2015)得到的最大倾滑位错5.7 m.比较两者的滑动模型误差(表 3)，方差分量估计定权时的滑动误差均值有略微的改善(0.002 m)，滑动误差最大值和标准差稍有增大的现象.但从两者对观测数据的拟合来看，方差分量估计定权下给出的GPS、InSAR数据拟合残差均值和标准差总体上均优于验前方差定权的结果(具体数值见表 3).残差统计直方图(图 6)显示前者残差分布更集中.其中，GPS E、N分量拟合残差标准差分别减小了0.4 cm和0.9 cm，InSAR LOS位移残差标准差减小了2.5 cm.并且值得注意的是，观测到的InSAR LOS位移场显示M_W7.3级最大余震位于沉降区域，且其西侧附近(图 5d方框)有较为显著的抬升(约0.5 m)，方差分量估计定权的拟合结果能很好地显示出这种位移特征(图 5b)，而验前方差定权的模型位移场对上述抬升形变(约0.5 m)细节的拟合稍有平滑的特点(图 5e)，相应的拟合残差较大(图 5f)，且M_W7.3余震显示落在InSAR LOS位移抬升区域，与观测值位移特征相反，证实方差分量估计定权得到的滑动分布模型能更好地解释观测到的InSAR同震形变场.

McNamara等(2017)基于三角位错的断层自动剖分方法，通过施加最小矩正则化平滑约束(minimum moment regularization)，利用GPS和ALOS-2/Radarsat-2 InSAR数据反演此地震的滑动分布，给出的最大滑动量与多数研究结果相比偏大，约~8 m，位于深度~16 km，深部滑移主要集中在12.5~22.5 km，深度超过20 km时仍有2 m左右的滑动.且从对ALOS-2 InSAR同震位移场的拟合来看，对于前文所述的约0.5 m抬升形变细节的拟合结果仍然偏小，约0.1 m.本文基于断层自动剖分方法通过施加拉普拉斯平滑约束，得到最大滑动量5.5 m，对应深度约12 km处，且深部滑移被限制在18 km以内.两者的最大走滑量比较一致，约1.2 m.我们也测试了最小矩正则化平滑约束下的滑动分布，结果显示最大滑动量为5.8 m，大于拉普拉斯平滑约束下的结果，但小于McNamara等(2017)的最大滑动量.为了便于比较，我们测试了增加Radarsat-2 InSAR数据(Feng et al., 2016)进行反演，结果显示，增加Radarsat-2 InSAR数据前后，无论是滑动特征还是最大滑动量值均没有明显差别.Diao等(2015)结果也表明，来自不同卫星(ALOS-2，Radarsat-2，Sentinel-1)的升降轨资料提供了相似的形变特征，增加或去除某种卫星影像，滑动分布结果基本不受影响.因此本文基于断层自动剖分方法的结果与McNamara等(2017)结果之间的差异主要是由定权的差别造成的.

Wang和Fialko(2018)利用基于三角位错的有限元模型(Finite Element Models，FEM)反演了本次2015年尼泊尔M_W7.8地震同震滑动分布.使用分裂节点(split-node)方法实现断层面上的滑动：(1)将整个模型空间网格化以使得地震破裂的两侧具有相同的节点分布；(2)为断层面上具有相同坐标的每个节点对分配一个额外的(“虚拟”)节点，该虚拟节点不属于任何有限元；(3)为每个节点对和相应的虚拟节点定义线性约束；(4)通过沿断层上的滑移方向应用位移边界条件来分裂断层两侧的节点.该模型沿走向，倾向和垂向的维度为1200 km×1200 km×400 km，包含超过一百万个四面体元，网格尺寸朝地震破裂方向递减，断层面被剖分成~3000个三角元，尺寸沿倾向增大，最小尺寸(~3 km)位于断层面上边界，略大于本文基于三角位错的断层自动剖分方法给出的最小三角元尺寸(~2.8 km).得到最大滑动量~7 m，大于本文基于断层自动剖分方法的结果，而两者的滑动特征非常相似.FEM方法考虑了地形效应以及介质结构的横向不均匀性，由于实际的断层几何结构和断层面上摩擦本构关系的复杂性，在有限元处理中如何考虑这些实际情况仍需大量深入的探索，此外，FEM方法比本文的断层自动剖分方法需要更多的计算资源(计算速度、数据存贮、传输等)，如何解决也需继续研究(袁杰和朱守彪，2014).总而言之，本文的断层自动剖分方法在耗费相对较少的计算资源的情况下，在确定模型分辨率矩阵与断层片最佳尺寸的关系的基础上，获得了可靠的滑动分布模型，且不失相应的滑动细节，而考虑地形效应和介质结构的不均匀性，也是今后研究的重要方面.

为了比较基于三角位错和矩形位错的滑动分布模型的差别，我们利用PSOKINV(Feng et al., 2013)和SDM程序(Wang et al., 2013)，基于矩形位错的方法反演此地震的断层滑动分布，设置相同的断层尺寸(220 km×160 km)并均匀剖分成10 km×10 km的子断层、走向(290°)、倾角(10.5°)，并允许滑动角在逆冲为主兼少量右旋走滑的约束下自由变化.同样地，采用蒙特卡洛方法评估滑动分布模型的精度，结果如图 7和图 8.

图 8

Fig. 8

图 8 基于SDM方法的均匀半空间(左)和考虑地壳分层结构(右)的同震滑动分布 (a)(b)为总滑动量; (c)(d)为滑动分布误差; (e)(f)和(g)(h)分别为走滑和倾滑分量. Fig. 8 The left and right panel represent total slip (a, b), slip error (c, d), strike slip (e, f) and dip slip (g, h) based on the SDM method considering the homogeneous and layered crust structure, respectively

与基于断层自动剖分的滑动分布结果相比，它们大体的滑动特征相一致，都表现为破裂向东南单侧传播，滑动优势方向与断层走向接近平行，主要滑动区域在加德满都北部、M_W7.8主震与M_W7.3最大余震之间，没有明显的地表破裂，但在细节上存在差别.基于PSOKINV方法(图 7c)和SDM方法(图 8)的结果相比，两者的最大滑动量分别为6.0 m和~5 m，前者地震矩7.43×10²⁰ Nm，接近于Feng等(2015)的结果7.55×10²⁰ Nm，大于SDM方法给出的结果(6.36×10²⁰ Nm和5.94×10²⁰ Nm).基于SDM方法得到的最大走滑位错为1.2 m，最大倾滑位错~5 m，与Wu等(2016)结果具有较好的一致性.但基于PSOKINV方法的反演结果显示具有更大的倾滑分量，最大为6.0 m，与Jiang等(2018)的结果6.19 m接近，而走滑分量较小，在0.7 m以下.基于SDM方法的深部滑移主要集中在8~18 km，其中无论是基于均匀半空间的结果，还是顾及地壳分层结构(速度结构采用Galetzka et al., 2015)的结果，两者没有显著的差别.其滑动分布结果相对PSOKINV方法的结果而言更加平滑，Wu等(2016)、Jiang等(2018)采用SDM方法的结果均显示出相对平滑的滑动特征.PSOKINV方法的结果显示在20 km深度时仍表现出约1~2 m的滑动.基于PSOKINV和SDM方法结果的上述差异主要是由两者采用的反演算法不同导致的.

本文利用方差分量估计定权的基于三角位错的断层自动剖分方法，结果显示滑动主要集中在8~18 km深度.而基于PSOKINV方法得到的滑动模型在20 km深度时仍表现出滑动特征，滑移量达1~2 m.三角位错的断层自动剖分方法得到最大滑动量5.5 m，接近于基于均匀半空间的SDM方法的最大滑动量5.3 m，小于本文基于PSOKINV方法的最大滑动量6.0 m.最大走滑位错为1.2 m，与Feng等(2015)、McNamara等(2017)、Jiang等(2018)以及本文基于SDM方法的结果相近，略小于Wu等(2016)的结果1.5 m.最大倾滑位错为5.5 m，与本文基于SDM方法的结果及Feng等(2015)、Wu等(2016)、李永生等(2016)结果比较接近，小于Jiang等(2018)结果6.19 m.地震矩6.83×10²⁰ Nm，与本文基于均匀半空间的SDM方法的结果及Avouac等(2015)、李永生等(2016)的结果比较一致，小于本文基于PSOKINV方法的结果.

对于此次2015年尼泊尔M_W7.8地震的研究中，已有研究大多基于矩形位错模型，将断层面均匀剖分成4 km×4 km(Feng et al., 2016；Zha and Dai, 2017)、5 km×5 km(Feng et al., 2015；Wu et al., 2016；Yang et al., 2018)、8km×8km(Avouac et al., 2015；Wei et al., 2018)或10 km×10 km(Galetzka et al., 2015；Yue et al., 2016；Yadav et al., 2017；Jiang et al., 2018)等尺寸的子断层片，或断层片尺寸随深度而增大(Wang and Fialko, 2015).这种基于矩形位错模型的断层滑动分布结果通常都不是对应观测值约束下的最优结果，表现为伪滑动或过度平滑特征，一般情况下二者兼而有之.而基于断层自动剖分方法的结果考虑了模型分辨率矩阵，它可以确定反演得到的断层模型是否适定，在这个基础上，通过高斯曲线拟合确定模型分辨率矩阵与断层片最佳尺寸之间的关系(Barnhart and Lohman, 2010；温扬茂等，2012)，即对断层片的尺寸进行迭代调整，直到得到最佳的断层剖分.本文基于断层自动剖分方法给出的滑动分布模型显示，断层近地表位置的三角元尺寸更小，表明这部分的滑动模型得到了更多的数据约束，具有更好的模型分辨率，而尺寸较大的三角元则分布在断层深部位置，观测数据量相对较少且空间分布远离断层地表迹线，模型分辨率相对较差.这种断层自动剖分方式得到的滑动分布被认为是模型分辨率矩阵约束下的最优滑动几何模型，即每个断层片上的滑移都可以通过足够的观测数据来确定(温扬茂等，2012).

比较三种方法得到的滑动误差分布(图 7和图 8，表 3)，PSOKINV方法的滑动误差平均为0.040 m，最大为0.299 m，标准差为0.055 m.均匀半空间的SDM方法滑动误差平均为0.048 m，最大为0.139 m，标准差为0.030 m.顾及地壳分层结构的SDM方法，其滑动分布误差均值较均匀半空间情况下的误差均值略大，为0.053 m，标准差为0.029 m，最大为0.128 m.与上述两种基于矩形位错的方法(PSOKINV和SDM方法)相比，本文利用方差分量估计的基于三角位错的断层自动剖分方法得到的滑动误差均值和标准差相对最小，均值为0.037 m，标准差为0.028 m，表明本方法得到的滑动模型无偏性和有效性相对更好，从而也表明该模型是更加可靠和稳定的.最大滑动误差0.141 m，不在主要的滑动区域，并且与其最大滑动量5.5 m相比不显著，信噪比较大.这为后续计算地震复发周期(刘静等，2015；谭凯等，2016)、库仑应力分析及地震危险性评估(熊维等，2015；Yang et al., 2018; Rahman and Bai, 2018)、震后黏弹性松弛(Zhao et al., 2017)等研究提供了重要的模型支撑.

从模型位移值与观测值的拟合残差统计直方图(图 6)来看(具体数值见表 3)，基于PSOKINV方法的残差分布相对较分散(图 6蓝色直方图)，GPS E、N分量拟合残差标准差分别达到1.0 cm和1.4 cm，InSAR LOS位移残差标准差为7.5 cm.顾及地壳分层结构的SDM方法的残差分布如图 6中白色直方图，残差分布相对PSOKINV方法的结果而言变得集中，GPS E、N分量拟合残差标准差为0.3 cm和0.5 cm，InSAR LOS位移残差标准差为6.7 cm.比较而言，基于方差分量估计的断层自动剖分方法给出的GPS E、N分量及InSAR LOS位移的拟合残差均值与上述顾及地壳分层结构的SDM方法给出的结果保持在基本相当的水平，略优于后者的结果，其改善主要体现在InSAR LOS位移拟合残差标准差减小了1.2 cm，残差分布相对更集中，也证实该方法得到了相对更优的滑动模型结果.

本文采用基于三角位错的断层自动剖分技术划分断层单元，联合InSAR和GPS同震位移数据反演了2015年尼泊尔M_W7.8地震的发震断层参数及同震滑动分布，反演中基于赫尔默特方差分量估计方法确定InSAR、GPS水平向与垂直向3种观测位移数据的相对权比.得到了如下几点结论：

(1) 2015年尼泊尔M_W7.8地震发震断层以低倾角(10.5°)逆冲为主，兼少量右旋走滑.由方差分量估计定权的均匀滑动反演给出的InSAR、GPS水平向与垂直向拟合残差RMS值，较数据中误差定权的结果分别减小0.8 cm，0.8 cm和0.5 cm.

(2) 由方差分量估计定权得到的三角位错滑动分布模型的最大滑动量5.5 m，大于验前方差定权给出的5.1 m的结果.通过方差分量估计定权得到的反演结果显示在主震破裂东北角形成一个滑动空区(slip gap).而由验前方差定权的反演结果没有表现出这种显著的滑动空区特征.

(3) 方差分量估计定权得到的滑动模型能更好地拟合出观测到的InSAR同震形变细节，其拟合残差标准差较验前方差定权的结果减小2.5 cm，GPS E、N分量残差标准差减小0.4 cm和0.9 cm.

(4) 与两种矩形位错滑动模型结果相比，方差分量估计给出的三角位错滑动分布误差均值和标准差相对最小，均值为0.037 m，标准差为0.028 m，表明该滑动模型的无偏性和有效性相对更好.