地质构造的边界识别是指对地下断裂构造、块体边界等线性分布的地质构造在地面上的展布位置进行定位，判断地下不同物性差异（如密度、磁化率和电阻率等）的构造带、不同岩性差异地质块体等构造的线性特征位置分布和变化范围。在不同密度或磁化率变化差异的地质构造边界位置，重磁数据异常变化较为明显，因此位场数据被广泛应用到地下地质构造的线性特征增强和提取中，且有较好的应用效果^[1-2]。目前，位场边界检测方法主要包括数值计算类和数理统计类两大类，均需导数或高阶导数运算，观测数据若存有干扰噪声，会放大计算过程中的数据误差，影响线性边界的准确定位。此外，小波模极大值（wavelet modulus maximum，WMM）方法是具有抗噪声干扰能力较强的位场边界检测方法，该方法即使在存有噪声异常干扰的情况下也可对线性构造进行清晰、准确的定位^[3]。

地表观测得到的位场异常通常是地下不同埋深的构造或地质体在观测面产生的综合叠加响应，具有多源性、非线性、非稳态性的特征。在位场数据资料处理过程中，可从位场数据异常中提取研究对象所产生的异常信息^[4-5]。二维经验模态分解（bi-dimensional empirical mode decomposition，BEMD）^[6-7]方法能够自适应地将数据中的各种频率成分有效分离出来，具有多分辨率的特性，其多尺度分解过程可以从高频到低频进行逐次尺度滤波，各尺度分解分量都反映了不同频次的模态特征。文献[8-10]将BEMD方法应用到位场多尺度分解中，并结合径向对数功率谱曲线对分解各尺度异常进行了分析和解释，取得了较好的应用效果。

本文首先通过建立模型重力位场数据，验证WMM方法在存有噪声干扰的情况下对线性构造进行定位和提取的有效性；其次利用BEMD方法对研究地区布格重力异常进行自适应多尺度分解，获得不同深度上构造或岩体线性分布特征；然后充分利用BEMD和WMM方法各自的数据处理优势特征，联合上述两种方法实现对分解各尺度异常进行场源边界检测，获得了不同深度上的场源边界信息；最后结合研究区域地质构造变化特征对各尺度下的边界异常信息进行分析和解释。

1 基本原理

1.1   小波模极大值

文献[1]提出了利用小波模极大值进行边缘检测的基本原理，并在图像领域得到成功应用。假设存在函数$\theta (x,y)$在整个二维平面的积分为1，且x、y在无限远处收敛到零，则称$\theta (x,y)$为平滑函数。那么，平滑函数$\theta (x,y)$在x、y方向上的一阶偏导数分别表示为：

因此，可以将函数${\psi }^x(x,y)$和${\psi }^y(x,y)$看作两个二维小波函数。

在尺度因子$s={\mathrm{2}}^j(j\in Z)$的情况下，对于任意函数$f(x,y)\in L^{\mathrm{2}}(R)$，则由${\psi }^x(x,y)$和${\psi }^y(x,y)$定义的小波变换具有两个分量，分别为：

梯度矢量可表示为：

梯度的模和相角分别表示为：

由二维网格数据结构可知，每个数据点周围有8个相邻数据点，分别用数字1~8标注，如图 1所示。由于梯度方向的对称性，只需考虑1~4区域，此时数据点（x₀，y₀）相角的正切值$\mathrm{t}\mathrm{a}\mathrm{n}{A}_{s}f(x_{\mathrm{0}},y_{\mathrm{0}})$必定落在$\left[\mathrm{1}-\sqrt[]{\mathrm{2}},\sqrt[]{\mathrm{2}}-\mathrm{1}\right)$、$\left[\sqrt[]{\mathrm{2}}-\mathrm{1},\sqrt[]{\mathrm{2}}+\mathrm{1}\right)$、$\left[\sqrt[]{\mathrm{2}}+\mathrm{1},+\mathrm{\infty }\right)$$\left(-\mathrm{\infty },-\mathrm{1}-\sqrt[]{\mathrm{2}}\right]$、$\left(-\mathrm{1}-\sqrt[]{\mathrm{2}},-\sqrt[]{\mathrm{2}}\right]$ 4个区间之一，与位于1~4区域中某个数据的正切值存在一一对应的关系。

图  1  网格数据点离散梯度示意图

Figure  1.  Diagram of the Discrete Gradient of the Grid Data Points

下载: 全尺寸图片 幻灯片

首先计算$\mathrm{t}\mathrm{a}\mathrm{n}{A}_{s}f(x_{\mathrm{0}},y_{\mathrm{0}})$数值的大小，然后判断其落入的区间编号N，将$M_{s}f(x_{\mathrm{0}},y_{\mathrm{0}})$与N区间以及N区间对称梯度方向上相邻两个数据点进行比较，若$M_{s}f(x_{\mathrm{0}},y_{\mathrm{0}})$大于这两个点，则数据点（x₀，y₀）为局部极大值点，否则就不是局部极大值点。因此，小波模局部极大值点刻画了二维数据边缘的特征点，获取到的局部极大值点就是边界数据点。

1.2   二维经验模态分解

文献[11-12]提出了一种自适应性地处理非线性、非平稳信号的分析方法——经验模态分解（empirical mode decomposition，EMD）方法，文献[13-14]将一维EMD拓展应用到二维数据中，利用BEMD方法对二维图像数据进行筛分，得到若干个不同尺度下的二维固有模态函数（bi-dimensional intrinsic mode function，BIMF）和一个剩余分量（Residue，Res）。具体分解流程见图 2。

图  2  BEMD分解流程图

Figure  2.  Flowchart of BEMD Method

下载: 全尺寸图片 幻灯片

设存在一个二维信号$s(x,y)$，其中x=1，2…m，y=1，2…n（m，n为正整数），利用BEMD方法对该信号进行多尺度分解，共得到K个二维固有模态函数BIMF和一个剩余分量Res，即：

2 算例与分析

文献[10]中通过建立球体组合模型试验，利用BEMD方法实现了对组合模型异常的多尺度分解，并结合径向对数功率谱方法实现对分解异常的定性、半定量分析，证明BEMD方法对位场进行多尺度分解是有效的，且分解各尺度异常具有实际地质含义。

为了验证利用WMM方法进行边界检测的有效性，建立了长方体组合模型进行实验，该组合模型是由3个不同形态、不同埋深、不同剩余密度的长方体组成，其空间分布见图 3，相关参数见表 1。组合模型正演重力异常见图 4（a），为模拟组合模型重力异常中混有噪声信号，且噪声信号幅值变化相对较小，因此，向组合模型重力异常中添加信噪比为1.8%的高斯白噪声（见图 4（b）），其中蓝色线框为水平长方体在平面位置投影。

图  3  组合模型空间分布示意图

Figure  3.  Schematic Diagram of the Space Distribution of Combined Model

下载: 全尺寸图片 幻灯片

表  1  组合模型相关参数

Table  1.  Related Parameters of the Combined Models

模型编号	（长, 宽, 高）/m	上顶埋深/m	下底埋深/m	剩余密度/(g·cm-3)
模型1	(600, 300, 300)	150	450	0.8
模型2	(650, 300, 250)	100	350	0.6
模型3	(500, 350, 400)	200	600	1.0

下载: 导出CSV 

| 显示表格

图  4  组合模型正演重力异常

Figure  4.  Gravity Anomalies of Combined Model

下载: 全尺寸图片 幻灯片

本文分别利用总水平梯度倾斜角（the tilt angle of the horizontal gradient，TAHG）^[15]和WMM方法对未加噪声和加有噪声的组合模型异常数据进行处理，图 5为TAHG边界检测结果，图 6为WMM方法边界检测结果。从图 5、图 6中可以看出，对未添加噪声的重力异常，两种方法均可以有效地提取长方体线性边界位置，与长方体实际边界位置能够较好地吻合，TAHG方法计算结果存在边界效应，并没有对内部数据造成较大影响；在添加噪声的情况下，由于TAHG方法属于数值类计算方法，在计算过程中需多次求导运算，高次求导运算会把噪声放大进而对整个数据造成污染，计算结果出现了较大偏差，获取的线性构造分布特征不明显，而WMM方法的计算结果仍能够准确、清晰地展现长方体边界位置，数据中噪声异常对计算结果未产生较大影响，具有较强的抗噪声干扰能力。

图  5  TAHG边界检测结果

Figure  5.  Results of Edge Detection by TAHG Method

下载: 全尺寸图片 幻灯片

图  6  WMM方法边界检测结果

Figure  6.  Results of Edge Detection by WMM Method

下载: 全尺寸图片 幻灯片

3 实际应用

三峡地区西邻青藏高原高山峡谷区域，东临长江中下游平原丘陵区域，区内断裂和褶皱构造发育，地质地貌特征复杂多样。在早期燕山运动中，该地区形成了巫山、大巴山等一系列褶皱山脉，这些山脉由西南-东北走向转为东西走向，地势由南向北逐渐降低。三峡地区断裂构造较发育，断裂带构造控制着该区域的构造演化过程，地质及其相关地球物理背景比较复杂，近几十年以来，该区域受到国内外地质工作者^[16-18]的密切关注。

