CN112748471A - 一种非结构化等效源的重磁数据延拓与转换方法 - Google Patents

Abstract

本发明公开了一种非结构化等效源的重磁数据延拓与转换方法，包括：通过观测得到物理场数据d₀；给定反演初始模型和反演网格边界，利用现有的四面体网格生成方法生成非结构化四面体网格，作为等效源反演网格空间；根据地磁场参数，对物理场数据d₀进行带规整化项的PDE三维反演计算，同时对反演等效源模型周围的反演网格空间，进行加密剖分动态调整，在满足迭代条件后得到异常体的最终等效源模型；根据实验需要设定不同的正演参数，利用步骤S3得到的等效源模型进行相应的三维正演计算，得到等效源产生的上延拓数据或下延拓数据、化极数据、异常总场数据、异常分量数据以及梯度张量数据。本发明实现对复杂环境中的异常数据进行高精度延拓和数据转换。

Description

一种非结构化等效源的重磁数据延拓与转换方法

技术领域

本发明涉及地球物理勘测技术领域，尤其涉及一种非结构化等效源的重磁数据延拓与转换方法。

背景技术

在重力或磁法探测中，测量分布在起伏观测曲面上，由于测量位置与异常体的距离不同，重磁数据所反映的主要异常信息内容也有所差异。在实际应用中，又往往需要将这些数据转换成所需要的参数和类型，比如不同观测高度(延拓)、异常分量、张量数据的换算等，重磁数据类型转换的主要挑战在于航磁(地面磁测)的起伏测量曲面以及数据的不规则测量位置(非网格化测量)，传统的数据类型转换以及化极计算无法处理这些数据，而采用传统单层、双层、三层等效源的数据转换方法，处理结果精度较低，无法满足应用需求。

现有文献1“Thomas Günther,Carsten Rücker,Klaus Spitzer.Three-dimensional modelling and inversion of dc resistivity data incorporatingtopography–II.Inversion[J].Geophysical Journal International,2006,166(2).”提出了基于非结构化网格的电法反演方法，未见基于非结构化网格的重磁反演；文献2“Dampney,C.N.G.THE EQUIVALENT SOURCE TECHNIQUE[J].geophysics,1969,34(1):39.”提出了等效源方法，利用单层等效场源模拟实测异常，可用于位场数据的空间延拓(包括曲面延拓)、梯度计算以及分量转换等；文献3“管志宁,安玉林,陈维雄.曲线与曲面上磁场向上延拓和分量转换[J].地球物理学报,1985,000(004):419.”提出了基于单层位和双层位理论的曲面上磁场的上延拓方法；文献4“Li Y.Reduction to the pole usingequivalent sources[J].SEG Technical Program Expanded Abstracts,1999,19(1):2484.”通过在有限的深度内设置单层的等效源实现了低纬度磁场异常的化极方法。文献5“Zuo B,Hu X,Marcelo Leo-Santos,et al.Downward Continuation and Transformationof Total-Field Magnetic Anomalies Into Magnetic Gradient Tensors BetweenArbitrary Surfaces Using Multilayer Equivalent Sources[J].GeophysicalResearch Letters,2020,47(16).”提出了一种通过带深度规整化项的多层等效源技术进行延拓与数据转换的方法。文献“Farquharson C.G.,C.R.W.Mosher,2009,Three-dimensional modelling of gravity data using finite differences,Journal ofApplied Geophysics,68,417-422.”实现了六面体结构化网格的三维差分正演，未实现反演。文献“Jahandari H.and C.G.Farquharson,2013,Forward modeling of gravity datausing finite-volume and finite-element methods on unstructured grids,Geophysics,78(3),G69–G80”实现了四面体非结构化网格的有限体积、有限元的重力三维正演，未实现反演。文献“Galley C.G.,P.G.Lelièvre_and C.G.Farquharson,2020,Geophysical inversion for 3D contact surface geometry,Geophysics,85,K27-K45.”采用积分方程实现了非结构化网格的曲面积分反演，与PDE方法在基础理论与方程上显然不同。

然而，这些技术存在以下问题：1)结构化六面体(立方体)网格对模型空间和异常体的剖分精度与效率存在条件限制：一方面，采用立方体单元在对复杂三维结构进行剖分时，容易导致模型离散化误差过大，拟合精度低；另一方面，若采用大量精细立方体单元对复杂三维结构进行剖分，则会导致三维正反演计算量过大，对硬件性能及计算时间要求较高；2)对自由空间和模型空间进行结构化六面体的网格剖分，难以在提高观测点、地形、场源的模型剖分精度的同时减少模型单元个数；3)难以在航空数据处理中对观测点以下至地形范围的空间采用非均匀的稀疏网格进行表示。

