CN111007571B - 一种基于三维构造张量的航磁数据地质体边界识别方法 - Google Patents

Abstract

本发明涉及一种基于三维构造张量的航磁数据地质体边界识别方法，包括：获取实测航磁异常数据并预处理；对预处理后的数据进行空间域计算；对预处理后的数据进行频率域计算；利用空间域计算的数据和频率域计算的数据得到水平导数、多阶垂向导数以及构造张量矩阵得到边界识别滤波器NSED1和边界识别滤波器NSED2；本发明突破了以往边界识别滤波器采用单梯度和同阶梯度进行滤波器定义的思想，定义的多梯度分量参与的边界识别滤波器在深浅叠加异常、临近叠加异常、正负叠加异常等复杂地质情况下均可获得更清晰、可靠的地质体水平位置圈定结果，避免叠加异常产生的干扰信息，提高了弱异常、深部异常对应地质体的边界位置增强和提取能力。

Description

一种基于三维构造张量的航磁数据地质体边界识别方法

技术领域

本发明涉及地球科学技术领域，具体为一种基于三维构造张量的航磁数据地质体边界识别方法。

背景技术

随着基于有人机、无人机、无人飞艇等移动平台航磁探测系统的不断发展，航磁探测获取的大面积、高精度数据的快速解释是现在研究的关键问题。而边界检测方法可圈定出地质体的水平位置，进而指导下一步的解释工作，被广泛地应用到航磁数据的处理中。目前，有很多方法用于地质构造边界的检测和拾取中。航磁异常的总水平导数最大值和垂向导数零值线位置与地质体的边界位置对应。为了进一步提高边界识别的分辨率，后来发展了多种基于不同阶导数比值的边界识别滤波器。但常规方法在深浅叠加异常、临近叠加异常等复杂地质情况下的边界识别方面存在精度低、分辨率差的问题。

发明内容

本发明所要解决的技术问题在于提供一种基于三维构造张量的航磁数据地质体边界识别方法，解决在深浅叠加异常、临近叠加异常等复杂地质情况下的边界识别方面存在精度低、分辨率差的问题。

本发明是这样实现的，一种基于三维构造张量的航磁数据地质体边界识别方法，该方法包括：

获取实测航磁异常数据；

对数据进行预处理；

对预处理后的数据进行空间域计算；

对预处理后的数据进行频率域计算；

利用空间域计算的数据和频率域计算的数据得到水平导数、多阶垂向导数以及构造张量矩阵；

利用张量矩阵得到边界识别滤波器NSED1和边界识别滤波器NSED2。

进一步地，对预处理后的数据进行空间域计算包括：

采用空间域求导方式得到异常的水平导数

表达式为：

其中，f为实测航磁异常数据，x、y为空间坐标的两个方向,x₁与y₁以及x₂与y₂为不同时间点航次异常数据对应的空间坐标方向。

进一步地，对预处理后的数据进行频率域计算包括：

通过频率域求导方式得到垂向导数

f为实测航磁异常数据；x、y、z为空间坐标的三个方向。频率域求导先将原始异常f做傅氏变换，用它的谱乘以波谱算子r(径向圆波数)，得到的结果再做反傅氏变换便得到垂向导数

进一步地，利用空间域和频率域计算方法得到水平导数、垂向导数以及构造张量矩阵：

根据实测或计算获得的梯度值

定义构造张量矩阵T，表达式为：

进一步地，采用边界识别滤波器NSED₁进行航磁数据的边界识别，经验参数k选择为0.0001，NSED₁表达式为：

进一步地，采用边界识别滤波器NSED₂进行航磁数据的边界识别，进一步提升边界识别结果的分辨率，通过试验确定参数k为0.0001，NSED₂表达式为：

本发明与现有技术相比，有益效果在于：

本发明突破了以往边界识别滤波器采用单梯度和同阶梯度进行滤波器定义的思想，引入用于图像处理的构造张量矩阵，并定义了基于航磁梯度数据的构造张量矩阵，进而建立了基于构造张量矩阵的边界识别滤波器，在深浅叠加异常、临近叠加异常、正负叠加异常等复杂地质情况下均可获得更清晰、可靠的地质体水平位置圈定结果，避免了产生叠加异常产生的干扰信息。多梯度分量参与边界识别滤波器的定义大大地丰富了数据中包含的地质体信息，提高了弱异常、深部异常对应地质体的边界位置增强和提取能力，具有更高的分辨率和精度。

附图说明

图1为一种基于三维构造张量的航磁数据地质目标体的边界识别方法流程图；

图2a为朱日和地区航磁异常图；

图2b为NSED1边界识别结果；

图2c为NSED2边界识别结果；

图2d为Theta图边界识别结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

请参阅图1，本发明提供一种基于三维构造张量的航磁数据地质目标体的边界识别方法，包括以下步骤：

获取实测航磁异常数据；

对数据进行预处理；

对预处理后的数据进行空间域计算；

对预处理后的数据进行频率域计算；

利用空间域计算的数据和频率域计算的数据得到水平导数、垂向导数以及构造张量矩阵；

如何利用上述的步骤得到NSED1和NSED2；

利用NSED1和NSED2得到边界位置。

其中，针对实测航磁异常数据，采用数据预处理手段进行有用异常信号的提取，进而计算对应的梯度分量，为后续地质目标体的边界识别提供数据源；

对预处理后的数据进行空间域计算包括：

采用空间域求导方式得到异常的水平导数

其中，f为实测航磁异常数据，x、y为空间坐标的两个方向,x₁与y₁以及x₂与y₂为不同时间点航次异常数据对应的空间坐标方向。

进一步地，对预处理后的数据进行频率域计算包括：

通过频率域求导方式得到垂向导数

f为实测航磁异常数据；

利用空间域计算的数据和频率域计算的数据定义航磁构造张量矩阵，表达式为：

综合利用航磁构造张量矩阵的全部分量定义边界识别滤波器：

式中，k取值0-0.1，由经验获得。

对于实测航磁异常数据中未获取的梯度分量，采用空间域求导和频率域计算的两种策略进行梯度分量的计算，进而定义出航磁数据对应的构造张量矩阵，采用两种形式的边界识别滤波器进行场源水平位置的圈定，k经验参数在选择时进行不同取值的尝试，通过解释结果的优劣确定最优的参数。

实施例：

以朱日和地区的实测航磁数据为例，通过基于构造张量的边界识别滤波器识别该地区的断裂分布。

具体步骤如下：

A、根据朱日和地区的实测航磁异常f见图2(a)，采用数据预处理技术进行数据处理，包括磁补偿、日变校正、正常场校正、高度校正、滞后校正、调平，并采用空间域求导方式得到异常的多阶水平导数(梯度值)

通过频率域求导方式得到多阶垂向导数(梯度分量)

B、根据实测或计算获得的梯度值

定义构造张量矩阵，表达式为：

C、采用边界识别滤波器NSED1进行航磁数据的边界识别，经验参数k选择为0.0001，图2(b)为该方法的解释结果，NSED1表达式为：

D、采用边界识别滤波器NSED2进行航磁数据的边界识别，进一步提升边界识别结果的分辨率，通过试验确定参数k为0.0001，图2(c)为该方法的解释结果，NSED2表达式为：：

为了对比本方法的分辨能力，选择现有的Theta图(Wijns，2005)进行对比，其表达式为：

图2(d)为Theta图边界识别结果，相比基于构造张量的边界识别方法具有更差的分辨率，圈定的构造边界位置模糊，而本发明提出方法可更清晰的显示局部弱异常对应断裂的边界位置，提升了航磁弱异常的解释能力。

在实际应用中，选择多方法进行航磁异常边界位置的拾取，将其结果进行相互印证，可更有助于后期的构造划分和断裂识别，获得可靠的解释结果。

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

Claims (1)

1.一种基于三维构造张量的航磁数据地质体边界识别方法，其特征在于，该方法包括：

获取实测航磁异常数据；

对数据进行预处理；

对预处理后的数据进行空间域计算；

对预处理后的数据进行频率域计算；

利用空间域计算的数据和频率域计算的数据得到水平导数、多阶垂向导数以及构造张量矩阵；

利用张量矩阵得到两种具有不同分辨率的边界识别滤波器NSED₁和边界识别滤波器NSED₂；

对预处理后的数据进行空间域计算包括：

采用空间域求导方式得到异常的水平导数

表达式为：

其中，f为实测航磁异常数据，x、y为空间坐标的两个方向,x₁与y₁以及x₂与y₂为不同时间点航磁异常数据对应的空间坐标方向；

对预处理后的数据进行频率域计算包括：

通过频率域求导方式得到垂向导数

f为实测航磁异常数据；x、y、z为空间坐标的三个方向，频率域求导先将实测航磁异常数据f做傅氏变换，用它的谱乘以波谱算子r，得到的结果再做反傅氏变换便得到垂向导数

利用空间域和频率域计算方法得到水平导数、垂向导数以及构造张量矩阵：

根据实测或计算获得的梯度值

定义构造张量矩阵T，表达式为：

采用边界识别滤波器NSED₁进行航磁数据的边界识别，经验参数k选择为0.0001，NSED₁表达式为：

采用边界识别滤波器NSED₂进行航磁数据的边界识别，进一步提升边界识别结果的分辨率，通过试验确定参数k为0.0001，NSED₂表达式为：

Images