布格重力异常的产生与地球内部物质不均匀分布有关，一般可以作为静态场来研究地质构造特征、地壳物性界面结构及地壳均衡状态等地质时间尺度的地学问题(陈石等， 2011a，2001b; 陈石和王谦身，2015)，但是对于较短时间尺度的重力场变化，通常需要借助更高精度的定点重复重力测量来获得.在地震监测和地震前兆学科领域，研究人员采用流动重力测量技术来获取整个测网覆盖区域范围内的重力场变化信号.重力场变化作为地下介质过程在地表的场兆反映之一，具有明确的地球物理意义，并可以认为是构造运动的一种表现形式.根据以往的研究结果表明，区域性的重力场变化梯度带对应潜在地震危险区(祝意青等，2009；申重阳等，2010).而在孕震期的重力场变化，在地下对应一个质量变化区域，这个区域可以称为“震质源”，其在地表的投影称为“震质中”(Kuo and Sun， 1993;顾功叙等，1997；郑金涵等，2003).通过重力场变化梯度带等信息，可以定性判断区域构造活动特征，为判断地震危险区提供参考.

但是如何应用地表观测的重力变化信号，定量开展对异常源的研究，即“以场求源”问题，一直是大家关注的问题.通常思路是采用地球物理反演方法来研究引起重力场变化的场源特征.一般认为，我们在地表观测的重力场变化，通常包含测量误差、近地表水循环(季节效应)、测点高程变化和地下介质密度变化及可能的物质运移等信息(陈运泰等，1980;陈石等，2011c).但是想精确分离这些地球物理过程在重力场变化中的贡献量，常常十分困难，这主要是由于与这些过程或深部介质变化的很多先验信息及参数无法得到有效观测或精确获得.缺少对地下介质和过程的先验认识，已经成为了进一步开展重力场变化反演和解释的主要难点之一.另外，就是流动重力观测获得到的重力场变化信号量与观测误差相比，信噪比并不高，如何从这种重力微变化信号中，提取与场源相关的物理信号，也是比较困难的问题.

在地震前兆学科研究中，流动重力观测方法是形变学科手段之一.在以往重力场变化问题研究中，数据计算和解释多采用大地测量学的技术手段，而很少应用勘探重力学方法.本文尝试使用在勘探重力学领域和位场研究中广泛采用的反演和解释方法，即位场数据的三维欧拉反褶积方法(Reid et al., 1990，Reid，1995)，来研究重力场变化数据.该方法适合在先验信息缺乏的情况下开展位场反演和解释.

本文将以川滇交界地区的流动重力测网实测数据为研究对象，通过对该区域2012-2014年期间实测的流动重力观测数据分析，以2012年10月期间的重力测量数据为基准，研究该时间范围内的多个累积重力场变化的场源特征.通过设计了一个与实际流动重力观测相似的理论模型模拟地表观测的重力场变化，并利用正演手段模拟观测到的重力场变化，再应用三维欧拉反褶积方法对场源参数进行反演，优选与模型最相近的构造指数及主要反演参数；最后，在模型优选参数的基础上，对川滇交界地区的流动重力测网5次重复测量获得的数据进行了反演和解释.本研究的意义在于，尝试将三维欧拉反褶积方法用于反演重力变化信号的场源位置参数，并开展了模型实验和实际应用研究.依据本文的研究方法和初步结论，可以在已知先验信息较少的情况下，开展“以场求源”的定量研究，可以为区域重力场变化信号的反演和解释提供一条新思路和方法.

1 方法原理

欧拉反褶积方法适合在已知先验信息较少的情况下，开展自动或半自动的场源位置反演，圈定构造范围和解释异常源的特征参数等.该方法最早由Reid提出(Reid et al., 1990)，主要是基于欧拉(Euler)齐次方程(公式1)建立的，通过求解欧拉方程可以确定场源位置的三维空间位置，可以用于对不同类型地质体给出参数估计，公式为

在应用三维欧拉反褶积方法进行反演时，需要首先对离散数据进行网格化，然后计算位场异常在x、y、z三个方向的梯度，并选择适当的构造指数N.通过对观测数据网格指定一个大小的数据窗口，采用窗口滑动的方法进行反演.窗口每个观测点数据带入公式1可以组成线性方程组，选择合适的窗口大小，要求解1式所示方程组窗口大小至少应该大于等于3×3.每次移动一次窗口，都可以求解得到一个场源位置参数，最后可以根据这些解的汇聚位置，来确定最佳反演结果，或者调整反演参数.

在选择合适的窗口大小时，应该使窗口大小覆盖异常范围，在一个合适窗口内可以覆盖一个场源的异常范围.在实际的应用中，常常需要通过计算网格尺寸来估计异常规模形态，通过对参数进行反复调整，最终得到一个最佳的反演结果.

2 模型反演

欧拉反褶积方法不需要已知场源先验信息，即可以对场源参数进行反演和解释，因此在勘探地球物理中被广泛应用.这主要是由于在欧拉齐次方程中场源物性参数在方程两边被相互化简，方程给出了位场异常、梯度和场源位置三者之间的相互关系.因此，反演计算前只要指定构造指数N和求出异常在X、Y、Z三个方向的导数即可组成线性方程组进行求解.但是，对于不同的地质问题和模型需要选择不同的反演参数，即构造指数N和滑动窗口大小.

构造指数是与不同的场源几何形体有关，通常对于岩脉、断层、台阶等问题，在重力反演问题中都对应一个相应的理论构造指数N.但由于实际地质构造情况复杂，场源形态也不确定，因此，如何确定和选择构造指数，通常还要通过对不同模型进行反复测试来确定.而滑动窗口大小的选择，一般要考虑测点间距、异常形态特征等，原则是在一个窗口内应该有足够多的测点可以满足方程组求解要求，并可以覆盖一个场源异常形态等.对得到反演结果评价方法，需要通过观察计算求得场源位置参数的水平误差和垂直误差大小，反演场源位置解的三维空间聚集程度和与地质认识是否符合等标准判断.

基于以上对该方法的认识和解释震前重力变化场的需求.本文在将三维欧拉反褶积方法应用于解释实际流动重力资料之前，根据现有的重力观测网络形态和规模，及孕震区空间特点进行了估算.设计了一个比较符合现有观测能力的理论模型，并根据重力异常正演计算方法，预先设计了在上地壳、中地壳和下地壳不同深度位置的场源模型，正演了在地表测网上的理论异常值，用于模拟观测到的重力变化信号和优选反演参数.模型异常如图 1所示，在图中矩形阴影位置为三个场源模型在地表观测网上的水平投影，地表测网点间距为5 km，测网长度在X、Y方向都为180 km.三个不同深度的场源模型参数如表 1所示.

图 1

Fig. 1

图 1 不同深度场源组合模型及正演异常 Fig. 1 The combined model of multi-depth bodies and their forward anomaly

表 1(Table 1)

表 1 理论场源模型参数 Table 1 Model parameters of theory field source 场源模型 直立六面体角点坐标(km) 密度变化(g·cm-3) x1 x2 y1 y2 z1 z2 M1 80 110 150 155 5 6 0.003 M2 50 60 100 110 10 13 0.003 M3 80 110 40 70 20 25 0.003

表 1 理论场源模型参数 Table 1 Model parameters of theory field source

用于模拟上地壳异常源的M1模型单元尺度为20 km×5 km×1 km，中心埋深5.5 km；用于模拟中地壳异常源的M2模型单元尺度为10 km×10 km×3 km，中心埋深11.5 km；用于模拟下地壳异常源的M3模型单元尺度为30 km×30 km×5 km，中心埋深22.5 km，三个模型的密度变化统一取为3.0×10^-3 g·cm^-3.三个模型的尺度特征各具差异，埋藏深度也不同，在地表在这组模型计算的理论异常异常范围在0-120×10^-8m·s^-2之间.

通过反复试验计算，我们选择了构造指数N为1时，结果与理论模型更加接近，反演结果的聚集程度更符合模型的几何特征.而滑动窗口大小在5~8个测点距离之间反演效果较理想.整个模型的反演深度误差范围限制在15%以内.图 2是我们根据以上原则得到一组最佳反演结果，图中黑点为反演的模型结果，灰色矩形框为实际模型位置，反演得到的216组解，可以明显看出结果的聚集程度较好，在三个深度范围内可以明显地区分出场源特征.在水平位置上，M2模型反演结果存在一定的误差，但M1和M3模型较为理想.因此，基于上述模型测试，我们选择该参数用于反演在川滇交界地区的流动重力测网实测重力变化资料.