由于布格重力异常数据在横向上具有较高的分辨率，为了研究三峡地区线性构造变化特征，利用该地区布格重力异常进行相关数据处理、分析与解释，结果见图 7。由图 7可知，研究区域异常整体呈现负异常，异常数值自西向东逐渐增大，存在近北北东向重力变化梯度带，三峡大坝正位于这个梯度带上，重力高值主要集中在东部江汉平原一带，重力低值覆盖了西部整个秭归盆地地区。本文利用BEMD方法提取了研究区域不同深度层次上布格重力异常特征信息，在进行BEMD分解过程中，需要预先设置阈值参数ε判断是否已经筛分出分解分量，ε可设定为较小的数，数值越小，获得的BIMF分量异常越多，但可能导致分解耗时过长。因此，为了在分解过程中获得较多的异常细节信息，设置分解阈值ε = 0.001，运用BEMD方法对三峡地区布格重力异常进行了多尺度分解，共得到5个BIMF和一个Res，分解结果见图 8。

图  7  研究区域布格重力异常

Figure  7.  Bouguer Gravity Anomalies in Study Area

下载: 全尺寸图片 幻灯片

图  8  研究区域布格重力异常的BEMD分解结果

Figure  8.  Decomposition Results of Bouguer Gravity Anomalies by BEMD Method in Study Area

下载: 全尺寸图片 幻灯片

由图 8可知，随着分解阶次增加，不同尺度下的异常呈现不同的形态变化特征，表现为从高频到低频变化。BIMF₁异常呈现零星状或条带状分布，且出现正负伴生的现象，异常幅值在-6.4 ~6.6 mGal之间变化，其主要是地表出露或埋藏较浅的不同密度和岩性地质体产生的；BIMF₂异常多呈现条带状并伴有零星状分布，异常展布范围相对较小，异常幅值在-8.2 ~7.4 mGal之间变化，其反映了浅部岩体或线性构造变化特征；BIMF₃异常幅值在-12.4~12.6 mGal之间变化，在研究区域中部位置存在明显的重力变化高值异常的区域，变化异常走向由北东走向转为北西走向，再转为北东走向，最后变成北西走向，形成一个整体近似菱形重力变化高值异常，异常幅值变化达12.6 mGal，三峡大坝就位于这个高异常条带的中部位置，且三峡大坝附近存有明显的低重力异常区域，在西北部神农架地区和西南部的恩施地区都存在重力高值；BIMF₄异常呈现区域性变化异常，异常变化分布范围较广，在东南部江汉平原地区有条带状正异常，秭归盆地为负异常，西北部神农架地区为高值异常，其为扬子板块出露部分，异常幅值变化较大，反映了区域深部岩性密度或构造特征变化；BIMF₅异常显示研究区域有正负伴生的区域异常，西部秭归盆地地区为负异常，东部江汉平原为正异常，异常幅值变化较大，反映了两个区域地层起伏变化特性；Res异常反映了该区域异常的变化趋势特征，区域异常表现为东高西低，显示了研究区域基底或者莫霍面西低东高的起伏变化情况，在大尺度上体现了该地区岩石圈变化特征。

另外，本文利用径向对数功率谱方法^[19-20]计算了上述各级固有模态函数的功率谱曲线，结果见图 9，其中红色虚线为拟合直线。将图9（a）~ 9（f）中曲线拟合出的直线斜率代入计算公式，得到各固有模态函数对应的场源深度，分别为7.3 km、16.7 km、25.8 km、45.3 km、52.7 km、71.3 km。由此可知，随着分解的增加，各级固有模态函数的波长和其所对应的场源深度逐渐增大。

图  9  BEMD分解各分量的径向对数功率谱曲线

Figure  9.  Radial Logarithmic Power Spectrum Curve of Bi-dimensional Intrinsic Mode Function at All Levels by BEMD Method

下载: 全尺寸图片 幻灯片

BIMF₁、BIMF₂和BIMF₃异常分布比较零散，异常强度和面积相对较小，可能反映了岩体或隐伏岩体、构造、地质异常体的物性和产状变化等；BIMF₄、BIMF₅和Res异常分布较集中，分布范围较广，反映区域内岩性密度变化以及地层或基底界面起伏变化。

