CN104658037A - 一种位场构造格架自动提取方法 - Google Patents

Abstract

本发明提供一种位场构造格架自动提取方法。该方法包括以下步骤：对来自待研究区域的重力位场数据或磁法位场数据进行预处理；对经预处理的重力位场数据或磁法位场数据进行多个尺度下多方向的边缘检测，分别得到各尺度的边缘；采用形态学骨骼算法将计算得到的各尺度的边缘细化为单像素宽度，并且每个点具有深度和强度属性，得到重力或磁综合构造格架图和综合构造强度格架图。根据本发明的方法，可以实现对控制矿床形成的地质构造的识别和定性解释，根据研究区先验知识确定潜在的矿床类型和控制矿床形成的构造属性，对不同类型构造格架进行筛选，从而实现金属矿床的靶区定位。

Description

一种位场构造格架自动提取方法

技术领域

本发明涉及一种重力位场和磁位场构造格架自动提取方法，是一项基于重力和磁异常数据进行地质构造探测的技术。更具体地，本发明涉及小波分析、图像处理、地球物理、地质学、矿产勘查等领域，根据本发明的方法能直接应用于矿产勘探和相关地质调查领域等。

背景技术

随着航空物探技术的发展，重力和磁法(下文中也简称为重磁)测量方法具有经济、快速、可以覆盖很多难以到达的景观地区等优点，具有广阔的应用前景。目前在矿产勘查领域，重力和磁法勘探方法通常被应用于具有强磁性或高密度的矿化体的直接探测，以及具有强异常的地质构造解释和反演。大部分有色金属和稀有金属矿床不能直接形成明显重力和磁异常，但是控制矿床形成的地质构造可以产生异常差异。所以，基于重力位场数据和磁位场数据进行控制矿床形成的地质构造探测，在金属矿床靶区定位中具有重要意义。

目前重磁位场构造信息自动识别和提取方法主要包括解析信号法、欧拉反褶积法、位场多尺度边缘检测方法等，其共同存在的问题是对方向性的信息不敏感，不能得到完整、准确的异常边界位置。在本发明申请人的发明名称为“位场多方向多尺度边缘检测方法”，申请号为200810006676.6的专利申请中，公开了一种位场多方向多尺度边缘检测方法，通过方向小波变换增强了方向性信息，获得了不同方向的异常源边界信息，实现了构造格架的自动提取，克服了解析信号法、欧拉反褶积法、位场多尺度边缘检测方法等对方向性信息不敏感的问题，是一项基于位场数据进行浅部地壳三维结构快速反演的技术方案。但该方法仍存在以下问题：(1)计算得到的边缘不是单像素点宽度，边缘及不同方向边缘交点所对应的实际地理范围较大，导致分析结果精确度不高；(2)需要基于不同尺度边缘进行人工矢量化，取边缘中心连线作为异常边界，从而生成浅部地壳三维结构图，导致工作效率低；(3)该浅部地壳三维结构图只能显示不同深度的构造信息，不能反映构造带及两侧岩性变化，不能指示构造的变形和活动强度；(4)对尺度没有明确定义，也就是没有对尺度这一特征赋予具体的地球物理属性。

因此，需要提供一种能够基于已采集的位场数据高精度地提取位场构造格架的方法。

发明内容

本发明的目的在于针对现有技术的上述不足，提供一种重磁位场构造格架自动提取方法，能够快速得到控制矿床形成的地质构造信息，从而实现金属矿床靶区定位。

