CN110888176A - 一种利用地面高精度重力测量的找矿方法 - Google Patents

一种利用地面高精度重力测量的找矿方法 Download PDF

Info

Publication number
CN110888176A
CN110888176A CN201911027784.6A CN201911027784A CN110888176A CN 110888176 A CN110888176 A CN 110888176A CN 201911027784 A CN201911027784 A CN 201911027784A CN 110888176 A CN110888176 A CN 110888176A
Authority
CN
China
Prior art keywords
gravity
field
anomaly
derivative
local
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911027784.6A
Other languages
English (en)
Other versions
CN110888176B (zh
Inventor
聂逢君
王彦国
邓居智
严兆彬
李满根
夏菲
何剑锋
张成勇
封志兵
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
East China Institute of Technology
Original Assignee
East China Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by East China Institute of Technology filed Critical East China Institute of Technology
Priority to CN201911027784.6A priority Critical patent/CN110888176B/zh
Publication of CN110888176A publication Critical patent/CN110888176A/zh
Application granted granted Critical
Publication of CN110888176B publication Critical patent/CN110888176B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G01V7/02Details
    • G01V7/06Analysis or interpretation of gravimetric records

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

一种利用地面高精度重力测量的找矿方法,主要包括六个方法,分别为位场高阶导数的波数域迭代法、基于迭代滤波法的位场分离方法、识别地质体边界的TASD法、位场数据处理的欧拉反褶积法、位场异常互相关成像技术方法和重力人机交互剖面正演模拟方法。本发明的优点是:利用地面高精度重力测量的找矿方法,对多条重力剖面的数据处理和解释,解释方法涵盖波数域导数迭代法、场分离的迭代滤波法,以及基于剩余异常的互相关成像技术应用于重力数据中,取得了较好的处理结果。

Description

一种利用地面高精度重力测量的找矿方法
技术领域
本发明涉及一种找矿技术,具体为一种利用地面高精度重力测量的找矿方法。
背景技术
中国北方砂岩型铀矿的成矿地质背景和条件与国外存在较大的差异,因此在找矿方法的探索上也具有 不同的思路。
中亚地区是盛产层间渗入氧化带砂岩型铀矿,层间氧化带矿床成矿模型主要来自这里,中国核工业1980 年代末引进的砂岩型铀矿成矿理论和找矿方法主要就是参考前苏联在中亚的成功经验。但是,中亚特殊的 成矿地质条件是中国不具备的,比如,中亚主要含矿目的层沉积环境是海陆过渡环境,砂体在侧向上的稳 定性不言而喻。又如美国中西部的怀俄明盆地,盆地属于挤压山间盆地,但盆地中目的层靠近隆起区元古 代富铀花岗岩,蚀源区铀源非常充足。中国北方的情况完全不同,目的层多属于陆相盆地沉积环境,特点 是物源近,相变大,堆积快,侧向连续性较差,不如中亚的稳定。另外,中国北方产铀盆地均分布在中亚 造山带两侧,或中间,盆地隆起区虽然有较多的花岗岩,但多数铀含量并不高,且风化程度也不如美国中 西部的怀俄明。
中国北方砂岩型铀矿成矿环境的复杂多变性决定了我们的找矿思路必须在中亚层间氧化带理论和美国 卷状砂岩型铀矿理论基础上有创新。前苏联在中亚地区的找矿方法主要是“次造山”理论指导下,利用钻 探技术,加上部分后生蚀变地球化学分带方法。由于海陆过渡相环境下形成的砂体分布规律容易掌握,利 用钻孔的成功率比较高。
1980年代末至1990年代初,中国主要是引进中亚的成矿理论与找矿方法,曾经在中国西北和东北许 多地区部署过大量的钻探工作,但收效甚微,主要原因是未能充分考虑中国成矿背景条件特征。1990年末 至2000年代初,我们加强了成矿地质条件的研究,尤其是目的层沉积体系和砂体分布规律的研究,配合钻 孔技术,在二连盆地、开鲁盆地和巴音戈壁盆地收到了很好的效果。
古河道砂岩型是二连盆地最为重要的类型,针对铀矿化发育在古河道砂体中,研究组开创性的利用地 面高精度重力测量与反演方法,准备定位河道的位置与深度,同时反演出砂体及附近的断层发育情况,为 在河道中找寻铀矿体及解释铀矿化的形成提供依据。
发明内容
本发明的目的在于提供一种利用地面高精度重力测量的找矿方法,该找矿方法准备定位河道的位置与 深度,同时反演出砂体及附近的断层发育情况,为在河道中找寻铀矿体及解释铀矿化的形成提供依据。
本发明采用的技术方案如下:一种利用地面高精度重力测量的找矿方法,主要包括六个方法,分别为 位场高阶导数的波数域迭代法、基于迭代滤波法的位场分离方法、识别地质体边界的TASD法、位场数据处 理的欧拉反褶积法、位场异常互相关成像技术方法和重力人机交互剖面正演模拟方法。
进一步的,所述位场高阶导数的波数域迭代法;其特征在于:
基本原理为:在波数域中,位场数据的任意阶任意方向导数波谱可以表示为:
Figure BDA0002248147660000021
其中S(m,p,q)(u,v)和S(u,v)分别为所求导数波谱和原始异常波谱;
Figure BDA0002248147660000022
Figure BDA0002248147660000023
分别为x方向m阶导数算子,y方向p阶导数算子和z方向q阶导数算子;u、v分别是x、y方向的波数;m、p、q均为大于等于零的实数。
由上式可知,求导阶次l和理论导数算子
Figure BDA0002248147660000024
分别为:
Figure BDA0002248147660000025
Figure BDA0002248147660000026
带入原始公式中,则有:
Figure BDA0002248147660000027
一般来说,求导阶次l越大,导数算子
Figure BDA0002248147660000028
引起的高频振荡越明显,则需对高频成分压制越强; 且
Figure BDA0002248147660000029
对不同方向的频率放大作用不同,为此在这里对S(m,p,q)乘以一个包含求导阶次l和导数算子
Figure BDA00022481476600000210
相关的低通滤波器P,即:
Figure BDA00022481476600000211
α≥1,β>0(一般可选α=1,β=1),
Figure BDA00022481476600000212
Figure BDA00022481476600000213
的绝对值数值形式,为了简便起 见,文中省略波数坐标(u,v)。那么经P滤波后的S(m,p,q)近似值(在此称为一阶近似谱)可表示为:
Figure BDA00022481476600000214
Figure BDA00022481476600000215
带入原式中,得S的一阶近似值:
Figure BDA00022481476600000216
再用S和S(1)的差值对
Figure BDA00022481476600000217
进行校正,得S(m,p,q)的二阶近似值:
Figure BDA00022481476600000218
重复以上步骤,得S(m,p,q)的n阶近似值:
Figure BDA00022481476600000219
且有:
Figure BDA00022481476600000220
则得:
Figure BDA0002248147660000031
Figure BDA0002248147660000032
Figure BDA0002248147660000033
分别是第n、n-1次迭代得到的波谱
Figure BDA0002248147660000034
Figure BDA0002248147660000035
的傅里叶反变换,ε为 给定的较小正实数。当迭代过程中满足
Figure BDA0002248147660000036
时终止迭代,此时可认为所求导数
Figure BDA0002248147660000037
模型检验为:采用两个二度体直立柱体的叠加异常进行数值检验。如图7b所示:模型体1长2km,上 顶埋深2km,下底埋深2.5km,上顶面中心坐标(11km,2km);模型体2长2km,上顶埋深2.5km,下底埋 深3.5km,上顶面中心坐标(14km,2.5km);两个模型体的剩余密度均为1.0g/cm3。组合模型体产生的重 力异常如图7a所示。
图8和图9是利用常规方法(未经向上延拓的及向上延拓一定高度的导数换 算)和迭代法得到的垂向二阶、三阶导数的对比结果。从图中可以看出:未经向 上延拓的导数换算得到的异常值(图8a、9a)与理论值相比,除产生明显的边 界效应外,其余数据还存在着明显的波动性,尤其垂向三阶导数(图9a)已无 法识别出有效异常;虽然向上延拓0.5倍点距的垂向二阶导数(图8b)和向上 延拓1倍点距的垂向三阶导数(图9b)分别相对于图8a和图9a而言,异常值 有了较好的改善,但数据的波动性仍较明显;为了进一步压制高频振荡,将垂向
基于迭代法的导数换算得到的垂向二阶、三阶导数(图8d和图9d)与常规方法计算的导数(图8a、b、 c和图9a、b、c)相比较,除了在边部有限数据产生较明显的边界效应外,其余数据点数值与理论值拟合 较好。
上述模型检验结果表明,未经向上延拓的导数换算计算结果稳定性较差,且导数阶次越高,数据波动 性越强,而这也正是由于理论导数算子对高频成分的放大所引起的;虽然向上延拓一定高度的导数异常稳 定性有所提高,但是存在着明显的不保幅性,这则是向上延拓一定高度的滤波算子对低频成分的消减所致; 而迭代法不仅计算结果稳定性强,且具有较强的保幅能力,这也恰是迭代法的导数算子在低频段更能逼近 理论导数因子而对高频成分更能有效压制的结果。
顺便指出的是,随着异常导数阶次的增加,其异常所反映的地质体边界特征也越明显。这种效果在图 8和图9虚线所指示的地质体边界位置可以得到充分的说明。
进一步的,所述基于迭代滤波法的位场分离方法;其特征在于:
基本原理为:令叠加异常u(x,y)的波谱为U(u,v),假设算子H(u,v,α)是一低通滤波器且满足0≤H≤1,那么将H(u,v,α)作用于U(u,v)可以达到场分离的目的,即有:
Ureg(u,v,α)=U(u,v)·H(u,v,α)
这里u、v分别是x、y方向上的波数,α是算子H的滤波参数。按照滤波参数α的选取可以将场 分离情况分为以下三种:①滤波参数α选取不当使得分离出的区域场ureg(ureg=F-1[Ureg],F-1是Fourier反变换)中含有局部场成分;②滤波参数α选取恰当使得叠加异常u分离较为彻底;③滤波参 数α选取不当使得分离出的局部场ured=u-ureg中含有区域场成分。
在实际应用中,由于滤波参数α不易确定而使得场的分离容易出现第一种或第三种情况,使得位场 分离效果不理想。
针对第一种情况,提出这样的改进方案:固定滤波参数α使得区域场ureg含有较多的局部场成分, 那么将H再次作用于Ureg以进一步剥离出局部异常,即有
Figure BDA0002248147660000041
倘若
Figure BDA0002248147660000042
还含 有局部场成分,则可以继续将H作用于
Figure BDA0002248147660000043
依次类推,那么H作用n次的区域场波谱为:
Figure BDA0002248147660000044
只要选合适的作用次数n,即可以使得区域场中的局部场成分较好的剥离出去。
针对第三种情况,改进方案为:固定参数α使得局部场ured含有较多的区域场成分,那么此时需要 将局部异常波谱中的区域场成分剥离出来“返还”给区域场,即
Figure BDA0002248147660000045
此时局部异常波谱变为:
Figure BDA0002248147660000046
若局部场中仍有区域 场信息,则继续对局部场进行剥离并“返还”给区域场,如此一直进行下去,至到剥离“返还”n次时, 满足局部场中没有或几乎没有区域场成分时终止,此时,区域异常波谱
Figure BDA0002248147660000047
和局部异常谱
Figure BDA0002248147660000048
分别为:
Figure BDA0002248147660000049
Figure BDA00022481476600000410
显然,上述第二种改进方案中,局部场剥离和“返还”过程与迭代滤波的思想模式是一致的,因此可 以将这种改进方案称为分离场的迭代滤波法。该方法的原理是:随着迭代次数的增加逐渐将局部场中的区 域异常剥离出来“返还”给区域场,当迭代达到一定次数时,可认为局部场中不包含(或者包含少量)区 域场成分,此时迭代次数为最佳,而随着迭代次数的继续增加,则会导致局部场成分被剥离过多而导致区 域场中含有局部场信息。
对于第一类型改进方案文中暂不予分析和模拟,但该方案同样有着自身的物理意义,例如,给定 H=exp(-r(u,v)·Δd)(r为圆波数,Δd为网格点距;一般可认为向上延拓1个点距的异常值中包 含着较多的局部场信息),则H作用n次的区域场波谱为
Figure BDA00022481476600000411
也就变成向上延拓进行异 常分离了,当然H还可以是其他形式的滤波算子。
文中着重研究第二种改进方案,这里给出低通滤波H形式为:
Figure BDA0002248147660000051
为了满足第二类情况的条件,这里设定α≥1,h≥50Δd,β≥1,而迭代次数(即剥离“返还” 次数)的选取方案采用的是文献[14]的差值互相关系数之差(反映的是分离出的区域场和局部场随迭代次 数的变化)。本文模型试验中选取参数α=β=1,h=100Δd。
模型试验:为了验证迭代滤波法进行场分离的正确性,进行了二维叠加异常试验,并将计算结果与常 规的向上延拓进行对分分析。选取的区域场源为一球体,其参数为:埋深10km,半径3km,质心平面点 位(7km,7km)处,产生的区域重力异常见图10。选取的局部场源为两个球体叠加,球体1的参数为:埋 深1km,半径0.5km,质心平面点位(5km,5km)处;球体2的参数为:埋深2km,半径0.8km,质心 平面点位(10km,10km)处,由浅部地质体1和2产生的局部重力异常见图11。图12是区域异常和局部 异常的叠加重力异常。
采用向上延拓值与理论区域异常的最小均方误差(图13)对应的延拓高度作为最佳高度(0.7km),对 应该延拓高度的延拓值如图14a所示,可以看出,向上延拓得到的区域异常(26a)受局部场影响较大,在地 质体2上方区域场与理论异常值相差高达1.382mGal(23.48%),这将导致分离出的剩余异常(26b)在幅值 上远小于理论异常值(24)。在地质体1和地质体2上方,向上延拓的剩余异常值分别为3.028mGal和2.195 mGal,相对于这两点的理论值3.566mGal和3.587mGal,误差达15.09%和38.81%。
为了进一步消除局部场对区域异常的影响,将向上延拓高度提升为1.2km,得到的延拓值(即区域异 常)和剩余异常分别见图15a和图15b。可以看出,相对于向上延拓0.7km的延拓异常,延拓高度为1.2km 的区域异常受局部异常影响变小,且分离出的剩余异常也明显更加接近于理论值,在地质体1、2上方,剩 余异常值分别变为3.972mGal和3.126mGal,误差减小为11.39%和12.85%。但是区域异常(图15a)不仅形 态上与理论异常(图10)仍相差较大(这说明区域场仍然受局部场的影响),而且幅值也明显小于理论值 (这又说明区域场的一部分信息包含在了剩余场中)。
图17a、图17b是根据差值互相关系数与迭代次数的关系(图16)选取最佳迭代次数(n=17)时得到 的区域异常和剩余异常。可以看出,迭代滤波分离出的区域异常在形态上与理论异常相近,在数值上也相 差不大,最大误差为0.526mGal(6.62%),小于向上延拓0.7km和1.2km的误差;得到的剩余场在地质体1 和地质体2上方异常值分别为3.875mGal和3.264mGal,误差百分比分别为8.67%和9%,也同样小于上 面两种延拓模式得到的剩余异常在地质体1、2上方的误差。
上述模型试验对比分析表明,向上延拓方法得到区域异常和剩余异常存在相互影响的现象。而迭代滤 波法得到区域场受局部场影响较小,可以较好地将叠加异常分离出来,提高了场分离的精度。
进一步的,所述识别地质体边界的TASD法;其特征在于:
基本原理:2008年,Cooper和Cowan另辟蹊径,利用导数标准差的归一化形式进行了场源边界的识 别,获得较好的应用效果,其归一化形式为
Figure BDA0002248147660000052
但该方法存在着识别地质体边界连续性差和易受噪声干扰的缺点,针对于此,笔者提出了导数标准差 的倾角法(Tilt angle of standard deviations),其形式如下:
Figure BDA0002248147660000061
其中改造因子m≥1。
模型检验:在此以上顶界面埋深为0.1km、0.5km和1km的三个边长均为1km的正方体组合为试验模型, 并设定三个模型体剩余密度均为1.0g/cm3,模型体平面位置分别对应于理论重力异常图(图18)中的三 条虚线,以下采用改进型识别方法和常规方法对理论重力异常进行相应处理,并对其计算结果进行对比分 析。
从图18易看出,归一化标准差法(图18b)较清晰地反映出了浅部地质体的边界,且识别出的边界 位置与实际吻合较好;但受干扰影响,对地质体2、3的边界反映模糊;而导数标准差的tilt angle(图 18c)可以清晰地反映出所有地质体的边界。表明了TASD法不仅具有较强的计算稳定性,同时也有着较好 的边界识别效果。
进一步的,所述位场数据处理的欧拉反褶积法;其特征在于:
方法原理:Thompson给出了Euler齐次方程如下:
Figure BDA0002248147660000062
其中,f(x,y,z)是位场异常;
Figure BDA0002248147660000063
分别是异常f在x、y、z三个方向上的一阶导数 值;N为构造指数(又称形状因子),不同的数值代表不同的地质体形状。当N=0时,上式可用来反演断裂 分布点。
模型验证:以一个上顶、下底分别为1km、2km,剩余密度为1.0g/cm3的直立柱体为例进行试验。在 反演中取构造指数N=0,取1×5的滑动窗口。重力异常(图19a)欧拉反褶积结果见图19b。可以看出, 基本上全部的反演点均落在地质体边界附近,仅有少量反演点位于地表处。这说明,尽管在反演中出现了 少许干扰,但该方法从整体上仍可以获得高精度的反演结果。
进一步的,所述位场异常互相关成像技术方法;其特征在于:
方法原理:在地形起伏的测区,笛卡尔坐标的(x,y)平面在基准面上,z取垂直向下为正,假设测区 地下任意坐标为(xq,yq,zq)、体积为vq的第q个点质量的剩余密度为Δσq,则它在测区上任意点 (x,y,z)处的重力异常Δgq(x,y,z)可表示为
Figure BDA0002248147660000064
其中,万有引力常数G=6.67×10-11m3/(kg·s2)。
借鉴Abdelrahman等的相邻最小二乘剩余异常的归一化互相关公式,定义测区上实测重力异常与第q个 点质量在测区上的重力异常的归一化互相关为
Figure BDA0002248147660000071
其中,(xi,yi,zi)为测区第i个观测点的坐标,Δg(xi,yi,zi)为该观测点的实测重力异常, Δgq(xi,yi,zi)为地下第q个点质量在该观测点的重力异常,N为测区观测点总数。
模型验证:以添加了随机干扰的边长为1km的正方形直立柱体的重力异常(图20a)为试验对象,对 互相关概率成像法进行检验。从成像图20b中可以看出,该方法不仅计算稳定性强,同时反演结果可以较 好地识别地质体的质心位置。这表明采用互相关概率成像方法反演地质体赋存状态是有效的。
进一步的,所述重力人机交互剖面正演模拟方法;该项技术的优点是便于将重力异常的处理、转换方 法得出的结果和其它地质、地球物理方法获取的先验信息输入到模型里,形成初始模型。并且根据计算结 果和实际重力异常的差异,随时方便地修改模型,直观地监督和指导正反演过程。
重力人机交互剖面正演模拟流程见图21。重力人机交互正反演技术,其特征在于:重力人机交互正反 演技术主要是依据横截面为多边形的二度体重力异常计算方法来实现的。通过对初始模型计算出的重力效 应与测线上的剩余重力异常进行对比,不断修正模型,直至达到计算出的重力效应与测线上的剩余重力异 常之差满足预定精度。
本发明的优点是:利用地面高精度重力测量的找矿方法,对多条重力剖面的数据处理和解释,解释方 法涵盖波数域导数迭代法、场分离的迭代滤波法,以及基于剩余异常的互相关成像技术应用于重力数据中, 取得了较好的处理结果。同时,依据重力异常导数分析、TASD和欧拉反褶积反演结果,对重力剖面的断裂 位置和产状进行了反演,获得了NL剖面上5条断裂、MD剖面上6条断裂、E400线9条断裂和SH线4条断 裂。另外,以重力资料为依据,结合研究区的地质资料、钻孔资料和电性资料,建立了多条剖面的地质- 地球物理模型;并进行了次一级构造单元的划分,指出了NL剖面和MD剖面的格日勒敖都凹陷东部次级凹 陷均呈现出明显的河道沉积特征,这为古河道的存在提供了地球物理证据。
附图说明
图1为本发明的NL线高程测量结果图。
图2为本发明的MD线高程测量结果图。
图3为本发明的E400线高程测量结果图。
图4为本发明的SH线高程测量结果图。
图5为本发明的CG-5重力仪静态变化曲线图。
图6为本发明的重力仪各点的动态零点位移曲线图。
图7为本发明的模型重力异常及模型示意图,其中a为模型体重力异常,b为模型体示意图。
图8为本发明的垂向二阶导数模型检验图,其中a为未进行向上延拓;b为向上延拓0.5倍点距;c为 向上延拓1倍点距;d为迭代法。
图9为本发明的垂向三阶导数模型检验图,其中a为未进行向上延拓;b为向上延拓1倍点距;c为向 上延拓2倍点距;d为迭代法。实线是理论值;点线是计算值。
图10为本发明的迭代滤波法中二维叠加异常试验的区域重力异常图。
图11为本发明的迭代滤波法中二维叠加异常试验的理论局部异常图。
图12为本发明的迭代滤波法中二维叠加异常试验的叠加重力异常图。
图13为本发明的向上延拓值和理论区域异常误差与延拓高度的关系曲线图。
图14为本发明的向上延拓0.7km得到的区域异常和剩余异常图;其中(a)区域异常;(b)剩余异常。
图15为本发明的向上延拓1.2km得到的区域异常和剩余异常图;其中(a)区域异常;(b)剩余异常。
图16为本发明的差值互相关系数与迭代次数关系曲线图。
图17为本发明的迭代滤波法迭代17次得到的区域异常和剩余异常图;其中(a)区域异常;(b)剩余 异常。
图18为本发明的TASD法与NSTD法识别场源边界结果图,其中(a)重力异常;(b)NSTD;(c)TASD。
图19为本发明的欧拉反褶积反演结果图。
图20为本发明的含有随机噪声的重力异常的互相关概率成像图,其中(a)含30%噪声的重力异常;(b) 互相关成像结果。
图21为本发明的重力人机交互正演模拟流程图。
图22为本发明的NL线相对布格重力异常图。
图23为本发明的MD线相对布格重力异常图。
图24为本发明的E400线相对布格重力异常图。
图25为本发明的SH线相对布格重力异常图。
图26为本发明的NL线一阶导数处理结果图。
图27为本发明的MD线一阶导数处理结果图。
图28为本发明的E400线一阶导数处理结果图。
图29为本发明的SH线一阶导数处理结果图。
图30为本发明的NL线重力异常TASD结果图。
图31为本发明的MD线重力异常TASD结果图。
图32为本发明的E400线重力异常TASD结果图。
图33为本发明的SH线重力异常TASD结果图。
图34为本发明的NL线欧拉反褶积反演结果图。
图35为本发明的MD线欧拉反褶积反演结果图。
图36为本发明的E400线欧拉反褶积反演结果图。
图37为本发明的SH线欧拉反褶积反演结果图。
图38为本发明的NL线断裂处理结果综合分析图。
图39为本发明的MD线断裂处理结果综合分析图。
图40为本发明的E400线断裂处理结果综合分析图。
图41为本发明的SH线断裂处理结果综合分析图。
图42为本发明的NL线重力异常分离结果图。
图43为本发明的MD线重力异常分离结果图。
图44为本发明的E400线重力异常分离结果图。
图45为本发明的SH线重力异常分离结果图。
图46为本发明的NL线剩余重力异常概率成像结果图。
图47为本发明的MD线剩余重力异常概率成像结果图。
图48为本发明的E400线剩余重力异常概率成像(a)与MT电性反演结果(b)图。
图49为本发明的SH线剩余重力异常概率成像结果图。
图50为本发明的图4-31NL线剩余重力异常与构造单元关系图。
图51为本发明的图4-32NL线反演断裂与构造单元关系图。
图52为本发明的图4-33NL线地质-地球物理综合模型图。
图53为本发明的图4-34MD线剩余重力异常与构造单元关系图。
图54为本发明的图4-35MD线反演断裂与构造单元关系图。
图55为本发明的图4-36MD线地质-地球物理综合模型图。
图56为本发明的图4-37E400线剩余重力异常与构造单元关系图。
图57为本发明的E400线反演断裂与构造单元关系图。
图58为本发明的E400剖面地质-地球物理综合模型图。
图59为本发明的SH线剩余重力场与构造单元关系图。
图60为本发明的SH线反演断裂与构造单元关系图。
图61为本发明的SH剖面地质-地球物理综合模型图。
具体实施方式
本发明从该技术的前期准备一直到项目成果出来的一些具体实施附上。
一、野外重力勘探工作;首先,本项目的GPS测地和重力测量组完成了四条剖线的快速静态GPS定位 测量和重力测量,其中包括两条剖线的数据补测工作。野外工作量-完成了计划设计的685个重力测点的数 据采集,完成率达154.60%。共完成了NL线、MD线、E400线和SH线四条重力剖面。其中剖面总长度为105.4km, 总物理点数1172个。
二、工作进展及资料整理;根据项目的设计要求在苏尼特右旗-二连浩特之间开展重力测量工作。由 于测点点位(经、纬度及高程)准确度直接影响到(相对)布格重力异常的精度,因此这里GPS测地工作 完全按照《全球定位系统(GPS)测量规范》(GB/T 18314—2001)要求进行。重力测量按照《大比例尺重力勘 查规范》(DZ/T 0171-1997)要求进行工作。
1、快速静态GPS测地,测区范围及自然地理概况;
1.1测区范围;测区位于内蒙古锡林浩特盟苏尼特右旗和二连浩特市之间。测区范围为:E 111°48.25′~112°24.36′,N 42°55.31′~43°19.56′。
测区坐标系:WGS84坐标系。测区在1个19带的六度带内,其中央子午线为111度。
1.2测区自然地理状况;二连浩特市位于中国正北方内蒙古自治区中部,与蒙古国扎门乌德市隔界 相望。二连浩特地势平坦,由西南向东北缓缓倾斜,平均海拔为932.2m。
苏尼特右旗位于内蒙古自治区中部,锡林郭勒盟西部,是锡盟的西大门,东邻苏尼特左旗、镶黄旗; 南靠乌兰察布市察哈尔右翼后旗、商都县;西接乌兰察布市的四子王旗;东北与本盟二连浩特市接壤。苏 尼特右旗属古湖盆地上升而形成剥蚀高原,平均海拔高度为1000m~1400m。测区主干道为国道208,其他 辅道均为草原沙石路。
2、观测方案与数据采集;
2.1观测基准,数据采集的坐标系统为WGS-84坐标系,其参考椭球几何参数为:长半径: a=6378137m;短半径:b=6356752.3142m;扁率:α=1/298.257223563。
2.2观测模式及数据采集,本次采用快速静态观测模式,采样率为1s,观测过程中GPS接收机不关 机,在测量点位上观测3~5分钟,然后直接移动接收机到下一个测点,为了控制解算精度,每隔10个测点 进行一次15分钟的加密观测;为了避免因时差导致的数据分割致使无法解算,在每天上午8点时刻同样需 要进行加密观测;同时在电池更换后的测点上同样也要进行加密15分钟观测,直至一天工作日的测量工作 完成后再关闭接收机。为了避免大功率发射塔、高压线等影响,测地人员在遇到高压线时适当的避开,若 无法避免或遇到发射塔时则同样延长测点的观测时间。
2.3检查点观测,为了检查普通测点的观测精度,每天对一部分的测点进行一次重复观测,检查点 一般不少于当天工作量的5%。
3、GPS内业计算;野外观测完毕后,当天需将数据从接收机中下载出来,利用测绘专业软件,将接 收机格式数据转换成标准格式。并利用软件对生成的文件进行必要的修改,如天线高,天线类似,初始坐 标,接收机等信息。修改完毕后再进行数据的质量检查,并进行粗略的单点定位计算,从而确保收集的数 据没有问题。
GPS数据处理流程:
(1)整理原始观测数据:将原始观测数据转换成标准RINEX,利用软件检查数据质量并编辑完善 RINEX头文件信息。
(2)获取IGS的高精度已有产品:在IGS网站上获取相应观测日的精密卫星轨道、精密卫星时钟改 正、卫星掩蔽信息及有关地球的时间、极移、自转和极轴方向等资料。
(3)选择处理策略:采用高效的精密单点处理(PPP)模式,固定JPL精密卫星轨道和钟差,估计 接收机钟偏(白噪声)和每5分钟估计对流层天顶梯度延迟,对流层模型采用GMF,高度截止角15°,卫星 和接收机天线采用绝对相位中心变化和偏移模型,海洋潮汐负载模型采用FES2004。
4、完成的工作量及精度评价;GPS测量导线总长度105.4km,有效物理点位1058个,实测物理点 1165个,检查点107个。
精度评定采用公式
Figure BDA0002248147660000111
其中p1、p2分别是普通点和检查点的测量 点位,n代表点的个数。经平差后平面点位误差为0.18m;高程点位误差为0.14m。本次静态GPS测量工 作,依靠先进精密的测量仪器、科学施工、严细管理,圆满完成了任务。精度可靠,各项指标均达到规范 中的相关要求。高程测量结果见图1-图4。
5、重力测量;
1重力仪性能,按照中华人民共和国地质矿产部行业标准《大比例尺重力勘探规范》要求,在进行 野外数据采集之前要对仪器各种性能指标进行测试。本次实验包括的内容有静态试验和动态实验。
1.1静态试验,本次静态试验地点为江西省南昌市东华理工大学地学楼,该试验室内有专门设置的加 固水泥墩来避免楼层震动对仪器的干扰。观测时间超过了规定的24小时;重力仪每隔100秒观测一次,仪 器本身带有固体潮改正软件,根据改正后的数据绘制出重力仪的静态漂移曲线如图5。从图中可以看出, 仪器所测曲线大致呈线性变化,掉格系数K约为-0.016mGal/h。
1.2动态试验,试验地点选在高度不同的9个点之间进行,重复观测13次,连续观测时间约为4小 时,试验点之间重力差大于0.5×10-5/s2,两点间单程观测时间均小于20分钟。由于仪器内置固体潮改正, 因此可直接绘制出各点的动态零点位移曲线图6。图中可以看出,仪器在各点的动态零点位移曲线表现出 明显的相似性,均大致呈线性递增变化。
仪器动态掉格系数公式为:
Figure BDA0002248147660000112
(单位:mGal/min)
其中:Δgi为第i点前后观测值之差;Δti为第i点前后观测时间之差,单位:分钟;n为重复观测 点数。根据动态观测数据结合上式可计算出重力仪动态零点位移率k=-0.0004mGal/min。
对动态观测精度进行统计,按照独立增量的形式计算出仪器在9个不同高度点之间连续观测的增量。 仪器动态观测精度评价公式为:
Figure BDA0002248147660000121
(单位:10-5m/s2)
其中:Vi为第i个相邻两点间单个独立增量与平均增量之差;m为增量总数;n为联测的边数。本次 试验m=104,n=8,由此计算得到重力仪的动态试验均方差为0.0066mGal。
由CG-5型重力仪的静态试验和动态试验结果可以看出,本次野外重力勘探采用的重力仪完全符合高 精度测量要求。
2重力数据采集及整理;
2.1重力基点及观测,由于本次重力测量工作采样间隔为0.1km,四条测线长度均小于100km(最长测 线为39.9km),因此可以采用相对重力测量模式进行操作。本次野外共设立了两个重力基点,两个基点均 选择在交通方便,标志明显的公路旁边,并用约1m长的钢筋竖直置于地下以便易于保存。为保证基点的精 度,本次测量中早基和晚基均采用“基-辅-基”的模式进行,并使两次在基点上测量误差不大于0.005mGal。
2.2普通点与检查点观测;研究区普通点重力测量采用单次观测模式,每个点观测时间为60s(每秒 仪器自动观测一次数据)。重力仪操作者在60s需观察数据变化,若数据变化过大(一般认定连续五个数据 的变化均大于0.01mGal),则需重新进行观测以保证测点的数据精度。
为了评价测量的准确性和精度,需对测量过的普通点进行一定比例的检查测量。本次野外四条重力勘 探线的检查点数均满足5%~10%的要求,且较短测线的检查点数大于10%,检查点基本均匀分布在整条测线 上。
3重力数据整理;由于研究区属于地势平坦的高原区,因此可以忽略地形对重力值的影响,局部地势 变化略大的地区采用目估地形改正值模式。因此,本工区的重力数据整理主要包括正常场改正(纬度改正)、 中间层改正和高度改正三项。
3.1正常场改正及精度,在我国,正常场的计算采用的是1901-1909年赫尔墨特正常重力公式:
Figure BDA0002248147660000122
(单位:mGal),式中
Figure BDA0002248147660000123
为测点纬度。
正常场改正均方误差由下式计算:
Figure BDA0002248147660000124
(单位:mGal),其中
Figure BDA0002248147660000125
为测点平均 纬度,εD为测点纵坐标均方误差(单位:km)。
取测点的平均纬度为43.1°,εD=±0.001km(由测量数据得到),经计算εD=±0.0008mGal,完全 满足中华人民共和国地质矿产部行业标准《大比例尺重力勘探规范》规定0.01mGal的要求。
3.2布格改正及精度,高度改正采用如下公式:
Figure BDA0002248147660000131
(单位:mGal),其中
Figure BDA0002248147660000132
为测点的纬度,h为测点海拔高程,高程数据由测地测量工作提供。
中间层改正按下式计算:Δgσ=-0.0419{σ}g/cm 3{h}m(单位:mGal);其中σ为中间层密度, 依据研究区岩石密度资料选σ=2.65g/cm3
布格改正误差(不考虑中间层密度误差)用下式计算:
Figure BDA0002248147660000133
式中εh为测点高程均方误差(由测量提供,εh=0.14m), 经过计算得:εb=±0.030mGal。
3.3相对布格重力异常值的计算,布格重力异常值采用下式计算:
Figure BDA0002248147660000138
式中:gk为测点相对重力值(单位:mGal);Δgh为高度改正 值(单位:mGal);Δgσ为中间层改正值(单位:mGal);
Figure BDA0002248147660000137
为相对正常重力值(单位:mGal)。
3.4,重力测量质量评价;按《重力勘探技术规程》的要求,要采用等精度、重复观测的方法进行质量 检查和评价。本次野外测量总计完成115个检查点,占总测点数的10.9%。测点重复观测精度评价公式为:
Figure BDA0002248147660000134
其中:di为第i点前后观测值之差;n为重复观测点数;经计算得εg=0.014mGal,同样满足规范要求 的0.03mGal。
对所获得的布格重力异常的精度评价要考虑到观测误差、布格改正误差、正常重力场改正误差。因此 布格重力异常精度要依下式进行评价:
Figure BDA0002248147660000135
前述已给出εg=±0.014mGal;εb=±0.030mGal;
Figure BDA0002248147660000136
经计算得ε=±0.036 mGal。满足《大比例尺重力勘查规范》要求的0.05mGal。
4重力测量完成的工作量情况及观测点精度评价,重力测量结果无任何跳点,且数据变化特征与测点 高程之间的镜像关系明显,这也相互验证了本次快速静态GPS测地和重力测量结果都是准确的。通过正常 场改正、布格改正得到的相对重力异常相对于相对重力值来说,异常形态与测点高程之间无明显相关性, 说明重力资料整理是正确的;同时异常曲线更为圆滑、简洁,所反映的“重力高”与“重力低”特征更为 明显。
三、采用的方法技术,为说明书内容。
四、重力剖面数据处理及资料解释;
第一:四条重力剖线的重力场特征;
1.1NL线重力剖面,NL线剖面主体走向EW,剖面主要位于格日勒敖都凹陷内,仅在剖面末端(东 部)处于东方红凸起上。从NL线相对布格重力异常(图22)中可以看出,重力异常呈现“西低东高”的 特征,在剖面西部,重力异常变化平缓,且呈现为相对弱负异常;在10-19km处异常呈明显的上升变化, 异常梯度变化可达每公里1.48mGal,格日勒敖都凹陷和东方红凸起的界限恰处于这个异常梯度带中。
1.2MD线重力剖面,MD线剖面主体走向为NE向,剖面主体位于东方红凸起上,NE段位于格日勒敖 都凹陷中。从相对布格重力异常(图23)可以看出,重力异常呈现出西南高、中部低、东北高的2个重力 高和1个重力低的特征,且重力异常的高低与剖面经过的凸起和凹陷关系密切。
1.3E400线重力剖面,E400线走向为ES,剖面西北段处于东方红凸起上,中段处于齐哈日格图凹陷 中,东南段处于苏尼特隆起上。从重力异常图(图24)中可以看出,该线异常在-1.43~20.71mGal间变化, 异常幅度变化大;整体上可将重力异常分为3个重力高和2个重力低,其中西北段的重力高很可能是由东 方红凸起引起的;东南端的重力高则与苏尼特隆起相关;而中部的两个重力低之间的重力高则可能说明齐 哈日格图凹陷中部还存在着局部的次级凸起。
1.4SH线重力剖面,SH线走向为135°,与E400线平行。从图25可以看出该剖线全部处于脑木根 凹陷之中,但从相对重力异常中可知该剖线两端重力低,中段重力高,整体上呈现一个重力高的特征。这 说明在闹木根凹陷中还存在着次一级的凸起。
第二:四条剖线断裂构造识别的重力异常数据处理:
2.1一阶导数分析,利用重力异常水平导数的极值点和垂向导数的梯度带(或零点)可以确定断层 的水平位置。本次采用前面提及的波数域迭代法计算了4条剖面重力异常的水平和垂向一阶导数,结果见 图26-41。图中可以看出,相对于原始重力异常(图22-37)来说,导数计算显著地提高了异常的横向分辨 率,同时结果表现出了较强的计算稳定性,说明利用迭代法进行实际数据的导数转换是可行的。
2.2导数标准差的tilt angle法(TASD),导数标准差的tilt angle法是编写人博士阶段的一个成 果,该方法是通过计算结果的极大值来识别地质体边界的,不仅具备一定的计算稳定性,同时还可以同时 识别出不同深度场源的边界。四条重力剖面的TASD结果见图30-45。
2.3基于欧拉反褶积的断裂反演,欧拉反褶积是以场位和位的导数满足欧拉齐次方程为理论基础的 一种反演方法。该方法能在先验信息较少的情况下,自动或半自动地反演出场源质心位置、边界位置、场 源深度等几何参数;由于该方法具有较强的灵活性和实用性,目前已成为了地球物理工作者常用的数据处 理反演方法和重磁资料解释工具。图34-49是基于波数域迭代求导基础上的欧拉反褶积反演结果。图中可 以看出,欧拉反演点可以较好地识别出断裂的位置和反演出断裂的产状。
2.4 4条剖面断层处理结果综合分析,为了进一步分析各种方法处理的结果反映的断层断点及产状 的准确性,将导数分析结果、TASD处理结果和欧拉反褶积反演结果叠置在一起,可获得4条剖面的综合分 析结果(见图38-53)。图中可以看出,不同方法反映的断层断点位置具有较好的一致性。由于欧拉反演点 是滑动窗口下最小二乘反演解,反演点分布较多的地方可认为是断裂的走向位置,因此笔者认为当三种方 法对断层断点及倾向判定存在差异性时,应该以欧拉反褶积反演结果为标准。
通过对重力异常的一阶导数分析、断裂识别的TASD法及欧拉反褶积反演结果的综合分析,可在NL 重力剖面上划分出5条断裂,MD剖面上可划分出6条断裂,在E400剖面上可划分出9条不同规模的断裂, 在SH线剖面上划分出4条断裂。
第三,4条剖面的重力场分离,由布格重力异常的地质地球物理意义可知,在地面上测量的重力值 经过各项改正后,所获得的布格重力异常是纵向和横向各种地质体与围岩的密度差所引起的各类异常的叠 加。即布格重力异常具有一定的复杂性,异常的分离也是一个较为复杂的问题。
若要研究沉积地层分布问题,就需要采用某种滤波方法获取这部分地质体因密度差所引起的重力异 常。这里采用迭代滤波法对四条剖面进行重力场的分离,结果见图42-57。图中可以看出,区域重力异常 反映的是原始重力异常的整体变化趋势,表明了分离出的区域场是合理的;局部(剩余)重力场正、负异 常相间出现,正异常和负异常幅度相当,同时正异常与负异常之和趋于零值,这证实了利用迭代滤波法分 离出的局部重力异常恰是浅部地质体产生的。
第四、4条剖线剩余重力异常相关概率成像,重力数据的相关成像是根据概率成像的概念发展而来 的,并经过了前人的模型和实际数据的试验证明了该方法在一定程度上能够勾画出重磁异常的场源分布, 并且具有计算效率高,计算过程稳定,容易实现等特点。
为了揭示4条剖面地下地质体的赋存状态,对4条剖面的剩余重力数据进行了归一化互相关概率成 像处理,结果见图46-61。图中可以看出,NL线呈现“两凸一凹”特征,“凹陷区”大致10km宽,根据 EZK1248-1903#(,工业铀矿井,未到基底)和EZK1420-2031#(无矿井,到基底)两个钻口资料,可证实 NL剖线的概率成像反演结果是合理的;MD线呈现“两凸两凹”特征,左侧凹陷区处于该线西段边部,中部 隆起区范围达16km,中部凹陷宽度约10km,根据EZK1376-1903#(矿化井,未到基底)和EZK1420-2031# (到基底)两个钻口资料,同样可证实MD剖线的概率成像反演结果是合理的;E400线概率成像结果呈现 “三凸两凹”的特征,虽然该条剖面编写者未能收集到相应地区的钻孔资料,但是从电性资料中可以验证 反演结果的可靠性;SH线主体处于一个凸起之上,两侧呈现出凹陷的迹象,DHZK-1#(无矿井,到基底) 的资料表明了重力资料反演结果的有效性。
第五、四条重力剖面地质-地球物理模型的建立及综合解释,在进行重力剖 面人机交互反演之前,需要收集研究区的岩石密度资料、过测线钻孔资料、地质 资料以及其他地球物理资料等相关内容。通过与其他单位的沟通,四口钻孔资料、 以及E400线CSAMT电性资料,这些重要的资料为重力正演模拟的可靠性提供了 有利的保证。依据相关资料,将正演模拟分为四层,分别是第四纪地层、第三纪 地层、古生代-中生代地层以及基底,分别对应的平均密度为2.05g/cm3、 2.15g/cm3、2.45g/cm3、2.75g/cm3,并取岩石平均密度为2.65g/cm3。在正演模 拟过程中,地层分布厚度变化则主要依据剩余重力异常相关成像结果。四条重力 剖面的地质-地球物理模型分别见图52、55、58和图61。
第六,NL剖面地质-地球物理综合解释,该线布设于格日勒敖都凹陷中,测线经过EZK1248-1903工 业铀矿井和EZK1420-2031无矿井。该剖面的布设目的在于验证重力勘探是否可用于寻找古河道,并利用重 力资料研究格日勒敖都凹陷的次级构造和沉积特征。
结合剩余重力异常(图50)、断裂分布(图51)和构建的地质-地球物理模型(图52),将NL剖线从 西至东依次划分为格日勒敖都凹陷中部次级凸起、格日勒敖都凹陷东部次级凹陷和东方红凸起。
6.1.1 NL剖面重力场与断裂构造关系,图51是依据重力异常导数分析、TASD、欧拉反褶积反演结 果勾勒出的5条不同规模的断裂。其中F1是格日勒敖都凹陷次级凸起中的一条断裂,在(剩余)重力场中 反映不明显,仅以重力梯级带的突然变化为表现特征,但导数分析、TASD和欧拉反褶积反演结果(图16) 对该断裂均有较为明显的反映;F2为格日勒敖都凹陷中部次级凸起与东部次级凹陷之间的界限断裂,在重 力场中表现为明显的梯级带特征;F3处于格日勒敖都凹陷东部次级凹陷中,在重力场中表现为缓慢变化的 梯级带特征,即表明了该条断裂倾角较小(与反演结果吻合);F4同样处于格日勒敖都凹陷东部次级凹陷中, 在重力场中无明显显示,但三种断裂识别方法同样表明了该断裂是存在的,在此需要指出的是,在重力野 外工作过程中,由于项目组成员观察仔细,发现该条断裂出露于地表;F5断裂是格日勒敖都凹陷与东方红 凸起的界限断裂,在重力场中表现为大型的重力异常梯级带,该接触带断裂已被EZK1420-2031钻孔予以证 实。
6.1.2 NL剖面重力场与沉积地层的分布关系;图52是依据剩余重力异常互相关成像结果(图46), 利用正演模拟得到的地层分布关系图。从研究区地质图中可以看出,该剖面出露于地表的地层为第三系, 无第四纪地层。从地质-地球物理模型中可以看出:第三纪地层厚度有限,整体上在格日勒敖都凹陷中部次 级凸起和东方红凸起上厚度较大,这也恰于不同年代的地层分布位置吻合,第三纪地层最大可解释深度为 200m,位于剖线的14.5km处;在新生代地层之下有着一定规模分布的中生代地层(古生代地层存在的可能 性较低,其原因有二:1、二连盆地古生代地层仅在局部地区有分布;2、EZK1420-2031钻孔见基底但未见 古生代地层),该时期的地层主要分布在格日勒敖都凹陷东部次级凹陷之中,具有一定的厚度和规模,且分 布具有明显的河道沉积特性,可推断“河道”宽度约为10km,底界面平均深度约为730m;本条剖面的基底 呈现为东、西两侧埋深浅、中部埋深大的特点,模拟的基底最浅深度仅为230m,处于测线末端。
6.2 MD剖面地质-地球物理综合解释,该线从西南到东北依次布设于东方红凸起上、格日勒敖都凹 陷中和东方红凸起上,测线经过EZK1376-1903矿化井和EZK1420-2031无矿井。该剖面的布设目的与NL线 相同。
依据剩余重力场异常形态(图53)、断裂展布特征(图34)和构建的地质-地球物理模型(图55), 将MD剖面从西南至东北依次划分为东方红凸起、格日勒敖都凹陷东部次级凹陷和东方红凸起。
6.2.1 MD剖面重力场与断裂构造关系,图54是依据重力异常导数分析、TASD、欧拉反褶积反演结 果勾勒出的6条不同规模的断裂。其中F1断裂位于东方红凸起中,在重力场中表现为明显的异常梯级带特 征;F2、F3同样位于东方红凸起中,在重力场中同样表现为异常梯级带特征,只是异常梯度大小不同;F4断裂为东方红凸起和格日勒敖都凹陷之间的界限断裂,该断裂在重力场中同样呈现为异常梯度变化特点; F5处于格日勒敖都凹陷之中,在重力场中无明显显示;F6同样是东方红凸起和格日勒敖都凹陷之间的界限 断裂,在重力场中以异常梯级带展示。
6.2.2 MD剖面重力场与地层分布关系,图55是依据剩余重力异常互相关成像结果(图47),利用 正演模拟得到的地层分布图。从研究区地质图中可以看出,该剖面出露于地表的地层为第三系,同样无第 四纪地层。第三纪地层沉积厚度有限,整体分布较为均匀,平均厚度约为130m,在东方红凸起上沉积较厚, 最大可解释厚度为190m;在新生代地层之下同样有着一定规模分布的中生代地层(古生代地层存在的可能 性较低,其原因与NL剖面相同),该时代的地层在东方红凸起中分布厚度变化不大,平均厚度约为210m, 在格日勒敖都凹陷东部次级凹陷中厚度较大,且具备一定的规模,分布同样具有河道沉积特性,可推断“河 道”宽度同样为10km,底界面平均厚度为700m;本条剖面的基底呈现为西南、东北两侧埋深浅、中部埋深 大的特点,模拟的基底最浅深度仅为260m,处于测线末端。本剖面由于与NL剖面在地理位置上相距较近, 因此岩石分布特征具有明显的可对比性。
6.3 E400线地质-地球物理综合解释,该线从西北到东南依次布设于东方红凸起上、齐哈日格图凹 陷中和苏尼特隆起之上。该剖面布设位置与2012年电法E400线位置完全重合。该剖面的布设目的在于探 究重力、电法的反演结果是否具备一致性或相似性,并进行重、电联合反演以提高地质模型的可靠性。
依据乌兰察布坳陷地质图,结合剩余重场异常特征(图56)、断裂展布(图57)和构建的地质-地 球物理模型(图58),将E400剖线从西北至东南依次划分为红方红凸起、齐哈日格图凹陷(次级构造依次 为:西部次级凹陷、中部次级凸起和东部次级凹陷)和苏尼特隆起。
6.3.1 E400剖面重力场与断裂构造关系,图57是依据重力异常导数分析、TASD、欧拉反褶积反演 结果勾勒出的9条不同规模的断裂。其中F1是东方红凸起与齐哈日格图凹陷的界限断裂,在(剩余)重力 场中表现为大型异常梯级带特征;F2为齐哈日格图凹陷东部次级凹陷中的一条断裂,在重力场无明显反映; 但三种断裂识别方法(导数分析、TASD、欧拉反褶积)均有较明显的反演结果,该断裂可能是沉积地层之 间的接触带或者地层变化较明显的位置;F3是齐哈日格图凹陷东部次级凹陷与中部次级凸起的界限,在重 力场中同样以明显的梯级带为反映特征;F4为齐哈日格图凹陷中部次级凸起与东部次级凹陷的界限断裂, 重力场中同样以明显的异常梯级带为反映;F5断裂是齐哈日格图凹陷东部次级凹陷内部的一条断裂,该断 裂在重力场中以异常梯级带为反映,由于梯级带范围较宽,可断定该条断裂具有一定的走向和规模;F6和 F8分别处于齐哈日格图凹陷和苏尼特隆起中,在重力场中均无明显反映;F7则是齐哈日格图凹陷与苏尼特 隆起的界限断裂,在重力场中表现为大型异常梯级带;F9是苏尼特隆起中的一条小型断裂,但在重力异常 中也表现出了明显的异常梯级带。
6.3.2 E400剖面重力场与沉积地层的分布关系,图58是依据剩余重力异 常互相关成像结果和MT视电阻率反演结果(图47),利用正演模拟得到的地层 分布图。从研究区地质图中可以看出,该剖面出露于地表的地层为第三系和第四 系。从地质-地球物理模型中可以看出:第四纪地层分布有限,且厚度很小,仅 在3.5~11.2km之间分布,平均厚度为65m,最大厚度为150m,主要分布在齐哈 日格图凹陷西部次级凹陷中;第三纪地层仍是近似水平层状均匀分布,平均厚度 为75m;在新生代地层之下有着一定规模分布的古生代和中生代地层(由于该测 线南侧有局部的早二叠世地层出露,因此古生代地层可能存在于本剖面之下), 古、中生代地层主要分布在齐哈日格图凹陷西部、东部次级凹陷之中,具有相当的厚度和规模,这两个次级凹陷的沉积地层可解释深度均大于1000m;本剖面基 底起伏特征明显,埋深较浅的位置位于东方红凸起、齐哈日格图凹陷中部次级凸 起和苏尼特隆起上,在这三个位置的最浅深度分别为80m、170m和150m。
6.4 SH剖面地质-地球物理综合解释,该测线是项目组在野外工作过程中临时添加的一条重力剖面, 剖面布设于闹木根凹陷中,其布设目的主要在于验证古河道是否存在于此,同时也为了研究闹木根凹陷地 下岩石分布特征。
依据乌兰察布坳陷地质图,结合剩余重场异常特征(图59)、断裂展布(图60)和构建的地质-地 球物理模型(图61),将SH剖线从西北至东南依次划分为闹木根凹陷西部次级凹陷、中部次级凸起和东部 次级凹陷。
6.4.1 SH剖面重力场与断裂构造关系,图60是依据重力异常导数分析、TASD、欧拉反褶积反演结果 勾勒出的4条断裂。其中F1、F2处于闹木根凹陷西部次级凹陷之中,在重力场中均表现为缓慢变化的异常 梯级带,在导数分析、TASD和欧拉反褶积中反映明显;F3是闹木根凹陷西部次级凹陷与中部次级凸起的界 限断裂,以异常梯级带在重力场中予以展示;F4是闹木根凹陷中部次级凸起与东部次级凹陷的界限,在重 力场中以明显的异常梯级带为反映特征。
6.4.2 SH剖面重力场与岩石分布关系;图61是依据剩余重力异常互相关成像结果(图48)正演模 拟得到的0~1500m的岩石分布图。从研究区地质图中可以看出,该剖面出露于地表的地层为第三系。从地 质-地球物理模型中可以看出:第三纪地层同样是近似水平层状均匀分布,平均厚度为80m;在新生代地层 之下有着平均厚度为420m的中生代地层(古生代地层存在的可能性较低,其原因与NL剖面大致相同),该 时代的地层并不像其他三条剖面那样起伏较为明显,地层底界面在闹木根凹陷东部次级凹陷中埋深较深, 最深处为760m,在闹木根凹陷中部次级凸起中埋深较浅,最浅处为370m;基底起伏变化同样不明显,只是 在闹木根凹陷中部次级凸起中有一定的隆起特征。
基于剩余重力异常相关成像所建立的地质-地球物理模型产生的重力异常与已知的剩余重力异常拟 合较好,且模拟地层分布与已知的地质资料和钻口资料吻合。这表明:构建的地质模型具备一定的可信度, 这为研究二连盆地沉积地层分布、古河道位置和隐伏隆起等重要的地质信息提供了有力的依据。

Claims (7)

1.一种利用地面高精度重力测量的找矿方法,其特征在于:主要包括六个方法,分别为位场高阶导数的波数域迭代法、基于迭代滤波法的位场分离方法、识别地质体边界的TASD法、位场数据处理的欧拉反褶积法、位场异常互相关成像技术方法和重力人机交互剖面正演模拟方法。
2.根据权利要求1所述的一种利用地面高精度重力测量的找矿方法,其特征在于:所述位场高阶导数的波数域迭代法;具体方法为:
在波数域中,位场数据的任意阶任意方向导数波谱可以表示为:
Figure FDA0002248147650000011
其中S(m,p,q)(u,v)和S(u,v)分别为所求导数波谱和原始异常波谱;
Figure FDA0002248147650000012
分别为x方向m阶导数算子,y方向p阶导数算子和z方向q阶导数算子;u、v分别是x、y方向的波数;m、p、q均为大于等于零的实数;
由上式可知,求导阶次l和理论导数算子
Figure FDA0002248147650000013
分别为:
l=m+p+q,
Figure FDA0002248147650000014
Figure FDA0002248147650000015
带入原始公式中,则有:
Figure FDA0002248147650000016
对S(m,p,q)乘以一个包含求导阶次l和导数算子
Figure FDA0002248147650000017
相关的低通滤波器P,即:
Figure FDA0002248147650000018
α≥1,β>0,
Figure FDA0002248147650000019
Figure FDA00022481476500000110
的绝对值数值形式,经P滤波后的S(m,p,q)近似值可表示为:
Figure FDA00022481476500000111
Figure FDA00022481476500000112
带入原式中,得S的一阶近似值:
Figure FDA00022481476500000113
再用S和S(1)的差值对
Figure FDA00022481476500000114
进行校正,得S(m,p,q)的二阶近似值:
Figure FDA00022481476500000115
重复以上步骤,得S(m,p,q)的n阶近似值:
Figure FDA0002248147650000021
且有:
Figure FDA0002248147650000022
则得:
Figure FDA0002248147650000023
Figure FDA0002248147650000024
分别是第n、n-1次迭代得到的波谱
Figure FDA0002248147650000025
的傅里叶反变换,ε为给定的较小正实数;当迭代过程中满足
Figure FDA0002248147650000026
时终止迭代,此时可认为所求导数
Figure FDA0002248147650000027
3.根据权利要求1所述的一种利用地面高精度重力测量的找矿方法,其特征在于:所述基于迭代滤波法的位场分离方法;具体方法为:
令叠加异常u(x,y)的波谱为U(u,v),假设算子H(u,v,α)是一低通滤波器且满足0≤H≤1,那么将H(u,v,α)作用于U(u,v)可以达到场分离的目的,即有:
Ureg(u,v,α)=U(u,v)·H(u,v,α)
这里u、v分别是x、y方向上的波数,α是算子H的滤波参数;
按照滤波参数α的选取可以将场分离情况分为以下三种:
①滤波参数α选取不当使得分离出的区域场ureg中含有局部场成分,其中ureg=F-1[Ureg],F-1是Fourier反变换;
②滤波参数α选取恰当使得叠加异常u分离较为彻底;
③滤波参数α选取不当使得分离出的局部场ured=u-ureg中含有区域场成分;
在实际应用中,由于滤波参数α不易确定而使得场的分离容易出现第一种或第三种情况,使得位场分离效果不理想;
针对第一种情况,提出改进方案:固定滤波参数α使得区域场ureg含有较多的局部场成分,那么将H再次作用于Ureg以进一步剥离出局部异常,即有
Figure FDA0002248147650000028
倘若
Figure FDA0002248147650000029
还含有局部场成分,则可以继续将H作用于
Figure FDA00022481476500000210
依次类推,那么H作用n次的区域场波谱为:
Figure FDA0002248147650000031
只要选合适的作用次数n,即可以使得区域场中的局部场成分较好的剥离出去;
针对第三种情况,改进方案为:固定参数α使得局部场ured含有较多的区域场成分,那么此时需要将局部异常波谱中的区域场成分剥离出来“返还”给区域场,即
Figure FDA0002248147650000032
此时局部异常波谱变为:
Figure FDA0002248147650000033
若局部场中仍有区域场信息,则继续对局部场进行剥离并“返还”给区域场,如此一直进行下去,至到剥离“返还”n次时,满足局部场中没有或几乎没有区域场成分时终止,此时,区域异常波谱
Figure FDA0002248147650000034
和局部异常谱
Figure FDA0002248147650000035
分别为:
Figure FDA0002248147650000036
Figure FDA0002248147650000037
上述第二种改进方案中,局部场剥离和“返还”过程与迭代滤波的思想模式是一致的,这种改进方案为分离场的迭代滤波法;
第二种改进方案中给出低通滤波H形式为:
Figure FDA0002248147650000038
为了满足第二类情况的条件,设定α≥1,h≥50Δd,β≥1,而迭代次数即剥离“返还”次数的选取方案采用的是差值互相关系数之差。
4.根据权利要求1所述的一种利用地面高精度重力测量的找矿方法,其特征在于:所述识别地质体边界的TASD法;具体方法为:
导数标准差的归一化形式进行场源边界的识别,其归一化形式为
Figure FDA0002248147650000039
上述方法存在着识别地质体边界连续性差和易受噪声干扰的缺点,提出导数标准差的倾角法,其形式如下:
Figure FDA0002248147650000041
其中改造因子m≥1。
5.根据权利要求1所述的一种利用地面高精度重力测量的找矿方法,其特征在于:所述位场数据处理的欧拉反褶积法;具体方法为:
Euler齐次方程如下:
Figure FDA0002248147650000042
其中,f(x,y,z)是位场异常;
Figure FDA0002248147650000043
分别是异常f在x、y、z三个方向上的一阶导数值;N为构造指数,不同的数值代表不同的地质体形状;当N=0时,上式可用来反演断裂分布点。
6.根据权利要求1所述的一种利用地面高精度重力测量的找矿方法,其特征在于:所述位场异常互相关成像技术方法;具体方法为:
在地形起伏的测区,笛卡尔坐标的(x,y)平面在基准面上,z取垂直向下为正,假设测区地下任意坐标为(xq,yq,zq)、体积为vq的第q个点质量的剩余密度为Δσq,则它在测区上任意点(x,y,z)处的重力异常Δgq(x,y,z)可表示为
Figure FDA0002248147650000044
其中,万有引力常数G=6.67×10-11m3/(kg·s2);
定义测区上实测重力异常与第q个点质量在测区上的重力异常的归一化互相关为:
Figure FDA0002248147650000045
其中,(xi,yi,zi)为测区第i个观测点的坐标,Δg(xi,yi,zi)为该观测点的实测重力异常,Δgq(xi,yi,zi)为地下第q个点质量在该观测点的重力异常,N为测区观测点总数。
7.根据权利要求1所述的一种利用地面高精度重力测量的找矿方法,其特征在于:所述重力人机交互剖面正演模拟方法;具体方法为:
(1)将重力异常的处理、转换方法得出的结果和其它地质、地球物理方法获取的先验信息输入到模型里,形成初始模型;
(2)根据计算结果和实际重力异常的差异,随时方便地修改模型,直观地监督和指导正反演过程;
(3)通过对初始模型计算出的重力效应与测线上的剩余重力异常进行对比,不断修正模型,直至达到计算出的重力效应与测线上的剩余重力异常之差满足预定精度。
CN201911027784.6A 2019-10-25 2019-10-25 一种利用地面高精度重力测量的找矿方法 Active CN110888176B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911027784.6A CN110888176B (zh) 2019-10-25 2019-10-25 一种利用地面高精度重力测量的找矿方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911027784.6A CN110888176B (zh) 2019-10-25 2019-10-25 一种利用地面高精度重力测量的找矿方法

Publications (2)

Publication Number Publication Date
CN110888176A true CN110888176A (zh) 2020-03-17
CN110888176B CN110888176B (zh) 2021-05-07

Family

ID=69746489

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911027784.6A Active CN110888176B (zh) 2019-10-25 2019-10-25 一种利用地面高精度重力测量的找矿方法

Country Status (1)

Country Link
CN (1) CN110888176B (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111490829A (zh) * 2020-04-08 2020-08-04 威海北洋电气集团股份有限公司 信号动态调节方法和系统及其在光纤传感中的应用
CN111580174A (zh) * 2020-05-29 2020-08-25 中国地质科学院 一种基于帕德近似的重磁数据向下延拓方法
CN111856615A (zh) * 2020-07-30 2020-10-30 中国人民解放军61540部队 重磁异常场重构方法及场源点物性参数解算方法、系统
CN112651102A (zh) * 2020-11-27 2021-04-13 中国地质科学院地球物理地球化学勘查研究所 一种起伏地形下的位场全自动极值点深度估算方法
CN113421194A (zh) * 2021-06-04 2021-09-21 贵州省地质矿产勘查开发局 一种根据布格重力异常图像提取隐伏断层的方法
CN115292660A (zh) * 2022-08-04 2022-11-04 中国自然资源航空物探遥感中心 一种自适应迭代的位场分离方法、系统、设备及计算机可读存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2539838C1 (ru) * 2013-07-16 2015-01-27 Федеральное Государственное Унитарное Предприятие "Сибирский Научно-Исследовательский Институт Геологии, Геофизики И Минерального Сырья" Способ прогноза скрытого оруденения, связанного с гранитоидами
CN104866653A (zh) * 2015-04-29 2015-08-26 中国地质科学院矿产资源研究所 一种获取地下三维密度结构的方法
CN106646660A (zh) * 2015-10-30 2017-05-10 核工业北京地质研究院 一种砂岩型铀矿综合地球物理勘探方法
CN107748399A (zh) * 2017-09-12 2018-03-02 中国石油化工股份有限公司 利用重力界面反演识别山前带深部构造层方法
CN109839670A (zh) * 2017-11-29 2019-06-04 核工业北京地质研究院 一种热液型铀矿床基底界面反演方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2539838C1 (ru) * 2013-07-16 2015-01-27 Федеральное Государственное Унитарное Предприятие "Сибирский Научно-Исследовательский Институт Геологии, Геофизики И Минерального Сырья" Способ прогноза скрытого оруденения, связанного с гранитоидами
CN104866653A (zh) * 2015-04-29 2015-08-26 中国地质科学院矿产资源研究所 一种获取地下三维密度结构的方法
CN106646660A (zh) * 2015-10-30 2017-05-10 核工业北京地质研究院 一种砂岩型铀矿综合地球物理勘探方法
CN107748399A (zh) * 2017-09-12 2018-03-02 中国石油化工股份有限公司 利用重力界面反演识别山前带深部构造层方法
CN109839670A (zh) * 2017-11-29 2019-06-04 核工业北京地质研究院 一种热液型铀矿床基底界面反演方法

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
KUNPENG GE ET AL.: "Rock magnetic investigation and its geological significance for vein-type uranium deposits in Southern China", 《GEOCHEMISTRY, GEOPHYSICS, GEOSYSTEMS》 *
王宝仁等: "《重力磁法实验学习教学指导书》", 28 February 1994, 地质出版社 *
王彦国: "位场数据处理的高精度方法研究及应用", 《中国博士学位论文全文数据库 基础科学辑》 *
王彦国等: "位场高阶导数的波数域迭代法", 《物探与化探》 *
王彦国等: "重磁勘探方法在砂岩型铀矿中的应用", 《东华理工大学学报(自然科学版)》 *
王明等: "Tilt-Euler方法在位场数据处理及解释中的应用", 《物探与化探》 *
郭良辉等: "重力和重力梯度数据三维相关成像", 《地球物理学报》 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111490829A (zh) * 2020-04-08 2020-08-04 威海北洋电气集团股份有限公司 信号动态调节方法和系统及其在光纤传感中的应用
CN111490829B (zh) * 2020-04-08 2022-09-20 威海北洋电气集团股份有限公司 信号动态调节方法和系统及其在光纤传感中的应用
CN111580174A (zh) * 2020-05-29 2020-08-25 中国地质科学院 一种基于帕德近似的重磁数据向下延拓方法
CN111580174B (zh) * 2020-05-29 2023-08-11 中国地质科学院 一种基于帕德近似的重磁数据向下延拓方法
CN111856615A (zh) * 2020-07-30 2020-10-30 中国人民解放军61540部队 重磁异常场重构方法及场源点物性参数解算方法、系统
CN111856615B (zh) * 2020-07-30 2023-01-24 中国人民解放军61540部队 重磁异常场重构方法及场源点物性参数解算方法、系统
CN112651102A (zh) * 2020-11-27 2021-04-13 中国地质科学院地球物理地球化学勘查研究所 一种起伏地形下的位场全自动极值点深度估算方法
CN113421194A (zh) * 2021-06-04 2021-09-21 贵州省地质矿产勘查开发局 一种根据布格重力异常图像提取隐伏断层的方法
CN113421194B (zh) * 2021-06-04 2022-07-15 贵州省地质矿产勘查开发局 一种根据布格重力异常图像提取隐伏断层的方法
CN115292660A (zh) * 2022-08-04 2022-11-04 中国自然资源航空物探遥感中心 一种自适应迭代的位场分离方法、系统、设备及计算机可读存储介质

Also Published As

Publication number Publication date
CN110888176B (zh) 2021-05-07

Similar Documents

Publication Publication Date Title
CN110888176B (zh) 一种利用地面高精度重力测量的找矿方法
CN109828316B (zh) 一种钙结岩型铀矿找矿勘查方法
CN106875471A (zh) 煤系含或隔水层三维可视化建模方法
CN104267429A (zh) 确定地层压力的方法及装置
Li et al. Onshore–offshore structure and hydrocarbon potential of the South Yellow Sea
CN112711078B (zh) 沉积盆地深部有利砂岩型铀成矿砂体识别方法
CN111045091B (zh) 一种玄武岩覆盖下古河道的识别定位方法
CN107272081A (zh) 一种山前地区沉积相带展布范围预测方法
CN105137482A (zh) 一种沉积体古坡度的计算方法
CN113325486A (zh) 覆盖区下的构造蚀变岩型矿物勘查方法、系统及装置
CN103954996A (zh) 一种基于旅行时法确定地层裂隙裂缝走向的方法及装置
Jaffal et al. Gravity study of the Western Bahira Basin and the Gantour Phosphatic Plateau, central Morocco: Interpretation and hydrogeological implications
Qi et al. Hierarchy and subsurface correlation of muddy baffles in lacustrine delta fronts: a case study in the X Oilfield, Subei Basin, China
CN100349014C (zh) 重力勘探数据处理变密度地形校正方法
Arinze et al. A scalar-geometric approach for the probable estimation of the reserve of some Pb-Zn deposits in Ameri, southeastern Nigeria
Hanna et al. Gravity, magnetics, and geology of the San Andreas fault area near Cholame, California
Haefner et al. Evaluation of the horizontal-to-vertical spectral ratio (HVSR) seismic method to determine sediment thickness in the vicinity of the South Well Field, Franklin County, OH
Volkov et al. Use of the program and goal-oriented approach to observe the vertical displacements of the earth’s surface in Russia
Beshr et al. Using modified inverse distance weight and principal component analysis for spatial interpolation of foundation settlement based on geodetic observations
Indriani et al. Geological structure model for recharge area in Patuha Geothermal Field
Lutter 3D Gravitational Inversion Modelling of the South Dunedin Sub Basin, Otago, New Zealand
CN117572509B (zh) 一种与斑岩活动有关热液脉型矿产的找矿方法
Du et al. The stability analysis method of leveling datum points in mining areas of western china based on SBAS-InSAR technology
DuRoss et al. Paleoseismic investigation of the northern strand of the Nephi segment of the Wasatch fault zone at Santaquin, Utah
CN114814982B (zh) 预测花岗岩体铀矿有利成矿部位的方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant