CN112363236B - 一种基于pde的重力场数据等效源延拓与数据类型转换方法 - Google Patents

Abstract

本发明提供一种基于PDE的重力场数据等效源延拓与数据类型转换方法，包括：输入重力数据，建立地形起伏曲面；根据起伏观测曲面的高程信息和预设反演最大深度，确定网格剖分的空间范围，并进行连续的结构化非均匀网格剖分，确定等效源反演网格空间；根据重力背景场，基于等效源反演网格空间对重力场数据进行带深度规整化因子、正值约束项以及规整化项的PDE三维反演计算，得到多层等效源模型；给定不同的观测参数，利用模型进行非线性PDE的重力场三维正演计算，得到重力异常铅垂数据、重力异常场数据以及重力梯度张量数据中的一种或多种。本发明的有益效果：能对复杂环境中的地下重力异常数据进行自适应、高精度的延拓与数据类型转换计算。

Description

一种基于PDE的重力场数据等效源延拓与数据类型转换方法

技术领域

本发明涉及地球物理勘测技术领域，尤其涉及一种基于PDE的重力场数据等效源延拓与数据类型转换方法。

背景技术

在重力探测中，人们通常对重力场的铅垂一次导数进行探测，但在实际的数据解释中，又往往需要将这些数据转换成所需要的参数和类型，比如不同观测高度(延拓)、梯度张量数据的换算等，重力场数据类型转换的主要挑战在于航测的起伏测量曲面以及数据的不规则测量位置(非网格化测量)，传统的数据类型转换以及延拓计算难以处理这些数据，而采用传统单层、双层、三层等效源的数据转换方法，处理结果精度较低，无法满足应用需求。

现有技术主要包括以下几种：(1)等效源方法，利用虚拟场源模拟实测异常，可用于位场数据的空间延拓(包括曲面延拓)、梯度计算以及分量转换等；选择单层等效源且将其布置于近地表是等效源方法的主要特点。

(2)使用了单层等效源进行重力场数据的分析，使用单层等效源对重力梯度张量数据进行了去噪；为了在保证计算效率的同时高精度重构位场，多层等效源方法是一个相对合理的选择，

(3)文献“李端；陈超；杜劲松；梁青；黎海龙；基于多层等效源的重力数据处理方法[C]//2017中国地球科学联合学术年会.0.”提出了一种将地下等效源分为浅、中、深三层来实现重力场数据的转换，

(4)文献“陈涛,张贵宾.利用重力异常计算重力梯度的等效源技术[J].地球物理学进展.”提出了一种多层等效源进行重力梯度计算的方法，其在垂向方向上将剖分区域分为三个子区域。

这些技术的应用使得人们可以通过等效源技术来进行地面重力数据的化极与数据类型转换运算，但存在以下问题：1)等效源网格的构建最大只有三层，且网格剖分不连续，现有技术由于等效源的设置层数有限、地下空间的网格剖分不合理，影响了重力场延拓与数据类型转换的精度及速度；2)需单独对每一个等效层的深度位置进行估算，然后单独放置；3)采用线性正演方法和线性模型反演方法计算等效源，抗干扰能力和计算精度较差。

发明内容

有鉴于此，本发明实际要解决的技术问题是：如何对非规则测量位置d的重力数据进行等效源延拓与数据类型转换，并提高计算精度；

为此，本发明提供了一种基于PDE的重力数据等效源延拓与数据类型转换方法，采用连续的结构化非均匀网格剖分，并基于非线性PDE(Partial Differential Equations，偏微分方程，还可包括有限元、有限体积、有限差分方法)正反演理论框架，引入深度规整化因子，在反演过程中直接确定等效源的深度和分布，可一次性计算整个三维自由空间的重力场数据。

本发明提供一种基于PDE的重力数据等效源延拓与数据类型转换方法，包括以下步骤：

S1、获取起伏观测曲面上的重力场数据d₀，并根据所述重力场数据所在区域的地形高度信息，建立地形起伏曲面；所述起伏观测曲面为观测时已记录好的已知值；

S2、根据所述起伏观测曲面的高程信息以及预先设定的反演最大深度，确定网格剖分的空间范围，并根据地形起伏曲面，对所述空间范围进行连续的结构化非均匀网格剖分，确定等效源反演网格空间；

S3、根据重力背景场U₀，基于等效源反演网格空间对重力场数据d₀进行带深度规整化因子、正值约束项以及规整化项的PDE三维反演计算，得到重力异常体的多层等效源模型；

S4、给定不同的观测参数，利用所述多层等效源模型进行非线性PDE的重力场三维正演计算，得到重力异常体在不同观测位置进行延拓后产生的重力异常铅垂数据、重力异常场数据U_s′以及重力梯度张量数据。

