CN109254327B - 三维强磁性体的勘探方法及勘探系统 - Google Patents

Abstract

本发明提出了一种三维强磁性体的勘探方法和勘探系统，所述三维强磁性体的勘探方法包括：设置目标区域、三维强磁性体的展布范围和三维强磁性体的磁化率分布数据；根据所述目标区域、所述展布范围、所述磁化率分布数据和三维强磁性体的磁场计算方法，计算所述三维强磁性体的磁场；以及若所述磁场和利用磁力仪测量得到的所述三维强磁性体的实际磁场相同，将所述磁化率分布数据作为所述三维强磁性体的实际磁化率分布数据以用于勘探所述三维强磁性体。本发明实施例解决了目前三维强磁性体的勘探方法和勘探系统难以平衡计算效率和计算精度、无法满足大规模磁场数据精细反演成像需求等问题。

Description

三维强磁性体的勘探方法及勘探系统

技术领域

本发明涉及矿物勘探技术领域，尤其涉及一种用于金属矿物勘探、且适用于任意几何形状和任意磁化率分布的三维强磁性体的勘探方法和勘探系统。

背景技术

随着社会经济的飞速发展，我国对矿产资源的需求也日益增大。金属矿类是我国主要的矿产资源之一，由于其磁性性质与围岩存在较大差异，磁法成为矿产资源勘探的重要手段，被广泛运用于寻找金属矿物。

在磁场数据的定量和定性解释中，往往会涉及不同形状磁性体的磁异常计算问题。传统的解析法只能求解少数规则地质体的异常磁场。对于任意复杂形状的地质体的磁异常，通常需要借助数值方法。此外，已有的绝大多数研究仅考虑弱磁性体的磁异常问题，而对强磁性体的磁异常问题的关注相对较少。随着磁铁矿等金属矿物需求的日益增加，以及磁场数据在地质填图过程中的重要性，以上问题已经引起了越来越多学者的重视。文献(Sharma P V.Rapid computation of magnetic anomalies and demagnetizationeffects caused by bodies of arbitrary shape.Pure and Applied Geophysics,1966,64(1):89-109.)通过将复杂磁性体剖分成许多小棱柱体，假设每个小棱柱体内磁性均匀，利用积分方程求解每个小棱柱体单元的有效磁化强度，以此根据体积分法模拟三维复杂地质体的磁异常。但对于规模较大的磁性体，这种方法需要大量的小棱柱体单元来逼近，将导致其计算效率大大降低。文献(方华竹.任意形状强磁性三度体(三维磁性体)磁异常的计算.地质学报,1978,1:63-78.)提出了一种计算稳定磁化场中多个任意形状强磁性均质三度体表面有效磁荷面密度磁异常的方法。为了解决高磁化率问题，文献(Eskola L,TervoT.Solving the magnetostatic field problem(a case of high susceptibility)bymeans of the method of subsections.Geoexploration,1980,18(2):79-95.)和文献(Furness P.A versatile integral equation technique for magneticmodeling.Journal of Appllied Geophysics,1999,41(4):345-359.)提出了一种面积分数值方法，并同时考虑了感应磁化和剩余磁化的影响。文献(Traynin P,Hansen RO.Magnetic modeling for highly permeable bodies with remanentmagnetization.Seg Technical Program Expanded Abstracts,1993,12(12):1396.)采用有限元法求解积分方程，给出了一种求解三维任意形状地质体静磁问题的数值方法。文献(Purss M B J,Cull J P.A new iterative method for computing the magnetic fieldat high magnetic susceptibilities.Geophysics,2005,70(5):53-62.)提出了一种计算任意复杂形状的高磁化率地质体的正演模拟方法，该技术采用由均匀的任意直径的球形体素定义分段模型，这些球形体素用于描述磁偶极子，其中每个磁偶极子的磁矩是等效体积的球形体素的磁矩，磁场计算包括对自由空间中每个球形体素的中心磁场的迭代计算，然后对周围的球形体素迭代计算，从而考虑了各个球形体素之间的磁相互作用。文献(Wallace Y.3D modeling of banded iron formation incorporatingdemagnetization—A case study at the Musselwhite Mine,Ontario,Canada.Exploration Geophysics,2007,38(4):1-5.)采用边界元法模拟了三维带状铁矿的退磁效应。

