基于欧拉反褶积的基底断裂研究

Research on Basement Fault Using Euler Deconvolution

摘要: 欧拉反褶积法作为无需先验信息的位场自动反演方法，因其计算高效、适用性广等特点，被广泛应用于场源体定位与边界识别。然而，噪声干扰引起的虚假解(如高频噪声使反演结果偏移真实场源)及复杂场源耦合效应(如多尺度异常叠加降低反演精度)等问题，严重制约其在深部构造精细解析中的应用。针对上述不足，文章提出“小波多尺度分析 + 欧拉反褶积联合反演”框架，通过理论建模、对比实验与实例应用验证其技术优势，最终在松辽盆地北部划分出8条主要断裂。为复杂地质区域断裂划分提供新的技术路径。

Abstract: The Euler deconvolution method, as a prior-information-free automatic potential field inversion technique, has been widely applied in source body positioning and boundary identification due to its computational efficiency and broad applicability. However, issues such as spurious solutions caused by noise interference (e.g., high-frequency noise deflecting inversion results from true field sources) and complex source coupling effects (e.g., multi-scale anomaly superposition reducing inversion accuracy) significantly limit its application in detailed deep structural analysis. To address these shortcomings, this paper proposes a “wavelet multiscale analysis + Euler deconvolution joint inversion” framework, verifying its technical advantages through theoretical modeling, comparative experiments, and practical applications, and ultimately delineating 8 major faults in the northern Songliao Basin. This provides a novel technical approach for fault mapping in complex geological regions.

1. 引言

重力数据处理方法的研究一直是地球物理勘探学者研究的热点[1]-[3]。重力具有勘探深度大、数据可信度高、不受电磁干扰的特点，重力异常能反映地下岩体的分布及区域断裂特征[4]。

欧拉反褶积法自提出以来，在国内外得到了广泛的应用与改进。Thompson et al. (1982)最初提出的欧拉反褶积法，通过利用局部重力异常数据反演地下异常源的位置和形态，成为勘探中不可或缺的重要工具[5]。在此基础上，国内外学者不断优化其算法，以提高反演精度。例如，李银飞(2019)通过引进分类算法中的K邻近算法，提高欧拉反演结果的准确性[6]；马国庆等(2020)提出了一种基于欧拉反褶积法的优化算法，重力梯度张量联合欧拉反褶积能够有效抑制噪声对结果的干扰[7]；金文博(2024年)基于相关成像的等效源法能够提高起伏地形下的重磁异常导数计算精度，提高了欧拉反褶积方法的精度[8]。

小波分析的理论基础形成于20世纪70~90年代。1974年法国地质物理学家Morlet首次提出小波概念，通过窗函数伸缩平移构造基函数，成功应用于非稳态地震信号分析[9]。21世纪以来，小波技术拓展至地球物理深度应用，2006年Lelievre P G等人结合连续小波变换与欧拉反演提升场源定位精度[10]。2011年Kaftana等利用离散小波优化磁异常反演，显著提升基底构造识别能力[11]，2000年高德章等人采用二维小波技术解析东海重力异常，成功分离沉积基底与莫霍面信号，验证了方法的地质解释可靠性[12]。2010年宁津生等结合离散与连续小波变换实现干扰源分离与场源深度精准测定[13]。2012年卢海创新性地将小波尺度分析与相关系数成像结合，提升重力场参数反演精度[14]。近年，2018年蒋勇平通过数值模拟与胶东实例验证了小波多尺度法在区域-局部异常分离中的优越性，推动方法走向工程化应用[15]。本文利用小波变换可以提取不同深度信息与良好降噪效果的优势，将小波变换的思想引入欧拉反褶积的计算，为解决欧拉反褶积解的发散问题提供新思路。

2. 基本原理

2.1. 欧拉反褶积

Thompson (1982)提出了该方法的基本理论框架，基于欧拉齐次方程(Euler Homogeneous Relation)来解决这一问题[5]。

考虑一个三维空间中的函数$f\left(x,y,z\right)$ ，如果它满足以下关系：

