一种磁全张量梯度反演地下三维磁性分布的波数域方法与流程

本发明涉及地球物理勘探技术，具体是一种磁全张量梯度反演地下三维磁性分布的波数域方法。

背景技术：

磁力方法是地球物理方法的主要分支之一，长期在资源勘查、工程勘察与区域构造研究中发挥着重要作用。磁力位场的正演和反演是磁力数据处理及其相应地质解释的主要途径。所谓正演即已知地下场源体的几何参数和物性参数(磁化强度或磁化率)，通过数学计算获得观测面上的理论磁场异常。而反演与之相反，即已知观测面磁场异常，通过计算进而获取地下磁源体的空间三维展布以及物性分布。磁力位场的正演是反演的基础，磁力位场的反演则是解决地下地质异常体相关构造解释的重要目的。

磁全张量梯度与磁三分量、磁总场异常相比，其受环境磁场影响小、能够更高分辨率地刻画磁性体的几何形态以及地下半空间三维分布情况。通过磁全张量梯度异常数据能够获得更多的反映场源特点及细节信息，甚至可以直接反应某些测量对象的磁矩矢量信息，此外，磁全张量梯度数据处理可以很好的描述目标体的磁化方向、走向等信息，可以促进精细地质解释。因此,近年来磁全张量梯度测量越来越流行，应用渐广。

三维物性反演是当前重磁数据处理解释研究的热点和发展趋势，用于给出地下空间的三维物性分布。三维物性反演技术通过利用地面或航空等实测数据，推断出磁性体(金属矿体、火成岩体、隐伏爆炸物、隐伏金属文物等)物性参数、几何参数及其在地下半空间的分布情况，从而达到寻找目标地质体的目的。该技术定位能力强，分辨率高，已成为地质填图、区域地质构造研究、金属矿勘探、工程勘察、考古等领域的关键技术。该技术应用三维网格模型来模拟近似地下的三维全空间，即将地下三维全空间剖分为众多规则排列、大小相等的长方体体积元，但各长方体的物性值(密度、磁化强度或磁化率)并不相同且是未知的、待反演计算的。三维物性反演主要通过反复的三维网格模型正演计算及迭代拟合来获取地下三维物性分布。

三维网格模型的正演和反演一般采用空间域算法，正演即利用空间域的数学解析式正演计算三维网格模型每个长方体在观测面上任意测点引起的重磁场，而观测面上每个测点的重磁场是由正演计算得到的所有立方体在该点引起的重磁场的累加求和，而反演则是正演的逆过程，即根据已知异常数据，通过反演，进而得到地下半空间磁性体分布情况。空间域算法具有计算简单、容易实现等优点，但同时也存在计算量大、效率低下、内存占用率高、存在边缘效应等缺点，直接影响反演效率和效果。波数域算法是地球物理数据处理常用的快速算法，与空间域算法相比，它具有计算简单、速度快的优点，在磁力数据处理解释中得到广泛应用。

总体来讲，传统三维磁力反演采用空间域算法，计算简单，但计算量大，对于大数据体模型反演效率低下。磁全张量梯度三维反演的波数域方法对提高反演精度和效率、促进精细地质解释具有重要意义。

技术实现要素：

为了快速重建地下空间的三维物性分布，通过利用地面或航空等实测磁全张量梯度异常数据，推断出磁性体(金属矿体、火成岩体、隐伏爆炸物、隐伏金属文物等)物性参数、几何参数及其在地下半空间的分布情况，从而达到寻找目标地质体的目的。该技术定位能力强，分辨率高，已成为地质填图、区域地质构造研究、金属矿勘探、工程勘察、考古等领域的关键技术。本发明提供一种磁全张量梯度反演地下三维磁性分布的波数域方法；本发明利用研究区内实测磁全张量梯度异常数据以及研究区内先验信息，通过对已知观测异常数据进行波数域迭代反演，不断更新磁模型数据，最终获得满足误差限的磁三维网格模型数据。

