CN112949133B - 一种基于pde的重磁联合反演方法 - Google Patents

Abstract

本发明公开了一种基于PDE的重磁联合反演方法，包括：给定反演初始模型和反演网格边界，利用网格生成方法生成非结构化四面体网格或结构化六面体网格，作为反演网格空间；根据反演初始模型，以及物理场的背景场参数进行PDE重力场与磁场三维正演计算，得到物理场数据d；根据物理场数据d,以及与S2相对应的反演计算，进行反演模型的PDE三维联合反演迭代计算，迭代过程中，根据迭代结果，对当前反演模型周围的反演网格空间进行加密剖分动态调整；当反演目标函数求解误差达到预先设定误差要求时，输出最终反演模型，否则返回S3继续进行三维联合反演迭代计算。本发明对于计算空间中物理场分布、描述几何形状具有良好的适用性，更高的稳定性和精度。

Description

一种基于PDE的重磁联合反演方法

技术领域

本发明涉及地球物理探测技术领域，尤其涉及一种基于PDE的重磁联合反演方法。

背景技术

根据物理场理论，运用定量描述，由已知的异常体性质求出重力场/磁场分布的过程称为正演问题；反之，由重力场/磁场异常求异常体密度、磁性参数、几何参数的过程被称为反演问题。对于反演问题进行求解时，必须建立在正演问题给出的势场理论数学表达式基础上才能进行。

文献“Li Y G,Oldenburg D W.3-D inversion of magnetic data[J].Geophysics,1996,61(2):394-408.”和文献“Li Y,Oldenburg D W.3-D inversion ofgravity data[J].Geophysics,1998,63(1):109-119.”提出了一种通过积分方程对结构化六面体(立方体)网格单元分别进行磁场反演与重力场反演的方法。文献“Li,Yaoguo.3-Dinversion of gravity gradiometer data[J].Seg Technical Program ExpandedAbstracts,2001:2135.”提出了一种通过重力梯度张量数据进行磁场三维反演的方法。“Peter G.Lelièvre,Oldenburg D W.Magnetic forward modelling and inversion forhigh susceptibility[J].Geophysical Journal International,2006.”提出了一种通过PDE(偏微分方程)进行有限体积的结构化立方体网格的正反演方法。文献“Thomas Günther,Carsten Rücker,Klaus Spitzer.Three-dimensional modelling and inversionof dc resistivity data incorporating topography–II.Inversion[J].GeophysicalJournal International,2006,166(2).”提出了基于非结构化网格的电法反演方法，未见基于非结构化网格的三维重磁反演；文献“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方法在基础理论与方程上显然不同。文献“Luis A.Gallardo,M.A.Pérez-Flores,E.Gómez-Refinement ofthree-dimensional multilayer models of basins and crustal environments byinversion of gravity and magnetic data.2004,397(1):37-54”实现了重磁数据联合反演。文献“Gallardo L A.Multiple cross-gradient joint inversion for geospectralimaging[J].Geophysical Research Letters,2007,34(19):5.”实现了电法跟磁法的交叉梯度联合反演，未见采用PDE方法的结构化和非结构的交叉梯度联合反演。

然而这些方法存在以下问题：1)结构化六面体(立方体)网格对联合反演模型空间和异常体的剖分精度与效率存在条件限制：一方面，采用立方体单元在对复杂三维结构进行剖分时，容易导致模型离散化误差过大，拟合精度低；另一方面，若采用大量精细立方体单元对复杂三维结构进行剖分，则会导致三维正反演计算量过大，对硬件性能及计算时间要求较高；2)积分方程方法所构建的计算矩阵为非稀疏矩阵，因此在联合反演中求解复杂度高，难以实现大规模并行化计算。

发明内容

针对以上技术问题，本发明提供了一种基于PDE的重磁联合反演方法，利用现有的网格生成方法构建非结构化四面体网格或结构化六面体网格，并基于非线性PDE(PartialDifferential Equations，偏微分方程)正反演理论框架，进行三维有限体积、三维标量有限元或三维矢量有限元方法的三维联合正反演计算。本方法所涉及的各个计算矩阵全部为稀疏矩阵，具备并行化计算条件，为大规模重磁数据处理分析提供了基础。

本发明提出的一种基于PDE的重磁联合反演方法，具体包括以下步骤：

S1、给定反演初始模型和反演网格边界，利用已知的网格生成方法生成非结构化四面体网格或结构化六面体网格，作为反演网格空间；

S2、根据反演初始模型，以及物理场的背景场参数进行PDE重力场与磁场三维正演计算，得到物理场数据d；所述物理场包括重力场和磁场；

S3、根据物理场数据d,以及与S2中三维正演计算相对应的反演计算，进行反演磁化率模型的PDE三维联合反演迭代计算，迭代过程中，根据迭代结果，对当前反演模型周围的反演网格空间进行加密剖分动态调整；

S4、当反演目标函数求解误差达到预先设定误差要求时，输出最终反演模型，否则返回S3继续进行三维联合反演迭代计算。

进一步地，步骤S2中，具体为：根据反演初始模型，采用三维有限体积、三维标量有限元或三维矢量有限元方法，基于偏微分方程PDE的物理场的正演公式，计算物理场数据d。

步骤S3中，重磁数据的PDE三维联合反演的目标函数如式(1)：

其中Φ_d可以表示为如式(2)：

其中为包含深度权重的权重矩阵；d^obs为实际观测数据，参数β₁和β₂为常数值；

Φ_cg为交叉梯度函数，可以表示为如式(3)：

