1. 引言
重力勘探是以研究对象与围岩存在密度差异为条件，通过观测和研究重力场变化来达到勘探目的的地球物理方法。上世纪中叶以来我国重力勘探持续发展，已逐渐成为地球物理勘探的一种重要手段，并在经济建设领域发挥重要作用 [1] [2]。重力勘探是一种快速、轻便而且有效的获取地下信息的地球物理勘探方法，且伴随着计算机技术的飞速发展和应用，重力数据处理与解释的水平有了极大的提高 [3]。现在，重力勘探被广泛应用于油气勘探，固体矿产资源勘探，大地及区域地质构造研究，水文及工程地质方面等。
重力数据处理方法的研究一直是地球物理勘探学者研究的热点：Dean在1958年提出了线性滤波理论，将解析延拓、二阶导数、圆滑等重力数据处理方法视为滤波 [4]。Gupta和Ramani在重力异常分离中首次对谱分解、向上延拓和图形圆滑三种方法进行了比较 [5]。曾华霖对垂向二阶导数、趋势分析方法和解析延拓法进行了对比，并在2002年提出了最佳向上延拓高度的求取方法 [6] [7]。魏献齐提出了椭圆向滤波器提升了重力资料处理中重力异常分离的效果 [8]。刘祥重提出了根据不同埋深的异常地质体对应不同的重力异常频率分离重力异常的方法 [9]。张凤旭在研究高精度重力异常处理方法时提出了组合滤波的方法，在一定程度上提高了重力归一化总梯度的分辨率，为寻找小型油气田提供了技术支持 [10]。肖锋研究了重力数据处理方法在钾盐矿中的应用，改进了传统小子域滤波法，并首次将Canny算子引入到重力数据处理中 [11]。陈曦研究了重力数据处理与解释在刚果钾盐矿区的应用，取得了很好的效果 [12]。刘璎等通过对云南省江城地区进行高精度重力数据位场分离、边缘检测等方法的处理与解释，验证了该区域深部找钾的可行性 [13]。郇恒飞等在金羊盆地重力数据处理中比较了单方法与组合多种重力数据处理方法划分断裂构造的效果，证明了组合多种重力数据处理方法划分断裂构造的有效性 [14]。
钾在农作物生长中具有促进酶的活化，促进光合作用，改善能量代谢，增强抗寒、抗寒性等多种作用，是农作物生长的必需营养元素。中国作为人口和农业大国，对钾盐资源需求量巨大，虽然全球钾盐资源极为丰富，但中国钾盐资源短缺且分布不均 [15] [16]。钾盐的有效供给对农业发展，粮食安全和民生质量等方面有重大的影响 [17]。因此，钾盐矿的勘探具有重要的现实意义。
钾盐是低密度岩体，钾盐矿厚度达到一定规模时，能在一定范围内产生相应的重力负异常 [18]。因此，通过重力勘探技术研究钾盐矿具有现实意义，且钾盐矿重力勘探数据的处理与解释有利于钾盐矿富集区域的圈定。
由于地面观测到的布格重力异常是地下物质密度分布不均匀的反映，在重力异常解释过程中，为了提取目标场源产生的异常，必须对重力数据进行处理和转换，压制其中的干扰信息，突出所需信息 [11]。本文通过使用不同方法对重力数据进行处理和分析，取得了一定的成果。
2. 重力勘探方法原理
2.1. 解析延拓
根据观测平面的观测重力异常计算高于或者低于它的平面上的重力异常称为解析延拓 [3] [19]。向上延拓即计算观测平面之上的平面，向下延拓即计算观测平面之下的平面。
解析延拓计算公式如下
(1)
式中：
为延拓高度，
为延拓后异常体波数谱，
是实测平面异常体波数谱。
向上延拓相当于一个低通滤波器，可以将由深部场源引起的区域重力异常从观测重力异常中分离出来，即压制浅部小地质体引起的小范围高频的局部异常，突出深部异常特征。向下延拓相当于一个高通滤波器，它起到压制以低频成分为主的深层区域异常，相对突出浅层局部场的作用 [11]。
本文使用Geosoft公司的Oasis montaj软件对原始数据进行二维滤波处理，完成解析延拓，得到了不同深度的重力异常数据，详见本文3.2章节。
2.2. 小波分析
小波分析 [20] 进行位场分离的基本原理是通过离散小波变换将重力异常数据分解为小波逼近和各阶小波细节 [21]，将各阶小波细节之和作为局部场，用总的布格重力异常减去各阶小波细节之和就可以得到区域场。
小波变换 [22] 分为连续小波变换和离散小波变换，在处理实际问题时，为了避免连续积分等复杂过程，通常需要将连续小波转换为离散小波 [23]。
本文使用美国MathWorks公司出品的Matlab软件对重力异常数据进行小波分析，选用coif2小波包分解出4阶和5阶小波细节，然后进行位场分离，详见本文3.3.1章节。
2.3. 归一化标准差
归一化标准差法 [24] 是一种通过对设定的窗口内数据的不同方向导数进行统计，最后对垂向一阶导数的标准偏差归一化的边界识别方法，即NSTD (normalized standerd deviations)。其计算公式如下：
(2)
式中
为重力场；
为标准差。
用归一化标准差法处理重力异常数据时，设定的窗口依次移动，数据光滑时，NSTD的值相对较小，在异常数据处，NSTD值相对较大 [25]。因此，NSTD极大值代表了构造边界位置。
本文使用美国MathWorks公司出品的Matlab软件对重力异常数据进行处理，选用3 × 3的窗口依次计算窗口内三个方向位场梯度的标准差，归一化为窗口中心的NSTD值，然后依次移动窗口得到各个窗口中心的NSTD值，详见本文4.2章节。
3. 重力异常处理与解释
3.1. 原始资料的预处理
数据网格化要求我们所测得的重力异常数据呈规则分布，但是在实际测量中，由于一些区域无法进行测量工作(如湖泊，断崖等)，会导致测点分布不均匀。由实测数据插值计算野外勘探工作无法测得的测点数据，使测点数据呈规律性分布，此过程即为数据网格化。经sufer软件网格化后的研究区布格重力异常如图1所示。
Figure 1. Bouguer gravity anomaly in the study area
图1. 研究区布格重力异常图
3.2. 解析延拓及初步解释
本文采取对研究区重力异常向上延拓100 m、200 m、300 m、500 m、800 m和1000 m六个不同高度来分析研究区的基底起伏情况，向上延拓不同高度的重力异常如图2所示。
向上延拓的结果是得到向上延拓高度一半深度的重力异常，图2各图大致对应地下50 m、100 m、150 m、250 m、400 m和500 m深度的密度不均匀特征。对比图2中各图可知，研究区基底较浅，由西向东基底埋深逐步下降，地下0 m到100 m地质体情况较为复杂，负异常出现在研究区的东南部，是重点找矿远景区。
图2(a)中由北西向东南存在三个阶梯状负异常圈闭(中心分别位于y = 247500，x = 203900；y = 249000、x = 203300与y = 251000、x = 2028000)，可能是最有利的找矿靶区。对比图(a)-(f)可知：中间负异常圈闭区域(中心为y = 249000、x = 203300)的强度与面积是三个负异常圈闭区域中最大的，其异常中心非常稳定，根据向上延拓时负异常圈闭的变化速率，判断该处低密度矿体产状较陡，向上延拓800 m以后逐渐消失，说明该处是寻找钾盐矿的有利位置，层底深在400 m以内；研究区最南侧负异常圈闭区域(中心为y = 251000、x = 2028000)也存在相当规模负异常圈闭，上延500 m以后逐渐消失，说明该处地下可能存在低密度矿体，层底深在250 m以内，但该区域处于研究区边缘部位，还需进一步勘探；研究区最北部的负异常圈闭区域(中心为y = 247500，x = 203900)强度较弱，范围也较小，上延300 m以后逐渐消失，说明该处地下可能存在低密度矿体，但埋藏较浅，层底深在150 m以内。
3.3. 位场分离和沉积盆地推测
3.3.1. 小波分解异常
本文采用四阶小波和五阶小波分解研究区的区域重力异常，研究区区域异常由布格重力异常减去剩余异常所得，见图3(a)和图4(a)，将各阶小波细节之和相加作为研究区的剩余异常，见图3(b)和图4(b)，通过分析研究区基底情况，以及结合小波分解得出的剩余重力异常给出5个靶区，见图3(b)和图4(b)中的红色圆圈符号及其标注A、B、C、D、E。对比发现四阶小波分解所得到的剩余异常相比五阶小波分解所得到的剩余异常缺失了众多的局部细节，两者中五阶小波分解所得到的剩余异常相较于四阶小波分解所得到的剩余异常更适合作为研究区的局部异常。
(a) 
(b)
Figure 3. Fourth order wavelet decomposition anomaly. (a) Regional anomaly; (b) Residual anomaly
图3. 四阶小波分解异常图。(a) 区域异常；(b) 剩余异常
(a) 
(b)
Figure 4. Fifth order wavelet decomposition anomaly. (a) Regional anomaly; (b) Residual anomaly
图4. 五阶小波分解异常图四阶小波分解异常图。(a) 区域异常；(b) 剩余异常
3.3.2. 向上延拓分解异常
从布格重力异常中去掉区域重力异常，就可得到剩余重力异常，研究区采用向上延拓500 m和800 m分离的区域重力异常见图5(a)和图6(a)，分离的剩余重力异常见图5(b)和图6(b)。
对比图5和图6可知，从区域异常来看，向上延拓800 m的区域重力异常保留的局部异常更少，更能反映出研究区的基底起伏情况；从剩余异常来看，向上延拓800 m的局部异常比向上延拓500 m局部异常突出，对于寻找低密度体更为清晰容易。
(a) 
(b)
Figure 5. Upward continuation of 500 m to decompose gravity anomaly. (a) Regional anomaly; (b) Residual anomaly
图5. 向上延拓500 m分解重力异常图。(a) 区域异常；(b) 剩余异常
3.3.3. 利用剩余重力异常推测钾盐沉积盆地的位置和范围
通过对比小波分析分解异常和向上延拓分解异常两种方法所得的剩余异常，可以发现向上延拓800 m分离的剩余异常相较于五阶小波分解所得到的剩余异常保留更为完整，所以向上延拓800 m分离的剩余异常是作为确定沉积盆地位置的基础资料的最佳选择，其所得的剩余异常的差值较小。将剩余异常图中小于零的布格重力异常提取出来，即可得到用剩余重力异常预测的低密度钾盐靶区范围，见图7。
从图7中可见，研究区中的小于零的重力异常圈闭主要分布在研究区的中东部，其他区域仅有小面积存在。这些大面积的低重力异常区域以负等值线圈闭，可作为寻找钾盐矿的有利部位。面积最大的低值区域位于研究区的中部偏南，呈不规则多中心形态；其余低值异常圈闭比较分散且圈闭面积相对较小。
(a) 
(b)
Figure 6. Upward continuation of 800 m to decompose gravity anomaly. (a) Regional anomaly; (b) Residual anomaly
图6. 向上延拓800 m分解重力异常图。(a) 区域异常；(b) 剩余异常
Figure 7. Low density Sylvite Mine prospect predicted by residual gravity anomaly
图7. 用剩余重力异常预测的低密度钾盐矿远景区
4. 断裂或构造边界识别
4.1. 断裂或构造边界的识别标志
在进行断裂或构造边界的识别时，断裂或构造边界通常具有以下几种标志 [26]：
1) 断裂导致重力异常梯度带的轴线有明显移动；
2) 重力异常等值线同向扭曲；
3) 重力异常在断裂两侧存在明显的差异；
4) 断裂的存在导致高值重力异常突然降低或低值重力异常突然升高。
4.2. 边界识别与钾盐矿成矿远景区预测
了解地下地质体的形态和规模是重力勘探的一项重要任务，而地下地质体的边界信息很难从布格重力异常中直接提取出来，所以选择一种方便且高效的边界识别方法是十分必要的。本文选用了归一化标准差法对研究区断裂或构造的构造边界进行识别。对研究区使用归一化标准差分进行边界识别如图8所示。
从图8中可以看出研究区断裂以北西向为主，而北东向断裂很少。识别结果与五阶小波细节较为吻合。综合对比研究区边界识别结果与向上延拓800 m剩余重力异常和五阶小波分解剩余重力异常，为研究区划分了如下6个钾盐矿成矿远景区，并以红色矩形圈闭(见图9中B1~B6)。将钾盐矿远景区预测(图7)结果与解析延拓结果(图6)对比可知，图9中B1~B6均是由负等值线圈闭的寻找钾盐矿的有利区域。其中B1、B2、B3埋深较浅，B6次之，B4与B5埋深最大。B1位于研究区北部，成矿远景区中心坐标为y = 246500，x = 2038700，主要由北西向断裂控制，底部埋深在250 m以内，东西向长度约2000 m，南北向长度约1000 m；B2位于研究区西北部，成矿远景区中心坐标为y = 245000，x = 2037800，主要由北西向和东西向断裂控制，底部埋深在150 m以内，成矿远景区形状不规则；B3位于研究区东北部，相对于其他远景区，面积最小，成矿远景区中心坐标为y = 248200，x = 2038000，主要由北西向和北东向断裂控制，底部埋深约150 m，储量较少；B4位于研究区中南部，相对于其他远景区，面积最大，但形状不规则呈串珠状分布，成矿远景区中心坐标为y = 246000，x = 2034000，主要由北西向和北东向断裂控制，底部埋深约250 m，与B1、B2、B3区域相比储量较多，是寻找钾盐矿的有利区域；B5位于研究区中东部，相对于其他远景区，面积较大，负等值线圈闭区域较为集中，成矿远景区中心坐标为y = 249000，x = 2034000，主要由北西向和北东向断裂控制，底部埋深约500 m，储量较多，是寻找钾盐矿的最有利区域；B6位于研究区南部，成矿远景区面积仅次于B4区域，负等值线圈闭区域较为集中，成矿远景区
(a) 
(b)
Figure 8. Study area boundary identification. (a) Normalized standard deviation; (b) Fracture distribution
图8. 研究区边界识别。(a) 归一化标准差；(b) 断裂分布
中心坐标为y = 249000，x = 2031000，主要由北西向和北东向断裂控制，底部埋深约400 m，储量较多，是寻找钾盐矿较为有利的区域。在今后的钾盐矿勘探工作中，还应配合钻孔工作来提升远景区划分的精确度。
Figure 9. Sylvite Mine prospect predicted in the Study area
图9. 研究区边界识别
5. 结论
不同的重力处理数据方法在钾盐矿勘探过程中都有其各自的优势与不足，本文采用了多种处理方法对研究区重力数据进行处理，主要得到了如下结论：
1) 通过解析延拓对研究区基底情况做出了初步解释，识别出研究区基底较浅且主体大致为北西向的斜坡，呈东南低西北高的梯度带状。
2) 通过分析对比小波分析和向上延拓分解重力异常的结果，最终选定向上延拓800 m求得的局部异常作为研究区的剩余异常，并根据研究区剩余重力异常初步推测了钾盐沉积盆地的位置和范围。
3) 通过使用归一化标准差法对研究区重力异常数据进行处理，对研究区断裂或构造的边界进行了识别，综合对比研究区边界识别结果与向上延拓800 m剩余重力异常和五阶小波分解剩余重力异常后圈定了研究区内寻找钾盐矿的六个有利位置。
综上所述，重力勘探是一种有效的钾盐矿勘探地球物理勘探方法，选取合适的重力数据处理方法在钾盐矿的重力勘探中可以取得很好的勘探效果。
基金项目
本论文受生态智慧矿山联合基金项目(项目编号：D2020402013)；河北工程大学教育教学改革研究与实践项目(项目编号：JGSZ2022012、JG2021041、JG2022034)；河北省高等教育教学改革研究与实践项目(项目编号：2021GJJG256)联合资助。