为实现这样的目的，本发明对测量得到的重力位场数据和磁位场数据采用多方向多尺度边缘检测方法提取重力构造格架和磁构造格架，采用形态学骨骼算法将得到的每一尺度构造格架细化处理成单像素宽度，将不同尺度的构造格架用不同颜色渐变叠加显示突出不同深度构造信息，生成综合构造格架图；将基于不同尺度构造格架各边缘点对应的梯度的模提取的构造强度值进行不同渐变颜色的叠加显示，突出密度变化和磁性变化强度信息，生成综合构造强度格架图。从而可以获得该研究区范围不同深度的重磁异常信息、表征不同深度地质构造的构造格架分布信息、不同深度构造格架的密度变化和磁性变化强度信息，实现对控制矿床形成的地质构造的识别和定性解释，根据研究分析区先验知识确定潜在的矿床类型和控制矿床形成的构造属性，对不同类型构造格架进行筛选，从而实现金属矿床的靶区定位。

本发明的位场构造格架自动提取方法，包括如下的具体步骤：

1)对从研究区域重磁测量获得的重力数据或磁法数据进行预处理

该预处理包括对磁法数据进行化极计算得到化极磁异常或进行伪重力计算得到伪重力异常；对重力数据进行预处理得到布格重力异常数据。

2)对经预处理的重力位场数据或磁法位场数据进行多个尺度下位场多方向边缘检测，包括将位场数据上延多个预定高度后得到多个尺度的重力位场数据或磁法位场数据，并分别对所得到的每一尺度下的位场数据进行多方向边缘检测，得到各尺度下位场边缘。

针对上延多个预定高度后的每一尺度，选择不同的方向α进行边缘检测计算，可突出不同方向的边缘信息。为了能够达到完整覆盖，方向α取值为kπ/(2^n-1)，其中k＝0,1,2,…,(2ⁿ-1)，n为大于或等于2的整数。针对每一方向，经预处理的重力数据或磁法数据的二维方向小波变换的模极大值点沿梯度的垂直方向连接得到的曲线构成边缘。获取同一尺度以不同方向计算得到的边缘，对所获取得到的边缘求并集作为该尺度的边缘，从而可实现各尺度下位场多方向的边缘检测。

3)采用形态学骨骼算法将计算得到的各尺度的边缘分别细化处理为单像素宽度，得到各尺度下的构造格架。

4)将计算得到的各尺度构造格架叠置生成综合构造格架图

对得到的各尺度的构造格架进行叠置生成反映不同深度信息的综合构造格架图，图上不同尺度边缘的横向偏移反应了构造格架的产状信息。

将位场上延不同高度后提取的边缘对应于不同深度的构造，深度为向上延拓后高度的一半(见作者为Jacobsen,B.H.，题目为A case for upwardcontinuation as a standard separation filter for potential-field maps，期刊名Geophysics，刊号v.52no.8，时间1987)，可以得到表征不同深度的构造格架图。

5)将计算得到的边缘上的每一边缘点处的梯度的模作为构造格架中该边缘点处强度值。赋有强度值的构造格架上边缘点强度值的高低反映了构造带及两侧岩性变化大小，反映了构造的变形和活动强度。建立反映不同深度信息的不同尺度的强度格架图，以及不同尺度叠加成综合构造强度格架图，揭示了区域地质构造格局。

具体地，本发明提供一种位场构造格架自动提取方法，包括以下步骤：

对来自待研究区域的位场数据进行预处理；

将经预处理的位场数据向上延拓多个预定高度得到多个相应尺度的位场数据；

分别针对每一尺度的位场数据进行多方向边缘检测，得到多个相应尺度的位场边缘；

采用形态学骨骼算法将计算得到的各尺度的位场边缘分别细化为单像素宽度，得到多个相应尺度的构造格架图。

优选地，该方法进一步包括将计算得到的所述多个相应尺度的构造格架图叠置生成综合构造格架图。

优选地，该方法进一步包括将各尺度构造格架图上的每一边缘点处的梯度的模作为该尺度构造格架图中该边缘点处的强度值，得到多个相应尺度的构造强度格架图。

优选地，该方法进一步包括将所述多个相应尺度的构造强度格架图叠置生成综合构造强度格架图。