本发明是通过以下技术方案实现的：

一种磁全张量梯度反演地下三维磁性分布的波数域方法,所述方法首先对该研究区域观测磁全张量梯度异常进行反演，进而获得磁三维网格模型数据；然后，对模型数据进行正演，得到理论磁全张量梯度异常；将其与实际观测数据相减得到理论异常残差，若异常残差的最小范数解不满足误差限，则对异常残差进行反演得到异常残差模型，重复上述步骤，不断更新磁模型数据，最终获得满足误差限的磁三维网格模型数据，重建地下半空间磁性体的三维分布；

进一步地，所述方法包括如下步骤：

s1：根据该研究区的先验信息，获得该区磁倾角、磁偏角、磁性体物性范围及上下界面埋深；

s2：设置反演滤波参数、反衰减因子、反演深度范围及步长；

s3：对观测全张量梯度异常tobs进行波数域三维快速约束反演，得到磁三维网格模型数据m，并根据先验信息对其进行物性约束；

s4：再对得到的三维网格数据m进行波数域三维快速正演，进而获得磁全张量梯度异常tcal；

s5：将理论异常与观测异常相减，得到理论异常残差δt；

s6：判断异常残差δt的最小范数解是否满足误差限，若不满足，通过对异常残差δt反演得到异常残差模型δm，并更新模型m＝m+δm，再重复前述步骤，直至满足误差限，最终得到磁三维网格模型数据；

进一步地，所述s3中波数域三维快速反演的过程具体为：

s311:对该研究区观测磁全张量梯度异常tobs进行等间距网格化；

s312：根据两个方向上测点数量，选取接近最大值的2的整数次幂进行扩边；

s313：对扩边后的观测异常数据进行二维傅里叶变换，得到观测异常频谱；

s314：将观测异常频谱与反演滤波因子以及界面约束因子相乘，然后对其进行二维反傅里叶变换，并同时进行消边，最终获得三维网格模型数据tcal；

进一步地，所述s314中反演滤波因子可由如下公式计算得到：

其中：为磁场强度分布方向矢量；t为地磁场强度；kx，ky分别为x，y的波数；n为滤波函数斜率，n∈(0,10)；z为反演深度；n为反衰减因子，初始为1，在每次迭代中增加0.5；为全张量梯度因子；

进一步地，所述s314中界面约束因子可由如下公式计算得到：

其中：α＝0.001为界面约束因子最佳参数值；r为界面约束因子斜率；zrange为反演的深度范围；ztop，zbottom分别为磁性体上下界面埋深；z为反演深度；

进一步地，所述s3中物性约束过程具体为：

s321：根据该研究区内的先验信息，获得磁性体的磁化率范围(κ1～κ2)；

s322：对于磁三维网格模型m，遍历每个网格单元，当磁化率值小于κ1时，修改反演的磁化率值，使其磁化率强制等于κ1，当磁化率值大于κ2时，修改反演的磁化率值，使其磁化率强制等于κ2；

进一步地，所述s4中波数域三维快速正演的过程具体为：

s41：根据三维网格模型数据m三个方向上的测点数量，选取接近最大值的2的整数次幂进行扩边；

s42：对扩边后的三维网格模型进行三维傅里叶变换，得到三维网格模型频谱；

s43：将三维网格模型频谱与相应的正演滤波因子相乘，得到三维全张量梯度异常频谱；

s44：对s43中获得的三维全张量梯度异常频谱进行三维反傅里叶变换，并对其进行消边，进而获得空间域三维全张量梯度异常，提取相应观测面的理论三维全张量梯度异常tcal；

进一步地，所述s43中正演滤波因子可由如下公式计算得到：

其中：为磁场强度分布方向矢量；t为地磁场强度；kx，ky，kz分别为x，y，z的波数；为全张量梯度因子；

本发明的有益效果如下：