图 2

Fig. 2

图 2 基于欧拉反褶积的反演结果三维空间分布 Fig. 2 The 3D inversion results using Euler Deconvolution method

一般认为，地震是构造运动的一种表现形式，地震的空间位置和震源机制是岩石圈动力学过程的反映，通过地震可以了解地下深部的地质构造和应力状态.川滇交界地区地处我国大陆活动地块边界带，构造背景复杂，中强地震频发.而文本的研究区位置如图 3所示，该区也多年来一直被划为地震重点危险区，区内的安宁河、则木河、小江断裂等在新生代以来都处于活动状态.由于研究区范围地形多以山地为主，原有流动重力测量多沿平行于断裂的沟谷进行，为了更好地观测活动断裂周边的地球物理场变化特征，我们在原有四川地震流动重力测点的基础上对重力测网进行了适当的加密和拓展，构成了由71个有效测点组成，覆盖研究区内主要活动断裂的观测系统.从2012年10月起，开展了多期的流动重力观测.图 3中黑色四边形为实际的流动重力测点位置分布，测点的东西向分布横跨了多条南北向断裂系统，这对断裂构造活动引起的形变异常监测是十分有利的条件之一.在研究区内，断裂带走向主要以北东向和南北向为主.比较重要的断裂系统，在图 3中都已经标出.特别是其中的安宁河、则木河及小江断裂，在历史上都发生过M_S 7.0以上的大地震.

图 3

Fig. 3

图 3 川滇交界地区的流动重力测点和断裂构造分布 Fig. 3 The mobile gravity mesh，fault system and history earthquake in the Fujian province

一般认为，重力位场异常可以体现地质构造运动过程，常被用于划分构造单元和识别断裂带形态、盆地边界等地质特征.自由空气重力异常可以被认为是表征整个地下物质盈亏程度的标志；区域布格重力异常可以反映深部构造形态，被称为构造运动的“活化石”.在研究区范围内，本文使用了EGM2008模型给出的全球自由空气重力异常数据(Pavlis et al., 2008，2012)，该数据空间形态特征如图 4a所示.在该数据基础上，使用了公开的全球海陆联合高程数据模型(Smith and Sandwell， 1997)进行了布格校正，再通过位场异常分离技术，剔除了浅部高频异常后，得到了如图 4b所示的区域布格重力异常模型.

图 4

Fig. 4

图 4 川滇交界地区的重力异常与历史地震分布 (a)自由空气重力异常；(b)区域布格重力异常. Fig. 4 Gravity anomaly and Earthquakes in the boundary of Sichuan and Yunnan

在图 4a给出的自由空气重力异常图中，异常范围在-220~260 mGal(10^-5m·s^-2)之间.该异常与地形具有较强的相关性，在山地地区异常多为正值，平原或盆地地区异常表现为负值.图 4b给出的区域布格重力异常，范围在-360~200×10^-5 m·s^-2之间，整体为负值异常，其中西北侧靠 近青藏高原地区为布格异常最小区，而向东南异常值逐渐增加，异常整体形态受断裂构造控制.在图 4中，我们标注了2013年以来的4级以上地震震中位置，其中，4至5级以上地震10次，5至6级以上地震1次，6级以上地震1次，从这些地震的震中位置分布上看用现在已知的地表断裂尚不能完全解释全部地震活动，而重力异常场的趋势变化特征与地震活动也具有一定的相关性.这表明这些重力异常背景场，尤其是布格重力异常场中也包含了一些深部构造分界信息或者可能隐伏在深部地壳内断裂的地球物理场特征.

根据川滇交界地区的流动重力测网在2012-2014年的5次测量成果，本文给出了如图 5所示的四期累积重力场变化结果.图 5a-5d都是以2012年10月为起算的累积重力异常，图中黑色等值线为重力场变化，其中虚线为变化负值区，实线为变化正值区，等值线间距为5×10^-8 m·s^-2.由于重力测点的空间分布不均匀性较大，本文采用了最小曲率方法，进行了空间插值，得到的等值线形态可以较好地反映实际测点的空间分布，等值线范围与测网覆盖区域统一，但对于没有实际重力测点覆盖区域的重力场变化，我们应该谨慎的解释，本文的分析和认识也仅对有实际测点覆盖的区域进行.

图 5

Fig. 5

图 5 川滇交界测网2012-2014年间重力场变化 (a)2013/03-2012/10年重力变化；(b)2013/08-2012/10年重力变化；(c)2014/03-2012/10年重力变化; (d)2014/06-2012/10年重力变化;图中红色实线为断裂构造位置，黑色线为重力变化等值线，实线为重力变化 正值区，虚线为重力变化负值区，等值线间距统一为5×10-8 m·s-2，图中蓝色圆圈为2013年 以来的历史地震震中. Fig. 5 The gravity variations in the boundary of Sichuan and Yunnan between 2012 to 2014

从图 5给出的四幅重力场变化图像形态特征看，容易发现，其中，重力场变化的信号成份较复杂，既有小尺度的局部异常，也有大范围的趋势性变化.图 5a中2013/03-2012/10期间的重力场变化在以则木河断裂和小江断裂为分界，东西的重力场变化模式存在一定差异，断裂以东正变化为主，以西则出现北正南负的变化特征.图 5b中2013/08-2012/10期间的重力变化整体幅度小于图 5a，空间特征上主要以安宁河断裂为分界，东侧正负变化相间，西侧以正变化为主.图 5c中2014/03-2012/10期间的重力变化整体呈现正变化特征，垂直于昭通断裂的北西向重力变化梯度带呈现一定的规模，另外在则木河断裂、小江断裂和莲峰断裂的交汇处出现较复杂的局部变化特征.图 5d中2014/06-2012/10期间的重力变化整体呈现负变化特征，垂直于昭通断裂的北西向出现了一个正负相间的重力变化梯度区，形态上与图 5b和图 5c相似，但变化量值不同.

从图 5a-d的四个测量期间的累积重力变化特征可以看出，以同一个时间起点为基准的不同时间尺度的累积重力场变化特征，存在一定一致性，但也具有明显的差异性.一致性特征在于垂直于昭通断裂的北西向重力变化梯度带呈现的位置基本一致，差异性在于不同时间尺度的累积重力变化，量值差异较大，在一定程度上与地表断裂位置相关.基于以前的研究表明，在一定规模尺度上的重力年变化存在正负波动的变化特征，这种波动具有规模上相近但反向变化的特点.仅就目前的观测数据还不足以有效地分离各种场源引起的重力场变化.如果我们将近地表的局部影响、数据测量的不确定性，作为一种随机过程来考虑，那么不同时间尺度的累积重力场变化相互对比，在地壳深部场源的重力变化特征应该具有一致性和相似性.而通过对不同时间段内的反演结果进行相似性分析，那些一致性较好的区域，应该是较可靠的异常源位置.

因此，如果可以从反演得到场源特性上去认识重力变化的机理问题，应该更多的关注那些多期重力变化反演结果场源重复性好的区域，对于那些聚集性好的位置，应该是存在物质运移较快，构造反复活动的区域，应该予以重点关注.

基于以上的原理，我们对不同累积时间尺度的重力场变化进行反演，目的就是通过对比结果之间的共同特性去发现较可靠的异常源，这也是目前处理此类微重力变化或低信噪比数据的有效思路之一.

本文采用三维欧拉反褶积方法，反演得到了如图 6所示的四期重力变化的场源形态特征.图中用圆圈表示反演的场源位置，大小与场源深度相对应，分别采用了绿、蓝、粉、桔黄四种颜色表示不同的累积时间反演解.在欧拉反褶积计算过程中，数据网格化间距为7 km，窗体大小选择为8倍间距.得到的反演解个数略有不同，2013/03-2012/10期间的得到323个反演解，2013/08-2012/10期间的得到387个反演解，2014/03-2012/10期间的得到363个反演解，2014/06-2012/10期间的得到390个反演解，反演结果的场源深度范围在5~70 km范围，与该区地壳厚度相关.根据上述对近20个月的重力变化场定性认识，我们可以从这些反演结果上进一步定量地分析场源特征，为进一步对孕震期震质源特征描述和分析解释提供依据.

图 6

Fig. 6

图 6 川滇交界地区2012-2014年重力变化的场源特征空间分布 图中黑色实线为断裂构造位置，灰色四边形为流动重力测点，蓝色三角形为主要城市位置，红色实心圆圈为地震 震中分布，不同颜色的空心圆圈表示不同时间段重力场变化的反演结果. Fig. 6 The spatial distribution of inversed field source using the gravity changes from 2012 to 2014 in the boundary of Sichuan and Yunnan

从图 6反演的场源位置分布整体特征上看，引起重力场变化的场源体位置在四次结果中，重合度较好的地区集中在东部以鲁甸地震震中以南位置和西部金河-箐河断裂至攀枝花之间的地区.但是，西部集中度较好的区域范围内缺乏有效的实际重力测点控制，因此，本文认为比较可靠的异常源应该是鲁甸地震震中以南呈北西向垂直于昭通断裂的条带状区域内.通过统计该区域的反演结果显示，场源深度主要集中在30±10 km范围内，这与鲁甸M_S 6.5地震的主震及余震分布也是比较一致的.

另外，我们除了对比不同时间段反演结果的空间集中度外，还应该考虑反演结果的深度一致性，反演结果的深度分布一致性较好的异常位置应该更可靠.而对于集中度不好的位置，我们认为这些场源的变化具有一定的随机特性，主要将其解释为这些地方虽然在短期内可能存在一定的场源(物质)运动，但是不具有一定的持续性，应该属于一种较短期的地壳活动现象或者与观测误差的不确定性有关.

4 结论和讨论

本文采用三维欧拉反褶积方法，从理论模型测试到实际资料反演和解释.尝试将该方法应用于解释重力变化数据.由于该方法不需要先验信息，即可以用于估算位场异常的场源参数，因此，适用于在不进行重力变化场源分离的情况下，开展对重力变化场源特征的反演和定量解释工作.根据本文的研究结果，主要结论如下：

(1)欧拉反褶积方法广泛应用于位场资料的解释中，也适用于反演重力场变化异常的场源三维空间位置特征，对于多期重力变化反演结果持续聚集的位置可能揭示了区域构造活动性较强.

(2)根据本文构建的理论模型和实验结果，欧拉反褶积可以通过选择合适的构造指数和滑动窗口大小，不需要给出准确场源的物性参数即可以得到较理想的反演结果.

(3)2014年8月3日鲁甸M_S 6.5地震的余震分布和发震构造(徐锡伟等，2014)为北西向，震源机制为走滑型，与近震中的北东向昭通逆冲断裂不一致，但是，在区域布格重力异常上存在明显的北西向异常梯度带，这说明布格重力异常可以揭示一定的深部壳内隐伏断裂或地壳横向不均匀性特征.

(4)而根据对川滇交界地区的流动重力实际资料处理和反演发现，在地震前的2年时间尺度范围内，出现了与布格重力异常较一致的北西向异常源活动特征，多期累积结构具有类似特性；这与已有鲁甸地震震中及周边的密度结构反演结果(陈石等，2014)和震前的近震中区域等效密度变化反演研究结果(石磊等，2014)相一致；但是，区域内其他小于6级地震前并没有发现类似的场源活动特征，这可能与本测网的场源分辨能力有关.

(5)对2012年10月以来多期的重力累积变化反演结果表明，在鲁甸地震震中以南的欧拉反演解重合度较好，其场源的深度分布主要集中于30±10 km范围，在该区属于中上地壳深度，可以认为是一种孕震深度范围的场源特征.

根据本文的研究结果，从方法原理上看三维欧拉反褶积方法合适于解释重力场变化问题.由于引起重力场变化的场源本身十分复杂，很难估算场源的合理密度参数范围及深度分布，这就决定采用常规的三维密度结构反演技术，进行密度结构变化的计算存在很大技术难度.而欧拉方法的主要优点就是不需要知道太多先验信息，这就使得该方法在解释重力场变化场源特征上具备独到的优点.通过对每期重力变化资料的跟踪和反演，重点可以关注反演结果中场源解位置聚合程度较高的地区，对不同测网和地质构造情况下，应该先通过模型实验，优选反演参数，多次反演综合评判反演结果的可靠性.本文的研究方法和结果，对于认识重力场变化信号的场源特征可以提供有效的帮助，适合于开展定量研究，可以为流动重力变化信号的反演和解释提供一条新思路.