改进Tilt-Euler反褶积法在重力数据反演中的适用性研究

Study on the Applicability of Improved Tilt-Euler Deconvolution Method in Gravity Data Inversion

摘要: 为解决传统欧拉反褶积法中构造指数选择困难和反演解发散的问题，本文提出并应用了一种改进的倾斜角欧拉(iTilt-Euler)反褶积方法。通过与传统欧拉法和倾斜角欧拉(Tilt-Euler)法的对比分析，本文评估了改进方法在不同埋深、多个场源体以及加噪声情况下的适用性。在单个立方体模型试算中，改进Tilt-Euler法能够准确识别浅部埋深地质体的边界，随着埋深增加，边界识别略有偏差。在组合模型试算中，该方法能够有效识别多个场源体及其接触区域的边界。通过加入2%的高斯噪声进行的试算显示，改进Tilt-Euler法能够有效抑制噪声对反演结果的影响，提高了识别精度。结果表明，改进Tilt-Euler法在处理复杂地质体和噪声干扰时表现出较传统方法更强的稳定性和准确性，具有较好的实际应用前景。

Abstract: To solve the problems of difficult index selection and divergent inversion solutions in traditional Euler deconvolution methods, this paper proposes and applies an improved tilt angle Euler (iTilt-Euler) deconvolution method. By comparing and analyzing with the traditional Euler method and Tilt-Euler method, this paper evaluates the applicability of the improved method under different burial depths, multiple field sources, and noisy conditions. In the trial calculation of a single cube model, the improved Tilt-Euler method can accurately identify the boundaries of shallow buried geological bodies. As the burial depth increases, there is a slight deviation in boundary identification. In the combination model calculation, this method can effectively identify the boundaries of multiple field source bodies and their contact areas. The trial calculation by adding 2% Gaussian noise shows that the improved Tilt-Euler method can effectively suppress the influence of noise on the inversion results and improve the recognition accuracy. The results show that the improved Tilt-Euler method exhibits stronger stability and accuracy than traditional methods when dealing with complex geological bodies and noise interference, and has good practical application prospects.

1. 引言

随着边缘识别技术的进步，准确推断地质体的边界和深度已成为位场数据处理的关键。重力和磁场数据在识别地质体边界方面具有独特优势。欧拉反褶积法作为一种有效的位场识别方法，能够快速识别场源体，在重磁数据解释中展现了广阔的应用前景[1]。然而，该方法仍面临构造指数选择困难和反演解发散等问题，影响了其实际应用。与传统欧拉法不同，倾斜角欧拉(Tilt-Euler)法无需构造指数，能更有效地反演地质体的边界和深度。但在计算倾斜角梯度的水平导数时，存在“解析奇点”，导致结果不稳定。近年来，针对这一问题的改进倾斜角欧拉(iTilt-Euler)法取得了显著进展，逐渐成为位场反演方法研究的热点。

欧拉反褶积方法最初由Peters提出，其核心思想是通过欧拉齐次方程中的构造指数来反演磁测数据的场源深度和位置[2]。Thompson将该方法扩展至平面，结合计算机辅助技术，通过磁性剖面数据估算自然场源的深度[3]。随后，英国物理学家Reid将其推广到三维空间，推动了该方法的发展[4]。近年来，针对构造指数的选择，研究者们通过改进计算技术加以优化。姚长利等采用变化构造指数的扫描技术，提高了场源解的准确性[5]。曹书锦等则通过核密度估算和误差分析，确定了最佳构造指数[6]。马国庆等提出的梯度欧拉反褶积法结合了场源参数和构造指数，取得了更好的反演效果[7]。周文月研究了小窗口多资料的影响，有效解决了场源混叠干扰带来的识别问题[8]。随着技术的进步，越来越多的地球物理学者将欧拉反褶积方法应用于实际勘探中。马国庆等扩展该方法至张量数据，显著提高了反演结果的分辨率并降低了噪声[9]。周文月等则在文顿丘盐区应用全张量航空重力欧拉反褶积法，成功划分了地下异常区[10]。

为解决欧拉反褶积法中构造指数选择困难和反演解发散等问题，本文采用了不依赖构造指数的Tilt-Euler法和改进Tilt-Euler法。通过构建单一立方体和组合立方体的正演模型，进行欧拉反褶积计算，成功克服了由于构造指数选择不当导致的反演结果发散问题，并提高了反演效果。同时，采用改进Tilt-Euler法，通过正反演模型试验，并与传统方法计算结果进行对比，分析改进Tilt-Euler反褶积法的适用性。本研究使用MATLAB软件进行计算和分析。

