伽师6.4级地震前后震源区视密度变化及其构造意义

刘代芹 陈 石 王晓强 张 贝 李 杰 吴传勇 卢红艳

1)中国科学技术大学地球和空间科学学院，合肥 230026 2)新疆帕米尔陆内俯冲国家野外科学观测研究站，乌鲁木齐 830011 3)中国地震局乌鲁木齐中亚地震研究所，乌鲁木齐 830011 4)新疆维吾尔自治区地震局，乌鲁木齐 830011 5)中国地震局地球物理研究所，北京 100081 6)北京白家疃地球科学国家野外观测科学研究站，北京 100095 7)中山大学地球科学与工程学院，广东省地球动力作用与地质灾害重点实验室，广州 510275

0 引言

通过研究重力场的静态异常特征(布格重力异常、自由空气重力异常、剩余均衡重力异常)和动态时间变化(随时变化的内部物质交换、质量迁移等)，可获取地壳内部结构、形态和特性的物理场变化信息(傅容珊，1988；滕吉文等，2006；Crossleyetal.，2013)。地震孕育和发生的整个过程是地壳内部、壳幔物质变动乃至核幔相互作用的结果，可使地壳出现物质运移、应力聚集和释放等活动，这些运动特征将导致地震断裂带内部出现因质量迁移引起的不均衡变化，使得某些活动断裂带具备孕育中强地震的能力及背景。因此，认识地震的孕育及发生的全过程，离不开对地壳构造运动及物质迁移的系统性研究。地下介质结构变化的异常信息可以通过地表重力观测进行捕捉。地球内部由于质量分布不均衡导致的状态变化，常会引起地表观测的重力异常发生变化。研究表明，地震活动区在震前较长的时间范围内将会出现区域重力异常变化信息(祝意青等，2008；申重阳等，2009)。时变重力场异常是与地震孕育发生相关的前兆异常之一。以往的震例，如1997年3月伊豆半岛地震群(Yoshidaetal.，1999)、2001年昆仑山口8.1级地震(祝意青等，2003)、2004年12月印尼9.1级大地震(Hanetal.，2006)、2008年汶川8.0级地震(Zhuetal.，2010；Zhangetal.，2020)、2013年芦山7.0级地震(祝意青等，2013，2015)以及2015年尼泊尔MW7.8地震(Chenetal.，2016)等均表明，地壳内部的物质运移与地震孕育环境的变化有关。在强震发生前数年，地壳能量的汇集或驱散将在一定程度上导致地壳内部区域性的物质运移，特别是对于6级以上地震的孕震过程而言，在地表通常可观测到μGal量级的重力场变化过程。因此，通过开展高精度重力及地球物理探测工作研究地壳深部特征构造的物性结构、监测震前地球物理场异常变化动态，对深入认识地壳深部的孕震环境和识别前兆异常信息具有重要的科学研究价值和现实意义。

本文以 “以场求源、场源结合”的思想为指导，基于地表重复重力观测数据，引入贝叶斯重力平差方法，将仪器的漂移率、气压导纳、潮汐因子、仪器的格值系数等作为未知量，通过贝叶斯原理以及ABIC评价准则选取最优的参数值，在得到平差值的同时可以获取每台重力仪的漂移特性及格值系数等，进而得到时变重力点值序列；在此基础上，采用时变重力信号的等效场源反演方法，最小化局部性浅源高频干扰，并获得研究区地壳内部近10a的区域性重力场源动态变化特征。具体而言，可分为以下几个步骤： 首先，对塔里木盆地西缘的重力观测网络系统、等效源反演方法的原理和研究区的场源监测能力进行介绍；之后，对实际重复重力观测数据进行了测试，并反演得到重力时空场源变化特征；最后，结合2020年伽师6.4级地震的发震构造区重力场变化过程进行分析与讨论。

本文较全面地获得了2020年1月19日伽师6.4级地震(39.83°N，77.21°E)前后柯坪断裂及其周边测网的重力场变化特征，为深入研究该区域深部的地壳物质运移过程和中强地震的孕育机理提供了地球物理依据。所得研究成果对于认识孕震构造背景和断裂带周边构造运动模式具有重要意义。