优选地，所述位场数据为重力位场数据或磁法位场数据，所述预处理进一步包括对磁法数据进行化极计算得到化极磁异常或进行伪重力计算得到伪重力异常；或对重力数据进行预处理得到布格重力异常。

优选地，针对每一尺度的位场数据进行多方向边缘检测，包括以下步骤：

设尺度参数s＝z/z₀，且z>z₀，z₀代表测量高度，z代表向上延拓后的高度，定义高度为零的位置(x,y)处的重力异常或磁异常为f₀(x,y)，

尺度s的平滑函数定义为：

${\mathrm{\θ}}_s(x,y)=s^{-2}\mathrm{\θ}(x/s+y/s)={\mathrm{\γ}}_{s{z}_0}(x,y),$

其中 ${\mathrm{\γ}}_{s{z}_0}(x,y)=\frac{1}{2\mathrm{\π}}k(x,y,s{z}_0),$ k(x,y,z)为格林函数，

在方向α的小波函数定义为：

${\mathrm{\ψ}}^{\mathrm{\α}}=(\cos \mathrm{\α}\frac{\∂}{\∂x}+\sin \mathrm{\α}\frac{\∂}{\∂y})\mathrm{\θ}(x,y)=D^{\mathrm{\α}}\mathrm{\θ}(x,y),$

其中，D表示一阶导数；

在尺度s和位置(x,y)的情况下，重力异常或磁异常f₀(x,y)在方向α的小波变换定义为：

$\begin{array}{c}{W}^{\mathrm{\α}}\left[f_0\right](x,y,s)=[f_0*{\mathrm{\ψ}}_s^{\mathrm{\α}}](x,y)\\ =f_0*\left(s{D}^{\mathrm{\α}}{\mathrm{\θ}}_s\right)\\ =s{D}^{\mathrm{\α}}[f_0*{\mathrm{\θ}}_s]\\ =s{D}^{\mathrm{\α}}[f_0*{\mathrm{\γ}}_{s{z}_0}]\end{array}$

其中，*表示卷积运算，

由位场向上延拓公式已知 $f_z(x,y)=f_0(x,y)*{\mathrm{\γ}}_{s{z}_0}(x,y),$

f_z(x,y)为f₀(x,y)从高度零上延到高度z＝sz₀处的重力异常或磁异常，通过将在测量高度z₀测量得到的重力异常或磁异常向上延拓高度z-z₀获得，

因此，W^α[f₀](x,y,s)＝sD^αf_z(x,y)

           ＝(z/z₀)D^αf_z(x,y)；

进一步，在尺度s和位置(x,y)，重力异常或磁异常f₀(x,y)在方向的小波变换定义为：

$W^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}\left[f_0\right](x,y,s)=(z/z_0)D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)$

则，f₀(x,y)的二维方向小波变换用梯度表示为：

$W\left[f_0\right](x,y,s,\mathrm{\α})=\left(\begin{array}{c}{W}^{\mathrm{\α}}\left[f_0\right](x,y,s)\\ W^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}\left[f_0\right](x,y,s)\end{array}\right)=(z/z_0)\left(\begin{array}{c}{D}^{\mathrm{\α}}{f}_z(x,y)\\ D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)\end{array}\right)$

$W\left[f_0\right](x,y,s,\mathrm{\α})=(z/z_0)\▿f_z(x,y,\mathrm{\α})$

其中为二维梯度。

对于位置(x,y)、尺度s和方向α，f₀(x,y)的二维方向小波变换W[f₀](x,y,s,α)与f_z(x,y)的梯度成正比，f₀(x,y)的二维方向小波变换W[f₀](x,y,s,α)可以用f_z(x,y)的梯度来表征。

对于高度z,定义梯度的模为：

$M\left[f_z\right](x,y,\mathrm{\α})=\sqrt{{\left|D^{\mathrm{\α}}{f}_z(x,y)\right|}^2+{\left|D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)\right|}^2},$

该梯度相应的沿水平方向的辐角为：

