CN113139289B - 一种基于积分方程的退磁效应下磁测数据的正反演方法 - Google Patents

Abstract

本发明提供了一种基于积分方程的退磁效应下磁测数据的正反演方法，包括以下步骤：根据地下目标区域，确定网格剖分的空间范围，使用正交网格对所述空间范围进行网格剖分；设定地磁场的磁倾角、磁偏角和磁感应强度；设定观测位置，利用基于积分方程的在退磁效应下的磁测数据正演计算方法，能够获得考虑退磁效应下在观测位置的磁测数据；输入磁场观测数据，通过基于积分方程在退磁效应下的磁测数据反演方法，结合正则化项进行模型最优化计算，对地下目标进行反演计算，得到反演结果。本发明能计算考虑退磁效应下的磁测数据正演，并能够根据磁测数据进行考虑退磁效应的反演。

Description

一种基于积分方程的退磁效应下磁测数据的正反演方法

技术领域

本发明涉及地球物理探测领域，尤其涉及一种基于积分方程的退磁效应下磁测数据的正反演方法。

背景技术

在地磁探测中，人们通常为了减小计算的时间，将异常体的磁化强度进行近似，这样探测出来的异常体的磁异常数据往往与实际的数据解释存在一定的误差。在低磁背景下，这种误差可以忽略不计，但在高磁背景，该误差对磁场异常探测有很大的影响，不考虑退磁效应的处理结果精度较低，无法满足实际应用的需求。因此考虑退磁效应下的磁梯度张量数据不仅能很好的描述异常体的磁化方向、空间几何结构、地质体的走向等信息，而且在高磁背景下精度准确，能很好地解释地下地质信息。

现有文献1“Sharma P V.Rapid computation of magnetic anomalies anddemagnetization effects caused by bodies of arbitrary shape[J].Pure&AppliedGeophysics,1966,64(1):89-109.”提出了积分方程考虑退磁效应下的磁场正演方法，但不包括反演方法。文献2“Lelièvre,P.G.and Oldenburg,D.W.Magnetic forward modellingand inversion for high susceptibility.2006,166:76-90.”提出了偏微分方程下的磁总场退磁方法。文献3“Lelievre,Peter.(2021).Forward modeling and inversion ofgeophysical magnetic data.”使用有限体积、有限元和磁偶极子等方法比较了退磁效应下磁总场和磁三分量之间的磁异常特征。文献4“Krahenbuhl R A,Li Y.Investigation ofmagnetic inversion methods in highly magnetic environments under strong self-demagnetization effect[J].2017,82(6):1-60.”提出了积分方程的退磁效应下的磁总场磁异常正反演方法。文献5“Zhdanov M S,Cai H, Wilson G A.3D Inversion of SQUIDMagnetic Tensor Data[J].2012,1(1).”提出了磁梯度张量数据的三维正反演方法，显著提高了反演的分辨率，但未考虑退磁效应。文献6“Zhdanov M S,Wilson G A,Technoimaging.3D inversion of full tensor magnetic gradiometry(FTMG)data[C].2011.”利用全张量磁梯度测量来反演澳大利亚塔拉旺地区的一个磁铁矿，反演结果与当地的地质情况非常吻合。目前在全球范围内，地球物理行业内已有偏微分方程和积分方程的退磁效应下的磁总场和磁三分量三维正演研究，但是积分方程下考虑退磁效应的磁梯度张量数据三维正反演和磁总场、磁三分量数据三维反演方法未见提出。

发明内容

有鉴于此，针对以上技术问题，本发明提供的一种基于积分方程的退磁效应下磁测数据的正反演方法，该方法能够进行磁测数据为磁梯度张量数据的正反演计算以及磁测数据为磁总场、磁三分量数据的反演计算，具体包括以下步骤：

S1、根据地下目标区域，确定网格剖分的空间范围，使用正交网格对所述空间范围进行网格剖分，得到剖分后的网格；

S2、设定地磁场的磁倾角、磁偏角和磁感应强度；

S3、在剖分后的网格中设定观测位置，利用基于积分方程的退磁效应下的磁测数据正演计算方法，获得在退磁效应下观测位置的磁测数据；

S4、根据磁测数据，通过基于积分方程的退磁效应下磁测数据的反演方法，结合正则化项进行模型最优化计算，对地下目标进行反演计算，得到反演结果。

进一步地，步骤S3中，基于积分方程的退磁效应下的磁测数据正演计算方法，其计算公式为：

M＝K(H₀+TM)

AM＝b

A＝I-KT

b＝KH₀

B_ij＝μ₀T_ijM,i∈{x,y,z},j∈{x,y,z}