1 塔里木盆地西缘的重力测网

2005年，新疆维吾尔自治区地震局在西南天山地区初步建立了由40个流动监测点构成的流动重力测网(喀伽网)，为了准确捕捉南天山地区的地震异常信息，分别于2012年、2013年对该重力网进行了优化改造，在原有的喀什-伽师监测网的基础上，将监测范围扩展到西昆仑、塔里木盆地、和田以及库车一带，形成了覆盖阿克苏、库车、喀什、乌恰、塔什库尔干、和田地区，共包括90余个重力监测点的陆地监测网(图1)。

图1 塔里木盆地西缘及周边的地震重力观测网络Fig.1 Seismic gravity observation network in the western margin of Tarim Basin and its surrounding areas.白色四边形为陆地重力观测点，蓝色五角星为绝对重力控制点，黄色圆圈为地震震中，红色实线为活动断裂构造

利用 2 台高精度石英弹簧相对重力仪对监测网进行定期常规复测。测网内由库车、塔什库尔干和乌什3个绝对点提供时空重力基准约束(图1 中的蓝色五角星)。从塔里木盆地西缘的流动重力分布点位(图1)可以发现，测点间距相对不均，测网西北部分布的测点较密。但塔里木盆地内部测点—和田地区一带仅有1条测线通过，监测能力较弱，存在局部监控盲区。

据中国地震台网测定，2020年1月19日伽师MS6.4地震的震源深度为16km，发震构造属于柯坪塔格褶皱-逆断裂带的前缘断裂(柯坪断裂)(温少妍等，2020)。该震区的主体位于柯坪逆冲推覆构造带，由5、6排褶皱-逆断裂带及其之间的新生代盆地或谷地构成(冯先岳，1991；Burchfieletal.，1999；邓起东等，2000；Fuetal.，2003；曲国胜等，2003；杨晓平等，2006；冉勇康等，2006；管树巍等，2007；李安等，2016)。GPS数据研究表明南天山山前-盆地结合带的变形明显大于天山内部，在变形梯度最大的部位历史上曾发生多次7～8级大地震(杨少敏等，2008；乔学军等，2011；李杰等，2015；刘代芹等，2016)。天山地区的逆冲断裂活动引发的地震破坏性极强(张培震等，1996；徐锡伟等，2006)，因此南天山山前盆-山结合部位的构造活动应该是我们需要重点关注的目标(陈杰等，2001；程建武等，2006；冉勇康等，2006)。柯坪逆冲推覆构造系位于塔里木盆地的西北缘(冯先岳，1991；Burchfieletal.，1999；邓起东等，2000；陈杰等，2006)，其北以喀拉铁克断裂与南天山晚古生代层系为界，南以柯坪塔格南缘逆断裂与塔里木盆地内NWW向的塔里木西南坳陷、巴楚断隆和阿瓦提坳陷为界，西界位于八盘水磨一带，与阿图什-八盘水磨逆冲构造带相接，东端位于阿克苏附近，与库车逆冲构造系逐渐过渡(邓起东等，2000；何文渊等，2002；李安，2013；吴传勇等，2014)。柯坪塔格逆冲推覆系自1990年以来共发生6级以上地震16次，2020年1月19日发生的伽师MS6.4地震是继2003年巴楚-伽师6.8级地震后再次发生的6级以上地震，该地震打破了该地区长达近20a的6级地震平静期(王筱荣，2001；高国英等，2007；温少妍等，2020)。

另外，2015年7月3日在研究区重力监测的有效时间范围内还发生了皮山MS6.5地震(37.6°N，78.2°E)，震源深度为10km。该地震位于青藏高原西北端和塔里木盆地的交界带，震中距泽普断裂较近，地震发生在推覆体前缘的皮山背斜上(李金等，2016；吴传勇等，2017)。但该地震周边的有效重力监测点较少，仅有1条测线通过(图1)。

2 场源反演方法

基于陆地时变重力观测研究地壳深部场源视密度结构特征时，通常需要考虑测量误差和场源干扰2种不确定性。由于重力场随时间变化的量级很小，通常在几十μGal的水平，测量误差主要体现在相对重力仪的非线性漂移和格值系数误差2方面。本文基于贝叶斯平差方法(Chenetal.，2019)重新处理了2013年第2期(2013-C2)—2020年第1期(2020-C1)原始重力观测资料(共14期)，获得了重力点值变化结果。

此外，要应用这些点值变化结果研究地壳内部场源的视密度变化特征，还需要分离局部浅源干扰，常见的方法有低通滤波方法、小波分离方法、最小二乘配置、遗传算法和神经网络等优化插值类方法(Chackoetal.，1980；侯遵泽等，1997；Raoetal.，1999；楼海等，2005；柯小平等，2009；陈国雄等，2014；李长波等，2014；许闯，2014)，然而以上方法对于单期资料的处理效果较为明显。

本文基于经典的重力位场反演方法理论(Lietal.，1998)，引入时变重力点值序列光滑的正则约束，以实现多期重力场变化的优化模型计算。采用 “等效源”近似描述时变重力场的场源视密度特征，这种等效类似于卫星重力数据反演的 “等效水厚度”(Schramaetal.，2007)概念。本文反演方法的优点在于，可以通过控制等效源的单元尺度压制高频干扰，通过场源的时空平滑正则约束项剔除局部性、脉冲式的场源干扰因素，进而更有效地获得地壳深部的场源物性特征。

2.1 基本原理

对于时变重力场的等效源反演问题，目标函数可以采用式(1)(Lietal.，1998)描述：

Φ=W0Gm-g(x，y，t)2+Φs(m)+ΦT(m)

(1)

其中，等式右端第1项为观测部分，第2和第3项分别为空间光滑约束和时间光滑约束，具体形式与模型的先验假设有关。式(1)中，函数g表示不同期次的重力观测时空变化点值；x、y、t分别为经度、纬度、时间；G为核函数，与等效源几何参数和观测点位置相关，可预先计算得出；m为待反演的等效源视密度参数；W0为权参数，与观测数据的噪声方差成反比。

下面将具体说明式(1)右端约束部分的构成。

2.1.1 空间约束

式(1)右端第2项Φs可以表示为

(2)

其中，W1、W2和W3是与模型在空间各方向光滑程度相关的权系数(超参数)，其具体的物理含义是期望待求的场源模型在X、XY和Y各个水平方向上都具有二阶光滑特征，具体光滑程度与W1、W2和W3这3个待优化的未知超参数有关。

2.1.2 时间约束

式(1)右端第3项ΦT可以表示为

(3)

其中，W4是与模型在时间上的光滑程度相关的权系数(超参数)，其物理含义是期望待求的场源模型在时间上具有二阶光滑特征，光滑程度由超参数W4控制。

在反演时，需要确定式(1)—(3)中的W0—W4共5个超参数后才可求解。为了合理地确定这些超参数，引入贝叶斯方法中的ABIC准则来实现超参数的优化计算(Akaike，1980)。具体而言，是通过ABIC最小化实现：

ABIC=-2log(L)+2N

(4)

式中，L是模型的似然函数，N是超参数个数。ABIC最小化问题是非线性问题，可以采用单纯形方法或者牛顿法等具有多参数非线性优化能力的方法求解。

2.2 检测板实验

在反演实际重力场观测数据之前，需要对不均匀分布的重力测网进行场源分辨能力测试。本文采用地震层析成像中常用的检测板方法开展场源分辨率测试。由于单独的重力异常不具有垂向分辨能力，因此，本文只进行水平分辨率的检测板实验。具体方法是设计一组正负相间的同尺度场源体，通过重力正演方法得到实际测点位置的理论重力异常，然后利用反演方法获得场源模型参数，并与已知场源参数进行对比以评价由于实际测网分布不均匀导致的空间非均匀场源分辨能力的差异性。通常地震重复重力测量的空间覆盖范围可达几百km，地球的曲率变化不能忽略，故采用球面六面体(Tesseroid单元)作为等效场源的离散化单元。

首先，采用0.5°×0.5°的场源体模型进行离散化，设定正负相间的场源体大小为1°×1°，每个场源体的视密度为±1×10-3g/cm3，场源体埋深为10km，每个场源模型的等效厚度为1km。通过理论计算得到每个等效场源体在地表可观测到的最大重力异常范围约为±32μGal。

图2 是1°×1°分辨率的检测板在地表可观测的理论重力异常与地表实际测点的空间分布情况。由于地表实际的重力测点分布不均匀，对于已知的检测板模型，实际仅能利用有限的陆地重力测网获得观测结果。本文在理论检测板的重力异常值(图2)中再加入5%的噪声，采用上节所述的反演方法分别进行无噪声和有噪声干扰条件下的等效源视密度特征反演，得到的反演结果如图3 所示。

图2 研究区重力测网分辨率检测板Fig.2 Gravity network resolution test board of the study area.

图3 是基于图2 所示的地表实际测网采集的重力异常反演后的场源参数结果。研究区测点间的平均距离为40～60km，因此，每个球坐标系下的六面等效场源体的大小设计为0.5°×0.5°。依据采样定理，采样频率应大于信号最高频率的2倍，因此，选择1°×1°分辨率的检测板进行实验。本文采用的1°×1°实验结果是相对合理的，图3 的相关结果主要用于解释不同区域的重力变化是否可靠提供参考。在图3a中，距离测点较远的区域无法恢复场源参数，在塔里木盆地西端的和田—阿克苏之间也存在大面积的监测空白区。由于每个场源体模型的大小为0.5°×0.5°，而检测板分辨率为1°×1°，故在例如和田西侧的无测点覆盖区内可以看到相邻场源体之间的密度值存在一定的 “渐变”特征，这就是空间平滑正则化约束的结果。在图3b的含噪声反演结果中采用了与图3a一致的颜色标注场源视密度特征，与图3a进行对比可以看出整体上两者区别不大，但在一些局部区域存在一定差异，主要表现为视密度值偏小。产生这种结果的原因是由于加入了随机误差，使得式(1)中与数据相关的估计权重参数W0减小。为了便于之后的分析，在图3b中用黑色虚线标示出具备1°场源分辨能力的区域位置，这对于后续实际数据反演结果的分析和解释具有指导意义。

另外在图3b中可以看出，对于塔里木盆地南部的皮山MS6.5地震，震中周边的场源检测板模型无法恢复；而对于伽师MS6.4地震，地表重力测网对震中周边的场源参数具有较好的分辨能力。因此，本文将重点分析北部的伽师MS6.4地震震中周边区域的场源参数变化特征。虽然本研究资料的时间范围也覆盖2015年的皮山MS6.5地震的时段，但是由于测网对该地震震中及周边监测能力较弱，故不再过多地尝试提取该地震的震前异常特征。

图3 测网反演的检测板图像Fig.3 Image of test board of the survey network obtained by inversion.a 无噪声模型；b 含噪声模型。白色四边形为重力测点的位置；黑色实线为断裂构造；黑色虚线为有场源监测能力范围；黄色实心圆圈为地震震中

3 实际数据测试

在针对实际重力测网进行理论场源模型的分辨率测试基础上，本文进一步对实际观测到的重力数据进行平差和反演。陆地时变重力场一般每年观测2次，分别在3—4月和8—9月进行，本文将上半年和下半年的观测分别定义为C1和C2。

3.1 贝叶斯平差

本文共处理了2013-C2—2020-C1共14期测量数据，通过贝叶斯重力平差方法得到了每期的重力点值和误差估计。首先，以2013-C2期的结果为基准，计算了全部14期的结果中具备同点位重复观测重力点的累积变化图像，选择2015-C1—2020-C1共6期重力时间序列变化数据绘制了图像，如图4 所示。

图4 贝叶斯平差后的重力场变化图像(2015—2020年)Fig.4 Image of gravity field change after Bayesian adjustment(2015—2020).