其中，m₁和m₂表示密度和磁化率两种不同的模型参数，D表示联合反演空间；Φ(m₁,m₂)表示三维联合反演目标函数，采用非线性优化求解；η为动态规整系数，在每一次迭代过程中不断调整η的值，当目标函数最小化达到阈值时，完成非线性优化求解过程。

本发明提供的有益效果是：基于PDE框架进行三维重磁联合正反演计算，利用重力梯度张量数据和磁梯度张量数据两种不同的数据联合对反演目标函数进行约束，具有更高的精度。

附图说明

图1是一种基于PDE的重磁联合反演方法的流程图；

图2是本发明实施例提供的结构化网格反演模型效果示意图；

图3是本发明实施例提供的非结构化网格反演模型效果示意图。

具体实施方式

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

请参阅图1，图1是本发明一种基于PDE的重磁联合反演方法的流程图，具体包括下面步骤：

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

S2、根据反演初始模型，以及物理场的背景场参数进行PDE重力场与磁场三维正演计算，得到物理场数据d；所述物理场为重力场和磁场；对于联合反演中的重力场，则所述重力场数据为重力梯度张量数据；对于联合反演中的磁场，所述磁场数据为磁梯度张量数据；

S3、根据物理场数据d,以及与S2中正演计算方法相对应的反演计算，进行反演磁化率模型的PDE三维联合反演迭代计算，迭代过程中，根据迭代结果，对当前反演模型周围的反演网格空间进行加密剖分动态调整；

S4、当三维联合反演目标函数求解误差达到预先设定误差要求时，输出最终反演模型，否则返回S3继续进行三维联合反演迭代计算；

本实施例中，所述三维联合正演计算采用三维有限体积的PDE方法并以重磁数据为例，需要说明的是，为满足有限体积求解条件，需对正反演网格空间进行扩展，请参考图2、图3，对设定的反演最大深度以下、观测面以上、以及反演网格空间周围的水平空间，根据有限体积法进行网格扩充。

磁场正演过程的PDE形式为：

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

为了保证正演过程的精度，对观测面和地形的网格进行了加密剖分。

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

其中，μ＝μ₀(1+χ)，μ₀表示真空磁导率，χ表示异常体非结构化网格模型的磁化率，B₀表示地磁场(背景场)，和/>分别为地磁场三个分量的数据，φ表示磁势，和/>分别表示磁三分量异常数据。

磁梯度张量数据表示为

对于重力场，正演公式为：γ表示万有引力常数，ρ表示非结构化网格模型的密度；

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

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

所述重力梯度张量数据表示为：

优选地，基于PDE的三维联合反演计算的目标函数为：

其中为包含深度权重的权重矩阵；d^obs为实际观测数据；参数β₁和β₂为常数值；Φ_cg为交叉梯度函数，可以表示为：

其中为包含深度权重的权重矩阵；参数β₁和β₂为常数值；Φ_cg为交叉梯度函数，可以表示为：

其中，m₁和m₂表示密度和磁化率两种不同的模型参数，D表示联合反演空间；Φ(m₁,m₂)表示三维联合反演目标函数，采用非线性优化求解；η为动态规整系数，在每一次迭代过程中不断调整η的值，当目标函数最小化达到阈值时，完成非线性优化求解过程。

通过迭代求解目标函数，即最小化误差，其中，每次迭代完成后得到新的m₁和m₂，即密度矩阵和磁化率矩阵，m₁和m₂用以拟合测量得到的磁场数据和重力场数据，最终得到最优化的异常体结构化或非结构化四面体网格模型的磁化率矩阵和密度矩阵m₁和m₂。

本发明提供的有益效果是：基于PDE框架进行三维重磁联合反演计算，利用重力梯度张量数据和磁梯度张量数据两种不同的数据联合对反演目标函数进行约束，具有更高的精度。

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

Claims (1)

1.一种基于PDE的重磁联合反演方法，其特征在于：具体包括以下步骤：

S1、给定反演初始模型和反演网格边界，利用已知的网格生成方法生成非结构化四面体网格或结构化六面体网格，作为反演网格空间；

S2、根据反演初始模型，以及物理场的背景场参数进行PDE重力场与磁场三维正演计算，得到物理场数据d；所述物理场包括重力场和磁场；

S3、根据物理场数据d,以及与S2中三维正演计算相对应的反演计算，进行反演磁化率模型的PDE三维联合反演迭代计算，迭代过程中，根据迭代结果，对当前反演模型周围的反演网格空间进行加密剖分动态调整；

S4、当三维联合反演目标函数求解误差达到预先设定误差要求时，输出最终反演模型，否则返回S3继续进行三维联合反演迭代计算；

步骤S2中，具体为：根据反演初始模型，采用三维有限体积、三维标量有限元或三维矢量有限元方法，基于偏微分方程PDE的物理场的三维正演公式，计算物理场数据d；其中物理场包括重力场和磁场；

对于联合反演中的重力场，则所述重力场数据为重力梯度张量数据；对于联合反演中的磁场，所述磁场数据为磁梯度张量数据；

步骤S3中，重磁数据的PDE三维联合反演的目标函数如式(1)：

其中Φ_d可以表示为如式(2)：

其中为包含深度权重的权重矩阵；d^obs为实际观测数据，参数β₁和β₂为常数值；Φ_cg为交叉梯度函数，可以表示为如式(3)：

其中，m₁和m₂表示密度和磁化率两种不同的模型参数，D表示联合反演空间；Φ(m₁,m₂)表示三维联合反演目标函数，采用非线性优化求解；η为动态规整系数，在每一次迭代过程中不断调整η的值，当目标函数最小化达到阈值时，完成非线性优化求解过程。