目前，三维强磁性体的勘探方法和勘探系统中所用到的磁异常数值模拟方法往往存在计算效率和计算精度难以平衡的问题。因此，本发明急需提出一种三维强磁性体的勘探方法和勘探系统，以解决现有技术中存在的问题。

发明内容

考虑到目前大部分磁异常正演模拟问题仅研究弱磁性体的影响，并且现有三维强磁性体的勘探方法和勘探系统大多存在计算效率和计算精度低等问题，本发明实施例旨在提供一种高效、高精度的三维强磁性体的勘探方法和勘探系统，以满足磁场数据精细反演成像的需求，提高磁测资料及解释精度。

一方面，本发明实施例提供一种三维强磁性体的勘探方法，包括以下步骤：

设置目标区域、三维强磁性体的展布范围和三维强磁性体的磁化率分布数据；

根据所述目标区域、所述展布范围、所述磁化率分布数据和三维强磁性体的磁场计算方法，计算所述三维强磁性体的磁场；以及

若所述磁场和利用磁力仪测量得到的所述三维强磁性体的实际磁场相同，将所述磁化率分布数据作为所述三维强磁性体的实际磁化率分布数据以用于勘探所述三维强磁性体；

其中，所述三维强磁性体的磁场计算方法包括：

(a)根据所述目标区域和所述展布范围建立初始三维棱柱体模型；

(b)将所述初始三维棱柱体模型划分成多个棱柱体；

(c)根据所述磁化率分布数据对每个所述棱柱体的磁化率进行赋值，得到对应所述三维强磁性体的目标三维棱柱体模型；

(d)设置高斯参数；

(e)根据所述目标三维棱柱体模型和所述高斯参数计算得到高斯偏移波数；

(f)根据所述目标三维棱柱体模型和所述高斯偏移波数计算格林函数单元积分、所述格林函数单元积分的导数以及波数域系数：

(g)根据地球主磁场模型，计算每个所述棱柱体中心处的地球主磁场；以及

(h)将每个所述棱柱体中心处的所述地球主磁场作为其磁场初始值；

(i)根据所述目标三维棱柱体模型、所述磁场初始值、所述高斯偏移波数、所述格林函数单元积分、所述格林函数单元积分的导数以及波数域系数计算得到空间域异常磁场；

(j)根据所述磁场初始值和所述空间域异常磁场计算得到总磁场；

(k)若所述总磁场满足给定迭代收敛条件，将所述总磁场作为所述三维强磁性体的所述磁场；

(l)若所述总磁场不满足所述给定迭代收敛条件，将所述总磁场作为所述磁场初始值并重复执行步骤(i)至步骤(l)。

在本发明其中一个实施例中，所述棱柱体为规则棱柱体。

在本发明其中一个实施例中，在每个所述棱柱体内的所述磁化率为常量。

在本发明其中一个实施例中，所述高斯参数包括高斯点的数目、每个所述高斯点及所述高斯点对应的高斯系数。

在本发明其中一个实施例中，所述高斯点的数目为4，相应的高斯点t_g和高斯系数a_g分别为：

在本发明其中一个实施例中，所述给定迭代收敛条件为：

其中，N_x,N_y和N_z分别为所述三维棱柱体模型中x,y和z方向上的所述棱柱体的数量，H⁽⁰⁾为所述磁场初始值，H⁽¹⁾为所述总磁场，ε为相对均方根误差，ε₀为期望的数值精度。