1)本发明使用实测磁全张量梯度异常进行反演，通过磁全张量梯度异常数据能够获得更多的反映场源特点及细节信息，提高计算精度；

2)本发明应用三维网格模型来模拟近似地下磁性体三维分布情况；

3)本发明采用波数域算法，它具有计算简单、速度快的优点，为解决大数据体磁三维反演提供有效方法；

4)本发明在反演迭代不断增加反衰减因子，减少反演滤波因子在深度方向的衰减性，同时降低“趋肤效应”，使反演结果更加接近地下磁性体真实分布情况；

5)本发明在波数域反演过程中，直接在反演滤波因子中引入界面约束因子，约束磁性体的上下界，并且大大提高了反演效率及反演精度；

6)本发明在迭代中引入直接物性约束法，限制物性值上下限，并在不断迭代中更新磁三维网格模型数据，使反演精度进一步提高；

7)本发明采用波数域的正、反演方法进行迭代，有利于更高效、更精确地获得地下半空间磁性体三维分布情况。

附图说明

图1中a和b为本发明实施例1中两种情况下组合直立长方体在地下半空间的三维分布情况图；

图2中a和b为本发明实施例1中两种情况下组合直立长方体的某一张量梯度异常图(uzz)；

图3中a和b为本发明实施例1中两种情况的波数域磁三维反演结果图；

图4为本发明磁全张量梯度快速三维反演的波数域方法流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细描述。应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。相反，本发明涵盖任何由权利要求定义的在本发明的精髓和范围上做的替代、修改、等效方法以及方案。

进一步，为了使公众对本发明有更好的了解，在下文对本发明的细节描述中，详尽描述了一些特定的细节部分。对本领域技术人员来说没有这些细节部分的描述也可以完全理解本发明。

下面结合附图和具体实施例对本发明作进一步说明，但不作为对本发明的限定。下面为本发明的举出最佳实施例：

实施例1

三维波数域反演是当前重磁数据处理解释研究的热点和发展趋势，用于快速重建地下半空间的三维物性分布。通过磁全张量梯度异常数据能够获得更多的反映场源特点及细节信息，有利于地下精细地质解释。本实施以两个组合直立长方体为例：两个直立长方体(1)物性相同、埋深相同，如图1中a所示；(2)物性不同，埋深不同，如图1中b所示。采用一种磁全张量梯度反演地下三维磁性分布的波数域方法，如图4所示，依次进行。所述方法包括如下步骤：

(1)根据模型数据先验信息，获得模型磁倾角、磁偏角、磁性体物性范围(以磁化率为例)及上下界面埋深；

(2)设置反演滤波参数、反衰减因子、反演深度范围及步长；

(3)对观测全张量梯度异常tobs进行波数域三维快速反演，得到磁三维网格模型数据m；

(4)再对得到的三维网格数据m进行波数域三维快速正演，进而获得理论磁全张量梯度异常tcal；

(5)将理论异常与观测异常相减，得到理论异常残差δt；

(6)判断异常残差δt的最小范数解是否满足误差限，若不满足，对异常残差进行反演得到异常残差模型δm，更新模型m＝m+δm；

(7)重复前述步骤，直至满足误差限，最终得到磁三维网格模型数据。

其中，步骤(3)中，所述波数域三维快速反演的过程具体为：

对该研究区观测磁全张量梯度异常tobs进行等间距网格化；

根据两个方向上测点数量，选取接近最大值的2的整数次幂进行扩边；

对扩边后的观测异常数据进行二维傅里叶变换，得到观测异常频谱；

将观测异常频谱与反演滤波因子以及界面约束因子相乘，并对其进行二维反傅里叶变换，同时进行消边，最终获得空间域三维网格模型数据m。

其中，步骤(3)中，所述反演滤波因子可由如下公式计算得到：

