传统的分场方法包括趋势分析、滑动平均和匹配滤波等，它们是一种划分区域场和局部场的“二分”形式，但效果往往不佳。Huang等^[1-2]提出一种自适应处理非线性、非平稳信号的分析方法——经验模态分解(EMD)，它将复杂信号分解成若干个不同频率的固有模态函数(IMF)和一个剩余分量(Res)。该方法较好地刻画了信号的局部特征，它不仅具有多分辨率的特性，且具有分解的自适应性，能够较好地保留信号的原有各种特征。Nunes等^[3-4]在一维EMD基础上，提出一种用于提取图像多尺度特征或空间频率算法——二维经验模态分解(BEMD)，并获得较好的应用效果。

近年来，国内外学者^[5-8]将BEMD方法用于重磁异常分离和滤波，对分解得到的分量异常进行定量解释。然而，目前鲜有利用BEMD方法对重力场动态变化的异常特征进行分离的相关研究。由于BEMD方法在处理分析非线性、非平稳信号方面具有明显的优势，因此本文尝试将该方法应用到川滇地区地震重力资料处理和弱孕震异常特征提取中，获得不同深度层次上的异常信息，这比单纯地分析重力差分和累积动态变化具有较高的分辨率。该方法为地震重力数据处理提供一种新的思路，将其分离得到的异常信息进行深入探讨分析，最终为震情跟踪、地震趋势研判提供重力学依据。

1 二维经验模态分解基本原理

BEMD原理与EMD类似，它将一维EMD中对线的筛分拓展为面，最终得到若干二维固有模态函数(BIMF)和一个剩余分量(Res)。

设存在二维信号s_ori(x, y)(x=1, 2, …, m; y=1, 2, …, n)(m、n为正整数)，该信号经过BEMD方法分解的具体步骤如下：

1) 初始化：Res₀=s_ori(x, y);

2) 提取第j级BIMF，初始值j=1，其中j为分解得到的BIMF个数，其具体分解步骤为：

① f_i－1(x, y)=Res_j－1(x, y)，初始值i=1;

② 查找f_i－1(x, y)的局部极大值点和极小值点;

③ 利用局部极大值点和极小值点分别进行插值计算，得到信号的上、下包络面E_{up, i－1}(x, y)和E_{low, i－1}(x, y);

④ 计算上、下包络面的均值：

$ {E_{{\rm{mean, }}i - 1}}\left( {x, y} \right) = \frac{{{E_{{\rm{up}}, i - 1}}\left( {x, y} \right) - {E_{{\rm{low}}, i - 1}}\left( {x, y} \right)}}{{{\rm{ }}2}}{\rm{;}} $

⑤ f_i-1(x, y)曲面减去上、下包络曲面的均值得到一个新的曲面，记作f_i(x, y)：

$ {f_i}\left( {x, y} \right):{f_i}\left( {x, y} \right) = {f_{i - 1}}\left( {x, y} \right) - {E_{{\rm{mean}}, i - 1}}\left( {x, y} \right); $

⑥ 计算终止条件——标准偏差(SD)，若SD ε(ε为阈值)，则有BIMF_j(x, y)=f_i(x, y)，否则i=i+1，重复上述步骤，其中筛分过程的停止条件SD计算公式为：

$ {\rm{S}}{{\rm{D}}_{ji}} = \sum\limits_{x = 1}^m {\sum\limits_{y = 1}^n {\frac{{|{f_{j(i - 1)}}\left( {x, y} \right) - {f_{ji}}\left( {x, y} \right)|{^2}}}{{{f_{j(i - 1)}}{{\left( {x, y} \right)}^2}}}} } ; $

3) 计算剩余值：Res_j(x, y)=Res_j－1(x, y)－BIMF_j(x, y);

4) j=j+1，重复步骤2)、3)直到Res_j(x, y)极大值点和极小值点个数小于K(K为较小整数)或极大值和极小值的数值近似相等。

最终，二维信号s_ori(x, y)(x=1, 2, …, m; y=1, 2, …, n)通过BEMD方法分解得到K个BIMF和一个剩余分量Res，即

$ {s_{{\rm{ori}}}}\left( {x, y} \right) = \sum\limits_{{\rm{ }}j = 1}^K {{\rm{BIM}}{{\rm{F}}_j}\left( {x, y} \right) + {\rm{Re}}{{\rm{s}}_K}\left( {x, y} \right)} $

