一种重力数据多层等效源延拓与数据转换方法转让专利
申请号 : CN202011104563.7
文献号 : CN112346139B
文献日 : 2021-12-21
发明人 : 左博新
申请人 : 中国地质大学(武汉)
摘要 :
本发明提供一种重力数据等效源延拓与数据转换方法，包括：输入起伏观测曲面上的重力场数据d0，根据观测面所在区域的地形信息，建立地形起伏曲面；根据起伏观测曲面以及设定的反演最大深度，确定网格剖分范围，并对所述范围进行网格剖分以确定模型空间；根据重力背景场，基于模型空间对重力场数据d0进行带约束项及规整化项的积分方程三维反演计算，得到重力异常体的等效源模型；根据延拓面的位置确定灵敏度矩阵G′，然后利用等效源模型进行基于积分方程的重力场正演计算，得到异常体产生的延拓后的重力场数据Us；对Us进行梯度张量转换。本发明能对重力异常体进行自适应且快速高效准确地生成所需要延拓与数据转换后的数据。
权利要求 :
1.一种重力数据多层等效源延拓与数据转换方法，其特征在于，包括以下步骤：S1、输入起伏观测曲面上的重力场数据d0，并根据起伏观测曲面所在区域的地形高度信息，建立地形起伏曲面；
S2、根据起伏观测曲面以及设定的反演最大深度，确定网格剖分的空间范围，并根据S1建立的地形起伏曲面对所述空间范围进行连续的结构化非均匀网格剖分，进一步确定反演模型求解空间；根据地形起伏曲面的最低点对网格剖分的空间范围进行划分，其中，对所述最低点以上的空间范围进行均匀网格剖分得到精细网格，对所述最低点以下的空间范围进行非均匀网格剖分得到扩展网格；所述地形起伏曲面至下底面之间的空间范围构成反演模型求解空间；
S3、根据重力场的背景场，基于S2反演模型求解空间对重力场数据d0进行带深度规整化因子、正值约束项以及规整化项的积分方程三维反演计算，得到重力异常体的多层等效源模型；所述多层等效源模型的模型深度面的数目大于3层；
积分方程三维反演计算的目标函数为：其中，φ表示优化目标，即误差，m表示待求解的多层等效源模型的密度矩阵，m≥0；F(·)表示基于积分方程的重力场正演操作，G表示灵敏度矩阵，与起伏观测曲面的位置以及网格剖分相关；β表示根据实际需求添加的规整化因子；mref表示参考模型的密度矩阵，Wr表示深度规整化因子；所述深度规整化因子Wr为：其中，z表示等效源到地形起伏曲面的距离，z0表示起伏观测曲面的高度，r表示深度系数；
S4、设定延拓面，并根据延拓面得到灵敏度矩阵G′，然后利用S3得到的多层等效源模型进行基于积分方程的重力场正演计算，得到重力异常体产生的延拓后的重力场数据Us；重力场数据Us＝F(m,G′)＝G′m，其中，m为密度矩阵，F(·)表示基于积分方程的重力场正演操作，G′为由延拓面确定的灵敏度矩阵；
S5、基于S4所求的重力场数据Us，利用重力场梯度张量转换公式，进行重力场数据的转换；重力场数据的转换如下：
其中， 为梯度算子，矩阵[*]中的每项因子均为重力场的不同张量，Ux、Uy以及Uz均为Us的分量数据。
2.根据权利要求1所述的一种重力数据多层等效源延拓与数据转换方法，其特征在于，S2中，所述网格剖分的空间范围包括上顶面和下底面，其中，所述上顶面为地形起伏曲面的最大高度所确定的平面，所述下底面为设定的反演最大深度所确定的平面。
说明书 :
一种重力数据多层等效源延拓与数据转换方法
技术领域
[0001] 本发明涉及重力勘探技术领域，尤其涉及一种重力数据多层等效源延拓与数据转换方法。
背景技术
[0002] 在重力探测中，人们通常对重力场的铅垂一次导数进行探测，但在实际的数据解释中，又往往需要将这些数据转换成所需要的参数和类型，比如不同观测高度(延拓)、梯度
张量数据的换算等，重力场数据类型转换的主要挑战在于航测的起伏测量曲面以及数据的
不规则测量位置(非网格化测量)，传统的数据类型转换以及延拓计算难以处理这些数据，
而采用传统单层、双层、三层等效源的数据转换方法，处理结果精度较低，无法满足应用需
求。
张量数据的换算等，重力场数据类型转换的主要挑战在于航测的起伏测量曲面以及数据的
不规则测量位置(非网格化测量)，传统的数据类型转换以及延拓计算难以处理这些数据，
而采用传统单层、双层、三层等效源的数据转换方法，处理结果精度较低，无法满足应用需
求。
[0003] 现有文献1“Dampney,C.N.G.THE EQUIVALENT SOURCE TECHNIQUE[J].geophysics,1969,34(1):39.”首次提出了等效源方法，利用虚拟场源模拟实测异常，可用
于位场数据的空间延拓(包括曲面延拓)、梯度计算以及分量转换等；选择单层等效源且将
其布置于近地表是等效源方法的主要特点，比如,文献2“Lyrio J C S D
O.2011.Equivalent source:A natural choice for gridding scatter gravity data
[J].In 12th International Congress of the Brazilian Geophysical Society，
pp.661‑665，doi:10.1190/sbgf2011‑136.”使用了单层等效源进行重力场数据的分析，文
献3“Martinez C，Li Y G.2016.Denoising of gravity gradient data using an
equivalent source technique[J].Geophysics，81(4):G67‑G79，doi:10.1190/geo2015‑
0379.1.”使用单层等效源对重力梯度张量数据进行了去噪；为了在保证计算效率的同时高
精度重构位场，多层等效源方法是一个相对合理的选择，文献4“李端；陈超；杜劲松；梁青；
黎海龙；基于多层等效源的重力数据处理方法[C]//2017中国地球科学联合学术年会.0.”
提出了一种将地下等效源分为浅、中、深三层来实现重力场数据的转换，文献5“陈涛,张贵
宾.利用重力异常计算重力梯度的等效源技术[J].地球物理学进展.”提出了一种多层等效
源进行重力梯度计算的方法，其在垂向方向上将剖分区域分为三个子区域。
于位场数据的空间延拓(包括曲面延拓)、梯度计算以及分量转换等；选择单层等效源且将
其布置于近地表是等效源方法的主要特点，比如,文献2“Lyrio J C S D
O.2011.Equivalent source:A natural choice for gridding scatter gravity data
[J].In 12th International Congress of the Brazilian Geophysical Society，
pp.661‑665，doi:10.1190/sbgf2011‑136.”使用了单层等效源进行重力场数据的分析，文
献3“Martinez C，Li Y G.2016.Denoising of gravity gradient data using an
equivalent source technique[J].Geophysics，81(4):G67‑G79，doi:10.1190/geo2015‑
0379.1.”使用单层等效源对重力梯度张量数据进行了去噪；为了在保证计算效率的同时高
精度重构位场，多层等效源方法是一个相对合理的选择，文献4“李端；陈超；杜劲松；梁青；
黎海龙；基于多层等效源的重力数据处理方法[C]//2017中国地球科学联合学术年会.0.”
提出了一种将地下等效源分为浅、中、深三层来实现重力场数据的转换，文献5“陈涛,张贵
宾.利用重力异常计算重力梯度的等效源技术[J].地球物理学进展.”提出了一种多层等效
源进行重力梯度计算的方法，其在垂向方向上将剖分区域分为三个子区域。
[0004] 这些技术的应用使得人们可以通过等效源技术来进行地面重力数据的化极与数据类型转换运算，但存在以下问题：1)上述方法构建的等效源层数通常小于等于三层，且网
格剖分不连续；2)需单独对每一个等效源的深度位置进行估算，然后单独放置；3)传统延拓
方法对于数据的延拓距离有限制(一般小于6倍的数据点距)，大于限制距离的延拓计算精
度较差。
格剖分不连续；2)需单独对每一个等效源的深度位置进行估算，然后单独放置；3)传统延拓
方法对于数据的延拓距离有限制(一般小于6倍的数据点距)，大于限制距离的延拓计算精
度较差。
发明内容
[0005] 有鉴于此，本发明提供了一种重力数据多层等效源延拓与数据转换方法，采用连续的结构化非均匀网格剖分，并基于积分方程正反演理论框架，引入深度规整化因子，在反
演过程中直接确定多层等效源的深度和分布，包括以下步骤：
演过程中直接确定多层等效源的深度和分布，包括以下步骤：
[0006] S1、输入起伏观测曲面上的重力场数据d0，并根据起伏观测曲面所在区域的地形高度信息，建立地形起伏曲面；
[0007] S2、根据起伏观测曲面以及设定的反演最大深度，确定网格剖分的空间范围，并根据S1建立的地形起伏曲面对所述空间范围进行连续的结构化非均匀网格剖分，进一步确定
反演模型求解空间；
反演模型求解空间；
[0008] S3、根据重力场的背景场，基于S2反演模型求解空间对重力场数据d0进行带深度规整化因子、正值约束项以及规整化项的积分方程三维反演计算，得到重力异常体的多层
等效源模型；
等效源模型；
[0009] S4、设定延拓面，并根据延拓面得到灵敏度矩阵G′，然后利用S3得到的多层等效源模型进行基于积分方程的重力场正演计算，得到重力异常体产生的延拓后的重力场数据
Us；
Us；
[0010] S5、基于S4所求的重力场数据Us，利用重力场梯度张量转换公式，进行重力场数据的转换。
[0011] 进一步地，S2中，所述网格剖分的空间范围包括上顶面和下底面，其中，所述上顶面为地形起伏曲面的最大高度所确定的平面，所述下底面为设定的反演最大深度所确定的
平面。
平面。
[0012] 进一步地，S2中，根据地形起伏曲面的最低点对网格剖分的空间范围进行划分，其中，对所述最低点以上的空间范围进行均匀网格剖分得到精细网格，对所述最低点以下的
空间范围进行非均匀网格剖分得到扩展网格；所述地形起伏曲面至下底面之间的空间范围
构成反演模型求解空间。
空间范围进行非均匀网格剖分得到扩展网格；所述地形起伏曲面至下底面之间的空间范围
构成反演模型求解空间。
[0013] 进一步地，所述扩展网格的垂直边以精细网格的垂直边的α1倍速度增长，且设定其最大增速为α2，其中，α2>α1>1。
[0014] 进一步地，S3中所述多层等效源模型的模型深度面的数目大于3层。
[0015] 进一步地，S3中，积分方程三维反演计算的目标函数为：
[0016]
[0017] 其中，φ表示优化目标即误差，m表示待求解的多层等效源模型的密度矩阵，m≥0；F(·)表示基于积分方程的重力场正演操作，G表示灵敏度矩阵，与起伏观测曲面的位置以
及网格剖分相关；β表示根据实际需求添加的规整化因子；mref表示参考模型的密度矩阵，Wr
表示深度规整化因子。
及网格剖分相关；β表示根据实际需求添加的规整化因子；mref表示参考模型的密度矩阵，Wr
表示深度规整化因子。
[0018] 进一步地，所述深度规整化因子Wr为：
[0019]
[0020] 其中，z表示等效源到地形起伏曲面的距离，z0表示起伏观测曲面的高度，r表示深度系数。
[0021] 进一步地，S4得到重力异常体产生的延拓后的重力场数据Us＝F(m,G′)＝G′m，其中，m为密度矩阵，F(·)表示基于积分方程的重力场正演操作，G′为由延拓面确定的灵敏度
矩阵。
矩阵。
[0022] 进一步地，S5中重力场数据的转换如下：
[0023]
[0024] 其中， 为梯度算子，矩阵[*]中的每项因子均为重力场的不同张量，Ux、Uy以及Uz均为Us的分量数据。
[0025] 本发明提供的技术方案带来的有益效果是：本发明提出的技术方案能对复杂环境中的地下重力异常体产生的重力场进行基于积分方程的非均匀网格的多层等效源延拓与
数据类型转换运算，采用连续网格剖分，数目通常在3层以上，精度更高；同时加入深度规整
化因子，不需要单独估计等效层的深度和范围，能对重力异常体进行自适应且快速高效准
确地生成所需要的延拓与数据类型转换后的数据，具有更高的稳定性和精度。
数据类型转换运算，采用连续网格剖分，数目通常在3层以上，精度更高；同时加入深度规整
化因子，不需要单独估计等效层的深度和范围，能对重力异常体进行自适应且快速高效准
确地生成所需要的延拓与数据类型转换后的数据，具有更高的稳定性和精度。
附图说明
[0026] 图1是本发明一种重力数据多层等效源延拓与数据转换方法的流程图；
[0027] 图2是本发明一种重力数据多层等效源延拓与数据转换方法的非均匀网格剖分的示意图。
具体实施方式
[0028] 为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。
[0029] 请参考图1，本发明的实施例提供了一种重力数据多层等效源延拓与数据转换方法，包括以下步骤：
[0030] S1、输入起伏观测曲面上的重力场数据d0，并根据观测面所在区域的地形高度信息，建立地形起伏曲面；所述重力场数据d0可以是重力异常铅垂数据或重力梯度张量数据，
本实施例以重力异常铅垂数据为例；
本实施例以重力异常铅垂数据为例；
[0031] S2、根据地形起伏曲面以及设定的反演最大深度，确定网格剖分的空间范围，并根据地形起伏曲面对所述空间范围进行连续的结构化非均匀网格剖分，进一步确定反演模型
求解空间；
求解空间；
[0032] 网格剖分的空间范围包括上顶面和下底面，其中，根据起伏观测曲面的最大高度确定上顶面，然后基于现有探测技术或实际经验估计出重力异常体可能存在的最大深度
(即反演最大深度)，并据此确定下底面；
(即反演最大深度)，并据此确定下底面；
[0033] 在确定网格剖分的空间范围之后，以地形起伏曲面的最低点对所述空间范围进行划分，对最低点以上的空间范围进行均匀网格剖分得到精细网格，对最低点以下的空间范
围进行非均匀的扩展网格剖分得到扩展网格，所述地形起伏曲面至下底面之间的空间范围
构成反演模型求解空间；优选地，若所述精细网格的垂直边为1长度单位(该长度单位的具
体数值可根据观测区域空间大小进行设定，比如100m为1长度单位)，则所述扩展网格的垂
直边以精细网格垂直边的1.2倍速增长，且设定最大增速为1.5倍，由此在保证一定反演精
度的基础上降低计算量。请参考图2，其为非均匀网格剖分的结果示意图，图中坐标轴分别
为东向(Easting)、北向(Northing)、垂向(Depth)。
围进行非均匀的扩展网格剖分得到扩展网格，所述地形起伏曲面至下底面之间的空间范围
构成反演模型求解空间；优选地，若所述精细网格的垂直边为1长度单位(该长度单位的具
体数值可根据观测区域空间大小进行设定，比如100m为1长度单位)，则所述扩展网格的垂
直边以精细网格垂直边的1.2倍速增长，且设定最大增速为1.5倍，由此在保证一定反演精
度的基础上降低计算量。请参考图2，其为非均匀网格剖分的结果示意图，图中坐标轴分别
为东向(Easting)、北向(Northing)、垂向(Depth)。
[0034] S3、根据重力场的背景场，基于S2反演模型求解空间对重力场数据d0进行带深度规整化因子、正值约束项以及规整化项的积分方程三维反演计算，得到重力异常体的多层
等效源模型，所述多层等效源模型具体是指反演模型求解空间中包含多个模型深度面，基
于上述方案，可求解得到的模型深度面数目通常大于3层。
等效源模型，所述多层等效源模型具体是指反演模型求解空间中包含多个模型深度面，基
于上述方案，可求解得到的模型深度面数目通常大于3层。
[0035] 基于积分方程的三维反演计算的目标函数为：
[0036]
[0037] 式中，φ表示优化目标，即误差， 表示目标函数的数值约束，表示目标函数的模型约束，m表示待求解的多层等效源模型的密度矩
阵，考虑到物体的物理性质导致密度必须为正值，故进行正值约束，即m≥0；F(·)表示基于
积分方程的重力场正演操作，G表示灵敏度矩阵，与起伏观测曲面的位置以及网格剖分相
关；β表示根据实际需求添加的规整化因子，若不需要可令β＝1；mref表示参考模型的密度矩
阵，Wr表示深度规整化因子，计算公式如下：
阵，考虑到物体的物理性质导致密度必须为正值，故进行正值约束，即m≥0；F(·)表示基于
积分方程的重力场正演操作，G表示灵敏度矩阵，与起伏观测曲面的位置以及网格剖分相
关；β表示根据实际需求添加的规整化因子，若不需要可令β＝1；mref表示参考模型的密度矩
阵，Wr表示深度规整化因子，计算公式如下：
[0038]
[0039] 其中，z表示等效源到地形起伏曲面的距离，z0表示起伏观测曲面的高度，r表示深度系数，本实施例中取r＝3。
[0040] S4、设定延拓面，并根据延拓面得到灵敏度矩阵G′，然后利用步骤S3得到的多层等效源模型进行基于积分方程的重力场正演计算，得到重力异常体产生的延拓后的重力场数
据Us，即：
据Us，即：
[0041] Us＝F(m,G′)＝G′m，
[0042] 其中，m表示S3输出的密度矩阵；
[0043] 通过下述公式从重力场数据Us转换为重力异常数据：
[0044]
[0045] 其中，d为铅垂方向的重力异常数据；
[0046] S5、基于所求的重力场数据Us，利用下述的重力场梯度张量转换公式，进行重力场数据的转换：
[0047]
[0048] 其中， 为梯度算子，矩阵[*]中的每项因子均为重力场的不同张量，Ux、Uy以及Uz均为Us的分量数据。
[0049] 上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员
在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多
形式，这些均属于本发明的保护之内。
在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多
形式，这些均属于本发明的保护之内。