发明内容

针对以上技术问题，本发明提供了一种基于非结构化等效源的重磁数据延拓与转换方法，利用现有的四面体网格生成方法生成非结构化网格，并基于非线性PDE(PartialDifferential Equations，偏微分方程)正反演理论框架，进行三维有限体积、三维标量有限元或三维矢量有限元方法三维反演计算得到异常体的非结构化四面体等效源模型，对该模型进行相应的基于PDE的三维有限体积、三维标量有限元或三维矢量有限元正演计算，能对复杂环境中异常数据进行自适应、高精度的延拓与数据类型转换计算。

本发明提出的一种非结构化等效源的重磁数据延拓与转换方法，具体包括以下步骤：

S1、在重力或磁法探测中，通过观测得到物理场数据d₀；所述物理场为重力场或磁场；

S2、给定反演初始模型和反演网格边界，利用现有的四面体网格生成方法生成非结构化四面体网格，作为等效源反演网格空间；

S3、根据地磁场参数，对物理场数据d₀进行带规整化项的PDE三维反演计算，在反演迭代的同时对反演等效源模型周围的反演网格空间，进行加密剖分动态调整，在满足迭代条件后得到异常体的最终等效源模型；

S4、根据实验需要设定不同的正演参数，利用步骤S3得到的等效源模型进行相应的三维正演计算，得到等效源产生的上延拓数据或下延拓数据、化极数据、异常总场数据、异常分量数据以及梯度张量数据。

进一步地，步骤S3中，所述反演计算采用三维有限体积、三维标量有限元或三维矢量有限元方法。

进一步地，步骤S3中，反演计算的目标函数为：

其中，

表示反演目标函数，采用非线性优化求解；d＝F(B₀,m),m≥0，表示正演计算得到的与d₀相应物理场数据；F(·)表示对等效源模型的PDE三维正演计算；m表示等效源模型；φ_reg表示规整化函数；B₀表示背景场，即地磁场，由所在地理位置确定。

进一步地，其特征在于：步骤S4中，异常总场数据、异常分量数据以及梯度张量数据，根据物理场不同而不同，具体为：

若物理场为磁场，则异常总场数据、异常分量数据以及梯度张量数据分别为磁总场异常数据、磁三分量异常数据和磁梯度张量数据，如下：

所述磁三分量异常数据，表示为式(2)：

式(2)中，B_s为磁三分量异常数据的向量形式；

为磁三分量异常数据的三个分量；

表示地磁场的三个分量；φ表示重力势或磁势；μ＝μ₀(1+χ)，μ₀表示真空磁导率，χ表示非结构化网格模型的磁化率；

所述磁总场异常数据，表示为式(3)：

所述磁梯度张量数据表示为式(4)：

若物理场为重力场，则异常总场数据、异常分量数据以及梯度张量数据分别为重力异常数据、重力三分量异常数据和重力梯度张量数据，如下：

所述重力三分量异常数据，表示为式(5)：

其中，g^x、g^y和g^z分别表示重力异常的三个分量的数据；

所述重力异常数据表示为式(6)：

所述重力梯度张量数据表示为式(7)：

进一步地，步骤S4中，利用步骤S3得到的等效源模型进行相应的三维正演计算，对于重力场，正演公式为：

对于磁场，正演公式为：

其中，γ表示万有引力常数，ρ表示非结构化网格模型的密度；

进一步地，步骤S4中，根据原始观测位置空间坐标p(x₀,y₀,z₀)，以及实验需要设定延拓后的观测面位置空间坐标p_up(x,y,z)/p_down(x,y,z)，进一步得到上延拓和/或下延拓后的物理场数据为：B′＝CB_s，其中，B′表示上延拓或下延拓后的数据，C表示根据延拓后观测面的位置与原始观测位置得到的插值函数。

步骤S4中根据实验需要进行磁场数据化极计算时，利用步骤S3中确定的等效源模型m，给定B′₀，进行PDE正演计算，得到化极数据B′_s＝F(B′₀,m)。

本发明提供的有益效果是：本发明提出的技术方案基于PDE框架进行三维有限体积、三维标量有限元或三维矢量有限元方法，能对异常体产生的异常数据进行上下延拓运算和数据类型转换，对于计算空间中物理场分布具有良好的适用性；本发明利用现有的四面体网格生成方法生成非结构化四面体网格，用非结构化四面体网格描述任意模型几何形状具有更高的灵活性，能对异常体进行自适应且快速高效准确地生成所需要的上延拓或下延拓数据和完成数据类型转换，具有更高的稳定性和精度。