2. 方法原理

2.1. 欧拉反褶积基本原理

欧拉方程描述了重力异常数据的几何参数与场源之间的关系。通过求解欧拉齐次方程，可以推算场源体的水平位置和深度，从而实现对断层、台阶等地质构造的识别与估算[11] [12]。

欧拉反褶积的基本原理是：在笛卡尔坐标系中，选取z轴向下为正，假设存在一个重力异常源$O\left(x_0,y_0,z_0\right)$ ，在观测平面上任意一点$Q\left(x,y,z\right)$ 处的重力异常形式可写成：

$f\left(x,y,z\right)=\frac{c}{{\left(\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2+{\left(z-z_0\right)}^2}\right)}^N}$ (1)

式中$O\left(x_0,y_0,z_0\right)$ 为场源中心位置坐标；$Q\left(x,y,z\right)$ 为观测点的坐标；为场源到观测点Q之间的距离[13]；c是一个常数；$N=1,2,3,\cdots$ ，是构造指数，与场源的几何构造、构造地质体的特性有关[14]。实际地下地质构造情况复杂，选取合适的构造指数尤为重要[15]。常见的地质构造体所对应的构造指数如表1所示。

Table 1. Structural indices selected for different geological structures

表1. 不同地质构造所选取的构造指数

构造指数	磁场反映的构造	重力场反映的构造
0	无限薄板盖层	岩墙、岩阶
0.5	垂直接触带	薄板
1.0	断层、水平窄岩脉	直立岩脉
2.0	直立岩脉、水平宽岩脉	球体
3.0	均匀质量球体	无

而对于任意环境的地形，欧拉齐次方程的一般表达式为：

$\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}=-N\left[f\left(x-x_0,y-y_0,z-z_0\right)\right]$ (2)

即方程是N阶其次的欧拉方程；为重力异常在$x,y,z$ 三个方向上的偏导数；为物性均匀的场源体的重力异常值。

当研究区存在背景干扰时，在上述方程中引入了一个未知常数B来消除背景场异常的影响，其表达式可改写为：

$\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}=N\left[B-f\left(x-x_0,y-y_0,z-z_0\right)\right]$ (3)

欧拉方程的求解过程依赖于场源构造指数的确定，而这一构造指数并非唯一且往往难以精确确定。在实际应用中，构造指数的选择直接影响反演结果的精度，尤其是在复杂地质条件下。如果构造指数选取不当，可能导致异常位置的推断出现较大的误差，甚至完全偏离实际地质体的边界。这种不准确的推断在处理具有多重构造和复杂地质体的区域时尤为明显，错误的构造指数可能会使反演结果产生显著偏差，从而影响后续的地质解释和资源评估。因此，如何改进欧拉反演解的收敛性，确保其在各种地质条件下都能稳定有效地收敛到正确的解，是该方法在地质勘探中的推广应用面临的关键问题之一。

2.2. 不依赖于构造指数的欧拉方程

由于此类方程仅含未知场源坐标，不含构造指数N，所以避免了由于构造指数的选择所导致的计算误差，相较于传统欧拉方程，其反演解具有更好的稳定性和收敛性。

2.2.1. Tilt-Euler法

Miller和Singh将解析信号相位的概念引入了边界识别[16]，称倾斜角，又称为斜导数。倾斜角的计算公式为：

$Tilt={\tan }^{-1}\left(\frac{VDR}{THDR}\right)$ (4)

倾斜角方法通过分析场源体的倾斜角变化特性来识别不同埋深的地质体边界。倾斜角在接近地质体边界时会发生明显变化，而在深部区域，地质体的影响相对较小，能够更准确地反演浅至深部的边界信息。该方法的优势在于无需依赖于复杂的构造指数，便能有效识别和定位场源体边界。Salem等人基于倾斜角公式进一步推导了Tilt-Euler法，该方法在传统欧拉反褶积法的基础上进行了改进，消除了对场源构造指数的依赖[17]。通过这一创新，Tilt-Euler法能够更精确地估计场源位置，并在复杂地质环境下提高了反演结果的稳定性和可靠性，极大地提升了欧拉反褶积法在实际应用中的效果和适用性。Tilt-Euler方程为：

$k_{x}x+k_{y}y+k_{z}z=k_x{x}_0+k_y{x}_0+k_z{x}_0$ (5)

其中$k_x,k_y,k_z$ 分别是倾斜角梯度沿$x,y,z$ 方向的导数：

$k_x=\frac{\partial Tilt}{\partial x}=\frac{1}{A^2}\left[\frac{\partial f}{\partial h}\frac{{\partial }^{2}f}{\partial x\partial z}-\left(\left(\frac{\partial f}{\partial z}\right)/\left(\frac{\partial f}{\partial h}\right)\right)\cdot \left(\frac{\partial f}{\partial x}\frac{{\partial }^{2}f}{\partial x^2}+\frac{\partial f}{\partial y}\frac{{\partial }^{2}f}{\partial y\partial x}\right)\right]$ (6)

$k_y=\frac{\partial Tilt}{\partial y}=\frac{1}{A^2}\left[\frac{\partial f}{\partial h}\frac{{\partial }^{2}f}{\partial y\partial z}-\left(\left(\frac{\partial f}{\partial z}\right)/\left(\frac{\partial f}{\partial h}\right)\right)\cdot \left(\frac{\partial f}{\partial x}\frac{{\partial }^{2}f}{\partial x\partial y}+\frac{\partial f}{\partial y}\frac{{\partial }^{2}f}{\partial y^2}\right)\right]$ (7)

$k_z=\frac{\partial Tilt}{\partial z}=\frac{1}{A^2}\left[\frac{\partial f}{\partial h}\frac{{\partial }^{2}f}{\partial z^2}-\left(\left(\frac{\partial f}{\partial z}\right)/\left(\frac{\partial f}{\partial h}\right)\right)\cdot \left(\frac{\partial f}{\partial x}\frac{{\partial }^{2}f}{\partial x\partial z}+\frac{\partial f}{\partial y}\frac{{\partial }^{2}f}{\partial y\partial z}\right)\right]$ (8)

其中$A=\sqrt{{\left(\frac{\partial f}{\partial x}\right)}^{\text{2}}+{\left(\frac{\partial f}{\partial y}\right)}^{\text{2}}+{\left(\frac{\partial f}{\partial z}\right)}^{\text{2}}}$ 为该场的解析信号振幅。公式(5)与传统欧拉反褶积方程类似，但不包含构造指数。因此，该方法克服了传统欧拉反褶积法由于构造指数选择不当而引起解的发散性，从而为实际数据的处理和解释提供了一种行之有效的途径。

2.2.2. 改进Tilt-Euler法

在区域背景场下，位于倾斜角分母上的总水平导数(THDR)可能趋近于零，使得倾斜角的计算出现“解析奇点”[18]。因此，刘鹏飞等对倾斜角进行了改进[19]，其数学定义式为：

$iTilt={\tan }^{-1}\left(\frac{\partial f/\partial z}{\sqrt{{\left(\partial f/\partial x\right)}^2+{\left(\partial f/\partial y\right)}^2+{\left(\partial f/\partial z\right)}^2}}\right)$ (9)

在该式中，解析信号振幅被用来替代分母中的总水平导数模，这一改进有效提高了计算结果的稳定性。同时，改进后的倾斜角方法充分发挥了其在识别微弱异常方面的优势，能够更敏锐地捕捉到不同埋藏深度的隐伏场源信息。在水平总梯度趋于零时，$k_x,k_y,k_z$ 的结果中存在“解析奇点”。针对这一问题，Huang等基于改进倾斜角(即公式9)推导了其欧拉反褶积形式[20]，记为iTilt-Euler，表达式如下：

${k^{\prime }}_x{x}_0+{k^{\prime }}_y{y}_0+{k^{\prime }}_z{z}_0={k^{\prime }}_{x}x+{k^{\prime }}_{y}y+{k^{\prime }}_{z}z$ (10)

其中${k^{\prime }}_x,{k^{\prime }}_y,{k^{\prime }}_z$ 分别为改进倾斜角沿$x,y,z$ 方向的导数，表示为：

${k^{\prime }}_x=\frac{\partial iTili}{\partial x}=\frac{1}{A^2}\left[A\frac{{\partial }^{2}f}{\partial x\partial z}-\left(\left(\frac{\partial f}{\partial z}\right)/A\right)\cdot \left(\frac{\partial f}{\partial x}\frac{{\partial }^{2}f}{\partial x^2}+\frac{\partial f}{\partial y}\frac{{\partial }^{2}f}{\partial y\partial x}+\frac{\partial f}{\partial z}\frac{{\partial }^{2}f}{\partial z\partial x}\right)\right]$ (11)

${k^{\prime }}_y=\frac{\partial iTili}{\partial y}=\frac{1}{A^2}\left[A\frac{{\partial }^{2}f}{\partial y\partial z}-\left(\left(\frac{\partial f}{\partial z}\right)/A\right)\cdot \left(\frac{\partial f}{\partial x}\frac{{\partial }^{2}f}{\partial x\partial y}+\frac{\partial f}{\partial y}\frac{{\partial }^{2}f}{\partial y^2}+\frac{\partial f}{\partial z}\frac{{\partial }^{2}f}{\partial z\partial y}\right)\right]$ (12)

${k^{\prime }}_z=\frac{\partial iTili}{\partial z}=\frac{1}{A^2}\left[A\frac{{\partial }^{2}f}{\partial z^2}-\left(\left(\frac{\partial f}{\partial z}\right)/A\right)\cdot \left(\frac{\partial f}{\partial x}\frac{{\partial }^{2}f}{\partial x\partial z}+\frac{\partial f}{\partial y}\frac{{\partial }^{2}f}{\partial y\partial z}+\frac{\partial f}{\partial z}\frac{{\partial }^{2}f}{\partial z^2}\right)\right]$ (13)

其中$A^{\prime }=\sqrt{{\left(\frac{\partial f}{\partial z}\right)}^2+A^2}$ 。在求解改进倾斜角各方向导数的过程中，仅有解析信号振幅A位于分母，即使在A趋于0的情况下，其各方向导数公式(11)~(13)中均不存在“解析奇点”。因此，仍然可以计算${k^{\prime }}_x,{k^{\prime }}_y,{k^{\prime }}_z$ 的值。改进Tilt-Euler既保留了Tilt-Euler法不依赖于构造指数的优点，又能较好地提高反演结果的稳定性。

3. 模型试算

对理论模型进行正反演计算是地球物理勘探中的一种重要手段，它能够反映源与场分布特征之间的联系，对于认识问题、深化解释工作具有十分重要的作用。本文设计了不同埋深的单一立方体模型，并使用改进的Tilt-Euler法进行处理，分析该方法在不同埋深条件下对地质体的识别效果。同时，考虑到实际情况中场源通常会叠加，本文还进行了组合模型试验，测试了欧拉反褶积法在存在其他场源干扰时的应用表现。通过对比传统欧拉法、Tilt-Euler法和改进Tilt-Euler法在处理加噪声和未加噪声的重力异常数据时的效果，重点分析了改进Tilt-Euler法在场源干扰情况下的识别性能。

3.1. 单个模型试算

3.1.1. 单个模型改进Tilt-Euler法反演结果

为了检验改进Tilt-Euler法在位场数据处理中的实际应用效果，本文选取了顶部埋深为2 km的单个立方体模型进行理论计算。图1(a)为单个立方体模型实际平面位置，设置网格间距为1 km × 1 km，剩余密度为0.3 g/cm³。

从重力异常结果(图1(b))可以看出，立方体模型的重力异常范围在0至35 mGal之间。通过改进Tilt方法的计算结果(图1(c))，可以发现该方法通过零值准确反映了边界位置，极大值对应于模型体中心，且所得到的模型大小与实际情况相符。改进的Tilt方法相较于重力异常图，能够更精确地收敛模型边界，同时有效避免了背景场区域高幅值的畸变现象，确保了计算过程的稳定性。从改进Tilt-Euler法的反演结果(图1(d))来看，反演效果整体较好，精度较高，能够准确识别模型的大小及边界位置。反演得到的深度为2.1 km，与实际深度非常接近。

3.1.2. 不同埋深计算结果

为了探讨改进Tilt-Euler法对于不同地质体埋深的反演效果，本文选取6个厚度和密度相同，但顶面埋深不同的单个立方体进行理论模型计算。设置100 km × 100 km的网格区域，模型1、2、3、4、5、6的顶面埋深分别为2、4、6、8、10、12 km，剩余密度均为0.3 g/cm³，具体参数设置见表2。

Figure 1. (a) Plan position diagram; (b) Gravity anomaly results; (c) Improve Tilt calculation results; (d) Improve the inversion results of Tilt-Euler method

图1. (a) 平面位置图；(b) 重力异常结果；(c) 改进Tilt计算结果；(d) 改进Tilt-Euler法反演结果

Table 2. Single cube model data

表2. 单个立方体模型数据

模型编号	顶面埋深 (km)	厚度 (km)	x轴范围 (km)	y轴范围 (km)	剩余密度 (g/cm^3)
1	2	3	(−30, 30)	(−30, 30)	0.3
2	4	3	(−30, 30)	(−30, 30)	0.3
3	6	3	(−30, 30)	(−30, 30)	0.3
4	8	3	(−30, 30)	(−30, 30)	0.3
5	10	3	(−30, 30)	(−30, 30)	0.3
6	12	3	(−30, 30)	(−30, 30)	0.3

从改进Tilt结果(图2)可以看出，该方法通过零值准确反映了模型的边界位置，极大值对应于模型体的中心，且能够增强埋深较深区域的弱异常信息。对于埋深较浅(2~6 km)的模型，所得到的模型大小较为准确。从图2(a)~(f)中的结果可以观察到，单一立方体模型的幅度异常范围在−0.8至0.8弧度(Radians)之间。随着埋深的增加，异常幅值逐渐偏向场源的边界外侧，且其位置略大于实际深度，并且出现了轻微的扭曲。总体而言，改进Tilt方法有效避免了背景场区域产生高幅值的畸变现象，确保了计算结果的稳定性。

Figure 2. Improved Tilt calculation results (a) Top burial depth of 2 km; (b) Top burial depth of 4 km; (c) Top burial depth of 6 km; (d) Top burial depth of 8 km; (e) Top burial depth of 10 km; (f) Top burial depth of 12 km

图2. 改进Tilt计算结果图(a) 顶部埋深2 km；(b) 顶部埋深4 km；(c) 顶部埋深6 km；(d) 顶部埋深8 km；(e) 顶部埋深10 km；(f) 顶部埋深12 km

以改进Tilt为基础数据，结合欧拉方程进行求解，设定构造指数N = 0，滑动窗口大小为11 × 11，误差范围限定为15%。反演结果如图3(a)~(f)所示。从图3可以看出，该方法整体效果良好，能够准确识别模型的大小和边界位置。对于实际埋深分别为2、4、6、8、10和12 km的立方体模型，改进Tilt-Euler法识别的平均深度分别为2.2、4.5、6.9、9.1、11.6和14.1 km。可以看出，该方法在识别埋深较浅的地质体时表现尤为准确。

Figure 3. Improved Tilt-Euler method inversion results (a) Top burial depth of 2 km; (b) Top burial depth of 4 km; (c) Top burial depth of 6 km; (d) Top burial depth of 8 km; (e) Top burial depth of 10 km; (f) Top burial depth of 12 km

图3. 改进Tilt-Euler法反演结果图(a) 顶部埋深2 km；(b) 顶部埋深4 km；(c) 顶部埋深6 km；(d) 顶部埋深8 km；(e) 顶部埋深10 km；(f) 顶部埋深12 km

3.1.3. 不同数据类型反演结果对比

为了探讨改进Tilt-Euler法相较于传统欧拉法和Tilt-Euler法的优势，本文选择顶部埋深为12 km的单个立方体模型，使用不同的欧拉反演方法进行计算，结果如图4所示。从单个立方体的重力异常图(图4(a))可以看出，模型中心的重力强度最大，约在模型的中心位置，随着距离中心的增加，重力值逐渐减小，虽然可以看到一些边界特征，但具体的边界位置仍不够清晰。

本文选取传统欧拉法、Tilt-Euler法和改进Tilt-Euler法进行欧拉反褶积反演，结果如图4(b)~(d)所示。

Figure 4. (a) Gravity anomaly results; (b) The traditional Euler method inversion result (N = 1); (c) Tilt-Euler method inversion results; (d) Improved Tilt-Euler method inversion results

图4. (a) 重力异常结果；(b) 传统欧拉法反演结果(N = 1)；(c) Tilt-Euler法反演结果；(d) 改进Tilt-Euler法反演结果

以重力异常为基础数据，选择构造指数为1，求解欧拉方程，得到的反演结果如图4(b)所示。从图中可以看出，传统欧拉法在计算水平位置时存在较大误差，且在深度方向上出现了大量的发散解，导致识别的模型深度远大于实际深度。尽管传统欧拉法能够在有限的先验信息下识别和判断场源体的轮廓，但其反演解的发散性和不稳定性问题显而易见，影响了反演结果的准确性。

以Tilt为基础数据，联合建立欧拉方程并进行求解，得到的反演结果如图4(c)所示。图中可以看出，Tilt-Euler法的反演结果较为理想，边界轮廓较为清晰，但在边界处仍然存在一些发散解，导致识别的深度为14.4 km，明显大于实际的深度。以改进Tilt为基础数据，结合欧拉方程进行求解时，反演结果如图4(d)所示。图中显示，改进Tilt-Euler法有效抑制了发散的欧拉解，反演效果较为稳定，识别的深度为12.2 km，略大于实际深度，表明该方法在深度识别上表现较为精确。

3.2. 多个模型试算

在实际的地球物理工作中，地质体都是以比较复杂的形态存在，因此从理论上对复杂形态的地质模型进行试算是必要且有意义的，同时可以用以检验改进Tilt-Euler法的可行性。

设置100 km × 100 km的网格区域，模型1、2、3的顶部埋深分别为2、4、6 km，剩余密度分别为0.3、0.4、0.5 g/cm³，其中在x轴(−50, −30)，y轴(30, 50)范围内，模型1底部与模型2顶部相互接触，具体参数见表3。

Table 3. Three cube combination model data

表3. 三个立方体组合模型数据

模型编号	顶面埋深 (km)	厚度 (km)	x轴范围 (km)	y轴范围 (km)	剩余密度 (g/cm^3)
1	2	2	(−50, 50)	(30, 50)	0.3
2	4	2	(−50, −30)	(−50, 50)	0.4
3	6	3	(20, 50)	(−50, −20)	0.5

组合模型的平面位置如图5(a)所示。从重力异常图(图5(b))可以看出，模型的异常范围在0至41 mGal之间。由于模型1埋深较浅且密度较小，其反映的重力异常强度最弱；而模型3的密度最大，因此其中心的重力异常较强，约为34 mGal。在模型1与模型2接触的区域，由于场源异常的叠加效应，产生了较大的重力异常，甚至超过了模型3中心处的最大值，达到约41 mGal。

Figure 5. (a) Combination model plane position diagram; (b) Composite model gravity anomaly map

图5. (a) 组合模型平面位置图；(b) 组合模型重力异常图

3.3. 不同欧拉反演结果对比

由于实际地层中存在各种噪声干扰，这些噪声可能导致边界识别效果变得模糊。为模拟实际情况，本文在模型中添加了2%的高斯噪声，以便对比分析上述反演方法在噪声干扰条件下的应用效果。向重力异常数据中添加了2%的高斯噪声，这导致边界识别结果变得模糊。为了减少高斯噪声对边界识别效果的影响，本文采用了向上延拓的方法处理加入噪声的数据，具体是将数据向上延拓3 km。

为了对比分析模型加噪后的传统欧拉反褶积法、Tilt-Euler法和改进Tilt-Euler法的计算结果，本文分别以添加2%高斯噪声并向上延拓3 km后计算得到的重力异常、倾斜角和改进倾斜角为基础数据，联合建立欧拉方程，求解欧拉方程，得到的计算结果如图6~8所示。

3.3.1. 传统欧拉反褶积反演结果

图6(a)展示了重力异常加噪并向上延拓3 km后的结果图。可以看到，经过噪声加入和延拓处理后，识别出的边界发生了严重的扭曲，已不再呈现立方体形状，且无法准确识别模型1与模型2接触区域的边界。

Figure 6. (a) Results of gravity anomaly with noise extension of 3 km; (b) Vzx; (c) Vzy; (d) Vzz; (e) THDR; (f) Add noise and extend 3 km to the traditional Euler inversion result (N = 1)

图6. (a) 重力异常加噪上延3 km结果图；(b) Vzx；(c) Vzy；(d) Vzz；(e) THDR；(f) 加噪上延3 km传统欧拉反演结果(N = 1)

通过联合重力异常数据与Vzx (图6(b))、Vzy (图6(c))以及Vzz (图6(d))，建立了相应的欧拉方程，并使用滑动窗口大小为11 × 11，误差范围设定为15%。最终得出的构造指数N = 1时的反演结果如图6(f)所示。从图中可以看出，加入噪声后的传统欧拉法反演效果较差，模型边界出现了大量的发散解，导致无法准确识别模型1和模型2接触区域的边界。

3.3.2. Tilt-Euler法反演结果

图7(a)展示了重力异常加噪并向上延拓3 km后的Tilt计算结果。从图中可以看出，经过噪声加入和延拓后，Tilt方法识别的模型边界发生了严重扭曲，最大异常值位于模型中心，但无法准确识别模型1与模型2接触区域的边界。同时，模型3的识别效果较差，识别出的边界范围大于实际范围。以Tilt为基础数据，构造指数设为0，滑动窗口大小为11 × 11，误差范围限制为15%，反演结果如图7(f)所示。从图中可以看出，经过加噪后的Tilt-Euler法能够较准确地识别模型的平面位置，相比传统欧拉法，其深度方向的反演解更加收敛，且识别出的深度更接近实际深度。然而，在模型边界的外侧仍然出现了大量虚假解，尽管如此，Tilt-Euler法仍然能够识别出模型1和模型2的接触区域边界。

Figure 7. (a) Tilt result after adding noise and extending by 3 km; (b) Tilt-kx; (c) Tilt-ky; (d) Tilt-kz; (e) Tilt-THDR; (f) Results of Tilt-Euler inversion with added noise and extension

图7. (a) 加噪上延3 km后Tilt结果图；(b) Tilt-kx；(c) Tilt-ky；(d) Tilt-kz；(e) Tilt-THDR；(f) 加噪上延Tilt-Euler法反演结果图

3.3.3. 改进Tilt-Euler法反演结果

图8(a)展示了重力异常加噪并向上延拓3 km后的改进Tilt计算结果。从图中可以看出，经过噪声加入和延拓后，改进Tilt能够识别模型1和模型2接触区域的边界，但出现了边界位置偏大的问题。虽然识别了接触区域，但识别结果存在一定的偏差。从图8(f)中的改进Tilt-Euler法反演效果图可以看出，尽管该方法能够识别模型的边界，边界外侧仍然存在大量发散解，导致识别效果不够精确。虽然能够勉强识别出模型1和模型2的接触区域，但对于埋深较深的模型3，识别效果较差，反演得到的深度大于实际深度，表明该方法在深埋地质体的识别上仍存在一定局限性。

Figure 8. (a) iTilt results after adding noise and extending by 3 km; (b) iTilt-kx; (c) iTilt-ky; (d) iTilt-kz; (e) iTilt-THDR; (f) Image of iTilt-Euler inversion results with added noise and extension

图8. (a) 加噪上延3 km后iTilt结果图；(b) iTilt-kx；(c) iTilt-ky；(d) iTilt-kz；(e) iTilt-THDR；(f) 加噪上延iTilt-Euler法反演结果图

4. 结论与认识

通过以上单个立方体模型、不加噪声和加噪声的组合模型试算，总结如下：

在单个立方体模型试算中，改进Tilt-Euler法对于浅部埋深的识别结果准确，随着场源埋深的增加，识别的结果偏向于模型边界外侧；在加噪的组合模型试算中，改进Tilt-Euler法使发散的欧拉解得到了较好的压制，让场源解更集中地分布在模型边界处，对于存在多个场源体的情况能有较好的识别。

在本研究中，提出的改进Tilt-Euler法在重力异常数据分析中表现出较好的适用性，尤其在处理高噪声数据和复杂地质结构时具有显著优势。然而，该方法也存在一定的局限性。例如，对于极端地质环境中深层异常的处理精度仍有待提高，且在面对大规模数据时可能会出现计算效率瓶颈。此外，当数据质量较差或异常特征极为复杂时，方法的稳定性可能受到影响。

未来的研究可以集中在以下几个方向：首先，结合先进的机器学习技术，进一步提升方法在复杂地质条件下的精度与稳定性。其次，探索该方法在其他类型地球物理数据(如磁力异常、地震波速度等)中的应用，推动跨领域的技术融合。最后，针对大规模数据集的处理，优化算法的计算效率，以满足更高时效性和更大数据量的需求。通过这些改进，能够进一步拓宽该方法的应用范围，并提高其在实际地质勘探中的实用价值。