$\begin{array}{c}f\left(tx,ty,tz\right)=t^{n}f\left(x,y,z\right)\end{array}$ (1)

其中，t是一个常数，n是一个整数，那么该函数$f\left(x,y,z\right)$ 就是一个齐次函数，阶数为n，可以得出该函数满足以下欧拉方程：

$\begin{array}{c}x\frac{\partial f}{\partial x}+y\frac{\partial f}{\partial y}+z\frac{\partial f}{\partial z}=nf\left(x,y,z\right)\end{array}$ (2)

则该方程被称为欧拉齐次方程。

假设在一个观测平面$z=0$ 上存在一个点源，其源点位置为$\left(x_{0,}{y}_0,z_0\right)$ ，该点源在该平面上的磁场异常$f\left(x,y\right)$ 可由以下公式表示：

$\begin{array}{c}f\left(x,y\right)=\frac{c}{\sqrt{{\left[{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2\right]}^N}}\end{array}$ (3)

其中，C是常数，N是构造指数，且C与坐标$x,y,z$ 无关。通过对该异常函数$f\left(x,y\right)$ 进行偏导数计算，并代入欧拉方程(2)，我们可以得到该点源的欧拉方程：

$\begin{array}{c}\left(x-x_0\right)\frac{\partial f}{\partial x}+\left(y-y_0\right)\frac{\partial f}{\partial y}+\left(z-z_0\right)\frac{\partial f}{\partial z}=-Nf\left(x,y,z\right)\end{array}$ (4)

方程(4)表明，通过建立线性方程，我们可以推导出点源的位置$\left(x0,y0,z0\right)$ 。其中N称为构造指数，不同的地质构造对应不同的构造指数，表1列出了构造指数对应的地质构造。

Table 1. Structural index

表1.常见场源的构造指数

构造指数	重力场反应构造
0	岩墙，台阶
0.5	薄板
1.0	岩柱
2.0	球体
3.0	无

2.2. 小波分析

对于信号$x\left(t\right)$ 和母小波函数$\psi \left(t\right)$ ，连续小波变换(Continuous Wavelet Transform, CWT)定义为：

$\begin{array}{c}W\left(a,b\right)=\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}x\left(t\right)\cdot \psi *\left(\frac{t-b}{a}\right)\text{d}t\end{array}$ (5)

式中：$W\left(a,b\right)$ 表示在尺度a和平移b处的变换系数；a为尺度参数(控制小波的伸缩)；b为平移参数(决定小波的位置) $\psi \left(t\right)$ 表示母小波的共轭复数。

将小波伸缩或平移得：

$\begin{array}{c}{\psi }_{a,b}\left(t\right)=\frac{1}{\sqrt{\left|a\right|}}\psi \left(\frac{t-b}{a}\right)\end{array}$ (6)

其中归一化因子$\frac{1}{\sqrt{\left|a\right|}}$ 保证了在不同尺度下小波能量的一致性。

则连续小波变换定义为：

$\begin{array}{c}{W}_f\left(a,b\right)=\frac{1}{\sqrt{\left|a\right|}}\underset{R}{{\displaystyle \int }}f\left(t\right)\overline{\psi }\left(\frac{t-b}{a}\right)\text{d}t\end{array}$ (7)

则对于连续小波变换，$f\left(t\right)$ 可重构变换为：

$\begin{array}{c}f\left(t\right)=\frac{1}{C_{\psi }}\underset{R^2}{{\displaystyle \iint }}\frac{W_{\psi }f\left(a,b\right){\psi }_{a,b}\left(t\right)}{a^2}\text{d}a\text{d}b\end{array}$ (8)

其中常数$C_{\psi }$ 必须满足小波的可逆性条件：

$\begin{array}{c}{C}_{\psi }=\underset{0}{\overset{+\infty }{{\displaystyle \int }}}\frac{{\left|\widehat{\psi }\left(\omega \right)\right|}^2}{\omega }\text{d}\omega <\infty \end{array}$ (9)