在图4a中，我们将图3b所示的有效场源监测范围用黑色虚线进行标记。图4a—e为2020年伽师MS6.4地震前的重力场变化图像，图4f为震后的变化图像。对比图4a—e可以看出，在地震前的2017—2018年(图4c，d)出现了较为明显的重力场增加变化，且从2015—2018年的演化过程中可以看出，重力场的时空变化逐年增加，能量不断积累，构造应力不断增强，特别是在伽师6.4级地震之前，2017—2018年(图4c，d)该地区重力值呈快速增加趋势，表明震中附近的区域构造应力突然增强，物质可能不断向该区域运移。但在2019年下半年，重力场开始趋于逐渐减弱的过程，在重力场正负交替阶段即发生了伽师6.4级地震。地震发生之后3个月针对该重力网再次进行了复测，获得了震后该地区的重力数据。从震后重力场变化图像(图4f)可以看出，重力场的时空演变特征继续延续伽师6.4级地震之前(图4e)的变化趋势，即重力场在该网NW侧出现大面积负值变化区域，特别是在震中出现了负值高值区，说明该地区在地震之后能量得到释放，从而导致地壳物质迁移，构造应力处于逐渐减小或地壳密度不断亏损的趋势。从图4 整体可以看出，2020年地震之后的重力场相对于2015年震前的重力场变化结果呈现区域性的减小态势，这与2017—2018年震前区域重力场突增的趋势特征明显不同，这表明伽师6.4级地震发生之前地壳构造应力的积累较为集中，该特征与本次地震的孕育过程密切相关。

3.2 视密度反演

在上节的平差结果的基础上，我们进一步应用第2部分的场源反演方法对全部期次的重力点值平差结果进行了模型反演。一般认为，重力变化与测点周边的环境密切相关，如测点的高程变化、地下水位变化、土壤含水、河流、降雨、蒸发、径流等。如果将这些不确定性统一归类为场源因素，则这些场源因素非常复杂，若要定量地扣除每一种场源干扰，实施难度又极大。因此，本研究将采用等效源反演技术，通过在反演模型引入时空平滑先验条件以压制局部高频干扰、周期性起伏波动变化等特征信号，进而实现深部时变重力场异常信号的提取和场源参数估计。

将与图4 中每个期次对应的反演结果绘制于图5 中，同样在图5a中用黑色虚线标注测网分辨能力较好的位置，且统一采用图5a中的颜色标尺来可视化图5b—f的结果。整体上，图5 所示的重力场源视密度变化特征与图4 的重力点值累积变化的过程相似，图5 中的场源视密度变化范围为±3kg/m3(地壳的平均密度约为2i700kg/m3)，反演得到的等效源层的密度变化约为正常地壳密度的1‰。从图5 的结果可以看出，由于反演模型中的时空正则化约束，反演得到的场源模型时空演化过程更加平滑。

对比图5d和5f在柯坪断裂区域范围内的场源视密度变化可以看出，地震前后出现了明显的与构造走向NEE较一致的区域性场源视密度变化过程。那么，从反演结果与构造形态相似性的角度，在一定程度上可以认为本文得到的反演结果反映了与断裂构造相关的深部重力场源信号。在地震前后(图5d，f)，场源的视密度增加过程由正转负，而且在地震后的视密度变化特征与柯坪构造高度相关，说明这些由重力监测获得场源变化信号与此次地震事件及构造控制的场源环境变化相关。

对比图5a—c可以看出等效源视密度变化由弱转强。首先，在阿图什及周边区域出现近EW向的场源视密度变化特征，自图5d开始变化方向与柯坪断裂系统高度相关，在图5d和5e中阿图什—和田区域范围内出现了1个小幅度的视密度减小区域，这个区域逐渐变大，发震后扩展至整个柯坪断裂系范围(图5f)。

图5 反演得到的地震前后震源区的视密度变化特征(2015—2020年)Fig.5 Variation characteristics of apparent density before and after earthquakes obtained by inversion(2015—2020).

综上所述，MS6.4伽师地震前，场源视密度的变化过程具有继承性发展特征，如果进一步从反演得到的等效源视密度变化的量级上估算，则本文研究区10km深度范围内、1km厚度等效层内约1‰的密度变化可以被陆地时变微重力测量很好地观测到。对于柯坪块体而言，此次伽师地震前后的场源视密度演化特征，未来可作为典型震例辅助判断强震风险源的位置。

3.3 残差分析

按照前述的平差和反演过程，图5 所示的场源视密度反演结果是基于图4 的重力点值平差结果计算得到的，反演的意义在于可以对一些局部高频的场源干扰和周期性的时变场源信号进行压制。两者之间的残差，即通过视密度模型正演得到的重力点值与反演输入的平差重力点值之差，就是反演模型剔除掉的局部场源干扰信号，该信号的直方图特征如图6 所示。

图6 平差重力变化(a)、反演模型重力变化(b)和残差重力变化直方图(c)(2015—2020年)Fig.6 The histogram of adjustment gravity change(a)，gravity change obtained by model inversion(b)and residual gravity change(c)(2015—2020).

从图6a中可以看出，全部14期参与反演的重力平差点值变化范围为±120μGal，而反演视密度模型解释的重力点值变化量约为±80μGal(图6b)，其残差部分的±40μGal重力变化是反演模型剔除掉的局部场源干扰和周期性高频变化。对于本文研究的塔里木盆地西缘地区多年期的累积重力场变化信号而言，这个量级的干扰基本合理。

4 结论与讨论

本文应用最新的时变重力贝叶斯平差方法和时变重力场视密度反演方法，系统性地跟踪分析了2020年伽师6.4级地震前后的重力场变化特征，并反演获得了场源视密度演化特征，得到的主要研究结论如下：

(1)塔里木盆地西缘及周边地区现有的地震重力测网的场源监测能力对于柯坪断裂系统及周边范围具有1°×1°的分辨能力，但对和田西侧构造体系及至阿克苏之间的塔里木盆地内部区域的监测能力相对较弱。

(2)伽师6.4地震前呈现的区域性视密度增加量级约为2kg/cm3，等效源模型的视密度变化约为正常地壳密度的0.74‰，该量级的场源多年累积的重力变化可以很好地被地表重力测网监测到。

(3)伽师6.4地震前的显著重力变化自2017年开始，视密度变化整体上呈现区域性增加的趋势，而在形态上先呈现EW向形态，后逐渐转为NEE向，与柯坪断裂系的构造方向趋于一致。2019年视密度变化趋势减弱，发震后出现NEE向的视密度减小。

(4)伽师6.4地震后出现了与柯坪构造走向较一致的场源视密度减小特征，其主要视密度较小的位置位于地震震中—阿图什一带，这可能与震后断层附近的构造应力得以迅速均衡调整，进而导致的地壳内部流体物质的再分布过程相关。另外，在伽师6.4级地震发生后3个月对塔里木盆地西缘及周边进行重力复测，所得结果表明，由于与震后间隔的时间较短，震中附近区域场源的视密度变化仍延续震前减小的态势，说明2020年4月观测到的流动重力数据有可能包含着此次地震的同震效应信息。

综上所述，地壳内部地质构造形变和地壳密度变化是地壳运动的2种基本形式，且两者相互耦合。地壳介质因形变引起的密度变化产生了重力效应(陈运泰等，1979；申重阳等，2009；陈石等，2014)，地壳物质密度变化是地壳运动的基本形式之一，探测地壳物质质量的重新分布(或物质运移)引起的密度变化，对了解深部物质的分异、调整和运动以及构造分区深层过程的动力学机制具有重要意义。

本文围绕伽师6.4地震前后的时空重力变化特征和等效场源视密度参数开展了一系列分析，所得结果可为新疆地区强震孕育、发生环境的研究提供地球物理场参考。此外，本文研究过程中采用的场源检测板实验方法可以为将来的区域重力测网设计、改进和优化布局提供指导。

致谢本研究基于所有参与天山地区流动重力观测的同事们的辛勤付出才得以开展；中国地震局重力观测技术管理部提供了相关基础资料；GEOIST开源Python软件包(1)https: ∥cea2020.gitee.io/geoistdoc为本文的模型测试和反演提供了支持；审稿专家对本文提出了建设性意见和建议。在此一并表示感谢!