一种提取局部重力异常的方法、系统、装置及介质与流程

1.本发明涉及航空重力异常数据处理以及矿产勘探技术领域，特别涉及一种基于二维曲率滤波法提取局部重力异常的方法、系统、装置及介质。


背景技术：

2.一般获得的航空自由空间重力数据经各种校正可得到航空布格重力异常数据，它是当代航空重力研究中获取最为简单且研究最为广泛的重力异常数据。航空布格重力异常是由区域重力异常和局部重力异常叠加组合而成。其中，由分布范围较大且影响较深的地质因素(例如：区域地质构造、莫霍面等)所引起的重力异常称为区域重力异常，由分布范围较小且影响较浅的地质因素(例如：金属矿产、盐丘等)所引起的重力异常称为局部重力异常。由此可见，局部重力异常的提取不仅是矿产资源勘探的重要研究内容，同时也有助于深部地质构造的研究。在此基础下，如何准确地提取出局部重力异常也就成为重要问题。
3.空间域非线性曲率滤波法于1987年从美国西方地球物理公司引进，在提取微弱位场异常中得到广泛应用，因是基于测线的位场异常分离方法，故本研究中将其称为一维曲率滤波法。
4.一维曲率滤波法是提取测线上尺度小、微弱的局部重力异常的方法，在航空重力异常数据处理和矿产资源勘探中得到实际应用。但是当获取较大尺度的局部重力异常时，该方法随着滤波窗口的增大会产生畸变现象。
5.现有技术实现了基于平面网格位场数据的空间域非线性曲率滤波法，思路是依次从八个方向进行一维曲率滤波，实现了对网格重力异常数据的局部重力异常的提取，但是这种方法在某些情况下可能存在计算偏差，而且该计算方式没有改变以测线为处理单元的思路，缺乏邻近测线之间联动。如果当第一次滤波出现畸变后，会把这种误差传递下去，进而影响到最终分离效果，导致获得局部重力异常存在偏差，造成数据的不准确，影响了在航空重力异常数据处理和矿产资源勘探中应用效果。


技术实现要素：

6.本发明针对现有技术普适性和准确性有待进一步提升的问题，提供了一种提取局部重力异常的方法、系统、装置及介质，根据平面重力异常的曲率变化非线性迭代逼近，逐步拟合出平面区域重力异常，从总场异常中减去区域重力异常得到局部重力异常。
7.为了实现以上发明目的，本发明采取的技术方案如下：
8.一种提取局部重力异常的方法，包括以下步骤：
9.s1：基于矿产勘探或地质构造研究中获得的航空空间重力异常，经布格校正和网格化处理得到航空布格重力异常网格数据。
10.s2：将航空布格重力异常网格数据作为第一次迭代计算的初始区域重力异常gf，估算给定滤波窗口宽度2wd内异常弯曲方向的特征值；
11.s3：判断异常弯曲方向是否受到叠加的局部重力异常的影响而发生交替变化，如
果是则更新区域重力异常值，如果不是则区域重力异常值保留原场值不变。
12.s4：在网格上依次从左至右、从上到下移动窗口，重复s2～s4步骤，直至整个网格计算完毕为止，计算得到的gf值为基于格网的区域重力异常的第一次迭代计算值；
13.s5：逐步减小滤波半窗口，每次减小滤波半窗口，就是在前一次区域重力异常迭代计算值gf的基础上，重复s2～s5步骤；以此逐次迭代逼近，直到滤波半窗口宽度减小为网格间距，最终迭代得到gf值就是基于网格的曲率滤波法计算获得的滤波窗口宽度为2wd的区域重力异常值；
14.s6：航空空间重力异常格网与最终迭代计算得到的区域重力异常格网做差，便得到局部重力异常。
15.进一步地，s2中用下式估算给定滤波窗口宽度2wd内异常弯曲方向的特征值：
16.skx1＝gf(kx2)-0.5[gf(kx1)+gf(k)]
[0017]
skx2＝gf(k)-0.5[gf(kx2)+gf(kx3)]
[0018]
skx3＝gf(kx3)-0.5[gf(k)+gf(kx4)]
[0019]
sky1＝gf(ky2)-0.5[gf(ky1)+gf(k)]
[0020]
sky2＝gf(k)-0.5[gf(ky2)+gf(ky3)]
[0021]
sky3＝gf(ky3)-0.5[gf(k)+gf(ky4)]
[0022]
式中，kx1、kx2、k、kx3、kx4为x方向上间隔距离为滤波窗口半宽度wd的测点位置，ky1、ky2、k、ky3、ky4为y方向上间隔距离为滤波窗口半宽度wd的测点位置，skx1、skx2、skx3为x方向上异常弯曲方向的特征值，sky1、sky2、sky3为y方向上异常弯曲方向的特征值。
[0023]
进一步地，s3中如果skx1、skx2、skx3同时满足skx1
·
skx2《0，skx3
·
skx2《0，或者sky1、sky2、sky3同时满足sky1
·
sky2《0，sky3
·
sky2《0，则可认为k点周围间距滤波窗口半宽度wd范围内的异常弯曲方向由于受到叠加的局部重力异常的影响发生了交替变化，由下面的式子更新k点的区域重力异常：
[0024]
gf(k)＝gf(k)-0.25(skx2+sky2)
[0025]
进一步地，s3中如果skx1、skx2、skx3不能同时满足skx1
·
skx2《0，skx3
·
skx2《0，或者sky1、sky2、sky3不能同时满足sky1
·
sky2《0，sky3
·
sky2《0，则认为k点周围间距滤波窗口半宽度wd范围内的异常弯曲方向没有受到叠加的局部重力异常的影响而发生交替变化，因此k点的区域重力异常值保留原场值不变；
[0026]
本发明还公开了一种提取局部重力异常的系统，包括：数据获取模块和局部重力异常提取模块；
[0027]
数据获取模块：获取航空空间重力异常数据处理或者矿产勘探中得到了航空布格重力异常网格数据，并输入至局部重力异常提取模块；
[0028]
局部重力异常提取模块：计算异常弯曲方向的特征值，判断异常弯曲方向是否受到叠加的局部异常的影响而发生交替变化，如果是则更新区域异常值，如果不是则区域异常值保留原场值不变。重复计算直至整个网格计算完毕为止，计算得到的值为格网区域重力异常的第一次迭代计算值；
[0029]
根据第一次迭代计算值，逐步减小滤波半窗口，重复计算上述直至滤波窗口减小到网格间距为止，获得区域重力异常值，数据获取模块中输入的航空重力异常格网数据与区域重力异常值做差，得到局部重力异常。
[0030]
本发明还公开了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述提取局部重力异常的方法。
[0031]
本发明还公开了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现上述提取局部重力异常的方法。
[0032]
与现有技术相比，本发明的优点在于：
[0033]
提高了在航空布格重力异常中提取局部重力异常的普适性和准确性，准确的分离结果可以提高矿产资源勘探中应用效果，也为地球深部构造反演研究提供更好的数据基础。
附图说明
[0034]
图1是本发明实施例二维曲率滤波法原理示意图；
[0035]
图中曲线为重力异常总场(航空布格重力异常)，虚线为区域重力异常，直线为重力异常总场曲线上横坐标为kx1、kx2、k等点的连线；
[0036]
图2是本发明实施例理论模型1，wd＝13km(130倍点距)，基于测线和网格的曲率滤波法分离效果对比图；
[0037]
图3是本发明实施例理论模型2，wd＝4km(40倍点距)，基于测线和网格的曲率滤波法分离效果对比图。
具体实施方式
[0038]
为使本发明的目的、技术方案及优点更加清楚明白，以下根据附图并列举实施例，对本发明做进一步详细说明。
[0039]
本发明的思路需要以同时考虑计算点周围测点的影响为前提。基于此提出了基于格网数据的改进方案，本发明的实质为根据平面位场异常的曲率变化非线性迭代逼近，逐步拟合出平面位场区域异常，本发明中将其称为二维曲率滤波法，其基本原理见图1所示，一种提取局部重力异常的方法，具体操作步骤如下：
[0040]
(1)基于矿产勘探或地质构造研究中获得的航空空间重力异常，经布格校正和网格化处理得到航空布格重力异常网格数据。
[0041]
(2)将原始重力异常网格数据作为第一次迭代计算的初始区域异常gf，用下式估算给定滤波窗口宽度2wd内异常弯曲方向的特征值：
[0042]
skx1＝gf(kx2)-0.5[gf(kx1)+gf(k)]
[0043]
skx2＝gf(k)-0.5[gf(kx2)+gf(kx3)]
[0044]
skx3＝gf(kx3)-0.5[gf(k)+gf(kx4)]
[0045]
sky1＝gf(ky2)-0.5[gf(ky1)+gf(k)]
[0046]
sky2＝gf(k)-0.5[gf(ky2)+gf(ky3)]
[0047]
sky3＝gf(ky3)-0.5[gf(k)+gf(ky4)]
[0048]
式中，kx1、kx2、k、kx3、kx4为x方向上间隔距离为滤波窗口半宽度wd的测点位置，ky1、ky2、k、ky3、ky4为y方向上间隔距离为滤波窗口半宽度wd的测点位置，skx1、skx2、skx3、sky1、sky2、sky3为图1(右)箭头所指的线段长度(对应的纵坐标长度)。
unit，cpu)，还可以是其他通用处理器、数字信号处理器(digital signal processor、dsp)、专用集成电路(application specific integrated circuit，asic)、现成可编程门阵列(field-programmable gatearray，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等，其是终端的计算核心以及控制核心，其适于实现一条或一条以上指令，具体适于加载并执行一条或一条以上指令从而实现相应方法流程或相应功能；本发明实施例所述的处理器可以用于上述提取局部重力异常的方法的操作，包括以下步骤：
[0062]
s1：基于矿产勘探或地质构造研究中获得的航空空间重力异常，经布格校正和网格化处理得到航空布格重力异常网格数据。
[0063]
s2：将原始重力异常网格数据作为第一次迭代计算的初始区域异常gf，估算给定滤波窗口宽度2wd内异常弯曲方向的特征值；
[0064]
s3：判断异常弯曲方向是否受到叠加的局部异常的影响而发生交替变化，如果是则更新区域异常值，如果不是则区域异常值保留原场值不变。
[0065]
s4：在网格上依次从左至右、从上到下移动窗口，重复s2～s4步骤，直至整个网格计算完毕为止，计算得到的gf值为格网区域异常的第一次迭代计算值；
[0066]
s5：逐步减小滤波半窗口，每次减小滤波半窗口，就是在前一次区域异常迭代计算值gf的基础上，重复s2～s5步骤；以此逐次迭代逼近，直到滤波半窗口宽度减小为网格间距，最终迭代得到gf值就是基于网格的曲率滤波法计算获得的滤波窗口宽度为2wd的位场区域异常值；
[0067]
s6：原始重力异常格网与最终迭代计算得到的区域异常格网做差，便得到局部异常。
[0068]
本发明再一个实施例中，本发明还提供了一种存储介质，具体为计算机可读存储介质(memory)，所述计算机可读存储介质是终端设备中的记忆设备，用于存放程序和数据。可以理解的是，此处的计算机可读存储介质既可以包括终端设备中的内置存储介质，当然也可以包括终端设备所支持的扩展存储介质。计算机可读存储介质提供存储空间，该存储空间存储了终端的操作系统。并且，在该存储空间中还存放了适于被处理器加载并执行的一条或一条以上的指令，这些指令可以是一个或一个以上的计算机程序(包括程序代码)。需要说明的是，此处的计算机可读存储介质可以是高速ram存储器，也可以是非不稳定的存储器(non-volatile memory)，例如至少一个磁盘存储器。
[0069]
可由处理器加载并执行计算机可读存储介质中存放的一条或一条以上指令，以实现上述实施例中有关提取局部重力异常的方法的相应步骤；计算机可读存储介质中的一条或一条以上指令由处理器加载并执行如下步骤：
[0070]
s1：基于矿产勘探或地质构造研究中获得的航空空间重力异常，经布格校正和网格化处理得到航空布格重力异常网格数据。
[0071]
s2：将原始重力异常网格数据作为第一次迭代计算的初始区域异常gf，估算给定滤波窗口宽度2wd内异常弯曲方向的特征值；
[0072]
s3：判断异常弯曲方向是否受到叠加的局部异常的影响而发生交替变化，如果是则更新区域异常值，如果不是则区域异常值保留原场值不变。
[0073]
s4：在网格上依次从左至右、从上到下移动窗口，重复s2～s4步骤，直至整个网格
计算完毕为止，计算得到的gf值为格网区域异常的第一次迭代计算值；
[0074]
s5：逐步减小滤波半窗口，每次减小滤波半窗口，就是在前一次区域异常迭代计算值gf的基础上，重复s2～s5步骤；以此逐次迭代逼近，直到滤波半窗口宽度减小为网格间距，最终迭代得到gf值就是基于网格的曲率滤波法计算获得的滤波窗口宽度为2wd的位场区域异常值；
[0075]
s6：原始重力异常格网与最终迭代计算得到的区域异常格网做差，便得到局部异常。
[0076]
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。
[0077]
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0078]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0079]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0080]
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的实施方法，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。