本发明提供的技术方案带来的有益效果是：本发明提出的技术方案能对起伏观测曲面上非规则测量位置处所获得的地下重力异常体产生的重力异常数据，进行高精度化极与数据类型转换计算；本发明采用连续网格剖分，等效源层数通常在3层以上，计算精度更高，同时加入深度规整化因子，不需要单独估计等效层的深度和范围，能对重力异常体进行自适应且快速高效准确地生成所需要的类型转换后的数据与延拓数据，具有更高的稳定性和精度。

附图说明

图1是本发明实施例提供的基于PDE的重力数据等效源延拓与数据类型转换方法的流程图；

图2是本发明实施例提供的非均匀网格剖分的效果示意图。

具体实施方式

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

请参考图1，本发明的实施例提供了一种基于PDE的重力数据等效源延拓与数据类型转换方法，包括以下步骤：

S1、获取起伏观测曲面上的重力场数据d₀，并根据重力场数据所在区域的地形高度信息，建立地形起伏曲面；所述重力场数据d₀可以是重力异常铅垂数据或重力梯度张量数据，本实施例以重力异常铅垂数据为例。

S2、根据航测的起伏观测曲面的高程信息以及设定的反演最大深度，确定网格剖分的空间范围，并根据地形起伏曲面，对所述空间范围进行连续的结构化非均匀网格剖分，进一步确定等效源反演网格空间。

具体地，网格剖分的空间范围包括上顶面和下底面，其中，根据起伏观测曲面的最大高度确定上顶面，然后基于现有探测技术或实际经验估计出重力异常体可能存在的最大深度并以此设定反演最大深度，从而确定下底面；

在确定网格剖分的空间范围之后，以地形起伏曲面的最低点对所述空间范围进行划分，对最低点以上的空间范围进行均匀网格剖分得到精细网格，对最低点以下的空间范围进行非均匀的扩展网格剖分得到扩展网格；优选地，若所述精细网格的垂直边为1长度单位(该长度单位的具体数值可根据观测区域空间大小进行设定，比如100m为1长度单位)，则所述扩展网格的垂直边以精细网格垂直边的1.2倍速增长，且设定最大增速为1.5倍，由此在保证一定反演精度的基础上降低计算量。请参考图2，其为本实施进行非均匀网格剖分的结果示意图，其中地形起伏曲面至下底面之间的空间范围，构成等效源反演网格空间。

S3、根据重力背景场U₀，基于等效源反演网格空间对重力场数据d₀进行带深度规整化因子、正值约束项以及规整化项的PDE三维反演计算，得到重力异常体的多层等效源模型，所述多层等效源模型具体是指反演模型求解空间中包含多个模型深度面，基于上述方案，可求解得到的模型深度面数目通常在3层以上。

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

其中，

U_s＝F(U₀,m)

m≥0

式中，φ表示优化目标(即误差)，

表示目标函数的数值约束，

表示目标函数的模型约束，m表示待求解的多层等效源模型的密度矩阵，考虑到物体的物理性质导致密度必须为正值，故进行正值约束，即m≥0；F(·)表示多层等效源模型的三维PDE正演计算，可采用有限体积或有限元方法，U_s表示正演操作得到的重力异常场数据，T(·)表示重力异常场数据到重力异常铅垂数据的转换函数；Q_x、Q_y、Q_z分别表示北向、东向及垂向上的插值函数，所述插值函数包含观测位置信息，可选用克里金插值函数等；U₀表示重力背景场强度矢量，U_0x、U_0y、U_0z分别表示其北向、东向及垂向分量；β表示根据实际需求添加的规整化因子，若不需要可令β＝1；m_ref表示参考模型的密度矩阵，W_r表示深度规整化因子：

其中，z表示等效源到地形起伏曲面的距离，z₀表示起伏观测曲面高度，r表示深度系数，一般取3。

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

正演计算中，重力场PDE方程为：

式中，U＝U₀+U_s，γ表示引力常量，ρ表示多层等效源模型的密度；基于上述方程，得到重力异常场数据U_s的正演计算公式：

U_s＝F(U₀,m)；

迭代求解目标函数，即最小化误差φ，其中，每次迭代完成后得到新的密度矩阵m(即密度ρ的分布)用以拟合测量得到的重力场数据d₀，最终得到多层等效源模型的密度矩阵m。所述正演计算得到三分量数据，在反演过程中可直接转化为铅垂数据或梯度张量数据，以满足不同观测数据类型的需要。

S4、给定不同的观测参数，利用多层等效源模型进行非线性PDE的重力场三维正演计算，可以得到重力异常体在不同观测位置进行延拓后产生的重力异常铅垂数据、重力异常场数据U_s′以及重力梯度张量数据中的一种或多种。

根据步骤S3中确定的密度矩阵m，进行正演计算，得到延拓后的重力场数据U_s′：

U_s′＝F(U₀′,m)；

重力场张量转换的计算公式为：

其中，

表示梯度算子，U＝[U_x,U_y,U_z]^T表示进行重力场张量转换的重力场数据，等式最右边矩阵中的因子为该重力场数据转换得到的不同张量。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (7)

1.一种基于PDE的重力数据等效源延拓与数据类型转换方法，其特征在于，包括以下步骤：

S1、获取起伏观测曲面上的重力场数据d₀，并根据所述重力场数据所在区域的地形高度信息，建立地形起伏曲面；所述起伏观测曲面为观测时已记录好的已知值；

S2、根据所述起伏观测曲面的高程信息以及预先设定的反演最大深度，确定网格剖分的空间范围，并根据地形起伏曲面，对所述空间范围进行连续的结构化非均匀网格剖分，确定等效源反演网格空间；

S3、根据重力背景场U₀，基于等效源反演网格空间对重力场数据d₀进行带深度规整化因子、正值约束项以及规整化项的PDE三维反演计算，得到重力异常体的多层等效源模型；

S4、给定不同的观测参数，利用所述多层等效源模型进行非线性PDE的重力场三维正演计算，得到重力异常体在不同观测位置进行延拓后产生的重力异常铅垂数据、重力异常场数据U_s′以及重力梯度张量数据；

所述多层等效源模型的模型深度面的层数大于3层；

步骤S3中，所述PDE三维反演计算的目标函数如式(1)所示：

式(1)中，

U_s＝F(U₀,m)

m≥0

φ表示优化目标，即误差；

表示目标函数的数值约束；

表示目标函数的模型约束，m表示待迭代优化求解的多层等效源模型的密度矩阵；m_ref表示参考模型的密度矩阵；F(·)表示多层等效源模型的三维PDE正演计算，U_s表示正演操作得到的重力异常场数据，T(·)表示重力异常场数据到重力异常铅垂数据的转换函数；Q_x、Q_y、Q_z分别表示北向、东向及垂向上的插值函数，所述插值函数包含观测位置信息；U₀表示重力背景场强度矢量，U_0x、U_0y、U_0z分别表示其北向、东向及垂向分量；β表示预设的权重参数；

所述深度规整化因子如式(2)所示：

式(2)中，z表示等效源到地形起伏曲面的距离，z₀表示起伏观测曲面高度，r表示深度系数。

2.根据权利要求1所述的基于PDE的重力数据等效源延拓与数据类型转换方法，其特征在于，步骤S2中，所述网格剖分的空间范围包括上顶面和下底面，其中，所述上顶面为起伏观测曲面的最大高度所确定的平面，所述下底面为预先设定的反演最大深度所确定的平面。

3.根据权利要求2所述的基于PDE的重力数据等效源延拓与数据类型转换方法，其特征在于，步骤S2中，根据地形起伏曲面，对所述空间范围进行连续的结构化非均匀网格剖分，确定等效源反演网格空间，具体为：

根据地形起伏曲面的最低点对所述空间范围进行划分；

对所述最低点以上的空间范围进行均匀网格剖分得到精细网格；

对所述最低点以下的空间范围进行非均匀网格剖分得到扩展网格；

所述地形起伏曲面至下底面之间的空间范围，即所述精细网格和所述扩展网格构成等效源反演网格空间。

4.根据权利要求3所述的基于PDE的重力数据等效源延拓与数据类型转换方法，其特征在于，所述扩展网格的垂直边以精细网格的垂直边的α₁倍速度增长，且设定其最大增速为α₂，其中，α₂>α₁>1。

5.根据权利要求1所述的基于PDE的重力数据等效源延拓与数据类型转换方法，其特征在于，步骤S3中，根据优化时迭代产生的密度矩阵m，与重力背景场U₀进行PDE正演计算，得到重力异常体产生的重力异常场数据：

U_s＝F(U₀,m)。

6.根据权利要求1所述的基于PDE的重力数据等效源延拓与数据类型转换方法，其特征在于，步骤S4中，根据步骤S3中确定的密度矩阵m，进行正演计算，得到延拓后的重力异常场数据U_s′：

U_s′＝F(U₀′,m)。

7.根据权利要求1所述的基于PDE的重力数据等效源延拓与数据类型转换方法，其特征在于，重力场张量数据的计算如式(3)：

其中，

表示梯度算子，U＝[U_x,U_y,U_z]^T表示进行重力场张量转换的重力场数据，式(3)最右边矩阵中的因子为该重力场数据转换得到的不同张量数据。

Images