另一方面，本发明实施例提供一种三维强磁性体的勘探系统，包括：存储器和连接所述存储器的一个或多个处理器，所述存储器存储有计算机程序代码，所述处理器用于执行所述计算机程序代码以实现如前述任一实施例所述的三维强磁性体的勘探方法。

本发明实施例具有以下优点或有益效果：本发明实施例的三维强磁性体的勘探方法和勘探系统通过利用具有简单、灵活等优点的棱柱体单元对复杂的三维强磁性体进行精细描述，能够十分方便地运用于磁化率分布复杂的三维强磁性体以及起伏地形情况的离散化；本发明实施例的三维强磁性体的勘探方法和勘探系统不仅能够高效、高精度地适用于计算磁化率较大的也即强磁性的三维强磁性体的磁场，从而提高利用磁场数据进行精细反演和解释的精度，还能够适用于三维弱磁性体，；另外，由于各个高斯偏移波数是相互独立的，在进行大规模的三维强磁性体的勘探时，本发明实施例的三维强磁性体的勘探方法和勘探系统具有占用计算机内存少、并行性好等优点。

附图说明

图1为本发明实施例提供的三维强磁性体的磁场计算方法的流程示意图；

图2为本发明实施例中球形三度体模型的示意图；

图3为本发明实施例中磁场分量数值解和解析解的对比图；

图4为本发明实施例中磁场分量的绝对误差图和迭代收敛趋势图；

图中符号说明如下：

Bx：磁场x分量，单位为nT；

By：磁场y分量，单位为nT；

Bz：磁场z分量，单位为nT；

本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

本发明实施例提供了一种三维强磁性体的勘探方法，主要包括：

设置目标区域、三维强磁性体的展布范围和三维强磁性体的磁化率分布数据；

根据所述目标区域、所述展布范围、所述磁化率分布数据和三维强磁性体的磁场计算方法，计算所述三维强磁性体的磁场；以及

若所述磁场和利用磁力仪测量得到的所述三维强磁性体的实际磁场相同，将所述磁化率分布数据作为所述三维强磁性体的实际磁化率分布数据以用于勘探所述三维强磁性体。下面详细介绍其中的三维强磁性体的磁场计算方法的执行原理：

参照图1，其为本发明实施例提供的一种三维强磁性体的勘探方法中的三维强磁性体的磁场计算方法的流程示意图，主要包括以下步骤：

第一步，复杂三维强磁性体模型表示：

首先，根据设置的所述目标区域及所述展布范围，建立三维棱柱体模型。选取空间中任意一点为坐标原点建立三维直角坐标系Oxyz，并确定三维棱柱体模型在x,y和z方向的起始位置。

然后，利用一系列与x,y,z轴平行的直线将此三维棱柱体模型均匀地划分成若干规则小棱柱体，小棱柱体的几何中心坐标为(x_i,y_j,z_n)，在x,y和z方向的尺寸分别为Δx,Δy和Δz_n，这里x和y方向采用等间隔剖分，z方向单元长度可以是任意的，Δz_n为z方向上第n个小棱柱体的垂向长度。模型x,y和z方向上小棱柱体的数量分别为N_x,N_y和N_z。

最后，根据设置的所述磁化率分布数据，对每个小棱柱体的磁化率进行赋值，在每个小棱柱体内，磁化率为常值，因此可以用小棱柱体中心(x_i,y_j,z_n)处的磁化率代表整个小棱柱体的磁化率。

第二步，高斯傅里叶变换的偏移波数计算：

设置高斯点数量N_G＝4，以及相应的高斯点t_g及高斯系数a_g，其中g＝1,2,…,4，即

根据x,y方向的剖分参数和高斯点，计算高斯傅里叶变换的偏移波数：

式(2)和(3)中，k_x和k_y分别为x,y方向的偏移波数；Δk_x和Δk_y分别为x,y方向的基波数；若p和q为偶数，则有

若p和q为奇数，则有