因此，将BIMF₁、BIMF₂和BIMF₃相加重构出局部异常，将BIMF₄、BIMF₅和Res相加作为区域异常，重构结果见图 10。图 10（a）中正、负异常呈现零星或区域性变化，部分异常变化特征较明显，主要存在3处高值异常区域，其中位于中部位置多条条带状高值异常整体呈现近菱形状，分布范围较广，异常幅值变化达30 mGal。此外，BIMF₃异常中仍存在一个类似菱形状高值异常，两个异常特征变化存在对应关系，因此该菱形状高值异常应给予关注，可能包含着某类岩性密度变化或地质构造信息。图 10（b）中重构区域异常呈现了西低东高的趋势异常形态，测区西部秭归盆地异常幅值变化最小，东南部的江汉平原地区存有最大幅值变化异常，整体上反映了该研究区域较深结晶基底界面或者地层埋深自西向东逐渐减小的起伏变化情况。

图  10  研究区域布格重力异常BEMD分解重构异常

Figure  10.  Reconstruction Bouguer Anomalies in Study Area by BEMD Method

下载: 全尺寸图片 幻灯片

此外，研究区域存有的主要地质构造是黄陵背斜，在其周围存在一系列的断裂分布，主要形成于印支运动时期，主要包括新华断裂、仙女山断裂、渔阳关-土门断裂、远安断裂、雾渡河断裂和天阳坪断裂等，这些断裂形成一个整体近似菱形状的断裂带，控制着区域的地质构造演化过程，其中黄陵背斜的相对隆升以及当阳、秭归等盆地的相对沉降在某种程度上都受到该断裂带演化过程的影响。文献[21]利用布格重力异常进行三维反演，反演结果显示这些断裂构造的切割深度达到了中地壳上部，属于区域性断裂。

由于BEMD分解得到各级固有模态函数是具有实际地质含义的，反映了不同深度层次上场源特性变化信息，本文对三峡地区布格重力异常进行BEMD分解得到各级BIMF和Res，并利用WMM方法对各尺度异常进行边界检测，检测结果见图 11。

图  11  BEMD分解各级固有模态函数的小波模极大值边界识别结果

Figure  11.  Results of Edge Detection of Bi-dimensional Intrinsic Mode Function at All Levels by WMM Method

下载: 全尺寸图片 幻灯片

由图 11可知，各尺度下边界检测结果与BEMD分解各分量异常特征相关，反映了三峡地区不同深度上的岩体、构造空间展布特征。利用BIMF₃尺度异常获得的线性构造特征中，在中部位置菱形状的线性构造提取得到的部分线性构造的展布位置与菱形断裂带构造（如新华断裂、渔阳关-土门断裂和雾渡河断裂等）存在对应关系，说明BIMF₃异常所反映的场源似深度信息与菱形断裂带构造埋藏深度能够较好的吻合，具有一定实际地质含义。然而，利用其余分量异常获得的线性构造与菱形断裂带构造并不存在明显的对应关系，说明上述不同频率的分量异常是埋藏较浅的地质体或地质构造和埋藏较深的地层起伏或基底变化的响应，菱形断裂带场源异常在上述频段内分解分量异常上反映不明显。

同样地，本文利用WMM方法对BEMD分解重构局部和区域异常进行边界检测，检测结果见图 12。由图 12可以看出，利用重构和区域局部异常获得的线性构造特征中均包含重构各尺度异常识别出的边界信息，是不同深度层次上的场源边界信息的综合响应。

图  12  三峡地区布格重力异常BEMD分解重构异常小波模极大值边界识别结果

Figure  12.  Results of Edge Detection of the Reconstruction Bouguer Anomalies of BEMD by WMM Method

下载: 全尺寸图片 幻灯片

因此，联合BEMD和WMM方法获得了三峡地区不同深度层位上的岩体或构造线性分布特征，实现了位场数据异常多尺度边界检测，且BIMF₃异常边界检测结果与区域地质构造特征能够较好的吻合，这比单纯地分析二分异常对应的边界检测结果具有更高的分辨率。

4 讨论与结论

BEMD方法能够自适应地实现位场数据的多尺度分离，且分解各尺度异常具有实际地质含义。在位场数据中有噪声干扰的情况下，WMM方法仍能够清晰、准确地获取地下地质体的边界信息。本文联合BEMD和WMM方法实现了三峡地区布格重力异常多尺度边界检测，充分发挥了上述两种数据处理方法的优势，对分解各尺度异常对应的不同深度层位上的岩体边界、地质构造线性分布特征进行清晰、准确的刻画，有助于了解地下目标体分布特征，能够对区域构造单元划分、地质填图、矿产资源圈定和结晶基底起伏变化特征等方面研究提供重要参考信息。

此外，BEMD方法在位场数据处理中取得了一定的应用效果，但在拓展应用过程中依然存在着一些缺陷和不足，包括边界效应、模态混叠和分解效率较低等问题，如何提高BEMD方法计算效率的同时又能获得较为准确的分解结果需作进一步研究。