${\mathrm{Af}}_z(x,y,\mathrm{\α})=\arg \mathrm{ument}(D^{\mathrm{\α}}{f}_z(x,y)+D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)),$

边缘点为模M[f_z](x,y,α)沿辐角方向Af_z(x,y,α)有局部极大值的点，

针对每一方向α，梯度的模的局部极大值点沿梯度的垂直方向连接得到的曲线构成边缘，

针对同一高度，以多个不同的方向α计算边缘，对计算得到的各边缘求并集得到该高度的边缘。

优选地，针对同一高度以多个不同的方向α计算边缘的步骤进一步包括，各方向α取值为kπ/(2^n-1)，其中k＝0,1,2,…,(2ⁿ-1)，n为大于或等于2的整数，以完整覆盖二维平面。

优选地，位场向上延拓多个预定高度后提取的边缘对应于不同深度的构造，对得到的各尺度的构造格架图进行叠加得到反映不同深度信息的综合构造格架图。

优选地，以渐变颜色表征高度的高低，来形成所述综合构造格架图。

优选地，以渐变颜色表征强度值的大小，来形成所述综合构造强度格架图。

根据本发明的方法，可以获得该研究区范围不同深度的重磁异常信息、表征不同深度地质构造的构造格架分布信息、不同深度构造格架的磁性和密度变化强度信息，实现对控制矿床形成的地质构造的识别和定性解释，根据研究区先验知识确定潜在的矿床类型和控制矿床形成的构造属性，对不同类型构造格架进行筛选，从而实现金属矿床的靶区定位。

从参照附图对实例性实施例的以下描述，本发明的进一步特征将变得明显。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为位场构造格架自动提取方法流程图；

图2根据本发明实例的单像素宽度构造格架图；

图3为根据本发明实例的综合构造格架图；

图4为根据本发明实例的构造强度格架图；

图5根据本发明实例的综合构造强度格架图；

具体实施方式

为了更好的理解本发明的技术方案，以下结合附图和具体实施例对本发明的实施方式作进一步的描述，但不作为对本发明的限定。

图1为本发明的位场构造格架自动提取方法的流程图。如图1所示，该方法包括如下步骤：

步骤101，对测量得到的重力位场数据或磁位场数据预处理。

对磁法数据进行化极计算得到化极磁异常或进行伪重力计算得到伪重力异常。

对重力数据进行预处理得到布格重力异常。

步骤102，采用多个尺度下位场多方向边缘检测方法对经预处理的重力位场数据或磁位场数据进行处理得到多个尺度的边缘。

位场多方向多尺度边缘检测方法包括将经预处理的重力位场数据或磁法位场数据向上延拓多个预定高度得到多个相应尺度的位场数据和分别针对每一尺度的位场进行多方向边缘检测得到各尺度位场边缘的步骤。

分别针对每一尺度的位场进行多方向边缘检测的计算方法，包括以下步骤：

定义地球表面点(x,y)处高程为零时的重力异常或磁异常为f₀(x,y)。

设尺度参数s＝z/z₀，且z>z₀，z₀代表测量高度，z代表向上延后的高度。

尺度s的平滑函数定义为：

${\mathrm{\θ}}_s(x,y)=s^{-2}\mathrm{\θ}(x/s+y/s)={\mathrm{\γ}}_{s{z}_0}(x,y),$

其中 ${\mathrm{\γ}}_{s{z}_0}(x,y)=\frac{1}{2\mathrm{\π}}k(x,y,s{z}_0),$ k(x,y,z)为格林函数。

在方向α的小波函数可以定义为：

${\mathrm{\ψ}}^{\mathrm{\α}}=(\cos \mathrm{\α}\frac{\∂}{\∂x}+\sin \mathrm{\α}\frac{\∂}{\∂y})\mathrm{\θ}(x,y)=D^{\mathrm{\α}}\mathrm{\θ}(x,y),$