第三步，格林函数单元积分I_Gr(k_x,k_y,z_n,z_obs)及其导数I′_Gr(k_x,k_y,z_n,z_obs)和波数域系数W(k_x,k_y)的计算：

其中

i表示虚数单位；z_obs表示观测点其z方向坐标；

第四步，设置地球主磁场：

根据地球主磁场模型，计算每个小棱柱体中心(x_i,y_j,z_n)处的地球主磁场H_b，其中：i＝1,2,…,N_x,j＝1,2,…,N_y和n＝1,2,…,N_z，此处N_x,N_y和N_z分别为三维棱柱体模型x,y和z方向上小棱柱体的数量。

第五步，迭代计算三维强磁性体的磁场。具体步骤包括：

步骤(1)三维棱柱体模型的磁场初值设置：

将地球主磁场作为三维棱柱体模型的磁场初值，即

H⁽⁰⁾(x_i,y_j,z_n)＝H_b(x_i,y_j,z_n) 式(7)

步骤(2)小棱柱体的空间域磁化强度计算：

M(x_i,y_j,z_n)＝χ(x_i,y_j,z_n)H⁽⁰⁾(x_i,y_j,z_n) 式(8)

式(8)中，χ(x_i,y_j,z_n)表示以(x_i,y_j,z_n)为中心的小棱柱体的磁化率。

步骤(3)对空间域磁化强度进行二维傅里叶变换，即：

步骤(4)计算波数域异常磁场

的三个分量：

步骤(5)：对波数域异常磁场

进行二维傅里叶逆变换，得到空间域异常磁场H_a。然后，计算总磁场H⁽¹⁾＝H_a+H_b。

步骤(6)：判断H⁽¹⁾是否满足迭代收敛条件：

式(11)中，ε₀为期望的数值精度。若H⁽¹⁾满足收敛条件，则输出结果，否则取H⁽⁰⁾＝H⁽¹⁾，重复步骤(2)至(6)，直至满足迭代收敛条件为止。

下面对本发明实施例提供的三维强磁性体的磁场计算方法的效果进行检验。

为了说明本发明实施例所提出的三维强磁性体的磁场计算方法应用于计算任意几何形状、任意磁化率分布的复杂三度体磁性异常的效率和精度，设计了如图2所示的球形三度体模型，具体细节介绍如下：

三维棱柱体模型内有一个强磁性异常球体，三维棱柱体模型的范围为：x方向从0m到1000m，y方向从0m到1000m，z方向从0m到1000m(z轴向下为正)；强磁性异常球体的球心与模型中心重合，空间坐标为(500m,500m,500m)，球体半径为200m；磁化率为5；目标区域地球主磁场大小为50000nT，磁倾角为45度，磁偏角为30度。将三维棱柱体模型剖分成100*100*100个大小相同的小棱柱体，计算地面(z_obs＝0)处(图2中棱柱体的上部平面)的磁场，计算点数量为100*100。期望相对均方根误差大小|ε₀|＝10^-3％。

本发明实施例的三维强磁性体的磁场计算方法利用Fortran语言实现，运行平台的配置为：CPU为i7-2600，主频为3.4GHz，内存为32GB。对于100*100*100的模型，本发明实施例的三维强磁性体的磁场计算方法迭代一次所需时间约为20秒，相对均方根误差收敛至10^-3％需要迭代16次，由此可见效率较高。磁场x分量、y分量和z分量的数值解和解析解对比如图3所示，从形态上看，两者是一致的。考虑磁场零值点的影响，如图4所示，使用数值解与解析解的绝对误差衡量本发明实施例的数值精度。同时，图4还展示了三个磁场分量的迭代收敛趋势。

本发明实施例还提供了一种三维强磁性体的勘探系统，包括：存储器和连接所述存储器的一个或多个处理器，所述存储器存储有计算机程序代码，所述处理器用于执行所述计算机程序代码以实现如前述实施例所述的三维强磁性体的勘探方法。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (3)

