地震是地球构造运动的一种表现形式，是自然界最具破坏力的地质灾害之一。通常，强震和大震的孕育发生受到局部应力场以及周边活动断裂控制，导致应力高度积累，在这些应力积累部位及其附近地区会形成明显的构造差异运动，且往往伴有重力场的显著变化（Chen et al，1979；Kuo et al，1999；滕吉文等，2008；Zhu et al，2010）。在地震地球物理观测与分析的学科领域中，一般会采用定期流动重力观测和连续重力重复观测的测量技术，来捕捉整个设定测网覆盖区内与场源变化相关的震前重力场变化信息（陈运泰等，1980；祝意青等，2009；陈石等，2011）。在地下介质变化过程中，重力场的地表场兆反应表现较为明显，且其物理意义明确，是一种较好的构造运动表现形式（高倩等，2015；刘芳等，2016）。如：祝意青等(2001, 2008, 2009, 2013)、陈石等(2014, 2016)将区域重力时间、空间变化分析方法应用于我国川滇地区、新疆地区以及青藏高原东北部中长期地震危险分析研究中，捕捉到震前重力变化高梯度变化特征，为地震危险区的判定提供了一定参考。

对于如何将重力信号用于异常的定量化研究，一直是专家和学者关注的焦点，因此引入地球物理反演方法，基于“以场溯源，场源结合”的分析思路，对已知场源变化信息进行定量化解释，研究重力场变化的场源特征（焦新华等，2009）。在反演过程中需要精确分离地球物理变化过程中蕴含的复杂信息，以获取有效观测和精确的先验信息和参数，这一直是重力场变化反演和解释的难点（高倩等，2015；刘芳等，2016；翟丽娜等，2021）。

本文通过模拟计算理论模型来分析重力场变化，利用欧拉反褶积方法进行反演计算，从而获取适合重力场反演的构造指数和反演参数，最终将三维欧拉反褶积方法应用于川滇地区流动重力数据问题，得到该区场源特征参数，在定性解释基础上，给出定量分析结果，旨在为地震监测和预报提供可参考的思路和方法。

欧拉反褶积方法由英国地球物理学家Reid等（1990）结合Thompson（1982）对欧拉齐次方程的研究总结后提出。在先验信息较少的情况下，可使用该方法实现对重力位场数据的快速反演，通过自动或半自动扫描确定场源位置，找出构造范围并推算其异常源信息（范美宁，2006；翟丽娜等，2021）。计算公式如下

式中：（x, y, z）为观测点位置；（x₀, y₀, z₀）为场源位置；T为位场异常，其在x、y、z三个方向的导数分别为$\frac{\partial T}{\partial x}$、$\frac{\partial T}{\partial y}$、$\frac{\partial T}{\partial z}$；N为构造指数；B为背景场。

构造指数N反映了位场异常伴随场源深度变化的衰减率，参数N的选择是确定场源位置的关键，一般与场源几何结构有关。结合前人研究结果（范美宁等，2006；鲁宝亮等，2009；高倩等，2015；刘芳等，2016；翟丽娜等，2021），在重力异常分析中，N的取值范围一般在0—2。

结合勘探地球物理研究中基于欧拉反褶积建立的超定方程组，计算得到最佳构造指数N，选取该构造指数进行欧拉反演。具体步骤如下：①对离散数据进行网格化处理，②计算x、y、z三个方向的梯度位场异常；③根据异常体性质选择合适的构造指数N和滑动窗口，且窗口滑动尺度最小为3×3；④求解方程组得到异常场源等参数，反复调整，直至获取最佳反演结果。

欧拉反褶积方法在勘探学中应用较广，即使场源信息未知，也可进行场源参数反演和解释。其计算要点在于，在进行反演时，最佳构造指数N与x、y、z三个方向的导数构成线性方程组，求解该方程组即可。构造指数N在重力反演计算中与场源体几何形状相关，实际应用时要根据复杂的地质情况和场源形态进行多次试算，以确定合适的构造指数。构造指数N确定后，根据测点间距和异常形态特征等因素选取合适的滑动窗口，窗口大小能够覆盖场源异常，反复调整参数直至得到最佳反演结果（高倩等，2015；刘芳等，2016；翟丽娜等，2021）。

根据上述对欧拉反褶积方法的认识和应用需求，结合川滇地区重力观测网络形态和规模，设计一个符合研究区观测规模的理论模型，该模型覆盖上地壳、中地壳和下地壳不同深度位置场源的3个直面正六面体模型。该模型如图 1所示，矩形阴影位置为3个场源模型在地表观测网上的水平投影，测网点间距设为5 km，测网长度在x、y两方向上为180 km。其中：M1模型位于中地壳，参数设为：x（50，60）、y（100，110）、z（10，13），中心埋深11.5 km；M2模型位于下地壳，参数设为：x（80，110）、y（40，70）、z（20，25），中心埋深22.5 km；M3模型位于位于上地壳，参数设为：x（80，110）、y（150，155）、z（5，6），中心埋深5.5 km，密度变化幅度选取3.0×10^-3 g·cm^-3。3个场源模型深度不同，图中地表投影区域由矩形框表示，理论异常范围在(0—120)×10^-8 m·s^-2。

图 1(Fig.1)

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

基于上述模型进行反演，构造指数设为N = 1，滑动窗口尺度设为5—8个测点距离，则反演误差深度在15%以内，基本获得最佳反演结果，共计216组解，如图 2所示，实际模型位置由灰色矩形框表示，反演模型结果由黑色圆点表示，计算结果聚集程度较好，不同深度将场源特征明显区分开来。M2和M3模型在水平位置上误差较小，结果较为理想，而M1误差较大。鉴于上述模型的试算结果较好，决定采用以上参数对川滇地区流动重力测网实测重力变化数据进行反演计算。