其中，D表示一阶导数；

在尺度s和位置(x,y)的情况下，重力异常或磁异常f₀(x,y)在方向α的小波变换定义为：

$\begin{array}{c}{W}^{\mathrm{\α}}\left[f_0\right](x,y,s)=[f_0*{\mathrm{\ψ}}_s^{\mathrm{\α}}](x,y)\\ =f_0*\left(s{D}^{\mathrm{\α}}{\mathrm{\θ}}_s\right)\\ =s{D}^{\mathrm{\α}}[f_0*{\mathrm{\θ}}_s]\\ =s{D}^{\mathrm{\α}}[f_0*{\mathrm{\γ}}_{s{z}_0}]\end{array}$

其中，*表示卷积运算，

由位场向上延拓公式已知：

$f_z(x,y)=f_0(x,y)*{\mathrm{\γ}}_{s{z}_0}(x,y),$

f_z(x,y)为f₀(x,y)从高度零上延到高度z＝sz₀处的重力异常或磁异常，通过将在测量高度z₀测量得到的重力异常或磁异常f_z0(x,y)向上延拓高度z-z₀获得，

因此，W^α[f₀](x,y,s)＝sD^αf_z(x,y)

             ＝(z/z₀)D^αf_z(x,y)；

同样，在尺度s和位置(x,y)，重力异常或磁异常f₀(x,y)在方向的小波变换定义为：

$W^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}\left[f_0\right](x,y,s)=(z/z_0)D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)$

f₀(x,y)的二维方向小波变换可以用梯度来表示：

$W\left[f_0\right](x,y,s,\mathrm{\α})=\left(\begin{array}{c}{W}^{\mathrm{\α}}\left[f_0\right](x,y,s)\\ W^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}\left[f_0\right](x,y,s)\end{array}\right)=(z/z_0)\left(\begin{array}{c}{D}^{\mathrm{\α}}{f}_z(x,y)\\ D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)\end{array}\right)$

$W\left[f_0\right](x,y,s,\mathrm{\α})=(z/z_0)\▿f_z(x,y,\mathrm{\α})$

其中为二维梯度，

上式建立了任意高度z(z>z₀)的重力异常或磁异常f_z(x,y)的水平梯度与高度为零的重力异常或磁异常f₀(x,y)的二维方向小波变换的联系。

对于位置(x,y)、尺度s和方向α，f₀(x,y)的二维方向小波变换W[f₀](x,y,s,α)与f_z(x,y)的梯度成正比，f₀(x,y)的二维方向小波变换W[f₀](x,y,s,α)可以用f_z(x,y)的梯度来表征。

对于向上延拓后的高度z,定义梯度的模为：

$M\left[f_z\right](x,y,\mathrm{\α})=\sqrt{{\left|D^{\mathrm{\α}}{f}_z(x,y)\right|}^2+{\left|D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)\right|}^2},$

该梯度相应的沿水平方向的辐角为：

${\mathrm{Af}}_z(x,y,\mathrm{\α})=\arg \mathrm{ument}(D^{\mathrm{\α}}{f}_z(x,y)+D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)),$

边缘点就是模M[f_z](x,y,α)沿辐角方向Af_z(x,y,α)有局部极大值的点。

针对每一向上延后的高度z，选择不同的方向α进行计算，可突出不同方向的边缘信息。为了能够达到二维平面完整覆盖，各方向α取值为kπ/(2^n-1)，其中k＝0,1,2,…,(2ⁿ-1)，n为大于或等于2的整数。针对每一方向，梯度的模的局部极大值点沿梯度的垂直方向连接得到的曲线构成边缘。在同一上延后的高度下，实现完整覆盖二维平面的各方向α计算边缘，对计算得到的各边缘求并集得到该尺度的边缘。

对每一尺度的位场进行上述计算，可实现各尺度下位场多方向的边缘检测。