其中：为磁场强度分布方向矢量；t为地磁场强度；kx，ky分别为x，y的波数；n为滤波函数斜率，n∈(0,10)；z为反演深度；n为反衰减因子，初始为1，在每次迭代中增加0.5；为全张量梯度因子。

其中，步骤(3)中，所述界面约束因子可由如下公式计算得到：

其中：α＝0.001为界面约束因子最佳参数值；r为界面约束因子斜率；zrange为反演的深度范围；ztop，zbottom分别为磁性体上下界面埋深；z为反演深度。

其中，步骤(3)中，所述物性约束过程具体为：

根据该研究区内的先验信息，获得磁性体的磁化率范围(κ1～κ2)；对于更新后的m，遍历每个网格单元，当磁化率值小于κ1时，修改反演的磁化率值，使其磁化率强制等于κ1，当磁化率值大于κ2时，修改反演的磁化率值，使其磁化率强制等于κ2。

其中，步骤(4)中，所述波数域三维快速正演的过程具体为：

根据三维网格模型数据m三个方向上的测点数量，选取接近最大值的2的整数次幂进行扩边；

对扩边后的三维网格模型进行三维傅里叶变换，得到三维网格模型频谱；

将三维网格模型频谱与相应的正演滤波因子相乘，得到三维全张量梯度异常频谱；

对其进行三维反傅里叶变换，并对其进行消边，进而获得空间域三维全张量梯度异常，提取相应观测面的理论三维全张量梯度异常tcal。

其中，步骤(4)中，所述正演滤波因子可由如下公式计算得到：

其中：为磁场强度分布方向矢量；t为地磁场强度；kx，ky，kz分别为x，y，z的波数；为全张量梯度因子。

在本实施例所述方法中，在不断更新三维网格模型数据过程中，采用波数域的正、反演方法进行迭代，有利于更高效、更精确地获得地下半空间磁性体三维分布情况。

对上述方法进行模型测试：

步骤一、建立三维网格模型，网格大小为51×51×21,分别采用两种不同情况，进行反演：两个直立长方体(1)物性相同，磁化率均为0.05si；埋深相同，上表面埋深距观测面均为650m，下表面埋深1250m(图1-a)；(2)物性不同，磁化率分别为0.1si和0.05si；埋深不同，分别为上表面为650m，下表面为1250m，以及上表面为250m，下表面为750m(图1-b)。利用如下公式正演全张量梯度异常的其中一个分量(磁场延深度方向二阶导数)作为已知观测全张量梯度异常(图2-a、b)：

其中：为磁场强度分布方向矢量；t为地磁场强度；kx，ky，kz分别为x，y，z的波数；κ为磁化率，和分别表示傅里叶变换和反傅里叶变换。

步骤二、根据该组合模型已知参数，提取先验信息，即其物性范围(以磁化率为例)，以及上下界面埋深情况，将此先验信息加入反演过程中，对反演过程进行约束；

步骤三、设置反演滤波参数、反衰减因子、反演深度范围及步长；

步骤四、将观测异常进行全张量梯度快速三维波数域反演，得到初次磁三维网格模型m；

步骤五、再对得到的三维网格数据进行波数域三维快速正演，进而获得理论磁全张量梯度异常tcal；

步骤六、将理论异常与观测异常相减，得到理论异常残差δt；

步骤七、判断异常残差δt最小范数解是否满足误差限，若不满足，对异常残差进行反演得到异常残差模型δm，更新模型m＝m+δm；

步骤八、重复前述步骤，直至满足误差限，最终得到磁三维网格模型数据(图3-a、b)。如图3，可见通过磁全张量梯度快速三维反演的波数域方法，更高效、更精确地获得地下半空间磁性体三维分布情况。与模型设定值的误差在可接受范围内，因此，模型试验结果证明该方法有效。

以上所述的实施例，只是本发明较优选的具体实施方式的一种，本领域的技术人员在本发明技术方案范围内进行的通常变化和替换都应包含在本发明的保护范围内。