附图说明

图1是一种非结构化等效源的重磁数据延拓与转换方法的流程图；

图2是本发明实施例提供的反演空间示意图；

图3是本发明实施例提供的上/下延拓以及化极转换示意图。

具体实施方式

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

请参考图1，本发明的实施例提供了一种基于非结构化等效源的重磁数据延拓与转换方法，以非结构化四面体等效源磁场数据延拓与转换方法为例包括以下步骤包括以下步骤：

S1、输入已有的物理场数据(重力场或磁场)d₀；

S2、给定反演初始模型和反演网格边界，利用现有的四面体网格生成方法生成非结构化四面体网格，作为等效源反演网格空间；

S3、根据地磁场参数，对磁场数据d₀进行带规整化项的PDE三维反演计算，所述反演计算可使用三维有限体积、三维标量有限元或三维矢量有限元方法，在反演迭代的同时对反演等效源模型周围的反演网格空间，进行加密剖分动态调整，在满足迭代条件后得到异常体的最终等效源模型；

S4、根据实验需要设定不同的正演参数，利用步骤S3得到的等效源模型进行相应的三维正演计算，得到等效源产生的上延拓和/或下延拓数据、化极数据、异常总场数据、异常分量数据以及梯度张量数据中的一种或多种。

基于PDE的三维反演计算的目标函数为：

其中，

表示反演目标函数，采用非线性优化求解；d＝F(B₀,m),m≥0，表示正演计算得到的与d₀相应物理场数据；F(·)表示对等效源模型的PDE三维正演计算；m表示等效源模型；φ_reg表示规整化函数，包括深度规整化函数、最小体积、最紧支撑函数等；B₀表示背景场，即地磁场，由所在地理位置确定。同时，在反演迭代过程中根据每次迭代结果，对当前等效源模型周围的反演网格空间进行加密剖分动态调整。如图2所示是本发明实施例提供的反演空间示意图。

本实施例中，所述正演计算采用有限体积的PDE方法，需要说明的是，为满足有限体积求解条件，需对反演网格空间进行扩展，对设定的反演最大深度以下、起伏观测曲面以上、以及反演网格空间周围的水平空间，根据有限体积法进行网格扩充。

在根据实验需要进行数据延拓与转换时，利用步骤S3得到的反演最终模型进行相应的三维正演计算，得到异常体产生的上延拓或下延拓数据(也可以得到上延拓和下延拓数据，这里不作限定)、化极数据、异常总场数据、异常分量数据以及梯度张量数据中的一种或多种。磁异常体正演过程的PDE形式为：

其中，μ＝μ₀(1+χ)，μ₀表示真空磁导率，χ表示非结构化网格模型的磁化率，φ表示重力势或磁势。

在基于非结构化四面体等效源的正演过程中，每个四面体单元具有均匀(相同)的磁化率，通过有限体积或有限元方法，采用非结构化四面体二阶插值方法，求解PDE得到所有四面体单元节点处的磁势φ。通过每个四面体单元10个节点的φ值，结合有限元非结构化四面体二阶插值函数，计算得到该四面体内任意观测点的磁三分量异常数据。接下来可根据三分量异常数据进一步计算得到其它类型的磁正反演数据。

对于三维磁正演，三分量异常数据表示为

其中，μ＝μ₀(1+χ)，μ₀表示真空磁导率，χ表示异常体非结构化网格模型的磁化率，B₀表示地磁场(背景场)，

和

分别为地磁场三个分量的数据，φ表示磁势，

和

分别表示磁三分量异常数据。

如图3所示，为该实例的上延拓、下延拓以及化极转换的示意图。

根据步骤S3中确定的反演最终模型矩阵m以及设定的化极正演参数B′₀(B′₀的磁倾角为90°，磁偏角为0°)，进行PDE正演计算，得到化极数据B′_s：

步骤S4中，根据实验需要所转换不同类型的数据的计算过程为：

磁总场异常数据表示为

磁梯度张量数据表示为

其中，

表示梯度算子，x、y、z分别表示北向、东向及垂向的位置量。

同样，对于重力场，正演公式为：

γ表示万有引力常数，ρ表示非结构化网格模型的密度；

对应的，重力三分量异常数据，表示为式：