且$\widehat{\psi }\left(\omega \right)$ 为母小波的傅里叶变换。

2.3. 小波多尺度分析与欧拉反褶积的联合反演

定义重力异常的多尺度分解形式：

$\begin{array}{c}\Delta g\left(x,y\right)={\displaystyle \sum }_{j=1}^n{D}_j\left(x,y\right)+A_n\left(x,y\right)\end{array}$ (10)

其中，D_j为第j层细节信号(高频局部异常)，A_n为第n层逼近信号(低频区域场)。利用小波分解提取不同尺度的异常分量(D_j)，分别输入欧拉反褶积进行反演，避免场源耦合干扰。

构造指数的多尺度优化：

高频分量(浅部异常)：优先选择N = 0(接触带/岩墙)；

低频分量(深部异常)：优先选择N = 2 (球体/基底断裂)。

$\begin{array}{c}{N}_j=f\left({\lambda }_j,{\sigma }_j\right)\end{array}$ (11)

其中，λ_j为尺度波长，σ_j为噪声水平。

根据小波分解尺度动态调整欧拉反褶积的窗口大小：

$\begin{array}{c}{W}_j=k\cdot {\lambda }_j\left(k=2\sim 3\right)\end{array}$ (12)

例如，浅层高频异常(j = 2)对应小窗口(W₂ = 5 × 5)，深层低频异常(j = 4)对应大窗口(W₄ = 15 × 15)。

多尺度联合反演流程：

Algorithm1. Step inversion algorithm

算法1. 分步反演算法

输入：原始重力异常数据Δg 1. 小波分解：Δg → {D1, D2, D3, A3} %A3代表三层分解，根据需求选择层数 2. 对每层Dj进行阈值去噪→${D^{\prime }}_j$ 3. 对${D^{\prime }}_j$ 和A3分别进行欧拉反褶积，参数Nj和Wj按2.3节动态选择 4. 融合各尺度反演结果，剔除虚假解 输出：场源位置(x0, y0, z0)

2.4. 小波去噪模型的建立

基于信号处理理论，受污染观测信号可构建如下数学模型：

$\begin{array}{c}s \left(n\right)＝f\left(n\right)+\sigma \text{e}\left(n\right)\end{array}$ (13)

式中：$\text{e}\left(n\right)$ 表征随机干扰项；σ为噪声水平系数。在重力场数据应用中，目标信号$f\left(n\right)$ 特指重力场异常信号。当$\epsilon \left(n\right)$ 服从高斯白噪声过程时，建立的小波去噪模型可有效实现重力异常信号g的优化提取。该模型在统计学框架下可解释为基于正交基函数的非参数回归过程。本方法通过三级处理流程实现信号增强：

a) 多尺度分解：选择适配的小波基函数对含噪信号实施多分辨率分析；

b) 系数优化处理：构建阈值函数对细节系数进行非线性量化；

c) 信号重构：通过逆变换合成降噪后的时域信号。

优化后的信号需满足双重要求：

a) 重构信号与原始信号具有同阶连续可微性；

b) 误差最小化：重构信号满足$\min \left\{E\left[\left(f\left(n\right)-f^\left(n\right)\right)2\right]\right\}$ 约束条件其中阈值选择策略与量化方式共同决定最终降噪性能，需通过参数优化实现信噪比最大化。

其中信噪比(SNR)是用于衡量信号质量的一个重要指标，他表示信号与噪声的相对强度或功率之间的比值[16]。信噪比越高，表示信号相对于噪声更强，质量更好。在计算信噪比时，它可以通过以下公式计算：

$\begin{array}{c}\text{SNR}=10\ast log10\left(\frac{Ps}{Pn}\right)\end{array}$ (14)

式中：$Ps$ 为有效信号功率均值，$Pn$ 为噪声功率统计量。该对数化处理符合信号处理领域对功率比的标准化表达规范，分贝单位制可有效压缩动态范围。

3. 模型试验

3.1. 单体模型实验

长方体是一种非常重要的地质模型，其不同的形态可以模拟不同的地质体。当长宽高相同时，可以模拟等轴状地质体；当厚度较薄而水平展布面积较大时，可以模拟岩床；而当长宽较小高度较大时，可以模拟岩柱、岩脉及岩墙等。建立三维直角坐标系，向东1 km为x轴，向南1 km为y轴，垂直地面向下为z轴。下面的模型数据都是基于此坐标系，具体参数如表2所示，空间位置如图1所示，正演结果如图2所示。

以模型重力场及其x、y、z方向导数为基础数据(如图2)，建立欧拉方程，得到反演结果如图3所示。可以看出，欧拉解很好地反映地质体投影所在的位置，说明欧拉反褶积方法能很好地估算出长方体边界。

Table 2. Cuboid model information

表2. 长方体模型参数

	正方体
相对密度(kg/m^3)	0.4
中心位置	(500, 500, 100)
X方向长度(m)	100
Y方向长度(m)	100
Z方向长度(m)	100

(a)正方体空间位置图 (b)正方体水平位置图

Figure 1. Map of Cuboid model location

图1. 正方体位置图

(a)理论重力异常 (b) x方向一阶导

(c) y方向一阶导 (d) z方向一阶导

Figure 2. Map of upward continuation

图2.正方体正演图

Figure 3. Map of Euler deconvolution

图3. 正方体理论重力欧拉解

在实际地球物理数据采集过程中，由于仪器精度、环境干扰等因素，噪声不可避免地会混入初始信号里。这些噪声可能会导致重力异常数据的失真，进而影响欧拉反褶积的性能。为了模拟真实环境并测试方法的适用性，在正演的结果中添加高斯噪声，并在此基础上应用小波多尺度分析，以期减弱噪声干扰并提升误差估计的准确性。在正方体模型的重力异常上加入0.05的高斯噪声，经过计算加入0.05高斯噪声之后的信噪比为16.2931，正演结果和欧拉解如图4所示。可以看出由于噪声的影响，导致解的聚集性很差，解集极为发散，又由于对噪声各方向导数的计算会加大这些干扰，严重影响欧拉反褶积对地质体的边界圈定。

(a) 正方体加噪正演 (b) 正方体加噪欧拉解

Figure 4. Map of Euler deconvolution for added noise

图4. 添加噪声欧拉解

在小波多尺度分析中，分解层数的选择对去噪效果至关重要。水平圆柱体重力异常信号具有显著的低频特性，添加的0.05高斯噪声主要表现为高频分量。考虑到实验数据网格为51 × 51，分解层数需满足$N\ge {\log }_2\left(信号主波长/噪声最小波长\right)$ ，根据模型参数(埋深$h=500\text{ m}$ ，半径$r=50\text{ m}$ )，信号主波长$\lambda \propto h$ ，噪声波长$\lambda \text{noise}\propto 网格间距\left(20\text{ m}\right)$ ，理论层数$N\approx 3\sim 5$ 。最大可分解层数为5层，因此我们选择将信号分解至第4层。这一层数能够在有效去除高频噪声的同时，保留信号的整体形态和重要细节。下面将小波去噪后的水平圆柱体进行欧拉反演，结果如图5所示。

Figure 5. Map of Euler deconvolution for denoised

图5. 小波去噪欧拉解

从图中可以看出，小波多尺度分析对噪声的去除有明显的效果。另外，信噪比也提升到了31.1876 dB，原本离散的虚假解基本上都能够有效地去除。对于正方体来说，小波逼近后的欧拉解对边界也能够很好地识别出来，深度方面也有很好的反映，整体来说使用小波多尺度分析去噪后的反演结果还是较为准确的，尤其是水平位置。

3.2. 组合模型试验

Table 3. Model data

表3. 球体组合模型参数

序号	中心点位置(x0, y0, z0)	半径(m)	剩余密度(g/cm^3)
1	(400, 800, 100)	70	0.5
2	(600, 700, 200)	120	0.5
3	(800, 600, 100)	80	0.5
4	(200, 300, 150)	80	0.5
5	(500, 200, 300)	65	0.5

(a)空间位置图 (b)水平位置图

Figure 6. Map of balls model location

图6. 球体组合位置图

在实际情况中，更多的是多个场源叠加在一起，一个计算窗口内可能包含多个异常。为了试验小波多尺度分析对于有其他场源干扰时欧拉反褶积的应用效果，本文进行了球体组合模型试验。具体参数见表3球体组合参数，具体位置如图6球体组合位置图所示。

在球体模型的重力异常加入0.01的高斯噪声，经过计算得到其信噪比为23.2779，加入噪声的正演结果以及欧拉解如图7所示。

(a) 球体组合加噪正演 (b) 球体组合加噪欧拉解

Figure 7. Map of Euler deconvolution for added noise

图7. 添加噪声欧拉解

从图7可以看出，欧拉解可以反映出模型大概的位置，但是缺乏很多细节方面的信息，同时存在很多虚假解，无法直接圈定出他们的边界。从埋深程度来看，识别结果比实际埋深较浅，与实际情况差距较大。综合来看，欧拉反褶积在噪声的干扰下，可以识别出埋藏较浅，体积较大的地质体，但是对埋藏较深，体积较小的地质体识别效果较差，说明欧拉反褶积受噪声干扰较大，抗干扰能力弱。

对添加了噪声的异常进行小波多尺度分析，实验数据网格为51 × 51，将结果分解到四层，分解之后的信噪比为29.2985 dB，对比加入0.01高斯噪声之后提升了6.0206 dB。得到的欧拉解如图8所示，从图中我们可以看出，经过小波多尺度分析的处理所得到的结果很好地将噪声所产生的异常进行了去除，对于模型1、2、3来说欧拉解很好地将其三个共同的边界给识别了出来，但是对于个体所给的信息较少，不能很好地将其识别出来。模型4的识别效果也很好，反演深度来看也接近实际埋深。而对于模型5，并没有得到很好的场源位置和深度信息，推测是由于模型5埋深较深，体积较小，从而在小波逼近的过程中，被当作区域异常分解掉。

Figure 8. Map of Euler deconvolution for denoised

图8. 小波去噪欧拉解

4. 实际资料应用

Figure 9. Bouguer gravity anomaly map of the research area

图9. 研究区布格重力异常图

本研究基于1:20万松辽盆地北部重力异常数据，采用Kriging插值法进行网格化处理，以保证数据精度及异常趋势的合理性，如图9所示。研究区布格重力异常整体分布在−400 g.u.至160 g.u.之间，呈现出明显的分区特征，主要受区域构造特征控制。

4.1. 小波多尺度分析

采用二维离散小波分解技术对松辽盆地北部地区重力异常进行分解，将信号分解为低频和高频部分。母函数采用coif4，分解层数为5。每一层的分解都会产生一个更低分辨率的逼近和更多的细节分量信息，由于一阶小波分解只对布格重力异常进行了一次重构，故本文不做一阶小波分解。研究区各阶小波逼近重力异常等值线如图10所示。

(a) 二阶小波逼近 (b) 三阶小波逼近

(c) 四阶小波逼近 (d) 五阶小波逼近

Figure 10. Maps of wavelet approaching

图10. 各阶小波逼近图

从各阶小波的逼近图我们可以看出：随分解尺度增大，重力异常等值线渐趋平滑。二阶分解结果与布格重力异常空间分布相似，幅值降低约30%，浅源地质体仍贡献40%场源特征，哈尔滨–大庆地区保留5~8个次级异常。三阶分解时场源分离效应增强，北安南缘等浅源异常区分辨率下降60%，仅余3~4个主异常中心，区域构造初显但松辽盆地中央凹陷带梯度特征仍显著。四阶分解基本滤除浅部干扰，齐齐哈尔基底隆起与哈尔滨幔坡带对应性增强。五阶分解虽细节损失严重，区域构造NNE向展布保持稳定，揭示松辽盆地北缘NNE向基底断裂体系控制白垩系沉积中心迁移及火山岩带分布。

