一种基于位场数据的鲁棒的高分辨率地质体边界识别方法

技术领域

本发明涉及一种基于位场数据的鲁棒的高分辨率地质体边界识别方法，属于地球物理勘探技术领域。

背景技术

位场边缘检测技术是一种历史悠久且应用广泛的地质边界识别手段。该技术能够突出原始位场数据中不易察觉的信息，在各种尺度的构造解释，如矿体探测、断层识别、不同地质体的圈定和区域构造格架的划分等中均发挥重要作用。基于过去的研究，位场边缘检测主要可分为三大类：（1）基于导数类，（2）基于相位或比值类和（3）基于统计/滑动窗口类。其中，水平导数和垂向导数（VDR法）是导数类方法的基础，例如，较早提出的总水平导数（THDR法）就是由不同方向的水平导数构成。倾斜角（TA法）是典型的基于相位或比值类的方法。最近还提出了ILTHG法和基于滑动窗口的NHF法等。

研究发现，同时利用顺序组合和比值的混合型边缘检测滤波器比传统的单一的基于导数、比值或滑动窗口的边缘检测器有更高的水平或垂向分辨能力，且更不容易产生假象。但是，伴有高阶导数的混合型滤波器，如ILTHG法，往往更容易受噪声干扰，在这种情况下，效果不如基于统计时窗的NHF法滤波器。总体来看，尽管位场边缘检测技术在不断发展，目前大多数滤波器都存在以下局限性：（1）对不同深度、不同尺寸和不同强度的场源的响应很难达到平衡；（2）低阶的边缘检测器分辨率有限，而对噪声的抵抗能力随着导数阶数的升高而迅速降低，即面临分辨率和抗噪性之间的折中。

发明内容

针对现有技术存在的上述缺陷，本发明提出了一种新的基于位场数据的鲁棒的高分辨率地质体边界识别方法，能够提高抗干扰能力，产生高分辨率的边缘检测结果，且不易产生虚假边界。

本发明是采用以下的技术方案实现的：本发明所述一种基于位场数据的鲁棒的高分辨率地质体边界识别方法，包括如下步骤：

步骤一：判断位场数据的信噪比是否满足要求，若信噪比太低，则对其做向上延拓或其他去噪处理，将去噪后的位场数据作为边缘检测的输入数据

步骤二：计算输入数据

S1：计算公式如下：

(1)

该步骤一般在频率-波数域中进行，然后经过傅里叶反变换转为空间域，从而得到

步骤三：对

S2：在任一3×3大小的窗口内，计算

(2)

其中，

S3：计算Harris滤波的响应函数，具体公式为

(3)

其中，Det(M)表示M的行列式，

步骤四：对

S4：二维高斯平滑具体表示为

(4)

其中，Gaussian()表示二维高斯分布函数，从而得到

步骤五：对

S5：分别计算

S6：将

(5)

其中，exp表示以自然常数e为底的指数函数，

步骤六：计算公式(5)从而得到

本发明的有益效果是：采用本发明所述的一种基于位场数据的鲁棒的高分辨率地质体边界识别方法，通过融合多种边缘检测滤波器的优势，从而能够得到高分辨率的地质体边界识别结果，且不容易产生虚假边界，在一定噪声强度下仍然具有鲁棒性，解决了现有方法不能兼顾高分辨率和强抗噪性的问题，抗噪能力和边界分辨率优于传统方法，边界识别结果可靠。

附图说明

图1为本发明的流程图；

图2为3D理论模型图；

图3为模型产生的重力异常图；

图4(a)为无噪情况下VDR法地质边界识别结果图；

图4(b)为无噪情况下THDR法地质边界识别结果图；

图4(c)为无噪情况下TA法地质边界识别结果图；

图4(d)为无噪情况下NHF法地质边界识别结果图；

图4(e)为无噪情况下ILTHG法地质边界识别结果图；

图4(f)为无噪情况下本发明方法地质边界识别结果图；

图5(a)为待添加的随机噪声图；

图5(b)为添加随机噪声后的重力异常图；

图6(a)为添加随机噪声后VDR法地质边界识别结果图；

图6(b)为添加随机噪声后THDR法地质边界识别结果图；

图6(c)为添加随机噪声后TA法地质边界识别结果图；

图6(d)为添加随机噪声后NHF法地质边界识别结果图；

图6(e)为添加随机噪声后ILTHG法地质边界识别结果图；

图6(f)为添加随机噪声后本发明方法地质边界识别结果图；

具体实施方式

为了使本发明目的、技术方案更加清楚明白，下面结合附图，对本发明作进一步详细说明。本发明的流程图，如图1所示，包括如下步骤：

步骤一：判断位场数据的信噪比是否满足要求，若信噪比太低，则对其做向上延拓或其他去噪处理，将去噪后的位场数据作为边缘检测的输入数据

步骤二：计算输入数据

S1：计算公式如下：

(1)

该步骤一般在频率-波数域中进行，然后经过傅里叶反变换转为空间域，从而得到

步骤三：对

S2：在任一3×3大小的窗口内，计算

(2)

其中，

S3：计算Harris滤波的响应函数，具体公式为

(3)

其中，Det(M)表示M的行列式，

步骤四：对

S4：二维高斯平滑具体表示为

(4)

其中，Gaussian()表示二维高斯分布函数，从而得到

步骤五：对

S5：分别计算

S6：将

(5)

其中，exp表示以自然常数e为底的指数函数，

步骤六：计算公式(5)从而得到

下面结合具体实施方式，对于本发明的模型测试进行解释和说明。

实施例一：

为了说明本方法的实现思路及实现过程并证明方法的有效性，建立一个三维模型进行测试，并和传统边缘检测方法的结果进行比较。

S1:建立如图2所示的理论三维模型。该模型由三个正密度和三个负密度的棱柱体组成，棱柱体①长为15 km，宽为50 km，高为5 km，密度差为-0.3 g/cm

S2:在此模型范围内，沿水平地表（z=0 km）布设重力观测点，测点间距为2 km。

S3:计算此模型产生的重力异常

S4:计算输入数据

S5:对

S6:对

S7:提取

为了说明本发明方法的效果，将本发明地质体边界识别结果与现有方法（包括VDR法、THDR法、TA法、NHF法和ILTHG法）结果进行对比。图4(a~e)分别为VDR法、THDR法、TA法、NHF法和ILTHG法的地质边界识别结果，图4(f)是本发明的地质边界识别结果。从图中可以明显看出，VDR法在零值点附近与真实的块体边界位置对应较好，然而，并不是所有的零值点都指示真实的边界位置，因此，VDR具有模糊性和虚假边界。THDR法的极大值指示了边界的位置，但具有较低的分辨率，对深源小目标体（棱柱体③）的敏感性较低。TA法的表现与VDR法大体相当，继承了VDR法的缺点，即存在明显的过剩的跨零点，导致了虚假边界的产生。NHF法比THDR法有更强的分辨率，但同样不能检测到深部的小块体。相比之下，ILTHG法和本发明方法更能够平衡不同大小、不同深度的场源的振幅强度，六个块体边界都被检测出来。相比与ILTHG法，本发明方法边缘线理更细，块体边界表现更尖锐，表明具有更高的分辨率。

实施例二：

为了进一步说明本方法的抗噪性和鲁棒性，对实施例一中模型产生的重力异常添加如图5(a)所示的随机噪声进行测试，并和现有方法（包括VDR法、THDR法、TA法、NHF法和ILTHG法）结果进行比较。添加随机噪声后的重力异常如图5所示。将加随机噪声后的重力异常做4km的向上延拓处理，作为边缘检测的输入数据

可以看到，VDR法和THDR法的结果相对模糊，且在棱柱体③和棱柱体⑤相接的位置存在边界的错移。TA法的结果与上两种方法的结果类似，但看起来受噪声的影响更大。NHF法受噪声的影响小，但分辨能力依然有限。ILTHG法的边界分辨力相对较高，但存在虚假的琐碎斑点，较大的噪声水平造成了ILTHG法的不稳定。相比于ILTHG法，本发明方法极大地降低了检测出离散虚假点作为边界的可能性。相较于NHF法，本发明方法显著地提高了边缘检测的水平分辨率，对边界的定位更加明确和可靠。可见，本发明方法优于现有的边缘检测方法。

以上所述仅为本发明的较佳实施例而己，并不以本发明为限制，凡在本发明的精神和原则之内所作的均等修改、等同替换和改进等，均应包含在本发明的专利涵盖范围内。