其中，g^x、g^y和g^z分别表示重力异常的三个分量的数据；

所述重力异常数据表示为式：

所述重力梯度张量数据表示为式(7)：

本发明提供的有益效果是：本发明提出的技术方案基于PDE框架进行三维有限体积、三维标量有限元或三维矢量有限元方法的正反演计算，对于计算空间中物理场分布具有良好的适用性；本发明利用现有的四面体网格生成方法生成非结构化网格，用非结构化四面体网格描述任意模型几何形状具有更高的灵活性，具有更高的稳定性和精度。

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

Claims (7)

1.一种非结构化等效源的重磁数据延拓与转换方法，其特征在于：具体包括以下步骤：

S1、在重力或磁法探测中，通过观测得到物理场数据d₀；所述物理场为重力场或磁场；

S2、给定反演初始模型和反演网格边界，利用现有的四面体网格生成方法生成非结构化四面体网格，作为等效源反演网格空间；

S3、根据地磁场参数，对物理场数据d₀进行带规整化项的PDE三维反演计算，在反演迭代的同时对反演等效源模型周围的反演网格空间，进行加密剖分动态调整，在满足迭代条件后得到异常体的最终等效源模型；

S4、根据实验需要设定不同的正演参数，利用步骤S3得到的等效源模型进行相应的三维正演计算，得到等效源产生的上延拓数据或下延拓数据、化极数据、异常总场数据、异常分量数据以及梯度张量数据。

2.如权利要求1所述的一种非结构化等效源的重磁数据延拓与转换方法，其特征在于：步骤S3中，所述反演计算采用三维有限体积、三维标量有限元或三维矢量有限元方法。

3.如权利要求1所述的一种非结构化等效源的重磁数据延拓与转换方法，其特征在于：步骤S3中，反演计算的目标函数为：

其中，

表示反演目标函数，采用非线性优化求解；d＝F(B₀,m),m≥0，表示正演计算得到的与d₀相应物理场数据；F(·)表示对等效源模型的PDE三维正演计算；m表示等效源模型；φ_reg表示规整化函数；B₀表示背景场，即地磁场，由所在地理位置确定。

4.如权利要求2所述的一种非结构化等效源的重磁数据延拓与转换方法，其特征在于：步骤S4中，异常总场数据、异常分量数据以及梯度张量数据，根据物理场不同而不同，具体为：

若物理场为磁场，则异常总场数据、异常分量数据以及梯度张量数据分别为磁总场异常数据、磁三分量异常数据和磁梯度张量数据，如下：

所述磁三分量异常数据，表示为式(2)：

式(2)中，B_s为磁三分量异常数据的向量形式；

为磁三分量异常数据的三个分量；

表示地磁场的三个分量；φ表示重力势或磁势；μ＝μ₀(1+χ)，μ₀表示真空磁导率，χ表示非结构化网格模型的磁化率；

所述磁总场异常数据，表示为式(3)：

所述磁梯度张量数据表示为式(4)：

若物理场为重力场，则异常总场数据、异常分量数据以及梯度张量数据分别为重力异常数据、重力三分量异常数据和重力梯度张量数据，如下：

所述重力三分量异常数据，表示为式(5)：

其中，g^x、g^y和g^z分别表示重力异常的三个分量的数据；

所述重力异常数据表示为式(6)：

所述重力梯度张量数据表示为式(7)：

5.如权利要求4所述的一种非结构化等效源的重磁数据延拓与转换方法，其特征在于：步骤S4中，利用步骤S3得到的等效源模型进行相应的三维正演计算，对于重力场，正演公式为：

对于磁场，正演公式为：

其中，γ表示万有引力常数，ρ表示非结构化网格模型的密度。

6.如权利要求4所述的一种非结构化等效源的重磁数据延拓与转换方法，其特征在于：步骤S4中，根据原始观测位置空间坐标p(x₀,y₀,z₀)，以及实验需要设定延拓后的观测面位置空间坐标p_up(x,y,z)/p_down(x,y,z)，进一步得到上延拓和/或下延拓后的物理场数据为：B′＝CB_s，其中，B′表示上延拓或下延拓后的数据，C表示根据延拓后观测面的位置与原始观测位置得到的插值函数。

7.如权利要求4所述的一种非结构化等效源的重磁数据延拓与转换方法，其特征在于：步骤S4中根据实验需要进行磁场数据化极计算时，利用步骤S3中确定的等效源模型m，给定B′₀，进行PDE正演计算，得到化极数据B′_s＝F(B′₀,m)。

Images