图 2(Fig.2)

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

川滇地区位于青藏高原东南部（图 3），地质构造复杂，地势起伏较大，地形地貌复杂多样，山谷、河流纵横交错，基岩主要由变质碎屑岩和位于弧后盆地的火山碎屑岩构成。在地质构造上，该区分布在滇南活动地块南段、川滇菱形块体南段和长江中下游断块构造亚区西部，断裂带集中分布在滇南地区，主要有澜沧—耿马断裂带、红河断裂带和小江断裂带等。区内地震活动性强，近年来先后发生汶川M_S 8.0地震、姚安M_S 6.0地震、芦山M_S 7.0地震、盈江M_S 6.1地震、鲁甸M_S 6.5地震、景谷M_S 6.6地震、康定M_S 6.3地震、九寨沟M_S 7.0地震、长宁M_S 6.0等破坏性强震。

图 3(Fig.3)

图 3 川滇地区断裂和流动重力测点分布 Fig.3 Distribution of fault and mobile gravity measurement points in the Sichuan-Yunnan area

沿主要构造带和地震危险区布设流动重力测点，采用绝对重力控制和相对重力联测方式，形成绝对重力控制、相对重力联测的高精度流动重力观测网，目前已持续开展数十期次测量，积累了丰富的重力观测资料，在地震危险性（祝意青等，2015；Zhu，2017）、同震效应（张培震等，2003）、地壳结构（付广裕等，2016）等研究领域得到广泛应用，成为探索和研究川滇地区地壳内部场源信号特征的有效工具和主要手段之一。选取研究区流动重力区域网2013—2017年11期流动重力数据进行整理，重新进行经典平差计算，以获取最佳解算结果。以2014年4月数据为参考基准，计算相邻两期区域重力场的差分动态变化和相对于首期的区域重力场累计动态变化，结果见图 4。

图 4(Fig.4)

图 4 流动重力数据动态变化异常 (a) 1309—1404期；(b) 1309—1504期；(c) 1404—1604期；(d)1404—1609期；(e)1404—1704期；(f)1404—1709期 Fig.4 The dynamic change anomaly of the mobile gravity data

由各期累积动态变化可知，测区重力场累积变化趋势和幅值具有较好的一致性，以罗平、兴义和六盘水等地为中心，周边地区呈显著正异常，周围则形成大范围负异常，正负异常之间形成重力变化梯度带和“零值线”，“零值线”位置和走向随着时间的推移会发生一定程度的转变。如图 4(a)所示为研究区1309—1404期流动重力数据动态变化，在研究时段内发生2014年鲁甸M_S 6.5、景谷M_S 6.6地震，分析2次地震震源区及周边发震前后重力场变化。由图可见，测区以香格里拉、大理、楚雄、罗平、兴义、六盘水、毕节等地为中心，周边为显著负异常区域，周围则形成大范围正异常区域，正负异常区域之间形成重力变化高梯度带和“零值线”，“零值线”的位置和走向发生一定程度的转变弯曲，呈现“中间低、四周高”的整体变化特征。鲁甸地震震中位于正异常区域，靠近正负异常变化高梯度带位置，景谷地震震中位于正异常区域边缘，异常变化特征不明显，见图中蓝色虚线框。

相比于震前异常，鲁甸M_S 6.5和景谷M_S 6.6地震发生后，鲁甸地震震源区正异常能量得到进一步释放，景谷地震震源区发生正负异常翻转，异常变化整体较平缓，见图 4(b)。2次地震主震位置处重力异常变化是地下不同场源、不同埋深和不同地质体的叠加响应，具有多源混合性、非线性、非稳态性特征，力求通过异常分离提取其中有用的异常信息，来探讨重力场变化特征与地下物质密度分布变化关系。

反演得到川滇地区重力变化场源空间特征，结果见图 5。由场源反演位置分布的整体特征可知，在引起重力场变化场源体位置的4次计算结果中，重合度较好区域位于研究区东部的鲁甸地震震中以南和研究区西部的金河—箐河断裂至攀枝花之间地区。与研究区东部和南部相比，西部地区重力分布控制点集中度较差。综合分析认为，异常源为鲁甸以南垂直于昭通断裂的NW向条带状区域，整体场源深度集中分布在20—40 km范围内，与鲁甸地震序列的主震和余震深度分布范围相吻合，也与前人分析结果相一致。

图 5(Fig.5)

图 5 川滇地区重力变化场源空间特征 Fig.5 Spatial characteristics of gravity change field sources in the Sichuan-Yunnan region

基于欧拉反褶积方法，在理论模型分析基础上，选用最优解参数对川滇地区流动重力数据进行重力场反演，得到与场源位置相关的几何参数，分析研究区所选震例发生前后重力场的变化，初步得到以下结论：

（1）欧拉反褶积方法在位场资料解释中应用较广，经试算和分析，发现通过该方法可反演得到重力场异常的场源三维空间位置。通过对研究区重力观测网多期重复观测重力资料进行反演，发现异常变化持续聚集的位置所在区与构造活动性增强区域一致。

（2）根据构建的理论模型和实验结果，欧拉反褶积可通过选择合适的构造指数和滑动窗口，获得较合理的反演结果。借助欧拉反褶积反演得到场源解位置集中区域，往往更接近真实解区域。以2014年鲁甸M_S 6.5地震为例，通过计算川滇地区重力场差分变化，发现重合度较好地区集中在鲁甸以南垂直于昭通断裂的NW向条带状区域，整体场源深度集中分布在20—40 km范围内，与鲁甸地震序列的主震和余震深度分布范围相对应。