步骤103，采用形态学骨骼算法对计算得到的各尺度的边缘进行细化处理成单像素宽度，得到各个尺度的格架图。

对上述计算得到的边缘图像采用Lam,L.,Seong-Whan Lee,and Ching Y.Suen,Thinning Methodologies-A Comprehensive Survey,IEEE Transactions onPattern Analysis and Machine Intelligence,v.14,no.9,September 1992的细化格架算法将边缘细化成单像素点宽度。

通过对计算得到的边缘进行细化，边缘及不同方向边缘交点所对应的实际地理范围较现有技术明显缩小，使得到的位场构造格架更接近于实际地质填图识别构造的特点，一方面图面信息清晰，增强了可读性，另一方面也便于进行地质解释。

步骤104，将计算得到的各尺度构造格架叠置生成综合构造格架图。

位场上延不同高度后提取的各尺度的边缘叠置，得到待研究区域的位场综合构造格架图。

可将位场上延后的不同高度对应于不同源区深度，见作者为Jacobsen,B.H.，题目为A case for upward continuation as a standard separation filter forpotential-field maps，Geophysics，v.52no.8，1987的文章，得到对应地球不同深度的构造。用渐变的颜色表示相应不同尺度或不同深度的格架，并将用不同的渐变颜色表征的各尺度构造格架叠加以突出显示不同深度的构造信息。通过如此叠加得到反映不同深度信息的综合构造格架图。

步骤105，用各尺度格架上各边缘点的梯度的模反映边缘强度。

针对不同的方向α，各边缘点处梯度的模M[f_z](x,y,α)大小不变，为一固定值，取边缘点处梯度的模代表该边缘点处构造强度值，建立不同尺度或不同深度构造强度格架图。

步骤106，将所述多个尺度的构造强度格架图叠置生成综合构造强度格架图。

对各尺度或深度边缘的强度值进行不同渐变颜色的叠加显示，可突出显示磁性和密度变化强度信息。

实例

下面以采用云南西部航磁数据进行位场构造格架自动提取为例，对本发明的技术方案进行解释。

首先，对航磁数据进行分块化极处理，将处理后的数据拼接成位场网格文件，网格大小为500米。本实例中的航磁数据为历史多次分别采集的航磁数据，测量高度在800-1200米范围内。

随后，将拼接的位场网格文件分别进行向上延拓处理得到高度上延后的位场数据，上延高度分别为1000、1500、2000、2500、3000、4000、5000、10000、15000、20000、25000、30000米，针对每一上延后高度进行32个方向的边缘检测，各方向α取值为kπ/(2^n-1)，其中k＝0,1,2,…,(2ⁿ-1)，n＝5，得到各尺度的边缘。

随后，采用骨骼算法对计算得到的各尺度的边缘进行细化处理，得到各个尺度的构造格架图。

图2为上延5000米上延高度边缘经骨骼算法细化提取的相应一个尺度的单像素宽度构造格架。可以看出，根据本发明的自动提取得到的单像素点宽度的位场构造格架更接近于实际地质填图识别构造特点，便于进行地质解释。另外图面信息更清晰，便于不同尺度格架的叠加分析。

图3为将位场分别向上延拓1000、1500、2000、2500、3000、4000、5000、10000、15000、20000、25000、30000米上延高度后经32个方向边缘检测得到的各尺度的格架进行颜色渐变叠加形成的综合构造格架图。将上延后的各个高度对应于相应的源区深度，可将该综合构造格架图用于表征研究区域不同深度的构造格架信息。用从灰白色到黑色渐变色代表了从低到高的向上延拓高度或从浅入深的深度，综合构造格架图反映了不同深度的构造信息。

图4为上延5000米上延高度格架的边缘点处的梯度的模所反应的构造强度，从灰白到黑色反应了梯度的模所代表的构造强度值逐渐增大。该图反映了相应一个尺度下构造强度变化。