1.一种三维强磁性体的勘探方法，其特征在于，包括：

设置目标区域、三维强磁性体的展布范围和三维强磁性体的磁化率分布数据；

根据所述目标区域、所述展布范围、所述磁化率分布数据和三维强磁性体的磁场计算方法，计算所述三维强磁性体的磁场；以及

若所述磁场和利用磁力仪测量得到的所述三维强磁性体的实际磁场相同，将所述磁化率分布数据作为所述三维强磁性体的实际磁化率分布数据以用于勘探所述三维强磁性体；

其中，所述三维强磁性体的磁场计算方法包括：

(a)根据所述目标区域和所述展布范围建立初始三维棱柱体模型；选取空间中任意一点为坐标原点建立三维直角坐标系Oxyz，并确定三维棱柱体模型在x,y和z方向的起始位置

(b)将所述初始三维棱柱体模型均匀地划分成多个规则棱柱体，小棱柱体的几何中心坐标为(x_i,y_j,z_n)，在x,y和z方向的尺寸分别为Δx,Δy和Δz_n，初始三维棱柱体模型在x,y和z方向上小棱柱体的数量分别为N_x,N_y和N_z；

(c)根据所述磁化率分布数据对每个所述棱柱体的磁化率进行赋值，得到对应所述三维强磁性体的目标三维棱柱体模型；

(d)设置高斯参数：高斯点数量N_G＝4，以及相应的高斯点t_g及高斯系数a_g，其中g＝1,2,…,4，即

(e)根据所述目标三维棱柱体模型和所述高斯参数计算得到高斯偏移波数；

式(2)和(3)中，k_x和k_y分别为x,y方向的偏移波数；Δk_x和Δk_y分别为x,y方向的基波数；

若p和q为偶数，则有

若p和q为奇数，则有

(f)根据所述目标三维棱柱体模型和所述高斯偏移波数计算格林函数单元积分I_Gr(k_x,k_y,z_n,z_obs)、所述格林函数单元积分的导数I′_Gr(k_x,k_y,z_n,z_obs)以及波数域系数W(k_x,k_y)：

其中

i表示虚数单位；z_obs表示观测点其z方向坐标；

(g)根据地球主磁场模型，计算每个所述棱柱体中心(x_i,y_j,z_n)处的地球主磁场H_b，其中：i＝1,2,…,N_x,j＝1,2,…,N_y和n＝1,2,…,N_z，此处N_x,N_y和N_z分别为三维棱柱体模型x,y和z方向上小棱柱体的数量；(h)将每个所述棱柱体中心处的所述地球主磁场作为其磁场初始值；

(i)根据所述目标三维棱柱体模型、所述磁场初始值、所述高斯偏移波数、所述格林函数单元积分、所述格林函数单元积分的导数以及波数域系数计算得到空间域异常磁场；

(j)根据所述磁场初始值和所述空间域异常磁场计算得到总磁场；

(k)若所述总磁场满足给定迭代收敛条件，将所述总磁场作为所述三维强磁性体的所述磁场；

(l)若所述总磁场不满足所述给定迭代收敛条件，将所述总磁场作为所述磁场初始值并重复执行步骤(i)至步骤(l)。

2.如权利要求1所述的三维强磁性体的勘探方法，其特征在于，所述给定迭代收敛条件为：

其中，N_x,N_y和N_z分别为所述三维棱柱体模型中x,y和z方向上的所述棱柱体的数量，H⁽⁰⁾为所述磁场初始值，H⁽¹⁾为所述总磁场，ε为相对均方根误差，ε₀为期望的数值精度。

3.一种三维强磁性体的勘探系统，其特征在于，包括：存储器和连接所述存储器的一个或多个处理器，所述存储器存储有计算机程序代码，所述处理器用于执行所述计算机程序代码以实现如权利要求1-2任意一项所述的三维强磁性体的勘探方法。

Images