另外，BEMD方法在进行信号分解的过程中受到很多因素的影响，其中影响分解过程的关键技术包括：局部极大值和极小值点的提取方式; 包络曲面的插值方法; 筛分过程中的停止准则; BEMD终止准则。为提高计算效率，获得较理想的分解结果，本文数据在进行筛分之前均进行扩边处理，采用邻近窗口法获取极值点、基于Delaunay三角剖分的三次样条插值进行包络曲面的插值计算、标准偏差SD作为筛分终止条件、极大值和极小值点个数小于3作为分解终止条件。

2 模型试验

为模拟野外观测的重力异常创建组合模型，其由5个埋藏较浅、体积较小的球体(剩余密度较小)和一个埋藏较深、体积较大的球体(剩余密度较大)组成，模型参数见表 1。大球体的重力异常分布范围广、幅值大，可视为区域异常; 5个小球体埋藏较浅，异常幅值较小，可视为局部异常。球体组合模型正演的局部、区域和叠加异常见图 1。

接下来，根据BEMD方法对叠加异常进行分解，由于在利用曲面插值方法进行边界拟合过程中，边界可能会出现发散和失真的情况，再经过多次插值迭代会向内部发散而污染整个数据。为减小在筛分过程中边界效应对分解结果的影响，对其进行数据扩边处理。通过对比分析原始异常的特征，异常特征变化相对简单，筛分条件和分解终止条件无需过于苛刻，为提高分解计算效率，设置筛分过程中阈值ε=0.01，对模型叠加异常进行分解共得到两组数据，分别为第一级固有模态函数(BIMF1)和剩余分量(Res)，对上述分解结果异常进行缩边处理，最终得到BEMD分解结果，见图 2。

从图 2(a)可以看出，分解得到的BIMF1中最明显的为5个小球体的异常，与模型理论局部异常对比，两者异常形态、位置一致，异常幅值较理论值偏小; 图 2(b)为分解得到的Res，其反映的是埋藏较深的大球体产生的异常，异常的形态、位置也与理论区域异常一致，幅值较理论幅值偏大些。定性来说，利用BEMD方法对重力数据进行分解是有效的，能够保持获得准确的异常形态且幅值不发生较大的畸变。

另外，计算模型叠加异常、BIMF1和Res的径向对数功率谱曲线，见图 3，其中虚线为拟合直线。图 3(a)可以看出，拟合出两个直线，说明叠加异常中存在两个不同频率的异常; 同样地，图 3(b)、3(c)中也拟合出两条直线，与图 3(a)中的拟合直线斜率大小相等，存在一一对应关系。利用图 3(b)、3(c)中拟合的直线斜率分别估算各分量反映地下场源的似深度，分别为10.85 m、74.64 m，这与模型理论值10 m、80 m基本吻合。因此，定量来说，BEMD方法可以将叠加异常中的两个不同频率分离出来，各分量具有一定的实际地质含义，主要体现在场源深度上。

3 数据处理与解释

川滇地区位于青藏高原东南部，地貌起伏较大，地形地势复杂，山谷和河流纵横交错。研究区域内断裂构造较发育，包括澜沧-耿马断裂带、无量山断裂带、红河断裂带和小江断裂带等。由于川滇地区位于欧亚板块与印度洋板块中国大陆碰撞带东缘、造山运动的突出位置，地质及其相关的地球物理背景比较复杂，近年来该区域受到国内外学者的密切关注，复杂的构造环境与强烈而频繁的地震活动规律成为该区域的研究热点^[9]。

当前往往通过对地震重力异常形态和特征进行分析研究，为地区震情跟踪、地震活动趋势分析和强震危险性提供意见。然而，川滇地区地壳构造运动和地震活动强度剧烈，区域重力场变化幅值较大，诱发地震“弱异常”信息可能湮没在背景场中，弱化了孕震异常，当然，强震孕震异常也可能存在重力趋势变化异常中^[10]。另外，随着陆态网络高精度绝对重力基准的建立，为相对重力联测提供了高精度起算基准点，可获得真实可靠的平差计算结果，避免重复重力观测获得的重力场动态变化出现畸变^[11]。故以测区内贵阳、泸州、水城等高精度绝对重力点作为起算基准点，对每期相对重力联测数据进行平差计算，再经过固体潮、零漂和一次项等改正，得到每期各测点的重力值，平差后的各期点值平均精度均优于12 μGal。由于1 a尺度差分的重力异常可以较好地消除季节性影响，较为客观地反映重力场变化特征，因此本文利用BEMD方法对1 a尺度2015-09~2016-09重力场差分动态变化(图 4)进行自适应性多尺度分解，然后对分离出来的各分量异常进行研究。

在利用BEMD方法处理数据之前，对异常进行扩边处理，为提高计算效率并获得较丰富的异常细节信息，筛分过程选取阈值ε= 0.03，共得到3组数据，见图 5，各级分量异常波长随着分解阶次的增加而增大。从分解结果来看，BIMF1表现为异常分布零散，分布范围相对较为集中，西鱼河-昭通断裂和威宁-水城断裂位置处异常幅值可达25 μGal，百色-合浦断裂位于正负异常的过渡位置，其异常幅值变化较大，可能反映地下较浅岩体、构造和地质体等的物性和产状变化特征; Res1为原始异常去除BIMF1后的剩余项，下同; BIMF2中异常分布范围较广，沿百色-合浦断裂和安宁河断裂呈现条带状正异常分布，异常幅值为17 μGal左右，而威宁-水城断裂正位于正负异常过渡带及梯度带附近; BIMF3为低频异常，异常呈北北东、西北北向条带状异常，整个测区异常幅值变化达30 μGal，反映地下较深构造变化特征; Res3为原始异常的趋势项，反映该地区的2015-09~2016-09重力差分动态变化的区域变化特征，异常幅值变化范围较大，鲜水河断裂、小江断裂和安宁河断裂位于负异常区，文山断裂、开远-蒙自断裂为正异常区，异常幅值最大为45 μGal，正负异常梯度带位置等值线发生一定程度的转折弯曲，从短时间尺度来看，该区域正处在重力正异常的能量积累。

另外，计算上述各级BIMF和剩余分量的径向对数功率谱曲线，见图 6，并对各谱曲线进行直线拟合，利用拟合出的直线斜率估算各级BIMF和Res3对应的场源似深度，分别为19 km、36.7 km、64.2 km和81.8 km。由此看出，重力数据经BEMD多尺度分解，各级BIMF反映的场源似深度随着分解阶次的增加而增大。因此，在BEMD分解结果中，BIMF1、BIMF2和BIMF3为高频局部重力变化，其大致反映中上地壳地下地质体、地质构造等密度变化情况; Res3为低频区域重力变化，其大致反映地层、基底界面起伏变化状况。所以，将BIMF1、BIMF2和BIMF3三者相加作为原始异常的局部场，Res3反映出异常的趋势项，将其作为异常的区域场，最终实现异常的“二分”，见图 7。从图 7看出，湮没在区域异常的弱异常在分离出的局部异常相对得到增强，异常特征更加明显; 同样，分离出的区域场正负异常过渡带或梯度带位置明显，更能反映出原始异常的区域特征。

另外，利用传统的匹配滤波方法对原始异常进行分场，结果见图 8。BEMD分场得到的局部和区域异常与匹配滤波得到的结果进行对比分析，两种方法得到异常形态和幅值变化都是一致的，能够较好地反映该地区重力差分动态变化的局部和趋势特征。

综上可知，BEMD方法不仅保留了传统“二分”法的优势，还可以将不同频率的异常从原始异常中分离出来，得到不同深度层次下的异常信息，然后再经过简单的组合叠加可得到局部异常和区域异常，实现分场。因此可将BEMD方法应用到地震重力异常分场中，然后通过对分解得到的异常进行分析研究，可使得弱化的重力异常信息得到增强，异常趋势变化特征更加明显。

4 结语

1) 通过模型试验验证，定性半定量地解释了BEMD方法分解得到的各分量的实际地质含义。

2) 在异常特征提取之前，对异常作特定的扩边处理，BEMD方法分解效果会得到进一步提升，不至于随着筛分过程的进行，内部数据被“污染”而导致分解结果“面目全非”。

3) BEMD算法可以提取不同尺度下流动重力异常特征信息，也可通过重构实现异常“二分”，这比纯粹的动态变化异常具有较高的分辨率，为流动重力数据处理、分析与解释提供了一种新的思路。