其中，M表示考虑退磁效应后的磁化强度矩阵；K表示磁化率矩阵，H₀表示磁场强度矢量,B_ij异常场的磁梯度张量数据,μ₀是真空中的磁导率；I为单位阵。

p的观测点位置位于每个正交网格单元的中心。

进一步地，步骤S4中，基于积分方程的退磁效应下磁测数据的反演方法的目标函数为：

minφ＝||F(κ)-d^obs||²+βΔκ

F(κ)＝Gκ

G＝T_ij(B₀+μ₀TM)

其中，φ表示优化目标，F(κ)表示正演过程计算出的磁测数据，磁测数据类型包括磁总场、磁三分量和磁梯度张量数据，d^obs表示实际观测磁测数据,β代表规整化因子，Δκ表示模型规整化项，G为灵敏度矩阵，κ表示模型磁化率矩阵， B₀为地磁场强度，μ₀表示真空磁导率；

本发明提出的技术方案能使用积分方程方法计算考虑退磁效应下的磁异常梯度张量数据，对高磁化率环境非常适用，同时本发明计算出来的磁正反演结果比不考虑退磁效应下计算出来的磁正反演结果误差相对较小，更加符合真实环境。

附图说明

图1是本发明一种基于积分方程的退磁效应下磁测数据的正反演方法的流程图；

图2是本发明的网格剖分原理示意图；

图3是本发明实施例提供的一个地下模型网格剖分示意图；

图4是本发明实施例提供的考虑不退磁效应和考虑退磁效应下的磁总场和磁三分量的磁异常数据示意图；

图5是本发明实施例提供的不考虑退磁效应和考虑退磁效应下磁梯度张量数据示意图；

图6是本发明实施例提供的不考虑退磁效应和考虑退磁效应下的反演结果示意图。

具体实施方式

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

请参考图1，一种基于积分方程的退磁效应下磁测数据的正反演方法，包括以下：

S1、根据地下目标区域，确定网格剖分的空间范围，使用正交网格对所述空间范围进行网格剖分，得到剖分后的网格。

网格剖分原理如图2所示，按照正交网格的形式进行剖分；实施例提供的模型剖分示意图如图3所示，选定X轴、Y轴、Z轴方向为200m×200m×200m 范围的区域作为网格剖分的空间范围，用来剖分的网格为10m×10m×10m的六面体形式，示意图中的不规则物体为磁化率为1SI的地下目标；

S2、设定地磁场的磁倾角、磁偏角和磁感应强度。常用的参数为磁倾角90 度，磁偏角0度，磁感应强度50000nT；

S3、在剖分后的网格中设定观测位置，利用基于积分方程的退磁效应下的磁测数据正演计算方法，获得在退磁效应下观测位置的磁测数据；

S4、根据磁测数据，通过基于积分方程的退磁效应下磁测数据的反演方法，结合正则化项进行模型最优化计算，对地下目标进行反演计算，得到反演结果。

步骤S3中，基于积分方程的退磁效应下的磁测数据正演计算方法，其计算公式为：

M＝K(H₀+TM)

AM＝b

A＝I-κT

b＝κH₀

B_ij＝μ₀T_ijM,i∈{x,y,z},j∈{x,y,z}

其中，M表示考虑退磁效应后的磁化强度矩阵,在反演算法的每次迭代中每次都会重新计算一次磁化强度；K表示磁化率矩阵，H₀表示磁场强度矢量,B_ij磁梯度张量异常数据,μ₀是真空中的磁导率,在SI制中μ₀＝4π×10^-7H/m；I为单位矩阵。T和T_ij的公式均为

区别在于T中p的观测点位置位于每个正交网格单元的中心，T_ij中p的观测点位置位于观测面。

为了便于理解，下面详细解释步骤S3中，基于积分方程的退磁效应下的磁测数据正演计算方法的推导过程。

由于剖分后每个网格单元的磁性体形状为正交网格，假定每个正交网格的磁化强度矢量为M＝(M₁,M₂,M₃)，则每个正交网格单元在观测点P的异常场磁场强度在i方向的上分量为：

公式(1)中

表示观测点p在i方向上异常场的磁场强度，p是观测面上观测点位置，

描述了场源和观测点磁场强度分量之间的关系，展开形式如公式(2)所示：

公式(2)中q是正交网格单元的中心点位置，同时，专利中考虑的磁性体为线性且不包含剩磁的，因此计算第p个正交网格单元的磁化强度三分量：

公式(3)中

为i方向上的总场磁场强度分量，H_0i为i方向上的地磁场强度分量，p代表第p个正交网格单元的测点位置，κ^p为第p个正交网格单元的磁化率。公式(3)写为矩阵方程的形式为：

M＝κ(H₀+TM) (4)

将公式(4)方程式移项可得

AM＝b

A＝I-κT

b＝κH₀

(5)

公式(5)中，M表示考虑退磁效应后计算出的磁化强度矩阵；κ表示磁化率矩阵，H₀表示地磁场强度矢量,I为单位矩阵。根据公式(5)可计算出正交网格单元的磁化强度矩阵M，同时，通过磁化场计算磁感应强度的公式为

公式(6)中B_ij表示磁梯度张量异常数据,即反演计算中的F(κ)，μ₀是真空中的磁导率；H_ij为磁梯度的磁场强度。

则根据公式(1)和公式(6)可得

B_ij＝μ₀T_ijM,i∈{x,y,z},j∈{x,y,z} (7)

步骤S4中，基于积分方程的退磁效应下磁测数据的反演方法的目标函数为：

minφ＝||F(κ)-d^obs||²+βΔκ

F(κ)＝Gκ

G＝T_ij(B₀+μ₀TM)

其中，F(κ)表示正演过程，d^obs表示实际观测磁测数据,磁测数据类型包括磁总场、磁三分量和磁梯度张量数据，β代表规整化因子，Δκ表示模型规整化项， G为灵敏度矩阵，κ表示模型磁化率矩阵，B₀为地磁场强度，μ₀表示真空磁导率。

在最优化计算的迭代过程中，每次迭代都会重新计算一次磁化强度M，进而重新计算预测的正演数据用于对d^obs进行拟合，M和磁梯度张量数据格式下的 F(κ)由以下的正演计算方法计算得来，即正演计算方法中的B_ij。

图4是本发明实施例提供的考虑不退磁效应和考虑退磁效应下的磁总场和磁三分量的磁异常数据示意图；

图4中B(nT)表示磁总场数据；Bx(nT)、By(nT)、Bz(nT)表示磁三分量的磁异常数据强度，坐标轴代表观测区域范围，纵坐标为X轴，横坐标为Y轴，单位为米，右侧色标表示磁异常数据的强度大小,单位为nT；

图5为发明实施例提供的不考虑退磁效应和考虑退磁效应下磁梯度张量数据示意图。

图5中，坐标轴代表观测区域范围，纵坐标为X轴，横坐标为Y轴，单位为米，图中从左至右，从上至下分别为磁梯度张量Bxx、Bxy、Bxz、Byy、Byz 和Bzz，右侧色标表示磁异常数据的磁梯度张量强度，单位为nT/m；

将计算得到的网格的磁化率强度进行成图，不考虑退磁效应和考虑退磁效应下的反演结果对比示意图如图6所示，可以看到考虑退磁效应下的反演结果更加收敛。

本发明的有益效果是：能使用积分方程法计算退磁效应下的磁梯度张量磁异常数据，对高磁化率环境非常适用；本发明计算出来的磁正反演结果比不考虑退磁效应下计算出来的磁正反演结果误差相对较小，更加符合真实环境。

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

Claims (1)

1.一种基于积分方程的退磁效应下磁测数据的正反演方法，其特征在于：

S1、根据地下目标区域，确定网格剖分的空间范围，使用正交网格对所述空间范围进行均匀网格剖分，得到剖分后的网格；

S2、设定地磁场的磁倾角、磁偏角和磁感应强度；

S3、在剖分后的网格中设定观测位置，利用基于积分方程的退磁效应下的磁测数据正演计算方法，获得在退磁效应下观测位置的磁测数据；

S4、根据磁测数据，通过基于积分方程的退磁效应下磁测数据的反演方法，结合正则化项进行模型最优化计算，对地下目标进行反演计算，得到反演结果；

步骤S3中，基于积分方程的退磁效应下的磁测数据正演计算方法，其计算公式为：

M＝K(H₀+TM)

AM＝b

A＝I-KT

b＝KH₀

B_ij＝μ₀T_ijM,i∈{x,y,z},j∈{x,y,z}

其中，M表示考虑退磁效应后的磁化强度矩阵；K表示磁化率矩阵，H₀表示磁场强度矢量,B_ij异常场的磁梯度张量数据,μ₀是真空中的磁导率；I为单位阵；

p的观测点位置位于每个正交网格单元的中心；

步骤S4中，基于积分方程的退磁效应下磁测数据的反演方法的目标函数为：

F(κ)＝Gκ

G＝T_ij(B₀+μ₀TM)

其中，F(κ)表示正演过程计算出的磁测数据，d^obs表示实际观测磁测数据,β代表规整化因子，Δκ表示模型规整化项，G为灵敏度矩阵，κ表示模型磁化率矩阵，B₀为地磁场强度，μ₀表示真空磁导率。

Images