CN110045432B - 基于3d-glq的球坐标系下重力场正演方法及三维反演方法 - Google Patents

基于3d-glq的球坐标系下重力场正演方法及三维反演方法 Download PDF

Info

Publication number
CN110045432B
CN110045432B CN201811481164.5A CN201811481164A CN110045432B CN 110045432 B CN110045432 B CN 110045432B CN 201811481164 A CN201811481164 A CN 201811481164A CN 110045432 B CN110045432 B CN 110045432B
Authority
CN
China
Prior art keywords
coordinate system
gravity
spherical coordinate
gravity field
dimensional
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811481164.5A
Other languages
English (en)
Other versions
CN110045432A (zh
Inventor
赵广东
柳建新
陈波
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN201811481164.5A priority Critical patent/CN110045432B/zh
Publication of CN110045432A publication Critical patent/CN110045432A/zh
Application granted granted Critical
Publication of CN110045432B publication Critical patent/CN110045432B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本申请提供了基于三维高斯勒让德数值积分公式(3D‑GLQ)的球坐标系下三维重力场快速高精度反演算法，核心内容包括以下两点：第一是在球坐标系下三维重力场正演方面，针对传统3D‑GLQ三维重力场正演算法计算效率低的问题，提出核矩阵等效存储策略，然后将核矩阵等效性与已有的自适应剖分策略相结合，在保证正演计算的最大相对误差低于0.1％的前提下，将计算效率提高两个数量级，同时大大减少了内存的占用；第二点是将上述快速高精度正演算法用于重力场三维反演当中，构造了球坐标系下重力场反演目标函数及深度加权函数，使得球坐标系下大规模重力场三维反演的计算效率和可信度大大提高。

Description

基于3D-GLQ的球坐标系下重力场正演方法及三维反演方法
技术领域
本申请涉及地球物理领域球坐标系下的重力场勘探计算，特别地，涉及一种基于3D-GLQ、提出了核矩阵等效存储策略、加入了成熟的自适应剖分技术的球坐标系下重力场正演方法及三维反演方法。
背景技术
重力勘探作为一种古老的地球物理方法，已被广泛运用于资源勘探(矿产、石油、天然气)、水文环境和地球大尺度研究。在勘探尺度，重力场反演常常在直角坐标系下进行，多种快速正演算法(包括FFT法，Gauss-FFT法，多级展开法等)使得反演的计算效率和精度得到了保证。然而当研究区域为跨度较大或者研究整个地球的重力场时，本申请必须考虑到地球曲率的影响，相应的数据处理、正演算法和反演都将在球坐标系下进行。球坐标系下的三维重力场正演方法有高斯勒让德积分法和泰勒级数展开法，这些传统方法计算效率低下，很难适应大规模三维反演的要求。在当今计算机计算性能难以有质的飞跃的背景下，想要大幅度提高计算效率，只能寻求新的正演方法。
另一方面，21世纪以来，卫星重力测量技术的发展不仅革命性地推动了地球重力场的研究，也推动了利用地球重力场信息探测地球内部组构及物质迁移的研究。日益丰富的卫星重力(包括卫星测高)、海洋重力、陆地重力及GPS测量数据使地球重力场的解算精度和分辨率得到极大提高。全球覆盖的高分辨率、高精度重力数据为获取区域及全球壳幔高分辨密度扰动模型创造了条件，也是固体地球科学研究重要的基础数据，同时为全球大尺度大规模三维卫星重力反演提供了新的契机。然而传统的勘探尺度的反演方法、目标函数已不再适用，需要重新构造球坐标系下重力场反演目标函数、深度加权函数。
发明内容
本申请目的在于提供基于3D-GLQ的球坐标系下重力场正演方法及三维反演方法，以解决目前地球的重力场正演方法效率太低，使得反演方法、目标函数已不适用现代科学发展的技术问题。
为实现上述目的，本申请提供了一种基于3D-GLQ的球坐标系下重力场正演方法，包括以下步骤：
A、将地下场源在径向、纬向、经向分别剖分成Nr、
Figure GDA0002492143660000021
Nλ段，记为
Figure GDA0002492143660000022
观测平面位于场源正上方，观测点数位
Figure GDA0002492143660000023
B、场源剖分成的tesseroid单元体是由三对曲面围成：一对同球心的曲面(r1,r2)、一对子午线断面(λ1,λ2)、一对同轴圆锥平面
Figure GDA0002492143660000024
将每一个tesseroid单元体的密度看作是常数，对于场源体外一计算点
Figure GDA0002492143660000025
由第q个tesseroid单元体
Figure GDA0002492143660000026
在该观测点处产生的重力位
Figure GDA0002492143660000027
重力加速度
Figure GDA0002492143660000028
和重力梯度张量
Figure GDA0002492143660000029
可以写为：
Figure GDA00024921436600000210
Figure GDA00024921436600000211
Figure GDA00024921436600000212
其中，
Figure GDA00024921436600000213
是密度分布；G是重力常数；δαβ是克罗内克符号(δαβ＝1，如果α＝β；δαβ＝0，如果α≠β)，
Figure GDA00024921436600000214
和
Figure GDA00024921436600000215
是第(l,n,m)个tesseroid单元体几何中心坐标和密度；并且
Figure GDA00024921436600000216
C、在重力位、重力加速度和重力梯度张量的核矩阵为如下公式前提下：
Figure GDA00024921436600000217
Figure GDA00024921436600000218
Figure GDA00024921436600000219
单个观测点
Figure GDA00024921436600000220
上的重力位、加速度三分量、梯度张量则写为：
Figure GDA0002492143660000031
Figure GDA0002492143660000032
Figure GDA0002492143660000033
其中
Figure GDA0002492143660000034
是高度为r0上的测点的坐标；p＝1,2,…,Nλand
Figure GDA00024921436600000313
D、当地下三维场源体剖分成均匀的等间隔的网格单元，并且观测面位于一个平面，观测点与tesseroid单元体的中心点的位置一一对应，则公式(13)～(15)中的核矩阵写为：
Figure GDA0002492143660000035
Figure GDA0002492143660000036
Figure GDA0002492143660000037
E、场源和观测面的剖分是等间隔的，观测点的位置与tesseroid单元体中心点位置一一对应；则λp和λm′中的下标p,m可以直接用于加减运算，以Kv为例，有
Figure GDA0002492143660000038
并且，当m≤p时，
Figure GDA0002492143660000039
写为
Figure GDA00024921436600000310
F、综合考虑交叉等效性公式(22)和平移等效性公式(23)，重力位的核矩阵等效性关系可以简化为
Figure GDA00024921436600000311
Ky,Kxy和Kyz是关于(λ′-λ)的奇函数，因此相比于(24)式，需添加一个符号函数
Figure GDA00024921436600000312
其中
Figure GDA0002492143660000041
优选的，公式(13)～(15)更进一步写成矩阵相乘的形式
Kv·ρ＝V, (16)
Kα·ρ＝gα, (17)
Kαβ·ρ＝gαβ, (18)
其中ρ表示tesseroid单元体密度向量，维度为Ntess；V、gα和gαβ正演结果向量，维度为Nobs；Kv、Kα和Kαβ重力场正演核矩阵从公式(10)～(12)得到。
优选的，并且，以一层tesseroid单元模型为例，因此(22)式中变量r0和rl′被省略时，
用K(q,p；n,m)来表示tesseroid单元体
Figure GDA0002492143660000042
在观测点
Figure GDA0002492143660000043
上产生的重力响应，则Kv,Kx,Kz,Kxx,Kxz,Kyy,Kzz的核矩阵元素相等。
优选的，当m≥p时，先运用(22)式所示的交叉等效性，然后再运用(23)式所示的平移等效性。
优选的，球坐标系下模型目标函数写为：
Figure GDA0002492143660000044
其中αr,αs,
Figure GDA0002492143660000045
和αλ是衡量(29)式中四项相对重要性的权重因子；
Figure GDA0002492143660000046
ρref是参考模型；
w(r)是深度加权函数，球坐标系下三维重力场反演的深度加权函数写为
Figure GDA0002492143660000047
本申请具有以下有益效果：
本申请提出了三维重力场正演核矩阵等效性算法，并将该算法与自适应细剖分策略相结合，在保证正演精度的前提下大幅度提高正演计算效率和减小内存占用。在计算球坐标系下三维重力场正演时，不需要计算所有核矩阵元素，只需计算沿着第一个经度圈由
Figure GDA0002492143660000051
单元体所产生的所有核矩阵元素
Figure GDA0002492143660000052
Figure GDA0002492143660000053
其余核矩阵元素都可以通过等效关系求得，这就大大提高了计算效率。与此同时，存储核矩阵所需的内存也减小到传统方法所需内存的1/Nλ。对于更为一般的情况，场源可能剖分为多层，那么对于多层的核矩阵的计算方法与单层方法相同，只需多次重复单层的流程即可。
为了提高正演精度，本申请运用已成熟的自适应剖分策略。本申请指出，本申请所提出的核矩阵等效存储策略和自适应剖分策略可以结合起来。因为核矩阵等效性是用于提高计算效率，而自适应细剖分策略主要用于提高正演精度，两者是互不影响的。自适应剖分策略只需运用于第一个经度圈上的tesseroid单元体即可。
然后将该快速高精度正演算法用于三维重力场反演，重新构造了球坐标系下三维重力场反演目标函数和深度加权函数，以提高反演结果的可信度和分辨率。
为了证明本申请所提出方法的正确性，本申请给出的一个合成反演模型算例模型包含280*140*20个单元体，观测点数位280*140，本申请所提出的反演方法用时10小时，所占内存为0.88Gb，如果用传统方法反演如此大规模数据，在当前看来是不可能的。从反演结果来看，反演得到的异常体位置和真实模型位置大体一致，证明了本申请所提出方法的有效性。然后将本方法运用于月球质量瘤的反演，成功反演出雨海和澄海地区地下三维密度分布，反演结果和前人得出的莫霍面深度信息基本一致，进一步证明了本申请所提出方法的正确性。因此本方法是一种快速高精度的三维重力场反演方法。
本申请所提出的方法既适用于勘探尺度的资源勘查，水文环境监测，又适用于大尺度三维密度成像、研究地球壳幔结构、构造运动以及岩石圈变迁等，应用范围广泛。
除了上面所描述的目的、特征和优点之外，本申请还有其它的目的、特征和优点。下面将参照图，对本申请作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：
图1是局部球坐标系统及tesseroid单元体几何示意图；
图2是地下三维场源以及观测面剖分示意图；
图3是本申请优选实施例的核矩阵交叉等效性示意图；
图4是本申请优选实施例的核矩阵平移等效性示意图；
图5是本申请优选实施例的合成反演模型示意图；其中(a)为合成反演模型在z＝27.5km切面处的密度分布；(b)加入噪声的观测数据；(c)预测数据(反演结果正演重构数据)；(d)观测数据和预测数据之间的差异；
图6是图5(a)上的黑色虚线A-A′,B-B′,C-C′,D-D′所表示的4条断面的反演结果；其中(a)为A-A′,B-B′,C-C′,D-D′四个断面的真实密度分布断面图；(b)为A-A′,B-B′,C-C′,D-D′四个断面的反演所得密度分布断面图；
图7是本申请优选实施例用于处理月球局部(雨海盆地和澄海盆地)重力数据的过程图片，其中(a)雨海和澄海盆地地形起伏；(b)自由空气重力异常；(c)地形改正；(d)布格重力异常；图(a)中的AM为ApennineMountains的缩写；CM为Caucasus Mountains的缩写；
图8是图7(d)上的黑色虚线A-A′,B-B′,C-C′,D-D′所表示的4条断面的反演密度分布；图中黑色加粗实线表示前人得出的莫霍面深度；(a)为A-A′断面的反演密度分布，(b)为B-B′断面的反演密度分布,(c)为C-C′断面的反演密度分布,(d)为D-D′断面的反演密度分布。
具体实施方式
以下结合附图对本申请的实施例进行详细说明，但是本申请可以根据权利要求限定和覆盖的多种不同方式实施。
参见图1至图4，本申请一种基于3D-GLQ的球坐标系下重力场快速高精度反演方法，提出了三维重力场正演核矩阵等效性算法，并将该算法与自适应细剖分策略相结合，该方法的核心点有两处：
第一点：基于3D-GLQ的球坐标系下三维重力场快速高精度正演算法。
对勘探尺度而言，重力场正反演常在直角坐标系下进行，然而当研究区域较大或者研究整个地球岩石圈时，地球曲率的影响不得不考虑，因此需要在球坐标系下进行。在球坐标系下，场源通常被剖分为如图1所示的tesseroid单元体。由图1可以看出，tesseroid单元体是由三对曲面围成的：一对同球心的曲面(r1,r2)；一对子午线断面(λ1,λ2)；一对同轴圆锥平面
Figure GDA0002492143660000061
对于场源体外任一观测点P(坐标为
Figure GDA0002492143660000062
)，由第q个tesseroid单元体(其坐标为
Figure GDA0002492143660000063
)在该观测点处产生的重力位
Figure GDA0002492143660000064
重力加速度
Figure GDA0002492143660000065
和重力梯度张量
Figure GDA0002492143660000066
可以写为：
Figure GDA0002492143660000067
Figure GDA0002492143660000071
Figure GDA0002492143660000072
其中与r有关的量(r、r′、r1、r2)都表示半径，r为观测点P的半径，r′为第q个tesseroid单元体的半径，r1为第q个tesseroid单元体的下底面半径，r2为第q个tesseroid单元体上界面半径；与
Figure GDA0002492143660000073
有关的量
Figure GDA0002492143660000074
都表示纬度度，
Figure GDA0002492143660000075
为观测点P的维度，
Figure GDA0002492143660000076
为第q个tesseroid单元体的维度，
Figure GDA0002492143660000077
为第q个tesseroid单元体积分下限纬度，
Figure GDA0002492143660000078
为第q个tesseroid单元体积分上限纬度，单位为°，与λ相关的量(λ、λ′、λ1、λ2)都表示经度，λ为观测点P的经度，λ′为第q个tesseroid单元体的经度，λ1为第q个tesseroid单元体积分下限经度，r2为第q个tesseroid单元体积分上限经度，单位为°；
Figure GDA0002492143660000079
是第q个tesseroid单元体密度分布；G是重力常数；δαβ是克罗内克符号(当α＝β时，δαβ＝1；当α≠β时，δαβ＝0)，并且
Figure GDA00024921436600000710
其中Δx,Δy,Δz表示公式(1)-(3)中Δα和Δβ当α,β∈{x,y,z}时的表达式；κ,ι分别为公式(1)-(3)中的量，cosΨ为ι中的量。
由于公式(1)～(3)中的包含对λ'和
Figure GDA00024921436600000713
的椭圆积分，因此他们没有解析公式，需要通过数值求解的方法去离散。上述公式(1)-(3)可以用如下三维高斯勒让德算法求解：
Figure GDA00024921436600000711
其中
Figure GDA00024921436600000712
I、J和K分别为径向、经向和纬向上的高斯节点数，以及
Figure GDA0002492143660000081
是变换到每个tesseroid单元体上的高斯节点；
Figure GDA0002492143660000082
和ri'、λj'、
Figure GDA0002492143660000083
分别是在区间[-1,1]上第i、j、k个高斯系数和高斯节点。
使用上述3D GLQ算法计算球坐标系下三维重力场将非常耗时和消耗内存，因此本申请提出核矩阵等效存储策略。
首先将地下场源在径向、纬向、经向分别剖分成Nr、
Figure GDA0002492143660000084
Nλ段，记为
Figure GDA0002492143660000085
观测平面位于场源正上方，观测点数位
Figure GDA0002492143660000086
其剖分示意图如图2所示。图2中的Δr、
Figure GDA0002492143660000087
和Δλ分别表示在径向、纬向和经向的剖分间隔。
在每一个tesseroid单元体中，其密度都看作是常数，因此可以从三维体积分号中提取出来，公式(1)～(3)可以写成单个tesseroid累加求和的形式
Figure GDA0002492143660000088
Figure GDA0002492143660000089
Figure GDA00024921436600000810
其中
Figure GDA00024921436600000811
和
Figure GDA00024921436600000812
是第(l,n,m)个tesseroid单元体几何中心坐标和密度；
本申请记
Figure GDA00024921436600000813
Figure GDA00024921436600000814
Figure GDA00024921436600000815
为重力位、重力加速度和重力梯度张量的核矩阵，这些核矩阵只与tesseroid单元体的大小和该单元体到观测点之间的几何距离有关，而与tesseroid单元体的密度分布情况无关。
因此单个观测点
Figure GDA0002492143660000091
上的重力位、加速度三分量、梯度张量可以写为
Figure GDA0002492143660000092
Figure GDA0002492143660000093
Figure GDA0002492143660000094
其中
Figure GDA0002492143660000095
是高度为r0上的测点的坐标；p＝1,2,…,Nλand
Figure GDA00024921436600000910
公式(13)～(15)可以更进一步写成矩阵相乘的形式
Kv·ρ＝V, (16)
Kα·ρ＝gα, (17)
Kαβ·ρ＝gαβ, (18)
其中ρ表示tesseroid单元体密度向量，维度为Ntess；V、gα和gαβ正演结果向量，维度为Nobs；Kv、Kα和Kαβ重力场正演核矩阵，可以从公式(10)～(12)得到。
当地下三维场源体剖分成均匀的等间隔的网格单元，并且观测面位于一个平面，观测点与tesseroid单元体的中心点的位置一一对应(如图2所示)，那么在正演核矩阵中将存在一系列的等效关系，本申请称这些等效关系为核矩阵等效性。考虑到公式(4)都是和(λ′-λ)有关的量，因此在公式(13)～(15)中的核矩阵可以写为
Figure GDA0002492143660000096
Figure GDA0002492143660000097
Figure GDA0002492143660000098
根据公式(1)～(4)，本申请发现Kv,Kx,Kz,Kxx,Kxz,Kyy和Kzz是的(λ′-λ)偶函数。另一方面，由于场源和观测面的剖分都是等间隔的，观测点的位置与tesseroid单元体中心点位置一一对应。因此λp和λm′中的下标p,m可以直接用于加减运算，以Kv为例，有
Figure GDA0002492143660000099
本申请称这种关系为交叉等效性，图3是这种关系的几何示意图。为了更为清晰，图3仅以一层模型为例，因此(22)式中变量r0和rl′可以被省略。将该层剖分为6*6的tesseroid单元体，观测点个数为6*6，位于单元体正上方。本申请使用K(q,p；n,m)来表示tesseroid单元体
Figure GDA0002492143660000101
在观测点
Figure GDA0002492143660000102
上产生的重力响应。例如，核矩阵元素K(3,2；2,5)表示tesseroid单元体
Figure GDA0002492143660000103
在测点
Figure GDA0002492143660000104
上产生的重力场响应。那么，对于Kv,Kx,Kz,Kxx,Kxz,Kyy和Kzz而言，从公式(1)～(4)可以看出，这些分量是(λ′-λ)的偶函数，因此由交叉的虚线(黑色或者灰色)所表示的重力场核矩阵元素大小相等；但对于Ky,Kxy和Kyz而言，由于这三个分量是(λ′-λ)的奇函数，因此由交叉的虚线(黑色或者灰色)所表示的重力场核矩阵元素相差一个负号。
上述所及的是交叉等效性，另一方面，本申请考虑另一等效关系，分两种情况讨论：
当m≤p时，
Figure GDA0002492143660000105
还可以写为
Figure GDA0002492143660000106
本申请称这种关系为平移等效性，图4是其几何示意图。图4模型和图3中一样，沿着平行的虚线(黑色或者灰色)所表示的核矩阵元素大小相等。而当m≥p时，由于无法直接运用公式(23)，本申请可以先运用(22)式所示的交叉等效性，转换成m≤p的情况，然后在使用平移等效性公式(23)。
综合考虑交叉等效性公式(22)和平移等效性公式(23)，重力位的核矩阵等效性关系可以简化为
Figure GDA0002492143660000107
对于Ky,Kxy和Kyz而言，他们是关于(λ′-λ)的奇函数，因此相比于(24)式，本申请需要额外添加一个符号函数
Figure GDA0002492143660000108
其中
Figure GDA0002492143660000109
该符号函数为人为定义的符号函数，当该符号函数的取值为1或者-1，用于综合Ky,Kxy和Kyz分量中关于(λ′-λ)奇函数所引起的符号差异。
至此，在计算球坐标系下三维重力场正演时，本申请不需要计算所有核矩阵元素，只需计算沿着第一个经度圈由
Figure GDA0002492143660000111
单元体所产生的所有核矩阵元素
Figure GDA0002492143660000112
Figure GDA0002492143660000113
其余核矩阵元素都可以通过等效关系求得，这就大大提高了计算效率。与此同时，存储核矩阵所需的内存也减小到传统方法所需内存的1/Nλ。对于更为一般的情况，场源可能剖分为多层，那么对于多层的核矩阵的计算方法与单层方法相同，只需多次重复单层的流程即可。再者，如果场源剖分与观测点位置也满足如图2中所述关系的话，核矩阵等效性也可以用于泰勒级数展开法重力场正演计算。
为了提高正演精度，本申请运用已成熟的自适应剖分策略。本申请所提出的核矩阵等效存储策略和自适应剖分策略可以结合起来。因为核矩阵等效性是用于提高计算效率，而自适应细剖分策略主要用于提高正演精度，两者是互不影响的。自适应剖分策略只需运用于第一个经度圈上的tesseroid单元体即可。
第二点：球坐标系下三维重力场反演方法的建立。
这里本申请将上面提出的快速高精度三维重力场正演方法运用到三维反演中去，运用吉洪诺夫正则化反演方法，通过最小化数据拟合差与模型拟合差，来寻求最小构造解：
minimizeΦ(ρ)＝Φd+λΦm, (27)
其中Φd是数据拟合差函数；Φm是模型目标函数；λ是正则化因子，用于平衡Φd与Φm之间的权重。正则化因子的选择有多重方法，本申请使用L曲线法。
地球物理界常用的数据拟合函数可以写为
Figure GDA0002492143660000114
其中dobs＝(d1,…,dN)T观测数据向量；dpre是预测数据；Wd＝diag{1/σ1,…,1/σN}；σi是第i个数据的误差标准差。
球坐标系下模型目标函数可以写为
Figure GDA0002492143660000121
其中αr,αs,
Figure GDA0002492143660000122
和αλ是衡量(29)式中四项相对重要性的权重因子；
Figure GDA0002492143660000123
ρref是参考模型；w(r)是深度加权函数，球坐标系下三维重力场反演的深度加权函数可以写为
Figure GDA0002492143660000124
在数据拟合函数(28)中，本申请所提出的快速算法可以用于计算dpre。通过最小化公式(27)，反演问题转化为一个线性方程组的求解问题，在本申请中，本申请使用成熟的共轭梯度法求解该方程组。
综上，本申请基于3D GLQ，提出正演核矩阵等效存储策略，并结合自适应剖分策略，然后将快速高精度正演方法用于球坐标系下三维重力场反演。实施的具体步骤如下：
步骤1、从网格化grid文件中输入观测数据，程序自动判断观测数据个数，起始范围，并计算出剖分间隔，以及每一个测点的坐标信息等；
步骤2、设置观测高度，设置观测数据误差均方差；
步骤3、设置反演参数，程序需要输入反演深度的起始范围，该范围可根据该地区的地质、区域构造、前人研究结果得到；
步骤4、运行程序，程序按照前述公式进行运算，得到最终反演结果。
实施例一、
本实施例给出一个合成反演模型，该模型包括两个具有相同密度600kg/m3的块体，如图5a和6a所示。用本申请所提出的方法正演计算了月球距参考半径(1738km)10km高度处的重力异常，然后再加入1％的随机噪声，作为观测数据(图5b)。模型精度方向起始范围为-35°到35°；纬度方向起始范围为15°到50°；深度方向上从0km到100km。这一模型剖分成280×140×20个以0.25°,0.25°,5km为间隔的tesseroid单元体。
用本申请所提出的球坐标系下三维重力场快速高精度正演算法反演该合成重力场数据，预测数据和其残差分别如图5c和5d所示。反演得到的密度分布如图6b所示。如图5d所示，观测数据和预测数据拟合很好，这就证明了本申请所提出的正反演方法的有效性。
在这个反演测试中，基于本申请所提出的方法的计算时间为10小时，而用传统的方法反演的话，计算时间预计1000小时左右，可见计算效率的提升是明显的。本反演所占内存仅为0.88Gb，而传统反演方法预计消耗内存245.86Gb。
实施例二、
为了进一步说明本申请所提出的方法的正确性和高效性，接下来本申请给出一个真实反演案例，将所提出方法运用于月球局部重力数据。本申请选择雨海盆地和澄海盆地作为研究区域，经度范围是35°W到35°E，纬度范围是15°N到50°N，如图7a所示。
这里本申请使用月球最新重力场模型GL1500E以及最新的月球地形模型LRO_LTM05_2050来计算布格重力异常。自由空气重力异常如图7b所示，地形改正结果如图7c所示。那么布格重力异常可以通过自由空气重力异常(图7b)减去地形改正(图7c)得到，如图7d所示。在这里，地形校正所用的密度是2560kg/m3，观测高度是10km。只用了模型的前720阶，对应于空间分辨率为0.25°。
本申请将地下场源剖分成280×140个等间隔的tesseroid单元体，根据前人研究基础，在深度方向上模型范围设置为0km到100km，剖分成20段，间隔为5km。本申请对布格重力异常数据做三维重力场反演，反演结果如图8所示。沿着4条断面的反演结果(图8)表明，高密度异常体主要分布在雨海和澄海的中心位置，其顶界面与前人研究的莫霍面起伏深度一致。在雨海地区，高密度异常体在深度分布范围大概为18-40km，在澄海地区，高密度异常深度分布范围大概为13-60km。
综上，本申请在正演方面提出了核矩阵等效存储策略，同时加入了成熟的自适应剖分技术，使得在保证正演精度的同时，大幅度提高正演计算效率，以及减小内存占用，使得大规模三维重力场反演成为可能。在反演方面，给出了球坐标系下三维重力场反演目标函数以及深度加权函数，以适应地球曲率的影响，本发明所提出的快速高精度正演方法直接用于反演，大大提高了大规模反演的可行性和可信度。
以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

Claims (5)

1.一种基于3D-GLQ的球坐标系下重力场正演方法，其特征在于，包括以下步骤：
A、将地下场源在径向、纬向、经向分别剖分成Nr、
Figure FDA0002436718230000011
Nλ段，记为
Figure FDA0002436718230000012
观测平面位于场源正上方，观测点数位
Figure FDA0002436718230000013
B、场源剖分成的tesseroid单元体是由三对曲面围成：一对同球心的曲面(r1,r2)、一对子午线断面(λ1,λ2)、一对同轴圆锥平面
Figure FDA0002436718230000014
将每一个tesseroid单元体的密度看作是常数，对于场源体外一计算点
Figure FDA0002436718230000015
由第q个tesseroid单元体
Figure FDA0002436718230000016
在该观测点处产生的重力位
Figure FDA0002436718230000017
重力加速度
Figure FDA0002436718230000018
和重力梯度张量
Figure FDA0002436718230000019
可以写为：
Figure FDA00024367182300000110
Figure FDA00024367182300000111
Figure FDA00024367182300000112
其中，α,β∈{x,y,z}；
Figure FDA00024367182300000113
是密度分布；G是重力常数；δαβ是克罗内克符号，
Figure FDA00024367182300000114
Figure FDA00024367182300000115
和
Figure FDA00024367182300000116
是第(l,n,m)个tesseroid单元体几何中心坐标和密度；并且
Figure FDA00024367182300000117
C、在重力位、重力加速度和重力梯度张量的核矩阵为如下公式前提下：
Figure FDA00024367182300000118
Figure FDA0002436718230000021
Figure FDA0002436718230000022
单个观测点
Figure FDA0002436718230000023
上的重力位、加速度三分量、梯度张量则写为：
Figure FDA0002436718230000024
Figure FDA0002436718230000025
Figure FDA0002436718230000026
其中
Figure FDA0002436718230000027
是高度为r0上的测点的坐标；
Figure FDA00024367182300000214
D、当地下三维场源体剖分成均匀的等间隔的网格单元，并且观测面位于一个平面，观测点与tesseroid单元体的中心点的位置一一对应，则公式(13)～(15)中的核矩阵写为：
Figure FDA0002436718230000028
Figure FDA0002436718230000029
Figure FDA00024367182300000210
E、场源和观测面的剖分是等间隔的，观测点的位置与tesseroid单元体中心点位置一一对应；则λp和λm′中的下标p,m能直接用于加减运算，对于重力位核矩阵Kv，有
Figure FDA00024367182300000211
并且，当m≤p时，
Figure FDA00024367182300000212
写为
Figure FDA00024367182300000213
F、综合考虑交叉等效性公式(22)和平移等效性公式(23)，重力位的核矩阵等效性关系可以简化为
Figure FDA0002436718230000031
Ky,Kxy和Kyz是关于(λ′-λ)的奇函数，因此相比于(24)式，需添加一个符号函数
Figure FDA0002436718230000032
其中
Figure FDA0002436718230000033
2.根据权利要求1所述的基于3D-GLQ的球坐标系下重力场正演方法，其特征在于，公式(13)～(15)更进一步写成矩阵相乘的形式
Kv·ρ＝V, (16)
Kα·ρ＝gα, (17)
Kαβ·ρ＝gαβ, (18)
其中ρ表示tesseroid单元体密度向量，维度为Ntess；V、gα和gαβ正演结果向量，维度为Nobs；Kv、Kα和Kαβ重力场正演核矩阵从公式(10)～(12)得到。
3.根据权利要求1所述的基于3D-GLQ的球坐标系下重力场正演方法，其特征在于，对于每一层tesseroid单元体，省略(22)式中的r0和r′l，
用K(q,p；n,m)来表示tesseroid单元体
Figure FDA0002436718230000034
在观测点
Figure FDA0002436718230000035
上产生的重力响应，则Kv,Kx,Kz,Kxx,Kxz,Kyy,Kzz的核矩阵元素相等。
4.根据权利要求1所述的基于3D-GLQ的球坐标系下重力场正演方法，其特征在于，当m≥p时，先运用(22)式所示的交叉等效性，然后再运用(23)式所示的平移等效性。
5.基于权利要求1所述的基于3D-GLQ的球坐标系下重力场正演方法的三维反演方法，其特征在于，球坐标系下模型目标函数写为：
Figure FDA0002436718230000041
其中αr,αs,
Figure FDA0002436718230000042
和αλ是衡量(29)式中四项相对重要性的权重因子；
Figure FDA0002436718230000043
ρref是参考模型；
w(r)是深度加权函数，球坐标系下三维重力场反演的深度加权函数写为
Figure FDA0002436718230000044
CN201811481164.5A 2018-12-05 2018-12-05 基于3d-glq的球坐标系下重力场正演方法及三维反演方法 Active CN110045432B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811481164.5A CN110045432B (zh) 2018-12-05 2018-12-05 基于3d-glq的球坐标系下重力场正演方法及三维反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811481164.5A CN110045432B (zh) 2018-12-05 2018-12-05 基于3d-glq的球坐标系下重力场正演方法及三维反演方法

Publications (2)

Publication Number Publication Date
CN110045432A CN110045432A (zh) 2019-07-23
CN110045432B true CN110045432B (zh) 2020-06-30

Family

ID=67273256

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811481164.5A Active CN110045432B (zh) 2018-12-05 2018-12-05 基于3d-glq的球坐标系下重力场正演方法及三维反演方法

Country Status (1)

Country Link
CN (1) CN110045432B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110501751B (zh) * 2019-08-23 2021-03-16 东北大学 一种基于多分量梯度数据联合和深度加权的相关成像方法
CN111400654B (zh) * 2020-03-13 2022-03-18 中南大学 一种基于Toplitze核矩阵的重力场快速正演方法及反演方法
CN111628822B (zh) * 2020-06-04 2021-03-26 清华大学 一种紫外光通信非视距链路中单次散射路径损耗的近似计算方法
CN112346139B (zh) * 2020-10-15 2021-12-21 中国地质大学（武汉） 一种重力数据多层等效源延拓与数据转换方法
CN112596106B (zh) * 2020-11-09 2024-04-26 中国人民武装警察部队黄金第五支队 一种球坐标系下重震联合反演密度界面分布的方法
CN113109883B (zh) * 2021-05-26 2023-01-03 湖南科技大学 基于等参变换全球离散网格球坐标下卫星重力场正演方法
CN113311494B (zh) * 2021-05-26 2022-04-26 中南大学 一种卫星重力场反演方法
CN113740915B (zh) * 2021-07-28 2024-04-19 中国人民武装警察部队黄金第五支队 一种球坐标系重力和接收函数联合反演地壳结构参数的方法
CN113985490B (zh) * 2021-09-22 2023-05-05 中国人民解放军战略支援部队信息工程大学 利用地形和地壳密度数据进行地表重力仿真的方法及装置
CN114966878B (zh) * 2022-04-12 2023-04-14 成都理工大学 一种基于混合范数和互相关约束的三维重力场反演方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106897425A (zh) * 2017-02-24 2017-06-27 中国科学院电子学研究所 一种地球重力场数据的三维可视化方法
CN107766691B (zh) * 2017-09-13 2019-12-31 中国地震局地震预测研究所 Grace卫星重力场球谐系数去相关的方法及电子设备
CN108490496B (zh) * 2018-03-26 2019-11-08 中国石油化工股份有限公司 基于拟径向基函数神经网络的重力场密度反演方法

Also Published As

Publication number Publication date
CN110045432A (zh) 2019-07-23

Similar Documents

Publication Publication Date Title
CN110045432B (zh) 基于3d-glq的球坐标系下重力场正演方法及三维反演方法
CN111400654B (zh) 一种基于Toplitze核矩阵的重力场快速正演方法及反演方法
CN109375280B (zh) 一种球坐标系下重力场快速高精度正演方法
Yun et al. Constraints on magma chamber geometry at Sierra Negra Volcano, Galápagos Islands, based on InSAR observations
Denker Regional gravity field modeling: theory and practical results
Hill et al. Combination of geodetic observations and models for glacial isostatic adjustment fields in Fennoscandia
CN113283802B (zh) 一种复杂艰险山区滑坡危险性评估方法
Song et al. Bridging the gap between geophysics and geology with generative adversarial networks
Zhao et al. Efficient 3‐D large‐scale forward modeling and inversion of gravitational fields in spherical coordinates with application to lunar mascons
CN113341476B (zh) 基于海底地形-重力联合提高海洋重力空间分辨率的方法
Bevilacqua et al. Radial interpolation of GPS and leveling data of ground deformation in a resurgent caldera: application to Campi Flegrei (Italy)
CN113514900B (zh) 基于密度约束的球坐标系重力和重力梯度联合反演方法
CN114943178A (zh) 一种三维地质模型建模方法、装置及计算机设备
Agoshkov et al. Variational assimilation of observation data in the mathematical model of the Black Sea taking into account the tide-generating forces
Gangui et al. Cosmic microwave background bispectrum from active models of structure formation
Yin et al. Estimating gravity changes caused by crustal strain: application to the Tibetan Plateau
Yin et al. A fast 3D gravity forward algorithm based on circular convolution
Din et al. Assessment of global geopotential models for modelling Malaysia Marine Geoid
Gvishiani et al. Gis-oriented database for the system analysis and prediction of the geodynamic stability of the Nizhne-Kansky massif
Zygmunt et al. Assessment of continental hydrosphere loading using GNSS measurements
Petrov et al. Spatial-temporal three-dimensional GIS modeling
CN114966878B (zh) 一种基于混合范数和互相关约束的三维重力场反演方法
Yan et al. Research and construction of a global hexagonal marine gravity gradient reference map for navigation
Radanović et al. Accuracy assessment and comparison of interpolation methods on geoid models
Yang et al. A novel phase unwrapping method for low coherence interferograms in coal mining areas based on a fully convolutional neural network

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant