摘要
重力勘探在矿产勘探中被广泛用于探测地下密度的变化，这些变化可能指示潜在的矿化带或成矿构造。本文对内蒙古下都伦胡东地区布格重力异常数据进行了系统处理与分析，采用位场数据处理技术，包括向上解析延拓、剩余重力异常提取和垂向二阶导数计算等方法。重点研究了欧拉反褶积法的理论基础，详细推导了其数学模型，并探讨了该方法的应用优势及局限性。通过单体模型试验，对比分析了不同构造指数对反演结果的影响；组合模型试验验证在有干扰场源的情况下，欧拉反褶积的适用性。最终，基于欧拉反褶积法对研究区断裂构造进行了精细解释，识别出一条区域性深大断裂和数条次级断裂，并据此圈定了6个具有找矿潜力的成矿远景区。综上所述，本研究不仅系统验证了欧拉反褶积法在复杂地质条件下的适用性与有效性，也为内蒙古下都伦胡东地区断裂构造识别与成矿远景预测提供了可靠的地球物理依据，对区域矿产勘查部署具有实际指导意义。
Abstract
Gravity exploration is widely utilized in mineral prospecting to detect subsurface density variations, which may indicate potential mineralization zones or ore-forming structures. This study systematically processed and analyzed the Bouguer gravity anomaly data from the Xiadulunhudong area of Inner Mongolia, employing poten‐tial field data processing techniques including upward continuation analysis, residual gravity anomaly extraction, and vertical second derivative calculation. The research focused on the theoretical foundations of Euler deconvolution, deriving its mathematical model in detail while discussing the method's advantages and limitations. Through single-source model tests, the influence of different structural indices on inversion results was comparatively analyzed. Combined-source model experiments further validated the applicability of Euler deconvolution under scenarios involving interfering field sources. Ultimately, based on Euler deconvolution, the fracture structures in the study area were meticulously interpreted, identifying a regional deep-seated fault and several secondary faults. These findings led to the delineation of six prospective metallogenic zones with mineral exploration potential Therefore, this work demonstrates the utility of Euler deconvolution in complex geology and provides a robust geophysical basis for delineating fault structures and identifying promising mineralization targets in the study area, with direct implications for guiding future mineral exploration in Inner Mongolia.
Keywords
0 引言
内蒙古是中国重要的矿产资源大省之一，其中金属矿产资源储量丰富（许立权等，2016；梅燕雄等，2024）。在金属矿勘探过程中，寻找矿体的位置和分布是一项重要的任务。重力勘探是以地质体与围岩的密度不同为基础来推断地下构造和矿产分布情况的一种地球物理勘探方法（曾华霖， 2005）。同时，作为一种无损性检测的勘探方法，具有较高水平高分辨率的特点，是矿产资源勘探中不可或缺的技术手段之一（柳建新等，2016）。研究区地质环境独特，重力场异常分布较为复杂（祁民等， 2009），因此对重力勘探的处理与解释提出了更高的要求。
重力数据处理方法的研究一直是地球物理勘探学者研究的热点（邹博伦，2019；代宏睿 2021；陈茂山等，2023）。重力具有勘探深度大、数据可信度高、不受电磁干扰的特点，重力异常能反演地下岩体的分布及区域断裂特征（马健等，2024）。重磁数据处理的核心目标在于从复杂叠加场中提取目标地质体的有效信息，其流程通常包含“位场分离”与 “边界识别”两个关键阶段。位场分离通过匹配滤波、小波多尺度分解等方法（Fedi and Quarta，1997），将区域场与局部场进行剥离，能够消除深部背景场或浅表干扰源的影响，从而增强目标异常的显著性；而边界识别则通过导数类算子对异常梯度带进行空间增强（如 THDR、VDR、ASM 等），其本质是通过数学变换凸显地质体的物性突变界面（Hinze et al.，2013）。这两类预处理技术为后续定量反演提供了高信噪比的输入数据，尤其对欧拉反褶积法 （Thompson，1982）这类依赖于场源几何特征的反演算法具有特殊意义——位场分离可降低深部场源对欧拉解算的干扰，而边界识别结果则能有效约束反演窗口的滑动范围（Reid et al.，2014）。例如，在澳大利亚 Yilgarn 克拉通的案例中，Cooper（2016）通过 NSTD 边界识别限定了欧拉反演的搜索半径，使场源深度解算误差从±18% 降至±7%。曾华霖和许德树（2002）提出了估计最佳向上延拓的方法；徐世浙（2006）提出的“积分-迭代法”在重力数据延拓中取得了非常好的效果。该方法通过对重力数据进行积分与迭代计算，不仅能够有效地提高数据处理精度，而且在复杂地质背景下具有较强的稳定性； 孙争等（2023）提出了一种确定最佳延拓高度的方法，以使向上延拓分离深部重力异常效果最佳。在边界识别方面，Nabighian（1972）提出的解析信号振幅法（Analytic Signal Amplitude，ASM）通过将磁场水平导数与垂直导数进行平方和开方运算，有效增强了地质体边界的垂向分辨率，该方法在磁异常边界识别中展现出独特优势；Hood and Teskey（1989）用总水平导数垂向导数的零值线，勾勒出了地下地质体的边界；Miller and Singh（1994）进一步提出倾斜角（Tilt angle）概念，定义为垂直导数与总水平导数的反正切，其对微弱异常信号具有显著放大作用，尤其适用于深部构造识别。Verduzco et al（.2004）发展的 Theta map 技术通过组合不同阶次导数实现了地质体边界的三维空间定位，在墨西哥湾盐丘构造识别中验证了其复杂构造解析能力；孙鹏飞和吴燕冈（2007）利用重磁的水平和垂直二阶导数划分了东北地区梯度带；Cooper and Cowan（2006）提出了归一化标准差法，并利用澳大利亚的航磁数据和南非的重力数据对NSTD进行了演示，结果表明NSTD能够更细致的描绘出地质体的边界；王万银等（2010） 利用二阶导数和梯度分析技术识别地质体边界方面取得了显著进展，提出的基于梯度带划分的算法，能够在复杂地质条件下精确地分辨不同类型的异常源边界，为矿产勘探提供了重要的技术支持； 针对传统方法的噪声敏感性问题，马国庆等（2013） 提出增强型均衡滤波器，通过引入自适应滤波算法将信噪比提升了40%以上，在中国四川盆地的应用中成功识别出区域地质构造的水平位置。
欧拉反褶积法自提出以来，在国内外得到了广泛的应用与改进。Thomphon（1982）最初提出的欧拉反褶积法，通过利用局部重力异常数据反演地下异常源的位置和形态，成为勘探中不可或缺的重要工具。在此基础上，国内外学者不断优化其算法，以提高反演精度。例如，李银飞（2019）通过通过引进分类算法中的 K 邻近算法，提高欧拉反演结果的准确性；马国庆等（2020）提出了一种基于欧拉反褶积法的优化算法，重力梯度张量联合欧拉反褶积能够有效抑制噪声对结果的干扰；金文博（2024）基于相关成像的等效源法能够提高起伏地形下的重磁异常导数计算精度，提高了欧拉反褶积方法的精度。
不同的构造指数对应着不同的场源，为了计算方便，Barbosa et al（.1999）提出了自动估算构造指数的方法来计算总场值与背景场之间的关系去获得最佳解；也有学者指出一些简单场源的构造指数范围是1~3（Reid et al.，2014）；姚长利等（2004）采用变构造指数对数据滑动窗口重复扫描，从而使得场源解更集中；周文纳（2014）通过联立不同高度的欧拉方程组，建立超定方程组来消除构造指数的影响，从而提升了反演结果的准确性。人工选定构造指数不足之处是在不能准确地判断地质情况时，特别是复杂地质构造情况下，容易引入人为干扰等不确定因素，无法快速准确地选取构造指数，需要重复试验，降低了反演效率，并且在不同的研究中，构造指数的取值也是不同的。综合上述研究，欧拉反褶积方法经过多次关键性改进，不仅在应用方便性上得到了提升，而且在反演结果的稳定性和准确性方面也有显著提高。
本研究的目的是在常规方法的基础上，使用欧拉反褶积快速找到成矿远景区，将有助于更好地了解该研究区的地质结构和矿产潜力，并为进一步的勘探活动提供指导。
1 区域地质概况
1.1 地质背景
研究区位于内蒙古自治区中部，构造位置处于华北板块北缘中段，位于西拉木伦河断裂以南的包尔汉图—白乃庙带岛弧带内（图1）。区内地层、岩浆活动及构造演化过程复杂，地质演化历史跨越古生代—新生代。岩性组合以沉积岩、火山岩和深成岩为主，其中燕山期岩浆-火山活动最为强烈，广泛发育花岗岩、闪长岩及石英岩等侵入岩体，为区内主要赋矿岩石类型（图2）。表层覆盖物以红色黏土及灰绿色、黄色砂土为主。区域矿产资源类型丰富，主要包括铁、铜、锌、金、硼、镁、稀土和煤等（朱卫平，2007）。成矿作用具有明显的多期性特征，元古宙为区内重要的成矿期，晚古生代—中生代亦形成了一些热液型矿床（沈存利等，2004）。不同成因、不同期次岩石单元在空间上相互叠置，砂砾岩与花岗岩广泛接触，反映了研究区在长期构造-岩浆作用影响下形成了复杂的地质结构格局及多类型矿床的叠加分布特征。
1.2 布格重力异常
基于收集到的下都伦胡东地区1∶10000布格重力异常图（图3）显示，研究区重力场特征具有明显的分区性，整体呈现西北—东南向的分带规律。根据重力异常形态特征及空间分布规律，可将研究区划分为3个显著的异常带：
（1）西北部重力高异常区：位于研究区的西北部，表现为显著的正异常特征，异常幅值较大，圈闭等值线形态复杂，出现多处同型扭曲现象，地质图中显示该处由燕山期中细粒—中粗粒闪长岩以及部分花岗细晶岩构成。表明该区存在高密度地质体的局部富集，且受后期构造改造作用明显，具有复杂的构造-岩浆活动特征。
（2）中部重力梯级带：该带呈北东向展布，贯穿整个研究区，表现为密集排列的重力等值线，反映明显的密度界面特征。这种线性梯度带通常指示深大断裂构造的存在，可能为重要的控岩控矿构造，对成矿物质的运移和富集具有重要控制作用。
1 —华北地台边缘带；2—阴山隆起带；3—白云鄂博—赤峰带；4—包尔汉图—白乃庙带岛弧带；5—温都尔庙—林西带；6—锡林浩特构造带； 7—混杂堆积；8—板块俯冲带；9—华北板块与西伯利亚板块缝合线；10—超基性岩体；11—蓝闪片岩；12—研究区
图2研究区地质及矿权分布图
1 —第四系更新统砂砾岩；2—燕山期中细粒—中粗粒黑云母花岗岩；3—燕山期细粒—中细粒石英闪长岩、石英二长岩；4—花岗细晶岩脉；5— 已有矿权位置
图3研究区布格重力异常图
（3）东南部重力低异常区：该区域呈现明显的负异常特征，可能反映低密度地质体。该异常带与西北部重力高异常区形成鲜明对比，地质图显示此处出露石英闪长岩、石英二长岩，暗示研究区深部构造可能存在显著的不均一性。
上述重力异常特征为研究区构造格架划分和成矿远景预测提供了重要的地球物理依据，特别是重力梯级带的识别对确定控矿构造具有重要意义。
2 重力数据的处理解释
2.1 位场分离
2.1.1 解析延拓
向上延拓相当于一个低通滤波器，起到压制以高频成分为主的浅层局部异常，相对突出深层区域场的作用（蒋甫玉，2008）。为了了解深部地质情况，本文对数据分别进行向上延拓 100 m、200 m、 300 m、500 m、800 m的处理，由于篇幅有限，本文仅展示100 m和800 m的结果（图4）。
通过分析延拓结果可知，该研究区基底较浅，布格重力异常值呈西北高、东南低之势，并且上下层继承良好；地表 0~100 m 之间重力异常的形态比较复杂，但是当延拓到800 m时，此时异常值变得较为简单，呈斜坡状，说明从此处开始地质体引起的异常可以当做区域场来进行处理。
2.1.2 剩余重力
剩余重力异常不仅能揭示浅部地质体的影响，还能反映不同深度地质体的作用，分别用向上延拓 0 m、100 m、200 m、300 m 的重力异常对应减去向上延拓 100 m、200 m、300 m、400 m 的重力异常，得到向上延拓 0～100 m、100～200 m、300～400 m 各个层位的重力异常，它们大致反映了地下 0～50 m、 50~100 m、100~150 m、150~200 m 深度各个层位密度不均匀体引起的重力异常变化情况（图5）。
对比不同深度的剩余重力异常结果（图5）可知，浅部异常主要集中分布于图5a～c所示深度范围内，异常呈现出西北、东南向重力值偏高，且以正异常为主。该类正异常表明研究区局部存在相对高密度体，其重力值明显高于周围背景场，推测与浅部高密度矿化地质体有关。综合异常分布特征及深度分解结果判断，该区潜在矿体主要赋存于 0～300 m深度范围内。
2.2 边界识别
通过解析延拓可知，该研究区的异常主要在浅部，而垂向二阶导的主要目的是突出浅部小的地质体所产生的异常，压制区域异常，并且其零值线大概反映出异常体的边界，可对布格重力异常向上延拓到不同高度，并求取垂向二阶导，图6为上延 100 m、200 m、300 m的垂向二阶导数图。
图4向上延拓成果图
a—向上延拓100 m的布格重力异常；b—向上延拓800 m的布格重力异常
图5分层剩余重力异常图
a—向上延拓0~100 m层位剩余布格重力异常；b—向上延拓100~200 m层位剩余布格重力异常；c—向上延拓200~300 m层位剩余布格重力异常；d—向上延拓300~400 m层位剩余布格重力异常
图6垂向二阶导数图
a—布格重力异常垂向二阶导；b—向上延拓100 m后的垂向二阶导数图；c—向上延拓200 m后的垂向二阶导数图；d—向上延拓300 m后的垂向二阶导数图
通过对布格重力异常数据及不同延拓高度下的垂向二阶导数法处理，成功提取出了主要地质体的边界特征。这些边界信息揭示了地下不同深度处的重力异常分布，有助于对深部地质体的结构进行初步推断。然而，从图6a、b中也可以观察到几个较小的重力低值圈闭，它们在空间分布上呈现出局部异常特征。这些异常可能代表了较浅源的地质体或次生构造体的影响，提示了可能存在的浅层地质构造或物质分布的变化。通过对比地质图，发现花岗岩晶脉群与布格重力异常的垂向二阶导对应良好，表明晶脉沿断层或裂隙分布，导数图也揭示了这些构造的特征。
3 欧拉反褶积
欧拉反褶积方法是一种基于位场异常数据的反演方法，广泛应用于地球物理探测，特别是在重力与磁力异常分析中。其基本思想是通过分析位场异常数据的梯度信息，推导出场源体的位置及其构造指数。Thompson（1982）提出了该方法的基本理论框架，基于欧拉齐次方程（Euler Homogeneous Relation）来解决这一问题。
考虑一个三维空间中的函数 f（ x，y，z），如果它满足以下关系：
(1)
式（1）中，t 是一个常数，n 是一个整数，那么该函数 f （ x，y，z）就是一个齐次函数，阶数为 n，可以得出该函数满足以下欧拉方程：
(2)
则该方程被称为欧拉齐次方程。
假设在一个观测平面 z=0 上存在一个点源，其源点位置为（ x0，y0，z0），该点源在该平面上的磁场异常 f（ x，y）可由以下公式表示：
(3)
式（3）中，C 是常数，N 是构造指数，且 C 与坐标 x，y，z 无关。通过对该异常函数 f（ x，y）进行偏导数计算，并代入欧拉方程（2），可以得到该点源的欧拉方程：
(4)
式（4）表明，通过建立线性方程，可以推导出点源的位置（ x0，y0，z0）。其中 N 称为构造指数，不同的地质构造对应不同的构造指数，表1为构造指数对应的地质构造。
表1常见的场源构造指数
3.1 理论模型试验
3.1.1 单体模型
为了了解不同的构造指数对欧拉解的影响，本文设计了球体理论模型，选取不同的构造指数对其进行欧拉反演。建立三维直角坐标系，向东1 km为 x 轴，向南 1 km 为 y 轴，垂直地面向下为 z 轴。下面的模型数据都是基于此坐标系，模型一参数如表2所示，位置如图7所示。
表2模型一参数
图7模型一位置图
以模型重力场及其 x、y、z 方向导数为基础数据，如图8所示，建立欧拉方程，构造指数分别选择 0、0.5、1和2对模型进行欧拉反演，结果如图9所示。
反演结果可以看出，当构造指数取值在0~1时，解的发散现象较为明显，边界识别效果较差。而当构造指数为 2 时，反演对球状地质体的解析效果最佳，能够清晰地刻画其边界，并较为准确地估算埋深。这一结果表明，欧拉反褶积反演过程中，对构
图8模型一正演图
a—理论重力异常；b—垂直导数；c—x方向一阶导数；d—y方向一阶导数
图9模型一欧拉解
a—构造指数0；b—构造指数0.5；c—构造指数1；d—构造指数2
3.1.2 组合模型
然而，在实际情况中，更多的是多个场源叠加在一起，一个计算窗口内可能包含多个异常。上文中使用单体模型探讨了构造指数对欧拉解的影响，下面将使用组合模型进行探究。长方体也是一种非常重要的地质模型，其不同的形态可以模拟不同的地质体。当长宽高相同时，可以模拟等轴状地质体；当厚度较薄而水平展布面积较大时，可以模拟岩床；而当长宽较小高度较大时，可以模拟岩柱、岩脉及岩墙等。为了实验对于有其他场源干扰时欧拉反褶积的应用效果，本文又设计了长方体组合模型试验，模型由 4 个长方体组成，其中一个正方体，一个垂直柱体和一个水平柱体模拟局部异常，一个较深的薄板横跨整个模型研究区模拟区域异常，位置如图10所示，具体参数见表3。
以模型重力场及其 x、y、z方向导数（图11）为基础数据，建立欧拉方程，得到反演结果如图12所示。
图10长方体组合位置图
表3模型二参数
图11模型二正演图
a—理论重力异常；b—垂直导数；c—x方向一阶导数；d—y方向一阶导数
图12模型二欧拉解
从图12中可以看出，欧拉反褶积能够将薄板模型所产生的区域异常很好地分离出去，对于地质体边界识别的效果也更好，也验证了欧拉反褶积的有效性。
3.2 欧拉解
为了深入探究研究区域内的断裂构造及边界特征，基于理论模型实验成果与欧拉反褶积快速反演原理，选定构造指数1作为参数，对研究区的重力异常数据进行处理，反演出由密度差异所指示的地质边界信息，并将其结果投影叠加至地质简图上，如图13所示。
图13研究区欧拉解
欧拉反褶积结果显示，研究区内存在多组不同深度的场源响应，反映出密度不均匀体在空间上的复杂分布特征。反演结果表明，区内断裂及构造边界发育，异常线性特征清晰，主要集中分布于400～500 m 深度范围内，推测与隐伏断裂或构造界面有关（图13）。
图14钻孔位置
表4部分钻孔资料
结合区域地质资料、布格重力异常及欧拉反褶积结果，对研究区主要断裂构造及成矿远景区进行了综合划分。将反演成果与已获取的钻孔资料进行对比分析（图14），除 ZK1-1 外，其余钻孔揭示的地质体位置与重力异常及反演结果具有较好对应关系（表4），验证了研究方法在本区应用的有效性。同时，欧拉反褶积结果与已有矿权分布位置具有一定一致性，进一步表明该方法在研究区找矿预测中的适用性。
4 成矿远景区预测
断裂对矿产形成起着控制和破坏两重作用。由于地质体具有地球物理特性，当断裂产生后，使地质体在三度空间发生位移和错断时，地层间的物理性质产生变化；断裂反映的是物性界面及地质体分布的陡变和扭曲错断位置，断裂的规模越大，两种物性界面的陡变和扭曲错断位置表现越清晰，重磁异常形态表现的分界性质也越明显。图15为本文所识别出的断裂，可以看出，研究区有一条大的断裂和若干条小型断裂，大型断裂走向为NE—WS，在重力异常图中，反映为重力异常梯度带；在重力异常剩余场分布图中，反映为正负异常分界，此条大的断裂对矿产资源的形成有着决定性作用。
图15本文识别的断裂
综合布格重力异常、向上延拓、垂向二阶导数及欧拉反褶积结果，并结合区域地质资料与钻孔揭示情况，对研究区成矿有利部位进行了综合分析与圈定。结果表明，研究区可划分出6个成矿远景区，分别编号为 A～F，其空间分布特征及范围如图16所示。
在重力剩余场中，这 6 处异常远景区表现为重力负值背景场上的高值圈闭区，通过重力异常垂向二阶导数分布图，排除了实测数据局部场异常的分离尚不彻底的因素，再结合已有矿权的实际开采情况，说明其密度大于围岩密度，确定为高密度矿体。并且图16中的 A~F 与图5a垂向二阶导中的负值等值线对应良好。其中 B 和 F 埋深较浅，在 200~300m；A、C、D 次之，埋深在 300~400 m；E 的埋深最大。从研究区地质图中可以看出，圈定的靶区地表主要出露燕山期中粗粒—中细、中粗粒黑云母花岗岩。而且出露的花岗岩密度变化范围大致为 2.4～3.1 g/ cm3。这种局部高重力异常圈闭现象只能是岩石内部含有较多的相对围岩的高密度矿物。但是物理勘探具有不确定性，具体还需要后期进行钻探工作来提高其精确度。
图16成矿远景区预测
5 结论
本文基于 1∶10000 布格重力数据，采用解析延拓法、分解剩余重力、垂向二阶导数法以及欧拉反褶积等方法对研究区重力数据进行了处理，取得了以下主要成果：
（1）对研究区重力异常进行了系统处理与分析。解析延拓及剩余重力分解结果表明，研究区基底埋藏较浅，基底埋深主要集中在 200～500 m 范围内，为浅部构造及成矿作用发育提供了有利地质条件；
（2）垂向二阶导数处理有效突出浅部异常特征，清晰刻画了主要地质体及构造边界的空间展布特征，为隐伏断裂及控矿构造的识别提供了重要依据；
（3）欧拉反褶积结果定量约束了主要异常体的空间位置及埋深特征，其解算结果与导数分析及剩余重力异常具有较好的一致性，表明该方法在研究区隐伏构造与赋矿地质体识别中具有较好的应用效果。
（4）综合重力异常特征与构造解释成果，推测研究区发育一条大型断裂及多条次级断裂构造，断裂与重力异常梯度带及异常中心关系密切，在此基础上圈定了 6 处成矿远景区，为下一步找矿靶区优选及钻探部署提供了地球物理依据。