图5为将位场分别上延1000、1500、2000、2500、3000、4000、5000、10000、15000、20000、25000、30000米上延高度后得到的边缘强度进行叠加得到的综合构造强度格架图。从图中可以看到，在不同尺度上，反映出主要构造带均表现为构造强度较大，对应的磁异常变换明显，指示了主要构造带均为磁异常的突变带，同时对应的深度非常大，反映了构造带控制了深部岩浆-成矿作用。

另外，图3和5中，具有深度和强度信息的构造格架图能够清楚的反映出区域中深度大、延展长的主干构造，和深度较浅、延展相对短的次级构造，以及互相的交割关系，因此根据本发明的方法得到的构造格架图可以帮助本领域技术人员认识所探测的区域构造格局。

将图3和图5与该区域的地表填图地质构造对比，空间位置和范围非常吻合，说明了根据本发明的自动提取方法的准确性和有效性。此外，与区域的地表填图地质构造对比，本发明的图3和图5增加了构造三维延伸、强度、结构等信息，可以帮助识别隐伏构造带。

沿大深度、高构造强度的构造带附近，不同方向构造的交汇部位以及构造带转折弯曲部位是发现潜在金属矿床的重要位置，根据本发明的方法利用航磁数据和重力数据可以快速准确提取得到具有深度和强度信息的构造格架图，帮助勘探人员准确快速地发现潜在的金属矿床。

以上借助优选实施例对本发明进行了详细说明，但是本发明不限于此。本技术领域技术人员可以根据本发明的原理进行各种修改。因此，凡按照本发明原理所作的修改，都应当理解为落入本发明的保护范围。

Claims (10)

1.一种位场构造格架自动提取方法，包括以下步骤：

对来自待研究区域的位场数据进行预处理；

将经预处理的位场数据向上延拓多个预定高度得到多个相应尺度的位场数据；

分别针对每一尺度的位场数据进行多方向边缘检测，得到多个相应尺度的位场边缘；

采用形态学骨骼算法将计算得到的各尺度的位场边缘分别细化为单像素宽度，得到多个相应尺度的构造格架图。

2.如权利要求1所述的位场构造格架自动提取方法，其特征在于，该方法进一步包括将计算得到的所述多个相应尺度的构造格架图叠置生成综合构造格架图。

3.如权利要求1所述的位场构造格架自动提取方法，其特征在于，该方法进一步包括：

将各尺度构造格架图上的每一边缘点处的梯度的模作为该尺度构造格架图中该边缘点处的强度值，得到多个相应尺度的构造强度格架图。

4.如权利要求3所述的位场构造格架自动提取方法，其特征在于，该方法进一步包括将所述多个相应尺度的构造强度格架图叠置生成综合构造强度格架图。

5.如权利要求1所述的位场构造格架自动提取方法，其特征在于，所述位场数据为重力位场数据或磁法位场数据，所述预处理进一步包括

对重力数据进行预处理得到布格重力异常；或

对磁法数据进行化极计算得到化极磁异常或进行伪重力计算得到伪重力异常。

6.如权利要求1所述的位场构造格架自动提取方法，其特征在于，

针对每一尺度的位场数据进行多方向边缘检测，包括以下步骤：

设尺度s＝z/z₀，且z>z₀，z₀代表测量高度，z代表向上延拓后的高度，定义高度为零的位置(x,y)处的重力异常或磁异常为f₀(x,y)，

尺度s的平滑函数定义为：

${\mathrm{\θ}}_s(x,y)=s^{-2}\mathrm{\θ}(x/s+y/s)={\mathrm{\γ}}_{{\mathrm{sz}}_0}(x,y),$

其中k(x,y,z)为格林函数，

在方向α的小波函数定义为：

${\mathrm{\ψ}}^{\mathrm{\α}}=(\cos \mathrm{\α}\frac{\∂}{\∂x}+\sin \mathrm{\α}\frac{\∂}{\∂y})\mathrm{\θ}(x,y)=D^{\mathrm{\α}}\mathrm{\θ}(x,y),$

其中，D表示一阶导数；

在尺度s和位置(x,y)的情况下，重力异常或磁异常f₀(x,y)在方向α的小波变换定义为：

$\begin{array}{c}{W}^{\mathrm{\α}}\left[f_0\right](x,y,s)=[f_0*{\mathrm{\ψ}}_s^{\mathrm{\α}}](x,y)\\ =f_0*\left(s{D}^{\mathrm{\α}}{\mathrm{\θ}}_s\right)\\ =s{D}^{\mathrm{\α}}[f_0*{\mathrm{\θ}}_s]\\ =s{D}^{\mathrm{\α}}[f_0*{\mathrm{\γ}}_{{\mathrm{sz}}_0}]\end{array}$

其中，*表示卷积运算，

由位场向上延拓公式已知

f_z(x,y)为f₀(x,y)从高度零上延到高度z＝sz₀处的重力异常或磁异常，通过将在测量高度z₀测量得到的重力异常或磁异常向上延拓高度z-z₀获得，

因此，W^α[f₀](x,y,s)＝sD^αf_z(x,y)

＝(z/z₀)D^αf_z(x,y)；

进一步，在尺度s和位置(x,y)，重力异常或磁异常f₀(x,y)在方向的小波变换定义为：

$W^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}\left[f_0\right](x,y,s)=(z/z_0)D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)$

则，f₀(x,y)的二维方向小波变换用梯度表示为：

$W\left[f_0\right](x,y,s,\mathrm{\α})=\left(\begin{array}{c}{W}^{\mathrm{\α}}\left[f_0\right](x,y,s)\\ W^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}\left[f_0\right](x,\mathrm{ys})\end{array}\right)=(z/z_0)\left(\begin{array}{c}{D}^{\mathrm{\α}}{f}_z(x,y)\\ D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)\end{array}\right)$

$W\left[f_0\right](x,y,s,\mathrm{\α})=(z/z_0)\▿f_z(x,y,\mathrm{\α})$

其中为二维梯度，

对于高度z,定义梯度的模为：

$M\left[f_z\right](x,y,\mathrm{\α})=\sqrt{{\left|D^{\mathrm{\α}}{f}_z(x,y)\right|}^2+{\left|D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x.y)\right|}^2},$

该梯度相应的沿水平方向的辐角为：

${\mathrm{Af}}_z(x,y,\mathrm{\α})=\arg \mathrm{ument}(D^{\mathrm{\α}}{f}_z(x,y)+D^{\mathrm{\α}+\frac{\mathrm{\π}}{2}}{f}_z(x,y)),$

边缘点为模M[f_z](x,y,α)沿辐角方向Af_z(x,y,α)有局部极大值的点，

针对每一方向α，梯度的模的局部极大值点沿梯度的垂直方向连接得到的曲线构成边缘，

针对同一高度，以多个不同的方向α计算边缘，对计算得到的各边缘求并集得到相应尺度的位场边缘。

7.如权利要求6所述的位场构造格架自动提取方法，其特征在于，针对同一高度以多个不同的方向α计算边缘的步骤进一步包括，各方向α取值为kπ/(2^n-1)，其中k＝0,1,2,…,(2ⁿ-1)，n为大于或等于2的整数，以完整覆盖二维平面。

8.如权利要求6所述的位场构造格架自动提取方法，其特征在于，位场向上延拓多个预定高度后提取的边缘对应于不同深度的构造，对得到的各尺度的构造格架图进行叠加得到反映不同深度信息的综合构造格架图。

9.如权利要求2或8所述的位场构造格架自动提取方法，其特征在于，以渐变颜色表征高度的高低，来形成所述综合构造格架图。

10.如权利要求4所述的位场构造格架自动提取方法，其特征在于，以渐变颜色表征强度值的大小，来形成所述综合构造强度格架图。