如图11是小波多尺度分解的细节图，从图中我们可以看出：二阶细节图(a)在北向、西南及东南向显示多个浅源低/高异常体，反映浅层小范围密度非均质性。四阶细节图(c)中齐齐哈尔–大庆–哈尔滨一带呈现北北东向重力异常带，与研究区整体重力场方向高度一致。五阶细节图(d)显示局部异常结构复杂化，多处高值响应区可能指示基底岩性界面或深部构造叠加。多阶次对比解析表明，区域地质体密度呈分层非均质分布，齐齐哈尔至哈尔滨基底埋藏深度呈现“浅–深–浅”波状变化，反映构造作用非均匀性。

(a) 二阶小波细节 (b) 三阶小波细节

(c) 四阶小波细节 (d) 五阶小波细节

Figure 11. Maps of wavelet details

图11. 各阶小波细节图

4.2. 欧拉反演

本研究聚焦于场源定位领域的欧拉反演技术，通过系统对比小波欧拉解算方法与常规布格重力异常欧拉反演算法的差异，揭示出联合小波变换的欧拉反演方程组在噪声抑制方面的显著优势。模型实验表明，该方法通过多尺度分解有效提升了欧拉反演的精确性，进而实现了对场源空间分布的精准定位，为发展多方法协同反演技术提供了重要理论支撑，对研究区进行小波多尺度分析与欧拉反褶积的联合反演，得到的欧拉解并将其以散点图的形式投影到小波四阶细节重力异常，如图12所示。

Figure 12. Maps of Euler deconvolution

图12. 研究区欧拉解(底图为四阶小波细节)

基于松辽盆地多源勘探数据，本研究开展重力场特征解析。通过小波分解技术提取异常细节信息，成功识别出研究区深部断裂构造特征。进一步对比小波域反演结果与布格重力异常欧拉解褶成果，验证了前者在解析复杂地质体边界及构造形态方面的可靠性。该研究不仅深化了对区域重力场特性的认识，更为联合反演技术的应用推广提供了实证依据。

Figure 13. Fracture tectonic delineation of research area

图13. 研究区断裂构造划分图

本文根据欧拉反演结果的聚集程度推断测区的断裂分布，并与常规方法划分断裂的结果相比较。为断裂构造的精细化研究提供了新方法。最终结合地质资料综合对比分析研究区布格重力异常特征，研究区区域重力异常特征以及不同深度的欧拉反褶积反演。可以较为精确地划分出研究区内深大断裂的分布和走向，经过综合分析对比本文为松辽盆地北部地区划分出主要断裂8条，分别编号为F1至F8，如图13所示。

5. 结论

本研究通过欧拉反褶积算法并融合小波多尺度分析技术，构建了适用于地质体反演的“小波多尺度分析–欧拉反褶积联合反演”框架，系统验证了方法在噪声抑制中的有效性，并成功应用于松辽盆地北部基底断裂的划分。得到以下主要结论：

1) 正演模型：构建了长方体正演模型，模拟地下地质体重力。使用“小波多尺度分析–欧拉反褶积联合反演”可以显著抑制高频噪声干扰，解决了传统欧拉反褶积解的发散性问题。

2) 松辽盆地北部深部构造解析：结合小波多尺度分析与欧拉反褶积算法，成功反演松辽盆地北部断裂分布，划分出8条主要断裂。

断裂的划分作为地质研究的基础框架，其建立需要多学科资料的交叉验证，包括地质填图、地球物理探测、地球化学分析、钻井岩心等多源数据的系统集成。受限于个人学术积累与研究经验，在本研究中，对断裂样式的深度、边界条件确定及演化序列重建等方面，可能存在局限性或认知偏差。客观而言，要实现更精准合理的断裂划分，尚需依托后续更系统的多学科数据支撑和更深入的区域构造演化研究。

NOTES

^*通讯作者。